蔡耀通 林 輝 孫 華 張 猛 龍江平
( 1. 中南林業科技大學林業遙感信息工程研究中心,湖南 長沙 410004;2. 林業遙感大數據與生態安全湖南省重點實驗室,湖南 長沙 410004;3. 南方森林資源經營與監測國家林業局重點實驗室,湖南 長沙 410004)
林分平均高是森林的重要參數,目前主要采用人工方法獲取,通過遙感技術獲取森林樹高主要有三種途徑:傳統光學遙感、激光雷達(Lidar)和合成孔徑雷達(SAR)[1-4]。近年來,隨著遙感對地觀測技術的發展,在軌SAR衛星越來越多,如ALOS-2,RADARSAT-2,COSMO-SkyMed,TanDEM-X(TDX)等。TDX作為工作在X波段的SAR衛星,高時間和空間分辨率的特點使其具有獲取高質量測繪產品的能力。“一發雙收”技術的運用,使得TDX只需要單軌飛行便能獲得干涉數據,同時能消除時間去相干對數據質量的影響。回顧現有的研究,相對于其他SAR衛星數據,TDX數據被廣泛用于森林樹高的估測研究,并在研究中取得了不錯的效果[5-7]。在植被覆蓋的林區,高頻波段(如X波段)難以穿透植被到達地表,此時利用TDX能獲得涵蓋森林樹高信息的數字表面模型(DSM),之后結合DEM能準確提取森林樹高,因此TDX數據具備了獲取高精度森林樹高的能力。
2011年,TDX單軌干涉數據開始提供,國內外學者開始利用TDX數據采用InSAR和POLIn-SAR方法對森林樹高進行反演研究[8-12]。相對而言,DSM-DEM差分法處理流程簡明,原理淺顯且算法簡單。隨著TDX數據的積累,該法在森林樹高的反演研究中也被大量應用。然而SAR微波具有較強的穿透性,容易穿透林區中植被間的間隙(低郁閉度林分)到達地表使得相位中心降低而無法得到準確的植被層相位中心高度。且基于DSM-DEM差分法的樹高反演通常以林分為單位進行,反演過程容易混合了林隙地面高程與林分冠層高度的影響而使林分平均高反演結果趨于偏低。因此,傳統的DSM-DEM差分法并不能取得很好的反演結果。為了解決這一問題,大多研究利用差分法獲得的冠層高度模型(CHM)進行回歸建模反演林分平均高[5]。但是這類解決方法局限性較大,回歸模型地域性強,普適性差,無法進行大范圍推廣。基于此,本研究以2對經配準的TanDEM-X /TerraSAR-X HH單極化干涉對和Lidar DEM為數據源,首先對干涉對進行干涉處理獲得DSM,其次運用DSM-DEM差分法獲取試驗區CHM,以完全約束最小二乘法混合像元分解獲得試驗區植被豐度并對CHM進行校正,從而實現林分平均高反演,為結合多源遙感數據提高混合像元效應下的林分平均高反演精度提供了行之有效的方案。
利用2對TDX干涉對數據通過InSAR處理獲得試驗區DSM,采用DSM-DEM差分法結合Lidar DEM得到CHM。為了降低傳統的DSMDEM差分法由于混合像元效應及SAR微波對植被的穿透性對CHM估算結果的影響,本研究以TDX干涉對、GF-2影像和Lidar點云數據為基礎,在極化干涉技術下利用TDX干涉對生成覆蓋研究區的DSM,并聯合Lidar DEM數據獲取研究區CHM。隨后基于高分辨率GF-2影像采用完全約束最小二乘混合像元分解方法獲得試驗區植被豐度,以此對樣地內CHM估算結果進行提取與校正,得到樣地林分平均高,結合地面調查數據進行精度驗證。具體技術路線見圖1。

