基于序貫高斯條件模擬的GEDI數據聯合Landsat8反演森林地上生物量

摘要：[目的]單一遙感技術估測森林生物量存在較大局限性，本研究旨在利用多源遙感協同技術互補激光雷達和光學遙感的優勢，提高生物量估測精度。[方法]以星載激光雷達GEDI和光學遙感Landsat8數據為主要信息源，采用序貫高斯條件模擬方法實現GEDI光斑數據由“點”到“面”的空間擴展，結合地面1 38塊生物量調查樣地，利用隨機森林回歸方法估測云南省香格里拉云冷杉林的地上生物量。[結果]（1）采用序貫高斯條件模擬方法對GEDI光斑點進行空間插值，模擬的12個生物物理指標在空間分布上呈現出隨機性、破碎化的特征，這與森林的空間分布聚集性非常相似，參與建模的9個指標OEC均大于0.90；（2）利用單一Landsat8光學遙感數據和地形因子構建的隨機森林模型精度為：R2=0.82，RMSE=35.51 t·hm-2，P-0.77；Landsat8數據協同星載激光雷達GEDI數據構建的隨機森林模型精度為：R2=0.86，RMSE=32.11 t·hm-2，P=0.80，模型精度明顯提升；（3）利用多源遙感技術估測的香格里拉2019年云冷杉林地上的生物量總量為37 042 605.68 t，平均生物量為123.28 t·hm-2。[結論]基于地統計學的序貫高斯條件模擬方法考慮到研究對象的空間異質性、能克服一定的平滑效應，用于實現激光點由“點”到“面”的空間擴展是可行的。星載激光雷達GEDI與光學遙感Landsat8協同的多源遙感數據可有效填補單一遙感數據源的缺陷，提高森林生物量的估測精度，能為激光雷達聯合光學遙感估測大范圍、全覆蓋的森林生物量提供參考。

關鍵詞：GEDI;Landsat 8；序貫高斯條件模擬；隨機森林；生物量

中圖分類號：X87; S757.2 文獻標識碼：A 文章編號：1001-1498（2024）03-0049-12

森林生物量是森林生態系統長期物質循環與能量流動的結果，它是反映森林生產力和森林生態系統功能的重要參數，在森林的經營、監測與評價中起到重要作用。森林生物量的估測對合理保護、管理和利用森林資源有重要意義。傳統的森林生物量估測主要采取實地測量法，雖然精度較高，但受到時間和空間差異的限制，需要耗費大量的人力和物力，并且對森林有一定的破壞作用。隨著3S技術的迅速發展，結合遙感數據估測生物量成為一種有效的技術手段。目前，森林生物量估測大多是利用野外實測的生物量樣地數據，通過分析樣地數據與光學遙感因子之間的相關關系建立估測方程推廣到區域尺度的生物量分布。然而，光學遙感易受天氣影響，對植被的穿透性較弱，只能獲取森林冠層表面的信息，而且信號容易飽和、數據存在混合像元等問題，導致模型估算精度不高。

激光雷達（LiDAR）能快速獲取目標的空間三維信息、對植被的穿透力強、可實現對森林冠層垂直結構參數的精準監測，在森林參數反演、生物量估測方面具有很大優勢。由于激光雷達是離散采樣，很難提供全覆蓋的數據，盡管它突破了光學遙感在森林應用領域的諸多缺陷，但依然存在數據獲取和處理較難等問題。此外，任何單一的遙感技術都存在局限性，激光雷達技術也不例外，多傳感器協同應用能夠填補激光雷達在林業應用領域存在的不足，因此，激光雷達協同光學遙感以獲取更高精度的森林參數迅速得到諸多研究者的應用與推廣。池弘等利用GLAS星載激光雷達結合Landsat/ETM+數據估測了長白山地區森林地上生物量；邊瑞利用無人機激光雷達結合Landsat數據估測了祁連山國家公園的森林生物量；巨一琳等聯合機載LiDAR數據和多光譜數據估算了根河市大興安嶺生態觀測站寒溫帶天然林的生物量。

