An Atmospheric Delay Correction Method Based on GNSS Higher-Order Horizontal Gradient and Particle Filtering

Abstract

This invention provides an atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering, comprising: GNSS data acquisition; obtaining GNSS ZTD using parameter estimation; determining a mapping function and correcting the mapping function using the ETILTING method to obtain anisotropic delay; correcting the GNSS ZTD using the anisotropic delay to obtain a corrected GNSS ZTD; performing particle filtering on the corrected GNSS ZTD and gradient observation data to construct linear relationship models between the observed ZTD value and the observed GNSS delay gradient parameters and altitude; projecting the regularly gridded sea-level ZTD onto InSAR coordinates after convolution interpolation to obtain a neutral atmospheric delay correction model; subtracting the neutral atmospheric delay correction model from the estimated delay models at different times to obtain the delay difference; and correcting the interferogram using the delay difference to obtain a corrected interferogram. The method described in this invention demonstrates significant advantages in reducing the impact of atmospheric delay and improving the accuracy of monitoring minute surface displacements.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing technology and atmospheric science, specifically to a neutral atmospheric delay correction model for synthetic aperture radar interferometry (InSAR) based on the zenith total delay (ZTD) and its horizontal gradient of the Global Navigation Satellite System (GNSS). Background Technology

[0002] Interferometric synthetic aperture radar (InSAR) technology, as an advanced space geodesy tool, has been widely used to monitor small deformations on the Earth's surface. Despite its powerful capabilities, InSAR measurement results may still be affected by interference from the Earth's atmosphere, a phenomenon commonly referred to as atmospheric propagation delay. This effect limits the accuracy of InSAR technology in displacement detection, typically to the centimeter level or higher [1].

[0003] Atmospheric delay effects are mainly divided into two parts: one part is caused by ionospheric disturbances. It is worth noting that microwave signals actually arrive at the receiver earlier in the ionosphere, rather than being delayed; the other part is caused by the delay caused by the neutral atmosphere. In this field, recent technological advances have successfully developed some robust correction methods to handle ionospheric phase changes in InSAR. In particular, the range-segmented spectrum method (SSM)[2], which utilizes the frequency modulation characteristics of SAR, has shown its effectiveness in correcting InSAR ionospheric delay in high-frequency and low-frequency SAR systems[3-5]. However, neutral atmospheric delay remains a major noise source in InSAR measurements because it is closely related to the total amount of water vapor along the microwave path. Water vapor in the atmosphere, especially in the lower troposphere, exhibits significant fluctuations in both space and time. In some extreme cases, as pointed out by Kinoshita et al.[6] and Kinoshita and Furuya[7], when cumulonimbus clouds develop, atmospheric delay can affect InSAR line-of-sight (LOS) changes by more than 20 cm.

[0004] To improve the accuracy of InSAR in measuring surface displacement down to the millimeter level, establishing an accurate atmospheric delay model is essential. This scheme aims to utilize GNSS atmospheric observation data, including the total tropospheric delay (ZTD) and horizontal gradient. A new model is developed to mitigate the impact of InSAR atmospheric delay by using a mapping function determined by the third-order continued fraction MFlsmcom and adding higher-order horizontal gradients via the ETILTING method. This model is expected to improve the accuracy and reliability of InSAR technology in monitoring minute surface deformations.

[0005] References:

[0006] [1]Zebker H A,Rosen P A,Hensley S.Atmospheric effects ininterferometric synthetic aperture radar surface deformation and topographicmaps[J].Journal of geophysical research:solid earth,1997,102(B4):7547-7563.

[0007] [2]Gomba G,Parizzi A,De Zan F,et al.Toward operational compensationof ionospheric effects in SAR interferograms:The split-spectrum method[J].IEEE Transactions on Geoscience and Remote Sensing,2015,54(3):1446-1461.