圖 1 技術路線Fig. 1 Technical route
遙感影像通過純凈端元提取和豐度估計,能在像元尺度上表達各地物的光譜反射貢獻率,得到各類端元所占比例(豐度),從而獲得像元內林分冠層和其他成分所占比例[13],進而對CHM估算結果進行校正。基于GF-2影像采用完全約束最小二乘混合像元分解方法獲得試驗區植被豐度圖,并對CHM估算結果進行校正。首先通過最小噪聲分離進行主成分分析并估計影像噪聲,根據噪聲估計結果計算純凈像元指數(PPI),其中PPI計算的迭代次數、迭代單位、閾值系數分別設置為10 000、250、2.50。隨后,調用n維可視化工具(n-D Visualizer)獲取各地類純凈像元,采用完全約束最小二乘法進行混合像元分解反演植被豐度。最后,在試驗區中生成150個隨機點,利用同期Google Earth影像(分辨率為0.3 m)混合像元分解結果對GF-2影像植被豐度進行精度檢驗。
DSM-DEM差分法通過帶有植被高度的DSM和DEM差分獲得CHM[14],本研究中DSM來自于TDX數據干涉處理結果,DEM通過Lidar點云數據獲取[15-16]。傳統差分法通常利用式(1)獲得林分冠層高度(H),然而由于混合像元效應的存在,H常常包含了地表高度的貢獻,因此林分平均高會有不同程度的過低估測。式(2)~(3)利用植被豐度值校正CHM估計值后,得到去除地表高度貢獻后的林分冠層高度(H'),根據實際情況剔除CHM中大于35 m和小于0 m數據,并按林分位置矢量邊界提取校正值作為林分平均高。

式中:H和H'分別表示差分法獲得的林分冠層高度和校正后的林分冠層高度,DNi和DN′i分別表示CHM、校正后CHM的第i個像元 表示林分范圍內所包含的像元個數,VAi表示第i個像元的植被豐度。
反演結果主要通過計算其估測值(林分平均高反演值)與觀測值(外業數據)之間的決定系數(R2),均方根誤差(RMSE)以及估測精度(EA)3個指標進行精度驗證。
式中: 和yi分別表示第i個樣地的估測與實測林分平均高; 為所有實測林分平均高的算術平均值;n為樣地數量。
選擇位于內蒙古赤峰市喀喇沁旗西南部的旺業甸國有林場(118°15′E,41°40′N)為試驗區。外業數據采集開始日期為2017年9月20日,歷時10 d。研究區內共設置樣地38塊,其中油松(Pinus tabulaeformis)樣地21塊,落葉松(Larix gmelinii)樣地17塊,具體分布見圖2。外業調查在無人機飛行區域內按隨機抽樣方法設置集中樣地,樣地大小為25 m×25 m,每塊樣地內林木均為相同起源、林齡相近的人工純林。在樣地內開展每木檢尺,主要觀測因子有:樹種、樹高、胸徑、年齡、枝下高、冠幅(東西向、南北向)。樣地林分平均高分布情況見圖3。由圖3可知,樣地林分平均高大多為10~20 m,油松樣地林分平均高總體大于落葉松樣地林分平均高。其中油松樣地林分平均高分布在11.7~19.7 m,平均值16.3 m,標準差2.0 m;落葉松樣地林分平均高分布在6.7~17.3 m,平均值12.1 m,標準差3.1 m。樣地林分平均高分布符合林場內的林木樹高分布情況,因此本研究樣本具有較好的代表性。
無人機機載Lidar點云數據于2017年9月22—24日獲取,分布在林場3個區域共計面積約5 km2,覆蓋全部38個精細調查樣地。使用Terrasolid軟件對點云數據進行噪聲剔除及地面點、非地面點分類等預處理,生成1 m分辨率的DEM。采用與InSAR DSM一致的WGS84基準面,將Lidar DEM重采樣至5 m分辨率后與DSM進行配準,以便后續的差分處理。
SAR數據中包括2對條帶成像模式下獲得的TDX和TerraSAR影像構成的HH單極化干涉對,干涉對覆蓋范圍約為32 km×56 km,距離向分辨率為2.4 m,方位向分辨率為3.3 m。2干涉對分別于2014年1月9日和2014年1月31日獲取,數據獲取當天天氣狀況良好,中心入射角分別為34.8°和37.1°,基線適合反演森林結構參數(高程模糊度大于試驗區最大樹高)。干涉對具體參數見表1。

表 1 TanDEM-X 干涉對參數Table 1 Parameters of TanDEM-X interference pairs

圖 2 試驗區及樣地位置Fig. 2 The study area and location of samples site

圖 3 樣地林分平均高分布Fig. 3 Stand mean height distribution of sample plots
TDX數據主要用于DSM的提取。DSM的獲取基于InSAR測高原理,當天線向地面發射并接收電磁波時,若2衛星天線接收的回波信號具有高相干性,那么信號的傳播路徑差可以通過信號間的相位差與干涉測量間的幾何關系獲得[17-18]。
對TDX數據進行干涉處理,主要流程包括數據導入生成SLC數據對、基線估計、干涉圖生成與去平地效應、Goldstein濾波與相干性生成、相位解纏、軌道精煉與重去平、相位轉高程與地理編碼。所有處理流程均在SARscape 5.3.1中完成。
在InSAR處理中,參考DEM在相位解纏、地理編碼和地形補償三部分被應用,因此參考DEM的好壞將直接影響InSAR處理的質量。研究中使用的30 m分辨率參考DEM(ASTER GDEM V2)數據來源于中國科學院計算機網絡信息中心地理空間數據云平臺(http://www.gscloud.cn)。其垂直精度約20 m,能滿足InSAR處理的要求。
高分二號(GF-2)具有2 m空間分辨率及45 km成像幅寬能力,是我國迄今為止空間分辨率最高的民用衛星,其具有的4個多光譜波段能夠用于準確、有效地估計林區植被豐度。本研究利用覆蓋試驗區的4景GF-2影像,其獲取日期為2017年9月5日。對GF-2影像的預處理包括輻射定標、FLAASH大氣校正、正射校正、鑲嵌、裁剪與配準,最后將試驗區影像重采樣至5 m分辨率,所有預處理均在ENVI 5.3軟件中完成。
InSAR處理獲得的DSM見圖4。由圖4可知,2景InSAR DSM地表細節顯示清楚,處理結果良好,其范圍覆蓋整個旺業甸林場,可用于精確反演該區域林分平均高等森林結構參數。其中,DSM1高程分布在710~1 864 m,DSM2高程分布在650~1 829 m。

圖 4 干涉處理獲得的DSMFig. 4 DSM obtained by interference processing
干涉相干性統計圖及相應的InSAR DSM理論誤差見圖5。圖5a、c橫軸表示干涉對相干性系數,縱軸表示像元個數。TDX數據沒有受到時間去相干的影響,干涉對相干性估計情況較好。統計干涉對相干性,干涉對1(2014-01-09)和干涉對2(2014-01-31)相干性系數均值分別為0.81、0.82。圖4b、d,橫軸為相干性,縱軸為InSAR DSM理論精度。從圖5可知,2景InSAR DSM在誤差允許范圍內,其理論精度為2 m,基本滿足反演條件。

圖 5 TDX數據相干性及與InSAR理論精度關系Fig. 5 TDX data coherence and its relationship with InSAR theoretical accuracy
本研究中獲得的Lidar DEM原始分辨率為1 m,而后為便于DSM與DEM的差分,將DEM重采樣至5 m分辨率,圖6展示了A區、B區及C區的5 m分辨率Lidar DEM。去除點云密度不足區域后,Lidar DEM置信度在90%以上,滿足差分法反演林分平均高要求。

圖 6 重采樣后Lidar DEMFig. 6 Lidar DEM after resampling
本研究基于GF-2影像采用完全約束最小二乘混合像元分解方法獲得試驗區植被豐度,見圖7。由圖7可知,植被豐度圖整體性良好,影像無云覆蓋區域符合林場內植被分布實際情況。利用Google Earth同期影像反演的植被豐度對試驗區內150個隨機點的GF-2植被豐度進行精度驗證,其精度較好,總體精度為84%。GF-2植被豐度能在像元尺度上表達植被與其他地類的光譜貢獻比例,并用于校正差分法得到的CHM。

圖 7 植被豐度圖Fig. 7 Vegetation abundance map
分別對傳統差分法和本研究方法對林分平均高實測值與估測值進行統計分析,使用3個評價指標對2個樹種樣地林分平均高反演結果進行精度驗證,結果見表2。由表2可知,對傳統DSM-DEM差分法林分平均高估測值與實測值進行線性擬合,R2為0.24。考慮到偏高與偏低部分估測值偏差較大,可能受到了林分低郁閉度所產生的混合像元效應影響[19]。按本研究方法校正冠層高度模型后,林分平均高估測值與實測值間偏差相比剔除前減小明顯。

表 2 林分平均高反演精度驗證Table 2 Predicted accuracy validation of stand mean heights
可見,傳統差分法中油松、落葉松樣地林分平均高反演結果并不太理想,兩者RMSE分別為6.8 m和4.2 m,EA分別為53.2%和63.3%。CHM經過校正后,R2有所升高,表明了估測值和實測值間有更好的相關性。以植被豐度為校正因子能有效剔除林地中由于郁閉度所產生的混合像元作用和噪聲等因素對林分平均高反演的影響,經過校正后的冠層高度模型估測精度有了較大的提高。由此表明,植被豐度能表達像元尺度上的植被光譜貢獻比例,在一定程度上能夠減小由于SAR微波對植被的穿透性、林分低郁閉度和噪聲對差分法反演林分平均高的影響,從而提高林分平均高反演精度。
本研究以內蒙古旺業甸國有林場為試驗區,結合TDX和Lidar DEM數據使用DSM-DEM差分法估計CHM,利用完全約束最小二乘混合像元分解反演的植被豐度作為校正因子,以樣地為單位校正并統計林分CHM平均值,實現林分平均高反演。
1)本研究利用TDX數據,采用DSM-DEM差分法反演了林分平均高,并獲得了較好的反演結果,顯示了基于TDX數據的DSM-DEM差分法在林分平均高反演上具有較大的潛力。
2)使用植被豐度校正CHM后,林分平均高的估測精度大幅度提高,反映了以像元植被豐度為校正因子能有效剔除林地中由于郁閉度所產生的混合像元作用和噪聲等因素對林分平均高反演的影響。
研究中由于TDX數據(2014年1月)與外業數據(2017年9月)獲取時間不同步,忽略了數據獲取時間差間所產生的樹高生長;而X波段電磁波對森林冠層具有一定的穿透作用,也將影響了對林分平均高的估測[20]。因此,R2并沒有達到更高的水平。在像元尺度上利用植被豐度校正林分冠層高度模型,在一定程度上能提高林分平均高的估測精度。下一步研究將結合面向對象的混合像元分解技術,強化本研究方法在景觀異質性較大的天然林中的魯棒性。
[ 參 考 文 獻 ]
[1]廖明生, 王騰. 時間序列InSAR技術與應用[M]. 北京: 科學出版社, 2014.
[2]王超, 張紅, 劉智. 星載合成孔徑雷達干涉測量[M].北京: 科學出版社, 2002.
[3]張曉莉, 趙鵬祥, 高凌寒, 等. 基于ArboLiDAR的林分自動分割研究與應用 [J]. 中南林業科技大學學報,2017, 37(11): 76-83.
[4]孫華, 鞠洪波, 張懷清, 等. 基于 Worldview-2影像的林木冠幅提取與樹高反演 [J]. 中南林業科技大學學報, 2014, 34(10): 45-50.
[5]張王菲, 陳爾學, 李增元, 等. 干涉、極化干涉SAR技術森林高度估測算法研究進展 [J]. 遙感技術與應用, 2017, 32(6): 983-997.
[6]岳彩榮, 肖虹雁, 曹霸. 基于PolInSAR森林高度反演研究 [J]. 西南林業大學學報, 2016, 36(3): 137-143.
[7]章皖秋, 岳彩榮, 劉曉英. 基于TerraSAR-X/TanDEMX干涉DEM的森林冠層高度估測 [J]. 西南林業大學學報, 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]. 測繪學報, 2011,40(3): 296-300.
[11]杜凱, 林輝, 龍江平, 等. 基于GVB模型改進算法的PolInSAR林分高度反演 [J]. 中南林業科技大學學報, 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]藍金輝, 鄒金霖, 郝彥爽, 等. 高光譜遙感影像混合像元分解研究進展 [J]. 遙感學報, 2018, 22(1): 13-27.
[14]龐勇, 李增元, 陳爾學, 等. 干涉雷達技術用于林分高估測 [J]. 遙感學報, 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. 極化建模與雷達遙感應用[M]. 北京: 電子工業出版社, 2015.
[18]Woodhouse I H. 微波遙感導論[M]. 董曉龍, 徐星歐,徐曦煜, 譯. 北京: 科學出版社, 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.