CN102116633B - 深空光学导航图像处理算法仿真验证方法 - Google Patents

Abstract

本发明提供一种借助仿真的功能来完成深空光学导航图像处理算法的仿真验证方法。该方法将获得的星点轨迹仿真图像、随机生成的一系列姿态矩阵以及带有误差的节点预测值进行图像处理；所述的图像处理是先对带有误差的节点预测值进行修正；然后以每一个修正的节点像素坐标为中心，在其周围取匹配模板，进行交叉相关，得到多个恒星目标之间的相对像素坐标，再与真实值进行比较评定图像处理的像素误差；然后计算相对角矩，根据相对角矩采用双三角形恒星识别方法进行背景恒星的识别，根据背景恒星的识别结果，由恒星向量形成计算矩阵，通过矩阵运算，计算出目标小行星方位；最后将计算出的目标小行星方位与真实的目标小行星方位比较即可知偏差。

Description

深空光学导航图像处理算法仿真验证方法

技术领域

本发明涉及光学自主导航的图像处理领域，具体涉及一种深空光学导航图像处理算法的仿真验证方法。

背景技术

深空飞行段不同于近地球轨道，光学导航所采用的导航星目标除了视场中的恒星之外，还选择了巡航段周围分时出现的小行星，由于探测星等很低，为9～12等星，需要敏感器长时间的曝光，曝光时间内探测器姿态的晃动影响明显，导航天体成像不再是单个星点，而是随着扰动改变在图像中的位置，形成轨迹线，对图像处理方法和仿真方法带来新的挑战，需要重新进行研究。

目前，星点图像仿真方法有很多，但是星迹图像仿真的验证系统研究未见相关报导。

发明内容

本发明的目的在于提供一种借助仿真的功能来完成深空光学导航图像处理算法的仿真验证方法。

实现本发明目的的技术方案：一种深空光学导航图像处理算法仿真验证方法，其先在小行星和探测器设定轨道位置已知的前提下，查阅星表数据，获得图像中星点的初始位置和相对角矩；然后随机生成一系列的姿态矩阵表征曝光时间内的敏感器姿态扰动，进行图像仿真获得星点轨迹仿真图像；其特征在于：该方法将获得的星点轨迹仿真图像、随机生成的一系列姿态矩阵以及带有误差的导航节点预估值进行图像处理；所述的图像处理是先对带有误差的导航节点预估值进行修正；然后以每一个修正的节点像素坐标为中心，在其周围取匹配模板，进行交叉相关，得到多个恒星目标之间的相对像素坐标，再与真实值进行比较评定图像处理的像素误差；然后计算相对角矩，根据相对角矩采用双三角形恒星识别方法进行背景恒星的识别，根据背景恒星的识别结果，由恒星向量形成计算矩阵，通过矩阵运算，计算出目标小行星方位；最后将计算出的目标小行星方位与真实的目标小行星方位比较即可知偏差。

如上所述的一种深空光学导航图像处理算法仿真验证方法，其所述的节点为星点轨迹仿真图像上运动变化剧烈、特征突出的点；根据星表数据和光学导航敏感器的初始姿态矩阵及姿态扰动矩阵得到节点像素坐标真值，节点像素坐标真值加上0～10之间的随机浮点数得到带有误差的节点预测值。

本发明的效果在于：

本发明是在目前没有成熟的敏感器可以提供图像的情况下，借助仿真的功能来完成图像算法的验证。统一采用导航系统提供的测量信息，无误差的信息提供给图像仿真部分，加上测量误差的作为预测信息提供给图像处理部分，然后图像处理算法在预测信息的前提下对仿真得到的图像进行处理，得到的计算结果和真实值比较评定算法误差。

本发明为星迹图像的仿真与处理提供了研究基础，并给出了验证模式。对图像处理算法的前期开发研究工作提出了借鉴模式，在整个工程研制的过程中保证了软件算法的先行研究，不受制于硬件研制的约束，将算法的整体性能评价从整个系统中剥离，有利于整机调试的误差分析。

附图说明

图1为仿真图像；

图2为交互相关匹配原理图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步描述。

