Network optimization InSAR large gradient deformation phase unwrapping method based on edge detection

Abstract

The invention discloses a network optimization InSAR large gradient deformation phase unwrapping method based on edge detection, and the method specifically comprises the steps: 1, carrying out the differential interference processing of an original image, and obtaining a differential interference pattern; 2, Gaussian filtering is carried out on the differential interferogram, and the gradient magnitude and the gradient direction of each pixel in the image after Gaussian filtering are calculated; a gradient magnitude image is obtained; step 3, performing non-maximum suppression on the gradient magnitude image, and completing preliminary extraction of edges; 4, setting double thresholds to further screen edge points, and finally completing edge extraction; 5, generating a coherence graph based on the differential interferogram, extracting coherence points, and generating an unwrapping network based on the extracted coherence points and the edges extracted in the step 4; and step 6, phase unwrapping is carried out based on the unwrapping network. The method has certain universality for monitoring application scenes with large-gradient characteristic deformation, such as landslide, mining and the like, and can provide important technical support for geological disaster prevention and control.

Description

Technical Field

[0001] This invention belongs to the field of synthetic aperture radar interferometry technology, and particularly relates to a network-optimized InSAR large gradient deformation phase unwrapping method based on edge detection. Background Technology

[0002] InSAR phase unwrapping is crucial to the accuracy of deformation results. For phase unwrapping in large-gradient deformation scenarios, most methods exhibit errors due to the presence of arcs that do not satisfy the phase continuity assumption. Some studies utilize external information constraints to reduce the connection of large-gradient deformation arcs, ensuring overall unwrapping accuracy. Other studies focus on using offset tracking or model inversion deformation to provide constraints or prior information for unwrapping, thereby obtaining high-precision unwrapped phases. However, single-baseline phase unwrapping is still limited by the phase continuity assumption. Furthermore, for large-gradient deformation regions, phase jumps in the unwrapped phase often occur because of edges that do not satisfy the phase continuity assumption. Conventionally constructed phase unwrapping networks do not fully consider this issue; directly connecting high-coherence points behind the mask generates arcs that cross these edges. This inevitably leads to these arcs that do not satisfy phase continuity still being considered in subsequent solutions, resulting in unwrapping errors and ultimately, obvious deformation jump regions in the unwrapped results. Summary of the Invention

[0003] Purpose of the invention: In order to solve the problems existing in the prior art, the present invention provides a network-optimized InSAR large gradient deformation phase unwrapping method based on edge detection.

[0004] Technical solution: This invention provides a network-optimized InSAR large gradient deformation phase unwrapping method based on edge detection, specifically including the following steps:

[0005] Step 1: Perform differential interferometry on the original image to obtain a differential interferogram;

[0006] Step 2: Perform Gaussian filtering on the differential interferogram and calculate the gradient magnitude and gradient direction of each pixel in the Gaussian-filtered image; obtain the gradient magnitude image;

[0007] Step 3: Perform non-maximum suppression on the gradient magnitude image and complete the initial edge extraction;

[0008] Step 4: Set dual thresholds to further filter edge points, and finally complete the edge extraction;

[0009] Step 5: Generate a coherence map based on the differential interferogram and extract coherence points. Generate an unwrapped network based on the extracted coherence points and the edges extracted in Step 4.

[0010] Step 6: Perform phase unwrapping based on the unwrapped network.

[0011] Furthermore, in step 2, the gradient magnitude and gradient direction of each pixel are calculated using the following formula:

[0012]

[0013] Where (x,y) are the coordinates of the pixel, M(x,y) represents the gradient magnitude of the pixel, θ(x,y) is the gradient direction of the pixel, and d x (x,y) represents the gradient value of the image along the x-axis at pixel location (x,y), d y (x,y) represents the gradient value of the image along the y-axis at the pixel location (x,y).

[0014] Furthermore, step 3 specifically involves: for any pixel P, along the gradient direction of the pixel, obtaining two adjacent pixels of the pixel, comparing the gradient magnitude of pixel P with that of its two adjacent pixels, and if the gradient magnitude of pixel P is greater than the gradient magnitude of its two adjacent pixels, then pixel P is assigned the value 1; otherwise, pixel P is assigned the value 0.

[0015] Furthermore, step 4 specifically involves: setting a first threshold and a second threshold, with the first threshold being greater than the second threshold; classifying pixels with gradient magnitudes greater than or equal to the first threshold as strong edge pixels and assigning them a value of 1; assigning pixels with gradient magnitudes less than the second threshold a value of 0; and assigning a value of 1 to pixels that are less than the first threshold but greater than or equal to the second threshold, if the pixel is connected to a strong edge pixel, then assigning that pixel a value of 1, otherwise assigning it a value of 0.

[0016] Furthermore, step 5 specifically involves: setting a coherence threshold, using pixels with coherence values ​​higher than the coherence threshold as coherent points, employing the Delaunay triangulation method to generate arc segments connecting each coherent point to generate an initial triangulation, expanding outwards from the edge line obtained in step 4 as the center, and using the expanded edge as a buffer. For any triangulation, if a corner pixel exists within the buffer, the pixel in the triangulation is extracted, and all arc segments connected to that pixel are deleted, thereby obtaining the optimized untangled network.

[0017] Furthermore, step 6 uses the Edgelist method to untangle the entanglement.

[0018] A computer device includes a memory, a processor, and a computer program stored in the memory and capable of running on the processor, wherein the processor executes the computer program to implement the method.

[0019] A computer-readable storage medium storing a computer program that, when executed by a processor, implements the method.

[0020] Beneficial Effects: This invention proposes a network-optimized InSAR large gradient deformation phase unwrapping method based on edge detection. During InSAR phase unwrapping, this method uses edge detection to detect large gradient deformation edges in the interferometric pair, thereby eliminating arc segments that may not satisfy the phase continuity assumption. Based on the optimized network, L... p Norm-based solutions are used to recover the absolute phase of large-gradient deformation. Master-slave SAR images are processed by differential interferometry to obtain a wrapped phase map and a coherence map. The gradient is calculated on the wrapped interferogram, and the original phase gradient map is masked by setting high and low thresholds. Large-gradient deformation edge detection is performed using the Canny operator. Simultaneously, to effectively avoid the influence of low-coherence regions, a coherence threshold is set to extract high-coherence points. An initial unwrapped network is constructed using a Delaunay triangulation, and the initial network is optimized using large-gradient deformation phase edges. Finally, based on the optimized network, L... p The method employs phase unwrapping. Experiments show that the proposed method significantly improves performance in both single interferograms and deformation rates and time series. It has broad applicability to monitoring scenarios with large-gradient deformation characteristics, such as landslides and mining, and can provide important technical support for geological disaster prevention and control. Attached Figure Description

[0021] Figure 1 This is an overall flowchart of the present invention;

[0022] Figure 2 This is a schematic diagram of the 5×5 Gaussian filter of the present invention;

[0023] Figure 3 Figure (a) shows a schematic diagram of the convolution kernel calculated by the Sobel operator for gradient calculation, and Figure (b) shows a schematic diagram of the convolution calculation along the x-axis of the image.

[0024] Figure 4 This is a diagram illustrating the principle of magnitude nonmaximum suppression using the Canny operator based on the gradient direction.

[0025] Figure 5 The following are example images of the Canny method used for interferogram fringe edge detection, where Figure (a) is the interferogram and Figure (b) is the edge map;

[0026] Figure 6 Figure (a) shows the network diagram obtained using the Delaunay triangulation method on the interferometric image, and Figure (b) shows the optimized network diagram.

[0027] Figure 7 The figure shows the results obtained by applying the method of the present invention in Gya La Zong.

[0028] Figure 8 Figure (a) shows the line-of-sight deformation rate of the Jialazong landslide area obtained by the method of the present invention and the existing method, and Figure (b) shows the line-of-sight deformation rate obtained by the existing method.

