A radar distributed scatterer interferometry method and apparatus for high-relief areas

Abstract

This invention discloses a radar distributed scatterer interferometry method and device for steep terrain. The method involves selecting a main image and non-main images; calculating differential interferometric data; screening PS candidate pixels with amplitude deviation indices less than the amplitude deviation threshold; selecting highly coherent DS candidate pixels; constructing a distance-constrained Delaunay triangulation network by fusing PS candidate pixels and highly coherent DS candidate pixels; calculating differential residual elevation, differential deformation rate, ensemble coherence coefficient of connected edges, and spatial differential residual phase; calculating residual terrain, linear deformation, deformation rate, and unwrapped residual phase; calculating atmospheric orbit phase, nonlinear deformation, and noise phase; and obtaining temporal deformation based on linear and nonlinear deformation. This invention provides a high number of monitoring points and the widest coverage of InSAR deformation monitoring results, with low MAE and RMSE.

Description

Technical Field

[0001] This invention belongs to the field of interferometric synthetic aperture radar data processing technology, specifically relating to a radar distributed scatterer interferometry method and equipment for steep regions. Background Technology

[0002] DS-InSAR is a rapidly developing new geodetic technology in recent years. It fully combines the high temporal coherence of DS targets after spatial adaptive filtering with the high temporal coherence of PS targets under single-look conditions for interferometric point target analysis. It can effectively overcome the shortcomings of PSInSAR in monitoring natural scenes, while inheriting the advantages of PSI in high-precision and high-resolution deformation monitoring. Since its introduction in 2011, DS-InSAR has become one of the key development technologies of temporal InSAR, providing an important means for precise deformation monitoring in complex scenes such as volcanoes, mining areas, landslides, and urban areas.

[0003] Identification of homogeneous points and phase optimization are key technologies in DS-InSAR. The first-generation DS-InSAR technology, SqueeSAR, was developed in this field. TM The method first employs an amplitude-based two-sample Kolmogorov–Smirnov (KS) hypothesis testing method for homogeneous point identification; then, based on these homogeneous pixels, the covariance matrix is ​​estimated; finally, a 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 mechanisms. While this increases the spatial density and coverage of monitoring points, it also suffers from insufficient robustness with small samples and low computational efficiency, limiting the widespread application of DS-InSAR technology. Summary of the Invention

[0004] The purpose of this invention is to address the problems existing in the prior art by providing a radar distributed scatterer interferometry method and device for high and steep regions. This aims to improve the monitoring coverage, density, and reliability in high mountain and canyon areas.

[0005] The above-mentioned objectives of the present invention are achieved through the following technical means:

[0006] A radar distributed scatterer interferometry method for steep terrain includes the following steps:

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

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

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

[0010] Step 4: Based on differential interferometric data, select high-coherence DS candidate pixels;

[0011] Step 5: Construct a distance-constrained Delaunay triangulation network for fusing PS candidate pixels and highly coherent DS candidate pixels;

[0012] Step 6: Maximize the set coherence model between two points connected by the edge of the distance-constrained Delaunay triangular network to 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 difference residual phase obtained in Step 6 to obtain the residual terrain, linear deformation, deformation rate and unwrapped residual phase;

[0014] Step 8: Perform bandpass filtering on the unwrapped residual phase to obtain the atmospheric orbit phase, nonlinear deformation, and noise phase;

[0015] Step 9: If the temporal coherence coefficients of both the PS candidate pixel and the highly coherent DS candidate pixel are greater than the set coherence threshold, or the maximum number of iterations is reached, then the linear deformation and nonlinear deformation of the last iteration are combined to obtain the temporal deformation; otherwise, PS candidate pixels and highly coherent DS candidate pixels with temporal coherence coefficients less than the set coherence threshold are removed, and steps 5 to 9 are repeated.

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

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

[0018] As described above in step 4, 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 homogeneous pixels, and calculate the weighted coherence matrix T of DS candidate pixels;

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

