企业官网建设流程全解析

简介：基于离散元方法（DEM）的岩石物理建模MATLAB代码包，面向地质力学、油气开发、采矿工程等相关领域的研究人员与工程师。代码利用颗粒离散元模拟，分析岩石在应力、流体饱和度等条件下的弹性模量、波速变化等力学特性，可用于地震响应解释和岩土工程问题预测。压缩包内共8个文件，全部为.m格式脚本，总大小仅5KB，核心脚本搭建DEM主框架，其余脚本分别针对拉梅参数、弹性模量、P波/S波速度以及含水饱和度与速度关系等物性参数进行建模，结构简洁、便于修改和功能扩展，也适合教学演示和二次开发。已有652人学习/下载，作者自编代码、思路直观，适合快速搭建数值实验；通过组合不同脚本，可获得颗粒尺度相互作用对宏观物性的影响规律，为储层评价、地基处理等场景提供参考，使用中遇到问题可在下载页留言交流，有助于缩短理解颗粒材料建模的入门时间。

1. 岩石物理 DEM 在解决什么问题：一个关于孔隙形状的真实例子

同样是 20% 孔隙度的灰岩，溶孔为主时纵波速度能到 5500 m/s 以上，裂缝主导时可能只有 4000 m/s 出头。孔隙度一样，速度差这么多，光靠孔隙度-速度经验公式解释不了。岩石物理 DEM（微分等效介质理论）就是为解决这类问题来的：它把孔隙形状作为输入，通过向矿物基质中逐步加入孔隙的方式计算干岩石骨架的弹性模量，再用流体替换公式得到饱和岩石的速度。这个方向适合做测井岩石物理解释、叠前反演和储层建模的从业者，也适合那些发现孔隙度-速度交会图散成一团、想搞清楚散点来源的研究生。下面从方程、代码、参数标定到常见翻车点，按可以直接落地的路径往下走。

2. DEM 的物理图像与核心方程：为什么要“微分”地加入孔隙

2.1 从 Voigt-Reuss-Hill 平均值说起：体积分数不够用

传统等效介质做法是拿矿物组分做平均。Voigt 上限假设所有矿物承受相同应变，Reuss 下限假设所有矿物承受相同应力，Hill 平均把两者取算术平均。问题是这类平均只吃“体积分数”，不吃“孔隙形状”。同一块 15% 孔隙度的砂岩，粒间孔为主和含裂缝为主，用 VRH 算出来的骨架模量完全一样，这显然不符合观测。

真实岩石里，孔隙不是球形气孔，而是从粒间孔、铸模孔到微裂缝的一系列旋转椭球体。椭球越扁，同样的孔隙体积对骨架的“削软”能力越强。要表达这种差异，模型里就必须多一个几何参数：孔隙纵横比 α，也就是椭球短轴与长轴的比值。有了这个参数，等效介质理论才能解释“孔隙度相同、速度不同”的散点带。

DEM 的微分思想很直接：不要一次性把全部孔隙丢进矿物里，而是把总孔隙度切成成千上万个小增量，每一次只往“当前等效介质”中加入极小一份孔隙。加完一小份后，把新的等效模量当作下一次迭代的宿主介质，继续加下一份。

写成标准 DEM 方程组就是下面这两个递推关系：

∂K*/∂φ = (Kf − K*) / (1 − φ) × P(α, ν*)

∂G*/∂φ = (Gf − G*) / (1 − φ) × Q(α, ν*)

其中 K* 和 G* 是当前等效介质的体积模量与剪切模量，Kf 和 Gf 是孔隙填充物的模量，φ 是从 0 开始累计的孔隙度。P 和 Q 是形状因子，依赖当前孔隙纵横比 α 和宿主介质的泊松比 ν*。每迭代一步，ν* 都要用最新的 K*、G* 重新算一次，这就是“宿主不断更新”的含义。

这种逐步迭代与 Kuster-Toksöz 一类的低浓度近似有本质区别。Kuster-Toksöz 假设夹杂物稀疏、互不影响，只适合孔隙度很低的场景；DEM 把孔隙间的相互作用通过宿主模量的递进演化“隐式”吸收进去了，所以在中高孔隙度下仍然可用。和自洽近似（SCA）相比，DEM 有明确的宿主概念——矿物基质是被替换的连续相，孔隙是后来加入的夹杂，这个主客关系更贴近成岩过程。代价是结果对“添加顺序”敏感，这点放到第 4 章重点说。

