Minimum-Cost Flow Phase Unwrapping Method and Device Based on Phase Gradient Rate Aiding

Abstract

The present disclosure provides a minimum-cost flow phase unwrapping method and apparatus based on phase gradient rate assistance, including: obtaining an SAR data set and DEM data of an area to be monitored; and generating a differential interferogram and a coherence map based thereon; calculating a shortest-time baseline phase gradient rate and a full-time baseline phase gradient rate respectively based on the differential interferogram; constructing a Delaunay triangulation network based on the differential interferogram and the coherence map, determining an irrotational constraint and a design matrix; weighting each triangular arc segment in the triangulation network based on the coherence map, the shortest-time baseline phase gradient rate and the full-time baseline phase gradient rate to obtain a weight value corresponding to each triangular arc segment; solving the integer ambiguity based on the irrotational constraint, the design matrix and the weight value corresponding to each triangular arc segment, and completing the phase unwrapping task based on the integer ambiguity. In this embodiment, the absolute phase of a large-gradient deformation area is obtained, thereby realizing high-precision surface deformation monitoring.

Description

Technical Field

[0001] The present disclosure relates to the technical field of surface large-gradient deformation monitoring. Specifically, it relates to a minimum-cost flow phase unwrapping method and device assisted by phase gradient rate. Background Art

[0002] Remote sensing technology based on SAR (Synthetic Aperture Radar) data has shown great potential in the investigation, monitoring, and early warning of geological disasters such as landslides, mining area collapses, volcanoes, and earthquakes. Currently, SAR technology is mainly divided into two categories: phase-based SAR interferometry (InSAR) and intensity-based offset tracking. The InSAR technology can achieve precise observation of surface micro-deformations by acquiring two images of the same area and performing phase processing. The intensity-based offset tracking technology, on the other hand, estimates the offset using the intensity information between images and is suitable for large-gradient deformation monitoring, but the accuracy is relatively low.

[0003] In the InSAR technology, phase unwrapping is a key step that directly affects the accuracy of deformation observation. The minimum-cost flow (MCF) phase unwrapping algorithm is widely used due to its high efficiency and relative accuracy. This algorithm transforms the phase unwrapping problem into a network flow problem of calculating the minimum cost and realizes phase recovery through optimization. However, the MCF algorithm depends on the Itoh condition, that is, the phase gradient difference between adjacent pixels needs to satisfy a certain range. In practical applications, due to the limitation of the SAR satellite sampling resolution, large-gradient surface deformations often cause the Itoh condition to fail, making the MCF algorithm unable to effectively recover the absolute phase of the large-gradient deformation area, thus limiting the application effect of the InSAR technology in complex deformation scenarios. Summary of the Invention

[0004] The embodiments of the present disclosure at least provide a minimum-cost flow phase unwrapping method and device assisted by phase gradient rate, which can achieve high-precision surface deformation monitoring by obtaining the absolute phase of the large-gradient deformation area.

[0005] The embodiments of the present disclosure provide a minimum-cost flow phase unwrapping method assisted by phase gradient rate, including:

[0006] Obtain the SAR dataset and DEM data of the area to be monitored; and determine the initial SAR image set based on the SAR dataset; wherein, the initial SAR image set includes multiple SAR images of the area to be monitored; determine the SAR main image in the initial SAR image set, and based on the SAR main image, register the other SAR images in the initial SAR image set except the SAR main 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 the baseline set based on the target SAR image set according to the preset spatio-temporal baseline threshold; wherein, the baseline set includes multiple interference pairs, and each interference pair includes two target SAR images; and 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;

[0007] For the differential interferogram corresponding to each interference pair, calculate the phase gradient value corresponding to each pixel in the differential interferogram based on the phase gradient calculation formula, and determine the first phase gradient map based on the phase gradient value corresponding to each pixel; wherein, the first phase gradient map includes the first horizontal direction phase gradient map, the first vertical direction phase gradient map, the first diagonal direction phase gradient map, and the second diagonal direction phase gradient map; for each interference pair, determine the gradient wrapping result according to the first phase gradient map corresponding to the differential interferogram based on the wrapping formula, and mask the gradient wrapping result based on the coherence map corresponding to the differential interferogram to obtain the second phase gradient map; wherein, 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;

[0008] 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; and, calculate the first-direction shortest-time baseline stack map and the first-direction full-time baseline stack map respectively according to the first-direction phase gradient map set based on the stacking formula; calculate the second-direction shortest-time baseline stack map and the second-direction full-time baseline stack map respectively according to the second-direction phase gradient map set based on the stacking formula; calculate the third-direction shortest-time baseline stack map and the third-direction full-time baseline stack map respectively according to the third-direction phase gradient map set based on the stacking formula; calculate the fourth-direction shortest-time baseline stack map and the fourth-direction full-time baseline stack map respectively according to the fourth-direction phase gradient map set based on the stacking formula;

[0009] Calculate the shortest-time baseline phase gradient rate based on the gradient fusion formula according to the first-direction shortest-time baseline stack map, the second-direction shortest-time baseline stack map, the third-direction shortest-time baseline stack map, and the fourth-direction shortest-time baseline stack map; calculate the full-time baseline phase gradient rate based on the gradient fusion formula according to the first-direction full-time baseline stack map, the second-direction full-time baseline stack map, the third-direction full-time baseline stack map, and the fourth-direction full-time baseline stack map;

[0010] Construct a Delaunay triangulation based on the differential interferogram and the coherence map, and determine the irrotational constraint and the design matrix based on the Delaunay triangulation; weight each triangular arc segment in the Delaunay triangulation based on the coherence map, the shortest-time baseline phase gradient rate, and the full-time baseline phase gradient rate to obtain the weight values corresponding to each triangular arc segment in the Delaunay triangulation;

[0011] Solve the integer ambiguity based on the irrotational constraint, the design matrix, and the weight values corresponding to each triangular arc segment in the Delaunay triangulation, and complete the phase unwrapping task of the area to be monitored based on the integer ambiguity.

[0012] An embodiment of the present disclosure provides a minimum-cost flow phase unwrapping device assisted by a phase gradient rate, including:

[0013] A data processing module, configured to obtain an SAR data set and DEM data of an area to be monitored; and 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 main image in the initial SAR image set, and based on the SAR main image, register the other SAR images in the initial SAR image set except the SAR main image according to the DEM data, and determine a target SAR image set based on the registration result; wherein the target SAR image set includes 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 interference pairs, and each interference pair includes two target SAR images; and 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;

[0014] A gradient map determination module, configured to, for the differential interferogram corresponding to each interference pair, calculate the phase gradient value corresponding to each pixel in the differential interferogram based on a phase gradient calculation formula, 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; for each interference pair, determine a gradient wrapping result according to the first phase gradient map corresponding to the differential interferogram based on a wrapping formula, 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;

