马兰黄土孔隙分形特征与渗透性关系

鲁 拓，唐亚明，李喜安，任永彪，任宏玉

(1.长安大学地质工程与测绘学院, 西安 710064; 2.中国地质调查局西安地质调查中心, 西安 710054)

黄土是干旱与半干旱地区第四纪历史时期的风积沉积物，由于其特定的沉积环境和沉积时期，往往使其具有特殊的柱状节理和大孔隙架空结构,而这种疏松多孔的特质为水在土体中的运移提供了便利。随着黄土地区的工程建设日益扩大，其渗透性对工程稳定性影响愈发突出，如土坝坝基的渗透问题、土坡坡面的渗流问题、地下水运移问题等等。而孔隙空间则是作为渗透作用发生的主要场所，是影响渗透性的客观基础，其结构组成也往往与渗透性密切相关。

目前研究多孔介质孔隙结构组成的方法较多，间接方法如压汞法、低温氮吸附法、小角度X射线散射法等。直接方法如扫描电镜、CT扫描、核磁共振等。其中压汞法最为常用且已广泛应用于多孔介质的孔隙分析中。由于复杂的多孔介质孔隙结构很难依据传统几何方法准确量化其几何特征，而分形理论为解决多孔介质微观孔隙复杂程度提供了一种新的方法思路[1-2]。Pfeifer等[1]则在用分子吸附技术研究岩石孔隙结构后提出多孔介质具有分形结构的概念。随后，Rieu等[2]基于“Sierpinski”地毯模型建立了孔隙分形模型。由此，众多学者纷纷基于这种可间接测得孔隙分布含量的实验数据，运用各种分形模型计算了多种多孔介质的分形维数。Zhang等[3]基于压汞试验和热力学理论建立了孔隙表面分形维数计算模型。郁伯铭等[4]在研究大量分形多孔介质输运的特征后提出了基于孔隙率的孔隙质量分形模型及迂曲度分形模型，此模型对于评价多孔介质渗透率有着重要参考。郝兆栋等[5]基于几何模型利用压汞试验及低温氮吸附试验对腐泥煤多孔介质计算了综合分形维数。王聪乐等[6]基于三维球模型、三维毛管模型、几何模型、体积分形维数模型等利用压汞试验讨论了油气储层岩石各种分维数的关系及与渗透系数、孔隙率的关系。王艳艳等[7]综述了孔隙和颗粒分形方法对量化研究土壤结构及其水力性质具备可行性。Huang等[8]利用低温N2吸附数据对贵州龙潭组海陆过渡页岩的孔隙结构进行了分形表征，并讨论了其与渗透率的关系。Kong等[9]基于孔隙分形特征对多孔岩土材料建立了一种孔径分布的特征分形模型，以此来研究其水理性质。白耀文等[10]基于体积分形维数模型利用压汞试验对延长组致密砂岩油气储层岩样进行了分维计算，并探讨了与渗透率关系。Dou等[11]基于压汞试验及图像分析等技术讨论了孔隙不同分形模型对评价砂岩渗透率的影响，发现不同模型计算的分形维数不同，且分形维数越大，渗透系数越低。不仅如此，郁伯铭等[4]还认为孔隙分形的自相似区间最小值与最大值应至少相差两个数量级才有综合反应结构的意义。

综上，利用间接测量多孔介质孔隙分布含量，或者说得到其孔隙分布函数，并以此来计算孔隙分形维数，并用分形维数表征孔隙结构的复杂性，屡见不鲜。但是，基于各种分形模型利用压汞试验主要对马兰黄土孔隙分形维数的探讨及其与渗透性的关系方面的研究鲜有报道，以及哪一种模型在什么样的尺度范围内更加适用于表征马兰黄土孔隙结构特征都值得讨论。用渗透系数表征渗透性大小，基于各种分形模型利用压汞试验探讨马兰黄土分形维数特征与渗透性关系，这对马兰黄土孔隙结构量化研究及其对渗透性的影响具有理论参考价值，对实际黄土工程地区渗流破坏及其防治有实际参考意义。

