★ FEATURED ARTICLE 简介面向石油工程、环境工程、材料科学等领域的科研人员与工程师这份基于Matlab的多孔介质孔隙网络建模软件包提供了一套从微观图像处理到流体传输模拟的完整工具链。包内共398个文件压缩后约10.11MB核心为54个m脚本配合196个vtk数据可视化文件与72个dat网络/节点数据文件可支持CT、SEM图像分割、孔隙网络自动生成、孔隙率/渗透率/比表面积计算以及油水气多相流动模拟另有html、md、rst等格式文档便于查阅使用说明。资源还包含图片、表格和MAT数据文件便于对照运行和结果分析。目前已有86人学习适合需要开展PNM建模或快速搭建Matlab仿真环境的中高级使用者。通过该软件包可跳过从零编写算法的大量工作直接获得可修改、可扩展的建模脚本与示例数据用于油气储层分析、地下水污染物传输研究及多孔材料性能预测等场景。1. 多孔介质孔隙网络建模Matlab 软件包到底解决什么问题拿到一张 micro-CT 扫描出来的岩心切片很多人第一反应是先把孔隙识别出来于是花一整周调分割阈值、做连通域分析最后画出一张花花绿绿的孔隙分布图。可一旦导师追问渗透率是多少、毛管压力曲线能不能出来这套流程就卡住了——因为图像分割只回答了孔隙长什么样回答不了流体在里面怎么流动。基于 Matlab 的多孔介质孔隙网络建模软件包做的正是这后半件事把三维体数据里的孔隙空间抽象成由孔隙体和喉道组成的网络再用这个网络去算渗透率、相对渗透率、毛管压力这些工程关心的量。它适合手里有 CT 或 FIB-SEM 图像、却没有精力从零写孔网提取算法的岩土、能源、材料方向从业者也适合想先跑通流程再决定要不要上商业软件的研究生。做这行的人有个共识孔网建模的成败七成不在图像分割而在把孔隙空间抽象成网络时你做了哪些几何假设。同一个体数据假设不同算出来的渗透率能差一个数量级。这套软件包的价值就是把那些容易翻车的抽象步骤封装成可复现的流程让使用者把精力放在参数标定和结果解释上。下面从原理讲到落地把这条路走一遍。2. 把孔隙空间变成网络孔网建模的核心原理与软件包工作流2.1 两种主流抽象算法最大球法与中轴线法的取舍孔隙网络建模的第一步是把连续的三维孔隙空间离散成一堆孔隙体和连接它们的喉道。常见做法有两大流派最大球法Maximal Ball和中轴线法Medial Axis。最大球法的思路很直白对每个孔隙体素找一个以它为中心、且完全落在孔隙内部的最大内切球这个球的半径就是该位置的局部孔径。把半径在局部邻域内最大的那些球挑出来作为孔隙体两个孔隙体之间最窄的通道就是喉道。Matlab 里做这一步通常用bwdist做距离变换再用imregionalmax找局部极大值。它的优点是对噪声相对不敏感实现直观适合形状不规则的孔隙空间缺点是计算量大而且最大球在细长喉道处容易给出偏大的半径估计。中轴线法则是先把孔隙空间细化成骨架线再沿骨架线扫描截面截面极小处定义为喉道。这种方法在几何上更贴合管道流的直觉对喉道定位更准但骨架线对表面噪声非常敏感——一个微小的凸起就可能让骨架长出伪分支后面还要花力气做剪枝。我一般会优先用最大球法跑第一版结果因为它的数值稳定性好只在喉道形态特别关键的使用场景比如要精确模拟水锁效应才会切到中轴线法做交叉验证。两种方法在 Matlab 软件包里的最终输出是一样的一组孔隙节点一组喉道连接每个元素都带几何属性。理解这一点很重要——不管底层算法是哪派下游的渗流模拟都只吃这组拓扑和几何所以选算法的本质是选几何误差分布。2.2 孔网描述里的关键几何量与拓扑量孔网抽取完每个孔隙体和喉道都会附带一串几何量这些量直接决定后续流动模拟的准确性。孔隙体这边最重要的是等效半径、体积和形状因子。等效半径通常取最大内切球半径体积则直接从体素积分算。形状因子Shape Factor是截面面积与周长的无量纲组合定义为 G A / P²其中 A 是喉道截面面积P 是截面周长。圆管的形状因子约 0.0796三角形截面在 0.02 到 0.05 之间狭缝则更低。这个量看着不起眼但它直接进到毛管压力公式里——用圆管假设算出来的毛管压力跟真实三角形孔道能差出几倍。拓扑量这边最重要的是配位数连通数即每个孔隙体连接了多少个喉道。真实岩心的配位数分布通常在 2 到 6 之间如果提取结果里大量孔隙配位数是 1 或 2基本可以断定提取参数太激进把喉道切断了不少。另一个常用拓扑量是孔喉比——孔隙半径与喉道半径之比它控制着捕集和滞后效应在模拟驱替过程时是必调的参数。软件包通常会把这堆量组织成一个pore结构体和一个throat结构体前者存半径、体积、坐标、形状因子后者存两端孔隙编号、半径、长度、形状因子。这套数据格式同时也是行业里交换孔网模型的通用格式掌握之后迁移到其他平台也顺。2.3 软件包的常见工作流从 TIFF 堆栈到可模拟的 PNM这类软件包的典型工作流是五步。第一步读入二值化的三维体数据通常是一组 TIFF 切片或单个三维 MAT 数组第二步孔隙体识别用距离变换加局部最大找种子点第三步喉道划分在相邻孔隙之间找收缩截面第四步简化与清理包括合并过度分割的孔隙、削掉孤立死端、重算配位数第五步输出孔网结构供下游的稳态渗流或驱替模拟调用。这套流程里最容易出问题的不是最后一步而是第二步和第四步之间的衔接——种子点挑多了后面孔喉比失真挑少了小孔隙被合并成大孔隙渗透率偏高。所以成熟的软件包都会在这一段暴露种子最小半径和局部最大邻域大小两个参数给你调而不是黑匣子一把梭。3. 在 Matlab 里把孔网提取跑通三个可直接上手的最小骨架3.1 加载体数据并清理孤立噪声TIFF 堆栈读取函数不管软件包封装得多高级入口一定是体数据加载。我通常先写一个加载函数把目录下的 TIFF 切片组装成三维逻辑数组顺手做掉两步脏活按文件名正确排序以及用最小连通域阈值清掉打在孔隙相里的孤立噪声点。function bw loadTiffStack(tiffDir, minVol) % 读取按文件名排序的 TIFF 切片组装成三维逻辑体 % tiffDir : 存放切片的目录路径 % minVol : 小于该体素数的连通域视为噪声置为固体 files dir(fullfile(tiffDir, *.tif)); nz numel(files); info imfinfo(fullfile(tiffDir, files(1).name)); bw false(info.Height, info.Width, nz); for k 1:nz bw(:,:,k) imread(fullfile(tiffDir, files(k).name)) 0; end % 删除孤立噪声团块并保留最大的孔隙连通域 bw bwareaopen(bw, minVol); CC bwconncomp(bw, 26); numPixels cellfun(numel, CC.PixelIdxList); [~, idx] max(numPixels); bw false(size(bw)); bw(CC.PixelIdxList{idx}) true; end这段代码的逻辑是先把每张切片按大于 0 的阈值二值化然后用bwareaopen清掉小于minVol体素的连通域最后只保留最大连通域。minVol的常见取值是 50 到 200 体素具体取决于你的体素尺寸——如果体素边长是 1 微米一个等效半径 3 微米的孔隙约 113 体素那minVol取 50 比较合适避免误删真实小孔。注意这里保留最大连通域是有代价的它会删掉与主孔隙体不连通的孤立孔洞如果研究对象本身就是不连通的孔比如某些封闭空腔这段要改写成按体积阈值过滤而不是只取最大域。做连通性分析时统一用 26 邻域否则斜向通道会被误判为断开。3.2 种子点检测与孔隙识别距离变换加局部最大这一步是孔网提取的心脏。先对孔隙相做距离变换让每个孔隙体素的值等于它到最近固体边界的距离再乘上体素尺寸换成物理长度。然后把局部最大点作为孔隙种子候选并用最小半径阈值过滤掉那些由噪声产生的微小种子。function poreSeeds findPoreSeeds(bw, voxelSize, minR) % 基于距离变换的局部最大找到孔隙种子 % bw : 三维逻辑数组true 为孔隙 % voxelSize : 体素边长单位微米 % minR : 种子最小等效半径单位微米 D bwdist(~bw) * voxelSize; % 距固体边界的距离场 % 26邻域局部最大半径小于 minR 的候选直接丢弃 seedMask imregionalmax(D) (D minR); % 分水岭约束只保留离固体边界足够远的种子 seedMask seedMask (D minR); S bwlabeln(seedMask, 26); props regionprops3(S, D, MaxIntensity, Centroid); keep props.MaxIntensity minR; poreSeeds.radius props.MaxIntensity(keep); poreSeeds.centroid props.Centroid(keep, :); poreSeeds.label find(keep); end核心参数是minR。取太小噪声和表面起伏会被当作孔隙体导致配位数虚高取太大小孔隙被合并渗透率偏高。我一般先取2 * voxelSize跑一版再看孔隙数分布曲线做调整。imregionalmax默认用 26 邻域对于细长孔道这个邻域设置会让相邻种子容易被合并必要时可改成 6 邻域配合minR细化。这段代码给的是找种子真正的喉道识别在简化工程实现里经常用球间相交判定两粒种子中心距小于两者半径之和乘以一个相交系数beta就认为它们之间有喉道。beta通常取 1.1 到 1.3越大喉道越多、配位数越高是后续渗透率标定的主要旋钮。3.3 组织 PNM 结构体为下游模拟留好标准字段提取完种子的下一步是把孔喉信息组装成统一结构体并落盘。下面这段代码完成两件事补算孔隙体积以及以-v7.3格式保存以支持超过 2GB 的大文件。function writePNM(pnm, outfile) % 补算缺失的孔隙体积并保存为 PNM 标准格式 % 输入 pnm 至少包含 pore.radius 与 throat.radius、throat.connectivity if ~isfield(pnm.pore, volume) pnm.pore.volume 4/3 * pi * pnm.pore.radius.^3; end % 打印关键诊断量保存前最后确认一遍网络合理性 fprintf(孔隙数: %d\n, numel(pnm.pore.radius)); fprintf(喉道数: %d\n, size(pnm.throat.connectivity, 1)); fprintf(平均配位数: %.2f\n, ... mean(histcounts(pnm.throat.connectivity(:), ... (1:max(pnm.throat.connectivity(:))1)))); % 统一元数据 pnm.meta.savedBy PNM pipeline; save(outfile, pnm, -v7.3); endthroat.connectivity是 N×2 数组每行存一条喉道两端的孔隙编号这是孔网数据格式的命脉下游模拟全靠它组装传导率矩阵。保存前打印平均配位数是个好习惯——如果平均值落在 1.5 以下说明网络已经碎得不像样别急着存盘先回头调参数。-v7.3用 HDF5 格式存储读取速度比旧格式快还能绕开 2GB 文件限制。注意histcounts统计配位数时要先把connectivity展平否则会按列统计出错误分布。4. 孔网建模的必调参数分辨率、种子阈值与标定旋钮4.1 体素分辨率与 REV 尺寸最少像素数决定几何可信度孔网提取对体素分辨率的敏感程度远超想象。经验法则是等效孔径至少要有 8 到 10 个体素覆盖否则距离变换给出的最大内切球半径会严重偏小。比如孔径 20 微米的孔隙体素边长超过 2 微米时提取出的半径可能只有真实值的七成随后渗透率模拟值跟着跳水。做这一步前先做一个表征单元体REV收敛测试取边长为 L 的子块分别算孔隙度、比表面积和平均配位数逐步增加 L 看结果是否收敛。一般来说L 达到平均孔径的 10 倍以上这三个量才会稳定。很多翻车案例都发生在拿一个边长不足 5 倍孔径的裁剪块直接建模结果换了裁剪位置结果就大变源头就是 REV 没过关。我习惯把 REV 测试写进流程的第一步先取 100³、150³、200³ 三个子块各跑一轮提取对比孔隙度和配位数均值三者偏差小于 3% 才继续往下走。这一步看起来费时间实际上能省掉后面所有结果解释不了的排查时间。4.2 种子最小半径与相交系数两个直接影响孔隙数量的旋钮软件包暴露的参数里minR种子最小半径和beta喉道相交系数是出现频率最高的一对。minR控制小孔隙的去留beta控制喉道的密度两者的失衡会制造出两种典型的坏结果孔隙数爆炸或网络过稀。调参顺序上先把minR定住再动beta。判断标准看平均配位数目标值参考真实岩心的 2 到 6如果你的样品是颗粒堆积物配位数可以到 5 以上如果是致密碳酸盐岩2 到 3 才算正常。配位数太高多半是beta偏大把不该相连的孔隙接上了配位数太低先怀疑minR偏大导致小孔隙被吞并。工程上还有个低成本校验办法把提取出的孔网按配位数上色叠回原始 CT 切片上目检。如果配位数高的孔隙集中在某一层或某条裂缝带上说明算法很可能是按虚假连通路径把它们接在一起的这时调整局部邻域设置比全局改参数更有效。4.3 用实验数据标定形状因子与收缩角参数表与标定顺序几何抽象参数定完之后最后一道标定是针对流动阻力的。两个最关键的标定量是形状因子 G 和收缩角corner angle。下表给出常见默认值和调参方向照着做至少不会让首轮模拟偏得离谱。参数常见默认影响对象典型调参依据形状因子 G圆管 0.0796毛管压力、传导率用压汞毛管压力曲线反标定收缩角30°~60°三角形截面相对渗透率、残余饱和度用油水驱替实验数据标定孔喉比2~5捕集与滞后用离心或驱替效率数据约束喉道长两孔隙中心距减去两半径绝对渗透率模拟渗透率对拍实验值标定的顺序有讲究先用绝对渗透率实验数据锁喉道长度和形状因子因为绝对渗透率只对这两者敏感确认绝对渗透率对上之后再用毛管压力曲线调收缩角。如果一上来就同时调四个参数很容易出现刚好凑出正确答案的假象——换一个边界条件就打回原形。这个先几何、后拓扑、再流动的标定顺序是让模型具备预测能力的底线。5. 孔网建模避坑指南五个高频坑位与排查办法5.1 提取出的孔隙数忽多忽少局部最大邻域设置的锅现象同一块数据改imregionalmax的邻域参数孔隙体数量从几千跳到几万配位数分布完全变样。原因默认 26 邻域的局部最大对细长孔道太敏感——孔道里每个弯曲段都可能产生一个局部极大值被当成独立孔隙体。噪声像素在距离变换场里也会形成虚假尖峰。解决先用imgaussfilt3对距离场做轻度高斯平滑sigma 取 1 到 2 个体素削掉表面噪声产生的尖峰再把邻域改成 6 邻域让种子更稀疏、更符合局部最大的物理直觉。改完之后重跑配位数分布正常应该落在 2~6 的窄区间。5.2 孔隙率与实验值差五个百分点分割阈值的边界效应现象用氦测孔隙度实验值是 22%软件包从二值体里统计出的孔隙率只有 17%差得整条模拟曲线都站不住。原因原始 CT 灰度图分割时阈值取得太保守把大量微喉道和边界上的半体素判成了固体另一个常见来源是载荷函数里bwareaopen把小于阈值的真实小孔团块当噪声清掉了而这些小孔合计能贡献几个百分点的孔隙率。解决先不改孔网参数回到二值化环节把阈值往孔隙方向偏移 5% 到 10%重算孔隙率直到和实验值对上同时把minVol降到 20 体素以下看孔隙率是否有明显抬升。如果这两个动作都不奏效再检查是不是 REV 选取本身就落在了非均质区。5.3 模拟渗透率比实验高一个数量级形状因子默认值太天真现象用圆管假设G0.0796模拟的绝对渗透率是实验值的 5~10 倍油藏工程里这是最经典的翻车现场。原因真实孔喉截面几乎从不是圆形的多为三角形或不规则狭缝。圆管是水力半径最大的截面形态流动阻力最小所以模拟渗透率高得离谱。形状因子这个量差 3 倍渗透率就差一个数量级。解决先用经验公式把形状因子从圆管值改到 0.02~0.05 区间跑一遍看渗透率是否回落到合理范围然后用压汞毛管压力曲线的进汞段做一次形状因子反演把每个喉道的 G 换成样品实测值而不是全网络统一一个常数。5.4 数组一上 512³ 就 Out of Memory浮点中间结果是元凶现象300³ 体数据能跑换 512³ 直接报内存错误Matlab 卡死只能强杀。原因imregionalmax输入的是 double 类型的距离场一个 512³ 的 double 数组就是 1GB再加上bwlabeln和regionprops3的中间变量峰值内存能冲到 6~8GB。你以为自己买了 16GB 内存机器就够其实中间过程吃掉了一大半。解决第一步把距离场转成single精度对于毫米级孔径建模完全够用第二步在bwdist之前把bw保持为 logical 类型不要让任何中间变量被自动提升为 double第三步如果数据更大按 z 方向分块做距离变换再拼接。省下内存的另一个好处是速度也上去了double 数组的访存开销在三维场上非常可观。5.5 换个重建分辨率结果面目全非REV 与体素尺寸没做收敛检查现象同一块样品用 1 微米和 2 微米体素分别重建提取出的孔隙数和绝对渗透率差出两倍不知道信哪个。原因体素尺寸超过了平均孔径的十分之一几何信息在大尺寸体素下被平均掉了顶部细小的孔结构直接消失。这不是算法 bug而是采样定理在三维空间里的硬约束。解决对同一块子块做分辨率敏感性扫描——把 1 微米体素的数据先用imresize3降采样到 2 微米和 3 微米各跑一遍完整流程看孔隙度和渗透率的变化曲线。如果渗透率随分辨率增大呈单调下降说明你的特征尺度没有充分采样要么换更高分辨率的扫描要么把目标孔径范围从分析中剔除。6. 从建模到模拟把 PNM 接到渗流计算上的进阶用法孔网建完之后最有价值的落地动作是直接算绝对渗透率。用孔网算渗透率的思路很简单每个喉道看成一根圆管用泊肃叶定律给每根管算传导率再把所有喉道组装成稀疏矩阵固定入口出口压力解一次线性系统得到总流量最后由达西定律反推渗透率。这整条流程在 Matlab 里只需要几十行代码。function K_abs computePermeability(pnm, mu, L_core, A_core) % 基于孔网结构组装传导率矩阵并计算绝对渗透率 % pnm.throat.radius : 喉道半径数组, pnm.throat.connectivity : Nx2 喉道两端孔隙编号 nPore numel(pnm.pore.radius); g pi * pnm.throat.radius.^4 ./ (8 * mu .* pnm.throat.length); % 圆管传导率 A sparse(nPore, nPore); for k 1:size(pnm.throat.connectivity, 1) i pnm.throat.connectivity(k, 1); j pnm.throat.connectivity(k, 2); A(i,i) A(i,i) g(k); A(j,j) A(j,j) g(k); A(i,j) A(i,j) - g(k); A(j,i) A(j,i) - g(k); end % 固定入口出口压力入口孔隙压力 1 Pa出口 0 Pa % 这里用移除方程法处理 Dirichlet 边界 inlet find(pnm.pore.centroid(:,3) 50); % 示例按 z 坐标划入口 outlet find(pnm.pore.centroid(:,3) 200); free setdiff(1:nPore, [inlet; outlet]); P zeros(nPore, 1); P(inlet) 1; P(outlet) 0; Q zeros(nPore, 1); Q(free) A(free, free) \ (-A(free, [inlet; outlet]) * P([inlet; outlet])); % 总流量 入口孔的总流入 q_total sum(A(inlet, :) * P); K_abs q_total * mu * L_core / (A_core * 1); % 达西定律压差为 1Pa end这段代码有两点要说明传导率公式用的是圆管假设实际使用应该换成第 4 章标定过的形状因子修正版本边界处理用了最简单的移除方程法Matlab 的稀疏矩阵反斜杠对几千阶的系统解起来毫秒级完全不需要迭代法。跑通这条计算链路之后孔网建模才真正闭环——你可以把模拟渗透率跟实验室岩心实测值对拍对不齐就回头调第 4 章那些参数而不是对着 CT 图空转。我现在的习惯是每次接一个新样品先花半天把 4.1 的 REV 收敛测试和 4.3 的标定顺序走一遍把几何参数和流动参数分开标定绝不混着调。这习惯救过我很多次——有回换了个致密样品直接沿用上一批松散砂岩的参数模拟渗透率偏高 30 倍排查了半天才发现是形状因子没重新标定。孔网建模这行的玄学感大多来自跳过这一步。把这条先收敛、后标定、再模拟的路径走顺这套 Matlab 软件包就能从玩具变成趁手的生产工具。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
