Method for retrieving tropospheric wet delay and atmospheric water vapor content over polar sea ice with techdemosat-1 satellite grazing angle spaceborne global navigation satellite system reflectometry

Abstract

A method for retrieving tropospheric wet delay and atmospheric water vapor content over polar sea ice with TDS-1 satellite grazing angle spaceborne GNSS-R is provided, including: Si, obtaining TDS-1 GNSS-R raw intermediate frequency signal data, VMF3 grid data, GPT3 grid data and ERA5 data; S2, correcting an error of tropospheric wet delay estimation of grazing angle spaceborne GNSS-R; S3, constructing a grazing angle spaceborne GNSS-R tropospheric wet delay estimation model; S4, calculating grazing angle spaceborne GNSS-R ZWD; S5, calculating a T.sub.m value of a target point based on GPT3 model, substituting the T.sub.m value into a conversion factor II, and combining calculated GNSS-R ZWD to obtain a GNSS-R IWV estimated value; and S6, verifying inversion performance of GNSS-R ZWD and integrated water vapor (IWV) by using reference data.

Claims

1. A method for retrieving tropospheric wet delay and atmospheric water vapor content over polar sea ice with TDS-1 satellite grazing angle spaceborne GNSS-R, comprising the following steps: Step 1, obtaining TDS-1 GNSS-R raw intermediate frequency signal data, VMF3 grid data, GPT3 grid data and ERA5 data; Step 2, correcting an error of a tropospheric wet delay estimation of the grazing angle spaceborne GNSS-R, comprising: preprocessing raw intermediate frequency signals of the TDS-1 satellite by performing closed-loop tracking of direct signals and computing reflected signal carrier phases using an open-loop tracking model; removing dynamic phase components comprises receiver clock deviation, orbit error, ionospheric delay, and carrier phase ambiguity by approximating orbit error through a time-dependent correction function and performing ionospheric delay correction using total electron content information obtained from a global ionospheric map and ionospheric reference model data; reconstructing reflected signal phases by segmenting the residual carrier phase into different frequency components using wavelet transform, mitigating signal discontinuities caused by cycle slips through interpolation, and suppressing phase fluctuations using a filtering process; and correcting residual phase discontinuities by applying a phase unwrapping procedure based on an extended Kalman filter, implemented with a nonlinear signal model configured for in-phase and quadrature data, to produce a continuous and stabilized carrier phase measurement; Step 3, constructing a grazing angle spaceborne GNSS-R tropospheric wet delay estimation model by combining a grazing angle spaceborne GNSS-R phase measurement theoretical model and a slant tropospheric delay theoretical model considering a mapping function; specifically comprises following sub-steps: Step 3.1, constructing the grazing angle spaceborne GNSS-R phase measurement theoretical model, wherein an expression is:
.sub.R(t)=|{right arrow over (r.sub.tx)}(tt.sub.R){right arrow over (r.sub.sp)}(t)|+|{right arrow over (r.sub.sp)}(t){right arrow over (r.sub.tx)}(t)|+I.sub.R(t)+T.sub.R(t)b.sub.tx(tt.sub.R)b.sub.rx(t)(9) wherein {right arrow over (r.sub.tx)}, {right arrow over (r.sub.sp)} and {right arrow over (r.sub.rx)} respectively represent a position vector of a GNSS satellite, a mirror point and a position vector of a TDS-1 satellite in Earth-Centered, Earth-Fixed (ECEF) coordinates, and b.sub.tx and b.sub.rx are clock deviations of a GNSS transmitter and a TDS-1 satellite receiver respectively, I.sub.R is influence of ionosphere on a carrier phase, and T.sub.R is tropospheric delay; after correcting errors, noises and cycle slips caused by ionospheric propagation, a GNSS satellite orbit and clock error, a TDS-1 orbit and clock error, residual carrier phases are mainly affected by the tropospheric delay in atmosphere, and a tropospheric information measurement equation is expressed as:
{circumflex over (T)}.sub.R(t)={circumflex over ()}.sub.R(t)(|{right arrow over (r.sub.tx)}(tt.sub.R){right arrow over (r.sub.sp)}(t)|+|{right arrow over (r.sub.sp)}(t){right arrow over (r.sub.tx)}(t)|)+b.sub.tx(tt.sub.R)b.sub.rx(t)I.sub.R(t)(10) wherein {circumflex over ()}.sub.R is a reflected signal phase distance; Step 3.2, estimating GNSS-R slant tropospheric delay by using residual phases of reflected signals:
{circumflex over (T)}.sub.R(t)={circumflex over ()}.sub.R(t).sub.R(t).sub.R(t)+M(t)+(t)(11) wherein represents an estimation of variables in equation (9), .sub.R represents geometric and clock components in the equation (9), M represents an unknown integer of a phase ambiguity, and represents influence of various estimation errors and noises; the slant tropospheric delay is modeled as: $\begin{array}{cc}{T}_R\left(t\right)=2\left(\mathrm{ZHD}\left(t\right)m_{\mathrm{dry}}\left(\left(t\right)\right)+\mathrm{ZWD}\left(t\right)m_{\mathrm{wet}}\left(\left(t\right)\right)+L_g\left(e,\right)\right)& \left(12\right)\end{array}$ $\begin{array}{cc}\left(t\right)=\frac{\frac{1}{2}{\hat{T}}_R\left(\left(t\right)\right)-\left(t\right)m_{\mathrm{dry}}\left(\left(t\right)\right)-L_g\left(e,\right)-M\left(t\right)}{m_{\mathrm{wet}}\left(\left(t\right)\right)}& \left(13\right)\end{array}$ wherein m.sub.dry and m.sub.wet are mapping functions of dry delay and wet delay in troposphere, L.sub.a(e,) is an error caused by a gradient change of the atmosphere in a horizontal direction, e is a cut-off altitude angle of a satellite, a is an azimuth angle of the satellite, and are obtained by interpolation based on grid data, and then obtaining a following expression: $\begin{array}{cc}{\hat{T}}_R\left(t\right)=\frac{1}{2}{\hat{T}}_R\left(\left(t\right)\right)-\left(t\right)m_{\mathrm{dry}}\left(\left(t\right)\right)-\left(t\right)m_{\mathrm{wet}}\left(\left(t\right)\right)-L_g\left(e,\right);& \left(14\right)\end{array}$ Step 4, calculating grazing angle spaceborne GNSS-R ZWD; Step 5, calculating a T.sub.m value of a target point based on a GPT3 model, substituting the T.sub.m value into a conversion factor II, and combining calculated GNSS-R ZWD to obtain a GNSS-R IWV estimated value; and Step 6, verifying inversion performance of GNSS-R ZWD and IWV by using reference data.