[0015] An atlas determination module, configured to determine a first-direction phase gradient atlas based on the second horizontal-direction phase gradient atlases corresponding to the interference pairs in the baseline set, determine a second-direction phase gradient atlas based on the second vertical-direction phase gradient atlases corresponding to the interference pairs in the baseline set, determine a third-direction phase gradient atlas based on the third diagonal-direction phase gradient atlases corresponding to the interference pairs in the baseline set, and determine a fourth-direction phase gradient atlas based on the fourth diagonal-direction phase gradient atlases corresponding to the interference pairs in the baseline set; and, calculate a first-direction shortest-time baseline stack map and a first-direction full-time baseline stack map respectively according to the first-direction phase gradient atlas based on a stacking formula; calculate a second-direction shortest-time baseline stack map and a second-direction full-time baseline stack map respectively according to the second-direction phase gradient atlas based on the stacking formula; calculate a third-direction shortest-time baseline stack map and a third-direction full-time baseline stack map respectively according to the third-direction phase gradient atlas based on the stacking formula; calculate a fourth-direction shortest-time baseline stack map and a fourth-direction full-time baseline stack map respectively according to the fourth-direction phase gradient atlas based on the stacking formula;

[0016] A rate calculation module, configured to calculate a shortest-time baseline phase gradient rate based on a gradient fusion formula according to the first-direction shortest-time baseline stack map, the second-direction shortest-time baseline stack map, the third-direction shortest-time baseline stack map, and the fourth-direction shortest-time baseline stack map; calculate a full-time baseline phase gradient rate based on the gradient fusion formula according to the first-direction full-time baseline stack map, the second-direction full-time baseline stack map, the third-direction full-time baseline stack map, and the fourth-direction full-time baseline stack map;

[0017] A weight determination module, configured to construct a Delaunay triangulation network based on the differential interferogram and the coherence map, and determine an irrotational constraint and a design matrix based on the Delaunay triangulation network; weight each triangular arc segment in the Delaunay triangulation network based on the coherence map, the shortest-time baseline phase gradient rate, and the full-time baseline phase gradient rate to obtain weights corresponding to the triangular arc segments in the Delaunay triangulation network;

[0018] A phase unwrapping module, configured to solve for the integer ambiguity based on the irrotational constraint, the design matrix, and the weights corresponding to the triangular arc segments in the Delaunay triangulation network, and complete the phase unwrapping task for the area to be monitored based on the integer ambiguity.

[0019] 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 minimum-cost flow phase unwrapping method based on phase gradient rate assistance described in any of the above possible embodiments is executed.

[0020] 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 minimum-cost flow phase unwrapping method based on phase gradient rate assistance described in any of the above possible embodiments is implemented.

[0021] In the embodiment of the present disclosure, the minimum-cost flow phase unwrapping method and device based on phase gradient rate assistance utilize the phase gradient information in the SAR image interferogram that is not affected by phase unwrapping errors to determine the boundary of large-gradient deformation on the ground, optimize the Delaunay triangulation network, and control the unwrapping path to enter the large-gradient deformation field from a position with a smaller deformation gradient for unwrapping, so as to correctly recover the absolute phase of the large-gradient deformation field on the ground, and further achieve high-precision deformation monitoring of the ground.

[0022] To make the above objects, features, and advantages of the present disclosure more obvious and understandable, the following specific preferred embodiments are given, and in conjunction with the accompanying drawings, the detailed description is as follows. Description of the Drawings

[0023] To more clearly illustrate the technical solutions of the embodiments of the present disclosure, the drawings required to be cited in the embodiments will be briefly introduced below. The drawings here are incorporated into the specification and constitute a part of this specification. These drawings show the embodiments that conform to the present disclosure and are used together with the specification to illustrate the technical solutions of the present disclosure. It should be understood that the following drawings only show some embodiments of the present disclosure, and therefore should not be regarded as limiting the scope. For those of ordinary skill in the art, other related drawings can be obtained based on these drawings without creative efforts.

[0024] Figure 1 Shows a flowchart of a minimum-cost flow phase unwrapping method based on phase gradient rate assistance provided by an embodiment of the present disclosure;

[0025] Figure 2 Shows a spatio-temporal baseline distribution diagram of a small baseline set of data used in a simulation experiment and an example provided by an embodiment of the present disclosure;

[0026] Figure 3 Shows a schematic diagram of a ground deformation field of a simulation experiment provided by an embodiment of the present disclosure;

[0027] Figure 4 Shows the original interference pattern and the wrapped interference pattern generated by a simulation experiment provided by an embodiment of the present disclosure;

[0028] Figure 5 Shows a comparison diagram of the MCF phase unwrapping results of a simulation experiment provided by an embodiment of the present disclosure;

[0029] Figure 6 Shows a comparison diagram of the phase unwrapping results of the method provided by the present disclosure in a simulation experiment provided by an embodiment of the present disclosure;

[0030] Figure 7 Shows a schematic structural diagram of a minimum cost flow phase unwrapping device based on phase gradient rate assistance provided by an embodiment of the present disclosure;

[0031] Figure 8 Shows a schematic structural diagram of a computer device provided by an embodiment of the present disclosure. Detailed implementation manners

[0032] 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. Usually, the components of the embodiments of the present disclosure described and illustrated herein can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present disclosure provided in the accompanying drawings is not intended to limit the scope of the present disclosure to be protected, but merely represents selected embodiments of the present disclosure. Based on the embodiments of the present disclosure, all other embodiments obtained by those skilled in the art without creative efforts fall within the scope of protection of the present disclosure.

[0033] 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.

[0034] The term "and / or" in this article merely describes an association relationship and indicates that three relationships may exist. For example, A and / or B may represent: A exists alone, both 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.

[0035] To facilitate the understanding of this embodiment, the execution subject of the minimum cost flow phase unwrapping method based on phase gradient rate assistance provided by the embodiments of the present disclosure will be introduced in detail first. The execution subject of the minimum cost flow phase unwrapping method based on phase gradient rate assistance provided by the embodiments of the present disclosure is a computer device. This computer device can be a terminal device or a server. Among them, the terminal device can also be a mobile device, a user terminal, a terminal, a handheld device, a computing device, etc. The server can be an independent physical server, a server cluster or a distributed system composed of multiple physical servers, or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud storage, big data and artificial intelligence platforms. Optionally, this method can also be applied to an implementation environment composed of a computer device and a server.

[0036] The following will detail the minimum cost flow phase unwrapping method based on phase gradient rate assistance provided by the embodiments of the present application with reference to the accompanying drawings. Refer to Figure 1 As shown, this method includes the following S101 to S106:

[0037] S101, obtain the SAR data set and DEM data of the area to be monitored; and determine the initial SAR image set based on the SAR data set; determine the SAR master image in the initial SAR image set, and 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; determine the baseline set based on the target SAR image set according to the preset spatio-temporal baseline threshold; and 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.

[0038] It is understandable that SAR (Synthetic Aperture Radar) is a technology for imaging the ground through electromagnetic waves, capable of acquiring surface information all day and all weather, and is widely used in fields such as surface change monitoring and disaster warning. The SAR dataset is image data containing electromagnetic wave reflection information of the area to be monitored collected by a radar sensor, which can provide precise surface topography and deformation information. DEM (Digital Elevation Model) is a digital model describing the ground elevation distribution obtained through remote sensing technology or ground measurement. DEM data is used to assist in understanding the terrain undulation and is combined with SAR data through a differential interferogram to derive a more accurate phase unwrapping result. Among them, the initial SAR image set includes multiple SAR images of the area to be monitored; 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.

[0039] In SAR interferometric 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.

