An InSAR Turbulent Atmosphere Correction Method Based on a 3D Stochastic Model

Abstract

This application is applicable to the field of geodetic surveying technology, and provides an InSAR turbulent atmosphere correction method based on a three-dimensional random model, including: obtaining M interferograms corresponding to N SAR images of a research area, and dividing each interferogram into Y windows by using a window of a preset size, and taking the phase information within the y-th window of all interferograms as the y-th observation window; respectively for each observation window, constructing a spatio-temporal three-dimensional random model within the observation window according to M interferograms, constructing a functional model for all pixels within the observation window, and using the weighted least squares method to solve to obtain the deformation parameters corresponding to the observation window; splicing the deformation parameters corresponding to all observation windows, and calculating the deformation phase of the research area based on the splicing result. This application can effectively suppress the influence of turbulent atmosphere.

Description

Technical Field

[0001] This application belongs to the technical field of geodetic surveying, and particularly relates to an InSAR turbulent atmosphere correction method based on a three-dimensional stochastic model. Background Art

[0002] Interferometric Synthetic Aperture Radar (InSAR) is a new generation of space geodetic surveying technology. With its advantages such as all-weather observation ability and high spatial resolution, it shows unique value in the field of surface deformation monitoring. This technology has been successfully applied in research fields such as earthquakes, volcanoes, and landslides, accumulating rich scientific achievements. However, InSAR technology is interfered by tropospheric delay, which not only significantly reduces the accuracy of surface deformation monitoring but also leads to great uncertainty in the interpretation of the deformation field, becoming a bottleneck restricting the wide application of this technology in high-precision surface deformation monitoring.

[0003] Tropospheric delay can be divided into vertical stratified delay and turbulent atmosphere delay according to its physical properties. The vertical stratified delay shows a strong correlation with topography in space and is generally corrected by modeling it as a function of topography. The spatio-temporal variation of turbulent atmosphere delay is relatively complex, and it is difficult to establish an exact function model for it. Generally, three methods, namely spatio-temporal filtering, interferogram stacking, and stochastic model, are used to suppress turbulent atmosphere delay. The spatio-temporal filtering method designs a corresponding filter according to the significant differences in spatio-temporal characteristics between turbulent atmosphere delay and deformation signals to effectively separate the two. However, the core parameter of this method, the filtering window, is subjectively determined depending on user experience. The interferogram stacking method believes that turbulent atmosphere is randomly distributed in time. By differentiating and then stacking interferograms with the same time interval, the turbulent atmosphere of SAR images is obtained and the atmospheric influence is suppressed. However, this method is only applicable to linear deformation. The stochastic model method suppresses the atmospheric influence by establishing the variance-covariance matrix of turbulent atmosphere. However, the existing stochastic model only considers the temporal correlation of turbulent atmosphere in the interferogram and establishes the variance-covariance matrix of the temporal dimension of turbulent atmosphere pixel by pixel, without considering the spatial correlation of turbulent atmosphere.

[0004] In summary, the current methods for suppressing turbulent atmosphere delay are difficult to effectively suppress the influence of turbulent atmosphere. Summary of the Invention

[0005] The embodiments of this application provide an InSAR turbulent atmosphere correction method based on a three-dimensional stochastic model, which can solve the problem of being difficult to effectively suppress the influence of turbulent atmosphere.

[0006] The embodiments of this application provide an InSAR turbulent atmosphere correction method based on a three-dimensional stochastic model, including:

[0007] Obtain M interferograms corresponding to N SAR images of the research area, and divide each interferogram into Y windows using a window of a preset size. Take the phase information within the y-th window of all interferograms as the y-th observation window; ;

[0008] For each observation window respectively, construct a spatio-temporal three-dimensional stochastic model within the observation window based on the M interferograms; the spatio-temporal three-dimensional stochastic model is used to describe the spatio-temporal characteristics of turbulent atmosphere and decorrelation noise within the observation window;

[0009] For each observation window respectively, construct a functional model for all pixels within the observation window; the functional model is used to describe the relationship between the phase information within the observation window and the deformation parameters to be solved;

[0010] For each observation window respectively, solve the deformation parameters using the weighted least squares method based on the functional model and the spatio-temporal three-dimensional stochastic model corresponding to the observation window, and obtain the numerical values of the deformation parameters under the observation window;

[0011] Stitch the numerical values of the deformation parameters under all observation windows, and calculate the deformation phase of the research area based on the stitching result.

[0012] Optionally, constructing a spatio-temporal three-dimensional stochastic model within the observation window based on the M interferograms includes:

[0013] Based on the M interferograms, construct the spatial covariance matrix of turbulent atmosphere within the observation window for each interferogram;

[0014] Based on the constructed spatial covariance matrix of turbulent atmosphere, determine the spatial covariance matrix of turbulent atmosphere within the observation window for each SAR image scene;

[0015] Based on the spatial covariance matrices of turbulent atmosphere within the observation window for N SAR image scenes, determine the spatio-temporal three-dimensional covariance matrix of turbulent atmosphere within the observation window for the M interferograms;

[0016] Based on the spatio-temporal three-dimensional covariance matrix of turbulent atmosphere, determine the spatio-temporal three-dimensional stochastic model within the observation window.

