简介：这份资源面向地质勘探、地球物理与油藏工程方向的研究人员和学生，聚焦岩石物理性质与Gassman公式的MATLAB实现。包内共6个m文件，压缩包约3KB，涵盖Gassman公式计算、Reuss下限与Voigt上限估算、Whyliamson平均弹性参数方法、Biot-Gassman扩展算法以及流体敏感性参数K的测量函数，构成一套轻量但完整的岩石物理分析工具箱。已有366人学习下载，说明其在含流体岩石弹性响应模拟中具有一定参考价值。读者可借助这些脚本理解孔隙流体对岩石弹性模量的影响，掌握多相材料弹性界限的估算思路，并将其应用于地层评价与石油勘探场景，适合具备一定MATLAB基础、希望深入Gassman理论与岩石物理建模的进阶学习者。

1. 从一组 MATLAB 函数说起：岩石物理里那套 Gassman 公式到底怎么落地

做地震正演或者测井解释的人，迟早会撞上同一个问题：同一块砂岩，含水饱和度和含气饱和度下的纵波速度能差出百分之二三十，可岩芯测试又不可能把每种流体都替一遍。这时候绕不开的就是 Gassman 公式——它把干燥骨架的弹性模量、孔隙流体的体积模量、孔隙度这几个量串起来，让你能从一组已知状态推算另一组。但公式写在纸上是一回事，真拿 MATLAB 敲出来又是另一回事：Voigt 上限和 Reuss 下限到底该在什么场景下用、Biot 系数和 Gassman 公式怎么耦合、KMeasure 里的 K 参数物理意义是什么，这些细节没人给你拆，光看公式推不出来。这份源码包就是干这个的——六个 .m 文件，把从上下界估算到 Gassman 流体替换的链路全串起来了，适合做岩石物理正演、流体替换或者地震响应模拟的工程师直接拿去改。

2. 六个 .m 文件的分工：先搞清楚每个函数在算什么

2.1 Voigt 上限与 Reuss 下限：为什么不能只用其中一个

多相介质的等效弹性模量估算，最朴素的做法就是 Voigt 和 Reuss 两个极端。Voigt 假设各相应变一致，给出的是体积模量和剪切模量的上限；Reuss 假设各相应力一致，给出的是下限。voigt_upper_bounds.m和reuss_lower_bounds.m这两个文件干的就是这件事。

实际用的时候有个血泪经验：很多人只拿 Voigt 上限当等效模量往 Gassman 公式里代，结果算出来的速度比实测高出一大截。原因很简单——Voigt 是等应变假设，对应的是各相平行排列的理想情况，真实岩石的孔隙流体和骨架根本不是这种几何关系。Reuss 下限虽然偏低，但它反映的是流体和骨架串联的极端。真实值永远落在两者之间，所以这两个函数通常是一起用的，用来框定一个合理区间。

看一个典型的调用方式：

% 假设有两相：石英骨架和孔隙流体 K1 = 36e9; % 石英体积模量，单位 Pa K2 = 2.2e9; % 水的体积模量 G1 = 44e9; % 石英剪切模量 G2 = 0; % 流体剪切模量为零 f1 = 0.85; % 石英体积分数 f2 = 0.15; % 流体体积分数 % Voigt 上限 [K_voigt, G_voigt] = voigt_upper_bounds(K1, K2, G1, G2, f1, f2); % Reuss 下限 [K_reuss, G_reuss] = reuss_lower_bounds(K1, K2, G1, G2, f1, f2); % 输出对比 fprintf('Voigt 上限: K = %.2f GPa, G = %.2f GPa\n', K_voigt/1e9, G_voigt/1e9); fprintf('Reuss 下限: K = %.2f GPa, G = %.2f GPa\n', K_reuss/1e9, G_reuss/1e9);

这里的关键参数是体积分数f1和f2，两者之和必须为 1。如果传进去的分数不归一化，算出来的模量会直接偏掉。另外注意流体剪切模量G2通常设为零，因为液体和气体不承受剪切应力，这个值如果误设成非零，Reuss 下限会算出偏高的剪切模量，后面 Gassman 替换就全乱了。

2.2 whyli.m：上下界之间的那根“锚”

光有上下界还不够，实际做流体替换时需要的是一个确定的等效模量值，而不是一个区间。whyli.m这个文件名对应的应该是 Williamson 方法——一种在 Voigt 和 Reuss 之间做加权平均的经典方案。它的核心思路是引入一个权重因子，根据孔隙几何或者经验关系来调节上下界的贡献比例。

