林文科，陸亞剛，蔣先蝶，李桂英，李登秋，陸燈盛

1.福建師范大學 濕潤亞熱帶山地生態(tài)國家重點實驗室培育基地,福州 350007;2.福建師范大學 地理研究所,福州 350007;3.國家林業(yè)和草原局華東調(diào)查規(guī)劃設計院,杭州 310000;4.福建省減災中心,福州 350001

1 引 言

森林是減緩氣候變化和提供多種生態(tài)服務功能的重要陸地生態(tài)系統(tǒng)，是國家重要的戰(zhàn)略資源。森林蓄積量是反映森林資源質(zhì)量狀況和實現(xiàn)2060碳中和的重要基礎數(shù)據(jù)，精準估測森林蓄積量也是國家制定科學合理的碳中和方案的迫切需求。亞熱帶地區(qū)具有獨特的自然特征和森林分布特征，其復雜的地表增加了樣地調(diào)查的難度和蓄積量估測的不確定性。而中國具有廣袤的亞熱帶森林分布，如何提高亞熱帶地區(qū)森林蓄積量估測精度，實現(xiàn)快速、準確的森林資源變化監(jiān)測，對實現(xiàn)森林資源可持續(xù)發(fā)展具有重要意義。

遙感手段估測森林蓄積量主要通過光譜信息或垂直結構信息與樣地實測數(shù)據(jù)建立估測模型（李德仁 等，2012；Dos Reis 等，2019）。光學遙感數(shù)據(jù)，特別是Landsat 系列，是森林蓄積量估測最常用的遙感數(shù)據(jù)源（Barrett 等，2016；Babcock等，2018）。Sentinel-2 數(shù)據(jù)比Landsat 數(shù)據(jù)具有更高的空間、光譜和時間分辨率，在森林蓄積量估測中更具優(yōu)勢（Chrysafis 等，2017；Mura 等，2018）。 然而，包括Sentinel-2 在內(nèi)的各種光學遙感都存在著穿透性差，只能記錄植被水平結構特征，難以反映林分垂直結構信息等問題；另外，光學遙感數(shù)據(jù)的飽和問題，特別是在森林茂密的情況下，導致森林蓄積量的嚴重低估，限制了實際中的廣泛應用（Zhao 等，2016；Nuthammachot等，2022）。

激光雷達數(shù)據(jù)能獲得森林垂直結構信息，可以有效的提高蓄積量估測精度，已成為森林資源估測研究的重要數(shù)據(jù)源（李旺等，2015；Jayathunga 等，2018），但數(shù)據(jù)獲取成本昂貴，覆蓋面積小，難以進行大范圍的森林蓄積量反演研究。衛(wèi)星立體成像數(shù)據(jù)覆蓋范圍廣，重返周期短，數(shù)據(jù)成本低，可獲取大范圍的地面高程信息。Li等（2019）利用生長季和落葉季的資源三號（ZY-3）立體成像數(shù)據(jù)提取北方落葉松林的冠層高度模型（CHM），進而估測森林地上生物量，顯著地提高了森林地上生物量的估測精度。但這種方法僅僅適用于地形平緩的落葉林區(qū)，并不適合地形復雜的山區(qū)及亞熱帶常綠闊葉林區(qū)。

森林蓄積量估測通常采用傳統(tǒng)的線性回歸或機器學習算法（羅環(huán)敏等，2011；Lu等，2016）。機器學習算法能有效地擬合復雜的非線性關系，在森林蓄積量估測研究中得到廣泛的應用（蒙詩櫟等，2017；Feng 等，2017；Gao 等，2018）。但機器學習方法中的參數(shù)優(yōu)化和模型過擬合等問題，常導致預測結果精度低。貝葉斯分層模型以低樣本需求和良好的防過擬合能力成為當前研究的熱點（Junttila和Laine，2017；Ver Planck 等，2018）。貝葉斯方法能夠?qū)崿F(xiàn)復雜的、具有層次結構的參數(shù)關聯(lián)的有效估計，可以在較小訓練樣本集和高度相關的多維數(shù)據(jù)上取得良好的效果（Junttila等，2015）。

