蘇迪，高心丹

（東北林業(yè)大學信息與計算機工程學院，哈爾濱 150040）

森林是地球生態(tài)系統(tǒng)中重要的組成部分，其碳匯功能在維持全球碳平衡方面發(fā)揮著重要作用。森林蓄積量是生物量和碳儲量研究的重要參考因素，是評價森林固碳能力的重要指標，也是森林資源調(diào)查的重要因素。傳統(tǒng)的一、二類森林資源調(diào)查是測算森林蓄積量的重要方法，這兩種方法雖然準確度較高，但在時間與空間上存在較大的局限性，調(diào)查過程會耗費大量的時間和人力［1-2］。如何在準確獲取森林蓄積量信息的同時，快速、高效地為區(qū)域生態(tài)狀況評估提供參考［3-5］已成為研究熱點。

近年來，隨著遙感技術的發(fā)展，基于遙感因子的森林蓄積量估測方法因其可以快速、高效獲得蓄積量模型而引起了國內(nèi)外學者的廣泛關注［6-8］?；貧w建模估測法因其工作量相對較少、估測頻率高、估測覆蓋面積大等優(yōu)點逐步成為基于遙感因子蓄積量估測研究的主要方向［9-11］。按照自變量數(shù)量可以將蓄積量回歸模型分為一元和多元回歸，按照方程表現(xiàn)形式可以分為線性和非線性回歸［12-13］。研究森林蓄積量的線性回歸方法主要有最小二乘法、逐步回歸法、主成分回歸法、偏最小二乘法［14-16］。其中，偏最小二乘法結合了多種回歸方法的功能和優(yōu)點，有效地改善了模型自變量間的多重共線性問題，目前已廣泛地應用于各領域的估測研究中［17-20］。帶有點云數(shù)據(jù)的航測數(shù)據(jù)與其他遙感數(shù)據(jù)相比，具有分辨率高、細節(jié)更完善、時效性強、成本低等優(yōu)點，可以更準確地提取蓄積量的相關特征因子［21-23］，從而提高蓄積量估測精度。

本研究利用無人機航攝影像點云數(shù)據(jù)估測樹高和胸徑因子，使用正射影像通過分水嶺算法提取樹冠個數(shù)，結合提取的坡向等因子，用主成分回歸方法估測郁閉度，結合全部因子通過偏最小二乘回歸方法估測森林蓄積量。

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域

以帽兒山林場老山施業(yè)區(qū)森林為研究區(qū)域（127°18′0″～127°41′6″E，45°2′20″～45°18′16″N），該林場是東北林業(yè)大學實驗林場，位于尚志市西北部，場址距市區(qū)40 km，隸屬于尚志國有林場管理局。林業(yè)用地面積26 067 hm2，其中有林地面積23 204 m3，森林總蓄積1 879 380 hm2，森林覆蓋率83.29%。老山施業(yè)區(qū)屬于帽兒山林場的一部分，原始地帶植被為闊葉紅松混交林，現(xiàn)存植被以紅松、樟子松、落葉松人工林為主，該區(qū)域具有豐富的森林資源適合進行蓄積量研究。

1.2 研究數(shù)據(jù)

研究數(shù)據(jù)1 來源于航測數(shù)據(jù)，原始圖像數(shù)據(jù)來源于型號為DJI X3 的無人機拍攝而成的影像，機上搭載Canon?EOS?60D 數(shù)碼相機拍攝，有效像素1 240 萬，最大分辨率可達到4 000×3 000 ppi，成像時間為2015 年9 月14 日，利用航空影像中獲得DOM（digital orthophoto map，數(shù)字正射影像圖），DOM 影像面積達到4 723.169 7 hm2，其地理坐標為xian80，研究區(qū)域內(nèi)的DOM 與小班矢量圖疊加圖像見圖1。用lasmerge 軟件將多個點云數(shù)據(jù)合成為一個點云數(shù)據(jù)，按照研究數(shù)據(jù)2 實測范圍對點云數(shù)據(jù)進行裁剪，其地理坐標為xian80，其面積為3 129.276 hm2，平均點云密度為49 point/m2，按高程顯示的點云數(shù)據(jù)圖像見圖2。