本发明所述的一种深空光学导航图像算法处理仿真验证方法步骤如下：先在小行星和探测器设定轨道位置已知的前提下，查阅星表数据，获得图像中星点的初始位置和相对角矩；然后随机生成一系列的姿态矩阵表征曝光时间内的敏感器姿态扰动，进行图像仿真系获得星点轨迹仿真图像；再将获得的星点轨迹仿真图像、随机生成的一系列姿态矩阵以及带有误差的导航节点预估值进行图像处理；所述的图像处理是先对带有误差的导航节点预估值进行修正；然后以每一个修正的节点像素坐标为中心，在其周围取匹配模板，进行交叉相关，得到多个恒星目标之间的相对像素坐标，再与真实值进行比较评定图像处理的像素误差；然后计算相对角矩，根据相对角矩采用双三角形恒星识别方法进行背景恒星的识别，根据背景恒星的识别结果，由恒星向量形成计算矩阵，通过矩阵运算，计算出目标小行星方位；最后将计算出的目标小行星方位与真实的目标小行星方位比较即可知算法偏差。

下面结合具体实例进行描述。

步骤一：获得星点轨迹仿真图像

设定飞行器深空探测的轨道为已知，在日心黄道坐标系下，当前探测器的位置为：(1.02565885315093×10¹¹，-1.11396451349936×10¹¹，-1.67302798002×10⁸)，日心黄道坐标系到探测器坐标系的姿态旋转阵为

$\left[\begin{array}{ccc}0.75078599084622& -0.64689074293123& 0.13361423075012\\ 0& -0.20227860892816& -0.97932801673907\\ 0.66054552905084& 0.73526575541091& -0.15186794583113\end{array}\right],$

小行星的位置为(3.79101473705228×10¹¹，1.96420556491809×10¹¹，-6.3746405086921×10¹⁰)，星等为9.029。同时探测器在成像曝光时间内的姿态扰动共有10次(第一次为初始位置)，其旋转变换阵分别为：

姿态2

0.750708916 -0.646893753 0.134032070

-0.000024210 -0.202911210 -0.979197141

0.660633122 0.735088779 -0.152342915

姿态3

0.751079762 -0.646403472 0.134319550

0.000092203 -0.203346580 -0.979106817

0.660211468 0.735399700 -0.152669901

姿态4

0.751459820 -0.645959646 0.134328984

0.000211810 -0.203360782 -0.979103849

0.659778822 0.735785654 -0.152680636

姿态5

0.751685363 -0.645687443 0.134375745

0.000282835 -0.203431597 -0.979089120

0.659521823 0.736004967 -0.152733931

姿态6

0.751503818 -0.646012253 0.133828927

0.000225843 -0.202602551 -0.979261025

0.659728703 0.735948623 -0.152110693

姿态7

0.751547218 -0.646017081 0.133561637

0.000239436 -0.202197352 -0.979344768

0.659679257 0.736055815 -0.151806177

姿态8

0.751597019 -0.645915991 0.133770155

0.000255012 -0.202513481 -0.979279442

0.659622510 0.736057622 -0.152043813

姿态9

0.751668989 -0.645869250 0.133591329

0.000277725 -0.202242133 -0.979335511

0.659540487 0.736173235 -0.151839766

姿态10

0.751809798 -0.645749746 0.133376514

0.000321913 -0.201916158 -0.979402757

0.659379955 0.736367524 -0.151594671

统一采用依巴谷星表数据，读取星表中恒星的赤经、赤纬，换算恒星矢量，再转换到日心黄道坐标系下。假设小行星的空间矢量即为敏感器的光轴矢量，求取恒星和小行星矢量间的夹角，小于0.8度视场的恒星矢量，认为会在相机的靶面上成像；为了模拟曝光时间内的星点轨迹，随机生成一系列的姿态矩阵表征曝光时间内的敏感器姿态扰动，进行图像仿真，获得星点轨迹仿真图像，如图1所示。

步骤二：获得带有误差的导航节点预估值