2.2 孔隙纵横比 α：决定骨架软硬的关键参数

孔隙纵横比定义为旋转椭球短轴/长轴。α=1 是球形孔，对应溶孔或粒间孔的理想形态；α=0.2 左右是偏扁的粒间孔；α=0.01 到 0.05 接近微裂缝、粒缘缝和黏土收缩缝。实际岩样里通常同时存在多种形态，所以工程上都用“多孔组”来描述：每个孔组有一个 α 和对应的体积占比。

获取 α 的第一手资料是铸体薄片和微 CT。薄片统计能看到粒间孔、溶孔的形态范围，微 CT 可以给出三维的孔喉分布，再拟合成若干个旋转椭球。没有岩心时，就只能用测井或实验室速度做反演，也就是第 5 章的工作流。需要特别强调的是，α 对骨架模量的影响远超直觉。固定孔隙度 15%，α 从 1 降到 0.05，干岩石体积模量可能下降一半以上，纵波速度差接近 20%。这解释了为什么裂缝发育段的声波速度会“莫名其妙”低一截。

α 和孔隙度 φ 是相互独立的两个输入。φ 控制孔隙“有多少”，α 控制孔隙“有多软”。建模时不要试图用一个等效 α 去吸收所有误差——那样拟合出来的 α 没有物理意义，换个井就失效。第 4 章的坑 5 会专门展开这个问题。

2.3 DEM 与 Gassmann 公式的分工：干骨架归 DEM，流体归经典公式

DEM 的自热场景是计算干岩石骨架。干岩样不是任何时候都能测到：取心不全、洗油不彻底、实验室测试频率与测井不一致，都会让直接测量的干骨架不可靠。这时用 DEM 从矿物模量和孔隙形态出发算 K_dry、G_dry，是一条稳定路径。

得到干骨架后，加流体要用经典流体替换公式，也就是常说的 Gassmann 公式：

K_sat = K_dry + (1 − K_dry/K_m)² / (φ/K_f + (1−φ)/K_m − K_dry/K_m²)

G_sat = G_dry

K_m 是矿物基质体积模量，K_f 是孔隙流体体积模量。这个公式的假设是低频、孔隙流体不承担剪切，测井频率和实验室低频测试基本满足。对于超声波兆赫兹级别的测量，频散效应会让实测 K_sat 偏高，标定时要单独处理。

注意不要把已经算完的干骨架再拿去做一次 DEM 迭代来“加流体”。原因有三：一是 DEM 对路径依赖，流体步骤会改变整个迭代历史；二是流体在孔隙中的力学响应是“整体压力支撑”，不是“夹杂弹性体”，形状因子 P/Q 并不描述这种机制；三是 Gassmann 公式本身就建立在孔隙几何无关的假设上，既然用了 Gassmann，就不需要再引入 α 对流体相的影响。换句话说，DEM 管骨架形态，Gassmann 管流体替换，两者是上下游关系，不是叠加关系。

2.4 输入参数准备：矿物模量与流体模量的常用参考区间

DEM 的输入参数并不多，但每一个都对结果有直接影响。矿物模量决定初始宿主，孔隙纵横比分布决定干骨架走势，流体模量决定饱和速度的流体项敏感性。下面是建模时常用的参考区间，具体项目里应以岩心 XRD 和实测弹性测试为准：

其中黏土矿物是最容易踩坑的输入。伊利石、蒙脱石、高岭石的模量差异很大，混在泥质砂岩里时，直接用某个“默认黏土值”会让基质模量偏高一截。我一般会先做敏感性测试：黏土模量上下浮动 20%，看最终 Vp 预测差多少。如果差得超过实测残差，说明这个输入必须依赖实测，不能拍脑袋。

3. 用 Python 实现 DEM 干骨架计算：可运行代码与三个关键参数

3.1 第一步：用 VRH 平均准备矿物基质模量

DEM 的宿主必须是“除孔隙以外的全部固相”的等效模量，而不是某一种矿物的模量。常见做法是先做 Voigt-Reuss-Hill 平均。下面这段函数接受矿物模量列表和体积分数，返回 Hill 平均：

