Method for analyzing facility deformation based on PSInSAR and SqueeSAR

Abstract

The application discloses a facility deformation analysis method based on PSInSAR and SqueeSAR, and relates to the field of time-series remote sensing image processing; specifically: first, for the to-be-measured region, SAR image is collected, and after cutting and registration, an interferogram is generated; then, the orbit state vector of the satellite is interpolated and fitted; the DEM data of the target S in the WGS84 coordinate system is projected into the SAR coordinate system to obtain a reference position and calculate the absolute value of the Doppler frequency of each pixel of the auxiliary image, so that the satellite position S i is the reference position P' of the target S in each SAR auxiliary image; the star-ground simulation phase is removed and a differential interference phase is generated; PSInSAR is performed, and PS point selection is completed through four iterations; on the basis of PSInSAR, SqueeSAR is performed, Newton iteration method is used for geographic coding of the target S, and the actual position is obtained; corresponding to the deformation rate, the complete ground deformation range and the time-series deformation graph of the corresponding monitoring point are obtained. The application effectively improves the monitoring and subsequent evaluation of facility deformation.

Description

Technical Field

[0001] This invention relates to the field of time-series remote sensing image processing, specifically to a facility deformation analysis method based on PSInSAR (Permanent ScattererInSAR) and SqueeSAR (time-series InSAR). Background Technology

[0002] Currently, Interferometry Synthetic Aperture Radar (InSAR) technology, as a novel remote sensing measurement technique, is widely used in fields such as surface deformation measurement, natural disaster monitoring, earthquake monitoring, and mine subsidence monitoring. Its successful application marks the transformation of radar remote sensing from a qualitative scientific research field to a quantitative tool. InSAR technology and its extension, D-InSAR (Differential Interferometry Synthetic Aperture Radar), are beneficial for revealing the mechanisms of surface deformation and earthquakes, providing a novel and efficient method for dynamic research in geology and geophysics. D-InSAR technology has unparalleled advantages over other technologies, such as enabling high-precision monitoring of surface changes over wide areas, and its data processing procedures are relatively simple.

[0003] While D-InSAR technology has unparalleled advantages over other technologies, changes in the scattering properties of ground objects in SAR images between two observations, excessively long spatial baselines between two images, or discontinuous atmospheric changes over the observed surface can reduce or even completely decohere the interferograms acquired by D-InSAR. This decoherence manifests as noise errors or even disordered interferometric fringes in the interferograms. These errors severely affect the accuracy of D-InSAR measurement results and limit the application of this technology. Summary of the Invention

[0004] To overcome the problem that discontinuous atmospheric changes observed above the Earth's surface in D-InSAR technology can lead to reduced or even complete loss of coherence in interferograms, severely affecting the accuracy of measurement results, this invention proposes a facility deformation analysis method based on PSInSAR and SqueeSAR. It applies a new InSAR method using permanent coherent points, namely PSInSAR, to generate an average deformation rate map by identifying PS points with strong reflection signals and good coherence in radar images, thus calculating ground deformation with millimeter-level precision. Simultaneously, the DS-InSAR (Distributed Satellite System) method is used to screen out homogeneous points. Based on the statistical characteristics of these homogeneous points, the phase of distributed targets is optimized to reduce the impact of noise on distributed targets. The combined use of DS points and time-series InSAR methods for deformation calculation increases the point target density in non-urban areas, allowing for the acquisition of the complete surface deformation range.

[0005] The facility deformation analysis method based on PSInSAR and SqueeSAR has the following specific steps:

[0006] Step 1: For the administrative region to be measured, use SAR to collect images of the location at different times to form a multi-view image set;

[0007] Step 2: Crop and register the image set, and use the coherence coefficient of the interferogram to select the main image and auxiliary image that meet the interferometric quality requirements from the SAR image clipping;

[0008] The specific process is as follows:

[0009] First, the SAR images of each frame in the image set at the current time point are divided into blocks of 1024*1024 size with 128 overlapping pixels at the boundary, so as to obtain the number of rows and columns of the image.

[0010] Then, the 1024*1024 image sub-blocks are divided into non-overlapping 256*256 sub-blocks, and interpolation is performed by a factor of 256 or 512 to form coarse registration, thus obtaining the row and column registration offset of any sub-block i.

[0011] Next, based on the row and column registration offset of each image sub-block, the sub-block corresponding to the main image sub-block is found. Corresponding set of all auxiliary image sub-blocks And based on the interpolation kernel function of size kernelsize, for all auxiliary image subsets By performing boundary extension, a new set of extended auxiliary image blocks is obtained.

