Preprocessing method of laser ranging data based on array detection technology

卫星激光测距（Satellite Laser Ranging, SLR）技术是精度最高的现代空间大地测量技术之一^[1]，空间碎片激光测距（Space Debris Laser Ranging, DLR）和月球激光测距（Lunar Laser Ranging, LLR）是在卫星激光测距技术上发展而来。激光测距涉及多个学科领域，各学科的迅速发展推动着激光测距技术不断更新换代，越来越多的新技术被应用到激光测距技术中，美国APOLLO站曾将4×4APD阵列用于月球激光测距技术中并取得了十分满意的结果^[2-3]，2017年，中国科学院云南天文台成功将阵列探测技术用于卫星/空间碎片激光测距试验中。

阵列探测技术的应用，使得每个激光脉冲中被探测并记录的回波光子数增加，每个主波时刻可能会对应多个回波时刻；系统中使用的阵列超导探测器工作在自由状态，无距离门控，一直处于接收回波信号状态，除在探测器处于恢复时间的情况外，其它时间每来一个光子信号或噪声信号皆可能响应并输出信号，因此，数据量将大大增加；探测器的响应时间抖动^[4-5]、探测器各像元响应不一致、多通道数据采集系统各通道间测量误差等也会对激光测距数据产生影响。以上这些问题对数据预处理方法提出了新的要求。

文中结合阵列超导探测器特点，在常规卫星激光测距数据预处理方法基础上，利用多高斯拟合和互相关分析的方法修正多通道数据间的偏差，利用泊松统计分析、循环拟合滤波的方法进行信号识别、数据拟合，实现了基于阵列探测技术的激光测距数据通道间的时间偏差修正、自动识别、数据拟合等。

文中所使用的数据来自中国科学院云南天文台基于阵列探测技术的异地收/发激光测距试验平台，如图1所示。测距开始后，激光器发射激光脉冲并通过53 cm双筒望远镜发射光路指向卫星/空间碎片，主波探测器探测发射激光并输出主波脉冲信号，并被多通道激光测距采集与处理系统响应并记录此主波信号时刻。激光脉冲里被卫星/空间碎片反射回地面测站的少部分光子信号通过1.2 m望远镜接收光路传输至阵列超导探测器产生回波脉冲信号，然后被多通道激光测距数据采集与处理系统响应并记录各通道回波信号时刻。多通道激光测距数据采集与处理系统包含10个GT668PXIe-1模块，每个模块有2个采集通道；其中9个工作模块，1个热备份模块，最多支持18个通道进行数据采集；GT668PXIe-1模块时间分辨率为0.9 ps，多模块间的时间同步偏差小于50 ps，可由2.4节方法一同修正。

Figure 1.  Diagram of laser ranging system based on array detection

激光测距是测量激光脉冲在测站与目标之间往返飞行的时间，从而计算得到目标的距离值，测量得到的是时间数据。为了使测到的数据能够被用于科学研究，需要对原始数据进行一系列的处理，这个处理过程称之为数据预处理^[6]。

基于阵列探测技术的激光测距数据预处理方法流程：（1）主回波匹配，将各通道每个回波时刻与对应的主波时刻匹配起来；（2）信号识别，从大量噪声信号中识别激光测距信号；（3）循环拟合滤波，利用多次多项式拟合的方法，进一步滤波剔除野点；（4）通道间的时间偏差修正，利用地靶数据或互相关分析，修正通道间的时间偏差；（5）修正后数据处理，进一步提高数据精度。

该节数据预处理过程以2019-04-15日测量的38355空间碎片处理过程为例，多通道激光测距数据概况如图2所示，图2（a）为激光器工作频率，图中可以看出激光脉冲发射频率为100 Hz，中间的缺口表示激光器被暂停过；图2（b）为目标在测量时间内的理论距离，根据主波时刻和距离轨道预报数据进行插值计算得到；图2（c）为各通道采集数据计数，即探测器响应计数；图2（d）为各通道主回波匹配上的数据计数，匹配门限−60 ~60 μs；从图中可以看出，测量数据中1、2、8通道未采集到数据；16通道采集到数据，但未匹配到数据。

Figure 2.  Multi-channel laser ranging data statistic

得到激光测距测量数据后，首先，需要从这些数据中将每个回波时刻与对应的主波时刻匹配起来，这样才能得到激光脉冲的往返飞行时间。匹配的流程是，对于各通道每个回波时刻，与每一个主波时刻做减法，其差值再减去此主波时刻所对应的目标预报距离时间，得到距离残差（O-C）数据^[6]，图3（a）为各通道主回波匹配的O-C数据图。

Figure 3.  Main echo matching O-C plot. (a) Each channel, (b) all channels

