基于高分六號衛星紅邊波段的森林蓄積量遙感反演
——以西寧市針葉林為例*

任 楓，王 琦，楊 佳，任慶福**

（1.西安綠環林業技術服務有限責任公司，西安 710048；2.廣東省科學院生態環境與土壤研究所，廣州 510650；3.航天宏圖信息技術股份有限公司，北京 100089）

森林蓄積量是指森林中全部樹木材積的總和，它是反映一個地區森林資源的豐富程度、衡量森林生態環境優劣、進行森林經營和采伐利用的重要依據[1]。傳統森林蓄積量的估算是通過林分調查，建立材積方程，進而推算區域森林蓄積量，往往存在工作量大、效率低下、成本高等問題。隨著遙感技術的快速發展，森林蓄積量的遙感反演已成為獲取區域森林蓄積量的主要手段[2]。目前，光學遙感的森林蓄積量反演，多數研究采用Landsat TM/ETM/OLI[3-4]、MODIS[5]、SPOT5[6]、資源三號[7]、GF1/2[8-10]等衛星載荷，通過提取與蓄積量相關的因子（如地形因子、光譜信息、植被指數、影像紋理等）建立模型進行反演，因此，在建模過程中，敏感因子的選擇就成為影響蓄積量反演精度的關鍵因素。

一直以來，紅光波段和近紅外波段常被作為指示植被生長狀況的敏感波段，但有研究指出，介于紅光與近紅外波段之間的紅邊波段，對植被葉綠素更加敏感[11-13]。受衛星載荷的波段限制，紅邊波段在植被監測方面的研究局限于基于地物光譜儀[14-15]。近年來，有較多的衛星載荷（如美國的WorldView-2/3、歐空局的Sentinel-2、德國的RapidEye等）增加了紅邊波段，基于紅邊波段的森林蓄積量遙感反演的研究隨之也不斷問世[16-18]。這些研究多以Sentinel-2為遙感數據源，以模型的適用性評價為主要研究目標，極大豐富了蓄積量反演研究的數據源和模型方法，而關于紅邊波段對蓄積量反演精度的影響缺乏深入研究，加之光學遙感影像數據的有效性易受云霧影響，因此，僅僅以Sentinel-2為數據源極大限制了紅邊波段在森林蓄積量反演方面的研究與應用。國產高分六號衛星的寬幅影像（GF6-WFV）是中國首顆帶有紅邊波段的中高分辨率影像。但目前，基于GF6-WFV紅邊波段的研究主要集中在植被識別及分類[19-24]、森林擾動[25]等方面，在森林蓄積量反演方面的研究未見報道。因此，基于GF6-WFV衛星紅邊波段的森林蓄積量反演研究，對于明確紅邊波段對森林蓄積量反演的作用以及拓展國產高分衛星的應用具有重要的理論指導意義。

西寧市地處黃土高原與青藏高原交界地帶，在生態環境方面具有極其重要的地位。但該地區存在生態防護功能差、森林生態系統脆弱、林地質量差等問題[26-27]，迫切需要對全區的森林蓄積量進行動態監測。本研究以GF6-WFV為遙感數據源，結合西寧市DEM及2014年森林資源二類調查數據，從光譜特征、植被指數、地形因子、影像紋理4個方面選取了蓄積量反演的自變量，并在此基礎上采用多元線性回歸、隨機森林模型，研究紅邊波段對西寧市針葉林蓄積量遙感反演精度的影響，以期為西寧市森林蓄積量的區域動態監測提供方法借鑒，為紅邊波段在森林蓄積量遙感反演中的應用提供理論依據，為拓展國產高分六號遙感衛星的應用提供新的途徑。

1 資料與方法

1.1 研究區概況

