Navigation satellite dual base inSAR phase unwrapping method based on multi-frequency joint triangulation network

Abstract

This invention discloses a method for unwrapping phases in bistatic InSAR navigation satellites based on a multi-frequency joint triangulation. The method includes: establishing a Delaunay triangulation based on the position and corresponding phase of the PS point on the dominant frequency; optimizing the Delaunay triangulation based on the distance between PS points; constructing a multi-frequency joint triangulation based on the position and original phase of the PS points on multiple frequencies; calculating the residuals within each triangle of the joint triangulation to establish a dual network; calculating the flow on the dual network using a minimum cost flow solution method; determining the branch path based on the flow on the dual network, integrating along the branch path to obtain the unwrapped phase at the PS point.

Description

Technical Field

[0001] This invention belongs to the field of bistatic synthetic aperture radar technology, specifically relating to a navigation satellite bistatic InSAR phase unwrapping method based on a multi-frequency joint triangulation network. Background Technology

[0002] Bistatic Synthetic Aperture Radar Interferometric SAR (GNSS-based InBSAR) is a bistatic SAR system that uses in-orbit navigation satellites as transmitters and receivers located near the ground. It employs Persistent Scatters-InSAR (PS-InSAR) technology, and the receivers can be airborne, vehicle-mounted, or even fixed. Due to the abundance of in-orbit navigation satellite resources, the lack of a transmitter requirement, and the ability to utilize BeiDou L-band signals, GNSS-based bistatic InSAR systems can achieve millimeter-level deformation measurements of the observation area with relatively low hardware costs. It is a widely used radar remote sensing technology for landslide early warning in areas prone to geological disasters.

[0003] However, in areas with significant deformation, the deformation phase in the received scene reflection signal may fold due to the large deformation, affecting the accuracy of deformation inversion. Furthermore, when using MEO satellites with a 7-day heavy orbit period (MOP) compared to IGSO satellites with a 1-day MOP, the deformation phase is more prone to phase folding due to the increased cumulative time, posing a challenge to deformation measurement in geologically hazardous areas and severely limiting the deformation measurement accuracy of the navigation satellite bistatic InSAR system in areas of large deformation. In addition, the navigation satellite bistatic InSAR system can extract fewer PS points compared to traditional ground-based InSAR systems, making it more difficult to calculate the true phase. Therefore, a phase unwrapping method suitable for navigation satellite bistatic InSAR systems is urgently needed to unwrap the differential phase extracted by the navigation satellite bistatic InSAR system in large deformation scenarios, thereby obtaining accurate and reliable deformation inversion results. Summary of the Invention

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

[0005] A method for phase unwrapping of bistatic InSAR for navigation satellites based on multi-frequency joint triangulation includes:

[0006] Step 1: Based on the position of the PS point and its corresponding phase at the main frequency point, establish a Delaunay triangulation network;

[0007] Step 2: Optimize the Delaunay triangulation from the previous step based on the distance between PS points;

[0008] Step 3: Based on the location and original phase of the multi-frequency PS points, construct a multi-frequency joint triangulation network on the basis of the original triangulation network;

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

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

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

[0012] Beneficial effects:

[0013] This invention provides a phase unwrapping method for bistatic InSAR navigation satellites based on multi-frequency joint triangulation. Based on Delaunay triangulation theory and the minimum cost flow principle, it fully utilizes the effective information on high-precision PS points at multiple frequencies to perform phase unwrapping between discrete points, thereby obtaining reliable deformation inversion results. First, this method overcomes the disadvantage of traditional phase unwrapping methods in unwrapping two-dimensional images, enabling phase unwrapping even with a small number and sparse distribution of PS points. Second, this method introduces joint processing of PS points at multiple frequencies, increasing the input of effective information. Third, this method uses the minimum cost flow method for network flow solving, which reduces computational load and improves runtime. Attached Figure Description

[0014] Figure 1 The above is a flowchart of the algorithm of this invention;

[0015] Figure 2 A schematic diagram of multi-frequency triangular network construction.

[0016] Figure 3 1. Schematic diagram of Delaunary triangular mesh and dual mesh;

[0017] Figure 4 The images in the examples and the multi-frequency triangulation results constructed using the method proposed in this invention are:

[0018] Figure 5 The given embodiments are PS point phase results without using the method proposed in this invention;

[0019] Figure 6The PS point phase unwrapping results are shown in the examples using the method proposed in this invention.

[0020] Figure 7 The cumulative deformation results of the present invention are used in the examples described. Detailed Implementation