[0040] Specifically, the SAR master image refers to the reference image selected from the entire image set. Usually, the image that is closest in time or has the best quality is selected as the benchmark. Registration refers to 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, DEM data is used to correct the geometric errors of the images to ensure 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.

[0041] It is understandable that the spatio-temporal baseline threshold refers to the limitation on the validity between image pairs in terms of time and space, which is usually a parameter value used to screen out representative and relatively consistent interferometric pairs. In the present disclosure, the spatio-temporal baseline threshold is set such that the time is less than 36 days and the spatial distance is less than 150 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 a 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.

[0042] Specifically, an interferometric pair is formed by two radar observations (or images). The time and location of each observation may vary slightly. 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 surface changes. 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.

[0043] It is understandable 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 better signal quality, and then determine which areas have higher 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 caused by atmospheric factors.

[0044] Among them, determining the coherence map corresponding to each interferometric pair based on the differential interferogram corresponding to each interferometric pair can be achieved through a coherence calculation algorithm. This algorithm usually 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 can 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. The value of each pixel in the map represents the degree of signal coherence of that position in the differential interferogram.

[0045] S102. For each differential interferogram corresponding to an interference pair, calculate the phase gradient value corresponding to each pixel in the differential interferogram based on the phase gradient calculation formula, and determine the first phase gradient map based on the phase gradient values corresponding to each pixel; for each interference pair, determine the gradient wrapping result based on the first phase gradient map corresponding to the differential interferogram according to the wrapping formula, and mask the gradient wrapping result based on the coherence map corresponding to the differential interferogram to obtain the second phase gradient map.

[0046] Specifically, for each differential interferogram corresponding to an interference pair, based on the phase gradient calculation formula, 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 the first horizontal direction phase gradient map, the first vertical direction phase gradient map, the first diagonal direction phase gradient map, and the second diagonal direction phase gradient map.

[0047] Here, the phase gradient calculation formula can be expressed as:

[0048] ;

[0049] Among them, represents the phase gradient in a specific direction of the pixel; 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.

[0050] It can be understood that for each interference pair, first determine the gradient wrapping result based on the first phase gradient map according to the wrapping formula. The phase wrapping phenomenon occurs when the phase gradient value exceeds a certain range. Therefore, it is necessary to use the wrapping formula to correct the phase gradient to avoid calculation errors caused by wrapping. At this time, there will be random phase gradients caused by random noise and local micro gradients caused by atmospheric errors in the phase gradient. To solve this phenomenon, the gradient wrapping result can be masked in combination with the coherence map to obtain the second phase gradient map. In this way, by retaining effective gradient information in areas with high coherence and removing invalid or incorrect gradient information in areas with low 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;

[0051] Here, the wrapping formula can be expressed as:

[0052] ;

[0053] Among them, is expressed as the gradient winding result; is expressed as the phase gradient.

[0054] S103, determine the first-direction phase gradient atlas based on the second horizontal-direction phase gradient maps corresponding to each interference pair in the baseline set, determine the second-direction phase gradient atlas based on the second vertical-direction phase gradient maps corresponding to each interference pair in the baseline set, determine the third-direction phase gradient atlas based on the third diagonal-direction phase gradient maps corresponding to each interference pair in the baseline set, determine the fourth-direction phase gradient atlas based on the fourth diagonal-direction phase gradient maps corresponding to each interference pair in the baseline set; and, calculate the first-direction shortest-time baseline stack map and the first-direction full-time baseline stack map respectively according to the first-direction phase gradient atlas based on the stacking formula; calculate the second-direction shortest-time baseline stack map and the second-direction full-time baseline stack map respectively according to the second-direction phase gradient atlas based on the stacking formula; calculate the third-direction shortest-time baseline stack map and the third-direction full-time baseline stack map respectively according to the third-direction phase gradient atlas based on the stacking formula; calculate the fourth-direction shortest-time baseline stack map and the fourth-direction full-time baseline stack map respectively according to the fourth-direction phase gradient atlas based on the stacking formula.

[0055] Specifically, for different directional phase gradient atlases, respectively determine the phase gradient atlases in different directions based on the second phase gradient maps corresponding to each interference pair. For example, based on the phase gradient map in the horizontal direction, obtain the first-direction phase gradient atlas of all interference pairs; based on the phase gradient map in the vertical direction, obtain the second-direction phase gradient atlas, and so on. In this way, it is possible to accumulate data of multiple interference pairs for each direction, so as to better capture the phase change characteristics in different directions.

[0056] It can be understood that through the stacking formula, stack the phase gradient atlases in each direction, and calculate the shortest-time baseline stack map and the full-time baseline stack map corresponding to each direction. Stacking can effectively fuse the phase gradient information from different interference pairs, so as to synthesize the phase gradient characteristics under different time baselines in time. The shortest-time baseline stack map focuses on the phase changes in a short time, while the full-time baseline stack map synthesizes the phase changes in a long time.

[0057] Here, the stacking formula can be expressed as:

[0058] ;

[0059] Among them, It is expressed as the stacked result of the phase gradient in the k direction, where the value of k is 0, 45, 90, 135; m represents the number of interference pairs participating in the calculation of the stacked phase gradient; It is expressed as the time baseline of the differential interferogram corresponding to the nth interference pair.

[0060] S104, according to the shortest time baseline stacking map in the first direction, the shortest time baseline stacking map in the second direction, the shortest time baseline stacking map in the third direction, and the shortest time baseline stacking map in the fourth direction, calculate the shortest time baseline phase gradient rate based on the gradient fusion formula; according to the full time baseline stacking map in the first direction, the full time baseline stacking map in the second direction, the full time baseline stacking map in the third direction, and the full time baseline stacking map in the fourth direction, calculate the full time baseline phase gradient rate based on the gradient fusion formula.

[0061] It can be understood that the phase gradient rate refers to the change rate of the phase difference between adjacent pixels and is used to describe the spatial change characteristics of the phase. In the present disclosure, for different time baselines, the shortest time baseline phase gradient rate and the full time baseline phase gradient rate can be calculated respectively.

[0062] Specifically, through the gradient fusion formula, combined with the shortest time baseline stacking maps in the first to fourth directions, the shortest time baseline phase gradient rate is calculated. In this way, by fusing the phase gradient information in different directions, a comprehensive phase gradient rate can be obtained to reflect the rate of phase change in a short time.

[0063] Here, the gradient fusion formula can be expressed as:

[0064] ;

[0065] Among them, G represents the fusion result of the phase gradient stacking maps in four directions; the phase gradient rate obtained by fusing the full time baseline stacking maps in the first direction, the second direction, the third direction, and the fourth direction is called the full time baseline phase gradient rate, denoted as The phase gradient rate obtained by fusing the shortest time baseline stacking maps in the first direction, the second direction, the third direction, and the fourth direction is called the shortest time baseline phase gradient rate, denoted as .

[0066] Similarly, based on the full time baseline stacking maps in the first to fourth directions, the full time baseline phase gradient rate is calculated using the above gradient fusion formula to reflect the rate of phase change on a long-term time scale.

