基于CT掃描技術(shù)的土壤孔隙定量表達優(yōu)化

趙 玥 韓巧玲 趙燕東

(北京林業(yè)大學工學院，北京 100083)

基于CT掃描技術(shù)的土壤孔隙定量表達優(yōu)化

趙 玥 韓巧玲 趙燕東

(北京林業(yè)大學工學院，北京 100083)

現(xiàn)有土壤孔隙量化方法主要通過圖像處理軟件實現(xiàn)孔隙結(jié)構(gòu)的辨識與分析，此類通用的圖像處理軟件或醫(yī)學處理軟件未考慮土壤內(nèi)部物質(zhì)的復雜多變性以及孔隙結(jié)構(gòu)的細小和不規(guī)則性，從而導致孔隙分割精度低進而量化誤差大，為解決這一問題，本文針對土壤CT圖像的特點提出了一種孔隙量化方法。該方法主要包括圖像處理和量化分析兩部分：選用自適應中值濾波算法去除噪聲對孔隙邊緣的影響，并采用迭代最佳閾值法與Canny邊緣檢測算子相結(jié)合的方法，準確識別出土壤孔隙結(jié)構(gòu)及輪廓線；運用數(shù)學統(tǒng)計方法定量研究土壤孔隙率、孔隙數(shù)目、分形維數(shù)、成圓率等幾何指標，用以揭示孔隙結(jié)構(gòu)的復雜性和不規(guī)則性，實現(xiàn)對土壤孔隙的量化分析。最后，以凍融循環(huán)作用下的土壤為應用對象驗證該方法性能。結(jié)果表明，本文方法能精確地定位孔隙輪廓，有效地分割孔隙結(jié)構(gòu)，而且通過多種孔隙幾何指標的量化可揭示出凍融循環(huán)作用對土壤結(jié)構(gòu)的影響，為孔隙幾何特征和空間特征的量化表達奠定了基礎。

土壤斷層掃描圖像； 孔隙結(jié)構(gòu)； 圖像處理技術(shù)； 孔隙量化

引言

土壤孔隙是指土壤顆粒之間、團聚體之間或團聚體內(nèi)部的空隙。孔隙結(jié)構(gòu)包括孔隙數(shù)目、大小等幾何形態(tài)，其特征決定了土壤水分保持能力與傳導能力，對土壤水分和養(yǎng)分遷移等生態(tài)過程產(chǎn)生重大影響[1-4]。因此，對土壤孔隙結(jié)構(gòu)進行量化分析，是從根本上認識土壤內(nèi)部結(jié)構(gòu)的必要前提。

CT掃描技術(shù)是孔隙量化研究的有效技術(shù)手段[5-8]，而孔隙的辨識與分割則是其重要前提。目前，研究者主要采用圖像處理軟件完成土壤孔隙結(jié)構(gòu)的辨識，如文獻[9-11]等利用與CT機配套的醫(yī)學圖像處理軟件獲得土壤孔隙結(jié)構(gòu)，文獻[12-14]采用Image J和Photoshop軟件實現(xiàn)土壤CT圖像的二值化，獲得孔隙結(jié)構(gòu)圖像。其中，醫(yī)學圖像處理軟件集成了針對人體骨骼和器官特點的處理程序，適用于識別具有較大連通域的病灶；而公共圖像處理軟件主要針對前景和背景具有較大差異的圖像，實現(xiàn)目標的識別。而土壤CT圖像中孔隙細小不規(guī)則(孔徑多為毫米級)，與固相物質(zhì)的灰度值差異也較小，受部分容積效應的影響在圖像不同區(qū)域的灰度值不等，因此，采用現(xiàn)有圖像處理軟件實現(xiàn)土壤CT圖像中孔隙的精確辨識具有較大難度，容易錯誤判別部分區(qū)域的孔隙結(jié)構(gòu)，導致出現(xiàn)孔隙過分割和欠分割現(xiàn)象，使得實驗結(jié)果無法真實反映孔隙的性質(zhì)，從而對后續(xù)孔隙量化研究產(chǎn)生不可避免的影響。

為解決上述缺陷，本文將通過理論分析與實驗比較兩方面的研究，探討現(xiàn)有的圖像處理手段在土壤CT圖像中的表現(xiàn)，并結(jié)合孔隙率、孔隙分形維數(shù)等量化指標提出一種針對土壤CT圖像的孔隙量化方法，該方法針對土壤CT圖像前景孔隙細小、連通性低與背景固相物質(zhì)差異不一致等特點完成孔隙結(jié)構(gòu)的高效分離，將為土壤孔隙的量化研究奠定技術(shù)基礎。為驗證所選方法的有效性，以凍融循環(huán)土壤為應用對象，通過孔隙量化結(jié)果的分析，揭示出土壤凍融循環(huán)對孔隙結(jié)構(gòu)產(chǎn)生的影響。

1 孔隙量化系統(tǒng)

