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

Abstract

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

Description

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

技术领域

本发明涉及一种机群链路的前馈/反馈组合式载波跟踪方法，属于航空数据链、无线电导航技术领域。 

背景技术

目前，在航空数据链、无线电导航技术领域内，机群链路在恶劣条件下信号捕获、跟踪与测量存在一定困难，各种因素引起的恶劣条件对扩频信号捕获跟踪都会带来影响。这里从三个角度定义恶劣条件： 

1)运动动态恶劣：①两运动体相对运动位移(零阶导数)对时间高阶可导，且各阶导数如：速度(一阶导数)、加速度(二阶导数)、加加速度/冲击加速度(三阶导数)等变化较大；②位移n阶不可导/n+1阶不连续(n＝2，3，…)，即存在跳变点，如脉冲式推力控制启动/停止时加速度量最可能出现跳变点；③位移的n阶导数(n＝2，3，…)在较小的时间尺度内(如：毫秒、亚毫秒量级)发生幅度变化将导致n+1阶导数对应时刻的时间导数幅度很大； 

2)信号动态条件恶劣：①由于姿态变化引起的天线指向变化剧烈，导致馈入天线的信号功率大范围变化；②距离变化、发射机内部因素变化、接收机内部因素引起接收信号功率的短时间内突然变化，反映到时间导数较大；③外界干扰因素，诸如日凌干扰、敌对干扰、自然/非自然干扰等因素引起带内噪声突然增大； 

3)归入恶劣条件的其他因素(新体制引起的捕获/跟踪困难)：①采用加密序列控制扩频码型切换的加密跳码扩频调制，码型切换、弱信号低信噪比(匿踪需求)导致捕获困难、失锁概率增大，对捕获速度/准确度和跟踪精度/稳定性、失锁概率等提出苛刻的要求；②扩频测控/非扩频数传等时复用体制中模式切换引起捕获、跟踪困难、失锁概率增大，给载波变化参数和扩频码变化参数的预报带来很大困难。 

对于给定的捕获方法和跟踪环路结构、参数，从前面定义的恶劣条件出发分析恶劣条件对捕获和跟踪的影响如下： 

1)对于大动态变化范围的多普勒频移，捕获带加宽，引起频点搜索时间变长； 

2)对于较短时间内多普勒频移发生突变，产生大幅值的多普勒频移高阶导数，载波跟踪系统对高阶量跟踪能力不足会引起对多普勒频移高阶导数量的跟踪误差，导致跟踪误差增大、失锁概率高； 

3)姿态机动或者突发干扰引起的短时间信噪比突变，使得接收单元环路滤波器无法适应窄带接收条件下的高动态、宽带条件下的低信噪比，导致跟踪误差增大、失锁概率增高； 

4)采用加密跳码扩频调制体制和扩频测控/非扩频数传等时复用体制对载波频点、扩频码相位的捕获带来很大困难，一旦失锁、失捕则很难直接恢复捕获跟踪，采用引导捕获的模式不利于隐蔽或者提高通信质量，因此需要高精度的载波跟踪/码跟踪效果以获得高精度的测量结果用于降低失捕/失锁概率和实施载波变化参数和扩频码变化参数的准确预报用于辅助捕获和辅助跟踪。 

在机群组网任务中，将有可能面临这些恶劣条件的影响，如：机群编队构型变换时的编队机动、姿态机动、天线对准、信道质量恶劣、弱信号低信噪比，等等。因此要求跟踪环路能够适应大范围动态应力和信号应力、瞬间短时中断时能够利用运动参数外推控制跟踪环路开环跟踪直到信号恢复。这里还是需要重点考虑窄带跟踪的优势、困难与关键技术，然而高动态、低信噪比、高精度跟踪(RMS方差)条件下，跟踪环路的带宽大、方差小、窄带跟踪三者互为矛盾，引入载波辅助码环措施、开环频率跟踪估计和前馈/反馈跟踪环设计获得窄带高阶跟踪能力，引入实时载波平滑伪距算法抑制测量方差获得准确的无差测量结果。 

可见，提高载波多普勒频率估计与跟踪精度极为重要，除了决定跟踪精度、失锁概率，更有利于计算机群成员节点之间的相对速度、加速度等测量值，以及在加密跳码扩频调制体制和扩频测量/非扩频数传等时复用体制中用于外推预报相关参数对载波频点、扩频码相位的捕获。 

发明内容

本发明的目的在于提供一种机群链路的前馈/反馈组合式载波跟踪方法，本发明提出高阶统计估计模型对环路多普勒频率参数进行跟踪估计，直接用于载波跟踪的迭代式前馈控制；同时，保留三阶PLL的为闭环载波跟踪结构，实现精确跟踪。可以如此理解：基于跟踪误差统计滤波(无偏或者有偏估计)的迭代式前馈补偿获得高动态的跟踪能力，补偿了信号的高动态；基于三阶PLL的反馈式闭环载波跟踪消除补偿高动态信号后的控制残差(低动态、 低偏差)，获得高精度的跟踪载波跟踪。扩频码精密跟踪则依赖载波辅助码环解决。 

本发明一种机群链路的前馈/反馈组合式载波跟踪方法，可以在电路板的数字信号处理器DSP和FPGA上实现整个系统结构和算法。 

