Analysis on ranging error of terrain targets for ICESat-2 single-photon lidar

星载单光子激光雷达是近几年发展起来的新一代对地观测载荷，它采用低能量、高重频和多波束激光脉冲^[1]，通过采集地面目标的单光子事件来获取高精度和高分辨率目标信息，从而使星载单光子激光雷达的使用性能得到大幅提升，它也成为未来对地观测的一个重要发展方向。

单光子激光雷达只能探测到信号的有无，不能探测到信号的大小，无法像传统线性激光雷达一样能够采集到完整的回波波形后通过求时间重心来确定激光测距值^[2-3]，因此，其测距精度会受到光子事件随机性的影响。Wang和Zhou通过减小发射脉冲宽度，增加发射激光能量以及使用窄带滤光片等硬件设计措施来提升激光测距精度^[4-5]。Oh和黄科针对静态目标，通过累积几千次的光子事件产生光子直方图，并根据光子直方图的时间重心来估算激光测距值，从而修正光子事件随机性引入的误差^[6-7]。然而，星载单光子激光雷达始终处于高速运动状态，每次测量的目标都是在动态变化的，各种噪声会影响激光测距值的解算。因而，目前提升星载单光子激光雷达测距精度的主要手段还是集中在其点云数据的精细化处理方式上^[8]。但是，经过单光子数据处理后得到的有效信号的点云事件仍然存在一定的随机性。因此，ICESat-2研究团队在牺牲空间分辨率的情况下，针对陆地和植被目标，选择100 m范围内目点云高程的中值作为目标高程值^[9]，而海冰高程产品则是基于累积150个左右的光子事件所产生的光子直方图来获取初始海冰面高程值^[10]，不过没有提供相关的选取依据。

针对光子事件随机性对测距精度影响的问题，在不牺牲测量空间分辨率的情况下，通过累积邻近有限数量的光子事件产生光子直方图，并通过理论推导反算出目标响应函数获取激光测距结果。同时，以测距均方根误差和平均绝对误差为性能指标，构建一种综合考虑陆地地形、累积光斑数目和回波光子数对测距误差影响的分析方法。采用ICESat-2单光子激光雷达的实测数据和产品结果，选择随机起伏的陆地条带，验证激光测距值解算及其误差分析方法，这为单光子激光雷达测距值的获取提供技术支撑。

星载激光雷达所接收的脉冲回波信号可以近似为发射脉冲信号与目标响应函数的卷积。对于常见的线性陆地目标而言，其脉冲回波信号的理论表达式为：

式中：$ g\left(t\right) $和$ m\left(t\right) $分别表示发射激光脉冲和目标响应函数；$ \gamma $、$ {t}_{h} $和$ \delta $则是目标响应函数的面积、时间重心和均方根脉宽，其中$ {t}_{h} $直接反映激光测距值的大小。

单光子激光雷达仅能捕获到离散的光子事件，无法获取脉冲回波信号。然而，单光子激光雷达的重复频率很高，相邻光斑的间距很小，使得一段沿轨距离范围内的目标具有一定的相关性，因而，可以将当前光斑作为中心，按照一定的时间（高程）分辨率，通过累积前后相邻的多个光斑的光子事件来形成光子计数直方图。

将光子直方图进行归一化运算，归一化的光子直方图直接反映出光子事件的概率分布，基于单光子探测模型^[11]，该概率近似满足以下关系式：

式中：$ {P}_{r}\left(i\right) $和$ K\left(i\right) $表示在第$ i $个时间片的光子事件概率和光子数；$ d $表示死区时间所对应的时间片数目。因此，由光子直方图反算可以得到脉冲回波分布为：

考虑到累积光斑数是有限的，则反算得到的脉冲回波分布与理论的脉冲回波之间会存在一定的区别。

基于公式（1），可以反算得到目标响应函数为：

式中：$ f $和$ {f}^{-1} $分别表示傅里叶变换和傅里叶逆变换。

由于常见的线性目标所对应的目标响应函数在理论上满足高斯分布，则采用高斯模型对反算的目标响应函数$ M $进行非线性最小二乘拟合得到高斯函数$ {L}_{n} $。通过解算拟合高斯函数的时间重心^[12]，即可得到激光测距值的表达形式：

式中：$ n $表示累积光斑的个数；$ c $表示光速；$ \tau $表示单光子激光雷达的时间片长度；$ {L}_{n}\left(i\right) $表示拟合高斯函数在第$ i $个时间片对应的数值；$ m $表示时间片的总数。

当然，激光光斑内的实际陆地地形可能会偏离线性分布，导致目标响应函数无法采用高斯函数来表示，因此，在解算激光测距值时，可以将公式（5）中的$ {L}_{n} $替换为公式（4）中的$ M $，这样得到光斑内的陆地目标的平均距离。

星载单光子激光雷达的地形起伏、累积光斑数和平均回波光子数是影响陆地目标测距误差的三个重要因素。考虑到单光子激光雷达的测距误差无法采用解析形式来描述，则采用如图1所示的测距误差评估流程。

