基于X-ray CT的古土壤孔裂隙識別與表征

李 鑫, 盧玉東, 張曉周, 盧陽春, 楊亞慧

(1.長安大學 環境科學與工程學院, 陜西 西安 710054; 2.旱區地下水文與生態效應教育部重點實驗室, 陜西 西安 710054)

古土壤是黃土高原漫長的形成過程中在濕熱環境中形成的棕褐色沉積物[1]。由于其黏粒含量高[2]，往往在黃土地區滑坡地質災害研究中被當作隔水層[3-5]。按此觀點，上覆黃土水分下滲過程中，當遇到古土壤層時會形成局部飽水帶，使得黃土軟化，強度降低[6-7]，形成滑動面，致使上覆黃土沿著古土壤接觸面下滑[8-9]。但在對已經發生的黃土滑坡調查時，發現古土壤層往往被切穿，多數的解釋是因為滑坡發生過程中，由于重力作用，使得古土壤層被切穿。然而，由于古土壤的黏性較高，在外界干濕交替、凍融循環、震動等過程中容易形成裂隙。如果是因為古土壤中存在的孔隙、裂隙優勢滲流通道，使得水分穿過古土壤層，在其下形成更深層次的滑動面，也不無可能。因此，定量化研究古土壤孔隙、裂隙結構，有助于揭示古土壤中水分賦存及運移特征，對進一步開展黃土斜坡地質災害機理研究具有重要的實踐意義。目前在黃土—古土壤層序空隙微觀結構研究方面，有大量針對馬蘭黃土、離石黃土的研究成果。高國瑞[10]利SEM和X射線衍射分析了黃土骨架顆粒形態、連接與排列方式，并提出了影響濕陷性的12種組合類型；雷祥義[11]利用壓汞法將西安市附近黃土按照成因分為原生和次生孔隙，按尺寸分為特大孔隙(>250 μm)、大孔隙(16～250 μm)、中孔隙(4～16 μm)、小孔隙(1～4 μm)和微孔隙(<1 μm)；李曉軍[12]在國內較早利用CT對壓實黃土的微結構進行初步探討；谷天峰[13]通過計算馬蘭黃土荷載前后SEM圖像中孔隙微觀結構參數的變化，揭示了土體變形過程中孔隙及顆粒變化特征；陳瓊[14]利用氮氣吸附法研究了滑帶土微孔隙吸附特性；程亞楠[15]利用CT掃描的方法對黃土水力學參數進行了提取，利用編程建立了孔隙網絡模型，并進行了水力學預測；潘網生[16-17]利用SEM對涇陽南塬黃土滑坡上樣品進行掃描，獲取黃土微觀孔隙參數及非均質性和分形分維特征，并利用逾滲理論和網絡分析方法對黃土微觀滲流進行模擬；Li[18]利用SEM對馬蘭黃土進行掃描，并分析了黃土微觀孔隙特征與濕陷性的關系；Li[19]利用CT對馬蘭黃土大孔隙特征進行研究，并對黃土大孔隙滲流進行了模擬；延愷[20]利用CT對黃土二三維微觀結構進行表征，得到了團聚體中土顆粒三維微觀形態特征和接觸方式。另外，FIB-SEM[21]、核磁共振[22]等方法也開始應用于土體微結構的研究。上述方法基本都是針對黃土層序中馬蘭和離石黃土進行的研究，而針對古土壤孔隙裂隙特征和孔隙尺度水力特性的研究卻少見報道。為此，本研究利用工業微焦點X射線計算機斷層掃描技術(X-ray ICT)對涇陽南塬原狀S5古土壤樣品進行掃描，對古土壤中發育的孔隙、裂隙進行識別提取和定量化表征，為進一步進行古土壤水力學特性預測提供依據。

1 材料與方法

1.1 樣品采集與處理