光学导航敏感器的初始姿态矩阵已知，若其视轴在天体坐标系中的指向为赤经赤纬(αz，δz)，其X轴指向为(αx，δx)，则从天球坐标系O-UVW(坐标原点O位于地球中心，U轴指向春分点，W轴为地球自旋轴)到光学导航敏感器坐标系O-XYZ的变换矩阵M为：

$\left(\begin{array}{c}u\\ v\\ w\end{array}\right)=M\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)$        式1)

其中

a11＝cosαxcosδx，

a12＝sinαxcosδx，

a13＝sinδx，

a31＝cosαzcosδz，

a32＝sinαzcosδz，

a33＝sinδz。

又：

$y=z\×x=\left[\begin{array}{ccc}0& -a_{33}& a_{32}\\ a_{33}& 0& -a_{31}\\ -a_{32}& a_{31}& 0\end{array}\right]\left[\begin{array}{c}{a}_{11}\\ a_{12}\\ a_{13}\end{array}\right]=\left[\begin{array}{c}{a}_{13}{a}_{32}-a_{12}{a}_{33}\\ a_{11}{a}_{33}-a_{13}{a}_{31}\\ a_{12}{a}_{31}-a_{11}{a}_{32}\end{array}\right]$      式2)

则可得：

a21＝sinαzcosδzsinδx-sinδzsinαxcosδx；

a22＝sinδzcosαxcosδx-cosαzcosδzsinδx；

a23＝sinαxcosδxcosαxcosδz-sinαzcosδzcosαxcosδx；

此时矩阵M可确定下来；对视场中出现的星点R，假定它的赤经赤纬为(αi，δi)，在惯性坐标系中可表示为ri＝(cosαxcosδx，sinαxcosδx，sinδx)^T，由变换矩阵，则可知它在光学导航敏感器中的坐标为rs＝M^Tri；

假设星点在光学导航敏感器坐标系中的坐标为rs＝(X，Y，Z)，则星点在光学导航敏感器焦平面的位置坐标(x，y)为：x＝f×X/Z，y＝f×Y/Z，f为光学导航敏感器光学系统的焦距；对上式进行归一化处理，得x＝f/s×X/Z，y＝f/s′×Y/Z，式中s，s′分别为焦平面X Y方向上像元的大小。

在敏感器的每个姿态变化时刻，通过上述方法可得到节点位置在焦平面的像素坐标真实值，节点像素坐标真值加上0～10之间的随机浮点数得到带有误差的导航节点预估值；所述的节点为星点轨迹图像上运动变化剧烈、特征突出的点；如下表所示。

可以看出节点预测值的误差已经很大，有的已经超过了10个像素。

步骤三：将获得的星点轨迹仿真图像、随机生成的一系列姿态矩阵以及带有误差的导航节点预估值进行图像处理，其包括节点预测值修正、匹配和交叉相关以及角矩计算和星图识别。

(a)节点预测值修正(具体可参考申请人同期申请专利“星点轨迹图像的节点预测修正方法”)

根据星表数据和光学导航敏感器的初始姿态矩阵及姿态扰动矩阵，将导航星点从惯性坐标系转换为光学导航敏感器坐标系，得到导航星点在光学导航敏感器焦平面的位置坐标；再根据导航星点在光学导航敏感器焦平面的位置坐标模拟星点的静态成像；对上述静态成像的星点按照时间序列连接，得到带有节点的星点动态轨迹图像。

对得到的星点动态轨迹图像的每个节点进行运动是否突出的判断；当判断某节点时，设定该时刻节点为中心节点、其上一时刻和下一时刻的节点为相邻节点，利用这些节点连线之间的夹角进行判断：

假设中心节点坐标为(u₀，v₀)，相邻两节点分别为(u₁，v₁)和(u₂，v₂)，则以中心点为起始点，两个矢量描述为：a：(u₁-u₀)i+(v₁-v₀)j，b：(u₂-u₀)i+(v₂-v₀)j，则夹角为

当

时，该中心节点认为是满足要求的节点，thred阈值取值范围为[-1，1]；

以满足要求的节点为中心，取一定大小的窗口为节点模板。

根据根据导航测量值得到带有误差的节点预估位置，在真实星点轨迹图像上的节点预估位置处取窗口，对窗口内轨迹进行判断，将满足亮度及宽度要求的节点进行保留；

