A method and system for monitoring ground deformation

Abstract

The present application relates to a method and system for monitoring surface deformation. The method comprises: extracting homogeneous points from the SAR sequence image data of the monitoring area using the Gamma confidence interval discrimination method to generate a homogeneous point set of the SAR sequence image data; calculating the interference phase coherence coefficient between the secondary image and the main image in the SAR sequence image data based on the homogeneous point set; performing adaptive multi-view processing on the interference pattern of the SAR sequence image data according to the interference coherence coefficient to obtain the multi-view interference phase; screening PS candidate points based on the amplitude discreteness of the SAR sequence image to obtain PS candidate points; calculating the adaptive threshold of the interference coherence coefficient based on the homogeneous point set, screening the PS candidate points using the adaptive threshold to obtain a valid PS point mask; performing atmospheric correction on the SAR sequence image data based on the valid PS point mask and the multi-view interference phase. The present application improves the monitoring accuracy of surface deformation.

Description

Technical Field

[0001] The present application belongs to the field of geological disaster monitoring technology, and in particular relates to a surface deformation monitoring method and system. Background Art

[0002] To meet the demands of rapid economic and social development, many cities have undertaken large-scale municipal engineering projects. Due to coastal reclamation, groundwater extraction, and railway and subway construction, surface deformation has become a major geological disaster in many cities, seriously endangering public safety. Currently, technologies for monitoring urban surface deformation primarily rely on field measurements (such as levels and total stations) and remote sensing. Field measurements are unsuitable for large-scale, high-frequency dynamic monitoring due to high labor costs, high deployment risks, sparse monitoring points, and limited coverage. Among remote sensing monitoring methods, optical sensors are affected by clouds, making it difficult to obtain long sequences of effective optical images for monitoring in the cloudy and rainy climates of coastal cities in South China.

[0003] Spaceborne synthetic aperture radar (SAR) boasts all-day, all-weather operation, wide coverage, strong penetration (capable of obtaining valid data through cloud cover), and a short satellite revisit period. Furthermore, much SAR data is freely available for global download. Therefore, it enables large-scale, high-frequency dynamic monitoring, offering irreplaceable advantages over traditional surface deformation monitoring. Synthetic space radar interferometry (InSAR) has become the most effective and feasible technology for measuring surface deformation. Among these, time-series SAR interferometry has achieved significant results in surface deformation monitoring. These technologies can be broadly divided into two categories: Persistent Scatterer InSAR (PS-InSAR) and Small Baseline Subset InSAR (SBAS-InSAR). In 2008, British scholar Hooper proposed a new time series InSAR algorithm, StaMPS (Stanford Method for Persistent Scatterers). Depending on application requirements, it employs either PS-InSAR (PS-InSAR) or SBAS-InSAR (SBAS-InSAR) optimization models, either alone or in combination. Currently, it is one of the most widely used technologies for surface deformation monitoring. In densely built-up urban areas, PS-InSAR within StaMPS is often used for surface deformation monitoring. The main problems with existing StaMPS PS-InSAR technology include:

[0004] 1) To maintain the original resolution of the SAR image and thus maintain a high number of PS (Persistent Scatter) candidate points, no multi-look processing was performed on the interferometric phase, and no spatially uncorrelated noise was suppressed before PS-InSAR processing.

[0005] 2) During the secondary screening of PS points, the parameter sensitivity is high. Fine-tuning the parameters at the critical point can easily lead to too few or too many PS points. Too few PS points cannot fully present the spatial variation trend of deformation, while too many PS points often incorrectly retain noise points on the water body and vegetation surface.

[0006] 3) In the atmospheric phase correction link, the existing methods of StaMPS either use the spatiotemporal filtering method alone or use the external atmospheric products alone to estimate the atmospheric delayed phase, which cannot handle the large spatiotemporal variations of water vapor in the coastal areas of South China. Summary of the Invention

[0007] The present application provides a surface deformation monitoring method and system, which aims to solve at least one of the above-mentioned technical problems in the prior art to a certain extent.

[0008] In order to solve the above problems, this application provides the following technical solutions:

[0009] A method for monitoring ground deformation, comprising:

[0010] Based on the StaMPS algorithm, the Gamma confidence interval discrimination method is used to extract homogeneous points from the SAR sequence image data of the monitoring area to generate a homogeneous point set of the SAR sequence image data;

[0011] Calculating the interference phase coherence coefficient between the secondary image and the primary image in the SAR sequence image data based on the homogeneous point set;

[0012] performing adaptive multi-look processing on the interferogram of the SAR sequence image data according to the calculation result of the interference coherence coefficient to obtain the multi-look interferometric phase;

[0013] Calculating the amplitude dispersion of the SAR sequence image, and screening PS candidate points according to the amplitude dispersion to obtain PS candidate points;

[0014] Calculating an adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and screening PS candidate points using the adaptive threshold to obtain a valid PS point mask;

[0015] Based on the effective PS point mask and the multi-look interferometric phase, the atmospheric correction of the SAR sequence image data is performed using the GACOS model and the InSAR spatiotemporal filtering algorithm;

[0016] The deformation rate and cumulative deformation amount of the atmospheric-corrected SAR image sequence data are subjected to time series inversion based on the StaMPS algorithm to obtain the surface deformation monitoring results of the monitoring area.

