一种快速脉冲星导航整周模糊度解算方法

摘要

一种快速脉冲星导航整周模糊度解算方法，分为地面实现部分与星上处理部分，地面实现部分对应于步骤1～步骤7，通在地面实现相关模型与模板的建立以节省星上计算时间，星上处理部分对应于步骤8～步骤12，称为基于粒子群优化的压缩模板匹配搜索方法。本发明满足了未来脉冲星导航系统快速并准确地建立航天器初始位置的需求，利用本发明可以用显著性水平参数来控制模糊度求解成功概率，并使用压缩模板指数来提高了模糊度解算效率，最终有效保证了模糊度能够成功与快速解算，能够正确并快速确定航天器的初始位置。

说明书

技术领域

本发明属于航天器自主导航技术领域，涉及一种快速脉冲星导航整周模糊度解算方法。 

背景技术

X射线脉冲星导航（简称脉冲星导航）是实现航天器长时间高精度自主导航最有希望取得突破的技术，具有重要的工程应用价值，备受国际航天机构关注。2004年，美国国防部提出X射线脉冲星导航研究计划（XNAV），目前已完成可行论证、关键技术攻关与地面验证，即将在国际空间站和高轨道卫星上开展空间飞行试验。此外，欧洲空间局（ESA）、俄罗斯、德国、日本、英国、印度和澳大利亚等国家或组织也启动了X射线脉冲星导航理论方法和试验验证研究。 

整周模糊度解算问题是脉冲星导航中的一个重要问题，因为航天器直接测量的是脉冲信号的小数相位，只有获得整数相位，才能得到航天器的绝对位置信息。脉冲信号整数相位是有待解算的未知量，称为整周模糊度（以下简称模糊度）。当导航系统发生故障而重启时，航天器先验位置信息丢失，需要解出正确的模糊度才能重建航天器的初始位置信息。 

由于相位量测是有误差的，模糊度解算问题是一个概率问题，需要通过搜索的方法来解算。S.I.Sheikh在其博士论文“The use of variable celestial X-ray sources for spacecraft navigation.Ph.D.Dissertation,Dept Aero Eng,Maryland Univ,College Park,MD,2005”中提出了搜索空间的概念，对候选模糊度进行逐个测式。J.Sara等在报告“Sala J,Urruela A,Villares X,et al.Feasibility study for a spacecraft navigation system relying on pulsar timing information.ARIADNA study03/4202,2004”中建立了基本搜索结构，即由3颗导航脉冲星与1颗检验脉冲星构成一个模糊度解算单元，不断地更换解算单元中的检验星进行搜索，直至得到3颗导航星的唯一模糊度解。国内学者通过相应研究，一定程度上提高了模糊度解算速度：谢振华等在文献“一种新的XPNAV系统解脉冲周期模糊算法.电子与信息学报,2008,30(9):2124–2127”中给出了一种线性形式的检验方程来缩短搜索时间；谢强等在文献“谢强,许录平,张华,等.脉冲星导航解周期模糊匹配搜索算法.系统工程与电子技术,2011,33(11):2498–2500”中开发了一种匹配搜索技术，使用匹配搜索模板来减小候选模糊度数量，进而加快搜索进程。但是，现有的模糊度解算方法的解算速度仍有待进一步提高，况且现有算 法检验阈值的选取比较经验化，不能有效控制解算成功概率。因此，本发明从实际工程应用角度出发，提出一种新型快速脉冲星导航整周模糊度解算方法，以满足未来脉冲星导航系统快速并准确地建立航天器初始位置的需求。 

发明内容

本发明的目的就在于：克服现有技术的不足，提供一种快速脉冲星导航整周模糊度解算方法，以满足未来脉冲星导航系统快速并准确地建立航天器初始位置的需求，利用本发明可以用显著性水平参数来控制模糊度求解成功概率，并使用压缩模板指数来提高了模糊度解算效率，最终有效保证了模糊度能够成功与快速解算，能够快速确定航天器的初始位置。 

本发明的技术解决方案：一种快速脉冲星导航整周模糊度解算方法，分为地面实现部分与星上处理部分，地面实现部分对应于步骤1～步骤7，通在地面实现相关模型与模板的建立以节省星上计算时间，星上处理部分对应于步骤8～步骤12，称为基于粒子群优化的压缩模板匹配搜索方法；步骤如下： 

