MT-InSAR atmospheric correction method based on independent principal component analysis and non-stationarity detection

Abstract

The present invention discloses an MT-InSAR atmospheric correction method based on independent principal component analysis and non-stationary detection. It includes: using M differential interferograms formed by time-series N-scene SAR images, after selecting high coherence points, correcting the orbit error phase and the terrain-related tropospheric delay phase, and combining the least squares estimation to obtain the differential interferometric phase accumulated by N scenes in the time series. Then, the study area is gridded based on the position information of the high coherence points, and the differential interferometric phase of each grid is processed by the spatial independent component analysis (sICA) algorithm and the non-stationary detection algorithm, and the atmospheric error signal and the deformation source signal are identified and reconstructed, so as to calculate the time-series deformation results of the grid. The deformation results of high coherence points in the overlapping areas of adjacent grids are used to achieve the splicing of the deformation results of the overall study area, which effectively improves the accuracy of time-series InSAR deformation monitoring and avoids dependence on high-temporal and spatial resolution meteorological data or subjective filtering methods.

Description

Technical Field

[0001] The present invention relates to the technical field of time series synthetic aperture radar interferometry (TS-InSAR), and in particular to an MT-InSAR atmospheric correction method without independent principal component analysis and non-stationarity detection. Background Art

[0002] Surface deformation refers to the phenomenon of instantaneous deformation or slow and continuous deformation in a local area caused by natural factors or human activities. Among them, instantaneous deformation mainly includes deformation caused by natural disasters such as volcanoes, earthquakes, landslides, and mudslides; slow and continuous deformation is mainly caused by large-scale artificial construction in cities, large-scale exploitation of groundwater, oil, natural gas and other resources, and ground subsidence induced by other human life and industrial activities. The surface deformation caused by these natural factors and human activities will directly threaten the property and life safety of the general public. Therefore, accurate and effective long-term monitoring of surface deformation will help to fully grasp the inducing factors and evolution mechanisms of geological disasters and deformation accidents of buildings and structures, evaluate the distribution and development dynamics of geological disaster risks, and then carry out timely and effective governance, which is an effective way to avoid various casualties and losses.

[0003] Traditional surface deformation monitoring is mainly carried out by means of leveling total stations and global satellite navigation systems. This type of monitoring method is simple to operate and has high accuracy. However, this method has low spatial resolution, high labor costs, slow data updates, and is easily restricted by the measurement environment. It is difficult to describe the real changes of the surface in a long-term, comprehensive, and high-resolution manner. Synthetic aperture radar interferometry (InSAR) technology is considered to be a space geodetic technology with great potential due to its advantages of all-weather, all-day, high accuracy, large area, and high spatial resolution. With the development of technology and the enrichment of SAR images, InSAR technology plays an increasingly important role in topographic mapping and geological disaster monitoring. However, the accuracy of deformation monitoring using traditional differential interferometry technology is greatly affected by factors such as atmospheric errors. Based on this, time-series InSAR technologies represented by permanent scatterer synthetic aperture radar interferometry technology and differential interferometry short baseline set time series analysis technology have been proposed to improve the accuracy of deformation monitoring.

[0004] It is worth noting that the above two types of time-series InSAR technologies both need to correct the atmospheric delay phase in the differential interferogram when estimating the surface deformation sequence. In response to this problem, many scholars at home and abroad have conducted a lot of research to suppress or correct the differential atmospheric delay phase in the interferogram. It can be mainly divided into two categories: correcting the atmospheric delay phase in the interferogram based on various external meteorological data or atmospheric water vapor delay products and numerical meteorological models, and correcting it based on the combination relationship of the SAR image itself.

[0005] Atmospheric correction methods that rely on interferograms of external meteorological data are largely completely dependent on the temporal and spatial resolution of externally acquired meteorological data or water vapor data, which limits the improvement of deformation accuracy of such methods. In addition, it is difficult to ensure that external meteorological data can be accurately obtained in all regions and at all times around the world. Based on this, a method that relies on interference to perform atmospheric correction on itself was proposed. The logical matching method was proposed in 1995. This method is only applicable to interferograms with public SAR images, so it cannot reduce the impact of atmospheric delay phase in non-public SAR images in the interferogram. The time-series interferogram overlay method achieves the theoretical reduction of the atmospheric delay phase in N interferograms to 0.01% by averaging N independent interferograms. Obviously, the correction accuracy of this method mainly depends on the number of time-series InSAR interferograms. The spatiotemporal filtering method is currently the most commonly used time-series InSAR atmospheric correction method. This method can better process the non-modeled deformation signal of the test area. However, for different experimental areas, the selection of filters and filter windows is still a difficult problem.

[0006] Based on the above analysis of the existing atmospheric delay phase correction algorithms in time-series InSAR surface deformation monitoring, it can be seen that the existing methods are limited by the spatial resolution of external meteorological data, meteorological conditions, the number of SAR images and the optimization selection algorithm of the filter, which restricts the application and promotion of the current InSAR atmospheric correction algorithm, thereby affecting the accuracy of time-series InSAR technology for surface deformation monitoring. Summary of the invention

[0007] In view of the defects of the prior art, the present invention provides an MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection.

[0008] This method introduces the interferogram information and spatial independent component analysis knowledge of synthetic aperture radar data sets, and uses the coherence map obtained by differential interferometry processing to complete the selection of high coherence points; then, the first SAR image is used as the reference image to obtain the differential interferometry phase accumulated in N scenes in time series; then, the position information and coordinate information of high coherence points in the interferogram are fully utilized to realize the interferogram grid division; for each grid, the spatial independent component analysis (sICA) algorithm is applied to obtain the mixing matrix A of the corresponding grid; then, the non-stationary detection algorithm is applied to each column of the mixing matrix A of the grid to identify the deformation source signal, and the deformation phase is reconstructed using the identified columns in the mixing matrix A, and the temporal deformation result of the grid is calculated according to the reconstructed phase component; finally, the deformation results of the pixels in the overlapping area of ​​adjacent grids are used to realize the splicing of the deformation of the study area.

[0009] In order to achieve the above invention object, the technical solution adopted by the present invention is as follows:

[0010] An MT-InSAR atmospheric correction method based on independent principal component analysis and non-stationarity detection comprises the following steps:

[0011] S1: For M differential interference images and M coherence images formed by differential interference of N SAR images covering the same study area, after selecting high coherence points using the coherence images, the first SAR image is used as the reference image to obtain the differential interference phase accumulated by N images in time sequence;

[0012] S2: Using the position information of the high coherence points in the interference pattern, the interference pattern is gridded;

[0013] S3: For each grid, the spatial independent component analysis (sICA) algorithm is applied to the N-scene accumulated differential interference phase to obtain the mixing matrix A of the corresponding grid;

[0014] S4: applying a non-stationary detection algorithm to each column of the grid's mixing matrix A to identify the atmospheric error signal and the deformation source signal, then using the identified columns in the mixing matrix A to reconstruct the deformation phase, and calculating the grid's temporal deformation result based on the reconstructed phase component;

[0015] S5: The deformation results of the study area are spliced ​​through the deformation results of the high coherence points in the overlapping areas of adjacent grids.

[0016] Furthermore, the S1 includes:

[0017] S11: Obtain N SAR images covering the same study area in time series, use differential interferometry to form M differential interferograms and M coherence maps, and after phase unwrapping, obtain the differential interferometry phases corresponding to the M differential interferograms. and the coherence matrix C = [C1, C2, ..., C i , …, CM ] T .in With C i They represent the differential interference phase and coherence value corresponding to the i-th differential interference pattern respectively;

[0018] S12: Calculate the maximum and minimum values ​​of the coherence matrix C, compare the temporal coherence value of any pixel in the differential interferogram with the average value of the maximum and minimum values ​​of the matrix C, and realize the adaptive selection of high coherence points in the interferogram. Based on the high coherence points, the orbit error phase and the terrain-related tropospheric delay phase correction are completed. The corrected differential interferometer phase is φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T ,φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern.

[0019] S13: Taking the first SAR image as the reference image, establish a function model between the differential interferometry phase of N SAR images and the differential interferometry phase of M interferograms, and obtain the cumulative differential interferometry phase corresponding to the imaging moment of N SAR images in time series by combining the least squares estimation algorithm. in It represents the cumulative differential interferometry phase corresponding to the imaging moment of the i-th SAR image.

[0020] Further, the S12 includes:

[0021] S121: Calculate the maximum and minimum values ​​of the coherence matrix C, and calculate the average value thereof using the following formula (1):

[0022]

[0023] Among them, c max With c min Respectively represent the maximum and minimum values ​​of the coherence matrix C, Indicates c max With c min The average value of

[0024] S122: For any pixel p in the interference graph, compare its corresponding coherence value in the M interference graphs with If the coherence values ​​corresponding to pixel p in the M interference patterns are all greater than Then pixel p is selected as a high coherence point, otherwise it is a non-high coherence point.

[0025] S123: For the highly coherent points selected above, a function model of the differential observation phase, the corresponding pixel and the terrain is established to estimate the orbit error phase and the terrain-related tropospheric delay phase in each interferogram, as shown in formula (2):

[0026] μ=a0+a1x+a2y+a3xy+a4x 2 +a5y 2 +a6h (2)

