CN117422742A - 一种基于卡尔曼滤波的星敏感器星点质心提取方法 - Google Patents
一种基于卡尔曼滤波的星敏感器星点质心提取方法 Download PDFInfo
- Publication number
- CN117422742A CN117422742A CN202311603531.5A CN202311603531A CN117422742A CN 117422742 A CN117422742 A CN 117422742A CN 202311603531 A CN202311603531 A CN 202311603531A CN 117422742 A CN117422742 A CN 117422742A
- Authority
- CN
- China
- Prior art keywords
- star
- centroid
- star point
- time
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 49
- 238000000605 extraction Methods 0.000 title claims abstract description 30
- 238000003384 imaging method Methods 0.000 claims abstract description 55
- 238000000034 method Methods 0.000 claims abstract description 40
- 238000005259 measurement Methods 0.000 claims abstract description 30
- 238000012545 processing Methods 0.000 claims abstract description 13
- 239000013598 vector Substances 0.000 claims description 98
- 239000011159 matrix material Substances 0.000 claims description 67
- 230000001133 acceleration Effects 0.000 claims description 11
- 238000004364 calculation method Methods 0.000 claims description 10
- 230000007704 transition Effects 0.000 claims description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical group C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 2
- 238000004088 simulation Methods 0.000 description 11
- 230000008569 process Effects 0.000 description 7
- 238000002474 experimental method Methods 0.000 description 5
- 230000005540 biological transmission Effects 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 230000035945 sensitivity Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000011045 prefiltration Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/277—Analysis of motion involving stochastic approaches, e.g. using Kalman filters
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
- G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Multimedia (AREA)
- Navigation (AREA)
Abstract
本发明公开了一种基于卡尔曼滤波的星敏感器星点质心提取方法，涉及星点质心提取技术领域，包括步骤：首先，基于星点运动学模型生成星点预测窗口，其次，采用经典的质心方法计算预测窗口内星点的含噪声质心，最后，用每个卡尔曼滤波器分别对上述质心进行处理，对图像噪声进行滤波，得到高精度质心，仿真验证了基于卡尔曼滤波的估计模型。本发明利用当前的去噪质心和全运动参数来生成下一时刻的星点预测窗口，在这种循环中可以不断地实现高精度的星点去噪质心提取，提高了多曝光成像星敏感器的测量精度。
Description
技术领域
本发明涉及星点质心提取技术领域，特别涉及一种基于卡尔曼滤波的星敏感器星点质心提取方法。
背景技术
星敏感器目前被认为是最精确的姿态测量仪器，而且它还具有无漂移和自主性的特点。因此，星敏感器在航空航天导航领域得到了广泛的应用。然而，传统星敏感器的姿态更新率较低，通常仅为4-10Hz，无法独立完成导航任务。为了提高姿态更新率，星敏感器和惯性测量单元的组合导航系统已经成为一种经典的解决方案。
星敏感器的姿态数据没有漂移、但更新率低，而惯性测量单元的更新率高、但姿态数据随时间漂移。当上述两个导航装置结合时，可以定期使用星敏感器的姿态数据来校正惯性测量单元的数据漂移。组合导航系统可以输出无漂移的姿态数据，且具有高更新率。此外，当采用三视场星敏感器时，姿态更新率可达到单视场星敏感器的3倍。这是因为三个视场可以用来交替输出姿态数据。然而，上述两种方案也存在一些明显的缺点，如体积大、重量重、功耗高等缺点，不利于姿态测控系统向轻小型化方向的整体发展。
相比之下，对单一星敏感器姿态更新率提升方法的研究具有更广泛的应用价值。有关学者研究了星敏感器的工作特性，并将其工作过程分为三个阶段。星图曝光成像和星图像素数据的传输和处理，以及星跟踪和姿态计算。随后，他们提出了一种并行处理方法，从而可以将上述三个工作阶段流水线化，从而提高星敏感器的姿态更新率。在这种情况下，姿态更新率完全由上述三个工作阶段中最耗时的一个决定。对于传统的星敏感器，由于图像探测器的灵敏度较低，必须进行较长的曝光时间以获得足够的恒星探测灵敏度。因此，曝光时间是传统星敏感器姿态更新率的主要瓶颈。为了减少曝光时间，一些学者将高灵敏度图像探测器引入星敏感器领域，如电子倍增电荷耦合器件和增强型电荷耦合器件等。曝光时间的缩短在一定程度上提高了星敏感器的姿态更新率。此时，星图像素数据的传输和处理时间成为新瓶颈。为了进一步提高，Wang等人提出了一种基于分布式并行超块的星点质心提取方法，Ding等人还提出了一种多通道星点质心提取方法。上述两种方法都利用了并行处理的思想来提取星点质心，从而有效地减少了像素数据的处理时间。
然而，由于图像探测器本身的输出容量限制，像素数据的传输时间仍然不能显著减少。同时，Liebe的研究指出，在不改变其他参数的情况下，随着图像探测器的像素分辨率提高，星敏感器的姿态测量精度得到了相应提高。因此，大阵列图像探测器已逐渐应用于星敏感器领域。在这种情况下，星图像素数据的传输时间消耗更为严重，这极大地限制了姿态更新率的进一步提升。
为了打破上述限制，一些学者提出了一种卷帘快门曝光成像方法。在这种方法中，单帧星图中的每个星点实际上是不同时刻的成像结果。当姿态跟踪和单星点测量相结合时，每个星点都可以获得一组姿态数据，从而提高了姿态更新率。但是，单星质心的随机误差严重影响了姿态测量的精度，从而难以保证星敏感器的姿态测量精度。
为了显著提高姿态更新率，并克服上述方法的缺点，作者曾提出一种多曝光成像方法。根据载体的角速度，可以将N次的短曝光自适应地插入到一个曝光周期中，从而将N个时刻的星点记录在单帧星图中。当分别根据N个时刻进行星跟踪和姿态计算时，可以得到N个时刻对应的姿态信息，相当于将姿态更新率提升N倍，然而，由于多次曝光，在星图中也引入了严重的噪声和虚假星点，严重不利于实际应用。
发明内容
本发明的目的在于针对上述现有技术的不足，提供一种基于卡尔曼滤波的星敏感器星点质心提取方法，以解决现有技术中星点噪声大和虚假星点多等严重影响星敏感器姿态测量精度的问题，从而保证多曝光成像方法在星敏感器中的应用。
本发明具体提供如下技术方案：一种基于卡尔曼滤波的星敏感器星点质心提取方法，包括以下步骤：
结合前一帧星图中恒星星点的质心，以及星敏感器的全运动参数X(t)，生成星点预测窗口，并利用星点预测窗口对当前帧星图进行处理，获得星点图像；
采用质心法计算所述星点图像的质心，获得含噪声质心；
使用卡尔曼滤波器对所述含噪声质心中包含的图像噪声进行滤波，获得去噪质心；
通过所述去噪质心估计当前的姿态矩阵和全运动参数，作为最终输出信息；
使用所述最终输出信息生成下一时刻星点预测窗口，通过所述下一时刻星点预测窗口进行下一时刻星点质心的提取，通过循环往复所述下一时刻星点质心的提取过程，输出星点符合阈值的质心。
优选的，使用所述最终输出信息生成下一时刻星点预测窗口，包括如下步骤：
构建多曝光成像星敏感器的坐标系；
在所述坐标系中对恒星i进行处理，获得其在坐标系中的观测单位矢量；
获取星敏感器在时刻tk的姿态矩阵，并通过姿态矩阵获取所述观测单位矢量与参考单位矢量之间的关系；
通过所述观测单位矢量与参考单位矢量之间的关系，获得星敏感器在不同时刻的旋转矩阵；
在时间间隔小于阈值时，将旋转矩阵近似值、观测单位矢量与时间之间的关系结合，生成下一时刻星点预测窗口。
优选的，在所述坐标系中对恒星i进行处理，获得其在坐标系中的观测单位矢量，包括如下步骤：
获取恒星i的赤经αi和赤纬δi，并获得赤经αi和赤纬δi在坐标系中的参考单位矢量Vi；
将恒星i进行多曝光成像，设置曝光次数为N，获得N个时刻的恒星i的成像星点图像；
通过所述星点图像获取坐标系中恒星i的观测单位矢量；
其中，xi(tk)与yi(tk)分别为星点在时刻tk的坐标，f为星敏感器的焦距。
优选的，所述获取星敏感器在时刻tk的姿态矩阵，并通过姿态矩阵获取所述观测单位矢量与参考单位矢量之间的关系，包括如下步骤：
令A(tk)为星敏感器在tk时刻的姿态矩阵；
通过所述姿态矩阵表示所述观测单位矢量与参考单位矢量之间的关系Wi(tk)，具体表达式为：
Wi(tk)＝A(tk)Vi(k＝1,2,...,N；i-1,2,...,M) (3)
其中，M表示在多曝光成像星敏感器的视场中导航星的总数；
以同样方式获取tk+1时刻的观测单位矢量与参考单位矢量之间的关系Wi(tk+1)；
Wi(tk+1)＝A(tk+1)Vi (4)
将tk时刻的观测单位矢量与参考单位矢量之间的关系表达式代入tk+1时刻的表达式，获得时刻tk、时刻tk+1中观测单位矢量之间的关系表达式为：
Wi(tk+1)＝A(tk+1)Vi＝A(tk+1)A(tk)TWi(tk) (5)
其中，A(tk+1)为星敏感器在tk+1时刻的姿态矩阵。
优选的，通过所述观测单位矢量与参考单位矢量之间的关系，获得星敏感器在不同时刻之间的旋转矩阵，包括如下步骤：
定义R(tk,tk+1)作为多曝光成像星敏感器从时刻tk到tk+1的旋转矩阵，具体表达式为：
R(tk,tk+1)＝A(tk+1)A(tk)T (6)。
优选的，所述在时间间隔小于阈值时，将旋转矩阵近似值、观测单位矢量与时间之间的关系结合，生成下一时刻星点预测窗口，包括如下步骤：
令时刻tk和时刻tk+1之间的时间间隔小于阈值，此时姿态矩阵满足：
R(tk,tk+1)＝I+[[Δθ]]+O(|Δθ|2) (7)
其中，I为3×3阶单位矩阵，Δθ＝[ΔθxΔθyΔθz]T是最小旋转角矢量，[[Δθ]]和O(|Δθ|2)分别为Δθ的交叉乘积矩阵和高阶无穷小矩阵；θx、θy、θz分别为x、y、z轴上的旋转角；
在时间间隔Δt＝tk+1-tk小于阈值时，将多曝光成像星敏感器的运动看作是均匀角加速度运动；
令ω(tk)＝[ωx(tk)ωy(tk)ωz(tk)]T和α(tk)＝[αx(tk)αy(tk)αz(tk)]T分别表示星敏感器在时刻tk时的角速度和角加速度，且角速度和角加速度为全运动参数，则瞬时角速度ω(t)＝[ωx(tk)ωy(tk)ωz(tk)]T表示为：
ω(t)＝ω(tk)+α(tk)(t-tk)(tk≤t≤tk+1) (8)
对公式(8)进行积分，获得最小旋转角矢量；
忽略高阶无穷小矩阵O(|Δθ|2)时，旋转矩阵R(tk,tk+1)近似为：
将公式(10)与观测单位矢量同时代入时刻tk、时刻tk+1中观测单位矢量之间的关系表达式，获得：
通过公式(11)，将观测单位矢量的第三分量之间关系表示为：
并通过公式(12)推导出：
将公式(13)带入单位观测矢量之间关系表达式中，转换为：
由于旋转角度小于阈值，最大计算结果小于0.001时进行忽略，公式(14)近似为：
考虑到焦距f是常数，并且(xi(tk),yi(tk))是成像平面中星点的坐标，公式(15)相应地重写为：
并简化为二维公式：
通过最终的二维公式生成下一时刻星点预测窗口。
优选的，所述采用质心法计算所述星点的质心，获得含噪声质心，包括如下步骤：
获取星点的有效像素数据，使用每个预测窗口内星点的像素数据来计算真实星点的含噪声质心；
其中，计算真实星点的含噪声质心时采用质心方法，具体表达式为：
其中，m、n、f(m,n)为预测窗口中每个有效像素的x坐标、y坐标和灰度值，(xic(tk+1),yic(tk+1))是计算出的含噪声质心。
优选的，使用卡尔曼滤波器对所述含噪声质心中包含的图像噪声进行滤波，获得去噪质心，包括如下步骤：
Xk+1＝Φk+1,kXk+BkUk+Wk (19)
其中，Xk+1＝[xi(tk+1),yi(tk+1)]T和Xk＝[xi(tk),yi(tk)]T分别是时刻tk和时刻tk+1的状态向量，Φk+1,k为状态转移矩阵，其表示为：
其中BkUk是确定性驱动项，其表示为：
其中，f为星敏感器的焦距。
优选的，通过质心方法的计算公式描述了tk+1时刻的量测向量，其中量测方程表达式为：
Zk+1＝Hk+1Xk+1+Vk+1, (22)
其中，Zk+1＝[xic(tk+1),yic(tk+1)]T和Xk+1＝[xi(tk+1),yi(tk+1)]T分别是量测向量和状态向量，Hk+1＝I2×2是量测矩阵；Vk+1表示量测噪声，其协方差矩阵为：
Rk+1＝μI2×2(μ>0) (23)
其中μ与星图的噪声水平有关。
优选的，通过状态方程和量测方程构建星点质心滤波估计模型，具体表达式为：
其中，k和k+1分别表示当前和下一次曝光时刻；和/>是状态向量的估计，它对应于两个曝光时刻k和k+1；/>和P(k+1)k分别为一步预测向量和一步预测误差矩阵；P(k+1)和K(k+1)为估计误差矩阵和滤波器增益矩阵；当给定初始值/>和P(0)时，星点质心矢量/>根据上述公式估计。
与现有技术相比，本发明具有如下显著优点：
本发明通过星点预测窗口对当前时刻的星图进行处理，然后只将预测窗口内的星点像素数据传输到下一步，采用经典的质心法计算上述星点的质心，这些质心不可避免地包含许多图像噪声，使用卡尔曼滤波器分别对含噪声质心进行处理，对图像噪声进行滤波，得到初步的去噪质心，降低星图中的部分噪声和虚假星点，且利用去噪质心估计姿态矩阵和全运动参数作为最终输出信息，生成下一时刻的星点预测窗口，再次进行星图处理，在这种循环中不断地实现更高精度的星点去噪质心提取，提高了多曝光成像星敏感器的测量精度。
附图说明
图1为本发明提供的星点质心提取方法的流程图；
图2为本发明提供的多曝光成像星敏感器的基本工作原理示意图；
图3为本发明提供的±0.4像素噪声下的坐标估计误差图，其中(a)为x坐标结果，(b)为y坐标结果；
图4为本发明提供的±0.4像素噪声下的质心最大绝对距离误差图；
图5为本发明提供的±1.0像素噪声下的坐标估计误差图，其中(a)为x坐标结果，(b)为y坐标结果；
图6为本发明提供的±1.0像素噪声下的质心最大绝对距离误差图；
图7为本发明提供的±2.5像素噪声下的坐标估计误差图，其中(a)为x坐标结果，(b)为y坐标结果；
图8为本发明提供的±2.5像素噪声下的质心最大绝对距离误差图；
图9为本发明提供的夜场实验图；
图10为本发明提供的一帧多曝光成像星图。
具体实施方式
下面结合本发明中的附图，对本发明实施例的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。
为了显著提高姿态更新率，并克服现有技术的缺点，发明人在之前的研究中提出了一种多曝光成像方法。根据载体的角速度，可以将N次的短曝光自适应地插入到一个正常的曝光周期中，从而将N个时刻的星点记录在单帧星图中。当分别根据N个时刻进行星跟踪和姿态计算时，可以得到N个时刻对应的姿态信息，相当于将星敏感器的姿态更新率提升N倍。然而，由于多次曝光，在整个星图中也引入了更严重的成像噪声和更多的假恒星点(虚假目标)，这不利于真实星点的精确质心提取，从而影响了多曝光成像星敏感器的精度，不利于多曝光成像的实际应用。为了保证多曝光成像星敏感器的有效应用，本发明提出了一种基于卡尔曼滤波的星敏感器星点质心提取方法。它对去除假星点目标和提高质心精度具有显著效果，从而有利于多曝光成像星敏感器的实际应用。
如图1所示，所提出的星点质心提取方法可分为三个部分，星点预测窗口，质心计算和卡尔曼滤波。本发明通过结合前一时刻星图中星点的高精度质心，以及全运动参数，生成星点预测窗口，通过星点预测窗口对当前时刻的星图进行处理，然后只将预测窗口内的星点像素数据传输到下一步，采用经典的质心法计算上述星点的质心，这些质心不可避免地包含许多图像噪声，因此只是含噪声质心；最后，用每个卡尔曼滤波器分别对含噪声质心进行处理，对图像噪声进行滤波，得到高精度质心；利用高精度质心，可以估计姿态矩阵和全运动参数作为最终输出信息。同时，可以利用当前的高精度质心和全运动参数来生成下一时刻的星点预测窗口，因此，在这种循环中可以不断地实现高精度的星点质心提取。
下面对本发明的内容进行详细说明，具体包括如下步骤：
步骤S1：结合前一时刻星图中星点的高精度质心，以及全运动参数，生成星点预测窗口，利用星点预测窗口对当前时刻的星图进行处理，获得星点图像。
其中，通过星点预测窗口对当前时刻的星图进行处理。本步骤中，结合前一帧星图中恒星星点的质心，以及全运动参数，生成下一时刻星点预测窗口，包括如下步骤：
步骤S11：构建多曝光成像星敏感器的坐标系。图2显示了多曝光成像星敏感器的基本工作原理。Os-XsYsZs是星敏感器的坐标系。其坐标原点Os位于图像探测器的成像平面的中心位置。Xs-Os-Ys平面与成像平面重合，Xs轴和Ys轴分别平行于成像平面的水平方向和垂直方向。Zs轴指向视轴的方向。Oe-XeYeZe是天球坐标系。它是一个惯性系，它的坐标原点Oe位于地球的中心。这个平面与天球赤道平面重合。Xe轴指向春分点，Ze轴指向天北极，Ye轴由右手坐标系决定。
步骤S12：在坐标系中对恒星i进行处理，获得其在坐标系中的观测单位矢量。在整个天球中，恒星随时间移动得非常缓慢。在修正了岁差和章动后，恒星的位置可以看作是与时间无关的常数。以恒星i为例，获取其赤经(从Xe轴到恒星i在Xe-Oe-Ye平面投影的旋转角度)和赤纬(恒星i在Xe-Oe-Ye平面投影到恒星i的旋转角度)分别为(αi、δi)，并获得赤经αi和赤纬δi天体坐标系Oe-XeYeZe中的参考单位矢量Vi，可表示为：
将恒星i进行多曝光成像，设置曝光次数为N，获得恒星i的N个成像星点的星图。
通过星图获取坐标系中恒星i的观测单位向量；
其中，xi(tk)与yi(tk)分别为星图平面在时刻tk星点的坐标，f为多曝光成像星敏感器的焦距。
步骤S13：获取(多曝光成像，下面均以多曝光成像对星敏感器进行限定)星敏感器在时刻tk的姿态矩阵，并通过姿态矩阵获取所述观测单位矢量与参考单位矢量之间的关系。具体为：
令A(tk)为多曝光成像星敏感器在tk时刻的姿态矩阵。
通过姿态矩阵表示观测单位矢量与参考单位矢量之间的关系Wi(tk)，具体表达式为：
Wi(tk)＝A(tk)Vi(k＝1,2,...,N；i-1,2,...,M)
其中，M表示在多曝光成像星敏感器的视场中导航星的总数。
以同样方式获取时刻tk+1的观测单位矢量与参考单位矢量之间的关系Wi(tk+1)。
Wi(tk+1)＝A(tk+1)Vi
将tk时刻的观测单位矢量与参考矢量之间的关系表达式代入tk+1时刻的表达式，获得时刻tk到tk+1中观测单位矢量与参考单位矢量之间的关系的表达式为：
Wi(tk+1)＝A(tk+1)Vi＝A(tk+1)A(tk)TWi(tk)
其中，A(tk+1)为星敏感器在tk+1时刻星敏感器的姿态矩阵。
步骤S14：通过观测单位矢量与参考单位矢量之间的关系，获得多曝光成像星敏感器在不同时刻的旋转矩阵。具体为：
定义R(tk,tk+1)作为多曝光成像星敏感器从时刻tk到tk+1的旋转矩阵，具体表达式为：
R(tk,tk+1)＝A(tk+1)A(tk)T。
步骤S15：令时刻tk和时刻tk+1之间的时间间隔小于阈值，将旋转矩阵近似值、观测单位矢量与时间之间的关系结合，生成下一时刻星点预测窗口，此时旋转矩阵满足：
R(tk,tk+1)＝I+[[Δθ]]+O(|Δθ|2)
其中，I为3×3单位矩阵，Δθ＝[ΔθxΔθyΔθz]T是最小的旋转角度向量，[[Δθ]]为Δθ和O(|Δθ|2)的高阶无穷小矩阵。θx、θy、θz分别为x、y、z轴上的旋转角。
在时间间隔小于阈值时，Δt＝tk+1-tk，将多曝光成像星敏感器的运动看作是均匀角加速度运动。
令ω(tk)＝[ωx(tk)ωy(tk)ωz(tk)]T和α(tk)＝[αx(tk)αy(tk)αz(tk)]T分别表示星敏感器在时刻tk时的角度速度和角度加速度，且角速度和角加速度为全运动参数，则瞬时角速度ω(t)＝[ωx(tk)ωy(tk)ωz(tk)]T表示为：
ω(t)＝ω(tk)+α(tk)(t-tk)(tk≤t≤tk+1)
对上式进行积分，获得最小旋转角度矢量。
忽略高阶无穷小矩阵O(|Δθ|2)时，旋转矩阵R(tk,tk+1)近似为：
将上式与观测单位矢量同时代入时刻tk到tk+1中观测单位矢量与参考矢量之间的关系的表达式，获得：
通过上式将观测单位矢量的第三分量的关系表示为：
并通过上式推导出：
对基于上式带入单位矢量与时间之间关系的表达式中，转换为：
由于旋转角度小于阈值，最大计算结果小于0.001时进行忽略，上式近似为：
考虑到焦距f是常数，并且(xi(tk),yi(tk))是图像平面中星点的坐标，上式相应地重写为：
并简化为二维公式：
通过最终的二维公式生成下一时刻星点预测窗口。
步骤S2：采用质心法计算星点图像的质心，获得含噪声质心。
获取星点图像的有效像素数据，使用每个预测窗口内星点的像素数据来计算真实的星点的含噪声质心。
其中，计算真实星点的含噪声质心时采用质心方法，具体表达式为：
其中，m、n、f(m,n)为预测窗口中每个有效像素的x坐标、y坐标和灰度值，(xic(tk+1),yic(tk+1))是计算出的含噪声质心。
步骤S3：使用卡尔曼滤波器对含噪声质心中包含的图像噪声进行滤波，获得去噪(高精度)质心。
获取卡尔曼滤波的状态方程；
Xk+1＝Φk+1,kXk+BkUk+Wk
其中，Xk+1＝[xi(tk+1),yi(tk+1)]T和Xk＝[xi(tk),yi(tk)]T分别是时刻tk和时刻tk+1的状态向量，Φk+1,k为状态转移矩阵，其表示为：
其中BkUk是确定性驱动项，其表示为：
其中，f为星敏感器的焦距。
通过质心方法的计算公式描述了tk+1时刻的量测向量，其中量测方程表达式为：
Zk+1＝Hk+1Xk+1+Vk+1,
其中，Zk+1＝[xic(tk+1),yic(tk+1)]T和Xk+1＝[xi(tk+1),yi(tk+1)]T分别是量测向量和状态向量，Hk+1＝I2×2是量测矩阵；Vk+1表示量测噪声，其协方差矩阵为：
Rk+1＝μI2×2(μ>0)
其中μ与星图的噪声水平有关。
步骤S4：通过去噪质心估计当前的姿态矩阵和全运动参数，作为最终输出信息。
通过状态方程和量测方程构建基于卡尔曼滤波的星点质心滤波估计模型，具体表达式为：
其中，k和k+1分别表示当前和下一次曝光时刻；和/>是状态向量的估计，它对应于两个曝光时刻k和k+1；/>和P(k+1)k分别为一步预测向量和一步预测误差矩阵；P(k+1)和K(k+1)为估计误差矩阵和滤波器增益矩阵；当给定初始值/>和P(0)时，星点质心矢量/>根据上述公式估计。
步骤S5：使用最终输出信息生成下一时刻星点预测窗口，通过下一时刻星点预测窗口进行下一时刻星点质心的提取，通过循环往复下一时刻星点质心的提取过程，输出星点符合阈值的质心。
仿真和实验：
通过仿真验证了上述基于卡尔曼滤波的估计模型的可行性和有效性。多曝光成像星敏感器的仿真参数如表1所示。多曝光成像星敏感器的总周期T和多曝光次数N分别设置为100ms和10，因此每次曝光之间的采样间隔为10ms。
表1多曝光成像星敏感器的仿真参数
| 参数 | 值 |
| 周期T(ms) | 100 |
| 多曝光次数N | 10 |
| 探测器分辨率S×S(像素) | 2,048×2,048 |
| 像素大小：a×a(μm) | 5.5×5.5 |
| 焦距(mm) | 31.94 |
| 视场(°) | 20 |
整个仿真都运行在仿真计算机上，其具有Windows 10操作系统和主频率为2.60GHz的CPU。仿真计算软件选择MATLAB(版本号：R2014a)。仿真过程可分为几个步骤。首先，根据史密森天体物理天文台(SAO)的恒星星表，只选择星等小于等于6.0的恒星来建立多曝光成像星敏感器的导航恒星星表。然后，随机生成一组姿态数据作为多曝光成像星敏感器的初始姿态数据，在该姿态下，共获得视场内的M＝30颗导航星。第三，角速度的初始值ω(t0)和角加速度初始值α(t0)，假定为[10°/s 00]T和[1.0°/s200]T。随后，在这30颗导航星中，随机选择第10颗导航星作为测试对象。当多曝光成像星敏感器在上述假定角速度下移动时，可以得到第10颗导航星在连续采样时刻St下的质心数据，其总采样次数设定为300，并作为真实值。最后，在真实质心数据中加入三种不同大小的高斯白噪声，并采用基于卡尔曼滤波的估计模型对质心数据进行滤波。首先将噪声大小设置为零均值和±0.4像素的最大偏差，质心坐标的估计误差如图3和图4所示。在图3中，x坐标的估计误差的标准差分别为滤波前0.24像素和滤波后0.09像素；y坐标的估计误差的标准差分别为滤波前0.23像素和滤波后0.08像素。x坐标和y坐标的估计误差均减小，减小量可以分别计算为：
为了进行合理的评价，我们选择了较大的37.5％的值作为评价结果。也就是说，经过卡尔曼滤波估计后，坐标误差减少到原始坐标误差的37.5％。在图4中，计算了质心的真值和估计值之间的绝对距离误差，可以看出，经过滤波后的绝对距离误差显著降低。随后，将噪声大小放大为零均值和±1.0像素的最大偏差，质心坐标的估计误差如图5和图6所示。在图5中，x坐标估计误差的标准差分别为滤波前0.57像素和滤波后0.14像素；y坐标估计误差的标准差分别为0.57像素和滤波后0.15像素。x坐标和y坐标的估计误差均减小，减小量可以分别计算为：
同样，选择较大的值26.3％作为评价结果。因此，经过卡尔曼滤波估计后，坐标误差降低到原始坐标误差的26.3％。在图6中，滤波后的绝对距离误差也显著降低。最后，为了模拟基于卡尔曼滤波的估计模型在很大误差下的性能，将噪声大小进一步扩大为零均值和±2.5像素偏差，估计误差如图7和图8所示。在图7中，x坐标估计误差的标准差分别为滤波前1.46像素和滤波后0.29像素；y坐标估计误差的标准差分别为滤波前1.45像素和滤波后0.30像素。x
坐标和y坐标的估计误差的标准差均减小，减小量可以分别计算为：
同样，经过卡尔曼滤波估计后，坐标误差降低到原始坐标误差的20.7％。在图8中，滤波后的绝对距离误差也显著降低。根据图4、图6、图8的结果，得到了不同大小噪声作用下的绝对距离误差的均值和标准差见表2。在表2中，卡尔曼滤波估计的绝对距离误差随着噪声的增大而增加，但增加幅度小于噪声的增加幅度。当噪声大小增加为零均值和±2.5像素最大偏差时，滤波前均值为1.92像素，滤波后为0.40像素，滤波前标准差为0.76像素，滤波后为0.20像素。可以看出，经过卡尔曼滤波估计后，绝对距离误差的均值和标准差显著降低。
表2绝对距离误差的均值和标准差表
实验：
通过夜场观星实验进一步验证了星点预测窗口的可行性和有效性。图9为夜场观星实验设备。主要设备包括多曝光成像星敏感器、便携式转台、三脚架、计算机和电源。图9中的多曝光成像星敏感器所采用的主要设计参数，如焦距、视场等，与表1中所述完全相同，一个便携式转台被牢固地安装在一个三脚架上。多曝光成像星敏感器已安装在转台上，并随之旋转。所用计算机与在仿真中使用的计算机相同。电源用于提供220V交流和5V直流电压源。
转台的角速度和角加速度设为5°/s和0.5°/s2，多曝光次数N设定为5。在这些条件下，多曝光成像星敏感器开始工作，多曝光成像星敏感器获得的其中一帧星图如图10所示。
在图10中，星图中总共有8颗导航星，编号为1-8，其余为带有噪声的星空背景。由于多曝光次数N设为5，每颗导航星有5个星点。此外，还可以看出，星图的噪声水平相对较高。当使用全局阈值分割(全局阈值设定为128)和经典的质心方法时，共有4985个星点对象。其中，单一像素对象3664个，占总数的73.5％，双像素对象749个，占总数的15.0％，三像素及以上对象572个，占总数的11.5％。虽然大多数对象都是假目标(1像素和2像素对象)，并且可以消除，但候选星点对象的数量仍然很大(572个超过3像素的对象)，这严重影响了真实星点对象的精确提取。相比之下，当使用星点预测窗口进行质心提取时，只提取了8×5＝40个星点对象，仅占总数的0.8％。同理，在相同的实验条件下，又获得了99张星图，以进一步验证星点预测窗口。在100张星图中，有效星点对象数在每张星图总目标数中的平均占比仅为0.95％，这极大地证明了该方法的有效性。
总之，在仿真过程中，在真实质心数据中加入了不同大小的噪声，验证了基于卡尔曼滤波的估计模型的性能。在均值为零、最大偏差分别为±0.4、±1.0和±2.5像素的噪声条件下，滤波后的坐标误差分别降低到原始坐标误差的37.5％、26.3％和20.7％左右。卡尔曼滤波估计的绝对距离误差随噪声增加而增大，但增大幅度小于滤波前增加的幅度。在均值为零、最大偏差为±2.5像素的噪声条件下，滤波前的绝对距离误差的均值和标准差分别为1.92像素和0.76像素，误差太大，无法进行正确的星图识别。而滤波后的平均值和标准差分别为0.40像素和0.20像素，仍然保持在亚像素级精度水平，从而能够实现正确的星图识别。随后，进行了夜场观星实验，验证了星点预测窗口的性能。以单帧星图为例，它的噪声水平相对较高。当采用经典的质心方法时，共获得了4985个目标对象。其中，候选星点对象数为572(占总数的11.5％)，严重阻碍了真实星点对象的准确提取。此外，由于噪声水平高，必须将全局阈值设置为128，这将过滤掉真实星点的弱像素信号，导致质心精度下降。相比之下，当使用星点预测窗口进行质心提取时，可以在一个非常小的预测窗口内自适应地设置局部阈值，共提取了8×5＝40个星点对象，仅占总数的0.8％。同理，在相同的实验条件下，共获得100张星图，有效星点对象数在每帧星图总目标对象数的平均占比仅为0.95％。仿真结果和实验结果均验证了该方法的可行性和有效性。
以上内容是结合具体优选实施方式对本发明做进一步详细说明，对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。
Claims (10)
1.一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，包括以下步骤：
结合前一帧星图中恒星星点的质心，以及星敏感器的全运动参数，生成星点预测窗口，并利用星点预测窗口对当前帧星图进行处理，获得星点图像；
采用质心法计算所述星点图像的质心，获得含噪声质心；
使用卡尔曼滤波器对所述含噪声质心中包含的图像噪声进行滤波，获得去噪质心；
通过所述去噪质心估计当前的姿态矩阵和全运动参数，作为最终输出信息；
使用所述最终输出信息生成下一时刻星点预测窗口，通过所述下一时刻星点预测窗口进行下一时刻星点质心的提取，通过循环往复所述下一时刻星点质心的提取过程，输出星点符合阈值的质心。
2.如权利要求1所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，结合前一帧星图中恒星星点的质心，以及全运动参数，生成下一时刻星点预测窗口，包括如下步骤：
构建多曝光成像星敏感器的坐标系；
在所述坐标系中对恒星i进行处理，获得其在坐标系中的观测单位矢量；
获取星敏感器在时刻tk的姿态矩阵，并通过姿态矩阵获取所述观测单位矢量与参考单位矢量之间的关系；
通过所述观测单位矢量与参考单位矢量之间的关系，获得星敏感器在不同时刻的旋转矩阵；
在时间间隔小于阈值时，将旋转矩阵近似值、观测单位矢量与时间之间的关系结合，生成下一时刻星点预测窗口。
3.如权利要求2所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，在所述坐标系中对恒星i进行处理，获得其在坐标系中的观测单位矢量，包括如下步骤：
获取恒星i的赤经αi和赤纬δi，并获得赤经αi和赤纬δi在坐标系中的参考单位矢量Vi；
将恒星i进行多曝光成像，设置曝光次数为N，获得N个时刻的恒星i的成像星点图像；
通过所述星点图像获取坐标系中恒星i的观测单位矢量；
其中，xi(tk)与yi(tk)分别为星点在时刻tk的坐标，f为星敏感器的焦距。
4.如权利要求3所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，所述获取星敏感器在时刻tk的姿态矩阵，并通过姿态矩阵获取所述观测单位矢量与参考单位矢量之间的关系，包括如下步骤：
令A(tk)为星敏感器在tk时刻的姿态矩阵；
通过所述姿态矩阵表示所述观测单位矢量与参考单位矢量之间的关系Wi(tk)，具体表达式为：
Wi(tk)＝A(tk)Vi(k＝1,2,...,N；i-1,2,...,M) (3)
其中，M表示在多曝光成像星敏感器的视场中导航星的总数；
以同样方式获取tk+1时刻的观测单位矢量与参考单位矢量之间的关系Wi(tk+1)；
Wi(tk+1)＝A(tk+1)Vi (4)
将tk时刻的观测单位矢量与参考单位矢量之间的关系表达式代入tk+1时刻的表达式，获得时刻tk、时刻tk+1中观测单位矢量之间的关系表达式为：
Wi(tk+1)＝A(tk+1)Vi＝A(tk+1)A(tk)TWi(tk) (5)
其中，A(tk+1)为星敏感器在tk+1时刻的姿态矩阵。
5.如权利要求4所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，通过所述观测单位矢量与参考单位矢量之间的关系，获得星敏感器在不同时刻之间的旋转矩阵，包括如下步骤：
定义R(tk,tk+1)作为多曝光成像星敏感器从时刻tk到tk+1的旋转矩阵，具体表达式为：
R(tk,tk+1)＝A(tk+1)A(tk)T (6)。
6.如权利要求5所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，所述在时间间隔小于阈值时，将旋转矩阵近似值、观测单位矢量与时间之间的关系结合，生成下一时刻星点预测窗口，包括如下步骤：
令时刻tk和时刻tk+1之间的时间间隔小于阈值，此时姿态矩阵满足：
R(tk,tk+1)＝I+[[Δθ]]+O(|Δθ|2) (7)
其中，I为3×3阶单位矩阵，Δθ＝[ΔθxΔθyΔθz]T是最小旋转角矢量，[[Δθ]]和O(|Δθ|2)分别为Δθ的交叉乘积矩阵和高阶无穷小矩阵；θx、θy、θz分别为x、y、z轴上的旋转角；
在时间间隔Δt＝tk+1-tk小于阈值时，将多曝光成像星敏感器的运动看作是均匀角加速度运动；
令ω(tk)＝[ωx(tk) ωy(tk) ωz(tk)]T和α(tk)＝[αx(tk) αy(tk) αz(tk)]T分别表示星敏感器在时刻tk时的角速度和角加速度，且角速度和角加速度为全运动参数，则瞬时角速度ω(t)＝[ωx(tk) ωy(tk) ωz(tk)]T表示为：
ω(t)＝ω(tk)+α(tk)(t-tk)(tk≤t≤tk+1) (8)
对公式(8)进行积分，获得最小旋转角矢量；
忽略高阶无穷小矩阵O(|Δθ|2)时，旋转矩阵R(tk,tk+1)近似为：
将公式(10)与观测单位矢量同时代入时刻tk、时刻tk+1中观测单位矢量之间的关系表达式，获得：
通过公式(11)，将观测单位矢量的第三分量之间关系表示为：
并通过公式(12)推导出：
将公式(13)带入单位观测矢量之间关系表达式中，转换为：
由于旋转角度小于阈值，最大计算结果小于0.001时进行忽略，公式(14)近似为：
考虑到焦距f是常数，并且(xi(tk),yi(tk))是成像平面中星点的坐标，公式(15)相应地重写为：
并简化为二维公式：
通过最终的二维公式生成下一时刻星点预测窗口。
7.如权利要求1所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，所述采用质心法计算所述星点的质心，获得含噪声质心，包括如下步骤：
获取星点图像的有效像素数据，使用每个预测窗口内星点的像素数据来计算真实星点的含噪声质心；
其中，计算真实星点的含噪声质心时采用质心方法，具体表达式为：
其中，m、n、f(m,n)为预测窗口中每个有效像素的x坐标、y坐标和灰度值，(xic(tk+1),yic(tk+1))是计算出的含噪声质心。
8.如权利要求7所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，使用卡尔曼滤波器对所述含噪声质心中包含的图像噪声进行滤波，获得去噪质心，包括如下步骤：
获取卡尔曼滤波的状态方程；
Xk+1＝Dk+1，kXk+BkUk+Wk (19)
其中，Xk+1＝[xi(tk+1),yi(tk+1)]T和Xk＝[xi(tk),yi(tk)]T分别是时刻tk和时刻tk+1的状态向量，Φk+1,k为状态转移矩阵，其表示为：
其中BkUk是确定性驱动项，其表示为：
其中，f为星敏感器的焦距。
9.如权利要求8所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，通过质心方法的计算公式描述了tk+1时刻的量测向量，其中量测方程表达式为：
zk+1＝Hk+1Xk+1+Vk+1 (22)
其中，Zk+1＝[xic(tk+1),yic(tk+1)]T和Xk+1＝[xi(tk+1),yi(tk+1)]T分别是量测向量和状态向量，Hk+1＝I2×2是量测矩阵；Vk+1表示量测噪声，其协方差矩阵为：
Rk+1＝μI2×2(μ>0) (23)
其中μ与星图的噪声水平有关。
10.如权利要求9所述的一种基于卡尔曼滤波的星敏感器星点质心提取方法，其特征在于，通过状态方程和量测方程构建星点质心滤波估计模型，具体表达式为：
其中，k和k+1分别表示当前和下一次曝光时刻；和/>是状态向量的估计，它对应于两个曝光时刻k和k+1；/>和P(k+1)k分别为一步预测向量和一步预测误差矩阵；P(k+1)和K(k+1)为估计误差矩阵和滤波器增益矩阵；当给定初始值/>和P(0)时，星点质心矢量/>根据上述公式估计。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311603531.5A CN117422742A (zh) | 2023-11-28 | 2023-11-28 | 一种基于卡尔曼滤波的星敏感器星点质心提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311603531.5A CN117422742A (zh) | 2023-11-28 | 2023-11-28 | 一种基于卡尔曼滤波的星敏感器星点质心提取方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117422742A true CN117422742A (zh) | 2024-01-19 |
Family
ID=89530999
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311603531.5A Pending CN117422742A (zh) | 2023-11-28 | 2023-11-28 | 一种基于卡尔曼滤波的星敏感器星点质心提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117422742A (zh) |
-
2023
- 2023-11-28 CN CN202311603531.5A patent/CN117422742A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109029433B (zh) | 一种移动平台上基于视觉和惯导融合slam的标定外参和时序的方法 | |
| CN109631887B (zh) | 基于双目、加速度与陀螺仪的惯性导航高精度定位方法 | |
| CN107590777A (zh) | 一种星敏感器星点图像增强方法 | |
| CN106780699B (zh) | 一种基于sins/gps和里程计辅助的视觉slam方法 | |
| CN102607526B (zh) | 双介质下基于双目视觉的目标姿态测量方法 | |
| CN112082545B (zh) | 一种基于imu和激光雷达的地图生成方法、装置及系统 | |
| CN107796391A (zh) | 一种捷联惯性导航系统/视觉里程计组合导航方法 | |
| CN109724586B (zh) | 一种融合深度图和点云的航天器相对位姿测量方法 | |
| CN109029425A (zh) | 一种采用区域滤波的模糊星图复原方法 | |
| CN111798523A (zh) | 星相机在轨定标定姿及遥感影像几何定位方法、系统 | |
| CN106250898B (zh) | 一种基于尺度预测的图像局部区域特征提取方法 | |
| Wei et al. | Restoration of motion-blurred star image based on Wiener filter | |
| CN107942090B (zh) | 一种基于模糊星图提取光流信息的航天器角速度估计方法 | |
| CN104154932B (zh) | 基于emccd和cmos的高动态星敏感器的实现方法 | |
| CN112179373A (zh) | 一种视觉里程计的测量方法及视觉里程计 | |
| CN113483748B (zh) | 小天体柔性附着多节点相对位姿估计方法 | |
| CN111260736A (zh) | 一种空间相机内参在轨实时标定方法 | |
| CN105403886B (zh) | 一种机载sar定标器图像位置自动提取方法 | |
| CN110361001A (zh) | 一种用于空间碎片运动测量系统及标定方法 | |
| Ning et al. | Angular velocity estimation using characteristics of star trails obtained by star sensor for spacecraft | |
| Ning et al. | Spacecraft angular velocity estimation method using optical flow of stars | |
| CN117422742A (zh) | 一种基于卡尔曼滤波的星敏感器星点质心提取方法 | |
| CN116295363A (zh) | 一种星点快速提取与高精度定位方法 | |
| Yu et al. | A novel inertial-aided feature detection model for autonomous navigation in planetary landing | |
| Chen | Low-cost star tracker development with a laboratory simulation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication |