A co-seismic interference pattern unwrapping method, device and medium

Abstract

The present application relates to the field of surface deformation monitoring technology, and provides a method, device and medium for unwrapping a coseismic interferogram. The method includes: obtaining an interferogram and a distance offset map using SAR images; performing median filtering on the distance offset map to obtain the initial offset and offset standard deviation of each pixel; determining a special central pixel based on all offset standard deviations, and calculating the final offset of the special central pixel based on the initial offset; obtaining the initial prior gradient and the significance of the pixel based on the final offset; determining a significant pixel based on the significance, and obtaining multiple connected regions based on the significant pixel; determining the near field and the far field from the connected region and obtaining the final prior gradient based on the near field and the far field; constructing an objective function using the final prior gradient and solving to obtain the unwrapping phase of the interferogram gradient; performing path integration on the unwrapping phase to obtain the interferogram unwrapping result. The method of the present application can improve the accuracy of the unwrapping result.

Description

Technical Field

[0001] The present application relates to the technical field of surface deformation monitoring, and in particular to a co-seismic interference pattern unwrapping method, equipment and medium. Background Art

[0002] Spaceborne synthetic aperture radar interferometry (InSAR) has the advantages of wide coverage, high spatial resolution, and high accuracy in monitoring surface deformation. It is a common space geodetic technology for obtaining earthquake deformation. The InSAR co-seismic deformation field of medium and large earthquakes can provide important data support for identifying surface ruptures, inverting fault source parameters and slip distribution, studying seismogenic mechanisms, and assessing regional earthquake risks. However, the imaging data of the commonly used C-band synthetic aperture radar (SAR) satellites (such as Sentinel-1) in the world are very sensitive to large gradient deformation and surface scattering changes, resulting in a large number of aliasing fringes and incoherence noise in the deformation near field of the InSAR co-seismic interferogram, which brings great difficulties to phase unwrapping. Phase unwrapping directly affects the accuracy of the InSAR deformation field, restricting the interpretation and application research of InSAR technology in earthquake events.

[0003] InSAR phase unwrapping refers to the process of restoring the wrapped phase from (-π,π) to the absolute unwrapped phase, that is, adding several integer multiples of 2π cycles to the wrapped phase of each pixel. Traditional unwrapping methods take into account the continuity of surface deformation and introduce the phase continuity assumption (that is, assuming that the absolute value of the difference in the unwrapped phase of adjacent pixels does not exceed π) as a constraint to complete unwrapping. Affected by factors such as SAR imaging geometry, spatial sampling rate and decoherence noise, not all pixels meet this assumption, which will cause unwrapping errors and their propagation. Minimum Cost Flow (MCF) and Statistical-Cost Network-Flow Algorithm for Phase Unwrapping (SNAPHU) are the two most commonly used robust unwrapping methods. They complete global phase unwrapping by minimizing the L1 norm of the difference between the wrapped and unwrapped gradients and maximizing the conditional probability density function based on statistical information, respectively, which can effectively suppress the large-scale propagation of unwrapping errors. However, the near-field deformation gradient of the C-band InSAR coseismic interferogram of medium and large earthquakes is large, the changes are drastic, and the decoherence phenomenon is serious. The MCF and SNAPHU methods, which essentially rely on the phase continuity assumption, are difficult to complete the correct unwrapping, resulting in the absolute value of the near-field deformation being underestimated, and the far-field deformation is also disturbed by the propagation of the unwrapping error.

[0004] In the production of actual InSAR deformation fields, masking the near-field decoherence area in the interferogram so that it does not participate in the unwrapping is a common avoidance solution. The resulting lack of near-field deformation data will affect the characterization and analysis of the slip behavior of shallow faults in large strike-slip earthquakes, and may even introduce the illusion of slip loss during the inversion of fault models. How to choose the appropriate coherence mask threshold or draw the mask area is also a problem that requires repeated attempts and considerations in this solution.

[0005] In addition, in order to solve the problem of unwrapping large deformation gradient interferograms, some scholars preprocessed the data based on the classical unwrapping method. Before unwrapping, the co-seismic interferogram was conjugate multiplied with the prior rough deformation field of the entanglement (derived from the range offset, the initial unwrapping result or the forward result of the physical model) to weaken the fringe gradient and fringe aliasing, so that the remaining entangled phase is closer to the phase continuity assumption. Other scholars have tried to post-process the unwrapping results to improve the accuracy, such as manually correcting some errors based on the range offset, and fusing the range offset with the unwrapping results by principal component analysis. However, such data pre-processing or post-processing strategies do not fundamentally change the dependence on the phase continuity assumption in the unwrapping process, so the improvement in the accuracy of the final unwrapping result is limited. It can be seen that there is a problem of low accuracy in unwrapping the interferogram of seismic events. Summary of the invention

[0006] The present application provides a co-seismic interference pattern unwrapping method, device and medium, which can solve the problem of low accuracy in unwrapping the interference pattern of seismic events.

[0007] In a first aspect, an embodiment of the present application provides a co-seismic interference pattern unwrapping method, the co-seismic interference pattern unwrapping method comprising:

[0008] The first SAR image and the second SAR image of the target area are obtained, and the interference map and the range offset map of the target area are obtained by using the first SAR image and the second SAR image; the time point corresponding to the first SAR image is the time when the satellite last passes before the earthquake event occurs in the target area, and the time point corresponding to the second SAR image is the time when the satellite first passes after the earthquake event occurs;

[0009] Perform median filtering on the range offset map to obtain the initial offset of each pixel in the range offset map, and calculate the offset standard deviation of each pixel based on all the initial offsets;

[0010] According to all the offset standard deviations, a plurality of special center pixels are determined from all the pixels in the distance offset map, and a final offset of each special center pixel is calculated based on all the initial offsets; the offset standard deviation of the special center pixel is greater than a standard deviation threshold;

[0011] Based on the final offsets of all special center pixels and the initial offsets of other pixels, multiple initial prior gradients and the significance of each pixel in the distance offset map are obtained; the initial prior gradient is used to describe the horizontal or vertical deformation gradient of the pixel, and the significance is used to describe the degree of difference between the pixel and all other pixels;

[0012] According to all the saliencies, multiple salient pixels with saliency greater than a saliency threshold are determined from the distance-to-offset map, and the regions formed by the adjacent multiple salient pixels are all regarded as connected regions;

[0013] The connected area with the largest area is taken as the near-field area, and each other connected area is taken as the far-field area. Based on all the initial prior gradients, all the far-field areas and the near-field areas, multiple final prior gradients are obtained; the final prior gradients are used to describe the deformation of pixels in the near-field area or the far-field area;

