[发明专利]一种基于自适应二阶卡尔曼滤波器的载波跟踪环路及其滤波方法

权利要求书

1.一种基于自适应二阶卡尔曼滤波器的载波跟踪环路，其特征在于：包括本地载波发生器、积分清除器、载波环鉴相器、环路滤波器和自适应二阶卡尔曼滤波器；

数字中频信号分别与本地载波发生器复现的载波及90°相移后的载波相乘后，再分别通过积分清除器提高载噪比，然后分别将每个积分清除器的相干积分结果均输出至载波环鉴相器，利用积分清除器输出的相干积分结果估算当前跟踪环路的载波相位差异，再将载波环鉴相器的输出载波相位差异作为自适二阶卡尔曼滤波器的观测值，经附带偏差的载波相位动态系统模型的自适应二阶卡尔曼滤波器递推出相位差异的估计值，相位差异的估计值输出至环路滤波器，经环路滤波器输出的控制信号使本地载波发生器实时调整。

2.根据权利要求1所述的一种基于自适应二阶卡尔曼滤波器的载波跟踪环路，其特征在于：所述的环路滤波器采用二阶环路滤波器，环路带宽为7Hz～10Hz。

3.根据权利要求1所述的一种基于自适应二阶卡尔曼滤波器的载波跟踪环路，其特征在于：所述的附带偏差的载波相位动态系统模型的状态方程具体为：

θe,kfd,kNCO=1-Δt01θe,k-1fd,k-1NCO+ΔtΔt2200fd,k-1fa,k-1+Wθ000---(1)]]>

fd,kfa,k=1Δt01fd,k-1fa,k-1+Wd00Wa---(2)]]>

附带偏差的载波相位动态系统模型的观测方程具体为：

θe,k-1mea=1-Δt2θe,k-1fd,k-1NCO+Δt2Δt26fd,k-1fa,k-1+Vk-1---(3)]]>

其中，θ_e为本地载波发生器NCO的输出与输入信号的载波相位误差，f_d为输入信号的多普勒频移，f_a是由GPS卫星与接收机沿着视线速度方向的加速度引起的频率移动变化率，为本地载波发生器复制的多普勒频率，Δt为更新周期，表示估计的载波相位误差，θ_e，k为t_k时刻本地载波发生器的输出与输入信号的载波相位差，为t_k时刻本地载波发生器复制的多普勒频率，为t_k-1时刻本地载波发生器复制的多普勒频率，f_d，k-1为t_k-1时刻输入信号的多普勒频移，f_a，k-1为t_k-1时刻由GPS卫星与接收机沿着视线速度方向的加速度引起的频率移动变化率，f_d，k为t_k时刻输入信号的多普勒频移，f_a，k为t_k时刻由GPS卫星与接收机沿着视线速度方向的加速度引起的频率移动变化率，W_n＝[W_θ，W_d，W_a]^T为噪声向量，W_θ、W_a、W_a为离散的噪声协方差分别是由接收机时钟引起的相位偏差、由时钟频率漂移引起的相位偏差、由卫星与接收机之间的加速度引起的相位偏差，表为t_k-1时刻估计的载波相位误差，θ_e，k-1为t_k-1时刻本地载波发生器NCO的输出与输入信号的载波相位差，V_k-1为观测噪声矢量。

4.一种基于自适应二阶卡尔曼滤波器的载波跟踪环路的滤波方法，其特征在于：具体包括以下几个步骤：

步骤一：定义自适应二阶卡尔曼滤波的初始状态以及计算过程中的中间变量：

定义初始状态为：

V0=P0xb(P0b)-1,]]>x0(+)=x0*-V0b0*---(4)]]>

b0(+)=b0*,]]>P0b(+)=P0b---(5)]]>

P0x(+)=P0x-V0P0bV0T---(6)]]>

