CN102540216A - 一种自适应跟踪环路及实现方法 - Google Patents

Abstract

本发明公开了一种自适应跟踪环路，包括UKF滤波器、观测噪声方差阵检测模块、模糊推理系统、UT变换比例因子调节模块、状态补偿器、载波NCO、比例因子、码NCO、积分和清零模块、码环鉴相器和二阶码环滤波器；还公开了一种自适应跟踪环路的实现方法，包括步骤一：信号相关及积分、清零处理；步骤二：码相位跟踪；步骤三：UKF滤波器建模；步骤四：估计观测噪声方差阵；步骤五：估计过程噪声方差阵；步骤六：调节UT变换比例因子；步骤七：补偿状态估计偏差；步骤八：载波NCO辅助码NCO。本发明在载波跟踪环路中设计了UKF滤波器、观测噪声方差阵检测模块以及模糊推理系统，不仅能够解决跟踪环路中热噪声颤动与动态应力误差之间的矛盾，而且能够根据外界环境的变化对过程噪声方差阵与观测噪声方差阵进行自适应调节，有效提高了跟踪环路在高动态、强干扰等复杂多变环境下的自适应能力。

Description

一种自适应跟踪环路及实现方法

技术领域

本发明属于导航技术领域，涉及一种自适应跟踪环路及实现方法。

背景技术

GPS系统能够全天候地在全球范围内提供精确的PVT信息(位置、速度、时间)，在军用以及民用领域具有十分广阔的用途和发展前景。GPS接收机不仅可用于弹道导弹、巡航导航、超视距空空导弹、制导炸弹以及炮弹等武器的精确制导，而且在载人航天、再入飞行器、太空探索等应用领域也发挥了重要作用。

在上述应用领域，GPS接收机通常工作在高动态、强干扰等复杂多变的环境下，这将给接收机跟踪GPS信号带来以下问题：(1)载体的高动态特性使GPS载波产生较大的多普勒频移，对于普通的载波跟踪环路，较大的载波多普勒频移会超出跟踪环路的带宽而造成跟踪环路失锁。(2)外界的噪声干扰将导致GPS信号载噪比降低，进而加剧环路的热噪声颤动，引起测量误差增加，当测量误差超过跟踪门限时将造成跟踪环路失锁而无法正常工作。

为了改善GPS接收机跟踪环路在高动态、强干扰环境下的跟踪性能，国内外相关学者围绕跟踪环路开展了大量的研究工作。W.J.Hurd等利用准开环、近似最大似然估计方法(MLE)，通过估计多普勒频率和码延迟达到对伪距及其变化率的估计的目的，此方案增强了接收机在高动态下保持跟踪的能力。Miao Jianfeng等提出了一种基于改进的Kalman滤波的跟踪算法。这种算法通过对动态模型噪声和测量噪声进行解相关，可实现在强干扰环境下更稳定、更精确地跟踪载波相位和多普勒频移。Matthew Lashley与David M.Bevly研究了基于向量的跟踪算法，它不仅利用扩展卡尔曼滤波器来跟踪GPS信号，同时对接收机的位置、速度和时钟偏差等信息进行估计，适用于在高动态环境下对弱信号进行精确跟踪。Wei-Lung Mao等对基于模糊控制理论的跟踪环路进行了研究，它利用模糊带宽控制器和三阶锁相环对载波相位进行跟踪，可在高动态环境下实现更快速的跟踪和更宽的频率跟踪范围。Zhu Yunlong等提出了一种基于改进的无迹卡尔曼滤波器(UKF)的跟踪算法，该算法根据线性状态方程对UKF滤波器进行简化运算，不仅显著降低了UKF滤波器的计算复杂度，而且在高动态环境下可有效地改善跟踪环路的载噪比失锁门限，提高了跟踪环路的抗干扰性能。

然而，当GPS接收机工作在复杂多变的高动态、强干扰环境下时，上述所设计的跟踪环路无法随着外界环境的变化而自适应地对跟踪环路参数进行调节，这时将导致跟踪环路的跟踪精度降低，甚至会造成跟踪环路失锁而无法正常工作。

发明内容

本发明的目的是为了解决现有GPS接收机跟踪环路在高动态、强干扰等复杂多变环境下无法正常跟踪卫星信号的问题，提供了一种自适应跟踪环路及实现方法，利用改进的UKF滤波器来替代传统载波环路中的鉴相器与低通滤波器，直接估计载波相位与载波频率并能够有效滤除噪声信号，不仅克服了鉴相器工作范围小的缺点，而且消除了热噪声颤动与动态应力误差之间的矛盾，有效地改善了跟踪环路在高动态、强干扰等复杂多变环境下的跟踪性能。

一种自适应跟踪环路的实现方法，具体步骤如下：

步骤一：信号相关及积分、清零处理；

步骤二：码相位跟踪；

步骤三：UKF滤波器建模；

步骤四：估计观测噪声方差阵；

步骤五：估计过程噪声方差阵；

步骤六：调节UT变换比例因子；

步骤七：补偿状态估计偏差；

步骤八：载波NCO辅助码NCO。

一种自适应跟踪环路，包括UKF滤波器、观测噪声方差阵检测模块、模糊推理系统、UT变换(无迹变换)比例因子调节模块、状态补偿器、载波NCO(数控振荡器)、比例因子、码NCO、积分和清零模块、码环鉴相器和二阶码环滤波器；

载波NCO生成的本地同相载波cos信号和本地正交载波sin信号，输入的数字中频信号分别与本地同相载波cos信号、本地正交载波sin信号进行相关，得到同相支路数据A、正交支路数据A；码NCO生成的本地C/A码，分别同相支路数据A、正交支路数据A进行相关，得到同相支路数据B和正交支路数据B；积分和清零模块对同相支路数据B、正交支路数据B进行积分和累加，得到同相支路相关能量与正交支路相关能量，同相支路相关能量输出至码环鉴相器、观测噪声方差阵检测模块和UKF滤波器，正交支路相关能量也输出至码环鉴相器、观测噪声方差阵检测模块和UKF滤波器；码环鉴相器对同相支路相关能量和正交支路相关能量进行鉴相处理，得到数字中频信号中的C/A码与码NCO所生成本地C/A码的码相位误差；二阶码环滤波器滤除码相位误差中的噪声信号，然后将码相位误差输出至码NCO，码NCO根据码相位误差调整本地C/A码，使本地C/A码相位与输入数字中频信号中的C/A码相位对齐；载波NCO将输入数字中频信号中载波多普勒频移输出至比例因子，比例因子进行比例运算，得到输入数字中频信号中码多普勒频移，输出至码NCO，调整码NCO内的本地C/A码频率，使本地C/A码频率与输入数字中频信号中C/A码频率对齐；

UKF滤波器包括初始化环节、Sigma点采样环节、时间更新环节和量测更新环节；初始化环节将初始状态估计值

与初始状态协方差阵P₀均设置为零，并赋值给初始观测噪声方差阵R₀和初始过程噪声方差阵Q₀，然后将初始值传输给Sigma点采样环节；Sigma点采样环节在跟踪环路首个循环周期内，根据初始状态估计值及初始状态协方差P₀、UT变换比例因子α选取2n+1个Sigma采样点，n为状态向量的维数；如果不在首个循环周期，则Sigma点采样环节根据量测更新环节提供的状态估计值

状态协方差P与UT变换比例因子α来选取Sigma采样点，并使这些采样点的均值和协方差等于状态向量的均值和协方差，然后将Sigma采样点输出给时间更新环节；时间更新环节将Sigma采样点导入到状态方程和观测方程，结合解模糊环节输出的Q阵(过程噪声方差阵)，得到观测向量预测估计值、状态预测估计值及状态协方差预测估计值，并输出给量测更新环节；量测更新环节根据同相支路相关能量、正交支路相关能量，并结合观测噪声方差阵计算环节输出的R阵(观测噪声方差阵)，对状态预测估计值及状态协方差预测估计值进行滤波更新，得到状态估计值

状态协方差P、残差e、输入数字中频信号中载波与本地同相载波相位差的三阶导数ω₂和残差协方差理论值

将状态估计值