本发明提出一种全新的载波跟踪方式：前馈+反馈的载波跟踪控制策略。本方法的开环前馈控制环节基于统计滤波算法估计载波相位和多普勒频率跟踪残差，对当前(第k步)载波跟踪残差估计值迭代积分累加获得下一步(第k+1步)的多普勒频率值预报输出，获得高动态和全局线性特性；闭环反馈控制环节与为三阶Costas PLL，利用载波相位和频率跟踪残差估计结果用于闭环控制实现精密跟踪。可以看出，这种策略利用了高精度的载波跟踪误差估计结果和递推预报计算实现比较严酷条件下的开环跟踪，前馈补偿多普勒频移中高动态的成分，保留了三阶PLL用于窄带条件下快速消除稳态误差成分。 

本发明一种机群链路的前馈/反馈组合式载波跟踪方法，具体包括开环前馈控制环节和闭环反馈控制环节；其原理图及结构框图如图2、图1所示； 

前馈+反馈的载波跟踪控制策略，有如下恒等式： 

f_d(k)＝[f_FF(k)+f_FB(k)]+f_e(k) 

     ＝f_u(k)+f_e(k)                                  (1) 

$f_d(k+1)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{(n+1)}\left(\mathrm{\ϵ}\right)\·T^{(n+1)}}{(n+1)!}(\mathrm{\ϵ}\∈[\mathrm{kT},(k+1)T\left]\right)---\left(2\right)$

公式(2)是f_d(k+1)在f_d(k)的泰勒展开式，f_d(k)与f_d(k+1)采样间隔为积分-清除周期T_o其中，多普勒频移f_d(k)；多普勒频移跟踪残差f_e(k)；开环前馈补偿输出f_FF(k)；锁相环闭环控制输出f_FB(k)； 

令： ${\hat{f}}_d^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 为多普勒频移量的0阶～2阶导数的估计值，利用 

通过合适的统计估计模型和滤波算法获得。 

$f_{\mathrm{FF}}(k+1)={\hat{f}}_d(k+1|k)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\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步的补偿控制输出量为 $f_{\mathrm{FF}}(k+1)={\hat{f}}_d(k+1|k)$ (这里需要特别注意： 

是在第k步对f_d(k+1)的预报值，而不是在第k+1步对f_d(k+1)的估计值 

切勿混淆)。公式(3)的三个积分累加量需要通过状态估计获得，且舍去了泰勒展开式的余项 

因此预报值 

存在预报误差，根据公式(1)～公式(3)推导预报误差的解析形式： 

$\mathrm{\δ}{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}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{(n+1)}\left(\mathrm{\ϵ}\right)\·T^{n+1}}{(n+1)!}-\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{{\hat{f}}_d^{\left(n\right)}\left(k\right)\·T^n}{n!}$

$=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\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(\mathrm{\ϵ}\right)\·T^{n+1}}{(n+1)!}$

$=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{\δ}{f}_d^{\left(n\right)}\left(k\right)}{n!}\·T^n+\frac{f_d^{(n+1)}\left(\mathrm{\ϵ}\right)\·T^{n+1}}{(n+1)!}(N=2)---\left(4\right)$

公式(4)看出，从第k步到第k+1步的预报误差δf_d(k+1|k)包括两类成分：第k步载波多普勒频率参数(各阶导数) ${\hat{f}}_d^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 的估计误差 $\mathrm{\δ}{f}_d^{\left(n\right)}\left(k\right)=f_d^{\left(n\right)}\left(k\right)-{\hat{f}}_d^{\left(n\right)}\left(k\right);$ 泰勒展开式截断误差 

预报误差δ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)$

$=\mathrm{\δ}{f}_d(k+1|k)-\mathrm{\δ}{f}_d(k+1|k)$

$=0---\left(7\right)$

当公式(6)、公式(7)成立时，f_e(k+1)能够被三阶PLL环路通过闭环反馈控制过程消除。 

多普勒频移量的0阶～2阶导数的估计值 ${\hat{f}}_d^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,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)\≈A\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\cos {\mathrm{\θ}}_k+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\≈A\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\sin {\mathrm{\θ}}_k+n_Q\left(k\right)\end{array}\right.---\left(8\right)$

公式(8)中：A为信号幅度；Δω_d(k)为多普勒频移估计残差， $\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)={\mathrm{\ω}}_d\left(k\right)-{\widehat{\mathrm{\ω}}}_d\left(k\right);$ ε(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支路信号作为观测量，建立四阶非线性状态估计模型。一般情况，可选择EKF、UKF这类非线性滤波方法实现载波参数的估计。 

(一)两种新型统计估计模型 

(i)建模基础 

当采样间隔(积分-清除周期)T足够小时，令θ_e(k)为积分-清除器输出信号的载波相位、  ${\mathrm{\ω}}_e^{\left(n\right)}\left(k\right)=2\mathrm{\π}{f}_e^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 为载波相位θ_e(k)的n+1阶导数(分别表示残余多普勒频率f_e(k)、残余多普勒频率变化率f_e ⁽¹⁾(k)、残余多普勒频率二阶变化率f_e ⁽²⁾(k))，它们在采样间隔T内的迭代关系可以按照泰勒级数展开得： 

$\left\{\begin{array}{c}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_e^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_e^{\left(3\right)}\left({\mathrm{\ϵ}}_1\right)\·T^4}{24}\\ f_e(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_e^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_2\right)\·T^3}{6}\\ f_e^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_e^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_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({\mathrm{\ϵ}}_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=\mathrm{0,1,2,3})---\left(10\right)$

