Radar distributed scatterer interferometry method and device for high and steep regions

Abstract

The invention discloses a radar distributed scatterer interference measurement method and device for a high and steep area. The method comprises the following steps: selecting a main image and a non-main image; calculating differential interference data; screening PS candidate pixels of which the amplitude deviation indexes are smaller than an amplitude deviation threshold value; high-coherence DS candidate pixels are screened out; constructing a distance-constrained Delaunay triangulation network for fusing the PS candidate pixels and the high-coherence DS candidate pixels; resolving a differential residual elevation, a differential deformation rate, a set coherence coefficient of a connected edge and a spatial differential residual phase; calculating residual terrain, linear deformation, deformation rate and residual phase after unwrapping; calculating an atmospheric orbit phase, nonlinear deformation and a noise phase; and obtaining time sequence deformation based on the linear deformation and the nonlinear deformation. According to the method, the number of monitoring points of the orbit descending InSAR deformation monitoring result is large, the coverage range of the monitoring result is widest, and MAE and RMSE are small.

Description

Technical Field

[0001] The invention belongs to the technical field of interferometric synthetic aperture radar data processing, and in particular relates to a radar distributed scatterer interferometric measurement method and equipment for high-steep areas. Background Art

[0002] DS-InSAR is a new geodetic technology that has developed rapidly in recent years. It fully combines the DS target that maintains high temporal coherence characteristics after spatial adaptive filtering and the PS target that maintains high temporal coherence in single-view conditions to perform interferometric point target analysis, which can effectively overcome the shortcomings of PSInSAR in natural scene monitoring, but inherits the advantages of high-precision and high-resolution deformation monitoring of PSI. Since it was proposed in 2011, DS-InSAR has become one of the key development technologies of time-series InSAR, providing an important way for precise deformation monitoring of complex scenes such as volcanoes, mining areas, landslides, and urban areas.

[0003] Homogeneous point identification and phase optimization are the key technologies of DS-InSAR. The first generation of DS-InSAR technology SqueeSAR TM The method first uses the amplitude-based two-sample Kolmogorov–Smirnov (KS) hypothesis test method to identify homogeneous points; then estimates the covariance matrix based on these homogeneous pixels; finally, the phase triangulation method based on maximum likelihood estimation (ML) is used to filter the extracted DS target phase and separate it from the main scattering mechanism. Although the spatial density and coverage of monitoring points are increased, there are also shortcomings such as insufficient robustness of small samples and low computational efficiency, which limit the promotion and application of DS-InSAR technology. Summary of the invention

[0004] The purpose of the present invention is to address the problems existing in the prior art and provide a radar distributed scatterer interferometry method and device for high and steep areas, so as to improve the monitoring coverage, density and reliability of high mountain valley areas.

[0005] The above-mentioned purpose of the present invention is achieved by the following technical means:

[0006] A radar distributed scatterer interferometry method for high and steep areas comprises the following steps:

[0007] Step 1: Select the main image and non-main image from the multi-temporal SAR image data;

[0008] Step 2: Perform conjugate operation between the main image and the non-main image to remove the terrain phase introduced by the digital elevation model and generate differential interferometry data;

[0009] Step 3: Calculate the amplitude deviation index of the differential interferometer data, and select the PS candidate pixels whose amplitude deviation index is less than the amplitude deviation threshold;

[0010] Step 4: Screen out high coherence DS candidate pixels based on differential interferometry data;

[0011] Step 5: Construct a distance-constrained Delaunay triangulation network that fuses PS candidate pixels with high-coherence DS candidate pixels;

[0012] Step 6: Maximize the set coherence model between two points connected by the distance-constrained Delaunay triangulation network edge, and obtain the differential residual elevation, differential deformation rate, set coherence coefficient of the connected edge, and spatial differential residual phase;

[0013] Step 7: Perform three-dimensional phase unwrapping on the spatial differential residual phase obtained in step 6 to obtain residual topography, linear deformation, deformation rate and unwrapped residual phase;

