-
在星敏感器工作过程中,杂散光进入视场的过程如图1所示。恒星距离星敏感器非常远,其相对于星敏感器的张角有限,在星敏感器中成像成直径不超过一定值的光斑。太阳、大气中的月光等杂散光通常距离星敏感器近,张角大,进入星敏感器中会成像成直径较大的光斑。
在杂散光干扰下,星敏感器的成像模型可表示为:
$$ {I_i} = {S_i} + {B_i} + {N_i} $$ (1) 式中:
${I_i}$ 、${S_i}$ 、${B_i}$ 和${N_i}$ 分别表示第$i$ 个像素上的测量值、恒星能量强度、杂散光能量强度和噪声强度。杂散光图像中,星点质心可通过公式(2)获得:
$$ \arg \mathop {\min }\limits_{{\mathbf{\theta }},{\mathbf{v}}} \sum\nolimits_i {{{\left[ {{I_i} - B({x_i},{y_i}|{\mathbf{\theta }}) - S({x_i},{y_i}|{\mathbf{v}})} \right]}^2}} $$ (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})$ 存在有如下关系:$$\begin{split} & B(x,y) = B({x_0},{y_0}) + \left( {h\frac{\partial }{{\partial x}} + k\frac{\partial }{{\partial y}}} \right)B({x_0},{y_0}) \cdots {\text{ + }}\\ & \frac{1}{2}{\left( {h\frac{\partial }{{\partial x}} + k\frac{\partial }{{\partial y}}} \right)^2}B({x_0},{y_0}){\text{ + }}o\left( {{{(h + k)}^2}} \right) \end{split}$$ (3) 式中:
$h = x - {x_0}$ ;$k = y - {y_0}$ ;$ o\left( {{{(h + k)}^2}} \right) $ 表示无穷小量。在此忽略,公式(3)重写为公式(4):$$ B(x,y) \approx a \cdot {x^2} + b \cdot {y^2} + c \cdot xy + d \cdot x + e \cdot y + f $$ (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),表示在对噪声的最大似然估计意义下对杂散光曲面的最优估计。
$$ \arg \mathop {\min }\limits_{\mathbf{\theta }} {\mathbf{H}}({x_i},{y_i}|{\mathbf{\theta }}) $$ (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)即可得到杂散光的曲面模型。$$ \frac{{\partial {\mathbf{H}}(x,y|{\mathbf{\theta }})}}{{\partial {\mathbf{\theta }}}} = 0 $$ (6) -
文中构建的杂散光干扰下星点检测方法的基本思路是首先利用星点图像的背景区域像素估计出杂散光在像素上的能量强度,然后通过图像减运算得到去杂散光的图像,最后在此图像中提取星点。算法分为两个过程:杂散光背景估计和星点提取过程。计算流程图如图2所示。
-
根据第1节描述的杂散光曲面模型,杂散光在某像素上的能量强度估计过程可描述为:首先确定像素的邻域,然后根据邻域像素估计杂散光曲面描述函数的参数,最后根据曲面参数估计杂散光在该像素上的能量强度。
某像素邻域为以该像素为中心,取一定长度的矩形区域内的像素作为该像素的邻域,表示为:
$$ \varOmega {\text{ = }}\left\{ {({x_i},{y_i})|{x_i} \in [{x_0} - l,{x_0} + l],{y_i} \in [{y_0} - l,{y_0} + l]} \right\} $$ (7) 式中:
$({x_0},{y_0})$ 表示该像素的位置;$2 l + 1$ 为矩形区域的长度;$l$ 为算法参数。考虑到计算的方便性,采用图像滤波的计算方式,滤波矩阵可表示为:
$$ \begin{array}{*{20}{c}} {\left\{ {\begin{array}{*{20}{c}} {{{\left[ {{M_x}} \right]}_{ij}} = i - l - 1} \\ {{{\left[ {{M_y}} \right]}_{ij}} = j - l - 1} \end{array}} \right.}&{\left\{ {\begin{array}{*{20}{c}} {i = 1,2, \cdots ,2l + 1} \\ {j = 1,2, \cdots ,2l + 1} \end{array}} \right.} \end{array} $$ (8) 相对于杂散光来说,恒星入射的能量是干扰因素。由于星点为直径不超过一定值的光斑,且小于杂散光形成的光斑区域。在公式(8)基础上,对滤波矩阵做进一步的改变,增加了掩膜,可表示为:
$$ \begin{array}{*{20}{c}} {\left\{ {\begin{array}{*{20}{c}} {{{\left[ {{M_x}} \right]}_{ij}} = 0} \\ {{{\left[ {{M_y}} \right]}_{ij}} = 0} \end{array}} \right.}&{\left\{ {\begin{array}{*{20}{c}} {i = [l + 1 - r,l + 1 + r]} \\ {j = [l + 1 - r,l + 1 + r]} \end{array}} \right.} \end{array} $$ (9) 式中:掩膜长度
$r$ 为算法参数。图3所示为文中构建的一个滤波矩阵,其中滤波矩阵的长度为9,掩膜长度为3。
构建的这种滤波矩阵的物理意义为假定当前像素属于星点,取星点周边的背景区域像素来估计当前像素的值。这种方式可排除星点能量对估计杂散光能量的干扰。
在估计杂散光曲面过程时,由于噪声的影响,可能会得到存在较大误差的曲面。为了避免这种情况,文中使用最小均方误差来确定最优曲面模型。
对于像素
$({x_0},{y_0})$ ,该方法可得到${(2 l + 1)^2}$ 个曲面,每个曲面在该像素位置上的均方根误差为${E_j}$ ,其中最小均方误差所对应的曲面即为该方法确定的在该像素上的曲面。$$ \mathop {\min }\limits_{j \in [1,{{(2l + 1)}^2}]} ({E_j}) \Rightarrow \hat \theta $$ (10) 式中:
${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)所示。相比之下,这种方法得到的曲面能更准确描述杂散光能量。
图 4 确定最优杂散光曲面。(a)某一个曲面;(b)最优曲面
Figure 4. Determine the optimal surface of stray light. (a) One of the estimated surfaces; (b) The optimal surface
在获得杂散光曲面参数后,通过公式(11)即可获得杂散光在某像素上的能量强度,则:
$$ \hat B({x_i},{y_i}) \approx {\hat a_i} \cdot {x_M}^2 + {\hat b_i} \cdot {y_M}^2 + {\hat c_i} \cdot {x_M}{y_M} \cdots + $$ $${\hat d_i} \cdot {x_M} + {\hat e_i} \cdot {y_M} + {\hat f_i} $$ (11) 式中:
$({x_M},{y_M})$ 表示像素在公式(10)中所确定的曲面中相对位置。由于星点可能位于图像中的任何位置,该方法要估计图像中每一个像素上的杂散光能量强度,这构成了图像中的杂散光背景。
-
在估计出杂散光背景后,通过图像减运算得到去杂散光的图像,然后采用快速高斯拟合法来提取星点质心,计算流程图如图5所示。
快速高斯拟合法是高斯拟合法的一种改进方法,具有高精度和高计算效率的特点,其星点质心计算精度接近高斯拟合法,高于重心法和加权重心法,计算效率接近质心法,可满足星敏感器实时性需求[11]。另外,采用快速高斯拟合法可构建基于信噪比的筛选像素方法,这有利于获得更多更准确的星点像素。
星点提取过程可分为3个过程,详述如下:
(1) 在去杂散光图像中,利用阈值法检测星点,得到星点的感兴趣区域。
去杂散光的图像可用公式(12)表示:
$$ {G_i} = {I_i} - {\hat B_i} $$ (12) 计算过程采用的阈值由公式(13)获得:
$$ {T_i} = 3\cdot \sqrt {\frac{{\displaystyle\sum\nolimits_{i \in \Omega } {{{({I_i} - {{\hat B}_i})}^2}} }}{{\displaystyle\sum\nolimits_{i \in \Omega } 1 }}} $$ (13) (2)在感兴趣区域内,首先选择那些高亮度的像素并利用快速高斯拟合法估计星点的描述函数参数,然后根据描述函数参数反演星点能量并计算像素信噪比。
星点描述函数定义如公式(14)所示:
$$ S({x_i},{y_i}|{\mathbf{v}}){\text{ = }}A\exp \left( { - \frac{{{{({x_i} - {x_c})}^2}}}{{2\sigma _x^2}} - \frac{{{{({y_i} - {y_c})}^2}}}{{2\sigma _y^2}}} \right) $$ (14) 式中:
${\mathbf{v}} = (A,{\sigma _x},{\sigma _y},{x_c},{y_c})$ ,$({x_c},{y_c})$ 即为星点质心。根据公式(2)、公式(12)和公式(14),快速高斯拟合法的目标函数可以表示为:
$$ \arg \mathop {\min }\limits_{{\mathbf{\theta }},{\mathbf{v}}} \sum\nolimits_i {{{\left[ {{W_i}\left( {\ln {G_i} - \ln S({x_i},{y_i}|{\mathbf{v}})} \right)} \right]}^2}} $$ (15) 在公式(15)中,对数运算可将高斯拟合法中的非线性优化过程转换成线性优化过程,提高了计算效率。在式中添加权重
$W$ 可使快速高斯拟合法的计算结果接近高斯拟合法的计算结果。在不牺牲计算精度的同时,快速高斯拟合法提高了计算效率。根据参考文献[11]的描述,权重可取${W_i} = {G_i}$ 。对目标函数公式(15)求偏导即可获得星点质心。像素信噪比可以表示为:
$$ SN{R}_{i}=\frac{\widehat{S}({x}_{i},{y}_{i}|v)}{{G}_{i}-\widehat{S}({x}_{i},{y}_{i}|v)} $$ (16) (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)所示。图 6 实验使用的杂散光干扰图像。(a)模拟杂散光图像;(b)月光图像;(c)人工光源图像
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)生成,可表示为:
$$ {S_s}({x_i},{y_i}) = \frac{A}{{\sqrt {2\pi {\sigma ^2}} }}\exp \left( { - \frac{{{{({x_i} - {x_c})}^2}{\text{ + }}{{({y_i} - {y_c})}^2}}}{{2{\sigma ^2}}}} \right) + {N_i} $$ (17) 式中:
$ A $ 表示星点的亮度水平;$\sigma $ 表示星点在星敏感器中的离散程度;$ ({x}_{c},{y}_{c}) $ 表示添加到图像中的星点位置。 -
在此次实验中,星点的亮度水平
$A$ 为0.2~0.5之间的随机值。星点的离散程度$\sigma $ 为0.6~1.1之间的随机值。星点的位置在图像内随机选择。图像中添加独立同分布的高斯噪声,噪声的强度依次从0增大到0.1,并且每次噪声强度变化产生1000幅图像,在此次实验中每幅图像中添加400个星点。实验结果如图7和8所示。检测率定义为提取的真星和所有真星之间的比值,误检率定义为提取到的假星与所有提取的星点之间的比值。随着添加的噪声强度增大,背景法确定的阈值会偏大,这使得一些低信噪比的星点不易被检测到,从而导致星点检测率下降。当噪声强度增大时,星点形状特征信息会变得模糊,特别是低信噪比的星点,这使得基于形态学方法检测到的真实星点少,星点检测率降低。对于位于杂散光背景变化显著区域的像素来说,这些像素的灰度值与周围区域的像素存在一定的差异,随着噪声强度增大,这些差异会变的更加显著。如果多个这样的像素聚集在一起,就形成了一个和星点具有类似特征的“星点”——假星,这是算法提取到假星的原因。由于背景法和基于形态学的方法无法识别这类假星,导致星点误检率随噪声强度增大而上升。文中的方法在一定程度上对这类假星不敏感,误检率较低。
-
在提取能力实验的基础上,取那些提取出来的星点,计算星点的质心误差,如公式(18)所示,并使用星点的平均质心误差分析算法性能。结果如图9所示。
$$ Er = \sqrt {{{({{\hat x}_c} - {x_c})}^2} + {{({{\hat y}_c} - {y_c})}^2}} $$ (18) 在杂散光图像中,影响星点质心精度的因素主要有两种:星点图像中的噪声和杂散光。噪声强度增大会引起星点质心误差增大[11]。在此次实验中,随着添加的噪声强度增大,同一算法得到的质心误差曲线会随着增大。由于杂散光能量的存在,公式(2)所描述的星点质心估计过程不再是无偏估计,质心估计结果中会存在一定的估计偏差。基于形态学的方法对杂散光能量的估计较弱,相应的估计结果中存在较大的估计偏差。背景法相当于将杂散光能量建模成一个等值的平面,而文中是将杂散光能量建模成不等值的曲面。相比之下,文中对杂散光能量估计的更准确些,质心结果中存在较小的估计偏差,星点质心误差会更小些。
-
该实验利用月光、人工光源等构建了杂散光实验环境,星敏感器在此环境中采集图像,统计和对比各个算法的提取结果。
在月光和人工光源干扰下得到的星点图像如图5(b)和(c)所示,取杂散光实验环境中获得的10幅图像。在这10幅图像中,各个算法提取的星点统计结果如表1所示。
表 1 杂散光干扰下星点提取结果
Table 1. Star spots extraction results under stray light interference
Evalucation
indicatorMoonlight Artificial light OSFE NSTS TTE OSFE NSTS TTE Detection 203 147 162 173 105 117 False-detetion 7 19 18 4 9 9 在月光干扰实验中,各个算法得到的星点检测个数和误检个数相近,主要原因是月光形成的杂散光背景比较平缓,星点的信噪比较高,这有利于星点的检测。而在仿真图像中,添加了高强度的噪声,这增加了星点检测的难度。
相比于月光干扰,人工光源干扰下得到的杂光背景要更复杂,这使得提取结果变差。从检测到的星点数目和误检数目上看,在人工干扰实验中,文中方法的提取结果要优于背景法和基于形态学方法提取到的结果。
Star spot extraction based on optimal background estimation for star sensor anti stray light
-
摘要: 研究了杂散光干扰下星敏感器星点提取问题。在杂散光干扰下,星点提取的效果会变差,这会使解算的姿态不准确和不可靠。为了抵抗杂散光干扰,提出了一种基于最优背景估计的抗杂散光星点提取方法。首先,分析了杂散光在星敏感器中的成像特点,构建了杂散光的具有闭合解的曲面模型;然后,设计了基于最优曲面的杂散光背景估计方法和星点提取方法;最后,采用仿真的杂散光图像和星敏感器采集的杂散光图像验证了算法的性能。实验结果表明:该方法得到的星点检测率、误检率和质心定位精度均优于现有基于背景和形态学的方法,该方法对杂散光有较好的抑制能力。Abstract: In this article, the problem of extracting star spots for star sensor in stray light is considered. Under this setting, the effect of extracting star spots becomes worse, which affects the accuracy and reliability of the estimated attitude. An extraction method based on optimal surface fitting is proposed for this task. First, the imaging characteristics of stray light in star sensor are analyzed, and the surface model with a closed-form solution is built. Then, the method of estimating optimal stray light background and extracting star spots is proposed. The performance of the proposed method is verified by the simulated star images and star sensor images. Experimental results show that the detection rate, false detection rate and centroid accuracy obtained by the proposed method are better than that obtained by the methods based on threshold and morphological, which shows that the proposed method can resist the disturbance of stray light.
-
Key words:
- star sensor /
- anti stray light /
- optimal background estimation /
- surface fitting /
- star spot extraction
-
表 1 杂散光干扰下星点提取结果
Table 1. Star spots extraction results under stray light interference
Evalucation
indicatorMoonlight Artificial light OSFE NSTS TTE OSFE NSTS TTE Detection 203 147 162 173 105 117 False-detetion 7 19 18 4 9 9 -
[1] Liebe C C. Accuracy performance of star trackers - a tutorial [J]. IEEE Transactions on Aerospace Electronic Systems, 2002, 38(2): 587-599. doi: 10.1109/TAES.2002.1008988 [2] Guo Y, Yang H, Fu W, et al. Temperature control of star sensor baffle using 3 D printing and PCM thermal energy storage technology [J]. International Journal of Heat and Mass Transfer, 2021, 165: 120644. doi: 10.1016/j.ijheatmasstransfer.2020.120644 [3] Yang L I, Liao Z, Shengbo M U, et al. Stray light suppressing technique and simulation for star sensor [J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(12): 2620-2624. (in Chinese) [4] Xu Z, Liu D, Yan C, et al. Stray light elimination method based on recursion multi-scale gray-scale morphology for wide-field surveillance [J]. IEEE Access, 2021, 9(1): 16928-16936. [5] Yu L, Mao X, Jin H, et al. Study on image process method of star tracker for stray lights resistance filtering based on background [J]. Aerospace Shanghai, 2016, 33(4): 26-31. (in Chinese) [6] Yu L, Mao X, Zhou Q, et al. Image processing of star sensor based on maximum background estimation [J]. Laser & Infrared, 2017, 47(7): 889-895. (in Chinese) [7] He Y, Wang H, Feng L, et al. A novel method of eliminating stray light interference for star sensor [J]. IEEE Sensors Journal, 2020, 20(15): 8586-8596. doi: 10.1109/JSEN.2020.2984001 [8] Jiang J, Liu L, Zhang G. Robust and accurate star segmentation algorithm based on morphology [J]. Optical Engineering, 2016, 55(6): 063101. [9] Anderson D R. Taylor’s Formula and Integral Inequalities for Conformable Fractional Derivatives[M]. Berlin: Springer International Publishing, 2016: 25-43. [10] Björck A. Linear Least Squares Problems [M]. Berlin: Springer International Publishing, 2015: 211-430. [11] Wan X, Wang G, Wei X, et al. Star centroiding based on fast Gaussian fitting for star sensors [J]. Sensors, 2018, 18(9): 2836.