林业科学  2017, Vol. 53 Issue (12): 62-72   PDF    
DOI: 10.11707/j.1001-7488.20171207
0

文章信息

王新闯, 吴金汝, 陆凤连, 焦海明, 张合兵
Wang Xinchuang, Wu Jinru, Lu Fenglian, Jiao Haiming, Zhang Hebing
基于ICESat-GLAS数据的波形去噪算法优选与森林冠顶高反演模型建立
Optimal Selection of Algorisms for Denoising ICESat-GLAS Waveform Data and Development of a Forest Crown Height Estimation Model
林业科学, 2017, 53(12): 62-72.
Scientia Silvae Sinicae, 2017, 53(12): 62-72.
DOI: 10.11707/j.1001-7488.20171207

文章历史

收稿日期:2017-03-27
修回日期:2017-06-30

作者相关文章

王新闯
吴金汝
陆凤连
焦海明
张合兵

基于ICESat-GLAS数据的波形去噪算法优选与森林冠顶高反演模型建立
王新闯, 吴金汝, 陆凤连, 焦海明, 张合兵    
河南理工大学测绘与国土信息工程学院 焦作 454000
摘要:【目的】比较基于不同窗函数的GLAS数据去噪算法和多种森林冠顶高反演模型的精度,优选波形去噪算法并确立对森林冠顶高估算精度较高的反演模型,为森林生物量估测等研究奠定数据基础。【方法】首先基于布莱克曼窗函数和高斯窗函数对GLAS数据进行去噪处理,采用RMSE和SNR定量比较2种波形去噪方法的去噪效果;然后对去噪效果最好的窗函数去噪后的波形提取波形参数,并分针叶林、阔叶林、针阔混交林和不分林型4种情况,采用线性回归方法,以波形长度为参数建立波形参数模型,以波形长度、地形指数为参数建立地形因子模型,在地形因子模型基础上,逐步引入波形半能量高、波形前缘长度和波形后缘长度等参数建立全模型,比较3类模型的模拟效果。【结果】高斯窗函数去噪后的RMSE较低、SNR较高,去噪效果较优;冠顶高反演模型中,分林型和不分林型情况下,全模型模拟效果均优于其他2类模型。其中地形因子模型中针叶林效果较好:R2=0.853,RMSE=2.519 7 m;全模型中针阔混交林效果最好:R2=0.972,RMSE=1.001 4 m。【结论】高斯窗函数对GLAS波形去噪能力较强,且在复杂地形情况下,当引入多种波形参数结合地形因子建立多元线性回归模型时,模型对各林型最大冠顶高的解释能力显著提高,可在一定程度上克服地形因子模型在坡度较大地区对冠顶高解释困难的问题,实现复杂地形情况下森林冠顶高的精确估算。
关键词:大光斑激光雷达    布莱克曼窗去噪    高斯窗去噪    森林冠顶高反演模型    
Optimal Selection of Algorisms for Denoising ICESat-GLAS Waveform Data and Development of a Forest Crown Height Estimation Model
Wang Xinchuang, Wu Jinru, Lu Fenglian, Jiao Haiming, Zhang Hebing    
School of Surveying and Land Information Engineering, Henan Polytechnic University Jiaozuo 454000
Abstract: 【Objective】By comparing the GLAS data denoising algorithm based on different window functions and the accuracy of different inversion models of forest crowns, the waveform denoising algorithm is optimized and the inversion model with high estimation accuracy is obtained, which can lay the foundation for the study of forest biomass estimation.【Method】Firstly, the GLAS data are denoised by the Blackman window function and the Gaussian window function, and the effects of two kinds of waveform denoising method is quantitatively compared with the root mean square error(RMSE)and the signal to noise ratio(SNR). Then, the waveform data denoised by the window function that produced a better denoising effect were used to extract waveform data parameters. Four forest types (coniferous forests, broadleaf forests, mixed coniferous and broadleaf forests, and all forests)were included in this study. The linear regression method was used to develop a waveform-parameter model with the waveform length and a topographic-index model with the waveform length and the topographic-index. A full-parameter model was based on the topographic-index model with the height of median energy, the waveform leading edge, and the waveform trailing edge. Finally, the result of the crown height estimations of these three models were compared.【Result】The result show that the Gaussian window function gave a lower RMSE and a higher SNR than the Blackman window function. This indicates that use of the Gaussian window function for denoising provided more accurate data. For all forest types, the predictions of crown height given by the full-parameter model were more accurate than those given by other two models. The predictions of crown height for the coniferous forest given by the topographic-index model(R2=0.853, RMSE=2.519 7 m) were more accurate than the predictions of that model for other forest types. The predictions given for the mixed coniferous and broadleaf forest by the full-parameter model(R2=0.972, RMSE=1.001 4 m)were more accurate than the predictions of that model for other forest types.【Conclusion】In conclusion, the Gaussian window function performed better than the Blackman window function in denoising the GLAS waveform data. For the waveforms reflected from complex terrain, when the new multiple linear regression model was developed that included several waveform data parameters and the topographic-index, the performance of the model in interpreting the maximum crown height was significantly improved. This makes our model overcome the difficulty of interpreting the maximum crown height of forest on the highly sloping terrain, and thus an accurate estimation of forest crown height on complex terrains can be achieved.
Key words: large footprint laser radar    Blackman window denoising    Gaussian window denoising    forest crown height inversion model    

森林冠顶高是评价森林生态系统功能的核心参数之一,全球碳循环、森林生产力、气候变化等研究均需要量化森林冠顶高。传统的森林冠顶高参数测定方法需要耗费大量的人力、物力和时间,不利于大范围地测定森林参数,而且在实际应用中还存在测量仪器误差、样地空间分布代表性不足等问题(娄雪婷等,2011)。近年来,被动光学遥感技术在森林参数估测领域应用广泛,但多用于获取森林水平结构信息,很少能准确反映森林垂直结构信息,限制了遥感在估测森林冠顶高等方面的精度(邢艳秋等,2008李立存等,2011孙国清等,2006)。激光雷达技术具有很强的穿透能力,可以高精度地获取森林垂直结构信息,搭载于冰、云和陆地高程卫星(the ice, cloud and land elevation satellite,ICESat)上的地学激光测高系统(the geoscience laser altimetry system, GLAS)星载激光雷达ICESat-GLAS数据在森林垂直结构参数中得到了广泛应用。

