Atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering

Abstract

The invention provides an atmospheric delay correction method based on GNSS (Global Navigation Satellite System) high-order horizontal gradient and particle filtering. The atmospheric delay correction method comprises the following steps: acquiring GNSS data; gNSS ZTD is obtained by using a parameter estimation method; determining a mapping function, and correcting the mapping function by using an ETILTING method to obtain anisotropic delay; the GNSS ZTD is corrected by using the anisotropic delay, and the corrected GNSS ZTD is obtained; performing particle filtering on the corrected GNSS ZTD and gradient observation data, and constructing a linear relation model between the observed ZTD value and the height and a linear relation model between the observed GNSS delay gradient parameter and the height; after convolution interpolation is carried out on the regularly-gridded sea level ZTD, the sea level ZTD is projected to an InSAR coordinate, and a neutral atmospheric delay correction model is obtained; subtracting the neutral atmospheric delay correction model from the estimation delay models at different times to obtain a delay difference; and correcting the interference pattern through the delay difference to obtain a corrected interference pattern. The method disclosed by the invention has obvious advantages in the aspects of reducing the atmospheric delay influence and improving the ground surface infinitesimal displacement monitoring precision.

Description

Technical Field

[0001] The present invention relates to the fields of remote sensing technology and atmospheric science, and in particular to a neutral atmosphere delay correction model of synthetic aperture radar interferometry (InSAR) based on global navigation satellite system (GNSS) zenith total delay (ZTD) and its horizontal gradient. Background Art

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

[0003] The atmospheric delay effect can be divided into two parts: one is caused by ionospheric disturbances. It is worth noting that microwave signals actually arrive at the receiver earlier than delayed in the ionosphere; the other is the delay caused by the neutral atmosphere. In this area, recent technological advances have successfully developed some robust correction methods to deal with ionospheric phase changes in InSAR. In particular, the range-sliced ​​spectrum method (SSM) [2], which exploits 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 is still 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, the atmospheric delay may affect the InSAR line-of-sight (LOS) variation by more than 20 cm.

[0004] In order to improve the accuracy of InSAR in surface displacement measurement to the millimeter level, it is essential to establish an accurate atmospheric delay model. This scheme aims to use GNSS atmospheric observation data, including the total tropospheric delay (ZTD) and horizontal gradient, to determine the mapping function based on the third-order continued fraction MFlsmcom, and to add high-order horizontal gradients using the ETILTING method to develop a new model to mitigate the impact of InSAR atmospheric delay. Through this model, we expect to improve the accuracy and reliability of InSAR technology in monitoring small 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] In view of the above technical problems, the present invention provides an atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering. The method can establish an accurate atmospheric delay correction model, thereby improving the accuracy and reliability of InSAR technology in monitoring small surface deformations.

[0014] The technical solution provided by the present 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 collection;

[0017] Use parameter estimation method to estimate GNSS data and obtain GNSS ZTD;

[0018] The mapping function is determined by using a method for jointly determining the mapping function coefficients based on a third-order continued fraction; the mapping function is modified by using the ETILTING method to obtain anisotropic delay;

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

[0020] Particle filtering is performed on the corrected GNSS ZTD and gradient observation data to obtain the state estimation value of the particle filter and the coefficient of the linear relationship model; a linear relationship model between the observed ZTD value and the observed GNSS delay gradient parameter and the altitude is constructed;

[0021] A regular gridded sea level ZTD model is constructed based on the linear relationship model; the regular gridded sea level ZTD model is projected onto the InSAR coordinates after convolution interpolation to obtain a neutral atmosphere delay correction model;

[0022] Subtract the neutral atmosphere delay correction model from the estimated delay model at different times to obtain the delay difference suitable for the interferogram;

[0023] The interference pattern is corrected by the delay difference to obtain a corrected interference pattern.

[0024] In a possible implementation, the calculation formula of the method for estimating GNSS data using the parameter estimation method 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 zenith tropospheric dry component, E is the satellite altitude angle; ZWD is the delay of the zenith tropospheric wet component; M wet (E) is the projection function of the zenith tropospheric wet component.

