InSAR turbulent atmosphere correction method based on three-dimensional random model

Abstract

The invention is suitable for the technical field of geodetic survey, and provides an InSAR turbulent atmosphere correction method based on a three-dimensional random model, and the method comprises the steps: obtaining M interferograms corresponding to N-scene SAR images of a research region, dividing each interferogram into Y windows through a window with a preset size, taking the phase information in the yth window of all interferograms as the yth observation window; for each observation window, constructing a space-time three-dimensional random model in the observation window according to the M interferograms, constructing a function model of all pixels in the observation window, and solving by using a weighted least square method to obtain a deformation parameter corresponding to the observation window; the deformation parameters corresponding to all the observation windows are spliced, and the deformation phase of the research area is obtained through calculation based on the splicing result. The influence of turbulent atmosphere can be effectively inhibited.

Description

Technical Field

[0001] The present application belongs to the field of geodetic technology, and in particular to an InSAR turbulent atmosphere correction method based on a three-dimensional random model. Background Art

[0002] As a new generation of space geodetic technology, Interferometric Synthetic Aperture Radar (InSAR) has shown unique value in the field of surface deformation monitoring due to its advantages such as all-weather observation capability and high spatial resolution. This technology has been successfully applied to research fields such as earthquakes, volcanoes and landslides, and has accumulated rich scientific results. However, InSAR technology is affected by tropospheric delay, which not only significantly reduces the accuracy of surface deformation monitoring, but also leads to greater uncertainty in the deformation field interpretation process, becoming a bottleneck restricting the widespread application of this technology in high-precision surface deformation monitoring.

[0003] Tropospheric delay can be divided into vertical layer delay and turbulent atmosphere delay according to its physical properties. Vertical layer delay shows strong correlation with terrain in space, and is generally corrected by modeling it as a function of terrain. The spatiotemporal variation of turbulent atmosphere delay is relatively complex, and it is difficult to establish an exact functional model for it. Generally, three methods are used to suppress turbulent atmosphere delay: spatiotemporal filtering, interference pattern superposition, and random model. The spatiotemporal filtering method designs corresponding filters to achieve effective separation of turbulent atmosphere delay and deformation signal based on the significant difference in spatiotemporal characteristics. However, the core parameter of this method, the filter window, depends on the subjective determination of user experience. The interference pattern superposition method assumes that the turbulent atmosphere is randomly distributed in time. By differentiating and superimposing interference patterns with the same time interval, the turbulent atmosphere of the SAR image is obtained and the atmospheric influence is suppressed. However, this method is only applicable to linear deformation. The random model method suppresses the atmospheric influence by establishing the variance-covariance matrix of the turbulent atmosphere. However, the existing random model only considers the temporal correlation of the turbulent atmosphere in the interference pattern, and establishes the variance-covariance matrix of the time dimension of the turbulent atmosphere pixel by pixel, without considering the spatial correlation of the turbulent atmosphere.

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

[0005] The embodiment of the present application provides an InSAR turbulent atmosphere correction method based on a three-dimensional random model, which can solve the problem of difficulty in effectively suppressing the influence of the turbulent atmosphere.

[0006] The present application embodiment provides an InSAR turbulent atmosphere correction method based on a three-dimensional random model, comprising: Obtain M interferograms corresponding to N SAR images of the study area, and divide each interferogram into Y windows using a window of preset size, and use the phase information in the yth window of all interferograms as the yth observation window; ; For each observation window, a three-dimensional space-time random model is constructed according to M interference patterns. The three-dimensional space-time random model is used to describe the space-time characteristics of the turbulent atmosphere and incoherent noise in the observation window. For each observation window, a function model of all pixels in the observation window is constructed; the function model is used to describe the relationship between the phase information in the observation window and the deformation parameters to be solved; For each observation window, the deformation parameters are solved by weighted least square method based on the function model corresponding to the observation window and the three-dimensional random model of space and time, and the value of the deformation parameters under the observation window is obtained; The values ​​of the deformation parameters in all observation windows are spliced, and the deformation phase of the study area is calculated based on the splicing results.

[0007] Optionally, a spatiotemporal three-dimensional random model within the observation window is constructed according to the M interference patterns, including: According to the M interference patterns, the spatial covariance matrix of the turbulent atmosphere in the observation window is constructed for each interference pattern; Based on the constructed turbulent atmosphere spatial covariance matrix, determine the turbulent atmosphere spatial covariance matrix of each SAR image within the observation window; Based on the spatial covariance matrix of the turbulent atmosphere of N SAR images in the observation window, the spatial and temporal three-dimensional covariance matrix of the turbulent atmosphere of M interferograms in the observation window is determined; Based on the three-dimensional space-time covariance matrix of the turbulent atmosphere, the three-dimensional space-time random model within the observation window is determined.

