国产星载激光雷达森林回波波形模拟仿真

蔡龙涛 岳春宇 黄 缙 贺 涛 国爱燕 邢艳秋

(1.东北林业大学工程技术学院，哈尔滨 150040；2.北京空间机电研究所，北京 100094；3.中国空间技术研究院遥感卫星总体部，北京 100094)

0 引言

森林生态系统是陆地生物圈中占地面积最大、结构组成最为复杂以及物质资源最为丰富的生态系统。森林面积约3.815×109hm2，占全球陆地总面积的25.60%[1]，生物量约占地球陆地生态系统总生物量的90%[2]，碳储量约占全球碳储量的45%[3]。森林生态系统在抑制全球变暖、降低碳排放和促进碳循环方面有着重要作用。为确定2030年和2060年是否能完成碳达峰和碳中和任务，需对国内森林资源进行高精度动态监测。

森林资源调查是了解森林资源动态变化的有效手段，传统森林资源调查费时费力，且难以实现大区域尺度森林生物量观测研究[4-6]，目前多采用激光雷达进行森林资源调查。而激光雷达作为一种主动遥感技术，其发射器发射的激光脉冲具有较强的穿透力，可穿透森林冠层以探测林下结构和地表信息，实现森林植被信息的动态观测。目前，已接收星载激光雷达回波波形的对地观测卫星有ICESat1(Ice,cloud,and land elevation satellite)、GF-7和GEDI(Global ecosystem dynamics investigation)，后续即将发射的星载激光雷达波形类对地测高卫星为陆地生态系统碳监测卫星[7]和MOLI(Multi-footprint observation LiDAR and imager)卫星[8]。国外星载激光雷达波形类卫星如ICESat-GLAS(Geoscience Laser Altimeter System)波形数据已用于冰川探测[9-11]、水位探测[12]、地物分类[13-15]和估测森林资源[16-17]；GEDI卫星波形数据同样在树高估测[18-20]、植被覆盖度[21]和森林生物量估测[22-24]方面取得了一定进展。而国产星载激光雷达波形类卫星还未发射，目前该卫星尚处于技术攻关阶段，其波形数据应用潜力有待挖掘。

为探究国产星载激光雷达波形数据在森林资源调查方面应用潜力，对其回波波形进行模拟仿真是其中重要的环节。李松等[25]对星载激光雷达回波波形模拟仿真时把地形分为斜坡地形和阶梯地形。潘浩等[26]对回波波形模拟仿真时发现冠层回波与地面回波重叠程度随地形坡度增大而增大。庞勇等[27]发现随着地形坡度增大，仿真波形中地面波峰和植被波峰值随之降低,且冠层回波与地面回波发生信息混叠。其他研究多基于机载点云数据对星载激光雷达回波波形进行模拟仿真[28-29]。然而，上述研究并未从地形起伏特点方面对地表进行模拟仿真，且回波仿真过程中未考虑激光脉冲在大气传输过程中的能量衰减问题，以及探究多种因素(如地形坡度、郁闭度和森林类型)对回波波形仿真精度的影响。

本文依据有限元原理，基于林地地形随机分布特点，建立随机地形；考虑到激光脉冲在大气传输过程中能量的衰减，波形仿真过程中加入激光雷达辐射传输模型；利用回波仿真原理[30]，分别对GLAS发射波，不同地形坡度、郁闭度和森林类型条件下回波波形进行模拟仿真。依据GLAS实测波形与仿真波形相关性分析结果，对所建回波仿真系统有效性进行验证。然后利用本文建立的回波仿真系统，对国产星载激光雷达回波波形进行模拟仿真。

1 研究区概况和研究方法

1.1 研究区概况

研究区为吉林省汪清林业局经营区(图1)，该区域属于长白山系中低山区(43°5′～43°40′N,129°56′～131°4′E)，地处寒温带，总面积3.04×105hm2，南北长约60 km，东西长约85 km，地面高程变化范围为360～1 477 m，地形坡度变化范围为0°～45°。

图1 研究区位置及野外样地分布Fig.1 Location of study area and distribution of field sampling plots

林区内森林覆盖率达到95.95%，深山区林分以针叶林、阔叶林和混交林为主，带状分布于海拔500～1 100 m之间。针叶树主要有红松(PinuskoraiensisSieboldetZuccarini)、云杉(PiceaasperataMast.)和臭冷杉(Abiesnephrolepis(Trautv.)Maxim.)，阔叶树多为椴树(TiliatuanSzyszyl.)、蒙古栎(QuercusmongolicaFischerexLedebour.)、枫桦(BetulacostaTrautv.)、色木槭(AcermonoMaxim.)和白桦(BetulaplatyphyllaSuk.)等。