[0027] In a possible implementation, the method for determining a mapping function by using a method for jointly determining mapping function coefficients based on a third-order continued fraction comprises the following steps:

[0028] The highly accurate MFlsmabc algorithm is used to estimate the coefficients a, b and c simultaneously. If the MFlsmabc algorithm reaches convergence, the modeling process ends. If it does not converge, the MFlsmab algorithm is used to build the model.

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

[0030] In a possible implementation, the method of correcting the mapping function using the ETILTING method comprises the following steps:

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

[0032] Construct a design matrix that contains functions of altitude and azimuth angles to relate the observed values ​​to the horizontal gradient parameters, and use the least squares method to estimate the horizontal gradient parameters;

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

[0034] The expression for anisotropic retardation is:

[0035]

[0036] Where ΔG represents 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, G N and G E Represents the GNSS delay gradient parameters in the north-south and east-west directions.

[0037] In one possible implementation, the linear relationship models of the observed ZTD value and the observed GNSS delay gradient parameter and the altitude are as follows:

[0038] The linear relationship model between the observed ZTD value and altitude is as follows:

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

[0040] Where ZTD obs is the observed ZTD value, ZTD0 is the sea level at the regular grid, a is the height dependence coefficient, h is the GNSS station height, and i and j represent the grid index in the east and north directions, respectively;

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

[0042]

[0043]

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

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

[0046] A regular gridded sea level ZTD model is constructed based on the 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 and is denoted as ZTD model ;

[0049] Using simple trigonometric functions of the InSAR incident angle, projecting from the zenith direction to the LOS direction, the InSAR slant path delay is calculated to obtain the neutral atmosphere delay correction model.

[0050]

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

[0052] In a second aspect, 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] An estimation module, used to estimate GNSS data using a parameter estimation method to obtain GNSS ZTD;

[0055] An anisotropic delay acquisition module is used to determine the mapping function by using a third-order continued fraction-based mapping function coefficient joint determination method; the mapping function is corrected by using the ETILTING method to obtain the anisotropic delay;

[0056] A correction module, used for correcting the GNSS ZTD by using anisotropic delay to obtain a 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 estimation value of the particle filter and the coefficient of the linear relationship model; construct the linear relationship model between the observed ZTD value and the observed GNSS delay gradient parameter and the height respectively;

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

[0059] A delay difference obtaining module is used to subtract the neutral atmosphere delay correction model from the estimated delay models at different times to obtain a delay difference suitable for the interferogram;

[0060] The correction module is used to correct the interference pattern through the delay difference to obtain the corrected interference pattern.

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

[0062] In a fourth aspect, the present invention provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, 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] In a fifth aspect, the present invention provides a computer program product, including a computer program, which, 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 the present invention are as follows:

[0065] (1) The method of the present invention shows obvious advantages in reducing the impact of atmospheric delay and improving the accuracy of monitoring tiny surface displacements.

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

[0067] (3) The method described in the present invention uses particle filtering to extract linear height dependence from the sea level ZTD distribution and projects it to InSAR coordinates for neutral atmospheric delay correction; and uses cubic convolution to interpolate GNSS ZTD data on a grid to improve the atmospheric delay correction accuracy of each pixel in the InSAR image. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 is an overall flow chart of the method of the present invention;

[0069] Figure 2Modeling flow chart for MFlsmcom method;

[0070] Figure 3 A schematic diagram of the structure of an 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 structural diagram of an electronic device provided by the present invention. DETAILED DESCRIPTION

[0072] The content of the present invention is further described below in conjunction with specific embodiments, but the content of the present invention is not limited thereto at all.

