一种卫星激光测高仪测距精度仿真计算方法

贺　涛　张新伟（北京空间飞行器总体设计部，北京 100094）

一种卫星激光测高仪测距精度仿真计算方法

贺涛张新伟
（北京空间飞行器总体设计部，北京 100094）

针对卫星激光测高仪在测绘应用中的测距精度计算问题，提出一种测距精度全链路仿真计算方法。方法包括回波仿真、回波处理与精度计算，模拟了测距结果从获取到使用的整个过程，考虑了包括回波处理算法在内的所有影响测距精度的环节。在经典回波模型的基础上，增加了对数字表面模型（DSM）采样间隔的限制条件，以避免选用采样间隔过大的DSM进行仿真导致仿真准确度下降的问题。利用对可见光遥感图像信息进行人工多次选点并取均值的方法，解决真实地形中测量点与卫星之间距离真值的确定问题。将仿真得到的回波数据及精度计算结果与“冰卫星”（ICESat）回波数据及实际测距精度进行比对，结果表明：仿真得到的回波波峰位置和测距随机误差与“冰卫星”的一致，从而验证了仿真计算方法的正确性。

卫星激光测高仪；测距精度；全波形；复杂地形；人工选点

1　引言

测绘卫星的几何分辨率和水平测绘精度在不断提高，但高程精度仍很低，如世界观测-2（World-View-2）卫星的高程精度为5m，我国资源三号卫星无控制点高程精度约为10m［1］。研究表明，激光复合测绘方法利用卫星激光测高仪获取的具有一定精度的激光测距数据参与到立体平差后，测绘卫星的高程精度有明显的提高，从而实现1∶10 000比例尺测绘［2］。因为立体平差算法要求激光测距值达到一定精度，因此，在使用测距数据之前，要根据仪器配置参数及使用方式对测距精度进行仿真计算，预估出测距精度。

目前，针对卫星激光测高仪测距精度的研究包括理论分析方法［3］和仿真计算方法。理论分析方法对地表的描述只有平均高程、倾斜度和粗糙度3个量，在地形复杂的情况下，这3个量无法充分描述地形特点，因此难以实现对城市楼群等复杂地形情况下的测距精度进行分析。仿真计算方法可以通过对整个测距过程的模拟，将复杂地形的影响通过回波仿真体现在回波波形中，不需要地形的描述参量，也就避免了对地形特点描述不充分的问题。文献［4］中提出针对模拟生成的简单地形（平面、斜坡、台阶）条件下的卫星激光测高仪测距精度仿真计算方法，但此方法不适于对激光复合测绘中的卫星激光测高仪在真实地形条件下的测距精度仿真计算。

本文提出一种仿真计算方法，包括回波仿真、回波处理与误差计算3个步骤。通过对不同数字表面模型（DSM）采样间隔下仿真回波的相关性进行分析，评价使用的DSM对仿真准确度的影响，从而对DSM进行筛选，保证仿真的可信度；考虑到回波的信噪比较低且模拟的发射波形为理想高斯波，采用高斯分解法实现波形分解，从而得到被剥离波峰的均值实现重心确定，使仿真计算结果中包含回波处理引入的误差；采用由地面可见光遥感图像信息进行人工选点的真值确定策略，使真值确定策略能在地面高程分布状态未知的情况下使用。因此，本文提出的方法充分考虑了影响测距精度的环节及相关因素，可应用于卫星激光测高仪在测绘复杂地形时的测距精度计算。

2　回波仿真

卫星激光测高仪典型的工作过程为:激光测高仪对待测区域发射一个激光脉冲，经过大气传输后一小部分激光被目标反射回测高仪。光电探测器件将探测到的激光回波转变为电信号，经采样得到激光脉冲回波的数字信号，通过测量激光脉冲飞行时间可计算出距离值［4］。

根据卫星激光测高仪的工作过程，将激光发射到回波采样的各个环节进行分解，将回波仿真分为激光发射脉冲仿真、地表反射仿真、激光能量传输仿真、背景光噪声仿真、光电转换仿真和电路噪声仿真6个部分，并分别建立相应的数学模型，仿真流程如图1所示。

图1　激光测距回波仿真流程Fig.1　Simulation process of laser ranging record waveform

