CN104977594A - 基于非合作导航卫星信号的无码多普勒定轨方法 - Google Patents

Abstract

本发明公开了基于非合作导航卫星信号的无码多普勒定轨方法，通过接收机的射频前端将导航卫星信号转换为数字中频信号，对中频信号进行消码处理后得到连续载波信号，并估计载波信号的初始频率；将该频率值作为锁相环路跟踪的初值，恢复载波并测量提取出各个导航星信号的多普勒频移；然后采用基于多普勒频移特征的匹配算法寻找出某频移对应的导航星；建立高精度测量模型，结合轨道动力学模型，采用卡尔曼滤波器估计出高精度轨道参数。本发明能在明码信号不可用的情况下，利用卫星导航系统发送的加密导航信号，实现高精度定轨。

Description

基于非合作导航卫星信号的无码多普勒定轨方法

技术领域

本发明属于卫星导航技术领域，特别涉及了基于非合作导航卫星信号的无码多普勒定轨方法。

背景技术

卫星轨道精密测量是卫星应用的基础。目前，国内外采用的定轨方式主要是依赖地基测控网，结合星上应答机获得测距、测速信息后精密校正轨道并上注，但是这种定轨方式不但在恶劣环境时生存保障较为困难，而且定轨精度约为百米到十米量级并随着星上轨道预报后精度逐渐下降，有时难以满足精度要求。

法国提出并研制的DORIS定轨系统，是一套星上自主实时定轨软件，其解算周期为10s，能够快速、准确地确定卫星轨道。中国首颗海洋动力环境卫星海洋二号(HY-2)主要用于观测全球海洋动力环境参数。

卫星激光测距系统观测精度可达到毫米级，可以满足精度要求，但是全球观测台站少、分布不均，以及受观测区域和气候条件的限制，部分台站还无法实现常规白天测距，使观测数据较为有限，而且SLR设备较为复杂庞大，价格不菲。

目前，GPS、GLONASS等国外的卫星导航系统是当前天基跟踪测量的主流手段，是地基定轨系统的重要补充，其中以GPS为主。但是GPS卫星导航系统是国外军方控制，在恶劣环境下，明码信号面临关闭的风险，这样就无法通过普通的接收机得到卫星伪距、伪距率，也就无法实现定位和定轨。我国的北斗系统目前是一个区域导航系统，尚不能实现全轨道覆盖。

发明内容

为了解决上述背景技术提出的技术问题，本发明旨在提供基于非合作导航卫星信号的无码多普勒定轨方法，在明码信号不可用的情况下，利用卫星导航系统发送的加密导航信号，实现高精度定轨。

为了实现上述技术目的，本发明的技术方案为：

基于非合作导航卫星信号的无码多普勒定轨方法，包括以下步骤：

(1)对接收到的导航卫星信号经数字下变频后，采样得到数字中频信号，消除数字中频信号的保密调制码，根据消码后的信号计算得到载波频率的估计值；

(2)将步骤(1)得到的载波频率的估计值作为载波跟踪环跟踪的初始频率，产生载波跟踪环的本振信号，将步骤(1)中消码后的信号作为载波跟踪环的输入信号，恢复出载波信号，并测量提取出各个导航卫星信号的多普勒频移；

(3)计算各个导航卫星信号的先验值，结合步骤(2)得到的多普勒频移的实测值，建立含有接收机频差的多普勒频移实测值和先验值的识别匹配算法，完成导航星的识别与匹配；

(4)利用统计定轨算法实现轨道测量，利用轨道动力学方程建立状态方程，利用载波相位测量值建立相位差测量方程，再利用卡尔曼滤波器估计用户卫星的轨道参数。

其中，步骤(1)中所述根据消码后的信号得到载波频率的估计值的具体步骤；

(A)选取合适的搜索步长，产生若干候选频率，对每一个候选频率产生一个载波相位调整信号；

(B)将步骤(A)产生的载波相位调整信号分别与消码后的信号相乘，对每一个相乘后的信号按照一定单位长度进行分块，累加分块数据，从而得到信号增强的新数据；

