A quasi-real-time monitoring method for ionospheric TEC at a single station

Abstract

The invention belongs to the technical field of ionosphere monitoring, and in particular relates to a quasi-real-time monitoring method for ionosphere TEC of a single station. The method is continuously performed based on a GNSS dual-frequency receiver deployed at an ionosphere monitoring station, and comprises: receiving GNSS dual-frequency pseudo-range observations, carrier phase observations, receiver position information, satellite pitch angles, azimuth angles and satellite ephemeris data according to a time length set in a time window; pre-processing the received data to obtain oblique TEC of carrier phase smoothed pseudo-ranges in the line of sight of a satellite receiver; establishing a nonlinear model of vertical TEC directly above the ionosphere monitoring station; and according to the oblique TEC value of the carrier phase smoothed pseudo-range of each satellite, using a Kalman filter iterative solution in a time window sliding manner, obtaining the nonlinear model coefficient of the vertical TEC above the measuring station and the hardware deviation of the receiver and the satellite, thereby realizing quasi-real-time monitoring of the vertical TEC value above the measuring station.

Description

Technical Field

[0001] The invention belongs to the technical field of ionosphere monitoring, and in particular relates to a quasi-real-time monitoring method of single-station ionosphere TEC. Background Art

[0002] The upper atmosphere of the Earth absorbs the sun's extreme ultraviolet rays and X-rays and becomes partially ionized, forming the ionosphere about 60 to 1,000 kilometers above the ground. The typical value of the electron density in the ionosphere is three orders of magnitude smaller than the density of the neutral atmosphere, but because it is sufficient to affect the propagation of radio waves involved in radio astronomy, radar, satellite navigation and communication, shortwave communication and other systems, the ionosphere is an important component of the space (electromagnetic) environment and an important forecast object for space weather.

[0003] The ionosphere has regional characteristics, daily changes, seasonal changes, annual changes and fluctuation structures, which are mainly affected by many factors such as the sun, the geomagnetic field, and the neutral atmosphere. Therefore, there are complex ionospheric disturbances, which will affect the performance of ground-to-air radio systems whose signals pass through or are reflected by the ionosphere. It is very important to carry out ionospheric monitoring and then implement space weather warnings and forecasts.

[0004] For global navigation satellite system (GNSS) positioning, the ionosphere is the largest source of error. By using the influence of the ionosphere on GNSS measurements and inverting the ionosphere total electron content (TEC), an ionosphere correction model can be constructed and applied to GNSS to improve positioning accuracy. On the other hand, using GNSS receiver observations to conduct all-day, all-weather, wide-area ionosphere TEC monitoring has become an important means of studying the ionosphere and conducting ionospheric disturbance warning and forecasting.

[0005] GNSS dual-frequency signals are usually used to measure and invert the ionospheric total electron content (TEC), which can be used to conduct large-scale ionospheric observations and research around the clock. However, at present, the GNSS dual-frequency signal measurement data is generally measured in units of days and processed post-processed by computers, which makes it impossible to achieve real-time monitoring of the ionosphere.

[0006] There are two common ways to monitor ionospheric TEC: network monitoring and single-station monitoring. For active or abnormal ionospheric areas, especially remote areas and oceans, due to geographical restrictions and data transmission restrictions, it is usually necessary to carry out real-time monitoring of the ionosphere at a single station. Single-station monitoring is a flexible and convenient means. For areas where networking is difficult, single-station monitoring is the only option. However, whether it is solved jointly with data from other stations or solved by a single station, ionospheric TEC monitoring is highly dependent on algorithms. The existing technology usually uses post-processing to calculate the hardware deviation of GNSS satellites through the measurement data of the previous day, which requires a large amount of data and is not real-time enough. For scenarios with limited network bandwidth, such as offshore stations, the data received from dozens of satellites is transmitted back in real time, and the amount of data is very large. If the TEC data above the station at one moment is counted as one copy, the corresponding navigation ephemeris, pseudorange and carrier phase measurement data will be as much as two hundred times. For single-station monitoring, it is urgent to develop a quasi-real-time method for GNSS ionospheric TEC. Summary of the invention

