Phase Unwrapping Method and Device Based on Improved Minimum Cost Flow and Quasi-Criterion Identification

Abstract

The present disclosure 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-accuracy identification, including: obtaining SAR and DEM data of an area to be monitored, and generating a differential interferogram and a coherence map; calculating a set of time baseline phase gradient rates based on the differential interferogram and the coherence map; constructing a Delaunay triangulation network, and weighting each triangular arc segment according to the coherence map and the set of time baseline phase gradient rates; secondly, using the minimum cost flow method and these weights to complete the initial phase unwrapping of the area to be monitored; finally, performing quasi-accuracy identification on the unwrapped phases of each time baseline, and adjusting the initial phase unwrapping result according to the quasi-accuracy identification result to obtain the target phase unwrapping result. In this embodiment, by combining the improved minimum cost flow algorithm with the quasi-accuracy identification technology, the accuracy and reliability of phase unwrapping are improved, thereby realizing high-precision surface deformation monitoring.

Description

Technical Field

[0001] The present disclosure relates to the technical field of surface deformation monitoring. Specifically, it relates 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, and is widely used in the monitoring and early warning of natural disasters such as landslides, mining area subsidence, volcanoes, and earthquakes. Compared with traditional remote sensing technology, SAR technology has the advantages of all-weather and all-time, can penetrate clouds and monitor at night, and is particularly suitable for remote sensing monitoring under complex terrains and extreme climates.

[0003] With the development of technology, the application of SAR technology has shifted from traditional disaster response to more precise deformation monitoring. Especially, the phase-based Synthetic Aperture Radar Interferometry (InSAR) technology has become one of the core technologies in research due to its advantages in micro-surface deformation monitoring. InSAR can monitor surface deformation with centimeter-level or even millimeter-level accuracy by calculating the interference phase of two SAR images with different viewing angles, and is suitable for long-term monitoring of geological disasters.

[0004] In related technologies, the time-series monitoring method based on InSAR technology has become the mainstream technology. By selecting SAR images with tight spatial baselines and different time baselines, it reduces the decoherence problem caused by large baseline differences, improves the spatial resolution and time sampling rate of surface deformation monitoring, and provides accurate data for disaster early warning and risk assessment. However, although InSAR technology can provide accurate deformation information, during the phase unwrapping process, especially in large-gradient deformation areas, the accuracy is often challenged. At the same time, as a commonly used unwrapping technology, the minimum cost flow unwrapping method, although improved, is still limited by SAR resolution and is difficult to effectively address the phase recovery problem in large-gradient deformation areas. 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 and quasi-identification technology, the accuracy and reliability of phase unwrapping are improved, and thus high-precision surface deformation monitoring is achieved.

[0006] The embodiments of the present disclosure provide a phase unwrapping method based on improved minimum cost flow and quasi-identification, including:

[0007] Obtain the SAR data set and 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;

[0008] Calculate a set of time - baseline phase gradient rates based on the differential interferogram and the coherence map;

[0009] Construct a Delaunay triangulation based on the differential interferogram and the coherence map, and assign weights to each triangular arc segment in the Delaunay triangulation based on the coherence map and the set of time - baseline phase gradient rates to obtain weights corresponding to each triangular arc segment in the Delaunay triangulation;

[0010] Complete the initial phase unwrapping task of the area to be monitored based on the minimum - cost flow method and the weights corresponding to each triangular arc segment in the Delaunay triangulation to obtain an initial phase unwrapping result; wherein, the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines;

[0011] Perform quasi - accuracy identification on the unwrapped phases corresponding to each time baseline respectively, and adjust the initial phase unwrapping result based on the quasi - accuracy identification result to obtain a target phase unwrapping result.

[0012] An embodiment of the present disclosure provides a phase unwrapping device based on improved minimum - cost flow and quasi - accuracy identification, including:

[0013] A data acquisition module, configured to acquire an SAR data set and 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;

[0014] A rate calculation module, configured to calculate a set of time - baseline phase gradient rates based on the differential interferogram and the coherence map;

[0015] A weight determination module, configured to construct a Delaunay triangulation based on the differential interferogram and the coherence map, and assign weights to each triangular arc segment in the Delaunay triangulation based on the coherence map and the set of time - baseline phase gradient rates to obtain weights corresponding to each triangular arc segment in the Delaunay triangulation;

[0016] A phase unwrapping module, configured to complete the initial phase unwrapping task of the area to be monitored based on the minimum - cost flow method and the weights corresponding to each triangular arc segment in the Delaunay triangulation to obtain an initial phase unwrapping result; wherein, the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines;

[0017] A phase quasi - accuracy module, configured to perform quasi - accuracy identification on the unwrapped phases corresponding to each time baseline respectively, and adjust the initial phase unwrapping result based on the quasi - accuracy identification result to obtain a target phase unwrapping result.

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

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

[0020] In the phase unwrapping method and device based on improved minimum cost flow and quasi-identification provided in the embodiments of the present disclosure, first, a synthetic aperture radar (SAR) data set and a digital elevation model (DEM) data of the area to be monitored are obtained; subsequently, based on the SAR data set and the DEM data, a differential interferogram and a coherence map are generated; next, a set of time baseline phase gradient rates is calculated using the differential interferogram and the coherence map, and on this basis, a Delaunay triangulation network is constructed. Then, according to the coherence map and the set of time baseline phase gradient rates, weights are assigned to each triangular arc segment in the Delaunay triangulation network to obtain the weights corresponding to each triangular arc segment; then, using the minimum cost flow method and these weights, the initial phase unwrapping task of the area to be monitored is completed to obtain an initial phase unwrapping result. To further improve the accuracy of phase unwrapping, the present disclosure further proposes to perform quasi-identification on the unwrapped phases corresponding to each time baseline, and finally, based on the results of the quasi-identification, the initial phase unwrapping result is adjusted to finally obtain the target phase unwrapping result.

[0021] In the embodiments of the present disclosure, in the process of phase unwrapping using the improved minimum cost flow method, by optimizing the calculation of the phase gradient rate set and the weight assignment of the triangular arc segments, the error accumulation is effectively reduced, thereby ensuring a high-precision unwrapping result and accelerating the unwrapping process. At the same time, the present disclosure also uses the quasi-identification technology to post-process the initial phase unwrapping result, further improving the phase unwrapping accuracy, and then realizing high-precision surface deformation monitoring.

[0022] To make the above objects, features, and advantages of the present disclosure more obvious and understandable, the following specifically enumerates preferred embodiments and, in conjunction with the accompanying drawings, makes a detailed description as follows. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] To more clearly illustrate the technical solutions of the embodiments of the present disclosure, the accompanying drawings required to be cited in the embodiments will be briefly introduced below. The accompanying drawings here are incorporated into the specification and constitute a part of this specification. These accompanying drawings show embodiments that conform to 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 accompanying drawings only show some embodiments of the present disclosure and should not be regarded as limiting the scope. For those of ordinary skill in the art, without creative efforts, other related accompanying drawings can also be obtained based on these accompanying drawings.

[0024] Figure 1 Shows a flowchart of a phase unwrapping method based on improved minimum cost flow and quasi-criterion identification provided by the embodiments of the present disclosure;

[0025] Figure 2 Shows a flowchart of a method for determining a differential interferogram and a coherence map provided by the embodiments of the present disclosure;

[0026] Figure 3 Shows a flowchart of a method for calculating a time baseline phase gradient rate set provided by the embodiments of the present disclosure;

[0027] Figure 4 Shows a flowchart of a Delaunay triangulation construction method provided by the embodiments of the present disclosure;

[0028] Figure 5 Shows a flowchart of a method for weighting Delaunay triangulation arcs provided by the embodiments of the present disclosure;