状态协方差P反馈给Sigma点采样环节，并将状态协方差P输出至UT变换比例因子调节模块，残差e输出至残差协方差实际值检测环节，输入数字中频信号中载波与本地同相载波相位差的三阶导数ω₂输出至加加速度变化率检测环节，残差协方差理论值

输出至残差比较环节；

观测噪声方差阵检测模块包括载波功率检测环节、载噪比检测环节和观测噪声方差阵计算环节；载波功率检测环节检测同相支路相关能量，获得数字中频信号的载波功率C，将载波功率C传输给观测噪声方差阵计算环节；同时，载噪比检测环节对同相支路相关能量、正交支路相关能量分别进行相干积分和非相干积分，得到窄带功率和宽带功率，并将窄带功率与宽带功率作比例运算后得到载噪比C/N₀，然后传输给观测噪声方差阵计算环节；观测噪声方差阵计算环节利用获得的载波功率C和载噪比C/N₀，并结合设置的预检测积分时间估计出观测噪声方差阵，进而输出R阵对量测更新环节中的观测噪声方差阵进行更新；

模糊推理系统包括残差协方差实际值检测环节、残差比较环节、加加速度变化率检测环节、模糊化环节、模糊规则和解模糊环节；残差协方差实际值检测环节通过检测量测更新环节中的残差e得到残差协方差实际值C_e，输出至残差比较环节，残差比较环节将C_e与量测更新环节中的残差协方差理论值

进行比较，将二者的比值r′输出至模糊化环节，同时，加加速度变化率检测环节根据量测更新环节中的状态估计值

实时提取数字中频信号中载波与本地同相载波相位差的三阶导数ω₂，进而得到载体的加加速度变化率δ′，输出至模糊化环节，模糊化环节将输入量r′、δ′转化为模糊量r、δ，并输出给模糊规则，模糊规则基于协方差匹配原理，将模糊量r与δ推理生成模糊量μ，解模糊环节利用重心法把模糊规则输出的模糊量μ转化为精确的输出变量μ′，得到精确的输出量μ′后，通过式Q＝10^μ′Q₀来对Q阵进行估计，并输出给时间更新环节，完成对Q阵的更新；

UT变换比例因子调节模块通过量测更新环节中的状态协方差P得到第k步状态估计值与状态真实值之间的距离以及第k+1步未作比例变换的Sigma点到中心点的距离，并将第k步状态真实值的邻域作为第k+1步Sigma点的采样范围，获得UT比例变换因子α，进而输出给Sigma点采样环节；

状态补偿器用来修正UKF滤波器的滞后性所造成的状态估计偏差，它将量测更新环节所得到的前一时刻的状态估计值

作为输入量，根据状态估计值

中的载波相位差θ及其一阶导数ω₀、二阶导数ω₁和三阶导数ω₂，利用泰勒级数展开式来预测当前时刻的载波信息U，将其输出给载波NCO，载波NCO根据状态补偿器输出的当前时刻载波信息U对本地载波频率与相位进行调节，使之与数字中频信号的载波频率、相位精确对齐；

跟踪环路在完成C/A码与载波的精确同步后，由同相支路相关能量解调出接收信号中的导航电文。

本发明的优点在于：

(1)本发明在载波跟踪环路中设计了具有非线性跟踪能力的UKF滤波器，能够解决高动态、强干扰环境下跟踪环路中热噪声颤动与动态应力误差之间的矛盾，实现在高动态、强干扰环境下精确地跟踪GPS信号；

(2)通过检测载波功率与载噪比以获得噪声功率谱密度，可实现对观测噪声方差阵的准确估计，从而使观测噪声方差阵能够随着外界干扰噪声的变化而进行自适应调节，有效提高了跟踪环路在时变噪声干扰环境下的适应能力；

(3)基于协方差匹配的原则，采用模糊推理系统使残差协方差实际值与其理论值保持一致，可实现对过程噪声方差阵的自适应调节，增强了跟踪环路在高动态环境下对载体运动状态的跟踪能力。

附图说明

图1为本发明一种自适应跟踪环路的结构示意图；

图2为本发明的UKF滤波器的具体工作原理图；

图3a为本发明模糊推理系统中输入变量r的隶属度函数；

图3b为本发明模糊推理系统中输入变量δ的隶属度函数；

图3c为本发明模糊推理系统中输出变量μ的隶属度函数；

图4为本发明的模糊推理系统中生成的模糊控制量；

图5为本发明一种自适应跟踪环路的实现方法的流程图；

图中：

1-UKF滤波器        2-观测噪声方差阵检测模块   3-模糊推理系统

4-UT变换比例因子调节模块 5-状态补偿器       6-载波NCO

7-比例因子               8-码NCO            9-积分和清零模块

10-码环鉴相器            11-二阶码环滤波器  101-初始化环节

102-Sigma点采样环节      103-时间更新环节   104-量测更新环节

201-载波功率检测环节     202-载噪比检测环节 203-观测噪声方差阵计算环节

301-残差协方差实际值检测 302-残差比较环节   303-加加速度变化率检测环节环节

304-模糊化环节           305-模糊规则       306-解模糊环节

具体实施方式

下面结合附图对本发明进行详细说明。

本发明的一种适用于高动态、强干扰环境的新型自适应跟踪环路，如图1所示，包括UKF滤波器1、观测噪声方差阵检测模块2、模糊推理系统3、UT变换比例因子调节模块4、状态补偿器5、载波NCO(数控振荡器)6、比例因子7、码NCO8、积分和清零模块9、码环鉴相器10和二阶码环滤波器11。

本发明在传统跟踪环路的基础上，利用UKF滤波器1来替代传统载波环路中的鉴相器与低通滤波器，并设计了观测噪声方差阵检测模块2、模糊推理系统3、UT变换比例因子调节模块4和状态补偿器5；同时，载波NCO6将调整后的本地载波频率通过比例因子7来对码NCO8进行频率辅助；

如图1所示，输入的数字中频信号分别与载波NCO6生成本地同相载波cos信号和本地正交载波sin信号进行相关，得到同相支路数据A和正交支路数据A。然后，同相支路数据A和正交支路数据A分别与码NCO8生成的本地C/A码进行相关，得到同相支路数据B和正交支路数据B。

同相支路数据B和正交支路数据B分别输入至积分和清零模块9，经过积分和累加，得到同相支路相关能量和正交支路相关能量，同相支路相关能量输出至码环鉴相器10、观测噪声方差阵检测模块2和UKF滤波器1，正交支路相关能量也输出至码环鉴相器10、观测噪声方差阵检测模块2和UKF滤波器1。

码环鉴相器10对同相支路相关能量和正交支路相关能量进行鉴相处理，得到数字中频信号C/A码与码NCO8所生成本地C/A码的码相位误差。

二阶码环滤波器11滤除码相位误差中的噪声信号，然后将码相位误差输出至码NCO8，码NCO8根据码相位误差调整本地C/A码，使本地C/A码相位与输入的数字中频信号中C/A码相位对齐。

此外，载波NCO6将输入数字中频信号中载波多普勒频移输出至比例因子7，比例因子7进行比例运算，得到输入数字中频信号中码多普勒频移，输出至码NCO8，调整码NCO8内的C/A码频率，使本地C/A码频率与输入数字中频信号中C/A码频率对齐。

如图2所示，UKF滤波器1包括初始化环节101、Sigma点采样环节102、时间更新环节103和量测更新环节104；

初始化环节101将初始状态估计值

与初始状态协方差阵P₀均设置为零，并赋值给初始观测噪声方差阵R₀和初始过程噪声方差阵Q₀，使其能够满足跟踪环路在信噪比为-20dB、加加速度为5g的正常环境下完成对卫星信号的跟踪功能，然后将初始值传输给Sigma点采样环节102。

如果在跟踪环路首个循环周期内，Sigma点采样环节102根据初始状态估计值

及初始状态协方差P₀、UT变换比例因子α选取2n+1个Sigma采样点(n为状态向量的维数)；如果不在首个循环周期，Sigma点采样环节102则根据量测更新环节104提供的状态估计值

状态协方差P与UT变换比例因子α来选取Sigma采样点，并使这些采样点的均值和协方差等于状态向量的均值和协方差，然后将Sigma采样点输出给时间更新环节103。