[0014] Step 8: Bandpass filter the unwrapped residual phase to obtain the atmospheric orbit phase, nonlinear deformation and noise phase;

[0015] Step 9: If the temporal coherence coefficients of the PS candidate pixels and the high coherence DS candidate pixels are both greater than the set coherence threshold, or the maximum number of iterative cycles is reached, the linear deformation and nonlinear deformation of the last iterative cycle are combined to obtain the temporal deformation; otherwise, the PS candidate pixels and high coherence DS candidate pixels whose temporal coherence coefficients are less than the set coherence threshold are eliminated, and steps 5 to 9 are repeated.

[0016] In step 1 above, the main image and the non-main image are obtained based on the following steps:

[0017] The multi-temporal SAR image data with the largest calculated coherence coefficient with other multi-temporal SAR image data is selected as the main image, and other multi-temporal SAR image data are selected as non-main images.

[0018] As described above in step 4, the high coherence DS candidate pixels are obtained based on the following steps:

[0019] Step 4.1: Select reference pixels in the differential interferometric data, determine homogeneous pixels of the reference pixels, determine reference pixels as DS candidate pixels based on the homogeneous pixels, and calculate the weighted coherence matrix T of the DS candidate pixels;

[0020] Step 4.2: Perform singular value decomposition on the weighted coherence matrix T of the DS candidate pixel to obtain the optimized phase vector of the DS candidate pixel;

[0021] Step 4.3: Substitute the optimized phase vector of the DS candidate pixel into the phase fit goodness test model to obtain the phase fit goodness result, and select high coherence DS candidate pixels whose phase fit goodness result is lower than the test threshold.

[0022] As described above in step 4.1, determining the homogeneous pixels of the reference pixels includes the following steps:

[0023] Calculate the image block structure similarity measure D(x,y)

[0024]

[0025] Where x represents the reference pixel in the differential interferometry data, y represents the neighborhood pixel in the 11×11 window neighborhood centered on the reference pixel x in the differential interferometry data, a 3×3 window of the reference pixel is constructed with the reference pixel as the center, and the pixels in the 3×3 window centered on the reference pixel are the reference extended pixels, a 3×3 window of the neighborhood pixel is constructed with the neighborhood pixel as the center, and the pixels in the 3×3 window centered on the neighborhood pixel are the neighborhood extended pixels, i is the number of the pixel in the 3×3 window, x is the number of the pixel in the 3×3 window, i and i are the differential interferometry data vector of the reference extended pixel corresponding to pixel number i and the differential interferometry data vector of the neighborhood extended pixel, E() represents the expected calculation of the sample, * represents conjugate, and || represents the modulus of the complex number.

[0026] If D(x,y) is greater than the threshold D thresh , the corresponding neighborhood pixel y is used as the homogeneous pixel of the reference pixel x; otherwise, the corresponding neighborhood pixel y is used as the non-homogeneous pixel of the reference pixel x.

[0027] As described above in step 4.1, the DS candidate pixels are determined based on the following steps:

[0028] When the number of neighborhood pixels adjacent to a reference pixel that are identified as homogeneous pixels exceeds a set number, the reference pixel is a DS candidate pixel.

[0029] As described above in step 4.1, the weighted coherence matrix T is obtained based on the following steps:

[0030] Calculate the similarity weight coefficient w between DS candidate pixels and homogeneous pixels:

[0031]

[0032] Calculate the weighted coherence matrix T of the DS candidate pixels:

[0033]

[0034] Among them, w jis the similarity weight coefficient of the homogeneous pixel numbered j of the DS candidate pixel, y j represents the homogeneous pixel numbered j of the DS candidate pixel, Ω is the homogeneous pixel set of the DS candidate pixel, N p represents the number of homogeneous pixels of DS candidate pixels, and H represents the conjugate operation of the vector.

[0035] A computer device comprises a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps of the measurement method according to any one of claims 1 to 6 when executing the computer program.

[0036] A computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the measurement method according to any one of claims 1 to 6.

