基于三維X-CT圖像的結(jié)皮土壤孔隙結(jié)構(gòu)特征與滲透率

許智隼，胡五龍,2※

（1.武漢理工大學(xué)理學(xué)院，武漢430070；2.武漢理工大學(xué)綠色智能江海直達(dá)船舶與郵輪游艇研究中心，武漢430070）

0 引 言

在干旱半干旱地區(qū)，表層土壤會(huì)在外部物理夯實(shí)或內(nèi)部微生物作用下形成一種低入滲率、高抗剪強(qiáng)度的土壤結(jié)皮[1]，即土壤物理結(jié)皮或土壤生物結(jié)皮。根據(jù)不同的形成機(jī)理和發(fā)育過(guò)程，土壤物理結(jié)皮可以分為由雨水等外力打擊而形成的結(jié)構(gòu)性結(jié)皮和由徑流沖刷引起的沉積性結(jié)皮[2-3]。土壤結(jié)皮中的微生物可以有效地富集并吸收土壤中的重金屬，對(duì)治理環(huán)境污染具有一定的積極意義[4]。同時(shí)，土壤結(jié)皮也會(huì)顯著降低土壤水分向地下入滲，增加地表徑流[5-6]，造成水土流失[7]，進(jìn)而可能引發(fā)一系列嚴(yán)重的生態(tài)環(huán)境災(zāi)害。目前國(guó)內(nèi)外對(duì)土壤結(jié)皮的研究主要聚焦于其成因和功能方面，而對(duì)土壤結(jié)皮的微觀孔隙結(jié)構(gòu)及其對(duì)水力特性的影響研究較少[8-9]。

土壤孔隙是氣體、水和營(yíng)養(yǎng)物質(zhì)輸運(yùn)的通道，同時(shí)也是土壤微生物的棲息地[10-12]。因此，土壤孔隙結(jié)構(gòu)在很大程度上影響或決定了土壤肥力，對(duì)生態(tài)環(huán)境、地下水系統(tǒng)及大氣循環(huán)有著重要影響。同時(shí)，土壤中的所有物理過(guò)程及生化反應(yīng)都發(fā)生在孔隙中，而目前的試驗(yàn)技術(shù)手段很難對(duì)孔隙中的這些過(guò)程和反應(yīng)進(jìn)行準(zhǔn)確觀測(cè)，因此基于真實(shí)土樣三維孔隙結(jié)構(gòu)從孔隙尺度模擬土壤中水的輸運(yùn)，有助于精確描述土壤中水的輸運(yùn)過(guò)程和相關(guān)機(jī)制。

