Spaceborne photon-counting LiDAR on-orbit calibration based on natural surface

星载激光雷达是一种主动式的激光测量设备，它根据激光脉冲的渡越时间(Time of Flight, ToF)获得卫星与地表目标间的精确距离值，结合卫星平台的精确姿态、位置信息以及激光指向信息后可以获得目标的精确三维坐标。这种对全球尺度的地表高精度三维信息的获取能力，使得星载激光测高在极地冰盖及海冰检测、陆地高程测量、植被高度及生物量估计以及全球三维控制点获取等方面能够发挥重要作用^[1]。

对于星载激光雷达，测距精度^[2-3]以及激光光束指向的精度是影响目标的三维坐标定位精度的两个重要因素，其中光束指向精度尤为重要^[4-5]；以600 km的轨道高度以及1°坡度的地表为例，在激光光束近似沿天底方向发射的情况下，30″的指向误差将会使激光光斑坐标在平面方向产生87 m的误差，并在高程方向产生1.5 m的误差^[6]。

可见，激光光束指向信息对地表光斑坐标的精确解算有着极为重要的意义。一般来说，卫星发射前会在地面对包括激光光束指向、测距零点在内的各项参数进行精确测量。然而，在发射过程中以及在空间中运行时的恶劣环境将不可避免地使激光光束指向以及测距零点产生一定的系统误差。因此，在这种情况下，对激光光束指向的系统误差以及测距系统误差进行在轨标定就显得十分有必要。美国国家航空航天局(National Aeronautics and Space Administration, NASA)在2003年发射ICESat卫星(Ice, Cloud, and Land Elevation Satellite)后，对其各项系统误差开展了细致的在轨标定及检校工作，从而实现了GLAS (Geoscience Laser Altimeter System)数据测距5 cm、指向(包含横滚和俯仰项)1.5″(1σ)的精度^[7-8]。我国的在资源三号02星发射后，也使用地面探测器对其搭载的激光测高实验载荷开展了在轨误差标定和检校工作，使得激光测高实验载荷的高程测量精度从标定前的6.3 m提高到了标定后的1.6 m^[9-12]。

近年来星载激光遥感技术得到了快速发展，光子计数星载激光雷达有微脉冲、高重复频率、多波束、小光斑等显著优点，成为星载激光雷达发展的新趋势^[13-14]。目前我国针对线性体制激光雷达提出的使用地面探测器进行在轨标定和检校的方法^[9-11]，在光子计数星载激光雷达的在轨标定中存在诸多不便之处。例如，光子计数星载激光雷达极高的重复频率要求与之配合的地面探测器要具有极高的时间同步精度，极低的能量密度要求大幅提升探测器的灵敏度，多波束和小光斑的特点会使得地面探测器的布设范围和布设密度大大增加^[15]。因此，对于光子计数星载激光雷达，有必要寻找一种新的在轨标定方法。文中将结合NASA所采用的海面姿态机动扫描的在轨标定技术，提出一种基于自然地表的在轨标定新方法，根据对两个实验场地的标定试验结果，对标定场地形对在轨标定结果的影响进行分析和讨论。

星载激光雷达在轨标定的目的是剔除测量过程中的系统误差，从而提高获得的地面光斑坐标的精度。在星载激光测高中，地面激光脚点如公式(1)所示。

式中：P_g^ECEF为地面激光脚点在地心地固坐标系(Earth-Centered Earth-Fixed, ECEF)下的坐标；r_s^ECEF为卫星位置数据；R_SBF2ECEF为由卫星本体坐标系(Satellite Body Frame, SBF)到ECEF坐标系的变换矩阵，这一矩阵表示卫星姿态；ΔX ^SBF为激光器相对于SBF原点的偏移量；R_laser2SBF为激光器在SBF下的安置姿态，一般情况下可以用俯仰(pitch)和横滚(roll)两个欧拉角所对应的旋转矩阵表示R_laser2SBF；ρ′为星载激光雷达实际测距数据；Δρ为测距零点的系统误差。

