基于CT掃描技術的土壤孔隙定量表達優化

趙 玥 韓巧玲 趙燕東

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

基于CT掃描技術的土壤孔隙定量表達優化

趙 玥 韓巧玲 趙燕東

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

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

土壤斷層掃描圖像； 孔隙結構； 圖像處理技術； 孔隙量化

引言

土壤孔隙是指土壤顆粒之間、團聚體之間或團聚體內部的空隙?？紫督Y構包括孔隙數目、大小等幾何形態，其特征決定了土壤水分保持能力與傳導能力，對土壤水分和養分遷移等生態過程產生重大影響[1-4]。因此，對土壤孔隙結構進行量化分析，是從根本上認識土壤內部結構的必要前提。

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

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

1 孔隙量化系統

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

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

圖1 孔隙量化系統流程圖Fig.1 Flow chart of pore quantitative system

1.1 土壤CT圖像處理

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

1.1.1土壤CT圖像預處理

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

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

1.1.2自適應中值濾波

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

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

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

1.1.3迭代最佳閾值法

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

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

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

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

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

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

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

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

1.1.4Canny邊緣檢測

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

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

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

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

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

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

1.2 孔隙參數量化

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

1.2.1孔隙率

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

1.2.2孔隙數量

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

1.2.3孔隙成圓率

孔隙成圓率是指孔隙形態與標準圓形的接近程度，表示的是不規則孔隙的形態特征，其值越接近1，孔隙結構越接近標準圓形?？紫兜男螒B特征對土壤的水分運輸能力和透氣性能具有較大影響。由于成圓率無法直接計算，因此，本實驗采用間接法來求得孔隙成圓率

C=4πA/P2×100%

(1)

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

P——孔隙周長，像素

1.2.4孔隙分形維數

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

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

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

初始情況下，令

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

(2)

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

(ε=1,2,…)

(3)

(ε=1,2,…)

(4)

則“毯子”的體積為

(5)

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

l——圖像長度，像素

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

(6)

根據分形定義，分形表面積滿足關系

A(ε)=Fε2-D

(7)

lgA(ε)=c1lgε+lgF

(8)

式中F——常數D——分形維數

c1——直線斜率

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

c1=2-D

(9)

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

2 實例應用與驗證

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

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

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

2.1 土壤CT圖像濾波

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

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

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

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

2.2 土壤CT圖像二值化

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

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

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

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

2.3 土壤CT圖像邊緣檢測

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

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

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

2.4 土壤孔隙量化分析

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

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

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

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

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

3 討論

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

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

4 結束語

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

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 張慧娟，孫宇瑞,林劍輝,等.不同粗糙度尺度下預測表層土壤孔隙率量化指數比較研究[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 程亞南,劉建立,張佳寶.土壤孔隙結構定量化研究進展[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 楊永輝,武繼承,毛永萍,等.利用計算機斷層掃描技術研究土壤改良措施下土壤孔隙[J].農業工程學報,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 王恩姮,盧倩倩,陳祥偉.模擬凍融循環對黑土剖面大孔隙特征的影響[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掃描技術在我國土工試驗中的應用研究進展[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 華珊,陳研,梁露燾,等.利用基于偏微分方程的圖像濾波技術研究土壤孔隙結構[J].農業工程學報,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].農業機械學報，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 王聰穎,張慧娟,孫宇瑞,等.基于分形理論的土壤粗糙指數與孔隙率映射規律研究[J].農業機械學報,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

中央高?；究蒲袠I務費專項(BLX2015-36)和國家自然科學基金青年基金項目(41501283)

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

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