森林冠层高度是指森林冠顶与地面之间的垂直距离，是度量森林结构变化的重要参数。目前，森林冠层高度已应用于森林生物量和森林碳储量研究中，因此，准确估测森林冠层高度变得越来越重要。

激光雷达(light direction and ranging，LiDAR)是一种新型的主动遥感技术，其以高精度、高分辨率和主动工作的优点，在估测森林结构参数尤其是森林垂直结构参数方面具有其他技术无可比拟的优势(Van Leeuwen et al.，2010； Chopping et al.，2009； Hyde et al.，2007； 庞勇等，2005； 马利群等，2011； 李丹等，2011)。搭载在冰、云和陆地高程卫星(the ice，cloud，and land elevation satellite，ICESat)上的地学激光测高系统(geoscience laser altimeter system，GLAS)是目前唯一的极地轨道星载大光斑激光雷达设备，具有全天时、全天候观测的能力，并且能够免费提供完整波形数据，在估测森林冠层高度方面独具优势(Iqbal，2010； Gwenzi et al.，2014； Schutz et al.，2005)。Lefsky等(2005)利用GLAS波形长度及地形指数建立了森林冠层高度估测模型，模型的R²在0.59~0.68之间。Xing等(2010)针对吉林省长白山汪清林区，对Lefsky等(2005)提出的模型进行了改进，建立了森林冠层高度对数模型，在0~5°坡度范围内，模型的R²为0.92，在0~10°坡度范围内，模型的R²为0.80，随着坡度范围增大模型的R²逐渐降低，当坡度范围为0~30°时，模型的R²为0.56。

虽然ICESat-GLAS波形数据在估测森林冠层高度方面独具优势，但由于其光斑呈离散条带状分布，若要实现区域森林冠层高度的高精度估测，则需联合其他连续成像的光学遥感数据(吴迪等，2014； 杨婷等，2014； Anderson et al.，2008； Yu et al.，2015)。现有研究多是将ICESat-GLAS波形数据与多光谱数据联合，但由于多光谱数据的光谱分辨率低，像元内光谱信息易发生混合，导致像元内光谱特征提取受到限制，进而影响森林冠层高度估测精度(Simard et al.，2011； 徐定，2013； 董立新等，2011)。相比于多光谱数据，高光谱遥感数据包含了更加丰富的空间、辐射和光谱三重信息，能同时获取数十至上百个波段的连续光谱信息，且具有较高的光谱分辨率，在森林制图、森林资源调查等方面具有应用前景(谭炳香，2006)。

因此，本研究将ICESat-GLAS波形数据与HJ-1A/HSI高光谱数据联合，借助HSI高光谱数据提供的连续高分辨率光谱信息，实现区域森林冠层高度的估测，降低由于GLAS光斑呈离散条带状分布无法覆盖整个研究区造成的估测误差。首先从ICESat-GLAS 波形数据中提取波形参数(波形长度W和地形坡度参数TS)，建模反演GLAS光斑内的森林冠层高度； 然后基于支持向量回归机(support vector regression，SVR)算法，利用GLAS估测的森林冠层高度与HJ-1A/HSI高光谱数据的光谱信息建立森林冠层高度SVR估测模型，得到区域森林冠层高度分布图。

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

选择吉林省汪清林业局经营区(43°05′—43°40′N，129°56′—131°04′E)为研究区(图 1)。该研究区内平均海拔806 m，地面高程360~1 477 m，坡度0~45°。区内植物种类繁多，植被结构复杂，呈带状分布在海拔500~1 100 m范围内。针叶树以红松(Pinus koraiensis)、云杉(Picea asperata)和臭冷杉(Abies nephrolepis)为主，阔叶树以紫椴(Tilia amurensis)、蒙古栎(Quercus mongolica)、枫桦(Betula costata)、色木槭(Acer mono)和白桦(Betula platyphylla)为主。

1.2 ICESat-GLAS波形数据