[0029] Figure 6 Shows a flowchart of a method for quasi-criterion identification of unwrapped phases provided by the embodiments of the present disclosure;

[0030] Figure 7 Shows a schematic diagram of a ground deformation field of a simulation experiment provided by the embodiments of the present disclosure;

[0031] Figure 8 Shows a schematic diagram of an original interferogram and a wrapped interferogram generated by a simulation experiment provided by the embodiments of the present disclosure;

[0032] Figure 9 Shows a schematic diagram of the phase unwrapping result of the traditional MCF method of a simulation experiment provided by the embodiments of the present disclosure;

[0033] Figure 10 Shows a schematic diagram of the phase unwrapping result of the improved MCF method guided by the phase gradient rate provided by the present disclosure in a simulation experiment;

[0034] Figure 11A schematic diagram showing the phase unwrapping result of the unwrapping phase error correction method provided by the present disclosure for an analog experiment;

[0035] Figure 12 A schematic diagram showing the structure of a phase unwrapping device based on improved minimum cost flow and quasi-identification provided by the embodiments of the present disclosure;

[0036] Figure 13 A schematic diagram showing the structure of a computer device provided by the embodiments of the present disclosure. Detailed implementation manners

[0037] To make the objectives, 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 with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are only some of the embodiments of the present disclosure, rather than all the embodiments. Components of the embodiments of the present disclosure described and illustrated in the accompanying drawings herein can be arranged and designed in various different configurations. Therefore, the detailed description of the embodiments of the present disclosure provided in the accompanying drawings is not intended to limit the scope of the present disclosure claimed, but merely represents selected embodiments of the present disclosure. All other embodiments obtained by those skilled in the art based on the embodiments of the present disclosure without creative efforts fall within the scope of protection of the present disclosure.

[0038] It should be noted that similar reference numerals and letters denote similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.

[0039] The term "and / or" in this article merely describes an association relationship and means that three relationships may exist. For example, A and / or B may represent: A exists alone, A and B exist simultaneously, and B exists alone. In addition, the term "at least one" in this article means any one of multiple or any combination of at least two of multiple. For example, including at least one of A, B, and C may represent including any one or more elements selected from the set composed of A, B, and C.

[0040] To facilitate the understanding of this embodiment, the execution subject of the phase unwrapping method based on improved minimum cost flow and quasi-identification provided by the embodiments of the present disclosure will be introduced in detail first. The execution subject of the phase unwrapping method based on improved minimum cost flow and quasi-identification provided by the embodiments of the present disclosure is a computer device. This computer device can be a server. Among them, the server can be an independent physical server, or a server cluster or distributed system composed 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.

[0041] The following will describe in detail the phase unwrapping method based on improved minimum cost flow and quasi-identification provided by the embodiments of the present application with reference to the accompanying drawings. Refer to Figure 1 As shown, it is a flowchart of a phase unwrapping method based on improved minimum cost flow and quasi-identification provided by the embodiments of the present disclosure. The method includes the following S101 to S105:

[0042] S101, obtain the SAR data set and 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.

[0043] It can be understood that SAR (Synthetic Aperture Radar) is a technology that images the ground through electromagnetic waves and has the ability to work all day and all weather, without being affected by factors such as weather and light. SAR forms a high-resolution surface image by transmitting radar signals and receiving their reflected signals, and is widely used in fields such as surface change monitoring and disaster warning. Here, the SAR data set is image data containing electromagnetic wave reflection information of the area to be monitored collected by the SAR sensor, which can provide accurate surface topography and deformation information. DEM (Digital Elevation Model) is a digital model obtained by using remote sensing technology or ground measurement technology, representing the ground elevation information of a certain area. Here, the DEM data is used to assist in understanding the terrain undulation and combined with the SAR data through the differential interferogram to further derive a more accurate phase unwrapping result.

[0044] Specifically, in SAR interferometry processing, the differential interferogram is obtained by calculating the phase difference between two SAR images, and it is an image reflecting surface deformation. It can show the changes that occur on the surface between two time points and is often used to monitor phenomena such as earthquakes, volcanic activities, and urban subsidence. The coherence map can be used to represent the similarity degree between different pixels, that is, coherence. It is obtained by comparing the phase information of adjacent pixel points at different times or different orbits and usually reflects the signal quality, noise, and data reliability.

[0045] Exemplarily, with reference to Figure 2 shown in the figure, determining the differential interferogram and coherence map based on the SAR dataset and DEM data may include the following steps S201 to S204:

[0046] S201, determining an initial SAR image set based on the SAR dataset.

[0047] Among them, the initial SAR image set includes multiple SAR images of the area to be monitored. These SAR images usually contain observation data of the same area at different times and different angles. The area to be monitored usually refers to a specific area where ground changes need to be monitored, such as high-risk areas like earthquakes, landslides, or urban changes.

[0048] S202, determining the SAR master image in the initial SAR image set, using the SAR master image as a reference, registering the other SAR images in the initial SAR image set except the SAR master image according to the DEM data, and determining a target SAR image set based on the registration result.

[0049] Specifically, the SAR master image refers to the reference image selected from the entire image set. Multiple factors such as the temporal baseline, spatial baseline, and change in Doppler center frequency need to be comprehensively considered. That is, the image that is closest in time or has the best quality is selected from the orbit dataset as the master image, while ensuring that the selected master image has a good match with other images in terms of spatial baseline and Doppler center frequency. Registration means aligning other images with the master image so that SAR images at different times and angles can be accurately matched in space. During the registration process, the DEM data is used to correct the geometric errors of the images, thereby ensuring the spatial consistency of the images. Finally, through these registered SAR images and the SAR master image, a target SAR image set is generated for subsequent analysis. Among them, the target SAR image set includes multiple target SAR images.

[0050] S203, determining a baseline set based on the target SAR image set according to a preset spatio-temporal baseline threshold.

[0051] It can be understood that the spatio-temporal baseline threshold refers to the limitation on the validity between image pairs in terms of time and space. Usually, it is a parameter value used to screen out representative and relatively consistent interferometric pairs and exclude those image pairs with unsatisfactory interferometric effects due to too long time intervals or too large spatial distances. In the present disclosure, the spatio-temporal baseline threshold is set such that the time is less than 72 days and the spatial distance is less than 100 m; in some other embodiments, it can be set according to specific needs and will not be specifically limited herein. Based on the target SAR image set, a baseline set is determined through a preset spatio-temporal baseline threshold. The baseline set includes multiple interferometric pairs, and each interferometric pair includes two target SAR images. Here, an interferometric pair refers to the pairing of two images required for interferometric analysis through SAR images. These paired images are usually separated by a certain time and space and can reflect ground changes.

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

[0053] Specifically, an interferometric pair is formed by two radar observations (or images). The time and location of each observation may be slightly different. By comparing these two observation images, an interferogram reflecting surface changes (such as subsidence, displacement, etc.) can be obtained. The differential interferogram is obtained by subtracting the topographic error (usually corrected using DEM data) and can display fine information on the relative changes of the surface. The differential interferogram of each interferometric pair can reveal the phase difference between the two observations, and this phase difference reflects information on ground displacement or deformation.

[0054] It can be understood that the coherence map reflects the signal quality and correlation of each pixel in the interferogram. Coherence essentially describes the consistency of radar signals in time and space. In radar imaging, different surface features (such as vegetation, buildings, mountains, etc.) have different reflection characteristics for radar signals, which will affect the coherence of the signals. Areas with high coherence usually represent stable surface features and consistent reflections, while areas with low coherence may be due to reasons such as surface dynamic changes, scatterer differences, or atmospheric effects. By generating a coherence map based on the differential interferogram of each interferometric pair, it is possible to effectively identify which areas have good signal quality, and then determine which areas have high data credibility and which areas may require further correction or exclusion. For example, if the coherence of a certain area is low, it may indicate that there are significant surface changes in that area or signal distortion due to atmospheric factors.

