ZY3立體像對和機載LiDAR抽樣數據協同估測森林平均高*

趙俊鵬 趙 磊 陳爾學 萬祥星 徐昆鵬

(中國林業科學研究院資源信息研究所 國家林業和草原局林業遙感與信息技術重點實驗室 北京 100091)

森林是陸地生態系統的主體和重要資源，是人類社會生存發展的重要生態屏障，在保持水土、涵養水源、固碳釋氧、凈化空氣、應對氣候變化、美化環境等方面發揮著不可替代的作用(Heseetal.， 2005； Le Toanetal.， 2011； 王效科等， 2001)。作為森林生態系統最基本的數量表征，森林參數表達了森林與環境間物質循環與能量流動的關系，是評估森林碳源和碳匯的重要標志。森林高度是重要的森林參數，能夠反映森林垂直結構信息，對森林蓄積量、地上生物量等參數估測頗為重要。

激光雷達(light direction and ranging，LiDAR)是一種可主動獲取觀測對象信息的現代光學遙感技術，其以高精度、高分辨率等特點，在獲取森林垂直結構參數方面具有明顯優勢(龐勇等， 2005； Choppingetal.， 2009)，通過機載LiDAR系統獲取高密度激光點云數據能夠提取出高精度的森林高度信息(Popescuetal.， 2002； 龐勇等， 2008)； 但由于大范圍獲取同一區域的多期機載LiDAR數據實現森林資源動態監測成本很高，目前在我國大部分地區實施起來還十分困難。將機載LiDAR與星載遙感數據相結合實現大面積森林參數估測，可能是最有效且成本最低的方法(Tsuietal.， 2013)； 但該方法需要采用抽樣式LiDAR數據獲取模式以降低LiDAR數據獲取成本，同時要求衛星遙感數據完全覆蓋整個目標制圖區域，且衛星遙感特征對LiDAR提取的森林參數敏感。近年來，國外學者利用從Landsat、MODIS等多光譜影像中提取的波段反射率、植被指數、圖像紋理等特征結合機載LiDAR數據，在森林高度估測方面開展了大量探索(Hudaketal.， 2002； Chenetal.， 2011； Yuetal.， 2015； Wangetal.， 2018； Zhuetal.， 2020)，但采用的中低空間分辨率遙感光譜類紋理特征用于森林高度探測的機理并不明確，在森林資源規劃設計調查業務中應用受到限制。

高空間分辨率光學立體攝影測量遙感技術具有獲取高精度森林冠層頂部高度信息[相當于數字表面模型(digital surface model,DSM)]的能力，若再已知林下地形信息[相當于數字高程模型(digital elevation model,DEM)]，就可以通過二者差值得到森林冠層高度信息(canopy height model,CHM)。倪文儉等(2018)由GF-2異軌立體數據結合機載LiDAR獲取的DEM數據提取得到CHM，其與LiDAR提取的CHM存在明顯相關性，可用于森林高度估測。邢艷秋等(2018)以吉林省長春市凈月潭國家森林公園為研究區，通過ZY3不同立體像對組合提取的DSM與機載LiDAR提取的DEM作差得到CHM，進而以CHM、ZY3多光譜特征為遙感特征(自變量)，與實測樣地林分高(因變量)數據建立線性回歸模型對森林高度進行估測，既保證了較高的樹高模型估測精度，又實現了較大連續區域快速的林分平均高估測，亦成功試驗了將測繪數據應用于林分測高。以上研究初步展示了高空間分辨率光學立體攝影測量結合機載LiDAR數據提取的CHM具有估測森林高度的潛力，但并沒有針對森林資源調查業務情景開展應用方法探討。