时间更新环节103将Sigma采样点导入到状态方程和观测方程，结合解模糊环节306输出的Q阵，得到观测向量预测估计值、状态预测估计值及状态协方差预测估计值，并输出给量测更新环节104。

量测更新环节104根据同相支路相关能量、正交支路相关能量，并结合观测噪声方差阵计算环节203输出的R阵，对状态预测估计值及状态协方差预测估计值进行滤波更新，得到状态估计值

状态协方差P、残差e、输入数字中频信号中载波与本地同相载波相位差的三阶导数ω₂和残差协方差理论值

将状态估计值

状态协方差P反馈给Sigma点采样环节102，并将状态协方差P输出至UT变换比例因子调节模块4，残差e输出至残差协方差实际值检测环节301，输入数字中频信号中载波与本地同相载波相位差的三阶导数ω₂输出至加加速度变化率检测环节303，残差协方差理论值

输出至残差比较环节302。

观测噪声方差阵检测模块2包括载波功率检测环节201、载噪比检测环节202和观测噪声方差阵计算环节203；

载波功率检测环节201检测同相支路相关能量，获得数字中频信号的载波功率C，将载波功率C传输给观测噪声方差阵计算环节203；同时，载噪比检测环节202对同相支路相关能量、正交支路相关能量分别进行相干积分和非相干积分，得到窄带功率和宽带功率，并将窄带功率与宽带功率作比例运算后得到载噪比C/N₀，然后传输给观测噪声方差阵计算环节203。

观测噪声方差阵计算环节203利用获得的载波功率C和载噪比C/N₀，并结合设置的预检测积分时间估计出观测噪声方差阵(R阵)，进而输出R阵对量测更新环节104中的观测噪声方差阵进行更新；

模糊推理系统3包括残差协方差实际值检测环节301、残差比较环节302、加加速度变化率检测环节303、模糊化环节304、模糊规则305和解模糊环节306；

残差协方差实际值检测环节301通过检测量测更新环节104中的残差e得到残差协方差实际值C_e，输出至残差比较环节302，残差比较环节302将C_e与量测更新环节104中的残差协方差理论值

进行比较，将二者的比值r′输出至模糊化环节304，同时，加加速度变化率检测环节303根据量测更新环节104中的状态估计值实时提取数字中频信号中载波与本地同相载波相位差的三阶导数ω₂，进而得到载体的加加速度变化率δ′，输出至模糊化环节304，模糊化环节304将输入量r′、δ′转化为模糊量r、δ，并输出给模糊规则305，模糊规则305基于协方差匹配原理，将模糊量r与δ推理生成模糊量μ，输入变量r、输入变量δ与输出变量μ的隶属度函数分别如图3a、图3b、图3c所示，模糊规则305中的模糊控制规则如表1所示：

表1控制规则

解模糊环节306利用重心法把模糊规则305输出的模糊量μ转化为精确的输出变量μ′，图4为采用重心法离线产生的控制量表，得到精确的输出量μ′后，通过式Q＝10^μ′Q₀来对Q阵进行估计，并输出给时间更新环节103，完成对Q阵的更新。

UT变换比例因子调节模块4通过量测更新环节104中的状态协方差P得到第k步状态估计值与状态真实值之间的距离以及第k+1步未作比例变换的Sigma点到中心点的距离，并将第k步状态真实值的邻域作为第k+1步Sigma点的采样范围，可获得UT比例变换因子α，进而输出给Sigma点采样环节102。

状态补偿器5用来修正UKF滤波器1的滞后性所造成的状态估计偏差，它将量测更新环节104所得到的前一时刻的状态估计值

作为输入量，根据状态估计值

中的载波相位差θ及其一阶导数ω₀、二阶导数ω₁和三阶导数ω₂，利用泰勒级数展开式来预测当前时刻的载波信息U(包括载波相位与角频率)，将其输出给载波NCO6，载波NCO6根据状态补偿器5输出的当前时刻载波信息U对本地载波频率与相位进行调节，使之与数字中频信号的载波频率、相位精确对齐，从而有效地补偿滤波器滞后性所造成的跟踪误差；

跟踪环路在完成C/A码与载波的精确同步后，由同相支路相关能量能够解调出接收信号中的导航电文，从而可为后面的导航解算模块完成定位功能提供GPS卫星星历数据，解算出GPS接收机的三维位置、速度和时间。

本发明的一种自适应跟踪环路的实现方法，流程如图5所示，具体包括以下步骤：

步骤一：信号相关及积分、清零处理；

输入的数字中频信号分别与本地同相载波和本地正交载波相关后进入C/A码跟踪环路，经过本地C/A码解扩以及积分和清零模块9后，得到同相支路相关能量I_P(i)和正交支路相关能量Q_P(i)，它们表示为：

$\left.\begin{array}{c}{I}_P\left(i\right)=0.5\·K_C\·A\·D_i\·R\left[{\mathrm{\ϵ}}_i\right]\·\frac{\sin (\mathrm{\π}{T}_I\·\mathrm{\Δf})}{\mathrm{\π}{T}_I\·\mathrm{\Δf}}\·\cos \left(\mathrm{\θ}\right)+N_I\left(i\right)\\ Q_P\left(i\right)=0.5\·K_C\·A\·D_i\·R\left[{\mathrm{\ϵ}}_i\right]\·\frac{\sin (\mathrm{\π}{T}_I\·\mathrm{\Δf})}{\mathrm{\π}{T}_I\·\mathrm{\Δf}}\·\sin \left(\mathrm{\θ}\right)+N_Q\left(i\right)\end{array}\right\}---\left(1\right)$

式中，K_C为一个预检测积分时间内累加的次数，K_C＝T_I/T_S，T_I为预检测积分时间，T_S为数字中频信号的采样时间间隔，A为输入中频信号的幅值，D_i为导航电文，R[ε_i]是归一化的C/A码自相关函数，ε_i为本地C/A码与输入C/A码的相位差，

表示本地同相载波与输入信号中载波的相位差，Δf、

分别为积分间隔起始时刻本地参考信号与输入信号之间的载波频率差和相位差，N_I(i)、N_Q(i)分别为同相支路高斯白噪声和正交支路高斯白噪声。

步骤二：码相位跟踪；

码环鉴相器10根据积分和清零模块9输出的数据得到码相位误差，将得到的码相位误差输出至二阶码环滤波器11，码相位误差经过二阶码环滤波器11滤波处理后，对码NCO8输出控制信号，码NCO8根据控制信号和比例因子7输出的辅助参数调整本地C/A码相位和频率，使得本地C/A码与数字中频信号中的C/A码精确对齐。

步骤三：UKF滤波器1建模；

将积分和清零模块9输出的I_P(i)和Q_P(i)作为UKF滤波器1的输入观测向量Z_k，对UKF滤波器1中的状态估计值进行更新，下面对UKF滤波器1进行精确建模。

(1)状态方程与观测方程

UKF滤波器1采用四维状态向量对载波信号进行跟踪，其状态向量X分别为本地同相载波与输入载波相位差θ及其一阶导数ω₀、二阶导数ω₁和三阶导数ω₂，将状态向量根据泰勒级数展开公式可得系统状态方程为：

X_k＝Φ_k，k-1X_k-1+G_k，k-1U_k-1+W_k     (2)

