弹道导弹作为一种高技术的精确打击武器，具有飞行时间短、速度快、精度要求高等特点，将星光制导应用于弹道导弹有利于提高其制导精度和自主导航能力[1-2]。星敏感器是星光制导的关键设备，但星敏感器的低数据输出率，成为制约星敏感器在弹道导弹上应用的关键问题之一[3-5]。提高星图识别效率有利于提高星敏感器的数据输出率，而减小导航星表的存储容量对提高星图识别效率具有重要作用，因此开展弹载导航星表研究对星敏感器在弹道导弹上的应用具有重要意义[6-7]。
由于弹道导弹的特殊应用环境，星敏感器视场只需扫描天球的局部区域，所以针对扫描的局部区域，构建符合弹道导弹实际应用所需的弹载导航星表，能减少导航星表存储量，有效提高星图识别效率[8-11]。目前，适用于弹载导航星表构建的方法有椭圆投影法[12]、“伪天球”法[13-14]等。椭圆投影法根据导弹飞行轨迹在地心惯性坐标系平面上的椭圆投影，在投影区域内选取导航星。该方法所确定的局部选星区域远大于视场扫描的实际区域，选取的导航星数量较多，且需要处理视场扫描区域坐标越界的问题，计算较为繁琐。“伪天球”法根据导弹飞行参数构建“伪天球”坐标系，在“伪天球”坐标系中确定视场扫描的局部区域，采用基于近似均匀分布的基准点及“距离—星等”加权法，在全天域范围内均匀选出候选星集合，并将所有候选星投影至“伪天球”坐标系，在视场扫描区域内选取导航星。该方法避免了处理扫描区域坐标越界的问题，但要将大量候选星投影至“伪天球”坐标系，计算量较大。
本文针对椭圆投影法和“伪天球”法的缺点，以弹载导航星表构建的完备性和快速性为原则，提出一种基于空间平面选取弹载导航星的方法，有效减少了弹载导航星表的构建时间。
1 空间平面的构建
为了确定星敏感器视场扫描区域，需构建四个空间平面。如图1所示，OxCyCzC为天球坐标系，G为关机点，G′点是关机点G在天球上的投影，R为再入点，R′点是再入点R在天球上的投影。假设星敏感器在导弹关机点开始工作，在再入点停止工作。星敏感器视场扫描区域的起始边界中点应在G′点附近，定义T1、T2为起始边界点，G″为起始边界T1T2的中点。星敏感器视场扫描区域的结束边界中点应在R′点附近，定义T3、T4为结束边界点，R″为结束边界T3T4的中点。为了确定空间平面，通过点G′与点G″的夹角关系，点R′与点R″的夹角关系，解算出点G″及点R″坐标；通过点T1、T2与点G″的夹角关系，点T3、0T4与点R″的夹角关系，解算出点T1、T2、T3、T4坐标。根据边界4点T1、T2、T3、T4及相应法线向量，可沿视场扫描区域边界构建四个空间平面，确定视场扫描区域。
图1 天球坐标系下星敏感器视场扫描区域
1.1 投影点G′、R′的确定
根据标准弹道，可预先计算出关机点G在地球坐标系下的地理位置(αgj,δgj)，为便于坐标转换计算，根据式(1)，将关机点G由经纬度表示形式转换成三维坐标表示形式(xgj,ygj,zgj)：
(1)
其中，αgj、δgj为地球坐标系下地理经、纬度。
为便于在天球坐标系下开展研究，通过地球坐标系与天球坐标系的转换关系，可解算出关机点G在天球坐标系中的投影点
(2)
其中:
为地球坐标系到天球坐标系的转换矩阵；Sgp为发射瞬时春分点时角。Sgp可根据式(3)计算得到：
Sgp=Sgp0+T-8h
(3)
其中，Sgp0为当天零时0h的春分点时角(Sgp0=6h40m+d×3m56s，d为从元旦起至当天的天数，元旦子夜时的春分点时角为6h40m)； T为发射瞬时北京时间。
同理，根据再入点R的地理位置(αzr,δzr)，利用式(1)、式(2)、式(3)的方法，可解算出再入点R的天球投影点R′的三维坐标
1.2 边界中点G″、R″的确定
星敏感器视场扫描区域的起始边界中点G″与投影点G′的夹角为
(4)
其中，φcx(tg)为关机时刻tg时的飞行程序角；Δφ为飞行程序角偏差；θFOV为星敏感器视场大小。
由图1可知：点G″可看作是点G′绕G′OR′平面单位法线向量
正向旋转-φG得到的，即
(5)
其中，(G′)T为G′的转置矩阵；
为绕单位法线向量正向旋转-φG的坐标旋转矩阵。
星敏感器视场扫描区域的结束边界中点R″与投影点R′的夹角为
(6)
其中，Θ为再入角；θR为攻角；Δ为Θ和θR的总体误差。
同理，由图1可知，点R″可看作是点R′绕正向旋转φR得到的，即
(7)
1.3 边界点T1、T2、T3、T4的确定
由图1可知，点T1、T2与G″的夹角及点T3、T4与R″的夹角同为边界点与边界中点的夹角，令其为Ω，根据标准弹道可得到：
(8)
其中，ψ为偏航角； γ为滚动角。
根据式(4)，则
(9)
其中，
为绕T1OT2平面的单位法线向量
正向旋转-Ω的坐标旋转矩阵；
为绕T1OT2平面的单位法线向量正向旋转Ω的坐标旋转矩阵；
为绕T4OT3平面的单位法线向量
正向旋转-Ω的坐标旋转矩阵；
为绕T4OT3平面的单位法线向量正向旋转Ω的坐标旋转矩阵。
1.4 空间平面的确定
为了便于在星敏感器视场的扫描区域T1T2T3T4内选取导航星，需构建T1OT2、T2O2T4、T4OT3、T1O1T3四个平面。根据已知坐标点和法线向量，并以便于选星为原则，令平面法线向量方向均指向扫描区域T1T2T3T4一侧，构建4个空间平面方程：
(10)
其中，x、y、z为恒星的三维坐标；xT1、yT1、zT1为点T1的三维坐标；xT4、yT4、zT4为点T4的三维坐标；xn、yn、zn为T2O2T4、T3O1T1平面法线向量的三维坐标；xg、yg、zg为T1OT2平面法线向量的三维坐标；xr、yr、zr为T4OT3平面法线向量的三维坐标。
在天球坐标系中，所构建的4个空间平面分布如图2所示。
图2 空间平面分布
其中，阴影部分是星敏感器视场扫描区域T1T2T3T4。由图2可知：采用空间平面，能准确确定视场扫描区域T1T2T3T4。
2 导航星的选取
在星敏感器视场扫描区域T1T2T3T4内选取导航星，可构建满足导弹发射任务需要的弹载导航星表。利用上述构建的4个空间平面方程，能保证构建弹载导航星表的快速性，具体步骤如下：
1) 在候选星中，选取位于T1O1T3、T2O2T4两空间平面间的恒星，并将选取出的恒星设为星表H，如图3所示，星表H中恒星在天球坐标系中的分布呈环状。选取公式为
(11)
其中，x、y、z为恒星的三维坐标。
图3 星表H分布
2) 在星表H中，选取位于T1OT2平面法线向量指向一侧的恒星，并将选取出的恒星设为星表H12，如图4所示，星表H12中恒星在天球坐标系中的分布呈弧状。选取公式为
xg·x+yg·y+zg·z≥0
(12)
图4 星表H12分布
3) 在星表H中，选取位于T4OT3平面法线向量
指向一侧的恒星，并将选取出的恒星设为星表H43.，如图5所示，星表H43.中恒星在天球坐标系中的分布呈弧状。选取公式为
-xr·x-yr·y-zr·z≥0
(13)
图5 星表H43分布
4) 讨论∠G″OR″的大小。∠G″OR″是视场扫描的起始边界与结束边界的夹角，即T1OT2平面与T4OT3平面的夹角。由图2、图4、图5可知，如果∠G″OR″≤180°，位于扫描区域T1T2T3T4内的导航星即属于星表H12且属于星表H43；如果∠G″OR″>180°，位于扫描区域T1T2T3T4内的导航星属于星表H12或属于星表H43.。因此，弹载导航星表S为
(14)
3 仿真验证与分析
为了验证本文提出的基于空间平面构建弹载导航星表的完备性和快速性，设计仿真实验，并与椭圆投影法、“伪天球”法进行比较[15]。仿真条件如下：由于在全天域范围内均匀化选取候选星集合的方法多样且各有利弊，为简化实验步骤，增加实验可靠性，本文采用SKYMAP2000星表，从中提取角距大于0.5°，星等小于等于6.0的星，并剔除双星和变星，共计 4 908颗作为候选星。星敏感器采用圆形视场， θFOV=12°。由弹道发生器生成一条标准弹道，计算出关机点投影点G′(90.22°,20.98°)，再入点投影点R′(123.16°,16.26°)，采用Monte Carlo方法，在星敏感器的扫描区域内随机生成 10 000个视场，仿真结果如下。
图6～图8分别是采用椭圆投影法、“伪天球”法和本文方法对该标准弹道选取弹载导航星的分布图，图中T1、T2、T3、T4是视场扫描区域的边界点。对图6、图7、图8中的弹载导航星数进行统计，统计结果如表1。
图6 椭圆投影法选取弹载导航星分布
图7 “伪天球”法选取弹载导航星分布
图8 本文方法选取弹载导航星分布
表1 三种方法的仿真实验结果
方法出现3颗以上导航星概率/%导航星选取数/颗椭圆投影法100751“伪天球”法100410本文方法100410
对表1中数据进行分析可知：3种方法识别导航星的成功率均为100%，表明本文方法满足文献[12-14]对弹载导航星表完备性的要求。由图6～8可直观看出，椭圆投影法所确定的选星范围明显超出星敏感器的视场扫描区域，“伪天球”法与本文方法所确定的选星范围正好在视场扫描区域的边界之内。椭圆投影法选取的弹载导航星为751颗，比“伪天球”法和本文方法要多出341颗，多出的弹载导航星正是由于椭圆投影法确定的选星范围超出星敏感器视场区域造成的。
根据上述实验结果，可知随着弹载导航星的数量增加，计算量也随之增大，从而影响弹载导航星表的构建速度。此外，“伪天球”法需进行坐标转换，将全天域内选出的侯选导航星投影至“伪天球”坐标系，在视场扫描区域内选取导航星，计算量仍然较大。本文方法利用所构建的四个空间平面，直接在天球坐标系下对候选星进行筛选，相比于椭圆投影法、“伪天球”法，本文方法避免了扫描区域坐标越界的问题，且无需对候选星进行坐标转换，能有效减少计算量，大大缩短弹载导航星表构建时间。
为进一步验证本文方法的快速性，由弹道发生器随机生成100条标准弹道，分别采用椭圆投影法、“伪天球”法和本文方法对100条标准弹道进行弹载导航星表构建，统计不同方法构建弹载导航星表的算法耗时，仿真结果如图9所示，数据统计结果见表2。
图9 三种方法构建弹载导航星表耗时对比
表2 三种方法的平均法耗时
方法平均耗时/ms标准差椭圆投影法46.050.44“伪天球”法31.240.61本文方法4.270.16
由图9和表2数据可知，在对100条标准弹道进行弹载导航星表构建时，采用椭圆投影法的选星平均耗时最长，其次为“伪天球”法，本文方法耗时最短，平均耗时为4.27 ms。相比于椭圆投影法，本文方法选取弹载导航星的速度提高了90.7%，相比于“伪天球”法提高了86.3%，表明本文方法构建弹载导航星表的速度更快。同时，从表2三种方法选取弹载导航星的耗时标准差可以看出，本文方法的耗时标准差最低，说明本文方法构建弹载导航星表的速度快且稳定。
4 结论
本文以星敏感器在弹道导弹上的应用为背景，通过构建空间平面对弹载导航星进行选取，避免对候选星进行坐标转换计算以及在扫描区域内选星时考虑扫描区域坐标越界，保证了导航星表的完备性，大大缩短了弹载导航星表构建时间，具有工程实际应用价值。
[1] 李欣璐,杨进华,张刘,等.空间立体角“准均匀分布”导航星表划分[J].光学精密工程,2014,22(8):2242-2246.
[2] 陆敬辉,王宏力,孙渊,等.应用相近模式向量的主星星图识别算法[J].红外与激光工程,2011,40(1):164-167.
[3] 房建成,宁晓琳.航天器自主导航原理与方法[M].北京:国防工业出版社,2006.
[4] 朱长征,居永忠,杜晓辉.导航星库制定方法研究[J].宇航学报,2016,31(5):1327-1330.
[5] MORTARI D,JUNKINS J L,SAMAAN M A.Lost-in-space Pyramid Algorithm for Robust Star Pattern Recognition[C]//The 24th Annual AAs Guidance and Control Conference,Breckenridge,Colorado,2001.
[6] 张磊,何昕,魏仲慧,等.三角形星图识别算法的改进[J].光学精密工程,2010,18(2):459-463.
[7] DAVID D N,JAMES P A.Fast Access and Low Memory Star Pair Catalog for Star Pattern Identification[J].Journal of Guidance,Control,and Dynamics,2010,33(5):1396-1403.
[8] 蒋明,喻明艳,王进祥,等.一种新的全天自主星图识别算法[J].宇航学报,2017,28(4):1020-1024.
[9] BENJAMIN B S,DANIELE M.A Survey on Star Identification Algorithms[J].Algorithms,2009(2):93-107.
[10] 徐帆.惯性/天文/卫星组合导航方法研究[D].北京:北京航空航天大学,2008.
[11] 王常虹,吴志华,李葆华.一种大角加速度下的星跟踪算法[J].中国惯性技术学报,2010,18(6):706-710.
[12] 陆敬辉,王宏力,文涛,等.弹载导航星表的设计与实现[J].宇航学报,2011,32(9):1895-1901.
[13] 王宏力,陈聪,崔祥祥,等.适用于弹道导弹的导航星选取方法[J].北京航空航天大学学报,2013,39(2):143-147.
[14] 陈聪.弹载导航星表构建方法研究[D].西安:第二炮兵工程大学,2012.
[15] 王宏力,崔祥祥,陆敬辉.基于惯导姿态信息的高鲁棒性星图识别算法[J].中国惯性技术学报,2010,18(6):729-732.