高精度土壤CT圖像是獲取孔隙幾何特征和分布狀態(tài)的數(shù)據(jù)基礎，由于人眼對于灰度值的不敏感性，僅通過目視法無法得到孔隙的有效信息；而且原始土壤CT圖像在獲取和傳輸過程中易受到噪聲、偽影等干擾，無法直接用于孔隙結(jié)構(gòu)的量化研究[15]。因此，借助計算機圖像處理技術(shù)加深對CT圖像深層信息的理解是獲得精確實驗結(jié)果的必要前提。

基于上述原因，本文提出一種針對土壤CT圖像的孔隙量化系統(tǒng)，用以深入探索孔隙幾何特性。該系統(tǒng)主要分為2層：第1層為土壤CT圖像的處理，首先掃描土壤樣本獲取土壤CT圖像，這是提取孔隙有效信息的前提；其次，針對土壤CT圖像特征進行相應的處理，主要選用自適應中值濾波、迭代最佳閾值和Canny邊緣檢測等算法實現(xiàn)，從而得到用于孔隙量化分析的土壤二值圖像，為后續(xù)孔隙量化分析提供精準的數(shù)據(jù)。第2層為孔隙參數(shù)的量化，主要運用數(shù)學統(tǒng)計的方法獲得土壤孔隙率、孔隙數(shù)目、成圓率、分形維數(shù)等指標，完成對土壤孔隙的量化分析。兩層是逐漸遞進的，前者是后者的基礎，后者是前者的目的。如圖1所示。

圖1 孔隙量化系統(tǒng)流程圖Fig.1 Flow chart of pore quantitative system

1.1 土壤CT圖像處理

由于CT機的技術(shù)限制，原始土壤CT圖像邊界處易產(chǎn)生畸變，因此，需對其進行剪裁、存儲等預處理，以提取用于計算機處理分析的土壤有效區(qū)域。但預處理階段并沒有去除噪聲和偽影對圖像的干擾，仍需通過自適應中值濾波的方法來抑制無用信息，突出有效信息，以便于對土壤孔隙結(jié)構(gòu)進行精確分割。本文孔隙結(jié)構(gòu)分割主要涉及到選用迭代最佳閾值法將土壤固體物質(zhì)與孔隙結(jié)構(gòu)分離，并通過Canny邊緣檢測算子準確提取孔隙輪廓線，從而得到只有孔隙結(jié)構(gòu)的二值圖像，為后續(xù)孔隙參數(shù)的量化分析提供精確的數(shù)據(jù)基礎。

1.1.1土壤CT圖像預處理

原始土壤CT圖像包含兩部分信息：第一部分，如圖2a周邊黑色區(qū)域所示，為CT機器的掃描時間、電壓、窗寬和窗位等基本信息。另一部分為圖2a圓圈部分的土壤有效圖像，直徑約為100 mm。由于掃描基本信息與后續(xù)研究無關，且矩形圖像更易于計算機的處理，本文基于圓的內(nèi)切正方形算法對原始土壤CT圖像進行剪裁處理，保留正方形目標區(qū)域(圖2a正方形)，并存儲為如圖2b所示的易于計算機理解的圖像格式，用于土壤孔隙的量化研究。

圖2 土壤CT圖像預處理Fig.2 Preprocessing of soil CT images

1.1.2自適應中值濾波

由于CT機本身放射束硬化、部分容積效應等的約束，導致土壤CT圖像在成像、傳輸?shù)冗^程中存在噪聲，在圖像上表現(xiàn)為與孔隙邊緣有一定的相似性，會影響孔隙輪廓線的提取，因此，對圖像進行濾波處理是準確獲取孔隙結(jié)構(gòu)的前提。

現(xiàn)有研究中，濾波方法主要有均值濾波、維納濾波、中值濾波等，其中均值濾波是以鄰域內(nèi)均值代替原像素值，會模糊孔隙與固相物質(zhì)的邊界，破壞土壤圖像中原始細節(jié)信息；而維納濾波是基于最小均方誤差準則的線性濾波方法，無法去除土壤CT圖像中非平穩(wěn)的隨機噪聲，因此，這兩種方法均不適用于土壤CT圖像的濾波處理。而中值濾波采用模板中各像素點的中值代替原始像素值，能有效消除孤立噪聲點，較好保留了土壤大孔隙邊緣[16]，但由于恒定模板使得濾波效果和細節(jié)信息保護相矛盾，因此，中值濾波會去除圖像中土壤小孔隙，不適用于細小孔隙研究。在此基礎上，本文采用自適應窗寬的濾波算法，實現(xiàn)在消除圖像噪聲的同時保留孔隙的邊緣和細節(jié)信息。

