企业官网建设流程全解析

简介：本资源面向地质勘探、地球物理及数值计算方向的学习者与研究人员，围绕地震岩石物理基本理论，提供利用MATLAB算法计算干燥岩石与流体饱和岩石弹性参数的实现方案，重点求解纵波模量、剪切模量和泊松比，并涉及云计算环境下的大规模数值模拟思路。压缩包内共1个文件，为m脚本文件，整体约2KB，体积轻量，便于直接阅读与二次修改，适合作为岩石物理建模的入门参考或课程实验素材。资源从弹性参数定义出发，串联体积模量、剪切模量与泊松比对地震波传播速度和衰减的影响，帮助读者理解地壳物理特性与储层评估之间的关联。目前已有545人学习下载，具备一定参考热度。通过该脚本，读者可掌握岩石弹性参数的计算流程与MATLAB实现方式，并据此拓展到流体替换、横波速度预测及地下储层特性分析等场景，为后续科研与工程应用提供可复用的代码基础。

1. 干燥与饱和：为什么同一块岩石测出两套弹性参数

同一块砂岩岩心，干燥状态下测得的纵波速度可能是 3800 m/s，饱和水之后却跳到 4200 m/s 以上，泊松比也从 0.15 附近抬到 0.30 左右。这不是仪器漂移，而是孔隙流体改变了岩石的有效弹性响应。地震岩石物理基本理论要解决的核心问题之一，就是把这套变化用可计算的模型描述出来：给定矿物模量、孔隙度、流体类型，算出干燥岩石和流体饱和岩石的纵波模量、剪切模量和泊松比。

这套计算在实操中主要服务三类人：做测井解释的需要用饱和岩石模量标定横波速度；做地震正演的需要干燥骨架模量作为输入；做储层预测的需要对比饱和与干燥差异来识别流体。Gassmann 方程是连接干燥和饱和两端的桥梁，但它有明确的适用边界，用错了比不用更危险。下面从理论到代码，把这条链路走通。

2. 从矿物模量到干燥骨架：Gassmann 方程的输入到底怎么准备

2.1 四个必须分清的模量概念

初学者最容易翻车的地方，是把矿物模量、骨架模量、干燥岩石模量和饱和岩石模量混为一谈。这四个量在 Gassmann 方程里各有位置，代错了结果会偏得离谱。

矿物模量（K_min、μ_min）描述的是构成岩石的固体颗粒本身的弹性性质，石英的 K 约 36.6 GPa、μ 约 45 GPa，方解石 K 约 76.8 GPa、μ 约 32 GPa。骨架模量（K_frame、μ_frame）描述的是矿物颗粒在干燥状态下通过接触和胶结形成的固体框架的弹性响应，它一定小于矿物模量，因为框架里有孔隙。干燥岩石模量在多数语境下和骨架模量同义，指孔隙中只有空气且空气模量近似为零时的岩石等效模量。饱和岩石模量则是孔隙被流体填充后的等效模量。

剪切模量在 Gassmann 框架下有一个重要结论：流体饱和不改变岩石的剪切模量。因为理想流体不能承受剪切应力，孔隙流体的存在不影响岩石对剪切变形的抵抗能力。所以 μ_sat = μ_dry，这个等式在后续计算泊松比时会直接用到。

2.2 用 Hashin-Shtrikman 下界估算骨架模量

实际工作中，干燥骨架模量往往没有直接测量值，需要用矿物模量和孔隙度估算。常见做法是 Hashin-Shtrikman 下界，它给出的骨架模量是给定孔隙度下理论上的下限值，对于孔隙形状比较均匀的碎屑岩比较适用。

