RVE孔隙模型細觀結構特征分析與對比

楊鵬飛，李 顯，闕 云，蔡沛辰

(福州大學 土木工程學院，福州 350108)

1 研究背景

我國東南沿海地區(qū)，花崗巖殘積土分布非常廣泛，其不同于一般黏性土，屬于典型大孔隙結構[1]。除花崗巖殘積土自身孔隙結構外，在自然氣候等條件影響下，斜坡表面土壤內分布著許多孔隙，這將進一步改變土體孔隙結構特性[2]。在降雨等作用下水流可通過孔隙繞過大部分基質土壤，產生孔隙流，導致邊坡失穩(wěn)[3]。因此，如何準確刻畫土壤內孔隙分布情況是揭示土壤邊坡失穩(wěn)機制的基礎與核心內容之一。

不同學者對于界定大孔隙的標準存在一定差異，石輝等[4]定義大孔隙為等效直徑在0.2～3.0 mm范圍內的孔隙；時忠杰等[5]定義大孔隙直徑為0.4～2.3 mm；Marshall[6]將直徑>30 μm的孔隙定義為大孔隙；Jarvis[7]將直徑在0.3～0.5 mm的孔隙定義為大孔隙。本研究在前人的基礎上，考慮工業(yè)CT掃描儀的精度，將等效直徑大于最低分辨率0.15 mm的孔隙定義為大孔隙。

目前，土壤大孔隙的研究方法多為染色示蹤法[8]、穿透曲線法[9]等，但這些方法對觀測的土體試樣破壞性較大，無法得到大孔隙在土壤中的三維結構參數。相比之下，CT掃描技術對原狀土大孔隙研究更加便捷、精確，可在不破壞土體原始結構的條件下，提取大孔隙的三維結構特征參數。Anderson等[10]最早利用CT掃描儀對干土和濕土進行掃描；馮杰等[11]采用CT掃描得到了大孔隙在截面上的分布情況；薛華慶等[12]采用微米CT掃描表征了頁巖的發(fā)育情況；戚楠[13]采用不同的孔隙結構表征技術對儲層孔喉進行了定量表征和二維平面刻畫。但通常，對大孔隙空間結構參數的計算往往還需將掃描圖像作進一步處理，常用的圖像處理軟件有AVIZO、Image J等?，F(xiàn)有的大孔隙結構研究已取得了一定進展，但由于量化技術、地域性差異等原因還存在一定問題。很多學者對大孔隙的研究多集中在二維結構參數上，且研究多針對農業(yè)土壤成巖石結構，很少有系統(tǒng)討論三維結構特征參數的影響，并且在不同的研究條件、研究區(qū)域的影響下得出的結果也不盡相同。同時，現(xiàn)有的大孔隙結構參數獲取過程中存在計算機運算速度慢、所需內存容量大等特點。

鑒于此，本文從三維層面分析原狀花崗巖殘積土大孔隙分布特征，基于工業(yè)CT掃描圖像和AVIZO軟件構建了孔隙網絡模型，選取代表性體積單元(Representative Volume Element，RVE)，對花崗巖殘積土大孔隙三維細觀結構進行了定量表征，并對比了不同孔隙處理方法下的孔隙結構特征。研究成果可為短時強降雨下含大孔隙斜坡的水分運移機制研究奠定基礎，并可進一步豐富短時暴雨型滑坡理論基礎及支撐災害防治。

2 試樣與方法

2.1 原狀土取樣

原狀土取自福州某地山坡，分別選取兩個試驗區(qū)A和B，沿深度范圍進行取樣(A區(qū)試樣1和B區(qū)試樣3在深度20 cm處取樣，A區(qū)試樣2和B區(qū)試樣4在深度70 cm處取樣)，最終獲取4個原狀土樣，如圖1[14]所示。 測試土樣的基本物理參數，如表1[15]所示。

圖1 原狀土取樣[14]Fig.1 Sampling of undisturbed soil

表1 基本物理參數[15]

2.2 CT掃描試驗

