基于激光雷達(dá)的大別山毛竹冠層高度反演研究

李玲娜，孫 悅，萬(wàn)芳芳，徐小牛，黃興召

(安徽農(nóng)業(yè)大學(xué)林學(xué)與園林學(xué)院，安徽 合肥 230031)

毛竹(Phyllostachysedulis)適應(yīng)性強(qiáng)、生長(zhǎng)快、經(jīng)濟(jì)價(jià)值高，也是我國(guó)竹類資源中分布最廣的竹種?，F(xiàn)有的毛竹林主要屬于筍材兼用林，是山村林農(nóng)重要的經(jīng)濟(jì)支柱[1]。根據(jù)第九次森林資源清查數(shù)據(jù)，我國(guó)竹林面積為6.41×106hm2，其中毛竹林面積4.67×106hm2，占竹林面積的72.9%；毛竹株數(shù)1.41×1011株，胸徑主要分布在7～11 cm，占61.9%[2-4]。如何進(jìn)一步精準(zhǔn)獲取毛竹林的冠層高度，是掌握毛竹資源結(jié)構(gòu)特征及可持續(xù)管理利用的關(guān)鍵。

激光雷達(dá)是近年來(lái)發(fā)展較快的一種主動(dòng)遙感技術(shù)，為準(zhǔn)確提取森林結(jié)構(gòu)參數(shù)提供了有力的技術(shù)支撐[5-6]。搭載在有人飛機(jī)或無(wú)人機(jī)上的小光斑激光雷達(dá)統(tǒng)稱為機(jī)載激光雷達(dá)(Airborne LiDAR)[7]。利用機(jī)載激光雷達(dá)數(shù)據(jù)生成數(shù)字高程模型(Digital Elevation Model，DEM)和數(shù)字表面模型(Digital Surface Model，DSM)，計(jì)算出冠層高度模型(Canopy Height Model，CHM)[8]。機(jī)載激光雷達(dá)獲取高精度的點(diǎn)云數(shù)據(jù)，能夠更加精確的估算毛竹冠層高度[9]，但是機(jī)載激光雷達(dá)做不到大尺度連續(xù)的觀測(cè)，需要結(jié)合其他大尺度連續(xù)的數(shù)據(jù)共同反演毛竹冠層高度，且機(jī)載激光雷達(dá)數(shù)據(jù)獲取成本大，很難進(jìn)行大規(guī)模的推廣[10]。星載激光雷達(dá)以衛(wèi)星為搭載平臺(tái)，運(yùn)行軌道高，觀測(cè)范圍廣，全天時(shí)全天候?qū)Φ剡M(jìn)行觀測(cè)[11]。高精度的機(jī)載激光雷達(dá)數(shù)據(jù)和大尺度的星載激光雷達(dá)數(shù)據(jù)相結(jié)合，進(jìn)一步提升了該方法的應(yīng)用范圍[12]。

全球生態(tài)系統(tǒng)動(dòng)態(tài)監(jiān)測(cè)(Global Ecosystem Dynamics Investigation，GEDI)是目前唯一專門為測(cè)量植被的垂直結(jié)構(gòu)而設(shè)計(jì)的全波形激光雷達(dá)，其核心目標(biāo)是估算熱帶和溫帶森林的地上生物量，模擬地表碳收支，了解植被結(jié)構(gòu)對(duì)生物多樣性的影響[13-15]。近年來(lái)，有關(guān)GEDI反演森林冠層高度已有較多研究，但冠層高度受樹種、氣候環(huán)境、地形等因子的影響，導(dǎo)致模型預(yù)估效果存在差異[16-18]。因此，利用GEDI數(shù)據(jù)，開展不同樹種的冠層高度反演研究，是提升模型預(yù)估精度的基礎(chǔ)[19-22]。為此，以毛竹為研究對(duì)象，建立GEDI反演模型，探討不同因子對(duì)模型的影響，并找出影響模型預(yù)估的關(guān)鍵因子，以期為毛竹林的結(jié)構(gòu)監(jiān)測(cè)和可持續(xù)利用提供技術(shù)支撐。