2. The method for retrieving the tropospheric wet delay and the atmospheric water vapor content over the polar sea ice with the TDS-1 satellite grazing angle spaceborne GNSS-R according to claim 1, wherein a raw intermediate frequency signal data set of TDS-1 comprises raw intermediate frequency signals reflected by single-frequency multi-constellation GNSS (GPS/Galileo/BeiDou-3); the VMF3 grid data is 1 tropospheric wet delay data estimated by a VMF model, GPT3 grid data is atmospheric weighted average temperature data estimated by the GPT3 model, and ERA5 data is atmospheric water vapor content data.

3. The method for retrieving the tropospheric wet delay and the atmospheric water vapor content over the polar sea ice with the TDS-1 satellite grazing angle spaceborne GNSS-R according to claim 1, wherein the carrier phase {circumflex over ()}.sub.R,i.sup.OL of the reflected signals calculated by the open-loop tracking model is expressed by a following formula:
{circumflex over ()}.sub.R,i.sup.OL={tilde over ()}.sub.R+{circumflex over ()}.sub.R,i+.sub.1N.sub.R,i.sup.OL(1) wherein N.sub.R,i.sup.OL is phase integer ambiguity, and {tilde over ()}.sub.R is a reflected signal geometric model comprising cycle slips and signal dynamic deviation; modeling and eliminating most of the dynamic phase components before the cycle slip correction; wherein a receiver clock deviation is estimated from a direct signal measurement, and an ionospheric correction formula is derived by using vertical total electron content (vTEC) model provided by global ionospheric map (GIM) of international GNSS service (IGS) as follows: $\begin{array}{cc}{I}_R=-\frac{40.3}{f_{L1}^2}\left({\mathrm{sTEC}}_1+{\mathrm{sTEC}}_2-{\mathrm{sTEC}}_D\right)& \left(2\right)\end{array}$ wherein for GPS L1 signals, f.sub.L1=0.57542 GHz, sTEC.sub.D represents a slant total electron content (sTEC) along a direct signal path, and sTEC.sub.1 and sTEC.sub.2 represent sTEC along incident light and reflected light respectively; wherein for reflection over the polar sea ice, a sea ice surface elevation is modeled by using a digital elevation model (DEM) or IceSat-2 data products, and a precise GNSS satellite orbit and a clock deviation are obtained from a final orbit of Centre for Orbit Determination in Europe (CODE); the orbit error is approximated by a linear trend, and a GNSS-R satellite orbit error is parameterized as a linear function of time:
T.sub.Orb.sup.Fit(t)=a.sub.0+a.sub.1t(3) wherein a.sub.0 and a.sub.1 are estimated by least square fitting as follows: $\begin{array}{cc}\left(a_0,a_1\right)=\underset{a_0^{*},a_1^{*}}{\arg \min }{.Math.}_{k=1}^M{\left({\hat{}}_R\left(t_k\right)-a_0^{*}-a_1^{*}{t}_k\right)}^2& \left(4\right)\end{array}$ wherein t.sub.k is a time stamp and M is a number of carrier phase samples of GNSS reflected signals; then the reflected signal geometric model removing the signal dynamic deviation is expressed as:
{circumflex over ()}R=(|{right arrow over (r.sub.tx)}(tt.sub.R){right arrow over (r.sub.sp)}(t)|+|{right arrow over (r.sub.sp)}(t)r.sub.tx(t)|)+b.sub.tx(tt.sub.R)b.sub.rx(t)(5) wherein {right arrow over (r.sub.tx)}, {right arrow over (r.sub.sp)}, and {right arrow over (r.sub.rx)} respectively represent the position vector of the GNSS satellite, the mirror point and the position vector of the TDS-1 satellite in the Earth-Centered, Earth-Fixed (ECEF) coordinates, and b.sub.tx and b.sub.rx are the clock deviations of the GNSS transmitter and the TDS-1 satellite receiver respectively; obtaining the residual phases of the reflected signals with reduced dynamic components by using a following formula:
{circumflex over ()}.sub.R,i.sup.OL={circumflex over ()}.sub.R,i.sup.OL{tilde over ()}.sub.R+{circumflex over ()}.sub.R(6) wherein the {circumflex over ()}.sub.R,i.sup.OL is processed by the wavelet transform, and the signals are decomposed into components with different scales and frequencies by the wavelet transform to obtain time-frequency characteristics of the signals, and cycle slip detection is realized by detecting sudden changes or abnormal phenomena in the signals; when the cycle slips occur, data are interpolated or extrapolated to fill missing values or expand a data range, and then phase smoothing filter is applied to deal with noise or fluctuation in the data; in a phase unwrapping process insensitive to the noise, an extended Kalman filter (EKF) algorithm is used to process residual phases of spaceborne GNSS-R measurement, and an observation model is expressed as an in-phase/quadrature (I/Q) nonlinear function; the algorithm is capable of estimating and reducing integer period deviation, reducing measurement noise and realizing cycle slip mitigation; finally, a carrier phase measurement result after correcting or reducing cycle slip deviation is as follows:
{circumflex over ()}.sub.R,i.sup.Corr={circumflex over ()}R(t)+{circumflex over ()}.sub.R,i.sup.Corr(7) wherein {circumflex over ()}.sub.R,i.sup.Corr is a residual phase of a reflected signal after the cycle slip correction.

