Phase unwrapping method and device based on improved minimum cost flow and quasi-identification

Abstract

The invention relates to the technical field of surface deformation monitoring, and provides a phase unwrapping method and device based on improved minimum cost flow and quasi identification, and the method comprises the steps: obtaining SAR and DEM data of a to-be-monitored region, and generating a differential interferogram and a coherence graph; calculating a time baseline phase gradient rate set based on the differential interferogram and the coherence map; constructing a Delaunay triangulation network, and determining the weight of each triangular arc section according to the coherence graph and the time baseline phase gradient rate set; secondly, completing initial phase unwrapping of the to-be-monitored area by using a minimum cost flow method and the weights; and finally, performing quasi-identification on the unwrapping phase of each time baseline, and adjusting an initial phase unwrapping result according to a quasi-identification result to obtain a target phase unwrapping result. According to the embodiment of the invention, the improved minimum cost flow algorithm and the quasi-identification technology are combined, so that the precision and reliability of phase unwrapping are improved, and the high-precision deformation monitoring of the earth surface is realized.

Description

Technical Field

[0001] The present disclosure relates to the technical field of surface deformation monitoring, and in particular to a phase unwrapping method and device based on improved minimum cost flow and quasi-identification. Background Art

[0002] Remote sensing technology based on synthetic aperture radar (SAR) data has become an important tool for geological disaster monitoring, widely used for monitoring and early warning of natural disasters such as landslides, mine collapses, volcanoes, and earthquakes. Compared with traditional remote sensing technology, SAR technology has the advantages of all-weather and all-day operation, can penetrate clouds and monitor at night, and is particularly suitable for remote sensing monitoring in complex terrain and extreme climates.

[0003] With the development of technology, the application of SAR technology has shifted from traditional disaster response to more accurate deformation monitoring, especially phase-based SAR interferometry (InSAR) technology, which has become one of the core technologies of research due to its advantages in monitoring tiny surface deformation. InSAR can monitor surface deformation with centimeter or even millimeter-level accuracy by calculating the interference phase of two SAR images with different perspectives, making it suitable for long-term monitoring of geological disasters.

[0004] Among related technologies, time-series monitoring methods based on InSAR technology have become mainstream. By selecting SAR images with tight spatial baselines and differentiated temporal baselines, they reduce the incoherence problem caused by large baseline differences, thereby improving the spatial resolution and temporal sampling rate of surface deformation monitoring and providing accurate data for disaster warning and risk assessment. However, although InSAR technology can provide accurate deformation information, its accuracy is often challenged during the phase unwrapping process, especially in areas with large gradient deformation. At the same time, the minimum cost flow unwrapping method, as a commonly used unwrapping technique, has improved, but is still limited by SAR resolution and cannot effectively address the phase recovery problem in areas with large gradient deformation. Summary of the Invention

[0005] The embodiments of the present disclosure at least provide a phase unwrapping method and device based on improved minimum cost flow and quasi-identification. By combining the improved minimum cost flow algorithm with quasi-identification technology, the accuracy and reliability of phase unwrapping are improved, thereby realizing high-precision deformation monitoring of the surface.

[0006] The present disclosure provides a phase unwrapping method based on improved minimum cost flow and quasi-identification, including: Acquiring a SAR data set and a DEM data of the area to be monitored; and generating a differential interferogram and a coherence map based on the SAR data set and the DEM data; calculating a time baseline phase gradient rate set based on the differential interferogram and the coherence map; Constructing a Delaunay triangulation based on the differential interference pattern and the coherence map, and weighting each triangular arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set to obtain a weight corresponding to each triangular arc segment in the Delaunay triangulation; Based on the minimum cost flow method and the weights corresponding to the triangular arc segments in the Delaunay triangulation, an initial phase unwrapping task of the monitored area is completed to obtain an initial phase unwrapping result; wherein the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines; The unwrapped phase corresponding to each time baseline is respectively subjected to a pseudo-accuracy identification, and the initial phase unwrapping result is adjusted based on the pseudo-accuracy identification result to obtain a target phase unwrapping result.

[0007] The embodiment of the present disclosure provides a phase unwrapping device based on improved minimum cost flow and quasi-identification, including: A data acquisition module is used to acquire SAR data sets and DEM data of the area to be monitored; and generate differential interferograms and coherence maps based on the SAR data sets and DEM data; a rate calculation module, configured to calculate a time baseline phase gradient rate set based on the differential interference pattern and the coherence pattern; a weight determination module, configured to construct a Delaunay triangulation based on the differential interference pattern and the coherence map, and to weight each triangular arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set, thereby obtaining a weight corresponding to each triangular arc segment in the Delaunay triangulation; a phase unwrapping module, configured to complete an initial phase unwrapping task of the monitored area based on a minimum cost flow method and weights corresponding to each triangular arc segment in the Delaunay triangulation, and obtain an initial phase unwrapping result; wherein the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines; The phase alignment module is used to perform alignment identification on the unwrapped phase corresponding to each time baseline respectively, and adjust the initial phase unwrapping result based on the alignment identification result to obtain the target phase unwrapping result.

[0008] An embodiment of the present disclosure provides a computer device, including: a processor, a memory, and a bus, wherein the memory stores machine-readable instructions executable by the processor. When the computer device is running, the processor and the memory communicate through the bus, and when the machine-readable instructions are executed by the processor, a phase unwrapping method based on improved minimum cost flow and quasi-identification as described in any possible embodiment described above is performed.

[0009] An embodiment of the present disclosure provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the phase unwrapping method based on improved minimum cost flow and quasi-identification as described in any of the above possible implementations is implemented.

[0010] The phase unwrapping method and apparatus based on improved minimum cost flow and quasi-identification provided in the embodiments of the present disclosure first acquire a synthetic aperture radar (SAR) dataset and digital elevation model (DEM) data for the area to be monitored. Subsequently, a differential interferogram and a coherence map are generated based on the SAR dataset and DEM data. Next, the differential interferogram and coherence map are used to calculate a time-base phase gradient rate set, and based on this, a Delaunay triangulation is constructed. Based on the coherence map and the time-base phase gradient rate set, each triangular arc segment in the Delaunay triangulation is weighted to obtain a weight corresponding to each triangular arc segment. The minimum cost flow method and these weights are then used to complete the initial phase unwrapping task for the area to be monitored, obtaining an initial phase unwrapping result. To further improve the accuracy of phase unwrapping, the present disclosure further proposes performing quasi-identification on the unwrapped phases corresponding to each time baseline. Finally, based on the results of the quasi-identification, the initial phase unwrapping result is adjusted to ultimately obtain the target phase unwrapping result.

[0011] In the disclosed embodiments, an improved minimum cost flow method is used during phase unwrapping to optimize the calculation of phase gradient rate sets and the weighting of triangular arc segments, effectively reducing error accumulation, thereby ensuring high-precision unwrapping results and accelerating the unwrapping process. Furthermore, the disclosed embodiments utilize quasi-identification techniques to post-process the initial phase unwrapping results, further improving phase unwrapping accuracy and enabling high-precision surface deformation monitoring.

[0012] In order to make the above-mentioned objectives, features and advantages of the present disclosure more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the following briefly introduces the drawings that need to be cited in the embodiments. The drawings herein are incorporated into and constitute a part of the specification. These drawings illustrate embodiments consistent with the present disclosure and, together with the specification, are used to illustrate the technical solutions of the present disclosure. It should be understood that the following drawings only illustrate certain embodiments of the present disclosure and should not be regarded as limiting the scope. For those of ordinary skill in the art, other relevant drawings can be obtained based on these drawings without inventive effort.