試驗所用古土壤原狀樣品采集自涇陽南塬高莊鎮修石渡村取土場剖面上比較明顯的紅3條S5古土壤層，為減少采樣過程中對樣品的干擾，首先利用刀和鋸切割30 cm ×30 cm×30 cm的樣品，用保鮮膜包裝好后，利用模板箱搬運回室內。然后利用切割機進行加工，最終試驗掃描樣品尺寸為5 cm×5 cm×10 cm。樣品干密度ρd為1.62 g/cm3，比重Gs為2.76，孔隙比e為0.7，飽和稱重法測得孔隙度為41.18%，塑限16.18，液限32.39，黏粒含量33.80%，粉粒含量63.88%，砂粒含量2.32%，粒度不均勻系數Cu=30.11,曲率系數Cs=1.93，該古土壤為級配良好的不均勻土。

1.2 X射線計算機斷層掃描(X-ray CT)

利用X-ray工業CT設備(圖1)進行古土壤樣品掃描，原理是利用X射線穿過樣品中不同物質時具有不同的衰減系數，即致密組分衰減系數大，疏松組分衰減系數小，X射線穿過樣品，由探測器接收到衰減程度不同的X射線，將其轉換為電信號，在計算機上呈現為灰度值(0～255)不同的CT切片圖像[23]，其中孔隙的灰度值最小，在圖像上表現最暗，高密度物質灰度值最大，表現為最亮，中間灰度階則為致密程度不同的固體土顆粒，通過沿樣品豎向掃描可以獲得一系列連續切片(圖1)。本次CT掃描，試樣尺寸為5 cm×5 cm×10 cm，空間分辨率為102 μm，共獲得1 728張像素為1 024×1 024 pixel的二維圖像。

圖1 X-ray CT實驗裝置及掃描切片示意圖

1.3 CT圖像濾波處理與二值化分割

由于古土壤樣品在CT掃描過程中儀器系統運行和環境噪聲的影響，所得到的圖像上存在噪點。為了獲取高質量的圖像，圖像處理的第一步就是要對連續CT切片進行濾波處理，達到既消除噪點，又保留孔隙邊緣的目的。AVIZO軟件中提供了20多種濾波算法，按功能主要分為平滑濾波(如Median Filter，Bilateral Filter等)、邊緣檢測濾波(如Sobel Filter， Gaussian Filter等)，銳化濾波、頻域轉換濾波和灰度轉換濾波[24]5大類，其中前3類在CT圖像處理中較為常用。

本文經過對原始圖像(以第330張切片為例)進行眾多濾波算法(圖2只展示部分濾波結果，包括中值濾波、高斯濾波、盒式濾波、雙邊濾波)計算比對發現，雙邊濾波器(Bilateral Filter)在降噪和邊緣保留方面比其他濾波算法效果好，主要是因為雙邊濾波在對周邊像素灰度值進行加權平均時，所采用的權重值綜合考慮了像素歐幾里得距離和像素局部灰度差異，使其同時具備降噪平滑和邊緣保留的功能[25]。利用該濾波得到了既可以有效降噪，又能很好地凸顯孔隙邊緣的圖像，為下一步進行圖像準確分割提供了良好的基礎。

圖像閾值分割的主要目的是將樣品中空隙空間和黃土顆粒分割開。本文利用閾值分割法(Threshold tool)和頂帽分割(Top hat)相結合，較為精準的區分出空隙和固體顆粒。其中，Threshold tool主要是基于影像灰度進行分割，但是分割過程中會出現閾值選取過小時，某些閾值邊緣灰度較淺的孔隙被分割為顆粒，閾值過大時，又出現過度分割的現象；Top hat分割方法是基于黑頂帽或白頂帽算法計算相鄰像素間灰度對比值來找出更“黑”或者更“白”的聚集像素，以識別出孔隙和土壤顆粒。Tophat識別出的空隙可以加入或者減除Threshold tool分割得到的孔隙或者顆粒中，以彌補或減輕因孔隙識別不足或者過度分割造成的偏差。通過兩種工具的聯合使用，對雙邊濾波法得到的灰度圖像進行分割，小于設定閾值的部分為孔隙，賦值1以黑色表示，大于該閾值的部分為顆粒，賦值0以白色表示，由此獲取了較為精準的二值化圖像，如圖2所示。

