MATLAB岩石物理计算：干燥与饱和弹性参数及Gassmann方程实现-北京尧图网络

简介这份资源面向地质勘探、地球物理及数值计算方向的学习者与研究人员围绕地震岩石物理基本理论用MATLAB实现干燥岩石与流体饱和岩石的弹性参数计算涵盖纵波模量、剪切模量和泊松比三类核心指标帮助读者理解岩石在压力与剪切力作用下的变形特性及其对地震波传播的影响。压缩包为rar格式仅含1个m脚本文件体积约2KB属于轻量级算法源码便于直接阅读、调试与二次开发。资源已有545人学习下载说明其在相关课程与科研入门中具有一定参考价值。通过该脚本读者可掌握弹性参数的数学建模思路与MATLAB编程实现方法并可将计算任务迁移至云计算环境进行批量模拟为储层评估、地震活动预测等场景提供基础工具支撑适合具备一定岩石物理与编程基础的中级学习者使用。1. 从一组干燥砂岩的模量算错说起这份 MATLAB 资源到底解决什么前阵子帮一个做储层预测的朋友看数据他拿一组干燥砂岩的纵波模量去反演流体因子结果饱和后的泊松比算出来是负的。翻回代码一看他把 Gassmann 方程里的干岩石骨架模量直接填成了矿物模量流体项整个串了。这类错误在地震岩石物理入门阶段特别常见——公式背得滚瓜烂熟一到干燥、饱和两种状态的参数传递就翻车。这份资源就是冲着这个痛点来的一个cementedsand.m脚本配合一份压缩包把干燥岩石和流体饱和岩石的纵波模量、剪切模量、泊松比从头算一遍不依赖任何商业岩石物理软件MATLAB 环境里跑通即可。适合刚接触岩石物理数值计算的学生、做叠前反演需要正演模量的工程师以及想把云计算资源接进 MATLAB 做批量模拟的开发者。它不教你地震波动力学全套只把弹性参数这条线拉直。2. 弹性参数的计算链路从矿物模量到饱和岩石的 Gassmann 替换2.1 干燥岩石的纵波模量与剪切模量怎么来的干燥岩石的弹性参数不是直接测出来的通常是从矿物模量、孔隙度和孔隙纵横比通过等效介质模型推出来的。这份资源里cementedsand.m走的是胶结砂岩模型的路子核心思路是先算干岩石骨架的体积模量 ( K_{dry} ) 和剪切模量 ( G_{dry} )再往里面填流体。干岩石体积模量的常见形式是[ K_{dry} K_{mineral} \cdot (1 - \phi)^{a} ]其中 ( \phi ) 是孔隙度( a ) 是胶结指数跟颗粒接触方式有关。剪切模量类似但指数项通常不同。脚本里把这两个指数做成了可调参数方便你按不同固结程度去试。实际操作时你需要在脚本开头把矿物模量、孔隙度、胶结指数这几个量填进去。常见做法是查文献取石英或方解石的 ( K ) 和 ( G )孔隙度从测井曲线读胶结指数先按 1.5 到 2.5 之间试。脚本里对应的变量名是K_min、G_min、phi、a_cement改完直接运行就能得到干燥状态的 ( K_{dry} ) 和 ( G_{dry} )。提示胶结指数对结果影响很大同一组孔隙度下指数从 1.5 变到 2.5纵波模量能差出 30% 以上。建议先固定其他参数单独扫这个指数看敏感度。2.2 流体饱和后的 Gassmann 方程为什么不能直接替换干燥岩石算完下一步是饱和流体。Gassmann 方程是低频极限下的标准做法它把干岩石模量、矿物模量、流体模量和孔隙度串起来[ K_{sat} K_{dry} \frac{(1 - K_{dry}/K_{mineral})^2}{\phi/K_{fluid} (1-\phi)/K_{mineral} - K_{dry}/K_{mineral}^2} ]剪切模量在 Gassmann 假设下不变( G_{sat} G_{dry} )。泊松比则从 ( K_{sat} ) 和 ( G_{sat} ) 反算[ \nu \frac{3K - 2G}{2(3K G)} ]脚本里把这段写成了一个独立函数gassmann_sat输入干岩石模量、矿物模量、流体模量和孔隙度输出饱和模量。流体模量根据你是水、油还是气查表或按经验公式给。常见做法是水取 2.2 GPa油取 1.0 GPa 左右气取 0.1 GPa 量级。这里有个容易忽略的点Gassmann 方程要求流体压力均衡所以只适用于低频。如果你拿它去算超声波频段的数据饱和模量会偏高泊松比也会偏。脚本里没有做频散校正这是它的边界不是 bug。2.3 用 MATLAB 把整条链路串起来脚本的结构很直白参数区、干岩石计算区、Gassmann 替换区、泊松比输出区。下面这段是核心逻辑的简化版你可以对照自己的cementedsand.m看% 输入参数 K_min 36.6; % 石英矿物体积模量, GPa G_min 45.0; % 石英矿物剪切模量, GPa phi 0.15; % 孔隙度 a_cement 2.0; % 胶结指数 K_fluid 2.2; % 水体积模量, GPa % 干燥岩石模量 K_dry K_min * (1 - phi)^a_cement; G_dry G_min * (1 - phi)^a_cement; % Gassmann 饱和替换 K_sat K_dry (1 - K_dry/K_min)^2 / ... (phi/K_fluid (1-phi)/K_min - K_dry/K_min^2); G_sat G_dry; % 泊松比 nu_dry (3*K_dry - 2*G_dry) / (2*(3*K_dry G_dry)); nu_sat (3*K_sat - 2*G_sat) / (2*(3*K_sat G_sat)); fprintf(干燥: K%.2f GPa, G%.2f GPa, nu%.3f\n, K_dry, G_dry, nu_dry); fprintf(饱和: K%.2f GPa, G%.2f GPa, nu%.3f\n, K_sat, G_sat, nu_sat);逻辑说明先按胶结模型算干岩石再走 Gassmann 得到饱和模量最后统一转泊松比。参数说明K_min和G_min按矿物类型查phi从测井或岩心数据来a_cement控制固结程度K_fluid决定流体类型。跑完对比干燥和饱和两行输出如果饱和泊松比反而比干燥低很多大概率是流体模量填错了。3. 把计算任务丢到云端MATLAB 并行与批量参数扫描3.1 为什么单机跑参数扫描会卡岩石物理正演本身不慢一次计算毫秒级。但如果你要做不确定性分析比如孔隙度从 0.05 到 0.30、胶结指数从 1.0 到 3.0、流体模量分水油气三档组合起来就是几千到几万次计算。单机 for 循环跑下来MATLAB 的解释器开销会让总时间从几秒变成几分钟甚至更久。更别说后面还要接反演每次迭代都调一遍正演。云计算在这里的价值不是“算得快”而是“同时算很多组”。把参数空间切成块每块丢一个 worker总时间就压下来了。这份资源本身不包含云平台配置但cementedsand.m的写法留了并行接口——参数区是独立的计算区没有全局依赖可以直接套parfor。3.2 用 parfor 做参数扫描的改法把原来的单点计算改成扫描核心是把参数组合展开成矩阵然后用parfor分发。下面是一个可抄的骨架% 参数网格 phi_vec 0.05:0.01:0.30; a_vec 1.0:0.1:3.0; fluid_vec [2.2, 1.0, 0.1]; % 水, 油, 气 % 展开所有组合 [PHI, A, FLUID] ndgrid(phi_vec, a_vec, fluid_vec); params [PHI(:), A(:), FLUID(:)]; n size(params, 1); % 预分配输出 K_dry_out zeros(n,1); K_sat_out zeros(n,1); nu_sat_out zeros(n,1); % 并行扫描 parfor i 1:n phi params(i,1); a params(i,2); Kf params(i,3); Kd K_min * (1 - phi)^a; Gd G_min * (1 - phi)^a; Ks Kd (1 - Kd/K_min)^2 / ... (phi/Kf (1-phi)/K_min - Kd/K_min^2); Gs Gd; K_dry_out(i) Kd; K_sat_out(i) Ks; nu_sat_out(i) (3*Ks - 2*Gs) / (2*(3*Ks Gs)); end逻辑说明ndgrid把三个参数展成全组合parfor把循环体分发到多个 worker。参数说明phi_vec和a_vec的步长决定扫描密度fluid_vec对应不同流体。注意parfor里不能有依赖上一次迭代的变量所以所有输出都预分配成数组按索引写。注意parfor启动并行池有开销如果总计算量不到几百次单机 for 循环反而更快。一般组合数超过 1000 再考虑并行。3.3 云端部署时容易忽略的路径与许可问题把 MATLAB 放到云服务器上跑第一道坎不是代码是许可。MATLAB 的并行计算工具箱需要额外的 license云厂商的镜像里不一定带。常见做法是先用ver命令确认Parallel Computing Toolbox在不在不在的话要么换镜像要么改用parfeval配合本地 worker。第二道坎是路径。cementedsand.m里如果写了相对路径读数据云端工作目录和本地不一样会报文件找不到。稳妥的做法是在脚本开头用fullfile拼绝对路径或者把数据文件跟脚本放同一目录用mfilename(fullpath)取脚本位置再拼。第三道坎是结果回传。批量扫描的输出可能有几十万行直接save成 mat 文件再下载比在云端做可视化更省事。如果非要画图建议只回传抽样后的散点数据别把整个网格传回来。4. 避坑与排查弹性参数计算里最容易翻车的五个点4.1 饱和泊松比出现负值现象跑完脚本饱和岩石的泊松比小于 0甚至到 -0.5。原因Gassmann 方程里 ( K_{sat} ) 算出来比 ( K_{dry} ) 还小或者流体模量填成了负数。最常见的是把K_fluid和K_min的变量名写反了或者孔隙度填成了百分数15 而不是 0.15。解决先检查phi是否在 0 到 1 之间再检查K_fluid是否为正。如果都正常打印K_dry和K_sat对比正常情况 ( K_{sat} K_{dry} )。如果反了说明 Gassmann 分母里的符号错了对照公式逐项核对。4.2 干燥和饱和剪切模量不一致现象脚本输出里 ( G_{dry} ) 和 ( G_{sat} ) 不相等。原因Gassmann 理论假设剪切模量不受流体影响所以 ( G_{sat} ) 必须等于 ( G_{dry} )。如果代码里给 ( G_{sat} ) 单独套了一个公式或者并行扫描时索引串了就会出现不一致。解决在脚本里强制G_sat G_dry;不要给剪切模量做任何流体替换。如果用了并行检查G_sat_out(i)的赋值是否误写成了G_dry_out(i1)之类的越界索引。4.3 胶结指数取值的玄学现象同一组孔隙度别人算出来的纵波模量跟你差一倍。原因胶结指数没有唯一标准不同文献给的取值范围从 1.0 到 3.0 都有。取 1.0 接近临界孔隙度模型取 3.0 接近高固结砂岩。如果没说明取值依据结果没法对比。解决在脚本注释里写清楚a_cement的来源比如“按某地区砂岩岩心拟合值取 1.8”。做敏感度分析时把a_cement作为扫描变量输出模量随它的变化曲线而不是只报一个数。4.4 parfor 里变量切片报错现象parfor循环报“无法分类变量”或“无效索引”。原因parfor要求循环内的变量要么是广播变量不随迭代变要么是切片变量按索引写。如果代码里出现了K_dry_out [K_dry_out; Kd];这种动态增长就会报错。解决所有输出数组提前用zeros(n,1)预分配循环内只写K_dry_out(i) Kd;。如果需要在循环内累加改用parfeval或把累加拆成两阶段。4.5 云端跑完结果对不上本地现象同样的参数本地跑出来 ( \nu_{sat} 0.28 )云端跑出来 0.31。原因最常见的是 MATLAB 版本不同parfor的归约顺序有差异导致浮点累加结果有微小偏差。如果偏差大到 0.03那更可能是云端脚本读到了旧版参数文件或者并行池里某个 worker 的路径指向了另一个副本。解决在脚本开头打印version和pwd确认云端和本地一致。批量扫描时把参数矩阵和结果一起存成 mat 文件回本地用isequal对比关键列。如果只是浮点末位差异可以忽略如果量级不对先查文件路径。5. 进阶技巧用干燥模量反推流体因子与敏感度矩阵5.1 从饱和模量反推流体模量的迭代写法正演是已知流体算饱和模量反演是已知饱和模量反推流体。Gassmann 方程对 ( K_{fluid} ) 没有解析解得用数值迭代。常见做法是二分法或牛顿迭代把 ( K_{fluid} ) 限制在 0.01 到 3.0 GPa 之间。function Kf invert_fluid(K_sat_obs, K_dry, K_min, phi) % 二分法反推流体模量 lo 0.01; hi 3.0; for iter 1:50 mid (lo hi) / 2; Ks K_dry (1 - K_dry/K_min)^2 / ... (phi/mid (1-phi)/K_min - K_dry/K_min^2); if Ks K_sat_obs lo mid; else hi mid; end end Kf (lo hi) / 2; end逻辑说明二分法每次把区间折半50 次迭代后精度足够。参数说明K_sat_obs是观测到的饱和模量K_dry和K_min从干燥计算和矿物查表来phi是孔隙度。这个函数可以嵌在parfor里对每个观测点单独反推。5.2 敏感度矩阵哪个参数对泊松比影响最大做正演的人经常被问“泊松比变了到底是孔隙度变了还是流体变了”。回答这个问题需要敏感度矩阵也就是 ( \partial \nu / \partial \phi )、( \partial \nu / \partial K_{fluid} ) 这些偏导。不用解析求导数值差分就行参数扰动方式泊松比变化量典型值孔隙度0.15 → 0.160.008流体模量2.2 → 1.0 GPa-0.015胶结指数2.0 → 2.1-0.003矿物模量36.6 → 37.0 GPa0.001这张表是用脚本跑出来的每行只改一个参数其他固定。从量级看流体模量的影响最大孔隙度次之胶结指数再次矿物模量最小。这意味着如果你要反演流体泊松比是比纵波模量更敏感的量。5.3 我自己的习惯每次改参数前先存一份基准从那以后我每次动cementedsand.m的参数区之前都强制走一遍基准测试把默认参数跑一遍把K_dry、K_sat、nu_dry、nu_sat四个数存到一个baseline.mat里。改完参数再跑用isequal或max(abs(diff))对比确认只有目标参数在变其他没被误改。这个习惯帮我逮住过好几次变量名手滑——有一次把K_fluid改成了K_fluid_new结果脚本里另一处还在引用旧名跑出来饱和模量根本没变差点把错误结论写进报告。希望帮到你。本文还有配套的精品资源点击获取