然后利用得到的节点模板与上述保留的位置窗口进行粗匹配，修正节点位置，R_ij出现最大值的地方就是最佳节点；

$R_{\mathrm{ij}}=\frac{\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}[S_{\mathrm{ij}}(m,n)\×B(m,n)]}{\sqrt{\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left[S_{\mathrm{ij}}(m,n)\right]}^2}\sqrt{\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\left[B(m,n)\right]}^2}}$       式4)

式中：B表示节点模板矩阵，S_ij表示保留的位置窗口区域，i、j表示保留的位置窗口的行和列。

节点预测值修正得到的是修正后的节点坐标，如下表所示：

(b)匹配和交叉相关

以每一个修正的节点坐标为中心，在其周围取匹配模板，进行交叉相关，原理如图2所示。图中A、B、C分别表示三个不同轨迹上选取的待匹配位置窗口，r^A表示A轨迹上选取的修正节点，即匹配节点模板的中心，r^B表示B轨迹上选取的修正节点，r^C表示C轨迹上选取的修正节点。C^A _B表示B轨迹上的节点模板在A上进行相关匹配得到的匹配中心。C^A _C表示C轨迹上的节点模板在A上进行相关匹配得到的匹配中心，其余类推。AC、BC、CC表示该位置实际的理想对应中心。a，b，c表示这三个中心偏离理想中心的偏移量。

由上述图2示可以建立如下关系：

$\left\{\begin{array}{c}{C^A}_B-r^A=a-b\\ {C^A}_C-r^A=a-c\\ {C^B}_A-r^B=b-a\\ {C^B}_C-r^B=b-c\\ {C^C}_A-r^C=c-a\\ {C^C}_B-r^C=c-b\end{array}\right.$        式5)

上述方程组中独立的方程为3个，可以求解a，b，c三个偏移量。得到偏移量后，在r^A、r^B、r^C上修正偏移向量，利用新得到的位置点作为新的特征节点，建立新的节点模板，重复上述交叉相关的过程，得到新的修正偏移量，如此迭代循环，直到位置偏移量小于设定的残差时，将该次的匹配中心作为该条轨迹图像上节点的最终位置输出。同时为了提高节点的精度，保留最大相关的位置值和其左右次相关的位置值，利用三点进行二次曲线插值，求取曲线的峰值点坐标，作为最佳匹配位置输出。

经过节点预测修正、匹配和交叉相关，最后得到的多个恒星目标之间的相对像素坐标如下，与真实值进行比较可以评定图像处理的像素误差。

相对位置	计算得到的结果(单位：像素)	仿真的正确结果(单位：像素)	误差
				目标1与目标2	354.0594，-194.6563	354.08714，194.57557	0.02778，0.08072
目标1与目标3	-215.1823，-248.005	-215.07764，248.14098	0.10468，-0.13583
				目标1与目标4	225.8439，-474.7233	226.02030，474.54321	0.17641，0.18010
目标1与目标5	-122.6649，-117.906	-122.68744，118.03647	-0.02249，-0.13019
				目标1与目标6	无效	-18.9223，236.3419	无效
目标1与目标7	147.5030，-186.2485	147.56000，-186.15533	0.05704，0.09314
				目标2与目标3	-569.2417，-53.3489	-569.18938，-53.42043	0.05230，-0.07157
目标2与目标4	-128.2155，-280.0670	-128.17489，-280.1404	0.04058，-0.07334
				目标2与目标5	-476.7243，76.7500	-476.73751，76.52936	-0.01319，-0.22065
目标2与目标6	-372.9816，430.9982	-373.08854，430.8791	-0.10689，-0.11914
				目标2与目标7	-206.5564，8.4078	-206.5104，8.4418	0.045980.03395
目标3与目标4	441.0262，-226.7182	441.01994，-226.62522	-0.00627，0.09294
				目标3与目标5	92.5174，130.0989	92.45187，129.94978	-0.06550，-0.14908
目标3与目标6	196.2600，484.3471	196.09939，484.32699	-0.16064，-0.02007
				目标3与目标7	362.6853，61.7567	362.63353，61.95848	-0.05175，0.20180

