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

趙 玥 韓巧玲 趙燕東

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

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

趙 玥 韓巧玲 趙燕東

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

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

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

引言

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

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

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

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

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

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

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

1.1 土壤CT圖像處理

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

1.1.1土壤CT圖像預(yù)處理

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

圖2 土壤CT圖像預(yù)處理Fig.2 Preprocessing of soil CT images

1.1.2自適應(yīng)中值濾波

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

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

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

1.1.3迭代最佳閾值法

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

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

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

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

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

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

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

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

1.1.4Canny邊緣檢測(cè)

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

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

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

(3)對(duì)梯度幅值進(jìn)行非極大值抑制，剔除非極大值像素點(diǎn)。

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

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

1.2 孔隙參數(shù)量化

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

1.2.1孔隙率

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

1.2.2孔隙數(shù)量

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

1.2.3孔隙成圓率

孔隙成圓率是指孔隙形態(tài)與標(biāo)準(zhǔn)圓形的接近程度，表示的是不規(guī)則孔隙的形態(tài)特征，其值越接近1，孔隙結(jié)構(gòu)越接近標(biāo)準(zhǔn)圓形?？紫兜男螒B(tài)特征對(duì)土壤的水分運(yùn)輸能力和透氣性能具有較大影響。由于成圓率無(wú)法直接計(jì)算，因此，本實(shí)驗(yàn)采用間接法來(lái)求得孔隙成圓率

C=4πA/P2×100%

(1)

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

P——孔隙周長(zhǎng)，像素

1.2.4孔隙分形維數(shù)

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

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

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

初始情況下，令

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

(2)

上下兩張“毯子”分別按如下原則生長(zhǎng)

(ε=1,2,…)

(3)

(ε=1,2,…)

(4)

則“毯子”的體積為

(5)

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

l——圖像長(zhǎng)度，像素

對(duì)于二維曲面，當(dāng)毯子高度由ε-1增至ε時(shí)，曲面表面積為

(6)

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

A(ε)=Fε2-D

(7)

lgA(ε)=c1lgε+lgF

(8)

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

c1——直線斜率

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

c1=2-D

(9)

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

2 實(shí)例應(yīng)用與驗(yàn)證

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

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

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

2.1 土壤CT圖像濾波

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

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

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

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

2.2 土壤CT圖像二值化

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

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

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

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

2.3 土壤CT圖像邊緣檢測(cè)

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

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

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

2.4 土壤孔隙量化分析

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

基于上述土壤CT圖像數(shù)據(jù)庫(kù)，本實(shí)驗(yàn)所得孔隙參數(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ù)的增加，孔隙率逐漸減小。孔隙率對(duì)凍融循環(huán)次數(shù)的響應(yīng)主要受兩方面影響：孔隙數(shù)量，如凍融循環(huán)0次所示，孔隙數(shù)量最大，說(shuō)明黑色像素點(diǎn)覆蓋的區(qū)域較多，像素點(diǎn)占圖像的比例也就大，孔隙率也較高；孔隙面積，如凍融循環(huán)1次所示，雖然孔隙數(shù)量較凍融循環(huán)3次少，但由于孔隙平均面積大，單個(gè)孔隙所占像素?cái)?shù)多，其孔隙率仍然較大。土壤孔隙的分形維數(shù)則反映了孔隙的復(fù)雜程度，其變化與凍融循環(huán)次數(shù)沒有直接的正負(fù)相關(guān)關(guān)系，但其受孔隙率與孔隙成圓率綜合因素的影響。如比較表1中凍融循環(huán)3次和6次的參數(shù)，雖然前者的孔隙率與孔隙數(shù)量都較大，但是由于其孔隙成圓率較小，兩參數(shù)綜合作用下，使得前者的孔隙分形維數(shù)略小于后者。凍融循環(huán)對(duì)孔隙數(shù)量、孔隙面積和孔隙成圓率雖然產(chǎn)生顯著影響，但并沒有明顯的變化規(guī)律。綜上所述，各參數(shù)對(duì)凍融循環(huán)次數(shù)不是單一響應(yīng)的，而是通過(guò)互相影響與制約決定了土壤結(jié)構(gòu)與物理性質(zhì)的變化。

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

3 討論

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

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

4 結(jié)束語(yǔ)

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

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 張慧娟，孫宇瑞,林劍輝,等.不同粗糙度尺度下預(yù)測(cè)表層土壤孔隙率量化指數(shù)比較研究[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào),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ìn)展[J].土壤通報(bào),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ì)算機(jī)斷層掃描技術(shù)研究土壤改良措施下土壤孔隙[J].農(nóng)業(yè)工程學(xué)報(bào),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].中國(guó)科學(xué)：地球科學(xué),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)對(duì)黑土剖面大孔隙特征的影響[J].土壤學(xué)報(bào),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ù)在我國(guó)土工試驗(yàn)中的應(yīng)用研究進(jìn)展[J].地震工程學(xué)報(bào),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è)工程學(xué)報(bào),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è)機(jī)械學(xué)報(bào)，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è)機(jī)械學(xué)報(bào),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(yè)務(wù)費(fèi)專項(xiàng)(BLX2015-36)和國(guó)家自然科學(xué)基金青年基金項(xiàng)目(41501283)

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

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