大光斑激光雷达信号是入射激光脉冲与森林、地面相互作用的结果,激光回波信号是由发射脉冲和激光光斑内的森林、地表参数共同决定的,前者包括工作的波长、脉冲宽度、脉冲能量、光斑尺寸和记录回波脉冲的时间间隔,后者包括激光光斑内每株树的位置、高度、树种、树冠大小和形状、冠层的反射率以及地表的反射率、坡度、坡向等,这些因素都会对基于GLAS数据的森林冠顶高反演精度产生影响(庞勇等,2006)。此外,以森林分布为主要特征的激光脚点处的陡坡、枯落木等地物信息造成波形数据失真、叠加,也增加了森林冠顶高的反演难度(刘经南等,2005)。因此,在基于GLAS数据提取森林冠顶高的研究中,首先要对原始波形进行去噪处理。近年来,国内外学者在GLAS回波波形去噪方面做了很多努力。Iqbal等(2013)基于傅里叶变换对GLAS波形进行去噪并建立森林冠顶高估算模型,该方法能解释研究区79%的森林冠顶高变化。邱赛等(2012)在基于小波变换的基础上,选择常用的小波基函数对GLAS波形进行去噪,结果发现sym7小波基比db1小波基去噪效果更好,均方根误差为0.912 325,证明sym基函数更适合GLAS波形去噪。王爱娟等(2016)以吉林汪清区为例,提出了基于窗函数的林区GLAS数据去噪方法,结果表明布莱克曼窗函数对林区波形数据去噪效果最好,预测冠顶高度与实测冠顶高度的回归精度由0.725增至0.820。以上研究表明,这些去噪算法在一定程度上降低了背景噪声对GLAS回波波形的影响,但在地形复杂的地区适应性不强,因此有必要对不同去噪算法进行研究,选取去噪效果更好的算法。目前关于高斯窗函数对波形去噪的研究较少,由于GLAS激光发射器发射的信号为高斯脉冲信号,因此本研究在前人研究的基础上,对高斯窗函数的波形去噪效果进行研究,以进一步提高森林冠顶高反演精度。

由于GLAS按照较大尺度离散采样,地面的坡度和粗糙度混淆了来自植被冠层的返回信号,且地面的坡度和粗糙度会使脉冲宽度变宽,从而影响GLAS估算植被冠顶高度的精度(Pang et al., 2006)。邢艳秋等(2009)研究显示,在0°~15°地形坡度范围内,GLAS完整波形数据能够合理地反演森林冠顶高度,解释能力在51%~90%之间,当坡度超过15°时,地形对GLAS波形的影响很大,冠顶高估算模型对变量的解释能力显著下降。为了减弱地形的影响,Lefsky等(2005)开发了“地形指数法”,以波形长度和地形指数为自变量,提高了坡度较大地区森林最大冠顶高的提取精度,最终所建估测模型能解释59%~68%的森林最大冠顶高变化,并进一步利用森林冠顶高估测了森林生物量,所建模型R2达到0.73。Xing等(2010)在长白山地区应用全波形GLAS数据,以脉冲宽度的对数和地形指数为独立变量,对Lefsky等(2005)所建模型进行改进,改善了倾斜地面最大冠顶高估算精度,此模型能解释在0°~30°坡面上56%~92%的最大冠顶高变化。Hayashi等(2013)分析了不同因素对冠顶高估算的影响,认为GLAS数据的获取时间和信噪比以及坡度均会对森林冠顶高的估算结果产生显著影响。在剔除因获取时间和信噪比低导致估算误差较大的数据后,采用地形指数法分坡度区间(小于10°和大于等于10°)进行建模,提高了最大冠顶高的估算精度。董立新(2008)改进了“地形指数法”,提出“质心-地形指数法”,对坡度大于20°的地区应用“质心-地形指数法”,提高了坡度较大地区的森林最大冠顶高提取精度。汤旭光(2013)在地形指数法的基础上,通过引入波形半能量高这一参数,提高了最大冠顶高的估算精度。在上述研究中,基于GLAS估算森林冠顶高均需要区域的DEM支持,为了消除对DEM的依赖,缓解地形的影响,从而得到更加准确的冠顶高,Lefsky等(2007)在7个研究区基础上利用3个波形指数(波形宽度、下端边缘纠正系数和上端边缘纠正系数)修正了光斑内森林冠顶高的变化和地形的影响,定量估算得到了平均冠顶高,并且估算精度较以往研究有所提高。前人研究表明,森林冠顶高模型参数的选取对森林冠顶高反演精度具有重要影响。尽管基于地形指数和仅利用波形参数建立冠顶高反演模型都在一定程度上减弱了坡度对冠顶高估算精度的影响,但目前结合地形指数和多种波形参数建立冠顶高模型估算冠顶高的研究较少。鉴于此,本研究尝试基于地形数据和GLAS全波形数据提取的多种波形参数,对各林型分别建立波形参数模型、地形因子模型和全模型,检验不同模型对森林冠顶高的解释能力,确立反演精度较高的冠顶高反演模型,提高冠顶高估算精度,为森林生物量估测等研究奠定数据基础。

1 研究区概况和数据 1.1 研究区概况