在原狀土大孔隙結構特征觀測方法中，CT掃描屬于無損檢測方法，可在不破壞土體結構的情況下獲得土體內部孔隙結構圖像。其他學者研究的結果表明，在掃描土壤橫斷面圖像質量方面，工業(yè)CT的圖像質量要明顯高于醫(yī)學CT圖像[16]。在三維重構方面，工業(yè)CT解析出的三維圖像能更好地體現(xiàn)出土壤孔隙分布。鑒于此，對獲取的原狀土樣進行工業(yè)CT掃描，設備名稱為C450 kV高能量工業(yè)CT，將掃描試樣放置在該CT系統(tǒng)工作臺的中央，工作電壓為450 kV，電流為63 mA，掃描最低分辨率為150 μm。得到典型CT掃描圖像如圖2所示。

圖2 典型CT掃描圖像Fig.2 Typical CT scanning images

2.3 三維重構預處理

本次研究采用VG studio MAX 3.0軟件對CT圖像進行三維重構，并結合AVIZO軟件進而提取孔隙面積、孔隙等效半徑、孔隙配位數、喉道面積、喉道等效半徑、喉道長度等參數。

由于收集到的CT圖像與真實情況有一定的差距，因此在三維重構前需要對CT圖像進行歸一化和環(huán)狀偽影去除等操作。其中歸一化主要用來糾正掃描過程中的光束強度變化，能對不同時間收集的CT投影圖像的強度進行補償，使強度的變化充分均勻。環(huán)狀偽影是CT投影成像中較為常見的一種現(xiàn)象，使用平面場校正功能可以有效去除環(huán)狀偽影。

CT掃描得到的原始切片需后續(xù)進一步處理才能進行細觀結構識別。VG studio是目前最常用的圖像三維重構軟件之一，主要功能為顯示、處理和分析工業(yè)CT掃描后得到的圖像數據。由于本次試驗著重于研究花崗巖殘積土大孔隙的三維形態(tài)，在較高的分辨率情況下，細微的擾動對試驗結果可能造成較大的影響，因此需要選擇適當的分析區(qū)域，去除各類干擾，盡可能還原花崗巖殘積土大孔隙的三維形態(tài)。樣品尺寸為15 cm×15 cm×40 cm，選取的分析區(qū)域為13 cm×13 cm×36 cm。

2.4 三維模型構建

CT圖像三維重構是采用連續(xù)的CT圖像通過插值運算生成三維立體結構的一個過程[17]。本次試驗根據文獻[14]、文獻[18]、文獻[19]中的方法，選取了300×300×300體素的代表性體積單元(RVE)，其中體素是數字數據于三維空間分割上的最小單位。相鄰切片之間的間距與CT圖像的分辨率相同，即每隔0.15 mm掃描一張CT圖像。

AVIZO可視化軟件是一套全方位的三維數據處理軟件，可對三維圖像進行閾值分割、曲面重建、網格劃分、孔隙網絡模型建模等操作[14]。對圖像進行三維分析時，采用擬合法[20]，在每個試樣中隨機挑選10張CT圖像導入Image pro-plus軟件，然后利用Image Histogram功能提取出圖像灰度值，繪制灰度值分布圖，并通過軟件進行擬合，得到每張圖的閾值，對閾值求平均值，即可得到4個試樣的閾值分別為15 470、16 505、19 920和16 528。

圖3為模型中間位置處4個試樣三維圖像采用中值濾波器降噪后的閾值分割圖，其中藍色部分代表孔隙域，灰色部分代表基質域。將分割后的三維圖像進行對象分離后，對孔隙結構進行分析，即可得到三維孔隙的特征參數。圖4為4個原狀試樣的代表性體積單元三維重構結果，圖5為刪除孤立孔隙后的簡化孔隙網絡模型。

圖3 原狀土閾值分割Fig.3 Threshold segmentation of undisturbed soil samples

圖4 原狀土三維重構模型Fig.4 Three-dimensional reconstruction model of undisturbed soil

圖5 簡化孔隙網絡模型Fig.5 Simplified pore network model

3 結果與討論

3.1 孔喉特征參數分析