土壤孔隙結(jié)構(gòu)的傳統(tǒng)研究方法有切片法[13]、壓汞法[14]和氮吸附法[15]等，這些方法不僅會(huì)破壞土壤原有的孔隙結(jié)構(gòu)，而且難以重構(gòu)土壤內(nèi)部的三維孔隙結(jié)構(gòu)。顯微成像技術(shù)，尤其是 X射線斷層掃描成像技術(shù)（X-ray Computed Tomography，X-CT）的快速發(fā)展，使人們能夠獲取真實(shí)土壤的三維孔隙結(jié)構(gòu)并對(duì)其進(jìn)行分析。與傳統(tǒng)的研究方法相比，X-CT是一種不損害土壤孔隙結(jié)構(gòu)的三維成像技術(shù)，其分辨像素尺寸甚至可以小于1μm，已被廣泛應(yīng)用于油氣開(kāi)采、地下水鹽運(yùn)移、污染物擴(kuò)散等與土壤相關(guān)的工程實(shí)踐與科學(xué)研究中[16]。Katuwal等[17]和Hu等[18]利用X-CT研究了土壤孔隙的幾何結(jié)構(gòu)與其水氣特性間的關(guān)系。程亞南等[19]利用CT掃描技術(shù)重構(gòu)了土壤的三維孔隙結(jié)構(gòu)，并建立了孔隙網(wǎng)絡(luò)模型來(lái)預(yù)測(cè)孔隙結(jié)構(gòu)的水力學(xué)性質(zhì)。楊永輝等[20]利用X-CT研究了不同土壤結(jié)構(gòu)改良措施下的土壤孔隙特征，為合理地應(yīng)用耕種方式提供了科學(xué)依據(jù)。利用 X-CT技術(shù)對(duì)土壤結(jié)皮進(jìn)行研究，一方面可以在真實(shí)結(jié)皮土樣三維孔隙結(jié)構(gòu)基礎(chǔ)上分析其孔隙特征及分布規(guī)律，另一方面可以基于真實(shí)土樣結(jié)構(gòu)圖像孔隙尺度模擬結(jié)皮土樣中水的滲透過(guò)程。然而，迄今為止，尚未有利用X-CT研究土壤結(jié)皮三維孔隙結(jié)構(gòu)及其分層特征的報(bào)道，對(duì)土壤結(jié)皮的結(jié)構(gòu)特征和滲流特點(diǎn)的孔隙尺度認(rèn)識(shí)也尚不明確。格子玻爾茲曼方法（Lattice Boltzmann Method，LBM）是近三十年發(fā)展起來(lái)的流體數(shù)值計(jì)算方法，它采用粒子分布函數(shù)直接模擬計(jì)算域中流體粒子的運(yùn)動(dòng)，具有成像清晰、算法簡(jiǎn)單、易于并行和適應(yīng)復(fù)雜幾何邊界條件的特點(diǎn)，因而非常適合用來(lái)模擬多孔介質(zhì)中傳熱和傳質(zhì)過(guò)程[21]。結(jié)合X-CT技術(shù)，利用LBM可以更逼真地模擬孔隙介質(zhì)中流體運(yùn)動(dòng)[22]，為探究多孔介質(zhì)的滲流機(jī)理提供了一種高效精確的方法[23]。

本文擬對(duì)3個(gè)有表層結(jié)皮的土樣孔隙結(jié)構(gòu)及滲透率進(jìn)行分層分析，先采用X-CT和圖像處理技術(shù)獲取三個(gè)土樣的三維孔隙結(jié)構(gòu)圖像，然后將這些圖像沿深度方向分層，分析各深度層內(nèi)土樣的孔隙結(jié)構(gòu)特征。最后結(jié)合格子玻爾茲曼方法計(jì)算各深度層的滲透率，探討土樣的孔隙特征對(duì)滲透率的影響，以期為干旱區(qū)半干旱區(qū)合理管理利用地下水資源和保護(hù)生態(tài)環(huán)境提供參考。

1 材料與方法

1.1 三維X-CT圖像獲取

本文對(duì) 3個(gè)表層有結(jié)皮的土樣進(jìn)行孔隙結(jié)構(gòu)和滲透率分析。原土樣采自英國(guó)貝德福德郡一處田地，提取地表10 cm厚土層，在實(shí)驗(yàn)室烘干搗碎去除植物根莖和粗顆粒后做成半徑和高均為10 cm的土樣，三個(gè)土樣均為粉黏土。通過(guò)噴淋方式模擬降雨，給定三個(gè)土樣不同的降雨量。然后將土樣至于置于實(shí)驗(yàn)室環(huán)境試驗(yàn)箱中模擬半干旱地區(qū)蒸發(fā)環(huán)境條件，待表層出現(xiàn)明顯結(jié)皮后，用環(huán)刀從土樣上部重新提取小柱狀土樣，將提取的小柱狀土樣通過(guò)Phoenix Nanotom掃描系統(tǒng)獲得一系列的X射線投影圖像，然后使用掃描系統(tǒng)的配套軟件 Phoenix datos|x2和重建算法將投影圖片重構(gòu)為三維X-CT圖像，并導(dǎo)出三維圖像的二維連續(xù)切片圖。X-CT圖像的顯示分辨率與掃描視野相關(guān)，本文為了顯示土樣的全貌，采用了較低顯示分辨率的掃描方式，切片圖像的分辨率為40μm。掃描的柱狀土樣周邊孔隙結(jié)構(gòu)會(huì)受取樣時(shí)環(huán)刀影響發(fā)生改變，為此截取三個(gè)土樣中心區(qū)域的掃描圖像用于后文的分析和計(jì)算。由于原始的切片圖像是灰度圖，本文僅關(guān)注土樣的孔隙特征，因此利用開(kāi)源軟件 ImageJ將這些連續(xù)切片圖像進(jìn)行二值化處理，方便進(jìn)一步的可視化和計(jì)算[24]。最終的三維圖像如圖1所示，三個(gè)土樣的 尺 寸 分 別 為 937×937×1000 、 912×912×1000 和933×933×1000體素（voxel），其中淺色部分為固相，深色部分為孔隙。為了討論土樣孔隙結(jié)構(gòu)沿深度方向（圖1中z方向）的變化規(guī)律，由上表面開(kāi)始每100 voxel（即4 mm）取一土層，并依次編號(hào)為土層1～9。

