CN111780745B - 一种面向深空探测光学导航的短弧椭圆拟合优化方法 - Google Patents

Abstract

本发明公开了一种面向深空探测光学导航的短弧椭圆拟合优化方法，步骤如下：通过星历查询获得整个行星目标边缘的投影椭圆和模拟理论成像弧段；利用最小二乘椭圆拟合算法拟合上述获得的整个行星目标边缘的投影椭圆以获取拟合方程，从拟合方程的参数中提取椭圆的形状信息；获得带形状先验的椭圆标准方程，用列文伯格‑马夸尔特方法对该带形状先验的椭圆标准方程进行优化，以实现对模拟的相机实际成像弧段的拟合。本发明通过引入星历查询的方法，将行星目标边缘投影椭圆的形状先验信息引入椭圆的标准方程，通过带有先验信息的优化方法实现对带有误差的短弧的高精度拟合。

Description

一种面向深空探测光学导航的短弧椭圆拟合优化方法

技术领域

本发明属于导航制导与控制技术领域，具体涉及一种面向深空探测光学导航的短弧椭圆拟合优化方法。

背景技术

在深空探测任务中，探测器的导航系统是决定探测任务能否成功的关键因素之一。目前的深空探测任务，大都以地面测控的方式对探测器进行测控与导航，但是该方法成本高且通讯延时大。随着未来的深空探测任务变得更加复杂和具有挑战性，地面测控的方式将难以满足探测器的导航和控制需求。而自主导航对于深空探测而言是一种独立于地面测控的理想导航方法，探测器的自主导航不依赖地面深空网测控系统的支持，利用探测器自身搭载的测量手段和导航控制计算系统测定探测器自身位置速度等状态信息。而其中基于光学成像的自主导航技术因为其成本低、易实现等优点，成为了自主导航中最有希望的一种，并且已经在许多深空探测任务中做了尝试。例如，1998年，美国深空一号(Deep Space1)探测器利用搭载的相机对小行星进行拍照，并结合得到图像中的恒星背景实现了深空探测器巡航阶段中的自主导航。2005年，日本隼鸟号小行星探测器根据预先投掷在小行星上的导航信标，融合使用光学敏感器、雷达测距仪及激光测距仪进行自主导航，成功实现了与Itokawa小行星的交会和自主着陆。

在光学自主导航中，一步很重要的工作是对提取的目标边缘做椭圆拟合。绝大部分目前已有的研究中都采用了最小二乘算法来实现对目标边缘的拟合，在一些工作中甚至直接采用了圆拟合。最小二乘算法计算简单而且容易在工程上实现，但是在深空探测任务中，提取出的目标的边缘通常不是完整的，甚至只有一小段弧度，并且提取的边缘数据还存在误差。最小二乘算法在拟合弧度较小并且拟合数据带误差时会产生较大的椭圆拟合误差甚至会直接失效，这极大地影响了光学自主导航的可用范围和精度。

发明内容

针对于上述现有技术的不足，本发明的目的在于提供一种面向深空探测光学导航的短弧椭圆拟合优化方法，以解决现有技术在深空探测光学导航中，因提取的目标边缘弧度过小且带有误差而导致的椭圆拟合精度较低甚至直接失效的问题。本发明通过引入星历查询的方法，将行星目标边缘投影椭圆的形状先验信息引入椭圆的标准方程，通过带有先验信息的优化方法实现对带有误差的短弧的高精度拟合。

为达到上述目的，本发明采用的技术方案如下：

本发明的一种面向深空探测光学导航的短弧椭圆拟合优化方法，步骤如下：

1)通过星历查询结合航天器的自身位姿，获得整个行星目标边缘的投影椭圆和模拟理论成像弧段；

2)利用最小二乘椭圆拟合算法拟合上述步骤1)中获得的整个行星目标边缘的投影椭圆以获取拟合方程，从拟合方程的参数中提取椭圆的形状信息；

3)将上述形状信息引入椭圆的标准方程，获得带形状先验的椭圆标准方程，用列文伯格-马夸尔特方法对该带形状先验的椭圆标准方程进行优化，以实现对模拟的相机实际成像弧段的拟合。

进一步的，所述步骤1)中整个行星目标边缘的投影椭圆获得过程如下：

通过查询星历获得行星目标在惯性坐标系下表面的坐标P＝[X_w,Y_w,Z_w]，通过投影成像原理把行星目标表面坐标投影到相机的像素坐标系上，得到点P＝[u,v]^T；相机的投影原理如下式所示：

