一种基于单星敏感器的星光折射卫星自主导航方法

摘要

本发明公开了一种基于单星敏感器的星光折射卫星自主导航方法，包括以下几个步骤：步骤一，按照最佳安装角度将星敏感器安装在卫星上；步骤二，星敏感器拍摄星图后，使用三角形算法识别星图中的正常星；步骤三，利用识别的正常星计算星敏感器光轴指向和卫星姿态；步骤四，根据星敏感器光轴指向从星表中选星生成模拟折射星图；步骤五，利用模拟折射星图识别折射星，根据识别结果计算星光折射角；步骤六，将星光折射角代入系统模型，星载计算机利用最优估计方法得到卫星的导航信息。本发明提高了星光折射卫星自主导航的精度、降低了设计成本。

说明书

技术领域

本发明涉及一种基于单星敏感器的星光折射卫星自主导航方法，属于卫星自主导航和星 图识别的技术领域。

背景技术

由于星光折射间接敏感地平的方法是一种低成本、高精度的自主导航方法。美国对于星 光折射法自主导航的研究工作可以追溯到60年代。在实施Apollo计划的过程中，就对利用 天体掩星、星光在大气中的折射、星光穿越大气时的衰减等实现自主导航的方案进行了研 究.1975年由美国海军研究局和美国国防部高级研究计划局共同投资，麻省理工学院Draper 实验室对星光折射/星光色散自主导航方案进行了研究和论证，结果表明在一个轨道周期可观 测40颗折射星的理想条件下，导航精度可以达到100m。90年代初投入使用的 MADAN(multi-mission attitude determination and autonomous navigation)导航系统(多 任务姿态确定和自主导航系统)便利用了星光折射原理。二十世纪80年代初期，法国也进行 了星光折射法自主导航的研究。1985年和1986年，CNES多次释放平流层气球对星光折射进 行了实际测量，在此基础上，对大气折射的精确模型、测量方案、自然环境对系统观测的约 束、误差分配和系统性能优化等方面进行了深入的分析和仿真试验，当时预计该系统导航精 度为300m。

基于星光折射间接敏感地平的方法是一种低成本、高精度的自主导航方法。传统的星光 折射方法中使用了两个星敏感器，一个用来敏感非折射星，另一个用来观测折射星，多个星 敏感器的使用不仅增加了设计成本，而且增加了初始校准中安装矩阵的校准负担。针对这种 情况，本发明公开了一种基于单星敏感器的星光折射卫星自主导航方法。由于星敏感器的安 装角度决定了卫星在一个周期内观测到的折射星数目，从而影响到导航的精度。因此，本发 明利用球面几何原理给出了一种计算星敏感器最佳安装角度范围的方法。在量测信息的获取 过程中，折射角的求取是最重要的一环，而折射星的识别是折射角计算的前提。本发明依据 折射星和非折射星的区别，给出了一种只使用一个星敏感器就可以进行折射星识别的方法。

发明内容

本发明的目的是为了提高星光折射卫星自主导航的精度、降低设计成本，提出了一种基 于单星敏感器的星光折射卫星自主导航方法

一种基于单星敏感器的星光折射卫星自主导航方法，包括以下几个步骤：

步骤一：按照最佳安装角度将星敏感器安装在卫星上；

步骤二：星敏感器拍摄星图后，使用三角形算法识别星图中的正常星；

步骤三：利用识别的正常星计算星敏感器光轴指向和卫星姿态；

步骤四：根据星敏感器光轴指向从星表中选星生成模拟折射星图；

步骤五：利用模拟折射星图识别折射星，根据识别结果计算星光折射角；

步骤六：将星光折射角代入系统模型，星载计算机利用最优估计方法得到卫星的导航信 息。

一种基于单星敏感器的星光折射卫星自主导航方法，还包括：

（1）星敏感器的最佳安装角度的范围是：

θ∈(a₄,a₃)

其中，

$a_4=\arccos \left[\frac{\cos \alpha }{\cos \left(\frac{{\theta }_{\mathrm{FOV}}}{2}\right)}\right],$$a_3=\arccos \left[\frac{\cos \beta }{\cos \left(\frac{{\theta }_{\mathrm{FOV}}}{2}\right)}\right]$

其中，θ为星敏感器的最佳安装角度，θ_FOV为星敏 感器视场大小，ha＝20km，hb＝50km，R_e为地球半径，r为卫星到地心的距离。