即： $X_k={\left[\begin{array}{c}\mathrm{\θ}\\ {\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ {\mathrm{\ω}}_2\end{array}\right]}_k=\left[\begin{array}{cccc}1& T_I& T_I^2/2& T_I^2/6\\ 0& 1& T_S& T_I^2/2\\ 0& 0& 1& T_I\\ 0& 0& 0& 1\end{array}\right]{\left[\begin{array}{c}\mathrm{\θ}\\ {\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ {\mathrm{\ω}}_2\end{array}\right]}_{k-1}+\left[\begin{array}{cc}-1& 0\\ 0& -1\\ 0& 0\\ 0& 0\end{array}\right]{\left[\begin{array}{c}\mathrm{\Δ\θ}\\ \mathrm{\Δ}{\mathrm{\ω}}_0\end{array}\right]}_{k-1}+{\left[\begin{array}{c}{\mathrm{\η}}_1\\ {\mathrm{\η}}_2\\ {\mathrm{\η}}_3\\ {\mathrm{\η}}_4\end{array}\right]}_k$

式中，X_k、X_k-1分别为第k步状态向量和第k-1步状态向量，Φ_k，k-1为系统的状态转移矩阵，T_I为预检测积分时间，G_k，k-1为输入转移矩阵，U_k-1为控制输入量，Δθ、Δω₀分别为载波NCO6相位调整量和角频率调整量，W_k为动态模型噪声，

i＝1，…，4，为泰勒级数展开式的余项，用来描述上述模型受到某些随机干扰以及模型不准确所造成的影响，Y(t)表示连续相位过程的四阶变化率，被视为单边谱密度为N_y的高斯白噪声。

进而由系统的观测向量Z_k与状态向量可得到系统的观测方程：

$Z_k={\left[\begin{array}{c}{I}_P\left(i\right)\\ Q_P\left(i\right)\end{array}\right]}_k=\left[\begin{array}{c}{K}^{\′}\cos (L\·X_k)\\ K^{\′}\sin (L\·X_k)\end{array}\right]+V_k---\left(3\right)$

式中，

L＝[1，0，0，0]，V_k＝[N_I(i)，N_Q(i)]^T，为观测噪声向量，其方差阵为：

$R_k=E\left[V_k{V}_k^T\right]=\frac{N_0}{2T_I}\·I_{2\×2}---\left(4\right)$

式中，N₀为噪声功率谱密度，I_2×2为二阶单位阵。

(2)初始化UKF滤波器1

在UKF滤波器1工作的初始阶段，初始化环节101将初始状态估计值

与初始状态协方差阵P₀均设置为零，并赋值给初始观测噪声方差阵R₀和初始过程噪声方差阵Q₀，使其能够满足跟踪环路在信噪比为-20dB、加加速度为5g的正常环境下完成对卫星信号的跟踪功能，然后将初始值传输给Sigma点采样环节102。

(3)Sigma点采样

Sigma点采样环节102根据状态估计值及状态协方差选取2n+1个Sigma采样点，Sigma采样点用矩阵χ_k-1来表示，χ_k-1可表示如下：

$\left.\begin{array}{c}{\mathrm{\χ}}_{i,k-1}={\hat{X}}_{k-1},i=0\\ {\mathrm{\χ}}_{i,k-1}={\hat{X}}_{k-1}+{\left(\sqrt{(n+\mathrm{\λ})P_{k-1}}\right)}_i,i=1,...,n\\ {\mathrm{\χ}}_{i,k-1}={\hat{X}}_{k-1}-{\left(\sqrt{(n+\mathrm{\λ})P_{k-1}}\right)}_{i-n},i=n+1,...,2n\end{array}\right\}---\left(5\right)$

式中，χ_i，k-1为x_k-1的第i列，

表示第k-1步状态估计值，n为状态向量的维数，λ＝α²(n+τ)-n是一个尺度参数，α为UT变换(无迹变换)比例因子，它决定了采样点在其均值的分布情况，τ＝3-n为次级尺度调节因子，P_k-1为第k-1步状态协方差，

表示矩阵(n+λ)P_k-1平方根的第i列。

Sigma采样点的一阶统计特性权系数

和二阶统计权系数

为：

$\left.\begin{array}{c}{W}_i^{\left(m\right)}=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},i=0\\ W_i^{\left(c\right)}=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β}),i=0\\ W_i^{\left(m\right)}=W_i^{\left(c\right)}=\frac{1}{2(n+\mathrm{\λ})},i=1,...,2n\end{array}\right\}---\left(6\right)$

式中，

为用来计算均值的权值，为用来计算协方差的权值，β为验前分布因子(对于高斯分布，β的最优值为2)。

(4)时间更新

时间更新环节103根据Sigma点采样环节102输出的采样点和解模糊环节306输出的Q阵，得到观测向量预测估计值

状态预测估计值

及状态协方差预测估计值P_k，k-1。

${\mathrm{\χ}}_{k,k-1}^{*}={\mathrm{\Φ}}_{k,k-1}{\mathrm{\χ}}_{k-1}+G_{k,k-1}{U}_{k-1}---\left(7\right)$

${\hat{X}}_{k,k-1}=\underset{i=0}{\overset{i=2n}{\mathrm{\Σ}}}{W}_i^{\left(m\right)}{\mathrm{\χ}}_{i,k,k-1}^{*}---\left(8\right)$

$P_{k,k-1}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[{\mathrm{\χ}}_{i,k,k-1}^{*}-{\hat{X}}_{k,k-1}]{[{\mathrm{\χ}}_{i,k,k-1}^{*}-{\hat{X}}_{k,k-1}]}^T+Q---\left(9\right)$

$Z_{k,k-1}^{*}=\left[\begin{array}{c}{K}^{\′}\cos (L\·{\mathrm{\χ}}_{k,k-1}^{*})\\ K^{\′}\sin (L\·{\mathrm{\χ}}_{k,k-1}^{*})\end{array}\right]---\left(10\right)$

${\hat{Z}}_{k,k-1}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}{Z}_{i,k,k-1}^{*}---\left(11\right)$

式中，

为χ_k-1经状态方程变换后得到的点集，

为

经观测方程变换后得到的预测点集，

与

分别为矩阵与矩阵

的第i列。

(5)量测更新

量测更新环节104根据时间更新环节103的输出和观测向量Z_k，并结合观测噪声方差阵计算环节203输出的R阵，获得残差协方差理论值状态估计值

残差e以及状态协方差P。

$P_{\tilde{Z}\tilde{Z}}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[Z_{i,k,k-1}^{*}-{\hat{Z}}_{k,k-1}]{[Z_{i,k,k-1}^{*}-{\hat{Z}}_{k,k-1}]}^T+R---\left(12\right)$

$P_{\tilde{X}\tilde{Z}}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[{\mathrm{\χ}}_{i,k,k-1}^{*}-{\hat{X}}_{k,k-1}][Z_{i,k,k-1}^{*}-{\hat{Z}}_{k,k-1}]---\left(13\right)$

$K_k=P_{\tilde{X}\tilde{Z}}{P}_{\tilde{Z}\tilde{Z}}^{-1}---\left(14\right)$

$e_k=Z_k-{\hat{Z}}_{k,k-1}---\left(15\right)$

${\hat{X}}_k={\hat{X}}_{k,k-1}+K_k{e}_k---\left(16\right)$

$P_k=P_{k,k-1}-K_k{P}_{\tilde{Z}\tilde{Z}}{K}_k^T---\left(17\right)$

式中，

为观测向量预测估计值与状态预测估计值的互协方差，K_k为UKF滤波器1的增益，e_k为第k步的残差，

为第k步的状态估计值，P_k为第k步的状态协方差。

步骤四：估计观测噪声方差阵；

载波功率检测环节201通过对式(1)中的同相支路相关能量I_P(i)进行计算可直接获得输入中频载波信号的功率。在式(1)中，为了消除噪声N_I(i)和导航电文翻转对计算载波功率C的影响，这里用K_A个I_P(i)采样点的绝对值来求取载波功率：

$C=\frac{1}{2}{A}^2=\frac{1}{2}\·{\left(\frac{\underset{i=1}{\overset{K_A}{\mathrm{\Σ}}}\left|I_P\left(i\right)\right|}{0.5K_C\·K_A}\right)}^2---\left(18\right)$

同时，载噪比检测环节202在M次预检测积分时间内对同相支路相关能量I_P(i)与正交支路相关能量Q_P(i)分别进行相干积分和非相干积分，得到窄带功率NBP和宽带功率WBP：

$\mathrm{NBP}={\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{I}_P\left(i\right)\right)}_j^2+{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{Q}_P\left(i\right)\right)}^2---\left(19\right)$

$\mathrm{WBP}=\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}(I_P^2\left(i\right)+Q_P^2\left(i\right))\right)---\left(20\right)$

式中，M表示噪声带宽所改变的倍数，WBP的带宽为1/T_I，NBP的带宽为1/(MT_I)。

