基于CT掃描的草炭土孔隙結構分析及滲流模擬*

楊圣濤 呂 巖 賀元源 劉婷婷 馬曉禎

(吉林大學建設工程學院， 長春 130026， 中國)

0 引 言

草炭土主要分布于我國的東北、西南、長江中下游地區(佴磊等， 2012)，主要發育在地表低洼，積水較多的沼澤類環境。草炭土是植物殘體經氧化和不完全分解作用堆積而成的一種特殊的腐殖質土，具有特殊的物質組成與結構特征，在自然環境狀態下，草炭土表現出高滲透性、高孔隙比、高含水量等不良工程地質特性(徐燕等， 2011)。近年來，修建在草炭土分布區的路基大多出現了不均勻沉降的現象，這些現象的存在影響了當地建筑工程的安全性、線路工程的交通營運(劉飛， 2011)。引起不均勻沉降的因素之一便是草炭土的高滲透特性，與其他黏性軟土相比，草炭土中的有機質含量高，土體內部布滿了植物纖維，因而孔隙發育，從外部表現出疏松多孔的形態。草炭土的分解度與有機質含量與其高滲透性緊密相關，毛文飛(2015)、汪之凡等(2017)通過室內試驗，對吉林省靖宇縣采取的草炭土樣的滲透性的影響因素以及影響規律進行了較為全面的分析，認為草炭土的分解度和有機質含量對草炭土的滲透性產生影響，且分解度的影響更為顯著，滲透率隨著分解度的增大而降低。

近年來，CT掃描技術在巖土體的微觀結構研究方面越來越廣泛。CT掃描技術是一種無損的成像技術，能夠探測材料內部的幾何形態，相比于其他試驗方法，CT技術適用于在不破壞巖土體孔隙結構的情況下獲取到巖土體的細觀結構圖像，進而通過三維重構技術實現對巖土體的三維形態的研究。諸多學者運用CT掃描技術來開展對巖土體結構的三維重構進而深入研究。張標等(2021)通過CT掃描技術并使用Avizo重構了珊瑚骨架灰巖的三維數字孔隙模型，毛偉澤等(2020)基于CT掃描技術實現了對花崗巖礦物顆粒三值化的三維重構。Mimics為比利時Materialise公司開發的面向醫學影像技術的專業3D建模軟件(蘇奎等， 2020)。作為專業的醫用影像建模軟件，Mimics支持多種無損壓縮的格式圖片的導入，且Mimics具有成熟的三維尺度上的體積計算功能，適用于CT掃描序列圖像上的信息識別以及重構三維模型(王嬌等， 2015)。

此外，隨著計算機技術的普及和發展，越來越多的學者基于CT掃描的圖像在土體滲流數值模擬上完成了許多卓有成效的研究，宋帥兵(2020)分別基于有限元數值理論模型、孔隙網絡模型計算得到高廟子膨潤土的滲透率，并加以滲透率試驗進行驗證，認為CT 微米級模型和FIB/SEM納米級模型分別具有最高和最低的滲透率數值。

Lattice Blotzmann方法是20世紀80年代發展而來的介于微觀分子動力學和宏觀連續介質模型之間的介觀模擬方法，在多孔介質滲流的模擬方面具有獨特的理論性。昆明理工大學申林方、王志良等人(申林方等， 2014，2015； 侯汝幾等， 2015； 崔冠哲等， 2016； 王志良等， 2018))采用D2Q9模型運用格子Boltzmann方法對二維平面的二值圖像進行數值模擬并研究了土體滲流場的變化機制。此外，國外的研究人員以LBM理論研發了PALABOS代碼庫(Latt et al.，2021)，能夠模擬解決各種不同環境下的流體流動問題(Huang et al.，2015)。在代碼庫的基礎上修改其多孔介質的流動模擬程序能夠實現對巖、土體內部孔隙的滲流模擬。