(C)对步骤(B)得到的新数据进行FFT运算，检查频谱结果，寻找超过门限的尖峰，对于每一个尖峰，记录其对应的频率值，即为载波频率的估计值。

其中，步骤(3)的具体步骤：

(a)根据步骤(2)恢复出的载波信号，由用户卫星的初始轨道和周期上注的导航卫星星历计算出各导航卫星信号的多普勒频移先验值；

(b)将步骤(2)得到的多普勒频移实测值和步骤(a)得到的多普勒频移先验值分别按照高低顺序排列成两组数据样本，引入待求解的整体平移量，该整体平移量用于将多普勒频移先验值作整体平移；

(c)对各导航卫星信号的多普勒频移实测值与平移后的多普勒频移先验值的频差绝对值求和得到匹配指标；

(d)求解在匹配指标取最小值时的整体平移量，将多普勒频移先验值作整体平移，实现多普勒频移实测值与多普勒频移先验值的匹配，实现导航星识别。

其中，在步骤(d)中，采用三分法求解在匹配指标取最小值时的整体平移量，其具体步骤：

(Ⅰ)假设某导航卫星信号的多普勒频移先验值作平移的频率区间是选取该区间的1/3、2/3处的坐标 把整个区间分成三个小区间[m₁,m₂]，

(Ⅱ)分别计算在m₁、m₂处匹配指标e(m₁)、e(m₂)，若e(m₁)＞e(m₂)，则舍弃区间令否则舍弃令

(Ⅲ)重复步骤(Ⅰ)-(Ⅱ)，直到小于某一预设阈值时，即为整体平移量。

其中，在步骤(1)中，采用平方法消除数字中频信号的保密调制码。

其中，在步骤(1)中，采用交叉相关法消除数字中频信号的保密调制码。

采用上述技术方案带来的有益效果：

本发明可以提高利用卫星导航系统实现卫星定轨的可靠性和安全性，可用作卫星定轨接收机的备份工作模式。同时本发明可望达到较高定轨精度。

附图说明

图1是本发明的流程图；

图2是射频信号下变频原理框图；

图3是平方法消码示意图；

图4是交叉相关法消码示意图；

图5是锁相环结构示意图；

图6是本发明三分法的3种情况示意图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

如图1所示本发明的流程图，用户卫星射频前端将导航卫星信号转换为数字中频信号，对得到的中频信号进行消码处理得到连续载波信号；基于相位补偿、数据块累加及频谱分析来精确估计消码后的载波的频率；将频率估计值作为载波跟踪环路的频率初始值，载波环路对载波信号动态跟踪，测量出各个导航星信号的多普勒频移；然后采用基于多普勒频移特征的匹配算法寻找出某频移对应的导航星；将导航卫星当作一个频率已知的无线电信标，建立高精度测量模型，结合轨道动力学模型，采用扩展卡尔曼滤波算法估计出高精度轨道参数。其具体过程如下。

1、中频信号获取及其消码处理

在GPS接收机中，通常采用混频方式将射频信号从高频下变频得到频率较低的中频模拟信号，继而通过模数转换器转换为数字中频信号。射频信号下变频原理如图2所示。

在图2中，射频信号r₁(t)和本地载波信号r₂(t)经过混频器混频，得到下变频的中频信号。便于信号处理，射频信号采用正交下变频得到正交的中频信号。混频后的信号分为两路，分别采样得到正交的卫星中频信号r(t)。导航卫星信号经正交采样后可表示为：

r(t)＝ACDe^j2πft+n(t)  (1)式中，A为信号幅度，C为未知的伪随机码，D为电文数据位，e^j2πft是用复指数形式描述的正弦信号，其中f为中频载波频率，t是时间，π是圆周率，j表示复数的虚部，n(t)为白噪声。

对中频信号r(t)进行平方或者交叉相关法进行消码处理，平方法和交叉相关法消码原理分别如图3、4所示。