则宽带功率WBP和窄带功率NBP的数学期望分别为：

E(NBP)＝8M²(C+N₀/(MT_I))       (21)

E(WBP)＝8M(C+N₀/T_I)           (22)

为了提高载噪比的检测精度，用N个NBP和WBP采样点的平均值来分别表示E(NBP)和E(WBP)，则E(NBP)与E(WBP)的比值E(NP)可表示为：

$E\left(\mathrm{NP}\right)=\frac{E\left(\mathrm{NBP}\right)}{E\left(\mathrm{WBP}\right)}=\frac{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{NBP}}_j}{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{WBP}}_j}=\frac{{\mathrm{MT}}_I\·C/N_0+1}{T_I\·C/N_0+1}---\left(23\right)$

对式(23)进行运算便可得到信号载噪比C/N₀：

$C/N_0=\frac{E\left(\mathrm{NP}\right)-1}{M-E\left(\mathrm{NP}\right)}\·\frac{1}{T_I}---\left(24\right)$

综合式(18)和式(24)可知噪声功率谱密度为N₀：

$N_0=\frac{T_I}{2}\·\frac{{\left(\underset{i=1}{\overset{K_A}{\mathrm{\Σ}}}\right|I_P\left(i\right)\left|\right)}^2\·(M-E\left(\mathrm{NP}\right))}{K_A^2\·(E\left(\mathrm{NP}\right)-1)}---\left(25\right)$

将式(25)代入式(4)中以对观测噪声方差阵进行估计并输出给测更新环节104，从而可对其中的观测噪声方差阵进行更新。

步骤五：估计过程噪声方差阵；

本发明利用模糊推理系统3对过程噪声方差阵进行自适应调节，以满足跟踪环路高机动条件下对跟踪精度的要求。它的具体实现步骤为：

(1)量测更新环节104输出的状态估计值

包含了载波相位差三阶导数ω₂，加加速度变化率检测环节303利用M_j个采样点ω₂可计算载体加加速度变化率δ′：

${\mathrm{\δ}}^{\′}=\frac{1}{\frac{1}{2}\·2\mathrm{\π}{M}_j{T}_I}(\underset{i=\frac{M_j}{2}+1}{\overset{M_j}{\mathrm{\Σ}}}{\mathrm{\ω}}_2\left(i\right)-\underset{i=1}{\overset{M_j/2}{\mathrm{\Σ}}}{\mathrm{\ω}}_2\left(i\right))---\left(26\right)$

(2)残差协方差实际值检测环节301利用量测更新环节104中的M_e个残差采样点e计算残差协方差实际值：

$C_e=\frac{1}{M_e}\underset{i=1}{\overset{M_r}{\mathrm{\Σ}}}{e}_i{e}_i^T---\left(27\right)$

式中，M_e为估计残差协方差的窗口宽度。

然后残差比较环节302将C_e与量测更新环节104中的残差协方差理论值

进行比较，令：

$r^{\′}=\frac{\mathrm{tr}\left(C_e\right)}{\mathrm{tr}\left(P_{\tilde{Z}\tilde{Z}}\right)}---\left(28\right)$

式中，tr(C_e)、

分别表示对C_e与

求迹。

r′反映了残差协方差实际值C_e与残差协方差理论值

的匹配程度。

(3)模糊化环节304分别将δ′和r′转化为模糊量δ和r，并输出至作为模糊规则305，模糊规则305基于协方差匹配原理，通过判断δ与r的大小推理生成输出的模糊量μ，模糊规则如表1所示。

(4)解模糊环节306利用重心法把输出的模糊量μ转化为精确的输出变量μ′，通过式Q＝10^μ′Q₀来对Q进行自适应调节，从而输出给时间更新环节103对过程噪声方差阵进行更新。

步骤六：调节UT变换比例因子；

UT变换比例因子α用来调节Sigma采样点到中心点的距离，由于第k步状态协方差Pk反映了第k步状态估计值

与状态真实值

之间的关系，即：

$P_k=E\left[({\stackrel{\‾}{X}}_k-{\hat{X}}_k){({\stackrel{\‾}{X}}_k-{\hat{X}}_k)}^T\right]---\left(29\right)$

可利用量测更新环节104中的第k步状态协方差P_k近似地计算出第k步状态估计值与状态真实值

之间的距离d_k：

$d_k\≈\sqrt{\mathrm{tr}\left(P_k\right)}---\left(30\right)$

同理，根据第k+1步的状态协方差P_k+1可得到第k+1步未作UT变换的Sigma点到中心点的距离d_k+1：

$d_{k+1}\≈\sqrt{\mathrm{tr}\left(P_{k+1}\right)}---\left(31\right)$

然后将第k步状态真实值的邻域

作为第k+1步Sigma点的采样范围，从而得到UT变换比例因子α，即：

$\mathrm{\α}=\frac{d_k}{d_{k+1}}---\left(32\right)$

将上式所得到的结果通过Sigma点采样环节102对UT变换比例因子α进行更新，可实现控制采样范围的目标。

步骤七：补偿状态估计偏差；

状态补偿器5用来解决UKF滤波器1滞后性所带来的状态估计偏差的问题。量测更新环节104中的状态估计值

是基于前一时刻的输入载波信号得到的，因而状态估计值

中的载波信息与当前时刻输入载波信息在时间上相差一个预检测积分时间T_I。状态补偿器5将量测更新环节104中状态估计值

作为输入量，根据泰勒级数展开公式对当前时刻载波信息U(包括载波NCO6相位调整量与角频率调整量)进行预测估计，当前时刻载波信息U与前一时刻状态估计值

之间的传递关系为：

$U=\left[\begin{array}{c}\mathrm{\Δ\θ}\\ \mathrm{\Δ}{\mathrm{\ω}}_0\end{array}\right]=\left[\begin{array}{cccc}1& T_I& T_I^2/2& T_I^3/6\\ 0& 1& T_I& T_I^2/2\end{array}\right]\·\hat{X}---\left(33\right)$

进而当前时刻载波信息U对载波NCO6中的载波频率与相位进行调节，实现本地同相载波与输入数字中频信号中载波频率、相位精确对齐。

步骤八：载波NCO6辅助码NCO8。

载波环的颤动噪声比码环的颤动噪声小三个数量级，因而可以精确地为码环提供载波辅助。载波NCO6输出载波多普勒频移通过比例因子7转换成C/A码多普勒频移，进而对码NCO8进行频率辅助，用于辅助C/A码多普勒频移的比例因子为：

$\mathrm{\ρ}=\frac{f_{\mathrm{dc}}}{f_{\mathrm{dl}}}=\frac{f_c}{f_l}---\left(34\right)$

式中，ρ为比例因子，f_dl表示比例因子7输入的载波多普勒频移，f_dc表示比例因子7输出的C/A码多普勒频移，f_c为码片速率，f_l为载波频率。

比例因子7通过调整码NCO8内的C/A码频率，使本地C/A码频率与输入数字中频信号中C/A码频率对齐。

跟踪环路在实现对输入信号中载波与C/A码的精确同步后，根据同相支路相关能量可解调出输入信号中的导航电文，输出至后面的导航解算模块，以便解算出GPS接收机的三维位置、速度和时间。

本发明的一种GPS接收机自适应跟踪环路能够根据外界干扰噪声的变化自适应调节观测噪声方差阵，有效地提高了跟踪环路在时变噪声干扰环境下的自适应能力；同时，利用模糊推理系统使残差协方差理论值与其实际值保持一致，可实现对过程噪声方差阵的自适应调节，改善了跟踪环路在高动态环境下的跟踪性能，使跟踪环路更好地满足高动态、强干扰等复杂多变环境中精确跟踪GPS信号的要求。

Claims (3)

1.一种自适应跟踪环路，其特征在于，包括UKF滤波器、观测噪声方差阵检测模块、模糊推理系统、UT变换比例因子调节模块、状态补偿器、载波NCO、比例因子、码NCO、积分和清零模块、码环鉴相器和二阶码环滤波器；