1.2 結(jié)構(gòu)特征分析方法

在經(jīng)過(guò)二值化處理后，土樣的結(jié)構(gòu)信息被存儲(chǔ)在一個(gè)三維數(shù)組中，數(shù)組中的每個(gè)元素都為0或1，分別表示土樣的孔隙相和固相，因此土壤孔隙率可由數(shù)組中“0”元素的數(shù)量簡(jiǎn)單得到。

式中φ表示孔隙率；N0表示數(shù)組中“0”元素的數(shù)量；N0&1表示數(shù)組中“0”元素和“1”元素的總數(shù)。

孔隙變異度是指當(dāng)前土層與其相鄰?fù)翆拥目紫堵什町惓潭萚12]，可用下式來(lái)計(jì)算：

式中E表示孔隙變異度，%；φu和φb分別表示當(dāng)前土層和下層土的孔隙率。

在計(jì)算孔隙直徑時(shí)，采用形態(tài)模型方法[25]，即以不同直徑的球填充孔隙空間，對(duì)于每個(gè)孔隙體素，取包含該孔隙體素的最大球直徑為其對(duì)應(yīng)的孔隙直徑。由于三維圖像是由一個(gè)個(gè)體素組成的三維矩陣，而體素為一個(gè)邊長(zhǎng)為 1的小立方體，因此孔隙之間可以通過(guò)立方體的角點(diǎn)、棱線和面來(lái)連接。由于不連通孔隙和微小孔隙（孔隙直徑小于分辨率）中的水分幾乎不發(fā)生遷移，故在本文研究中忽略這部分孔隙體積，并假設(shè)流體僅在連通孔隙中傳輸，僅當(dāng)兩個(gè)孔隙體素以面連接才視為連通，孔隙連通率可由下式給出：

式中Rp為孔隙連通率，%；φ和φc分別為總孔隙率和連通孔隙率。

土壤孔隙復(fù)雜度可由分形維數(shù)[13]來(lái)進(jìn)行定量評(píng)價(jià)。在二維情況下，規(guī)則圖形的周長(zhǎng)-面積滿(mǎn)足以下關(guān)系：

式中C為二維圖形的周長(zhǎng)，mm；A為二維圖形的面積，mm2。

對(duì)于不規(guī)則圖形，可以使用分形周長(zhǎng)[26]來(lái)代替上式中的實(shí)際周長(zhǎng)。分形周長(zhǎng)和分形面積之間滿(mǎn)足

式中ε表示測(cè)量尺度，C(ε)和A(ε)表示在該尺度下測(cè)得的周長(zhǎng)和面積。因此，分形維數(shù)D可以通過(guò)采集不同測(cè)量尺度ε的C和A，然后繪制lgC-lgA曲線來(lái)獲得。由于實(shí)際中往往只應(yīng)用了一個(gè)尺度來(lái)測(cè)量圖形的相關(guān)信息，缺少其他尺度下的C和A值，因此常用下式來(lái)近似地求取D：