自適應中值濾波采用可變化的濾波窗，通過判斷濾窗中心像素是否為脈沖噪聲，進行不同的輸出操作，實現(xiàn)對圖像的濾波處理。算法主要包括噪聲檢測、濾窗選擇和濾波輸出3部分。實現(xiàn)過程為：若一定大小的濾波窗口Sij內(nèi)中值濾波值Zmed與該點像素值均介于最大濾波值Zmax與最小濾波值Zmin之間，則表明該像素點與中值點都不是噪聲，直接輸出該點灰度Zij，以最大程度保留原始信息；若該像素點是噪聲，則選擇濾窗內(nèi)中值濾波值Zmed輸出。在這一算法中，采用尺寸變化的濾窗是為了針對濾窗中不同類型的像素點分別進行處理，從而實現(xiàn)在保留孔隙細節(jié)信息的同時去除高低脈沖噪聲的干擾。

1.1.3迭代最佳閾值法

圖像濾波后，孔隙邊界較為清晰，與固體物質(zhì)的區(qū)別更為明顯，便于孔隙的提取。研究者普遍通過手動設定固定閾值的方法提取土壤孔隙，如PIERRET等[17]、 MONGA等[18]均通過該方法對孔隙的大小、數(shù)量等信息進行量化分析。手動設置閾值的方法需要根據(jù)具體的圖像設定相應閾值，使得圖像處理的效率較低，也不具備普遍適用性。另外，該方法中閾值的選擇是通過人眼觀察決定的，而人眼對于灰度的敏感性要遠遠低于計算機[19]，因此，針對土壤CT圖像完成閾值的自動選取，是實現(xiàn)對土壤孔隙信息自動判讀的關鍵。考慮到算法的實時性與精準性，本文選用基于逼近思想的迭代最佳閾值法，實現(xiàn)計算機自動選擇最佳閾值，以完成孔隙的精確分割，為孔隙輪廓線的提取奠定技術(shù)基礎。

迭代最佳閾值法主要思想：設定初始閾值，然后按照一定方式對其進行改進，直到滿足設定的準則，迭代結(jié)束，所得數(shù)值即為完成孔隙分割的最佳閾值。其具體實現(xiàn)過程為：

(1)由圖像最大灰度與最小灰度求得初始閾值T0。

(2)根據(jù)初始閾值T0將圖像分割為固體區(qū)域與孔隙區(qū)域，分別求得兩區(qū)域的平均灰度T1和T2。

(3)通過計算步驟(2)中T1和T2均值的方式確定迭代閾值T。

(4)比較迭代閾值T與初始閾值T0的差值d是否在允許范圍內(nèi)，若是，則認為該迭代閾值可以將圖像中固體顆粒與孔隙目標有效地分離；相反，則將迭代閾值T賦予初始閾值T0，繼續(xù)上述步驟(2)和(3)的迭代計算，直到新舊閾值差滿足迭代條件。

(5)最終迭代閾值則為圖像進行二值化處理的最佳閾值。

迭代最佳閾值法原理簡單，易于實現(xiàn)，能自動找到每幅圖像的精確閾值，對噪聲具有一定的魯棒性，可以準確地將土壤孔隙結(jié)構(gòu)分離，適用于大批量的土壤CT圖像。

1.1.4Canny邊緣檢測

基于土壤CT圖像，提取孔隙的輪廓線，有助于精確計算孔隙結(jié)構(gòu)的周長、面積等參數(shù)。邊緣檢測算子主要為基于微分法的一階微分和二階微分算子，其中Robert、Prewitt、Sobel等一階算子采用代表梯度和方向的2個模板對像素點進行卷積運算，將最大值作為邊緣點輸出，而Log等二階算子將二階導數(shù)的過零點判為邊緣。Robert算子適用于邊緣陡峭的低噪聲圖像，Prewitt和Sobel算子的平滑操作會模糊圖像邊緣，適用于圖像噪聲小的圖像，Log算子在邊緣檢測和抑制噪聲方面存在矛盾。而在土壤CT圖像的應用中，需要一種既能有效抑制圖像中的噪聲，又能精確檢測出孔隙細小邊緣的算法。而Canny邊緣檢測算子根據(jù)對信噪比與定位乘積進行測度，得到最優(yōu)化逼近算子，能滿足孔隙輪廓線提取對定位精度和邊緣響應的要求，實現(xiàn)對孔隙邊緣的精準檢測，是最適用于土壤CT圖像特點的算子，其實現(xiàn)過程如下：

(1)利用高斯濾波器對原圖像進行濾波。

(2)通過在鄰域內(nèi)求有限差分來計算圖像各點的梯度幅值和梯度方向。

(3)對梯度幅值進行非極大值抑制，剔除非極大值像素點。

(4)設置高低閾值，并結(jié)合連接分析方法，確定圖像邊緣。

Canny算子通過設定高低閾值，能夠避免邊緣信息丟失或者偽邊緣的情況，從而準確檢測出孔隙邊緣。綜上所述，采用Canny算子進行孔隙邊緣的檢測，能夠得到精確的孔隙圖像，不僅為孔隙結(jié)構(gòu)的深入研究提供技術(shù)指導，也為孔隙參數(shù)的精確量化奠定數(shù)據(jù)基礎。

1.2 孔隙參數(shù)量化