[0029] Figure 9 Figure 1 shows the time-series deformation analysis diagrams of typical points of the Jialazong landslide obtained by the method of the present invention and existing methods. Figure 2 shows the time-series deformation diagram at point P1 obtained by the existing method, Figure 3 shows the time-series deformation diagram at point P2 obtained by the existing method, Figure 4 shows the time-series deformation diagram at point P1 obtained by the method of the present invention, and Figure 5 shows the time-series deformation diagram at point P2 obtained by the method of the present invention. Detailed Implementation

[0030] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an improper limitation of the invention.

[0031] This invention provides a network-optimized InSAR large gradient deformation phase unwrapping method based on edge detection, the specific process of which is as follows: Figure 1 As shown;

[0032] In this embodiment, the original image is first subjected to differential interferometry processing before the phase unwrapping method is performed to obtain a differential interferogram.

[0033] (1) Edge detection using the Canny operator

[0034] 1) Gaussian filtering

[0035] The presence of noisy phases can lead to the detection of many erroneous edges. Therefore, Gaussian filtering is first applied to the wrapped phase gradient map. A weighted average is then performed using a convolution kernel generated according to the Gaussian formula, based on the neighborhood pixel values ​​within the center pixel window. This effectively filters out Gaussian noise in the image, ensuring a smooth input image for subsequent detection.

[0036] The formula for a two-dimensional Gaussian distribution is:

[0037]

[0038] In the formula, (x, y) are the point coordinates, and σ is the standard deviation, which controls the degree of dispersion; the larger the value, the stronger the filtering effect. In actual processing, a Gaussian filter template, i.e., a convolution kernel, needs to be given, and filtering is achieved by traversing the image pixel by pixel. To determine the Gaussian filter template, the Gaussian function can usually be discretized. Specifically, taking the center position of the template as the origin, the coordinates of each neighboring pixel are substituted into Equation 1 to obtain the coefficients of the Gaussian filter template at the corresponding position. Common Gaussian filters include... Figure 2 As shown. By Figure 2 It can be seen that the Gaussian filter template has a pyramid structure. Its main advantage is that as the distance from the center pixel increases, the weight assigned to the pixel decreases. That is, the closer the pixel is to the center pixel, the greater its contribution to the filtering of the center pixel. The weighted average of the neighboring pixels effectively eliminates noise.

[0039] 2) Calculate the image gradient magnitude and direction

[0040] After filtering the original image, the next step is to calculate the gradient magnitude and direction. First, the Sobel operator is used to calculate the image gradient, with convolution kernels used in the horizontal and vertical directions as follows: Figure 3 As shown, where G x With G y Convolution is performed along the x and y directions of the image respectively, resulting in gradient images d along the two directions. x and d y Then, the magnitude and direction of the gradient are calculated using equations 2 and 3, respectively.

[0041]

[0042] In the formula, d x With d y Let M be the gradient along the x and y directions, respectively, where (x, y) is the pixel position, M is the gradient magnitude, and θ is the gradient direction.

[0043] 3) Gradient image nonmaximum suppression

[0044] The gradient magnitude image obtained from the above steps contains coarse and weak edges. To further refine the edges and make them closer to true edges, while removing weak edges caused by noise points, it is necessary to find the maximum gradient magnitude locally, that is, to suppress non-maximum values. The method used by Canny edge detection uses the gradient direction as an aid to determine whether the current pixel is a maximum value in its region. Its principle is as follows: Figure 4 As shown.

[0045] As shown in the blue box, for pixel P0, its gradient direction has been determined based on the angle θ. Then, the amplitude values ​​of its two neighboring pixels P1 and P2 are taken along this gradient direction, and the amplitude value at pixel P0 is compared to the maximum of the three. Pixel P0 is assigned a value of 1. Similarly, for the red box, the amplitude value of pixel P0' is less than that of its neighboring pixel P1', therefore this pixel cannot be considered an edge point and is assigned a value of 0. This suppresses non-maximum values ​​and retains the point with the largest local gradient, resulting in a refined edge. This completes the initial edge extraction.

[0046] 4) Double threshold edge connection

