基于CT的黄土大孔隙形态三维分形研究

王 凤，邓念东，马逢清，王 超，江星辰

(西安科技大学地质与环境学院，西安 710054)

1 概述

黄土自身物理力学性质加之外部诱发因素而导致的地质灾害已经严重制约了黄土地区经济的发展[1]，而我国黄土分布广泛，有效解决这些黄土地质问题就显得至关重要。黄土的物理力学性质很大程度上取决于土体结构，因此，土力学纵深发展的关键在于结构性问题的解决。要探索黄土体结构与黄土灾变及宏观力学响应的关系，先要对土体结构进行定量量测与分析，从而进一步研究黄土体结构是如何影响土体内部应力和变形分布，以及如何影响与改变土体的破坏过程和规律[2-3]。

基于新型测量技术与图像处理技术的发展，国内外许多学者在土体微结构领域开展了很多研究工作[4-5]。这些研究多借助扫描电镜(Scanning Electron Microscope,SEM)技术获取一系列二维黄土断面扫描图像(即剖面)，然后利用分形理论对黄土内部结构进行精细量测和数字表述，进一步求解结构参数。

分形理论作为研究不规则几何对象的计算机图形技术的理论正不断深入到岩土工程的各个研究领域[6]。例如：运用分形理论对高铁地基黄土的分形维数与其孔隙比之间的关系进行了研究[7]，利用分形理论获得了土壤颗粒分布数量的多重分形参数，更细致地刻画了土壤颗粒的数量分布特征[8]，用分形理论揭示了岩体数据的内在关系，深层次地反映节理岩体的工程地质底蕴[9]。

这些方法一定程度上揭示了黄土微结构定量参数与其分形维数的关系，但仍然局限于二维，无法做到真正地以三维形式揭示黄土微结构定量分析的方法。以往的理论认为三维表面的分形维数DS=D+1(D为二维轮廓的分形维数)，这种描述适用于各向同性的表面结构，虽然不少研究者得出的一致结论认为天然土层是横观各向同性的[10]，但实际取样过程中不可避免的误差会对实验的结果造成不可预测的影响。

计算机断层扫描技术(Computed Tomography,CT)作为一种无损检测技术，近年来已推广至岩土工程界。相比SEM技术只能观测样品表面结构形态并在制样过程中可能对样品表面的损坏等方面的缺陷，CT技术可以无损且连续地观测黄土内部的微结构特征。但到目前为止，诸多学者多是利用CT技术观测黄土样品在受压或加湿等条件下内部结构的变化，关于利用CT技术进行黄土内部结构定量化的工作鲜有报道。尽管CT照片的空间分辨率低于SEM照片[11-12]，且由于黄土内部结构复杂，各介质之间具有含糊性和不确定性，在生成的CT照片中无法准确区分黄土的微观介质与颗粒。但较大孔隙或裂隙却很清晰，而其对黄土的物理力学性质及工程都有较大影响，因此利用CT技术对黄土孔隙进行定量化研究具有积极的理论与实际意义。同时考虑到黄土内部结构图像中的三维信息与拓扑关系，通常的三维分形维数计算方法如计盒法、方格法、结构函数法等计算过程比较复杂，不便于实践应用。

本文利用CT技术，采用周长-面积法求取黄土样品中孔隙形态的分形维数，简化三维分形的计算过程，提高三维分形维数计算的工作效率。

2 研究方法

2.1 基本原理

面积-周长法起源于Mandelbrot提出的小岛法。对于规则图形来说，其面积A与测量单位尺寸ε成二次方比例，周长P与ε成一次方比例，周长和面积存在如下关系：

而对于二维空间中的不规则图形，周长和面积的关系要复杂一些，Mandelbort提出应该用分形周长曲线代替原来的光滑周长，从而得到下列关系式：

其中，a0是和小岛形状有关的常数；ε是测量尺寸。对上式取对数后得：

2.2 研究思路

黄土大孔隙的表面形态呈现为一个具有粗糙表面的不规则状复杂体，只具有统计意义下的自相似性，经过图像处理和人工选定的黄土骨架颗粒或孔隙可以理解为边界闭合的图形。因此，可以通过不断变换测量尺寸ε来绘制以为横坐标，为纵坐标的图，将求解分形维数的问题变为在特定的测量尺寸下绘制不同大小孔隙的周长和面积的双对数图来求解分形维数的问题。

