1 引言

森林是一种重要的自然资源，在全球碳循环中发挥着重要作用。树木生物量是森林生态系统功能及生物圈、大气和海洋之间碳通量作用的关键，而树木的高度是生物量估计以及碳储量的重要参数^[1]。星载激光测高仪系统是一种通过接收和处理激光发射脉冲经大气传输在目标表面反射后的回波波形数据，具备测距和反演目标表面特性功能的遥感系统^[2]。该测高系统可以主动获取全球地表及目标三维信息，弥补了光学影像高程精度低的缺点，在获取植被高度、地面生物量特征以及其他植被结构参数方面具有绝对的优势^[3]。

基于激光测高系统采集的回波信号估测植被冠层高度时一般采用直接估算法^[4-9]，根据植被冠层高度的物理定义，计算回波信号中冠层顶部到地面返回位置之间的高程差值来获取。这类方法是通过波形数据本身包含的信息提取相关波形参数来反演冠层高度的，因此适用于各种森林类型，在估算植被地区冠层高度时得到了广泛应用。其中，植被冠层顶部由信号开始位置定义，信号开始位置由背景噪声阈值决定。然而，由于星载激光测高系统大光斑的特性和地表环境复杂等因素，从回波中提取准确地面返回位置具有一定的困难。

目前还没有广泛适用的地面返回位置提取方法。Sun等人^[10]基于波形的高斯分解，将波形分解后的最后一个高斯回波波峰所在位置视为地面返回。但最后一个高斯波也可能由地面低矮植被或粗糙度等造成，并不能准确代表光斑内地面高程。Chen等人^[11]将分解后的最后两个高斯回波中波峰振幅更大的回波波峰所在的位置确认为地面返回位置，取得了更好的反演结果。但在地形环境复杂的区域，GLAS回波复杂的情况下，单个高斯回波波峰也无法准确代表地面，进而对植被冠层高度反演带来更多误差。针对以上问题，本文提出了一种新的地面返回位置提取算法，首先利用广义高斯分解确定地面对应的回波位置范围，然后在地面范围波形内检测所有波峰，并根据设置的中值阈值，计算有效候选峰的加权高程，实现复杂波形下地面返回位置的准确提取和植被冠层高度反演。

1 研究方法

2003年美国国家宇航局发射了星载激光测高卫星ICESat（ice, cloud, and land elevation satellite），它搭载的GLAS（geoscience laser altimeter system）系统是第一个实现全球连续观测的星载激光雷达测距系统，目前已被广泛应用于全球或区域尺度的地学研究中。该系统以40 Hz频率向地表发射近红外波段（1 064 nm）的激光光束，用于地面测高。激光光束在地表形成直径约为65 m近圆形的激光足印，同一轨道的激光足印间隔为175 m^[12]。GLAS系统提供了15种标准数据产品，其中GLA01数据产品记录了完整原始波形数据，该波形数据反映了对应光斑内的地物信息，可以用于植被冠层高度估测。文中在进行植被高度反演时，首先对GLAS回波数据进行滤波，减少噪声对回波的影响；然后根据广义高斯分解方法，提取地面位置对应的波形范围以及提取信号的起始位置；根据地面位置对应的波形范围，进一步确定地面高程，从而确定冠层高度，研究过程主要步骤，如图1所示。

图 1. 高度反演及精度验证流程图

Fig. 1. Flow chart of canopy height inversion and accuracy verification

下载图片 查看所有图片

1.1 GLAS波形数据及预处理

由于激光测高仪在地面上的光斑直径往往在几十米的量级，因此，回波信号通常为多模式复杂曲线，由混入探测器噪声、背景噪声、量子噪声和热噪声的若干个高斯分量组成，除背景噪声为非零均值高斯白噪声外，其余都被认为是零均值高斯白噪声。高斯滤波是一类根据高斯函数的性状确定权值的线性平滑滤波方法，对于一直服从正态分布的噪声非常有效；高斯滤波器的宽度越大，对噪声的抑制效果越好，但同时会减小回波信号峰值，展宽高斯波形以及改变波形重心。文中采用自适应的高斯滤波方法^[9]对回波信号进行滤波，通过原始信号和滤波信号确定背景噪声的均值和方差，并根据公式（1）确定信号的噪声阈值：

式中：Noise_th表示噪声阈值；Noise_mean表示噪声均值；Noise_std表示噪声标准偏差；N表示倍数。

1.2 基于GLAS波形的地面返回位置提取

要实现GLAS波形地面返回位置的精确提取，首先要确定含有地面信息的波形范围。Allouis等人^[13]的研究表明GLAS的回波波形可以用2个广义高斯函数来近似，其中一个代表植被信号部分；另一个代表地面信号部分。因此，本文设置了2个步骤来正确识别地面返回位置。

