全波形星载激光测距误差抑制的滑动窗口高斯拟合算法

谢俊峰,刘 仁

1. 自然资源部国土卫星遥感应用中心, 北京 100048； 2. 河海大学地球科学与工程学院, 江苏 南京 211100

星载激光测高仪作为主动式测量设备，因其较高的测距精度，已广泛被用于深空探测和对地观测[1]。在对地观测方面，2003年1月美国首次发射了GLAS(geoscience laser altimeter system)激光测高系统[2-3]，而后我国于2016年5月发射资源三号02星，该星搭载了一台激光测高仪作为试验性载荷用于对地观测[4-6]。随着我国激光技术的发展，我国于2019年11月发射了高分七号(GF-7)星载激光测高系统[7-8]，2020年7月发射了资源三号03星激光测高仪[9]，同年12月我国发射的高分十四号卫星(天绘三号)也搭载了一套激光测高系统，2021年还将发射陆地生态碳卫星激光测高系统。

星载激光测高仪在轨后，研究人员主要关注于修正卫星在轨后激光指向角及测距系统误差[1,4,10-14]，却鲜有对由激光时间测量引起的测距随机误差进行改正的研究。基于阈值时刻鉴别体制的星载激光测高仪，一般利用前后缘阈值时刻估算波形重心位置作为激光出光与返回时刻，进行激光渡越时间测量[15-17]，计算激光测距值。这种模式下测距随机误差难以被发现，测距精度相对较低。对于全波形星载激光测高仪(如GLAS，GF-7、陆地生态碳卫星、天绘三号卫星)，其测距方式通过发射与返回波形峰值时间差计算得到[18]，其测时原理如图1所示。由于激光模拟回波信号经数字化处理后的量化误差存在，导致基于波形获取的时间存在一个微小的随机误差，这将产生厘米级的测距随机误差。而该测距误差，对于如高分七号搭载的厘米级别测距精度的全波形激光测高仪而言，往往不应被忽略。

图1 全波形星载激光测高仪测时原理示意Fig.1 Schematic of time measurement of full-waveform spaceborne laser

为了减小全波形星载激光测高仪测距随机误差，本文针对其发射与返回波形，提出了一种依赖峰值初始位置的全波形星载激光测距误差抑制的滑动窗口高斯拟合算法。该方法在确定峰值位置后，应用直线判断原理，抑制了波形中噪声，随后利用有效波形点开展波形拟合，确定最优峰值时刻，从而提取精确的激光测距值。由于GF-7卫星在轨运行期间激光指向整体精度优于1.5″，使得由激光指向抖动引起激光在平坦地形的测距误差可以忽略不计。故选取平静的内陆湖面、冰面及江苏平地为试验区域，以GF-7星载激光数据为试验对象，利用一般峰值法与本文方法提取的测距值，通过计算湖面与冰面的激光相对高程精度，以及利用高精度LiDAR[19]计算平地的绝对高程精度，从而验证了本文方法对测距随机误差的改进程度。

1 全波形星载激光测距提取算法

1.1 基于一般峰值法提取的测距误差分析

通常情况下，波形采样间隔非常小，大多数大气和地形条件下激光获取的回波近乎平滑。一般全波形激光发射和返回波形的最大值对应时刻即为激光发射或者返回时刻，根据它们的时间差计算渡越时间，称为一般峰值法。然而，实际在轨工作中，星载激光模拟回波信号经数字化处理后产生量化误差，从而导致了波形峰值处离散采样点存在似噪声的波动，使得部分情况下发射与返回波形出现多个最大值，或是峰值点并不是真实的波形峰值点(图2、图3)等情况，最终导致基于一般峰值法计算的激光测距存在明显的随机误差。

图2 GF-7星载激光测高仪3种典型非标准高斯发射波形峰值Fig.2 Three typical non-standard Gaussian peak values of GF-7 spaceborne laser emission waveform

图3 GF-7星载激光测高仪三种典型非标准高斯返回波形峰值Fig.3 Three typical non-standard Gaussian peak graphs of GF-7 spaceborne laser echo

