基于CT的黃土孔隙尺度優先流特性

寧瑞浩， 冷艷秋*， 何芝遠， 李澤坤， 馬哲

(1. 長安大學地質工程與測繪學院，西安 710054； 2.西部礦產資源與地質環境教育部重點實驗室，西安 710054)

黃土在中國主要分布在山西、陜西、甘肅、寧夏的黃土高原一帶，是一種疏松多孔、豎向節理、裂隙發育，具有很強的濕陷性和水敏性的結構性土；因其特殊的濕陷性和水敏性，水便成為了黃土地區誘發地質災害的關鍵因素。常規的滲水試驗表明降雨或者灌溉的入滲深度只有2～3 m[1]，但是黑方臺和涇陽等地的滑坡實例卻顯示水運移到坡體以下數十米，優先流便是這類地質現象的合理解釋，因此探索黃土孔隙、裂隙中優先流的滲流特性，對防控黃土地質災害具有重大意義。

優勢流又稱為優先流，是一種非平衡流，主要特征為部分水流繞過大部分土體基質沿著優勢通道迅速入滲到土體深部[2]。目前對優先流的研究主要集中在兩個尺度，即宏觀尺度和孔隙尺度[3]，Watanabe等[4]通過構建局部非孔洞模型，研究了碳酸鹽巖的孔隙度與滲透率和內部優勢流的關系，為碳酸鹽巖儲層的開發利用提供了依據。趙寬耀等[5]通過物探探測與數值模擬相結合的方法，分析了黃土邊坡在優勢流與基質流共同作用下不同灌溉強度水的滲流過程。陳思婕等[6]將雙重滲透模型與無限邊坡穩定分析法相耦合，根據土壤含水量及孔隙水壓力觀測數據，模擬了強降雨條件下滑體內的水動力過程，以此分析了優先流對滑坡觸發機理的影響。然而這些研究都是集中在宏觀層面的。潘網生等[7]認為定量化描述黃土微細觀孔隙、裂隙結構模型,建立從微細觀到宏觀不同尺度下研究黃土滑坡機理的系統理論方法，對認識優先流對滑坡的觸發機理尤為重要，僅僅開展宏觀尺度的研究并不能揭示黃土中優先流的本質規律。

研究孔隙尺度內的優先流特性首先要獲取真實的黃土內部孔隙結構，在眾多方法中計算機斷層掃描(computed tomography,CT)技術作為一種無損檢測技術正被廣泛用于巖土材料的孔隙定量化表征。Li等[8]利用CT掃描技術定量化研究了馬蘭黃土的大孔隙結構特征，結果表明馬蘭黃土在垂直和水平方向上的孔隙形狀、聯通特征有很大差異，這與室內試驗與數值模擬結果較為吻合。王超等[9]采用微米CT建立了黃土巖的孔隙網絡模型，研究了黃土巖中礦物成分、形狀對黃土巖的滲透率各向異性的影響。Wang等[10]通過CT掃描技術研究了黃土區露天煤礦排土場重構土壤孔隙三維分布的多重分形表征，為礦區土壤重構提供了依據。上述成果說明利用CT技術獲取黃土內部真實孔隙結構并進行定量化的研究是可行的。

目前已有眾多學者在真實孔隙結構的基礎上開展優先流滲流特性研究，Li等[11]、Larsbo等[12]在孔隙尺度上研究了孔隙與孔隙網絡特征(平均直徑、配位數等)對土體內部優勢流滲流規律的影響。Kim等[13]通過數值模擬的方法研究了孔隙尺度內河床附近優勢流及湍流對溶質異常輸運特性的影響。Zhang等[14]采用低場核磁共振技術對納米顆粒在非均質多孔介質中的運移方式實時監測，觀察到了明顯的優勢流現象。魯拓等[15]研究了基于多種分形模型，計算了不同尺度范圍孔隙的分形維數，并討論了二者之間的聯系。Sarah等[16]利用高分辨率的微米CT研究了土壤內部的細微觀結構并利用貝葉斯統計方法探討了土壤的微觀特征與其飽和導水率等參數的相關性。Li等[17]利用AVIZO軟件研究了馬蘭黃土的大孔隙特征并基于形狀因子對其進行分類，為研究孔隙尺度內優勢流特征提供了依據。

