森林是人类赖以生存的自然资源之一，与气候变暖、碳排放、碳循环等热点话题有着千丝万缕的联系。森林资源调查是世界各国了解森林植被变化必不可少的手段。传统的森林资源调查对相关森林参数如树种、树高、胸径、冠幅、郁闭度、地形坡度等进行野外实测，再计算蓄积量、森林生物量等森林物理生态参数，但实测调查不仅耗时费力投入高，而且缺乏时间和空间尺度上对森林参数的连续监测。由于激光脉冲可以穿透森林冠层间隙反映出林下结构及地表信息，激光雷达遥感技术在探测植被的垂直结构信息等森林植被参数估算及森林生态监测中的应用日益引起关注^[1-2]。目前，机载激光雷达和星载激光雷达数据在森林结构参数估测中有着广泛的应用^[3]。国内外学者利用激光雷达数据在反演森林冠层高度，冠层垂直结构分布、森林分类、郁闭度、森林生物量^[4-7]等方面作出了许多有益的研究，但现有的研究仅使用了波形数据的部分信息，更多提取信息的方法还有待探讨。

机载激光雷达多为小光斑系统，光斑尺寸小于林木冠幅，主要是利用激光点云估算单木的结构参数，数据获取成本高，难以实现大范围的数据获取，制约了其在大区域的广泛应用。星载激光雷达采用卫星平台，例如搭载于冰、云与陆地抬升卫星(Ice, Cloud, and Land Elevation Satellite, ICESat)卫星上的地球科学激光测高仪系统(Geoscience Laser Altimeter System, GLAS)载荷是目前唯一对地观测的大光斑星载激光雷达系统，具有运行轨道高、观测范围广等优点，可实现全球范围覆盖，在估算大尺度森林参数方面具有无可比拟的优势。中国2016年下半年立项的陆地生态系统碳监测卫星将搭载面向森林目标的全波形记录激光雷达载荷，标示着星载激光雷达技术在中国森林资源调查中将得到更深入、更全面的应用。

建立面向森林植被的大光斑激光雷达回波仿真器，可以更好地理解森林植被空间结构和光学特性对激光雷达回波波形的影响，研究目标植被结构参数与回波波形参数之间的联系，为波形信息提取、参数反演及激光雷达载荷设计提供理论依据。现有的面向森林目标的激光雷达回波仿真模型主要有3类:(1)文献[8]的模型针对稠密森林，只考虑植被表层起伏分布和单次反射，并假设植被反射率相同，没有考虑复杂的森林空间结构；(2)文献[9]的模型结合森林场景中林木冠层的空间结构和辐射传输模型而构建; (3)文献[10-11]的模型是利用机载激光雷达提供的森林冠层几何参数构建间隙率函数，采用混合几何光学辐射传输模型(geometric optical and radiative transfer, GORT)模拟激光雷达的回波。后两类模型都考虑了激光在植被冠层内部的多次散射，但是对地表只做了无起伏有坡度的假设，并且没有考虑林下植被覆盖的情况。

本文构建了基于统计规律的参数化森林植被三维空间分布模型，利用植被冠层的辐射传输理论得到适用于星载激光雷达的森林植被反射模型，再根据星载激光雷达的回波理论，得到面向植被目标的激光雷达回波仿真器。本文在仿真过程中讨论了单层林和具有林下植被覆盖的森林结构，考虑了地形坡度及地表起伏变化的影响；并通过仿真大光斑的激光雷达回波对坡度变化的敏感性，发现适当减小光斑尺寸，可以降低坡度对回波的影响，更有利于反演植被目标信息。

1 星载激光雷达的植被回波仿真理论

不考虑大气环境影响，星载激光系统的发射脉冲经过一次菲涅耳衍射到达目标植被，与目标作用发生二向反射，再一次经过菲涅耳衍射传回接收望远镜视场，经过光电转换和滤波，得到输出回波脉冲信号。经目标表面反射后到达望远镜的回波信号的功率p(t)表示为^[12-13]：

