拓十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

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

伪谱法弹性波正演模拟:从原理到避坑实战指南

伪谱法弹性波正演模拟:从原理到避坑实战指南

简介:这是一套面向地球物理、工程波动模拟初学者的初步虚谱法(伪谱法)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、σxz:

rho ∂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”是模型数组还是标量,以及有没有做速度到拉梅参数的换算,很多结果怪异的问题都出在这一步。

时间递推用跳蛙格式,即速度在n+1/2时刻、应力在n时刻交错更新。它是二阶精度的,空间误差由伪谱法控制在几乎为零,时间误差就成了总误差的主要来源。如果要做长时间模拟,可以换四阶Runge-Kutta,但每步要算四次导数场,成本高很多,初步程序保持二阶中心差分即可。

2.3 波数域求导算子:整个伪谱法程序的核心就这一段

把空间导数封装成一个函数,后续所有递推都复用它。Python实现如下:

import numpy as np def spectral_derivative(field, dx, axis=0): """ 沿指定轴对场做波数域一阶求导。 以二维波场形状 (nz, nx) 为准: axis=0 对应 z 方向,间距为 dz;axis=1 对应 x 方向,间距为 dx。 """ nx = field.shape[axis] # 角波数向量:fftfreq 返回频率索引,乘 2*pi 后是角波数,单位 rad/m k = 2.0 * np.pi * np.fft.fftfreq(nx, d=dx) # 把波数向量广播到 field 的目标轴 shape = [1] * field.ndim shape[axis] = nx k = k.reshape(shape) # 正变换、在波数域乘 i*k、反变换取实部 derivative = np.fft.ifft( np.fft.fft(field, axis=axis) * (1j * k), axis=axis ).real return derivative

这段的要点有三个。第一,fftfreq(nx, d=dx)返回的频率索引从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和scipy;Fortran代码则用gfortran编译,确保编译器和运行时库都是64位。32位和64位混用,链接阶段大概率会报“无法定位程序输入点getcurrentpackagefullname”之类的动态库错误,这类报错基本都和位数不匹配有关。

3.2 最小两层模型:一套立刻能用的参数

为了验证程序能跑,不用上来就上一个真模型,我用一个两层介质模型:上层2000m/s,下层3000m/s,横波速度按根号三比例对应。网格200×200,网格间距10米,震源用20Hz的Ricker子波、垂直集中力,放在深度500米处。记录时长1.5秒,时间步长0.5毫秒。

参数值选取理由
网格 nx×nz200×200两层模型只验证物理过程,够用即可
dx=dz10 mS波最短波长约57.8m,约5.8点/波长
上层 Vp/Vs/ρ2000 / 1155 / 2000 kg/m³Vp/Vs=√3,接近真实沉积岩比例
下层 Vp/Vs/ρ3000 / 1732 / 2200 kg/m³界面反射系数适中,便于观察
界面深度1000 m给反射波留出清晰的走时窗口
震源Ricker,20 Hz,垂直集中力集中力同时激发P波和S波
震源位置x=1000 m,z=500 m离顶面和边界都足够远
dt0.5 ms约为二维稳定极限的1/3,偏保守
记录长度1.5 s反射波有足够时间回到地表

这里的关键是网格间距和震源主频的匹配。20Hz主频对应上层横波波长约57.8m,10m网格每波长约5.8个点,满足伪谱法4到5点的经验要求。如果把主频提到40Hz,最短波长降一半,网格间距就要缩到5m左右,计算量翻四倍,这个权衡在第4章还会展开。

3.3 主循环:跳蛙递推的顺序不能写反

拿到程序后,主循环通常是这样的结构,我把它重写成一个尽量贴近各类初步程序的Python版本:

# 伪谱法弹性波模拟主循环(跳蛙格式,二阶时间差分) # 数组形状统一为 (nz, nx),axis=0 是深度 z,axis=1 是水平 x for it in range(nt): # 第一步:由应力更新速度分量 vx += dt / rho * ( spectral_derivative(sxx, dx, axis=1) + # ∂σxx/∂x spectral_derivative(sxz, dz, axis=0) # ∂σxz/∂z ) vz += dt / rho * ( spectral_derivative(sxz, dx, axis=1) + # ∂σxz/∂x spectral_derivative(szz, dz, axis=0) # ∂σzz/∂z ) # 在震源位置加载垂直集中力源,只加在 vz 分量 vz[nsz, nsx] += dt / rho[nsz, nsx] * wavelet[it] # 第二步:由速度更新应力分量 sxx += dt * ( (lam + 2.0 * mu) * spectral_derivative(vx, dx, axis=1) + lam * spectral_derivative(vz, dz, axis=0) ) szz += dt * ( lam * spectral_derivative(vx, dx, axis=1) + (lam + 2.0 * mu) * spectral_derivative(vz, dz, axis=0) ) sxz += dt * mu * ( spectral_derivative(vx, dz, axis=0) + # ∂vx/∂z spectral_derivative(vz, dx, axis=1) # ∂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,按上层Vp=2000m/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/s,dx=10m,二维稳定极限约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的取值至少要覆盖一个中心波长,中心波长用震源主频对应的波长来算:λ_c=vmax/f0。在20Hz主频、3000m/s最大速度的模型里,中心波长150米,L建议取15到20个网格(dx=10m时)。L太薄时,波在阻尼带内还没衰减到位就撞到硬边界,反射能量依旧可观。

给一段阻尼带实现,可以直接替换第3.3节主循环里的“第三步”:

# 生成二维阻尼衰减系数场,四个边界各加 L 个网格 def build_damper(nz, nx, L, vmax, dt): sig_max = 3.0 * vmax / (L * dx) # 单位 1/s,L*dx 是带的总长度(米) damp = np.ones((nz, nx), dtype=np.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个主频周期,让子波初始时刻接近零,避免在t=0时刻给波场一个阶跃激励。实现如下:

# 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章模型为例,上层Vs=1155m/s,f₀=20Hz时最短波长约57.8m,dx=10m相当于每波长约5.8个点,处于安全区间。如果把主频从20Hz提到40Hz,最短波长降一半,dx就必须缩到5m左右,计算量涨四倍,这就是主频和网格步长的直接权衡。

如果压缩包默认震源是爆炸源,而你需要同时看P波和S波,换成垂直集中力源即可。爆炸源只会辐射纯纵波,无论后来怎么调吸收边界和网格,横波分量始终是零,这一点在验证环节最容易把人带偏。震源加载位置建议离边界至少10个网格,否则即使有阻尼带,源与人工边界之间的多次反射也会干扰早期波场。

5. 伪谱法程序避坑指南:5个最常见的翻车现场与排查方法

5.1 波场图上一片棋盘格噪声:高频Nyquist分量在作怪

现象:模拟几步后,波场图出现颗粒状交替亮暗的棋盘格,尤其在震源附近最明显,振幅随步数增长。

原因:单点加载震源在空间上是一个极窄的尖峰,它的频谱在Nyquist波数附近仍然有可观的能量。伪谱法对这个分量是全精度放大的,不像差分法有天然的抑制,于是波场里出现以单个网格为周期的交替扰动,视觉上就是棋盘格。

解决:把震源先做空间平滑,再乘子波。常见做法是给震源区一个高斯半径,比如σ_source=1.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编译,别直接用网上别人编好的exe;Python部分统一到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章的公式重新选择网格间距,而不是程序逻辑有问题。

第三步是能量监测,这个前面提过。在没有阻尼带和震源持续加载的均匀模型中,总能量应该守恒;在带阻尼带的模型中,能量应单调衰减而不是振荡上升。把每步总能量画出来,曲线形状正常,程序才算真正通过验收。

我拿到的每一个伪谱法程序,都会先跑这三步再做物理实验。走时对不上先查震源类型,能量发散了先查时间步长,波形不收敛先查网格间距,顺序不要倒过来。这个习惯帮我挡掉了大量“看起来正常其实参数错位”的翻车现场。希望帮到你。

本文还有配套的精品资源,点击获取

返回列表