[0017] Optionally, the spatial covariance matrix of turbulent atmosphere within the observation window for the -th interferogram is:

[0018] ;

[0019] wherein, , the element in the spatial covariance matrix of turbulent atmosphere , represents the position within the observation window​ and the distance between two pixels at the position is , , the size of the observation window is , representing the number of rows of the observation window, representing the number of columns of the observation window;

[0020]

[0021] represents the sill value of the th interferogram, represents the range value of the th interferogram.

[0022] Optionally, before the step of constructing the turbulent atmospheric spatial covariance matrix of each interferogram within the observation window based on M interferograms, the InSAR turbulent atmospheric correction method further includes:

[0023] Performing theoretical variogram fitting on each interferogram using the spherical model to obtain the sill value and range value corresponding to each interferogram.

[0024] Optionally, based on the constructed turbulent atmospheric spatial covariance matrix, determining the turbulent atmospheric spatial covariance matrix of each SAR image within the observation window includes:

[0025] Calculating the turbulent atmospheric spatial covariance matrix of the th SAR image within the observation window through the following formula :

[0026] ;

[0027]

[0028]

[0029]

[0030]

[0031] where, represents the coefficient matrix. The coefficient matrix B is an M×N coefficient matrix, which is determined by the relationship between the interferogram and the SAR image. For each row of it, the position of the master image is 1, the position of the slave image is also 1, and the remaining positions are 0. .

[0032] Optionally, based on the turbulent atmosphere spatial covariance matrix of N scene SAR images within the observation window, determine the three-dimensional spatio-temporal covariance matrix of turbulent atmosphere for M interferograms within the observation window, including:

[0033] Calculate the three-dimensional spatio-temporal covariance matrix of turbulent atmosphere for M interferograms within the observation window through the following formula :

[0034]

[0035]

[0036] Wherein, represents the coefficient matrix, , represents the coefficient matrix connecting the interferogram and the SAR image. The coefficient matrix F is an M×N coefficient matrix, which is determined by the relationship between the interferogram and the SAR image. For each row of it, the position of the master image is -1, the position of the slave image is 1, and the remaining positions are 0, represents the Kronecker product operation, represents a unit matrix of size, , represents the turbulent atmosphere spatial covariance matrix of all SAR images within the observation window;

[0037] .

[0038] Optionally, based on the three-dimensional spatio-temporal covariance matrix of turbulent atmosphere, determine the three-dimensional spatio-temporal random model within the observation window, including:

[0039] Calculate the three-dimensional spatio-temporal random model through the following formula :

[0040]

[0041]

[0042]

[0043] Wherein, represents the decorrelation noise variance matrix of all pixels of the interferograms within the observation window, represents the decorrelation noise variance matrix of the pixels of the th interferogram within the observation window, The element in is the variance of the decorrelation noise of the pixel , represents the number of pixels within the observation window, = , , Represents pixel The coherence coefficient of , L represents the number of multi-looks.

[0044] Optionally, the function model is:

[0045]

[0046]

[0047]

[0048]

[0049]

[0050] in, Indicates that M interference patterns are within the observation window The phase vector of pixels, represents the coefficient matrix, represents the deformation parameter to be solved within the observation window, represents the residual phase, represents an M×4 coefficient matrix, represents the Kronecker product operation, Indicates a The identity matrix of size, , , represents the number of pixels in the observation window, Represents the phase information corresponding to pixel p in all interference patterns, It represents the residual phase of pixel p in all interference patterns. Indicates The cumulative time of the first scene SAR image relative to the first scene SAR image, represents the angle of incidence, represents the slant distance from the satellite to the surface of the earth, Indicates Scene SAR image and the The vertical baseline length of the interferogram generated by the SAR image, , , .

[0051] Optionally, the number of deformation parameters is multiple, and the values of the deformation parameters in all observation windows are concatenated, including:

[0052] For each deformation parameter, calculate the mean value of the deformation parameter values under all observation windows, and use the calculated mean value as the final value of the deformation parameter;

[0053] Use the final values of all deformation parameters as the splicing result.

[0054] The above solution of this application has the following beneficial effects:

[0055] In the embodiment of this application, since the constructed spatio-temporal three-dimensional random model can express the spatio-temporal characteristics of turbulent atmosphere, when solving the deformation parameters for deformation phase calculation based on the spatio-temporal three-dimensional random model and the function model of pixels, the influence of turbulent atmosphere can be effectively suppressed, and thus the InSAR surface deformation monitoring accuracy can be effectively improved.

[0056] In addition, since the construction of the spatio-temporal three-dimensional random model and the solution of deformation parameters are carried out window by window, while greatly reducing the operation difficulty, the accuracy of the spatio-temporal three-dimensional random model and deformation parameters is improved, which is convenient for effectively suppressing the influence of turbulent atmosphere.

[0057] Other beneficial effects of this application will be described in detail in the subsequent specific implementation part. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] In order to more clearly illustrate the technical solutions in the embodiments of this application, the following will briefly introduce the drawings required for the embodiments or the description of the prior art. Obviously, the following drawings are only some embodiments of this application. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0059] Figure 1 It is a flowchart of the InSAR turbulent atmosphere correction method based on a three-dimensional random model provided by an embodiment of this application;