Figure 1.  Analysis process of ranging error for single-photon lidar

单光子激光雷达测距误差的评估分析主要包括四个过程：单光子点云仿真^[13]、激光测距值解算、参考距离提取和测距误差分析。其中，单光子点云仿真是对机载激光雷达点云数据所生成的伪回波进行随机抽样后，再经过死区时间剔除所产生的。测距值是由反算的目标响应函数的时间重心解算得到的，参考距离是光斑内目标的加权平均高程。测距误差即为解算测距值与参考距离之间的差异。

参考距离是在考虑激光光斑的高斯空间分布和机载激光点云的强度情况下得到的。以目标的相对高程来表征参考距离，则其表达式为：

式中：$ e\left(j\right) $、$ h\left(j\right) $和$ d\left(j\right) $和分别为第$ j $个机载激光点的强度、相对高程及其与光斑中心的水平距离；$ l $表示当前光斑中机载激光点云的总个数，$ r $表示当前光斑的均方根半径。

采用均方根误差（RMSE）和平均绝对误差（MAE）形式来表征某个机载条带目标的测距误差^[14]，其中，RMSE和MAE分别代表条带目标测距值的随机误差和系统误差。根据累积多光斑事件得到的测距值以及参考距离的表达式，可以得到测距值的RMSE和MAE的形式为：

式中：$ p $表示沿轨光斑点的总个数。

ICESat-2星载单光子激光雷达是目前国际上唯一一个对地观测的光子计数激光雷达，其主要目的在于利用单光子激光雷达的高精度测距结果监测极区冰层的消融情况。单光子激光雷达发射系统和接收系统的主要技术参数如表1所示^[15]。

不同地形起伏会影响激光测距结果，因而选取平坡、缓坡和陡坡三种目标类型^[16]作为分析对象。三种目标所对应的机载激光点云数据分别位于新西兰金海湾区域，美国科罗拉达州西南部滑坡区域和美国爱达荷州岩基泥石流区域，其选取的区域面积大小均为150 m×750 m，三种区域的点密度已经足够描述伪回波波形^[17]，相关数据的主要参数如表2所示。

以ICESat-2单光子激光雷达的强激光作为测量光束，预估陆地目标的平均回波光子数为3，累积光斑数目选取为11，仿真出三种地形的测距误差分布图。

从图2中可以看出，累积光斑前后不同坡度的地形都会对测距误差造成影响。在累积多光斑前，三种不同坡度区域的测距误差值波动性很大，平坡、缓坡和陡坡的沿轨光斑点的最大测距误差值分别达到了100 cm、200 cm和400 cm；在累积多光斑后，误差值的波动趋近于平稳，测距误差值得到明显减小，其沿轨光斑点的最大测距误差值也分别减小到50 cm、100 cm和200 cm。通过统计分析得到，三种坡度区域所对应的均方根误差和平均绝对误差分别见表3所示。

Figure 2.  Range errors before and after accumulating multiple spots at three slope areas

从表3中的结果可以得出，累积光斑前后，随着坡度的增加，条带的测距误差值均会相应增加。且无论平坡、缓坡还是陡坡区域，累积光斑后的RMSE和MAE均能得到改善，三种地形的RMSE值分别减小10.55 cm、28.61 cm和83.91 cm，MAE值分别减小7.17 cm、19.56 cm和64.37 cm。

理想情况下，对于短死区的单光子激光雷达系统而言，当光斑大小固定不变且累积光斑数目趋近于无穷大时，目标响应的时间重心对应于实际的激光测距值。然而，对于有限的累积数，累积直方图分布与回波光子数分布之间出现偏差，其分布形式与累积数有关，这使得不同累积数条件下的激光测距值也会发生改变。

以$\Delta {\rm{\rho }}\left( {{n}} \right) = {{RMSE}}\left( {{n}} \right) + {{MAE}}\left( {{n}} \right)$来表征激光测距误差的大小，模拟出三种目标区域所对应的测距误差与累积光斑数目之间的关系曲线，如图3所示。

Figure 3.  Distribution of the sum values of MAE and RMSE at three slope areas

从图3的仿真结果中可以看出，随着累积光斑数目的增加，三种目标区域的测距误差值都会逐渐减小，但是当累积光斑数目超过某个值时，其绝对误差值又会逐渐增大。三种目标区域所对应的最优累积光斑数目取值不同，对于平坡、缓坡和陡坡，最优的累积光斑数分别为19、23和21。考虑到星载单光子激光雷达测量过程中，无法预知被测陆地区域目标的坡度先验知识，因此，采用以下原则来确定最终的累积光斑数目优化值：

式中：$ \Delta {{\rm{\rho }}_\mu }\left( {\mu = 1,2,3} \right)$分别表示平坡、缓坡和陡坡目标的测距误差。基于该原则得到，对于陆地目标，当累积光斑数为21个时，ICESat-2单光子激光雷达的测距误差达到最小。