二值化土壤CT圖像顯示出孔隙的幾何形態(tài)及分布，但由于人眼觀察獲取到的孔隙信息有限，因此，需借助計算機對孔隙信息進行深入挖掘，從而實現(xiàn)參數(shù)的量化分析。由于圖像在計算機中是以矩陣的形式存儲的，基于二值化圖像對孔隙參數(shù)進行量化分析的實質(zhì)就是對矩陣中不同像素點運用數(shù)學的方法進行分類與統(tǒng)計，而孔隙率、孔隙數(shù)目、孔隙周長、孔隙面積、成圓率以及分形維數(shù)這6個指標分別從不同方面表現(xiàn)了土壤中孔隙的拓撲結(jié)構(gòu)及分布，對土壤的物理性質(zhì)具有重要的影響。

1.2.1孔隙率

土壤孔隙率是指土顆粒間的孔隙體積占土壤總體積的比率，表征了土壤中自由水和空氣的體積，是評價土壤質(zhì)量和農(nóng)作物生長環(huán)境的重要指標。由于在圖像中，二值化圖像的灰度值只有2個狀態(tài)：代表土壤孔隙的黑色區(qū)域(像素為0)和代表土壤固態(tài)物質(zhì)的白色區(qū)域(像素為1)，因此，黑色區(qū)域(像素值為0)的像素點占總像素點的比重，即為該土壤圖像的孔隙率。

1.2.2孔隙數(shù)量

孔隙數(shù)量即圖像中所包括的孔隙數(shù)目，在二值圖像上表現(xiàn)為像素值為0的集合個數(shù)。由于孔隙的大小、形態(tài)各異，每個孔隙所包含的像素點數(shù)目不同，因此，本文基于連通域的方法來計算孔隙的數(shù)量。其基本思想為若中心點像素值與周圍區(qū)域的8個像素值相同，則稱這些像素點是連通的，將其判斷為一個連通域，圖像中連通域的數(shù)量即為孔隙數(shù)量。

1.2.3孔隙成圓率

孔隙成圓率是指孔隙形態(tài)與標準圓形的接近程度，表示的是不規(guī)則孔隙的形態(tài)特征，其值越接近1，孔隙結(jié)構(gòu)越接近標準圓形。孔隙的形態(tài)特征對土壤的水分運輸能力和透氣性能具有較大影響。由于成圓率無法直接計算，因此，本實驗采用間接法來求得孔隙成圓率

C=4πA/P2×100%

(1)

式中A——孔隙面積，像素×像素

P——孔隙周長，像素

1.2.4孔隙分形維數(shù)

土壤孔隙是細小的不規(guī)則結(jié)構(gòu)，要對其復雜程度進行準確的描述，找到定量的指標是必要的。孔隙分形維數(shù)描述的是土壤孔隙的自相似特性，是土壤孔隙不規(guī)則性的綜合體現(xiàn)[20]，可用于表征土壤孔隙結(jié)構(gòu)的復雜性。

由于土壤CT圖像是由像素點組成的灰度圖像，每個像素點都有相應的灰度值，表示為灰度函數(shù)Z=f(i,j)，在三維空間中，灰度函數(shù)可等價為一個曲面，因此，本文利用雙毯法計算孔隙分形維數(shù)。雙毯法通過求取以灰度曲面為中心的具有一定高度的立方體體積，來估計灰度曲面的分形維數(shù)，并且只適用于求取圖像表面分形維數(shù)。其基本思想如下：

將土壤CT圖像的灰度函數(shù)f(i,j)視為1個曲面，其中(i,j)為圖像上的坐標，選定灰度曲面上的高度ε處，分別建立覆蓋該曲面的毯子，上毯子表示為u0(i,j)，下毯子表示為b0(i,j)，則構(gòu)成1個厚度為2ε的“毯子”，毯子的表面積為體積除以2ε。對于不同的高度ε，可分別計算出分形表面積。

初始情況下，令

u0(i,j)=b0(i,j)=f(i,j)

(2)

上下兩張“毯子”分別按如下原則生長

(ε=1,2,…)

(3)

(ε=1,2,…)

(4)

則“毯子”的體積為

(5)

式中w——圖像寬度，像素

l——圖像長度，像素

對于二維曲面，當毯子高度由ε-1增至ε時，曲面表面積為

(6)

根據(jù)分形定義，分形表面積滿足關系

A(ε)=Fε2-D

(7)

lgA(ε)=c1lgε+lgF

(8)

式中F——常數(shù)D——分形維數(shù)

c1——直線斜率

對于不同的ε，可以計算出不同的lgA(ε)，通過最小二乘法擬合直線斜率，可以計算出c1，由

c1=2-D

(9)

即可計算出分形維數(shù)D，作為表征土壤圖像表面紋理不規(guī)則性和復雜度的參數(shù)。

2 實例應用與驗證

為測試本文選用的基于CT掃描技術(shù)的孔隙量化方法在孔隙結(jié)構(gòu)提取與孔隙參數(shù)量化上的應用效果，以對孔隙結(jié)構(gòu)影響較為強烈的凍融循環(huán)土壤為研究對象進行結(jié)果分析，并通過不同方法的對比實驗，證明本文方法的優(yōu)越性。