[0047] The above steps have yielded relatively refined edges, but some false edges still exist. Therefore, to further improve the quality of the extracted edges, the Canny method uses a dual threshold to further filter edge points. Specifically, two thresholds, high and low, are used to filter the edge points initially selected in the previous step. Pixels with gradient magnitude values ​​below the low threshold are set as non-edge pixels with a value of 0; pixels with gradient magnitude values ​​above the high threshold are set as strong edge pixels with a value of 1. Pixels with gradient magnitude values ​​above the low threshold but below the high threshold are further examined: if the pixel is connected to a strong edge pixel, the pixel is set to 1; otherwise, it is set to 0. This completes the main steps of edge detection using the Canny method. For example, consider a real interferogram of a landslide, as shown in the image. Figure 5 As shown in (a), edge detection is performed using the Canny method described above, and the resulting edges are as follows: Figure 5 (b) As shown by the white line, the results show that the Canny method can effectively detect the edges of the interferogram, especially for areas with densely varied stripes and large deformation gradients. Based on this method, relatively reliable and fine edge information can be obtained.

[0048] (2) Network optimization based on edge constraints

[0049] Regarding the generation of the initial unwrapped network, this invention generates a coherence map from the interferogram and extracts coherent points (coherent points are pixels with coherence values ​​higher than the coherence threshold). Using the Delaunay triangulation method, arc segments are generated to connect each coherent point, forming a triangular network as shown below. Figure 6 As shown in (a), arc segment removal across large gradient deformation edges is achieved based on edge constraints. First, after extracting the edge line, a buffer zone is established along both sides of the edge line by setting a certain number of pixels. That is, the edge region is expanded based on the set number of pixels. For example, if the number of pixels is set to 2, the edge line expands outward by two pixels. Figure 6 (a) In the initial triangular network constructed, any cell whose corner point falls within the buffer is extracted, and all arc segments connected to these cells are deleted, thereby achieving the removal of large gradient arc segments. The optimized network is as follows: Figure 6As shown in (b), it can be seen from the figure that arc segments crossing the edge line were effectively eliminated.

[0050] (3) Based on optimized network phase unwrapping

[0051] After optimizing the network, this embodiment uses the Edgelist method for untangling.

[0052] The specific form of its optimization model is shown in the following equation:

[0053] min ∑f(k ij )

[0054]

[0055] In the formula, the variable k is constructed. ij The difference between the fuzzy numbers at the two endpoints i and j of the arc segment is an integer, n. j and n i Let be the fuzzy number of the absolute phase at both endpoints of the arc segment, and let be an integer. and These are the winding phases at pixels i and j, respectively. Based on this model, two additional variables p are introduced. ij and q ij , and k ij The relationship is:

[0056] k ij =p ij -q ij

[0057] The optimization model can then be expressed as:

[0058] min ∑c ij (p ij +q ij )

[0059]

[0060] In the formula, p ij With q ij All are non-negative real numbers, c ij The cost function is used, and the transformed model also uses L. 1 The norm is given, and the constraint coefficient matrix is ​​a fully monomodular matrix; therefore, it can be solved using linear programming. The advantage of the Edgelist method is that it does not require closed triangles to calculate residuals, but solves directly on the arc segments; furthermore, it can directly solve for the absolute phase ambiguity number n at each pixel. j and n iThis eliminates the need to perform integration after calculating the fuzzy gradient on the arc segment. This method is well-suited for untangling triangular meshes that are not closed due to arc segment removal after network optimization.

[0061] Figure 7 Images obtained by phase unwrapping of the Xalazonghua area using the method of this invention at different time periods are shown. Each column in the image corresponds to an interference pair. Figure 7 Figures (a), (c), (e), and (g) in the figure are unwrapped diagrams obtained by existing methods. Figure 7 Figures (b), (d), (f), and (h) in the figure are the unwrapping diagrams obtained by the method proposed in this invention. It can be seen that the proposed unwrapping method can effectively eliminate most of the phase jumps in the unwrapping of the landslide area and successfully recover the large gradient deformation phase of the landslide.

