★ FEATURED ARTICLE 先想清楚一个事你在网上一搜“波束形成”出来的资料十有八九是通信方向讲阵列天线怎么对信号加权、怎么调零陷。但如果你要做的是医学超声成像、无损检测、声呐侧扫这类东西你遇到的“波束形成”大概率是另一副面孔——它不关心赋形和抗干扰它只关心一件事把每个阵元收到的回波信号对整齐然后加起来。这就是延迟求和业内缩写叫DASDelay And Sum。这篇文章我想把基于MATLAB的DAS波束形成实现从头到尾拆一遍。你可以把它看成一份“把原理翻译成可运行代码”的实操笔记涉及阵元时延怎么算、采样点怎么插值、图像怎么从一堆A线拼成B超帧以及我实际调试中踩过的那些坑。适合正在做超声成像仿真作业、刚进实验室要写波束形成代码、或者想从原理图过渡到MATLAB跑通第一版成像的同学。先说清楚它能干什么DAS是最简单、最稳健的波束形成器你把通道数据按几何时延对齐后求和就得到一条聚焦后的射频线再沿着成像区域扫描一轮就拼出一帧图像。虽然它牺牲了分辨率上限但计算效率高、实现门槛低而且它是后续DMAS、自适应波束、相干系数加权这些花活儿的基准线。把这个基础版本彻底吃透后面做任何改进都有底气。1. 为什么是DAS从阵列信号到一帧超声图像1.1 DAS到底在做什么很多同学第一次接触DAS会把它误认为“把N个通道的数据直接相加”。真这么做得到的东西基本没用因为超声回波不是同时到达每个阵元的。想象你在一个水池边扔了一块石头水面波纹传到池壁不同位置的时间不一样接收端如果想“听清”石头落水那一下就不能简单地把各位置的波形叠在一起而是要把每个通道的波形按到达时间差挪一挪让它们对齐后再叠加。数学上也很好理解。第i个阵元在t时刻收到的由成像点p散射回来的信号可以写成( s_i(t) a_i \cdot x(p, t-\tau_i(p)) n_i(t) )其中τ_i(p)就是声波从发射阵元到散射点p、再返回到第i个接收阵元的总传播时间。DAS做的就是估计出这个τ_i然后把s_i朝零时刻反向搬移( y(p, t) \sum_{i1}^{N} w_i \cdot s_i(t \tau_i(p)) )取t0那一下就是该像素点的波束输出。w_i是窗函数用来压低旁瓣常见的有汉宁窗、海明窗。为什么对齐后再相加有用因为目标点的回波是相干的对齐后同相叠加信号幅度涨N倍而噪声是随机的、部分相消的幅度只涨√N倍。这一涨一消就是DAS能提升信噪比和对比度的根源。1.2 时延计算的几何本质DAS的核心工作量不在“求和”在“算时延”。你得知道每个像素点对应到每个阵元信号要挪多少个采样点。这个时延只跟几何关系有关传导路径分两段发射段和接收段。如果阵列是相控阵发射时用的是聚焦延迟发射那你需要知道发射焦点坐标如果是平面波发射那发射波前是一条线或一个面时延计算方式又不一样。最常用的仿真和教学场景是单阵元发射、全阵元接收如合成孔径或者全阵元同时发射、全阵元接收如平面波。我下面讲代码时默认用的是“逐点动态聚焦接收 单焦点发射”的常规扫描模式这最容易对应到一条条A线A-line的生成过程。接收时延的计算公式是( \tau_{i,rec}(p) \frac{ \sqrt{(x_p - x_i)^2 y_p^2} }{ c } )这里(x_i,0)是第i个阵元的坐标p是成像区域内某个像素点c是声速。如果是发射焦点在(x_f, y_f)发射时延则算从发射阵元到焦点的距离再返程总时延是发射时延加接收时延。说得更直白一点每个像素点和每个阵元的连线的长度除以声速就是你该对齐的时间偏移量。这里面最容易被忽视的是“距离单位”和“数组索引”的换算。MATLAB里你算出来的是秒而射频数据的横轴是采样点序号必须除以采样率fs。这步一旦搞错轻则图像糊成一片重则整幅图全是斜纹。1.3 为什么MATLAB适合做这件事我没有贬低C或Python的意思但DAS这个算法确实是为MATLAB量身定做的场景。原因有三。第一DAS本质上是“循环累加”在MATLAB里可以用矩阵运算直接做广播一行代码把N个通道的时延插值全算完换成纯Python循环同样的规模可能要慢一个数量级。当然Python配NumPy也行但MATLAB的语法在高维数组广播上更顺手。第二MATLAB自带的插值函数和图像显示工具太省事了。延迟求和需要在每个通道的采样序列上做非整数点搬移sinc插值、线性插值、三次插值都有现成函数成像输出还方便直接调imagesc、imshow。第三调试体验。DAS这种算法诡异得很理论上全逻辑都对但图像就是不干净。MATLAB的变量可视化能一眼看出时延矩阵是不是出现了负值、插值后数组有没有NaN、是不是索引越界这些在命令行窗口里全都能实时观察。2. 仿真参数与成像区域怎么定2.1 阵元与孔径的选择逻辑做仿真不是参数越大越好。阵元数量N、阵元间距d、孔径D(N-1)*d这三个量决定了系统的横向分辨率、栅瓣水平还有计算量。先说阵元间距。这是一条硬规则阵元间距如果大于半波长λ/2方向响应会出现栅瓣也就是除了真实目标在别的方向也会“看到”一个虚影。真要说起来栅瓣位置满足( \sin\theta \sin\theta_0 m\frac{\lambda}{d} )当dλ/2时m取±1都可能落在物理区域内(−1到1之间)于是图像里就多出假目标。所以我的习惯是中心频率定下来先算波长然后把d取在0.5λ附近宁可密一点也别省钱。仿真里声速常取1540 m/s软组织声速中心频率如果是7.5 MHz波长远一点算λ c / f 1540 / (7.5e6) ≈ 0.205 mm所以阵元间距取0.1 mm到0.12 mm左右就非常安全。阵元数量N决定了横向分辨率。阵列的远场波束宽度大约正比于λ/D。D越大波束越窄但D受限于探头物理尺寸。仿真相控阵常用的64阵元、间距0.1 mm孔径也只有6.3 mm在深度40 mm附近横向分辨率是毫米级——不要指望它和光学图像一样锐利超声图像本来就要接受它的“涂抹感”。如果初始调试只是为了看算法是否正确16到32阵元完全够用跑得快才能快速迭代。2.2 声速、采样率与空间分辨率的换算做超声仿真的人最容易漏掉的一件事是把声速假设、采样率和图像深度绑在一起算。射频数据的时间轴长度取决于你要看多深成像深度10 cm那么单程传播时间是depth/c往返是双程所以回波最晚到达时间为( T_{max} 2 \times \frac{depth}{c} )采样率决定时间量化精度。如果fs 50 MHz而中心频率是7.5 MHz一个采样周期对应的时间是20 ns换算成距离就是( \Delta z c \times \Delta t 1540 \times 20e-9 \approx 0.031 \text{ mm} )也就是说时延计算时四舍五入到一个采样点轴向位置就会偏差30多微米。你说这对一帧图像重不重要单看似乎无所谓但叠加几百条扫描线后误差会在耗时延矩阵上累积出条纹。所以代码里时延矩阵一定要用浮点数存储做插值而不是取整后直接索引。轴向分辨率还有个经验值理论上约等于λ/2也就是约0.1 mm但实际受系统带宽影响往往只能做到0.20.3 mm。这个数字不用背知道数量级就行方便你验证仿真结果是不是“物理上合理”。2.3 关注动态范围从包络到对数压缩直接对DAS输出画图会得到一个噪点很多的波形因为射频信号是正负振荡的灰度图上看起来像毛刺。所以成图前必须做两件事包络检测和动态范围压缩。包络检测最常用的方法是希尔伯特变换对每一条A线做abs(hilbert(aline))。物理意义是把射频信号变成解析信号取模得到包络。这个操作会把射频振荡抹平留下慢变的幅度轮廓对应的就是组织反射强度。动态范围压缩是为了让人眼能看清。超声信号幅度动态范围经常超过60 dB而屏幕灰度只有8 bit约48 dB线性动态范围。所以我会用对数或者分贝压缩( I_{dB}(z) 20 \log_{10}( I(z) \epsilon ) )然后再映射到0~255的灰度区间。这步不做图像会变成“只有一个亮点四周全黑”——不是没有回波是被动态范围压扁了。3. 从零实现DASMATLAB代码与逐段解析3.1 第一步生成仿真射频数据没有实测数据时仿真生成RF数据是验证算法最可靠的路径。我习惯的做法是在成像区域里撒几个点目标每个点给一个反射系数然后让每个阵元接收由这些点散射回来的信号。每个阵元的回波是多个子波的叠加子波的包络可以用高斯脉冲近似或者更严谨地用带通高斯调制脉冲。核心代码如下注释里写了每个变量的物理含义% ---- 参数设置 ---- c 1540; % 声速单位 m/s fs 50e6; % 采样率单位 Hz depth 0.04; % 成像深度40 mm t 0:1/fs:2*depth/c;% 双程回波时间轴 N 64; % 阵元数量 pitch 0.1e-3; % 阵元间距 xele (0:N-1) * pitch;% 阵元x坐标 fc 7.5e6; % 中心频率 t_axis reshape(t, [], 1); % ---- 放置点目标 ---- pts [0.00, 0.015; % 深度15mm中心处 -0.003, 0.025; % 偏左3mm深度25mm 0.004, 0.035]; % 偏右4mm深度35mm amp [1, 0.8, 0.6]; % 反射系数 % ---- 脉冲形状 ---- impulse (tau) exp(-(tau/2e-7).^2) .* cos(2*pi*fc*tau);生成RF数据时要注意一点把每个点目标到达该阵元的时间tau 2*dist/c算出来后要用脉冲函数在对应时间位置采样并把每个点的回波累加在一起。rf zeros(length(t), N); for p 1:size(pts,1) xp pts(p,1); yp pts(p,2); for i 1:N dist sqrt((xp - xele(i))^2 yp^2); tau 2 * dist / c; % 双程往返时延 rf(:, i) rf(:, i) amp(p) * impulse(t - tau); end end这里我用了两层循环初学版本可读性优先。实际跑大数据集时这段可以改写成矩阵广播后面专门讲优化。3.2 第二步计算逐点时延接下来是DAS的核心。你要为成像区域里的每一个像素算它到每一个阵元的双程时延。先把成像网格建出来。假设轴向从0到depth横向从-6mm到6mm步长分别为0.2 mm和0.1 mm。x_range -0.006:0.0001:0.006; % 横向范围 z_range 0.0002:0.0002:0.04; % 轴向范围 [X, Z] meshgrid(x_range, z_range); X reshape(X, [], 1); Z reshape(Z, [], 1);时延矩阵的尺寸是像素数 × 阵元数。对每个阵元计算它与所有像素点的欧氏距离然后除以声速得到单程时间双程要乘以2numPix length(X); tauMatrix zeros(numPix, N); for i 1:N dists sqrt((X - xele(i)).^2 Z.^2); tauMatrix(:, i) 2 * dists / c; end这个矩阵每行代表“当前像素第i个阵元收到的回波延迟了多少秒”行与行之间代表的是不同的成像点。之后做延迟求和时要按这个矩阵将射频数据的采样点搬移。如果发射时是有聚焦的比如发射焦点在(0, F)那么发射时延不能忽略。一般公式( \tau_{total} \tau_{transmit} \tau_{receive} )但合成孔径或者平面波成像可以把发射时延做统一化处理。我建议初学时先做“全阵元同时发射、平面波前”的模型省掉发射时延把注意力放在接收DAS本身等接收过程跑通了再加发射聚焦否则两个变量在校验时互相干扰难排查问题。3.3 第三步延迟求和与插值取舍有了时延矩阵理论上只要把每列信号按对应时间偏移后累加。但坑在时间偏移量不一定是采样周期的整数倍。如果取整就引入量化误差。假设采样率50MHz时间量化误差为20ns换算下来是约0.031mm的误差。单次误差还无所谓但相邻像素点的误差不一致会产生垂直方向的条纹这正是很多初版DAS图噪难看的来源。所以这里一定要做插值。MATLAB有很多插值选项。最稳妥的、接近理论模型的是sinc插值。理想带限信号重建嘛。但sinc插值计算量大实际工程中线性插值就能得到很好的效果三次插值更好一点。我给出一种高效做法先对时间轴归一化到位移采样点序号然后用interp1做线性插值。为了避免interp1在边界产生NaN需要把目标索引限制在有效范围内。idx tauMatrix * fs; % 把秒换算成采样点索引 % rf矩阵尺寸为 LT x NLT是时间轴长度 lineCount length(x_range); nPixPerLine length(z_range); image zeros(nPixPerLine, lineCount); % 延迟求和主循环 for k 1:numPix row 0; for i 1:N targetIndex 1 idx(k, i); % MATLAB索引从1开始 if targetIndex 1 targetIndex length(t) % interp1需要单调递增的x这里直接用线性插值近似 weight 1.0; % 如果想加窗在这里乘上窗函数 row row weight * interp1(0:length(t)-1, rf(:, i), targetIndex-1, linear, 0); end end image(k) row; end这个三重循环写起来最直观但性能很差。优化思路是在循环外面把整列射频数据一次性插值到目标索引序列上aligned zeros(numPix, N); for i 1:N aligned(:, i) interp1(0:length(t)-1, rf(:, i), idx(:, i), linear, 0); end image sum(aligned, 2);这样瞬间清爽而且运行速度提升了不止一个量级。注意linear, 0后面的0表示超出边界时补零这比让interp1返回NaN再手动处理安全得多。3.4 第四步成像显示与量化评价拿到延迟求和的结果后直接reshape成图像矩阵还不行因为输出仍然是射频振荡信号。需要包络检测和对数压缩。imageRF reshape(image, nPixPerLine, lineCount); env abs(hilbert(imageRF)); % 对每一列做解析信号求包络 env env / max(env(:)); % 归一化 img 20 * log10(env 1e-12); % 分贝压缩 img img - max(img(:)); % 让最大值为0 dB img(img -60) -60; % 截断到60dB动态范围 imgGray mat2gray(img); % 映射到0~1 imagesc(x_range*1000, z_range*1000, imgGray); colormap(gray); xlabel(横向位置 (mm)); ylabel(深度 (mm)); axis image;观察结果时不要只看“有没有亮斑”。要判断点目标的位置是否准、横向涂抹是否严重、有没有栅瓣形成的假目标、深场区域的旁瓣有没有异常抬高。更严格的评价可以用点扩散函数PSF来测量-6dB宽度比如在图像里找到点目标位置画一条经过该点的横向剖面线量它的半高全宽FWHM。这是我每次仿真完必做的步骤。还有一种快速自检方法只放一个位于中心轴上的点目标然后看图像有没有左右对称性。DAS的几何模型本身是对称的如果图像不对称多半是阵列坐标或者时延索引写错了趁早排查。3.5 性能优化向量化、并行和数据降维仿真从32阵元换到128阵元、成像网格细化之后DAS的“像素数×阵元数×插值”三重压力会瞬间暴露。我给几个亲测有效的加速手段。第一用griddata或者纯向量化重写时延矩阵生成。前面for i 1:N里面用sqrt函数其实可以整块做矩阵广播让MATLAB自己处理内存对齐。第二把interp1换成自写的线性插值公式因为对等间隔数据做线性插值可以用floor和fraction两步完成省去函数调用开销。第三如果计算资源允许把最内层的阵元循环改成parfor或者用gpuArray把矩阵搬到GPU上算。不过插值类操作在GPU上未必稳定建议先确认结果和CPU完全一致再接GPU。我实际中受益最大的优化是不做全像素插值而是把时延索引的四舍五入误差控制在半个采样点内。具体做法是先对时延取整到最近的采样点然后用二次插值抛物线拟合补正剩余的小数部分。这个方式省掉了逐点调用插值函数但效果几乎和interp1持平。代价是代码复杂一点初学可以先不碰。4. 实战踩坑记录与问题排查速查表4.1 栅瓣是怎么出现的怎么避免栅瓣这个问题我之前调参时就踩过一次。当时为了压缩计算量把阵元间距从0.1mm改成了0.25mm心想着反正64个阵元反正仿真用点目标结果图像里每个真实点目标的左边和右边各多出两个对称的虚点。我当时第一反应是坐标反了排查了半天最后才发现是间距超过λ/2才会出现的栅瓣。记住一句话仿真里看不出“实测噪声掩盖栅瓣”的效应栅瓣会更明显因为点目标回波能量极集中。检查时看到对称虚假目标先算λ再算pitch与λ/2的关系通常一看就破案。4.2 时延变负、索引越界的处理还有一个经典问题在深度较浅、像素点位于阵列上方时像素到阵元的距离很小时延对应的索引可能是0甚至负数。interp1如果被问负索引要么报错要么返回NaN直接把整幅图像弄出一个“黑洞区”。解决办法是在计算时延后加一句idx(idx 0) 0;并统一用外插值为0的interp1保护数据。但要分清楚是“像素点在阵元上方”导致的无回波区还是“坐标系设置错误”导致的负时延。如果是前者应直接将该区域的像素标记为无数据而不是强行补零否则会在探头表面附近出现人为的暗带。4.3 图像条纹/网格伪影的根源如果你做出来的图像不是干净的点扩散而是满屏竖条纹或者斜纹90%出在时延计算和采样量化的匹配度上。我曾因为图省事把时延直接四舍五入成了整数索引结果图像边缘出现了明显的“阶梯感”。换成浮点索引加插值后立刻消失。另一种条纹来源很隐蔽时间轴和距离轴不一致。比如射频数据的深度范围和成像网格的z范围起点没对齐差那么几十个采样点看起来就像图像整体错位加条纹。排查方法很简单在代码里把“时间轴对应的物理深度”和“z_range”打印出来对比确保两者起点一致。4.4 计算太慢从矩阵化到分块处理仿真点数一多DAS最内层的循环可能跑几分钟。这时别再扣代码细节了先做分块。成像区域很大时我把z方向切成若干块每块独立算时延、插值和求和最后拼起来。好处是每块都能在内存里放得下也不会因为整个矩阵太大导致interp1内存爆炸。另外MATLAB的interp1在多列输入时维度处理不一样。如果是沿行插值需要转置或者用squeeze保证维度一致。这个坑我踩过好多次报错时看维度信息一眼能看出来。4.5 常见问题速查表现象可能原因快速排查与修复图像全黑或亮点极小对数压缩前未做包络检测abs(hilbert(...))后再log压缩点目标位置偏移严重声速设置不对或双程/单程忘记区分检查tau 2*dist/c中的系数2对称虚假目标阵元间距大于半波长计算λc/fc检查pitch图像斜纹/条纹时延取整、网格与时间轴不对齐改用浮点索引插值对齐深度起点图像左亮右暗窗函数方向搞反确认窗函数是按阵元序号排列不是按像素顺序生成图像颠倒图像矩阵reshape顺序不对记得先按z外层、x内层顺序reshape再用axis xy修正y轴方向5. 从DAS出发还能往哪走5.1 DMAS与自适应波束的简单扩展DAS跑通之后下一步最常见的改进是DMAS也就是延迟乘法求和。思路很简单延迟对齐之后不再是直接求和而是先两两相乘开方再求和。这样做的好处是增强相干信号抑制不相关的杂波和旁瓣图像对比度有明显提升。坏处是计算复杂度从O(N)变成O(N²)64个阵元就要计算64×644096次乘法仿真还行实时系统就要认真权衡了。用MATLAB实现DMAS时有个技巧可以绕过显式双重循环。延迟对齐后的矩阵aligned尺寸是numPix×N你只要构造一个“所有成对组合”的索引矩阵再用矩阵点乘。不过更简洁的方式是套现成的公式[ y_{DMAS}(p) \sum_{i1}^{N} \sum_{ji1}^{N} \text{sign}(s_i s_j) \sqrt{|s_i s_j|} ]这个操作在数值上很稳定而且保留符号可以避免负信号乘积的相位错乱。我第一次写DMAS的时候忘了开方之后的正负号处理出来的图像像是黑白底片折腾了好久才意识到DMAS的乘积累加是要带符号的。5.2 频域波束形成与相干补偿如果你觉得时域插值太麻烦另一个方向是把通道数据做FFT搬到频域用相移代替时移。延迟τ在频域里就是乘以exp(-j2πfτ)这样就不需要插值了而且可以用矩阵乘法一次性处理所有像素。但频域做法对信号带宽和参数一致性要求更高采样率低或距离轴短时误差反而比时域插值大。相干补偿则是另一个常用的增强手段。DAS把所有回波求和后可以计算每通道信号和总信号的相关系数然后用这个相关系数作为权重加权再求和。这个思想被用在不少低复杂度增强算法里效果明显、实现也不算复杂。我第一次试的时候几个不明显的小旁瓣确实被压下去了代价是图像整体变暗一点需要在对比度和亮度之间做平衡。5.3 最后的实战心得从原理到最后跑通一帧图像我最大的体会是DAS这种算法一小时能看懂一天能写出第一版但调试细节可能要花上一周。它不像深度学习那样有那么多指标可以刷它拼的就是对几何、时间、采样这三者关系的理解是否准确。你在MATLAB里看到的每一行代码背后都是物理模型如果图像不对先怀疑自己画坐标和算时间时偷的懒再去怀疑算法本身。做仿真时我建议你把验证拆成三个阶段先验证单通道的脉冲位置是否对再验证单像素点延迟求和后的时域波形是否有正确的峰值最后再拼整幅图。跳步直接跑到最后一步很容易被全局图像的假象带偏。这个习惯帮我避开了无数个“看起来差不多其实差很多”的bug。从DAS起步后面还有自适应波束、最小方差波束形成、频域相干叠加等一大票方向等着你去探索。但先把这篇里的时延、插值、包络、对数压缩吃透你会发现后面那些高级算法本质上都是在回答同一个问题怎么“对齐得更准、加得更聪明”。
阅读完成 · 觉得有帮助？
