Navigation star bistatic InSAR phase unwrapping method based on multi-frequency-point combined triangulation network

Abstract

The invention discloses a navigational star bistatic InSAR phase unwrapping method based on a multi-frequency-point combined triangulation network, and the method comprises the steps: building a Delaunay triangulation network according to the position of a PS point on a main frequency point and a corresponding phase; optimizing the Delaunay triangulation network in the previous step on the basis of the PS inter-point distance; constructing a multi-frequency-point joint triangulation network on the basis of an original triangulation network based on the position and the original phase of the multi-frequency-point PS point; calculating residual errors in triangles in the combined triangulation network, and establishing a dual network; calculating the flow on the dual network through a minimum cost flow calculation method; and determining a branch cutting path according to the flow on the dual network, and carrying out integration along the branch cutting path to obtain an unwrapping phase on the PS point.

Description

Technical Field

[0001] The present invention belongs to the technical field of bistatic synthetic aperture radar, and in particular relates to a navigation satellite bistatic InSAR phase unwrapping method based on a multi-frequency point joint triangulation network. Background Art

[0002] Global Navigation Satellite System-based Bistatic Interferometric SAR (GNSS-based InBSAR) utilizes in-orbit navigation satellites as transmitters and ground-based receivers. This system employs Persistent Scatter InSAR (PS-InSAR) technology, and the receivers can be airborne, vehicle-mounted, or even fixed. Due to the abundance of in-orbit navigation satellites, the lack of a transmitter, and the use of Beidou satellite L-band signals, the Navigation Satellite Bistatic InSAR system can measure millimeter-level deformations across the observation area at a low hardware cost. It is a radar remote sensing technology widely used for landslide early warning in areas prone to geological hazards.

[0003] However, for areas experiencing large deformation, the deformation phase in the received scene reflection signal will fold due to the large deformation, affecting the accuracy of deformation inversion. At the same time, when using a MEO satellite with a reorbit period of 7 days for observation, compared with an IGSO satellite with a reorbit period of 1 day, the deformation phase will be more prone to phase folding due to the increase in accumulation time, which poses a challenge to deformation measurement in areas prone to geological disasters. As a result, the deformation measurement accuracy of the navigation satellite bistatic InSAR system in large deformation areas is seriously restricted. In addition, the number of PS points that can be extracted by the navigation satellite bistatic InSAR system is smaller than that of traditional ground-based InSAR systems, making it more difficult to calculate the true phase. Therefore, a phase unwrapping method suitable for the navigation satellite bistatic InSAR system is urgently needed to perform phase unwrapping on the differential phase extracted by the navigation satellite bistatic InSAR system in large deformation scenes, so as to obtain accurate and reliable deformation inversion results. Summary of the Invention

[0004] In view of this, the present invention provides a navigation satellite bistatic InSAR phase unwrapping method based on a multi-frequency point joint triangulation network. The technical solution of this method is:

[0005] The phase unwrapping method of bistatic InSAR based on multi-frequency point joint triangulation network includes:

[0006] Step 1: Establish a Delaunay triangulation based on the position and corresponding phase of the PS point on the main frequency point;

[0007] Step 2: Based on the distance between PS points, optimize the Delaunay triangulation in the previous step;

[0008] Step 3: Based on the positions and original phases of the multi-frequency PS points, a multi-frequency joint triangulation network is constructed on the basis of the original triangulation network.

[0009] Step 4: Calculate the residuals of each triangle in the joint triangulation network and establish a dual network;

[0010] Step 5: Calculate the flow on the dual network using the minimum cost flow solution method;

[0011] Step 6: Determine the branch cutting path according to the flow on the dual network, integrate along the branch cutting path, and obtain the unwrapping phase at the PS point.

[0012] Beneficial effects:

[0013] This paper provides a phase unwrapping method for bistatic InSAR (InSAR) using a multi-frequency joint triangulation network. Based on Delaunay triangulation theory and the principle of minimum cost flow, this method fully utilizes the effective information at high-precision PS points across multiple frequencies to perform phase unwrapping between discrete points, thereby obtaining reliable deformation inversion results. First, this method overcomes the disadvantages of traditional phase unwrapping methods for unwrapping two-dimensional images, enabling phase unwrapping even when the number of PS points is small and sparsely distributed. Second, this method introduces joint processing of PS points across multiple frequencies, increasing the input of effective information. Third, this method utilizes the minimum cost flow method for network flow analysis, reducing computational complexity and improving runtime. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 , is a flow chart of the algorithm of the present invention;

