Minimum cost flow phase unwrapping method and device based on phase gradient rate assistance

Abstract

The invention provides a minimum cost flow phase unwrapping method and device based on phase gradient rate assistance. The method comprises the steps that an SAR data set and DEM data of a to-be-monitored area are acquired; based on this, a differential interferogram and a coherence map are generated; respectively calculating a shortest time baseline phase gradient rate and a full time baseline phase gradient rate based on the differential interferogram; constructing a Delaunay triangulation network based on the differential interferogram and the coherence graph, and determining a rotation-free constraint and a design matrix; on the basis of the coherence graph, the shortest time baseline phase gradient rate and the full time baseline phase gradient rate, performing weight determination on each triangular arc segment in the triangulation network to obtain a weight corresponding to each triangular arc segment; and based on the non-rotation constraint, designing a matrix and weights corresponding to the triangular arc segments to solve the integer ambiguity, and completing a phase unwrapping task based on the integer ambiguity. According to the embodiment, the absolute phase of the large-gradient deformation area is obtained, so that high-precision deformation monitoring of the earth surface is realized.

Description

Technical Field

[0001] The present disclosure relates to the technical field of large-gradient deformation monitoring of the surface, and in particular to a minimum cost flow phase unwrapping method and device based on phase gradient rate assistance. 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, mine collapses, volcanoes, earthquakes, etc. At present, SAR technology is mainly divided into two categories: phase-based SAR interferometry (InSAR) and intensity-based offset tracking. InSAR technology can accurately observe small surface deformations by acquiring two images of the same area and performing phase processing; while intensity-based offset tracking technology uses the intensity information between images to estimate the offset, which is suitable for large gradient deformation monitoring, but the accuracy is relatively low.

[0003] In 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. The algorithm transforms the phase unwrapping problem into a network flow problem with the minimum computational cost, and achieves phase recovery through optimization. However, the MCF algorithm relies on the Itoh condition, that is, the phase gradient difference between adjacent pixels must meet a certain range. However, in practical applications, due to the limitation of SAR satellite sampling resolution, large gradient deformation of the surface often causes the Itoh condition to fail, making the MCF algorithm unable to effectively recover the absolute phase of the large gradient deformation area, thereby limiting the application effect of InSAR technology in complex deformation scenes. Summary of the invention

[0004] The embodiments of the present disclosure at least provide a minimum cost flow phase unwrapping method and device based on phase gradient rate assistance, which obtains the absolute phase of the large gradient deformation area, thereby achieving high-precision deformation monitoring of the surface.

[0005] The embodiment of the present disclosure provides a minimum cost flow phase unwrapping method based on phase gradient rate assistance, including: Acquire a SAR data set and DEM data of the area to be monitored; and determine an initial SAR image set based on the SAR data set; wherein the initial SAR image set includes a plurality of SAR images of the area to be monitored; determine a SAR main image in the initial SAR image set, and register other SAR images in the initial SAR image set except the SAR main image based on the DEM data, taking the SAR main image as a reference, and determine a target SAR image set based on the registration result; wherein the target SAR image set includes a plurality of target SAR images; determine a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; wherein the baseline set includes a plurality of interference pairs, each of which includes two target SAR images; and, generate a differential interference map corresponding to each interference pair based on the DEM data; and generate 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.

[0006] The embodiment of the present disclosure provides a minimum cost flow phase unwrapping device based on phase gradient rate assistance, comprising: A data processing module, used for acquiring a SAR data set and a DEM data of an 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.

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

[0008] An embodiment of the present disclosure provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the minimum cost flow phase unwrapping method based on phase gradient rate assistance as described in any possible implementation manner described above is implemented.

