★ FEATURED ARTICLE 简介本资源面向地质勘探、地球物理及数值计算方向的学习者与研究人员围绕地震岩石物理基本理论提供利用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.30K_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.22M_sat 可能变化 5% 以上K_fl 从 2.2 变到 0.1M_sat 可能下降 30%。这意味着在流体替换研究中流体模量的准确性比骨架模量的准确性更关键。5.3 一个实用检查习惯每次算完一组参数我会把干燥和饱和的 Vp、Vs、ν 列成一张小表然后问自己三个问题Vs 在饱和前后是否相等ν_sat 是否大于 ν_dryVp_sat/Vp_dry 是否在 1.05 到 1.30 之间如果三个答案都是肯定的结果基本可信如果有任何一个不对先回去查输入参数不要急着调模型。这个习惯帮我省了很多返工时间希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