研究區位于青海省東部的西寧市，地理位置介于100°52'7″-101°54'58″E、36°13'40″-37°28'9″N，平均海拔3127.8m。研究區地處黃河一級支流湟水河中游，是黃土高原向青藏高原的過渡地帶，地形復雜，以川水、淺山、腦山3種地形為主，具有典型的干旱高原大陸性氣候，年平均氣溫為5.8℃，極端最低氣溫為-26℃， 極端最高氣溫為41℃，無霜期約160d。年平均降水量368mm，集中于7-9月；年蒸發量為1100mm，年日照時數2600h。研究區分布的主要喬木樹種有青海云杉（Picea crassifolia）、落葉松（Larix principis-rupprechtii Mayr）、油松（Pinus tabulaeformis Carr.）、祁連圓柏（Sabina przewalskii）、白樺（Betula platyphylla Suk）、青楊（Populus cathayana Rehd.）、旱柳（Salix matsudana Koidz）等。

1.2 數據來源及預處理

1.2.1 遙感影像數據及預處理

遙感影像數據來源于自然資源衛星遙感云服務平臺（http://www.sasclouds.com/chinese/normal/），選用2019年8月15日覆蓋西寧市的高分六號衛星寬幅影像（簡稱GF6-WFV），其衛星具體參數見表1。采用國產遙感處理軟件PIE-Basic6.0（https://www.piesat.cn/product/PIE-Basic/index.htMLR）對GF6-WFV影像進行輻射定標、幾何校正、大氣校正等預處理，最終得到各波段的地表反射率數據。其中，影像的幾何校正處理所需的高程數據來源于地理空間數據云（http://www.gscloud.cn/）的GDEMV2 30M分辨率數字高程數據集。從該數據集中下載覆蓋GF6-WFV整 幅 影 像（93°00'00″- 103°00'00″E、35°00'00″-40°00'00″N）的DEM數據，DEM的分幅數據拼接采用ArcGIS10.2軟件完成。

表1 GF6-WFV衛星參數Table 1 The satellite parameters of GF6-WFV

1.2.2 針葉林蓄積量

針葉林蓄積量的地面驗證數據來源于西寧市2014年森林資源二類調查數據。由于研究區地處西北干旱區，針葉樹種生長速度較慢，且與影像獲取時間滯后在一個齡級（10a）內，因此，忽略影像獲取時間與蓄積量調查時間的滯后效應。利用三倍標準差分析方法剔除針葉林蓄積量異常的調查小班，共獲得503個調查小班的針葉林蓄積量數據，得到西寧市針葉林調查小班（平均胸徑在5cm及以上）的分布（圖1）。

圖1 GF6-WFV影像及針葉林調查小班分布Fig.1 GF6-WFV imagery and the sub-compartments of coniferous forest in study area

1.3 自變量選取與確定

1.3.1 植被指數

以GF6-WFV影像前6個波段（即表1中的B1、B2、B3、B4、B5和B6）的地表反射率作為單波段反射率變量。選取歸一化植被指數（NDVI）、比值植被指數（RVI）、增強植被指數（EVI）和土壤調節植被指數（SAVI）4個常規植被指數，在此基礎上選取基于紅邊波段的2個植被指數即MTCI指數[28]和NDREI指數[29]，作為針葉林蓄積量反演的植被指數變量，具體計算式見表2。

表2 植被指數計算公式Table 2 The formula of different vegetation index

1.3.2 地形及紋理特征變量

選取坡度（Slope）、坡向（Aspect）作為地形因子變量。影像紋理特征的描述采用灰度共生矩陣方法（Gray-Level Co-occurrence Matrix, GLCM），選取均勻性（H）、異質性（D）、對比度（CT）、能量（ENG）、相關性（CR）、信息熵（ENT）6個指標來刻畫GF6-WFV影像的紋理特征。

采用滑動窗口法（窗口大小3×3）分別計算GF6-WFV影像前6個波段的6個紋理特征，得到研究區內36張紋理特征變量的分布圖，然后分別計算各調查小班范圍內的紋理特征變量的平均值，得到503個調查小班對應的36個紋理特征變量值。

1.3.3 自變量分組與確定

（1）將選取的6個單波段反射率變量、6個植被指數變量、2個地形因子變量以及36個紋理特征變量有重復地分為無紅邊波段處理（No Red-Edge）和加入紅邊波段處理（Red-Edge Added）兩組。