ICESat-GLAS 激光雷达是美国国家航天局戈达德宇航中心于2003年发射的，其上搭载了3个激光器，每个激光器每秒发射40个波长为1 064 nm的激光脉冲，用于地物探测。激光脉冲的脉冲宽度为4 ns，每个激光脉冲在目标地物上形成一个直径约为70 m的椭圆形光斑，相邻光斑中心点间距为170 m，相邻轨道间距在赤道附近为15 km，在纬度达到80°处为2.5 km。

ICESat-GLAS提供15种数据产品(GLA01，GLA02，GLA03，… GLA15)，为了尽可能与野外数据获取时间一致，减少数据时间差异造成的误差，本研究共获取了L3E(2006年3月)、L3F(2006年6月)、L3G(2006年11月)、L3H(2007年3月)和L3I(2007年10月)5期GLA01和GLA14数据产品，数据获取网址为http://nsidc.org/data/icesat/。GLA01产品中包含完整的波形数据，对于陆地而言，每个波形数据包含544帧数据，每一帧对应的真实距离为0.15 m。GLA14产品中包含地面高程及经纬度等信息。2种数据产品间使用头文件中的索引号和光斑号进行匹配。

1.3 HJ-1A/HSI高光谱数据

环境减灾小卫星星座A/B星上搭载的超光谱成像仪(hyper spectral imaging radiometer，HSI)，是世界上第1个基于空间调制干涉成像原理设计的民用高光谱成像仪(Zhao et al.，2010)。HSI共包含115个波段，光谱范围为0.45~0.95 μm，空间分辨率为100 m，平均光谱分辨率为4.32 nm，幅宽为50 km，重访周期为96 h，且具备±30°的侧视能力。本研究共获取4景HSI高光谱遥感影像，网址为http://www.cresda.com/n16/index.html。由于云层及多景拼接的影响，数据获取时间存在差异，其中2景获取时间为2010年6月5日，1景获取时间为2010年6月9日，另外1景获取时间为2010年10月1日。

1.4 野外实测数据

本研究于2006年9月和2007年9月分别进行了野外数据调查。沿着ICESat-GLAS轨道方向随机选取183个GLAS光斑作为野外数据调查样地，每个样地的水平投影面积为500 m²。利用手持超声波测高仪(Vertex IV)测量样方内每株树的树高，取最大值作为此样方的森林冠层高度。使用坡度测量仪(HG4DQY-1)测量每个样地的地形坡度。

2 数据处理与研究方法 2.1 ICESat-GLAS波形数据处理及参数提取

利用IDL_Reader程序包先将二进制格式的GLA01和GLA14数据产品转换成ASCII格式，为了便于后续分析和计算，将其转换为由电压(单位： V)表示的格式。由于数据获取时间不同，为了统一分析计算，对数据进行标准化，即将每一时刻的接收能量除以回波总能量。对标准化后的数据进行滤波、高斯分解拟合等平滑处理，利用平滑后的波形数据提取波形参数。具体处理内容如下：

1) 波形平滑处理。ICESat-GLAS发射的激光脉冲是高斯脉冲，因此本研究使用高斯平滑滤波器对波形进行平滑滤波(Xing et al.，2010)，平滑后的波形如图 2所示。

2) 提取波形长度W。波形长度W定义为波形起波点与止波点之间的垂直距离。起波点与止波点由背景噪声阈值确定，背景噪声阈值的设定参考以往研究，采用波形前100帧数据的平均值与其4倍的标准差之和(Lefsky et al.，2005)。波形长度W计算公式如下：

$ W=(SigEnd-SigStart)\times BinSize。 $

(1)

式中：SigStart为波形起波点； SigEnd为波形止波点； BinSize为每帧数据代表的真实距离，对陆地而言，每帧代表真实距离0.15 m，因此BinSize取值为0.15 m。

