Madagascar下CGFWI全波形反演实践：从RSF到梯度更新

简介：面向地球物理勘探与地震数据处理研究者，这份资源聚焦Madagascar开源平台上的全波形反演（Waveform Inversion），围绕初至波、多次反射与复杂波动信息，解决地下速度模型高精度反演问题。包内共24个文件，以rsf格式的速度/梯度模型数据、vpl可视化文件、Python与C源程序及SConstruct构建脚本为主，压缩包仅380KB，结构紧凑。已有445人学习下载。资源呈现了完整的FWI实验流程：包含Mztzfwi2d.c核心反演程序、bldutil.py与configure.py等辅助脚本、SConstruct自动化编译配置，以及smvel、vel等速度模型和shotcur快照文件，覆盖从数据预处理、初始模型建立、波场模拟到误差函数最小化与参数更新的关键环节。用户可按文档运行并复现二维全波形反演迭代过程，观察模型更新与波形拟合变化，rsf、vpl文件则便于直接查看速度场和反演结果，对系统理解Madagascar下的FWI参数设置、梯度计算与优化策略具有直接参考价值。

1. 为什么一个 rar 里的全波形反演代码值得跑一遍

全波形反演（waveform inversion）这几年几乎成了高分辨率速度建模的代名词，常规走时层析只能给出光滑背景，FWI 则利用地震波的相位、幅值和全频带信息把速度模型推到可以分辨薄层和断块的尺度。Madagascar 平台是一套以 RSF 为统一数据格式的开源地球物理软件包，CGFWI 正好是跑在这个平台上的 2D 全波形反演实现。压缩包不大，没有 GUI、没有一键安装，只有 SConstruct、configure.py、一个 C 源程序和一堆 .rsf 文件，但这套结构恰好把一个最小可复现的 FWI 流程拆得很清楚。适合有 Madagascar 基础的工程师，想用一套能跑通的例子去理解梯度计算和模型更新细节。

2. 在 Madagascar 里解锁 CGFWI：从 configure 到 scons

2.1 解包并检查 RSF 环境

拿到 CGFWI.rar 后，第一件事不是急着打开 README，而是先确认 Madagascar 的 RSF 环境变量已经可用。因为 Mztzfwi2d.c 大量调用 rsf.h，编译时要用 Madagascar 安装目录下的头文件和动态库，如果 RSFROOT 没有设置，scons 会在第一轮依赖扫描阶段就把整个构建停住。

我一般会先做一次环境检查：

unrar x CGFWI.rar cd CGFWI echo "RSFROOT=$RSFROOT" which scons which sfheader

这里unrar x解压并保留内部目录结构；echo看 RSF 根路径；which scons检查构建工具，which sfheader确认 Madagascar 基础命令在当前 PATH 里。如果四个步骤里有哪一个输出为空，说明你还没有把 Madagascar 的环境变量加载进来，常见做法是先执行安装目录下的source env.sh（或你安装时写进 shell 配置的那段脚本），再重新打开终端重复上面的检查。

2.2 configure.py 与 SConstruct 的分工

CGFWI 目录里同时有 configure.py 和 SConstruct，容易被人忽视。configure.py 只做系统探测，把编译器、RSF 路径、编译选项写进配置；SConstruct 才是真正定义编译规则的文件。顺序一旦搞反，scons 会因为缺少配置变量而报一堆 undefined。

这组文件在构建里的角色可以排成一张表：

文件	构建角色	运行角色
configure.py	生成编译配置	一般只在构建前运行
SConstruct	定义 target/dependency	可以通过 scons 触发数据流
Mztzfwi2d.c	编译进主程序	反演主体代码
bldutil.py	scons 辅助函数	不需要直接调用
rsf.h	RSF C 接口头文件	被源文件 include

常规的构建命令是：

python configure.py --rsfroot=$RSFROOT scons -Q -j4

--rsfroot指定 Madagascar 安装位置，configure 把路径写入配置，后续 Mztzfwi2d.c 里的#include "rsf.h"才能找到真实头文件；-Q是 quiet 模式，只显示关键进度，出错时建议去掉-Q看完整 gcc 命令行；-j4指定 4 个并行编译任务。