[0014] Figure 1A flow chart of a phase unwrapping method based on improved minimum cost flow and quasi-identification provided by an embodiment of the present disclosure is shown; Figure 2 A flowchart of a method for determining a differential interference pattern and a coherence pattern provided by an embodiment of the present disclosure is shown; Figure 3 A flow chart of a method for calculating a time baseline phase gradient rate set provided by an embodiment of the present disclosure is shown; Figure 4 A flowchart of a Delaunay triangulation construction method provided by an embodiment of the present disclosure is shown; Figure 5 A flowchart of a method for determining the weight of arc segments of a Delaunay triangulation provided by an embodiment of the present disclosure is shown; Figure 6 A flow chart of an unwrapped phase alignment identification method provided by an embodiment of the present disclosure is shown; Figure 7 A schematic diagram showing a surface deformation field of a simulation experiment provided by an embodiment of the present disclosure is shown; Figure 8 A schematic diagram showing an original interference pattern and a tangled interference pattern generated by a simulation experiment provided by an embodiment of the present disclosure is shown; Figure 9 A schematic diagram showing a phase unwrapping result of a traditional MCF method in a simulation experiment provided by an embodiment of the present disclosure; Figure 10 A schematic diagram showing a phase unwrapping result of an improved MCF method based on phase gradient rate guidance provided by the present disclosure in a simulation experiment provided by an embodiment of the present disclosure is shown; Figure 11 A schematic diagram showing a phase unwrapping result of the unwrapping phase error correction method provided by the present disclosure in a simulation experiment provided by an embodiment of the present disclosure; Figure 12 A schematic structural diagram of a phase unwrapping device based on improved minimum cost flow and quasi-identification provided by an embodiment of the present disclosure is shown; Figure 13 A schematic structural diagram of a computer device provided by an embodiment of the present disclosure is shown. DETAILED DESCRIPTION

[0015] In order to make the purpose, technical solutions and advantages of the embodiments of the present disclosure clearer, the technical solutions in the embodiments of the present disclosure will be clearly and completely described below in conjunction with the drawings in the embodiments of the present disclosure. Obviously, the described embodiments are only part of the embodiments of the present disclosure, not all of the embodiments. The components of the embodiments of the present disclosure generally described and shown in the drawings herein can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present disclosure provided in the drawings is not intended to limit the scope of the disclosure for which protection is sought, but merely represents selected embodiments of the present disclosure. Based on the embodiments of the present disclosure, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of the present disclosure.

[0016] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not require further definition or explanation in subsequent drawings.

[0017] The term "and / or" herein simply describes an association relationship, indicating that three relationships can exist. For example, A and / or B can represent the existence of A alone, the simultaneous existence of A and B, and the existence of B alone. In addition, the term "at least one" herein refers to any combination of at least two of any one or more of a plurality of items. For example, "at least one of A, B, and C" can represent any one or more elements selected from the set consisting of A, B, and C.

[0018] To facilitate understanding of this embodiment, we first provide a detailed introduction to the implementation of the phase unwrapping method based on improved minimum cost flow and quasi-identification provided in this embodiment. The implementation of the phase unwrapping method based on improved minimum cost flow and quasi-identification provided in this embodiment is a computer device. This computer device can be a server. The server can be a standalone physical server, a server cluster or distributed system consisting of multiple physical servers, or a cloud server that provides basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud storage, big data, and artificial intelligence platforms.

[0019] The following describes in detail the phase unwrapping method based on improved minimum cost flow and quasi-identification provided by the embodiment of the present application in conjunction with the accompanying drawings. Figure 1 FIG. 1 is a flow chart of a phase unwrapping method based on improved minimum cost flow and quasi-identification according to an embodiment of the present disclosure. The method includes the following steps S101 to S105: S101, acquiring a SAR data set and a DEM data of a region to be monitored; and generating a differential interferogram and a coherence map based on the SAR data set and the DEM data.

[0020] As you can understand, SAR (Synthetic Aperture Radar) is a technology that uses electromagnetic waves to image the ground. It has all-day, all-weather capabilities and is unaffected by factors such as weather and lighting. SAR transmits radar signals and receives their reflected signals, forming high-resolution surface images. It is widely used in areas such as surface change monitoring and disaster warning. Here, a SAR dataset is image data collected by a SAR sensor that contains electromagnetic wave reflection information from the monitored area, providing accurate information on surface topography and deformation. A DEM (Digital Elevation Model) is a digital model obtained using remote sensing or ground-based measurement techniques that represents the ground elevation of a specific area. Here, DEM data is used to assist in understanding topographic relief and is combined with SAR data through differential interferograms to derive more accurate phase unwrapping results.

[0021] Specifically, in SAR interferometry, a differential interferogram (DIF) is an image of surface deformation obtained by calculating the phase difference between two SAR images. It can show changes in the surface between two points in time and is often used to monitor phenomena such as earthquakes, volcanic activity, and urban subsidence. A coherence map can be used to indicate the degree of similarity between different pixels, or coherence. This map is derived by comparing the phase information of adjacent pixels at different times or orbits, and generally reflects signal quality, noise, and data reliability.

[0022] For example, referring to Figure 2 As shown, determining the differential interferogram and the coherence map based on the SAR data set and the DEM data may include the following steps S201 to S204: S201: Determine an initial SAR image set based on the SAR data set.

[0023] The initial SAR image set consists of multiple SAR images of the area to be monitored. These SAR images typically contain observation data of the same area at different times and angles. The area to be monitored typically refers to a specific area where ground change monitoring is required, such as a high-risk area for earthquakes, landslides, or urban change.

[0024] S202: Determine a SAR main image in the initial SAR image set, use the SAR main image as a reference, register other SAR images in the initial SAR image set except the SAR main image according to the DEM data, and determine a target SAR image set based on the registration results.

[0025] Specifically, the SAR master image refers to the reference image selected from the entire image set. This requires comprehensive consideration of multiple factors, including temporal baseline, spatial baseline, and changes in Doppler center frequency. This means selecting the image closest in time or with the best quality from the orbital dataset as the master image, while ensuring that the selected master image and other images have good matching in terms of spatial baseline and Doppler center frequency. Registration involves aligning other images with the master image so that SAR images from different times and angles can be accurately spatially aligned. During the registration process, DEM data is used to correct the geometric errors of the images, thereby ensuring their spatial consistency. Ultimately, these registered SAR images and the SAR master image are used to generate a target SAR image set for subsequent analysis. The target SAR image set includes multiple target SAR images.

[0026] S203: Determine a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold.

[0027] It can be understood that the spatiotemporal baseline threshold refers to the limitation on the validity between image pairs in time and space, which is usually a parameter value used to screen out representative and relatively consistent interference pairs, and exclude those image pairs whose interference effects are not ideal due to too long time intervals or too large spatial distances. In the present disclosure, the spatiotemporal baseline threshold is set to a time of less than 72 days and a spatial distance of less than 100m; in some other embodiments, it can be set according to specific needs, and no specific limitation is made here. Based on the target SAR image set, the baseline set is determined by a preset spatiotemporal baseline threshold. The baseline set includes multiple interference pairs, and each interference pair includes two target SAR images. Here, the interference pair refers to the pairing of two images required for interference analysis through SAR images. These paired images are usually separated by a certain time and space and can reflect ground changes.

[0028] S204 , generating a differential interference map corresponding to each interference pair based on the DEM data; and generating a coherence map corresponding to each interference pair based on the differential interference map corresponding to each interference pair.

[0029] Specifically, an interferogram pair is formed by two radar observations (or images), each taken at slightly different times and locations. By comparing these two observations, an interferogram can be generated that reflects surface changes (such as subsidence and displacement). A differential interferogram is generated by subtracting topographic errors (usually corrected using DEM data) and can reveal detailed information about relative surface changes. The differential interferogram of each interferogram pair reveals the phase difference between the two observations, which reflects ground displacement or deformation.

[0030] It can be understood that the coherence map reflects the signal quality and correlation of each pixel in the interferogram. Coherence essentially describes the temporal and spatial consistency of the radar signal. In radar imaging, different surface features (such as vegetation, buildings, and mountains) reflect radar signals differently, which affects the signal coherence. Areas with high coherence typically indicate stable surface features and consistent reflections, while areas with low coherence may be due to dynamic surface changes, differences in scatterers, or atmospheric influences. By generating a coherence map based on the differential interferogram of each interferometer pair, it is possible to effectively identify areas with good signal quality, thereby determining which areas have high data confidence and which areas may require further correction or exclusion. For example, if the coherence in a certain area is low, it may indicate significant surface changes or signal distortion caused by atmospheric factors.

[0031] Determining the coherence map corresponding to each interference pair based on the differential interferogram corresponding to each interference pair can be achieved through a coherence calculation algorithm. This algorithm generally involves statistical analysis of the phase information in the differential interferogram to quantify the phase stability of adjacent pixels or the same pixel at different times. Specifically, a coherence formula that calculates the standard deviation of the phase difference or the change in the phase gradient can be used to calculate the coherence value of each pixel, thereby obtaining a coherence map in which the value of each pixel represents the degree of signal coherence at that position in the differential interferogram.

[0032] S102: Calculate a time baseline phase gradient rate set based on the differential interference map and the coherence map.

