先验信息辅助的运动模糊星图星点提取方法

何贻洋，王宏力，冯 磊，由四海，刘 珂，高迎东

(火箭军工程大学 导弹工程学院，西安 710025)

星敏感器是目前公认的精度最高的姿态敏感器之一[1].在航空、航天领域，无论是近地卫星或者深空探测，还是弹道导弹或者小卫星上，星敏感器都得到了广泛的应用.星敏感器通常捷联安装于上述的载体上，若载体受到外界环境的冲击或者进行大姿态角机动时，星敏感器拍摄的星图会出现运动模糊现象，导致星图信噪比降低，可能严重影响星点提取成功率和精度，进而降低星图识别成功率[2-3].

目前对运动模糊星图处理的主流思路是先复原星图再进行后续的星点提取、识别等过程.根据上述思路，运动模糊星图处理的核心在于星图复原，故而有不少学者开展了星图复原的研究.目前星图复原方法主要是非盲图像复原方法，非盲图像复原方法应用最为广泛的是逆滤波方法[4-5]、维纳滤波方法[6-7]和Lucy-Richardson(L-R)方法[8-9].逆滤波方法是应用最早的图像复原方法，但该方法缺点很明显，即对图像噪声非常敏感；维纳滤波方法可以实现最小均方复原效果，但该方法需已知图像噪声特性的峰值信噪比较高时才能获得较好的复原效果；L-R方法是一种统计学方法，如果图像的先验概率模型和条件概率模型与实际不符时，图像复原效果会受到极大影响.另外，非盲图像复原方法受图像先验信息的准确性影响很大，故而对图像的模糊核函数估计也是研究的重点，Radon变换方法[10-11]、基于星像运动轨迹的估计方法[12-14]是当前模糊核函数估计的常用方法.传统的Radon变换方法的估计精度受星图噪声影响很大，Jiang等[15]提出了改进的Radon变换方法，该方法基于模糊星图噪声特性并且结合了Z函数和双线性插值方法，取得了较好的估计效果；吴小娟等[12]基于绕星敏感器光轴或者非光轴转动时的星像运动轨迹模型进而得到相应条件下星图模糊核函数的解析表达式；廖育富等[14]根据星像运动轨迹模型提出了一种操作简单利于工程化实现的模糊核函数估计方法，但估计精度不高且容易受到噪声影响.上述都是基于星图自身信息估计模糊核函数，但星图本身的先验信息有限，所以近几年有学者提出基于惯性导航系统(inertial navigation system，INS)辅助对模糊核函数进行估计[16-18]，取得了较好的估计效果.

综上所述，对于运动模糊星图的研究聚焦于模糊核函数的估计和复原方法，期望获得质量较好的复原星图以保证星点提取精度和星图识别成功率，但星图复原过程会占据星敏感器一定的计算资源，影响星敏感器的数据更新率.因此，本文提出了一种基于INS和星敏感器的先验信息直接提取运动模糊星图上星点的方法.该方法首先利用INS和星敏感器的姿态先验信息确定局部区域生长的“种子”点，从而得到运动模糊星图上的星像区域；然后基于先验信息建立了星点提取误差修正模型，补偿利用质心法直接提取运动模糊星图星点的坐标误差.该方法可以跳过星图复原操作直接对运动模糊星图进行星点提取，这对于保证高动态环境下星敏感器的数据更新率具有重要意义.

1 方法分析

高动态条件下星敏感器成像受到影响，星图会出现模糊现象，图1是模拟的一幅运动模糊星图，模拟星图的灰阶为8位，灰度范围为0～255，添加的高斯白噪声均值为0，均方差为6(灰度值)，载体的角运动速度设置为ωx=-2°/s，ωy=-3°/s，ωz=2°/s.分别利用L-R方法[8-9]和维纳滤波方法[6-7]进行复原处理后，得到局部放大的复原星图以及相应的灰度分布如图2所示.

从上述星图复原的结果来看，利用维纳滤波方法和L-R方法复原星图后，并不能将星点恢复到静态条件下的高斯弥散斑，但星点“拖尾”轨迹相对缩短，利用质心法对复原星图上星点提取的误差小于直接对原始星图上星点提取的误差.如果能够估计出直接利用质心法提取模糊星图上星点的误差，然后对星点坐标进行补偿修正，这就能跳过星图复原这个过程，从而大大节省了时间，保证星敏感器数据更新率，同时仍能获得较好的星点提取精度.基于上述思想，本文提出了基于INS辅助的直接对运动模糊星图星点提取的方法.

图2 不同方法的星图复原结果局部放大和灰度分布图

星敏感器常与INS组合使用，虽然INS存在姿态输出误差随时间累积的缺点，但其短时间的姿态输出精度很高，文献[16-18]中将INS的姿态输出引入到模糊核函数估计中取得了较好的效果，本文进一步利用INS提供的姿态先验信息，提出一种直接对运动模糊星图进行星点提取的方法，具体思路如下：