可见在ECEF坐标系下的激光指向受r_s^ECEF、R_SBF2ECEF、ΔX ^SBF、R_laser2SBF的影响；由于经后处理后的卫星姿态与位置数据（R_SBF2ECEF、r_s^ECEF）能够达到很高的精度，并且由于ΔX ^SBF远小于测距值ρ，由ΔX ^SBF的误差所导致的光斑定位误差极小，因此在对星载激光雷达进行在轨标定时，主要关注的是对R_laser2SBF这一旋转矩阵以及测距值ρ的标定。其中，R_laser2SBF可进一步分解为以俯仰和横滚两个方向的欧拉角所表示的旋转矩阵之积，ρ可以分解为实测数据ρ′与测距零点的系统误差Δρ之和。因此，在星载激光雷达的在轨标定中，激光器在卫星平台中的俯仰角、横滚角以及测距系统误差(Δρ)就是需要计算的3个基本参数。

基于地面探测器的在轨标定技术直接通过地面探测器捕获激光光束在地表形成的光斑，然后根据公式(1)做间接平差得到系统误差参数。基于海面姿态机动扫描的在轨标定技术是利用海面平均高程这一先验信息，根据测距残差解算系统误差参数。基于自然地表地形匹配的在轨标定技术则是以地表高精度的高程数据为先验信息，根据回波波形相关系数或高程残差解算系统误差参数。

2003年美国ICESat卫星发射后，NASA使用基于姿态机动和海洋表面扫描的方法对激光载荷进行了在轨标定。这一方法将海洋表面假设为具有一定粗糙度的理想平面，在已知卫星精确轨道位置以及海洋表面高程后，利用实际星载激光雷达的测距值与理论测距值之间的残差对系统误差进行解算。

图1所示为海洋表面激光测距几何原理的简单示意图，图中θ表示光束指向，σ_θ为指向角变化量，s表示表面坡度，σ_ρ为测距值变化量。测距值关于指向的导数σ_ρ/ σ_θ是一个关于θ和s的函数。在卫星进行姿态机动时σ_ρ/σ_θ会随卫星姿态发生改变，从而使得分别由指向系统误差和测距系统误差导致的测距值的误差彼此分离；再利用由海面高程计算得到的理论测距值与星载激光雷达的实际测距值之间的残差对系统误差参数进行迭代逼近。

Figure 1.  Ocean surface laser ranging geometry principle

以姿态机动的方式对海面进行扫描的在轨标定方法，对标定场的要求较低、可选区域范围广，因而可以在卫星运行轨道的较大弧段内进行标定工作，这使得系统参数迭代解算中的样本数量可以达到很大的量级，对于有效估计系统误差十分有利。并且，这一标定方案只需要低分辨率的海面平均高程和海面平均浪高数据，这些数据大都已有成熟的数据库，获取难度和获取成本也很低。

该方法不仅适用于以GLAS为代表的线性体制星载激光雷达，也同样适用于星载光子计数激光雷达。2018年美国发射的ICESat-2卫星中搭载了世界首个光子计数星载激光雷达ATLAS (Advanced Topographic Laser Altimeter System)，该系统也采用了基于姿态机动和海洋表面扫描的方法对指向与测距系统误差进行在轨标定^[16]。然而该方案中，卫星在进行姿态机动时可能会引入额外的误差从而对系统误差的计算产生影响，这就要求卫星平台的姿态测量和姿态控制必须达到极高的精度^[17]。这一特点导致基于海面姿态机动扫描的在轨标定方案不适用于我国的大平台、三轴姿态稳定卫星平台。

基于自然地表的在轨标定方法是一种利用先验的地表地形数据，如数字高程模型(Digital Elevation Model, DEM)或数字表面模型(Digital Surface Model, DSM)进行在轨标定的方法。其基本原理是通过寻找地表DEM或DSM与卫星下传测量结果之间的高程残差、或波形残差（针对线性体制激光雷达）的最小值点，求解星载激光雷达的各项系统误差参数。