[0009] The phase gradient rate-assisted minimum cost flow phase unwrapping method and device provided in the embodiments of the present invention determine the boundary of large gradient deformation of the surface by utilizing the phase gradient information in the SAR image interferogram that is not affected by the phase unwrapping error, optimize the Delaunay triangulation, and control the unwrapping path to enter the large gradient deformation field from the position with a smaller deformation gradient for unwrapping, so as to correctly restore the absolute phase of the large gradient deformation field of the surface, thereby realizing high-precision deformation monitoring of the surface.

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

[0011] In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the following briefly introduces the drawings required to be cited in the embodiments. The drawings herein are incorporated into the specification and constitute a part of the specification. These drawings illustrate embodiments consistent with 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 illustrate certain embodiments of the present disclosure and should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without creative work.

[0012] Figure 1 A flow chart of a minimum cost flow phase unwrapping method based on phase gradient rate assistance provided by an embodiment of the present disclosure is shown; Figure 2 A time-space 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 is shown; Figure 3 A schematic diagram showing a surface deformation field of a simulation experiment provided by an embodiment of the present disclosure is shown; Figure 4An original interference pattern and a tangled interference pattern generated by a simulation experiment provided by an embodiment of the present disclosure are shown; Figure 5 A comparison diagram of MCF phase unwrapping results of a simulation experiment provided by an embodiment of the present disclosure is shown; Figure 6 A comparison diagram of phase unwrapping results of a method provided by the present disclosure in a simulation experiment provided by an embodiment of the present disclosure is shown; Figure 7 A schematic structural diagram of a phase gradient rate-assisted minimum cost flow phase unwrapping device provided in an embodiment of the present disclosure is shown; Figure 8 A schematic diagram of the structure of a computer device provided by an embodiment of the present disclosure is shown. DETAILED DESCRIPTION

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

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

[0015] The term "and / or" herein only describes an association relationship, indicating that three relationships may exist. For example, A and / or B may represent the following three situations: A exists alone, A and B exist at the same time, and B exists alone. In addition, the term "at least one" herein represents any combination of at least two of any one or more of a plurality of. For example, including at least one of A, B, and C may represent including any one or more elements selected from the set consisting of A, B, and C.

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

[0017] The following is a detailed description of the phase gradient rate-assisted minimum cost flow phase unwrapping method provided in the embodiment of the present application in conjunction with the accompanying drawings. Figure 1 As shown, the method includes the following S101~S106: S101, acquiring a SAR data set and a DEM data of a to-be-monitored area; and determining an initial SAR image set based on the SAR data set; determining a SAR main image in the initial SAR image set, taking the SAR main image as a reference, and 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; determining a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; and, generating a differential interferogram corresponding to each interference pair based on the DEM data; and generating a coherence map corresponding to each interference pair based on the differential interferogram corresponding to each interference pair.

[0018] It is understandable that SAR (Synthetic Aperture Radar) is a technology that uses electromagnetic waves to image the ground. It can obtain surface information all day and all weather, and is widely used in the fields of surface change monitoring and disaster warning. SAR datasets are image data collected by radar sensors that contain electromagnetic wave reflection information of the area to be monitored, which can provide accurate surface morphology and deformation information. DEM (Digital Elevation Model) is a digital model that describes the ground elevation distribution obtained through remote sensing technology or ground measurement. DEM data is used to assist in understanding the terrain undulations, and is combined with SAR data through differential interferograms to derive more accurate phase unwrapping results. 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 change monitoring is required, such as high-risk areas such as earthquakes, landslides or urban changes.

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

[0020] 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 accurately match in space. In the registration process, DEM data is used to correct the geometric errors of the image to ensure the consistency of the image in space. Finally, a target SAR image set is generated through these registered SAR images and SAR master images for subsequent analysis. Among them, the target SAR image set includes multiple target SAR images.