為解決基于單一數(shù)據(jù)集、有限的建模樣本、亞熱帶地區(qū)復雜的地形和多樣的森林類型導致森林蓄積量模型估測精度低的問題，本研究將聯(lián)合Sentinel-2 多光譜數(shù)據(jù)，ZY-3 立體像對和機載Lidar 數(shù)據(jù)，結合森林類型和地形等因素，采用貝葉斯分層方法構建森林蓄積量估測模型，解決建模對樣本量的需求以及亞熱帶地區(qū)樣本采集困難的矛盾，提高森林蓄積量估測制圖精度，為大范圍的森林資源遙感調(diào)查提供關鍵技術支撐。

2 材料與方法

2.1 研究區(qū)

金寨縣位于安徽省西部（圖1），屬于北亞熱帶濕潤季風區(qū)，年均降水量1500 mm 左右，梅雨顯著、夏雨集中。夏季平均氣溫在22℃以上，冬季平均氣溫在10℃以下。該縣總面積為3814 km2，居安徽省各縣之首。全縣平均海拔500 m，大別山山脈由西南向東北貫穿全境，形成相對高差1600余米的垂直地勢分布特征，使得金寨縣南北的物候相差半月左右。全縣森林資源豐富，森林覆蓋率高達72.75% （http：//www.ahjinzhai.gov.cn/［2021-08-09］）。典型植被屬亞熱帶常綠闊葉林，但由于溫度、降水量和光照等環(huán)境因子受海拔高度的影響，植被類型具有明顯的垂直地帶性，形成了杉木、馬尾松、板栗、竹類、油桐等為主的森林植被分布特征。

圖1 研究區(qū)位置Fig.1 Location of the study area

2.2 研究數(shù)據(jù)及技術路線

本研究利用了Sentinel-2多光譜、ZY-3立體成像、機載Lidar 和樣地調(diào)查等多種數(shù)據(jù)（表1）。為了評估不同遙感數(shù)據(jù)在森林蓄積量估測中的作用，制定了圖2的技術路線。主要步驟包括（1）數(shù)據(jù)的預處理：對機載Lidar、ZY-3、Sentinel-2 和外業(yè)樣地調(diào)查數(shù)據(jù)的收集和處理，分別得到數(shù)字地形模型（DTM）、數(shù)字地面模型（DSM）、多光譜波段和樣地蓄積量；（2）變量的提取：基于Lidar-DTM 和ZY-3 DSM 生成的林分冠層高度模型數(shù)據(jù)（CHM），Sentinel-2 多光譜波段以及各種植被指數(shù)；（3）變量的提取和篩選：利用統(tǒng)計方法提取基于CHM 和Sentinel-2 多光譜的變量，并將二者合并作為多源數(shù)據(jù)變量，形成3個數(shù)據(jù)變量集，然后分別使用逐步回歸方法篩選變量，篩選的變量分為整體變量（對所有樣本建模選中的變量）和分層變量（分層樣本單獨建模選中的變量）；（4）貝葉斯分層建模：基于3種變量數(shù)據(jù)集篩選出來的變量分別利用雙因素貝葉斯分層方法進行蓄積量估測建模，這里的雙因素是指坡向和森林類型兩個分層因素；（5）結果分析：比較分析3種數(shù)據(jù)方案的建模效果，選出最優(yōu)的數(shù)據(jù)方案結合森林分布圖反演整個研究區(qū)的蓄積量空間分布。

圖2 基于多源遙感數(shù)據(jù)的森林蓄積量估測技術路線Fig.2 Framework of developing forest growing stock volume using multi-sensor remotely sensed data

表1 研究所需數(shù)據(jù)Table 1 Datasets used in research

2.3 數(shù)據(jù)收集與整理

以研究區(qū)的森林類型和年齡為分層依據(jù)，運用分層抽樣法布設樣地。樣地大小均為25.82 m ×25.82 m（1 畝）。樣地調(diào)查于2019-09—2019-10 進行，共調(diào)查樣地71個，包括闊葉林（麻櫟、板栗、楊樹、楓香等）29個、馬尾松26個、和杉木16個。野外作業(yè)時，對樣地進行每木檢尺，起測胸徑為5 cm，詳細記錄單木的樹種、樹高、胸徑等因子，并記錄樣地的優(yōu)勢樹種、平均樹高、齡組、郁閉度、坡度、坡向等因子。依據(jù)樹種、胸徑查找安徽省立木材積表計算單木材積（安徽省地方標準：杉木一元立木材積表（DB 34/T 1724—2012），馬尾松立木材積表（DB34/T 3345-2019），硬闊立木材積表（DB34/T 3907-2021），楊樹立木材積表（DB34/T 3620-2020）），對樣地內(nèi)的單木材積累加得到樣地蓄積量（m3/畝），進而換算為立方米/公頃（m3/ha）。表2為樣地數(shù)據(jù)的統(tǒng)計結果。

