阅读说明：本技术 一种实时计算电离层相位闪烁指数计算方法及系统 (Method and system for calculating ionospheric phase flicker index in real time ) 是由 田羽森 王先毅 孙越强 杜起飞 刘黎军 王冬伟 李伟 白伟华 蔡跃荣 柳聪亮 孟 于 2019-11-22 设计创作，主要内容包括：本发明公开了一种实时计算电离层相位闪烁指数计算方法及系统,所述方法包括：从接收机中获取载波相位的测量值,计算信号的多普勒频率测量值；基于信号的多普勒频率测量值,得到去除卫星和接收机相对运动产生的相位变化的载波相位；用六阶Butterworth高通滤波器对载波相位进行滤波,去除卫星钟差、接收机钟差、对流程误差、以及拟合的残差；得到滤波后的载波相位；基于滤波后的载波相位,结合热噪声修正值,计算电离层相位闪烁指数。本发明的方法能够提供高精度的电离层相位闪烁指数测量值。(The invention discloses a method and a system for calculating an ionospheric phase flicker index in real time, wherein the method comprises the following steps: obtaining a measured value of a carrier phase from a receiver, and calculating a Doppler frequency measured value of a signal; based on the Doppler frequency measurement value of the signal, obtaining a carrier phase for removing phase change generated by relative motion of the satellite and the receiver; filtering the carrier phase by using a six-order Butterworth high-pass filter to remove satellite clock error, receiver clock error, flow error and fitted residual error; obtaining a filtered carrier phase; and calculating the ionospheric phase flicker index based on the filtered carrier phase and by combining the thermal noise correction value. The method can provide high-precision ionospheric phase scintillation index measurement values.)

一种实时计算电离层相位闪烁指数计算方法及系统

技术领域

本发明涉及卫星通信、导航定位领域，具体涉及一种实时计算电离层相位闪烁指数计算方法及系统。

背景技术

电离层闪烁指电离层不均匀结构引起的电磁波幅度的衰落和相位抖动，给卫星通信、导航定位等带来严重危害。因此实时监测电离层闪烁具有重要意义。电离层相位闪烁是衡量电离层闪烁的重要观测量之一。接收机测量的载波相位中除了闪烁相位，还包含卫星接收机相对运动引起的变化，接收机和卫星的时钟误差、对流层以及热噪声引起的相位抖动。

在计算电离层闪烁相位时需要去除由卫星运动带来的相位趋势。常见的一种方法是用4阶多项式拟合出相位趋势并减去，但是这种方法只应用于事后数据处理中。

发明内容

本发明的目的在于克服上述技术缺陷，提出了一种实时计算电离层相位闪烁指数计算方法。该方法可以有效去除其他因素的影响，准确测量电离层闪烁相位。

为了实现上述目的，本发明提供了一种实时计算电离层相位闪烁指数计算方法，所述方法包括：

从接收机中获取载波相位的测量值，计算信号的多普勒频率测量值；

基于信号的多普勒频率测量值，得到去除卫星和接收机相对运动产生的相位变化的载波相位；

用六阶Butterworth高通滤波器对载波相位进行滤波，去除卫星钟差、接收机钟差、对流程误差、以及拟合的残差；得到滤波后的载波相位；

基于滤波后的载波相位，结合热噪声修正值，计算电离层相位闪烁指数。

作为上述方法的一种改进，所述基于信号的多普勒频率测量值，得到去除卫星和接收机相对运动产生的相位变化的载波相位；具体包括：

根据多普勒频率记录的时间和多普勒频率计算多项式拟合的观测量；

利用最小二乘法求出多项式的系数；

根据当前时间计算出拟合的载波多普勒大小；

将拟合的载波多普勒和测量的载波多普勒相减得到了去除运动趋势的载波频率值；

对去除运动趋势的载波频率值进行积分得到载波相位。

作为上述方法的一种改进，所述根据多普勒频率记录的时间和多普勒频率计算多项式拟合的观测量；具体包括：