如图3所示，一般峰值法依据发射与返回波形提取的激光渡越时间将产生1～5个采样间隔的随机误差，对于目前较高的波形0.5 ns采样频率(2 GHz)的GF-7星载激光会导致0.05～2.5 ns时间测量误差，即0.075～0.375 m的测距误差。对于其他常用的低采样频率星载激光测高仪来说，一般峰值法带来的随机测距误差将翻倍。综上，对于全波形星载激光测高仪，一般峰值法提取的激光渡越时间直接影响了激光的测距精度，继而影响星载激光测高精度。

1.2 基于滑动窗口的高斯曲线拟合峰值的测距误差抑制

1.2.1 波形噪声抑制与激光渡越时间提取

针对全波形激光模拟回波信号经数字化处理后产生量化误差，引起的激光测距随机误差，本文提出了一种全波形星载激光测距误差抑制的滑动窗口高斯拟合方法，适用于星载激光各种情况下的发射与返回波形峰值拟合，方法基本流程如图 4所示。

图4 基于滑动窗口的高斯曲线拟合方法流程Fig.4 Flow chart of Gaussian curve fitting method based on sliding window

依据上述流程，本文方法核心步骤如下：

(1) 波形初始峰值检索。

根据输入波形查找到波形最大值作为初始峰值，以它为中心将波形分为左侧(上升沿)与右侧(下降沿)波形进行算法实现。

(2) 基于滑动窗口波峰左侧(右侧)波形有效点提取。

算法设置1×3窗口，由初始峰值点沿波形下降方向(作为正向)进行滑动检索，沿正向方向单个窗口内第1个点定义为P1，第2、3个点分别定义为P2、P3。若P2=P1，剔除P1，窗口起点滑动至P2重新获取3个波形采样点计算，直至P2≠P1。此时利用窗口内前两个点P1(x1,y1)和P2(x2,y2)构建的正向直线L为

L12(x,y)=(x2-x1)(y-y1)-

(y2-y1)(x-x1)

(1)

对于已构建的正向直线，窗口内第3个波形点所在位置可定义为3类。①直线上点；②内点：沿正向直线L前进方向右侧点；③外点：沿正向直线L前进方向左侧点。将窗口内波形第3个点P3(x3,y3)，代入式(1)，进行内外点判断[20]，计算公式为

(2)

对于标准高斯波形，波峰右侧波形沿下降方向窗口内第3个点始终为内点；波峰左侧恰好相反，窗口内第3个点始终为外点，如图5(a)所示。

图5 高斯波形采样点分布与波峰右侧波形噪声剔除Fig.5 Gaussian waveform sampling point and noise eliminating on right waveform

实际在轨星载激光波形，1×3的窗口内第3个点P3仅存在两类情况：①波峰左侧波形P3为外点或者右侧波形P3为内点；②P3为内点(波峰左侧波形)或为外点(波峰右侧波形)。当在第1种情况下，P1、P2、P3满足标准高斯曲线分布，直接保留P1、P2、P3，窗口起点滑动至P2重新计算。当在第2种情况下，P3为波形凸点；其中，若P3y=P1y，P1,P2为噪声点，直接剔除且保留P3，窗口起点滑动至P3再次计算；若P3y≠P1y，P2为噪声点直接剔除，保留P1,P3且窗口起点滑动至P3再次计算(以波峰右侧波形为例，示意图如图5(b)所示)。上述两种情况下，直至波峰左侧和右侧波形同时各保留N个凸点(N≥2)时停止试验，对于GF-7星载激光，本文设置N为6。

(3) 高斯曲线拟合。

星载激光测高仪发射与返回波形均满足标准高斯分布，利用高斯曲线方程对选定的N个波形拟合点进行波形拟合，曲线拟合方程为[21]

(3)

式中，y为波形拟合点幅值；x为波形拟合点时刻；x0为拟合波形峰值时刻；δ为拟合波形脉宽；A为拟合波形幅值。

1.2.2 激光测距值计算

根据激光计时系统记录的发射与返回波形起始时间，加上它们各自到波峰位置的时间差即可计算出激光渡越时间，并转换为激光测距值。以GF-7全波形星载激光测高仪为例，其测距公式为[15]

τR=a*((T2+Techo)-(T1+Twf))+b

(4)

(5)

式中，Rlaser为星载激光精确测距值；τR为激光渡越时间；T1为发射波形起点时间；T2为返回波形起点时间；Twf为发射波形峰值至发射波形起点时间；Techo为返回波形峰值至返回起点时间；a为激光计时系统常数因子；b为激光计时偏移改正量；c为光速，c=299 792 458 m/s。

