Star spot extraction based on optimal background estimation for star sensor anti stray light

作为飞行器姿态控制单元的核心设备之一，星敏感器通过敏感天球上的恒星输出三轴姿态，具有精度高、自主性强和可靠性高等优点，是飞行器自主航行期间最可靠的设备^[1]。在星敏感器工作期间，太阳、月光、地气光、空间碎片反射光及测试环境中的光等杂散光照射进星敏感器，会造成星点提取数目少和星点质心定位精度差等问题，影响姿态解算结果。为了解决杂散光干扰下星点提取效果差的问题，文中开展了去杂散光干扰的星点提取研究，提出了基于杂散光最优背景估计的星点提取方法。

目前，星敏感器在抗杂散光技术上主要采用两种方法：遮光罩和抗杂光算法^[2-3]。遮光罩是一种可吸收和反射杂散光的挡光设备，这种方法虽然在一定程度上减少了杂散光的影响，但是不能完全排除杂散光的干扰。抗杂光算法是一类从杂散光图像中提取星点的算法，具有使用灵活、可移植性强及可面向多种干扰等特点，是目前星敏感器中抗杂光干扰的研究热点^[4]。抗杂光星点提取方法可以分为两类：基于杂光背景估计的方法和基于星点分布的方法。基于杂光背景估计的方法通常先估计出杂光背景然后在残差图像上检测星点，典型的方法有基于背景滤波的方法和最大背景估计的方法^[5-6]。前一种方法使用星点区域的边缘像素来估计星点图像的背景，后一种方法利用目标点附近亮区域像素估计背景。从构造背景的过程上看这两种方法是将杂散光背景建模成一个平面，然后通过目标附近区域的像素确定平面的值，这个平面是等值面。在实际中，杂散光在遮光罩内反射会形成多种形态的面，包括但不限于平面。在使用等值面提取星点时会产生提取星点中包含一定的杂散光能量或者星点能量被过多去除的情况，这会引起星点质心估计误差增大的情况。基于星点在图像中表现出的“凸起”特点，研究者提出了利用灰度前向差分方法和数学形态学等基于星点形状的提取方法^[7-8]。前者利用了星点和背景在灰度值上的阶跃特征，参考文献[7]中提出了基于“削顶”特点的提取方法。后者利用星点边缘结构元素进行检测和提取星点，参考文献[8]提出了一种改进的top-hat变换来提取星点。由于杂散光背景的复杂性，同样的星点在不同的背景下可能会表现出不同的特点。另外，在噪声干扰下，低信噪比星点的形态特点会表现的不显著。这些会影响这种基于星点形状方法的提取效果。这类基于星点特征的提取方法更多关注杂散光干扰下得星点检测问题，对星点质心估计的准确研究较少。

针对杂散光干扰下星点提取问题，文中提出了一种基于最优背景估计的抗杂散光星点提取方法，和参考文献[5-6]不同的是文中将杂散光成像建模成曲面。文中首先根据杂散光在星敏感器中的成像特点，构建了杂散光的曲面成像模型，然后提出了利用邻域像素估计杂散光最优背景的方法和星点提取方法。最后构建了仿真和实物实验，对比和分析了算法的性能。

在星敏感器工作过程中，杂散光进入视场的过程如图1所示。恒星距离星敏感器非常远，其相对于星敏感器的张角有限，在星敏感器中成像成直径不超过一定值的光斑。太阳、大气中的月光等杂散光通常距离星敏感器近，张角大，进入星敏感器中会成像成直径较大的光斑。

Figure 1.  Diagram of star sensor in stray light

在杂散光干扰下，星敏感器的成像模型可表示为：

式中：${I_i}$、${S_i}$、${B_i}$和${N_i}$分别表示第$i$个像素上的测量值、恒星能量强度、杂散光能量强度和噪声强度。

杂散光图像中，星点质心可通过公式（2）获得：

式中：$S({x_i},{y_i}|{\mathbf{v}})$表示参数为${\mathbf{v}}$的星点描述函数；$B({x_i},{y_i}|{\mathbf{\theta }})$表示参数为${\mathbf{\theta }}$的杂散光斑描述函数。公式（2）表示对噪声的最大似然估计意义下对恒星和杂散光能量的最优估计。由此可以看出，准确确定杂散光能量强度有利于准确估计星点质心。考虑到杂散光来源种类繁多、在遮光罩内的光线反射会增加成像复杂性、星敏感器测量能力有限等情况，很难用一个描述函数来表示图像中的杂散光光斑。

在星敏感器成像过程中，相邻位置的背景光亮度具有强相关性，在位置横纵坐标及亮度构成的三维空间中通常呈光滑、平缓的曲面。基于这个特点，文中将杂散光成像光斑在局部区域内建模成曲面。