[0014] An objective function is constructed using all final prior gradients, and the objective function is maximized to obtain an unwrapped phase of each interferogram gradient; the objective function is used to describe the probability of the unwrapped phase of all interferogram gradients, and the interferogram gradient is a vector obtained by subtracting the entangled phase of each pixel in the interferogram from that of adjacent pixels in the horizontal or vertical direction;

[0015] Path integration is performed on all unwrapped phases to obtain an interference pattern unwrapping result of the target area; the interference pattern unwrapping result is used to describe the seismic deformation of the target area.

[0016] Optionally, obtaining an interference map and a range offset map of the target area using the first SAR image and the second SAR image includes:

[0017] Performing differential interferometry on the first SAR image and the second SAR image to obtain an interferogram of the target area;

[0018] The offset of the first SAR image and the second SAR image is tracked to obtain a range offset map of the target area.

[0019] Optionally, calculate the final offset for each special center pixel based on all initial offsets, including:

[0020] For each special central pixel, perform the following steps:

[0021] Using the maximum inter-class variance threshold segmentation method, other pixels in a window of a preset scale centered on the special central pixel are divided into two categories;

[0022] The difference between the average of the initial offsets of the two types of pixels and the initial offset of the special central pixel is calculated respectively, and each pixel included in the pixel type with the smaller difference is used as a homogeneous pixel corresponding to the special central pixel;

[0023] The median of the initial offsets of all homogeneous pixels is calculated and used as the final offset of the special central pixel.

[0024] Optionally, based on the final offsets of all special center pixels and the initial offsets of other pixels, multiple initial prior gradients and the saliency of each pixel are obtained, including:

[0025] Take the initial offset of each other pixel as the final offset of the other pixels;

[0026] For each pixel, the final offset of the pixel is subtracted from the final offset of the next pixel in the horizontal direction to obtain the initial prior gradient of the pixel in the horizontal direction, and the final offset of the pixel is subtracted from the final offset of the next pixel in the vertical direction to obtain the initial prior gradient of the pixel in the vertical direction;

[0027] Based on the final offsets of all pixels, the saliency of each pixel is calculated.

[0028] Optionally, based on all initial prior gradients, all far-field regions, and near-field regions, multiple final prior gradients are obtained, including:

[0029] For each far-field region, a plurality of candidate prior gradients whose absolute values ​​are greater than π are determined from all initial prior gradients corresponding to all pixels in the far-field region;

[0030] For each candidate prior gradient, if the candidate prior gradient is negative, -π is used as the final prior gradient corresponding to the candidate prior gradient; if the candidate prior gradient is positive, π is used as the final prior gradient corresponding to the candidate prior gradient;

[0031] Each initial prior gradient of the non-candidate prior gradient in each far-field region is taken as the final prior gradient, and each initial prior gradient in the near-field region is taken as the final prior gradient.

[0032] Optionally, construct the objective function using all final prior gradients, including:

[0033] Constructing the unwrapping probability density function of each interferogram gradient using all final prior gradients;

[0034] For each interference pattern gradient, based on the unwrapping probability density function corresponding to the interference pattern gradient, a conditional probability density function corresponding to the interference pattern gradient is obtained;

[0035] Construct the objective function based on all conditional probability density functions.

[0036] Optionally, the unwrapping probability density function is:

[0037]

[0038] in, represents the unwrapping probability density of the j-th interference pattern gradient, represents the unwrapped phase to be solved for the jth interferogram gradient, j = 1, 2, ..., k, k represents the total number of interferogram gradients, σ ​​represents the coherence parameter, G RO represents all final prior gradients, G INF represents the winding gradient corrected mean, T coh represents the coherence threshold, and coherence represents the coherence of the j-th interference pattern gradient;

[0039] The conditional probability density function is:

[0040]

[0041] in, represents the conditional probability density of the j-th interference pattern gradient, represents the winding phase of the j-th interference pattern gradient, and a and b represent arbitrary integers.

[0042] Optionally, the objective function is:

[0043]

[0044] Wherein, F(ΔΦ|ΔΨ) represents the objective function value, ΔΦ represents the unwrapped phase of all interferogram gradients, and ΔΨ represents the wrapped phase of all interferogram gradients.

[0045] In a second aspect, an embodiment of the present application provides a co-seismic interference pattern unwrapping device, comprising:

[0046] An acquisition module acquires a first SAR image and a second SAR image of a target area, and uses the first SAR image and the second SAR image to obtain an interference diagram and a range offset diagram of the target area; the time point corresponding to the first SAR image is the time when the satellite last passes before an earthquake event occurs in the target area, and the time point corresponding to the second SAR image is the time when the satellite first passes after the earthquake event occurs;

[0047] The calculation module performs median filtering on the distance offset map to obtain the initial offset of each pixel in the distance offset map, and calculates the offset standard deviation of each pixel based on all the initial offsets;

[0048] A first determination module determines a plurality of special center pixels from the distance offset map according to all offset standard deviations, and calculates a final offset of each special center pixel based on all initial offsets; the offset standard deviation of the special center pixel is greater than a standard deviation threshold;

[0049] The saliency acquisition module obtains the saliency of each pixel in multiple initial prior gradients and the distance offset map based on the final offset of all special center pixels and the initial offset of other pixels; the initial prior gradient is used to describe the horizontal or vertical deformation gradient of the pixel, and the saliency is used to describe the degree of difference between the pixel and all other pixels;

[0050] The second determination module determines, according to all the saliencies, a plurality of salient pixels having a saliency greater than a saliency threshold from all the pixels in the distance-to-offset map, and regards the regions formed by the plurality of adjacent salient pixels as connected regions;

[0051] The final prior gradient acquisition module takes the connected area with the largest area as the near field area and each other connected area as the far field area. Based on all the initial prior gradients, all the far field areas and the near field areas, multiple final prior gradients are acquired. The final prior gradient is used to describe the deformation of pixels in the near field area or the far field area.

[0052] A construction module is provided, wherein an objective function is constructed using all final prior gradients, and the objective function is maximized to obtain an unwrapped phase of each interferogram gradient; the objective function is used to describe the probability of the unwrapped phase of all interferogram gradients, and the interferogram gradient is a vector obtained by subtracting the wrapped phase of each pixel in the interferogram from that of adjacent pixels in the horizontal or vertical direction;

[0053] The path integration module performs path integration on all unwrapped phases to obtain the interference pattern unwrapping results of the target area; the interference pattern unwrapping results are used to describe the seismic deformation of the target area.

[0054] In a third aspect, an embodiment of the present application provides a terminal device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the above-mentioned co-seismic interference pattern unwrapping method when executing the above-mentioned computer program.