载波NCO生成的本地同相载波cos信号和本地正交载波sin信号，输入的数字中频信号分别与本地同相载波cos信号、本地正交载波sin信号进行相关，得到同相支路数据A、正交支路数据A；码NCO生成的本地C/A码，分别同相支路数据A、正交支路数据A进行相关，得到同相支路数据B和正交支路数据B；积分和清零模块对同相支路数据B、正交支路数据B进行积分和清零，得到同相支路相关能量与正交支路相关能量，同相支路相关能量输出至码环鉴相器、观测噪声方差阵检测模块和UKF滤波器，正交支路相关能量也输出至码环鉴相器、观测噪声方差阵检测模块和UKF滤波器；码环鉴相器对同相支路相关能量和正交支路相关能量进行鉴相处理，得到数字中频信号C/A码与码NCO所生成本地C/A码的码相位误差；二阶码环滤波器滤除码相位误差中的噪声信号，然后将码相位误差输出至码NCO，码NCO根据码相位误差调整本地C/A码，使本地C/A码相位与输入的数字中频信号中C/A码相位对齐；载波NCO将输入数字中频信号中载波多普勒频移输出至比例因子，比例因子进行比例运算，得到输入数字中频信号中码多普勒频移，输出至码NCO，调整码NCO内的本地C/A码频率，使本地C/A码频率与输入数字中频信号中C/A码频率对齐；

UKF滤波器包括初始化环节、Sigma点采样环节、时间更新环节和量测更新环节；初始化环节将初始状态估计值

与初始状态协方差阵P₀均设置为零，并赋值给初始观测噪声方差阵R₀和初始过程噪声方差阵Q₀，然后将初始值传输给Sigma点采样环节；Sigma点采样环节在跟踪环路首个循环周期内，根据初始状态估计值

及初始状态协方差P₀、UT变换比例因子α选取2n+1个Sigma采样点，n为状态向量的维数；如果不在首个循环周期，则Sigma点采样环节根据量测更新环节提供的状态估计值

状态协方差P与UT变换比例因子α来选取Sigma采样点，并使这些采样点的均值和协方差等于状态向量的均值和协方差，然后将Sigma采样点输出给时间更新环节；时间更新环节将Sigma采样点导入到状态方程和观测方程，结合解模糊环节输出的Q阵，所述的Q阵为过程噪声方差阵，得到观测向量预测估计值、状态预测估计值及状态协方差预测估计值，并输出给量测更新环节；量测更新环节根据同相支路相关能量、正交支路相关能量，并结合观测噪声方差阵计算环节输出的R阵，所述的R阵为观测噪声方差阵，对状态预测估计值及状态协方差预测估计值进行滤波更新，得到状态估计值

状态协方差P、残差e、输入数字中频信号中载波与本地同相载波相位差的三阶导数ω₂和残差协方差理论值

将状态估计值

状态协方差P反馈给Sigma点采样环节，并将状态协方差P输出至UT变换比例因子调节模块，残差e输出至残差协方差实际值检测环节，输入数字中频信号中载波与本地同相载波相位差的三阶导数ω₂输出至加加速度变化率检测环节，残差协方差理论值

输出至残差比较环节；

观测噪声方差阵检测模块包括载波功率检测环节、载噪比检测环节和观测噪声方差阵计算环节；载波功率检测环节检测同相支路相关能量，获得数字中频信号的载波功率C，将载波功率C传输给观测噪声方差阵计算环节；同时，载噪比检测环节对同相支路相关能量与正交支路相关能量分别进行相干积分和非相干积分，得到窄带功率和宽带功率，并将窄带功率与宽带功率作比例运算后得到载噪比C/N₀，然后传输给观测噪声方差阵计算环节；观测噪声方差阵计算环节利用获得的载波功率C和载噪比C/N₀，并结合设置的预检测积分时间估计出观测噪声方差阵，进而输出R阵对量测更新环节中的观测噪声方差阵进行更新；

模糊推理系统包括残差协方差实际值检测环节、残差比较环节、加加速度变化率检测环节、模糊化环节、模糊规则和解模糊环节；残差协方差实际值检测环节通过检测量测更新环节中的残差e得到残差协方差实际值C_e，输出至残差比较环节，残差比较环节将C_e与量测更新环节中的残差协方差理论值

进行比较，将二者的比值r′输出至模糊化环节，同时，加加速度变化率检测环节根据量测更新环节中的状态估计值实时提取数字中频信号中载波与本地同相载波相位差的三阶导数ω₂，进而得到载体的加加速度变化率δ′，输出至模糊化环节，模糊化环节将输入量r′、δ′转化为模糊量r、δ，并输出给模糊规则，模糊规则基于协方差匹配原理，将模糊量r与δ推理生成模糊量μ，解模糊环节利用重心法把模糊规则输出的模糊量μ转化为精确的输出变量μ′，得到精确的输出量μ′后，通过式Q＝10^μ′Q₀来对Q阵进行估计，并输出给时间更新环节，完成对Q阵的更新；

UT变换比例因子调节模块通过量测更新环节中的状态协方差P得到第k步状态估计值与状态真实值之间的距离以及第k+1步未作比例变换的Sigma点到中心点的距离，并将第k步状态真实值的邻域作为第k+1步Sigma点的采样范围，获得UT比例变换因子α，进而输出给Sigma点采样环节；

状态补偿器用来修正UKF滤波器的滞后性所造成的状态估计偏差，它将量测更新环节所得到的前一时刻的状态估计值

作为输入量，根据状态估计值

中的载波相位差θ及其一阶导数ω₀、二阶导数ω₁和三阶导数ω₂，利用泰勒级数展开式来预测当前时刻的载波信息U，将其输出给载波NCO，载波NCO根据状态补偿器输出的当前时刻载波信息U对本地载波频率与相位进行调节，使之与数字中频信号的载波频率、相位精确对齐；

跟踪环路在完成C/A码与载波的精确同步后，由同相支路相关能量解调出接收信号中的导航电文。

2.根据权利要求1所述的一种自适应跟踪环路，其特征在于，所述的模糊规则中输出模糊量μ与输入模糊量r、输入模糊量δ的关系如表1所示：

表1控制规则

3.一种自适应跟踪环路的实现方法，其特征在于，具体包括以下步骤：

步骤一：信号相关及积分、清零处理；

输入的数字中频信号分别与本地同相载波和本地正交载波相关后进入C/A码跟踪环路，经过本地C/A码解扩以及积分和清零模块后，得到同相支路相关能量I_P(i)和正交支路相关能量Q_P(i)，表示为：

$\left.\begin{array}{c}{I}_P\left(i\right)=0.5\·K_C\·A\·D_i\·R\left[{\mathrm{\ϵ}}_i\right]\·\frac{\sin (\mathrm{\π}{T}_I\·\mathrm{\Δf})}{\mathrm{\π}{T}_I\·\mathrm{\Δf}}\·\cos \left(\mathrm{\θ}\right)+N_I\left(i\right)\\ Q_P\left(i\right)=0.5\·K_C\·A\·D_i\·R\left[{\mathrm{\ϵ}}_i\right]\·\frac{\sin (\mathrm{\π}{T}_I\·\mathrm{\Δf})}{\mathrm{\π}{T}_I\·\mathrm{\Δf}}\·\sin \left(\mathrm{\θ}\right)+N_Q\left(i\right)\end{array}\right\}---\left(1\right)$

式中，K_C为一个预检测积分时间内累加的次数，K_C＝T_I/T_S，T_I为预检测积分时间，T_S为数字中频信号的采样时间间隔，A为输入中频信号的幅值，D_i为导航电文，R[ε_i]是归一化的C/A码自相关函数，ε_i为本地C/A码与输入C/A码的相位差，

表示本地同相载波与输入信号中载波的相位差，Δf、

分别为积分间隔起始时刻本地参考信号与输入信号之间的载波频率差和相位差，N_I(i)、N_Q(i)分别为同相支路高斯白噪声和正交支路的高斯白噪声；

步骤二：码相位跟踪；

码环鉴相器根据积分和清零模块输出的数据得到码相位误差，将得到的码相位误差输出至二阶码环滤波器，码相位误差经过二阶码环滤波器滤波处理后，对码NCO输出控制信号，码NCO根据控制信号和比例因子输出的辅助参数调整本地C/A码相位和频率，使得本地C/A码与数字中频信号中的C/A码精确对齐；

