当前位置：首页 > news >正文

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

news 2026/6/14 18:25:52

基本原理

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

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);
endfunction [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';endend% 添加噪声和系统增益变化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);
endfunction 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]);
endfunction 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;
endfunction 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);
endfunction 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)); % 假设最大深度40mmset(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超图像');endfunction apply_smoothing(~, ~)h = fspecial('gaussian', [5 5], 1.5);img_smoothed = imfilter(img, h);imagesc(img_smoothed);colormap(gray);title('平滑处理后的B超图像');endfunction 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扫成像的方法 www.youwenfan.com/contentcno/65820.html

重建过程关键步骤详解

带通滤波：
- 目的：去除低频噪声和高频干扰
- 方法：设计Butterworth带通滤波器
- 实现：bandpass_filter函数
包络检测：
- 目的：提取RF信号的幅度信息
- 方法：希尔伯特变换获取解析信号
- 数学表达：
```
analytic_signal = hilbert(rf_signal);
envelope = abs(analytic_signal);
```
对数压缩：
- 目的：将大动态范围压缩到适合显示的0-255范围
- 方法：取对数并缩放
- 公式：
```
log_compressed = 20 * log10(envelope + eps);
```
时间增益补偿(TGC)：
- 目的：补偿深度增加导致的信号衰减
- 方法：应用与深度成正比的增益曲线
- 实现：
```
tgc = 1.0 + 1.5 * depth / max_depth;
tgc_corrected = log_compressed .* tgc';
```
扫描转换：
- 目的：将数据映射到图像坐标系
- 注意：本代码模拟线性阵列探头，直接显示数据矩阵
- 对于扇形扫描，需要插值到笛卡尔坐标系