基于ICESat-2數據及TanDEM-X DEM的林下地形反演

張 晨,朱建軍,付海強

(中南大學 地球科學與信息物理學院，湖南 長沙 410083)

數字高程模型(Digital Elevation Model, DEM)作為最重要的基礎地理信息產品之一[1]，已經在基礎工程建設、水資源管理、地質學、冰川研究、氣象以及軍事制導等研究領域發揮著重要的作用[2]。隨著對地觀測技術的發展，目前獲取的全球DEM向高精度、高分辨率、多樣性的方向發展。然而受現有對地觀測技術穿透性的限制，獲取的DEM數據中包含植被、建筑物等高度信息。尤其在植被覆蓋區域，傳統光學遙感及微波遙感手段獲取的是植被冠層的高程信息，并非準確的林下地表高程。為此，如何獲得準確的植被高，將其從上述DEM產品中移除，對于獲取大范圍、高精度地表高程信息至關重要。

現階段，主要的DEM產品有SRTM DEM[3]、ASTER GDEM[4]、AW3D30[5]及TanDEM-X DEM[6]。其中SRTM及TanDEM-X DEM分別通過C、X 波段InSAR (Synthetic Aperture Radar Interferometry)技術獲取，ASTER GDEM及AW3D30產品通過光學遙感手段獲取。在森林植被覆蓋區，這些DEM受全部或部分植被高度信號的干擾，無法反映林下地形真實情況。因此，要想得到林下地形，需從上述DEM產品中扣除植被高信號。目前，現有研究主要是對SRTM DEM中包含的植被高信號進行扣除。2007年，Wilson等通過實地調查植被高結合Landsat-7遙感影像估算C波段SAR信號的穿透深度為植被高的50%。而后，從SRTM剔除相應的植被高[7]。2012年，Gallant等利用Landsat TM數據進行地塊分析，通過最小二乘估算出地塊邊緣的植被高偏移量，并將其從SRTM DEM中扣除。但是，以上兩種方法僅限于小范圍，無法滿足大范圍林下地形制圖的需求[8]。2013年，Baugh等人利用ICESat-1GLAS數據測量結果發現：C波段SAR信號在亞馬遜熱帶雨林區的穿透深度為植被高的40%，隨后，結合全球植被高產品，從SRTM中扣除植被高度。盡管該方法適用于大范圍林下地形制圖，但對C波段SAR信號穿透深度的計算方法過于簡略，沒有充分考慮穿透深度在空間上的顯著差異[9]。針對該問題，2016年，O'Loughlin等利用ICESat-1GLAS、MODIS數據以及全球植被高產品，建立SRTM DEM中C波段SAR信號穿透深度與MODIS植被連續場(Vegetation Continuous Field, VCF)的關聯，之后，利用VCF預測穿透深度的空間變化，最終實現SRTM DEM中植被高信號扣除[2]。

上述研究工作主要利用Simard等[10]發布的全球植被高產品。但是，該全球植被高產品所用數據為2007年之前獲取，受植被高時空變化影響，該全球植被高產品無法用于TanDEM-X DEM(2018年發布)等新一代全球DEM中植被高信號的扣除。2019年5月NASA發布了ICESat-2/ATLAS(Advanced Topographic Laser Altimeter System)數據，其采用先進的光子點云數據采集技術，與ICESat-1/GLAS數據相比，ICESat-2/ATLAS具有更高的空間分辨率，對坡度具有更強的抵抗力，更易于測量山區地形以及植被高[12]。因此，如何利用ICESat-2/ATLAS數據獲取大范圍植被高并將其從TanDEM-X DEM中扣除，對提升DEM質量具有重要意義?；诖耍疚奶岢鲆环N基于ICESat-2/ATLAS數據和TanDEM-X DEM的林下地形生成方法，該方法主要包含三個步驟：①融合ICESat-2/ATLAS、 MODIS、溫度、降雨及高程輔助數據并結合隨機森林機器學習法獲取植被高；②基于得到的植被高數據、TanDEM-X DEM及VCF數據，利用隨機森林方法對X波段SAR數據的穿透深度進行建模及預測；③從TanDEM-X DEM扣除有效的植被高，得到大范圍、高精確的林下地形。最后，選用我國東北區域作為實驗區對所提算法的有效性進行驗證。

1 扣除TanDEM-X DEM中植被信號的林下地形反演

1.1 基于ICESat-2 ATLAS植被高度反演

