基于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.