式中，$\psi = \frac{{2R}}{c} + \frac{{{\rho ^2}}}{{cR}}{\rm{ - }}\frac{{2h\left( \rho \right)}}{c} $；∑表示光斑的区域范围；A_R表示接收望远镜的面积；T_a表示大气衰减系数；η表示接受系统能量利用率；R表示测高系统的卫星轨道高度；c表示光波在真空中的传播速度；β(α, ρ)表示光斑内目标的表面反射率；α表示光束的入射角；ρ表示被测目标的横截面坐标半径；h(ρ)表示光斑内目标表面的轮廓分布；f(t)和a(ρ, R)分别表示发射脉冲的时间分布和空间分布。

一般认为，发射激光脉冲在时域和空域都近似满足高斯分布，可表示为：

式中，F(t)、I(ρ, R)分别表示发射激光脉冲的时域分布和空域分布；θ_T为光束发散角；E_t为发射激光的脉冲能量；δ_t为发射激光脉冲的均方根脉宽。根据高斯函数的分布特性，若要更完整地记录发射激光脉冲，应使发射脉冲在空间分布上的光斑半径和时间脉宽分别取3Rtanθ_T和3δ_t，即可包含其99%以上的能量。

对于森林植被目标，在系统参数不变和环境因素恒定的情况下，植被仿真只关心植被参数与回波波形特性的相互作用，ψ中的第一、二项可忽略，则式(1)可归一化为：

森林植被目标轮廓特性由h(ρ)给出，反射率由β(α, ρ)表示，积分号中的部分是森林目标响应函数，由发射脉冲的空间分布和目标的轮廓分布及反射特性决定。回波信号由发射脉冲的时间分布与目标响应函数的卷积给出。森林植被的目标响应函数很难用解析的方式表达，只能用离散的数值来描述。

2 森林植被的目标模型

由于林木空间结构复杂多变，森林区域的精确建模非常困难。对于以GLAS为例的星载激光雷达的大尺度遥感观测，植被目标的反射模型建模时，不考虑具体到一枝一叶这样的微观尺度，通过分析林木空间结构的统计规律建立森林植被的三维空间分布数据模型，进而得到与实际森林样地数据无偏的森林空间结构数据。

2.1   林木的空间几何模型

影响激光雷达回波波形的森林植被的空间结构参数主要有林木的位置、树高、树种、树冠形状和大小。首先，按照林木的空间格局生成水平位置分布。由于自然环境影响和植被种群生长竞争，不同种群的分布具有不同的空间格局，大致可分成均匀分布、聚集分布和随机分布3种基本的分布类型^[14]。依据不同分布格局生成光斑面积内的森林植株水平位置坐标。以阈值距离d为参数，控制林木的最小间距。其次，在生成林木水平位置坐标的基础上，根据植株距离与高度的统计规律得到相应的树高分布。一般认为，林分中林木的树高是正态分布，单株林木的高度与邻近林木相关性强，与距离远的林木相关极小^[15]。因此，先对实测的样地数据进行统计分析，建立反映林木间相互关系的参数模型，林木高度分布的相关性由协方差估计公式来描述；再由Cholesky分解方法建立树高空间分布的仿真模型，获得一组对应水平位置坐标上与样地数据无偏的虚拟树高数据。改变树高均值、方差、树间距等参数，可以获得不同分布的森林植被的垂直树高分布数据。利用这一方法获得的森林树高分布既可以反映真实森林的树高分布特点，又减少了外业实测的工作量。

通常情况下，林木的冠幅和冠层高度都与林木的胸径正相关。通过拟合林业实测数据中的树高与胸径的关系，由树高分布得到对应的林木胸径值，进而确定冠幅和冠层高度。对不同树种采用椭球体、圆锥体、球体或倒圆锥体来对树冠的形状进行近似处理，最终形成不同树种按照树高、冠幅形状和大小变化的森林植被参数化模型，得到目标轮廓分布的离散表示h(ρ)。

2.2   森林植被反射特性