前人所做研究采用CT掃描的分辨率多數大于50 μm，雷祥義[18]根據壓汞實驗的結果將大孔隙劃分為孔隙半徑為16～250 μm，中孔隙為孔隙半徑為4～16 μm，采用50 μm以上的分辨率顯然并不能滿足精度需求，因此在提取黃土內部真實孔隙結構的基礎上開展微細觀尺度的優先流滲流特性研究。采用高精度的微米CT，掃描分辨率為4 μm，以此構建真實孔隙結構，研究延安黃土孔隙尺度內的優先流滲流特性。

1 實驗材料與方法

1.1 實驗材料

實驗所用土樣取自陜西延安邊坡坡腳，取樣深度約為75 m，所取原狀土樣為圓柱體，直徑和高度分別為25 cm和50 cm，在密封和標記后運至實驗室，運輸過程中采取合理的避震措施從而避免擾動。通過對所取試樣隨機測量，研究區黃土試樣的基本物理性質見表1。圖1為試驗所采用黃土試樣的粒度分布曲線，其中粉粒(2～75 μm)含量81.72%，黏粒(<2 μm)含量為6.65%。結合粒度分布曲線及塑性指數，按照文獻[19]對黃土的分類，研究區黃土屬于黃土質砂黏土。

表1 土樣基本物理性質Table 1 Basic physical properties of soil samples

圖1 粒度分布曲線Fig.1 Particle size distribution curve

1.2 微米CT掃描

實驗所采用CT掃描設備為Phoenix v|tome|x s工業CT(圖2)，該CT系統組合有180 kV/15 W 高功率nanofocus X射線管和240 kV/320 W 的微焦點管，是一種靈活可靠的掃描分析工具。本次測試試樣為圓柱體，直徑7 mm，高15 mm(圖3)，掃描分辨率為4 μm，共獲取X、Y、Z方向切片14 200張，每張二維切片像素點的灰度值代表著該點的X射線衰減系數，密度越大衰減系數越大，借助這一原理再結合適當的分割方法，三維重建等手段，便可以把試樣內部的孔隙裂隙識別重建出來以用于二維、三位參數的提取，具體技術路線見圖4。

圖2 CT掃描設備Fig.2 CT scanning equipment

圖3 掃描土樣Fig.3 Scanned soil sample

圖4 技術流程圖Fig.4 Technical flow chart

1.3 室內滲透實驗

為研究黃土微細觀與宏細觀滲流之間的聯系并驗證計算結果可靠性，進行室內滲透實驗與數值計算結果進行對比。本次室內變水頭滲透實驗采用土樣為延安黃土原狀土樣與重塑土樣(圖5)，各自設置5組平行試驗，重塑土樣采用分層靜壓法制備，將原狀黃土進行碾碎、過篩后放入涂好凡士林的壓樣筒內，控制重塑土樣干密度與原狀樣相等，分五層壓實，層間接觸部位進行刮毛，每層壓實時間不少于45 min。試樣為標準尺寸，直徑、高度分別為61.8、40 mm，在真空缸中飽和后進行滲透實驗。滲透實驗流程按照《土工試驗方法標準》(GB/T 50123—2019)規定，待水頭管中水頭穩定后在起始時刻標定水頭高度，每隔20～40 s測定水頭和時間的變化，如此連續測記10次后再使水頭回升至適當高度，重復試驗5次，所得結果如表2所示。

圖5 實驗土樣Fig.5 Experimental soil sample

表2 室內滲透實驗所得滲透系數Table 2 Permeability coefficient obtained from indoor permeability test

2 CT圖像預處理與三維重建

2.1 建模軟件選取

在目前常用的多孔介質二、三維孔隙結構及顆粒定量化分析平臺中，AVIZO軟件可以滿足圖像降噪處理-圖像二值分割-三維體重建-數據分析-滲流模擬-絕對滲透率計算等各方面的需求，在科學和工業領域、生命科學和生物醫學、材料科學、土壤學、巖土領域有著廣泛的應用，能夠滿足本研究的相關需求，因此選用AVIZO平臺進行優先流的相關研究。

2.2 圖像預處理

在CT圖像的采集和獲取過程中，會因為采集設備和儀器的少量干擾和周圍環境因素的影響，使得采集的圖像或多或少的存在部分噪聲，為了進一步提高掃描圖像的質量，突出圖像中孔隙的邊緣輪廓特征，以便能夠準確分割目標孔隙和背景，提取出感興趣的孔隙部分，通常采用圖像噪聲濾波方法和邊緣增強進行預處理。

