CN105699992A - 高动态gnss载波信号跟踪方法及系统 - Google Patents

Abstract

本发明提出一种高动态GNSS载波信号跟踪方法，时频图像可以同时结合时间域和频率域来分析信号，在原理上平衡频率分辨率与时间分辨率之间的矛盾，将高动态引起的信号不确定性限制在数个像素内，因此，与传统跟踪方法不同，本发明的动态跟踪能力不受环路带宽的限制，使得接收机可以对极高动态的载波信号进行跟踪，扩展了基于北斗/GPS的GNSS接收机的使用范围。

Description

高动态GNSS载波信号跟踪方法及系统

技术领域

本发明涉及卫星导航技术领域，特别涉及的是高动态GNSS载波信号跟踪方法及系统。

背景技术

GNSS(GlobalNavigationSatelliteSystem)即全球导航卫星系统。卫星导航接收机捕获并跟踪多个GNSS卫星的信号，然后解调其中调制的导航数据。卫星导航接收机利用测距码计算GNSS卫星与用户的相对距离，利用导航数据中的星历数据解算卫星位置和时间模型，进而计算出用户的位置。

利用GNSS信号(如GPS(全球定位系统)、北斗、GLONASS(全球卫星导航系统)等)进行导航已经广泛应用于空间低轨卫星、导弹、以及地面目标。当载体运行在高速、高加速度与高加加速度的环境中时，信号的多普勒效应将极大的影响接收机对信号的接收，因此对接收机中的信号跟踪方法和其跟踪性能均提出了极高的要求。

美国喷气动力实验(JPL)较早地在高动态GPS信号跟踪技术方面作了深入的研究并于1987年提出采用最大似然估计(MLE)方法对伪码延时和载波多普勒频率进行了估计。自此之后，高动态环境下的GNSS信号接收技术主要可以归为两类：一是研究适合高动态环境下的GNSS信号跟踪的频率估计算法，将算法嵌入接收机的载波环路内，以适应高动态环境下GNSS信号的接收。显然它具有体积小、成本低等特点。另一种是将接收机与惯导系统结合，给GNSS接收机提供IMU(惯性测量单元)的速度辅助(即提供多普勒频移的先验知识)，组成捷联惯导系统。在这种导航系统中，惯性导航系统为GNSS提供速率信息，以适应高动态定位导航；GNSS为惯导系统提供时间标准等信息，以消去惯性导航器件因时间而累积的误差。

高动态GNSS信号与普通GNSS信号的最大区别在于载波多普勒的变化率很大。具体来说，跟踪高动态GNSS信号有以下几个难点：

(1)信号的不确定性：根据信号的不确定性原理，信号的频率分辨率与时间分辨率不能同时任意高。即给定观测时间，在该时间内的载波多普勒频率估计精度越高，则该频率在观测时间内出现的具体时间就越不确定，反之亦然。因此，无论是依靠频域的传统跟踪方法，还是依靠时域的矢量跟踪方法，因为只考虑了单域下的跟踪，信号估计的不确定性将影响跟踪性能；

(2)动态水平精确检测困难：无论是依靠环路中的频率估计，还是捷联惯导系统中的IMU，都存在不能准确实时的对载波多普勒变化率进行估计的问题，影响了接收机对GNSS信号的跟踪准确性；

(3)跟踪环路设计复杂：现有的高动态跟踪方法，或是跟踪环路阶数高，运算量大，或是受IMU误差影响大，需要随时进行误差补偿，因此使得环路的设计十分复杂。

发明内容

本发明所要解决的技术问题是提供一种高动态GNSS载波信号跟踪方法，动态跟踪能力不受环路带宽的限制，可以实现极高动态的载波信号的跟踪。

为解决上述问题，本发明提出一种高动态GNSS载波信号跟踪方法，包括以下步骤：

S1：确定本地稳定载波f_L1和f_L2，对中频信号进行解调产生解调信号x₁和x₂，对解调信号x₁和x₂分段进行相关积分，产生两个N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数，其中，所述本地稳定载波f_L1根据当前多普勒频率确定，所述本地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B；

S2：根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像；

S3：若在时频图像中的瞬时频率曲线不完整，则拼接两幅时频图像，若在一幅时频图像中的瞬时频率曲线完整则无需拼接，根据上一次检测到的多普勒频率变化率确定瞬时频率搜索范围；

S4：对时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波；

S5：将瞬时频率搜索范围的时频图像进行空间变换，在变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率f′_d；

S6：根据当前载波多普勒频率f_d确定本地变动载波f_LH，选择所述解调信号x₁和x₂中包含变动载波的信号进行二次解调，以剥离载波多普勒；