[0055] Among them, 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. Such an algorithm typically 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, it may include calculating the standard deviation of the phase difference or the coherence formula of the change in the phase gradient to calculate the coherence value of each pixel, and then obtaining a coherence map, where the value of each pixel in the map represents the degree of signal coherence at that position in the differential interferogram.

[0056] S102, calculate a set of time - baseline phase - gradient rates based on the differential interferogram and the coherence map.

[0057] It can be understood that the set of time - baseline phase - gradient rates is a set that contains the phase - gradient rates corresponding to different time baselines. Here, the "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 of the shortest time baseline, the first time baseline, and the second time baseline, which respectively represent the information of surface deformation within different time periods (i.e., the time intervals are 12 days, 24 days, and 36 days). By analyzing these different time baselines, the change characteristics of surface deformation at different time scales can be captured.

[0058] Here, the phase - gradient rate reflects the change rate of the phase difference between adjacent pixels. In the differential interferogram, the phase difference is directly related to the surface deformation. Therefore, the phase - gradient rate can indirectly reveal the distribution and rate of surface deformation. When determining the set of time - baseline phase - gradient rates, the differential interferogram can be pre - processed 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 the low - coherence region on the result; secondly, by calculating the phase difference between adjacent pixels and combining the time - baseline information, the phase - gradient rates under different time baselines can be obtained.

[0059] Specifically, the set of time - baseline phase - gradient rates will contain the phase - gradient rates corresponding to different time baselines. In this disclosure, the set of time - baseline phase - gradient rates includes the shortest - time - baseline phase - gradient rate, the first - time - baseline phase - gradient rate, and the second - time - baseline phase - gradient rate.

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

[0061] Exemplarily, referring to Figure 3As shown, when calculating the time - baseline phase - gradient rate set based on the differential interferogram and coherence map, the following steps S301 to S304 can be included:

[0062] S301. For each differential interferogram corresponding to an 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 values corresponding to each pixel.

[0063] Specifically, for each differential interferogram corresponding to an interference pair, determine the phase - gradient value of each pixel in the image. These gradient values reflect the change rate of the phase difference between adjacent pixels. The first phase - gradient map can be determined through the calculated phase - gradient values. Among them, the first phase - gradient map includes a first horizontal - direction phase - gradient map, a first vertical - direction phase - gradient map, a first diagonal - direction phase - gradient map, and a second diagonal - direction phase - gradient map.

[0064] Here, the present disclosure uses a phase - gradient calculation formula to calculate the phase - gradient value corresponding to each pixel. The formula can be expressed as:

[0065] ;

[0066] Among them, represents the phase - gradient in a specific direction of the pixel point; represents the vertical - direction phase - gradient; represents the first diagonal - direction phase - gradient; represents the horizontal - direction phase - gradient; represents the second diagonal - direction phase - gradient.

[0067] In some other embodiments, other methods can also be used to calculate the phase - gradient value corresponding to each pixel, such as methods based on frequency - domain analysis, methods based on edge detection, etc., which are not specifically limited herein.

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

[0069] 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 minute gradients caused by atmospheric errors in the phase gradient. Therefore, it is necessary to correct the phase gradient to avoid calculation errors caused by wrapping. To solve this phenomenon, the gradient wrapping result can be masked by combining with the coherence map to obtain the second phase gradient map. In this way, by retaining the effective gradient information in the regions with higher coherence and removing the invalid or incorrect gradient information in the regions with lower coherence, the calculation accuracy can be improved. Among them, the second phase gradient map includes the second horizontal direction phase gradient map, the second vertical direction phase gradient map, the third diagonal direction phase gradient map, and the fourth diagonal direction phase gradient map.

[0070] Here, the present disclosure uses a wrapping formula to implement the task of determining the gradient wrapping result. The wrapping formula can be expressed as:

[0071] ;

[0072] Among them, is expressed as the gradient wrapping result; is expressed as the phase gradient.

[0073] In some other embodiments, when determining the gradient wrapping result of each interference pair, methods such as the quality-guided method and the path-tracking method can also be used to directly process the phase map, so as to obtain a continuous phase distribution, and then calculate a phase gradient without jumps. In addition, machine learning or deep learning techniques can also be combined to predict and correct the phase wrapping problem by training a model, which is not specifically limited here.

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

[0075] Here, due to the inevitable presence of noise and errors in interferometric measurements, the phase gradient map of a single interference pair may not accurately reflect the true surface deformation information. Therefore, by stacking the phase gradient maps of multiple interference pairs, their information can be effectively fused, thereby suppressing noise, improving the signal-to-noise ratio of the data, and then more accurately reflecting the phase gradient characteristics of the surface under different time baselines.

[0076] Specifically, for different directional phase gradient atlases, based on the second phase gradient maps corresponding to each interference pair, phase gradient atlases in different directions are determined. Among them, the first-direction phase gradient atlas includes first-direction phase gradient stack maps corresponding to each time baseline, the second-direction phase gradient atlas includes second-direction phase gradient stack maps corresponding to each time baseline, the third-direction phase gradient atlas includes third-direction phase gradient stack maps corresponding to each time baseline, and the fourth-direction phase gradient atlas includes fourth-direction phase gradient stack maps corresponding to each time baseline. These stack maps can intuitively display the phase gradient change along a certain direction at different time baselines, which helps to analyze and identify the phase change characteristics at different time baselines in this direction.

[0077] Here, when the present disclosure uses a stacking formula to stack the phase gradient maps corresponding to each interference pair, the formula can be expressed as:

[0078] ;

[0079] Among them, represents the phase gradient stacking result in the k direction, where k takes values of 0, 45, 90, 135; m represents the number of interference pairs participating in the phase gradient stacking calculation; represents the time baseline of the differential interferogram corresponding to the nth interference pair.

[0080] S304. Based on the first-direction phase gradient stack maps, second-direction phase gradient stack maps, third-direction phase gradient stack maps, and fourth-direction phase gradient stack maps corresponding to each time baseline, determine the time baseline phase gradient rate set.

[0081] It can be understood that by stacking the phase gradient atlases in each direction corresponding to each time baseline and calculating the change rate of these stack maps 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 the change speed of the phase gradient over time.

[0082] In some possible embodiments, in order to more accurately determine the time baseline phase gradient rate set, it can be further refined on the original basis by introducing steps (1) to (3) for determining the phase gradient rate for a specific time baseline (such as the shortest time baseline, the first time baseline, the second time baseline, etc.):

[0083] (1) 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 stack diagram corresponding to the shortest time baseline, the second - direction phase - gradient stack diagram, the third - direction phase - gradient stack diagram, and the fourth - direction phase - gradient stack diagram;

[0084] (2) Determine the first - time - baseline phase - gradient rate corresponding to the first time baseline based on the gradient fusion formula, the first - direction phase - gradient stack diagram corresponding to the first time baseline, the second - direction phase - gradient stack diagram, the third - direction phase - gradient stack diagram, and the fourth - direction phase - gradient stack diagram;