[0017] The technical solution adopted in the embodiment of the present application also includes: extracting homogeneous points from the SAR sequence image data of the monitoring area using the Gamma confidence interval discrimination method includes:

[0018] For each pixel in the SAR sequence image data, a search window is determined, and the pixel is used as the reference pixel. The reference sample mean is calculated based on the time dimension sampling of the pixel. by The Gamma distribution confidence interval is calculated using the set α value; the Gamma distribution confidence interval calculation formula is:

[0019]

[0020] Where 1-α is the confidence level, Represents the standard GammaG(NL,1) distribution Percentile, N is the number of images in the SAR image sequence data, L is the number of views in the SAR image sequence data, a=L; σ 2 is the backscatter intensity of SAR sequence image data;

[0021] Determine whether the time dimension sample mean of the neighborhood pixels in the search window falls into the Gamma distribution confidence interval. If so, determine that the neighborhood pixels and the reference pixel belong to the same statistical distribution, and obtain the initial homogeneous pixel subset Ω of the pixel. init ;

[0022] Using the initial homogeneous pixel subset Ω init The average value of the sample means of the K pixels in the reference sample mean Update and update according to the updated Recalculate the Gamma distribution confidence interval;

[0023] According to the new Gamma distribution confidence interval, the initial homogeneous pixel subset Ω init The pixels in the θ are judged twice, and the neighboring pixels that fall into the new Gamma distribution confidence interval and are connected to the reference pixel space are retained to generate a selected homogeneous pixel subset Ω of the pixel.

[0024] The technical solution adopted in the embodiment of the present application further includes: the calculation of the interference coherence coefficient between the secondary image and the primary image in the SAR sequence image data based on the homogeneous point set includes:

[0025] For each pixel in the SAR sequence image data, a homogeneous pixel subset Ω is taken to estimate the interference coherence coefficient to obtain a biased estimate of the interference coherence coefficient of the pixel; the interference coherence coefficient calculation formula is:

[0026]

[0027] Where Z1(l) and Z2(l) represent the complex signals corresponding to pixel l in the two SAR images Z1 and Z2 used for interference, and K is the number of pixels in the homogeneous pixel subset Ω where the pixel is located;

[0028] The double Bootstrapping method is used to correct the biased estimation of the interference coherence coefficient to obtain the unbiased interference coherence coefficient of the SAR sequence image data; the interference coherence coefficient of the approximate unbiased estimation is for:

[0029]

[0030] In which, all points in the homogeneous pixel subset Ω are represented as X = (x1, x2, ..., x K ), K is the number of homogeneous points, R is the number of times Bootstrapping repeatedly samples X, is the sample of Bootstrapping replication, and γ is the biased estimate of the interference coherence coefficient.

[0031] The technical solution adopted by the embodiment of the present application also includes: the adaptive multi-look processing of the interferogram of the SAR sequence image data according to the interference coherence coefficient calculation result includes:

[0032] Preprocessing the SAR sequence image data to generate an interference pattern of the SAR sequence image data; the preprocessing includes:

[0033] Performing a spatiotemporal baseline analysis on the SAR image sequence data, selecting an image with the shortest spatial and temporal baselines in the SAR image sequence data as the primary image, and the other images as secondary images;

[0034] According to the pairing criterion of the PS-InSAR algorithm, the secondary image is registered with the primary image using orbital parameters and a DEM model, and the secondary image and the primary image are combined into an interference pair to generate an interference map.

[0035] The technical solution adopted by the embodiment of the present application also includes: the adaptive multi-look processing of the interferogram of the SAR sequence image data according to the interference coherence coefficient calculation result includes:

[0036] For each pixel in the SAR sequence image data, taking the homogeneous pixel subset where the pixel is located as a unit, weighting the interference pattern complex signal within the homogeneous pixel subset according to the interference coherence coefficient, taking the phase thereof, and obtaining a multi-look interferometric phase;

[0037] For any pixel p, the calculation formula of the multi-view interference phase is:

[0038]

[0039] Where k is the number of homogeneous pixels in the homogeneous pixel subset where pixel P is located, Coh i is the interference coherence coefficient of the i-th pixel, Ifg i is the complex interference pattern signal of the i-th pixel.

[0040] The technical solution adopted in the embodiment of the present application further includes: the calculation of the amplitude dispersion of the SAR sequence image data is specifically as follows:

[0041]

[0042] Among them, D A represents the amplitude dispersion, σ A and μ A are the standard deviation and mean of SAR sequence image data calculated along the time dimension.

[0043] The technical solution adopted in the embodiment of the present application further includes: the adaptive threshold of the interference coherence coefficient is calculated based on the homogeneous point set, and the adaptive threshold is used to screen the PS candidate points, including:

[0044] According to the actual threshold value set And the number of homogeneous pixels in the homogeneous pixel subset, calculate the adaptive threshold:

[0045]

[0046] where σ co h is the Cramer-Rao lower bound standard deviation of the unbiased interference coherence coefficient;

[0047] For each pixel in the SAR sequence image data, the sequence of the unbiased interference coherence coefficient in the time dimension is statistically analyzed according to the adaptive threshold of the pixel. If the precise estimate of the interference coherence coefficient in the time series is If the number of is greater than the set ratio, the point is determined to be a valid PS point;

[0048] The intersection of all valid PS points and the PS candidate points is calculated to obtain a valid PS point mask.