import numpy as np def hsw_lower_bound(K_min, mu_min, K_inc, mu_inc, phi): """ Hashin-Shtrikman 下界，用于估算干燥骨架模量。 K_min, mu_min: 矿物模量 (GPa) K_inc, mu_inc: 孔隙相模量 (GPa)，干燥时取 0 phi: 孔隙度 (小数) """ # 干燥条件下孔隙相模量为零，直接代入 K_hs = 1.0 / ((1 - phi) / (K_min + 4/3 * mu_min) + phi / (K_inc + 4/3 * mu_min)) - 4/3 * mu_min mu_hs = 1.0 / ((1 - phi) / (mu_min + mu_min * (9*K_min + 8*mu_min) / (6*(K_min + 2*mu_min))) + phi / (mu_inc + mu_min * (9*K_min + 8*mu_min) / (6*(K_min + 2*mu_min)))) - \ mu_min * (9*K_min + 8*mu_min) / (6*(K_min + 2*mu_min)) return K_hs, mu_hs # 石英为主的高孔隙砂岩示例 K_min, mu_min = 36.6, 45.0 phi = 0.25 K_dry, mu_dry = hsw_lower_bound(K_min, mu_min, 0.0, 0.0, phi) print(f"干燥骨架模量: K_dry={K_dry:.2f} GPa, mu_dry={mu_dry:.2f} GPa")

这段代码的关键参数是 K_inc 和 mu_inc，干燥时取零。如果孔隙中已经有流体，这两个值要换成流体的体积模量和剪切模量。孔隙度 phi 直接决定骨架模量的下降幅度，phi 从 0.20 变到 0.30，K_dry 可能下降 30% 以上。实际使用时要注意，HS 下界假设孔隙形状比较规则，对于裂缝型储层会低估骨架模量，这时候需要换用自洽模型或差分有效介质模型。

2.3 干燥岩石的纵波模量和泊松比

拿到 K_dry 和 μ_dry 之后，干燥岩石的纵波模量 M_dry 和泊松比 ν_dry 可以直接算：

def dry_elastic_params(K_dry, mu_dry): """计算干燥岩石的纵波模量 M 和泊松比 nu""" M_dry = K_dry + 4/3 * mu_dry nu_dry = (3*K_dry - 2*mu_dry) / (2*(3*K_dry + mu_dry)) return M_dry, nu_dry M_dry, nu_dry = dry_elastic_params(K_dry, mu_dry) print(f"干燥纵波模量 M_dry={M_dry:.2f} GPa, 泊松比 nu_dry={nu_dry:.3f}")

纵波模量 M 有时也叫 P-wave modulus，它和纵波速度的关系是 Vp = sqrt(M/ρ)。泊松比的物理意义是横向应变与纵向应变之比的负值，干燥砂岩通常在 0.10 到 0.20 之间，比饱和砂岩明显偏低。这个差异正是流体识别的基础。

3. 流体饱和后的模量：Gassmann 方程怎么代、参数怎么选

3.1 Gassmann 方程的标准形式与代入顺序

Gassmann 方程描述的是低频条件下，孔隙流体完全连通且与固体框架达到压力平衡时的饱和岩石模量。标准形式如下：

K_sat = K_dry + (1 - K_dry/K_min)^2 / (phi/K_fl + (1-phi)/K_min - K_dry/K_min^2)

其中 K_fl 是孔隙流体的有效体积模量。这个公式的代入顺序很重要：先算分母，再算分子，最后加到 K_dry 上。分母如果接近零，说明输入的参数组合在物理上不合理，计算结果会发散。

def gassmann_saturated(K_dry, mu_dry, K_min, K_fl, phi): """ Gassmann 方程计算饱和岩石体积模量。 K_fl: 流体体积模量 (GPa)，水约 2.2，油约 1.0，气约 0.1 """ numerator = (1 - K_dry / K_min) ** 2 denominator = phi / K_fl + (1 - phi) / K_min - K_dry / K_min ** 2 if denominator <= 0: raise ValueError("Gassmann 分母非正，检查输入参数是否物理合理") K_sat = K_dry + numerator / denominator mu_sat = mu_dry # 剪切模量不变 return K_sat, mu_sat K_sat, mu_sat = gassmann_saturated(K_dry, mu_dry, K_min, 2.2, phi) print(f"饱和体积模量 K_sat={K_sat:.2f} GPa, 剪切模量 mu_sat={mu_sat:.2f} GPa")

