AVO正演从理论到实践：Zoeppritz方程、Aki-Richards近似与Python实现

简介：这份资源面向石油物探方向的研究生及地震数据处理初学者，聚焦AVO正演模型实验与地震数据正演这一核心课题。包内共4个cpp源码文件，压缩包约12KB，均为C++实现的正演程序，涵盖加噪音条件下的AVO正演模型实验、角度区域处理以及直接形成CDP道集等典型环节，部分代码已在中石油相关软件中投入实际应用，具备工程参考价值。目前已有197人学习下载。读者可借助这些源码理解AVO正演的基本流程与实现思路，掌握加噪音对正演结果的影响、角度域数据的组织方式以及CDP道集的直接生成方法，从而为地震数据正演建模、算法复现与后续处理研究提供可运行的代码基础，适合作为课程实验、课题入门与算法对照的参考资料。

1. AVO 正演到底在算什么：从一份研究生作业说起

如果你手上拿到一个叫AVOANDAVAForward.rar的压缩包，打开一看是几段 Fortran 或 MATLAB 脚本，注释里写着「avo正演」「地震数据正演」，那大概率是石油物探方向研究生的课程作业或课题起步代码。它要干的事其实很具体：给定一套水平层状介质模型，每层的纵波速度、横波速度、密度已知，用 Zoeppritz 方程或其近似式，算出不同入射角下反射界面的反射系数，再和子波褶积，合成出一张 CDP 道集。这张道集就是后续做 AVO 属性分析、烃类检测的输入。换句话说，AVO 正演是「已知地下模型，正着推地震响应」，和反演的方向正好相反。它适合三类人：刚进课题组要跑通第一个合成道集的研究生、需要造标签数据做深度学习的算法工程师、以及想验证自己 AVO 属性提取流程对不对的物探从业者。这一章先把「地震正演」这件事的边界划清楚，后面几章再动手。

2. 从 Zoeppritz 到 Aki-Richards：选哪个公式决定你的道集长什么样

2.1 精确解和近似解的取舍

Zoeppritz 方程是 AVO 正演的理论基石，它给出平面波在两种弹性介质分界面上，反射系数随入射角变化的精确解。四个边界条件——位移连续、应力连续——联立出四个方程，解出反射和透射的 P 波、S 波系数。问题是这个解太复杂，写出来是一堆含三角函数的有理式，物理意义不直观，而且对速度密度的小扰动不敏感，做属性分析时不好用。

所以实际正演里，绝大多数人用的是近似式。Aki-Richards 公式是最常见的一个，它把反射系数写成入射角 θ 的函数：

R(θ) ≈ (1/2)(ΔVp/Vp)(1/cos²θ) - 4(Vs/Vp)²(ΔVs/Vs)sin²θ + (1/2)(Δρ/ρ)(1 - 4(Vs/Vp)²sin²θ)

其中 ΔVp、ΔVs、Δρ 是界面两侧的速度密度差，Vp、Vs、ρ 是平均值。这个式子把 AVO 响应拆成了三个部分：零偏移距项、梯度项、曲率项。Shuey 进一步化简，把 R(θ) 写成截距 P 加梯度 G 乘 sin²θ 的形式，这就是后来 AVO 属性分析里 P、G 交会图的来源。

选哪个？我的经验是：做方法验证、写论文对比精确解和近似解误差时，用 Zoeppritz；做合成道集喂给反演或神经网络时，用 Aki-Richards 或 Shuey，因为快，而且和后续属性提取的假设一致。如果你用精确解生成道集，再用近似式去反演，误差里混着近似误差，说不清楚是谁的锅。

2.2 用 Python 实现 Aki-Richards 正演的最小代码

下面这段代码是我一般用来快速验证模型响应的最小实现，输入是上下两层介质参数和入射角数组，输出反射系数曲线。

import numpy as np def aki_richards(vp1, vs1, rho1, vp2, vs2, rho2, angles_deg): """ Aki-Richards 近似计算反射系数 vp1, vs1, rho1: 上层纵波速度(m/s)、横波速度(m/s)、密度(g/cc) vp2, vs2, rho2: 下层参数 angles_deg: 入射角数组，单位度 返回: 反射系数数组 """ theta = np.radians(angles_deg) vp = (vp1 + vp2) / 2.0 vs = (vs1 + vs2) / 2.0 rho = (rho1 + rho2) / 2.0 dvp = vp2 - vp1 dvs = vs2 - vs1 drho = rho2 - rho1 term1 = 0.5 * (dvp / vp) / (np.cos(theta)**2) term2 = -4.0 * (vs / vp)**2 * (dvs / vs) * (np.sin(theta)**2) term3 = 0.5 * (drho / rho) * (1 - 4.0 * (vs / vp)**2 * np.sin(theta)**2) return term1 + term2 + term3 # 示例：砂岩页岩界面，含气砂岩速度降低 angles = np.arange(0, 45, 5) R = aki_richards(3000, 1500, 2.4, 2600, 1300, 2.2, angles) for a, r in zip(angles, R): print(f"入射角 {a:2d}° 反射系数 {r:+.4f}")