[0067] S105. Construct a Delaunay triangulation based on the differential interferogram and coherence map, and determine the irrotational constraint and design matrix based on the Delaunay triangulation; weight each triangular arc segment in the Delaunay triangulation based on the coherence map, the shortest-time baseline phase gradient rate, and the full-time baseline phase gradient rate to obtain the weight values corresponding to each triangular arc segment in the Delaunay triangulation.

[0068] It can be understood that the Delaunay triangulation is a triangular grid structure generated by a set of scattered points, with good mathematical properties, and is widely used in the fields of geographic information system (GIS) and computational geometry. This triangulation requires that no other points are contained within the circumcircle of any triangle, thus avoiding the appearance of "long and narrow" triangles. In this solution, the Delaunay triangulation is used to divide the area to be monitored into multiple triangular regions and further used to define the constraint conditions for phase unwrapping.

[0069] Specifically, when constructing the Delaunay triangulation based on the differential interferogram and coherence map, the following (1) to (2) may be included:

[0070] (1) For each differential interferogram corresponding to an interferometric pair, extract high-quality points that meet the preset coherence threshold based on the coherence map corresponding to the differential interferogram, and determine the high-quality point set corresponding to the interferometric pair based on the extraction result.

[0071] (2) For each high-quality point set corresponding to an interferometric pair, use the point-by-point insertion method to generate the Delaunay triangulation corresponding to the interferometric pair; wherein, each Delaunay triangulation includes multiple triangles.

[0072] It can be understood that high-quality points that meet the preset coherence threshold are extracted through the coherence map, and these points represent the ground object information with high coherence in the monitoring area. By selecting high-quality points, it can be ensured that in the subsequent triangulation construction process, the point set used has high accuracy and representativeness, ensuring the quality of the points in the Delaunay triangulation. Then, based on these extracted high-quality points, the Delaunay triangulation can be generated using the point-by-point insertion method, where each Delaunay triangulation includes multiple triangles. The point-by-point insertion method is a classic mesh generation method that adjusts the existing triangular structure by adding new points one by one, and finally forms a grid that meets the Delaunay conditions.

[0073] It is understandable that the irrotationality constraint means that during the unwrapping process, the phase change between adjacent regions must meet certain physical consistency, that is, the unwrapping result cannot have rotation or unreasonable phase changes, ensuring the continuity and physical rationality of the phase during the unwrapping process. The design matrix is a mathematical tool used to transform various constraints and known conditions into computable mathematical forms. In the phase unwrapping problem, the design matrix is used to describe the mutual relationships and their constraint conditions between each triangle, ultimately helping to solve the overall phase.

[0074] Specifically, for the Delaunay triangulation corresponding to each interference pair, construct a design matrix corresponding to the Delaunay triangulation; and, calculate the residuals of each triangle in the Delaunay triangulation based on the residual calculation formula according to the second phase gradient map, and construct the irrotationality constraint based on the residuals of each triangle in the Delaunay triangulation.

[0075] Here, the design matrix can be expressed as:

[0076] ,

[0077] ;

[0078] Among them, ; M represents the number of triangles in the Delaunay triangulation; N represents N arc segments; for the triangular closed loop , the side with a counterclockwise rotation direction The corresponding matrix element , the side with a clockwise rotation direction in the closed loop The corresponding matrix element , if is not in , then the matrix element .

[0079] Here, the residual calculation formula can be expressed as:

[0080] ;

[0081] Among them, represents the residual of the th triangle; respectively represent the three vertices of the th triangle; represents the wrapped phase gradient between vertex and vertex ; represents the wrapped phase gradient between vertex and vertex ; represents the vertex and the vertex between the winding phase gradients.

[0082] Exemplarily, in the Delaunay triangulation, each arc segment of a triangle corresponds to a weight value. The weight value is determined based on the phase gradient rate in the differential interferogram and the information in the coherence map. The phase gradient rate and the coherence map provide the signal strength and change speed of the phase change between different regions, thereby affecting the connection weights between adjacent regions. During the unwrapping process, reasonable weight allocation helps to balance the contributions of different regions to the final unwrapping result and improve the unwrapping accuracy.

[0083] Specifically, after obtaining the coherence map, the phase gradient rate of the shortest time baseline, and the phase gradient rate of the full time baseline, for the Delaunay triangulation corresponding to each interferometric pair, based on the arc segment weight calculation formula, the weights of the arc segments of each triangle in the Delaunay triangulation corresponding to each interferometric pair are determined according to the coherence map, the phase gradient rate of the shortest time baseline, and the phase gradient rate of the full time baseline.

[0084] Here, the arc segment weight calculation formula includes:

[0085] ;

[0086] wherein, 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 shortest time baseline at the starting position of the arc segment ; represents the phase gradient rate of the shortest time baseline at the ending position of the arc segment ; represents the phase gradient rate threshold, and its value is , represents the phase gradient rate of the full time baseline; represents the coherence value at the starting position of the arc segment ; represents the coherence value at the ending position of the arc segment .

[0087] S106. Solve the integer ambiguity based on the irrotational constraint, the design matrix, and the weights corresponding to the arc segments of each triangle in the Delaunay triangulation, and complete the phase unwrapping task of the area to be monitored based on the integer ambiguity.

[0088] It is understandable that ambiguity resolution (AR) is a common technique in systems such as satellite positioning (e.g., GNSS) and radar (e.g., SAR). Especially in the process of unwrapping phase measurements, the ambiguity refers to the uncertainty caused by periodicity in the phase data, and how to accurately determine the ambiguity (i.e., integer multiples of the wavelength). In SAR, due to the periodicity of the phase values, multiple ambiguity resolution problems may occur during the unwrapping process. By solving the ambiguity based on the curl-free constraint, the design matrix, and the calculated weights, inconsistent phase ambiguities can be eliminated, making the phase unwrapping more accurate. Finally, by determining the ambiguity, the phase of the area to be monitored can be effectively unwrapped, completing the phase unwrapping task, thereby providing accurate terrain deformation information or relevant data of other monitoring targets.

[0089] Specifically, when solving the ambiguity, the following (a) to (b) may be included:

[0090] (a) Construct an objective function;

[0091] (b) For each Delaunay triangulation corresponding to an interference pair, solve the ambiguity based on the objective function according to the curl-free constraint, the design matrix, and the weights corresponding to each triangular arc segment in the Delaunay triangulation corresponding to the interference pair.

[0092] Here, the objective function can be expressed as:

[0093] ;

[0094] where W represents a vector including the weights corresponding to each triangular arc segment in the Delaunay triangulation; X represents the ambiguity; A represents the design matrix; and U represents the curl-free constraint.

[0095] It is understandable that after obtaining the ambiguity corresponding to each interference pair, the ambiguity corresponding to each interference pair can be integrated to obtain the unwrapping result corresponding to each interference pair; and the phase unwrapping task of the area to be monitored can be completed based on the unwrapping result corresponding to each interference pair.

[0096] In the method and device for minimum cost flow phase unwrapping assisted by phase gradient rate provided in the embodiments of the present disclosure, by using the phase gradient information in the SAR image interferogram that is not affected by phase unwrapping errors, the boundary of large gradient deformations on the ground surface is judged, the Delaunay triangulation is optimized, and the unwrapping path is controlled to enter the large gradient deformation field from a position with a smaller deformation gradient for unwrapping, so as to correctly restore the absolute phase of the large gradient deformation field on the ground surface, and further achieve high-precision deformation monitoring of the ground surface.