表2 外業(yè)調(diào)查樣地數(shù)據(jù)統(tǒng)計Table 2 Statistic data of collected sample plots

兩景Sentinel-2 L1C 級多光譜數(shù)據(jù)由歐洲空間局網(wǎng)站下載，獲取時間為2020-04-09。使用Sen2cor 插件進行輻射定標和大氣校正，將L1C 級產(chǎn)品處理成L2A 級（地表反射率數(shù)據(jù)）。根據(jù)研究目的，去除3 個60 m 波段（氣溶膠、水蒸氣、短波卷云），使用SNAP 融合工具將4 個10 m 和6 個20 m 分辨率的波段融合成10 m 分辨率（Brodu，2017）。為了消除地形對光譜反射率的影響，使用SCS+C 校正模型和機載Lidar 生成的2 m 空間分辨率的DTM 數(shù)據(jù)對融合后的Sentinel-2影像進行地形校正（高永年和張萬昌，2008），然后計算植被指數(shù)，如歸一化植被指數(shù)（NDVI），比值植被指數(shù)（RVI）、土壤調(diào)整植被指數(shù)（SAVI）、增強型植被指數(shù)（EVI））。矢量格式的金寨縣森林地圖提取自2019年安徽省林地“一張圖”?！耙粡垐D”中的每個多邊形都標有詳細的土地覆蓋或森林類型。通過疊加Sentinel-2 多光譜影像、CHM 和土地覆蓋邊界，對森林區(qū)域進行了修改，將土地覆蓋分為4類：闊葉、馬尾松、杉木和其他地類，用于森林蓄積量空間分布的反演，經(jīng)過調(diào)整后的金寨縣森林類型分布如圖1（c）。

2.4 聯(lián)合機載Lidar 和ZY-3 立體成像數(shù)據(jù)提取林分冠層高度模型（CHM）數(shù)據(jù)

機載Lidar 點云數(shù)據(jù)使用RIEGL-VQ-1560i 激光雷達航攝儀獲取，是具有雙通道的機載激光掃描儀系統(tǒng)，系統(tǒng)的內(nèi)部無縫集成高性能IMU/GNSS系統(tǒng)以及一億像素的數(shù)碼相機。點云獲取時間是2019-06，密度為2點/m2。對點云數(shù)據(jù)進行去噪處理后，將離散點云回波點分成地面點和非地面點兩部分，利用不規(guī)則三角網(wǎng)插值算法對地面點進行插值運算生成數(shù)字地形模型（DTM），以2 m 的空間分辨率導出研究區(qū)的DTM數(shù)據(jù)。

ZY-3 立體成像數(shù)據(jù)采集時間為2020-04-09。從立體成像數(shù)據(jù)提取DSM 的主要步驟包括輸入數(shù)據(jù)、選擇地面控制點和連接點、計算模型、生成核線和導出DSM（Xie 等，2019）?？刂泣c所用的高程信息以Lidar DTM 高程數(shù)據(jù)為參考。根據(jù)立體成像原理提取DSM 的過程分為4 步：（1）相對定向：確定兩個不同影像的相對位置，獲得成像瞬時兩影像之間的位置關系和姿態(tài)參數(shù)；（2）絕對定向：通過地面選取的控制點（GCP）確定立體模型在地面測量坐標系中的位置；（3）連接點：通過連接左右影像建立對應關系生成核線影像；（4）輸出DSM：根據(jù)核線影像提取DSM，以2 m空間分辨率輸出。經(jīng)分析發(fā)現(xiàn)正視和后視影像生成的DSM比正視和前視獲取的DSM 噪聲少，精度高。因此，將前者生成的DSM和Lidar數(shù)據(jù)提取的DTM進行差值運算得到空間分辨率為2 m 的CHM。對該CHM進行后處理，檢查高于50 m 和低于0 m 數(shù)值的像元，對照Lidar 影像確定是否對該像元進行濾波處理或替代處理。

2.5 蓄積量建模變量的提取和篩選

