简介面向地球物理勘探与地震数据处理研究者这份资源聚焦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 到 scons2.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.pyscons 辅助函数不需要直接调用rsf.hRSF 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 数据流结果写到 stdoutCGFWI 的 C 源码也按这个约定写。可以这样试探# 先看构建目标 scons -n # 试运行如果缺少输入文件会报 usage ./Mztzfwi2d vel.rsf out.rsfscons -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.rsfCGFWI 文件清单里的反演流程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|n3RSF 文件头是 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.rsfsyn 由当前速度模型正演而来。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 的包之后默认会检查的几组参数参数经验范围/做法影响震源主频先 35 Hz再逐级提到 15 Hz频率越高能恢复的细节越细但非线性越强时间采样间隔 d1满足 CFL 条件差分稳定因子控制在 0.8 以下时间步长太大会数值发散步长/阻尼系数先看梯度 rms取 1e-3 ~ 1e-2过大发散过小停滞正则化权重从 0.01 倍最大梯度算起平衡模型平滑度与数据拟合迭代次数2060 次为常见区间太多过拟合太少欠拟合频率从低到高几乎是 FWI 的标准做法。低频对初始模型误差不敏感先恢复大尺度构造再用高频细节刻薄边界。如果一开始就用 15 Hz 主频很容易卡在局部极小值而且 CFL 条件会迫使时间采样间隔变小算力成本成倍上升。正则化方面如果 CGFWI 的程序内部没有自动做模型平滑常见做法是在更新前对梯度做一次空间平滑# 对梯度做高斯或矩形平滑减少高频噪声 sfgsmooth grads.rsf grad_smooth.rsf rect15 rect25rect1 和 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 vplvel.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%这个频段的反演基本做透了再去提高震源主频继续下一轮。这时候继续硬跑只是把噪声变成模型里的随机层回退到梯度开始稳定前的那次模型往往会更可信。本文还有配套的精品资源点击获取