CN106646543A - 基于主从式aukf算法的高动态卫星导航信号载波追踪方法 - Google Patents

Abstract

本发明公开了一种基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，包括：a、将接收的卫星信号经过伪码剥离后得到信号y，信号y分别与本地压控振荡器NCO产生的相位差为90°的两个复现载波混频得到信号I和信号Q；b、信号I和信号Q分别经过积分清零电路，得到信号y_I和信号y_Q作为观测量，并计算其向量形式Y(k)；c、Y(k)通过主从式AUKF算法处理，得到多普勒频移变化的相关状态向量d、通过环路滤波电路对频率估计值进行修正，得到修正值M(k)；e、用M(k)值控制本地压控振荡器NCO，使其产生的载波频率和接收信号载波频率保持一致。其能抑制滤波发散，使得在高动态下噪声统计特性发生变化时，载波追踪仍然具有很好的稳定性。

Description

基于主从式AUKF算法的高动态卫星导航信号载波追踪方法

技术领域

本发明涉及无线通信技术领域，特别涉及一种基于主从式AUKF算法的高动态卫星导航信号载波追踪方法。

背景技术

目前，有关高动态接收机载波跟踪研究可以归结为两个方面：外部其它传感器辅助载波跟踪和载波跟踪算法的改进。

一方面，外部其他传感器辅助载波跟踪典型的应用是惯导辅助。因为载波环本身是一种跟踪较为紧密的环路，所以，来自外界惯导系统的速度辅助信息必须及时和精确，否则会误导跟踪环路并最终导致信号失锁。然而，精确的惯性系统通常意味着昂贵的价格，并且外界对其参数校正的这一过程有可能相当复杂。所以，在没有惯导辅助时，研究高动态接收机的载波环路的跟踪算法也是当今热点研究课题之一。

另一方面，研究人员提出了多种载波环路跟踪算法。有研究者提出了FLL辅助PLL的载波追踪方法，即通过设定一个经验阈值让环路在两者之间切换；也有研究者提出了小波变换的载波追踪方法，即在PLL的鉴相器和环路滤波器之间加入小波变换，小波变换可以对鉴相器的输出进行噪声消除；还有人提出了一种模糊PLL载波跟踪结构，并利用梯度下降法和遗传算法优化模糊参数的载波追踪方法；以及提出了用无迹卡尔曼滤波(UKF)对载波频率进行估计的跟踪算法和基于平方根无迹卡尔曼滤波(SR-UKF)的载波跟踪算法。

载波跟踪算法中，UKF滤波算法是经典卡尔曼滤波(KF)的扩充，它是以无迹变换(UT变换)为基础，以卡尔曼线性滤波为框架，逼近非线性分布，在一步预测方程中，使用UT变换来处理均值和协方差的非线性传递的一种非线性滤波方法。相比扩展卡尔曼滤波(EKF)不需要线性化过程，不需要计算雅克比矩阵，具有比EKF更高的估计精度。

但是，UKF算法同KF和EKF—样，通常先验假设过程噪声和观测嗓声为统计特性已知不变的白噪声，这与高动态条件下卫星导航接收机的工作环境不符，因此不能保证载波跟踪系统的收敛性和稳定性，即在先验的噪声统计特性与实际的噪声统计特性不相符时，状态估计性能将变差甚至发散。

发明内容

基于此，针对上述问题，本发明的目的在于提供一种基于主从式自适应无迹卡尔曼滤波算法AUKF的高动态卫星导航信号载波追踪方法，其能自适应调整过程噪声方差，从而达到减小模型估计误差，抑制滤波发散的目的，使得在高动态下噪声统计特性发生变化时，载波追踪仍然具有很好的稳定性。

本发明的技术方案是：一种基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，包括以下步骤：

a、将接收的卫星信号经过伪码剥离后得到信号y，信号y分别与本地压控振荡器NCO产生的相位差为90°的两个复现载波混频得到信号I和信号Q；

b、所述信号I和信号Q分别经过积分清零电路，得到信号y_I和信号y_Q作为观测量，并计算信号y_I和信号y_Q的向量形式Y(k)；

c、所述向量Y(k)通过主从式AUKF滤波算法处理，得到多普勒频移变化的相关状态向量

d、所述状态向量通过环路滤波电路对频率估计值进行修正，得到修正值M(k)；

e、用所述M(k)值控制所述本地压控振荡器NCO，使其产生的载波频率和接收信号载波频率保持一致；