本文圖像為由飛利浦Brilliance 64排螺旋CT機掃描土壤樣本所得的CT圖像。掃描樣本分別為經(jīng)歷0、1、3、6、9次凍融循環(huán)的土壤樣本，每次凍融循環(huán)選取3個重復樣本，單個樣本可得220張掃描圖像，單次凍融循環(huán)共有660幅掃描圖像，因此，本實驗的土壤圖像數(shù)據(jù)庫共包括3 300幅土壤CT圖像。

根據(jù)圖像中土壤有效面積的位置, 采用最大內(nèi)切正方形的方法，將每張土壤CT圖像剪切為211像素×211像素后進行實驗。

2.1 土壤CT圖像濾波

為避免土壤CT圖像中的噪聲干擾，精確后續(xù)的邊緣化檢測與定量分析，本實驗采取常用的中值濾波、均值濾波和維納濾波3種濾波方法對土壤CT圖像進行濾波處理，與本文自適應中值濾波效果進行對比，并將實驗結(jié)果直觀地顯示在圖3中。由于篇幅有限，圖3為隨機選取土壤圖像數(shù)據(jù)庫中一幅典型圖像。

圖3 土壤CT圖像的濾波處理效果Fig.3 Filtering effects of soil CT image

由圖3可以看出，相較于圖3a所示的原圖，圖3b和圖3d所示的中值濾波算法和維納濾波算法較好地保留了土壤孔隙的邊緣，但丟失了圖像的細節(jié)信息；圖3c所示的均值濾波算法平滑了土壤孔隙的邊緣，使圖像整體信息變得模糊。上述3種濾波算法雖然能夠有效去除噪聲，但是對于孔隙邊緣的檢測和細節(jié)信息的保護不盡人意，不適用于土壤CT圖像孔隙結(jié)構(gòu)的精確提取。而本文的自適應中值濾波算法，在去除噪聲污染的基礎上，完整保留了土壤孔隙的弱邊緣和細節(jié)信息，如圖3e所示。

經(jīng)實驗證明，本文采用的自適應中值濾波算法，能最大程度地保證孔隙結(jié)構(gòu)，有利于進行孔隙結(jié)構(gòu)的分割處理。

2.2 土壤CT圖像二值化

目前常用于土壤孔隙研究的方法為全局固定閾值法，該方法主要通過閾值的調(diào)節(jié)，選取使得孔隙直徑與真值最相近的數(shù)值作為二值化的閾值，本文采用此方法提取出的土壤孔隙結(jié)構(gòu)圖如圖4b所示。由圖4b明顯觀察到提取的孔隙數(shù)目偏多，孔隙連通域明顯偏大，這是因為圖像不同區(qū)域的灰度差不同，適合特定區(qū)域的閾值并不適用于整幅圖像，以此為依據(jù)的孔隙結(jié)構(gòu)提取易出現(xiàn)過分割或欠分割的現(xiàn)象，從而使得后續(xù)孔隙參數(shù)量化分析存在偏差。

由上文可知，閾值選取不理想，易導致孔隙連通域失真，而本文選取的迭代最佳閾值法則克服了這一難題。該方法選取的閾值綜合了圖像整體信息，是像素點共同作用的結(jié)果，適用于描述土壤CT圖像各方面信息的變化，采用該方法提取的土壤孔隙結(jié)構(gòu)如圖4c所示。通過與圖4a對比可知，采用迭代最佳閾值法提取孔隙結(jié)構(gòu)的形態(tài)、數(shù)目更加符合實際情況，分割效果更為理想。

圖4 土壤CT圖像二值化結(jié)果Fig.4 Soil CT image binarization results

圖5 土壤CT圖像邊緣檢測結(jié)果Fig.5 Soil CT image edge detection results

2.3 土壤CT圖像邊緣檢測

對孔隙參數(shù)進行量化分析，不僅需要提取孔隙結(jié)構(gòu)，還需明確孔隙輪廓線。本實驗采用經(jīng)典的Roberts、Prewitt、Sobel、Log、Canny 5種邊緣檢測算子對孔隙結(jié)構(gòu)進行邊緣檢測，并隨機選取一幅圖像用于檢測效果的比較分析，檢測結(jié)果如圖5所示。

由圖5可以看出，相較于圖5a，圖5b和圖5c所示的Roberts算子和Prewitt算子對孔隙邊緣定位較準，但提取的孔隙輪廓線不具備完好的連通性。圖5d所示的Sobel算子對孔隙邊緣定位不準，得到的孔隙結(jié)構(gòu)有一定程度的失真。圖5e所示的Log算子平滑過程中使得邊緣具有延展性，導致孔隙結(jié)構(gòu)被明顯放大，孔隙形態(tài)失真。上述4種邊緣檢測算子都能夠提取孔隙輪廓線，但是對于孔隙邊緣的定位和形態(tài)的判斷存在一定缺陷，不符合精確提取土壤CT圖像孔隙輪廓線的要求。而圖5f所示的Canny邊緣檢測算子具有定位精準和單邊緣響應強的優(yōu)點，提取的孔隙邊緣連續(xù)性和清晰度比較理想，還可以檢測出圖像中的弱邊緣，更適合土壤CT圖像的細節(jié)刻畫。

