1. 项目概述聊聊我今天想说的这个项目——“基于MATLAB的延迟求和DAS波束形成算法实现”。很多人第一次看到“波束形成”四个字会觉得是不是军工或者雷达才用的东西其实它在医学超声、相控阵雷达、声呐、麦克风阵列语音增强里都有极其广泛的应用。本质上它就是让一组传感器阵元通过“对齐信号到达时间差”后叠加在一起从而在某个方向上获得更强的有效信号同时把杂散的干扰压下去。你手机上的麦克风降噪、医院里的超声探头发射接收背后都离不开这个逻辑。这个项目适合三类人一类是刚接触阵列信号处理的学生想搞懂原理但教材太拗口另一类是正在做超声成像仿真或者雷达信号处理仿真的工程师需要一个可靠的基线方案去做对比还有一类就是想在MATLAB里把一个信号处理管线完整跑一遍、拿到能出图的成果用来写报告或者做毕业设计的人。这个项目最大的价值是算法本身不复杂但它是整个成像系统的心脏把DAS写成代码、调出清晰图像的过程能让你把信号处理里“延时、插值、叠加、加窗”这四件事彻底吃透。2. DAS波束形成的核心设计与思路拆解2.1 为什么需要阵列单换能器不是也能成像吗我这么说吧单阵元超声探头在医学成像里确实是存在的但它只能沿一条线发射和接收扫出来的图像就像用一根针去探测世界既慢又看不清深部结构。阵列探头则像一排“眼睛”同时从不同位置观察目标区域再把这些数据用算法合成出一张图像。合成的过程就叫波束形成。对于超声成像来说核心需求是“聚焦”发射的时候我们希望声能在目标点附近汇聚得集中一点接收的时候我们希望回波信号里来自目标点的成分被强调而来自其他位置的杂波被抑制。DAS就是完成这两件事的最基础算法。有一个非常直观的类比。一群人在操场上喊一句话如果有人站在远处听他会听到混乱的声音叠在一起但如果每个人都事先安排好自己的“喊话时间”让启动信号到达听者耳朵的那一瞬间是同时的那听者会听到一个清晰响亮的声音。DAS干的就是这个“安排时间”的工作——它给每个阵元信号做不同的时间偏移让它们在某一个焦点上“时钟对齐”再把对齐后的信号求和。2.2 DAS的数学本质和工程实现路径我们用一维线性阵列来理解这个算法这也是最经典的情形。假设阵列有 (N) 个阵元间距是固定的 (d)阵元位置为 (x_n)目标聚焦点位于 (\vec{r_f})。那么第 (n) 个阵元到焦点位置的往返时间可以写成[ t_n(\vec{r_f}) \frac{|\vec{r_f} - \vec{r_n}|}{c} \times 2 ]这里 (c) 是声速超声软组织一般取1540 m/s乘2是因为声波走了一个来回发射到目标再反射回来。每个阵元记录下来的接收信号是 (s_n(t))我们要对第 (n) 个阵元的信号做一个时间移位让它们都在 (t_n) 这个时刻对齐到焦点上然后累加所有阵元的输出。[ S_{DAS}(\vec{r_f}) \sum_{n1}^{N} w_n \cdot s_n(t - t_n(\vec{r_f})) ](w_n) 是权重系数最常见的就是等权重矩形窗或者按幅度锥削Hanning、Hamming窗来降低旁瓣。MATLAB实现的时候最核心的步骤有三个延时计算、信号采样点对齐插值、加权叠加。在实际工程里我们会发现 (t_n) 乘以采样率 (f_s) 后并不一定是整数。比如某个阵元的延时是 (3.72) 个采样周期那信号这个点就不落在整数时间坐标上就需要插值。MATLAB里常用的做法是用interp1做线性插值或者用resample和delayseq来做更平滑的延迟。精度要求再高的话可以通过频域移相来实现分数延时。2.3 为什么选MATLAB而不是其他语言来做这个项目很多人在做算法仿真时都会纠结到底用C、Python还是MATLAB我的经验是如果目标是“理解并验证算法”那MATLAB在这类阵列信号处理项目里的优势是压倒性的——矩阵操作是语言内置的xcorr、interp1、fir1、phasedelay这类函数拿来就能用画图不需要额外引库二维图像直接imagesc出图调参可以即时可视化。Python也不是不行但你需要自己组装numpy、scipy、matplotlib并且很多时候信号处理的细节要自己撸。对于入门验证和快速迭代MATLAB的交互式体验确实更好。当然如果你将来要把代码跑到实时系统上MATLAB的代码生成Coder也能把DAS算法转到C/C这个项目写的核心循环逻辑对代码生成还是比较友好的。3. 核心算法实现与关键参数精讲3.1 从数据流的角度理解整个实现我建议把整个DAS流程看作是“接收信号 → 多通道同步 → 延时补偿 → 幅度加权 → 相干累加 → 图像重建”的数据流。早期做实验中我犯过一个错误上来就写循环去一帧帧做延时结果代码慢到怀疑人生。后来把数据组织成一个大矩阵用向量化一次处理所有通道速度提升是数量级的。具体的矩阵化思路是这样把每个阵元的回波信号存成矩阵的行行数是阵元数 (N)列数是采样点数 (M)。对某个聚焦点我们计算所有通道的整型延时下标向量再用“查找索引线性插值”的方式一次性取出对应采样点。我通常用interp1(t_received, signal_matrix, t_aligned)来操作其中t_aligned是每个通道经过延时换算后的采样时间点向量。这里的核心角度要分清楚我们不是“在时间轴上截取一段”而是将焦点处的期望回波时间点映射到各个通道的接收信号中去“读出一个值”。一旦想通了这一点后面的实现就水到渠成了。3.2 关键的延时计算细节延时计算是整个DAS的灵魂。第一件事是确定焦点坐标和阵元坐标的参考系。我习惯把线阵阵元中心放在原点阵元沿 (x) 轴等距排开焦点的极坐标用 ((r, \theta)) 表示。这样一来单程传播距离就是[ d_n(r, \theta) \sqrt{(x_n - r\sin\theta)^2 (r\cos\theta)^2} ]我在仿真平面上把成像区域划分成横向若干个像素点纵向是深度。对每一个成像点都计算一次所有通道的延时。对了这里有个容易忽略的点如果目标点是焦点那针对每个像素点都计算一套延时这个过程的计算量是 (N_x \times N_z \times N)也就是图像宽 × 图像深 × 阵元数。所以写代码时一定要避免在像素循环内部重复计算那些跟通道位置无关的常量。我项目中常用的参数如下参数数值说明阵元数 (N)128模拟医学超声线阵的常见尺寸阵元间距 (d)0.3 mm换能器中心频率决定半波长左右中心频率 (f_0)5 MHz常规腹部超声频率声速 (c)1540 m/s人体软组织近似声速采样率 (f_s)50 MHz奈奎斯特条件至少要两倍于中心频率这里留足余量成像深度范围0~50 mm根据目标组织特性设定横向像素数200成像宽度的离散网格点数纵向像素数300成像深度的离散层数从 (t_n) 换算到采样下标时一定要记得是“往返时间”。我见过很多初学者直接用单程距离除以声速去做延时出来的图焦点全错位了。具体做法dt 1 / fs; t_n 2 * sqrt((x_n - x_focus)^2 z_focus^2) / c; sample_index round(t_n / dt) 1; % 1是因为MATLAB索引从1开始如果后续用interp1就不需要这个取整操作直接把t_aligned t - t_n作为查询点传入即可。3.3 插值方法的选择线性插值够不够关于插值我直接用结论绝大数情况下线性插值够用但要分场景。使用线性插值interp1(..., linear)需要保证采样率远高于奈奎斯特率否则会产生插值不平滑的问题。以 (f_s50) MHz、信号带宽 (B2) MHz 为例过采样倍率达到25倍线性插值的误差已经足够小实践中完全可用。如果信号带宽高、采样率紧贴奈奎斯特频率那就要考虑更高阶的方法比如立方样条插值或频域零填充实现分数延时。频域方法在MATLAB里可以用ifft(fft(signal) .* exp(-1j*2*pi*f*tau))实现任意整数、分数延时而且相位精度最好。代价是计算量更大。这个项目里我先用线性插值后续对比尺寸再换成频域方法二者做对比图放在报告里效果极佳。3.4 加权窗与图像质量的权衡我之前做过一组对比实验同一个仿真数据不加窗、加Hamming窗、加Chebyshev窗结果差别非常直观。不加窗时主瓣窄、成像分辨率最高但旁瓣很高——在医学超声里高旁瓣意味着浅层组织的能量会泄漏到深层去形成伪像。加窗本质上是牺牲一定的主瓣宽度换取旁瓣的显著压低。权重 (w_n) 用的是阵元位置的函数。矩形窗权重全为1Hamming窗权重边缘小、中心大。它的作用是对空间采样做幅度锥削。在实际成像中强烈建议用Hamming或Hann窗除非你要专门研究旁瓣指标。我代码里把窗函数做成一个参数可以随时切换来对比效果。switch window_type case rect w ones(1, N); case hamming w hamming(N); case hann w hann(N, periodic); case cheb w chebwin(N, 40); % 40dB旁瓣抑制 end刚才去对比不同窗函数的“分辨率-旁瓣”折中时可以参考一个量化指标主瓣宽度-3dB宽度和最强旁瓣电平。以128阵元、半波长间距的线阵为例矩形窗主瓣约为 (0.89 \cdot \lambda)最强旁瓣约为-13.3dBHamming窗主瓣大概 (1.33 \cdot \lambda)最强旁瓣能压到-42dB左右。这就是分辨率换低旁瓣的经典代价做实验报告时把这两个数字写进去会很有说服力。4. 实操过程从仿真信号到最终图像4.1 生成模拟超声回波信号在没有实测数据的情况下先用仿真数据验证算法是正规路子。我采用“点散射子模型”——强反射小目标模拟血管壁、组织边界等按已知位置放好每个散射点看作一个回波源接收信号就是这些回波经过延时、衰减和噪声叠加后的总和。核心思路是对每个散射点计算它到每个阵元的往返时间再把发射信号比如一个带通调制的脉冲延时到对应位置并累加。这里的发射脉冲我常设置为中心频率5 MHz、带宽2 MHz的高斯调制正弦波。% 生成发射脉冲 pulse_freq 5e6; pulse_bandwidth 2e6; t_pulse -1e-6:dt:1e-6; pulse cos(2 * pi * pulse_freq * t_pulse) .* exp(-t_pulse.^2 * pulse_bandwidth^2); % 放置目标散射点 targets [0, 10e-3, 1.0, 0 5e-3, 20e-3, 0.8, 1 -5e-3, 20e-3, 0.7, 1 2e-3, 35e-3, 0.5, 2]; % 每一行x位置, z位置, 幅度, 散射强弱循环里就是计算每个阵列通道的往返时间然后把散射幅度加到对应时刻的样点上。不要忘了在最后的接收信号上叠加高斯白噪声来模拟电子噪声。散射点的分布设计是有讲究的我故意放了三个离得较近的散射点相距1~2 mm用来测试系统的分辨率能力又放了不同深度的点用来观察聚焦效果随深度的变化。这种仿真布置使我能一眼看出图像的轴向差异和横向分辨情况比只放一个点全面得多。4.2 核心DAS函数编写把DAS封装成一个函数是让整个项目可复用、可调参的正确做法。函数输入接收矩阵、阵列几何参数、成像范围输出重建图像。代码核心部分可以直接这样写function image das_bf(received_signal, x_elem, z_start, z_end, N_x, N_z, c, fs) % received_signal: N_elem x N_samples 矩阵 % x_elem: N_elem x 1阵元x坐标 % z_start/z_end: 成像深度范围 N_elem length(x_elem); z_axis linspace(z_start, z_end, N_z); x_axis linspace(min(x_elem)*0.8, max(x_elem)*0.8, N_x); image zeros(N_x, N_z); t_axis (0:size(received_signal,2)-1) / fs; for ix 1:N_x for iz 1:N_z x_focus x_axis(ix); z_focus z_axis(iz); integrated 0; for ne 1:N_elem dist sqrt((x_elem(ne) - x_focus)^2 z_focus^2); t_delay 2 * dist / c; t_aligned t_axis t_delay; % 若采样点超出接收时间范围则跳过 if t_aligned(end) t_axis(end) t_aligned(1) t_axis(1) sample interp1(t_axis, received_signal(ne,:), t_aligned, linear, 0); integrated integrated w(ne) * sample; end end image(ix, iz) integrated; end end end注意几个细节interp1的最后一个参数填0作用是防止查询点超出范围报NaNt_aligned(end) t_axis(end)这个判断是防止焦点太深导致延时超出采样时间延时的查询方向是t_axis t_delay而不是t_axis - t_delay因为接收信号里的目标回波本来就比参考时间起点晚出现我们要“往后找”。这里我一开始老搞错方向总感觉对信号做延时应该减掉某个值后来从物理上推了一遍才彻底搞清楚阵列中某个阵元接收到的回波和通道0的接收回波相比如果距离焦点更远那它的回波到达时间就更靠后。要在叠加时把这个靠后的时间“补偿回去”让各通道在某时刻对齐到同一个焦点就要把信号往前面拉也就是在时间轴上去查询更晚时间点的值。所以查询点是t_axis t_delay。反过来做就全错了。4.3 图像重建与后处理DAS输出的是一个二维矩阵但直接用原始矩阵画图的话动态范围太大近场强反射会把深层细节全压成黑色。因此需要做动态范围压缩常用的做法是取对数包络并设置显示窗口。流程是这样先对重建矩阵取绝对值实信号可以用Hilbert变换得到包络再转成分贝最后减去一个门槛来设置显示上下限。env abs(hilbert(image)); % 包络提取 log_image 20 * log10(env / max(env(:))); % 归一化分贝图 disp_range 40; % 显示动态范围40dB log_image(log_image -disp_range) -disp_range; figure; imagesc(x_axis*1e3, z_axis*1e3, log_image); colormap(gray); axis tight; xlabel(横向位置 (mm)); ylabel(深度 (mm));这里用的Hilbert包络在实际超声图像形成中非常重要因为射频信号是振荡的直接取绝对值会产生振铃状纹理取包络后才能显示干净的组织结构。要是忘了做包络图像会显得特别“噪”这也是很多新手代码出图后感觉不对劲的主要原因。4.4 针对不同参数组合做实验的注意事项写DAS代码时每换一组参数就需要重新计算整套延时。我建议把延时表针对每个聚焦点、每个阵元的延时下标先算好并缓存起来不要在每次重建图像时反复算。由于聚焦点数量 (N_x \times N_z) 乘以阵元数很容易到百万级循环计算虽然耗时但MATLAB矩阵化能一次处理所有通道。一次优化创造的性能差距很大。以200×300像素 × 128阵元为例单纯三层循环计算要跑几秒钟甚至更久而矩阵方式预先用meshgrid生成坐标矩阵再一次性计算距离和延时通常能把单帧降到零点几秒。这在调试时感受特别明显矩阵方式能做到改一个参数立即看到效果然后我再对图像改动做微调整个过程顺滑得多。5. 常见问题与排查技巧实录5.1 图像出现横向错位或目标位置偏移这个问题的根源九成是延时方向搞反了。你要时刻记住焦点的距离越远回波到各阵元的延迟差越大补偿是要“后对齐”不是“前对齐”。我用一个简单测试来区分在正中央放置一个散射点重建后图像点如果落在了中央偏左或偏右的位置那多半就是延时符号反了。另一种可能是阵元坐标和焦点坐标用的参考系不一致比如一个用毫米一个用米导致距离计算完全错误。5.2 图像整体模糊不清、分辨率很差排除参数问题后最可能的原因是阵元间距太大。阵元间距如果超过半波长阵列方向图会出现栅瓣像图像里出现周期性的伪目标。超声探头的阵元间距一般设计为中心频率的半波长也就是 (d \lambda/2 c/(2f_0))。我们的项目里 (f_05) MHz半波长约0.154 mm。如果仿真里为了省事设成1 mm那就直接遇到栅瓣效应了。同时还有个细节如果阵元数太少比如只有8阵元DAS的主瓣会非常宽两个点很难分辨开。经验上要分辨某个横向距离 (\Delta x)阵列的孔径长度 (L Nd) 至少要接近目标距离处可分辨的宽度要求。比如孔径10 mm在30 mm深处横向分辨率大概在1 mm量级这是由衍射极限决定的。做定量分析时可以用瑞利判据估算一下别指望有限孔径能突破物理极限。5.3 图像动态范围拉不开背景噪声没压住这种情况首先要确认是不是噪声加得太猛。信号噪声比SNR低于10 dB的时候即使DAS本身有阵列增益重建图像也可能被淹没。阵元数128时DAS对不相关噪声的理论增益是 (10\log_{10}(128) \approx 21) dB这个数字可以拿来预估需要多高的输入SNR才能稳住图像。其次是没有做包络提取就直接imagesc导致图像显示为振荡条纹看起来就像“噪声重”。还有就是动态范围设置得太小或者太大设太小时对比度不足设太大时弱散射点细节会被近场强反射遮盖。一般静态仿真取40~60 dB比较好。5.4 运算速度太慢循环代码半天跑不完如果代码确实是按最直观的方式写的——三重循环每个像素每个阵元都独立计算那性能瓶颈就是必然的。解决三步走第一把散射点的累加放到矩阵运算里用sum(..., 2)第二对每个焦点点的延时做向量化第三对于工程要求较高的场景可以用parfor替代for来并行处理成像点循环。MATLAB的并行计算工具箱在“成像点之间互不依赖”的任务里加速比很可观。5.5 焦点离阵列太近时出现伪影近场区的波形其实不再是简单的球面波而更像柱面波或更复杂的波场。标准DAS假设的是理想的球面波传播路径当焦点离阵列几个波长以内时这种近似就失效了。这也是为什么超声探头前面要做匹配层、声透镜的原因之一。我在仿真里会做一个判断焦点距阵列表面至少大于10倍波长否则对结果要去打一个预期较低的标记具体项目中这一条不太会引起混淆但如果要做近场成像比如光声显微镜就要换分层介质模型来处理了。6. 性能对比验证与进阶扩展方向6.1 用点散射目标检验DAS重建质量回到我那个模拟数据我在中央深度20 mm处放了两个横向距离1.5 mm的点目标。用50 MHz采样率、128阵元、Hamming窗重建之后可以从图像上直观看清两个点。我建议在这个测试里记录两项指标峰值旁瓣电平PSLR和横向分辨率-6dB宽度。用MATLAB在图像峰值附近取横向切片数据找到主瓣下降6dB和旁瓣峰值的幅值记录下来这个就是最能拿得出手的定量结果。对比矩形窗和Hamming窗的结果矩形窗旁瓣在-13.3 dB左右Hamming窗能到-40多dB。做报告时下这种定量结论说服力比“图像看起来更干净了”强得多。6.2 自适应波束形成的下一步扩展DAS做到能出图以后自然的进阶方向是自适应波束形成。最常见的是MVDR最小方差无失真响应它根据采样协方差矩阵动态计算权重来抑制干扰。MVDR能把分辨率再往上推但对模型误差也非常敏感协方差矩阵估计不准时性能反而会崩。所以实际应用里常做对角加载来增强稳健性。如果你把这个项目当课程设计或者毕设我建议的路线是先在仿真里实现DAS并把图像定量指标测好再实现MVDR并对比两者的分辨率、旁瓣、抗噪性。对比表一放“为什么现在医用超声大多还在用DAS类的方法而不是MVDR”这个问题就能引出一个相当有深度的讨论——因为实时系统对计算量和稳健性的要求极高自适应方法虽然理论好但落地难。这个论点特别能体现你对工程问题的理解。6.3 光声/超声联合成像延伸另一个有意思的扩展方向是把DAS移植到光声成像中。光声成像的物理过程是短脉冲激光打到组织上组织吸收光能产生超声信号再用超声阵列接收重建图像。它的回波模型和超声不太一样光声信号是单程传播从目标到阵元所以延时计算去掉2倍距离其余逻辑基本不变。改一行核心代码就能跑光声的重建管线。做过超声DAS之后再去做光声DAS几乎是无缝切换这种能力会很受实验室导师喜欢。延时计算从2*d/c改成d/c声速换成对应介质速度再把光源的分布考虑进去一个基础的光声重建仿真就出来了。基于这个框架还能调整光通量分布和超声换能器的频率响应进一步与真实系统比对。6.4 关于声速不均匀问题的坦白最后说一个所有超声成像越不过去的话题人体组织声速并不是均匀的1540 m/s脂肪组织约1450 m/s骨骼约3000 m/s以上DAS建立在均匀声速假设上遇到组织不均匀图像就会变形、模糊。现在临床超声里好多伪像都源于此。你在项目总结里若把“DAS假设均匀声速”这个局限单独提出来并说清楚下一步可能的补偿方法是“逐层声速校正”那文章的专业层次立刻就不一样了。我个人在实际操作中的体会是DAS这个算法看起来简单但这套“先对齐后叠加”的逻辑是整个阵列信号处理的主心骨。把这个项目完整做下来你对延迟、加窗、插值、阵列几何空间关系的理解会有一个质的飞跃。最后再分享一个小技巧——连续调参数对比图像时一定要把窗口函数、动态范围、声速这些参数固化下来只改要研究的那一项否则很容易被视觉上的差异误导。这只是一个入门级波束形成器真正吃透它的人再去看自适应波束形成、平面波成像、合成孔径成像都会觉得一路畅通。 
长河君 · 2026/10/16 13:27:07 · 阅读 1.2k 基于MATLAB的DAS波束形成算法实现与超声成像仿真
1. 项目概述聊聊我今天想说的这个项目——“基于MATLAB的延迟求和DAS波束形成算法实现”。很多人第一次看到“波束形成”四个字会觉得是不是军工或者雷达才用的东西其实它在医学超声、相控阵雷达、声呐、麦克风阵列语音增强里都有极其广泛的应用。本质上它就是让一组传感器阵元通过“对齐信号到达时间差”后叠加在一起从而在某个方向上获得更强的有效信号同时把杂散的干扰压下去。你手机上的麦克风降噪、医院里的超声探头发射接收背后都离不开这个逻辑。这个项目适合三类人一类是刚接触阵列信号处理的学生想搞懂原理但教材太拗口另一类是正在做超声成像仿真或者雷达信号处理仿真的工程师需要一个可靠的基线方案去做对比还有一类就是想在MATLAB里把一个信号处理管线完整跑一遍、拿到能出图的成果用来写报告或者做毕业设计的人。这个项目最大的价值是算法本身不复杂但它是整个成像系统的心脏把DAS写成代码、调出清晰图像的过程能让你把信号处理里“延时、插值、叠加、加窗”这四件事彻底吃透。2. DAS波束形成的核心设计与思路拆解2.1 为什么需要阵列单换能器不是也能成像吗我这么说吧单阵元超声探头在医学成像里确实是存在的但它只能沿一条线发射和接收扫出来的图像就像用一根针去探测世界既慢又看不清深部结构。阵列探头则像一排“眼睛”同时从不同位置观察目标区域再把这些数据用算法合成出一张图像。合成的过程就叫波束形成。对于超声成像来说核心需求是“聚焦”发射的时候我们希望声能在目标点附近汇聚得集中一点接收的时候我们希望回波信号里来自目标点的成分被强调而来自其他位置的杂波被抑制。DAS就是完成这两件事的最基础算法。有一个非常直观的类比。一群人在操场上喊一句话如果有人站在远处听他会听到混乱的声音叠在一起但如果每个人都事先安排好自己的“喊话时间”让启动信号到达听者耳朵的那一瞬间是同时的那听者会听到一个清晰响亮的声音。DAS干的就是这个“安排时间”的工作——它给每个阵元信号做不同的时间偏移让它们在某一个焦点上“时钟对齐”再把对齐后的信号求和。2.2 DAS的数学本质和工程实现路径我们用一维线性阵列来理解这个算法这也是最经典的情形。假设阵列有 (N) 个阵元间距是固定的 (d)阵元位置为 (x_n)目标聚焦点位于 (\vec{r_f})。那么第 (n) 个阵元到焦点位置的往返时间可以写成[ t_n(\vec{r_f}) \frac{|\vec{r_f} - \vec{r_n}|}{c} \times 2 ]这里 (c) 是声速超声软组织一般取1540 m/s乘2是因为声波走了一个来回发射到目标再反射回来。每个阵元记录下来的接收信号是 (s_n(t))我们要对第 (n) 个阵元的信号做一个时间移位让它们都在 (t_n) 这个时刻对齐到焦点上然后累加所有阵元的输出。[ S_{DAS}(\vec{r_f}) \sum_{n1}^{N} w_n \cdot s_n(t - t_n(\vec{r_f})) ](w_n) 是权重系数最常见的就是等权重矩形窗或者按幅度锥削Hanning、Hamming窗来降低旁瓣。MATLAB实现的时候最核心的步骤有三个延时计算、信号采样点对齐插值、加权叠加。在实际工程里我们会发现 (t_n) 乘以采样率 (f_s) 后并不一定是整数。比如某个阵元的延时是 (3.72) 个采样周期那信号这个点就不落在整数时间坐标上就需要插值。MATLAB里常用的做法是用interp1做线性插值或者用resample和delayseq来做更平滑的延迟。精度要求再高的话可以通过频域移相来实现分数延时。2.3 为什么选MATLAB而不是其他语言来做这个项目很多人在做算法仿真时都会纠结到底用C、Python还是MATLAB我的经验是如果目标是“理解并验证算法”那MATLAB在这类阵列信号处理项目里的优势是压倒性的——矩阵操作是语言内置的xcorr、interp1、fir1、phasedelay这类函数拿来就能用画图不需要额外引库二维图像直接imagesc出图调参可以即时可视化。Python也不是不行但你需要自己组装numpy、scipy、matplotlib并且很多时候信号处理的细节要自己撸。对于入门验证和快速迭代MATLAB的交互式体验确实更好。当然如果你将来要把代码跑到实时系统上MATLAB的代码生成Coder也能把DAS算法转到C/C这个项目写的核心循环逻辑对代码生成还是比较友好的。3. 核心算法实现与关键参数精讲3.1 从数据流的角度理解整个实现我建议把整个DAS流程看作是“接收信号 → 多通道同步 → 延时补偿 → 幅度加权 → 相干累加 → 图像重建”的数据流。早期做实验中我犯过一个错误上来就写循环去一帧帧做延时结果代码慢到怀疑人生。后来把数据组织成一个大矩阵用向量化一次处理所有通道速度提升是数量级的。具体的矩阵化思路是这样把每个阵元的回波信号存成矩阵的行行数是阵元数 (N)列数是采样点数 (M)。对某个聚焦点我们计算所有通道的整型延时下标向量再用“查找索引线性插值”的方式一次性取出对应采样点。我通常用interp1(t_received, signal_matrix, t_aligned)来操作其中t_aligned是每个通道经过延时换算后的采样时间点向量。这里的核心角度要分清楚我们不是“在时间轴上截取一段”而是将焦点处的期望回波时间点映射到各个通道的接收信号中去“读出一个值”。一旦想通了这一点后面的实现就水到渠成了。3.2 关键的延时计算细节延时计算是整个DAS的灵魂。第一件事是确定焦点坐标和阵元坐标的参考系。我习惯把线阵阵元中心放在原点阵元沿 (x) 轴等距排开焦点的极坐标用 ((r, \theta)) 表示。这样一来单程传播距离就是[ d_n(r, \theta) \sqrt{(x_n - r\sin\theta)^2 (r\cos\theta)^2} ]我在仿真平面上把成像区域划分成横向若干个像素点纵向是深度。对每一个成像点都计算一次所有通道的延时。对了这里有个容易忽略的点如果目标点是焦点那针对每个像素点都计算一套延时这个过程的计算量是 (N_x \times N_z \times N)也就是图像宽 × 图像深 × 阵元数。所以写代码时一定要避免在像素循环内部重复计算那些跟通道位置无关的常量。我项目中常用的参数如下参数数值说明阵元数 (N)128模拟医学超声线阵的常见尺寸阵元间距 (d)0.3 mm换能器中心频率决定半波长左右中心频率 (f_0)5 MHz常规腹部超声频率声速 (c)1540 m/s人体软组织近似声速采样率 (f_s)50 MHz奈奎斯特条件至少要两倍于中心频率这里留足余量成像深度范围0~50 mm根据目标组织特性设定横向像素数200成像宽度的离散网格点数纵向像素数300成像深度的离散层数从 (t_n) 换算到采样下标时一定要记得是“往返时间”。我见过很多初学者直接用单程距离除以声速去做延时出来的图焦点全错位了。具体做法dt 1 / fs; t_n 2 * sqrt((x_n - x_focus)^2 z_focus^2) / c; sample_index round(t_n / dt) 1; % 1是因为MATLAB索引从1开始如果后续用interp1就不需要这个取整操作直接把t_aligned t - t_n作为查询点传入即可。3.3 插值方法的选择线性插值够不够关于插值我直接用结论绝大数情况下线性插值够用但要分场景。使用线性插值interp1(..., linear)需要保证采样率远高于奈奎斯特率否则会产生插值不平滑的问题。以 (f_s50) MHz、信号带宽 (B2) MHz 为例过采样倍率达到25倍线性插值的误差已经足够小实践中完全可用。如果信号带宽高、采样率紧贴奈奎斯特频率那就要考虑更高阶的方法比如立方样条插值或频域零填充实现分数延时。频域方法在MATLAB里可以用ifft(fft(signal) .* exp(-1j*2*pi*f*tau))实现任意整数、分数延时而且相位精度最好。代价是计算量更大。这个项目里我先用线性插值后续对比尺寸再换成频域方法二者做对比图放在报告里效果极佳。3.4 加权窗与图像质量的权衡我之前做过一组对比实验同一个仿真数据不加窗、加Hamming窗、加Chebyshev窗结果差别非常直观。不加窗时主瓣窄、成像分辨率最高但旁瓣很高——在医学超声里高旁瓣意味着浅层组织的能量会泄漏到深层去形成伪像。加窗本质上是牺牲一定的主瓣宽度换取旁瓣的显著压低。权重 (w_n) 用的是阵元位置的函数。矩形窗权重全为1Hamming窗权重边缘小、中心大。它的作用是对空间采样做幅度锥削。在实际成像中强烈建议用Hamming或Hann窗除非你要专门研究旁瓣指标。我代码里把窗函数做成一个参数可以随时切换来对比效果。switch window_type case rect w ones(1, N); case hamming w hamming(N); case hann w hann(N, periodic); case cheb w chebwin(N, 40); % 40dB旁瓣抑制 end刚才去对比不同窗函数的“分辨率-旁瓣”折中时可以参考一个量化指标主瓣宽度-3dB宽度和最强旁瓣电平。以128阵元、半波长间距的线阵为例矩形窗主瓣约为 (0.89 \cdot \lambda)最强旁瓣约为-13.3dBHamming窗主瓣大概 (1.33 \cdot \lambda)最强旁瓣能压到-42dB左右。这就是分辨率换低旁瓣的经典代价做实验报告时把这两个数字写进去会很有说服力。4. 实操过程从仿真信号到最终图像4.1 生成模拟超声回波信号在没有实测数据的情况下先用仿真数据验证算法是正规路子。我采用“点散射子模型”——强反射小目标模拟血管壁、组织边界等按已知位置放好每个散射点看作一个回波源接收信号就是这些回波经过延时、衰减和噪声叠加后的总和。核心思路是对每个散射点计算它到每个阵元的往返时间再把发射信号比如一个带通调制的脉冲延时到对应位置并累加。这里的发射脉冲我常设置为中心频率5 MHz、带宽2 MHz的高斯调制正弦波。% 生成发射脉冲 pulse_freq 5e6; pulse_bandwidth 2e6; t_pulse -1e-6:dt:1e-6; pulse cos(2 * pi * pulse_freq * t_pulse) .* exp(-t_pulse.^2 * pulse_bandwidth^2); % 放置目标散射点 targets [0, 10e-3, 1.0, 0 5e-3, 20e-3, 0.8, 1 -5e-3, 20e-3, 0.7, 1 2e-3, 35e-3, 0.5, 2]; % 每一行x位置, z位置, 幅度, 散射强弱循环里就是计算每个阵列通道的往返时间然后把散射幅度加到对应时刻的样点上。不要忘了在最后的接收信号上叠加高斯白噪声来模拟电子噪声。散射点的分布设计是有讲究的我故意放了三个离得较近的散射点相距1~2 mm用来测试系统的分辨率能力又放了不同深度的点用来观察聚焦效果随深度的变化。这种仿真布置使我能一眼看出图像的轴向差异和横向分辨情况比只放一个点全面得多。4.2 核心DAS函数编写把DAS封装成一个函数是让整个项目可复用、可调参的正确做法。函数输入接收矩阵、阵列几何参数、成像范围输出重建图像。代码核心部分可以直接这样写function image das_bf(received_signal, x_elem, z_start, z_end, N_x, N_z, c, fs) % received_signal: N_elem x N_samples 矩阵 % x_elem: N_elem x 1阵元x坐标 % z_start/z_end: 成像深度范围 N_elem length(x_elem); z_axis linspace(z_start, z_end, N_z); x_axis linspace(min(x_elem)*0.8, max(x_elem)*0.8, N_x); image zeros(N_x, N_z); t_axis (0:size(received_signal,2)-1) / fs; for ix 1:N_x for iz 1:N_z x_focus x_axis(ix); z_focus z_axis(iz); integrated 0; for ne 1:N_elem dist sqrt((x_elem(ne) - x_focus)^2 z_focus^2); t_delay 2 * dist / c; t_aligned t_axis t_delay; % 若采样点超出接收时间范围则跳过 if t_aligned(end) t_axis(end) t_aligned(1) t_axis(1) sample interp1(t_axis, received_signal(ne,:), t_aligned, linear, 0); integrated integrated w(ne) * sample; end end image(ix, iz) integrated; end end end注意几个细节interp1的最后一个参数填0作用是防止查询点超出范围报NaNt_aligned(end) t_axis(end)这个判断是防止焦点太深导致延时超出采样时间延时的查询方向是t_axis t_delay而不是t_axis - t_delay因为接收信号里的目标回波本来就比参考时间起点晚出现我们要“往后找”。这里我一开始老搞错方向总感觉对信号做延时应该减掉某个值后来从物理上推了一遍才彻底搞清楚阵列中某个阵元接收到的回波和通道0的接收回波相比如果距离焦点更远那它的回波到达时间就更靠后。要在叠加时把这个靠后的时间“补偿回去”让各通道在某时刻对齐到同一个焦点就要把信号往前面拉也就是在时间轴上去查询更晚时间点的值。所以查询点是t_axis t_delay。反过来做就全错了。4.3 图像重建与后处理DAS输出的是一个二维矩阵但直接用原始矩阵画图的话动态范围太大近场强反射会把深层细节全压成黑色。因此需要做动态范围压缩常用的做法是取对数包络并设置显示窗口。流程是这样先对重建矩阵取绝对值实信号可以用Hilbert变换得到包络再转成分贝最后减去一个门槛来设置显示上下限。env abs(hilbert(image)); % 包络提取 log_image 20 * log10(env / max(env(:))); % 归一化分贝图 disp_range 40; % 显示动态范围40dB log_image(log_image -disp_range) -disp_range; figure; imagesc(x_axis*1e3, z_axis*1e3, log_image); colormap(gray); axis tight; xlabel(横向位置 (mm)); ylabel(深度 (mm));这里用的Hilbert包络在实际超声图像形成中非常重要因为射频信号是振荡的直接取绝对值会产生振铃状纹理取包络后才能显示干净的组织结构。要是忘了做包络图像会显得特别“噪”这也是很多新手代码出图后感觉不对劲的主要原因。4.4 针对不同参数组合做实验的注意事项写DAS代码时每换一组参数就需要重新计算整套延时。我建议把延时表针对每个聚焦点、每个阵元的延时下标先算好并缓存起来不要在每次重建图像时反复算。由于聚焦点数量 (N_x \times N_z) 乘以阵元数很容易到百万级循环计算虽然耗时但MATLAB矩阵化能一次处理所有通道。一次优化创造的性能差距很大。以200×300像素 × 128阵元为例单纯三层循环计算要跑几秒钟甚至更久而矩阵方式预先用meshgrid生成坐标矩阵再一次性计算距离和延时通常能把单帧降到零点几秒。这在调试时感受特别明显矩阵方式能做到改一个参数立即看到效果然后我再对图像改动做微调整个过程顺滑得多。5. 常见问题与排查技巧实录5.1 图像出现横向错位或目标位置偏移这个问题的根源九成是延时方向搞反了。你要时刻记住焦点的距离越远回波到各阵元的延迟差越大补偿是要“后对齐”不是“前对齐”。我用一个简单测试来区分在正中央放置一个散射点重建后图像点如果落在了中央偏左或偏右的位置那多半就是延时符号反了。另一种可能是阵元坐标和焦点坐标用的参考系不一致比如一个用毫米一个用米导致距离计算完全错误。5.2 图像整体模糊不清、分辨率很差排除参数问题后最可能的原因是阵元间距太大。阵元间距如果超过半波长阵列方向图会出现栅瓣像图像里出现周期性的伪目标。超声探头的阵元间距一般设计为中心频率的半波长也就是 (d \lambda/2 c/(2f_0))。我们的项目里 (f_05) MHz半波长约0.154 mm。如果仿真里为了省事设成1 mm那就直接遇到栅瓣效应了。同时还有个细节如果阵元数太少比如只有8阵元DAS的主瓣会非常宽两个点很难分辨开。经验上要分辨某个横向距离 (\Delta x)阵列的孔径长度 (L Nd) 至少要接近目标距离处可分辨的宽度要求。比如孔径10 mm在30 mm深处横向分辨率大概在1 mm量级这是由衍射极限决定的。做定量分析时可以用瑞利判据估算一下别指望有限孔径能突破物理极限。5.3 图像动态范围拉不开背景噪声没压住这种情况首先要确认是不是噪声加得太猛。信号噪声比SNR低于10 dB的时候即使DAS本身有阵列增益重建图像也可能被淹没。阵元数128时DAS对不相关噪声的理论增益是 (10\log_{10}(128) \approx 21) dB这个数字可以拿来预估需要多高的输入SNR才能稳住图像。其次是没有做包络提取就直接imagesc导致图像显示为振荡条纹看起来就像“噪声重”。还有就是动态范围设置得太小或者太大设太小时对比度不足设太大时弱散射点细节会被近场强反射遮盖。一般静态仿真取40~60 dB比较好。5.4 运算速度太慢循环代码半天跑不完如果代码确实是按最直观的方式写的——三重循环每个像素每个阵元都独立计算那性能瓶颈就是必然的。解决三步走第一把散射点的累加放到矩阵运算里用sum(..., 2)第二对每个焦点点的延时做向量化第三对于工程要求较高的场景可以用parfor替代for来并行处理成像点循环。MATLAB的并行计算工具箱在“成像点之间互不依赖”的任务里加速比很可观。5.5 焦点离阵列太近时出现伪影近场区的波形其实不再是简单的球面波而更像柱面波或更复杂的波场。标准DAS假设的是理想的球面波传播路径当焦点离阵列几个波长以内时这种近似就失效了。这也是为什么超声探头前面要做匹配层、声透镜的原因之一。我在仿真里会做一个判断焦点距阵列表面至少大于10倍波长否则对结果要去打一个预期较低的标记具体项目中这一条不太会引起混淆但如果要做近场成像比如光声显微镜就要换分层介质模型来处理了。6. 性能对比验证与进阶扩展方向6.1 用点散射目标检验DAS重建质量回到我那个模拟数据我在中央深度20 mm处放了两个横向距离1.5 mm的点目标。用50 MHz采样率、128阵元、Hamming窗重建之后可以从图像上直观看清两个点。我建议在这个测试里记录两项指标峰值旁瓣电平PSLR和横向分辨率-6dB宽度。用MATLAB在图像峰值附近取横向切片数据找到主瓣下降6dB和旁瓣峰值的幅值记录下来这个就是最能拿得出手的定量结果。对比矩形窗和Hamming窗的结果矩形窗旁瓣在-13.3 dB左右Hamming窗能到-40多dB。做报告时下这种定量结论说服力比“图像看起来更干净了”强得多。6.2 自适应波束形成的下一步扩展DAS做到能出图以后自然的进阶方向是自适应波束形成。最常见的是MVDR最小方差无失真响应它根据采样协方差矩阵动态计算权重来抑制干扰。MVDR能把分辨率再往上推但对模型误差也非常敏感协方差矩阵估计不准时性能反而会崩。所以实际应用里常做对角加载来增强稳健性。如果你把这个项目当课程设计或者毕设我建议的路线是先在仿真里实现DAS并把图像定量指标测好再实现MVDR并对比两者的分辨率、旁瓣、抗噪性。对比表一放“为什么现在医用超声大多还在用DAS类的方法而不是MVDR”这个问题就能引出一个相当有深度的讨论——因为实时系统对计算量和稳健性的要求极高自适应方法虽然理论好但落地难。这个论点特别能体现你对工程问题的理解。6.3 光声/超声联合成像延伸另一个有意思的扩展方向是把DAS移植到光声成像中。光声成像的物理过程是短脉冲激光打到组织上组织吸收光能产生超声信号再用超声阵列接收重建图像。它的回波模型和超声不太一样光声信号是单程传播从目标到阵元所以延时计算去掉2倍距离其余逻辑基本不变。改一行核心代码就能跑光声的重建管线。做过超声DAS之后再去做光声DAS几乎是无缝切换这种能力会很受实验室导师喜欢。延时计算从2*d/c改成d/c声速换成对应介质速度再把光源的分布考虑进去一个基础的光声重建仿真就出来了。基于这个框架还能调整光通量分布和超声换能器的频率响应进一步与真实系统比对。6.4 关于声速不均匀问题的坦白最后说一个所有超声成像越不过去的话题人体组织声速并不是均匀的1540 m/s脂肪组织约1450 m/s骨骼约3000 m/s以上DAS建立在均匀声速假设上遇到组织不均匀图像就会变形、模糊。现在临床超声里好多伪像都源于此。你在项目总结里若把“DAS假设均匀声速”这个局限单独提出来并说清楚下一步可能的补偿方法是“逐层声速校正”那文章的专业层次立刻就不一样了。我个人在实际操作中的体会是DAS这个算法看起来简单但这套“先对齐后叠加”的逻辑是整个阵列信号处理的主心骨。把这个项目完整做下来你对延迟、加窗、插值、阵列几何空间关系的理解会有一个质的飞跃。最后再分享一个小技巧——连续调参数对比图像时一定要把窗口函数、动态范围、声速这些参数固化下来只改要研究的那一项否则很容易被视觉上的差异误导。这只是一个入门级波束形成器真正吃透它的人再去看自适应波束形成、平面波成像、合成孔径成像都会觉得一路畅通。