蔡耀通 林 輝 孫 華 張 猛 龍江平

( 1. 中南林業(yè)科技大學(xué)林業(yè)遙感信息工程研究中心，湖南 長沙 410004；2. 林業(yè)遙感大數(shù)據(jù)與生態(tài)安全湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410004；3. 南方森林資源經(jīng)營與監(jiān)測國家林業(yè)局重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410004)

林分平均高是森林的重要參數(shù)，目前主要采用人工方法獲取，通過遙感技術(shù)獲取森林樹高主要有三種途徑：傳統(tǒng)光學(xué)遙感、激光雷達(dá)（Lidar）和合成孔徑雷達(dá)（SAR）[1-4]。近年來，隨著遙感對地觀測技術(shù)的發(fā)展，在軌SAR衛(wèi)星越來越多，如ALOS-2，RADARSAT-2，COSMO-SkyMed，TanDEM-X（TDX）等。TDX作為工作在X波段的SAR衛(wèi)星，高時間和空間分辨率的特點(diǎn)使其具有獲取高質(zhì)量測繪產(chǎn)品的能力?！耙话l(fā)雙收”技術(shù)的運(yùn)用，使得TDX只需要單軌飛行便能獲得干涉數(shù)據(jù)，同時能消除時間去相干對數(shù)據(jù)質(zhì)量的影響。回顧現(xiàn)有的研究，相對于其他SAR衛(wèi)星數(shù)據(jù)，TDX數(shù)據(jù)被廣泛用于森林樹高的估測研究，并在研究中取得了不錯的效果[5-7]。在植被覆蓋的林區(qū)，高頻波段（如X波段）難以穿透植被到達(dá)地表，此時利用TDX能獲得涵蓋森林樹高信息的數(shù)字表面模型（DSM），之后結(jié)合DEM能準(zhǔn)確提取森林樹高，因此TDX數(shù)據(jù)具備了獲取高精度森林樹高的能力。

2011年，TDX單軌干涉數(shù)據(jù)開始提供，國內(nèi)外學(xué)者開始利用TDX數(shù)據(jù)采用InSAR和POLIn-SAR方法對森林樹高進(jìn)行反演研究[8-12]。相對而言，DSM-DEM差分法處理流程簡明，原理淺顯且算法簡單。隨著TDX數(shù)據(jù)的積累，該法在森林樹高的反演研究中也被大量應(yīng)用。然而SAR微波具有較強(qiáng)的穿透性，容易穿透林區(qū)中植被間的間隙（低郁閉度林分）到達(dá)地表使得相位中心降低而無法得到準(zhǔn)確的植被層相位中心高度。且基于DSM-DEM差分法的樹高反演通常以林分為單位進(jìn)行，反演過程容易混合了林隙地面高程與林分冠層高度的影響而使林分平均高反演結(jié)果趨于偏低。因此，傳統(tǒng)的DSM-DEM差分法并不能取得很好的反演結(jié)果。為了解決這一問題，大多研究利用差分法獲得的冠層高度模型（CHM）進(jìn)行回歸建模反演林分平均高[5]。但是這類解決方法局限性較大，回歸模型地域性強(qiáng)，普適性差，無法進(jìn)行大范圍推廣?；诖耍狙芯恳?對經(jīng)配準(zhǔn)的TanDEM-X /TerraSAR-X HH單極化干涉對和Lidar DEM為數(shù)據(jù)源，首先對干涉對進(jìn)行干涉處理獲得DSM，其次運(yùn)用DSM-DEM差分法獲取試驗(yàn)區(qū)CHM，以完全約束最小二乘法混合像元分解獲得試驗(yàn)區(qū)植被豐度并對CHM進(jìn)行校正，從而實(shí)現(xiàn)林分平均高反演，為結(jié)合多源遙感數(shù)據(jù)提高混合像元效應(yīng)下的林分平均高反演精度提供了行之有效的方案。

1 基于TanDEM-X數(shù)據(jù)的林分平均高反演方法

1.1 技術(shù)路線