其中，k为离散时间，k＝1，2，3…。

作为对本发明的进一步改进，所述步骤b中的积分清零电路是以累加清零的方式实现的，每次累加的点数N为中频采样频率与载波跟踪环路的更新频率之比；当完成一次N点累加后，输出累加结果，并将累加寄存器清零；累加时钟为中频采样时钟，清零周期为环路的更新周期。

作为对本发明的进一步改进，所述步骤b中，对所述观测向量Y(k)，有：

其中，L^T＝[1 0 0 0]；X(k)＝[θ(k) ω₀(k) ω₁(k) ω₂(k)]^T，θ(k)为k时刻载波运动参数，ω₀(k)为k时刻载波多普勒角频率、ω₁(k)和ω₂(k)为k时刻多普勒角频率变换率和角频率变化率的变化率；A为接收信号的幅度，由于不是需要估计的参数，可设其为单位值；测量噪声向量n^T(k)＝[n_I(k) n_Q(k)]为零均值高斯白噪声，其协方差矩阵：

其中，n_I(k)表示所述信号I在k时刻的零均值高斯白噪声，n_Q(k)表示信号Q在k时刻的零均值高斯白噪声；n_I(l)表示所述信号I在l时刻的零均值高斯白噪声，n_Q(l)表示信号Q在l时刻的零均值高斯白噪声；I′为单位矩阵，σ²表示加加速度方差。

作为对本发明的进一步改进，所述接收信号离散形式下的状态向量X(k)是根据高机动目标跟踪跃动模型Jerk Model模型确定的，具体过程如下：

b11、定义在Jerk Model模型下，高动态运动目标加加速度满足如下指数衰减规律：

其中，α为衰减速率，γ为比例常数；w(t)为t时刻的激励白噪声，其方差即策动噪声强度为 为加加速度方差；j(t)为加加速度即加速度在单位时间内的变化量，t为连续时间，t＞0；

b12、根据所述步骤b11中的衰减规律，计算载波运动参数θ(t)分别与载波多普勒角频率ω₀(t)、多普勒角频率变换率ω₁(t)和角频率变化率的变化率ω₂(t)的关系式：

其中，θ(t)为载波信号中多普勒频移引起的t时刻相位变化；

b13、定义步骤b12中所述关系式在连续时间形式下的向量形式为：

其中,X(t)为连续形式下t时刻的状态向量，表示单位时间内X(t)的变化量，A和B为系数矩阵,并且有：

X(t)＝[θ(t) ω₀(t) ω₁(t) ω₂(t)]^T

B＝[0 0 0 1]^T

b14、根据步骤b13中所述的向量形式推导离散时间形式下的状态向量满足如下关系：

X(k+1)＝ΦX(k)+u(k)

其中，Φ为状态转移矩阵，u(k)为k时刻的策动噪声向量，并且有：

u(k)的协方差矩阵为：

其中，T为采样时间间隔，Φ(t_k+1-τ)表示(t_k+1-τ)时刻的状态转移矩阵，B(τ)表示τ时刻的系数矩阵，ω(τ)表示τ时刻的激励白噪声，q为策动噪声强度。

作为对本发明的进一步改进，所述步骤c中的主从式AUKF滤波算法具体包括如下循环过程：

c1、主UKF算法过程：用k-1时刻的输出状态向量以及步骤c2估计得出的k时刻的噪声强度来调整其时间更新式，得到状态向量状态向量和所述观测向量Y(k)用于测量更新式中，最后得到k时刻输出状态向量并同时计算新息协方差ν_k，进入步骤c2；

c2、从UKF算法过程：用k-1时刻的噪声强度来调整其时间更新式，得到噪声强度噪声强度和所述步骤c1产生的新息协方差ν_k用于其测量更新式中，最后得到k时刻的噪声强度再进入步骤c1。

作为对本发明的进一步改进，所述步骤c1中的时间更新式包括：

χ^* _k|k-1＝f(χ_k-1)

γ_k|k-1＝h(χ_k|k-1)

所述步骤c1中的测量更新式包括：

各参数存在如下关系：

λ＝n(ε²-1),

其中，x₀是初始状态向量，初始化：P₀为初始状态协方差矩阵，初始化：λ是比例因子；n表示扩维后的系统状态向量维数；ε是散步程度因子，通常取正值；β为验前分布因子，对于高斯分布，取2最优；为一阶统计特性时的加权系数；为二阶统计特性时的加权系数；Q^x是先验过程噪声协方差矩阵Q的修正值，其噪声强度q^x由从UKF估计得到，以逼近Q的噪声强度；R^x是假设已知的观测噪声协方差矩阵R。