[0007] The purpose of the present invention is to overcome the defects of the prior art and propose a quasi-real-time monitoring method of ionospheric TEC at a single station.

[0008] In order to achieve the above object, the present invention proposes a quasi-real-time monitoring method of ionospheric TEC at a single station, which is continuously performed based on a GNSS dual-frequency receiver deployed at an ionospheric monitoring station, comprising:

[0009] Step 1) receiving GNSS dual-frequency pseudo-range observations, carrier phase observations, receiver position information, satellite elevation angles, azimuth angles and satellite ephemeris data according to the duration set in the time window;

[0010] Step 2) preprocessing the received data to obtain the oblique TEC of the carrier phase smoothed pseudorange in the line of sight direction of the satellite receiver;

[0011] Step 3) Establish a nonlinear model of the vertical TEC directly above the ionospheric monitoring station;

[0012] Step 4) According to the oblique TEC value of the carrier phase smoothed pseudorange of each satellite, the Kalman filter is iteratively solved by sliding the time window to obtain the nonlinear model coefficient of the vertical TEC above the station and the hardware deviation of the receiver and the satellite, thereby realizing quasi-real-time monitoring of the vertical TEC value above the station.

[0013] Preferably, the step 1) comprises:

[0014] Receive and decode the ranging information of dual-frequency pseudorange and carrier phase measurements of multiple GNSS satellites according to the duration set in the time window;

[0015] Receive the elevation and azimuth information of each satellite;

[0016] The number, observation time, pitch angle, azimuth angle, GNSS dual-frequency pseudorange and ranging information of carrier phase measurement of each satellite are cached separately.

[0017] Preferably, the step 2) comprises:

[0018] Step 2-1) Obtain the oblique TEC value of the jth satellite based on the dual-frequency pseudorange observation according to the data within the cached time window setting time length and the slant TEC value of the jth satellite based on dual-frequency carrier phase observations

[0019]

[0020] Among them, f1 and f2 are the two frequencies of GNSS satellite signals. are the code pseudorange observations of the jth satellite at two frequencies f1 and f2, are the carrier phase observations of the jth satellite at two frequencies f1 and f2, respectively, and λ1 and λ2 are the wavelengths corresponding to the two frequencies f1 and f2, respectively;

[0021] Step 2-2) Calculate the jth satellite according to the following formula: and Baseline value between

[0022]

[0023] Among them, i represents the i-th data measurement moment in the data received during the time window setting duration, and Num represents the total number of moments corresponding to the data during the time window setting duration; Indicates the elevation angle of the jth satellite relative to the receiver at the i-th data measurement moment;

[0024] Step 2-3) Calculate the oblique TEC value TEC of the jth satellite carrier phase after smoothing pseudorange according to the following formula: s,sys j :

[0025]

[0026] Wherein, the subscript s represents the satellite, and the subscript sys represents the GNSS satellite system, where sys=1 is the GPS system signal, and sys=2 is the Beidou system signal;

[0027] Step 2-4) Calculate the vertical TEC value TEC directly above the ionospheric monitoring station according to the following formula:

[0028]

[0029] Among them, b r,sys To correspond to the receiver hardware bias of each GNSS satellite system, is the satellite hardware bias of the jth satellite; j is the zenith angle of the line-of-sight path between the jth satellite and the receiver at the ionospheric piercing point.

[0030] Preferably, the nonlinear model of step 3) is a spherical harmonic function model, which satisfies the following formula:

[0031]

[0032] Where TEC represents the vertical TEC value directly above the ionospheric monitoring station, M and N represent the order and series of the spherical harmonic function, m and n represent the mth order and nth order, respectively, and θ j is the geographical latitude of the point where the jth satellite signal penetrates the ionosphere, φ j is the solar hour angle of the jth satellite signal puncture point, is the normalized Legendre function, A nm and B nm are the spherical harmonic coefficients to be determined.

[0033] Preferably, the step 4) comprises:

[0034] Step 4-1) construct the state equation and measurement equation of Kalman filtering;

[0035] Step 4-2) Calculate the zenith angle x of the jth satellite signal at the ionosphere penetration point j , latitude θ j and longitude φ j ;