3)提取地形坡度参数TS。一般在平坦地区，波形长度可近似为森林冠层高度；但在林地中，受地形坡度和光斑直径的影响，波形展宽现象明显(Duncanson et al.，2010； 庞勇等，2005； Yang et al.，2011)。目前，多是利用从DEM数据中提取的地形指数来校正地形坡度对波形的影响，但是DEM数据分辨率的高低也会造成一些不必要的误差。因此本研究尝试从ICESat-GLAS波形数据中提取地形坡度参数TS来校正地形坡度对波形的影响。

本研究中地形坡度参数TS定义为地面回波对应高斯分量的波形长度W_gf与光斑直径D比值的反正切值(Mahoney et al.，2014)，由于大气衰减、信号饱和等影响导致地面回波波形展宽W_m，因此用W_m对W_gf进行修正，最终得到TS的计算公式如下：

$ \text{TS}=\arctan \left( \frac{{{W}_{\text{gf}}}-{{W}_{\text{m}}}}{D} \right)。 $

(2)

W_m的计算采用Mahoney等(2014)对全球大量数据进行统计分析得到的公式：

$ {{W}_{\text{m}}}=4.689+0.759A。 $

(3)

式中： A为地面回波峰值能量值。

2.2 HJ-1A/HSI高光谱数据处理

HSI数据是二级产品，已经经过了几何校正和辐射定标，本研究在此基础上首先利用地形图对数据进行精校正，然后利用ENVI软件中的FLAASH大气校正模块对HSI高光谱数据进行大气校正。HSI高光谱数据由115个波段组成，相邻波段间相关性强、数据冗余度高，因此本研究利用最小噪声分离法(MNF)对数据进行降噪降维处理。MNF是在主成分分析基础上发展起来的一种以信噪比为衡量标准的线性变换算法，可以将信息集中到尽可能少的几个分量中，同时将噪声分离，被广泛应用于高光谱数据处理中(李海涛等，2007)。对HSI高光谱数据进行MNF处理后，发现前3个MNF(MNF1，MNF2，MNF3)分量包含的信息占信息总量的98%，能够满足研究要求，因此，本研究采用前3个MNF分量进行分析。为了匹配数据，将GLAS光斑中心点坐标加入到MNF分量图中，以光斑中心点为圆心建立直径为70 m的缓冲区，将缓冲区所覆盖像元的平均值作为光斑的MNF值。

2.3 GLAS森林冠层高度估测模型建立

本研究利用GLAS波形参数建立GLAS森林冠层高度估测模型。众多研究表明，地形坡度会造成GLAS波形展宽，是影响森林冠层高度估测精度的重要因素(Pang et al.，2011； Lee et al.，2011)。如图 3所示，由于地形坡度及光斑大小的作用，波形长度W不完全等于真实的森林冠层高度H，存在一定的误差，当地形坡度为α时，波形的最大展宽即最大估测误差为Dtanα。由此可见，地形坡度越大，GLAS波形展宽越明显，波形长度直接用于估测森林冠层高度的误差越大。

因此，本研究从ICESat-GLAS波形数据中提取地形坡度参数TS，尽量降低地形坡度造成的影响。从图 4可以看出，地形坡度参数TS与地形坡度实测值间存在着很强的线性关系(R²=0.78)，从ICESat-GLAS波形数据中提取地形坡度参数TS校正地形坡度对波形造成的影响是可行的。基于以上分析，拟建立模型如下:

$ H=\alpha W+bD\tan \text{TS}+c。 $

(4)

式中： a，b，c为相应系数。

2.4 基于SVR算法的森林冠层高度分布图

支持向量机是在统计学理论上发展起来的一种针对有限样本实现结构风险最小化思想的机器学习方法，具有较高的泛化能力。支持向量回归机(SVR)是支持向量机在函数回归领域的应用，能解决神经网络方法中无法避免的局部最优问题(曾绍华，2006)。该算法的基本思想是：将数据集X[(x₁，y₁)，(x₂，y₂)，···(x_n，y_n)]通过核函数K(x，x_i)转换到高维特征空间，然后在高维特征空间中进行线性回归，从而取得在原空间非线性回归的效果，并基于结构风险最小化原则，引入lagrange函数将回归问题转换为求最优解问题，如下所示：