研究区设置在露水河林区及其附近区域,该林区(127°29′—128°02′E, 42°24′—42°49′N)位于吉林省抚松县境内,长白山西北麓。全局经营区东西宽约40 km,南北长约50 km,总经营面积121 295 hm2。东南部地势起伏不大,比较平坦,西北部地势起伏较大,南部位于起伏的熔岩高台地,坡度为0°~53°。林区属中温带大陆性气候,气温较低,降水充沛,年平均气温0.9~1.5 ℃。研究区森林属长白山顶级植物群落区系,主要为阔叶林和红松(Pinus koraiensis)阔叶混交林,大部分为天然林和天然次生林,其中成熟林居多,森林高大茂密,郁闭度较高,也分布有部分人工林。该区森林类型有阔叶林、针阔混交林和针叶林,三者面积比例大致为3:1:1。本研究样地在不同坡度区间(0°~5°、5°~10°、10°~15°、15°~20°、20°~25°、25°~30°、> 30°)均有分布,样地林型及不同坡度区间的分布情况如表 1所示。

表 1 各林型样地不同坡度区间分布 Tab.1 Distribution of different slope intervals for each forest type

表 1可知,各林型样地在不同坡度区间均匀分布,阔叶林在研究区分布范围最广,地形最为复杂。总体来说,所有实测数据在平缓地区和坡度较大地区均匀分布,在10°~25°坡度范围内分布最多,样地地形复杂。

1.2 ICESat-GLAS数据