1 試驗(yàn)地概況與研究方法

1.1 試驗(yàn)地概括

實(shí)驗(yàn)選擇安徽省金寨縣關(guān)廟鄉(xiāng)竹園村(115°25′E，31°30′N)為研究地，該地平均海拔450 m，屬于北亞熱帶季風(fēng)氣候，年平均溫度14.6 ℃，極端高溫和低溫分別為35 ℃和-14 ℃。地帶性植被為常綠闊葉林，毛竹林為該地主要的森林植被類型之一。土壤主要是黃棕壤，土層厚度在0.5～1.2 m，土壤pH值在5.5～6.5之間。2021年7月，選擇代表性立地條件及竹林結(jié)構(gòu)，建立25 m×25 m的標(biāo)準(zhǔn)樣地28塊(見圖1)。在樣地內(nèi)進(jìn)行每木檢尺，利用超聲波樹高測(cè)定儀(LD6172)測(cè)定樹高，作為模型預(yù)估效果驗(yàn)證數(shù)據(jù)。

圖1 毛竹林臨時(shí)樣地分布圖Fig.1 The distribution of Moso bamboo sampling plots

1.2 研究數(shù)據(jù)獲取

1.2.1星載激光雷達(dá)數(shù)據(jù) GEDI星載全波形激光雷達(dá)數(shù)據(jù)于2021年7月在NASA/USGS (National Aeronautics and Space Administration/United States Geological Survey)土地處理分布式活動(dòng)檔案中心發(fā)布，GEDI二級(jí)產(chǎn)品包括地面高程和垂直結(jié)構(gòu)參數(shù)[23-25]。提取所需地區(qū)GEDI光斑內(nèi)的經(jīng)緯度、不同百分位樹高值(RH)、地面高程、覆蓋度、波形參數(shù)、信號(hào)強(qiáng)度、平均海平面等參數(shù)。

1.2.2機(jī)載激光雷達(dá)數(shù)據(jù) 機(jī)載激光雷達(dá)數(shù)據(jù)(2021年7月)通過(guò)規(guī)劃的航線，使用Riegl·VUX-1無(wú)人機(jī)激光雷達(dá)系統(tǒng)獲取平均點(diǎn)密度42 pts·m-2，飛行高度 90 m的數(shù)據(jù)(飛行區(qū)包含設(shè)置的28塊樣地)。獲取到機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)后，利用點(diǎn)云魔方軟件自帶的功能處理，獲得數(shù)字地形模型(Digital Terrain Model，DTM)和數(shù)字地表模型(Digital Surface Model，DSM)，通過(guò)計(jì)算可獲得冠層高度模型(Canopy Height Model，CHM)[26-28]。

1.3 研究方法

1.3.1星載激光雷達(dá)數(shù)據(jù)的提取與篩選 根據(jù)獲取的機(jī)載激光雷達(dá)數(shù)據(jù)的經(jīng)緯度參數(shù)，我們?cè)诰W(wǎng)站進(jìn)行篩選，并下載所選28塊樣地星載激光雷達(dá)數(shù)據(jù)[28-29]。提取的參數(shù)如表1所示。

表1 提取參數(shù)介紹

為了保證獲取的星載激光雷達(dá)數(shù)據(jù)更具有，我們除了篩選經(jīng)緯度，還對(duì)其他參數(shù)進(jìn)行了處理。

(1)degrade_flag=0，若該值為“1”，表示指示指向和/定位信息的狀態(tài)下降，因此刪除degrade_flag=1的光斑，保留degrade_flag=0的光斑。

(2)quality_flag=1，質(zhì)量標(biāo)志為“1”表示波形滿足基于能量、靈敏度、振幅和實(shí)時(shí)表面跟蹤質(zhì)量的特定標(biāo)準(zhǔn)，可以表示為有效波形。

(3)surface_flag=0，表面標(biāo)志為“0”表示elev_lowestmode在 DEM 或 MSS 的 300 m范圍內(nèi)，可以表示該參數(shù)有效。