[0008] Optional, The spatial covariance matrix of the turbulent atmosphere in the observation window for: ; in, , the spatial covariance matrix of the turbulent atmosphere Elements in , Indicates the position within the observation window and location The distance between two pixels on , , the size of the observation window is , Indicates the number of rows in the observation window, Indicates the number of columns in the observation window; Indicates The sill value of the interferogram, Indicates The range value of the interference pattern.

[0009] Optionally, before the step of constructing a turbulent atmosphere spatial covariance matrix of each interferogram within the observation window according to the M interferograms, the InSAR turbulent atmosphere correction method further includes: The spherical model is used to fit the theoretical variation function of each interference pattern to obtain the base value and range value corresponding to each interference pattern.

[0010] Optionally, based on the constructed turbulent atmosphere spatial covariance matrix, the turbulent atmosphere spatial covariance matrix of each SAR image in the observation window is determined, including: The following formula is used to calculate the The spatial covariance matrix of the turbulent atmosphere of the scene SAR image in the observation window : ; in, Represents the coefficient matrix. The coefficient matrix B is an M×N coefficient matrix determined by the relationship between the interferogram and the SAR image. For each row, the position of the main image is 1, the position of the slave image is also 1, and the rest of the positions are 0. .

[0011] Optionally, based on the spatial covariance matrix of the turbulent atmosphere of N SAR images in the observation window, the temporal and spatial three-dimensional covariance matrix of the turbulent atmosphere of M interference patterns in the observation window is determined, including: The three-dimensional covariance matrix of the turbulent atmosphere in space and time within the observation window of M interference patterns is calculated by the following formula: : in, represents the coefficient matrix, , It represents the coefficient matrix connecting the interferogram and the SAR image. The coefficient matrix F is an M×N coefficient matrix determined by the relationship between the interferogram and the SAR image. For each row, the position of the main image is -1, the position of the slave image is 1, and the rest of the positions are 0. represents the Kronecker product operation, Indicates a The identity matrix of size, , represents the spatial covariance matrix of the turbulent atmosphere of all SAR images within the observation window; .

[0012] Optionally, based on the turbulent atmosphere three-dimensional space-time covariance matrix, a three-dimensional space-time random model within the observation window is determined, including: The three-dimensional random model of space and time is calculated by the following formula : in, represents the incoherence noise variance matrix of all interferogram pixels in the observation window, Indicates The incoherence noise variance matrix of the pixels in the observation window of the interference pattern is Elements in is the pixel in the observation window The variance of the decoherent noise, , represents the number of pixels in the observation window, = , , Represents pixel The coherence coefficient of , L represents the number of multi-looks.

[0013] Optionally, the function model is: 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, , , .

[0014] Optionally, the number of deformation parameters is multiple, and the values ​​of the deformation parameters in all observation windows are spliced, including: For each deformation parameter, calculate the mean value of the deformation parameter in all observation windows, and use the calculated mean value as the final value of the deformation parameter; The final values ​​of all deformation parameters are taken as the stitching result.

[0015] The above solution of the present application has the following beneficial effects: In the embodiments of the present application, since the constructed spatiotemporal three-dimensional random model can express the spatiotemporal characteristics of the turbulent atmosphere, when solving the deformation parameters for deformation phase calculation based on the spatiotemporal three-dimensional random model and the function model of pixels, the influence of the turbulent atmosphere can be effectively suppressed, thereby effectively improving the accuracy of InSAR surface deformation monitoring.

[0016] In addition, since the construction of the three-dimensional space-time random model and the solution of the deformation parameters are carried out window by window, the accuracy of the three-dimensional space-time random model and deformation parameters is improved while greatly reducing the difficulty of calculation, which facilitates the effective suppression of the influence of the turbulent atmosphere.