1.3 基于激光脚点高程的测距精度验证

本文采用控制变量法，在仅改变激光测距值条件下进行激光脚点高程计算，对比在不同测距值下的高程精度。首先，根据星载激光在轨几何定位原理，考虑大气[22-23]、潮汐[24-25]引起的测距误差，构建严密几何定位模型[4]；然后，利用一般峰值法与本文方法提取的测距值分别计算激光点高程；最后，利用平静湖面、冰面及高精度LiDAR数据，进行两类测距值的激光点高程相对精度与绝对精度对比。其流程如图6所示。

图6 基于激光脚点高程的测距精度验证流程Fig.6 The flow of ranging accuracy verification based on laser footprint elevation

(1) 基于激光测高相对精度的测距验证。

借助内陆平静的湖面和冰面，计算激光在湖面或者冰面上各点的高程，随后统计该湖面或冰面上的激光高程标准偏差

(6)

理想条件下，在同一平静的湖面或者冰面激光脚点高程基本一致。由于激光测距随机误差存在，同一平静湖面和冰面激光高程会存在一定波动，可采用激光点高程标准偏差衡量测距精度，即同一湖面或冰面激光高程标差偏差越小，测距精度越高[4]。

(2) 基于激光测高绝对精度的测距验证。

将激光测距值作为影响激光脚点高程唯一变量，选取高精度的地面控制数据，如：RTK控制点，高精度LiDAR点云，验证激光脚点高程绝对精度。通过计算激光点与地面控制点高程差值(绝对精度)的均值hmean与中误差hRMSE，用于对比不同测距值下的激光高程绝对精度，其均值与中误差计算公式为[4]

(7)

(8)

式中，hli为第i激光点高程；hgi为第i激光点地面实际高程；n为试验中的激光点个数。

2 试验与验证

2.1 试验数据

本文以GF-7星载激光测高仪为试验对象，其采用双波束激光同时对地观测，每条波束发射能量均为100～180 mJ，脉宽为4～8 ns，发散角为30 μrad[8]。激光接收系统采用数字化回波采样设备，同时记录下发射与返回高增益和低增益波形，波形采样间隔为0.5 ns，回波最大采样长度为1200个采样间隔[7]。

由于GF-7卫星过顶时间短，激光地面点数据较少，地面密度小，故而在基于激光测高相对精度的测距验证试验中，选取湖泊较大的瑞典维纳恩湖和我国江苏省太湖作为试验区域。维纳恩湖位于瑞典南部，最北边纬度约为北纬59.3°，GF-7星载激光经过该湖泊的数据为2019年12月19日的第691轨，此时维纳恩湖已经进入寒冷冰期，湖面已变成厚厚的冰面。激光在该湖泊分布如图7所示，由于天气因素，波束1点有效点8个，波束2有有效点5个。太湖常年无冰期，风浪较小，适于进行试验，其中GF-7过太湖的激光数据为2020年5月8日第2844轨，如图8中东南角所示，但由于激光测量时地表云层较厚导致波束1部分数据无返回信号，波束2湖面无激光点。故采用波束1位于太湖湖面8个有效激光点(南北跨越50 km)。

图7 维纳恩湖冰面试验区域及其试验数据分布Fig.7 The experimental area and its experimental data map of the ice surface of Lake Vanern

图8 江苏省试验区域与试验数据分布Fig.8 The experimental area and its experimental data map of Jiangsu province

基于激光测高绝对精度的测距验证试验，本文选取GF-7星载激光过江苏省西北部宿迁地区的激光点，如图8中西北角所示，该数据为2019年12月12日GF-7卫星第595轨。该区域范围为：33.152 46°N—34.683 58°N，118.104 62°E—118.619 99°E，其中波束1共55个激光点，波束2共67个激光点。用于进行激光点高程绝对精度验证数据，为采用徕卡ALS70机载激光扫描系统获取区域内的LiDAR点云数据。该设备最大脉冲频率为500 kHz，最大扫描频率为200 Hz，最大视场角为75°，可接收无限次回波记录。现场测量时飞行航高为2400 m，地面光斑大小约为60 cm，获取的点云密度约为1.53/m2，高程精度中误差为0.12 m，数据采集时间为2018年，获取的数据面积约3432 km2。