按照辐射传输理论和平均冠层透射理论，面向森林植被的星载激光雷达遥感属于大尺度遥感，森林植被冠层可看做是浑浊散射体，假定冠层可用水平无限均匀介质近似构成，整个冠层又可垂直分成若干薄层，每层中植被的结构特性参数和光学特性参数均可近似为常数。图 1是森林植被目标的空间结构模型示意图。对于图 1(a)所示的森林模型，沿林木高度方向分成若干等厚均匀薄层，每层中的元素再进行网格划分。如图 1(b)所示，在每一个小网格中可能是树叶，激光脉冲分量被反射，部分透射；可能是林木间隙，则激光脉冲分量透过，不发生反射，直到在下一层中遇到树叶部分再计算反射能量；或者是地表，按漫反射来计算。由于星载激光雷达是近天顶入射和接收的测量机制，故只考虑发射脉冲垂直入射的情况。在星载激光雷达光束近天顶入射时，激光脉冲的时间变量t与植被高度变量z满足t=2(R－z)/c，其中，R为卫星轨道高度，c为光速。所以在仿真时林地的垂直树高分辨率与发射脉冲的时间分辨率要一致。以GLAS为例，发射脉冲时间分辨率为1 ns时，树高的垂直分辨率为0.15 m。

图  1  森林植被的模型示意图

Figure  1.  3D Model of Forest Vegetation

下载: 全尺寸图片 幻灯片

把植被冠层当作浑浊散射体，则其光学性质要用二向性反射特性来描述。由于星载激光雷达的激光发射装置和接收装置同在天顶方向，满足热点条件^[16]。对于植被冠层，二向间隙率(bidirectional gap probability, BDGP) p(z)可以表示为：

式中，z是冠层高度；u_L(z′)是深度为z′处的单侧叶面积体密度；G(z′)是Ross-NilsonG函数，即在入射方向上单位叶面积在z′深度投影面积的平均值。

假设树冠在一个水平薄层里是均匀的，故冠层厚度为Δz、表面积为Δs的单层树冠的一次散射的二向反射率表示为：

式中，Γ=GR_L是树冠的面散射相位函数。将单个叶片看作朗伯反射体，在树冠某高度的水平薄层里，其反射率R_L、G函数和单侧叶面密度u_L在树冠体内均为常数。按照树高和树冠形状的设定计算每棵树对应网格点的反射率，高度为z处的冠层反射率是厚度为Δz的薄层内各个树冠反射率分布的集合。在仿真多树种混交林时，利用椭球体、圆锥体、倒圆锥体、球体等几何形状对树冠形状进行近似，将不同的Г、u_L、G代入式(6)，ds的积分域变为树冠在该层上的投影面积，可计算相应的等效反射率β_i(z)。

由于到达林下地表的直射和散射光只占总辐射的一小部分，地表散射的大部分光还要经过植被的多次散射才能逸出植被。对于不同类型的林下地表，包括有部分枯枝落叶的土壤或有小绿叶灌木覆盖的土壤，其反射特性则更加复杂。故模型中设置林下植被层，通过调节林木高度和方差来模拟二层林或小绿叶灌木覆盖层，输入覆盖率参数表示地表覆盖物的比例。当表示二层林时反射率的计算同前，当表示小灌木等覆盖物时，因为可能存在枯枝落叶，反射率适当调低。研究结果表明，土壤的反射特性由二向反射来描述，当入射光垂直入射或近垂直入射时，土壤的反射因子近似为常数。本文中设置抛物面面型和地表粗糙度来模拟地表的随机起伏变化，根据坡度的设定把土壤反射率赋值到对应的网格点，并计算每株林木的冠层等效反射率，赋值到三维网格的对应格点，就可以建立起完整的具有反射特性的植被模型，包括植株位置、高度、树高、树冠形状及冠层反射率和地表反射率等参数。

3 回波信号波形仿真及对比

