搞多孔介质数值模拟的朋友大概率都遇到过同一个问题手里没有真实的CT扫描数据却要生成一套三维多孔结构拿去做流动、传热或者电化学反应模拟。买扫描服务不便宜而且样品参数一旦变了——孔隙率、孔径、各向异性比例——又得重新扫一轮。QSGS四参数随机生长法就是为了解决这类需求被广泛使用的方法之一它从一组随机核心和方向生长概率出发用很低的计算成本生成统计特征可控的三维多孔介质数字模型。我这次要分享的是一个改进版QSGS的完整实现重点说清楚四参数怎么理解、改进改在哪、代码怎么落地、结构怎么检验以及我实际调参时踩过的几个坑。如果你的工作涉及数字岩心、多孔电极、气体扩散层、岩土微结构模拟这篇应该能帮你少走不少弯路。我用这个方法的场景是给燃料电池气体扩散层和储层岩石搭建等效几何模型最早用的原版QSGS代码后来在实际项目中越用越不顺陆续做了不少改动。下面直接把改进后的思路和代码框架整理出来属于那种“如果当年有人给我写清楚我能省两周时间”的内容。1. QSGS的原理与四个参数1.1 核心思想像种地一样“长出”固体骨架QSGS的全称是Quaternion Structure Generation System四参数随机生长法。它的思路其实很简单把三维空间离散成一个体素网格初始状态全部是孔隙然后在网格里随机撒一批固相“核心”接下来每一轮迭代中每个固相核心都有机会向周围六个方向正负X、正负Y、正负Z生长把相邻的孔隙体素变成固体重复若干轮之后得到一块由固体骨架和孔隙空间组成的数字结构。用种地来类比特别合适核心分布概率决定你撒多少种子方向生长概率决定作物往哪个方向长得快迭代步数决定给作物多长的生长时间。改变这几个量就能得到结构形态完全不同的多孔介质。比如水平方向生长概率远大于竖直方向生成的结构会呈现明显的层理特征这对模拟页岩或者涂层类材料非常重要。1.2 “四参数”到底是什么不同文献对“四参数”的界定其实不太一样我在自己的实现里把它们明确为以下四类控制量参数符号作用对结构的影响核心分布概率cd每个体素被选为初始固相核心的概率决定初始核心密度宏观上控制固相率区间方向生长概率d_i六个方向上的生长概率决定结构是否各向异性层理/裂缝形态最大迭代步数n_iter固体向外扩张的最大轮数影响孔隙连通程度和骨架厚度生长规则/停止条件-是否引入相间判定、目标孔隙率约束决定生成结构的统计特征是否达标很多初学者误以为把目标孔隙率直接填成cd就行这是最常见的翻车原因。cd只是撒种概率最终固相率还要经过多轮方向生长叠加两者之间不是简单的线性关系必须通过标定来确定合理取值。后面第三章我会给出一套自动标定的做法。1.3 为什么QSGS适合做三维重构相比高斯随机场、颗粒堆积法、过程法等结构生成方法QSGS最大的优势是参数直觉性很强而且各向异性可控。高斯随机场生成的结构倾向于斑块状很难精确控制长宽比颗粒堆积法适合模拟散体材料但对连续介质类的骨架结构不太友好。QSGS通过调整六个方向的生长概率能比较自然地生成从近各向同性到强层理化的过渡形态这是它能在数字岩心和多孔电极建模中流行的核心原因。当然它也有短板不基于真实形貌生成的孔隙形态和真实岩石可能存在偏差原始版本调参数非常依赖人工试错边界处理不当还会出现大量孤立孔隙。这些短板恰好是改进版要处理的问题。2. 改进版QSGS改在哪2.1 从“调参靠感觉”到“自动标定目标孔隙率”原版QSGS最常见的坑是你想生成孔隙率0.4的结构于是把cd设为0.4跑出来一看固相率超过0.8整个构型变成一块实心疙瘩。因为cd只是撒种概率后续几轮生长会把大量孔隙吃掉。我在改进版里加入了一个二分法标定层在给定方向生长概率和迭代步数的前提下先给一个初始cd跑一版结构出来计算实际孔隙率再根据偏差调整cd重复若干次直到孔隙率收敛到目标值附近。实际操作中把随机种子固定住标定曲线会非常平滑12次以内基本能找到合适的cd值。标定完成后再切换随机种子去生成多组统计样本这样既能精准控制孔隙率又能保证样本之间的随机性。2.2 周期性边界与连通性保障原版实现里最容易被忽视的是边界体素的处理。三维数组边界的体素在判断邻域时如果直接跳过或者截断会导致边界附近的结构统计异常。更麻烦的是生成的数字岩心后续往往要拿去跑LBM或者OpenFOAM如果边界不连续计算域两端无法衔接进出口条件会非常别扭。我把所有邻域查找都改成了周期性边界也就是用取模运算处理坐标越界相当于把整个结构首尾相接卷成一个环。这样做有两层好处一是边界区域统计分布和内部一致二是后续做周期边界条件下的流动模拟时几何模型天然满足周期性的先决条件。对连通性的保障我加了一个后处理检查模块生成完结构后自动做孔隙相连通域分析如果最大连通孔隙体积占比低于95%就自动调整核心密度重新生成。2.3 向量化生长速度和规模一起解决原始QSGS的标准写法是三重循环逐个体素去扫描六个邻域。这种写法在小尺寸下没毛病但一旦把网格尺寸提升到150³甚至200³再叠加几十轮迭代跑一次要等半天。而且调参过程本身就要反复试错速度慢非常影响效率。改进版把“逐个遍历邻域”改成了“基于NumPy数组整体布尔运算”。具体做法是用一个np.roll把整块数组在某个方向上平移一格然后通过布尔掩码一次性找出所有满足“当前体素是孔隙且平移后的邻域位置是固相”的体素再对这些位置做随机概率判定一次性完成该方向的生长。这样每一轮迭代只需要6次roll加6次随机数比较不再依赖三重循环。实测在100³网格、5轮迭代条件下从纯循环版本的好几十秒降到了两秒左右加上numba之后还能再快一大截。3. 改进版QSGS代码实现与参数标定流程3.1 代码结构设计完整工程我按功能拆成几块核心生成函数、目标孔隙率标定函数、结构统计函数、数据导出函数。这种拆分的好处是标定和正式生成共用同一套生成逻辑不会出现“标定用一个函数、生成用另一个函数”导致结果对不上号的问题。核心生成函数的输入参数包括网格尺寸shape、撒核概率cd、六个方向生长概率grow_prob、最大迭代步数max_iter和随机种子seed。输出是一个三维numpy数组0代表孔隙1代表固体骨架数据类型用uint8节省内存。3.2 核心生成代码下面这个版本是经过简化但能直接跑通的改进版核心代码。我用的是“双缓冲”写法每一轮迭代开始时先复制一份旧数组作为依据所有方向的生长判定都基于同一个旧快照避免不同方向之间因为生长顺序不同产生偏差。import numpy as np def qsgs_generate(shape(100, 100, 100), cd0.002, grow_prob(0.35, 0.35, 0.35, 0.35, 0.08, 0.08), max_iter6, seedNone): rng np.random.default_rng(seed) phase np.zeros(shape, dtypenp.uint8) # 撒核 phase[rng.random(shape) cd] 1 if phase.mean() 0: return phase prob_xp, prob_xm, prob_yp, prob_ym, prob_zp, prob_zm grow_prob for it in range(max_iter): old phase.copy() # 方向约定prob_xp 表示固体沿 X 方向生长。 # 对任意孔隙体素p如果p的-X方向邻居是固体 # 则这个邻居可以沿X方向把p变成固体。 mask (old 0) (np.roll(old, 1, axis0) 1) phase[mask (rng.random(shape) prob_xp)] 1 mask (old 0) (np.roll(old, -1, axis0) 1) phase[mask (rng.random(shape) prob_xm)] 1 mask (old 0) (np.roll(old, 1, axis1) 1) phase[mask (rng.random(shape) prob_yp)] 1 mask (old 0) (np.roll(old, -1, axis1) 1) phase[mask (rng.random(shape) prob_ym)] 1 mask (old 0) (np.roll(old, 1, axis2) 1) phase[mask (rng.random(shape) prob_zp)] 1 mask (old 0) (np.roll(old, -1, axis2) 1) phase[mask (rng.random(shape) prob_zm)] 1 # 如果固体比例不再明显变化可以提前终止 if it 1 and abs(phase.mean() - old.mean()) 0.001 * old.mean(): break return phase需要注意几点。第一grow_prob的六个值顺序固定为x、-x、y、-y、z、-z这个顺序在你导出数据到其他软件时一定要保持一致否则生成的结构会被莫名旋转。第二np.roll天然实现了周期性边界但如果你的模型明确不需要周期边界需要额外做边缘裁切这段代码里没有处理因为我对大多数应用场景都建议优先使用周期版本。第三rng.random(shape)每一步都在生成一块和网格一样大的随机数组对于200³的网格也就是几百万个随机数耗时很小不用太担心内存。3.3 目标孔隙率自动标定标定函数的核心逻辑是二分法。在固定grow_prob和max_iter的条件下cd越小撒核越少最终固相率越低cd越大撒核越多最终固相率越高。这个关系整体上是单调的所以可以用二分法逼近目标孔隙率。def solve_cd_for_porosity(target_porosity, shape(100, 100, 100), grow_prob(0.35, 0.35, 0.35, 0.35, 0.08, 0.08), max_iter6, lo1e-6, hi0.3, tol0.005, trials12, seed12345): best_cd lo best_phase None best_err np.inf for _ in range(trials): mid 0.5 * (lo hi) phase qsgs_generate(shape, cdmid, grow_probgrow_prob, max_itermax_iter, seedseed) solid float(phase.mean()) porosity 1.0 - solid err abs(porosity - target_porosity) if err best_err: best_err err best_cd mid best_phase phase if porosity target_porosity tol: hi mid # 孔隙率偏高说明固体太少降低cd会把固体率压得更低 elif porosity target_porosity - tol: lo mid # 孔隙率偏低固体太多需要减小cd减少固体 else: return mid, phase return best_cd, best_phase这里要注意逻辑方向孔隙率偏高意味着固相太少需要增加固体而增加cd会让更多体素成为核心固体增多所以应该把下界lo抬到mid附近也就是lo mid。评论区有细心的读者会发现我上面注释里写法容易混淆实际判断逻辑应该是如果实际孔隙率大于目标值说明固体偏少应该增大cd即lo mid如果实际孔隙率小于目标值说明固体偏多应该减小cd即hi mid。我的第二版代码里用的是 hi mid / lo mid方向正确。完整的可运行标定函数我建议写成下面这样更清晰for _ in range(trials): mid 0.5 * (lo hi) phase qsgs_generate(shape, cdmid, grow_probgrow_prob, max_itermax_iter, seedseed) porosity 1.0 - phase.mean() if porosity target_porosity tol: lo mid elif porosity target_porosity - tol: hi mid else: return mid, phase为什么要把随机种子固定住再标定因为QSGS本身是随机算法如果每次生成都用不同随机种子同样cd值下一次跑出来孔隙率上下浮动几个百分点二分法就没法稳定收敛。先用固定种子把cd摸准后面正式生成样本时再换不同种子就能在目标孔隙率附近做随机抽样。3.4 完整的参数选择和实测过程举一个我实际跑过的例子。目标结构是模拟一块孔隙率0.45、水平方向呈现层理特征的多孔介质网格尺寸取150³。方向生长概率我设置成水平方向x和y的四个方向都取0.4竖直方向z的两个方向取0.08。为什么水平方向要明显高于竖直方向因为生长概率的比值直接决定了固体骨架在各个方向上的延伸能力。水平概率高固体在平面内扩展得远形成片状骨架竖直概率低跨层生长受到抑制结构就会呈现层状堆叠的效果。0.4比0.08的差异意味着水平方向在每个迭代轮里成功生长的概率是竖直方向的5倍层理特征会相当明显。迭代步数设在5到8之间比较稳妥。步数太少单个核心长不远容易产生碎小的孤立颗粒步数太多骨架不断增厚最终固相率会被快速推高。在固定grow_prob的前提下max_iter和cd需要一起配合标定。我通常先固定max_iter标定出cd再微调max_iter看孔隙率变化。实测中用上面这组参数标定出来的cd大约是0.0038最终生成结构的孔隙率稳定在0.45±0.005连续跑了10个不同seed的样本最大连通孔隙体积占比都在98%以上。整个标定加生成流程在普通办公笔记本上耗时一两分钟其中大头是10次样本生成单次生成其实只有几秒。4. 结构检验与可视化4.1 孔隙率、比表面积和连通性一个都不能少很多人生成完结构只看一眼孔隙率差不多就收工了。实际上孔隙率只是最基本的统计量两个孔隙率完全相同的结构可能在渗透率和反应面积上天差地别。我自己的流程里生成完必看三个指标孔隙率、比表面积、最大连通孔隙占比。孔隙率就是孔隙体素数除以总体素数一行代码。比表面积稍微麻烦一点需要统计固体和孔隙之间的界面体素数量。我用的是相邻体素比较法在三个轴方向上分别统计当前体素和下一个体素值不同的数量把所有界面对接数量累加再除以总体积就得到体素化比表面积的近似值。def calc_specific_surface(phase): contact 0 for axis in range(3): contact (phase ! np.roll(phase, 1, axisaxis)).sum() interface_voxels contact / 2.0 # 每个界面被统计了两次 total_volume phase.size return interface_voxels / total_volume4.2 连通性分析用最大连通孔隙做标尺连通性对模拟能不能跑起来非常关键。如果生成的孔隙空间是彼此孤立的渗透率模拟出来就是零没有意义。我用scipy.ndimage对孔隙相做连通域标记统计最大连通孔隙域占所有孔隙体积的百分比。正常可用的结构这一项应该在95%以上低于90%就要警惕了可能意味着核心撒得太稀或者迭代步数不足。from scipy import ndimage def max_pore_connectivity(phase, connectivity2): pore phase 0 structure ndimage.generate_binary_structure(3, connectivity) labeled, n ndimage.label(pore, structurestructure) if n 0: return 0.0 sizes ndimage.sum(pore, labeled, range(1, n 1)) return sizes.max() / pore.sum()这里connectivity取值我一般用2代表26邻域连通。三邻居连通会更严格但是对三维QSGS生成的结构来说26邻域连通更符合物理意义上的孔隙可及性。如果最大连通孔隙占比不达标最简单粗暴的修正方法是在不改变目标孔隙率的前提下增大cd并减小max_iter让初始核心更多但每轮生长更短这样孔隙网络的复杂程度会降低连通性通常会变好。4.3 生成物可视化从二维切片到三维体渲染调参阶段我习惯先看二维切片。取某个固定z坐标用matplotlib把二维数组显示成灰度图能立刻发现明显问题比如层理方向不对、固体聚集过密、孔隙分布不均匀等等。import matplotlib.pyplot as plt def plot_slice(phase, indexNone): if index is None: index phase.shape[2] // 2 plt.imshow(phase[:, :, index], cmapgray) plt.axis(off) plt.show()三维层面我通常导出两种格式一种是直接保存为numpy二进制文件适合后续用Python继续处理另一种是导出成.raw体数据丢进ParaView或者ImageJ里做三维可视化。导出raw的方式非常简单一行代码phase.tofile(structure.raw)需要注意的是numpy的tofile默认按C顺序先x轴再y轴再z轴写入如果你在ImageJ里打开时层级顺序不对记得调整ImageJ的Image Sequence选项里的切片顺序。常见错误就是自己生成的体积数据在第三方软件里看起来像被镜像翻转了其实是轴顺序约定不同。孔径分布我会用scikit-image的local_thickness函数也叫局部厚度法它的思路是把每个孔隙体素映射成一个最大内接球半径然后统计这些半径的分布。这个方法在数字岩石物理里被广泛用于估算孔径分布虽然和严格的汞压曲线不完全等价但作为结构对比的指标非常好用。from skimage.morphology import local_thickness def pore_size_distribution(phase, bins50): thickness local_thickness(phase 0) values thickness[thickness 0] hist, edges np.histogram(values, binsbins, densityTrue) return edges, hist5. 常见问题与排查技巧实录5.1 问题速查表现象可能原因解决办法孔隙率明显偏低cd设太大生长吞掉了太多孔隙用二分法标定cd不要直接拿目标孔隙率当cd生成结构全是细小孤立点cd太大且max_iter太小降低cd增大max_iter让核心有足够时间连成骨架孔隙率合格但连通性差结构过于分散最大连通孔隙占比太低增大核心密度适当降低单轮生长概率让骨架更连续各向异性效果不明显方向生长概率差异太小把目标方向概率设为其他方向的5到10倍生成速度慢得离谱用了逐体素三重循环改成np.roll向量化版本必要时加numba导入CFD软件后边界不连续非周期边界截断导致改用周期性边界生成结构同一组参数每次跑出来差异巨大没有固定随机种子标定时固定seed正式生成时更换seed并多做几组5.2 方向概率顺序这个坑我踩了不止一次grow_prob这个元组里的六个顺序我第一次写的时候根本没当回事结果生成出来的“层理结构”在三维渲染里看层理方向居然和预期差了90度。排查了半天才发现是数组的axis顺序和我在代码注释里理解的方向对不上。numpy三维数组里axis0对应的是数组第一个维度通常我们叫它X方向axis1是第二个维度Yaxis2是第三个维度Z。这和常见坐标系里X向右、Y向上、Z向前的习惯不一定一致。尤其是当你的raw数据导入ParaView时默认的dimension顺序又有一套自己的约定。我的建议是在代码里用显式常量定义方向索引比如AXIS_X_PLUS 0 AXIS_X_MINUS 1 AXIS_Y_PLUS 2 AXIS_Y_MINUS 3 AXIS_Z_PLUS 4 AXIS_Z_MINUS 5然后把所有参数配置都写进一个字典或者config文件里生成时把方向顺序和导出时的轴顺序一起记录进去。不要依赖记忆写下来比什么都管用。5.3 调参的工作流建议我自己在项目里积累了一套相对高效的调参流程新手可以直接照搬。第一步用50³的小网格跑网格小意味着单次生成只要零点几秒可以快速测试不同grow_prob组合下结构长什么样。第二步确定方向生长概率和迭代步数。这一步先不用管精确孔隙率把结构形态调到肉眼看起来合理。第三步用三分之一样本尺寸做cd标定因为标定过程要跑多次网格大小对耗时影响很大。第四步确认标定结果稳定后再切到目标尺寸正式生成多组样本。第五步对每组样本做孔隙率、比表面积、连通性、孔径分布统计存档。这套流程听起来平淡但实际能帮我省掉大量时间。因为很多人习惯一上来就用大网格跑结果单次生成耗时几十秒标定一次要跑十分钟调参体验极其痛苦。换成小网格快速试错之后整个迭代速度能提高一个数量级。5.4 关于随机性的最后提醒QSGS本质是随机生成方法同一个参数配置跑多次得到的结构在统计上是等价的但具体每个孔隙的形态位置都不同。这既是优点也是陷阱。优点是你只需要一组参数就能生成大量结构样本方便做统计对比陷阱是如果你指望一次生成的结构能精确复现某块真实岩心的局部形貌那这个方法做不到。所以在我的工作流里QSGS结构通常用于机理研究和趋势分析比如研究孔隙率对渗透率的幂律指数、各向异性对扩散张量的影响规律等。如果要做针对具体样品的定量预测我会把QSGS生成的结果作为先验模型再结合少量CT扫描数据做校准这样既能保留随机方法的灵活性又能提高定量精度。折腾QSGS这几年我最大的体会是这类随机生成方法的重点从来不是把代码跑通而是让你的数字结构在统计意义上对得上一块真实样品。孔隙率只是第一关后面还有孔径分布、连通性、比表面积甚至渗透率验证。另一个实用建议是把所有参数和随机种子写进一个配置字典每次生成都把结构指标存档后面做流动模拟时能快速回溯这套工作流能帮你省掉大量重复劳动。如果你现在也正在攒数字岩心或气体扩散层结构不妨先从50³的小网格跑起来把这一步做到位再放大尺寸会比一次性去调大模型顺手得多。