作为对本发明的进一步改进，所述步骤c2中的时间更新式包括：

P_k|k-1＝Q^q

所述步骤c2中的测量更新式包括：

各参数存在如下关系：

其中，Q^q和R^q分别是过程噪声协方差矩阵和测量噪声协方差矩阵；为一阶统计特性时的加权系数；为二阶统计特性时的加权系数。

作为对本发明的进一步改进，所述步骤c1中的时间更新式也可为以下简化式：

所述步骤c1中的测量更新式也可为以下简化式：

其中，各参数存在如下关系：

作为对本发明的进一步改进，所述步骤d中的环路滤波电路满足如下特性：

其中，H(z)是环路滤波电路的转移函数，K₀K_d是环路增益，w_n是自然频，ξ是阻尼比，T_S是预检测积分时间。

本发明的有益效果是：

(1)观测向量Y(k)进入主UKF模块，以主UKF模块产生的新息协方差ν_k为依据，采用从UKF模块实时估计策动噪声强度q并自适应调整主UKF的策动噪声矢量协方差矩阵Q，使其减小模型估计误差，抑制滤波发散；

(2)即便在先验的噪声统计特性与实际的噪声统计特性不相符时，也能通过自适应调整过程噪声方差，以保持良好的状态估计性能，保证载波跟踪系统的收敛性和稳定性；

(3)其能自适应调整过程噪声方差，从而达到减小模型估计误差，抑制滤波发散的目的，使得在高动态下噪声统计特性发生变化时，载波追踪仍然具有很好的稳定性；

(4)通过对主UKF算法过程中各计算式进行合理的简化，以适应在FPGA中的应用，极大地提升了运行速率并降低了硬件成本。

附图说明

图1是将本发明实施例所述基于主从式AUKF算法的高动态卫星导航信号载波追踪方法应用于高动态追踪环路的电路原理图；

图2是主从式AUKF滤波算法的流程框图；

图3是采用传统UKF滤波算法和本发明主从式AUKF滤波算法追踪多普勒频移对比图。

具体实施方式

下面结合附图对发明的实施例进行详细说明。

实施例

参考图1和图2，一种基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，包括以下步骤：

a、将接收的卫星信号经过伪码剥离后得到信号y，信号y分别与本地压控振荡器NCO产生的相位差为90°的两个复现载波混频得到信号I和信号Q；

b、所述信号I和信号Q分别经过积分清零电路，得到信号y_I和信号y_Q作为观测量，并计算信号y_I和信号y_Q的向量形式Y(k)；

c、所述向量Y(k)通过主从式AUKF算法处理，得到多普勒频移变化的相关状态向量

d、所述状态向量通过环路滤波电路对频率估计值进行修正，得到修正值M(k)；

e、用所述M(k)值控制所述本地压控振荡器NCO，使其产生的载波频率和接收信号载波频率保持一致；

其中，k为离散时间，k＝1，2，3…。

基于上述方法的高动态追踪环路如图1，主要包括本地压控振荡器NCO、积分清零电路、主从式AUKF算法模块和环路滤波电路。

在另一个实施例中，所述步骤b中的积分清零电路是以累加清零的方式实现的，每次累加的点数N为中频采样频率与载波跟踪环路的更新频率之比；当完成一次N点累加后，输出累加结果，并将累加寄存器清零；累加时钟为中频采样时钟，清零周期为环路的更新周期。

在另一个实施例中，所述步骤b中，对所述观测向量Y(k)，有：

其中，A为接收信号的幅度，测量噪声向量n^T(k)＝「n_I(k) n_Q(k)]为零均值高斯白噪声，其协方差矩阵：

设定L^T＝[1 0 0 0],θ(k)＝L^TX(k)，

则测量方程为：

其中，A为接收信号的幅度，由于不是需要估计的参数，可设其为单位值；X(k)＝[θ(k) ω₀(k) ω₁(k) ω₂(k)]^T，θ(k)为k时刻载波运动参数，ω₀(k)为k时刻载波多普勒角频率、ω₁(k)和ω₂(k)为k时刻多普勒角频率变换率和角频率变化率的变化率；n_I(k)表示所述信号I在k时刻的零均值高斯白噪声，n_Q(k)表示信号Q在k时刻的零均值高斯白噪声；n_I(l)表示所述信号I在l时刻的零均值高斯白噪声，n_Q(l)表示信号Q在l时刻的零均值高斯白噪声；I′为单位矩阵，σ²表示加加速度方差。