def vrh_average(moduli, fractions): # moduli: [(K1, G1), (K2, G2), ...] 与 fractions 一一对应 kv = sum(f * k for (k, g), f in zip(moduli, fractions)) kr = 1.0 / sum(f / k for (k, g), f in zip(moduli, fractions)) gv = sum(f * g for (k, g), f in zip(moduli, fractions)) gr = 1.0 / sum(f / g for (k, g), f in zip(moduli, fractions)) return 0.5 * (kv + kr), 0.5 * (gv + gr)

逻辑上，Voigt 上限是体积加权算术平均，Reuss 下限是体积加权调和平均，Hill 取两者中值。这段代码的用途是把石英、方解石、黏土等多个固相合并成一个等效基质。如果你处理的是白云岩这类致密岩石，可以再把“只用 Reuss 下限”作为对比方案跑一遍，看宿主选取对结果的影响幅度。

3.2 最小可运行 DEM：离散增量迭代实现

下面是完整的 DEM 干骨架计算函数。核心思路是把总孔隙度切成 n_steps 个小份，按孔组顺序逐份加入，每加入一份就用最新的 K、G 更新泊松比和形状因子：

import numpy as np def shape_factors(alpha, nu): # 球形孔隙（alpha=1）的 Eshelby 形状因子解析式 P0 = (1.0 + nu) / (3.0 * (1.0 - nu)) Q0 = 2.0 * (4.0 - 5.0 * nu) / (15.0 * (1.0 - nu)) if alpha >= 1.0: return P0, Q0 # 扁球孔隙的工程代理：越扁，P/Q 越小，骨架越软 # 生产项目中请用完整 Eshelby 解或查表替换此代理 soften = 1.0 / (1.0 + 2.0 * (1.0 - alpha) / alpha) return P0 * soften, Q0 * soften def dem_dry(Km, Gm, Kf, Gf, phi_total, alpha_list, frac_list, n_steps=2000): K, G = Km, Gm phi = 0.0 dphi = phi_total / n_steps # 按体积占比分配步数，末组补齐，保证总步数严格等于 n_steps steps = [max(1, int(round(f * n_steps))) for f in frac_list] steps[-1] += n_steps - sum(steps) # 按 alpha_list 顺序逐孔组添加：顺序会影响结果，固定“先硬后软” for alpha, n in zip(alpha_list, steps): for _ in range(n): phi += dphi nu = (3.0 * K - 2.0 * G) / (2.0 * (3.0 * K + G)) P, Q = shape_factors(alpha, nu) K += (Kf - K) * P * dphi / (1.0 - phi) G += (Gf - G) * Q * dphi / (1.0 - phi) return K, G

逻辑说明：每一步迭代的宿主都是上一步更新后的等效介质，因此孔隙度增大到 20%、30% 时，模型仍然能跑，不会出现低浓度近似那种发散问题。泊松比 ν* 必须用最新 K、G 重算，因为形状因子对宿主的力学状态敏感，这一步做错会让 α 的影响被低估。参数说明：n_steps 是精度与速度的折中，建议不低于 1000，孔隙度超过 25% 时我一般加到 3000；alpha_list 与 frac_list 长度必须一致，frac_list 之和应为 1，函数内部会自动补齐末组步数；Kf、Gf 是孔隙流体模量，干骨架计算时 Kf 取气体或实验环境流体值，Gf 恒为 0。

注意：shape_factors 里的 soften 是一个工程代理形式，用来演示 DEM 的趋势行为。做定量标定时，请把这里替换成完整 Eshelby 解或查表得到的扁球 P/Q 系数，否则不同 α 之间的数值间距会有偏差。

3.3 用趋势曲线验证：α 对速度的影响一眼可见

跑完 DEM 后，最直接的验证方式是把干骨架模量换算成纵波速度，画出不同 α 下的孔隙度-速度趋势带：