在無紅邊波段的處理中，自變量包括ρ1、ρ2、ρ3、ρ4、NDVI、RVI、EVI、SAVI、Slope、Aspect以及影像前4個波段對應6個紋理指標的24個紋理特征變量；在加入紅邊波段的處理中，其自變量是在無紅邊波段處理的自變量基礎上，增加GF6-WFV影像第5、第6紅邊波段相關的自變量，即ρ1、ρ2、ρ3、ρ4、ρ5、ρ6、NDVI、RVI、EVI、SAVI、Slope、Aspect以及影像前6個波段對應6個紋理指標的24個紋理特征變量。

（2）采用主成分分析法（Principal Component Analysis）對兩組處理的紋理特征變量進行降維處理，獲得新的變量，作為紋理特征的備選變量。

（3）分別將兩組處理的備選自變量與針葉林蓄積量建立逐步回歸方程，根據AIC最小準則，最終確定兩組處理的自變量。

1.4 模型構建和檢驗

1.4.1 模型選取

從總體樣本（503個調查小班的蓄積量及對應自變量）中隨機抽取75%的數據作為模型的訓練樣本，其余25%作預測樣本。選取線性、非線性模型分別對無紅邊波段（No Red-Edge）、加入紅邊波段（Red-Edge Added）的兩組處理進行蓄積量反演，模型的構建及精度評價通過Python語言實現。

線性反演模型采用多元線性回歸模型（MLR）,非線性反演模型采用隨機森林模型（RF）。

隨機森林模型是一種基于決策樹的機器學習算法[30]，它是在決策樹的基礎上，利用bootstrap方法從訓練集 X = { x1, x2,… , xn}中有放回的多次重復采樣，然后結合對應的目標值對這些樣本訓練模型，在訓練結束之后，對未知樣本的預測可以通過簡單多數原則確定檢驗樣本的反演結果。

1.4.2 模型構建及精度評價

（1）基于兩組處理的訓練樣本，分別對MLR、RF模型的參數進行訓練，以決定系數為代價函數，得到各模型的最優參數。

（2）將得到的最優參數帶入MLR、RF模型，對兩組處理的預測樣本進行蓄積量的預測。

（3）通過十折交叉驗證方法，分別將兩組處理的MLR、RF模型預測的蓄積量與預測樣本的蓄積量進行對比，計算兩者之間的決定系數（R2）和均方根誤差（RMSE），以衡量模型反演的性能及精度。計算式為

式中，n為預測樣本個數，y?i為模型預測的蓄積量（m3·hm-2）；yi為預測樣本的蓄積量（m3·hm-2）。

2 結果與分析

2.1 紋理特征變量降維

無紅邊波段處理（No Red-Edge）的24個紋理特征變量包括GF6-WFV影像前4個波段（B1、B2、B3、B4）對應的6個紋理指標（H、D、CT、ENG、CR、ENT）；加入紅邊波段處理（Red-Edge dded）的36個紋理特征變量包括GF6-WFV影像前6個波段（B1、B2、B3、B4、B5、B6）對應的6個紋理指標（H、D、CT、ENG、CR、ENT），采用主成分分析法分別對兩組處理的紋理特征變量進行降維處理，選擇方差累計貢獻率大于95%的主成分變量代替原有的紋理特征變量。

經計算，無紅邊波段處理的紋理特征變量中提取到9個主成分（記為P1N, P2N, …, P9N），其方差累計貢獻率為95.59%（表3）。對主成分中24個變量的系數進行排序，選取系數最大的前3個變量，對9個主成分逐個分析，從波段方面來看，第1、2、3個主成分的方差貢獻率分別為45.44%、19.78%、10.06%，分別解釋了B2、B3、B4波段的紋理特征；紋理指標以異質性（D）、均勻性（H）、能量（ENG）以及信息熵（ENT）為主。

表3 無紅邊波段處理的紋理特征變量主成分分析結果Table 3 The PCA result of the texture feature variables for the No Red-Edge treatment