在另一个实施例中，所述接收信号离散形式下的状态向量X(k)是根据高机动目标跟踪跃动模型Jerk Model模型确定的，具体过程如下：

b11、定义在Jerk Model模型下，高动态运动目标加加速度满足如下指数衰减规律：

其中，α为衰减速率，γ为比例常数；w(t)为t时刻的激励白噪声，其方差即策动噪声强度为 为加加速度方差；j(t)为加加速度即加速度在单位时间内的变化量，t为连续时间，t＞0；

b12、根据所述步骤b11中的衰减规律，计算载波运动参数θ(t)分别与载波多普勒角频率ω₀(t)、多普勒角频率变换率ω₁(t)和角频率变化率的变化率ω₂(t)的关系式：

其中，θ(t)为载波信号中多普勒频移引起的t时刻相位变化；

b13、定义步骤b12中所述关系式在连续时间形式下的向量形式为：

其中,X(t)为连续形式下t时刻的状态向量，表示单位时间内X(t)的变化量，A和B为系数矩阵,并且有：

X(t)＝[θ(t) ω₀(t) ω₁(t) ω₂(t)]^T

B＝[0 0 0 1]^T

b14、根据步骤b13中所述的向量形式推导离散时间形式下的状态向量满足如下关系：

X(k+1)＝ΦX(k)+u(k)

其中，Φ为状态转移矩阵，u(k)为k时刻的策动噪声向量，并且有：

u(k)的协方差矩阵为：

其中，T为采样时间间隔，Φ(t_k+1-τ)表示(t_k+1-τ)时刻的状态转移矩阵，B(τ)表示τ时刻的系数矩阵，ω(τ)表示τ时刻的激励白噪声，q为策动噪声强度。

在另一个实施例中，所述步骤c具体包括如下循环过程：

c1、主UKF算法过程：用k-1时刻的输出状态向量以及步骤c2估计得出的k时刻的噪声强度来调整其时间更新式，得到状态向量状态向量和所述观测向量Y(k)用于测量更新式中，最后得到k时刻输出状态向量并同时计算新息协方差ν_k，进入步骤c2；

c2、从UKF算法过程：用k-1时刻的噪声强度来调整其时间更新式，得到噪声强度噪声强度和所述步骤c1产生的新息协方差ν_k用于其测量更新式中，最后得到k时刻的噪声强度再进入步骤c1。

在另一个实施例中，所述步骤c1主UKF算法过程中，Sigma采样和时间更新的计算公式如下：

χ^* _k|k-1＝f(χ_k-1) (4.2)

γ_k|k-1＝h(χ_k|k-1) (4.6)

所述步骤c1中的测量更新式包括：

各参数存在如下关系：

其中，x₀是初始状态向量，初始化：P₀为初始状态协方差矩阵，初始化：λ是比例因子；n表示扩维后的系统状态向量维数；ε是散步程度因子，通常取正值；β为验前分布因子，对于高斯分布，取2最优；为一阶统计特性时的加权系数；为二阶统计特性时的加权系数；Q^x是先验过程噪声协方差矩阵Q的修正值，其噪声强度q^x由从UKF估计得到，以逼近Q的噪声强度；R^x是假设已知的观测噪声协方差矩阵R。

在另一个实施例中，所述步骤c2从UKF算法过程中，Sigma采样和时间更新的计算公式如下：

P_k|k-1＝Q^q (4.14)

所述步骤c2中的测量更新式如下：

各参数存在如下关系：

其中，Q^q和R^q分别是过程噪声协方差矩阵和测量噪声协方差矩阵；为一阶统计特性时的加权系数；为二阶统计特性时的加权系数。

在另一个实施例中，还提供对步骤c1中的时间更新式和测量更新式进行简化的方法，以期能实现在FPGA中的应用，简化过程如下：

对χ_k|k-1中的协方差矩阵进行Cholesky分解，设定：

联合式(4.5)-(4.6)、(7.1)-(7.3)对式(4.7)进行简化，得到结果如下：

其中，

a＝W₁+18W₂+2cW₂,

由于，

另外，设：