4. The method for retrieving the tropospheric wet delay and the atmospheric water vapor content over the polar sea ice with the TDS-1 satellite grazing angle spaceborne GNSS-R according to claim 3, wherein the sTEC is calculated by the GIM by a following formula: $\begin{array}{cc}{\mathrm{sTEC}}_{1,2,D}={}_{1,2,D}\frac{{\mathrm{vTEC}}_{1,2,D}}{\cos \left({}_{1,2,D}\right)}& \left(8\right)\end{array}$ wherein is a proportional factor of ionospheric density above the TDS-1 orbit, is an elevation angle of incident ray at an ionospheric penetration point (IPP), vTEC is vertical TEC at the IPP obtained from the GIM, and subscripts 1, 2 and D represent incident light, reflected light and direct path respectively; values of .sub.1, .sub.2 and .sub.D are determined based on an international reference ionosphere (IRI) 2020 model.

5. The method for retrieving the tropospheric wet delay and the atmospheric water vapor content over the polar sea ice with the TDS-1 satellite grazing angle spaceborne GNSS-R according to claim 1, wherein a calculation expression of the grazing angle spaceborne GNSS-R ZWD in the step 4 is: $\begin{array}{cc}\left(t\right)=\frac{{\hat{T}}_R\left(t\right)-\stackrel{\_}{{\hat{T}}_R}}{m_{\mathrm{wet}}\left(\left(t\right)\right)}+\left(t\right)& \left(15\right)\end{array}$ wherein - represents an average value of time series.

