程曉菲，武 剛

(北京林業(yè)大學(xué)信息學(xué)院，北京 100089)

單木是組成陸地生態(tài)系統(tǒng)的最小喬木實(shí)體。隨著新一代信息技術(shù)特別是遙感影像處理技術(shù)的發(fā)展與普及，高效采集海量單木信息建立單木數(shù)據(jù)庫(kù)，是實(shí)現(xiàn)精準(zhǔn)林業(yè)特別是城市林業(yè)精準(zhǔn)集約化管理的基礎(chǔ)，也是新一代智慧林業(yè)發(fā)展的目標(biāo)[1-4]。事實(shí)上，國(guó)內(nèi)一些個(gè)別的林業(yè)領(lǐng)域已實(shí)現(xiàn)單木集約管理，如古樹名木、北京市楊柳雌株精準(zhǔn)管理、果樹管理等[5-6]，但這些主要通過(guò)傳統(tǒng)的地面調(diào)查實(shí)現(xiàn)，效率低且工作量巨大?；谶b感影像的單木提取是實(shí)現(xiàn)單木數(shù)據(jù)庫(kù)高效建設(shè)的關(guān)鍵技術(shù)，其中基于遙感影像的單木定位(探測(cè))是基礎(chǔ)。

近年來(lái)遙感平臺(tái)和相應(yīng)的處理技術(shù)發(fā)展很快[7]，雖然主動(dòng)遙感技術(shù)(如激光雷達(dá)LiDAR等)已用于單木檢測(cè)(包括單木位置、樹高等)和反演模型建立等方面研究[8]，但其成本較高且不適用于大規(guī)模獲取。而基于高分辨率光學(xué)遙感影像的單木提取方法，因?yàn)槠渚哂屑夹g(shù)基礎(chǔ)深厚、影像尺度大且易得性高、便于推廣應(yīng)用等特點(diǎn)，依然是研究熱點(diǎn)[9-10]。大量關(guān)于被動(dòng)遙感的研究表明：遙感數(shù)據(jù)源[11]、季節(jié)及樹木生長(zhǎng)狀況[12]、定位方法選擇[13]等均為影響單木定位效果的重要因素。針對(duì)單木定位問題，學(xué)者們提出了以局部最大值法[9]為代表的探測(cè)單木樹尖的直接定位方法，以及單木樹冠探測(cè)與繪制(ITCD，Individual Tree Crown Detection and Delineation)等間接定位方法[14-15]，如各種改進(jìn)的分水嶺法[16-17]、谷底跟蹤法[18]、區(qū)域生長(zhǎng)法[19]和標(biāo)記點(diǎn)過(guò)程[3]。近年來(lái)學(xué)者們也嘗試各種定位新方法，如爬峰法[20]、Snake模型[21-22]、機(jī)器學(xué)習(xí)方法[23]等。但目前單木定位方法精度難以達(dá)到實(shí)用要求，還沒有一種方法能夠適應(yīng)各種林分[24-25]。對(duì)于ITCD而言，其主要問題是單木排列緊密及冠層交錯(cuò)均會(huì)導(dǎo)致單個(gè)樹冠內(nèi)、重疊樹冠間出現(xiàn)過(guò)分割或欠分割問題，導(dǎo)致單木定位產(chǎn)生誤判和漏判問題[26-27]，因此需要探究更高效且穩(wěn)定的單木定位新方法。

Chan等[28]提出了一種基于區(qū)域活動(dòng)輪廓模型的圖像分割方法，即Chan-Vese模型(簡(jiǎn)稱CV模型)。CV模型利用像元灰度區(qū)域的統(tǒng)計(jì)信息，能夠較好處理模糊邊界或者離散狀邊界的目標(biāo)[29]，已廣泛應(yīng)用于醫(yī)學(xué)、農(nóng)業(yè)等領(lǐng)域[30-33]。但是，CV模型應(yīng)用時(shí)，初始輪廓的大小、位置及形狀，均影響其分割的效率和結(jié)果[34]。實(shí)際應(yīng)用中，通常采用人工介入的方法將包裹目標(biāo)物體的圓或矩形作為初始輪廓，其效率低且在一定程度上忽略了圖片的視覺特征[35]。