[0037] A computer program product comprises a computer program, which implements the steps of the measurement method according to any one of claims 1 to 6 when executed by a processor.

[0038] Compared with the prior art, the present invention has the following advantages:

[0039] 1. High density. Compared with the traditional PSI, KS-DSI, FaSHP-DSI and other classic time-series InSAR analysis methods, the number of monitoring points of the InSAR deformation monitoring results of the ascending and descending orbits of the present invention is the highest, which is 1.7 and 2.8 times that of the currently widely used FaSHP-DSI method.

[0040] 2. Wide coverage. Compared with the traditional PSI, KS-DSI, FaSHP-DSI and other classic time-series InSAR analysis methods, the InSAR deformation monitoring results of the ascending and descending orbits of the present invention have the widest coverage.

[0041] 3. High precision. The MAE and RMSE calculated by the difference between the InSAR deformation monitoring results of the ascending and descending orbits and the GNSS deformation results are the smallest, which are better than 1.5 mm / yr and 3.5 mm / yr respectively. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a schematic diagram of the process of the present invention.

[0043] Figure 2 Figures 2 and 3 are deformation rate result diagrams, where (a), (b), (c), and (d) are the PSI, KS-DSI, FaSHP-DSI, and deformation rate results monitored by the present invention of the ascending InSAR, respectively.

[0044] Figure 3Figures 2 and 3 are deformation rate result diagrams, where (a), (b), (c), and (d) are the deformation rate results monitored by the PSI, KS-DSI, FaSHP-DSI, and the present invention of descending InSAR, respectively. DETAILED DESCRIPTION

[0045] In order to facilitate those skilled in the art to understand and implement the present invention, the present invention is further described in detail below in conjunction with embodiments. The embodiments described herein are only used to illustrate and explain the present invention and are not intended to limit the present invention.

[0046] Example:

[0047] This embodiment proposes a radar distributed scatterer interferometry method for high and steep areas, and the input data includes a multi-phase SAR image data set and a known rough-precision digital elevation model (DEM). Figure 1 As shown:

[0048] Step 1: Based on the comprehensive coherence coefficient theory, the multi-temporal SAR image data with the largest coherence coefficient calculation result with other multi-temporal SAR image data is selected as the main image, and other multi-temporal SAR image data are selected as non-main images;

[0049] Step 2: Perform conjugate operation between the main image and the non-main image, and remove the terrain phase introduced by the rough precision digital elevation model to generate differential interferometry data;

[0050] Step 3: Calculate the amplitude deviation index of the differential interferometer data, set the amplitude deviation index threshold, and select the PS candidate pixels whose amplitude deviation index is less than the amplitude deviation threshold;

[0051] Step 4: Use the differential interference data generated in step 2 to screen out high coherence DS candidate pixels. The specific steps are as follows:

[0052] Step 4.1: Use the image block structure similarity measure D(x,y) in formula (1) to identify homogeneous pixels;

[0053]

[0054] Where x represents the reference pixel in the differential interferometry data, y represents the neighborhood pixel in the 11×11 window neighborhood centered on the reference pixel x in the differential interferometry data, a 3×3 window of the reference pixel is constructed with the reference pixel as the center, and the pixels in the 3×3 window centered on the reference pixel are the reference extended pixels, a 3×3 window of the neighborhood pixel is constructed with the neighborhood pixel as the center, and the pixels in the 3×3 window centered on the neighborhood pixel are the neighborhood extended pixels, i is the number of the pixel in the 3×3 window, x is the number of the pixel in the 3×3 window, i and iare the differential interferometric data vector of the reference extended pixel corresponding to pixel number i and the differential interferometric data vector of the neighborhood extended pixel, E() represents the expected calculation of the sample, * represents conjugate, and || represents the modulus of the complex number. If D(x,y) is greater than a certain threshold D thresh , the corresponding neighborhood pixel y is used as the homogeneous pixel of the reference pixel x; otherwise, the corresponding neighborhood pixel y is used as the non-homogeneous pixel of the reference pixel x.