1.4 三維重建及可視化

基于移動立方體(MC， marching cubes)算法[26]，對上述閾值分割的系列二值化圖片進行三維重建，得到孔隙和古土壤顆粒三維模型，為了更加直觀地觀察孔隙形態和三維展布情況，利用分水嶺算法對三維孔隙體進行分割。基于光的發射和吸收原理，模擬預設光源通過孔隙體素并產生不同程度的發射和吸收，從而呈現出不同明暗程度的映射色系，達到可視化的目的。通過三維渲染效果圖(圖3)，可清晰地看到孔隙、裂隙的幾何形態、空間展布、連接方式、連通程度等信息。

圖2 圖像濾波處理及二值化分割

圖3 古土壤孔裂隙結構及基質顆粒

1.5 古土壤空隙參數計算及分類

為了對古土壤空隙結構進行定量表征，在圖像分割及三維重建的基礎上，首先計算了孔隙度(n)和形狀因子(SF)，并利用三維空隙形狀因子對古土壤中各類空隙進行分類，在分類的基礎上，計算了不同類型空隙的等效孔徑(d)。

孔隙度n為空隙體積(某切片內空隙面積)與土體(某切片)總體積(總面積)之比，見公式(1)：

(1)

式中：AV，AT——二維層面空隙面積和總面積(μm2)；VV，VT——三維空隙體積和總體積(μm3)。

形狀因子(SF, shape factor)是表征三維空間幾何體相對于球體褒癟長短的形狀特性，見公式(2)：

(2)

式中：A3D——三維空隙表面積(μm2)，V3D——三維空隙體積(μm3)。

等效孔徑d是表征古土壤中空隙尺寸的參數，由于樣品中空隙類型多樣，而AVIZO軟件提供的等效孔徑計算方法是利用球形直徑計算的，對非球形空隙 (如裂隙、枝杈狀孔隙、長柱狀孔隙)進行計算，難以真實準確表征所有空隙孔徑特性。實際上，對于裂隙來說，等效孔徑代表裂隙開度，在計算時應利用垂直于最大費雷特徑的寬度計算；對于枝杈狀孔隙來說，實際上是由若干長柱狀孔隙交叉連通形成的，而長柱狀孔隙局部為圓柱體，因此枝杈狀孔隙和長柱狀孔隙的等效直徑應利用圓柱體體積公式來推求；對于橢球和似球體孔隙則可按照球形直徑進行計算。本文對AVIZO中算法進行改進，針對不同類型空隙，應用不同算法來求取等效直徑〔公式(3)〕，以更加準確地表征空隙特征：

(3)

式中：B3D——裂隙開度(μm)；V3D——孔隙體積(μm3)；L3D——孔隙長度(μm)。

2 結果與分析

2.1 二維層面結構定性分析

從圖2中古土壤原始二維切片圖像可以看出，圖像以0～255的不同灰度組成，可以明顯辨識出3類組分，其中黑色為孔隙、白色為高密度物質(主要為鈣質成分)、中間灰色部分代表土顆粒。從圖2中二值化分割結果可以看出，古土壤中孔隙和裂隙共存，呈現隨機分散分布，形狀各異，有狹長延伸型、似橢圓或圓型、不規則型孔隙，局部連通性較好，但是從整體上看，大部分都是孤立存在。

2.2 三維孔裂隙結構宏觀定性分析