[0055] In a fourth aspect, an embodiment of the present application provides a computer-readable storage medium, which stores a computer program, and when the computer program is executed by a processor, the above-mentioned co-seismic interference pattern unwrapping method is implemented.

[0056] The above solution of the present application has the following beneficial effects:

[0057] In some embodiments of the present application, by determining multiple special center pixels from all pixels in the distance offset map according to all offset standard deviations, and calculating the final offset of each special center pixel based on all distance offsets, pixels with large offset standard deviations can be processed separately, so that the accuracy of the pixel offset is improved; each significant pixel is connected to adjacent significant pixels to obtain multiple connected regions, and then the connected region with the largest area is used as the near-field region, and each of the other connected regions is used as the far-field region, and multiple final prior gradients are obtained based on all initial prior gradients, all far-field regions, and near-field regions. The near-field region and the far-field region are determined according to the significance of the pixels, and the final prior gradient is obtained according to the near-field region and the far-field region, so that different processing is performed on the near-field region and the far-field region, so that the final prior gradient conforms to the different actual conditions of the near-field region and the far-field region; the objective function is constructed by utilizing all the final prior gradients without relying on the traditional phase continuity assumption, which can improve the accuracy of the objective function and reduce the error, so that the unwrapping phase of the interference pattern gradient obtained by maximizing the objective function is highly accurate; finally, the interference pattern unwrapping result is obtained by utilizing the accurate unwrapping phase of the interference pattern gradient, which can effectively improve the accuracy of the unwrapping result.

[0058] Other beneficial effects of the present application will be described in detail in the subsequent specific implementation section. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0060] Figure 1 A flowchart of a co-seismic interference pattern unwrapping method provided in an embodiment of the present application;

[0061] Figure 2 A schematic diagram of a first disentanglement result provided in an embodiment of the present application;

[0062] Figure 3 A schematic diagram of statistical results of unwrapping residuals provided in an embodiment of the present application;

[0063] Figure 4 A schematic diagram of filtering results of the distance offset provided by an embodiment of the present application;

[0064] Figure 5 A schematic diagram of a second disentanglement result provided in an embodiment of the present application;

[0065] Figure 6A schematic diagram of LOS deformation comparison along a vertical fault section provided by an embodiment of the present application;

[0066] Figure 7 A schematic diagram of the structure of a co-seismic interference pattern unwrapping device provided in one embodiment of the present application;

[0067] Figure 8 A schematic diagram of the structure of a terminal device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0068] In the following description, specific details such as specific system structures, technologies, etc. are provided for the purpose of illustration rather than limitation, so as to provide a thorough understanding of the embodiments of the present application. However, it should be clear to those skilled in the art that the present application may also be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to prevent unnecessary details from obstructing the description of the present application.

[0069] It should be understood that when used in the present specification and the appended claims, the term "comprising" indicates the presence of described features, wholes, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, wholes, steps, operations, elements, components and / or combinations thereof.

[0070] It should also be understood that the term “and / or” used in the specification and appended claims refers to any and all possible combinations of one or more of the associated listed items, and includes these combinations.

[0071] As used in the specification and appended claims of this application, the term "if" can be interpreted as "when" or "uponce" or "in response to determining" or "in response to detecting", depending on the context. Similarly, the phrase "if it is determined" or "if [described condition or event] is detected" can be interpreted as meaning "uponce it is determined" or "in response to determining" or "uponce [described condition or event] is detected" or "in response to detecting [described condition or event]", depending on the context.

[0072] In addition, in the description of the present application specification and the appended claims, the terms "first", "second", "third", etc. are only used to distinguish the descriptions and cannot be understood as indicating or implying relative importance.

[0073] References to "one embodiment" or "some embodiments" etc. described in the specification of this application mean that one or more embodiments of the present application include specific features, structures or characteristics described in conjunction with the embodiment. Therefore, the statements "in one embodiment", "in some embodiments", "in some other embodiments", "in some other embodiments", etc. that appear in different places in this specification do not necessarily refer to the same embodiment, but mean "one or more but not all embodiments", unless otherwise specifically emphasized in other ways. The terms "including", "comprising", "having" and their variations all mean "including but not limited to", unless otherwise specifically emphasized in other ways.

[0074] In view of the low accuracy of the existing unwrapping of the interference graph of earthquake events, an embodiment of the present application provides a co-seismic interference graph unwrapping method, which determines multiple special center pixels from all pixels of the distance offset graph according to all offset standard deviations, and calculates the final offset of each special center pixel based on all distance offsets, so that pixels with large offset standard deviations can be processed separately, so that the accuracy of the pixel offset is improved; each significant pixel is connected with adjacent significant pixels to obtain multiple connected regions, and then the connected region with the largest area is used as the near field region, and each of the other connected regions is used as the far field region, based on all initial prior gradients, There are far-field regions and near-field regions, multiple final prior gradients are obtained, the near-field region and the far-field region are determined according to the significance of pixels, and the final prior gradient is obtained according to the near-field region and the far-field region, so that different processing is performed on the near-field region and the far-field region, so that the final prior gradient conforms to the different actual conditions of the near-field region and the far-field region; the objective function is constructed by using all final prior gradients without relying on the traditional phase continuity assumption, which can improve the accuracy of the objective function and reduce errors, so that the unwrapping phase of the interference pattern gradient obtained by maximizing the objective function is highly accurate; finally, the interference pattern unwrapping result is obtained by using the accurate unwrapping phase of the interference pattern gradient, which can effectively improve the accuracy of the unwrapping result.

[0075] Next, the co-seismic interference pattern unwrapping method provided in this application is exemplified.

[0076] like Figure 1 As shown, the co-seismic interference pattern unwrapping method provided by the present application comprises the following steps:

[0077] Step 11: Acquire a first SAR image and a second SAR image of the target area, and use the first SAR image and the second SAR image to obtain an interference map and a range offset map of the target area.

[0078] The time point corresponding to the first SAR image is the last satellite transit time before the earthquake event in the target area, and the time point corresponding to the second SAR image is the first satellite transit time after the earthquake event. For example, if the main shock of earthquake event A occurs at 5:00 on August 3, the first SAR image can be the SAR image of the target area at 9:00 on August 1, and the second SAR image can be the SAR image of the target area at 9:00 on August 13. The target area is the area where the co-seismic surface deformation study of the earthquake event is required, such as the area affected by earthquake event A.

[0079] In some embodiments of the present application, the first SAR image and the second SAR image of the target area can be obtained by accessing a website that discloses satellite image data, such as the Sentinel Science Data Center. The above-mentioned step of obtaining the interference map and the range offset map of the target area using the first SAR image and the second SAR image specifically includes:

[0080] In the first step, differential interferometry is performed on the first SAR image and the second SAR image to obtain an interference map of the target area.