圖1 DOM 與小班矢量圖疊加圖像Fig.1 DOM and sub?compartment vector overlay image

圖2 點云數(shù)據(jù)圖像Fig.2 Point cloud data image

研究數(shù)據(jù)2 是實測數(shù)據(jù)，為帽兒山林場老山實驗區(qū)森林資源二類調(diào)查數(shù)據(jù)。其中包括2015 年調(diào)查的1∶10 000 地形圖，618 個以矢量地理信息數(shù)據(jù)形式儲存的小班二類調(diào)查數(shù)據(jù)。

研究數(shù)據(jù)3 為實測數(shù)據(jù)，為東北林業(yè)大學帽兒山實驗林場森林固定樣地數(shù)據(jù)，2004 年調(diào)查更新的比例尺1∶25 000 的固定樣地分布圖，調(diào)查總面積為26 496 hm2，包含263 個樣地的二類調(diào)查信息。

1.3 樣地篩選

不同類型的GIS 因子數(shù)據(jù)對蓄積量模型精度有嚴重的影響，為了減少該影響，對研究數(shù)據(jù)2、3中的針闊混交林小班采用標準差分析法剔除樣本中離群較大的數(shù)據(jù)，即剔除樣本中的數(shù)據(jù)。重復上述步驟，最終選擇到120 個針闊混交林樣本數(shù)據(jù)。將所有數(shù)據(jù)進行中心標準化處理，去除因各個因子間量綱不同帶來的影響，統(tǒng)一所有因子的量綱，如式（1）所示：

2 特征參數(shù)提取與結果

利用GIS 技術提取與森林蓄積量相關的樹冠個數(shù)、平均樹高、平均胸徑、坡度、坡向、海拔、小班面積、郁閉度這8 個GIS 因子建立蓄積量模型，通過該方法提高蓄積量的估測精度。本研究先使用ArcGIS10.2 軟件將DOM 圖像和點云數(shù)據(jù)圖像按照小班邊界進行裁剪，為進行特征因子提取做準備。本研究結合點云數(shù)據(jù)、DOM 圖像和研究區(qū)小班矢量數(shù)據(jù)，以二類調(diào)查小班為單位對樹冠個數(shù)、平均樹高等因子進行估測，對海拔、坡向等因子進行提取，獲得建立蓄積量模型所需的特征因子，特征因子提取流程圖見圖3。

圖3 特征因子提取流程圖Fig.3 Feature factor extraction flow chart

2.1 特征參數(shù)的提取

2.1.1 分水嶺樹冠個數(shù)提取

分水嶺算法適合提取邊緣微弱較敏感、分辨率高、細節(jié)多的影像輪廓，但傳統(tǒng)分水嶺算法存在過分割問題，為了解決該問題，本研究將120 個小班正射影像作為分水嶺分割的原始圖像，首先使用支持向量機的分類方法將其分為針葉林和闊葉林；其次采用最大熵法分別獲得針葉林和闊葉林的最優(yōu)閾值，比較針、闊葉林的最大熵最優(yōu)閾值，發(fā)現(xiàn)針葉林最優(yōu)閾值大于闊葉林最優(yōu)閾值；最后對形態(tài)學開閉重建處理后的針闊混交林圖像進行基于閾值標記的分水嶺分割，第1 次將小于闊葉林閾值區(qū)域進行標記和分水嶺分割，第2 次將大于闊葉林最優(yōu)閾值，小于針葉林最優(yōu)閾值進行分割得到樹冠輪廓圖像。

2.1.2 樹高和胸徑的估測

使用LIDAR360 工具對點云數(shù)據(jù)進行去噪處理，生成0.5 m 分辨率的DEM（digital elevation model，數(shù)字高程模型）和DSM（digital surface model，數(shù)字地表模型）。由于DEM 模型中只有地形的高程信息，不包括DSM 模型中的森林樹木高度等地表信息的高程，因此從DSM 中減去DEM 即可獲得CHM（canopy height model，冠層高度模型）。