由泰勒公式可知，图像中某像素的信息可以由该像素附近的像素信息描述^[9]。对于图像中的像素$({x_0},{y_0})$存在有如下关系：

式中：$h = x - {x_0}$；$k = y - {y_0}$；$ o\left( {{{(h + k)}^2}} \right) $表示无穷小量。在此忽略，公式（3）重写为公式（4）：

式中：$ a = {B_{xx}} $；$ b = {B_{yy}} $；$ c = 2{B_{xy}} $； $ d = {B_x} - 2{x_0}{B_{xx}} - 2{y_0}{B_{xy}} $；$ e = {B_y} - 2{y_0}{B_{yy}} - 2{x_0}{B_{xy}} $；$f = B - {x_0}{B_x} - {y_0}{B_y} + x_0^2{B_{xx}} + 2{x_0} {y_0}{B_{xy}} + y_0^2{B_{yy}}$；${B_{xx}}{\text{ = }}\dfrac{{{\partial ^2}B({x_0},{y_0})}}{{\partial {x^2}}}； {B_{yy}}{\text{ = }}\dfrac{{{\partial ^2}B({x_0},{y_0})}}{{\partial {y^2}}}； {B_{xy}}{\text{ = }} \dfrac{{{\partial ^2}B({x_0},{y_0})}}{{\partial x\partial y}}$；${B_x} = \dfrac{{\partial B({x_0},{y_0})}}{{\partial x}}； {B_y} = \dfrac{{\partial B({x_0},{y_0})}}{{\partial y}}$。

公式（4）即为文中定义的杂散光曲面模型。在星点图像的背景区域，恒星能量为0，那么，公式（2）可重写为公式（5），表示在对噪声的最大似然估计意义下对杂散光曲面的最优估计。

式中：$ {\mathbf{H}}({x_i},{y_i}|{\mathbf{\theta }}){\text{ = }}\sum\limits_i {{{\left[ {{I_i} - B({x_i},{y_i}|{\mathbf{\theta }})} \right]}^2}} $；参数${\mathbf{\theta }}{\text{ = (}}a,b,c, d,e,f)$。

公式（5）是一个线性最小二乘问题，可通过对$ {\mathbf{H}}(x,y|{\mathbf{\theta }}) $求偏导获得参数的闭合形式的最优解^[10]，如公式（6）所示。得到的参数${\mathbf{\hat \theta }}$代入公式（4）即可得到杂散光的曲面模型。

文中构建的杂散光干扰下星点检测方法的基本思路是首先利用星点图像的背景区域像素估计出杂散光在像素上的能量强度，然后通过图像减运算得到去杂散光的图像，最后在此图像中提取星点。算法分为两个过程：杂散光背景估计和星点提取过程。计算流程图如图2所示。

Figure 2.  Process of extracting star spot based on stray light background estimation

根据第1节描述的杂散光曲面模型，杂散光在某像素上的能量强度估计过程可描述为：首先确定像素的邻域，然后根据邻域像素估计杂散光曲面描述函数的参数，最后根据曲面参数估计杂散光在该像素上的能量强度。

某像素邻域为以该像素为中心，取一定长度的矩形区域内的像素作为该像素的邻域，表示为：

式中：$({x_0},{y_0})$表示该像素的位置；$2 l + 1$为矩形区域的长度；$l$为算法参数。

考虑到计算的方便性，采用图像滤波的计算方式，滤波矩阵可表示为：

相对于杂散光来说，恒星入射的能量是干扰因素。由于星点为直径不超过一定值的光斑，且小于杂散光形成的光斑区域。在公式（8）基础上，对滤波矩阵做进一步的改变，增加了掩膜，可表示为：

式中：掩膜长度$r$为算法参数。

图3所示为文中构建的一个滤波矩阵，其中滤波矩阵的长度为9，掩膜长度为3。

Figure 3.  An example of a constructed filter matrix

构建的这种滤波矩阵的物理意义为假定当前像素属于星点，取星点周边的背景区域像素来估计当前像素的值。这种方式可排除星点能量对估计杂散光能量的干扰。

在估计杂散光曲面过程时，由于噪声的影响，可能会得到存在较大误差的曲面。为了避免这种情况，文中使用最小均方误差来确定最优曲面模型。

对于像素$({x_0},{y_0})$，该方法可得到${(2 l + 1)^2}$个曲面，每个曲面在该像素位置上的均方根误差为${E_j}$，其中最小均方误差所对应的曲面即为该方法确定的在该像素上的曲面。

式中：${E_j} = E\left\{ {\left| {{I_i} - B({x_i},{y_i}|{{\mathbf{\theta }}_j})} \right|_2^2} \right\}$；$j$表示第$j$个曲面。