[0097] To verify the performance of the proposed method for phase unwrapping, the present disclosure obtained a total of 91 scenes of Sentinel-1 descending orbit data from January 2, 2019 to December 29, 2021 for simulation experiments and real experiments. This satellite is a C-band SAR satellite developed by the European Space Agency, which can achieve full coverage and periodic observation of the global land surface. And due to the high temporal resolution (12 / 6 days) and spatial resolution (~20 meters) of Sentinel-1 SAR satellite data and the free and open access policy, it has become a commonly used data source in the InSAR field. Then, 267 interferometric pairs with a maximum temporal baseline of 36 days and a maximum spatial baseline of 150 meters were generated, and the spatio-temporal baseline map is as shown in Figure 2 shown. Then, simulation experiments and real data experiments were carried out based on the obtained Sentinel-1 SAR data. Among them, the deformation field in the simulation experiment was simulated, and the spatio-temporal baseline was real data. In the real experiment, both the surface deformation field and the spatio-temporal baseline were real data.

[0098] Specifically, as shown in Figure 3 shown, the present disclosure simulated three common types of surface deformation fields: circular deformation field, trapezoidal deformation field, and rectangular deformation field. Among them, the deformation of the circular deformation field slowly increases from the boundary of the circle to the center of the circle, and there is no obvious deformation boundary compared with the stable area; the deformation of the trapezoidal deformation field slowly increases from top to bottom, and there are obvious deformation boundaries on both sides and the lower side compared with the stable area, and there is no obvious deformation boundary on the upper side; there is an obvious deformation boundary between the rectangular deformation field and the stable area. Then, 264 original interferograms were generated according to the spatio-temporal baseline, and the phase of the original interferogram was wrapped to obtain the wrapped interferogram. When generating the simulated interferogram, on the basis of linear deformation, a periodic deformation with an annual cycle was added. The first three interferometric phase diagrams generated by simulation are as shown in Figure 4 shown, where (a)-(c) are the simulated original interferometric phase diagrams, and their temporal baselines are 12 days, 24 days, and 36 days respectively. (d)-(f) are the interferograms obtained by wrapping (a)-(c). As can be seen from Figure 2 it, the phase of the wrapped interferogram is between, and when the temporal baseline is 24 days and 36 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.

[0099] After generating the simulated wrapped phase diagram, phase unwrapping was carried out respectively by the traditional MCF method and the method provided by the embodiment of the present disclosure. The unwrapping result of the traditional MCF method is as shown in Figure 5 shown. Figure 5 (a) is the unwrapped phase diagram obtained from the wrapped phase diagram with a temporal baseline of 12 days by the MCF method. Similarly, Figure 5 (b) and Figure 4(c) Unwrapped phase diagrams obtained by the MCF method from wrapped phase diagrams with 24-day and 36-day winding phases respectively. Figure 5 (d)-(f) Differences between the unwrapped phase diagrams and the original phase diagrams with the same time baselines respectively. As can be seen from Figure 5 (a), when the time baseline is 12 days, the phase is correctly restored without unwrapping error, and only the error caused by the computer storage precision exists in Figure 5 (d). When the time baseline is 24 days, the circular deformation area is correctly unwrapped, and the correct unwrapping results are obtained for some areas in the trapezoidal deformation area, while all areas in the rectangular deformation area are incorrectly unwrapped. When the time baseline is extended to 36 days, the circular deformation area still obtains the correct phase unwrapping result, while the area with correct results in the trapezoidal deformation area decreases, and all areas in the rectangular deformation area still obtain incorrect phase unwrapping results.

[0100] The unwrapping results of the method provided by the embodiments of the present disclosure are as shown in Figure 6 the figure Figure 6 (a) is the unwrapped phase diagram obtained from the wrapped phase diagram with a 12-day time baseline by the proposed phase unwrapping method. Similarly, Figure 6 (b) and Figure 6 (c) are the unwrapped phase diagrams obtained by phase unwrapping the wrapped phase diagrams with 24-day and 36-day winding phases respectively by the proposed method. Figure 6 (d)-(f) Differences between the unwrapped phase diagrams and the simulated original phase diagrams with the same time baselines respectively. As can be seen from the figure, when the time baseline is 12 days, the phase is correctly restored without unwrapping error, and only the error caused by the computer storage precision exists in Figure 6 (d). When the time baselines are 24 days and 36 days, the circular deformation area obtains the correct unwrapping result, and the correct unwrapping results are obtained for most areas in the trapezoidal deformation area, and only scattered points with incorrect phase unwrapping exist at the edge of the deformation field. While all areas in the rectangular deformation area still obtain incorrect unwrapping results, because the wrapped phase gradient disappears at the edge of the rectangular deformation field when the time baseline is greater than 12 days, which leads to the unwrapping network being unable to enter the interior of the deformation field from the boundary with a smaller deformation gradient for unwrapping.

[0101] 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 on the implementation process. The specific execution order of each step should be determined according to its function and possible internal logic.

[0102] Based on the same inventive concept, an apparatus for minimum-cost flow phase unwrapping assisted by phase gradient rate corresponding to the minimum-cost flow phase unwrapping method assisted by phase gradient rate is further provided in the embodiments of the present disclosure. Since the principle of problem-solving of the apparatus in the embodiments of the present disclosure is similar to that of the above-mentioned method in the embodiments of the present disclosure, the implementation of the apparatus can refer to the implementation of the method, and the repeated parts will not be elaborated.

[0103] Referring Figure 7 As shown, it is a schematic diagram of an apparatus 700 for minimum-cost flow phase unwrapping assisted by phase gradient rate provided in the embodiments of the present disclosure. The apparatus includes:

[0104] A data processing module 701, configured to obtain an SAR data set and DEM data of a to-be-monitored area; and determine an initial SAR image set based on the SAR data set; wherein, the initial SAR image set includes multiple SAR images of the to-be-monitored area; determine the SAR master image in the initial SAR image set, and based on the SAR master image, register other SAR images in the initial SAR image set except the SAR master image according to the DEM data, and determine a target SAR image set based on the registration result; wherein, the target SAR image set includes 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 interference pairs, and each interference pair includes two target SAR images; and 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;

[0105] A gradient map determination module 702, configured to, for the differential interferogram corresponding to each interference pair, calculate the phase gradient value corresponding to each pixel in the differential interferogram based on a phase gradient calculation formula, 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; for each interference pair, determine a gradient wrapping result based on the first phase gradient map corresponding to the differential interferogram according to a wrapping formula, 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;