[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 a possible implementation, the GNSS data originates from the Nevada Geological Laboratory at the University of Nevada, Reno.

[0076] S200, estimating the GNSS data using a parameter estimation method to obtain a GNSS ZTD (Zenith Tropospheric Delay).

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

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

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

[0080] S300, determining a mapping function using a third-order continued fraction-based mapping function coefficient joint determination method (MFlsmcom); and modifying the mapping function using an ETILTING method to obtain anisotropic delay.

[0081] In a possible implementation, in S300, a method for determining a mapping function by jointly determining a mapping function coefficient based on a third-order continued fraction is used, such as Figure 2 As shown, the following steps are included:

[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 reaches convergence, the mapping function construction is terminated; if it does not converge, the MFlsmab algorithm is used to construct the mapping function;

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

[0085] In a possible implementation, in S300, the method of correcting the mapping function using the ETILTING method (Extended TILTING, a horizontal gradient multi-parameter modeling method taking into account high-order variation terms) 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 consists of actual observation data and contains the values ​​of anisotropic delays observed at different altitudes and azimuths), and 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, and a second-order horizontal gradient parameter is added to better fit the anisotropic delayed high-order variation term caused by local terrain disturbances and weather events. The expression of the ETILTING method is: Among them, AD(e,φ) represents 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, G N and G E 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 above-mentioned method for jointly determining the mapping function coefficients based on the third-order continued fraction (MFlsmcom), where the mathematical model of the MFlsmabc method is expressed as:

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

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

[0091]

[0092] In the formula, n represents the number of mapping factors involved in modeling, and the matrix represents that the partial differential of the third-order continued fraction on the a, b, c coefficients contains the a, b, c coefficients themselves, and the initial values ​​of the a, b, c coefficients need to be given:

[0093]

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

[0095] x is the parameter matrix to be estimated:

[0096]

[0097] l is a constant matrix:

[0098]

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

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

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

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

[0103] The construction process of the MFlsma model is similar to that of the MFlsmabc model, with the main difference being that the MFlsma model sets the values ​​of b and c to known a priori values ​​and uses the least squares method to estimate the coefficient a. If a mapping function for an elevation angle of 3° or a starting elevation angle of 3.3° is included in the application of the MFlsma model, the model is simplified to the MFfast model. The MFlsmab method only fixes the coefficient c and uses the least squares method to estimate the coefficients a and b.

[0104] On the basis of the above, the appropriate mapping function is selected according to the convergence.

[0105] S400, the GNSS ZTD is corrected using anisotropic delay to obtain a corrected GNSS ZTD. The correction formula is:

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

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

[0108] S500, performing particle filtering on the corrected GNSS ZTD and gradient observation data to obtain a state estimation value of the particle filter and obtain coefficients of a linear relationship model; constructing linear relationship models between the observed ZTD value and the observed GNSS delay gradient parameter and the altitude respectively;

[0109] The linear relationship model between the observed ZTD value and altitude is as follows:

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

[0111] Where ZTD obs is the observed ZTD value, ZTD0 is the sea level at the regular grid, a is the height dependence coefficient, h is the GNSS station height, and i and j represent the grid index in the east and north directions, respectively;

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

[0113]

[0114]

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

[0116] In a possible implementation, the state estimation value M' of the particle filter is:

[0117]

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

[0119] The state estimation value M' includes the coefficients ZTD0,a (height dependence coefficient) in the linear relationship model, H ZTD , D e and D n , these coefficients are obtained through particle filtering to construct the linear relationship between the observed ZTD value and the observed GNSS delay gradient parameter and the altitude, respectively, and then used for the correction of atmospheric delay.

[0120] In a possible implementation, the particle filter includes the following steps:

[0121] Weight calculation: For each particle, calculate its weight, which is obtained by comparing the particle-predicted observation value with the actual GNSS observation data.

[0122] Resampling: Resample according to the weight of the particles to reduce the diversity of particles and focus on more likely state areas. This step allows particles with higher weights to have more opportunities to be replicated, while particles with lower weights may be discarded.

[0123] S600, constructing a regular grid sea level ZTD model based on the linear relationship model; projecting the regular grid sea level ZTD onto the InSAR coordinates after convolution interpolation to obtain a neutral atmosphere delay correction model.

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

[0125] S610, constructing a regular gridded sea level ZTD model: constructing a regular grid to cover the entire GNSS station distribution area; each grid point represents a sea level ZTD value ZTD0, which is consistent with ZTD0 in step S500;

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

[0127] S630, adjust the gradient parameter: apply the observed GNSS delay gradient parameter G obtained in S500 N and G E The linear relationship model with the height is used to refine the model and obtain the adjusted ZTD model. These gradient parameters take into account the changes in ZTD values ​​in different directions, thereby improving the spatial resolution and accuracy of the model.

[0128] S640, convolution interpolation: performing convolution interpolation on the adjusted ZTD model to obtain the sea level ZTD (ZTD0) to smooth the discontinuity of the ZTD value and improve the continuity of the model. In this embodiment, cubic convolution interpolation is used.

[0129] S650, using ZTD0, calculates the estimated ZTD for each InSAR pixel, denoted as ZTD model ;

[0130] S660 uses simple trigonometric functions of the InSAR incident angle, projecting from the zenith direction to the LOS direction, to calculate the InSAR (STD InSAR ) Slant path delay, and obtain the neutral atmosphere delay correction model;

[0131]

[0132] Among them, STD InSAR InSAR slant path delay, ZTD model represents the estimated ZTD of each InSAR pixel and θ represents the radar incidence angle.

[0133] S700, subtracting the atmospheric delay correction model from the estimated delays at different times to obtain a delay difference suitable for the interferogram.

[0134] It should be noted that the estimated delays at different times are obtained by collecting GNSS data sources at different times. The data source of GNSS (Global Navigation Satellite System) refers to the satellite signal data used to determine the position and time of a certain point on the earth. These data sources have different times, which means that they are collected at different times, reflecting the satellite positions, signal propagation conditions, and the conditions of the earth's atmosphere (especially the troposphere) at different points in time. Therefore, the GNSS data sources are at different times, and the atmospheric delays obtained are also at different times.

[0135] S800, correcting the interference pattern by using the delay difference to obtain a corrected interference pattern.

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

[0137] The atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering provided by the present invention is described below. The atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering described below and the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering described above can be referenced to each other.

[0138] Figure 3 is a schematic diagram of the structure of an atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering provided by an embodiment of the present invention, such as 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] A collection module 31, used for GNSS data collection;