[0021] It can be understood that the spatiotemporal baseline threshold refers to the limitation on the validity between image pairs in time and space, and is usually a parameter value used to screen out interference pairs that are representative and relatively consistent. In the present disclosure, the spatiotemporal baseline threshold is set to a time of less than 36 days and a spatial distance of less than 150m; in some other embodiments, it can be set according to specific needs and is not specifically limited here. Based on the target SAR image set, a baseline set is determined by a preset spatiotemporal baseline threshold. The baseline set includes multiple interference pairs, and each interference pair includes two target SAR images. Here, an interference pair refers to the pairing of two images required for interference analysis using SAR images. These paired images are usually separated by a certain time and space and can reflect ground changes.

[0022] Specifically, an interferogram is formed by two radar observations (or images), and the time and location of each observation may be slightly different. By comparing the two observation images, an interferogram reflecting surface changes (such as settlement, displacement, etc.) can be obtained. The differential interferogram is obtained by subtracting the terrain error (usually corrected using DEM data), which can show detailed information about the relative changes in the surface. The differential interferogram of each interferogram pair can reveal the phase difference between the two observations, which reflects the information of ground displacement or deformation.

[0023] It can be understood that the coherence map reflects the signal quality and correlation of each pixel in the interference map. Coherence essentially describes the consistency of the radar signal in time and space. In radar imaging, different features of the surface (such as vegetation, buildings, mountains, etc.) have different reflection characteristics for radar signals, which will affect the coherence of the signal. Areas with high coherence usually represent stable surface features and consistent reflections, while areas with low coherence may be due to dynamic changes in the surface, differences in scatterers, or atmospheric effects. By generating a coherence map based on the differential interferogram of each interference 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 need further correction or exclusion. For example, if the coherence of a certain area is low, it may indicate that there are large surface changes in the area, or signal distortion due to atmospheric factors.

[0024] Among them, determining the coherence map corresponding to each interference pair based on the differential interference map corresponding to each interference pair can be achieved by a coherence calculation algorithm. This algorithm usually involves statistical analysis of the phase information in the differential interference map to quantify the phase stability of adjacent pixels or the same pixel at different times. Specifically, the coherence formula for calculating the standard deviation of the phase difference or the change of the phase gradient can be included to calculate the coherence value of each pixel, thereby obtaining a coherence map, in which the value of each pixel represents the degree of signal coherence at that position in the differential interference map.

[0025] S102, for the differential interference pattern corresponding to each interference pair, calculate the phase gradient value corresponding to each pixel in the differential interference pattern based on the phase gradient calculation formula, and determine the first phase gradient pattern based on the phase gradient value corresponding to each pixel; for each interference pair, determine the gradient winding result according to the first phase gradient pattern corresponding to the differential interference pattern based on the winding formula, and mask the gradient winding result based on the coherence map corresponding to the differential interference pattern to obtain the second phase gradient pattern.

[0026] Specifically, for the differential interference image corresponding to each interference pair, based on the calculation formula of the phase gradient, the phase gradient value of each pixel in the image is determined, and these gradient values ​​reflect the rate of change of the phase difference between adjacent pixels. The first phase gradient image can be determined by the calculated phase gradient value; 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.

[0027] Here, the phase gradient calculation formula can be expressed as: ; 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; Expressed as the phase gradient in the second diagonal direction.

[0028] It can be understood that, for each interference pair, the gradient winding result is first determined based on the winding formula and the first phase gradient map. The phase winding phenomenon occurs when the phase gradient value exceeds a certain range, so it is necessary to use the winding formula to correct the phase gradient to avoid calculation errors caused by winding. At this time, there will be random phase gradients caused by random noise and local tiny gradients caused by atmospheric errors in the phase gradient. In order to solve this phenomenon, the gradient winding result can be masked in combination with the coherence map to obtain a second phase gradient map. In this way, the calculation accuracy is improved by retaining valid gradient information in areas with higher coherence and removing invalid or erroneous gradient information in areas with lower coherence. Among them, the second phase gradient map includes a second horizontal phase gradient map, a second vertical phase gradient map, a third diagonal phase gradient map, and a fourth diagonal phase gradient map; Here, the winding formula can be expressed as: ; in, It is represented as the gradient winding result; It is expressed as phase gradient.

