GNSS troposphere wet delay estimation method, system, equipment and medium
阅读说明：本技术 一种gnss对流层湿延迟估计方法、系统、设备及介质 (GNSS troposphere wet delay estimation method, system, equipment and medium ) 是由 张兴刚 贺成艳 于 2021-09-14 设计创作，主要内容包括：本发明公开了一种GNSS对流层湿延迟估计方法、系统、设备及介质,所述方法包括以下步骤：获取不同频点的观测值,在每个历元形成无电离层组合观测方程；根据先验ZWD随机游走参数进行Kalman滤波状态更新和观测更新,得到未知变量预报值和其方差协方差阵P-(t|t-1)；计算给定先验ZWD随机游走噪声和观测噪声时的数学期望；以所述数学期望作为伪观测值,应用滑动窗口里的数据估计新的ZWD随机游走参数；用新的ZWD随机游走参数更新先验ZWD随机游走参数。本发明的方法可以实时估计随时间变化的随机模型参数,可改善对流层湿延迟估计值。(The invention discloses a GNSS troposphere wet delay estimation method, a GNSS troposphere wet delay estimation system, GNSS troposphere wet delay estimation equipment and a GNSS troposphere wet delay estimation medium, wherein the GNSS troposphere wet delay estimation method comprises the following steps: obtaining observed values of different frequency points, and forming an ionosphere-free combined observation equation in each epoch; performing Kalman filtering state updating and observation updating according to prior ZWD random walk parameters to obtain an unknown variable prediction value And its variance covariance matrix P t|t‑1 (ii) a Calculating a mathematical expectation for a given prior ZWD random walk noise and observation noise; estimating new ZWD random walk parameters by using data in a sliding window with the mathematical expectation as a pseudo-observed value; the a priori ZWD random walk parameters are updated with the new ZWD random walk parameters. The method can estimate the random model parameters changing along with time in real time and can improve the wet delay of the troposphereAnd (6) estimating the value.)
1. A GNSS tropospheric wet delay estimation method, comprising:
obtaining observed values of different frequency points, and forming an ionosphere-free combined observation equation in each epoch;
performing Kalman filtering state updating and observation updating according to prior ZWD random walk parameters to obtain an unknown variable prediction valueAnd its variance covariance matrix Pt|t-1(ii) a The unknown variables comprise station coordinates, a receiver clock, an ambiguity parameter and a ZWD; calculating to obtain a new parameter filtering value based on the non-ionized layer combined observation equationAnd its variance covariance matrix Pt|tAnd calculating to obtain the after-the-fact residual errorPerforming one backward smoothing to obtain a parameter smoothing estimation valueAnd its variance covariance matrix Pt-1|tAnd anAndcovariance matrix P oft,t-1|t；
According toPt|t、Pt-1|t、Pt,t-1|tAnd post-mortem residualCalculating a mathematical expectation for a given prior ZWD random walk noise and observation noise; taking the mathematical expectation as a pseudo observed value, inserting the mathematical expectation into the tail part of a sliding window queue, maximizing, and estimating new ZWD random walk parameters by using data in the sliding window; deleting doorAnd (4) removing the head data of the sliding window, and updating the prior ZWD random walk parameters by using the new ZWD random walk parameters to finish the GNSS troposphere wet delay estimation.
2. The GNSS troposphere wet delay estimation method of claim 1, wherein the step of obtaining observations at different frequency points and forming an ionosphere-free combined observation equation at each epoch specifically includes:
extracting observed values of different frequency points from the code pseudo range and the carrier observation file, and forming the following non-ionosphere combination in each epoch:
wherein the content of the first and second substances,representing the satellite position; c represents the speed of light;representing the coordinate of the survey station; dtiRepresenting the satellite clock error; dtrRepresenting the receiver clock error;representing a tropospheric mapping function;a tropospheric delay correction value is represented, is the dry delay mapping function, ZHD is the zenith delay statics delay,is the wet delay mapping function, ZWD is the wet zenith delay as previously described; c δ trelIs a correction value for the effects of relativity,whereinWhich represents the constant of the earth's gravity,andwavelength and ambiguity for the ionospheric-free combination;andpseudorange and carrier noise, respectively; deltaPCOAnd deltaPCVRespectively an antenna phase center and a phase change correction value; brcvAnd biReceiver and satellite code hardware delays, respectively; deltaWindupCorrecting for phase winding; the coordinates of the survey station after the correction,is the coordinate of the receiver monument,is a corrected value of the solid tide at the measuring station,is the corrected value of the tide,is the extreme tidal position correction value.
3. The GNSS troposphere wet delay estimation method of claim 1, further comprising, before obtaining observations at different frequency points and forming an ionosphere-free combined observation equation for each epoch, the following steps:
data preparation, comprising: specifying a priori ZWD random walk parameter; downloading code pseudo range and carrier wave observation files, code difference DCB files, precision orbit files, clock error files, planet ephemeris files, tide files, earth orientation parameter files, receiver antenna files, satellites and receiver hardware delay files.
4. The GNSS tropospheric wet delay estimation method of claim 1, characterized by being used for monitoring and forecasting extreme weather events.
5. A GNSS tropospheric wet delay estimation system, comprising:
the equation acquisition module is used for acquiring observed values of different frequency points and forming an ionosphere-free combined observation equation in each epoch;
a parameter obtaining module for performing Kalman filtering state update and observation update according to prior ZWD random walk parameter to obtain unknown variable prediction valueAnd its variance covariance matrix Pt|t-1(ii) a The unknown variables comprise station coordinates, a receiver clock, an ambiguity parameter and a ZWD; calculating to obtain a new parameter filtering value based on the non-ionized layer combined observation equationAnd its variance covariance matrix Pt|tAnd calculating to obtain the after-the-fact residual errorPerforming one backward smoothing to obtain a parameter smoothing estimation valueAnd its variance covariance matrix Pt-1|tAnd anAndcovariance matrix P oft,t-1|t；
A parameter update module for updating the parameters based onPt|t、Pt-1|t、Pt,t-1|tAnd post-mortem residualCalculating a mathematical expectation for a given prior ZWD random walk noise and observation noise; taking the mathematical expectation as a pseudo observed value, inserting the mathematical expectation into the tail part of a sliding window queue, maximizing, and estimating new ZWD random walk parameters by using data in the sliding window; deleting the head of the sliding windowAnd data, updating the prior ZWD random walk parameters by using the new ZWD random walk parameters, and finishing the GNSS troposphere wet delay estimation.
6. The GNSS troposphere wet delay estimation system of claim 5, wherein the equation acquisition module acquires observations at different frequency points, and the step of forming an ionosphere-free combined observation equation at each epoch specifically comprises:
extracting observed values of different frequency points from the code pseudo range and the carrier observation file, and forming the following non-ionosphere combination in each epoch:
wherein the content of the first and second substances,representing the satellite position; c represents the speed of light;representing the coordinate of the survey station; dtiRepresenting the satellite clock error; dtrRepresenting the receiver clock error;representing a tropospheric mapping function;a tropospheric delay correction value is represented, is the dry delay mapping function, ZHD is the zenith delay statics delay,is the wet delay mapping function, ZWD is the wet zenith delay as previously described; c δ trelIs a correction value for the effects of relativity, which represents the constant of the earth's gravity,andwavelength and ambiguity for the ionospheric-free combination;andpseudorange and carrier noise, respectively; deltaPCOAnd deltaPCVRespectively an antenna phase center and a phase change correction value; brcvAnd biReceiver and satellite code hardware delays, respectively; deltaWindupCorrecting for phase winding; the coordinates of the survey station after the correction,is the coordinate of the receiver monument,is a corrected value of the solid tide at the measuring station,is the corrected value of the tide,is the extreme tidal position correction value.
7. The GNSS tropospheric wet delay estimation system of claim 5 further comprising:
the data preparation module is used for appointing prior ZWD random walk parameters; downloading code pseudo range and carrier wave observation files, code difference DCB files, precision orbit files, clock error files, planet ephemeris files, tide files, earth orientation parameter files, receiver antenna files, satellites and receiver hardware delay files.
8. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that the processor when executing the computer program realizes the steps of the GNSS tropospheric wet delay estimation method of any of claims 1 to 3.
9. A computer-readable storage medium having a computer program stored thereon, which, when being executed by a processor, carries out the steps of the GNSS tropospheric wet delay estimation method of any one of claims 1 to 3.
Technical Field
The invention belongs to the technical field of satellite positioning, relates to the field of troposphere wet delay estimation, and particularly relates to a GNSS troposphere wet delay estimation method, a GNSS troposphere wet delay estimation system, GNSS troposphere wet delay estimation equipment and a GNSS troposphere wet delay estimation medium.
Background
The troposphere is one of the main error sources of space geodetic measurement technologies such as a Global Navigation Satellite System (GNSS), Very Long Baseline Interferometry (VLBI) and Doppler orbit determination and radio positioning system (DORIS), and the troposphere has become a ground reference system for realizing millimeter-level high precision at present, reliably measures and explains springback after ice, and monitors restriction factors of climate change influences such as sea level, global warming and the like.
Currently, tropospheric Delay is generally divided into dry Delay (ZHD) and Wet Delay (ZWD).
(1) The dry delay can be accurately modeled by empirical models, which can be divided into two categories: the first type uses ground meteorological data as input, and calculates tropospheric Delay (ZTD) in Zenith direction according to empirical meteorological formulas such as Saastamoinen (1972) and Hopfield (1969); the second category requires time and approximate coordinates to calculate tropospheric delay using averages of parameters such as anniversary, semianniversary, and rate of temperature change for Numerical Weather Models (NWM), including GPT and GPT2, among others.
(2) In contrast to dry delay ZHD, because tropospheric wet delay ZWD depends on atmospheric moisture content, it varies dramatically over time and space, making ZWD difficult to predict and model accurately and can only be estimated as an unknown parameter. Generally, the wet delay of a GNSS signal in the inclination direction is mapped to the zenith direction through a mapping function, and the wet delay ZWD of the zenith, the coordinates of a receiver and clock error are estimated epoch by epoch; mapping functions that map oblique tropospheric delays to zenith directions include Neill (1996) and VMF1, among others. In addition, it is sometimes desirable to use a three-dimensional mapping function that estimates tropospheric zenith delays while simultaneously estimating tropospheric gradients, taking into account the effects of atmospheric azimuthal asymmetry.
With the development of new GNSS data processing techniques, (near) real-time high precision products of GNSS tracks and clocks enable one to estimate ZWD in real-time with a precision of about 1 cm, compared to post-processed ZWD products, and there is still a need to improve the (near) real-time ZWD precision. The currently predominant methods include:
(1) a move from single GPS constellations to multi-GNSS constellations so that more observations and better constellation geometry can be utilized.
(2) The discrete empirical mapping function is refined so that a portion of the irregular variation ZWD enters the empirical model and can perform better at a particular elevation angle. Both of these approaches have been developed relatively well.
(3) The real-time ZWD estimate is improved by a real-time estimation ZWD stochastic model. When processing GNSS data, the ZWD is usually modeled as Random Walk Noise (RWN), and different software currently adopts different fixed RWN values, which range from 0.01 toHowever, this method of fixing RWN sometimes cannot distinguish between normal weather events and extreme weather events, resulting in inaccurate estimation of ZWD under extreme weather conditions, and therefore it is necessary to determine the optimal RWN according to the GNSS station location and local weather conditions.
For the above prior art method (3), RWN can be obtained from the numerical weather model NWM, as can be known from the RWN definition; however, the calculation of RWN by NWM still suffers from the following disadvantages:
(1) because the NWM is inaccurate, and the size and trend of the ZWD depends on the mapping function used and the GPS observation stochastic model, the RWN derived by the NWM may not be consistent with the background model of the software.
(2) The NWM does not guarantee the accuracy of the RWN estimation.
(3) The products of NWM are not publicly available.
(4) Computing the RWN by the NWM requires computing using ray tracing, resulting in complex and time consuming computations.
(5) The NWM is usually issued once every 6 hours, and the time resolution is low; for example, strong convection weather such as storm with a duration of only 2 hours suddenly comes temporarily, and cannot meet the requirement.
Disclosure of Invention
The present invention is directed to a GNSS tropospheric wet delay estimation method, system, device and medium, so as to solve one or more of the above technical problems. The method of the invention is based on Bayes principle and maximum likelihood estimation method, fully considers the geographical position of the GNSS observation station and the time-varying characteristic of rapid change of the troposphere, can estimate the random model parameter changing along with the time in real time, and can improve the estimated value of the troposphere wet delay.
In order to achieve the purpose, the invention adopts the following technical scheme:
the invention provides a GNSS troposphere wet delay estimation method, which comprises the following steps:
obtaining observed values of different frequency points, and forming an ionosphere-free combined observation equation in each epoch;
performing Kalman filtering state updating and observation updating according to prior ZWD random walk parameters to obtain an unknown variable prediction valueAnd its variance covariance matrix Pt|t-1(ii) a The unknown variables comprise station coordinates, a receiver clock, an ambiguity parameter and a ZWD; calculating to obtain a new parameter filtering value based on the non-ionized layer combined observation equationAnd its variance covariance matrix Pt|tAnd calculating to obtain the after-the-fact residual errorPerforming one backward smoothing to obtain a parameter smoothing estimation valueAnd its variance covariance matrix Pt-1|tAnd anAndcovariance matrix P oft,t-1|t；
According toPt|t、Pt-1|t、Pt,t-1|tAnd post-mortem residualCalculating a mathematical expectation for a given prior ZWD random walk noise and observation noise; taking the mathematical expectation as a pseudo observed value, inserting the mathematical expectation into the tail part of a sliding window queue, maximizing, and estimating new ZWD random walk parameters by using data in the sliding window; deleting the data at the head of the sliding window, and updating the prior ZWD random walk parameter by using the new ZWD random walk parameter to finish the GNSS troposphere wet delay estimation.
The method of the invention is further improved in that the step of obtaining the observed values of different frequency points and forming the ionosphere-free combined observation equation in each epoch specifically comprises the following steps:
extracting observed values of different frequency points from the code pseudo range and the carrier observation file, and forming the following non-ionosphere combination in each epoch:
wherein the content of the first and second substances,representing the satellite position; c represents the speed of light;representing the coordinate of the survey station; dtiRepresenting the satellite clock error; dtrRepresenting the receiver clock error;representing a tropospheric mapping function;a tropospheric delay correction value is represented, is the dry delay mapping function, ZHD is the zenith delay statics delay,is the wet delay mapping function, ZWD is the wet zenith delay as previously described; c δ trelIs a correction value for the effects of relativity,whereinWhich represents the constant of the earth's gravity,andwavelength and ambiguity for the ionospheric-free combination;andpseudorange and carrier noise, respectively; deltaPCOAnd deltaPCVRespectively an antenna phase center and a phase change correction value; brcvAnd biReceiver and satellite code hardware delays, respectively; deltaWindupCorrecting for phase winding; the coordinates of the survey station after the correction,is the coordinate of the receiver monument,is a corrected value of the solid tide at the measuring station,is the corrected value of the tide,is the extreme tidal position correction value.
The method is further improved in that before the observation values of different frequency points are obtained and the ionosphere-free combined observation equation is formed in each epoch, the method further comprises the following steps:
data preparation, comprising: specifying a priori ZWD random walk parameter; downloading code pseudo range and carrier wave observation files, code difference DCB files, precision orbit files, clock error files, planet ephemeris files, tide files, earth orientation parameter files, receiver antenna files, satellites and receiver hardware delay files.
A further improvement of the method of the present invention is for monitoring and forecasting extreme weather events.
The invention discloses a GNSS troposphere wet delay estimation system, which comprises:
the equation acquisition module is used for acquiring observed values of different frequency points and forming an ionosphere-free combined observation equation in each epoch;
a parameter obtaining module for performing Kalman filtering state update and observation update according to prior ZWD random walk parameter to obtain unknown variable prediction valueAnd its variance covariance matrix Pt|t-1(ii) a The unknown variables comprise station coordinates, a receiver clock, an ambiguity parameter and a ZWD; calculating to obtain a new parameter filtering value based on the non-ionized layer combined observation equationAnd its variance covariance matrix Pt|tAnd calculating to obtain the after-the-fact residual errorPerforming one backward smoothing to obtain a parameter smoothing estimation valueAnd its variance covariance matrix Pt-1|tAnd anAndcovariance matrix P oft,t-1|t；
A parameter update module for updating the parameters based onPt|t、Pt-1|t、Pt,t-1|tAnd post-mortem residualCalculating a mathematical expectation for a given prior ZWD random walk noise and observation noise; taking the mathematical expectation as a pseudo observed value, inserting the mathematical expectation into the tail part of a sliding window queue, maximizing, and estimating new ZWD random walk parameters by using data in the sliding window; deleting the data at the head of the sliding window, and updating the prior ZWD random walk parameter by using the new ZWD random walk parameter to finish the GNSS troposphere wet delay estimation.
The system of the invention is further improved in that the equation acquisition module acquires observed values of different frequency points, and the step of forming the ionosphere-free combined observation equation in each epoch specifically comprises the following steps:
extracting observed values of different frequency points from the code pseudo range and the carrier observation file, and forming the following non-ionosphere combination in each epoch:
wherein the content of the first and second substances,representing the satellite position; c represents the speed of light;representing the coordinate of the survey station; dtiRepresenting the satellite clock error; dtrRepresenting the receiver clock error;representing a tropospheric mapping function;a tropospheric delay correction value is represented, is the dry delay mapping function, ZHD is the zenith delay statics delay,is the wet delay mapping function, ZWD is the wet zenith delay as previously described; c δ trelIs a correction value for the effects of relativity,andis a non-ionized layer combinationWavelength and ambiguity of;andpseudorange and carrier noise, respectively; deltaPCOAnd deltaPCVRespectively an antenna phase center and a phase change correction value; brcvAnd biReceiver and satellite code hardware delays, respectively; deltaWindupCorrecting for phase winding;the coordinates of the survey station after the correction,is the coordinate of the receiver monument,is a corrected value of the solid tide at the measuring station,is the corrected value of the tide,is the extreme tidal position correction value.
The system of the invention is further improved in that the system further comprises:
the data preparation module is used for appointing prior ZWD random walk parameters; downloading code pseudo range and carrier wave observation files, code difference DCB files, precision orbit files, clock error files, planet ephemeris files, tide files, earth orientation parameter files, receiver antenna files, satellites and receiver hardware delay files.
A computer device of the invention comprises a memory, a processor and a computer program stored in the memory and executable on the processor, the processor implementing the steps of any of the above mentioned GNSS tropospheric wet delay estimation methods of the invention when executing the computer program.
A computer-readable storage medium of the invention stores a computer program which, when executed by a processor, implements the steps of any of the above-described GNSS tropospheric wet delay estimation methods of the invention.
Compared with the prior art, the invention has the following beneficial effects:
when extreme weather occurs, different random walk noise constraints and GNSS observation weights can cause that GNSS troposphere estimation values have great difference, and the current method for setting GNSS troposphere random walk parameters according to experience cannot meet the troposphere estimation requirements under the extreme weather. The method disclosed by the invention is independent of a numerical weather model, adopts a variance component estimation algorithm suitable for estimating the GNSS troposphere random walk noise in real time, and can provide a more mature and reliable GNSS troposphere product for monitoring extreme weather and climate change. Specifically, the method takes the Bayesian principle and the maximum likelihood estimation method as the basis, fully considers the geographical position of the GNSS observation station and the time-varying characteristic of rapid change of the troposphere, can estimate the random model parameters changing along with time in real time, can improve the troposphere zenith delay estimation value, and is beneficial to monitoring and tracking extreme weather events.
Further explained, the method of the present invention is independent of the numerical weather model, and does not require purchasing real-time numerical weather products; compared with the traditional ray tracing algorithm depending on a numerical weather model, the method is easy to implement and high in calculation efficiency. The method has the advantages that the characteristic that GNSS observation noise and multipath depend on altitude angles is taken into consideration, the geographical position difference of a GNSS observation station and the time-varying characteristic of rapid change of a troposphere are fully considered, the optimal random model parameter changing along with time can be estimated in real time, the ZWD can be improved (exemplarily, the ZWD can be improved by about 11%), the local peak value of the ZWD is more obvious, and severe weather events can be monitored and forecasted more reliably.
Further, when the random walk noise is directly estimated by applying the conventional variance component estimation technique, convergence is difficult and only post-processing is possible. The method carries out reasonable algorithm assumption on the basis of the expected maximization algorithm, restrains the phase variance component of the GNSS carrier, ensures the convergence of the algorithm and can carry out real-time estimation.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art are briefly introduced below; it is obvious that the drawings in the following description are some embodiments of the invention, and that for a person skilled in the art, other drawings can be derived from them without inventive effort.
Fig. 1 is a flowchart illustrating a GNSS tropospheric wet delay estimation method according to an embodiment of the present invention;
FIG. 2 is a diagram illustrating a conditional mathematical expectation sliding window queue for random noise in an embodiment of the present invention;
FIG. 3 is an IGS test station distribution diagram for performing algorithm value test validation, in accordance with an embodiment of the present invention;
FIG. 4 is a schematic diagram of a random walk noise time sequence of each IGS station estimated by applying the method of the present invention in the embodiment of the present invention;
FIG. 5 is a schematic view of the troposphere ZTD time series in extreme weather in the embodiment of the invention; where ZWD is converted to ZTD for comparison with JPL post hoc tropospheric products.
Detailed Description
In order to make the purpose, technical effect and technical solution of the embodiments of the present invention clearer, the following clearly and completely describes the technical solution of the embodiments of the present invention with reference to the drawings in the embodiments of the present invention; it is to be understood that the described embodiments are only some of the embodiments of the present invention. Other embodiments, which can be derived by one of ordinary skill in the art from the disclosed embodiments without inventive faculty, are intended to be within the scope of the invention.
Referring to fig. 1, the improved GNSS troposphere wet delay estimation method according to the embodiment of the present invention can be decomposed into the following three steps when estimating ZWD per epoch, wherein Step1 and Step2 are conventional ZWD estimation algorithms; step3 is an algorithm module newly added in the embodiment of the invention, and is an innovative core of the method in the embodiment of the invention.
Step1, data preparation, including specifying prior random parameters, download code pseudo range and carrier observation files and inter-code difference DCB files, precision orbit and clock error files, planet ephemeris files and tide files, earth orientation parameter files, receiver antenna files, satellite and receiver hardware delay files. After the files are prepared, the next step is carried out to extract the observed value and correct the observed value.
Step 2. precision point-and-point positioning (PPP) estimation ZWD. Hereinafter with GPSL1And L2The observed values of the frequency points are briefly introduced for example.
Firstly, extracting observed values of different frequency points from a code pseudo range and a carrier observed file, and forming the following non-ionosphere combination in each epoch:
wherein the content of the first and second substances,the satellite position can be obtained by interpolation calculation of a precise orbit file; c represents the speed of light;is a coordinate of a survey station and can be directly extracted by SINEX; dtiThe satellite clock error can be obtained by calculating a precise satellite clock error file; dtrIs the receiver clock error;is a function of the tropospheric mapping,is tropospheric delayCorrection values, divided into dry retardation and wet retardation:
is the dry delay mapping function, ZHD is the zenith delay statics delay,is the wet delay mapping function and the ZWD is the wet zenith delay as previously described.
cδtrelIs a relativistic effect correction value and can be calculated by a receiver, satellite positions and velocities:
wherein the content of the first and second substances,is the constant of earth's gravity.
Andwavelength and ambiguity for the ionospheric-free combination;andthe pseudorange and the noise of the carrier wave, respectively.
δPCOAnd deltaPCVRespectively are an antenna phase center and a phase change correction value, and can be calculated and corrected by an antenna file; brcvAnd biAre respectively provided withThe hardware delay of the receiver and the satellite code can be extracted and corrected by a hardware delay file; deltaWindupFor phase winding correction, it can be calculated from the satellite position and the handset position, and the specific calculation method is Wu et al (1993).
When the observation equations (1) and (2) are formed, a precise orbit and a clock difference file are required to be applied to preprocess data, wherein the preprocessing comprises cycle slip detection, receiver millisecond slip repair, cycle slip detection and receiver millisecond slip repair. Meanwhile, ephemeris files, tide files, earth orientation parameters and antenna files are required to be applied to correct the GPS position:
wherein the content of the first and second substances,the coordinates of the survey station after the correction,is the coordinate of the receiver monument,is a corrected value of the solid tide at the measuring station,is the corrected value of the tide,is the extreme tidal position correction value. Other modifications such as calculating phase center corrections, phase unwrapping corrections, relativistic effect corrections, prior tropospheric corrections, etc. may refer to the IERS2020 protocol and are not described in detail herein.
Then, according to the prior troposphere random parameter information, Kalman filtering state updating and observation updating are carried out to obtain new prediction values of unknown variables (including station coordinates, a receiver clock, an ambiguity parameter and a ZWD)Covariance matrix P with predicted valuest|t-1The subscript represents the estimation quantity of the epoch t moment obtained from the observation value (including the t-1 moment) in front of the epoch t-1, and the subscript in the following means and so on; calculating to obtain new parameter filtering valueAnd its variance Pt|tThe subscript represents the estimator of epoch t time calculated from the observation values (including t time) in front of epoch t; and calculating the after-the-fact residual errorThe subscript represents the estimation quantity calculated by the observation value (including the time t) in front of the epoch t, wherein the estimation quantity is obtained by the time t of the epoch; finally, performing one-way smoothing (Lag-one smoothing) to obtain a parameter smoothing estimation valueAnd its variance Pt-1|tAnd anAndcovariance matrix P oft,t-1|t. To be provided withPt|t、Pt-1|t、Pt，t-1|tAnd post-mortem residualAs an input value of Step3, used to calculate random walk noise in the next Step, and the third Step is entered.
Step3. estimate new ZWD random walk parameters. First, according toPt|t、Pt-1|t、Pt,t-1|tAnd post-mortem residualThe mathematical expectation given the current a priori ZWD random walk noise and observation noise is calculated.
Then, the mathematical expectation is taken as a pseudo observed value, the pseudo observed value is inserted into the tail of a sliding window queue, and then maximization is carried out, namely, the data in the sliding window is applied to estimate new ZWD random walk parameters. And finally deleting the head data of the sliding window, substituting the new estimation parameters into Step1 to update the prior ZWD random walk parameters, and repeating the steps from Step1 to Step3.
As mentioned above, Step3 is the main innovation point of the algorithm of the present invention, and compared with the prior art, the method of the present invention has the following advantages:
1) in the existing algorithm, the random walk noise of the ZWD is set as a fixed value, and the space-time difference of the ZWD in different areas at different times is not considered, so that the estimation of the real-time ZWD estimation value is inaccurate, and the phase of the ZWD time sequence lags. The method of the invention estimates the ZWD random walk noise and observation noise variance component in real time mainly according to the Bayes principle, and improves the ZWD by about 11 percent.
2) The algorithm of the invention has high time resolution, and particularly, when extreme weather occurs, the method of the invention can make local peak values of ZWD more obvious, and can more reliably monitor and forecast severe weather events.
3) When the existing square difference component estimation technology is used for directly estimating the random walk noise, convergence is difficult, and generally only after-treatment can be carried out. The algorithm of the invention not only can estimate the random walk noise model in real time, but also is more stable and more reliable in result compared with the similar maximum likelihood estimation method.
In order to better understand the technical scheme of the application, the implementation steps of the technology are described in detail below by combining a specific example. This example forms an ionospheric-free combination LC and PC to simulate real-time estimation ZWD with GPS constellations L1 and L1 carriers, P1 and P1 pseudo-random noise codes as basic observations.
Referring to fig. 1 to fig. 5, an improved GNSS tropospheric wet delay estimation method according to an embodiment of the present invention includes the following steps:
step 1: the method comprises the steps of downloading a rinex observation file in advance, and preparing a GPS precision satellite orbit and clock error file, an earth rotation parameter ERP matched with the clock error and the orbit, a GPS code difference file DCB, an IGS antenna file, a JPL planet ephemeris file and a tide FES2004 file.
Step 2: to elaborate on this step, it is necessary to first give a definition of ZWD time-varying random walk noise and then introduce the state equations and observation equations corresponding to precise point positioning.
1) Definition of time-varying random walk noise
Considering tropospheric random noise sequence ZWDtI t is 0, 1, …, and satisfies
ZWDt＝ZWDt-1+wt (6)
Wherein, wtIndependent and zero normal distribution from the mean, whose covariance depends linearly on the scaling factor q and ZWDtAnd ZWDt-1Sampling interval Δ t, i.e.
E{wt}＝0，Cov{wt，wt}＝qt·Δt (7)
Then { ZWDtI t 0, 1, … is defined as random walk noise, parameter qtTime-varying is the Spectral Density (SPD) that describes the characteristics of the time-varying random walk noise.
2) State equation definition
Suppose that in the precise point positioning PPP, Kalman filtering is applied to process data, and the state equation is as follows:
wherein, the right arrow represents a vector, tandt-1 represents two epoch time before and after, t is 1, …, n;is an r X1 dimensional state vector, X, Y and Z are geocentric coordinate components of the receiver, respectively; the cdt represents the receiver clock difference and,representing the ambiguity parameter of the ionospheric-free combination, mtIndicating the number of visible satellites. PhitRepresenting a state transition matrix;is a process noise vector, follows a gaussian distribution,the mean value of (a) is zero, and the variance matrix satisfies the condition:
wherein Q isk,t(k ═ 1, …, p) is a diagonal matrix or scalar; alpha is alphak(k-1, …, p) represents process noise QtThe variance component of (a) is a parameter to be solved. It can be seen that the random walk parameter q of the convective wet zenith delay ZWDtThe requirement of the formula (4) is met.
3) Observation equation definition
Reading an epoch observation data, and performing data preprocessing to form the following ionosphere-free combination:
wherein the notation is the same as in formulas (1) and (2). Equations (8) and (9) are non-linear and can be estimated using extended kalman filtering. For convenience of description, equations (8) and (9) are written in a general matrix form:
whereinIs an unknown vector containing the receiver coordinates, receiver clock, ambiguities and troposphere ZWD.Is composed of an ionosphere-free layerThe constructed observation vector is used as a vector for observing,is the observation noise vector.RtIs the observed noise matrix at time t:
wherein R isk,t(k ═ 1, …, q) is a diagonal matrix or scalar; beta is ak(k-1, …, q) represents process noise RtThe variance component of (a) is a parameter to be solved.
According to the formula (6-11), the predicted value of the state vector can be obtained according to the extended Kalman filtering theoryAnd its predicted variance Pt|t-1(ii) a Further observation and update are carried out to obtain the estimation value prediction value of the state vectorSum covariance variance Pt|t(ii) a Finally, one-step smoothing is carried out to obtain a smooth state vectorSum covariance variance Pt-1|tAnd further turning to Step3 to estimate random noise or variance components.
Step 3：
Step 3.1 calculating the mathematical expectation
Suppose thatSatisfy the requirement oft＝1，…，n。
WhereinIs the observation data that is to be observed,is an unknown variable, becauseNot directly observable, is called hidden variable.
Will be provided withIt is treated as complete data and its likelihood function is assumedDependent on the unknown parameter vector Θ ═ αk(k＝1,…,p),βj(j ═ 1, …, q) }, assuming Θ(t-1)Representing a priori unknown parameters estimated using historical data, a mathematical expectation can be obtained:
wherein the content of the first and second substances,is the after-the-fact residual error,and Pt|tFiltering the parameter and its variance;and Pt-1|tSmoothing the estimate and its variance P for the parametert-1|t，Pt,t-1|tIs thatAndthe covariance matrix of (a) is obtained,
step 3.2 reaction of XtAnd YtAnd respectively inserting the tail parts of the sliding window queues.
Step3.3 maximize the M-Step.
From the queue elements in the sliding window, the variance component is calculated according to the following formula:
step 3.3.1 sliding window variance component estimation
Formula (16) and formula (17) use historical data in a certain time to recurrently update the current epoch parameter, the parameter estimation value of the current epoch uses the prior value of the next iteration, and the sliding window can not only track the rapidly changing time-varying parameter, but also can rapidly filter the outlier.
The variables L and N and the sliding window sizes used to emphasize process noise and observation noise are not the same because they may have different rates of change. For example, since the variance component is almost constant for the observed values, N should take a large value, e.g., greater than 2000 epochs, and the random walk noise L should be relatively small. According to practical experience, L should exceed 10 minutes, and adjusting the size of the sliding window will change the temporal behavior of the unknowns.
Step3.3.2 boundary constraint
In order to strengthen the numerical stability and ensure the convergence of the algorithm, it is necessary to set upper and lower limits on the parameters to be estimated. The precision of GPS carrier phase measurement is about 2 to 3 cm, which is much higher than random walk noiseRandom walk noise is easily absorbed into GPS carrier measurement, so that the estimation value of the random walk noise is not reduced, and the variance component of observation noise is increased continuously. Without the boundary constraint, the random walk noise can be underestimated, which in turn can result in the tropospheric zenith total delay being underestimated. Generally, the upper and lower bounds of a parameter can be obtained by analyzing the physical characteristics of such parameter. However, it is difficult to specifically quantify model errors due to geographical differences and the complexity of weather changes. However, observing the GNSS phase residuals gives a hint to indirectly constrain the unknown parameters. When strong convection weather occurs, the 15 minute RMS average of the GPS phase residuals may fluctuate between 11 mm and 25 mm. Thus, if the phase residual is less than a certain valueThe carrier phase of a value (e.g. 15 mm) to update the parameter for the purpose of (indirectly) constraining the parameter change.
The embodiment of the invention provides a GNSS troposphere time-varying random walk noise estimation method, and relates to the field of extreme weather prediction by applying GNSS; the GNSS troposphere real-time random walk parameter estimation technology under extreme weather can be independent of a numerical weather model to estimate the GNSS troposphere random walk parameters in real time on the basis of an expectation maximization algorithm and Kalman filtering. The method is independent of a numerical weather model, and random walk parameters of troposphere Wet-sky top Delay (ZWD) can be estimated in real time only by applying a GNSS precision single-point positioning technology PPP; extracting random walk noise parameters by adopting sliding window lengths with different lengths for noise and observation noise components in the Kalman filtering process in a recursion mode; and boundary constraint conditions of the GNSS carrier variance component are specified to ensure the convergence and reliability of the whole algorithm.
Step 3.4 deletes the sliding window queue head element in preparation for the next evaluation as shown in fig. 2.
The present invention uses precision point-positioning (PPP) to calculate zenith wet delay. The PPP mapping function is a global mapping function (Boehm et al, 2006a) and the PPP is implemented based on a Kalman filter, equation (16,17) is used to estimate the parameters of random walk noise. To ensure convergence, the size of the sliding window used for observation in (17) should be greater than 2000 epochs, and the length of the sliding window for random walk noise in equation (16) is about 10 minutes. Estimation of random walk noise parameter q using only carrier phase observations with residuals less than 15 mmt。
The test scheme provided by the embodiment of the invention specifically comprises the following steps: 13 IGS stations are selected, randomly and uniformly distributed over different locations on earth with different weather (see fig. 3). Wherein the atmospheric river passes twice over the survey station SCIP, resulting in the SCIP undergoing two heavy rainfall events. The PPP observation equations are standard GPS dual-ionosphere-free combinations, expressed in PL and LC as code pseudorange and carrier phase observations, respectively (units: meters), that have corrected satellite clocks, relativistic effects, solid earth tide, polar tide, ocean tide, phase wrap, a priori tropospheric delay and phase closure. In (13), the stochastic nature of PC and LC observations is considered to be dependent on altitude angle, but with different variance components, namely:
wherein the content of the first and second substances,is the observation noise of the satellite s at epoch t,is the elevation angle, beta, from the receiver to the satellite1andβ2Respectively representing variance components PC and LC variance components, mtIs the number of satellites in view.
For all GPS stations except SCIP, the final GFZ Multi-GNSS (GBM) product (https:// cddis. nasa. gov/archive/gnss/products/mgex /) that uses satellite orbits and clocks. For the survey station SCIP, the orbits and satellite clocks of the ephemeris are downloaded from the center CODE (http:// ftp. aiub. unebe. ch/CODE/2017 /). A sample interval of 15 minutes for accurate orbit correction and 30 minutes for accurate clock, are generated at the GFZ. The present invention uses CODE Differential CODE Bias (DCB) product to align C1 to P1, using frequency bands of L1 and L2 for GPS.
Figure 4 shows the estimates of the random walk parameters (excluding the station SCIP) over a week. The horizontal axis is the number of epochs, the spectral density qtSquare root of squareThe results show that the estimated values of IGS stations at different geographical locations are significantly different, roughly to the extentThe smallest and most stable value is the HRAG station, and its time series is much smoother than other time series, because HRAG usesThe hydrogen is the primary clock and the HRAG is located in desert areas, where the upper troposphere is relatively stable. Other sites show greater interference compared to the HRAG, with the magnitude of korr being approximately three times that of the HRAG. In addition, at ZIM2 stations, the ZWD improvement may reach 11%.
Fig. 5 shows the station SCIP zenith troposphere delay ZTD time series with two peaks indicating extreme weather-atmospheric river two passes through the station SCIP. O denotes the ZTD time series estimated using the normal PPP algorithm, o denotes the results obtained using the algorithm of the invention, post ZTD products published by the black + curve JPL. SCIP experiences two extreme rainfall days due to atmospheric river landing on the west coast of north america twice in 2017 at 1 month. It can be seen that the ZTD and JPL time series estimated by the algorithm of the present invention exhibit similar behavior with an average shift of 3.9mm in the vertical axis direction and an RMS of about 7.8mm, both solutions clearly capturing the abrupt changes in the troposphere. From 2017, day 1, month 18 to 1 month 19, the ZTD first rose sharply by about 10cm, then the ZTD value decreased immediately thereafter as atmospheric river flow exited. Likewise, ZTD again increases dramatically as the second atmospheric river arrives between 2017 at 1/22 and 2017 at 1/23. It can be clearly seen that the ZTD time series of the invention are more pronounced than the peaks of the JPL ZTD strong precipitation periods, which can differ by up to about 2 cm in peak value, and therefore, the algorithm of the invention is more advantageous in monitoring and forecasting severe weather events. Furthermore, it can be seen that the ZTD value of normal PPP creates a time lag when extreme weather occurs, which is also less than the ZTD estimate of the present invention. Of the three solutions, the solution of the embodiment of the present invention is significantly better.
In summary, the core innovation point of the invention is that 1) according to the expectation-maximization algorithm and the Kalman filtering theory, the invention develops an algorithm framework capable of estimating the GNSS troposphere random walk noise in real time. 2) The algorithm is independent of a numerical weather model, and the calculation is efficient. 3) The algorithm adopts a sliding window and an indirect constraint carrier phase variance component to estimate GNSS troposphere random walk noise, thereby ensuring the convergence and numerical reliability of the algorithm. 4) Compared with the ZTD value estimated by the traditional GNSS precise single-point positioning technology, the algorithm can improve the ZTD by about 11 percent. 5) Particularly, the algorithm can remarkably improve the ZTD estimated value in extreme weather, and can more effectively utilize GNSS to monitor and forecast the extreme weather.
The invention discloses a new algorithm framework capable of estimating Random Walk Noise (RWN) in real time when GNSS is applied to extreme weather prediction. In current GNSS data processing software, RWN parameters are typically given empirically as constants. However, different RWN parameters and corresponding observation weights can cause huge differences in ZWD estimated values under extreme weather, and real-time estimation of the ZWD time-varying optimal RWN parameters is the key for accurately estimating the ZWD values under extreme weather and monitoring and forecasting extreme weather events. The invention takes Kalman filtering and expectation maximization algorithm as the basis, considers the characteristic that GNSS observation noise and multipath depend on altitude angle, fully considers the geographical position difference of GNSS observation stations and the time-varying characteristic of rapid change of troposphere, can estimate the optimal random model parameter changing along with time in real time, improves the ZWD by about 11 percent, makes the local peak value of the ZWD more obvious, and can more reliably monitor and forecast severe weather events.
As will be appreciated by one skilled in the art, embodiments of the present application may be provided as a method, system, or computer program product. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.
The present application is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the application. It will be understood that each flow and/or block of the flow diagrams and/or block diagrams, and combinations of flows and/or blocks in the flow diagrams and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art can make modifications and equivalents to the embodiments of the present invention without departing from the spirit and scope of the present invention, which is set forth in the claims of the present application.