[0049] The technical solution adopted in the embodiment of the present application further includes: performing atmospheric correction on the SAR sequence image data based on the effective PS point mask and the multi-look interferometric phase using the GACOS model and the InSAR spatiotemporal filtering algorithm, including:

[0050] Using the StaMPS algorithm to estimate the spatial correlation phase component of the SAR sequence image data and to estimate the phase noise of the effective PS point mask;

[0051] The effective PS point mask is screened in combination with the amplitude dispersion and phase noise estimation, and three-dimensional phase unwrapping is performed on the screened PS points:

[0052] Based on the phase unwrapping results of the PS points, the GACOS model and the InSAR spatiotemporal filtering algorithm are used to filter out the atmospheric delay and orbit error phase from the interferometric phase, and perform atmospheric correction on the SAR sequence image data.

[0053] The technical solution adopted in the embodiment of the present application also includes: the atmospheric correction of the SAR sequence image data using the GACOS model and the InSAR spatiotemporal filtering algorithm includes:

[0054] The GACOS model is used to estimate the atmospheric delay phase of each interferometer pair, and the atmospheric phase component related to terrain is removed from the corresponding interferometer phase.

[0055] Using the InSAR spatiotemporal filtering algorithm to perform time-dimensional low-pass filtering on the SAR sequence image data to obtain the atmospheric delay and orbit error of the main image;

[0056] The InSAR spatiotemporal filtering algorithm is used to perform high-pass filtering in the time dimension and low-pass filtering in the space dimension on the SAR data to obtain the atmospheric delay, orbit error and DEM error of the secondary image;

[0057] The atmospheric delay and orbit error of the primary image, the atmospheric delay, orbit error and DEM error of the secondary image are removed from the interferometric phase, and refined atmospheric correction is performed on the SAR sequence image data.

[0058] Another technical solution adopted in the embodiment of the present application is: a surface deformation monitoring system, comprising:

[0059] Homogeneous point extraction module: It is used to extract homogeneous points from the SAR sequence image data of the monitoring area based on the StaMPS algorithm and the Gamma confidence interval discrimination method to generate a homogeneous point set of the SAR sequence image data;

[0060] A multi-look processing module is used to calculate the interference coherence coefficient between the secondary image and the primary image in the SAR sequence image data based on the homogeneous point set, and perform adaptive multi-look processing on the interferogram of the SAR sequence image data according to the calculation result of the interference coherence coefficient to obtain the multi-look interference phase;

[0061] Initial PS point screening module: used to calculate the amplitude dispersion of the SAR sequence image, and screen PS candidate points according to the amplitude dispersion to obtain PS candidate points;

[0062] Valid PS point screening module: used to calculate the adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and use the adaptive threshold to screen the PS candidate points to obtain a valid PS point mask;

[0063] Atmospheric correction module: used for performing atmospheric correction on the SAR sequence image data based on the effective PS point mask and the multi-look interferometric phase using the GACOS model and the InSAR spatiotemporal filtering algorithm;

[0064] Sequence inversion module: used to perform time series inversion on the deformation rate and cumulative deformation of the atmospherically corrected SAR sequence image data based on the StaMPS algorithm to obtain the surface deformation monitoring results of the monitoring area.

[0065] Compared with the prior art, the beneficial effects produced by the embodiments of the present application are as follows: the surface deformation monitoring method and system of the embodiments of the present application improve the classic StaMPS algorithm based on homogeneous point extraction. Before the PS point screening, the interference phase is subjected to adaptive multi-view processing based on the homogeneous point set. By distinguishing between targets with two different scattering mechanisms, permanent scattering points PS and distributed scattering points DS, the phase noise is effectively suppressed and the high coherence signal of the PS point target is avoided from being confused by the phase noise of other points, thereby improving the subsequent solution accuracy. The interference coherence coefficient after deviation correction and the adaptive threshold are used to assist in PS point screening. While effectively removing low coherence points on the surface of water bodies and vegetation, the integrity of the spatial deformation trend is maintained as much as possible, avoiding excessive screening of PS points. By combining the GACOS atmospheric model and InSAR spatiotemporal filtering to perform fine corrections on the atmospheric delay phase, the accuracy of the atmospheric correction is improved, thereby improving the monitoring accuracy of surface deformation. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 is a flow chart of a method for monitoring ground deformation according to an embodiment of the present application;

[0067] Figure 2 This is a flowchart of atmospheric fine correction using the GACOS atmospheric model and spatiotemporal filtering in accordance with an embodiment of the present application.

[0068] Figure 3Schematic diagram of the structure of the surface deformation monitoring system according to an embodiment of the present application. DETAILED DESCRIPTION

[0069] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0070] See also Figure 1 , is a flow chart of a surface deformation monitoring method according to an embodiment of the present application. The surface deformation monitoring method according to an embodiment of the present application comprises the following steps:

[0071] S10: Acquire SAR sequence image data of the monitoring area;

[0072] In this step, the acquired SAR image sequence data is the SLC time series data of Sentinel-1IW TOPS mode. The parameters of the SLC time series data include time range, orbit direction, revisit period, polarization mode, incident angle, wavelength, slant range sampling interval, and azimuth sampling interval, as shown in Table 1:

[0073] Table 1. Sentinel-1IW SLC time series data parameters for the monitoring area

[0074]

[0075]

[0076] S20: preprocessing the SAR sequence image data to generate an interferogram of the SAR sequence image data;