1)基于INS的短时姿态输出和星敏感器在上一帧星图的姿态先验信息，预测出恒星在当前帧星图的粗略星点坐标，将其作为种子点进行局部区域生长，从而获得星像区域；

2)利用质心法[3]直接在星像区域计算星点质心坐标；

3)根据INS得到在曝光时间内星敏感器的转动角速度，计算星图上星像的运动轨迹长度，根据星像运动轨迹与直接质心法的星点提取误差之间的关系，修正星点坐标. 本文方法的实施步骤如图3所示.

1.1 基于先验信息的局部区域生长

局部区域生长法[19]的三要素是种子点、生长准则和终止条件，而种子点是其中最为关键的一环.利用局部区域生长可以快速而准确获取运动模糊星图上各个星像的范围，但必须首先解决如何获取满足条件的种子点的问题.单帧星图的先验信息有限，本文基于INS从t0→t0+Δt时刻的姿态输出并结合t0时刻星敏感器输出的姿态信息，得到t0+Δt时刻的粗略姿态，解得t0+Δt时刻星敏感器光轴指向，进而预测出t0+Δt时刻恒星在成像平面上的粗略坐标，将其作为种子点，利用局部区域生长法获得运动模糊星图上星像范围.为质心法提取星点坐标奠定基础，具体思路如下：

假设星敏感器根据当前这一帧星图获得t0时刻的姿态四元数为q，利用INS的姿态输出得到从t0→t0+Δt的姿态机动四元数为q′，那么根据四元数的乘法顺序和规则，t0+Δt时刻星敏感器的姿态四元数可以由下式得到:

q″=q′q.

(1)

文献[20]中关于姿态四元数与星敏感器光轴指向(A,D)和星敏感器的滚动角θ之间的关系式：

(2)

(3)

式中:A∈[0,2π)，D∈[-π/2,π/2]，θ∈[0,2π).

得到星敏感器光轴指向(A,D)之后，预测星敏感器在t0+Δt时刻星图上星像坐标就是建立恒星在星图上质心坐标与在天球坐标系中的坐标(α,δ)之间的数学模型，利用t0时刻的星敏感器的先验信息可以得到(α,δ).根据星敏感器成像原理以及星敏感器坐标系与天球坐标系之间的坐标变换关系，在星敏感器光轴指向(A,D)、恒星在天球坐标系中的坐标(α,δ)和星敏感器视场确定的情况下，理论上出现在星图上的恒星位置也随之确定.得到此时恒星的星像坐标为

(4)

如果考虑到在实际中星敏感器的滚动角θ一般不为0，那么则有

(5)

根据式(1)～(5)，可以预测得到t0+Δt时刻星图上星像坐标，将这些星像坐标记为{Si}，其中Si=(xi,yi)，i=1,…,N，N为星图上星像个数. 以{Si}作为局部区域生长的种子点，从而获得模糊星图上各个星像的实际范围关于区域生长法的原理可具体详见文献[19]，在此不再赘述.

1.2 星点提取误差修正模型

质心法具有操作简单、鲁棒性好的优点，是目前常用的星点提取方法.星像的能量分布对质心法的精度有很大影响，根据文献[3]，星像能量分布服从二维高斯分布时，质心法精度较高.在高动态环境下，星像不再是静态时的高斯弥散斑，而是形成了一条“拖尾”.在曝光时间内Te内，星像区域的灰度也不再服从二维高斯分布，因此直接用质心法提取星点的误差会大大增加.为了对该误差进行修正，根据文献[21]的研究，首先应对高动态条件下的运动模糊星图的灰度分布进行分析建模，在此基础上分析推导直接用质心法提取星点的误差.

将星敏感器的曝光时间Te等分为N个微小的区间[ti,ti+1]，i=0,1,…,N-1，Δt=ti+1-ti=Te/N，t0=0，tN=Te，由于[ti,ti+1]非常小，因而可近似认为在[ti,ti+1]内，恒星在成像平面上所形成星像的灰度仍服从二维高斯分布，设在[ti,ti+1]的二维高斯分布的中心点坐标为(xi,yi)，在[ti,ti+1]内，恒星在成像平面上所形成星像的灰度分布如下：

式中:Gi=G0/N，G0是与曝光时间Te内星像总能量相关的常值；Gi(x,y)为[ti,ti+1]内成像平面上的坐标(x,y)处的灰度值;σ为高斯弥散斑半径，代表了弥散程度，根据二维高斯分布的特点以及3σ准则，星像99%以上的能量集中于距离星点的3σ范围以内.