首先，对仿真区域的DSM进行三角格网化处理［5］，结合激光指向角得到发射激光束中心轴线与每片三角网格平面的夹角后，设置卫星激光测高仪性能参数并仿真得到激光发射脉冲，读取三角格网平面与激光束夹角；然后，根据回波全波形采样频率及卫星高度，设置仿真时间步长及仿真时间段，结合DSM数据与发射激光数据，得到激光光斑内部地面上每个三角格网内反射的光能量；根据地反射模型，结合滤光片的带宽及激光回波的能量，考虑接收望远镜的接收口径，将地面所有三角格网在某一时刻的反射光相加，得到该时刻到达光电探测器的光功率；最后，根据光电探测器件的光电转换、滤波特性及电路噪声特性，按照设定的量化位数，得到量化后的回波数据。

主要仿真步骤的数学模型如下。

1）激光发射脉冲模型

激光发射脉冲一般采用理想高斯模型，其发射激光强度如式（1）所示［6］。

式中:r为激光束截面上点与发射光轴的距离；L为激光束传播的距离，对于激光测距的目标，L= Ltarget，Ltarget为目标到激光测高仪的距离；t为激光脉冲传播的时间；θtr为激光发射脉冲光斑中心能量e—1/2处的光束发散角；Etr为激光发射脉冲的能量；ktr为激光发射脉冲的均方根脉宽。

2）地表反射模型

地表高程信息由DSM给出，采用三角格网将DSM中离散的点连成面，计算每个三角格网对入射其上的激光束的反射能量。由于地表反射特性数据难以获取且易发生变化，为简化计算，仿真中将地表按各点反射率相同的朗伯体处理。在激光收发同轴的测距过程中，其反射光束与入射光束满足式（2）［7］。

Iref=Itr·cosθdem·ηsurface（2）

式中:Iref为反射激光强度；Itr为入射激光强度，即激光发射脉冲强度；θdem为入射光束与三角格网面的夹角；ηsurface为激光垂直地面入射时的地表反射率，仿真中各点取值相同。

3）激光能量传输模型

激光在大气中传输会产生延迟，若激光斜入射则有光束弯曲，其根本原因是大气折射率不为1，存在空间差异，并且随大气条件变化而变化，因此要根据气象参数对测距结果进行校正［8］。目前，常用的校正方法是天顶角延迟与映射函数的乘积模型［9］，较为复杂。为提高仿真效率，本文在仿真中不计算产生的延迟，在处理时也不进行校正，只加入校正的残差（厘米量级）作为测距的随机误差［9-10］，可获得与先仿真计算延迟量再校正一样的效果。不考虑大气影响，激光回波接收望远镜接收到的光强Irec满足菲涅尔衍射传播规律，如式（3）所示。

式中:α为能量衰减系数；c为真空中的光速。

4）背景光噪声模型［11］

接收望远镜可接收到处于滤光片带宽范围内的光能量，其中一部分为地表反射的太阳光，即为背景光噪声，其功率见式（4）。

式中:Eλ为窄带滤光片中心波长λ附近波长的太阳光在目标表面产生的光谱辐照度；Δλ为窄带滤光片带宽；ηre为接收望远镜光学系统光学效率；θre为接收视场角；ε为窄带滤光片中心波长附近的地表反射率；Dre为接收望远镜的接收口径。

5）光电转换模型

雪崩光电二极管（APD）在线性工作区内，其接收光功率Pr（t）与输出电压Vd（t）之间的关系如式（5）所示。

Vd（t）=γGw Pr（t）（5）

式中:γ为APD中电子和空穴电离概率之比；G为APD的增益系数；w为电子电量。

接收光功率Pr（t）是背景光噪声功率PB与激光反射光功率之和，如式（6）所示。

式中:（x，y）为地面三角格网中心点的坐标；r（x，y）为点（x，y）与激光发射光轴的距离；L（x，y）为点（x，y）与激光测高仪的距离；S为地面每个三角格网的面积；代表求激光光斑内部所有点反射光功率之和。

6）电路噪声模型［11］

电路噪声包含光电探测器噪声、光子噪声及热噪声，这3类噪声均可认为是均值为0的高斯噪声，其标准差代表噪声的强度。

卫星激光测高仪一般采用APD作为光电探测器，其噪声主要包括暗电流噪声id和前置放大电流噪声ia，产生的噪声电压up方差如式（7）所示。