步骤三：UKF滤波器建模；

将积分和清零模块输出的I_P(i)和Q_P(i)作为UKF滤波器的输入观测向量Z_k，对UKF滤波器中的状态估计值进行更新，下面对UKF滤波器进行精确建模；

(1)状态方程与观测方程

UKF滤波器采用四维状态向量对载波信号进行跟踪，其状态向量X分别为本地同相载波与输入载波相位差θ及其一阶导数ω₀、二阶导数ω₁和三阶导数ω₂，将状态向量根据泰勒级数展开公式可得系统状态方程为：

X_k＝Φ_k，k-1X_k-1+G_k，k-1U_k-1+W_k    (2)

即： $X_k={\left[\begin{array}{c}\mathrm{\θ}\\ {\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ {\mathrm{\ω}}_2\end{array}\right]}_k=\left[\begin{array}{cccc}1& T_I& T_I^2/2& T_I^2/6\\ 0& 1& T_S& T_I^2/2\\ 0& 0& 1& T_I\\ 0& 0& 0& 1\end{array}\right]{\left[\begin{array}{c}\mathrm{\θ}\\ {\mathrm{\ω}}_0\\ {\mathrm{\ω}}_1\\ {\mathrm{\ω}}_2\end{array}\right]}_{k-1}+\left[\begin{array}{cc}-1& 0\\ 0& -1\\ 0& 0\\ 0& 0\end{array}\right]{\left[\begin{array}{c}\mathrm{\Δ\θ}\\ \mathrm{\Δ}{\mathrm{\ω}}_0\end{array}\right]}_{k-1}+{\left[\begin{array}{c}{\mathrm{\η}}_1\\ {\mathrm{\η}}_2\\ {\mathrm{\η}}_3\\ {\mathrm{\η}}_4\end{array}\right]}_k$

式中，X_k、X_k-1分别为第k步状态向量和第k-1步状态向量，Φ_k，k-1为系统的状态转移矩阵，T_I为预检测积分时间，G_k，k-1为输入转移矩阵，U_k-1为控制输入量，Δθ、Δω₀分别为载波NCO6相位调整量和角频率调整量，W_k为动态模型噪声，i＝1，…，4，为泰勒级数展开式的余项，用来描述上述模型受到某些随机干扰以及模型不准确所造成的影响，Y(t)表示连续相位过程的四阶变化率，被视为单边谱密度为N_y的高斯白噪声；

进而由系统的观测向量Z_k与状态向量可得到系统的观测方程：

$Z_k={\left[\begin{array}{c}{I}_P\left(i\right)\\ Q_P\left(i\right)\end{array}\right]}_k=\left[\begin{array}{c}{K}^{\′}\cos (L\·X_k)\\ K^{\′}\sin (L\·X_k)\end{array}\right]+V_k---\left(3\right)$

式中，

L＝[1，0，0，0]，V_k＝[N_I(i)，N_Q(i)]^T，为观测噪声向量，其方差阵为：

$R_k=E\left[V_k{V}_k^T\right]=\frac{N_0}{2T_I}\·I_{2\×2}---\left(4\right)$

式中，N₀为噪声功率谱密度，I_2×2为二阶单位阵；

(2)初始化UKF滤波器

在UKF滤波器工作的初始阶段，初始化环节将初始状态估计值

与初始状态协方差阵P₀均设置为零，并赋值给初始观测噪声方差阵R₀和初始过程噪声方差阵Q₀，使其能够满足跟踪环路在信噪比为-20dB、加加速度为5g的正常环境下完成对卫星信号的跟踪功能，然后将初始值传输给Sigma点采样环节；

(3)Sigma点采样

Sigma点采样环节根据状态估计值及状态协方差选取2n+1个Sigma采样点，Sigma采样点用矩阵χ_k-1来表示，χ_k-1表示如下：

$\left.\begin{array}{c}{\mathrm{\χ}}_{i,k-1}={\hat{X}}_{k-1},i=0\\ {\mathrm{\χ}}_{i,k-1}={\hat{X}}_{k-1}+{\left(\sqrt{(n+\mathrm{\λ})P_{k-1}}\right)}_i,i=1,...,n\\ {\mathrm{\χ}}_{i,k-1}={\hat{X}}_{k-1}-{\left(\sqrt{(n+\mathrm{\λ})P_{k-1}}\right)}_{i-n},i=n+1,...,2n\end{array}\right\}---\left(5\right)$

式中，χ_i，k-1为χ_k-1的第i列，

表示第k-1步状态估计值，n为状态向量的维数，λ＝α²(n+τ)-n是一个尺度参数，α为UT变换比例因子，即无迹变换比例因子，它决定了采样点在其均值的分布情况，τ＝3-n为次级尺度调节因子，P_k-1为第k-1步状态协方差，

表示矩阵(n+λ)P_k-1平方根的第i列；

Sigma采样点的一阶统计特性权系数

和二阶统计权系数

为：

$\left.\begin{array}{c}{W}_i^{\left(m\right)}=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},i=0\\ W_i^{\left(c\right)}=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β}),i=0\\ W_i^{\left(m\right)}=W_i^{\left(c\right)}=\frac{1}{2(n+\mathrm{\λ})},i=1,...,2n\end{array}\right\}---\left(6\right)$

式中，

为用来计算均值的权值，为用来计算协方差的权值，β为验前分布因子；

(4)时间更新

时间更新环节根据Sigma点采样环节输出的采样点和解模糊环节输出的Q阵，得到观测向量预测估计值

状态预测估计值

及状态协方差预测估计值P_k，k-1；

${\mathrm{\χ}}_{k,k-1}^{*}={\mathrm{\Φ}}_{k,k-1}{\mathrm{\χ}}_{k-1}+G_{k,k-1}{U}_{k-1}---\left(7\right)$

${\hat{X}}_{k,k-1}=\underset{i=0}{\overset{i=2n}{\mathrm{\Σ}}}{W}_i^{\left(m\right)}{\mathrm{\χ}}_{i,k,k-1}^{*}---\left(8\right)$

$P_{k,k-1}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[{\mathrm{\χ}}_{i,k,k-1}^{*}-{\hat{X}}_{k,k-1}]{[{\mathrm{\χ}}_{i,k,k-1}^{*}-{\hat{X}}_{k,k-1}]}^T+Q---\left(9\right)$

$Z_{k,k-1}^{*}=\left[\begin{array}{c}{K}^{\′}\cos (L\·{\mathrm{\χ}}_{k,k-1}^{*})\\ K^{\′}\sin (L\·{\mathrm{\χ}}_{k,k-1}^{*})\end{array}\right]---\left(10\right)$

${\hat{Z}}_{k,k-1}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}{Z}_{i,k,k-1}^{*}---\left(11\right)$

式中，

为χ_k-1经状态方程变换后得到的点集，

为

经观测方程变换后得到的预测点集，

与

分别为矩阵

与矩阵

的第i列；

(5)量测更新

量测更新环节根据时间更新环节的输出和观测向量Z_k，并结合观测噪声方差阵计算环节输出的R阵，获得残差协方差理论值

状态估计值

残差e以及状态协方差P；

$P_{\tilde{Z}\tilde{Z}}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[Z_{i,k,k-1}^{*}-{\hat{Z}}_{k,k-1}]{[Z_{i,k,k-1}^{*}-{\hat{Z}}_{k,k-1}]}^T+R---\left(12\right)$

$P_{\tilde{X}\tilde{Z}}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[{\mathrm{\χ}}_{i,k,k-1}^{*}-{\hat{X}}_{k,k-1}][Z_{i,k,k-1}^{*}-{\hat{Z}}_{k,k-1}]---\left(13\right)$

$K_k=P_{\tilde{X}\tilde{Z}}{P}_{\tilde{Z}\tilde{Z}}^{-1}---\left(14\right)$

$e_k=Z_k-{\hat{Z}}_{k,k-1}---\left(15\right)$

${\hat{X}}_k={\hat{X}}_{k,k-1}+K_k{e}_k---\left(16\right)$

$P_k=P_{k,k-1}-K_k{P}_{\tilde{Z}\tilde{Z}}{K}_k^T---\left(17\right)$

