Ground surface deformation inversion method considering atmospheric effect

Abstract

The invention discloses a surface deformation inversion method considering an atmospheric effect, and relates to the technical field of remote sensing. The method comprises the following steps: screening out differential interferograms influenced by atmospheric pollution from a plurality of initial differential interferograms to obtain a plurality of differential interferograms; performing principal component analysis on the evaluation factors of the differential interferogram to obtain the weight of each evaluation factor; performing weighted summation on the values of all the evaluation factors and the weights of the evaluation factors to obtain a comprehensive score of the differential interferogram, performing screening according to the comprehensive score to obtain a target differential interferogram, and constructing an interferogram network; performing interference net inversion to obtain a phase value of each pixel in the interference pattern network, and sequentially performing earth tide phase correction, atmospheric delay phase correction and terrain phase residual correction on the phase values of the pixels to obtain corrected phase time sequence data; and obtaining earth surface deformation information of the research area based on the corrected phase time sequence data. According to the method, the precision of surface deformation inversion is improved.

Description

Technical Field

[0001] This invention relates to the field of remote sensing technology, and in particular to a method for inverting surface deformation that takes into account atmospheric effects. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR) technology has the advantages of being available in all weather conditions, covering a wide area, and providing high precision. Combined with space geodesy and active microwave remote sensing technologies, it is widely used in the field of earth sciences and has become an important tool for monitoring ground deformation caused by earthquakes, volcanoes, landslides, permafrost, man-made structures, groundwater extraction, underground oil and mining. It has shown great potential in areas such as surface deformation monitoring and geological disaster early warning.

[0003] In the practical application of InSAR technology, obtaining high-precision surface deformation fields is one of the core objectives, and the achievement of this objective highly depends on the selection of high-quality interferometric image pairs. Interferometric image pairs often contain various non-deformation phase errors. If low-quality interferometric image pairs containing such errors are directly used for analysis, it will seriously interfere with the accuracy of deformation signal extraction. Therefore, constructing a high-quality interferogram network and screening out effective, low-error interferograms from a massive number of interferometric image pairs has become a key step in improving the accuracy of InSAR surface deformation inversion.

[0004] To achieve the construction of high-quality interferometric networks, various interferometric network construction models have been proposed in existing technologies. Among them, Short Baseline Assembled Aperture Radar Interferometry (SBAS-InSAR) has become one of the mainstream technologies for temporal deformation monitoring due to its ability to effectively reduce the impact of temporal and spatial decorrelation. However, the overall quality of the interferometric networks generated by existing SBAS-InSAR interferometric pair construction methods is relatively low, which further leads to low inversion accuracy of surface deformation. Summary of the Invention

[0005] Therefore, it is necessary to provide a surface deformation inversion method that takes into account atmospheric effects to address the aforementioned technical problems. This method improves the accuracy of surface deformation inversion.

[0006] The present invention adopts the following technical solution: This invention provides a method for inverting surface deformation that takes into account atmospheric effects, comprising: Acquire temporal SAR images of the study area within a preset time period, and combine multiple SAR images in the temporal SAR images by interferometric image pairing to generate multiple initial differential interferograms; Multiple differential interferograms are obtained by filtering out differential interferograms affected by atmospheric pollution from the multiple initial differential interferograms generated; For each differential interferogram, principal component analysis is performed on the evaluation factors of the differential interferogram to obtain the weight of each evaluation factor; the evaluation factors include phase standard deviation, elevation phase correlation coefficient and spatial correlation vanishing distance; The values ​​of all evaluation factors are weighted and summed to obtain the comprehensive score of the differential interferogram. The differential interferograms are sorted in descending order of comprehensive score. A predetermined number of differential interferograms are selected as target differential interferograms, and an interferogram network is constructed based on the target differential interferograms. Interferometric inversion is performed based on the interferogram network to obtain the phase value of each pixel in the interferogram network. For each pixel, the phase value of the pixel is successively corrected by solid tide phase correction, atmospheric delay phase correction, and topographic phase residual correction to obtain the corrected phase time series data. Based on the corrected phase time series data, the surface deformation information of the study area is obtained.

[0007] Preferably, multiple SAR images in a time-series SAR imagery are combined by interferometric image pairs to generate multiple initial differential interferograms, specifically including: For each SAR image in the time-series SAR imagery, the SAR image is registered using the sequential network enhanced spectral diversity method; the number of connections during the registration process is set to a preset number of connections. An initial differential interferogram is generated using the registered SAR imagery.

[0008] Preferably, differential interferograms affected by atmospheric pollution are filtered out from the generated multiple initial differential interferograms to obtain multiple differential interferograms, specifically including: The initial differential interferogram is sequentially processed by multi-view processing, phase correction, phase filtering, and phase unwrapping; The short-range scale component of the initial differential interferogram is determined by the phase standard deviation of the initial differential interferogram after unwrapping. The topographic correlation components of the initial differential interferogram are determined using the elevation phase correlation coefficient; The wide-scale components of the initial differential interferogram are determined based on the spatial correlation vanishing distance of the optimal fitting model of the semi-variogram. Based on short-range scale components, terrain-related components, and wide-area scale components, the differential interferograms affected by atmospheric pollution are screened out from the multiple initial differential interferograms generated using the three-standard-deviation principle, resulting in multiple differential interferograms.

