阅读说明：本技术 一种实时计算电离层相位闪烁指数计算方法及系统 (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.)

1. A method of calculating ionospheric phase scintillation indices in real-time, the method comprising:

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.

2. The method of claim 1, wherein the carrier phase is derived based on doppler frequency measurements of the signal, with phase changes resulting from relative motion between the satellite and the receiver removed; the method specifically comprises the following steps:

calculating observation quantity of polynomial fitting according to the time recorded by the Doppler frequency and the Doppler frequency;

solving coefficients of the polynomial by using a least square method;

calculating the fitted carrier Doppler size according to the current time;

subtracting the fitted carrier Doppler and the measured carrier Doppler to obtain a carrier frequency value with the movement trend removed;

and integrating the carrier frequency values without the motion tendency to obtain carrier phases.

3. The method of claim 2, wherein the calculating of the observation of the polynomial fit from the time of the doppler frequency recordings and the doppler frequency; the method specifically comprises the following steps:

establishing a relation between the time of Doppler frequency recording and the Doppler frequency:

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

wherein T (N) is a time matrix with a variable at the power of time; p (N) is the coefficient of the polynomial; c (N) is a carrier Doppler vector; n represents the number of data used to perform the fitting;

the time matrix T (N) is:

wherein t (i) represents the time corresponding to the ith sampling point, and i is more than or equal to 1 and less than or equal to N;

the carrier doppler vector c (n) is:

wherein, carph (i) represents the carrier phase doppler corresponding to the ith sampling point;

the polynomial coefficient vector P (N) is:

P_ithe ith coefficient of the polynomial.

4. The method according to claim 3, wherein the coefficients of the polynomial are found by a least squares method; the method specifically comprises the following steps:

normalizing the curve equation by multiplying the two sides of T (N) P (N) C (N) by the transpose of T (N) ((N))^TObtaining:

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

wherein U (N) is:

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

wherein D (N) is:

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

the increment Δ U (n +1) from U (n) to U (n +1) is:

wherein N is more than or equal to 1 and less than or equal to N-1;

the recurrence formula for U (n) is therefore:

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

calculating by using the recursion formula and the initial value U (1) to obtain U (N);

the increments of D (n) to D (n +1) are:

the recurrence formula for D (N) is therefore:

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

calculating by using the recursion formula and the initial value D (1) to obtain D (N);

p (n) is obtained by the least square method based on u (n) p (n) ═ d (n).

5. The method of claim 4, wherein calculating a fitted carrier Doppler magnitude from a current time; the method specifically comprises the following steps:

the fitted carrier doppler magnitudes are:

wherein t is the current time.

6. The method of claim 5, wherein calculating the ionospheric phase flicker index based on the filtered carrier phase in combination with the thermal noise correction value comprises:

phase sigma due to thermal noise_pllComprises the following steps:

wherein, B_nIs the loop bandwidth, C/N, of the receiver₀Represents the reconstruction ratio, T is the loop integration time;

ionospheric phase flicker index corrected for thermal noise error

wherein the content of the first and second substances,

7. A system for calculating an ionospheric phase scintillation index in real time, the system comprising:

the Doppler frequency measurement value calculation module is used for acquiring the measurement value of the carrier phase from the receiver and calculating the Doppler frequency measurement value of the signal;

the carrier phase correction module is used for obtaining a carrier phase for removing phase change generated by relative motion of the satellite and the receiver based on the Doppler frequency measurement value of the signal;

the carrier phase filtering module is used for 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 the ionospheric phase flicker index machine loss module is used for calculating the ionospheric phase flicker index based on the filtered carrier phase and combined with the thermal noise correction value.

Technical Field

The invention relates to the field of satellite communication and navigation positioning, in particular to a method and a system for calculating an ionospheric phase scintillation index in real time.

Background

Ionospheric scintillation refers to the electromagnetic wave amplitude fading and phase jitter caused by the uneven structure of the ionosphere, which brings serious harm to satellite communication, navigation positioning and the like. Therefore, the real-time monitoring of ionospheric scintillation is of great significance. Ionospheric phase scintillation is one of the important observables to measure ionospheric scintillation. Besides the flicker phase, the carrier phase measured by the receiver also contains the change caused by the relative motion of the satellite receiver, the clock error of the receiver and the satellite, the phase jitter caused by the troposphere and the thermal noise.