式中:σ2（id）为暗电流噪声id的方差；σ2（ia）为前置放大电流噪声ia的方差；Pd和Pa分别为暗电流噪声和前置放大电流噪声的功率谱密度，A2/Hz；RL为APD等效电阻；B为等效噪声带宽。

光子噪声产生的电压方差如式（8）所示。

式中:keff为空穴电离系数与价带电子电离系数之比；hf为电路的冲击响应；ηd为衰减系数；h为普朗克常数；υ0为光频率；Pg为入射至APD表面的光功率；τ为参与积分的时间变量。

热噪声产生的电压方差如式（9）所示。

式中:K为玻尔兹曼常数；Tdet为探测器温度。

为保证仿真结果的准确性，须要对仿真中使用的DSM进行评价，分析其采样间隔对仿真准确度的影响程度。对不同采样间隔的DSM分别进行回波仿真，通过多点仿真计算各采样间隔条件下波形相关系数的均值，给出回波相关性与DSM采样间隔的关系如图2所示。目前，可获取最小采样间隔的真实地形DSM为德国法伊欣根城市楼群的DSM，其间隔为0.09m，故以此为参考波形，每次将DSM进行降采样，从中可得到相关系数随DSM采样间隔的变化趋势。进行降采样试验后，得到如图2中的曲线，其中相邻两点的相关系数差值小于0.05时，两者中采样间隔较小的DSM仿真波形与无穷小采样间隔仿真回波的相关系数大于0.9，对测距精度计算结果影响普遍小于3%。

图2　相关系数与采样间隔的关系Fig.2　Relation between relevance and sampling distance

3　回波处理与误差计算

激光测距的全波形回波处理算法已十分成熟，常用的有高斯分解法、小波变换法和反卷积法［12］。在实际应用中，一般会根据具体情况选择理论上引入误差最小的算法。发射波为理想高斯波时，高斯分解法与反卷积法均能得到较好的测距精度，但高斯分解算法中的噪声对测距精度的影响更小。本文中，由于回波仿真中的发射激光为理想高斯波，且回波信号信噪比较低，因此选用理论上引入误差最小的高斯分解法作为回波处理算法，流程如图3所示。

图3　高斯分解法的流程Fig.3　Process of Gauss-decomposition method

测距精度为测量值与真实值之差。通过上述仿真，已可以得到回波波形中每个波峰对应的测距值，但由于激光光斑覆盖范围较大，从回波数据中无法得到波峰与地面点的对应关系，即无法得到真实值。为解决以上问题，采取人工选点的做法，下面以一个测量点的遥感影像及仿真回波为例，进行说明。

（1）计算出光斑中心位置，得到光斑覆盖区域的遥感影像，如图4所示。

（2）在仿真得到的回波波形（如图5所示）中，可找到时间轴上的最后一个峰值。

（3）考虑到地面或屋顶存在一定倾斜，就算找到了该面与波峰的对应，也只知道该波峰代表的距离值对应的是卫星与该面间的平均距离。为克服随机性，采用多次选点的方式。选点时，选点人利用所有经验，尽可能使选取点与卫星距离与测量值一致。取这些点真实高程的均值作为高程真值，从而得到与卫星距离真值，与测量值做差，即为该点测距误差。

图4　光斑位置示意Fig.4　Position of laser footprint

图5　激光测距回波仿真结果Fig.5　Simulation result of laser ranging record waveform

4　实例验证

下面给出仿真计算ICESat测距精度的实例，用于说明本文方法的用途并验证方法的正确性。按照ICESat的激光测高仪设置仿真计算参数［11，13］，如表1所示。

基于 MATLAB软件编写回波仿真程序，选取ICESat在2005年2月24日和2008年10月9日对加拿大多伦多城市区域进行测量的位置［14］。对该区域DSM进行评价，判断其采样间隔是否符合要求。目前可获取的为国际摄影测量与遥感学会（ISPRS）网站提供的采样间隔为25cm的DSM［15］，对该采样间隔进行多次加倍后，根据ICESat卫星的参数配置进行回波仿真，得到采样间隔与波形之间的相关系数，见表2。从表2中可知，0.25m与0.5m采样间隔DSM仿真所得波形的相关系数差值均小于0.05，可判断0.25m采样间隔DSM用于回波仿真时，仿真回波与采样间隔无穷小情况下的仿真回波的相关系数大于0.9，对测距精度计算结果影响小于3%。

表1　ICEEat的激光测高仪测距精度仿真计算参数Table1　Eimulative computation parameters of ranging accuracy for ICEEat laser altimeter