[0077] In this step, the preprocessing of SAR sequence image data specifically includes:

[0078] S21: Perform spatiotemporal baseline analysis on the SAR image sequence data, select the image with the shortest spatial and temporal baseline in the SAR image sequence data as the main image, and the other images as the secondary images;

[0079] S22: According to the pairing criteria of the PS-InSAR algorithm, the secondary image is registered with the primary image using orbital parameters and a DEM (Digital Elevation Model) model, and the secondary image and the primary image are combined into an interferogram pair to generate an interferogram.

[0080] S23: After removing the flat ground phase and terrain phase of the interferogram based on the DEM model, the interferometric phase file after registration is exported.

[0081] S30: In the StaMPS algorithm, the Gamma confidence interval discriminant method is used to extract the Statistically Homogeneous Pixel (SHP) of the SAR sequence image data to generate a homogeneous point set of the SAR sequence image data;

[0082] In this step, the backscatter intensity of a pixel (i.e., resolution unit) in the SAR image is the superposition of the backscatter intensities of a large number of scatterers belonging to the same pixel. If a pixel is entirely composed of distributed scatterers, then the pixel is a DS (Distributed Scatterer) point target, which can be regarded as a zero-mean complex Gaussian random variable. Given a SAR image sequence data with N periods of SLC (Single Look Complex) format, for a DS point target, its time dimension sampling {Z1, Z2, ..., Z N} T It can be regarded as a zero-mean complex Gaussian random variable, and its covariance matrix can be defined as:

[0083]

[0084] in is the amplitude of the i-th image, the diagonal elements of the matrix is the backscatter intensity of the i-th SAR image. In the case of multiple views, assuming the number of views is L, the backscatter intensity expression is It obeys the Gamma distribution, and the distribution parameters are and L:

[0085]

[0086] For a DS point target, the time dimension vector of its backscatter intensity is expressed as The purpose of this application is to find a confidence interval for σ. According to SAR statistical theory, σ follows the Gamma distribution σ~G(a,b), where a=L, Its mean is μ=σ 2 , the variance is According to the central limit theorem, the sample mean Normal distribution The corresponding The confidence interval for is:

[0087]

[0088] in is the standard overall distribution probability density function Percentile, 1-α is the confidence level (for example, if the confidence level is 95%, then α = 5%). The neighboring pixels of the pixel are judged one by one. If the sample mean of the neighboring pixel falls within the confidence interval, the neighboring pixel is considered to belong to the same statistical distribution as the pixel, and the neighboring pixel and the pixel are classified into a homogeneous point subset.

[0089] When the time series of SAR image data is short, that is, the number of image periods N is less than a certain value, the sample mean of the Gamma distribution does not obey the normal distribution. In this case, the confidence interval is asymmetric and the central limit theorem cannot be used for inference. In this case, if a random variable follows the Gamma distribution σ~G(a,b), where a=L, but make Y follows the standard Gamma distribution Y~G(Na,1), then the corresponding Gamma distribution confidence interval is:

[0090]

[0091] in Represents the standard Gamma distribution Percentile. Because b=σ 2 / L, then a=L, substitute into (4), and after transformation, we can get The 1-α confidence interval for is:

[0092]

[0093] Where 1-α is the confidence level (for example, if the confidence level is 95%, then α = 5%), Represents the standard Gamma distribution G(NL,1) Percentile, N is the number of images in the SAR image sequence data, L is the number of views in the SAR image sequence data, σ 2 is the backscatter intensity of SAR sequence image data, is the sample mean of σ. Therefore, regardless of the length of the time series of SAR image data, formula (7) can be uniformly used as the confidence interval to identify homogeneous points.

[0094] Based on the above-mentioned principle of homogeneous point identification based on Gamma distribution confidence intervals, in the embodiment of the present application, the process of extracting homogeneous points from SAR sequence image data using the Gamma confidence interval identification method specifically includes:

[0095] S31: Using an adjusted boxplot, we detect and remove outliers in SAR image data by taking into account the medcouple (MC) indicator of data distribution skewness.

[0096] S32: using the one-sided Anderson–Darling (AD) test method to determine whether the time dimension vector of each pixel in the SAR sequence image data after removing the singular value obeys the Gamma distribution. If it obeys the Gamma distribution, execute S33;

[0097] S33: Use formula (7) as the confidence interval to identify homogeneous points of SAR sequence image data and extract the homogeneous point set of SAR sequence image data;

[0098] The specific extraction process is as follows: for each pixel in the SAR sequence image data, a search window is determined, the pixel is used as the reference pixel, and the reference sample mean is calculated based on the time dimension sampling of the pixel. by And the set α value is used to calculate the Gamma distribution confidence interval; all the neighborhood pixels in the search window are judged one by one. If the time dimension sample mean of the neighborhood pixel falls into the Gamma distribution confidence interval, it is determined that the neighborhood pixel and the reference pixel belong to the same statistical distribution, and the initial homogeneous pixel subset Ω of the pixel is obtained. init Since the length N of the time series is finite, the sample mean Therefore, the embodiment of the present application adopts the initial homogeneous pixel subset Ω init The average of the sample means of the K pixels in the reference sample mean Update, and then according to the updated Recalculate the Gamma distribution confidence interval, and calculate the initial homogeneous pixel subset Ω according to the new Gamma distribution confidence interval init The pixels in the SAR image data are judged twice, and the neighboring pixels that fall into the new Gamma distribution confidence interval and are connected to the reference pixel space are retained, thereby generating a selected homogeneous pixel subset Ω of the pixel. The above processing is performed on each pixel in the SAR image data sequence to generate a homogeneous point set for all pixels. In the embodiment of the present application, the initial homogeneous pixel subset Ω is generated. init When generating a selected homogeneous pixel subset Ω, set α = 50%; when generating a selected homogeneous pixel subset Ω, set α = 5%. The specific α value can be set according to the actual application.