目标4与目标5	-348.5088，356.8170	-348.52193，356.8050	-0.01308，-0.01205
				目标4与目标6	-244.7662，711.0652	-244.8912，711.12070	-0.12502，0.05548
目标4与目标7	-78.3409，288.4748	-78.36226，288.50712	-0.02132，0.03228
				目标5与目标6	103.7427，354.2482	103.63073，354.31572	-0.11194，0.06753
目标5与目标7	无效	270.1679，-68.3422	无效
				目标6与目标7	166.4252，-422.5904	166.58547，-422.45374	0.16022，0.13664

(c)计算相对角矩

根据上述结果计算相对角矩，已知某一星点在图像中的坐标为u1，v1(中心坐标表示)，另一星点在图像中的坐标为u2，v2，相机焦距为f。则两点之间的角矩(夹角的余弦值)计算公式如下：

   式6)

取f＝73337个像素长度，得到的角矩值如下表所示：

目标之间	角矩值
		目标1相对目标2	0.9999848220
目标1相对目标3	无效
		目标1相对目标4	0.9999743274
目标1相对目标5	0.9999973170
		目标1相对目标6	0.9999947647
目标1相对目标7	0.9999947592
		目标2相对目标3	0.9999696575
目标2相对目标4	0.9999911834
		目标2相对目标5	0.9999783117
目标2相对目标6	0.9999698013
		目标2相对目标7	0.9999960293
目标3相对目标4	0.9999771962
		目标3相对目标5	0.9999976306
目标3相对目标6	0.9999746009
		目标3相对目标7	0.9999874385

目标4相对目标5	0.9999768613
		目标4相对目标6	0.9999474449
目标4相对目标7	0.9999916913
		目标5相对目标6	0.9999873401
目标5相对目标7	无效
		目标6相对目标7	0.9999808301

(d)计算出目标小行星在日心黄道坐标系下的方位

根据相对角矩采用双三角形恒星识别方法进行背景恒星的识别，根据背景恒星的识别结果，由恒星向量形成计算矩阵，通过矩阵运算，计算出目标小行星在日心黄道坐标系下的方位。

设有m≥3颗恒星被识别，这些恒星在日心黄道坐标系中的单位矢量分别为U_i＝(x′_i，y′_i，z′_i)^T，它们与目标小行星的角距为cosθ_i，记目标小行星在日心黄道坐标系中的单位矢量为A＝(X，Y，Z)^T

则有：

由最小二乘法求解A。

得到的小行星空间指向为：(0.6605456697，0.7352658062，-0.1518671045)，其真实指向为：(0.66054552905084，0.73526575541091，-0.15186794583113)，计算值与真值的偏差为0.176189095267495850”。

Claims (1)

1.一种深空光学导航图像处理算法仿真验证方法，其先在小行星和探测器设定轨道位置已知的前提下，查阅星表数据，获得图像中星点的初始位置和相对角矩；然后随机生成一系列的姿态矩阵表征曝光时间内的敏感器姿态扰动，进行图像仿真获得星点轨迹仿真图像；根据星表数据和光学导航敏感器的初始姿态矩阵及姿态扰动矩阵得到节点像素坐标真值，节点像素坐标真值加上0～10之间的随机浮点数得到带有误差的节点预估值；其特征在于：该方法将获得的星点轨迹仿真图像、随机生成的一系列姿态矩阵以及带有误差的导航节点预估值进行图像处理；所述的图像处理是先对带有误差的导航节点预估值进行修正；然后以每一个修正的节点像素坐标为中心，在其周围取匹配模板，进行交叉相关，得到多个恒星目标之间的相对像素坐标，再与真实值进行比较评定图像处理的像素误差；然后计算相对角矩，根据相对角矩采用双三角形恒星识别方法进行背景恒星的识别，根据背景恒星的识别结果，由恒星向量形成计算矩阵，通过矩阵运算，计算出目标小行星方位；最后将计算出的目标小行星方位与真实的目标小行星方位比较即可知算法偏差；所述的节点为星点轨迹仿真图像上运动变化剧烈、特征突出的点。

Images