[0081] It should be noted that the band of the above interference pattern is the C band, which is used to describe the surface deformation information of the target area after the earthquake event.

[0082] In the second step, the offset of the first SAR image and the second SAR image is tracked to obtain a range offset map of the target area.

[0083] For example, software such as GAMMA and SARscape may be used to perform differential interferometry and offset tracking on the first SAR image and the second SAR image to obtain an interference pattern and a range offset.

[0084] It is worth mentioning that by obtaining the interference pattern of the target area and the distance offset of each pixel, the surface deformation condition of the target area can be described, providing prior information of the surface deformation for the subsequent unwrapping step.

[0085] Step 12, performing median filtering on the range offset map to obtain the initial offset of each pixel in the range offset map, and calculating the offset standard deviation of each pixel based on all the initial offsets.

[0086] Specifically, the preset scale of the window is set according to the range size and noise level of the earthquake event, and then for each distance offset, the pixel corresponding to the distance offset is taken as the center, and the median of all distance offsets in the spatial filtering window is obtained. The median is used as the initial offset of the pixel, and then the standard deviation of all initial offsets in the spatial filtering window is calculated to obtain the offset standard deviation of the pixel.

[0087] Exemplarily, the above steps can be run using computer software such as MATLAB to obtain the initial offset and offset standard deviation of each pixel.

[0088] It is worth mentioning that the obtained distance offset often has data phenomena such as salt and pepper noise. By performing median filtering on the distance offset, data phenomena such as salt and pepper noise can be greatly reduced.

[0089] Step 13, according to all the offset standard deviations, a plurality of special central pixels are determined from all the pixels in the distance offset map, and a final offset of each special central pixel is calculated based on all the initial offsets.

[0090] The standard deviation of the offset of the special central pixel is greater than the standard deviation threshold.

[0091] In some embodiments of the present application, the step of calculating the final offset of each special central pixel based on all initial offsets is specifically to perform the following steps for each special central pixel:

[0092] In the first step, the maximum inter-class variance threshold segmentation method is used to divide other pixels in a window of a preset scale centered on the special central pixel into two categories.

[0093] In the second step, the difference between the initial offset means of the two types of pixels and the initial offset of the special central pixel is calculated respectively, and each pixel included in the pixel category with the smaller difference is taken as a homogeneous pixel corresponding to the special central pixel.

[0094] In the third step, the median of the initial offsets of all homogeneous pixels is calculated and used as the final offset of the special central pixel.

[0095] Exemplarily, the standard deviation threshold is 0.5, the standard deviation of the offset of the first pixel is 0.7, then the first pixel is the special center pixel, the preset scale is 3, then there are 9 pixels in the 3*3 window centered on the first pixel, and the maximum inter-class variance threshold segmentation method is used to divide the other 8 pixels into two categories, wherein the difference between the mean of these two categories of pixels and the initial offset of the special center pixel is 0.2 and 1.1 respectively, the pixels contained in the first category of pixels with the smaller difference are taken as homogeneous pixels, namely the second pixel, the fourth pixel, and the seventh pixel, and the initial offsets of the second pixel, the fourth pixel, and the seventh pixel are 1.04, 0.95, and 1.35 respectively, and the median is 1.04, then the final offset of the special center pixel is 1.04.

[0096] It is worth mentioning that for pixels whose offset standard deviation is greater than the standard deviation threshold, it can be considered that there are cross-surface rupture traces in the window centered on the pixel and the final offset difference is too large. By determining the homogeneous pixels and then taking the median, the influence of pixels that are too different from the central pixel can be avoided. The special filtering of the central pixel that only considers homogeneous pixels can improve the accuracy of the final offset of the central pixel.

[0097] Step 14, based on the final offsets of all special central pixels and the initial offsets of other pixels, obtain the saliency of each pixel in multiple initial prior gradients and the distance offset map.

[0098] The above initial prior gradient is used to describe the horizontal or vertical deformation gradient of the pixel, and the significance is used to describe the difference between the pixel and all other pixels.

[0099] In some embodiments of the present application, the above-mentioned method of obtaining the significance of each pixel in the multiple initial prior gradients and distance offset maps based on the final offsets of all special central pixels and the initial offsets of other pixels includes:

[0100] In the first step, the initial offset of each other pixel is used as the final offset of the other pixels.

[0101] Specifically, if the other pixels are pixels other than the special center pixel, for the pixel whose offset standard deviation is less than or equal to the standard deviation threshold in step 13, the initial offset of the pixel is used as the final offset of the pixel. That is, among all the pixels, for the pixel marked as the special center pixel in step 13, its final offset is the initial offset obtained after the processing of step 13, and for the pixel not marked as the special center pixel in step 13, its final offset is the initial offset obtained through step 12.

[0102] In the second step, for each pixel, the final offset of the pixel is subtracted from the final offset of the next pixel in the horizontal direction to obtain the horizontal initial prior gradient of the pixel, and the final offset of the pixel is subtracted from the final offset of the next pixel in the vertical direction to obtain the vertical initial prior gradient of the pixel.

[0103] It should be noted that if the pixel is the last pixel in the horizontal direction, no difference is made for the pixel in the horizontal direction, and if the pixel is the last pixel in the vertical direction, no difference is made for the pixel in the vertical direction. That is, the number of initial prior gradients is k = (m-1) × n + (n-1) × m, where m represents the number of rows of the distance offset map, and n represents the number of columns of the distance offset map.

[0104] In the third step, the saliency of each pixel is calculated based on the final offsets of all pixels.

[0105] Exemplarily, a dual-window based multi-scale feature extraction method may be used, and all final offsets may be used as inputs of the method to obtain the saliency of each pixel.

[0106] It is worth mentioning that the initial prior gradient obtained according to the final offset with high precision has high accuracy.

[0107] Step 15: According to all the saliencies, multiple salient pixels having saliency greater than a saliency threshold are determined from all the pixels in the distance-to-offset map, and regions formed by adjacent salient pixels are all regarded as connected regions.

[0108] Specifically, for each pixel, if the significance of the pixel is greater than the significance threshold, the pixel is marked as a significant pixel.

[0109] For example, if the significance threshold is 0.08, the significance of the first pixel is 0.04, and the significance of the second pixel is 0.09, the second pixel is marked as a significant pixel. If the second significant pixel, the third significant pixel, and the fourth significant pixel are adjacent, the region formed by the second significant pixel, the third significant pixel, and the fourth significant pixel is regarded as a connected region.

[0110] It is worth mentioning that for pixels whose significance is greater than the significance threshold, it can be considered that the difference between the pixel and all other pixels is high, and it is easier to pay attention to the pixel. The connected area obtained based on the significant pixel can reflect the area where significant deformation occurs in the target area.