因此，本研究基于樹木的形態(tài)學(xué)及其影像光譜特征，在自動(dòng)提取初始輪廓的基礎(chǔ)上，通過(guò)對(duì)多個(gè)實(shí)驗(yàn)區(qū)高分辨率衛(wèi)星影像進(jìn)行多方法比較實(shí)驗(yàn)分析，探討基于CV模型的單木定位方法應(yīng)用效果。

1 材料與方法

1.1 研究區(qū)域及其遙感影像

首先選擇了北京林業(yè)大學(xué)及其附近的2個(gè)實(shí)驗(yàn)區(qū)，通過(guò)目視和實(shí)地測(cè)量，進(jìn)行了初步試驗(yàn)，并與河北文安界圍農(nóng)場(chǎng)、深圳寶安荔枝園、山東棲霞市蘋果園、浙江臨安青山湖綠道的4個(gè)實(shí)驗(yàn)區(qū)實(shí)驗(yàn)影像進(jìn)行比較，具體的實(shí)驗(yàn)影像數(shù)據(jù)及其參數(shù)如表1所示。試驗(yàn)區(qū)1為北京市海淀區(qū)的東升八家郊野公園，主要樹種為針葉圓柏(Sabinachinensis)，該地區(qū)樹木密度較高，樹冠重疊現(xiàn)象較為嚴(yán)重。裁剪的實(shí)驗(yàn)影像1如圖1A中藍(lán)框區(qū)域所示，其中目視定位且實(shí)地驗(yàn)證共有397棵參考單木。試驗(yàn)區(qū)2為北京林業(yè)大學(xué)校園，裁剪的實(shí)驗(yàn)影像6(如圖1F中藍(lán)框區(qū)域)，該影像含有籃球場(chǎng)及其周邊，目視定位且實(shí)地驗(yàn)證了24棵參考單木，其樹冠茂盛，重疊嚴(yán)重，樹種主要為闊葉的國(guó)槐(Sophorajaponica)和饅頭柳(Salixmatsudanavar.matsudanaf.umbraculifera)。它們分別代表現(xiàn)實(shí)單木存在的林分和非林分(四旁樹)兩種形式。

后續(xù)為了對(duì)比考察CV模型單木定位法對(duì)不同類型影像源的適用性及其定位效果，選擇了試驗(yàn)區(qū)2的”北京二號(hào)”影像7(圖1G)，以及目視定位可靠且代表不同類型(針葉林、闊葉林、經(jīng)濟(jì)林)的4幅林分影像2～5(圖1B—1E)[19,36-39]。

圖1 研究區(qū)實(shí)驗(yàn)影像Fig.1 Experimental image of study area

1.2 圖像預(yù)處理

一般單木定位(探測(cè))的圖像預(yù)處理包括校正、配準(zhǔn)、增強(qiáng)、樹木區(qū)域提取等環(huán)節(jié)，它們是準(zhǔn)確提取單木位置信息的基礎(chǔ)。本研究遙感影像已進(jìn)行過(guò)校正和配準(zhǔn)，實(shí)驗(yàn)只涉及圖像增強(qiáng)和樹木區(qū)域提取2個(gè)環(huán)節(jié)。

為提高圖像的層次感，增強(qiáng)樹冠和非樹冠之間的差異，減少同一樹冠內(nèi)由于小樹枝及陰影引起的噪聲，采用高斯濾波方法增強(qiáng)圖像[14]。

樹木區(qū)域提取的目的是盡可能地剔除遙感影像中非樹木區(qū)，如土壤、街道、房屋、操場(chǎng)、草地等，以保證定位處理方法僅作用于樹木區(qū)域。通常在高分辨率遙感影像中，采用歸一化植被指數(shù)(NDVI)用于掩膜運(yùn)算。

采用綠度分割方法[24,40]：

dgreen=(R-G)/(R+G)。

(1)

其中：dgreen為綠度，R和G分別代表紅色波段和綠色波段光譜反射值?；诰G度統(tǒng)計(jì)直方圖，使用最大類間方差法(OTSU)確認(rèn)閾值[41]，去除小于閾值的像元點(diǎn)，以此將圖像中的道路、草地等無(wú)關(guān)信息剔除，獲取樹木區(qū)域圖像I。