1.2 研究数据

1.2.1ICESat-GLAS波形数据

为验证本文建立的回波仿真系统有效性，本文选用应用较为成熟的ICESat-GLAS波形数据，对GLAS发射波仿真波形和GLAS回波仿真波形进行相关性分析。结合相关性分析结果验证本文所建回波仿真系统有效性。

ICESat-1/GLAS是第1个极地轨道大光斑激光雷达卫星，该卫星共提供15个数据产品：GLA01～GLA15。其中，GLA01数据产品记录了GLAS发射波和回波波形数据；GLA14数据产品记录了GLAS波形数据对应的地面光斑地理位置和高程数据。GLAS波形数据可从美国国家冰雪数据中心(http:∥nsidc.org/data/ice-sat/)下载，美国国家冰雪数据中心拥有2003年至2009年采集的所有ICESat-GLAS回波波形数据。GLAS载荷参数如表1所示。

表1 ICESat-GLAS回波仿真相关载荷参数Tab.1 Load parameters of ICESat-GLAS for echo simulation

1.2.2国产星载激光雷达载荷参数

目前，国产星载激光雷达所搭载的卫星还未发射，无法获取国产星载激光雷达实测波形。为探究国产星载激光雷达回波波形在森林结构参数估测方面应用潜力，本文基于国产星载激光雷达载荷参数和实地调查数据，结合回波仿真原理[30]对国产星载激光雷达回波波形进行模拟仿真，以获取国产星载激光雷达仿真波形。国产星载激光雷达回波波形模拟仿真过程涉及到的载荷参数如表2所示。

表2 国产星载激光雷达载荷参数Tab.2 Load parameters of domestic spaceborne LiDAR

1.2.3实地调查数据

本文利用分层随机采样法，选取286组光斑数据，其中森林样地251组，其他类型样地如水地、裸地和草地等共计35组。分别于2006年9月、2007年9月和2010年9月3次采集获取。其中，实地调查样地点位分布图如图1所示。在实地勘测过程中，以针叶林、阔叶林和针阔混交林作为研究对象，使用GPS对已选定的激光光斑采样点进行定位，把验证数据对应光斑中心点作为地面调查样地的圆心，依据公式πR2cosθ=500 m2，建立水平投影面积为500 m2的圆形样地，记录样地内植被分布情况、植被类型和植被覆盖度、样地坡度与对应样地半径。另外，对针叶林、阔叶林和混交林分类时，主要结合我国森林资源调查主要技术规定，将针叶林蓄积量占总蓄积量65%以上的样地定义为针叶林，阔叶林蓄积量占总蓄积量65%以上的样地定义为阔叶林，任何一个树种蓄积量占总蓄积量不到65%的样地定义为混交林。

1.3 回波仿真相关原理

1.3.1基于有限元原理构建地表响应函数模型

本文所选林分为天然林，林分内树种多样，林层结构复杂，难以实现光斑内林分信息三维模拟仿真。本文依据有限元原理，以光斑中心为原点，自原点出发，0.05 m为间隔把正圆划分为不同直径的同心圆。然后以原点为中心，从极坐标0°开始，以3°为间隔画直线把光斑等分成120个扇形。直线和同心圆共同把光斑划分为72 000个小区域，如图2所示。

图2 光斑划分示意图Fig.2 Schematic of spot division

星载激光雷达光斑划分为若干扇形小区域后，在垂直方向上以0.15 m为间隔对光斑内地物进行垂直分层。假定在垂直方向上可分为m层，每层有k个小区域，统计每层(第j层)小区域内地物在地面的投影面积，并将所有层地物投影面积按时间序列排列起来，可构成地表响应函数模型，计算式为

(1)

G(t)={N(t(1),t(2),…,t(m))}

(2)

式中t(j)——发射波与第j层地面目标物接触时的时刻

N(j)——第j层地面目标物在地面的投影面积

Ni——第j层第i个小区域在地面投影面积

G(t)——地表响应函数

1.3.2回波仿真模型

依据回波仿真原理[30]，星载激光雷达回波波形为发射波函数与地表响应函数(后向散射截面的集合)的卷积。若不考虑大气影响，回波仿真模型定义为

E(t)=F(t)*G(t)

(3)

式中E(t)——回波波形

F(t)——激光发射脉冲函数

*——卷积运算符