以平方法为例，对中频信号r(t)平方可得：

r²(t)＝A²C²D²e^j2π2ft+n²(t)+2ACDe^j2πft×n(t)

                                   (2)

＝A²e^j2π2ft+n′(t)

处于±1状态的调制码，经平方后均为+1，而+1不改变信号相位和载波频率。平方后的载波信号变为连续的二倍频载波，其他两项频谱是噪声，n′(t)＝n²(t)+2ACDe^j2πft×n(t)。

2、无码载波频率估计技术

传统的捕获环节利用伪随机码良好的相关特性，完成对载波频率和码相位的粗略估计，实现跟踪通道初始化以及导航卫星识别。但是，消码处理后的信号削弱了信号强度，所以对该信号进行FFT变换，几乎检测不到峰值；而且在轨高动态条件下的多普勒频移搜索范围为几十kHz，若采用逐步搜索法，会耗费大量时间。所以载波频率初值的获取是需要解决的关键问题。

由于信号叠加可增强信号强度，以平方消码后的信号为例，以T为单位分成L个数据块进行叠加：

$\underset{i=0}{\overset{L}{\mathrm{\Σ}}}{r}^2(t+\mathrm{iT})=A^2{e}^{j2\mathrm{\π}2\mathrm{ft}}\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{e}^{j2\mathrm{\π}2\mathrm{fiT}}+\underset{i=1}{\overset{L-1}{\mathrm{\Σ}}}{n}^{\′}(t+\mathrm{iT})---\left(3\right)$

式中，i表示第i个数据块。

上式中在e^j2π2fiT＝1时可以取得最大值。但e^j2π2fiT与频率f有关，每个数据块的起始载波相位不同步，则信号直接累加不能保证每个数据块都能取到最大值。

让本地产生一个频率为f₁的正弦载波调整信号p(t)，对r²(t)进行相位补偿。相位补偿后的信号：

$\begin{array}{c}s\left(t\right)=r^2\left(t\right)p\left(t\right)\\ =r^2\left(t\right)e^{j2\mathrm{\π}{f}_{1}t}\\ =A^2{e}^{j2\mathrm{\π}2\mathrm{ft}}\·e^{j2\mathrm{\π}{f}_{1}t}+n^{\′}\left(t\right)e^{j2\mathrm{\π}{f}_{1}t}\\ =A^{2e^{j2\mathrm{\π}(2f+f_1)t}+n^{\′}\left(t\right)e^{j2\mathrm{\π}{f}_{1}t}}\end{array}---\left(4\right)$

其中，调整信号的引入只改变了载波频率，与噪声相乘依然是噪声。

仍然将相位补偿后的信号以T为单位分为L个数据块，并对L个数据块累加：

$\begin{array}{c}S\left(t\right)=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}s(t+\mathrm{iT})\\ =A^2{e}^{j2\mathrm{\π}(2f+f_1)t}\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{e}^{j2\mathrm{\π}(2f+f_1)\mathrm{iT}}+\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{n}^{\′}(t+\mathrm{iT})e^{j2\mathrm{\π}{f}_1(t+\mathrm{iT})}\\ =A^2{e}^{j2\mathrm{\π}{f}^{\′}t}\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{e}^{\mathrm{j\π}{f}^{\′}\mathrm{iT}}+\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{n}^{\′}(t+\mathrm{iT})e^{j2\mathrm{\π}{f}_1(t+\mathrm{iT})}\end{array}---\left(5\right)$

式中，S(t)为累加后的数据，f′＝2f+f₁。

当f′T＝M是一个整数时，才有

e^j2πf′iT＝e^j2πiM＝1       (6)

把使(6)式成立的f₁称作优化频率。若数据块的长度为1ms，即T＝1ms，则f′T是整数的条件就等同于f′是1kHz的整数倍，f₁的搜索范围是[0,1kHz]。

当(6)式成立时，累加后的数据可表示为：