图4展示了一个位于图像背景变化剧烈附近位置的杂散光曲面筛选过程。图4（a）为确定的在该像素上的一个曲面，由于背景变化剧烈，存在确定的曲面不能准确描述背景的可能性。在此，文中通过对比残差图像的均方误差，选择具有最小的均方误差的曲面作为该位置上的最优曲面，如图4（b）所示。相比之下，这种方法得到的曲面能更准确描述杂散光能量。

Figure 4.  Determine the optimal surface of stray light. (a) One of the estimated surfaces; (b) The optimal surface

在获得杂散光曲面参数后，通过公式（11）即可获得杂散光在某像素上的能量强度，则：

式中：$({x_M},{y_M})$表示像素在公式（10）中所确定的曲面中相对位置。

由于星点可能位于图像中的任何位置，该方法要估计图像中每一个像素上的杂散光能量强度，这构成了图像中的杂散光背景。

在估计出杂散光背景后，通过图像减运算得到去杂散光的图像，然后采用快速高斯拟合法来提取星点质心，计算流程图如图5所示。

Figure 5.  The process of calculating the centroid of a star spot

快速高斯拟合法是高斯拟合法的一种改进方法，具有高精度和高计算效率的特点，其星点质心计算精度接近高斯拟合法，高于重心法和加权重心法，计算效率接近质心法，可满足星敏感器实时性需求^[11]。另外，采用快速高斯拟合法可构建基于信噪比的筛选像素方法，这有利于获得更多更准确的星点像素。

星点提取过程可分为3个过程，详述如下：

（1） 在去杂散光图像中，利用阈值法检测星点，得到星点的感兴趣区域。

去杂散光的图像可用公式（12）表示：

计算过程采用的阈值由公式（13）获得：

（2）在感兴趣区域内，首先选择那些高亮度的像素并利用快速高斯拟合法估计星点的描述函数参数，然后根据描述函数参数反演星点能量并计算像素信噪比。

星点描述函数定义如公式（14）所示：

式中：${\mathbf{v}} = (A,{\sigma _x},{\sigma _y},{x_c},{y_c})$，$({x_c},{y_c})$即为星点质心。

根据公式（2）、公式（12）和公式（14），快速高斯拟合法的目标函数可以表示为：

在公式（15）中，对数运算可将高斯拟合法中的非线性优化过程转换成线性优化过程，提高了计算效率。在式中添加权重$W$可使快速高斯拟合法的计算结果接近高斯拟合法的计算结果。在不牺牲计算精度的同时，快速高斯拟合法提高了计算效率。根据参考文献[11]的描述，权重可取${W_i} = {G_i}$。对目标函数公式（15）求偏导即可获得星点质心。

像素信噪比可以表示为：

（3）在获得包含星点的一定区域内的像素信噪比后，根据信噪比来筛选属于星点的像素，然后再次使用快速高斯拟合法获得星点质心。相比上一过程，此次参与计算星点质心的像素更多，得到的结果会更准确。在此次计算筛选像素过程中，信噪比阈值可设为1。

文中分别构建了仿真和实物实验来验证算法的性能，并与现有方法进行了对比。

基于最大背景估计的方法是基于背景滤波方法的改进，这里作为背景法的代表与文中方法对比，标记为“TTE”；相比于利用灰度前向差分技术的星点提取方法，基于形态学方法采用的星点形状特征更全面，面向的星点更多，这种方法在文中作为基于星点形状特征的提取方法的代表，标记为“NSTS”；文中提出的方法标记为“OSFE”。

在文中方法中有两个参数：决定滤波矩阵长度的$l$和决定掩膜长度的$r$。这两个参数会影响杂散光曲面模型估计的准确性，也会影响算法的计算效率。

掩膜长度由星点的尺寸决定。在星敏感器中，星点尺寸通常不超过$7 \times 7$ pixel。以此作为约束，文中$r$取为3个像素。

滤波矩阵长度决定了像素的邻域范围。对于估计杂散光的曲面模型来说，$l$越大估计的曲面参数会越准确，但是$l$太大，参与计算的像素会变多，会影响算法计算效率。此外邻域过大，杂散光的变化趋势可能会发生改变，会引入额外的误差。考虑到滤波矩阵与星点直径之间的关系，文中取约3倍星点直径作为滤波矩阵的长度，即$l$为10个像素。这样，一方面可保证选择的邻域中有足够多的背景像素，另一方面又不会引入过多的像素参与计算。

在仿真实验中，仿真图像由不同参数下生成的杂散光背景图像和星点光斑叠加合成。