目前常用的濾波方法有高斯濾波、中值濾波、非局部均值濾波法等，圖6分別為采用上述濾波方法的預處理結果。由圖6可知，非局部均值濾波法幾乎能完全去除圖像中的高斯噪聲、脈沖噪聲等，同時還能夠很好地保留圖像細節，很大程度地減輕邊緣模糊[圖6(d)]；其他方法雖然能夠消除噪聲但是不能很好的保留圖片信息，致使邊緣模糊，因此選用非局部均值濾波法對所獲取的CT圖像進行降噪處理。非局部均值濾波法的基本思想為根據圖像的自相似性來計算鄰域像素權重，通過加權的形式將最相近的幾個像素塊中的中心點結合起來估計真實值，缺點是對計算機的計算能力要求較高；雖然采用的濾波方法在有效去除噪聲的同時能夠很好地保留圖像細節，但仍然存在一定的邊緣模糊現象，即黃土固相和孔隙之間交界處的過渡帶[圖6(e)]。采用AVIZO軟件中的unsharp mask模塊對圖像進行邊緣增強處理[圖6(f)]，消減固相和孔隙之間交界處的過渡帶，大大增加分割精度。

2.3 圖像分割

在黃土孔隙結構的三維重建及要素提取研究中，針對黃土CT掃描圖像的分割技術是整個研究的核心與關鍵點之一，其分割結果直接影響到后續孔隙結構三維模型重建的精確性。到目前為止還沒有任何一種完備的分割方法能夠適用各種應用環境，常用的圖像分割方法可以劃分為:基于閾值的分割方法、基于區域的分割方法、基于邊緣的分割方法以及基于某些特定理論的分割方法等。本研究采用基于閾值的分割方法。閾值分割法的基本原理為給定一個恰當的灰度閾值，圖像某點處的灰度值大于該閾值認為是土顆粒，小于該閾值則認為是孔隙，因此閾值的選取尤為關鍵，過大會導致過分割，過小則會導致欠分割；常用的閾值分割方法有最大熵法、Ostu法、矩量法等，對同一張灰度圖像采用不同的分割方法的局部放大圖如圖7所示，最大熵法與矩量法確定的閾值偏大，存在明顯的過分割現象，Ostu分割法的分割效果與實際情況最為接近，然而還存在輕微的過分割現象，需手動調試閾值。基于此，分割閾值采用Ostu法確定參考閾值，在AVIZO軟件中通過手動調整確定最佳閾值[圖7(d)]，對部分難以識別小孔隙、微小裂隙通過AVIZO軟件中的Top-hat算法進行識別[圖7(e)]，最終取二者并集[圖7(f)]作為最終孔隙分割結果。經計算分割后的微觀孔隙度約為19%，小于實際孔隙度，這是由于小于4 μm的孔隙無法被識別出來。

圖6 圖像預處理結果Fig.6 Image preprocessing results

2.4 孔隙結構三維重建

2.4.1 代表性體元確定

目前，基于CT掃描圖像三維重構的基本方法主要包括體繪制和面繪制兩種方法，采用可以避免圖像信息丟失，能夠保留數據完整性的體繪制方法對土樣中的孔隙網絡進行三維重構。考慮到計算機的計算能力有限并且必須兼顧樣品具有夠的代表性，需確定代表性體積單元(representative elementary volume，REM)[20]。本文研究確定代表性體積單元的方法為在土樣任意一點處取一系列不同邊長的立方體(圖8)，計算其孔隙度，當孔隙度趨于穩定時即認為該邊長的立方體為最小代表性體積單元。如圖9所示，在立方體邊長達到1 250 μm時，代表性立方體孔隙度已達到穩定值，即立方體邊長達到1 250 μm以上時即可認為該子體積具有代表性。

圖7 圖像分割結果Fig.7 Image segmentation results

圖8 REV選取示意圖Fig.8 Schematic diagram of REV selection

目前對優先流的常見分類為：裂隙流、管流、指流、側向流，指流和側向流分別發育在上細下粗和含有較大阻水塊體的土體中，而本次所取土樣土質較均一且不含巖石、樹根等阻水塊體。基于此，在掃描土樣中分別提取含裂隙(C1、C2，圖10)、含管狀通道(C3、C4，圖10)與不含優勢通道的(C5、C6，圖10)2 000 μm×2 000 μm×2 000 μm的立方體，排除孤立孔隙，在各立方體連通孔隙基礎上對比分析裂隙流、管流的滲流特性。

圖9 孔隙度與立方體邊長關系曲線Fig.9 Relation curve between porosity and cube size s

圖10 各代表性立方體Fig.10 Representative elementary volume

2.4.2 隙網絡模型的構建

