★ FEATURED ARTICLE 简介这份资源面向石油物探方向的研究生及地震数据处理初学者聚焦AVO正演模型实验与地震数据正演这一核心技能点帮助读者理解含噪条件下AVO响应特征及CDP道集的直接生成过程。压缩包内共4个文件全部为cpp源代码整体约12KB涵盖加噪音正演、角度区域分析、b5-25参数实验以及直接形成CDP等不同实验分支代码结构紧凑便于逐文件对照阅读与二次修改。其中部分代码已在中石油相关软件中实际应用具备一定的工程参考价值。目前已有203人学习下载说明该方向存在稳定的学习需求。读者可从中获得AVO正演的实现思路、噪音添加与角度域处理的具体写法以及从模型到CDP道集的完整代码路径适合作为课程实验、课题预研或算法复现的起点也可用于对比不同参数设置下的正演结果差异为后续反演与解释工作打基础。1. 从一份正演代码包说起地震波场模拟到底在算什么第一次拿到“AVO 正演”这个方向的任务时很多人以为不过是套个 Zoeppritz 公式算几条反射系数曲线。真跑起来才发现从弹性参数模型到叠前角道集中间隔着一整套波场模拟的工程链路。所谓地震数据正演本质是给定地下介质的纵波速度、横波速度、密度通过求解波动方程或近似公式算出检波器位置能记录到的地震响应。它既是理解 AVO 响应特征的“黑匣子”拆解工具也是做反演、做储层预测绕不开的第一步。石油物探方向的研究生几乎都要在这上面花掉几个月。这篇笔记就按我实际跑通一套 AVO 正演流程的顺序把选型、参数、代码和踩过的坑讲清楚让后来的人少走弯路。2. 正演方法怎么选从 Zoeppritz 到波动方程的三条路线2.1 三种正演粒度的适用边界做 AVO 正演第一件事不是写代码而是想清楚你要的“正演”到底精细到什么程度。常见做法分三档。第一档是解析解路线直接用 Zoeppritz 方程或其近似式Aki-Richards、Shuey 等算反射系数随入射角的变化。它假设界面上下是均匀各向同性介质、平面波入射算得飞快适合做 AVO 响应规律分析、岩性替换、流体替换这类需要成千上万次试算的场景。缺点是它只给反射系数不给波形不含调谐、不含绕射、不含层间多次波。第二档是射线追踪加褶积路线。先用射线追踪算出各层的走时和入射角再把反射系数序列与子波褶积得到合成道集。它比纯解析解多了子波和走时信息能模拟薄层调谐计算量仍然可控是很多研究生做叠前正演的主力方案。第三档是全波场数值模拟用有限差分求解声波或弹性波方程。它能自然包含绕射、多次波、转换波、频散等所有波动现象最接近真实地震记录但计算量大、边界处理麻烦、参数敏感。做 AVO 机理研究时前两档通常够用要验证复杂构造或强横向变速下的响应才上第三档。选型判断很简单如果你的目标是“看反射系数随角度的变化规律”走第一档如果要“合成出能拿去和实际道集对比的叠前数据”走第二档如果研究的是波场传播本身走第三档。别一上来就上有限差分那是拿大炮打蚊子调试成本会让你怀疑人生。2.2 用 Python 跑通 Zoeppritz 正演的最小实现下面这段代码是我常用的 Zoeppritz 精确解实现输入上下介质的纵波速度、横波速度、密度和入射角输出反射系数。它不依赖任何地球物理专用库只用 NumPy方便你嵌到自己的流程里。import numpy as np def zoeppritz(vp1, vs1, rho1, vp2, vs2, rho2, theta1): 精确 Zoeppritz 方程求解 PP 波反射系数 theta1: 入射角单位弧度可为标量或数组 返回: PP 反射系数 theta1 np.asarray(theta1, dtypefloat) # 由 Snell 定律求透射角 p np.sin(theta1) / vp1 # 射线参数 theta2 np.arcsin(np.clip(p * vp2, -1, 1)) phi1 np.arcsin(np.clip(p * vs1, -1, 1)) phi2 np.arcsin(np.clip(p * vs2, -1, 1)) # 角度相关项 a rho2 * (1 - 2 * np.sin(phi2)**2) - rho1 * (1 - 2 * np.sin(phi1)**2) b rho2 * (1 - 2 * np.sin(phi2)**2) 2 * rho1 * np.sin(phi1)**2 c rho1 * (1 - 2 * np.sin(phi1)**2) 2 * rho2 * np.sin(phi2)**2 d 2 * (rho2 * vs2**2 - rho1 * vs1**2) # 构造矩阵并求解 M np.zeros((4, 4) theta1.shape) # 这里省略完整矩阵填充实际使用时按标准 Zoeppritz 矩阵逐项赋值 # 关键矩阵元素含 sin/cos 角度项注意广播维度 # 求解后取反射系数 Rpp # Rpp ... return None # 占位完整实现见下文说明上面这段是骨架重点在参数含义和维度处理。vp1/vs1/rho1是上层介质参数vp2/vs2/rho2是下层theta1是入射角数组。实际填充矩阵时四个方程的系数分别对应位移和应力的连续性条件矩阵元素里全是sin、cos和密度、速度的组合。新手最容易翻车的地方是角度单位——NumPy 的三角函数吃弧度你从 Excel 或文本读进来的是度忘了转就会得到一条完全错误的曲线而且它不报错只是数值离谱属于典型的“玄学 bug”。参数设置上速度单位统一用 m/s密度用 g/cm³ 或 kg/m³ 都行只要上下层一致。入射角一般取 0 到 40 度步长 1 度足够覆盖常规叠前角道集范围。如果你要做宽角反射可以到 60 度但 Zoeppritz 在临界角附近会出现复数需要额外处理建议先避开。2.3 从反射系数到合成角道集子波与褶积拿到反射系数曲线只是半成品要变成“地震数据”还得乘子波、加走时。常见做法是选一个零相位雷克子波主频根据目标层深度和分辨率定浅层用 40 到 50 Hz深层用 20 到 30 Hz。def ricker(freq, length, dt): 生成雷克子波 freq: 主频 Hz length: 采样点数 dt: 采样间隔 s t np.arange(length) * dt - (length - 1) * dt / 2 pi2f2t2 (np.pi * freq * t) ** 2 return (1 - 2 * pi2f2t2) * np.exp(-pi2f2t2) def build_gather(rc_series, angles, wavelet): 把反射系数序列与子波褶积生成角道集 rc_series: 形状 (n_layers, n_angles) n_angles len(angles) gather np.zeros((len(wavelet) rc_series.shape[0] - 1, n_angles)) for i in range(n_angles): gather[:, i] np.convolve(rc_series[:, i], wavelet, modefull) return gatherricker里freq决定分辨率length一般取 2 到 3 个周期dt是时间采样率常用 1 ms 或 2 ms。build_gather对每个角度道单独褶积输出就是叠前角道集。注意np.convolve的modefull会让输出变长后续做显示时要裁掉两端否则道集上下会多出无意义的拖尾。这一步的坑在于子波长度和采样率不匹配——子波太长会拖慢计算太短则主瓣没采全合成记录会出现旁瓣假象。3. 弹性参数模型怎么建从测井曲线到层状模型3.1 测井数据的整理与去噪正演的输入是弹性参数模型最常用的来源是测井曲线。实际拿到的纵波速度、横波速度、密度曲线往往采样率不一致、有异常值、有缺失段。我一般先做三件事统一采样率到 0.5 m 或 1 m用中值滤波去掉明显毛刺对缺失段做线性插值或按岩性分段赋值。import numpy as np from scipy.signal import medfilt def clean_curve(curve, kernel5): 中值滤波去毛刺kernel 必须为奇数 return medfilt(curve, kernel_sizekernel) def resample_curve(depth_old, curve_old, depth_new): 按新深度轴重采样 return np.interp(depth_new, depth_old, curve_old)medfilt的kernel别开太大5 到 7 点足够开大了会把薄层信息抹掉正演出来的道集会“糊”。重采样用np.interp是线性插值对速度曲线够用如果要做精细的层序分析可以换样条插值但要注意过冲。3.2 层状模型的抽层原则测井曲线是连续采样的但 Zoeppritz 和射线追踪需要层状模型。抽层就是决定哪些深度段合并成一层。原则是速度或密度发生明显台阶的地方必须分层薄层如果厚度小于子波主波长的一半可以合并因为它对波形的影响本来就有限。我通常用速度曲线的梯度做辅助判断梯度超过阈值的深度点作为候选界面再人工检查一遍去掉明显由噪声引起的假界面。抽完层后每层赋一个平均速度、平均密度界面深度取中点。这一步没有标准答案同一个测井曲线不同人抽出来的层数可能差一倍正演结果也会有差异。所以做对比实验时一定要固定抽层方案否则你分不清结果变化是来自参数还是来自分层。3.3 流体替换与岩性替换的快速实现AVO 正演经常要做流体替换比如把含水砂岩换成含气砂岩看反射系数怎么变。常用 Gassmann 方程。def gassmann(vp_sat, vs_sat, rho_sat, rho_fl_new, k_fl_new, k_min, k_fl_old, phi): Gassmann 流体替换返回新的 vp, vs, rho 参数均为国际单位制速度 m/s密度 kg/m3模量 Pa # 由饱和岩石速度反算干岩石模量 mu rho_sat * vs_sat**2 k_sat rho_sat * vp_sat**2 - 4/3 * mu # Gassmann 方程求干岩石体积模量 k_dry (k_sat * (phi * k_min / k_fl_old 1 - phi) - k_min) / \ (phi * k_min / k_fl_old k_sat / k_min - 1 - phi) # 新流体下的饱和模量 k_sat_new k_dry (1 - k_dry / k_min)**2 / \ (phi / k_fl_new (1 - phi) / k_min - k_dry / k_min**2) # 新速度与密度 rho_new rho_sat phi * (rho_fl_new - rho_fl_old) vp_new np.sqrt((k_sat_new 4/3 * mu) / rho_new) vs_new np.sqrt(mu / rho_new) return vp_new, vs_new, rho_new这里k_min是矿物模量砂岩常用石英的 36 GPak_fl_old和k_fl_new是替换前后流体模量水约 2.2 GPa气约 0.1 GPa。phi是孔隙度。坑在于单位——测井速度是 m/s密度是 g/cm³直接代进去会差几个数量级。我习惯在函数入口统一转成 SI出口再转回来。另一个坑是 Gassmann 假设低频对高频子波并不严格成立所以流体替换结果更适合做定性趋势分析别拿它当定量标定。4. 正演流程里的避坑清单五个让我返工的血泪经验4.1 角度单位不统一导致曲线全错现象算出来的反射系数曲线在 0 度附近就不对和理论值对不上。原因入射角用度传入函数内部直接调np.sin把度当弧度算了。解决在函数入口强制theta np.deg2rad(theta)或者约定所有接口只用弧度并在文档里写死。这个坑我踩过两次第二次是因为换了一个人写调用脚本忘了约定。4.2 子波主频与采样率不匹配产生假频现象合成道集上出现规则的高频振荡像干扰波。原因子波主频设成 60 Hz但采样率只有 2 ms奈奎斯特频率才 250 Hz理论上够但子波长度取了 5 个周期离散化后主瓣采样点太少褶积时产生假频。解决保证每个周期至少 8 到 10 个采样点即dt 1/(10*freq)。60 Hz 子波dt至少 1 ms最好 0.5 ms。4.3 抽层过细导致层间多次波假象现象合成道集上出现很多规则的低幅同相轴像是多次波。原因抽层时把测井曲线的每个小波动都当成一层层数上百层间反射叠加出假多次波。解决设定最小层厚阈值比如 5 m小于它的层合并到相邻层。同时检查速度曲线如果层间速度差小于 5%也可以合并。4.4 密度曲线缺失时用 Gardner 公式的陷阱现象密度曲线缺失用 Gardner 公式rho 0.23 * vp^0.25补结果 AVO 梯度完全不对。原因Gardner 公式是统计关系对含气砂岩偏差极大而含气正是 AVO 最关心的场景。解决如果研究目标涉及流体密度必须实测或从其他可靠来源获取别用 Gardner 硬补。实在没有至少分岩性用不同的 Gardner 系数并做敏感性分析。4.5 边界处理不当让有限差分结果全是边界反射现象全波场模拟的道集上边界处出现强能量盖住了有效信号。原因有限差分没有加吸收边界波碰到模型边缘全反射回来。解决用 PML 或海绵边界PML 层数取 10 到 20 个网格衰减系数要调太小吸不干净太大数值不稳定。我一般先用海绵边界快速验证正式跑再换 PML。5. 进阶技巧用批量正演做 AVO 敏感性分析单次正演只能看一个模型真正有价值的是批量跑参数组合看哪些参数对 AVO 响应最敏感。我常用的做法是固定岩石物理背景对孔隙度、含水饱和度、泥质含量各取若干水平交叉组合批量做流体替换和正演最后统计不同角度段的反射系数变化。import itertools import numpy as np def batch_forward(phi_list, sw_list, vsh_list, base_params): 批量正演返回每个参数组合的 AVO 曲线 results [] for phi, sw, vsh in itertools.product(phi_list, sw_list, vsh_list): # 根据 phi, sw, vsh 计算新的 vp, vs, rho # 这里调用前面的 gassmann 和岩石物理模型 vp, vs, rho rock_physics_model(phi, sw, vsh, base_params) # 计算反射系数 rc zoeppritz(base_params[vp1], base_params[vs1], base_params[rho1], vp, vs, rho, np.deg2rad(np.arange(0, 41, 1))) results.append({phi: phi, sw: sw, vsh: vsh, rc: rc}) return resultsitertools.product做笛卡尔积参数水平别取太多每个 3 到 4 个水平组合数就上百了再乘角度数计算量不小。跑完后我一般画两种图一种是固定两个参数看第三个参数变化时 AVO 曲线的分异程度另一种是算不同角度段的反射系数差值做敏感性排序。经验上孔隙度和含水饱和度对 AVO 梯度影响最大泥质含量主要影响截距。这个结论不是绝对的取决于你的岩石物理模型和背景趋势所以批量正演的意义就在于用你自己的数据验证一遍。验证正演结果是否可靠最直接的办法是拿实际井旁道做对比。把正演合成的角道集和实际道集并排显示看主要同相轴的时间、相位、振幅随角度的变化趋势是否一致。如果趋势对不上先检查抽层和子波再检查流体替换参数。别指望完全吻合实际数据有噪声、有衰减、有各向异性正演是简化模型能对上趋势就算成功。我自己在这个方向上最大的教训是别急着写代码先把模型和参数在纸上理清楚。我见过太多人包括我自己一上来就调库跑结果卡在单位、角度、采样率这些低级问题上浪费好几天。后来我养成了一个习惯任何正演任务先手算一个两层模型的解析解和代码结果对一遍对上了再往下做。这个“后悔药”帮我省了无数次返工。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
