企业官网建设流程全解析

简介：本资源是一个面向地球物理勘探与地震资料处理初学者的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 正演例程，至少应该包含：

如果你的 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)/(Z2+Z1)，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 是有偏的。

所以例程包里如果同时有精确解和近似解，我建议你在主程序里同时计算两条曲线，并输出相对误差或直接叠图。常规的界限是：

这个“先看角度范围再选公式”的习惯，能帮你避开很多解读阶段的假象。正演里算错的反射系数，到了反演阶段就是地震资料上的假亮点。

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 度到多少度）
• 输出方式（曲线、道集、交会图）

以最常见的含气砂岩模型为例，参数大致是这样的：

注意含气砂岩的 Vp 明显比含水砂岩低，但 Vs 几乎不变，这就是气层导致的“纵波速度下降、横波速度基本不变”的经典流体响应，也是 AVO 能够识别流体的底层逻辑。如果你在例程里把这些参数替换进去，直接就能看到含水砂岩顶面的反射振幅随角度变化缓慢，而含气砂岩顶面的反射振幅随角度明显变负。

3.3 运行与验证：得到的道集合理吗

主脚本运行后，你通常会看到两类图：一类是反射系数曲线 R(θ)，另一类是合成的 AVO 道集。道集怎么看？横轴是入射角或偏移距，纵轴是时间或深度，颜色代表振幅。在某个反射界面上，如果振幅从左到右（小角度到大角度）越来越“亮”或越来越“暗”，说明这个界面的 AVO 响应强烈。

跑完第一步，先做三件验证工作：

1. 零角度处的反射系数，用手算一下波阻抗差，确认和曲线起点一致
2. 看大角度方向的曲线是否出现异常跳动，如果有，考虑临界角效应（后面专门讲）
3. 把精确解和 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 Williams（1989）提出的含气砂岩 AVO 分类框架。这个分类虽然老，但现在工业界解释叠前道集时仍然天天在用：

用上面那组含气砂岩参数算出来的是典型的第 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 越界导致的复数传播。现在我的习惯是每次修改参数后，都固定输出一组与解析解对比的验证数值，一旦结果偏离预期，立刻回溯是参数问题还是代码问题，而不是埋头在图上找原因。希望这篇拆解能帮你省掉那些我已经替你踩过的坑。

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