CN106646430B - 一种基于地面探测器的激光足印中心确定方法 - Google Patents

Abstract

本发明公开了一种基于地面探测器的足印中心确定方法，所述方法包括以下步骤：S1：对地面探测器捕获的激光足印进行预处理；S2：利用高斯曲面拟合法确定激光足印中心。本发明通过对激光足印中心中的孤立点、缺失点和异常点进行预处理，减小了由于外界环境及探测器不一致造成的噪声误差，使激光足印阵列恢复均匀平滑形状；利用足印范围内覆盖85％的能量进行高斯曲面拟合，求取曲面的中心位置，即为激光足印的中心，以此作为激光测高仪在轨几何检校的地面参考。

Description

一种基于地面探测器的激光足印中心确定方法

技术领域

本发明涉及激光测高仪在轨检校技术领域，特别涉及一种基于地面探测器的激光足印中心确定方法。

背景技术

激光测高因具有方向性好、测距精度高等特点，在深空探测和地球科学领域中体现了巨大的应用潜力，将星载激光测高技术应用于高分辨率测绘卫星，辅助光学相机等载荷提高高程精度是一种新颖的手段。激光测高仪上包含激光发射器和激光接收器，各部件按照设计的安装角安装在卫星上，但各角在卫星在发射过程中受到机械的震动及在运行过程中复杂的太空环境会发生变化，致使激光发射器与激光发射器与激光测高仪、卫星平台之间的实际安装角与设计值有微小偏差，使得激光打在地面上的水平位置不确定。因此，高精度在轨几何检校是激光数据应用的前提。

激光以一定的发散角从卫星平台发射到地面并非一个点，而是一个有一定半径的圆形光斑，称为激光足印。对于基于地面探测器的星载激光测高仪在轨几何检校而言，激光足印中心(即，激光发射器指向地面点的确切位置)的确定是在轨几何检校的第一步，其精度高低直接影响到检校的精度。

目前，激光足印中心确定的方法不多，从公开的资料来看，只有美国GLAS卫星通过地面探测器接收卫星过境时的激光足印的实验取得了成功，实现了在轨几何检校。其激光发射器是以40HZ的频率在地球表面“照”出一个直径为70m的足印光斑，激光能量探测器只有0和1两个档位，“0”代表没有响应到激光能量，“1”代表响应到激光能量。通过求取显示为“1”的探测器能量阵列的形心来作为足印中心，得到中心点的地面坐标，实现发射器安装角的反算。

2016年5月30日，我国首颗载有激光测高仪的卫星——ZY-302(资源三号02星)在太原卫星发射中心成功发射升空，其激光发射器是以2HZ的频率发射激光脉冲，试验人员利用地面探测器成功捕获激光足印光斑，该探测器分为三个档位，每个档位采用8个能量等级，可以准确确定激光足印中心，顺利地完成了在轨几何检校。

首先，激光在大气中传播经过两次衰减才被星上激光接收器接收，而地面探测器捕获的足印是经过1次衰减得到；其次，激光从卫星平台发射至地面，还会受到风、云、雾、空中颗粒等的影响发生吸收、前向散射等现象；第三，由于探测器的一致性问题，使得地面探测器接收到的激光能量等级分布并不像实验室里呈现类高斯分布的特性。因此如何准确确定激光足印的中心位置，是本发明主要要解决的技术问题。

发明内容

针对基于地面探测器的激光足印中心确定需求，本发明提出了一种基于地面探测器的激光足印中心确定方法，可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。

本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明目的和优点。

本发明提供了一种基于地面探测器的足印中心确定方法，其特征在于，所述方法包括以下步骤：

S1：对地面探测器捕获的激光足印进行预处理；

S2：利用高斯曲面拟合法确定激光足印中心；其中，

步骤S1具体包括以下子步骤：

S11：计算激光足印在地面上的半径r；

S12：检测所述激光足印内的孤立点、缺失点和异常点；

S13：对检测出来的孤立点、缺失点和异常点进行校正；

步骤S2具体包括以下子步骤：