将式(7.5)、式(7.6)代入到式(4.8)中，得到观测值协方差矩阵如下：

其中，

W₄＝W₆(1-a)²,

W₅＝2W₂,

W₆＝W₃+18W₂,

W₇＝W₁+18W₂,

m＝W₅sin²(c),

n＝2W₂[cos(c)-a]²,

l＝m-n,

关于状态值与观测值协方差矩阵的简化结果如下：

其中，

综合上述，整个主UKF滤波算法简化后可以表示如下，加权系数和初始化：

Sigma采样和时间更新：

测量更新：

由式(7.10)和式(7.11)可以可知，要使能在FPGA中用相应模块实现，主要包括4×4维矩阵的加法、减法，正弦余弦的运算，4×4维矩阵与4×1维矩阵的乘法运算，4×1维矩阵与1×2维矩阵的乘法运算，4×2维矩阵与2×2维矩阵的乘法运算，4×2维矩阵与2×1维矩阵的乘法运算，4×2维矩阵与2×4维矩阵的乘法运算，Cholesky分解等。

在另一个实施例中，所述步骤d中的环路滤波电路满足如下特性：

其中，H(z)是环路滤波电路的转移函数，K₀K_d是环路增益，w_n是自然频，ξ是阻尼比，T_S是预检测积分时间。

仿真过程中，采用传统UKF滤波算法和本发明主从式AUKF滤波算法追踪多普勒频移对比图如图3，关于仿真部分参数如下。

主UKF中参数设置：

载波信号的中频频率：l575.42MHz；

采样时间：T＝lms，载噪比：CNR＝32dB-Hz；

状态初始值：x₀＝[0 32995.519 323.356 0]；

初始状态协方差和观测噪声协方差矩阵分别为：

P₀＝diag([0.001 0.001 0.001 0.001])

R＝diag{1×10^-9 1×10^-9}；

UT变换中，取ε＝1.027，κ＝0，β＝2。

从UKF中参数设置：

初始噪声强度：q₀＝5；

初始噪声强度协方差：P_q0＝1e-8；

噪声强度过程协方差：Q^q＝1e-10；

观测噪声协方差：R^q＝diag{[1e-8 1e-8]}；

UT变换中，取ε＝1.027，κ＝0，β＝2。

以上所述实施例仅表达了发明的具体实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离发明构思的前提下，还可以做出若干变形和改进，这些都属于发明的保护范围。

Claims (9)

1.基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，包括以下步骤：

a、将接收的卫星信号经过伪码剥离后得到信号y，信号y分别与本地压控振荡器NCO产生的相位差为90°的两个复现载波混频得到信号I和信号Q；

b、所述信号I和信号Q分别经过积分清零电路，得到信号y_I和信号y_Q作为观测量，并计算信号y_I和信号y_Q的向量形式Y(k)；

c、所述向量Y(k)通过主从式AUKF滤波算法处理，得到多普勒频移变化的相关状态向量

d、所述状态向量通过环路滤波电路对频率估计值进行修正，得到修正值M(k)；

e、用所述M(k)值控制所述本地压控振荡器NCO，使其产生的载波频率和接收信号载波频率保持一致；

其中，k为离散时间，k＝1，2，3…。

2.根据权利要求1所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于：所述步骤b中的积分清零电路是以累加清零的方式实现的，每次累加的点数N为中频采样频率与载波跟踪环路的更新频率之比；当完成一次N点累加后，输出累加结果，并将累加寄存器清零；累加时钟为中频采样时钟，清零周期为环路的更新周期。

3.根据权利要求2所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于：所述步骤b中，对所述观测向量Y(k)，有：

$Y\left(k\right)=\left[\begin{array}{c}{y}_I\left(k\right)\\ y_Q\left(k\right)\end{array}\right]=\left[\begin{array}{c}A\sin \left(L^{T}X\right(k\left)\right)\\ A\cos \left(L^{T}X\right(k\left)\right)\end{array}\right]+n\left(k\right)$

其中，L^T＝[1 0 0 0]；X(k)＝[θ(k) ω₀(k) ω₁(k) ω₂(k)]^T，θ(k)为k时刻载波运动参数，ω₀(k)为k时刻载波多普勒角频率、ω₁(k)和ω₂(k)为k时刻多普勒角频率变换率和角频率变化率的变化率；A为接收信号的幅度，由于不是需要估计的参数，可设其为单位值；测量噪声向量为零均值高斯白噪声，其协方差矩阵：