S7：检测二次解调的载波相位残差以实现载波相位的跟踪；

S8：根据当前多普勒频率切换本地稳定载波f_L1和f_L2，返回执行步骤S1，产生两个新的N点离散信号O₁(k)与O₂(k),k＝2,3,...,N+1，以使环路更新并保持对载波信号的跟踪。

根据本发明的一个实施例，所述步骤S1包括以下步骤：

S11：根据相关积分的积分时间Δt，确定检测带宽B；

S12：在初始化状态下或跟踪状态下，根据当前载波多普勒频率，确定本地稳定载波f_L1，使得解调后的载波多普勒频率处于频带[f_L1，f_L1+B]的大致中间位置；在跟踪状态下，根据载波多普勒频率变化率的检测结果，设置f_L2＝f_L1+B或f_L2＝f_L1-B；

S13：在接到本地频率切换命令之前保持本地稳定载波f_L1和f_L2不变，分别使用本地稳定载波f_L1和f_L2对中频信号进行解调产生解调信号x₁和x₂，对解调信号x₁和x₂分段进行相关积分，产生两个N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N。

根据本发明的一个实施例，所述步骤S2包括以下步骤：

S21：在跟踪状态下，根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，在初始化状态下，则离散信号O₁(k)与O₂(k)都选取；

S22：对选取的离散信号进行平方计算，消除离散信号中的导航比特翻转，此时离散信号内的多普勒与多普勒变化率将相应增倍；

S23：将消除导航比特翻转后的离散信号进行Wigner-Ville变换，得到离散信号的Wigner-Ville分布；

S24：将Wigner-Ville分布中的负值全部设为0，并根据Wigner-Ville分布中的最大值将Wigner-Ville分布的值归一化为在[0,预设值]范围内的数值，从而将Wigner-Ville分布转换为时频图像，若信号中存在变动载波，则在时频图像中将生成载波信号的瞬时频率曲线，根据步骤S21中选取的离散信号，而相应生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像。

根据本发明的一个实施例，所述步骤S3包括以下步骤：

S31：根据载波多普勒频率变化率f′_d，确定变动载波信号的瞬时频率曲线在时频图像中的最大可移动频率p＝2f′_dT_s，其中，T_s为跟踪环的周期，根据p设定瞬时频率搜索范围

S32：若瞬时频率搜索范围未超出当前时频图像的频带，则取出该范围内图像，否则，将当前时频图像和与其频带紧邻的下一时频图像进行拼接，使得瞬时频率曲线完整准确，并取出瞬时频率搜索范围内图像；

S33：若有拼接，则在瞬时频率搜索范围离开两幅时频图像的拼接处之后，结束时频图像的拼接。

根据本发明的一个实施例，所述步骤S4包括以下步骤：

S41：将时频图像表示为由不同频率分量与噪声n(t)相加而成的信号：

S42：对I(t)进行调频调制：其中μ为调频调制的尺度参数；

S43：计算调制结果z(t)的Wigner-Ville变换，根据变换结果确定峰值的估计值：以滤除噪声。

根据本发明的一个实施例，所述步骤S5包括以下步骤：

S51：对瞬时频率搜索范围内的时频图像进行Radon变换或Hough变换，通过该空间变换，将时频图像由时频坐标系(t,f)投影至Radon或Hough空间坐标系(θ,λ)，其中tcosθ+fsinθ＝λ；

S52：搜索Radon或Hough空间内最大值，得到该最大值的坐标

S53：计算当前载波多普勒频率变化率与载波多普勒频率其中，Δt为积分时间。

根据本发明的一个实施例，所述步骤S8包括以下步骤：

S81：根据当前载波多普勒频率变化率确定瞬时频率曲线的移动方向；

S82：若瞬时频率曲线处于频带[f_L1，f_L1+B]内，且其移动方向为正，则设置f_L2＝f_L1+B，否则f_L2＝f_L1-B，以使图像拼接时，用于拼接的时频图像频带连续；

S83：若瞬时频率曲线移动至频带[f_L2，f_L2+B]内，则根据瞬时频率曲线的移动方向，设置f_L1＝f_L2+B或f_L1＝f_L2-B；

S84：返回执行步骤S1以使环路更新并保持对载波信号的跟踪。

本发明还提供一种高动态GNSS载波信号跟踪系统，包括：

第一载波压控振荡器和第二载波压控振荡器，分别用于产生本地稳定载波f_L1和f_L2；

初级解调器，用于根据本地稳定载波f_L1和f_L2对中频信号进行解调产生解调信号x₁和x₂；