[0009] Preferably, the initial differential interferogram is subjected to multi-view processing, phase correction, phase filtering, and phase unwrapping in sequence, specifically including: Perform multi-view processing on the initial differential interferogram; Based on digital elevation model data, phase components caused by flat terrain effect and topographic relief are removed from the initial interferometric phase map after multi-view processing in order to perform phase correction; Phase filtering is applied to the phase correction result; The phase filtering result is unwrapped using the minimum cost flow algorithm to obtain the initial differential interferogram after unwrapping.

[0010] Preferably, based on the corrected phase time series data, the surface deformation information of the study area is obtained, specifically including: The terrain phase residuals in the corrected phase time series data are evaluated to obtain the evaluation results. The evaluation results include the terrain phase residual values ​​of each phase cell, the statistical distribution characteristics of the terrain phase residuals, the spatial dispersion of the terrain phase residuals, and the results of the terrain phase residual exceeding the limit. Based on the evaluation results, the phase time series data were optimized by trend surface fitting correction, spatial interpolation correction and residual compensation correction. The optimized phase time series data were then subjected to secondary residual verification to remove abnormal phase data that failed the verification, and the target phase time series data of each phase element was obtained. Based on the target phase time series data of all phase elements, the surface deformation information of the study area is inverted.

[0011] Preferably, principal component analysis is performed on the evaluation factors of the difference interferogram to obtain the weight of each evaluation factor, specifically including: Calculate the Pearson correlation coefficient between any two evaluation factors in the differential interferogram, and construct a correlation coefficient matrix based on the calculated Pearson correlation coefficient; Based on the correlation coefficient matrix, a characteristic equation is constructed, and the characteristic equation is solved to obtain the characteristic value corresponding to each evaluation index. Based on multiple eigenvalues, the cumulative variance contribution rate of each evaluation factor is calculated, and multiple principal components with a cumulative variance contribution rate greater than or equal to a preset percentage are selected. Calculate the linear combination coefficient of each evaluation index in different principal components; Calculate the cumulative variance contribution rate of each principal component; The initial weight of each evaluation factor is calculated based on the cumulative variance contribution rate of the principal components and the linear combination coefficient of each evaluation factor in different principal components. The initial weights are normalized to obtain the weight of each evaluation factor.

[0012] Preferably, the correlation coefficient matrix is ​​as follows: ; in, R The correlation coefficient matrix, For the first differential interferogram The evaluation factor and the first Pearson correlation coefficients among the evaluation factors p The number of evaluation factors; The characteristic equation is: ; in, For eigenvalues, E for p An identity matrix of order 1; No. k The formula for calculating the cumulative variance contribution rate of each evaluation factor is as follows: ; in, For the first k The cumulative variance contribution rate of each evaluation factor For the first k The characteristic values ​​of each evaluation factor; No. j The evaluation factor is at the _ i The formula for calculating the linear combination coefficients in each principal component is as follows: ; in, For the first j The evaluation factor is at the ... i The linear combination coefficients of the principal components Represented as the first i Principal component pairs of evaluation factors The load factor, , , m The number of principal components, p The number of evaluation factors; No. i The formula for calculating the cumulative variance contribution rate of each principal component is: ; in, For the first i Cumulative variance contribution rate of each principal component; No. j The formula for calculating the initial weights of each evaluation factor is: ; in, For the first j The initial weights of each evaluation factor, ; No. j The formula for calculating the weight of each evaluation factor is: ; in, For the first j The weights of each evaluation factor For the first k The initial weights of each evaluation factor. This invention provides a surface deformation inversion device that takes atmospheric effects into account, comprising:

[0013] The acquisition module is used to acquire time-series SAR images of the study area within a preset time period, and to combine multiple SAR images in the time-series SAR images by interferometric image pairing to generate multiple initial differential interferograms. The filtering module is used to filter out differential interferograms affected by atmospheric pollution from the multiple generated initial differential interferograms, resulting in multiple differential interferograms; The principal component analysis module is used to perform principal component analysis on the evaluation factors of each difference interferogram to obtain the weight of each evaluation factor. The evaluation factors include phase standard deviation, elevation phase correlation coefficient, and spatial correlation vanishing distance. The module is used to perform a weighted summation of the values ​​of all evaluation factors and their weights to obtain the comprehensive score of the differential interferogram. The differential interferograms are sorted in descending order of comprehensive score. A preset number of differential interferograms are selected as target differential interferograms, and an interferogram network is constructed based on the target differential interferograms. The inversion module is used to perform interferometric network inversion based on the interferogram network to obtain the phase value of each pixel in the interferogram network. For each pixel, the phase value of the pixel is sequentially corrected by solid tide phase, atmospheric delay phase, and topographic phase residual to obtain the corrected phase time series data. The determination module is used to obtain surface deformation information of the study area based on the corrected phase time series data.

