简介这是一份面向声学、石油勘探及地球物理研究人员的MATLAB数值模拟案例围绕二维井孔声场问题结合有限差分法与完全匹配层PML技术实现边界吸收与高精度求解。资源包共1个文件为MATLAB脚本体量仅2KB代码精简但完整覆盖井孔几何参数与声学属性设置、差分网格划分、差分步长选取、PML边界构造及迭代求解流程并包含初始化、边界条件设定和结果可视化环节适合用于学习井声波传播特性、地层参数反演或作为有限差分数值模拟的入门参考。已有159人学习下载通过阅读作者提供的MATLAB源码可快速理解井孔声场模拟的核心步骤包括PML的吸收机制、时间步进策略以及声速分布与声场强度等结果的可视化输出对科研和工程实践均有实用价值。1. 井声场模拟中的边界反射与 PS008.m 的解决路径在二维井孔声场模拟里最让人头疼的不是波动方程本身而是计算域截断处的反射。井壁反射波、地层界面反射波需要传播到远处但有限网格必然截断空间如果截断边界处理不干净虚假反射会折返回井孔混在首波和后续波群里直接影响声速和衰减系数的提取。PS008.m 是这类问题的一个完整 MATLAB 实现它在网格外缘铺设完全匹配层Perfectly Matched Layer, PML把边界反射压到可忽略水平。对做声波测井、孔内声波探测或实验室超声模拟的人来说这份代码的价值在于能直接看到 PML 参数和差分数值模拟如何在二维井孔模型里配合工作而不是停留在公式层面。2. 波动方程离散化有限差分与完全匹配层的衔接2.1 二维声波方程的位移-应力形式井孔声场虽然在物理上是柱坐标下的轴对称传播但 PS008.m 这类平面切片模型通常直接按二维直角坐标处理对应 r-z 平面径向-深度上的声压场。控制方程用一阶速度-应力声波方程组表示∂p/∂t -ρc²(∂ux/∂x ∂uz/∂z)∂ux/∂t -(1/ρ) ∂p/∂x∂uz/∂t -(1/ρ) ∂p/∂z其中 p 是声压ux 和 uz 分别是 x、z 方向的质点振动速度ρ 是介质密度c 是声速。之所以用一阶方程组而不是直接对二阶声压方程做差分是因为交错网格上速度与声压天然错开半个网格步长空间二阶精度更容易保持且后续加入 PML 衰减项时每个分量方程都能单独修正。PS008.m 里的几个主数组就是按这种交错布局划分的p 数组定义在网格单元中心ux 数组在 x 方向与 p 错开半格uz 数组在 z 方向错开半格。代码里对 p、ux、uz 三个数组的更新顺序、差分长度都会遵循这一套布局改动网格尺寸时需要同步检查数组维度是否仍然匹配否则运行中会出现矩阵尺寸不兼容的报错。对于井孔问题井内是泥浆通常声速 1500 m/s 左右、密度 1000~1500 kg/m³井外是地层声速 3000~5000 m/s、密度 2000~3000 kg/m³两种介质的声阻抗差异很大。差分求解时交界面的参数赋值方式会直接影响波场的反射系数稍后在 3.3 节展开。2.2 完全匹配层的复数坐标伸缩原理PML 的理论基础是复数坐标伸缩。在吸收域内把空间坐标 x 替换成带复数的形式x x - (i/ω) ∫σx dx其中 ω 是角频率σx 是 x 方向上的衰减系数。平面波进入 PML 后振幅按 exp(-∫σx dx) 衰减。这里的关键是 σx 从 PML 内边界到外边界逐渐增大理想情况下不会产生数值反射。二维问题要在计算域四边都设置 PML角落区域还得同时做 x、z 两个方向的伸缩。实际代码里PML 的微分方程与内部区域共用同一套交错网格只是方程右侧多出一项衰减项。以速度分量为例∂ux/∂t σx ux -(1/ρ) ∂p/∂x离散化时这一项可以直接写成ux(n1) (ux(n) - (dt/ρ) ∂p/∂x) / (1 dt·σx)这种隐式形式的好处是保持差分格式稳定性PS008.m 中 PML 区域的速度更新几乎都能看到类似的分母项。内部区域 σx 0公式自动退化为普通差分因此完全可以用同一段代码处理两种区域靠衰减系数矩阵做区分。与传统吸收边界条件ABC相比PML 对宽入射角的波吸收效果稳定得多。ABC 在掠射角较大时反射明显而 PML 对大部分入射角都能把反射系数压制到 1% 以下这也是它在井孔声场这类需要长时间传播、多次反射的应用中成为首选的原因。2.3 PS008.m 中差分系数的组织方式把以上公式转成 MATLAB 代码时PS008.m 通常会把与位置相关的系数预先算好存放在与声压场尺寸一致的矩阵里。这样主循环中不需要逐网格判断介质类型直接查表索引即可。下面这段是系数组装的常见写法% 系数预计算常见实现方式 c2 c.^2; % 声速平方逐网格点 rho_face_x 0.5*(rho(1:end-1,:) rho(2:end,:)); % 界面密度平均 rho_face_z 0.5*(rho(:,1:end-1) rho(:,2:end)); Cu_x dt ./ rho_face_x; % ux 更新系数 Cu_z dt ./ rho_face_z; % uz 更新系数 Cp dt * c2 .* rho; % 声压更新系数 pmlx zeros(nx, nz); % x方向 PML 衰减系数矩阵 pmlz zeros(nx, nz); % z方向 PML 衰减系数矩阵这里的 Cu_x、Cu_z 与交错网格中 ux、uz 的实际维度对应。密度在界面处取算术平均是为了避免声阻抗突变导致差分项产生数值振荡如果直接给每个网格点赋单一介质密度跨界面差分会把界面当成分离的反射源模拟结果里会出现明显的虚假散射。pmlx、pmlz 先预留为全零矩阵只在第 3 章提到的 PML 层区域内赋值。这样做既保证内部区域不受衰减影响又能让所有数组维度在初始化阶段就对齐。调试时打印 pmlx 的部分行能直观确认衰减系数是否从内到外递增。下表总结了四个核心数组的约定格式数组常规尺寸更新方程中承担的角色pnx × nz声压场中心节点ux(nx1) × nzx 方向速度交错半格uznx × (nz1)z 方向速度交错半格pmlx / pmlznx × nzPML 衰减系数边界区域非零这类差分格式在 PS008.m 里体现得比较直观。对 401×401 网格四个数组的内存占用在 MATLAB 中约为 401×401×8×6 字节也就是不到 30 MB普通机器跑起来没有压力如果网格提升到 801×801内存占用会翻 4 倍就需要在初始化时用单精度或者压缩边界缓存。3. 网格步长与 PML 参数让 PS008.m 的骨架跑起来3.1 空间步长与 CFL 稳定性约束有限差分不是凭空选参数的。空间步长 dx 必须满足一个波长内至少 8~10 个网格点否则波动场的数值频散会掩盖真实信号。对井孔声场模拟dx 按下式估算dx c_min / (10 * fmax)其中 c_min 是模型中的最低声速fmax 是激励源频谱中的有效最高频率。比如井内泥浆声速 1500 m/s、声源中心频率 20 kHz频谱最高成分按 50 kHz 计算dx 就是 3 mm。时间步长再由 CFL 条件约束dt ≤ dx / (sqrt(2) * c_max)sqrt(2) 来自二维差分格式的稳定性上限c_max 取模型中最大声速。若地层声速 4500 m/s、dx 3 mmdt 约 4.7×10⁻⁷ s。PS008.m 里已经给出的 dx、dt 和网格数通常是按这组公式配好的改动 c_min 或 fmax 时需要重新核对这两个条件。一个常用的参数组合如下表参数示例取值说明nx, nz401 × 401总网格数含 PMLdx, dz0.003 m空间步长dt4×10⁻⁷ s时间步长npml20 层PML 厚度c_mud1500 m/s井内流体声速c_form3500~4500 m/s地层纵波声速fmax50 kHz激励源有效最高频率注意PML 层本身也要占网格点。20 层 PML 意味着每条边向内缩 20 个网格点四个角还要同时占用纵横两个方向的层区。如果模型的有效物理尺寸是 100 mm × 100 mm、dx 3 mm不加 PML 是 34×34 网格加完每边 20 层 PML 后总网格大约是 74×74。很多人在设置 nx、nz 时忘掉 PML 厚度导致实际有效区域比预期小后续分析波形时会发现直达波到达时间对不上。我一般先将有效计算域尺寸算好再加 2*npml 得到总网格数。3.2 PML 厚度与衰减系数梯度PML 厚度 npml 的经验范围是 10~30 层。少于 10 层时即使衰减系数很大波场在层内也来不及衰减就碰到截断边界多于 30 层对大多数井孔问题只是增加计算量。衰减系数的空间分布最常用二次抛物线σx(k) σmax · (k / npml)²其中 k 是距离 PML 内边界的层序号从 1 到 npml。σmax 的取值决定整体吸收强度理论反射系数约束下常用的经验公式为σmax -3 · c_max · ln(0.001) / (2 · npml · dx)这个公式把理论反射系数控制在 10⁻³ 量级。如果波形对比显示边界反射仍明显可以把抛物线指数从 2 改为 3或者把 σmax 乘 1.2重新跑对比实验。下面这段 MATLAB 代码展示了 PML 参数初始化过程% 完全匹配层参数初始化 npml 20; % PML 层数 sigma_max 3 * cmax / (2 * npml * dx) * log(1/1e-3); pmlx zeros(nx, nz); % x方向衰减系数矩阵 for k 0:npml-1 val sigma_max * ((k 0.5) / npml)^2; pmlx(k1, :) val; % 左边界区域 pmlx(nx-k, :) val; % 右边界区域 end代码里 (k 0.5)/npml 而不是 k/npml是为了把衰减系数放到交错网格的速度节点位置与声压节点错开半格。这样在差分离散时可以避免奇偶节点上的系数不匹配导致的高频振荡。如果你同时设置 z 方向的 pmlz逻辑完全一样只是沿第二维赋值循环。PML 区域内的 σx 与 σz 可以设置成不同大小。对于井孔模型如果深度方向z的波场传播距离比径向x更远或者想要重点观察沿井轴传播的导波可以适当加大 pmlz 的值让沿 z 方向漏出的能量更快衰减而 x 方向保持原有设置。3.3 井孔几何的网格映射井孔在网格模型里表现为一条竖直的低速通道。PS008.m 里常见的做法是设定井孔半径 r_well、井心位置 (xc, zc)计算每个网格点到井心的距离落在半径范围内的点标记为流体介质。半径在离散网格上用距离场判断可以避免逐网格写 if 条件的低效循环% 井孔几何映射矩阵化写法 x (0:nx-1) * dx; z (0:nz-1) * dz; [X, Z] meshgrid(x, z); r sqrt((X - xc).^2 (Z - zc).^2); inside_well r r_well; % 逻辑掩码 rho rho_form * ones(nx, nz); rho(inside_well) rho_mud; c(inside_well) c_mud;r_well 这里按实际物理半径除以 dx 后的网格数换算。若井眼半径 0.1 mdx 3 mm约 33 个网格点跨过 33 个点来模拟井壁的曲率精度足够。井壁处声阻抗差异大速度更新时的密度应做界面平均。我在 2.3 节给出的 rho_face_x 就是为跨界面差分准备的。如果 PS008.m 直接按网格点赋密度而没有界面平均模拟结果里井壁附近会出现一层周期性的数值振荡频率接近网格截断频率。判断方法很简单把井壁附近的波场快照放大若看到沿井壁排列的高频同心波纹多半是密度插值方式不对。边界条件方面外部截断边界由 PML 吸收不需要额外处理井壁如果是裸眼井流体-地层界面自然由阻抗差决定反射系数如果是套管井还需要在井壁网格上设置钢套管的声速和密度。PS008.m 的代码结构里通常会预留介质参数矩阵只需要把套管参数插入到 rho 和 c 的对应环带不需要改动差分主循环。4. 声源激励与迭代求解把 PS008.m 主循环跑通4.1 Ricker 子波激励源的加载井孔声场模拟中声源最常用的是 Ricker 子波时间域表达式为s(t) (1 - 2π²f₀²(t - 1/f₀)²) · exp(-π²f₀²(t - 1/f₀)²)f₀ 是中心频率。加 1/f₀ 的时间偏移是为了让子波在 t0 时接近零避免源值突变激发高频振荡。加载方式有两种把子波加到声压数组的源点位置或加到速度数组的源点位置。使用压力加载时模拟结果直接对应压力接收器测到的波形使用速度加载时首波到达时间会稍有不同。源点位置要选在内部区域且离开 PML 至少 5 个网格点否则子波的初始波前会被 PML 层内的异常系数污染。加载过程通常在时间循环内实现% 声源加载Ricker 子波中心频率 f0 f0 20000; % 20 kHz for n 1:nt t (n-1) * dt; tau pi * f0 * (t - 1/f0); src (1 - 2*tau^2) * exp(-tau^2); p(isrc, jsrc) p(isrc, jsrc) src * dt; % ... 下面是差分更新 p, ux, uz ... end每个时间步重新计算 src 是可读性较好的写法但 nt 到几万步时会增加不少耗时。预先分配 src_all zeros(nt,1)在循环前把整列子波算好循环内直接 p(isrc,jsrc) p(isrc,jsrc) src_all(n) * dt速度会快 20% 左右。对 PS008.m 这种课程设计级的代码保留逐次计算版本更容易对照公式理解调试完成后再优化即可。4.2 井壁流体-固体边界的处理井壁界面的处理方式决定了模拟结果能否反映真实的反射和透射。纯声学模型中井内流体和地层都当作声学介质处理井壁两侧只需要满足法向速度连续和声压连续。差分网格上这对应着速度分量的密度插值ux 和 uz 的更新系数在跨越井壁时使用界面两侧密度的平均值。这里有个常见的坑如果沿 x 方向跨越井壁时rho_face_x 用的是算术平均而两侧密度差超过 3 倍比如泥浆 1200 kg/m³石灰岩 2700 kg/m³建议改用调和平均 2ρ1ρ2/(ρ1ρ2)。调和平均对阻抗突变的保守程度更高能减少界面上的速度不连续。地层与井内流体的阻抗比越大两种平均方式的差异越明显。很多井声场模拟的实际问题是套管井。套管壁厚 5~10 mm在 3 mm 网格下只有 1~3 个网格点厚度无法准确模拟钢套管的板波和共振效应。常见做法是把套管的网格做细化或者把套管简化为一个高阻抗薄层用等效声速和密度填充这一圈网格。若 PS008.m 的目标是观察裸眼井地层波传播直接按流体-地层两种介质处理即可。4.3 PML 内部区域的时间推进循环时间循环是 PS008.m 的主体部分。以交错网格的一阶声波方程为例完整的推进过程如下% 主循环交替更新速度与声压交错网格 for n 1:nt % 更新 ux中心差分近似 dp/dx ux(2:end-1,:) ux(2:end-1,:) - Cu_x .* ... (p(2:end,:) - p(1:end-1,:)) / dx; % 更新 uz uz(:,2:end-1) uz(:,2:end-1) - Cu_z .* ... (p(:,2:end) - p(:,1:end-1)) / dz; % PML 区域附加衰减隐式形式 ux(2:end-1,:) ux(2:end-1,:) ./ ... (1 dt * pmlx(2:end-1,:)); uz(:,2:end-1) uz(:,2:end-1) ./ ... (1 dt * pmlz(:,2:end-1)); % 更新声压 p p p - Cp .* (... (ux(2:end,:) - ux(1:end-1,:)) / dx ... (uz(:,2:end) - uz(:,1:end-1)) / dz); % 源加载与接收点记录略 p(isrc, jsrc) p(isrc, jsrc) src_all(n) * dt; rec(n) p(irec, jrec); end这段代码的顺序是先在速度场上做空间差分再施加 PML 衰减最后更新声压。注意 ux(2:end-1,:) 的索引方式速度场定义在半个网格点上所以 p(2:end,:)-p(1:end-1,:) 实际对应 p 在相邻整点间的差分结果落在半网格点上。pmlx(2:end-1,:) 与 ux 的内部区域维度对齐加了衰减后如果速度异常减小说明 pmlx 的数值偏大需要调小 σmax。这个循环的主体是矩阵运算没有任何逐点 for 嵌套在 401×401 网格、6000 步迭代下 MATLAB 运行约几分钟。如果把 PML 更新单独抽出来、跳过内部区域的全零矩阵运算可以减少约 10%~20% 的耗时。对 PS008.m 这种带有教学性质的代码保持可读性优先并不建议一开始就做碎片化优化。5. 用反射系数定量验证 PML 吸收效果验证 PML 是否正常工作的办法不是看波场快照里的颜色深浅而是跑一组对照一个算例正常启用 PML另一个把 PML 区域直接置零当刚性边界其余参数完全一致。在井孔附近固定一个接收点记录两种设置下的波形。PML 有效时接收波形在初至之后会迅速衰减到接近背景噪声无 PML 时波形尾部会拖出一串周期性强的反射震相周期约等于波从接收点到边界再返回的时间。把两个波形叠加对比反射尾巴的差异一目了然。定量验证通常用反射系数。公式为R A_reflected / A_directA_direct 取直达波峰值A_reflected 取边界反射波到达时间窗内的最大幅值。下面这段脚本可以完成该计算% 反射系数两次模拟结果对比 load(wave_with_pml.mat); % 含 PML 的接收记录 load(wave_no_pml.mat); % 刚性边界记录 t_direct 50:150; % 直达波所在时间窗 t_refle 200:400; % 边界反射波时间窗 A_direct max(abs(wave_with_pml(t_direct))); A_reflected max(abs(wave_with_pml(t_refle))); R A_reflected / A_direct; fprintf(反射系数 R %.4f\n, R);如果 R 大于 0.05优先检查两处一是 σmax 是否偏小把抛物线指数从 2 提高到 3二是 npml 是否不够从 20 增到 30。修改后重新跑对照通常两次以内就能收敛到可接受量级。快速调试时还有个技巧先把模型缩小到 100×100 网格、300 步观察 PML 层内部的波场快照。如果 PML 层内出现明显强于内部区域的波动说明衰减系数没有正确施加。最常见的原因是 pmlx 与 ux 的维度不匹配——pmlx 是 nx×nz而 ux 是 (nx1)×nz索引 ux(2:end-1,:) 时与 pmlx 没有对齐会漏掉个别 PML 节点。沿 x 方向逐行打印 pmlx(1:5,:)确认系数从内到外递增就能定位到具体的赋值遗漏点。本文还有配套的精品资源点击获取