1 试验材料和方法

1.1 试样基本物理指标

试样选自吕梁地区柳林县、隰县、大宁县、吉县、乡宁县、汾西县、孝义市、汾阳市、交城县、岚县10个地区的典型原状Q3马兰黄土(图1)，采取深挖取样方式获得原状黄土试样，表1为试样基本物理指标及采用室内变水头渗透试验测试的渗透系数。

图1 Q3马兰黄土取样地点及现场

表1 基本物理指标及渗透系数

1.2 孔隙分布特征与渗透系数关系

孔隙分析采用美国麦克仪器AutoPore IV 9500高性能全自动压汞仪测试，测试压力范围0.03～419 MPa，汞表面张力485 mN/m，汞接触角 130°。颗粒平均粒径采用Bettersize 2000激光粒度仪测试然后再求取平均粒径，为后文用Kozeny-Carman方程计算理论渗透率提供试验参数。

通过压汞法所测孔隙结构参数排驱压力Pcd的定义为进汞曲线中平缓处的拐点[图2(a)，试样11]，若进汞曲线拐点不明显则将第i次与i+1次进汞的饱和度SHg差值首次大于3%的点作为排驱压力点[图2(a)，试样9]。雷祥义[12]首次应用压汞法定量揭示了黄土中不同大小孔隙的分布，并按照所测孔隙直径d的大小将孔隙划分为微孔隙n1(2 μm≥d)、小孔隙n2(8 μm≥d>2 μm)、中孔隙n3(32 μm≥d>8 μm)、大孔隙n4(32 μm

P为进汞压力，MPa；SHg为进汞饱和度，%；n1、n2、n3、n4分别为微孔隙、小孔隙、中孔隙、大孔隙的百分含量

2 分形特征

2.1 分形模型及计算结果

分形理论的一个最基本关系式就是分形物体的度量M(ε)与测量尺度ε服从标度关系[1,5]

M(ε)∝εD

(1)

式(1)中：D为分形维数；M(ε)为一个物体的质量、体积、面积、长度等参数；ε为测量尺度。

当测量尺度ε在自相似区间范围内改变时，分形物体度量M(ε)也随之改变，但分形维数D基本不变，对式(1)两边取对数则得到一条直线，其斜率为分形维数。虽然分形理论已广泛应用于多孔介质的孔隙量化分析中，但是不同的多孔介质在什么样的自相似区间内应用何种模型仍值得探究。由“Sierpinski”地毯模型和基于实测孔隙分布含量衍生的计算孔隙分形维数模型有:三维球模型[9]、三维毛管模型[9]、几何模型[9]、体积分形模型(润湿相模型)[4,9]、质量分形(面积分形)模型[4]，由热力学理论和实测孔隙分布函数衍生建立的分形模型有孔隙表面分形模型[3]，值得一提的是质量分形模型强调孔隙率，体积分形模型强调孔径搭配。本文中整理了这些计算模型的表达函数，如表2所示。并且，尺度范围、模型类型、参数选取更是分形理论正确应用的前提。尺度范围问题，即孔径上下限选取问题(自相似区间问题)以及选取该区间是否具有综合代表性问题。

表2 分形模型

但是，由于实际测量试验的误差，上下限还受控于测试手段的缺陷，如当压力超过10 MPa时会使得疏松多孔介质发生压缩甚至破裂从而改变原生孔隙系统[13]，故经过筛选选取孔径大于0.2 μm以上的孔隙参与分形计算[图3(a)]。依据分形理论，分形自相似区间应该至少相差两个数量级dmin/dmax≤10-2才具有分形计算意义[14]，所算维数才具有一定的代表性，区间数量级往往记为dmin/dmax≤10-2和dmin/dmax≤10-3，d为所测直径。结合本文压汞试验实际，试样最大孔径为429.59 μm，故选取的尺度范围为三种：范围Ⅰ：430～2 μm；范围Ⅱ：32～0.2 μm；范围Ⅲ：430～0.2 μm。Ⅰ、Ⅱ区间数量级为dmin/dmax≤10-2，Ⅰ区间涵盖大、中、小孔隙，Ⅱ区间涵盖中、小、部分微孔隙。Ⅲ区间数量级为dmin/dmax≤10-3，涵盖大、中、小、部分微孔隙。模型选取问题，考虑该尺度范围内采取该模型计算时是否具有线性回归的特征。我们采用线性回归的拟合度R2(>0.75)来判定是否具有线性特征，再者通过模型计算的分形维数值理应介于2～3，最后通过该模型计算的分形维数与渗透系数是否具有相关性来决定是否选取该模型。参数选取问题，主要为各模型的因变量与自变量特征及进一步的选取孔径数值d的问题。第i次进汞和第i+1次进汞区间内，代表此孔隙体积区间的孔径一般为调和平均数dh、几何平均数dg、算术平均数da、平方平均数ds这四种，四种孔径平均数对计算分形维数是否有影响及影响多大值得探讨。

因此，对于Q3马兰黄土的分形模型也主要基于渗透系数并从尺度范围、模型类型、参数选取三个方面进行讨论(表2、图3)。选取的参与分形维数计算的孔径范围下限(0.2 μm)为大于10 MPa所对应的孔径，此范围的进汞过程中，土样未发生压缩，故具有合理性[图3(a)]。由图3(b)知通过四种孔径平均数计算的分形维数相差甚小，并应用其他模型计算得到类似结果，故涉及孔径参与的分形维数的计算本文采用算术平均数da代表d。由图3(c)知，相同试样的相同尺度范围内不同模型计算的分形维数相差很大，其中三维球模型计算的分形维数最大(>3)。相同试样不同尺度范围内的分形维数相差很大，其中范围Ⅰ与范围Ⅱ相差较近，范围Ⅰ计算的分形维数普遍大于范围Ⅱ、Ⅲ[图3(b)、图3(c)]，这说明孔径越大，孔隙的空间展布越是复杂，分形维数越大。通过图3(c)，筛选出合理的分形维数，分别是：范围Ⅰ质量分形模型计算的Ⅰ-D5、孔隙表面分形模型计算的Ⅰ-D6，范围Ⅱ内质量分形模型计算的Ⅱ-D5、孔隙表面分形模型计算的Ⅱ-D6；范围Ⅲ三维毛管模型计算的Ⅲ-D2、体积分形模型计算的Ⅲ-D4、质量分形模型计算的Ⅲ-D5、孔隙表面分形模型计算的Ⅲ-D6，部分计算过程如图3(d)所示。

图3 Q3马兰黄土分形维数计算

2.2 分形与渗透系数关系

理论渗透率，亦称固有渗透率，采用半经验Kozeny-Carman方程计算，记为k′，表达式为

(2)

式(2)中：平均粒径采用几何平均粒径来等效[15]。假定室内试验的渗流过程符合达西定律，得

(3)

(4)

式中:φ为孔隙率；CKC为K-C常数，通常通过试验测定；ρw为水的密度，kg/m3；g为重力加速度，m/s2；μw为水的动力黏滞系数，Pa·s；C为常数，经过与实测数据进行校正，取C=1.592×10-6，由式(4)可得理论渗透系数。

将拟合程度大于0.7以上的分形维数认为是与渗透系数K具有相关性的，由图4知，基于孔隙率的质量分形模型计算的分形维数在三个尺度范围内均与渗透系数呈正相关，相关性均达到了0.79以上，这说明，随着质量分形维数增大渗透性在增强。并且无论是大、中、小孔隙组合，或中、小、部分微孔隙组合、还是大、中、小、部分微孔隙组合，其孔隙率始终是渗透性的主要影响因素。基于孔径搭配的体积分形模型计算的分形维数只在尺度范围Ⅲ:430～0.2 μm内与渗透系数呈负相关，相关性为0.70以上，这意味着随着体积分形维数增大渗透性在下降。并且仅当大、中、小、部分微孔隙组合时才能综合反映孔隙的搭配复杂程度，且这种搭配结构复杂性对渗透性有重要影响。

图4 Q3马兰黄土孔隙分形维数与渗透系数相关性分析

孔隙结构迂曲程度对渗透作用的发生同样影响深远，且迂曲度已作为渗透性评价的重要指标[4]，而迂曲度分形维数则是刻画这种迂曲程度的度量。郁伯铭等[4]基于孔隙分形维数对迂曲度分形维数进行了一般表达，即

(5)

(6)

式中：DT为迂曲度分形维数；τav为孔隙平均迂曲度；D为孔隙三维分形维数减去1，或者为孔隙二维分形维数。

图5为通过质量分形维数D5、体积分形维数D4计算的迂曲度分形维数，即将其减去1认为是式(6)中的D。结果表明，通过此两者分形维数计算的迂曲度分形维数均与渗透系数呈良好的负相关，相关性均大于0.79。这说明随着迂曲度分形维数的增加，孔隙迂曲程度也在增加，Q3马兰黄土的渗透能力则在减小，同时也说明所筛选的分形维数具有合理性。

图5 Q3马兰黄土孔隙迂曲度分形维数与渗透系数相关性分析

图6为质量分形维数和体积分形维数与其他参数的相关性分析。排驱压力Pcd、微孔隙含量n1与质量分形维数呈负相关，相关性均大于0.70，从孔隙演化的角度看，随着沉积环境、颗粒组分、密度等的改变，大的孔隙在被压缩转化成较小孔隙的过程中，总孔隙率在减小的同时而较小孔隙含量却在增加。故微孔隙含量会影响孔隙率再而影响质量分形维数。中值孔径与质量分形维数呈指数正相关，相关性0.74以上，甚至高达0.90以上，孔径越大孔隙率就越大，导致质量分形维数越大。而前述变化规律与体积分形维数则恰好相反，微孔隙含量增加导致排驱压力增大，中值孔径也会增大，进而增加孔隙搭配结构的复杂性，导致体积分形维数增大。

We为退汞效率，%;Sp为分选系数；Dr为变异系数；dg为中值粒径，μm

综上，质量分形维数强调的是尺度范围内的孔隙率，孔隙率越大，其值也越大。而体积分形维数强调的是孔径搭配，其值越大表示孔径搭配结构越是复杂且较小孔隙含量越多。两者均作为影响Q3马兰黄土渗透性的重要因素。

3 结论

基于压汞试验和分形理论,在不同尺度范围应用多种模型计算分形维数并与实测及理论渗透系数进行相关性分析，取得如下认识。

(1)基于渗透性的Q3马兰黄土孔隙适用的分形模型为质量分形模型和体积分形模型，质量分形维数在尺度范围430～2 μm、32～0.2 μm、430～0.2 μm内均与渗透系数呈正相关，相关性0.79以上，说明孔隙率始终是左右渗透性的关键因素。在尺度范围430～0.2 μm体积分形维数与渗透系数呈负相关，相关性0.70以上，说明大、中、小、部分微孔隙组合时才能综合反映孔隙的搭配复杂程度。同时通过它们计算的迂曲度分形维数均有渗透系数呈负相关，相关性0.79以上，说明了这两种模型计算的分形维数具有合理性。

(2)排驱压力、微孔隙含量与质量分形维数呈负相关、中值孔径与其呈正相关，相关性0.74以上，而这种变化规律与体积分形维数恰好相反。大的孔隙在转化成较小孔隙的过程中，总孔隙率在减小的同时而较小孔隙含量却在增加，增加孔隙搭配结构的复杂性，进而影响渗透性，造成质量分形维数减小的同时体积分形维数增大。因此基于孔隙率的质量分形维数与基于孔径搭配的体积分形维数对马兰黄土孔隙结构量化有理论参考价值。对黄土工程地区渗流破坏及其预防有实际参考意义。