[0055] When the number of neighboring pixels adjacent to the reference pixel that are identified as homogeneous pixels exceeds the set number (such as 20), the reference pixel becomes a DS candidate pixel, and the similarity weight coefficient w between the DS candidate pixel and the homogeneous pixel is calculated:

[0056]

[0057] It can be seen from this formula that when D(x, y) = 0, w = 1; when D(x, y) ≥ D thresh , w = 0. A weight coefficient is assigned to each homogeneous pixel of the screened DS candidate pixels, and the weighted coherence matrix T of the DS candidate pixels is calculated using the following formula:

[0058]

[0059] Among them, w j is the similarity weight coefficient of the homogeneous pixel numbered j of the DS candidate pixel, y j represents the homogeneous pixel numbered j of the DS candidate pixel, Ω is the homogeneous pixel set of the DS candidate pixel, N p represents the number of homogeneous pixels of DS candidate pixels, and H represents the conjugate operation of the vector.

[0060] Step 4.2: Perform singular value decomposition on the weighted coherence matrix T of the DS candidate pixel, use the strongest scattering mechanism obtained by the singular value decomposition as the optimal phase vector, and obtain the optimized phase vector of the DS candidate pixel. The details are as follows:

[0061] The weighted coherence matrix T is subjected to singular value decomposition:

[0062]

[0063] Among them, λ i ~λ M is the eigenvalue of the weighted coherence matrix T, sorted in descending order λ1≥λ2≥L≥λ M , μ k is the eigenvalue λ k The corresponding eigenvector, k is the sequence number of the eigenvalue. At this time, the phase vector of the DS candidate pixel after optimization |||| is the 2-norm operator of a vector.

[0064] Step 4.3: Substitute the optimized phase vector of the DS candidate pixel in step 4.2 into the phase fit goodness test model to obtain the phase fit goodness result. According to the test threshold, select high coherence DS candidate pixels whose phase fit goodness result is lower than the test threshold.

[0065] Step 5: Using the PS candidate pixels and the high coherence DS candidate pixels, a distance-constrained Delaunay triangulation network is constructed that fuses the PS candidate pixels and the high coherence DS candidate pixels;

[0066] Step 6: Maximize the set coherence model between two points connected by the distance-constrained Delaunay triangulation network edge, and obtain the differential residual elevation, differential deformation rate, set coherence coefficient of the connected edge, and spatial differential residual phase;

[0067] Step 7: Perform three-dimensional phase unwrapping on the spatial differential residual phase obtained in step 6 to obtain residual topography, linear deformation, deformation rate and unwrapped residual phase;

[0068] Step 8: Using the different spatiotemporal characteristics of the atmospheric orbit phase, nonlinear deformation phase, and noise phase, the unwrapped residual phase is bandpass filtered to obtain the atmospheric orbit phase, nonlinear deformation phase, and noise phase;

[0069] Step 9: Calculate the temporal coherence coefficients of the PS candidate pixels and the high coherence DS candidate pixels according to the noise phase. If the temporal coherence coefficients of the PS candidate pixels and the high coherence DS candidate pixels are both greater than the set coherence threshold, or the maximum number of iterative cycles is reached, then according to steps 7 and 8 of the last iterative cycle, linear deformation and nonlinear deformation are combined to obtain temporal deformation; otherwise, remove the PS candidate pixels and the high coherence DS candidate pixels whose temporal coherence coefficients are less than the set coherence threshold, and repeat steps 5 to 9.

[0070] A person of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing related hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods.

[0071] In this embodiment, the multi-phase SAR image data of the ascending and descending orbits of the European Space Agency Sentinel-1 that coincide with the start time of the water storage period in the area are used to verify the method of the present invention. The data time span is more than 2 years, starting from the first SAR image data when the Baihetan Hydropower Station began to store water after its completion. The ascending orbit time span is from April 9, 2021 to May 29, 2023, including a total of 66 images; the descending orbit time span is from April 17, 2021 to May 31, 2023, including a total of 77 images. In order to reduce the impact of terrain on deformation, the DEM data product of TanDEM-X is selected for reference terrain phase removal. Verify the effectiveness of the present invention.