利用2對TDX干涉對數(shù)據(jù)通過InSAR處理獲得試驗(yàn)區(qū)DSM，采用DSM-DEM差分法結(jié)合Lidar DEM得到CHM。為了降低傳統(tǒng)的DSMDEM差分法由于混合像元效應(yīng)及SAR微波對植被的穿透性對CHM估算結(jié)果的影響，本研究以TDX干涉對、GF-2影像和Lidar點(diǎn)云數(shù)據(jù)為基礎(chǔ)，在極化干涉技術(shù)下利用TDX干涉對生成覆蓋研究區(qū)的DSM，并聯(lián)合Lidar DEM數(shù)據(jù)獲取研究區(qū)CHM。隨后基于高分辨率GF-2影像采用完全約束最小二乘混合像元分解方法獲得試驗(yàn)區(qū)植被豐度，以此對樣地內(nèi)CHM估算結(jié)果進(jìn)行提取與校正，得到樣地林分平均高，結(jié)合地面調(diào)查數(shù)據(jù)進(jìn)行精度驗(yàn)證。具體技術(shù)路線見圖1。

圖 1 技術(shù)路線Fig. 1 Technical route

1.2 植被豐度反演

遙感影像通過純凈端元提取和豐度估計(jì)，能在像元尺度上表達(dá)各地物的光譜反射貢獻(xiàn)率，得到各類端元所占比例（豐度），從而獲得像元內(nèi)林分冠層和其他成分所占比例[13]，進(jìn)而對CHM估算結(jié)果進(jìn)行校正?；贕F-2影像采用完全約束最小二乘混合像元分解方法獲得試驗(yàn)區(qū)植被豐度圖，并對CHM估算結(jié)果進(jìn)行校正。首先通過最小噪聲分離進(jìn)行主成分分析并估計(jì)影像噪聲，根據(jù)噪聲估計(jì)結(jié)果計(jì)算純凈像元指數(shù)（PPI），其中PPI計(jì)算的迭代次數(shù)、迭代單位、閾值系數(shù)分別設(shè)置為10 000、250、2.50。隨后，調(diào)用n維可視化工具（n-D Visualizer）獲取各地類純凈像元，采用完全約束最小二乘法進(jìn)行混合像元分解反演植被豐度。最后，在試驗(yàn)區(qū)中生成150個隨機(jī)點(diǎn)，利用同期Google Earth影像（分辨率為0.3 m）混合像元分解結(jié)果對GF-2影像植被豐度進(jìn)行精度檢驗(yàn)。

1.3 林分平均高估測方法

DSM-DEM差分法通過帶有植被高度的DSM和DEM差分獲得CHM[14]，本研究中DSM來自于TDX數(shù)據(jù)干涉處理結(jié)果，DEM通過Lidar點(diǎn)云數(shù)據(jù)獲取[15-16]。傳統(tǒng)差分法通常利用式（1）獲得林分冠層高度（H），然而由于混合像元效應(yīng)的存在，H常常包含了地表高度的貢獻(xiàn)，因此林分平均高會有不同程度的過低估測。式（2）～（3）利用植被豐度值校正CHM估計(jì)值后，得到去除地表高度貢獻(xiàn)后的林分冠層高度（H'），根據(jù)實(shí)際情況剔除CHM中大于35 m和小于0 m數(shù)據(jù)，并按林分位置矢量邊界提取校正值作為林分平均高。

式中：H和H'分別表示差分法獲得的林分冠層高度和校正后的林分冠層高度，DNi和DN′i分別表示CHM、校正后CHM的第i個像元 表示林分范圍內(nèi)所包含的像元個數(shù)，VAi表示第i個像元的植被豐度。

1.4 精度驗(yàn)證

反演結(jié)果主要通過計(jì)算其估測值（林分平均高反演值）與觀測值（外業(yè)數(shù)據(jù)）之間的決定系數(shù)（R2），均方根誤差（RMSE）以及估測精度（EA）3個指標(biāo)進(jìn)行精度驗(yàn)證。

