简介基于Matlab开发的多孔介质孔隙网络建模软件包PNM为科研与工程人员提供完整的孔隙尺度模拟工具适用于油气储层、地下水污染迁移、生物流体传输及多孔材料设计等场景。资源共398个文件压缩包约10.11MB内含m脚本核心算法与GUI、dat数据文件孔隙/节点及链接参数、vtk可视化文件三维网络模型、html/md文档及示例图像模块划分清晰便于按需调用与二次开发。目前已有86人学习使用。软件包涵盖图像预处理与分割、孔隙网络自动生成、物理特性分析、流体传输模拟及参数化研究等完整流程可读取CT、SEM图像提取结构信息构建节点-边网络模型计算孔隙率、渗透率、比表面积并模拟单相或多相流动的路径分配与不均匀性。用户既可修改脚本实现定制分析也可在GUI环境下交互操作。附带的Berea砂岩等示例数据和详细文档能帮助用户快速掌握从图像处理到结果分析的完整工作流显著降低专业编程门槛适合多孔介质领域的科研工作者与工程师。1. 基于Matlab的多孔介质孔隙网络建模软件包把岩心CT体数据变成能算数的孔喉结构钻井取出来的岩心、电池隔膜、土壤团聚体、陶瓷滤芯只要关心“流体在这些材料里怎么流动”就离不开孔隙结构定量化。基于Matlab的多孔介质孔隙网络建模软件包做的是这样一件事拿CT扫描或其他数值重建得到的三维体数据通过最大球算法识别出一个个孔隙和连接它们的喉道输出一张可计算的网络图。这张网络把复杂的孔隙空间浓缩成“孔喉”的拓扑集合孔隙度、配位数、孔径分布、渗透率估算都可以在这张网上完成。不少初学者看到三维体数据就蒙了以为拿到图像就等于拿到结果实际上图像本身是黑匣子真正能参与计算的是网络。这个软件包的价值在于把“体素级几何”翻译成“图论级结构”让岩心分析、电池材料研究、地下水渗流模拟都有统一的中间表示。我见过太多人卡在“会看图但不会算数”这一步明明有体数据却只能靠估算孔隙度交差。适合两类人看一类手里已经有CT体数据或数字岩心想快速算出孔隙网络参数另一类做后续流动模拟需要一个可靠的上游建模工具。接下来按理论、实操、验证、避坑的顺序展开最后补一个我常用的进阶自检习惯。2. 建模原理先立住最大球算法、中轴法和孔喉图结构的取舍2.1 最大球算法为什么成为提取孔隙网络的事实标准孔隙空间的真实形态极其不规则直接把每个体素都当作节点来计算三维网格动辄上亿体素任何流动模拟都扛不住。最大球算法的思路很直观对每个孔隙体素向外扩展一个半径不断增大的球直到球面碰到固体骨架为止这个球就是该位置的最大内切球。把所有最大内切球按半径排序半径大的局部极大球就是孔隙被较大球覆盖但半径较小的球则降级为喉道的组成部分。这个算法成为事实标准有两个原因。第一它直接对应“孔隙”和“喉道”的直觉定义大球是孔连接两个大球的收缩通道是喉道。第二它对二值体数据的要求很低只要掩膜质量过关就能稳定运行不需要网格剖分和额外的几何拟合。常见做法是把它作为软件包的第一级提取器输出结果再交给后续流动模块。一个值得注意的边界条件是这里的“球”不是真实物理颗粒而是孔隙空间的局部内切球。算法输出的每个球都有半径、质心坐标和连通关系这些信息共同构成孔隙网络的几何描述。很多刚接触的人误以为球就是真实砂粒这是理解上最普遍的偏差。2.2 中轴法适合做什么骨架提取作为补充手段中轴法先把孔隙空间细化为单体素宽度的骨架线再在骨架线上标注半径距离值分支点和端点天然构成网络拓扑。它的优点是直接给出联通路径不需要复杂的合并规则缺点是对体素噪声非常敏感表面一个凸起就会在骨架上长出一个小分枝导致网络出现大量虚假喉道。我在多孔介质建模中走了几年弯路后得出的经验是中轴法更适合做二次校验而不是主提取器。比如最大球算法给出配位数异常偏高的结果时我会单独抽取中轴骨架看一眼拓扑形态确认是否存在过度分割。从这个角度上说软件包里同时保留两条路径是有意义的但默认流程应优先选择最大球。中轴法的另一个用途是计算孔隙曲折度。骨架路径长度和直线距离之比能反映流体实际行走路径的弯曲程度这个量在后续渗透率模拟里非常有用。主提取器用最大球中轴骨架只做几何诊断两者互不替代。2.3 从局部最大球到喉道图结构生成与合并规则最大球算法生成一堆半径不等的球之后关键步骤是判断哪些球属于同一个孔隙。常用做法是设定一个合并阈值当两个球的质心距离小于二者半径之和乘以某个系数且重叠体积占比较小时就把它们合并为同一个孔隙。滑移半径差异较大的球之间会留下一个狭窄的公共截面这个截面的特征半径就是喉道半径。这里的参数选择直接影响网络质量。阈值设得过高会把真正独立的两个孔隙粘成一个配位数虚高阈值设得过低又会产生大量只有喉道没有孔隙体积的碎片。实际项目中我会先把阈值放在0.5附近跑一遍统计孔隙数量、配位数分布、喉道长度分布三个指标。如果孔隙数量与岩心薄片统计值偏离超过30%基本可以判定阈值不合适。最终输出的图结构包含两个核心对象孔隙节点记录质心、半径、体积、所属区域编号喉道边记录连接的两个孔隙编号、喉道半径、长度、横截面积。这个结构与后续单相和两相模拟程序完全解耦软件包内部怎么存储不重要导出时统一成这种边-点结构就够用了。3. 用软件包跑通一次孔隙网络建模体素输入、参数给定和结果初筛3.1 第一步把三维体数据规范化成二值矩阵孔隙网络建模的输入不是原始灰度CT图像而是经过分割的语义标签体数据。常规做法是输出一个三维数组元素值为1表示孔隙体素0表示固体骨架。获得这个数组的方式可以是阈值分割、机器学习分割或数值重建对软件包而言这层接口是统一的。原始CT灰度值通常是16位整型直接作为输入会让内存占用翻倍且没有意义。我一般会在读入后立刻做二值化并转成uint8这样既能保留全部结构信息又能把内存占用控制到最小。下面这段是典型的前处理脚本% 读取重建体数据这里假设已有三维灰度数组 raw % raw 的尺寸例如 [1024 1024 1024]uint16 类型 raw load(ct_reconstruction.mat, raw); seg zeros(size(raw), uint8); % 阈值分割灰度值大于 solid 阈值即认为是孔隙 solidThreshold 16000; seg(raw solidThreshold) 1; % 清理孤立体素避免噪声点影响后续最大球提取 seg bwareaopen(seg, 1000, 6); % 去掉小于1000体素的连通域 clear raw这段代码先把数据二值化再用形态学连通域过滤清掉小孤立团块。参数1000是体积下限建议根据体素尺寸调整体素分辨率越高这个下限可以设得越小。第4个参数6表示6连通邻域三维数据下取6连通是逐体素层面的标准做法。这里有个容易踩的点bwareaopen会删掉孤立的孔隙团块如果岩心本身含有孤立孔洞这种预处理会让孔隙度偏低。取舍原则是保留尺寸大于最小孔隙半径对应体积的连通域具体数值需要通过孔径分布预估。3.2 第二步最小孔隙半径、合并阈值等关键参数怎么给提取参数放在一个结构体里传给核心函数。下面我给出常用字段的参考范围这不是某个固定仓库的接口规范而是不同实现里最通用的参数命名和取值逻辑opts.removeHalo false; % 是否裁剪边界外圈体素 opts.minPoreRadius 3; % 最小孔隙半径单位体素 opts.minThroatRadius 2; % 最小喉道半径单位体素 opts.mergeThreshold 0.5; % 孔隙合并系数通常取0.3~0.8 opts.maxIterations 500; % 最大球膨胀迭代次数过大无收益 opts.voxelSize 2.75; % 体素尺寸单位微米参数的具体含义如下。minPoreRadius是低于该半径的孔隙团块会被直接忽略符合实际物理意义数字岩心分辨率有限过小的孔隙体素不可信。mergeThreshold决定大小球合并的保守程度0.5是一个相对中庸的起点。voxelSize必须和真实物理尺寸对应它决定所有半径和长度的单位换算这个参数设错会让后续渗透率计算差一个数量级。一个最常见的参数坑是先随便给一组值跑通再看结果调参数导致每次结果都不一样。我会预先打印体数据的基本统计量包括总孔隙体素数、连通域个数、各向异性程度把参数范围固定下来后再做批量敏感性测试。3.3 第三步运行提取并解读返回的网络结构调用核心提取函数通常是单行语句返回的数据结构一般把孔隙、喉道和全局信息打包在一起。这里用示意名extractPoreNetwork表示软件包主函数实际使用时替换成对应版本的实际函数名即可。network extractPoreNetwork(seg, opts); % 核心输出结构 p network.pores; % 结构体数组radius, coord, volume, label t network.throats; % 结构体数组pore1, pore2, radius, length g network.global; % 全局信息孔隙度、体素尺寸、边界类型 % 打印基础统计量 fprintf(孔隙总数%d\n, length(p)); fprintf(喉道总数%d\n, length(t)); fprintf(平均配位数%.2f\n, 2 * length(t) / length(p));返回结果的命名习惯因软件包版本不同会有差异但pores和throats这两个字段是行业通用约定。network.global里通常还会记录提取参数的历史版本方便追溯。打印平均配位数后如果数值小于2或大于8基本可以判断参数不合理正常砂岩数字岩心大多落在2到6之间。这里再强调一个易错行为不要直接拿网络里的孔隙半径直方图和粒度分布对比孔隙半径是内切球半径不是颗粒半径两者概念不同。正确对比对象是压汞孔径分布曲线后面第4章会展开讲。3.4 第四步孔隙度、配位数、孔径分布初筛脚本拿到网络之后的第一件事是验证全局孔隙度是否和体数据分割结果一致。计算公式是全部孔隙体积加喉道体积之和除以总体积。这里有一个常见工程近似把孔隙近似为球体喉道近似为圆柱体体积可以解析计算。% 累计孔隙体积球体近似 poreVolume sum(4/3 * pi * [p.radius].^3); % 累计喉道体积圆柱近似 throatVolume sum(pi * [t.radius].^2 .* [t.length]); % 全局孔隙体积 totalVoxelVolume numel(seg) * opts.voxelSize.^3; porosity (poreVolume throatVolume) / totalVoxelVolume; fprintf(网络计算孔隙度%.3f\n, porosity); fprintf(体数据直接统计孔隙度%.3f\n, sum(seg(:)) / numel(seg));两行孔隙度的差异不应超过1个百分点。如果网络孔隙度明显偏低通常是边界处孔隙被当成了固体背景如果明显偏高大概率是分割阈值设低了把固体微孔也算进了孔隙空间。把这一步作为每次建模仿真的必做检查可以挡住绝大多数低级错误。孔径分布的初筛可以用直方图。绘制时横轴建议用对数坐标因为砂岩孔径跨度从纳米到微米线性坐标会把小孔全部挤在左侧看不见。这个直方图后续要和压汞毛管压力曲线做交叉验证所以提前就把坐标尺度统一。4. 验证网络算不算数孔隙度、绝对渗透率与入口出口边界的对照实验4.1 从网络直接算孔隙度手算验证和代码验证孔隙网络建模完成后第一件事就是独立核算孔隙度不要信任软件包自动输出的数值。手算思路是把每个孔隙球和喉道圆柱体积累加再和体数据总骨架体积做比值。这个计算和上一章的初筛脚本类似但要注意两点区别。第一喉道体积如果并入孔隙体积会虚增大约2%到5%的孔隙度取决于喉道长径比分布。标准做法是在文档里写明“体积统计包含喉道”还是“孔隙单独统计”否则后面算相对渗透率时容易产生可比性问题。第二边界孔隙的体积被截断了球体近似公式不再成立。如果边界孔隙占比超过10%应该使用逐体素统计孔隙度作为主结果网络体积仅作参考。我实际项目中一般会保留两条计算路径一条基于体素直统计一条基于网络几何累加二者在交付报告里同时保留。网络体积值的意义更多在于验证提取质量而不是替换体素统计。4.2 用Poiseuille近似估计绝对渗透率并与实测对比绝对渗透率是孔隙网络模型最常见的验证指标。把每个喉道近似为半径均匀的毛细管泊肃叶定律给出单根喉道的流量表达式再按网络拓扑组装成全局流动问题就能算出等效渗透率。% 设定入口和出口边界孔隙编号例如x方向第1层和最后一层 inlet [p.label] min([p.label]); outlet [p.label] max([p.label]); % 每个喉道的传导率忽略孔本身阻力这是常用近似 cond pi * [t.radius].^4 ./ (8 * [t.length]); % 构建稀疏线性方程组求解孔隙压力场 % 这里示意性地只显示装配过程完整求解需引入有限体积法 A sparse(length(p), length(p)); for k 1:length(t) i t(k).pore1; j t(k).pore2; A(i,i) A(i,i) cond(k); A(j,j) A(j,j) cond(k); A(i,j) A(i,j) - cond(k); A(j,i) A(j,i) - cond(k); end这段代码的物理含义是喉道看作电阻孔隙看作节点线性方程组求解结果就是各节点压力。把入口压力设为1出口设为0求解后统计入口总流量再由达西公式反推渗透率。省略号的后续步骤就是常规稀疏线性解法实际项目中我会用pcg或者直接A\b求解。这个验证只在喉道为层流且单相不可压缩时成立实际岩心中微裂缝和大孔隙同时存在时误差很明显。所以我通常把计算值和时间实测值做对数坐标对比如果偏差在一个数量级以内网络拓扑和几何基本可信偏差超过一个数量级优先检查边界条件和参数而不是怀疑算法本身。4.3 开放边界和周期性边界的稳定度测试边界处理是孔隙网络建模里最容易被忽略但也最能影响结果的部分。常见做法有三种封闭边界把所有边界孔隙都当实体开放边界保留边界孔隙并可自由进出周期性边界把对边体素连接成环形。三种方式的孔隙度、配位数计算结果差异在均匀介质里不显著但在非均质岩心里差异可能达到20%以上。实际操作中我会取同一块体数据分别用开放边界和周期性边界各跑一遍对比孔隙度和绝对渗透率。如果两者差异小于5%边界选择对结果影响不大可以放心用计算量更小的封闭边界如果差异明显说明孔隙网络在边缘被截断了需要先扩展计算域。这里有一个值得保留的习惯每次交付网络参数时把边界类型、体素尺寸、分割阈值三项一起附上。原因很简单网络本身只是一张图离开这些元数据就无法复现。软件包输出里通常会留一个字段记录这些信息应手动检查一遍不要盲信。5. 避坑孔隙网络建模常见的 5 个翻车点与排查路径5.1 边界处孔隙被当成固体孔隙度直接打五折现象体数据统计的孔隙度明明是25%提取出来的网络孔隙度只有12%左右而且配位数偏低。原因最大球算法从体数据内部膨胀遇到边界时被截断边界处大量浅层孔隙被判为无效背景。解决打开opts.removeHalofalse同时将opts.boundaryopen或将边界外圈额外标记为孔隙填充层。更稳妥的做法是把体数据在边界方向外扩若干体素填充成孔隙空间提取后再删除外扩区域对应的网络节点。这条血泪经验我反复提给同事很多所谓“孔隙网络提取效果差”的反馈其实都是边界裁切问题不是算法问题。5.2 合并阈值设错配位数暴高暴低现象平均配位数算出来接近12明显超出常见砂岩的4到6范围。原因mergeThreshold设得太高大量相邻球被错误合并成一个大孔导致每个孔连接的喉道数量虚增。解决先降到0.3重新提取观察配位数分布曲线是否从高耸单峰变成双峰或多峰。多峰出现不一定是错误如果孔隙结构本身是双模态孔径分布配位数出现两个聚集区间是正常的。判断标准是查看两组峰分别对应的孔径范围是否和薄片统计一致一致就不需要强行调参。5.3 全程用16位或双精度数组内存消耗翻倍现象1000立方体素的体数据内存占用超过6GB提取过程频繁内存交换速度极慢甚至崩溃。原因灰度CT数据、分割标签、中间球数据全部用默认double存储。解决读入后立刻转uint8保存分割结果中间最大球半径用单精度数组计算完成即时清理。一条可用的记忆方法1000立方体素的uint8数组约1GB如果内存规划超过这个规模的4倍就应该考虑分块提取。分块会引入跨块连通性问题所以内存仍充足时尽量不要主动分块。5.4 体素尺寸单位写错所有半径和渗透率一起偏离现象网络几何看起来完全正常孔隙半径也都合理但渗透率计算值和实测差出两个数量级。原因voxelSize设置从微米写成了毫米或者没有统一单位直接参与了泊肃叶公式计算。解决提取前后各打印一次voxelSize以及最小喉道半径对应的物理尺寸和压汞曲线最小可测孔径对照。这类错误最容易出现在从某个历史脚本复制参数的项目里。我的习惯是让opts.voxelSize只在主配置文件中出现一次任何子模块都从network.global.voxelSize读取杜绝参数漂移。5.5 输出字段命名不一致后处理脚本张冠李戴现象换了一个提取版本之后后处理脚本报错说找不到字段或者静默用错了半径字段最终排列图全乱。原因不同版本把喉道半径的字段名从r改成radius把coord从行向量改成列向量。解决在正式计算前加一个字段兼容适配层统一映射到内部命名再进入后续计算。这个适配层只有几行代码但值得为它单独写一个函数。命名混乱在开源生态里太常见了短视地省掉这一步后期排查会非常痛苦。6. 把孔隙网络接到单相流动模拟上一个让我少走弯路的自检习惯孔隙网络提取不是终点多数场景下最终目标是流动参数。我习惯把单相模拟分为三步走从提取结果读入孔隙喉道列表按第4章的线性系统装配渗透率再和体素级格子玻尔兹曼或有限体积结果做交叉验证。这一步能把网络本身的误差和流动求解器的误差分离开。% 读取网络并转化为标准边列表 % 标准边列表每行: [pore_i, pore_j, radius, length] edgeList zeros(length(t), 4); for k 1:length(t) edgeList(k, 1) t(k).pore1; edgeList(k, 2) t(k).pore2; edgeList(k, 3) t(k).radius; edgeList(k, 4) t(k).length; end save(network_edges.mat, edgeList, opts, porosity);我先承认一个翻车教训以前我总跳过这个验证步骤直接用网络结果推相对渗透率直到有一次两相模拟结果和实验差异巨大排查后发现是孔隙网络提取时的边界条件设错了。从那以后每次做流动模拟前我都会花10分钟检查三件事孔隙度是否一致、平均配位数是否落在合理区间、绝对渗透率是否在同一个数量级内。这个习惯也适用于软件包升级。升级提取核心后我不直接相信新版本给出的漂亮统计图而是先跑一遍相同体数据的渗透率对比确认结果一致后才会把新版本用于生产级项目。工具是拿来用的它的输出必须服务于下游的定量结论而不是停在“跑通了”这一步。如果你正准备给自己的体数据建网络模型建议从第3章的最小脚本开始把孔隙度和配位数两个指标作为质量门禁稳定之后再往流动模拟走。希望帮到你。本文还有配套的精品资源点击获取
简介基于Matlab开发的多孔介质孔隙网络建模软件包PNM为科研与工程人员提供完整的孔隙尺度模拟工具适用于油气储层、地下水污染迁移、生物流体传输及多孔材料设计等场景。资源共398个文件压缩包约10.11MB内含m脚本核心算法与GUI、dat数据文件孔隙/节点及链接参数、vtk可视化文件三维网络模型、html/md文档及示例图像模块划分清晰便于按需调用与二次开发。目前已有86人学习使用。软件包涵盖图像预处理与分割、孔隙网络自动生成、物理特性分析、流体传输模拟及参数化研究等完整流程可读取CT、SEM图像提取结构信息构建节点-边网络模型计算孔隙率、渗透率、比表面积并模拟单相或多相流动的路径分配与不均匀性。用户既可修改脚本实现定制分析也可在GUI环境下交互操作。附带的Berea砂岩等示例数据和详细文档能帮助用户快速掌握从图像处理到结果分析的完整工作流显著降低专业编程门槛适合多孔介质领域的科研工作者与工程师。1. 基于Matlab的多孔介质孔隙网络建模软件包把岩心CT体数据变成能算数的孔喉结构钻井取出来的岩心、电池隔膜、土壤团聚体、陶瓷滤芯只要关心“流体在这些材料里怎么流动”就离不开孔隙结构定量化。基于Matlab的多孔介质孔隙网络建模软件包做的是这样一件事拿CT扫描或其他数值重建得到的三维体数据通过最大球算法识别出一个个孔隙和连接它们的喉道输出一张可计算的网络图。这张网络把复杂的孔隙空间浓缩成“孔喉”的拓扑集合孔隙度、配位数、孔径分布、渗透率估算都可以在这张网上完成。不少初学者看到三维体数据就蒙了以为拿到图像就等于拿到结果实际上图像本身是黑匣子真正能参与计算的是网络。这个软件包的价值在于把“体素级几何”翻译成“图论级结构”让岩心分析、电池材料研究、地下水渗流模拟都有统一的中间表示。我见过太多人卡在“会看图但不会算数”这一步明明有体数据却只能靠估算孔隙度交差。适合两类人看一类手里已经有CT体数据或数字岩心想快速算出孔隙网络参数另一类做后续流动模拟需要一个可靠的上游建模工具。接下来按理论、实操、验证、避坑的顺序展开最后补一个我常用的进阶自检习惯。2. 建模原理先立住最大球算法、中轴法和孔喉图结构的取舍2.1 最大球算法为什么成为提取孔隙网络的事实标准孔隙空间的真实形态极其不规则直接把每个体素都当作节点来计算三维网格动辄上亿体素任何流动模拟都扛不住。最大球算法的思路很直观对每个孔隙体素向外扩展一个半径不断增大的球直到球面碰到固体骨架为止这个球就是该位置的最大内切球。把所有最大内切球按半径排序半径大的局部极大球就是孔隙被较大球覆盖但半径较小的球则降级为喉道的组成部分。这个算法成为事实标准有两个原因。第一它直接对应“孔隙”和“喉道”的直觉定义大球是孔连接两个大球的收缩通道是喉道。第二它对二值体数据的要求很低只要掩膜质量过关就能稳定运行不需要网格剖分和额外的几何拟合。常见做法是把它作为软件包的第一级提取器输出结果再交给后续流动模块。一个值得注意的边界条件是这里的“球”不是真实物理颗粒而是孔隙空间的局部内切球。算法输出的每个球都有半径、质心坐标和连通关系这些信息共同构成孔隙网络的几何描述。很多刚接触的人误以为球就是真实砂粒这是理解上最普遍的偏差。2.2 中轴法适合做什么骨架提取作为补充手段中轴法先把孔隙空间细化为单体素宽度的骨架线再在骨架线上标注半径距离值分支点和端点天然构成网络拓扑。它的优点是直接给出联通路径不需要复杂的合并规则缺点是对体素噪声非常敏感表面一个凸起就会在骨架上长出一个小分枝导致网络出现大量虚假喉道。我在多孔介质建模中走了几年弯路后得出的经验是中轴法更适合做二次校验而不是主提取器。比如最大球算法给出配位数异常偏高的结果时我会单独抽取中轴骨架看一眼拓扑形态确认是否存在过度分割。从这个角度上说软件包里同时保留两条路径是有意义的但默认流程应优先选择最大球。中轴法的另一个用途是计算孔隙曲折度。骨架路径长度和直线距离之比能反映流体实际行走路径的弯曲程度这个量在后续渗透率模拟里非常有用。主提取器用最大球中轴骨架只做几何诊断两者互不替代。2.3 从局部最大球到喉道图结构生成与合并规则最大球算法生成一堆半径不等的球之后关键步骤是判断哪些球属于同一个孔隙。常用做法是设定一个合并阈值当两个球的质心距离小于二者半径之和乘以某个系数且重叠体积占比较小时就把它们合并为同一个孔隙。滑移半径差异较大的球之间会留下一个狭窄的公共截面这个截面的特征半径就是喉道半径。这里的参数选择直接影响网络质量。阈值设得过高会把真正独立的两个孔隙粘成一个配位数虚高阈值设得过低又会产生大量只有喉道没有孔隙体积的碎片。实际项目中我会先把阈值放在0.5附近跑一遍统计孔隙数量、配位数分布、喉道长度分布三个指标。如果孔隙数量与岩心薄片统计值偏离超过30%基本可以判定阈值不合适。最终输出的图结构包含两个核心对象孔隙节点记录质心、半径、体积、所属区域编号喉道边记录连接的两个孔隙编号、喉道半径、长度、横截面积。这个结构与后续单相和两相模拟程序完全解耦软件包内部怎么存储不重要导出时统一成这种边-点结构就够用了。3. 用软件包跑通一次孔隙网络建模体素输入、参数给定和结果初筛3.1 第一步把三维体数据规范化成二值矩阵孔隙网络建模的输入不是原始灰度CT图像而是经过分割的语义标签体数据。常规做法是输出一个三维数组元素值为1表示孔隙体素0表示固体骨架。获得这个数组的方式可以是阈值分割、机器学习分割或数值重建对软件包而言这层接口是统一的。原始CT灰度值通常是16位整型直接作为输入会让内存占用翻倍且没有意义。我一般会在读入后立刻做二值化并转成uint8这样既能保留全部结构信息又能把内存占用控制到最小。下面这段是典型的前处理脚本% 读取重建体数据这里假设已有三维灰度数组 raw % raw 的尺寸例如 [1024 1024 1024]uint16 类型 raw load(ct_reconstruction.mat, raw); seg zeros(size(raw), uint8); % 阈值分割灰度值大于 solid 阈值即认为是孔隙 solidThreshold 16000; seg(raw solidThreshold) 1; % 清理孤立体素避免噪声点影响后续最大球提取 seg bwareaopen(seg, 1000, 6); % 去掉小于1000体素的连通域 clear raw这段代码先把数据二值化再用形态学连通域过滤清掉小孤立团块。参数1000是体积下限建议根据体素尺寸调整体素分辨率越高这个下限可以设得越小。第4个参数6表示6连通邻域三维数据下取6连通是逐体素层面的标准做法。这里有个容易踩的点bwareaopen会删掉孤立的孔隙团块如果岩心本身含有孤立孔洞这种预处理会让孔隙度偏低。取舍原则是保留尺寸大于最小孔隙半径对应体积的连通域具体数值需要通过孔径分布预估。3.2 第二步最小孔隙半径、合并阈值等关键参数怎么给提取参数放在一个结构体里传给核心函数。下面我给出常用字段的参考范围这不是某个固定仓库的接口规范而是不同实现里最通用的参数命名和取值逻辑opts.removeHalo false; % 是否裁剪边界外圈体素 opts.minPoreRadius 3; % 最小孔隙半径单位体素 opts.minThroatRadius 2; % 最小喉道半径单位体素 opts.mergeThreshold 0.5; % 孔隙合并系数通常取0.3~0.8 opts.maxIterations 500; % 最大球膨胀迭代次数过大无收益 opts.voxelSize 2.75; % 体素尺寸单位微米参数的具体含义如下。minPoreRadius是低于该半径的孔隙团块会被直接忽略符合实际物理意义数字岩心分辨率有限过小的孔隙体素不可信。mergeThreshold决定大小球合并的保守程度0.5是一个相对中庸的起点。voxelSize必须和真实物理尺寸对应它决定所有半径和长度的单位换算这个参数设错会让后续渗透率计算差一个数量级。一个最常见的参数坑是先随便给一组值跑通再看结果调参数导致每次结果都不一样。我会预先打印体数据的基本统计量包括总孔隙体素数、连通域个数、各向异性程度把参数范围固定下来后再做批量敏感性测试。3.3 第三步运行提取并解读返回的网络结构调用核心提取函数通常是单行语句返回的数据结构一般把孔隙、喉道和全局信息打包在一起。这里用示意名extractPoreNetwork表示软件包主函数实际使用时替换成对应版本的实际函数名即可。network extractPoreNetwork(seg, opts); % 核心输出结构 p network.pores; % 结构体数组radius, coord, volume, label t network.throats; % 结构体数组pore1, pore2, radius, length g network.global; % 全局信息孔隙度、体素尺寸、边界类型 % 打印基础统计量 fprintf(孔隙总数%d\n, length(p)); fprintf(喉道总数%d\n, length(t)); fprintf(平均配位数%.2f\n, 2 * length(t) / length(p));返回结果的命名习惯因软件包版本不同会有差异但pores和throats这两个字段是行业通用约定。network.global里通常还会记录提取参数的历史版本方便追溯。打印平均配位数后如果数值小于2或大于8基本可以判断参数不合理正常砂岩数字岩心大多落在2到6之间。这里再强调一个易错行为不要直接拿网络里的孔隙半径直方图和粒度分布对比孔隙半径是内切球半径不是颗粒半径两者概念不同。正确对比对象是压汞孔径分布曲线后面第4章会展开讲。3.4 第四步孔隙度、配位数、孔径分布初筛脚本拿到网络之后的第一件事是验证全局孔隙度是否和体数据分割结果一致。计算公式是全部孔隙体积加喉道体积之和除以总体积。这里有一个常见工程近似把孔隙近似为球体喉道近似为圆柱体体积可以解析计算。% 累计孔隙体积球体近似 poreVolume sum(4/3 * pi * [p.radius].^3); % 累计喉道体积圆柱近似 throatVolume sum(pi * [t.radius].^2 .* [t.length]); % 全局孔隙体积 totalVoxelVolume numel(seg) * opts.voxelSize.^3; porosity (poreVolume throatVolume) / totalVoxelVolume; fprintf(网络计算孔隙度%.3f\n, porosity); fprintf(体数据直接统计孔隙度%.3f\n, sum(seg(:)) / numel(seg));两行孔隙度的差异不应超过1个百分点。如果网络孔隙度明显偏低通常是边界处孔隙被当成了固体背景如果明显偏高大概率是分割阈值设低了把固体微孔也算进了孔隙空间。把这一步作为每次建模仿真的必做检查可以挡住绝大多数低级错误。孔径分布的初筛可以用直方图。绘制时横轴建议用对数坐标因为砂岩孔径跨度从纳米到微米线性坐标会把小孔全部挤在左侧看不见。这个直方图后续要和压汞毛管压力曲线做交叉验证所以提前就把坐标尺度统一。4. 验证网络算不算数孔隙度、绝对渗透率与入口出口边界的对照实验4.1 从网络直接算孔隙度手算验证和代码验证孔隙网络建模完成后第一件事就是独立核算孔隙度不要信任软件包自动输出的数值。手算思路是把每个孔隙球和喉道圆柱体积累加再和体数据总骨架体积做比值。这个计算和上一章的初筛脚本类似但要注意两点区别。第一喉道体积如果并入孔隙体积会虚增大约2%到5%的孔隙度取决于喉道长径比分布。标准做法是在文档里写明“体积统计包含喉道”还是“孔隙单独统计”否则后面算相对渗透率时容易产生可比性问题。第二边界孔隙的体积被截断了球体近似公式不再成立。如果边界孔隙占比超过10%应该使用逐体素统计孔隙度作为主结果网络体积仅作参考。我实际项目中一般会保留两条计算路径一条基于体素直统计一条基于网络几何累加二者在交付报告里同时保留。网络体积值的意义更多在于验证提取质量而不是替换体素统计。4.2 用Poiseuille近似估计绝对渗透率并与实测对比绝对渗透率是孔隙网络模型最常见的验证指标。把每个喉道近似为半径均匀的毛细管泊肃叶定律给出单根喉道的流量表达式再按网络拓扑组装成全局流动问题就能算出等效渗透率。% 设定入口和出口边界孔隙编号例如x方向第1层和最后一层 inlet [p.label] min([p.label]); outlet [p.label] max([p.label]); % 每个喉道的传导率忽略孔本身阻力这是常用近似 cond pi * [t.radius].^4 ./ (8 * [t.length]); % 构建稀疏线性方程组求解孔隙压力场 % 这里示意性地只显示装配过程完整求解需引入有限体积法 A sparse(length(p), length(p)); for k 1:length(t) i t(k).pore1; j t(k).pore2; A(i,i) A(i,i) cond(k); A(j,j) A(j,j) cond(k); A(i,j) A(i,j) - cond(k); A(j,i) A(j,i) - cond(k); end这段代码的物理含义是喉道看作电阻孔隙看作节点线性方程组求解结果就是各节点压力。把入口压力设为1出口设为0求解后统计入口总流量再由达西公式反推渗透率。省略号的后续步骤就是常规稀疏线性解法实际项目中我会用pcg或者直接A\b求解。这个验证只在喉道为层流且单相不可压缩时成立实际岩心中微裂缝和大孔隙同时存在时误差很明显。所以我通常把计算值和时间实测值做对数坐标对比如果偏差在一个数量级以内网络拓扑和几何基本可信偏差超过一个数量级优先检查边界条件和参数而不是怀疑算法本身。4.3 开放边界和周期性边界的稳定度测试边界处理是孔隙网络建模里最容易被忽略但也最能影响结果的部分。常见做法有三种封闭边界把所有边界孔隙都当实体开放边界保留边界孔隙并可自由进出周期性边界把对边体素连接成环形。三种方式的孔隙度、配位数计算结果差异在均匀介质里不显著但在非均质岩心里差异可能达到20%以上。实际操作中我会取同一块体数据分别用开放边界和周期性边界各跑一遍对比孔隙度和绝对渗透率。如果两者差异小于5%边界选择对结果影响不大可以放心用计算量更小的封闭边界如果差异明显说明孔隙网络在边缘被截断了需要先扩展计算域。这里有一个值得保留的习惯每次交付网络参数时把边界类型、体素尺寸、分割阈值三项一起附上。原因很简单网络本身只是一张图离开这些元数据就无法复现。软件包输出里通常会留一个字段记录这些信息应手动检查一遍不要盲信。5. 避坑孔隙网络建模常见的 5 个翻车点与排查路径5.1 边界处孔隙被当成固体孔隙度直接打五折现象体数据统计的孔隙度明明是25%提取出来的网络孔隙度只有12%左右而且配位数偏低。原因最大球算法从体数据内部膨胀遇到边界时被截断边界处大量浅层孔隙被判为无效背景。解决打开opts.removeHalofalse同时将opts.boundaryopen或将边界外圈额外标记为孔隙填充层。更稳妥的做法是把体数据在边界方向外扩若干体素填充成孔隙空间提取后再删除外扩区域对应的网络节点。这条血泪经验我反复提给同事很多所谓“孔隙网络提取效果差”的反馈其实都是边界裁切问题不是算法问题。5.2 合并阈值设错配位数暴高暴低现象平均配位数算出来接近12明显超出常见砂岩的4到6范围。原因mergeThreshold设得太高大量相邻球被错误合并成一个大孔导致每个孔连接的喉道数量虚增。解决先降到0.3重新提取观察配位数分布曲线是否从高耸单峰变成双峰或多峰。多峰出现不一定是错误如果孔隙结构本身是双模态孔径分布配位数出现两个聚集区间是正常的。判断标准是查看两组峰分别对应的孔径范围是否和薄片统计一致一致就不需要强行调参。5.3 全程用16位或双精度数组内存消耗翻倍现象1000立方体素的体数据内存占用超过6GB提取过程频繁内存交换速度极慢甚至崩溃。原因灰度CT数据、分割标签、中间球数据全部用默认double存储。解决读入后立刻转uint8保存分割结果中间最大球半径用单精度数组计算完成即时清理。一条可用的记忆方法1000立方体素的uint8数组约1GB如果内存规划超过这个规模的4倍就应该考虑分块提取。分块会引入跨块连通性问题所以内存仍充足时尽量不要主动分块。5.4 体素尺寸单位写错所有半径和渗透率一起偏离现象网络几何看起来完全正常孔隙半径也都合理但渗透率计算值和实测差出两个数量级。原因voxelSize设置从微米写成了毫米或者没有统一单位直接参与了泊肃叶公式计算。解决提取前后各打印一次voxelSize以及最小喉道半径对应的物理尺寸和压汞曲线最小可测孔径对照。这类错误最容易出现在从某个历史脚本复制参数的项目里。我的习惯是让opts.voxelSize只在主配置文件中出现一次任何子模块都从network.global.voxelSize读取杜绝参数漂移。5.5 输出字段命名不一致后处理脚本张冠李戴现象换了一个提取版本之后后处理脚本报错说找不到字段或者静默用错了半径字段最终排列图全乱。原因不同版本把喉道半径的字段名从r改成radius把coord从行向量改成列向量。解决在正式计算前加一个字段兼容适配层统一映射到内部命名再进入后续计算。这个适配层只有几行代码但值得为它单独写一个函数。命名混乱在开源生态里太常见了短视地省掉这一步后期排查会非常痛苦。6. 把孔隙网络接到单相流动模拟上一个让我少走弯路的自检习惯孔隙网络提取不是终点多数场景下最终目标是流动参数。我习惯把单相模拟分为三步走从提取结果读入孔隙喉道列表按第4章的线性系统装配渗透率再和体素级格子玻尔兹曼或有限体积结果做交叉验证。这一步能把网络本身的误差和流动求解器的误差分离开。% 读取网络并转化为标准边列表 % 标准边列表每行: [pore_i, pore_j, radius, length] edgeList zeros(length(t), 4); for k 1:length(t) edgeList(k, 1) t(k).pore1; edgeList(k, 2) t(k).pore2; edgeList(k, 3) t(k).radius; edgeList(k, 4) t(k).length; end save(network_edges.mat, edgeList, opts, porosity);我先承认一个翻车教训以前我总跳过这个验证步骤直接用网络结果推相对渗透率直到有一次两相模拟结果和实验差异巨大排查后发现是孔隙网络提取时的边界条件设错了。从那以后每次做流动模拟前我都会花10分钟检查三件事孔隙度是否一致、平均配位数是否落在合理区间、绝对渗透率是否在同一个数量级内。这个习惯也适用于软件包升级。升级提取核心后我不直接相信新版本给出的漂亮统计图而是先跑一遍相同体数据的渗透率对比确认结果一致后才会把新版本用于生产级项目。工具是拿来用的它的输出必须服务于下游的定量结论而不是停在“跑通了”这一步。如果你正准备给自己的体数据建网络模型建议从第3章的最小脚本开始把孔隙度和配位数两个指标作为质量门禁稳定之后再往流动模拟走。希望帮到你。本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