这段代码里几个参数需要说明。vp1/vs1/rho1是上层，通常代表页岩盖层；vp2/vs2/rho2是下层储层。含气砂岩的典型特征是 Vp 明显降低、Vs 变化小、密度略降，所以你会看到反射系数随入射角增大而变得更负——这就是所谓的第三类 AVO 异常。angles_deg一般取 0 到 40 度，超过 40 度近似误差会变大。np.cos(theta)**2在零角度时为 1，大角度时迅速增大，这是近似的固有特性，实际处理中远角道集信噪比低，通常截断在 35 到 40 度。

跑完这段，你会得到一条 R-θ 曲线。如果曲线从正变负、或者负值越来越负，说明模型有 AVO 异常。但这只是单个界面，真正的地震道集需要把多个界面的反射系数和子波褶积，再按角度排列成道集。

2.3 从反射系数到合成道集：子波褶积和角度道集排列

单个界面的反射系数只是一个数，地震记录是一个时间序列。要把反射系数变成地震道，需要和地震子波做褶积。常用的是 Ricker 子波，主频一般取 30 到 40 Hz，对应常规地震资料的主频范围。

def ricker_wavelet(freq, length, dt): """生成 Ricker 子波 freq: 主频(Hz) length: 采样点数 dt: 采样间隔(s) """ t = np.arange(length) * dt - (length * dt) / 2 pi2f2t2 = (np.pi * freq * t) ** 2 return (1 - 2 * pi2f2t2) * np.exp(-pi2f2t2) def build_angle_gather(layer_vp, layer_vs, layer_rho, layer_t, angles_deg, freq=35, dt=0.001): """ 多层模型合成角度道集 layer_vp/vs/rho: 每层参数列表，长度 n layer_t: 每层顶界面的双程旅行时列表，长度 n angles_deg: 角度数组 返回: 道集矩阵 (n_angles, n_samples) """ n_layers = len(layer_vp) n_samples = int(layer_t[-1] / dt) + 200 wavelet = ricker_wavelet(freq, 81, dt) gather = np.zeros((len(angles_deg), n_samples)) for i in range(n_layers - 1): R = aki_richards(layer_vp[i], layer_vs[i], layer_rho[i], layer_vp[i+1], layer_vs[i+1], layer_rho[i+1], angles_deg) idx = int(layer_t[i+1] / dt) for j, r in enumerate(R): gather[j, idx:idx+len(wavelet)] += r * wavelet return gather

这里layer_t是每个界面反射波的双程旅行时，需要根据层厚度和速度算出来。n_samples留了 200 个点的尾巴防止子波被截断。idx是界面在时间轴上的位置，每个角度的反射系数乘上同一个子波，叠加到对应位置。最终gather的每一行是一个角度的道，列是时间采样。把 gather 画出来就是一张角度道集图，横轴角度、纵轴时间、颜色表示振幅。

提示：子波长度一般取 81 或 101 个点，太短会截断旁瓣，太长会拖尾干扰深层反射。主频根据你的目标层深度和分辨率需求调，浅层用高频，深层用低频。

3. 模型参数怎么设：速度密度从哪来、角度范围怎么定

3.1 用测井曲线还是经验公式

做 AVO 正演，模型参数是命根子。最理想的情况是有一口井的纵波速度、横波速度、密度曲线，直接读出来做层状简化。但很多研究生手里没有实测横波曲线，这时候就得用经验公式估算。常见的有 Castagna 泥岩线：

Vp = 1.16 * Vs + 1360 (m/s)

或者 Gardner 公式从速度估密度：

ρ = 0.23 * Vp^0.25 (g/cc, Vp 单位 ft/s)

这些公式有适用条件，Castagna 适用于水饱和碎屑岩，Gardner 适用于常规沉积岩。如果你做的是碳酸盐岩或者含气层，经验公式误差会很大，这时候宁可用岩石物理模型（比如 Gassmann 流体替换）去算，也不要硬套。