[0106] An atlas determination module 703 is configured to determine a first-direction phase gradient atlas based on the second horizontal-direction phase gradient maps corresponding to the interference pairs in the baseline set, determine a second-direction phase gradient atlas based on the second vertical-direction phase gradient maps corresponding to the interference pairs in the baseline set, determine a third-direction phase gradient atlas based on the third diagonal-direction phase gradient maps corresponding to the interference pairs in the baseline set, and determine a fourth-direction phase gradient atlas based on the fourth diagonal-direction phase gradient maps corresponding to the interference pairs in the baseline set; and calculate a first-direction shortest-time baseline stack map and a first-direction full-time baseline stack map respectively according to the first-direction phase gradient atlas based on a stacking formula; calculate a second-direction shortest-time baseline stack map and a second-direction full-time baseline stack map respectively according to the second-direction phase gradient atlas based on the stacking formula; calculate a third-direction shortest-time baseline stack map and a third-direction full-time baseline stack map respectively according to the third-direction phase gradient atlas based on the stacking formula; calculate a fourth-direction shortest-time baseline stack map and a fourth-direction full-time baseline stack map respectively according to the fourth-direction phase gradient atlas based on the stacking formula;

[0107] A rate calculation module 704 is configured to calculate a shortest-time baseline phase gradient rate based on a gradient fusion formula according to the first-direction shortest-time baseline stack map, the second-direction shortest-time baseline stack map, the third-direction shortest-time baseline stack map, and the fourth-direction shortest-time baseline stack map; calculate a full-time baseline phase gradient rate based on the gradient fusion formula according to the first-direction full-time baseline stack map, the second-direction full-time baseline stack map, the third-direction full-time baseline stack map, and the fourth-direction full-time baseline stack map;

[0108] A weight determination module 705 is configured to construct a Delaunay triangulation network based on the differential interferogram and the coherence map, and determine an irrotational constraint and a design matrix based on the Delaunay triangulation network; determine weights for each triangular arc segment in the Delaunay triangulation network based on the coherence map, the shortest-time baseline phase gradient rate, and the full-time baseline phase gradient rate, to obtain weights corresponding to each triangular arc segment in the Delaunay triangulation network;

[0109] A phase unwrapping module 706 is configured to solve for the integer ambiguity based on the irrotational constraint, the design matrix, and the weights corresponding to each triangular arc segment in the Delaunay triangulation network, and complete the phase unwrapping task of the area to be monitored based on the integer ambiguity.

[0110] In some possible embodiments, the phase gradient calculation formula includes:

[0111] ;

[0112] Among them, is expressed as the phase gradient in a specific direction of the pixel points; is expressed as the phase gradient in the vertical direction; is expressed as the phase gradient in the first diagonal direction; is expressed as the phase gradient in the horizontal direction; is expressed as the phase gradient in the second diagonal direction;

[0113] The winding formula includes:

[0114] ;

[0115] Among them, is expressed as the gradient winding result; is expressed as the phase gradient;

[0116] The stacking formula includes:

[0117] ;

[0118] Among them, is expressed as the stacking result of the phase gradient in the k direction, where the value of k is 0, 45, 90, 135; m represents the number of interference pairs participating in the phase gradient stacking calculation; is expressed as the time baseline of the differential interferogram corresponding to the nth interference pair;

[0119] The gradient fusion formula includes:

[0120] ;

[0121] Among them, G is expressed as the fusion result of the phase gradient stacking diagrams in four directions; the phase gradient rate obtained by fusing the full-time baseline stacking diagrams in the first direction, the full-time baseline stacking diagrams in the second direction, the full-time baseline stacking diagrams in the third direction, and the full-time baseline stacking diagrams in the fourth direction is called the full-time baseline phase gradient rate, denoted as , and the phase gradient rate obtained by fusing the shortest-time baseline stacking diagrams in the first direction, the shortest-time baseline stacking diagrams in the second direction, the shortest-time baseline stacking diagrams in the third direction, and the shortest-time baseline stacking diagrams in the fourth direction is called the shortest-time baseline phase gradient rate, denoted as .

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

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

[0124] 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.

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

[0126] For each Delaunay triangulation corresponding to an interference pair, construct a design matrix corresponding to the Delaunay triangulation; and, calculate the residuals of each triangle in the Delaunay triangulation based on the residual calculation formula according to the second phase gradient map, and construct an irrotational constraint based on the residuals of each triangle in the Delaunay triangulation;

[0127] The design matrix includes:

[0128] ,

[0129] ;

[0130] Wherein, ; M represents the number of triangles in the Delaunay triangulation; N represents N arcs; for a triangular closed loop , the side with a counterclockwise rotation direction The corresponding matrix element , the side with a clockwise rotation direction in the closed loop The corresponding matrix element , if Is not in , then the matrix element ;

[0131] The residual calculation formula includes:

[0132] ;

[0133] Wherein, Represents the residual of the rd triangle; Respectively represent the three vertices of the th triangle; Represents the wrapped phase gradient between vertex and vertex ; Represents the wrapped phase gradient between vertex and vertex ; Represents the wrapped phase gradient between vertex and vertex ;

[0134] The weight determination module 705 is specifically configured to: for each Delaunay triangulation corresponding to an interference pair, based on the arc weight calculation formula, determine weights for each triangular arc in the Delaunay triangulation corresponding to each interference pair according to the coherence map, the shortest time baseline phase gradient rate, and the full time baseline phase gradient rate;

[0135] The arc weight calculation formula includes:

[0136] ;

[0137] wherein, represents the weight of arc ; the starting position of arc is and the ending position is ; represents the shortest time baseline phase gradient rate at the starting position of arc ; represents the shortest time baseline phase gradient rate at the ending position of arc ; represents the phase gradient rate threshold, and its value is , represents the full time baseline phase gradient rate; represents the coherence value at the starting position of arc ; represents the coherence value at the ending position of arc .

[0138] In some possible embodiments, the phase unwrapping module 706 is specifically configured to:

[0139] Construct an objective function;

[0140] For each Delaunay triangulation corresponding to an interference pair, solve for the integer ambiguity based on the objective function according to the irrotational constraint, the design matrix, and the weights corresponding to each triangular arc in the Delaunay triangulation corresponding to the interference pair;

[0141] The objective function includes:

[0142] ;

[0143] wherein, W represents a vector including the weights corresponding to each triangular arc in the Delaunay triangulation; X represents the integer ambiguity; A represents the design matrix; U represents the irrotational constraint;

[0144] The phase unwrapping module 706 is specifically configured to: integrate the integer ambiguity corresponding to each interference pair to obtain an unwrapping result corresponding to each interference pair; and complete the phase unwrapping task of the area to be monitored based on the unwrapping result corresponding to each interference pair.

[0145] Based on the same inventive concept, an embodiment of the present disclosure also provides a computer device. Referring to Figure 8 As shown, it is a schematic structural diagram of a computer device 800 provided by an embodiment of the present disclosure, including a processor 801, a memory 802, and a bus 803. Among them, the memory 802 is used to store execution instructions, including an internal memory 8021 and an external memory 8022; here, the internal memory 8021 is also called the main memory, which is used to temporarily store the operation data in the processor 801 and the data exchanged with the external memory 8022 such as a hard disk. The processor 801 exchanges data with the external memory 8022 through the internal memory 8021.

[0146] In an embodiment of the present application, the memory 802 is specifically configured to store the application program code for implementing the solution of the present application and is controlled by the processor 801 to execute. That is, when the computer device 800 runs, the processor 801 communicates with the memory 802 through the bus 803, so that the processor 801 executes the application program code stored in the memory 802, and further executes the method described in any of the foregoing embodiments. Among them, the memory 802 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 801 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.