[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 screen high coherence DS candidate pixels with phase fit goodness results lower than the test threshold.

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

[0023] Calculate the image patch structural similarity measure D(x,y)

[0024]

[0025] Where x represents the reference pixel in the differential interferometric data, and y represents the neighboring pixels within an 11×11 window centered on the reference pixel x. A 3×3 window is constructed centered on the reference pixel, and the pixels within this 3×3 window are called reference extended pixels. Similarly, a 3×3 window is constructed centered on the neighboring pixels, and the pixels within this 3×3 window are called neighbor extended pixels. Let i be the pixel number within the 3×3 window, and x... i and y i These are the differential interferometric data vectors of the reference extended pixel and the neighboring extended pixel, respectively. E() represents the expected value 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 If the corresponding neighboring pixel y is considered a homogeneous pixel of the reference pixel x, then the corresponding neighboring pixel y is considered a non-homogeneous pixel of the reference pixel x.

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

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

[0029] As described in step 4.1 above, 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 candidate pixels in DS:

[0033]

[0034] Among them, w jy is the similarity weight coefficient of homogeneous pixels with the same number j among the DS candidate pixels. j Let Ω represent the homogeneous pixel numbered j of the DS candidate pixel, Ω be the set of homogeneous pixels of the DS candidate pixel, and N be the set of homogeneous pixels of the DS candidate pixel. p H represents the number of homogeneous pixels among the candidate pixels in DS, and H represents the conjugate operation of the vector.

[0035] A computer device includes a memory and a processor, the memory storing a computer program, the processor executing the computer program to implement the steps of the measurement method according to any one of claims 1 to 6.

[0036] A computer-readable storage medium having a computer program stored thereon, which, 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 includes a computer program that, when executed by a processor, implements the steps of the measurement method according to any one of claims 1 to 6.

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

[0039] 1. High density. Compared with the traditional time-series InSAR analysis methods such as PSI, KS-DSI, and FaSHP-DSI, the number of monitoring points for the rising and falling orbit InSAR deformation monitoring results of this 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 traditional time-series InSAR analysis methods such as PSI, KS-DSI, and FaSHP-DSI, the InSAR deformation monitoring results of the rising and falling orbits of this invention have the widest coverage.

[0041] 3. High precision. The MAE and RMSE calculated by comparing the InSAR deformation monitoring results and the GNSS deformation results of the lifting rail in this invention are the smallest, which are better than 1.5mm / yr and 3.5mm / yr, respectively. Attached Figure Description

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

[0043] Figure 2 The figure shows the deformation rate results. (a), (b), (c), and (d) are the PSI, KS-DSI, FaSHP-DSI, and deformation rate results monitored by the present invention for ascending InSAR, respectively.

[0044] Figure 3The figure shows the deformation rate results. (a), (b), (c), and (d) are the PSI, KS-DSI, FaSHP-DSI, and deformation rate results monitored by the present invention for the down-orbit InSAR, respectively. Detailed Implementation

[0045] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to embodiments. The embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0046] Example:

[0047] This embodiment proposes a radar distributed scatterer interferometry method for steep terrain. The input data includes a multi-temporal SAR image dataset and a known coarse-precision digital elevation model (DEM). The specific flowchart is shown below. Figure 1 As shown:

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

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

[0050] Step 3: Calculate the amplitude deviation index of the differential interferometric data, set the amplitude deviation index threshold, and screen PS candidate pixels with amplitude deviation index less than the amplitude deviation threshold;

[0051] Step 4: Using the differential interferometric data generated in Step 2, select 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) of equation (1) to identify homogeneous pixels;

[0053]

[0054] Where x represents the reference pixel in the differential interferometric data, and y represents the neighboring pixels within an 11×11 window centered on the reference pixel x. A 3×3 window is constructed centered on the reference pixel, and the pixels within this 3×3 window are called reference extended pixels. Similarly, a 3×3 window is constructed centered on the neighboring pixels, and the pixels within this 3×3 window are called neighbor extended pixels. Let i be the pixel number within the 3×3 window, and x... i and y iLet E(x,y) be the differential interferometric data vector of the reference extended pixel and the differential interferometric data vector of the neighboring extended pixel, respectively, and let E() represent the expected value of the sample, * denote conjugate, and || denote the modulus of the complex number. If D(x,y) is greater than a certain threshold D thresh If the neighboring pixel y is homogeneous with the reference pixel x, then the corresponding neighboring pixel y is considered a homogeneous pixel of the reference pixel x; otherwise, the corresponding neighboring pixel y is considered a 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 a set number (e.g., 20), the reference pixel becomes a DS candidate pixel. The similarity weight coefficient w between the DS candidate pixel and the homogeneous pixels is then calculated.

[0056]