加入紅邊波段處理的紋理特征變量中提取到13個主成分（記為P1A, P2A, …, P13A），其方差累計貢獻率為95.45%（表4）。加入紅邊波段之后，第1主成分的前3個變量與無紅邊波段的處理一致，但其方差貢獻率降為36.09%；第2主成分的方差貢獻率為18.38%，且增加了B6波段的均勻性（H）；第3主成分則反映了B5波段的紋理特征，以能量（ENG）、 均勻性（H）、異質性（D）為主要指標，其方差貢獻率為10.27%。因此，對于無紅邊波段處理的紋理特征變量，以9個主成分（即P1N, P2N, …, P9N）替代原有的24個特征變量；對于加入紅邊波段處理的紋理特征變量，則以13個主成分（即P1A, P2A, …, P13A）替代原有的36個特征變量。

表4 加入紅邊波段處理的紋理特征變量主成分分析結果表Table 4 The PCA result of the texture feature variables for the Red-Edge Added treatment

2.2 自變量篩選

無紅邊波段處理包括4個單波段地表反射率變量（1ρ、ρ2、ρ3、ρ4）、4個植被指數變量（NDVI、RVI、EVI、SAVI）、2個地形因子變量（Slope、Aspect）以及9個紋理主成分變量（P1N, P2N, …, P9N）共19個備選自變量；加入紅邊波段處理包括6個單波段地表反射率變量（1ρ、ρ2、ρ3、ρ4、ρ5、ρ6）、6個植被指數變量（NDVI、RVI、EVI、SAVI、MTCI、NDREI）、2個地形因子變量（Slope、Aspect）以及13個紋理主成分變量（P1A, P2A, …, P13A）共27個備選自變量。分別將各備選自變量與針葉林蓄積量建立逐步回歸方程，依據AIC最小準則，得到兩組處理最終篩選出的自變量（表5）。表5中，無紅邊波段處理篩選出11個變量（ρ2、ρ3、ρ4、NDVI、SAVI、Slope、P1N、P2N、P4N、P5N、P9N），而加入紅邊波段處理篩選出13個變量（1ρ、ρ2、ρ5、MTCI、Slope、P1A、P2A、P3A、P5A、P6A、P7A、P9A、P11A）。

表5 無紅邊波段和加入紅邊波段處理下篩選出的反演變量Table 5 The screened variables for retrieval under the No Red-Edge and the Red-Edge Added treatment

從單波段地表反射率來看，兩組處理均包括B2波段的地表反射率；植被指數方面，No Red-Edge處理篩選出了NDVI、SAVI，而Red-Edge Added處理只有MTCI；紋理方面，No Red-Edge處理的紋理主成分包括了第1、2、4、5、9主成分，結合表3可知，紋理指標以異質性（D）、均勻性（H）、能量（ENG）以及信息熵（ENT）為主，波段集中在B3、B4；Red-Edge Added處理的紋理成分包括了第1、2、3、5、6、7、9、11主成分，結合表4可知，加入紅邊波段后，B5波段的紋理特征出現，且以異質性（D）、均勻性（H）和信息熵（ENT）為主要指標。

2.3 模型構建及其反演結果

2.3.1 多元線性模型

利用75%的訓練樣本，根據無紅邊波段（No Red-Edge）、加入紅邊波段（Red-Edge Added）篩選出的自變量，分別與針葉林蓄積量構建多元線性模型，得到兩組處理的多元線性回歸方程為

式（3）、式（4）中，Vno、Vadd分別為No Red-Edge、Red-Edge Added處理反演的蓄積量（m3·hm-2）。

根據訓練樣本建立的多元線性回歸模型，兩組處理的顯著性均小于0.001，方程的決定系數（R2）均大于0.55，且Red-Edge Added處理的決定系數有所提高。

2.3.2 隨機森林模型

利用75%的訓練樣本，根據無紅邊波段（No Red-Edge）、加入紅邊波段（Red-Edge Added）兩組處理篩選出的自變量，分別與針葉林蓄積量構建隨機森林模型。關鍵的模型參數主要包括決策樹個數、決策樹最大深度、決策樹的最大特征數、葉子節點的最小樣本個數和最大葉子節點數5個參數。

以決定系數為代價函數，經過多次試驗，得到各參數的最優值（表6）。結果表明對于訓練樣本，兩組處理的R2均在0.85以上，且Red-Edge Added處理的決定系數略高于No Red-Edge處理。