CHM 中每個點的像素值即為該點的高程值，通過MATLAB 2013b 軟件編程，讀取各個小班樣本的CHM 模型平均高程，只保留高程值不為0 的點進行平均高程的計算。由于CHM 為冠層模型不能等同于樹高，因此需要結合讀取的CHM 平均高程利用一元線性回歸方法建立平均樹高模型和平均胸徑模型。

2.1.3 海拔、坡向、坡度、小班面積因子的提取

使用LIDAR360 軟件在DEM 基礎上生成等高線圖（圖4a）。其中比例尺為1∶1 000、三角形最大邊為30 m、間曲線等高距為2.5 m、首曲線等高距為5 m、計曲線等高線為25 m。利用DEM 數(shù)據(jù)生成坡向圖像和以10°為間隔進行處理生成坡度圖像，如圖4b、c 所示。

使用Python2.0 編程，用雙線性插值法提取等高線、坡向和坡度平均值。使用ArcGIS10.2 軟件結合小班矢量圖計算小班面積，提取各項因子為接下來建立郁閉度和蓄積量模型提供特征因子。

2.1.4 郁閉度估測

圖4 等高線、坡向、坡度圖像Fig.4 Contour line，aspect，slope image

主成分回歸是主成分分析和回歸估計結合的方法，主成分分析是一種變量降維技術，是將多個相關變量重新組合成幾個主成分變量，減少變量個數(shù)的變量分析方法。該方法能夠簡化分析過程，提高分析效率。同理主成分回歸是用合成后的主成分變量進行回歸分析，不僅保留了原自變量對因變量的影響，還消除了各個變量相互間的共線性問題。本研究通過主成分回歸將樹冠個數(shù)、平均樹高、平均胸徑、坡度、坡向、海拔、小班面積這7 個特征因子組成若干相互獨立的主成分，根據(jù)各個主成分的貢獻率大小選取模型最終的主成分變量，建立郁閉度回歸模型。使用MATLAB 2013b 編程，分析1～70 號小班數(shù)據(jù)，得出了各主成分的特征根、主成分得分、主成分累積貢獻率。按照貢獻率累積80%以上的選取標準，根據(jù)貢獻率的大小順序選取主成分，主成分與特征因子間的關系方程，再根據(jù)主成分的構成，還原郁閉度方程，得到郁閉度與7個特征因子的關系模型。

2.2 特征參數(shù)的提取結果

2.2.1 樹冠個數(shù)提取結果

實驗使用2.1.1 部分改進后分水嶺分割算法提取樹冠輪廓信息，將樹冠分割結果圖像進行二值化處理計算連通區(qū)域得到樹冠個數(shù)，通過公式（2）計算分割準確率和平均分割準確率來顯示分割效果。

式中：Ad為準確率；為平均準確率；Nc為分割的樹冠個數(shù)；Nd為小班樹木株數(shù)總數(shù)（實測小班株數(shù)）；n為圖像個數(shù)。

通過計算得到120 個樣地DOM 圖像的傳統(tǒng)分水嶺平均分割精度為52.89%，改進后分水嶺平均分割精度為80.03%，能夠達到基本分割要求，改善了傳統(tǒng)分水嶺存在的過分割現(xiàn)象，分割后得到的樹冠個數(shù)可以用于建立蓄積量模型。經(jīng)計算，23 號樣地樹冠分割精度與平均分割精度最相近，所以選擇該樣地圖像進行顯示（圖5）。

圖5 樣地23 圖像Fig.5 Sample plot 23 image

2.2.2 平均樹高、平均胸徑估測結果

根據(jù)2.1.2 部分得到平均樹高和平均胸徑模型結果如下：

由圖6 和平均樹高、平均胸徑模型的平均精度可知，估測值與實測值的擬合程度較好，且都是由CHM 直接進行估測而成，充分利用點云數(shù)據(jù)所生成的CHM 數(shù)據(jù)。