$R=E\[n\left(k\right)n^T\left(l\right)\]=\left\{\begin{array}{cc}0& k\≠l\\ {\mathrm{\σ}}^2{I}^{\′}& k=l\end{array}\right.$

其中，n_I(k)表示所述信号I在k时刻的零均值高斯白噪声，n_Q(k)表示信号Q在k时刻的零均值高斯白噪声；n_I(l)表示所述信号I在l时刻的零均值高斯白噪声，n_Q(l)表示信号Q在l时刻的零均值高斯白噪声；I′为单位矩阵，σ²表示加加速度方差。

4.根据权利要求3所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，所述接收信号离散形式下的状态向量X(k)是根据高机动目标跟踪跃动模型Jerk Model模型确定的，具体过程如下：

b11、定义在Jerk Model模型下，高动态运动目标加加速度满足如下指数衰减规律：

$\frac{d}{dt}j\left(t\right)=\mathrm{\γ}\×(-\mathrm{\α}j(t)+w(t\left)\right)$

其中，α为衰减速率，γ为比例常数；w(t)为t时刻的激励白噪声，其方差即策动噪声强度为 为加加速度方差；j(t)为加加速度即加速度在单位时间内的变化量，t为连续时间，t＞0；

b12、根据所述步骤b11中的衰减规律，计算载波运动参数θ(t)分别与载波多普勒角频率ω₀(t)、多普勒角频率变换率ω₁(t)和角频率变化率的变化率ω₂(t)的关系式：

$\frac{d}{dt}\mathrm{\θ}\left(t\right)={\mathrm{\ω}}_0\left(t\right)$

$\frac{d}{dt}{\mathrm{\ω}}_0\left(t\right)={\mathrm{\ω}}_1\left(t\right)$

$\frac{d}{dt}{\mathrm{\ω}}_1\left(t\right)={\mathrm{\ω}}_2\left(t\right)$

$\frac{d}{dt}{\mathrm{\ω}}_2\left(t\right)=-{\mathrm{\α\ω}}_2\left(t\right)+w\left(t\right)$

其中，θ(t)为载波信号中多普勒频移引起的t时刻相位变化；

b13、定义步骤b12中所述关系式在连续时间形式下的向量形式为：

$\stackrel{\·}{X}\left(t\right)=AX\left(t\right)+Bw\left(t\right)$

其中,X(t)为连续形式下t时刻的状态向量，表示单位时间内X(t)的变化量，A和B为系数矩阵,并且有：

X(t)＝[θ(t) ω₀(t) ω₁(t) ω₂(t)]^T

$A=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 0& -\mathrm{\α}\end{array}\right]$

B＝[0 0 0 1]^T

b14、根据步骤b13中所述的向量形式推导离散时间形式下的状态向量满足如下关系：

X(k+1)＝ΦX(k)+u(k)

其中，Φ为状态转移矩阵，u(k)为k时刻的策动噪声向量，并且有：

$\mathrm{\Φ}=\left[\begin{array}{cccc}1& T& \frac{T^2}{2}& \frac{T^3}{6}\\ 0& 1& T& \frac{T^2}{2}\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]$

$u\left(k\right)={\∫}_{t_k}^{t_{k+1}}\mathrm{\Φ}(t_{k+1}-\mathrm{\τ})B\left(\mathrm{\τ}\right)w\left(\mathrm{\τ}\right)d\mathrm{\τ},$

u(k)的协方差矩阵为：

$Q\left(k\right)=q\left[\begin{array}{cccc}\frac{T^7}{252}& \frac{T^6}{72}& \frac{T^5}{30}& \frac{T^4}{24}\\ \frac{T^6}{72}& \frac{T^5}{20}& \frac{T^4}{8}& \frac{T^3}{6}\\ \frac{T^5}{30}& \frac{T^4}{8}& \frac{T^3}{3}& \frac{T^2}{2}\\ \frac{T^4}{24}& \frac{T^3}{6}& \frac{T^2}{2}& T\end{array}\right]$

其中，T为采样时间间隔，Φ(t_k+1-τ)表示(t_k+1-τ)时刻的状态转移矩阵，B(τ)表示τ时刻的系数矩阵，ω(τ)表示τ时刻的激励白噪声，q为策动噪声强度。

5.根据权利要求3或4所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，所述步骤c中的主从式AUKF滤波算法具体包括如下循环过程：