预检测积分器，用于对所述初级解调器产生的解调信号x₁和x₂，分段进行相关积分，产生两个由N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数；

时频图像生成器，用于根据当前跟踪频率从预检测积分器产生的离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像；

时频图像拼接器，用于接收所述时频图像生成器的时频图像，在时频图像中的瞬时频率曲线不完整，拼接两幅时频图像，并用于确定瞬时频率搜索范围；

时频图像去噪器，用于对时频图像拼接器输出的时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波；

空间变换器，用于将时频图像去噪器输出的瞬时频率搜索范围的时频图像进行空间变换；

信号参数估计器，用于在空间变换器的变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率f′_d，且用于根据当前载波多普勒频率f_d与载波多普勒频率变化率f′_d，控制第一载波压控振荡器和第二载波压控振荡器对本地稳定载波f_L或f_L2进行更新，使用更新后的本地稳定载波fL或fL2对中频信号进行解调；

第三载波压控振荡器，用于根据信号参数估计器输出的当前载波多普勒频率f_d确定本地变动载波f_LH；

二次解调器，用于选择所述初级解调器产生的解调信号x₁和x₂中包含变动载波的信号，使用本地变动载波f_LH进行二次解调，以剥离载波多普勒；

鉴相器，用于检测二次解调的载波相位残差，将载波相位残差反馈给第三载波压控振荡器，以跟踪本地变动载波f_LH；

其中，所述本地稳定载波f_L1根据当前多普勒频率确定，所述本地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B。

根据本发明的一个实施例，所述空间变换器为Radon变换器或Hough变换器。

采用上述技术方案后，本发明相比现有技术具有以下有益效果：

(1)和单域方式不同，本发明的时频图像可以同时结合时间域和频率域来分析信号，在原理上平衡频率分辨率与时间分辨率之间的矛盾，将高动态引起的信号不确定性限制在数个像素内，因此，与传统跟踪方法不同，本发明的动态跟踪能力不受环路带宽的限制，使得接收机可以对极高动态的载波信号进行跟踪，扩展了基于北斗/GPS的GNSS接收机的使用范围；

(2)Wigner-Ville变换(信号中心协方差函数的傅立叶变换)在所有的时频变换中具有最佳的能量聚集性，且本发明中的信号参数估计精度与动态程度无关，高动态的大小不会影响信号的检测概率与估计精度，提高了信号跟踪的鲁棒性；

(3)实现了对载波多普勒变化率的实时估计，不再需要依靠IMU的辅助，极大的简化了跟踪环路的设计。

附图说明

图1是本发明实施例的高动态GNSS载波信号跟踪方法的流程示意图；

图2是本发明实施例的对x₁与x₂其中一个信号，从相关积分到时频图像生成的过程示意图；

图3是本发明实施例的确定瞬时频率搜索范围与时频图像拼接的过程示意图；

图4是本发明实施例的时频图像去噪前后的效果对比示意图；

图5是本发明实施例的时频图像空间变换的示意图；

图6是本发明实施例的高动态GNSS载波信号跟踪系统的结构框图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。

在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施的限制。

图1示出本发明实施例的高动态GNSS载波信号跟踪方法，包括以下步骤：

S1：确定本地稳定载波f_L1和f_L2，对中频信号进行解调产生解调信号x₁和x₂，对解调信号x₁和x₂分段进行相关积分，产生两个由N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数，其中，所述本地稳定载波f_L1根据当前多普勒频率确定，所述本地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B；

S2：根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像；

S3：若在时频图像中的瞬时频率曲线不完整，则拼接两幅时频图像，若在一幅时频图像中的瞬时频率曲线完整则无需拼接，确定瞬时频率搜索范围；

S4：对时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波；

S5：将瞬时频率搜索范围的时频图像进行空间变换，在变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率f′_d；

S6：根据当前载波多普勒频率f_d确定本地变动载波f_LH，选择所述解调信号x₁和x₂中包含变动载波的信号进行二次解调，以剥离载波多普勒；

S7：检测二次解调的载波相位残差以跟踪本地变动载波f_LH；

S8：根据当前多普勒频率切换本地稳定载波f_L1和f_L2，返回执行步骤S1，产生两个新的N点离散信号O₁(k)与O₂(k),k＝2,3,...,N+1，以使环路更新并保持对载波信号的跟踪。

下面结合附图1-5对本发明的方法进行详细的描述，可以通过步骤S1-S8实现高动态GNSS载波信号的跟踪初始化和跟踪。