采用孔隙模型分析大孔隙結構特征時，主要將孔隙空間劃分為孔隙結構和連接孔隙之間的喉道結構。因此孔喉的尺寸特征分析主要包括孔隙數量、喉道數量、孔隙結構參數和喉道結構參數。孔隙特征參數見表2，喉道特征參數見表3。

從表2可以看出：不同試樣的連通孔隙數量差值較大，試樣3的孔隙數量是試樣1的2倍左右，而不同試樣的孔隙半徑差值較小。圖6為4個試樣的孔隙平均半徑分布，可以看出4個試樣連通孔隙的半徑多介于1～4 mm之間。

從表3可以看出：4個試樣的喉道數量差距較大，其中試樣3的喉道總數為3 765個，而試樣1的喉道總數僅有819個，說明不同試樣之間的孔隙連通情況差異性較大。4個試樣的喉道等效半徑、喉道長度和喉道表面積平均值差異性較小，表明在花崗巖殘積土中雖然孔隙數量差異性較大，但由于孔隙形成原因類似，孔隙結構之間差異性較小。

由2.4節(jié)可知，本文REV模型尺寸為300×300×300體素，而文獻[14]中為大尺寸三維重構模型，其相比于小尺寸模型，能夠更真實地反映土體的原始結構，但同時計算量更大，運算也更加耗時。再對比二者的孔隙特征參數結果發(fā)現(xiàn)，4個試樣的孔隙平均半徑、喉道長度和喉道表面積相差較小(平均誤差分別為4%、0.8%和2.9%)，孔隙表面積和喉道等效半徑平均誤差在5%和17%左右，分析原因可能是原狀土孔隙三維結構較為復雜，選取區(qū)域與整體區(qū)域結構特征存在不同。此外，由于REV選取的自身特點使得二者在孔隙配位數方面差異性較大[14]，但總體REV選取的方法仍可在一定程度上對孔吼結構參數進行分析。

圖6 孔隙半徑分布Fig.6 Distribution of pore radius

3.2 孔隙拓撲結構量化表征

土壤孔隙空間的拓撲結構是指土壤內部孔隙與喉道連接關系，即連通率[21]?？紫锻負浣Y構量化可由孔隙配位數來表征。由表2可知，4個試樣的孔隙配位數最大值分別為24、31、47、38，經過統(tǒng)計發(fā)現(xiàn)花崗巖殘積土中孔隙配位數主要集中在20以內，超過20的孔隙配位數分布呈離散點狀，且數量皆為個位數，其與大尺寸模型統(tǒng)計結果較為相近[14]。因此本文在研究孔隙配位數分布情況時，只統(tǒng)計20以內的數值，并將其繪于圖7。

表2 孔隙特征參數

表3 喉道特征參數

從圖7可以看出：4個試樣的孔隙配位數分布差異性較大，其中試樣1中配位數為1～4的孔隙構成了土壤孔隙的主體，占總孔隙數量的65.52%；試樣2中配位數為1～4的孔隙構成了土壤孔隙的主體，占總數量的69.93%；試樣3中配位數為1～7的孔隙構成了土壤孔隙的主體，占總數量的91.73%，試樣4中配位數為1～6的孔隙構成了土壤孔隙的主體，占總孔隙數量的89.22%。上述現(xiàn)象表明：在花崗巖殘積土中，不同位置的土壤內部孔隙拓撲空間結構差異性較大，而深度對其影響較小，如試驗區(qū)B的不同深度試樣3和4的孔隙連通性都較好，而試驗區(qū)A中試樣1和試樣2的孔隙連通性都較差。

圖7 孔隙配位數分布Fig.7 Distribution of coordination number of pores

同時對比文獻[14]中大尺寸模型，小尺寸試樣1、2、4中構成土壤孔隙主體的配位數差異不大，而試樣3中獲得的配位數1～7的孔隙占孔隙總數量比例偏大，分析原因可能是由于REV模型試樣3相比大尺寸模型孔隙數量較少，整體分布不均勻所致。綜上，REV模型可一定程度上表征孔隙的空間拓撲結構。

3.3 不同處理方法孔隙參數的差異性分析