[0017] Other beneficial effects of the present application will be described in detail in the subsequent specific implementation section. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0019] Figure 1 A flowchart of an InSAR turbulent atmosphere correction method based on a three-dimensional random model provided in one embodiment of the present application; Figure 2 This is a comparison diagram of deformation time series calculated by the ordinary least squares method and the method of the present application at the GPS station SACY in an example of the present application; Figure 3 This is a comparison diagram of deformation time series calculated by the ordinary least squares method and the method of the present application at the GPS station AZU1 in an example of the present application; Figure 4 A comparison diagram of deformation time series calculated by the ordinary least squares method and the method of the present application at the GPS station BLSA in an example of the present application; Figure 5 This is a comparison diagram of the deformation time series calculated by the ordinary least squares method and the method of the present application at the GPS station CAS4 in an example of the present application. DETAILED DESCRIPTION

[0020] In the following description, specific details such as specific system structures, technologies, etc. are provided for the purpose of illustration rather than limitation, so as to provide a thorough understanding of the embodiments of the present application. However, it should be clear to those skilled in the art that the present application may 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 prevent unnecessary details from obstructing the description of the present application.

[0021] It should be understood that when used in the present specification and the appended claims, the term "comprising" indicates the presence of 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 combinations thereof.

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

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

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

[0025] References to "one embodiment" or "some embodiments" etc. described in the specification of this application mean that one or more embodiments of the present application include specific features, structures or characteristics described in conjunction with the embodiment. Therefore, the statements "in one embodiment", "in some embodiments", "in some other embodiments", "in some 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 "including", "comprising", "having" and their variations all mean "including but not limited to", unless otherwise specifically emphasized in other ways.

[0026] In view of the current difficulty in effectively suppressing the influence of turbulent atmosphere, the embodiment of the present application provides an InSAR turbulent atmosphere correction method based on a three-dimensional random model, which determines the deformation phase of the study area by constructing a spatiotemporal three-dimensional random model window by window and solving the deformation parameters. Since the constructed spatiotemporal three-dimensional random model can express the spatiotemporal characteristics of the turbulent atmosphere, when solving the deformation parameters used for deformation phase calculation based on the spatiotemporal three-dimensional random model and the function model of pixels, the influence of the turbulent atmosphere can be effectively suppressed, thereby effectively improving the accuracy of InSAR surface deformation monitoring.

[0027] In addition, since the construction of the three-dimensional space-time random model and the solution of the deformation parameters are carried out window by window, the accuracy of the three-dimensional space-time random model and deformation parameters is improved while greatly reducing the difficulty of calculation, which facilitates the effective suppression of the influence of the turbulent atmosphere.

[0028] The InSAR turbulent atmosphere correction method based on a three-dimensional random model provided by the present application is exemplarily described below in conjunction with specific embodiments.

[0029] like Figure 1 As shown, the InSAR turbulent atmosphere correction method based on a three-dimensional random model provided in an embodiment of the present application includes the following steps: Step 11, 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, and use the phase information in the yth window of all interferograms as the yth observation window.

[0030] Above ; The above-mentioned study area is an area where surface deformation monitoring is required, and synthetic aperture radar (SAR) images can be collected by SAR sensors carried by satellites. After collecting N SAR images of the study area, the SAR images can be processed by the commonly used interferogram acquisition method to obtain an interferogram. That is, after selecting a main image, the remaining SAR images are registered, and then the time-space baseline threshold is set to perform differential processing to obtain M interferograms, and the flat ground phase and terrain phase are removed according to the external digital elevation model (DEM) data, and then the phase is unwrapped, and finally M unwrapped interferograms are obtained. It can be understood that the interferogram in the above step 11 is an unwrapped interferogram.

[0031] The size of the above window is , Indicates the number of rows in the window. Indicates the number of columns of the window. After dividing each interference pattern, the windows at the same position in all interference patterns can be used as an observation window, that is, the phase information in the yth window of all interference patterns is used as the yth observation window. It can be understood that in order to smoothly transition between windows, a certain overlap rate can be set between windows. Generally, the window size can be set to , the window moving step can be set to 24.

[0032] Step 12: for each observation window, construct a spatiotemporal three-dimensional random model within the observation window according to the M interference patterns.

[0033] The above-mentioned three-dimensional space-time random model is used to describe the space-time characteristics of the turbulent atmosphere and incoherent noise within the observation window.

[0034] Step 13, for each observation window, construct a function model of all pixels in the observation window.

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

[0036] Step 14, for each observation window, respectively, based on the function model corresponding to the observation window and the three-dimensional space-time random model, the deformation parameters are solved by using the weighted least squares method to obtain the value of the deformation parameters under the observation window.