式中： 和yi分別表示第i個樣地的估測與實(shí)測林分平均高； 為所有實(shí)測林分平均高的算術(shù)平均值；n為樣地?cái)?shù)量。

2 數(shù)據(jù)來源及處理

2.1 外業(yè)數(shù)據(jù)

選擇位于內(nèi)蒙古赤峰市喀喇沁旗西南部的旺業(yè)甸國有林場（118°15′E，41°40′N）為試驗(yàn)區(qū)。外業(yè)數(shù)據(jù)采集開始日期為2017年9月20日，歷時10 d。研究區(qū)內(nèi)共設(shè)置樣地38塊，其中油松（Pinus tabulaeformis）樣地21塊，落葉松（Larix gmelinii）樣地17塊，具體分布見圖2。外業(yè)調(diào)查在無人機(jī)飛行區(qū)域內(nèi)按隨機(jī)抽樣方法設(shè)置集中樣地，樣地大小為25 m×25 m，每塊樣地內(nèi)林木均為相同起源、林齡相近的人工純林。在樣地內(nèi)開展每木檢尺，主要觀測因子有：樹種、樹高、胸徑、年齡、枝下高、冠幅（東西向、南北向）。樣地林分平均高分布情況見圖3。由圖3可知，樣地林分平均高大多為10～20 m，油松樣地林分平均高總體大于落葉松樣地林分平均高。其中油松樣地林分平均高分布在11.7～19.7 m，平均值16.3 m，標(biāo)準(zhǔn)差2.0 m；落葉松樣地林分平均高分布在6.7～17.3 m，平均值12.1 m，標(biāo)準(zhǔn)差3.1 m。樣地林分平均高分布符合林場內(nèi)的林木樹高分布情況，因此本研究樣本具有較好的代表性。

2.2 機(jī)載 Lidar 數(shù)據(jù)

無人機(jī)機(jī)載Lidar點(diǎn)云數(shù)據(jù)于2017年9月22—24日獲取，分布在林場3個區(qū)域共計(jì)面積約5 km2，覆蓋全部38個精細(xì)調(diào)查樣地。使用Terrasolid軟件對點(diǎn)云數(shù)據(jù)進(jìn)行噪聲剔除及地面點(diǎn)、非地面點(diǎn)分類等預(yù)處理，生成1 m分辨率的DEM。采用與InSAR DSM一致的WGS84基準(zhǔn)面，將Lidar DEM重采樣至5 m分辨率后與DSM進(jìn)行配準(zhǔn)，以便后續(xù)的差分處理。

2.3 TDX數(shù)據(jù)

SAR數(shù)據(jù)中包括2對條帶成像模式下獲得的TDX和TerraSAR影像構(gòu)成的HH單極化干涉對，干涉對覆蓋范圍約為32 km×56 km，距離向分辨率為2.4 m，方位向分辨率為3.3 m。2干涉對分別于2014年1月9日和2014年1月31日獲取，數(shù)據(jù)獲取當(dāng)天天氣狀況良好，中心入射角分別為34.8°和37.1°，基線適合反演森林結(jié)構(gòu)參數(shù)（高程模糊度大于試驗(yàn)區(qū)最大樹高）。干涉對具體參數(shù)見表1。

表 1 TanDEM-X 干涉對參數(shù)Table 1 Parameters of TanDEM-X interference pairs

圖 2 試驗(yàn)區(qū)及樣地位置Fig. 2 The study area and location of samples site

圖 3 樣地林分平均高分布Fig. 3 Stand mean height distribution of sample plots

TDX數(shù)據(jù)主要用于DSM的提取。DSM的獲取基于InSAR測高原理，當(dāng)天線向地面發(fā)射并接收電磁波時，若2衛(wèi)星天線接收的回波信號具有高相干性，那么信號的傳播路徑差可以通過信號間的相位差與干涉測量間的幾何關(guān)系獲得[17-18]。

