基于細(xì)化法的土壤孔隙骨架提取算法研究

韓巧玲 趙 玥 趙燕東 潘賢君 彭 涌 鄭一力

(1.北京林業(yè)大學(xué)工學(xué)院， 北京 100083； 2.城鄉(xiāng)生態(tài)環(huán)境北京實(shí)驗(yàn)室， 北京 100083；3.林業(yè)裝備與自動(dòng)化國家林業(yè)局重點(diǎn)實(shí)驗(yàn)室， 北京 100083)

0 引言

土壤孔隙結(jié)構(gòu)是土壤固相顆粒、土壤有機(jī)物、微生物等各相物質(zhì)在空間排列后形成的空隙[1-3]。其幾何形狀、空間分布、連通性能等三維特征決定了養(yǎng)分運(yùn)移能力、透氣性等土壤功能，對(duì)于土壤生態(tài)過程具有重要影響。因此，對(duì)孔隙三維特征的描述和量化分析是研究土壤功能多樣性的一個(gè)重要方向。

近年來，計(jì)算機(jī)斷層掃描技術(shù)(Computer tomography, CT)在土壤領(lǐng)域的成功應(yīng)用，能夠?qū)崿F(xiàn)土壤孔隙三維結(jié)構(gòu)的可視化[4-5]。由于土壤組成成分的復(fù)雜性，使得孔隙結(jié)構(gòu)具有多種幾何形狀和多尺度特征：①孔隙呈不規(guī)則的幾何形狀，并向多個(gè)方向發(fā)育延伸。②具有不同數(shù)量級(jí)的幾何尺寸。③多尺度特征的孔隙相互交錯(cuò)。復(fù)雜的拓?fù)浣Y(jié)構(gòu)，使得土壤孔隙的研究難度大。為了探索復(fù)雜孔隙結(jié)構(gòu)的體積、曲折度、連通性等特征，需要結(jié)合重構(gòu)技術(shù)對(duì)其空間結(jié)構(gòu)和拓?fù)涮匦赃M(jìn)行分析，在這一過程中，孔隙骨架模型的構(gòu)建方法至關(guān)重要。

現(xiàn)有的骨架化方法中，應(yīng)用較多的主要有廣義勢(shì)場(chǎng)法、Vonoroi圖法、Reeb圖法、距離變換法和細(xì)化法。廣義勢(shì)場(chǎng)法的構(gòu)建精度較高，但算法復(fù)雜度高，運(yùn)行時(shí)間長，不利于實(shí)際應(yīng)用[6]；Voronoi圖和Reeb圖是基于幾何分析獲取的骨架模型，適用于結(jié)構(gòu)相對(duì)簡(jiǎn)單的幾何圖形[7-9]。應(yīng)用范圍最廣泛的距離變換法和細(xì)化法易操作，精度高，適用于不規(guī)則的土壤孔隙，有利于保持孔隙骨架模型的中心性和連通性[10-11]。

本文將采用距離變換法和細(xì)化法構(gòu)建土壤孔隙骨架模型，并結(jié)合土壤孔隙結(jié)構(gòu)的特性增加細(xì)化法的限制條件。試驗(yàn)分別以自定義規(guī)則結(jié)構(gòu)和土壤孔隙結(jié)構(gòu)為應(yīng)用對(duì)象，在連通性、細(xì)化性和中心性3個(gè)層面比較2種算法的構(gòu)建效果，并通過對(duì)中心性的量化分析，證明細(xì)化法對(duì)土壤孔隙骨架模型構(gòu)建的優(yōu)越性，為從孔隙尺度理解土壤功能提供一種較為先進(jìn)的技術(shù)手段。

1 土壤樣本采集與圖像預(yù)處理

1.1 土壤樣本的采集

研究所用土壤樣本采自黑龍江省克山農(nóng)場(chǎng)(125°8′～125°37′E，48°12′～48°23′N)，該地區(qū)屬于典型黑土區(qū)，年平均氣溫1.3℃左右，年降水量約502.5 mm，其土壤類型以黏化濕潤均腐土為主[12]。于2017年10月初，采用自制的內(nèi)徑和高度分別為10 cm的圓柱形有機(jī)玻璃管于土壤深度為0～40 cm的侵蝕溝壁進(jìn)行原狀土取樣。通過機(jī)械分層法，每10 cm取一次土樣，每層重復(fù)取樣3次，共得到12個(gè)土壤樣本，以供后續(xù)CT掃描用。