[0037] That is, by using the weighted least squares method to solve the deformation parameters based on the function model corresponding to the observation window and the three-dimensional random model of space and time, the values ​​of the deformation parameters to be solved in step 13 are obtained, and the values ​​of the deformation parameters obtained by solving the deformation parameters are used as the values ​​of the deformation parameters under the observation window.

[0038] Step 15, stitching the values ​​of the deformation parameters in all observation windows, and calculating the deformation phase of the study area based on the stitching results.

[0039] In some embodiments of the present application, there are multiple deformation parameters. When stitching, the mean of the numerical values ​​of the deformation parameters under all observation windows can be calculated for each deformation parameter respectively, and the calculated mean value can be used as the final value of the deformation parameter; then the final values ​​of all deformation parameters are used as the stitching results.

[0040] Exemplarily, for the average deformation rate, the average of the values ​​of the average deformation rate in all observation windows may be used as the final value of the average deformation rate.

[0041] It should be noted that after solving the final values ​​of all deformation parameters, the corresponding temporal low-frequency deformation phase can be obtained, and then the residual phase is subjected to temporal low-pass filtering. The filtered phase is added to the temporal low-frequency deformation phase to obtain the deformation phase of the study area, which is convenient for deformation monitoring.

[0042] It is worth mentioning that since the constructed spatiotemporal three-dimensional random model can express the spatiotemporal characteristics of the turbulent atmosphere, when solving the deformation parameters used for deformation phase calculation based on the spatiotemporal three-dimensional random model and the pixel function model, it can effectively suppress the influence of the turbulent atmosphere, thereby effectively improving the accuracy of InSAR surface deformation monitoring.

[0043] In addition, since the construction of the three-dimensional space-time random model and the solution of the deformation parameters are carried out window by window, the accuracy of the three-dimensional space-time random model and deformation parameters is improved while greatly reducing the difficulty of calculation, which facilitates the effective suppression of the influence of the turbulent atmosphere.

[0044] The specific implementation method of constructing the spatiotemporal three-dimensional random model in the observation window according to the M interference patterns in step 12 is exemplarily described below in conjunction with a specific embodiment.

[0045] Specifically, the specific implementation method of constructing the spatiotemporal three-dimensional random model in the observation window according to the M interference patterns in the above step 12 includes the following steps: Step 12.1, based on the M interference patterns, construct the turbulent atmosphere spatial covariance matrix of each interference pattern within the observation window.

[0046] The above-mentioned turbulent atmosphere spatial covariance matrix can also be called the turbulent atmosphere spatial variance-covariance matrix.

[0047] In some embodiments of the present application, The spatial covariance matrix of the turbulent atmosphere in the observation window for: ; in, , the spatial covariance matrix of the turbulent atmosphere The elements on the main diagonal represent the variance of a pixel in the observation window, and the elements on the non-main diagonal represent the covariance between two pixels in the observation window. The spatial covariance matrix of the turbulent atmosphere is Elements in , Indicates the position within the observation window and location The distance between two pixels on , , the size of the observation window is , Indicates the number of rows in the observation window, Indicates the number of columns in the observation window.

[0048] in, Indicates The sill value of the interferogram, Indicates The range value of the interference pattern.

[0049] It can be understood that, before the step of constructing the turbulent atmosphere spatial covariance matrix of each interference pattern within the observation window based on the M interference patterns, the above-mentioned InSAR turbulent atmosphere correction method also includes the following step of determining the base value and the range value: using a spherical model to fit the theoretical variation function of each interference pattern to obtain the base value and the range value corresponding to each interference pattern.

[0050] In some embodiments of the present application, Sill value of the interferogram and The range of the interference pattern The process of obtaining is as follows: InSAR phase components are composed of: (1) is the interference phase after unwrapping, , , , and They are deformation phase, terrain residual phase, vertical stratified atmosphere phase, turbulent atmosphere phase and noise. Assuming vertical stratified atmosphere The phase-elevation model has been well removed. For the interferogram of the short time-space baseline, the deformation phase and terrain residual phase The influence of is small, and the noise is smaller than the turbulent atmosphere, so the main component of the phase is the turbulent atmosphere. The turbulent atmosphere has significant spatial correlation. As a core tool in geostatistics, the variogram is often used to analyze the spatial correlation of data. , the variation function Defined as: (2) It is the distance The variogram value at The distance is The number of data point pairs. The phase of the interference pattern is used as spatial data Substitute into formula (2) and set different distances Obtain the corresponding variance function value , It is the distance The number of Since the number of samples is always limited, the experimental variogram (i.e., variogram value) obtained above is not continuous in space. In order to obtain a spatially continuous variogram, it is necessary to fit the experimental variogram using a theoretical variogram. In this application, the following spherical model is used for fitting: (3) is the maximum variance function value, called the sill value, is the maximum correlation distance, called the range value. Fit the theoretical variation function for each interference pattern (i.e. And the corresponding variance function value After fitting) obtain the corresponding sill value and range value , M is the number of interference patterns.

