简介面向地震岩石物理与地质勘探学习者提供基于岩石物理基本理论计算干燥岩石与流体饱和岩石弹性参数的MATLAB算法脚本。资源聚焦纵波模量、剪切模量与泊松比三类核心参数其中纵波模量反映岩石抗体积变形能力剪切模量控制横波传播速度泊松比揭示应力作用下的横向变形特征可帮助研究人员快速掌握数值计算流程适用于储层预测、地震活动分析等场景。包体为1个rar压缩包内含1个m文件大小仅2KB脚本为MATLAB可执行程序包含完整的计算逻辑与参数定义结构紧凑便于直接运行、调试或二次开发。已有545人学习下载。通过该脚本可直观理解岩石弹性参数的物理含义与计算实现方式为后续结合云计算开展大规模数值模拟、评估地下储层特性提供基础工具适合科研人员、学生及工程师学习参考。1. 干燥岩石与流体饱和岩石的弹性参数为什么值得亲手算一遍在地震岩石物理里我们经常要计算干燥岩石、流体饱和岩石的弹性参数比如纵波模量、剪切模量和泊松比。表面上看这只是三个公式套数值但真正做过储层预测的人都知道直接把反演得到的 Vp、Vs 代入ρVp²只能得到一个“流体和骨架混在一起”的表观模量既说不清孔隙里到底是油还是水也判断不了岩性变化。更麻烦的是同一套速度数据干燥骨架和流体饱和骨架对应的泊松比往往差出 0.05 到 0.15这个差量恰好是流体识别的敏感区间。把干燥岩石和流体饱和岩石分开计算本质上是给地下介质做一个“解耦”先得到不受流体影响的骨架参数再通过流体替换回到目标饱和度。这条路一旦走通测井解释里的含水饱和度、地震反演里的弹性参数体、以及油藏工程里对压力衰竭的预测就都有了共同的底层输入。适合正在做岩石物理建模、测井评价或者叠前反演标定的工程师。2. 先厘清纵波模量、剪切模量与泊松比三个参数各自听谁的“指挥”2.1 纵波模量 M、剪切模量 μ、泊松比 ν 的定义与换算纵波模量 M 在中文文献里也常叫 P 波模量定义为纵波速度的平方乘以密度M ρ × Vp²。它和体积模量 K 不是一回事M 包含了对体积变化和形状变化的共同抵抗能力数学上写作M K 4/3 μ。剪切模量 μ 则只描述形状变化表达式是μ ρ × Vs²几乎所有岩石物理公式都对 μ 格外“偏爱”因为流体不承担剪切应力。泊松比 ν 描述横向压缩与纵向拉伸的比例用速度表达时最常用这个形式ν (Vp² - 2Vs²) / (2(Vp² - Vs²))Vp/Vs 越高ν 越大。水饱和砂岩的 ν 通常在 0.30 到 0.35气饱和或干燥砂岩往往降到 0.20 到 0.25因此泊松比常被当成流体指示器。需要注意当速度单位用 km/s 时上面的公式不受影响因为单位会被约掉但一旦把速度代入 M 和 μ 的表达式密度单位必须和速度单位匹配速度用 m/s 时密度用 kg/m³模量单位是 Pa速度用 km/s 时密度用 g/cm³模量单位是 GPa。很多初学翻车就是在这里速度用 km/s密度用 g/cm³算出来的 M 成了 GPa 的量级却把它当作 Pa 去和文献对比。三个参数之间还可以互相换算。如果已经算出了 M 和 μ泊松比可以用ν (M - 2μ)/(2(M - μ))计算体积模量则是K M - 4/3 μ。实际项目里我一般先用 Vp、Vs、密度算出 M 和 μ再用速度形式算泊松比因为后一步可以直接和测井的 Vp/Vs 曲线对比便于发现输入数据是否错位。2.2 为什么干燥岩石和饱和岩石必须分开算Gassmann 方程的“锚点”是干骨架Gassmann 方程是地震岩石物理流体替换的基石它回答了一个关键问题已知干燥岩石的弹性模量孔隙中充满流体后体积模量会变成多少。方程的低频形式是K_sat K_dry (1 - K_dry/K_min)² / (φ/K_fluid (1 - φ)/K_min - K_dry/K_min²)其中K_sat是流体饱和岩石的体积模量K_dry是干燥岩石骨架的体积模量K_min是固体矿物的体积模量K_fluid是孔隙流体的体积模量φ是孔隙度。方程里有一个非常关键的结果剪切模量在流体替换前后不变即μ_sat μ_dry。这是因为 Gassmann 假设孔隙压力平衡、低频条件下流体不承担剪切力流体只改变体积模量不改变剪切模量。这个性质决定了计算顺序必须先确定一套“与流体无关”的干骨架模量再代入流体。很多反演项目遇到的问题是井上能直接测得的是饱和岩石的速度和密度从中能反推的是K_sat可真正对含水饱和度敏感的是K_dry。如果把K_sat当成骨架参数去做岩性判别含油层和含水层的差异会被流体效应放大误判率会明显上升。Gassmann 方程的价值就是允许我们把测量到的饱和岩石参数“剥离流体”还原成一个虚构的干燥岩石状态然后再替换到任意目标流体。2.3 一次完整的干燥→饱和计算应该经过哪几步从实际操作的顺序看一条完整的计算链在程序里是线性展开的每步的结果都要作为下一步的输入。第一步是确定固体矿物的弹性模量通常用 Voigt-Reuss-Hill 平均把石英、黏土、长石等矿物混成一套等效矿物。第二步是选择干骨架模型比如临界孔隙度模型或 Krief 经验公式由矿物模量和孔隙度算出K_dry与μ_dry。第三步是准备流体参数按目标流体水、油、气得到体积模量和密度。第四步用 Gassmann 方程算出K_sat同时保持μ_sat μ_dry再用流体替换后的混合密度计算新的 Vp、Vs最后换算泊松比。看起来很线性但每一步都有参数选择风险矿物模量取错了后面全偏干骨架模型选得不合适流体替换的幅度会失真流体参数用错温压条件气体等效体积模量能差出一个数量级。所以后面几章我会把每步的典型参数、公式和代码分别展开并把最常见的踩坑点单独列出来。3. 计算干燥岩石的弹性参数矿物平均与干骨架模型怎么选3.1 用 Voigt-Reuss-Hill 把石英、黏土、长石混成“一块矿物”地下岩石不是单一矿物砂岩是石英、长石、黏土矿物的混合体泥岩则以黏土为主。Gassmann 方程里的K_min和μ_min必须是等效矿物模量直接用某种纯矿物模量代替会带来明显偏差。最常用的混合方法是 Voigt-Reuss-Hill 平均Voigt 平均假设应变均匀实际是等应变上限Reuss 平均假设应力均匀是等应力下限两者取算术平均就是 Hill 平均真实矿物骨架通常落在这个上下限之间。实现代码不需要很复杂关键是同时处理体积模量和剪切模量两个序列。一个可复现的最小函数如下def vrh_average(fractions, k_list, mu_list): # fractions: 各矿物体积分数例如 [0.65, 0.25, 0.10] # k_list: 对应矿物体积模量单位 GPa # mu_list: 对应矿物剪切模量单位 GPa # Voigt 上限体积分数加权求和 kv sum(f * k for f, k in zip(fractions, k_list)) gv sum(f * g for f, g in zip(fractions, mu_list)) # Reuss 下限体积分数除以模量后再求倒数 kr 1.0 / sum(f / k for f, k in zip(fractions, k_list)) gr 1.0 / sum(f / g for f, g in zip(fractions, mu_list)) # Hill 平均上下限取算术平均 k_hill (kv kr) / 2.0 mu_hill (gv gr) / 2.0 return k_hill, mu_hill调用时只需要传入矿物分数和模量表。以常见的石英砂岩混入少量黏土为例石英K37 GPaμ44 GPa黏土K21 GPaμ7 GPa长石K38 GPaμ15 GPa。如果砂岩中含 70% 石英、20% 长石、10% 黏土那么fracs [0.70, 0.20, 0.10] k_values [37.0, 38.0, 21.0] mu_values [44.0, 15.0, 7.0] k_min, mu_min vrh_average(fracs, k_values, mu_values) print(K_min , k_min, GPa, mu_min , mu_min, GPa)这段代码里最重要的一点是分数之和必须等于 1否则上下限都会失真。实际测井中矿物分数来自元素俘获谱或 XRD 分析这些数据往往有误差因此我建议对黏土含量做敏感性测试黏土含量从 5% 变到 25%等效K_min可能下降 4-8 GPa足够让后续流体替换的 Vp 偏差超过 50 m/s。矿物模量表中的参数不要照抄所有文献优先选与本区成岩环境相近的数值。3.2 干骨架模量临界孔隙度模型与 Krief 经验式二选一有了等效矿物模量下一步是算干骨架。这里不存在从地震数据直接测得的干骨架必须通过模型或实验室测量约束。最常用的是临界孔隙度模型它假设骨架模量随孔隙度线性降低在孔隙度到达临界值φ_c时模量降为 0固体颗粒不再有承力结构。砂岩的临界孔隙度通常取 0.36-0.40。表达式是K_dry K_min × (1 - φ/φ_c) μ_dry μ_min × (1 - φ/φ_c)这个模型的优点是参数少、物理含义直观缺点是只适合中高孔隙度区间。孔隙度小于 0.05 时线性公式会把骨架模量推向接近矿物模量实际致密岩石可能会因为微裂缝提前软化所以致密储层往往不适用。另一条常用路线是 Krief 经验式它在低孔隙度下衰减更慢在高孔隙度下和临界孔隙度模型趋于一致n 3 / (1 - φ) K_dry K_min × (1 - φ)^n μ_dry μ_min × (1 - φ)^n两条模型放一起时会发现一个典型差异在中孔10%-20%区间Krief 给出的K_dry通常比临界孔隙度模型高 10%-20%这会造成流体替换后的泊松比差 0.02 左右。选择时不能只看拟合优度更要看目标是流体识别还是岩性识别流体识别需要尽量保守的干骨架刚度岩性识别则要贴近真实岩石骨架。下面是两个函数的最小实现def dry_moduli_critical(k_min, mu_min, phi, phi_c0.4): # 临界孔隙度模型φ φ_c 时骨架失去刚度 if phi phi_c: return 0.0, 0.0 k_dry k_min * (1.0 - phi / phi_c) mu_dry mu_min * (1.0 - phi / phi_c) return k_dry, mu_dry def dry_moduli_krief(k_min, mu_min, phi): # Krief 经验式指数与孔隙度相关 if phi 1.0 or phi 0.0: return None n 3.0 / (1.0 - phi) factor (1.0 - phi) ** n return k_min * factor, mu_min * factor用临界孔隙度模型时有一个必须处理的边界情况当φ ≥ φ_c代码直接返回 0 模量。如果不做这个判断1 - phi/phi_c变成负数后续 Gassmann 方程分母会出现非物理值。很多批量处理脚本没有这个保护结果会在高孔隙度段输出异常的负模量后面算速度时直接开根号报错。3.3 干燥岩石的泊松比不是额外输入而是 K_dry/μ_dry 的结果计算干燥岩石参数时常见的误会是“先算干燥泊松比”。其实干燥岩石的泊松比是派生的由K_dry和μ_dry共同决定而不是一个独立输入。根据弹性参数关系干燥岩石体积模量与剪切模量之比K_dry/μ_dry一旦确定泊松比就唯一确定。这意味着如果你在代码里给干燥岩石单独指定一个泊松比又同时指定K_dry、μ_dry三者的关系很可能是矛盾的后续 Gassmann 替换后 Vp/Vs 也会表现出人为倾向。处理干骨架的合理做法是只挑一个主控制参数通常选μ_dry或K_dry另一个由模型给出。临界孔隙度模型同时给出两个天然自洽。如果是手工作业里只有一条干燥泊松比曲线那么要先用ν_dry反推K_dry/μ_dry再结合其中一个模量确定另一个千万不要把两个模量和泊松比全部平行输入。4. 流体饱和岩石的参数计算把 Gassmann 流体替换跑起来4.1 水、油、气的体积模量和密度取数时别只看“一个数”流体参数是整个计算链里最容易让人低估动态变化的一环。纯水的体积模量在常温常压下大约是 2.2 GPa密度 1.0 g/cm³但地层水含盐后密度会升到 1.05-1.10 g/cm³体积模量也升高到 2.4-2.8 GPa。油的体积模量受溶解气和温度影响很大轻质油一般在 0.8-1.5 GPa稠油可能低于 0.6 GPa。天然气最极端地面条件下体积模量几乎为 0但在 3000 米深处、温度 90°C、压力 30 MPa 时甲烷的体积模量可以达到 0.03-0.08 GPa密度 0.15-0.30 g/cm³。如果拿地面条件下的气体参数去算流体替换会把气层响应压到几乎不可见。工程上常见的取值方式是使用 Batzle-Wang 经验公式输入温度、压力、气油比和矿化度来计算每种流体的体积模量与密度。手工作业或者早期评价中也可以用简化参数表盐水K2.3 GPa, ρ1.05 g/cm³轻质油K1.0 GPa, ρ0.80 g/cm³天然气储层条件K0.05 GPa, ρ0.20 g/cm³。关键是统一单位GPa配合g/cm³时速度单位必须用 km/s密度不要混用 kg/m³。4.2 从干燥到饱和的 Python 最小实现Gassmann 公式与密度更新下面这段代码把第 3 章算出的干骨架模量带入 Gassmann 方程同时更新密度并输出饱和岩石的速度与泊松比。这是整套计算的核心建议直接保留到项目脚本里import math def fluid_substitution(k_dry, mu_dry, k_min, mu_min, k_fluid, phi, rho_min, rho_fluid): # 输入干骨架体积模量/剪切模量单位 GPa # 矿物模量、流体模量单位 GPa # 孔隙度小数矿物密度、流体密度单位 g/cm3 # Gassmann 方程求饱和岩石体积模量 denominator (phi / k_fluid (1.0 - phi) / k_min - k_dry / (k_min * k_min)) k_sat k_dry (1.0 - k_dry / k_min) ** 2 / denominator # 剪切模量不受流体影响 mu_sat mu_dry # 饱和岩石密度骨架 孔隙流体 rho_sat (1.0 - phi) * rho_min phi * rho_fluid # 从模量换算 Vp、Vs注意模量 GPa 与密度 g/cm3 对应速度 km/s vp math.sqrt((k_sat 4.0 / 3.0 * mu_sat) / rho_sat) vs math.sqrt(mu_sat / rho_sat) # 泊松比用 Vp/Vs 表示 ratio_vp_vs (vp / vs) ** 2 nu (ratio_vp_vs - 2.0) / (2.0 * (ratio_vp_vs - 1.0)) return k_sat, mu_sat, rho_sat, vp, vs, nu调用时先算矿物等效模量再算干骨架然后选择一个目标流体饱和度。比如一个中孔砂岩孔隙度 0.25矿物模量K_min36.5 GPa, μ_min40 GPa矿物密度 2.65含水饱和k_min, mu_min 36.5, 40.0 rho_min 2.65 phi 0.25 k_dry, mu_dry dry_moduli_critical(k_min, mu_min, phi, phi_c0.4) result fluid_substitution(k_dry, mu_dry, k_min, mu_min, 2.3, phi, rho_min, 1.05) print(Vp%.3f km/s Vs%.3f km/s nu%.3f % (result[3], result[4], result[5]))注意代码里denominator的构成第一项φ/K_fluid是流体所占的柔度第二项(1-φ)/K_min是矿物柔度第三项K_dry/K_min²是骨架柔度的修正。这个分母数值上通常很小所以K_sat对浮点误差敏感。如果使用 32 位浮点运算在高孔隙度低流体模量组合下可能出现分母为零的异常建议在项目里统一转成 64 位浮点。4.3 从饱和模量反算 Vp、Vs、泊松比顺序错了就全错很多人会在 Gassmann 替换之后直接拿原来的速度曲线套泊松比这是一个方向性错误。流体替换改变的不仅是体积模量孔隙流体密度变化还会让饱和岩石整体密度改变。举例来说同一套干骨架含水后密度比含气后高 0.15-0.25 g/cm³这个差异会直接进入 Vp、Vs 的计算速度等于模量除以密度的平方根密度变化 10%速度就要变化 5% 左右。如果沿用原密度等于把替换后的模量配到错误的惯性项上结果速度和泊松比都会向错误方向偏移。正确顺序永远是这样先算K_sat保持μ_sat不变再计算ρ_sat最后通过Vp sqrt((K_sat 4/3 μ_sat)/ρ_sat)和Vs sqrt(μ_sat/ρ_sat)得到速度再转泊松比。这个顺序在上一节代码里已经实现。实践里我还会额外输出中间量K_sat和ρ_sat直接和测井反演的动态体积模量对比这一步经常能暴露出密度曲线质量问题比只对比 Vp 更敏感。5. 避坑清单流体替换计算结果不物理的几种常见翻车现场5.1 现象泊松比大于 0.5 或变成负数原因最常见的是 Vp/Vs 输入不物理。如果 Vp 和 Vs 来自不同深度采样点或者 Vs 曲线有异常尖峰速度比的平方就可能小于 2泊松比公式分子变负。另一种情况是单位换算错位速度用 m/s、密度用 g/cm³算出的模量混夹着不同量纲最终泊松比完全失真。解决先对 Vp、Vs 做深度匹配和异常值剔除再用一个已知水层的 Vp/Vs 区间做合理性检查。如果水层泊松比不在 0.30-0.36优先怀疑输入速度曲线而不是修改公式。另外在代码里做一个边界保护任何ν 0.5或ν 0的样本直接标红输出不要静默处理。5.2 现象流体替换前后剪切模量变了原因Gassmann 方程要求μ_sat μ_dry但很多人会在写代码时把μ_sat写成由矿物模量和孔隙度重新计算或者套用体积模量替换公式把剪切模量也做了流体替换。剪切模量对流体的响应在地震频带里确实可以忽略流体替换后它变了就说明代码路径有误。解决检查替换函数里是否有mu_sat mu_dry这行赋值。如果没有直接补上。同时打印替换前后的剪切模量曲线两条线理论上应完全重合任何偏差都是实现错误。5.3 现象孔隙度大于临界孔隙度时干骨架模量变成负数原因临界孔隙度模型线性公式在φ φ_c时没有物理意义但批量脚本往往扫过全井段未固结高孔段就会出现负模量。负模量进入 Gassmann 方程后甚至可能生成虚数速度。解决临界孔隙度函数里增加if phi phi_c: return 0.0, 0.0把这段孔隙度标记为骨架失效。如果实际目标储层孔隙度确实超过临界孔隙度说明应该换用更复杂的未固结砂岩模型比如接触理论模型而不是继续死磕临界孔隙度公式。5.4 现象气体参数用的是地面条件导致流体替换失真原因天然气的体积模量和密度对温度压力极度敏感地面条件近似 0.0001 GPa 与储层条件下 0.03-0.08 GPa 相差两个数量级。用错参数后气层的 Gassmann 流体响应会非常弱甚至出现含气层和水层的 Vp 几乎一致的现象。解决储层流体参数一定要输入温压条件用 Batzle-Wang 公式或 PVT 实验数据计算。如果项目早期没有 PVT 资料我一般至少会做一次敏感性分析气体体积模量取 0.02、0.05、0.10 GPa 三个值看目标储层泊松比变化是否超过 0.02。如果超过说明该项目不能靠简化参数必须补 PVT 数据。5.5 现象速度、密度、孔隙度曲线深度不对齐原因声波测井、密度测井和中子孔隙度测井来自不同仪器本身就有深度误差和采样率差异。没有做深度匹配和重采样就直接进入流体替换会在地层边界上产生抛物线假象特别是薄储层中Vp 可能来自相邻泥岩层孔隙度却来自储层。解决流体替换前先做多曲线深度对齐再统一重采样到同一采样率。做完这一步后用交会图质量控制孔隙度和密度交会应该沿干净趋势分布出现明显离群点说明深度匹配还没到位。6. 用实际井点做一次标定反推临界孔隙度并检验流体替换计算流程跑通之后不要直接信任默认参数。我一般会用一口有可靠测井解释的目标井做标定已知含水饱和度、孔隙度、矿物体积分数以及实测 Vp、Vs然后反推干骨架模型里的关键参数比如临界孔隙度φ_c。目标函数很简单用选定的φ_c做流体替换算出的 Vs 或 Vp/Vs 与实测曲线最接近的那个值就是本区比较可信的骨架标定参数。def calibrate_phi_c(vp_obs, vs_obs, phi, k_min, mu_min, k_fluid, rho_min, rho_fluid): # vp_obs/vs_obs: 实测纵波/横波速度km/s # phi: 孔隙度数组rho_fluid: 目标流体密度 best_phi_c 0.4 best_error 1e9 for phi_c in [0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.45]: total_error 0.0 for i in range(len(phi)): k_dry, mu_dry dry_moduli_critical(k_min, mu_min, phi[i], phi_c) result fluid_substitution(k_dry, mu_dry, k_min, mu_min, k_fluid, phi[i], rho_min, rho_fluid) vp_calc, vs_calc result[3], result[4] total_error (vp_calc - vp_obs[i]) ** 2 total_error ((vs_calc - vs_obs[i]) * 2) ** 2 if total_error best_error: best_error total_error best_phi_c phi_c return best_phi_c, best_error代码里给 Vs 误差乘了权重 2因为 Vs 对流体不敏感能更好约束骨架参数如果只拟合 Vp很可能得到多个等价的φ_c。标定结果出来后再看替换后泊松比曲线和实测泊松比曲线的残差残差小于 0.02 说明参数组合可靠残差系统偏移则说明矿物模量或流体模量仍然有偏差而不是骨架模型的问题。这套标定方法不需要额外软件适合在储层评价早期把干骨架模型固定下来。多年做岩石物理项目的习惯让我坚持一点任何弹性参数计算结果都必须先接受实际测井曲线的交叉验证否则参数再合理也只是“看起来合理”而真正能经受井点检验的计算链后期做流体预测时才会少返工。希望这套从干燥岩石到流体饱和岩石的计算流程能帮你在自己的数据上少踩几个坑。本文还有配套的精品资源点击获取
从干燥到饱和：Gassmann流体替换与弹性参数计算指南
简介面向地震岩石物理与地质勘探学习者提供基于岩石物理基本理论计算干燥岩石与流体饱和岩石弹性参数的MATLAB算法脚本。资源聚焦纵波模量、剪切模量与泊松比三类核心参数其中纵波模量反映岩石抗体积变形能力剪切模量控制横波传播速度泊松比揭示应力作用下的横向变形特征可帮助研究人员快速掌握数值计算流程适用于储层预测、地震活动分析等场景。包体为1个rar压缩包内含1个m文件大小仅2KB脚本为MATLAB可执行程序包含完整的计算逻辑与参数定义结构紧凑便于直接运行、调试或二次开发。已有545人学习下载。通过该脚本可直观理解岩石弹性参数的物理含义与计算实现方式为后续结合云计算开展大规模数值模拟、评估地下储层特性提供基础工具适合科研人员、学生及工程师学习参考。1. 干燥岩石与流体饱和岩石的弹性参数为什么值得亲手算一遍在地震岩石物理里我们经常要计算干燥岩石、流体饱和岩石的弹性参数比如纵波模量、剪切模量和泊松比。表面上看这只是三个公式套数值但真正做过储层预测的人都知道直接把反演得到的 Vp、Vs 代入ρVp²只能得到一个“流体和骨架混在一起”的表观模量既说不清孔隙里到底是油还是水也判断不了岩性变化。更麻烦的是同一套速度数据干燥骨架和流体饱和骨架对应的泊松比往往差出 0.05 到 0.15这个差量恰好是流体识别的敏感区间。把干燥岩石和流体饱和岩石分开计算本质上是给地下介质做一个“解耦”先得到不受流体影响的骨架参数再通过流体替换回到目标饱和度。这条路一旦走通测井解释里的含水饱和度、地震反演里的弹性参数体、以及油藏工程里对压力衰竭的预测就都有了共同的底层输入。适合正在做岩石物理建模、测井评价或者叠前反演标定的工程师。2. 先厘清纵波模量、剪切模量与泊松比三个参数各自听谁的“指挥”2.1 纵波模量 M、剪切模量 μ、泊松比 ν 的定义与换算纵波模量 M 在中文文献里也常叫 P 波模量定义为纵波速度的平方乘以密度M ρ × Vp²。它和体积模量 K 不是一回事M 包含了对体积变化和形状变化的共同抵抗能力数学上写作M K 4/3 μ。剪切模量 μ 则只描述形状变化表达式是μ ρ × Vs²几乎所有岩石物理公式都对 μ 格外“偏爱”因为流体不承担剪切应力。泊松比 ν 描述横向压缩与纵向拉伸的比例用速度表达时最常用这个形式ν (Vp² - 2Vs²) / (2(Vp² - Vs²))Vp/Vs 越高ν 越大。水饱和砂岩的 ν 通常在 0.30 到 0.35气饱和或干燥砂岩往往降到 0.20 到 0.25因此泊松比常被当成流体指示器。需要注意当速度单位用 km/s 时上面的公式不受影响因为单位会被约掉但一旦把速度代入 M 和 μ 的表达式密度单位必须和速度单位匹配速度用 m/s 时密度用 kg/m³模量单位是 Pa速度用 km/s 时密度用 g/cm³模量单位是 GPa。很多初学翻车就是在这里速度用 km/s密度用 g/cm³算出来的 M 成了 GPa 的量级却把它当作 Pa 去和文献对比。三个参数之间还可以互相换算。如果已经算出了 M 和 μ泊松比可以用ν (M - 2μ)/(2(M - μ))计算体积模量则是K M - 4/3 μ。实际项目里我一般先用 Vp、Vs、密度算出 M 和 μ再用速度形式算泊松比因为后一步可以直接和测井的 Vp/Vs 曲线对比便于发现输入数据是否错位。2.2 为什么干燥岩石和饱和岩石必须分开算Gassmann 方程的“锚点”是干骨架Gassmann 方程是地震岩石物理流体替换的基石它回答了一个关键问题已知干燥岩石的弹性模量孔隙中充满流体后体积模量会变成多少。方程的低频形式是K_sat K_dry (1 - K_dry/K_min)² / (φ/K_fluid (1 - φ)/K_min - K_dry/K_min²)其中K_sat是流体饱和岩石的体积模量K_dry是干燥岩石骨架的体积模量K_min是固体矿物的体积模量K_fluid是孔隙流体的体积模量φ是孔隙度。方程里有一个非常关键的结果剪切模量在流体替换前后不变即μ_sat μ_dry。这是因为 Gassmann 假设孔隙压力平衡、低频条件下流体不承担剪切力流体只改变体积模量不改变剪切模量。这个性质决定了计算顺序必须先确定一套“与流体无关”的干骨架模量再代入流体。很多反演项目遇到的问题是井上能直接测得的是饱和岩石的速度和密度从中能反推的是K_sat可真正对含水饱和度敏感的是K_dry。如果把K_sat当成骨架参数去做岩性判别含油层和含水层的差异会被流体效应放大误判率会明显上升。Gassmann 方程的价值就是允许我们把测量到的饱和岩石参数“剥离流体”还原成一个虚构的干燥岩石状态然后再替换到任意目标流体。2.3 一次完整的干燥→饱和计算应该经过哪几步从实际操作的顺序看一条完整的计算链在程序里是线性展开的每步的结果都要作为下一步的输入。第一步是确定固体矿物的弹性模量通常用 Voigt-Reuss-Hill 平均把石英、黏土、长石等矿物混成一套等效矿物。第二步是选择干骨架模型比如临界孔隙度模型或 Krief 经验公式由矿物模量和孔隙度算出K_dry与μ_dry。第三步是准备流体参数按目标流体水、油、气得到体积模量和密度。第四步用 Gassmann 方程算出K_sat同时保持μ_sat μ_dry再用流体替换后的混合密度计算新的 Vp、Vs最后换算泊松比。看起来很线性但每一步都有参数选择风险矿物模量取错了后面全偏干骨架模型选得不合适流体替换的幅度会失真流体参数用错温压条件气体等效体积模量能差出一个数量级。所以后面几章我会把每步的典型参数、公式和代码分别展开并把最常见的踩坑点单独列出来。3. 计算干燥岩石的弹性参数矿物平均与干骨架模型怎么选3.1 用 Voigt-Reuss-Hill 把石英、黏土、长石混成“一块矿物”地下岩石不是单一矿物砂岩是石英、长石、黏土矿物的混合体泥岩则以黏土为主。Gassmann 方程里的K_min和μ_min必须是等效矿物模量直接用某种纯矿物模量代替会带来明显偏差。最常用的混合方法是 Voigt-Reuss-Hill 平均Voigt 平均假设应变均匀实际是等应变上限Reuss 平均假设应力均匀是等应力下限两者取算术平均就是 Hill 平均真实矿物骨架通常落在这个上下限之间。实现代码不需要很复杂关键是同时处理体积模量和剪切模量两个序列。一个可复现的最小函数如下def vrh_average(fractions, k_list, mu_list): # fractions: 各矿物体积分数例如 [0.65, 0.25, 0.10] # k_list: 对应矿物体积模量单位 GPa # mu_list: 对应矿物剪切模量单位 GPa # Voigt 上限体积分数加权求和 kv sum(f * k for f, k in zip(fractions, k_list)) gv sum(f * g for f, g in zip(fractions, mu_list)) # Reuss 下限体积分数除以模量后再求倒数 kr 1.0 / sum(f / k for f, k in zip(fractions, k_list)) gr 1.0 / sum(f / g for f, g in zip(fractions, mu_list)) # Hill 平均上下限取算术平均 k_hill (kv kr) / 2.0 mu_hill (gv gr) / 2.0 return k_hill, mu_hill调用时只需要传入矿物分数和模量表。以常见的石英砂岩混入少量黏土为例石英K37 GPaμ44 GPa黏土K21 GPaμ7 GPa长石K38 GPaμ15 GPa。如果砂岩中含 70% 石英、20% 长石、10% 黏土那么fracs [0.70, 0.20, 0.10] k_values [37.0, 38.0, 21.0] mu_values [44.0, 15.0, 7.0] k_min, mu_min vrh_average(fracs, k_values, mu_values) print(K_min , k_min, GPa, mu_min , mu_min, GPa)这段代码里最重要的一点是分数之和必须等于 1否则上下限都会失真。实际测井中矿物分数来自元素俘获谱或 XRD 分析这些数据往往有误差因此我建议对黏土含量做敏感性测试黏土含量从 5% 变到 25%等效K_min可能下降 4-8 GPa足够让后续流体替换的 Vp 偏差超过 50 m/s。矿物模量表中的参数不要照抄所有文献优先选与本区成岩环境相近的数值。3.2 干骨架模量临界孔隙度模型与 Krief 经验式二选一有了等效矿物模量下一步是算干骨架。这里不存在从地震数据直接测得的干骨架必须通过模型或实验室测量约束。最常用的是临界孔隙度模型它假设骨架模量随孔隙度线性降低在孔隙度到达临界值φ_c时模量降为 0固体颗粒不再有承力结构。砂岩的临界孔隙度通常取 0.36-0.40。表达式是K_dry K_min × (1 - φ/φ_c) μ_dry μ_min × (1 - φ/φ_c)这个模型的优点是参数少、物理含义直观缺点是只适合中高孔隙度区间。孔隙度小于 0.05 时线性公式会把骨架模量推向接近矿物模量实际致密岩石可能会因为微裂缝提前软化所以致密储层往往不适用。另一条常用路线是 Krief 经验式它在低孔隙度下衰减更慢在高孔隙度下和临界孔隙度模型趋于一致n 3 / (1 - φ) K_dry K_min × (1 - φ)^n μ_dry μ_min × (1 - φ)^n两条模型放一起时会发现一个典型差异在中孔10%-20%区间Krief 给出的K_dry通常比临界孔隙度模型高 10%-20%这会造成流体替换后的泊松比差 0.02 左右。选择时不能只看拟合优度更要看目标是流体识别还是岩性识别流体识别需要尽量保守的干骨架刚度岩性识别则要贴近真实岩石骨架。下面是两个函数的最小实现def dry_moduli_critical(k_min, mu_min, phi, phi_c0.4): # 临界孔隙度模型φ φ_c 时骨架失去刚度 if phi phi_c: return 0.0, 0.0 k_dry k_min * (1.0 - phi / phi_c) mu_dry mu_min * (1.0 - phi / phi_c) return k_dry, mu_dry def dry_moduli_krief(k_min, mu_min, phi): # Krief 经验式指数与孔隙度相关 if phi 1.0 or phi 0.0: return None n 3.0 / (1.0 - phi) factor (1.0 - phi) ** n return k_min * factor, mu_min * factor用临界孔隙度模型时有一个必须处理的边界情况当φ ≥ φ_c代码直接返回 0 模量。如果不做这个判断1 - phi/phi_c变成负数后续 Gassmann 方程分母会出现非物理值。很多批量处理脚本没有这个保护结果会在高孔隙度段输出异常的负模量后面算速度时直接开根号报错。3.3 干燥岩石的泊松比不是额外输入而是 K_dry/μ_dry 的结果计算干燥岩石参数时常见的误会是“先算干燥泊松比”。其实干燥岩石的泊松比是派生的由K_dry和μ_dry共同决定而不是一个独立输入。根据弹性参数关系干燥岩石体积模量与剪切模量之比K_dry/μ_dry一旦确定泊松比就唯一确定。这意味着如果你在代码里给干燥岩石单独指定一个泊松比又同时指定K_dry、μ_dry三者的关系很可能是矛盾的后续 Gassmann 替换后 Vp/Vs 也会表现出人为倾向。处理干骨架的合理做法是只挑一个主控制参数通常选μ_dry或K_dry另一个由模型给出。临界孔隙度模型同时给出两个天然自洽。如果是手工作业里只有一条干燥泊松比曲线那么要先用ν_dry反推K_dry/μ_dry再结合其中一个模量确定另一个千万不要把两个模量和泊松比全部平行输入。4. 流体饱和岩石的参数计算把 Gassmann 流体替换跑起来4.1 水、油、气的体积模量和密度取数时别只看“一个数”流体参数是整个计算链里最容易让人低估动态变化的一环。纯水的体积模量在常温常压下大约是 2.2 GPa密度 1.0 g/cm³但地层水含盐后密度会升到 1.05-1.10 g/cm³体积模量也升高到 2.4-2.8 GPa。油的体积模量受溶解气和温度影响很大轻质油一般在 0.8-1.5 GPa稠油可能低于 0.6 GPa。天然气最极端地面条件下体积模量几乎为 0但在 3000 米深处、温度 90°C、压力 30 MPa 时甲烷的体积模量可以达到 0.03-0.08 GPa密度 0.15-0.30 g/cm³。如果拿地面条件下的气体参数去算流体替换会把气层响应压到几乎不可见。工程上常见的取值方式是使用 Batzle-Wang 经验公式输入温度、压力、气油比和矿化度来计算每种流体的体积模量与密度。手工作业或者早期评价中也可以用简化参数表盐水K2.3 GPa, ρ1.05 g/cm³轻质油K1.0 GPa, ρ0.80 g/cm³天然气储层条件K0.05 GPa, ρ0.20 g/cm³。关键是统一单位GPa配合g/cm³时速度单位必须用 km/s密度不要混用 kg/m³。4.2 从干燥到饱和的 Python 最小实现Gassmann 公式与密度更新下面这段代码把第 3 章算出的干骨架模量带入 Gassmann 方程同时更新密度并输出饱和岩石的速度与泊松比。这是整套计算的核心建议直接保留到项目脚本里import math def fluid_substitution(k_dry, mu_dry, k_min, mu_min, k_fluid, phi, rho_min, rho_fluid): # 输入干骨架体积模量/剪切模量单位 GPa # 矿物模量、流体模量单位 GPa # 孔隙度小数矿物密度、流体密度单位 g/cm3 # Gassmann 方程求饱和岩石体积模量 denominator (phi / k_fluid (1.0 - phi) / k_min - k_dry / (k_min * k_min)) k_sat k_dry (1.0 - k_dry / k_min) ** 2 / denominator # 剪切模量不受流体影响 mu_sat mu_dry # 饱和岩石密度骨架 孔隙流体 rho_sat (1.0 - phi) * rho_min phi * rho_fluid # 从模量换算 Vp、Vs注意模量 GPa 与密度 g/cm3 对应速度 km/s vp math.sqrt((k_sat 4.0 / 3.0 * mu_sat) / rho_sat) vs math.sqrt(mu_sat / rho_sat) # 泊松比用 Vp/Vs 表示 ratio_vp_vs (vp / vs) ** 2 nu (ratio_vp_vs - 2.0) / (2.0 * (ratio_vp_vs - 1.0)) return k_sat, mu_sat, rho_sat, vp, vs, nu调用时先算矿物等效模量再算干骨架然后选择一个目标流体饱和度。比如一个中孔砂岩孔隙度 0.25矿物模量K_min36.5 GPa, μ_min40 GPa矿物密度 2.65含水饱和k_min, mu_min 36.5, 40.0 rho_min 2.65 phi 0.25 k_dry, mu_dry dry_moduli_critical(k_min, mu_min, phi, phi_c0.4) result fluid_substitution(k_dry, mu_dry, k_min, mu_min, 2.3, phi, rho_min, 1.05) print(Vp%.3f km/s Vs%.3f km/s nu%.3f % (result[3], result[4], result[5]))注意代码里denominator的构成第一项φ/K_fluid是流体所占的柔度第二项(1-φ)/K_min是矿物柔度第三项K_dry/K_min²是骨架柔度的修正。这个分母数值上通常很小所以K_sat对浮点误差敏感。如果使用 32 位浮点运算在高孔隙度低流体模量组合下可能出现分母为零的异常建议在项目里统一转成 64 位浮点。4.3 从饱和模量反算 Vp、Vs、泊松比顺序错了就全错很多人会在 Gassmann 替换之后直接拿原来的速度曲线套泊松比这是一个方向性错误。流体替换改变的不仅是体积模量孔隙流体密度变化还会让饱和岩石整体密度改变。举例来说同一套干骨架含水后密度比含气后高 0.15-0.25 g/cm³这个差异会直接进入 Vp、Vs 的计算速度等于模量除以密度的平方根密度变化 10%速度就要变化 5% 左右。如果沿用原密度等于把替换后的模量配到错误的惯性项上结果速度和泊松比都会向错误方向偏移。正确顺序永远是这样先算K_sat保持μ_sat不变再计算ρ_sat最后通过Vp sqrt((K_sat 4/3 μ_sat)/ρ_sat)和Vs sqrt(μ_sat/ρ_sat)得到速度再转泊松比。这个顺序在上一节代码里已经实现。实践里我还会额外输出中间量K_sat和ρ_sat直接和测井反演的动态体积模量对比这一步经常能暴露出密度曲线质量问题比只对比 Vp 更敏感。5. 避坑清单流体替换计算结果不物理的几种常见翻车现场5.1 现象泊松比大于 0.5 或变成负数原因最常见的是 Vp/Vs 输入不物理。如果 Vp 和 Vs 来自不同深度采样点或者 Vs 曲线有异常尖峰速度比的平方就可能小于 2泊松比公式分子变负。另一种情况是单位换算错位速度用 m/s、密度用 g/cm³算出的模量混夹着不同量纲最终泊松比完全失真。解决先对 Vp、Vs 做深度匹配和异常值剔除再用一个已知水层的 Vp/Vs 区间做合理性检查。如果水层泊松比不在 0.30-0.36优先怀疑输入速度曲线而不是修改公式。另外在代码里做一个边界保护任何ν 0.5或ν 0的样本直接标红输出不要静默处理。5.2 现象流体替换前后剪切模量变了原因Gassmann 方程要求μ_sat μ_dry但很多人会在写代码时把μ_sat写成由矿物模量和孔隙度重新计算或者套用体积模量替换公式把剪切模量也做了流体替换。剪切模量对流体的响应在地震频带里确实可以忽略流体替换后它变了就说明代码路径有误。解决检查替换函数里是否有mu_sat mu_dry这行赋值。如果没有直接补上。同时打印替换前后的剪切模量曲线两条线理论上应完全重合任何偏差都是实现错误。5.3 现象孔隙度大于临界孔隙度时干骨架模量变成负数原因临界孔隙度模型线性公式在φ φ_c时没有物理意义但批量脚本往往扫过全井段未固结高孔段就会出现负模量。负模量进入 Gassmann 方程后甚至可能生成虚数速度。解决临界孔隙度函数里增加if phi phi_c: return 0.0, 0.0把这段孔隙度标记为骨架失效。如果实际目标储层孔隙度确实超过临界孔隙度说明应该换用更复杂的未固结砂岩模型比如接触理论模型而不是继续死磕临界孔隙度公式。5.4 现象气体参数用的是地面条件导致流体替换失真原因天然气的体积模量和密度对温度压力极度敏感地面条件近似 0.0001 GPa 与储层条件下 0.03-0.08 GPa 相差两个数量级。用错参数后气层的 Gassmann 流体响应会非常弱甚至出现含气层和水层的 Vp 几乎一致的现象。解决储层流体参数一定要输入温压条件用 Batzle-Wang 公式或 PVT 实验数据计算。如果项目早期没有 PVT 资料我一般至少会做一次敏感性分析气体体积模量取 0.02、0.05、0.10 GPa 三个值看目标储层泊松比变化是否超过 0.02。如果超过说明该项目不能靠简化参数必须补 PVT 数据。5.5 现象速度、密度、孔隙度曲线深度不对齐原因声波测井、密度测井和中子孔隙度测井来自不同仪器本身就有深度误差和采样率差异。没有做深度匹配和重采样就直接进入流体替换会在地层边界上产生抛物线假象特别是薄储层中Vp 可能来自相邻泥岩层孔隙度却来自储层。解决流体替换前先做多曲线深度对齐再统一重采样到同一采样率。做完这一步后用交会图质量控制孔隙度和密度交会应该沿干净趋势分布出现明显离群点说明深度匹配还没到位。6. 用实际井点做一次标定反推临界孔隙度并检验流体替换计算流程跑通之后不要直接信任默认参数。我一般会用一口有可靠测井解释的目标井做标定已知含水饱和度、孔隙度、矿物体积分数以及实测 Vp、Vs然后反推干骨架模型里的关键参数比如临界孔隙度φ_c。目标函数很简单用选定的φ_c做流体替换算出的 Vs 或 Vp/Vs 与实测曲线最接近的那个值就是本区比较可信的骨架标定参数。def calibrate_phi_c(vp_obs, vs_obs, phi, k_min, mu_min, k_fluid, rho_min, rho_fluid): # vp_obs/vs_obs: 实测纵波/横波速度km/s # phi: 孔隙度数组rho_fluid: 目标流体密度 best_phi_c 0.4 best_error 1e9 for phi_c in [0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.45]: total_error 0.0 for i in range(len(phi)): k_dry, mu_dry dry_moduli_critical(k_min, mu_min, phi[i], phi_c) result fluid_substitution(k_dry, mu_dry, k_min, mu_min, k_fluid, phi[i], rho_min, rho_fluid) vp_calc, vs_calc result[3], result[4] total_error (vp_calc - vp_obs[i]) ** 2 total_error ((vs_calc - vs_obs[i]) * 2) ** 2 if total_error best_error: best_error total_error best_phi_c phi_c return best_phi_c, best_error代码里给 Vs 误差乘了权重 2因为 Vs 对流体不敏感能更好约束骨架参数如果只拟合 Vp很可能得到多个等价的φ_c。标定结果出来后再看替换后泊松比曲线和实测泊松比曲线的残差残差小于 0.02 说明参数组合可靠残差系统偏移则说明矿物模量或流体模量仍然有偏差而不是骨架模型的问题。这套标定方法不需要额外软件适合在储层评价早期把干骨架模型固定下来。多年做岩石物理项目的习惯让我坚持一点任何弹性参数计算结果都必须先接受实际测井曲线的交叉验证否则参数再合理也只是“看起来合理”而真正能经受井点检验的计算链后期做流体预测时才会少返工。希望这套从干燥岩石到流体饱和岩石的计算流程能帮你在自己的数据上少踩几个坑。本文还有配套的精品资源点击获取