S21：提取所述激光足印内能量等级近似服从高斯分布的数据；

S22：基于步骤S21提取的数据，计算激光足印中心点的坐标。

优选的，所述步骤S11具体为：

根据下式计算激光足印在地面上的半径r：

其中，r为激光足印在地面上的半径，R₁为足印中心与星下点之间的水平距离，R₂为星下点与足印边缘的最小距离，H为卫星轨道高度，θ为激光发射器的发射角，δ为激光能量发散角。

优选的，所述步骤S12具体包括：

S121：在所述激光足印中选定一个3×3的窗口，设所述3×3的窗口的中心点(也称为窗口中心元素)的坐标为I₅(x₅,y₅)，足印中心元素的坐标为P(x_p,y_p)，集合X为所述3×3的窗口内所有元素的集合，集合Y为除了窗口中心元素外的邻域内所有元素的集合，优选的，集合X和Y可表示为：

X＝{I_i},i＝1,2...9，Y＝{I_i}，i＝1,2...9,i≠5

S122：计算窗口中心元素I₅(x₅,y₅)和足印中心元素P(x_p,y_p)之间的距离d：

S123：利用窗口中心元素I₅(x₅,y₅)和足印中心元素P(x_p,y_p)的距离d、足印半径r和窗口中心元素的能量等级I₅来确定所述窗口中心元素是否属于孤立点、缺失点或异常点。

优选的，所述步骤S123具体包括：

(i)，如果d＞r，I₅≠0，并且集合Y中的元素为0的个数大于等于第一阈值，则判断所述窗口中心元素I₅(x₅,y₅)为孤立点；优选的，第一阈值可以为7或者8；

(ii)，如果d＜r，I₅＝0，并且集合Y中的元素不为0的个数大于第二阈值，则判断所述窗口中心元素为缺失点；优选的，第二阈值可以为7或8；

(iii)，如果d＜r，I₅≠0，|I₅-M|＞ε，并且集合Y中的元素不为0的个数大于等于第三阈值，则判断所述窗口中心元素为异常点，其中，M为点I₅(x₅,y₅)的不含有误差的能量等级，ε为噪声误差，则判断所述窗口中心元素为异常点；优选的，第三阈值可以为7或8。

优选的，在“判断所述窗口中心元素为异常点”的步骤中，M为集合X中的所有元素的平均值mean(X)，ε为集合X中的所有元素的标准差std(X)，其中，mean(X)和std(X)的计算公式如下：

优选的，步骤S13具体包括：

(i)，对于孤立点，将孤立点删除；

(ii)，对于缺失点，采用邻域点均值内插的方法进行补偿，具体为，根据邻域内的元素内插出窗口中心元素I₅(x₅,y₅)的能量等级I₅：

I₅＝(∑I_i)/8，I_i∈Y

其中，8为集合Y的元素个数；

(iii)，对于异常点，也采用邻域点均值内插的方法进行补偿，具体为，根据邻域内的元素内插出窗口中心元素I₅(x₅,y₅)的能量等级I₅：

I₅＝(∑I_i)/8，I_i∈Y

其中，8为集合Y的元素个数。

优选的，所述步骤S22具体包括以下子步骤：

S221，根据下式(1)计算激光足印中心点的坐标(x_u,y_u)：

其中，A为高斯曲面的峰值，x_i,y_i表示步骤S21提取的第i个数据的坐标，P_i为坐标x_i,y_i点的能量等级，σ_x为高斯曲面在X方向上的标准差，σ_y为高斯曲面在Y方向上的标准差，x_μ、y_μ为待求的激光足印中心点的坐标；

S222，对上式(1)两边取自然对数：

令

将式(3)代入式(2)可得：

当有m个数据参与曲面拟合时，写成矩阵的形式为：

简记为：

Z＝XB (6)

S223，根据最小二乘原理，选择单位权阵P，将式(6)改写为：

(X^TPX)B＝(X^TP)XB＝(X^TP)Z＝X^TPZ (7)

由式(7)可得，矩阵B的广义最小二乘解为：

B＝(X^TPX)^-1X^TPZ (8)

