閻夢晴，于司杭，王培亮，王逸飛，李明寶

(東北林業(yè)大學土木工程學院，黑龍江哈爾濱150040)

我國是世界上水土流失最嚴重的國家之一［1］。水土流失會破壞土壤的肥力，造成土壤干層，淤積水庫、阻塞河道、抬高河床，惡化生態(tài)環(huán)境。雖然自新中國成立以來，我國的水土流失治理已經(jīng)逐漸成為國建的生態(tài)建設重點工程，但是由于我國國土面積大，自然地理條件特殊，以及長期以來對水土資源的過度利用，我國水土流失仍然面積大、分布廣，防治任務艱巨［2］。防止水土流失，采用科學、環(huán)保的植物固土措施，引入數(shù)字圖像處理技術研究植物根系固土機制，能夠為我國農(nóng)業(yè)水土保持增添新的技術手段，促進我國生態(tài)文明建設的發(fā)展。

水土流失的防治措施有工程措施和植物措施，其中植物措施不僅可以防治水土流失，還可改善局部的生態(tài)環(huán)境，節(jié)約經(jīng)濟成本，具有生態(tài)可持續(xù)性。近年來，由于生態(tài)環(huán)境建設的需要及各學科科學技術的發(fā)展，植物根系固定土體機理的研究成為根系研究的焦點［3］。根系與土壤形成復合結構體，起到固定植物地上部分以及固定土壤、防止水土流失的作用。對植物根系固定土體的機理進行綜合定量研究，對于搭建優(yōu)良的人工水土保持生態(tài)系統(tǒng)不僅具有廣闊的應用前景，而且具有重要的學術意義［1］。對于根系固土力學機理的研究，已經(jīng)得到越來越多學者的重視，植物根系防治水土流失已經(jīng)成為水土保持領域的一個研究熱點。

近年來，國內外科研專家運用各種科學技術方法對根系固定土體的力學機理進行了研究，測定植物根系本身的抗拉力［4－5］，采用原位剪切試驗得到的有根土和無根土的抗擠強度值作為衡量根系固土性質的指標［6］。隨著科學技術手段的不斷進步，土壤中根系的幾何形態(tài)和分布情況能夠通過各種試驗方法得到，在研究植物根系在土壤中的分布以及根系幾何形態(tài)的提取與仿真方面，也不可避免地需要對根土斷層圖像、土壤裂縫圖像進行數(shù)字圖像處理，為科研人員提供根系的形態(tài)參數(shù)數(shù)據(jù)。

運用數(shù)字圖象處理技術，對根系固土抗裂性試驗中的試驗圖像進行特征提取，將裂縫的圖像特征(如裂縫寬度、深度等)用以衡量根系固土的抗裂有效性，補充完善了宏觀的試驗指標，通過這些指標可以直觀地反映根系固土的程度。通過對圖像特征的提取，建立更加全面的參數(shù)，這些參數(shù)可作為研究抗裂性及其影響因素之間的橋梁。運用數(shù)字圖像處理，能夠使根系固土抗裂性這一性質轉化為更加具體且全面的數(shù)學指標，通過這些量化的指標可建立抗裂性與其影響因素之間的關系，探究它們的影響機理，這將為研究植物根系固土防治水土流失奠定很好的理論基礎。數(shù)字圖像處理技術的引入，將成為水土保持研究領域的一個重要環(huán)節(jié)，將為分析根系固土機制提供有效的技術手段。

1 基于MATLAB的數(shù)字圖像處理

數(shù)字圖像處理，把一幅圖像定義為一個二維函數(shù)f(x，y)，其中x、y是空間坐標，f是任意坐標(x，y)處的幅度(稱為亮度或灰度)。MATLAB(矩陣實驗室)是處理矩陣和線性代數(shù)的強大的數(shù)值計算軟件，將圖像作為矩陣來處理，選用MATLAB軟件能有效地提高圖像處理的效率，準確地提取出數(shù)字圖像的特征。將坐標值(x，y)數(shù)字化稱為采樣，將幅值8數(shù)字化稱為量化，采樣和量化后得到一個實數(shù)矩陣。假設，對圖像f(x，y)采樣后得到M×N的圖像，坐標值為離散量。將基于MATLAB的數(shù)字圖像處理技術應用于根系固土抗裂性試驗，旨在得到土壤裂縫的圖像特征。

1.1 邊緣檢測算法 從輸入的圖像中提取屬性，進行圖像分割。分割算法的基本原理是基于灰度的突變來分割圖像，基于檢測灰度值的不連續(xù)來檢測邊緣?；叶鹊牟贿B續(xù)性的檢測采用一階導數(shù)和二階導數(shù)。圖像處理中的二階導數(shù)由于對噪聲過于敏感而很少直接用于邊緣檢測，僅在與其他邊緣檢測技術結合時使用，在此不作討論。

圖像二維函數(shù)f(x，y)的梯度為

向量?f的幅值為