[0051] Step 12.2, based on the constructed turbulent atmosphere spatial covariance matrix, determine the turbulent atmosphere spatial covariance matrix of each SAR image within the observation window.

[0052] The above-mentioned turbulent atmosphere spatial covariance matrix can also be called the turbulent atmosphere spatial variance-covariance matrix.

[0053] In some embodiments of the present application, the first The spatial covariance matrix of the turbulent atmosphere of the scene SAR image in the observation window : ; in, , Represents the coefficient matrix. The coefficient matrix B is an M×N coefficient matrix determined by the relationship between the interferogram and the SAR image. For each row, the position of the main image is 1, the position of the slave image is also 1, and the rest of the positions are 0.

[0054] in, The process of obtaining is as follows: Assume The interferogram is composed of SAR images j and k (i.e. Scene SAR image and The difference of the scene SAR image) is obtained, and each pixel in the observation window has the following propagation law according to the variance-covariance: It is The position of the interference pattern in the observation window and location The covariance between two pixels on (when b=d, c=e, represents variance). It is The variance-covariance corresponding to the scene SAR image, It is is the variance-covariance corresponding to the scene SAR image, and N is the number of SAR images.

[0055] Next, Rewritten into the following matrix form: ; is the interference pattern variance-covariance vector corresponding to a pixel, is the variance-covariance vector of the SAR image corresponding to a pixel. The variance-covariance of the SAR image corresponding to the pixel is obtained by least squares solution: Step 12.3, based on the spatial covariance matrix of the turbulent atmosphere of the N SAR images in the observation window, determine the spatiotemporal three-dimensional covariance matrix of the turbulent atmosphere of the M interference patterns in the observation window.

[0056] In some embodiments of the present application, the three-dimensional spatiotemporal covariance matrix of the turbulent atmosphere of the M interference patterns within the observation window can be calculated by the following formula: : in, represents the coefficient matrix, , It represents the coefficient matrix connecting the interferogram and the SAR image. The coefficient matrix F is an M×N coefficient matrix determined by the relationship between the interferogram and the SAR image. For each row, the position of the main image is -1, the position of the slave image is 1, and the rest of the positions are 0. represents the Kronecker product operation, Indicates a The identity matrix of size, , represents the spatial covariance matrix of the turbulent atmosphere of all SAR images within the observation window; .

[0057] Step 12.4, based on the three-dimensional space-time covariance matrix of the turbulent atmosphere, determine the three-dimensional space-time random model within the observation window.

[0058] In some embodiments of the present application, the spatiotemporal three-dimensional random model can be calculated by the following formula: : in, represents the incoherence noise variance matrix of all interferogram pixels in the observation window, Indicates The incoherence noise variance matrix of the pixels in the observation window of the interference pattern is Elements in is the pixel in the observation window The variance of the decoherent noise, , represents the number of pixels in the observation window, = , , Represents pixel The coherence coefficient of , L represents the number of multi-looks.

[0059] In some embodiments of the present application, the function model constructed in step 13 is: Among them, G is an M×4 coefficient matrix, The elements of the row are The parameters of the interference pattern are composed of 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 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, , Indicates that pixel p is The phase information corresponding to the interference pattern, 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, is 0, represents the angle of incidence, It represents the slant distance from the satellite (i.e. the satellite collecting SAR images) to the surface of the earth. Indicates Scene SAR image and the The vertical baseline length of the interferogram generated by the SAR image, , , .

[0060] The following is an exemplary description of the construction process of the function model: The small baseline set interferometric synthetic aperture radar (SBAS-InSAR) method often models the temporal low-frequency deformation and terrain residual phase. For the interferogram composed of the jth and kth SAR images, the temporal low-frequency deformation and terrain residual of pixel p in the observation window are as follows:

[0061] is the modeled temporal low-frequency deformation and terrain residual phase, is the radar wavelength, , and are the average deformation rate, average acceleration and average rate of change of acceleration, respectively. is the vertical baseline corresponding to the interference pattern, is the slant distance from the satellite to the Earth's surface, is the angle of incidence, is the terrain residual. and are the accumulated time of the jth and kth SAR images relative to the first SAR image. Considering that there are M interferograms, they can be written in matrix form as: ; is the phase of pixel p in all interference patterns, is an M×4 coefficient matrix, is the temporal low-frequency deformation parameter and terrain residual to be obtained for pixel p, is the residual phase, which is mainly composed of unmodeled deformation and turbulent atmospheric phase. The above only establishes a function model for one pixel. Since the method of this application is solved window by window, it is also necessary to establish a function model for all pixels in the window: It should be noted that for the new generation of satellites (such as Sentinel-1), the spatial baseline is controlled to be very short (usually less than 200 m), and the influence of terrain residuals is greatly weakened. Therefore, the influence of terrain residuals on Sentinel-1 satellite data can often be ignored, and the above function model can be further simplified.

[0062] 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 spatiotemporal three-dimensional random model in step 14 is: By using the weighted least squares method to solve the function model corresponding to the observation window and the three-dimensional random model of space and time, we can get: 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 terrain residual within the observation window.

[0063] The temporal low-frequency deformation phase can be obtained according to the solved temporal low-frequency deformation parameters. Since the residual phase contains a small amount of unmodeled deformation phase, it can be extracted by temporal low-pass filtering. At this time, the deformation phase in the residual phase accounts for a low proportion, so a large window filter can be used to extract the deformation phase. The final deformation phase within an observation window can be obtained by adding the temporal low-frequency deformation phase to the non-model deformation phase extracted by filtering from the residual phase.

[0064] The InSAR turbulent atmosphere correction method of the present application is exemplarily described below with reference to specific examples.

[0065] In this example, experiments were conducted using ESA's Sentinel-1 satellite data in an area with densely distributed Global Positioning System (GPS) stations, which can be used to evaluate InSAR deformation results. The Sentinel data covers 30 SAR images from January 10, 2018 to December 24, 2018. A total of 109 interferograms were generated by setting a spatial baseline threshold of 150 m and a temporal baseline threshold of 60 days. In order to suppress noise, the interferograms were multi-looked in 20×5 in range and azimuth directions with a spatial resolution of about 80 m, and the minimum cost flow method was used for phase unwrapping.

[0066] 1) Construct a three-dimensional random model of space and time within the observation window In order to reduce the influence of deformation phase and terrain residual phase, the short-time and space baseline interferogram generated above is used to obtain the experimental variogram. Since the method proposed in this application proposal is mainly aimed at turbulent atmosphere, in order to avoid the influence of vertical stratified atmosphere, the area with elevation higher than 250 m is masked. The step size is set to 160 m (corresponding to 2 pixels), and a series of experimental variogram values ​​are solved. The theoretical variogram of the spherical model is used to fit the experimental variogram to obtain the corresponding base value. and range value The spatial window size was set to 25 × 25 pixels, and a spatiotemporal three-dimensional random model was generated.

[0067] 2) Solving within the observation window Since the orbit control of the Sentinel satellite is good, the spatial baseline of all interferograms is less than 150 m. Therefore, the influence of terrain residuals can be ignored compared with deformation and turbulent atmosphere. Therefore, only the temporal low-frequency deformation phase is modeled as a function, and the temporal low-frequency deformation parameters within an observation window are solved using the weighted least squares method. Then, the moving step is set to 24 pixels, and the observation window is moved for solution. The overlapping area between the observation windows is averaged, and then the temporal low-frequency deformation phase is solved according to the temporal low-frequency deformation parameters. Finally, the non-model deformation phase is extracted from the residual phase using a temporal low-pass filter. The temporal low-pass filter window is set to 180 days. The temporal low-frequency deformation and the non-model deformation are added to obtain the final deformation phase.

[0068] The deformation rate of InSAR line of sight (LOS) is solved using the ordinary least squares method and the method of this application, and the deformation field accuracy of the two methods is evaluated using a series of indicators.

[0069] In order to evaluate the uncertainty of the deformation field of the two methods, the following indicators are used for evaluation: ; It is the epoch The corresponding displacement, is the predicted linear displacement, is the number of SAR images, is the cumulative time of the uth SAR image relative to the first SAR image, represents the uncertainty in the time series displacement due to temporal random residual noise, is the mean of the cumulative time.

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