S224，将根据式(8)得到的B的值代入到式(3)，得到激光足印中心的坐标x_μ、y_μ。

本发明还提供了一种计算机程序，当执行所述计算机程序时，执行上述任一种方法。

本发明通过对激光足印中心中的孤立点、缺失点和异常点进行预处理，减小了由于外界环境及探测器不一致造成的噪声误差，使激光足印阵列恢复均匀平滑形状；利用足印范围内覆盖85％的能量进行高斯曲面拟合，求取曲面的中心位置，即为激光足印的中心，以此作为激光测高仪在轨几何检校的地面参考。

附图说明

图1为根据本发明实施例的、基于地面探测器的激光足印中心确定方法的流程图；

图2为根据本发明实施例的、计算激光足印在地面上的半径的示意图；

图3A-3C分别为预处理前的、地面探测器阵列能量等级、平面示意图和三维示意图；

图4为根据本发明实施例的、孤立点、缺失点或异常点检测示意图；

图5A-5C分别示出了孤立点、缺失点或异常点的示意图；

图6A-6C分别为预处理后的、地面探测器阵列能量等级、平面示意图和三维示意图；

图7示出了小于86.5％能量等级的最大区间的区域及中心点位置；

图8示出了高斯曲面拟合后的地面探测器阵列值、平面示意图和三维示意图；

图9示出了高斯曲面拟合前后各点叠加结果的示意图。

具体实施方式

激光测高是通过量测激光出光时刻到接收到激光的时间间隔来准确测算距离的，但由于激光在到达地面的途中，受到风、云、雾、空中颗粒等的影响会发生吸收、前向散射等一系列现象，并且由于地面探测器的一致性问题，使得地面上探测器阵列显示的能量等级不符合地面实验室里的高斯分布特性，因此，要得到足印的中心，首先需要对获得的探测器阵列中的能量等级进行预处理，即孤立点、缺失点、异常点检测和处理，再对足印范围内小于86.5％的能量等级的最大区间的数据进行高斯曲面拟合，以确定足印的中心位置，用于激光测高仪在轨几何检校。

下面对本发明进行更全面的描述，其中说明本发明的示例性实施例。

如图1所示，本发明提出的基于地面探测器的足印中心确定方法，包括以下步骤：S1：对地面探测器捕获的激光足印进行预处理；S2：利用高斯曲面拟合法确定足印中心。

下面结合其他附图对上述步骤S1和S2进行详细说明。

S1：对地面探测器捕获的激光足印进行预处理(即，对探测器阵列中的能量等级进行预处理)

步骤S1具体包括以下子步骤：

S11：计算激光足印在地面上的半径r；

激光发射器不可能完全垂直地安装在卫星平台上，且卫星在飞行过程中的姿态也时刻发生着变化，为了获得回波，激光发射器以发射角θ(包含安装角和姿态角)和发散角δ发射激光，地面探测器在地面上捕获的激光并非一个点，而是一个偏离星下点且半径为r的光斑，称为激光足印，在本发明中，激光足印由探测器阵列中的能量等级来表示，即，由对应于探测器阵列中的能量等级值的各个点来表示激光足印。由于在轨几何检校时，为了排除一些干扰因素对于探测器阵列能量等级的影响，以计算所得的激光足印半径来约束一些较远的无效点或者误亮点，即后续孤立点、缺失点、异常点的约束条件。

如图2所示，根据下式计算激光足印在地面上的半径r：

其中，r为激光足印在地面上的半径，R₁为足印中心与星下点之间的水平距离，R₂为星下点与足印边缘的最小距离，H为卫星轨道高度，θ为激光发射器的发射角，δ为激光能量发散角。由于H、δ为设计值，都是已知量，θ为激光发射器与卫星平台的实际安装角与卫星飞行过程中姿态角之和。检校前以设计值为安装角，飞行中得姿态角可以通过陀螺和星敏获得，因此为已知量。所以根据上述公式可以计算出足印半径r。

S12：检测激光足印内的孤立点、缺失点和异常点；