利用CT技术能以相同的距离连续拍摄断面图像这一特性，可以得到一组一个孔隙从小变大再变小的连续照片，以Q3黄土为例，由于照片数量较多，不能一一呈现，因此选取有较大差异变化的照片进行说明。将这些等间距的断面图叠加在一起就生成了一幅三维模拟图，而将测量后得到的一组连续断面图中同一孔隙周长和面积的数据在双对数坐标下绘制成离散点数据，若这些数据存在直线部分，则其斜率的倒数即为该孔隙的三维分形维数D。应用此法求取多个孔隙的三维分形维数并取均值即可视为该区黄土的三维分形维数。

2.3 实现方法

由以上分析知，黄土大孔隙三维分形维数计算过程为：

Step1取得典型黄土样品的CT照片。

Step2在一组连续的CT照片中选取具有代表性的孔隙。要求孔隙的边缘清晰并在照片中呈现出从无到有再到无、由小到大再到小的变化趋势。

Step3选取特征孔隙后，计算孔隙的周长与面积。由于通过图像的灰度值进行编程计算算法较复杂，利用CAD软件(如南方CASS软件)可以容易地计算出指定孔隙的周长和面积，将这些离散数据绘制于双对数坐标下，若存在直线部分，取直线部分斜率的倒数即为该孔隙的分形维数。

Step4重复上述步骤，选取多个孔隙进行计算，计算结果的均值即该黄土的三维分形维数。

3 应用实例

3.1 实验设备与实验概况

实验所用的Q2和Q3黄土土样取自陕西西安临潼区骊山山麓凤凰池施工现场，由于该现场存在典型的新鲜黄土剖面，因此使用环刀在不大力作用下可在施工现场取样。为了使孔隙保持其原始状态，取样过程严格按照规范《原状土取样技术标准JGJ89-92》执行。在制作拍摄土样时采用刀切法制作立方体土样，便于进行三维扫描。同时由于本实验借助CT机观测土样内部结构，因此刀切法对土样表面颗粒孔隙形态的破坏不影响最终的分形维数统计结果。

通过土工实验，得到该区Q2和Q3黄土物理力学性质，如表1所示。由表1可以看出，该区的Q2黄土含水率、初始密度和c值均高于Q3黄土，Q2黄土φ值与Q3黄土相当。Q2黄土的孔隙比e小于Q3黄土，这主要是因为土体长期受上覆荷载，大孔隙被压实挤密所致。

表1 黄土物理力学性质与分维数对照

实验所用的CT机借助长沙理工大学MG325型CT设备，该设备实验时的加速电压是300 kV，电流是2 mA，电流焦斑尺寸是1.9 mm，射线焦斑尺寸是1.9 mm，CT机扫描层间距为0.2 mm。

3.2 分维值确定

根据以上分析，要计算如图1所示的Q3黄土图像中孔隙形态的三维分形维数。

图1 黄土不同层孔隙变化的CT照片

分维值的确定需要进行以下操作：(1)选取合适的孔隙，即在连续的CT照片中选取连续变化的孔隙。由图1可以发现，孔隙图像基本完整，边缘明显。(2)利用CAD软件的功能(以南方CASS7.1为例)，将照片导入软件中，计算选定孔隙的周长和面积。(3)将步骤(2)中得到的周长和面积数据统计之后拟合双对数图，求取拟合后直线段部分斜率的倒数即可得到该孔隙的三维分形维数。经过以上3步得到了Q3黄土一个孔隙在连续变化条件下的周长和面积，如表2所示。

表2 孔隙周长和面积计算值

将表2中的周长和面积在Excel中转换成双对数数据，并在Matlab程序中采用最小二乘法进行线性拟合，如图2所示。由图2可以看出，表2中的周长和面积数据在双对数坐标下显示出了良好的线性关系，这表明土体中的孔隙形态具有分形特征。计算直线部分的斜率并取其倒数即得到了该区Q3黄土孔隙形态的一个三维分形维数，D=1.327。

图2 孔隙周长和面积的双对数关系

重复上述步骤，分别选取并计算Q2和Q3黄土CT图像中多个孔隙的三维分形维数，将得到的数据求取均值即为该区Q2和Q3黄土的三维分形维数：Dq2=1.151,Dq3=1.352。

3.3 结果分析

发现该区Q2黄土孔隙三维分形维数小于Q3黄土。由于基于周长-面积算法的分形维数反映了土壤孔隙轮廓边界的曲折程度，孔隙轮廓越不规则则其分维数越大，因此由实验结果可以发现该区黄土Q3内部大孔隙较Q2黄土要复杂，孔隙轮廓更不规则。同时对比表1中的物理力学性质，发现土体大孔隙的三维分形维数与含水率、初始密度和内聚力c成反比例相关，与孔隙比e成正比例相关。