自遙感數據提取出對森林參數敏感的遙感特征(自變量)且獲得若干樣地的待估測森林參數(因變量)實測值后，可以通過建立因變量與自變量的統計模型實現每個遙感分辨單元的森林參數估測。國內外最常用的統計模型包括多元線性回歸、逐步回歸、k-鄰近(KNN)、支持向量機回歸、隨機森林回歸等，但這些模型通常沒有考慮變量的空間統計特性，具有空間屬性的森林，在生長過程中會受各種因素影響而存在空間效應(馮益明等， 2004)。空間統計學則提供了有效分析和解釋空間數據關系的方法(Rossietal.， 1992)，在應用空間統計學結合遙感數據進行森林參數反演的相關工作中，國外已有學者進行了大量探索，結果發現利用克里格(Kriging)、外部漂移克里格(Kriging with external drift,KED)、協同克里格(co-Kriging，CK)、回歸克里格(regression Kriging,RK)等方法(Salesetal.， 2007； Hernndezetal.， 2011； Fayadetal.， 2016； Pouladietal.， 2019)進行森林參數估測可以得到較高的估測精度。國內主要相關研究，如馮益明等(2006)應用空間統計學半方差理論對QuickBird影像的全色波段進行分析，實現了人工白皮松(Pinusbungeana)林的冠幅估測； 郭含茹等(2015)以臺州市仙居縣為研究區，基于Landsat TM影像和森林資源二類調查數據，應用地理加權回歸方法對森林碳儲量進行了估測； 王紫君等(2016)利用Landsat ETM+數據，在逐步回歸分析基礎上，引入模型殘差的反距離插值結果，對上海城市森林生物量分布特征進行了估測； 閭妍宇等(2017)也應用地理加權回歸方法，結合森林資源二類調查數據和Landsat TM影像，對景谷縣的思茅松(Pinuskesiyavar.langbianensis)生物量進行了估測。但國內這些研究并沒有或很少采用回歸克里格方法，所用遙感特征也只是常規的光譜類特征，更不是針對森林平均高估測的應用目標。

廣西壯族自治區2018年實現了南寧市高峰林場機載LiDAR數據和地面樣地數據的獲取，完成了森林平均高、蓄積量等森林資源專題信息的提取，同時也得到了相應區域的高精度DEM信息。在當前數據條件下(已具備高精度的DEM)，如何更經濟(不依賴全覆蓋的機載LiDAR數據)、更有效地更新森林平均高專題圖？本研究針對該問題，探索一種適用于已具備林下地形，可協同利用少量實測樣地數據、抽樣式采集的機載激光雷達(LiDAR)條帶數據和區域全覆蓋的資源三(ZY3)立體像對數據有效估測區域森林平均高的方法，以期為提高森林資源調查效率和精度提供技術支撐。

1 研究區概況

以廣西高峰林場界牌和東升分場為研究區。圖1b中黃色區域表示高峰林場轄區，圖1c和d中2個紅色多邊形區域表示2個分場邊界。研究區總面積約52 km2，地理坐標為108°19′52″—108°26′10″E，22°56′50″—23°0′35″N。高程77.88～463.65 m，坡度0°～74.9°，地形十分陡峭。屬亞熱帶季風氣候，年均氣溫21.6 ℃，年均降水量1 200～1 500 mm。主要樹種包括桉樹(Eucalyptus)、杉木(Cunninghamialanceolata)、馬尾松(Pinusmassoniana)以及其他闊葉樹種。

圖1 研究區位置和覆蓋范圍示意Fig.1 Location and coverage map of study area

2 數據與方法

2.1 數據與處理

2.1.1 實測數據 實測數據來源于2018年1月13—25日開展的外業樣地調查，共調查54塊20 m×20 m方形樣地，其中37塊樣地用于模型構建、17塊樣地用于精度評價。對樣地內所有胸徑大于5 cm的樹木每木檢尺，測量因子包括胸徑、樹高和樹種等。采用斷面積加權法計算森林平均高，公式如下：

(1)

式中：H為樣地內森林平均高；hi為第i株林木高度；gi為第i株林木胸高斷面積；n為樣地內林木株數。