[0008] [3]Zhang B,Wang C,Ding X,et al.Correction of ionospheric artifacts inSAR data:Application to fault slip inversion of 2009southern Sumatraearthquake[J].IEEE Geoscience and Remote Sensing Letters,2018,15(9):1327-1331.

[0009] [4]Furuya M,Suzuki T,Maeda J,et al.Midlatitude sporadic-E episodesviewed by L-band split-spectrum InSAR[J].Earth,Planets and Space,2017,69:1-10.

[0010] [5]Gomba G, González FR, De Zan F.Ionospheric phase screencompensation for the Sentinel-1TOPS and ALOS-2ScanSAR modes[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 55(1):223-235.

[0011] [6]Kinoshita Y, Furuya M, Hobiger T, et al. Are numerical weather model outputs helpful to reduce tropospheric delay signals in InSAR data? [J].Journal of Geodesy,2013,87:267-277.

[0012] [7]Kinoshita Y,Furuya M.Localized delay signals detected by syntheticaperture radar interferometry and their simulation by WRF 4DVAR[J].SOLA,2017,13:79-84. Summary of the Invention

[0013] To address the aforementioned technical problems, this invention provides an atmospheric delay correction method based on GNSS high-order horizontal gradients and particle filtering. This method can establish an accurate atmospheric delay correction model, thereby improving the accuracy and reliability of InSAR technology in monitoring minute surface deformations.

[0014] The technical solution provided by this invention is as follows.

[0015] In a first aspect, the present invention provides an atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering, comprising the following steps:

[0016] GNSS data acquisition;

[0017] The GNSS ZTD is obtained by estimating GNSS data using parameter estimation methods;

[0018] The mapping function is determined by a method based on the joint determination of the coefficients of the mapping function using the third-order continued fraction; the mapping function is then modified using the ETILTING method to obtain the anisotropic delay.

[0019] The GNSS ZTD is corrected by using anisotropic delay to obtain the corrected GNSS ZTD;

[0020] Particle filtering is applied to the corrected GNSS ZTD and gradient observation data to obtain the state estimate of the particle filter and the coefficients of the linear relationship model; linear relationship models of the observed ZTD value and the observed GNSS delay gradient parameter with altitude are constructed respectively.

[0021] A regular gridded sea level ZTD model is constructed based on a linear relationship model; the regular gridded sea level ZTD model is then convolved and interpolated and projected onto InSAR coordinates to obtain a neutral atmospheric delay correction model.

[0022] The delay difference applicable to the interferogram is obtained by subtracting the neutral atmospheric delay correction model from the estimated delay models at different times.

[0023] The corrected interferogram is obtained by correcting the interferogram using the delay difference.

[0024] In one possible implementation, the calculation formula for the method of estimating GNSS data using parameter estimation is as follows:

[0025] ZTD = ZHD·M dry (E)+ZWD·M wet (E)

[0026] Where ZTD is the total tropospheric delay; ZHD is the zenith tropospheric dry component delay; M dry (E) is the projection function of the dry component of the zenith troposphere, E is the satellite elevation angle; ZWD is the delay of the wet component of the zenith troposphere; M wet (E) is the projection function of the wet component of the zenith troposphere.

[0027] In one possible implementation, the method for determining the mapping function using a joint determination method for the coefficients of the mapping function based on a third-order continued fraction includes the following steps:

[0028] The highly accurate MFlsmabc algorithm is used to simultaneously estimate the coefficients a, b, and c. If the MFlsmabc algorithm converges, the modeling process ends; if it does not converge, the MFlsmab algorithm is used instead for model construction.

[0029] Evaluate the convergence status of the MFlsmab and MFlsma algorithms; if neither converges, use the MFfast algorithm for computation.

[0030] In one possible implementation, the method of modifying the mapping function using the ETILTING method includes the following steps:

[0031] Calculate the anisotropic delay at different azimuth and elevation angles;

[0032] Construct a design matrix that contains functions of elevation and azimuth to correlate the observations with the horizontal gradient parameters, and use the least squares method to estimate the horizontal gradient parameters.

[0033] Determine the model expression for the ETILTING method to modify the mapping function, solve for the horizontal gradient parameters, including the first and second order horizontal gradient parameters in the north-south and east-west directions, and obtain the anisotropic delay.

[0034] The expression for anisotropic delay is:

[0035]

[0036] Where ΔG represents the anisotropic delay, MF is the mapping function, e is the satellite elevation angle, i is the i-th order, M is the total order, φ represents the azimuth angle, and G N and G E This represents the GNSS delay gradient parameters in the north-south and east-west directions.

[0037] In one possible implementation, the linear relationship between the observed ZTD value and the observed GNSS delay gradient parameter and altitude is modeled as follows:

[0038] The observed linear relationship between ZTD values ​​and altitude is shown in the following model:

[0039] ZTD obs =ZTD0(i,j)+a·h(i,j)

[0040] In the formula, ZTD obs The observed ZTD values ​​are: ZTD0 is the sea level at the regular grid points, a is the height dependence coefficient, h is the GNSS station height, and i and j represent the grid indices in the east and north directions, respectively.

[0041] Observed GNSS delay gradient parameter G N and G E The linear relationship model between altitude and height is as follows:

[0042]

[0043]

[0044] Among them, H ZTD D represents the scale height of the ZTD. e D represents the grid spacing in the east-west direction. n This indicates the grid spacing in the north-south direction.

[0045] In one possible implementation, the method of constructing a regular gridded sea level ZTD model based on a linear relationship model; and projecting the regular gridded sea level ZTD onto InSAR coordinates after convolutional interpolation to obtain a neutral atmospheric delay correction model includes:

[0046] A regular gridded sea level ZTD model is constructed based on a linear relationship model to obtain an adjusted ZTD model;

[0047] The sea level ZTD is obtained by using cubic convolution interpolation on the adjusted ZTD model, denoted as ZTD0;

[0048] The estimated ZTD for each InSAR pixel is calculated using ZTD0, denoted as ZTD. model ;

[0049] Using a simple trigonometric function of the InSAR incident angle, the InSAR slant path delay is calculated by projecting from the zenith direction to the LOS direction, thus obtaining a neutral atmospheric delay correction model.

[0050]

[0051] Among them, STD InSAR θ represents the InSAR slant path delay, and θ represents the radar incident angle.

[0052] Secondly, the present invention provides an atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering, comprising:

[0053] Acquisition module, used for GNSS data acquisition;

[0054] The estimation module is used to estimate GNSS data using parameter estimation methods to obtain the GNSS ZTD;

[0055] The anisotropic delay acquisition module is used to determine the mapping function using a joint determination method of mapping function coefficients based on the third-order continued fraction; and to modify the mapping function using the ETILTING method to obtain the anisotropic delay.

[0056] The correction module is used to correct the GNSS ZTD using anisotropic delay, so as to obtain the corrected GNSS ZTD.

[0057] The linear relationship determination module is used to perform particle filtering on the corrected GNSS ZTD and gradient observation data to obtain the state estimate of the particle filter and obtain the coefficients of the linear relationship model; and to construct linear relationship models between the observed ZTD value and the observed GNSS delay gradient parameter and altitude, respectively.

[0058] The model acquisition module is used to construct a regular gridded sea level ZTD model based on a linear relationship model; after convolution interpolation, the regular gridded sea level ZTD model is projected onto InSAR coordinates to obtain a neutral atmospheric delay correction model.

[0059] The delay difference acquisition module is used to subtract the neutral atmospheric delay correction model from the estimated delay models at different times to obtain the delay difference suitable for the interferogram.

[0060] The correction module is used to correct the interferogram using the delay difference to obtain the corrected interferogram.

[0061] Thirdly, the present invention provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in the first aspect.

[0062] Fourthly, the present invention provides a non-transitory computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in the first aspect.

[0063] Fifthly, the present invention provides a computer program product, including a computer program that, when executed by a processor, implements the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in the first aspect.

[0064] The beneficial effects of this invention are as follows:

[0065] (1) The method described in this invention has shown significant advantages in reducing the impact of atmospheric delay and improving the accuracy of monitoring minute surface displacements.

[0066] (2) The method described in this invention determines the mapping function based on the third-order continuous fraction MFlsmcom, and the ETILTING method determines the anisotropic delay to improve the correction accuracy.

[0067] (3) The method described in this invention uses particle filtering to extract the linear height dependence from the sea level ZTD distribution and projects it onto InSAR coordinates for neutral atmospheric delay correction; by using cubic convolution interpolation on the grid to interpolate GNSS ZTD data, the atmospheric delay correction accuracy of each pixel in the InSAR image is improved. Attached Figure Description

[0068] Figure 1 This is an overall flowchart of the method described in this invention;

[0069] Figure 2Flowchart for modeling the MFlsmcom method;

[0070] Figure 3 A schematic diagram of the atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering provided by the present invention;

[0071] Figure 4 This is a schematic diagram of the structure of the electronic device provided by the present invention. Detailed Implementation

[0072] The present invention will be further described below with reference to specific embodiments, but the content of the present invention is not limited thereto.

[0073] like Figure 1 As shown, an atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering includes the following steps:

[0074] S100, GNSS (Global Navigation Satellite System) data acquisition.

[0075] In one possible implementation, the GNSS data originates from the Nevada Geological Laboratory at the University of Nevada, Reno.

[0076] S200 uses parameter estimation to estimate GNSS data and obtains GNSS ZTD (Zenith Tropospheric Delay).

[0077] In one possible implementation, the calculation formula for the method of estimating GNSS data using parameter estimation is as follows:

[0078] ZTD = ZHD·M dry (E)+ZWD·M wet (E)

[0079] Where ZTD is the total zenith delay in the troposphere; ZHD is the zenith tropospheric dry component delay; M dry (E) is the projection function of the dry component of the zenith troposphere, E is the satellite elevation angle; ZWD is the delay of the wet component of the zenith troposphere; M wet (E) is the projection function of the wet component of the zenith troposphere.

[0080] S300 uses the joint determination method of mapping function coefficients based on third-order continued fractions (MFlsmcom) to determine the mapping function; the ETILTING method is used to modify the mapping function to obtain the anisotropic delay.

[0081] In one possible implementation, S300 utilizes a method for determining the mapping function based on a joint determination method of the mapping function coefficients of a third-order continued fraction, such as... Figure 2 As shown, it includes the following steps:

[0082] (1) The highly accurate MFlsmabc algorithm is used to simultaneously estimate the coefficients a, b and c;

[0083] (2) Convergence evaluation: If the MFlsmabc algorithm converges, the mapping function construction ends; if it does not converge, the MFlsmab algorithm is used to construct the mapping function.

[0084] (3) Evaluate the convergence status of the MFlsmab algorithm and the MFlsma algorithm in turn. If neither converges, the more efficient MFfast algorithm is used to construct the mapping function.

[0085] In one possible implementation, the method of modifying the mapping function using the ETILTING method (ExtendedTILTING, a multi-parameter modeling method for horizontal gradients that takes into account higher-order variation terms) in S300 includes the following steps:

[0086] (1) Parameter estimation: Construct a design matrix (the design matrix contains partial derivatives related to the horizontal gradient parameters to be estimated) and an observation vector (the observation vector is composed of actual observation data, containing the values ​​of anisotropic delay observed at different elevation angles and azimuth angles), solve the linear least squares problem to obtain the estimated values ​​of the horizontal gradient parameters, including the first-order and second-order horizontal gradient parameters.

[0087] (2) The TILTING method is extended to the ETILTING method by adding a second-order horizontal gradient parameter to better fit the anisotropic delay higher-order variation term caused by local topographic disturbances and weather events. The expression for the ETILTING method is: Where AD(e,φ) represents the anisotropic delay, MF is the mapping function, e is the satellite elevation angle, i is the i-th order, M is the total order, φ represents the azimuth angle, and G... N and G E This represents the i-th order horizontal gradient parameter in the north-south and east-west directions.

[0088] It should be noted that the mapping function is determined according to the aforementioned method for jointly determining the coefficients of the mapping function based on the third-order continued fraction (MFlsmcom), where the mathematical model of the MFlsmcom method is expressed as:

[0089] v n×1 =Ax-l,P

[0090] Where A is the design matrix, and its expression is:

[0091]

[0092] In the formula, n represents the number of mapping factors involved in the modeling, and the matrix represents the partial derivative of the third-order continued fraction with respect to the coefficients a, b, and c, which includes the coefficients a, b, and c themselves. Initial values ​​for the coefficients a, b, and c need to be given.

[0093]

[0094] h0 represents the initial value in height, and w0 represents the initial value in width.

[0095] x is the matrix of parameters to be estimated:

[0096]

[0097] l is a constant matrix:

[0098]

[0099] In the formula, The mapping factor, representing the input modeled by the mapping function, is obtained through ray tracing. This represents the prior mapping factor, which is calculated from the initial values ​​of the coefficients a, b, and c.

[0100] P is the weight matrix of the mapping factor:

[0101] P n×n =E n×n

[0102] Where E represents the prior mapping factor matrix, which is calculated from the initial values ​​of the coefficients a, b, and c.

[0103] The construction process of the MFlsma model is similar to that of the MFlsmabc model. The main difference is that the MFlsma model sets the values ​​of b and c to known prior values ​​and uses the least squares method to estimate the coefficient of a. If a mapping function with an elevation angle of 3° or an initial elevation angle of 3.3° is incorporated into the application of the MFlsma model, then the model simplifies to the MFfast model. The MFlsmab method only fixes the coefficient of c and uses the least squares method to estimate the coefficients of a and b.

[0104] Based on the above, a suitable mapping function is selected according to convergence.

[0105] S400 uses anisotropic delay to correct the GNSS ZTD, obtaining the corrected GNSS ZTD. The correction formula is:

[0106] ZTD(e, φ) = ZHD·MF h(e)+ZWD·MF w (e)+AD(e, φ)

[0107] In the formula, ZTD represents the total tropospheric delay, ZHD represents the dry tropospheric delay, ZWD represents the wet tropospheric delay, and MF represents the total tropospheric delay. h MF is the tropospheric dry delay gradient mapping function. w φ is the tropospheric wet delay gradient mapping function, e is the satellite elevation angle, and φ represents the azimuth angle.

[0108] S500 performs particle filtering on the corrected GNSS ZTD and gradient observation data to obtain the state estimate of the particle filter and the coefficients of the linear relationship model; and constructs linear relationship models between the observed ZTD value and the observed GNSS delay gradient parameter and altitude.

[0109] The observed linear relationship between ZTD values ​​and altitude is shown in the following model:

[0110] ZTD obs =ZTD0(i,j)+a·h(i,j)

[0111] In the formula, ZTD obs The observed ZTD values ​​are: ZTD0 is the sea level at the regular grid points, a is the height dependence coefficient, h is the GNSS station height, and i and j represent the grid indices in the east and north directions, respectively.

[0112] Observed GNSS delay gradient parameter G N and G E The linear relationship model between altitude and height is as follows:

[0113]

[0114]

[0115] Among them, H ZTD D represents the scale height of the ZTD. e D represents the grid spacing in the east-west direction. n G represents the grid spacing in the north-south direction. N and G E These represent the gradient parameters in the north-south and east-west directions, respectively.

[0116] In one possible implementation, the state estimate M' of the particle filter is:

[0117]

[0118] in, For sampling, denoted as sample weight, and N as the number of particles.

[0119] The state estimate M' includes the coefficients ZTD0,a (high dependence coefficient) from the linear relationship model, H ZTD D e and D n These coefficients are obtained through particle filtering to construct linear relationships between the observed ZTD values ​​and the observed GNSS delay gradient parameters and altitude, which are then used for atmospheric delay correction.

[0120] In one possible implementation, particle filtering includes the following steps:

[0121] Weight Calculation: For each particle, its weight is calculated. The weight is obtained by comparing the predicted observations of the particle with the actual GNSS observation data.

[0122] Resampling: Resampling is performed based on the particle weights to reduce particle diversity and concentrate on more likely state regions. This step allows particles with higher weights to have a greater chance of being copied, while particles with lower weights may be discarded.

[0123] S600 constructs a regular gridded sea level ZTD model based on a linear relationship model; after convolution interpolation, the regular gridded sea level ZTD is projected onto InSAR coordinates to obtain a neutral atmospheric delay correction model.

[0124] In one possible implementation, S600 includes the following steps:

[0125] S610, Construction of a Regular Gridded Sea Level ZTD Model: Construct a regular grid covering the entire GNSS station distribution area; each grid point represents a sea level ZTD value ZTD0, consistent with ZTD0 in step S500;

[0126] S620, Adjust ZTD value: Adjust the ZTD value of each grid point using the linear relationship model between the observed ZTD value obtained from S500 and the height;

[0127] S630, Adjust gradient parameters: Apply the observed GNSS delay gradient parameter G obtained from S500. N and G E The model is refined by using a linear relationship model with height to obtain the adjusted ZTD model. These gradient parameters take into account the variation of ZTD values ​​in different directions, thereby improving the spatial resolution and accuracy of the model.

[0128] S640, Convolutional Interpolation: Perform convolutional interpolation on the adjusted ZTD model to obtain the sea level ZTD (ZTD0), in order to smooth the discontinuity of ZTD values ​​and improve the continuity of the model. In this embodiment, cubic convolutional interpolation is used.

[0129] S650 uses ZTD0 to calculate the estimated ZTD for each InSAR pixel, denoted as ZTD. model ;

[0130] S660, using a simple trigonometric function of the InSAR incident angle, projects from the zenith direction to the LOS direction to calculate the InSAR (STD) value. InSAR Obtain a neutral atmospheric delay correction model by slanting path delay;

[0131]

[0132] Among them, STD InSAR Indicates InSAR slant path delay, ZTD model Let θ represent the estimated ZTD for each InSAR pixel, and θ represent the radar incident angle.

[0133] S700 subtracts the estimated delays at different times from the atmospheric delay correction model to obtain the delay difference applicable to the interferogram.

[0134] It should be noted that the estimated delays at different times are obtained by collecting GNSS data sources from different times. GNSS (Global Navigation Satellite System) data sources refer to satellite signal data used to determine the location and time of a point on Earth. These data sources are collected at different times, meaning they are collected at different times and reflect the satellite positions, signal propagation conditions, and the state of the Earth's atmosphere (especially the troposphere) at different points in time. Therefore, since the GNSS data sources are from different times, the resulting atmospheric delays are also from different times.

[0135] S800, by correcting the interferogram through delay difference, obtains the corrected interferogram.

[0136] It should be noted that the interferogram is obtained by synthetic aperture radar interferometry (InSAR) technology through image registration, removal of flattening effects, atmospheric correction, and phase unwrapping. The corrected interferogram can be used in fields such as topographic mapping, surface deformation monitoring, precipitation prediction, and landslide monitoring.

[0137] The atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering provided by the present invention will be described below. The atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering described below can be referred to in correspondence with the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering described above.

[0138] Figure 3 This is a schematic diagram of the atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering provided in an embodiment of the present invention, as shown below. Figure 3As shown, it includes: an acquisition module 31, an estimation module 32, an anisotropic delay acquisition module 33, a correction module 34, a linear relationship determination module 35, a model acquisition module 36, a delay difference acquisition module 37, and a correction module 38, wherein:

[0139] Acquisition module 31 is used for GNSS data acquisition;

[0140] Estimation module 32 is used to estimate GNSS data using parameter estimation methods to obtain GNSS ZTD;

[0141] The anisotropic delay acquisition module 33 is used to determine the mapping function using a joint determination method of mapping function coefficients based on the third-order continuous fraction; and to modify the mapping function using the ETILTING method to obtain the anisotropic delay.

[0142] Correction module 34 is used to correct the GNSS ZTD using anisotropic delay to obtain the corrected GNSS ZTD;

[0143] The linear relationship determination module 35 is used to perform particle filtering on the corrected GNSS ZTD and gradient observation data to obtain the state estimate of the particle filter and obtain the coefficients of the linear relationship model; and to construct linear relationship models between the observed ZTD value and the observed GNSS delay gradient parameter and the altitude, respectively.

[0144] Model acquisition module 36 is used to construct a regular gridded sea level ZTD model based on a linear relationship model; after convolution interpolation, the regular gridded sea level ZTD is projected onto InSAR coordinates to obtain a neutral atmospheric delay correction model.

[0145] The delay difference acquisition module 37 is used to subtract the neutral atmospheric delay correction model from the estimated delay models at different times to obtain the delay difference suitable for the interferogram.

[0146] The correction module 38 is used to correct the interferogram through the delay difference to obtain the corrected interferogram.

[0147] Figure 4 An example is a schematic diagram of the physical structure of an electronic device, such as... Figure 4 As shown, the electronic device may include a processor 410, a communications interface 420, a memory 430, and a communication bus 440. The processor 410, communications interface 420, and memory 430 communicate with each other via the communication bus 440. The processor 410 can call logical instructions from the memory 430 to execute an atmospheric delay correction method based on GNSS high-order horizontal gradients and particle filtering.

[0148] Furthermore, the logical instructions in the aforementioned memory 430 can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, essentially, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0149] On the other hand, the present invention also provides a computer program product, which includes a computer program that can be stored on a non-transitory computer-readable storage medium. When the computer program is executed by a processor, the computer is able to execute the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering provided by the above methods.