The phase trend caused by the satellite motion needs to be removed when calculating the ionospheric scintillation phase. A common method is to fit a 4 th order polynomial to the phase trend and subtract, but this method is only applied in post-hoc data processing.

Disclosure of Invention

The invention aims to overcome the technical defects and provides a method for calculating the ionospheric phase scintillation index in real time. The method can effectively remove the influence of other factors and accurately measure the ionospheric scintillation phase.

In order to achieve the above object, the present invention provides a method for calculating an ionospheric phase flicker index in real time, the method comprising:

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.

As an improvement of the above method, the carrier phase is obtained by removing phase changes caused by relative motion between the satellite and the receiver based on doppler frequency measurements of the signals; the method specifically comprises the following steps:

calculating observation quantity of polynomial fitting according to the time recorded by the Doppler frequency and the Doppler frequency;

solving coefficients of the polynomial by using a least square method;

calculating the fitted carrier Doppler size according to the current time;

subtracting the fitted carrier Doppler and the measured carrier Doppler to obtain a carrier frequency value with the movement trend removed;

and integrating the carrier frequency values without the motion tendency to obtain carrier phases.

As an improvement of the above method, the calculating of the observation of polynomial fitting from the time of doppler frequency recording and the doppler frequency; the method specifically comprises the following steps:

establishing a relation between the time of Doppler frequency recording and the Doppler frequency:

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

wherein T (N) is a time matrix with a variable at the power of time; p (N) is the coefficient of the polynomial; c (N) is a carrier Doppler vector; n represents the number of data used to perform the fitting;

the time matrix T (N) is:

wherein t (i) represents the time corresponding to the ith sampling point, and i is more than or equal to 1 and less than or equal to N;

the carrier doppler vector c (n) is:

wherein, carph (i) represents the carrier phase doppler corresponding to the ith sampling point;

the polynomial coefficient vector P (N) is:

P_ithe ith coefficient of the polynomial.

As an improvement of the above method, the coefficients of the polynomial are found by a least square method; the method specifically comprises the following steps:

normalizing the curve equation by multiplying the two sides of T (N) P (N) C (N) by the transpose of T (N) ((N))^TObtaining:

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

wherein U (N) is:

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

wherein D (N) is:

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

the increment Δ U (n +1) from U (n) to U (n +1) is:

wherein N is more than or equal to 1 and less than or equal to N-1;

the recurrence formula for U (n) is therefore:

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

calculating by using the recursion formula and the initial value U (1) to obtain U (N);

the increments of D (n) to D (n +1) are:

the recurrence formula for D (N) is therefore:

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

calculating by using the recursion formula and the initial value D (1) to obtain D (N);

p (n) is obtained by the least square method based on u (n) p (n) ═ d (n).

As an improvement of the above method, the fitted carrier doppler magnitude is calculated according to the current time; the method specifically comprises the following steps:

the fitted carrier doppler magnitudes are:

wherein t is the current time.

As an improvement of the above method, calculating an ionospheric phase flicker index based on the filtered carrier phase in combination with a thermal noise correction value specifically includes:

phase sigma due to thermal noise_pllComprises the following steps:

wherein, B_nIs the loop bandwidth, C/N, of the receiver₀Represents the reconstruction ratio, T is the loop integration time;

ionospheric phase flicker index corrected for thermal noise error

Comprises the following steps:

wherein the content of the first and second substances,

is the filtered carrier phase.

The invention also provides a system for calculating the ionospheric phase scintillation index in real time, which comprises:

the Doppler frequency measurement value calculation module is used for acquiring the measurement value of the carrier phase from the receiver and calculating the Doppler frequency measurement value of the signal;

the carrier phase correction module is used for obtaining a carrier phase for removing phase change generated by relative motion of the satellite and the receiver based on the Doppler frequency measurement value of the signal;

the carrier phase filtering module is used for 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 the ionospheric phase flicker index machine loss module is used for calculating the ionospheric phase flicker index based on the filtered carrier phase and combined with the thermal noise correction value.

The invention has the advantages that:

1. the method of the invention removes the phase trend by using a polynomial fitting method and can run on the GNSS ionosphere scintillator in real time, and in addition, the method also considers the measurement error brought by the loop and corrects the measurement error;

2. the method has stable operation and accurate result, and can provide high-precision ionospheric phase scintillation index measurement values.

Drawings

Fig. 1 is a flowchart of a method for calculating an ionospheric phase scintillation index in real time according to the present invention.

Detailed Description

The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

In order to calculate the phase jitter caused by ionospheric flicker, the phase jitter caused by other factors needs to be filtered. The phase change due to the relative motion of the satellite and the receiver is removed first. The phase change caused by the relative motion between the satellite and the receiver can be calculated by fitting a polynomial to the carrier phase, and subtracted from the measured carrier phase. In order to reduce the calculation amount, the method fits the carrier Doppler of the signal through a second-order polynomial, and then subtracts the fitted carrier Doppler from the measured carrier Doppler so as to obtain the carrier Doppler with the movement trend removed. Therefore, the carrier phase to remove the relative motion between the satellite and the receiver can be calculated by the following formula:

coefficient p of polynomial in the formula_iCalculated by the least square method, see formula (2), t is time,

is the carrier doppler measured by the receiver;

UP＝D (2)

where U is the input time series normalization matrix, P is the column vector of polynomial coefficients, and D is the product of carrier doppler and time series. The polynomial coefficient P can be found by multiplying both sides of the equation by the inverse of the U matrix.

The satellite clock error, the receiver clock error and the phase change caused by neutral atmosphere can be filtered by a 6-order Butterworth high-pass filter, and the cut-off frequency of the filter is 0.1 Hz. The parameters of the filter are shown in table 1,

table 1: butterworth filter coefficient

The phase due to thermal noise can be estimated using the carrier-to-noise ratio, the loop bandwidth and the integration time, and is given by the following formula,

in the formula B_nIs the loop bandwidth, C/N, of the receiver₀Representing the carrier-to-noise ratio, T is the loop integration time; sigma_pllPhase due to thermal noise; the error caused by the thermal noise can be corrected by using the estimated noise jitter and combining the following formula:

is the ionospheric phase flicker index corrected for errors.

As shown in fig. 1, embodiment 1 of the present invention provides a method for calculating an ionospheric phase scintillation index in real time, including:

step 1) obtaining a measured value of a carrier phase from a receiver;

step 2) carrying out differential operation on the carrier phase to obtain a Doppler frequency measurement value of the signal;

step 3) calculating the observed quantity of polynomial fitting according to the time recorded by the Doppler frequency and the Doppler frequency;

the specific calculation method is as follows;

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

wherein T (N) is a matrix with a variable at the power of time;

p (N) is the coefficient of the polynomial;

c (N) is the measured carrier Doppler;

n represents the number of data used to perform the fit.

The time matrix T (N) is:

where t (i) represents the time corresponding to the ith sampling point,

the carrier measurement vector c (n) is:

carph (i) represents the sampled value of the ith carrier phase doppler;

the polynomial coefficient vector P (N) is:

P_irepresents the ith coefficient of the polynomial;

the polynomial coefficients are solved by using the least square method, firstly, the curve equation is normalized, namely coordinates on two sides of the equation are multiplied by the transposition of T (N),

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

wherein U (N) is a group of,

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

wherein D (N) is a group of,

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

since the receiver resources are limited and all past values cannot be stored, the recursion formulas of u (n) and d (n) need to be found, and u (n) and d (n) need to be found in an iterative manner.

The increments of U (N) to U (N +1) are:

the recurrence formula of u (n) is,

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

the increments of D (N) to D (N +1) are:

the recurrence formula for D (N) is therefore:

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

step 4) solving coefficients of the polynomial by using a least square method and calculating the fitted carrier Doppler size according to the current time;

step 5) subtracting the fitted carrier Doppler and the measured carrier Doppler to obtain a carrier frequency value without the movement trend;

step 6) integrating the carrier frequency values without the motion tendency to obtain carrier phases;

step 7) filtering the carrier phase by using a six-order Butterworth high-pass filter, and removing satellite clock error, receiver clock error, flow error and fitted residual error;

and 8) calculating the ionospheric phase flicker index by combining the thermal noise correction value.