由于古土壤濕熱氣候中受環境改造明顯，且黏性較高，質地較黃土脆硬，容易產生裂隙及孔道，孔裂隙相連接共同構成空隙網絡。從三維數字模型中(圖3)可以看出，古土壤中空隙在形態上主要由裂隙、枝杈狀孔道、柱狀孤立孔隙、橢球及似球狀孤立孔隙組成，其中裂隙和枝杈狀孔道連通性較好，為古土壤中水分運移提供優先通道，是古土壤中產生優先流的重要原因，很大程度上決定了古土壤的宏觀滲透性質；長柱狀孔隙雖然分散孤立存在，難以導水，但是在水分入滲時，裂隙和連通孔隙拓展發育過程中，可能會將該部分孔隙貫通，形成更大的范圍的水力聯系；孤立的橢球或球狀孔隙通常為儲氣空間，一般起到阻水的作用，在極限條件下可能會以逾滲的方式發生水分交換。從發展演化的觀點來看，裂隙發育部位在樣品中上部，總體沿垂向延伸，開度隨深度減小，表明裂隙與特大孔隙是孤立大孔隙拓展的驅動條件，后者是前者產生更大水力貢獻的基礎。通過以上定性描述，大致了解了古土壤中孔隙裂隙形態、賦存方式，為認識孔隙裂隙結構提供了良好的基礎。但是要深入研究古土壤內部孔裂隙結構，還需要對其結構參數進行定量化表征。

2.3 孔隙度定量表征

在圖像分割的基礎上，由式(1)分別計算了Z，Y，X三個方向上的二維孔隙度值，得到了3個方向上連續切片的孔隙度分布曲線(圖4)，以及三向孔隙度基本統計表(表1)。從圖4和表1中可以看出，總體上，3個方向上孔隙度平均值為9.68%，與三維孔隙度相等；XZ和YZ平面上的孔隙度曲線形狀和分布趨勢相近，部分切片段內YZ孔隙度波動趨勢比XZ滯后，表明有傾斜孔隙存在；XY平面上的孔隙度值與前兩者差異較大，說明孔隙展布在三維空間上具有明顯的垂向性或近垂向性。由表1知，XZ和YZ平面變異系數較接近，都小于10%，XY平面變異系數是前兩者的2倍之多，由于樣品中孔隙和裂隙并存，沿Z方向(XY平面)的孔隙多寡變化強烈，造成XY平面變異系數較大，以上現象都說明該古土壤空隙具有較強的垂向性和空間異質性。

圖4 X，Y，Z三向切片孔隙度分布曲線

平面孔隙度/%平均標準差最小值最大值變異系數/%XY9.681.866.6714.1619.20XZ9.680.908.1112.059.27YZ9.680.756.8111.987.78三維9.68

2.4 空隙類型定量表征與分析

由圖3直觀地看出古土壤中空隙的形態，但是要將每一類提取出來，分類表征，需要借助三維形狀因子加以區分，因此，依據公式(2)計算了所有空隙的形狀因子，并利用自然斷點法，將系列值分成了5類(圖5)，即將前述五類空隙分別提取出來。當形狀因子SF>40.18為裂隙，7.01

圖5 基于形狀因子的不同類型空隙分布頻率及累積頻率

2.5 孔徑分布定量表征

在空隙類型劃分的基礎上，利用公式(3)計算了不同類型空隙的等效直徑，得到古土壤孔徑分布信息，并繪制了各孔徑區間孔隙數量及其體積百分比柱狀圖(圖6)。統計知,通過CT掃描識別出的最小孔徑為0.109 mm，最大為5.059 mm，平均值為0.477 mm。從圖6可以看出，樣品中孔徑主要分布在0.1～0.7 mm之間，該部分孔隙占到總孔隙數目的82.90%，0.7～1.0 mm占13.98%，1.0～5.0 mm占3.09%，大于5.0 mm的占0.02%，然而體積百分比方面，0.1～0.7 mm區間上孔隙體積占總孔隙體積的26.9%，0.7～1.0 mm占42.6%，1.0～5.0 mm占26.6%，大于5.0 mm的占3.9%。該結果與前述基于形狀因子分類得到的孔隙分布和體積占比結果一致，說明基于形狀因子對空隙類型進行分類表征可以達到對古土壤孔隙定量表征的目的。

圖6 孔徑分布與對應體積百分比直方圖

3 討 論

3.1 古土壤孔裂隙識別的水力學意義