1） 潜在的地面返回位置应该在GLAS波形后端的地面范围波形中，波形的地面范围代表了光斑内因地面地形变化引起的高程变化，由代表地面部分的广义高斯函数来确定，如图2（a）蓝色波形。使用广义高斯函数作为波形分解的建模模型时，回波可表达为

式中：A表示振幅值；σ表示波宽；μ表示时间位置；α表示形状参数；A_canopy、σ_canopy、μ_canopy、α_canopy表示植被部分波形的高斯参数；A_ground、σ_ground、μ_ground、α_ground表示地面部分波形的高斯参数。

图 2. GLAS波形地面返回位置提取方法

Fig. 2. Ground return position extraction method based on GLAS waveform

下载图片 查看所有图片

2） 检测出地面范围波形内所有的候选峰，将其中最大峰值的一半处设为中值阈值，超过中值阈值的候选峰为有效候选峰。若有效候选峰只有一个，则该有效峰为地面回波，其最大峰值处确认为地面返回位置；若有效候选峰不只一个，计算所有有效候选峰的加权高程，其中权值为有效候选峰的峰值，并将其作为估测的地面返回位置，如图2（b）所示。这里中值阈值的设置是因为GLAS波形仍然包含一些振荡，包括非常低的峰值，显然不是指地形信号，通过设置中值阈值可以排除这些无效的候选峰，有助于提取更准确的地面返回位置。图2给出了提取地面返回位置高程的2个过程。

1.3 基于GLAS回波植被冠层高度计算

本文应用的GLAS大光斑植被冠层高度的反演模型是直接估算法，根据植被最大冠层高度的物理定义，通过计算从GLAS回波信号中冠层顶部到地面返回位置之间的高程差值即可得到。在GLAS波形中，植被冠层顶部由信号开始位置定义，而信号的开始位置由信号的噪声阈值的位置确定。在用公式（1）确定阈值时，在不同区域森林参数的研究中，N的取值有4.5、4、3或2。根据Enßle等人^[9]的研究，在平坦地形下，将噪声阈值中N取3时通常有最佳的效果，故本文中回波开始位置取N=3的噪声阈值位置。本文植被冠层高度计算公式如下：

式中：H表示GLAS光斑内植被最大冠层高度；Wf_start表示植被冠层顶部位置高程；μ_ground表示地面返回位置高程。

2 结果与讨论

2.1 研究区域与数据

本文研究区域选取的是美国佛蒙特州西北方向植被地区，该地处于中高纬度地区（44ºN～45ºN，72ºW～73ºW），植被属于温带阔叶林和混交林的生物群落。这里森林分布广泛，有针叶树和北方阔叶林，糖枫、黄桦、美洲山毛榉和白灰是北方阔叶林中常见的主要树木，其他物种包括东部铁杉和东部白松。

本研究采用了佛蒙特州西北地区GLAS第3个激光器的L3i、L3j、L3k期激光数据，激光脉冲发射和接收过程中数据质量易受云和系统噪声的影响，所以首先需要对整体数据进行质量筛选^[14-15]。通过GLAS数据产品中提供的质量指标来完成对整体数据的过滤，筛选条件如下：1） 波形信噪比要大于10；2） GLAS参考高程和SRTM高程差小于100 m，除去因低云污染而失真的GLAS数据；3） 波形信号为非饱和信号；4） 为避免地形坡度过大造成的波形展宽对植被高度反演的影响，本研究使用地形坡度低于10°的波形数据，坡度值由SRTM数据计算得到。

机载激光雷达对植被空间结构和地形的探测能力较强，在植被参数定量测量和反演方面具有明显的优势，且机载Lidar数据估测的树高与实际测量树高具有高度的一致性^[16]。为了验证本文估测植被冠层高度的准确性，使用佛蒙特州地区同时期的机载Lidar点云数据为验证数据。这些数据是由徕卡ALS-50激光雷达系统在2008年2月11日至2009年4月16日之间秋冬季节收集的，且当时地面没有积雪。采集得到的激光雷达点云数据高程精度约为15.0 cm，水平精度约为0.5 m，激光点云间隔约为1.4 m，获取了2 424 km²的机载激光雷达点云数据，激光雷达点云密度为 2.02 points/m²。

2.2 GLAS回波地面返回位置提取结果

本文基于在地面范围波形内检测所有波峰，并设置中值阈值，进而实现复杂波形下植被区域的地面返回位置的准确提取，将本文使用的方法表示为方法1。为了验证本文方法提取地面返回位置的可靠性，将本文所用的方法与常用的地面返回位置的提取方法进行比较。一种常用的地面返回位置的提取方法是高斯分解法^[11]，将高斯分解后高程最低的2个高斯回波中，波峰振幅更大的回波波峰所在的位置确认为地面返回位置，这里表示为方法2。图3显示了方法1和方法2提取的地面返回位置结果，其中方法2高斯子波1和高斯子波2表示高斯分解后高程最低的2个高斯波，波峰振幅更大的回波波峰位置表示方法2提取的地面返回位置。