在对线性体制的星载激光雷达在轨标定方面，唐新明等人提出了一种针对资源三号02星的利用大范围DSM进行指向角粗标定的方法^[12]，该方法利用经地表DSM(ALOS全球数字表面模型ALOS World 3D-30 m)，将星载激光雷达的高程测量精度由几十米提高到了3 m以内。唐新明、谢俊峰等人还提出了一种利用金字塔搜索的地形匹配在轨标定方法^[18]，该方法根据卫星实测数据与DSM之间的高程残差，利用金字塔搜索的方法，成功地将高程误差由100 m降低到了3 m左右。张文豪、李松等人提出了一种针对具备波形记录功能的线性体制星载激光雷达的在轨标定方法^[19]，该方法根据由DSM得到的仿真波形与卫星实测波形之间的相关系数求解系统误差参数，在使用GLAS数据验证时，实现了约2 m的平面定位精度。易洪、李松等人提出了一种基于自然地表测距残差的在轨标定法思路^[20]，并使用GLAS数据进行了仿真验证。

对于星载光子计数激光雷达，其更小的光斑间距使得在同样的标定场内能够获得更多的数据，这有利于系统误差参数的精确计算。Yaming Nan等人根据卫星实测数据与DEM之间的高程残差，通过最小二乘法求解其最小值点，实现对系统误差参数的解算^[21]。在利用人为加入坐标偏移后的ICESat-2点云做验证时，该方法得到的点云坐标与ICESat-2 ATL03数据中提供的官方点云坐标之间的平面误差约为3 m，这符合ICESat-2卫星6.4 m的平面定位精度。然而，该方法还存在一定问题：在对表示激光光束指向的两个角量的标定中，只有一个角量标定的较为精确，而另一个角量仍存在较大的误差。

根据前人的工作经验，结合光子计数星载激光雷达技术的工作特点，笔者提出了以测距残差为代价函数，通过最小二乘法求解星载激光雷达系统误差的新方法。

笔者使用公式(2)所示函数近似表示星载激光雷达对光斑范围(X±X₀ , Y±Y₀)内目标的测距结果，其中，$ \overrightarrow{\alpha }$表示待估的系统误差参数，L(X+u, Y+v)表示星载激光雷达发射天线与地表微面元之间的距离，公式(3)所示H(X+u, Y+v)表示激光雷达接收孔径接收到的能量大小，κ为一个受地表反射率、大气透过率、激光雷达光学系统效率影响的常数，n_g和n_l分别表示地表面元的法向量和激光指向向量，σ_B表示地表光斑的1σ半径。考虑到各种因素导致的测距随机误差后，构建了公式(4)所示的星载激光雷达测距值的观测函数；然后利用最小二乘估计的得到系统误差参数的估计结果，最小二乘递推公式如公式(5)所示。

受篇幅限制，该方法详细的理论模型和实现算法将另文介绍。目前，使用仿真点云数据和ICESat-2数据进行验证试验并取得了初步的结果。

首先，使用仿真点云数据对标定算法进行验证。验证过程中，首先根据卫星姿态、轨道数据求得指向$ \overrightarrow{n}$，然后在$ \overrightarrow{n}$的基础上人为地加入指向误差后（包括系统误差及随机误差）的得到$ {\overrightarrow{n}}_{0}$，并根据$ {\overrightarrow{n}}_{0}$计算出实际地表光斑中心坐标，然后依据对应区域的地表DEM数据，按照参考文献[22]中所示方法对激光雷达所接收的光子进行仿真，得到每一个接收光子所对应的飞行时间(Time of Flight, ToF)真值t₀，并在t₀的基础上加入测距随机误差与测距系统误差后得到飞行时间t，作为接收光子飞行时间的实测值。将结合$ \overrightarrow{n}$、t以及卫星的姿态、轨道数据计算得到的单光子点云坐标$ \overrightarrow{P}$作为仿真得到的带有误差的“实测结果”。仿真所使用的相关参数及误差大小见表1，在仿真过程中，笔者在激光指向角中加入了200″量级系统误差与2″的随机误差，并在测距结果中加入了10 m的系统误差与1 m的随机误差。