[0027] Among them, μ represents the original differential interference phase corresponding to the high coherence point in the differential interference diagram, x and y represent the coordinate position corresponding to the high coherence point, h represents the elevation value corresponding to this high coherence point, a0,…,a5 and a6 represent the unknown parameters corresponding to the orbit error phase and the terrain-related tropospheric delay phase, and the estimation of the unknown parameters can be obtained by least squares estimation. Then it is removed from the original differential interference phase φ, and the differential interference phase after removal is φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T . Where φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern.

[0028] Further, the S13 includes:

[0029] S131: Taking the first SAR image as the reference image, a function model between the differential interferometric phases of the N SAR images and the differential interferometric phases of the M interferograms is established, as shown in formula (3):

[0030]

[0031] Among them, φ o represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the M differential interference patterns, and φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T . Where φ o irepresents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern; represents the cumulative differential interferometric phase corresponding to the imaging moment of N SAR images in the time series, and in represents the cumulative differential interferometry phase corresponding to the imaging moment of the i-th SAR image;

[0032] The coefficient matrix B is shown in formula (4):

[0033]

[0034] S132: Using the least squares estimation algorithm to obtain the cumulative differential interferometric phase corresponding to the imaging time of N SAR images in the time series That is, as shown in formula (5):

[0035]

[0036] Among them, B T Represents the transpose of matrix B.

[0037] Further, the S2 includes:

[0038] Using the position information of high coherence points in the interference pattern, the study area is evenly divided into regular rectangular grids with a side length of 1 km, while ensuring that there are 10% common high coherence points between adjacent grids.

[0039] Further, the S3 includes:

[0040] The accumulated differential interference phase corresponding to the imaging time of N SAR images in the time sequence In each grid, the spatial independent component analysis (sICA) algorithm is used to obtain the mixing matrix A and the source signal matrix S. Phase The functional relationship with the mixing matrix A is:

[0041]

[0042] Among them, the phase It is a two-dimensional matrix of N×q, where q is the number of highly coherent points contained in each grid. The mixing matrix A is N×n s A two-dimensional matrix, N represents the number of time-series SAR images, n s Indicates the number of independent signal sources. S is composed of n s The source signals consist of n s ×q two-dimensional source matrix.

[0043] Further, the S4 includes:

[0044] S41: Apply the non-stationarity detection algorithm to each column of the mixing matrix A. 0 With the alternative hypothesis H 1 They are:

[0045]

[0046] Where t represents the imaging time of the SAR image, y(t) and y(t-1) represent the observation data of the time series at time t and time t-1, respectively, c, d, and e(t) represent the drift coefficient, trend term coefficient, and residual term of the regression model, respectively;

[0047] S42: According to the Perron-Phillips unit root test method, under the condition of significance level α, if we accept H 0 , if the column is considered non-stationary, delete the column corresponding to the mixing matrix A and the row corresponding to the source matrix S; otherwise, if the column is considered stable, keep the column corresponding to the mixing matrix A and the row corresponding to the source matrix S;

[0048] S43: The matrices after removing the non-stationary parts of the mixing matrix A and the source matrix S are respectively recorded as matrices And the matrix The deformation phase It is expressed as:

[0049]

[0050] Among them, the deformation phase It is a two-dimensional matrix of N×q, where q is the number of highly coherent points contained in each grid, and the mixing matrix yes The two-dimensional matrix, Represents the number of columns remaining after the non-stationary part of the mixing matrix A is deleted, and the source signal matrix yes A two-dimensional matrix.

[0051] Further, the S5 includes:

[0052] S51: obtaining common high coherence points in the overlapping area, and calculating the difference between the deformation results of high coherence points in the overlapping area of ​​adjacent grids;

[0053] S52: According to the positions of the highly coherent points of each grid, a spatial adjustment network is constructed, and the difference between the deformation results is combined to perform an overall correction of the grid blocking effect using the least squares network adjustment, so as to obtain a spatially continuous and smooth deformation result.

[0054] The present invention also discloses an MT-InSAR atmosphere correction system, which can be used to implement the above-mentioned MT-InSAR atmosphere correction method, specifically comprising:

[0055] Data preprocessing module: obtain N SAR images covering the same study area, generate M differential interference images and M coherence images, complete phase unwrapping, and obtain the differential interference phase φ and coherence matrix C.

[0056] High coherence point selection module: calculates the maximum and minimum values ​​of the coherence matrix C, compares the temporal coherence value of the pixel with the average value of C, selects high coherence points, and performs orbit error phase and terrain-related tropospheric delay phase correction based on the high coherence points.

[0057] Differential interferometry phase calculation module: Taking the first SAR image as a reference, a function model between the differential interferometry phase and the differential interferometry phase of the M interferograms is established, and the least squares estimation algorithm is used to obtain the cumulative differential interferometry phase corresponding to the imaging moment of the N SAR images in the time series.

[0058] Grid division module: Utilize the position information of high coherence points in the interference pattern to divide the interference pattern into grids.

