简介这份资源面向石油物探方向的研究生及地震数据处理初学者聚焦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, freq35, dt0.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[i1], layer_vs[i1], layer_rho[i1], angles_deg) idx int(layer_t[i1] / dt) for j, r in enumerate(R): gather[j, idx:idxlen(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 (m/s)Vs (m/s)密度 (g/cc)Vp/Vs页岩盖层300015002.402.00含水砂岩280016002.301.75含气砂岩240015502.151.55致密灰岩550030002.651.83这张表是我做正演时的起手模板你可以根据实际工区调整。注意含气砂岩的 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/sf 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/sVs 从 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_samples5000): sampler qmc.LatinHypercube(d5) samples sampler.random(nn_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% 是物理上不可能的模型网络学了一堆垃圾。希望帮到你。本文还有配套的精品资源点击获取