全球生態系統動力學調查（GEDI）多波束激光雷達相比于GLAS，采樣密度有較大提高，實現對森林冠層高度、生物量等參數的測量具有重要意義。本研究以典型高寒山區香格里拉市為研究區，在對GEDI波形數據處理的基礎上提取大量的森林參數，利用基于地統計學的序貫高斯條件模擬方法對GEDI光斑點森林參數變量實現由“點”到“面”的空間擴展，獲取全覆蓋的GEDI數據，結合Landsat8-OLI數據和2016年調查的138個野外樣地數據，利用隨機森林方法估算香格里拉2019年云冷杉林地上生物量，評估多源遙感數據協同集成應用的優勢。

1 材料與方法

1.1 研究區概況

香格里拉市位于云南省西北部“三江并流”區域，地理坐標為26°52＇～28°52＇N，99°20＇～100°19＇E，總面積11 613 km2，最高海拔為5545 m，最低海拔為1 503 m，屬于山地寒溫帶季風氣候，全年氣溫偏低，年平均氣溫僅為5.4 ℃，四季交替不明顯，夏秋時節雨水較多，冬春時節較為干旱，年均降水量為268～945 mm。境內生態系統完好，森林資源豐富，森林覆蓋率達74.99%，優勢樹種有云冷杉（Picea asperataMast. and Abies fabri （Mast.） Craib）、高山松（Pinus densata Mast.）、云南松 （Pinusyunnanensis Franch.）、高山櫟（Quercussemicarpifolia Smith）等。研究區地理位置見圖1。

1.2 數據來源與處理

1.2.1 星載激光雷達GEDI數據

全球生態系統動力學調查（GEDI）雷達由美國國家航天航空局（NASA）2018年12月5日發射成功，自2019年4月開始收集全球51.6°N～51.6°S之間基于足印的測量數據，軌道高度約400 km，光斑直徑為25m，沿軌光斑點間隔距離為60 m，相鄰軌道間距約600 m，提供了較高的采樣密度和精度。GEDI包含的數據產品種類豐富，根據數據處理階段的不同分為L1、L2、L3、L4四個級。本研究所使用的GEDI數據為2019年4-12月獲取的L2B產品，它是光斑尺度的生物物理指標，包含陸地表面上每個激光足印點精確的經緯度坐標、高程、海拔、覆蓋率和垂直剖面指標等，數據從NASA官方平臺獲取（https：//search.earthdata.nasa.gov）。

提取37個GEDI光斑數據指標，包括quality_flag、degrade_flag、sensitivity 3個質量指標，各指標的描述見表1。GEDI為波形數據，采集過程中會受到云層和其他環境因素的影響，存在部分地理定位精度差以及信號質量差的光斑，因此根據3個質量指標對初始的3 858個GEDI光斑點進行篩選。quality_flag有0和1兩個值，值為O時表示光斑信號質量差，需剔除；degrade_flag有0和1兩個值，值為1時表示衛星處于降軌狀態，所獲取的光斑點無效，需剔除；Sensitivity是靈敏度參數，閾值為0到1，在陸地上使用的GEDI光斑點靈敏度閾值一般在0.95以上，剔除閾值小于0.95的光斑點。篩選后共有1 302個光斑點，其中80%作為訓練數據集，用于獲取各變量的空間分布結果，剩余的20%作為驗證數據集，用于檢驗空間插值的精度。篩選前后的研究區GEDI光斑點分布見圖2。

1.2.2 Landsat8遙感影像數據

使用分辨率為30m的Landsat8-OLI數據，共需3景影像覆蓋整個研究區，數據獲取時間為2019年12月。遙感影像廣泛被用于森林生物量的估測研究，在前期對影像預處理的基礎上，提取NDVI、 RVI、 DVI、PVI、EVI、ARVI、RDVI、V13、SARVI、MSAVI、SLAVI、VIS234、B53、B345、ND43、ND563等16個植被指數用于云冷杉林地上生物量的估測。

1.2.3 野外調查樣地數據