c1、主UKF算法过程：用k-1时刻的输出状态向量以及步骤c2估计得出的k时刻的噪声强度来调整其时间更新式，得到状态向量状态向量和所述观测向量Y(k)用于测量更新式中，最后得到k时刻输出状态向量并同时计算新息协方差ν_k，进入步骤c2；

c2、从UKF算法过程：用k-1时刻的噪声强度来调整其时间更新式，得到噪声强度噪声强度和所述步骤c1产生的新息协方差ν_k用于其测量更新式中，最后得到k时刻的噪声强度再进入步骤c1。

6.根据权利要求5所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，所述步骤c1中的时间更新式包括：

${\mathrm{\χ}}_{k-1}=\[{\hat{x}}_{k-1},{\hat{x}}_{k-1}+\sqrt{n+\mathrm{\λ}}{P}_{x_{k-1}},{\hat{x}}_{k-1}-\sqrt{n+\mathrm{\λ}}{P}_{x_{k-1}}\]$

χ^* _k|k-1＝f(χ_k-1)

${\widehat{\mathrm{\χ}}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\mathrm{\χ}}_{i,k|k-1}^{*}$

$P_{x_k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{i,k|k-1}){({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{i,k|k-1})}^T+Q^x$

${\mathrm{\χ}}_{k|k-1}=\[{\hat{x}}_{k|k-1},{\hat{x}}_{k|k-1}+\sqrt{n+\mathrm{\λ}}{P}_{x_{k|k-1}},{\hat{x}}_{k|k-1}-\sqrt{n+\mathrm{\λ}}{P}_{x_{k|k-1}}\]$

γ_k|k-1＝h(χ_k|k-1)

${\hat{y}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\mathrm{\γ}}_{i,k|k-1}$

所述步骤c1中的测量更新式包括：

$P_{y_k{y}_k}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\mathrm{\γ}}_{i,k|k-1}-{\hat{y}}_{k|k-1}){({\mathrm{\γ}}_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T+R^x$

$P_{x_k{y}_k}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\mathrm{\χ}}_{i,k|k-1}-{\hat{x}}_{k|k-1}){({\mathrm{\γ}}_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$K_{x_k}=P_{x_k{y}_k}{P}_{y_k{y}_k}^{-1}$

${\hat{x}}_k={\hat{x}}_{k|k-1}+K_{x_k}(y_k-{\hat{y}}_{k|k-1})$

$P_{x_k}=P_{x_k|k-1}-K_{x_k}{P}_{y_k{y}_k}{K}_{x_k}^T$

各参数存在如下关系：

$w_0^m=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},$

$w_0^c=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(n-{\mathrm{\ϵ}}^2+\mathrm{\β}),$

$w_i^m=w_i^c=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},i=1,...2n,$

$\mathrm{\η}=\sqrt{n+\mathrm{\λ}},$

λ＝n(ε²-1),

其中，x₀是初始状态向量，初始化：P₀为初始状态协方差矩阵，初始化：λ是比例因子；n表示扩维后的系统状态向量维数；ε是散步程度因子，通常取正值；β为验前分布因子，对于高斯分布，取2最优；为一阶统计特性时的加权系数；为二阶统计特性时的加权系数；Q^x是先验过程噪声协方差矩阵Q的修正值，其噪声强度q^x由从UKF估计得到，以逼近Q的噪声强度；R^x是假设已知的观测噪声协方差矩阵R。

7.根据权利要求6所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，所述步骤c2中的时间更新式包括：

${\hat{q}}_{k|k-1}={\hat{q}}_{k-1}$

P_k|k-1＝Q^q

${\hat{S}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{w}_{qi}^c{\mathrm{\ψ}}_{i,k|k-1}$

所述步骤c2中的测量更新式包括：

$P_{S_k{S}_k}=\underset{i=0}{\overset{2n}{\Σ}}{w}_{qi}^c({\mathrm{\ψ}}_{i,k|k-1}-{\hat{S}}_{k|k-1}){({\mathrm{\ψ}}_{i,k|k-1}-{\hat{S}}_{k|k-1})}^T+R^q$

$K_{q_k}=P_{q_k{S}_k}{P}_{S_k{S}_k}^{-1}$

${\hat{q}}_k={\hat{q}}_{k|k-1}-K_{q_k}({\hat{S}}_k-{\hat{S}}_{k|k-1})$

$P_{q_k}=P_{q_k|k-1}-K_{q_k}{P}_{S_k{S}_k}{K}_{q_k}^{-1}$

各参数存在如下关系：

$w_{qi}^m=w_{qi}^c=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},i=1,...2n$