常见失败还有一种：configure.py 走完了，scons 却在编译某个.c文件时报告找不到 rsf.h。这时别急着改 C 代码，先看 scons 打印出的编译命令里有没有-I$RSFROOT/include。我一般用scons -n干跑一遍，找出哪一步少了 include 路径，再重新执行 configure。

2.3 常见运行方式与 stdin/stdout 约定

构建完成后，目录里会多出一个可执行文件，名字由 SConstruct 控制，常见的是直接叫Mztzfwi2d。Madagascar 程序通常从 stdin 读取 RSF 数据流，结果写到 stdout，CGFWI 的 C 源码也按这个约定写。可以这样试探：

# 先看构建目标 scons -n # 试运行，如果缺少输入文件会报 usage ./Mztzfwi2d < vel.rsf > out.rsf

scons -n是 dry-run，只列出会执行的构建动作；./Mztzfwi2d < vel.rsf把当前速度模型作为输入喂给程序，输出重定向到 out.rsf。如果程序需要额外参数，一般会在 stderr 打出一段 usage。不要直接拿自己的野外数据来跑，先用包内的 vel.rsf 做冒烟测试，确认流程通了再换真实数据。

文件列表里还出现 vel.rsf、shotcur.rsf、grads.rsf，说明 SConstruct 可能把整个流程串成了数据依赖。这样 scons 不只是编译器，还承担了一部分工作流引擎的职责：当输入模型变化，后续目标会自动重算。CGFWI 的价值就在这里，它把一次反演变成可复现的数据流关系，而不是靠人肉记录中间文件。

3. 从 vel.rsf 到 shotcur.rsf：CGFWI 文件清单里的反演流程

3.1 两个速度模型：初始模型和真实模型的角色

FWI 本质是个优化问题：给定初始速度模型，用波动方程正演得到模拟波形，再和观测炮集做差，把残差反投影回模型空间求梯度，更新速度。vel.rsf 是当前速度模型，smvel.rsf 从命名看更接近经过平滑的启动模型。两者在迭代中的作用不同，smvel.rsf 提供平滑起步，避免一开始就引入过多小尺度假象；vel.rsf 负责承接梯度加步长后的更新，每次迭代后都会变化。

这组文件在反演链条里的位置可以这么看：

文件	判断	在反演链条里的位置
smvel.rsf	平滑/初始化模型	反演起点，提供低频框架
vel.rsf	当前速度模型	正演输入，每次迭代被更新
shotcur.rsf	观测炮集记录	目标函数中的观测资料
grads.rsf	梯度场	模型更新方向
addgrads.rsf	累加梯度	用于线搜索或优化器累积
deltav.rsf	扰动模型	反映迭代步长与阻尼

跑反演前要先确认这些 RSF 文件的维度是否一致：

# 检查速度模型的空间维度和采样间隔 sfheader vel.rsf | grep -E "n1|n2|n3|d1|d2|d3" sfheader smvel.rsf | grep -E "n1|n2|n3"

RSF 文件头是 ASCII 文本，n1/n2/n3表示每个维度的采样点数，d1/d2/d3表示采样间隔。两个速度模型 n1、n2 必须一致；不一致时 scons 跑数据流会直接报 dimension mismatch。看头文件时还要注意 label 字段：如果里面写的是 s/km 或 us/ft，那这些文件存的可能是慢度或换算后的速度平方，直接当速度用会把正演主频全部带偏。

3.2 shotcur.rsf 如何驱动残差

shotcur.rsf 是整个反演的目标，也就是观测数据。全波形反演把地震记录里每道整段波形都当作信息源，所以它对模型细节的敏感度远高于初至走时层析。数据残差越接近零，说明当前模型与地下真实速度越匹配。

FWI 的目标函数通常写成 L2 形式，示意如下：

# L2 目标函数示意，obs 为观测炮集，syn 为模拟炮集 objective = 0.5 * np.sum((obs - syn) ** 2)

