简介这份资源面向地质勘探、地球物理及数值计算方向的学习者与研究人员围绕地震岩石物理基本理论用MATLAB实现干燥岩石与流体饱和岩石的弹性参数计算涵盖纵波模量、剪切模量和泊松比三类核心指标。压缩包内共1个文件为m脚本文件整体约2KB体积轻量便于直接阅读源码、理解公式到代码的映射关系并可在本地或云端环境中快速运行验证。资源结合云计算思路适合需要处理大规模岩体状态模拟、批量计算不同条件下弹性参数的场景帮助读者掌握从理论模型到数值实现的完整链路。目前已有545人学习下载可作为地震波传播分析、地下储层特性评估与地球资源探测相关课题的入门参考与算法练习素材。1. 地震岩石物理入门从干燥骨架到流体饱和的弹性参数计算很多刚接触地震岩石物理的同行第一次拿到测井曲线或者岩心测试报告时都会卡在同一个地方实验室报出来的是干燥岩石的纵波模量、剪切模量和泊松比但地下储层里孔隙流体是客观存在的这两套数据到底怎么换算、能不能直接对比我当年做第一个碳酸盐岩项目时就因为忽略了流体替换这一步把干燥骨架的泊松比直接套到饱和层段结果横波预测和实测对不上被审核人一眼看穿。这个标题要解决的就是这个核心问题以地震岩石物理基本理论为出发点把干燥岩石和流体饱和岩石的基本弹性参数——纵波模量、剪切模量和泊松比——算清楚、算对、算到能用于正演和反演。适合谁看做过一点测井解释、想往岩石物理建模方向走的工程师或者正在做AVO正演、需要自己构造弹性参数曲线的研究生。下面我从理论框架开始一步步把公式、参数和代码落地。2. 理论框架与参数定义为什么干燥和饱和要分开算2.1 从胡克定律到各向同性弹性参数各向同性线弹性介质中应力与应变通过广义胡克定律联系独立弹性常数只有两个。工程上常用的三件套是纵波模量 ( M )、剪切模量 ( \mu ) 和泊松比 ( u )。这里要特别注意一个容易混淆的点纵波模量 ( M ) 和体积模量 ( K ) 不是一回事。( M ) 是轴向应力与轴向应变之比侧向受限而 ( K ) 是静水压力与体积应变之比。很多文献里把 ( M ) 写成 ( K \frac{4}{3}\mu )这个关系只在各向同性条件下成立用的时候要确认前提。干燥岩石的弹性参数直接由矿物基质和孔隙结构决定流体饱和后孔隙流体的存在会改变岩石的有效弹性响应。Gassmann理论给出了低频条件下流体饱和岩石体积模量的解析表达式这是整个流体替换流程的基石。剪切模量在Gassmann理论中不随流体变化因为理想流体不能承受剪切应力这一点是后续所有计算的锚点。2.2 干燥岩石弹性参数的计算路径干燥岩石的纵波模量 ( M_{dry} )、剪切模量 ( \mu_{dry} ) 和泊松比 ( u_{dry} ) 通常来自实验室超声测量或测井数据。如果拿到的是纵波速度 ( V_p ) 和横波速度 ( V_s ) 以及密度 ( \rho )计算路径非常直接[ \mu_{dry} \rho_{dry} V_{s,dry}^2 ][ M_{dry} \rho_{dry} V_{p,dry}^2 ][ u_{dry} \frac{M_{dry} - 2\mu_{dry}}{2(M_{dry} - \mu_{dry})} ]这里 ( \rho_{dry} ) 是干燥岩石的体密度单位kg/m³速度单位m/s模量单位Pa。实际工作中我一般统一用GPa输出方便和测井曲线量纲对齐。注意泊松比的表达式是从 ( M K \frac{4}{3}\mu ) 和 ( u \frac{3K - 2\mu}{2(3K \mu)} ) 联立推出来的代入时确保 ( M \mu )否则会出现负泊松比这种在常规岩石中不常见的值。2.3 流体饱和岩石的Gassmann方程与参数映射Gassmann方程的核心思想是低频条件下孔隙流体有足够时间在孔隙空间中达到压力平衡。饱和岩石的体积模量 ( K_{sat} ) 由下式给出[ K_{sat} K_{dry} \frac{\left(1 - \frac{K_{dry}}{K_{min}}\right)^2}{\frac{\phi}{K_{fl}} \frac{1 - \phi}{K_{min}} - \frac{K_{dry}}{K_{min}^2}} ]其中 ( K_{dry} ) 是干燥岩石体积模量( K_{min} ) 是矿物基质体积模量( K_{fl} ) 是孔隙流体体积模量( \phi ) 是孔隙度。干燥岩石体积模量 ( K_{dry} ) 由 ( M_{dry} - \frac{4}{3}\mu_{dry} ) 得到。饱和岩石剪切模量 ( \mu_{sat} \mu_{dry} )这是Gassmann理论的直接结论。饱和岩石纵波模量 ( M_{sat} K_{sat} \frac{4}{3}\mu_{sat} )泊松比 ( u_{sat} \frac{3K_{sat} - 2\mu_{sat}}{2(3K_{sat} \mu_{sat})} )。整个计算链条中( K_{min} )、( K_{fl} ) 和 ( \phi ) 是三个必须准确输入的参数任何一个偏差都会传递到最终结果。3. 动手实现用Python计算干燥与饱和弹性参数3.1 输入参数准备与单位统一先明确输入数据。假设我们有一组干燥砂岩样品数据纵波速度 ( V_p 3500 ) m/s横波速度 ( V_s 2100 ) m/s干燥密度 ( \rho_{dry} 2350 ) kg/m³孔隙度 ( \phi 0.18 )矿物基质体积模量 ( K_{min} 36 ) GPa石英为主孔隙流体为盐水体积模量 ( K_{fl} 2.2 ) GPa。这些值在常规砂岩储层中具有代表性实际使用时替换成自己的测量值即可。单位统一是第一个坑。速度用m/s密度用kg/m³算出来的模量是Pa除以 ( 10^9 ) 转成GPa。Gassmann方程中所有模量必须同量纲我习惯全部转成GPa再代入避免中间步骤出现 ( 10^9 ) 量级的系数错误。3.2 干燥岩石弹性参数的计算代码import numpy as np # 输入参数 Vp_dry 3500.0 # 干燥纵波速度, m/s Vs_dry 2100.0 # 干燥横波速度, m/s rho_dry 2350.0 # 干燥密度, kg/m³ phi 0.18 # 孔隙度, 小数 K_min 36.0 # 矿物基质体积模量, GPa K_fl 2.2 # 流体体积模量, GPa # 干燥岩石模量计算 mu_dry rho_dry * Vs_dry**2 / 1e9 # 剪切模量, GPa M_dry rho_dry * Vp_dry**2 / 1e9 # 纵波模量, GPa K_dry M_dry - 4.0/3.0 * mu_dry # 体积模量, GPa # 干燥岩石泊松比 nu_dry (M_dry - 2*mu_dry) / (2*(M_dry - mu_dry)) print(f干燥岩石: mu_dry{mu_dry:.3f} GPa, M_dry{M_dry:.3f} GPa, K_dry{K_dry:.3f} GPa, nu_dry{nu_dry:.4f})这段代码的逻辑很直接从速度和密度出发先算 ( \mu_{dry} ) 和 ( M_{dry} )再反推 ( K_{dry} )。注意 ( K_{dry} ) 必须为正如果出现负值说明输入的 ( V_p/V_s ) 比值不合理通常 ( V_p/V_s ) 小于1.5时就要警惕。泊松比的计算中分母 ( M_{dry} - \mu_{dry} ) 不能为零实际岩石中这个条件自然满足。3.3 流体饱和岩石弹性参数的Gassmann替换# Gassmann流体替换 # 饱和岩石体积模量 numerator (1 - K_dry / K_min)**2 denominator phi / K_fl (1 - phi) / K_min - K_dry / K_min**2 K_sat K_dry numerator / denominator # 饱和岩石剪切模量不变 mu_sat mu_dry # 饱和岩石纵波模量 M_sat K_sat 4.0/3.0 * mu_sat # 饱和岩石泊松比 nu_sat (3*K_sat - 2*mu_sat) / (2*(3*K_sat mu_sat)) print(f饱和岩石: K_sat{K_sat:.3f} GPa, mu_sat{mu_sat:.3f} GPa, M_sat{M_sat:.3f} GPa, nu_sat{nu_sat:.4f})Gassmann方程的分母是数值稳定性的关键。当 ( \phi ) 很小且 ( K_{fl} ) 远小于 ( K_{min} ) 时( \phi/K_{fl} ) 这一项会很大导致分母整体偏大( K_{sat} ) 接近 ( K_{dry} )。这是物理上合理的低孔硬岩中流体对整体刚度的影响有限。反过来高孔软岩中 ( K_{sat} ) 会明显大于 ( K_{dry} )。如果分母出现负值或零检查 ( K_{dry} ) 是否超过了 ( K_{min} )这在物理上不可能说明输入数据有误。3.4 结果验证与中间量检查算完之后不要直接扔进正演流程先做三个检查。第一( \mu_{sat} ) 必须等于 ( \mu_{dry} )如果不等代码里哪里改了剪切模量。第二( K_{sat} ) 必须大于 ( K_{dry} )因为流体填充孔隙会提高整体体积模量。第三( u_{sat} ) 通常大于 ( u_{dry} )因为流体饱和会提高 ( V_p/V_s ) 比值。如果出现相反的趋势优先怀疑 ( K_{fl} ) 给错了比如把气体的体积模量约0.1 GPa当成水的代入。我一般会把干燥和饱和的 ( V_p )、( V_s ) 也反算出来做对比# 反算饱和岩石速度 rho_sat rho_dry phi * 1000.0 # 假设流体密度1000 kg/m³ Vp_sat np.sqrt(M_sat * 1e9 / rho_sat) Vs_sat np.sqrt(mu_sat * 1e9 / rho_sat) print(f饱和速度: Vp_sat{Vp_sat:.1f} m/s, Vs_sat{Vs_sat:.1f} m/s) print(f速度变化: dVp{Vp_sat-Vp_dry:.1f} m/s, dVs{Vs_sat-Vs_dry:.1f} m/s)饱和密度 ( \rho_{sat} \rho_{dry} \phi \rho_{fl} )这里假设流体密度为1000 kg/m³。速度变化量是判断流体替换是否合理的直观指标水饱和通常使 ( V_p ) 略微下降或基本不变( V_s ) 基本不变或略降因为密度增加抵消了模量增加的部分效果。如果 ( V_p ) 大幅上升检查 ( K_{fl} ) 是否给得过高。4. 避坑与排查干燥/饱和计算中的五个血泪教训4.1 现象泊松比算出负值原因干燥岩石的 ( V_p/V_s ) 比值低于1.45导致 ( M_{dry} 2\mu_{dry} )泊松比公式分子为负。这种情况在含气砂岩或高裂隙岩石中可能出现但更常见的是横波速度拾取错误。解决先检查 ( V_s ) 是否被低估。如果数据确实如此不要强行修正负泊松比在特定条件下是物理真实的但需要在报告中注明。我一般会画 ( V_p/V_s ) 交会图把异常点标出来单独处理。4.2 现象Gassmann替换后 ( K_{sat} ) 小于 ( K_{dry} )原因分母中 ( K_{dry}/K_{min}^2 ) 项过大通常是因为 ( K_{dry} ) 被高估或者 ( K_{min} ) 被低估。另一个可能是 ( \phi ) 给成了小数但代码里当成了百分数。解决确认 ( \phi ) 的输入格式0.18就是18%。检查 ( K_{dry} ) 是否小于 ( K_{min} )如果 ( K_{dry} K_{min} )说明干燥岩石比矿物基质还硬物理上不可能。回到速度数据检查 ( V_p ) 是否异常高。4.3 现象饱和岩石 ( V_p ) 大幅高于干燥岩石原因( K_{fl} ) 给成了水的值但实际孔隙流体是气水混合或者 ( \phi ) 给得过大。水饱和砂岩的 ( V_p ) 通常比干燥状态低或持平因为密度增加的效果往往超过模量增加。解决用Wood公式计算混合流体的等效体积模量( 1/K_{fl} S_w/K_w (1-S_w)/K_g )其中 ( S_w ) 是含水饱和度。不要直接用纯水的 ( K_{fl} ) 代替混合流体。4.4 现象剪切模量在饱和后发生变化原因代码中误用了饱和体积模量去反算剪切模量或者把 ( M_{sat} ) 和 ( K_{sat} ) 的关系搞混了。解决Gassmann理论的核心结论之一就是 ( \mu_{sat} \mu_{dry} )。在代码里把这一行单独写出来并加注释任何情况下都不要让剪切模量参与流体替换的迭代。4.5 现象干燥和饱和参数在交会图上完全重叠原因孔隙度输入为0或者 ( K_{fl} ) 等于 ( K_{min} )。前者说明数据里根本没有孔隙后者说明流体和矿物一样硬这两种情况在实际储层中都不存在。解决检查 ( \phi ) 是否从测井曲线正确读取检查 ( K_{fl} ) 是否误用了矿物模量。如果确实是非孔隙岩石如致密石英岩那干燥和饱和参数本来就应该一样不需要做流体替换。5. 进阶技巧从单点计算到曲线批量处理与敏感性分析单点计算跑通之后实际项目里面对的是整条测井曲线。我一般会把上面的逻辑封装成函数输入 ( V_p )、( V_s )、( \rho )、( \phi ) 的数组输出干燥和饱和的 ( M )、( \mu )、( u ) 曲线。这里有一个效率上的技巧Gassmann方程的分母在每条曲线上都要重新计算但 ( K_{min} ) 和 ( K_{fl} ) 通常是常数可以提前算好中间项。def gassmann_batch(Vp_dry, Vs_dry, rho_dry, phi, K_min, K_fl): 批量计算干燥和饱和弹性参数 输入均为numpy数组模量单位GPa mu_dry rho_dry * Vs_dry**2 / 1e9 M_dry rho_dry * Vp_dry**2 / 1e9 K_dry M_dry - 4.0/3.0 * mu_dry nu_dry (M_dry - 2*mu_dry) / (2*(M_dry - mu_dry)) # 预计算不随深度变化的项 K_ratio K_dry / K_min numerator (1 - K_ratio)**2 denominator phi / K_fl (1 - phi) / K_min - K_dry / K_min**2 K_sat K_dry numerator / denominator mu_sat mu_dry.copy() M_sat K_sat 4.0/3.0 * mu_sat nu_sat (3*K_sat - 2*mu_sat) / (2*(3*K_sat mu_sat)) return { M_dry: M_dry, mu_dry: mu_dry, nu_dry: nu_dry, M_sat: M_sat, mu_sat: mu_sat, nu_sat: nu_sat, K_dry: K_dry, K_sat: K_sat }批量处理时最常见的翻车点是数组形状不匹配。测井曲线的深度采样间隔可能不同( \phi ) 曲线和速度曲线的采样点必须对齐。我一般先用pandas的reindex对齐深度再转成numpy数组传入。另一个坑是曲线中可能含有空值或负值需要在计算前用掩码过滤否则Gassmann分母会出现除零或负数。敏感性分析是判断结果可靠性的重要手段。( K_{min} ) 的取值对 ( K_{sat} ) 影响很大砂岩中石英含量从80%变到95%( K_{min} ) 可能从32 GPa变到36 GPa最终 ( u_{sat} ) 的差异可以达到0.02到0.03。我习惯对 ( K_{min} )、( K_{fl} ) 和 ( \phi ) 各做±10%的扰动画出 ( u_{sat} ) 的包络线。如果包络线宽度超过0.05说明输入参数的不确定性已经主导了结果这时候需要回到岩石物理测量本身去提高精度而不是在计算环节反复调参。还有一个实用技巧把干燥和饱和的 ( M )、( \mu )、( u ) 曲线与实测饱和曲线对比。如果实测数据来自水饱和岩心而计算曲线在泥质含量高的层段偏差明显优先检查 ( K_{min} ) 是否应该用混合矿物的等效模量而不是纯石英的值。黏土矿物的体积模量通常在20到25 GPa之间比石英低不少忽略黏土的影响会让 ( K_{min} ) 偏高进而低估 ( K_{sat} )。最后说一个我自己的习惯每次做完流体替换都会把干燥和饱和的 ( V_p/V_s ) 比值曲线画在一起。如果两条曲线在非储层段几乎重合、在储层段明显分离说明计算逻辑是对的。如果全井段都分离或者全井段都重合那大概率是 ( \phi ) 曲线或者 ( K_{fl} ) 出了问题。这个检查花不了两分钟但能省下后面正演阶段大量的返工时间。希望帮到你。本文还有配套的精品资源点击获取