孔隙網絡模型(pore network model,PNM)的基本假定為一個相互連通的孔隙網絡有多個孔隙組成，不同孔隙之間通過孔道連接，孔喉是孔道中橫截面積最小的地方。在AVIZO中將孔隙網絡模型用球棍模型來表示，分別代表孔隙與孔隙之間的聯通關系，通過構建孔隙網絡模型可以獲取孔隙之間的孔道長度、孔喉面積以及配位數等關鍵信息，孔隙網絡模型如圖11所示。

圖11 孔隙網絡模型Fig.11 Pore network model

2.5 重建結果分析

2.5.1 孔隙與孔喉分布特征分析

根據所建立孔隙網絡模型所得數據，孔喉和孔隙等效直徑的分布范圍基本一致，對孔隙和孔喉的等效直徑以10 μm為間隔進行直方圖統計，發現孔隙與孔喉近似服從對數正態分布。根據圖12，孔隙主要分布在孔徑為40～140 μm的區間內，占總數目的91.83%，孔喉主要分布在孔徑為10～50 μm的區間內，占總數目的85.07%。

圖12 孔隙與孔喉分布特征Fig.12 Distribution characteristics of pores and pore throats

2.5.2 孔隙連通特征分析

孔隙的連通性可以用配位數表征，一個孔隙的配位數是指與其他孔隙連接的孔喉的數量，配位數越多則代表著與周圍的孔隙連通性越好，根據構建的孔隙網絡模型，對獲取的孔隙配位數進行統計，發現配位數近似服從對數正態分布(圖13)，其值介于1～28之間；含優勢通道的立方體中平均值為7.17，眾數為6，不含優勢通道的立方體中平均值為6.15，眾數為5，這說明在含優勢通道的立方體中，孔隙的整體連通性要好于不含優勢通道的立方體。

圖13 孔隙配位數分布特征Fig.13 Distribution characteristics of pore coordination number

3 滲流模擬

3.1 基本原理

Navier-Stokes方程是用于描述不可壓縮黏性流體的動量守恒方程，其矢量形式為

(1)

式(1)中:ρ流體密度，kg/m3；V為速度矢量；P為流體壓力,kPa；f為單位體積流體所受外力，kN/m3；μ為流體動力黏度,Pa·s。對V有

(2)

根據式(2)將式(1)展開，得

(3)

一般條件下很難得到其精確解，考慮以下條件將Stokes方程簡化：流體為不可壓縮流體；流體為牛頓流體，動力黏度μ為常數；流體為穩定流，流速不隨時間改變；流體流動方式為層流；對黃土孔隙尺度內流動的水，上述條件均可滿足。

簡化后的Stokes方程為

(4)

為保證方程在整個體積的有效性，將其轉化為體積平均形式，使方程在體積空間內更加平滑，即

(5)

式(5)中:D為速度擾動場張量；d壓力擾動場向量；I單位張量。

通過求解系統體積v上速度擾動場張量的平均值來求得滲透率張量k，即

(6)

在AVIZO中用有限體積法求解方程組，方程在交錯網格上離散，假定體素為立方體，待求解壓力位于體素中心，而速度在體素面上分解，待求解量的時間導數趨于零時認為收斂，誤差計算方法為

(7)

式(7)中:n為當前迭代次數；?t為時間增量；c2為虛擬壓縮系數[21]。

對于本次計算，當誤差小于10-4時即認為收斂，當求得絕對滲透率后，滲透系數計算公式為

(8)

3.2 模擬結果

模擬滲流方向為自上而下，上下邊界為定水頭邊界，四周為不透水邊界，上下邊界的水頭差為0.002 m，即代表性立方體的高度，各立方體內部流速分布如圖14～圖16所示。根據圖14(a)(e)～圖16(a)(e)，在含優勢通道的代表性立方體中(C1、C2、C3、C4)高流速的滲流路徑沿著裂隙展布方向集中，裂隙、管狀通道內的流線更加順直平滑，滲流壓力大，而其周圍的大孔隙內滲流路徑則較為分散，滲流壓力也更小，壓力分布更加不均勻；而在不含優勢通道立方體中(C5、C6)，流線及壓力的分布較為均勻，流速也更慢，表現為均勻入滲,各代表性立方體滲透系數及絕對滲透率如表3所示。

表3 滲流模擬所得滲透系數Table 3 Permeability coefficients obtained from seepage simulations

圖14 含裂隙代表性體元模擬結果Fig.14 Representative voxel simulation results with fractures