[0085] (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 stack diagram corresponding to the second time baseline, the second - direction phase - gradient stack diagram, the third - direction phase - gradient stack diagram, and the fourth - direction phase - gradient stack diagram.

[0086] Here, due to the shape of the surface deformation field, the stacked results of the phase gradients in four directions have different shapes and need to be fused to obtain the complete gradient - rate information of the surface deformation. When performing phase - gradient fusion, to avoid the gradients in different directions from canceling each other out, the method of taking the mean of absolute values is adopted. That is, the gradient fusion formula can be expressed as:

[0087] ;

[0088] where \(G\) represents the fusion result of the stacked diagrams of the phase gradients in four directions; represents the stacked diagram of the phase gradient in the \(k\) - direction, and the value of \(k\) is 0, 45, 90, 135, represents the first direction, represents the second direction, represents the third direction, represents the fourth direction.

[0089] It can be understood that the shortest time - baseline phase - gradient rate provides an initial reference for subsequent analysis regarding surface deformation or terrain changes. By calculating the phase - gradient rates corresponding to each time baseline respectively, the rate of surface deformation or terrain changes within a specific time period can be obtained.

[0090] Here, the shortest time baseline mentioned in the present 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 some other embodiments, if specific requirements are met, further calculations can be extended to consider longer time baselines (such as 48 days, 60 days, etc.) to provide a wider analysis scope. By flexibly adjusting the length of the time baseline, different analysis requirements can be adapted to ensure more targeted analysis results can be provided.

[0091] S103. Construct a Delaunay triangulation network based on the differential interferogram and the coherence map, and weight each triangular arc segment in the Delaunay triangulation network based on the coherence map and the time baseline phase gradient rate set to obtain the weight values corresponding to each triangular arc segment in the Delaunay triangulation network.

[0092] It can be understood that the Delaunay triangulation network is a method of connecting a set of points on a plane into a triangular grid, and these triangles satisfy the empty circumcircle property, that is, no other points are contained within the circumcircle of any triangle. In geographic information systems and remote sensing data analysis, the Delaunay triangulation network 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 network is used to divide the area to be monitored into multiple triangular regions and is further used to define the constraints for phase unwrapping.

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

[0094] Exemplarily, with reference to Figure 4 as shown, when constructing the Delaunay triangulation network, the following steps S401 to S402 may be included:

[0095] S401. For each differential interferogram corresponding to an interference pair, extract high-quality points that meet a preset coherence threshold based on the coherence map corresponding to the differential interferogram, and determine a high-quality point set corresponding to the interference pair based on the extraction result.

[0096] Specifically, high-quality points that meet a preset coherence threshold are extracted through a coherence map. These points represent ground object information with high coherence within the monitoring area. By selecting high-quality points, it is possible to ensure that the point set used in the subsequent Delaunay triangulation process has high accuracy and representativeness, guaranteeing the quality of points in the Delaunay triangulation. In this disclosure, the preset coherence threshold is set to 0.4. In some other embodiments, it can be adjusted according to the actual situation and is not specifically limited herein.

[0097] S402. For each set of high-quality points corresponding to an interference pair, use the point-by-point insertion method to generate a Delaunay triangulation corresponding to the interference pair.

[0098] It can be understood that after extracting the set of high-quality points corresponding to each interference pair, a Delaunay triangulation can be generated based on these extracted high-quality points using the point-by-point insertion method. Each Delaunay triangulation includes multiple triangles. The point-by-point insertion method is a classic mesh generation method that forms a mesh meeting the Delaunay condition by adding new points one by one and adjusting the existing triangle structure.

[0099] In some possible embodiments, to improve the accuracy and efficiency of surface deformation analysis, especially when dealing with complex terrains and various interference factors, during the process of determining weights for each triangular arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set, the weight allocation process can be further optimized. Referring to Figure 5 as shown, it can include the following steps S501 to S503:

[0100] S501. For each Delaunay triangulation corresponding to an interference pair, determine whether there is a triangular arc segment corresponding to the shortest time baseline in the Delaunay triangulation.

[0101] Specifically, check each Delaunay triangulation generated by different interference pairs to identify whether there is a triangular arc segment defined by the shortest time baseline (i.e., the interference pair with the shortest time interval between two observations). The shortest time baseline usually means higher quality of the interferogram because a shorter time interval reduces the influence of factors other than surface changes (such as atmospheric interference, vegetation growth, etc.) on interferometry.

[0102] S502. If it exists, set the weight value of the triangular arc segment corresponding to the shortest time baseline to a preset coherence value.

[0103] Here, if the triangular arcs defined by the shortest time baseline are found, these arcs 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 arcs in deformation analysis. In this disclosure, the preset coherence value is set to 1.

[0104] S503. If not, weight the triangular arcs in the Delaunay triangulation based on the arc weight calculation formula, the coherence map, and the set of time baseline phase gradient rates.

[0105] Here, if there are no triangular arcs defined by the shortest time baseline in the Delaunay triangulation, or for a more comprehensive weight assignment to all arcs, an arc weight calculation formula can be used to weight them. Among them, the arc weight calculation formula can be expressed as:

[0106] ;

[0107] Among them, represents the weight of arc ; The starting position of arc is and the ending position is ; represents the phase gradient rate at the starting position of arc ; represents the phase gradient rate at the ending position of arc ; represents the phase gradient rate threshold, and its value is , represents the phase gradient rate; represents the coherence value at the starting position of arc ; represents the coherence value at the ending position of arc ;

[0108] S104. Complete 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 arcs in the Delaunay triangulation, and obtain the initial phase unwrapping result.

[0109] It can be understood that the Minimum Cost Flow Method (MCF) is an optimization algorithm in graph theory. In a network, each edge has a weight or cost, and the goal of the algorithm is to find the minimum-cost path or flow allocation scheme from the source node to the sink node in a given network. In the context of phase unwrapping, the minimum cost flow method is used to find an allocation scheme of phase values such that the cumulative cost of the phase difference between adjacent pixels (i.e., the phase gradient of the arc segment) is minimized.

[0110] Specifically, after weighting each triangle arc segment in the Delaunay triangulation network, the minimum cost flow algorithm is used to optimize the flow allocation in this network, decouple the phase information, and thus effectively complete phase unwrapping to obtain the initial phase unwrapping result, where the initial phase unwrapping result includes the unwrapped phases corresponding to multiple temporal baselines. Here, by guiding the unwrapping network through the phase gradient rate, the phase unwrapping results of most interferograms can be accurately restored.

[0111] S105. Quasi-identification is respectively performed on the unwrapped phases corresponding to each temporal baseline, and the initial phase unwrapping result is adjusted based on the quasi-identification result to obtain the target phase unwrapping result.

[0112] It can be understood that during the process of phase unwrapping, in order to ensure the accurate recovery of deformation information, especially for regions where the gradients at the boundaries of the deformation field are all wrapped, the complexity of the deformation field boundaries may cause phase gradients to be confused in these regions. Simply relying on the minimum cost flow method often cannot provide accurate unwrapping results.

[0113] In response to this, the present disclosure proposes to respectively perform quasi-identification on the unwrapped phases corresponding to each temporal baseline, and adjust the initial phase unwrapping result based on the quasi-identification result to obtain the target phase unwrapping result. Here, quasi-identification refers to a process of verifying and calibrating the unwrapped phase data, aiming to evaluate its accuracy and reliability. By performing quasi-identification on the unwrapped phases corresponding to each temporal baseline, the phase changes of the entire deformation field at different time points can be comprehensively evaluated, so as to more accurately capture the deformation information. Based on the quasi-identification, the initial phase unwrapping result is adjusted, and the adjustment process may involve correcting the phase values, re-unwrapping the error regions, or smoothing the overall phase field to ensure that the deformation information in all regions can be accurately restored to obtain the target phase unwrapping result.

[0114] Exemplarily, referring to Figure 6 As shown, when respectively performing quasi-identification on the unwrapped phases corresponding to each temporal baseline and adjusting the initial phase unwrapping result based on the quasi-identification result, the following steps S601~S605 may be included:

[0115] S601. Determine the first matrix to be quasi-fitted based on the initial phase unwrapping result.

[0116] 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 result, and together they form the first matrix to be quasi-fitted. The shortest time baseline usually refers to the shortest interval between two adjacent observation times in the time series. Since the time interval is short and the deformation change is small, the unwrapped phase corresponding to it can be considered relatively stable, that is, the "quasi-stable phase". The unwrapped phase corresponding to the first time baseline may have certain uncertainties due to the complexity of the deformation field or observation errors, that is, the "non-quasi-stable phase".

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

[0118] Specifically, for the first matrix to be quasi-fitted, 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, use the gross error observation equation to calculate the first gross error estimate, which reflects the 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:

[0119] ;

[0120] ;

[0121] Among them, V represents the observation residual; A represents an n×m-dimensional coefficient matrix, n represents the number of re-unwrapped phases of the matrix to be quasi-fitted, and m represents the number of SAR images in the initial SAR image set; represents the hyperparameter; represents the gross error estimate; L represents the matrix to be quasi-fitted; represents the matrix composed of multiple non-quasi-stable phases; P represents the identity matrix; R represents the adjustment factor.

[0122] 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 the second matrix to be quasi-fitted based on the adjusted initial phase unwrapping result.

[0123] It can be understood that after obtaining 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 (such as rounding) the obtained first gross error estimate to an integer multiple of the whole cycle, and then removing the constrained gross error from the corresponding phase.

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

[0125] S604. For the second matrix to be quasi-aligned, mark the unwrapped phase corresponding to the shortest time baseline and the adjusted unwrapped phase corresponding to the first time baseline as the second quasi-stable phase, and mark the unwrapped phase corresponding to the second time baseline as the 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.

[0126] Specifically, for the second matrix to be quasi-aligned, the unwrapped phase corresponding to the shortest time baseline and the adjusted unwrapped phase corresponding to the 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 used again to calculate the second gross error estimate, which reflects the degree of difference between the unwrapped phase corresponding to the second time baseline and the quasi-stable phase.

[0127] S605. Adjust the unwrapped 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.

[0128] It can be understood 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 to obtain a more accurate target phase unwrapping result.

[0129] 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 have been quasi-aligned and adjusted, so as to ensure that the deformation information of the entire deformation field can be accurately restored.

[0130] In the phase unwrapping method and device based on improved minimum cost flow and quasi-identification provided in the embodiments of the present disclosure, during the phase unwrapping process using the improved minimum cost flow method, by optimizing the calculation of the phase gradient rate set and the weight determination of the triangular arc segments, the error accumulation is effectively reduced, thereby ensuring a high-precision unwrapping result and accelerating the unwrapping process. At the same time, the present disclosure also uses the quasi-identification technology to post-process the initial phase unwrapping result, further improving the phase unwrapping accuracy, and then realizing high-precision surface deformation monitoring.

[0131] To verify the performance of the phase unwrapping method proposed in the present disclosure, the present disclosure simulates three common types of surface deformation rate fields (as shown in Figure 7 (a) and (b)) to analyze the influence 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 of the three are different. Among them, the deformation of deformation field one increases slowly from 0 to 25 cm / a from top to bottom, and there are obvious deformation boundaries between the two sides and the lower side and the stable area, while there is no obvious deformation boundary on the upper side; there is a significant deformation boundary between the deformation of deformation field two and the stable area, and the deformation rate within the area is 25 cm / a. The deformation of deformation field three increases slowly from 0 to 25 cm / a from the boundary to the center, and there is no obvious deformation boundary compared with the surrounding stable area; for the convenience of distinguishing this type of deformation rate field, their shapes are respectively set as trapezoidal deformation field, rectangular deformation field and circular deformation field in this paper. Then, the present disclosure obtains the orbital information of 30 scenes of Sentinel-1 descending orbit images from January 2, 2019 to December 28, 2019. And according to the time baseline threshold of 72 days and the space baseline threshold of 100 meters, 153 interference pairs are generated. 153 original interferograms are generated according to the spatio-temporal baseline, and then the phase of the original interferogram is wrapped to to obtain the wrapped interferogram. When generating the simulated interferogram, periodic deformation is added on the basis of linear deformation. The first six interferometric phase diagrams generated by simulation are as shown in Figure 8 (a)-(f) are the simulated original interferometric 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 interferograms obtained by wrapping (a)-(f). It can be seen from Figure 8 that the phase of the wrapped interferogram is between , and when the time baseline is greater than 24 days, phase jumps occur inside the circular deformation field, inside and at the boundary of the trapezoidal deformation field, and at the boundary of the rectangular deformation field.

[0132] After generating the simulated wrapped phase diagram, the present disclosure performs phase unwrapping through the traditional MCF method and the MCF method based on the phase gradient rate respectively. The unwrapping result of the traditional MCF method is as Figure 9 shown, Figure 9 . (a) is the unwrapped phase diagram obtained from the wrapped phase diagram with a time baseline of 12 days through the MCF method. Similarly, Figure 9 . (b)-9.(f) are the unwrapped phase diagrams obtained from the wrapped phase diagrams of 24 days, 36 days, 48 days, 60 days, and 72 days respectively through the MCF method. Figure 9 . (g)-(l) are the differences between the unwrapped phase diagrams and the original phase diagrams with the same time baseline. From Figure 9 . (a), it can be seen that when the time baseline is 12 days, the phase is correctly restored without unwrapping errors. In Figure 9 . (g), there are only errors caused by the computer storage precision. When the time baseline is greater than 12 days, the circular deformation area obtains the correct unwrapping result, while the trapezoidal deformation area and the rectangular deformation area obtain incorrect unwrapping results.

[0133] Specifically, the unwrapping result of the MCF method guided by the phase gradient rate is as Figure 10 shown, Figure 10 . (a) is the unwrapped phase diagram obtained from the wrapped phase diagram with a time baseline of 12 days through the proposed phase unwrapping method. Similarly, Figure 10 . (b)-10.(f) are the unwrapped phase diagrams obtained by phase unwrapping of the wrapped phase diagrams of 24 days, 36 days, 48 days, 60 days, and 72 days respectively; Figure 10 . (g)-(l) are the differences between the unwrapped phase diagrams and the simulated original phase diagrams with the same time baseline. When there are unwrapped regions at the deformation field boundary, the correct unwrapping result is obtained. However, with the time cumulative effect of the surface deformation, the deformation gradient at the deformation boundary gradually becomes larger, and the MCF method based on the phase gradient rate obtains an incorrect phase unwrapping result. Then, the result after correcting the unwrapping phase error (i.e., quasi-identification adjustment) is as Figure 11 shown, Figure 11 . (a)-11.(f) are the results after correcting the unwrapped phases of the wrapped phase diagrams of 12 days, 24 days, 36 days, 48 days, 60 days, and 72 days respectively; Figure 11 . (g)-11.(l) are the differences between the unwrapped phase diagrams and the simulated original phase diagrams with the same time baseline. It can be seen from the figure that the unwrapping errors are basically all corrected.

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

[0135] Based on the same inventive concept, an embodiment of the present disclosure further provides a phase unwrapping device based on improved minimum cost flow and quasi-criterion identification corresponding to the phase unwrapping method based on improved minimum cost flow and quasi-criterion identification. Since the principle of solving problems by the device in the embodiment of the present disclosure is similar to that of the phase unwrapping method based on improved minimum cost flow and quasi-criterion identification in the above 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 elaborated.

[0136] Referring to Figure 12 As shown, it is a schematic diagram of a phase unwrapping device 1200 based on improved minimum cost flow and quasi-criterion identification provided by an embodiment of the present disclosure. The device includes:

[0137] A data acquisition module 1201, configured to acquire an SAR data set and DEM data of an area to be monitored; and generate a differential interferogram and a coherence map based on the SAR data set and the DEM data;

[0138] A rate calculation module 1202, configured to calculate a time baseline phase gradient rate set based on the differential interferogram and the coherence map;

[0139] A weight determination module 1203, configured to construct a Delaunay triangulation network based on the differential interferogram and the coherence map, and determine weights for each triangular arc segment in the Delaunay triangulation network based on the coherence map and the time baseline phase gradient rate set, to obtain weights corresponding to each triangular arc segment in the Delaunay triangulation network;

[0140] A phase unwrapping module 1204, configured to complete the initial phase unwrapping task of the area to be monitored based on the minimum cost flow method and the weights corresponding to each triangular arc segment in the Delaunay triangulation network, to obtain an initial phase unwrapping result; wherein, the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines;

[0141] A phase quasi-criterion module 1205, configured to perform quasi-criterion identification on the unwrapped phases corresponding to each time baseline respectively, and adjust the initial phase unwrapping result based on the quasi-criterion identification result, to obtain a target phase unwrapping result.

[0142] In some possible embodiments, the data acquisition module 1201 is specifically configured to:

[0143] Determine an initial SAR image set based on the SAR data set; wherein, the initial SAR image set includes multiple SAR images of the area to be monitored;

[0144] Determine the SAR master image in the initial SAR image set. Based on the SAR master image, register the other SAR images in the initial SAR image set except the SAR master image according to the DEM data, and determine the target SAR image set based on the registration result; wherein, the target SAR image set includes multiple target SAR images.

[0145] Determine a baseline set based on the target SAR image set according to a preset spatio-temporal baseline threshold; wherein, the baseline set includes multiple interference pairs, and each interference pair includes two target SAR images.

[0146] Generate a differential interferogram corresponding to each interference pair based on the DEM data; and generate a coherence map corresponding to each interference pair based on the differential interferogram corresponding to each interference pair.

[0147] In some possible embodiments, the rate calculation module 1202 is specifically configured to:

[0148] 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 direction phase gradient map, a first vertical direction phase gradient map, a first diagonal direction phase gradient map, and a second diagonal direction phase gradient map.

[0149] For each interference pair, determine a gradient wrapping result according to the first phase gradient map corresponding to the differential interferogram, and mask the gradient wrapping result based on the coherence map corresponding to the differential interferogram 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.

[0150] Determine the first-direction phase gradient map set based on the second horizontal-direction phase gradient maps corresponding to the interference pairs in the baseline set, determine the second-direction phase gradient map set based on the second vertical-direction phase gradient maps corresponding to the interference pairs in the baseline set, determine the third-direction phase gradient map set based on the third diagonal-direction phase gradient maps corresponding to the interference pairs in the baseline set, and determine the fourth-direction phase gradient map set based on the fourth diagonal-direction phase gradient maps corresponding to the interference pairs in the baseline set; wherein, the first-direction phase gradient map set includes first-direction phase gradient stack maps corresponding to each time baseline, the second-direction phase gradient map set includes second-direction phase gradient stack maps corresponding to each time baseline, the third-direction phase gradient map set includes third-direction phase gradient stack maps corresponding to each time baseline, and the fourth-direction phase gradient map set includes fourth-direction phase gradient stack maps corresponding to each time baseline;

[0151] Determine the time baseline phase gradient rate set based on the first-direction phase gradient stack maps, second-direction phase gradient stack maps, third-direction phase gradient stack maps, and fourth-direction phase gradient stack maps corresponding to each time baseline.

[0152] In some possible embodiments, the time baseline includes the shortest time baseline, the first time baseline, and the second time baseline, and 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; the rate calculation module 1202 is specifically configured to:

[0153] 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 stack map, second-direction phase gradient stack map, third-direction phase gradient stack map, and fourth-direction phase gradient stack map corresponding to the shortest time baseline;

[0154] Determine the first time baseline phase gradient rate corresponding to the first time baseline based on the gradient fusion formula, the first-direction phase gradient stack map, second-direction phase gradient stack map, third-direction phase gradient stack map, and fourth-direction phase gradient stack map corresponding to the first time baseline;

[0155] 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 stack map, second-direction phase gradient stack map, third-direction phase gradient stack map, and fourth-direction phase gradient stack map corresponding to the second time baseline;

[0156] The gradient fusion formula includes:

[0157] ;

[0158] Among them, G represents the fusion result of the four-direction phase gradient stacked images; represents the phase gradient stacked image in the k direction, and the value of k is 0, 45, 90, 135, represents the first direction, represents the second direction, represents the third direction, represents the fourth direction.

[0159] In some possible embodiments, the weight determination module 1203 is specifically configured to:

[0160] For each differential interferogram corresponding to an interference pair, extract high-quality points that meet a preset coherence threshold based on the coherence map corresponding to the differential interferogram, and determine a high-quality point set corresponding to the interference pair based on the extraction result;

[0161] For each high-quality point set corresponding to an interference pair, use the point-by-point insertion method to generate a Delaunay triangulation corresponding to the interference pair; wherein, each Delaunay triangulation includes a plurality of triangles;

[0162] The weight determination module 1203 is specifically configured to:

[0163] For each Delaunay triangulation corresponding to an interference pair, determine whether there is a triangle arc segment corresponding to the shortest temporal baseline in the Delaunay triangulation;

[0164] If it exists, set the weight of the triangle arc segment corresponding to the shortest temporal baseline to a preset coherence value;

[0165] If it does not exist, weight each triangle arc segment in the Delaunay triangulation based on the arc segment weight calculation formula, the coherence map, and the temporal baseline phase gradient rate set;

[0166] The arc segment weight calculation formula includes:

[0167] ;

[0168] Among them, represents the weight of the arc segment ; the starting position of the arc segment is and the ending position is ; represents the phase gradient rate at the starting position of the arc segment ; represents the arc segment The phase gradient rate at the end position; Denoted as the phase gradient rate threshold, whose value is , Denoted as the phase gradient rate; Denoted as the coherence value at the starting position of the arc segment ; Denoted as the arc segment The coherence value at the end position.

[0169] In some possible embodiments, the phase alignment module 1205 is specifically configured to:

[0170] Determine a first matrix to be aligned based on the initial phase unwrapping result; wherein, the first matrix to be aligned includes the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the first time baseline;

[0171] For the first matrix to be aligned, mark the unwrapped phase corresponding to the shortest time baseline as the first stable phase, and mark the unwrapped phase corresponding to the first time baseline as the first non-stable phase; and determine a first gross error estimate based on the gross error observation equation, the first stable phase, and the first non-stable phase;

[0172] 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 aligned based on the adjusted initial phase unwrapping result; wherein, the second matrix to be aligned includes the unwrapped phase corresponding to the shortest time baseline, the adjusted unwrapped phase corresponding to the first time baseline, and the unwrapped phase corresponding to the second time baseline;

[0173] For the second matrix to be aligned, mark the unwrapped phase corresponding to the shortest time baseline and the adjusted unwrapped phase corresponding to the first time baseline as the second stable phase, and mark the unwrapped phase corresponding to the second time baseline as the second non-stable phase; and determine a second gross error estimate based on the gross error observation equation, the second stable phase, and the second non-stable phase;

[0174] Adjust the unwrapped 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.

[0175] In some possible embodiments, the gross error observation equation includes:

[0176] ;

[0177] ;

[0178] Wherein, V represents the observation residual; A represents an n×m dimensional coefficient matrix, n represents the number of unwrapped phases of the matrix to be quasi-determined, and m represents the number of SAR images in the initial SAR image set; is represented as a hyperparameter; is represented as the gross error estimate; L represents the matrix to be quasi-determined; is represented as a matrix composed of multiple non-stable phases; P represents the identity matrix; R represents the adjustment factor.

[0179] Based on the same inventive concept, an embodiment of the present disclosure also provides a computer device. Referring to Figure 13 As shown, it is a schematic structural diagram of a computer device 1300 provided by an embodiment of the present disclosure, including a processor 1301, a memory 1302, and a bus 1303. Among them, the memory 1302 is used to store execution instructions, including an internal memory 13021 and an external memory 13022; here, the internal memory 13021 is also called the main memory, which is used to temporarily store the operation data in the processor 1301 and the data exchanged with the external memory 13022 such as a hard disk. The processor 1301 exchanges data with the external memory 13022 through the internal memory 13021.

[0180] In an embodiment of the present application, the memory 1302 is specifically used to store the application program code for implementing the solution of the present application, and is controlled by the processor 1301 to execute. That is, when the computer device 1300 runs, the processor 1301 communicates with the memory 1302 through the bus 1303, so that the processor 1301 executes the application program code stored in the memory 1302, and further executes the method described in any of the foregoing embodiments. Among them, the memory 1302 may be, but is not limited to, a random access memory, a read-only memory, a programmable read-only memory, an erasable read-only memory, an electrically erasable read-only memory, etc. The processor 1301 may be an integrated circuit chip with signal processing capabilities. The above-mentioned processor may be a general-purpose processor, including a central processing unit, a network processor, etc.; it may also be a digital signal processor, an application-specific integrated circuit, a field-programmable gate array, or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components. It can implement or execute the various methods, steps, and logic block diagrams disclosed in the embodiments of the present invention. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc.

[0181] It can be understood that the structure schematically shown in the embodiment 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 combine certain components, or split certain components, or different component arrangements. The components shown may be implemented in hardware, software, or a combination of software and hardware.

[0182] An embodiment of the present disclosure also provides a computer-readable storage medium, on which a computer program is stored. When the computer program is run by a processor, it executes the steps of the phase unwrapping method based on improved minimum cost flow and quasi-identification described in the above method embodiment. Among them, the storage medium may be a volatile or non-volatile computer-readable storage medium.

[0183] An embodiment of the present disclosure also provides a computer program product, which carries 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 and will not be elaborated here.

[0184] Among them, the above computer program product can be specifically implemented by means of hardware, software, or a combination thereof. In an optional embodiment, the computer program product is specifically embodied as a computer storage medium. In another optional embodiment, the computer program product is specifically embodied as a software product, such as a software development kit, etc.

[0185] Those skilled in the art can clearly understand that for the convenience and brevity of description, the specific working processes of the systems and devices described above can refer to the corresponding processes in the foregoing method embodiments and will not be elaborated herein. In several embodiments provided in the present 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 is only a logical function division, and there may be other division methods in actual implementation. For another example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the couplings or direct couplings or communication connections shown or discussed with each other can be through some communication interfaces. The indirect couplings or communication connections of the devices or units can be in electrical, mechanical or other forms. The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they can be located in one place or distributed to multiple network units. Some or all of the units can be selected according to actual needs to achieve the purpose of the solution of this embodiment. In addition, in each embodiment of the present disclosure, the functional units can be integrated into one processing unit, or each unit can exist physically alone, or two or more units can be integrated into one unit. If the functions are implemented in the form of software function units and sold or used as independent products, they can be stored in a non-volatile computer-readable storage medium executable by a processor. Based on such an understanding, the technical solution of the present disclosure, in essence, or the part that contributes to the prior art or a part of this 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 causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in each embodiment of the present disclosure. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories, random access memories, magnetic disks, or optical discs that can store program codes.

[0186] Finally, it should be noted that the above-described embodiments are only specific implementation manners of the present disclosure, used to illustrate the technical solutions of the present disclosure, rather than limiting it. The protection scope of the present disclosure is not limited thereto. Although the present disclosure has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that any person skilled in the art within the technical scope disclosed by the present disclosure can still modify the technical solutions recorded in the foregoing embodiments or can easily think of changes, or perform equivalent replacements on some of the technical features; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present disclosure, and should all be covered within the protection scope of the present disclosure. Therefore, the protection scope of the present disclosure should be subject to the protection scope of the claims.

Claims

1. A phase unwrapping method based on improved minimum cost flow and quasi-criterion identification, characterized in that Including: Obtain the SAR data set and DEM data of the area to be monitored; Generate a differential interferogram and a coherence map based on the SAR data set and the DEM data; Calculate a set of time - baseline phase gradient rates based on the differential interferogram and the coherence map; Construct a Delaunay triangulation network based on the differential interferogram and the coherence map, and weight each triangular arc segment in the Delaunay triangulation network based on the coherence map and the set of time - baseline phase gradient rates to obtain the weights corresponding to each triangular arc segment in the Delaunay triangulation network; Complete the initial phase unwrapping task of the area to be monitored based on the minimum - cost flow method and the weights corresponding to each triangular arc segment in the Delaunay triangulation network to obtain an initial phase unwrapping result; wherein, the initial phase unwrapping result includes unwrapped phases corresponding to multiple time baselines, and the time baselines include the shortest time baseline, the first time baseline, and the second time baseline; Determine a first matrix to be quasi - calibrated based on the initial phase unwrapping result; wherein, the first matrix to be quasi - calibrated includes the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the first time baseline; For the first matrix to be quasi - calibrated, mark the unwrapped phase corresponding to the shortest time baseline as the first quasi - stable phase, and mark the unwrapped phase corresponding to the first time baseline as the 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; 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 quasi - calibrated based on the adjusted initial phase unwrapping result; wherein, the second matrix to be quasi - calibrated includes the unwrapped phase corresponding to the shortest time baseline, the adjusted unwrapped phase corresponding to the first time baseline, and the unwrapped phase corresponding to the second time baseline; For the second matrix to be quasi - calibrated, mark the unwrapped phase corresponding to the shortest time baseline and the adjusted unwrapped phase corresponding to the first time baseline as the second quasi - stable phase, and mark the unwrapped phase corresponding to the second time baseline as the second non - quasi - stable phase; and determine a second gross error estimate based on the gross error observation equation, the second quasi - stable phase, and the second non - quasi - stable phase; Adjust the unwrapped phase corresponding to the second time baseline in the adjusted initial phase unwrapping result based on the second gross error estimate to obtain a target phase unwrapping result.

2. The phase unwrapping method based on improved minimum cost flow and quasi-criterion identification according to claim 1, wherein The generating a differential interferogram and a coherence map based on the SAR data set and the DEM data includes: Determine an initial SAR image set based on the SAR data set; wherein, the initial SAR image set includes multiple SAR images of the area to be monitored; Determine the SAR master image in the initial SAR image set. Based on the SAR master image, register the other SAR images in the initial SAR image set except the SAR master image according to the DEM data, and determine the target SAR image set based on the registration result; wherein, the target SAR image set includes multiple target SAR images. Determine a baseline set based on the target SAR image set according to a preset spatio-temporal baseline threshold; wherein, the baseline set includes multiple interferometric pairs, and each interferometric pair includes two target SAR images. Generate a differential interferogram corresponding to each interferometric pair based on the DEM data; and generate a coherence map corresponding to each interferometric pair based on the differential interferogram corresponding to each interferometric pair.

3. The phase unwrapping method based on improved minimum cost flow and quasi-criterion identification according to claim 2, wherein The calculating the time baseline phase gradient rate set based on the differential interferogram and the coherence map includes: For the differential interferogram corresponding to each interferometric 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 values corresponding to each pixel; wherein, the first phase gradient map includes a first horizontal direction phase gradient map, a first vertical direction phase gradient map, a first diagonal direction phase gradient map, and a second diagonal direction phase gradient map. For each interferometric pair, determine a gradient wrapping result according to the first phase gradient map corresponding to the differential interferogram, and mask the gradient wrapping result based on the coherence map corresponding to the differential interferogram 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. Determine a first direction phase gradient map set based on the second horizontal direction phase gradient maps corresponding to the interferometric pairs in the baseline set, determine a second direction phase gradient map set based on the second vertical direction phase gradient maps corresponding to the interferometric pairs in the baseline set, determine a third direction phase gradient map set based on the third diagonal direction phase gradient maps corresponding to the interferometric pairs in the baseline set, and determine a fourth direction phase gradient map set based on the fourth diagonal direction phase gradient maps corresponding to the interferometric pairs in the baseline set; wherein, the first direction phase gradient map set includes a first direction phase gradient stack map corresponding to each time baseline, the second direction phase gradient map set includes a second direction phase gradient stack map corresponding to each time baseline, the third direction phase gradient map set includes a third direction phase gradient stack map corresponding to each time baseline, and the fourth direction phase gradient map set includes a fourth direction phase gradient stack map corresponding to each time baseline. Determine the time baseline phase gradient rate set based on the first direction phase gradient stack map, the second direction phase gradient stack map, the third direction phase gradient stack map, and the fourth direction phase gradient stack map corresponding to each time baseline.

4. The phase unwrapping method based on improved minimum cost flow and quasi-criterion identification according to claim 3, characterized in that 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; determining the time - baseline phase - gradient rate set based on the first - direction phase - gradient stack map, the second - direction phase - gradient stack map, the third - direction phase - gradient stack map, and the fourth - direction phase - gradient stack map corresponding to each time - baseline includes: Determining the shortest - time - baseline phase - gradient rate corresponding to the shortest time - baseline based on the gradient - fusion formula, the first - direction phase - gradient stack map, the second - direction phase - gradient stack map, the third - direction phase - gradient stack map, and the fourth - direction phase - gradient stack map corresponding to the shortest time - baseline; Determining the first - time - baseline phase - gradient rate corresponding to the first time - baseline based on the gradient - fusion formula, the first - direction phase - gradient stack map, the second - direction phase - gradient stack map, the third - direction phase - gradient stack map, and the fourth - direction phase - gradient stack map corresponding to the first time - baseline; Determining the second - time - baseline phase - gradient rate corresponding to the second time - baseline based on the gradient - fusion formula, the first - direction phase - gradient stack map, the second - direction phase - gradient stack map, the third - direction phase - gradient stack map, and the fourth - direction phase - gradient stack map corresponding to the second time - baseline; The gradient - fusion formula includes: ； Among them, G represents the fusion result of the four-direction phase gradient stacked images; represents the phase gradient stacked image in the k direction, where k takes values of 0, 45, 90, and 135, represents the first direction, represents the second direction, represents the third direction, represents the fourth direction.

5. The phase unwrapping method based on improved minimum cost flow and quasi-criterion identification according to claim 4, wherein Constructing a Delaunay triangulation based on the differential interferogram and the coherence map includes: For the differential interferogram corresponding to each interferometric pair, extracting high - quality points that meet a preset coherence threshold based on the coherence map corresponding to the differential interferogram, and determining a high - quality point set corresponding to the interferometric pair based on the extraction result; For the high - quality point set corresponding to each interferometric pair, generating a Delaunay triangulation corresponding to the interferometric pair using the point - by - point insertion method; wherein each Delaunay triangulation includes a plurality of triangles; Weight - setting for each triangular arc segment in the Delaunay triangulation based on the coherence map and the time - baseline phase - gradient rate set includes: For the Delaunay triangulation corresponding to each interferometric pair, determining whether there is a triangular arc segment corresponding to the shortest time - baseline in the Delaunay triangulation; If it exists, setting the weight value of the triangular arc segment corresponding to the shortest time - baseline to a preset coherence value; If it does not exist, weighting each triangular arc segment in the Delaunay triangulation based on the arc - segment weight - calculation formula, the coherence map, and the time - baseline phase - gradient rate set; The arc - segment weight - calculation formula includes: ； Among them, represents the weight value of the arc segment ; the starting position of the arc segment is , and the ending position is ; represents the phase gradient rate of the starting position of the arc segment ; represents the phase gradient rate of the ending position of the arc segment ; represents the phase gradient rate threshold, and its value is , represents the phase gradient rate; represents the coherence value of the starting position of the arc segment ; represents the coherence value of the ending position of the arc segment .

6. The phase unwrapping method based on improved minimum cost flow and quasi-criterion identification according to claim 1, wherein The gross - error observation equation includes: ； ； Among them, V represents the observation residual; A represents an n×m dimensional coefficient matrix, where n represents the number of phases to be re-unwrapped for the matrix to be quasi-determined, and m represents the number of SAR images in the initial SAR image set; is represented as a hyperparameter; is represented as the gross error estimate; L represents the matrix to be quasi-determined; is represented as a matrix composed of multiple non-stable phases; P represents the identity matrix; R represents the adjustment factor.

7. A phase unwrapping device based on improved minimum cost flow and quasi-criterion identification, characterized in that, Includes: A data acquisition module for acquiring an 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; A rate calculation module for calculating a time - baseline phase - gradient rate set based on the differential interferogram and the coherence map; A weight determination module, configured to construct a Delaunay triangulation based on the differential interferogram and the coherence map, and determine weights for each triangular arc segment in the Delaunay triangulation based on the coherence map and the time baseline phase gradient rate set, so as to obtain weights corresponding to each triangular arc segment in the Delaunay triangulation; A phase unwrapping module, configured to complete the initial phase unwrapping task of the area to be monitored based on the minimum cost flow method and the 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, and the time baselines include the shortest time baseline, the first time baseline, and the second time baseline; A first matrix determination module, configured to determine a first matrix to be quasi-stabilized based on the initial phase unwrapping result; wherein, the first matrix to be quasi-stabilized includes the unwrapped phase corresponding to the shortest time baseline and the unwrapped phase corresponding to the first time baseline; A first matrix quasi-stability module, configured to, for the first matrix to be quasi-stabilized, mark the unwrapped phase corresponding to the shortest time baseline as the first quasi-stable phase, and mark the unwrapped phase corresponding to the first time baseline as the first non-quasi-stable phase; and determine a first gross error estimation value based on the gross error observation equation, the first quasi-stable phase, and the first non-quasi-stable phase; A second matrix determination module, configured to adjust the unwrapped phase corresponding to the first time baseline in the initial phase unwrapping result based on the first gross error estimation value, and determine a second matrix to be quasi-stabilized based on the adjusted initial phase unwrapping result; wherein, the second matrix to be quasi-stabilized includes the unwrapped phase corresponding to the shortest time baseline, the adjusted unwrapped phase corresponding to the first time baseline, and the unwrapped phase corresponding to the second time baseline; A second matrix quasi-stability module, configured to, for the second matrix to be quasi-stabilized, mark the unwrapped phase corresponding to the shortest time baseline and the adjusted unwrapped phase corresponding to the first time baseline as the second quasi-stable phase, and mark the unwrapped phase corresponding to the second time baseline as the second non-quasi-stable phase; and determine a second gross error estimation value based on the gross error observation equation, the second quasi-stable phase, and the second non-quasi-stable phase; A phase unwrapping module, configured to adjust the unwrapped phase corresponding to the second time baseline in the adjusted initial phase unwrapping result based on the second gross error estimation value, so as to obtain a target phase unwrapping result.

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

9. A computer device, comprising a storage medium, a processor, and a computer program stored on the storage medium and executable on the processor, wherein, When the processor executes the computer program, it implements the phase unwrapping method based on improved minimum cost flow and quasi-stability identification according to any one of claims 1 to 6.

Images