表6 無紅邊波段和加入紅邊波段處理下隨機森林模型參數優化結果Table 6 The parameter optimization results of rand forest model under the No Red-Edge and the Red-Edge Added treatment

2.3.3 模型驗證

利用25%的預測樣本，將No Red-Edge、Red-Edge Added兩組處理篩選出的自變量，代入到構建的MLR、RF模型中，得到兩組處理各模型的蓄積量反演值。將十折交叉驗證結果平均，得到西寧市針葉林蓄積量的遙感反演結果與觀測值之間的關系圖（圖2）。由圖2可見，兩組處理的MLR、RF模型反演的蓄積量結果與觀測值之間的相關性，其顯著性均小于0.001；從散點圖的分布來看，MLR模型的反演結果較離散，而RF模型的反演結果更集中；且Red-Edge Added處理RF模型反演結果與觀測值更接近1:1線，其精確度達到0.9186。另外，MLR、RF模型反演的蓄積量在較低水平（小于50m3·hm-2）的均高于觀測值，而在較高水平（大于100m3·hm-2）的均低于觀測值。

圖2 無紅邊波段(a)和加入紅邊波段(b)后預測樣本在多元線性回歸(1)和隨機森林(2)模型蓄積量的反演結果Fig.2 The forest stock volume retrieved results of the multiple linear regression model (1) and the random forest regression model(2) based on the predicting samples for the No Red-Edge treatment (a) and the Red-Edge Added treatment (b) respectively

2.3.4 精度分析

由表7可見，相對于No Red-Edge處理，在加入紅邊波段后，MLR模型反演結果與觀測值擬合的決定系數R2表現出一定的提高（由0.5172增大到0.5487），但其均方根誤差幾乎沒有發生變化（由26.4m3·hm-2到26.3m3·hm-2），表明MLR模型不具備檢測GF6-WFV的紅邊波段對西寧市針葉林蓄積量反演影響的能力。而對于RF模型，在加入紅邊波段之后，其決定系數R2則表現出了明顯的增高（由0.6120增大到0.6719），且均方根誤差降低了10.3%（由23.2m3·hm-2降至20.8m3·hm-2）。對表7中的數據進一步分析可知，與MLR模型相比，RF模型在No Red-Edge、Red-Edge Added處理的決定系數R2分別提高了18.3%、22.5%，均方根誤差分別降低了12.1%、20.9%。表明RF模型可以更好地檢測出紅邊波段對蓄積量反演影響，在去除模型對反演精度的影響后，相對于No Red-Edge處理，Red-Edge Added處理的反演結果與觀測值擬合的決定系數R2提高了11.6%，均方根誤差降低了9.1%。

表7 無紅邊波段和加入紅邊波段處理下不同模型的森林蓄積量反演精度結果對比Table 7 Comparison of the retrieved precision of the forest volume for different models under the No Red-Edge and the Red-Edge Added treatment

3 結論與討論

3.1 討論

（1）遙感影像的紋理特征增強了森林冠層與灌草的差別，因而在森林蓄積量的遙感反演研究中得到了廣泛的應用[32-33]。有研究指出紋理特征是物體的空間結構特征的顯示，與影像的光譜特性無關，因此可以任意選擇一個波段[34]，王月婷等[4,35]的研究就采用了Landsat8的全色波段。但也有研究表明不同波段的紋理特征對森林蓄積量的影響不同，劉俊等[32]研究得到近紅外波段與柞樹蓄積量具有較強的相關性，曹霖等[18]也將不同波段的紋理指標作為不同的紋理特征變量。本研究No Red-Edge處理中選取了GF6-WFV前4個波段的紋理指標，分析處理后，篩選出的主成分信息集中在紅波段與近紅外波段，這與張沁雨等[24,32]的研究結論一致。同時本研究在Red-Edge Added處理中選取了GF6-WFV前6個波段的紋理指標，且在篩選出的主成分信息中捕捉到了紅邊波段的信息，進一步表明在森林蓄積量的遙感反演研究中，需要考慮不同波段的紋理特征。