假设各个高斯分布的弥散斑半径相等，将各个微小曝光时间内形成的星像灰度分布函数进行叠加，得到在曝光时间Te内恒星在成像平面上星像的灰度分布函数如下：

利用质心法可以得到[ti,ti+1]内对应的中心点(xi,yi)的计算公式如下：

将式(6)、(7)分别代入式(8)、(9)，则有

根据文献[14,21-22]研究可知绕光轴的角速度对星图的影响很小，因而星图上星像运动轨迹近似为匀速直线运动上的一段，假设星像在成像平面的OX及OY轴上的运动速度分别为vx≈-fωy/dCCD，vy≈fωx/dCCD[12]，ωx、ωy为星敏感器转动角速度大小，逆时针转动为正.则在曝光时间Te内，星像运动轨迹在在成像平面的OX及OY轴上的分量为Lx=vxTe，Ly=vyTe，且有

式中i=0,1,…,N-1，Δt=ti+1-ti.

同理，可得：

(10)

可见，星图运动模糊引起的星点坐标误差等于星像运动轨迹长度的一半.利用INS可得到星敏感器转动的角速度，进而估算出星像运动轨迹长度，根据式(10)可对星点坐标进行误差补偿，从而获得较为准确的星点坐标.

2 方法评价

对于运动模糊星图处理问题，本文提出一种基于INS辅助直接对运动模糊星图进行星点提取的方法.INS的辅助作用主要体现在以下两个方面：1)根据INS的姿态输出可以得到相邻两帧星图之间星敏感器的姿态角变化量，结合星敏感器在上一帧星图的姿态输出，可以预测出当前这帧星图上星点的粗略坐标；2)基于INS的角速度输出可估算曝光时间内的星像运动轨迹长度.

本文方法需要利用恒星在上一帧星图的质心坐标和INS系统的角速度输出等先验信息来预测恒星在当前帧星图中的坐标，作为局部区域生长所需的“种子”点.在实际环境中，先验信息也不可避免的存在误差，那么预测的星点坐标也会产生误差，但若“种子”点仍在实际星像区域范围内，那么局部区域生长就能准确的获得实际星像范围.由此可见，局部区域生长法对“种子”点的精度要求不高，这在一定程度上增强了本文方法对先验信息误差的容错能力；局部区域生长只需关心存在星像的局部星图，星图其他区域不需做进一步处理，这相当于对星图做了一次去噪处理；本文方法避免了扫描整幅星图，而且不需对运动模糊星图进行复原操作，相比于先复原星图再提取星点的方法，本文方法能够大大节省星点提取时间，故本文方法还具有快速性好、抗噪能力强的优点.

值得注意的是本文方法要求已知恒星的先验信息，也就是需要已知恒星在上一帧星图的星点坐标，对应的星敏感器工作模式是星跟踪模式.如果星敏感器处于星捕获模式下，也就是有新的恒星进入视场，而新进入视场内的恒星难以获得其在上一帧星图中的先验信息，那么这种情况下本文方法可能将不再适用.

3 仿真实验及分析

3.1 仿真条件设置

采用变形星敏感器视场设置为9°×9°，焦距为65.76 mm，像元大小为20 um，面阵大小为512×512，曝光时间为100 ms；从SKY2000主星表筛选出星等于小于6的恒星，并剔除其中的双星与变星，剩余4 908颗组成本文的星表；星敏感器角速度设置为(-2°/s -3°/s 2°/s)，角速度方向记顺时针转动为正；模拟星图的灰阶为8位，灰度范围为0～255，添加的高斯白噪声的均值为0，均方差为6(灰度值)；仿真实验工具是MATLAB(R2015a版)，仿真实验的电脑配置为2.5 GHz Core(TM)i7-6500U，8.0 GB RAM.

利用Monte-Carlo方法随机产生100组星敏感器在t0时刻的姿态四元数{qi}和100组t1时刻的光轴指向{(αi,δi)}i=1,…,100，从t0→t1时刻的星敏感器姿态变化四元数{q′i}可利用INS系统的输出而得到，根据式(2)～(5)可得出t1时刻星敏感器成像平面上星像的粗略坐标，将此作为本文方法中局部区域生长法的种子点.

3.2 仿真实验分析

3.2.1 仿真实验1

为了考察“种子”点的误差对于本文方法结果的影响大小，也即考察本文方法对先验信息误差的容错能力，设计了本仿真实验.

为了便于对一幅星图上的星点提取精度进行描述，以星图上提取到的星点坐标与相应的理想星点坐标距离的平均值作为衡量星点提取精度高低的指标，定义如下：