1.3 基于CV模型的定位處理

1.3.1 初始輪廓提取

CV模型常用于獲取目標(biāo)物體的精確輪廓，首先沿著物體邊緣大致設(shè)置一條閉合曲線，使其在約束條件下，朝著內(nèi)法線的方向移動(dòng)，最終在目標(biāo)邊界停止。在傳統(tǒng)方法中，通常采用基于人工交互的方式將其設(shè)置在目標(biāo)邊緣。單木定位中，因樹冠數(shù)量諸多，手工設(shè)置初始輪廓工作繁重且不實(shí)用，宜采用自動(dòng)化方法實(shí)現(xiàn)初始輪廓提取。根據(jù)樹冠頂點(diǎn)受太陽(yáng)光照強(qiáng)度較大，在影像中通常具有較高的光譜反射值，而樹冠邊緣反射值較低[42]，本研究通過(guò)計(jì)算影像中多個(gè)局部最大值區(qū)域并進(jìn)行連通域標(biāo)記，作為初始輪廓。其中，局部最大值區(qū)域是指連通且值相同的像元區(qū)，其周圍像元值均小于區(qū)域的像元值。具體提取過(guò)程如下：

1)設(shè)樹木區(qū)域圖像I的像元(x,y)灰度值為G(x,y)，其8鄰域內(nèi)像元灰度值為Gneigh(x,y)，遍歷每個(gè)像元，則由公式(2)可得到一個(gè)二值圖A。

(2)

2)對(duì)二值圖A中的局部最大值區(qū)域采用4鄰域方法進(jìn)行形態(tài)學(xué)連通域查找并標(biāo)記[43]，獲取初始輪廓標(biāo)記圖C=C1∪C2∪…∪CL。其中，Ci為第i個(gè)連通域，Ci(x,y)=i，若記Ci在樹木區(qū)域圖像I中對(duì)應(yīng)的區(qū)域?yàn)棣竔,i=1,2,…,L，則Ω=Ω1∪Ω2∪…∪ΩL，即為CV模型演化的初始輪廓。

1.3.2 樹冠輪廓分割

基于CV模型對(duì)初始輪廓進(jìn)行演化迭代，主要是通過(guò)輪廓曲線的能量函數(shù)值的不斷優(yōu)化完成。樹冠輪廓曲線由水平集函數(shù)φ(x,y,t)表達(dá)，曲線把圖像劃分為內(nèi)外兩個(gè)部分，c1和c2分別記為曲線內(nèi)部和外部的圖像像元平均值，將圖像的能量函數(shù)定義為：

(3)

其中：μ≥0、v≥0、λ1≥0、λ2≥0是各個(gè)能量項(xiàng)的權(quán)重系數(shù)；x、y分別為圖像橫縱軸坐標(biāo)；當(dāng)曲線逼近樹冠邊緣時(shí)取值接近0，且c1=c2，此時(shí)曲線ci位于目標(biāo)樹冠邊界。H(φ)和δ(φ)分別為Heaviside函數(shù)和Dirac函數(shù)的正則化表示，用于劃分演化區(qū)域，限定輪廓在演變過(guò)程中在零水平集函數(shù)附近取值，其表達(dá)式為：

(4)

(5)

式中，ε為數(shù)字運(yùn)算選取的一個(gè)小正數(shù)(閾值)。

依據(jù)變分原理和梯度下降原理，對(duì)水平集函數(shù)φ求式(3)能量泛函極小化，得到其偏微分方程：

(6)

式中：v為能量項(xiàng)的權(quán)重系數(shù)，φ0(x,y)為初始水平集函數(shù)。

式(6)可通過(guò)有限差分法求解，即設(shè)空間步長(zhǎng)為h，時(shí)間步長(zhǎng)為Δt，則得到迭代方程為[28]：

(7)

樹冠輪廓具體迭代計(jì)算步驟為：

①n=0，根據(jù)待分割圖Ω，建立初始水平集函數(shù)φ0(x,y)，并根據(jù)式(3)構(gòu)造能量函數(shù)；