本文以吉林省敦化市的草炭土為研究對象，通過CT掃描試驗獲取草炭土的切片斷層圖像。通過Mimics21.0體積計算的孔隙率結合真實孔隙率對比確定草炭土CT圖像的二值化的最優閾值?；诓萏客恋闹貥嬋S模型以及序列CT圖像的橫截面與縱截面分析了草炭土的內部的孔隙形態。從LBM的基本理論出發，通過改編程序對草炭土CT序列圖像進行單相BGK滲流模擬，研究了LBM理論下草炭土的滲透率以及其內部流線與流速的分布?；谝陨线@些研究內容，本研究能夠作為草炭土或軟土滲流模型研究的理論基礎，可以指導這些地區的工程建設。

1 試驗方法

1.1 樣品的采集

草炭土樣品采集自吉林省敦化市江源鎮地區地表下的典型草炭土層。為了獲取不受擾動的草炭土試樣，試驗坑被開挖至地面1.5m左右的深度，用內徑為0.1m的PVC圓筒推切草炭土并將其包裹后取出(圖1)。將取出的樣品用黑色塑料袋與膠帶進行二次密封處理，并標記土樣編號，直接運送至實驗室進行冷藏。冷藏期間保證冷藏室溫度為3～5℃，防止樣品中的植物纖維取出后再度被微生物分解。

圖1 草炭土的現場樣品

1.2 土體的基本物理性質指標

本文選取地表下1.05～1.20m深度內的3份草炭土樣品，分別通過室內試驗測得其土體的天然密度ρ，土粒密度ρs、天然含水率ω等基本物理指標，經式(1)、式(2)換算得到3份草炭土樣的孔隙率n，如表1：

表1 草炭土的基本物理性質

(1)

(2)

1.3 CT掃描試驗

本文進行CT掃描試驗所用的儀器為比利時BRUKER-MICROCT公司生產的 Skyscan 1172 型桌面X射線CT系統，試驗在中國科學院長春應用化學研究所完成。CT掃描的草炭土樣品選擇的是1.1中進行室內試驗的同一批冷藏封存的3份樣品，試樣編號一一對應。試樣尺寸是直徑約8mm，高約10mm的圓柱體。將樣品放置在真空冷凍干燥機內，在-45℃的條件下進行抽真空24h，使其內部的水分全部被蒸發。通過X射線穿透樣品進行逐層掃描，便可以得到在不同角度方向該層面的“切片”投影圖像(鞏林賢， 2020)，并通過 Micro-CT掃描儀自帶的容積重構軟件Nrecon(V1.6.1.0)中的FDK重建算法處理“投影”圖像，將投影圖像轉換成二維斷層切片圖像。

掃描后得到的二維斷層切片圖像的像素單元尺寸為5 μm，其能反映的最小孔隙直徑為5μm。最小掃描層厚為7.5μm，經掃描及投影后得到草炭土3份土樣的斷層切片圖像數目分別為1332張， 1324張， 1323張。二維斷層圖像的分辨率為2000×2000，圖片格式為8位位圖的BMP格式，能夠無損且顏色準確地保留圖像信息。

2 草炭土CT序列圖像的二值化最優閾值選取及三維模型分析

2.1 CT圖像的預處理

在采集草炭土樣品的CT圖像的過程中，在草炭土的側邊加上了固定環圈防止掃描過程中草炭土自身發生擾動，因而在二維斷層切片圖像中草炭土的周邊顯現出固定環圈的痕跡。同時，由于設備運行過程中，因電壓不穩造成圖像局部灰度失衡，出現細小的斑點，使得圖像信息失真(趙玥等， 2018)。為了消除以上原因對圖像信息產生的干擾，本文采用MATLAB分別對原圖像進行裁剪； 利用中值濾波法進行圖像去噪處理，處理前后效果見圖2。

圖2 圖像的預處理

經過程序處理的圖像，像素大小由原圖像的2000×2000裁剪為1100×1100，消除了圖像中交織的噪聲，使得肉眼能夠更清晰地分辨出孔隙與土顆粒的信息。

2.2 CT圖像的二值化的最優閾值選取

2.2.1 二值化最優閾值選取的方法