1.2 土壤CT圖像的獲取

所用的土壤CT圖像由黑龍江省中醫(yī)藥大學(xué)的Philips Brilliance 128排64層螺旋CT機(jī)掃描土壤樣本獲得。在土壤樣本掃描前，為了保證土壤成像的精度，采用標(biāo)準(zhǔn)水膜對(duì)CT機(jī)進(jìn)行校正。掃描過程的參數(shù)選擇為：電壓120 kV，電流196 mA，掃描層厚0.625 mm，窗寬和窗位分別為2 000和800[13]。經(jīng)過螺旋掃描，每個(gè)土壤樣本得到236幅土壤CT圖像。

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

由于醫(yī)療CT機(jī)掃描土壤樣本后得到的是包含256個(gè)灰度級(jí)的DICOM格式圖像，不便于進(jìn)行后續(xù)的圖像處理操作，因此，為保證土壤CT圖像的可操作性，將原始土壤CT圖像存儲(chǔ)為BMP格式。圖1a所示為原始土壤CT圖像，圓圈內(nèi)為有效的土壤圖像，為了減少算法占用內(nèi)存和提高圖像處理的執(zhí)行效率，第1步需要對(duì)原始土壤CT圖像進(jìn)行剪裁處理。另外，考慮到土壤圓柱與PVC玻璃管邊界處會(huì)存在一定晃動(dòng)，可能會(huì)導(dǎo)致土壤邊界存在一定畸變的情況，因此，本試驗(yàn)采用最大內(nèi)切正方形的方法，將土壤CT圖像剪裁成尺寸為289像素×289像素的圖像，如圖1b所示。

第2步，基于剪裁后的土壤CT圖像進(jìn)行孔隙結(jié)構(gòu)辨識(shí)。而由于CT掃描儀受電路板等硬件的限制和掃描層厚等參數(shù)的約束，會(huì)導(dǎo)致土壤CT圖像受到噪聲的污染。因此，選用自適應(yīng)中值濾波算法去除圖像噪聲的干擾，增強(qiáng)孔隙結(jié)構(gòu)信息[14-15]?；趫D1c所示的濾波圖像，本文采用自適應(yīng)模糊C均值算法完成孔隙結(jié)構(gòu)的二值化[16]，結(jié)果如圖1d所示，黑色結(jié)構(gòu)為土壤孔隙。

圖1 土壤CT圖像Fig.1 Soil CT images

基于孔隙二值化圖像進(jìn)行三維重構(gòu)，是構(gòu)建孔隙骨架模型的最后一個(gè)預(yù)處理環(huán)節(jié)。為了保證三維模型的真實(shí)性和消除土壤樣本由于移動(dòng)、碰撞造成的畸變影響，本文選取土壤樣本中間部位的216幅圖像，采用光線投影法完成孔隙三維模型的重構(gòu)[17]。如圖1e所示，重構(gòu)后模型的尺寸為289像素×289像素×216幅，可為后續(xù)土壤孔隙骨架模型的構(gòu)建提供精準(zhǔn)的數(shù)據(jù)基礎(chǔ)。

2 研究方法

2.1 距離變換法

距離變換法是一種基于局部區(qū)域進(jìn)行的賦值變換，其主要思想是計(jì)算并標(biāo)記空間目標(biāo)體素點(diǎn)到模型邊界體素點(diǎn)的距離，并依據(jù)此距離完成骨架模型的構(gòu)建[18]。對(duì)于三維離散圖像，體素點(diǎn)(k,j,i)的距離變換值Dv的表達(dá)式為

(1)

(2)

式中d——體素點(diǎn)(k,j,i)到體素點(diǎn)(x,y,z)的歐氏距離

B——所有邊界體素點(diǎn)的集合

v1、v2——體素點(diǎn)

距離變換法主要包括3個(gè)步驟：①采用式(1)和式(2)計(jì)算孔隙模型中所有體素點(diǎn)的距離變換值，組成距離變換矩陣。②選取距離變換矩陣中數(shù)值最大的體素點(diǎn)作為骨架模型的種子點(diǎn)，由式(1)、(2)可知，三維目標(biāo)的中心體素點(diǎn)具有最大的距離變換值。③以種子點(diǎn)為起點(diǎn)，確定其26 鄰域內(nèi)距離變換具有局部極大值的體素點(diǎn)，該體素點(diǎn)即為下一個(gè)骨架點(diǎn)。重復(fù)步驟③，直至符合條件的體素點(diǎn)不再發(fā)生變化，連接所有體素點(diǎn)，則為孔隙骨架模型。