對TDX數(shù)據(jù)進(jìn)行干涉處理，主要流程包括數(shù)據(jù)導(dǎo)入生成SLC數(shù)據(jù)對、基線估計(jì)、干涉圖生成與去平地效應(yīng)、Goldstein濾波與相干性生成、相位解纏、軌道精煉與重去平、相位轉(zhuǎn)高程與地理編碼。所有處理流程均在SARscape 5.3.1中完成。

2.4 DEM數(shù)據(jù)

在InSAR處理中，參考DEM在相位解纏、地理編碼和地形補(bǔ)償三部分被應(yīng)用，因此參考DEM的好壞將直接影響InSAR處理的質(zhì)量。研究中使用的30 m分辨率參考DEM（ASTER GDEM V2）數(shù)據(jù)來源于中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺（http：//www.gscloud.cn）。其垂直精度約20 m，能滿足InSAR處理的要求。

2.5 高分二號影像數(shù)據(jù)

高分二號（GF-2）具有2 m空間分辨率及45 km成像幅寬能力，是我國迄今為止空間分辨率最高的民用衛(wèi)星，其具有的4個多光譜波段能夠用于準(zhǔn)確、有效地估計(jì)林區(qū)植被豐度。本研究利用覆蓋試驗(yàn)區(qū)的4景GF-2影像，其獲取日期為2017年9月5日。對GF-2影像的預(yù)處理包括輻射定標(biāo)、FLAASH大氣校正、正射校正、鑲嵌、裁剪與配準(zhǔn)，最后將試驗(yàn)區(qū)影像重采樣至5 m分辨率，所有預(yù)處理均在ENVI 5.3軟件中完成。

3 結(jié)果與分析

3.1 提取DSM

InSAR處理獲得的DSM見圖4。由圖4可知，2景InSAR DSM地表細(xì)節(jié)顯示清楚，處理結(jié)果良好，其范圍覆蓋整個旺業(yè)甸林場，可用于精確反演該區(qū)域林分平均高等森林結(jié)構(gòu)參數(shù)。其中，DSM1高程分布在710～1 864 m，DSM2高程分布在650～1 829 m。

圖 4 干涉處理獲得的DSMFig. 4 DSM obtained by interference processing

干涉相干性統(tǒng)計(jì)圖及相應(yīng)的InSAR DSM理論誤差見圖5。圖5a、c橫軸表示干涉對相干性系數(shù)，縱軸表示像元個數(shù)。TDX數(shù)據(jù)沒有受到時間去相干的影響，干涉對相干性估計(jì)情況較好。統(tǒng)計(jì)干涉對相干性，干涉對1（2014-01-09）和干涉對2（2014-01-31）相干性系數(shù)均值分別為0.81、0.82。圖4b、d，橫軸為相干性，縱軸為InSAR DSM理論精度。從圖5可知，2景InSAR DSM在誤差允許范圍內(nèi)，其理論精度為2 m，基本滿足反演條件。

圖 5 TDX數(shù)據(jù)相干性及與InSAR理論精度關(guān)系Fig. 5 TDX data coherence and its relationship with InSAR theoretical accuracy

3.2 提取 DEM

本研究中獲得的Lidar DEM原始分辨率為1 m，而后為便于DSM與DEM的差分，將DEM重采樣至5 m分辨率，圖6展示了A區(qū)、B區(qū)及C區(qū)的5 m分辨率Lidar DEM。去除點(diǎn)云密度不足區(qū)域后，Lidar DEM置信度在90%以上，滿足差分法反演林分平均高要求。

圖 6 重采樣后Lidar DEMFig. 6 Lidar DEM after resampling

3.3 植被豐度反演結(jié)果

本研究基于GF-2影像采用完全約束最小二乘混合像元分解方法獲得試驗(yàn)區(qū)植被豐度，見圖7。由圖7可知，植被豐度圖整體性良好，影像無云覆蓋區(qū)域符合林場內(nèi)植被分布實(shí)際情況。利用Google Earth同期影像反演的植被豐度對試驗(yàn)區(qū)內(nèi)150個隨機(jī)點(diǎn)的GF-2植被豐度進(jìn)行精度驗(yàn)證，其精度較好，總體精度為84%。GF-2植被豐度能在像元尺度上表達(dá)植被與其他地類的光譜貢獻(xiàn)比例，并用于校正差分法得到的CHM。