步骤1:确定导航脉冲星与检验脉冲星。导航脉冲星为3颗，选取自转周期较大，相位量测精度较高的脉冲星，检验脉冲星为1颗，按照相位量测精度由高至低排序。3颗导航脉冲星与1颗检验脉冲星构成1个模糊度解算单元。 

步骤2:建立模糊度解算单元的模糊度量测方程。所建立模糊度量测方程表达式为：y=m+Bx+e 

式中：m=[m₁,m₂,m₃,m₄]^T，为模糊度解算单元各脉冲星的模糊度；y=[y₁,y₂,y₃,y₄]^T为模糊度解算单元各脉冲星的相位量测；e=[e₁,e₂,e₃,e₄]^T为模糊度解算单元各脉冲星的相位量测误差，设其服从零均值的4维正态分布，即e～N₄(0,R)，R为对称正定的方差矩阵，满足E(ee^T)=R>0；x=[x₁,x₂,x₃]^T为航天器3维位置矢量；B=[b₁,b₂,b₃,b₄]^T为脉冲星方向矢量矩阵，有 

b_i=c^-1f_i[cosδ_icosα_i,cosδ_isinα_i,sinδ_i]^T(i=1,2,3,4) 

其中，c为真空中的光速，f_i(i=1,2,3,4)为模糊度解算单元各脉冲星的自转频率，α_i与δ_i(i=1,2,3,4)分别为模糊度解算单元各脉冲星的赤径与赤纬。 

步骤3:建立线性模糊度接受域模型，具体步骤如下： 

步骤3.1:定义C≡I-BX，其中，X=(B^TR^-1B)^-1B^TR^-1，I为4阶单位阵。通过奇异值分解（SVD），将C分解为其中，σ₁为C的唯一非零奇异值，u₁与v₁分 别为C对应于σ₁的左奇异向量与右奇异向量。 

步骤3.2:求得方差矩阵R的逆的开方矩阵W，W满足W≥0，且W²=R^-1。 

步骤3.3:给定显著性水平α，求得自由度为1的χ²分布的上侧分位数(1)。 

步骤3.4:建立线性模糊度接受域模型，表达式为 

$|l^T\tilde{m}-\tilde{M}|\le D$

式中：为模糊度解算单元各脉冲星的模糊度估值，即为待检验模糊度；$l={\sigma }_1{v}_1^T,$记为l=[l₁,l₂,l₃,l₄]^T；${D=\left|\right|k\left|\right|}^{-1}\sqrt{{\chi }_{\alpha }^2\left(1\right)},$有k=Wu₁，表示向量的模值； $\tilde{M}=z^{T}k{\left|\right|k\left|\right|}^{-2},$z=WCy。 

步骤4:根据航天器任务确定航天器参考原点以及航天器至参考原点的最大可能距离R_max。以参考原点为球心，以R_max为半径的球应能覆盖航天器任何可能位置。 

步骤5:依据步骤4选定的R_max确定模糊度搜索空间。模糊度搜索空间用模糊度的下界和上界来描述，即对于一个模糊度解算单元，b_Li与b_Ui(i=1,2,3,4)分别为模糊度的下界与上界，即m_i∈[b_Li,b_Ui]。b_Li与b_Ui的计算公式为： 

$\left(\begin{array}{c}{b}_{\mathrm{Li}}=\mathrm{round}\left({-c}^{-1}{R}_{\max }{f}_i\right)\\ b_{\mathrm{Ui}}=\mathrm{round}\left(c^{-1}{R}_{\max }{f}_i\right)+1\end{array}\right)$

其中，round(·)表示四舍五入; 

步骤6:依据步骤5确定的模糊度搜索空间，再根据步骤4选定的R_max，确定压缩模板指数γ_m。γ_m的确定遵循以下三个原则：(1)γ_m>1；(2)要确保步骤8中粒子群算法在迭代次数N_I内以尽可能高的概率找到步骤8中定义的初始接受域的一个解；(3)如果可以满足第二条原则，γ_m值尽可能取大。基于上述三个原则，对于给定的模糊度搜索空间与选定的R_max，通过数学仿真对γ_m进行试配。 