[0021] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0022] This invention is a phase unwrapping method for bistatic InSAR navigation satellites based on a multi-frequency joint triangulation. First, a Delaunay triangulation is established based on the position and corresponding phase of the PS point at the dominant frequency. Then, the Delaunay triangulation from the previous step is optimized based on the distance between PS points. Next, a multi-frequency joint triangulation is constructed based on the position and original phase of the PS points at multiple frequencies. Then, the residuals within each triangle of the joint triangulation are calculated to establish a dual network. Next, the flow on the dual network is calculated using a minimum cost flow solution method. Finally, the branch paths are determined based on the flow on the dual network, and integration is performed along the branch paths to obtain the unwrapped phase at the PS point.

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

[0024] Step 1: Based on the position and corresponding phase of the PS point at the dominant frequency, construct a Delaunay triangulation; specifically:

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

[0026] 1) There are no other points inside the circle between any two vertices;

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

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

[0029] 1) Constructing a super triangle: To ensure that all given points are inside the triangulation, a super triangle containing all points needs to be constructed first. 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: Inserts the given points one by one into the current triangulation. Each time a point is inserted, the topology of the triangulation is updated to satisfy the Delaunay condition.

[0031] Step 2: Optimize the Delaunay triangulation from the previous step based on the distance between PS points; specifically:

[0032] When constructing a triangulation mesh, the connected points should be spatially continuous. However, due to differences in surface scattering characteristics at different locations in the image, it is necessary to limit the side lengths of the triangulation mesh and remove excessively long edges. The length of each edge in the triangulation mesh N = (V, E) is calculated, and a distance threshold dis for the edge length is set based on the scene characteristics. men For example, dis men =150m. Based on this, edges in the triangulation network N that are greater than the distance threshold are removed. Then, the connectivity of each point in the triangulation network is checked, and disconnected points are deleted to obtain the optimized triangulation network N. ad =(V con E men ), where V con To remove disconnected PS points, E men It is the set of edges of the triangulation network after the distance threshold selection.

[0033] Step 3: Based on the location and original phase of the multi-frequency PS points, construct a multi-frequency joint triangulation network on the basis of the original triangulation network; specifically:

[0034] Since the interferograms of the same satellite at multiple frequencies have the same configuration, they are theoretically already registered images. Therefore, based on the triangulation generated from the PS points on the main frequency band, the PS points on the other two frequencies are used to refine the triangulation. A denser triangulation means more nodes and edges, which can more strictly enforce the consistency of phase in overlapping regions of the multiple images through flow constraints on the edges. Let the set of PS points at multiple frequencies be V. fu Its corresponding original phase is First, the phase of the PS point on other frequency points. Phase shifted to the main frequency point It can be represented as:

[0035]

[0036] Where, λ fu λ represents the wavelength at multiple frequency points, while λ represents the wavelength at the dominant frequency point.

[0037] Then, iterate through the multi-frequency point PS set V. fu For each point p∈V fu It can be divided into points p inside the triangular network. in and external point p out Two situations, such as Figure 2 As shown. For the triangular network N... ad External multi-frequency PS point p out Directly combine p out and V con Construct an extended triangulation network N extend =(V con +p out E extend ). For the triangular network N ad The nth triangle T n (i,j,k),i,j,k∈V con Multi-frequency PS point p within T ∈p in Based on point sets Constructing a triangular mesh subset Then merged into the extended triangular network N extend In the meantime, the final multi-frequency joint triangular network can be represented as:

[0038]

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

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

[0041]

[0042] Where m, j, k represent the three vertices of the triangle. and Let be the entanglement phase gradient of the three sides, and round be the rounding operator.

[0043] Below, we establish N. ass dual network N * =(V * E * First, let's introduce the concept of dual networks: For a given planar graph N = (V, E), it has surfaces (i.e., closed polygons) F1, F2, ..., F... n If there exists N * =(V * E * The following conditions must be met:

[0044] 1) For any face F of graph N i N * It has one and only one node inside.

[0045] 2) For surface F of graph N i ,F j public boundary e k N * There exists and there exists an edge. Make and With e k intersect.

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

[0047] In addition, grounding nodes need to be constructed. Based on the points with non-zero residual values, super source nodes and super sink nodes are constructed respectively, with residual values ​​of -z and f, where z and f represent the number of positive and negative residual points, respectively.

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

[0049] In the dual network N * In the dual network N, the minimum cost flow method is applied to connect the pairs of positive and negative residual points, and the minimum cost flow set is calculated. The calculated flow is essentially a directional branch cut, and the phase matrix can be integrated based on the magnitude and direction of the flow to obtain the untangled result. * In this context, the mathematical description of the minimum cost flow problem is as follows:

[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 it.

[0055] Currently, algorithms for solving the minimum cost flow problem mainly include the original dual algorithm, the flawed algorithm, the relaxation algorithm, and the network simplex algorithm. The relaxation algorithm, which is relatively efficient, will be used for computation below. The relaxation algorithm is a type of algorithm that can solve multi-source, multi-sink minimization problems; the minimization problem in phase untangling can be directly solved using this algorithm.

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

[0057] After the above steps, the minimum cost flow *x* can be calculated, and the branch cut can be directly obtained from the flow *x*. Based on the location of the branch cut, the entangled phase is solved in the Delaunay triangulation using the method of adding and subtracting 2xπ through the branch cut. Thus, the unentangled phase result of the navigation satellite bistatic InSAR system is obtained.

[0058] Table 1 Parameter Table for Embodiments

[0059]

[0060] The processing results of the embodiment are described below. In this embodiment, based on the parameters in Table 1, Longxigou in Chongqing was used as the experimental scenario, which is a natural landslide. The experiment was conducted using echo data from BeiDou MEO satellites from August 6, 2024 to September 3, 2024. Taking PRN22 on August 14, 2024 as an example, the imaging results and the results of the multi-frequency triangulation network constructed using the method proposed in this invention are as follows: Figure 4 As shown. The original phase and unwrapped phase results at point PS are respectively as follows: Figure 5 and Figure 6 As shown, the method proposed in this invention can effectively untangle the abruptly entangled phase and obtain continuous untangled phase results.

[0061] Based on the unwrapped differential phase, a three-dimensional deformation inversion is performed. The deformation in the east, north, and sky directions is projected onto the one-dimensional deformation along the radar receiver's line of sight, yielding the cumulative deformation results from August 6, 2024 to September 3, 2024, as shown below. 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-orbiting time.

[0062] The specific embodiments described above only illustrate the design principles of the present invention. The shapes and names of the components in this description may differ and are not limited. Therefore, those skilled in the art can modify or make equivalent substitutions to the technical solutions described in the foregoing embodiments; and these modifications and substitutions do not depart from the inventive spirit and technical solutions of the present invention, and should all fall within the protection scope of the present invention.

Claims

1. A method for phase unwrapping of bistatic InSAR for navigation satellites based on multi-frequency joint triangulation, including: Step 1: Based on the position of the PS point and its corresponding phase at the main frequency point, establish a Delaunay triangulation network; Step 2: Optimize the Delaunay triangulation from the previous step based on the distance between PS points; Step three, based on the position of multi-frequency PS point and the original phase, the multi-frequency joint triangulation network is constructed on the basis of the original triangulation network: first, the phase of PS point on other frequency is calculated based on the phase of the main frequency point Turn to the phase on the main frequency point , expressed as: ； wherein, is the wavelength at the primary frequency point, is the wavelength at the primary frequency point; Secondly, regarding the triangular network The nth triangle Multi-frequency PS points within Based on point sets Constructing a triangular mesh subset Then merged into the extended triangular network. In the end, the final multi-frequency joint triangular network is represented as: ； Step 4: Calculate the residuals within each triangle of the joint triangulation network and establish the dual network; Step 5: Calculate the traffic on the dual network using the minimum cost flow solution method; Step 6: Determine the branch cutting path based on the flow on the dual network, integrate along the branch cutting path, and obtain the untangled phase at point PS.

2. The multi-finger joint triangulation network based navigation satellite bistatic InSAR phase unwrapping method according to claim 1, characterized in that, In step one, a triangular network is constructed using the point-by-point insertion method.

3. The multi-finger joint triangulation network based navigation satellite bistatic InSAR phase unwrapping method according to claim 1, characterized in that, In step one, 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, a super triangle containing all points needs to be constructed first; the vertex coordinates of the super triangle must be large enough to ensure that all points are inside it; 2) Point-by-point insertion: Inserts the given points one by one into the current triangulation; Each time a point is inserted, the topology of the triangular network is updated to satisfy the Delaunay condition.

4. The multi-facet joint-triangulation network based navigation star bistatic InSAR phase unwrapping method according to claim 1, characterized in that, In step two, the triangular network is calculated. The length of each side in the triangular mesh is set according to the characteristics of the scene, with a distance threshold for the side length. Based on this, triangular meshes are eliminated. Find edges with a distance greater than a threshold, then check the connectivity of each point in the triangulation, delete disconnected points, and obtain the optimized triangulation. ,in To the point set after removing disconnected PS points, It is the set of edges of the triangulation network after the distance threshold selection.

5. The multi-facet joint triangulation network based navigation satellite bistatic InSAR phase unwrapping method according to claim 1, characterized in that, In step four, calculate the value of each triangle in the triangulation network. The residual value is expressed as: ； in Represents the three vertices of a triangle. , and The phase gradient of the three sides is the winding phase gradient. This is the floor operator.

Images