[0060] Figure 2 It is a comparison diagram of the deformation time series solved by the GPS site SACY, the ordinary least squares method and the method of this application in an example of this application;

[0061] Figure 3 It is a comparison diagram of the deformation time series solved by the GPS site AZU1, the ordinary least squares method and the method of this application in an example of this application;

[0062] Figure 4 It is a comparison diagram of the deformation time series solved by the GPS site BLSA, the ordinary least squares method and the method of this application in an example of this application;

[0063] Figure 5In an example of this application, it is a comparison chart of deformation time series solved at the GPS station CAS4, by the ordinary least squares method, and by the method of this application. Detailed implementation manners

[0064] In the following description, specific details such as specific system structures and technologies are presented for the purpose of illustration rather than limitation, so as to thoroughly understand the embodiments of this application. However, those skilled in the art should clearly understand that this application can also be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to avoid unnecessary details from interfering with the description of this application.

[0065] It should be understood that when used in the specification of this application and the appended claims, the term "comprising" indicates the presence of the described features, wholes, steps, operations, elements, and / or components, but does not exclude the presence or addition of one or more other features, wholes, steps, operations, elements, components, and / or their combinations.

[0066] It should also be understood that the term "and / or" used in the specification of this application and the appended claims refers to any combination and all possible combinations of one or more of the associated listed items, and includes these combinations.

[0067] As used in the specification of this application and the appended claims, the term "if" can be interpreted as "when", "once", "in response to determining", or "in response to detecting" according to the context. Similarly, the phrase "if determined" or "if [the described condition or event] is detected" can be interpreted as meaning "once determined", "in response to determining", "once [the described condition or event] is detected", or "in response to detecting [the described condition or event]" according to the context.

[0068] In addition, in the description of the specification of this application and the appended claims, the terms "first", "second", "third", etc. are only used for differentiating descriptions and cannot be understood as indicating or implying relative importance.

[0069] The reference to "one embodiment" or "some embodiments" etc. described in the specification of this application means that specific features, structures, or characteristics described in combination with that embodiment are included in one or more embodiments of this application. Thus, statements such as "in one embodiment", "in some embodiments", "in other some embodiments", "in still other embodiments", etc. that appear in different places in this specification do not necessarily refer to the same embodiment, but mean "one or more but not all embodiments", unless otherwise specifically emphasized in other ways. The terms "comprising", "including", "having", and their variants all mean "including but not limited to", unless otherwise specifically emphasized in other ways.

[0070] In view of the problem that it is currently difficult to effectively suppress the influence of turbulent atmosphere, an embodiment of the present application provides an InSAR turbulent atmosphere correction method based on a three-dimensional random model. By constructing a spatio-temporal three-dimensional random model window by window and solving deformation parameters, the deformation phase of the study area is determined. Since the constructed spatio-temporal three-dimensional random model can express the spatio-temporal characteristics of the turbulent atmosphere, when solving the deformation parameters for deformation phase calculation based on the spatio-temporal three-dimensional random model and the functional model of pixels, the influence of the turbulent atmosphere can be effectively suppressed, thereby effectively improving the accuracy of InSAR surface deformation monitoring.

[0071] In addition, since the construction of the spatio-temporal three-dimensional random model and the solution of deformation parameters are carried out window by window, while greatly reducing the computational difficulty, the accuracy of the spatio-temporal three-dimensional random model and deformation parameters is improved, which is convenient for effectively suppressing the influence of the turbulent atmosphere.

[0072] The following uses specific embodiments to exemplarily illustrate the InSAR turbulent atmosphere correction method based on a three-dimensional random model provided by the present application.

[0073] As Figure 1 shown, the InSAR turbulent atmosphere correction method based on a three-dimensional random model provided by an embodiment of the present application includes the following steps:

[0074] Step 11, obtain M interferograms corresponding to N SAR images of the study area, and divide each interferogram into Y windows by using a window of a preset size. The phase information within the y-th window of all interferograms is used as the y-th observation window.

[0075] The above ; the above study area is the area where surface deformation monitoring is required. The synthetic aperture radar (SAR) images can be collected by the SAR sensor carried by the satellite. After collecting N SAR images of the study area, the interferograms can be obtained by processing the SAR images through the common method of obtaining interferograms. That is, after selecting a master image, register the remaining SAR images, and then set the spatio-temporal baseline threshold for differential processing to obtain M interferograms. After removing the flat-earth phase and terrain phase according to the external digital elevation model (DEM) data and performing phase unwrapping, finally obtain M unwrapped interferograms. It can be understood that the interferograms in the above step 11 are unwrapped interferograms.

[0076] The size of the above window is , represents the number of rows of the window, Indicates the number of columns of the window. After dividing each interferogram, the windows at the same position in all interferograms can be used as an observation window. That is, the phase information within the y-th window of all interferograms is used as the y-th observation window. It can be understood that in order to achieve a smooth transition between windows, a certain overlap rate can be set between windows. Generally, the window size can be set to , and the moving step size of the window can be set to 24.

[0077] Step 12: For each observation window respectively, construct a spatio-temporal three-dimensional random model within the observation window based on M interferograms.