（2）模拟折射星图的具体生成方法是：

利用星敏感器的光轴指向和星敏感器参数从星表中选择落在星敏感器视场中星光矢量为 s的折射星，将折射星投影到星敏感器像平面上；

$\alpha -{\theta }_R<\arccos (\frac{r}{\left|r\right|}·s)<\beta$

其中，ha＝20km，R_e为地球半径，s为折射星的星光矢量，r为卫星 相对于地心的位置矢量，|r|为r模长，代表卫星到地心的距离，θ_R为折射高度为20km时对 应的星光折射角。

（3）星光折射角是：

$R=\arccos \left(S_{c1}^T{S}_{c2}\right)$

其中，

$S_{\mathrm{ci}}=\frac{1}{\sqrt{P_{\mathrm{xi}}^2+P_{\mathrm{yi}}^2+f^2}}\left(\begin{array}{c}-P_{\mathrm{xi}}\\ -P_{\mathrm{yi}}\\ f\end{array}\right),$i＝1或2

其中，f为星敏感器光学系统的角距，S_c1为折射星的星光在星敏感器坐标系下折射前的 单位矢量，S_c2为折射星的星光在星敏感器坐标系下折射后的单位矢量，(P_x1,P_y1)为折射星折 射前在像平面的位置坐标，(P_x2,P_y2)为折射星折射后在像平面的位置坐标。

（4）系统模型包括状态模型和量测方程，系统的状态模型是：

$\left(\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mu \frac{x}{r^3}[1-J_2\left(\frac{R_e}{r}\right)(7.5\frac{z^2}{r^2}-1.5)]+\Delta F_x\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mu \frac{y}{r^3}[1-J_2\left(\frac{R_e}{r}\right)(7.5\frac{z^2}{r^2}-1.5)]+\Delta F_y\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mu \frac{z}{r^3}[1-J_2\left(\frac{R_e}{r}\right)(7.5\frac{z^2}{r^2}-1.5)]+\Delta F_z\end{array}\right)$

$r=\sqrt{x^2+y^2+z^2}$

其中，系统状态为X＝[x y z v_x v_y v_z]^T，[x y z]^T代表卫星的位置矢量，[v_x v_y v_z]^T代 表卫星的速度矢量，μ＝3.986×10¹⁴m³/s²为地心引力常数，J₂＝0.00108263为地球引力系数， ΔF_x,ΔF_y,ΔF_z为地球非球形摄动的高阶摄动项和日、月摄动以及太阳光压摄动和大气摄动的 影响；

系统的量测方程是：

$Z=h_a+v=\left(\begin{array}{c}\sqrt{r^2-u_1^2}+u_1\tan R_1-R_e+v_1\\ \sqrt{r^2-u_2^2}+u_2\tan R_2-R_e+v_2\\ ·\\ ·\\ ·\\ \sqrt{r^2-u_i^2}+u_i\tan R_i-R_e+v_i\end{array}\right)$

其中，i代表观测到的折射星数目，R_e为地球半径，R₁,R₂...R_i均为星光折射角。

本发明的优点在于：

（1）本发明根据球面几何原理能够计算星敏感器的最佳安装角度范围，根据计算结果安 装星敏感器后使得卫星上安装的星敏感器在轨道周期中能够观测到更多的折射星，提高了导 航的精度；

（2）本发明提出了一种单星敏感器的折射星识别方法，相较于传统方法既降低了设计成 本又减轻了星敏感器在初始校准时的负担。

附图说明

图1是折射星的选择图；

图2是折射区域图；

图3是正常星识别星图；

图4是折射星识别星图；

图5是速度误差实验结果图；

图6是位置误差实验结果图；

图7是位置误差随轨道高度变化的实验结果图；

图8是一个轨道周期内观测到的折射星数目随轨道高度的变化结果图；

图9是一个轨道周期内观测到的折射星数目随安装角度的变化结果图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。本发明是一种单星敏感器的星 光折射卫星自主导航方法，包括以下几个步骤：

步骤一：按照最佳安装角度将星敏感器安装在卫星上；

卫星自主导航中，根据平流层的厚度折射高度一般选取20km-50km，即图1中ha=20km， hb=50km；假设星敏感器安装在卫星轨道平面内，星光的矢量为s，则由图1得符合方程（1） 的恒星星光产生折射后可以被卫星接收到。星光未发生折射的恒星为正常星，星光发生折射 的恒星为折射星。