[0111] Step 16, taking the connected region with the largest area as the near-field region, and each of the other connected regions as the far-field region, and obtaining multiple final prior gradients based on all initial prior gradients, all far-field regions, and near-field regions.

[0112] The above final prior gradient is used to describe the deformation of pixels in the near-field region or the far-field region.

[0113] In some embodiments of the present application, the above steps of obtaining multiple final prior gradients based on all initial prior gradients, all far-field regions, and near-field regions are specifically as follows:

[0114] In the first step, for each far-field region, multiple candidate prior gradients with absolute values ​​greater than π are determined from all initial prior gradients corresponding to all pixels in the far-field region.

[0115] In the second step, for each candidate prior gradient, if the candidate prior gradient is negative, -π is used as the final prior gradient corresponding to the candidate prior gradient; if the candidate prior gradient is positive, π is used as the final prior gradient corresponding to the candidate prior gradient.

[0116] Exemplarily, the first initial prior gradient in the first far-field region is -2π, then the initial prior gradient is used as a candidate prior gradient, and the candidate prior gradient is truncated to a value of -π to obtain the final prior gradient corresponding to the first initial prior gradient in the first far-field region.

[0117] It should be noted that after obtaining the final prior gradient corresponding to each candidate prior gradient, all final prior gradients in the far field area need to be smoothed by mean filtering.

[0118] In the third step, each initial prior gradient of the non-candidate prior gradient in each far-field region is taken as the final prior gradient, and each initial prior gradient in the near-field region is taken as the final prior gradient.

[0119] Exemplarily, the above steps can be performed using computer software such as Matlab to obtain multiple final prior gradients.

[0120] It should be noted that compared with the near-field area, the seismic deformation gradient in the far-field area is small in magnitude and smooth and continuous, and should meet the characteristics of the phase continuity assumption. If the double earthquake surface deformation of the target area needs to be studied, the connected area with the largest area and the connected area with the second largest area are both regarded as the near-field area, and the other connected areas are all regarded as the far-field area.

[0121] It is worth mentioning that by determining the near-field area and the far-field area, and truncating and smoothing the initial prior gradient of the far-field area, the gradient of the far-field area can be made to conform to the actual situation that the magnitude of the surface deformation gradient in the far-field area is small, thereby improving the accuracy of all final prior gradients.

[0122] Step 17, constructing an objective function using all final prior gradients, and maximizing the objective function to obtain the unwrapped phase of each interference pattern gradient.

[0123] The above-mentioned objective function is used to describe the probability of the unwrapped phase of all interference pattern gradients. The interference pattern gradient is a vector obtained by subtracting the wrapped phase of each pixel in the interference pattern from that of adjacent pixels in the horizontal or vertical direction.

[0124] In some embodiments of the present application, the above steps of constructing the objective function using all final a priori gradients and maximizing the objective function to obtain the unwrapped phase of each interference pattern gradient are specifically as follows:

[0125] In the first step, the disentangled probability density function of each interferogram gradient is constructed using all final prior gradients.

[0126] Specifically, the above disentanglement probability density function is:

[0127]

[0128] in, represents the unwrapping probability density of the j-th interference pattern gradient, represents the unwrapped phase to be solved for the jth interferogram gradient, j = 1, 2, ..., k, k represents the total number of interferogram gradients, σ ​​represents the coherence parameter, G RO represents all final prior gradients, G INF represents the winding gradient corrected mean, T coh represents the coherence threshold, and coherence represents the coherence of the j-th interference pattern gradient.

[0129] It should be noted that the above winding gradient correction mean is obtained based on the phase continuity assumption.

[0130] In the second step, for each interference pattern gradient, based on the unwrapping probability density function corresponding to the interference pattern gradient, the conditional probability density function corresponding to the interference pattern gradient is obtained.

[0131] Specifically, the above conditional probability density function is:

[0132]

[0133] in, represents the conditional probability density of the j-th interference pattern gradient, represents the winding phase of the j-th interference pattern gradient, and a and b represent arbitrary integers.

[0134] The third step is to construct the objective function based on all conditional probability density functions.

[0135] Specifically, the above objective function is:

[0136]

[0137] Wherein, F(ΔΦ|ΔΨ) represents the objective function value, ΔΦ represents the unwrapped phase of all interferogram gradients, and ΔΨ represents the wrapped phase of all interferogram gradients.

[0138] The fourth step is to maximize the objective function and obtain the unwrapped phase of the interference pattern gradient.

[0139] Exemplarily, the objective function can be maximized by using a dynamic cost cycle cancellation algorithm (DCC) to obtain the optimal value of all interferogram gradient unwrapping phases ΔΦ. The objective function can be constructed using computer software such as matlab, and the dynamic cost cycle cancellation algorithm in SNAPHU software can be used to maximize the objective function to obtain the unwrapping phase.

[0140] It is worth mentioning that using the final prior gradient to construct the objective function without relying on the traditional phase continuity assumption can improve the accuracy of the objective function and reduce the error, so that the unwrapped phase obtained by maximizing the objective function is highly accurate.

[0141] Step 18: perform path integration on all unwrapped phases to obtain an interference pattern unwrapping result of the target area.

[0142] The above interferogram unwrapping results are used to describe the seismic deformation in the target area.

[0143] Specifically, path integration may be performed on the unwrapping phases of all interference pattern gradients to obtain a plurality of unwrapping phase values ​​of the target region, and all the unwrapping phase values ​​are used as unwrapping results of the target region.

[0144] Exemplarily, software such as Matlab and SNAPHU can be used to perform path integration on the unwrapped phase to obtain the unwrapped result.

[0145] It is worth mentioning that the unwrapped phase of the interference pattern gradient solved in step 17 already has a rotation-free characteristic, and a simple arbitrary path integral can be performed to obtain the unwrapped result.

[0146] The co-seismic interference pattern unwrapping method of the present application is exemplarily described below with reference to a specific example.

[0147] A double-fault moment dislocation model is constructed based on the fault geometry of the 2023 Türkiye large double earthquake event, and its InSAR observed deformation field is simulated as follows: Figure 2 The unwrapped truth value is wrapped and decoherent noise is added based on the true coherence of the earthquake event to obtain Figure 2 (b) shows the simulated C-band co-seismic InSAR interferogram. Perlin noise is used to simulate the speckle noise of the range offset, and the following is obtained: Figure 2 (c) shows a simulated range offset with a noise standard deviation (NoiseSTD) of 12.7 cm. The simulated C-band coseismic InSAR interferograms were unwrapped using MCF, SNAPHU, and the coseismic interferogram unwrapping method of this application. To ensure a fair comparison, the three methods use the same surface rupture trace mask and all use coherence as auxiliary data.