Km, Gm = vrh_average([(37.0, 44.0), (76.8, 32.0)], [0.9, 0.1]) phi_list = np.linspace(0.02, 0.25, 12) for alpha in [1.0, 0.2, 0.03]: vp_trend = [] for phi in phi_list: K_dry, G_dry = dem_dry(Km, Gm, 2.3, 0.0, phi, alpha_list=[alpha], frac_list=[1.0]) rho = 2.65 * (1.0 - phi) + 1.05 * phi Vp = np.sqrt((K_dry + 4.0/3.0 * G_dry) / rho) vp_trend.append(Vp) print("alpha =", alpha, np.round(vp_trend, 2))

密度这里用基质密度 2.65 g/cc 与流体密度 1.05 g/cc 做简单加权，正式项目里用测井密度曲线更准。运行后你会发现三条趋势带在孔隙度 15% 时已经拉开明显差距，α=1 的趋势在最上方，α=0.03 在最下方。这就是 DEM 对“孔隙形状”表达能力的最直观验证。

实测数据对拍的标准动作是：把岩心或测井段的 Vp、Vs 投到同一张孔隙度-速度交会图上，看数据点落在哪两条趋势带之间。如果数据点整体高于 α=1 的趋势，往往不是模型错，而是基质矿物模量取低了，或泥质含量没扣除干净。

3.4 三个容易被忽略的边界条件：频散、压力与残余流体

DEM 参数标定不是纯粹的数学拟合，数据本身的物理条件必须对齐。第一是频散：超声波岩心测试在兆赫兹量级，测井在几十千赫兹，两者测到的模量本身就有差异，标定 DEM 时不能把超声波数据直接与测井速度混在一个目标函数里。第二是压力：有效压力升高会让微裂缝闭合，软孔的有效 α 从 0.01 升到 0.05 甚至更大。所以“同一块岩心、不同围压标定出的 α 分布不同”不是程序 bug，而是岩石的真实响应。第三是残余流体：干岩样洗油不彻底时，孔隙里残留的水会让 Kf 偏大，掩盖软孔的削软作用。判断方法很简单——如果 G_dry 异常偏高，基本可以怀疑含水没抽干净。

4. DEM 建模避坑清单：五个让结果翻车的习惯

4.1 坑 1：孔隙添加顺序不同，干骨架结果能差 10%

现象：同一套 α 分布，先加硬孔再加软孔与先软后硬算出来的 K_dry 相差 5%~10%，孔隙度越大差得越多。

原因：DEM 方程组本身是路径依赖的。每一步都把上一步的等效介质当宿主，所以“软孔出现时宿主多硬”会影响软孔对最终模量的贡献。先加软孔时，软孔在较硬的纯矿物宿主里被“放大”，削软作用更强；先加硬孔时，宿主已经被削弱，软孔的相对影响反而变小。

解决：固定一个添加顺序并写进建模说明里，我习惯用“先硬后软”，也就是 α 从大到小。同时做一次顺序敏感性测试：如果两种顺序的结果差异超过你的可接受误差，那这套 DEM 参数只能做定性解释，不能报定量结果。

4.2 坑 2：把孔隙度一次性代入，DEM 退化成线性近似

现象：算出来的 K_dry 随孔隙度接近直线下降，α 几乎不起作用，结果和低浓度近似差不多。

原因：迭代步数太少，或者直接把 φ_total 代进了微分方程的静态形式。“微分”两个字被丢了，孔隙间相互作用自然也没了。

解决：n_steps 至少取 1000，孔隙度超过 20% 时取 3000，然后做收敛性检查——把步数翻倍，看结果变化是否小于 0.5%。另外一个自检技巧是用 α=1 跑一组结果，和 Gassmann 干岩石极限趋势做对比，球形孔隙的 DEM 干骨架应该与 Gassmann 无机几何假设下的结果接近，如果差太多说明离散格式有问题。

4.3 坑 3：多矿物岩石直接用单一矿物当宿主

现象：碳酸盐岩建模时拿方解石模量直接当 Km，算出来 Vp 系统性偏高；泥质砂岩拿石英当宿主，低孔隙度段拟合很好，高孔隙度段越偏越大。

原因：宿主必须代表除孔隙以外所有固相的等效介质。方解石、白云石、黏土混合在一起时，物理上的“基质”已经不是纯方解石，单矿物模量会高估基质的刚性。