得到各通道O-C数据后，先将各通道O-C数据进行简单叠加，得到所有通道O-C数据，简单叠加可以增加有效信号点数并平均各通道噪声，方便下一步的测距信号识别，图3（b）为所有通道O-C数据。

在激光测距过程中，由于系统噪声、环境噪声等的影响，测量的数据中有大部分噪声，需要从噪声中识别出信号数据^[6]。此过程中，处理的数据对象是2.1节中获得的测距O-C数据。

回波识别可以凭借人眼经验判断识别，也可以采取自动识别算法^[7]。泊松（Poisson）统计滤波算法能对高噪声、低信号水平的激光测距数据进行有效滤除噪声。美国McDonald站将Poisson统计滤波算法应用到月球激光测距中并取得了良好的效果，澳大利亚Stromlo站将Poisson统计滤波算法应用到卫星激光测距中并实现了无人职守的全自动观测^[8-11]。

Poisson统计滤波算法认为O-C数据中的噪声点服从Poisson分布，而有效信号则在短时间内近似呈线性分布，通过判断数据点是否服从泊松分布来滤波去噪，保留有效信号^[7-10]。

Poisson分布给出了在给定的单位时间间隔内，如果在总体样本中平均有$\lambda $个随机事件发生，那么单位时间内发生$n$个随机事件的概率如公式（1）所示：

图4（a）为分段进行Poisson统计滤波算法进行信号识别的过程，图4（b）为信号识别结果。将整段数据按照同样的间隔划分为8段，图4（a）为其中第3段。

Figure 4.  Segmental signal identification (3/8) and signal identification results (b)

图4（a）中左上图为O-C数据的直方图统计结果；右上图为根据直方图统计结果得到的频数分布和Poisson分布统计结果，从图中可以看出，数据中噪声频数分布统计结果和Poisson分布统计结果基本一致；左下图为累计概率分布统计结果，当累计概率大于相应阈值时，认为该区域信号为疑似信号或确定信号，待进一步判断，其余区域信号为噪声；右下图为分段信号识别结果。

从图4（b）可以看出，使用泊松滤波的方法可以很好地识别出激光测距信号，并且不受O-C数据图中竖条纹干扰，根据试验系统情况和测距数据分析，该竖条纹可能是激光后向散射被超导探测器探测，也可能是由超导探测器自身干扰引起。

为了进一步滤除噪声，并且得到数据点的分布特征以及计算数据的符合精度，需要对数据进行多项式拟合滤波^[6]。

对识别后的O-C使用最小二乘拟合法，进行低阶多项式拟合，去除系统趋势项，避免引进高阶震荡项，获得拟合残差数据。

将O-C拟合成：

式中：$t$为主波时刻；${f_{{\rm{o} }- {\rm{c}}}}(t)$为多项式拟合的O-C值；${a_i}$为拟合得到的系数。

拟合后，计算距离残差（O-C）$\Delta {R_i}$在时间${T_i}$时候的拟合残差为：

计算拟合残差序列的标准偏差，删除拟合残差值大于2.5倍标准偏差的数据，然后计算剩余数据的内符合精度（标准偏差RMS），若精度没达到要求，拟合阶数增加1，直至达到相应要求。

图5（a）为循环拟合滤波结果，图5（b）为拟合残差图，由多项式拟合去除系统趋势项得到。

Figure 5.  Cyclic fitting filtering results (a) and fitting residual plot (b)

由于发射的激光脉冲服从高斯分布，所以认为回波信号也服从高斯（Gauss）分布，回波光子分布示意如图6（a）所示，由于激光飞行过程中目标反射影响，同时受探测器和数据采集设备影响，回波光子分布会有展宽。

Figure 6.  Echo photon distribution diagram of laser ranging data (a) and time offset diagram between channels (b)

由于阵列探测器各像元响应不一致、探测器到数据采集终端连接线长短不一致、多个GT668PXIe-1模块时间同步不一致等原因，导致多通道激光测距数据各通道间存在一定时间偏差，如图6（b）所示，因此需要对通道间的时间偏差进行修正。

通过地靶测量的方法可以在修正系统延迟的同时，修正各通道间的时间偏差。实现方法是经过上述处理后，分通道计算系统延迟；然后在探测目标数据处理时，分通道减去系统延迟，就可以在修正系统延迟的同时修正通道间的时间偏差。

这种方法简单易行，但需要频繁测量地靶；在地靶测量时，难以保证所有通道都能采集到数据，实际应用中有可能需要舍弃一些通道的数据，只保留地靶测量和目标测量都采集到数据的通道。所以通过地靶测量方法来修正通道间的时间偏差存在一定问题。

采用互相关分析和多高斯拟合的方法也可修正各通道间的时间偏差。如图6（a）所示，回波光子分布可以通过多高斯拟合的方法进行拟合。文中使用两个高斯函数叠加来拟合回波信号分布情况，如公式（4）所示：

