孫 華，鞠洪波，張懷清，凌成星

（1. 中南林業(yè)科技大學 林業(yè)遙感信息工程研究中心，湖南 長沙 410004；2. 中國林業(yè)科學研究院資源信息所，北京 100091）

基于Worldview-2影像的林木冠幅提取與樹高反演

孫 華1,2，鞠洪波2，張懷清2，凌成星2

（1. 中南林業(yè)科技大學 林業(yè)遙感信息工程研究中心，湖南 長沙 410004；2. 中國林業(yè)科學研究院資源信息所，北京 100091）

以湖南省攸縣黃豐橋國有林場杉木人工林為例，探討林木冠幅提取與樹高反演方法研究?；赪orldview-2影像，采用均值漂移分割算法開展樣地內杉木冠幅信息提取。通過設置不同分割尺度確定最佳的冠幅分割參數為hs=10，hr=6，M=20。對提取的冠幅邊界進行平滑處理，利用平滑后的影像冠幅與實測樹高，分別建立了冠幅樹高曲線估計模型和非線性聯立方程組反演模型。其中以樹高作為啞變量，建立的影像冠幅樹高非線性聯立方程組模型的擬合效果最佳，模型決定系數R2為0.899，模型的變動系數（CV），平均百分標準誤差（EMPSE）均在10%以內，是樹高反演的一種有效手段。

林業(yè)遙感；林木冠幅；杉木人工林；圖像分割；均值漂移算法；非線性聯立方程組

進入21世紀，隨著高空間分辨率遙感影像不斷涌現，遙感技術在林業(yè)中的應用由最初的森林類型識別和信息提取，向更加精細的方向發(fā)展，已經開始用于森林參數提取研究。樹高是一個重要的森林參數，是森林蓄積量和地上生物量估測的一個重要指標。對于多光譜影像而言，通過提取影像各波段反射率、植被指數、紋理因子及其衍生變量的散點圖，采用空間統(tǒng)計學或逐步回歸技術提取與樹高相關性高的波段或變量，建立回歸模型來反演樹高是比較常見的方法[1-4]。近年來，冠幅胸徑模型以及冠幅樹高模型的研究受到了國內外的關注，已有的文獻研究表明冠幅與樹高、胸徑有著良好的數學關系，林木位置和樹冠信息提取成為高分辨率遙感影像研究的焦點[5-6]。為了從影像上快速識別和分析冠幅特征，并將冠幅信息從遙感影像數據中分離并提取出來，面向對象的分割技術被認為是一種有效的方式[7-9]。近年來隨著商業(yè)遙感分析軟件eCognition的應用，面向對象的遙感影像分析方法越來越受到研究者的關注。國內外學者在圖像分割方面開展了大量的研究，借助各種理論提出了上千種分割算法?，F已提出的分割算法都是針對具體的問題，并沒有一種適用于所有圖像的通用分割算法。對于樹冠信息提取而言，研究的算法主要側重于局部最大值法，基于輪廓的方法，模板匹配法，區(qū)域增長法，分水嶺算法等[10-13]。

論文以Worldview-2影像為數據源，采用均值漂移算法（Mean Shift algorithm）開展杉木樹冠信息提取與樹高反演研究。首先從影像中提取冠幅的邊緣特征，通過圖像分割的方式確定最佳的冠幅提取尺度。同時開展地面樣地調查，進行每木檢尺，記錄每棵樹的位置、胸徑、樹高、冠幅長度（東西和南北）、以及林分郁閉度等信息。分析調查數據中實測冠幅與影像提取冠幅之間的相關性。結合實測樹高與冠幅的關系，最終建立杉木人工林影像冠幅樹高聯立方程組反演模型。

1 材料與方法

1.1 研究區(qū)概況