[0059] Independent principal component analysis module: For each grid, the spatial independent component analysis (sICA) algorithm is applied to the N-scene accumulated differential interference phase to obtain the mixing matrix A of the corresponding grid.

[0060] Non-stationarity detection module: The non-stationarity detection algorithm is applied to each column of the grid's mixing matrix A to identify the atmospheric error signal and the deformation source signal, and then the deformation phase is reconstructed using the identified columns in the mixing matrix A. The temporal deformation result of the grid is calculated based on the reconstructed phase component.

[0061] Deformation result splicing module: The deformation results of the study area are spliced ​​through the deformation results of high coherence points in the overlapping areas of adjacent grids.

[0062] The present invention also discloses a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the above-mentioned MT-InSAR atmospheric correction method when executing the program.

[0063] The present invention also discloses a computer-readable storage medium, on which a computer program is stored. When the program is executed by a processor, the above-mentioned MT-InSAR atmosphere correction method is implemented.

[0064] Compared with the prior art, the advantages of the present invention are:

[0065] 1. It does not rely on any external data, but makes full use of the spatial independent component analysis algorithm and the interferogram's own temporal and spatial distribution characteristics of atmospheric errors for correction.

[0066] 2. Avoid the errors caused by unreasonable selection of filters and filter windows in the current spatiotemporal filtering algorithm.

[0067] 3. The method is simple and efficient, easy to implement and widely used.

[0068] 4. It will help promote the development and application of time-series InSAR interferometry technology in engineering applications.

[0069] 5. Improved the precision and accuracy of deformation monitoring by effectively suppressing the influence of turbulent atmospheric phase in the interference pattern.

[0070] 6. It can better handle deformation monitoring problems under complex terrain and environmental conditions, such as mountainous areas, cities and other areas.

[0071] 7. The intelligent selection of high coherence points in the interference pattern and the efficient stitching of deformation results are realized, making the results more comprehensive and reliable.

[0072] 8. By reducing the reliance on external meteorological data, the complexity of data acquisition and processing is reduced, while the autonomy and independence of processing are improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 The figure is a flow chart of the MT-InSAR atmospheric correction method according to an embodiment of the present invention.

[0074] Figure 2 This is a spatiotemporal baseline distribution diagram of the InSAR interferogram simulated in an embodiment of the present invention.

[0075] Figure 3 Result diagrams of simulation experiments of embodiments of the present invention; (a) is a simulated deformation velocity diagram, (b) and (c) are two simulated differential interference diagrams affected by atmospheric errors; (df) are deformation velocity diagrams of the method proposed in the present invention, the Kalman-based method, and the traditional SBAS estimation, respectively. DETAILED DESCRIPTION

[0076] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and examples.

[0077] Reference Figure 1 , an MT-InSAR atmospheric correction method based on independent principal component analysis and non-stationarity detection, including:

[0078] S1: For the M differential interference images and M coherence images formed by differential interference of N SAR images covering the same study area, after completing the selection of high coherence points, the first SAR image is used as the reference image to obtain the differential interference phase accumulated by N scenes in time sequence.

[0079] In this embodiment, the sub-steps of step S1 are as follows: S11: Obtain N SAR images covering the study area and on the same track, complete image registration and resampling, perform interference processing on the master and slave SAR images to generate M interference patterns, and then use external digital elevation model (DEM) data to simulate the terrain phase of the study area to obtain M differential interference patterns covering the study area. After performing phase unwrapping and coherence calculation on the M differential interference patterns, obtain the differential interference phase corresponding to the M differential interference patterns and the coherence matrix C = [C1, C2, ..., C i , …, C M ] T .in With C i They respectively represent the differential interference phase and coherence value corresponding to the i-th differential interference pattern.

[0080] S12: Calculate the maximum and minimum values ​​of the coherence matrix C, and calculate the average value thereof using the following formula (1):

[0081]

[0082] Among them, c max With c min Respectively represent the maximum and minimum values ​​of the coherence matrix C, Indicates c max With c min For any pixel p in the interference graph, compare its corresponding coherence value in the M interference graphs with the coherence value in formula (1) If the coherence values ​​corresponding to pixel p in the M interference patterns are all greater than c, then pixel p is selected as a high coherence point, otherwise it is a non-high coherence point.

[0083] S13: Using the highly coherent points selected above, a function model of the differential observation phase and the corresponding pixels and terrain is established to estimate the orbit error phase and the terrain-related tropospheric delay phase in each interferogram, as shown in formula (2).

[0084] μ=a0+a1x+a2y+a3xy+a4x 2 +a5y 2 +a6h (2)