我一般会这么用：

% 在 Voigt 和 Reuss 之间做加权 % w 是权重，0 对应纯 Reuss，1 对应纯 Voigt w = 0.5; % 初始猜测，通常根据岩芯标定调整 K_eff = whyli(K_voigt, K_reuss, w); G_eff = whyli(G_voigt, G_reuss, w); fprintf('等效体积模量: %.2f GPa\n', K_eff/1e9); fprintf('等效剪切模量: %.2f GPa\n', G_eff/1e9);

权重w是这里最玄学的参数。没有岩芯标定数据的时候，我一般从 0.5 开始试，然后拿已知流体饱和状态的测井速度反推——如果含水饱和度的纵波速度对不上，就调w直到匹配。有岩芯数据的话，直接拿干燥和饱和两种状态的模量反算w，比拍脑袋靠谱得多。

2.3 GassmanFormula.m 与 Boit_Gassman.m：流体替换的核心

这两个文件是整个包的核心。GassmanFormula.m实现的是经典 Gassman 方程：

[ K_{sat} = K_{dry} + \frac{(1 - K_{dry}/K_{min})^2}{\phi/K_f + (1-\phi)/K_{min} - K_{dry}/K_{min}^2} ]

其中 (K_{sat}) 是饱和岩石体积模量，(K_{dry}) 是干燥骨架体积模量，(K_{min}) 是矿物基质体积模量，(K_f) 是流体体积模量，(\phi) 是孔隙度。剪切模量在 Gassman 公式里不受流体影响，所以 (G_{sat} = G_{dry})。

Boit_Gassman.m则是在 Gassman 基础上引入了 Biot 系数 (\alpha = 1 - K_{dry}/K_{min})，用来描述孔隙压力变化对有效应力的贡献。Biot 系数在流体替换和孔隙压力预测里非常关键，尤其是做 4D 地震或者压力监测的时候，光用 Gassman 公式不够，得把 Biot 系数单独拎出来算。

调用逻辑大致是这样：

% 已知干燥骨架模量和矿物模量 K_dry = 15e9; % 干燥骨架体积模量 G_dry = 12e9; % 干燥骨架剪切模量 K_min = 36e9; % 石英矿物模量 phi = 0.2; % 孔隙度 % 第一种流体：水 K_water = 2.2e9; [K_sat_water, G_sat_water] = GassmanFormula(K_dry, G_dry, K_min, K_water, phi); % 第二种流体：气 K_gas = 0.1e9; [K_sat_gas, G_sat_gas] = GassmanFormula(K_dry, G_dry, K_min, K_gas, phi); % 计算 Biot 系数 alpha = Boit_Gassman(K_dry, K_min); fprintf('含水饱和: K = %.2f GPa, G = %.2f GPa\n', K_sat_water/1e9, G_sat_water/1e9); fprintf('含气饱和: K = %.2f GPa, G = %.2f GPa\n', K_sat_gas/1e9, G_sat_gas/1e9); fprintf('Biot 系数: %.3f\n', alpha);

参数说明：K_dry和G_dry通常来自干燥岩芯的超声波测试或者经验关系估算；K_min根据岩性选取，石英约 36 GPa，方解石约 70 GPa；phi是孔隙度，小数形式。注意K_f对气体非常敏感，0.1 GPa 和 0.01 GPa 算出来的饱和模量能差不少，所以气体体积模量最好用状态方程算，别随便估。

2.4 KMeasure.m：那个容易被忽略的 K 参数

KMeasure.m从文件名看是计算流体敏感性参数 K 的。在 Gassman 框架里，K 参数描述的是岩石骨架对孔隙流体压力变化的响应程度，本质上和 Biot 系数是同一类东西，但有些文献里会单独定义。这个函数在实际流体替换流程里通常用来做质量控制——如果算出来的 K 值超出物理合理范围（比如大于 1 或者小于 0），说明输入的模量或者孔隙度有问题。

我一般会在跑完 Gassman 替换后加一步验证：

% 计算 K 参数做质控 K_param = KMeasure(K_dry, K_min, phi); if K_param < 0 || K_param > 1 warning('K 参数超出合理范围: %.3f，检查输入模量和孔隙度', K_param); end

这个检查看着简单，但能拦住不少低级错误——比如把矿物模量和骨架模量搞反了，或者孔隙度输成了百分数而不是小数。

3. 从干燥骨架到饱和岩石：完整流体替换流程与参数标定