野外調查數據為2016年香格里拉市森林資源二類調查的角規控制樣地數據，共有138塊云冷杉樣地，每塊樣地面積為1 hm2，包含樣地點坐標、小班地類、優勢樹種、平均胸徑、樹高等參數，利用二元材積表計算樣地的林分蓄積量。由于Landsat8-OLI和GEDI的數據采集時間均為2019年，所以根據《云南省主要樹種材積生長率表》的材積凈生長率（云杉0.85%，冷杉0.88%）推算2019年的林分蓄積量，然后利用林木生物量擴展因子法計算2019年的樣地生物量。蓄積量一生物量轉換模型見式（1），云冷杉生物量擴展因子及木材密度見表2，138塊樣地的2019年生物量統計參數見表3。為與遙感影像像元尺度相匹配，將每hm2樣地生物量數據轉換為900m2的樣地生物量數據，樣地點的空間分布見圖1。

B= V·SVD·BEF（1）

式中，B為森林地上生物量，V為蓄積量，SVD為某一樹種或樹種組的木材密度，BEF為生物量擴展因子。

1.2.4 數字高程模型數據

數字高程模型（DEM）數據為美國國家航空航天局（NASA）和日本宇宙航空研究開發機構（JAXA）共同開發的ASTER GDEMV2數字高程數據產品，分辨率為30 m。該數據用于Landsat8影像的地形輻射校正，以及基于柵格表面算法提取坡度、坡向和高程3個地形因子（圖3）。根據圖3可知，香格里拉海拔較高，總體呈現出西北高東南低的趨勢，海拔落差較大，境內地勢陡峭，多為斜坡、陡坡和峭破，地形結構復雜。

1.3 研究方法

1.3.1 變異函數

變異函數的理論模型可通過球狀模型、指數模型和高斯模型進行擬合，用于描述區域化變量的變化規律，變異函數的模型由模型類型、基臺值（Sill）、變程（Range）和塊金值（Nugget）唯一確定。本研究利用GS+9軟件對GEDI光斑變量進行變異函數分析，采用3個模型進行預測，以模型的決定系數（R2），殘差平方和（RSS）對模型進行評價。在分析前，對所有變量進行歸一化處理以提高模型的收斂速度，即將所有變量的數值調至[0，1]或[-1，1]，計算方法見式（2），并對所有變量的數據結構進行正態分布性檢驗，將不符合正態分布的變量進行開立方變換，符合條件后方可進行變異函數分析。在前期對GEDI光斑變量進行變異函數分析和篩選的基礎上，選擇最優變異函數模型用于實現序貫高斯條件模擬。

Y=y-ymin/ymax-ymin（2）

式中，y為歸一化后的值，y為初始值，ymin為初始數據的最小值，ymax為初始數據的最大值。

1.3.2 序貫高斯條件模擬

序貫高斯條件模擬（SGCS）以變異函數理論為基礎，結合區域化變量，用隨機模擬方法估計研究變量。根據已知數據構造高斯函數，將區域化隨機變量Z（x）的每一個取值看作是符合高斯函數（即正態分布函數）F（x）的一次隨機實現。在每一個模擬位置X處，F（x）是以n個已知數據Z （xi）（i=1，2，…，n）和此前的m-1個模擬值z （xj）（j=1，2…，m-1）為條件的累積條件概率密度函數，在模擬過程通過構建局部累計條件概率分布來實現最終結果。該方法對所預測的空間數據可能的取值結果及其概率進行度量，能在最大程度上正確反映區域化變量的空間波動性，再現真實資源特性空間變異曲線的波動。在大多數情況下，原始數據并不滿足模擬過程的要求，需通過正態變換后再進行模擬，最后通過逆變換將模擬結果轉換為初始值。

1.3.3 隨機森林回歸

隨機森林（RF）是以決策樹為基礎的機器學習集成算法，最早由Breiman等提出。它利用Bootsrap重抽樣方法從訓練數據中隨機抽取一部分樣本，對每個Bootsrap樣本進行決策樹建模，然后組合多棵決策樹的預測，通過投票決策的方式得出最終預測結果。在隨機森林中，每棵決策樹都是獨立并在隨機選擇的子樣本上進行訓練的，從而有效降低過擬合的風險。該方法能夠進行特征重要性分析，相比于神經網絡和支持向量機等其它算法，它在分析變量關系上具有顯著優勢。