一般将发射激光脉冲近似为高斯脉冲，按GLAS系统参数设置发射脉冲的参数^[17]，即轨道高度z=600 km，激光波长λ=1 064 nm，光束发散角θ_T=110 mrad，发射能量E_t=75 mJ，发射脉冲半峰全宽为4 ns，接收系统的时间分辨率为1 ns，生成发射激光脉冲的时空分布模型。根据仿真森林区域位置设置森林的分布类型，如人工林一般为均匀分布，自然林一般为聚集分布。按照表 1中的参数，设定林木棵树、阈值距离、森林区域边长生成林木的水平位置分布。设置树高和树高方差依据植株距离与高度的统计规律，生成林木的垂直空间分布。通过拟合得到树高、胸径和冠层冠幅的关系，计算冠幅几何尺寸，根据仿真区域的树种特点，可以用椭球体、圆锥体、球体或倒圆锥体来近似树冠形状。把生成的林木三维空间作网格化划分，如图 1所示，对每一棵树从冠顶向下按式(6)计算每个格点的反射率，遍历所有树，生成冠层反射率集合。若有二层林或地面灌木覆盖物，就按前述生成林下植被模型。地形面型可选平面、带坡度的斜面或抛物面，地表的起伏用目标粗糙度近似，满足均值为0、方差可调的正态分布，假设地表反射率为常数(见表 1)，生成地形模型。冠层向下投影，从上往下覆盖，无阴影处即保留地表信息，如此3个模型对应叠加，最后得到图 1所示的三维网格中每个格点赋予反射率值的林木冠层及林下结构和地表的3D森林植被反射模型。

表  1  植被仿真模型输入参数

Table  1.  Input Parameters of Vegetation Model

参数	值	参数	值	参数	值
样地直径/m	66	冠层反射率(RL)	0.55	水平分辨率/m	0.2
平均树高/m	15.7	地表反射率(RG)	0.29	垂直分辨率/m	0.15
树高方差/m2	3~9	单侧叶面积体密度(uL)	0.8	G函数	0.5
阈值距离/m	1	树冠形状	椭圆体，圆锥体	地表斜率/(°)	0~15
灌木层平均树高/m	0.5~6	灌木层方差/m2	0.3~3	地表方差/m2	0.3~1

下载: 导出CSV 

| 显示表格

以大兴安岭地区(50.68°N~51.32°N，120.71°E~121.42°E)样地的实测数据分析统计规律，生成森林模型。由于森林调查样地的选取本身也是在林分中按照比例进行抽样设置，由样地测量结果近似代表林分的真值，因此该区域的森林分布特点可以由样地结构特征表示。按统计规律生成多个与实际森林样地数据无偏的森林空间结构，一个激光光斑覆盖的区域模拟一个实测样地，仿真回波波形，与该区域的GLAS回波做相关性分析。森林仿真输入参数见表 1。

图 2和图 3所示分别为GLAS的L3I工作时段经过大兴安岭上述区域的部分脚点轨迹和回波波形。以样地参数为依据，样地附近的森林结构是类似的，每个光斑内的坡度会有一定差异，调整光斑内的坡度，经大量仿真波形计算，仿真波形与GLAS回波的相关系数R²均达到0.91以上，高度相关。图 4所示为森林植被目标回波模拟器仿真得到的波形与GLAS回波波形的对比。GLAS数据是以UTC-J2000为时间标准按时间序列记录的，即以2000年1月1日午时为时间轴起点，图 3和图 4中的UTC表示GLAS回波数据的时间节点。

图  2  GLAS脚点与实测样地的位置示意图

Figure  2.  GLAS Footprint Track and Position of Field Measured Sample Area

下载: 全尺寸图片 幻灯片

图  3  大兴安岭地区样地附近的部分GLAS回波波形

Figure  3.  Return Waveforms over Greater Khingan Mountains Region from GLAS

下载: 全尺寸图片 幻灯片

图 4中所示的GLAS回波波形有展宽的单波峰，也有多波峰的波形。模型中根据实际目标森林结构，考虑了二层林的存在和地面灌木覆盖等因素，使得仿真回波更接近真实回波。对应目标区域森林结构的回波，本文建立的植被波形仿真器模拟的回波波形都与之符合得很好，图 4中所示的仿真波形与真实回波的相关系数R²达到0.97以上。这说明根据实际目标森林的先验知识，经过树高、方差、坡度、树冠形状、地表覆盖等参数的调整，本文建立的仿真器可以很好地仿真星载激光雷达的回波波形。

图  4  大兴安岭地区的仿真波形与GLAS回波波形对比图