[0033] It can be understood that the time baseline phase gradient rate set is a collection of phase gradient rates corresponding to different time baselines. Here, "time baseline" refers to the time interval between two SAR images used to generate the differential interferogram. In this disclosure, special attention is paid to three cases: the shortest time baseline, the first time baseline, and the second time baseline. These three cases represent information on surface deformation over different time periods (i.e., time intervals of 12 days, 24 days, and 36 days). By analyzing these different time baselines, the changing characteristics of surface deformation at different time scales can be captured.

[0034] Here, the phase gradient rate reflects the rate of change of the phase difference between adjacent pixels. In the differential interferogram, the phase difference is directly related to the surface deformation, so the phase gradient rate can indirectly reveal the distribution and rate of surface deformation. When determining the time baseline phase gradient rate set, the differential interferogram can be preprocessed first, including steps such as noise removal and filtering, to improve the accuracy and reliability of the phase information; then, using the coherence map as a weighting factor, the phase information in the differential interferogram is weighted to weaken the interference of low coherence areas on the results; secondly, by calculating the phase difference between adjacent pixels and combining it with the time baseline information, the phase gradient rate under different time baselines can be obtained.

[0035] Specifically, the time baseline phase gradient rate set will include phase gradient rates corresponding to different time baselines. In the present disclosure, the time baseline phase gradient rate set includes the shortest time baseline phase gradient rate, the first time baseline phase gradient rate, and the second time baseline phase gradient rate.

[0036] In some other embodiments, the time interval may also include 48 days, 60 days, 72 days, etc., and the corresponding time baselines are expressed as the third time baseline, the fourth time baseline, and the fifth time baseline, etc., which are not specifically limited here.

[0037] For example, referring to Figure 3 As shown, when calculating the time baseline phase gradient rate set based on the differential interference map and the coherence map, the following steps S301 to S304 may be included: S301 : For a differential interferogram corresponding to each interference pair, calculate a phase gradient value corresponding to each pixel in the differential interferogram, and determine a first phase gradient map based on the phase gradient values ​​corresponding to each pixel.

[0038] Specifically, for the differential interferogram corresponding to each interference pair, the phase gradient values ​​of each pixel in the image are determined. These gradient values ​​reflect the rate of change of the phase difference between adjacent pixels. The calculated phase gradient values ​​can be used to determine a first phase gradient map; the first phase gradient map includes a first horizontal phase gradient map, a first vertical phase gradient map, a first diagonal phase gradient map, and a second diagonal phase gradient map.

[0039] Here, the present disclosure uses a phase gradient calculation formula to calculate the phase gradient value corresponding to each pixel, which can be expressed as: ; in, Expressed as the phase gradient in a specific direction of the pixel; Expressed as the vertical phase gradient; It is represented as the phase gradient in the first diagonal direction; Expressed as horizontal phase gradient; It is represented as the second diagonal phase gradient.

[0040] In some other embodiments, other methods may be used to calculate the phase gradient value corresponding to each pixel, such as a method based on frequency domain analysis, a method based on edge detection, etc., which are not specifically limited here.

[0041] S302 : For each interference pair, determine a gradient wrapping result according to a first phase gradient map corresponding to the differential interference map, and mask the gradient wrapping result based on a coherence map corresponding to the differential interference map to obtain a second phase gradient map.

[0042] It can be understood that for each interference pair, the gradient wrapping result can be determined based on the first phase gradient map. The phase wrapping phenomenon occurs when the phase gradient value exceeds a certain range. At this time, there will be random phase gradients caused by random noise and local micro gradients caused by atmospheric errors in the phase gradient. Therefore, the phase gradient needs to be corrected to avoid calculation errors caused by wrapping. In order to solve this phenomenon, the gradient wrapping result can be masked in combination with the coherence map to obtain a second phase gradient map. In this way, by retaining valid gradient information in areas with higher coherence and removing invalid or erroneous gradient information in areas with lower coherence, the calculation accuracy is improved. Among them, the second phase gradient map includes a second horizontal phase gradient map, a second vertical phase gradient map, a third diagonal phase gradient map, and a fourth diagonal phase gradient map.

[0043] Here, the present disclosure uses a winding formula to achieve the determination task of the gradient winding result. The winding formula can be expressed as: ; in, Represented as gradient winding result; It is expressed as a phase gradient.

[0044] In some other embodiments, when determining the gradient wrapping results for each interferometer pair, methods such as mass guidance and path tracking can be used to directly process the phase image to obtain a continuous phase distribution, thereby calculating a jump-free phase gradient. Furthermore, machine learning or deep learning techniques can be combined to train models to predict and correct phase wrapping issues, which are not specifically limited here.

[0045] S303, determine the first directional phase gradient atlas based on the second horizontal directional phase gradient map corresponding to each interference pair in the baseline set, determine the second directional phase gradient atlas based on the second vertical directional phase gradient map corresponding to each interference pair in the baseline set, determine the third directional phase gradient atlas based on the third diagonal directional phase gradient map corresponding to each interference pair in the baseline set, and determine the fourth directional phase gradient atlas based on the fourth diagonal directional phase gradient map corresponding to each interference pair in the baseline set.

[0046] Because noise and errors are inevitable in interferometry, the phase gradient map from a single interferometer pair may not accurately reflect the true surface deformation. Therefore, stacking the phase gradient maps from multiple interferometer pairs effectively fuses their information, suppressing noise and improving the signal-to-noise ratio (SNR) of the data, thereby more accurately reflecting the phase gradient characteristics of the surface under different time baselines.

[0047] Specifically, for different directional phase gradient atlases, phase gradient atlases in different directions are determined based on the second phase gradient maps corresponding to each interference pair. The first directional phase gradient atlas includes a first directional phase gradient stacked map corresponding to each time baseline, the second directional phase gradient atlas includes a second directional phase gradient stacked map corresponding to each time baseline, the third directional phase gradient atlas includes a third directional phase gradient stacked map corresponding to each time baseline, and the fourth directional phase gradient atlas includes a fourth directional phase gradient stacked map corresponding to each time baseline. These stacked maps can intuitively display the phase gradient changes along a certain direction on different time baselines, helping to analyze and identify the phase change characteristics under different time baselines in that direction.

[0048] Here, when the present disclosure uses a stacking formula to stack the phase gradient images corresponding to each interference pair, the formula can be expressed as: ; in, It is represented by the phase gradient stacking result in the k direction, where the values ​​of k are 0, 45, 90, and 135; m is represented by the number of interferometer pairs participating in the phase gradient stacking calculation; It is represented as the time basis of the differential interferogram corresponding to the nth interference pair.

[0049] S304 , determining the time baseline phase gradient rate set based on the first direction phase gradient stacked graph, the second direction phase gradient stacked graph, the third direction phase gradient stacked graph, and the fourth direction phase gradient stacked graph corresponding to each time baseline.

[0050] It can be understood that by stacking the phase gradient atlases in each direction corresponding to each time baseline, and calculating the rate of change of these stacked atlases in the time dimension, that is, the phase gradient rate, the phase gradient rate corresponding to each time baseline can be obtained, and then the time baseline phase gradient rate set can be obtained. The phase gradient rate reflects how fast the phase gradient changes over time.

[0051] In some possible embodiments, in order to more accurately determine the time baseline phase gradient rate set, the original basis can be further refined by introducing phase gradient rate determination steps (1) to (3) for a specific time baseline (such as the shortest time baseline, the first time baseline, the second time baseline, etc.): (1) determining the shortest time baseline phase gradient rate corresponding to the shortest time baseline based on a gradient fusion formula, a first direction phase gradient stacked graph, a second direction phase gradient stacked graph, a third direction phase gradient stacked graph, and a fourth direction phase gradient stacked graph corresponding to the shortest time baseline; (2) determining a first time baseline phase gradient rate corresponding to the first time baseline based on a gradient fusion formula, a first direction phase gradient stacked graph, a second direction phase gradient stacked graph, a third direction phase gradient stacked graph, and a fourth direction phase gradient stacked graph corresponding to the first time baseline; (3) Determine the second time baseline phase gradient rate corresponding to the second time baseline based on the gradient fusion formula, the first direction phase gradient stacked graph, the second direction phase gradient stacked graph, the third direction phase gradient stacked graph, and the fourth direction phase gradient stacked graph corresponding to the second time baseline.