$\alpha -{\theta }_R<\arccos (\frac{r}{\left|r\right|}·s)<\beta ---\left(1\right)$

其中，r为卫星相对于地心的位置矢量，|r|为r模长并且代表卫星到地心的距离，α和β 如图1所示，可表示为：

$\alpha =\arcsin \left(\frac{\mathrm{ha}+R_e}{r}\right);$$\beta =\arcsin \left(\frac{\mathrm{hb}+R_e}{r}\right);$

θ_R为折射高度为20km时对应的星光折射角，当卫星处在轨道上的某一点时，所有符合上 述公式的恒星投影到天球上就变成了一个圆环，如图2折射区域所示。圆环与星敏感器视场 的相交部分(如图2中阴影部分所示)就是能够被星敏感器观测到的折射星所在区域。其中，C 为星敏感器视轴在天球上的投影点，求取星敏感器的最佳安装角度θ实际上就是求C在天球 上的最佳位置。适当调整星敏感器的安装角度就能够使得相交区域的面积达到最大，这就是 求取星敏感器最佳安装角度的目的。

如图2所示A点是卫星沿着其位置矢量反方向在天球上的投影点；B,F,D,E分别为星 敏感器视场投影与折射圆环的交点。为相应的大圆弧,为视场 与外圆相交形成的圆弧，为视场与内圆相交形成的圆弧。由图1得，圆环的宽度对应的 角距离为θ_d＝(β-α+θ_R)，若取轨道高度为686km，得θ_d≈0.67°，相对于视场来说非常小， 在这种情况下可以通过如下的方法来得到最佳视轴投影范围：先求取得最大值时视轴位 置，再求取得最大值时视轴位置，将最佳视轴投影范围定位在二者之间。

在球面三角形ABE中，对应的球面角大小为2a，因此当取得最大的弧长时，角度 a也是最大的，我们可以通过求取角度a的最大值来得到相交弧长的最大值。对应角距离 为a₁，对应角距离为a₂，则由球面几何的余弦定理，在球面三角形ABC中有

cosa₂＝cosa₁cosa₃+sina₁sina₃cosa    (2)

其中，a₁,a₂,a₃为三角形三边的角距离，由图4得，a₁＝β，其中θ_FOV为星 敏感器视场大小，带入方程（2）得

$a=\arccos \left[\frac{\cos \left(\frac{{\theta }_{\mathrm{FOV}}}{2}\right)-\cos \beta \cos a_3}{\sin \beta \sin a_3}\right]---\left(3\right)$

对方程（3）求导，可得a取极大值时

$a_3=\arccos \left[\frac{\cos \beta }{\cos \left(\frac{{\theta }_{\mathrm{FOV}}}{2}\right)}\right]---\left(4\right)$

同理可以求得得到最大弧长时A点到视场中心的角距离为

$a_4=\arccos \left[\frac{\cos \alpha }{\cos \left(\frac{{\theta }_{\mathrm{FOV}}}{2}\right)}\right]---\left(5\right)$

则星敏感器的最佳安装角度的范围为θ∈(a₄,a₃)；以轨道高度686km，星敏感视场10°×10° 为例，可得θ∈(64.82°,65.40°)。而卫星的轨道高度一般都高于686km，随着轨道高度的增加 这个角度范围会更小，因此可以近似的取中点作为最佳的安装角度，此时令 ha+R_e＝m₁，hb+R_e＝m₂，如图2所示，则可得相交区域的面积为对应球 面角为2a的一部分圆环

$S=r_s^2·\frac{\sqrt{r^2-m_1^2}-\sqrt{r^2-m_2^2}}{r}·[\arccos \frac{\sqrt{m_2^2-r^2{m}_3^2}}{m_2}+\arccos \frac{\sqrt{m_1^2-r^2{m}_3^2}}{m_1}]---\left(6\right)$

其中，r_s为卫星到天球表面的半径，一般将整个星空等效为单位天球，在这种情况下r_s＝1。 假设天球上恒星分布均匀，全天范围内星敏感器能够观测到恒星数目的公式

N＝6.57e_1.08M    (7)

其中，M为星敏感器的星等敏感极限，利用此公式我们可以得到星敏感器能够观测到的 折射星数目