（2）本研究采用逐步回歸法，篩選出了No Red-Edge、Red-Edge Added兩組處理的自變量。對于Red-Edge Added處理，紅邊1波段（B5）的地表反射率入選，這與Hu等[16]的研究結論一致。Red-Edge Added處理中，植被指數只有MTCI指數入選，而MTCI指數主要用于反映植被冠層的綠素含量，與NDVI相比，MTCI指數減少了大氣、土壤背景的影響[36]，削弱反射率的“碗邊效應”[37]，更適合遙感定量反演與植被生長的參數[38-39]。另外，本研究中的NDREI指數未入選，從計算公式來看，盡管MTCI與NDREI中的兩個紅邊波段均參與計算，但MTCI涉及紅波段，而NDREI只是兩個紅邊波段的運算，進一步表明即使加入紅邊波段，紅波段對蓄積量的貢獻不可忽視，曹霖等[18]研究證實了這一點。

（3）本研究采用不同模型對有無紅邊波段處理的蓄積量進行反演，結果表明加入紅邊波段以后，針葉林蓄積量的反演精度得到了提高，均方根誤差降低了9.1%，這與曹霖等[18]在2018年基于Sentinel-2對吉林省中東部的森林蓄積量的反演結果一致。與此同時，本研究發現蓄積量在較低水平（小于50m3·hm-2）的反演結果均高于觀測值，而在較高水平（大于100m3·hm-2）的反演結果卻低于觀測值，這在曹霖等[18]的研究中同樣有此“飽和”現象，曹霖等[18]解釋這可能與植被指數的“飽和”現象有關。

（4）本研究引入紅邊波段后提高了針葉林蓄積量的反演精度，但與觀測值相比，還存在著諸多不確定性及誤差，主要表現在：①數據在時間上的匹配度：本研究忽略了蓄積量獲取時間與影像數據獲取時間的滯后性，盡管針葉樹種在一個齡級內生長緩慢，但在初期（蓄積量較低水平）由于生長速率較快，可能會出現反演結果高于觀測值，而在后期（蓄積量較高水平）生長速率降低，反演結果低于觀測值的情況；②蓄積量的觀測誤差：本研究的蓄積量數據來源于西寧市2014年森林資源二類調查數據，盡管對蓄積量的原始數據進行了數據清洗及預處理，但也會存在因調查小班數目多、分布廣、參與觀測人員多等原因而導致的觀測誤差；③遙感影像的選擇：本研究采用光學遙感對蓄積量進行反演，但光學遙感并不能穿透樹冠，缺乏更多森林內部結構等信息，進而影響蓄積量的反演精度；④模型的選擇：本研究只選取了多元線性模型與隨機森林模型來反演西寧市針葉林的蓄積量，缺乏更多模型（如支持向量機、KNN模型、BP神經網絡等）的對比印證。因此，在后續蓄積量的遙感反演研究中，需要引入多源數據包括光學、高光譜、雷達遙感以及LiDAR等，結合同期的蓄積量實測數據，兼顧反演模型的適用性，從森林蓄積量的形成機理出發，在宏觀與微觀尺度上對森林蓄積量的遙感反演作綜合研究。本研究揭示了紅邊波段在針葉林蓄積量的遙感反演方面的作用，為森林蓄積量的遙感反演及紅邊波段的應用提供了更多的科學參考。

3.2 結論

（1）不同波段的紋理特征與針葉林蓄積量的相關性不同。No Red-Edge與Red-Edge Added兩組處理的紋理特征變量降維后，其主成分信息主要解釋了紅、近紅外、紅邊1波段的紋理特征。

（2）與蓄積量相關的光譜特征變量No Red-Edge處理主要包括紅、近紅外波段的地表反射率，而Red-Edge Added處理主要為紅邊1波段的地表反射率；與蓄積量相關的植被指數變量No Red-Edge處理主要包括NDVI與SAVI，而Red-Edge Added處理為MTCI。

（3）在去除模型對反演精度的影響之后，相對于No Red-Edge處理，Red-Edge Added處理的反演結果與觀測值擬合的決定系數R2提高了11.6%，均方根誤差降低了9.1%，表明WF6-WFV的紅邊波段提高了西寧市針葉林的蓄積量反演精度。