我一般会建一个表格，把每层的 Vp、Vs、密度、厚度列清楚，再检查一下 Vp/Vs 比值是否合理。砂岩的 Vp/Vs 一般在 1.6 到 1.8，页岩在 1.8 到 2.0，含气砂岩可能低到 1.5。如果算出来 Vp/Vs 是 1.2 或者 2.5，那八成是参数填错了。

这张表是我做正演时的起手模板，你可以根据实际工区调整。注意含气砂岩的 Vp 比含水砂岩低了 400 m/s，但 Vs 只降了 50 m/s，密度降了 0.15，这就是 AVO 异常的来源。

3.2 角度范围、子波主频和采样率

角度范围不是随便定的。常规海上拖缆最大入射角能到 40 到 45 度，陆上可控震源可能只有 30 到 35 度。你做正演时如果取到 50 度，合成道集在远角部分会失真，因为 Aki-Richards 近似在大角度误差急剧增大。我的习惯是最大取 40 度，步长 5 度，这样得到 9 个角度的道集，足够做 AVO 属性拟合。

子波主频决定分辨率。主频 35 Hz、采样率 1 ms 是常规配置。如果你要模拟薄层调谐，主频可以提到 50 Hz，采样率 0.5 ms。但要注意，主频越高，子波旁瓣越明显，合成道集上会出现假的同相轴，别把它当成真实反射。

采样率的选择要满足 Nyquist 定理，1 ms 采样对应 500 Hz Nyquist 频率，远高于地震信号带宽，没问题。但如果你做的是高频正演，比如 100 Hz 主频，采样率至少 0.5 ms。

3.3 层厚和调谐效应

层厚小于子波波长四分之一时，顶底反射会干涉，形成调谐。调谐效应会让振幅和 AVO 梯度都发生变化，如果你用调谐后的道集去反演，得到的阻抗和真实值有偏差。做正演时，如果目标层很薄，要么把层厚设得足够大避开调谐，要么就专门研究调谐对 AVO 的影响。

我一般会先算一下子波的主波长：λ = Vp / f。比如 Vp 3000 m/s，f 35 Hz，λ 约 86 m，四分之一波长约 21 m。如果储层厚度小于 21 m，就要小心调谐。这时候可以做一个层厚扫描，从 5 m 到 50 m 变化，看 AVO 梯度的变化趋势，找到稳定区间。

4. 跑通正演后怎么验证：三个检查点和两个对比实验

4.1 检查点一：零角度反射系数是否等于波阻抗差

零角度时，Aki-Richards 退化为：

R(0) = (ρ2Vp2 - ρ1Vp1) / (ρ2Vp2 + ρ1Vp1)

也就是波阻抗差除以波阻抗和。你可以在代码里加一行，把 angles 设为 0，看输出是否等于手算的波阻抗反射系数。如果不等，检查公式实现有没有漏项或者符号错误。这是最基本的自检，但很多人跳过，结果后面道集极性反了都不知道。

4.2 检查点二：道集同相轴是否随角度变化

合成道集画出来后，看目标层对应的同相轴。如果振幅随角度不变，说明你的反射系数计算里角度项没起作用，可能是np.radians忘了加，或者sin²θ写成了sinθ。如果振幅随角度变化但趋势不对，比如含水砂岩应该振幅减小，结果反而增大，那可能是 Vp/Vs 比值设反了。

4.3 检查点三：和精确 Zoeppritz 解对比

找一个公开的 Zoeppritz 实现，或者自己写一个，把同一组模型参数输入，对比近似解和精确解在 0 到 40 度的差异。一般来说，Aki-Richards 在 30 度以内误差小于 5%，40 度时可能到 10%。如果你发现误差超过 20%，检查一下速度对比度是不是太大——近似式假设速度差远小于平均速度，如果上下层速度差了一倍，近似就失效了。

4.4 对比实验：含气与含水砂岩的 AVO 响应

这是最直观的验证。用同一套骨架参数，只把孔隙流体从水换成气，看道集变化。含水砂岩的反射系数随角度可能变化不大，含气砂岩则会出现明显的振幅增大或极性反转。如果你做出来的含气道集和含水道集几乎一样，那说明流体替换没做对，或者参数里 Vs 没跟着变。

4.5 对比实验：不同子波主频对 AVO 梯度的影响

用 25 Hz、35 Hz、45 Hz 三个主频分别合成道集，提取 AVO 梯度，看梯度值是否稳定。如果主频变化导致梯度大幅波动，说明调谐效应严重，你的层厚可能太薄，或者子波旁瓣干扰了反射系数提取。这时候要么加厚层，要么在提取属性前做谱白化。

5. 避坑与排查：AVO 正演里最容易翻车的五个地方

5.1 道集极性反转但没发现