${\hat{q}}_0=E\[q_0\]$

$P_{{\mathrm{qx}}_0}=E\[(q_0-{\stackrel{\‾}{q}}_0){(q_0-{\stackrel{\‾}{q}}_0)}^T\]$

其中，Q^q和R^q分别是过程噪声协方差矩阵和测量噪声协方差矩阵；为一阶统计特性时的加权系数；为二阶统计特性时的加权系数。

8.根据权利要求6或7所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，所述步骤c1中的时间更新式也可为以下简化式：

$\left\{\begin{array}{c}{\hat{x}}_{k|k-1}=\mathrm{\Φ}{\hat{x}}_{k-1}\\ P_{x_k|k-1}={\mathrm{\ΦP}}_{x_{k-1}}{\mathrm{\Φ}}^T+Q^x\\ S_{x_k|k-1}=chol\left(P_{x_k|k-1}\right)\\ {\hat{y}}_{k|k-1}=\left[\begin{array}{c}\sin \left(b\right)\\ \cos \left(b\right)\end{array}\right]\end{array}\right.$

所述步骤c1中的测量更新式也可为以下简化式：

$\left\{\begin{array}{c}{P}_{y_k{y}_k}=\left[\begin{array}{cc}(W_4+l){\sin }^2\left({\hat{x}}_{1,k|k-1}\right)+m+R_{11}& (W_4+l)\sin \left({\hat{x}}_{1,k|k-1}\right)\cos \left({\hat{x}}_{1,k|k-1}\right)\\ (W_4+l)\sin \left({\hat{x}}_{1,k|k-1}\right)\cos \left({\hat{x}}_{1,k|k-1}\right)& (W_4+l){\cos }^2\left({\hat{x}}_{1,k|k-1}\right)+m+R_{22}\end{array}\right]\\ P_{x_k{y}_k}=2\sqrt{n+\mathrm{\λ}}\×W_2\sin \left(c\right)\left[\begin{array}{c}{s}_{11,x_k|k-1}\\ s_{21,x_k|k-1}\\ s_{31,x_k|k-1}\\ s_{41,x_k|k-1}\end{array}\right]\[\begin{array}{cc}\cos \left(b\right)& -\sin \left(b\right)\end{array}\]\\ K_{x_k}=P_{x_k{y}_k}{P}_{y_k{y}_k}^{-1}\\ {\hat{x}}_k={\hat{x}}_{k|k-1}+K_{x_k}(y_k-{\hat{y}}_{k|k-1})\end{array}\right.$

其中，各参数存在如下关系：

$\left\{\begin{array}{c}{W}_1=w_0^m=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},\\ W_2=w_i^m=w_i^c=\frac{\mathrm{\λ}}{2(n+\mathrm{\λ})},i=1,...2n\\ W_3=w_0^c=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(n-{\mathrm{\ϵ}}^2+\mathrm{\β})\\ b={\hat{x}}_{1,k|k-1}\\ {\hat{x}}_0=E\[x_0\]\\ P_{x_0}=E\[(x_0-{\stackrel{\‾}{x}}_0){(x_0-{\stackrel{\‾}{x}}_0)}^T\]\end{array}\right.$

9.根据权利要求1所述的基于主从式AUKF算法的高动态卫星导航信号载波追踪方法，其特征在于，所述步骤d中的环路滤波电路满足如下特性：

$G_1=\frac{1}{K_0{K}_d}\frac{4{\left(w_n{T}_S\right)}^2}{4+4{\mathrm{\ξw}}_n{T}_S+{\left(w_n{T}_S\right)}^2}$

$G_2=\frac{1}{K_0{K}_d}\frac{8{\mathrm{\ξw}}_n{T}_S}{4+4{\mathrm{\ξw}}_n{T}_S+{\left(w_n{T}_S\right)}^2}$

$H\left(z\right)=\frac{K_0{K}_d(G_1+G_2)z^{-1}-K_0{K}_d{G}_2{z}^{-2}}{1+K_0{K}_d(G_1+G_2)z^{-1}-(1-K_0{K}_d{G}_2)z^{-2}}$

其中，H(z)是环路滤波电路的转移函数，K₀K_d是环路增益，w_n是自然频，ξ是阻尼比，T_S是预检测积分时间。