建立多普勒频率记录的时间和多普勒频率之间的关系式：

T(N)P(N)＝C(N)

其中，T(N)是以时间的幂次作为变量的时间矩阵；P(N)是多项式的系数；C(N)是载波多普勒向量；N表示用来进行拟合的数据的个数；

时间矩阵T(N)为：

其中，t(i)表示第i个采样点对应的时间，1≤i≤N；

载波多普勒向量C(N)为：

其中，carph(i)表示第i个采样点对应的载波相位多普勒；

多项式系数向量P(N)为：

P_i表示多项式的第i个系数。

作为上述方法的一种改进，所述利用最小二乘法求出多项式的系数；具体包括：

对曲线方程进行正规化，即T(N)P(N)＝C(N)的两侧乘上T(N)的转置T(N)^T，得到：

U(N)P(N)＝D(N)

其中U(N)为：

U(N)＝T(N)^TT(N)

其中D(N)为：

D(N)＝T(N)^TC(N)

U(n)到U(n+1)的增量ΔU(n+1)为：

其中，1≤n≤N-1；

所以U(n)的递推公式为：

U(n+1)＝U(n)+ΔU(n+1)

利用该递推公式和初始值U(1)计算得到U(N)；

D(n)到D(n+1)的增量为：

所以D(N)的递推公式为：

D(N+1)＝D(N)+ΔD(N+1)

利用该递推公式和初始值D(1)计算得到D(N)；

基于U(N)P(N)＝D(N)，利用最小二乘法求出P(N)。

作为上述方法的一种改进，所述根据当前时间计算出拟合的载波多普勒大小；具体为：

拟合的载波多普勒大小为：

其中，t为当前时间。

作为上述方法的一种改进，基于滤波后的载波相位，结合热噪声修正值，计算电离层相位闪烁指数，具体包括：

热噪声引起的相位σ_pll为：

其中，B_n是接收机的环路带宽，C/N₀表示再造比，T是环路积分时间；

经过热噪声误差修正的电离层相位闪烁指数

为：

其中，

为滤波后的载波相位。

本发明还提供了一种实时计算电离层相位闪烁指数计算系统，所述系统包括：

多普勒频率测量值计算模块，用于从接收机中获取载波相位的测量值，计算信号的多普勒频率测量值；

载波相位修正模块，用于基于信号的多普勒频率测量值，得到去除卫星和接收机相对运动产生的相位变化的载波相位；

载波相位滤波模块，用于使用六阶Butterworth高通滤波器对载波相位进行滤波，去除卫星钟差、接收机钟差、对流程误差、以及拟合的残差；得到滤波后的载波相位；

电离层相位闪烁指数机损模块，用于基于滤波后的载波相位，结合热噪声修正值，计算电离层相位闪烁指数。

本发明的优势在于：

1、本发明的方法利用多项式拟合的方法去除相位趋势并且能够在GNSS电离层闪烁仪上实时运行，此外，该方法还考虑了环路带来的测量误差，并对其修正；

2、本发明的方法运算稳定，结果准确，能够提供高精度的电离层相位闪烁指数测量值。

附图说明

图1为本发明的实时计算电离层相位闪烁指数计算方法的流程图。

具体实施方式

下面结合附图对本发明的技术方案进行详细的说明。

为了计算电离层闪烁引起的相位抖动需要滤除其他因素引起的相位抖动。首先要去除由于卫星和接收机相对运动产生的相位变化。卫星与接收机相对运动产生的相位变化一般可以用多项式拟合载波相位的方式将运动产生的相位变化趋势计算出来，并从测量得到的载波相位中减去。为了减小计算量，本方法通过二阶多项式拟合信号的载波多普勒，然后用测量的载波多普勒减去拟合的载波多普勒从而得到去除运动趋势的载波多普勒。因此，去除卫星与接收机相对运动的载波相位可以由下式计算得到：