圖6 平均數(shù)高和平均胸徑擬合曲線Fig.6 Average height and average DBH curve

2.2.3 郁閉度估測結果

由2.1.4 部分將7 個變量進行主成分分析，其中第1 主成分貢獻率32.087%，第2 主成分貢獻率24.375%，第3 主成分貢獻率18.005%，第4 主成分貢獻率13.769%。前4 個主成分貢獻率已達到88.235%，符合主成分選取要求，選取前4 個主成分因子為新變量。則公式（5）～（8）為4 個主成分與特征因子間的關系方程。

第1 主成分為：

第2 主成分為：

第3 主成分為：

第4 主成分為：

公式（5）～（8）中：X1 為小班平均樹高；X2 為平均胸徑；X3 為坡度；X4 為坡向；X5 為海拔；X6 為小班面積；X7 為分水嶺提取的樹冠個數(shù)。

主成分回歸在減輕原始變量的多重共線性問題后，得到以4 個主成分為新自變量的郁閉度模型，回歸方程為：

根據(jù)4 個主成分的構成，將回歸方程還原，得到關于自變量X1—X7 的回歸方程為：

方程（9）和（10）中X8 為郁閉度，通過71～120號50 個檢驗數(shù)據(jù)進行精度檢驗，郁閉度模型（公式10）精度為83.18%，可以達到精度要求，再結合前文所有提取和估測的特征因子建立蓄積量模型，進一步研究各個特征因子與蓄積量之間的關系。

3 蓄積量估測模型

3.1 偏最小二乘回歸方法

偏最小二乘回歸方法是一種集典型相關分析和主成分分析優(yōu)點以及功能于一身的多元線性分析方法。該方法由于精度高，穩(wěn)健性實用性好而被廣泛地應用于各種估測研究中。

設已知單個因變量Y和自變量X＝［x1，x2，…，xp］，樣本個數(shù)n，偏最小二乘法從X和Y的相關矩陣中提取主成分t，使用Y和X對t進行回歸，對于回歸分析的需要，偏最小二乘回歸在提取成分時添加以下兩個目標：1）盡可能多地在t中攜帶X矩陣數(shù)據(jù)表中的變異信息；2）t和Y的相關程度能夠達到最大。在對X和Y進行一次成分提取后，分別進行X對t的回歸和Y對t的回歸。如果回歸方程達到了滿意的精度，就停止計算，否則將X關于t回歸后的殘差矩陣與Y關于t回歸后的殘差矩陣進行新一輪的成分提取。反復重復這個過程，直至滿足交叉有效性原則規(guī)定的條件，最終確定提取主變量成分的個數(shù)，建立偏最小二乘回歸方程［19］。

3.2 蓄積量模型的建立

借助SPSS 軟件，通過相關性分析得到與蓄積量存在顯著關系的8 個自變量如表1 所示，表中X1—X8 同前文。

表1 表明各因子間存在多重相關性，一般回歸模型已不適用，偏最小二乘法對自變量的選擇門檻不高，不需要選擇最優(yōu)因子，且較多變量有益于對提取的主成分進行累計解釋分析，因此本研究采用偏最小二乘法建立森林蓄積量估測模型。

在MATLAB 2013b 軟件中按照偏最小二乘原理編寫PLS 回歸代碼，對1～120 號樣本用偏最小二乘法建立蓄積量模型，對模型進行交叉有效檢驗，得到原始變量的回歸方程系數(shù)，據(jù)此可以擬合出基于偏最小二乘回歸方法蓄積量方程。

表1 特征因子相關系數(shù)Table 1 Characteristic factor correlation relationship table

式中：V為蓄積量，方程包含的X1—X8 與上文一致，將30 個固定樣地檢驗樣本數(shù)據(jù)代入方程（公式11），得到蓄積量估測值，為模型評價和精度檢驗做準備。

3.3 模型評價和精度檢驗