其中，是相机的内参数矩阵，/>是相机的外参数矩阵，其中R是旋转矩阵，T是平移向量。

进一步的，所述步骤1)具体还包括：根据相机视角在获得的整个行星目标的投影椭圆上截取可见的弧段，即为模拟理论成像弧段。

进一步地，所述步骤2)具体包括：

椭圆的一般方程为：

其中，(x_i,y_i)是椭圆上的点，为椭圆一般方程的自变量组成的矩阵，a，b，c，d，f_e为椭圆一般方程的系数，为方便表示将其写为矩阵形式a＝[a,b,c,d,e,f_e]^T；且方程的系数满足约束：

4ac-b²＝1

该约束写成矩阵的形式：

引入了一个拉格朗日乘数λ₁的最小二乘拟合的约束椭圆拟合问题的目标函数为：

min J＝||D·a||²+λ₁(1-a^TCa) (1)

其中：是椭圆弧段上所有点组成的自变量的矩阵；将拟合问题视为秩不足的广义特征值问题：

(D^TD)a＝λ₁Ca (2)

把上式(2)代入式(1)，得：

J＝a^TD^TDa＝λ₁a^TCa＝λ₁

避免D^TD的奇异，修改a、D和C如下：

a＝[a₁,a₂]^T，a₁＝[a,b,c]^T，a₂＝[d,e,f_e]^T

定义矩阵：

得到：

S₁a₁+S₂a₂＝λ₁C₁a₁

最后得到a₁的解如下：

通过以下给出的标准椭圆参数来计算完整行星目标的投影椭圆的形状信息：

其中，A和B分别是椭圆的长轴和短轴。

进一步地，所述步骤3)具体包括：

31)将上述步骤2)中拟合获得的椭圆的形状信息引入椭圆的标准方程，椭圆的形状由椭圆的长轴和短轴之比确定，即k＝A/B；椭圆的一般方程中五个系数均与长轴、短轴有关，无法将椭圆的形状信息代入，故选择对椭圆的标准方程做相应的变换，变换后的椭圆标准方程如下：

其中，C为椭圆中心的x坐标，D为椭圆中心的y坐标，θ为椭圆的倾角；

上述变换后的椭圆标准方程无法构造最小二乘，故构造一个均方误差函数，通过优化的方式去拟合求解C、D、A、θ四个参数；

32)优化椭圆的标准方程，其目标函数为：

优化问题的本质是一个非线性最小二乘问题：

其中，x＝(A,C,D,θ)，是由椭圆标准方程四个参数组成的矩阵；

将f(x)进行一阶泰勒展开：

f(x+Δx)≈f(x)+J(x)Δx

式中，J(x)是f(x)关于x的导数，为一个雅可比矩阵；当前的目标是寻找一个下降矢量Δx，使得||f(x+Δx)||²达到最小；为求Δx，需要求解一个线性的最小二乘问题：

在该线性最小二乘问题中引入信赖区域，则问题变为求解：

式中，μ是信赖区域的半径，D取成单位阵I；上式(3)为一个带不等式约束的优化问题，用拉格朗日乘子将其转化为一个无约束优化问题；

式中，λ₂为拉格朗日乘子，将上式(4)展开，求其关于Δx的导数，并令其为零，得到如下方程：

(J(x)^TJ(x)+λ₂I^TI)Δx＝-J(x)^Tf(x)

令：

H＝J(x)^TJ(x)，g＝-J(x)^Tf(x)

可得：

(H+λ₂I)Δx＝g

由此求得下降矢量Δx；

33)确定信赖区域的范围，使用如下公式(5)来判断泰勒近似是否够好，公式如下：

式中，f(x+Δx)-f(x)实际函数下降值，J(x)Δx是近似函数下降值；若ρ接近于1，则近似是好的；若则认为近似差，需要缩小近似范围；反之，若/>则说明实际下降的比预计下降的大，则需要放大近似范围。

本发明的有益效果：

1、本发明利用星历查询获得行星目标的投影椭圆，通过最小二乘椭圆拟合算法能够获得相机获取的行星目标实际成像弧段的形状先验信息；

2、本发明将获得的形状先验信息引入椭圆的标准方程，然后用列文伯格-马夸尔特方法实现对带先验信息的椭圆标准方程的优化，能够在不显著增加计算代价的前提下提高椭圆拟合的精度，并且在传统最小二乘算法不能使用的小弧段拟合中依然可以使用，扩大了光学导航在深空探测中的适用范围。