星载激光雷达ICESat-GLAS数据来自搭载于美国2003年发射的科学实验卫星ICESat上的第1台星载激光雷达传感器GLAS。星载激光雷达ICESat-GLAS数据免费获取,可用来测量冰盖高及其随时间的变化、云层和气溶胶的高度、植被的高度以及海冰的厚度等,是第1个用于连续全球观测的星载激光测高系统(邢艳秋等,2009)。GLAS搭载3个激光仪器,以91天重复轨道方案来执行观测任务,相邻轨道间距在赤道附近为15 km,在纬度80°处间距为2.5 km,激光光斑直径为60~70 m,光斑间隔为170 m。GLAS通过测量脉冲从卫星到地面的往返时间来测量地面高度,标准的脉冲宽度为4 ns,相当于地面高度0.6 m,解压后的1 000帧数据采样间隔为1 ns,对应地面高度为0.15 m,可以通过数据帧数来计算地面高度。本文使用研究区及周边附近地区的2006年的L3F、L3G和2007年的L3I的GLA01数据进行研究,数据从美国国家冰雪数据中心(http://nsidc.org/data/ice-sat/)下载。

1.3 样地实测数据及其处理 1.3.1 样地设置与调查

利用已经搜集到的林区GLAS数据、数字林相图、DEM和森林资源二类调查资料进行样地设置。露水河林区的有效GLAS光斑数据主要分布在研究区西北部,所有林型和坡度区间均有光斑分布。考虑林型、坡度对估算冠顶高的影响和样地的代表性,所设样地涵盖了3种不同林型(阔叶林、针叶林和针阔混交林)和不同坡度区间(0°~5°、5°~10°、10°~15°、15°~20°、20°~25°、25°~30°、> 30°)的森林。由于GLAS光斑为60 m×70 m的椭圆形状,因此在实际样地测量中,为与GLAS光斑相对应,以光斑中心为圆心,利用全站仪设置半径为30 m的圆形样地。在样地内,以样地最大冠顶高树木为中心设置半径为7.5 m的圆形样方,并在样地不同冠顶高区间(以10 m为间距)设置半径为7.5 m的圆形样方。为避免样方光谱信息互相干扰,样方之间距离尽可能大。本文所选用的47个样地数据于2015年8—9月实测得到,样地空间分布如图 1所示。

图 1 野外样地分布 Figure 1 Distribution of sample plots

样地调查内容主要包括GPS经纬度、高程、森林类型、坡度、坡向、样地情况、树种、胸径、树高、郁闭度、叶面积指数和林龄等。部分森林参数详细测量方法如下:用高精度GPS定位样地和样方中心坐标及其高程。对样地和样地内的样方中胸径大于2 cm的乔木进行每木检尺,利用激光测高仪获得树高,通过胸径尺测量距地面1.3 m处的立木直径。测量胸径时,胸径尺与树干垂直且紧贴树干周围。本文利用样地内与最高树木相差2 m的所有树高的平均值作为该样地树木实测最大冠顶高。

1.3.2 样地数据处理

由于样地调查时间和获取的GLAS数据时间不一致,为了减弱该期间因林木自然生长、森林采伐或自然灾害等因素导致的森林变化,首先与林业局管理人员进行沟通,查阅其2006—2014年有GLAS数据分布区域的采伐记录,对比与GLAS数据时相一致的SPOT5影像并实地调查期间的SPOT5影像,剔除那些位于采伐区或受自然灾害影响区域的GLAS数据;对于因林木自然生长引起的冠顶高变化,根据当地林龄与树高和胸径的关系予以校正,获取与GLAS数据时间一致的森林胸径、树高等基本数据。然后根据冠顶高数据,获取各样地和样方内的最大冠顶高。

2 数据处理与研究方法 2.1 技术路线

整个研究的开展主要包括遥感数据获取与处理、样地数据获取与处理以及森林冠顶高模型建立3大部分。首先对GLAS波形进行预处理,将GLA01产品中记录的激光雷达回波原始数据提取出来,由于激光脉冲受大气环境和激光器本身性能的影响,其接收能量会发生变化,且地形也会改变波形能量,因此需要将波形数据进行标准化,即用每刻的接收能量除以回波总能量,标准化波形对应的面积都为1,便于对比分析不同脉冲的波形;其次采用布莱克曼窗和高斯窗2种窗函数对GLAS波形进行去噪处理;然后利用SNR和RMSE 2种精度衡量指标,对2种波形去噪方法进行定量比较,优选出去噪效果最优算法;最后对去噪后的波形进行波形参数提取,并根据DEM和提取后的波形参数建立冠顶高反演模型。具体方案及路线如图 2所示。

图 2 技术路线 Figure 2 Technical roadmap
2.2 GLAS波形去噪算法 2.2.1 高斯窗函数

高斯窗函数是一种指数窗,主瓣较宽,因此其频率分辨率低,无负的旁瓣。第一旁瓣衰减达-55 dB,常用来截短一些非周期信号,如指数衰减信号。

假设f(t)∈L2(R),以g(t)作为窗函数的窗口傅里叶变换定义为:

$ {\rm{W}}{{\rm{F}}_{\rm{R}}}\left({\omega, b} \right) = \int_{ - \infty }^{ + \infty } {f\left(t \right)g\left({t - b} \right){{\rm{e}}^{ - j\omega t}}{\rm{d}}t} 。$ (1)

式中:b为时间位移;ω为频率位移;f(t)为在离散时间t的信号;WFR(ω, b)为t=b附近的频谱特性;j为虚数单位。

令:

$ {g_{\omega, b}}\left(t \right) = g\left({t - b} \right){{\rm{e}}^{ - j\omega t}}, $ (2)

则:

$ {\rm{W}}{{\rm{F}}_{\rm{R}}}(\omega, b) = \int_{ - \infty }^{ + \infty } {f\left(t \right){{\overline g }_{\omega, b}}\left(t \right){\rm{d}}t = \left\langle {f\left(t \right), {g_{\omega, b}}\left(t \right)} \right\rangle } 。$ (3)

g(t)的有效窗口宽度为Dt,则WFR(ω, b)给出的是f(t)在[b-Dt/2, b+Dt/2]时间区间内的频谱信息。Dt越小,对信号的定位能力越强。

假设f(t)的傅里叶变换为F(η), gω, b(t)的傅里叶变换为Gω, b(η), 根据Parseval定理:

$ {\rm{W}}{{\rm{F}}_{\rm{R}}}(\omega, b) = \left\langle {F\left(\eta \right), {G_{\omega, b}}\left(\eta \right)} \right\rangle /2{\rm{ \mathit{ π} }}。$ (4)

G(η)的有效窗口宽度为Dω,则WFR(ω, b)给出的是f(η)在[ω-Dω/2, ω+Dω/2]时间区间内的频谱信息。Dω越小,对信号的定位能力越强。

从公式可以看出,根据不确定原理,傅里叶变换的时域窗口和频域窗口宽度不可能同时提高。在所有窗函数中,只有高斯函数DtDω的乘积最小,因此为了尽可能减小信号通过滤波器失真的同时保证滤波器有效地进行滤波,选择高斯窗函数对信号进行去噪是不错的选择。

2.2.2 布莱克曼窗函数

布莱克曼窗的时域形式可以表示为:

$ \begin{array}{l} w\left(t \right) = 0.42 - 0.5{\rm{cos}}\left({2{\rm{ \mathit{ π} }}\frac{{t - 1}}{{N - 1}} } \right){\rm{ }}+\\ \;\;\;0.08\left({4{\rm{ \mathit{ π} }}\frac{{t - 1}}{{N - 1}}} \right), t = 1, 2, \cdots, N。\end{array} $ (5)

频域特性为:

$ \begin{array}{*{20}{c}} {w\left(\omega \right) = 0.42{W_{\rm{R}}} + 0.25[{W_{\rm{R}}}\left({\omega - \frac{{2{\rm{ \mathit{ π} }}}}{{N - 1}}} \right) + {\rm{ }}}\\ {{W_{\rm{R}}}\left({\omega + \frac{{2{\rm{ \mathit{ π} }}}}{{N - 1}}} \right)] + }\\ {0.04\left[ {{W_{\rm{R}}}\left({\omega - \frac{{4{\rm{ \mathit{ π} }}}}{{N - 1}}} \right) + {W_{\rm{R}}}\left({\omega + \frac{{4{\rm{ \mathit{ π} }}}}{{N - 1}}} \right)} \right]}。\end{array} $ (6)

式中:t为离散时间;ω为频率;WR(ω)为矩形窗函数的幅度频率特性函数;N为窗函数长度。

布莱克曼窗函数的最大旁瓣值比主瓣值低57 dB,但主瓣宽度是矩形窗函数主瓣宽度的3倍,为12 π/N,虽其频率识别精度低,但幅值识别精度高,有更好的选择性。

2.3 窗函数设计FIR滤波器去噪

本文基于Matlab软件,用2种窗函数设计FIR滤波器对波形进行去噪处理,具体过程如下。

1) 根据过渡带宽和阻带衰减要求,选择窗函数的类型并估计窗口长度N。本文选择高斯窗函数和布莱克曼窗函数。经多次对比试验,窗函数长度N取值为21。然后,调用Matlab中的窗函数求出中心值归一化为1的窗函数序列w(n), n为离散时间。

2) 根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)。根据线性叠加原理,对于其理想单位脉冲响应hd(n),可以根据给定的窗函数长度N由Matlab语句得到。

3) 计算滤波器的单位脉冲响应h(n),即滤波器的系数向量。h(n)是理想单位脉冲响应和窗函数的乘积, 在Matlab中用点乘命令表示为h(n)=hd(nw(n)。

4) 将原始波形数据分别与各窗函数相对应的滤波器的系数向量相互作用,进行波形去噪处理。

2.4 GLAS波形参数提取

本研究采用4组回波波形参数尝试建立和改进森林冠顶高反演模型:波形长度L、波形半能量高(HOME,height of median energy)H0、波形前缘长度l1和波形后缘长度l2。在噪声估计的基础上,利用噪声均值与其4倍的标准差之和作为波形信号的始末阈值。对于波峰位置的确定,首先计算信号的一阶导数,一阶导数由正变负的零点位置即为波峰,以此检测波形中的所有波峰。由于信号穿过冠层时,穿透能力减小,能量产生衰减,因此以信号最大强度的1/5作为检测波峰的阈值,从信号起始位置向后搜索,第1个超过波峰阈值的波峰视为波形的第1个波峰位置,波形前缘长度l1即为第1个波峰位置与信号起始位置的差值;从信号结束位置向前搜索,第1个超过波峰阈值的波峰即为地面波峰位置,波形长度L为地面波峰位置与信号起始位置的差值,波形后缘长度l2为地面波峰位置与信号结束位置的差值。波形半能量高是回波中植被部分能量一半的高度,波形的质心位置即波形面积一半处的位置,该位置与地面回波位置的差值即为波形半能量高H0。各波形参数如图 3所示。