上述求取分形維數(shù)D的方法稱(chēng)為“小島法”[12,27]，其討論的是復(fù)雜二維圖形的周長(zhǎng)-面積關(guān)系，本文研究土壤三維孔隙結(jié)構(gòu)，需在上式基礎(chǔ)上進(jìn)行改進(jìn)。假設(shè)復(fù)雜三維圖形的表面積S（mm2）和體積V（mm3）滿(mǎn)足

那么，三維圖形的分形維數(shù)D可以近似地表示為

式中分形維數(shù)D表征孔隙結(jié)構(gòu)復(fù)雜度（后文簡(jiǎn)稱(chēng)“孔隙復(fù)雜度”）；S表示三維圖形的分形表面積，而V表示其分形體積，二者的單位需要相互對(duì)應(yīng)，在本文中分別為mm2和 mm3。

從式（8）來(lái)看，孔隙復(fù)雜度D與比表面積具有相似的描述能力，即相同表面積下，體積越大的孔隙，D值越小，其孔隙結(jié)構(gòu)越簡(jiǎn)單。

1.3 滲透率計(jì)算

基于 Bhatnagar-Gross-Krook（BGK）算子的格子Boltzmann基本方程[28]如下式所示：

式中i表示離散速度方向；ei為無(wú)量綱離散速度；fi(x,t)和fieq(x,t)分別為t時(shí)刻x處的分布函數(shù)和平衡態(tài)分布函數(shù)；c=Δx/Δt表示離散格子速度，Δx和Δt分別為格子空間步長(zhǎng)和格子時(shí)間步長(zhǎng)，在 LB模擬中一般取為 1；τ為無(wú)量綱松弛時(shí)間。

平衡態(tài)分布函數(shù)fieq(r,t)由Maxwell-Boltzmann分布的Hermite展開(kāi)[29]來(lái)確定

式中cs=c/代表無(wú)量綱格子聲速；ρ和u分別表示x處流體粒子的宏觀密度ρ（g/cm3）和宏觀速度u（mm/s）；wi是權(quán)重系數(shù)，由離散速度模型所確定。

本文采用 D3Q19模型[30]，即將式（9）和式（10）在三維速度空間按照 19個(gè)方向進(jìn)行分解，格子模型如圖2所示。

那么，i方向上的離散速度矢量ei即圖2中的各空間矢量，如e0=（0, 0, 0），e1=（1, 0, 0），e7=（1, 1, 0），以此類(lèi)推。與之相對(duì)應(yīng)的權(quán)重系數(shù)分別取w0=1/3，w1-6=1/18，w7-18=1/36。

流體的宏觀密度和宏觀速度別由分布函數(shù)的一階矩和二階矩給出

式（11）和（12）得到的流體密度和速度為格子單位，可換算為實(shí)際物理單位。

流體的運(yùn)動(dòng)黏度υ（mm2/s）為

作用在流體粒子上的壓力P0（Pa）由狀態(tài)方程[21]給出

在模擬過(guò)程中，假設(shè)水是不可壓縮的，其密度ρ=1.0 g/cm3。在LB M模擬中，流體的壓力與密度成正比，在水密度恒定的條件下流場(chǎng)中的壓力處處相等，與流體力學(xué)中由壓力梯度來(lái)驅(qū)動(dòng)流體流動(dòng)的客觀事實(shí)相悖。為模擬土壤中水滲透過(guò)程，本文將土樣頂部的流體密度設(shè)置為ρ1=1.0+0.000 2Lzg/cm3，Lz是圖像的高度，土樣的底部則設(shè)置為ρ2=1.0 g/cm3，四周為周期性邊界條件。壓力邊界條件的具體格式參考Zou等[31]的非平衡態(tài)外推邊界。