[0029] S103, determining a first directional phase gradient atlas based on the second horizontal directional phase gradient map corresponding to each interference pair in the baseline set, determining a second directional phase gradient atlas based on the second vertical directional phase gradient map corresponding to each interference pair in the baseline set, determining a third directional phase gradient atlas based on the third diagonal directional phase gradient map corresponding to each interference pair in the baseline set, and determining a fourth directional phase gradient atlas based on the fourth diagonal directional phase gradient map corresponding to each interference pair in the baseline set; and, calculating a first direction shortest time baseline stacking map and a first direction full time baseline stacking map respectively according to the first directional phase gradient atlas based on a stacking formula; calculating a second direction shortest time baseline stacking map and a second direction full time baseline stacking map respectively according to the second directional phase gradient atlas based on a stacking formula; calculating a third direction shortest time baseline stacking map and a third direction full time baseline stacking map respectively according to the third directional phase gradient atlas based on a stacking formula; calculating a fourth direction shortest time baseline stacking map and a fourth direction full time baseline stacking map respectively according to the fourth directional phase gradient atlas based on a stacking formula.

[0030] Specifically, for different directional phase gradient atlases, the phase gradient atlases in different directions are determined based on the second phase gradient maps corresponding to each interference pair. For example, based on the phase gradient map in the horizontal direction, the first direction phase gradient atlas of all interference pairs is obtained; based on the phase gradient map in the vertical direction, the second direction phase gradient atlas is obtained, and so on. In this way, the data of multiple interference pairs can be accumulated for each direction, so as to better capture the phase change characteristics in different directions.

[0031] It can be understood that the phase gradient atlas in each direction is stacked through the stacking formula, and the shortest time baseline stacking diagram and the full time baseline stacking diagram corresponding to each direction are calculated. 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 stacking diagram focuses on the phase changes in a short time, while the time baseline stacking diagram integrates the phase changes in a long time.

[0032] Here, the stacking formula can be expressed as: ; 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.

[0033] S104, calculating the shortest time baseline phase gradient rate based on the gradient fusion formula according to the shortest time baseline stacked graph in the first direction, the shortest time baseline stacked graph in the second direction, the shortest time baseline stacked graph in the third direction, and the shortest time baseline stacked graph in the fourth direction; calculating the full time baseline phase gradient rate based on the gradient fusion formula according to the full time baseline stacked graph in the first direction, the full time baseline stacked graph in the second direction, the full time baseline stacked graph in the third direction, and the full time baseline stacked graph in the fourth direction.

[0034] It is understood that the phase gradient rate refers to the rate of change of the phase difference between adjacent pixels, which is used to describe the spatial variation 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.

[0035] Specifically, the shortest time baseline phase gradient rate is calculated by combining the shortest time baseline stacking diagram from the first direction to the fourth direction through the gradient fusion formula. 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.

[0036] Here, the gradient fusion formula can be expressed as: ; 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 .

[0037] Similarly, based on the full-time baseline stacking diagram from the first direction to the fourth direction, the full-time baseline phase gradient rate is calculated using the above-mentioned gradient fusion formula to reflect the rate of phase change on a long-term time scale.

[0038] S105, constructing a Delaunay triangulation network based on the differential interference graph and the coherence graph, and determining the irrotational constraint and the design matrix based on the Delaunay triangulation network; weighting each triangular 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, and obtaining the weight corresponding to each triangular arc segment in the Delaunay triangulation network.

[0039] It is understandable that the Delaunay triangulation is a triangular mesh structure generated by a set of scattered points. It has good mathematical properties and is widely used in the field of geographic information systems (GIS) and computational geometry. The triangulation requires that the circumscribed circle of any triangle does not contain other points, thereby avoiding the appearance of "long strip" triangles. In this scheme, the Delaunay triangulation is used to divide the area to be monitored into multiple triangular areas, and further used to define the constraints of phase unwrapping.