(4)sensitivity≥0.9，大于等于0.9的靈敏度閾值表示光斑質(zhì)量較好，光斑能夠穿透毛竹冠層打到地面，因此刪除靈敏度小于0.9的光斑。

(5)|elev_lowestmode-digital_elevation_model-mean_sea_surface|<50，如果星載激光雷達(dá)數(shù)據(jù)獲取的(地面高程-平均海平面高度)與TanDEM-X的絕對(duì)值差值大于50 m，則被認(rèn)為是受云層影響的錯(cuò)誤數(shù)據(jù)，即刪除該光斑。

(6)RH100≤30，通常來(lái)說(shuō)，毛竹不會(huì)超過(guò)30 m，如果出現(xiàn)了大于30 m的數(shù)據(jù)則被認(rèn)為是受云霧影響或非毛竹冠層高度，數(shù)據(jù)不正確。

(7)landsat_treecover>0，modis_treecover>0，由于毛竹與闊葉樹相比冠幅小，激光更易打到地面以造成誤差，故刪除覆蓋度不是大于零的錯(cuò)誤值。

1.3.2機(jī)載激光雷達(dá)數(shù)據(jù)重采樣 機(jī)載激光雷達(dá)數(shù)據(jù)較為精準(zhǔn)，具有1 m的高分辨率，GEDI光斑直徑約為25 m，因此需要對(duì)機(jī)載激光雷達(dá)進(jìn)行重采樣，降低其分辨率至25 m[30-31]。首先，研究使用了ArcMap最鄰近法對(duì)DTM進(jìn)行重采樣；其次，我們用Matlab軟件編寫重采樣代碼，提取CHM不同的百分位樹高值(CHM 90～100)；最后，使用ArcMap表面分析工具提取DTM中的坡度、坡向等地形指數(shù)數(shù)據(jù)。

1.3.3星載激光雷達(dá)與機(jī)載激光雷達(dá)數(shù)據(jù)融合 星載激光雷達(dá)與機(jī)載激光雷達(dá)數(shù)據(jù)的融合，既能解決機(jī)載激光雷達(dá)數(shù)據(jù)獲取成本高的問(wèn)題，又能解決星載激光雷達(dá)數(shù)據(jù)分布不連續(xù)的問(wèn)題。首先運(yùn)用ArcMap的投影與變換將重采樣過(guò)的星載激光雷達(dá)數(shù)據(jù)與機(jī)載激光雷達(dá)數(shù)據(jù)投影到同一坐標(biāo)系；其次根據(jù)經(jīng)緯度坐標(biāo)，將CHM 90～100、坡度和坡向與星載激光雷達(dá)數(shù)據(jù)融合；最后剔除融合數(shù)據(jù)的重復(fù)值與異常值。研究使用70%的大別山毛竹林機(jī)載激光雷達(dá)數(shù)據(jù)進(jìn)行建模，30%的數(shù)據(jù)進(jìn)行模型驗(yàn)證。

1.4 模型構(gòu)建與驗(yàn)證

首先對(duì)星載數(shù)據(jù)RH 90～100與CHM 90～100進(jìn)行擬合分析，找出擬合程度最優(yōu)的一種組合，即為機(jī)載數(shù)據(jù)最匹配的百分位樹高值，以避免機(jī)載激光雷達(dá)數(shù)據(jù)對(duì)模型估算產(chǎn)生誤差。然后進(jìn)一步分析星載激光雷達(dá)的不同特征參數(shù)對(duì)毛竹冠層高度精確度的影響，我們共提取5個(gè)變量，分別為植被覆蓋度、太陽(yáng)高度角、太陽(yáng)方位角、坡度與坡向，利用隨機(jī)森林對(duì)5個(gè)特征參數(shù)進(jìn)行判別，找出關(guān)鍵影響因子。再利用逐步回歸構(gòu)建的毛竹冠層高度估算模型，根據(jù)均方根誤差(RMSE)和決定系數(shù)(R2)評(píng)價(jià)建模效果，選出最優(yōu)模型。最后，利用樣地實(shí)測(cè)的毛竹數(shù)據(jù)驗(yàn)證模型預(yù)估效果。