[0071] After evaluating the intrinsic accuracy of the deformation rate, the external accuracy of the deformation rate is evaluated using external GPS data. First, the GPS three-dimensional deformation is projected to the LOS deformation using the imaging geometry of the SAR satellite. Then, the InSAR pixels within 200 m near each GPS station are searched, and the average deformation of these pixels is taken as the InSAR deformation corresponding to the station. The comparison results of the GPS deformation rate, the deformation rate calculated by the ordinary least squares method, and the deformation rate calculated by the method of this application are shown in Table 1. At most GPS stations, the deformation rate of the method of this application is closer to the deformation rate of GPS. The RMSE of the deformation rate of the ordinary least squares method and the method of this application are 8.01 mm / yr and 6.59 mm / yr, respectively, which indicates that the external accuracy of the deformation rate has been improved by about 17.73%.

[0072] Table 1 Deformation rate comparison results

[0073] It is also possible to verify that the deformation time series obtained by the method of the present 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 GPS deformation time series, and the deformation time series at the GPS station is plotted as follows: Figures 2 to 5 As shown in the figure, 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 and has higher consistency. Figures 2 to 5 The horizontal axis represents time, and the vertical axis represents displacement. Figures 2 to 5 The deformation time series at GPS stations SACY, AZU1, BLSA, and CAS4 are shown respectively.

[0074] Finally, the root mean square error of the deformation time series of each GPS station is counted. The root mean square error comparison results of the deformation time series solved by the ordinary least squares method and the deformation time series solved by the method of the present application are shown in Table 2. The RMSE mean values ​​of the deformation time series of the ordinary least squares method and the method of the present application are 7.83 mm and 4.26 mm, respectively, which indicates that the external symbol accuracy of the deformation time series is improved by about 45.59%.

[0075] Table 2 Comparison results of the root mean square error of deformation time series

[0076] From the above experimental data, it can be seen that the average uncertainties of the deformation rate calculated by the ordinary least squares method and the method proposed in this application proposal are 1.7 mm / yr and 0.8 mm / yr, respectively, and the intrinsic accuracy of the deformation rate is improved by about 52.9%.

[0077] In the example verification, taking the GPS deformation rate as the true value, the RMSE of the deformation rate solved by the ordinary least squares method and the method proposed in this application proposal are 8.6 mm / yr and 6.9 mm / yr respectively, and the external accuracy of the deformation rate is improved by about 19.7%.

[0078] In the example verification, taking the GPS deformation time series as the true value, the average RMSE of the deformation time series solved by the ordinary least squares method and the method proposed in this application proposal are 7.8 mm and 4.4 mm respectively, and the external symbol accuracy of the deformation time series is improved by about 43.6%.

[0079] In summary, the InSAR turbulent atmosphere correction method provided in the embodiment of the present application has the following advantages: 1) Extend the random model from one-dimensional time to three-dimensional space and time, making up for the deficiency of the traditional one-dimensional random model in time that does not consider the spatial correlation of turbulent atmosphere; 2) Introducing windows into the construction of three-dimensional random models in space and time solves the problem of random models being too large and memory overflowing when taking into account the spatial correlation of the turbulent atmosphere; 3) The spatiotemporal three-dimensional random model effectively describes the spatiotemporal characteristics of the turbulent atmosphere, improves the accuracy of InSAR surface deformation monitoring, and solves the problem that traditional methods rely on empirical parameters (such as filter window size) and are difficult to effectively suppress the turbulent atmosphere.

[0080] The above is a preferred embodiment of the present application. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles described in the present application. These improvements and modifications should also be regarded as the scope of protection of the present application.

Claims

1. A turbulent atmosphere correction method for InSAR based on a three-dimensional random model, characterized in that: include: Obtain M interferograms corresponding to N SAR images of the study area, and divide each interferogram into Y windows using a window of preset size, and use the phase information in the yth window of all interferograms as the yth observation window: ; For each observation window, a spatiotemporal three-dimensional random model is constructed according to the M interference patterns in the observation window; the spatiotemporal three-dimensional random model is used to describe the spatiotemporal characteristics of the turbulent atmosphere and the incoherent noise in the observation window; For each observation window, a function model of all pixels in the observation window is constructed; the function model is used to describe the relationship between the phase information in the observation window and the deformation parameter to be solved; For each observation window, respectively, the deformation parameter is solved by using a weighted least square method based on a function model corresponding to the observation window and a three-dimensional space-time random model to obtain a value of the deformation parameter under the observation window; The values ​​of the deformation parameters in all observation windows are spliced, and the deformation phase of the study area is calculated based on the splicing results.