[0085] Among them, μ represents the original differential interference phase corresponding to the high coherence point in the differential interference diagram, x and y represent the coordinate position corresponding to the high coherence point, h represents the elevation value corresponding to this high coherence point, a0,…,a5 and a6 represent the unknown parameters corresponding to the orbit error phase and the terrain-related tropospheric delay phase, and the estimation of the unknown parameters can be obtained by least squares estimation. Then it is removed from the original differential interference phase φ, and the differential interference phase after removal is φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T . Where Φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern.

[0086] S14: Taking the first SAR image as the reference image, a function model between the differential interferometric phase of the N SAR images and the differential interferometric phase of the M interferograms is established, as shown in formula (3):

[0087]

[0088] Among them, φ o represents the differential interference phase corresponding to the M-amplitude differential interference pattern, and φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T ,φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern; represents the cumulative differential interferometric phase corresponding to the imaging moment of N SAR images in the time series, and in represents the cumulative differential interferometry phase corresponding to the imaging time of the i-th SAR image; the coefficient matrix B is shown in formula (4):

[0089]

[0090] S15: Use the least squares estimation algorithm to obtain the cumulative differential interferometric phase corresponding to the imaging time of N SAR images in the time series That is, as shown in formula (5):

[0091]

[0092] S2: Using the position information and coordinate information of the high coherence points in the interference pattern, the interference pattern is divided into grids.

[0093] In this embodiment, the sub-steps of step S2 are as follows:

[0094] S21: Using the position information of high coherence points in the interference pattern, the study area is evenly divided into regular rectangular grids with a side length of 1 km, while ensuring that there are 10% common high coherence points between adjacent grids.

[0095] S3: For each grid, the spatial independent component analysis (sICA) algorithm is applied to the N-scene accumulated differential interference phase to obtain the mixing matrix A of the corresponding grid.

[0096] In some specific embodiments, step S3 may further include:

[0097] S31: accumulating differential interferometric phases corresponding to the imaging moments of N SAR images in the time sequence The spatial independent component analysis (sICA) algorithm is applied in each grid to obtain the mixing matrix A and the source signal matrix S. Phase The functional relationship with the mixing matrix A is:

[0098]

[0099] Among them, the phase It is a two-dimensional matrix of N×q, where q is the number of highly coherent points contained in each grid. The mixing matrix A is N×n s A two-dimensional matrix, N represents the number of time-series SAR images, n s Indicates the number of independent signal sources. S is composed of n s The source signals consist of n s ×q two-dimensional source matrix.

[0100] S4: Apply the non-stationary detection algorithm to each column of the grid's mixing matrix A to identify the atmospheric error signal and the deformation source signal, then use the identified columns in the mixing matrix A to reconstruct the deformation phase, and calculate the grid's temporal deformation result based on the reconstructed phase component.

[0101] In this embodiment, the sub-steps of step S4 are as follows: S41: Apply the non-stationarity detection algorithm to each column of the above mixing matrix A. 0 With the alternative hypothesis H 1 They are:

[0102]

[0103] Where t represents the time of the SAR image, y(t) and y(t-1) represent the observation data of the time series at time t and time t-1, respectively, c, d, and e(t) represent the drift coefficient, trend term coefficient, and residual term of the regression model, respectively;

[0104] S42: According to the Perron-Phillips unit root test method, under the condition of significance level α, if we accept H 0 , if the column is considered non-stationary, delete the column corresponding to the mixing matrix A and the row corresponding to the source matrix S; otherwise, if the column is considered stable, retain the column corresponding to the mixing matrix A and the row corresponding to the source matrix S;

[0105] S43: The matrices after removing the non-stationary parts of the mixing matrix A and the source matrix S are respectively recorded as matrices And the matrix The deformation phase It is expressed as:

[0106]

[0107] Among them, the deformation phase It is a two-dimensional matrix of N×q, where q is the number of highly coherent points contained in each grid, and the mixing matrix yes The two-dimensional matrix, Represents the number of columns remaining after the non-stationary part of the mixing matrix A is deleted, and the source signal matrix yes A two-dimensional matrix.

[0108] S5: The deformation results of the study area are spliced ​​through the deformation results of the high coherence points in the overlapping areas of adjacent grids.

[0109] In this embodiment, the sub-steps of step S5 are as follows: S51: obtaining common high coherence points in the overlapping area, and calculating the difference in deformation results of high coherence points in the overlapping areas of adjacent grids;

[0110] S52: According to the positions of the highly coherent points of each grid, a spatial adjustment network is constructed, and the difference between the deformation results is combined to perform an overall correction of the grid blocking effect using the least squares network adjustment, so as to obtain a spatially continuous and smooth deformation result.

[0111] Through the above S1-S5, the MT-InSAR high-precision surface deformation solution based on independent principal component analysis and non-stationarity detection is realized, which avoids the need to set different filters and filter windows for different research areas in the traditional time-series InSAR method, and avoids the dependence on external high-temporal and spatial resolution external meteorological data. The method proposed in the present invention is easy to implement and helps to promote and promote the development of the current time-series InSAR interferometry technology and its application in engineering.