考虑到单光子事件仿真随机性会引起激光测距误差的起伏，在固定光斑累积数为21，回波光子数为3的情况下，针对平坡、缓坡和陡坡三种地形，重复100次仿真模拟，得到其测距误差的起伏分布，如图4所示。

Figure 4.  Fluctuation distributions of range error at three slope areas

从图4可以看出，随着仿真次数的变化，三种地形所对应的测距误差会发生微小起伏，其平均值分别为10.86 cm、16.75 cm和66.40 cm，其标准差分别为0.70 cm、1.25 cm和4.22 cm，对应的变异系数（标准差/平均值）分别为6.45%，7.46%和6.36%。上述结果表明，测距误差的起伏远小于其均值水平，即累积光斑数目后的测距误差具备收敛性。

若陆地目标的回波光子数发生改变，则最终的优化累积光斑数目和对应的测距误差也会相应地出现变化，图5仿真模拟了不同回波光电子数条件下的最优累积光斑数目的分布。

从图5的结果可以看出，随着平均回波光子数的增加，最优累积光斑数会逐渐减小，当平均回波光子数超过三个时，累加光斑数趋于稳定。当平均回波光子数为10时，对应的累积光斑数目为15，这与ICESat-2提供的ATL07产品中海冰累积150个光子是相吻合的，这也证明了文中累积光斑数优化方法的正确性。

Figure 5.  Relationship between the optimal number of accumulative spots and the average echo photons

在最优累积光斑数目的情况下，不同回波光子数对应的最小测距误差也是不同的，图6给出了三种区域平均回波光子数与测距误差的关系图。

Figure 6.  Relationship between the average echo photons and range error

图6的结果表明，随着回波光子数的增加，相应的测距误差值会逐渐减小，当回波光子数超过四个时，三种地形的测距误差会趋于稳定，因此，适当地增加回波光子数量有利于提高测距精度。

以ICESat-2卫星过境美国犹他州西瓦利城附近陆地区域（40°42′ N, 112°7′ W）的一个条带数据作为分析对象，其条带的总长度为420 m，对应的光斑数为600个。其中，滤噪后的光子事件数据及其对应的目标高程数据分别存储在ATL03和ATL08产品中。根据ALT03产品中光子事件的经纬度信息，选取相应的机载激光雷达点云数据作为参考源，其点密度为11.93 pts/m²，如图7所示。

Figure 7.  Study area with corresponding true elevation and slope distribution

通过对图7沿轨光斑的坡度统计分析，光斑内目标坡度为0°~5°、5°~20°和20°~35°的光斑数分别占总光斑数的3%、70%和27%，其地形存在不规则的随机起伏特点。

为减少噪声对测距误差的影响，仅选取ALT03产品中置信度标记为4的光子事件进行处理分析，该条带上的平均回波光子数为2.2个，对应的累积光斑数为21。采用累积21个光斑的光子事件所得到的光子直方图，解算得到沿轨条带上的目标高程。将该高程与累加前的目标高程以及ICESat-2提供的ATL08产品中的目标高程进行比较，如图8所示。

Figure 8.  Comparison of elevation before and after the accumulative spots and the stripe elevation label of ATL08 products

通过图8显示的结果计算分析，累积光斑前后RMSE值由114.25 cm减小到了63.84 cm，MAE值由70.97 cm减小到了48.52 cm，两者之和减小的幅度达到72.86 cm。将ATL08产品提供的间隔100 m范围的A、B、C、D和E这五个高程值点与目标参考高程值比较，计算得到其RMSE值和MAE值分别是137.96 cm和97.24 cm，两者的值均大于累积21个光斑后该条带测出的结果值。

相比于仿真的误差分析结果，实际地形的测距误差结果略偏大，这主要是由于ICESat-2过境的这片条带地形的相关性不足所导致的，使得某些位置处的目标回波偏离理想的单高斯回波信号，存在多个模式回波的混叠。

文中从星载单光子激光雷达测距模型出发，提出了一种通过累积直方图反算出目标响应函数的时间重心来确定激光测距值的方法，分析了陆地地形、累积光斑数目和回波光子数对测距误差的影响。从数值模拟和实际数据验证的角度，论证了在不降低采样分辨率的条件下，该方法的正确性与可行性。

通过仿真分析表明，测距精度的提升与目标起伏程度有关，目标越平坦，测距精度提升效果就越好，反之亦然；累积光斑数的优化结果与目标平均回波光子数存在反比关系，当平均回波光子数超过三个时，累积光斑数的优化值趋近于平稳。随着平均回波光子数的增加，测距误差会逐渐减小。同时，多次重复仿真表明，所提出的激光测距解算方法具有较好的收敛性。然而，激光测距精度提升效果受地形起伏的制约。若地形无法采用单一的平面形式来描述，其测距精度的提升效果就会下降。

所提出的激光测距解算方法有利于提高单光子激光雷达在陆地目标的测距精度，这对我国未来单光子激光雷达数据处理具有一定的参考意义。