Mimics中能夠通過閾值(Threshold value)劃分得到該灰度值區間下的蒙罩(Mask)，蒙罩部分能夠快速且準確地計算蒙罩體積(Mask volume)。通過調整閾值形成各個部分的灰度值區間，生成并查看mask部分的屬性(Properties)能夠得到各個部分體積計算的結果，從而確定各個部分的體積大小。通過計算V孔隙mask/(V土顆粒mask+V孔隙mask)能夠得到Mimics計算的CT序列圖像的孔隙率。

本文使用閾值分割的方法進行CT圖像的二值化。采用試湊法，結合室內試驗換算得到的草炭土樣品的真實孔隙率與Mimics的體積計算功能，通過調整土顆粒與孔隙的閾值形成各自的灰度值區間并生成蒙罩，分別計算不同閾值下的孔隙率。當樣品的真實孔隙率與Mimics計算得到的孔隙率一致時，以該閾值作為最優閾值，進行圖像的二值化，以保證后續圖像信息轉換的準確性。

2.2.2 二值最優閾值的確定

將預處理完畢后的草炭土樣品CT序列圖像導入至Mimics 21.0的工作區間上，在軸斷面(Axial)窗口中，使用某條直徑做剖面線(Profile Line)進行灰度值提取，見圖3。

圖3 草炭土CT圖像剖面線

不同的介質由于密度不同，灰度值存在差異，能夠通過灰度值的大小來進行定量識別。從圖4可見，當剖面線經過孔隙部分時，灰度值位于0～80之間，波谷處呈鋸齒狀波動； 當剖面線經過土顆粒部分時灰度值明顯上升，位于60～200之間，曲線波動幅度小。分析得出，由于孔隙的邊緣不明顯，在灰度值位于60～80之間時，像素點的位置既有可能是孔隙部分也有可能是土顆粒部分。此外，當剖面線經過土顆粒中夾帶的礦物質時，出現明顯的波峰，灰度位于約150以上，礦物質顆粒越大越明顯，灰度波峰越高。

圖4 剖面線土顆粒、孔隙、礦物質灰度值分布圖

本文分別計算了土顆粒與孔隙的閾值為60～80下的孔隙率，得到閾值大小與土體孔隙率關系曲線，并結合室內試驗換算得到的同一深度下的3份草炭土樣品的真實孔隙率，投射至關系曲線上，如圖5所示。

圖5 土體真實孔隙率對應的閾值

由圖5可知，本文將76作為土樣1與土樣2 CT圖像的最優閾值， 78作為土樣3的最優閾值。

2.2.3 二值化處理

在得到圖像的二值最優閾值后，采用MATLAB編程進行圖像的二值化，效果如圖6所示。

圖6 二值化效果

經過選取最優的閾值進行圖像分割后，草炭土的CT圖像分成了孔隙與土顆粒兩部分，分割的效果較好，能夠完整地提取圖像中的信息，為后續草炭土的滲流模擬研究提供精確信息。

2.3 草炭土CT序列圖像的三維模型與分析

在Mimics 21.0的工作區間上將閾值區間劃分為土顆粒，礦物質顆粒和孔隙3個部分。將每部分存放于各自的Mask內之后利用腔隙填充(Cavity fill)、編輯蒙罩(Edit masks)等方法修改蒙罩中的像素，編輯處理完成后經3D Calculate得到各個部分較為精確的三維結構模型，組合后便得到了完整的草炭土三維重構模型，如圖7所示。

圖7 草炭土三維重構模型

基于三維重構模型，能夠獲得土顆粒、孔隙、礦物質的空間分布規律。從圖上可以看出，草炭土中的土顆粒與孔隙之間交互分布，礦物質顆粒均勻地分布在土顆粒中?？紫渡y地分布于模型的各個區域，孔徑大小不一，大孔徑與細孔徑的孔隙形態各異，大孔徑孔隙的貫通性良好，而細孔徑的孔隙則散布于各個區域。CT序列圖像的橫截面與縱截面圖像上能夠清楚地看見大孔徑孔隙的形態，如圖8所示。草炭土的土顆粒部分為有機質并夾雜有植物纖維，植物纖維在有機質中呈條狀分布，從圖中也能看出從植物纖維夾有長條狀或圓管狀的植物孔隙。結合三維模型與截面圖像能夠發現，植物纖維在土中呈現架空狀的結構也會形成大孔徑孔隙，草炭土孔隙的形成受到植物纖維分布的影響，因此孔隙結構與植物纖維的含量和分解度緊密相關。