2.2 基于平静湖面的激光高程相对精度验证

2.2.1 波形噪声抑制试验与分析

采用本文基于滑动窗口高斯拟合方法，分别对GF-7星载激光第2844轨波束1湖面上8个激光点，以及第691轨冰面上波束1的8个激光点和波束2的5个激光点，共计21个激光点的发射与返回波形进行试验。剔除波形中噪声点，利用自动挑选出波峰周围12个波形点进行波形峰值拟合，3类典型的发射与返回波形拟合结果如图9和图10所示。卫星实际在轨后，激光波形峰值形状主要为6种情况，其中发射波形3种，返回波形3种，以下利用本文方法对这6种情况下的波形拟合结果进行详细分析。

图9 发射波形拟合结果Fig.9 Fitting results of emission waveform

图10 返回波形拟合结果Fig.10 Fitting results of echo waveform

针对以上6种不同的情况下的波形，由整体结果可以看出，本文算法均能准确地拟合出回波波形，其中每类波形脉宽吻合度极高，拟合波形峰值与实际回波最大值点幅值相当。从拟合波形细节观察，本文算法避开了各种情况下的原始波形中噪声点；且在上述6种典型波形情况下，拟合后波形峰值均能位于原始波形幅值最大的两个采样点之间。充分说明本文算法能够有效地拟合出实际波形的峰值，抑制激光发射与回波波形中噪声，从而可根据拟合后的波峰即可精确获取激光发射或返回波形时刻。

根据一般峰值法与本文方法，计算第691轨冰面上波束1的8个激光点和波束2的5个激光点的发射脉冲时刻(Twf)与返回时刻(Techo)，结果见表1。两种方法计算的发射和返回波形峰值时间差分别为ΔTwf、ΔTecho；由两种方法计算的测距差值为ΔTrange。

表1 一般峰值法与本文方法提取的激光渡越时间差值Tab.1 The time difference of laser transit time between general peak method and this paper method

由表1结果可分析，无论是星载激光测高仪发射波形还是返回波形，经本文方法拟合后波形峰值时刻发生了一定偏移。其中，相对于一般峰值法，本文方法提取的发射波形峰值点位置偏移了0.235 ns；本文方法提取的返回波形峰值点位置偏移了0.49 ns。将发射波形与返回波形峰值偏移量应用于GF-7星载激光测高仪测距计算式(4)之中，本文方法提取的激光渡越时间与一般峰值法提取的渡越时间偏差为0.50 ns，相当于7.5 cm的测距误差。根据上述分析可知，本文方法能够明显减小激光测距随机误差。

2.2.2 平静湖面激光高程相对精度验证

为分析本文方法提取的激光测距值，能否有效提升星载激光测高仪相对测高精度。根据在轨检校后激光指向角，分别利用一般峰值法与本文方法提取的测距值，对第2844轨波束1太湖湖面上8个激光点，以及第691轨冰面上波束1的8个激光点和波束2的5个激光点进行激光脚点高程解算，并转换至正常高，结果见表2、表3。试验中，仅有激光测距唯一变量，其余参数，如激光指向、大气改正/潮汐改正值、高程异常等完全一致。

根据表2、表3结果可以看出，在维纳恩湖冰面上，GF-7星载激光波束1与波束2高程相对精度基本一致，其中利用一般峰值法提取的激光测距值计算的激光高程标准偏差为0.125 m，基于本文方法提取的激光测距值计算相同激光点高程的标准偏差为0.082 m。对于内陆太湖，利用一般峰值法提取的激光测距值计算太湖湖面8个激光点正常高的标准偏差为0.123 m；而基于本文方法提取的激光测距值计算相同激光点正常高的偏差为0.065 m。其中本文方法对于太湖湖面的激光点精度提升更高，然而实际情况下维纳恩湖冰面结果更加准确、可信，因为太湖湖面或多或少存在一定风浪，产生了该现象。综上，本文方法提取的激光测距值对星载激光测高仪相对高程精度提升了近4.2 cm，由原始的12.5 cm提升至8.3 cm。

表3 GF-7激光点在瑞典维纳恩湖冰面高程相对精度Tab.3 The relative elevation accuracy of GF-7 laser footprint on the ice surface of Vanern Lake in Sweden