圖 7 植被豐度圖Fig. 7 Vegetation abundance map

3.4 林分平均高反演結(jié)果

分別對傳統(tǒng)差分法和本研究方法對林分平均高實(shí)測值與估測值進(jìn)行統(tǒng)計(jì)分析，使用3個評價(jià)指標(biāo)對2個樹種樣地林分平均高反演結(jié)果進(jìn)行精度驗(yàn)證，結(jié)果見表2。由表2可知，對傳統(tǒng)DSM-DEM差分法林分平均高估測值與實(shí)測值進(jìn)行線性擬合，R2為0.24?？紤]到偏高與偏低部分估測值偏差較大，可能受到了林分低郁閉度所產(chǎn)生的混合像元效應(yīng)影響[19]。按本研究方法校正冠層高度模型后，林分平均高估測值與實(shí)測值間偏差相比剔除前減小明顯。

表 2 林分平均高反演精度驗(yàn)證Table 2 Predicted accuracy validation of stand mean heights

可見，傳統(tǒng)差分法中油松、落葉松樣地林分平均高反演結(jié)果并不太理想，兩者RMSE分別為6.8 m和4.2 m，EA分別為53.2%和63.3%。CHM經(jīng)過校正后，R2有所升高，表明了估測值和實(shí)測值間有更好的相關(guān)性。以植被豐度為校正因子能有效剔除林地中由于郁閉度所產(chǎn)生的混合像元作用和噪聲等因素對林分平均高反演的影響，經(jīng)過校正后的冠層高度模型估測精度有了較大的提高。由此表明，植被豐度能表達(dá)像元尺度上的植被光譜貢獻(xiàn)比例，在一定程度上能夠減小由于SAR微波對植被的穿透性、林分低郁閉度和噪聲對差分法反演林分平均高的影響，從而提高林分平均高反演精度。

4 結(jié)論與討論

本研究以內(nèi)蒙古旺業(yè)甸國有林場為試驗(yàn)區(qū)，結(jié)合TDX和Lidar DEM數(shù)據(jù)使用DSM-DEM差分法估計(jì)CHM，利用完全約束最小二乘混合像元分解反演的植被豐度作為校正因子，以樣地為單位校正并統(tǒng)計(jì)林分CHM平均值，實(shí)現(xiàn)林分平均高反演。

1）本研究利用TDX數(shù)據(jù)，采用DSM-DEM差分法反演了林分平均高，并獲得了較好的反演結(jié)果，顯示了基于TDX數(shù)據(jù)的DSM-DEM差分法在林分平均高反演上具有較大的潛力。

2）使用植被豐度校正CHM后，林分平均高的估測精度大幅度提高，反映了以像元植被豐度為校正因子能有效剔除林地中由于郁閉度所產(chǎn)生的混合像元作用和噪聲等因素對林分平均高反演的影響。

研究中由于TDX數(shù)據(jù)（2014年1月）與外業(yè)數(shù)據(jù)（2017年9月）獲取時間不同步，忽略了數(shù)據(jù)獲取時間差間所產(chǎn)生的樹高生長；而X波段電磁波對森林冠層具有一定的穿透作用，也將影響了對林分平均高的估測[20]。因此，R2并沒有達(dá)到更高的水平。在像元尺度上利用植被豐度校正林分冠層高度模型，在一定程度上能提高林分平均高的估測精度。下一步研究將結(jié)合面向?qū)ο蟮幕旌舷裨纸饧夹g(shù)，強(qiáng)化本研究方法在景觀異質(zhì)性較大的天然林中的魯棒性。

[ 參 考 文 獻(xiàn) ]

[1]廖明生, 王騰. 時間序列InSAR技術(shù)與應(yīng)用[M]. 北京： 科學(xué)出版社, 2014.

[2]王超, 張紅, 劉智. 星載合成孔徑雷達(dá)干涉測量[M].北京： 科學(xué)出版社, 2002.