现象：合成道集上目标层振幅随角度从负变正，你以为这是第三类 AVO，结果检查发现是反射系数符号搞反了。 原因：Aki-Richards 公式里 ΔVp 定义为下层减上层，如果你写成上层减下层，整个曲线极性就反了。 解决：在代码里固定dvp = vp2 - vp1，并在零角度检查波阻抗差符号。如果上层波阻抗大于下层，反射系数应为负，道集上表现为波峰还是波谷取决于子波极性，但相对关系要对。

5.2 角度单位混用

现象：反射系数曲线形状怪异，大角度时数值爆炸。 原因：np.sin和np.cos接受弧度，但你传进去的是角度值，35 度当成 35 弧度算，结果完全不对。 解决：在函数入口统一用np.radians转换，或者在参数名里写明angles_deg，调用时检查。

5.3 子波采样率和道集采样率不一致

现象：褶积后道集同相轴变宽或变窄，时间厚度对不上。 原因：子波的dt和道集的dt不一致，比如子波用 1 ms 生成，道集用 2 ms 采样，褶积时没有重采样。 解决：生成子波和构建道集用同一个dt，或者在褶积前用scipy.signal.resample统一采样率。

5.4 层厚设得太薄导致调谐

现象：目标层顶底反射分不开，AVO 梯度随层厚剧烈变化。 原因：层厚小于四分之一波长，顶底反射干涉。 解决：先算主波长，确保层厚大于四分之一波长。如果实际储层就是薄，那就把正演目的改成研究调谐效应，而不是提取真实 AVO 属性。

5.5 忽略横波速度的流体敏感性

现象：含水换含气后，Vp 降了，但 Vs 没变，AVO 异常不明显。 原因：Gassmann 流体替换中，Vs 对流体不敏感，但 Vp 和密度敏感。如果你只改 Vp 不改密度，反射系数变化不够。 解决：用 Gassmann 方程同时计算 Vp、Vs、密度变化，或者至少按经验把密度也调低 0.1 到 0.2 g/cc。

6. 进阶技巧：用 AVO 正演造深度学习训练集

如果你跑通了单个模型的正演，下一步很可能是批量生成道集，用来训练神经网络做 AVO 反演或流体识别。这时候单条曲线的手工操作就不够了，需要参数化扫描。

我一般会定义一个参数空间：Vp 从 2200 到 3200 m/s，Vs 从 1200 到 1800 m/s，密度从 2.0 到 2.5 g/cc，层厚从 10 到 50 m，子波主频从 25 到 45 Hz。用拉丁超立方采样抽 5000 组，每组生成一个角度道集，标签是对应的 Vp、Vs、密度、流体类型。生成脚本的核心循环和前面一样，只是外面套一层采样。

from scipy.stats import qmc def generate_dataset(n_samples=5000): sampler = qmc.LatinHypercube(d=5) samples = sampler.random(n=n_samples) # 映射到参数范围 vp_shale = 3000 + samples[:, 0] * 200 vp_sand = 2200 + samples[:, 1] * 1000 vs_sand = 1200 + samples[:, 2] * 600 rho_sand = 2.0 + samples[:, 3] * 0.5 thickness = 10 + samples[:, 4] * 40 # 对每组参数调用 build_angle_gather # 保存道集和标签

这里用拉丁超立方而不是均匀网格，是因为 5 维均匀网格点数会爆炸，拉丁超立方能用更少样本覆盖更均匀。生成 5000 组大概需要几分钟到十几分钟，取决于采样点数和角度数。保存成 npy 或 hdf5 格式，训练时直接读。

注意：生成训练集时，角度范围要和实际资料一致。如果你用 0 到 40 度训练，实际资料只有 0 到 30 度，网络在远角部分会外推，误差不可控。另外，子波主频也要和实际资料匹配，否则网络学到的是子波特征而不是 AVO 特征。

还有一个技巧是加噪声。合成道集太干净，网络会过拟合。我一般加 5% 到 10% 的高斯噪声，或者按实际资料的信噪比加。加噪后再做 AVO 属性提取，看梯度是否稳定，如果噪声一加梯度就乱飞，说明你的正演参数太理想，实际资料更差。

最后说一个我自己的习惯：每次生成完数据集，随机抽 10 个道集画出来，肉眼扫一遍。如果看到某个道集同相轴断裂、振幅异常大、或者时间轴对不齐，大概率是某组参数越界了。比如 Vp 小于 Vs，或者密度为负，这些在采样时就要卡住边界。别小看这一步，我见过有人生成了几万条道集，训练完才发现里面有 30% 是物理上不可能的模型，网络学了一堆垃圾。希望帮到你。

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