本文所用样地在地面投影为直径25.2 m的正圆，而验证数据对应光斑直径为70 m，两者覆盖区域面积相差较大，不同区域内林木个数同样存在较大差异。为高精度模拟仿真国产星载激光雷达回波波形，需对直径25.2 m样地林木个数进行扩展，在保证林木密度、森林类型和郁闭度相同的条件下把光斑直径扩展到70 m，近似模拟70 m光斑内林分地面三维信息。

1.3.3激光雷达辐射传输模型

激光脉冲在大气中传输过程中，受到大气分子、水蒸气和气溶胶等因素影响，产生一系列物理反应(如大气折射、后向散射和大气分子吸收等)，这些反应使得激光脉冲能量值产生一定程度衰减。其中，激光在大气中的透过率决定了激光脉冲回波波形的振幅值，该投过率可由朗伯-比尔定律[31]表示，公式为

(4)

式中τatm(λ)——波长为λ时激光脉冲在大气中的单程投过率

L——激光脉冲发射器与目标物之间的距离

β(λ)——波长为λ时总衰减系数

其中β(λ)由4部分构成[32]

β(λ)=σm+σa+sm+sa

(5)

式中，σm、σa分别为分子、气溶胶吸收系数，sm、sa分别为分子、汽溶胶散射系数。另外，大气总衰减系数值一般为0.5[26]。

1.4 星载激光雷达回波波形信噪比

星载激光雷达实测回波波形中存在大量噪声数据。为高精度模拟仿真国产星载激光雷达回波波形，需对国产星载激光雷达理论仿真波形进行添加噪声处理。本文对验证数据回波波形进行信噪比分析，以确定国产星载激光雷达回波波形添加混合噪声信噪比设定值。其中，星载激光雷达回波波形信噪比定义公式[33]为

(6)

式中RSNR——信噪比

s(i)——去噪前验证数据回波波形

f(i)——去噪后验证数据回波波形

1.5 混合噪声模拟

星载激光雷达回波波形为大尺度遥感数据[34]，激光脉冲在大气传输过程中经由两次菲涅耳衍射后，得到的回波波形为多模式复杂曲线，其中混入了探测器噪声、背景噪声、量子噪声和热噪声等若干高斯分量,这些噪声数据中除背景噪声为非零均值高斯白噪声外，其余噪声均认为是零均值高斯白噪声[35-37]。

由于星载激光雷达回波波形中大部分噪声数据可视为零均值高斯白噪声，因此，本文加入高斯白噪声以代替星载激光雷达回波波形中的噪声数据。另外，加入的噪声数据可通过改变信噪比的方式进行控制。混合噪声模拟公式为

(7)

其中

DRN=DRawRN-mean(DRawRN)

(8)

式中Dnoise——混合噪声模拟数据

DSNRV——信噪比设定值

DRawRN——随机波形数据，波形长度为544

W——国产星载激光雷达理论仿真波形数据

n——国产星载激光雷达回波波形帧数，取544

std()——标准差函数

mean()——求平均值函数

2 实验与结果分析

按照统计学标准，样本抽样个数一般不低于30。本文对仿真波形进行相关性分析时，从样本总个数较多(超过30)的样本中随机抽取30个样本数据作为调查数据，进行相关性分析；对样本总个数较少(低于30)的样本整体进行相关性分析。

2.1 发射波波形

基于验证数据对应测高系统载荷参数，结合发射波函数模型[30]，对其发射波波形以及发射波激光脉冲能量分布进行模拟仿真，仿真结果如图3所示。

图3 发射波仿真波形与激光脉冲能量分布示意图Fig.3 Schematics of emission wave simulation waveform and laser pulse energy distribution

对比分析发射波仿真波形与实测波形，发现发射波实测波形波峰点与发射波仿真波形波峰点之间存在一定距离偏移，如图4a所示。据激光雷达测高原理与发射波函数可知，回波波形测距结果不会因发射波波峰位置变化而改变。故本文在对发射波仿真波形进行验证时，可对其波形波峰点进行左右平移，以获取发射波仿真波形与发射波实测波形相关系数最大值，并把该值作为发射波仿真波形与发射波实测波形相关系数。平移后发射波仿真波形与实测波形示意图如图4b所示。

图4 波峰点平移前后发射波仿真波形与实测波波形示意图Fig.4 Schematics of transmitted wave simulated waveforms and measured waveforms before and after crest point was shifted

平移前发射波仿真波形与发射波实测波形相关系数为0.32。平移后发射波仿真波形与发射波实测波形相关系数为0.96。研究结果表明发射波仿真波形波峰点平移对相关性分析结果会产生较大影响，而不会影响测距结果。因此，有必要对发射波仿真波形波峰点进行平移，以获取发射波仿真波形与发射波实测波形相关系数最大值，并把该值作为两者相关系数实际值。

