简介本资源是一个面向地球物理勘探与地震资料处理初学者的MATLAB入门级AVO正演建模工具包聚焦于振幅随偏移距变化AVO理论的编程实现与可视化验证。压缩包为RAR格式仅含1个核心文件——avoMODING.m脚本大小925B结构精简但功能完整涵盖AVO参数输入、Shuey或Aki-Richards等经典正解模型调用、角度域振幅计算及AVO响应曲线绘制等关键环节。已有135人学习下载适用于高校地质工程/地球物理学专业课程实践、科研入门训练或地震解释方法自学。读者可直接运行该脚本通过修改速度模型、密度、泊松比等岩石物理参数实时观察不同岩性组合下的AVO响应特征快速建立理论公式与实际地震表现之间的映射关系并为后续流体识别与储层预测打下编程与建模基础。 你有没有过这种经历从某个网盘或者U盘里扒下一个名为 avoMODING.rar 的压缩包解压出来一坨 MATLAB 的 .m 文件文件名倒是挺规整——zoeppritz_avo.m、shuey_approx.m、fluid_replace.m——但当你双击运行主脚本屏幕上要么报出一串矩阵维度的红色错误要么画出来的图和你在地震教科书的PPT里看到的那张AVO道集完全对不上。这个压缩包我在实验室帮人排过好几次了今天干脆把它彻底讲明白AVO正演模拟到底在模拟什么、里面的 MATLAB 例程每行代码在干嘛、以及你拿到这种来路不明的 rar 之后应该怎么最快跑出第一张可用的道集。这篇文章适合三类人刚接手叠前道集解释的勘探地球物理方向学生、需要在项目中快速搭一套AVO正演验证流程的工程师、以及单纯想搞懂“一个反射系数是怎么随入射角变化的”的MATLAB使用者。我会从物理原理讲到代码实现再讲到排查经验和扩展思路保证你读完能自己动手改参数、画出有意义的结果。1. avoMODING到底是干什么的一次AVO正演模拟的完整需求拆解1.1 压缩包背后的物理问题AVO 全称 Amplitude Variation with Offset中文一般叫“振幅随偏移距变化”或“振幅随入射角变化”。它要回答的问题非常直接当一束地震波以不同角度打到地下某个岩性分界面上时反射回来的能量大小会不会变如果会变变化的方式和岩层里的流体油、气、水有什么关系这个问题的工程背景是传统的地震剖面只能看到反射界面的“亮点”或“暗点”但亮点不一定是油气可能是煤层、火成岩或者钙质夹层。而 AVO 引入了一个额外的维度——入射角。你可以把它想象成用不同角度的灯光照同一个物体如果这个物体是哑光的正面照和斜着照亮度差别不会太大如果它是镜面的稍微换个角度反射光强度就剧烈变化。地下岩层里的流体种类恰恰会改变反射系数对入射角的“敏感度”。avoMODING 这个 rar 包里的 MATLAB 例程核心任务就是把这种“反射系数随入射角变化”的曲线、道集、交会图给算出来。它是叠前地震解释的最前端工具后面接的 AVO 属性分析、流体因子反演、弹性波阻抗反演全都建立在这个正演模拟的基础上。1.2 一个例程至少应该包含哪些模块拿到一个 avoMODING.rar我建议你先别急着运行打开文件夹看看它有没有这几类文件。一个像样的 AVO 正演例程至少应该包含模块对应文件常见命名作用精确反射系数计算zoeppritz_avo.m / solve_zoeppritz.m用 Zoeppritz 方程求解四个反射/透射系数近似公式计算shuey_approx.m / aki_richards_approx.m用线性近似公式快速计算 R(θ)用于对比和属性分析模型参数设置model_parameters.m / define_model.m定义上下层的 Vp、Vs、密度和入射角范围道集生成与绘图plot_avo_gather.m / wiggle_trace.m把反射系数显示成道集或曲线流体替换fluid_replace.m / gassmann.m利用 Gassmann 方程在含水、含油、含气之间切换看 AVO 响应差异如果你的 rar 里只有前两个文件那多半是个阉割版建议自己补一个参数设置脚本和绘图函数不然没法直观看到结果。如果文件特别多而且互相乱调用也别慌先用matlab的依赖分析工具或者手动grep一下函数名理清调用关系。我在实际折腾这个例程包的时候发现一个规律大部分人拿它跑不出结果不是代码本身的问题而是他根本不知道“正演的前提是先定义模型”。AVO 正演不是从地震数据里提取什么东西而是先假设“地下有一个含气砂岩它的 Vp、Vs、密度是这样”然后基于弹性波动理论算出这个模型应该产生什么样的反射振幅。所以参数的合理性直接决定结果的可用性。2. Zoeppritz方程和它的三个近似avomod的数学骨架2.1 精确解到底怎么求Zoeppritz 方程是 1919 年提出的它基于界面两侧位移连续和应力连续的边界条件联立四个方程同时求解入射纵波在界面上产生的反射纵波PP、反射横波PS、透射纵波TP、透射横波TS四个振幅系数。在 MATLAB 里实现这个方程核心就是一个 4×4 矩阵的求解问题。我贴一段在例程包里最常见的实现方式单位统一用 m/s 和 kg/m³入射角用度内部转弧度function [Rpp, Rps, Tpp, Tps] zoeppritz_avo(vp1, vs1, rho1, vp2, vs2, rho2, theta1) th1 theta1 * pi / 180; p sin(th1) / vp1; % 射线参数Snell 定理 th2 asin(p * vp2); % 透射纵波角 ph1 asin(p * vs1); % 反射横波角 ph2 asin(p * vs2); % 透射横波角 M [ sin(th1) cos(ph1) -sin(th2) cos(ph2) cos(th1) -sin(ph1) cos(th2) sin(ph2) sin(2*th1) (vp1/vs1)*cos(2*ph1) (rho2*vp2*vs2)/(rho1*vp1*vs1)*sin(2*th2) -(rho2*vp2*vs2)/(rho1*vp1*vs1)*cos(2*ph2) cos(2*ph1) -(vs1/vp1)*sin(2*ph1) -(rho2*vp2)/(rho1*vp1)*cos(2*ph2) -(rho2*vs2)/(rho1*vp1)*sin(2*ph2) ]; B [ -sin(th1) cos(th1) sin(2*th1) -cos(2*ph1) ]; X M \ B; Rpp X(1); Rps X(2); Tpp X(3); Tps X(4); end这里最容易踩的坑是不同教材对 Zoeppritz 矩阵的符号约定不一样有的把应力的正方向定义成朝下有的把位移分量取正方向定义成朝上导致最终结果看起来差一个负号。所以写完矩阵先别急着往下接先用垂直入射θ0验证此时反射系数应该约等于 (Z2-Z1)/(Z2Z1)ZρVp 是波阻抗。如果对不上优先检查第三行、第四行的符号而不是去改入射角。2.2 Shuey近似为什么是实际项目里最常用的Zoeppritz 的精确解虽然理论完整但公式复杂物理直觉差。1985 年 Shuey 在 Aki-Richards 线性近似的基础上把反射系数改写成关于入射角的显式表达式R(θ) R0 G·sin²θ K·(tan²θ - sin²θ)其中R0 是法向入射反射系数也叫 AVO 截距Intercept反映垂直入射时的振幅强度G 是 AVO 梯度Gradient控制振幅随入射角变化的速度是整个 AVO 分析里最核心的属性K 与纵波速度相对变化率有关在入射角小于 30 度时第三项贡献很小通常省略例程包里的 shuey_approx.m 实现通常长这样function R shuey_approx(vp1, vs1, rho1, vp2, vs2, rho2, theta) th theta * pi / 180; dvp vp2 - vp1; drho rho2 - rho1; vp (vp1 vp2) / 2; rho (rho1 rho2) / 2; vs (vs1 vs2) / 2; R0 0.5 * (dvp/vp drho/rho); G R0 - (dvp/vp) * 4*(vs/vp)^2 - (drho/rho) * 2*(vs/vp)^2; K 0.5 * dvp/vp; R R0 G * sin(th).^2 K * (tan(th).^2 - sin(th).^2); end别看这公式简单它把复杂的弹性波传播问题压缩成了三个参数和两个三角函数项直接让后续的截距-梯度分析成为可能。实际解释的流程是把实际地震道集上每个反射界面的振幅随角度的变化趋势拟合出来得到截距 P 和梯度 G然后看 P×G 的异常。含气砂岩的 P×G 通常会出现明显负异常而含水砂岩虽然有负的 P但 G 不会显著变负。2.3 误差边界什么时候不能再用两项近似我见过不少人把 Shuey 近似当万能公式用入射角都采到 45 度了还在拿两项近似做拟合结果梯度 G 被严重污染。实测下来当入射角超过 30 度以后公式里的第三项 K·(tan²θ - sin²θ) 的贡献会迅速增大如果你只取前两项拟合出来的 R0 和 G 是有偏的。所以例程包里如果同时有精确解和近似解我建议你在主程序里同时计算两条曲线并输出相对误差或直接叠图。常规的界限是最大入射角推荐方法小于 20 度两项 Shuey 近似完全够用20~30 度三项 Shuey 近似注意密度项精度大于 30 度优先用 Zoeppritz 精确解Shuey 只用于趋势分析这个“先看角度范围再选公式”的习惯能帮你避开很多解读阶段的假象。正演里算错的反射系数到了反演阶段就是地震资料上的假亮点。3. 跑通例程的完整流程从解压到画出第一张AVO道集3.1 解压后第一件事检查文件结构与依赖关系把 avoMODING.rar 解压到本地之后我强烈建议第一件事不是双击运行而是把文件夹放到一个纯英文路径下比如D:\codes\avoModing\。Windows 下 MATLAB 对中文路径的支持时好时坏尤其是当你后面要调用 MEX 文件、第三方工具箱或者写入文件时中文路径会带来一堆莫名其妙的报错。这个习惯花十秒钟就能养成但它能帮你省下一整晚的排查时间。然后打开 MATLAB用cd切到该目录运行depfun(main_avo.m) % 查看主脚本依赖的所有函数或者直接在编辑器里打开主脚本逐个点一下函数名看能否跳转到对应文件。如果发现有函数名标红找不到定义优先检查是不是子文件夹没有加进路径。用addpath(genpath(pwd))一次性把当前目录及所有子目录加进搜索路径是解决“函数未定义”最粗暴也最有效的办法。3.2 主程序参数表哪些参数必须提前想清楚跑正演之前先把这个模型的“地质身份”定好。一个 AVO 正演模型至少需要四组参数上覆泥岩的 Vp、Vs、密度下伏砂岩的 Vp、Vs、密度入射角范围从 0 度到多少度输出方式曲线、道集、交会图以最常见的含气砂岩模型为例参数大致是这样的参数上覆泥岩含水砂岩含气砂岩Vpm/s280030002600Vsm/s120015001500ρkg/m³235023502050Vp/Vs 比2.332.001.73注意含气砂岩的 Vp 明显比含水砂岩低但 Vs 几乎不变这就是气层导致的“纵波速度下降、横波速度基本不变”的经典流体响应也是 AVO 能够识别流体的底层逻辑。如果你在例程里把这些参数替换进去直接就能看到含水砂岩顶面的反射振幅随角度变化缓慢而含气砂岩顶面的反射振幅随角度明显变负。3.3 运行与验证得到的道集合理吗主脚本运行后你通常会看到两类图一类是反射系数曲线 R(θ)另一类是合成的 AVO 道集。道集怎么看横轴是入射角或偏移距纵轴是时间或深度颜色代表振幅。在某个反射界面上如果振幅从左到右小角度到大角度越来越“亮”或越来越“暗”说明这个界面的 AVO 响应强烈。跑完第一步先做三件验证工作零角度处的反射系数用手算一下波阻抗差确认和曲线起点一致看大角度方向的曲线是否出现异常跳动如果有考虑临界角效应后面专门讲把精确解和 Shuey 近似的曲线叠在一起看偏差是否在可接受范围内我自己的习惯是直接在命令行里打几个关键值对比一下[R0_zoe] zoeppritz_avo(2800,1200,2350,2600,1500,2050,0); [R0_shu] shuey_approx(2800,1200,2350,2600,1500,2050,0); fprintf(Zoeppritz R0 %.4f, Shuey R0 %.4f\n, R0_zoe, R0_shu);如果这两个值差超过 0.005说明某个函数的参数顺序或者符号约定有问题先修这个再往下走。4. 结果解读4类AVO异常和截距-梯度交会图4.1 含气砂岩在道集上长什么样跑出第一张 AVO 道集之后最想知道的当然是这个结果到底能不能说明地下含气这里需要引入 Rutherford and Williams1989提出的含气砂岩 AVO 分类框架。这个分类虽然老但现在工业界解释叠前道集时仍然天天在用类型含气砂岩阻抗法向入射反射系数振幅随角度变化特征1类高阻抗比泥岩硬正值振幅先减后增可能出现极性反转2类近零阻抗接近零反射很弱极性反转常见3类低阻抗比泥岩软负值振幅绝对值随角度增大4类低阻抗更特殊负值振幅绝对值随角度减小用上面那组含气砂岩参数算出来的是典型的第 3 类法向反射系数为负并且随入射角增大振幅的绝对值越来越大。对应的图形特征是道集上这个反射轴的“亮度”从左到右越来越强而且是负极性先负后正或先黑后白取决于显示约定。为什么第 3 类最常见因为绝大多数浅层、中深层含气砂岩都比围岩泥岩更“软”——纵波速度低、密度低导致阻抗差本来就很大再加上泊松比降低横波速度差异相对小于是角度项进一步把负振幅拉大。你如果看到自己的正演结果居然在 20 度以后振幅往回缩那要看是不是参数里给出了异常的 Vs 值或者密度压得太低。4.2 从正演到AVO属性P-G交会图怎么用例程包如果够完整里面多半还有一个函数用来拟合法向入射截距 P 和梯度 G。做法很简单对反射系数序列做最小二乘拟合theta_deg 0:0.5:30; Rpp zoeppritz_avo(vp1,vs1,rho1,vp2,vs2,rho2,theta_deg); A [ones(length(theta_deg),1), sin(theta_deg*pi/180).^2]; coef A \ Rpp(:); P coef(1); % 截距 G coef(2); % 梯度得到 P 和 G 之后把不同模型含水、含油、含气的正演结果放到同一个 P-G 交会图里你会看到它们分布在不同的象限或区域。典型含气砂岩的 P×G 为正的负值区域第三象限或沿着负 P 负 G 方向含水砂岩则更靠近坐标原点或正向区域。这个交会图是 AVO 解释里最有名的“甜点探测器”正演的意义就在于你知道一个真实气藏对应的 P、G 应该在哪个位置再看实际数据的 P、G 点是否落进来。5. 实际跑代码时最容易翻车的三个地方5.1 DLL初始化失败可能是路径和运行库的问题很多人在 MATLAB 里调用外部代码或 MEX 文件时会撞见类似这样的报错OSError: [WinError 1114] 动态链接库(DLL)初始化例程失败。Error loading D:...\xxx.dll这个错误我见到太多次了它在 Windows MATLAB 环境下高发原因通常不是代码逻辑而是系统层面的 DLL 加载问题。最常见的诱因有三个路径里有中文或空格导致 DLL 依赖的本地资源找不到目标 DLL 依赖的 Visual C 运行库缺失需要装 vc_redist.x64.exe杀毒软件把 DLL 隔离或拦截了加载时初始化函数无法执行排查建议按顺序来先把整个工程目录挪到D:\codes\这种纯英文路径再确认 MATLAB 的位数matlab -arch和你调用的 DLL 位数一致最后用Dependencies之类的工具打开 DLL看缺失的依赖项。不要一上来就怀疑 MATLAB 安装坏了大多数 WinError 1114 都是环境问题。5.2 矩阵维度报错与复数结果临界角没有处理Zoeppritz 求解里asin(p * vp2)可能算出复数因为射线参数 p sin(θ1)/vp1 是固定值当入射角增大到一定程度时p * vp2 1导致反正弦函数的定义域越界。这个入射角就是临界角。超过临界角后透射波会变成非均匀波折射回介质内部反射系数在临界角附近会出现剧烈的振幅变化。如果你不加处理直接把这个复数结果拿去画道集图里就会出现一撮“毛刺”或者 NaN 空洞。解决办法是在循环里检查abs(p * vp2)超过 1 就做截断或直接丢弃该角度同时在道集绘制时限制最大显示角度。经验法则是最大入射角取临界角的 80% 左右既能保证信息量又不会让临界角附近的噪声干扰注意力。5.3 符号约定不统一先和解析解比对再往下走我前面提到过 Zoeppritz 矩阵符号乱的问题这里再展开。不同代码库、不同论文里对“位移正方向”“反射系数极性”的定义经常不一致。最典型的例子是同一个地质模型你在 A 例程里算出的 3 类 AVO 道集是“负黑正白”在 B 例程里可能完全反过来。如果你拿自己的结果和别人的图对比发现极性反了先别怀疑地质参数先检查是不是符号约定不同。怎么快速自查用最简单的地质界面上覆是高速高密度下伏是低速低密度计算垂直入射反射系数。正常约定下它应该是负值即反射波与入射波相位相反。如果你的代码算出正号那么整个道集的颜色约定就要整体取反或者你在绘图时有意反转了极性。把这个验证脚本写进例程的头部注释里能救很多人的命。6. 从单道正演到合成道集扩展例程的进阶思路6.1 用Ricker子波做褶积生成合成地震道只画反射系数曲线在正演层面虽然够用但地震解释人员看的是“地震道集”也就是反射系数经过子波褶积后的结果。扩展例程很自然的下一步就是把 Ricker 子波和反射系数做褶积生成更接近真实地震记录的道集。t 0:0.001:0.8; w ricker_wavelet(30, 0.001); % 30Hz Ricker 子波需要自备或自己写 r_trace zeros(size(t)); r_trace(200) Rpp(1); % 假设一个界面在 0.2s for i 2:length(theta_deg) r_trace_i zeros(size(t)); r_trace_i(200) Rpp(i); synth(:,i) conv(w, r_trace_i, same); end这里最关键的是“道集上每个角度的子波波形要保持一致”否则你观察到的振幅变化可能只是子波旁瓣的干涉结果而不是真实的 AVO 响应。实际地震道集在近偏移距和远偏移距上的子波会因为动校正拉伸而有差异这是另一个处理环节的问题正演阶段可以暂时忽略但心里要有数。6.2 从正演走向反演AVO属性提取与流体识别正演例程跑熟之后你会自然想到一个应用如果我从合成道集或实际道集里拟合出 P 和 G能不能反推下伏岩层的弹性参数这一步就是从正演到反演的桥梁。常用的套路是利用 P、G 组合出流体因子比如流体因子 F P G某些物性条件下与含气饱和度相关性好泊松比变化率 Δσ 的近似公式λρ、μρ 弹性参数反演我在实际项目里比较常用的是把 P-G 交会图和流体替换结果结合先用 Gassmann 方程把同一个砂岩分别替换成含水、含油、含气三种状态正演出三组 P-G 点再把实测数据的 P-G 点投影到图上看落在哪个流体附近。这样一来正演就不再是单纯画几条曲线而是直接参与储层流体判别的决策链。Gassmann 流体替换的简化实现并不复杂核心是把岩石骨架的体积模量从含水状态换算到目标流体状态再重新算纵波速度K_sat1 rho1 * (vp1^2 - 4/3 * vs1^2); K_sat2 ...; % 带入目标流体参数 vp2 sqrt((K_sat2 4/3 * vs2^2) / rho2);注意这里的单位要统一密度用 kg/m³速度用 m/s模量单位就是 Pa。我踩过一次坑密度用 g/cm³、速度用 km/s算出来的 K 小了 10 的 6 次方倍所有速度更新全错。建议在脚本开头强制做单位转换把所有参数统一成国际单位制再计算。最后再分享一个我在实际使用中的体会拿到 avoMODING 这类例程包别急着贪多求全先把 Zoeppritz 精确解、Shuey 近似、单界面道集这三样东西跑明白比下载十个扩展包都管用。我刚接触 AVO 那会儿曾在临界角处理上栽过跟头画出过一条“振幅先增后减又暴增”的道集后来发现就是 p×vs2 越界导致的复数传播。现在我的习惯是每次修改参数后都固定输出一组与解析解对比的验证数值一旦结果偏离预期立刻回溯是参数问题还是代码问题而不是埋头在图上找原因。希望这篇拆解能帮你省掉那些我已经替你踩过的坑。本文还有配套的精品资源点击获取 AVO正演模拟入门：Zoeppritz方程与MATLAB实现全解析
简介本资源是一个面向地球物理勘探与地震资料处理初学者的MATLAB入门级AVO正演建模工具包聚焦于振幅随偏移距变化AVO理论的编程实现与可视化验证。压缩包为RAR格式仅含1个核心文件——avoMODING.m脚本大小925B结构精简但功能完整涵盖AVO参数输入、Shuey或Aki-Richards等经典正解模型调用、角度域振幅计算及AVO响应曲线绘制等关键环节。已有135人学习下载适用于高校地质工程/地球物理学专业课程实践、科研入门训练或地震解释方法自学。读者可直接运行该脚本通过修改速度模型、密度、泊松比等岩石物理参数实时观察不同岩性组合下的AVO响应特征快速建立理论公式与实际地震表现之间的映射关系并为后续流体识别与储层预测打下编程与建模基础。 你有没有过这种经历从某个网盘或者U盘里扒下一个名为 avoMODING.rar 的压缩包解压出来一坨 MATLAB 的 .m 文件文件名倒是挺规整——zoeppritz_avo.m、shuey_approx.m、fluid_replace.m——但当你双击运行主脚本屏幕上要么报出一串矩阵维度的红色错误要么画出来的图和你在地震教科书的PPT里看到的那张AVO道集完全对不上。这个压缩包我在实验室帮人排过好几次了今天干脆把它彻底讲明白AVO正演模拟到底在模拟什么、里面的 MATLAB 例程每行代码在干嘛、以及你拿到这种来路不明的 rar 之后应该怎么最快跑出第一张可用的道集。这篇文章适合三类人刚接手叠前道集解释的勘探地球物理方向学生、需要在项目中快速搭一套AVO正演验证流程的工程师、以及单纯想搞懂“一个反射系数是怎么随入射角变化的”的MATLAB使用者。我会从物理原理讲到代码实现再讲到排查经验和扩展思路保证你读完能自己动手改参数、画出有意义的结果。1. avoMODING到底是干什么的一次AVO正演模拟的完整需求拆解1.1 压缩包背后的物理问题AVO 全称 Amplitude Variation with Offset中文一般叫“振幅随偏移距变化”或“振幅随入射角变化”。它要回答的问题非常直接当一束地震波以不同角度打到地下某个岩性分界面上时反射回来的能量大小会不会变如果会变变化的方式和岩层里的流体油、气、水有什么关系这个问题的工程背景是传统的地震剖面只能看到反射界面的“亮点”或“暗点”但亮点不一定是油气可能是煤层、火成岩或者钙质夹层。而 AVO 引入了一个额外的维度——入射角。你可以把它想象成用不同角度的灯光照同一个物体如果这个物体是哑光的正面照和斜着照亮度差别不会太大如果它是镜面的稍微换个角度反射光强度就剧烈变化。地下岩层里的流体种类恰恰会改变反射系数对入射角的“敏感度”。avoMODING 这个 rar 包里的 MATLAB 例程核心任务就是把这种“反射系数随入射角变化”的曲线、道集、交会图给算出来。它是叠前地震解释的最前端工具后面接的 AVO 属性分析、流体因子反演、弹性波阻抗反演全都建立在这个正演模拟的基础上。1.2 一个例程至少应该包含哪些模块拿到一个 avoMODING.rar我建议你先别急着运行打开文件夹看看它有没有这几类文件。一个像样的 AVO 正演例程至少应该包含模块对应文件常见命名作用精确反射系数计算zoeppritz_avo.m / solve_zoeppritz.m用 Zoeppritz 方程求解四个反射/透射系数近似公式计算shuey_approx.m / aki_richards_approx.m用线性近似公式快速计算 R(θ)用于对比和属性分析模型参数设置model_parameters.m / define_model.m定义上下层的 Vp、Vs、密度和入射角范围道集生成与绘图plot_avo_gather.m / wiggle_trace.m把反射系数显示成道集或曲线流体替换fluid_replace.m / gassmann.m利用 Gassmann 方程在含水、含油、含气之间切换看 AVO 响应差异如果你的 rar 里只有前两个文件那多半是个阉割版建议自己补一个参数设置脚本和绘图函数不然没法直观看到结果。如果文件特别多而且互相乱调用也别慌先用matlab的依赖分析工具或者手动grep一下函数名理清调用关系。我在实际折腾这个例程包的时候发现一个规律大部分人拿它跑不出结果不是代码本身的问题而是他根本不知道“正演的前提是先定义模型”。AVO 正演不是从地震数据里提取什么东西而是先假设“地下有一个含气砂岩它的 Vp、Vs、密度是这样”然后基于弹性波动理论算出这个模型应该产生什么样的反射振幅。所以参数的合理性直接决定结果的可用性。2. Zoeppritz方程和它的三个近似avomod的数学骨架2.1 精确解到底怎么求Zoeppritz 方程是 1919 年提出的它基于界面两侧位移连续和应力连续的边界条件联立四个方程同时求解入射纵波在界面上产生的反射纵波PP、反射横波PS、透射纵波TP、透射横波TS四个振幅系数。在 MATLAB 里实现这个方程核心就是一个 4×4 矩阵的求解问题。我贴一段在例程包里最常见的实现方式单位统一用 m/s 和 kg/m³入射角用度内部转弧度function [Rpp, Rps, Tpp, Tps] zoeppritz_avo(vp1, vs1, rho1, vp2, vs2, rho2, theta1) th1 theta1 * pi / 180; p sin(th1) / vp1; % 射线参数Snell 定理 th2 asin(p * vp2); % 透射纵波角 ph1 asin(p * vs1); % 反射横波角 ph2 asin(p * vs2); % 透射横波角 M [ sin(th1) cos(ph1) -sin(th2) cos(ph2) cos(th1) -sin(ph1) cos(th2) sin(ph2) sin(2*th1) (vp1/vs1)*cos(2*ph1) (rho2*vp2*vs2)/(rho1*vp1*vs1)*sin(2*th2) -(rho2*vp2*vs2)/(rho1*vp1*vs1)*cos(2*ph2) cos(2*ph1) -(vs1/vp1)*sin(2*ph1) -(rho2*vp2)/(rho1*vp1)*cos(2*ph2) -(rho2*vs2)/(rho1*vp1)*sin(2*ph2) ]; B [ -sin(th1) cos(th1) sin(2*th1) -cos(2*ph1) ]; X M \ B; Rpp X(1); Rps X(2); Tpp X(3); Tps X(4); end这里最容易踩的坑是不同教材对 Zoeppritz 矩阵的符号约定不一样有的把应力的正方向定义成朝下有的把位移分量取正方向定义成朝上导致最终结果看起来差一个负号。所以写完矩阵先别急着往下接先用垂直入射θ0验证此时反射系数应该约等于 (Z2-Z1)/(Z2Z1)ZρVp 是波阻抗。如果对不上优先检查第三行、第四行的符号而不是去改入射角。2.2 Shuey近似为什么是实际项目里最常用的Zoeppritz 的精确解虽然理论完整但公式复杂物理直觉差。1985 年 Shuey 在 Aki-Richards 线性近似的基础上把反射系数改写成关于入射角的显式表达式R(θ) R0 G·sin²θ K·(tan²θ - sin²θ)其中R0 是法向入射反射系数也叫 AVO 截距Intercept反映垂直入射时的振幅强度G 是 AVO 梯度Gradient控制振幅随入射角变化的速度是整个 AVO 分析里最核心的属性K 与纵波速度相对变化率有关在入射角小于 30 度时第三项贡献很小通常省略例程包里的 shuey_approx.m 实现通常长这样function R shuey_approx(vp1, vs1, rho1, vp2, vs2, rho2, theta) th theta * pi / 180; dvp vp2 - vp1; drho rho2 - rho1; vp (vp1 vp2) / 2; rho (rho1 rho2) / 2; vs (vs1 vs2) / 2; R0 0.5 * (dvp/vp drho/rho); G R0 - (dvp/vp) * 4*(vs/vp)^2 - (drho/rho) * 2*(vs/vp)^2; K 0.5 * dvp/vp; R R0 G * sin(th).^2 K * (tan(th).^2 - sin(th).^2); end别看这公式简单它把复杂的弹性波传播问题压缩成了三个参数和两个三角函数项直接让后续的截距-梯度分析成为可能。实际解释的流程是把实际地震道集上每个反射界面的振幅随角度的变化趋势拟合出来得到截距 P 和梯度 G然后看 P×G 的异常。含气砂岩的 P×G 通常会出现明显负异常而含水砂岩虽然有负的 P但 G 不会显著变负。2.3 误差边界什么时候不能再用两项近似我见过不少人把 Shuey 近似当万能公式用入射角都采到 45 度了还在拿两项近似做拟合结果梯度 G 被严重污染。实测下来当入射角超过 30 度以后公式里的第三项 K·(tan²θ - sin²θ) 的贡献会迅速增大如果你只取前两项拟合出来的 R0 和 G 是有偏的。所以例程包里如果同时有精确解和近似解我建议你在主程序里同时计算两条曲线并输出相对误差或直接叠图。常规的界限是最大入射角推荐方法小于 20 度两项 Shuey 近似完全够用20~30 度三项 Shuey 近似注意密度项精度大于 30 度优先用 Zoeppritz 精确解Shuey 只用于趋势分析这个“先看角度范围再选公式”的习惯能帮你避开很多解读阶段的假象。正演里算错的反射系数到了反演阶段就是地震资料上的假亮点。3. 跑通例程的完整流程从解压到画出第一张AVO道集3.1 解压后第一件事检查文件结构与依赖关系把 avoMODING.rar 解压到本地之后我强烈建议第一件事不是双击运行而是把文件夹放到一个纯英文路径下比如D:\codes\avoModing\。Windows 下 MATLAB 对中文路径的支持时好时坏尤其是当你后面要调用 MEX 文件、第三方工具箱或者写入文件时中文路径会带来一堆莫名其妙的报错。这个习惯花十秒钟就能养成但它能帮你省下一整晚的排查时间。然后打开 MATLAB用cd切到该目录运行depfun(main_avo.m) % 查看主脚本依赖的所有函数或者直接在编辑器里打开主脚本逐个点一下函数名看能否跳转到对应文件。如果发现有函数名标红找不到定义优先检查是不是子文件夹没有加进路径。用addpath(genpath(pwd))一次性把当前目录及所有子目录加进搜索路径是解决“函数未定义”最粗暴也最有效的办法。3.2 主程序参数表哪些参数必须提前想清楚跑正演之前先把这个模型的“地质身份”定好。一个 AVO 正演模型至少需要四组参数上覆泥岩的 Vp、Vs、密度下伏砂岩的 Vp、Vs、密度入射角范围从 0 度到多少度输出方式曲线、道集、交会图以最常见的含气砂岩模型为例参数大致是这样的参数上覆泥岩含水砂岩含气砂岩Vpm/s280030002600Vsm/s120015001500ρkg/m³235023502050Vp/Vs 比2.332.001.73注意含气砂岩的 Vp 明显比含水砂岩低但 Vs 几乎不变这就是气层导致的“纵波速度下降、横波速度基本不变”的经典流体响应也是 AVO 能够识别流体的底层逻辑。如果你在例程里把这些参数替换进去直接就能看到含水砂岩顶面的反射振幅随角度变化缓慢而含气砂岩顶面的反射振幅随角度明显变负。3.3 运行与验证得到的道集合理吗主脚本运行后你通常会看到两类图一类是反射系数曲线 R(θ)另一类是合成的 AVO 道集。道集怎么看横轴是入射角或偏移距纵轴是时间或深度颜色代表振幅。在某个反射界面上如果振幅从左到右小角度到大角度越来越“亮”或越来越“暗”说明这个界面的 AVO 响应强烈。跑完第一步先做三件验证工作零角度处的反射系数用手算一下波阻抗差确认和曲线起点一致看大角度方向的曲线是否出现异常跳动如果有考虑临界角效应后面专门讲把精确解和 Shuey 近似的曲线叠在一起看偏差是否在可接受范围内我自己的习惯是直接在命令行里打几个关键值对比一下[R0_zoe] zoeppritz_avo(2800,1200,2350,2600,1500,2050,0); [R0_shu] shuey_approx(2800,1200,2350,2600,1500,2050,0); fprintf(Zoeppritz R0 %.4f, Shuey R0 %.4f\n, R0_zoe, R0_shu);如果这两个值差超过 0.005说明某个函数的参数顺序或者符号约定有问题先修这个再往下走。4. 结果解读4类AVO异常和截距-梯度交会图4.1 含气砂岩在道集上长什么样跑出第一张 AVO 道集之后最想知道的当然是这个结果到底能不能说明地下含气这里需要引入 Rutherford and Williams1989提出的含气砂岩 AVO 分类框架。这个分类虽然老但现在工业界解释叠前道集时仍然天天在用类型含气砂岩阻抗法向入射反射系数振幅随角度变化特征1类高阻抗比泥岩硬正值振幅先减后增可能出现极性反转2类近零阻抗接近零反射很弱极性反转常见3类低阻抗比泥岩软负值振幅绝对值随角度增大4类低阻抗更特殊负值振幅绝对值随角度减小用上面那组含气砂岩参数算出来的是典型的第 3 类法向反射系数为负并且随入射角增大振幅的绝对值越来越大。对应的图形特征是道集上这个反射轴的“亮度”从左到右越来越强而且是负极性先负后正或先黑后白取决于显示约定。为什么第 3 类最常见因为绝大多数浅层、中深层含气砂岩都比围岩泥岩更“软”——纵波速度低、密度低导致阻抗差本来就很大再加上泊松比降低横波速度差异相对小于是角度项进一步把负振幅拉大。你如果看到自己的正演结果居然在 20 度以后振幅往回缩那要看是不是参数里给出了异常的 Vs 值或者密度压得太低。4.2 从正演到AVO属性P-G交会图怎么用例程包如果够完整里面多半还有一个函数用来拟合法向入射截距 P 和梯度 G。做法很简单对反射系数序列做最小二乘拟合theta_deg 0:0.5:30; Rpp zoeppritz_avo(vp1,vs1,rho1,vp2,vs2,rho2,theta_deg); A [ones(length(theta_deg),1), sin(theta_deg*pi/180).^2]; coef A \ Rpp(:); P coef(1); % 截距 G coef(2); % 梯度得到 P 和 G 之后把不同模型含水、含油、含气的正演结果放到同一个 P-G 交会图里你会看到它们分布在不同的象限或区域。典型含气砂岩的 P×G 为正的负值区域第三象限或沿着负 P 负 G 方向含水砂岩则更靠近坐标原点或正向区域。这个交会图是 AVO 解释里最有名的“甜点探测器”正演的意义就在于你知道一个真实气藏对应的 P、G 应该在哪个位置再看实际数据的 P、G 点是否落进来。5. 实际跑代码时最容易翻车的三个地方5.1 DLL初始化失败可能是路径和运行库的问题很多人在 MATLAB 里调用外部代码或 MEX 文件时会撞见类似这样的报错OSError: [WinError 1114] 动态链接库(DLL)初始化例程失败。Error loading D:...\xxx.dll这个错误我见到太多次了它在 Windows MATLAB 环境下高发原因通常不是代码逻辑而是系统层面的 DLL 加载问题。最常见的诱因有三个路径里有中文或空格导致 DLL 依赖的本地资源找不到目标 DLL 依赖的 Visual C 运行库缺失需要装 vc_redist.x64.exe杀毒软件把 DLL 隔离或拦截了加载时初始化函数无法执行排查建议按顺序来先把整个工程目录挪到D:\codes\这种纯英文路径再确认 MATLAB 的位数matlab -arch和你调用的 DLL 位数一致最后用Dependencies之类的工具打开 DLL看缺失的依赖项。不要一上来就怀疑 MATLAB 安装坏了大多数 WinError 1114 都是环境问题。5.2 矩阵维度报错与复数结果临界角没有处理Zoeppritz 求解里asin(p * vp2)可能算出复数因为射线参数 p sin(θ1)/vp1 是固定值当入射角增大到一定程度时p * vp2 1导致反正弦函数的定义域越界。这个入射角就是临界角。超过临界角后透射波会变成非均匀波折射回介质内部反射系数在临界角附近会出现剧烈的振幅变化。如果你不加处理直接把这个复数结果拿去画道集图里就会出现一撮“毛刺”或者 NaN 空洞。解决办法是在循环里检查abs(p * vp2)超过 1 就做截断或直接丢弃该角度同时在道集绘制时限制最大显示角度。经验法则是最大入射角取临界角的 80% 左右既能保证信息量又不会让临界角附近的噪声干扰注意力。5.3 符号约定不统一先和解析解比对再往下走我前面提到过 Zoeppritz 矩阵符号乱的问题这里再展开。不同代码库、不同论文里对“位移正方向”“反射系数极性”的定义经常不一致。最典型的例子是同一个地质模型你在 A 例程里算出的 3 类 AVO 道集是“负黑正白”在 B 例程里可能完全反过来。如果你拿自己的结果和别人的图对比发现极性反了先别怀疑地质参数先检查是不是符号约定不同。怎么快速自查用最简单的地质界面上覆是高速高密度下伏是低速低密度计算垂直入射反射系数。正常约定下它应该是负值即反射波与入射波相位相反。如果你的代码算出正号那么整个道集的颜色约定就要整体取反或者你在绘图时有意反转了极性。把这个验证脚本写进例程的头部注释里能救很多人的命。6. 从单道正演到合成道集扩展例程的进阶思路6.1 用Ricker子波做褶积生成合成地震道只画反射系数曲线在正演层面虽然够用但地震解释人员看的是“地震道集”也就是反射系数经过子波褶积后的结果。扩展例程很自然的下一步就是把 Ricker 子波和反射系数做褶积生成更接近真实地震记录的道集。t 0:0.001:0.8; w ricker_wavelet(30, 0.001); % 30Hz Ricker 子波需要自备或自己写 r_trace zeros(size(t)); r_trace(200) Rpp(1); % 假设一个界面在 0.2s for i 2:length(theta_deg) r_trace_i zeros(size(t)); r_trace_i(200) Rpp(i); synth(:,i) conv(w, r_trace_i, same); end这里最关键的是“道集上每个角度的子波波形要保持一致”否则你观察到的振幅变化可能只是子波旁瓣的干涉结果而不是真实的 AVO 响应。实际地震道集在近偏移距和远偏移距上的子波会因为动校正拉伸而有差异这是另一个处理环节的问题正演阶段可以暂时忽略但心里要有数。6.2 从正演走向反演AVO属性提取与流体识别正演例程跑熟之后你会自然想到一个应用如果我从合成道集或实际道集里拟合出 P 和 G能不能反推下伏岩层的弹性参数这一步就是从正演到反演的桥梁。常用的套路是利用 P、G 组合出流体因子比如流体因子 F P G某些物性条件下与含气饱和度相关性好泊松比变化率 Δσ 的近似公式λρ、μρ 弹性参数反演我在实际项目里比较常用的是把 P-G 交会图和流体替换结果结合先用 Gassmann 方程把同一个砂岩分别替换成含水、含油、含气三种状态正演出三组 P-G 点再把实测数据的 P-G 点投影到图上看落在哪个流体附近。这样一来正演就不再是单纯画几条曲线而是直接参与储层流体判别的决策链。Gassmann 流体替换的简化实现并不复杂核心是把岩石骨架的体积模量从含水状态换算到目标流体状态再重新算纵波速度K_sat1 rho1 * (vp1^2 - 4/3 * vs1^2); K_sat2 ...; % 带入目标流体参数 vp2 sqrt((K_sat2 4/3 * vs2^2) / rho2);注意这里的单位要统一密度用 kg/m³速度用 m/s模量单位就是 Pa。我踩过一次坑密度用 g/cm³、速度用 km/s算出来的 K 小了 10 的 6 次方倍所有速度更新全错。建议在脚本开头强制做单位转换把所有参数统一成国际单位制再计算。最后再分享一个我在实际使用中的体会拿到 avoMODING 这类例程包别急着贪多求全先把 Zoeppritz 精确解、Shuey 近似、单界面道集这三样东西跑明白比下载十个扩展包都管用。我刚接触 AVO 那会儿曾在临界角处理上栽过跟头画出过一条“振幅先增后减又暴增”的道集后来发现就是 p×vs2 越界导致的复数传播。现在我的习惯是每次修改参数后都固定输出一组与解析解对比的验证数值一旦结果偏离预期立刻回溯是参数问题还是代码问题而不是埋头在图上找原因。希望这篇拆解能帮你省掉那些我已经替你踩过的坑。本文还有配套的精品资源点击获取