定义第k步到第k+1步的载波NCO积分多普勒值(第k步到第k+1步的载波相位差分值)为： 

$\mathrm{\Δ}{\mathrm{\θ}}_u(k+1)={\mathrm{\θ}}_u(k+1)-{\mathrm{\θ}}_u\left(k\right)$

$={\∫}_{\mathrm{kT}}^{(k+1)T}2\mathrm{\π}{f}_u\left(t\right)\mathrm{dt}$

$=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_u^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_u^{\left(3\right)}\left({\mathrm{\ξ}}_1\right)\·T^4}{24}---\left(11\right)$

将公式(10)、公式(11)代入公式(8)得： 

$\left\{\begin{array}{c}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_d^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_d^{\left(3\right)}\left({\mathrm{\ξ}}_1\right)\·T^4}{24}-\mathrm{\Δ}{\mathrm{\θ}}_u(k+1)\\ f_d(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ξ}}_2\right)\·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ξ}}_3\right)\·T^2}{2}\\ f_d^{\left(2\right)}(k+1)=\underset{n=2}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-2}}{(n-2)!}+f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_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}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_d^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}-\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_u^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_e^{\left(3\right)}\left({\mathrm{\ϵ}}_1\right)\·T^4}{24}\\ f_d(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}-\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_2\right)\·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}-\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_n^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ϵ}}_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_d^{\left(3\right)}\left({\mathrm{\ϵ}}_4\right)\·T\end{array}\right.---\left(13\right)$

