简介fd.zip是一份面向地球物理勘探与数值模拟初学者的四阶声波有限差分正演程序包。压缩包共三个文件包含一个C语言源码和两个数据文件整体仅51KB。源码以四阶有限差分方法求解声波波动方程覆盖网格离散化、时间步进、边界条件与初始条件设置并涉及雷克子波加载、稳定性条件判断等实现细节方便对照代码理解正演模拟中波动方程、有限差分格式与时间递推的核心流程两个数据文件分别保存模拟记录和波前信息可用于绘制波场快照、分析声波在不同介质中的传播路径、反射与折射特征。已有125人学习这一资源适合正在入门地震正演模拟或需要参考紧凑C语言实现的开发者研读通过阅读和运行该程序可快速搭建自己的声波正演实验。1. 四阶声波有限差分地震正演在解决什么在地下速度模型相对简单的场景里地震正演通常先拿声波方程做骨架只考虑纵波标量压力场省去横波和转换波。四阶有限差分是这个骨架里最常用的空间离散方案——五点中心差分把空间导数精度从传统的 O(h²) 提到 O(h⁴)让网格尺度可以比二阶格式放得更粗又不像交错网格那样需要定义额外的速度/应力分量。fd.zip 这类基本程序解决的核心问题是给一个速度模型和一个 Ricker 子波如何用可复现的四阶 FD 算法生成一套合成地震记录。它适合打算自己搭正演基线、校对商业软件模板、或者为后续弹性波正演先跑通流程的工程师。2. 四阶精度的声波离散格式与稳定性条件2.1 二维声波方程与五点差分模板怎样做到四阶起点是标量声波方程∂²P/∂t² c² ( ∂²P/∂x² ∂²P/∂z² ) S其中 P 是声压c 是介质速度S 是震源项。基本程序里把 P 作为唯一求解量不引入粒子速度分量这样内存占用和代码逻辑都最直接。对空间二阶导采用五点中心差分替换常规三点差分∂²P/∂x² ≈ [ -P(x-2h) 16P(x-h) - 30P(x) 16P(xh) - P(x2h) ] / (12h²)截断误差为 O(h⁴)。把 P(x±h) 与 P(x±2h) 在 x 处同时做 Taylor 展开要求 h⁻²、h⁻¹、h⁰、h¹、h² 各阶项分别抵消再让展开系数满足二阶导的单一条件就会解出这组系数。实际写程序时建议把系数记成下面这个表格手写循环时不容易敲错。加权系数u(x-2h)u(x-h)u(x)u(xh)u(x2h)截断误差二阶中心差分01-210O(h²)四阶中心差分-1/124/3-5/24/3-1/12O(h⁴)用 -1/12、4/3、-5/2、4/3、-1/12 去乘相邻数组值再统一除以 12就是在代码里常见的(-u[i-2] 16*u[i-1] - 30*u[i] 16*u[i1] - u[i2]) / 12.0。这个写法只关心 x 方向z 方向完全同理。提示这里的“四阶”只指空间导数精度不是时间积分精度。四阶龙格库塔是另一套时间推进算法不能与空间四阶差分混用对基本正演程序时间方向保留二阶即可实现简单且稳定性条件明确。2.2 时间二阶推进与 CFL 稳定参数完整时间推进采用P^(n1)(i,j) 2P^n(i,j) - P^(n-1)(i,j) (c·Δt)² · [Lxx Lzz]P^n这就是教科书里常写的 O(Δt², Δx⁴) 格式时间中心差分二阶空间五点差分四阶。它只需保存 p0、p1、p2 三个时间层内存开销对二维模型非常友好网格达到百万量级时三层 double 数组也只要 24 MB 左右。稳定性条件由 von Neumann 分析导出。五点模板在二维网格上的最大特征值出现在最高可分辨波数处由此得到(c·Δt/h)² ≤ 3/8即 c·Δt/h ≤ 0.612工程上通常直接取 0.6写成Δt 0.6·h / c_maxc_max 是整个模型的最快速度不是震源附近速度。如果模型里有 3500 m/s 的高速层而震源布在 1800 m/s 区域按 1800 算出的 Δt 会在高速层快速溢出正演跑几十步就变成 NaN。反过来按高速层取 Δt 只是多花一点计算量不会错。下面这段小型函数可以独立验证五点模板的写法/* 一维五点二阶导要求调用时 i 至少大于等于 2 且小于 n-2 */ static double d2x(const double *u, int i) { return (-u[i-2] 16.0*u[i-1] - 30.0*u[i] 16.0*u[i1] - u[i2]) / 12.0; }调用前要注意数组边界i 从 2 到 n-3否则会越界读取。实际二维程序里x 和 z 两个方向各算一次再把结果相加这也是第 5 章完整代码里 Laplacian 的计算方式。3. 边界吸收与网格裁剪Cerjan 阻尼带的核心设置3.1 为什么零边界和简单衰减会毁掉合成记录正演网格必然是有限范围。把四周边界直接当成零位移会在边界和四角持续产生强反射因为四阶模板在靠近边界处无法照常使用边界反射会带着明显的锯齿形态向内部传播。地震正演中这类伪反射能量往往比真实弱反射层信号还强直接掩盖有效同相轴。基本程序里最常见的处理是 Cerjan 阻尼带也叫海绵吸收边界。原理很简单在计算区域四周包裹一圈衰减带每个时间步推进后让波场乘上一个随距离边界衰减的系数波传进吸收带后能量逐步消失到达硬边界时已经小到可以忽略。3.2 Cerjan 阻尼带公式与参数组合衰减系数定义为到吸收带内壁的距离 d 的函数g(d) exp( -α² · d² )内部区域 d 0g 1波场不受影响越靠近最外层边界d 越大衰减越强。实现时把 x、z 两个方向分别做一维系数再用 gx[i] · gz[j] 相乘。四角区域同时被两个方向的系数压制反射抑制效果比单方向叠加好得多。参数经验取值说明NB吸收带厚度2040 格至少要覆盖几个网格波长α衰减强度0.050.12 / 格过大会在带内形成硬边界过小吸收不足最外层衰减量exp( -(α·NB)² )建议小于 10⁻²稳定时取 10⁻³这两个参数是耦合的。NB 取 20 时α 取 0.08最外层衰减约为 exp(-2.56) ≈ 0.077余量偏大NB 取 30、α 取 0.08 时最外层约为 0.003能压住大部分反射。更稳妥的做法是先设 NB30跑一个均匀模型观察边界反射残差再在 0.050.12 之间调 α。以下是可直接搬用的系数生成代码static void make_taper(double *g, int n, int nb, double alpha) { for (int i 0; i n; i) { double d 0.0; if (i nb) d nb - i; else if (i n - nb) d i - (n - nb - 1); g[i] exp(-alpha * alpha * d * d); } }调用方式make_taper(gx, NX, NB, 0.08); make_taper(gz, NZ, NB, 0.08);。每个时间步推进完成后对全部网格点执行p2[j*NXi] * gx[i] * gz[j]。由于带内 d0内部点乘数为 1不需要为内部循环特判边界代码分支更少。注意如果吸收带边缘出现二次波动往往不是 NB 太薄而是 α 太大。α 过大时 g(d) 在带内下降过快形同硬边界会反向激发新波。反之 α 太小直达波穿过吸收带后仍在边界反弹表现为从边界方向斜插进来的弧形伪同相轴在单道记录上很难和真实反射区分。如果后续要提升模拟保真度可以把 Cerjan 带替换成 PML。PML 对广角入射更透明却能大幅增加内存和实现复杂度需要额外维护辅助波场变量。对 fd.zip 这类基本程序Cerjan 阻尼带是性价比最高的缺省方案。4. 震源子波、空间步长与时间步长的匹配计算4.1 Ricker 子波的时域表达式与频段含义地震正演里震源时间函数基本都用 Ricker 子波w(t) [ 1 - 2π² f₀² (t - t0)² ] · exp( -π² f₀² (t - t0)² )f₀ 是峰值频率t0 是主峰时刻。Ricker 子波零相位、低频干净是合成记录里最容易识别的波形。有效频率范围可粗略按 2.5 倍 f₀ 估算空间采样必须按这个高频端来约束否则高频分量会先出现数值色散。4.2 空间步长、时间步长与每波长网格点数关系四阶格式对每波长采样点数 Nppw 的经验要求大约在 68 个点。按此设计网格步长h ≤ c_min / ( 2.5 · f₀ · Nppw )这里用低速介质的最小速度 c_min是为了保证最短波长有足够网格分辨率。之后再用全局最大速度 c_max 约束时间步Δt ≤ 0.6 · h / c_max顺序不能反。如果先用 c_max 定 h低速区波长变短很可能在浅部低速层产生明显色散。举一个可复算的例子介质情况速度选取参数低速目标层1800 m/sf₀25 HzNppw8h ≤ 1800/(2.5×25×8) 3.6 m取 4 m高速基底3500 m/sc_max3500由 h4 得 Δt ≤ 0.6×4/3500 ≈ 0.686 ms实际取 Δt0.6 ms留出稳定余量。这个组合下Ricker 子波在 1800 m/s 低速层的每波长采样约为 1800/(25×4)18 点在高速层约为 35 点色散可以压得比较干净。子波时长建议覆盖 35 个主周期t0 约取 3/f₀。t0 太短子波起始段被截断会产生高频毛刺t0 太长前段空转时间浪费。可以用下面代码先输出子波序列检查形态for (int n 0; n 200; n) { double t n * dt; double w ricker(t, f0, t0); printf(%10.6f %13.6e\n, t, w); }观察输出波形主峰两侧应当各有一个对称的旁瓣且首尾趋于零。如果首尾在零值附近有明显跳变就需要加大 t0 或增加时间采样长度。4.3 一个常见误区很多人以为四阶精度足够高可以把 Nppw 压到 5 以下。实际在速度模型存在尖锐界面时四阶模板跨越界面会引入非物理振荡这个误差并不会因为空间精度升高而消失。基本程序尽量让速度场平滑过渡尖锐界面留给专门的界面条件或可变网格方法处理。5. 基本程序落地二维均匀介质四阶 FD 的 C 实现5.1 数组分工与时间层滚动完整程序围绕三组数组展开数组作用说明p0、p1、p2前时刻、当前时刻、下一时刻波场三指针循环滚动c2速度平方每格一个 double可预计算gx、gz吸收带衰减系数每个时间步乘性应用时间层更新不复制数据只交换指针推进时从 p1 读写入 p2结束后把 p0 指向旧 p1p1 指向旧 p2p2 指向旧 p0。这样下一时间步仍从 p1 读取当前波场。5.2 带吸收边界的二维四阶 FD 全部代码下面代码把前面所有参数设置集中在一个可编译的 C 程序里。模型为 400×400 网格均匀介质 2500 m/s源在网格中心接收器在水平方向 3/4 位置每 25 步打印一个采样点。#include stdio.h #include math.h #define NX 400 #define NZ 400 #define NT 1200 #define NB 30 static double gx[NX], gz[NZ]; static double ricker(double t, double f0, double t0) { double u M_PI * f0 * (t - t0); u u * u; return (1.0 - 2.0*u) * exp(-u); } static void make_taper(double *g, int n, int nb, double alpha) { for (int i 0; i n; i) { double d 0.0; if (i nb) d nb - i; else if (i n - nb) d i - (n - nb - 1); g[i] exp(-alpha * alpha * d * d); } } int main(void) { static double p0[NX*NZ], p1[NX*NZ], p2[NX*NZ]; static double c2[NX*NZ]; double c0 2500.0, dx 4.0, dt 6.0e-4; double f0 25.0, t0 0.09; double sr dt * dt / (dx * dx); double alpha 0.08; int isrc NX/2, jsrc NZ/2; int iRec 3*NX/4, jRec NZ/2; for (int j 0; j NZ; j) for (int i 0; i NX; i) c2[j*NX i] c0 * c0; make_taper(gx, NX, NB, alpha); make_taper(gz, NZ, NB, alpha); for (int n 0; n NT; n) { double t n * dt; /* 内部区域空间四阶、时间二阶推进 */ for (int j 2; j NZ-2; j) for (int i 2; i NX-2; i) { int ix j*NX i; double lap_x (-p1[ix-2] 16.0*p1[ix-1] - 30.0*p1[ix] 16.0*p1[ix1] - p1[ix2]) / 12.0; double lap_z (-p1[ix-2*NX] 16.0*p1[ix-NX] - 30.0*p1[ix] 16.0*p1[ixNX] - p1[ix2*NX]) / 12.0; p2[ix] 2.0*p1[ix] - p0[ix] c2[ix]*sr*(lap_x lap_z); } /* 震源项按 (c*dt/dx)^2 * w(t) 注入 */ p2[jsrc*NX isrc] c2[jsrc*NX isrc]*sr * ricker(t, f0, t0); /* 整场乘吸收衰减带内系数为 1.0 */ for (int j 0; j NZ; j) for (int i 0; i NX; i) p2[j*NX i] * gx[i]*gz[j]; if (n % 25 0) printf(%d %.10f\n, n, p1[jRec*NX iRec]); /* 滚动时间层不拷贝整块波场 */ double *tmp p0; p0 p1; p1 p2; p2 tmp; } return 0; }这里sr就是 Δt²/h²。空间导数算完后乘c2[ix]*sr等价于把波动方程两边同时乘 Δt² 后的显式更新。震源加载放在推进之后表示当前时间步汇入压力增量。放在推进之前会把子波反向并产生一个时间步的相位偏差。如果发现波形形状完全正确但到时差一个 Δt首先检查震源注入位置。数组全部声明为static避免大块栈内存导致程序崩溃。网格进一步扩大后应改成calloc动态分配同时把 NX、NZ 作为参数传入推进函数。5.3 编译与运行编译命令gcc -O2 -o fdfwd fdfwd.c -lm ./fdfwd输出两列时间步序号和接收点波场值。首波到时与接收距离、介质速度的对应关系是立即可以检查的指标源到接收器水平距离 100 格×4 m400 mc2500 m/s理论到时约 0.16 s对应第 267 步。输出会在第 275 步前后出现第一个明显极值。6. 震源到时、边界伪影与数值色散的排查6.1 均匀模型到时验证上面代码的接收器位置已经给出一个天然验证点。理论到时 t r/c 400/2500 0.16 s对应 n267。程序每 25 步输出一次第一个极值应当在 n275 附近出现。如果要做更精确的到时验证把输出改成每步打印再用一个小脚本完成自动比对./fdfwd | awk NR1{first$1} $20{next} {print $1, $2} | head -20首波到达前波场应接近零到达后开始出现完整的三峰波形。首峰位置偏差超过 2Δt 时优先怀疑震源加载顺序或时间层滚动顺序。6.2 三个高频问题快速排查症状最可能原因检查步骤波前出现梳齿状高频振荡空间步长相对有效频率过粗缩小 h 后看波形是否明显变干净记录尾部持续有半正弦隆起吸收带 α 太小或 NB 太薄打印 gx 在最外层的值确认小于 10⁻³到时偏差大于 3Δt 且不收敛震源项顺序或指针滚动有误手推前三个时间步与程序逐值对照6.3 用多次数值导验证网格色散一个容易操作的分辨率验证技巧把接收道写成文件对时间序列做一次二阶导再做一次四阶导和原始波形对比主峰位置。在色散可控的范围内导数波形主峰应与原始波形对齐如果导数波形整体拖尾、旁瓣被拉开说明网格对高频成分的相速度已经明显偏离真实值。这个技巧适合在放大网格时快速找到允许的最大 h。对基本正演程序来说先跑通均匀介质线再逐步加入层状模型始终用理论到时作为基线之后的真实构造差异才可信而不是数值格式变化带来的假象。本文还有配套的精品资源点击获取 四阶有限差分声波方程正演：从空间离散到合成地震记录的实现要点
简介fd.zip是一份面向地球物理勘探与数值模拟初学者的四阶声波有限差分正演程序包。压缩包共三个文件包含一个C语言源码和两个数据文件整体仅51KB。源码以四阶有限差分方法求解声波波动方程覆盖网格离散化、时间步进、边界条件与初始条件设置并涉及雷克子波加载、稳定性条件判断等实现细节方便对照代码理解正演模拟中波动方程、有限差分格式与时间递推的核心流程两个数据文件分别保存模拟记录和波前信息可用于绘制波场快照、分析声波在不同介质中的传播路径、反射与折射特征。已有125人学习这一资源适合正在入门地震正演模拟或需要参考紧凑C语言实现的开发者研读通过阅读和运行该程序可快速搭建自己的声波正演实验。1. 四阶声波有限差分地震正演在解决什么在地下速度模型相对简单的场景里地震正演通常先拿声波方程做骨架只考虑纵波标量压力场省去横波和转换波。四阶有限差分是这个骨架里最常用的空间离散方案——五点中心差分把空间导数精度从传统的 O(h²) 提到 O(h⁴)让网格尺度可以比二阶格式放得更粗又不像交错网格那样需要定义额外的速度/应力分量。fd.zip 这类基本程序解决的核心问题是给一个速度模型和一个 Ricker 子波如何用可复现的四阶 FD 算法生成一套合成地震记录。它适合打算自己搭正演基线、校对商业软件模板、或者为后续弹性波正演先跑通流程的工程师。2. 四阶精度的声波离散格式与稳定性条件2.1 二维声波方程与五点差分模板怎样做到四阶起点是标量声波方程∂²P/∂t² c² ( ∂²P/∂x² ∂²P/∂z² ) S其中 P 是声压c 是介质速度S 是震源项。基本程序里把 P 作为唯一求解量不引入粒子速度分量这样内存占用和代码逻辑都最直接。对空间二阶导采用五点中心差分替换常规三点差分∂²P/∂x² ≈ [ -P(x-2h) 16P(x-h) - 30P(x) 16P(xh) - P(x2h) ] / (12h²)截断误差为 O(h⁴)。把 P(x±h) 与 P(x±2h) 在 x 处同时做 Taylor 展开要求 h⁻²、h⁻¹、h⁰、h¹、h² 各阶项分别抵消再让展开系数满足二阶导的单一条件就会解出这组系数。实际写程序时建议把系数记成下面这个表格手写循环时不容易敲错。加权系数u(x-2h)u(x-h)u(x)u(xh)u(x2h)截断误差二阶中心差分01-210O(h²)四阶中心差分-1/124/3-5/24/3-1/12O(h⁴)用 -1/12、4/3、-5/2、4/3、-1/12 去乘相邻数组值再统一除以 12就是在代码里常见的(-u[i-2] 16*u[i-1] - 30*u[i] 16*u[i1] - u[i2]) / 12.0。这个写法只关心 x 方向z 方向完全同理。提示这里的“四阶”只指空间导数精度不是时间积分精度。四阶龙格库塔是另一套时间推进算法不能与空间四阶差分混用对基本正演程序时间方向保留二阶即可实现简单且稳定性条件明确。2.2 时间二阶推进与 CFL 稳定参数完整时间推进采用P^(n1)(i,j) 2P^n(i,j) - P^(n-1)(i,j) (c·Δt)² · [Lxx Lzz]P^n这就是教科书里常写的 O(Δt², Δx⁴) 格式时间中心差分二阶空间五点差分四阶。它只需保存 p0、p1、p2 三个时间层内存开销对二维模型非常友好网格达到百万量级时三层 double 数组也只要 24 MB 左右。稳定性条件由 von Neumann 分析导出。五点模板在二维网格上的最大特征值出现在最高可分辨波数处由此得到(c·Δt/h)² ≤ 3/8即 c·Δt/h ≤ 0.612工程上通常直接取 0.6写成Δt 0.6·h / c_maxc_max 是整个模型的最快速度不是震源附近速度。如果模型里有 3500 m/s 的高速层而震源布在 1800 m/s 区域按 1800 算出的 Δt 会在高速层快速溢出正演跑几十步就变成 NaN。反过来按高速层取 Δt 只是多花一点计算量不会错。下面这段小型函数可以独立验证五点模板的写法/* 一维五点二阶导要求调用时 i 至少大于等于 2 且小于 n-2 */ static double d2x(const double *u, int i) { return (-u[i-2] 16.0*u[i-1] - 30.0*u[i] 16.0*u[i1] - u[i2]) / 12.0; }调用前要注意数组边界i 从 2 到 n-3否则会越界读取。实际二维程序里x 和 z 两个方向各算一次再把结果相加这也是第 5 章完整代码里 Laplacian 的计算方式。3. 边界吸收与网格裁剪Cerjan 阻尼带的核心设置3.1 为什么零边界和简单衰减会毁掉合成记录正演网格必然是有限范围。把四周边界直接当成零位移会在边界和四角持续产生强反射因为四阶模板在靠近边界处无法照常使用边界反射会带着明显的锯齿形态向内部传播。地震正演中这类伪反射能量往往比真实弱反射层信号还强直接掩盖有效同相轴。基本程序里最常见的处理是 Cerjan 阻尼带也叫海绵吸收边界。原理很简单在计算区域四周包裹一圈衰减带每个时间步推进后让波场乘上一个随距离边界衰减的系数波传进吸收带后能量逐步消失到达硬边界时已经小到可以忽略。3.2 Cerjan 阻尼带公式与参数组合衰减系数定义为到吸收带内壁的距离 d 的函数g(d) exp( -α² · d² )内部区域 d 0g 1波场不受影响越靠近最外层边界d 越大衰减越强。实现时把 x、z 两个方向分别做一维系数再用 gx[i] · gz[j] 相乘。四角区域同时被两个方向的系数压制反射抑制效果比单方向叠加好得多。参数经验取值说明NB吸收带厚度2040 格至少要覆盖几个网格波长α衰减强度0.050.12 / 格过大会在带内形成硬边界过小吸收不足最外层衰减量exp( -(α·NB)² )建议小于 10⁻²稳定时取 10⁻³这两个参数是耦合的。NB 取 20 时α 取 0.08最外层衰减约为 exp(-2.56) ≈ 0.077余量偏大NB 取 30、α 取 0.08 时最外层约为 0.003能压住大部分反射。更稳妥的做法是先设 NB30跑一个均匀模型观察边界反射残差再在 0.050.12 之间调 α。以下是可直接搬用的系数生成代码static void make_taper(double *g, int n, int nb, double alpha) { for (int i 0; i n; i) { double d 0.0; if (i nb) d nb - i; else if (i n - nb) d i - (n - nb - 1); g[i] exp(-alpha * alpha * d * d); } }调用方式make_taper(gx, NX, NB, 0.08); make_taper(gz, NZ, NB, 0.08);。每个时间步推进完成后对全部网格点执行p2[j*NXi] * gx[i] * gz[j]。由于带内 d0内部点乘数为 1不需要为内部循环特判边界代码分支更少。注意如果吸收带边缘出现二次波动往往不是 NB 太薄而是 α 太大。α 过大时 g(d) 在带内下降过快形同硬边界会反向激发新波。反之 α 太小直达波穿过吸收带后仍在边界反弹表现为从边界方向斜插进来的弧形伪同相轴在单道记录上很难和真实反射区分。如果后续要提升模拟保真度可以把 Cerjan 带替换成 PML。PML 对广角入射更透明却能大幅增加内存和实现复杂度需要额外维护辅助波场变量。对 fd.zip 这类基本程序Cerjan 阻尼带是性价比最高的缺省方案。4. 震源子波、空间步长与时间步长的匹配计算4.1 Ricker 子波的时域表达式与频段含义地震正演里震源时间函数基本都用 Ricker 子波w(t) [ 1 - 2π² f₀² (t - t0)² ] · exp( -π² f₀² (t - t0)² )f₀ 是峰值频率t0 是主峰时刻。Ricker 子波零相位、低频干净是合成记录里最容易识别的波形。有效频率范围可粗略按 2.5 倍 f₀ 估算空间采样必须按这个高频端来约束否则高频分量会先出现数值色散。4.2 空间步长、时间步长与每波长网格点数关系四阶格式对每波长采样点数 Nppw 的经验要求大约在 68 个点。按此设计网格步长h ≤ c_min / ( 2.5 · f₀ · Nppw )这里用低速介质的最小速度 c_min是为了保证最短波长有足够网格分辨率。之后再用全局最大速度 c_max 约束时间步Δt ≤ 0.6 · h / c_max顺序不能反。如果先用 c_max 定 h低速区波长变短很可能在浅部低速层产生明显色散。举一个可复算的例子介质情况速度选取参数低速目标层1800 m/sf₀25 HzNppw8h ≤ 1800/(2.5×25×8) 3.6 m取 4 m高速基底3500 m/sc_max3500由 h4 得 Δt ≤ 0.6×4/3500 ≈ 0.686 ms实际取 Δt0.6 ms留出稳定余量。这个组合下Ricker 子波在 1800 m/s 低速层的每波长采样约为 1800/(25×4)18 点在高速层约为 35 点色散可以压得比较干净。子波时长建议覆盖 35 个主周期t0 约取 3/f₀。t0 太短子波起始段被截断会产生高频毛刺t0 太长前段空转时间浪费。可以用下面代码先输出子波序列检查形态for (int n 0; n 200; n) { double t n * dt; double w ricker(t, f0, t0); printf(%10.6f %13.6e\n, t, w); }观察输出波形主峰两侧应当各有一个对称的旁瓣且首尾趋于零。如果首尾在零值附近有明显跳变就需要加大 t0 或增加时间采样长度。4.3 一个常见误区很多人以为四阶精度足够高可以把 Nppw 压到 5 以下。实际在速度模型存在尖锐界面时四阶模板跨越界面会引入非物理振荡这个误差并不会因为空间精度升高而消失。基本程序尽量让速度场平滑过渡尖锐界面留给专门的界面条件或可变网格方法处理。5. 基本程序落地二维均匀介质四阶 FD 的 C 实现5.1 数组分工与时间层滚动完整程序围绕三组数组展开数组作用说明p0、p1、p2前时刻、当前时刻、下一时刻波场三指针循环滚动c2速度平方每格一个 double可预计算gx、gz吸收带衰减系数每个时间步乘性应用时间层更新不复制数据只交换指针推进时从 p1 读写入 p2结束后把 p0 指向旧 p1p1 指向旧 p2p2 指向旧 p0。这样下一时间步仍从 p1 读取当前波场。5.2 带吸收边界的二维四阶 FD 全部代码下面代码把前面所有参数设置集中在一个可编译的 C 程序里。模型为 400×400 网格均匀介质 2500 m/s源在网格中心接收器在水平方向 3/4 位置每 25 步打印一个采样点。#include stdio.h #include math.h #define NX 400 #define NZ 400 #define NT 1200 #define NB 30 static double gx[NX], gz[NZ]; static double ricker(double t, double f0, double t0) { double u M_PI * f0 * (t - t0); u u * u; return (1.0 - 2.0*u) * exp(-u); } static void make_taper(double *g, int n, int nb, double alpha) { for (int i 0; i n; i) { double d 0.0; if (i nb) d nb - i; else if (i n - nb) d i - (n - nb - 1); g[i] exp(-alpha * alpha * d * d); } } int main(void) { static double p0[NX*NZ], p1[NX*NZ], p2[NX*NZ]; static double c2[NX*NZ]; double c0 2500.0, dx 4.0, dt 6.0e-4; double f0 25.0, t0 0.09; double sr dt * dt / (dx * dx); double alpha 0.08; int isrc NX/2, jsrc NZ/2; int iRec 3*NX/4, jRec NZ/2; for (int j 0; j NZ; j) for (int i 0; i NX; i) c2[j*NX i] c0 * c0; make_taper(gx, NX, NB, alpha); make_taper(gz, NZ, NB, alpha); for (int n 0; n NT; n) { double t n * dt; /* 内部区域空间四阶、时间二阶推进 */ for (int j 2; j NZ-2; j) for (int i 2; i NX-2; i) { int ix j*NX i; double lap_x (-p1[ix-2] 16.0*p1[ix-1] - 30.0*p1[ix] 16.0*p1[ix1] - p1[ix2]) / 12.0; double lap_z (-p1[ix-2*NX] 16.0*p1[ix-NX] - 30.0*p1[ix] 16.0*p1[ixNX] - p1[ix2*NX]) / 12.0; p2[ix] 2.0*p1[ix] - p0[ix] c2[ix]*sr*(lap_x lap_z); } /* 震源项按 (c*dt/dx)^2 * w(t) 注入 */ p2[jsrc*NX isrc] c2[jsrc*NX isrc]*sr * ricker(t, f0, t0); /* 整场乘吸收衰减带内系数为 1.0 */ for (int j 0; j NZ; j) for (int i 0; i NX; i) p2[j*NX i] * gx[i]*gz[j]; if (n % 25 0) printf(%d %.10f\n, n, p1[jRec*NX iRec]); /* 滚动时间层不拷贝整块波场 */ double *tmp p0; p0 p1; p1 p2; p2 tmp; } return 0; }这里sr就是 Δt²/h²。空间导数算完后乘c2[ix]*sr等价于把波动方程两边同时乘 Δt² 后的显式更新。震源加载放在推进之后表示当前时间步汇入压力增量。放在推进之前会把子波反向并产生一个时间步的相位偏差。如果发现波形形状完全正确但到时差一个 Δt首先检查震源注入位置。数组全部声明为static避免大块栈内存导致程序崩溃。网格进一步扩大后应改成calloc动态分配同时把 NX、NZ 作为参数传入推进函数。5.3 编译与运行编译命令gcc -O2 -o fdfwd fdfwd.c -lm ./fdfwd输出两列时间步序号和接收点波场值。首波到时与接收距离、介质速度的对应关系是立即可以检查的指标源到接收器水平距离 100 格×4 m400 mc2500 m/s理论到时约 0.16 s对应第 267 步。输出会在第 275 步前后出现第一个明显极值。6. 震源到时、边界伪影与数值色散的排查6.1 均匀模型到时验证上面代码的接收器位置已经给出一个天然验证点。理论到时 t r/c 400/2500 0.16 s对应 n267。程序每 25 步输出一次第一个极值应当在 n275 附近出现。如果要做更精确的到时验证把输出改成每步打印再用一个小脚本完成自动比对./fdfwd | awk NR1{first$1} $20{next} {print $1, $2} | head -20首波到达前波场应接近零到达后开始出现完整的三峰波形。首峰位置偏差超过 2Δt 时优先怀疑震源加载顺序或时间层滚动顺序。6.2 三个高频问题快速排查症状最可能原因检查步骤波前出现梳齿状高频振荡空间步长相对有效频率过粗缩小 h 后看波形是否明显变干净记录尾部持续有半正弦隆起吸收带 α 太小或 NB 太薄打印 gx 在最外层的值确认小于 10⁻³到时偏差大于 3Δt 且不收敛震源项顺序或指针滚动有误手推前三个时间步与程序逐值对照6.3 用多次数值导验证网格色散一个容易操作的分辨率验证技巧把接收道写成文件对时间序列做一次二阶导再做一次四阶导和原始波形对比主峰位置。在色散可控的范围内导数波形主峰应与原始波形对齐如果导数波形整体拖尾、旁瓣被拉开说明网格对高频成分的相速度已经明显偏离真实值。这个技巧适合在放大网格时快速找到允许的最大 h。对基本正演程序来说先跑通均匀介质线再逐步加入层状模型始终用理论到时作为基线之后的真实构造差异才可信而不是数值格式变化带来的假象。本文还有配套的精品资源点击获取