綜合室內滲透實驗及數值計算結果可知重塑黃土由于內部結構被完全破壞導致孔隙之間連通性極差，滲透系數最小，要比原狀黃土小近一個數量級。不含優勢通道的代表性立方體滲透系數與原狀黃土室內滲透實驗所得滲透系數在同一個數量級上，含優勢通道的代表性立方體滲透系數要比原狀黃土大近一個數量級，這與前人所得結果一致[22]，說明優勢通道對黃土的滲透性起著主導作用，也間接證明了黃土微細觀尺度與宏細觀尺度滲流的統一性。

為詳細研究黃土內部優先流的分布及流速特征，借助前述對黃土內部孔隙進行概化的孔隙網絡模型將孔喉按流量大小進行展示[圖14(b)(e)～圖16(b)(e)]，通過孔喉截面面積及流量換算流速并分組統計，依據數值計算及室內滲透實驗所得滲透系數將孔喉按流速劃分為三部分，不流動-低流速孔喉(T1)：流速小于1 μm/s，多位于孔隙連通性較差部位，流速與重塑黃土滲透系數相當，處于同一數量級；均勻流速孔喉(T2)：流速介于1～100 μm/s，距離優勢通道較遠，但孔隙連通性一般部位，與原狀黃土及含優勢通道的代表性立方體滲透系數處于同一數量級；優先流發育孔喉(T3)：流速大于100 μm/s，多位于優勢通道及孔隙連通性較好部位。

圖15 含管狀通道代表性體元模擬結果Fig.15 Representative voxel simulation results with pipes

對三個流速段的孔喉數量進行統計，統計結果見圖17。三種立方體中均以均勻流速孔喉為主，其中無優勢通道立方體中的T2占比最高，含裂隙立方體與含管狀通道立方體中的T2通道數量幾乎相同。將含裂隙立方體與含管狀通道立方體進行對比分析發現前者T3占比為30%，其含量超過T1，同時也比后者T3占比高，因此含裂隙立方體的絕對滲透率比含管狀通道立方體大。將含裂隙與不含優勢通道的代表性立方體孔喉流速分布進行對比，發現裂隙發育的立方體中T1所占百分比明顯降低，T2數量約減少10%，而T3所占比例顯著增加，約增加15%，說明裂隙對整體的滲透特性影響為均勻促進型，即裂隙的存在會整體地提升立方體的滲透性；而管狀通道發育的立方體與不含優勢通道的立方體相比T1占比約增加5%，T2數量約減少10%，T3所占比例增加幅度較裂隙發育的立方體小，但管狀通道內的最大流速顯著高于裂隙內部，這說明管狀優勢通道對整體的滲透特性影響為滲流集中型，即滲流路徑沿著管狀通道集中，對非優勢通道在滲流過程中所發揮的作用有所削弱。在不含優勢通道的代表性立方體中仍有相當一部分孔喉流速大于100 μm/s，這說明局部聯通性較好的大孔隙也可形成優勢通道。對三種流速段孔喉直徑進行統計發現T1部分孔喉平均直徑為11.2 μm, T2部分孔喉平均直徑為19.4 μm, T3部分孔喉平均直徑為29.1 μm,由此可見孔喉流速與孔喉直徑關系為正相關。

圖16 不含優勢通道代表性體元模擬結果Fig.16 Representative voxel simulation results without dominant channels

圖17 各流速段孔喉百分比Fig.17 Percentage of orifice throat at each flow rate section

4 結論與展望

在對黃土內部孔隙進行三位重建的基礎上進行滲流模擬，結合室內滲透實驗，可得出以下結論。

(1)基于真實孔隙結構數值計算得到的滲透系數與室內滲透實驗所得原狀黃土的滲透系數在同一數量級上，證明了計算結果的可靠性。

(2)從壓力分布圖可知，裂隙中不僅有著更高的流速而且有著更高的壓力，這種壓力差會促進聯通孔隙向著微裂隙轉變，即優先流對黃土內部的微觀孔隙結構存在著一定的改造作用，但這種作用的大小依賴于邊界壓差。

(3)根據原狀黃土、重塑黃土滲透系數將孔喉按流速分組統計，發現裂隙與管狀通道對整體滲透特性的影響分別為均勻促進型與滲流集中型。

更高的分辨率意味著更小的試樣尺寸及更大的計算機能力要求，隨著未來CT掃描技術與計算機計算能力的提升，可以進行更大尺度更高分辨率的研究，從而更好的研究黃土微觀結構和宏、細觀結構之間的聯系。