K_fl 的选择直接决定饱和模量的高低。水的体积模量在常温常压下约 2.2 GPa，但在地层条件下要考虑温度压力校正，高温下可能降到 1.8 GPa 左右。如果孔隙中是油气混合，需要用 Wood 公式先算混合流体的等效模量：1/K_fl = S_oil/K_oil + S_gas/K_gas + S_water/K_water。这一步在流体替换研究中是必须做的，直接用水模量代替混合流体模量会高估饱和模量。

3.2 饱和岩石的纵波模量和泊松比

拿到 K_sat 和 μ_sat 之后，计算方式和干燥岩石完全一致：

M_sat = K_sat + 4/3 * mu_sat nu_sat = (3*K_sat - 2*mu_sat) / (2*(3*K_sat + mu_sat)) print(f"饱和纵波模量 M_sat={M_sat:.2f} GPa, 泊松比 nu_sat={nu_sat:.3f}")

对比干燥和饱和的结果，通常会看到 M_sat 明显大于 M_dry，而 μ 不变，ν_sat 大于 ν_dry。这个规律在砂岩中非常稳定，如果算出来 ν_sat 反而小于 ν_dry，几乎可以肯定是输入参数出了问题，优先检查 K_fl 是否代错或者 phi 是否超出了合理范围。

3.3 从模量到速度：验证计算结果是否合理

模量本身看不见摸不着，最终要转成速度才能和测井或实验数据对比：

def moduli_to_velocity(M, mu, rho): """将纵波模量和剪切模量转换为纵波、横波速度 (m/s)""" Vp = np.sqrt(M * 1e9 / rho) Vs = np.sqrt(mu * 1e9 / rho) return Vp, Vs # 假设干燥和饱和密度 rho_dry = 2.65 * (1 - phi) * 1000 # kg/m3 rho_sat = (2.65 * (1 - phi) + 1.0 * phi) * 1000 Vp_dry, Vs_dry = moduli_to_velocity(M_dry, mu_dry, rho_dry) Vp_sat, Vs_sat = moduli_to_velocity(M_sat, mu_sat, rho_sat) print(f"干燥: Vp={Vp_dry:.0f} m/s, Vs={Vs_dry:.0f} m/s") print(f"饱和: Vp={Vp_sat:.0f} m/s, Vs={Vs_sat:.0f} m/s")

注意单位换算：模量用 GPa，密度用 kg/m³，速度出来是 m/s。干燥和饱和的密度不同，饱和密度要加上孔隙流体的质量。如果算出来的 Vs 在饱和前后完全一样，说明剪切模量确实没变，这是正确的；如果 Vp_sat 比 Vp_dry 高出 20% 以上，对于高孔隙砂岩是正常的，但如果是低孔隙致密岩石，这个增幅应该很小。

4. 避坑与排查：Gassmann 计算中最容易翻车的五个地方

4.1 现象：饱和模量比干燥模量还小

原因：K_fl 代入了负值或者分母计算时符号搞反。Gassmann 方程要求 K_dry < K_min 且 K_fl > 0，如果 K_dry 被高估到接近甚至超过 K_min，分子会变成负数。

解决：先检查 K_dry 是否小于 K_min。对于孔隙度 25% 的砂岩，K_dry 通常在 10 到 20 GPa 之间，如果算出来 30 GPa 以上，说明骨架模量估算有问题，回到 HS 下界重新检查孔隙度输入。

4.2 现象：泊松比算出来大于 0.5 或者小于 0

原因：K 或 μ 出现了负值，或者单位不统一。泊松比的物理上限是 0.5（对应不可压缩材料），超过 0.5 说明 K 为负或者 μ 为负。

解决：在每一步计算后加断言检查。K_dry > 0、μ_dry > 0、K_sat > 0、μ_sat > 0，这四个条件必须同时满足。如果 μ_dry 为负，通常是 HS 下界公式中孔隙相模量代入错误，干燥条件下 mu_inc 应该取 0 而不是取矿物剪切模量。