6. The method for retrieving the tropospheric wet delay and the atmospheric water vapor content over the polar sea ice with the TDS-1 satellite grazing angle spaceborne GNSS-R according to claim 1, wherein a calculation formula of the step 5 is as follows: $\begin{array}{cc}\mathrm{IWV}=\mathrm{II}\mathrm{ZWD}& \left(16\right)\end{array}$ $\begin{array}{cc}\mathrm{II}=\frac{{10}^6}{R_V.Math.\left(k_2^{}+k_3/T_m\right)}& \left(17\right)\end{array}$ $\begin{array}{cc}{k}_2^{}=k_2-k_1\frac{R_d}{R_v},R_d=R/M_d& \left(18\right)\end{array}$ wherein R.sub.v=461.522 J kg.sup.1, k.sub.1=77.604 K.Math.hPa.sup.1, k.sub.2=64.79 K.Math.hPa.sup.1, k.sub.3=373900 K.sup.2.Math.hPa.sup.1, R=8.31434 J/(mol.Math.K), is an ideal universal constant, M.sub.a=28.964410.sup.3 kg/mol, is a molar mass of dry air, and T.sub.m is an atmospheric weighted average temperature calculated by the GPT3 model.

7. The method for retrieving the tropospheric wet delay and the atmospheric water vapor content over the polar sea ice with the TDS-1 satellite grazing angle spaceborne GNSS-R according to claim 1, wherein the step 6 is: verifying retrieval performance of the GNSS-R ZWD and the IWV by using reference ZWD data and reference IWV data; wherein the reference ZWD data is the VMF3 grid data and the reference IWV data is the ERA5 data.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The accompanying drawings, which are a part of this application, are used to provide a further understanding of the disclosure. The illustrative embodiments of the disclosure and their descriptions are used to explain the disclosure, but they do not constitute an improper limitation of the disclosure. Obviously, the attached drawings in the following description are only some embodiments, and other drawings may be obtained according to these drawings without creative work for ordinary people in the field. In the attached drawings:

(2) FIG. 1 is a flow chart of the present disclosure.

(3) FIG. 2A is a comparison result of tropospheric wet delay estimated in the VMF3 grid model and retrieved from data of trajectory 1 in GNSS-R of TDS-1 satellite in the embodiment of the present disclosure.