在步骤S1中，可以通过载波压控震荡器产生本地稳定载波f_L1和f_L2，本地稳定载波f_L1可以根据信号捕获阶段确定的多普勒频率确定，本地稳定载波f_L2可以根据本地稳定载波f_L1和检测带宽B确定，检测带宽B根据对解调信号x₁和x₂分段进行相关积分的积分时间而定。使用本地稳定载波f_L1和f_L2对中频信号进行解调，产生解调信号x₁和x₂，可以使用预检测积分器对解调信号x₁和x₂分段进行相关积分，产生两个由N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数。

更具体的，步骤S1进一步包括以下步骤：

S11：根据相关积分的积分时间Δt，确定检测带宽B；

S12：在信号接收机成功捕获信号后，包括初始化状态和跟踪状态，根据当前载波多普勒频率使载波压控震荡器产生本地稳定载波f_L1，使得解调后的载波多普勒频率处于频带[f_L1，f_L1+B]的大致中间位置，并在接到本地频率切换命令之前保持不变，进入到后续处理流程，使得环路进入跟踪状态；在跟踪状态下，根据载波多普勒频率变化率的检测结果，设置f_L2＝f_L1+B或f_L2＝f_L1-B，使得O₁(k)与O₂(k)所处的频带紧邻；

S13：在接到本地频率切换命令之前保持本地稳定载波f_L1和f_L2不变，跟踪锁定，分别使用本地稳定载波f_L1和f_L2对解扩后的中频信号进行初始解调产生解调信号x₁和x₂，去除信号中的固定载波，使用预检测积分器对解调信号x₁和x₂分别分段进行相关积分，产生两个N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，执行后续的处理流程。

在步骤S2中，根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，在初始状态下则可以不用选取，将离散信号O₁(k)与O₂(k)均做消除导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像，若瞬时频率曲线完全处于一幅时频图像中时，则无需进行拼接，若瞬时频率曲线一半处于一幅时频图像中，而另一半处于另一幅时频图像中，则对两幅时频图像进行拼接，由于地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B，因而频带可接合。

更具体的，步骤S2进一步包括以下步骤：

S21：在跟踪状态下，根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，若两个离散信号中均存在载波则两个都选取，在初始化状态下，则离散信号O₁(k)与O₂(k)都选取；

S22：对选取的离散信号进行平方计算，消除离散信号中的导航比特翻转，此时离散信号内的多普勒与多普勒变化率将加倍；

S23：将消除导航比特翻转后的离散信号进行Wigner-Ville变换，得到离散信号的Wigner-Ville分布；

S24：将Wigner-Ville分布中的负值全部设为0，并根据Wigner-Ville分布中的最大值将Wigner-Ville分布的值归一化为在[0,预设值]范围内的数值，例如归一化为在[0,255]范围内的整型数，从而将Wigner-Ville分布转换为时频图像，若信号中存在变动载波，则在时频图像中将生成载波信号的瞬时频率曲线，根据步骤S21中选取的离散信号，而相应生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像。

参看图2，给出了本发明一个实施例中对x₁与x₂其中一个信号，从相关积分到时频图像生成的过程。

首先，随着时间流逝，中频信号经解调后持续进入由N段组成，每段积分时间为Δt的预检测积分器进行相关积分，得到离散信号O。

其次，通过计算其平方消除信号中的导航电文比特翻转，得到O²。

然后，对O²计算其Wigner-Ville变换。Wigner-Ville变换由以下公式给出：

$W_X(t,f)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}x(t+\mathrm{\τ}/2)x^{*}(t+\mathrm{\τ}/2)e^{-j2\mathrm{\π}f\mathrm{\τ}}d\mathrm{\τ};$

若给定高动态载波信号为其中，P为信号功率，a与b分别为卫星视线方向上的加速度与加加速度，f₀为初始多普勒频率，为初始载波相位，n(t)为带限白噪声，则将其带入Wigner-Ville变换公式可得：

其中，N_X为噪声n(t)的时频变换。由该结果可知，在Wigner-Ville分布中，沿着信号的瞬时频率所在坐标，将出现一系列冲击函数。

最后，对信号的Wigner-Ville分布，将负值全部设为0，并根据分布中的最大值将Wigner-Ville分布的值归一化为在[0,255]范围内的整型数，从而将分布转换为图像。

这样，沿着信号的瞬时频率所在坐标处的像素，将出现峰值，在图像中观察，信号的瞬时频率将呈现为一条曲线。此时，对信号的参数估计就转化为了图像中的曲线检测，从而可以采用图像的模式识别技术进行检测。

与传统的方法不同，信号的不确定度不再分布于整个信号带宽中，而是被限制在瞬时频率所处的时频图像像素点的周围，因此本发明可以直接提高多普勒频率、频率变化率的实时估计精度。