[0015] Figure 2 , schematic diagram for constructing a multi-frequency point triangulation network

[0016] Figure 3 , is a schematic diagram of Delaunary triangulation and dual network;

[0017] Figure 4 , are the imaging results in the examples and the multi-frequency point triangulation results constructed using the method proposed in the present invention;

[0018] Figure 5 , is the PS point phase result of the embodiment without using the method proposed by the present invention;

[0019] Figure 6, is the PS point phase unwrapping result using the method proposed in the present invention in the embodiment;

[0020] Figure 7 , are the cumulative deformation results of the embodiments using the present invention. DETAILED DESCRIPTION

[0021] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0022] The present invention is a phase unwrapping method for a navigation satellite bistatic InSAR based on a multi-frequency joint triangulation network. First, a Delaunay triangulation network is established based on the positions and corresponding phases of PS points on the main frequency points. Then, the Delaunay triangulation network from the previous step is optimized based on the distances between PS points. Next, a multi-frequency joint triangulation network is constructed based on the original triangulation network based on the positions and original phases of the multi-frequency PS points. The residuals within each triangle within the joint triangulation network are calculated to establish a dual network. Next, the flow on the dual network is calculated using a minimum cost flow solution. Finally, a branch cut path is determined based on the flow on the dual network, and integration is performed along the branch cut path to obtain the unwrapped phase at the PS point.

[0023] The algorithm flow chart of the present invention is as follows Figure 1 As shown, the specific steps are:

[0024] Step 1: Establish a Delaunay triangulation based on the position and corresponding phase of the PS point on the main frequency point; specifically:

[0025] A Delaunay triangulation is a set of connected but non-overlapping triangles whose circumscribed circles do not contain any other points in the domain. It satisfies the following two conditions:

[0026] 1) There is no circle between any two vertices that contains other points;

[0027] 2) The circumcircle of a triangle formed by any three vertices does not contain any other points.

[0028] Triangulation algorithms can be divided into three categories: divide-and-conquer, point-by-point interpolation, and triangulation growth. Since the phase unwrapped original image typically occupies a large amount of space, a low-memory algorithm should be selected. Therefore, this method uses point-by-point interpolation to construct the triangulation. Let the selected PS point set be V, and a Delaunay triangulation N = (V, E) be established, where E is the edge set. The Delaunay triangulation algorithm can be divided into two main steps:

[0029] 1) Constructing a super triangle: To ensure that all given points are inside the triangulation network, you need to first construct a super triangle that contains all the points. The vertex coordinates of the super triangle must be large enough to ensure that all points are inside it.

[0030] 2) Point-by-point insertion: Insert the given points into the current triangulation one by one. Each time a point is inserted, the topology of the triangulation is updated to meet the Delaunay condition.

[0031] Step 2: Based on the distance between PS points, optimize the Delaunay triangulation in the previous step; specifically:

[0032] When constructing a triangulated network, the connected points should be continuous in space. However, due to the differences in surface scattering characteristics at different locations in the image, it is necessary to limit the length of the triangulated network edges and remove excessively long edges. Calculate the length of each edge in the triangulated network N = (V, E) and set the distance threshold of the triangulated network edge length according to the scene characteristics. men , for example, dis men = 150m. Based on this, remove the edges in the triangulated network N that are greater than the distance threshold, then check the connectivity of each point in the triangulated network, delete the disconnected points, and obtain the optimized triangulated network N ad =(V con ,E men ), where V con is the point set after removing disconnected PS points, E men It is the edge set of the triangulated network after being selected by the distance threshold.

[0033] Step 3: Based on the positions and original phases of the multi-frequency PS points, a multi-frequency joint triangulation network is constructed on the basis of the original triangulation network. Specifically:

[0034] Since the interferograms of the same satellite at multiple frequency points have the same configuration, they are theoretically registered images. Therefore, based on the triangulated network generated by the PS points on the main frequency band, the PS points on the other two frequency points are used to refine the triangulated network. A dense triangulated network means more nodes and edges, which can more strictly enforce the consistency of the phase of multiple images in the overlapping area through the flow constraints on the edges. Let the multi-frequency PS point set be V fu , and its corresponding original phase is First, the phase of the PS point at other frequencies Switch to the phase at the main frequency It can be expressed as:

[0035]

[0036] Among them, λ fu is the wavelength at multiple frequency points, and λ is the wavelength at the main frequency point.

[0037] Then, traverse the multi-frequency PS point set V fu Every point p∈V fu , can be divided into triangulated internal points p in and external point p out Two situations, such as Figure 2 As shown. For the triangulated network N ad External multi-frequency point PS point p out , directly combined with p out and V con Build extended triangulation N extend =(V con +p out ,E extend ). For the triangulated network N ad The nth triangle T n (i,j,k),i,j,k∈V con Multiple frequency points PS point p T ∈p in , based on point set Build a triangulated subset Then merge into the extended triangulation N extend , the final multi-frequency point joint triangulation network can be expressed as:

[0038]

[0039] Step 4: Calculate the residuals of each triangle in the joint triangulation network and establish a dual network; specifically:

[0040] Calculate each triangle Tri in the triangulation network i The residual value is expressed as:

[0041]

[0042] Where m, j, k represent the three vertices in the triangle. and is the winding phase gradient of the three edges, and round is the rounding operator.

[0043] Next, we build N ass The dual network N * =(V * ,E * ). First, we introduce the concept of dual network: for a given planar graph N = (V, E), it has faces (i.e. closed polygons) F1, F2, ..., F n , if there is N * =(V * ,E * ) meets the following conditions:

[0044] 1) For any face F of graph N i , N * There is only one node inside

[0045] 2) For face F of figure N i ,F j The common boundary k , N * There exists an edge in Make and With e k intersect.

[0046] Then it is called graph N * is the dual graph of graph N. The schematic diagram of triangulated network and dual network is as follows Figure 3 As shown in the figure, black nodes represent extracted high-quality PS points, connected by solid lines to form a triangulated network. The dashed lines in the figure form the dual network, and the white nodes represent residual points. The values ​​of the residual points are obtained by using the solid triangles in the triangulated network and can be +1, 0, or -1. For any triangle in the solid triangulated network, there is exactly one node inside it in the dashed network; for any edge of a triangle in the triangulated network, there is exactly one edge in the dashed network that intersects it. Therefore, the dashed network is called the dual network of the triangulated network.

[0047] In addition, it is necessary to construct grounding nodes. According to the points where the residual values ​​are not zero, super source nodes and super sink nodes are constructed respectively. Their residual values ​​are -z and f respectively, where z and f represent the number of positive and negative residual points.

[0048] Step 5: Calculate the flow on the dual network using the minimum cost flow solution method; specifically:

[0049] In the dual network N * In the process, the minimum cost flow method is applied to connect the positive and negative residual point pairs and calculate the minimum cost flow set. The calculated flow is the branch cut with direction, and then the phase matrix can be integrated according to the size and direction of the flow to obtain the disentanglement result. In the dual network N * In the mathematical description of the minimum cost flow problem, the following is the mathematical description:

[0050]

[0051] The constraints are:

[0052]

[0053] 0≤x ij ≤u ij ,(i,j)∈E * (6)

[0054] where c ij For the edge The cost on x ij For the edge Traffic on u ij For the edge The capacity on.

[0055] Currently, the main algorithms for solving the minimum cost flow problem include the original dual algorithm, the flawed algorithm, the relaxation algorithm, and the network simplex algorithm. The following calculation uses the more efficient relaxation algorithm. The relaxation algorithm can solve multi-source and multi-sink minimization problems, and the minimization problem in phase unwrapping can be directly solved using this algorithm.

[0056] Step 6: Determine the branch-cut path based on the flow on the dual network, integrate along the branch-cut path, and obtain the unwrapping phase at the PS point; specifically:

[0057] After these steps, the minimum cost flow x can be calculated, and the branch cut can be directly obtained from the flow x. Based on the branch cut location, the wrapped phase is solved in the Delaunay triangulation by integrating the branch cut with the addition or subtraction method of 2xπ. This provides the unwrapped phase results for the navigation satellite bistatic InSAR system.

[0058] Table 1 Parameters of the embodiment

[0059]

[0060] The following is an explanation of the processing results of the embodiment. In this embodiment, according to the parameters in Table 1, Longxigou in Chongqing is used as the experimental scene, which is a natural landslide. The experiment was carried out using the echo data of Beidou MEO satellite data from August 6, 2024 to September 3, 2024. Taking PRN22 on August 14, 2024 as an example, the imaging results and the multi-frequency point triangulation results constructed using the method proposed in the present invention are as follows: Figure 4 The original phase and unwrapped phase results at PS point are shown as follows: Figure 5 and Figure 6 As shown, the method proposed in the present invention can effectively untangle the jumping wrapped phase and obtain a continuous untangle phase result.