$\begin{array}{c}S\left(t\right)=A^2{e}^{j2\mathrm{\π}{f}^{\′}t}\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{e}^{j2f^{\′}\mathrm{iT}}+\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{n}^{\′}(t-\mathrm{iT})e^{j2\mathrm{\π}{f}_1(t+\mathrm{iT})}\\ =A^2{e}^{j2\mathrm{\π}{f}^{\′}t}L+\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{n}^{\′}(t+\mathrm{iT})e^{j2\mathrm{\π}{f}_1(t+\mathrm{iT})}\end{array}---\left(7\right)$

平方信号经相位补偿和数据块累加得到了一个被放大了L倍的信号。对累加后的数据进行FFT运算，得到频谱分量的峰值，对平方法消码的信号而言，将尖峰频率值减去f₁再除以2，结果就对应某颗卫星信号的实际载波频率。

3、载波跟踪测量与恢复

中频数据消除了调制码，作为锁相环(即载波跟踪环)的输入信号，经过锁相环跟踪得到与输入信号同频的输出信号，由该输出信号的相位微分即可得到其频率。

锁相环技术的核心就是反馈控制环路，主要包括三部分：鉴相器、环路滤波器、数字压控振荡器。典型的锁相环结构如图5所示。

图5中，u_i(t)为输入信号，u_o(t)为压控振荡器输出的本地复现信号。通常情况下，环路滤波器为低通滤波器。锁相环的原理是通过鉴相器对比输入信号与本地信号的相位差值，然后由压控振荡器根据该差值实时闭环反馈调整本地复现信号，从而最终实现两者频率和相位的精确同步。当两者精确同步时，本地复现信号u_o(t)的频率和相位即为输入信号u_i(t)的频率和相位，即完成了对目标载波的跟踪测量。

4、导航卫星识别技术

在通常的卫星导航中，是通过伪码匹配的方式来识别导航卫星的，但是无码载波恢复技术仅能得到不同频移的单载波信号。由于用户卫星的初始轨道可认为是已知的，而导航卫星的星历可以通过地面上注，所以可以计算得到一系列多普勒频移的先验值。本发明将多普勒频移的实测值与先验值分别按照高低顺序排列成两组样本，分析了引起二者不一致的因素，设计了一种不受这些因素影响的匹配策略。

1)基于多普勒频移的导航星识别原理

用户卫星接收机测量的导航卫星的多普勒频移可以表示为

$f_d=\frac{1}{\mathrm{\λ}}\stackrel{\·}{R}+\mathrm{\δf}+\mathrm{\ϵ}---\left(8\right)$

式中，f_d为测量的多普勒频移，λ为载波波长；为用户卫星至导航卫星的伪距率真值；δf为用户卫星接收机时钟频率误差；ε为多普勒频移的测量误差，一般是白噪声。

已知用户卫星的初始轨道及地面周期上注的导航卫星的星历，可计算得到一系列先验多普勒频移：

${\hat{f}}_d=\frac{1}{\mathrm{\λ}}\stackrel{\·}{\hat{R}}=\frac{{(X-\hat{x})}^T(V-\hat{v})}{\mathrm{\λ}\hat{R}}---\left(9\right)$

式中，为多普勒频移先验值，和分别为用户卫星至导航卫星的含有用户星初轨误差的距离和距离变化率；X和V分别为导航卫星的位置和速度矢量，由上注的导航卫星星历计算获得；和分别为用户卫星的位置和速度矢量，由不精确的初始轨道计算得到。

由式(8)、(9)可知，不考虑先验根数及测量误差的影响，实测的多普勒频移f_d与先验估计的多普勒频移之间的频率差就是用户卫星接收机时钟频率误差。

假设某时刻有N颗可见导航卫星，将第n颗可见卫星信号的多普勒频移实测值f_d,n(n∈[1,N])及先验值分别从低到高排列。由于接收机钟差的存在，二者有固定的频率差，引入Δf作为的整体平移量，实现两组频率的一一对应。计算f_d,n与平移后的的频差绝对值，并求和得到匹配指标e，当e取最小值时认为两组频率已一一对应，该平移量就是接收机频差的估计值。e取最小值可以表示为：