[0147] 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 800. In other embodiments of the present application, the computer device 800 may include more or fewer components than shown in the figure, or combine certain components, or split certain components, or have different component arrangements. The illustrated components may be implemented in hardware, software, or a combination of software and hardware.

[0148] The embodiments of the present disclosure also provide 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 minimum-cost flow phase unwrapping method based on phase gradient rate assistance described in the above method embodiments. Among them, the storage medium can be a volatile or non-volatile computer-readable storage medium.

[0149] The embodiments of the present disclosure also provide a computer program product, which carries program codes. The instructions included in the program codes can be used to execute the steps of the minimum-cost flow phase unwrapping method based on phase gradient rate assistance described in the above method embodiments. For details, please refer to the above method embodiments and will not be elaborated here. Among them, the above computer program product can be specifically implemented in the form 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. Those skilled in the art can clearly understand that for the convenience and conciseness of description, the specific working processes of the above-described systems and devices can refer to the corresponding processes in the foregoing method embodiments and will not be elaborated here. In several embodiments provided by 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. Also, for 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 displayed or discussed mutual coupling or direct coupling or communication connection can be through some communication interfaces. The indirect coupling or communication connection of the devices or units can be in an electrical, mechanical, or other form.

[0150] The unit described as a separation component may or may not be physically separated. The component shown as a unit may or may not be a physical unit, that is, it may 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, each functional unit can be integrated in a processing unit, or each unit can exist physically alone, or two or more units can be integrated in one unit. If the function is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a non-volatile computer-readable storage medium 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 part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device 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 that can store program codes such as USB flash drives, mobile hard disks, read-only memories, random access memories, magnetic disks, or optical discs.

[0151] Finally, it should be noted that the above embodiments are only specific implementation manners of the present disclosure, used to illustrate the technical solutions of the present disclosure, rather than limiting them. 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 easily conceive of changes, or perform equivalent replacements on some of the technical features; and these modifications, changes or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present disclosure, and should all be covered by 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 minimum cost flow phase unwrapping method based on phase gradient rate assistance, characterized in that: include: Obtain SAR data set and DEM data of the area to be monitored; and determining an initial SAR image set based on the SAR data set; wherein the initial SAR image set includes a plurality of SAR images of the area to be monitored; determining a SAR main image in the initial SAR image set, taking the SAR main image as a reference, registering other SAR images in the initial SAR image set except the SAR main image according to the DEM data, and determining a target SAR image set based on the registration result; wherein the target SAR image set includes a plurality of target SAR images; determining a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; wherein the baseline set includes a plurality of interference pairs, each of which includes two target SAR images; and generating a differential interference map corresponding to each interference pair based on the DEM data; and generating a coherence map corresponding to each interference pair based on the differential interference map corresponding to each interference pair; For the differential interference image corresponding to each interference pair, the phase gradient value corresponding to each pixel in the differential interference image is calculated based on the phase gradient calculation formula, and the first phase gradient image is determined based on the phase gradient value corresponding to each pixel; wherein the first phase gradient image includes a first horizontal direction phase gradient image, a first vertical direction phase gradient image, a first diagonal direction phase gradient image, and a second diagonal direction phase gradient image; for each interference pair, the gradient winding result is determined according to the first phase gradient image corresponding to the differential interference image based on the winding formula, and the gradient winding result is masked based on the coherence image corresponding to the differential interference image to obtain a second phase gradient image; wherein the second phase gradient image includes a second horizontal direction phase gradient image, a second vertical direction phase gradient image, a third diagonal direction phase gradient image, and a fourth diagonal direction phase gradient image; Determine a first direction phase gradient atlas based on the second horizontal direction phase gradient map corresponding to each interference pair in the baseline set, determine a second direction phase gradient atlas based on the second vertical direction phase gradient map corresponding to each interference pair in the baseline set, determine a third direction phase gradient atlas based on the third diagonal direction phase gradient map corresponding to each interference pair in the baseline set, and determine a fourth direction phase gradient atlas based on the fourth diagonal direction phase gradient map corresponding to each interference pair in the baseline set; and, calculate the first direction shortest time baseline stacking map and the first direction full time baseline stacking map respectively based on the stacking formula according to the first direction phase gradient atlas; calculate the second direction shortest time baseline stacking map and the second direction full time baseline stacking map respectively based on the stacking formula according to the second direction phase gradient atlas; calculate the third direction shortest time baseline stacking map and the third direction full time baseline stacking map respectively based on the stacking formula according to the third direction phase gradient atlas; calculate the fourth direction shortest time baseline stacking map and the fourth direction full time baseline stacking map respectively based on the stacking formula according to the fourth direction phase gradient atlas; According to the shortest time baseline stacking diagram in the first direction, the shortest time baseline stacking diagram in the second direction, the shortest time baseline stacking diagram in the third direction, and the shortest time baseline stacking diagram in the fourth direction, the shortest time baseline phase gradient rate is calculated based on the gradient fusion formula; according to the full time baseline stacking diagram in the first direction, the full time baseline stacking diagram in the second direction, the full time baseline stacking diagram in the third direction, and the full time baseline stacking diagram in the fourth direction, the full time baseline phase gradient rate is calculated based on the gradient fusion formula; A Delaunay triangulation network is constructed based on the differential interference graph and the coherence graph, and a curl-free constraint and a design matrix are determined based on the Delaunay triangulation network; each triangular arc segment in the Delaunay triangulation network is weighted based on the coherence graph, the shortest time baseline phase gradient rate, and the full time baseline phase gradient rate to obtain a weight corresponding to each triangular arc segment in the Delaunay triangulation network; Integer ambiguity is solved based on the irrotational constraint, the design matrix and the weights corresponding to each triangular arc segment in the Delaunay triangulation network, and the phase unwrapping task of the monitored area is completed based on the integer ambiguity.

2. The method according to claim 1, characterized in that: The phase gradient calculation formula includes: ； in, It is expressed as the phase gradient in a specific direction of the pixel; It is expressed as the phase gradient in the vertical direction; It is represented as the phase gradient in the first diagonal direction; It is expressed as the horizontal phase gradient; It is expressed as the phase gradient in the second diagonal direction; The winding formula includes: ； in, Represented as gradient winding result; It is expressed as phase gradient; The stacking formula includes: ； in, It is represented by the phase gradient stacking result in the k direction, and the values ​​of k are 0, 45, 90, and 135; m is represented by the number of interference pairs participating in the phase gradient stacking calculation; It is represented as the time basis of the differential interferogram corresponding to the nth interferometer pair; The gradient fusion formula includes: ； Where G represents the fusion result of the four-directional phase gradient stacking images; the phase gradient rate obtained by fusion of the first direction full-time baseline stacking image, the second direction full-time baseline stacking image, the third direction full-time baseline stacking image, and the fourth direction full-time baseline stacking image is called the full-time baseline phase gradient rate, which is recorded as The phase gradient rate obtained by fusing the shortest time baseline stacking diagram in the first direction, the shortest time baseline stacking diagram in the second direction, the shortest time baseline stacking diagram in the third direction, and the shortest time baseline stacking diagram in the fourth direction is called the shortest time baseline phase gradient rate, which is recorded as .