按照一定的间距在地面上铺设探测器，当卫星过境时，激光发射器发射激光脉冲，地面探测器捕获激光，将光信号转换为电信号，再根据电流的大小设置若干等级，激光能量高的电流大，能量等级越高，亮起的灯数目越多，激光能量低的电流小，能量等级越低，亮起的灯数目越少。将受激光而亮灯的探测器构成的矩阵称为探测器阵列，每个探测器亮灯的数量称为该点的能量等级。由于实验室中86.5％的激光能量的服从类高斯分布，呈现中间高，边缘低的趋势，因此地面探测器获得的能量等级也应该服从类高斯分布，从中心往外呈现递减的趋势。但阵列周围有个别未复位或由于地面探测器一致性等原因引起的孤立点，还有一些由于地面探测器一致性或者大气等原因引起的缺失点、异常点，使得能量等级形成的曲面复杂多样，无法通过直接拟合的方法求取足印中心。下面对孤立点、缺失点和异常点作出具体说明。

孤立点：位于实际足印有效范围之外，由于地面探测器未复位或受探测器一致性、大气等影响而误亮的点。

缺失点：位于实际足印有效范围之内，距离足印中心比较近，由于地面探测器一致性或者大气等原因而缺省的点。

异常点：位于实际足印有效范围之内，能量等级明显低于或高于周围值的点，这类点也可能是由于未复位、地面探测器一致性、大气等原因引起的。

下面结合图3A，对孤立点、缺失点和异常点作进一步的解释说明。如图3A所示，其中示出了根据本发明实施例的、地面探测器阵列能量等级示意图。图3A中第一行中的8和3因为距离足印中心点较远，并且周围没有探测器响应，当属于孤立点；足印内部两个0点位于足印内部，且周围点都有值，因此是缺失点；由于边缘能量等级为8的点虽位于足印内部，但并不位于足印形心处，其值明显高于周围的点，因此可以确定为异常点。

步骤S12具体包括：

S121：如图4所示，在所述激光足印中选定一个3×3的窗口，设所述3×3的窗口的中心点(下文称为窗口中心元素)的坐标为I₅(x₅,y₅)，足印中心元素的坐标为P(x_p,y_p)，集合X为所述3×3的窗口内所有元素的集合，集合Y为除了窗口中心元素外的邻域内所有元素的集合(即，集合Y为所述3×3的窗口内除了中心元素外的其他所有元素的集合)，集合X和Y可表示为：

X＝{I_i},i＝1,2...9，Y＝{I_i}，i＝1,2...9,i≠5

S122：计算窗口中心元素I₅(x₅,y₅)和足印中心元素P(x_p,y_p)之间的距离d：

S123：利用窗口中心元素I₅(x₅,y₅)和足印中心元素P(x_p,y_p)的距离d、足印半径r和窗口中心元素的能量等级I₅来确定所述窗口中心元素是否属于孤立点、缺失点或异常点。具体包括：

(i)，如果d＞r，I₅≠0，并且集合Y中的元素为0的个数大于等于第一阈值，则判断所述窗口中心元素I₅(x₅,y₅)为孤立点；优选的，第一阈值可以为7或者8(优选的，第一阈值为8，因为在整个检校场中，在没有激光打下来的情况下，几百个探测器仅有1-2个会发生响应，因此若无激光数据，激光探测器发生响应的概率很小，因此，孤立点周围的点应当属于没有接收到激光，即能量等级为0，因此第一阈值定为8)，即，集合Y中的元素几乎为0或者全部为0；

孤立点，主要是由于未复位或者探测器一致性的原因而引起的误亮的点，在曲面拟合中易造成多值或无解的情况。如图5A所示，窗口中心点能量等级为8，而周围能量等级为0，中心元素I₅(x₅,y₅)与足印中心元素P(x_p,y_p)点的距离d大于足印半径r，那么此点便为孤立点。