$\left\{\begin{array}{c}{f}_u(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_u^{\left(3\right)}\left({\mathrm{\ζ}}_1\right)\·T^3}{6}\\ f_u^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_u^{\left(3\right)}\left({\mathrm{\ζ}}_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({\mathrm{\ζ}}_3\right)\·T\end{array}\right.---\left(14\right)$

根据公式(14)给出的关于f_u(k+1)、f_u ⁽¹⁾(k+1)、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)\≈A\left(k\right)\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\·\cos {\mathrm{\θ}}_e\left(k\right)+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\≈A\left(k\right)\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\·\sin {\mathrm{\θ}}_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 {\mathrm{\θ}}_e\left(k\right)\\ \sin {\mathrm{\θ}}_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为归一化量测噪声向量，一般认为是高斯白噪声，其协方差阵为 $R=E\left[n\left(k\right)n^T\left(k\right)\right]={\mathrm{\σ}}_n^{2}I,$ ${\mathrm{\σ}}_n^2={10}^{-\mathrm{SNR}}$ (SNR为积分-清除器输出信号的信噪比)；I为二阶单位阵。 

(iv)建立两种载波参数估计模型 

①估计模型1 

根据公式(11)、公式(12)、公式(17)建立估计模型1。状态向量取为x(k)＝[θ_e(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{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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\left(k\right)={[{\mathrm{\θ}}_e\left(k\right),f_d\left(k\right),{\stackrel{\·}{f}}_d\left(k\right),{\stackrel{\·\·}{f}}_d\left(k\right)]}^T,$ 控制输入量为 $u\left(k\right)={[{\hat{f}}_u\left(k\right),{\hat{f}}_u^{\′}\left(k\right),{\hat{f}}_u^{\′\′}\left(k\right)]}^T.$ 状态方程可表示为： 

x(k+1)＝Fx(k)+Gu(k)+w(k)        (19) 

公式(19)中：状态转移阵 $F=\left[\begin{array}{cccc}1& 2\mathrm{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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)定义。 

控制输入量为 $u\left(k\right)={[{\hat{f}}_u\left(k\right),{\hat{f}}_u^{\′}\left(k\right),{\hat{f}}_u^{\′\′}\left(k\right)]}^T$ 由公式(13)建立的线性估计模型通过线性kalman滤波获得。令状态向量为：x_u(k)＝[f_u(k)，f_u′(k)，f_u″(k)]^T，根据公式(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)为量测噪声，与频标的相位噪声(晶振的短期稳定度或铷原子频标Allan方差)和载波NCO位数(对于32bit以上的计数器，可以忽略)有关，是微小量。 

两种估计模型中，由公式(18)和公式(17)定义的估计模型1是含标量控制输入的4阶非线性统计估计模型；由公式(19)和公式(17)定义的估计模型2是含3阶控制输入向量的4阶非线性统计估计模型，3阶控制输入向量可由公式(20)和公式(11)定义的线性估计模型通过线性滤波实时计算。 

本发明一种机群链路的前馈/反馈组合式载波跟踪方法，其优点在于： 

首先，公式(18)和公式(17)构成的载波参数的非线性统计估计模型1、由公式(19)和公式(17)构成的载波参数的非线性统计估计模型2，配合相应的非线性滤波算法(如：UKF滤波)构成载波参数的非线性统计估计方法具有明显的优点： 

①状态空间信息量丰富，提供载波相位的0～3阶导数； 

②对残余多普勒频移高阶导数建模，模型阶数更高，描述运动更准确，保留高动态成分信息； 

③滤波算法可实时递推获得全状态估计。 

④具有4阶无差估计精度，且跟踪偏差为泰勒展开式余项，与多普勒频移的3阶导数成正比、与积分-清除周期的三次方成正比，量级很小，能够适用于比较极端的恶劣动态条件； 

⑤应用约束小，只需完成捕获的情况下残余多普勒频移过大、扩频码相位尚未完全对准的低相关输出、低信噪比输出时即可快速进入跟踪状态，基于开环前馈控制过程迅速引导进入精密跟踪过程实施闭环反馈跟踪过程； 

⑥估计模型和统计滤波算法对大信号动态(信噪比SNR短时间大范围变化)不敏感、对强运动动态(多普勒频移短时间大范围变化)不敏感，估计模型鲁棒性强，能支持窄带跟踪，恶劣条件下能保持优良的接收灵敏度、接收信噪比、动态跟踪精度、失锁概率； 

⑦采用前馈控制大大补偿了接收信号的动态，基本完成解调环节载波剥离工作，使得进入积分-清除器的多普勒频移残余成分大大降低，在滤波递推周期内积分-清除器输出信号的相位变化很小(例如：积分-清除周期T＝1ms时采样频率为f＝1kHz，对于≤10Hz的残余多普勒频移，远离Nequist频率，在步长范围内载波相位变化≤3.6°)，由于非线性动力学过程变化范围很小时呈线性动力学特性，非线性滤波容易获得更高的性能和稳定度； 

⑧计算复杂度不高，实用性好，实施成本低，在现有平台上(DSP、FPGA)容易编程(C语言、VHDL语言、Verilog语言)实现 

常规的载波参数估计算法建立的非线性估计模型对真实物理过程进行了简化，忽略了一些因素，但这些因素可能在恶劣条件下凸显出来影响估计效果甚至导致结果的不可信不可用，甚至滤波发散；然而本发明给出的新型模型推导过程完全基于恒等式推导，对信号的动力学过程没有进行任何近似，因此性能相对得到提高。特别指出：公式(18)和公式(17)构成的载波参数的非线性统计估计模型1、公式(19)和公式(17)构成的载波参数的非线性统计估计模型2利用非线性滤波算法精确地估计出了载波多普勒频移参数和载波相位参数，能够直接 作为载波相位和载波多普勒参数测量值用于其他用途，如：测速、测角、测姿态、基线变化量测量，等等。 

附图说明

图1所示为本发明前馈+反馈的载波跟踪控制方法结构框图 

图2所示为本发明前馈+反馈的载波跟踪环路原理框图 

具体实施方式

本发明一种机群链路突发扩频信号的载波跟踪与参数估计方法，可以在电路板的数字信号处理器DSP和FPGA上实现整个系统结构和算法。 

一、前馈/反馈的载波跟踪控制方法 

本发明提出一种全新的载波跟踪方式：前馈+反馈的载波跟踪控制策略。本方法的开环前馈控制环节基于统计滤波算法估计载波相位和多普勒频率跟踪残差，对当前(第k步)载波跟踪残差估计值迭代积分累加获得下一步(第k+1步)的多普勒频率值预报输出，获得高动态和全局线性特性；闭环反馈控制环节与为三阶Costas PLL，利用载波相位和频率跟踪残差估计结果用于闭环控制实现精密跟踪。可以看出，这种策略利用了高精度的载波跟踪误差估计结果和递推预报计算实现比较严酷条件下的开环跟踪，前馈补偿多普勒频移中高动态的成分，保留了三阶PLL用于窄带条件下快速消除稳态误差成分。 

图1给出了前馈+反馈的载波跟踪控制策略结构框图；图2给出了前馈+反馈的载波跟踪环路原理框图。重点讨论开环前馈控制环节和闭环反馈控制环节的关系、开环前馈控制环节的算法原理及输入输出特性。图中：①接收信号中频频率f_lFd(k)＝f_IF+f_d(k)；②中频频率标称值f_IF；③多普勒频移f_d(k)；④多普勒频移跟踪残差f_e(k)；⑤积分-清除器输出的信号噪声n_e(k)；⑥多普勒频移跟踪残差量测值 ${\tilde{f}}_e\left(k\right)=f_e\left(k\right)+n_e\left(k\right);$ ⑦开环前馈补偿输出f_FF(k)；⑧锁相环闭环控制输出f_FB(k)；⑨载波跟踪环频率控制输出f_u(k)。 

根据图1、图2描述的前馈+反馈的载波跟踪控制策略，有如下恒等式： 

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}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{(n+1)}\left(\mathrm{\ϵ}\right)\·T^{(n+1)}}{(n+1)!}(\mathrm{\ϵ}\∈[\mathrm{kT},(k+1)T\left]\right)---\left(2\right)$

公式(2)是f_d(k+1)在f_d(k)的泰勒展开式，f_d(k)与f_d(k+1)采样间隔为积分-清除周期T。 

令： ${\hat{f}}_d^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 为多普勒频移量的0阶～2阶导数的估计值，利用 

通过合适的统计估计模型和滤波算法获得。 

$f_{\mathrm{FF}}(k+1)={\hat{f}}_d(k+1|k)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\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步的补偿控制输出量为 $f_{\mathrm{FF}}(k+1)={\hat{f}}_d(k+1|k)$ (这里需要特别注意： 

是在第k步对f_d(k+1)的预报值，而不是在第k+1步对f_d(k+1)的估计值 

切勿混淆)。公式(3)的三个积分累加量需要通过状态估计获得，且舍去了泰勒展开式的余项 

因此预报值 

存在预报误差，根据公式(1)～公式(3)推导预报误差的解析形式： 

$\mathrm{\δ}{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}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{(n+1)}\left(\mathrm{\ϵ}\right)\·T^{n+1}}{(n+1)!}-\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{{\hat{f}}_d^{\left(n\right)}\left(k\right)\·T^n}{n!}$

$=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\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(\mathrm{\ϵ}\right)\·T^{n+1}}{(n+1)!}$

$=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{\δ}{f}_d^{\left(n\right)}\left(k\right)}{n!}\·T^n+\frac{f_d^{(n+1)}\left(\mathrm{\ϵ}\right)\·T^{n+1}}{(n+1)!}(N=2)---\left(4\right)$