$ \begin{matrix} \underset{{{\alpha }^{*}}}{\mathop{\max }}\,-\frac{1}{2}\sum\limits_{i=1,j=1}^{n}{\left( {{\alpha }_{i}}-\alpha _{i}^{*} \right)\left( {{\alpha }_{j}}-\alpha _{j}^{*} \right)}\cdot K\left( x-{{x}_{i}} \right)+ \\ \sum\limits_{i=1}^{n}{\left( {{\alpha }_{i}}-\alpha _{i}^{*} \right){{y}_{i}}-\sum\limits_{i=1}^{n}{\left( {{\alpha }_{i}}+\alpha _{i}^{*} \right)\varepsilon ;}} \\ \end{matrix} $

$ \text{s}\text{.t}\text{.}\left\{ \begin{matrix} \sum\limits_{i=1}^{n}{\left( {{\alpha }_{i}}-\alpha _{i}^{*} \right)\left( {{\alpha }_{j}}-\alpha _{j}^{*} \right)=0,} \\ 0\le \alpha _{i}^{*}\le C,0\le {{\alpha }_{i}}\le C, \\ 0\le \alpha _{j}^{*}\le C,0\le {{\alpha }_{J}}\le C。 \\ \end{matrix} \right. $

(5)

式中：ε为误差； α_i，α_j，α_i^*，α_j^*(i=1，2，…，n；j=1，2，…，n)为lagrange乘子； C为惩罚参数，用于平衡模型复杂度和损失误差； K(x，x_i)为核函数，x_i为核函数中心，x为因变量；α^*=[α₁^*，α₂^*，…，α_n^*]^T为最优解。

本研究采用LIBSVM工具包实现SVR算法在区域森林冠层高度估测中的应用，具体步骤如下：

1) 利用已建立的GLAS森林冠层高度估测模型[式(4)]计算研究区所有GLAS 光斑内的森林冠层高度。

2) 利用径向基核函数K(x，x_i)将GLAS估测的森林冠层高度与相应的3个MNF分量构成的样本数据集映射到高维特征空间中，并将样本数据随机分成训练集和测试集。K(x，x_i)表达式如下：

$ K\left( x-{{x}_{i}} \right)=\exp \left( -\frac{\left\| x-{{x}_{i}} \right\|}{2{{\sigma }^{2}}} \right)。 $

(6)

式中： σ为函数的宽度参数，用于控制函数的径向作用范围。

3)应用LIBSVM工具包实现SVR算法的关键是确定惩罚参数C和核函数参数σ。本研究采用遗传算法寻找最优的C和σ。遗传算法是通过模拟生物进化过程中自然选择和遗传变异的一种随机优化方法。由于遗传算法的运算对象是表示个体的符号串，所以首先通过编码将实际问题转换到遗传空间，用染色体代表实际问题； 然后，给定参数的初始寻优范围(本研究设定C和σ 的寻优范围分别为[0.01，1 000]和[0.01，1 000]，设置种群最大数量为200，种群交叉概率为0.4，交叉验证次数为4)进行参数寻优，根据结果调整参数范围继续寻优直至找到最优结果。

4)利用得到的最优参数C和σ 及训练集数据建模，并用测试集数据对所建模型的泛化能力进行测试，得到最优SVR模型。

5)将研究区MNF1，MNF2和MNF3栅格数据导入到SVR模型中，计算整个研究区各像元内的森林冠层高度，生成森林冠层高度分布图。

2.5 模型精度评价指标

为了评价模型的性能，验证最终得到的森林冠层高度分布图的准确度，本研究采用决定系数(R²)和均方根误差(RMSE)作为评价指标。各评价指标表达式如下：

