MATLAB數(shù)字圖像處理在根系固土中的應(yīng)用研究

閻夢(mèng)晴，于司杭，王培亮，王逸飛，李明寶

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

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

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

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

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

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

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

1.1 邊緣檢測(cè)算法 從輸入的圖像中提取屬性，進(jìn)行圖像分割。分割算法的基本原理是基于灰度的突變來(lái)分割圖像，基于檢測(cè)灰度值的不連續(xù)來(lái)檢測(cè)邊緣。灰度的不連續(xù)性的檢測(cè)采用一階導(dǎo)數(shù)和二階導(dǎo)數(shù)。圖像處理中的二階導(dǎo)數(shù)由于對(duì)噪聲過(guò)于敏感而很少直接用于邊緣檢測(cè)，僅在與其他邊緣檢測(cè)技術(shù)結(jié)合時(shí)使用，在此不作討論。

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

向量?f的幅值為

梯度的幅值通常簡(jiǎn)稱為梯度，由式(1)可知梯度具有微分的性質(zhì)，當(dāng)灰度恒定時(shí)梯度值為零，當(dāng)灰度變化時(shí)梯度值相應(yīng)變化，由此可通過(guò)計(jì)算梯度值來(lái)檢測(cè)灰度的變化程度，從而檢測(cè)邊緣。確定圖像中灰度快速變化的位置，基本準(zhǔn)則是尋求灰度的一階導(dǎo)數(shù)的幅值大于某個(gè)指定閾值的位置。在MATLAB中，采用edge函數(shù)進(jìn)行邊緣檢測(cè)。

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

1.2 算法比較

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

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

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

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

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

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

高斯拉普拉斯算子為

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

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

2 結(jié)果與分析

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

2.1 邊緣檢測(cè)器的選用 土壤的裂縫圖像來(lái)自于牛國(guó)權(quán)［8］所做的根系固土抗裂有效性試驗(yàn)，見(jiàn)圖2(a)。

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

圖2 不同邊緣檢測(cè)器的處理結(jié)果

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

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

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

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

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

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

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

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

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

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

圖3 使用Otsu方法進(jìn)行閾值處理后的結(jié)果

圖4 輸入圖像的直方圖

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

圖5 基于梯度的邊緣信息改進(jìn)全局閾值處理結(jié)果

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

圖6 最終結(jié)果圖

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

3 結(jié)論

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

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

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

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

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

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

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

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

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

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