通常認為，古土壤相對于黃土具有較小的滲透性能[4-5]，對垂向入滲的水分產生明顯的阻滯效應，并引起水分的側向滲透，形成局部飽和帶，改變滲流場，降低坡體強度[16]。上述認識是基于野外和室內滲透試驗建立的，滲流特性也是基于宏觀達西定律或Richard方程描述的，如何在孔裂隙尺度上對古土壤水分入滲規律進行精細刻畫，對于揭示古土壤在斜坡地質災害中水力學貢獻至關重要。本文無論從二維切片分割(圖2)還是三維孔裂隙模型結果(圖3)，都可以看出古土壤中發育規模不同、產狀各異、形態多樣的復雜空隙，具有較強的非均質性和各向異性，說明利用適用于均質假定的達西定律描述其滲流規律是有偏差的，因此有必要在孔裂隙尺度(微米尺度)上開展研究。本文采用新技術新方法對空隙類型進行了有效劃分，并提出各類空隙的度量算法，得到準確的孔徑分布結果。有助于進一步開展孔隙網絡建模、孔隙裂隙尺度滲流模型的建立與滲流模擬，從而從本質上揭示古土壤在斜坡地質災害中的重要作用。

3.2 分辨率對CT法測定孔隙度的影響

本文從CT掃描連續切片三維重建得到的孔裂隙結構，計算得到古土壤的孔隙度為9.68%，而利用飽和稱重法計算的孔隙度為41.18%，壓汞法計算得到的孔隙度為30.2%，3種方法測得的孔隙度值差異較大，主要原因是:飽和稱重法測定時，由于土顆粒的親水性，水分子直徑0.4 nm，可以進入所有孔隙，因此得到的孔隙度最大，壓汞法最佳測量范圍為0.003～400 μm的孔隙[27]，對于特大孔隙測量不準確，因此比稱重法偏小。而本文CT法測得的空隙受本次掃描分辨率影響，測得的最小孔隙為109 μm，最大孔隙達到5 mm，得到的孔隙僅僅為大孔隙和特大孔隙，因此可以看出，古土壤中等效直徑小于100 μm的空隙占總空隙的比例較大。

在空隙識別與分類中，由于本次掃描分辨率僅為102 μm，對于微小孔隙難以識別，但是可以識別出絕大多數有水力學意義的大孔隙(d≥100 μm)[28]，為了系統性的對古土壤空隙特征精確識別，在今后應開展多尺度、多分辨率CT掃描，從微—介—宏觀系統精確表征古土壤的內部空隙特性，為進一步研究古土壤孔裂隙雙重介質滲流及其水力學特性研究提供科學基礎。

4 結 論

本文為了進行古土壤孔隙裂隙結構識別與表征，優選雙邊濾波法對古土壤CT掃描圖像進行濾波處理，進一步結合閾值分割和頂帽分割較為精準地區分了空隙和固體顆粒，利用移動立方體(MC)算法重建了三維孔裂隙模型，經過光學渲染得到了古土壤內部三維可視化結構，通過定性分析與定量表征，對于空隙二三維形態特征、連接方式、空隙類型與孔徑分布規律等進行了分析。主要研究結論為：

(1) 依據二三維定性分析可知，古土壤孔隙裂隙結構并存，整體上呈分散分布，局部孔裂隙連通性較好；在形態上識別為裂隙、枝杈狀孔隙、長柱狀孔隙、橢球及球狀孔隙；空間展布上具有明顯的垂向性和空間異質性。

(2) 根據定量表征結果可知，古土壤大孔隙度(>100 μm的孔隙度)為9.86%，僅占總孔隙度(飽和稱重法孔隙度41.18%)的23.4%，說明古土壤以等效直徑在100 μm以下的孔裂隙為主；形狀因子是有效提取與識別孔裂隙結構的關鍵參數，依其可對空隙進行精確分類表征；可識別孔徑分布在0.1～5.0 mm之間，以0.1～0.7 mm孔隙居多，而0.7～1 mm孔隙體積貢獻率最大。

(3) CT無損掃描及三維可視化重建技術可以有效識別與表征古土壤內部孔裂隙結構，為進一步開展孔隙尺度滲流模擬提供可靠基礎。