Figure  4.  Simulated Waveforms vs. GLAS Echo Waveforms over Greater Khingan Mountains Region

下载: 全尺寸图片 幻灯片

4 光斑大小对坡度地形植被回波的影响

目前唯一的对地观测星载激光雷达GLAS的脚点光斑直径约为66 m，远大于单株林木的树冠冠幅，一般情况下，即使林木枝叶稠密，植株间距也会使得光斑内含有地面信息。然而，山区林地的回波很多表现为展宽的单峰波形，地形坡度增大会导致回波脉宽的展宽与叠加，甚至使波形中的地面回波和冠层回波很难清楚地辨别出来。

保持林木输入参数如表 1不变，不考虑林下植被，选取树冠形状为椭球体，将坡度作为变量，对不同坡度地形的林地的回波波形进行仿真，如图 5所示。分别以GLAS的系统参数为输入得到图 5(a)，以改变激光束发散角使地面激光光斑直径等于30 m为参数得到图 5(b)的仿真波形。

图  5  地形坡度对不同直径光斑回波的影响

Figure  5.  Influence of Terrain Slope on the Different Size of Echo Footprint

下载: 全尺寸图片 幻灯片

由图 5中的曲线可看出，随着坡度的增加波形脉宽展宽，脉冲峰值能量下降，波形的时间长度增加。对于直径66 m的光斑，当坡度小于5°时，回波波形的展宽和混叠不明显，可以明显区分冠层的回波波形与地表的回波波形；随着坡度增大，可以观察到明显的波形展宽，当坡度达到15°时，已经不能区分冠层回波和地面回波。而对于直径30 m的光斑，坡度为15°时，依然能分辨出冠层回波和地面回波。即对于同样的树高分布，缩小光斑尺寸，回波对于坡度的敏感性也降低。

对于直径为66 m光斑的系统，地形坡度较小时，当光斑内地面高程最高高度在h₂时，冠层回波和地面回波可以区分，参看图 6右侧的第2个波形示意图；当地形存在较大坡度、光斑内地面最高高度在h₁、发射脉冲入射到h₁高度时，该高度层反射的被探测器接收到的信号中既包含树冠顶层反射贡献的能量，也包含地面反射贡献能量，树冠和地面回波将融合在一起很难被区分，参看图 6右侧的第1个波形示意图。

图  6  坡度对回波的影响示意图

Figure  6.  Influence of Slope to Return Waveform

下载: 全尺寸图片 幻灯片

而对于30 m光斑的系统，同样坡度的地形，光斑内的地表最高高程点相对植被冠层高度下降了很多，地面回波与植被冠层回波不发生混叠，可以分辨植被冠层回波和地面回波，参见图 6右侧第2个波形示意图。由此可见，坡度导致的回波展宽混叠，不仅是由于坡度引起的，混叠的程度与光斑尺寸及树高都有关系。对于同样的森林树高，相对小的光斑尺寸对坡度的敏感度小，可以反演更大坡度的植被信息，即减小发射脉冲的发散角，将有利于更多山区大坡度森林植被回波波形信息的提取。

5 结语

对于森林植被目标，激光雷达回波模型的建立相比普通陆地目标更为复杂，森林植被目标的空间结构参数多，植被的反射特性与植被空间结构相关。本文根据现有的森林测量数据分析其统计规律，按统计规律建立森林植被反射模型，通过调整设置参数可以获得不同结构的林分仿真模型，适应性较强。在此基础上仿真了大量大兴安岭地区的星载激光雷达森林植被回波，与GLAS回波对比达到高度相关。本文通过仿真大光斑星载激光雷达的森林回波波形，分析了坡度对回波波形的影响，其与光斑尺寸和树高相关。在可能范围内，减小发射脉冲的发散角或降低轨道高度，都将提高回波波形对森林目标特性的表达，有利于目标参数的反演。本文所建立的植被回波仿真器能较好地反映星载激光雷达回波波形特性，利用回波仿真波形分析光斑大小对于地形坡度影响的敏感性，对中国未来研制和发射的激光雷达载荷系统设计有一定的参考意义。