[0012] Extended auxiliary image sub-blocks Parallel resampling is performed to output finely registered SAR auxiliary image sub-blocks of size 1024*1024, and geometric averaging is performed on the finely registered SAR auxiliary image sub-blocks in the overlapping areas to generate corresponding interferograms.

[0013] Finally, the SAR coherence coefficient in the interferogram is calculated;

[0014] The calculation formula is:

[0015]

[0016] Where ρ (i,j) This represents the coherence coefficient value of the element in the i-th row and j-th column; W1 and W2 are the window calculation sizes; * indicates the conjugate of the complex data; This represents the set of main image sub-blocks containing the element in the i-th row and j-th column; This represents the set of extended auxiliary image sub-blocks containing the elements in the i-th row and j-th column.

[0017] Step 3: Interpolate and fit the 12 orbital state vectors of the SAR satellite's position, velocity, and acceleration to obtain the satellite state vector at the imaging time of each resolution cell of the SAR.

[0018] Step 4: Select the target S in the area to be measured from the main image that meets the interferometric quality requirements. Project the DEM data in the WGS84 coordinate system onto the SAR coordinate system to obtain the reference DEM location.

[0019] Based on the SAR range equation and the Doppler equation:

[0020]

[0021]

[0022] Where |·| represents Euclidean distance, and <·,·> represents vector dot product. The position of the DEM of target S in the WGS84 coordinate system projected onto the SAR coordinate system; The current velocity of the satellite relative to the ground. The current spatial position of the satellite is given by λ; r is the distance between the target S and the satellite at the imaging time; λ refers to the wavelength of the SAR; and f is the distance between the target S and the satellite at the imaging time. dc The frequency of the Doppler center.

[0023] Then, using the imaging speed formula:

[0024]

[0025]

[0026] (row,col) represents the position. The corresponding row and column position of target S in the SAR image; t is the time when the satellite images target S; t min r is the first-phase imaging time of the SAR image; min ρ represents the distance between the satellite and the prime ministerial target. r f represents the range resolution of the SAR image. prf The pulse repetition frequency, This represents the floor operator.

[0027] Step 5: For each auxiliary image that meets the interference quality requirements, when the absolute value of the Doppler frequency of each pixel in the auxiliary image is the smallest, the corresponding satellite position S is... i The reference DEM location P' of target S in each SAR auxiliary image;

[0028] Step 6: Calculate the phase of the simulated interferometry between satellite and ground. Further calculation of the differential interference phase of the main image and the auxiliary image

[0029] Star-ground simulated interferometry phase The formula is as follows:

[0030]

[0031] The differential interference phase formula is as follows:

[0032]

[0033] S m To achieve the required main image quality for interference, S s* Auxiliary images that meet interference quality requirements;

[0034] Step 7: Utilizing differential interference phase PSinSAR phase unwrapping is performed on the pixel resolution units of the main image and the auxiliary image, and permanent scattering points (PS points) in each image are found after four iterations.

[0035] Phase unwrapping utilizes a multi-component temporal coherence model to obtain the temporal coherence coefficients of adjacent PS candidate points. The calculation formula is:

[0036] in, For any two adjacent PS candidate points x1 and x2, the phase difference operator is... Let Δφ represent the differential interference phase of the x-th PS candidate point in the i-th interferogram, M be the total number of interferograms, and Δφ be the differential interference phase. H,x,i and These represent the simulated residual terrain and the interferometric phase of the linear deformation of the x-th PS candidate point in the i-th interferogram, respectively.

[0037] A nonlinear blind source separation method based on spatiotemporal filtering is used to separate components such as orbital error, atmospheric delay phase, residual terrain phase, and deformation phase. Homogeneous points are removed, and the coherence coefficient, three-dimensional phase unwrapping value, and interference phase are iterated. The output of each iteration includes: temporal coherence coefficient, DEM after residual elevation correction, average deformation rate map, and time-series deformation map of monitoring points.

[0038] Three-dimensional phase unwrapping includes: construction of time-baseline plane constrained triangulation for one-dimensional time phase unwrapping, and construction of distance-azimuth plane constrained triangulation for two-dimensional spatial phase unwrapping;

[0039] The component separation adopts a nonlinear blind source separation method based on spatiotemporal filtering, and the obtained coherence coefficients are used for PS selection and uncertainty assessment.

[0040] Step 8: Select Ds points based on PSInSAR and perform DespecKS and fusion of DS and PS based on PSInSAR.

[0041] DespecKS includes: using FaSHP to quickly identify DSs targets, calculating the sample coherence matrix based on the identified DSs locations, and performing DSs phase filtering based on PTA.