(ii)，如果d＜r，I₅＝0，并且集合Y中的元素不为0的个数大于第二阈值，则判断所述窗口中心元素为缺失点；优选的，第二阈值可以为7或8，(优选的，第二阈值为7，理论上，激光的能量经过大气层等会发生衰减，但其能量等级分布与实验室中内高外地的特性不会发生太大变化，然而从激光在轨检校捕获的多个足印阵列来看，在足印内部，并不是所有的探测器都有能量等级，会出现有多个点能量等级缺失，激光探测器接收激光却不响应的概率高于没有激光而激光探测器响应的概率，但这种探测器被放置在同一区域的概率又很低，因此，第二阈值设置为7。因此，缺失点周围的点大部分不为0或者全部不为0)即，集合Y中的元素大部分不为0或者全部不为0；

在激光足印内，窗口中心点能量等级为0，而周围的值大部分不为0，以此来确定其为缺失点。如图5B中能量等级为0的元素I₅(x₅,y₅)。

(iii)，如果d＜r，I₅≠0，|I₅-M|＞ε，并且集合Y中的元素不为0的个数大于等于第三阈值，则判断所述窗口中心元素为异常点，其中，M为点I₅(x₅,y₅)的不含有误差的能量等级，ε为噪声误差，则判断所述窗口中心元素为异常点；优选的，第三阈值可以为7或8(这里暂定为8，从激光在轨检校捕获的多个足印阵列来看，在足印内部，会出现少量地面探测器获取的能量等级有异常，异常点周围的能量等级形成的曲面近似光滑，大部分不为0或者全部不为0。可以认为都是正常点，因此，第三阈值设置为8。)，即，集合Y中的元素大部分不为0或者全部不为0。优选的，该步骤中还可以规定I₅(x₅,y₅)不是距离所述足印中心最近的能量等级最高的点，用以排除将探测器阵列中距离足印中心最近的最大值点为异常点的情况。

在激光足印内部异常点会高于或低于阈值而使得曲面不平滑，一般为高于周围点，有多个峰值，会造成曲面拟合多解或无解的情况。优选的，M为集合X中的所有元素的平均值mean(X)，ε为标准差std(X)，若|I₅-mean(X)|＞std(X)，则I₅(x₅,y₅)点被判定为异常点。其中，mean(X)和std(X)的计算公式如下：

如图5C中，中心点I₅(x₅,y₅)虽位于激光足印内部，但不是靠近足印中心且能量等级最高的点，由于

满足|I₅-M|＞ε，因此确定中心点I₅(x₅,y₅)为异常点，将不参与中心点校正的计算。

S13：对检测出来的孤立点、缺失点和异常点进行校正；

对孤立点剔除时直接赋值为0，缺失点补偿和异常点校正主要有双线性内插法、曲面拟合法、邻域点均值内插法，都顾及了周围点的能量等级。由于能量等级本身就带有一定截断误差，所以本发明采用邻域点均值内插法来补偿缺失点、校正异常点。

步骤S13具体包括：

(i)，对于孤立点，将孤立点删除；具体的，通过将孤立点赋值为0来删除孤立点，如图5A中的孤立点直接赋值为0；

(ii)，对于缺失点，采用邻域点均值内插的方法进行补偿，具体为，根据邻域内的元素内插出窗口中心元素I₅(x₅,y₅)的能量等级I₅：

I₅＝(∑I_i)/8，I_i∈Y

其中，8为集合Y的元素个数。

如图5B示出的例子，计算的缺失点I₅(x₅,y₅)的能量等级I₅为：I₅＝(1+2+3+1+4+1+3+5)/8＝2.5。

(iii)，对于异常点，也采用邻域点均值内插的方法进行校正值I₅，即如上式所示。具体的，对于异常点，采用邻域点均值内插的方法进行补偿，具体为，根据邻域内的元素内插出窗口中心元素I₅(x₅,y₅)的能量等级I₅：

I₅＝(∑I_i)/8，I_i∈Y

其中，8为集合Y的元素个数。

如图5C所示的例子，I₅＝mean(Y)＝(4+2+1+5+1+2+1+0)/8＝16/8＝2。

滑动窗口，以第一行第一列作为中心点向右向下滑动开始进行运算，判断所有的点的属性，最终预处理后的结果如图6A所示，共有2个孤立点，2个缺失点，1个异常点，其余的都是正常点，不进行能量等级的校正和更改。