當(dāng)土樣中的體積流量（m3/s）在前后兩個(gè)時(shí)間步間相差小于10-6，則認(rèn)為滲流達(dá)到穩(wěn)定狀態(tài)。此時(shí)，可由達(dá)西定律求出滲透率[32]

式中q為體積平均流速，mm/s；k為滲透率，mm2；?P為壓力梯度，Pa/mm。

在z方向（豎直方向）上施加壓力梯度時(shí)，z方向上的滲透率分量kzz表示為[23]

式中的平均流速qz（mm/s）由下式給出：

式中∑ui,z表示圖像中各個(gè)流體粒子在z方向上的速度分量之和，mm/s；N=NxNyNz，為土樣圖像的體素（voxel）個(gè)數(shù)，Nx、Ny和Nz分別為土樣圖像x、y、z方向的體素?cái)?shù)。

模擬過(guò)程中，所有物理量包括長(zhǎng)度、密度等均直接采用格子單位，計(jì)算結(jié)果最后轉(zhuǎn)換成實(shí)際的物理單位，1 lu =1 voxel（lu為格子單位）。

2 結(jié)果與分析

2.1 孔隙結(jié)構(gòu)特征分析

2.1.1 孔隙率與孔隙連通性分析

表1顯示了3個(gè)結(jié)皮土樣的基本信息，從統(tǒng)計(jì)結(jié)果來(lái)看，盡管三個(gè)土樣的孔隙率具有較大的差異，但其平均孔隙直徑均在0.2 mm左右。土樣1具有較好的整體連通性，孔隙連通率（連通孔隙體積與總孔隙體積之比）達(dá)91.26%，而土樣2和3的整體連通性則較差。

表1 土樣的基本幾何參數(shù)Table 1 Basic geometric parameters of soil crusts

圖3a為三個(gè)土樣各水平切片上的孔隙率沿深度變化規(guī)律，從圖中可以看出，土樣表層具有較大的孔隙率，但與其緊鄰的下層土的孔隙率陡降，在1～2 mm深度達(dá)到一個(gè)極小值。從圖3b的孔隙變異情況來(lái)看，土樣上表面幾層切片在孔隙率上具有較大的差異，其變異系數(shù)E＞5%。這種差異可能會(huì)導(dǎo)致其滲透能力表現(xiàn)出截然不同的特點(diǎn)。

與圖3a相對(duì)應(yīng)，各土樣的1～9土層的最大深度分別為h=4～36 mm（每層4 mm）。圖4a顯示了三個(gè)土樣每個(gè)土層的孔隙率，從圖中可以看出，當(dāng)土層厚度取為4 mm時(shí)，三個(gè)土樣的連通孔隙率均表現(xiàn)出兩端孔隙率小、中間孔隙率大的整體變化趨勢(shì)。圖4a中每個(gè)土樣虛實(shí)線的相對(duì)位置關(guān)系顯示了土層孔隙連通性的差異。從計(jì)算結(jié)果來(lái)看，三個(gè)土樣各土層的孔隙連通率最大值分別為99.21%、94.64%和57.35%，處于土樣中部，而最小值分別為45.99%、27.30%，和11.74%，處于土樣的兩端，各土樣均表現(xiàn)出兩端連通性較差而中間連通性較好的整體趨勢(shì)。圖4b顯示了三個(gè)土樣各土層總孔隙率與孔隙連通率間的關(guān)系，對(duì)二者進(jìn)行相關(guān)性分析，顯示相關(guān)系數(shù)r=0.97，P＜0.01，表明總孔隙率與孔隙連通率存在極顯著的相關(guān)性。從圖中的趨勢(shì)來(lái)看，土壤孔隙率越高，其孔隙連通性越好。

2.1.2 孔徑分布分析