NASA(National Aeronautics and Space Administration)于2018年9月15日發射ICESat-2衛星，用于測量冰蓋、冰川、冰架、海冰以及植被高度。ICESat-2搭載先進的地形激光高度計系統(ATLAS)傳感器，通過采集光子點云數據對地表進行測量，具有小光斑、多波束、高采樣率的特性[11]。采集的數據具有較高的空間分辨率，能記錄豐富的地形信息，對坡度影響具有較強的抵御能力，在林下地形測繪及植被冠層探測等方面具有巨大潛力。本文主要采用ICESat-2/ATLAS ATL08產品中的地形及植被高數據，該數據在沿軌方向采樣間隔為100 m[12]。由于地面和植被有較低的反射率，返回信號大部分來源于強光波束，因此需要根據升降軌來判別強弱光，從而選取強光信號；其次，數據采集期間，會有大量的噪聲光子存在，但是夜間數據存在的背景噪聲明顯小于白天，這里需要根據夜間標志(night_flag)對數據進行篩選、剔除；最后，需要根據云標志(cloud_flag_asr)以及植被標志(canopy_flag)對數據進行篩選，從而得到可靠數據。但是，基于ICESat-2/ATLAS數據提取的植被高度仍為空間離散點，需考慮如何借助其他輔助數據填補其空白區。

從空間大尺度角度來看，影響植被分布的因素主要有緯度變化、海陸分布以及海拔高度三種情況。緯度影響主要體現在不同緯度溫度不同，進而影響植被高度分布，即植被的緯度地帶性；海陸分布影響主要是從海洋向內陸水汽含量逐漸變少，從而影響降雨量的變化，繼而影響植被高度分布，即植被的經度地帶性；海拔高度影響主要表現在隨著海拔的變化溫度也隨之發生改變，使得植被分布隨著海拔的增加呈現垂直分布，即植被的垂直地帶性[13-14]?；谏鲜隹紤]，本文選取全球降水測量(Global Precipitation Measurement，GPM)數據、陸地地表溫度產品MOD11A2、植被連續場(Vegetation Continuous Fields，VCF)數據MOD44B以及高程數據TanDEM-X DEM幾個輔助數據進行植被高的反演。主要原理為通過尋求植被高度與上述輔助數據之間的關聯，利用得到的關聯及輔助數據推演空間連續植被高度圖。

所有輔助數據均被采樣至1 km的分辨率，如果一個分辨單元內包含多個ICESat-2點，則計算這些點的均值，以此表征該像元內的植被高。由于ICESat-2 ATLAS數據的精度受坡度變化的影響較小，因此，根據地面坡度對數據進行篩選。此外，盡管ICESat-2相較于ICESat-1具備較高的沿軌方向分辨率，但作為一種“點”數據，其難達到“面”測量的需求。為解決這一問題，本文利用ICESat-2/ATLAS(ATL08)植被高數據及上述輔助數據，利用隨機森林(Random Forest， RF)進行訓練，尋求二者之間的關聯，在此基礎之上，以輔助數據為輸入，利用RF對植被高度進行預測[11]。

1.2 扣除TanDEM-X DEM中植被高信號的林下地形反演

當得到植被高度后，需要將反演的植被高產品及輔助數據進行重采樣，使其與TanDEM-X DEM具備相同的分辨率(90 m)。TanDEM-X DEM通過X波段InSAR技術獲得，由于X波段SAR信號具有一定穿透能力，測得的高程并非位于冠層頂部，不能從TanDEM-X DEM中直接扣除植被高(Hveg)，需要考慮穿透深度的影響。進而TanDEM-X DEM中扣除的植被高度Hrm可表示為：

Hrm=a×Hveg.

(1)

式中：a為扣除后植被高與估算的植被高度的比例系數。隨后，通過移除DEM中的植被高度得到林下地形的高度信息，可以表述為：

DEMbare-earth=DEM-Hrm.

(2)

聯合式(1)和式(2)，有效植被占比可以表示為：

(3)

植被覆蓋區，InSAR穿透深度的主要與植被高、植被密度、地形、含水量等因素有關。為了準確估計每一個像素的扣除植被比，本文首先依據ICESat-2 ATL08產品提供的“裸”地表高程信息，計算其對應的扣除植被比；然后，通過RF方法構建扣除植被比a與植被高(Hveg)、VCF、高程(TanDEM-X DEM)之間的關系，依據所得到的模型，估算其余像素的扣除植被占比；最后，計算每一個像素對應的有效植被高度，并將其從TanDEM-X DEM數據中去除，即可得到林下地形信息。具體流程如圖1所示。

圖1 扣除TanDEM-X DEM中植被高信號的林下地形反演流程圖

2 實驗結果與分析

為了驗證算法有效性，本文選取我國東北區域作為實驗區，主要包括黑龍江省，吉林省，遼寧省，內蒙古東部(呼倫貝爾市，興安盟，通遼市和赤峰市)以及河北北部(承德市，秦皇島市，唐山市以及張家口市)。該區域屬于溫帶季風氣候，夏季高溫多雨、冬季寒冷干燥。植被覆蓋面積廣闊，以針葉林和針闊混交林為主。ICESat-2數據分布及DEM圖2所示。