$ {{R}^{2}}=\frac{\sum\limits_{i=1}^{N}{{{\left( {{f}_{i}}-\overline{y} \right)}^{2}}}}{\sum\limits_{i=1}^{N}{{{\left( {{y}_{i}}-\overline{y} \right)}^{2}}}}; $

(7)

$ \text{RMSE}=\sqrt{\frac{1}{N}\sum\limits_{i=1}^{N}{{{\left( {{y}_{i}}-{{f}_{i}} \right)}^{2}}}} $

(8)

式中： f_i为森林冠层高度估测值； y_i为森林冠层高度实测值； y为森林冠层高度实测值的平均值； N为样本数。

3 结果与讨论 3.1 GLAS森林冠层高度估测模型

本研究基于GLAS 波形参数(W和TS)建立森林冠层高度估测模型，并利用此模型计算研究区所有GLAS 光斑内的森林冠层高度。基于W和TS所建森林冠层高度估测模型的R²=0.78，RMSE=2.51 m，模型验证的R²=0.85，RMSE=1.67 m。模型表达式如下：

$ H=0.51W+0.04D\tan \text{TS}+4.45 $

(9)

图 5所示为模型验证散点图，可见实测值与估测值之间存在偏差。本研究对偏差较大的样本数据中对应的ICESat-GLAS波形数据及野外数据进行对比分析，发现这些样地植被覆盖度较高，林下植被丰富，地形较复杂，由于枝叶的遮挡，可能会导致野外实测森林冠层高度误差增大；此外，林下植被对GLAS波形也会造成一定的影响，进而影响森林冠层高度估测。

图 6所示为考虑地形坡度与否森林冠层高度估测误差随坡度的变化情况。分析图 6发现，不考虑地形坡度，直接将GLAS波形长度W作为森林冠层高度的估测误差明显高于本研究所建模型的估测误差，而且随着地形坡度的增加，估测误差随之增大。此外，本研究提出的模型估测误差随地形坡度变化的趋势不明显，说明地形坡度参数TS的加入削弱了地形坡度对森林冠层高度估测的影响。

针对此研究区，Xing等(2010)利用地形指数修正地形坡度对波形的影响，建立了森林冠层高度对数模型，该模型在0~5°范围内的R²=0.92，RMSE=2.56 m，随着坡度范围增大，其精度逐渐下降，当坡度范围扩大到0~30°时，其模型的R²=0.56，RMSE=4.23 m。该模型虽然在0~5°坡度范围内精度很高，但模型中森林冠层高度与波形长度间为对数关系，而理论上在平坦地区，波形长度近似等于森林冠层高度，是线性关系，所以该模型在平坦地区不好解释。而本研究所建模型为线性模型，克服了模型解释困难的问题；且模型是基于各个坡度范围样方数据建立的，可应用于不同坡度范围。此外，本研究从ICESat-GLAS波形数据中提取地形坡度参数TS，减少了借助DEM数据计算地形指数所造成的误差。

3.2 森林冠层高度分布图

针对GLAS数据在森林冠层高度区域估测方面的不足，本研究将ICESat-GLAS波形数据与HJ-1A/HSI高光谱数据联合，基于SVR算法，利用GLAS估测的森林冠层高度与从HSI 高光谱数据中得到的3个MNF分量(MNF1，MNF2，MNF3)建立了森林冠层高度SVR估测模型，此模型建模的R²=0.70，RMSE=3.62 m，模型验证的R²=0.67，RMSE=4.42 m，模型最优参数C=1.976 6，σ =0.805 1。图 7所示为最终得到的森林冠层高度分布图，表 1为森林冠层高度估测误差统计分析结果。

分析图 7可以看出，研究区内森林冠层高度在0~35 m范围内，研究区的北部和西南部森林冠层高度较高，而在研究区的南部、东北部和东南部森林冠层高度较低。结合野外调查数据分析发现，研究区的北部和西南部，生长的多是常绿乔木，而在研究区的南部、东北部和东南部，生长的多是落叶松(Larix gmelinii)。本研究采用的高光谱数据获取时间为6月和10月，在10月落叶松已经开始变黄落叶，因此导致10月高光谱数据获取的光谱信息存在误差，使得最终得到的森林冠层高度偏低。