图 3 GLAS波形参数 Figure 3 GLAS waveform parameters
2.5 GLAS森林冠顶高反演模型建立

在基于GLAS数据估算森林冠顶高的研究中,直接法和统计法较为常用(Iqbal et al., 2013Lefsky et al., 20052007Rosette et al., 2008Popescu et al., 2011Chen,2010)。直接法利用GLAS回波信号开始位置和地面回波位置的垂直高度差来估算冠顶高,统计法则基于统计分析来估算冠顶高。在地形较为平坦的区域,基于GLAS完整波形,利用直接法提取波形长度即可得到较为准确的森林冠顶高,如图 4a中,冠顶高即为信号起始位置和地面回波位置的差值;但当坡度增大时,回波信号受到影响,植被冠层部分会与地面回波信号混杂,且地面回波位置及最大冠顶位置会发生展宽,精确提取冠顶高较为困难。如图 4bf为GLAS光斑大小,θ为地面坡度,ftanθ为地面回波信号向上延伸的高度,即坡度对信号影响的垂直距离。地面回波位置随着坡度的增大其延伸也越大,仅仅依靠波形长度对冠顶高进行估算受地形影响较大,因此,采取一定方法减弱地形影响,对冠顶高的精确提取具有重大意义。

图 4 坡度对GLAS回波信号的影响 Figure 4 Effects of slope on GLAS echo signal

Chen(2010)研究表明,与直接法相比,统计法可以提高森林冠顶高估算精度,通过引入地形指数和GLAS波形参数可在一定程度上减弱地形的影响(Xing et al., 2010Lefsky et al., 2007Nie et al., 2015)。前人研究中普遍采用地形起伏度(高程最大和最小值之差)(Iqbal et al., 2013Lefsky et al., 2005Hayashi et al., 2013Chen,2010)和坡度对回波信号的垂直影响距离(汤旭光,2013)作为地形指数。由表 1可知,本研究区并非平坦区域,且各林型在不同坡度区间(0°~5°、5°~10°、10°~15°、15°~20°、20°~25°、25°~30°、>30°)均有分布,因此本研究引入地形指数,基于统计法建立森林冠顶高反演模型。

本研究地形指数的具体选取流程为:首先提取GLAS激光光斑所处地理位置多种采样窗口(3×3,5×5,7×7,9×9)下的地形起伏度并实测样地坡度对回波信号的垂直影响距离,然后计算波形长度与地面实测最大冠顶高的差异,通过建立地形起伏度和坡度对回波信号影响的垂直距离与高度差异之间的相关性,从而评估哪一个地形指数最适合估算因地形坡度引起的差异,最后选择其作为地形指数。本研究中地形起伏度从研究区10 m分辨率的DEM影像中提取,坡度对回波信号的垂直影响距离根据样地实测坡度计算得到。2种地形指数与高度差异之间的相关性分析如表 2所示。

表 2 高度差异与各地形指数的Pearson相关系数 Tab.2 Pearson correlation coefficients between forest crown height difference and terrain relief of different moving windows

表 2可知,5×5采样窗口下地形起伏度与高度差异之间的相关性最大,因此本文以5×5采样窗口下地形起伏度作为地形指数g参与森林冠顶高模型建立。

森林冠顶高模型参数选取对森林冠顶高反演精度有重要影响。相比于前人仅利用波形长度结合地形因子或仅利用波形参数建立的模型,本研究尝试建立3类模型:采用线性回归方法,首先建立为以波形长度为参数的波形参数模型及以波形长度和地形指数为参数的地形因子模型,然后在地形因子模型基础上,逐步引入波形半能量高、波形前缘长度、波形后缘长度等参数分林型建立全模型。引入波形半能量高,是利用其对冠层的垂直分布敏感而对坡度不敏感的特性作为冠顶高反演的重要参数(董立新,2008);波形前缘长度可以反映植被冠层和地形复杂起伏对回波信号的综合影响,引入该参数可有效抵消由于植被自身特征对回波造成的影响(Lefsky et al., 2007);波形后缘长度可以反映地形坡度和粗糙度对回波信号的影响(Lefsky et al., 2007),引入该参数可在一定程度上修正由于DEM对地形表达不准确对冠顶高提取的影响。所建冠顶高反演模型及参数如表 3所示。

表 3 森林冠顶高反演模型及参数 Tab.3 Forest crown height inversion models and parameters

通过各林型实测数据的回归分析确定不同模型的参数。

3 结果与分析 3.1 波形去噪结果对比与分析

本文选择研究区内47个脉冲点分别基于高斯窗函数和布莱克曼窗函数对GLAS波形进行去噪处理,其中1号点经过去噪处理后的波形如图 5所示。

图 5 去噪前后波形对比 Figure 5 Comparison of waveform before and after denoising

图 5可知,1号点经高斯低通滤波和小波去噪后,波形均得到了一定程度的平滑,保留了原始波形的大部分细节信息。为定量检验2种方法的去噪效果,本文选用均方根误差(root mean square error, RMSE)和信噪比(signal to noise ratio, SNR)2个指标对去噪效果进行评价。计算公式如下:

$ {\rm{RMSE}} = \sqrt {\frac{{{{\sum\limits_{i = 1}^L {\left[ {s\left( i \right) - f\left( i \right)} \right]} }^2}}}{L}} ; $ (7)
$ {\rm{SNR}} = 10 \times {\rm{lg}}\left[ {\frac{{\sum\limits_{i = 1}^L {{f^2}\left(i \right)} }}{{\sum\limits_{i = 1}^L {{{\left[ {s\left(i \right) - f\left(i \right)} \right]}^2}} }}} \right]。$ (8)

式中:s(i)表示原始信号强度(V);f(i)表示去噪后信号强度(V);i表示信号位置;L表示信号长度。

根据式(7)和(8),对所选47个脉冲点的布莱克曼窗去噪和高斯窗去噪效果分别进行定量评价,结果表明:由高斯窗去噪得到的信号RMSE大部分都小于布莱克曼窗去噪后的RMSE,SNR大部分都大于布莱克曼窗去噪后的SNR。经布莱克曼窗去噪处理后波形的信噪比平均值为20.958 81,均方根误差0.000 285;经高斯窗去噪处理后波形的信噪比平均值为23.360 704,均方根误差0.000 233。信噪比越高,均方根误差越低,说明信噪分离效果越好。由此可得,高斯窗去噪方法信噪分离效果优于布莱克曼窗去噪方法,能够较好地保留原始信号的有用信息。

3.2 冠顶高模型精度评价与分析

本研究共利用47组实测冠顶高数据及地形指数和去噪前后的不同波形参数进行回归分析,回归分析在SPSS软件中进行。各林型在不同模型下的去噪前后模型精度如表 4所示。

表 4 冠顶高反演模型精度比较 Tab.4 Accuracy comparison of forest crown height inversion models

采用R2、调整后的R2和RMSE 3个指标对波形去噪前后所建模型对各林型冠顶高拟合效果进行评价。由表 4可知,去噪前所建模型的R2范围为0.533~0.848,RMSE范围为2.475 0~3.836 4 m;而去噪后各林型所建模型的R2显著提高,范围为0.686~0.972,RMSE普遍降低,范围为1.004 1~2.756 2 m。这说明,对波形进行去噪,可以提取较为精确的波形参数,进而可以提高冠顶高估算精度。

表 4可知,去噪后5种模型均能较好地对复杂地形条件下的森林冠顶高进行估算。与波形参数模型相比,引入地形指数所建地形因子模型中,RMSE有所降低。针叶林和针阔混交林模型模拟效果最佳,R2高达0.853,阔叶林相对较差,R2为0.697。

在引入波形半能量高、波形前缘长度和波形后缘长度等波形参数后,所建模型改善了各林型的地形因子模型。董立新(2008)研究发现,在坡度较大地区,引入波形半能量高对坡度的敏感性较低,引入该参数对森林冠顶高的提取结果比较稳定。由于本研究区各林型在各坡度区间均有分布,因此本文引入波形半能量高建立模型3,减弱坡度较大地区对冠顶高估算的影响。模型3在各林型中总体较为稳定,该模型在地形指数校正模型的基础上引入波形半能量高,R2范围为0.721~0.915,RMSE范围为2.031 0~2.554 9 m;模型4在引入波形前缘长度后,冠顶高估算精度在不同林型中总体最高,但调整后的R2普遍降低。Lefsky等(2007)研究指出,波形后缘长度与地形最为相关,波形前缘长度与冠顶高变化有关,适用于估计平均树高而不是最大树高,这可以解释本文在估算最大冠顶高时引入波形前缘长度模型调整后的R2普遍降低的问题。进一步引入波形后缘长度建立模型5,可以发现,模型的R2有了进一步提高,范围为0.728~0.972,RMSE范围为1.001 4~2.216 4 m;效果最好的仍为针叶林和针阔混交林,其中针阔混交林的R2提高到0.972,RMSE为1.001 4 m,调整后的R2为0.902。对于各森林类型,全模型的RMSE范围在1.001 4~2.554 9 m之间,地形因子模型的RMSE范围在2.337 3~2.689 8 m之间。其中地形因子模型中针叶林效果最好:R2=0.853,RMSE=2.519 7 m;全模型中针阔混交林效果最好:R2=0.972,RMSE=1.001 4 m。当不分森林类型时,全模型R2范围为0.761~0.790,RMSE范围为2.559 7~2.734 9 m,总体效果较好,这说明当缺乏森林类型数据时,利用统一的森林冠顶高,基于GLAS数据建立冠顶高反演模型也能准确地估算森林冠顶高。

由以上分析可知,基于去噪后的GLAS波形数据估算森林冠顶高时,全模型总体优于波形参数模型和地形因子模型,针叶林优于针阔混交林,阔叶林最差。

4 讨论

基于GLAS波形数据提取森林冠顶高时,对GLAS波形的去噪效果在一定程度上影响波形参数的准确提取,而建立冠顶高反演模型所选取的模型参数可直接影响到所建模型对冠顶高的模拟效果,因此选择合适的波形去噪算法和模型参数至关重要。

王爱娟等(2016)探讨了汉宁窗、布莱克曼窗、矩形窗、三角窗和凯撒窗5种窗函数对林区GLAS波形数据去噪的效果,得出布莱克曼窗去噪效果最好。本文考虑GLAS激光发射器发射的信号为高斯脉冲信号,可假设激光回波由高斯分量构成(邢艳秋等,2009),因此采用高斯窗函数对波形进行去噪处理,结果发现高斯窗去噪效果优于布莱克曼窗,且从基于高斯窗去噪后的波形对各林型建立冠顶高反演模型的模拟结果来看,所建模型R2最大为0.972,高于王爱娟等(2016)所建模型R2最大为0.91的研究结果。这说明相比其他窗函数,高斯窗函数在由高斯分量构成的GLAS激光回波去噪处理中更有优势。在基于去噪后的GLAS波形直接估算森林冠顶高时,Iqbal等(2013)通过对GLAS波形进行去噪处理,利用直接法,可以解释79%的冠顶高变化,RMSE为3.18 m,冠顶高估算效果较好。本文在基于去噪后的GLAS波形建立模型估算森林冠顶高时,去噪后的各模型精度显著提高。以上研究表明,在基于GLAS数据估算冠顶高的过程中,对波形去噪方法的优选很有必要。