[0150] In another aspect, the present invention also provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, is implemented to perform the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering provided by the above methods.

[0151] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0152] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0153] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. An atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering, characterized in that, Includes the following steps: GNSS data acquisition; The GNSS ZTD is obtained by estimating GNSS data using parameter estimation methods; A mapping function is determined using a joint determination method based on the coefficients of a third-order continued fraction; the mapping function is then corrected using the ETILTING method to obtain anisotropic delays. The method for determining the mapping function using the joint determination method based on the coefficients of a third-order continued fraction includes the following steps: coefficients a, b, and c are simultaneously estimated using the highly accurate MFlsmabc algorithm; if the MFlsmabc algorithm converges, the modeling process ends; if it does not converge, the MFlsmab algorithm is used for model construction; the convergence status of the MFlsmab and MFlsma algorithms is evaluated; if both converge... If convergence fails, the MFfast algorithm is used for calculation. The method of using the ETILTING method to correct the mapping function and obtain the anisotropic delay includes the following steps: calculating the anisotropic delay for different azimuth and elevation angles; constructing a design matrix containing functions of elevation and azimuth angles to correlate the observations with the horizontal gradient parameters, and using the least squares method to estimate the horizontal gradient parameters; determining the model expression of the ETILTING method to correct the mapping function, solving for the horizontal gradient parameters, including the first and second order horizontal gradient parameters in the north-south and east-west directions, to obtain the anisotropic delay. The expression for anisotropic delay is: in, The function represents the anisotropic delay, MF is the mapping function, e is the satellite elevation angle, i is the i-th order, and M is the total order. Indicates azimuth. and GNSS delay gradient parameters representing the north-south and east-west directions; The GNSS ZTD is corrected by using anisotropic delay to obtain the corrected GNSS ZTD; Particle filtering is applied to the corrected GNSS ZTD and gradient observation data to obtain the state estimate of the particle filter and the coefficients of the linear relationship model; linear relationship models of the observed ZTD value and the observed GNSS delay gradient parameter with altitude are constructed respectively. A regular gridded sea level ZTD model is constructed based on a linear relationship model; the regular gridded sea level ZTD model is then convolved and interpolated and projected onto InSAR coordinates to obtain a neutral atmospheric delay correction model. The delay difference applicable to the interferogram is obtained by subtracting the neutral atmospheric delay correction model from the estimated delay models at different times. The corrected interferogram is obtained by correcting the interferogram using the delay difference.