表2　多伦多城市区域DEM采样间隔与对应回波的相关系数Table2　Relevance index and DEM sampling distance in Toronto

利用该DSM进行回波仿真，并与ICESat实测回波进行比对，结果如图6所示。可见，ICESat实测回波与仿真回波的主要波峰位置均相同，证明了回波仿真的正确性。值得注意的是，仿真回波与实测回波存在一定的差异，可能的原因包括DSM与真实地形差异、反射率因素、激光光斑位置误差、发散角外激光反射。

每个测量点的回波经处理后得到距离测量值，并换算为测量点的高程值；再对激光地面光斑内区域进行人工选点，得到高程真值。在ICESat数据产品中，GLA06级产品包括经过校正的数字高程模型（DEM）数据，就是回波的最后一个峰对应的高程。数据产品的高程值中包含测距误差、姿态轨道误差及坐标系误差。由于卫星的天顶角仅有0.3º，可以认为测距精度与高程测量精度是一致的，因此统一写为测距精度（见表3）。定轨误差为3cm，姿态轨道造成的误差较小且为低频变化量，短时间内表现为系统误差。坐标系带来的误差也为系统误差。因此，高程中包含的随机误差与测距的随机误差是相同的。将仿真得到的测量结果与DSM点的高程值进行比较，可得到测距精度仿真计算结果；将ICESat的产品与DSM点的高程值进行比较，可以得到ICESat测距精度。具体比较结果如表3所示。

图6　仿真回波与ICESat实测回波对比Fig.6　Comparison between simulation record waveform and real record waveform by ICESat

表3　仿真得到高程值、ICEEat产品中给出高程值与DEM图中高程值的比较Table3　Elevation comparison of simulation result，ICEEat data and DEM

对表3中的仿真计算测距精度和ICESat的测距精度进行比较，可以验证本文仿真计算方法的正确性。可以看出:仿真计算测距精度均值十分接近0；仿真计算测距精度的标准差与ICESat数据测距精度的标准差基本一致。由于没有足够的数据，目前仅能用这两组数据对仿真计算方法的正确性进行初步验证。从验证结果来看，仿真计算方法得到的精度标准差与ICESat真实测距精度标准差的差异均小于0.1m，初步证明了仿真计算方法是正确的。

5　结束语

本文针对卫星激光测高仪的测距精度，提出了一种从测距数据获取到应用的仿真计算方法，通过回波仿真、回波数据处理、人工选点评价测距精度的方式，使此方法适用于测绘应用中的真实复杂地形。回波仿真结果及测距精度的仿真计算结果与ICESat实际的回波及测距精度的比对表明:使用本文方法进行激光测高仪测距回波的仿真及随机误差的计算，具有一定的准确度。利用该方法可以进行激光测高仪测距精度的仿真计算，以预先剔除低精度测距数据，避免其参与到测绘平差中；也可以通过改变某一因素对应的仿真参数值，研究各因素对测距精度的影响规律，为激光测高仪的设计提供一定的参考。

（References）

［1］李德仁，王密.资源三号卫星在轨几何定标及精度评估［J］.航天返回与遥感，2012，33（3）:1-6 Li Deren，Wang Mi.On-orbit geometric calibration and accuracy assessment of ZY-3［J］.Spacecraft Recovery& Remote Sensing，2012，33（3）:1-6（in Chinese）

［2］梁德印.应用激光测高仪提高测绘卫星定位精度的研究［J］.航天器工程，2014，23（6）:17-23 Liang Deyin.Study on improving positioning accuracy of mapping satellite using laser altimeter［J］.Spacecraft Engineering，2014，23（6）:17-23（in Chinese）

［3］马跃，阳凡林，卢秀山，等.对地观测星载激光测高系统高程误差分析［J］.红外与激光工程，2015，44（3）:1043-1047 Ma Yue，Yang Fanlin，Lu Xiushan，et al.Elevation error analysis of spaceborne laser altimeter for earth observation［J］.Infrared and Laser Engineering，2015，44（3）: 1043-1047（in Chinese）

［4］朱近，孙世君.星载激光测高仪回波信号仿真分析［J］.航天返回与遥感，2013，34（1）:67-78 Zhu Jin，Sun Shijun.Simulation analysis of space borne laser altimeter echo signal［J］.Spacecraft Recovery& Remote Sensing，2013，34（1）:67-78（in Chinese）