[0072] A comparative analysis of the InSAR deformation rate of the entire Baihetan reservoir area was conducted. Figure 2 Figure 2 is the deformation rate result diagram. (a), (b), (c), and (d) are the deformation rate results monitored by PSI, KS-DSI, FaSHP-DSI, and the present invention respectively. By comparison, it can be found that the spatial distribution of deformation results of several methods is similar, but the monitoring point density of PSI and KS-DSI is low. FaSHP has a higher monitoring point density, but the monitoring coverage is low in the high and steep terrain areas along the banks of the Jinsha River. The proposed method has the highest monitoring point density and the best coverage in the high and steep terrain areas.

[0073] The deformation inversion performance of several methods is quantitatively evaluated and compared by comparing the GNSS monitoring results of three locations of the Eighth Hydropower Bureau with the InSAR deformation rate results of each method. The evaluation indicators are mean absolute error (MAE) and root mean square error (RMSE). Table 1 shows the quantitative evaluation results of the comparison between the deformation rate of the ascending and descending orbit InSAR of PSI, KS-DSI, FaSHP-DSI, and the present invention and the GNSS deformation rate, as well as the statistics of the number of InSAR phase change monitoring points in the area. In Table 1, the uplink data corresponding to the same method is the data corresponding to the ascending orbit, and the downlink data is the data corresponding to the descending orbit. It can be seen from the comparison that the ascending and descending orbit results of the present invention are the best, with the deformation accuracy of the ascending orbit being better than 1.5 mm / year and the accuracy of the descending orbit being better than 3.3 mm / year; the number of ascending and descending orbit monitoring points of the proposed method is 1.7 and 2.8 times that of the currently most widely used FaSHP-DSI method.

[0074] Table 1

[0075]

[0076] Embodiment 2:

[0077] This embodiment provides a computer device, including a memory and a processor. The memory stores a computer program. When the processor executes the computer program, each step in the above-mentioned embodiment 1 is implemented.

[0078] Embodiment 3:

[0079] This embodiment provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, each step in the above-mentioned embodiment 1 is implemented.

[0080] Embodiment 4:

[0081] This embodiment provides a computer program product, including a computer program. When the computer program is executed by a processor, each step in the above-mentioned embodiment 1 is implemented.

[0082] The above are only preferred embodiments of the present invention, and the protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, some improvements and modifications without departing from the principle of the present invention should be regarded as the protection scope of the present invention.

Claims

1. A radar distributed scatterer interferometry method for high and steep areas, characterized in that: The following steps are involved: Step 1: Select the main image and non-main image from the multi-temporal SAR image data; Step 2: Perform conjugate operation between the main image and the non-main image to remove the terrain phase introduced by the digital elevation model and generate differential interferometry data; Step 3: Calculate the amplitude deviation index of the differential interferometer data, and select the PS candidate pixels whose amplitude deviation index is less than the amplitude deviation threshold; Step 4: Screen out high coherence DS candidate pixels based on differential interferometry data; Step 5: Construct a distance-constrained Delaunay triangulation network that fuses PS candidate pixels with high-coherence DS candidate pixels; Step 6: Maximize the set coherence model between two points connected by the distance-constrained Delaunay triangulation network edge, and obtain the differential residual elevation, differential deformation rate, set coherence coefficient of the connected edge, and spatial differential residual phase; Step 7: Perform three-dimensional phase unwrapping on the spatial differential residual phase obtained in step 6 to obtain residual topography, linear deformation, deformation rate and unwrapped residual phase; Step 8: Bandpass filter the unwrapped residual phase to obtain the atmospheric orbit phase, nonlinear deformation and noise phase; Step 9: If the temporal coherence coefficients of the PS candidate pixels and the high coherence DS candidate pixels are both greater than the set coherence threshold, or the maximum number of iterative cycles is reached, the linear deformation and nonlinear deformation of the last iterative cycle are combined to obtain the temporal deformation; otherwise, the PS candidate pixels and high coherence DS candidate pixels whose temporal coherence coefficients are less than the set coherence threshold are eliminated, and steps 5 to 9 are repeated.

