简介：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(x+h) - P(x+2h) ] / (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(x+h) | u(x+2h) | 截断误差 |
|---|---|---|---|---|---|---|
| 二阶中心差分 | 0 | 1 | -2 | 1 | 0 | O(h²) |
| 四阶中心差分 | -1/12 | 4/3 | -5/2 | 4/3 | -1/12 | O(h⁴) |
用 -1/12、4/3、-5/2、4/3、-1/12 去乘相邻数组值，再统一除以 12，就是在代码里常见的(-u[i-2] + 16*u[i-1] - 30*u[i] + 16*u[i+1] - u[i+2]) / 12.0。这个写法只关心 x 方向，z 方向完全同理。
提示：这里的“四阶”只指空间导数精度，不是时间积分精度。四阶龙格库塔是另一套时间推进算法，不能与空间四阶差分混用；对基本正演程序，时间方向保留二阶即可，实现简单且稳定性条件明确。
2.2 时间二阶推进与 CFL 稳定参数
完整时间推进采用：
P^(n+1)(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_max
c_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[i+1] - u[i+2]) / 12.0; }调用前要注意数组边界：i 从 2 到 n-3，否则会越界读取。实际二维程序里，x 和 z 两个方向各算一次，再把结果相加；这也是第 5 章完整代码里 Laplacian 的计算方式。
3. 边界吸收与网格裁剪：Cerjan 阻尼带的核心设置
3.1 为什么零边界和简单衰减会毁掉合成记录
正演网格必然是有限范围。把四周边界直接当成零位移，会在边界和四角持续产生强反射；因为四阶模板在靠近边界处无法照常使用，边界反射会带着明显的锯齿形态向内部传播。地震正演中这类伪反射能量往往比真实弱反射层信号还强，直接掩盖有效同相轴。
基本程序里最常见的处理是 Cerjan 阻尼带，也叫海绵吸收边界。原理很简单：在计算区域四周包裹一圈衰减带，每个时间步推进后，让波场乘上一个随距离边界衰减的系数，波传进吸收带后能量逐步消失，到达硬边界时已经小到可以忽略。
3.2 Cerjan 阻尼带公式与参数组合
衰减系数定义为到吸收带内壁的距离 d 的函数：
g(d) = exp( -α² · d² )
内部区域 d = 0，g = 1，波场不受影响；越靠近最外层边界，d 越大，衰减越强。实现时把 x、z 两个方向分别做一维系数，再用 gx[i] · gz[j] 相乘。四角区域同时被两个方向的系数压制，反射抑制效果比单方向叠加好得多。
| 参数 | 经验取值 | 说明 |
|---|---|---|
| NB，吸收带厚度 | 20～40 格 | 至少要覆盖几个网格波长 |
| α，衰减强度 | 0.05～0.12 / 格 | 过大会在带内形成硬边界，过小吸收不足 |
| 最外层衰减量 | exp( -(α·NB)² ) | 建议小于 10⁻²，稳定时取 10⁻³ |
这两个参数是耦合的。NB 取 20 时，α 取 0.08，最外层衰减约为 exp(-2.56) ≈ 0.077，余量偏大；NB 取 30、α 取 0.08 时，最外层约为 0.003，能压住大部分反射。更稳妥的做法是先设 NB=30，跑一个均匀模型观察边界反射残差，再在 0.05～0.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*NX+i] *= gx[i] * gz[j]。由于带内 d=0，内部点乘数为 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 的经验要求大约在 6～8 个点。按此设计网格步长：
h ≤ c_min / ( 2.5 · f₀ · Nppw )
这里用低速介质的最小速度 c_min，是为了保证最短波长有足够网格分辨率。之后再用全局最大速度 c_max 约束时间步：
Δt ≤ 0.6 · h / c_max
顺序不能反。如果先用 c_max 定 h，低速区波长变短，很可能在浅部低速层产生明显色散。
举一个可复算的例子：
| 介质情况 | 速度 | 选取参数 |
|---|---|---|
| 低速目标层 | 1800 m/s | f₀=25 Hz，Nppw=8，h ≤ 1800/(2.5×25×8) = 3.6 m，取 4 m |
| 高速基底 | 3500 m/s | c_max=3500，由 h=4 得 Δt ≤ 0.6×4/3500 ≈ 0.686 ms |
实际取 Δt=0.6 ms，留出稳定余量。这个组合下，Ricker 子波在 1800 m/s 低速层的每波长采样约为 1800/(25×4)=18 点，在高速层约为 35 点，色散可以压得比较干净。
子波时长建议覆盖 3～5 个主周期，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 指向旧 p1，p1 指向旧 p2，p2 指向旧 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[ix+1] - p1[ix+2]) / 12.0; double lap_z = (-p1[ix-2*NX] + 16.0*p1[ix-NX] - 30.0*p1[ix] + 16.0*p1[ix+NX] - p1[ix+2*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 m=400 m，c=2500 m/s，理论到时约 0.16 s，对应第 267 步。输出会在第 275 步前后出现第一个明显极值。
6. 震源到时、边界伪影与数值色散的排查
6.1 均匀模型到时验证
上面代码的接收器位置已经给出一个天然验证点。理论到时 t = r/c = 400/2500 = 0.16 s，对应 n=267。程序每 25 步输出一次，第一个极值应当在 n=275 附近出现。
如果要做更精确的到时验证，把输出改成每步打印，再用一个小脚本完成自动比对：
./fdfwd | awk 'NR==1{first=$1} $2<0{next} {print $1, $2}' | head -20首波到达前波场应接近零，到达后开始出现完整的三峰波形。首峰位置偏差超过 2Δt 时，优先怀疑震源加载顺序或时间层滚动顺序。
6.2 三个高频问题快速排查
| 症状 | 最可能原因 | 检查步骤 |
|---|---|---|
| 波前出现梳齿状高频振荡 | 空间步长相对有效频率过粗 | 缩小 h 后看波形是否明显变干净 |
| 记录尾部持续有半正弦隆起 | 吸收带 α 太小或 NB 太薄 | 打印 gx 在最外层的值，确认小于 10⁻³ |
| 到时偏差大于 3Δt 且不收敛 | 震源项顺序或指针滚动有误 | 手推前三个时间步，与程序逐值对照 |
6.3 用多次数值导验证网格色散
一个容易操作的分辨率验证技巧：把接收道写成文件，对时间序列做一次二阶导，再做一次四阶导，和原始波形对比主峰位置。在色散可控的范围内，导数波形主峰应与原始波形对齐；如果导数波形整体拖尾、旁瓣被拉开，说明网格对高频成分的相速度已经明显偏离真实值。
这个技巧适合在放大网格时快速找到允许的最大 h。对基本正演程序来说，先跑通均匀介质线，再逐步加入层状模型，始终用理论到时作为基线；之后的真实构造差异才可信，而不是数值格式变化带来的假象。
本文还有配套的精品资源，点击获取