其中和分别表示自适应二阶卡尔曼滤波初始状态向量x₀的均值和协方差，和分别表示自适应二阶卡尔曼滤波初始偏差向量b₀的均值和协方差，表示x₀和b₀的互协方差矩阵，V₀，为定义的混合矩阵；x₀(+)、b₀(+)、分别表示自适应无偏差滤波的初始状态估计值、自适应偏差滤波的初始状态估计值、自适应偏差滤波的初始估计均方误差矩阵、自适应无偏差滤波的初始估计均方误差矩阵；

定义中间变量：

N_k＝H_kU_k+D_k，Vk=Uk-KkxNk---(7)]]>

Uk=U‾k[I-λkbQk-1b[Pkb(-)]-1]---(8)]]>

U‾k=(Φk-1Vk-1+Bk-1)Ak-1-1---(9)]]>

uk=(U‾k+1-Uk+1)Akbk(+)---(10)]]>

Q‾kx=Qkx+Uk+1QkbU‾k+1---(11)]]>

其中H_k为观测矩阵，为自适应偏差滤波的滤波增益矩阵，I为单位矩阵，为自适应偏差滤波的渐消因子，表示系统噪声协方差，和表示偏差系统噪声协方差，为自适应偏差滤波的预测均方误差矩阵，Φ_k-1为状态转移矩阵，A_k-1，B_k-1，D_k为随机偏差系数矩阵，u_k，U_k，V_k，N_k，为定义的混合矩阵；b_k(+)表示自适应偏差滤波的状态估计值；

步骤二：计算自适应无偏差滤波过程中的状态一步预测值x_k(-)：

                        x_k(-)＝Φ_k-1x_k-1(+)+u_k-1                 (12)

其中，Φ_k-1为状态转移矩阵，x_k-1(+)为状态估计值，u_k-1为步骤一中定义的混合矩阵；

步骤三：计算自适应无偏差滤波方法中的预测均方误差矩阵

Pkx(-)=λkx(Φk-1Pk-1x(+)Φk-1T+Qk-1x)---(13)]]>

其中，为估计均方误差矩阵，为系统噪声协方差，为自适应无偏差滤波的渐消因子；

步骤四：计算自适应无偏差滤波的滤波增益矩阵

Kkx=Pkx(-)HkT[HkPkx(-)HKT+Rk]-1---(14)]]>

其中，R_k为观测噪声协方差；

步骤五：计算自适应无偏差滤波的估计均方误差矩阵

Pkx(+)=(I-KkxHk)Pkx(-)---(15)]]>

步骤六：计算自适应无偏差滤波的渐消因子

(1)计算k时刻观测量z_k的新息序列

ηkx=Zk-Hkxk(-)---(16)]]>

其中，H_k为观测矩阵，x_k(-)为自适应无偏差滤波的状态一步预测值；

(2)计算k时刻新息序列协方差的理论值：

t_k时刻自适应无偏差滤波的预测均方误差协方差矩阵为

Pkx(-)=λkx(Φk-1Pk-1x(+)Φk-1T+Qk-1x)---(17)]]>

其中Φ_k-1为状态转移矩阵，自适应无偏差滤波的为估计均方误差矩阵，为系统噪声的协方差，t_k为时间量，新息序列协方差的理论值为：

Ckx=E[ηkx(ηkx)T]=HkPkx(-)HkT+Rk---(18)]]>

其中E[·]表示求解协方差矩阵，R_k为量测噪声方差阵；

(3)计算k时刻新息序列协方差的估计值

当滤波增益矩阵为最优增益矩阵时，新息序列是白噪声序列，则新息序列的自相关函数等于零：

E[ηk+jx(ηkx)T]=Pkx(-)HkT-KkxCkx=0---(19)]]>

计算新息序列协方差的估计值

C‾kx=λk-1xηkx(ηkx)T1+λk-1x,k>112η0x(η0x)T,k=1---(20)]]>

其中，为k＝0时的新息序列，为k时刻的新息序列，为k-1时刻的渐消因子；

(4)计算自适应无偏差滤波的渐消因子

设预测均方误差协方差的估计值和新息序列协方差的估计值分别为：

P‾kx(-)=λkxPkx(-),]]>C‾kx=αkxCkx---(21)]]>