2. The InSAR turbulent atmosphere correction method according to claim 1, characterized in that: The step of constructing the spatiotemporal three-dimensional random model within the observation window according to the M interference patterns comprises: According to the M interference patterns, constructing a turbulent atmosphere spatial covariance matrix of each interference pattern within the observation window; Based on the constructed turbulent atmosphere spatial covariance matrix, determining the turbulent atmosphere spatial covariance matrix of each SAR image within the observation window; Determine the spatiotemporal three-dimensional covariance matrix of the turbulent atmosphere of the M interference patterns in the observation window based on the spatial covariance matrix of the turbulent atmosphere of the N SAR images in the observation window; Based on the turbulent atmosphere spatiotemporal three-dimensional covariance matrix, a spatiotemporal three-dimensional random model within the observation window is determined.

3. The InSAR turbulent atmosphere correction method according to claim 2, characterized in that: No. The spatial covariance matrix of the turbulent atmosphere in the observation window is for: ； in, , the spatial covariance matrix of the turbulent atmosphere Elements in , Indicates the position within the observation window and location The distance between two pixels on , , the size of the observation window is , Indicates the number of rows in the observation window, Indicates the number of columns in the observation window; Indicates The sill value of the interferogram, Indicates The range value of the interference pattern.

4. The InSAR turbulent atmosphere correction method according to claim 3, characterized in that: Before the step of constructing a turbulent atmosphere spatial covariance matrix of each interference pattern within the observation window according to the M interference patterns, the InSAR turbulent atmosphere correction method further includes: The spherical model is used to fit the theoretical variation function of each interference pattern to obtain the base value and range value corresponding to each interference pattern.

5. The InSAR turbulent atmosphere correction method according to claim 3, characterized in that: Determining the turbulent atmosphere spatial covariance matrix of each SAR image within the observation window based on the constructed turbulent atmosphere spatial covariance matrix includes: The following formula is used to calculate the The spatial covariance matrix of the turbulent atmosphere of the SAR image in the observation window : ； in, Represents the coefficient matrix. The coefficient matrix B is an M×N coefficient matrix determined by the relationship between the interferogram and the SAR image. For each row, the position of the main image is 1, the position of the slave image is also 1, and the rest of the positions are 0. .

6. The InSAR turbulent atmosphere correction method according to claim 5, characterized in that: The determining of the turbulent atmosphere spatiotemporal three-dimensional covariance matrix of the M interference patterns in the observation window based on the turbulent atmosphere spatial covariance matrix of the N SAR images in the observation window comprises: The three-dimensional spatiotemporal covariance matrix of the turbulent atmosphere in the observation window of M interference patterns is calculated by the following formula: : in, represents the coefficient matrix, , It represents the coefficient matrix connecting the interferogram and the SAR image. The coefficient matrix F is an M×N coefficient matrix determined by the relationship between the interferogram and the SAR image. For each row, the position of the main image is -1, the position of the slave image is 1, and the rest of the positions are 0. represents the Kronecker product operation, Indicates a The identity matrix of size, , represents the spatial covariance matrix of the turbulent atmosphere of all SAR images within the observation window; 。 7. The InSAR turbulent atmosphere correction method according to claim 6, characterized in that: The step of determining the spatiotemporal three-dimensional random model within the observation window based on the spatiotemporal three-dimensional covariance matrix of the turbulent atmosphere comprises: The three-dimensional random model of space and time is calculated by the following formula : in, represents the incoherence noise variance matrix of all pixels in the observation window of the interference pattern, Indicates The incoherence noise variance matrix of the pixels in the observation window of the interferogram is Elements in is the pixel in the observation window The variance of the decoherent noise, , represents the number of pixels in the observation window, = , , Represents pixel The coherence coefficient of , L represents the number of multi-looks.

8. The InSAR turbulent atmosphere correction method according to claim 7, characterized in that: The function model is: 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, , , .

9. The InSAR turbulent atmosphere correction method according to claim 1, characterized in that: There are multiple deformation parameters, and the values ​​of the deformation parameters in all observation windows are spliced, including: For each deformation parameter, calculate the mean of the values ​​of the deformation parameter in all observation windows, and use the calculated mean as the final value of the deformation parameter; The final values ​​of all deformation parameters are taken as the stitching result.

Images