[0036] Step 4-3) Determine whether it is the first Kalman filter calculation. If it is, assign an initial value of 0 to the estimated parameter x0 and an initial value of P0 to the cross-covariance of the estimated parameters. If it is not, inherit the previous estimated parameter x1. k-1 and the cross-covariance P of the estimated parameters k-1 ;

[0037] Step 4-4) Calculate the estimated value of the parameter to be estimated and the covariance matrix P k ;

[0038] Step 4-5) Cache the estimated values ​​of the parameters to be estimated and the covariance matrix P k ; and P k If the hardware bias of satellite j is relatively and Pk-1 middle If there is a change, it is updated; if the jth satellite appears for the first time, then and P k Add the hardware bias of satellite j in variable;

[0039] Step 4-6) The time sliding window slides backward by a fixed length, and according to the next set of preprocessed observations, step 4-2) is re-entered to start the next round of Kalman filter iteration.

[0040] Preferably, the state equation of the Kalman filter in step 4-1) is:

[0041] x k =Φ k,k-1 x k-1 +W k

[0042] Among them, Φ k,k-1 is the state transition matrix from the k-1th moment to the kth moment, W k is the system noise matrix at time k, x k is the parameter vector to be estimated at the kth moment, satisfying the following formula:

[0043]

[0044] Wherein, subscript k represents the time; A 00,k ,A 10,k ,…A 33,k ,B 33,k is the third-order spherical harmonic coefficient to be calculated at the kth moment; b r，1,k is the receiver hardware deviation corresponding to the GPS signal at the kth moment; b r，2,k is the receiver hardware deviation corresponding to the Beidou system signal at the kth moment; represents the hardware bias of the satellite numbered 1, 2, …, J at the kth time. The matrix superscript T represents the matrix transpose.

[0045] Preferably, the measurement equation of Kalman filtering in step 4-1) is:

[0046] z k =H k x k +V k

[0047] Among them, z k is the oblique TEC measurement value after carrier phase smoothing pseudorange at time k The set of , j∈[1,J],

[0048]

[0049] H k is the measurement matrix, satisfying the following formula:

[0050]

[0051] V k is the satellite's The measurement noise.

[0052] Preferably, the step 4-4) comprises:

[0053] Step 4-4-1) Fusion of the predicted values ​​at time k and the measured value z k , and obtain the result of the parameter to be estimated at time k

[0054]

[0055] in, To use the state equation to predict the estimated value of the parameter vector at time k, the following equation is satisfied:

[0056]

[0057] in, is the Kalman filter output result of the parameter vector to be estimated at time k-1;

[0058] K k is the Kalman gain at time k, satisfying the following formula:

[0059]

[0060] R k = Cov(V k )

[0061]

[0062] Among them, R k To measure the noise V k The covariance matrix of is the estimated error covariance matrix of the estimated parameter at time k, P k-1 is the error covariance matrix of the estimated parameter at time k-1, Q k is the system noise W at time k k The covariance matrix satisfies the following formula:

[0063] Q k = Cov(W k )

[0064] Step 4-4-2) Update the error covariance matrix P at time k according to the following formula: k, used to estimate the error covariance matrix of the estimated parameter at time k+1:

[0065]

[0066] Where I is the identity matrix.

[0067] Compared with the prior art, the advantages of the present invention are:

[0068] 1. The quasi-real-time monitoring method of the single-station ionospheric TEC proposed in the present invention is based on the GNSS dual-frequency receiver deployed at the ionospheric monitoring station to continuously solve the ionospheric monitoring problem in areas where networking and data transmission are difficult;

[0069] 2. The method of the present invention adopts a sliding time window data preprocessing method, which greatly reduces the demand for data volume and well solves the problem of limited equipment development environment and computing resources;

[0070] 3. The nonlinear modeling proposed in the present invention and the quasi-real-time TEC inversion algorithm for a single station through Klaman filtering make it possible to monitor the vertical TEC above the measuring station with a relatively small amount of data and achieve monitoring accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 This is a flow chart of the sliding time window data preprocessing method in the present invention;

[0072] Figure 2 This is a single-station Kalman filter TEC inversion calculation flow chart in the present invention;

[0073] Figure 3 It is the output result of the method of the present invention. DETAILED DESCRIPTION