結合預留的121～150 號30 個驗證樣本實測數(shù)據(jù)和蓄積量方程（公式11）估算的森林蓄積量得到圖7 實測值與估測值折線圖。

圖7 實測值與估測值折線圖Fig.7 Line chart of stock volume measured value and estimated value

從圖7 中可以看出估測的蓄積量值整體偏大，通過對蓄積量實測值和估測值進行配對T檢驗，得到T＝-0.374，P＝0.657＞0.05，說明實測值和估測值相比雖然整體偏小，但不存在顯著性差異。

使用實測值與估測值對模型進行評價，根據(jù)模型的決定系數(shù)R2，均方根誤差RMSE，總體相對誤差RS，平均相對誤差MRE，預估精度Pr為評價指標。根據(jù)公式（12）～（16）計算：

式中：yi為蓄積量實測值；為蓄積量模型預測值；為蓄積量實測樣本平均值；n為檢驗樣本數(shù)；i表示第i個樣本。本研究蓄積量模型（公式11）的R2＝0.738 6，RMSE＝5.135 3 m3/hm2，說明模型擬合效果較好。利用模型偏差統(tǒng)計量進行比較并評價模型的預測能力，RS＝-1.285 8、MRE＝-0.263 0，說明估測蓄積量與實測蓄積量偏離不大，可以很好地進行蓄積量估測；精度達到88.43%，估測能力較好，并且有一定的可推廣性。

以蓄積量實測值為橫坐標，估測值為縱坐標，建立實測值和估測值散點圖，建立的散點圖趨勢線擬合效果較好，實測值和估測值比較吻合（圖8）。

圖8 蓄積量實測值與估測值的散點圖Fig.8 Scatter plot of measured and estimated values of stock volume

4 結論

利用無人機影像生成的點云數(shù)據(jù)和正射影像估測蓄積量的方法，能夠在充分利用無人機影像的同時最大限度地保留無人機影像的細節(jié)和各種因子的特征，提高特征參數(shù)的提取精度。該方法不僅提高了實施效率、降低了估測成本，還為代替人工野外實測森林調(diào)查信息提供了可能。無人機數(shù)據(jù)生成的正射影像與高分辨率衛(wèi)星影像相比具有重疊度大、分辨率高等優(yōu)點；其生成的點云數(shù)據(jù)不僅可以代替機載雷達獲得高程信息，還具有航線靈活、操作方便、時效性好等優(yōu)點。無人機航測數(shù)據(jù)彌補了高分辨率衛(wèi)星影像無法同時獲得影像和高程信息的缺點，在估測森林生物量和蓄積量方面應用前景廣闊。

本研究利用無人機航攝方法采集了黑龍江帽兒山林場老山實驗區(qū)航攝影像，利用三維點云數(shù)據(jù)提取CHM 模型高程，估測了平均樹高和平均胸徑；用改進的分水嶺算法獲得了樹冠個數(shù)；用DEM數(shù)據(jù)獲得坡度等特征信息，估測了森林蓄積量。使用帽兒山林場固定樣地數(shù)據(jù)作為檢驗樣本，結果表明，平均樹高提取精度達到97.34%，平均胸徑提取精度達到91.27%，樹冠個數(shù)提取精度達到80.03%，郁閉度模型精度為83.18%，森林蓄積量估測精度達到88.43%，滿足森林調(diào)查精度要求。研究得到如下結論：1）在提取CHM 平均高程值時不計算高程值為0 的點，減少了點云數(shù)據(jù)空白點對平均樹高、胸徑模型的影響；2）通過遙感技術提取特征因子，用主成分回歸的方法可有效地估測郁閉度；3）模型中所有因子都是通過遙感與回歸方法提取或估測生成，使用偏最小二乘回歸方法建立蓄積量模型，為信息化、自動化提取特征因子估測蓄積量提出可能。由于時間和試驗條件限制，仍存在很多不足和改進的地方，如波段信息對蓄積量估測的影響等方面有待進一步研究。