考虑到很难有准确的杂散光光斑描述函数，文中采用数据插值方法生成杂散光背景图像。在图像中随机选择$n$个像素，像素灰度值随机赋值，然后根据这些像素的灰度值利用随机分布数据的方式插值其他像素的灰度值，得到图像即为杂散光背景图像。在此次实验中，图像包含$1\;024 \times 1\;024$ pixel，$n$为1~1024之间的随机数，一个生成的杂散光背景如图6（a）所示。

Figure 6.  Star images with different stray light. (a) Star image with the simulated stray light; (b) Star sensor image under the condition of moonlight; (c) Star sensor image under the condition of an artificial light

星点光斑图像由公式（17）生成，可表示为：

式中：$ A $表示星点的亮度水平；$\sigma $表示星点在星敏感器中的离散程度；$ （{x}_{c},{y}_{c}） $表示添加到图像中的星点位置。

在此次实验中，星点的亮度水平$A$为0.2~0.5之间的随机值。星点的离散程度$\sigma $为0.6~1.1之间的随机值。星点的位置在图像内随机选择。图像中添加独立同分布的高斯噪声，噪声的强度依次从0增大到0.1，并且每次噪声强度变化产生1000幅图像，在此次实验中每幅图像中添加400个星点。实验结果如图7和8所示。

Figure 7.  Star spots detection rate

Figure 8.  Star spots false detection rate

检测率定义为提取的真星和所有真星之间的比值，误检率定义为提取到的假星与所有提取的星点之间的比值。随着添加的噪声强度增大，背景法确定的阈值会偏大，这使得一些低信噪比的星点不易被检测到，从而导致星点检测率下降。当噪声强度增大时，星点形状特征信息会变得模糊，特别是低信噪比的星点，这使得基于形态学方法检测到的真实星点少，星点检测率降低。对于位于杂散光背景变化显著区域的像素来说，这些像素的灰度值与周围区域的像素存在一定的差异，随着噪声强度增大，这些差异会变的更加显著。如果多个这样的像素聚集在一起，就形成了一个和星点具有类似特征的“星点”——假星，这是算法提取到假星的原因。由于背景法和基于形态学的方法无法识别这类假星，导致星点误检率随噪声强度增大而上升。文中的方法在一定程度上对这类假星不敏感，误检率较低。

在提取能力实验的基础上，取那些提取出来的星点，计算星点的质心误差，如公式（18）所示，并使用星点的平均质心误差分析算法性能。结果如图9所示。

Figure 9.  Centroid error of star spots extracted

在杂散光图像中，影响星点质心精度的因素主要有两种：星点图像中的噪声和杂散光。噪声强度增大会引起星点质心误差增大^[11]。在此次实验中，随着添加的噪声强度增大，同一算法得到的质心误差曲线会随着增大。由于杂散光能量的存在，公式（2）所描述的星点质心估计过程不再是无偏估计，质心估计结果中会存在一定的估计偏差。基于形态学的方法对杂散光能量的估计较弱，相应的估计结果中存在较大的估计偏差。背景法相当于将杂散光能量建模成一个等值的平面，而文中是将杂散光能量建模成不等值的曲面。相比之下，文中对杂散光能量估计的更准确些，质心结果中存在较小的估计偏差，星点质心误差会更小些。

该实验利用月光、人工光源等构建了杂散光实验环境，星敏感器在此环境中采集图像，统计和对比各个算法的提取结果。

在月光和人工光源干扰下得到的星点图像如图5（b）和（c）所示，取杂散光实验环境中获得的10幅图像。在这10幅图像中，各个算法提取的星点统计结果如表1所示。

在月光干扰实验中，各个算法得到的星点检测个数和误检个数相近，主要原因是月光形成的杂散光背景比较平缓，星点的信噪比较高，这有利于星点的检测。而在仿真图像中，添加了高强度的噪声，这增加了星点检测的难度。

相比于月光干扰，人工光源干扰下得到的杂光背景要更复杂，这使得提取结果变差。从检测到的星点数目和误检数目上看，在人工干扰实验中，文中方法的提取结果要优于背景法和基于形态学方法提取到的结果。

针对星敏感器中杂散光干扰下提取到的星点数目少和质心精度差等问题，文中提出了一种基于杂散光最优背景估计的抗杂散光星点提取方法。该方法分析和构建了杂散光的曲面模型，并设计了估计杂散光最优背景的方法和星点提取方法。经过仿真实验和实物实验测试，在检测率、误检率和质心精度上，该方法的表现优于现有的基于背景的提取法和基于形态学的方法。这显示该方法具有较好的抗杂散光能力。在实际中，杂散光的成像特点多种多样，后续将构建更多测试环境对算法做进一步的测试和分析。