3.1 数据准备：你需要哪些输入量

跑通这套代码，最少需要以下几类数据：

如果只有测井速度而没有模量，得先做反演：

% 从纵波速度 Vp、横波速度 Vs 和密度 rho 反演模量 Vp = 3500; % m/s Vs = 2000; % m/s rho = 2300; % kg/m3 G = rho * Vs^2; % 剪切模量 K = rho * (Vp^2 - 4/3*Vs^2); % 体积模量 fprintf('反演体积模量: %.2f GPa\n', K/1e9); fprintf('反演剪切模量: %.2f GPa\n', G/1e9);

这一步的坑在于：测井速度是饱和状态的，反演出来的是饱和模量，要先用 Gassman 公式反推干燥模量，才能做后续的流体替换。反推的时候把 Gassman 公式里的 K_sat 当已知，K_dry 当未知，解一个二次方程就行。

3.2 分步执行：从干燥到饱和的完整链路

假设你已经有了干燥骨架模量，现在要算含水饱和和含气饱和两种状态下的弹性参数，完整流程如下：

第一步，用 Voigt 和 Reuss 估算混合矿物的等效模量。如果岩石是纯石英，跳过这步；如果是石英加黏土，必须做。

% 石英和黏土两相混合 K_quartz = 36e9; G_quartz = 44e9; K_clay = 21e9; G_clay = 7e9; f_quartz = 0.8; f_clay = 0.2; [K_min_voigt, G_min_voigt] = voigt_upper_bounds(K_quartz, K_clay, G_quartz, G_clay, f_quartz, f_clay); [K_min_reuss, G_min_reuss] = reuss_lower_bounds(K_quartz, K_clay, G_quartz, G_clay, f_quartz, f_clay); % 用 whyli 取加权平均 w_min = 0.5; K_min = whyli(K_min_voigt, K_min_reuss, w_min); G_min = whyli(G_min_voigt, G_min_reuss, w_min);

第二步，用 Gassman 公式做流体替换。注意这里输入的 K_dry 和 G_dry 是干燥骨架的，不是矿物的。

% 干燥骨架参数 K_dry = 15e9; G_dry = 12e9; phi = 0.2; % 水饱和 K_water = 2.2e9; [K_sat_w, G_sat_w] = GassmanFormula(K_dry, G_dry, K_min, K_water, phi); % 气饱和 K_gas = 0.1e9; [K_sat_g, G_sat_g] = GassmanFormula(K_dry, G_dry, K_min, K_gas, phi);

第三步，从模量算速度，做最终对比。

% 密度计算 rho_min = 2650; % 矿物密度 rho_water = 1000; rho_gas = 200; rho_sat_w = rho_min * (1 - phi) + rho_water * phi; rho_sat_g = rho_min * (1 - phi) + rho_gas * phi; % 速度计算 Vp_w = sqrt((K_sat_w + 4/3*G_sat_w) / rho_sat_w); Vs_w = sqrt(G_sat_w / rho_sat_w); Vp_g = sqrt((K_sat_g + 4/3*G_sat_g) / rho_sat_g); Vs_g = sqrt(G_sat_g / rho_sat_g); fprintf('水饱和: Vp = %.0f m/s, Vs = %.0f m/s\n', Vp_w, Vs_w); fprintf('气饱和: Vp = %.0f m/s, Vs = %.0f m/s\n', Vp_g, Vs_g);

这套流程跑下来，你能清楚看到流体替换对纵波速度的影响远大于横波——因为 Gassman 公式里流体只影响体积模量，不影响剪切模量。这也是为什么做流体识别的时候，Vp/Vs 比值比单纯 Vp 更管用。

3.3 参数敏感性：哪些输入最要命

跑通流程只是第一步，真正做项目的时候，输入参数的误差会直接传到结果上。我做过一轮简单的敏感性测试，结果如下：

孔隙度是最敏感的，所以测井解释的孔隙度精度直接决定流体替换的可靠性。干燥骨架模量排第二，这也是为什么岩芯测试那么重要——没有实测 K_dry，光靠经验关系估，误差可能到 20% 以上。流体模量对气体特别敏感，但好在气体模量可以通过状态方程算得比较准，只要温压数据可靠。

4. 避坑与排查：跑这套代码时最容易翻车的五个地方

4.1 现象：算出来的饱和模量比干燥模量还小

原因：Gassman 公式里 K_sat 必须大于等于 K_dry，因为流体填充孔隙会让岩石更难压缩。如果算出来反了，大概率是 K_f 输错了——比如把气体体积模量输成了 2.2e9（水的值），或者孔隙度输成了负值。