从表 1可知，森林冠层高度估测误差的最大值为7.10 m，最小值为0.07 m，平均值为1.78 m，估测误差的标准差为1.49 m，中值为1.45 m，Q1为0.75 m，Q3为2.31 m，说明有25%的估测误差≤0.75 m，50%的估测误差在0.75~2.31 m之间。综上分析可知，本研究提出的方法对于区域森林冠层高度估测是可行的，能够满足林业应用需求。

本研究中联合ICESat-GLAS波形数据和HSI高光谱数据估测森林冠层高度的误差来源主要有2方面： 1) 客观原因导致的数据获取时间上的差异。在利用ICESat-GLAS波形数据估测森林冠层高度时，采用的都是2006和2007年的数据，虽然不是同一个月份的数据，但在同一年内森林冠层高度之间的差异不明显，数据获取时间差异对GLAS估测森林冠层高度的结果影响不大。但引入HSI高光谱数据进行区域推广时，采用的HSI数据是2010年的，与ICESat-GLAS波形数据和野外数据间存在3~4年的时间差，这期间树木会有所增长，因此导致最终得到的区域森林冠层高度比实际值偏低。2) 地形因素。GLAS 光斑直径约为70 m，光斑覆盖面积大，光斑内地形坡度变化会造成波形展宽、地面回波与植被回波重叠等现象。虽然本研究从ICESat-GLAS波形数据中提取了TS参数对地形坡度进行了修正，但是不可能完全消除其影响。

吴迪等(2014)联合GLAS波形数据和多角度光学遥感数据，基于随机森林法建立了塔河林业局营林区的森林冠层高度分布图(R²=0.72，RMSE=1.83 m)，精度稍高于本研究，但是其只是利用坡度小于10°范围内的数据进行森林冠层高度制图，没有考虑坡度较大的区域；而本研究利用研究区0~35°坡度范围内的数据进行分析建模，最后生成森林冠层高度分布图，因此本研究结果对于地形复杂的区域更具可行性。

4 结论

为了提高ICESat-GLAS数据估测森林冠层高度的精度，弥补ICESat-GLAS数据估测区域森林冠层高度的不足，本研究从ICESat-GLAS波形数据中提取波形参数(波形长度W和地形坡度参数TS)建立GLAS森林冠层高度估测模型，并基于SVR算法建立了GLAS估测的森林冠层高度与从HSI高光谱数据提取的3个MNF分量的森林冠层高度SVR估测模型，估测了无GLAS光斑覆盖区域的森林冠层高度，生成森林冠层高度分布图。通过分析可得到以下结论：

1) 从ICESat-GLAS波形数据中提取的地形坡度参数TS与地形坡度间存在很强的线性关系。基于波形长度W和地形坡度参数TS建立的森林冠层高度估测模型，减少了借助DEM数据计算地形指数所造成的误差，模型建模的R²=0.78，RMSE=2.51 m，模型验证的R²=0.85，RMSE=1.67 m。

2) 本研究借助HJ-1A/HSI高光谱数据提供的连续高分辨率光谱信息，基于SVR算法，建立GLAS估测的森林冠层高度与HSI高光谱前3个MNF分量森林冠层高度的SVR估测模型，实现了森林冠层高度的区域制图，提高了森林冠层高度估测的精度，模型建模的R²=0.70，RMSE=3.62 m，模型验证的R²=0.67，RMSE=4.42 m。

由于客观原因的存在，导致获取的数据存在时间差异，影响了结果的准确性。本研究虽然考虑了地形坡度对森林冠层高度估测模型的影响，但是并不能完全消除其影响，如何彻底解决这些因素的影响尚需进一步研究。