MATLAB岩石物理计算：干燥与饱和弹性参数及Gassmann方程实现

相关推荐

Widar3.0实战：用Matlab实现WiFi CSI行为识别全流程

会议纪要哪个软件比较准确？2026年7款主流工具统一维度实测横评

MegaSR C105 RAID驱动兼容指南：Linux内核与固件匹配实战

更多资讯

LLC闭环仿真全流程：从控制策略到环路补偿的防发散实践

Cesium全球体积云渲染实战：光线步进与性能优化

跌倒检测数据集全解析：格式转换、清洗与训练避坑指南

Mistral AI 发布第二代开源编码模型：Devstral2 和 Devstral Small2 实测与 TaoToken 接入

Stata 17 模型代码版本兼容与回归结果复现指南

OpenWebUI部署实战：从Ollama到vLLM的模型对接与避坑指南

今日更新

金融AI落地指南：从调研报告到技术选型与工程实践

大模型上下文成本与注意力优化技术深度解析——为何O(n²)平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测

基于GMM的语音识别系统原理与MATLAB实现详解

本周精选

金融AI落地指南：从调研报告到技术选型与工程实践

大模型上下文成本与注意力优化技术深度解析——为何O(n²)平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测

基于GMM的语音识别系统原理与MATLAB实现详解

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？