机群链路的前馈+反馈组合式载波跟踪方法

摘要

本发明涉及一种机群链路的前馈/反馈组合式载波跟踪方法，属于航空数据链、无线电导航技术领域。本发明提供了一种可以在电路板的数字信号处理器DSP和FPGA上实现机群链路高动态信号的精密跟踪与测量方法的体系构架，本发明提出高阶统计估计模型对环路多普勒频率参数进行跟踪估计，直接用于载波跟踪的迭代式前馈控制；同时，保留三阶PLL的为闭环载波跟踪结构，实现精确跟踪；扩频码精密跟踪则依赖载波辅助码环解决。本发明解决了传统意义上恶劣动态条件下传统高动态接收机频繁失锁、精度不佳的缺陷。本发明公开的方法能够广泛应用于基于抑制载波调制直接序列扩频体制的卫星导航接收机、测距系统和通信系统。

主权项

1.一种机群链路的前馈+反馈组合式载波跟踪方法，是在电路板的数字信号处理器DSP和FPGA上实现整个系统结构和方法，其特征在于：该方法是一种载波跟踪方式：前馈+反馈的载波跟踪控制策略，本方法的开环前馈控制环节基于统计滤波算法估计载波相位和多普勒频移跟踪残差，对当前第k步载波跟踪残差估计值迭代积分累加获得下一步第k+1步的多普勒频移值预报输出，获得高动态和全局线性特性；闭环反馈控制环节为三阶Costas PLL，将载波相位和频率跟踪残差估计结果用于闭环控制实现精密跟踪；前馈+反馈的载波跟踪控制策略，有如下恒等式：f_d(k)＝[f_FF(k)+f_FB(k)]+f_e(k)(1)＝f_u(k)+f_e(k)$f_d(k+1)=\underset{n=0}{\overset{N}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^n}{n!}+\frac{f_d^{(N+1)}\left(ϵ\right)·T^{N+1}}{(N+1)!}(ϵ\in [\mathrm{kT},(k+1)T\left]\right)---\left(2\right)$公式(2)是f_d(k+1)在第k点的泰勒展开式，f_d(k)与f_d(k+1)采样间隔为积分-清除周期T；其中，多普勒频移f_d(k)；多普勒频移跟踪残差f_e(k)；开环前馈补偿输出f_FF(k)；锁相环闭环控制输出f_FB(k)；令：为多普勒频移量的0阶～2阶导数的估计值，利用通过合适的统计估计模型和滤波算法获得；其中n＝0，1，2；$f_{\mathrm{FF}}(k+1)={\hat{f}}_d(k+1|k)=\underset{n=0}{\overset{N}{\Sigma }}\frac{{\hat{f}}_d^{\left(n\right)}\left(k\right)·T^n}{n!}---\left(3\right)$公式(3)是关于f_d(k+1)的一步递推预报公式，预报值作为第k+1步开环前馈控制补偿量，即：第k+1步的补偿控制输出量为其中，这里需要特别注意：是在第k步对f_d(k+1)的预报值，而不是在第k+1步对f_d(k+1)的估计值公式(3)的三个积分累加量需要通过状态估计获得，且舍去了泰勒展开式的余项因此预报值存在预报误差，根据公式(1)～公式(3)推导预报误差的解析形式：$\delta f_d(k+1|k)=f_d(k+1)-f_{\mathrm{FF}}(k+1)$$=f_d(k+1)-{\hat{f}}_d(k+1|k)$$=\underset{n=0}{\overset{N}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^n}{n!}+\frac{f_d^{(N+1)}\left(ϵ\right)·T^{N+1}}{(N+1)!}-\underset{n=0}{\overset{N}{\Sigma }}\frac{{\hat{f}}_d^{\left(n\right)}\left(k\right)·T^n}{n!}---\left(4\right)$$=\underset{n=0}{\overset{N}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)-{\hat{f}}_d^{\left(n\right)}\left(k\right)}{n!}·T^n+\frac{f_d^{(N+1)}\left(ϵ\right)·T^{N+1}}{(N+1)!}$$=\underset{n=0}{\overset{N}{\Sigma }}\frac{{\mathrm{\delta f}}_d^{\left(n\right)}\left(k\right)}{n!}·T^n+\frac{f_d^{(N+1)}\left(ϵ\right)·T^{N+1}}{(N+1)!}(N=2)$公式(4)看出，从第k步到第k+1步的预报误差δf_d(k+1|k)包括两类成分：第k步载波多普勒频移参数各阶导数(n＝0，1，2)的估计误差泰勒展开式截断误差预报误差δf_d(k+1|k)也就是第k+1步的开环前馈补偿控制的误差，如果没有PLL闭环反馈控制过程，则此预报误差就是第k+1步的多普勒频移跟踪残差，即：f_e(k+1)＝δf_d(k+1|k)                (5)若三阶PLL环路的闭环反馈控制作用为：f_FB(k+1)＝-δf_d(k+1|k)              (6)则由公式(4)得：$f_e(k+1)=f_d(k+1)-[f_{\mathrm{FF}}(k+1)+f_{\mathrm{FB}}(k+1)]$$=[f_d(k+1)-{\hat{f}}_d(k+1|k)]+f_{\mathrm{FB}}(k+1)$$\left(7\right)$$={\mathrm{\delta f}}_d(k+1|k)-{\mathrm{\delta f}}_d(k+1|k)$$=0$当公式(6)、公式(7)成立时，f_e(k+1)能够被三阶PLL环路通过闭环反馈控制过程消除；多普勒频移量的0阶～2阶导数的估计值作为多普勒频移参数测量值，即：载波跟踪环路的本地再生载波相对中频频点标称值的偏差成分，用于计算并输出开环前馈控制f_FF(k+1)，最大程度补偿了多普勒频移的恶劣动态，则锁相环承担的多普勒动态将被补偿干净而相当微弱，为预报误差δf_d(k+1|k)，通过反馈控制量f_FB(k)消除；所述的开环前馈补偿，是应用一种适用于恶劣条件下的载波参数统计估计方法实现，其中载波参数估计方法将积分-清除器的原始输出同相I支路和正交Q支路信号作为统计估计模型和滤波算法的输入量，该信号是接收通道A/D采样输出的数字中频信号与本地载波NCO正交混频解调、与本地码NCO驱动的再生伪码相关解扩，载波剥离、伪码剥离之后的信号再进入积分-清除器处理得到的输出信号；所述机群链路的前馈+反馈组合式载波跟踪方法实施前提是：【条件1】已完成载波频点捕获和码相位粗捕获，此时多普勒频移跟踪残差在-500Hz～+500Hz之内，码相位误差在1/4码片以内，但载波环路尚未进入锁定状态；或者跟踪失锁后刚刚完成信号重捕；【条件2】正常跟踪状态；此时，积分-清除器原始输出的同相I支路和正交Q支路信号在相关间隔末输出结果为：$\left\{\begin{array}{c}{I}_{\mathrm{ps}}\left(k\right)\approx A·R\left[ϵ\left(k\right)\right]·\sin c[{\mathrm{\Delta \omega }}_d\left(k\right)·N/2]\cos {\theta }_k+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\approx A·R\left[ϵ\left(k\right)\right]·\sin c[{\mathrm{\Delta \omega }}_d\left(k\right)·N/2]\sin {\theta }_k+n_Q\left(k\right)\end{array}\right.---\left(8\right)$公式(8)中：A为信号幅度；Δω_d(k)为多普勒频移估计残差，其中ε(k)为码相位(延时)估计偏差——真实延时和估计延时的差，ε(k)＝Δτ；R(·)为伪随机码理想的二电平自相关函数，均为时间的函数；N为积分清除器的积分点数；θ_k为载波相位误差，θ_k＝k·N·Δw_d(k)-Δw_d(k)·N/2+Δφ；n_I(k)，n_Q(k)为随机噪声；选择载波相位的0、1、2、3阶时间导数即载波相位、载波多普勒频移、载波多普勒一阶变化率、载波多普勒二阶变化率作为状态变量，系统控制输入为载波跟踪环频率控制输出值f_u(k)的积分多普勒相位差，选择积分-清除器输出的I、Q支路信号作为观测量，建立四阶非线性状态估计模型；(一)两种新型统计估计模型(i)建模基础当采样间隔积分-清除周期T足够小时，令θ_e(k)为积分-清除器输出信号的载波相位、(n＝0，1，2)为载波相位θ_e(k)的n+1阶导数，分别表示残余多普勒频移f_e(k)、残余多普勒频移变化率残余多普勒频移二阶变化率它们在采样间隔T内的迭代关系按照泰勒级数展开得：$\left\{\begin{array}{c}{\theta }_e(k+1)={\theta }_e\left(k\right)+\underset{n=0}{\overset{2}{\Sigma }}\frac{2\pi f_e^{\left(n\right)}\left(k\right)·T^{n+1}}{(n+1)!}+\frac{{2\mathrm{\pi f}}_e^{\left(3\right)}\left({ϵ}_1\right)·T^4}{24}\\ f_e(k+1)=\underset{n=0}{\overset{2}{\Sigma }}\frac{f_e^{\left(n\right)}\left(k\right)·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({ϵ}_2\right)·T^3}{6}\\ f_e^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\Sigma }}\frac{f_e^{\left(n\right)}\left(k\right)·T^{n-1}}{(n-1)!}+\frac{f_e^{\left(3\right)}\left({ϵ}_3\right)·T^2}{2}\\ f_e^{\left(2\right)}(k+1)=f_e^{\left(2\right)}\left(k\right)+f_e^{\left(3\right)}\left({ϵ}_4\right)·T\end{array}\right.---\left(9\right)$公式(9)中：(ε₁，ε₂，ε₃，ε₄∈[kT，(k+1)T])；(ii)两种系统动力学描述方式①描述方式1：根据公式(4)得：$f_e^{\left(n\right)}\left(k\right)=f_d^{\left(n\right)}\left(k\right)-f_u^{\left(n\right)}\left(k\right),$其中，n＝0，1，2，3                (10)定义第k步到第k+1步的载波NCO积分多普勒值第k步到第k+1步的载波相位差分值为：${\mathrm{\Delta \theta }}_u(k+1)={\theta }_u(k+1)-{\theta }_u\left(k\right)$$={\int }_{\mathrm{kT}}^{(k+1)T}2\pi f_u\left(t\right)\mathrm{dt}$$\left(11\right)$$=\underset{n=0}{\overset{2}{\Sigma }}\frac{{2\mathrm{\pi f}}_u^{\left(n\right)}\left(k\right)·T^{n+1}}{(n+1)!}+\frac{{2\mathrm{\pi f}}_u^{\left(3\right)}\left({\xi }_1\right)·T^4}{24}$将公式(10)、公式(11)代入公式(8)得：$\left\{\begin{array}{c}{\theta }_e(k+1)={\theta }_e\left(k\right)+\underset{n=0}{\overset{2}{\Sigma }}\frac{{2\mathrm{\pi f}}_d^{\left(n\right)}\left(k\right)·T^{n+1}}{(n+1)!}+\frac{{2\mathrm{\pi f}}_d^{\left(3\right)}\left({\xi }_1\right)·T^4}{24}-{\mathrm{\Delta \theta }}_u(k+1)\\ f_d(k+1)=\underset{n=0}{\overset{2}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^n}{n!}+\frac{f_d^{\left(3\right)}\left({\xi }_2\right)·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^{n-1}}{(n-1)!}+\frac{f_d^{\left(3\right)}\left({\xi }_3\right)·T^2}{2}\\ f_d^{\left(2\right)}(k+1)=\underset{n=2}{\overset{2}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^{n-2}}{(n-2)!}+f_d^{\left(3\right)}\left({\xi }_4\right)·T\end{array}\right.---\left(12\right)$公式(11)、公式(12)中：(ξ₁，ξ₂，ξ₃，ξ₄∈[kT，(k+1)T]；第k步到第k+1步的载波NCO积分多普勒值Δθ_u(k+1)作为系统状态方程的控制输入量，直接通过载波周数计数器和载波NCO相位寄存器计算获得；基于描述方式1建立的载波参数估计模型的系统方程为带控制信号的4阶系统；②描述方式2：将公式(10)、公式(12)代入公式(9)得：$\left\{\begin{array}{c}{\theta }_e(k+1)={\theta }_e\left(k\right)+\underset{n=0}{\overset{2}{\Sigma }}\frac{{2\mathrm{\pi f}}_d^{\left(n\right)}\left(k\right)·T^{n+1}}{(n+1)!}-\underset{n=0}{\overset{2}{\Sigma }}\frac{{2\mathrm{\pi f}}_u^{\left(n\right)}\left(k\right)·T^{n+1}}{(n+1)!}+\frac{2\pi f_e^{\left(3\right)}\left({ϵ}_1\right)·T^4}{24}\\ f_d(k+1)=\underset{n=0}{\overset{2}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^n}{n!}-\underset{n=0}{\overset{2}{\Sigma }}\frac{f_u^{\left(n\right)}\left(k\right)·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({ϵ}_2\right)·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\Sigma }}\frac{f_d^{\left(n\right)}\left(k\right)·T^{n-1}}{(n-1)!}-\underset{n=1}{\overset{2}{\Sigma }}\frac{f_u^{\left(n\right)}\left(k\right)·T^{n-1}}{(n-1)!}+\frac{f_e^{\left(3\right)}\left({ϵ}_3\right)·T^2}{2}\\ f_d^{\left(2\right)}(k+1)=f_d^{\left(2\right)}\left(k\right)-f_u^{\left(2\right)}\left(k\right)+f_e^{\left(3\right)}\left({ϵ}_4\right)·T\end{array}\right.---\left(13\right)$$\left\{\begin{array}{c}{f}_u(k+1)=\underset{n=0}{\overset{2}{\Sigma }}\frac{f_u^{\left(n\right)}\left(k\right)·T^n}{n!}+\frac{f_u^{\left(3\right)}\left({\zeta }_1\right)·T^3}{6}\\ f_u^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\Sigma }}\frac{f_u^{\left(n\right)}\left(k\right)·T^{n-1}}{(n-1)!}+\frac{f_u^{\left(3\right)}\left({\zeta }_2\right)·T^2}{2}\\ f_u^{\left(2\right)}(k+1)=f_u^{\left(2\right)}\left(k\right)+f_u^{\left(3\right)}\left({\zeta }_3\right)·T\end{array}\right.---\left(14\right)$根据公式(14)给出的关于f_u(k+1)、的状态方程，将第k+1步输出f_u(k+1)为观测量，建立估计模型并获得实时估计值作为控制输入带入公式(15)；(iii)量测方程积分-清除器原始输出的同相I支路和正交Q支路信号为量测向量，引用公式(8)：$\left\{\begin{array}{c}{I}_{\mathrm{ps}}\left(k\right)\approx A\left(k\right)·R\left[ϵ\left(k\right)\right]·\sin c[{\mathrm{\Delta \omega }}_d\left(k\right)·N/2]·\cos {\theta }_e\left(k\right)+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\approx A\left(k\right)·R\left[ϵ\left(k\right)\right]·\sin c[{\mathrm{\Delta \omega }}_d\left(k\right)·N/2]·\sin {\theta }_e\left(k\right)+n_Q\left(k\right)\end{array}\right.---\left(15\right)$z(k)＝h[x(k)]+n(k)＝[I_ps(k)，Q_ps(k)]^T+[n_I(k)，n_Q(k)]^T        (16)一般采用归一化的形式进行计算：$z\left(k\right)=h\left[x\left(k\right)\right]+n\left(k\right)=\left[\begin{array}{c}\cos {\theta }_e\left(k\right)\\ \sin {\theta }_e\left(k\right)\end{array}\right]+\left[\begin{array}{c}{n}_c\left(k\right)\\ n_s\left(k\right)\end{array}\right]=\left[\begin{array}{c}\cos (C·x\left(k\right))\\ \sin (C·x\left(k\right))\end{array}\right]+\left[\begin{array}{c}{n}_c\left(k\right)\\ n_s\left(k\right)\end{array}\right]---\left(17\right)$公式(17)中：量测向量为$z\left(k\right)={[\frac{I_{\mathrm{ps}}\left(k\right)}{\sqrt{I_{\mathrm{ps}}^2\left(k\right)+Q_{\mathrm{ps}}^2\left(k\right)}},\frac{Q_{\mathrm{ps}}\left(k\right)}{\sqrt{I_{\mathrm{ps}}^2\left(k\right)+Q_{\mathrm{ps}}^2\left(k\right)}}]}^T;$h[x(k)]＝[cosθ_e(k)，sinθ_e(k)]^T＝[cos(C·x(k))，sin(C·x(k))]^T；C＝[1，0，0，0]；n(k)＝[n_c(k)，n_s(k)]^T为归一化量测噪声向量，一般认为是高斯白噪声，其协方差阵为其中，SNR为积分-清除器输出信号的信噪比；I为二阶单位阵；(iv)建立两种载波参数估计模型①估计模型1根据公式(11)、公式(12)、公式(17)建立估计模型1；状态向量取为x(k)＝[θ_c(k)，f_d(k)，f′_d(k)，f″_d(k)]^T，控制输入量为u(k)＝Δθ_u(k+1)；状态方程表示为：x(k+1)＝Fx(k)+Gu(k)+w(k)                    (18)公式(18)中：状态转移矩阵$F=\left[\begin{array}{cccc}1& 2\mathrm{\pi T}& \pi T^2& {\mathrm{\pi T}}^3/3\\ 0& 1& T& T^2/2\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right];$控制输入耦合矩阵为G＝[-1，0，0，0]^T；将泰勒展开式余项看作系统噪声向量w(k)＝[w₁(k)，w₂(k)，w₃(k)，w₄(k)]^T；量测方程由公式(17)定义；②估计模型2根据公式(13)、公式(14)、公式(17)建立估计模型2；状态向量取为控制输入量为状态方程表示为：x(k+1)＝Fx(k)+Gu(k)+w(k)                    (19)公式(19)中：状态转移阵$F=\left[\begin{array}{cccc}1& 2\mathrm{\pi T}& \pi T^2& {\mathrm{\pi T}}^3/3\\ 0& 1& T& T^2/2\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right];$控制输入耦合矩阵为$G=\left[\begin{array}{ccc}2\mathrm{\pi T}& \pi T^2& \pi T^3/3\\ 1& T& T^2/2\\ 0& 1& T\\ 0& 0& 1\end{array}\right];$将泰勒展开式余项看作系统噪声向量w(k)＝[w₁(k)，w₂(k)，w₃(k)，w₄(k)]^T；量测方程由公式(17)定义；控制输入量为由公式(13)建立的线性估计模型通过线性kalman滤波获得；令状态向量为：根据公式(13)状态方程表示为：x_u(k+1)＝F_ux_u(k)+w_u(k)                (20)公式(20)中：状态转移阵$F_u=\left[\begin{array}{ccc}1& T& T^2/2\\ 0& 1& T\\ 0& 0& 1\end{array}\right];$将泰勒展开式余项看作系统噪声向量w_u(k)＝[w_u1(k)，w_u2(k)，w_u3(k)]^T；量测方程由公式(21)定义：z_u(k)＝C_ux_u(k)+n_u(k)                  (21)公式(21)中：量测耦合矩阵C_u＝[1，0，0，0]；n_u(k)为量测噪声，与频标的相位噪声和载波NCO位数有关，是微小量；两种估计模型中，由公式(18)和公式(17)定义的估计模型1是含标量控制输入的4阶非线性统计估计模型；由公式(19)和公式(17)定义的估计模型2是含3阶控制输入向量的4阶非线性统计估计模型，3阶控制输入向量由公式(20)和公式(11)定义的线性估计模型通过线性滤波实时计算。