解决：先用 vrh_average 把全岩矿物合并成 Km、Gm，再把这个值作为 DEM 的初始宿主。更严格的做法是对“VRH 平均”和“Reuss 下限”两种宿主各跑一遍，把结果作为骨架模量的不确定性范围。

4.4 坑 4：干骨架和流体替换两头打架，速度预测在高孔隙度段失真

现象：低孔隙度段 Vp 预测很准，到高孔隙度段系统性偏低或偏高；或者流体替换后的速度与实测差得离谱。

原因：最常见的有两类。一是把流体当成夹杂又在 DEM 里迭代了一遍，流体项被路径依赖放大，等于算了两次“孔隙”；二是 Gassmann 公式里用了错误的骨架模量，比如拿饱和岩石的模量反推干骨架时单位没统一，或者 K_dry/K_m 算出来接近甚至超过 1，已经违背了 Gassmann 的基本约束。

解决：坚持“DEM 算干骨架，Gassmann 加流体”的单向流程。每次跑完检查 K_dry/K_m 是否在 0~1 区间内，这个比值超出正常范围（一般 0.3~0.8）时，先回头检查输入而不是调流体模量。

4.5 坑 5：单纵横比拟合就是玄学现场，换个井全部重来

现象：用单一 α 拟合 Vp 时残差很好看，但 Vs 或者泊松比趋势对不上；换一口井就得重新调 α，且调出来的值和薄片统计完全对不上。

原因：单一 α 是一个强非线性参数，它能够吸收矿物模量误差、孔隙度误差、甚至密度误差，导致 Vp 拟合“看起来成功”。但这种成功没有物理约束，一旦换数据就失效。更麻烦的是，Vp 对 α 和 φ 存在强耦合，单独用 Vp 约束时解不唯一。

解决：换成双孔组模型，固定硬孔 α1=0.2 左右，只反演软孔占比 f_soft 和软孔 α2，并且同时拟合 Vp、Vs 两条曲线。标定结果还要和薄片/CT 统计做交叉验证，如果反演出 60% 的软孔，而薄片上微孔只有 20%，那数学上再“最优”的结果也不可信。

5. 从岩心和测井标定 DEM 输入：双孔组软孔占比反演工作流

5.1 标定前的数据准备与质量检查

双孔组反演前，先把数据整理成一个深度对齐的表格，至少包含以下几个来源：

质量检查有三个重点。深度对齐是前提，孔隙度与速度曲线相差半米反演结果都会乱。压力要一致，超声波岩心数据必须在模拟地层压力的围压下测试，否则微裂缝状态与地下不同。泥质校正要做，黏土含量高的层段矿物模量要先算进 VRH，不能事后当误差处理。

5.2 固定硬孔形态，只反演软孔占比

双孔组模型的设定思路是：硬孔 α1 代表粒间孔和溶孔，形态相对稳定，取 0.2 是砂岩和碳酸盐岩里都比较稳的经验值；软孔 α2 代表微孔和裂缝，变化大，是反演的主角。真正需要拟合的未知量只有软孔体积占比 f_soft。

这样做的理由有两个。第一是降低病态性：同时反演 α1、α2、f_soft 三个参数时，目标函数面上存在一条近似的“深层谷”，很多组合都能给出同样好的 Vp 拟合。第二是与薄片可对应：薄片统计能直接给出粒间孔和微孔的占比，f_soft 反演结果可以和这个占比互相验证，多一层物理约束。

5.3 用实测 Vp、Vs 拟合软孔占比的最小脚本

下面这个脚本把 DEM 干骨架和 Gassmann 流体替换串成完整的前向模型，再用 scipy 优化软孔占比：