[0057] From this formula, we can see that: when D(x,y)=0, w=1; when D(x,y)≥D thresh w = 0. Assign a weight coefficient to each homogeneous pixel in the selected DS candidate pixels, and calculate the weighted coherence matrix T of the DS candidate pixels using the following formula:

[0058]

[0059] Among them, w j y is the similarity weight coefficient of homogeneous pixels with the same number j among the DS candidate pixels. j Let Ω represent the homogeneous pixel numbered j of the DS candidate pixel, Ω be the set of homogeneous pixels of the DS candidate pixel, and N be the set of homogeneous pixels of the DS candidate pixel. p H represents the number of homogeneous pixels among the candidate pixels in DS, 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 pixels. Use the strongest scattering mechanism obtained from the singular value decomposition as the optimal phase vector to obtain the optimized phase vector of the DS candidate pixels. Specifically:

[0061] Singular value decomposition of the weighted coherence matrix T:

[0062]

[0063] Where, λ i ~λ M The eigenvalues ​​of the weighted coherence matrix T are sorted in descending order: λ1≥λ2≥L≥λ M μ k For the eigenvalue λ k The corresponding feature vector, where k is the index of the feature value. At this point, the optimized phase vector of the DS candidate pixel... |||| is the 2-norm operator for vectors.

[0064] Step 4.3: Substitute the optimized phase vector of the DS candidate pixels from Step 4.2 into the phase fit goodness test model to obtain the phase fit goodness result. Based on the test threshold, select high-coherence DS candidate pixels with phase fit goodness results lower than the test threshold.

[0065] Step 5: Construct a distance-constrained Delaunay triangulation network by fusing PS candidate pixels and high-coherence DS candidate pixels using PS candidate pixels and high-coherence DS candidate pixels.

[0066] Step 6: Maximize the set coherence model between two points connected by the edge of the distance-constrained Delaunay triangular network to 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 difference residual phase obtained in Step 6 to obtain the residual terrain, linear deformation, deformation rate and unwrapped residual phase;

[0068] Step 8: Utilize the different spatiotemporal characteristics of atmospheric orbit phase, nonlinear deformation phase, and noise phase to perform bandpass filtering on the unwrapped residual phase to obtain the atmospheric orbit phase, nonlinear deformation phase, and noise phase;

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

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

[0071] In this embodiment, multi-temporal SAR imagery data from the ESA Sentinel-1 system, coinciding with the start of the impoundment period in the region, is used to verify the method of the present invention. The data spans more than two years, beginning with the first SAR imagery data taken after the Baihetan Hydropower Station began impounding water. The ascent period spans from April 9, 2021 to May 29, 2023, comprising 66 images; the descent period spans from April 17, 2021 to May 31, 2023, comprising 77 images. To reduce the impact of terrain on deformation, TanDEM-X DEM data products are selected for reference terrain phase removal. This verifies the effectiveness of the present invention.

[0072] A comparative analysis of the InSAR deformation rate across the entire Baihetan Reservoir area was conducted. Figure 2 The figures (a), (b), (c), and (d) show the deformation rate results for ascending-orbit InSAR (PSI), KS-DSI, FaSHP-DSI, and the deformation rate monitored by this invention, respectively. A comparison reveals that the spatial distribution of deformation results from these methods is similar, but PSI and KS-DSI have lower monitoring point densities. FaSHP has a higher monitoring point density, but its coverage is low in the steep terrain along both banks of the Jinsha River. The proposed method has the highest monitoring point density and the best coverage in steep terrain areas.

[0073] The deformation rate results of various InSAR methods were compared with the monitoring results of three GNSS stations at the Eighth Hydropower Bureau to quantitatively evaluate and compare the deformation inversion performance of several methods. The evaluation indicators were 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 PSI, KS-DSI, FaSHP-DSI, and the proposed method's ascending and descending orbit InSAR and GNSS deformation rates, as well as the statistics of the number of InSAR phase change monitoring points in the region. In Table 1, the uplink data corresponding to the same method is the data corresponding to ascending orbit, and the downlink data is the data corresponding to descending orbit. The comparison shows that the ascending and descending orbit results of the proposed method are the best, with ascending orbit deformation accuracy better than 1.5 mm / year and descending orbit accuracy 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, respectively.

[0074] Table 1

[0075]

[0076] Example 2:

[0077] This embodiment provides a computer device, including a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the various steps in Embodiment 1 above.