其中为自适应无偏差滤波的渐消因子，为自适应无偏差滤波的标量因子，和是根据标准卡尔曼滤波计算得到的理论值；最优增益阵的估计值为

K‾kx=P‾kx(-)HkTC‾k-1=[λkxPkx(-)]HkT[1αkxCk-1]=λkxαkx[Pkx(-)HkTCk-1]=λkxαkxKkx---(22)]]>

其中自适应无偏差滤波的标量因子为

αkx=max{1,trace(C‾kx)trace(Ckx)}---(23)]]>

其中trace(·)表示矩阵的迹运算；

新息序列协方差的估计值表示为

C‾kx=αkx(HkPkx(-)HkT+Rk)=λkxHkPkx(-)HkT+Rk---(24)]]>

由公式(24)得到

αkxHkPkx(-)HkT=λkxHkPkx(-)HkT+(1-αkx)Rk---(25)]]>

则自适应无偏差滤波的渐消因子为

λkx=max{1,trace(C‾kx)trace(Ckx)}---(26)]]>

步骤七：计算自适应无偏差滤波的状态估计值x_k(+)：

xk(+)=xk(-)+Kkxηkx---(27)]]>

步骤八：计算自适应偏差滤波的状态一步预测值b_k(-)：

                            b_k(-)＝A_k-1b_k-1(+)               (28)

其中，A_k-1为随机偏差系数矩阵，b_k-1(+)为自适应偏差滤波的状态估计值；

步骤九：计算自适应偏差滤波的预测均方误差矩阵

Pkb(-)=λkb(Ak-1Pk-1b(+)Ak-1T+Qk-1b)---(29)]]>

其中，为自适应偏差滤波的估计均方误差矩阵，为偏差系统噪声协方差，为自适应无偏差滤波的渐消因子；

步骤十：计算自适应偏差滤波的滤波增益矩阵

Kkb=Pkb(-)NkT[HkPkx(-)HxT+Rk+NkPkb(-)NKT]-1---(30)]]>

其中，R_k为观测噪声协方差，为自适应无偏差滤波的预测均方误差矩阵；

步骤十一：计算自适应偏差滤波的估计均方误差矩阵

Pkb(+)=[I-KkbNk]Pkb(-)---(31)]]>

步骤十二：计算自适应偏差滤波的渐消因子

自适应偏差滤波的渐消因子的计算方法和自适应无偏差滤波的渐消因子的计算方法相同，具体步骤为：

(1)计算k时刻观测量z_k的新息序列

ηkb=zk-Hkxk(-)-Nkbk(-)---(32)]]>

(2)计算k时刻新息序列协方差的理论值：

Ckb=HkPkx(-)HkT+Rk+NkPkb(-)NkT---(33)]]>

(3)计算k时刻新息序列协方差的估计值

C‾kb=λkbCkb=λk-1ηkb(ηkb)T1+λk-1,k>112η0b(η0b)T,k=1---(34)]]>

其中，为k＝0时的新息序列，为k时刻的新息序列，为k-1时刻的渐消因子；

(4)计算自适应偏差滤波的渐消因子

λkb=max{1,trace(C‾kb)trace(Ckb)}---(35)]]>

步骤十三：计算自适应偏差滤波的状态估计值b_k(+)：

bk(+)=bk(-)+Kkbηkb---(36)]]>步骤十四：计算自适应二阶卡尔曼滤波的状态一步预测值状态估计值预测均方误差矩阵和估计均方误差矩阵分别为：

其中x_k(-)表示自适应无偏差滤波的状态一步预测值；x_k(+)表示自适应无偏差滤波的状态估计值；b_k(-)表示自适应偏差滤波的状态一步预测值；b_k(+)表示自适应偏差滤波的状态估计值；自适应偏差滤波的预测均方误差矩阵；表示自适应无偏差滤波的估计均方误差矩阵；表示自适应偏差滤波的预测均方误差矩阵；表示自适应偏差滤波的估计均方误差矩阵；U_k和V_k为步骤一中定义的混合矩阵。