对于任意两个通道，理论上认为分布特征应该非常相似或一致，对任意两通道回波分布进行移位互相关分析，可以得到一组互相关序列。其中互相关数值最大的位置表示两通道回波分布相关程度最高，因此此时的位置移动量即为这两通道间时间偏差。

实际使用中，首先选取一个通道作为标准通道，目前通常选取1通道为标准通道（此数据采集通道和主波数据采集通道为同一块数据采集卡，偏差相对较小）或选取数据最多的通道为标准通道（数据分布相对更明显）。然后让每个通道数据和标准通道进行互相关分析。

不同通道间互相关分析相关系数为（$\tau $为通道间的延迟网格个数）：

归一化互相关分析相关系数为：

此过程中处理的数据为2.3节获得的拟合残差数据，如图5（b）所示。

图7（a）为各通道拟合残差直方图统计结果和双高斯拟合结果，图7（b）为其中通道13拟合残差直方图统计结果和双高斯拟合结果。

Figure 7.  Fitting residual histogram statistics and multi-peaks Gaussian fitting results. (a) Each channel, (b) channel 13

根据双高斯拟合结果对各通道与标准通道进行互相关分析，图8(a)为各通道与标准通道（14通道）直方图统计结果的互相关序列和高斯拟合结果（公式（7）），对公式（7）求导数得到公式（8），当拟合函数导数为0时，即$x = {b_1}$时，$\dot r\left( x \right) = 0$，$r\left( x \right)$最大，如图8(b)所示，此$x\left( {{b_1}} \right)$值即为当前通道与标准通道的时间偏差。

Figure 8.  Cross-correlation sequence of each channel and standard channel (channel 14) & Gaussian fitting results. (a) Cross-correlation sequence of each channel and standard channel & Gaussian fitting results, (b) cross-correlation sequence of channel 13 and its maximum value

总数据拟合残差直方图统计如图9所示，修正前内符合精度1.883 m （见图9（a）），修正后内符合精度1.809 m （见图9（b）），激光脉宽为5.7 ns。

Figure 9.  Fitting residual histogram statistics before (a) and after (b) correction（all channels）

互相关分析方法适用于所有的测量目标（包括地靶），但当识别信号点数较少时，因激光测距数据分布特征不明显而导致该方法失效。

在实际应用中，结合两种方法进行修正。在中国科学院云南天文台4×4阵列激光测距试验平台中，已经实现数据自动处理，在测量完成时，在10 min内可以输出数据预处理结果，根据2.2和2.3节泊松统计分析和循环拟合滤波的结果判断是否可以通过互相关分析方法进行修正，对于无法修正或者修正结果异常的情况，会提示需要使用地靶测量的方法进行修正。

同时，对于比较重要的目标，可以事后分别使用两种修正方法进行修正，再进行进一步的分析。

修正通道间的时间偏差后，再次进行数据拟合（与2.3节方法一致），滤除少量野点，进一步提高数据内符合精度。2.4节修正通道间的时间偏差后内符合精度为1.809 m，再次进行数据拟合后内符合精度为1.794 m。

经过以上匹配、识别、拟合、修正等过程后，修正并记录系统延迟等信息，可按照相应格式^[6]需求输出高精度数据，用于空间碎片的精密定轨等科学应用。

该节主要选取2018-12-01和2019-04-15的激光测距数据进行处理。

部分测量数据处理结果如表1所示，测量中激光脉宽5.7 ns，数据预处理为自动处理，过程中无需人为干预。

从表1中可以看出通过通道间时间偏差修正能有效提高精度，表中修正前后内符合精度（RMS）平均提高约7.2 cm。

其中2019-04-15日10571-1数据预处理结果如图10所示；2018-12-1日hy2a数据预处理结果如图11所示。

Figure 10.  10571 data preprocessing results. (a) Original data, (b) signal identification results, (c) histogram statistics of fitting residual after correction

Figure 11.  hy2a data preprocessing results. (a) Original data, (b) signal identification results, (c) histogram statistics of fitting residual after correction

修正通道间的时间偏差后，再次进行数据拟合，结果如表2所示，内符合精度（RMS）平均提高约3.2 cm。

文中研究了基于阵列探测技术的激光测距数据预处理方法，该方法已经成功应用于中国科学院云南天文台4×4阵列激光测距试验平台试验中，能很好的修正通道间的时间偏差，还促进了4×4阵列激光测距试验平台激光测距数据自动处理的实现，提高了系统自动化程度。

在目前试验系统中，受激光脉宽影响，数据内符合（RMS）精度偏低，未来在基于阵列探测技术的高精度激光测距测量中，该数据预处理方法中通道间的时间偏差修正将发挥出更大的作用。