[0112] Based on the simulation data, the specific implementation steps of the present invention are as follows:

[0113] According to the parameters of the Sentinel-1 satellite data set, 70 SAR images lasting two years were simulated. Using the temporal and spatial baseline threshold method, a total of 233 interferograms were generated. The corresponding temporal and spatial baseline distribution is shown in the figure. Figure 2 As shown in Figure 2. The simulated deformation speed is shown in Figure 2. Figure 3 (a) The atmospheric error in the interferogram is simulated by a numerical function, and the maximum delay phase range of the entire simulation process is 5 rad to 12 rad, as shown in Figure 3 As shown in (b) and (c), it is assumed that there is no phase component related to terrain and orbit information, and the 233 interferograms have been correctly phase-unwrapped.

[0114] In order to prove the superiority of the method proposed in this invention, Figure 3 (d), (e), and (f) respectively show the deformation velocity results estimated by the method proposed in this invention, the Kalman-based method, and the traditional small baseline subset (SBAS) method. Figure 3 (a) and 3(d), it can be seen that the deformation results estimated by the method proposed in the present invention have a high consistency with the simulated real deformation results in space. In contrast, the deformation results of the Kalman-based method 3(e) and the traditional SBAS method 3(f) show obvious atmospheric residual errors. Furthermore, the RMSE of the deformation results obtained by the method proposed in the present invention is 0.5mm, while the RMSE of the deformation results obtained by the traditional SBAS method is as high as 3.28mm. This once again illustrates the reliability of the deformation results estimated by the method of the present invention.

[0115] In yet another embodiment of the present invention, a MT-InSAR atmosphere correction system is provided. The system can be used to implement the above-mentioned MT-InSAR atmosphere correction method, specifically comprising:

[0116] Data preprocessing module: obtain N SAR images covering the same study area, generate M differential interference images and M coherence images, complete phase unwrapping, and obtain the differential interference phase φ and coherence matrix C.

[0117] High coherence point selection module: calculates the maximum and minimum values ​​of the coherence matrix C, compares the temporal coherence value of the pixel with the average value of C, selects high coherence points, and performs orbit error phase and terrain-related tropospheric delay phase correction based on the high coherence points.

[0118] Differential interferometry phase calculation module: Taking the first SAR image as a reference, a function model between the differential interferometry phase and the differential interferometry phase of the M interferograms is established, and the least squares estimation algorithm is used to obtain the cumulative differential interferometry phase corresponding to the imaging moment of the N SAR images in the time series.

[0119] Grid division module: Utilize the position information of high coherence points in the interference pattern to divide the interference pattern into grids.

[0120] Independent principal component analysis module: For each grid, the spatial independent component analysis (sICA) algorithm is applied to the N-scene accumulated differential interference phase to obtain the mixing matrix A of the corresponding grid.

[0121] Non-stationarity detection module: The non-stationarity detection algorithm is applied to each column of the grid's mixing matrix A to identify the atmospheric error signal and the deformation source signal, and then the deformation phase is reconstructed using the identified columns in the mixing matrix A. The temporal deformation result of the grid is calculated based on the reconstructed phase component.

[0122] Deformation result splicing module: The deformation results of the study area are spliced ​​through the deformation results of high coherence points in the overlapping areas of adjacent grids.

[0123] Result visualization module: used to display the processed deformation monitoring results, including the following:

[0124] Time series deformation result diagram: displays the deformation results at different time points or time periods, such as deformation rate, deformation amount, etc.

[0125] Spatial deformation result map: shows the deformation distribution of different locations in the study area.

[0126] Abnormal detection results: Mark possible abnormal deformation or atmospheric interference signals to help users analyze the accuracy and credibility of the results.

[0127] Data comparison and analysis: Compare and analyze the differences before and after processing to help users evaluate the effect of atmospheric correction and the accuracy of deformation monitoring.

[0128] In another embodiment of the present invention, a terminal device is provided, the terminal device includes a processor and a memory, the memory is used to store a computer program, the computer program includes program instructions, and the processor is used to execute the program instructions stored in the computer storage medium. The processor may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc., which are the computing core and control core of the terminal, which are suitable for implementing one or more instructions, and are specifically suitable for loading and executing one or more instructions to implement the corresponding method flow or corresponding functions; the processor described in the embodiment of the present invention can be used for the operation of the MT-InSAR atmospheric correction method.

[0129] In another embodiment of the present invention, a storage medium is provided, specifically a computer-readable storage medium (Memory), which is a memory device in a terminal device for storing programs and data. It is understandable that the computer-readable storage medium here can include both the built-in storage medium in the terminal device and the extended storage medium supported by the terminal device. The computer-readable storage medium provides a storage space, which stores the operating system of the terminal. In addition, one or more instructions suitable for being loaded and executed by the processor are also stored in the storage space, and these instructions can be one or more computer programs (including program codes). It should be noted that the computer-readable storage medium here can be a high-speed RAM memory or a non-volatile memory, such as at least one disk memory.