在通过仿真得到单光子点云坐标$ \overrightarrow{P}$后，将$ \overrightarrow{P}$输入标定算法，将经算法解算出的系统误差参数与表1中设置的真实系统误差参数进行对比，评估标定精度。

图2及图3是给定不同迭代初值时，标定结果与仿真预设参数之间的误差变化。从图中可以看到，指向角系统误差的迭代初值与表1所示的真实系统误差参数相差在−200″~+200″的范围内时(500 km轨道高度下对应约±500 m的平面误差)，由不同初值得到的标定结果与表1所示真值之间的残差基本稳定在0.2″(指向角)及1 cm (测距)以内，这表明，该标定算法对于不同的初始值都能得到较为稳定的标定结果。

Figure 2.  Calibration results with different initial Roll values

Figure 3.  Calibration results with different initial Pitch values

图4所示为使用不同长度的数据进行标定时，标定结果与表1所示的真实系统误差参数之间的残差的变化，用于验证不同样本大小对标定结果的影响。可以看到在数据长度达到2~3 km以上时，标定结果中指向角的精度基本可以达到0.5″，测距系统误差Δρ的标定精度可以达到分米量级。仿真验证试验证明了新标定算法的正确性。

Figure 4.  Calibration results with different length of datasets

除了使用仿真点云进行标定算法验证外，还使用了ICESat-2的实测数据进行了验证试验。

由于ICESat-2本身所提供的点云坐标数据是经过标定改正后的点云坐标，因此，笔者认为这一数据中包含了实际物理过程中产生的随机误差以及标定后残余的部分系统误差，为了尽可能地还原标定前数据中存在较大系统误差的情况，首先利用NASA所提供的点云坐标$ \overrightarrow{P}$以及卫星姿态、位置数据，根据公式(1)所示的激光脚点定位方程，反推出激光器在卫星本体坐标系下的安置姿态矩阵R_laser2SBF以及测距值ρ，并加入200″的指向系统误差以及10 m的测距系统误差，从而得到一组新的带有较大系统误差的点云坐标$ {\overrightarrow{P}}^{\prime }$。然后将$ {\overrightarrow{P}}^{\prime }$输入标定算法进行改正得到$ {{\overrightarrow{P}}^{\prime }}_{0}$，通过将$ {{\overrightarrow{P}}^{\prime }}_{0}$与NASA官方提供的标定改正后的点云坐标$ \overrightarrow{P}$进行对比，对标定算法进行验证和评估。

N折交叉验证是一种广泛应用于数学建模和机器学习中的一种模型评估方法，其基本原理是将数据集分为训练集和评估集两个相互独立的部分，使用训练集对模型中的关键参数进行求解，并使用求解出的参数进行预报，然后通过预报结果与评估集中的数据之间的对比对模型性能进行评估。

在此，笔者也是用交叉验证的方法对标定算法进行评估。交叉验证中，首先将数据集分为N段，任取其中1段作为验证数据，剩余N-1段作为输入标定算法的训练数据，这样便得到了N组相互独立的验证数据与训练数据。使用交叉验证的方式保证了验证数据与训练数据之间的独立性，使用训练数据通过标定算法得到系统误差参数后，将系统误差参数带入激光脚点定位方程中对验证数据进行改正，并将改正后的数据与NASA提供的ATL03数据中的点云坐标进行对比，从而对标定算法进行评估。

表2为使用南极麦克莫多干谷地区的DEM数据及对应的ICESat-2实测数据（ATL03数据中陆地类型置信度大于3的点云数据）对标定算法进行交叉验证试验的结果。试验中，将6 km的ICESat-2数据平均分为6段，然后从中取5段作为标定算法的输入数据，再利用标定算法输出的系统误差参数解算剩下1段数据的点云坐标并与NASA的ATL03数据中自带的点云坐标进行比较。