[0140] An estimation module 32, configured to estimate the GNSS data using a parameter estimation method to obtain a GNSS ZTD;

[0141] An anisotropic delay acquisition module 33 is used to determine the mapping function by using a third-order continued fraction-based mapping function coefficient joint determination method; and to modify the mapping function by using an ETILTING method to obtain the anisotropic delay;

[0142] A correction module 34, used to correct the GNSS ZTD using the anisotropic delay to obtain a 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 a state estimation value of the particle filter and obtain the coefficients of the linear relationship model; construct a linear relationship model between the observed ZTD value and the observed GNSS delay gradient parameter and the height respectively;

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

[0145] A delay difference obtaining module 37, used for subtracting the neutral atmosphere delay correction model from the estimated delay models at different times to obtain a delay difference suitable for the interferogram;

[0146] The correction module 38 is used to correct the interference pattern by using the delay difference to obtain a corrected interference pattern.

[0147] Figure 4 An example of a physical structure diagram of an electronic device is shown in FIG. Figure 4 As shown, the electronic device may include: a processor 410, a communication interface 420, a memory 430 and a communication bus 440, wherein the processor 410, the communication interface 420 and the memory 430 communicate with each other through the communication bus 440. The processor 410 may call the logic instructions in the memory 430 to execute the atmospheric delay correction method based on the GNSS high-order horizontal gradient and particle filtering.

[0148] In addition, the logic instructions in the above-mentioned memory 430 can be implemented in the form of a software functional unit and can be stored in a computer-readable storage medium when it is sold or used as an independent product. Based on such an understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product, and the computer software product is stored in a storage medium, including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk and other media that can store program codes.

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

[0150] On the other hand, 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 execute the atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering provided by the above-mentioned methods.

[0151] The device embodiments described above are merely illustrative, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. Ordinary technicians in this field can understand and implement it without paying creative labor.