在步骤S3中，若在时频图像中的瞬时频率曲线不完整，则拼接两幅相邻的时频图像，使得频带相接，若在一幅时频图像中的瞬时频率曲线完整则无需拼接，确定瞬时频率搜索范围。

具体的，步骤S3进一步包括以下步骤：

S31：根据载波多普勒频率变化率f′_d，确定变动载波信号的瞬时频率曲线在时频图像中的最大可移动频率p＝2f′_dT_s，其中，T_s为跟踪环的周期，根据p设定瞬时频率搜索范围

S32：若瞬时频率搜索范围未超出当前时频图像的频带，则取出该范围内图像，送入后续处理程序执行，否则，将当前时频图像和与其频带紧邻的下一时频图像进行拼接，使得瞬时频率曲线完整准确，并取出瞬时频率搜索范围内图像，送入后续处理程序执行；

S33：若有拼接，则在瞬时频率搜索范围离开两幅时频图像的拼接处之后，结束时频图像的拼接。

参看图3，示出了本发明一个实施例中，确定瞬时频率搜索范围过程，与根据本地频率切换控制的时频图像拼接过程(为更好的显示效果，图像进行了反相)。

设Δt＝0.2ms，则单个时频图像的频带B为2500Hz。本实例中多普勒频率变化率f′_d为正，因此设置f_L2＝f_L1+B，使解调信号x₁与x₂对应的时频图像的频带毗邻，且总的多普勒频带为0～5000Hz。

如图3a所示，随着时间的流逝，瞬时频率曲线将从低频位置向高频位置移动，而当瞬时频率曲线处于时频图像1中时，由于解调信号x₂的欠采样，在时频图像2中将出现一条伪曲线。另一方面，如图3c所示，当瞬时频率曲线处于时频图像2中时，由于对解调信号x₁解调产生的负频率，在时频图像1中将同样出现一条伪曲线。伪曲线会干扰对真实瞬时频率曲线的识别，因此需要限定瞬时频率搜索范围。

根据载波多普勒频率变化率的检测结果，可以确定载波信号的瞬时频率曲线在时频图像中的最大可移动频率p＝2f′_dT_s。其中，T_s为跟踪环的周期。因此，可以根据p设定瞬时频率搜索范围为

如图3b所示，为了保证瞬时频率曲线移动时的完整连续，当其搜索范围移动至两幅时频图像的交会频带时，需要对两幅时频图像进行拼接。如图3c所示，在瞬时频率搜索范围离开两幅时频图像的交会频带后，再对图像进行分割，将时频图像1释放出来，等待更新。若瞬时频率曲线继续向高频位置移动，则更新f_L1＝f_L2+B，使解调信号x₁对应生成新的频带为5000Hz～7500Hz的时频图像3(也即新的时频图像1)，等待与时频图像2拼接，以保证瞬时频率曲线移动时的完整连续。

在步骤S4中，对时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波，实现图像去噪。

具体的，步骤S4进一步包括以下步骤：

S41：将时频图像表示为由不同频率分量与噪声n(t)相加而成的信号：

S42：对I(t)进行调频调制：其中μ为调频调制的尺度参数；

S43：计算调制结果z(t)的Wigner-Ville变换，根据变换结果确定峰值的估计值：以滤除噪声。

参看图4，示出了本发明一个实施例中，时频图像去噪前后的效果对比。由于时频峰值滤波对噪声的抑制作用，图4a中的时频图像在经历了图像去噪之后，载噪比为38dB-Hz的信号的瞬时频率曲线能够在图4b中更加清晰地显示出来。

在步骤S5中，将瞬时频率搜索范围的时频图像进行空间变换，例如可以是Radon变换(拉东变换)或Hough变换(霍夫变换)，在空间域中观察，在变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率f′_d。

具体的，步骤S5进一步包括以下步骤：

S51：对瞬时频率搜索范围内的时频图像进行Radon变换或Hough变换，通过该空间变换，将时频图像由时频坐标系(t,f)投影至Radon或Hough空间坐标系(θ,λ)，其中tcosθ+fsinθ＝λ；

S52：搜索Radon或Hough空间内最大值，得到该最大值的坐标

S53：计算当前载波多普勒频率变化率与载波多普勒频率其中，Δt为积分时间。

参看图5，示出了本发明一个实施例中，图像空间变换的结果。图5b与图5c分别为图5a中的时频图像的Radon变换与Hough变换结果。

在实际应用中，只需使用Radon变换与Hough变换中的一种。下面以Radon变换为例，说明空间变换的作用与信号参数估计的方法。

