简介这份资源面向石油物探方向的研究生及地震数据处理初学者聚焦AVO正演模型实验与地震数据正演这一核心课题。包内共4个cpp源码文件压缩包约12KB均为C实现的正演程序涵盖加噪音条件下的AVO正演模型实验、角度区域处理以及直接形成CDP道集等典型环节代码结构紧凑便于直接编译调试与二次开发。其中部分代码已在中石油相关软件中投入实际应用具备一定的工程参考价值。目前已有203人学习下载说明其在物探教学与科研场景中具有一定认可度。读者可借助这些源码理解AVO正演的基本流程、噪音添加方式与CDP道集生成逻辑对照实验参数观察不同角度区域对正演结果的影响从而加深对地震波场模拟与AVO分析方法的认识适合作为课程实验、课题预研或算法验证的起步素材。1. 从一份“avo正演”压缩包说起地震数据正演到底在算什么拿到一个名为 AVOANDAVAForward 的压缩包很多人第一反应是“里面是不是有能直接跑的代码”但真正值得先弄明白的是avo正演、地震数据正演、AVO 分析、石油物探这几件事是怎么串起来的。简单说地震正演就是给定一个地下介质模型层速度、密度、厚度用物理公式算出“如果在这里放炮、在那里接收应该收到什么样的地震记录”。而 AVO 正演更进一步它关心的是振幅随偏移距的变化——同一反射界面近炮检距和远炮检距的反射振幅为什么不一样这个差异里藏着岩性和流体的信息。这份材料面向的典型读者是石油物探方向的研究生或者刚进入地震资料处理解释岗位、需要自己造数据验证算法的人。它解决的核心问题是你手上没有真实井旁地震道或者真实数据信噪比太低、控制变量太难于是先用正演造出一批“已知答案”的合成记录拿去做反演、做属性分析、做机器学习标签。适合谁适合那些不想只调库、想搞清楚 Zoeppritz 方程到底怎么落到代码里的人。下面我从模型参数怎么定、代码怎么搭、坑在哪一层层拆开讲。2. AVO 正演的理论骨架与最小可跑模型2.1 为什么 AVO 正演绕不开 Zoeppritz 与 Shuey 近似AVO 正演的物理内核是弹性波在层状介质分界面上的反射透射。完整解是 Zoeppritz 方程它给出反射系数随入射角的变化但形式复杂、非线性强直接拿来做反演不现实。工程上更常用的是 Shuey 近似它把反射系数写成截距 A、梯度 B 和曲率项的组合物理意义清晰A 对应垂直入射反射系数B 反映振幅随偏移距的变化率。做正演时如果你只是想快速生成一批 AVO 响应曲线Shuey 足够如果你要验证全波形或精确振幅就得回到 Zoeppritz。我一般会先用 Shuey 搭一个最小模型确认参数单位和角度范围没问题再换 Zoeppritz 做对比。这样做的原因是 Shuey 计算量小、可解析出问题容易定位Zoeppritz 数值求解涉及复数运算和临界角处理一旦量纲错了曲线会直接飞掉排查成本高。2.2 用 Python 搭一个两层 AVO 正演的最小脚本下面这段代码实现两层介质、基于 Shuey 近似的反射系数计算并画出振幅随入射角的变化曲线。它不依赖任何商业软件装好 numpy 和 matplotlib 就能跑。import numpy as np import matplotlib.pyplot as plt def shuey_reflectivity(vp1, vs1, rho1, vp2, vs2, rho2, theta_deg): 两层介质 Shuey 近似 AVO 反射系数 theta_deg: 入射角单位度 返回: 反射系数数组 theta np.deg2rad(theta_deg) # 垂直入射反射系数截距 A A (vp2 * rho2 - vp1 * rho1) / (vp2 * rho2 vp1 * rho1) # 梯度 B这里采用简化形式实际可用更完整的 Shuey 表达式 dvp vp2 - vp1 dvs vs2 - vs1 drho rho2 - rho1 vp_avg (vp1 vp2) / 2.0 vs_avg (vs1 vs2) / 2.0 rho_avg (rho1 rho2) / 2.0 B 0.5 * (dvp / vp_avg) - 2.0 * (vs_avg / vp_avg)**2 * (dvs / vs_avg) \ - 0.5 * (drho / rho_avg) # Shuey 近似R A B * sin^2(theta) R A B * np.sin(theta)**2 return R # 参数示例第一层为页岩第二层为含气砂岩数值为虚构示例仅用于演示 vp1, vs1, rho1 3000.0, 1500.0, 2.40 vp2, vs2, rho2 2600.0, 1400.0, 2.20 angles np.linspace(0, 40, 81) R shuey_reflectivity(vp1, vs1, rho1, vp2, vs2, rho2, angles) plt.figure(figsize(6, 4)) plt.plot(angles, R, b-, linewidth1.5) plt.xlabel(Incident angle (degree)) plt.ylabel(Reflection coefficient) plt.title(AVO response (Shuey approximation)) plt.grid(True) plt.tight_layout() plt.savefig(avo_curve.png, dpi150) plt.show()逻辑说明函数先算垂直入射反射系数 A再算梯度 B最后用 sin²θ 加权得到不同角度的反射系数。参数说明vp1/vs1/rho1 是上层纵波速度、横波速度、密度vp2/vs2/rho2 是下层对应参数theta_deg 是入射角数组。注意速度单位统一用 m/s密度用 g/cm³混用会导致 A 和 B 的量级完全错乱。跑完这张图你应该能看到一条随角度单调或非单调变化的曲线如果是一条水平直线说明 B 算出来接近零检查一下 vs 和 rho 的差异是否被正确传入。2.3 从单界面到多层卷积模型与子波的选择单界面 AVO 曲线只是反射系数真正的地震道是反射系数序列与子波的卷积。做地震数据正演时你需要把每一层的反射系数按时间采样排成序列再与一个零相位或最小相位子波做卷积。常见做法是用雷克子波主频根据目标层深度和分辨率需求选 2540 Hz时间采样率取 1 ms 或 2 ms然后加上随机噪声模拟实际资料。这里有个容易忽略的点AVO 正演通常按角度道集组织也就是说你要对每一个角度单独算反射系数序列再分别卷积。如果只算一个角度然后复制到所有角度那就退化成常规叠后正演失去了 AVO 的意义。我一般会先生成 0° 到 40°、间隔 2° 的角度道集每个角度对应一个反射系数序列最后拼成一个角度道集矩阵行是时间列是角度。3. 把正演数据做成可用的角度道集参数、噪声与格式3.1 角度道集的参数表与常见取值范围做 AVO 正演参数不是随便填的。下面这张表是我在多个模拟项目中总结的常用范围针对的是碎屑岩储层场景碳酸盐岩或深层致密层需要另行调整。参数含义常用范围备注入射角范围角度道集的起止角0°40°超过 40° 后 Shuey 近似误差增大角度间隔相邻角度差1°5°间隔太小计算量大太大丢失 AVO 梯度细节子波主频雷克子波峰值频率2540 Hz深层取低频浅层取高频时间采样率时间轴离散间隔12 ms与目标层厚度匹配避免假频信噪比信号与随机噪声能量比520 dB验证反演鲁棒性时取低值层数模型总层数310 层层数过多时反射相互干涉需正演后检查这张表不是死规定但如果你第一次做建议从 0°35°、间隔 2°、主频 30 Hz、采样率 1 ms 开始这个组合在大多数碎屑岩模型上都能得到清晰的 AVO 响应。3.2 加噪声与角度道集拼接的代码实现下面代码在上一节基础上生成多层反射系数序列按角度循环加噪声后拼成角度道集。import numpy as np def ricker_wavelet(freq, length, dt): 生成雷克子波 t np.arange(-length//2, length//2) * dt pi2f2t2 (np.pi * freq * t) ** 2 w (1 - 2 * pi2f2t2) * np.exp(-pi2f2t2) return w def build_angle_gather(layer_params, angles, dt0.001, wavelet_freq30, snr_db10): layer_params: 列表每个元素为 (vp, vs, rho, thickness_m) 返回: 角度道集矩阵 (时间采样点数, 角度数) n_layers len(layer_params) # 计算每层双程旅行时生成时间轴 times [] t_cum 0.0 for i in range(n_layers - 1): vp layer_params[i][0] thick layer_params[i][3] t_cum 2 * thick / vp times.append(t_cum) total_time t_cum 0.2 # 尾部留白 n_samples int(total_time / dt) reflect_series np.zeros(n_samples) # 在每个界面时间点放置反射系数垂直入射近似用于演示 for i, t in enumerate(times): idx int(t / dt) if idx n_samples: vp1, vs1, rho1 layer_params[i][0], layer_params[i][1], layer_params[i][2] vp2, vs2, rho2 layer_params[i1][0], layer_params[i1][1], layer_params[i1][2] A (vp2 * rho2 - vp1 * rho1) / (vp2 * rho2 vp1 * rho1) reflect_series[idx] A wavelet ricker_wavelet(wavelet_freq, 81, dt) gather np.zeros((n_samples, len(angles))) for j, ang in enumerate(angles): # 这里简化处理每个角度用同一反射序列卷积实际应按角度重算反射系数 trace np.convolve(reflect_series, wavelet, modesame) # 加噪声 noise np.random.randn(n_samples) * np.std(trace) / (10 ** (snr_db / 20.0)) gather[:, j] trace noise return gather, dt # 虚构三层模型页岩-含气砂岩-页岩 layers [ (3000.0, 1500.0, 2.40, 100.0), (2600.0, 1400.0, 2.20, 50.0), (3200.0, 1600.0, 2.50, 200.0), ] angles np.arange(0, 36, 2) gather, dt build_angle_gather(layers, angles, dt0.001, wavelet_freq30, snr_db10) print(角度道集形状:, gather.shape)逻辑说明先根据层厚度和速度算双程旅行时把反射系数放到对应时间采样点再用雷克子波卷积最后按信噪比加高斯噪声。参数说明layer_params 每层四个值顺序是 vp、vs、rho、厚度angles 是角度数组snr_db 控制噪声强度数值越小噪声越强。注意这里为了演示每个角度用了同一反射序列实际 AVO 正演中每个角度应重新计算反射系数否则角度道集上不会有振幅随角度变化的现象。如果你跑出来所有角度道长得一模一样就是这个原因。3.3 输出格式SEG-Y 还是 NumPy 数组正演结果最终要给别人用或者喂给反演算法。如果只是自己验证存成 NumPy 的 .npy 或 .npz 最方便读取快、无精度损失。如果要导入到主流解释软件就得写 SEG-Y。常见做法是用 segyio 或 obspy 写 SEG-Y 文件道头里填好 CDP、偏移距、角度等信息。这里不展开 SEG-Y 的字节级细节但提醒一点角度道集写 SEG-Y 时很多人把角度信息塞在自定义道头里结果对方软件读不出来。稳妥的做法是同时输出一份纯文本的角度-道号对照表别让格式成为协作的瓶颈。4. 避坑与排查AVO 正演里最容易翻车的五个地方4.1 现象AVO 曲线随角度几乎不变梯度接近零原因最常见的是 vs 和 rho 在两层之间差异太小或者干脆填了相同值。Shuey 近似里梯度 B 对 vs 和 rho 的变化敏感如果这两项没差异B 自然接近零。另一个可能是角度范围太小sin²θ 变化不明显。解决检查两层 vs 和 rho 是否确实不同把角度上限提到 35°40°再看曲线是否有斜率。4.2 现象反射系数绝对值大于 1原因Zoeppritz 或 Shuey 在临界角附近会出现全反射反射系数模可以到 1但如果你在垂直入射公式里算出了大于 1 的值多半是密度或速度单位不一致比如上层用 m/s、下层用 km/s。解决统一量纲速度全部用 m/s密度全部用 g/cm³ 或 kg/m³不要混。另外检查两层参数是否出现了“速度倒转”且密度也倒转的极端情况那可能是模型本身不合理。4.3 现象角度道集上同相轴随角度倾斜方向不对原因这通常是反射系数随角度的变化趋势与预期相反根源可能是 Shuey 近似中梯度项的符号写错或者 vs/vp 比值算反。解决拿一个已知的简单模型比如上层低速、下层高速手动算几个角度的反射系数和代码输出对比。如果符号相反检查 B 的表达式里各项的正负号。4.4 现象加噪后 AVO 梯度完全被淹没原因信噪比设得太低或者噪声是加在卷积之后但幅度没有按信号能量归一化。解决先算信号的标准差再按目标信噪比反推噪声标准差不要直接给一个固定噪声幅度。如果是为了测试反演鲁棒性可以从 20 dB 逐步降到 5 dB观察梯度估计在哪一级开始失稳。4.5 现象多层模型反射相互干涉分不清哪个是目标层原因层厚度太薄或者子波主频太低导致相邻反射的波峰波谷叠在一起。解决先算调谐厚度雷克子波主频 f 对应的调谐厚度约为 vp/(2.5f)。如果目标层厚度小于这个值要么提高主频要么接受干涉并做正演记录时标注清楚。我一般会在正演前把层厚度和主频代入这个公式估一下避免做完才发现分辨不开。5. 进阶用法用正演数据验证 AVO 反演与属性交会正演造出来的数据最大的价值不是画几条曲线而是给反演算法提供一个“答案已知”的测试集。你可以把正演得到的角度道集当成输入跑一遍 AVO 反演看反演出来的截距 A 和梯度 B 跟正演时设定的真值差多少。这个闭环能帮你判断反演算法在什么信噪比、什么角度范围下开始失效。具体做法是先生成一批随机层模型每个模型算角度道集同时记录真实的 A 和 B然后把道集喂给反演流程得到估计的 A 和 B最后画交会图横轴真实 A、纵轴估计 A看散点是否落在对角线上。如果偏离严重再按信噪比分组看是噪声问题还是角度覆盖问题。下面是一个简单的交会验证代码片段# 假设 true_A, true_B 是正演时记录的真值数组 # est_A, est_B 是反演估计值数组 import numpy as np import matplotlib.pyplot as plt fig, axes plt.subplots(1, 2, figsize(10, 4)) axes[0].scatter(true_A, est_A, s10, alpha0.6) axes[0].plot([-0.2, 0.2], [-0.2, 0.2], r--) axes[0].set_xlabel(True A) axes[0].set_ylabel(Estimated A) axes[0].set_title(Intercept crossplot) axes[1].scatter(true_B, est_B, s10, alpha0.6) axes[1].plot([-0.5, 0.5], [-0.5, 0.5], r--) axes[1].set_xlabel(True B) axes[1].set_ylabel(Estimated B) axes[1].set_title(Gradient crossplot) plt.tight_layout() plt.savefig(avo_crossplot.png, dpi150) plt.show()参数说明true_A/true_B 来自正演模型参数直接计算est_A/est_B 来自反演输出。红色虚线是对角参考线散点越贴近这条线反演越可靠。如果梯度 B 的散点明显偏离优先检查角度范围是否覆盖了足够大的入射角以及噪声是否在远角度把信号压垮了。还有一个进阶技巧用交会图区分岩性。在正演阶段你可以设定不同流体气、油、水的砂岩模型观察 A-B 交会点的分布区域。含气砂岩通常表现为低 A、负 B 或异常梯度油和水则落在另一片区域。这个分布规律可以直接迁移到实际数据的解释中作为流体识别的先验模板。我自己做这类正演时最大的教训是不要一上来就堆几十层复杂模型先把两层模型的角度道集跑通确认 AVO 曲线的趋势和量纲都对再逐步加层、加噪、加角度。很多翻车都发生在“模型太复杂出错了不知道是哪一层的问题”。另一个习惯是每次正演都保存一份参数日志记录速度、密度、厚度、子波主频、信噪比不然过两天回头看数据根本想不起来当时设了什么。希望帮到你。本文还有配套的精品资源点击获取