梯度的幅值通常簡稱為梯度，由式(1)可知梯度具有微分的性質，當灰度恒定時梯度值為零，當灰度變化時梯度值相應變化，由此可通過計算梯度值來檢測灰度的變化程度，從而檢測邊緣。確定圖像中灰度快速變化的位置，基本準則是尋求灰度的一階導數(shù)的幅值大于某個指定閾值的位置。在MATLAB中，采用edge函數(shù)進行邊緣檢測。

在MATLAB中，根據(jù)土壤裂縫的圖像以及邊緣檢測器的功能，選用 Sobel、Prewitt、Roberts、LoG、Canny 這 5 種邊緣檢測器，對比它們邊緣檢測的效果，以確定最適用于土壤裂縫圖像做分割的邊緣檢測器。

1.2 算法比較

1.2.1 Sobel邊緣檢測器。Sobel邊緣檢測器使用圖1(b)中的Sobel近似導數(shù)發(fā)現(xiàn)邊緣，即圖1(a)所示的3×3鄰域的行、列間的離散差來計算梯度，每行每列的中心像素以2加權提供平滑效果［7］。

圖1 邊緣檢測器模板及其實現(xiàn)的一階導數(shù)

式中，z為灰度。在(x，y)處，當?f≥T(T為指定閾值)，則(x，y)處的像素為邊緣像素。

1.2.2 Prewitt邊緣檢測器。Prewitt邊緣檢測器使用圖1(c)中的近似導數(shù)發(fā)現(xiàn)邊緣。與Sobel檢測器相比，Prewitt邊緣檢測器的計算更簡單，但是檢測結果中的噪聲較大。

1.2.3 Roberts邊緣檢測器。Roberts邊緣檢測器使用圖1(d)的近似導數(shù)發(fā)現(xiàn)邊緣。Roberts邊緣檢測器是最古老且最簡單的數(shù)字圖像處理邊緣檢測器之一，優(yōu)點在其速度快，但功能不全面。

1.2.4 LoG 檢測器。LoG(Laplacian of a Gaussian)檢測器使用了高斯拉普拉斯算子與圖像卷積(濾波)。對于高斯函數(shù)

高斯拉普拉斯算子為

使用?2G(x，y)卷積圖像，使得圖像平滑，降低了噪音，計算拉普拉斯算子得到雙邊緣圖像，查找雙邊緣之間的零交叉即可定位邊緣。

1.2.5 Canny邊緣檢測器。與LoG檢測器相比，Canny邊緣檢測器指定了高斯濾波器中的標準差σ，以此來平滑圖像，減少噪聲。在(x，y)處計算局部梯度［+］1/2和邊緣方向arctan(gx/gy)，通過尋找f(x，y)的梯度局部最大發(fā)現(xiàn)邊緣。Canny邊緣檢測器使用了兩個閾值來檢測強邊緣和弱邊緣，此方法更有可能檢測到真正的弱邊緣。

2 結果與分析

為得到土壤裂縫的形態(tài)特征參數(shù)，要選用適宜的邊緣檢測器以及閾值處理方法。通過5種不同邊緣檢測器的處理結果圖與原裂縫圖的對比，選擇出精度最高的邊緣檢測器。在此基礎上，采用兩種不同的閾值處理方法，對比所得圖像的精確度，確定最適宜的閾值處理方法。

2.1 邊緣檢測器的選用 土壤的裂縫圖像來自于牛國權［8］所做的根系固土抗裂有效性試驗，見圖2(a)。

圖2為5種不同邊緣檢測器處理土壤裂縫圖像后的結果圖。由圖2(b)、(c)可知，Prewitt邊緣檢測器和Sobel邊緣檢測器的效果均不佳，二者顯然無法檢測到圖像中的弱邊緣，無法有效分割出土壤裂縫的圖像。Roberts邊緣檢測器原理如圖1所示，與Prewitt和Sobel邊緣檢測器相比，Roberts邊緣檢測器模板是非對稱的，雖然簡單，但用于檢測土壤裂縫的圖像特征時，如圖2(d)所示，其效果不好。由圖2(e)可知，LoG檢測器結果比前3種檢測器的效果好，但不及Canny邊緣檢測器。由圖2(f)可知，Canny邊緣檢測器功能強大，能檢測到細微的圖像特征，但是在檢測土壤裂縫時，還需調整參數(shù)盡可能地減少不相關的細節(jié)。

圖2 不同邊緣檢測器的處理結果

2.2 圖像閾值處理 做分割后得到的裂縫圖像是初步處理結果。在初步處理圖像時，選擇的是各個檢測器自動選擇的閾值，得到的結果不精確。需根據(jù)裂縫圖像的圖形特征，選擇最優(yōu)閾值。

閾值處理后的圖像為二值圖像，它定義為

式中，a像素對應物體，b像素對應背景，按照慣例，a=1(白)，b=0(黑)。

2.2.1 使用Otsu方法進行最佳全局閾值處理。將圖像的直方圖的成分表示為

其中，n為圖像像素總和，nq為具有灰度級q的像素數(shù)量，L為圖像中可能的灰度級總數(shù)(灰度級為整數(shù))。

假設初始閾值為 k，C1 是灰度級為［0，1，2，…，k］的一組像素，C2是灰度級為［k+1，k+2，…，L －1］的一組像素，則最大類間方差為