由于探测器阵列间距比较大，在判定点的属性时选定的窗口不宜选择过大，因此，本发明的实施例选定一个3×3的窗口，当然其他窗口也可以根据实际需要来选择，本发明对此不作具体限定。

图6A-6C示出了预处理后的、地面探测器阵列的能量等级、平面示意图和三维示意图，与图3A-3C相比，可以看出，地面探测器阵列能量等级经过孤立点剔除、缺失点补偿和异常点校正之后，无论从平面视图还是从三维视图来看，能量等级曲面都基本恢复了平滑的形状，呈现内高外地的趋势，将噪声有效地控制在小范围内。

在执行了上述步骤S1之后，下面详细描述步骤S2：利用高斯曲面拟合法确定足印中心。

步骤S2具体包括以下子步骤：

S21：提取激光足印内能量等级近似服从高斯分布的数据；

发明人经过大量试验研究发现，激光足印86.5％的能量近似服从高斯分布，因此本发明优先选择探测器阵列内86.5％的能量等级作为限制，截取86.5％的能量等级中最大等级区间作为探测器阵列中心提取的数据，利用高斯曲面拟合法来实现激光足印中心的提取。

S22：基于步骤S21提取的数据，计算激光足印中心点的坐标。

步骤S22具体包括以下子步骤：

S221，根据下式(1)计算激光足印中心点的坐标(x_u,y_u)：

由于激光足印内每一个点P_i(x_i,y_i)的能量等级P_i近似服从以P_u(x_u,y_u)的高斯分布，所以可以表达成如式(1)所示的形式。

在该实施例中，假设激光能量等级符合以P_u(x_μ,y_μ)点为中心的高斯分布，上式中，A为高斯曲面的峰值，x_i,y_i表示步骤S21提取的第i个数据的坐标，P_i为坐标x_i,y_i点的能量等级(该能量等级为已知值)，σ_x、σ_y为高斯曲面的半宽信息，具体的，σ_x为高斯曲面在X方向上的标准差，σ_y为高斯曲面在Y方向上的标准差，σ_x或σ_y越大，则数据在x或y方向上的分布越分散，反之，数据分布就越集中，x_μ、y_μ为待求的激光足印中心点的坐标。

上式中，待估参数A为高斯曲面的峰值，σ_x、σ_y为高斯曲面的半宽信息，是在保证回波同时不损失高程精度的情况下的设计值，但由于受到大气各种因素的影响会有衰减，并且通过能量等级进行取舍换算，是一个可以通过曲面拟合计算出来的值，x_μ、y_μ为高斯曲面峰值点的坐标，即地面探测器阵列能量等级中心点坐标(待求的激光足印中心点的坐标)，通过曲面拟合求得。i＝1,2，......m，表示共有m个点参与曲面拟合计算(m为大于等于2的正整数)。

S222，对上式(1)两边取自然对数：

令

将式(3)代入式(2)可得：

当有m个数据参与曲面拟合时，写成矩阵的形式为：

简记为：

Z＝XB (6)

S223，根据最小二乘原理，选择单位权阵P，将式(6)改写为：

(X^TPX)B＝(X^TP)XB＝(X^TP)Z＝X^TPZ (7)

由式(7)可得，矩阵B的广义最小二乘解为：

B＝(X^TPX)^-1X^TPZ (8)

因为Z与各点的能量等级P_i有关，X与各点的坐标(x_i,y_i)有关，均为已知量，且P也是已知的单位权阵，因此，由公式(8)可计算出B。

S224，将根据式(8)得到的B的值代入到式(3)，得到激光足印中心的坐标x_μ、y_μ。

具体的，根据式(8)解算得到B，由于B＝[b₁ b₂ b₃ b₄ b₅]^T，所以得到b₁、b₂、b₃、b₄、b₅，将b₁、b₂、b₃、b₄、b₅代入式(3)，得到高斯曲面的5个参数：A、x_μ、y_μ、σ_x、σ_y。P为权阵，可取为单位阵。其中，x_μ、y_μ为探测器阵列的中心，即激光足印的中心位置，利用此点的地面点坐标用于激光测高仪在轨几何检校。