[0014] The present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described method for inverting surface deformation taking into account atmospheric effects.

[0015] The present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the above-mentioned method for inverting surface deformation taking into account atmospheric effects.

[0016] The above-mentioned at least one technical solution adopted in this invention can achieve the following beneficial effects: The study acquires temporal SAR images of the study area within a preset time period. Multiple SAR images are combined using interferometric pairs to generate multiple initial differential interferograms (DIAs). These initial DIAAs provide complete surface phase information for the study area within the preset time span, laying a solid foundation for subsequent deformation inversion. DIAAs affected by atmospheric pollution are filtered out from the generated initial DIAAs, resulting in multiple DIAAs. Removing these low-quality DIAAs allows for the early elimination of core interference sources, effectively reducing the impact of invalid atmospheric phase on subsequent inversion calculations. For each DIAAs, principal component analysis is performed on the evaluation factors to obtain the weight of each factor. The values ​​of all evaluation factors are weighted and summed to obtain a comprehensive score for the DIAAs. The DIAAs are then sorted in descending order of comprehensive score, and a preset number of DIAAs are selected as target DIAAs. An interferometric network is constructed based on the target differential interferogram (DII), transforming multi-dimensional evaluation factors into a single comprehensive quality score. This enables quantitative ranking of the quality of each DII. The interferometric network, built upon the selected DIIIs, exhibits higher phase consistency and data reliability compared to networking based on all initial interferograms, effectively reducing phase calculation errors during interferometric network inversion. Interferometric network inversion is performed based on this network to obtain the phase value of each pixel. For each pixel, the phase value is sequentially corrected for solid tide phase, atmospheric delay phase, and topographic phase residuals, yielding corrected phase time-series data. A step-by-step correction method is then applied, targeting the three core non-deformation phase errors in the interferometric network phase. Based on the corrected phase time-series data, surface deformation information for the study area is obtained. This method improves the accuracy of surface deformation inversion. Attached Figure Description

[0017] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this invention, illustrate exemplary embodiments of the invention and are used to explain the invention, but do not constitute an undue limitation of the invention. In the drawings:

[0018] Figure 1 A schematic diagram of the process for a surface deformation inversion method that takes atmospheric effects into account, provided by the present invention; Figure 2 The flowchart of the InSAR interferometric pair construction method provided by this invention is shown below; Figure 3 A schematic diagram of a surface deformation inversion device that takes atmospheric effects into account, provided by the present invention; Figure 4A schematic diagram of a computer device for implementing a surface deformation inversion method that takes atmospheric effects into account, as provided by the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0020] Devices such as desktop computers, servers, and laptops are capable of executing the present invention. For ease of explanation, the following description will focus on servers as the executing entity.

[0021] The technical solutions provided by the various embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0022] Figure 1 This is a schematic diagram of a surface deformation inversion method that takes atmospheric effects into account in this invention, specifically including the following steps: S101: Acquire temporal SAR images of the study area within a preset time period, and combine multiple SAR images in the temporal SAR images by interferometric image pairing to generate multiple initial differential interferograms.

[0023] In an exemplary embodiment, multiple SAR images in a time-series SAR imagery are combined by interferometric image pairs to generate multiple initial differential interferograms. Specifically, this includes: for each SAR image in the time-series SAR imagery, registering the SAR image using a sequential network enhancement spectrum diversity method; setting the number of connections during the registration process to a preset number of connections; and generating an initial differential interferogram using the registered SAR imagery.

[0024] Specifically, SAR data preprocessing was performed using the open-source software Interferometric Synthetic Aperture Radar Scientific Computing Environment (ISCE).

[0025] Step 1: Use the network-based enhanced spectral diversity (NESD) method to refine and register SAR images. Considering data processing quality and efficiency, the number of connections is set to 5.

[0026] Step 2: Using the registered images, generate differential interferograms through interferometry; perform interferometry on the registered images to generate original differential interferograms including information such as surface deformation, atmospheric delay, and topographic phase.

[0027] Step 3: Perform multi-view processing on the results obtained from differential interferometry; the purpose of multi-view processing is to reduce noise.

[0028] Step 4: Based on data such as the Digital Elevation Model (DEM), remove the phase components caused by flat terrain effect and topographic relief from the interferometric phase; use DEM data to separate the topographic phase while retaining the deformation and atmospheric phase components.

[0029] Step 5: Considering data processing quality and efficiency, set the number of connections and set the filter strength to 0.2 for phase filtering.

[0030] Step 6: Use the Minimum Cost Flow (MCF) algorithm for phase unwrapping, etc., to solve the phase transition problem and obtain a continuous phase field.

[0031] S102: From the generated multiple initial differential interferograms, the differential interferograms affected by atmospheric pollution are screened out to obtain multiple differential interferograms.