[0148] The unwrapping results obtained by using MCF are as follows Figure 2 As shown in (d), the untangling result obtained by using SNAPHU is as follows Figure 2 (e) shows that the unwrapping result obtained by using the co-seismic interferogram unwrapping method of the present application is as follows Figure 2 As shown in (f), the disentanglement residual between the disentanglement result corresponding to MCF and the disentanglement truth value is as follows: Figure 2 As shown in (h), the root mean square error (RMSE) is 7.6 cm, and the unwrapping residual between the unwrapping result corresponding to SNAPHU and the unwrapping true value is as follows: Figure 2 As shown in (i), the RMSE is 5.3 cm. The unwrapping residual between the unwrapping result and the unwrapping true value corresponding to the coseismic interferogram unwrapping method of the present application is as follows: Figure 2 As shown in (j), the RMSE is 1.9 cm.

[0149] It can be seen that the MCF and SNAPHU methods significantly underestimate the near-field deformation close to the fault rupture. The magnitude and spatial distribution of the unwrapping results corresponding to the coseismic interferogram unwrapping method of this application are relatively consistent with the simulated unwrapping true value. MCF not only has large residuals in the near-field area, but also has certain residuals in the far-field area of ​​the upper half of the interferogram, indicating that the unwrapping error has spread over a large range, and the residual RMSE has reached 7.6cm; the residuals of SNAPHU are concentrated in the near-field area, with an RMSE of 5.3cm; and the unwrapping error of the coseismic interferogram unwrapping method of this application is greatly reduced, with an RMSE of only 1.9cm. The small amount of residuals mainly comes from the noise of the distance offset and the incoherence noise of the interferogram itself.

[0150] In order to further explore the influence of the noise level of the range offset on the accuracy of the method proposed in this invention, 60 sets of range offsets with different noise levels (the standard deviation of the noise is between 5 and 25 cm) are simulated, and the unwrapping accuracy of the co-seismic interferogram unwrapping method of this application is calculated when different range offsets are used as prior information. The statistical results of the unwrapping residuals are shown in Figure 2. Figure 3 As shown, Figure 3 The horizontal axis is the noise standard deviation (Noise STD) of the range offset, and the vertical axis is the unwrapping residual RMSE, both in centimeters (cm). The points are simulated samples, and the straight lines are the results of linear fitting. The unwrapping accuracy of the co-seismic interferogram unwrapping method of the present application shows an obvious positive correlation with the noise level of the range offset. When the standard deviation of the range offset is in the normal range of practical application (i.e., 5-25cm), the overall error obtained by this method is relatively small, and the RMSE is less than 3.5cm. The standard deviation of the actual range offset is often between 10 and 15cm. Within this range, the average RMSE of this method is 1.9cm, which is 76% and 64% higher than that of MCF and SNAPHU, respectively.

[0151] Next, the co-seismic interference pattern unwrapping method of the present application is exemplified in conjunction with a specific experiment.

[0152] For the 2016 Mw7.0 Kumamoto earthquake, according to the co-seismic interferogram unwrapping method of this application, the two SAR images of 20160408 and 20160420 are first registered, and the original range offset obtained by the offset tracking technology is as follows: Figure 4 As shown in (a), the window size is set to 2.5 km and the standard deviation threshold is 0.3. The final offset obtained by ordinary median filtering is as follows: Figure 4 As shown in (b), the final offset obtained by filtering using the co-seismic interferogram unwrapping method of the present application is as follows: Figure 4 (c) (i.e., homogeneous pixel median filtering in the figure), the difference between the two filtering methods is as follows Figure 4 (d) as shown.

[0153] It can be seen that the homogeneous pixel median filtering method in the co-seismic interferogram unwrapping method of the present application effectively avoids the problem of excessive smoothing of large-scale deformations across the surface rupture area.

[0154] Figure 5 The unwrapping results of MCF, SNAPHU and the method of the present application on the coseismic interferogram are compared, as well as the comparison with several Global Navigation Satellite System (GNSS) stations and local LiDAR (Light Detection and Ranging) data in the region. Figure 5 The middle triangles represent GNSS stations, the thick solid lines represent the surface rupture traces, and the thin dashed lines represent the near-field contours determined by this method. Figure 5 (a) is the InSAR coseismic interferogram, rad is radians, and AA'-DD' represents four sections perpendicular to the fault trace, which are used for subsequent detailed comparative analysis. Figure 5 (d) is the prior information of the distance offset. The unwrapping result of MCF is as follows: Figure 5 As shown in (b), the untangling result of SNAPHU is as follows Figure 5 (c) shows, the untangling result of the method of the present application is as follows Figure 5 The scatter plot between the unwrapping results of the method of the present application and the slant range (LOS, Line Of Sight) deformation of 23 GNSS stations in the region is shown in (e). Figure 5(f) shows that the horizontal axis is the unwrapping value of the method of the present application (i.e. the unwrapping result mentioned above), in centimeters (cm), and the vertical axis is the LOS deformation data of the GNSS station, in cm, and the mean absolute deviation (MAE, Masked Autoencoders) is only 2.3cm, which is 34% more accurate than the traditional method. In order to highlight the advantages of the present application in the near field area, the unwrapping results of different methods are compared with the mean absolute deviation (N-MAE) of the three GNSS stations in the near field area, and the N-MAE of MCF, SNAPHU and the method of the present application are 24.0cm, 16.8cm and 4.8cm respectively. Figure 5 (a) Figure 5 (b) Figure 5 (c) Figure 5 (d) Figure 5 In (e), the horizontal axis represents longitude, and the vertical axis represents latitude.

[0155] The LOS deformation along the four sections AA'-DD' is as follows: Figure 6 As shown, Figure 6 The horizontal axis represents the distance along the profile in kilometers (km), and the vertical axis represents the LOS deformation in meters (m). Figure 6 (a) is the LOS deformation comparison of the AA' section, where the dot symbol represents the range offset curve, the five-pointed star represents the LOS deformation curve of the method of the present application, the square represents the LOS deformation curve of MCF, the triangle represents GNSS data, the diamond represents the LOS deformation curve of SNAPHU, and the plus sign represents LiDAR data. Figure 6 (b) is the LOS deformation comparison of the BB' section. Figure 6 (c) is the LOS deformation comparison of CC' section, Figure 6 (d) is the LOS deformation comparison of the DD' section. Figure 6 It is further shown that the traditional MCF and SNAPHU methods may have a deviation of more than 1 meter in the near-field area, which is invalid unwrapping; while the unwrapping results of the applied method in the near-field area are closer to the external verification data such as range offset, GNSS and LiDAR, and are less affected by the noise of the range offset.