[0099] S40: Estimate the interferometric phase coherence between the secondary image and the primary image in the SAR image sequence data based on a homogeneous point set, and use a double bootstrapping (with replacement sampling) method to correct the bias of the interferometric phase coherence to obtain the unbiased interferometric coherence coefficient of the SAR image sequence data;

[0100] In this step, interferometric phase coherence is a measure of the similarity between two SAR image signals, which is used to measure the accuracy of InSAR phase estimation and guide phase filtering and phase unwrapping. In the covariance matrix in formula (3), the non-diagonal elements That is coherence. When the number of views is L, the coherence calculation formula using the sliding window is:

[0101]

[0102] The amplitude value is the interference coherence coefficient, Z i (l) is the complex signal of pixel l. When using a sliding window for coherence calculation, the implicit assumption is that adjacent pixels within a window all belong to the same statistical distribution. However, in a scene with rich textures, the pixels within a window are likely to have different statistical distributions. Indifferent averaging of neighboring pixels with different scattering characteristics will lead to deviations in the estimation of the interference phase coherence and a decrease in image resolution. Therefore, the embodiment of the present invention uses a homogeneous point set instead of a traditional sliding window to estimate the interference coherence coefficient, and uses a double Bootstrapping (sampling with replacement) method to further correct the deviation, thereby improving the estimation accuracy of the interference coherence coefficient.

[0103] Specifically, interferometric phase coherence estimation and deviation correction of SAR sequence image data based on a homogeneous point set include:

[0104] S41: for each pixel, taking a subset of homogeneous pixels Ω to estimate the interference coherence coefficient, and obtaining a biased estimate of the interference coherence coefficient of the pixel;

[0105] Among them, if the sliding window in formula (8) is replaced by a set of homogeneous points, the calculation formula of the interference coherence coefficient is:

[0106]

[0107] where Z1(l) and Z2(l) represent the complex signal corresponding to pixel l in the two interferometric SAR images Z1 and Z2, and K is the number of pixels within the homogeneous pixel subset Ω in which the pixel resides. This step can avoid estimation bias caused by texture differences across the SAR images.

[0108] S42: A non-parametric double bootstrapping method is used to correct the biased estimation of the interferometric coherence coefficient and obtain the unbiased interferometric coherence coefficient;

[0109] Among them, for the two SAR images involved in the interference, the complex signal corresponding to each pixel in the homogeneous point set Ω can be regarded as a pair of complex observation values x i =(z 1i ,z 2i ); then all points in the homogeneous point set Ω can be expressed as X=(x1,x2,…,x K ), K is the number of homogeneous points, which can be regarded as a Bootstraping sample That is, we sample K times with replacement from the original observation X, which is the first re-sampling. After re-sampling X R times, we can get a set of double bootstrapping samples. After calculating the interference coherence coefficient for each sample using formula (9), we can get R bootstrapping replicated samples: The bootstrapping replicates provide a biased estimate of the interferometric coherence coefficient. Perform bias correction to obtain an approximately unbiased estimate of the interference coherence coefficient

[0110]

[0111] S50: using a homogeneous point set, performing adaptive multi-look processing on the interferogram of the SAR sequence image data according to the calculation result of the interference coherence coefficient, and obtaining the multi-look interferometric phase;

[0112] In this step, since there are surfaces with different scattering characteristics such as vegetation, water bodies, and artificial buildings in the study area, in order to avoid confusing scatterers with different characteristics and avoid the phase noise of scatterers with low coherence from interfering with scatterers with high coherence during multi-look processing, the embodiment of the present application uses a homogeneous point set instead of a regular sliding window to perform adaptive multi-look processing on all interference patterns in the time series. By distinguishing between targets with two different scattering mechanisms, permanent scattering points PS and distributed scattering points DS, the phase noise is effectively suppressed, and the high coherence signal of the PS point target is avoided from being confused by the phase noise of other points, thereby improving the subsequent solution accuracy.

[0113] Specifically, for each pixel, the interferogram complex signal is weighted according to the precisely estimated interference coherence coefficient within the homogeneous pixel subset, and its phase is taken. This operation is repeated for each pixel in the SAR sequence image data to obtain the interferometric phase after adaptive multi-look. For any pixel p, the multi-look interferometric phase is calculated as follows:

[0114]

[0115] Where K is the number of homogeneous pixels in the homogeneous pixel subset where the pixel is located, Coh i is the interference coherence coefficient of the i-th pixel, Ifg i is the complex interference pattern signal of the i-th pixel.

[0116] Based on the above, for each pixel, the interferometric phase is weighted by the coherence coefficient of the homogeneous pixel subset to which it belongs. The PS point is located within its own homogeneous pixel subset and is not affected by the phase noise of other target points. Meanwhile, the DS point performs phase averaging weighted by the coherence coefficient within its homogeneous pixel subset, significantly reducing phase noise. Compared with the prior art, the embodiment of the present application suppresses phase noise while preserving local detail and avoiding the need to reduce the number of PS points.