$e_{\min }=\min \left(\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\right|f_{d,n}-{\hat{f}}_{d,n}-\mathrm{\Δf}\left|\right)---\left(10\right)$ 2)匹配计算方法

接收机时钟频率误差是有界的，所以整体平移量Δf的取值范围也是有限制的。假设待平移的频率范围为搜索步长为f_step，那么覆盖平移范围的所有平移点为：

$\mathrm{\Δf}=-\mathrm{\Δ}\hat{f}+(j-1)f_{\mathrm{step}},1\≤j\≤\frac{2\mathrm{\Δ}\hat{f}}{f_{\mathrm{step}}}+1---\left(11\right)$

式中，j表示第j次搜索。

对于上述遍历穷举式的计算方法，需要平移个频率点，而且随着搜索步长变小，匹配速度的速度越慢。

因为在两组频率几乎对齐时e有一个最小值，所以提出采用三分法进行匹配计算。算法如下：

(1)待平移的频率范围是选取区间的1/3、2/3处的坐标 $m_1=\mathrm{\Δ}{\hat{f}}_1+[{\mathrm{\Δ}\hat{f}}_2-\mathrm{\Δ}{\hat{f}}_1]/3,m_2=\mathrm{\Δ}{\hat{f}}_2-[\mathrm{\Δ}{\hat{f}}_2-\mathrm{\Δ}{\hat{f}}_1]/3,$ 把整个区间分成三个小区间[m₁,m₂]，使e最小的频率肯定会落在其中一个区间，三种可能的情况如图6的(a)、(b)、(c)所示；

(2)根据式(10)计算e(m₁)、e(m₂)，若e(m₁)＞e(m₂)，则舍弃区间令否则舍弃令

(3)重复执行1)、2)，直到小于某一阈值即停止搜索，获得使e最小的接收机频差估计值。

采用三分法的执行次数与逐步改变Δf相比，计算量大大减小，提高了导航星识别和接收机频差估计的效率，还可以根据需要调整终止条件，获得更高的匹配精度。

5、基于多普勒频移的高精度实时定轨算法

性能优良的定轨算法是实现高精度定轨的又一重要方面，由于多普勒频移测量值本质上是伪距率信息，所以需要研究动力学定轨算法，以扩展卡尔曼滤波器算法实现精密定轨。

采用卫星的位置和速度作为状态量，即