[0078] The above spatio-temporal three-dimensional random model is used to describe the spatio-temporal characteristics of turbulent atmosphere and decoherence noise within the observation window.

[0079] Step 13: For each observation window respectively, construct a function model for all pixels within the observation window.

[0080] The above function model is used to describe the relationship between the phase information within the observation window and the deformation parameters to be solved. The deformation parameters include average deformation rate, average acceleration, average acceleration change rate, and terrain residual.

[0081] Step 14: For each observation window respectively, solve the deformation parameters using the weighted least squares method based on the function model and spatio-temporal three-dimensional random model corresponding to the observation window, and obtain the numerical values of the deformation parameters under the observation window.

[0082] That is, by solving the deformation parameters using the weighted least squares method based on the function model and spatio-temporal three-dimensional random model corresponding to the observation window, obtain the numerical values of the deformation parameters to be solved in Step 13, and use the obtained numerical values of the deformation parameters as the numerical values of the deformation parameters under the observation window.

[0083] Step 15: Stitch the numerical values of the deformation parameters under all observation windows, and calculate the deformation phase of the study area based on the stitching result.

[0084] In some embodiments of the present application, the number of deformation parameters is multiple. When stitching, for each deformation parameter respectively, calculate the mean value of the numerical values of the deformation parameters under all observation windows, and use the calculated mean value as the final value of the deformation parameter; then use the final values of all deformation parameters as the stitching result.

[0085] Exemplarily, for the average deformation rate, the mean value of the numerical values of the average deformation rate under all observation windows can be used as the final value of the average deformation rate.

[0086] It should be noted that after obtaining the final values of all deformation parameters, the corresponding low-frequency deformation phase in time domain can be obtained. Then, the residual phase is filtered by a low-pass filter in time domain, and the filtered phase is added to the low-frequency deformation phase in time domain to obtain the deformation phase of the study area, which is convenient for deformation monitoring.

[0087] It is worth mentioning that since the constructed spatio-temporal three-dimensional stochastic model can express the spatio-temporal characteristics of turbulent atmosphere, when solving the deformation parameters for calculating the deformation phase based on the spatio-temporal three-dimensional stochastic model and the functional model of pixels, the influence of turbulent atmosphere can be effectively suppressed, and thus the accuracy of InSAR surface deformation monitoring can be effectively improved.

[0088] In addition, since the construction of the spatio-temporal three-dimensional stochastic model and the solution of deformation parameters are carried out window by window, the operation difficulty is greatly reduced, and at the same time, the accuracy of the spatio-temporal three-dimensional stochastic model and deformation parameters is improved, which is convenient for effectively suppressing the influence of turbulent atmosphere.

[0089] The following will exemplarily illustrate the specific implementation method of constructing the spatio-temporal three-dimensional stochastic model within the observation window according to M interferograms in step 12 in combination with specific embodiments.

[0090] Specifically, the specific implementation method of constructing the spatio-temporal three-dimensional stochastic model within the observation window according to M interferograms in step 12 includes the following steps:

[0091] Step 12.1: According to M interferograms, construct the spatio-temporal covariance matrix of turbulent atmosphere for each interferogram within the observation window.

[0092] The above spatio-temporal covariance matrix of turbulent atmosphere can also be called the spatio-temporal variance-covariance matrix of turbulent atmosphere.

[0093] In some embodiments of the present application, the spatio-temporal covariance matrix of turbulent atmosphere for the th interferogram within the observation window is:

[0094] ;

[0095] where , the elements on the main diagonal of the spatio-temporal covariance matrix of turbulent atmosphere represent the variance of a certain pixel within the observation window, and the elements off the main diagonal represent the covariance between two pixels within the observation window. The element in the spatio-temporal covariance matrix of turbulent atmosphere , represents the distance between two pixels at positions and within the observation window, , , the size of the observation window is , represents the number of rows of the observation window, represents the number of columns of the observation window.

[0096]

[0097] Among them, represents the sill value of the th interferogram, represents the range value of the th interferogram.

[0098] It can be understood that before the step of constructing the turbulent atmosphere spatial covariance matrix of each interferogram within the observation window according to M interferograms, the above InSAR turbulent atmosphere correction method further includes the following steps of determining the sill value and the range value: fitting the theoretical variogram of each interferogram using the spherical model to obtain the sill value and the range value corresponding to each interferogram.

[0099] In some embodiments of the present application, the process of obtaining the sill value of the th interferogram and the range value of the th interferogram is as follows:

[0100] The InSAR phase components are:

[0101] (1)

[0102] is the unwrapped interferometric phase, , , , and are the deformation phase, the topographic residual phase, the vertically stratified atmosphere phase, the turbulent atmosphere phase and the noise respectively. Assuming that the vertically stratified atmosphere has been well removed through the phase-elevation model, for interferograms with short spatio-temporal baselines, the influence of the deformation phase and the topographic residual phase is very small, and the noise is of a smaller magnitude compared to the turbulent atmosphere. Therefore, the main component of the phase is the turbulent atmosphere. The turbulent atmosphere has significant spatial correlation. The variogram, as a core tool in geostatistics, is often used to analyze the spatial correlation of data. For spatial data , the variogram is defined as:

[0103] (2)

[0104] is the variogram value at the distance , is the number of data point pairs with a distance of . Substitute the phase of the th interferogram as the spatial data into Equation (2), set different distances to obtain the corresponding variogram values . is the number of at the distance

[0105] Since the number of samples is always limited, the experimental variogram (i.e., the variogram value) obtained above is not continuous in space. To obtain a spatially continuous variogram, it is also necessary to fit the experimental variogram with a theoretical variogram. In this application, the following spherical model is used for fitting:

[0106] (3)

[0107] is the maximum variogram value, called the sill value, is the maximum correlation distance, called the range value. After fitting the theoretical variogram for each interferogram (i.e., fitting the distance and the corresponding variogram value ), the corresponding sill value and range value are obtained, where M is the number of interferograms.

[0108] Step 12.2: Based on the constructed spatial covariance matrix of the turbulent atmosphere, determine the spatial covariance matrix of the turbulent atmosphere within the observation window for each SAR image scene.

[0109] The above spatial covariance matrix of the turbulent atmosphere can also be referred to as the spatial variance-covariance matrix of the turbulent atmosphere.

[0110] In some embodiments of this application, the spatial covariance matrix of the turbulent atmosphere within the observation window for the th SAR image scene can be calculated by the following formula :

[0111] ;

[0112]

[0113]

[0114]

[0115]

[0116] Among them, , represents the coefficient matrix. The coefficient matrix B is an M×N coefficient matrix, which is determined by the relationship between the interferogram and the SAR image. For each row of it, the position of the master image is 1, the position of the slave image is also 1, and the remaining positions are 0.

[0117] Among them, is obtained as follows:

[0118] Suppose the th interferogram is obtained by differencing the SAR images j and k (i.e., the th scene SAR image and the th scene SAR image). For each pixel in the observation window, according to the variance-covariance propagation law, there is:

[0119]

[0120] is the covariance between two pixels at positions and in the th interferogram in the observation window (when b = d and c = e, represents the variance). is the variance-covariance corresponding to the th scene SAR image, is the variance-covariance corresponding to the th scene SAR image, and N is the number of SAR images.

[0121] Next, rewrite into the following matrix form:

[0122] ;

[0123]

[0124]

[0125] is the variance-covariance vector of the interferogram corresponding to a pixel, is the variance-covariance vector of the SAR image corresponding to a pixel. The variance-covariance of the corresponding SAR image of this pixel is obtained by least squares solution:

[0126]

[0127] Step 12.3: Based on the turbulent atmosphere spatial covariance matrix of N scene SAR images in the observation window, determine the turbulent atmosphere spatio-temporal three-dimensional covariance matrix of M interferograms in the observation window.

[0128] In some embodiments of the present application, the three-dimensional spatio-temporal covariance matrix of the turbulent atmosphere within the observation window for M interferograms can be specifically calculated through the following formula :

[0129]

[0130]

[0131] where represents the coefficient matrix, , represents the coefficient matrix connecting the interferogram and the SAR image. The coefficient matrix F is an M×N coefficient matrix, which is determined by the relationship between the interferogram and the SAR image. For each row of it, the position of the master image is -1, the position of the slave image is 1, and the remaining positions are 0. represents the Kronecker product operation, represents an identity matrix of size, , represents the three-dimensional spatial covariance matrix of the turbulent atmosphere of all SAR images within the observation window;

[0132] .

[0133] Step 12.4: Determine the three-dimensional spatio-temporal random model within the observation window based on the three-dimensional spatio-temporal covariance matrix of the turbulent atmosphere.

[0134] In some embodiments of the present application, the three-dimensional spatio-temporal random model can be specifically calculated through the following formula :

[0135]

[0136]

[0137]

[0138] where represents the decorrelation noise variance matrix of all pixels of the interferograms within the observation window, represents the decorrelation noise variance matrix of the pixels of the th interferogram within the observation window, The element in is the variance of the decorrelation noise of the pixel , represents the number of pixels within the observation window, = , , represents the coherence coefficient of a pixel point , and L represents the number of looks.

[0139] In some embodiments of the present application, the function model constructed in step 13 is:

[0140]

[0141]

[0142]

[0143]

[0144]

[0145] where G is an M×4 coefficient matrix, and the elements of the th row are composed of the parameters of the th interferogram, represents the phase vector of M interferograms at pixels within the observation window, represents the coefficient matrix, represents the deformation parameter to be solved within the observation window, represents the residual phase, represents the Kronecker product operation, represents an identity matrix of size, , , represents the number of pixels within the observation window, represents the phase information corresponding to all interferograms of pixel p, , represents the phase information corresponding to pixel p in the th interferogram, represents the residual phase corresponding to pixel p in all interferograms, represents the cumulative time of the th scene SAR image relative to the first scene SAR image, is 0, represents the incidence angle, represents the th scene SAR image and the th scene SAR image generated by the vertical baseline length of the interferogram, , , .

[0146] An exemplary description of the construction process of the function model is as follows:

[0147] The Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) method often models the low-frequency deformation over time and the topographic residual phase. For the interferogram formed by the j-th and k-th SAR images, the low-frequency deformation over time and the topographic residual of pixel p in the observation window are as follows:

[0148]

[0149] is the low-frequency deformation over time and the topographic residual phase for modeling, is the radar wavelength, 、 and are the average deformation rate, average acceleration, and average acceleration change rate, respectively, is the vertical baseline corresponding to the interferogram, is the slant range from the satellite to the ground surface, is the incident angle, is the topographic residual. and are the cumulative times of the j-th and k-th SAR images relative to the first SAR image, respectively. Considering there are M interferograms, written in matrix form is:

[0150] ;

[0151] is the phase corresponding to pixel p in all interferograms, is an M×4 coefficient matrix, is the low-frequency deformation parameter over time and the topographic residual to be obtained for pixel p, is the residual phase, mainly composed of the unmodeled deformation and the turbulent atmosphere phase. The above only establishes a function model for one pixel. Since the method of this application is solved window by window, a function model for all pixels within a window also needs to be established:

[0152]

[0153]

[0154] It should be noted that for a new generation of satellites (such as Sentinel-1), the spatial baseline is controlled very short (usually less than 200 m), and the influence of the topographic residual is greatly weakened. Therefore, for Sentinel-1 satellite data, the influence of the topographic residual can often be ignored, and the above function model can be further simplified.

[0155] In some embodiments of the present application, the specific process of solving the deformation parameters by using the weighted least squares method based on the function model corresponding to the observation window and the spatio-temporal three-dimensional stochastic model in step 14 is as follows:

[0156] By using the weighted least squares method to solve the function model corresponding to the observation window and the spatio-temporal three-dimensional stochastic model, the following can be obtained:

[0157]

[0158] It should be noted that in practical applications, the results calculated by the above formula are the values of the average deformation rate, average acceleration, average acceleration change rate, and topographic residual within the observation window.

[0159] According to the solved low-frequency deformation parameters in time, the low-frequency deformation phase in time can be obtained. Since there are a small amount of unmodeled deformation phases in the residual phase, they can be extracted by time low-pass filtering. At this time, the proportion of the deformation phase in the residual phase is relatively low, so a large window filtering can be used to extract the deformation phase. Adding the low-frequency deformation phase in time to the non-model deformation phase filtered and extracted from the residual phase can obtain the final deformation phase within an observation window.

[0160] The InSAR turbulent atmosphere correction method of the present application will be exemplarily described below with specific examples.

[0161] In this example, experiments were carried out in a certain area using Sentinel-1 satellite data of the European Space Agency. There are dense Global Positioning System (GPS) stations in this area, which can be used to evaluate the InSAR deformation results. The coverage time of the Sentinel data is from January 10, 2018 to December 24, 2018, with a total of 30 SAR images. A spatial baseline threshold of 150 m and a temporal baseline threshold of 60 days were set, and a total of 109 interferograms were generated. In order to suppress noise, the interferograms were multi-looked in the range and azimuth directions with 20×5, and the spatial resolution is about 80 m. The minimum cost flow method was used for phase unwrapping.

[0162] 1) Construct the spatio-temporal three-dimensional stochastic model within the observation window

[0163] In order to reduce the influence of the deformation phase and the topographic residual phase, the experimental variogram was obtained by using the short spatio-temporal baseline interferograms generated above. Since the method proposed in the present application mainly targets the turbulent atmosphere, in order to avoid the influence of the vertically stratified atmosphere, the area with an elevation higher than 250 m was masked. The step size was set to 160 m (corresponding to 2 pixels), and a series of experimental variogram values were solved. The theoretical variogram of the spherical model was used to fit the experimental variogram to obtain the corresponding sill value and range value . The spatial window size was set to 25×25 pixels, and a spatio-temporal three-dimensional random model was generated.

[0164] 2) Solving within the observation window

[0165] Due to the good orbit control of Sentinel satellites, the spatial baselines of all interferograms are less than 150 m. Therefore, compared with deformation and turbulent atmosphere, the influence of topographic residuals can be ignored. Thus, only the time-low-frequency deformation phase is modeled functionally, and the time-low-frequency deformation parameters within an observation window are solved using the weighted least squares method. Then, the moving step size is set to 24 pixels, and the observation window is moved for solution. The overlapping area between observation windows is averaged. After that, the time-low-frequency deformation phase is solved based on the time-low-frequency deformation parameters. Finally, the non-model deformation phase is extracted from the residual phase using time low-pass filtering. The time low-pass filtering window is set to 180 days, and the final deformation phase can be obtained by adding the time-low-frequency deformation and the non-model deformation.

[0166] The ordinary least squares method and the method of this application are used to solve the InSAR line-of-sight (LOS) deformation rate, and the deformation field accuracies of the two methods are evaluated through a series of indicators.

[0167] To evaluate the uncertainties of the deformation fields of the two methods, the following indicators are used for evaluation:

[0168] ;

[0169] is the epoch corresponding displacement, is the predicted linear displacement, is the number of SAR images, is the cumulative time of the u-th SAR image relative to the first SAR image, represents the uncertainty caused by time random residual noise in the time series displacement, is the mean of the cumulative time.