經(jīng)實驗證明，Canny算子能夠準確定位和檢測土壤孔隙的輪廓線，其與迭代最佳閾值法的結(jié)合有利于提取孔隙的真實結(jié)構(gòu)，為孔隙的量化研究提供精確的數(shù)據(jù)基礎。

2.4 土壤孔隙量化分析

有研究表明，季節(jié)性凍融循環(huán)對土壤團聚體結(jié)構(gòu)與物理性質(zhì)的影響較為強烈，而團聚體結(jié)構(gòu)與物理性質(zhì)的改變直觀地表現(xiàn)為土壤孔隙數(shù)量與形態(tài)的變化，因此，為探明凍融循環(huán)對土壤結(jié)構(gòu)的影響，本文基于土壤CT圖像，實現(xiàn)了對孔隙結(jié)構(gòu)的量化研究。本文主要通過對孔隙率、孔隙數(shù)量、孔隙成圓率、分形維數(shù)等參數(shù)進行具體的分析，實現(xiàn)對土壤孔隙結(jié)構(gòu)細致深入的研究，從而加強對土壤孔隙幾何特性和空間特性的了解，進而判斷凍融循環(huán)對農(nóng)業(yè)系統(tǒng)、生態(tài)系統(tǒng)等方面的影響。

基于上述土壤CT圖像數(shù)據(jù)庫，本實驗所得孔隙參數(shù)如表1所示。由于篇幅有限，部分參數(shù)以多幅圖像均值形式表示。

表1 不同凍融次數(shù)孔隙參數(shù)的比較Tab.1 Comparison of pore parameters of different freezing-thawing cycles

由原始土壤CT圖像可知，土壤孔隙表現(xiàn)為CT圖像中不同范圍的黑色區(qū)域，土壤孔隙率則體現(xiàn)為黑色區(qū)域所占圖像比例，由表1比較得出，隨著凍融循環(huán)次數(shù)的增加，孔隙率逐漸減小。孔隙率對凍融循環(huán)次數(shù)的響應主要受兩方面影響：孔隙數(shù)量，如凍融循環(huán)0次所示，孔隙數(shù)量最大，說明黑色像素點覆蓋的區(qū)域較多，像素點占圖像的比例也就大，孔隙率也較高；孔隙面積，如凍融循環(huán)1次所示，雖然孔隙數(shù)量較凍融循環(huán)3次少，但由于孔隙平均面積大，單個孔隙所占像素數(shù)多，其孔隙率仍然較大。土壤孔隙的分形維數(shù)則反映了孔隙的復雜程度，其變化與凍融循環(huán)次數(shù)沒有直接的正負相關關系，但其受孔隙率與孔隙成圓率綜合因素的影響。如比較表1中凍融循環(huán)3次和6次的參數(shù)，雖然前者的孔隙率與孔隙數(shù)量都較大，但是由于其孔隙成圓率較小，兩參數(shù)綜合作用下，使得前者的孔隙分形維數(shù)略小于后者。凍融循環(huán)對孔隙數(shù)量、孔隙面積和孔隙成圓率雖然產(chǎn)生顯著影響，但并沒有明顯的變化規(guī)律。綜上所述，各參數(shù)對凍融循環(huán)次數(shù)不是單一響應的，而是通過互相影響與制約決定了土壤結(jié)構(gòu)與物理性質(zhì)的變化。

由表1數(shù)據(jù)可知，同一圖像的孔隙參數(shù)在不同凍融次數(shù)間的差異較小，因此，研究凍融循環(huán)對土壤孔隙結(jié)構(gòu)的影響需要保證孔隙量化過程中所產(chǎn)生的誤差。而本文方法針對孔隙量化表達的所有環(huán)節(jié)，采取不同處理算法確保高獲取數(shù)據(jù)的精度，從而得到較為真實的孔隙參數(shù)，完成土壤孔隙的定量表達。

3 討論

本文基于CT掃描技術(shù)的孔隙量化系統(tǒng)綜合了3種處理算法，以凍融循環(huán)土壤CT圖像為研究對象進行不同算法實驗結(jié)果的比較，從而選擇適用于土壤孔隙量化表達的算法，增強了孔隙形狀的還原度與孔隙輪廓的精確度。但由于算法的比較主要通過目視法定性分析實現(xiàn)，由此獲得的孔隙參數(shù)未與真值進行對比，缺乏一定的說服力。