圖8 植物纖維形成的根管狀孔隙

3 草炭土CT序列圖像的滲流模擬

3.1 LBM的基本原理

Lattice Boltzmann Method(LBM)基于Boltzmann于1872年從分子運動論和統計物理的理論推導出來Boltzmann方程，從微觀動力學角度出發，將流體連續介質看作大量位于網格節點上的離散流體質點粒子(何雅玲等， 2009)。粒子按碰撞和遷移規則在網格上運動，通過對各網格流體質點及運動特征進行統計分析，獲得流體宏觀運動規律。LBM完整的格子模型由分布函數的演化方程、離散速度模型、平衡態分布函數3部分組成(楊佳慶， 2010)。

本文使用Qian et al.(1992)提出的D3Q19模型，為了估算碰撞項，這里運用BGK(Bhatnagar-Gross-Krook)動力學對碰撞項進行線性化的處理，BGK單相松弛LBM模型中松弛時間τ由單個參數定義，以此來實現對Navier-Stokes方程的求解(Mandzhieva， 2017)。模型的粒子分布函數如式(3)所示：

fi(x，t)-fi(x+eiδt，t+Δt)=

(3)

式中：τ為松弛時間(粒子碰撞速率)，與流體的黏度有關； δt為離散時間步長；fi(x，t)為平衡粒子分布函數。方程左邊表示的是粒子的流動，粒子在每個方向將離散速度的傳遞至各自相鄰的節點； 方程右邊表示的是粒子間的碰撞出現的平衡粒子分布的一部分松弛效應。

D3Q19的空間離散速度為圖9所示。

圖9 D3Q19模型

D3Q19模型中平衡態分布函數表示為：

(4)

式中：ρ為密度；cs為格子聲速；ωi為權系數；u為流體的宏觀平均流速；ei為粒子的離散速度。

D3Q19模型中格子聲速cs以及權系數ωi分別為：

(5)

(6)

粒子的離散速度ei的大小由格子速度決定，格子速度c=δx/δt，δx為格子步長。

BGK單相松弛LBM模型存在著的滲透率會隨松弛頻率的增加幾乎呈現線性增加的狀態，但是當c=1時能夠避免這一非物理現象(Mandzhieva， 2017)。因此本研究選擇δx=δt，使得c=δx/δt=1，從而實現在細觀尺度下由Lattice Blotzmann方程來替代Navier-Stokes方程。

因而，流體的宏觀流動特性可以由細觀上的每個時間步的分布函數近似得到，如式(7)，式(8)所示：

(7)

(8)

式中：ρ為模型的宏觀密度；u為模型的宏觀速度。

3.2 模型邊界條件的處理

在滲流問題的模型的設置上，邊界條件一直都對數值模擬的結果起著至關重要的作用(周瀟， 2016)。結合本研究問題，采用周期性邊界與反彈邊界進行數值模型的邊界條件處理。設定草炭土模型的流體出入口采用周期性邊界格式，不透水邊界以及草炭土顆粒采用反彈格式。

3.2.1 反彈格式邊界

反彈格式假設粒子與土顆粒表面以及壁面發生碰撞，碰撞后粒子速度發生逆轉，速度大小不變，朝著相反的方向發生反彈，反彈格式邊界處粒子的運移方式如圖10所示。

圖10 反彈格式邊界處流體粒子運移的4個階段

流體粒子遷移到邊界處，速度發生逆轉，碰撞后粒子沿原路返回，在圖10中發生的碰撞分布函數為：

f2， 8， 9(i，l)=f4， 7， 10(i，l)

