ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

水浸探头与相控阵声场仿真:从瑞利积分到多元高斯叠加

水浸探头与相控阵声场仿真:从瑞利积分到多元高斯叠加 简介面向超声检测与换能器设计人员的声场仿真脚本包聚焦水浸探头、相控阵超声与聚焦探头的声场模拟。通过MATLAB源码可快速查看不同介质与结构参数下的声压分布、声束形状和聚焦效果理解频率、直径、焦点位置等对检测性能的影响尤其适合声学仿真入门者及探头设计工程师用于参数预研。压缩包内共3个m文件包体仅3KB均为轻量级可运行脚本覆盖相控阵探头仿真、水浸聚焦探头仿真等典型模块结构清晰便于二次修改。已有271人浏览学习适合通过阅读源码掌握从模型建立到结果分析的完整思路在此基础上调整边界条件或探头参数即可快速评估不同方案降低试错成本提升超声检测系统开发或课题研究效率。1. 水浸探头和相控阵声场仿真先用哪个模型决定成败做超声检测的人大多经历过这种场景探头买回来规格书上写着“聚焦声束”“分辨率高”可真到了工件检测时焦点位置偏了、声束宽度对不上缺陷尺寸返工成本极高。与其反复上机试错不如先在仿真里把声场看明白。本文拆解的“探头声场仿真.zip”包含相控阵探头仿真.m、shuangjiezhigaosi.m、jujiaotantoushui.m三个核心脚本覆盖水浸平面探头、水浸聚焦探头和相控阵探头三类典型场景。它的技术核心不是简单调用pfield之类现成函数而是从瑞利-索末菲积分和多元高斯叠加模型Multi-Gaussian Beam ModelMGBM出发自己实现声压分布计算。这套代码适合无损检测工程师、超声探头设计人员和做相控阵成像算法研究的学生能直接帮你算出声束聚焦位置、焦柱直径和偏转角度把“经验调参”变成“参数可算”。2. 水浸探头声场仿真的数学基础从瑞利积分到多元高斯叠加2.1 瑞利-索末菲积分为什么在水浸场景下会失效一半水浸检测的特点在于声波先从水进入固体工件两种介质的声速、密度差异导致声束在界面处发生折射。瑞利-索末菲积分描述的是声源表面每个微元辐射的球面波在空间某点叠加的结果理论上只要知道界面两侧的声压和法向振速的连续性条件就能严格求解。但实际计算中积分面对应的是换能器表面和固液界面两层边界离散网格一旦加密计算量马上爆炸。尤其是计算水浸聚焦探头的三维声场时双界面积分涉及大量数值振荡收敛速度慢不适合做参数扫描。这也是代码里选多元高斯叠加模型的原因。MGBM 的核心思想是把活塞探头的表面振动近似为一组高斯函数的线性组合然后用高斯声束在分层介质中的传播解析解代替数值积分。它的适用前提是声束在传播过程中保持接近高斯分布形态且传播距离远大于探头半径。水浸检测条件下探头到工件的距离通常有较多波长满足这个前提。2.2 多元高斯叠加模型参数表和使用边界高斯叠加模型的精度取决于系数选取。常用的是一组 10 项或 15 项复系数分别对应不同阶数的拉盖尔-高斯展开。工程中最常见的是 Wen 和 Breazeale 提出的 10 阶解析近似其系数为系数值说明B11.0150.373第一项复振幅B2-0.282-0.593第二项修正A13.005-3.406第一项高斯宽度下降参数A20.908-2.267第二项高斯宽度参数A3-3.4736.593第三项高阶修正截断项数10 或 11项数增加对近场描述更好但远离声源时差异开始缩小shuangjiezhigaosi.m双高斯叠加的命名和这个模型直接对应。代码里用两组高斯函数叠加来模拟圆形活塞探头的辐射场计算流程分三步function [P, x, z] shuangjiezhigaosi(freq, radius, c, rho, Nx, Nz, xmax, zmax) % freq: 探头中心频率 Hz % radius: 探头半径 m % c: 介质声速 m/s % rho: 介质密度 kg/m^3 lambda c / freq; k 2 * pi / lambda; x linspace(-xmax, xmax, Nx); z linspace(0.001, zmax, Nz); % 避免 z0 处奇点 [X, Z] meshgrid(x, z); P zeros(Nz, Nx); % 10 阶高斯系数B1, B2, A1, A2, A3, A4, A5 coeff [1.015 - 0.373i, 0.282 0.593i, ... % B1 B2 3.005 - 3.406i, 0.908 - 2.267i, ... % A1 A2 -3.473 6.593i, 2.545 - 4.851i, ... % A3 A4 1.280 0.359i, -1.476 1.884i]; % A5 for m 1:8 % 每个高斯项独立传播后在观察点叠加 P P coeff(m) ./ (1 1i * coeff(2*m) * Z / (k * radius^2)) ... .* exp(-(X.^2 ./ (radius^2 * (1 1i * coeff(2*m) * Z / (k * radius^2)))) ... - 1i * k * Z); end P P * rho * c * (1 / radius^2); % 量纲归一化 end这段代码的主体逻辑是对每个高斯项先计算高斯宽度在传播距离 Z 上的复扩展量1 1i * A * Z / (k * radius^2)再乘以高斯横向分布指数最后在观察点叠加所有项。linspace(0.001, zmax, Nz)避开源平面上的奇点meshgrid直接生成二维坐标网格。参数radius的平方出现在分母上意味着频率越高、半径越大近场距离radius^2/lambda越长高斯项的相位变化越剧烈。2.3 声压归一化和量纲处理仿真结果通常需要归一化显示C 或 Python 里很多人直接取绝对值最大值归一但在 MATLAB 脚本里需要先判断是否需要保留相位信息。shuangjiezhigaosi.m中返回复数声压场P在实际使用时再做abs(P)。如果需要显示声压级分贝值公式为20*log10(abs(P)/max(abs(P(:))))。量纲上代码里最后乘rho*c的作用是把振动速度势转换成声压否则高频时幅值会差好几个数量级。3. 相控阵探头的声束偏转与聚焦延迟法则怎么和声场叠加结合3.1 相控阵仿真为什么不能用“单探头声场平移”来处理相控阵探头仿真.m和前两个脚本最大的差异在于声源模型变成了多阵元离散阵列。常见误区是把单个阵元的声场平移后叠加这在阵元间距接近半波长时误差不大但阵元间距一旦达到 1 倍波长以上栅瓣会从大角度方向冒出来平移叠加会把栅瓣位置算错。正确做法是每个阵元都是一个独立的子声源先计算它在空间中的辐射声场可以用多元高斯模型也可以用瑞利积分再乘上对应的延时相位后叠加。延时法则Delay Law决定声束的偏转和聚焦位置。相控阵相控阵的聚焦公式为delay_n (sqrt((x_f - x_n)^2 z_f^2) - sqrt(x_f^2 z_f^2)) / c其中x_n是第 n 个阵元中心横坐标(x_f, z_f)是焦点坐标。代码中这一步直接在阵元循环里完成先用pdist2之类函数计算阵元到焦点的距离再转成时间延迟。3.2 相控阵声场叠加代码实现下面这段是相控阵延时聚焦计算的常见结构和相控阵探头仿真.m的流程一致function [Ptotal, x, z] xiangkongzhen_focus(freq, pitch, Nelem, width, c, focal_z, Nx, Nz) % pitch: 阵元中心间距 m % Nelem: 阵元数量 % width: 阵元宽度 m % focal_z: 聚焦深度 m lambda c / freq; k 2 * pi / lambda; x linspace(-Nelem * pitch / 2, Nelem * pitch / 2, Nx); z linspace(0.001, focal_z * 1.5, Nz); [X, Z] meshgrid(x, z); Ptotal zeros(Nz, Nx); x_elem ((1:Nelem) - (Nelem 1) / 2) * pitch; % 阵元中心坐标 for n 1:Nelem % 目标点方向角与距离 r sqrt((X - x_elem(n)).^2 Z.^2); % 计算该阵元到聚焦点的距离差 rf sqrt((focal_z * 0 - x_elem(n))^2 focal_z^2); tau_n (rf - sqrt((x_elem(n) - 0)^2 focal_z^2)) / c; % 阵元辐射以一维活塞源近似 % 每个阵元在目标点的复声压 p_n exp(1i * k * (r - (Z - focal_z)) ); % 基尔霍夫近似 p_n p_n ./ sqrt(r); % 柱面波衰减 Ptotal Ptotal p_n .* exp(1i * 2 * pi * freq * tau_n); end Ptotal Ptotal / Nelem; end每个阵元都被视为一个等效的活塞子源子源的辐射在远场近似下呈柱面波衰减幅度按1/sqrt(r)变化相位按exp(i*k*r)变化。tau_n就是聚焦延时的量化值从阵元到焦点所需时间与参考点阵列中心处到焦点时间的差。叠加后Ptotal在焦点处同相叠加产生最大值偏离焦点则相位散乱互相抵消。参数上pitch不能超过2*lambda否则栅瓣进入可见空间范围。工程经验是阵元间距取0.6 ~ 0.8倍波长既能保证主瓣增益又能把栅瓣压到-12 dB以下。Nelem越多主瓣越窄但计算量和硬件成本也会上升。代码里最耗时的部分是meshgrid生成的大矩阵X、Z每个阵元循环都要对两个Nz×Nx矩阵做运算阵元数超过 32 时会有明显卡顿。3.3 角偏转仿真的小技巧用旋转坐标代替延迟计算中的三角函数需要声束偏转时很多人习惯在每个目标点都计算theta atan2(Z, X - x_elem(n))再算延迟量。这个写法逻辑直观但速度很慢。更高效的做法是直接在代码里把聚焦点坐标改成(z_f * sin(theta), z_f * cos(theta))延时公式不变。这样只需修改聚焦点位置一个变量就能同时实现聚焦深度和偏转角度的控制而不用每个阵元重新计算三角函数。4. 聚焦探头声场仿真参数怎么设焦点距离、孔径和频率的互相制约4.1 jujiaotantoushui.m 中的几何聚焦模型jujiaotantoushui.m对应的是水浸聚焦探头它在结构上比平面探头多了一个声透镜或凹面晶片。几何聚焦的核心是焦点位置由探头曲率半径决定但实际声场焦点和几何焦点并不重合——由于衍射效应声束在几何焦点附近会形成一个焦柱最大声压位置通常在几何焦点和近场末端之间偏移。代码实现中聚焦探头可以等效成相位延迟在探头表面按抛物线分布。生成相位分布时核心参数是焦点距离F和探头半径a聚焦相位为phi(r) k * (sqrt(F^2 r^2) - F)其中r是距探头中心的径向距离。仿真时把这个相位分布加到平面活塞探头的高斯模型上就能得到聚焦声场。% 聚焦圆盘声源抛物线相位分布 声压扫描 function P focused_disk(freq, radius, F, c, rho, N, z_axis) % radius: 晶片半径 % F: 几何焦距 % N: 径向采样点数 lambda c / freq; k 2 * pi / lambda; r_axis linspace(0, radius, N); % 只在径向取半剖面 phase_profile k * (sqrt(F^2 r_axis.^2) - F); % 抛物线近似 P zeros(size(z_axis)); for iz 1:length(z_axis) z0 z_axis(iz); % 每个径向位置上的环形面积加权叠加 integrand r_axis .* exp(1i * phase_profile) ... .* exp(1i * k * sqrt(z0^2 r_axis.^2)) ./ sqrt(z0^2 r_axis.^2); % 以梯形法完成径向积分 P(iz) trapz(r_axis, integrand); end P P * 2 * pi / lambda; end这段代码用一维径向积分代替二维面积分利用了轴对称特性。exp(1i * phase_profile)等效于声透镜的相位调制trapz是 MATLAB 自带的梯形积分函数数值稳定性比sum好。观察点在近场时sqrt(z0^2 r_axis.^2)变化迅速轴向采样间隔不能太大否则会漏掉振荡细节。4.2 水浸聚焦的关键参数表和调整原则参数取值范围影响调整倾向F 数焦距/直径F/1 ~ F/3比值越小聚焦越强焦柱越短薄工件选 F/3厚工件选 F/1.5频率2.25 ~ 15 MHz决定焦柱直径和穿透力晶粒粗大声衰减大时降频晶片直径6 ~ 25 mm决定焦点处声压增益空间受限时缩直径水中焦距根据工件表面距调整保证焦点落在工件内目标深度加上工件内声速修正水浸聚焦有一个容易忽略的点声束经过水/钢界面后焦点的实际位置会比按几何光学计算的深因为钢中声速约 5900 m/s远大于水中声速约 1480 m/s。实际计算焦点深度时要按斯涅尔定律重构不能直接用空气中的焦距数值。这也是代码仿真比实物调试有优势的地方——直接在z轴上扫描声压最大值就能找出真正的声学焦点位置。5. 排错与验证近场振荡、计算速度和栅瓣的实用检查手段5.1 从归一化声压图判断聚焦位置拿到shuangjiezhigaosi.m或jujiaotantoushui.m的输出结果后第一步不是看三维图而是先画轴向声压一维曲线。具体做法是取z轴上的声压绝对值的最大值位置对应就是声学焦点。如果轴向曲线出现双峰或主峰离几何焦点过远常见原因是聚焦相位分布中F值没有考虑透镜声速修正。判断方法[~, idx] max(abs(P(:, Nx/2)))用这一行代码直接定位焦点所在网格索引再换算实际深度。5.2 计算时间太长时的三个优化方向相控阵代码里最影响速度的是每个阵元生成声压矩阵时都要对全场格点做距离计算。第一个方向是把高斯的解析式直接写成二进制心算加速用bsxfun或implicit expansion替代显式meshgrid空间直接少一个大矩阵。第二个方向是把阵元分组具有相同延时变化率的阵元合并计算再乘以各自的相位旋转。第三个方向是只计算二维切面而不是完整三维体聚焦探头通常把半径和轴向两个方向都降成二维线阵时间能缩短到原来的 1/10。5.3 栅瓣怎么在代码里提前识别偏转声束的相控阵仿真结果里如果除主瓣外出现强度异常高的次峰先检查阵元间距pitch是否超过波长。若间距超限栅瓣必然出现在theta_g asin(lambda/pitch - sin(theta_main))方向。仿真时直接在代码里计算这个角度对比声压极值图里峰的位置能快速判断是栅瓣还是其他杂散信号。水浸场景下横波在工件中传播时会额外产生表面波这部分声场的回波并不代表探头发射声场本身验证时要先做时域截断处理只取直达波对应的时间窗。5.4 把仿真声场和实验扫查结果对齐的两个技巧第一个技巧是用水浸探头对着平底孔试块做 C 扫得到不同深度处的反射幅值曲线再和仿真轴向声压归一化曲线对比通常误差在 10% 以内就算模型成立。第二个技巧是聚焦探头关注焦点处焦柱直径实验中通过横越平底孔扫查得到 -6 dB 声束宽度仿真中则以abs(P) max(abs(P(:)))/2的区域宽度为准。两者对比时注意实验测到的是双程响应接收路径也经过一次声场调制严格比较需要把仿真单程声场做平方处理。本文还有配套的精品资源点击获取
返回列表