对瞬时频率搜索范围内的时频图像进行Radon变换，Radon变换由下式给出：

G(r，θ)＝R(I(t，f)}＝∫∫I(t，f)δ(λ-tcosθ-fsinθ)dtdf

通过该变换，将时频图像由时频坐标系(t,f)投影至Radon空间坐标系(θ,λ)，其中tcosθ+fsinθ＝λ。若原时频图像中存在瞬时频率曲线，则在Radon空间中将产生特定的峰值。

搜索Radon空间内峰值，得到该峰值的坐标便可计算当前载波多普勒频率变化率与载波多普勒频率

由表达式可知，载波多普勒频率与载波多普勒频率变化率的估计精度仅取决于Radon空间内峰值位置的估计精度，而与接收机所处动态程度无关，因此高动态的大小不会影响信号的检测概率与估计精度，提高了信号跟踪的鲁棒性。

有了载波多普勒频率与载波多普勒频率变化率，便可控制产生本地变动载波，对载波信号进行精确的跟踪与解调。

在步骤S6中，根据当前载波多普勒频率f_d和多普勒频率变化率，使用载波压控振荡器产生本地变动载波f_LH，本地变动载波f_LH和当前载波多普勒频率f_d大致相同，选择解调信号x₁和x₂中包含变动载波的信号进行二次解调，以完全剥离载波多普勒。

在步骤S7中，使用鉴相器检测二次解调的载波相位残差，用于对相位进行精确跟踪，以实现载波相位的跟踪。

在步骤S8中，实现环路更新和本地频率切换控制，根据当前多普勒频率切换更新本地稳定载波f_L1和f_L2，使信号的Wigner-Ville变换产生的时频图像的频率范围能够及时扩展至新的频段，保证跟踪的持续，产生两个新的N点离散信号O₁(k)与O₂(k),k＝2,3,...,N+1，返回执行步骤S1以使环路更新并保持对载波信号的跟踪。

具体的，步骤S8进一步包括以下步骤：

S81：根据当前载波多普勒频率变化率确定瞬时频率曲线的移动方向；

S82：若瞬时频率曲线处于频带[f_L1，f_L1+B]内，且其移动方向为正，则设置f_L2＝f_L1+B，否则f_L2＝f_L1-B，以使图像拼接时，用于拼接的时频图像频带准确连续；

S83：若瞬时频率曲线移动至频带[f_L2，f_L2+B]内，则根据瞬时频率曲线的移动方向，设置f_L1＝f_L2+B或f_L1＝f_L2-B；

S84：返回执行步骤S1以使环路更新并保持对载波信号的跟踪。

参看图6，本实施例的高动态GNSS载波信号跟踪系统，包括：

第一载波压控振荡器1和第二载波压控振荡器2，分别用于产生本地稳定载波f_L1和f_L2；

初级解调器11，用于根据第一载波压控振荡器1和第二载波压控振荡器2产生的本地稳定载波f_L1和f_L2，对接收的中频信号进行解调产生解调信号x₁和x₂；较佳的，中频信号为解扩后的中频信号。

预检测积分器3，用于对所述初级解调器11产生的解调信号x₁和x₂，分段进行相关积分，产生两个由N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数；

时频图像生成器4，用于根据当前跟踪频率从预检测积分器3产生的离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像；

时频图像拼接器5，用于接收所述时频图像生成器4的时频图像，在时频图像中的瞬时频率曲线不完整，拼接两幅时频图像，并用于确定瞬时频率搜索范围；

时频图像去噪器6，用于对时频图像拼接器5输出的时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波；

空间变换器7，用于将时频图像去噪器6输出的瞬时频率搜索范围的时频图像进行空间变换；

信号参数估计器8，用于在空间变换器7的变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率f′_d，且用于根据当前载波多普勒频率f_d与载波多普勒频率变化率f′_d，控制第一载波压控振荡器1和第二载波压控振荡器2对本地稳定载波f_L或f_L2进行更新，使用更新后的本地稳定载波fL或fL2对中频信号进行解调；

第三载波压控振荡器9，用于根据信号参数估计器输出的当前载波多普勒频率f_d确定本地变动载波f_LH；

二次解调器12，用于选择所述初级解调器11产生的解调信号x₁和x₂中包含变动载波的信号，使用第三载波压控振荡器9产生的本地变动载波f_LH进行二次解调，以剥离载波多普勒；

鉴相器10，用于检测二次解调的载波相位残差，将载波相位残差反馈给第三载波压控振荡器9，以实现载波相位的跟踪，更精确的跟踪变动载波；

其中，所述本地稳定载波f_L1根据当前多普勒频率确定，所述本地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B。