[0170] By calculating the uncertainties of the two methods, the uncertainty of the method of this application has been significantly reduced compared with the ordinary least squares method. The average uncertainties of the deformation fields of the method of this application and the ordinary least squares method are 1.7 mm / yr and 0.8 mm / yr respectively, which indicates that the internal consistency accuracy of the deformation rate has been improved by about 52.9%.

[0171] After evaluating the internal symbol accuracy of the complete shape change rate, the external accuracy of the shape change rate is evaluated using external GPS data. First, the GPS three-dimensional deformation is projected onto the LOS direction deformation using the imaging geometry of the SAR satellite. Then, for each GPS site, InSAR pixel points within 200 m nearby are searched, and the average deformation of these pixel points is taken as the InSAR deformation corresponding to the site. The comparison results of the GPS shape change rate, the shape change rate solved by the ordinary least squares method, and the shape change rate solved by the method of this application are shown in Table 1. At most GPS sites, the shape change rate of the method of this application is closer to the GPS shape change rate. The RMSE values of the shape change rate of the ordinary least squares method and the method of this application are 8.01 mm / yr and 6.59 mm / yr respectively, indicating that the external symbol accuracy of the shape change rate has been improved by approximately 17.73%.

[0172] Table 1 Comparison results of shape change rates

[0173]

[0174] It is also possible to verify that the deformation time series obtained by the method of this application better resists the influence of turbulent atmosphere by analyzing the deformation time series solved by the two methods. Specifically, the deformation time series solved by the two methods are quantitatively evaluated using the GPS deformation time series, and the deformation time series at the GPS sites are plotted as Figures 2 to 5 shown. Compared with the ordinary least squares method, the deformation time series solved by the method of this application is more consistent with the GPS deformation time series, with higher consistency. Among them, Figures 2 to 5 the abscissas of both represent time, and the ordinates both represent displacement, Figures 2 to 5 respectively showing the deformation time series at GPS sites SACY, AZU1, BLSA, and CAS4.

[0175] Finally, the root mean square error of the deformation time series of each GPS site is statistically calculated. The comparison results of the root mean square error of the deformation time series solved by the ordinary least squares method and the method of this application are shown in Table 2. The average RMSE values of the deformation time series of the ordinary least squares method and the method of this application are 7.83 mm and 4.26 mm respectively, indicating that the external symbol accuracy of the deformation time series has been improved by approximately 45.59%.

[0176] Table 2 Comparison results of root mean square error values of deformation time series

[0177]

[0178] From the above experimental data, it can be seen that the average uncertainties of the shape change rates solved by the ordinary least squares method and the method proposed in this application are 1.7 mm / yr and 0.8 mm / yr respectively, and the internal symbol accuracy of the shape change rate has been improved by approximately 52.9%.

[0179] In the instance verification, taking the GPS deformation rate as the true value, the RMSEs of the deformation rates calculated by the ordinary least squares method and the method proposed in this application are 8.6 mm / yr and 6.9 mm / yr respectively, and the external consistency accuracy of the deformation rate is improved by about 19.7%.

[0180] In the instance verification, taking the GPS deformation time series as the true value, the average RMSEs of the deformation time series calculated by the ordinary least squares method and the method proposed in this application are 7.8 mm and 4.4 mm respectively, and the external consistency accuracy of the deformation time series is improved by about 43.6%.

[0181] In summary, the InSAR turbulent atmosphere correction method provided by the embodiments of this application has the following advantages:

[0182] 1) Expand the random model from one-dimensional time to three-dimensional space-time, making up for the deficiency that the one-dimensional time random model of the traditional method does not consider the spatial correlation of the turbulent atmosphere;

[0183] 2) Introduce a window into the construction of the three-dimensional space-time random model, solving the problems of the overly large random model and memory overflow when considering the spatial correlation of the turbulent atmosphere;

[0184] 3) The three-dimensional space-time random model effectively describes the spatio-temporal characteristics of the turbulent atmosphere, improves the accuracy of InSAR surface deformation monitoring, and solves the problems that the traditional method depends on empirical parameters (such as the size of the filtering window) and is difficult to effectively suppress the turbulent atmosphere.

[0185] The above is the preferred implementation manner of this application. It should be noted that for those of ordinary skill in the art in this technical field, without departing from the principle described in this application, several improvements and refinements can still be made, and these improvements and refinements should also be regarded as the protection scope of this application.

Claims

1. An InSAR turbulent atmosphere correction method based on a three-dimensional random model, characterized in that, Including: Obtain M interferograms corresponding to N SAR images of the study area, and divide each interferogram into Y windows using a window of a preset size. Take the phase information within the y-th window of all interferograms as the y-th observation window: ; For each observation window, constructing a spatio-temporal three-dimensional stochastic model within the observation window according to the M interferograms; the spatio-temporal three-dimensional stochastic model is used to describe the spatio-temporal characteristics of the turbulent atmosphere and decorrelation noise within the observation window; For each observation window, constructing a functional model for all pixels within the observation window; the functional model is used to describe the relationship between the phase information within the observation window and the deformation parameters to be solved; For each observation window, based on the functional model and the spatio-temporal three-dimensional stochastic model corresponding to the observation window, using the weighted least squares method to solve for the deformation parameters, and obtaining the numerical values of the deformation parameters under the observation window; Stitching the numerical values of the deformation parameters under all observation windows, and calculating the deformation phase of the study area based on the stitching result.