(9)

3.2.2 周期性邊界

設置上下邊界為周期性邊界條件，假設計算域的出口的流體再次從入口流入，保證在非穩態與周期性的出入口邊界條件下使得滲流模擬的過程趨于穩定，最終的計算結果得到收斂，計算得到草炭土的滲透率。

3.3 程序的編制

3.3.1 基于MATLAB的圖像信息轉換程序

通過MATLAB中的imread函數功能，自編程序讀取CT序列圖像的信息，提取圖像數據并轉化成數據結構。數據結構包含序列圖像各個像素點的索引值，儲存在代碼能夠讀取的二進制(.dat)文件中。程序設定的數據結構表示：

B(i，j，numFiles)=

(10)

式中：i，j表示的是圖像中像素單元的位置；numFiles表示的是圖像的序列號。

3.3.2 基于LBM的滲流模擬程序

滲流模擬程序基于PALABOS開源代碼庫(https：∥palabos.unige.ch/)，采用C++語言依照以上基本理論編制了草炭土滲流的計算程序。根據構建的數據結構表征的數值，在程序中規定：流體粒子在讀取到值為0的體素中正常遷移，值為1的體素速度發生反彈，值為2的體素中被設定為無法遷移。出入口兩端設置固定的壓差，模擬流體在孔隙內線性流動，當相鄰時間步長的速度變化量小于收斂值，則判定流體的運移達到穩定狀態，停止模擬。

3.3.3 程序中滲透率的計算與單位轉換

本研究中模擬是水這一不可壓縮流體在孔隙內的流動。程序基于達西定律由流體粒子的流速計算求解土體的滲透率，關系式如下：

(11)

在格子單位領域滲透率是以無量綱單位輸出的，使用式(12)將其轉換為物理單位。

(12)

式中：Lphysical等于圖像分辨率與像素單元大小的乘積；LPALABOS等于格子單位的長度。

3.4 草炭土CT序列圖像的細觀滲流模擬

3.4.1 CT序列圖像的選取

由于LBM模型無法完整地讀取圓形切片圖像上的信息，因此在滲流模擬的研究中重新從原始圖像中裁剪出方形區域，并對圖像進行預處理、二值化。選取處理后的草炭土CT切片圖像中連續的814張圖像，將序列圖像轉化為二進制的信息文件。草炭土整體的序列圖像如圖11所示：

圖11 草炭土方形CT序列圖像示意圖

3.4.2 模擬初始條件的設定

模型模擬設定的初始條件見表2。

表2 模擬的初始條件

3.4.3 模擬結果分析

3.4.3.1 數值模擬的草炭土的滲透率結果

基于上述模擬過程分別測得3份草炭土樣品的滲透率，結果如表3所示。

表3 格子單位與物理單位下計算的滲透率

本文得到基于LBM理論下3份土樣的模擬滲透率值。通過與文獻對比發現，基于LBM法計算得到的草炭土的滲透率與Avizo軟件的X-lab算法(鞏林賢， 2020)以及實驗室實測(汪之凡等， 2017)的滲透率值接近，說明基于LBM理論的模擬結果具有一定的可靠性。

3.4.3.2 流速云圖與流線圖的分析

本文結合Paraview可視化軟件，將模擬得到土樣1的VTI格式文件導入至Paraview中進行可視化。得到圖12草炭土滲流路徑云圖，圖13為不同數量位移積點的流速流線分布圖。從流速云圖上我們可以看出，草炭土的孔隙分布復雜，孔隙孔徑大小不一，大孔徑的貫通性優于散布的細孔徑孔隙，存在著大孔徑滲流通道。在土體滲流穩定后，流體流動的方向優先選擇連通性良好的孔隙所形成的通道。

圖12 土樣1孔隙內部的滲流路徑云圖(上方入，下方出)

圖13 土樣1孔隙內部的流線分布圖