2 結(jié)果與分析

2.1 激光雷達(dá)數(shù)據(jù)融合

我們將網(wǎng)站下載的星載激光雷達(dá)數(shù)據(jù)與無(wú)人機(jī)獲取的機(jī)載激光雷達(dá)數(shù)據(jù)在ArcMap中進(jìn)行融合，計(jì)算每個(gè)區(qū)域的CHM 90～100、高程、坡度和坡向。圖2a為研究地冠層高度融合的分布圖，由藍(lán)色到黃色依次代表冠層高度由低到高，其中紅色的點(diǎn)連成的線是星載激光雷達(dá)的激光器發(fā)射到地面上的光斑。我們發(fā)現(xiàn)，該地毛竹冠層高度在0～25.8 m之間。圖2b為星載與機(jī)載激光雷達(dá)融合的高程分布圖，該地地勢(shì)中間高，四周低，高程差為達(dá)到264 m，高程主要在580～640 m之間。圖2c為數(shù)據(jù)融合的坡度分布圖，該地中部地勢(shì)較為平坦，南部、東部及東南部坡度較大。圖2d為數(shù)據(jù)融合的坡向分布圖，主要朝向西北方、北方、南方和東南方，其中坡向?yàn)槲鞅焙蜄|南的占比達(dá)到了70%。

圖2 研究地多源數(shù)據(jù)分布圖(a)冠層高度分布圖 (b)高程分布圖 (c)坡度分布圖 (d)坡向分布圖Fig.2 Distribution maps of multi-source data in the field(a) Canopy height distribution map (b) Elevation distribution map(c) Slope gradient distribution map (d) Slope aspect distribution map

2.2 百分位冠層高度匹配

基于機(jī)載激光雷達(dá)與星載激光雷達(dá)融合的數(shù)據(jù)，對(duì)RH 90～100與CHM 90～100進(jìn)行相關(guān)性分析，結(jié)果如表2所示。星載與機(jī)載激光雷達(dá)冠層高度的數(shù)據(jù)匹配效果好，R2在0.75以上的達(dá)到92.6%。其中，RH 90與CHM 100的R2值最小，為0.63；RH 97與CHM 90的R2值最大，為0.90。對(duì)比分析發(fā)現(xiàn)，表格數(shù)據(jù)以RH 97為中心，呈現(xiàn)向兩邊遞減的趨勢(shì)。包含RH 90、CHM 100的組合匹配程度較低，包含RH 96、RH 97、RH 98的組合匹配程度較高。

2.3 特征參數(shù)重要性分析

利用隨機(jī)森林，分析植被覆蓋度、太陽(yáng)高度角、太陽(yáng)方位角、坡度與坡向這5個(gè)特征變量對(duì)毛竹冠層高度精確度的影響。特征變量重要性的大小排序如圖3所示，太陽(yáng)高度角和太陽(yáng)方位角對(duì)星載激光雷達(dá)數(shù)據(jù)精確度影響較大，植被覆蓋度以及坡度也是影響星載激光雷達(dá)數(shù)據(jù)精確度的重要原因，坡向的影響相對(duì)較小。

2.4 大別山毛竹冠層高度的預(yù)估與驗(yàn)證

根據(jù)表2，我們發(fā)現(xiàn)CHM 90和RH 97擬合程度最高，R2達(dá)到了0.90，所以將RH 97以及特征參數(shù)分析的植被覆蓋度(landsat_treecover)，TAN坡度，坡向，太陽(yáng)高度角(solar_elevation)，太陽(yáng)方位角(solar_azimuth)作為自變量，CHM 90作為因變量進(jìn)行逐步回歸分析，建立多個(gè)毛竹冠層高度預(yù)估模型，算出預(yù)估冠層高度與實(shí)測(cè)數(shù)據(jù)之間的決定系數(shù)(R2)和均方根誤差(RMSE)(表3)。結(jié)果表明，在加上植被覆蓋度(landsat_treecover)這個(gè)參數(shù)后，R2由0.63提高到0.75，RMSE由2.39 m減少到1.81 m。加上坡度參數(shù)后，模型方程的精確度又進(jìn)一步提高(R2增加了0.08，RMSE減少了0.19 m)。