鑒于外業(yè)調(diào)查樣地的大小為25.82 m×25.82 m，針對10 m空間分辨率的Sentinel-2光譜波段及植被指數(shù)數(shù)據(jù)，以樣地中心點提取3×3 個像素的統(tǒng)計值作為新的變量賦予對應的樣地，統(tǒng)計值包括最大值、最小值、平均值、標準差和方差。對于CHM 數(shù)據(jù)，則以13×13 個像素為統(tǒng)計單元，提取樣地范圍的最大值、平均值、標準差、方差和高度百分位（10th，20th，…，90th）作為變量。

對于這些以樣地范圍提取的統(tǒng)計變量，需要通過相關性分析進一步篩選出用于森林蓄積量建模的變量。通過將訓練樣本的蓄積量值與樣本所在范圍提取的變量建立逐步回歸方程，將被選入方程中的變量作為建模變量。具體方法為，對于3種數(shù)據(jù)（CHM，Sentinel-2 和多源數(shù)據(jù)），利用逐步回歸方法對3種數(shù)據(jù)集分別重復以下操作：（1）利用所有訓練樣本建立逐步回歸方程提取整體變量；（2）將訓練樣本按森林類型進行分組，每個組的樣本各自建立逐步回歸方程提取分層變量；（3）將整體變量和分層變量合并作為貝葉斯分層建模變量。

2.6 基于貝葉斯分層方法的森林蓄積量估測模型