[0117] S60: Calculate the amplitude dispersion of the SAR sequence image data, and screen PS candidate points of the SAR sequence image data according to the amplitude dispersion to obtain PS candidate points;

[0118] In this step, the amplitude dispersion calculation formula is as follows:

[0119]

[0120] Among them, D A represents the amplitude dispersion, σ A and μ A are the standard deviation and mean of the SAR sequence image data calculated along the time dimension. A The threshold is set to 0.35~0.4.

[0121] S70: Calculating an adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and screening PS candidate points using the adaptive threshold to obtain a valid PS point mask;

[0122] In this step, firstly, according to the set real threshold And the number of homogeneous pixels in the homogeneous point set, calculate an adaptive threshold:

[0123]

[0124] where σ coh is the Cramer-Rao lower bound standard deviation of the unbiased interference coherence coefficient, and is calculated as follows:

[0125]

[0126] For each pixel in the SAR sequence image data, the adaptive threshold of the pixel is calculated according to formula (13), and the sequence formed by the interference coherence coefficient calculated based on the homogeneous point set and corrected for deviation in the time dimension is statistically analyzed. If the accurate estimate of the interference coherence coefficient in the time series is If the number of pixels is greater than the set ratio (which can be adjusted according to the situation, such as 85%), the point is determined to be a valid PS point. After analyzing all pixels one by one, the intersection with the PS candidate points is calculated to obtain the valid PS point mask.

[0127] Based on the above, the embodiment of the present application uses the interference coherence coefficient after deviation correction and the adaptive threshold to assist in PS point screening. While effectively removing low coherence points on the surface of water bodies and vegetation, it tries to maintain the integrity of the spatial deformation trend and avoids excessive screening of PS points.

[0128] S80: Use the StaMPS algorithm to estimate the spatial correlation phase component of the SAR image sequence data, estimate the phase noise of the effective PS point mask, and obtain the terrain residual phase estimate;

[0129] S90: Combine amplitude dispersion and phase noise estimation to screen the effective PS point mask, and perform three-dimensional (spatial two-dimensional + temporal dimension) phase unwrapping on the screened PS points;

[0130] S100: Based on the phase unwrapping results of the PS point, the GACOS model and InSAR spatiotemporal filtering algorithm are used to filter out atmospheric delay, orbit error phase, etc. from the interferometric phase, and perform refined atmospheric correction on the SAR sequence image data;

[0131] In this step, the GACOS model is first used to estimate the atmospheric delay phase of each interferometer pair, and the atmospheric phase components related to the terrain are removed. The spatiotemporal filtering method in StaMPS is then used to filter out the residual atmospheric phase, thereby improving the accuracy of atmospheric correction and thus improving the monitoring accuracy of surface deformation. Figure 2 FIG. 1 is a flowchart of a refined atmospheric correction method combining the GACOS atmospheric model and the InSAR spatiotemporal filtering algorithm according to an embodiment of the present application, which specifically includes the following steps:

[0132] S101: First, based on the GACOS atmospheric model, the atmospheric delay phase of each SAR data acquisition date is estimated, and the atmospheric phase component related to terrain is removed from the corresponding interferometric phase to complete the preliminary atmospheric correction of the SAR image sequence data;

[0133] S102: Using the InSAR spatiotemporal filtering algorithm to perform time-dimensional low-pass filtering on the SAR data to obtain the atmospheric delay and orbit error of the main image;

[0134] S103: Using the InSAR spatiotemporal filtering algorithm, the SAR data is subjected to high-pass filtering in the time dimension and low-pass filtering in the space dimension to obtain the atmospheric delay, orbit error, and DEM error of other secondary images;

[0135] S104: removing the atmospheric delay and orbit error of the primary image, the atmospheric delay, orbit error and DEM error of the secondary image from the interferometric phase, and completing the refined atmospheric correction of the SAR sequence image data.

[0136] S110: Based on the StaMPS algorithm, the annual deformation rate and cumulative deformation of the SAR image sequence data after fine atmospheric correction are inverted to obtain the surface deformation monitoring results of the monitoring area.

[0137] Based on the above, the surface deformation monitoring method of the embodiment of the present application proposes the SHP StaMPS algorithm based on homogeneous point extraction, and improves the classic StaMPS algorithm. Before the PS point screening, the interference phase is subjected to adaptive multi-view processing based on the homogeneous point set. By distinguishing between targets with two different scattering mechanisms, permanent scattering points PS and distributed scattering points DS, it not only effectively suppresses the phase noise, but also avoids the high coherence signal of the PS point target being confused by the phase noise of other points, thereby improving the subsequent solution accuracy. The interference coherence coefficient after deviation correction and the adaptive threshold are used to assist in PS point screening. While effectively removing low coherence points on the surface of water bodies and vegetation, the integrity of the spatial deformation trend is maintained as much as possible, and excessive screening of PS points is avoided. By combining the GACOS atmospheric model and InSAR spatiotemporal filtering to perform fine correction of the atmospheric delay phase, the accuracy of the atmospheric correction is improved, thereby improving the monitoring accuracy of surface deformation.

[0138] See also Figure 3 , is a schematic diagram of the structure of a surface deformation monitoring system according to an embodiment of the present application. The surface deformation monitoring system 40 according to an embodiment of the present application includes:

[0139] Homogeneous point extraction module 41: used for extracting homogeneous points from the SAR sequence image data of the monitoring area based on the StaMPS algorithm and using the Gamma confidence interval discrimination method to generate a homogeneous point set of the SAR sequence image data;

[0140] Multi-look processing module 42: used to calculate the interference coherence coefficient of the SAR sequence image data based on the homogeneous point set, and perform adaptive multi-look processing on the interferogram of the SAR sequence image data according to the calculation result of the interference coherence coefficient to obtain the multi-look interferometric phase;

[0141] Initial PS point screening module 43: used to calculate the amplitude dispersion of the SAR sequence image data, and screen the PS candidate points of the SAR sequence image data according to the amplitude dispersion to obtain PS candidate points;

[0142] Valid PS point screening module 44: used for calculating the adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and screening the PS candidate points using the adaptive threshold to obtain a valid PS point mask;

[0143] Atmospheric correction module 45: configured to perform atmospheric correction on the SAR sequence image data based on the effective PS point mask and the multi-look interferometric phase using a GACOS model and an InSAR spatiotemporal filtering algorithm;

[0144] The sequence inversion module 46 is used to perform time series inversion on the deformation rate and cumulative deformation of the atmospherically corrected SAR sequence image data based on the StaMPS algorithm to obtain the surface deformation monitoring results of the monitoring area.

[0145] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application is not limited to the embodiments shown herein, but is intended to encompass the broadest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for monitoring ground deformation, characterized in that: include: Based on the StaMPS algorithm, the Gamma confidence interval discrimination method is used to extract homogeneous points from the SAR sequence image data of the monitoring area to generate a homogeneous point set of the SAR sequence image data; Calculating the interference coherence coefficient between the secondary image and the primary image in the SAR sequence image data based on the homogeneous point set; performing adaptive multi-look processing on the interferogram of the SAR sequence image data according to the calculation result of the interference coherence coefficient to obtain the multi-look interferometric phase; Calculating the amplitude dispersion of the SAR sequence image, and screening PS candidate points according to the amplitude dispersion to obtain PS candidate points; Calculating an adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and screening PS candidate points using the adaptive threshold to obtain a valid PS point mask; Based on the effective PS point mask and the multi-look interferometric phase, the atmospheric correction of the SAR sequence image data is performed using the GACOS model and the InSAR spatiotemporal filtering algorithm; The deformation rate and cumulative deformation of the atmospherically corrected SAR image sequence data are inverted in time series based on the StaMPS algorithm to obtain the surface deformation monitoring results of the monitoring area; wherein: The method of extracting homogeneous points from the SAR sequence image data of the monitoring area using the Gamma confidence interval discrimination method includes: For each pixel in the SAR sequence image data, a search window is determined, and the pixel is used as the reference pixel. The reference sample mean is calculated based on the time dimension sampling of the pixel. by The Gamma distribution confidence interval is calculated using the set α value; the Gamma distribution confidence interval calculation formula is: Where 1-α is the confidence level, Represents the standard GammaG(NL,1) distribution Percentile, N is the number of images in the SAR image sequence data, L is the number of views in the SAR image sequence data, a=L; σ 2 is the backscatter intensity of SAR sequence image data; Determine whether the time dimension sample mean of the neighborhood pixels in the search window falls into the Gamma distribution confidence interval. If so, determine that the neighborhood pixels and the reference pixel belong to the same statistical distribution, and obtain the initial homogeneous pixel subset Ω of the pixel. init ; Using the initial homogeneous pixel subset Ω init The average value of the sample means of the K pixels in the reference sample mean Update and update according to the updated Recalculate the Gamma distribution confidence interval; According to the new Gamma distribution confidence interval, the initial homogeneous pixel subset Ω init The pixels in the θ are judged twice, and the neighboring pixels that fall into the new Gamma distribution confidence interval and are connected to the reference pixel space are retained to generate a selected homogeneous pixel subset Ω of the pixel.

2. The surface deformation monitoring method according to claim 1, characterized in that: The step of calculating the interference coherence coefficient between the secondary image and the primary image in the SAR sequence image data based on the homogeneous point set includes: For each pixel in the SAR sequence image data, a homogeneous pixel subset Ω is taken to estimate the interference coherence coefficient to obtain a biased estimate of the pixel's interference coherence coefficient; the interference coherence coefficient calculation formula is: Where Z1(l) and Z2(l) represent the complex signals corresponding to pixel l in the two SAR images Z1 and Z2 used for interference, and K is the number of pixels in the homogeneous pixel subset Ω where the pixel is located; The double bootstrapping method is used to correct the biased estimation of the interference coherence coefficient to obtain the unbiased interference coherence coefficient of the SAR sequence image data; the interference coherence coefficient of the approximate unbiased estimation is obtained. for: In which, all points in the homogeneous pixel subset Ω are represented as X = (x1, x2, ..., x K ), K is the number of homogeneous points, R is the number of times Bootstrapping repeatedly samples X, is the sample of Bootstrapping replication, and γ is the biased estimate of the interference coherence coefficient.