同時(shí)，我們將回歸模型預(yù)估出的毛竹冠層高度進(jìn)行驗(yàn)證，如圖4所示。與實(shí)測(cè)數(shù)據(jù)相比較我們發(fā)現(xiàn)，均方根誤差(RMSE)達(dá)到1.33 m，決定系數(shù)(R2)達(dá)到0.96。

表2 星載激光雷達(dá)數(shù)據(jù)與機(jī)載激光雷達(dá)數(shù)據(jù)匹配結(jié)果

表3 毛竹冠層高度擬合方程

圖3 特征變量重要性排序Fig.3 Sorting characteristic parameters by their importance

圖4 實(shí)測(cè)數(shù)據(jù)驗(yàn)證毛竹冠層高度Fig.4 Estimated Moso bamboo canopy heights verified by measured data

3 討論與結(jié)論

激光雷達(dá)數(shù)據(jù)具有較強(qiáng)的觀測(cè)高度參數(shù)提取能力。毛竹是禾本科剛竹屬單軸散生型常綠喬木狀竹類植物，與闊葉樹、針葉樹不同，毛竹冠幅小。與冠幅大的闊葉樹和針葉樹相比，激光雷達(dá)發(fā)射的激光束更易打到地面[8,13]，如何獲取有效數(shù)據(jù)是研究的主要難點(diǎn)。

基于隨機(jī)森林算法分析各特征參數(shù)對(duì)GEDI數(shù)據(jù)反演冠層高度的影響時(shí)，我們發(fā)現(xiàn)除與太陽(yáng)輻射有關(guān)的2個(gè)參數(shù)外，毛竹的覆蓋度參數(shù)對(duì)激光雷達(dá)反演毛竹冠層高度的精準(zhǔn)度影響非常大[26，32]。影響闊葉樹或針葉樹等大冠幅樹種反演精度的最主要的參數(shù)，可能不是覆蓋度，反而是坡度等地形因素，這要求分析變量的重要性部分有所不同[33-35]。研究表明在復(fù)雜地形條件下需要考慮太陽(yáng)高度角、太陽(yáng)方位角、植被覆蓋度等因素，以提高毛竹林冠層高度反演模型的估算精度。機(jī)載激光雷達(dá)空間分辨率為1 m，GEDI星載激光雷達(dá)分辨率約為25 m。為了讓二者數(shù)據(jù)進(jìn)行融合，我們用ArcMap將機(jī)載激光雷達(dá)數(shù)據(jù)下采樣到25 m[35]。在匹配最優(yōu)百分位數(shù)值時(shí)，研究選擇了CHM 90和RH 97這一組合，R2和RMSE分別為0.96和1.33 m。Potapov等[36]利用GEDI數(shù)據(jù)構(gòu)建全球森林冠層高度模型RMSE為6.64 m；Lang等[37]利用GEDI數(shù)據(jù)構(gòu)建全球森林冠層高度模型，該模型產(chǎn)生的全球冠層高度估計(jì)值的預(yù)期RMSE為2.78 m。因此，利用GEDI數(shù)據(jù)可以大尺度精準(zhǔn)的反演冠層高度。

綜上所述，筆者將星載激光雷達(dá)數(shù)據(jù)與機(jī)載激光雷達(dá)數(shù)據(jù)進(jìn)行融合，以毛竹為研究對(duì)象對(duì)其冠層高度進(jìn)行反演，在評(píng)估特征重要性參數(shù)的基礎(chǔ)上，通過(guò)逐步回歸建立冠層高度預(yù)估模型。模型中充分考慮太陽(yáng)方位角、高度角、植被覆蓋度、坡度、坡向等因子，與實(shí)測(cè)數(shù)據(jù)相比較，模型精確度高(R2=0.96，RMSE=1.33 m)。研究結(jié)果為大尺度高精度高分辨率反演毛竹冠層高度提供了技術(shù)支撐。