公式(4)看出，从第k步到第k+1步的预报误差δf_d(k+1|k)包括两类成分：第k步载波多普勒频率参数(各阶导数) ${\hat{f}}_d^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 的估计误差 $\mathrm{\δ}{f}_d^{\left(n\right)}\left(k\right)=f_d^{\left(n\right)}\left(k\right)-{\hat{f}}_d^{\left(n\right)}\left(k\right);$ 泰勒展开式截断误差 

预报误差δ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)$

$=\mathrm{\δ}{f}_d(k+1|k)-\mathrm{\δ}{f}_d(k+1|k)$

$=0---\left(7\right)$

当公式(6)、公式(7)成立时，f_e(k+1)能够被三阶PLL环路通过闭环反馈控制过程消除。由于开环前馈补偿控制已充分补偿了恶劣条件引起的应力动态与信号动态，闭环反馈控制的PLL环路能够采用窄带跟踪技术将一定范围低频扰动完全消除。因此，即使前馈控制环节的载波参数估计存在较小的非零均值偏差，只要在PLL闭环系统的正常工作范围内，都能彻底被三阶PLL的闭环控制过程消除达到很高的精度。 

多普勒频移量的0阶～2阶导数的估计值 ${\hat{f}}_d^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 作为多普勒频移参数测量值，即：载波跟踪环路的本地再生载波相对中频频点标称值的偏差成分， 

用于计算 

并输出开环前馈控制f_FF(k+1)，最大程度补偿了多普勒频移的恶劣动态，则锁相环承担的多普勒动态将被补偿干净而相当微弱(为预报误差δf_d(k+1|k))，可通过反馈控制量f_FB(k)消除，由此可以得出本章最重要的结论： 

【结论6.1】开环前馈补偿+闭环反馈控制的载波跟踪模式中，选择合理的统计估计模型和滤波算法，不需要满足最优统计滤波及无偏估计所需的苛刻约束要求，在估计存在一定范围偏差情况下，仍然能够获得恶劣条件下的高精度载波跟踪与测量，满足误差均值为零和输出方差极小化。 

【结论6.2】载波跟踪三阶锁相环的闭环反馈控制环节与开环前馈补偿环节工作相互独立，PLL采用反正切鉴相器直接对积分-清除器输出I/Q支路信号鉴相，而不应利用任何前馈补偿环节统计滤波算法获得的载波参数估计值，达到利用反馈调节原理完全消除前馈环节估计/预报误差引起的开环跟踪偏差。 

该结论给出的命题指出了开环前馈补偿环节采用的统计估计模型和滤波算法的选择原 则： 

①优越的动态跟踪特性，估计误差动态足够小，能够确保窄带PLL闭环反馈跟踪环节正常工作和满足跟踪精度； 

②估计误差界能够确保PLL闭环反馈跟踪环节正常工作和满足跟踪精度； 

③估计结果不要求最优、无偏，甚至零均值； 

④统计估计建模和滤波算法带来的计算复杂度不高，能满足高更新率的估计输出(一般选择步长为积分-清除周期，即：扩频码周期)。 

二、载波参数的新型统计估计算法 

对于高动态环境下的直序扩频测距体制接收机载波跟踪环路频率估计算法，建立准确、合理的载波参数运动行为描述模型是至关重要的。 

本发明提出并验证了一种适用于恶劣条件下的载波参数统计估计方法，用于开环前馈补偿。本发明中，载波参数估计方法的将积分-清除器的原始输出同相I支路和正交Q支路信号作为统计估计模型和滤波算法的输入量，该信号是接收通道A/D采样输出的数字中频信号与本地载波NCO正交混频解调、与本地码NCO驱动的再生伪码相关解扩和将载波剥离、伪码剥离之后的进入积分-清除器处理得到的输出。本发明的方法实施前提是： 

【条件1】已完成载波频点捕获和码相位粗捕获，此时多普勒频移跟踪残差在-500Hz～+500Hz之内，码相位误差在1/4码片以内，但载波环路尚未进入锁定状态；或者跟踪失锁后刚刚完成信号重捕； 

【条件2】正常跟踪状态。 

此时，积分-清除器原始输出的同相I支路和正交Q支路信号在相关间隔末输出结果为： 