步骤7:根据步骤6确定的γ_m值，建立压缩模板V。记3颗导航脉冲量模糊度增量估值为$\delta {\tilde{m}}_n={[\delta {\tilde{m}}_1,\delta {\tilde{m}}_2\delta {\tilde{m}}_3]}^T,$V是由$\delta {\tilde{m}}_n$构成的向量集合，满足 

$V=\left\{\delta {\tilde{m}}_n\right|\left|{\delta }_r\left(\delta {\tilde{m}}_n\right)\right|\le (1+1/{\gamma }_m)D/\left|l_4\right|,\left|\delta {\tilde{m}}_i\right|<b_{\mathrm{Ui}}{-b}_{\mathrm{Li}},\delta {\tilde{m}}_i\in Z,i=\mathrm{1,2,3}\}$

式中：Z代表整数集合；记l_n=[l₁,l₂,l₃]，有 

${\delta }_r\left(\delta {\tilde{m}}_n\right)=-1/l_4{l}_4^T\delta {\tilde{m}}_n-\mathrm{round}(-1/l_4{l}_n^T\delta {\tilde{m}}_n)\in [-\mathrm{0.5,0.5})$

步骤8:使用粒子群优化方法搜寻满足初始接受域的初始解。根据步骤3建立的线性 模糊度接受域模型与步骤6确定的γ_m值，确定初始接受域为 

$|l^T{\tilde{m}}_0-\tilde{M}|\le D/{\gamma }_m$

式中：为模糊度求解单元初始模糊度估值，有m_0i∈[b_Li,b_Ui](i=1,2,3,4)。对第j（初始设j=1）颗检验星，使用粒子群优化方法获得满足初始接受域的初始解进而确定3颗导航脉冲星初始模糊度估值

步骤9:根据步骤7建立的压缩模板V与步骤8获得的满足初始接受域的初始解按下式构建3颗导航脉冲星模糊度估值构成的备选集合S： 

$S=\{{\tilde{m}}_n={\tilde{m}}_{n0}+\delta {\tilde{m}}_n|\delta {\tilde{m}}_n\in V,{\tilde{m}}_{\mathrm{ni}}\in [b_{\mathrm{Li}},b_{\mathrm{Ui}}],i=\mathrm{1,2,3}\}$

步骤10:如果|S|=1，其中指集合元素的个数，则完成搜索，S中的元素即为导航星的模糊度；如果|S|=0，本次搜索失败，返回步骤8重新搜索；如果|S|>1，进行步骤11。 

步骤11:根据步骤9或步骤12确定的3颗导航脉冲星模糊度估值备选集合S。使j=j+1，对于第j颗检验星，使用3维顺序遍历来构建进一步缩小范围的模糊度估值备选集合S： 

其中，U由步骤2得到的线性形式的接受域模型得到，对于第j颗检验星有： 

$U=[(\tilde{M}-D{l}_n^T{\tilde{m}}_n)/l_4,(\tilde{M}+D-l_n^T{\tilde{m}}_n){/l}_4]$

步骤12:如果|S|=1，搜索完成，S中元素即为导航星的模糊度；如果|S|=0，本次搜索失败，返回步骤8重新搜索；如果|S|>1，返回步骤11。 

本发明与现有技术相比的优点在于： 

(1)本发明基于假设检验及统计学原理建立了描述模糊度求解问题的接受域模型，基于SVD分解的方法获得了线性形式的接受域方程，因而可以便于模糊度搜索方法的实现，并使得可以用显著性水平参数α来控制模糊度求解成功概率，最终有效保证了模糊度能够成功解算，能够正确确定航天器的初始位置。 

(2)本发明设计了基于粒子群优化的压缩模板匹配搜索方法进行模糊度搜索，引入参数γ_m来压缩匹配搜索模板的尺寸，提高了模糊度解算效率，在航天器初始位置不确定程度较大情形下，本发明所述方法的计算效率约是当前先进解算方法（即文献“谢强,许录平,张华,等.脉冲星导航解周期模糊匹配搜索算法.系统工程与电子技术,2011,33(11):2498–2500”中所述方法）的两倍，最终有效保证了模糊度能够快速解算，能够快速确定航 天器的初始位置。 

附图说明

图1为本发明快速脉冲星导航整周模糊度解算方法流程图。 

具体实施方式

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