[0032] In an exemplary embodiment, differential interferograms affected by air pollution are filtered out from multiple generated initial differential interferograms to obtain multiple differential interferograms. Specifically, this includes: performing multi-view processing, phase correction, phase filtering, and phase unwrapping sequentially on the initial differential interferograms; determining the short-range scale component of the initial differential interferograms using the phase standard deviation of the unwrapped initial differential interferograms; determining the terrain-related component of the initial differential interferograms using the elevation phase correlation coefficient; determining the wide-area scale component of the initial differential interferograms based on the spatial correlation vanishing distance of the optimal fitting model of the semi-variogram; and filtering out differential interferograms affected by air pollution from the multiple generated initial differential interferograms based on the short-range scale component, terrain-related component, and wide-area scale component using the three-times-standard-deviation principle to obtain multiple differential interferograms.

[0033] In an exemplary embodiment, the initial differential interferogram is sequentially subjected to multi-view processing, phase correction, phase filtering, and phase unwrapping, specifically including: performing multi-view processing on the initial differential interferogram; based on digital elevation model data, removing phase components caused by flat terrain effect and topographic undulation from the multi-view processed initial interferometric phase map to perform phase correction; performing phase filtering on the phase correction result; and using the minimum cost flow algorithm to perform phase unwrapping on the phase filtering result to obtain the unwrapped initial differential interferogram.

[0034] Specifically, tropospheric atmospheric delay inversion and component evaluation are performed. The phase standard deviation of the unwrapped interferogram is used to evaluate the short-range scale components of the tropospheric atmosphere. The elevation phase correlation coefficient is used to evaluate the topographically related tropospheric atmospheric components, and the spatial correlation vanishing distance (effective distance) of the optimal fitting model of the semi-variogram is used to evaluate the wide-scale components. Using the three-standard-deviation principle, interferograms significantly affected by air pollution are screened out from each evaluation factor dimension, achieving preliminary filtering of the initial difference interferogram.

[0035] Specifically, 1. For each initial differential interferogram generated by the solution, extract the atmospheric phase values ​​of the corresponding short-range scale component, topographic-related component, and wide-area scale component to obtain the atmospheric phase datasets of the three scale components for each initial differential interferogram; 2. Perform statistical analysis on the atmospheric phase datasets of the three scale components to calculate the mean and standard deviation of the atmospheric phase for the short-range scale component, the topographic-related component, and the wide-area scale component; 3. Based on the principle of three times the standard deviation, set the atmospheric phase threshold interval for each scale component: the atmospheric phase threshold interval for each scale component is the mean of the atmospheric phase of each component - 3 × the standard deviation of the atmospheric phase of each component to... 4. Threshold judgment is performed on the atmospheric phase values ​​of the three scale components of each initial differential interferogram. If the atmospheric phase value of any scale component exceeds its corresponding atmospheric phase threshold interval, the initial differential interferogram is judged to be a differential interferogram affected by air pollution. If the atmospheric phase values ​​of the three scale components all fall within their respective atmospheric phase threshold intervals, the initial differential interferogram is judged to be a valid differential interferogram without significant air pollution. 5. Based on the above judgment results, differential interferograms judged to be affected by air pollution are removed from all initial differential interferograms, and all valid differential interferograms judged to be without significant air pollution are retained, finally obtaining multiple differential interferograms after screening.

[0036] S103: For each differential interferogram, perform principal component analysis on the evaluation factors of the differential interferogram to obtain the weight of each evaluation factor; the evaluation factors include phase standard deviation, elevation phase correlation coefficient and spatial correlation vanishing distance.

[0037] In an exemplary embodiment, principal component analysis is performed on the evaluation factors of the difference interferogram to obtain the weight of each evaluation factor. Specifically, this includes: calculating the Pearson correlation coefficient between any two evaluation factors in the difference interferogram, and constructing a correlation coefficient matrix based on the calculated Pearson correlation coefficient. A characteristic equation is constructed based on the correlation coefficient matrix, and the characteristic equation is solved to obtain the eigenvalues ​​corresponding to each evaluation index. Based on multiple eigenvalues, the cumulative variance contribution rate of each evaluation factor is calculated, and multiple principal components with cumulative variance contribution rates greater than or equal to a preset percentage are selected. The linear combination coefficients of each evaluation index in different principal components are calculated. The cumulative variance contribution rate of each principal component is calculated. Based on the cumulative variance contribution rate of the principal components and the linear combination coefficients of each evaluation factor in different principal components, the initial weight of each evaluation factor is calculated. The initial weights are normalized to obtain the weight of each evaluation factor.

[0038] The correlation coefficient matrix is ​​given by formula (1): (1); in, R The correlation coefficient matrix, For the first differential interferogram The evaluation factor and the first Pearson correlation coefficients among the evaluation factors p This represents the number of evaluation factors.

[0039] The characteristic equation is formula (2): (2); in, For eigenvalues, E for p An identity matrix of order 1.

[0040] No. k The formula for calculating the cumulative variance contribution rate of each evaluation factor is formula (3): (3); in, For the first k The cumulative variance contribution rate of each evaluation factor For the first k The eigenvalues ​​of each evaluation factor.