表2 GF-7激光点在太湖湖面高程相对精度

假设理想条件下，即激光测距无误差、同一无风浪湖面或冰面上的激光测高应该完全一致。因此，激光在冰面或无风浪湖面的高程相对偏差直接大致可以反映出星载激光测高仪测距精度。结合表3冰面的试验结果分析，由表1中可以看出，相对于一般峰值法，基于本文方法提取的激光测距精度提升了7.5 cm。

2.3 基于高精度机载LiDAR的激光高程绝对精度验证

分别对GF-7第595轨江苏境内的波束1的55个激光点和波束2的67个激光点进行波形拟合，并利用一般峰值法与本文方法提取出两套测距值。在其他参数完全一致的条件下，根据激光几何定位模型计算每个激光脚点大地坐标。从激光足印落点位置，发现GF-7卫星第595轨激光过宿迁市等多个城区，如图8所示，部分星载激光点落在房屋等建筑与林木上。剔除这些异常激光点后，波束1和波束2分别剩余48和61个激光点。利用机载LiDAR点云内插激光点高程，将激光高程与内插的LiDAR的高程差的均值和中误差作为GF-7星载激光高程绝对精度。GF-7星载激光平地地区两类测距计算的高程绝对精度如图11和图12所示。

由图11和图12结果显示，无论是GF-7激光波束1或波束2，本文方法提取的测距值计算的激光高程差曲线向0值附近收缩，可判断出本文方法能够明显提升星载激光高程绝对精度。为定量分析其对激光高程绝对精度提升空间，统计上述波束1与波束2所有点在两类测距值下的均值与中误差，结果见表4。

图11 两类测距值下第595轨波束1激光高程绝对精度Fig.11 The 595th track beam 1 laser elevation absolute precision under two sets of ranging value

图12 两类测距值下第595轨波束2激光高程绝对精度Fig.12 The 595th track beam 2 laser elevation absolute precision under two sets of ranging value

由表4可以看出，在平原地区，基于一般峰值法提取的测距值计算的GF-7星载激光波束1激光高程绝对精度为-0.056±0.192 m；基于本文方法提取的测距值计算的波束1激光高程绝对精度为-0.040±0.177 m。同理，GF-7星载激光波束2在两类测距值下的高程绝对精度分别为-0.038±0.238 m、-0.029±0.189 m。总体来看，本文方法可将GF-7星载激光高程绝对精度由初始-0.047±0.215 m提升至-0.035±0.183 m，提升了4.5 cm。从量级上看，高程绝对精度提升不大，但对于绝对高程精度本已非常高的GF-7星载激光测高仪而言，再提升4.5 cm精度具有十分重要的现实意义。

表4 两类测距值下平地地区激光测高绝对精度

3 结 论

本文在分析星载激光测高仪全波形数据过程中，发现波峰附近存在明显类似噪声的现象，采用一般峰值法提取的测距值存在较大随机误差的问题，故而提出了一种全波形星载激光测距误差抑制的滑动窗口高斯拟合算法，用于提取激光测距值，该方法明显提高了星载激光测距精度。利用一般峰值法与本文方法提取的测距值，对瑞典维纳恩湖冰面、江苏太湖湖面，以及江苏平地地区的GF-7星载激光数据，分别进行激光高程相对与绝对精度验证，得到相关结论如下：

(1) 本文方法对噪声不敏感，可适用于星载激光各类情况下波形，正确拟合出波峰位置。本文方法提取的测距值，较一般峰值法，测距精度提升了7.5 cm。

(2) 利用瑞典维纳恩湖冰面的GF-7激光点为试验对象，经验证结果表明，经本文方法提取的激光测距值，计算的高程相对精度提升了4.2 cm。

(3) 以平地地区高精度机载LiDAR点云为高程验证数据，本文方法相对一般峰值法提取的测距值，计算的激光高程绝对精度提升了4.5 cm。

综上，本文提出的方法能够抑制激光发射与回波波形噪声，有效地减小全波形星载激光测距随机误差，提升了激光脚点高程精度，目前已用于GF-7星载激光数据业务化处理中。针对不同地物区域的星载激光返回波形引起的测距误差，还需进一步开展深入的研究分析。