在以地形指数和波形长度为因变量建立的冠顶高估算模型研究中,Lefsky等(2005)所建估算模型能解释59%~68%的森林最大冠顶高变化,RMSE范围为4.85~12.66 m;Rosette等(2008)所建冠顶高估算模型R2为0.89,RMSE=2.99 m;仅利用波形自身参数所建冠顶高估算模型中,Lefsky等(2007)在7个研究区基础上利用3个波形指数(波形宽度、下端边缘纠正系数和上端边缘纠正系数)所建模型能解释83%的森林冠顶高变化,RMSE为5 m。相比于汤旭光(2013)利用地形指数和波形长度以及波形半能量高分林型建立多元线性回归模型,各林型调整后R2最大为0.904,RMSE最小为2.021 m,本研究在此基础上还引入了波形前缘长度和波形后缘长度,各林型所建模型R2最大为0.972,RMSE最小为1.001 4 m,模型的模拟效果相比前人研究有所改善。各林型中,针叶林和针阔混交林冠顶高模拟效果最好,这与汤旭光(2013)的结论一致。尽管在缺少林型数据时,本研究所建模型对冠顶高的反演能力相对较差,R2最大为0.79,RMSE最小为2.559 7 m,但总体来说,也能获取相对准确的森林冠顶高数据。以上分析表明,分林型引入多种波形参数建立冠顶高反演模型在地形复杂地区具有很强的适用性。本研究利用地形指数和波形参数虽在一定程度上修正了坡度对冠顶高提取的影响,但研究中影响冠顶高反演模型对森林冠顶高模拟精度的因素还包括时间差异、树种差异、样地地表状况等不确定性因素,因此在今后的研究中,基于地形数据和波形参数联合光学遥感数据,修正影响森林冠顶高反演精度的多种因素,对进一步提高森林垂直结构参数的提取精度有重大意义。

由于GLAS光斑呈条带形分布,本研究中样地布设原则是与GLAS光斑的位置大小相匹配,来估算GLAS光斑内森林冠顶高并顾及样地的可达性,所以本研究样地也呈条带形分布,这在一定程度上会影响模型的推广,为此,本研究在样地设置时,已综合考虑了研究区主要林型及地形情况,样地尽可能涵盖研究区森林类型及不同地形条件。为了进一步增强冠顶高反演模型的适用性,今后将通过补充野外调查,获取更多的样地数据,进一步优化模型,并结合光学遥感数据,综合考虑森林冠顶高空间扩展的可能有效因子,确定光斑内冠顶高空间外推的最佳尺度和模型,进行森林冠顶高的空间扩展研究。

5 结论

本文以吉林省露水河林区为研究区,基于大光斑激光雷达数据,采用布莱克曼窗函数和高斯窗函数对波形进行去噪,并定量比较2种方法的去噪效果。分别基于去噪后的GLAS波形,引入地形指数和波形参数建立不同森林类型冠顶高反演模型,通过对比2类冠顶高反演模型比较复杂地形中各反演参数对不同林型的冠顶高反演能力。结果表明:

1) 在对GLAS数据进行去噪处理中,与布莱克曼窗函数去噪相比,高斯窗函数去噪得到的RMSE较低、SNR较高,信噪分离效果较好,适应性较强。基于GLAS波形参数建立冠顶高反演模型时,对波形去噪后,所建模型显著提高了模型对森林冠顶高的模拟效果。因此在基于GLAS数据提取冠顶高时,对波形进行去噪处理很有必要。

2) 在地形复杂地区,基于GLAS数据分林型对冠顶高的反演能力比不分林型的反演能力强,对针叶林和针阔混交林的反演效果较好,所建模型的R2最高为0.972,RMSE最小为1.001 4 m;阔叶林反演效果较差,这可能是由于样地内最大冠顶高林木为落叶阔叶树种,GLAS获取时间为10月底,由于林木凋零,导致冠顶回波信号较弱而难以识别,从而导致估测值偏低。

3) 在所建森林冠顶高反演模型中,各林型的冠顶高估算效果有所不同。其中地形因子模型中针叶林效果最好:R2=0.853,RMSE=2.519 7 m;全模型中针阔混交林效果最好:R2=0.972,RMSE=1.001 4 m。总体来说,各模型中,针叶林和针阔混交林冠顶高模拟效果最好。

4) GLAS数据反演森林冠顶高时,引入多种波形参数显著改善了由地形指数建立的校正模型,其中效果最好的反演参数为波形长度、地形指数、波形半能量高和波形后缘长度,说明结合波形自身参数建立模型在森林冠顶高反演方面具有很大的应用潜力。