三個(gè)土樣的孔徑d分布如圖5所示，均近似服從泊松分布。盡管三個(gè)土樣孔徑分布區(qū)間各不同，但直徑d小于0.4 mm的孔隙均占據(jù)主要部分，分別占比94.02%，93.54%和 94.61%。而孔徑分布云圖顯示，直徑大于0.4 mm的孔隙主要分布在土樣的表層部分。

圖6顯示了各個(gè)土層的平均孔隙直徑，從圖中可以看出表層土（h≤4 mm）的平均孔徑與相鄰?fù)翆颖憩F(xiàn)出較大的差異，其余各層間的變化則較為平緩，總體上各土層的平均孔徑隨著深度的增加而減小，三個(gè)土樣各土層的平均孔徑分別由表層土（h≤4 mm）的0.43、0.37和0.50 mm降低至深層土（h≥32 mm）的 0.15、0.14和0.14 mm。這種變化規(guī)律一方面是由于表層土直接裸露在空氣中，與外界的物質(zhì)交換頻繁，且受雨水侵蝕、動(dòng)植物活動(dòng)的影響較大；另一方面，土壤水分?jǐn)y沙下滲時(shí)，細(xì)土粒在下部淤積，從而阻塞孔隙，使得下土層的孔徑變小。

2.1.3 孔隙復(fù)雜度分析

由于土壤中的水主要在連通孔隙中流動(dòng)，因此忽略孔徑小于圖像分辨率的微孔，在此僅討論土樣中連通孔隙的復(fù)雜程度。表2列出了三個(gè)土樣各土層中連通孔隙的表面積S與體積V值，在此基礎(chǔ)上根據(jù)式（7）計(jì)算各土層的孔隙復(fù)雜度D。

表2 各土層連通孔隙的表面積（S）與體積（V）Table 2 Surface area(S) and volume(V) of connected pores of each soil layer

最終的孔隙復(fù)雜度計(jì)算結(jié)果如圖7所示。從圖7 a可以看出，三個(gè)土樣間的孔隙復(fù)雜度有著明顯的不同，土樣1的D值最小，而土樣3的最大。從式（8）及孔隙復(fù)雜度的定義來(lái)看，孔隙復(fù)雜度與土壤孔隙的體積密切相關(guān)。圖7 b顯示了孔隙復(fù)雜度D與連通孔隙率間的關(guān)系，從總體上的趨勢(shì)來(lái)看孔隙率越高的土層，具有更小的D值，其孔隙結(jié)構(gòu)更加簡(jiǎn)單。對(duì)D和φ進(jìn)行Spearman相關(guān)性分析，二者的相關(guān)系數(shù)r=-0.90，P＜0.01，這表明孔隙復(fù)雜度與連通孔隙率呈極顯著負(fù)相關(guān)關(guān)系。圖7 b中，土層1-1和3-1偏離總體趨勢(shì)，這是由于土樣1和土樣2表層土壤因外部作用產(chǎn)生了大面積的凹陷，形成了不平整的表面，這些表面凹陷在算法中被當(dāng)成了孔隙統(tǒng)計(jì)進(jìn)來(lái)，從而使得表層土壤的孔隙復(fù)雜度大幅降低。

2.2 滲透率分析

2.2.1 孔隙結(jié)構(gòu)參數(shù)對(duì)滲透率的影響

通過(guò)Boltzmann方法和1.3節(jié)的邊界處理方式模擬了三個(gè)土樣各土層的滲流過(guò)程，并計(jì)算其滲透率，當(dāng)前后兩個(gè)時(shí)間步的滲透率值相差小于10-6，則認(rèn)為迭代收斂，滲流達(dá)到穩(wěn)定狀態(tài)[23]。最終的滲透率計(jì)算結(jié)果如圖8所示，從其沿深度變化的規(guī)律來(lái)看，各土層的局部滲透率沿深度方向先增加后減小，這與其連通孔隙率沿深度變化的規(guī)律基本一致。土層3-1滲透率kzz高達(dá)1093μm2，而土樣3其他土層的最大局部滲透率僅為23.522μm2，二者之間存在數(shù)量級(jí)的差異。分析土層3-1的孔隙發(fā)現(xiàn)在該土層存在一個(gè)上下通透的穿孔，如圖9所示，圖中白色透明的部分即為該土樣穿孔，該穿孔的存在是其滲透能力顯著強(qiáng)于其他土層的主要原因。