[0040] Specifically, when constructing a Delaunay triangulation based on the differential interference graph and the coherence graph, the following (1) to (2) can be included: (1) for each differential interference pattern corresponding to the interference pair, extracting high-quality points that meet a preset coherence threshold based on a coherence map corresponding to the differential interference pattern, and determining a high-quality point set corresponding to the interference pair based on the extraction result; (2) 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.

[0041] It can be understood that high-quality points that meet the preset coherence threshold are extracted through the coherence graph, 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 the point set used in the subsequent triangulation construction process has high accuracy and representativeness, and the quality of the points in the Delaunay triangulation is guaranteed. 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 network construction method, which adjusts the existing triangle structure by adding new points one by one, and finally forms a mesh that meets the Delaunay conditions.

[0042] It can be understood that the irrotation constraint means that the phase change between adjacent regions must meet certain physical consistency during the unwrapping process, that is, the unwrapping result cannot have rotation or unreasonable phase change, ensuring the continuity and physical rationality of the phase during the unwrapping process. The design matrix is ​​a mathematical tool used to convert various constraints and known conditions into computable mathematical forms. In the phase unwrapping problem, the design matrix is ​​used to describe the relationship between each triangle and its constraints, and ultimately help solve the overall phase.

[0043] Specifically, 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.

[0044] Here, the design matrix can be expressed as: , ; 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 .

[0045] Here, the residual calculation formula can be expressed as: ; 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 .

[0046] For example, in a Delaunay triangulation network, each arc segment of a triangle corresponds to a weight. The weight is determined based on the information in the phase gradient rate and coherence map in the differential interferogram. The phase gradient rate and coherence map provide the signal strength and speed of phase changes between different regions, thereby affecting the connection weight between adjacent regions. In the untangling process, reasonable weight distribution helps to balance the contribution of different regions to the final untangling result and improve the accuracy of untangling.

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

[0048] Here, 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; 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.

[0049] S106, solving integer ambiguity based on the irrotational constraint, the design matrix and the weights corresponding to each triangular arc segment in the Delaunay triangulation, and completing the phase unwrapping task of the area to be monitored based on the integer ambiguity.

[0050] It is understandable that integer ambiguity resolution (AR) is a common technology in satellite positioning (such as GNSS) and radar (such as SAR) systems. Especially in the unwrapping process of phase measurement, integer ambiguity refers to the uncertainty caused by periodicity in phase data, and how to accurately determine the ambiguity (i.e., integer multiples of wavelength). In SAR, due to the periodicity of phase values, multiple unambiguity problems may occur in the unwrapping process. By solving integer ambiguities based on irrotational constraints, design matrices, and calculated weights, inconsistent phase ambiguities can be eliminated, making phase unwrapping more accurate. Finally, by determining the integer ambiguity, the phase of the area to be monitored can be effectively untangled and the phase unwrapping task can be completed, thereby providing accurate terrain deformation information or other related data of other monitoring targets.

[0051] Specifically, when solving the integer ambiguity, the following (a) to (b) can be included: (a) Constructing the objective function; (b) for each interference pair corresponding to the Delaunay triangulation, solving the integer ambiguity based on the objective function 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.

[0052] Here, the objective function can be expressed as: ; Wherein, W represents a vector including weights corresponding to each triangular arc segment in the Delaunay triangulation; X represents the integer ambiguity; A represents the design matrix; and U represents the irrotational constraint.

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

[0054] The phase gradient rate-assisted minimum cost flow phase unwrapping method and device provided in the embodiments of the present invention determine the boundary of large gradient deformation of the surface by utilizing the phase gradient information in the SAR image interferogram that is not affected by the phase unwrapping error, optimize the Delaunay triangulation, and control the unwrapping path to enter the large gradient deformation field from the position with a smaller deformation gradient for unwrapping, so as to correctly restore the absolute phase of the large gradient deformation field of the surface, thereby realizing high-precision deformation monitoring of the surface.