湖南省攸縣黃豐橋林場介于113°04′～113°43′E，27°06′～ 27°24′N 之間。東北部與江西的蓮花、萍鄉(xiāng)交界，東南與茶陵縣接壤，西北部與株洲、醴陵毗鄰。全場地貌以中低山為主，境內最高海拔1 270 m，最低海拔115 m，坡度介于20～35°之間。林場地處中亞熱帶季風濕潤氣候區(qū)，年均氣溫17.8℃，平均無霜期為292 d，年均降水量1 410.8 mm。全場現有林地面積10 122.6 hm2，活立木蓄積879 705 m3，其中有林地蓄積879 688 m3，占活立木蓄積99.99%，四旁樹蓄積17 m3，僅占0.01%。林場的森林覆蓋率為86.24%。樹種主要以杉木Cunninghamia lanceolata(Lamb.) Hook、松類Pinusspp為主，其中杉木面積3 197.6 hm2，占用材林面積89.9%，蓄積593 738 m3，占96.56%，全部為人工造林[14]。

1.2 均值漂移算法

均值漂移算法（Mean Shift）是一種是非參數估計密度函數的方法，使每一個點通過有效的統(tǒng)計迭代“漂移”到密度函數的局部極大值點。算法思想最初由Fukunaga（1975）等人[15]在一篇關于概率密度梯度函數的估計論文中提出來的，Cheng（1995）對算法進行了推廣，定義了核函數以及設定了一個權重系數，擴大了算法的適用范圍。Comaniciu（2002）等人[16]證明了算法的收斂性，將Mean Shift算法運用到圖像的特征空間分析中。Mean Shift算法本質上是一個自適應的梯度上升搜索峰值的方法，近年來，Mean shift算法在圖像平滑、分割以及視頻跟蹤等領域獲得了非常成功的應用效果。并逐步應用于中高分辨率遙感影像的實現遙感影像的多尺度分割，在分割參數控制及穩(wěn)定性方面研究成果頗多，有關Meanshift算法的相關數學推導見參考文獻[16]。

1.3 建模方法

本研究擬采取曲線估計和近代統(tǒng)計學中聯立方程組法構建Worldview-2影像冠幅樹高模型[17]。為了表達方便，用HT代表實測樹高，CW代表實測冠幅，PCW代表提取的影像冠幅。聯立方程組模型構建具體步驟如下：

（1）通過地面實測樣木數據，構建HT-CW模型；

（2）建立基于分割算法提取的林木冠幅與實測冠幅（CW-PCW）模型；

（3）HT-CW模型與CW-PCW模型聯立，建立聯立方程組模型。

1.4 模型評價

利用確定系數R2、估計值的標準差（Standard Error of Estimate，SEE）、變動系數（Coeff i cient of Variation，CV）和平均百分標準誤差（Mean Percent Standard Error，MPSE）4個評價指標對模型進行檢驗。

式中，yi表示實測值，y?i表示預測值，y表示實測值的均值，p為模型參數個數。

2 結果與分析

2.1 冠幅提取分析

樹木的冠幅信息在圖像上所表示同質范圍相對較小，要想在遙感影像上準確提取冠幅信息，分割尺度是一個必須要考慮的問題。對于Mean Shift而言，通過控制核的帶寬參數（hs，hr），有效的將同質區(qū)域影像分割到同一對象塊中，最后利用最小區(qū)域面積（M）參數逐步實現影像的多尺度分割?？紤]到單棵樹的樹冠面積相對較小，在保持空間特征帶寬hs與顏色特征帶寬hr不變（hs=10，hr=6）的情況下，逐步減少最小區(qū)域面積參數M的值，發(fā)現隨著最小區(qū)域面積值的減少，分割出來的斑塊數量呈幾何級數增長（圖1），而且分割效果不佳，出現了過分割現象。

圖1 斑塊數量與最小區(qū)域面積尺度M的相關曲線Fig.1 Relation between plaque numbers and minimum area scale parameter