2. The InSAR turbulent atmosphere correction method according to claim 1, characterized in that The constructing of the spatio-temporal three-dimensional stochastic model within the observation window according to the M interferograms includes: According to the M interferograms, constructing the spatio-temporal covariance matrix of the turbulent atmosphere within the observation window for each interferogram; Based on the constructed spatio-temporal covariance matrix of the turbulent atmosphere, determining the spatio-temporal covariance matrix of the turbulent atmosphere within the observation window for each SAR image scene; Based on the spatio-temporal covariance matrices of the turbulent atmosphere within the observation window for the N SAR image scenes, determining the spatio-temporal three-dimensional covariance matrix of the turbulent atmosphere within the observation window for the M interferograms; Based on the spatio-temporal three-dimensional covariance matrix of the turbulent atmosphere, determining the spatio-temporal three-dimensional stochastic model within the observation window.

3. The InSAR turbulent atmosphere correction method according to claim 2, wherein The spatial covariance matrix of the turbulent atmosphere within the observation window for the first interferogram is: ； Among them, , the covariance matrix of the turbulent atmosphere space The elements in , Indicates the distance between two pixels at positions and position in the observation window, , , the size of the observation window is , Indicates the number of rows of the observation window, Indicates the number of columns of the observation window; represents the baseline value of the th interferogram, represents the range value of the th interferogram.

4. The InSAR turbulence atmospheric correction method according to claim 3, wherein Before the step of constructing the spatio-temporal covariance matrix of the turbulent atmosphere within the observation window for each interferogram according to the M interferograms, the InSAR turbulent atmosphere correction method further includes: Using a spherical model to perform theoretical variogram fitting on each interferogram, and obtaining the sill value and range value corresponding to each interferogram.

5. The InSAR turbulent atmosphere correction method according to claim 3, wherein The determining of the spatio-temporal covariance matrix of the turbulent atmosphere within the observation window for each SAR image scene based on the constructed spatio-temporal covariance matrix of the turbulent atmosphere includes: Obtained by calculating using the following formula, the spatial covariance matrix of the turbulent atmosphere of the scene SAR image within the observation window is: ； Among them, represents the coefficient matrix. The coefficient matrix B is an M×N coefficient matrix, which is determined by the relationship between the interferogram and the SAR image. For each row of it, the position of the master image is 1, the position of the slave image is also 1, and the remaining positions are 0. .

6. The InSAR turbulent atmosphere correction method according to claim 5, characterized in that The determining of the spatio-temporal three-dimensional covariance matrix of the turbulent atmosphere within the observation window for the M interferograms based on the spatio-temporal covariance matrices of the turbulent atmosphere within the observation window for the N SAR image scenes includes: The three-dimensional spatio-temporal covariance matrix of the turbulent atmosphere within the observation window for M interference patterns is calculated through the following formula :[[]]END]] Among them, represents the coefficient matrix, , represents the coefficient matrix connecting the interferogram and the SAR image. The coefficient matrix F is an M×N coefficient matrix, which is determined by the relationship between the interferogram and the SAR image. For each row of it, the position of the master image is -1, the position of the slave image is 1, and the rest of the positions are 0. represents the Kronecker product operation, represents an identity matrix of size, , represents the turbulent atmospheric spatial covariance matrix of all SAR images within the observation window; 。 7. The InSAR turbulent atmosphere correction method according to claim 6, wherein The determining of the spatio-temporal three-dimensional stochastic model within the observation window based on the spatio-temporal three-dimensional covariance matrix of the turbulent atmosphere includes: Calculate the three-dimensional spatio-temporal random model using the following formula : Among them, represents the decorrelation noise variance matrix of all pixels in the observation window, represents the -th interferogram's decorrelation noise variance matrix of pixels in the observation window, The element in the observation window is the variance of the decorrelation noise of the pixel , , represents the number of pixels in the observation window, = , , represents the coherence coefficient of the pixel point , and L represents the number of looks.

8. The InSAR turbulent atmosphere correction method according to claim 7, characterized in that, The functional model is: Among them, represents the radar wavelength, represents the phase vectors of M interferograms within the observation window for the pixels, represents the coefficient matrix, represents the deformation parameters to be solved within the observation window, represents the residual phase, represents a coefficient matrix of M×4, represents the Kronecker product operation, represents an identity matrix of the size, , , represents the number of pixels within the observation window, represents the phase information corresponding to all interferograms for pixel p, represents the residual phase corresponding to all interferograms for pixel p, represents the cumulative time of the th scene SAR image relative to the first scene SAR image, represents the incident angle, represents the th scene SAR image and the th scene SAR image for the vertical baseline length of the generated interferogram, , , .

9. The InSAR turbulent atmosphere correction method according to claim 1, characterized in that, The number of deformation parameters is multiple, and the stitching of the numerical values of the deformation parameters under all observation windows includes: For each deformation parameter, calculating the mean value of the numerical values of the deformation parameter under all observation windows, and using the calculated mean value as the final value of the deformation parameter; Taking the final values of all deformation parameters as the stitching result.

Images