2.1.2 機載LiDAR數據 機載LiDAR數據采集于2018年1月，使用CAF-LiCHy數據采集系統(Pangetal.， 2016)，其中LiDAR傳感器型號為Riegl LMS-Q680i，數據采集時傳感器參數配置見表1。將覆蓋整個研究區的LiDAR點云數據采用Kraus濾波算法得到地面激光點，通過計算每個像元內激光點高度平均值生成空間分辨率為1 m的DEM，為與其他數據空間分辨率保持一致，本研究將DEM的空間分辨率重采樣至20 m。基于預處理的點云數據，對其進行高度歸一化，提取分位數高度變量，結合37塊實測樣地的森林平均高與點云分位數高度變量構建反演模型，得到整個研究區LiDAR估測的森林平均高(HLi)，具體數據處理和信息提取方法參考曹林等(2015)。

表1 LiDAR傳感器參數配置Tab.1 Parameter configuration of the LiDAR sensor

2.1.3 ZY3數據 獲取覆蓋整個研究區的一景ZY3數據，成像時間為2018年3月，1臺2.1 m 分辨率和2臺3.5 m分辨率全色相機按照下視、前視、后視方式排列，各視影像組合可組成立體像對。ZY3立體像對數據通過密集點匹配技術產生點云，并由此生成DSM，將其與LiDAR數據提取的DEM相減，可得到覆蓋整個研究區的CHM(Lietal.， 2019a)，記為CHMZY3，作為反演森林平均高的特征用于模型構建。

2.2 研究方法

2.2.1 技術路線 高空間分辨率光學立體攝影測量遙感技術具有獲取高精度DSM的能力，在森林覆蓋區獲取的DSM代表森林冠層頂部高度信息，但要測量森林高度，還需要已知林下地形信息，即需要高精度的DEM。由于光學遙感波段電磁波對樹冠沒有穿透能力，只能靠樹冠間的空隙到達林下地表，因此光學立體攝影測量在高郁閉度森林很難獲取高精度的DEM。相比之下，機載LiDAR通過獲取冠層和地表的高密度激光點云數據，可以同時提取出高精度的DSM和DEM，進而得到CHM。ZY3是我國于2012年發射的第一顆民用三線陣立體測繪衛星，其所產生的DSM具有覆蓋范圍大、成本低的特點，只要目標區域已具備高精度的DEM，就很容易得到CHMZY3；但得到的CHMZY3不是森林平均高，只是一個與森林平均高有很好關系的遙感特征，通常需要采用若干地面實測樣地數據，建立一個由該遙感特征(自變量X)估測森林平均高(因變量Y，由樣地數據計算得到)的模型，基于該模型可將CHMZY3轉化為森林平均高。實測樣地數據的獲取是耗費人力、物力的，本研究希望探索通過抽樣式采集機載LiDAR數據，配合少量地面實測樣地數據，分2個階段建立估測模型的方法。第一階段先建立由LiDAR數據估測森林平均高的模型(Y為實測樣地森林平均高，X為LiDAR提取的點云特征變量)，得到LiDAR抽樣數據對應的森林平均高(HLi)； 第二階段以HLi為實測森林平均高(Y)，建立由CHMZY3特征(X)估測森林平均高的模型。上述針對森林資源調查業務應用情景的技術路線如圖2所示。在本研究中，已有林下地形為LiDAR數據采集得到的DEM。

圖2 針對業務應用情景的技術路線Fig.2 Technical route for practical application scenario

圖3 試驗數據分布Fig.3 Distribution of the experimental data

2.2.2 機載LiDAR抽樣數據模擬 從全覆蓋的LiDAR數據中抽取12條飛行條帶，如圖3所示，水平和垂直方向上分別有5條和7條條帶，用于模擬抽樣式采集的LiDAR數據。條帶寬度600 m、條帶中心線之間間隔1.5 km，數據面積占整個研究區面積的65%。在抽取的條帶內隨機選擇3 000個像素點用于森林平均高反演模型的構建。

2.2.3 估測模型 1) OLS模型 普通最小二乘(ordinary least squares,OLS)模型是一種最小化回歸模型殘差平方和的估計方法，其可以根據一組數據求出自變量與因變量之間的函數關系，并使估測值與實測值之間的殘差平方和最小。模型如式(2)所示：