[0078] Example 3:

[0079] This embodiment provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps in Embodiment 1 above.

[0080] Example 4:

[0081] This embodiment provides a computer program product, including a computer program that, when executed by a processor, implements the steps in Embodiment 1 above.

[0082] The above are merely preferred embodiments of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A radar distributed scatterer interferometry method for steep terrain, characterized in that, Includes the following steps: Step 1: Select the main image and non-main images from the multi-temporal SAR image data; Step 2: Perform conjugate operation between the main image and non-main images to remove the terrain phase introduced by the digital elevation model and generate differential interferometric data; Step 3: Calculate the amplitude deviation index of the differential interferometric data and screen PS candidate pixels whose amplitude deviation index is less than the amplitude deviation threshold; Step 4: Based on differential interferometric data, select high-coherence DS candidate pixels; Step 5: Construct a distance-constrained Delaunay triangulation network for fusing PS candidate pixels and highly coherent DS candidate pixels; Step 6: Maximize the set coherence model between two points connected by the edge of the distance-constrained Delaunay triangular network to 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 difference residual phase obtained in Step 6 to obtain the residual terrain, linear deformation, deformation rate and unwrapped residual phase; Step 8: Perform bandpass filtering on the unwrapped residual phase to obtain the atmospheric orbit phase, nonlinear deformation, and noise phase; Step 9: If the temporal coherence coefficients of both the PS candidate pixel and the highly coherent DS candidate pixel are greater than the set coherence threshold, or the maximum number of iterations is reached, then the linear deformation and nonlinear deformation of the last iteration are combined to obtain the temporal deformation; otherwise, PS candidate pixels and highly coherent DS candidate pixels with temporal coherence coefficients less than the set coherence threshold are removed, and steps 5 to 9 are repeated.

2. The radar distributed scatterer interferometry method for steep terrain as described in claim 1, characterized in that, In step 1, the main image and non-main images are obtained based on the following steps: The multi-temporal SAR image with the highest coherence coefficient among other multi-temporal SAR image data is selected as the main image, and the other multi-temporal SAR image data are selected as non-main images.

3. The radar distributed scatterer interferometry method for steep terrain as described in claim 1, characterized in that, In step 4, 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 homogeneous pixels, and calculate the weighted coherence matrix T of DS candidate pixels; Step 4.2: Perform singular value decomposition on the weighted coherence matrix T of the DS candidate pixels to obtain the optimized phase vector of the DS candidate pixels; 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 screen high coherence DS candidate pixels with phase fit goodness results lower than the test threshold.

4. The radar distributed scatterer interferometry method for steep terrain as described in claim 3, characterized in that, In step 4.1, determining the homogeneous pixels of the reference pixel includes the following steps: Calculate the image patch structural similarity measure D(x,y) Where x represents the reference pixel in the differential interferometric data, and y represents the neighboring pixels within an 11×11 window centered on the reference pixel x. A 3×3 window is constructed centered on the reference pixel, and the pixels within this 3×3 window are called reference extended pixels. Similarly, a 3×3 window is constructed centered on the neighboring pixels, and the pixels within this 3×3 window are called neighbor extended pixels. Let i be the pixel number within the 3×3 window, and x... i and y i These are the differential interferometric data vectors of the reference extended pixel and the neighboring extended pixel, respectively. E() represents the expected value calculation of the sample, * represents the conjugate, and || represents the modulus of the complex number. If D(x,y) is greater than the threshold D thresh If the corresponding neighboring pixel y is considered a homogeneous pixel of the reference pixel x, then the corresponding neighboring pixel y is considered a non-homogeneous pixel of the reference pixel x.

5. The radar distributed scatterer interferometry method for steep terrain as described in claim 4, characterized in that, In step 4.1, DS candidate pixels are determined based on the following steps: When the number of neighboring pixels adjacent to the 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 steep terrain as described in 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 y is the similarity weight coefficient of homogeneous pixels with the same number j among the DS candidate pixels. j Let Ω represent the homogeneous pixel numbered j of the DS candidate pixel, Ω be the set of homogeneous pixels of the DS candidate pixel, and N be the set of homogeneous pixels of the DS candidate pixel. p H represents the number of homogeneous pixels among the candidate pixels in DS, 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, it implements the steps of the measurement method according to any one of claims 1 to 6.

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

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

Images