可选的，空间变换器7为Radon变换器或Hough变换器。

关于本发明系统部分的内容可以参看本发明前述方法中的具体描述内容，在此不再赘述。

本发明虽然以较佳实施例公开如上，但其并不是用来限定权利要求，任何本领域技术人员在不脱离本发明的精神和范围内，都可以做出可能的变动和修改，因此本发明的保护范围应当以本发明权利要求所界定的范围为准。

Claims (9)

1.一种高动态GNSS载波信号跟踪方法，其特征在于，包括以下步骤：

S1：确定本地稳定载波f_L1和f_L2，对中频信号进行解调产生解调信号x₁和x₂，对解调信号x₁和x₂分段进行相关积分，产生两个N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数，其中，所述本地稳定载波f_L1根据当前多普勒频率确定，所述本地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B；

S2：根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像；

S3：若在时频图像中的瞬时频率曲线不完整，则拼接两幅时频图像，若在一幅时频图像中的瞬时频率曲线完整则无需拼接，根据上一次检测到的多普勒频率变化率确定瞬时频率搜索范围；

S4：对时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波；

S5：将瞬时频率搜索范围的时频图像进行空间变换，在变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率

S6：根据当前载波多普勒频率f_d确定本地变动载波f_LH，选择所述解调信号x₁和x₂中包含变动载波的信号进行二次解调，以剥离载波多普勒；

S7：检测二次解调的载波相位残差以实现载波相位的跟踪；

S8：根据当前多普勒频率切换本地稳定载波f_L1和f_L2，返回执行步骤S1，产生两个新的N点离散信号O₁(k)与O₂(k),k＝2,3,...,N+1，以使环路更新并保持对载波信号的跟踪。

2.如权利要求1所述的高动态GNSS载波信号跟踪方法，其特征在于，所述步骤S1包括以下步骤：

S11：根据相关积分的积分时间Δt，确定检测带宽B；

S12：在初始化状态下或跟踪状态下，根据当前载波多普勒频率，确定本地稳定载波f_L1，使得解调后的载波多普勒频率处于频带[f_L1，f_L1+B]的大致中间位置；在跟踪状态下，根据载波多普勒频率变化率的检测结果，设置f_L2＝f_L1+B或f_L2＝f_L1-B；

S13：在接到本地频率切换命令之前保持本地稳定载波f_L1和f_L2不变，分别使用本地稳定载波f_L1和f_L2对中频信号进行解调产生解调信号x₁和x₂，对解调信号x₁和x₂分段进行相关积分，产生两个N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N。

3.如权利要求1所述的高动态GNSS载波信号跟踪方法，其特征在于，所述步骤S2包括以下步骤：

S21：在跟踪状态下，根据当前跟踪频率从离散信号O₁(k)与O₂(k)中选取具有载波的信号，在初始化状态下，则离散信号O₁(k)与O₂(k)都选取；

S22：对选取的离散信号进行平方计算，消除离散信号中的导航比特翻转，此时离散信号内的多普勒与多普勒变化率将相应增倍；

S23：将消除导航比特翻转后的离散信号进行Wigner-Ville变换，得到离散信号的Wigner-Ville分布；

S24：将Wigner-Ville分布中的负值全部设为0，并根据Wigner-Ville分布中的最大值将Wigner-Ville分布的值归一化为在[0,预设值]范围内的数值，从而将Wigner-Ville分布转换为时频图像，若信号中存在变动载波，则在时频图像中将生成载波信号的瞬时频率曲线，根据步骤S21中选取的离散信号，而相应生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像。

4.如权利要求1所述的高动态GNSS载波信号跟踪方法，其特征在于，所述步骤S3包括以下步骤：

S31：根据载波多普勒频率变化率f′_d，确定变动载波信号的瞬时频率曲线在时频图像中的最大可移动频率p＝2f′_dT_s，其中，T_s为跟踪环的周期，根据p设定瞬时频率搜索范围

S32：若瞬时频率搜索范围未超出当前时频图像的频带，则取出该范围内图像，否则，将当前时频图像和与其频带紧邻的下一时频图像进行拼接，使得瞬时频率曲线完整准确，并取出瞬时频率搜索范围内图像；

S33：若有拼接，则在瞬时频率搜索范围离开两幅时频图像的拼接处之后，结束时频图像的拼接。

5.如权利要求1所述的高动态GNSS载波信号跟踪方法，其特征在于，所述步骤S4包括以下步骤：

S41：将时频图像表示为由不同频率分量与噪声n(t)相加而成的信号： $I\left(t\right)=\underset{k=1}{\overset{N}{\Σ}}{X}_k\left(t\right)+n\left(t\right);$