圖2 （hr=6，M=20）和圖3（hs=20，M=20）分別表示在兩個參數固定不變的情況下，斑塊數量隨空間尺度和光譜尺度變化的關系，從中可以得知，兩者的曲線的變化趨勢基本一致，尺度越大，所生成的斑塊數量越少；從兩條曲線的比較來看，在最小區(qū)域面積保持不變的情況下（M=20），顏色特征尺度hr比空間尺度hs對斑塊數量影響更敏感，隨著尺度的增加，斑塊的數量變化要大于空間尺度參數。

圖2 斑塊數量與尺度參數hs的相關曲線Fig.2 Relation between plaque numbers and scale parameter hs

圖3 斑塊數量與尺度參數hr相關Fig.3 Relation between plaque numbers and scale parameter hr

經過對上述3個分割參數的比較與分析之后，要達到準確提取杉木冠幅，確定最優(yōu)的分割尺度，除了能準確顯示冠幅邊界信息之外，還需考慮通過分割之后的斑塊不能太破碎，尤其是針對同質區(qū)域范圍內的影像，分割結果盡可能保持圖像的一致性，其次斑塊數量要控制在一定的范圍內，這樣可以避免產生過分割現象，經過反復試驗和比較，認為hs=10，hr=6，M=20對于研究區(qū)杉木冠幅提取具有較好的效果，共提取杉木有效冠幅樣本227個，使用邊界平滑算法對提取的冠幅進行平滑處理，效果如圖4所示。

圖4 Meanshift算法冠幅提取結果Fig.4 Extracted results of tree crows by using Meanshift

2.2 曲線估計模型

為了得到更為準確的模型，分別采用實測冠幅和影像提取的平均冠幅建立冠幅樹高模型。從提取的227個杉木冠幅樣本中，隨機抽取157個樣本作為建模數據，剩余70個樣本作為檢驗數據。采用曲線估計中的11種模型進行冠幅樹高建模，并選出最佳擬合模型，模型擬合結果如表1所示。從表中可以看出，利用實測冠幅信息建立的冠幅樹高模型的擬合最佳曲線為二次方程，決定系數R2為0.52。利用影像提取的冠幅信息建立的冠幅樹高擬合最佳曲線為冪函數模型，決定系數R2為0.36。

表 1 冠幅樹高最佳曲線估計模型Table 1 Optimal evaluation models for crown-height

圖5 實測冠幅樹高模型殘差分析Fig.5 Residual distribution of HT-CW model

圖6 影像冠幅樹高模型殘差分析Fig.6 Residual distribution of HT-PCW model

圖5 和圖6的殘差分布結果表明2個模型檢驗數據的殘差絕大部分都在置信帶內 [-2σ?，2σ?]，殘差點分布隨機性較好，以0為基準線上下對稱分布。從表2模型檢驗結果分析可知，利用影像提取的冠幅直接建立冠幅樹高模型的精度低于實測冠幅樹高建立的模型，因此利用影像冠幅直接建立樹高的遙感反演模型的效果并不理想。

表2 冠幅樹高模型檢驗結果Table 2 Examination results of regression models for crown-height

2.3 非線性聯立方程組模型

實測冠幅與樹高數據分析可知，冠幅與兩者之間存在較好的二次方程的關系。如何建立影像上的冠幅與實測樹高的反演模型呢？這需要從分析實測冠幅和影像提取冠幅之間的關系入手。首先，對多尺度均值漂移提取的227個平滑后的杉木冠幅與實測冠幅進行匹配，再建立曲線估計模型，從模型建立的情況來看，實測冠幅與遙感影像提取的冠幅存在良好的線性關系（圖7），其決定系數R2=0.646。

圖7 冠幅提取結果與實測值擬合效果Fig.7 Observed values of tree crown width against estimated values