草炭土的孔隙內部的流線分布圖中，Number of Points是指位移積點的數目。結合流速云圖我們可以看出，從Number of Points=1000～10000的過程中，位移積點的位置越來越集中于貫通的大孔徑孔隙。大孔徑孔隙內的流線數目越來越多，而部分連通的孔隙散布著細小的流線。說明大孔徑內的流線數目多于其余不貫通孔隙，進一步表明了大孔徑孔隙是草炭土滲流的主要通道，對草炭土的滲透性起著至關重要的作用。

3.4.3.3 流速的變化分析

本文分別從流速云圖XZ截面(Y=314處)上選取Z=381與Z=240兩條X向(豎向)截面截線(圖14)以及XY截面(Z=237處)上選取Y=229與Y=414兩條X向(豎向)截面截線(圖15)，提取截面截線上的流速信息并繪制流速分布曲線。

從流速云圖(圖14a、圖15a)中可以看出，大孔徑孔隙內流體滲流量大，流速明顯高于不貫通的細孔徑孔隙。結合流速分布曲線(圖14b、圖15b)可知大孔徑孔隙通道出現了整個流域內最大的滲流速度，這些大孔徑孔隙之間互相貫通，流體在其中的流通性非常好； 而部分連通的孔隙不存在較大的流速，流體在其中流動斷斷續續，這也是由于草炭土土體孔隙復雜，水在不貫通孔隙流動的過程中能量被消耗，流速逐漸減緩。

圖14 土樣1 XZ縱截面上的代表性截線流速

圖15 土樣1 XY縱截面上的代表性截線流速

由此可見土體的孔隙結構尤其是大孔徑孔隙的數量在一定程度上影響著草炭土的滲流特性。而大孔徑孔隙的分布與土中植物纖維的分布有關，因而草炭土的滲透率受到植物纖維的含量以及分解度影響，植物纖維的含量越多，越分解不完全，草炭土內部越架空，內部的孔隙貫通性越好，滲透率越高； 反之，則滲透率越低。

4 結 論

(1)本文采用試湊法，通過Mimics21.0的Mask體積計算功能得到草炭土CT序列圖像的孔隙率，結合土體真實孔隙率，得到了3份樣品圖像二值化處理的最優閾值為76， 76， 78。該方法確定的CT圖像二值最優閾值能夠精確地提供圖像提取的信息。

(2)基于三維重構模型并結合CT序列圖像的橫截面與縱截面發現草炭土的孔隙結構復雜，孔徑的大小不一，孔隙形態各異，植物纖維形成的圓管狀孔隙通道以及架空狀孔隙結構是大孔徑孔隙形成的主要原因。因此土體內部大孔徑孔隙的分布與植物纖維的分布緊密相關，植物纖維形成的大孔徑孔隙會成為草炭土滲流的優先路徑。

(3)基于LBM理論，采用C++算法改編程序，在細觀尺度上對草炭土的孔隙進行單相BGK滲流模擬，計算了草炭土的格子單位下的滲透率，并轉化為物理單位下的滲透率，計算得到3份土樣的滲透率分別為7.1358 μm2、5.5194μm2、6.713675μm2。通過與以往的文獻數據對比發現LBM算法的結果與Avizo的X-lab算法以及實驗室實測值的滲透率值接近，說明基于LBM理論的模擬結果具有一定的可靠性。

(4)通過Paraview可視化的滲流路徑云圖、流線分布圖、以及代表性截線流速曲線分析發現，在細觀尺度上草炭土內的流體優先選擇大孔徑的滲流通道，流線大都集中于這些大孔徑的孔隙中，而部分連通的孔隙散布著細小的流線。大孔徑孔隙通道內出現了整個流域內最大的滲流速度。細孔徑孔隙內不存在較大的流速，流體在其中流動斷斷續續，能量消耗快于大孔徑孔隙，流速逐漸減緩。說明草炭土的孔隙結構尤其是大孔徑孔隙的數量會影響土體的滲透性。而孔隙結構與植物纖維的含量以及分解程度有關，因而土體中的植物纖維的含量以及分解度會影響土體的滲透性，植物纖維的含量越多，分解度越低，草炭土的滲透率越高，反之則越低。