表2中，Δx、Δy表示由算法标定校正后的点云与NASA ATL03数据中自带的点云在局地坐标系(此处使用站心地平坐标系)的x轴、y轴方向上的差异，用来表示平面方向的偏差；Δz表示由算法标定校正后的点云与NASA ATL03数据中自带的点云在局地坐标系(此处使用站心地平坐标系)的z轴方向上的差异，用来表示高程偏差。可见平面方向相对于ATL03提供的点云坐标偏差大约在3 m左右，这符合NASA公布ICESat-2的6.4 m平面定位精度的设计指标，也与参考文献[21]中的得到的结论一致。文中的标定算法改正后的点云与NASA公布的点云坐标之间仍然存在一定偏差，这是由多种因素导致的，包括文中的标定算法改正后的系统误差残差、NASA所提供的官方点云数据中的系统误差残差以及部分的随机误差，但是二者之间3 m的偏差小于NASA设计指标中6.4 m的平面定位精度可以说明笔者所提出的标定算法的标定精度与NASA所使用的标定算法的标定精度是处于同一量级上的。

此外，为了更进一步地检验标定算法的效果，文中还使用了连云港附近的高程变化较大、地形特征较为丰富的地区进行标定实验；与南极麦克莫多干谷不同，连云港地区存在一定的植被以及人工建筑。因此，需要首先对解算系统误差参数过程中使用的ICESat-2数据和地表高程数据进行一定的处理。对于ICESat-2数据，结合ALT08数据中对光子点云的分类标签(分类标签包括：“噪声”、“地表”、“冠层”和“冠层顶”)，去除ATL03数据中来自噪声和植被冠层反射的光子点云，只保留地表反射的光子点云。对于地表高程数据，通过点云处理软件对机载LiDAR点云进行分类，去除植被冠层点云，只保留地面点云数据，并依据地面的机载点云生成不含植被冠层的DEM。标定中使用去除植被冠层影响后的ICESat-2数据和DEM数据作为输入，求解系统误差参数。利用标定得到的系统误差参数对ICESat-2点云改正后，将人工建筑作为特征点，通过对比改正后的ICESat-2数据中人工建筑反射的光子点云与地面DSM数据（为了确保DSM中包含人工建筑顶部的高程，该DSM数据与前文解算系统误差参数过程中使用的DEM数据不同，该DSM数据是直接由包含植被冠层数据在内的未分类的机载LiDAR数据生成的），可以直观地判断在轨标定的效果。

表3为使用连云港云台山地区的DSM数据及对应的ICESat-2实测数据对标定算法进行交叉验证试验的结果。图5与图6分别是文中标定算法改正后的点云以及NASA提供的官方点云与对应平面坐标处的DSM高程的对比；图5中，标定后的点云是在对ATL03点云数据使用DBSCAN算法滤波后（DBSCAN算法搜索直径为3 m，为了防止错误地去除建筑物顶部反射的光子点云，这里并没有像前文解算系统误差参数过程中所使用的数据一样利用ATL08数据对植被冠层反射的光子点云进行滤除），结合标定得到的系统误差参数后改正的结果。其中，图5所示数据中，经过文中的标定方法改正后的点云与对应平面坐标处的DSM高程之间的均方根误差为2.22 m；图6所示数据中，NASA官方提供的点云与对应平面坐标处的高程之间的均方根误差为2.94 m，相较之下，经文中标定方法改正的点云数据与用DSM表示的地表“真值”之间的误差更小。同时在图5、图6中，可以明显地看出,对于黑色圆圈所框的四处目标，文中标定算法改正后的点云与DSM的重合程度要好于NASA官方提供的点云与DSM的重合程度。这表明标定算法得到了一个正确的匹配结果。

Figure 5.  Comparion of calibrated point cloud with DSM elevation

Figure 6.  Comparion of NASA official point cloud with DSM elevation

