Field II 超声相控阵仿真系列：多角度平面波相干合成提升成像质量

1. 从“快”到“好”为什么单次平面波成像不够用大家好我是老张在超声成像仿真这个领域摸爬滚打了十来年用过不少工具Field II算是我的老朋友了。今天咱们不聊那些复杂的理论推导就说说一个非常实际的问题怎么让超声图像又快又好很多刚接触超声相控阵仿真的朋友可能都听说过或者尝试过平面波成像。它的魅力太大了——一次发射整个成像区域的回波数据就都拿到了帧频可以轻松飙到几千甚至上万帧每秒。这速度用来观察心脏跳动、血液流动简直是“神器”。我刚开始用的时候也特别兴奋觉得这技术太酷了。但很快现实就给我上了一课。当我兴冲冲地把仿真出来的图像和传统的聚焦扫描图像一对比心就凉了半截图像怎么这么“糊”目标的边缘不清晰对比度也差深一点的目标都快和背景噪声混在一起了。这其实就是单次平面波成像的“阿喀琉斯之踵”。你可以把它想象成用手电筒照明。传统的聚焦发射就像把手电筒的光调成一个很细、很亮的光束集中照在一个点上这个点当然被照得清清楚楚。而单次平面波发射相当于把手电筒的聚光罩拿掉让光均匀地洒向整个房间。整个房间是同时被照亮了成像速度快但每个地方的光线都变弱了也不够集中导致你看清书本上的小字高分辨率和分辨阴影里的细节高对比度就变得困难。在声学上这是因为声波能量没有在特定深度聚焦导致横向分辨率和对比度噪声比显著下降。对于要求精细诊断的场景比如要看清细微的血管壁、微小的钙化点这样的图像质量是远远不够的。那么有没有办法既能保留平面波“快”的优点又能把“好”给找回来呢当然有这就是我们今天要深入探讨的“复合平面波成像”的核心思路用多个角度的“光”去照亮同一个场景然后把它们巧妙地合成一张更清晰的照片。简单来说我们不满足于只从正面“看”一次我们还从左边一点点角度、右边一点点角度各“看”几次。每一次虽然单独看都不够清晰但当我们把所有这些从不同角度获取的信息按照严格的物理规律也就是相干合成叠加在一起时那些在所有角度里都稳定存在的真实信号会被增强而那些随机的噪声和伪影则会被抑制。这就好比多个人从不同角度描述同一个物体综合他们的描述你就能在脑海里构建出一个远比单人描述更准确、更立体的形象。2. 仿真基石如何用Field II搭建你的相控阵模型工欲善其事必先利其器。在开始玩转多角度合成之前我们得先把“舞台”——也就是超声相控阵探头模型——在Field II里搭建好。这一步的参数设置是后续所有仿真的基础一点都不能马虎。我见过不少新手在这里栽跟头仿真结果不对折腾半天才发现是探头参数设错了。Field II是一个基于声场计算的MATLAB工具箱它的强大之处在于可以非常物理地模拟超声波的发射、传播和接收过程。我们首先要定义一个符合实际的换能器阵列。下面我就以一个临床上常用的P7-4相控阵探头为例带大家一步步走通。别担心我会把每个参数代表什么、为什么这么设都讲清楚。%% 1. 换能器核心参数设置 trans.nameP7-4; % 探头型号方便自己记录 trans.fc6e6; % 中心频率6 MHz这是该探头的典型工作频率 trans.numele 64; % 阵元总数64个阵元 trans.width 136.9e-6; % 单个阵元宽度单位米约137微米 trans.pitch 171.1e-6; % 阵元中心间距节距单位米约171微米 trans.kerf trans.pitch - trans.width; % 阵元间的间隙kerf trans.heigh 14e-3; % 阵元高度Elevation方向单位米14毫米 trans.elevationFocus 60e-3; % 高程方向聚焦深度单位米60毫米 trans.focus [0 0 100e6]/1000; % 初始定义在平面波成像中通常不直接使用 trans.ElementPos trans.pitch * (-((trans.numele-1)/2):((trans.numele-1)/2)); % 计算每个阵元的中心位置坐标 trans.c 1540; % 生物软组织中的平均声速单位米/秒 xT trans.ElementPos; % 将阵元位置坐标赋值给变量xT后续计算会频繁用到这里有几个关键点需要拎出来说说。trans.pitch这个参数至关重要它直接决定了我们探头的最大可用视角和避免栅瓣的条件。trans.kerf是阵元间的缝隙虽然小但在计算声场时是有影响的。trans.ElementPos这行代码是在生成一个包含64个元素的数组每个元素值代表对应阵元在X轴上的位置。比如中间两个阵元的位置接近0两端的阵元位置则分别是正负最大值。这个坐标数组是我们后续计算发射和接收延时的基础。定义好几何和物理参数后我们需要在Field II中创建发射和接收孔径。你可以把它们理解成探头的工作模式开关。%% 2. 创建发射与接收孔径 emit xdc_linear_array(trans.numele, trans.width, trans.heigh, trans.kerf, 3, 10, trans.focus); rcv xdc_linear_array(trans.numele, trans.width, trans.heigh, trans.kerf, 3, 10, trans.focus);xdc_linear_array函数创建了一个线性阵列模型。里面的参数3和10是Field II内部用于控制声场计算精度的参数通常保持默认即可。接下来我们要给这个探头定义它的“声音特质”包括它发出什么样的脉冲以及它如何被电信号激励。%% 3. 设置脉冲响应与激励信号 impulse sin(2*pi*trans.fc * (0:1/userset.fs:2/trans.fc)); impulse_response impulse .* hanning(max(size(impulse))); xdc_impulse(emit, impulse_response); % 设置发射孔径的脉冲响应 xdc_impulse(rcv, impulse_response); % 设置接收孔径的脉冲响应通常与发射相同 excitation sin(2*pi*trans.fc * (0:1/userset.fs:2/trans.fc)); % 生成激励信号 xdc_excitation(emit, excitation); % 将激励信号加载到发射孔径上这里impulse_response是探头的固有特性可以理解为它受到一个理想冲击时产生的振动波形。我们用一个中心频率为trans.fc的正弦波加汉宁窗来模拟这样更接近实际探头的带通特性。excitation则是我们实际加在阵元上的电压信号。注意在平面波发射模式下我们不会在这里设置聚焦延时聚焦或者说偏转的延时会在下一步单独施加。最后别忘了设置整个系统的采样频率这决定了回波信号数字化的精细程度。set_sampling(userset.fs); % 例如 userset.fs 100e6; 即100 MHz采样率3. 发射的艺术计算多角度平面波的发射延时模型搭好了现在让我们来指挥这个探头“工作”。单次平面波发射就是让所有阵元在同一时刻被激发产生一个垂直于探头表面的平面波前。但为了实现复合平面波成像我们需要这个波前能“歪着”发射出去也就是进行角度偏转。这个“歪”是怎么实现的呢靠的就是给不同阵元施加微小的时间差即发射延时。原理很像小时候玩过的打水漂如果你同时扔出一排石子它们会同时入水但如果你从左到右依次延迟扔出这一排石子入水点就会连成一条斜线相当于水波前进的方向发生了偏转。对于我们的线性阵列计算这个延时公式非常简洁tx_delay pitch * [0:(N-1)] * sin(steer_angle) / c其中pitch是阵元间距[0:(N-1)]是阵元索引假设中心为参考点steer_angle是我们想要的偏转角度c是声速。这里有个非常重要的符号约定在Field II和多数仿真惯例中通常定义角度为正时波束向右偏转角度为负时向左偏转。这个公式的物理意义是为了形成一个朝向steer_angle方向的平面波位置靠右的阵元需要稍微提前发射以“追赶”上左侧阵元发出的波前。在代码中我们可以将其封装成一个函数方便多次调用function delays plane_wave_delayt(trans, steer_angle) % 计算平面波偏转发射的延时 % trans: 换能器参数结构体 % steer_angle: 偏转角度弧度制 element_index 0:(trans.numele - 1); delays trans.pitch * element_index * sin(steer_angle) / trans.c; % 通常我们会将延时归一化使最小延时为0 delays delays - min(delays); end在实际仿真中我们不会只发射一个角度。为了获得更好的合成效果我们会在一个对称的角度范围内均匀地发射多个不同角度的平面波。比如设定一个ang_range 20度即±10度然后在这个范围内均匀取ang_num 21个角度。角度越多合成效果理论上越好但计算量也越大需要权衡。%% 定义角度序列 userset.angrange 20; % 总角度范围度 userset.angnum 21; % 发射角度数量 Angles linspace(-deg2rad(userset.angrange/2), deg2rad(userset.angrange/2), userset.angnum); %% 循环进行多角度发射仿真 for i 1:userset.angnum txsteer Angles(i); % 当前发射角度 tx_d plane_wave_delayt(trans, txsteer); % 计算当前角度的发射延时 % 配置发射孔径 xdc_apodization(emit, 0, ones(1, trans.numele)); % 发射孔径变迹这里用均匀变迹即全开 xdc_center_focus(emit, [0 0 0]); % 设置孔径中心 xdc_focus_times(emit, 0, tx_d); % !!!关键步骤将计算好的延时施加到发射孔径上 % 配置接收孔径接收通常是全孔径、无延时聚焦等待回波 xdc_apodization(rcv, 0, ones(1, trans.numele)); xdc_focus_times(rcv, 0, zeros(1, trans.numele)); % ... 接下来调用 calc_scat 等函数计算散射回波 ... end踩坑提醒xdc_focus_times函数是施加延时的关键。它的第二个参数是焦距对于平面波我们设为0表示无限远焦点第三个参数就是我们的延时向量。这里一定要确保延时向量的长度和阵元数一致并且顺序对应。4. 核心拼图多角度回波数据的相干合成波束形成好了经过上一步的循环我们得到了21组回波数据每组数据对应一个特定发射角度。它们单独来看每一幅都是对比度不佳的“模糊”图像。现在我们要像一位高明的拼图师把这些碎片拼成一幅高清全景图。这个过程就是延时叠加波束形成也是复合平面波成像的精华所在。波束形成的目标是为图像中的每一个像素点(x, z)从所有阵元的所有采样时间中找到并累加与之对应的回波信号。这个“对应”关系就是声波从发射到散射点再回到接收阵元的总传播时间。对于复合平面波这个时间计算比单次发射复杂一点因为它包含了发射延时和接收延时两部分。总时间 (发射路径时间 接收路径时间) / 声速对于一次偏转角度为TXangle的发射声波从探头平面传播到像素点(x, z)的路径并不是简单的直线。由于是平面波发射我们可以将发射过程等效为一个从虚拟源点发出的、方向为TXangle的平面波到达该像素点所需的时间。这个虚拟源点的位置与探头孔径有关。一个常用的简化计算公式如下% 假设 xT 是阵元位置坐标数组xT(end)是最大正坐标值 halfaper sign(TXangle) * xT(end); % 孔径半长考虑偏转方向 dTX z * cos(TXangle) (x halfaper) * sin(TXangle); % 发射路径距离 dRX sqrt((xT - x).^2 z.^2); % 接收路径距离每个阵元到像素点的距离 tau (dTX dRX) / trans.c; % 总传播时间这里dTX的计算是理解关键。它利用了平面波的特性波前是直线。z*cos(TXangle)是轴向深度分量(xhalfaper)*sin(TXangle)是横向偏移分量。halfaper的引入是为了将计算参考点对齐确保公式在正负角度下都正确。但请注意Field II 函数calc_scat或calc_scat_multi在返回回波数据时还会返回一个tstart参数。这个参数代表了数据采集开始的绝对时间零点包含了硬件触发、电缆延时等。因此在实际用于索引射频数据的延时应该是delay_for_index round((tau - tstart) * sampling_frequency);现在我们可以编写波束形成的核心循环了。为了提高计算效率我们通常先对射频RF数据进行希尔伯特变换得到解析信号然后进行插值读取。% 假设 rf_data 是三维数据 [采样点数, 阵元数, 发射角度数] % 假设 x, z 是待成像像素点的坐标向量 rf_hilbert hilbert(rf_data); % 对射频数据做希尔伯特变换得到复信号 % 初始化一个三维矩阵存储每个角度单独形成的图像 dasdata zeros(length(z), length(x), userset.angnum); for ii 1:userset.angnum current_rf rf_hilbert(:,:,ii); % 取出第ii个角度的数据 TXangle Angles(ii); halfaper sign(TXangle) * xT(end); % 为所有像素点计算延时矩阵 tau [像素点数 阵元数] % 这里使用向量化计算提升速度避免嵌套循环 Z_mat z(:); % 列向量 X_mat x(:); % 列向量 dTX_all Z_mat * cos(TXangle) (X_mat halfaper) * sin(TXangle); % [像素数1] dRX_all sqrt((xT.^2) (Z_mat.^2) - 2 * X_mat * xT (X_mat.^2)); % 利用广播[像素数阵元数] tau_all (dTX_all dRX_all) ./ trans.c; % [像素数阵元数] % 将时间转换为采样点索引 sample_index round((tau_all - tstart(ii)) * userset.fs); sample_index max(1, min(size(current_rf,1), sample_index)); % 限制索引在有效范围内 % 延时叠加对每个像素点从每个阵元数据中取出对应时刻的信号并求和 for pix 1:length(X_mat) for ele 1:trans.numele idx sample_index(pix, ele); dasdata(pix, :, ii) dasdata(pix, :, ii) current_rf(idx, ele); end end end % 最终将所有角度合成的图像进行相干叠加复数求和 final_image sum(dasdata, 3); final_image_amplitude abs(final_image); % 取模值得幅度图像注意上面的双循环是为了清晰展示原理在实际编程中利用MATLAB的矩阵运算和interp1插值函数可以极大地优化这部分代码实现向量化操作速度能快几十倍甚至上百倍。核心思想是构建一个索引矩阵然后一次性从射频数据矩阵中“抓取”出所有需要的信号值。5. 效果对比与参数调优让你的图像更清晰经过前面一番“折腾”我们终于得到了合成后的图像final_image_amplitude。是骡子是马拉出来溜溜。最直观的方法就是把它和单次平面波比如0度发射的图像放在一起对比。你可以用imagesc函数显示图像并用dB刻度20*log10(图像/最大值)来更好地展示动态范围。通常你会观察到以下几个方面的显著改善分辨率提升点目标的尺寸点扩散函数会明显变小变得更“锐利”。旁瓣抑制目标周围的伪影旁瓣会减弱背景更干净。对比度提升弱散射体与背景的区分度更高。但是复合平面波成像的效果并不是简单地“角度越多越好”。这里面有几个关键参数需要你根据实际情况仔细调优这也是体现工程经验的地方角度范围 (ang_range)这个范围不能无限大。它受到探头孔径和阵元间距的物理限制。根据相控阵理论最大无栅瓣偏转角约为arcsin(λ / (2*pitch))。对于我们的P7-4探头λ≈0.256mm pitch0.171mm这个角度大概在±45度左右。但实际应用中通常使用±10度到±20度的范围因为角度太大时波束质量会下降边缘区域的信号强度太弱。角度数量 (ang_num)在确定的角度范围内发射多少个角度。数量越多合成增益越高图像信噪比改善越明显但仿真或数据采集时间也线性增加。我个人的经验是在±15度范围内取31或41个角度已经能获得非常好的效果再增加数量带来的边际效益就很小了。你可以做一个实验分别用5 11 21 41个角度合成图像观察图像质量的变化曲线找到那个“性价比”最高的点。合成策略我们上面做的是相干合成即对复数信号直接求和。这是最常用、效果也最好的方法。但还有一种叫非相干合成即先对每个角度形成的图像取幅度abs然后再求和。非相干合成对运动伪影的鲁棒性稍好但分辨率提升效果不如相干合成。在实际中如果目标如心脏运动很快可能需要考虑更复杂的运动补偿算法与合成策略的结合。为了更直观地展示参数影响我们可以设计一个小实验参数组合角度范围角度数量预期效果计算成本组合A±5°11分辨率略有提升旁瓣抑制一般低组合B±10°21分辨率、对比度显著提升效果均衡中组合C±15°41图像非常清晰背景干净接近理想聚焦高组合D±20°61边缘图像质量可能下降计算耗时很长极高调优建议我建议你从“组合B”开始。先跑通整个流程得到一幅比单次平面波好得多的图像建立信心。然后固定角度范围逐步增加角度数量观察图像质量的饱和点。再固定一个合适的角度数量去微调角度范围。这个过程可能需要反复几次但能让你深刻理解每个参数是如何影响最终成像结果的。仿真的一大优势就是可以低成本地做这些“参数扫描”实验找到最适合你当前仿真目标比如仿血管、仿囊肿的最佳参数集。最后别忘了处理图像显示。由于相控阵的成像区域是扇形的而我们的像素网格(x, z)可能是矩形的。在显示前需要将矩形网格中落在扇形区域外的像素点屏蔽掉通常设为NaN或0。可以使用inpolygon函数判断像素点是否在扇形内或者直接根据像素点的极坐标角度atan2(x, z)是否在探头最大偏转角度内来判断。处理好这个细节你的成像结果才会看起来专业又美观。