2. The atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering according to claim 1, characterized in that, The calculation formula for the method of estimating GNSS data using parameter estimation is as follows: Where ZTD is the total tropospheric delay; ZHD is the zenith tropospheric dry component delay; ZWD is the projection function of the dry component of the zenith troposphere, where E is the satellite elevation angle and ZWD is the delay of the wet component of the zenith troposphere. This is the projection function of the wet component of the zenith troposphere.

3. The atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering according to claim 1, characterized in that, The observed ZTD values ​​and observed GNSS delay gradient parameters have the following linear relationships with altitude: The observed linear relationship between ZTD values ​​and altitude: In the formula, ZTD obs These are the observed ZTD values, where ZTD0 represents the sea level at the regular grid points. a , where h is the GNSS station height, and i and j represent the grid indices in the east and north directions, respectively; Observed GNSS delay gradient parameter G N and G E Linear relationship with altitude: Among them, H ZTD Indicates the scale height of ZTD. D e Indicates the grid spacing in the east-west direction. D n This indicates the grid spacing in the north-south direction.

4. The atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering according to claim 1, characterized in that, The method for constructing a regular gridded sea level ZTD model based on a linear relationship model, and then projecting the regular gridded sea level ZTD onto InSAR coordinates after convolution interpolation to obtain a neutral atmospheric delay correction model includes: A regular gridded sea level ZTD model is constructed based on a linear relationship model to obtain an adjusted ZTD model; The sea level ZTD is obtained by using cubic convolution interpolation on the adjusted ZTD model, denoted as ZTD0; The estimated ZTD for each InSAR pixel is calculated using ZTD0, denoted as ZTD. model ; Using a simple trigonometric function of the InSAR incident angle, the InSAR slant path delay is calculated by projecting from the zenith direction to the LOS direction, thus obtaining a neutral atmospheric delay correction model. Among them, STD InSAR Indicates the InSAR slant path delay. θ Indicates the radar incident angle.