(4) FIG. 2B is a comparison result of tropospheric wet delay estimated in the VMF3 grid model and retrieved from data of trajectory 2 in GNSS-R of TDS-1 satellite in the embodiment of the present disclosure.

(5) FIG. 3A is a comparison result between the atmospheric water vapor content retrieved from the spaceborne GNSS-R data of trajectory 1 of TDS-1 satellite and the atmospheric water vapor content provided by ERA5 in the embodiment of the present disclosure.

(6) FIG. 3B is a comparison result between the atmospheric water vapor content retrieved from the spaceborne GNSS-R data of trajectory 2 of TDS-1 satellite and the atmospheric water vapor content provided by ERA5 in the embodiment of the present disclosure.

(7) It should be noted that these drawings and written descriptions are not intended to limit the conceptual scope of the disclosure in any way, but to explain the concept of the disclosure to those skilled in the art by referring to specific embodiments.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(8) In order to make the objective, technical scheme and advantages of the embodiment of the present disclosure more clear, the technical scheme in the embodiment of the present disclosure will be described clearly and completely with the accompanying drawings. The following embodiments are used to illustrate the present disclosure, but are not used to limit the scope of the present disclosure.

Embodiment 1

(9) In order to verify the effectiveness of the proposed method, the TDS-1 raw intermediate frequency signal observation data, VMF3 grid data, GPT3 grid data and ERA5 data from Sep. 1, 2014 to Mar. 25, 2019 are obtained for experiments, and the experimental results of the disclosure are compared with reference data.

(10) The disclosure relates to a method for retrieving tropospheric wet delay and atmospheric water vapor content over polar sea ice with TDS-1 satellite grazing angle spaceborne GNSS-R, including the following steps.

(11) S1, obtaining TDS-1 GNSS-R raw intermediate frequency signal data, VMF3 grid data, GPT3 grid data and ERA5 data; where the acquisition time of TDS-1 raw intermediate frequency signal data set is from Sep. 1, 2014 to Mar. 25, 2019, including the raw intermediate frequency signals reflected by single-frequency multi-constellation GNSS (GPS/Galileo/BeiDou-3). the VMF3 grid data is 1 tropospheric wet delay data estimated by a VMF model, GPT3 grid data is atmospheric weighted average temperature data estimated by the GPT3 model, and ERA5 data is atmospheric water vapor content data.

(12) S2, correcting an error of a tropospheric wet delay estimation of the grazing angle spaceborne GNSS-R; correcting the error of the tropospheric wet delay estimation of the grazing angle spaceborne GNSS-R in the S2 includes dynamic phase components and deviation removal and phase reconstruction; where the dynamic phase components and deviation removal includes clock deviation, orbit error, ionospheric error and phase ambiguity, and the phase reconstruction includes cycle slip identification based on wavelet transform and cycle slip correction based on extended Kalman filter, including the following sub-steps.

(13) S2.1, before correcting each error, processing raw intermediate frequency observation data of TDS-1 satellite, including closed-loop processing of direct signals and calculation of open-loop (OL) tracking model of reflected signals; where carrier phase {circumflex over ()}.sub.R,i.sup.OL of the reflected signals calculated by the open-loop tracking model is expressed by a following formula:
{circumflex over ()}.sub.R,i.sup.OL={tilde over ()}.sub.R+{circumflex over ()}.sub.R,i+.sub.1N.sub.R,i.sup.OL(1) where N.sub.R,i.sup.OL is phase integer ambiguity, {circumflex over ()}.sub.R,i is residual carrier phase estimation, .sub.1GNSS signal wavelength, and {tilde over ()}.sub.R is a reflected signal geometric model including cycle slips and signal dynamic deviation.

(14) S2.2, modeling and eliminating dynamic phase components before the cycle slip correction; where a receiver clock deviation is estimated from a direct signal measurement, the current TDS-1 satellite only generates a single frequency measurement result, so an ionospheric correction formula is derived by using vertical total electron content (vTEC) model provided by global ionospheric map (GIM) of international GNSS service (IGS) as follows:

(15) $\begin{array}{cc}{I}_R=-\frac{40.3}{f_{L1}^2}\left({\mathrm{sTEC}}_1+{\mathrm{sTEC}}_2-{\mathrm{sTEC}}_D\right)& \left(2\right)\end{array}$ where for GPS L1 signals, f.sub.L=1.57542 GHz, sTEC.sub.D represents a slant total electron content (sTEC) along the direct signal path, and sTEC.sub.1 and sTEC.sub.2 represent sTEC along incident light and the reflected light respectively. sTEC is calculated by the GIM and calculated by a following formula:

(16) $\begin{array}{cc}{\mathrm{sTEC}}_{1,2,D}={}_{1,2,D}\frac{{\mathrm{vTEC}}_{1,2,D}}{\cos \left({}_{1,2,D}\right)}& \left(3\right)\end{array}$ where is a proportional factor of ionospheric density above TDS-1 orbit, is an elevation angle of incident ray at an ionospheric penetration point (IPP), vTEC is vertical TEC at the IPP obtained from the GIM, and subscripts 1, 2 and D represent incident light, reflected light and direct path respectively; values of .sub.1, .sub.2 and .sub.D are determined based on international reference ionosphere (IRI) 2020 model.

(17) For the reflection over the polar sea ice, a sea ice surface elevation is modeled by using digital elevation model (DEM) or IceSat-2 data products, and a precise GNSS satellite orbit and a clock deviation are obtained from a final orbit of Centre for Orbit Determination in Europe (CODE). Due to the lack of precise orbit information of TDS-1, the orbit error is approximated by a linear trend, that is, the orbit error is approximated by a linear trend, and a GNSS-R satellite orbit error is parameterized as a linear function of time:
T.sub.Orb.sup.Fit(t)=a.sub.0+a.sub.1t(4) where .sub.0 and .sub.1 are estimated by least square fitting as follows:

(18) 0$\begin{array}{cc}\left(a_0,a_1\right)=\arg \min {.Math.}_{k=1}^M{\left({\hat{}}_R\left(t_k\right)-a_0^{*}-a_1^{*}{t}_k\right)}^2& \left(5\right)\end{array}$ where t.sub.k is a time stamp and M is a number of carrier phase samples of GNSS reflected signals. In addition to the linear components in the orbit error term, the linear components in other error terms such as carrier integer ambiguity and delay correction residual are removed.

(19) Then the reflected signal geometric model removing the signal dynamic deviation is expressed as:
{circumflex over ()}R=(|{right arrow over (r.sub.tx)}(tt.sub.R){right arrow over (r.sub.sp)}(t)|+|{right arrow over (r.sub.sp)}(t)r.sub.tx(t)|)+b.sub.tx(tt.sub.R)b.sub.rx(t)(6) where {right arrow over (r.sub.tx)}, {right arrow over (r.sub.sp)} and {right arrow over (r.sub.rx)} respectively represent position vector of GNSS satellite, mirror point and position vector of TDS-1 satellite in Earth-Centered, Earth-Fixed (ECEF) coordinates, and b.sub.tx and b.sub.rx are clock deviations of GNSS transmitter and TDS-1 satellite receiver respectively.

(20) S2.3, obtaining residual phases of the reflected signals with reduced dynamic components by using the following formula:
{circumflex over ()}.sub.R,i.sup.OL={circumflex over ()}.sub.R,i.sup.OL{tilde over ()}.sub.R+{circumflex over ()}.sub.R(7) herein the {circumflex over ()}.sub.R,i.sup.OL is processed by the wavelet transform, and the signals are decomposed into components with different scales and frequencies by the wavelet transform to obtain time-frequency characteristics of the signals, and cycle slip detection is realized by detecting sudden changes or abnormal phenomena in the signals; when cycle slips occur, data are interpolated or extrapolated to fill missing values or expand a data range, and then phase smoothing filter is applied to deal with noise or fluctuation in the data; in a phase unwrapping process insensitive to the noise, an extended Kalman filter (EKF) algorithm is used to process residual phases of spaceborne GNSS-R measurement, and the observation model is expressed as an in-phase/quadrature (I/Q) nonlinear function; the algorithm is capable of estimating and reducing integer period deviation, reducing measurement noise and realizing cycle slip mitigation; finally, carrier phase measurement result after correcting or reducing cycle slip deviation is as follows:
{circumflex over ()}.sub.R,i.sup.Corr={circumflex over ()}R(t)+{circumflex over ()}.sub.R,i.sup.Corr(8) where {circumflex over ()}.sub.R,i.sup.Corr is the residual phase of the reflected signal after the cycle slip correction.