[0074] The technical solution of the present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0075] Example 1

[0076] Embodiment 1 of the present invention discloses a quasi-real-time monitoring method of ionospheric TEC at a single station, which is continuously performed based on a GNSS dual-frequency receiver deployed at an ionospheric monitoring station, and includes:

[0077] Step 1) receiving GNSS dual-frequency pseudo-range observations, carrier phase observations, receiver position information, satellite elevation angles, azimuth angles and satellite ephemeris data according to the duration set in the time window;

[0078] Step 2) preprocessing the received data to obtain the oblique TEC of the carrier phase smoothed pseudorange in the line of sight direction of the satellite receiver;

[0079] Step 3) Establish a nonlinear model of the vertical TEC directly above the ionospheric monitoring station;

[0080] Step 4) According to the oblique TEC value of the carrier phase smoothed pseudorange of each satellite, the Kalman filter is iteratively solved by sliding the time window to obtain the nonlinear model coefficient of the vertical TEC above the station and the hardware deviation of the receiver and the satellite, thereby realizing quasi-real-time monitoring of the vertical TEC value above the station.

[0081] like Figure 1 As shown, the sliding time window data preprocessing method corresponding to step 2) in the present invention includes the following steps:

[0082] Step 2-1, receiving the ranging information of the full system dual-frequency pseudorange and carrier phase measurement output by the GNSS dual-frequency high-precision receiver, and decoding it;

[0083] Step 2-2, receiving the elevation and azimuth information of the GNSS satellite output by the GNSS high-precision receiver;

[0084] Step 2-3, cache the observation time, satellite system, satellite number, elevation angle, azimuth angle, code pseudorange 1, code pseudorange 2, carrier phase measurement 1, and carrier phase measurement 2 information data within 20 minutes;

[0085] Step 2-4, pre-process the 20 minutes of data cached in step 3, and use the 20 minutes of observations to calculate the oblique TEC of each satellite based on the dual-frequency pseudorange observations and the oblique TEC based on the dual-frequency phase observations:

[0086]

[0087] In the formula, the superscript j represents the satellite number; is the slant TEC of the jth satellite based on dual-frequency pseudorange observations, is the slant TEC of the jth satellite based on dual-frequency carrier observations; f1 and f2 are the two frequencies of the GNSS satellite signal, are the code pseudorange observations of the jth satellite at two frequencies f1 and f2 respectively; are the carrier phase observations of the jth satellite at the two frequencies f1 and f2 respectively; λ1 and λ2 are the wavelengths corresponding to the two frequencies f1 and f2 respectively.

[0088] Step 2-5, calculate for the jth satellite and Baseline value between

[0089]

[0090] in, is the jth satellite calculated using 20 minutes of data and The superscript j indicates the satellite number. Indicates the sum of 20 minutes of data, i indicates the i-th data measurement time in the 20 minutes of data, and Num indicates the total number of time points corresponding to the 20 minutes of data; Represents the i-th data measurement time Represents the i-th data measurement time Indicates the elevation angle of the jth satellite relative to the receiver at the i-th data measurement time in the 20-minute data;

[0091] It should be noted that 20 minutes is an embodiment of the present application. If cache conditions permit, data of several hours or longer can also be used as the TEC calculation time window. The TEC calculation accuracy may be better than the 20-minute data used in the present application.

[0092] Step 2-6, calculate the oblique TEC after carrier phase smoothing pseudorange, denoted as TEC s,sys j :

[0093]

[0094] Among them, the superscript j represents the satellite number, the subscript sys represents the GNSS satellite system, sys=1 is the GPS system, and sys=2 is the Beidou system.

[0095] Among them, the oblique TEC (TEC) after the j-th satellite carrier phase smoothing pseudorange obtained in step 2-6 is s,sys j ) contains the hardware bias of the receiver and the satellite. TEC s,sys j The vertical TEC (denoted as TEC) at the ionosphere penetration point of the satellite signal has the following relationship:

[0096]

[0097] That can also be expressed as:

[0098]

[0099] Among them, b r,sys is the receiver hardware bias corresponding to each GNSS satellite system, such as b r,1 is the receiver hardware bias of the GPS system, b r,2 is the hardware deviation of the Beidou system receiver. is the satellite hardware bias of the jth satellite;j is the zenith angle of the line-of-sight path between the jth satellite and the receiver at the ionospheric piercing point.