import numpy as np from scipy.optimize import minimize def forward(phi, f_soft, Vp_obs, Vs_obs, Km, Gm, Kf=2.3): # 双孔组：硬孔 alpha1=0.2，软孔 alpha2=0.03 alpha_list = [0.2, 0.03] frac_list = [1.0 - f_soft, f_soft] K_dry, G_dry = dem_dry(Km, Gm, Kf, 0.0, phi, alpha_list, frac_list, n_steps=2000) # Gassmann 流体替换，得到饱和岩石模量 K_sat = K_dry + (1.0 - K_dry / Km)**2 / ( phi / Kf + (1.0 - phi) / Km - K_dry / Km**2) G_sat = G_dry rho = 2.65 * (1.0 - phi) + 1.05 * phi Vp = np.sqrt((K_sat + 4.0/3.0 * G_sat) / rho) Vs = np.sqrt(G_sat / rho) return Vp, Vs def cost(f_soft, phi, Vp_obs, Vs_obs, Km, Gm): Vp_s, Vs_s = forward(phi, f_soft, Vp_obs, Vs_obs, Km, Gm) # 归一化残差，避免 Vs 量级小被 Vp 吃掉 return ((Vp_s - Vp_obs) / Vp_obs)**2 + ((Vs_s - Vs_obs) / Vs_obs)**2 # 示例：phi=0.18, Vp=3.2 km/s, Vs=1.9 km/s, 基质为纯石英示意 res = minimize(cost, x0=0.3, method='L-BFGS-B', bounds=[(0.05, 0.95)], args=(0.18, 3.2, 1.9, 37.0, 44.0)) print("f_soft =", res.x[0], "残差 =", res.fun)

逻辑说明：forward 先算干骨架，再用 Gassmann 把盐水效果加回去，最后换算 Vp、Vs。这样脚本既适用于实验室饱水岩样，也适用于测井条件。目标函数用归一化残差是因为 Vp 数值在 3~6 km/s，Vs 在 1.5~3.5 km/s，量级差异会掩盖 Vs 的约束作用。

参数说明：x0 是 f_soft 初值 0.3，表示先假设 30% 软孔；bounds 限制在 0.05~0.95，避免优化器跑到无物理意义的边界。Km、Gm 这里传了石英单矿物的示意值，正式计算一定要换成 vrh_average 的结果。Kf=2.3 GPa 是地层盐水体积模量的典型值，实际按温度压力查表。

这里有一个重要的边界提醒：如果你的资料只有 Vp，没有 Vs，那这个反演基本不可做。单靠 Vp 时 f_soft 与 α2 几乎完全相关，能拟合出无数个同样好的解。遇到这种情况，我建议把 α2 固定为 0.03，只反演 f_soft，并明确告诉项目方这是“低约束标定”，结果只能做趋势参考。

5.4 标定完成后的合理性诊断表

反演不是输出一个 f_soft 就结束，下面这几项检查每次都要过一遍：

6. 进阶用法：把反演的纵横比当孔隙结构指示器用

标定出来的 f_soft 和 α2 不要只停留在“建模参数”这个层面，它本身就是储层评价的线索。最实用的做法是把 f_soft 沿深度画成曲线，和裂缝发育段、溶孔段做对比。软孔占比突然升高的层段，通常对应微裂缝密集段或泥质微孔发育段；反过来，f_soft 很低而孔隙度也不低的位置，往往是优质粒间孔或溶孔储层。这个指标比单看孔隙度更接近“孔隙结构”的真实描述。

第二个用法是把标定好的双孔组参数带到全区井上，生成一组“结构归一化”的骨架速度趋势带。之后新井的实测速度点投到这张图上，偏离趋势带的幅度和方向，直接解读为孔隙结构异常。我在某区块做过一次：全区 40 口井用同一套双孔组参数，只有一口井在目的层段整体低于 α=0.03 趋势带，后来成像测井证实是构造裂缝发育段。这就是把 DEM 从“模量预测工具”升级成了“结构诊断工具”。

第三个用法是和统计建模结合。用 DEM 生成不同 φ、不同 α 分布下的弹性模量组合，作为合成样本库，再喂给机器学习模型做孔隙类型分类。DEM 在这里的价值是提供物理约束，让合成样本落在合理的模量空间内，而不是随机采样一个宽泛的参数范围。这样训练出的分类器在外推时不容易跑飞。

我自己的习惯是每次跑 DEM 前先固定“先硬后软”的添加顺序，并做一次顺序敏感性测试；如果换顺序结论就倒过来，那这次标定只能定性，不能定量。有一年我做过一个项目，Vp 拟合残差只有 1.5%，但 Vs 系统偏高 8%，查到最后是黏土模量取了石英参数。这种错怪不了 DEM，怪输入。把矿物模量、孔隙度、流体模量这三类输入按实测校准清楚，DEM 给出的干骨架和流体替换结果才真正能当作解释依据。希望帮到你。

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