[0156] It can be seen that the co-seismic interferogram unwrapping method of the present application can effectively improve the accuracy of the unwrapping results and reduce errors.

[0157] The co-seismic interference pattern unwrapping device provided in the present application is exemplarily described below.

[0158] like Figure 7 As shown, the embodiment of the present application provides a co-seismic interference pattern unwrapping device, and the co-seismic interference pattern unwrapping device 700 includes:

[0159] The acquisition module 701 acquires a first SAR image and a second SAR image of the target area, and uses the first SAR image and the second SAR image to obtain an interference diagram and a range offset diagram of the target area; the time point corresponding to the first SAR image is the time when the satellite last passes over the target area before the earthquake event occurs, and the time point corresponding to the second SAR image is the time when the satellite first passes over the target area after the earthquake event occurs;

[0160] The calculation module 702 performs median filtering on the distance offset map to obtain an initial offset of each pixel in the distance offset map, and calculates the offset standard deviation of each pixel based on all the initial offsets;

[0161] The first determination module 703 determines a plurality of special central pixels from the distance offset map according to all offset standard deviations, and calculates a final offset of each special central pixel based on all initial offsets; the offset standard deviation of the special central pixel is greater than a standard deviation threshold;

[0162] The saliency acquisition module 704 acquires the saliency of each pixel in the multiple initial prior gradients and the distance offset map based on the final offsets of all special center pixels and the initial offsets of other pixels; the initial prior gradient is used to describe the horizontal or vertical deformation gradient of the pixel, and the saliency is used to describe the difference between the pixel and all other pixels;

[0163] The second determination module 705 determines, according to all the saliencies, a plurality of salient pixels having saliencies greater than a saliency threshold from all the pixels in the distance-to-offset map, and regards the regions formed by the plurality of adjacent salient pixels as connected regions;

[0164] The final prior gradient acquisition module 706 takes the connected area with the largest area as the near field area, and each other connected area as the far field area, and acquires multiple final prior gradients based on all initial prior gradients, all far field areas and near field areas;

[0165] A construction module 707 is used to construct an objective function using all final prior gradients, and to maximize the objective function to obtain an unwrapped phase of each interferogram gradient; the objective function is used to describe the probability of the unwrapped phase of all interferogram gradients, and the interferogram gradient is a vector obtained by subtracting the wrapped phase of each pixel in the interferogram from that of adjacent pixels in the horizontal or vertical direction;

[0166] The path integration module 708 performs path integration on all the unwrapped phases to obtain an interference pattern unwrapping result of the target area; the interference pattern unwrapping result is used to describe the seismic deformation of the target area.

[0167] It should be noted that the information interaction, execution process, etc. between the above-mentioned devices / units are based on the same concept as the method embodiment of the present application. Their specific functions and technical effects can be found in the method embodiment part and will not be repeated here.

[0168] The technicians in the relevant field can clearly understand that for the convenience and simplicity of description, only the division of the above-mentioned functional units and modules is used as an example for illustration. In practical applications, the above-mentioned function allocation can be completed by different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiment can be integrated in a processing unit, or each unit can exist physically separately, or two or more units can be integrated in one unit. The above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units. In addition, the specific names of the functional units and modules are only for the convenience of distinguishing each other, and are not used to limit the scope of protection of this application. The specific working process of the units and modules in the above-mentioned system can refer to the corresponding process in the aforementioned method embodiment, which will not be repeated here.

[0169] like Figure 8 As shown, an embodiment of the present application provides a terminal device. The terminal device D10 of this embodiment includes: at least one processor D100 ( Figure 8 Only one processor is shown in the figure), a memory D101, and a computer program D102 stored in the memory D101 and executable on the at least one processor D100, wherein the processor D100 implements the steps of any of the above-mentioned method embodiments when executing the computer program D102.

[0170] Specifically, when the processor D100 executes the computer program D102, a plurality of special center pixels are determined from all pixels in the distance offset map according to all offset standard deviations, and the final offset of each special center pixel is calculated based on all distance offsets, so that pixels with large offset standard deviations can be processed separately, so that the accuracy of pixel offsets is improved; each significant pixel is connected to adjacent significant pixels to obtain a plurality of connected regions, and then the connected region with the largest area is used as the near-field region, and each of the other connected regions is used as the far-field region, and based on all initial prior gradients, all far-field regions and near-field regions, a plurality of connected regions are obtained. Multiple final prior gradients are used to determine the near-field region and the far-field region according to the significance of pixels, and the final prior gradient is obtained according to the near-field region and the far-field region, so as to realize different processing of the near-field region and the far-field region, so that the final prior gradient conforms to the different actual conditions of the near-field region and the far-field region; the objective function is constructed by using all final prior gradients without relying on the traditional phase continuity assumption, which can improve the accuracy of the objective function and reduce the error, so that the unwrapping phase of the interference pattern gradient obtained by maximizing the objective function is highly accurate; finally, the unwrapping result is obtained by using the accurate unwrapping phase of the interference pattern gradient, which can effectively improve the accuracy of the unwrapping result.

[0171] The processor D100 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc.

[0172] In some embodiments, the memory D101 may be an internal storage unit of the terminal device D10, such as a hard disk or memory of the terminal device D10. In other embodiments, the memory D101 may also be an external storage device of the terminal device D10, such as a plug-in hard disk, a smart memory card (SMC, SmartMedia Card), a secure digital (SD, Secure Digital) card, a flash card (Flash Card), etc. equipped on the terminal device D10. Further, the memory D101 may also include both an internal storage unit of the terminal device D10 and an external storage device. The memory D101 is used to store an operating system, an application program, a boot loader (BootLoader), data and other programs, such as the program code of the computer program, etc. The memory D101 may also be used to temporarily store data that has been output or is to be output.

[0173] An embodiment of the present application further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the steps in the above-mentioned method embodiments can be implemented.

[0174] An embodiment of the present application provides a computer program product. When the computer program product runs on a terminal device, the terminal device can implement the steps in the above-mentioned method embodiments when executing the computer program product.

[0175] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the present application implements all or part of the processes in the above-mentioned embodiment method, which can be completed by instructing the relevant hardware through a computer program. The computer program can be stored in a computer-readable storage medium. When the computer program is executed by the processor, the steps of the above-mentioned various method embodiments can be implemented. Among them, the computer program includes computer program code, and the computer program code can be in source code form, object code form, executable file or some intermediate form. The computer-readable medium may at least include: any entity or device that can carry the computer program code to the co-seismic interferogram unwrapping method device / terminal device, recording medium, computer memory, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), electric carrier signal, telecommunication signal and software distribution medium. For example, a USB flash drive, a mobile hard disk, a magnetic disk or an optical disk. In some jurisdictions, according to legislation and patent practice, computer-readable media cannot be electric carrier signals and telecommunication signals.