$\left\{\begin{array}{c}{I}_{\mathrm{ps}}\left(k\right)\≈A\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\cos {\mathrm{\θ}}_k+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\≈A\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\sin {\mathrm{\θ}}_k+n_Q\left(k\right)\end{array}\right.---\left(8\right)$

公式(8)中：A为信号幅度；Δω_d(k)为多普勒频移估计残差， $\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)={\mathrm{\ω}}_d\left(k\right)-{\widehat{\mathrm{\ω}}}_d\left(k\right);$ ε(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支路信号作为观测量，建立四阶非线性状态估计模型。一般情况，可选择EKF、UKF这类非线性滤波方法实现载波参数的估计。 

(一)两种新型统计估计模型 

(i)建模基础 

在估计频率的同时估计相位可以使频率估计误差减小，而为适合高动态环境的需要及提高相位和频率的估计精度，将频率的一、二阶导数一并估计是合适的。当采样间隔(积分-清除周期)T足够小时，令θ_e(k)为积分-清除器输出信号的载波相位、 ${\mathrm{\ω}}_e^{\left(n\right)}\left(k\right)=2\mathrm{\π}{f}_e^{\left(n\right)}\left(k\right)(n=\mathrm{0,1,2})$ 为载波相位θ_e(k)的n+1阶导数(分别表示残余多普勒频率f_e(k)、残余多普勒频率变化率f_e ⁽¹⁾(k)、残余多普勒频率二阶变化率f_e ⁽²⁾(k))，它们在采样间隔T内的迭代关系可以按照泰勒级数展开得： 

$\left\{\begin{array}{c}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_e^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_e^{\left(3\right)}\left({\mathrm{\ϵ}}_1\right)\·T^4}{24}\\ f_e(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_e^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_2\right)\·T^3}{6}\\ f_e^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_e^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_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({\mathrm{\ϵ}}_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=\mathrm{0,1,2,3})---\left(10\right)$

在此给出两种状态空间描述方式作为本发明提出并推荐的载波参数估计模型，后面进行仿真试验比较了两种估计模型的性能。 

定义第k步到第k+1步的载波NCO积分多普勒值(第k步到第k+1步的载波相位差分值)为： 

$\mathrm{\Δ}{\mathrm{\θ}}_u(k+1)={\mathrm{\θ}}_u(k+1)-{\mathrm{\θ}}_u\left(k\right)$

$={\∫}_{\mathrm{kT}}^{(k+1)T}2\mathrm{\π}{f}_u\left(t\right)\mathrm{dt}$

$=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_u^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_u^{\left(3\right)}\left({\mathrm{\ξ}}_1\right)\·T^4}{24}---\left(11\right)$

将公式(10)、公式(11)代入公式(8)得： 

$\left\{\begin{array}{c}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_d^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_d^{\left(3\right)}\left({\mathrm{\ξ}}_1\right)\·T^4}{24}-\mathrm{\Δ}{\mathrm{\θ}}_u(k+1)\\ f_d(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ξ}}_2\right)\·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ξ}}_3\right)\·T^2}{2}\\ f_d^{\left(2\right)}(k+1)=\underset{n=2}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-2}}{(n-2)!}+f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_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}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_d^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}-\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_u^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_e^{\left(3\right)}\left({\mathrm{\ϵ}}_1\right)\·T^4}{24}\\ f_d(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}-\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_2\right)\·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}-\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_n^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ϵ}}_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_d^{\left(3\right)}\left({\mathrm{\ϵ}}_4\right)\·T\end{array}\right.---\left(13\right)$