由上可知，實測冠幅與影像冠幅存在線性關系，而實測冠幅樹高之間是二次方程的關系，是一種非線性關系。假定遙感影像冠幅沒有度量誤差，則可以通過實測冠幅建立影像冠幅與樹高的非線性聯立方程組反演模型。非線性聯立方程組的嚴格數學模型和推導方法見唐守正等（2009），研究所建模型的分線性聯立方程組運算通過ForStat軟件得到。從提取的冠幅227個樣本中，抽取157個樣本進行構建非線性聯立方程組反演模型，剩余70個樣本用來做模型精度檢驗，所得冠幅/樹高的非線性聯立方程組如下：

式中，聯立方程組的誤差變量指的是CW和HT。CW的殘差平方和為33.159，確定系數R2為0.642，H的殘差平方和為619.904，確定系數R2為0.355。利用樣本驗證數據來檢驗非線性聯立方程組模型預測值與實際觀測值是否存在較好的線性擬合關系，分析殘差的分布是否在置信帶[-2σ?，2σ?]內。冠幅樹高模型的檢驗結果如圖8所示。

圖8 非線性聯立方程組模型殘差分布Fig.8 Residual distribution of non-linear simultaneous model

從圖8的殘差分布可知，大部分驗證數據的殘差分布在置信帶 [-2σ?，2σ?]內。從表3 驗證數據的結果來看，驗證數據所得的冠幅樹高模型的決定系數0.419，與之前的冠幅與樹高直接建立的反演方程相比，并沒有得到明顯改善。模型檢驗數據估計值的標準差（SEE）、變動系數（CV）、平均百分標準誤差（MPSE）的統(tǒng)計結果與表3統(tǒng)計結果不差上下。從上述分析可知，假定影像冠幅提取的結果沒有度量誤差，建立影像冠幅與實測冠幅以及樹高的聯立方程組反演模型，沒有達到預期的建模目標，需要對非線性聯立方程組的建模因子做進一步的分析，提高模型的精度，使得模型解釋更為合理。

表3 非線性聯立方程組模型檢驗結果Table 3 Examination results of non-linear simultaneous equations

2.4 基于啞變量的非線性聯立方程組模型

在考慮樹高以及冠幅測量存在人為誤差建立的聯立方程組沒有得到良好的效果時，應該考慮樣地內林木樹高的差異，將樹高視為啞變量，進行數量化分級，建立基于啞變量的非線性聯立方程組反演模型。樣地實地調查的過程中，樣地中林木樹高差異較大，樣地內最高的杉木為18.6 m，最低的僅為7.2 m，為了便于計算將樹高分為4個等級，如表4所示。在ForStat軟件中構建啞變量非線性聯立方程組，將模型中所有參數的初始值設置為1，采用牛頓——拉格朗日方法進行計算。建模結果如下：

式中，CW表示實測平均冠幅，PCW表示影像上提取的冠幅長度，HT表示樹高，樹高的4個數量化等級表示如下：（1）當E1=E2=E3=0表示樹高第一個等級（I），（2）當E1=1且E2=E3=0時表示第二個等級（II），（3）當E2=1且E1=E3=0表示第三個等級（III），（4）當E3=1且E1=E2=0表示第四個等級（IV）。模型中，CW的殘差平方和為33.159，確定系數R2為0.642，H的殘差平方和為97.535，確定系數R2為0.899。

表4 樹高數量化分級Table 4 Quantitative hierarchy for DBH and HT

表5中模型檢驗結果表明，冠幅樹高模型的擬合效果很好，決定系數R2為0.885。而且，模型檢驗數據估計值的標準差（SEE）、變動系數（CV）、平均百分標準誤差（MPSE）的統(tǒng)計結果較非線性聯立方程組有很大的改進，誤差基本控制在10%以內。圖9中殘差點的分布隨機性較好，絕大部分驗證數據的殘差分布在置信帶 [-2σ?，2σ?]內。從上述分析可知，說明基于啞變量的非線性聯立方程組建立影像冠幅與樹高的反演模型效果十分理想。

表5 啞變量非線性聯立方程組模型檢驗結果Table 5 Examination results of non-linear simultaneous equations when dummy variable was used