式中的多项式的系数p_i通过最小二乘法计算得到，见式(2)，t为时间，

是接收机测量得到的载波多普勒；

UP＝D (2)

其中，U为输入的时间序列正规化矩阵，P是多项式系数的列向量，D是载波多普勒和时间序列的乘积。等式两边同时乘以U矩阵的逆便可以求得多项式系数P。

卫星钟差，接收机钟差以及中性大气引起的相位变化可以用一个6阶的Butterworth高通滤波器进行滤除，滤波器的截止频率为0.1Hz。滤波器的参数如表1所示，

表1：Butterworth滤波器系数

热噪声引起的相位可以用载噪比、环路带宽和积分时间进行估计，公式如下，

式中B_n是接收机的环路带宽，C/N₀表示载噪比，T是环路积分时间；σ_pll为热噪声引起的相位；利用估计出来的噪声抖动结合下式可以修正热噪声带来的误差：

为经过误差修正的电离层相位闪烁指数。

如图1所示，本发明的实施例1提出了一种实时计算电离层相位闪烁指数计算方法，包括：

步骤1)从接收机中获取载波相位的测量值；

步骤2)对载波相位进行差分运算，求出信号的多普勒频率测量值；

步骤3)计算根据多普勒频率记录的时间和多普勒频率计算多项式拟合的观测量；

具体计算方法如下；

T(N)P(N)＝C(N) (5)

其中T(N)是以时间的幂次作为变量的矩阵；

P(N)是多项式的系数；

C(N)是测量的载波多普勒；

N表示用来进行拟合的数据的个数。

时间矩阵T(N)为：

其中t(i)表示第i个采样点对应的时间，

载波测量向量C(N)为：

carph(i)表示第i个载波相位多普勒的采样值；

多项式系数向量P(N)为：

P_i表示多项式第i个系数；

使用最小二乘法对多项式系数进行求解，首先要对曲线方程进行正规化，即方程两侧坐标乘上T(N)的转置，

U(N)P(N)＝D(N) (9)

其中U(N)为，

U(N)＝T(N)^TT(N) (10)

其中D(N)为，

D(N)＝T(N)^TC(N) (11)

由于接收机资源有限，无法将过去的所有值存储起来，所以需要求出U(N)和D(N)的递推公式，通过迭代的方式求出U(N)和D(N)。

U(N)到U(N+1)的增量为：

所以U(N)的递推公式为，

U(N+1)＝U(N)+ΔU(N+1) (13)

D(N)到D(N+1)的增量为：

所以D(N)的递推公式为：

D(N+1)＝D(N)+ΔD(N+1) (15)

步骤4)利用最小二乘法求出多项式的系数并根据当前时间计算出拟合的载波多普勒大小；

步骤5)将拟合的载波多普勒和测量的载波多普勒相减得到了去除运动趋势的载波频率值；

步骤6)对去除运动趋势的载波频率值进行积分得到载波相位；

步骤7)用六阶Butterworth高通滤波器对载波相位进行滤波，去除卫星钟差、接收机钟差、对流程误差、以及拟合的残差；

步骤8)结合热噪声修正值，计算电离层相位闪烁指数。

实施例2

本发明的实施例2提供了一种实时计算电离层相位闪烁指数计算系统，所述系统包括：

多普勒频率测量值计算模块，用于从接收机中获取载波相位的测量值，计算信号的多普勒频率测量值；

载波相位修正模块，用于基于信号的多普勒频率测量值，得到去除卫星和接收机相对运动产生的相位变化的载波相位；

载波相位滤波模块，用于使用六阶Butterworth高通滤波器对载波相位进行滤波，去除卫星钟差、接收机钟差、对流程误差、以及拟合的残差；得到滤波后的载波相位；

电离层相位闪烁指数机损模块，用于基于滤波后的载波相位，结合热噪声修正值，计算电离层相位闪烁指数。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。