[0100] For step 3), in low-latitude areas, especially the Equatorial Ionization Anomaly (EIA), the ionospheric vertical TEC has the characteristic of nonlinear variation with latitude, forming a TEC maximum near 20 degrees north latitude (the northern hemisphere hump of EIA), and the TEC decreases rapidly in the area north of the hump. Therefore, using nonlinear functions to model the ionospheric vertical TEC can obtain more accurate TEC inversion calculation results. When the nonlinear function of TEC modeling uses spherical harmonic functions, the specific expression of the ionospheric vertical TEC above the station is as follows:

[0101]

[0102] Where TEC represents the vertical TEC value directly above the ionospheric monitoring station, M and N represent the order and series of the spherical harmonic function, m and n represent the mth order and nth order, respectively, and θ j is the geographical latitude of the point where the jth satellite signal penetrates the ionosphere, φ j is the solar hour angle of the jth satellite signal puncture point, is the normalized Legendre function, A nm and B nm are the spherical harmonic coefficients to be determined.

[0103] It should be noted that this embodiment takes spherical harmonics as an example, but is not limited thereto, and may also be a polynomial function, etc. For spherical harmonics, when a third-order spherical harmonic function is used, M=N=3, but it is not limited to the third order, and may also be other orders.

[0104] A nm and B nm It can be calculated using the Kalman filter method. The number of spherical harmonic coefficients is:

[0105] N SHC =(N+1) 2 -(NM)(N-M+1) (7)

[0106] In the Kalman filter calculation process, define the estimated parameter vector x at the kth moment k for:

[0107]

[0108] Wherein, subscript k represents the time; A 00,k ,A 10,k ,…A 33,k ,B 33,kis the third-order spherical harmonic coefficient to be calculated at the kth moment; b r，1,k is the receiver hardware deviation corresponding to the GPS signal at the kth moment; b r，2,k is the receiver hardware deviation corresponding to the Beidou system signal at the kth moment; represents the hardware bias of the satellite numbered 1, 2, ..., J at the kth time. Kalman filter calculation x k The state equation is:

[0109] x k =Φ k,k-1 x k-1 +W k (9)

[0110] Among them, Φ k,k-1 is the state transfer matrix from the k-1th moment to the kth moment, as shown below:

[0111]

[0112] W k is the system noise matrix at time k:

[0113]

[0114] W k Each value in is the corresponding parameter vector x to be estimated k The noise variance at time k.

[0115] According to formulas (5) and (6), Kalman filtering calculates x k The measurement equation in the process is expressed as:

[0116] z k =H k x k +V k (12)

[0117] Among them, z k is the oblique TEC measurement value after carrier phase smoothing pseudorange at time k According to formula (5):

[0118]

[0119] According to formulas (5) and (6), the H in the measurement equation k as follows:

[0120]

[0121] Where V k is the satellite's The measurement noise.

[0122] Assume that the measurement noise V k and system noise W k are all zero-mean white noise sequences. k The covariance matrix of is taken as:

[0123] R k = Cov(V k )=0.01I J (15)

[0124] Among them, I J is the J×J identity matrix. k The covariance matrix of is taken as:

[0125]

[0126] in, is(N SHC +2+J)×(N SHC +2+J); is the total number of spherical harmonic coefficients of the Mth and Nth order spherical harmonic functions, refer to formula (7); 2 represents the use of satellite data of both GPS and BeiDou systems; J is the total number of satellites used.

[0127] The recursive process of Kalman filter calculation in step 4) can be divided into two steps. The first step is prediction, and the second step is to correct the prediction result using the measurement information.

[0128] The first step of prediction includes using the state equation to predict the estimated value of the parameter vector at time k, denoted as

[0129]

[0130] Among them, Φ k,k-1 is the state transfer matrix, refer to formula (9), (10); is the Kalman filter output result of the parameter vector to be estimated at time k-1. It also includes the estimation error covariance matrix at time k:

[0131]

[0132] P in the formula k-1 is the estimated error covariance matrix at time k-1, Q k Refer to formula (16).