圖2 實驗區及ICESat-2數據分布與DEM

2.1 植被高反演結果

基于圖1算法流程，本文通過對ICESat-2 ATLAS數據以及輔助數據的處理提取計算、激光光斑植被高建模，使用隨機森林算法外推出沒有ICESat-2 ATLAS數據區域的植被高，得到研究區域植被高分布圖，如圖3所示。其中，大部分植被高在13 m左右，高植被區域主要在吉林省的東部、遼寧省的南部以及黑龍江省伊春市。占據植被分布面積較多的大興安嶺區域植被高主要集中在13～20 m。

圖3 平均植被高圖

對上述反演結果，用TanDEM-X森林/非森林(TanDEM-X Forest/Non-Forest，TanDEM-X FNF)數據進行掩膜，最后采用ICESat-2 數據中預留的少部分數據進行最終的精度評定。植被高制圖結果一個分辨單元內若存在多個驗證數據點，則計算其平均值，用該平均值代表該像元的平均植被高。如圖4所示，本文反演的植被高結果與驗證數據具有較好的一致性。植被高結果均方根誤差RMSE為2.96 m，相關系數R2為0.87。誤差主要分布在±5 m之間。反演結果存在低矮植被區域高估，高植被區域略微低估的現象。造成上述誤差的原因，主要包含采樣尺度誤差，模型誤差以及地理定位誤差。采樣尺度誤差主要體現在：為了使得數據保持相同分辨率，從而進行采樣處理；模型誤差主要是隨機森林建模過程中的誤差；地理定位誤差主要體現在多幅遙感影像進行糾正時精度不一致導致的誤差。

圖4 植被高反演結果驗證

本文還根據全球土地覆蓋圖GlobCover 2009對植被類型分類進行精度評定。研究區域主要有闊葉林，針葉落葉常綠林以及闊葉針葉混交林幾種植被類型。由圖5可知，幾種林分反演結果與總體結果趨勢相近。該區域，優勢樹種是針葉林，主要集中在大興安嶺地區，其反演精度也相對較高，均值14.21 m，RMSE為2.55 m，R2為0.86，相對精度為17.96%。其次是闊葉林(>5m)占據較大面積分布，主要集中在黑龍江省伊春市以及遼寧省的東南部，均值為14.12 m，RMSE為3.0 m，R2為0.85，相對精度為21.21%。最后是混交林，其反演精度與針葉林相當，植被高均值為12.51m，RMSE為2.53 m，R2為0.87，相對精度為22.19%。

圖5 不同植被類型高度反演結果

2.2 林下地形計算結果與精度驗證

利用隨機森林建模計算出扣除植被的百分比，在植被高制圖的基礎上得到TanDEM-X DEM中包含的植被高信號，即需要扣除的植被高，如圖6所示。將得到的結果從TanDEM-X DEM中移除，得到林下地形。為了驗證算法的正確性，通過ICESat-2地面高程點(815 710個)，分別對改正前、后的DEM驗證進行精度驗證。其統計結果如圖7所示。

圖6 扣除的植被高圖

圖7 DEM結果驗證統計結果

從圖7(b)可以看出校正前DEM誤差整體表現出負的偏態，說明TanDEM-X DEM中包含了明顯的植被高信號。而經過改正后的DEM，由于移除了植被高的影響，其誤差分布在0附近。為了進一步驗證精度，分別計算了校正前、后的均方根誤差(RMSE)和誤差均值。校正后DEM的RMSE為9.14 m，相較于原始的DEM(RMSE=11.65 m)，其精度提高了約21.6%；且其誤差大多集中在±20 m之間。誤差均值也由-6.5 m改善至0.03 m。綜上結果表明：本文所提的方法可以較好的校正現有DEM中包含的植被高度信號。

3 結束語

本文以東北區域為例，基于輔助數據和ICESat-2/ATLAS數據，采用隨機森林建模方法，反演植被高。植被高結果RMSE為2.96 m。根據研究區域主要植被類型的不同，對不同類型植被進行精度評定，結果表明該區域的優勢樹種反演結果精度相對較好。其后，利用所得的植被高產品，采用隨機森林建模，推算出研究區域的林下地形，RMSE為9.14 m，相比原DEM精度提高了21.2%。結果表明本文方法可較好的移除現有DEM中植被高信號，為快速生成大范圍、高精度的DEM提供了參考。盡管星載激光雷達具備反演大范圍植被高的能力，但分辨率過低，適用于大尺度分析，難以滿足高分辨率林下地形制圖的需求。如何進一步提高植被高反演結果的分辨率將是后期研究工作的重點。