(2)

式中：YOLS為OLS模型估測值；β0為模型常數項；βi為變量的權重；p為自變量個數；Xi為自變量；ε為隨機誤差。

2) KNN模型k-近鄰(k-nearest neighbor, KNN)模型是一種非參數估測方法，其基于預測點與觀測數據特征之間的距離關系，從而對待預測點的值進行估測(McRobertsetal.， 2002)。模型如式(3)所示：

(3)

式中：YKNN為KNN模型估測值；Wi為不同觀測點的距離權重；Yi為第i個觀測數據的值；k為距離預測點最近觀測數據的個數。

3) RK模型 回歸克里格(regression Kriging, RK)模型是一種將線性模型或其他模型與模型殘差空間插值相結合的估計方法(Odehetal.， 1995)。本研究中，首先采用OLS模型和KNN模型估測整個研究區的森林高度； 然后采用簡單克里格方法對2種模型的殘差(建模數據實測值與模型估測值的差)分別進行空間插值； 最后將模型估測結果與殘差插值結果相加，得到最終估測結果，分別記為RKOLS和RKKNN。模型如式(4)所示：

RKOLS/KNN=YOLS/KNN+Krad。

(4)

式中： RKOLS/KNN為OLS或KNN模型對應RK模型的估測結果；YOLS/KNN為OLS或KNN模型估測值；Krad為OLS或KNN模型對應殘差的空間插值。

2.2.4 模型估測精度評價 為比較幾種模型的估測精度，以17塊未參與建模的樣地數據為驗證樣本對模型估測精度進行檢驗。選用相關系數(correlation coefficient,r)、平均誤差(mean orror, ME)、均方根誤差(root mean square error, RMSE)、估測精度(estimate accuracy, EA)作為評價指標：

(5)

(6)

(7)

(8)

3 結果與分析

3.1 模型擬合

在OLS模型中，自變量X僅有1個(p=1)，即CHMZY3，因變量Y為HLi，經模型擬合，得到如式(9)所示回歸方程：

HLi=6.008 9+0.568 8×CHMZY3。

(9)

模型F統計量為178(P<0.01)，表明回歸分析中的因變量與自變量極相關，模型擬合效果較好。

k的選擇對KNN模型最終估測結果影響很大，已有研究(鄭剛等， 2010； 向安民等， 2017)顯示，k在1～10范圍內，對不同參數的估測效果較好。為得到最優模型參數，本研究對建模數據采用交叉驗證方法，基于10折交叉驗證結果，選擇模型最優時的k值。由圖4可知，當k=6時，建模數據交叉驗證得到的RMSE最小，為2.21 m。因此，本研究選擇6作為森林平均高估測的最佳k值。

圖4 不同k值下模型的均方根誤差Fig.4 The RMSE of models with different k values

在RK模型中，分別采用指數、Stable模型擬合OLS、KNN模型殘差的變異函數，并對變異函數結構進行分析，從中選擇最優模型。塊金值反映區域化變量的隨機性大小， 基臺值反映系統內總的變異， 變程反映變量在空間上的影響距離，從表2可以看出，2種擬合模型的基臺值相近，但變程不同，OLS和KNN殘差的2種擬合模型中，Stable和指數模型分別能夠獲得更大的變程。從圖5可以看出，2種模型殘差的變異函數擬合結果均較好。

表2 模型殘差的不同變異函數系數Tab.2 Parameters of different variogram models for residuals

3.2 估測效果

使用未參與建模樣地數據對2類4種模型的估測精度進行檢驗，結果如圖6和表3所示。從圖6可以看出，RKOLS和RKKNN模型相比OLS和KNN模型得到的估測值與實測值之間的相關性更好； 從表3可以看出，考慮模型殘差空間關系的RKOLS和RKKNN模型能夠得到比OLS和KNN模型更高的估測精度。RKOLS模型的估測效果最好，得到了最高的估測精度(EA=87.42%)以及最低的平均誤差和均方根誤差(ME=-0.16 m，RMSE=1.84 m)； 其次是RKKNN模型(EA=87.32%)。在對比分析各種指標的基礎上，將估測效果最好的RKOLS模型應用于整個研究區，得到研究區內森林平均高估測結果如圖7所示。