[0052] Here, due to the morphology of the surface deformation field, the phase gradient stacking results in the four directions have different forms and need to be fused to obtain the complete gradient rate information of the surface deformation. When performing phase gradient fusion, in order to avoid the gradients in different directions from canceling each other out, the method of taking the absolute value average is adopted, that is, the gradient fusion formula can be expressed as: ; Among them, G represents the fusion result of the four-directional phase gradient stacking images; It is represented as a phase gradient stacking diagram in the k direction, where the values ​​of k are 0, 45, 90, and 135. Represented as the first direction, Represented as the second direction, Indicated as the third direction, Indicated as the fourth direction.

[0053] It can be understood that the shortest time baseline phase gradient rate provides an initial reference for subsequent analysis of surface deformation or topographic change. By calculating the phase gradient rate corresponding to each time baseline, the rate of surface deformation or topographic change within a specific time period can be obtained.

[0054] The minimum time baseline mentioned in this disclosure is 12 days, the first time baseline is 24 days, and the second time baseline is 36 days. These time baselines can be set as standard references. In other embodiments, if specific needs arise, the calculation can be further extended to consider longer time baselines (such as 48 days, 60 days, etc.), thereby providing a wider analysis scope. By flexibly adjusting the length of the time baseline, different analysis requirements can be met, ensuring the provision of more targeted analysis results.

[0055] S103: Construct a Delaunay triangulation based on the differential interference pattern and the coherence map, and weight each triangular arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set to obtain a weight corresponding to each triangular arc segment in the Delaunay triangulation.

[0056] It can be understood that the Delaunay triangulation is a method of connecting a set of points on a plane into a triangular mesh. These triangles satisfy the empty circumcircle property, that is, the circumcircle of any triangle does not contain other points. In geographic information systems and remote sensing data analysis, the Delaunay triangulation is often used to interpolate discrete surface deformation data points into a continuous deformation surface for easy analysis and visualization. In this solution, the Delaunay triangulation is used to divide the area to be monitored into multiple triangular regions and further used to define the constraints for phase unwrapping.

[0057] Specifically, after constructing the Delaunay triangulation, weights are assigned to each triangular arc segment in the triangulation based on the coherence map and the time-baseline phase gradient rate set. The weights reflect the reliability or importance of different arc segments in deformation analysis. For example, regions with high coherence and arc segments with stable deformation rates may be assigned higher weights because they provide more reliable deformation information. This process helps improve the accuracy and robustness of the deformation model.

[0058] For example, referring to Figure 4 As shown, when constructing a Delaunay triangulation, the following steps S401 to S402 may be included: S401 , for each differential interferogram corresponding to the interference pair, extract high-quality points that meet a preset coherence threshold based on a coherence map corresponding to the differential interferogram, and determine a high-quality point set corresponding to the interference pair based on the extraction result.

[0059] Specifically, the coherence graph extracts high-quality points that meet a preset coherence threshold. These points represent highly coherent ground features within the monitored area. Selecting high-quality points ensures that the point set used in the subsequent triangulation process is highly accurate and representative, guaranteeing the quality of the points in the Delaunay triangulation. In this disclosure, the preset coherence threshold is set to 0.4. In other embodiments, this threshold can be adjusted based on actual conditions and is not specifically limited here.

[0060] S402 : For each high-quality point set corresponding to an interference pair, generate a Delaunay triangulation corresponding to the interference pair using a point-by-point interpolation method.

[0061] It can be understood that after extracting the high-quality point set corresponding to each interference pair, a point-by-point interpolation method can be used to generate a Delaunay triangulation based on these extracted high-quality points. Each Delaunay triangulation consists of multiple triangles. Point-by-point interpolation is a classic meshing method that adjusts the existing triangle structure by adding new points one by one, ultimately forming a mesh that meets the Delaunay conditions.

[0062] In some possible embodiments, in order to improve the accuracy and efficiency of surface deformation analysis, especially when dealing with complex terrain and multiple interference factors, the weight allocation process can be further optimized in the process of weighting each triangle arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set. Figure 5 As shown, the following steps S501 to S503 may be included: S501 : for each interference pair corresponding to the Delaunay triangulation, determine whether there is a triangle arc segment corresponding to the shortest time baseline in the Delaunay triangulation.

[0063] Specifically, the Delaunay triangulation generated by each interferometer pair is examined to identify whether there are triangular arc segments defined by the shortest time baseline (i.e., the interferometer pair with the shortest time interval between the two observations). The shortest time baseline generally means a higher quality interferogram because the shorter time interval reduces the influence of factors other than surface changes (such as atmospheric interference and vegetation growth) on the interferometric measurement.

[0064] S502: If it exists, the weight of the triangle arc segment corresponding to the shortest time baseline is set to a preset coherence value.

[0065] Here, if triangular arc segments defined by the shortest time baseline are found, these arc segments will be assigned a preset, relatively high coherence value as a weight. This preset value is usually based on historical data, experimental verification or expert judgment, and is intended to reflect the high reliability and importance of these arc segments in deformation analysis. In this disclosure, the preset coherence value is set to 1.

[0066] S503: If not, weight each triangular arc segment in the Delaunay triangulation based on the arc segment weight calculation formula, the coherence graph, and the time baseline phase gradient rate set.

[0067] Here, if there is no triangle arc segment defined by the shortest time baseline in the Delaunay triangulation, or in order to more comprehensively assign weights to all arc segments, the arc segment weight calculation formula can be used to determine the weights. The arc segment weight calculation formula can be expressed as: ; in, Represented as arc segments The weight of arc The starting point is , the end point is ; Represented as arc segments The phase gradient rate at the starting position; Represented as arc segments The phase gradient rate at the endpoint position; It is expressed as the phase gradient rate threshold, which is , Expressed as the phase gradient rate; Represented as arc segments The coherence value of the starting position; Represented as arc segments The coherence value of the end position.

[0068] S104 , completing the initial phase unwrapping task of the area to be monitored based on the minimum cost flow method and the weights corresponding to the triangular arc segments in the Delaunay triangulation to obtain an initial phase unwrapping result.

[0069] The Minimum Cost Flow (MCF) method is an optimization algorithm in graph theory. In a network, each edge has a weight or cost. The algorithm's goal is to find the minimum-cost path or flow distribution scheme from a source to a sink in a given network. In the context of phase unwrapping, the MCF method is used to find a phase value distribution scheme that minimizes the cumulative cost of the phase difference between adjacent pixels (i.e., the phase gradient of an arc segment).

[0070] Specifically, after weighting each triangular arc segment in the Delaunay triangulation network, a minimum cost flow algorithm is used to optimize flow distribution within the network, decoupling phase information and effectively completing phase unwrapping. This initial phase unwrapping result includes the unwrapped phases corresponding to multiple time baselines. By guiding the unwrapping network with the phase gradient rate, the phase unwrapping results for most interferograms can be accurately recovered.

[0071] S105 , performing a pseudo-accuracy assessment on the unwrapped phase corresponding to each time baseline, and adjusting the initial phase unwrapping result based on the pseudo-accuracy assessment result to obtain a target phase unwrapping result.

[0072] It is understandable that in the process of phase unwrapping, in order to ensure the accurate recovery of deformation information, especially for areas where gradients are entangled at the boundaries of the deformation field, the complexity of the deformation field boundaries may cause phase gradients to be confused in these areas. Simply relying on the minimum cost flow method often cannot provide accurate unwrapping results.

[0073] In this regard, the present disclosure proposes to perform quasi-accurate identification on the unwrapped phase corresponding to each time baseline respectively, and adjust the initial phase unwrapping result based on the quasi-accurate identification result to obtain the target phase unwrapping result. Here, quasi-accurate identification refers to a process of verification and calibration of the unwrapped phase data, aiming to evaluate its accuracy and reliability. By performing quasi-accurate identification on the unwrapped phase corresponding to each time baseline, the phase changes of the entire deformation field at different time points can be comprehensively evaluated, thereby more accurately capturing the deformation information. Based on the quasi-accurate identification, the initial phase unwrapping result is adjusted. The adjustment process may involve correction of the phase value, re-unwrapping of the erroneous area, or smoothing of the overall phase field to ensure that the deformation information of all areas can be accurately restored to obtain the target phase unwrapping result.

[0074] For example, referring to Figure 6 As shown, when performing a pseudo-accuracy assessment on the unwrapped phase corresponding to each time baseline and adjusting the initial phase unwrapping result based on the pseudo-accuracy assessment result, the following steps S601 to S605 may be included: S601: Determine a first matrix to be approximated based on the initial phase unwrapping result.