$\mathrm{\σ}={\left[\begin{array}{cccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T---\left(12\right)$

其中，σ是待估计的状态量，x,y,z为卫星的三轴位置，为卫星的三轴速度。

σ的变化规律可以用如下的卫星轨道动力学方程来描述：

$\stackrel{\·}{\mathrm{\σ}}=f(\mathrm{\σ},F_{\mathrm{\ϵ}})---\left(13\right)$

其中，F_ε为摄动力。多普勒测量为

$f_d=\frac{1}{\mathrm{\λ}}\·\frac{r_G-r_s}{|r_G-r_s|}\·({\mathrm{\υ}}_G-{\mathrm{\υ}}_s)+f+f_{\mathrm{\ϵ}}---\left(14\right)$

其中，f_d是多普勒频移，λ是波长，r_G是某导航卫星的位置矢量，r_s是某用户卫星的位置矢量，即x,y,z，υ_G是某导航卫星的速度矢量，υ_s是某用户卫星的速度矢量，即f_ε是接收机测量噪声，f是接收机的频移，可以表示为

$\stackrel{\·}{f}=n---\left(15\right)$

其中，n为接收机频差的驱动噪声。

将式(13)和式(15)作为状态方程，将式(14)作为测量方程，利用扩展卡尔曼滤波器实现轨道参数的精密估计，即可完成最终的定轨计算。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (6)

1.基于非合作导航卫星信号的无码多普勒定轨方法，其特征在于，包括以下步骤：

(1)对接收到的导航卫星信号经数字下变频后，采样得到数字中频信号，消除数字中频信号的保密调制码，根据消码后的信号计算得到载波频率的估计值；

(2)将步骤(1)得到的载波频率的估计值作为载波跟踪环跟踪的初始频率，产生载波跟踪环的本振信号，将步骤(1)中消码后的信号作为载波跟踪环的输入信号，恢复出载波信号，并测量提取出各个导航卫星信号的多普勒频移；

(3)计算各个导航卫星信号的先验值，结合步骤(2)得到的多普勒频移的实测值，建立含有接收机频差的多普勒频移实测值和先验值的识别匹配算法，完成导航星的识别与匹配；

(4)利用统计定轨算法实现轨道测量，利用轨道动力学方程建立状态方程，利用载波相位测量值建立相位差测量方程，再利用卡尔曼滤波器估计用户卫星的轨道参数。

2.根据权利要求1所述基于非合作导航卫星信号的无码多普勒定轨方法，其特征在于，步骤(1)中所述根据消码后的信号得到载波频率的估计值的具体步骤；

(A)选取合适的搜索步长，产生若干候选频率，对每一个候选频率产生一个载波相位调整信号；

(B)将步骤(A)产生的载波相位调整信号分别与消码后的信号相乘，对每一个相乘后的信号按照一定单位长度进行分块，累加分块数据，从而得到信号增强的新数据；

(C)对步骤(B)得到的新数据进行FFT运算，检查频谱结果，寻找超过门限的尖峰，对于每一个尖峰，记录其对应的频率值，即为载波频率的估计值。

3.根据权利要求1所述基于非合作导航卫星信号的无码多普勒定轨方法，其特征在于，步骤(3)的具体步骤：

(a)根据步骤(2)恢复出的载波信号，由用户卫星的初始轨道和周期上注的导航卫星星历计算出各导航卫星信号的多普勒频移先验值；

(b)将步骤(2)得到的多普勒频移实测值和步骤(a)得到的多普勒频移先验值分别按照高低顺序排列成两组数据样本，引入待求解的整体平移量，该整体平移量用于将多普勒频移先验值作整体平移；

(c)对各导航卫星信号的多普勒频移实测值与平移后的多普勒频移先验值的频差绝对值求和得到匹配指标；

(d)求解在匹配指标取最小值时的整体平移量，将多普勒频移先验值作整体平移，实现多普勒频移实测值与多普勒频移先验值的匹配，实现导航星识别。

4.根据权利要求3所述基于非合作导航卫星信号的无码多普勒定轨方法，其特征在于：在步骤(d)中，采用三分法求解在匹配指标取最小值时的整体平移量，其具体步骤：

(Ⅰ)假设某导航卫星信号的多普勒频移先验值作平移的频率区间是选取该区间的1/3、2/3处的坐标 $m_1=\mathrm{\Δ}{\hat{f}}_1+[\mathrm{\Δ}{\hat{f}}_2-\mathrm{\Δ}{\hat{f}}_1]/3,$ $m_2=\mathrm{\Δ}{\hat{f}}_2+[\mathrm{\Δ}{\hat{f}}_2-\mathrm{\Δ}{\hat{f}}_1]/3,$ 把整个区间分成三个小区间[m₁,m₂]，

(Ⅱ)分别计算在m₁、m₂处匹配指标e(m₁)、e(m₂)，若e(m₁)＞e(m₂)，则舍弃区间令 $\mathrm{\Δ}{\hat{f}}_1=m_1,$ 否则舍弃令 $\mathrm{\Δ}{\hat{f}}_2=m_2;$

(Ⅲ)重复步骤(Ⅰ)-(Ⅱ)，直到小于某一预设阈值时，即为整体平移量。

5.根据权利要求1所述基于非合作导航卫星信号的无码多普勒定轨方法，其特征在于：步骤(1)中，采用平方法消除数字中频信号的保密调制码。

6.根据权利要求1所述基于非合作导航卫星信号的无码多普勒定轨方法，其特征在于：步骤(1)中，采用交叉相关法消除数字中频信号的保密调制码。