蔡耀通 林 辉 孙 华 张 猛 龙江平
( 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.