参考文献(References)
董立新. 2008. 基于多源遥感数据的三峡库区森林冠层高度与生物量估算方法研究. 北京: 中国科学院研究生院博士学位论文.
(Dong L X. 2008.Study on canopy height and biomass estimation method of Three Gorges Reservoir Area based on multi-source remote sensing data. Beijing:PhD thesis of Graduate University of Chinese Academy of Sciences.[in Chinese])
李立存, 张淑芬, 邢艳秋. 2011. 全站仪和测高仪在树高测定上的比较分析[J]. 森林工程, 27(4): 38-41.
(Li L C, Zhang S F, Xing Y Q. 2011. Comparative analysis on tree height measured by total station and vertex Ⅳ altimeter[J]. Forest Engineering, 27(4): 38-41. [in Chinese])
刘经南, 张小红. 2005. 利用激光强度信息分类激光扫描测高数据[J]. 武汉大学学报:信息科学版, 30(3): 189-193.
(Liu J N, Zhang X H. 2005. Classification of laser scanning altimetry data using laser intensity[J]. Geomatics and Information Science of Wuhan University, 30(3): 189-193. [in Chinese])
娄雪婷, 曾源, 吴炳方. 2011. 森林地上生物量遥感估测研究进展[J]. 国土资源遥感, 23(1): 1-8. DOI:10.6046/gtzyyg.2011.01.01
Lou X T, Zeng Y, Wu B F. Advances in the estimation of above-ground biomass of forest using remote sensing[J]. Remote Sensing For Land & Resources, 23(1): 1-8.
庞勇, 孙国清, 李增元. 2006. 林木空间格局对大光斑激光雷达波形的影响模拟[J]. 遥感学报, 10(1): 97-103.
(Pang Y, Sun G Q, Li Z Y. 2006. Large footprint lidar waveform modelling of forest spatial patterns[J]. Journal of Remote Sensing, 10(1): 97-103. DOI:10.11834/jrs.20060115 [in Chinese])
邱赛, 邢艳秋, 李立存, 等. 2012. 基于小波变换的ICESat-GLAS波形处理[J]. 森林工程, 28(5): 33-35, 59.
(Qiu S, Xing Y Q, Li L C, et al. 2012. ICESat-GLAS data processing based on wavelet transform[J]. Forest Engineering, 28(5): 33-35, 59. [in Chinese])
孙国清, RansonK J, 张钟军. 2006. 利用激光雷达和多角度频谱成像仪数据估测森林垂直参数[J]. 遥感学报, 10(4): 523-530.
(Sun G Q, Ranson K J, Zhang Z J. 2006. Forest vertical parameters from lidar and multi-angle imaging spectrometer data[J]. Journal of Remote Sensing, 10(4): 523-530. [in Chinese])
汤旭光. 2013. 基于激光雷达与多光谱遥感数据的森林地上生物量反演研究. 北京: 中国科学院大学博士学位论文.
(Tang X G.2013.Inversion of forest ground biomass based on LiDAR and multispectral remote sensing data. Beijing:PhD thesis of University of Chinese Academy of Sciences.[in Chinese]) http://cdmd.cnki.com.cn/Article/CDMD-80062-1013025748.htm
王爱娟, 邢艳秋, 邱赛, 等. 2016. 基于窗函数的林区ICESAT-GLAS波形数据消噪研究[J]. 西北林学院学报, 31(1): 214-220.
(Wang A J, Xing Y Q, Qiu S, et al. 2016. Denoising of forest ICESat-GLAS waveform data based on window function[J]. Journal of Northwest Forestry University, 31(1): 214-220. [in Chinese])
邢艳秋, 王立海. 2008. 基于森林生物量相容性模型长白山天然林生物量估测[J]. 森林工程, 24(2): 1-4.
(Xing Y Q, Wang L H. 2008. Biomass estimation of natural forest in Changbai Mountains based on forest biomass compatibility model[J]. Forest Engineering, 24(2): 1-4. [in Chinese])
邢艳秋, 王立海. 2009. 基于ICESat-GLAS完整波形的坡地森林冠层高度反演研究——以吉林长白山林区为例[J]. 武汉大学学报:信息科学版, 34(6): 696-700.
(Xing Y Q, Wang L H. 2009. ICESat-GLAS full waveform-based study on forest canopy height retrieval in sloped area-a case study of forests in Changbai Mountains, Jilin[J]. Geomatics and Information Science of Wuhan University, 34(6): 696-700. [in Chinese])
Chen Q. 2010. Retrieving vegetation height of forests and woodlands over mountainous areas in the pacific coast region using satellite laser altimetry[J]. Remote Sensing of Environment, 114(7): 1610-1627. DOI:10.1016/j.rse.2010.02.016
Hayashi M, Saigusa N, Oguma H, et al. 2013. Forest canopy height estimation using ICESat/GLAS data and error factor analysis in Hokkaido, Japan[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 81(7): 12-18.
Iqbal I A, Dash J, Ullah S, et al. 2013. A novel approach to estimate canopy height using ICESat/GLAS data:a case study in the New Forest National Park, UK[J]. International Journal of Applied Earth Observation and Geoinformation, 23(1): 109-118.
Lefsky M A, Harding D J, Keller M, et al. 2005. Estimates of forest canopy height and aboveground biomass using ICESat[J]. Geophysical Research Letters, 32(22): 441.
Lefsky M A, Keller M, Pang Y. 2007. Revised method for forest canopy height estimation from geoscience laser altimeter system waveforms[J]. Journal of Applied Remote Sensing, 1(1): 1-18.
Nie S, Wang C, Zeng H, et al. 2015. A revised terrain correction method for forest canopy height estimation using ICESat/GLAS data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 108: 183-190. DOI:10.1016/j.isprsjprs.2015.07.008
Pang Y, Li Z Y, Lefsky M, et al. 2006. Model based terrain effect analyses on ICESAT GLAS waveforms[J]. IEEE International Conference on Geoscience and Remote Sensing Symposium: 3232-3235.
Popescu S C, Zhao K, Neuenschwander A, et al. 2011. Satellite lidar vs. small footprint airborne LiDAR:comparing the accuracy of aboveground biomass estimates and forest structure metrics at footprint level[J]. Remote Sensing of Environment, 115(11): 2786-2797. DOI:10.1016/j.rse.2011.01.026
Rosette J A B, North P R J, Suarez J C. 2008. Vegetation height estimates for a mixed temperate forest using satellite laser altimetry[J]. International Journal of Remote Sensing, 29(5): 1475-1493. DOI:10.1080/01431160701736380
Xing Y, de Gier A, Zhang J, et al. 2010. An improved method for estimating forest canopy height using ICESat/GLAS full waveform data over sloping terrain:a case study in Changbai Mountains, China[J]. International Journal of Applied Earth Observation and Geoinformation, 12(5): 385-392. DOI:10.1016/j.jag.2010.04.010