[0133] The second step uses the measurement information to correct the prediction results, including using the measurement matrix, the measurement covariance matrix and equation (18) to define the Kalman gain K at time k k :

[0134]

[0135] In the formula, H k Refer to formula (14), R k Refer to formula (15).

[0136] Fusion k-time prediction value and the measured value z k , and obtain the result of the parameter to be estimated at time k:

[0137]

[0138] Update the estimated error covariance matrix for the next estimation, that is:

[0139]

[0140] Where I is the unit matrix.

[0141] If it is the first Kalman filter calculation, the estimated parameters are assigned an initial value of 0; or if the statistical characteristics of the initial state x0 are known, the estimated parameters are assigned the following initial values:

[0142]

[0143] Using equations (17) to (22), the Kalman filter process can converge quickly to estimate the vertical TEC and hardware bias.

[0144] like Figure 2 As shown, the GNSS satellite hardware bias information inheritance and Kalman filtering algorithm in step 4) of the present invention includes the following steps:

[0145] Step 4-1, construct the state equation and measurement equation of Kalman filtering;

[0146] Step 4-2, calculate the zenith angle χ at the ionospheric puncture point j , latitude θ j and longitude φ j ,j=1,2,…,J is the satellite number at the current epoch;

[0147] Step 4-3, determine whether it is the first Kalman filter calculation;

[0148] Step 4-4, if yes, then assign an initial value of 0 to the estimated parameter x0 and an initial value of the cross-covariance P0 of the estimated parameters; otherwise, inherit the previous estimated parameter x k-1 and the cross-covariance P of the estimated parameters k-1 ;

[0149] Step 4-5, calculate the covariance matrix R of measurement noise and system noise k and Q k ;

[0150] Step 4-6, calculate the state transfer matrix Φ k,k-1 and the measurement matrix H k ;

[0151] Step 4-7, calculate the covariance matrix of the prediction error and the Kalman filter gain matrix K k ;

[0152] Step 4-8, calculate the estimated value of the parameter to be estimated and the covariance matrix P k ;

[0153] Steps 4-9, cache estimates and the covariance matrix P k ; and P k If the hardware bias of satellite j is relatively and P k-1 middle If there is a change, it is updated; if the jth satellite appears for the first time, then and P k Add the hardware bias of satellite j in variable;

[0154] In step 4-10, when the time sliding window slides back 5 minutes, step 4-2 is re-entered to start the next round of Kalman filter iteration based on the next set of preprocessed observations.

[0155] Due to the movement of GNSS satellites, the received GNSS satellite signals may fade out of the receiver's reception range, and some satellites may re-enter the receiver's reception range. Satellite hardware bias in the parameters to be estimated The number of x varies with the number of visible satellites and is not a fixed vector. The estimated parameter x in the Kalman filtering process k And the intermediate calculation matrix must also change with the change of the observed quantity.

[0156] Parameters to be estimated is the hardware bias of the satellite, which generally changes very little. After the satellite fades out and then enters the observation area, the hardware bias of the satellite The result retained before fading out is inherited as the initial value. Correspondingly, the cross-correlation matrix of the observed quantity also inherits the result retained before fading out as the initial value of the satellite when it enters this time, thereby avoiding the instability problem caused by multiple initializations in the Kalman filtering process and ensuring the reliability of the results.

[0157] In one embodiment, an ionosphere monitoring device is provided, which is deployed in a GNSS dual-frequency receiver of an ionosphere monitoring station and is implemented based on the above method. The device includes:

[0158] The data cache module is used to cache GNSS dual-frequency pseudo-range observations, carrier phase observations, receiver position information, satellite pitch angles, azimuth angles and satellite ephemeris data according to the duration set in the time window, and transmit them to the TEC inversion calculation module and the storage module respectively; in one embodiment, we set the time window to 20 minutes.

[0159] The TEC inversion calculation module is implemented based on an embedded processor and is used to pre-process the received data using a sliding time window, establish a spherical harmonic function model, and use Kalman filtering to iteratively solve the nonlinear model coefficients and the hardware deviations of the receiver and the satellite, thereby realizing quasi-real-time monitoring of the TEC value and transmitting the TEC value to the storage module. In one embodiment, we set the time sliding window to slide back for 5 minutes to perform the next round of Kalman filtering iteration.