1.3.4 精度評估方法

（1）空間插值精度評估

序貫高斯條件模擬是基于Kriging估計的一種不確定性空間插值方法，會產生一定的估測誤差。本研究采用總體估計值一致性（OEC）對GEDI各變量的序貫高斯條件模擬結果進行精度評估，即以GEDI光斑點模擬估計的總體總值（TGs）與光斑點估計的總體總值置信區間中值（TGp）的比值來衡量，OEC越接近于1，表明空間插值精度越高。各指標的計算公式如下：

式中，A為研究區總面積，a為光斑面積，n為光斑點觀測數量，gci為光斑點實測值，m為像元總數，scj為第/個像元的模擬值。OEC代表模擬總值與光斑點估計總值的接近程度，二者相等時，OEC=1。

（2）生物量估測模型精度評估

模型精度的高低直接影響到生物量的估測結果，本研究基于138個云冷杉樣地生物量數據，采用留一交叉驗證法驗證云冷杉林生物量非參數估測模型的精度，以模型的決定系數（R2）、均方根誤差（RMSE）、以及總體估測精度（P）作為精度評價指標。

2 結果與分析

2.1 GEDI光斑點變量篩選及序貫高斯條件模擬

2.1.1 GEDI變量篩選

保留光斑點的quality_flag和degrade_flag兩個指標只有1和0單一固定值，參與構建生物量估測模型的意義不大，因此僅用于光斑點的篩選，最終采用剩余的35個變量進行變異函數分析，各變量的模型擬合精度見圖4。根據變異函數的分析結果，剔除R2小于0.5的變量，從而減小序貫高斯條件模擬方法對GEDI光斑點進行空間插值的誤差，獲取精度較高的空間分布結果。篩選后保留12個GEDI變量，其變異函數的擬合參數見表4。

2.1.2 GEDI變量序貫高斯條件模擬

基于GS+9軟件擬合的12個GEDI變量的最優變異函數，根據模型參數（變程、塊金值、偏基臺值）進行簡單Kriging插值，最后采用ArcGIS 10.5中“高斯地統計模擬”工具完成序貫高斯條件模擬。通過反復實驗，依次設定模擬次數為10次、25次、50次、75次、100次、125次，對比分析后，發現50次模擬后像元均值的方差變化趨于穩定，因此確定模擬的次數為50次。為與GEDI光斑點的面積尺度相匹配，輸出GEDI各變量模擬的柵格數據分辨率為22.15 m，并對經過歸一化處理和正態變換后的變量進行反歸一化處理和逆變換，最終得到真實的模擬結果。序貫高斯條件模擬的流程見圖5，12個GEDI變量模擬的結果見圖6。利用預留20%的實測光斑數據檢驗模擬結果的精度，其總體估計值一致性的檢驗結果見表5。

根據總體估計值一致性檢驗結果，pgap_theta_ error、rg_a2、rg_a6 3個變量的OEC較低，且經過檢驗總體估計值未在置信區間內，因此剔除OEC最低的3個變量，最終保留滿足精度要求的9個變量參與構建云冷杉林生物量估測模型。

2.2 隨機森林建模

2.2.1 利用Landsat8光學數據構建估測模型

利用Python 3.7軟件，基于Landsat8影像的16個植被指數和3個地形因子，采用隨機森林特征重要性分析方法對參與建模的19個特征進行特征重要性貢獻度分析，其中對建模貢獻度最大的變量是V13，貢獻率為12.97%；貢獻率最小的變量是VIS234，貢獻率為1.52%，各特征的建模貢獻度占比見圖7。采用所有特征構建隨機森林回歸模型，模型的R2=0.82、RMSE=35.51 t·hm-2、P=0.77，模型精度擬合的散點圖見圖8。

2.2.2 聯合GEDI和Landsat8數據構建估測模型