［5］潘浩，李国元，王华斌，等.一般地形大光斑激光测高回波模拟研究［J］.地理信息世界，2015（2）:82-87 Pan Hao，Li Guoyuan，Wang Huabin，et al.Research on large spot laser altimetry echo simulation of general terrain［J］.Geomatics World，2015（2）:82-87（in Chinese）

［6］Bahaa E A Saleh，Malvin Carl Teich.Fundamentals of photonics［M］.2nd ed.New Jersey:John Wiley& Sons，2007

［7］张雷洪，孙刘杰，马秀华.激光跟踪中目标卫星表面BRDF对回波信号的影响［J］.红外与激光工程，2012，41（8）:2049-2052 Zhang Leihong，Sun Liujie，Ma Xiuhua.Influence of bidirectional reflection distribution function of target satellite surface on echo in laser tracking［J］.Infrared and Laser Engineering，2012，41（8）:2049-2052（in Chinese）

［8］赵欣，张毅，赵平建，等.星载激光测高仪大气传输延迟对测距精度的影响［J］.红外与激光工程，2011，40（3）: 438-442 Zhao Xin，Zhang Yi，Zhao Pingjian，et al.Influence of atmospheric transmission delay of satellite laser altimeter on ranging precision［J］.Infrared and Laser Engineering，2011，40（3）:438-442（in Chinese）

［9］李松，马跃，周辉，等.星载激光测高系统大气折射延迟的影响分析［J］.应用光学，2014，35（3）:520-524 Li Song，Ma Yue，Zhou Hui，et al.Analysis on atmospheric refraction delay for satellite laser altimeter system ［J］.Journal of Applied Optics，2014，35（3）:520-524（in Chinese）

［10］马跃，李松，翁寅侃，等.星载激光测高仪大气干项延迟校正［J］.红外与激光工程，2013，42（4）:910-914 Ma Yue，Li Song，Weng Yinkan，et al.Hydrostatic delay correction for satellite laser altimeter［J］.Infrared and Laser Engineering，2013，42（4）:910-914（in Chinese）

［11］周辉.激光测量系统地表测量技术指标综合论证报告［R］.武汉:武汉大学，2015 Zhou Hui.The demonstration report on overall qualification of laser ranging system［R］.Wuhan:Wuhan University，2015（in Chinese）

［12］王成，习晓环，骆社周，等.星载激光雷达数据处理与应用［M］.北京:科学出版社，2015 Wang Cheng，Xi Xiaohuan，Luo Shezhou，et al.Data processing and application of satellite lidar［M］.Beijing:Science Press，2015（in Chinese）

［13］Peggy L J.The ICESat/GLAS instrument operations report，NASA/TM-2012-208641［R］.Washington D. C.:NASA，2012

［14］National Snow and Ice Data Center.GLAS/ICESat L1A global altimetry data，Version 33［EB/OL］. ［2016-01-27］.http://nsidc.org/data/GLA01

［15］International Society for Photogrammetry and Remote Sensing.WG III/4:3D scene analysis［EB/OL］. ［2016-01-27］.http://www2.isprs.org/commissions/ comm3/wg4.html

（编辑：夏光）

Simulative Computation Method of Ranging Accuracy for Satellite Laser Altimeter

HE Tao ZHANG Xinwei
（Beijing Institute of Spacecraft System Engineering，Beijing 100094，China）

To consider the ranging accuracy computation of satellite laser altimeter applied in mapping，a complete chain simulative computation method is established.The method includes record waveform simulation，record waveform processing and ranging accuracy computation，and simulates the over-all working stage of the altimeter.Based on the existing model，a constraint of DSM（digital surface model）sampling interval is adopted to avoid any possible simulation accuracy degradation.The true distance is decided by artificial point choice.After it is done several times，the mean value is obtained，which solves the problem of determining true value. The comparison between ICESat data and simulation result is performed.Compared with the waveform and true ranging accuracy of ICESat，the simulation result shows that the location of record waveform peak and the ranging random error of the simulative computation and ICESat data match well.

satellite laser altimeter；ranging accuracy；full waveform；complex terrain；artificial point choice

V447.1

A

10.3969/j.issn.1673-8748.2016.03.015

2016-02-01；

2016-03-03

贺涛，男，硕士研究生，研究方向为航天器总体设计。Email：zd09ht@163.com。