obs 来自 shotcur.rsf，syn 由当前速度模型正演而来。L2 对异常值非常敏感，所以实际数据里的强噪声会直接污染梯度。如果 CGFWI 内部没有做稳健范数替换，我一般会在喂给程序之前对 shotcur.rsf 做一次带通滤波，把明显不在有效频段的能量切掉，否则反演前几步就会被浅层强能量带偏。

这里要特别强调：shotcur.rsf 应该是经过解编、静校正、滤波后的炮集，而不是原始记录。CGFWI 文件夹里只放一个 shotcur.rsf，说明预处理已经完成。如果你把自己的数据放进来做反演，先要保证观测数据和模拟数据的时间采样点数、炮道几何完全一致，否则残差计算那一步就会因为数组长度对不上直接退出。

3.3 梯度文件与更新链路

grads.rsf 是梯度，但梯度不是直接加到模型上，而是需要先做尺度调整和正则化。addgrads.rsf 的出现说明 CGFWI 内部可能累积了多次梯度，常见用途有两个：一是为 L-BFGS 提供相邻两次迭代的梯度差，二是用来判断收敛——如果累积梯度整体量级不再下降，说明优化已经进入平台区。deltav.rsf 则是速度更新量，可以把它和速度模型相加来看实际修正的大小。

判断梯度是否健康，我用 sfattr 看统计量：

# 用 sfattr 看梯度最大值、最小值和 rms sfattr grads.rsf sfattr addgrads.rsf

只看一个数字还不够。梯度最大值和 rms 的比值如果超过正常范围，说明某个网格点的梯度异常突出，这种点往往是观测数据覆盖不足或震源位置附近有采样假象。遇到这种情况不要继续加大迭代，先回到模型参数化那一步，对梯度做平滑或裁掉覆盖极差的边缘区域。

4. 梯度、步长和正则化：CGFWI 的四组关键参数

4.1 从 L2 残差到速度更新

FWI 的梯度不是直接对目标函数求偏导得到的，而是通过伴随状态法，用一次正传波场和一次反传残差波场在时间轴上互相关求得。CGFWI 的主体是 2D 时间域有限差分正演，每一步迭代里有两笔主要开销：正演波场的存储，以及残差波场的伴随反传。内存不够时，最常见做法是 checkpointing，也就是每隔 N 时间步存一个波场快照，反传时再从快照里恢复。

迭代链路可以抽象成这段伪码：

# 示意性伪代码，对应 CGFWI 的迭代骨骼 for it in range(niter): syn = forward(vel) # 正演模拟波场 res = observed - syn # 观测与模拟残差 grad = backprop(res, vel) # 伴随反传得到梯度 step = line_search(res, grad, vel) # 一维搜索 vel += step * grad # 模型更新

正演负责把速度模型映射到数据空间，反传把数据残差映射回模型空间，步长把梯度变成合理的模型修改量。这个流程里最容易出问题的是 line_search：如果步长直接取固定值，前几步降损耗很快，后面会因为目标函数非线性和梯度量级变化而出现抖动。

在 CGFWI 这个规模下，我一般不用随机梯度。全波形反演对梯度噪声很敏感，L-BFGS 只需要少量内存就能近似海森矩阵的逆。实际跑的时候，前几次迭代用梯度下降热身，再切换到 L-BFGS，比全程 L-BFGS 更稳。原因是 FWI 目标函数非线性强，海森近似在前几步很容易被坏梯度带偏。

4.2 四组参数怎么设

下面这张表，是我拿到类似 CGFWI 的包之后默认会检查的几组参数：

参数	经验范围/做法	影响
震源主频	先 3～5 Hz，再逐级提到 15 Hz	频率越高，能恢复的细节越细，但非线性越强
时间采样间隔 d1	满足 CFL 条件，差分稳定因子控制在 0.8 以下	时间步长太大会数值发散
步长/阻尼系数	先看梯度 rms，取 1e-3 ~ 1e-2	过大发散，过小停滞
正则化权重	从 0.01 倍最大梯度算起	平衡模型平滑度与数据拟合
迭代次数	20～60 次为常见区间	太多过拟合，太少欠拟合

