简介光声峰峰值成像利用光吸收产生的超声信号重建组织内部光吸收分布在生物医学光学成像与病变识别中具有实用价值。面向光声成像研究者、生物医学工程相关专业学生以及需要快速构建成像算法的开发者这份MATLAB资源提供了一套完整的峰峰值幅值成像代码包。压缩包内共有4个m文件整体仅3KB分别对应信号预处理与滤波、主流程控制、三维曲面可视化以及数据翻转平滑等模块代码精简、分工明确便于直接阅读和二次修改。已有433人学习下载。通过研读和运行这些脚本读者能够掌握光声信号去噪、峰峰值特征提取、三维图像重建等关键操作理解从原始声波信号到光吸收分布图的完整处理链条还可根据自身数据调整参数以优化成像质量适合作为光声成像入门与算法验证的参考资料。1. 光声峰峰值成像为什么用峰峰值而不是最大值或均值做光声成像的同学应该遇到过这样一个选择从采集到的射频RF信号重建图像时到底取信号的哪个特征来代表这个像素有人直接取最大值有人计算均值还有人在做“幅值成像”时下意识地用了信号的峰峰值。区别在于光声信号本质上是激光激发产生的超声波瞬态波它通常是双极性振荡包含着正负极值。直接取最大值只关注了正向半波漏掉了负向半波的信息而取均值几乎不能反映信号的瞬态强度。峰峰值Peak-to-Peak, PP定义为信号在某一时间窗口内最大值与最小值之差。在光声峰峰值成像中这个值表达了该空间位置处光声信号的绝对振荡幅度对直流偏移不敏感也天然抑制了缓慢变化的背景漂移。相比最大值、均值或均方根峰峰值能在不依赖绝对基线的情况下刻画吸收体的浓度和尺寸而且计算简单、无需拟合。这篇文章就围绕这个主题展开从信号模型到MATLAB实现再到参数调优和三维扩展。适合正在写光声信号处理脚本、做图像重建或调试采集系统的工程师阅读哪怕是刚接触MATLAB光声数据分析的人也可以照着例子跑通。2. 光声信号模型与峰峰值特征提取原理2.1 光声信号的产生与探测光声成像的原理并不复杂短脉冲激光照射吸收体如血管中的血红蛋白吸收体快速热膨胀产生压力波这个压力波被超声换能器接收形成时域信号我们称之为一个A-line。典型的光声信号类似于带阻尼的振荡正半周和负半周幅度可能不对称但整体包络反映了吸收体的光吸收分布。由于换能器带宽有限信号并不是理想的狄拉克脉冲而是一个有限频宽的瞬态波。这种信号的特征参数很多。对于一个采集到的RF数据矩阵行代表时间采样列代表空间位置一维扫描时是一个A-line索引二维B-scan时是扫描角度或x位置三维体积数据则是x和y两个方向的排列。光声峰峰值成像要做的就是沿着时间轴每一列计算信号的峰值到峰值的幅度再把这个幅度映射到对应的空间坐标上形成强度图像。2.2 峰峰值Peak-to-Peak的定义与数学表达设光声信号为 ( s(t) )其采样序列为 ( s[n], n1,\dots,N )。对应一个空间点的峰峰值定义为[ PP \max_{n \in W} s[n] - \min_{n \in W} s[n] ]其中 ( W ) 是所选的时间窗口可以是整个A-line长度也可以是激光发射后的一个特定时间区间。窗口的选择直接决定了成像的深度范围窗口起点对应探测起始深度窗口长度对应成像深度厚度。在MATLAB中这个计算可以用最朴素的循环完成但更好的做法是使用max和min函数在指定维度上作用。对于二维矩阵data(Nt, Nx)一行代码就能完成pp max(data, [], 1) - min(data, [], 1);这里的max(data, [], 1)是沿第一个维度时间维度求最大值返回一个1×Nx的行向量。min同理。两个向量相减得到每个A-line的峰峰值。这个操作没有显式循环效率很高。不过要注意直接对原始RF信号做max和min有一个陷阱信号中的随机噪声会产生额外的尖峰。如果噪声是零均值高斯白噪声那么取最大和最小会受到噪声极值的影响且采样点越多这种影响越大。所以在实际工作中我更倾向于先对信号进行带通滤波滤除带外噪声再计算峰峰值。但滤波本身会引入振铃可能改变信号的极值这一点在后面的章节会详细讨论。2.3 与最大值、均方根等幅值特征对比如果需要向同事解释为什么要用峰峰值最直接的办法是把几种常用的幅值特征放在一起对比。下表列出了常见特征的定义、适用场景和局限性特征定义对直流偏移敏感性抗噪性光声适用性最大值(\max_{n \in W} s[n])敏感差只反映正向半波容易受基线漂移影响最小值(\min_{n \in W} s[n])敏感差只反映负向半波很少单独使用峰峰值(\max - \min)不敏感中等完整反映振荡幅度适合双极性信号均值(\frac{1}{N}\sum s[n])不敏感好正负抵消后接近零无法表征光声强度均方根RMS(\sqrt{\frac{1}{N}\sum s[n]^2})不敏感中等反映能量但计算量大且受直流偏置影响包络希尔伯特(| \mathcal{H}(s) |)不敏感较好反映瞬时幅度但需要额外变换计算量较大在光声成像中信号是一个近似零均值的衰减振荡所以均值几乎等于零没有区分度。最大值只捕捉了正向峰值如果一个吸收体产生的信号正向半波很小而负向半波很大这通常发生在信号相位反转时最大值就会错误地低估信号强度。均方根虽然能反映能量但对直流偏置没有峰峰值那种“自动抵消”的能力——不过光声信号直流偏置通常已被前期去基线处理。因此峰峰值在这种场景下是最简单、最稳健的幅值特征。另外很多超声成像里的“幅值”指的是包络检测后的峰值但那是经过了Hilbert变换或正交解调之后的事情这里我们讨论的是RF域直接计算的峰峰值两者概念不同写代码时别混淆。3. 用MATLAB实现峰峰值成像的完整流程3.1 数据准备光声RF数据的读取与预处理3.1.1 读入二维或三维数据矩阵光声数据常见的存储格式有.mat、.h5、.dat或厂商自定义格式。MATLAB处理这些并不复杂。假设你已经把数据读取到一个变量rf_data中它的大小是[Nt, Nx]或[Nt, Nx, Ny]。作为示例我们先生成一段模拟光声信号以便演示完整流程% 模拟光声RF数据一个类似阻尼振荡的信号 fs 50e6; % 采样率 50 MHz t (0:999)/fs; % 时间向量20微秒 f0 5e6; % 中心频率 5 MHz t0 5e-6; % 信号到达时间 % 双极性衰减振荡 signal sin(2*pi*f0*(t-t0)) .* exp(-abs(t-t0)/2e-6); % 构造两个空间位置信号一个强一个弱 rf_data zeros(1000, 4); rf_data(:,1) signal * 1.0 0.01*randn(1000,1); rf_data(:,2) signal * 0.5 0.01*randn(1000,1); rf_data(:,3) signal * 0.2 0.01*randn(1000,1); rf_data(:,4) 0.01*randn(1000,1); % 只有噪声的区域这里rf_data有1000个时间采样点4个空间位置。真实的系统里每一列就是一个A-line的原始采样数据。注意我们加了高斯噪声模拟真实采集中的电子噪声。读取真实数据时如果用的是.h5可以用h5read(file.h5, /dataset)如果是.dat且知道采样格式用fread加reshape。这里不展开但有一个通用原则先确认数据的采样率和时间轴长度因为它决定了后续频率滤波的参数。另外多通道数据可能按通道交织存储读取后要重排成[通道数, 采样点数, 帧数]再转置成我们习惯的[采样点数, 通道数]。3.1.2 去基线与带通滤波的取舍原始RF信号往往带有直流偏置这来自换能器放大电路或模拟前端的偏压。直流偏置如果不去除在计算峰峰值时会被max和min共同抵消掉因为峰峰值是差值所以实际上直流偏置对峰峰值没有影响。这是峰峰值的一个天然优势因此在做峰峰值成像时我一般不做专门去基线。带通滤波则是一个需要认真思考的步骤。光声信号的频谱集中在换能器带宽内典型范围从几百kHz到几十MHz。带外噪声特别是低频漂移和高频白噪声会影响max和min的取值。建议先做一个带通滤波但要用零相位滤波以避免相位失真导致峰值偏移% 带通滤波参数 f_low 1e6; % 低频截止 1 MHz f_high 15e6; % 高频截止 15 MHz % 设计巴特沃斯滤波器 [b, a] butter(4, [f_low, f_high]/(fs/2), bandpass); % 零相位滤波 rf_filtered filtfilt(b, a, rf_data);这里butter的阶数为4通带范围1~15 MHz。为什么用4阶而不是更高阶数越高频带边缘越陡峭但群延迟越明显filtfilt可以减轻相位失真但高阶滤波器对瞬态信号容易产生振铃。光声信号本身是宽带瞬态过高阶数的滤波器会在信号前后沿处产生明显振荡反而干扰峰峰值的计算。我的经验是2~4阶巴特沃斯就足够最常用的反而是2阶。如果不想滤波也可以直接计算峰峰值然后在成像前对峰峰值结果做平滑。但那样噪声尖峰的影响可能会残留在图像中。所以建议顺序是滤波 - 计算峰峰值。注意滤波后信号的正负峰值幅度会略微改变但不影响相对对比度。3.2 峰峰值计算的核心操作3.2.1 沿时间轴提取每个A-line的峰峰值有了预处理后的rf_filtered计算峰峰值非常简单。对于二维数据pp_image_1d max(rf_filtered, [], 1) - min(rf_filtered, [], 1);得到的pp_image_1d是一个长度为通道数的向量每个值就是该A-line的峰峰值。对于三维数据例如 B-scan 或多个角度采集的体积数据rf_data大小是[Nt, Nx, Ny]我们需要沿第一个维度计算保留另外两个维度pp_image_2d max(rf_filtered, [], 1) - min(rf_filtered, [], 1); % 结果维度为 [1, Nx, Ny]squeeze后变为 [Nx, Ny] pp_image_2d squeeze(pp_image_2d);这里的max(x, [], 1)在本版本MATLAB中沿第一维返回1×Nx×Ny数组。如果不写squeeze后面看图像时会有一个单一维度不易操作。这是新手常踩的坑。3.2.2 参数说明窗口选择、噪声阈值与输出映射计算峰峰值时不一定总要用整个时间长度。有时信号只出现在某个深度范围而其他时间段都是纯噪声。如果整段计算噪声会导致每个位置都有一个基础峰峰值降低了对比度。此时应该选择信号所在的时间窗口。最常见的做法是用激光触发延迟和换能器焦距来估计时间窗口。例如光声信号在 t5μs 到达持续约5μs那么可以只取 t3μs 到 t10μs 之间的采样点win_start 3e-6; % 窗口起始时间秒 win_end 10e-6; % 窗口结束时间 idx_start round(win_start * fs) 1; idx_end round(win_end * fs); windowed_segment rf_filtered(idx_start:idx_end, :, :); pp_windowed squeeze(max(windowed_segment, [], 1) - min(windowed_segment, [], 1));之后可设置一个噪声阈值把峰峰值低于某个值的像素置为背景值通常这个阈值取无信号区域的峰峰值平均加若干倍标准差。也可以直接用 noise_level作为掩膜。噪声水平可以通过测量采集一段无激光时的数据得到或者取每个A-line前几个采样点还没有接收到光声信号的峰峰值来估计。成像映射这一步对于二维B-scan直接imagesc(pp_image)对于一维线扫描用plot(pp_image)或image逐列填充。但要注意实际坐标轴映射还需要知道换能器的扫描位置与时间-深度转换关系。基本关系是深度 (d c \cdot t / 2)声速 (c) 在生物组织中约为1540 m/s所以时间轴可以换算为深度轴。如果imagesc时手动指定坐标depth t * 1540 / 2 * 1e3; % 单位 mm imagesc(pp_image_1d); colormap(hot); xlabel(通道位置);这里只是示意真实数据还需要根据采集几何定义x轴和y轴。3.3 用MATLAB把峰峰值结果显示为灰度图现在我们把上面的步骤整合到一个完整的、可运行的示例中并添加必要的注释% 完整峰峰值成像示例模拟数据 fs 50e6; t (0:999)/fs; % 生成两个不同强度的光声信号加上噪声 s sin(2*pi*5e6*(t-5e-6)) .* exp(-abs(t-5e-6)/1.5e-6); rf_data [s; 0.6*s; 0.3*s; zeros(size(s))] 0.02 * randn(1000,4); % 滤波 [b, a] butter(4, [1e6 15e6]/(fs/2), bandpass); rf_f filtfilt(b, a, rf_data); % 窗口选择2~12 us w_start 2e-6; w_end 12e-6; [~, i1] min(abs(t - w_start)); [~, i2] min(abs(t - w_end)); rf_w rf_f(i1:i2, :); % 峰峰值 pp max(rf_w, [], 1) - min(rf_w, [], 1); % 显示 figure; subplot(2,1,1); imagesc(rf_w); title(滤波后RF数据每列一个通道); subplot(2,1,2); plot(pp, o-); xlabel(通道序号); ylabel(峰峰值); grid on;代码中[~, i1] min(abs(t - w_start))是一种常见的按时间找索引方法求时间向量与目标值差的绝对值最小值所在位置就是最近的索引。这种写法比round(w_start*fs)1更直观尤其在时间向量有非均匀间隔时。运行后可以看到前三个通道的峰峰值明显高于只有噪声的通道而且近似呈1:0.6:0.3的关系符合我们设置信号的幅度比例。这说明峰峰值成像能定量反应光声信号幅值。如果数据是三维的只需要把rf_f设为[Nt, Nx, Ny]然后max和min沿第一维计算最后squeeze得到二维图像。后续显示时用imagesc即可。4. 峰峰值成像参数调优窗口、阈值与对数压缩4.1 时域窗口长度对峰峰值的影响上面提到窗口的选择直接决定了参与计算的信号长度。窗口太短可能只截取到半个振荡周期导致max和min无法捕捉完整的正负峰峰峰值偏小。窗口太长则会引入更多的噪声采样点噪声极值变大峰峰值升高。那么如何定量选择窗口最稳妥的方法是基于已知信号的带宽光声信号的振荡周期约为中心频率的倒数。以5 MHz中心频率为例一个周期是200 ns。5个周期的持续时间约为1μs。如果系统采集的信号长度有5μs那么信号振荡大约25个周期。为了捕捉完整的正负峰值窗口至少需要覆盖信号的主要能量段通常取2~3倍信号持续时间。但实际中我们总希望窗口尽量窄以提高轴向分辨率。一个折中做法是先观察几个A-line的波形手动估计信号从开始到衰减至10%的时间长度然后取这个时间窗口。窗口起点也很关键。如果起点在信号到达之前那么这段纯噪声会被纳入计算每个位置的峰峰值都会被噪声抬高。我一般先取整段信号计算一个粗略峰峰值图像找到信号覆盖的时间区间然后再用这个区间重新计算。这个过程可以做成半自动用包络检波Hilbert变换找到信号能量超过噪声阈值的区域然后用该区域作为窗口。4.2 噪声阈值与动态范围峰峰值成像的后处理中阈值设定的意义不仅在于美观更在于避免低幅值区域被噪声主导。光声信号幅度差异较大例如强吸收体血管与弱吸收体脂肪可能相差两个数量级。如果直接显示线性灰度图低幅值区域的噪声会非常明显。常见的做法是计算每个A-line在无信号时间段的峰峰值统计其均值和标准差得到噪声水平 (N_0)。设阈值为 (T N_0 k \cdot \sigma_N)k 通常取3~5。将小于 (T) 的峰峰值置为 (T) 或直接置为背景值。在MATLAB中实现如下% 假设已经得到峰峰值矩阵 pp_image (MxN) % 估计噪声取每个A-line最前面100个采样点 noise_seg rf_data(1:100, :); noise_pp max(noise_seg, [], 1) - min(noise_seg, [], 1); threshold mean(noise_pp) 3 * std(noise_pp); pp_thresholded pp_image; pp_thresholded(pp_thresholded threshold) 0;注意这种阈值处理是对每个A-line独立估计的可以有效应对通道间噪声水平不一致的情况。对于三维数据可以先把噪声段改成一个三维子块然后对时间维计算峰峰值得到一张二维噪声图再扩展为整个体积使用的阈值分布。阈值设得太高小吸收体会丢失太低则噪声背景明显。我一般在调试时使用滑块交互式调整figure; im imagesc(pp_image); thr mean(noise_pp) 3*std(noise_pp); caxis([thr max(pp_image(:))]);这虽然不改变数据但通过调节显示范围caxis可以快速评估阈值。后续再将同样的裁剪逻辑应用到最终保存的数据。4.3 对数压缩与动态范围显示光声图像往往需要显示很大的灰度范围。峰峰值从几十到几万直接imagesc会压缩低幅值信息。通常采用对数压缩将数据映射到人眼可感知的动态范围。MATLAB中可以直接用log(pp_image 1)或db函数。其中加1是为了避免log(0)产生无穷大。常见的动态范围压缩公式是[ I_{\text{display}} \frac{\log_{10}(1 \alpha \cdot I)}{\log_{10}(1 \alpha \cdot I_{\max})} ]其中 (\alpha) 是压缩系数控制压缩强度。(\alpha) 越大低幅值区域被分配更多灰度级。实现为alpha 100; I_disp log10(1 alpha * (pp_image / max(pp_image(:)))); imagesc(I_disp);也可以使用商用的显示方案先去掉底噪再取对数。很多人会把这个处理直接放进行业成熟的重建脚本里。但要注意对数压缩会改变图像之间的相对比例如果后续要做定量分析如血氧饱和度应该在原始峰峰值图像上进行而不是压缩后的显示数据。4.4 性能优化向量化与并行计算光声数据集往往很大例如三维体积数据有1000个时间点 × 256 × 256 的空间点数据量约6550万个采样点用双层循环计算峰峰值会非常慢。上面的max(..., [], 1)方法已经向量化但对三维数据需要小心处理。一个常见的错误是直接使用max(rf, [], 1)结果维度变成[1, Nx, Ny]然后错误地使用max(result, [], 2)等导致计算混乱。正确做法见前面代码。如果体积数据太大内存无法同时容纳可以采用分块处理。例如沿沿x方向每次读取一块256×16×Nt的子体积计算峰峰值后再拼接。MATLAB的tall数组和imageDatastore也能处理但更直接的还是for循环配合matfile对象。此外如果有Parallel Computing Toolbox可以对通道维使用parfor% 将数据按通道分给不同并行worker pp_all zeros(Nx, Ny); parfor i 1:Nx pp_all(i, :) max(rf_time(:, i, :), [], 1) - min(rf_time(:, i, :), [], 1); end注意这里rf_time是一个[Nt, Nx, Ny]的三维数组parfor切片时MATLAB会自动传输每个i需要的数据。不过并行开销大只有当数据维度较长时使用才划算。对于单次二维图像向量化操作已经足够快不需要并行。5. 验证与进阶三维峰峰值投影、多波长比值与滤波顺序技巧5.1 三维体积数据上的峰峰值投影当你有三维光声数据B-scan 堆叠成体积时除了生成每一层的二维图像往往还需要做一次最大强度投影MIP或平均强度投影展示整体结构。峰峰值成像可以在三维上直接计算体积数据得到一个三维幅值矩阵然后再沿深度方向做投影。假设rf_volume维度是[Nt, Nx, Ny]先求峰峰值体数据pp_volume max(rf_volume, [], 1) - min(rf_volume, [], 1); pp_volume squeeze(pp_volume); % 变成 Nx×Ny这样每个(x,y)位置都是一个深度上最大振荡幅度的投影。但这个投影丢失了深度信息。更常用的是最大强度投影沿深度轴时间轴做而不是先算峰峰值。两种做法目的不同如果你想要的是“这个位置不管多深都存在强吸收体”应该对原始RF信号取深度方向的峰峰值如果你想要的是“某个深度切面”应该对每个时间片分别做峰峰值。前者把三维体积压缩成二维地图适合观察整体血管形态后者用于逐层分析。MATLAB中实现最大强度投影的方式是mip_image max(pp_volume, [], 3); % 如果 pp_volume 是 Nx×Ny×Nz如果你的原始三维数据维度是[Nt, Nx, Ny]那么pp_volume恰好是[Nx, Ny]无法再沿深度投影。这种情况下你需要把体积按时间窗口分割成多个子体积每个子体积计算峰峰值得到多幅二维图像然后再沿深度轴即子体积序号取最大。这个过程中窗口的划分会影响投影结果的轴向分辨率。小窗口带来高分辨率但噪声变大大窗口则相反。实践中我通常将时间窗口设为信号振荡周期的3~5倍。5.2 多波长光声中的峰峰值比值光声成像的一个重要应用是多波长成像通过不同波长下吸收体的吸收系数差异来计算血氧饱和度等参数。峰峰值成像在这种场景下相当方便因为峰峰值与光声信号幅值成正比而信号幅值又正比于吸收系数和光通量。如果我们采集了两个波长 (\lambda_1) 和 (\lambda_2) 的峰峰值图像 (PP_1) 和 (PP_2)那么比值 (R PP_1 / PP_2) 可以抵消掉光通量不均和表面衰减的影响保留与血氧饱和度相关的吸收比信息。计算时要注意分母不能为零而且要使用同一时间窗口。MATLAB代码可以写成PP1 max(data_lambda1, [], 1) - min(data_lambda1, [], 1); PP2 max(data_lambda2, [], 1) - min(data_lambda2, [], 1); ratio (PP1 eps) ./ (PP2 eps);其中eps是为了防止分母为零。比值图像一般不需要再做大动态范围压缩直接显示即可。不过多波长成像更精细的做法是先用上述峰峰值方法得到强度图像再在空间上做一个中值滤波平滑掉噪声然后再算比值。经过平滑之后比值图会稳定很多。5.3 一个被忽略的细节滤波顺序对峰峰值的影响最后分享一个我踩过的坑。很多人会先对RF信号做带通滤波再计算峰峰值。这个顺序很自然但滤波器的振铃效应可能让信号产生额外的正负极值特别是在信号边缘处。filtfilt虽说是零相位但零相位不代表零振铃。对于短促的光声脉冲一个通带截止频率接近中心频率的滤波器会在脉冲前后形成对称的过冲这个过冲幅度有时甚至比原始信号峰值还大。解决的办法是先计算峰峰值然后再对峰峰值图像进行空间域平滑。或者如果一定要在时间域滤波那就选用 Butterworth 2阶且通带带宽足够宽的滤波器并且在计算峰峰值之前截去信号两端各几微秒的数据避开振铃的起始段。我个人的策略是原始RF信号先只做去除直流减均值用max和min计算峰峰值然后对峰峰值图像做3×3中值滤波和低通空间滤波。这样既保留了峰峰值对双极性振荡的真实度量又避免了滤波振铃的干扰。在实际的大光声数据集上这个策略对比度稳定、伪影少。这也是我在处理光声峰峰值成像时最常推荐的做法。本文还有配套的精品资源点击获取
光声峰峰值成像：MATLAB实现与参数调优指南
简介光声峰峰值成像利用光吸收产生的超声信号重建组织内部光吸收分布在生物医学光学成像与病变识别中具有实用价值。面向光声成像研究者、生物医学工程相关专业学生以及需要快速构建成像算法的开发者这份MATLAB资源提供了一套完整的峰峰值幅值成像代码包。压缩包内共有4个m文件整体仅3KB分别对应信号预处理与滤波、主流程控制、三维曲面可视化以及数据翻转平滑等模块代码精简、分工明确便于直接阅读和二次修改。已有433人学习下载。通过研读和运行这些脚本读者能够掌握光声信号去噪、峰峰值特征提取、三维图像重建等关键操作理解从原始声波信号到光吸收分布图的完整处理链条还可根据自身数据调整参数以优化成像质量适合作为光声成像入门与算法验证的参考资料。1. 光声峰峰值成像为什么用峰峰值而不是最大值或均值做光声成像的同学应该遇到过这样一个选择从采集到的射频RF信号重建图像时到底取信号的哪个特征来代表这个像素有人直接取最大值有人计算均值还有人在做“幅值成像”时下意识地用了信号的峰峰值。区别在于光声信号本质上是激光激发产生的超声波瞬态波它通常是双极性振荡包含着正负极值。直接取最大值只关注了正向半波漏掉了负向半波的信息而取均值几乎不能反映信号的瞬态强度。峰峰值Peak-to-Peak, PP定义为信号在某一时间窗口内最大值与最小值之差。在光声峰峰值成像中这个值表达了该空间位置处光声信号的绝对振荡幅度对直流偏移不敏感也天然抑制了缓慢变化的背景漂移。相比最大值、均值或均方根峰峰值能在不依赖绝对基线的情况下刻画吸收体的浓度和尺寸而且计算简单、无需拟合。这篇文章就围绕这个主题展开从信号模型到MATLAB实现再到参数调优和三维扩展。适合正在写光声信号处理脚本、做图像重建或调试采集系统的工程师阅读哪怕是刚接触MATLAB光声数据分析的人也可以照着例子跑通。2. 光声信号模型与峰峰值特征提取原理2.1 光声信号的产生与探测光声成像的原理并不复杂短脉冲激光照射吸收体如血管中的血红蛋白吸收体快速热膨胀产生压力波这个压力波被超声换能器接收形成时域信号我们称之为一个A-line。典型的光声信号类似于带阻尼的振荡正半周和负半周幅度可能不对称但整体包络反映了吸收体的光吸收分布。由于换能器带宽有限信号并不是理想的狄拉克脉冲而是一个有限频宽的瞬态波。这种信号的特征参数很多。对于一个采集到的RF数据矩阵行代表时间采样列代表空间位置一维扫描时是一个A-line索引二维B-scan时是扫描角度或x位置三维体积数据则是x和y两个方向的排列。光声峰峰值成像要做的就是沿着时间轴每一列计算信号的峰值到峰值的幅度再把这个幅度映射到对应的空间坐标上形成强度图像。2.2 峰峰值Peak-to-Peak的定义与数学表达设光声信号为 ( s(t) )其采样序列为 ( s[n], n1,\dots,N )。对应一个空间点的峰峰值定义为[ PP \max_{n \in W} s[n] - \min_{n \in W} s[n] ]其中 ( W ) 是所选的时间窗口可以是整个A-line长度也可以是激光发射后的一个特定时间区间。窗口的选择直接决定了成像的深度范围窗口起点对应探测起始深度窗口长度对应成像深度厚度。在MATLAB中这个计算可以用最朴素的循环完成但更好的做法是使用max和min函数在指定维度上作用。对于二维矩阵data(Nt, Nx)一行代码就能完成pp max(data, [], 1) - min(data, [], 1);这里的max(data, [], 1)是沿第一个维度时间维度求最大值返回一个1×Nx的行向量。min同理。两个向量相减得到每个A-line的峰峰值。这个操作没有显式循环效率很高。不过要注意直接对原始RF信号做max和min有一个陷阱信号中的随机噪声会产生额外的尖峰。如果噪声是零均值高斯白噪声那么取最大和最小会受到噪声极值的影响且采样点越多这种影响越大。所以在实际工作中我更倾向于先对信号进行带通滤波滤除带外噪声再计算峰峰值。但滤波本身会引入振铃可能改变信号的极值这一点在后面的章节会详细讨论。2.3 与最大值、均方根等幅值特征对比如果需要向同事解释为什么要用峰峰值最直接的办法是把几种常用的幅值特征放在一起对比。下表列出了常见特征的定义、适用场景和局限性特征定义对直流偏移敏感性抗噪性光声适用性最大值(\max_{n \in W} s[n])敏感差只反映正向半波容易受基线漂移影响最小值(\min_{n \in W} s[n])敏感差只反映负向半波很少单独使用峰峰值(\max - \min)不敏感中等完整反映振荡幅度适合双极性信号均值(\frac{1}{N}\sum s[n])不敏感好正负抵消后接近零无法表征光声强度均方根RMS(\sqrt{\frac{1}{N}\sum s[n]^2})不敏感中等反映能量但计算量大且受直流偏置影响包络希尔伯特(| \mathcal{H}(s) |)不敏感较好反映瞬时幅度但需要额外变换计算量较大在光声成像中信号是一个近似零均值的衰减振荡所以均值几乎等于零没有区分度。最大值只捕捉了正向峰值如果一个吸收体产生的信号正向半波很小而负向半波很大这通常发生在信号相位反转时最大值就会错误地低估信号强度。均方根虽然能反映能量但对直流偏置没有峰峰值那种“自动抵消”的能力——不过光声信号直流偏置通常已被前期去基线处理。因此峰峰值在这种场景下是最简单、最稳健的幅值特征。另外很多超声成像里的“幅值”指的是包络检测后的峰值但那是经过了Hilbert变换或正交解调之后的事情这里我们讨论的是RF域直接计算的峰峰值两者概念不同写代码时别混淆。3. 用MATLAB实现峰峰值成像的完整流程3.1 数据准备光声RF数据的读取与预处理3.1.1 读入二维或三维数据矩阵光声数据常见的存储格式有.mat、.h5、.dat或厂商自定义格式。MATLAB处理这些并不复杂。假设你已经把数据读取到一个变量rf_data中它的大小是[Nt, Nx]或[Nt, Nx, Ny]。作为示例我们先生成一段模拟光声信号以便演示完整流程% 模拟光声RF数据一个类似阻尼振荡的信号 fs 50e6; % 采样率 50 MHz t (0:999)/fs; % 时间向量20微秒 f0 5e6; % 中心频率 5 MHz t0 5e-6; % 信号到达时间 % 双极性衰减振荡 signal sin(2*pi*f0*(t-t0)) .* exp(-abs(t-t0)/2e-6); % 构造两个空间位置信号一个强一个弱 rf_data zeros(1000, 4); rf_data(:,1) signal * 1.0 0.01*randn(1000,1); rf_data(:,2) signal * 0.5 0.01*randn(1000,1); rf_data(:,3) signal * 0.2 0.01*randn(1000,1); rf_data(:,4) 0.01*randn(1000,1); % 只有噪声的区域这里rf_data有1000个时间采样点4个空间位置。真实的系统里每一列就是一个A-line的原始采样数据。注意我们加了高斯噪声模拟真实采集中的电子噪声。读取真实数据时如果用的是.h5可以用h5read(file.h5, /dataset)如果是.dat且知道采样格式用fread加reshape。这里不展开但有一个通用原则先确认数据的采样率和时间轴长度因为它决定了后续频率滤波的参数。另外多通道数据可能按通道交织存储读取后要重排成[通道数, 采样点数, 帧数]再转置成我们习惯的[采样点数, 通道数]。3.1.2 去基线与带通滤波的取舍原始RF信号往往带有直流偏置这来自换能器放大电路或模拟前端的偏压。直流偏置如果不去除在计算峰峰值时会被max和min共同抵消掉因为峰峰值是差值所以实际上直流偏置对峰峰值没有影响。这是峰峰值的一个天然优势因此在做峰峰值成像时我一般不做专门去基线。带通滤波则是一个需要认真思考的步骤。光声信号的频谱集中在换能器带宽内典型范围从几百kHz到几十MHz。带外噪声特别是低频漂移和高频白噪声会影响max和min的取值。建议先做一个带通滤波但要用零相位滤波以避免相位失真导致峰值偏移% 带通滤波参数 f_low 1e6; % 低频截止 1 MHz f_high 15e6; % 高频截止 15 MHz % 设计巴特沃斯滤波器 [b, a] butter(4, [f_low, f_high]/(fs/2), bandpass); % 零相位滤波 rf_filtered filtfilt(b, a, rf_data);这里butter的阶数为4通带范围1~15 MHz。为什么用4阶而不是更高阶数越高频带边缘越陡峭但群延迟越明显filtfilt可以减轻相位失真但高阶滤波器对瞬态信号容易产生振铃。光声信号本身是宽带瞬态过高阶数的滤波器会在信号前后沿处产生明显振荡反而干扰峰峰值的计算。我的经验是2~4阶巴特沃斯就足够最常用的反而是2阶。如果不想滤波也可以直接计算峰峰值然后在成像前对峰峰值结果做平滑。但那样噪声尖峰的影响可能会残留在图像中。所以建议顺序是滤波 - 计算峰峰值。注意滤波后信号的正负峰值幅度会略微改变但不影响相对对比度。3.2 峰峰值计算的核心操作3.2.1 沿时间轴提取每个A-line的峰峰值有了预处理后的rf_filtered计算峰峰值非常简单。对于二维数据pp_image_1d max(rf_filtered, [], 1) - min(rf_filtered, [], 1);得到的pp_image_1d是一个长度为通道数的向量每个值就是该A-line的峰峰值。对于三维数据例如 B-scan 或多个角度采集的体积数据rf_data大小是[Nt, Nx, Ny]我们需要沿第一个维度计算保留另外两个维度pp_image_2d max(rf_filtered, [], 1) - min(rf_filtered, [], 1); % 结果维度为 [1, Nx, Ny]squeeze后变为 [Nx, Ny] pp_image_2d squeeze(pp_image_2d);这里的max(x, [], 1)在本版本MATLAB中沿第一维返回1×Nx×Ny数组。如果不写squeeze后面看图像时会有一个单一维度不易操作。这是新手常踩的坑。3.2.2 参数说明窗口选择、噪声阈值与输出映射计算峰峰值时不一定总要用整个时间长度。有时信号只出现在某个深度范围而其他时间段都是纯噪声。如果整段计算噪声会导致每个位置都有一个基础峰峰值降低了对比度。此时应该选择信号所在的时间窗口。最常见的做法是用激光触发延迟和换能器焦距来估计时间窗口。例如光声信号在 t5μs 到达持续约5μs那么可以只取 t3μs 到 t10μs 之间的采样点win_start 3e-6; % 窗口起始时间秒 win_end 10e-6; % 窗口结束时间 idx_start round(win_start * fs) 1; idx_end round(win_end * fs); windowed_segment rf_filtered(idx_start:idx_end, :, :); pp_windowed squeeze(max(windowed_segment, [], 1) - min(windowed_segment, [], 1));之后可设置一个噪声阈值把峰峰值低于某个值的像素置为背景值通常这个阈值取无信号区域的峰峰值平均加若干倍标准差。也可以直接用 noise_level作为掩膜。噪声水平可以通过测量采集一段无激光时的数据得到或者取每个A-line前几个采样点还没有接收到光声信号的峰峰值来估计。成像映射这一步对于二维B-scan直接imagesc(pp_image)对于一维线扫描用plot(pp_image)或image逐列填充。但要注意实际坐标轴映射还需要知道换能器的扫描位置与时间-深度转换关系。基本关系是深度 (d c \cdot t / 2)声速 (c) 在生物组织中约为1540 m/s所以时间轴可以换算为深度轴。如果imagesc时手动指定坐标depth t * 1540 / 2 * 1e3; % 单位 mm imagesc(pp_image_1d); colormap(hot); xlabel(通道位置);这里只是示意真实数据还需要根据采集几何定义x轴和y轴。3.3 用MATLAB把峰峰值结果显示为灰度图现在我们把上面的步骤整合到一个完整的、可运行的示例中并添加必要的注释% 完整峰峰值成像示例模拟数据 fs 50e6; t (0:999)/fs; % 生成两个不同强度的光声信号加上噪声 s sin(2*pi*5e6*(t-5e-6)) .* exp(-abs(t-5e-6)/1.5e-6); rf_data [s; 0.6*s; 0.3*s; zeros(size(s))] 0.02 * randn(1000,4); % 滤波 [b, a] butter(4, [1e6 15e6]/(fs/2), bandpass); rf_f filtfilt(b, a, rf_data); % 窗口选择2~12 us w_start 2e-6; w_end 12e-6; [~, i1] min(abs(t - w_start)); [~, i2] min(abs(t - w_end)); rf_w rf_f(i1:i2, :); % 峰峰值 pp max(rf_w, [], 1) - min(rf_w, [], 1); % 显示 figure; subplot(2,1,1); imagesc(rf_w); title(滤波后RF数据每列一个通道); subplot(2,1,2); plot(pp, o-); xlabel(通道序号); ylabel(峰峰值); grid on;代码中[~, i1] min(abs(t - w_start))是一种常见的按时间找索引方法求时间向量与目标值差的绝对值最小值所在位置就是最近的索引。这种写法比round(w_start*fs)1更直观尤其在时间向量有非均匀间隔时。运行后可以看到前三个通道的峰峰值明显高于只有噪声的通道而且近似呈1:0.6:0.3的关系符合我们设置信号的幅度比例。这说明峰峰值成像能定量反应光声信号幅值。如果数据是三维的只需要把rf_f设为[Nt, Nx, Ny]然后max和min沿第一维计算最后squeeze得到二维图像。后续显示时用imagesc即可。4. 峰峰值成像参数调优窗口、阈值与对数压缩4.1 时域窗口长度对峰峰值的影响上面提到窗口的选择直接决定了参与计算的信号长度。窗口太短可能只截取到半个振荡周期导致max和min无法捕捉完整的正负峰峰峰值偏小。窗口太长则会引入更多的噪声采样点噪声极值变大峰峰值升高。那么如何定量选择窗口最稳妥的方法是基于已知信号的带宽光声信号的振荡周期约为中心频率的倒数。以5 MHz中心频率为例一个周期是200 ns。5个周期的持续时间约为1μs。如果系统采集的信号长度有5μs那么信号振荡大约25个周期。为了捕捉完整的正负峰值窗口至少需要覆盖信号的主要能量段通常取2~3倍信号持续时间。但实际中我们总希望窗口尽量窄以提高轴向分辨率。一个折中做法是先观察几个A-line的波形手动估计信号从开始到衰减至10%的时间长度然后取这个时间窗口。窗口起点也很关键。如果起点在信号到达之前那么这段纯噪声会被纳入计算每个位置的峰峰值都会被噪声抬高。我一般先取整段信号计算一个粗略峰峰值图像找到信号覆盖的时间区间然后再用这个区间重新计算。这个过程可以做成半自动用包络检波Hilbert变换找到信号能量超过噪声阈值的区域然后用该区域作为窗口。4.2 噪声阈值与动态范围峰峰值成像的后处理中阈值设定的意义不仅在于美观更在于避免低幅值区域被噪声主导。光声信号幅度差异较大例如强吸收体血管与弱吸收体脂肪可能相差两个数量级。如果直接显示线性灰度图低幅值区域的噪声会非常明显。常见的做法是计算每个A-line在无信号时间段的峰峰值统计其均值和标准差得到噪声水平 (N_0)。设阈值为 (T N_0 k \cdot \sigma_N)k 通常取3~5。将小于 (T) 的峰峰值置为 (T) 或直接置为背景值。在MATLAB中实现如下% 假设已经得到峰峰值矩阵 pp_image (MxN) % 估计噪声取每个A-line最前面100个采样点 noise_seg rf_data(1:100, :); noise_pp max(noise_seg, [], 1) - min(noise_seg, [], 1); threshold mean(noise_pp) 3 * std(noise_pp); pp_thresholded pp_image; pp_thresholded(pp_thresholded threshold) 0;注意这种阈值处理是对每个A-line独立估计的可以有效应对通道间噪声水平不一致的情况。对于三维数据可以先把噪声段改成一个三维子块然后对时间维计算峰峰值得到一张二维噪声图再扩展为整个体积使用的阈值分布。阈值设得太高小吸收体会丢失太低则噪声背景明显。我一般在调试时使用滑块交互式调整figure; im imagesc(pp_image); thr mean(noise_pp) 3*std(noise_pp); caxis([thr max(pp_image(:))]);这虽然不改变数据但通过调节显示范围caxis可以快速评估阈值。后续再将同样的裁剪逻辑应用到最终保存的数据。4.3 对数压缩与动态范围显示光声图像往往需要显示很大的灰度范围。峰峰值从几十到几万直接imagesc会压缩低幅值信息。通常采用对数压缩将数据映射到人眼可感知的动态范围。MATLAB中可以直接用log(pp_image 1)或db函数。其中加1是为了避免log(0)产生无穷大。常见的动态范围压缩公式是[ I_{\text{display}} \frac{\log_{10}(1 \alpha \cdot I)}{\log_{10}(1 \alpha \cdot I_{\max})} ]其中 (\alpha) 是压缩系数控制压缩强度。(\alpha) 越大低幅值区域被分配更多灰度级。实现为alpha 100; I_disp log10(1 alpha * (pp_image / max(pp_image(:)))); imagesc(I_disp);也可以使用商用的显示方案先去掉底噪再取对数。很多人会把这个处理直接放进行业成熟的重建脚本里。但要注意对数压缩会改变图像之间的相对比例如果后续要做定量分析如血氧饱和度应该在原始峰峰值图像上进行而不是压缩后的显示数据。4.4 性能优化向量化与并行计算光声数据集往往很大例如三维体积数据有1000个时间点 × 256 × 256 的空间点数据量约6550万个采样点用双层循环计算峰峰值会非常慢。上面的max(..., [], 1)方法已经向量化但对三维数据需要小心处理。一个常见的错误是直接使用max(rf, [], 1)结果维度变成[1, Nx, Ny]然后错误地使用max(result, [], 2)等导致计算混乱。正确做法见前面代码。如果体积数据太大内存无法同时容纳可以采用分块处理。例如沿沿x方向每次读取一块256×16×Nt的子体积计算峰峰值后再拼接。MATLAB的tall数组和imageDatastore也能处理但更直接的还是for循环配合matfile对象。此外如果有Parallel Computing Toolbox可以对通道维使用parfor% 将数据按通道分给不同并行worker pp_all zeros(Nx, Ny); parfor i 1:Nx pp_all(i, :) max(rf_time(:, i, :), [], 1) - min(rf_time(:, i, :), [], 1); end注意这里rf_time是一个[Nt, Nx, Ny]的三维数组parfor切片时MATLAB会自动传输每个i需要的数据。不过并行开销大只有当数据维度较长时使用才划算。对于单次二维图像向量化操作已经足够快不需要并行。5. 验证与进阶三维峰峰值投影、多波长比值与滤波顺序技巧5.1 三维体积数据上的峰峰值投影当你有三维光声数据B-scan 堆叠成体积时除了生成每一层的二维图像往往还需要做一次最大强度投影MIP或平均强度投影展示整体结构。峰峰值成像可以在三维上直接计算体积数据得到一个三维幅值矩阵然后再沿深度方向做投影。假设rf_volume维度是[Nt, Nx, Ny]先求峰峰值体数据pp_volume max(rf_volume, [], 1) - min(rf_volume, [], 1); pp_volume squeeze(pp_volume); % 变成 Nx×Ny这样每个(x,y)位置都是一个深度上最大振荡幅度的投影。但这个投影丢失了深度信息。更常用的是最大强度投影沿深度轴时间轴做而不是先算峰峰值。两种做法目的不同如果你想要的是“这个位置不管多深都存在强吸收体”应该对原始RF信号取深度方向的峰峰值如果你想要的是“某个深度切面”应该对每个时间片分别做峰峰值。前者把三维体积压缩成二维地图适合观察整体血管形态后者用于逐层分析。MATLAB中实现最大强度投影的方式是mip_image max(pp_volume, [], 3); % 如果 pp_volume 是 Nx×Ny×Nz如果你的原始三维数据维度是[Nt, Nx, Ny]那么pp_volume恰好是[Nx, Ny]无法再沿深度投影。这种情况下你需要把体积按时间窗口分割成多个子体积每个子体积计算峰峰值得到多幅二维图像然后再沿深度轴即子体积序号取最大。这个过程中窗口的划分会影响投影结果的轴向分辨率。小窗口带来高分辨率但噪声变大大窗口则相反。实践中我通常将时间窗口设为信号振荡周期的3~5倍。5.2 多波长光声中的峰峰值比值光声成像的一个重要应用是多波长成像通过不同波长下吸收体的吸收系数差异来计算血氧饱和度等参数。峰峰值成像在这种场景下相当方便因为峰峰值与光声信号幅值成正比而信号幅值又正比于吸收系数和光通量。如果我们采集了两个波长 (\lambda_1) 和 (\lambda_2) 的峰峰值图像 (PP_1) 和 (PP_2)那么比值 (R PP_1 / PP_2) 可以抵消掉光通量不均和表面衰减的影响保留与血氧饱和度相关的吸收比信息。计算时要注意分母不能为零而且要使用同一时间窗口。MATLAB代码可以写成PP1 max(data_lambda1, [], 1) - min(data_lambda1, [], 1); PP2 max(data_lambda2, [], 1) - min(data_lambda2, [], 1); ratio (PP1 eps) ./ (PP2 eps);其中eps是为了防止分母为零。比值图像一般不需要再做大动态范围压缩直接显示即可。不过多波长成像更精细的做法是先用上述峰峰值方法得到强度图像再在空间上做一个中值滤波平滑掉噪声然后再算比值。经过平滑之后比值图会稳定很多。5.3 一个被忽略的细节滤波顺序对峰峰值的影响最后分享一个我踩过的坑。很多人会先对RF信号做带通滤波再计算峰峰值。这个顺序很自然但滤波器的振铃效应可能让信号产生额外的正负极值特别是在信号边缘处。filtfilt虽说是零相位但零相位不代表零振铃。对于短促的光声脉冲一个通带截止频率接近中心频率的滤波器会在脉冲前后形成对称的过冲这个过冲幅度有时甚至比原始信号峰值还大。解决的办法是先计算峰峰值然后再对峰峰值图像进行空间域平滑。或者如果一定要在时间域滤波那就选用 Butterworth 2阶且通带带宽足够宽的滤波器并且在计算峰峰值之前截去信号两端各几微秒的数据避开振铃的起始段。我个人的策略是原始RF信号先只做去除直流减均值用max和min计算峰峰值然后对峰峰值图像做3×3中值滤波和低通空间滤波。这样既保留了峰峰值对双极性振荡的真实度量又避免了滤波振铃的干扰。在实际的大光声数据集上这个策略对比度稳定、伪影少。这也是我在处理光声峰峰值成像时最常推荐的做法。本文还有配套的精品资源点击获取