简介本资源是一套基于MATLAB实现的弹性波动方程有限差分法数值模拟代码包面向地球物理、声学仿真、结构动力学等领域的本科生、研究生及科研初学者用于理解波动方程离散化原理与数值求解实践。压缩包含30个文件以28个.m脚本为主涵盖初始化fdInitArray、边界处理fdInitBound、差分更新fdModUD、模型加载fdLoadModel、结果可视化plotcoloura/plotTraces等核心模块辅以2个.mat数据文件wncEdge.mat和wncUnity.mat提供预设介质参数与初始场整体仅46KB轻量易读便于逐行调试与教学演示。已有245人学习下载适合希望掌握时-空域中心差分离散、稳定性控制、位移场迭代更新及MATLAB高效数组运算的用户。读者可直接运行获得弹性波在均匀介质中的传播快照与波形记录配套函数命名规范、模块职责清晰是深入理解地震波正演建模与数值方法落地的优质入门范例。1. 这不是“跑个脚本就出图”的MATLAB练习——它是一套可复现、可调试、带物理边界的弹性波有限差分全链路实现你打开fdModUD.m发现它不依赖任何Toolbox连Signal Processing Toolbox都不调用却能稳定推进上千时间步你加载wncEdge.mat里面存的不是简单数组而是预计算的四阶精度边界插值权重矩阵你运行plotptwist2.m它输出的不是单帧位移快照而是带相位旋转校正的横波偏振轨迹动画。这不是教学示例而是一个面向地震波前建模真实需求打磨过的有限差分系统它把弹性波动方程的张量形式σ_ij C_ijkl ε_kl η ∂ε_kl/∂t拆解为纵波/横波耦合更新逻辑用fdInitMod2.m构建非均匀层状介质模型靠pedgeQuad2b.m实现二阶精度吸收边界最终通过fdSEGY2.m直接导出符合SEGY Rev1标准的二进制地震数据体。适合需要验证波场传播机理、调试各向异性参数敏感性、或为全波形反演FWI准备合成数据的地球物理建模者——尤其当你手头只有MATLAB基础环境R2018a及以上且必须避开GPU加速依赖时这套纯CPU向量化实现反而更可控、更易断点追踪。2. 弹性波动方程的离散化不是套公式从连续张量方程到MATLAB数组索引映射2.1 为什么必须用速度-应力双变量格式——避免数值频散与伪影的底层约束弹性波动方程在各向同性介质中常被简化为位移形式如摘要中给出的ρ∂²u/∂t² ∇·(c²∇u)但该形式在高波数区域会引入严重频散且难以自然嵌入自由表面、吸收边界等物理条件。本项目采用速度-应力交错网格staggered-grid velocity-stress FD将速度分量v_x、v_z和应力分量σ_xx、σ_zz、σ_xz分别定义在空间网格的不同偏移位置。其核心离散方程组为ρ ∂v_x/∂t ∂σ_xx/∂x ∂σ_xz/∂z ρ ∂v_z/∂t ∂σ_xz/∂x ∂σ_zz/∂z ∂σ_xx/∂t λ ∂v_x/∂x 2μ ∂v_x/∂x λ ∂v_z/∂z ∂σ_zz/∂t λ ∂v_z/∂z 2μ ∂v_z/∂z λ ∂v_x/∂x ∂σ_xz/∂t μ (∂v_x/∂z ∂v_z/∂x)其中λ、μ为Lamé常数由密度ρ和纵波速α、横波速β推导μ ρβ²λ ρα² − 2μ。这种格式天然满足应力-应变本构关系且中心差分在交错网格上能保持二阶精度显著抑制网格频散。MATLAB中不显式声明符号变量而是通过数组索引偏移实现物理量定位v_x存于(i,j)σ_xx存于(i0.5,j)σ_xz存于(i0.5,j0.5)—— 这直接决定了所有差分算子的索引偏移量。提示查看fdInitArray.m中vx,vz,sxx,szz,sxz的预分配尺寸。你会发现vx是Nx × Nz而sxx是(Nx1) × Nzsxz是(Nx1) × (Nz1)。这种尺寸设计正是交错网格的内存体现强行统一尺寸会导致边界越界或物理量错位。2.2 时间推进三步龙格-库塔法RK3替代显式欧拉——稳定性与精度的平衡点显式欧拉法u^{n1} u^n Δt·f(u^n)虽简单但CFL条件苛刻Δt ≤ 0.6·min(Δx,Δz)/max(α,β)且局部截断误差为O(Δt²)。本项目在fdModUD.m中采用三阶龙格-库塔法RK3其递推结构为% Step 1: k1 f(u^n) k1_vx fdComputeVxDeriv(sxx, sxz, rho, dx, dz); k1_vz fdComputeVzDeriv(sxz, szz, rho, dx, dz); k1_sxx fdComputeSxxDeriv(vx, vz, lam, mu, dx, dz); % ... 其他应力分量k1计算 % Step 2: u^(n1/2) u^n 0.5*Δt*k1; k2 f(u^(n1/2)) vx_half vx 0.5*dt*k1_vx; % ... 计算k2_vx, k2_vz等 % Step 3: u^(n1) u^n Δt*(2/3*k2 - 1/3*k1) vx vx dt*(2/3*k2_vx - 1/3*k1_vx); vz vz dt*(2/3*k2_vz - 1/3*k1_vz); sxx sxx dt*(2/3*k2_sxx - 1/3*k1_sxx); % ... 更新全部5个场变量该格式局部误差为O(Δt⁴)允许在相同CFL数下使用更大时间步长同时保持线性稳定性域足够覆盖典型地震波频带10–50 Hz。实际运行中fdInitGen.m会根据模型最大波速max_vel和最小网格间距min_dh自动计算推荐dt并写入结构体par.dt。2.2.1 RK3在MATLAB中的向量化实现关键避免for循环嵌套fdComputeVxDeriv函数不使用for i2:Nx-1, for j2:Nz-1而是利用MATLAB原生数组切片function dvxdt fdComputeVxDeriv(sxx, sxz, rho, dx, dz) % sxx: (Nx1) x Nz, sxz: (Nx1) x (Nz1), rho: Nx x Nz % 计算 ∂σ_xx/∂x ∂σ_xz/∂z → 结果尺寸 Nx x Nz d_sxx_dx (sxx(2:end,:) - sxx(1:end-1,:)) / dx; % 差分后尺寸 Nx x Nz d_sxz_dz (sxz(:,2:end) - sxz(:,1:end-1)) / dz; % 差分后尺寸 (Nx1) x Nz % 注意d_sxz_dz需在x方向取平均以匹配vx网格 d_sxz_dz_avg 0.5 * (d_sxz_dz(1:end-1,:) d_sxz_dz(2:end,:)); % → Nx x Nz dvxdt (d_sxx_dx d_sxz_dz_avg) ./ rho; end此处d_sxz_dz_avg的构造是关键因sxz定义在(i0.5,j0.5)其z方向差分结果位于(i0.5,j)需再沿x方向平均才能落到vx(i,j)网格点。这种索引对齐错误是初学者最常遇到的“波场发散”根源。3. 边界与源项从数学理想到物理可实现的工程落地3.1 吸收边界不是加个系数那么简单——pedgeQuad2b.m实现的二阶PML等效无限介质中的波传播需人工截断计算域传统海绵层sponge layer在大角度入射时反射率高。本项目采用二次多项式完美匹配层Quadratic PML其核心是将坐标拉伸为复数∂/∂x → (1 i·σ_x(x)/ω)⁻¹ ∂/∂x其中σ_x(x)在边界内按二次函数增长。pedgeQuad2b.m预计算了该PML在交错网格上的离散权重生成wncEdge.mat中的wedge_x,wedge_z矩阵尺寸与对应场变量一致。应用时只需在差分算子中乘以这些权重% 在fdModUD.m中更新应力时以sxx为例 d_vx_dx (vx(2:end,:) - vx(1:end-1,:)) / dx; sxx sxx dt * (lam .* d_vx_dx mu .* d_vx_dx) .* wedge_x; % wedge_x尺寸为(Nx1) x Nz自动广播匹配wedge_x在内部区域为1在PML层内从1平滑衰减至0.01衰减曲线由fdInitBound.m根据PML厚度par.pml_thick默认20格和最大衰减系数par.pml_max默认100生成。实测表明该PML在45°入射角下反射率低于−60 dB远优于简单指数衰减。3.2 震源注入力偶矩源moment tensor而非点脉冲——fdInitForce.m的物理建模逻辑地震震源本质是介质内力偶矩释放本项目支持六分量力偶矩张量M [Mxx Mzz Mxz Mzx Mzx Mzz]对称实际3个独立分量。fdInitForce.m将其离散为网格点上的等效力% Mxx分量作用于σ_xx方程∂σ_xx/∂t Mxx * δ(x-x0)δ(z-z0) * δ(t-t0) % 在MATLAB中用高斯型空间分布近似δ函数 [xg,zg] meshgrid(1:Nx,1:Nz); dist2 (xg - x0).^2 (zg - z0).^2; gauss exp(-dist2 / (2*sigma_x*sigma_z)); force_sxx Mxx * gauss * (1/(sigma_x*sigma_z*sqrt(pi))); % 注入到sxx场的对应位置考虑交错网格偏移 sxx(round(x0)1, round(z0)) sxx(round(x0)1, round(z0)) force_sxx;sigma_x,sigma_z由par.src_sigma控制默认0.5格确保源频谱主频与网格分辨率匹配避免混叠。对比简单点源力偶矩源能自然产生P波与S波的相对振幅比这对验证各向异性参数至关重要。3.2.1 检查源项是否生效三步验证法时域验证运行单步par.nt 1检查sxx,szz,sxz在源点附近是否出现预期符号模式如Mxx0时源点右侧sxx为正左侧为负频域验证对vx在源正上方接收点提取时序fft(vx_rec)应显示主频在1/(2*pi*par.src_sigma)附近能量守恒验证计算全域动能0.5*sum(rho.*(vx.^2vz.^2),all)与应变能0.5*sum((1/lam).*sxx.^2(1/lam).*szz.^2(1/mu).*sxz.^2,all)之和无耗散时应基本恒定允许浮点误差1e-12。4. 可视化与数据导出超越imagesc的物理意义还原4.1plotptwist2.m横波偏振轨迹的相位校正算法横波S波偏振方向携带介质各向异性信息但常规位移图imagesc(vx)或quiver无法直观呈现。plotptwist2.m实现瞬时相位旋转校正% 对vx, vz做Hilbert变换得到解析信号 vx_a hilbert(vx); vz_a hilbert(vz); % 计算瞬时相位角θ(t) atan2(imag(vz_a), imag(vx_a)) theta angle(vz_a ./ vx_a); % 将每个采样点的(vx,vz)绕原点旋转-θ得到径向分量vr和切向分量vt vr vx.*cos(theta) vz.*sin(theta); vt -vx.*sin(theta) vz.*cos(theta); % 绘制(vr, vt)随时间变化的轨迹Lissajous图 plot(vr(:), vt(:), .k, MarkerSize, 1); xlabel(Radial component); ylabel(Tangential component);该图中闭合椭圆表示线性偏振圆形表示圆偏振复杂形状表示非均匀各向异性。此方法无需预设参考相位直接从数据本身提取偏振演化是识别裂缝方向的关键步骤。4.2fdSEGY2.m生成工业级SEGY文件的字段填充逻辑SEGY格式要求严格字节对齐与字段语义。fdSEGY2.m不仅写入trace数据还填充关键二进制头Binary Header和扩展文本头Extended Textual Header字段名SEGY位置填充值来源物理意义SampleIntervalbytes 111–112round(par.dt*1e6)采样间隔微秒NumberOfSamplesbytes 115–116par.nt每道采样点数DataFormatCodebytes 117–1181IEEE浮点格式SourceXbytes 181–184par.src_x震源X坐标米GroupXbytes 189–192rec_x(i)第i个检波器X坐标特别注意fdSEGY2.m将vx和vz分别存为独立SEGY文件vx.segy,vz.segy并设置TraceHeader.TraceIdentificationCode 1CMP道和5垂直分量符合SEG标准。用户可用OpendTect或SeisSpace直接加载无需额外转换。4.2.1 快速验证SEGY完整性segy_info命令行工具在Linux/Mac终端执行需安装segyiopip install segyio python -c import segyio; fsegyio.open(vx.segy); print(f.bin[Samples]); print(len(f.trace[0]))输出应显示Samples字段值等于par.nt且首道长度与之完全一致。若不等说明fdSEGY2.m中fwrite的precision参数未设为float32导致字节错位。5. 调试与性能优化当波场“炸开”或“不动”时查什么5.1 波场爆炸数值不稳定的五级排查清单级别检查项命令/操作预期结果失败含义L1CFL数是否超限cfl par.dt * max(par.vel_p, par.vel_s) / min(par.dx, par.dz)cfl 0.6网格太粗或时间步太大重跑fdInitGen.mL2密度/波速是否为零或NaNany(isnan(par.rho(:)))any(par.rho0)L3边界权重是否全1max(abs(wedge_x(:)-1))1e-10wncEdge.mat未正确加载检查路径L4RK3中间步是否溢出在fdModUD.m中k1_vx计算后加assert(all(isfinite(k1_vx(:))))无报错某个差分算子除零如rho0L5内存是否碎片化memory命令查看PhysicalMemory.Available2GB多次运行未clear all重启MATLAB5.2 加速技巧从向量化到内存布局优化预分配所有中间数组fdInitArray1.m已完成但若修改模型尺寸务必重新运行该脚本避免MATLAB动态扩容禁用图形渲染运行前执行set(0,DefaultFigureVisible,off)plotcoloura.m等绘图函数将跳过屏幕绘制提速30%启用多线程BLAS在MATLAB命令行输入maxNumCompThreads(0)让底层线性代数库自动使用全部CPU核心避免parfor陷阱本项目未用parfor因FD更新存在严格时序依赖若强行并行会在sxx(i,j)更新时读取未完成的vx(i1,j)导致结果随机。注意fdCREWES.m包含一个隐藏的CPU亲和性设置feature(SetNumWorkerThreads, 4)适用于Linux服务器。Windows用户需注释此行否则MATLAB可能报错。最后若需快速生成测试快照直接运行fdInitGen; fdLoadModel; fdInitMod2; fdInitForce; for it1:100, fdModUD; end; plotcolourk(vx, vz, time, 100);这行命令组合将跳过所有初始化检查直奔第100步波场可视化——是验证安装完整性的黄金指令。本文还有配套的精品资源点击获取