频率从低到高几乎是 FWI 的标准做法。低频对初始模型误差不敏感，先恢复大尺度构造，再用高频细节刻薄边界。如果一开始就用 15 Hz 主频，很容易卡在局部极小值，而且 CFL 条件会迫使时间采样间隔变小，算力成本成倍上升。

正则化方面，如果 CGFWI 的程序内部没有自动做模型平滑，常见做法是在更新前对梯度做一次空间平滑：

# 对梯度做高斯或矩形平滑，减少高频噪声 sfgsmooth < grads.rsf > grad_smooth.rsf rect1=5 rect2=5

rect1 和 rect2 是矩形平滑半径，单位是网格点，一般取 3 到 8。平滑半径太大会把断层边界抹掉，太小又起不到压制假象的作用。还有一点要注意：平滑必须在梯度域做，而不是直接对速度模型做，因为平滑模型会同时影响正演波场和梯度，两边的误差容易被混在一起。

4.3 观察收敛的三个信号

判断 CGFWI 是否真的在收敛，我会同时看三个信号。第一是目标函数值还能不能继续下降；第二是梯度 rms 的变化速率；第三是 deltav.rsf 是否从有规律的层状特征变成随机花斑。如果前两项都在下降，但第三项出现高波数噪声，说明正则化权重偏小。

抓取这些信号不一定要写复杂脚本，直接用 sfattr 就够了：

# 若梯度 rms 数量级超过 1e3，通常先缩小步长再继续 sfattr grads.rsf | grep -E "max|rms"

梯度 rms 的量级取决于炮数、震源幅度和几何扩散，不同数据之间没有绝对可比性。我通常把第一次迭代的 rms 记为基准，后续迭代只要出现比基准大两个数量级的变化，就先停下来检查是不是步长设置出问题。还有一种情况：梯度 rms 一直在下降，但目标函数基本不动，这通常是当前频率已经到底，继续跑只是在拟合噪声。

5. 用 vel.vpl 和梯度曲线给反演结果做体检

5.1 把 vel.rsf 快速变成图

构建跑通之后，最直接的验证是打开 vel.vpl。这个文件是速度模型渲染成剖面图的属性模板，里面可能预先配置好色标、幅值范围或裁剪区间。就算没有这个文件，用 sfvelplot 也能直接出图：

# 纵向剖面方式显示速度模型 sfvelplot < vel.rsf | sfplot & # 或加载已有的属性模板 sfplot < vel.rsf vpl=vel.vpl &

第一行把 vel.rsf 转成 sfplot 能读的格式，再弹出窗口显示；第二行明确指定 vpl 模板。如果 vpl 里定义了速度显示范围，就用那个范围；没有定义时 sfplot 会自动取全局最大最小值，这会让颜色对比度拉满，反而看不清浅层细节。

5.2 一个续跑技巧：给长 FWI 加检查点

反演跑到十几步之后才发散是常有的事。我一般会在每次外部迭代结束后，把当前 vel.rsf 复制到带编号的 checkpoint 目录，避免最后一发不可收拾：

mkdir -p ckpt for i in {1..10} do ./Mztzfwi2d < vel.rsf >> fwi.log cp vel.rsf ckpt/vel_$(date +%H%M%S).rsf done

这段循环每调用一次程序，就保存一份当前模型。日期后缀能让你知道哪一步开始变差。如果后期检查发现某次迭代之后出现了明显随机花斑，就回退到上一个 checkpoint，把步长减半再继续。

5.3 看梯度是否真的稳定下来

判断反演能不能结束，不只看目标函数，也要看梯度 rms 的走势。如果连续三轮sfattr grads.rsf的结果变化不到 1%，这个频段的反演基本做透了，再去提高震源主频继续下一轮。这时候继续硬跑只是把噪声变成模型里的随机层，回退到梯度开始稳定前的那次模型往往会更可信。

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