3. The surface deformation monitoring method according to claim 2, characterized in that: The adaptive multi-view processing of the interference pattern of the SAR sequence image data according to the interference coherence coefficient calculation result includes: Preprocessing the SAR sequence image data to generate an interference pattern of the SAR sequence image data; the preprocessing includes: Performing a spatiotemporal baseline analysis on the SAR image sequence data, selecting an image with the shortest spatial and temporal baselines in the SAR image sequence data as the primary image, and the other images as secondary images; According to the pairing criterion of the PS-InSAR algorithm, the secondary image is registered with the primary image using orbital parameters and a DEM model, and the secondary image and the primary image are combined into an interference pair to generate an interference map.

4. The surface deformation monitoring method according to claim 3, characterized in that: The adaptive multi-view processing of the interference pattern of the SAR sequence image data according to the interference coherence coefficient calculation result includes: For each pixel in the SAR sequence image data, taking the homogeneous pixel subset where the pixel is located as a unit, weighting the interference pattern complex signal within the homogeneous pixel subset according to the interference coherence coefficient, taking the phase thereof, and obtaining a multi-look interferometric phase; For any pixel p, the calculation formula of the multi-view interference phase is: Where K is the number of homogeneous pixels in the homogeneous pixel subset where pixel p is located, Coh i is the interference coherence coefficient of the i-th pixel, Ifg i is the complex interference pattern signal of the i-th pixel.

5. The surface deformation monitoring method according to any one of claims 1 to 4, characterized in that: The calculation of the amplitude dispersion of the SAR sequence image data is specifically as follows: Among them, D A represents the amplitude dispersion, σ A and μ A are the standard deviation and mean of SAR sequence image data calculated along the time dimension.

6. The surface deformation monitoring method according to claim 5, characterized in that: Calculating the adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and screening the PS candidate points using the adaptive threshold comprises: According to the actual threshold value set And the number of homogeneous pixels in the homogeneous pixel subset, calculate the adaptive threshold: where σ coh is the Cramer-Rao lower bound standard deviation of the unbiased interference coherence coefficient; For each pixel in the SAR sequence image data, the sequence of the unbiased interference coherence coefficient in the time dimension is statistically analyzed according to the adaptive threshold of the pixel. If the precise estimate of the interference coherence coefficient in the time series is If the number of is greater than the set ratio, the point is determined to be a valid PS point; The intersection of all valid PS points and the PS candidate points is calculated to obtain a valid PS point mask.

7. The surface deformation monitoring method according to claim 6, characterized in that: The atmospheric correction of the SAR sequence image data using the GACOS model and the InSAR spatiotemporal filtering algorithm based on the effective PS point mask and the multi-look interferometric phase comprises: Using the StaMPS algorithm to estimate the spatial correlation phase component of the SAR sequence image data and to estimate the phase noise of the effective PS point mask; The effective PS point mask is screened in combination with the amplitude dispersion and phase noise estimation, and three-dimensional phase unwrapping is performed on the screened PS points: Based on the phase unwrapping results of the PS points, the GACOS model and the InSAR spatiotemporal filtering algorithm are used to filter out the atmospheric delay and orbit error phase from the interferometric phase, and perform atmospheric correction on the SAR sequence image data.

8. The surface deformation monitoring method according to claim 7, characterized in that: The atmospheric correction of the SAR sequence image data using the GACOS model and the InSAR spatiotemporal filtering algorithm comprises: The GACOS model is used to estimate the atmospheric delay phase of each interferometer pair, and the atmospheric phase component related to terrain is removed from the corresponding interferometer phase. Using the InSAR spatiotemporal filtering algorithm to perform time-dimensional low-pass filtering on the SAR sequence image data to obtain the atmospheric delay and orbit error of the main image; The InSAR spatiotemporal filtering algorithm is used to perform high-pass filtering in the time dimension and low-pass filtering in the space dimension on the SAR data to obtain the atmospheric delay, orbit error and DEM error of the secondary image; The atmospheric delay and orbit error of the primary image, the atmospheric delay, orbit error and DEM error of the secondary image are removed from the interferometric phase, and refined atmospheric correction is performed on the SAR sequence image data.

9. A surface deformation monitoring system using the surface deformation monitoring method according to claim 1, characterized in that: include: Homogeneous point extraction module: It is used to extract homogeneous points from the SAR sequence image data of the monitoring area based on the StaMPS algorithm and the Gamma confidence interval discrimination method to generate a homogeneous point set of the SAR sequence image data; A multi-look processing module is used to calculate the interference coherence coefficient between the secondary image and the primary image in the SAR sequence image data based on the homogeneous point set, and perform adaptive multi-look processing on the interferogram of the SAR sequence image data according to the calculation result of the interference coherence coefficient to obtain the multi-look interference phase; Initial PS point screening module: used to calculate the amplitude dispersion of the SAR sequence image, and screen PS candidate points according to the amplitude dispersion to obtain PS candidate points; Valid PS point screening module: used to calculate the adaptive threshold of the interference coherence coefficient based on the homogeneous point set, and use the adaptive threshold to screen the PS candidate points to obtain a valid PS point mask; Atmospheric correction module: used for performing atmospheric correction on the SAR sequence image data based on the effective PS point mask and the multi-look interferometric phase using the GACOS model and the InSAR spatiotemporal filtering algorithm; Sequence inversion module: used to perform time series inversion on the deformation rate and cumulative deformation of the atmospherically corrected SAR sequence image data based on the StaMPS algorithm to obtain the surface deformation monitoring results of the monitoring area.

Images