附图说明

图1为本发明方法的流程图；

图2为航天器对火星的模拟投影成像示意图；

图3a为拟合弧段高斯白噪声标准差为0.1时两种算法的拟合结果对比图；

图3b为拟合弧段高斯白噪声标准差为0.5时两种算法的拟合结果对比图；

图3c为拟合弧段高斯白噪声标准差为1时两种算法的拟合结果对比图；

图4a为椭圆长轴的迭代结果示意图；

图4b为椭圆中心x坐标的迭代结果示意图；

图4c为椭圆中心y坐标的迭代结果示意图；

图4d为椭圆倾斜角的迭代结果示意图；

图5a为不同拟合弧度时两种算法拟合结果椭圆中心x坐标误差对比图；

图5b为不同拟合弧度时两种算法拟合结果椭圆中心y坐标误差对比图；

图6a为不同干扰误差情况下本发明算法拟合结果椭圆中心x坐标误差示意图；

图6b为不同干扰误差情况下本发明算法拟合结果椭圆中心y坐标误差示意图。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

参照图1所示，本发明的一种面向深空探测光学导航的短弧椭圆拟合优化方法，步骤如下：

1)通过星历查询的方法，结合航天器的自身位姿，获得整个行星目标边缘的投影椭圆和模拟理论成像弧段；

所述步骤1)中整个行星目标边缘的投影椭圆获得过程如下：

通过查询星历获得行星目标在惯性坐标系下表面的坐标P＝[X_w,Y_w,Z_w]，通过投影成像原理把行星目标表面坐标投影到相机的像素坐标系上，得到点P＝[u,v]^T；相机的投影原理如下式所示：

其中，是相机的内参数矩阵，/>是相机的外参数矩阵，其中R是旋转矩阵，T是平移向量。

根据相机视角在获得的整个行星目标的投影椭圆上截取可见的弧段，即为模拟理论成像弧段。

2)利用最小二乘椭圆拟合算法拟合上述步骤1)中获得的整个行星目标边缘的投影椭圆以获取拟合方程，从拟合方程的参数中提取椭圆的形状信息；

椭圆的一般方程为：

其中，(x_i,y_i)是椭圆上的点，为椭圆一般方程的自变量组成的矩阵，a，b，c，d，f_e为椭圆一般方程的系数，为方便表示将其写为矩阵形式a＝[a,b,c,d,e,f_e]^T；且方程的系数满足约束：

4ac-b²＝1

该约束写成矩阵的形式：

引入了一个拉格朗日乘数λ₁的最小二乘拟合的约束椭圆拟合问题的目标函数为：

min J＝||D·a||²+λ₁(1-a^TCa) (1)

其中：是椭圆弧段上所有点组成的自变量的矩阵；

将拟合问题视为秩不足的广义特征值问题：

(D^TD)a＝λ₁Ca (2)

把上式(2)代入式(1)，可得：

J＝a^TD^TDa＝λ₁a^TCa＝λ₁

避免D^TD的奇异，修改a、D和C如下：

a＝[a₁,a₂]^T，a₁＝[a,b,c]^T，a₂＝[d,e,f_e]^T

定义矩阵：

得到：

S₁a₁+S₂a₂＝λ₁C₁a₁

最后得到a₁的解如下：

通过以下给出的标准椭圆参数来计算完整行星目标的投影椭圆的形状信息：

其中，A和B分别是椭圆的长轴和短轴。

3)将上述形状信息引入椭圆的标准方程，获得带形状先验的椭圆标准方程，用列文伯格-马夸尔特方法对该带形状先验的椭圆标准方程进行优化，以实现对模拟的相机实际成像弧段的拟合；

31)将上述步骤2)中拟合获得的椭圆的形状信息引入椭圆的标准方程，椭圆的形状由椭圆的长轴和短轴之比确定，即k＝A/B；椭圆的一般方程中五个系数均与长轴、短轴有关，无法将椭圆的形状信息代入，故选择对椭圆的标准方程做相应的变换，变换后的椭圆标准方程如下：

其中，C为椭圆中心的x坐标，D为椭圆中心的y坐标，θ为椭圆的倾角；

上述变换后的椭圆标准方程无法构造最小二乘，故构造一个均方误差函数，通过优化的方式去拟合求解C、D、A、θ四个参数；

32)优化椭圆的标准方程，其目标函数为：

优化问题的本质是一个非线性最小二乘问题：

其中，x＝(A,C,D,θ)，是由椭圆标准方程四个参数组成的矩阵；

