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

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

0 引言

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

近年來，計算機斷層掃描技術(Computer tomography, CT)在土壤領域的成功應用，能夠實現土壤孔隙三維結構的可視化[4-5]。由于土壤組成成分的復雜性，使得孔隙結構具有多種幾何形狀和多尺度特征：①孔隙呈不規(guī)則的幾何形狀，并向多個方向發(fā)育延伸。②具有不同數量級的幾何尺寸。③多尺度特征的孔隙相互交錯。復雜的拓撲結構，使得土壤孔隙的研究難度大。為了探索復雜孔隙結構的體積、曲折度、連通性等特征，需要結合重構技術對其空間結構和拓撲特性進行分析，在這一過程中，孔隙骨架模型的構建方法至關重要。

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

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

1 土壤樣本采集與圖像預處理

1.1 土壤樣本的采集

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

1.2 土壤CT圖像的獲取

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

1.3 土壤CT圖像的預處理

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

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

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

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

2 研究方法

2.1 距離變換法

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

(1)

(2)

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

B——所有邊界體素點的集合

v1、v2——體素點

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

2.2 細化法

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

2.2.1細化條件的建立

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

(1)為了保證骨架模型的連通性，所有刪除的目標點都應保證在其26鄰域內至少有2個相同的體素點。按照該條件進行模型遍歷，直到生成的骨架模型中目標體素點的數目不再改變。

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

2.2.2偽分支和毛刺的處理

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

去偽分支和毛刺處理主要包括毛刺的刪除和狹窄缺口的填補，其實現過程主要包括3個步驟：①依據某一準則判斷出組成孔隙骨架中軸線的體素點。②以這些體素點組成的中心路徑為標準，刪除路徑之外的冗余體素點，以去除骨架上的毛刺。③通過形態(tài)學操作，填補骨架路徑的狹窄缺口，從而保證孔隙骨架模型的拓撲不變性。

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

3 實例應用與驗證

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

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

3.1 規(guī)則結構的骨架模型

為詳細比較2種算法對規(guī)則結構骨架模型的構建效果，選取了規(guī)則曲線、規(guī)則圓環(huán)和規(guī)則圓柱體進行骨架模型的構建。如圖2所示，灰色區(qū)域為原始結構，紅色曲線為提取的骨架模型。

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

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

3.2 土壤孔隙骨架模型

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

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

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

3.3 骨架模型中心性評價

為了更精確評價2種骨架化算法的性能，本文采用中心性作為定量評價標準。由于該指標均需要基于標準骨架結構進行分析，而土壤孔隙的真實骨架模型難以確定，因此，試驗選用CT掃描的規(guī)則圓柱體模型進行量化試驗的測試對象。

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

(3)

式中K——骨架的數量

i——骨架結構的序號

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

圖4 圓柱體模型中心性評價結果Fig.4 Evaluation results of centrality

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

結合規(guī)則幾何結構和土壤孔隙結構的試驗結果可知，細化法構建的骨架模型具有良好的細化性、連通性和中心性，能夠更好地保留土壤孔隙結構的拓撲特性，可為土壤孔隙結構的空間分布研究提供一種先進的技術手段。

4 結論

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

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

(3)通過對骨架模型連通性、細化性和中心性3個指標的定性和定量分析發(fā)現，細化法具有優(yōu)越的描述土壤孔隙形狀及拓撲特征的能力，可進一步應用于土壤物理結構及水文特性的探索研究。