[3]張曉莉, 趙鵬祥, 高凌寒, 等. 基于ArboLiDAR的林分自動分割研究與應(yīng)用 [J]. 中南林業(yè)科技大學(xué)學(xué)報(bào),2017, 37(11)： 76-83.

[4]孫華, 鞠洪波, 張懷清, 等. 基于 Worldview-2影像的林木冠幅提取與樹高反演 [J]. 中南林業(yè)科技大學(xué)學(xué)報(bào), 2014, 34(10)： 45-50.

[5]張王菲, 陳爾學(xué), 李增元, 等. 干涉、極化干涉SAR技術(shù)森林高度估測算法研究進(jìn)展 [J]. 遙感技術(shù)與應(yīng)用, 2017, 32(6)： 983-997.

[6]岳彩榮, 肖虹雁, 曹霸. 基于PolInSAR森林高度反演研究 [J]. 西南林業(yè)大學(xué)學(xué)報(bào), 2016, 36(3)： 137-143.

[7]章皖秋, 岳彩榮, 劉曉英. 基于TerraSAR-X/TanDEMX干涉DEM的森林冠層高度估測 [J]. 西南林業(yè)大學(xué)學(xué)報(bào), 2016, 36(6)： 64-72.

[8]Kugler F, Schulze D, Hajnsek I, et al. TanDEM-X pol-InSAR performance for forest height estimation [J].IEEE Transactions on Geoscience and Remote Sensing,2014, 52(10)： 6404-6422.

[9]Neumann M, Ferro-Famil L, Reigber A. Estimation of forest structure, ground, and canopy layer characteristics from multibaseline polarimetric interferometric SAR data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(3)： 1086-1104.

[10]談璐璐, 楊立波, 楊汝良. 基于ESPRIT算法的極化干涉SAR植被高度反演研究 [J]. 測繪學(xué)報(bào), 2011,40(3)： 296-300.

[11]杜凱, 林輝, 龍江平, 等. 基于GVB模型改進(jìn)算法的PolInSAR林分高度反演 [J]. 中南林業(yè)科技大學(xué)學(xué)報(bào), 2018, 38(1)： 49-54.

[12]Karila K, Vastaranta M, Karjalainen M, et al. Tandem-X interferometry in the prediction of forest inventory attributes in managed boreal forests [J]. Remote Sensing of Environment, 2015, 159： 259-268.

[13]藍(lán)金輝, 鄒金霖, 郝彥爽, 等. 高光譜遙感影像混合像元分解研究進(jìn)展 [J]. 遙感學(xué)報(bào), 2018, 22(1)： 13-27.

[14]龐勇, 李增元, 陳爾學(xué), 等. 干涉雷達(dá)技術(shù)用于林分高估測 [J]. 遙感學(xué)報(bào), 2003, 7(1)： 8-13.

[15]Kenyi L W, Dubayah R, Hofton M, et al. Comparative analysis of SRTM-NED vegetation canopy height to LIDAR-derived vegetation canopy metrics [J]. International Journal of Remote Sensing, 2009, 30(11)：2797-2811.

[16]Rignot E. Dual-frequency interferometric SAR observations of a tropical rain-forest [J]. Geophysical Research Letters, 1996, 23(9)： 993-996.

[17]Cloude S R. 極化建模與雷達(dá)遙感應(yīng)用[M]. 北京： 電子工業(yè)出版社, 2015.

[18]Woodhouse I H. 微波遙感導(dǎo)論[M]. 董曉龍, 徐星歐,徐曦煜, 譯. 北京： 科學(xué)出版社, 2014.

[19]Miliaresis G, Delikaraoglou D. Effects of percent tree canopy density and DEM misregistration on SRTM/NED vegetation height estimates [J]. Remote Sensing, 2009, 1(2)： 36-49.

[20]Praks J, Kugler F, Papathanassiou K P, et al. Height estimation of boreal forest： interferometric model-based inversion at L- and X-band versus HUTSCAT profiling scatterometer [J]. IEEE Geoscience & Remote Sensing Letters, 2007, 4(3)： 466-470.