$\mathrm{starnum}=\frac{N·S}{4\pi r_s^2}$

将所有已知量带入得

$\mathrm{starnum}=\frac{6.57e^{1.08M}}{4\pi }·\frac{\sqrt{r^2-m_1^2}-\sqrt{r^2{m}_2^2}}{r}·[\arccos \frac{\sqrt{m_2^2-r^2{m}_3^2}}{m_2}+\arccos \frac{\sqrt{m_1^2-r^2{m}_3^2}}{m_1}]---\left(8\right)$

将式(6)对r求导得

$\frac{\mathrm{dS}}{\mathrm{dr}}=r_s^2·{\Omega }_1·{\Omega }_2<0---\left(9\right)$

其中，

${\Omega }_1=\frac{\sqrt{r^2-m_2^2}-\sqrt{r^2-m_1^2}}{r^2}+\frac{1}{\sqrt{r^2-m_1^2}}-\frac{1}{\sqrt{r^2-m_2^2}}<0$

${\Omega }_2=\frac{1}{m_3\sqrt{m_2^2-r^2{m}_3^2}}+\frac{1}{m_3\sqrt{m_1^2-r^2{m}_3^2}}>0$

说明，相交区域的面积会随着轨道高度的增大而减小，则卫星一个运行周期内观测到的 折射星数目也会随着轨道高度的增加而减小，其定位精度也会相应降低。

步骤二：星敏感器拍摄星图后，使用三角形算法识别星图中的正常星；

在本方法中，采用三角形星图识别方法，这种方法是星图识别中最常用的，在张光军《星 图识别》一书中有详细的介绍。此方法要求视场中的恒星数目不小于3颗，星敏感器同时还 用来观测折射星，因此为了有足够的非折射星来进行三角形识别应该选取较大视场的星敏感 器，取星敏感视场为12°×12°，像平面分辨率为1024×1024，依据式步骤一计算得到星敏感器 的最佳安装角度范围θ＝(64.77°,65.36°)，将安装角度定为中点θ＝65.06°。将拍摄星图中的任 意三颗星组成一个三角形，利用三角形的边角信息与存储在导航行库中信息进行匹配，若匹 配成功，说明这三颗星都没有发生折射，否则说明其中存在折射星，则重新选择三颗星进行 匹配，直到匹配成功为止。图3为使用三角形算法进行星图识别的结果。从图中我们发现， 有4颗星没有识别成功，初步判断为折射星。

步骤三：利用识别的正常星计算星敏感器光轴指向和卫星姿态；

利用匹配成功的正常星计算星敏感器光轴的指向和卫星的姿态，具体计算方法在房建成、 宁晓琳《天文导航原理及应用》中有详细的介绍。

步骤四：根据星敏感器光轴指向从星表中选星生成模拟折射星图；

利用星敏感器的光轴指向和星敏感器参数从星表中选择落在星敏感器视场中星光矢量为 s的折射星，将折射星投影到星敏感器像平面上；

$\alpha -{\theta }_R<\arccos (\frac{r}{\left|r\right|}·s)<\beta$

其中，ha＝20km，R_e为地球半径，s为折射星的星光矢量，r为 卫星相对于地心的位置矢量，|r|为r模长，代表卫星到地心的距离，θ_R为折射高度为20km 时对应的星光折射角。

步骤五：利用模拟折射星图识别折射星，根据识别结果计算星光折射角；

由图4可以看出折射星折射前、后在星敏感器像平面的位置。计算模拟星图中的折射星 和折射星投影之间的欧式距离，就可以识别出折射星。某颗星折射前后在像平面的位置坐标 为(P_x1,P_y1)，(P_x2,P_y2)，利用星敏感器的成像原理可以得到折射前、后星光在星敏感器坐标系 下的单位矢量S_c1，S_c2，如下式所述

$S_{\mathrm{ci}}=\frac{1}{\sqrt{P_{\mathrm{xi}}^2+P_{\mathrm{yi}}^2+f^2}}\left(\begin{array}{c}-P_{\mathrm{xi}}\\ -P_{\mathrm{yi}}\\ f\end{array}\right),(i=\mathrm{1,2})---\left(10\right)$

其中，f为星敏感器光学系统的角距

两个单位矢量的夹角即为当前时刻的星光折射角