圖9 啞變量非線性聯立方程組殘差分布Fig.9 Residual distribution of non-linear simultaneous equations when dummy variable was used

3 結論與討論

以Worldview-2影像對研究對象，以地面調查樣地為基礎，首先利面向對象的圖像分析方法開展杉木冠幅信息提取研究；其次，對提取的冠幅邊界信息進行平滑處理，結合實測冠幅和樹高建立聯立方程組反演模型。對模型預測值與實際觀測值是否存在較好的線性擬合關系進行了檢驗，分析了模型的殘差分布，得出以下幾點結論。

（1）面向對象的影像分割方法可以有效的提取杉木冠幅信息。研究運多尺度均值漂移分割方法對杉木冠幅信息進行提取。顏色特征尺度hr，空間尺度hs以及分割尺度（M）是多尺度均值漂移冠幅提取的關鍵，經過多次反復試驗，確定最合適的分割參數為hs=10，hr=6，M=20，多尺度均值漂移提取的有效斑塊數為227個，從斑塊提取的有效性比較來看，多尺度均值漂移算法是提取林木冠幅的一種有效方法。

（2）建立了影像冠幅樹高聯立方程組反演模型。影像分割平滑后的冠幅與實測冠幅存在良好的線性關系，決定系數R2=0.646。同時，實測冠幅與樹高最佳擬合曲線模型為二次方程。通過實測冠幅，建立了影像冠幅與樹高的非線性聯立方程組。其中，以考慮樹高以及實測冠幅存在測量誤差，且將樹高作為啞變量建立的影像冠幅樹高分線性誤差變量聯立方程組模型擬合效果最好，擬合效果比較理想，決定系數R2為0.899，模型的變動系數（CV）、平均百分標準誤差（MPSE）在10%以內。模型的適用性檢驗表明，基于啞變量的非線性聯立方程組模型的殘差點分布隨機性較好，以0為基準線上下對稱分布，說明利用從影像提取冠幅進行樹高的反演是可行的，可以為快速獲取杉木人工林平均胸徑、樹高等森林參數提供方法與技術參考。

[1] 曾 濤,琚存勇,蔡體久,等. 利用變量投影重要性準則篩選郁閉度估測參數[J]. 北京林業(yè)大學學報,2010,32(6):37-41.

[2] Franklin S., Wulder M., Gerylo G. Texture analysis of IKONOS panchromatic data for Douglas-f i r forest age class separability in British Columbia[J]. International Journal of Remote Sensing ,2002, 22(13): 2627-2632.

[3] Vander Sanden JJ and Hoekman DH. Review of relationships between grey-tone co-occurrence, semivariance, and autocorrelation based image texture analysis approaches[J].Canadian Journal of Remote Sensing, 2005, 31(3):207-213.

[4] Jon Pasher, Douglas J. King. Multivariate forest structure modeling and mapping using high resolution airborne imagery and topographic information[J]. Remote Sensing of Environment,2010, (114): 1718-1732.

[5] Marshall D. D.,Gregory P.J.,Hann DW. Crown prof i le equations for stand-grown western hemlock trees in northwestern Oregon[J]. Canadian Journal of Forest Research, 2003, (33): 2059-2066.

[6] Turan S?nmez. Diameter at breast height-crown diameter prediction models for Picea orientalis[J]. African Journal of Agricultural Research, 2009,4(3):215-219.

[7] 覃先林,李增元,易浩若.高空間分辨率衛(wèi)星遙感影像樹冠信息提取方法研究[J].遙感技術與應用,2005,20(2):228-232

[8] 馮益明,李增元,張 旭.基于高空間分辨率影像的林分冠幅估計[J].林業(yè)科學,2006,42(5):110-113

[9] Pouliot D A, King D J, Bell F W,et al.Automated tree crown detection and delineation in high-resolution digital camera imagery of coniferous forest regeneration[J]. Remote Sensing of Environment. 2002,(82):322-334