考慮到與樣地面積和Landsat8數據的柵格大小相匹配，將基于2.1節GEDI光斑數據9個變量的序貫高斯條件模擬結果重采樣成分辨率為30 m的柵格數據。針對保留的9個GEDI光斑變量、Landsat8-OLI影像的16個植被指數和3個地形因子，采用隨機森林特征重要性分析方法對參與建模的所有特征進行特征重要性貢獻度分析，其中對建模貢獻率最大的變量是V13，貢獻率為12.05%；貢獻率最小的變量是RVI，貢獻率為0.96%，各特征的建模貢獻度占比見圖9。

利用9個GEDI光斑變量聯合Landsat8影像的16個植被指數特征、以及3個地形因子建立隨機森林回歸模型，模型的R2=0.86、RMSE=32.11t·hm-2、P=0.80，兩種遙感數據聯合后構建的模型精度明顯提升。模型精度擬合的散點圖見圖10。

2.2.3 生物量估測模型地形因子貢獻度分析

根據圖7和圖9可知，兩個隨機森林模型中3個地形因子的建模貢獻度占比均較大，尤其以DEM較為顯著，前后兩個模型中的貢獻度分別達到11 .18%和7.21%。香格里拉地處高寒山區，海拔落差較大，云冷杉為喜陰樹種，主要分布在海拔較高的陰坡和半陰坡地帶，生長趨勢嚴重受到海拔高度的影響，垂直地帶差異性強，因此DEM對模型的貢獻度較大。

2.3 多源遙感數據協同的生物量估測

根據兩個隨機森林模型精度的對比分析結果，利用星載激光雷達GEDI數據協同Landast8光學遙感數據構建模型估測香格里拉2019年云冷杉林的地上生物量。估測的2019年云冷杉林地上生物量的總量為37 042 605.68 t，平均生物量為123.28 t·hm-2，生物量空間分布結果見圖11。香格里拉地處高寒山區，平均海拔約為3 500 m，總體呈現出西北高東南低的趨勢。云冷杉林主要集中分布于南部地區，根據圖11可知，高生物量值并非全部集中分布于南部地區，而是分布在海拔3 600 m以上的高寒地區，這是因為冷杉林主要分布于海拔3 500～3 800 m的陰坡、半陰坡地帶，云杉林主要分布于海拔3 200～3 700 m的陰坡、半陰坡以及溝谷地帶。此外，香格里拉人口較多的鄉鎮主要分布在南部的建塘鎮、虎跳峽鎮、金江鎮和三壩鄉，受人文因素和自然因素的影響，云冷杉林地上生物量總體分布不均，區域差異較大。

3 討論

3.1 序貫高斯條件模擬方法對激光點空間擴展的應用分析

地統計學是研究空間變異及格局的有效方法，最初被應用于地質研究領域，隨著統計學的發展，以及森林資源變化監測的重要性日益顯著，該方法也被廣泛應用于林業研究領域。星載激光雷達GEDI只能提供離散的采樣點，不具有成像性，所以研究采用基于地統計學的序貫高斯條件模擬方法實現了GEDI光斑由“點”到“面”的空間擴展，獲取與地面調查樣地對應的變量指標數據。相比于Kriging插值，該方法克服了一定的平滑效應，能夠在最大程度上真實地展現研究對象的空間分布格局及變化。根據1 2個GEDI變量的模擬結果（圖6）可知，每個生物物理指標都呈現出破碎化、隨機性的特征，這與香格里拉森林空間分布聚集性非常相似，且參與建模的9個變量OEC均在0.9以上，精度滿足要求，因此，序貫高斯條件模擬方法可用于激光點的空間擴展，為離散的“點”數據轉向大區域尺度的“面”數據的應用提供了可靠的技術支持。

3.2 云冷杉材積生長率不確定性及驅動因素分析