$\left\{\begin{array}{c}{f}_u(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_u^{\left(3\right)}\left({\mathrm{\ζ}}_1\right)\·T^3}{6}\\ f_u^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_u^{\left(3\right)}\left({\mathrm{\ζ}}_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({\mathrm{\ζ}}_3\right)\·T\end{array}\right.---\left(14\right)$

根据公式(14)给出的关于f_u(k+1)、f_u ⁽¹⁾(k+1)、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)\≈A\left(k\right)\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\·\cos {\mathrm{\θ}}_e\left(k\right)+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\≈A\left(k\right)\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[\mathrm{\Δ}{\mathrm{\ω}}_d\left(k\right)\·N/2]\·\sin {\mathrm{\θ}}_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 {\mathrm{\θ}}_e\left(k\right)\\ \sin {\mathrm{\θ}}_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为归一化量测噪声向量，一般认为是高斯白噪声，其协方差阵为 $R=E\left[n\left(k\right)n^T\left(k\right)\right]={\mathrm{\σ}}_n^{2}I,$ ${\mathrm{\σ}}_n^2={10}^{-\mathrm{SNR}}$ (SNR为积分-清除器输出信号的信噪比)；I为二阶单位阵。 

(iv)建立两种载波参数估计模型 

①估计模型1 

根据公式(11)、公式(12)、公式(17)建立估计模型1。状态向量取为x(k)＝[θ_e(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{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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\left(k\right)={[{\mathrm{\θ}}_e\left(k\right),f_d\left(k\right),{\stackrel{\·}{f}}_d\left(k\right),{\stackrel{\·\·}{f}}_d\left(k\right)]}^T,$ 控制输入量为 $u\left(k\right)={[{\hat{f}}_u\left(k\right),{\hat{f}}_u^{\′}\left(k\right),{\hat{f}}_u^{\′\′}\left(k\right)]}^T.$ 状态方程可表示为： 

x(k+1)＝Fx(k)+Gu(k)+w(k)         (19) 

公式(19)中：状态转移阵 $F=\left[\begin{array}{cccc}1& 2\mathrm{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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)定义。 

控制输入量为 $u\left(k\right)={[{\hat{f}}_u\left(k\right),{\hat{f}}_u^{\′}\left(k\right),{\hat{f}}_u^{\′\′}\left(k\right)]}^T$ 由公式(13)建立的线性估计模型通过线性kalman滤波获得。令状态向量为：x_u(k)＝[f_u(k)，f_u′(k)，f_u″(k)]^T，根据公式(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)为量测噪声，与频标的相位噪声(晶振的短期稳定度或铷原子频标Allan方差)和载波NCO位数(对于32bit以上的计数器，可以忽略)有关，是微小量。 

两种估计模型中，由公式(18)和公式(17)定义的估计模型1是含标量控制输入的4阶非线性统计估计模型；由公式(19)和公式(17)定义的估计模型2是含3阶控制输入向量的4阶非线性统计估计模型，3阶控制输入向量可由公式(20)和公式(11)定义的线性估计模型通过线性滤波实时计算。估计模型2是由两个子系统构成的分布式二级递推估计模型，总模型阶数和计算复杂度略高于估计模型1，但从后面的仿真实验可以看出，选择同样滤波算法的前提下，模型阶数的提高的确能够带来更高的精度，但计算量代价随之增加。 

载波参数的非线性估计模型一般采用非线性系统滤波算法如：EKF、UKF等实现状态估计。另外，也可尝试选择类似的非线性滤波算法如：粒子滤波、Stirling一阶/二阶插值滤波、二阶EKF、迭代EKF等等算法。由于本发明采用UKF作为两种估计模型的滤波算法，进行性能比较。 

(二)UKF非线性滤波算法简介 

EKF是常用于非线性动力学系统的统计滤波，应用广泛，这里不再赘述。EKF对于非线性的观测方程采用线性化近似方法，需要计算Jacobi矩阵，使得精度、稳定性受到影响、实现变得复杂。UKF是一种基于UT变换的线性滤波算法。UT变换是一种计算随机变量经过非线性变换后统计分布的方法，它采用确定性采样方法来捕获分布的高阶信息，能更好地逼近系统的非线性特征，从而较好地解决非线性问题。将UT变换应用到卡尔曼滤波的时间更新和测量更新中就得到了UKF算法。它仍然采用与EKF类似的一套递推公式，通过状态向量与估计协方差阵的递推以及利用测量时刻的观测信息来进行更新，得到状态估计值和估计协方差阵。与EKF不同的是，UKF利用一系列Sigma采样点，通过UT变换来进行状态与协方差阵的递推和更新，这样做使得估计结果能够以三阶精度逼近真实的状态值，远高于EKF[，而且使得稳定性也优于EKF。并且由于UKF不需要计算Jacobi矩阵，实现相对简单，计算量与EKF相当。 

另外由于频率估计的系统噪声与观测噪声均为加性噪声，因此可以只对状态进行Sigma点采样，而将系统噪声与观测噪声的信息提出来处理，从而使得算法得到简化，降低了计算量。UKF算法步骤如下： 

1)初始化 

设定初始状态估计值 与初始估计协方差阵P(0|0)。 

2)时间更新 

a.给定Sigma采样点与相应加权因子： 

$\mathrm{\χ}(k-1|k-1)=[\hat{X}(k-1|k-1),\hat{X}(k-1|k-1)+\sqrt{(m+\mathrm{\λ})P(k-1)(k-1)},\hat{X}(k-1|k-1)-\sqrt{(m+\mathrm{\λ})P(k-1|k-1)}]---\left(22\right)$

$W_0^{\left(m\right)}=\mathrm{\λ}/(m+\mathrm{\λ})$ $W_0^{\left(c\right)}=\mathrm{\λ}/(m+\mathrm{\λ})+(1-{\mathrm{\α}}^2+\mathrm{\β})$

$W_i^{\left(m\right)}=W_i^{\left(c\right)}=0.5/(m+\mathrm{\λ})(i=1,...,2m)---\left(23\right)$

λ＝α²(m+κ)-m，计算时κ＝0，β＝2，α＝0.01 

公式(23)中：m为状态的维数4，矩阵x(k-1|k-1)共有2m+1列，每列是一个Sigma采样点x_i(k-1|k-1)，(i＝0，…，2m)。 

b．状态预测 

χ(k|k-1)＝Φχ(k-1|k-1)， $\hat{X}\left(k\right|k-1)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i^{\left(m\right)}{\mathrm{\χ}}_i\left(k\right|k-1)---\left(24\right)$

$P\left(k\right|k-1)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}({\mathrm{\χ}}_i\left(k\right|k-1)-\hat{X}\left(k\right|k-1)){({\mathrm{\χ}}_i\left(k\right|k-1)-\hat{X}\left(k\right|k-1))}^T+Q---\left(25\right)$

Z_i(k|k-1)＝h(x_i(k|k-1))(i＝0，…，2m)， $\hat{Z}\left(k\right|k-1)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i^{\left(m\right)}{Z}_i\left(k\right|k-1)---\left(26\right)$

3)测量更新 

$P_{\tilde{Z}\left(k\right|k-1)}=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}(Z_i\left(k\right|k-1)-\hat{Z}\left(k\right|k-1)){(Z_i\left(k\right|k-1)-\hat{Z}\left(k\right|k-1))}^T+R---\left(27\right)$

$P_{\hat{X}\left(k\right|k-1)\hat{Z}\left(k\right|k-1)}=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}({\mathrm{\χ}}_i\left(k\right|k-1)-\hat{X}\left(k\right|k-1)){(Z_i\left(k\right|k-1)-\hat{Z}\left(k\right|k-1))}^T---\left(28\right)$

$K\left(k\right)=P_{\hat{X}\left(k\right|k-1)\hat{Z}\left(k\right|k-1)}{\left(P_{\hat{Z}\left(k\right|k-1)}\right)}^{-1}---\left(29\right)$