[0130] One or more instructions stored in a computer-readable storage medium may be loaded and executed by a processor to implement the corresponding steps of the MT-InSAR atmospheric correction method in the above embodiment; one or more instructions in the computer-readable storage medium may be loaded and executed by a processor.

[0131] It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0132] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0133] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0134] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0135] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the implementation methods of the present invention, and should be understood that the protection scope of the present invention is not limited to such special statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. A MT-InSAR atmospheric correction method based on independent principal component analysis and non-stationarity detection, characterized in that: The following steps are involved: S1: For M differential interference images and M coherence images formed by differential interference of N SAR images covering the same study area, after selecting high coherence points using the coherence images, the first SAR image is used as the reference image to obtain the differential interference phase accumulated by N images in time sequence; S2: Using the position information of the high coherence points in the interference pattern, the interference pattern is gridded; S3: For each grid, the spatial independent component analysis (sICA) algorithm is applied to the N-scene accumulated differential interference phase to obtain the mixing matrix A of the corresponding grid; S4: applying a non-stationary detection algorithm to each column of the grid's mixing matrix A to identify the atmospheric error signal and the deformation source signal, then using the identified columns in the mixing matrix A to reconstruct the deformation phase, and calculating the grid's temporal deformation result based on the reconstructed phase component; S4 includes the following sub-steps: S41: Apply the non-stationarity detection algorithm to each column of the mixing matrix A; where the original hypothesis H 0 With the alternative hypothesis H 1 They are: Where t represents the imaging time of the SAR image, y(t) and y(t-1) represent the observation data of the time series at time t and time t-1, respectively, c, d, and e(t) represent the drift coefficient, trend term coefficient, and residual term of the regression model, respectively; S42: According to the Perron-Phillips unit root test method, under the condition of significance level α, if we accept H 0 , if the column is considered non-stationary, delete the column corresponding to the mixing matrix A and the row corresponding to the source matrix S; otherwise, if the column is considered stable, retain the column corresponding to the mixing matrix A and the row corresponding to the source matrix S; S43: The matrices after removing the non-stationary parts of the mixing matrix A and the source matrix S are respectively recorded as matrices And the matrix The deformation phase It is expressed as: Among them, the deformation phase It is a two-dimensional matrix of N×q, where q is the number of highly coherent points contained in each grid, and the mixing matrix yes The two-dimensional matrix, Represents the number of columns remaining after the non-stationary part of the mixing matrix A is deleted, and the source signal matrix yes A two-dimensional matrix of S5: The deformation results of the study area are spliced ​​through the deformation results of the high coherence points in the overlapping areas of adjacent grids.

2. The MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection according to claim 1, characterized in that: The S1 includes: S11: Obtain N SAR images covering the same study area in time series, use differential interferometry to form M differential interferograms and M coherence maps, and after phase unwrapping, obtain the differential interferometry phases corresponding to the M differential interferograms. and the coherence matrix C = [C1, C2, ..., C i , …, C M ] T ;in With C i They represent the differential interference phase and coherence value corresponding to the i-th differential interference pattern respectively; S12: Calculate the maximum and minimum values ​​of the coherence matrix C, compare the temporal coherence value of any pixel in the differential interferogram with the average value of the maximum and minimum values ​​of the matrix C, and realize the adaptive selection of high coherence points of the interferogram; complete the orbit error phase and terrain-related tropospheric delay phase correction based on the high coherence points, and the differential interferometer phase after correction is φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T ,φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern; S13: Taking the first SAR image as the reference image, establish a function model between the differential interferometry phase of N SAR images and the differential interferometry phase of M interferograms, and obtain the cumulative differential interferometry phase corresponding to the imaging moment of N SAR images in time series by combining the least squares estimation algorithm. in It represents the cumulative differential interferometry phase corresponding to the imaging moment of the i-th SAR image.

3. The MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection according to claim 2, characterized in that: The S12 includes: S121: Calculate the maximum and minimum values ​​of the coherence matrix C, and calculate the average value thereof using the following formula (1): Among them, c max With c min Respectively represent the maximum and minimum values ​​of the coherence matrix C, Indicates c max With c min The average value of S122: For any pixel p in the above interference graph, compare its corresponding coherence value in the M interference graphs with If the coherence values ​​corresponding to pixel p in the M interference patterns are all greater than Then pixel p is selected as a high coherence point, otherwise it is a non-high coherence point; S123: for the selected high coherence points, a function model of the differential observation phase, the corresponding pixel and the terrain is established to estimate the orbit error phase and the terrain-related tropospheric delay phase in each interferogram, as shown in formula (2); μ=a0+a1x+a2y+a3xy+a4x 2 +a5y 2 +a6h (2) Among them, μ represents the original differential interference phase corresponding to the high coherence point in the differential interference diagram, x and y represent the coordinate positions corresponding to the high coherence point, h represents the elevation value corresponding to the high coherence point, a0,…,a5 and a6 represent the unknown parameters corresponding to the orbit error phase and the terrain-related tropospheric delay phase, respectively. The estimated value of the unknown parameters can be obtained by least squares estimation; then it is removed from the original differential interference phase φ, and the differential interference phase after removal is φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T ; where φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern.