2.2 細(xì)化法

細(xì)化法的基本思想是通過逐步均勻地剝蝕模型中邊界體素點(diǎn)，以獲得連通性不變的骨架模型[19]。該方法運(yùn)算過程由2部分組成：①細(xì)化條件的建立，這是構(gòu)建骨架模型的核心問題和必備條件。②偽分支和毛刺的處理，這是保證孔隙骨架模型連通性、中心性和細(xì)化性的重要因素。

2.2.1細(xì)化條件的建立

本文基于ROCKETT[20]提出的單像素寬度細(xì)化算法，建立孔隙結(jié)構(gòu)26鄰域內(nèi)的鄰接矩陣，通過判斷鄰域內(nèi)的中心體素點(diǎn)是否會(huì)影響骨架結(jié)構(gòu)的連通性，來決定是否從現(xiàn)有骨架模型中刪除該體素點(diǎn)，并通過不斷迭代判斷實(shí)現(xiàn)對(duì)骨架的剝蝕。在這一過程中，細(xì)化條件的建立需要保證骨架模型的連通性、細(xì)化性和中心性，考慮到土壤孔隙的特性，本文細(xì)化法共包含兩個(gè)條件：

(1)為了保證骨架模型的連通性，所有刪除的目標(biāo)點(diǎn)都應(yīng)保證在其26鄰域內(nèi)至少有2個(gè)相同的體素點(diǎn)。按照該條件進(jìn)行模型遍歷，直到生成的骨架模型中目標(biāo)體素點(diǎn)的數(shù)目不再改變。

(2)由于孔隙結(jié)構(gòu)細(xì)小不均勻的特性，圖像噪聲會(huì)引起骨架模型具有冗余的枝節(jié)，因此，在細(xì)化法中增加了一個(gè)基于閾值的限制條件。骨架閾值長度T設(shè)置為5個(gè)體素點(diǎn)。遍歷所有骨架結(jié)構(gòu)，將長度小于該閾值的骨架刪除，以保證模型的細(xì)化性。

2.2.2偽分支和毛刺的處理

骨架化過程中，細(xì)化法僅考慮了模型的局部連通性，易導(dǎo)致初始骨架模型存在冗余的分支、毛刺以及缺口的現(xiàn)象。因此，為了保證孔隙骨架模型的連通性和細(xì)化性，需要對(duì)骨架模型進(jìn)行去毛刺處理。

去偽分支和毛刺處理主要包括毛刺的刪除和狹窄缺口的填補(bǔ)，其實(shí)現(xiàn)過程主要包括3個(gè)步驟：①依據(jù)某一準(zhǔn)則判斷出組成孔隙骨架中軸線的體素點(diǎn)。②以這些體素點(diǎn)組成的中心路徑為標(biāo)準(zhǔn)，刪除路徑之外的冗余體素點(diǎn)，以去除骨架上的毛刺。③通過形態(tài)學(xué)操作，填補(bǔ)骨架路徑的狹窄缺口，從而保證孔隙骨架模型的拓?fù)洳蛔冃浴?/p>

在這一修正過程中，中心路徑的確定對(duì)于骨架模型去毛刺處理的精度具有重大影響。本文通過3個(gè)步驟建立中心路徑：①在連通域的兩端分別指定一個(gè)起始點(diǎn)和一個(gè)終止點(diǎn)。②由起始點(diǎn)進(jìn)行路徑漫游，判斷在26連通域內(nèi)能達(dá)到終點(diǎn)的所有路徑。③所有路徑中的最短路徑即為中心路徑。以此為標(biāo)準(zhǔn)，可以快速精確地完成骨架模型的去毛刺處理。

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

為測(cè)試距離變換法和細(xì)化法在土壤孔隙骨架模型構(gòu)建的可行性和精確性，選用2組自定義的規(guī)則結(jié)構(gòu)、1組基于3D打印技術(shù)和CT掃描技術(shù)獲得的規(guī)則圓柱體模型和12組基于二維土壤CT圖像生成的孔隙模型為研究對(duì)象。通過對(duì)規(guī)則結(jié)構(gòu)骨架模型和土壤孔隙骨架模型構(gòu)建結(jié)果的比較，從細(xì)化性、中心性和連通性3個(gè)層面分析2種算法的骨架化性能。