圖5 OLS和KNN模型殘差的變異函數Fig.5 Variograms and fitted models of residuals from OLS and KNN

圖6 實測數據與估測結果比較Fig.6 Scatter plots of predicted versus validation data

表3 估測精度比較Tab.3 The accuracy comparison of four kinds of models

圖7 研究區森林平均高估測結果Fig.7 Spatial distribution of the estimated MHF

4 討論

本研究基于研究區2018年獲取的機載LiDAR數據、ZY3立體像對數據和少量實測樣地數據，通過合理假設和數據“模擬”，對森林資源調查業務應用情景的森林平均高估測技術路線進行驗證。該技術路線實際應用時需要目標區域具備高精度的DEM，目前一般通過采集覆蓋目標區域的高密度機載LiDAR點云數據實現。由于DEM具有長期穩定不變的特點，因此只要目標區域獲取過一次全覆蓋的機載LiDAR數據，后續對目標區域森林平均高的估測或更新就可以通過本研究建議的技術路線(圖2)實現。除了機載LiDAR外，長波長極化干涉SAR也具有測量林下地形的能力，而無論是機載還是星載極化干涉SAR系統，最近5年均會有快速發展，如我國機載極化干涉SAR系統目前已在開展應用校飛試驗，歐洲空間局P波段雷達生物量探測衛星(BIOMASS)也計劃于2022年前后發射，我國也將P波段極化SAR衛星納入了未來研發規劃(Le Toanetal.， 2011； Lietal.， 2019b)。

圖2技術路線的重要輸入之一是ZY3立體像對數據，這對立體測圖衛星來說比較容易實現目標區域1年1次全覆蓋觀測，從而可大大提高森林平均高的監測效率。目前已在軌運行的GF-7是一顆雙線陣立體測圖衛星，DSM的精測精度相比ZY3有大幅度提升； 另外，我國正在研制的陸地生態系統碳監測衛星，將搭載多波束激光雷達、多角度多光譜等有效載荷，也具有高精度冠層頂部高程信息精準提取能力。

按照圖2技術路線，還需獲取機載LiDAR抽樣數據和少量實測樣地數據，主要是為森林平均高估測模型建立提供具有代表性的地面實況。以往該過程通常只依賴樣地調查數據，需要大量樣地調查才能保證模型具有足夠的代表性，費時費力； 但采用抽樣式獲取LiDAR數據、分兩階段建立估測模型的方法能夠大大降低人力、物力投入，同時也可增加地面實況數據的樣本量，有利于建立更具有代表性、更穩健的估測模型。而且，當前無人機LiDAR測樹技術已比較成熟，對林場尺度的業務應用來說，利用無人機實現若干LiDAR抽樣數據的獲取也是一種更靈活、高效和低成本的地面實況數據獲取方法，這無疑使本研究建議技術路線的具體實施具有更高的可行性。

5 結論

1) 本研究中2類4種模型均可有效估測森林平均高，回歸克里格模型(RKOLS、RKKNN)優于非空間模型(OLS、KNN)，RKOLS模型估測效果最好。OLS和KNN模型的RMSE分別為1.88和1.96 m，EA分別為87.18%和86.64%； RKKNN和RKOLS模型的RMSE分別為1.86和1.84 m，EA分別為87.32%和87.42%。

2) 在林下地形已知時，協同利用少量實測樣地數據、抽樣式采集的機載LiDAR條帶數據和區域全覆蓋的ZY3立體像對數據能夠實現區域森林平均高的高效、高精度估測。

3) 已在軌運行的GF-7衛星和正在研制的陸地生態系統碳監測衛星、P波段極化SAR衛星以及無人機LiDAR等遙感手段，可靈活應用于本研究建議的技術路線，對提高森林資源調查效率和精度具有重要支撐作用。