根据本发明的实施例，如图6A中，整个探测器阵列能量等级总和为81.5，则86.5％的能量等级总和为81.5×86.5％＝70.4975，而探测器阵列中能量等级大于等于2的点的能量等级总和为65.5，能量等级为1的点的总和为16，由于激光在实验室中呈现的是类高斯分布，用不同能量等级区间截断时，能量等级比较高的区域位置呈现震荡趋势，因此不宜选择过大的区域用于足印中心提取，并且激光在大气中传播经过两次衰减才被星上激光接收器接收，而地面探测器捕获的足印是经过1次衰减得到，因此在此舍弃能量等级为1的点，以能量等级大于等于2的点进行高斯曲面拟合计算探测器阵列中心点坐标。

将选定的各点的坐标(x_i,y_i)及其能量等级P_i代入式(4)中，得到矩阵X和矩阵Z，

利用最小二乘方法根据式(7)得(X^TPX)B＝X^TPZ为：

根据(8)，则B＝(X^TPX)^-1X^TPZ，在此，没有对X^TPX求逆，而是将(X^TPX)B＝X^TPZ写成增广矩阵G＝[X^TPX|X^TPZ]的形式，采用高斯消元法化简之后求得B的解，即：

也即：

于是

再将B代入式(3)中解出高斯曲面的5个参数分别为A＝5.48499，x_μ＝6.549994，y_μ＝4.73468，σ_x＝1.71984，σ_y＝1.57786。如图7所示，以行列号的值代表所求地面点的坐标，则中心点的坐标即图7的箭头所指点。

曲面拟合所得的中心点位于探测器阵列中能量等级为8的点的左上方，点的坐标为(6.55，4.73)，而在轨几何检校中若采用能量等级最高的8点作为足印中心，则中心点的坐标为(7，5)，两点之间的像素距离为：

若探测器阵列间隔为20m，则两点之间的平面距离为10.4m，假设卫星高度为500km，由于激光发射器得安装可近似认为垂直于地面，即认为θ＝0(包括垂轨和沿轨方向上的角度)，10.4m的水平距离引起的检校误差为：δ＝arctan(10.4/500000)×180×3600/π＝4.29"，说明了激光足印中心的精度对在轨几何检校精度的影响是很大的。

将高斯曲面的参数回代到式(1)中，离散化各点，保留两位小数，得到如图8平面视图。

如图8所示，拟合和和拟合前，能量等级分布都是上低下高的趋势，这是偏移所致，由内而外是内高外低的趋势。能量等级最高的8点在拟合后的离散结果为5.23，二者，相差最大，而其他点前后相差不大，也反过来验证了计算的正确性。从图9拟合前后能量等级的三维套合结果来看，二者只有最高点有较大误差，其余各点都套合紧密，为最小二乘拟合的最佳结果。

另外，本发明还提供了一种计算机程序，当执行所述计算机程序时，执行上述实施例所述的方法。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (4)

1.一种基于地面探测器的足印中心确定方法，激光以一定的发散角从卫星平台发射到地面并非一个点，而是一个有一定半径的圆形光斑，称为激光足印；其特征在于，所述方法包括以下步骤：

S1：对地面探测器捕获的安装在卫星平台上的激光发射器发射的激光足印进行预处理；

S2：利用高斯曲面拟合法，对足印范围内小于86.5％的能量等级的最大区间的数据进行高斯曲面拟合，确定激光足印中心；其中，

步骤S1具体包括以下子步骤：

S11：计算激光足印在地面上的半径r；

S12：检测所述激光足印内的孤立点、缺失点和异常点；

S13：对检测出来的孤立点、缺失点和异常点进行校正；

步骤S2具体包括以下子步骤：

S21：提取所述激光足印内能量等级近似服从高斯分布的数据；

S22：基于步骤S21提取的数据，计算激光足印中心点的坐标；

步骤S12具体包括：