由表2、表3所示结果可以看出：使用文中提出的标定方法校准后的点云坐标相较于NASA提供的ATL03数据中的点云坐标的平面偏差分别为3 m和3.5 m，对应的指向角偏差约为1.23″和1.44″。需要说明的是，对于ICESat-2卫星来说，笔者无法知道地表光斑的真实位置，并且NASA提供的ATL03数据中的点云坐标本身也存在一定误差(ICESat-2平面定位精度的设计指标为不低于6.4 m)，因此，对于表2、表3中所提供的结果，其相对于ATL03的点云坐标平面偏差不大于6.4 m，就可以认为文中的标定算法已经基本达到了ATL03数据的精度水平。

图5与图6使用点云高程与DSM的对比，定性地比较ATL03数据与使用文中标定算法校正后的点云数据之间的精确度。图5、图6中所示4个目标处，可以看到：由文中所提出的标定算法校正后的点云与DSM的重合程度要明显好于ATL03数据所提供的点云与DSM的重合程度。从这一点可定性的确认，在该标定场范围内，由文中所提出的标定算法校正后的点云坐标要比ATL03数据所提供的点云坐标有更好的平面精度。

另一方面，对于基于自然地表的标定方法来说，标定场本身的地形也是影响标定精度的重要因素。对于标定场地形对标定精度的影响，比较表2和表3，可以发现南极麦克莫多干谷地区的标定结果的稳定度(以交叉验证试验结果的标准差衡量)要好于连云港地区的标定结果。对于测距精度来说，图7与图8分别为连云港地区和麦克莫多干谷地区标定后测距值的残差（标定后的测距值与由DSM得到的测距值之差），可以看到南极麦克莫多干谷地区的测距残差要明显小于连云港地区的测距残差。

Figure 7.  Ranging residual after calibration in Lianyungang

Figure 8.  Ranging residual after calibration in Antarctic McMudro Dry Valleys

受模型的非线性性以及标定场复杂地形的影响，定量地描述标定精度与各因素之间的关系是一个很困难的问题。影响标定精度的因素不仅有作为输入数据的卫星姿态、轨道、测距值的误差大小，标定场地表类型（如植被、人工建筑等）以及坡度、坡向、地形的复杂程度等因素同样会影响基于自然地表的在轨标定方法的标定精度。笔者将在之后的工作中对标定精度与各影响因素（包括输入数据噪声、标定场地形地貌等）之间的关系进行更进一步的深入研究。

在轨标定是一项影响星载激光雷达光斑定位精度的重要技术，文中对星载激光雷达在轨标定的主要目标和核心理论模型进行了阐述，对现有主要在轨标定工作的基本原理、实现方法和技术特点进行了较为深入的分析。文中对基于自然地表的在轨标定方法进行了重点讨论。基于自然地表的在轨标定，利用地表已知先验地形数据与卫星激光雷达下传测量结果之间的高程残差、或波形残差求解星载激光雷达的各项系统误差结果。相较于其他类型的星载激光雷达在轨标定方法，基于自然地表的在轨标定方法实施成本较低，不需要在地面布设大规模探测器等辅助标定设备，同时对卫星平台本身没有特殊要求，是一种具有较好推广价值的在轨标定技术。

文中在目前已有在轨标定方法的基础上提出了一种新的基于自然地表的地形匹配在轨标定方法，其适用于光子计数星载激光雷达的高精度在轨标定。文中使用ICESat-2数据，以南极麦克莫多干谷地区和中国连云港地区为例，通过对标定后平面偏差和高程偏差的统计，证明了标定理论和方法的可行性和正确性。同时，对文中所示结果的精度进行了更进一步的分析，指出了标定场地形会对基于自然地表的在轨标定方法的精度产生影响。文中对现有星载激光雷达在轨标定方法进行了梳理和总结，并提出了一种新的基于自然地表的在轨标定方法，对我国未来光子计数型激光测高卫星的在轨标定工作具有重要参考意义。