[0042] Step 9: Utilize the reference DEM location in the main image Combined with the reference DEM location P' in each auxiliary image and the satellite state vector, the target S is geocoded using the Newton-Raphson iteration method to obtain the actual geographical location.

[0043] The calculation formula is:

[0044]

[0045] (x p ,y p ,z p R is the actual geographical coordinate of the target S. a R is the radius of the major semi-axis, Δh is the elevation difference, and R b Let x be the radius of the minor semi-axis, (x) m ,y m ,z m ) represents the coordinates of the geographic coding standard point m, |r mp | represents the absolute distance between the target S and the standard point m, (v mx ,v my ,v mz Let v be the velocity in each direction at the standard point m. m -v p Let r be the velocity difference between the target S and the reference point m. m -r p Let S be the radius difference between the target point S and the standard point m;

[0046] Step 10: Based on the actual geographical location of target S obtained from the geocoding, correspond to its deformation rate to obtain the complete surface deformation range; at the same time, obtain the average deformation rate map and the time series deformation map of the corresponding monitoring points for subsequent analysis and result conclusion.

[0047] The advantages of this invention are:

[0048] This invention is based on a facility deformation analysis method using PS-InSAR and SqueeSAR. It establishes a complete system for analyzing the results of PS-InSAR and SqueeSAR, realizing the analysis of facility deformation, effectively improving the monitoring and subsequent evaluation of facility deformation, and enhancing the application capabilities of PS-InSAR and SqueeSAR in geological disasters. Attached Figure Description

[0049] Figure 1 This is a flowchart of the facility deformation analysis method based on PSInSAR and SqueeSAR of the present invention;

[0050] Figure 2 This is the user interface of the Beijing Capital International Airport example used in this invention;

[0051] Figure 3 This is a schematic diagram illustrating the interpolation and fitting of 12 orbital state vectors of a satellite according to the present invention;

[0052] Figure 4 This is a detailed flowchart of the PSInsar invention;

[0053] Figure 5 This is a detailed flowchart of the SqueeSAR invention.

[0054] Figure 6 This is a comparison diagram of surface subsidence within a time interval, representing an example of the present invention.

[0055] Figure 7 The PS point selection results and key step calculation time for the four iterations of this invention;

[0056] Figure 8 The time-baseline constrained triangulation used for PS point selection in this invention;

[0057] Figure 9 This is the statistical result of the high coherence edge of the distance-azimuth plane-constrained triangulation network during the four iterations of this invention;

[0058] Figure 10 These are the results of coherence coefficient calculations during the four iterations of this invention.

[0059] Figure 11 The corrected DEM, average deformation rate, and time-series deformation diagram of a randomly selected PS point are shown in this invention.

[0060] Figure 12 This is a comparison chart of the existing SBAS deformation results and DEM products with the results of the algorithm of this invention;

[0061] Figure 13 The uncertainty of the existing SBAS deformation results and the PSInSAR results of the DEM product compared with the algorithm of this invention is given. Detailed Implementation

[0062] The present invention will be further described below with reference to the accompanying drawings and examples.

[0063] The core idea of ​​PS-InSAR technology is to statistically analyze the backscattering characteristics of ground objects in long-term SAR interferograms to find stable target models, i.e., target ground objects with high coherence in SAR images. This strategy will cause some ground objects with poor coherence in some interferograms to be completely discarded, thus limiting its application.

[0064] In urban surface monitoring with numerous man-made structures, the backscattering coefficient of each pixel in a SAR image is mainly provided by scatterers with relatively high backscattering coefficients, such as buildings and concrete pavements. These features are frequently detected in PS-InSAR and are referred to as PS ground point targets. However, in natural areas with few or no man-made structures, the surface mainly consists of bare land, rocks, or sparse vegetation. Their backscattering coefficients are basically the same. The observation value of each pixel in the SAR image is an integral within the surface pixel range weighted by the scattering coefficients of these features. These features are called distributed targets (DTs), and they have low coherence in PS-InSAR and cannot be effectively detected. In summary, PS-InSAR technology is suitable for urban surface monitoring with many PS point targets, but not for DT target monitoring in natural areas.

[0065] To address the shortcomings of D-InSAR technology, this invention utilizes PS-InSAR and SqueeSAR for facility deformation analysis, proposing a facility deformation analysis method based on PS-InSAR and SqueeSAR, such as... Figure 1 As shown, the specific steps are as follows:

[0066] Step 1: For the administrative region to be measured, use SAR to collect images of the location at different times to form a multi-view image set;

[0067] Step 2: Crop and register the image set, and use the coherence coefficient of the interferogram to select the main image and auxiliary image that meet the interferometric quality requirements from the SAR image clipping;

[0068] The specific process is as follows:

[0069] First, the SAR images of each frame in the image set at the current time point are divided into blocks of 1024*1024 size with 128 overlapping pixels at the boundary, so as to obtain the number of rows and columns of the image.

[0070] Then, the 1024*1024 image sub-blocks are divided into non-overlapping 256*256 sub-blocks, resulting in 4*4 grandchild image blocks. Interpolation is then performed at a factor of 256 or 512. Based on the maximum correlation coefficient criterion, coarse registration is performed in parallel on each image sub-block to obtain the row and column registration offset of any sub-block i.

[0071] Next, based on the row and column registration offset of each image sub-block, the sub-block corresponding to the main image sub-block is found. Corresponding set of all auxiliary image sub-blocks And based on the interpolation kernel function of size kernelsize in the fine registration, for all auxiliary image subsets By performing boundary extension, a new set of extended auxiliary image blocks is obtained.

[0072] Extended auxiliary image sub-blocks Parallel resampling is performed to output finely registered SAR auxiliary image sub-blocks of size 1024*1024, and geometric averaging is performed on the finely registered SAR auxiliary image sub-blocks in the overlapping areas to generate corresponding interferograms.

[0073] Finally, the SAR coherence coefficient in the interferogram is calculated;

[0074] The calculation formula is:

[0075]

[0076] Where ρ (i,j) This represents the coherence coefficient value of the element in the i-th row and j-th column; W1 and W2 are the window calculation sizes; * indicates the conjugate of the complex data; This represents the set of main image sub-blocks containing the element in the i-th row and j-th column; This represents the set of extended auxiliary image sub-blocks containing the elements in the i-th row and j-th column.

[0077] Step 3: Interpolate and fit the 12 orbital state vectors of the SAR satellite's position, velocity, and acceleration to obtain the satellite state vector at the imaging time of each resolution cell of the SAR.

[0078] This step provides high-precision satellite data for subsequent geocoding.

[0079] Step 4: Select the target S in the area to be measured from the main image that meets the interferometric quality requirements. Project the DEM data in the WGS84 coordinate system onto the SAR coordinate system to obtain the reference DEM location.

[0080] Based on the SAR range equation and the Doppler equation:

[0081]

[0082]

[0083] Where |·| represents Euclidean distance, and <·,·> represents vector dot product. The position of the DEM of target S in the WGS84 coordinate system projected onto the SAR coordinate system; The current velocity of the satellite relative to the ground. The current spatial position of the satellite is given by λ; r is the distance between the target S and the satellite at the imaging time; λ refers to the wavelength of the SAR; and f is the distance between the target S and the satellite at the imaging time. dc The frequency of the Doppler center.

[0084] Then, using the imaging speed formula:

[0085]

[0086]

[0087] (row,col) represents the position. The corresponding row and column position of target S in the SAR image; t is the time when the satellite images target S; t min r is the first-phase imaging time of the SAR image; min ρ represents the distance between the satellite and the prime ministerial target. r f represents the range resolution of the SAR image. prf The pulse repetition frequency, This represents the floor operator.

[0088] Determine the integer value as arbitrary The row and column positions (row, col) of a target in a SAR image are denoted as... The position in the WGS84 coordinate system is: the reference DEM position projected onto the SAR image at row and column position (row,col), i.e. It is the position of the DEM projection of target S in the WGS84 coordinate system onto the SAR coordinate system;

[0089] Step 5: For each auxiliary image that meets the interference quality requirements, when the absolute value of the Doppler frequency of each pixel in the auxiliary image is the smallest, the corresponding satellite position S is... i The reference DEM location P' of target S in each SAR auxiliary image;

[0090] P' refers to the DEM position of the target location determined according to the Doppler equation in the WGS84 coordinate system.

[0091] Step 6: Calculate the phase of the simulated interferometry between satellite and ground. Further calculation of the differential interference phase of the main image and the auxiliary image

[0092] Star-ground simulated interferometry phase The formula is as follows:

[0093]

[0094] The differential interference phase formula is as follows:

[0095]

[0096] S m To achieve the required main image quality for interference, S s* Auxiliary images that meet interference quality requirements;

[0097] Step 7: Utilizing differential interference phase PSinSAR phase unwrapping is performed on the pixel resolution units of the main image and the auxiliary image, and permanent scattering points (PS points) in each image are found after four iterations.

[0098] Points on the polarization field (PS) are generally considered the dominant scatterer within a pixel-resolution unit, and their extraction is a crucial prerequisite for obtaining coseismic deformation field results. The PS can be viewed as a natural corner reflector; therefore, radar signals exhibit long-term stability on this scatterer, almost unaffected by surrounding weak scatterers, and largely unaffected by geometric decoherence. For surface deformation measurements, the stability of the PS is primarily reflected in its phase stability.

[0099] Phase unwrapping utilizes a multi-component temporal coherence model to obtain the temporal coherence coefficients of adjacent PS candidate points. The calculation formula is:

[0100] in, For any two adjacent PS candidate points x1 and x2, the phase difference operator is... Let Δφ represent the differential interference phase of the x-th PS candidate point in the i-th interferogram, M be the total number of interferograms, and Δφ be the differential interference phase. H,x,i and These represent the simulated residual terrain and the interferometric phase of the linear deformation of the x-th PS candidate point in the i-th interferogram, respectively.

[0101] Its difference form can be expressed as:

[0102]

[0103]

[0104] Where b i With t i Let θ and λ represent the vertical baseline and observation time of the i-th image, respectively, and let θ and λ represent the radar incident angle and wavelength, respectively. t0 represents the initial deviation constant of the periodic motion. This optimization problem can be solved quickly using a grid search method under CPU parallelism and GPU computing strategies.

[0105] A nonlinear blind source separation method based on spatiotemporal filtering is used to separate components such as orbital error, atmospheric delay phase, residual terrain phase, and deformation phase. Homogeneous points are removed, and the coherence coefficient, three-dimensional phase unwrapping value, and interference phase are iterated. The output of each iteration includes: temporal coherence coefficient, DEM after residual elevation correction, average deformation rate map, and time-series deformation map of monitoring points.

[0106] Three-dimensional phase unwrapping based on time coherence coefficient mainly includes the following steps:

[0107] 1. Construct a spatial triangulation network based on the PS candidate point set;

[0108] 2. For the phase difference of two PS candidate points connected by the triangular mesh, maximize the multi-component time coherence model to obtain the time coherence coefficient of adjacent PS candidate points, as well as the residual phase gradient and the average deformation rate gradient.

[0109] 3. Construct a constrained triangulation in the baseline-time two-dimensional plane and remove triangles with longer baseline and time intervals;

[0110] 4. Based on the extended minimum cost flow method, phase unwrapping is performed on the phase gradient in the time dimension. The time unwrapped value is obtained by estimating the main image prior and integer linear programming.

[0111] 5. Remove the low-coherence triangulation edges from step 2, and construct a highly coherent constrained maximum connected triangulation in the azimuth-distance two-dimensional plane;

[0112] 6. Based on the phase gradient time unwrapping value in step 4, use the sparse minimum cost flow method to unwrap the residual wrapping phase gradient multiple in the spatial dimension and solve the optimization model;

[0113] 7. Integrate the unwrapped phase gradient obtained in step 6 along the optimal path obtained in step 5 to obtain the unwrapped value of the entangled interference phase. Based on the estimated value obtained from the expression of each side of the triangular mesh unfolded in the azimuth-range plane, use the time coherence coefficient in step 2 as the weight of the equation, and solve for the unwrapped value of the entangled interference phase by weighted least squares.

[0114] Step 8: Select point Ds based on PSInSAR and perform SqueeSAR;

[0115] SqueeSAR includes: DespecKS, a fusion of DS and PS based on PSInSAR;

[0116] Based on psinsar, the Ds point was selected, which increased the amount of reliable data compared to psinsar.

[0117] DespecKS includes: using FaSHP to quickly identify DSs targets, calculating the sample coherence matrix based on the identified DSs locations, and performing DSs phase filtering based on PTA.

[0118] The fusion method is mainly based on PSInSAR processing.

[0119] Step 9: Utilize the reference DEM location in the main image The reference DEM location P' in each auxiliary image is combined with the satellite state vector, and the Newton-Raphson iteration method is used to geocode the initial target S of the image to obtain the actual geographical location.

[0120] Using the SAR auxiliary image imaging position P' as the initial value, the actual geographical location's latitude, longitude, elevation, and average deformation rate are obtained. Based on the Earth ellipsoid model, Doppler equations, and distance equations:

[0121]

[0122] (x p ,y p ,z p R is the actual geographical coordinate of the target S. a R is the radius of the major semi-axis, Δh is the elevation difference, and R b Let x be the radius of the minor semi-axis, (x) m ,y m ,z m ) represents the coordinates of the geographic coding standard point m, |r mp | represents the absolute distance between the target S and the standard point m, (v mx ,v my ,v mz Let v be the velocity in each direction at the standard point m. m -v p Let r be the velocity difference between the target S and the reference point m. m -r p Let S be the radius difference between the target point S and the standard point m;

[0123] Step 10: Based on the actual geographical location of target S obtained from the geocoding, correspond to its deformation rate to obtain the complete surface deformation range; at the same time, obtain the average deformation rate map and the time series deformation map of the corresponding monitoring points for subsequent analysis and result conclusion.

[0124] In this invention, the identification of homogeneous points refers to determining whether two pixels within a certain spatial range in a SAR image dataset belong to the same ground feature. The main principle is to extract the intensity information of the same pixel in the time dimension, measure the similarity between two samples, and thus determine whether they are homogeneous points; homogeneous points are then selected for use in subsequent iterations to select ps points.

[0125] Example:

[0126] Settlement analysis was conducted using Beijing Capital International Airport in Shunyi District as an example. The main operating interface is as follows: Figure 2 As shown:

[0127] The dataset consists of 31 TerraSAR images, covering the period from 2012 to 2016. The image size is 4608*4608 pixels, and the azimuth and distance resolutions are 1.9 meters and 0.9 meters, respectively. See the example data list below:

[0128]

[0129]

[0130] The specific process is as follows:

[0131] Step 1: Crop and register the image, and generate an interferogram;

[0132] First, the SAR image is cropped into blocks of 1024×1024 size with 128 overlapping pixels, to obtain the number of rows and columns of the image. Then, the coarsely registered image blocks of 1024×1024 size are divided into non-overlapping 256×256 blocks. Next, these 4×4 grandchild image blocks are interpolated by a factor of 256 or 512.

[0133] The following phase-independent coherent estimation method is used for parallel computation. Based on the row and column registration offsets of all pixels in any fitted image sub-block, the main steps include coarse and fine registration of SAR images. Generally, in the registration process between the main image and the auxiliary image, coarse registration is performed first, followed by fine registration. Fine registration generally requires an accuracy of 1 / 8 pixel.

[0134] Parallel resampling of the extended auxiliary image sub-blocks is performed using an interpolation kernel function of size kernelsize; finely registered SAR auxiliary image sub-blocks of the same size as the image sub-blocks are output, and geometric averaging is performed on the auxiliary image sub-blocks in the overlapping regions and corresponding interferograms are generated.

[0135] Step 2: Generate coherence coefficients and obtain appropriate orbital state plots by interpolating and fitting the 12 orbital state vectors (position, velocity, and acceleration) of the satellite.

[0136] like Figure 3 As shown, high-precision satellite data is provided for subsequent geocoding based on the geometric relationships of SAR observations.

[0137] Step 3: Project the external DEM to the SAR coordinate system. Project the DEM in the WGS84 coordinate system within the observation area to the geocentric inertial SAR coordinate system.

[0138] This step is mainly to obtain the DEM of the SAR coordinate system, which will facilitate subsequent satellite-to-ground simulation interferometric phase calculation.

[0139] Step 4: Satellite-to-Ground Simulation Phase Removal and Differential Interferometric Phase Generation. The reference DEM in the determined SAR coordinate system is used to determine, according to the formula, the satellite position S corresponding to the minimum absolute value of the Doppler frequency of each pixel in the auxiliary image. i As the imaging location for SAR auxiliary images.

[0140] Step 5: Perform PSinSAR; PS point selection can be completed in 4 iterations.

[0141] like Figure 4 The diagram shown is a flowchart of PSinSAR. The key intermediate steps of each iteration are: three-dimensional phase unwrapping, separation of components such as orbital error, atmospheric delay phase, residual terrain phase, and deformation phase.

[0142] The key outputs of each iteration mainly include: time coherence coefficient, DEM after residual elevation correction, average deformation rate map, and time-series deformation map of monitoring points.

[0143] The tomography function outputs the corresponding TomSAR digital elevation model geocoded image TIFF plot and the TomSAR average deformation rate geocoded image TIFF plot.

[0144] Step 6: Perform SqueeSAR based on PSInSAR;

[0145] The flowchart of SqueeSAR is as follows: Figure 5 As shown.

[0146] Step 7: Use Newton's iteration method to geocode the actual location of target S to obtain the actual latitude, longitude, altitude and average deformation rate.

[0147] Step 8: Individual facility analysis function. Click on the corresponding point to view the settlement situation within the corresponding time interval.

[0148] The specific generated results are as follows: Figure 6 As shown.

[0149] In this invention, the PS point selection results and the calculation time of key steps are obtained through four iterations, such as... Figure 7 As shown, the total time required for PSInSAR to run more than 30,000 PS points is less than 9 hours. Figure 7 (a) The left side represents the number of PS candidate points before each iteration, and the right side represents the number of PS candidate points retained after each iteration; Figure 7 (b) The left side represents the time taken to calculate the time coherence coefficient for each iteration, and the right side represents the time required for time and space phase unwrapping.

[0150] A spatiotemporally constrained triangulation sparse 3D phase unwrapping method based on a time coherence model is adopted. Key intermediate results include: construction of a time-baseline plane constrained triangulation for one-dimensional temporal phase unwrapping, and construction of a distance-azimuth plane constrained triangulation for two-dimensional spatial phase unwrapping. The time-baseline constrained triangulation is shown below. Figure 8 As shown, the distance-azimuth plane constrained triangulation of the four-iteration update process is as follows: Figure 9 As shown in Table 1, the number of correct unwrapped edges and the error rate were selected as quantitative indicators for the three-dimensional unwrapping process in successive iterations. With increasing iteration count, the unwrapping error rate decreased; with increasing temporal coherence coefficient, the number of correct unwrapped edges increased; and the unwrapping accuracy in the fourth iteration reached as high as 99.163%.

[0151] Evaluation of three-dimensional phase unwrapping results

[0152]

[0153] A nonlinear blind source separation method based on spatiotemporal filtering is employed to separate components such as orbital error, atmospheric delay phase, residual terrain phase, and deformation phase, thereby obtaining the coherence coefficient for PS selection and uncertainty assessment. The calculation results of the coherence coefficient during the four-iteration update process are shown below. Figure 10 As shown, with increasing iterations, low-coherence pixels are successively eliminated, leaving only high-coherence PS pixels. The final corrected DEM, average deformation rate, and temporal deformation of a randomly selected PS point are shown in the figure. Figure 11 As shown, Figure 12 A comparison was made between the SBAS deformation results published in the multiview journal and the 90-meter resolution DEM product of TanDEM, and it was found that the distribution of the obtained results was similar to... Figure 10 The results are consistent, but differ in that this invention obtains results under a single view, resulting in higher resolution and richer detail. According to unbiased estimation theory, the uncertainty of the corrected DEM and the mean deformation rate results is as follows: Figure 13 As shown, the elevation uncertainty of point PS is better than 1.8 meters, and the average deformation rate uncertainty is better than 0.6 mm / yr.

[0154] This invention uses DInSAR technology to interferoscopically remove geographic terrain models and obtain qualitative information on ground deformation. However, it cannot achieve ground displacement measurement accuracy above the centimeter level, is severely affected by atmospheric interference and noise, and cannot obtain time-series data. Simply using PSInSAR technology to find stable ground targets, estimate and remove the influence of atmospheric factors, and obtain time-series ground deformation with millimeter-level ground displacement measurement accuracy is also problematic. However, stable ground target points are unevenly distributed regionally and are generally man-made targets.

Claims

1. A method for analyzing facility deformation based on PSInSAR and SqueeSAR, characterized in that, The specific steps are as follows: Step one, for the administrative region to be tested, SAR is used to collect images of the place at different times to form a multi-scene image set; Step two, the image set is cropped and image registration is performed, and the main image and auxiliary image meeting the interference quality requirements are selected from the SAR image clip using the coherence coefficient of the interferogram; Step three, 12 orbit state vectors of the SAR satellite, including position, speed and acceleration, are interpolated and fitted to obtain the satellite state vector at the imaging moment of each resolution unit of the SAR; Step four, select the target S of the region to be measured from the main image meeting the interference quality requirements, project the DEM data in the WGS84 coordinate system to the SAR coordinate system to obtain the reference DEM position Step five, for each sub-image meeting the interference quality requirements, when the absolute value of the Doppler frequency of each pixel in the sub-image is the smallest, the corresponding satellite position S i is the reference DEM position P' of the target S in each SAR sub-image; Step six, calculating the star-to-ground simulated interferometric phase Further calculating the difference interferometric phase of the primary and secondary images Step seven, using differential interferometric phase The phase unwrapping of PSinSAR is performed on the pixel resolution cells of the primary and secondary images, and the permanent scatterers PS points in each image are found after four iterations; In the phase unwrapping, the maximum multi-component time coherence model is used to obtain the time coherence coefficients of adjacent PS candidate points The calculation formula is: wherein, is the phase difference operator between two adjacent PS candidates x1 and x2, φx,i represents the differential interferometric phase of the xth PS candidate in the ith interferogram, M is the total number of interferograms, Δφ H,x,i and φx,i represents the simulated residual topography and the linear deformation interferometric phase of the xth PS candidate in the ith interferogram, respectively. Step eight, based on PSInSAR, Ds points are selected, DespecKS is performed, and DS and PS fusion based on PSInSAR is performed; DespecKS includes: selecting FaSHP to quickly identify DSs targets, calculating a sample coherence matrix according to the positions of the identified DSs, and performing DSs phase filtering based on PTA; Step nine, use the reference DEM position in the primary image and the reference DEM position P' in each secondary image, in combination with the satellite state vector, to geocode the target S using Newton's iteration method to obtain the actual geographic position; Step ten, according to the actual geographic location of the target S obtained through geographic coding, the corresponding deformation rate is obtained to obtain the complete surface deformation range; at the same time, the average deformation rate map and the time series deformation map corresponding to the monitoring point are obtained, which are used for subsequent analysis and result conclusion.

2. The PSInSAR and SqueeSAR based facility deformation analysis method of claim 1, wherein, The step two is specifically: Firstly, the SAR images in the image set at the current time point are divided into blocks according to 1024*1024 size and 128 pixel boundary overlap, respectively, to obtain the row and column block numbers of the images; Then, the 1024*1024 size image sub-block is not overlapped 256*256 divided, 256 or 512 times interpolation is carried out, coarse registration is formed, and the row and column registration offset of any sub-block i is obtained Then, according to the row and column registration offset of each image sub-block, find all the corresponding auxiliary image sub-block sets of the main image sub-block According to the size of the interpolation kernel function, kernelsize, all auxiliary image sub-block sets of the main image sub-block are extended to obtain a new extended auxiliary image sub-block set Extended auxiliary image sub-block Parallel resampling is performed, and a fine registration SAR auxiliary image sub-block with a size of 1024*1024 is output, and a geometric mean of the fine registration SAR auxiliary image sub-blocks in the overlapping area is calculated to generate an interference image. Finally, the SAR coherence coefficient in the interferogram is calculated; The calculation formula is: where ρ (i,j) represents the coherence coefficient value of the element in the ith row and jth column; W1, W2 are window calculation sizes; * represents the conjugate of complex data; represents the main image sub-block set of the element in the ith row and jth column; represents the extended auxiliary image sub-block set of the element in the ith row and jth column.

3. The PSInSAR and SqueeSAR based facility deformation analysis method of claim 1, wherein, In the step four, according to the SAR distance equation and the Doppler equation: Where |·| represents the Euclidean distance, <·,·> represents the vector inner product, P is the position of the DEM of the target S in the WGS84 coordinate system projected to the SAR coordinate system; is the velocity of the satellite relative to the ground at the current time, is the spatial position of the satellite at the current time; r is the distance between the target S and the satellite at the imaging time, λ refers to the wavelength of the SAR, f dc is the frequency of the Doppler center; Then, the imaging fast and slow time formula is used: (row, col) is the position corresponding target S in the SAR image; t is the time when the satellite images the target S; t min is the first imaging time of the SAR image; r min is the distance between the satellite and the first imaging target; p r is the range resolution of the SAR image; f prf is the pulse repetition frequency, denotes the floor operator.

4. The PSInSAR and SqueeSAR based facility deformation analysis method of claim 1, wherein, In step six, the satellite-ground simulated interference phase The formula is as follows: The differential interference phase formula is as follows: S m S s* S 5. The PSInSAR and SqueeSAR based facility deformation analysis method of claim 1, wherein, In the step seven, a nonlinear blind source separation method based on space-time filtering is used to separate the orbit error, atmospheric delay phase, residual terrain phase, deformation phase components, remove homogeneous points, and iterate the coherence coefficient, three-dimensional phase unwrapping value, and interference phase; the output of each iteration result includes: time coherence coefficient, residual elevation corrected DEM, average deformation rate map, and time series deformation map of monitoring points; Three-dimensional phase unwrapping includes: time-baseline plane constraint triangular network construction for one-dimensional phase unwrapping in time, and distance-azimuth plane constraint triangular network construction for two-dimensional phase unwrapping in space; The component separation adopts a nonlinear blind source separation method based on space-time filtering to obtain the coherence coefficient for PS selection and uncertainty evaluation.

6. The PSInSAR and SqueeSAR based facility deformation analysis method of claim 1, wherein, In the step nine, the calculation formula is: (x p ,y p ,z p are the actual geographic position coordinates of the target S, R a is the radius of the major semiaxis, Δh is the height difference, R b is the radius of the minor semiaxis, (x m ,y m ,z m ) are the coordinates of the geocoding standard point m, r mp is the absolute distance of the target S from the standard point m, (v mx ,v my ,v mz ) are the respective directional velocities of the standard point m, v m -v p is the velocity difference of the target S from the standard point m, r m -r p is the radius difference of the target S from the standard point m.

Images