(21) S3, constructing a grazing angle spaceborne GNSS-R tropospheric wet delay estimation model by combining a grazing angle spaceborne GNSS-R phase measurement theoretical model and a slant tropospheric delay theoretical model considering a mapping function, specifically including following sub-steps:

(22) S3.1, constructing a grazing angle spaceborne GNSS-R phase measurement theoretical model, and an expression is:
.sub.R(t)=|{right arrow over (r.sub.tx)}(tt.sub.R){right arrow over (r.sub.sp)}(t)|+|{right arrow over (r.sub.sp)}(t){right arrow over (r.sub.tx)}(t)|+I.sub.R(t)+T.sub.R(t)b.sub.tx(tt.sub.R)b.sub.rx(t)(9) where {right arrow over (r.sub.tx)}, {right arrow over (r.sub.sp)} and {right arrow over (r.sub.rx)} respectively represent the position vector of the GNSS satellite, the mirror point and the position vector of the TDS-1 satellite in the Earth-Centered, Earth-Fixed (ECEF) coordinates, and b.sub.tx and b.sub.rx are the clock deviations of the GNSS transmitter and the TDS-1 satellite receiver respectively, I.sub.R is influence of ionosphere on carrier phase, and T.sub.R is tropospheric delay.

(23) In order to obtain accurate carrier phase estimation for tropospheric delay estimation, it is necessary to correct the errors, noises and cycle slips caused by ionospheric propagation, GNSS satellite orbit and clock error, TDS-1 orbit and clock error. After correcting these errors, the residual carrier phase is mainly affected by tropospheric delay in the atmosphere. Furthermore, the tropospheric information measurement equation may be expressed as:
{circumflex over (T)}.sub.R(t)={circumflex over ()}.sub.R(t)(|{right arrow over (r.sub.tx)}(tt.sub.R){right arrow over (r.sub.sp)}(t)|+|{right arrow over (r.sub.sp)}(t){right arrow over (r.sub.tx)}(t)|)+b.sub.tx(tt.sub.R)b.sub.rx(t)I.sub.R(t)(10) where {circumflex over ()}.sub.R is a reflected signal phase distance.

(24) S3.2, estimating GNSS-R slant tropospheric delay by using the residual phases of the reflected signals:
{circumflex over (T)}.sub.R(t)={circumflex over ()}.sub.R(t).sub.R(t).sub.R(t)+M(t)+(t)(11) where represents an estimation of variables in equation (9), .sub.R represents geometric and clock components in the equation (9), M represents an unknown integer of the phase ambiguity, and represents influence of various estimation errors and noises. The slant tropospheric delay is modeled as:

(25) $\begin{array}{cc}{T}_R\left(t\right)=2\left(\mathrm{ZHD}\left(t\right)m_{\mathrm{dry}}\left(\left(t\right)\right)+\mathrm{ZWD}\left(t\right)m_{\mathrm{wet}}\left(\left(t\right)\right)+L_g\left(e,\right)\right)& \left(12\right)\end{array}$$\begin{array}{cc}\left(t\right)=\frac{\frac{1}{2}{\hat{T}}_R\left(\left(t\right)\right)-\left(t\right)m_{\mathrm{dry}}\left(\left(t\right)\right)-L_g\left(e,\right)-M\left(t\right)}{m_{\mathrm{wet}}\left(\left(t\right)\right)}& \left(13\right)\end{array}$ where m.sub.dry and m.sub.wet are mapping functions of dry delay and wet delay in the troposphere, L.sub.g(e,) is an error caused by a gradient change of the atmosphere in a horizontal direction, e is a cut-off altitude angle of a satellite, is an azimuth angle of the satellite, and are obtained by interpolation based on grid data, and then obtaining a following expression:

(26) $\begin{array}{cc}{\hat{T}}_R\left(t\right)=\frac{1}{2}{\hat{T}}_R\left(\left(t\right)\right)-\left(t\right)m_{\mathrm{dry}}\left(\left(t\right)\right)-\left(t\right)m_{\mathrm{wet}}\left(\left(t\right)\right)-L_g\left(e,\right).& \left(14\right)\end{array}$

(27) S4, calculating grazing angle spaceborne GNSS-R ZWD, and the calculation expression is:

(28) $\begin{array}{cc}\left(t\right)=\frac{{\hat{T}}_R\left(t\right)-\stackrel{\_}{{\hat{T}}_R}}{m_{\mathrm{wet}}\left(\left(t\right)\right)}+\left(t\right)& \left(15\right)\end{array}$ where - represents an average value of time series.

(29) S5, the T.sub.m value of the target point is calculated based on the GPT3 model, and the T.sub.m value is substituted into the conversion factor II, and then the estimated value of GNSS-R IWV is obtained by combining the calculated GNSS-R ZWD. The calculation formula is as follows:

(30) $\begin{array}{cc}\mathrm{IWV}=\mathrm{II}\mathrm{ZWD}& \left(16\right)\end{array}$$\begin{array}{cc}\mathrm{II}=\frac{{10}^6}{R_V.Math.\left(k_2^{}+k_3/T_m\right)}& \left(17\right)\end{array}$$\begin{array}{cc}{k}_2^{}=k_2-k_1\frac{R_d}{R_v},R_d=R/M_d& \left(18\right)\end{array}$ where R.sub.v=461.522 J kg.sup.1 K.sup.1, k.sub.1=77.604 K.Math.hPa.sup.1, k.sub.2=64.79 K.Math.hPa.sup.1, k.sub.3=373900 K.sup.2.Math.hPa.sup.1, R=8.31434 J/(mol.Math.K), R is an ideal universal constant, M.sub.d=28.964410.sup.3 kg/mol, M.sub.d is a molar mass of dry air, and T.sub.m is an atmospheric weighted average temperature.

(31) S6, using the reference ZWD data and the reference IWV data to verify the inversion performance of GNSS-R ZWD and IWV; where the reference ZWD data is the VMF3 grid data and the reference IWV data is the ERA5 data.

(32) FIG. 2A-FIG. 2B show the estimated ZWD of ZWD and VMF3 grid model retrieved from the two trajectory data of TDS-1 satellite spaceborne GNSS-R. FIG. 3A-FIG. 3B show the IWV estimated by GNSS-R ZWD and the IWV of the interpolated ERA5 data along two trajectories.

(33) As can be seen from FIG. 2A-FIG. 2B and FIG. 3A-FIG. 3B, the inversion results of GNSS-R ZWD and the estimation results of VMF3 grid model have good consistency and high accuracy, and gradients are observed from two GNSS-R orbits with similar steepness. In addition, the inversion results of GNSS-R IWV are in good agreement with those of ERA5 IWV. These preliminary results show that it is feasible to invert ZWD and IWV in high latitudes by using TDS-1 grazing angle spaceborne GNSS-R raw intermediate frequency signal data.

(34) The above is only the preferred embodiment of the present disclosure, and it does not limit the present disclosure in any form. Although the present disclosure has been disclosed in the preferred embodiment, it is not used to limit the present disclosure. Any person familiar with this patent may make some changes or modify it into an equivalent embodiment by using the technical content suggested above without departing from the technical scheme of the present disclosure. However, any simple modification to the above embodiment is made according to the technical essence of the present disclosure.