简介这是一套面向地球物理、工程波动模拟初学者的初步虚谱法伪谱法MATLAB程序用于在复杂介质中模拟弹性波传播兼顾谱方法的高精度与有限差分式的直接求解适合地震波、声波和地下结构探测等应用场景。压缩包内共2个m文件整体仅3KB均为可直接运行的MATLAB源代码包含计算网格建立、材料参数设置、初始波场与边界条件配置、波动方程求解及结果可视化等基础功能模块。程序基于快速傅里叶变换FFT实现用户可按需调整网格密度、时间步长与物性参数从而适配不同研究目标。目前已有179人学习下载适合需要快速入手弹性波数值模拟的科研人员和工程师通过阅读和修改源码可进一步结合具体模型开展地震波传播、地下探测等深入模拟研究。1. 初步虚谱法程序弹性波模拟选伪谱法而不是差分法的关键理由做弹性波正演模拟时大多数人第一步会想到有限差分成熟、资料多、随手就能找到全套代码。但模型稍微大一点差分法的代价立刻显形——每个最小波长要放10到15个网格点三维模型一跑就是几天起步。伪谱法也叫虚谱法改用FFT在波数域里对空间求导一个正弦分量理论上两个网格点就能表示实际取4到5个点波场干净程度就能超过八阶差分这是它在弹性波模拟里最值钱的地方。这个“初步虚谱法程序”压缩包就是一条伪谱法弹性波正演的完整落地路径。下面按“原理→跑通→调参→避坑→验证”的顺序把这条路线讲透适合想用粗网格换高精度、又不想反复调数值频散的从业者。2. 伪谱法原理与弹性波方程离散为什么粗网格能换来高精度2.1 有限差分的分辨率瓶颈与伪谱法的替代思路伪谱法的本质是把空间导数的计算从网格局部挪到波数域全局。有限差分算子无论阶数多高本质上是对Taylor展开的截断。八阶差分在波数较低时接近理想导数一旦波数逼近Nyquist它的振幅响应就会明显偏离理想的ik——体现到波场里就是数值频散高频分量速度变慢或变快波前面出现拖着尾巴的振荡。要压住这种频散只有加密网格这一条路而加密网格意味着内存和计算量按模型维度的次方增长。伪谱法绕开了这个限制。它的做法是对波场做FFT正变换在波数域把每个谱分量乘上ik或者所需的任意阶导数算子再反变换回空间域。FFT对正弦分量是全精度的最大可表示波数就是Nyquist波数π/dx所以理论上每个波长两个网格点就能精确表示一个正弦波。实际模拟中取4到5个点/波长是为了照顾震源附近的奇异性和时间离散误差但已经比差分法少一半以上的网格。弹性波模拟尤其吃这个红利。模型里P波和S波速度差异明显Vp/Vs通常在根号二到根号三之间S波波长只有P波的一半左右。差分法为了保证S波不出频散整个网格都要按S波最短波长加密而伪谱法在最稀疏的网格上也能同时分辨两种波这是它在弹性波模拟里一直被保留的原因。对只需要做二维两层模型验证的场景来说这个优势更直接网格从300×300降到150×150内存少了四倍单步耗时也大幅下降。空间离散方式每波长网格点最大精确波数频散特征单步计算量二阶差分20~30有限强频散需极密网格小八阶差分10~15较高轻微频散中伪谱法4~5Nyquist无空间频散每次求导两次FFT顺带说一个检索层面的坑伪谱法还有个别名叫虚谱法二者都是pseudo-spectral的不同译法代码结构完全一致。看到“虚谱”别以为是另一个技术家族在文献和程序包里两个词混用的情况非常普遍。2.2 弹性波方程用一阶速度-应力形式写比二阶位移形式更顺手伪谱法可以作用在二阶位移方程上但工程上我更推荐一阶速度-应力方程组。原因有三个二阶方程里出现对x和z的混合二阶偏导伪谱法虽然也能算但边界条件和震源加载的物理意义不如一阶直观一阶方程里每个空间导数都是对单轴的代码结构规整不容易写错时间上可以直接用二阶中心差分做跳蛙递推存储量只有五个变量。方程写出来是下面这样五个未知量分别是水平速度vx、垂直速度vz以及三个应力分量σxx、σzz、σxzrho ∂vx/∂t ∂σxx/∂x ∂σxz/∂z rho ∂vz/∂t ∂σxz/∂x ∂σzz/∂z ∂σxx/∂t (λ2μ) ∂vx/∂x λ ∂vz/∂z ∂σzz/∂t λ ∂vx/∂x (λ2μ) ∂vz/∂z ∂σxz/∂t μ ∂vx/∂z μ ∂vz/∂xλ和μ是拉梅参数由Vp、Vs和密度换算λρ(Vp²−2Vs²)μρVs²。网格模型只要给每个点填上Vp、Vs、ρ三个量再逐点换算成λ和μ递推里需要的所有系数就齐了。这里有个容易踩的换算细节有些初步程序直接以λ2μ和μ的形式存参数省去每步除法有的则是每步都算。前者快很多后者代码易读但耗时。模拟前先确认参数文件里的“vp”“vs”“rho”是模型数组还是标量以及有没有做速度到拉梅参数的换算很多结果怪异的问题都出在这一步。时间递推用跳蛙格式即速度在n1/2时刻、应力在n时刻交错更新。它是二阶精度的空间误差由伪谱法控制在几乎为零时间误差就成了总误差的主要来源。如果要做长时间模拟可以换四阶Runge-Kutta但每步要算四次导数场成本高很多初步程序保持二阶中心差分即可。2.3 波数域求导算子整个伪谱法程序的核心就这一段把空间导数封装成一个函数后续所有递推都复用它。Python实现如下import numpy as np def spectral_derivative(field, dx, axis0): 沿指定轴对场做波数域一阶求导。 以二维波场形状 (nz, nx) 为准 axis0 对应 z 方向间距为 dzaxis1 对应 x 方向间距为 dx。 nx field.shape[axis] # 角波数向量fftfreq 返回频率索引乘 2*pi 后是角波数单位 rad/m k 2.0 * np.pi * np.fft.fftfreq(nx, ddx) # 把波数向量广播到 field 的目标轴 shape [1] * field.ndim shape[axis] nx k k.reshape(shape) # 正变换、在波数域乘 i*k、反变换取实部 derivative np.fft.ifft( np.fft.fft(field, axisaxis) * (1j * k), axisaxis ).real return derivative这段的要点有三个。第一fftfreq(nx, ddx)返回的频率索引从0到nx/2再到负半轴乘2π之后正好是角波数如果程序里FFT库返回的是循环频率而非角频率乘的因子要相应调整。第二乘的是1jk这是频域求导的傅里叶变换性质如果要求二阶导改成(1jk)**2即可伪谱法求高阶导数就是一次FFT的事这也是它区别于差分法的重要特性。第三反变换后必须取实部——由于浮点误差ifft会带回极小的虚部直接参与递推会被逐时间步放大最终污染整个波场。如果你拿到的是Fortran版本核心逻辑一模一样先调用FFT库做正变换把实数组转成复数谱乘上虚数单位乘波数再逆变换取实部。区别只在于FFT库的布局约定比如某些库返回的是物理排列的实部虚部需要先做fftshift数值实现不复杂但移植时最容易在这些地方翻车。3. 把初步虚谱法程序跑起来文件确认、环境准备与最小两层算例3.1 解压之后先确认四类文件缺了别急着跑一个典型的初步伪谱法程序包解压后通常包含四类东西主程序源码可能是Fortran的.f90、Python的.py或Matlab的.m参数定义要么是独立的文本/配置块要么写在主程序开头的常量区输出与绘图脚本把模拟结果写成二进制或文本的地震记录以及一个模型/算例目录。如果压缩包里带README先看README的“运行方式”一节那里会写明预期的输出文件名和物理单位。没有README是常态。我拿到这类包一般先按文件大小排个序最大的多半是结果或模型数据文件最小且能直接读的才是可执行入口。用编辑器打开主程序先搜“main”或“program”找到时间递推主循环的位置再搜“parameter”或“const”把网格尺寸、时间步长、震源位置这几组常量抄出来。这一步花十分钟后面能省下几小时的翻车排查。环境方面最常出现的坑是终端直接报“gfortran不是内部或外部命令”“conda不是内部或外部命令”这类信息。它的本质是编译器或Python解释器的路径没加入系统PATH而不是程序本身有问题。Windows下我建议统一装Anaconda并创建一个专门环境装好numpy和scipyFortran代码则用gfortran编译确保编译器和运行时库都是64位。32位和64位混用链接阶段大概率会报“无法定位程序输入点getcurrentpackagefullname”之类的动态库错误这类报错基本都和位数不匹配有关。3.2 最小两层模型一套立刻能用的参数为了验证程序能跑不用上来就上一个真模型我用一个两层介质模型上层2000m/s下层3000m/s横波速度按根号三比例对应。网格200×200网格间距10米震源用20Hz的Ricker子波、垂直集中力放在深度500米处。记录时长1.5秒时间步长0.5毫秒。参数值选取理由网格 nx×nz200×200两层模型只验证物理过程够用即可dxdz10 mS波最短波长约57.8m约5.8点/波长上层 Vp/Vs/ρ2000 / 1155 / 2000 kg/m³Vp/Vs√3接近真实沉积岩比例下层 Vp/Vs/ρ3000 / 1732 / 2200 kg/m³界面反射系数适中便于观察界面深度1000 m给反射波留出清晰的走时窗口震源Ricker20 Hz垂直集中力集中力同时激发P波和S波震源位置x1000 mz500 m离顶面和边界都足够远dt0.5 ms约为二维稳定极限的1/3偏保守记录长度1.5 s反射波有足够时间回到地表这里的关键是网格间距和震源主频的匹配。20Hz主频对应上层横波波长约57.8m10m网格每波长约5.8个点满足伪谱法4到5点的经验要求。如果把主频提到40Hz最短波长降一半网格间距就要缩到5m左右计算量翻四倍这个权衡在第4章还会展开。3.3 主循环跳蛙递推的顺序不能写反拿到程序后主循环通常是这样的结构我把它重写成一个尽量贴近各类初步程序的Python版本# 伪谱法弹性波模拟主循环跳蛙格式二阶时间差分 # 数组形状统一为 (nz, nx)axis0 是深度 zaxis1 是水平 x for it in range(nt): # 第一步由应力更新速度分量 vx dt / rho * ( spectral_derivative(sxx, dx, axis1) # ∂σxx/∂x spectral_derivative(sxz, dz, axis0) # ∂σxz/∂z ) vz dt / rho * ( spectral_derivative(sxz, dx, axis1) # ∂σxz/∂x spectral_derivative(szz, dz, axis0) # ∂σzz/∂z ) # 在震源位置加载垂直集中力源只加在 vz 分量 vz[nsz, nsx] dt / rho[nsz, nsx] * wavelet[it] # 第二步由速度更新应力分量 sxx dt * ( (lam 2.0 * mu) * spectral_derivative(vx, dx, axis1) lam * spectral_derivative(vz, dz, axis0) ) szz dt * ( lam * spectral_derivative(vx, dx, axis1) (lam 2.0 * mu) * spectral_derivative(vz, dz, axis0) ) sxz dt * mu * ( spectral_derivative(vx, dz, axis0) # ∂vx/∂z spectral_derivative(vz, dx, axis1) # ∂vz/∂x ) # 第三步应用吸收边界第4章展开 # 第四步在接收点处把 vx/vz 写入记录道注意这里的存储细节。vx代表水平振动速度vz代表垂直振动速度nsz是深度索引nsx是水平索引。加载垂直集中力时改的是vz而不是vx否则辐射图会绕着一个错误的轴转。如果震源是爆炸源则应该同时往sxx、szz、sxz上加各向同性压力而不是直接改速度分量——很多初步程序把爆炸源实现成“往所有点加同一个速度扰动”得到的结果看着有波但波型比例完全错误。时间递推的顺序是先更新速度再更新应力还是反过来其实可以互换只要震源加在正确的位置、并保持交错时刻的一致性。但每个时间步内部顺序要统一先算完所有速度分量再算所有应力分量不能混着来否则时间同步被打破高频成分会迅速失稳。上面的写法重在清晰效率不是最优。spectral_derivative每调用一次就是一次FFT加一次逆FFT这个循环里一共调用了12次其中对vx的x方向导数和vz的z方向导数在速度更新和应力更新里重复算了。优化时可以先把六个一阶导数场一次性算好再组装应力更新整体能省掉约1/3的FFT开销。初步程序不追求性能但这个逻辑值得记着后续做三维扩展时会用到。3.4 跑通后的第一道验收直达波与反射波的到达时间跑完之后先看接收器输出的两组记录。vz记录上第一个到达的是直达P波初走时约等于震源到接收点的距离除以上层纵波速度随后会看到来自界面的反射P波和反射转换波。如果vz上和vx上除了直达波外什么都没有检查震源类型和界面两侧波阻抗差——速度差太小也会让反射系数低到看不见这时加大两层速度比再试。一个快速的手工验算是把震源到界面的垂直距离和接收点的水平距离代入初等几何关系算出反射P波的走时再与程序输出的记录道对比。以第3.2节的参数为例震源深500m、界面在1000m、接收点水平距离100m时反射P波路径长约1503m按上层Vp2000m/s算走时约0.75秒直达P波走时约0.255秒。误差在1到2毫秒以内说明程序核心逻辑基本正确超过这个量就要回去检查网格方向或介质参数是否装反了。4. 三个必调参数时间步长、吸收边界与震源子波改错了就翻车4.1 时间步长伪谱法的稳定极限不是差分法那个公式伪谱法的空间导数没有频散误差但这不意味着可以无脑用大时间步长。如果时间差分仍然是二阶中心差分稳定性条件来自最大可表示的波数k_maxπ/dx与介质最大波速vmax的乘积。一维情况下理论极限约为0.637·dx/vmax二维时波数向量可以沿对角方向叠加k_max变为π√2/dx极限步长缩到约0.45·dx/vmax三维更严约0.37·dx/vmax。伪谱法能精确表示到Nyquist波数而差分法在高波数部分的振幅响应实际上是衰减的相当于天然滤掉了一部分不稳定成分所以伪谱法对时间步长更敏感。我一般不会顶着极限值用而是取二维极限的一半左右dt 0.3·dx/vmax。这样既留出安全余量又不会因为步长太小让长时程模拟的步数猛增。以第3章那个两层模型为例vmax取下层纵波速3000m/sdx10m二维稳定极限约1.5毫秒取0.5毫秒是极限的1/3属于稳妥选择。如果压缩包代码里时间步长是写死的先按这个公式重新算一遍再跑。判断步长是否过大不一定要等波场爆炸。最快的诊断方法是打印每一时间步的总能量在均匀无吸收模型里总能量应当基本守恒。如果看到某个分量能量随步数单调上升比如从1e-2涨到1e0基本可以断定步长越过稳定极限。把dt缩小到原来的1/4再跑能量曲线趋于平稳就说明问题出在此处而非程序逻辑。提示步长的大小对伪谱法的影响是“全有或全无”的越界一步就会在几十步内爆掉。养成每个新模型先跑50步看能量的习惯比跑完整个记录才发现翻车要省时得多。4.2 吸收边界阻尼带的厚度和衰减系数要一起调初步程序很少带PML最常见的是在计算域四周加一层阻尼带也叫海绵边界或吸收层。它的原理很简单每时间步对边界区的波场乘一个小于1的衰减因子让波在到达人工边界前衰减到可忽略。实现不难但参数配不对时阻尼带本身就会变成反射源效果比不加还糟。阻尼系数一般取成空间位置的函数例如σ(x)σ_max·(x/L)²其中L是阻尼带的网格数x是该点到计算域边界的归一化距离。σ_max的经验范围是2到3倍的vmax/(L·dx)。L的取值至少要覆盖一个中心波长中心波长用震源主频对应的波长来算λ_cvmax/f0。在20Hz主频、3000m/s最大速度的模型里中心波长150米L建议取15到20个网格dx10m时。L太薄时波在阻尼带内还没衰减到位就撞到硬边界反射能量依旧可观。给一段阻尼带实现可以直接替换第3.3节主循环里的“第三步”# 生成二维阻尼衰减系数场四个边界各加 L 个网格 def build_damper(nz, nx, L, vmax, dt): sig_max 3.0 * vmax / (L * dx) # 单位 1/sL*dx 是带的总长度米 damp np.ones((nz, nx), dtypenp.float64) for i in range(L): factor sig_max * ((i 1) / L) ** 2 * dt damp[i, :] * np.exp(-factor) # 上边界 damp[-(i 1), :] * np.exp(-factor) # 下边界 damp[:, i] * np.exp(-factor) # 左边界 damp[:, -(i 1)] * np.exp(-factor) # 右边界 return damp # 每个时间步在递推之后执行 vx * damp vz * damp sxx * damp szz * damp sxz * damp注意角点区域会被重复衰减这个实现在角点的衰减系数比边上大一倍实际影响不大如果要严格处理需要按到最近边界的距离分别计算x和z方向的衰减因子再相乘。更重要的是阻尼带内最好保持常数速度模型不要放界面或强速度梯度否则波在带内产生反射这部分反射同样会污染内部波场。4.3 震源子波Ricker子波的主频和网格间距是配对关系震源子波最常用Ricker表达式是f(t)(1−2π²f₀²(t−t₀)²)exp(−π²f₀²(t−t₀)²)其中t₀一般取1.2到1.5个主频周期让子波初始时刻接近零避免在t0时刻给波场一个阶跃激励。实现如下# Ricker 子波f0 为主频dt 为时间步长 t np.arange(nt) * dt t0 1.2 / f0 wavelet (1.0 - 2.0 * (np.pi * f0 * (t - t0)) ** 2) * \ np.exp(-(np.pi * f0 * (t - t0)) ** 2)主频f₀越高波场分辨率越高能分辨更薄的层但代价是S波最短波长同步变短需要更细的网格。经验约束是每个最短波长至少要有4到5个网格点即dx ≤ v_s_min/(4·f₀)。这里速度取整个模型里最小的S波速度因为S波波长最短最容易频散。以第3章模型为例上层Vs1155m/sf₀20Hz时最短波长约57.8mdx10m相当于每波长约5.8个点处于安全区间。如果把主频从20Hz提到40Hz最短波长降一半dx就必须缩到5m左右计算量涨四倍这就是主频和网格步长的直接权衡。如果压缩包默认震源是爆炸源而你需要同时看P波和S波换成垂直集中力源即可。爆炸源只会辐射纯纵波无论后来怎么调吸收边界和网格横波分量始终是零这一点在验证环节最容易把人带偏。震源加载位置建议离边界至少10个网格否则即使有阻尼带源与人工边界之间的多次反射也会干扰早期波场。5. 伪谱法程序避坑指南5个最常见的翻车现场与排查方法5.1 波场图上一片棋盘格噪声高频Nyquist分量在作怪现象模拟几步后波场图出现颗粒状交替亮暗的棋盘格尤其在震源附近最明显振幅随步数增长。原因单点加载震源在空间上是一个极窄的尖峰它的频谱在Nyquist波数附近仍然有可观的能量。伪谱法对这个分量是全精度放大的不像差分法有天然的抑制于是波场里出现以单个网格为周期的交替扰动视觉上就是棋盘格。解决把震源先做空间平滑再乘子波。常见做法是给震源区一个高斯半径比如σ_source1.5倍的dx让源在空间上分布到8到10个网格点同时检查FFT后是否取了实部虚部残留也会产生类似的高频噪声。如果程序本身没有平滑函数可以在加载震源前对相邻网格按高斯权重分配能量。5.2 边界反射比预期早出现阻尼带没盖住最大波长现象波场图上在计算域边界附近出现强反射弧反射波到达内部接收点的时间明显早于模型里真实界面的理论走时。原因阻尼带厚度L没有按最大中心波长设计。L太薄时长波长成分在带内衰减不够振幅在到达硬边界时仍然可观边界反射自然回传。解决把L加大到至少一个中心波长。用vmax/f0算出中心波长后再换算成网格数如果程序里阻尼带厚度写死改参数或预处理速度模型时把边界区扩展。验证方法是给一个无反射界面的均匀模型跑一次把接收点能量画成时间曲线观察末段是否有明显长时间拖尾的反射能量。阻尼带的σ_max也要同步调到2到3倍vmax/(L·dx)薄带配大衰减、厚带配小衰减两种组合效果不同需要交叉验证。5.3 振幅随时间指数增长直到NaN时间步长越过稳定极限现象前面的波形看着正常到几百步之后某个应力分量量级从1e-2跳到1e20甚至直接变成NaN程序挂掉。原因按照4.1节算出的单方向稳定条件只是一维理论在二维模型里波动能量沿多个方向传播实际允许的步长通常更小。很多初步程序的dt是作者用他的模型试出来的换到你自己的网格尺寸和速度模型后稳定余量可能已经不够。解决把dt缩小到当前值的一半甚至1/4重跑看是否仍然发散。同时建议在时间循环里加一个能量检测每50步打印一次波场总能量看到指数上升就立即终止避免跑完整个记录长度才发现翻车、白烧算力。稳定步长与dx、vmax的具体取值参考4.1的公式但最终以你的模型能量曲线为准这是这类程序最不可省的一步基本功。5.4 横波分量离奇失踪震源类型和参数化把S波灭掉了现象接收记录上只有纵波初至之后全是微弱的低频尾巴理论上应当明显的反射转换波消失vx分量尤其干净。原因两类常见误操作。一是用爆炸源加载它只激发P波S波天然为零二是参数换算时把μ设成了0或很小的值导致S波速度接近0波场根本传播不出去。解决换成垂直集中力源加载在vz分量上同时检查拉梅参数换算μρVs²如果模型文件里Vs列填了0或没填μ就会变成0。一张快速自检图是把Vp、Vs画成按深度的曲线看Vs站点是否与Vp同步变化若Vs全程为0程序里再聪明也算不出S波。5.5 程序在Windows下报动态库或命令找不到环境没有对齐现象终端执行编译命令时报“gfortran不是内部或外部命令”运行Python时报“numpy模块不存在”或者程序启动直接报“无法定位程序输入点getcurrentpackagefullname于动态链接库…”运行就中断。原因三类问题混在一起——编译器或解释器的PATH没有配好、Python环境不对、以及32位/64位运行时库混用。后者在下载了旧版编译好的现成程序包时最容易出现因为动态链接库的位数和主程序不匹配系统加载时就报找不到入口点。解决Fortran源码重新用本地gfortran编译别直接用网上别人编好的exePython部分统一到Anaconda的64位环境建环境后执行conda install numpy scipy别用系统自带的Python。检查位数的方法是打开终端分别敲gfortran --version和python --version确认输出里有没有带32位字样。这一类报错的排查逻辑和网上常见的“conda不是内部或外部命令”完全一样先确认环境变量再确认位数最后才是代码问题。6. 验证伪谱法程序正确性解析解对比与网格收敛性检查写完代码、跑通模拟不等于程序是对的。我验证任何正演程序都走固定的三步解析解走时对比、网格收敛性检查和能量守恒检查。这三步能过滤掉九成以上的隐性错误。第一步用两层介质模型或均匀半空间模型把接收点的波场与解析走时对比。均匀半空间里直达P波走时是r/Vp直达S波走时是r/Vs两层模型里反射P波走时按镜像源法计算公式简单手算即可。把程序输出的单道记录拆成vx和vz两列找到初至时间误差在1到2毫秒内算通过。严格检查可以再加一个垂直自由表面边界对比Rayleigh波存在与否但初步程序一般不需要。第二步是网格收敛性检验。把dx、dz同时减半dt等比缩小重跑同一个模型对比同一接收点的波形。伪谱法如果实现正确两次结果的波形差异应该在1%以内且差值主要集中在高频尾部。如果减半网格后波形明显变化说明原网格本身就不满足分辨率要求需要按第4章的公式重新选择网格间距而不是程序逻辑有问题。第三步是能量监测这个前面提过。在没有阻尼带和震源持续加载的均匀模型中总能量应该守恒在带阻尼带的模型中能量应单调衰减而不是振荡上升。把每步总能量画出来曲线形状正常程序才算真正通过验收。我拿到的每一个伪谱法程序都会先跑这三步再做物理实验。走时对不上先查震源类型能量发散了先查时间步长波形不收敛先查网格间距顺序不要倒过来。这个习惯帮我挡掉了大量“看起来正常其实参数错位”的翻车现场。希望帮到你。本文还有配套的精品资源点击获取 伪谱法弹性波正演模拟：从原理到避坑实战指南
简介这是一套面向地球物理、工程波动模拟初学者的初步虚谱法伪谱法MATLAB程序用于在复杂介质中模拟弹性波传播兼顾谱方法的高精度与有限差分式的直接求解适合地震波、声波和地下结构探测等应用场景。压缩包内共2个m文件整体仅3KB均为可直接运行的MATLAB源代码包含计算网格建立、材料参数设置、初始波场与边界条件配置、波动方程求解及结果可视化等基础功能模块。程序基于快速傅里叶变换FFT实现用户可按需调整网格密度、时间步长与物性参数从而适配不同研究目标。目前已有179人学习下载适合需要快速入手弹性波数值模拟的科研人员和工程师通过阅读和修改源码可进一步结合具体模型开展地震波传播、地下探测等深入模拟研究。1. 初步虚谱法程序弹性波模拟选伪谱法而不是差分法的关键理由做弹性波正演模拟时大多数人第一步会想到有限差分成熟、资料多、随手就能找到全套代码。但模型稍微大一点差分法的代价立刻显形——每个最小波长要放10到15个网格点三维模型一跑就是几天起步。伪谱法也叫虚谱法改用FFT在波数域里对空间求导一个正弦分量理论上两个网格点就能表示实际取4到5个点波场干净程度就能超过八阶差分这是它在弹性波模拟里最值钱的地方。这个“初步虚谱法程序”压缩包就是一条伪谱法弹性波正演的完整落地路径。下面按“原理→跑通→调参→避坑→验证”的顺序把这条路线讲透适合想用粗网格换高精度、又不想反复调数值频散的从业者。2. 伪谱法原理与弹性波方程离散为什么粗网格能换来高精度2.1 有限差分的分辨率瓶颈与伪谱法的替代思路伪谱法的本质是把空间导数的计算从网格局部挪到波数域全局。有限差分算子无论阶数多高本质上是对Taylor展开的截断。八阶差分在波数较低时接近理想导数一旦波数逼近Nyquist它的振幅响应就会明显偏离理想的ik——体现到波场里就是数值频散高频分量速度变慢或变快波前面出现拖着尾巴的振荡。要压住这种频散只有加密网格这一条路而加密网格意味着内存和计算量按模型维度的次方增长。伪谱法绕开了这个限制。它的做法是对波场做FFT正变换在波数域把每个谱分量乘上ik或者所需的任意阶导数算子再反变换回空间域。FFT对正弦分量是全精度的最大可表示波数就是Nyquist波数π/dx所以理论上每个波长两个网格点就能精确表示一个正弦波。实际模拟中取4到5个点/波长是为了照顾震源附近的奇异性和时间离散误差但已经比差分法少一半以上的网格。弹性波模拟尤其吃这个红利。模型里P波和S波速度差异明显Vp/Vs通常在根号二到根号三之间S波波长只有P波的一半左右。差分法为了保证S波不出频散整个网格都要按S波最短波长加密而伪谱法在最稀疏的网格上也能同时分辨两种波这是它在弹性波模拟里一直被保留的原因。对只需要做二维两层模型验证的场景来说这个优势更直接网格从300×300降到150×150内存少了四倍单步耗时也大幅下降。空间离散方式每波长网格点最大精确波数频散特征单步计算量二阶差分20~30有限强频散需极密网格小八阶差分10~15较高轻微频散中伪谱法4~5Nyquist无空间频散每次求导两次FFT顺带说一个检索层面的坑伪谱法还有个别名叫虚谱法二者都是pseudo-spectral的不同译法代码结构完全一致。看到“虚谱”别以为是另一个技术家族在文献和程序包里两个词混用的情况非常普遍。2.2 弹性波方程用一阶速度-应力形式写比二阶位移形式更顺手伪谱法可以作用在二阶位移方程上但工程上我更推荐一阶速度-应力方程组。原因有三个二阶方程里出现对x和z的混合二阶偏导伪谱法虽然也能算但边界条件和震源加载的物理意义不如一阶直观一阶方程里每个空间导数都是对单轴的代码结构规整不容易写错时间上可以直接用二阶中心差分做跳蛙递推存储量只有五个变量。方程写出来是下面这样五个未知量分别是水平速度vx、垂直速度vz以及三个应力分量σxx、σzz、σxzrho ∂vx/∂t ∂σxx/∂x ∂σxz/∂z rho ∂vz/∂t ∂σxz/∂x ∂σzz/∂z ∂σxx/∂t (λ2μ) ∂vx/∂x λ ∂vz/∂z ∂σzz/∂t λ ∂vx/∂x (λ2μ) ∂vz/∂z ∂σxz/∂t μ ∂vx/∂z μ ∂vz/∂xλ和μ是拉梅参数由Vp、Vs和密度换算λρ(Vp²−2Vs²)μρVs²。网格模型只要给每个点填上Vp、Vs、ρ三个量再逐点换算成λ和μ递推里需要的所有系数就齐了。这里有个容易踩的换算细节有些初步程序直接以λ2μ和μ的形式存参数省去每步除法有的则是每步都算。前者快很多后者代码易读但耗时。模拟前先确认参数文件里的“vp”“vs”“rho”是模型数组还是标量以及有没有做速度到拉梅参数的换算很多结果怪异的问题都出在这一步。时间递推用跳蛙格式即速度在n1/2时刻、应力在n时刻交错更新。它是二阶精度的空间误差由伪谱法控制在几乎为零时间误差就成了总误差的主要来源。如果要做长时间模拟可以换四阶Runge-Kutta但每步要算四次导数场成本高很多初步程序保持二阶中心差分即可。2.3 波数域求导算子整个伪谱法程序的核心就这一段把空间导数封装成一个函数后续所有递推都复用它。Python实现如下import numpy as np def spectral_derivative(field, dx, axis0): 沿指定轴对场做波数域一阶求导。 以二维波场形状 (nz, nx) 为准 axis0 对应 z 方向间距为 dzaxis1 对应 x 方向间距为 dx。 nx field.shape[axis] # 角波数向量fftfreq 返回频率索引乘 2*pi 后是角波数单位 rad/m k 2.0 * np.pi * np.fft.fftfreq(nx, ddx) # 把波数向量广播到 field 的目标轴 shape [1] * field.ndim shape[axis] nx k k.reshape(shape) # 正变换、在波数域乘 i*k、反变换取实部 derivative np.fft.ifft( np.fft.fft(field, axisaxis) * (1j * k), axisaxis ).real return derivative这段的要点有三个。第一fftfreq(nx, ddx)返回的频率索引从0到nx/2再到负半轴乘2π之后正好是角波数如果程序里FFT库返回的是循环频率而非角频率乘的因子要相应调整。第二乘的是1jk这是频域求导的傅里叶变换性质如果要求二阶导改成(1jk)**2即可伪谱法求高阶导数就是一次FFT的事这也是它区别于差分法的重要特性。第三反变换后必须取实部——由于浮点误差ifft会带回极小的虚部直接参与递推会被逐时间步放大最终污染整个波场。如果你拿到的是Fortran版本核心逻辑一模一样先调用FFT库做正变换把实数组转成复数谱乘上虚数单位乘波数再逆变换取实部。区别只在于FFT库的布局约定比如某些库返回的是物理排列的实部虚部需要先做fftshift数值实现不复杂但移植时最容易在这些地方翻车。3. 把初步虚谱法程序跑起来文件确认、环境准备与最小两层算例3.1 解压之后先确认四类文件缺了别急着跑一个典型的初步伪谱法程序包解压后通常包含四类东西主程序源码可能是Fortran的.f90、Python的.py或Matlab的.m参数定义要么是独立的文本/配置块要么写在主程序开头的常量区输出与绘图脚本把模拟结果写成二进制或文本的地震记录以及一个模型/算例目录。如果压缩包里带README先看README的“运行方式”一节那里会写明预期的输出文件名和物理单位。没有README是常态。我拿到这类包一般先按文件大小排个序最大的多半是结果或模型数据文件最小且能直接读的才是可执行入口。用编辑器打开主程序先搜“main”或“program”找到时间递推主循环的位置再搜“parameter”或“const”把网格尺寸、时间步长、震源位置这几组常量抄出来。这一步花十分钟后面能省下几小时的翻车排查。环境方面最常出现的坑是终端直接报“gfortran不是内部或外部命令”“conda不是内部或外部命令”这类信息。它的本质是编译器或Python解释器的路径没加入系统PATH而不是程序本身有问题。Windows下我建议统一装Anaconda并创建一个专门环境装好numpy和scipyFortran代码则用gfortran编译确保编译器和运行时库都是64位。32位和64位混用链接阶段大概率会报“无法定位程序输入点getcurrentpackagefullname”之类的动态库错误这类报错基本都和位数不匹配有关。3.2 最小两层模型一套立刻能用的参数为了验证程序能跑不用上来就上一个真模型我用一个两层介质模型上层2000m/s下层3000m/s横波速度按根号三比例对应。网格200×200网格间距10米震源用20Hz的Ricker子波、垂直集中力放在深度500米处。记录时长1.5秒时间步长0.5毫秒。参数值选取理由网格 nx×nz200×200两层模型只验证物理过程够用即可dxdz10 mS波最短波长约57.8m约5.8点/波长上层 Vp/Vs/ρ2000 / 1155 / 2000 kg/m³Vp/Vs√3接近真实沉积岩比例下层 Vp/Vs/ρ3000 / 1732 / 2200 kg/m³界面反射系数适中便于观察界面深度1000 m给反射波留出清晰的走时窗口震源Ricker20 Hz垂直集中力集中力同时激发P波和S波震源位置x1000 mz500 m离顶面和边界都足够远dt0.5 ms约为二维稳定极限的1/3偏保守记录长度1.5 s反射波有足够时间回到地表这里的关键是网格间距和震源主频的匹配。20Hz主频对应上层横波波长约57.8m10m网格每波长约5.8个点满足伪谱法4到5点的经验要求。如果把主频提到40Hz最短波长降一半网格间距就要缩到5m左右计算量翻四倍这个权衡在第4章还会展开。3.3 主循环跳蛙递推的顺序不能写反拿到程序后主循环通常是这样的结构我把它重写成一个尽量贴近各类初步程序的Python版本# 伪谱法弹性波模拟主循环跳蛙格式二阶时间差分 # 数组形状统一为 (nz, nx)axis0 是深度 zaxis1 是水平 x for it in range(nt): # 第一步由应力更新速度分量 vx dt / rho * ( spectral_derivative(sxx, dx, axis1) # ∂σxx/∂x spectral_derivative(sxz, dz, axis0) # ∂σxz/∂z ) vz dt / rho * ( spectral_derivative(sxz, dx, axis1) # ∂σxz/∂x spectral_derivative(szz, dz, axis0) # ∂σzz/∂z ) # 在震源位置加载垂直集中力源只加在 vz 分量 vz[nsz, nsx] dt / rho[nsz, nsx] * wavelet[it] # 第二步由速度更新应力分量 sxx dt * ( (lam 2.0 * mu) * spectral_derivative(vx, dx, axis1) lam * spectral_derivative(vz, dz, axis0) ) szz dt * ( lam * spectral_derivative(vx, dx, axis1) (lam 2.0 * mu) * spectral_derivative(vz, dz, axis0) ) sxz dt * mu * ( spectral_derivative(vx, dz, axis0) # ∂vx/∂z spectral_derivative(vz, dx, axis1) # ∂vz/∂x ) # 第三步应用吸收边界第4章展开 # 第四步在接收点处把 vx/vz 写入记录道注意这里的存储细节。vx代表水平振动速度vz代表垂直振动速度nsz是深度索引nsx是水平索引。加载垂直集中力时改的是vz而不是vx否则辐射图会绕着一个错误的轴转。如果震源是爆炸源则应该同时往sxx、szz、sxz上加各向同性压力而不是直接改速度分量——很多初步程序把爆炸源实现成“往所有点加同一个速度扰动”得到的结果看着有波但波型比例完全错误。时间递推的顺序是先更新速度再更新应力还是反过来其实可以互换只要震源加在正确的位置、并保持交错时刻的一致性。但每个时间步内部顺序要统一先算完所有速度分量再算所有应力分量不能混着来否则时间同步被打破高频成分会迅速失稳。上面的写法重在清晰效率不是最优。spectral_derivative每调用一次就是一次FFT加一次逆FFT这个循环里一共调用了12次其中对vx的x方向导数和vz的z方向导数在速度更新和应力更新里重复算了。优化时可以先把六个一阶导数场一次性算好再组装应力更新整体能省掉约1/3的FFT开销。初步程序不追求性能但这个逻辑值得记着后续做三维扩展时会用到。3.4 跑通后的第一道验收直达波与反射波的到达时间跑完之后先看接收器输出的两组记录。vz记录上第一个到达的是直达P波初走时约等于震源到接收点的距离除以上层纵波速度随后会看到来自界面的反射P波和反射转换波。如果vz上和vx上除了直达波外什么都没有检查震源类型和界面两侧波阻抗差——速度差太小也会让反射系数低到看不见这时加大两层速度比再试。一个快速的手工验算是把震源到界面的垂直距离和接收点的水平距离代入初等几何关系算出反射P波的走时再与程序输出的记录道对比。以第3.2节的参数为例震源深500m、界面在1000m、接收点水平距离100m时反射P波路径长约1503m按上层Vp2000m/s算走时约0.75秒直达P波走时约0.255秒。误差在1到2毫秒以内说明程序核心逻辑基本正确超过这个量就要回去检查网格方向或介质参数是否装反了。4. 三个必调参数时间步长、吸收边界与震源子波改错了就翻车4.1 时间步长伪谱法的稳定极限不是差分法那个公式伪谱法的空间导数没有频散误差但这不意味着可以无脑用大时间步长。如果时间差分仍然是二阶中心差分稳定性条件来自最大可表示的波数k_maxπ/dx与介质最大波速vmax的乘积。一维情况下理论极限约为0.637·dx/vmax二维时波数向量可以沿对角方向叠加k_max变为π√2/dx极限步长缩到约0.45·dx/vmax三维更严约0.37·dx/vmax。伪谱法能精确表示到Nyquist波数而差分法在高波数部分的振幅响应实际上是衰减的相当于天然滤掉了一部分不稳定成分所以伪谱法对时间步长更敏感。我一般不会顶着极限值用而是取二维极限的一半左右dt 0.3·dx/vmax。这样既留出安全余量又不会因为步长太小让长时程模拟的步数猛增。以第3章那个两层模型为例vmax取下层纵波速3000m/sdx10m二维稳定极限约1.5毫秒取0.5毫秒是极限的1/3属于稳妥选择。如果压缩包代码里时间步长是写死的先按这个公式重新算一遍再跑。判断步长是否过大不一定要等波场爆炸。最快的诊断方法是打印每一时间步的总能量在均匀无吸收模型里总能量应当基本守恒。如果看到某个分量能量随步数单调上升比如从1e-2涨到1e0基本可以断定步长越过稳定极限。把dt缩小到原来的1/4再跑能量曲线趋于平稳就说明问题出在此处而非程序逻辑。提示步长的大小对伪谱法的影响是“全有或全无”的越界一步就会在几十步内爆掉。养成每个新模型先跑50步看能量的习惯比跑完整个记录才发现翻车要省时得多。4.2 吸收边界阻尼带的厚度和衰减系数要一起调初步程序很少带PML最常见的是在计算域四周加一层阻尼带也叫海绵边界或吸收层。它的原理很简单每时间步对边界区的波场乘一个小于1的衰减因子让波在到达人工边界前衰减到可忽略。实现不难但参数配不对时阻尼带本身就会变成反射源效果比不加还糟。阻尼系数一般取成空间位置的函数例如σ(x)σ_max·(x/L)²其中L是阻尼带的网格数x是该点到计算域边界的归一化距离。σ_max的经验范围是2到3倍的vmax/(L·dx)。L的取值至少要覆盖一个中心波长中心波长用震源主频对应的波长来算λ_cvmax/f0。在20Hz主频、3000m/s最大速度的模型里中心波长150米L建议取15到20个网格dx10m时。L太薄时波在阻尼带内还没衰减到位就撞到硬边界反射能量依旧可观。给一段阻尼带实现可以直接替换第3.3节主循环里的“第三步”# 生成二维阻尼衰减系数场四个边界各加 L 个网格 def build_damper(nz, nx, L, vmax, dt): sig_max 3.0 * vmax / (L * dx) # 单位 1/sL*dx 是带的总长度米 damp np.ones((nz, nx), dtypenp.float64) for i in range(L): factor sig_max * ((i 1) / L) ** 2 * dt damp[i, :] * np.exp(-factor) # 上边界 damp[-(i 1), :] * np.exp(-factor) # 下边界 damp[:, i] * np.exp(-factor) # 左边界 damp[:, -(i 1)] * np.exp(-factor) # 右边界 return damp # 每个时间步在递推之后执行 vx * damp vz * damp sxx * damp szz * damp sxz * damp注意角点区域会被重复衰减这个实现在角点的衰减系数比边上大一倍实际影响不大如果要严格处理需要按到最近边界的距离分别计算x和z方向的衰减因子再相乘。更重要的是阻尼带内最好保持常数速度模型不要放界面或强速度梯度否则波在带内产生反射这部分反射同样会污染内部波场。4.3 震源子波Ricker子波的主频和网格间距是配对关系震源子波最常用Ricker表达式是f(t)(1−2π²f₀²(t−t₀)²)exp(−π²f₀²(t−t₀)²)其中t₀一般取1.2到1.5个主频周期让子波初始时刻接近零避免在t0时刻给波场一个阶跃激励。实现如下# Ricker 子波f0 为主频dt 为时间步长 t np.arange(nt) * dt t0 1.2 / f0 wavelet (1.0 - 2.0 * (np.pi * f0 * (t - t0)) ** 2) * \ np.exp(-(np.pi * f0 * (t - t0)) ** 2)主频f₀越高波场分辨率越高能分辨更薄的层但代价是S波最短波长同步变短需要更细的网格。经验约束是每个最短波长至少要有4到5个网格点即dx ≤ v_s_min/(4·f₀)。这里速度取整个模型里最小的S波速度因为S波波长最短最容易频散。以第3章模型为例上层Vs1155m/sf₀20Hz时最短波长约57.8mdx10m相当于每波长约5.8个点处于安全区间。如果把主频从20Hz提到40Hz最短波长降一半dx就必须缩到5m左右计算量涨四倍这就是主频和网格步长的直接权衡。如果压缩包默认震源是爆炸源而你需要同时看P波和S波换成垂直集中力源即可。爆炸源只会辐射纯纵波无论后来怎么调吸收边界和网格横波分量始终是零这一点在验证环节最容易把人带偏。震源加载位置建议离边界至少10个网格否则即使有阻尼带源与人工边界之间的多次反射也会干扰早期波场。5. 伪谱法程序避坑指南5个最常见的翻车现场与排查方法5.1 波场图上一片棋盘格噪声高频Nyquist分量在作怪现象模拟几步后波场图出现颗粒状交替亮暗的棋盘格尤其在震源附近最明显振幅随步数增长。原因单点加载震源在空间上是一个极窄的尖峰它的频谱在Nyquist波数附近仍然有可观的能量。伪谱法对这个分量是全精度放大的不像差分法有天然的抑制于是波场里出现以单个网格为周期的交替扰动视觉上就是棋盘格。解决把震源先做空间平滑再乘子波。常见做法是给震源区一个高斯半径比如σ_source1.5倍的dx让源在空间上分布到8到10个网格点同时检查FFT后是否取了实部虚部残留也会产生类似的高频噪声。如果程序本身没有平滑函数可以在加载震源前对相邻网格按高斯权重分配能量。5.2 边界反射比预期早出现阻尼带没盖住最大波长现象波场图上在计算域边界附近出现强反射弧反射波到达内部接收点的时间明显早于模型里真实界面的理论走时。原因阻尼带厚度L没有按最大中心波长设计。L太薄时长波长成分在带内衰减不够振幅在到达硬边界时仍然可观边界反射自然回传。解决把L加大到至少一个中心波长。用vmax/f0算出中心波长后再换算成网格数如果程序里阻尼带厚度写死改参数或预处理速度模型时把边界区扩展。验证方法是给一个无反射界面的均匀模型跑一次把接收点能量画成时间曲线观察末段是否有明显长时间拖尾的反射能量。阻尼带的σ_max也要同步调到2到3倍vmax/(L·dx)薄带配大衰减、厚带配小衰减两种组合效果不同需要交叉验证。5.3 振幅随时间指数增长直到NaN时间步长越过稳定极限现象前面的波形看着正常到几百步之后某个应力分量量级从1e-2跳到1e20甚至直接变成NaN程序挂掉。原因按照4.1节算出的单方向稳定条件只是一维理论在二维模型里波动能量沿多个方向传播实际允许的步长通常更小。很多初步程序的dt是作者用他的模型试出来的换到你自己的网格尺寸和速度模型后稳定余量可能已经不够。解决把dt缩小到当前值的一半甚至1/4重跑看是否仍然发散。同时建议在时间循环里加一个能量检测每50步打印一次波场总能量看到指数上升就立即终止避免跑完整个记录长度才发现翻车、白烧算力。稳定步长与dx、vmax的具体取值参考4.1的公式但最终以你的模型能量曲线为准这是这类程序最不可省的一步基本功。5.4 横波分量离奇失踪震源类型和参数化把S波灭掉了现象接收记录上只有纵波初至之后全是微弱的低频尾巴理论上应当明显的反射转换波消失vx分量尤其干净。原因两类常见误操作。一是用爆炸源加载它只激发P波S波天然为零二是参数换算时把μ设成了0或很小的值导致S波速度接近0波场根本传播不出去。解决换成垂直集中力源加载在vz分量上同时检查拉梅参数换算μρVs²如果模型文件里Vs列填了0或没填μ就会变成0。一张快速自检图是把Vp、Vs画成按深度的曲线看Vs站点是否与Vp同步变化若Vs全程为0程序里再聪明也算不出S波。5.5 程序在Windows下报动态库或命令找不到环境没有对齐现象终端执行编译命令时报“gfortran不是内部或外部命令”运行Python时报“numpy模块不存在”或者程序启动直接报“无法定位程序输入点getcurrentpackagefullname于动态链接库…”运行就中断。原因三类问题混在一起——编译器或解释器的PATH没有配好、Python环境不对、以及32位/64位运行时库混用。后者在下载了旧版编译好的现成程序包时最容易出现因为动态链接库的位数和主程序不匹配系统加载时就报找不到入口点。解决Fortran源码重新用本地gfortran编译别直接用网上别人编好的exePython部分统一到Anaconda的64位环境建环境后执行conda install numpy scipy别用系统自带的Python。检查位数的方法是打开终端分别敲gfortran --version和python --version确认输出里有没有带32位字样。这一类报错的排查逻辑和网上常见的“conda不是内部或外部命令”完全一样先确认环境变量再确认位数最后才是代码问题。6. 验证伪谱法程序正确性解析解对比与网格收敛性检查写完代码、跑通模拟不等于程序是对的。我验证任何正演程序都走固定的三步解析解走时对比、网格收敛性检查和能量守恒检查。这三步能过滤掉九成以上的隐性错误。第一步用两层介质模型或均匀半空间模型把接收点的波场与解析走时对比。均匀半空间里直达P波走时是r/Vp直达S波走时是r/Vs两层模型里反射P波走时按镜像源法计算公式简单手算即可。把程序输出的单道记录拆成vx和vz两列找到初至时间误差在1到2毫秒内算通过。严格检查可以再加一个垂直自由表面边界对比Rayleigh波存在与否但初步程序一般不需要。第二步是网格收敛性检验。把dx、dz同时减半dt等比缩小重跑同一个模型对比同一接收点的波形。伪谱法如果实现正确两次结果的波形差异应该在1%以内且差值主要集中在高频尾部。如果减半网格后波形明显变化说明原网格本身就不满足分辨率要求需要按第4章的公式重新选择网格间距而不是程序逻辑有问题。第三步是能量监测这个前面提过。在没有阻尼带和震源持续加载的均匀模型中总能量应该守恒在带阻尼带的模型中能量应单调衰减而不是振荡上升。把每步总能量画出来曲线形状正常程序才算真正通过验收。我拿到的每一个伪谱法程序都会先跑这三步再做物理实验。走时对不上先查震源类型能量发散了先查时间步长波形不收敛先查网格间距顺序不要倒过来。这个习惯帮我挡掉了大量“看起来正常其实参数错位”的翻车现场。希望帮到你。本文还有配套的精品资源点击获取