[0055] In order to verify the performance of the phase unwrapping of the proposed method, a total of 91 scenes of Sentinel-1 down-orbit data from January 2, 2019 to December 29, 2021 were obtained for simulation experiments and real experiments. The 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 because of 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. Afterwards, a maximum time baseline of 36 days, a maximum spatial baseline of 150 meters, and a total of 267 interference pairs were generated. The time-space baseline diagram is shown as follows: Figure 2 As shown. Afterwards, simulation experiments and real data experiments were carried out based on the acquired Sentinel-1 SAR data. Among them, the deformation field of the simulation experiment is simulated, and the spatiotemporal baseline is real data. In the real experiment, the surface deformation field and spatiotemporal baseline are both real data.

[0056] Specifically, refer to Figure 3 As shown, the present disclosure simulates 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 increases slowly 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 increases slowly from top to bottom, and there are obvious deformation boundaries on the two sides and the bottom compared with the stable area, and there is no obvious deformation boundary on the top; there is an obvious deformation boundary between the rectangular deformation field and the stable area. Afterwards, 264 original interference patterns were generated according to the space-time baseline, and then the phase of the original interference pattern was twisted between them to obtain a twisted interference pattern. When generating the simulated interference pattern, a periodic deformation with an annual period is added on the basis of the linear deformation. The interference phase diagrams of the first three scenes generated by simulation are shown as follows. Figure 4 As shown, (a)-(c) are the simulated original interference phase diagrams, and their time baselines are 12 days, 24 days and 36 days respectively. (d)-(f) are the interference diagrams obtained by winding (a)-(c). Figure 2 It can be seen that the phase of the entanglement interference pattern is between , and when the time baseline is 24 days and 36 days, phase jumps appear inside the circular deformation field, inside and at the boundary of the trapezoidal deformation field, and at the boundary of the rectangular deformation field.

[0057] After generating the simulated winding phase diagram, phase unwrapping was performed using the traditional MCF method and the method provided by the embodiment of the present disclosure. The unwrapping results of the traditional MCF method are as follows: Figure 5 shown. Figure 5 (a) is the unwrapped phase image obtained by the MCF method from the wrapped phase image with a time baseline of 12 days. Similarly, Figure 5 (b) and Figure 4(c) The unwrapped phase images obtained by the MCF method from the 24-day and 36-day wrapped phase images, respectively. Figure 5 (d)-(f) are the differences between the unwrapped phase image and the original phase image of the same time baseline. Figure 5 (a) It can be seen that when the time baseline is 12 days, the phase is correctly restored and no unwrapping error occurs. Figure 5 (d) only contains errors caused by computer storage accuracy. When the time baseline is 24 days, the circular deformation area obtains correct unwrapping results, part of the trapezoidal deformation area obtains correct unwrapping results, and all the rectangular deformation areas obtain incorrect unwrapping results. When the time baseline is extended to 36 days, the circular deformation area still obtains correct phase unwrapping results, while the area with correct results in the trapezoidal deformation area decreases, and the rectangular deformation area still obtains all incorrect phase unwrapping results.

[0058] The untangling results of the method provided in the embodiment of the present disclosure are as follows: Figure 6 Shown Figure 6 (a) is the unwrapped phase image obtained by the proposed phase unwrapping method from the wrapped phase image with a time baseline of 12 days. Similarly, Figure 6 (b) and Figure 6 (c) Unwrapped phase images obtained by phase unwrapping the 24-day and 36-day wrapped phase images using the proposed method. Figure 6 (d)-(f) are the differences between the unwrapped phase image and the simulated original phase image of the same time baseline. It can be seen from the figure that when the time baseline is 12 days, the phase is correctly restored and no unwrapping error occurs. Figure 6 (d) only contains errors caused by computer storage accuracy. When the time baseline is 24 days and 36 days, the circular deformation area obtains correct unwrapping results, and most of the trapezoidal deformation area obtains correct unwrapping results, with only scattered points at the edge of the deformation field that are not correctly phase unwrapped. However, the rectangular deformation area still obtains incorrect unwrapping results. This is because the edge of the rectangular deformation field has the phenomenon of winding phase gradient disappearance when the time baseline is greater than 12 days, which makes it impossible for the unwrapping network to enter the deformation field from the boundary with a smaller deformation gradient to unwrap.