去掉土層3-1后，發(fā)現(xiàn)三份土樣中含結(jié)皮的表層具有較小的滲透率值，分別為7.472、0.006和6.960μm2，這表明結(jié)皮對(duì)土壤水的滲流有較強(qiáng)的阻滯作用。圖10顯示了土樣各結(jié)構(gòu)參數(shù)與滲透率間的關(guān)系，從圖中結(jié)果來(lái)看，滲透率隨孔隙率和孔隙連通性的提高而增大，隨孔隙復(fù)雜度D的提高而降低。滲透率與平均孔隙直徑間則沒(méi)有明顯的單調(diào)變化趨勢(shì)，從圖10 c來(lái)看，平均孔徑在0.16～0.20 mm間的土層具有較高的滲透率，這是由于這些土層均具有較高的孔隙連通率。

對(duì)這些結(jié)構(gòu)參數(shù)與滲透率進(jìn)行 Spearman相關(guān)性分析，結(jié)果如圖10所示。連通孔隙率與滲透率間極顯著正相關(guān)（相關(guān)系數(shù)為正，且P＜0.01），孔隙復(fù)雜度與滲透率極顯著負(fù)相關(guān)（相關(guān)系數(shù)為負(fù)，且P＜0.01），而平均孔徑和比表面積對(duì)滲透率大小的描述能力則較差（P＞0.05）。

應(yīng)用SPSS對(duì)圖10中的散點(diǎn)進(jìn)行曲線估計(jì)以選擇最優(yōu)的回歸模型，其結(jié)果如表3所示。在曲線估計(jì)時(shí)，選擇了SPSS中所有可選的回歸模型。從估計(jì)的結(jié)果來(lái)看，孔隙連通率φe和孔隙復(fù)雜度D能夠一定程度上對(duì)滲透率的數(shù)值結(jié)果做出預(yù)測(cè)，其決定系數(shù)較高且回歸分析具有顯著的統(tǒng)計(jì)學(xué)意義（P＜0.05），而應(yīng)用比表面積Ss和平均孔徑dave的回歸模型的誤差則較大。

表3 結(jié)構(gòu)參數(shù)與滲透率間的最優(yōu)回歸模型估計(jì)Table 3 Optimal regression model estimation between structural parameters and permeability

2.2.2 滲透率預(yù)測(cè)模型與改進(jìn)

從前文的結(jié)構(gòu)特征分析來(lái)看，孔隙率更高的土層具有更高的孔隙連通率和更小的孔隙復(fù)雜度，其相關(guān)性分析的P值均小于0.05。這些相關(guān)性分析表明，土壤孔隙率能夠在很大程度上描述土壤孔隙的結(jié)構(gòu)特征，并能通過(guò)孔隙率來(lái)預(yù)測(cè)土壤的滲透率。在已有的研究中，大多數(shù)的滲透率預(yù)測(cè)模型都是將孔隙率作為輸入?yún)?shù)。早在20世紀(jì)，Kozeny就將滲透率與孔隙率和孔隙形狀聯(lián)系起來(lái)，建立了水力傳導(dǎo)理論的基礎(chǔ)[33]，而后Carman在進(jìn)一步的工作中提出了廣為人知的Kozeny-Carman（K-C）方程[34-35]。然而原始K-C方程中有許多難以獲取的結(jié)構(gòu)參數(shù)，因此一些基于K-C方程的簡(jiǎn)化模型在實(shí)際應(yīng)用中更為廣泛[36]。利用兩種常用滲透率預(yù)測(cè)模型擬合LBM數(shù)值計(jì)算結(jié)果，如圖11所示，其中簡(jiǎn)化K-C方程的擬合參數(shù)C=（1.302±0.036）×104μm2，擬合程度R2=0.97；而Tomadakis的預(yù)測(cè)模型中的擬合參數(shù)C=（4.603±0.176）×104μm2，相應(yīng)的擬合程度R2=0.94。從擬合結(jié)果來(lái)看，簡(jiǎn)化K-C方程具有更高的精度，即k-φ之間滿(mǎn)足