将f(x)进行一阶泰勒展开：

f(x+Δx)≈f(x)+J(x)Δx

式中，J(x)是f(x)关于x的导数，为一个雅可比矩阵；当前的目标是寻找一个下降矢量Δx，使得||f(x+Δx)||²达到最小；为求Δx，需要求解一个线性的最小二乘问题：

在该线性最小二乘问题中引入信赖区域，则问题变为求解：

式中，μ是信赖区域的半径，D取成单位阵I；上式(3)为一个带不等式约束的优化问题，用拉格朗日乘子将其转化为一个无约束优化问题；

式中，λ₂为拉格朗日乘子，将上式(4)展开，求其关于Δx的导数，并令其为零，得到如下方程：

(J(x)^TJ(x)+λ₂I^TI)Δx＝-J(x)^Tf(x)

令：

H＝J(x)^TJ(x)，g＝-J(x)^Tf(x)

可得：

(H+λ₂I)Δx＝g

由此求得下降矢量Δx；

33)确定信赖区域的范围，使用如下公式(5)来判断泰勒近似是否够好，公式如下：

式中，f(x+Δx)-f(x)实际函数下降值，J(x)Δx是近似函数下降值；若ρ接近于1，则近似是好的；若则认为近似差，需要缩小近似范围；反之，若/>则说明实际下降的比预计下降的大，则需要放大近似范围。

以下通过具体实例说明算法的使用流程：

设定如下计算条件和技术参数：

(1)以火星探测为计算背景，考虑火星探测的环绕段任务，以STK数据作为火星的星历数据，并将火星假设为一个球；

(2)采用了一条航天器环绕火星的轨道，轨道的轨道六根数为：偏心率：0.61545319，升交点赤经：286.419°，近火点幅角：320.359°，真近点角：0°，半长轴：9501.9km，轨道倾角：86.9°；

(3)相机的视角设为6°×8°；

(4)在模拟理论成像弧段中添加旋转、平移、缩放以及高斯白噪声以模拟实际成像弧段；

采用Matlab软件进行仿真验证，可以获得图2-图6b的仿真结果。图2是航天器对火星的模拟投影成像示意图，模拟理论成像弧段的弧度为71.54°。图3a为拟合弧段高斯白噪声标准差为0.1时两种算法的拟合结果对比图；图3b为拟合弧段高斯白噪声标准差为0.5时两种算法的拟合结果对比图；图3c为拟合弧段高斯白噪声标准差为0.5时两种算法的拟合结果对比图；图4a为椭圆长轴的迭代结果示意图；图4b为椭圆中心x坐标的迭代结果示意图；图4c为椭圆中心y坐标的迭代结果示意图；图4d为椭圆倾斜角的迭代结果示意图；图5a为不同拟合弧度时两种算法拟合结果椭圆中心x坐标误差对比图；图5b为不同拟合弧度时两种算法拟合结果椭圆中心y坐标误差对比图；图6a为不同干扰误差情况下本发明算法拟合结果椭圆中心x坐标误差；图6b为不同干扰误差情况下本发明算法拟合结果椭圆中心y坐标误差。从图3a、图3b和图3c中可以看出，当拟合弧段的高斯白噪声变大时，最小二乘算法的精度急剧下降而本发明的优化方法的精度总体上可以保持较高的水平。从图4a-图4d可以看到参数很快收敛，迭代次数少，计算量不会特别大。从图5a和图5b中可以得出结论，本发明的优化方法可以在更小的弧度上实现椭圆拟合，相比于最小二乘拟合算法所需的最小的弧度，本发明的算法减少了将近40°，极大地扩大了深空探测光学导航算法的适用范围，体现了该方法的优势。图6a和图6b则验证了本发明算法的鲁棒性。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (5)

1.一种面向深空探测光学导航的短弧椭圆拟合优化方法，其特征在于，步骤如下：

1)通过星历查询结合航天器的自身位姿，获得整个行星目标边缘的投影椭圆和模拟理论成像弧段；

2)利用最小二乘椭圆拟合算法拟合上述步骤1)中获得的整个行星目标边缘的投影椭圆以获取拟合方程，从拟合方程的参数中提取椭圆的形状信息；

3)将上述形状信息引入椭圆的标准方程，获得带形状先验的椭圆标准方程，用列文伯格-马夸尔特方法对该带形状先验的椭圆标准方程进行优化，以实现对模拟的相机实际成像弧段的拟合。