由于土壤孔隙結(jié)構(gòu)的不規(guī)則性，目前尚未有獲取孔隙真值的標準方法，難以準確獲取土壤孔隙參數(shù)的真值。但已有研究人員采用不同方法進行孔隙結(jié)構(gòu)的量化分析，因此，研究者可以通過與現(xiàn)有量化方法的對比實驗證明所用方法在孔隙還原與輪廓提取方面的優(yōu)勢。

4 結(jié)束語

土壤孔隙狀況是土壤的物理特性之一，它決定了土壤中空氣和水分的含量，進而影響土壤各種變化進程及作物的生長，是目前評價土壤質(zhì)量優(yōu)劣的指標和研究熱點之一。本研究針對土壤CT圖像特點，提出了一種針對土壤CT圖像處理的孔隙量化系統(tǒng)，該系統(tǒng)通過自適應中值濾波、迭代最佳閾值法與Canny邊緣檢測算子的結(jié)合，提取清晰的土壤孔隙輪廓線，增強了孔隙形狀的還原度與孔隙輪廓的精確度。在此基礎上，通過對像素值的統(tǒng)計分析，獲得土壤孔隙率、孔隙數(shù)目、成圓率、孔隙分形維數(shù)等幾何指標，為后續(xù)土壤孔隙連通性、土壤水分運移狀況分析等研究奠定了數(shù)據(jù)基礎。并經(jīng)實驗證明，本文提出的將數(shù)字圖像處理技術(shù)與高精度土壤CT圖像相結(jié)合的孔隙量化方法，能夠有效解決無法準確區(qū)分具有相似性背景與目標的缺陷，較好地從土壤固相中分離出孔隙結(jié)構(gòu)，準確識別和定位孔隙輪廓線，為土壤孔隙研究提供了一種有效技術(shù)手段，對土壤孔隙幾何特征和空間特征的研究具有重要影響。

1 HILL R L, HORTON R, CRUSE R M.Tillage effects on soil water retention and pore size distribution of two mollisols [J].Soil Science Society of America Journal, 1984, 49(5):1264-1270.

2 張慧娟，孫宇瑞,林劍輝,等.不同粗糙度尺度下預測表層土壤孔隙率量化指數(shù)比較研究[J].應用基礎與工程科學學報,2009,19(1):69-76.

ZHANG Huijuan, SUN Yurui, LIN Jianhui, et al.Comparative study on predicting porosity by different roughness indices [J].Journal of Basic Science and Engineering, 2009, 17(1):69-76.(in Chinese)

3 TOKUMOTO I, NOBORIO K, KOGA K.Coupled water and heat flow in a grass field with aggregated Andisol during soil-freezing periods[J].Cold Regions Science and Technology, 2010, 62(2-3):98-106.

4 STEPHENS P R, HEWITT A E, SPARLING G P, et al.Assessing sustainability of land management using a risk identification model[J].Pedosphere, 2003, 13:41-48.

5 PETROVIC A M, SIEBERT J E, RIEKE P E.Soil bulk density analysis in three dimensions by computed tomographic scanning [J].Soil Science Society, 1982, 46:445-450.

6 程亞南,劉建立,張佳寶.土壤孔隙結(jié)構(gòu)定量化研究進展[J].土壤通報,2012(4):988-994.

CHENG Ya’nan, LIU Jianli, ZHANG Jiabao.Advance in the study on quantification of soil pore structure [J].Chinese Journal of Soil Science, 2012(4):988-994.(in Chinese)

7 ELLIOT T R, HECK R J.A comparison of optical and X- ray CT technique for void analysis in soil thin section [J].Geoderma, 2007, 141(1-2): 60-70.

8 楊永輝,武繼承,毛永萍,等.利用計算機斷層掃描技術(shù)研究土壤改良措施下土壤孔隙[J].農(nóng)業(yè)工程學報,2013,29(23):99-108.

YANG Yonghui, WU Jicheng, MAO Yongping, et al.Using computed tomography scanning to study soil pores under different soil structure improvement measures [J].Transactions of the CSAE, 2013,29(23):99-108.(in Chinese)

9 WAMER G S, NIEBER J L, MOORE I D, et al.Characterizing macropores in soil by computed tomography [J].Soil Science Society of America Journal, 1989, 53(3): 653-660.

10 PEYLON R L, GANTZER C J, ANDERSON S H.Fractal dimension to describe soil macropores structure using X-ray computed tomography [J].Water Resources Research, 1994,30(3): 691-700.

11 CORTINA-JANUCHS M G, QUINTANILLA-DAMINGUEZ J, VEGA-CORONA A, et al.Detection of pore space in CT soil images using artificial neural networks [J].Biogeosciences, 2011, 8(2):279-288.

12 RANJITH P U, STEPHEN H A, CLARK J G, et al.Influence of prairie restoration on CT-measured soil pore characteristics [J].Journal of Environmental Quality, 2008, 37(1):219-228.

13 趙世偉，趙勇鋼，吳金水.黃土高原植被演替下土壤孔隙的定量分析[J].中國科學：地球科學,2010，40(2)：223-231.