$R=\arccos \left(S_{c1}^T{S}_{c2}\right)---\left(11\right)$

从图4中我们发现，有2颗折射星没有投影星，实际上这两颗星也是折射星，但是其星 光是经过对流层而被星敏感器观测到的，即这两颗星的星光折射高度要低于20km。而根据式 (1)得到的折射星的折射高度是在20km-50km的平流层，因此这两颗折射星没有投影星。当确 定了星敏感器的安装角度后，如图1所示，由于β-α＜θ_FOV，因此按照上述方案安装星敏感 器后，星敏感器会观测到经过对流层的星光，即被图1中对应的视场区域观测到的光线。 由于没有对应投影星，这种情况不影响折射星的识别。

步骤六：将星光折射角代入系统模型，星载计算机利用最优估计方法得到卫星的导航信 息；

系统模型包括系统的状态模型和两侧方程，系统的状态模型如下：

$\left(\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mu \frac{x}{r^3}[1-J_2\left(\frac{R_e}{r}\right)(7.5\frac{z^2}{r^2}-1.5)]+\Delta F_x\\ \frac{d{v}_y}{\mathrm{dt}}=-\mu \frac{y}{r^3}[1-J_2\left(\frac{R_e}{r}\right)(7.5\frac{z^2}{r^2}-1.5)]+\Delta F_y\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mu \frac{z}{r^3}[1-J_2\left(\frac{R_e}{r}\right)(7.5\frac{z^2}{r^2}-1.5)]+\Delta F_z\end{array}\right)---\left(12\right)$

$r=\sqrt{x^2+y^2+z^2}$

其中，系统状态为X＝[x y z v_x v_y v_z]^T，[x y z]^T代表卫星的位置矢量，[v_x v_y v_z]^T代表卫 星的速度矢量；μ＝3.986×10¹⁴m³/s²为地心引力常数；J₂＝0.00108263为地球引力系数； ΔF_x,ΔF_y,ΔF_z为地球非球形摄动的高阶摄动项和日、月摄动以及太阳光压摄动和大气摄动等 摄动力的影响。

量测方程为

$Z=h_a+v=\left(\begin{array}{c}\sqrt{r^2-u_1^2}+u_1\tan R_1-R_e+v_1\\ \sqrt{r^2-u_2^2}+u_2\tan R_2-R_e+v_2\\ ·\\ ·\\ ·\\ \sqrt{r^2-u_i^2}+u_i\tan R_i-R_e+v_i\end{array}\right)---\left(13\right)$

其中，i代表观测到的折射星数目，即量测方程的维数随着星敏感器所观测到的折射星数 目而变化。由于所建立的系统模型是非线性的，因此采用UKF滤波算法进行导航解算。由于 星敏感器在某些采样周期中无法观测到折射星，因此在没有折射星是滤波器只进行状态更新， 有折射星出现时同时进行状态更新和量测更新。

本发明以轨道高度为686km的对地观测卫星为例，仿真中卫星保持三轴稳定、对地定向 模式，UKF滤波周期为2s。仿真中的轨道数据由STK生成。仿真参数为：轨道长半轴7064.14km， 偏心率0，轨道倾角98.1358°，升交点赤经254.145°，近升角距0，星敏感器精度3″，视 场大小12°×12°，像平面分辨率为1024×1024，依据式步骤一计算得到星敏感器的最佳安装角 度范围θ＝(64.77°,65.36°)，将安装角度定为中点θ＝65.06°。

在上述仿真条件下，新方法仿真结果如图5和图6所示。图5为速度误差实验结果，图 6为位置误差实验结果。滤波收敛后系统的平均位置误差为103.9m(RMS)，最大位置误差为 259.6m，平均速度误差为0.168m/s，最大速度误差为0.342m/s。

将轨道高度进行变化，图7为位置误差随轨道高度变化的实验结果，图8为一个轨道周 期内观测到的折射星数目随轨道高度的变化结果，由图7和图8可以看出，随着轨道高度的 增加，星敏感器在一个轨道周期内观测到的折射星数目减少，同时导航系统的导航误差增大， 与式(9)得到的结论相同。为了验证星敏感器安装方案合理性，对安装角度与一个轨道周期中 观测到的折射星数目进行了仿真，仿真曲线如图9所示，仿真中将星敏感器的安装角度设置 为60°～70°。