為消除運(yùn)行環(huán)境和計(jì)算機(jī)硬件對(duì)2種算法性能的影響，保證2種算法的可比性，所有試驗(yàn)均在同一臺(tái)計(jì)算機(jī)上采用Matlab R2018a實(shí)現(xiàn)，計(jì)算機(jī)處理器為4.00 GHz Intel Core i7 4790，內(nèi)存16 GB，操作系統(tǒng)為Windows 7。

3.1 規(guī)則結(jié)構(gòu)的骨架模型

為詳細(xì)比較2種算法對(duì)規(guī)則結(jié)構(gòu)骨架模型的構(gòu)建效果，選取了規(guī)則曲線、規(guī)則圓環(huán)和規(guī)則圓柱體進(jìn)行骨架模型的構(gòu)建。如圖2所示，灰色區(qū)域?yàn)樵冀Y(jié)構(gòu)，紅色曲線為提取的骨架模型。

由圖2a、2b可知，距離變換法針對(duì)規(guī)則曲線和規(guī)則圓環(huán)構(gòu)建的骨架模型具有明顯的斷裂，連通性較差；而且通過與原始結(jié)構(gòu)的對(duì)比發(fā)現(xiàn)，其構(gòu)建的骨架模型長度偏短，主要體現(xiàn)在骨架模型的兩端(藍(lán)色方框所示)。相較于距離變換法，細(xì)化法提取出的骨架結(jié)構(gòu)沒有出現(xiàn)明顯的斷裂，而且在首端和尾端也沒有出現(xiàn)骨架缺失的情況，具有較好的連通性。

由圖2c、2f藍(lán)色方框可知，對(duì)于規(guī)則圓柱體模型，兩種算法構(gòu)建的骨架結(jié)構(gòu)兩端均沒有出現(xiàn)明顯的缺失現(xiàn)象。但是，在圓柱體直徑較小時(shí)，距離變換法構(gòu)建的孔隙骨架存在多處斷裂的情況，而細(xì)化法構(gòu)建的骨架模型則在很大程度上改善了這一問題。綜上所述，細(xì)化法構(gòu)建的骨架模型與真實(shí)情況更為一致。

3.2 土壤孔隙骨架模型

圖3為2種算法針對(duì)土壤孔隙骨架模型的構(gòu)建結(jié)果。為詳細(xì)比較2種算法對(duì)土壤孔隙骨架模型的構(gòu)建細(xì)節(jié)，隨機(jī)選取了原狀土壤樣本中的兩個(gè)大孔隙結(jié)構(gòu)進(jìn)行比較分析。圖3a、3b、3d、3e所示為兩個(gè)大孔隙A、B骨架模型的放大圖。

由圖3a、3b、3d、3e可以看出，距離變換法和細(xì)化法都能得到較為清晰的土壤孔隙骨架模型，并分列于相應(yīng)原孔隙的中心。但是距離變換法出現(xiàn)斷裂的情況，破壞了骨架的連通性，基于該孔隙骨架模型不利于對(duì)孔隙特征的量化分析。這一原因主要是因?yàn)榫嚯x變換法在局部區(qū)域可能沒有絕對(duì)極大值，從而難以準(zhǔn)確判定骨架體素點(diǎn)的位置，導(dǎo)致體素點(diǎn)偏離中心。相比較而言，細(xì)化法能夠定義骨架模型的簡(jiǎn)單點(diǎn)并進(jìn)行循環(huán)判斷，保證骨架模型具有較好的連通性；并通過去偽分支和毛刺的處理避免了體素點(diǎn)重復(fù)和骨架存在狹窄缺口的情況，保證了骨架的體素點(diǎn)處于孔隙的中心位置。

如圖3c、3f所示，在真實(shí)土壤孔隙骨架模型中，單個(gè)體素點(diǎn)之間存在空隙，這是由于圖像尺度增大導(dǎo)致的視覺效果。由圖3c、3f的放大圖可知，距離變換法仍存在大量骨架結(jié)構(gòu)缺失，對(duì)于部分體素點(diǎn)也出現(xiàn)了偏移中心位置的現(xiàn)象；而且局部區(qū)域也出現(xiàn)了多個(gè)體素點(diǎn)重疊的情況。通過對(duì)整個(gè)骨架模型的分析，發(fā)現(xiàn)距離變換法的細(xì)化性和連通性較差。反之，細(xì)化法的骨架結(jié)構(gòu)完整，具有單體素點(diǎn)的特征，也處于土壤孔隙的中心位置，具有良好的連通性、細(xì)化性和中心性。