[0152] Through the description of the above implementation methods, those skilled in the art can clearly understand that each implementation method can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solution is essentially or the part that contributes to the prior art can be embodied in the form of a software product, and the computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a disk, an optical disk, etc., including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment 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, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions 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: The following steps are involved: GNSS data collection; Use parameter estimation method to estimate GNSS data and obtain GNSS ZTD; Determine the mapping function by using a third-order continued fraction-based joint determination method for mapping function coefficients; The mapping function is modified using the ETILTING method to obtain anisotropic delay; The GNSS ZTD is corrected by using the anisotropic delay to obtain the corrected GNSS ZTD; Particle filtering is performed on the corrected GNSS ZTD and gradient observation data to obtain the state estimation value of the particle filter and the coefficient of the linear relationship model; a linear relationship model between the observed ZTD value and the observed GNSS delay gradient parameter and the altitude is constructed; A regular gridded sea level ZTD model is constructed based on the linear relationship model; the regular gridded sea level ZTD model is projected onto the InSAR coordinates after convolution interpolation to obtain a neutral atmosphere delay correction model; Subtract the neutral atmosphere delay correction model from the estimated delay model at different times to obtain the delay difference suitable for the interferogram; The interference pattern is corrected by the delay difference to obtain a corrected interference pattern.

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

3. The atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering according to claim 1 is characterized in that: The method for determining a mapping function by using a method for jointly determining mapping function coefficients based on a third-order continued fraction comprises the following steps: The highly accurate MFlsmabc algorithm is used to estimate the coefficients a, b and c simultaneously. If the MFlsmabc algorithm reaches convergence, the modeling process ends. If it does not converge, the MFlsmab algorithm is used to build the model. Evaluate the convergence status of the MFlsmab algorithm and the MFlsma algorithm; if neither converges, use the MFfast algorithm for calculation.

4. The atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering according to claim 1 is characterized in that: The method of using the ETILTING method to correct the mapping function to obtain anisotropic delay comprises the following steps: Calculate anisotropic delay at different azimuth and elevation angles; Construct a design matrix that contains functions of altitude and azimuth angles to relate the observed values ​​to the horizontal gradient parameters, and use the least squares method to estimate the horizontal gradient parameters; The model expression of the ETILTING method is determined to modify the mapping function, solve the horizontal gradient parameters, including the first-order and second-order horizontal gradient parameters in the north-south and east-west directions, and obtain the anisotropic delay; The expression for anisotropic retardation is: Where ΔG represents 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, G N and G E Represents the GNSS delay gradient parameters in the north-south and east-west directions.

5. The atmospheric delay correction method based on GNSS high-order horizontal gradient and particle filtering according to claim 1, characterized in that: The linear relationships between the observed ZTD value and the observed GNSS delay gradient parameter and the altitude are as follows: The observed linear relationship between ZTD value and altitude is: ZTD obs =ZTD0(i,j)+a·h(i,j) Where ZTD obs is the observed ZTD value, ZTD0 is the sea level at the regular grid, a is the height dependence coefficient, h is the GNSS station height, and i and j represent the grid index in the east and north directions, respectively; Observed GNSS delay gradient parameter G N and G E Linear relationship with height: Among them, H ZTD Denotes the scale height of ZTD, D e Indicates the grid spacing in the east-west direction, D n Indicates the grid spacing in the north-south direction.

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

7. An atmospheric delay correction device based on GNSS high-order horizontal gradient and particle filtering, characterized in that: include: Acquisition module, used for GNSS data acquisition; An estimation module, used to estimate GNSS data using a parameter estimation method to obtain GNSS ZTD; An anisotropic delay acquisition module is used to determine the mapping function by using a third-order continued fraction-based mapping function coefficient joint determination method; the mapping function is corrected by using the ETILTING method to obtain the anisotropic delay; A correction module, used for correcting the GNSS ZTD by using anisotropic delay to obtain a 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 estimation value of the particle filter and the coefficient of the linear relationship model; construct the linear relationship model between the observed ZTD value and the observed GNSS delay gradient parameter and the height respectively; The model acquisition module is used to construct a regular grid sea level ZTD model based on the linear relationship model; the regular grid sea level ZTD model is projected onto the InSAR coordinates after convolution interpolation to obtain a neutral atmosphere delay correction model; A delay difference obtaining module is used to subtract the neutral atmosphere delay correction model from the estimated delay models at different times to obtain a delay difference suitable for the interferogram; The correction module is used to correct the interference pattern through the delay difference to obtain the corrected interference pattern.

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

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

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

Images