②根據(jù)當(dāng)前樹冠輪廓迭代結(jié)果φn，按照式(7)計(jì)算φn+1；

③若|φn+1-φn|≤ε(閾值)，則N=n，轉(zhuǎn)步驟④，否則n=n+1,轉(zhuǎn)向步驟②；

④輸出最終輪廓(樹冠曲線)φN(x,y)。

1.3.3 單木位置計(jì)算

對(duì)上述分割樹冠輪廓二值圖φN(x,y)進(jìn)行連通域查找，基于圖像幾何矩原理，用公式(8)計(jì)算其重心作為單木位置坐標(biāo)(x,y)[44]，連通域個(gè)數(shù)即為提取的單木數(shù)量。

(8)

其中：φN(x,y)為二值圖的像元值；C10和C01為圖像一階矩，分別表示樹冠輪廓x坐標(biāo)值與像元灰度值乘積的累加和以及樹冠輪廓y坐標(biāo)值與灰度值乘積的累加和；C00為圖像零階矩，表示樹冠輪廓像元灰度值的累加和。

1.4 定位結(jié)果評(píng)價(jià)指標(biāo)

影像1的局部影像機(jī)器定位與目視定位的對(duì)比結(jié)果見圖2，其中圖2B為CV模型定位結(jié)果與目視定位的對(duì)比，圖2C為標(biāo)記分水嶺定位結(jié)果與目視定位的對(duì)比。

圖2 影像機(jī)器定位與目視定位對(duì)比Fig.2 Comparison between image machine positioning and visual positioning

2 結(jié)果與分析

以配有Intel Core(TM)i5、2.60GHz、12GB內(nèi)存、64位Windows 10操作系統(tǒng)，利用PyCharm Community Edition 2018.3.5、Python編程語(yǔ)言及其OpenCV、MATLAB. Engine等組件開發(fā)的算法實(shí)驗(yàn)平臺(tái)，對(duì)7幅影像以目視解譯數(shù)據(jù)作為參考，采用CV模型和3種傳統(tǒng)的主流單木定位方法(梯度分水嶺法、標(biāo)記分水嶺法和局部最大值法)進(jìn)行了實(shí)驗(yàn)比較分析。

2.1 CV模型的單木位置圖像精度

影像2的各算法實(shí)驗(yàn)影像見圖3。其中圖3A—3D分別為圖像原圖、初始輪廓圖、CV模型樹冠輪廓圖以及基于CV模型的單木位置圖像。圖3E—3G為局部最大值法、梯度分水嶺法以及標(biāo)記分水嶺的單木定位圖像?？梢钥闯?，基于樹冠的光譜反射原理，初始輪廓大概標(biāo)記了單木位置，并在此基礎(chǔ)上，進(jìn)行水平集函數(shù)的構(gòu)建及迭代，使樹冠輪廓逐步收斂到目標(biāo)邊緣，能夠較好地識(shí)別單木。相較于局部最大值的漏判，梯度分水嶺多棵單木誤識(shí)別為1棵，CV模型能夠較好地進(jìn)行單木定位。7組實(shí)驗(yàn)的單木位置精度評(píng)價(jià)見表2，CV模型均取得了最高匹配率。

在影像1排列緊密、大小均勻且排列較為規(guī)則的針葉林中，因其樹冠輻射值較高且分布較均勻，局部最大值在3種傳統(tǒng)方法中匹配率(M)最高，梯度分水嶺及標(biāo)記分水嶺等傳統(tǒng)方法將多個(gè)小樹冠誤合并識(shí)別為1個(gè)樹冠，因此漏判率(Rom)較高；影像2～3中，樹冠大小不一且排列凌亂，局部最大值方法無(wú)法準(zhǔn)確獲取單木位置，導(dǎo)致了一定程度的漏判和誤判，樹冠間無(wú)明顯區(qū)分，梯度分水嶺及標(biāo)記分水嶺法易將樹冠誤分割為多部分， CV模型的準(zhǔn)確率(Rmat)較高，并且誤判率(Rcom)較低，能夠較好地識(shí)別較小單木，處理漏判現(xiàn)象，具有一定的泛化效果。在影像4～7中，闊葉樹樹冠較大且樹冠重疊嚴(yán)重，傳統(tǒng)方法將樹冠誤分割為多個(gè)樹冠，出現(xiàn)了一定的過(guò)分割現(xiàn)象，也因此導(dǎo)致了較高的誤判率，CV模型較好地緩解了這一狀況，但存在一部分相鄰樹冠難以分割的情況，而總體匹配率效果較好。