如图1所示，本发明实现快速脉冲星导航整周模糊度解算的步骤如下： 

步骤1:确定导航脉冲星与检验脉冲星。导航脉冲星为3颗，选取自转周期较大，相位量测精度较高的脉冲星，检验脉冲星为1颗，按照相位量测精度由高至低排序。相位量测精度可按文献“Sheikh S I.The use of variable celestial X-ray sources for spacecraft navigation.Ph.D.Dissertation,Dept Aero Eng,Maryland Univ,College Park,MD,2005”所述方法进行估算。3颗导航脉冲星与1颗检验脉冲星构成1个模糊度解算单元。本实施例选取3颗导航星为PSR J1846-0258，PSR B1509-58与PSR J1930+1852，选取第1颗检验星为PSR B1823-13，其他检验星依次为PSR J1124-5916，PSR B1706-44，PSR B0833-45，PSR J1617-5055，PSR J1420-6048，PSR J0205+6449。 

步骤2:建立模糊度解算单元的模糊度量测方程。所建立模糊度量测方程表达式为： 

y=m+Bx+e 

式中：m=[m₁,m₂,m₃,m₄]^T，为模糊度解算单元各脉冲星的模糊度；y=[y₁,y₂,y₃,y₄]^T为模糊度解算单元各脉冲星的相位量测；e=[e₁,e₂,e₃,e₄]^T为模糊度解算单元各脉冲星的相位量测误差，设其服从零均值的4维正态分布，即e～N₄(0,R)，R为对称正定的方差矩阵，满足E(ee^T)=R>0；x=[x₁,x₂,x₃]^T为航天器3维位置矢量，本实施例设定x为x=[2.6×10⁷,-2.6×10⁷,2.6×10⁷]^Tm；B=[b₁,b₂,b₃,b₄]^T为脉冲星方向矢量矩阵，有 

b_i=c^-1f_i[cosδ_icosα_i,cosδ_isinα_i,sinδ_i]^T(i=1,2,3,4) 

其中，c为真空中的光速，fi(i=1,2,3,4)为模糊度解算单元各脉冲星的自转频率，α_i与δ_i(i=1,2,3,4)分别为模糊度解算单元各脉冲星的赤径与赤纬。 

步骤3:建立线性模糊度接受域模型，具体步骤如下： 

步骤3.1:定义C≡I-BX，其中，X=(B^TR^-1B)^-1B^TR^-1，I为4阶单位阵。通过奇异值分解（SVD），将C分解为其中，σ₁为C的唯一非零奇异值，u₁与v₁分别为C对应于σ₁的左奇异向量与右奇异向量。 

步骤3.2:求得方差矩阵R的逆的开方矩阵W，W满足W≥0，且W²=R^-1。 

步骤3.3:给定显著性水平α，求得自由度为1的χ²分布的上侧分位数(1)。本实施例α设为0.05，α可以控制模糊度解算成功概率，若使用N颗检验星成功解算出模糊度，则模糊度解算成功概率大于等于(1-α)^N。 

步骤3.4:建立线性模糊度接受域模型，表达式为 

$|l^T\tilde{m}-\tilde{M}|\le D$

式中：为模糊度解算单元各脉冲星的模糊度估值，即为待检验模糊度；$l={\sigma }_1{v}_1^T,$记为l=[l₁,l₂,l₃,l₄]^T；${D=\left|\right|k\left|\right|}^{-1}\sqrt{{\chi }_{\alpha }^2\left(1\right)},$有k=Wu₁，表示向量的模值； $\tilde{M}=z^{T}k{\left|\right|k\left|\right|}^{-2},$z=WCy。 

步骤4:根据航天器任务确定航天器参考原点以及航天器至参考原点的最大可能距离R_max。以参考原点为球心，以R_max为半径的球应能覆盖航天器任何可能位置。本实施例描述对象选为探月航天器，故参考原点设为地心，R_max根据地月距离加上余量设为5×10⁸m。 

步骤5:依据步骤4选定的R_max确定模糊度搜索空间。模糊度搜索空间用模糊度的下界和上界来描述，即对于一个模糊度解算单元，b_Li与b_Ui(i=1,2,3,4)分别为模糊度的下界与上界，即m_i∈[b_Li,b_Ui]。b_Li与b_Ui的计算公式为： 