图 3. 结果比较

Fig. 3. Results comparison

下载图片 查看所有图片

由图3中提取结果可看到，总体上方法1提取到的地面返回位置与机载Lidar地面点云的分布更为接近，且方法1提取地面返回位置的准确性更高。尽管有一些情况下方法1和方法2提取的地面高程都是接近机载Lidar地面高程的，但方法2确认地面返回位置时表现得更不稳定，有2点可能的原因：1） 高程最低的2个高斯子波作为候选地面波峰时，其中可能包括很小的震荡，显然不是地形信号的主要部分作为地面候选峰，漏掉了真正的有效候选峰；2） 仅通过将最后2个高斯子波中波峰振幅更大的回波波峰位置确认为地面返回位置时，也可能除去了本属于地面信号部分且包含地面信号能量较多的子波，忽略了部分有效的地面信号。本文提出的方法1在解决以上问题时有很大的优势，中值阈值的设置可以有效地滤除一些无效候选峰，并且考虑了地面信号部分所有的有效候选峰后提取的地面返回位置，更能准确代表GLAS光斑下的地面高程。

2.3 GLAS波形植被冠层高度反演结果

基于本文提出的地面返回位置提取方法，进行进一步的植被冠层高度反演表示为方法1，另一种基于高斯分解的地面返回位置，反演植被冠层高度为方法2。选取了佛蒙特州地区148个GLAS回波数据，以该区域同时期的机载Lidar获取的植被冠层高度为验证数据。利用方法1和方法2的反演结果如图4所示。

图 4. GLAS植被冠层高度反演结果

Fig. 4. GLAS vegetation canopy height inversion results

下载图片 查看所有图片

图4中，分别给出了148个GLAS光斑内2种方法反演的植被冠层高度与机载Lidar植被冠层高度的对比。图中结果表明2种方法得到的数据点大部分位于1:1线两侧附近，但方法1得到的数据点在1:1线两侧分布更为集中，且方法1下估算的冠层高度与机载Lidar冠层高度差值绝大部分在10 m以内，只有少数几个异常值大于10 m，而方法2的异常值明显更多。以机载Lidar结果为验证数据，求得方法1下反演的GLAS植被冠层高度RMSE为2.82 m，相关系数R为0.81。方法2下反演的GLAS植被冠层高度RMSE为3.75 m，相关系数R为0.67。由此可见，方法1的GLAS冠层高度反演结果跟机载Lidar测量结果的一致性更高。

另外本文又进行了不同区域的研究，选取了美国俄怀明州Teton国家森林公园地区的72个GLAS回波数据，对光斑内植被冠层高度进行了反演，见图5。图中给出了Teton森林GLAS植被冠层高度与机载Lidar结果对比，结果显示本文方法下GLAS冠层高度反演结果与机载Lidar测量结果的匹配程度很好，RMSE为2.92 m，相关系数R为0.85，也说明本文方法在不同植被地区的适用性很好。

图 5. Teton区域GLAS冠层高度与机载Lidar冠层高度对比

Fig. 5. Comparison of GLAS canopy height with airborne Lidar canopy height

下载图片 查看所有图片

3 结论

本文是基于GLAS数据植被冠层参数反演方法的研究，设计了一种新的地面返回位置提取算法，进而用来反演植被冠层高度，与现有的方法比较，本文方法的处理结果与验证数据的一致性更好，并具有良好的地域实用性。本文主要研究结论如下：1） GLAS波形地面返回位置的确认是影响植被冠层高度反演精度的关键问题，由于GLAS光斑范围面积大，地形坡度、地表粗糙度、低矮植被等环境因素对波形影响较大，地面返回位置的确定存在一定的不确定性。本文通过用2个广义高斯函数近似表达GLAS回波的植被信号部分和地面信号部分，得到回波中地面信号的范围，通过在地面范围波形内检测所有波峰，并设置中值阈值，计算所有有效候选峰的加权高程，实现了复杂波形下的地面返回位置的准确提取。中值阈值的设置可以有效地防止一些因地表条件复杂等因素造成的尾边较小的回波被认为是地面信号的情况，提高了GLAS光斑下地面返回位置判别的准确性。2） 本文植被冠层高度反演结果RMSE为2.82 m，相关系数为0.81，与校验数据机载Lidar冠层高度一致性很好，结果优于基于高斯分解的冠层高度反演结果。3） 目前已有研究表明在坡度小于10°的植被区域，根据物理模型的直接估算法反演植被冠层高度结果是可靠的，本文的研究结果也侧面验证了这一结论。未来可以针对本文提出的GLAS植被冠层高度反演方法在不同坡度地形上的适用性和反演精度开展进一步的研究。