應(yīng)用比表面積SS的原始K-C方程如下所示：

為使預(yù)測(cè)模型能夠較為準(zhǔn)確地描述孔隙復(fù)雜度D對(duì)滲透率的影響作用，并提升模型的預(yù)測(cè)精度，這里將D引入簡(jiǎn)化K-C方程。由于D與k負(fù)相關(guān)，即k隨D的增加而減少。而從D的定義來(lái)看，D與比表面積SS具有相似的描述能力，且D與k的相關(guān)程度要顯著高于D與SS的相關(guān)程度。因此，可用D來(lái)替換式（20）中的比表面積SS，即

改進(jìn)模型擬合LBM數(shù)值計(jì)算結(jié)果如圖12（b）所示，其中擬合參數(shù)C=(2.480±0.068)×104μm2，R2=0.98，RMSE=18.31μm2，較改進(jìn)之前的模型（R2=0.96，RMSE=25.06μm2）具有更高的預(yù)測(cè)精度。從上式來(lái)看，k隨孔隙率φ的增加而增加，隨孔隙復(fù)雜度D的增加而減小，與圖10表現(xiàn)出的趨勢(shì)基本一致。

3 結(jié) 論

本文采用X-CT成像技術(shù)和圖像分析方法，獲取了三個(gè)表層有結(jié)皮的土樣，分析了其孔隙結(jié)構(gòu)特征及相關(guān)規(guī)律，結(jié)合單相滲流LBM（Lattice Boltzmann Method，格子玻爾茲曼方法）計(jì)算了各土樣不同深度土層的局部滲透率，在此基礎(chǔ)上對(duì)結(jié)構(gòu)特征參數(shù)與滲透率的相關(guān)性進(jìn)行了分析，并提出了根據(jù)土壤孔隙率和復(fù)雜度計(jì)算滲透率的分形維數(shù)模型，得到了以下結(jié)論：

1）結(jié)構(gòu)分析表明，土樣表層和底層的孔隙率及其孔隙連通性均小于中間土層，三個(gè)土樣各土層的孔隙連通率最大分別為99.21%、94.64%和57.35%，處于土樣中部，而最小值分別為45.99%，27.30%和11.74%，處于土樣的兩端；頂層土孔隙結(jié)構(gòu)較相鄰?fù)翆泳哂休^大的差異，變異程度更大，其變異系數(shù)E＞5%。各土層的平均孔隙直徑隨著深度的增加而減小，分別由0.43、0.37和0.50 mm降低至0.15、0.14和0.14 mm；孔隙率更高的土層具有更好的孔隙連通性和更低的孔隙復(fù)雜度，相關(guān)系數(shù)達(dá)0.97。

2）滲透率計(jì)算結(jié)果顯示，土樣表層具有較低的局部滲透率，分別為7.472、0.006和6.960μm2，表明結(jié)皮對(duì)土壤滲流具有阻滯作用；孔隙率越高、孔隙連通性越好、孔隙復(fù)雜度越低的土層滲透能力越強(qiáng)。

3）基于LBM數(shù)值模擬的結(jié)果，將原始K-C方程中的比表面積替換為孔隙復(fù)雜度D，改進(jìn)后的模型對(duì)土壤滲透率具有更高預(yù)測(cè)精度，其決定系數(shù)R2由0.96提升至0.98，而均方根誤差RMSE則由25.06降低至18.31μm2。