本研究參考《云南省主要樹種材積生長率表》云杉和冷杉的材積凈生長率推算2019年云冷杉林地上生物量，該表按15個齡級分別計算各齡級的材積凈生長率，其中云杉I齡級的材積凈生長率為8.84%，XV齡級的凈生長率為0.35%；冷杉I齡級的材積凈生長率為14.64%，XV齡級的凈生長率為0.12%，兩個樹種不同齡級的材積凈生長率差距較大，因減少成本和受其他因素限制，統一采用總平均凈生長率，相關研究表明，利用該表給出的凈生長率計算的材積生長量結果偏小。根據2016年二類調查的林木材積推算2019年的林木材積，期間有3年的時間跨度，近年來氣候變化多樣，時間跨度導致的氣候因素會影響林木的生長趨勢，進而影響生物量的生長量，賈勃等的研究發現，生物量的生長量受到年均溫度的正向影響；何瀟等的研究也表明氣候因子對林分生物量有顯著影響。此外，研究利用生物量擴展因子法計算2019年云冷杉林的地上生物量，三重不確定性因素會導致估測結果存在微弱的偏差。胥麗等利用該方法獲取了香格里拉2019年的櫟類森林地上生物量精度較高的估測結果，但后續研究中，若條件允許建議采用實地驗證法對估測結果進行檢驗，并對林木材積生長率變化氣候驅動因素進行分析，降低生長率不確定性因素導致的生物量估測偏差。

3.3 生物量估測精度分析及模型優化建議

針對生物量的估測結果，本研究與前人的結果做了一些對比。謝福明等基于KNN模型估測了香格里拉2016年云冷杉林地上生物量，RMSE為41.64 t·hm-2．而本研究基于隨機森林模型的RMSE為32.11 t·hm-2，估測精度明顯提升。根據云冷杉材積凈生長率，以及生物量擴展因子法，基于2016年香格里拉森林資源二類調查的小班數據計算了2019年云冷杉林地上生物量，其總生物量為39 961 196.39 t，138塊樣地的2019年生物量平均值為125.77 t·hm-2。相比于二者的結果，本研究的估測結果具有一定可靠性，但在精度方面依然有可提升的空間。非參數模型是一種黑箱操作原理，它是直接或間接地從實際系統的實驗分析中得到響應，機械式的系統化難以提高模型的預測精度，為提高模型的預測精度，獲取更加精準的生物量空間分布制圖，可對模型進行優化。宋涵玥等采用隨機搜索法（Random Searching）對模型參數進行優化，優化后的生物量估測精度提高2.01%。此外，還可采用遺傳優化算法、粒子群優化算法和貝葉斯優化方法對模型進行優化。

4 結論

本研究采用星載激光雷達GEDI數據協同Landsat8-OLI影像估測了滇西北生態系統較為完善的高寒山區香格里拉云冷杉林的地上生物量，評估多源遙感數據協同估測森林生物量的優勢。主要結論如下：

（1）利用地統計學GS+9軟件對GEDI觀測點做變異函數分析，經第一次篩選的12個GEDI變量變異函數的最優模型均為指數模型，模型決定系數均大于0.5。基于變異函數的分析結果，獲取12個變量的序貫高斯條件模擬空間分布，采用OEC檢驗其模擬精度，并進行第二次篩選，篩選后9個GEDI變量的OEC均高于0.9，較接近于1，且模擬估計的總體總值均在置信區間內，滿足精度要求，基于地統計學的序貫高斯條件模擬方法可實現GEDI光斑點由“點”到“面”的空間擴展。

（2）利用單一光學遙感Landsat8數據的1 6個植被指數和3個地形因子構建生物量估測模型，模型的R2=0.82，RMSE=35.51 t·hm-2，P-0.77；星載激光雷達GEDI數據協同光學遙感Landsat8數據構建的生物量估測模型R2=0.86，RMSE=32.11t·hm-2，P=0.80，精度明顯提升。多源遙感數據協同能夠互補不同的數據優勢，提高森林生物量的估測精度。

（3）利用2019年的兩種遙感數據，結合2016年野外調查數據，根據材積凈生長率和生物量擴展因子法獲取138塊云冷杉樣地2019年的生物量，利用隨機森林回歸方法估測的香格里拉2019年云冷杉林的地上生物量總量為37 042 605.68 t，平均生物量為123.28 t·hm-2，主要分布在38.54～282.19 t·hm-2之間。

（責任編輯：彭南軒）

基金項目：云南省農業聯合專項一重點項目（20230IBD070001-002）;云南省教育廳科學研究基金項目（2023Y0728），