研究中随机选取30组发射波实测波形，与发射波仿真波形进行相关性分析，得出发射波仿真波形波峰点平移之后与发射波实测波形相关系数均值为0.96。研究结果显示发射波仿真波形与发射波实测波形具有较高的相关性，表明本文建立的发射波仿真模型可对国产星载激光雷达发射波波形进行模拟仿真。

2.2 光斑内地面信息三维模型

为高精度模拟仿真国产星载激光雷达回波波形，需对林木和林地地形高精度模拟仿真。本文依据林木冠层结构把针叶树冠型定义为椭球体型，阔叶树冠型定义为椭球体型、上半球体型和下半球体型。本文在李松等[25]研究基础上，依据实测林地内地形无规律分布特点，通过有限元原理，获取了每个小区域内地物三维坐标，并对其进行曲面拟合，建立了随机地形，如图5a所示。该地形无规律起伏的特点，增大了模拟地形的地表粗糙度，提高了实际地形仿真精度。之后，联合林木冠型定义方式和随机地形，对光斑内林分进行了三维建模，建模结果如图5b所示。

图5 随机地形和林分三维模拟示意图Fig.5 Three-dimensional simulation diagrams of random terrain and forest stands

2.3 仿真波形添加混合噪声

结合1.5节混合噪声模拟原理，对国产星载激光雷达仿真波形进行添加噪声处理：首先，基于混合噪声模拟公式(7)获取噪声波形，噪声波形获取过程中可通过改变信噪比确定添加混合噪声数据的大小；然后，把混合噪声波形数据与国产星载激光雷达仿真波形数据对应帧数振幅相加，可得带有噪声数据的国产星载激光雷达仿真波形，如图6a所示。

图6 国产星载激光雷达仿真波形添加噪声示意图Fig.6 Schematics of adding noise to simulation waveform of domestic spaceborne LiDAR

2.4 仿真波形验证

2.4.1不同地形坡度回波波形

地形坡度较大区域，星载激光雷达回波波形存在波形展宽现象[38-39]，影响森林结构参数估测精度[40]。为探究本文所建回波仿真系统在不同地形坡度条件下是否有效，以10°为间隔把地形坡度分为0°～10°、10°～20°、20°～30°和30°以上共4个小组，对不同地形坡度条件下回波仿真波形与实测波形进行相关性分析，结果如表3所示。

由表3可知，回波仿真波形与实测波形相关系数整体随地形坡度增大而降低。分析其原因：地形坡度较大时，星载激光雷达回波波形出现波形展宽现象[38-39]，且据实地调查发现地形坡度较大时光斑内地形呈无规律起伏分布，以致回波仿真波形与实测波形存在较大差异，降低了仿真波形与实测波形相关性。

表3 不同地形坡度回波波形仿真精度Tab.3 Simulation accuracy of received waveform under different terrain slopes

2.4.2不同郁闭度回波波形

林分郁闭度为林分在地面的投影面积与林地面积之比，对森林每年碳增量存在较大影响[41]。研究中以0.2为间隔，将郁闭度分为0～0.2、0.2～0.4、0.4～0.6、0.6～0.8、0.8～1.0等5个范围，对不同郁闭度条件下回波仿真波形与实测波形进行相关性分析，结果如表4所示。

表4 不同郁闭度回波波形仿真精度Tab.4 Simulation accuracy of received waveform under different canopy covers

由表4可知，仿真波形与实测波形相关系数随郁闭度增大呈先减小后增大的趋势。分析其原因：对郁闭度范围为0～0.2和0.2～0.4条件下仿真波形与实测波形相关系数进行分析，发现郁闭度较小条件下(低于0.2)的林分，其冠层反映在回波波形上为回波波形振幅偏低，且冠层回波能量之和较小。因此，仿真波形中冠层回波波形与实测波形中冠层回波波形振幅值相差较小，且郁闭度在0～0.2条件下地形坡度平均值为9.24°，低于郁闭度0.2～0.4条件下地形坡度平均值10.50°，以至郁闭度在0～0.2条件下仿真波形与实测波形相关系数高于郁闭度在0.2～0.4条件下仿真波形与实测波形相关系数。