圖3 影像2各算法實(shí)驗(yàn)影像圖Fig.3 Image 2 experimental results of each algorithms

表2 單木位置精度評(píng)價(jià)

2.2 CV模型影像定位效果

所有實(shí)驗(yàn)影像的整體定位效果見表3，局部最大值方法因其受移動(dòng)窗口因素影響較大，同一樹冠(特別是闊葉樹)中易產(chǎn)生多個(gè)局部最大值點(diǎn)，因此導(dǎo)致較高的誤判率(Rcom)，平均為0.434。梯度分水嶺和標(biāo)記分水嶺通過(guò)在原始圖像上分別計(jì)算圖像梯度函數(shù)及標(biāo)記樹冠大概位置，在其基礎(chǔ)上進(jìn)行分水嶺分割，容易產(chǎn)生過(guò)分割現(xiàn)象，誤判率較高。相較于其他3種定位方法，CV模型在7組影像實(shí)驗(yàn)中均取得了最高的匹配率(M)，平均高出其他方法約23%，綜合定位效果最好。在同一地區(qū)不同影像中(影像6～7)，不同類型樹木中均具有較好的定位效果。但CV模型在樹木排列緊密區(qū)域，目視不清時(shí)，也容易存在一定的欠分割問題，導(dǎo)致漏判率較高(影像3～7)。影像6的局部區(qū)域?qū)嶒?yàn)圖像見圖4。該區(qū)域樹冠緊密，CV模型將排列緊密樹冠誤認(rèn)為同一個(gè)(圖4B)，僅識(shí)別出1棵單木(圖4C)。

表3 整體單木定位精度

圖4 CV模型局部實(shí)驗(yàn)圖像Fig.4 Local experimental image of CV model

3 結(jié) 論

本研究提出的CV模型單木定位法，首先利用樹木的形態(tài)學(xué)特征及其在影像中的光譜特征，自動(dòng)設(shè)置初始輪廓；其次，進(jìn)行基于水平集函數(shù)構(gòu)建及迭代，對(duì)樹冠輪廓進(jìn)行演變；最終獲取單木位置。該方法基于圖像全局信息，利用曲線內(nèi)外的灰度均值而不是梯度信息進(jìn)行分割，能夠在邊界模糊或梯度無(wú)意義的圖像中取得較好的分割效果，能夠快速準(zhǔn)確地收斂到目標(biāo)位置。實(shí)驗(yàn)表明：該方法比其他3種方法(梯度分水嶺法、標(biāo)記分水嶺法、局部最大值法)具有更高的匹配率，平均高出23%，且與影像源無(wú)關(guān)，顯示了良好的泛化應(yīng)用潛力。

更全面地評(píng)價(jià)單木信息還需要進(jìn)行更多方法(如區(qū)域增長(zhǎng)、模板匹配、機(jī)器學(xué)習(xí)等)和更多影像的對(duì)比實(shí)驗(yàn)。同時(shí)，實(shí)驗(yàn)也表明，雖然CV模型的單木定位法一定程度上解決了單木樹冠識(shí)別的過(guò)分割及欠分割問題，綜合定位效果有所提高，但在郁閉度較高、目視都難以識(shí)別單木的情況下，仍存在一定的欠分割問題；后續(xù)可針對(duì)郁閉度差異進(jìn)行單木定位實(shí)驗(yàn)研究，探討其適用性。隨著機(jī)器學(xué)習(xí)方法在目標(biāo)檢測(cè)方面的突破，新的單木定位方法值得期待。另外，結(jié)合多時(shí)像(春、夏、秋、冬)高分辨率影像的融合應(yīng)用技術(shù)的研究，或許對(duì)解決郁閉度高、單木難定位識(shí)別問題有所幫助。