超声原始数据重构成B扫成像的MATLAB实现

基本原理

1. 原始RF数据：超声探头接收的原始射频信号
2. 包络检测：提取RF信号的幅度信息
3. 对数压缩：压缩动态范围以适应显示
4. 时间增益补偿(TGC)：补偿深度相关的信号衰减
5. 扫描转换：将数据映射到图像坐标系

MATLAB代码

function ultrasound_bscan_reconstruction()
    % 超声B扫成像重建
    close all; clc;
    
    % 生成模拟超声RF数据
    [rf_data, fs, c] = simulate_ultrasound_data();
    
    % 显示原始RF数据
    display_raw_rf(rf_data, fs, c);
    
    % 重建B超图像
    b_mode_image = reconstruct_bscan(rf_data, fs, c);
    
    % 显示最终B超图像
    display_bscan_image(b_mode_image);
end

function [rf_data, fs, c] = simulate_ultrasound_data()
    % 模拟超声RF数据生成
    c = 1540;               % 声速 (m/s)
    fs = 100e6;             % 采样频率 (Hz)
    num_samples = 1500;     % 每条扫描线采样点数
    num_lines = 128;        % 扫描线数量
    depth = 0.04;           % 成像深度 (m)
    
    % 创建时间和空间轴
    t = (0:num_samples-1)/fs;   % 时间轴
    z = t * c / 2;              % 深度轴
    
    % 生成扫描线位置 (线性阵列)
    x_pos = linspace(-0.02, 0.02, num_lines);
    
    % 初始化RF数据矩阵
    rf_data = zeros(num_samples, num_lines);
    
    % 添加组织模型中的散射点
    scatterers = [
        % x(m), z(m), 反射强度
        0.000, 0.010, 1.0;
        0.005, 0.015, 0.8;
        -0.004, 0.020, 0.9;
        0.000, 0.030, 1.2;
        0.008, 0.025, 0.7;
        -0.007, 0.018, 0.6;
    ];
    
    % 添加高反射界面 (模拟器官边界)
    z_boundary = 0.012:0.002:0.038;
    for i = 1:length(z_boundary)
        scatterers = [scatterers; 
            linspace(-0.015, 0.015, 20)', repmat(z_boundary(i), 20, 1), repmat(1.5, 20, 1)];
    end
    
    % 为每条扫描线生成RF信号
    for line = 1:num_lines
        % 探头脉冲响应 (高斯调制正弦波)
        t_pulse = (-10:10)/fs * 5;
        pulse = exp(-(t_pulse*1e6).^2) .* sin(2*pi*5e6*t_pulse);
        
        % 初始化扫描线信号
        line_signal = zeros(num_samples, 1);
        
        % 添加散射点响应
        for s = 1:size(scatterers, 1)
            dx = scatterers(s, 1) - x_pos(line);
            dz = scatterers(s, 2);
            distance = sqrt(dx^2 + dz^2);  % 到散射点的距离
            
            % 计算往返时间
            time_of_flight = 2 * distance / c;
            
            % 找到最近的采样点
            sample_idx = round(time_of_flight * fs) + 1;
            
            % 确保在有效范围内
            if sample_idx > 10 && sample_idx < num_samples - 10
                % 添加脉冲响应
                start_idx = sample_idx - 10;
                end_idx = sample_idx + 10;
                line_signal(start_idx:end_idx) = line_signal(start_idx:end_idx) + ...
                    scatterers(s, 3) * pulse';
            end
        end
        
        % 添加噪声和系统增益变化
        line_signal = line_signal + 0.02 * randn(size(line_signal));
        rf_data(:, line) = line_signal;
    end
    
    % 添加系统增益随深度变化 (模拟TGC效果)
    tgc_gain = linspace(1, 5, num_samples)';
    rf_data = rf_data .* repmat(tgc_gain, 1, num_lines);
end

function display_raw_rf(rf_data, fs, c)
    % 显示原始RF数据
    figure('Name', '原始RF数据', 'Position', [100, 100, 900, 600]);
    
    % 计算深度轴
    num_samples = size(rf_data, 1);
    t = (0:num_samples-1)/fs;
    depth = t * c / 2;
    
    % 显示单条RF信号
    subplot(2, 2, [1, 3]);
    plot(t*1e6, rf_data(:, round(size(rf_data, 2)/2)));
    xlabel('时间 (\mus)');
    ylabel('幅度');
    title('单条扫描线RF信号');
    grid on;
    
    % 显示RF数据图像
    subplot(2, 2, 2);
    imagesc([], depth*1000, rf_data);
    colormap(gray);
    colorbar;
    title('RF数据矩阵');
    ylabel('深度 (mm)');
    xlabel('扫描线索引');
    axis tight;
    
    % 显示RF数据幅度谱
    subplot(2, 2, 4);
    [pxx, f] = pwelch(rf_data(:, round(size(rf_data, 2)/2)), [], [], [], fs);
    plot(f/1e6, 10*log10(pxx));
    xlabel('频率 (MHz)');
    ylabel('功率谱密度 (dB/Hz)');
    title('RF信号频谱');
    grid on;
    xlim([0, 20]);
end

function b_mode_image = reconstruct_bscan(rf_data, fs, c)
    % 重建B超图像
    % 参数设置
    f0 = 5e6;       % 中心频率 (Hz)
    bw = 0.7;       % 带宽系数
    
    % 获取数据尺寸
    [num_samples, num_lines] = size(rf_data);
    
    % 步骤1: 带通滤波
    rf_filtered = bandpass_filter(rf_data, fs, f0, bw);
    
    % 步骤2: 包络检测 (希尔伯特变换)
    analytic_signal = hilbert(rf_filtered);
    envelope = abs(analytic_signal);
    
    % 步骤3: 对数压缩
    log_compressed = 20 * log10(envelope + eps);  % 加eps避免log(0)
    
    % 步骤4: 时间增益补偿 (TGC)
    depth = ((0:num_samples-1)/fs) * c / 2;
    tgc = 1.0 + 1.5 * depth / max(depth);  % 线性增益补偿
    tgc_corrected = log_compressed .* tgc';
    
    % 步骤5: 降采样 (可选)
    dec_factor = 2;
    b_mode_downsampled = tgc_corrected(1:dec_factor:end, :);
    
    % 步骤6: 扫描转换 (此处为线性阵列，直接显示)
    b_mode_image = b_mode_downsampled;
end

function filtered_data = bandpass_filter(data, fs, f0, bw)
    % 带通滤波器设计
    nyquist = fs/2;
    f_low = max(0, f0*(1 - bw/2));    % 低频截止
    f_high = min(nyquist, f0*(1 + bw/2)); % 高频截止
    
    % 设计滤波器
    [b, a] = butter(4, [f_low, f_high]/nyquist, 'bandpass');
    
    % 应用滤波器
    filtered_data = filtfilt(b, a, data);
end

function display_bscan_image(b_mode_image)
    % 显示B超图像
    figure('Name', 'B超图像重建', 'Position', [100, 100, 1000, 800]);
    
    % 归一化图像
    img = b_mode_image - min(b_mode_image(:));
    img = img / max(img(:));
    
    % 显示图像
    imagesc(img);
    colormap(gray);
    title('重建的B超图像');
    xlabel('扫描线索引');
    ylabel('深度采样点');
    
    % 添加深度刻度 (示例值)
    depth_ticks = linspace(0, size(b_mode_image, 1), 6);
    depth_labels = round(linspace(0, 40, 6)); % 假设最大深度40mm
    set(gca, 'YTick', depth_ticks, 'YTickLabel', depth_labels);
    ylabel('深度 (mm)');
    
    % 添加颜色条
    colorbar;
    
    % 添加图像处理选项
    uicontrol('Style', 'pushbutton', 'String', '增强对比度', ...
        'Position', [20 20 100 30], 'Callback', @enhance_contrast);
    uicontrol('Style', 'pushbutton', 'String', '平滑处理', ...
        'Position', [140 20 100 30], 'Callback', @apply_smoothing);
    uicontrol('Style', 'pushbutton', 'String', '边缘增强', ...
        'Position', [260 20 100 30], 'Callback', @edge_enhancement);
    
    % 回调函数
    function enhance_contrast(~, ~)
        img_enhanced = imadjust(img);
        imagesc(img_enhanced);
        colormap(gray);
        title('对比度增强的B超图像');
    end

    function apply_smoothing(~, ~)
        h = fspecial('gaussian', [5 5], 1.5);
        img_smoothed = imfilter(img, h);
        imagesc(img_smoothed);
        colormap(gray);
        title('平滑处理后的B超图像');
    end

    function edge_enhancement(~, ~)
        img_edges = img - imgaussfilt(img, 2) + 0.5;
        img_edges(img_edges < 0) = 0;
        img_edges(img_edges > 1) = 1;
        imagesc(img_edges);
        colormap(gray);
        title('边缘增强的B超图像');
    end
end

1. 模拟超声原始数据生成

• 创建模拟组织环境（散射点模型）
• 生成超声脉冲响应（高斯调制正弦波）
• 模拟扫描线采集过程
• 添加时间增益补偿(TGC)和随机噪声

2. 原始RF数据显示

• 单条扫描线RF信号波形
• RF数据矩阵图像
• RF信号频谱分析

3. B扫成像重建处理

• 带通滤波：提取探头频带内的信号
• 包络检测：使用希尔伯特变换提取信号幅度
• 对数压缩：压缩动态范围（40-60dB）
• 时间增益补偿：补偿深度相关的信号衰减
• 降采样：减少数据量（可选步骤）

4. B超图像显示

• 显示重建的B超图像
• 添加深度刻度
• 交互式图像处理按钮：1对比度增强 2平滑处理 3边缘增强

参考 超声原始数据的重构成B扫成像的方法 youwenfan.com/contentcsa/65820.html

重建过程关键步骤详解

1. 带通滤波：

   • 目的：去除低频噪声和高频干扰
   • 方法：设计Butterworth带通滤波器
   • 实现：bandpass_filter函数

2. 包络检测：

   • 目的：提取RF信号的幅度信息

   • 方法：希尔伯特变换获取解析信号

   • 数学表达：

     analytic_signal = hilbert(rf_signal);
     envelope = abs(analytic_signal);
     

3. 对数压缩：

   • 目的：将大动态范围压缩到适合显示的0-255范围

   • 方法：取对数并缩放

   • 公式：

     log_compressed = 20 * log10(envelope + eps);
     

4. 时间增益补偿(TGC)：

   • 目的：补偿深度增加导致的信号衰减

   • 方法：应用与深度成正比的增益曲线

   • 实现：

     tgc = 1.0 + 1.5 * depth / max_depth;
     tgc_corrected = log_compressed .* tgc';
     

5. 扫描转换：

   • 目的：将数据映射到图像坐标系
   • 注意：本代码模拟线性阵列探头，直接显示数据矩阵
   • 对于扇形扫描，需要插值到笛卡尔坐标系

运行结果

运行此代码将生成四个图形窗口：

1. 原始RF数据显示（包含RF信号波形、数据矩阵和频谱）
2. 重建的B超图像（显示组织结构和散射点）
3. 交互式界面（提供图像增强选项）

B超图像将显示模拟的器官边界和散射点，您可以使用界面上的按钮进行图像增强处理。

这个实现提供了超声成像的基本原理演示，实际临床应用中的系统会更加复杂，包含更多优化和校正步骤。