$\left(\begin{array}{c}{b}_{\mathrm{Li}}=\mathrm{round}\left({-c}^{-1}{R}_{\max }{f}_i\right)\\ b_{\mathrm{Ui}}=\mathrm{round}\left(c^{-1}{R}_{\max }{f}_i\right)+1\end{array}\right)$

其中，round(·)表示四舍五入。 

步骤6:依据步骤5确定的模糊度搜索空间，再根据步骤4选定的R_max，确定压缩模板指数γ_m。γ_m的确定遵循以下三个原则：(1)γ_m>1；(2)要确保步骤8中粒子群算法在迭代次数N_I内以尽可能高的概率找到步骤8中定义的初始接受域的一个解；(3)如果可以满足第二条原则，γ_m值尽可能取大。基于上述三个原则，对于给定的模糊度搜索空间与选定的R_max，通过数学仿真对γ_m进行试配。对于本实施例，通过试配，本发明γ_m设为8。 

步骤7:根据步骤5选择的γ_m值，建立压缩模板V。记3颗导航脉冲量模糊度增量估值为$\delta {\tilde{m}}_n={[\delta {\tilde{m}}_1,\delta {\tilde{m}}_2\delta {\tilde{m}}_3]}^T,$V是由$\delta {\tilde{m}}_n$构成的向量集合，满足 

$V=\left\{\delta {\tilde{m}}_n\right|\left|{\delta }_r\left(\delta {\tilde{m}}_n\right)\right|\le (1+1/{\gamma }_m)D/\left|l_4\right|,\left|\delta {\tilde{m}}_i\right|<b_{\mathrm{Ui}}{-b}_{\mathrm{Li}},\delta {\tilde{m}}_i\in Z,i=\mathrm{1,2,3}\}$

式中：Z代表整数集合；记l_n=[l₁,l₂,l₃]，有 

${\delta }_r\left(\delta {\tilde{m}}_n\right)=-1/l_4{l}_4^T\delta {\tilde{m}}_n-\mathrm{round}(-1/l_4{l}_n^T\delta {\tilde{m}}_n)\in [-\mathrm{0.5,0.5})$

步骤8:使用粒子群优化方法搜寻满足初始接受域的初始解。根据步骤2建立的线性模糊度接受域模型与步骤6确定的γ_m值，确定初始接受域为 

$|l^T{\tilde{m}}_0-\tilde{M}|\le D/{\gamma }_m$

式中：为模糊度求解单元初始模糊度估值，有m_0i∈[b_Li,b_Ui](i=1,2,3,4)。对第j（初始设j=1）颗检验星，使用粒子群优化方法获得满足初始接受域的初始解进而确定3颗导航脉冲星初始模糊度估值粒子群优化方法的具体步骤如下。 

步骤8.1:在模糊度搜索空间内任意设定S个粒子的初始位置。粒子位置即为模糊度求解单元的模糊度估值(k=1,2,…,S)，满足m_ki∈[b_Li,b_Ui](i=1,2,3,4)。本实施例取S=40。 

步骤8.2:对于每个粒子（k=1,2,…,S），执行步骤8.2.1～步骤8.2.4： 

步骤8.2.1:对于i=1,2,3，生成随机数r_ki～U(b_Li,b_Ui)（U代表均匀分布），并使得 ${\tilde{m}}_{\mathrm{ki}}=\mathrm{round}\left(r_{\mathrm{ki}}\right).$

步骤8.2.2:使得$m_{k4}=\mathrm{round}\left[(\tilde{M}-l_n^T{\tilde{m}}_{\mathrm{nk}}){/l}_4\right],$式中${\tilde{m}}_{\mathrm{nk}}={[{\tilde{m}}_{1k},{\tilde{m}}_{2k},{\tilde{m}}_{3k}]}^T.$

步骤8.2.3:初始化粒子个体最佳位置为其初始位置：其中，p_k=[p_k1,p_k2,p_k3,p_k4]^T表示第k个粒子的个体最佳位置。 

步骤8.2.4:如果f(p_k)<f(g)，则更新g：g←p_k，其中，f为目标函数，定义为 $f\left(\tilde{m}\right)=|l^T\tilde{m}-\tilde{M}|,$g=[g₁,g₂,g₃,g₄]^T指群体最佳位置。 