4. The MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection according to claim 2, characterized in that: The S13 includes: S131: Taking the first SAR image as the reference image, a function model between the differential interferometric phases of the N SAR images and the differential interferometric phases of the M interferograms is established, as shown in formula (3): Among them, φ o represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the M differential interference patterns, and φ o =[φ o 1, φ o 2,…,φ o i ,…,φ o M ] T ; where φ o i represents the differential interference phase after removing the orbit error phase and the terrain-related tropospheric delay phase from the i-th differential interference pattern; represents the cumulative differential interferometric phase corresponding to the imaging moment of N SAR images in the time series, and in represents the cumulative differential interferometry phase corresponding to the imaging time of the i-th SAR image; the coefficient matrix B is shown in formula (4): S132: Using the least squares estimation algorithm to obtain the cumulative differential interferometric phase corresponding to the imaging time of N SAR images in the time series That is, as shown in formula (5): Among them, B T Represents the transpose of matrix B.

5. The MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection according to claim 1, characterized in that: The S2 includes: Using the position information of high coherence points in the interference pattern, the study area is evenly divided into regular rectangular grids with a side length of 1 km, while ensuring that there are 10% common high coherence points between adjacent grids.

6. The MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection according to claim 1, characterized in that: The S3 includes: The accumulated differential interference phase corresponding to the imaging time of N SAR images in the time sequence In each grid, the spatial independent component analysis (sICA) algorithm is used to obtain the mixing matrix A and the source signal matrix S; the phase The functional relationship with the mixing matrix A is: Among them, the phase It is a two-dimensional matrix of N×q, where q is the number of highly coherent points contained in each grid. The mixing matrix A is N×n s A two-dimensional matrix, N represents the number of time-series SAR images, n s Indicates the number of independent signal sources; S is composed of n s The source signals consist of n s ×q two-dimensional source matrix.

7. The MT-InSAR atmosphere correction method based on independent principal component analysis and non-stationarity detection according to claim 1, characterized in that: The S5 includes: S51: obtaining common high coherence points in the overlapping area, and calculating the difference between the deformation results of high coherence points in the overlapping area of ​​adjacent grids; S52: According to the positions of the highly coherent points of each grid, a spatial adjustment network is constructed, and the difference between the deformation results is combined to perform an overall correction of the grid blocking effect using the least squares network adjustment, so as to obtain a spatially continuous and smooth deformation result.

8. An MT-InSAR atmospheric correction system, characterized by: The system can be used to implement the MT-InSAR atmosphere correction method described in any one of claims 1 to 7, specifically comprising: Data preprocessing module: obtain N SAR images covering the same study area, generate M differential interferograms and M coherence maps, complete phase unwrapping, and obtain differential interferometry phase φ and coherence matrix C; High coherence point selection module: calculates the maximum and minimum values ​​of the coherence matrix C, compares the temporal coherence value of the pixel with the average value of C, selects high coherence points, and performs orbit error phase and terrain-related tropospheric delay phase correction based on the high coherence points; Differential interferometry phase calculation module: Taking the first SAR image as a reference, a function model between the differential interferometry phase and the differential interferometry phase of the M interferograms is established, and the cumulative differential interferometry phase corresponding to the imaging moment of the N SAR images in the time series is obtained using the least squares estimation algorithm; Grid division module: Grid division of the interference pattern is performed using the position information of the high coherence points in the interference pattern; Independent principal component analysis module: For each grid, the spatial independent component analysis (sICA) algorithm is applied to the N-scene accumulated differential interference phase to obtain the mixing matrix A of the corresponding grid; Non-stationarity detection module: Apply the non-stationarity detection algorithm to each column of the grid's mixing matrix A to identify the atmospheric error signal and the deformation source signal, then use the identified columns in the mixing matrix A to reconstruct the deformation phase, and calculate the grid's temporal deformation result based on the reconstructed phase component; Deformation result splicing module: The deformation results of the study area are spliced ​​through the deformation results of high coherence points in the overlapping areas of adjacent grids.

9. A computer device, characterized in that: The invention comprises a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor implements the MT-InSAR atmosphere correction method according to any one of claims 1 to 7 when executing the program.

Images