$\hat{X}\left(k\right|k)=\hat{X}\left(k\right|k-1)+K\left(k\right)(Z\left(k\right)-\hat{Z}\left(k\right|k-1)),$ $P\left(k\right|k)=P\left(k\right|k-1)-K\left(k\right)P_{\hat{Z}\left(k\right|k-1)}{K}^T\left(k\right)---\left(30\right)$

Claims (1)

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}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{(N+1)}\left(\mathrm{\ϵ}\right)\·T^{N+1}}{(N+1)!}(\mathrm{\ϵ}\∈[\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}{\mathrm{\Σ}}}\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)推导预报误差的解析形式：

$\mathrm{\δ}{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}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{(N+1)}\left(\mathrm{\ϵ}\right)\·T^{N+1}}{(N+1)!}-\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{{\hat{f}}_d^{\left(n\right)}\left(k\right)\·T^n}{n!}---\left(4\right)$

$=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\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(\mathrm{\ϵ}\right)\·T^{N+1}}{(N+1)!}$

$=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\frac{{\mathrm{\δf}}_d^{\left(n\right)}\left(k\right)}{n!}\·T^n+\frac{f_d^{(N+1)}\left(\mathrm{\ϵ}\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{\δf}}_d(k+1|k)-{\mathrm{\δ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)\≈A\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[{\mathrm{\Δ\ω}}_d\left(k\right)\·N/2]\cos {\mathrm{\θ}}_k+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\≈A\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[{\mathrm{\Δ\ω}}_d\left(k\right)\·N/2]\sin {\mathrm{\θ}}_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}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{2\mathrm{\π}{f}_e^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{{2\mathrm{\πf}}_e^{\left(3\right)}\left({\mathrm{\ϵ}}_1\right)\·T^4}{24}\\ f_e(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_e^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_2\right)\·T^3}{6}\\ f_e^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_e^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_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({\mathrm{\ϵ}}_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{\Δ\θ}}_u(k+1)={\mathrm{\θ}}_u(k+1)-{\mathrm{\θ}}_u\left(k\right)$

$={\∫}_{\mathrm{kT}}^{(k+1)T}2\mathrm{\π}{f}_u\left(t\right)\mathrm{dt}$ $\left(11\right)$

$=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{{2\mathrm{\πf}}_u^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{{2\mathrm{\πf}}_u^{\left(3\right)}\left({\mathrm{\ξ}}_1\right)\·T^4}{24}$

将公式(10)、公式(11)代入公式(8)得：

$\left\{\begin{array}{c}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{{2\mathrm{\πf}}_d^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{{2\mathrm{\πf}}_d^{\left(3\right)}\left({\mathrm{\ξ}}_1\right)\·T^4}{24}-{\mathrm{\Δ\θ}}_u(k+1)\\ f_d(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ξ}}_2\right)\·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_d^{\left(3\right)}\left({\mathrm{\ξ}}_3\right)\·T^2}{2}\\ f_d^{\left(2\right)}(k+1)=\underset{n=2}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-2}}{(n-2)!}+f_d^{\left(3\right)}\left({\mathrm{\ξ}}_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}{\mathrm{\θ}}_e(k+1)={\mathrm{\θ}}_e\left(k\right)+\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{{2\mathrm{\πf}}_d^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}-\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{{2\mathrm{\πf}}_u^{\left(n\right)}\left(k\right)\·T^{n+1}}{(n+1)!}+\frac{2\mathrm{\π}{f}_e^{\left(3\right)}\left({\mathrm{\ϵ}}_1\right)\·T^4}{24}\\ f_d(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^n}{n!}-\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_2\right)\·T^3}{6}\\ f_d^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_d^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}-\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_e^{\left(3\right)}\left({\mathrm{\ϵ}}_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({\mathrm{\ϵ}}_4\right)\·T\end{array}\right.---\left(13\right)$

$\left\{\begin{array}{c}{f}_u(k+1)=\underset{n=0}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^n}{n!}+\frac{f_u^{\left(3\right)}\left({\mathrm{\ζ}}_1\right)\·T^3}{6}\\ f_u^{\left(1\right)}(k+1)=\underset{n=1}{\overset{2}{\mathrm{\Σ}}}\frac{f_u^{\left(n\right)}\left(k\right)\·T^{n-1}}{(n-1)!}+\frac{f_u^{\left(3\right)}\left({\mathrm{\ζ}}_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({\mathrm{\ζ}}_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)\≈A\left(k\right)\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[{\mathrm{\Δ\ω}}_d\left(k\right)\·N/2]\·\cos {\mathrm{\θ}}_e\left(k\right)+n_I\left(k\right)\\ Q_{\mathrm{ps}}\left(k\right)\≈A\left(k\right)\·R\left[\mathrm{\ϵ}\left(k\right)\right]\·\sin c[{\mathrm{\Δ\ω}}_d\left(k\right)\·N/2]\·\sin {\mathrm{\θ}}_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 {\mathrm{\θ}}_e\left(k\right)\\ \sin {\mathrm{\θ}}_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{\πT}& \mathrm{\π}{T}^2& {\mathrm{\π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{\πT}& \mathrm{\π}{T}^2& {\mathrm{\π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{\πT}& \mathrm{\π}{T}^2& \mathrm{\π}{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)定义的线性估计模型通过线性滤波实时计算。

Images