[0061] According to the unwrapped differential phase, a three-dimensional deformation inversion is performed, and the deformation in the east, north, and sky directions is projected into a one-dimensional deformation in the radar receiver line of sight. The cumulative deformation results from August 6, 2024 to September 3, 2024 are obtained as follows: Figure 7 As shown, the effectiveness of the method proposed in this invention is verified, which can overcome the phase folding problem caused by large deformation and long re-orbit time.

[0062] The above specific embodiments merely illustrate the design principles of the present invention. The shapes and names of the components described herein may vary and are not limiting. Therefore, those skilled in the art may modify or substitute equivalents for the technical solutions described in the above embodiments. Such modifications and substitutions, without departing from the inventive spirit and technical solutions of the present invention, shall fall within the scope of protection of the present invention.

Claims

1. A phase unwrapping method for bistatic InSAR based on a multi-frequency joint triangulation network, including: Step 1: Establish a Delaunay triangulation based on the position and corresponding phase of the PS point on the main frequency point; Step 2: Based on the distance between PS points, optimize the Delaunay triangulation in the previous step; Step 3: Based on the positions and original phases of the multi-frequency PS points, a multi-frequency joint triangulation network is constructed on the basis of the original triangulation network. Step 4: Calculate the residuals of each triangle in the joint triangulation network and establish a dual network; Step 5: Calculate the flow on the dual network using the minimum cost flow solution method; Step 6: Determine the branch cutting path according to the flow on the dual network, integrate along the branch cutting path, and obtain the unwrapping phase at the PS point.

2. The method for phase unwrapping of a navigation satellite bistatic InSAR based on a multi-frequency point joint triangulation network according to claim 1, characterized in that: In step 1, the triangulated network is constructed using the point-by-point insertion method.

3. The method for phase unwrapping of a navigation satellite bistatic InSAR based on a multi-frequency point joint triangulation network according to claim 1, characterized in that: In step 1, the Delaunay triangulation generation algorithm can be divided into two main steps: 1) Constructing a super triangle: To ensure that all given points are inside the triangulation network, you need to first construct a super triangle that contains all points; the vertex coordinates of the super triangle must be large enough to ensure that all points are inside it; 2) Point-by-point insertion: insert the given points into the current triangulation network one by one; Each time a point is inserted, the topology of the triangulation network is updated to satisfy the Delaunay condition.

4. The method for phase unwrapping of a navigation satellite bistatic InSAR based on a multi-frequency point joint triangulation network according to claim 1, wherein: In step 2, the length of each edge in the triangulated network N = (V, E) is calculated, and the distance threshold of the triangulated network edge length is set according to the scene characteristics. men ; Based on this, remove the edges in the triangulated network N that are greater than the distance threshold, then check the connectivity of each point in the triangulated network, delete the disconnected points, and obtain the optimized triangulated network N ad =(V con ,E men ), where V con is the point set after removing disconnected PS points, E men It is the edge set of the triangulated network after being selected by the distance threshold.

5. The method for phase unwrapping of a navigation satellite bistatic InSAR based on a multi-frequency point joint triangulation network according to claim 1, wherein: In step 3, the phase of the PS point at other frequencies is Switch to the phase at the main frequency Expressed as: Among them, λ fu is the wavelength at multiple frequency points, and λ is the wavelength at the main frequency point.

6. The method for phase unwrapping of a navigation satellite bistatic InSAR based on a multi-frequency point joint triangulation network according to claim 1, wherein: In step 3, for the triangulated network N ad The nth triangle T n (i,j,k),i,j,k∈V con Multiple frequency points PS point p T ∈p in , based on point set Constructing a triangulated subset Then merge into the extended triangulation N extend , the final multi-frequency point joint triangulation network is expressed as:

7. The method for phase unwrapping of a navigation satellite bistatic InSAR based on a multi-frequency point joint triangulation network according to claim 1, wherein: In step 4, calculate each triangle Tri in the triangulation network i The residual value is expressed as: Where m, j, k represent the three vertices in the triangle. and is the winding phase gradient of the three edges, and round is the rounding operator.

Images