可以假设“种子”点的误差是服从均值为0的高斯白噪声，共设置10组不同大小的“种子”点的误差，分别为0.5,1.0,…,5.0像素.根据t1时刻的光轴指向{(αi,δi)}得到100幅模拟的运动模糊星图，其他仿真条件见仿真条件设置.得到在不同的“种子”点误差条件下，利用本文方法的星点提取误差仿真结果如图4、5所示.图4中的星点提取成功率定义为ratio=Nsuc/Ntotal，其中Ntotal为用于仿真实验的星图总数，Nsuc为星图上的星点被全部成功提取的星图数量.

图4 不同“种子”点误差条件下的星点提取成功率

图5 不同“种子”点误差条件下err的均值

从本仿真实验条件下，“种子”点误差不超过3.5像素时，本文方法的星点提取成功率为100%，“种子”点误差大于3.5像素时，星点提取成功率略有下降；仿真结果如图5所示，星点提取误差err的均值随着“种子”点误差的增大而略有增加，但增幅较小，上述仿真实验结果说明本文方法对“种子”点误差(即先验信息误差带来的影响)具有较好的容错能力.

3.2.2 仿真实验2

根据t1时刻的光轴指向{(αi,δi)}得到100幅模拟的运动模糊星图，其他仿真条件见仿真条件设置.利用本文方法直接对模糊星图上星像质心进行提取和补偿；作为比较，利用L-R方法和维纳滤波方法对模糊星图进行复原，然后对复原星图利用块扫描法[23]+质心法进行星点提取，分别简记为Method-1和Method-2，L-R方法的迭代次数设置为10.得到仿真结果如图6、7所示.

图6 3种方法err值的对比

利用本文方法对100幅运动模糊星图进行星点提取以及利用L-R方法以及维纳滤波方法先对运动模糊星图复原再进行星点提取，100幅星图的星点提取精度如图6所示，本文方法、Method-1以及Method-2这3种方法err的平均值分别为0.658、1.063、1.205像素，可以看出本文方法的星点提取误差小于Method-1和Method-2；3种方法的耗时情况如图7的所示，本文方法、Method-1和Method-2的平均耗时分别为0.063 2、0.555 3、0.112 0 s，可知本文方法的耗时小于Method-1和Method-2.从上述仿真实验结果可知，相比于L-R方法和维纳滤波方法，本文方法在快速性方面具有一定的优势，星点提取误差虽然比其他两种方法有所降低，但误差仍然较大.

3.2.3 仿真实验3

为了更全面的分析星图噪声大小对星点提取的影响，研究本文方法的抗噪性能，设计如下仿真实验.在模拟星图上加入不同强度的高斯白噪声，共设置10组不同强度的高斯白噪声，高斯白噪声均值为0，均方差分别为1,2,…,10(灰度值，星图灰度值范围为0～255)，其他仿真条件与仿真实验2相同.得到仿真结果如图8～10所示.

图8 不同噪声条件下3种方法的星点提取成功率对比

图9 不同噪声条件下3种方法耗时的对比

根据仿真结果图8可知，随着星图噪声增大，3种方法都受到了影响，Method-1的星点提取成功率下降幅度最大，Method-2次之，本文方法的星点提取成功率下降幅度最小，在仿真实验中，星点提取成功率最低为0.8，可见本文方法具有较好的抗噪性能.从仿真结果图9可知，3种方法的平均耗时随着噪声的变化，波动较小，相比于Method-1和Method-2，本文方法的速度更快.为了考察在不同噪声条件下，当星图上的星点被成功提取时3种方法的星点提取精度.在上述仿真实验条件下，得到在不同噪声条件下3种方法的星点提取误差的仿真结果如图10所示，图10中的星点提取平均误差为每组噪声条件下err的均值.

图10 不同噪声条件下3种方法的星点提取误差对比

从图10的仿真结果可知，随着噪声的增大，3种方法的星点提取误差均有逐渐增大的趋势，但本文方法的星点误差远比其他两种方法要小，以噪声均方差等于10为例，本文方法的星点提取平均误差为0.957像素，Method-1的星点提取平均误差为1.784像素，Method-2的星点提取平均误差为2.470像素.

根据仿真实验2的结果，无论从星点提取成功率还是从星点提取平均误差来讨论，本文方法在不同噪声水平下的表现都优于Method-1和Method-2，这说明本文方法具有较好的抗噪性能.

4 结 论

1)利用INS的输出和星敏感器自身已经获得的姿态先验信息预测恒星在当前帧星图上的粗略坐标，将其作为“种子”点，局部区域生长对“种子”点误差具有较好的容错能力，这在一定程度上减轻了先验信息误差对本文方法的影响.

2)建立了利用质心法提取运动模糊星图上星点质心坐标的误差修正模型.本文方法不需要对运动模糊星图进行复原且不用对整幅星图进行扫描，这大大节省了运算时间，有利于提高星敏感器数据更新率.

3)仿真实验表明，本文方法的快速性和抗噪能力均优于常用的L-R方法和维纳滤波方法，星点提取精度也得到了一定的提高.