5. An atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering, characterized in that, The apparatus is used to implement the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in any one of claims 1 to 4, comprising: Acquisition module, used for GNSS data acquisition; The estimation module is used to estimate GNSS data using parameter estimation methods to obtain the GNSS ZTD; The anisotropic delay acquisition module is used to determine the mapping function using a joint determination method of mapping function coefficients based on the third-order continued fraction; and to modify the mapping function using the ETILTING method to obtain the anisotropic delay. The correction module is used to correct the GNSS ZTD using anisotropic delay, so as to obtain the corrected GNSS ZTD. The linear relationship determination module is used to perform particle filtering on the corrected GNSS ZTD and gradient observation data to obtain the state estimate of the particle filter and obtain the coefficients of the linear relationship model; and to construct linear relationship models between the observed ZTD value and the observed GNSS delay gradient parameter and altitude, respectively. The model acquisition module is used to construct a regular gridded sea level ZTD model based on a linear relationship model; after convolution interpolation, the regular gridded sea level ZTD model is projected onto InSAR coordinates to obtain a neutral atmospheric delay correction model. The delay difference acquisition module is used to subtract the neutral atmospheric delay correction model from the estimated delay models at different times to obtain the delay difference suitable for the interferogram. The correction module is used to correct the interferogram using the delay difference to obtain the corrected interferogram.

6. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in any one of claims 1 to 4.

7. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in any one of claims 1 to 4.

8. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering as described in any one of claims 1 to 4.

Images