[0041] No. j The evaluation factor is at the ... i The formula for calculating the linear combination coefficients in each principal component is formula (4): (4); in, For the first j The evaluation factor is at the ... i The linear combination coefficients of the principal components Represented as the first i Principal component pairs of evaluation factors The load factor, , , m The number of principal components, p This represents the number of evaluation factors.

[0042] No. i The formula for calculating the cumulative variance contribution rate of each principal component is formula (5): (5); in, For the first i Cumulative variance contribution rate of each principal component; No. j The formula for calculating the initial weight of each evaluation factor is formula (6): (6); in, For the first j The initial weights of each evaluation factor, ; No. j The formula for calculating the weight of each evaluation factor is formula (7): (7); in, For the first j The weights of each evaluation factor For the first k The initial weights of each evaluation factor.

[0043] S104: Weighted sum of the values ​​of all evaluation factors and their weights to obtain the comprehensive score of the differential interferogram. Sort the differential interferograms in descending order of comprehensive score, select a preset number of differential interferograms as target differential interferograms, and construct an interferogram network based on the target differential interferograms.

[0044] For the initial differential interferograms after preliminary screening, the evaluation factors are first standardized and principal component analysis is performed. The weights of each factor are calculated based on the cumulative variance contribution rate, and a comprehensive scoring framework with fused weights is established. Then, based on the statistical distribution of the score values, the interferogram with the smallest atmospheric delay is selected, and a high-quality interferogram network is constructed to complete the ground deformation inversion.

[0045] Based on the constructed interferogram network, interferometric network inversion is performed. When constructing the first interferogram network, low spatial coherence is masked and high-quality pixel target points with temporal coherence greater than 0.6 are selected for interferometric network inversion.

[0046] Specifically, the first 148 interferograms are selected from smallest to largest to construct an interferogram network.

[0047] S105: Perform interferometric network inversion based on the interferogram network to obtain the phase value of each pixel in the interferogram network. For each pixel, perform solid tide phase correction, atmospheric delay phase correction, and topographic phase residual correction on the phase value of the pixel in sequence to obtain the corrected phase time series data.

[0048] Specifically, considering the influence of Earth's solid tides on the surface phase, solid tide phase correction is performed on the pixel phase based on the theoretical model and relevant parameters of solid tides. The atmospheric delayed phase is corrected using the fifth-generation meteorological reanalysis dataset method (ERA5 method) published by the European Centre for Medium-Range Weather Forecasts (ECMWF). Based on the previous removal of flat and topographic phases, the topographic phase residuals are further corrected. After a series of correction operations in the previous steps, the corrected phase time series data for each pixel is obtained. The residual phases in the corrected phase time series are evaluated. Based on the residual phase evaluation, the corrected phase time series is further optimized and corrected to obtain the target phase time series data.

[0049] S106: Based on the corrected phase time series data, the surface deformation information of the study area is obtained.

[0050] In an exemplary embodiment, the surface deformation information of the study area is obtained based on the corrected phase time series data. Specifically, this includes: evaluating the topographic phase residuals in the corrected phase time series data to obtain evaluation results; the evaluation results include the topographic phase residual values ​​of each phase cell, the statistical distribution characteristics of the topographic phase residuals, the spatial dispersion of the topographic phase residuals, and the results of determining if the topographic phase residuals exceed the limits; based on the evaluation results, the phase time series data are optimized by using trend surface fitting correction, spatial interpolation correction, and residual compensation correction, respectively, and a second residual verification is performed on the optimized phase time series data to remove abnormal phase data that fail the verification, thereby obtaining the target phase time series data for each phase cell; based on the target phase time series data of all phase cells, the surface deformation information of the study area is inverted.

[0051] In one exemplary embodiment, the present invention provides as follows Figure 2 The flowchart shown below illustrates the InSAR interferometric pair construction method. Figure 2 As shown, SAR data preprocessing is performed using the open-source software ISCE (InSAR Scientific Computing Environment), and time series analysis is conducted using the MintPy (Miami InSAR Time-series software in Python) platform. The implementation process includes the following steps:

[0052] ISCE Experimental Procedure: Basic Data Preparation 1. Use the network-based enhanced spectral diversity (NESD) method for refined registration of SAR images, setting the number of connections to 5 to balance data processing quality and efficiency; 2. Use the registered images to generate differential interferograms using interferometry techniques; 3. Perform multi-view processing on the results obtained from the differential interferometry; 4. Based on data such as the digital elevation model (DEM), remove phase components caused by flat terrain effects and topographic relief from the interferometric phase; 5. Set the number of connections to balance data processing quality and efficiency, and set the filter strength to 0.2 for phase filtering; 6. Use the minimum cost flow (MCF) algorithm for phase unwrapping, etc.