ZHAO Shiwei, ZHAO Yonggang, WU Jinshui.Quantitative analysis of soil pores under natural vegetation successions on the Loess Plateau [J].Science of China: Earth Science, 2010, 40(2):223-231.(in Chinese)

14 王恩姮,盧倩倩,陳祥偉.模擬凍融循環(huán)對黑土剖面大孔隙特征的影響[J].土壤學報,2014,51(3):490-496.

WANG Enheng, LU Qianqian, CHEN Xiangwei.Characterization of macro-pores in mollisol profile subjected to simulated freezing-thawing alternation [J].Journal of Soil Science, 2014, 51(3):490-496.(in Chinese)

15 王艷麗,程展林.CT掃描技術(shù)在我國土工試驗中的應用研究進展[J].地震工程學報,2015,37(增刊1):35-39.

WANG Yanli, CHENG Zhanlin.Progress in the application of CT scanning technology in chinese soil tests [J].China Earthquake Engineering Journal, 2015,37(Supp.1):35-39.(in Chinese)

16 華珊,陳研,梁露燾,等.利用基于偏微分方程的圖像濾波技術(shù)研究土壤孔隙結(jié)構(gòu)[J].農(nóng)業(yè)工程學報,2014,30(3):78-85.

HUA Shan, CHEN Yan, LIANG Lutao, et al.Studying soil pore structure by using image filtering technology based on partial differential equation model [J].Transactions of the CSAE, 2014, 30(3): 78-85.(in Chinese)

17 PIERRET A, CAPOWIEZ Y, BELZUNCES L, et al.3D reconstruction and quantification of macropores using computed tomography and image analysis [J].Geoderma, 2002, 106(3-4):247-271.

18 MONGA O, BOUSSOA M, GARNIERB P, et al.3D geometric structures and biological activity:application to microbial soil organic matter decomposition in pore space [J].Ecological Modelling, 2008, 216(3-4):291-302.

19 姬偉,陶云,趙德安,等.基于CLAHE的蘋果樹樹枝迭代閾值分割方法研究[J/OL].農(nóng)業(yè)機械學報，2014,45(4):69-75.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20140411&flag=1.DOI：10.6041/j.issn.1000-1298.2014.04.011.

JI Wei, TAO Yun, ZHAO Dean, et al.Iterative threshold segmentation of apple branch images based on CLAHE [J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2014,45(4):69-75.(in Chinese)

20 王聰穎,張慧娟,孫宇瑞,等.基于分形理論的土壤粗糙指數(shù)與孔隙率映射規(guī)律研究[J].農(nóng)業(yè)機械學報,2011,42(11):32-38.

WANG Congying, ZHANG Huijuan, SUN Yurui, et al.Relationship between soil surface porosity and roughness indices based on fractal theory [J].Transactions of the Chinese Society for Agricultural Machinery, 2011,42(11):32-38.(in Chinese)

OptimizationofSoilPoreQuantitativeExpressionBasedonComputedTomographyScanningTechnology

ZHAO Yue HAN Qiaoling ZHAO Yandong

(SchoolofTechnology,BeijingForestryUniversity,Beijing100083,China)

In recent years, image processing software was wisely applied to identify and analyze pore structure.However, these softwares, such as Photoshop and Image J, did not take into account the complexity of the internal material in the soil and the irregularity of pore structure, and they caused low pore segmentation precision.In order to solve the problem, a pore quantitative method based on the characteristics of soil computed tomography (CT) image was proposed.This method mainly included image processing and quantification analysis.Firstly, the adaptive median filtering algorithm was adopted to remove the effect of image noise on the edge of pore.Then, the method of iterative optimal threshold and canny edge detection was used to identify the pore structure in the soil and the contour line of the pore.Secondly, soil pore structure had evident spatial characteristics, which included soil porosity, pore number, pore radius, spore size distribution, circularity, fractal dimension, and so on.These pore geometry indicators were calculated by using the mathematical statistics method, and they could reveal the complexity and irregularity of pore structure.These geometry indicators were useful for realizing the quantitative analysis of the soil porosity.Finally, the method was applied to the soil under freeze-thaw cycle.The results showed that the method can accurately locate the pore profile, and segment the pore structure effectively.Furthermore, the effect of freezing and thawing cycles on the soil structure was revealed by quantifying the geometrical parameters of various soil pores, thus it proved the effectiveness of the method and laied foundation for quantification of soil pore geometry and spatial characteristics.

soil computed tomography image; pore structure; image processing technology; quantitation of pore

10.6041/j.issn.1000-1298.2017.10.031

S152

A

1000-1298(2017)10-0252-08

2017-01-12

2017-02-14

中央高校基本科研業(yè)務費專項(BLX2015-36)和國家自然科學基金青年基金項目(41501283)

趙玥(1986—)，女，講師，博士，主要從事圖像處理與模式識別等研究，E-mail: zhaoyue0609@126.com

趙燕東(1965—)，女，教授，博士生導師，主要從事生態(tài)信息智能檢測與控制研究，E-mail: yandongzh@bjfu.edu.cn