S121：在所述激光足印中选定一个3×3的窗口，设所述3×3的窗口的中心点，也称为窗口中心元素的坐标为I₅(x₅,y₅)，足印中心元素的坐标为P(x_p,y_p)，集合X为所述3×3的窗口内所有元素的集合，集合Y为除了窗口中心元素外的邻域内所有元素的集合，集合X和Y表示为：

X＝{I_i},i＝1,2...9，Y＝{I_i}，i＝1,2...9,i≠5

S122：计算窗口中心元素I₅(x₅,y₅)和足印中心元素P(x_p,y_p)之间的距离d：

S123：利用窗口中心元素I₅(x₅,y₅)和足印中心元素P(x_p,y_p)的距离d、足印半径r和窗口中心元素的能量等级I₅来确定所述窗口中心元素是否属于孤立点、缺失点或异常点；

其中，所述步骤S123具体包括：

(i)，如果d＞r，I₅≠0，并且集合Y中的元素为0的个数大于等于第一阈值，则判断所述窗口中心元素I₅(x₅,y₅)为孤立点；第一阈值为7或者8；

(ii)，如果d＜r，I₅＝0，并且集合Y中的元素不为0的个数大于第二阈值，则判断所述窗口中心元素为缺失点；第二阈值为7或8；

(iii)，如果d＜r，I₅≠0，|I₅-M|＞ε，并且集合Y中的元素不为0的个数大于等于第三阈值，则判断所述窗口中心元素为异常点，其中，M为点I₅(x₅,y₅)的不含有误差的能量等级，ε为噪声误差，则判断所述窗口中心元素为异常点；第三阈值为7或8；

步骤S13具体包括：

(i)，对于孤立点，将孤立点删除；

(ii)，对于缺失点，采用邻域点均值内插的方法进行补偿，具体为，根据邻域内的元素内插出窗口中心元素I₅(x₅,y₅)的能量等级I₅：

I₅＝(∑I_i)/8，I_i∈Y

其中，8为集合Y的元素个数；

(iii)，对于异常点，也采用邻域点均值内插的方法进行补偿，具体为，根据邻域内的元素内插出窗口中心元素I₅(x₅,y₅)的能量等级I₅：

I₅＝(∑I_i)/8，I_i∈Y

其中，8为集合Y的元素个数。

2.根据权利要求1所述的方法，其特征在于，所述步骤S11具体为：

根据下式计算激光足印在地面上的半径r：

其中，r为激光足印在地面上的半径，R₁为足印中心与星下点之间的水平距离，R₂为星下点与足印边缘的最小距离，H为卫星轨道高度，θ为激光发射器的发射角，δ为激光能量发散角。

3.根据权利要求1所述的方法，其特征在于，在“判断所述窗口中心元素为异常点”的步骤中，M为集合X中的所有元素的平均值mean(X)，ε为集合X中的所有元素的标准差std(X)，其中，mean(X)和std(X)的计算公式如下：

4.根据权利要求1所述的方法，其特征在于，所述步骤S22具体包括以下子步骤：

S221，根据下式(1)计算激光足印中心点的坐标(x_u,y_u)：

其中，A为高斯曲面的峰值，x_i,y_i表示步骤S21提取的第i个数据的坐标，P_i为坐标x_i,y_i点的能量等级，σ_x为高斯曲面在X方向上的标准差，σ_y为高斯曲面在Y方向上的标准差，x_μ、y_μ为待求的激光足印中心点的坐标；

S222，对上式(1)两边取自然对数：

令

将式(3)代入式(2)可得：

当有m个数据参与曲面拟合时，写成矩阵的形式为：

简记为：

Z＝XB (6)

S223，根据最小二乘原理，选择单位权阵P，将式(6)改写为：

(X^TPX)B＝(X^TP)XB＝(X^TP)Z＝X^TPZ (7)

由式(7)可得，矩阵B的广义最小二乘解为：

B＝(X^TPX)^-1X^TPZ (8)

S224，将根据式(8)得到的B的值代入到式(3)，得到激光足印中心的坐标x_μ、y_μ。

Images