[0075] Here, the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the first time baseline are extracted from the initial phase unwrapping results, and together they form the first matrix to be approximated. The shortest time baseline generally refers to the shortest interval between two adjacent observation moments in a time series. The unwrapped phase corresponding to this short time baseline is considered relatively stable due to its short time interval and small deformation change, i.e., a "quasi-stable phase." However, the unwrapped phase corresponding to the first time baseline may have some uncertainty due to the complexity of the deformation field or observation errors, i.e., a "non-quasi-stable phase."

[0076] S602: For the first matrix to be approximated, mark the unwrapped phase corresponding to the shortest time baseline as a first quasi-stable phase, and mark the unwrapped phase corresponding to the first time baseline as a first non-quasi-stable phase; and determine a first gross error estimate based on the gross error observation equation, the first quasi-stable phase, and the first non-quasi-stable phase.

[0077] Specifically, for the first matrix to be quasi-stable, the unwrapped phase corresponding to the shortest time baseline can be marked as the first quasi-stable phase, and the unwrapped phase corresponding to the first time baseline can be marked as the first non-quasi-stable phase. Then, the gross error observation equation is used to calculate the first gross error estimate, which reflects the degree of difference between the unwrapped phase corresponding to the first time baseline and the unwrapped phase corresponding to the shortest time baseline, that is, the possible error size. Here, the gross error observation equation can be expressed as: ; ; Where V represents the observation residual; A represents the n×m dimension coefficient matrix, n represents the number of phases to be re-unwrapped by the matrix to be estimated, and m represents the number of SAR images in the initial SAR image set; Expressed as a hyperparameter; It is represented as the gross error estimate; L is represented as the matrix to be estimated; It is represented as a matrix composed of multiple non-quasi-stable phases; P is represented as the unit matrix; R is represented as the adjustment factor.

[0078] S603: Adjust the unwrapped phase corresponding to the first time baseline in the initial phase unwrapping result based on the first gross error estimate, and determine a second matrix to be approximated based on the adjusted initial phase unwrapping result.

[0079] It is understandable that after calculating the first gross error estimate, the unwrapped phase corresponding to the first time baseline in the initial phase unwrapping result can be adjusted based on the first gross error estimate to reduce the error so that the adjusted unwrapped phase will be closer to the true value. Here, the adjustment of the initial phase unwrapping result can be specifically manifested as: constraining the obtained first gross error estimate (for example, rounding it off) to The constrained gross error is then removed from the corresponding phase.

[0080] Specifically, after completing the phase adjustment task of the first time baseline, a second matrix to be estimated is constructed based on the adjusted initial phase unwrapping result. The matrix includes the unwrapped phase corresponding to the shortest time baseline, the unwrapped phase corresponding to the adjusted first time baseline, and the unwrapped phase corresponding to the second time baseline.

[0081] S604, for the second matrix to be approximated, mark the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the adjusted first time baseline as a second quasi-stable phase, and mark the unwrapped phase corresponding to the second time baseline as a second non-quasi-stable phase; and determine the second gross error estimate based on the gross error observation equation, the second quasi-stable phase and the second non-quasi-stable phase.

[0082] Specifically, for the second matrix to be approximated, the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the adjusted first time baseline can be marked as the second quasi-stable phase, and the unwrapped phase corresponding to the second time baseline can be marked as the second non-quasi-stable phase. Then, similar to step S602, the gross error observation equation is again used to calculate a second gross error estimate. This estimate reflects the degree of difference between the unwrapped phase corresponding to the second time baseline and the quasi-stable phase.

[0083] S605: Adjust the unwrapping phase corresponding to the second time baseline in the adjusted initial phase unwrapping result based on the second gross error estimate to obtain the target phase unwrapping result.

[0084] It is understandable that, similar to step S603, the unwrapped phase corresponding to the second time baseline in the adjusted initial phase unwrapping result is adjusted based on the second gross error estimate. In this way, the error can be further reduced, and a more accurate target phase unwrapping result can be obtained.

[0085] In some other embodiments, if there are other longer time baselines (such as 48 days, 60 days, etc.), this process can be iterated until the unwrapped phases corresponding to all time baselines are accurately identified and adjusted, thereby ensuring that the deformation information of the entire deformation field can be accurately recovered.

[0086] The phase unwrapping method and apparatus based on improved minimum cost flow and quasi-identification provided in the embodiments of the present disclosure effectively reduce error accumulation during the phase unwrapping process by optimizing the calculation of phase gradient rate sets and the weighting of triangular arc segments using the improved minimum cost flow method, thereby ensuring high-precision unwrapping results and accelerating the unwrapping process. Furthermore, the present disclosure utilizes quasi-identification technology to post-process the initial phase unwrapping results, further improving phase unwrapping accuracy and enabling high-precision surface deformation monitoring.