步骤8.3:定义最大迭代次数为N_I，本实施例中取N_I=10，直至f(g)<D/γ_m或迭代次数达到N_I，对于每个粒子k=1,2,…,S，重复执行步骤8.3.1～步骤8.3.3： 

步骤8.3.1:对于i=1,2,3,4，执行步骤8.3.1.1～步骤8.3.1.3： 

步骤8.3.1.1:生成随机数r_p～U(0,1)，r_g～U(0,1)，及r_d～U(0,1)。 

步骤8.3.1.2:计算粒子速度$v_{\mathrm{ki}}=c_p{r}_p(p_{\mathrm{ki}}-{\tilde{m}}_{\mathrm{ki}})+c_g{r}_g(g_i-{\tilde{m}}_{\mathrm{ki}}),$其中，c_p=c_g=2 

步骤8.3.1.3:更新粒子位置： 

${\bar{m}}_{\mathrm{ki}}\leftarrow \left(\begin{array}{cc}{\bar{m}}_{\mathrm{ki}}+1,& v_{\mathrm{ki}}>0\\ {\bar{m}}_{\mathrm{ki}}-1,& v_{\mathrm{ki}}<0,\\ {\bar{m}}_{\mathrm{ki}}-d,& v_{\mathrm{ki}}=0\end{array}\right),d=\left(\begin{array}{cc}1,& r_d\le 1/3\\ 0,& 1/3<r_d\le 2/3\\ -1,& r_d>2/3\end{array}\right)$

步骤8.3.2:如果$f\left({\tilde{m}}_k\right)<f\left(p_k\right),$更新p_k：$p_k\leftarrow {\tilde{m}}_k.$

步骤8.3.3:如果f(p_k)<f(g)，更新g：g←p_k。 

步骤8.4:使得进而确定3颗导航脉冲星初始模糊度估值 ${\tilde{m}}_{n0}={[{\tilde{m}}_{01},{\tilde{m}}_{02},{\tilde{m}}_{03}]}^T.$

步骤9:根据步骤7建立的压缩模板V与步骤8获得的满足初始接受域的初始解按下式构建3颗导航脉冲星模糊度估值构成的备选集合S： 

$S=\{{\tilde{m}}_n={\tilde{m}}_{n0}+\delta {\tilde{m}}_n|\delta {\tilde{m}}_n\in V,{\tilde{m}}_{\mathrm{ni}}\in [b_{\mathrm{Li}},b_{\mathrm{Ui}}],i=\mathrm{1,2,3}\}$

步骤10:如果|S|=1，其中||指集合元素的个数，则完成搜索，S中的元素即为导航星的模糊度；如果|S|=0，本次搜索失败，返回步骤8重新搜索；如果|S|>1，进行步骤11。 

步骤11:根据步骤9或步骤12确定的3颗导航脉冲星模糊度估值备选集合S。使j=j+1，对于第j颗检验星，使用3维顺序遍历来构建进一步缩小范围的模糊度估值备选集合S： 

其中，U由步骤2得到的线性形式的接受域模型得到，对于第j颗检验星有： 

$U=[(\tilde{M}-D{l}_n^T{\tilde{m}}_n)/l_4,(\tilde{M}+D-l_n^T{\tilde{m}}_n){/l}_4]$

步骤12:如果|S|=1，搜索完成，S中元素即为导航星的模糊度；如果|S|=0，本次搜索失败，返回步骤8重新搜索；如果|S|>1，返回步骤11。 

通过500次的蒙特卡洛试验，本实施例的平均检验星个数N=2.354，预测成功概率(1-α)^N=0.886，实际成功概率为0.930；本实施例计算时间为文献“一种新的XPNAV系统解脉冲周期模糊算法.电子与信息学报,2008,30(9):2124–2127”中方法的39.3%，为文献“谢强,许录平,张华,等.脉冲星导航解周期模糊匹配搜索算法.系统工程与电子技术,2011,33(11):2498–2500”中方法的54.3%。可见本发明可以有效控制模糊度解算成功概率并有效提高解算速度。 

以上实施例为本发明的较佳实施方式之一，凡是在本发明的精神和原则之下进行的等同替换，局部改进都将视为在本发明的保护范围之内。