[0176] In the above embodiments, the description of each embodiment has its own emphasis. For parts that are not described or recorded in detail in a certain embodiment, reference can be made to the relevant descriptions of other embodiments.

[0177] Those of ordinary skill in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of this application.

[0178] The above is a preferred embodiment of the present application. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles described in the present application. These improvements and modifications should also be regarded as the scope of protection of the present application.

Claims

1. A co-seismic interference pattern unwrapping method, characterized in that: include: Acquire a first SAR image and a second SAR image of a target area, and use the first SAR image and the second SAR image to obtain an interference map and a range offset map of the target area; the time point corresponding to the first SAR image is the time when the satellite last passes over the target area before the earthquake event occurs, and the time point corresponding to the second SAR image is the time when the satellite first passes over the target area after the earthquake event occurs; Performing median filtering on the distance offset map to obtain an initial offset of each pixel in the distance offset map, and calculating an offset standard deviation of each pixel based on all the initial offsets; Determining a plurality of special center pixels from all pixels of the distance offset map according to all offset standard deviations, and calculating a final offset of each special center pixel based on all initial offsets; The standard deviation of the offset of the special central pixel is greater than the standard deviation threshold; Based on the final offsets of all special central pixels and the initial offsets of other pixels, a plurality of initial prior gradients and the significance of each pixel in the distance offset map are obtained; the initial prior gradient is used to describe the horizontal or vertical deformation gradient of the pixel, and the significance is used to describe the degree of difference between the pixel and all other pixels; According to all the saliencies, a plurality of salient pixels having saliencies greater than a saliency threshold are determined from all the pixels in the distance-to-offset map, and regions formed by adjacent salient pixels are all considered as connected regions; The connected region with the largest area is taken as the near-field region, and each of the other connected regions is taken as the far-field region, and multiple final prior gradients are obtained based on all initial prior gradients, all far-field regions and the near-field region; the final prior gradients are used to describe the deformation of pixels in the near-field region or the far-field region; An objective function is constructed using all final prior gradients, and the objective function is maximized to obtain an unwrapped phase of each interferogram gradient; the objective function is used to describe the probability of the unwrapped phase of all interferogram gradients, and the interferogram gradient is a vector obtained by subtracting the entangled phase of each pixel in the interferogram from that of adjacent pixels in the horizontal or vertical direction; Path integration is performed on all unwrapped phases to obtain an interference pattern unwrapping result of the target area; the interference pattern unwrapping result is used to describe the seismic deformation of the target area.

2. The co-seismic interference pattern unwrapping method according to claim 1, characterized in that: The step of obtaining the interference map and the range offset map of the target area by using the first SAR image and the second SAR image includes: Performing differential interferometry on the first SAR image and the second SAR image to obtain an interference map of the target area; The first SAR image and the second SAR image are subjected to offset tracking to obtain a range offset map of the target area.

3. The co-seismic interference pattern unwrapping method according to claim 1, characterized in that: The final offset of each special center pixel is calculated based on all initial offsets, including: For each of the special central pixels, the following steps are performed: Using the maximum inter-class variance threshold segmentation method, other pixels in a window of a preset scale centered on the special central pixel are divided into two categories; The difference between the average of the initial offsets of the two types of pixels and the initial offset of the special central pixel is calculated respectively, and each pixel included in the pixel type with the smaller difference is used as a homogeneous pixel corresponding to the special central pixel; The median of the initial offsets of all homogeneous pixels is calculated and used as the final offset of the special central pixel.

4. The co-seismic interference pattern unwrapping method according to claim 1, characterized in that: The method of obtaining a plurality of initial a priori gradients and the significance of each pixel based on the final offset of all special central pixels and the initial offsets of other pixels comprises: Using the initial offset of each of the other pixels as the final offset of the other pixels; For each pixel, respectively, subtract the final offset of the pixel from the final offset of the next pixel in the horizontal direction to obtain the initial prior gradient of the pixel in the horizontal direction, and subtract the final offset of the pixel from the final offset of the next pixel in the vertical direction to obtain the initial prior gradient of the pixel in the vertical direction; Based on the final offsets of all pixels, the saliency of each of said pixels is calculated.

5. The co-seismic interference pattern unwrapping method according to claim 1, characterized in that: The step of obtaining a plurality of final prior gradients based on all initial prior gradients, all far-field regions and the near-field region comprises: For each of the far-field regions, determine a plurality of candidate prior gradients whose absolute values ​​are greater than π from all initial prior gradients corresponding to all pixels in the far-field region; For each candidate prior gradient, if the candidate prior gradient is a negative number, -π is used as the final prior gradient corresponding to the candidate prior gradient; if the candidate prior gradient is a positive number, π is used as the final prior gradient corresponding to the candidate prior gradient; Each initial prior gradient of the non-candidate prior gradient in each of the far-field regions is taken as the final prior gradient, and each initial prior gradient in the near-field region is taken as the final prior gradient.

6. The co-seismic interference pattern unwrapping method according to claim 1, characterized in that: The objective function is constructed using all final prior gradients, including: Constructing the unwrapping probability density function of each interferogram gradient using all final prior gradients; For each of the interference pattern gradients, based on the unwrapping probability density function corresponding to the interference pattern gradient, obtaining a conditional probability density function corresponding to the interference pattern gradient; Construct the objective function based on all conditional probability density functions.

7. The co-seismic interference pattern unwrapping method according to claim 6, characterized in that: The unwrapping probability density function is: in, represents the unwrapping probability density of the j-th interference pattern gradient, represents the unwrapped phase to be solved for the jth interferogram gradient, j = 1, 2, ..., k, k represents the total number of interferogram gradients, σ ​​represents the coherence parameter, G RO represents all final prior gradients, G INF represents the winding gradient corrected mean, T coh represents the coherence threshold, and coherence represents the coherence of the j-th interference pattern gradient; The conditional probability density function is: in, represents the conditional probability density of the j-th interference pattern gradient, represents the winding phase of the j-th interference pattern gradient, and a and b both represent arbitrary integers.

8. The co-seismic interference pattern unwrapping method according to claim 7, characterized in that: The objective function is: Wherein, F(ΔΦ|ΔΨ) represents the objective function value, ΔΦ represents the unwrapped phase of all interferogram gradients, and ΔΨ represents the wrapped phase of all interferogram gradients.

9. A terminal device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the co-seismic interference pattern unwrapping method according to any one of claims 1 to 8 is implemented.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the co-seismic interference pattern unwrapping method according to any one of claims 1 to 8 is implemented.

Images