企业官网建设流程全解析

简介：面向石油工程、环境工程、材料科学等领域的科研人员与工程师，这份基于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'); end

throat.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）。下表给出常见默认值和调参方向，照着做至少不会让首轮模拟偏得离谱。

标定的顺序有讲究：先用绝对渗透率实验数据锁喉道长度和形状因子，因为绝对渗透率只对这两者敏感；确认绝对渗透率对上之后，再用毛管压力曲线调收缩角。如果一上来就同时调四个参数，很容易出现"刚好凑出正确答案"的假象——换一个边界条件就打回原形。这个"先几何、后拓扑、再流动"的标定顺序，是让模型具备预测能力的底线。

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 模拟渗透率比实验高一个数量级：形状因子默认值太天真

现象：用圆管假设（G=0.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 软件包就能从玩具变成趁手的生产工具。希望帮到你。

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