[0059] Those skilled in the art will appreciate that, in the above method of specific implementation, the order in which the steps are written does not imply a strict execution order and does not constitute any limitation on the implementation process. The specific execution order of the steps should be determined by their functions and possible internal logic.

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

[0061] Reference Figure 7 FIG. 1 is a schematic diagram of a phase gradient rate-assisted minimum cost flow phase unwrapping device 700 provided in an embodiment of the present disclosure, wherein the device comprises: The data processing module 701 is used to obtain a SAR data set and a DEM data of a region to be monitored; and determine an initial SAR image set based on the SAR data set; wherein the initial SAR image set includes a plurality of SAR images of the region to be monitored; determine a SAR main image in the initial SAR image set, and register other SAR images in the initial SAR image set except the SAR main image based on the DEM data, taking the SAR main image as a reference, and determine a target SAR image set based on the registration result; wherein the target SAR image set includes a plurality of target SAR images; determine a baseline set based on the target SAR image set according to a preset spatiotemporal baseline threshold; wherein the baseline set includes a plurality of interference pairs, each of which includes two target SAR images; and, generate a differential interference map corresponding to each interference pair based on the DEM data; and generate a coherence map corresponding to each interference pair based on the differential interference map corresponding to each interference pair; A gradient map determination module 702 is used to calculate the phase gradient value corresponding to each pixel in the differential interference map corresponding to each interference pair based on the 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, a gradient winding result is determined according to the first phase gradient map corresponding to the differential interference map based on the winding formula, and the gradient winding result is masked based on the 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; The atlas determination module 703 is used to determine the first direction phase gradient atlas based on the second horizontal direction phase gradient map corresponding to each interference pair in the baseline set, determine the second direction phase gradient atlas based on the second vertical direction phase gradient map corresponding to each interference pair in the baseline set, determine the 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 the fourth direction phase gradient atlas based on the fourth diagonal direction phase gradient map corresponding to each interference pair in the baseline set; and, based on the stacking formula, respectively calculate the first direction shortest time baseline stacking map and the first direction full time baseline stacking map according to the first direction phase gradient atlas; based on the stacking formula, respectively calculate the second direction shortest time baseline stacking map and the second direction full time baseline stacking map according to the second direction phase gradient atlas; based on the stacking formula, respectively calculate the third direction shortest time baseline stacking map and the third direction full time baseline stacking map according to the third direction phase gradient atlas; based on the stacking formula, respectively calculate the fourth direction shortest time baseline stacking map and the fourth direction full time baseline stacking map according to the fourth direction phase gradient atlas; The rate calculation module 704 is used 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; 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 705 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 706 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 monitored area based on the integer ambiguity.

[0062] In some possible embodiments, 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, It is represented as the 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 .

[0063] In some possible embodiments, the weight determination module 705 is specifically used to: 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 an 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.

[0064] In some possible embodiments, the weight determination module 705 is specifically used to: 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 weight determination module 705 is specifically used to: for the Delaunay triangulation corresponding to each interference pair, determine the weight of each triangular arc segment in the Delaunay triangulation corresponding to each interference pair based on the arc segment weight calculation formula 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; 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.

[0065] In some possible embodiments, the phase unwrapping module 706 is specifically configured to: 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; X represents the integer ambiguity; A represents the design matrix; U represents the irrotational constraint; The phase unwrapping module 706 is specifically used to: integrate the integer ambiguity corresponding to each interference pair to obtain the unwrapping result corresponding to each interference pair; and complete the phase unwrapping task of the monitored area based on the unwrapping result corresponding to each interference pair.