[0062] Figure 8 To obtain a line-of-sight deformation rate map of the Jialazong landslide area using existing methods and the method of this application, by Figure 8 It can be seen that this application has eliminated large gradient arcs.

[0063] Figure 9 To illustrate the deformation time series diagrams of typical points in the Jialazong landslide area obtained using existing methods and the method of this application, by... Figure 9 It is evident that the unwrapping method of the present invention can effectively reduce or even eliminate the impact of unwrapping errors in time-series InSAR analysis compared to existing methods.

[0064] It should also be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. An edge detection-based network-optimized InSAR large-gradient deformation phase unwrapping method, characterized in that, Specifically comprising the following steps: Step 1: differential interference processing is performed on the original image to obtain a differential interference image; Step 2: Gaussian filtering is performed on the differential interference image, and the gradient amplitude and gradient direction of each pixel in the Gaussian filtered image are calculated; a gradient amplitude image is obtained; Step 3: non-maximum suppression is performed on the gradient amplitude image, and preliminary extraction of the edge is completed; Step 4: a double threshold is set to further screen the edge points, and finally the extraction of the edge is completed; Step 5: a coherence map is generated based on the differential interference image, and coherent points are extracted, and an unwinding network is generated based on the extracted coherent points and the edge extracted in step 4; Step 6: phase unwinding is performed based on the unwinding network.

2. The network optimization InSAR large gradient deformation phase unwrapping method based on edge detection according to claim 1, characterized in that, The step 2 specifically calculates the gradient amplitude and gradient direction of each pixel by using the following formula: where (x, y) is the coordinate of the pixel, M(x, y) represents the gradient magnitude of the pixel, θ(x, y) is the gradient direction of the pixel, d x (x, y) represents the gradient value of the image along the x-axis at the pixel position (x, y), d y (x, y) represents the gradient value of the image along the y-axis at the pixel position (x, y).

3. The network optimization InSAR large gradient deformation phase unwrapping method based on edge detection according to claim 1, characterized in that, The step 3 is specifically: for any one pixel P, along the gradient direction of the pixel, two adjacent pixels of the pixel are obtained, the gradient amplitudes of the pixel P and its two adjacent pixels are compared, if the gradient amplitude of the pixel P is greater than the gradient amplitudes of the two adjacent pixels, the pixel P is assigned a value of 1, otherwise the pixel P is assigned a value of 0.

4. The network optimization InSAR large gradient deformation phase unwrapping method based on edge detection according to claim 1, characterized in that, The step 4 is specifically: a first threshold and a second threshold are preset, and the first threshold is greater than the second threshold, the pixel with a gradient amplitude greater than or equal to the first threshold is regarded as a strong edge pixel and is assigned a value of 1, the pixel with a gradient amplitude less than the second threshold is assigned a value of 0, for the pixel with a gradient amplitude less than the first threshold and greater than or equal to the second threshold, if the pixel is connected with the strong edge pixel, the pixel is assigned a value of 1, otherwise the pixel is assigned a value of 0.

5. The network optimization InSAR large gradient deformation phase unwrapping method based on edge detection according to claim 1, characterized in that, The step 5 is specifically: a coherence threshold is set, the pixel with a coherence value higher than the coherence threshold is regarded as a coherent point, a Delaunay triangular network method is used to generate an initial triangular network connecting each coherent point, the edge obtained in step 4 is expanded to both sides as a buffer zone, for any one triangular network, if there is a corner pixel in the buffer zone, the pixels in the triangular network are extracted, and all the arc segments connected with the pixel are deleted, so as to obtain an optimized unwinding network.

6. The network-optimized InSAR large-gradient deformation phase-unwrapping method based on edge detection according to claim 5, characterized in that, Step 6 uses the Edgelist method for unwinding.

7. A computer device comprising a memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that, The processor executes the computer program to realize the method of any one of claims 1 to 6.

8. A computer-readable storage medium storing a computer program, the computer-readable storage medium comprising instructions that, when executed by a computer, cause the computer to perform the method of any one of claims 1 to 7. The computer program is executed by the processor to realize the method of any one of claims 1 to 6.

Images