2.根据权利要求1所述的面向深空探测光学导航的短弧椭圆拟合优化方法，其特征在于，所述步骤1)中整个行星目标边缘的投影椭圆获得过程如下：

通过查询星历获得行星目标在惯性坐标系下表面的坐标P＝[X_w,Y_w,Z_w]^T，通过投影成像原理把行星目标表面坐标投影到相机的像素坐标系上，得到点P′＝[u,v]^T；相机的投影原理如下式所示：

其中，是相机的内参数矩阵，/>是相机的外参数矩阵，其中R是旋转矩阵，T是平移向量。

3.根据权利要求1所述的面向深空探测光学导航的短弧椭圆拟合优化方法，其特征在于，所述步骤1)具体还包括：根据相机视角在获得的整个行星目标的投影椭圆上截取可见的弧段，即为模拟理论成像弧段。

4.根据权利要求1所述的面向深空探测光学导航的短弧椭圆拟合优化方法，其特征在于，所述步骤2)具体包括：

椭圆的一般方程为：

其中，(x_i,y_i)是椭圆上的点，为椭圆一般方程的自变量组成的矩阵，a，b，c，d，f_e为椭圆一般方程的系数，将其写为矩阵形式a＝[a,b,c,d,e,f_e]^T；且方程的系数满足约束：

4ac-b²＝1

该约束写成矩阵的形式：

引入了一个拉格朗日乘数λ₁的最小二乘拟合的约束椭圆拟合问题的目标函数为：

min J＝||D·a||²+λ₁(1-a^TCa) (1)

其中：是椭圆弧段上所有点组成的自变量的矩阵；将拟合问题视为秩不足的广义特征值问题：

(D^TD)a＝λ₁Ca (2)

把上式(2)代入式(1)，得：

J＝a^TD^TDa＝λ₁a^TCa＝λ₁

避免D^TD的奇异，修改a、D和C如下：

a＝[a₁,a₂]^T，a₁＝[a,b,c]^T，a₂＝[d,e,f_e]^T

D＝[D₁,D₂],

定义矩阵：

得到：

S₁a₁+S₂a₂＝λ₁C₁a₁

最后得到a₁的解如下：

通过以下给出的标准椭圆参数来计算完整行星目标的投影椭圆的形状信息：

其中，A和B分别是椭圆的长轴和短轴。

5.根据权利要求4所述的面向深空探测光学导航的短弧椭圆拟合优化方法，其特征在于，所述步骤3)具体包括：

31)将上述步骤2)中拟合获得的椭圆的形状信息引入椭圆的标准方程，椭圆的形状由椭圆的长轴和短轴之比确定，即k＝A/B；选择对椭圆的标准方程做相应的变换，变换后的椭圆标准方程如下：

其中，C为椭圆中心的x坐标，D为椭圆中心的y坐标，θ为椭圆的倾角；

构造一个均方误差函数，通过优化的方式去拟合求解C、D、A、θ四个参数；

32)优化椭圆的标准方程，其目标函数为：

优化问题是一个非线性最小二乘问题：

其中，x＝(A,C,D,θ)，是由椭圆标准方程四个参数组成的矩阵；

将f(x)进行一阶泰勒展开：

f(x+Δx)≈f(x)+J(x)Δx

式中，J(x)是f(x)关于x的导数，为一个雅可比矩阵；当前的目标是寻找一个下降矢量Δx，使得||f(x+Δx)||²达到最小；为求Δx，需求解一个线性的最小二乘问题：

在该线性最小二乘问题中引入信赖区域，则问题变为求解：

式中，μ是信赖区域的半径，D取成单位阵I；上式(3)为一个带不等式约束的优化问题，用拉格朗日乘子将其转化为一个无约束优化问题；

式中，λ₂为拉格朗日乘子，将上式(4)展开，求其关于Δx的导数，并令其为零，得到如下方程：

(J(x)^TJ(x)+λ₂I^TI)Δx＝-J(x)^Tf(x)

令：

H＝J(x)^TJ(x)，g＝-J(x)^Tf(x)

得：

(H+λ₂I)Δx＝g

由此求得下降矢量Δx；

33)确定信赖区域的范围，使用如下公式(5)来判断泰勒近似是否够好，公式如下：

式中，f(x+Δx)-f(x)实际函数下降值，J(x)Δx是近似函数下降值；若ρ接近于1，则近似是好的；若则认为近似差，需要缩小近似范围；反之，若/>则说明实际下降的比预计下降的大，则需要放大近似范围。