[0160] The storage module is used to store the GNSS measurement data received by the GNSS high-precision OEM board, and is also used to store the TEC value output by the TEC inversion calculation module.

[0161] like Figure 3 The results of the Klaman filter single-station quasi-real-time TEC inversion algorithm for spherical harmonics modeling TEC implemented by the present invention are shown, and compared with the results of the TEC reference value and the Klaman filter single-station algorithm for linear function modeling TEC. The TEC reference value in the figure is the result of inverting TEC using spherical harmonics for a whole day of measurement data, and Klaman+linear is the result of the Klaman filter single-station algorithm for linear function modeling TEC. It can be seen from the result graph that the Klaman+spherical harmonics result using the present invention is closer to the reference result.

[0162] In summary, the single-station ionospheric TEC quasi-real-time monitoring method provided by the present invention uses a sliding time window data preprocessing method, nonlinear function TEC modeling, GNSS satellite hardware bias information inheritance and Kalman filtering algorithm to achieve quasi-real-time single-station inversion calculation of ionospheric TEC.

[0163] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention is described in detail with reference to the embodiments, it should be understood by those skilled in the art that any modification or equivalent replacement of the technical solutions of the present invention does not depart from the spirit and scope of the technical solutions of the present invention and should be included in the scope of the claims of the present invention.

Claims

1. A quasi-real-time monitoring method for ionospheric TEC at a single station, based on a GNSS dual-frequency receiver deployed at an ionospheric monitoring station, comprising: Step 1) receiving GNSS dual-frequency pseudo-range observations, carrier phase observations, receiver position information, satellite elevation angles, azimuth angles and satellite ephemeris data according to the duration set in the time window; Step 2) preprocessing the received data to obtain the oblique TEC of the carrier phase smoothed pseudorange in the direction of the satellite receiver's line of sight; Step 3) Establish a nonlinear model of the vertical TEC directly above the ionospheric monitoring station; Step 4) According to the oblique TEC value of the carrier phase smoothed pseudorange of each satellite, the nonlinear model coefficient of the vertical TEC above the measuring station and the hardware deviation of the receiver and the satellite are obtained by using the Kalman filter iteration solution through the time window sliding method, thereby realizing the quasi-real-time monitoring of the vertical TEC value above the measuring station; The nonlinear model of step 3) is a spherical harmonic function model, which satisfies the following formula: Where TEC represents the vertical TEC value directly above the ionospheric monitoring station, M and N represent the order and series of the spherical harmonic function, m and n represent the mth order and nth order, respectively, and θ j is the geographical latitude of the point where the jth satellite signal penetrates the ionosphere, φ j is the solar hour angle of the jth satellite signal puncture point, is the normalized Legendre function, A nm and B nm are the spherical harmonic coefficients to be determined; The step 4) comprises: Step 4-1) construct the state equation and measurement equation of Kalman filtering; Step 4-2) Calculate the zenith angle x of the jth satellite signal at the ionosphere penetration point j , latitude θ j and longitude φ j ; Step 4-3) Determine whether it is the first Kalman filter calculation. If it is, assign an initial value of 0 to the estimated parameter x0 and an initial value of P0 to the cross-covariance of the estimated parameters. If it is not, inherit the previous estimated parameter x1. k-1 and the cross-covariance P of the estimated parameters k-1 ; Step 4-4) Calculate the estimated value of the parameter to be estimated and the covariance matrix P k ; Step 4-5) Cache the estimated values ​​of the parameters to be estimated and the covariance matrix P k ; and P k If the hardware bias of satellite j is relatively and P k-1 middle If there is a change, it is updated; if the jth satellite appears for the first time, then and P k Add the hardware bias of satellite j in variable; Step 4-6) The time sliding window slides backward by a fixed length, and according to the next set of preprocessed observations, step 4-2) is re-entered to start the next round of Kalman filter iteration.