貝葉斯分層模型是考慮了先驗分布的多層模型，允許同時對不同層次的數(shù)據(jù)進行建模，從而有效地考慮數(shù)據(jù)的嵌套和依賴關系。多層模型的核心是通過對線性組合預測因子η進行變換來預測響應變量y，并假設y服從于一個通過反函數(shù)f假設確定的分布D，yi~D(f(ηi)，θ)，D是對模型中誤差分布和反函數(shù)的描述，參數(shù)θ表示不隨數(shù)據(jù)點變化的特定參數(shù)，i表示數(shù)據(jù)點的個數(shù)。線性組合的預測因子一般表達為η=Xβ+Zμ，式中，β表示整體水平上的系數(shù)，μ表示分層水平上的系數(shù)，β和μ組成模型的參數(shù)集，分別代表固定影響和隨機影響。X和Z為對應于系數(shù)β和μ組成的變量矩陣。因此y、X和Z組成了模型的變量數(shù)據(jù)集。獲得整體水平回歸參數(shù)β需要總體樣本的先驗分布，但先驗分布需要前期大量數(shù)據(jù)的支撐，因此假定整體水平的參數(shù)不受正態(tài)先驗的限制。分層水平的回歸參數(shù)μ假設來自均值為零且協(xié)方差矩陣Σ未知的多元正態(tài)分布：μ~N(0，Σ），當假設不同分組參數(shù)之間的協(xié)方差為零時，矩陣可以分解為：μk~N(0，Σk），k指示分組的因素。并假設同組別內(nèi)不同級別（由j索引）的參數(shù)是獨立的，那么可以得到：μkj~N(0，Vk），將協(xié)方差矩陣Vk作為參數(shù)進行建模，令：Vk=D(σk)ΩkD(σk)，這里D(σk)表示具有對角元素的對角矩陣σk，Ωk代表Vk參數(shù)化后的相關矩陣。為Ωk指定參數(shù)是通過LKJ-相關先驗參數(shù)ξ> 0 來定義先驗（Lewandowski 等，2009），Ωk~LKJ(ξ)，這種先驗通常比二分之一的柯西先驗更能使模型收斂，但仍然是相對薄弱的信息。通過LKJ模型來得到Ωk，并最終可以求出未知協(xié)方差矩陣Σk，從而得到分層水平的回歸參數(shù)μ。

上述過程由R 程序中的brms包完成，而brms通過調(diào)用Stan程序完成計算。brms通過擴展的lme4 語法提供了高效和可讀的Stan代碼，通過R調(diào)用Stan程序在后臺運行，運行結果以可視化的形式再返回R 中，關于brms包中模型的一般性概述可參考（Bates等，2015；Bürkner，2017）。brms中應用的公式語法建立在lme4 的語法之上，語法的基本模式與lme4 一致，但在分層和參數(shù)設置上具有更大的靈活性。lme4 基本語法的形式為response~pterms+(gterms|group)，response代表響應變量，pterms部分表示通過觀測這部分的變量對整體水平的影響是相同的，gterms部分包含了所謂的組級效應，這些效應假定隨著組中指定的組變量而變化，group表示分組的策略。而brms在此基礎上進一步擴展使其能使用多組級效應模型和多成員模型。

本研究同時考慮了兩種分組因子—森林類型（闊葉林、馬尾松和杉木）和坡向（陽坡、半陽坡、半陰坡和陰坡）。然后聯(lián)合兩種分組因子將樣本分成12 組建立模型，命名為雙因素貝葉斯分層模型。在貝葉斯分層模型中，模型變量由兩部分組成：整體變量和分層變量。根據(jù)先驗知識將一些具有整體效應的變量作為整體變量。為了獲得整體和分層變量，利用所有的樣本和變量進行逐步回歸建模，篩選出整體變量，將樣本依據(jù)分組進行逐步回歸建模，篩選出每個分層的變量，最后的變量組成用于貝葉斯分層建模的變量數(shù)據(jù)集。

2.7 森林蓄積量估測結果評價

基于樣本蓄積量和遙感變量建立的蓄積量估測模型需要進行精度評價。本文采用決定系數(shù)R2評價模型回歸系數(shù)擬合優(yōu)度，采用均方根誤差（RMSE）和相對均方根誤差（rRMSE）評價模型預測能力。R2越高，模型擬合效果越好；而RMSE和rRMSE 越低，模型的預測效果越好。由于樣本數(shù)量的限制，本研究采用留一交叉驗證法，也就是說，對于71 個樣本，保留其中一個單獨的樣本作為驗證，其余70個樣本作為訓練數(shù)據(jù)構建模型，重復操作直至每個樣本被驗證。每一次驗證，訓練數(shù)據(jù)幾乎使用了所有的樣本，可以保證產(chǎn)生可靠的預測，避免了隨機因素的影響。通過比較基于不同數(shù)據(jù)集建立的貝葉斯分層模型的性能，選擇最佳的數(shù)據(jù)方案，并將其應用于整個研究區(qū)的森林蓄積量預測。

3 結果分析

3.1 森林蓄積量估測模型分析

基于CHM、Sentinel-2 以及多源遙感數(shù)據(jù)3 種方案的建模變量篩選和模型精度驗證結果（表3）表明，多源遙感數(shù)據(jù)的結合具有較好的建模效果。多源數(shù)據(jù)集在整體變量和分層變量都選中了高度和光譜變量，兩者共同在模型中發(fā)揮作用，模型的決定系數(shù)R2從單一數(shù)據(jù)Sentinel-2 的0.81 和CHM 的0.83 提高到了多源數(shù)據(jù)的0.93。RMSE 和rRMSE 數(shù)據(jù)中也體現(xiàn)了多源數(shù)據(jù)在改善森林蓄積量估測精度的作用?；贑HM 數(shù)據(jù)的蓄積量估測誤差比Sentinel-2數(shù)據(jù)的估測誤差降低了4.3%，二者的結合再降低了5.2%。從森林蓄積量反演結果（圖3）上看，CHM 反演的結果主要集中在50—150 m3/ha，而Sentinel-2 數(shù)據(jù)反演的結果集中在50—200 m3/ha，更多的出現(xiàn)在150—200 m3/ha，很少出現(xiàn)低于50 m3/ha的預測?；诙嘣磾?shù)據(jù)的反演結果則出現(xiàn)較多低于50 m3/ha 的預測，以及高于250 m3/ha 的預測結果，蓄積量較高的地方主要集中在西南部。

表3 貝葉斯建模變量及建模精度Table 3 A summary of selected variables for modeling and evaluation results of modeling performance

圖3 不同數(shù)據(jù)方案反演的部分研究區(qū)森林蓄積量空間分布圖Fig.3 Spatial distribution of forest growing stock volume from different data scenarios

3.2 貝葉斯分層模型估算結果分析

基于3種數(shù)據(jù)集對不同森林類型的蓄積量估測的精度（表4）分析表明，3 種數(shù)據(jù)集對所有森林類型估測效果呈現(xiàn)相同的趨勢，估測精度從高到低都依次是多源數(shù)據(jù)、CHM 和Sentinel-2 數(shù)據(jù)。其中，多源數(shù)據(jù)對闊葉林蓄積量估測精度的作用最明顯，比CHM 數(shù)據(jù)的估測精度提高了6.7%，比Sentinel-2 數(shù)據(jù)精度提高了12.7%。3 種數(shù)據(jù)集的蓄積量預測值與實際值之間的關系圖和殘差（圖4）可以看出，對于多源數(shù)據(jù)的森林蓄積量估測，驗證樣本點多靠近1∶1 線上，殘差也從單一數(shù)據(jù)時（Sentinel-2 或CHM）的0—100 m3/ha 變?yōu)?—75 m3/ha，且殘差并沒有隨蓄積量的增大而出現(xiàn)明顯波動；然而，使用單一數(shù)據(jù)都不可避免的出現(xiàn)殘差值隨實測蓄積量的增大而增大。此外，比較兩種單一數(shù)據(jù)發(fā)現(xiàn)，總體上Sentinel-2 數(shù)據(jù)的預測結果不如CHM 數(shù)據(jù)，特別是對于闊葉林的預測和低蓄積量值的預測。

表4 根據(jù)不同森林類型評估蓄積量建模結果Table 4 Evaluation of forest growing stock volume modeling results according to forest types

圖4 樣地蓄積量觀測值與估測值之間的關系圖以及二者之間的殘差圖Fig.4 Relationships between the estimated forest growing stock volume and reference data and residual between estimates and reference data

不同坡向的森林蓄積量估測結果（表5）表明，CHM 數(shù)據(jù)集和多源數(shù)據(jù)集對蓄積量估測的精度受坡向的影響。在朝陽面（陽坡、半陽坡）兩個數(shù)據(jù)集的估測誤差都小于背陰面（半陰坡、陰坡），其中CHM 數(shù)據(jù)陽坡比陰坡精度高出5.5%，多源數(shù)據(jù)則高出6.9%。這與ZY-3 立體像對提取DSM 時不同坡向的精度差異有關（Rahlf 等，2014）。而在Sentinel-2 數(shù)據(jù)集中則沒有這種規(guī)律，陽坡比陰坡僅僅高出0.3%，這是因為前期對Sentinel-2 影像做過地形校正，很大程度上消除了不同坡向引起的光譜值差異。3 種數(shù)據(jù)的估測結果表明對于包含CHM 的數(shù)據(jù)集，進行坡向分層是必要的；而在Sentinel-2 數(shù)據(jù)中，除了半陽坡，其他各個坡向之間的估測精度沒有明顯差異。此外，將兩種數(shù)據(jù)聯(lián)合后，各個坡向的估測誤差都顯著下降，且在坡向上體現(xiàn)出差異，說明多源數(shù)據(jù)進行坡向分層建模有利于提高森林蓄積量的估測效果。

表5 基于不同坡向的蓄積量估測結果Table 5 Evaluation of accuracy assessment results according to different slope aspects

分析3 種數(shù)據(jù)集對不同森林蓄積量等級的預測效果（表6）可以看出：除區(qū)間200—250 m3/ha外，CHM 的估測誤差都小于Sentinel-2 數(shù)據(jù)，但兩種數(shù)據(jù)都不可避免的在高值區(qū)域出現(xiàn)了較高的估測誤差；多源數(shù)據(jù)則展現(xiàn)出良好的估測結果，除100—200 m3/ha 的區(qū)間外，多源數(shù)據(jù)都有最低的估測誤差。此外CHM 數(shù)據(jù)和Sentinel-2 數(shù)據(jù)的估測誤差都呈現(xiàn)出兩頭大中間小的特點，但這種現(xiàn)象在多源數(shù)據(jù)中則不明顯，各個蓄積量區(qū)間的RMSE 趨于相等，說明數(shù)據(jù)聯(lián)合有利于改善高值低估和低值高估的現(xiàn)象。比如，在0—50 m3/ha 和50—100 m3/ha 的蓄積量區(qū)間中，多源數(shù)據(jù)的估測誤差比Sentinel-2 數(shù)據(jù)分別降低了18.0 m3/ha 和16.4 m3/ha，在大于250 m3/ha 的蓄積量區(qū)間也降低了近21.8 m3/ha，rRMSE 只有14.1%，估測精度得到大幅度提升。

表6 3種數(shù)據(jù)在不同蓄積量范圍的估測結果精度分析Table 6 Evaluation of results based on three data sets according to forest stock volume ranges

4 討 論

4.1 多源數(shù)據(jù)聯(lián)合能提高森林蓄積量估測精度

Sentinel-2 等多光譜數(shù)據(jù)一般只能提供林分冠層水平信息，無法獲得垂直高度信息。因樹高與森林蓄積量密切相關，基于光譜數(shù)據(jù)的蓄積量估測普遍不如基于高度數(shù)據(jù)的估測精度（Feng 等，2017）。如利用Landsat TM 數(shù)據(jù)估測雪松蓄積量時，當雪松蓄積量高于200 m3/ha 時，飽和問題使得擬合效果大幅度下降（Lee 和Phua，2010），在亞馬遜熱帶雨林的森林生物量和中國亞熱帶森林生物量估測時遇到同樣的問題（Lu，2005；Zhao等，2016；Gao 等，2018），利用其他單一數(shù)據(jù)源，如RapidEye （Feng 等，2017）、SAR 數(shù) 據(jù)（Dos Reis 等，2019）也不可避免的因飽和問題影響估測結果。本研究基于每個樹種的Sentinel-2 數(shù)據(jù)中出現(xiàn)明顯的光譜飽和現(xiàn)象（圖5），但在蓄積量估測結果的散點圖中（圖4（b）），沒有出現(xiàn)明顯的飽和現(xiàn)象。原因在于使用的貝葉斯模型對于蓄積量的估測考慮了森林類型，同時本研究的樣本又具有特殊性，即高森林蓄積量樣本（>250 m3/ha）全部為杉木（圖5），且杉木的樣本也基本分布在高蓄積量區(qū)域，而杉木在250 m3/ha 之前的光譜還未出現(xiàn)光譜飽和，因此高值區(qū)域的估測很容易分布在大于250 m3/ha 的區(qū)域，使得原來在高值區(qū)域的蓄積量飽和現(xiàn)象不明顯。另外，馬尾松和闊葉在蓄積量大于200 m3/ha 已經(jīng)出現(xiàn)光譜飽和，但由于應用了分層建模，這兩種樹種的光譜飽和不會對杉木的估測造成影響，從側(cè)面也說明了分層建模在減少飽和問題上的作用。

圖5 樣地蓄積量與短波紅外波段的反射率之間的關系Fig.5 Relationships between the forest growing stock volume and the surface reflectance of shortwave spectral band

基于樹高變量的森林蓄積量估測，不存在數(shù)據(jù)飽和問題，但樹高變量與蓄積量之間的關系受其他因素，如樹種、年齡、立地條件等的影響。如何從不同維度獲取森林的信息，充分利用多源數(shù)據(jù)各自的優(yōu)勢，通過數(shù)據(jù)互補充分發(fā)揮數(shù)據(jù)潛力來彌補單一數(shù)據(jù)的不足，是提高森林蓄積量估測精度的重要手段（Simard 等，2011；Laurin 等，2018）。如在多光譜數(shù)據(jù)的基礎上加入年齡、地形因子和SAR 數(shù)據(jù)，桉樹蓄積量的估測誤差減少了近1/3（Dos Reis 等，2019）。利用Lidar 的冠層穿透力結合Landsat OLI 數(shù)據(jù)估測亞熱帶森林生物量，比單一數(shù)據(jù)明顯提升了估測效果（曹林等，2016）。本研究進一步證實了以上結果，Sentinel-2數(shù)據(jù)聯(lián)合CHM后，估測精度提高近10%。

4.2 多源數(shù)據(jù)對于高值和低值區(qū)域估測精度的提升

在蓄積量估測研究中，蓄積量的低值高估和高值低估常常是困擾研究者的難題（Lu等，2016；Gao 等，2018）。本研究發(fā)現(xiàn)Sentinel-2 光譜和CHM 數(shù)據(jù)的結合可以有效降低這種現(xiàn)象對估測精度的影響。本研究中總樣本的平均蓄積量值為141.6 m3/ha，而不同蓄積量范圍的估測精度驗證發(fā)現(xiàn)，當蓄積量位于100—200 m3/ha時，多源數(shù)據(jù)的估測精度不如CHM，特別是100—150 m3/ha 時，反而降低了5.5%，但當蓄積量偏離平均值附近時，即在0—100 m3/ha 和大于200 m3/ha 時，多源數(shù)據(jù)的估測誤差都小于其他兩種數(shù)據(jù)，說明數(shù)據(jù)聯(lián)合對于模型的改進主要集中在高值和低值區(qū)域，而對平均值附近的改進不明顯。因為當響應變量和預測變量相關性不高時，通過預測變量得到的響應變量的估測值更容易趨近于平均值，這在3種數(shù)據(jù)方案結果中都明顯存在，因此平均值附近的預測誤差往往最低。雖然在平均值附近多源數(shù)據(jù)的精度有所下降，但不影響整體精度的提高，因為蓄積量估測的問題主要出現(xiàn)高值區(qū)和低值區(qū)，而數(shù)據(jù)聯(lián)合方法在這兩個區(qū)域很好地實現(xiàn)了估測誤差的下降，說明通過數(shù)據(jù)聯(lián)合來提高蓄積量估測精度是切實可行的方法（Li等，2019）。

4.3 分層貝葉斯建模的依據(jù)和優(yōu)勢

研究中基于森林類型分層是基于不同森林類型的蓄積量和樹高的關系不同。另外在森林類型分層的基礎上，發(fā)現(xiàn)CHM 和蓄積量的關系與坡向相關，這種相關性可能來自3方面的原因，（1）樹木生長受坡向的影響；（2）資源三號DSM 與機載Lidar DTM 之間位置有錯位，使得二者相減沿錯位移動方向產(chǎn)生坡向相關特征；（3）光照條件隨坡向發(fā)生變化，導致點云匹配出現(xiàn)差異。本文的前期實驗中，曾對只使用Sentinel-2 數(shù)據(jù)的不同建模方法做比較，發(fā)現(xiàn)分坡向的建模結果并不明顯優(yōu)于不分坡向的建模結果，說明這種相關性更多的不在于坡向間樹木生長的差異，可能在于CHM 數(shù)據(jù)在不同坡向上生產(chǎn)的差異。Rahlf 等（2014）發(fā)現(xiàn)基于ZY-3立體像對提取的DSM 高程會受到坡向的影響，背光面往往出現(xiàn)DSM 數(shù)值偏大的現(xiàn)象，而向陽面與實際地表的狀況較為吻合。這種坡向引起的差異，說明分層建模的必要性。然而分層建模意味著更多的樣本需求，本研究只有71 個野外調(diào)查樣地，因此利用貝葉斯分層建模既可以減少不同森林類型和坡向差異對整體樣本的影響，又能有效利用整體樣本，從而提高估測精度，減少低值高估和高值低估問題。從分樹種以及分坡向精度驗證結果也可以發(fā)現(xiàn)，不同樹種和坡向的估測精度之間存在明顯的差異，為進一步提高森林蓄積量估測精度提供了方向。

5 結 論

森林蓄積量遙感估測中林分高度信息比光譜信息更有優(yōu)勢，CHM 數(shù)據(jù)的估測結果明顯優(yōu)于Sentinel-2 光譜數(shù)據(jù)，而通過聯(lián)合光譜信息和CHM數(shù)據(jù)的方法是提高森林蓄積量估測精度的有效手段。通過本研究的實驗結果，得到以下結論：

（1）3 種數(shù)據(jù)集中，多源數(shù)據(jù)在蓄積量估測上具有明顯的優(yōu)勢，其總體估測結果和分層結果評價都優(yōu)于單獨使用Sentinel-2 數(shù)據(jù)集或CHM 數(shù)據(jù)集。基于多源數(shù)據(jù)的森林蓄積量估測誤差比Sentinel-2 數(shù)據(jù)降低了13.6 m3/ha，比CHM 降低了7.4 m3/ha；

（2）多源數(shù)據(jù)可以有效改進高值低估和低值高估問題。多源數(shù)據(jù)在0—50 m3/ha 區(qū)域的估測誤差比單獨使用Sentinel-2 數(shù)據(jù)估測誤差降低約1/4；在大于250 m3/ha 區(qū)域的估測誤差為比單獨使用Sentinel-2 或CHM 數(shù)據(jù)估測誤差也降低了近1/3。因此多源數(shù)據(jù)對于提高森林蓄積量估測精度特別是高值和低值區(qū)域具有顯著作用；

（3）在亞熱帶地區(qū)使用多源遙感數(shù)據(jù)進行森林蓄積量估測的結果表明，當樣本數(shù)不足或數(shù)據(jù)隨某個屬性存在明顯差異時，使用貝葉斯分層建模方法能發(fā)揮不同數(shù)據(jù)的優(yōu)勢，從而提高遙感手段估測森林蓄積量的精度。