[0053] MintPy experimental steps: 1. Construct an interferogram network based on the processed single-pair interferograms; 1.1. Based on existing research and theory, tropospheric atmospheric delay inversion and component evaluation are carried out. The phase standard deviation of the unwrapped interferogram is used to evaluate the short-range scale components of the tropospheric atmosphere. The elevation phase correlation coefficient is used to evaluate the topographically related tropospheric atmospheric components. The spatial correlation disappearance distance (effective distance) of the semivariogram optimal fitting model is used to evaluate the wide-area scale components. 1.2. Using the three-standard-deviation principle, interferograms significantly affected by air pollution are screened out from each evaluation factor dimension to achieve preliminary filtering of interferograms; 1.3 For the interferograms after preliminary screening, the evaluation factors are first standardized and principal component analysis is performed. The weights of each factor are calculated based on the cumulative variance contribution rate, and a comprehensive scoring framework with fused weights is established. Then, based on the statistical distribution of the score values, the interferogram with the smallest atmospheric delay is selected, and a high-quality interferogram network is constructed to complete the ground deformation inversion.

[0054] 2. Based on the constructed interferogram network, perform interferometric network inversion. When constructing the first interferogram network, use a mask to reduce spatial coherence and select high-quality pixel target points with temporal coherence greater than 0.6 for interferometric network inversion;

[0055] 3. For each pixel in the interferogram network, solve for its corresponding phase value; 4. Considering the influence of Earth's solid tides on the surface phase, solid tide phase correction is performed on the pixel phase based on the theoretical model and relevant parameters of solid tides; 5. Atmospheric delayed phase was corrected using the fifth-generation meteorological reanalysis dataset method (ERA5 method) published by the European Centre for Medium-Range Weather Forecasts (ECMWF); 6. Based on the previous removal of flat land and terrain phases, further corrections are made to the terrain phase residuals; 7. After a series of correction operations in the previous steps, the corrected phase timing data of each pixel is obtained; 8. Evaluate the residual phase in the corrected phase timing; 9. Based on the residual phase assessment, the corrected phase time series is further optimized and corrected to obtain the final phase time series data, and finally high-precision surface deformation information is obtained.

[0056] In one exemplary embodiment, a comparative experiment was conducted using the present invention and the most commonly used spatiotemporal baseline thresholding method. The results showed that the present invention obtained a larger mean coherence coefficient of the interferometric image pair, increased the number of high-quality pixels calculated by 6.5%, had a stronger detection capability for mountainous vegetation-covered areas, increased the coefficient of determination R2 between the cumulative deformation inversion results and GPS deformation results by 36.9%, reduced the root mean square error (RMSE) by 46.1%, and reduced the mean absolute error (MAE) by 44.7%.

[0057] In one exemplary embodiment, the present invention uses a specific region as the study area. This region is a mixed plain and hilly area, extensively covered by urban and mountainous vegetation. Forty-two Sentinel-1A SAR images covering the study area from November 2021 to March 2023 are used as the experimental database. Simultaneously, the present invention selects a 30-meter resolution Digital Elevation Model (DEM) as the topographic data. The results compared with the spatiotemporal baseline threshold network method are as follows: the tropospheric atmospheric delayed principal component weight reconstruction network model has a greater number of high-quality pixels retrieved from time-series InSAR in mountainous vegetation-covered areas, and its detection capability in mountainous vegetation-covered areas is stronger. The results show that the tropospheric atmospheric delayed principal component weight reconstruction network model has a higher mean coherence coefficient of interferometric pairs, and the selected interferograms are of better quality. The tropospheric atmospheric delayed principal component weight reconstruction network method identifies and extracts more pixels, increasing the number of pixels by 6.8% compared with the spatiotemporal baseline threshold network model, and has a larger spatial coverage. In addition, the tropospheric atmospheric delayed principal component weight reconstruction network method has more pixels with a rate standard deviation of less than 7 mm / a, and the number of high-quality pixels is increased by 6.5% compared with the spatiotemporal baseline threshold network method.

[0058] The key feature of this invention is the proposed method for surface deformation inversion that takes atmospheric effects into account. Its core lies in constructing a principal component weighted interferogram network model that considers tropospheric atmospheric delay, achieving high-quality selection of InSAR interferometric pairs. Compared with existing technologies, this invention has the advantages of minimal tropospheric atmospheric delay and lowest pollution level in the interferometric pairs, resulting in higher interferogram quality.

[0059] When applying the surface deformation inversion method that takes atmospheric effects into account provided by this invention, it is not necessary to consider... Figure 1 The steps shown are executed in sequence. The specific execution order of each step can be determined as needed, and this invention does not impose any restrictions on it.

[0060] The above describes one or more embodiments of the present invention for a surface deformation inversion method that takes atmospheric effects into account. Based on the same idea, the present invention also provides a corresponding surface deformation inversion device that takes atmospheric effects into account, such as... Figure 3 As shown.

[0061] Figure 3 A schematic diagram of a surface deformation inversion device that takes atmospheric effects into account, provided by the present invention, includes: The acquisition module 301 is used to acquire time-series SAR images of the study area within a preset time period, and to combine multiple SAR images in the time-series SAR images by interferometric image pairing to generate multiple initial differential interferograms.

[0062] The filtering module 302 is used to filter out differential interferograms affected by atmospheric pollution from the multiple generated initial differential interferograms, thereby obtaining multiple differential interferograms.