3. The method according to claim 2, characterized in that The constructing of a Delaunay triangulation network based on the differential interference graph and the coherence graph comprises: For the differential interference graph corresponding to each interference pair, high-quality points satisfying a preset coherence threshold are extracted based on a coherence graph corresponding to the differential interference graph, and a high-quality point set corresponding to the interference pair is determined based on the extraction result; For each high-quality point set corresponding to the interference pair, a point-by-point interpolation method is used to generate a Delaunay triangulation corresponding to the interference pair; wherein each Delaunay triangulation includes a plurality of triangles.

4. The method according to claim 3, characterized in that The step of determining the irrotational constraint and the design matrix based on the Delaunay triangulation network includes: For the Delaunay triangulation corresponding to each interference pair, a design matrix corresponding to the Delaunay triangulation is constructed; and based on the residual calculation formula, the residual of each triangle in the Delaunay triangulation is calculated according to the second phase gradient map, and an irrotational constraint is constructed based on the residual of each triangle in the Delaunay triangulation; The design matrix includes: ， ； in, ; M represents the number of triangles in the Delaunay triangulation; N represents N arc segments; for a triangular closed loop , the counterclockwise side The corresponding matrix elements , the clockwise edge of the closed loop The corresponding matrix elements ,like Not Available If ; The residual calculation formula includes: ； in, Expressed as The residual of the triangle; Respectively expressed as The three vertices of a triangle; Represented as a vertex and vertices The winding phase gradient between Represented as a vertex and vertices The winding phase gradient between Represented as a vertex and vertices The winding phase gradient between The step of weighting each triangular arc segment in the Delaunay triangulation network based on the shortest time baseline phase gradient rate and the full time baseline phase gradient rate includes: For the Delaunay triangulation corresponding to each interference pair, based on the arc segment weight calculation formula, weight each triangular arc segment in the Delaunay triangulation corresponding to each interference pair according to the coherence graph, the shortest time baseline phase gradient rate and the full time baseline phase gradient rate; The arc weight calculation formula includes: ； in, Represented as an arc segment The weight of The starting point is The end point is ; Represented as an arc segment The shortest time baseline phase gradient rate at the starting position; Represented as an arc segment The shortest time baseline phase gradient rate at the endpoint position; It is represented as the phase gradient rate threshold, and its value is , It is expressed as the full-time baseline phase gradient rate; Represented as an arc segment The coherence value of the starting position of ; Represented as an arc segment The coherence value of the end position.

5. The method according to claim 4, characterized in that The step of solving integer ambiguity based on the irrotational constraint, the design matrix and the weights corresponding to each triangular arc segment in the Delaunay triangulation network comprises: Construct the objective function; For the Delaunay triangulation corresponding to each interference pair, solving integer ambiguity according to the irrotational constraint corresponding to the interference pair, the design matrix and the weights corresponding to each triangular arc segment in the Delaunay triangulation based on the objective function; The objective function includes: ； Wherein, W represents a vector including weights corresponding to each triangular arc segment in the Delaunay triangulation network; X represents the integer ambiguity; A represents the design matrix; U represents the irrotational constraint; The step of completing the phase unwrapping task of the area to be monitored based on the integer ambiguity includes: The integer ambiguity corresponding to each interference pair is integrated to obtain an unwrapping result corresponding to each interference pair; and the phase unwrapping task of the area to be monitored is completed based on the unwrapping result corresponding to each interference pair.

6. A minimum cost flow phase unwrapping device based on phase gradient rate assistance, characterized in that: include: Data processing module, used to obtain SAR data set and DEM data of the area to be monitored; and determining an initial SAR image set based on the SAR data set; wherein the initial SAR image set includes a plurality of SAR images of the area to be monitored; determining a SAR main image in the initial SAR image set, taking the SAR main image as a reference, registering other SAR images in the initial SAR image set except the SAR main image according to the DEM data, and determining a target SAR image set based on the registration result; wherein the target SAR image set includes a plurality of target SAR images; determining a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; wherein the baseline set includes a plurality of interference pairs, each of which includes two target SAR images; and generating a differential interference map corresponding to each interference pair based on the DEM data; and generating a coherence map corresponding to each interference pair based on the differential interference map corresponding to each interference pair; A gradient map determination module, for calculating the phase gradient values ​​corresponding to each pixel in the differential interference map corresponding to each interference pair based on a phase gradient calculation formula, and determining 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 interference pair, determining a gradient winding result according to the first phase gradient map corresponding to the differential interference map based on a winding formula, and masking the gradient winding result based on a coherence map corresponding to the differential interference map to obtain a second phase gradient map; wherein the second phase gradient map includes a second horizontal direction phase gradient map, a second vertical direction phase gradient map, a third diagonal direction phase gradient map, and a fourth diagonal direction phase gradient map; An atlas determination module is used to determine a first directional phase gradient atlas based on a second horizontal directional phase gradient map corresponding to each interference pair in the baseline set, determine a second directional phase gradient atlas based on a second vertical directional phase gradient map corresponding to each interference pair in the baseline set, determine a third directional phase gradient atlas based on a third diagonal directional phase gradient map corresponding to each interference pair in the baseline set, and determine a fourth directional phase gradient atlas based on a fourth diagonal directional phase gradient map corresponding to each interference pair in the baseline set; and, based on a stacking formula, respectively calculate a first directional shortest time baseline stacking map and a first directional full time baseline stacking map according to the first directional phase gradient atlas; based on a stacking formula, respectively calculate a second directional shortest time baseline stacking map and a second directional full time baseline stacking map according to the second directional phase gradient atlas; based on a stacking formula, respectively calculate a third directional shortest time baseline stacking map and a third directional full time baseline stacking map according to the third directional phase gradient atlas; and based on a stacking formula, respectively calculate a fourth directional shortest time baseline stacking map and a fourth directional full time baseline stacking map according to the fourth directional phase gradient atlas; a rate calculation module, configured to calculate the shortest time baseline phase gradient rate based on the gradient fusion formula according to the shortest time baseline stacking diagram in the first direction, the shortest time baseline stacking diagram in the second direction, the shortest time baseline stacking diagram in the third direction, and the shortest time baseline stacking diagram in the fourth direction; and calculate the full time baseline phase gradient rate based on the gradient fusion formula according to the full time baseline stacking diagram in the first direction, the full time baseline stacking diagram in the second direction, the full time baseline stacking diagram in the third direction, and the full time baseline stacking diagram in the fourth direction; A weight determination module is used to construct a Delaunay triangulation network based on the differential interference graph and the coherence graph, and determine the irrotational constraint and the design matrix based on the Delaunay triangulation network; weight each triangle arc segment in the Delaunay triangulation network based on the coherence graph, the shortest time baseline phase gradient rate and the full time baseline phase gradient rate to obtain the weight corresponding to each triangle arc segment in the Delaunay triangulation network; The phase unwrapping module is used to solve the integer ambiguity based on the irrotational constraint, the design matrix and the weights corresponding to each triangular arc segment in the Delaunay triangulation network, and complete the phase unwrapping task of the area to be monitored based on the integer ambiguity.

7. A storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 5 is implemented.

8. A computer device comprising a storage medium, a processor, and a computer program stored in the storage medium and executable on the processor, characterized in that: When the processor executes the computer program, the method according to any one of claims 1 to 5 is implemented.

Images