S42：对I(t)进行调频调制：其中μ为调频调制的尺度参数；

S43：计算调制结果z(t)的Wigner-Ville变换，根据变换结果确定峰值的估计值： $\hat{X}\left(t\right)={\hat{f}}_Z\left(t\right)=\frac{\arg \underset{f}{\max }\[W_Z\left(t,f\right)\]}{\mathrm{\μ}},$ 以滤除噪声。

6.如权利要求1所述的高动态GNSS载波信号跟踪方法，其特征在于，所述步骤S5包括以下步骤：

S51：对瞬时频率搜索范围内的时频图像进行Radon变换或Hough变换，通过该空间变换，将时频图像由时频坐标系(t,f)投影至Radon或Hough空间坐标系(θ,λ)，其中tcosθ+fsinθ＝λ；

S52：搜索Radon或Hough空间内最大值，得到该最大值的坐标

S53：计算当前载波多普勒频率变化率与载波多普勒频率 $f_d=\frac{\widehat{\mathrm{\λ}}}{\mathrm{\Δ}tN}+\frac{f_d^{\′}\mathrm{\Δ}tN}{2},$ 其中，Δt为积分时间。

7.如权利要求1所述的高动态GNSS载波信号跟踪方法，其特征在于，所述步骤S8包括以下步骤：

S81：根据当前载波多普勒频率变化率确定瞬时频率曲线的移动方向；

S82：若瞬时频率曲线处于频带[f_L1，f_L1+B]内，且其移动方向为正，则设置f_L2＝f_L1+B，否则f_L2＝f_L1-B，以使图像拼接时，用于拼接的时频图像频带连续；

S83：若瞬时频率曲线移动至频带[f_L2，f_L2+B]内，则根据瞬时频率曲线的移动方向，设置f_L1＝f_L2+B或f_L1＝f_L2-B；

S84：返回执行步骤S1以使环路更新并保持对载波信号的跟踪。

8.一种高动态GNSS载波信号跟踪系统，其特征在于，包括：

第一载波压控振荡器和第二载波压控振荡器，分别用于产生本地稳定载波f_L1和f_L2；

初级解调器，用于根据本地稳定载波f_L1和f_L2对中频信号进行解调产生解调信号x₁和x₂；

预检测积分器，用于对所述初级解调器产生的解调信号x₁和x₂，分段进行相关积分，产生两个由N点离散信号O₁(k)与O₂(k),k＝1,2,3,...,N，N为正整数；

时频图像生成器，用于根据当前跟踪频率从预检测积分器产生的离散信号O₁(k)与O₂(k)中选取具有载波的信号，消除离散信号中的导航比特翻转并计算Wigner-Ville变换，得到离散信号的时频分布，对时频分布进行归一化，生成一幅包含完整瞬时频率曲线的时频图像或两幅分别包含瞬时频率曲线的一部分的时频图像；

时频图像拼接器，用于接收所述时频图像生成器的时频图像，在时频图像中的瞬时频率曲线不完整，拼接两幅时频图像，并用于确定瞬时频率搜索范围；

时频图像去噪器，用于对时频图像拼接器输出的时频图像在瞬时频率搜索范围内的部分进行时频峰值滤波；

空间变换器，用于将时频图像去噪器输出的瞬时频率搜索范围的时频图像进行空间变换；

信号参数估计器，用于在空间变换器的变换结果中进行二维峰值搜索，根据峰值位置估计当前载波多普勒频率f_d与载波多普勒频率变化率f′_d，且用于根据当前载波多普勒频率f_d与载波多普勒频率变化率f′_d，控制第一载波压控振荡器和第二载波压控振荡器对本地稳定载波f_L或f_L2进行更新，使用更新后的本地稳定载波fL或fL2对中频信号进行解调；

第三载波压控振荡器，用于根据信号参数估计器输出的当前载波多普勒频率f_d确定本地变动载波f_LH；

二次解调器，用于选择所述初级解调器产生的解调信号x₁和x₂中包含变动载波的信号，使用本地变动载波f_LH进行二次解调，以剥离载波多普勒；

鉴相器，用于检测二次解调的载波相位残差，将载波相位残差反馈给第三载波压控振荡器，以跟踪本地变动载波f_LH；

其中，所述本地稳定载波f_L1根据当前多普勒频率确定，所述本地稳定载波f_L2和所述本地稳定载波f_L1之间相差一个检测带宽B。

9.如权利要求8所述的高动态GNSS载波信号跟踪系统，其特征在于，所述空间变换器为Radon变换器或Hough变换器。