[0063] Principal component analysis module 303 is used to perform principal component analysis on the evaluation factors of each differential interferogram to obtain the weight of each evaluation factor; the evaluation factors include phase standard deviation, elevation phase correlation coefficient and spatial correlation vanishing distance.

[0064] Module 304 is used to perform a weighted summation of the values ​​of all evaluation factors and their weights to obtain the comprehensive score of the differential interferogram. The differential interferograms are sorted in descending order of comprehensive score. A preset number of differential interferograms are selected as target differential interferograms, and an interferogram network is constructed based on the target differential interferograms.

[0065] The inversion module 305 is used to perform interferometric network inversion based on the interferogram network to obtain the phase value of each pixel in the interferogram network. For each pixel, the phase value of the pixel is sequentially corrected by solid tide phase correction, atmospheric delay phase correction, and topographic phase residual correction to obtain the corrected phase time series data.

[0066] The determination module 306 is used to obtain surface deformation information of the study area based on the corrected phase time series data.

[0067] Specific limitations regarding the surface deformation inversion device that takes atmospheric effects into account can be found in the limitations of the surface deformation inversion method that takes atmospheric effects into account, as described above, and will not be repeated here. Each module in the aforementioned surface deformation inversion device that takes atmospheric effects into account can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in the processor of a computer device in hardware form or independently of the processor, or they can be stored in the memory of a computer device in software form, so that the processor can call and execute the corresponding operations of each module.

[0068] The present invention also provides a computer-readable storage medium storing a computer program that can be used to execute the above-described... Figure 1 The provided method for inverting surface deformation taking into account atmospheric effects.

[0069] The present invention also provides Figure 4 The schematic diagram of the computer device shown is as follows: Figure 4 As shown, at the hardware level, this computer device includes a processor, internal bus, network interface, memory, and non-volatile memory, and may also include other hardware required for business operations. The processor reads the corresponding computer program from the non-volatile memory into memory and then executes it to achieve the above. Figure 1 The provided method for inverting surface deformation taking into account atmospheric effects.

[0070] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the methods described above. Any references to memory, storage, databases, or other media used in the embodiments provided by this invention can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical storage, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0071] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this invention.

Claims

1. A method for surface deformation inversion considering atmospheric effects, characterized in that, The application relates to a method for obtaining surface deformation information of a research area. The method comprises the following steps: acquiring time-series SAR images of the research area in a preset time period, and combining multiple SAR images in the time-series SAR images to generate multiple initial differential interferograms; screening out differential interferograms affected by atmospheric pollution from the generated multiple initial differential interferograms to obtain multiple differential interferograms; for each differential interferogram, performing principal component analysis on evaluation factors of the differential interferogram to obtain the weight of each evaluation factor; the evaluation factors include a phase standard deviation, a height phase correlation coefficient and a spatial correlation disappearance distance; weighting and summing the values of all the evaluation factors and the weights of the evaluation factors to obtain a comprehensive score of the differential interferogram, sorting the differential interferograms according to the comprehensive scores from large to small, screening out a preset number of differential interferograms as target differential interferograms, and constructing an interferogram network based on the target differential interferograms; performing interferogram network inversion based on the interferogram network to obtain the phase value of each pixel in the interferogram network, and sequentially performing solid tide phase correction, atmospheric delay phase correction and terrain phase residual correction on the phase value of each pixel to obtain corrected phase time-series data; 2. The method of claim 1, wherein, obtaining surface deformation information of the research area based on the corrected phase time-series data. The method comprises the following steps: for each SAR image in the time-series SAR images, performing registration on the SAR image by using a sequential network enhancement spectrum diversity method; the connection number in the registration process is set as a preset connection number; 3. The method of claim 1, wherein, generating an initial differential interferogram by using the registered SAR image. The method comprises the following steps: sequentially performing multi-view processing, phase correction, phase filtering and phase unwrapping on the initial differential interferogram; determining a short-distance scale component of the initial differential interferogram by using the phase standard deviation of the unwrapped initial differential interferogram; determining a terrain correlation component of the initial differential interferogram by using a height phase correlation coefficient; determining a wide-area scale component of the initial differential interferogram based on a spatial correlation disappearance distance of an optimal fitting model of a semi-variogram function; 4. The method of claim 3, wherein, based on the short-distance scale component, the terrain correlation component and the wide-area scale component, screening out differential interferograms affected by atmospheric pollution from the generated multiple initial differential interferograms by using a three-standard-deviation principle to obtain multiple differential interferograms. The method comprises the following steps: performing multi-view processing on the initial differential interferogram; based on digital elevation model data, removing phase components caused by flat ground effects and terrain undulations from the initial interferometric phase image after multi-view processing to perform phase correction; performing phase filtering on the result of the phase correction; 5. The method of claim 1, wherein, performing phase unwrapping on the result of the phase filtering by using a minimum cost flow algorithm to obtain the unwrapped initial differential interferogram. The method comprises the following steps: The terrain phase residual in the corrected phase time series data is evaluated to obtain an evaluation result; the evaluation result includes terrain phase residual values of each phase element, statistical distribution characteristics of the terrain phase residual, spatial dispersion of the terrain phase residual, and terrain phase residual out-of-limit determination results; Based on the evaluation result, the phase time series data is optimized by using trend surface fitting correction, spatial interpolation correction, and residual compensation correction, and the optimized phase time series data is subjected to secondary residual verification to eliminate abnormal phase element data that fails the verification, thereby obtaining target phase time series data of each phase element; Based on the target phase time series data of all phase elements, the surface deformation information of the research area is inverted.