针对仿真波形与实测波形相关系数随郁闭度(郁闭度大于0.2)增大而增大的问题，分析其原因：依据回波仿真原理[30]，林分郁闭度可直接影响星载激光雷达回波振幅，而郁闭度越大表明林分中林木个数越多。由于林木在林分中随机分布，光斑内林木个数越多，地面信息三维仿真模型与实际地面信息越相近，极端条件下如郁闭度为1时，林木随机分布这一因素对回波波形仿真精度的影响达到最小。因此，郁闭度大于0.2后，仿真波形与实测波形相关系数随郁闭度增大而增大。

2.4.3不同森林类型回波波形

不同森林类型条件下星载激光雷达回波波形森林结构参数估测精度存在较大差异[42]。为验证本文所建回波仿真系统在不同森林类型条件下是否有效，研究中分别对针叶林、阔叶林和混交林回波仿真波形与实测波形进行相关性分析，结果如表5所示。

表5 不同森林类型回波波形仿真精度Tab.5 Simulation accuracy of received waveform under different forest types

由表5可知不同森林类型条件下仿真波形与实测波形相关系数整体从小到大依次为针叶林、阔叶林和混交林。

依据林木冠形特征和地形分布规律，对不同森林类型条件下仿真波形和实测波形相关系数差异性进行分析：发现针叶树冠型多为圆锥体型，冠层枝叶多分布于冠层中下部位置，反映在回波波形上表现为针叶林冠层回波波形与地面回波波形之间距离较近。当地形存在一定坡度时，受波形展宽影响，冠层回波波形与地面回波波形存在一定程度的重叠；且受林木随机分布影响，地形坡度越大，冠层回波波形与地面回波波形重叠度越大。因此，针叶林仿真波形和实测波形相关系数较低。相对于针叶林，阔叶树冠型多为椭球体型、上半球体型和下半球体型，冠层枝叶多集中于冠层中上部位置，反映在回波波形上表现为阔叶林冠层回波波形与地面回波波形之间距离较远。一定地形坡度条件下，冠层回波波形与地面回波波形重叠度较低。因此阔叶林仿真波形和实测波形相关系数高于针叶林仿真波形和实测波形相关系数。相对于针叶林和阔叶林，混交林同时存在针叶树和阔叶树，反映在回波波形上表现为冠层回波波形与地面回波波形既存在重叠部分，也存在未重叠部分。重叠部分仿真波形多为针叶树冠层回波，该波形与光斑内针叶树冠层实测波形相关性较高(任意地形坡度下)；未重叠部分仿真波形多为阔叶树冠层回波，该波形与光斑内阔叶树实测回波波形相关性较高(任意地形坡度下)。因此，任意地形坡度下，仿真波形中冠层回波与地面回波重叠部分和未重叠部分都与实测波形具有较高相关性，表现在仿真波形与实测波形相关性上为混交林仿真波形与实测波形具有更高的相关性。

2.5 国产星载激光雷达回波仿真波形

对251组验证数据回波波形进行信噪比分析，发现验证数据回波波形信噪比均值为19。因此，本文对国产星载激光雷达理论仿真波形添加噪声时，把混合噪声信噪比设定为19。

2.4节已验证本文建立的回波仿真系统有效性，依据该回波仿真系统，结合国产星载激光雷达载荷参数和实地调查数据，可对国产星载激光雷达回波波形进行模拟仿真。星载激光雷达回波波形受地形坡度影响较大[38-39]，故本文以地形坡度为例，分别对坡度0°、10°、20°和30°条件下国产星载激光雷达回波波形进行模拟仿真，仿真结果如图7所示。

图7 不同地形坡度仿真波形示意图Fig.7 Schematics of simulation waveforms under different terrain slopes

对图7中不同地形坡度条件下国产星载激光雷达仿真波形进行分析，发现国产星载激光雷达仿真波形中地面回波波形长度随地形坡度增大而增大。这与星载激光雷达回波波形受地形坡度影响结果相同[38-39]，说明本文所建回波仿真系统适用于国产星载激光雷达回波波形模拟仿真。

3 结论

(1)验证数据发射波仿真波形与实测波形相关系数均值为0.96，表明本文所得发射波波形可用于星载激光雷达回波波形模拟仿真。

(2)利用有限元原理可实现林分三维信息的模拟仿真，且构建的随机地形与实际地形起伏规律更为相近，可用于国产星载激光雷达回波波形模拟仿真。

(3)不同地形坡度、郁闭度和森林类型条件下验证数据仿真波形与实测波形相关系数均值分别为0.87、0.85和0.87；且国产星载激光雷达仿真波形中地面回波波形长度随地形坡度增大而增大，与波形展宽现象一致。表明本文所建回波仿真系统可用于星载激光雷达回波波形模拟仿真。