2. The quasi-real-time monitoring method of single-station ionospheric TEC according to claim 1, characterized in that: The step 1) comprises: Receive and decode the ranging information of dual-frequency pseudorange and carrier phase measurements of multiple GNSS satellites according to the duration set in the time window; Receive the elevation and azimuth information of each satellite; The number, observation time, pitch angle, azimuth angle, GNSS dual-frequency pseudorange and ranging information of carrier phase measurement of each satellite are cached separately.

3. The quasi-real-time monitoring method of single-station ionospheric TEC according to claim 1, characterized in that: The step 2) comprises: Step 2-1) Obtain the oblique TEC value of the jth satellite based on the dual-frequency pseudorange observation according to the data within the cached time window setting time length and the slant TEC value of the jth satellite based on dual-frequency carrier phase measurement observations Among them, f1 and f2 are the two frequencies of GNSS satellite signals. are the code pseudorange observations of the jth satellite at two frequencies f1 and f2, are the carrier phase observations of the jth satellite at two frequencies f1 and f2, respectively, and λ1 and λ2 are the wavelengths corresponding to the two frequencies f1 and f2, respectively; Step 2-2) Calculate the jth satellite according to the following formula: and Baseline value between Among them, i represents the i-th data measurement moment in the data received during the time window setting duration, and Num represents the total number of moments corresponding to the data during the time window setting duration; Indicates the elevation angle of the jth satellite relative to the receiver at the i-th data measurement moment; Step 2-3) Calculate the oblique TEC value TEC of the jth satellite after phase smoothing pseudorange according to the following formula: s,sys j : Wherein, the subscript s represents the satellite, and the subscript sys represents the GNSS satellite system, where sys=1 is the GPS system, and sys=2 is the BeiDou system; Step 2-4) Calculate the vertical TEC value TEC directly above the ionospheric monitoring station according to the following formula: Among them, b r,sys To correspond to the receiver hardware bias of each GNSS satellite system, is the satellite hardware bias of the jth satellite; j is the zenith angle of the line-of-sight path between the jth satellite and the receiver at the ionospheric piercing point.

4. The quasi-real-time monitoring method of single-station ionospheric TEC according to claim 1, characterized in that: The state equation of Kalman filtering in step 4-1) is: x k =Φ k,k-1 x k-1 +W k Among them, Φ k,k-1 is the state transition matrix from the k-1th moment to the kth moment, W k is the system noise matrix at time k, x k is the parameter vector to be estimated at the kth moment, satisfying the following formula: Wherein, subscript k represents the time; A 00,k ,A 10,k ,…A 33,k ,B 33,k is the third-order spherical harmonic coefficient to be calculated at the kth moment; b r，1,k is the receiver hardware deviation corresponding to the GPS signal at the kth moment; b r，2,k is the receiver hardware deviation corresponding to the Beidou system signal at the kth moment; represents the hardware bias of the satellite numbered 1, 2, …, J at the kth time. The matrix superscript T represents the matrix transpose.

5. The quasi-real-time monitoring method of single-station ionospheric TEC according to claim 1, characterized in that: The measurement equation of Kalman filtering in step 4-1) is: z k =H k x k +V k Among them, z k is the oblique TEC measurement value after carrier phase smoothing pseudorange at time k A collection of j∈[1,J], H k is the measurement matrix, satisfying the following formula: V k is the satellite's The measurement noise.

6. The quasi-real-time monitoring method of single-station ionospheric TEC according to claim 1, characterized in that: The step 4-4) comprises: Step 4-4-1) Fusion of the predicted values ​​at time k and the measured value z k , and obtain the result of the parameter to be estimated at time k in, To use the state equation to predict the estimated value of the parameter vector at time k, the following equation is satisfied: in, is the Kalman filter output result of the parameter vector to be estimated at time k-1; K k is the Kalman gain at time k, satisfying the following formula: R k =Cov(V k ) Among them, R k To measure the noise V k The covariance matrix of is the estimated error covariance matrix of the estimated parameter at time k, P k-1 is the error covariance matrix of the estimated parameter at time k-1, Q k is the system noise W at time k k The covariance matrix satisfies the following formula: Q k =Cov(W k ) Step 4-4-2) Update the error covariance matrix P at time k according to the following formula: k , used to estimate the error covariance matrix of the estimated parameter at time k+1: Where I is the identity matrix.

Images