6. The method of claim 1, wherein, The evaluation factors of the differential interferogram are subjected to principal component analysis to obtain the weight of each evaluation factor, specifically including: The Pearson correlation coefficient between any two evaluation factors in the differential interferogram is calculated, and a correlation coefficient matrix is constructed based on the calculated Pearson correlation coefficient; A characteristic equation is constructed based on the correlation coefficient matrix, and the characteristic equation is solved to obtain the characteristic value corresponding to each evaluation index; Based on the plurality of characteristic values, the cumulative variance contribution rate of each evaluation factor is calculated, and a plurality of principal components with a cumulative variance contribution rate greater than or equal to a preset percentage are selected; The linear combination coefficients of each evaluation index in different principal components are calculated; The cumulative variance contribution rate of each principal component is calculated; Based on the cumulative variance contribution rate of the principal component and the linear combination coefficients of each evaluation factor in different principal components, the initial weight of each evaluation factor is calculated; The initial weight is normalized to obtain the weight of each evaluation factor.

7. The method of claim 6, wherein, The correlation coefficient matrix is: ； wherein R is a correlation coefficient matrix, is the Pearson correlation coefficient between the i-th evaluation factor and the j-th evaluation factor in the differential interferogram, is the Pearson correlation coefficient between the i-th evaluation factor and the j-th evaluation factor in the differential interferogram, is the Pearson correlation coefficient between the i-th evaluation factor and the j-th evaluation factor in the differential interferogram, p is the number of evaluation factors; The characteristic equation is: ； wherein is an eigenvalue, E is p is an identity matrix of order The formula for calculating the cumulative variance contribution rate of the first evaluation factor is: k ​ ； wherein, is the cumulative variance contribution rate of the first k evaluation factor, is the eigenvalue of the first k evaluation factor; The calculation formula of the linear combination coefficient of the first evaluation factor in the first principal component is: j The calculation formula of the linear combination coefficient of the first evaluation factor in the first principal component is: i The calculation formula of the linear combination coefficient of the first evaluation factor in the first principal component is: ； wherein is the linear combination coefficient of the jth evaluation factor in the ith principal component, j is the linear combination coefficient of the jth evaluation factor in the ith principal component, i is the linear combination coefficient of the jth evaluation factor in the ith principal component, is the loading coefficient of the ith principal component to the jth evaluation factor, i is the loading coefficient of the ith principal component to the jth evaluation factor, is the loading coefficient of the ith principal component to the jth evaluation factor, , , m is the number of principal components, p is the number of evaluation factors; The formula for calculating the cumulative variance contribution rate of the first principal component is: i ​ ； wherein is the cumulative variance contribution rate of the first i principal component; The formula for calculating the initial weight of the first evaluation factor is: j ​ ； wherein is the initial weight for the j th evaluation factor, ; The calculation formula of the weight of the first evaluation factor is: j the second evaluation factor is: ； in, For the first j The weights of each evaluation factor For the first k The initial weights of each evaluation factor.

8. A device for surface deformation inversion taking into account atmospheric effects, characterized in that, including: The acquisition module is configured to acquire time-series SAR images of a research area within a preset time period, and combine a plurality of SAR images in the time-series SAR images to generate a plurality of initial differential interferograms; The screening module is configured to screen out differential interferograms affected by atmospheric pollution from the generated plurality of initial differential interferograms to obtain a plurality of differential interferograms; The principal component analysis module is configured to, for each differential interferogram, perform principal component analysis on evaluation factors of the differential interferogram to obtain the weight of each evaluation factor; the evaluation factors include phase standard deviation, elevation phase correlation coefficient, and spatial correlation distance; The construction module is configured to weight and sum the values of all evaluation factors and the weights of the evaluation factors to obtain a comprehensive score of the differential interferogram, sort the differential interferograms in descending order of the comprehensive score, select a preset number of differential interferograms as target differential interferograms, and construct an interferogram network based on the target differential interferograms; The inversion module is configured to perform interferometric network inversion based on the interferogram network to obtain phase values of each pixel in the interferogram network, and sequentially perform solid tide phase correction, atmospheric delay phase correction, and terrain phase residual correction on the phase values of each pixel to obtain corrected phase time series data; The determination module is configured to obtain surface deformation information of the research area based on the corrected phase time series data.

9. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, and the computer program is executed by a processor to implement the method in any one of claims 1-7.

10. A computer device, comprising: A computer program product, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein the processor implements the method in any one of claims 1-7 when executing the program.

Images