4 与基于SEM方法的比较

基于CT的周长-面积方法与基于SEM的周长-面积方法是以黄土孔隙轮廓线的粗糙程度为研究对象并对其进行定量化的2种不同的方法。这2种方法在计算分形维数时各有优缺点，简述如下。

由于SEM技术的局限性，该方法只能计算样品表面孔隙的分形维数，尽管基于SEM技术的黄土孔隙分形计算方法发展迅速，但仍局限于二维；CT技术可以对样品任意地切取连续的剖面，对连续剖面照片上的同一孔隙应用周长-面积法求取其分形维数，从而真三维地揭示了黄土孔隙的空间分布和规律。

基于SEM的周长-面积方法对土样制备过程以及SEM图像效果要求都比较高，尤其要保证土样表面未经扰动；基于CT的周长-面积方法可连续对土样切取剖面，因此，土样表面的轻微扰动对于图像生成和实验结果并无影响。

SEM技术的分辨率高于CT技术，借助SEM技术可以观测到土体内部细小的孔隙结构。但由于黄土微结构图像的复杂性和不确定性，在使用周长-面积法计算孔隙的分形维数时，需要对图像进行二值化处理，选择不同的阈值，得到的二值图像亦不相同，因而得到的分形维数也不相同，即阈值与分形维数存在一定的关系，由于阈值选择不同或不恰当，易造成计算结果的差异进而影响分析结果。尽管CT图像空间分辨率低于SEM图像，但由图1可以看出，黄土CT图像中的孔隙与周边介质边缘明显，这有助于准确计算孔隙的周长面积比和分形维数。

5 结束语

本文基于CT技术的黄土大孔隙三维分形维数定量化方法，通过CT扫描图片，计算黄土大孔隙的周长和面积，得出黄土大孔隙的三维分形维数，可以定量地从真三维的角度准确描述黄土内部大孔隙在不同载荷作用下的变化情况。该方法不但避免了制备土样时对土样表面孔隙形态造成的破坏，而且孔隙分形维数的计算结果不随阈值的变化而变化，加之测量软件的发展和应用大大简化了孔隙周长和面积的计算过程，显著提高了计算分形维数的工作效率，因此，该方法具有较大的理论和现实意义。

由于CT技术空间分辨率的局限性而不能观测到土体内部更细微的结构，且没有定量地分析出黄土大孔隙三维分形维数与黄土宏观物理力学性质之间的关系，下一步工作将重点研究两者之间的关系。

[1]王念秦,张倬元.黄土滑坡灾害研究[M].兰州:兰州大学出版社,2005.

[2]姚 涛.基于三轴土样变形数字图像测量的黄土变形和强度研究[D].大连:大连理工大学,2008.

[3]Terzaghi K.Theoretical Soil Mechanics[M].[S.l.]:Weley,1943.

[4]岳中琦.岩土细观介质空间分布数字表述和相关力学数值分析的方法,应用和进展[J].岩石力学与工程学报,2006,25(5):875-888.

[5]邵生俊,周飞飞,龙吉勇.原状黄土结构性及其定量化参数研究[J].岩土工程学报,2004,26(4):531-536.

[6]Tovey N K.A Digital Computer Technique for Orientation Analysis of Micrographs of Soil Fabric[J].Journal of Microscopy,1990,120(3):303-305.

[7]任 权,王家鼎,袁中夏,等.高速铁路地基黄土微结构的分形研究[J].水文地质工程地质,2007,34(6):76-78.

[8]李 敏,李 毅.土壤颗粒数量分布的局部分形及多重分形特性[J].西北农林科技大学学报:自然科学版,2011,(11):224-230.

[9]范卫峰,林启太.基于分形理论的岩石质量评价[J].矿业研究与开发,2005,25(3):31-33.

[10]龚晓南.土塑性力学[M].杭州:浙江大学出版社,1990.

[11]朱元青,陈正汉.原状Q_3黄土在加载和湿陷过程中细观结构动态演化的CT-三轴试验研究[J].岩土工程学报,2009,31(8):1219-1228.

[12]王慧妮,倪万魁.基于计算机X射线断层术与扫描电镜图像的黄土微结构定量分析[J].岩土力学,2012,33(1):243-254.