解决：先检查 K_f 是否在合理范围。水的体积模量约 2.2 GPa，油约 1.0-1.5 GPa，气约 0.01-0.1 GPa。如果 K_f 大于 K_min，公式会出问题。另外确认 phi 是小数，0.2 不是 20。

4.2 现象：Voigt 上限和 Reuss 下限算出来一样

原因：两相的模量太接近，或者体积分数输成了 0 和 1。如果 f1=1、f2=0，上下界自然重合。另一种可能是调用函数时参数顺序传错了，把 K 和 G 搞混了。

解决：打印中间变量确认。在调用上下界函数之前，先输出各相的 K、G 和体积分数，确认没有零值或者顺序错误。

4.3 现象：Biot 系数算出大于 1 或小于 0

原因：Biot 系数 α = 1 - K_dry/K_min，物理上必须在 0 到 1 之间。如果 K_dry > K_min，α 就是负的；如果 K_dry 是负的，α 就大于 1。这通常意味着干燥骨架模量估算错了——K_dry 必须小于 K_min，因为骨架包含孔隙，肯定比纯矿物好压缩。

解决：检查 K_dry 的来源。如果是反演出来的，确认反演流程没有把饱和模量当干燥模量用。如果是经验关系估的，换一个更保守的估算公式。

4.4 现象：含气饱和度的速度算出来比含水还低很多，但实测没差那么多

原因：Gassman 公式假设孔隙流体压力均匀，低频假设。实际地震频带下，部分饱和岩石存在流体斑块效应，速度会比 Gassman 预测的高。这是 Gassman 公式的固有局限，不是代码 bug。

解决：如果工区存在明显的气水混合带，考虑用斑块饱和模型做修正。简单做法是在 Gassman 结果上加一个经验校正量，或者用 Biot-Gassman 的频散版本。这套代码里的Boit_Gassman.m可能已经包含了一些修正逻辑，值得细看。

4.5 现象：whyli.m 的权重 w 怎么调都匹配不上实测

原因：Williamson 方法是经验性的，权重 w 没有物理约束。如果岩芯标定数据本身有误差，或者岩石的孔隙结构不是简单的两相混合，光调 w 可能永远对不上。

解决：先确认岩芯数据质量。如果岩芯测试的干燥和饱和速度本身就不一致，先排查测试流程。如果数据没问题，尝试分段标定——不同孔隙度区间用不同的 w，而不是全井段一个值。

5. 进阶用法：用 KMeasure 做质控和不确定性分析

跑通基本流程之后，这套代码还有一个容易被忽略的用法：拿KMeasure.m做不确定性传播分析。具体做法是给每个输入参数加一个扰动，跑蒙特卡洛模拟，看输出速度的分布范围。

% 蒙特卡洛不确定性分析 N = 1000; Vp_samples = zeros(N, 1); for i = 1:N % 给输入加随机扰动 K_dry_i = K_dry * (1 + 0.1 * randn()); phi_i = phi * (1 + 0.05 * randn()); K_f_i = K_water * (1 + 0.08 * randn()); % 跑 Gassman 替换 [K_sat_i, G_sat_i] = GassmanFormula(K_dry_i, G_dry, K_min, K_f_i, phi_i); % 算速度 rho_i = rho_min * (1 - phi_i) + rho_water * phi_i; Vp_samples(i) = sqrt((K_sat_i + 4/3*G_sat_i) / rho_i); end % 统计结果 Vp_mean = mean(Vp_samples); Vp_std = std(Vp_samples); fprintf('Vp 均值: %.0f m/s, 标准差: %.0f m/s\n', Vp_mean, Vp_std); fprintf('90%% 置信区间: [%.0f, %.0f] m/s\n', prctile(Vp_samples, 5), prctile(Vp_samples, 95));

这个分析的价值在于：你能告诉决策者，基于当前输入数据的质量，流体替换预测的速度有多大的不确定范围。如果 90% 置信区间跨过了含水饱和和含气饱和的分界线，那这个预测就不能用来做流体识别，得先补数据。

我自己的习惯是，每次做流体替换之前，先跑一遍蒙特卡洛，看看不确定范围。如果范围太宽，先回去搞数据，别急着出结论。从那以后我每次做流体替换都强制走一遍 KMeasure 质控加蒙特卡洛，拦住了不少“看起来合理但实际不可靠”的结果。希望帮到你。

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