為對比不同處理方法下大孔隙特征參數的差異性，本文將二維圖像處理[22]和三維重構處理分別提取出的孔隙參數平均值匯總于表4。為減少工作量，二維圖像每隔1.5 cm選取一張進行孔隙參數分析，三維重構處理可分為孔隙模型和孔隙網絡模型兩種。

表4 2D和3D孔隙參數平均值

從表4可以看出：采用二維圖像處理得到的孔隙平均直徑最小，而且不同試樣之間差異性較大，而通過刪除孤立孔隙而建立的孔隙網絡模型得到的孔隙平均直徑最大。同時，對比三維孔隙模型和孔隙網絡模型處理結果可知，在孔隙網絡模型中通過等效半徑計算出4個試樣的平均孔隙體積分別為33.36 、42.68 、47.27 、49.18 mm3，通過孔隙模型計算出4個試樣的平均孔隙體積分別為6.28、6.05、5.57、8.34 mm3，說明在土壤內部有大量的孤立孔隙存在，由于孔隙體積越小，和其他孔隙產生連通的幾率也越小，因此孤立孔隙的體積一般遠小于連通孔隙。當采用孔隙網絡模型分析孔隙參數時，孤立孔隙被刪除，因此造成2種統(tǒng)計方式平均孔隙體積相差5倍以上?？紫兜钠骄砻娣e和平均孔隙體積類似，在三維孔隙模型中4個試樣的平均孔隙表面積分別為16.47、15.87 、14.80 、19.11 mm2，而通過孔隙網絡模型計算出4個試樣的平均孔隙表面積分別為136.38 、178.46 、249.98 、267.31 mm2，可以看出通過2種方式統(tǒng)計出的平均孔隙表面積差距比平均孔隙體積更大，原因在于隨著孔隙的增大，孔隙表面復雜程度也隨之增加。

本研究中對于大孔隙的分析采用的是三維孔隙網絡模型處理方法，原因在于：二維處理方法是對眾多CT掃描切片進行參數分析，只能從二維層面體現(xiàn)孔隙特征，而事實上，土體中的大孔隙是三維結構，二維表征的精確度有待商榷；三維孔隙模型處理方法雖是從三維層面對土體大孔隙結構進行表征，但統(tǒng)計時將各種孤立孔隙都計算在內，而這些孤立孔隙一般不參與滲流過程，對于后續(xù)滲流研究不能起到很好的指導作用；三維孔隙網絡模型處理方法是在三維孔隙模型的基礎之上，檢測連通孔隙和孤立孔隙，并將孤立孔隙刪除之后形成的連通孔隙網絡模型，這樣可以有針對性地統(tǒng)計有效孔隙結構特征參數，還可以提高后續(xù)滲流模擬研究中的效率。

4 結 論

(1)花崗巖殘積土RVE模型中不同試樣的連通孔隙數量差異性較大，而孔隙結構參數差值較小，且孔隙半徑多介于1～4 mm之間，表明在花崗巖殘積土中雖然孔隙數量差異性較大，但由于孔隙形成原因類似，孔隙結構之間差異性較小。

(2)不同區(qū)域的土壤內部孔隙拓撲空間結構差異性較大，部分孔隙的配位數可達40以上，但深度對其影響較小。其中孔隙配位數主要集中在20以內，超過20的呈離散點狀分布。

(3)對比不同處理方法知，2D圖像處理得到的孔隙平均直徑最小，孔隙差異性較大，而通過刪除孤立孔隙而建立的孔隙網絡模型，得到的孔隙平均直徑最大，孔隙分布較為均勻。對比3D處理結果表明，花崗巖殘積土內有大量的孤立孔隙存在，而孤立孔隙的體積一般遠小于連通孔隙。

(4)明確土壤內部孔隙特征分布是揭示土壤邊坡失穩(wěn)機制的基礎，本文對三維層面花崗巖殘積土內大孔隙分布特征進行了定量表征，考慮到土體原始三維結構計算復雜，采用REV孔隙模型代替大尺寸模型研究，其結果能較為精確地表征土壤孔隙結構，但由于土壤內部結構復雜，整體結構會與選取區(qū)域存在部分差異。REV模型可為小尺寸原狀土定量化研究孔隙內部水分運移機理提供一定借鑒。