4.3 现象：高频实验数据对不上 Gassmann 预测

原因：Gassmann 方程是低频极限，适用条件是波在孔隙流体中传播时压力有足够时间平衡。超声波实验频率在 MHz 量级，流体来不及平衡，会出现频散，实测速度高于 Gassmann 预测。

解决：如果对比的是超声波实验数据，不要直接用 Gassmann 结果做定量对比，只能做趋势判断。要做定量对比，需要用 Biot 理论或喷射流模型考虑频散效应。测井频率在 kHz 量级，地震在 Hz 量级，这两个频段 Gassmann 基本适用。

4.4 现象：含气砂岩的饱和模量算出来比水饱和还高

原因：气的体积模量约 0.1 GPa，远小于水的 2.2 GPa，含气砂岩的 K_sat 应该低于含水砂岩。如果算反了，通常是 K_fl 在含气时代入了水的值，或者混合流体模量计算时饱和度权重搞反。

解决：用 Wood 公式算混合流体模量时，1/K_fl = Σ(S_i/K_i)，含气饱和度越高，K_fl 越小。检查饱和度之和是否为 1，含气时 K_fl 应该在 0.1 到 1.0 GPa 之间，含水时在 2.0 到 2.5 GPa 之间。

4.5 现象：同一组参数换台机器结果不一样

原因：浮点数精度差异或者单位换算系数不同。GPa 到 Pa 是 1e9，密度用 g/cm³ 还是 kg/m³ 差 1000 倍，这些换算在代码里如果写死容易出错。

解决：把所有单位换算集中在一个函数里，输入输出都标注清楚单位。模量统一用 GPa，密度统一用 kg/m³，速度统一用 m/s。在代码开头定义常量 GPa_TO_PA = 1e9，所有换算引用这个常量，不要手写数字。

5. 进阶技巧：用流体替换做定量解释和敏感性分析

5.1 从饱和反推干燥：流体替换的逆向计算

实际测井得到的是饱和岩石模量，但做正演时需要干燥骨架模量。Gassmann 方程可以反解：

def gassmann_dry_from_sat(K_sat, mu_sat, K_min, K_fl, phi): """从饱和模量反推干燥骨架模量""" # 先算干燥体积模量 a = K_sat / K_min b = phi / K_fl + (1 - phi) / K_min # 解二次方程 K_dry = (K_sat * (b - 1) + K_min) / (a * (b - 1) + 1) # 更稳健的做法是用迭代，但解析解在多数情况下够用 mu_dry = mu_sat return K_dry, mu_dry

这个反推在流体替换中很常用：先反推干燥骨架，再代入新的流体模量正算新的饱和模量，就完成了从一种流体饱和状态到另一种的替换。比如把含水砂岩替换成含气砂岩，看速度变化量，这个变化量就是地震流体识别的理论依据。

5.2 参数敏感性：哪个输入对结果影响最大

做敏感性分析时，固定其他参数，逐个扰动 10%，看输出变化多少。根据我的经验，对饱和纵波模量影响从大到小的顺序通常是：孔隙度 > 流体模量 > 矿物模量 > 干燥骨架模量。孔隙度从 0.20 变到 0.22，M_sat 可能变化 5% 以上；K_fl 从 2.2 变到 0.1，M_sat 可能下降 30%。这意味着在流体替换研究中，流体模量的准确性比骨架模量的准确性更关键。

5.3 一个实用检查习惯

每次算完一组参数，我会把干燥和饱和的 Vp、Vs、ν 列成一张小表，然后问自己三个问题：Vs 在饱和前后是否相等？ν_sat 是否大于 ν_dry？Vp_sat/Vp_dry 是否在 1.05 到 1.30 之间？如果三个答案都是肯定的，结果基本可信；如果有任何一个不对，先回去查输入参数，不要急着调模型。这个习惯帮我省了很多返工时间，希望帮到你。

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