其中，P1(k)為集合C1發(fā)生的概率，P2(k)為集合C2發(fā)生的概率，mG為全局均值(平均灰度)。

令P2(k)=1－P1(k)，則式(8)可寫為

通過類間方差最大化思想即可確定最佳閾值，在MATLAB中，采用工具箱函數(shù)graythresh()可計算Otsu閾值T=0.431 4。據(jù)此閾值再采用Canny邊緣檢測器即可得到如圖3所示圖像。

2.2.2 使用基于梯度的邊緣信息改進全局閾值處理。在MATLAB中，將土壤裂縫圖像作為輸入圖像，采用imhist()函數(shù)，用plot方法繪制輸入圖像的直方圖。如圖4所示，該圖像直方圖是單峰的，通常，當物體遠小于背景時，物體對直方圖的貢獻可忽略不計，為改善這種情況，往往使用邊緣信息以得到梯度圖像。

圖3 使用Otsu方法進行閾值處理后的結果

圖4 輸入圖像的直方圖

自定義函數(shù)percentile2i()，用以設定在邊緣圖像中將較少的像素用于計算閾值，即使用高百分位(99.9)來估計梯度的閾值，因為梯度圖像中較大的值是出現(xiàn)在物體與背景的邊界處。處理結果如圖5所示，閾值T=0.368 6。

圖5 基于梯度的邊緣信息改進全局閾值處理結果

比較圖3和圖5，圖像閾值處理的結果顯然前者較好，這說明在處理土壤裂縫這類圖片時，使用Otsu方法進行閾值處理能得到最優(yōu)閾值。結合Canny邊緣檢測器和Otsu閾值處理方法，最后得到的裂縫圖像如圖6所示。

圖6 最終結果圖

土壤裂縫圖像不同于其他圖像，為了獲取裂縫的形態(tài)參數(shù)，通過對比研究了5種邊緣檢測器以及2種閾值處理方法，確定了最優(yōu)邊緣檢測器和閾值處理方法，即結合Canny邊緣檢測器和Otsu閾值處理的方法，能夠得到最精確的裂縫圖像。圖6所示的圖像能夠精確地反映原圖中土壤的裂縫形態(tài)。將二維圖像數(shù)字化，可得到裂縫邊緣的二值化圖像，在MATLAB中以矩陣形式存儲了圖像像素，再后期運用裂縫形態(tài)特征參數(shù)進行根系固土機制研究時，即可直接調用。

3 結論

該研究以根系固土試驗中土壤裂縫圖像為例，運用MATLAB軟件，對圖像進行數(shù)字圖像處理，對比Sobel、Prewitt、Roberts、LoG、Canny這5 種邊緣檢測器的效果，選擇最適合用于檢測土壤裂縫的Canny邊緣檢測器。比較基于梯度的邊緣信息改進全局閾值處理方法和Otsu方法進行最佳全局閾值處理的效果，確定使用Otsu方法進行閾值處理能得到最優(yōu)閾值。結合Canny邊緣檢測器和Otsu閾值處理方法，是最適用于土壤裂縫圖像的數(shù)字圖像處理方法。

為防止水土流失，對植物根系固定土體機制的研究中，引入基于MATLAB的數(shù)字圖像處理方法，能夠得到宏觀圖像中如根系、土壤裂縫等圖像的形態(tài)特征參數(shù)。處理后的根系、裂縫等邊緣像素的特征以及根系圖像像素尺寸與實際尺寸存在線性關系，根據(jù)各個參數(shù)的幾何特性可分析出根系的長度、表面積、平均直徑、體積以及根系間的夾角，或是裂縫的長度、寬度、面積等。通過該研究結果，可以為關于根系固土相關領域的研究人員提供可靠分析手段以及根系的形態(tài)參數(shù)數(shù)據(jù)。

［1］劉震.水土保持60年:成就·經(jīng)驗·發(fā)展對策［J］.中國水土保持科學，2009(4):1 －6.

［2］田衛(wèi)堂，胡維銀，李軍，等.我國水土流失現(xiàn)狀和防治對策分析［J］.水土保持研究，2008，15(4):204 －209.

［3］陳麗華，余新曉，宋維峰，等.林木根系固土力學機制［M］.北京:科學出版社，2008:87－115.

［4］張瑜.微根窗根系的圖像處理方法研究［D］.哈爾濱:東北林業(yè)大學，2012.

［5］史敏華，王棣，李任敏.石灰?guī)r區(qū)主要水保灌木根系分布特征與根抗拉力研究初報［J］.山西林業(yè)科技，1994(1):17－19.

［6］趙麗兵，張寶貴，蘇志珠.草本植物根系增強土壤抗剪切強度的量化研究［J］.中國生態(tài)農(nóng)業(yè)學報，2008，16(3):718 －722.

［7］GONZALEZ R C，WOODS R E.Digital Image Processing［M］.3rd ed.Prentice Hall，2007.

［8］牛國權.基于力學特性的根系固土抗裂性有效性研究［D］.呼和浩特:內蒙古農(nóng)業(yè)大學，2010.