[0066] Based on the same technical concept, the embodiment of the present disclosure also provides a computer device. Figure 8 , which is a schematic diagram of the structure of a computer device 800 provided in an embodiment of the present disclosure, including a processor 801, a memory 802, and a bus 803. The memory 802 is used to store execution instructions, including a memory 8021 and an external memory 8022; the memory 8021 is also called an internal 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 memory 8021.

[0067] In the embodiment of the present application, the memory 802 is specifically used to store the application code for executing the solution of the present application, and the execution is controlled by the processor 801. That is, when the computer device 800 is running, the processor 801 communicates with the memory 802 through the bus 803, so that the processor 801 executes the application code stored in the memory 802, and then executes the method described in any of the above embodiments. Among them, the memory 802 can be, but 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 device, a discrete gate or transistor logic device, a discrete hardware component. The disclosed methods, steps and logic block diagrams in the embodiments of the present invention can be implemented or executed. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc.

[0068] It is to be understood that the structure illustrated 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 some components, or split some components, or arrange the components differently. The components shown in the figure may be implemented in hardware, software, or a combination of software and hardware.

[0069] The embodiment of the present disclosure also provides a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the phase gradient rate-assisted minimum cost flow phase unwrapping method described in the above method embodiment are executed. The storage medium can be a volatile or non-volatile computer-readable storage medium.

[0070] The embodiment of the present disclosure also provides a computer program product, which carries a program code, and the instructions included in the program code 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 embodiment. For details, please refer to the above method embodiment, which will not be repeated here. Among them, the above computer program product can be implemented specifically by hardware, software or a combination thereof. In an optional embodiment, the computer program product is specifically embodied as a computer storage medium, and in another optional embodiment, the computer program product is specifically embodied as a software product, such as a software development kit, etc. A person skilled in the art can clearly understand that for the convenience and simplicity of description, the specific working process of the system and device described above can refer to the corresponding process in the above method embodiment, which will not be repeated here. In several embodiments provided in the present disclosure, it should be understood that the disclosed system and method can be implemented in other ways. The device embodiments described above are only schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation. 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 mutual coupling or direct coupling or communication connection shown or discussed may be an indirect coupling or communication connection through some communication interface, device or unit, which may be electrical, mechanical or other forms.

[0071] The unit described as a separate component may or may not be physically separated, and the component displayed as a unit may or may not be a physical unit, that is, it may be located in one place, or it may be distributed on multiple network units. Some or all of the units may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. In addition, the functional units in each embodiment of the present disclosure may be integrated into a processing unit, or each unit may exist physically alone, or two or more units may be integrated into one unit. If the function is implemented in the form of a software functional unit and sold or used as an independent product, it may be stored in a non-volatile computer-readable storage medium executable by a processor. Based on this understanding, the technical solution of the present disclosure is essentially or the part that contributes to the prior art or the part of the technical solution can be embodied in the form of a software product, which is stored in a storage medium, including several instructions for a computer device to perform all or part of the steps of the method described in each embodiment of the present disclosure. The aforementioned storage medium includes: various media that can store program codes, such as a U disk, a mobile hard disk, a read-only memory, a random access memory, a disk or an optical disk.

[0072] Finally, it should be noted that the above-described embodiments are only specific implementation methods of the present disclosure, which are used to illustrate the technical solutions of the present disclosure, rather than to limit them. The protection scope of the present disclosure is not limited thereto. Although the present disclosure is described in detail with reference to the above-described embodiments, ordinary technicians in the field should understand that any technician familiar with the technical field can still modify the technical solutions recorded in the above-described embodiments within the technical scope disclosed in the present disclosure, or can easily think of changes, or make equivalent replacements for some of the technical features therein; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present disclosure, and should be included in the protection scope of the present disclosure. Therefore, the protection scope of the present disclosure should be based on 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 an 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; 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; 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