2. According to claim 1, a radar distributed scatterer interferometry method for high and steep areas is characterized in that: In step 1, the main image and the non-main image are obtained based on the following steps: The multi-temporal SAR image data with the largest calculated coherence coefficient with other multi-temporal SAR image data is selected as the main image, and other multi-temporal SAR image data are selected as non-main images.

3. According to claim 1, a radar distributed scatterer interferometry method for high and steep areas is characterized in that: In step 4, the high coherence DS candidate pixels are obtained based on the following steps: Step 4.1: Select reference pixels in the differential interferometric data, determine homogeneous pixels of the reference pixels, determine reference pixels as DS candidate pixels based on the homogeneous pixels, and calculate the weighted coherence matrix T of the DS candidate pixels; Step 4.2: Perform singular value decomposition on the weighted coherence matrix T of the DS candidate pixel to obtain the optimized phase vector of the DS candidate pixel; Step 4.3: Substitute the optimized phase vector of the DS candidate pixel into the phase fit goodness test model to obtain the phase fit goodness result, and select high coherence DS candidate pixels whose phase fit goodness result is lower than the test threshold.

4. According to claim 3, a radar distributed scatterer interferometry method for high and steep areas is characterized in that: In step 4.1, determining the homogeneous pixels of the reference pixels comprises the following steps: Calculate the image block structure similarity measure D(x,y) Where x represents the reference pixel in the differential interferometry data, y represents the neighborhood pixel in the 11×11 window neighborhood centered on the reference pixel x in the differential interferometry data, a 3×3 window of the reference pixel is constructed with the reference pixel as the center, and the pixels in the 3×3 window centered on the reference pixel are the reference extended pixels, a 3×3 window of the neighborhood pixel is constructed with the neighborhood pixel as the center, and the pixels in the 3×3 window centered on the neighborhood pixel are the neighborhood extended pixels, i is the number of the pixel in the 3×3 window, x is the number of the pixel in the 3×3 window, i and i are the differential interferometry data vector of the reference extended pixel corresponding to pixel number i and the differential interferometry data vector of the neighborhood extended pixel, E() represents the expected calculation of the sample, * represents conjugate, and | | represents the modulus of the complex number. If D(x,y) is greater than the threshold D thresh , the corresponding neighborhood pixel y is used as the homogeneous pixel of the reference pixel x; otherwise, the corresponding neighborhood pixel y is used as the non-homogeneous pixel of the reference pixel x.

5. The radar distributed scatterer interferometry method for high and steep areas according to claim 4, characterized in that: In step 4.1, the DS candidate pixels are determined based on the following steps: When the number of neighborhood pixels adjacent to a reference pixel that are identified as homogeneous pixels exceeds a set number, the reference pixel is a DS candidate pixel.

6. The radar distributed scatterer interferometry method for high and steep areas according to claim 4, characterized in that: In step 4.1, the weighted coherence matrix T is obtained based on the following steps: Calculate the similarity weight coefficient w between DS candidate pixels and homogeneous pixels: Calculate the weighted coherence matrix T of the DS candidate pixels: Among them, w j is the similarity weight coefficient of the homogeneous pixel numbered j of the DS candidate pixel, y j represents the homogeneous pixel numbered j of the DS candidate pixel, Ω is the homogeneous pixel set of the DS candidate pixel, N p It represents the number of homogeneous pixels of DS candidate pixels, and H represents the conjugate operation of the vector.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, the steps of the measurement method according to any one of claims 1 to 6 are implemented.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the measurement method according to any one of claims 1 to 6 are implemented.

9. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the measurement method according to any one of claims 1 to 6 are implemented.

Images