3.3 骨架模型中心性評(píng)價(jià)

為了更精確評(píng)價(jià)2種骨架化算法的性能，本文采用中心性作為定量評(píng)價(jià)標(biāo)準(zhǔn)。由于該指標(biāo)均需要基于標(biāo)準(zhǔn)骨架結(jié)構(gòu)進(jìn)行分析，而土壤孔隙的真實(shí)骨架模型難以確定，因此，試驗(yàn)選用CT掃描的規(guī)則圓柱體模型進(jìn)行量化試驗(yàn)的測(cè)試對(duì)象。

中心性表示實(shí)際的骨架結(jié)構(gòu)S(i)與構(gòu)建的骨架模型S′(i)之間的偏差程度。根據(jù)歐幾里德理論，中心性Cd的計(jì)算公式為

(3)

式中K——骨架的數(shù)量

i——骨架結(jié)構(gòu)的序號(hào)

由式(3)可知，中心性表示構(gòu)建的骨架模型偏離原始結(jié)構(gòu)中心位置的距離，其值越小，構(gòu)建的骨架模型與真實(shí)的骨架模型越接近，則說明算法的中心性越好。圖4為2種骨架模型構(gòu)建算法的中心性評(píng)價(jià)結(jié)果。

圖4 圓柱體模型中心性評(píng)價(jià)結(jié)果Fig.4 Evaluation results of centrality

由圖4可以看出，兩種算法的骨架偏移距離和圓柱體的孔徑?jīng)]有明顯的關(guān)系。不論圓柱體的孔徑多大，細(xì)化法的骨架偏移距離一直小于距離變換法。細(xì)化法平均骨架偏移距離為0.10 mm，距離變換法則為0.15 mm。骨架偏移距離越接近0，表示中心性越好，因此，細(xì)化法的平均中心性優(yōu)于距離變換法。此外，細(xì)化法的方差為0.05 mm，僅為距離變換法(0.12 mm)的2/5倍，由此可知，細(xì)化法針對(duì)不同直徑的孔隙結(jié)構(gòu)具有更為穩(wěn)定的性能和更強(qiáng)的魯棒性。

結(jié)合規(guī)則幾何結(jié)構(gòu)和土壤孔隙結(jié)構(gòu)的試驗(yàn)結(jié)果可知，細(xì)化法構(gòu)建的骨架模型具有良好的細(xì)化性、連通性和中心性，能夠更好地保留土壤孔隙結(jié)構(gòu)的拓?fù)涮匦?，可為土壤孔隙結(jié)構(gòu)的空間分布研究提供一種先進(jìn)的技術(shù)手段。

4 結(jié)論

(1)針對(duì)自定義規(guī)則結(jié)構(gòu)和土壤孔隙結(jié)構(gòu)的骨架模型構(gòu)建試驗(yàn)可得，距離變換法和細(xì)化法的骨架模型構(gòu)建效果均不受模型類型的影響。其中，距離變換法構(gòu)建的骨架模型中存在部分體素點(diǎn)缺失和偏移的現(xiàn)象，使得骨架的實(shí)際長度小于真實(shí)情況。而細(xì)化法構(gòu)建的骨架結(jié)構(gòu)則避免了這一問題，骨架模型具有良好的連通性和細(xì)化性。

(2)通過對(duì)2種算法中心性的量化分析可得，細(xì)化法的平均骨架偏移距離為0.10 mm，比距離變換法(0.15 mm)減少了50%，說明細(xì)化法具有更強(qiáng)的魯棒性，更加適用于不同孔徑的土壤孔隙結(jié)構(gòu)。

(3)通過對(duì)骨架模型連通性、細(xì)化性和中心性3個(gè)指標(biāo)的定性和定量分析發(fā)現(xiàn)，細(xì)化法具有優(yōu)越的描述土壤孔隙形狀及拓?fù)涮卣鞯哪芰Γ蛇M(jìn)一步應(yīng)用于土壤物理結(jié)構(gòu)及水文特性的探索研究。