[0087] In order to test the performance of the phase unwrapping method proposed in this disclosure, this disclosure simulates three common types of surface deformation rate fields (such as Figure 7 (a) and (b) are shown to analyze the impact of different deformation characteristics on phase unwrapping. The maximum deformation rate of these three deformation fields is 25 cm / a, but the spatial distribution characteristics of the deformation rates are different. Deformation field 1 gradually increases from 0 to 25 cm / a from top to bottom, with clear deformation boundaries between the sides and bottom and the stable region, but no obvious deformation boundary on the top. Deformation field 2 has a significant deformation boundary between the deformation and the stable region, and the deformation rate within this region is 25 cm / a. Deformation field 3 gradually increases from 0 to 25 cm / a from the boundary to the center, with no obvious deformation boundary compared to the surrounding stable region. To facilitate the distinction between these types of deformation rate fields, this paper defines their shapes as trapezoidal, rectangular, and circular, respectively. Subsequently, this paper acquired orbital information for 30 Sentinel-1 descending images from January 2, 2019, to December 28, 2019. 153 interferograms were generated based on the time baseline threshold of 72 days and the space baseline threshold of 100 meters. 153 original interferograms were generated based on the time and space baselines, and the phases of the original interferograms were wrapped to When generating the simulated interference pattern, a periodic deformation is added on the basis of the linear deformation. The first six interference phase diagrams generated by simulation are as follows: Figure 8 As shown, (a)-(f) are the simulated original interference phase diagrams, and their time baselines are 12 days, 24 days, 36 days, 48 ​​days, 60 days and 72 days respectively. (g)-(l) are the interference diagrams obtained by wrapping (a)-(f). Figure 8 It can be seen that the phase of the entanglement interference pattern is When the time baseline is greater than 24 days, phase jumps appear inside the circular deformation field, inside and at the boundary of the trapezoidal deformation field, and at the boundary of the rectangular deformation field.

[0088] After generating the simulated entangled phase diagram, the present disclosure performs phase unwrapping by the traditional MCF method and the MCF method based on the phase gradient rate. The unwrapping results of the traditional MCF method are as follows: Figure 9 As shown, Figure 9 (a) is the unwrapped phase diagram obtained by the MCF method from the wrapped phase diagram with a time baseline of 12 days. Similarly, Figure 9 .(b)-9.(f) Unwrapped phase diagrams obtained by the MCF method from the wrapped phase diagrams of 24 days, 36 days, 48 ​​days, 60 days, and 72 days, respectively. Figure 9 (g)-(l) are the differences between the unwrapped phase image and the original phase image of the same time baseline. Figure 9 As can be seen in (a), when the time baseline is 12 days, the phase is correctly recovered and no unwrapping error occurs. Figure 9 (g) only has errors due to computer storage accuracy. When the time baseline is greater than 12 days, the circular deformation area obtains correct unwrapping results, while the trapezoidal and rectangular deformation areas obtain incorrect unwrapping results.

[0089] Specifically, the unwrapping results of the MCF method guided by the phase gradient rate are as follows: Figure 10 As shown, Figure 10 (a) is the unwrapped phase image obtained by the proposed phase unwrapping method from the wrapped phase image with a time baseline of 12 days. Similarly, Figure 10 (b)-10. (f) Unwrapped phase images obtained by phase unwrapping from the wrapped phase images of 24 days, 36 days, 48 ​​days, 60 days, and 72 days, respectively; Figure 10 (g)-(l) are the differences between the unwrapped phase map and the simulated original phase map of the same time baseline. When there is an unwrapped area at the boundary of the deformation field, the correct unwrapping result is obtained. However, with the time accumulation effect of the surface deformation, the deformation gradient at the deformation boundary gradually increases, and the MCF method based on the phase gradient rate obtains an incorrect phase unwrapping result. Then, after the unwrapping phase error correction (i.e., the quasi-identification adjustment), the result is as follows Figure 11 As shown, Figure 11 (a)-11. (f) The results of unwrapping phase correction of the wrapped phase diagrams of 12 days, 24 days, 36 days, 48 ​​days, 60 days, and 72 days respectively; Figure 11 (g) - (l) are the differences between the unwrapped phase image and the simulated original phase image for the same time baseline. As can be seen from the figures, the unwrapping error has been almost completely corrected.

[0090] Those skilled in the art will understand that in the above-mentioned method of the specific implementation method, the writing order of each step does not mean a strict execution order and does not constitute any limitation on the implementation process. The specific execution order of each step should be determined by its function and possible internal logic.

[0091] Based on the same inventive concept, the embodiment of the present disclosure also provides a phase unwrapping device based on improved minimum cost flow and quasi-identification corresponding to the phase unwrapping method based on improved minimum cost flow and quasi-identification. Since the principle of solving the problem by the device in the embodiment of the present disclosure is similar to the above-mentioned phase unwrapping method based on improved minimum cost flow and quasi-identification in the embodiment of the present disclosure, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be repeated.

[0092] Reference Figure 12 FIG. 1 is a schematic diagram of a phase unwrapping device 1200 based on improved minimum cost flow and quasi-identification according to an embodiment of the present disclosure, wherein the device includes: The data acquisition module 1201 is used to acquire a SAR data set and a DEM data of the area to be monitored; and generate a differential interferogram and a coherence map based on the SAR data set and the DEM data; A rate calculation module 1202 is configured to calculate a time baseline phase gradient rate set based on the differential interferogram and the coherence map; A weight determination module 1203 is configured to construct a Delaunay triangulation based on the differential interference pattern and the coherence map, and to weight each triangle arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set to obtain a weight corresponding to each triangle arc segment in the Delaunay triangulation; A phase unwrapping module 1204 is configured to complete an initial phase unwrapping task for the monitored area based on a minimum cost flow method and the weights corresponding to the triangular arc segments in the Delaunay triangulation, thereby obtaining an initial phase unwrapping result; wherein the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines; The phase alignment module 1205 is configured to perform alignment identification on the unwrapped phase corresponding to each time baseline, and adjust the initial phase unwrapping result based on the alignment identification result to obtain a target phase unwrapping result.

[0093] In some possible embodiments, the data acquisition module 1201 is specifically configured to: Determining an initial SAR image set based on the SAR data set; wherein the initial SAR image set includes a plurality of SAR images of the area to be monitored; Determining a SAR main image in the initial SAR image set, registering other SAR images in the initial SAR image set except the SAR main image based on the DEM data using the SAR main image as a reference, and determining a target SAR image set based on the registration results; wherein the target SAR image set includes multiple target SAR images; Determining a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; wherein the baseline set includes a plurality of interference pairs, each of the interference pairs including two target SAR images; A differential interference map corresponding to each interference pair is generated based on the DEM data; and a coherence map corresponding to each interference pair is generated based on the differential interference map corresponding to each interference pair.

[0094] In some possible embodiments, the rate calculation module 1202 is specifically configured to: For the differential interferogram corresponding to each interference pair, calculate the phase gradient value corresponding to each pixel in the differential interferogram, and determine a first phase gradient map based on the phase gradient value corresponding to each pixel; wherein the first phase gradient map includes a first horizontal phase gradient map, a first vertical phase gradient map, a first diagonal phase gradient map, and a second diagonal phase gradient map; For each interference pair, a gradient wrapping result is determined based on a first phase gradient map corresponding to the differential interference map, and the gradient wrapping result is masked based on a coherence map corresponding to the differential interference map to obtain a second phase gradient map; wherein the second phase gradient map includes a second horizontal phase gradient map, a second vertical phase gradient map, a third diagonal phase gradient map, and a fourth diagonal phase gradient map; A first directional phase gradient atlas is determined based on the second horizontal directional phase gradient map corresponding to each interference pair in the baseline set, a second directional phase gradient atlas is determined based on the second vertical directional phase gradient map corresponding to each interference pair in the baseline set, a third directional phase gradient atlas is determined based on the third diagonal directional phase gradient map corresponding to each interference pair in the baseline set, and a fourth directional phase gradient atlas is determined based on the fourth diagonal directional phase gradient map corresponding to each interference pair in the baseline set; wherein, the first directional phase gradient atlas includes a first directional phase gradient stacking map corresponding to each time baseline, the second directional phase gradient atlas includes a second directional phase gradient stacking map corresponding to each time baseline, the third directional phase gradient atlas includes a third directional phase gradient stacking map corresponding to each time baseline, and the fourth directional phase gradient atlas includes a fourth directional phase gradient stacking map corresponding to each time baseline; The time baseline phase gradient rate set is determined based on the first direction phase gradient stacking map, the second direction phase gradient stacking map, the third direction phase gradient stacking map, and the fourth direction phase gradient stacking map corresponding to each time baseline.

[0095] In some possible embodiments, the time baseline includes a shortest time baseline, a first time baseline, and a second time baseline; the time baseline phase gradient rate set includes a shortest time baseline phase gradient rate, a first time baseline phase gradient rate, and a second time baseline phase gradient rate; and the rate calculation module 1202 is specifically configured to: Determine the shortest time baseline phase gradient rate corresponding to the shortest time baseline based on a gradient fusion formula, a first direction phase gradient stacked map, a second direction phase gradient stacked map, a third direction phase gradient stacked map, and a fourth direction phase gradient stacked map corresponding to the shortest time baseline; Determine a first time baseline phase gradient rate corresponding to the first time baseline based on a gradient fusion formula, a first direction phase gradient stacked graph, a second direction phase gradient stacked graph, a third direction phase gradient stacked graph, and a fourth direction phase gradient stacked graph corresponding to the first time baseline; Determine a second time baseline phase gradient rate corresponding to the second time baseline based on a gradient fusion formula, a first direction phase gradient stacked map, a second direction phase gradient stacked map, a third direction phase gradient stacked map, and a fourth direction phase gradient stacked map corresponding to the second time baseline; The gradient fusion formula includes: ; Among them, G represents the fusion result of the four-directional phase gradient stacking images; It is represented as a phase gradient stacking diagram in the k direction, where the values ​​of k are 0, 45, 90, and 135. Represented as the first direction, Represented as the second direction, Indicated as the third direction, Indicated as the fourth direction.

[0096] In some possible embodiments, the weight determination module 1203 is specifically configured to: For each differential interferogram corresponding to the interference pair, extracting high-quality points that meet a preset coherence threshold based on a coherence map corresponding to the differential interferogram, and determining a high-quality point set corresponding to the interference pair based on the extraction result; For each high-quality point set corresponding to the interference pair, a point-by-point interpolation method is used to generate a Delaunay triangulation corresponding to the interference pair; wherein each Delaunay triangulation includes a plurality of triangles; The weight determination module 1203 is specifically used to: For each interference pair corresponding to the Delaunay triangulation, determining whether there is a triangle arc segment corresponding to the shortest time baseline in the Delaunay triangulation; If it exists, the weight of the triangle arc segment corresponding to the shortest time baseline is set to a preset coherence value; If not, weighting each triangular arc segment in the Delaunay triangulation based on an arc segment weight calculation formula, the coherence graph, and the time baseline phase gradient rate set; The arc weight calculation formula includes: ; in, Represented as arc segments The weight of arc The starting point is , the end point is ; Represented as arc segments The phase gradient rate at the starting position; Represented as arc segments The phase gradient rate at the endpoint position; It is expressed as the phase gradient rate threshold, and its value is , Expressed as the phase gradient rate; Represented as arc segments The coherence value of the starting position; Represented as arc segments The coherence value of the end position.

[0097] In some possible embodiments, the phase alignment module 1205 is specifically configured to: Determine a first matrix to be estimated based on the initial phase unwrapping result; wherein the first matrix to be estimated includes an unwrapped phase corresponding to the shortest time baseline and an unwrapped phase corresponding to the first time baseline; For the first matrix to be approximated, marking the unwrapped phase corresponding to the shortest time baseline as a first quasi-stable phase, and marking the unwrapped phase corresponding to the first time baseline as a first non-quasi-stable phase; and determining a first gross error estimate based on a gross error observation equation, the first quasi-stable phase, and the first non-quasi-stable phase; Adjusting the unwrapped phase corresponding to the first time baseline in the initial phase unwrapping result based on the first gross error estimate, and determining a second matrix to be accurately estimated based on the adjusted initial phase unwrapping result; wherein the second matrix to be accurately estimated includes the unwrapped phase corresponding to the shortest time baseline, the unwrapped phase corresponding to the adjusted first time baseline, and the unwrapped phase corresponding to the second time baseline; For the second matrix to be accurately approximated, marking the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the adjusted first time baseline as a second quasi-stable phase, and marking the unwrapped phase corresponding to the second time baseline as a second non-quasi-stable phase; and determining a second gross error estimate based on the gross error observation equation, the second quasi-stable phase, and the second non-quasi-stable phase; The unwrapping phase corresponding to the second time baseline in the adjusted initial phase unwrapping result is adjusted based on the second gross error estimate to obtain the target phase unwrapping result.

[0098] In some possible embodiments, the gross error observation equation includes: ; ; Where V represents the observation residual; A represents the n×m dimension coefficient matrix, n represents the number of phases to be re-unwrapped by the matrix to be estimated, and m represents the number of SAR images in the initial SAR image set; Expressed as a hyperparameter; It is represented as the gross error estimate; L is represented as the matrix to be estimated; It is represented as a matrix composed of multiple non-quasi-stable phases; P is represented as the unit matrix; R is represented as the adjustment factor.

[0099] Based on the same technical concept, the embodiment of the present disclosure also provides a computer device. Figure 13 1 is a schematic diagram of the structure of a computer device 1300 provided in an embodiment of the present disclosure, including a processor 1301, a memory 1302, and a bus 1303. The memory 1302 is used to store execution instructions and includes a memory 13021 and an external memory 13022. The memory 13021 is also referred to as internal memory and is used to temporarily store calculation data in the processor 1301 and data exchanged with an external memory 13022 such as a hard disk. The processor 1301 exchanges data with the external memory 13022 through the memory 13021.

[0100] In the embodiments of the present application, memory 1302 is specifically used to store application code for executing the solutions of the present application, and its execution is controlled by processor 1301. Specifically, when computer device 1300 is running, processor 1301 communicates with memory 1302 via bus 1303, enabling processor 1301 to execute the application code stored in memory 1302, thereby performing the methods described in any of the aforementioned embodiments. Memory 1302 may be, but is not limited to, random access memory, read-only memory, programmable read-only memory, erasable read-only memory, electrically erasable read-only memory, etc. Processor 1301 may be an integrated circuit chip with signal processing capabilities. Such processors may be general-purpose processors, including central processing units (CPUs) and network processors; they may also be digital signal processors, application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. These processors may implement or execute the various methods, steps, and logic block diagrams disclosed in the embodiments of the present invention. A general-purpose processor may be a microprocessor or any conventional processor.

[0101] It should be understood that the structure illustrated in the embodiments of the present application does not constitute a specific limitation on the computer device 1300. In other embodiments of the present application, the computer device 1300 may include more or fewer components than shown, or may combine or separate certain components, or arrange the components differently. The illustrated components may be implemented in hardware, software, or a combination of software and hardware.

[0102] The present disclosure also provides a computer-readable storage medium having a computer program stored thereon. When executed by a processor, the computer program executes the steps of the phase unwrapping method based on improved minimum cost flow and quasi-identification described in the above method embodiment. The storage medium may be a volatile or non-volatile computer-readable storage medium.

[0103] An embodiment of the present disclosure also provides a computer program product, which carries a program code. The instructions included in the program code can be used to execute the steps of the phase unwrapping method based on improved minimum cost flow and quasi-identification described in the above method embodiment. For details, please refer to the above method embodiment, which will not be repeated here.

[0104] The computer program product may be implemented in hardware, software, or a combination thereof. In one embodiment, the computer program product is implemented as a computer storage medium. In another embodiment, the computer program product is implemented as a software product, such as a software development kit.

[0105] Those skilled in the art will clearly understand that, for ease of description and brevity, the specific operating processes of the systems and devices described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here. In the several embodiments provided in this disclosure, it should be understood that the disclosed systems and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of the units described is merely a logical functional division, and actual implementations may employ other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be omitted or not implemented. Furthermore, the coupling or direct coupling or communication connection shown or discussed between each other may be through some communication interface, or the indirect coupling or communication connection between devices or units may be electrical, mechanical, or other forms. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, i.e., they may be located in one location or distributed across multiple network units. Some or all of the units may be selected to achieve the objectives of the present embodiment according to actual needs. In addition, the functional units in the various embodiments of the present disclosure can be integrated into a processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. If the function is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a non-volatile computer-readable storage medium that can be executed by a processor. Based on this understanding, the technical solution of the present disclosure, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the method described in the various embodiments of the present disclosure. The aforementioned storage medium includes: various media that can store program code, such as a USB flash drive, a mobile hard disk, a read-only memory, a random access memory, a magnetic disk, or an optical disk.

[0106] Finally, it should be noted that the above-described embodiments are only specific implementation methods of the present disclosure, which are used to illustrate the technical solutions of the present disclosure, rather than to limit them. The scope of protection of the present disclosure is not limited thereto. Although the present disclosure has been described in detail with reference to the above-described embodiments, those skilled in the art should understand that any person skilled in the art can modify or easily conceive of changes to the technical solutions described in the above-described embodiments within the technical scope disclosed in the present disclosure, or replace some of the technical features therein with equivalents. Such modifications, changes, or replacements do not deviate from the spirit and scope of the technical solutions of the embodiments of the present disclosure, and should be included in the scope of protection of the present disclosure. Therefore, the scope of protection of the present disclosure shall be subject to the scope of protection of the claims.

Claims

1. A phase unwrapping method based on improved minimum cost flow and quasi-identification, characterized in that: include: Obtain SAR data set and DEM data of the area to be monitored; and generating a differential interferogram and a coherence map based on the SAR data set and the DEM data; calculating a time baseline phase gradient rate set based on the differential interferogram and the coherence map; Constructing a Delaunay triangulation network based on the differential interference graph and the coherence graph, and weighting each triangular arc segment in the Delaunay triangulation network based on the coherence graph and the time baseline phase gradient rate set to obtain a weight corresponding to each triangular arc segment in the Delaunay triangulation network; Based on the minimum cost flow method and the weights corresponding to the triangular arc segments in the Delaunay triangulation network, the initial phase unwrapping task of the monitored area is completed to obtain an initial phase unwrapping result; wherein the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines; The unwrapping phase corresponding to each time baseline is respectively subjected to a quasi-accurate identification, and the initial phase unwrapping result is adjusted based on the quasi-accurate identification result to obtain a target phase unwrapping result.

2. The phase unwrapping method based on improved minimum cost flow and quasi-identification according to claim 1, characterized in that: The generating of a differential interference graph and a coherence graph based on the SAR data set and the DEM data comprises: Determine an initial SAR image set based on the SAR data set; wherein the initial SAR image set includes a plurality of SAR images of the area to be monitored; Determine a SAR main image in the initial SAR image set, take the SAR main image as a reference, register other SAR images in the initial SAR image set except the SAR main image according to the DEM data, and determine a target SAR image set based on the registration result; wherein the target SAR image set includes a plurality of target SAR images; Determine a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; wherein the baseline set includes a plurality of interference pairs, each of the interference pairs includes two target SAR images; A differential interference map corresponding to each interference pair is generated based on the DEM data; and a coherence map corresponding to each interference pair is generated based on the differential interference map corresponding to each interference pair.

3. The phase unwrapping method based on improved minimum cost flow and quasi-identification according to claim 2, characterized in that: The step of calculating a time baseline phase gradient rate set based on the differential interference graph and the coherence graph comprises: For the differential interference graph corresponding to each interference pair, the phase gradient value corresponding to each pixel in the differential interference graph is calculated, and a first phase gradient graph is determined based on the phase gradient value corresponding to each pixel; wherein the first phase gradient graph includes a first horizontal direction phase gradient graph, a first vertical direction phase gradient graph, a first diagonal direction phase gradient graph, and a second diagonal direction phase gradient graph; For each interference pair, a gradient winding result is determined according to a first phase gradient map corresponding to the differential interference map, and the gradient winding result is masked based on a coherence map corresponding to the differential interference map to obtain a second phase gradient map; wherein the second phase gradient map includes a second horizontal direction phase gradient map, a second vertical direction phase gradient map, a third diagonal direction phase gradient map, and a fourth diagonal direction phase gradient map; A first directional phase gradient atlas is determined based on a second horizontal directional phase gradient map corresponding to each interference pair in the baseline set, a second directional phase gradient atlas is determined based on a second vertical directional phase gradient map corresponding to each interference pair in the baseline set, a third directional phase gradient atlas is determined based on a third diagonal directional phase gradient map corresponding to each interference pair in the baseline set, and a fourth directional phase gradient atlas is determined based on a fourth diagonal directional phase gradient map corresponding to each interference pair in the baseline set; wherein the first directional phase gradient atlas includes a first directional phase gradient stacking map corresponding to each time baseline, the second directional phase gradient atlas includes a second directional phase gradient stacking map corresponding to each time baseline, the third directional phase gradient atlas includes a third directional phase gradient stacking map corresponding to each time baseline, and the fourth directional phase gradient atlas includes a fourth directional phase gradient stacking map corresponding to each time baseline; The time baseline phase gradient rate set is determined based on the first direction phase gradient stacking diagram, the second direction phase gradient stacking diagram, the third direction phase gradient stacking diagram, and the fourth direction phase gradient stacking diagram corresponding to each time baseline.

4. The phase unwrapping method based on improved minimum cost flow and quasi-identification according to claim 3, characterized in that: The time baseline includes a shortest time baseline, a first time baseline and a second time baseline, and the time baseline phase gradient rate set includes a shortest time baseline phase gradient rate, a first time baseline phase gradient rate and a second time baseline phase gradient rate; the determining of the time baseline phase gradient rate set based on a first direction phase gradient stacking diagram, a second direction phase gradient stacking diagram, a third direction phase gradient stacking diagram and a fourth direction phase gradient stacking diagram corresponding to each time baseline includes: Determine the shortest time baseline phase gradient rate corresponding to the shortest time baseline based on the gradient fusion formula, the first direction phase gradient stacking diagram, the second direction phase gradient stacking diagram, the third direction phase gradient stacking diagram, and the fourth direction phase gradient stacking diagram corresponding to the shortest time baseline; Determine a first time baseline phase gradient rate corresponding to the first time baseline based on a gradient fusion formula, a first direction phase gradient stacked graph, a second direction phase gradient stacked graph, a third direction phase gradient stacked graph, and a fourth direction phase gradient stacked graph corresponding to the first time baseline; Determine a second time baseline phase gradient rate corresponding to the second time baseline based on a gradient fusion formula, a first direction phase gradient stacked graph, a second direction phase gradient stacked graph, a third direction phase gradient stacked graph, and a fourth direction phase gradient stacked graph corresponding to the second time baseline; The gradient fusion formula includes: ； Among them, G represents the fusion result of the phase gradient stacking images in four directions; It is represented as a phase gradient stacking diagram in the k direction, where the values ​​of k are 0, 45, 90, and 135. Represents the first direction, Represented as the second direction, Indicated as the third direction, Indicated as the fourth direction.

5. The phase unwrapping method based on improved minimum cost flow and quasi-identification according to claim 4, characterized in that: The constructing of a Delaunay triangulation network based on the differential interference graph and the coherence graph comprises: For the differential interference graph corresponding to each interference pair, high-quality points satisfying a preset coherence threshold are extracted based on a coherence graph corresponding to the differential interference graph, and a high-quality point set corresponding to the interference pair is determined based on the extraction result; For each high-quality point set corresponding to each interference pair, a point-by-point interpolation method is used to generate a Delaunay triangulation corresponding to the interference pair; wherein each of the Delaunay triangulations includes a plurality of triangles; The step of weighting each triangular arc segment in the Delaunay triangulation network based on the coherence graph and the time baseline phase gradient rate set includes: For each interference pair corresponding to the Delaunay triangulation network, determine whether there is a triangle arc segment corresponding to the shortest time baseline in the Delaunay triangulation network; If it exists, the weight of the triangle arc segment corresponding to the shortest time baseline is set to a preset coherence value; If it does not exist, weighting each triangular arc segment in the Delaunay triangulation network based on the arc segment weight calculation formula, the coherence graph and the time baseline phase gradient rate set; The arc weight calculation formula includes: ； in, Represented as an arc segment The weight of The starting point is The end point is ; Represented as an arc segment The phase gradient rate at the starting position; Represented as an arc segment The phase gradient rate at the end position; It is represented as the phase gradient rate threshold, and its value is , Expressed as the phase gradient rate; Represented as an arc segment The coherence value of the starting position of ; Represented as an arc segment The coherence value of the end position.

6. The phase unwrapping method based on improved minimum cost flow and quasi-identification according to claim 4, characterized in that: The step of respectively performing a pseudo-identification on the unwrapped phase corresponding to each time baseline, and adjusting the initial phase unwrapping result based on the pseudo-identification result, comprises: Determine a first matrix to be estimated based on the initial phase unwrapping result; wherein the first matrix to be estimated includes an unwrapped phase corresponding to the shortest time baseline and an unwrapped phase corresponding to the first time baseline; For the first matrix to be quasi-stabilized, marking the unwrapped phase corresponding to the shortest time baseline as a first quasi-stable phase, and marking the unwrapped phase corresponding to the first time baseline as a first non-quasi-stable phase; and determining a first gross error estimate based on a gross error observation equation, the first quasi-stable phase, and the first non-quasi-stable phase; The unwrapped phase corresponding to the first time baseline in the initial phase unwrapping result is adjusted based on the first gross error estimate, and a second matrix to be estimated is determined based on the adjusted initial phase unwrapping result; wherein the second matrix to be estimated includes the unwrapped phase corresponding to the shortest time baseline, the unwrapped phase corresponding to the adjusted first time baseline, and the unwrapped phase corresponding to the second time baseline; For the second matrix to be quasi-accurate, the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the adjusted first time baseline are marked as second quasi-stable phases, and the unwrapped phase corresponding to the second time baseline is marked as a second non-quasi-stable phase; and a second gross error estimate is determined based on the gross error observation equation, the second quasi-stable phase and the second non-quasi-stable phase; The unwrapping phase corresponding to the second time baseline in the adjusted initial phase unwrapping result is adjusted based on the second gross error estimate to obtain the target phase unwrapping result.

7. The phase unwrapping method based on improved minimum cost flow and quasi-identification according to claim 6, characterized in that: The gross error observation equation includes: ； ； Where V represents the observation residual; A represents the n×m dimension coefficient matrix, n represents the number of phases to be re-unwrapped by the matrix to be estimated, and m represents the number of SAR images in the initial SAR image set; Expressed as a hyperparameter; It is represented as the gross error estimate; L is represented as the matrix to be estimated; It is represented as a matrix composed of multiple non-quasi-stable phases; P is represented as the unit matrix; R is represented as the adjustment factor.

8. A phase unwrapping device based on improved minimum cost flow and quasi-identification, characterized in that: include: A data acquisition module is used to obtain SAR data sets and DEM data of the area to be monitored; and generating a differential interferogram and a coherence map based on the SAR data set and the DEM data; A rate calculation module, configured to calculate a time baseline phase gradient rate set based on the differential interference graph and the coherence graph; A weight determination module, used for constructing a Delaunay triangulation network based on the differential interference graph and the coherence graph, and weighting each triangular arc segment in the Delaunay triangulation network based on the coherence graph and the time baseline phase gradient rate set to obtain a weight corresponding to each triangular arc segment in the Delaunay triangulation network; A phase unwrapping module, used to complete the initial phase unwrapping task of the monitored area based on the minimum cost flow method and the weights corresponding to each triangular arc segment in the Delaunay triangulation network, and obtain an initial phase unwrapping result; wherein the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines; The phase alignment module is used to respectively perform alignment identification on the unwrapping phase corresponding to each time baseline, and adjust the initial phase unwrapping result based on the alignment identification result to obtain the target phase unwrapping result.

9. A storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the phase unwrapping method based on improved minimum cost flow and quasi-identification according to any one of claims 1 to 7 is implemented.

10. A computer device comprising a storage medium, a processor, and a computer program stored in the storage medium and executable on the processor, characterized in that: When the processor executes the computer program, the phase unwrapping method based on improved minimum cost flow and quasi-identification according to any one of claims 1 to 7 is implemented.

Images