式中，为观测向量预测估计值与状态预测估计值的互协方差，K_k为UKF滤波器1的增益，e_k为第k步的残差，

为第k步的状态估计值，P_k为第k步的状态协方差。

步骤四：估计观测噪声方差阵；

载波功率检测环节通过对式(1)中的同相支路相关能量I_P(i)进行计算可直接获得输入中频载波信号的功率；在式(1)中，为了消除噪声N_I(i)和导航电文翻转对计算载波功率C的影响，用K_A个I_P(i)采样点的绝对值来求取载波功率：

$C=\frac{1}{2}{A}^2=\frac{1}{2}\·{\left(\frac{\underset{i=1}{\overset{K_A}{\mathrm{\Σ}}}\left|I_P\left(i\right)\right|}{0.5K_C\·K_A}\right)}^2---\left(18\right)$

同时，载噪比检测环节在M次预检测积分时间内对同相支路相关能量I_P(i)与正交支路相关能量Q_P(i)分别进行相干积分和非相干积分，得到窄带功率NBP和宽带功率WBP：

$\mathrm{NBP}={\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{I}_P\left(i\right)\right)}_j^2+{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{Q}_P\left(i\right)\right)}^2---\left(19\right)$

$\mathrm{WBP}=\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}(I_P^2\left(i\right)+Q_P^2\left(i\right))\right)---\left(20\right)$

式中，M表示噪声带宽所改变的倍数，WBP的带宽为1/T_I，NBP的带宽为1/(MT_I)；

则宽带功率WBP和窄带功率NBP的数学期望分别为：

E(NBP)＝8M²(C+N₀/(MT_I))    (21)

E(WBP)＝8M(C+N₀/T_I)        (22)

用N个NBP和WBP采样点的平均值来分别表示E(NBP)和E(WBP)，则E(NBP)与E(WBP)的比值E(NP)表示为：

$E\left(\mathrm{NP}\right)=\frac{E\left(\mathrm{NBP}\right)}{E\left(\mathrm{WBP}\right)}=\frac{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{NBP}}_j}{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{WBP}}_j}=\frac{{\mathrm{MT}}_I\·C/N_0+1}{T_I\·C/N_0+1}---\left(23\right)$

对式(23)进行运算便可得到信号载噪比C/N₀：

$C/N_0=\frac{E\left(\mathrm{NP}\right)-1}{M-E\left(\mathrm{NP}\right)}\·\frac{1}{T_I}---\left(24\right)$

综合式(18)和式(24)可知噪声功率谱密度为N₀：

$N_0=\frac{T_I}{2}\·\frac{{\left(\underset{i=1}{\overset{K_A}{\mathrm{\Σ}}}\right|I_P\left(i\right)\left|\right)}^2\·(M-E\left(\mathrm{NP}\right))}{K_A^2\·(E\left(\mathrm{NP}\right)-1)}---\left(25\right)$

将式(25)代入式(4)中以对观测噪声方差阵进行估计并输出给测更新环节，从而可对其中的观测噪声方差阵进行更新；

步骤五：估计过程噪声方差阵；

利用模糊推理系统对过程噪声方差阵进行自适应调节，以满足跟踪环路高机动条件下对跟踪精度的要求，它的具体实现步骤为：

(1)量测更新环节输出的状态估计值

包含了载波相位差三阶导数ω₂，加加速度变化率检测环节利用M_j个采样点ω₂计算载体加加速度变化率δ′：

${\mathrm{\δ}}^{\′}=\frac{1}{\frac{1}{2}\·2\mathrm{\π}{M}_j{T}_I}(\underset{i=\frac{M_j}{2}+1}{\overset{M_j}{\mathrm{\Σ}}}{\mathrm{\ω}}_2\left(i\right)-\underset{i=1}{\overset{M_j/2}{\mathrm{\Σ}}}{\mathrm{\ω}}_2\left(i\right))---\left(26\right)$

(2)残差协方差实际值检测环节利用量测更新环节中的M_e个残差采样点e计算残差协方差实际值：

$C_e=\frac{1}{M_e}\underset{i=1}{\overset{M_r}{\mathrm{\Σ}}}{e}_i{e}_i^T---\left(27\right)$

式中，M_e为估计残差协方差的窗口宽度；

然后残差比较环节将C_e与量测更新环节中的残差协方差理论值

进行比较，令：

$r^{\′}=\frac{\mathrm{tr}\left(C_e\right)}{\mathrm{tr}\left(P_{\tilde{Z}\tilde{Z}}\right)}---\left(28\right)$

式中，tr(C_e)、

分别表示对C_e与

求迹；r′反映了残差协方差实际值C_e与残差协方差理论值

的匹配程度；

(3)模糊化环节分别将δ′和r′转化为模糊量δ和r，并输出至作为模糊规则，模糊规则如表1所示；

(4)解模糊环节利用重心法把输出的模糊量μ转化为精确的输出变量μ′，通过式Q＝10^μ′Q₀来对Q进行自适应调节，从而输出给时间更新环节对过程噪声方差阵进行更新；

步骤六：调节UT变换比例因子；

UT变换比例因子α用来调节Sigma采样点到中心点的距离，由于第k步状态协方差Pk反映了第k步状态估计值

与状态真实值之间的关系，即：

$P_k=E\left[({\stackrel{\‾}{X}}_k-{\hat{X}}_k){({\stackrel{\‾}{X}}_k-{\hat{X}}_k)}^T\right]---\left(29\right)$

利用量测更新环节中的第k步状态协方差P_k计算出第k步状态估计值

与状态真实值

之间的距离d_k：

$d_k\≈\sqrt{\mathrm{tr}\left(P_k\right)}---\left(30\right)$

同理，根据第k+1步的状态协方差P_k+1可得到第k+1步未作UT变换的Sigma点到中心点的距离d_k+1：

$d_{k+1}\≈\sqrt{\mathrm{tr}\left(P_{k+1}\right)}---\left(31\right)$

然后将第k步状态真实值的邻域

作为第k+1步Sigma点的采样范围，从而得到UT变换比例因子α，即：

$\mathrm{\α}=\frac{d_k}{d_{k+1}}---\left(32\right)$

将上式所得到的结果通过Sigma点采样环节对UT变换比例因子α进行更新，可实现控制采样范围的目标；

步骤七：补偿状态估计偏差；

量测更新环节中的状态估计值

是基于前一时刻的输入载波信号得到的，因而状态估计值中的载波信息与当前时刻输入载波信息在时间上相差一个预检测积分时间T_I；状态补偿器将量测更新环节中状态估计值

作为输入量，根据泰勒级数展开公式对当前时刻载波信息U进行预测估计，当前时刻载波信息U与前一时刻状态估计值之间的传递关系为：

$U=\left[\begin{array}{c}\mathrm{\Δ\θ}\\ \mathrm{\Δ}{\mathrm{\ω}}_0\end{array}\right]=\left[\begin{array}{cccc}1& T_I& T_I^2/2& T_I^3/6\\ 0& 1& T_I& T_I^2/2\end{array}\right]\·\hat{X}---\left(33\right)$

利用当前时刻载波信息U对载波NCO中的载波频率与相位进行调节，实现本地同相载波与输入数字中频信号中载波频率、相位精确对齐；

步骤八：载波NCO辅助码NCO；

载波NCO输出载波多普勒频移通过比例因子转换成C/A码多普勒频移，进而对码NCO进行频率辅助，用于辅助C/A码多普勒频移的比例因子为：

$\mathrm{\ρ}=\frac{f_{\mathrm{dc}}}{f_{\mathrm{dl}}}=\frac{f_c}{f_l}---\left(34\right)$

式中，ρ为比例因子，f_dl表示比例因子输入的载波多普勒频移，f_dc表示比例因子输出的C/A码多普勒频移，f_c为码片速率，f_l为载波频率；

比例因子通过调整码NCO内的C/A码频率，使本地C/A码频率与输入数字中频信号中C/A码频率对齐；

跟踪环路在实现对输入信号中载波与C/A码的精确同步后，根据同相支路相关能量解调出输入信号中的导航电文，输出至后面的导航解算模块。

Images