[10] Leckie D G, Gougeon F A, Tinis S,et al. Automated tree recognition in old growth conifer stands with high resolution digital imagery[J]. Remote Sensing of Environment, 2005, (94):311-326

[11] 黃建文,鞠洪波,趙 峰,等. 利用遙感進行退耕還林成活率及長勢監(jiān)測方法的研究[J]. 遙感學報.2007,11(6):899-905

[12] 鄧 廣. 高空間分辨率遙感影像單株立木識別與樹冠分割算法研究[J]. 北京中國林業(yè)科學研究院博士論文,2009,8-10

[13] Wang Le. A Multi-scale approach for delineating individual tree crowns with very high resolution imagery[J]. Photogrammetric Engineering & Remote Sensing,2010,76(4):371-378

[14] 孫 華,鞠洪波,張懷清,等. 三種回歸分析方法在Hyperion影像LAI反演中的比較[J]. 生態(tài)學報,2012,32(24):7781-7790

[15] K. Fukunaga,L.D. Hostetler. The Estimation of the Gradient of a Density Function, with Applications in Pattern Recognition[J].IEEE Transactions on Information Theory,1975,21(1):32-40

[16] Dorin Comaniciu,Peter Meer. Mean shift:A robust approach toward feature space analysis[J]. IEEE Transaction on pattern analysis and machine intelligence,2002,24(5):601-619

[17] 唐守正,郎奎建,李?？? 統(tǒng)計和生物數學模型計算(ForStat教程)[M].北京:科學出版社,2009,197-290.

[18] 李曉冬, 王雪峰, 賀 鵬. 林木冠層圖像的分割方法研究[J].中南林業(yè)科技大學學報, 2013, 33(7): 40-44.

[19] 劉 剛, 張 雨, 臧 卓, 等. 多源遙感數據森林信息的提取和比較分析[J]. 中南林業(yè)科技大學學報, 2012, 32(10): 158-161.

Study on crown diameter extraction and tree height inversion based on worldview-2 data

SUN Hua1,2, JU Hong-bo2, ZHANG Huai-qing2, LING Cheng-xing2
(1. Research Center of Forestry Remote Sensing & Information Engineering, Central South University of forestry & Technology,Changsha 410004, Hunan, China; 2. Research Institute of Forest Resources Information Technique, Chinese Academy of Forestry,Beijing 100091, China)

By taking the Chinese fi r plantation in Huangfengqiao Forest Farm in Youxian county, Hunan province as the tested objective,the methods of tree crown width extraction and tree height inversion were studied. Based on worldview-2 data, the tree crown width information of the plantation was extracted by adopting mean shift segmentation algorithm. After setting a series of Image segmentation scale, the optimal crown width segmentation parameters’ values were determined, they arehs=10,hr=6,M=20.The extracted crown width border was smoothened, and the estimation model of tree crown=height curve and the inversion model of nonlinear simultaneous equations were respectively established based on the smoothed image’s crown width and measured tree height. Of them, the model of image tree crown-high nonlinear simultaneous equations set by taking tree high as the dummy argument had the optimal fi tting results,the model coeff i cient of determinationR2=0.899, the variation coeff i cient of modelCCVand the average percentage standard errorsEMPSEboth were less than 10%, so the last selected model is an effective means for tree height inversion of Chinese fi r plantation.

forest remote sensing; tree crown width; Chinese fir plantation; image segmentation; mean shift algorithm; non-linear simultaneous equations

S771.8

A

1673-923X(2014)10-0045-06

2013-12-21

國家“十二五”863課題：“數字化森林資源監(jiān)測關鍵技術研究”（2012AA102001）；林業(yè)公益性行業(yè)科研專項（201104028）

孫 華（1979-），男，湖南隆回人，博士，講師，研究方向為林業(yè)信息技術

張懷清（1973- ），男，湖南寧鄉(xiāng)人，研究員，碩士生導師，主要從事林業(yè)可視化模型技術與濕地監(jiān)測技術

[本文編校：吳 彬]