張 晨,朱建軍,付海強(qiáng)

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

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

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

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

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

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

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

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

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

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

當(dāng)?shù)玫街脖桓叨群螅枰獙⒎囱莸闹脖桓弋a(chǎn)品及輔助數(shù)據(jù)進(jìn)行重采樣，使其與TanDEM-X DEM具備相同的分辨率(90 m)。TanDEM-X DEM通過X波段InSAR技術(shù)獲得，由于X波段SAR信號具有一定穿透能力，測得的高程并非位于冠層頂部，不能從TanDEM-X DEM中直接扣除植被高(Hveg)，需要考慮穿透深度的影響。進(jìn)而TanDEM-X DEM中扣除的植被高度Hrm可表示為：

Hrm=a×Hveg.

(1)

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

DEMbare-earth=DEM-Hrm.

(2)

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

(3)

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

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

2 實(shí)驗(yàn)結(jié)果與分析

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

圖2 實(shí)驗(yàn)區(qū)及ICESat-2數(shù)據(jù)分布與DEM

2.1 植被高反演結(jié)果

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

圖3 平均植被高圖

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

圖4 植被高反演結(jié)果驗(yàn)證

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

圖5 不同植被類型高度反演結(jié)果

2.2 林下地形計(jì)算結(jié)果與精度驗(yàn)證

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

圖6 扣除的植被高圖

圖7 DEM結(jié)果驗(yàn)證統(tǒng)計(jì)結(jié)果

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

3 結(jié)束語

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