Landslide deformation estimation method and device, storage medium and computer equipment

Abstract

The present disclosure relates to the technical field of ground deformation detection, and provides a landslide deformation estimation method and device, a storage medium and computer equipment, which comprises: acquiring multiple differential interferograms of a target study area and a coherence coefficient map corresponding to each differential interferogram; for each differential interferogram, determining a deformation gradient field of the differential interferogram in multiple gradient extraction directions, and determining a weight field corresponding to each deformation gradient field based on the coherence coefficient map corresponding to the differential interferogram; based on the deformation gradient fields of all differential interferograms in the same gradient extraction direction and the corresponding weight fields, determining an average deformation gradient field of the target study area in the gradient extraction direction; and using Poisson equation to reconstruct a deformation rate field according to the average deformation gradient fields in each gradient extraction direction, so as to obtain a deformation rate field of the target study area and determine a deformation estimation result of the target study area. The present embodiment realizes deformation rate field reconstruction and improves estimation reliability and spatial continuity.

Description

Technical Field

[0001] This disclosure relates to the field of surface deformation detection technology, and more specifically, to a landslide deformation estimation method, apparatus, storage medium, and computer equipment. Background Technology

[0002] Interferometric Synthetic Aperture Radar (InSAR) technology possesses high-precision, all-weather, and wide-area Earth observation capabilities, enabling the acquisition of surface deformation information at the centimeter or even millimeter level. Temporal InSAR technology further suppresses the effects of atmospheric delay and spatiotemporal incoherence, and has been widely applied in deformation monitoring fields such as land subsidence, landslides, earthquakes, volcanic activity, and glacier movement, becoming an important technical means for geological disaster research and investigation.

[0003] In related technologies, phase gradient processing methods mainly avoid the phase unwrapping step by calculating the phase difference between adjacent pixels in the interferogram. However, the stacking results of this type of method only contain dimensionless values ​​and lack physical quantities corresponding to the surface deformation rate, thus failing to characterize the complete deformation field inside the landslide body. For uniform sliding or slowly creeping areas inside the landslide body, this type of method is difficult to extract effective signals. At the same time, interferograms often suffer from incoherence problems in areas with dense vegetation cover and fragmented terrain. The monitoring point density is low, and related processing methods rely on phase unwrapping. However, unwrapping errors are prone to propagate along the path, affecting all subsequent inversion results. Furthermore, large-scale matrix operations lead to long processing cycles, which contradicts the rapid response requirements of landslide monitoring. Summary of the Invention

[0004] This disclosure provides at least one landslide deformation estimation method, apparatus, storage medium, and computer equipment. By using multi-directional weighted averaging and Poisson integral reconstruction, the deformation rate field is reconstructed, improving the estimation reliability and spatial continuity.

[0005] This disclosure provides a landslide deformation estimation method, including: Obtain multiple differential interferograms of the target study area and a coherence coefficient map corresponding to each differential interferogram; For each differential interferogram, the deformation gradient field of the differential interferogram in multiple preset gradient extraction directions is determined, and the weight field corresponding to each deformation gradient field is determined based on the coherence coefficient map corresponding to the differential interferogram. Based on the deformation gradient field and corresponding weight field of all the differential interferograms in the same gradient extraction direction, the average deformation gradient field of the target study area in the gradient extraction direction is determined. Using the Poisson equation, the deformation rate field is reconstructed based on the average deformation gradient field in each of the gradient extraction directions to obtain the deformation rate field of the target study area; and the deformation estimation result of the target study area is determined based on the deformation rate field of the target study area.

[0006] This disclosure provides a landslide deformation estimation device, comprising: The data acquisition module is used to acquire multiple differential interferograms of the target study area and a coherence coefficient map corresponding to each differential interferogram; The gradient and weight field determination module is used to determine the deformation gradient field of each differential interferogram in multiple preset gradient extraction directions, and to determine the weight field corresponding to each deformation gradient field based on the coherence coefficient map corresponding to the differential interferogram. The average gradient field determination module is used to determine the average deformation gradient field of the target study area in the gradient extraction direction based on the deformation gradient field and the corresponding weight field of all the differential interferograms in the same gradient extraction direction. The deformation estimation module is used to reconstruct the deformation rate field based on the average deformation gradient field in each of the gradient extraction directions using the Poisson equation, thereby obtaining the deformation rate field of the target study area; and to determine the deformation estimation result of the target study area based on the deformation rate field of the target study area.

[0007] This disclosure provides a computer device including a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the computer device is running, the processor communicates with the memory via the bus. When the machine-readable instructions are executed by the processor, they perform the landslide deformation estimation method as described in any of the above possible embodiments.

[0008] This disclosure provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the landslide deformation estimation method as described in any of the possible embodiments above.

[0009] The landslide deformation estimation method, apparatus, storage medium, and computer equipment provided in this disclosure construct deformation gradient fields and their corresponding coherence coefficient weight fields in multiple gradient extraction directions, and perform weighted averaging based on all differential interferograms in the same direction. The contribution of high coherence regions to the average deformation gradient field is enhanced, and the noise influence of low coherence regions is effectively suppressed, thereby obtaining a more stable and reliable average deformation gradient field in each direction. Furthermore, the average deformation gradient field in multiple directions is reconstructed by integrating the Poisson equation, restoring the gradient domain information to a globally consistent deformation rate field. This avoids the error accumulation problem caused by pixel-by-pixel unwrapping or single-direction integration, and improves the spatial continuity and numerical accuracy of the deformation rate field. Based on the reconstructed deformation rate field, the deformation estimation results are determined, which can characterize the spatial distribution of surface deformation in the target study area and provide a quantitative basis for landslide disaster identification and monitoring.

[0010] Thus, this disclosure achieves an effective characterization of the spatial distribution of surface deformation in the target study area by performing weighted averaging based on the coherence coefficient weight field in the multi-gradient extraction direction and reconstructing the deformation rate field using the Poisson equation, thereby improving the reliability and spatial continuity of the deformation estimation results.

[0011] To make the above-mentioned objects, features and advantages of this disclosure more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0012] To more clearly illustrate the technical solutions of the embodiments of this disclosure, the accompanying drawings referenced in the embodiments will be briefly described below. These drawings are incorporated in and constitute a part of this specification. They illustrate embodiments conforming to this disclosure and, together with the specification, serve to explain the technical solutions of this disclosure. It should be understood that the following drawings only show some embodiments of this disclosure and should not be considered as limiting the scope. Those skilled in the art can obtain other related drawings based on these drawings without creative effort.

[0013] Figure 1 A flowchart of a landslide deformation estimation method provided by an embodiment of this disclosure is shown; Figure 2 A flowchart of a deformation gradient field determination method provided by an embodiment of this disclosure is shown; Figure 3 A flowchart of a weight field determination method provided by an embodiment of this disclosure is shown; Figure 4 A flowchart of a method for determining an average deformation gradient field provided by an embodiment of this disclosure is shown; Figure 5A flowchart of a deformation rate field reconstruction method provided by an embodiment of this disclosure is shown; Figure 6 This diagram illustrates a comparison of deformation estimation results using different methods provided in an embodiment of this disclosure. Figure 6 (a) is a schematic diagram of the landslide area in the target study area provided in an embodiment of this disclosure. Figure 6 (b) is a schematic diagram of the deformation result of the conventional phase gradient rate method provided in the embodiments of this disclosure. Figure 6 (c) is a schematic diagram of the deformation result of the improved phase gradient stacking method provided in the embodiments of this disclosure. Figure 6 (d) is a schematic diagram of the deformation result of the method of this disclosure provided in the embodiments of this disclosure; Figure 7 A schematic diagram of the structure of a landslide deformation estimation device provided in an embodiment of this disclosure is shown; Figure 8 A schematic diagram of the structure of a computer device provided in an embodiment of this disclosure is shown. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of the embodiments of this disclosure clearer, the technical solutions of the embodiments of this disclosure will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this disclosure, and not all of them. The components of the embodiments of this disclosure described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of this disclosure provided in the accompanying drawings is not intended to limit the scope of the claimed disclosure, but merely represents selected embodiments of this disclosure. All other embodiments obtained by those skilled in the art based on the embodiments of this disclosure without inventive effort are within the scope of protection of this disclosure.

[0015] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0016] In this document, the term "and / or" merely describes a relationship, indicating that three relationships can exist. For example, A and / or B can represent three cases: A alone, A and B simultaneously, and B alone. Furthermore, the term "at least one" in this document means any combination of at least two of any one or more elements. For example, including at least one of A, B, and C can mean including any one or more elements selected from the set consisting of A, B, and C.

[0017] Interferometric Synthetic Aperture Radar (InSAR) is a space-based Earth observation technology developed in recent years. It can acquire surface deformation information with centimeter-level or even millimeter-level accuracy, offering advantages such as high precision, all-weather operation, and wide-area monitoring. Temporal InSAR technology can effectively suppress some of the effects of atmospheric delay and spatiotemporal incoherence, and is widely used in various fields of surface deformation monitoring, including land subsidence, landslides, earthquakes, volcanic activity, and glacial movement, becoming an important tool for geological disaster research.

[0018] However, these methods heavily rely on phase unwrapping, which is their most vulnerable aspect. Interferograms record phase values ​​wrapped between -π and π, and unwrapping is essential to recover the true deformation. However, in densely vegetated, fragmented mountainous areas, the interferometric phase signals are highly fragmented, making the unwrapping algorithm prone to erroneous transitions. These errors propagate along the path to all subsequent deformation inversions, leading to unreliable monitoring results. Furthermore, factors such as vegetation growth and soil moisture changes can cause radar echo signals to lose coherence, resulting in large unusable areas in the interferogram and extremely sparse monitoring points, making it impossible to obtain sufficient deformation information. In addition, traditional time-series InSAR methods require the construction of complex spatiotemporal networks and large-scale matrix operations, resulting in enormous computational loads and long processing cycles, making it difficult to meet the timeliness requirements of wide-area rapid surveys.

[0019] Research has revealed that, in response to the aforementioned technical challenges, Phase Gradient Rate (PGR) was an early simplified processing strategy. Its basic idea is to calculate the phase difference, or phase gradient, between adjacent pixels in the interferogram instead of solving for the absolute deformation value of each pixel. Since real deformation often manifests as a continuous phase change in space, while atmospheric interference and noise have random or large-scale slowly varying characteristics, calculating the gradient can effectively highlight local deformation signals.

[0020] However, traditional PGR phase gradient stacking results do not have clearly corresponding physical quantities. The gradient field only reflects the boundary characteristics of the deformation zone to a certain extent and does not correspond to the complete deformation field. It cannot identify uniform or slowly deformed regions inside the landslide body, and the boundary range of the gradient stacking field becomes blurred as the gradient step size increases.

[0021] Based on the above research, this disclosure provides a landslide deformation estimation method, apparatus, storage medium, and computer equipment, including: constructing deformation gradient fields and their corresponding coherence coefficient weight fields in multiple gradient extraction directions, and performing weighted averaging based on all differential interferograms in the same direction. The contribution of high coherence regions to the average deformation gradient field is enhanced, and the noise influence of low coherence regions is effectively suppressed, thereby obtaining a more stable and reliable average deformation gradient field in each direction; further, the average deformation gradient field in multiple directions is reconstructed by integrating the Poisson equation, restoring the gradient domain information to a globally consistent deformation rate field, avoiding the error accumulation problem caused by pixel-by-pixel unwrapping or single-direction integration, and improving the spatial continuity and numerical accuracy of the deformation rate field; determining the deformation estimation result based on the reconstructed deformation rate field, which can characterize the spatial distribution of surface deformation in the target study area and provide a quantitative basis for landslide disaster identification and monitoring.

[0022] In this embodiment of the disclosure, by performing weighted averaging based on the coherence coefficient weight field in the multi-gradient extraction direction and reconstructing the deformation rate field using the Poisson equation, the spatial distribution of surface deformation in the target study area is effectively characterized, thereby improving the reliability and spatial continuity of the deformation estimation results.

[0023] To facilitate understanding of this embodiment, the executing entity of the landslide deformation estimation method provided in this disclosure embodiment will first be described in detail. The executing entity of the landslide deformation estimation method provided in this disclosure embodiment is a computer device. This computer device can be a terminal device or a server. The terminal device can also be a mobile device, user terminal, terminal, handheld device, computing device, vehicle-mounted device, wearable device, etc. The server can be an independent physical server, a server cluster or distributed system composed of multiple physical servers, or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud storage, big data, and artificial intelligence platforms. Optionally, this method can also be applied to an implementation environment composed of computer devices and servers.

[0024] The landslide deformation estimation method provided in the embodiments of this application will be described in detail below with reference to the accompanying drawings. See also Figure 1 The diagram shows a flowchart of a landslide deformation estimation method provided in an embodiment of this disclosure. The method includes the following steps S101 to S104: S101, acquire multiple differential interferograms of the target study area and the coherence coefficient diagram corresponding to each differential interferogram.

[0025] Understandably, the target study area refers to the geographical region where landslide deformation monitoring and estimation are required. This area can be a densely vegetated mountainous area, a hilly area with fractured terrain, or an area prone to landslide geological hazards, such as mining areas or the vicinity of reservoirs. The size of the target study area is determined according to actual monitoring needs, and can be a local area of ​​tens of square kilometers or a wide area covering hundreds or even thousands of square kilometers.

[0026] Here, the differential interferogram of the target study area is a data product generated by synthetic aperture radar interferometry. It records the remaining interferometric phase information between two synthetic aperture radar images of the same target study area acquired at different times, after removing the flat-field effect and terrain phase. This interferometric phase exists in the form of a wrapped phase, and its numerical range is limited to -π to π, reflecting the small displacement of the ground surface along the radar line of sight. The differential interferogram can be obtained by performing differential interferometric processing on multi-temporal synthetic aperture radar images covering the target study area. In the specific implementation process, suitable image pairs can be selected from the multi-temporal synthetic aperture radar images of the target study area. Then, each image pair can be precisely registered, interferometrically generated, and processed to remove the flat-field effect and terrain phase, thereby generating a differential interferogram corresponding to each image pair.

[0027] Here, the coherence coefficient map corresponding to each differential interferogram refers to a two-dimensional image generated by calculating the coherence between corresponding pixels in two synthetic aperture radar (SAR) images. The coherence coefficient map quantifies the reliability of the interference phase quality at each pixel location. Its value typically ranges from 0 to 1; a value closer to 1 indicates higher interference phase quality, while a value closer to 0 indicates lower interference phase quality. The coherence coefficient map is obtained during differential interferometry processing by calculating the coherence within a local window of the registered SAR images.

[0028] Among them, the surface observation data acquired through synthetic aperture radar (SAR) systems is an active microwave imaging radar that obtains backscattering information of surface targets by emitting microwave pulses towards the surface and receiving the echo signals. SAR images record the echo signal intensity and phase information of surface targets within each resolution cell, and can be acquired through spaceborne or airborne SAR systems.

[0029] For example, when acquiring a differential interferogram, one image can be selected as the master image from the multi-temporal synthetic aperture radar (SAR) images of the target study area, and the remaining images can be registered with the master image to ensure all images are in a unified geometric coordinate system. Then, interferometry is performed on each pair of registered master and slave images to generate an interferogram. Each pixel in the interferogram records the phase difference between the master and slave images at that pixel location. Next, the flat-land effect phase and terrain phase caused by the Earth's curvature and topographic undulations are removed from the interferogram to obtain the differential interferogram. This differential interferogram retains the phase contribution caused by surface deformation, as well as other non-deformation phase contributions such as atmospheric delay and noise. Correspondingly, the coherence coefficient map is calculated synchronously during the generation of the differential interferogram. That is, when interferometry is performed on the two registered SAR images, a rectangular window is selected around each pixel location, and the coherence coefficient at that pixel location is calculated using the complex signal values ​​of all pixels within the window. After traversing all pixel locations, the coherence coefficient map corresponding to the entire differential interferogram is obtained.

[0030] In some possible embodiments, a geocoded file of the target study area can also be obtained to convert the subsequently generated deformation rate field from the radar coordinate system to the geographic coordinate system, so as to perform spatial overlay analysis with geographic information data.

[0031] S102, for each of the differential interferograms, determine the deformation gradient field of the differential interferogram in multiple preset gradient extraction directions, and determine the weight field corresponding to each deformation gradient field based on the coherence coefficient map corresponding to the differential interferogram.

[0032] Understandably, in order to extract spatial variation information related to surface deformation from differential interferograms (DIAs) and avoid unwrapping the entangled phase, deformation gradient fields can be calculated separately along multiple gradient extraction directions of the DIA. The deformation gradient field reflects the rate of change of deformation rate in space, i.e., the difference in deformation rate between adjacent pixels. It can effectively highlight local deformation signals while suppressing large-scale, slowly changing atmospheric delays and random noise.

[0033] Here, the deformation gradient field refers to the physical field characterizing the rate of change of the surface deformation rate in the target study area in space. The value of each pixel reflects the difference in deformation rate between that pixel and its neighboring pixels along a specified direction. The preset multiple gradient extraction directions represent the directions selected on the two-dimensional plane for calculating the spatial rate of change. These are the basic spatial dimensions for subsequent deformation field reconstruction and can be specifically set according to actual needs, such as the range direction, azimuth direction, and the two diagonal directions between the range and azimuth directions. Specifically, the range direction corresponds to the column direction of the synthetic aperture radar image, i.e., the projection direction of the radar beam's line of sight onto the horizontal plane; the azimuth direction corresponds to the row direction of the synthetic aperture radar image, i.e., the satellite's flight direction; and the diagonal directions correspond to the two directions with an angle of 45 degrees between the range and azimuth directions.

[0034] For example, considering that surface deformation may change in any direction in space, calculating the gradient along only a single direction cannot fully describe the two-dimensional spatial distribution characteristics of the deformation. Therefore, when determining the deformation gradient field of each differential interferogram, reference is made to... Figure 2 As shown, for each gradient extraction direction, all pixel positions of the differential interferogram are traversed, and the following operations S201~S204 are performed for each pixel position: S201, in the complex phase field of the differential interferogram, determine the first adjacent pixel and the second adjacent pixel that are located on both sides of the pixel position along the current gradient extraction direction and are respectively separated from the pixel position by a preset step size.

[0035] Here, the complex phase field of the differential interferogram is represented as a complex signal field constructed by taking the wrapped phase value of each pixel in the differential interferogram as the argument. Each complex signal unit in this complex signal field is obtained by mapping the wrapped phase value of the pixel to the complex domain using Euler's formula. The complex phase field of the differential interferogram is the basic input data for subsequent conjugate multiplication operations, which can preserve complete phase information without unwrapping. The preset step size represents the offset used when shifting pixels along the current gradient extraction direction on the differential interferogram grid. The spatial scale of gradient calculation can be adjusted by setting the value of the step size; for example, the step size can be one pixel, two pixels, or more pixels.

[0036] Specifically, before determining the first and second adjacent pixels, the row and column coordinates of the current pixel position within the entire differential interferogram grid can be determined, along with the row and column offsets corresponding to the current gradient extraction direction. Using the current pixel position as the center, offsets are made along both the positive and negative directions of the current gradient extraction direction by a preset step size to obtain the first and second adjacent pixels located on either side of the current pixel position. Thus, the first and second adjacent pixels are located on either side of the current pixel position along the current gradient extraction direction, with distances from the current pixel position equal to the preset step size.

[0037] For example, in the distance direction, the column coordinates of the current pixel position can be subtracted and added by a preset step size, while the row coordinates remain unchanged, to obtain the first and second adjacent pixels; in the azimuth direction, the row coordinates of the current pixel position can be subtracted and added by a preset step size, while the column coordinates remain unchanged, to obtain the first and second adjacent pixels; in the northeast-southwest diagonal direction, the row coordinates of the current pixel position can be subtracted by a preset step size and the column coordinates can be added by a preset step size, and the row coordinates of the current pixel position can be added by a preset step size and the column coordinates can be subtracted by a preset step size, to obtain the first and second adjacent pixels; in the southeast-northwest diagonal direction, the row coordinates of the current pixel position can be added by a preset step size and the column coordinates can be added by a preset step size, and the row coordinates of the current pixel position can be subtracted by a preset step size and the column coordinates can be subtracted by a preset step size, to obtain the first and second adjacent pixels.

[0038] In some possible embodiments, if the current pixel position is located in the boundary region of the differential interferogram, causing it to exceed the image boundary after shifting by a preset step in a certain direction and thus failing to obtain effective adjacent pixels, a cyclic shift method can be used to obtain the corresponding pixel value from the opposite side boundary of the image, so that the boundary pixel position can also obtain the complete first and second adjacent pixels, ensuring that all pixel positions can participate in gradient calculation, thereby ensuring that the final output deformation gradient field has no missing data.

[0039] S202, perform conjugate multiplication of the complex phase value of the first adjacent pixel and the complex phase value of the second adjacent pixel, and determine the phase difference between the first adjacent pixel and the second adjacent pixel based on the conjugate multiplication result.

[0040] Specifically, conjugate multiplication is a mathematical operation used to eliminate common phase components between two complex signals. It involves multiplying the complex phase value of a first adjacent pixel by the conjugate of the complex phase value of a second adjacent pixel. This conjugate multiplication eliminates the common phase components between the two pixels caused by common phase contribution sources. The true phase difference between the first and second adjacent pixels can be determined based on the argument of the conjugate multiplication result. This process does not require unwrapping the entangled phases. Here, the phase difference refers to the difference in phase values ​​between the first and second adjacent pixels caused by factors such as surface deformation and atmospheric delay; it is the minimum difference between the entangled phase values ​​of these two pixels.

[0041] In some possible embodiments, the phase difference can be obtained by extracting the argument of the conjugate multiplication result. The argument can be extracted by the arctangent function or the complex argument function, without being specifically limited here.

[0042] For example, taking the east-west direction as an example, the result of the conjugate multiplication in the east-west direction can be expressed as: ; in, It is represented as the complex result of conjugate multiplication along the east-west direction; This is represented as the complex phase value obtained by shifting the complex phase field to the left along the column direction by a preset step size. It represents the complex phase value obtained by shifting the complex phase field to the right along the column direction by a preset step size; x represents the column coordinate in the differential interferogram grid, corresponding to the range direction; y represents the row coordinate in the differential interferogram grid, corresponding to the azimuth direction.

[0043] For example, the phase difference in the east-west direction can be expressed as: ; in, It represents the phase difference in the east-west direction; arg(·) represents the operation of taking the complex argument.

[0044] S203, based on the phase difference and the time baseline of the differential interferogram, determine the deformation rate gradient value between the first adjacent pixel and the second adjacent pixel, and use the deformation rate gradient value as the gradient value corresponding to the pixel position in the deformation gradient field corresponding to the current gradient extraction direction.

[0045] Here, the time baseline is the time difference between the acquisition times of the two synthetic aperture radar images that make up the differential interferogram. It can be represented by a set time unit, such as year, month or day. The specific time unit used can be determined according to the actual processing needs and the application scenario of the output results.

[0046] For example, when it is necessary to output the deformation rate gradient value in millimeters per year, the time baseline can be expressed in years. For example, the revisit period of a single Sentinel-1 satellite is 12 days, which is approximately 0.0329 years. When the monitoring time span is short and it is necessary to output the daily or monthly variation rate, the time baseline can also be expressed in days or months, without any specific limitation.

[0047] Specifically, based on the phase difference, the time baseline, and a preset conversion coefficient between phase and deformation, the deformation rate gradient value between the first and second adjacent pixels can be determined. The phase difference characterizes the phase difference between two pixels. Multiplying it by the preset conversion coefficient converts it into a deformation difference value. Dividing this by the time baseline converts the deformation difference value into a deformation rate difference value per unit time, i.e., the deformation rate gradient value. The deformation rate gradient value characterizes the rate of change of deformation rate between the first and second adjacent pixels in space. After obtaining this value, it can be used as the gradient value corresponding to the current pixel position in the deformation gradient field corresponding to the current gradient extraction direction, for subsequent multi-directional gradient field weighted averaging and deformation field reconstruction.

[0048] The preset conversion coefficient between phase and deformation is a proportional factor that converts the phase change in the radar line-of-sight direction into the surface deformation displacement. It can be expressed as wavelength / 4π and can establish a quantitative relationship between phase difference and surface deformation displacement. Specifically, it can be set according to the wavelength of the synthetic aperture radar sensor used. For example, the radar wavelength of Sentinel-1 satellite is 5.6 cm, and the corresponding conversion coefficient is about 45 mm per radian.

[0049] For example, taking the east-west direction as an example, the calculation of the east-west deformation rate gradient value can be expressed as: ; in, It is expressed as the deformation rate gradient value in the east-west direction, in millimeters per year; The time difference between the acquisition times of the two synthetic aperture radar images that constitute the differential interferogram is expressed in years. It represents the preset conversion coefficient between phase and deformation, used to convert phase difference into deformation displacement.

[0050] S204, After the traversal is completed, the deformation gradient field of the differential interferogram in the current gradient extraction direction is obtained.

[0051] Furthermore, after all pixel positions of the differential interferogram have been traversed, the gradient values ​​of all pixel positions of the differential interferogram in the current gradient extraction direction can be calculated, thereby obtaining the complete deformation gradient field corresponding to the differential interferogram and the gradient extraction direction. The deformation gradient field is a two-dimensional data array with the same spatial size as the differential interferogram, which can characterize the spatial variation distribution of the surface deformation rate reflected by the differential interferogram in the specified direction.

[0052] In this way, by performing gradient calculations along multiple gradient extraction directions on the complex phase field, the deformation gradient fields and corresponding weight fields of each differential interferogram in each direction can be obtained. This allows for the acquisition of spatial gradient information of surface deformation without relying on phase unwrapping, reducing phase unwrapping errors and their propagation risks.

[0053] Furthermore, after determining the deformation gradient field of each differential interferogram in each gradient extraction direction, since different pixel positions of each deformation gradient field have different phase quality reliability, the weight field corresponding to each deformation gradient field can also be determined based on the coherence coefficient map corresponding to the differential interferogram. The weight field is used to assign corresponding weights to the gradient values ​​at different positions in the deformation gradient field, which can reflect the reliability of the gradient values ​​at different pixel positions, thereby reducing the contribution of unreliable gradient data in the subsequent weighted average processing.

[0054] Specifically, when determining the weight field, refer to Figure 3 As shown, for each gradient extraction direction, all pixel positions of the differential interferogram are traversed, and the following operations S301~S304 are performed for each pixel position: S301, determine the first adjacent pixel and the second adjacent pixel corresponding to the pixel position.

[0055] Understandably, according to the above step S201, the first and second adjacent pixels corresponding to each pixel position can be determined. Since the weight value needs to be determined one-to-one with the pixel position in the deformation gradient field when determining the weight value, the first and second adjacent pixels corresponding to the current pixel position can be determined first before determining the weight value, so as to serve as the basis for finding the position of the coherence coefficient in the coherence coefficient map later.

[0056] S302, obtain the first coherence coefficient corresponding to the first adjacent pixel and the second coherence coefficient corresponding to the second adjacent pixel from the coherence coefficient map.

[0057] Here, the value of each pixel in the coherence coefficient graph reflects the reliability of the interference phase quality at that location, containing coherence distribution information for the entire target study area. Therefore, when determining the weight values, the coherence coefficients at the same positions as the first and second adjacent pixels in the coherence coefficient graph can be identified and designated as the first coherence coefficient and the second coherence coefficient, respectively. The first coherence coefficient is the coherence value corresponding to the first adjacent pixel position in the coherence coefficient graph, which can characterize the interference phase quality at the location of the first adjacent pixel; the second coherence coefficient is the coherence value corresponding to the second adjacent pixel position in the coherence coefficient graph, which can characterize the interference phase quality at the location of the second adjacent pixel.

[0058] S303, the minimum value between the first coherence coefficient and the second coherence coefficient is determined as the weight value corresponding to the pixel position in the weight field corresponding to the current gradient extraction direction.

[0059] It is understandable that the minimum value of the first coherence coefficient and the second coherence coefficient represents the coherence level of the poorer quality pixel among the first and second adjacent pixels. Considering that the calculation of the deformation gradient value depends on the complex phase values ​​of the first and second adjacent pixels, if the interference phase quality of either of these two pixels is low, the gradient value calculated based on these two pixels will be unreliable. Therefore, when determining the weight value, the minimum value of the first coherence coefficient and the second coherence coefficient can be used as the weight value of the current position to reduce the contribution of unreliable gradient data in subsequent processing.

[0060] For example, the value of the coherence coefficient can be determined through a minimization strategy. Taking the east-west direction mentioned above as an example, the weight value can be calculated as follows: ; in, This is represented as the weight value corresponding to the current pixel position in the weight field corresponding to the east-west direction; It is represented as the first coherence coefficient corresponding to the first adjacent pixel position in the coherence coefficient diagram; This is represented as the second coherence coefficient corresponding to the second adjacent pixel position in the coherence coefficient diagram; This means taking the minimum of the two values ​​within the parentheses.

[0061] S304, After the traversal is completed, the weight field of the differential interferogram in the current gradient extraction direction is obtained.

[0062] In this way, by traversing all pixel positions, the weight values ​​of all pixel positions in the differential interferogram along the current gradient extraction direction can be determined, resulting in a complete weight field corresponding to the differential interferogram and the gradient extraction direction. The weight field is a two-dimensional data array with the same spatial dimensions as the deformation gradient field, providing corresponding mass weights for the gradient value at each pixel position in the deformation gradient field.

[0063] In this disclosure, by determining the minimum value of the coherence coefficients of the first and second adjacent pixels as the weight of the corresponding gradient value, the influence of unreliable gradient data in low coherence regions on the average deformation gradient field can be effectively suppressed during subsequent weighted averaging, thereby improving the reliability of the gradient field data and the accuracy of the subsequent deformation field reconstruction results.

[0064] It should be noted that the deformation gradient field and weight field in other gradient extraction directions in this disclosure can be calculated according to the above formula. It is only necessary to substitute the positions of the first and second adjacent pixels in the corresponding directions into the above calculation process, which will not be elaborated here.

[0065] S103, based on the deformation gradient field and corresponding weight field of all the differential interferograms in the same gradient extraction direction, determine the average deformation gradient field of the target study area in the gradient extraction direction.

[0066] Specifically, the average deformation gradient field is a comprehensive gradient field obtained by weighted averaging of the deformation gradient fields of all differential interferograms along the same gradient extraction direction. This eliminates random errors introduced by atmospheric delay, noise, and other factors in a single differential interferogram. Specifically, the weighted average value at each pixel location can be calculated based on the deformation gradient fields of all differential interferograms along the same gradient extraction direction and their corresponding weight fields. Since each differential interferogram corresponds to a different time baseline, and the coherence coefficients of each pixel location differ, the gradient field of a single differential interferogram inevitably contains noise and local errors. Weighted averaging effectively suppresses these errors, making the average deformation gradient field more accurately reflect the true and stable deformation gradient distribution of the target study area in that direction.

[0067] For example, for each gradient extraction direction, when weighted and fused the deformation gradient fields and corresponding weight fields of all differential interferograms in that direction, refer to Figure 4 As shown, the steps S401~S403 may be included: S401, for each gradient extraction direction, the deformation gradient fields of all the differential interferograms in the gradient extraction direction are multiplied pixel by pixel with their respective weight fields, and the multiplication results of all the differential interferograms are accumulated pixel by pixel to obtain the gradient accumulation field.

[0068] Here, for each differential interferogram, the deformation gradient field and the weight field are multiplied pixel by pixel. This weight field weights the gradient value at each pixel location in the deformation gradient field. Furthermore, the weighted results for the same pixel location in all differential interferograms are summed to obtain the weighted gradient sum of that pixel location across all interferograms. Finally, the data field obtained after performing the same operation on all pixel locations can be used as the gradient accumulation field. The gradient accumulation field is a two-dimensional data array with the same spatial size as a single differential interferogram, and it is the weighted gradient sum of all differential interferograms along the same gradient extraction direction.

[0069] For example, the expression for calculating the gradient accumulation field in the east-west direction can be expressed as: ; in, This is represented as a gradient accumulation field along the east-west direction; It is represented as the deformation gradient field in the east-west direction of the i-th differential interferogram; Let N represent the weight field of the i-th differential interferogram in the east-west direction; N represents the total number of differential interferograms in the east-west direction; and i represents the index number of the differential interferogram.

[0070] S402, the weight fields of all the differential interferograms in the gradient extraction direction are accumulated pixel by pixel to obtain the weight accumulation field.

[0071] Specifically, to normalize each pixel position during subsequent calculation of the average deformation gradient field, the total weight of each pixel position across all differential interferograms can be determined. This involves pixel-by-pixel summing of the weight fields of all differential interferograms along the same gradient extraction direction, resulting in a weight accumulation field to determine the total weight of each pixel position. The weight accumulation field is a two-dimensional data array with the same spatial dimensions as the gradient accumulation field, recording the sum of the weights of each pixel position across all differential interferograms.

[0072] For example, for the east-west direction, the expression for calculating the weighted cumulative field can be expressed as: ; in, It is represented as a weighted cumulative field in the east-west direction.

[0073] S403, based on the gradient accumulation field and the weight accumulation field, determine the average deformation gradient field of the target study area in the gradient extraction direction.

[0074] Furthermore, after obtaining the gradient accumulation field and the weight accumulation field, the gradient accumulation field can be used as the sum of weighted gradients, and the weight accumulation field can be used as the normalized denominator of the weighted sum. By dividing the gradient accumulation field and the weight accumulation field pixel by pixel, the average deformation gradient field of the target study area in the gradient extraction direction can be determined. Each pixel value of this average deformation gradient field represents the weighted average gradient value at that location across all available differential interferograms.

[0075] For example, the calculation of the average deformation gradient field in the east-west direction can be expressed as: ; in, It is represented as the average deformation gradient field in the east-west direction.

[0076] In this way, the average deformation gradient field integrates the gradient information of all differential interferograms along the gradient extraction direction, has a higher signal-to-noise ratio and more stable statistical representativeness, and can provide an accurate and reliable data foundation for subsequent deformation rate field reconstruction.

[0077] In some possible embodiments, to improve the data quality of the average deformation gradient field, a filtering and smoothing process can be performed on the average deformation gradient field. Specifically, a two-dimensional median filter is applied to each average deformation gradient field and its corresponding average coherence coefficient field. The median filter uses a sliding window approach, taking the median value of the pixels within the window as the filtering result for the center pixel. The window size is typically set to 15×15 or 21×21 pixels, and the specific window size can be set according to the spatial resolution of the target study area and the spatial distribution characteristics of the deformation field. The average coherence coefficient field is a reliability evaluation field corresponding to the average deformation gradient field. It can be obtained by arithmetically averaging the weight fields of all differential interferograms along the same gradient extraction direction, and is used to measure the reliability of the average deformation gradient field at each pixel location. In this way, median filtering can effectively remove isolated noise points, eliminate gradient outliers caused by local noise, and maintain the sharpness of deformation edges, avoiding excessive smoothing that erases the true deformation boundaries, thus preserving the spatial detail features of the deformation field.

[0078] In some possible embodiments, the filtered data can also undergo quality screening, i.e., low-quality regions can be removed using a coherence coefficient mask. A coherence coefficient threshold is set, for example, 0.2, and the gradient values ​​in the average deformation gradient field corresponding to pixels with an average coherence coefficient below this threshold are set to invalid values. In this way, data from areas with severe phase decoherence, such as water bodies, steep terrain, or dense vegetation, can be removed, preventing unreliable gradient values ​​from these areas from introducing systematic errors in the subsequent Poisson equation solution, thereby improving the accuracy of deformation field reconstruction.

[0079] In some possible embodiments, invalid value gaps generated after masking can also be filled. For example, a neighborhood interpolation algorithm can be used to fill the invalid value locations. During interpolation, only the gradient information of surrounding valid pixels is used for estimation, so that the filling result conforms to the continuous change trend of local deformation. In this way, the global data continuity of the gradient field can be guaranteed, so as to meet the requirement that the input data has no gaps in the subsequent solution of the Poisson equation.

[0080] S104, using the Poisson equation, the deformation rate field is reconstructed based on the average deformation gradient field in each of the gradient extraction directions to obtain the deformation rate field of the target study area; and the deformation estimation result of the target study area is determined based on the deformation rate field of the target study area.

[0081] It is understandable that the Poisson equation is a partial differential equation used to reconstruct the antiderivative from a gradient field. Given a gradient field, solving the Poisson equation can yield the globally optimal solution that best matches the input gradient field in the least squares sense. This disclosure proposes using the Poisson equation to extract the average deformation gradient field along each gradient direction, and then performing spatial integration reconstruction on the gradient domain information. This allows the deformation spatial variation information contained in the gradient field to be restored to a complete deformation rate field, thus obtaining a spatially continuous deformation rate field for the target study area with definite physical units.

[0082] Specifically, the average deformation gradient fields along different gradient extraction directions describe the spatial variation characteristics of the deformation rate from different dimensions. When constructing the Poisson equation, one can refer to... Figure 5 As shown, in the process of reconstructing the deformation rate field, the average deformation gradient fields in two mutually perpendicular target directions can be selected as the basic inputs for solving the Poisson equation. Specifically, this can include the following steps S501~S505: S501, obtain the average deformation gradient field in the first target direction and the average deformation gradient field in the second target direction from the average deformation gradient fields in each of the gradient extraction directions.

[0083] Here, the first target direction is one of the two mutually perpendicular directions involved in the reconstruction of the Poisson equation, and the range direction can be chosen. The second target direction represents the other of the two mutually perpendicular directions involved in the reconstruction of the Poisson equation, which is the azimuth direction. The first and second target directions are perpendicular to each other, providing complete two-dimensional spatial gradient information for the Poisson equation, thus ensuring a unique solution to the reconstruction problem from the gradient field to the deformation field.

[0084] S502, take the partial derivative of the average deformation gradient field in the first target direction in the first target direction to obtain the first partial derivative field; and take the partial derivative of the average deformation gradient field in the second target direction in the second target direction to obtain the second partial derivative field.

[0085] Specifically, by taking the partial derivative of the average deformation gradient field along the first target direction, the spatial rate of change of the gradient field along the first direction can be calculated, yielding the first partial derivative field. Similarly, by taking the partial derivative of the average deformation gradient field along the second target direction, the second partial derivative field can be obtained, which characterizes the gradient rate of change in the second direction.

[0086] Here, taking the partial derivative of the average deformation gradient field along the first target direction is equivalent to calculating the spatial rate of change of this gradient field along the first direction, i.e., the second-order change information of the deformation rate field to be reconstructed along the first direction. Correspondingly, in the process of constructing the Poisson equation, the differential relationship between the deformation rate field to be reconstructed and the average deformation gradient field can be expressed as: ; ; Where u represents the deformation rate field of the target study area to be reconstructed; This is expressed as the average deformation gradient field along the first target direction; It is represented as the average deformation gradient field in the direction of the second target.

[0087] Specifically, the above formula describes the spatial gradient of the deformation rate field u and the average deformation gradient field. , The theoretical relationship between the two is that, ideally, the partial derivative of the deformation rate field along the first target direction should be equal to the average deformation gradient field along the first target direction, and the partial derivative along the second target direction should be equal to the average deformation gradient field along the second target direction. However, in actual data, due to noise and errors, the above equations do not always hold true. Therefore, it is necessary to obtain the deformation rate field reconstruction result that best matches the input gradient field using the Poisson equation in the least squares sense.

[0088] S503, Based on the first partial-derived field and the second partial-derived field, determine the target divergence field.

[0089] Here, the target divergence field is a two-dimensional data field obtained by adding the corresponding pixels of the first partial derivative field and the second partial derivative field. It can characterize the spatial variation intensity of the average deformation gradient field and is mainly used as the source term of the Poisson equation to drive the reconstruction of the deformation field.

[0090] For example, the calculation of the target divergence field can be expressed as: ; in, Represented as the target divergence field, it is a two-dimensional data field obtained by adding the corresponding pixels of the first partial derivative field and the second partial derivative field; Here, to calculate the above divergence on a discrete grid, a central difference approach can be used for approximation, as shown in the following equation: ; in, It is expressed as the partial derivative of the first partial derivative field at pixel position (i,j) along the first target direction; It is represented as the gradient value of the first partial derivative field at pixel position (i,j+1) at the adjacent pixel along the positive direction of the first target direction; It is represented as the gradient value of the first partial derivative field at pixel position (i,j-1) at the pixel adjacent to the first target direction in the negative direction; This represents the ground resolution corresponding to the pixel interval, used to convert the pixel coordinate difference into actual spatial distance. Similarly, the partial derivative of the second partial derivative field in the direction of the second target can be discretized and approximated in the same way, which will not be elaborated here.

[0091] S504, using the target divergence field as the source term of the Poisson equation and the deformation rate field to be reconstructed as the unknown quantity to be solved in the Poisson equation, construct the Poisson equation.

[0092] Specifically, the target divergence field can be used as the right-hand source term of the Poisson equation, and the Laplace operator of the deformation rate field to be reconstructed can be used as the left-hand term of the equation to construct the complete Poisson equation. This allows for the establishment of a quantitative constraint relationship between the deformation rate field and the average deformation gradient field, thereby obtaining a numerical solution for the deformation rate field. This enables the gradient domain information to be restored to a continuous deformation rate distribution in the spatial domain.

[0093] Here, the Poisson equation can be expressed as follows: the Laplace operator of the deformation rate field is equal to the target divergence field: ; in, It is expressed as the Laplace operator of the deformation rate field to be reconstructed, that is, the sum of the second-order partial derivatives of the deformation rate field in space; It is expressed as the second-order partial derivative of the deformation rate field along the direction of the first target; It is expressed as the second partial derivative of the deformation rate field along the direction of the second target; the right side of the equation is the target divergence field; It is expressed as the partial derivative of the average deformation gradient field along the first target direction, that is, the spatial rate of change of the first partial derivative field along the first target direction; This is expressed as the partial derivative of the average deformation gradient field along the second target direction, i.e., the spatial rate of change of the second partial derivative field along the second target direction. This equation establishes a quantitative constraint relationship between the deformation rate field to be reconstructed and the input gradient field.

[0094] S505, Solve the Poisson equation based on the preset boundary conditions to obtain the deformation rate field of the target study area.

[0095] Furthermore, the Poisson equation can be solved according to preset boundary conditions. Under given boundary constraints, the deformation rate field distribution satisfying the Poisson equation can be obtained, thus yielding the deformation rate field of the target study area. The preset boundary conditions are constraints applied to the boundary of the solution region when solving the Poisson equation, ensuring a unique solution. These preset boundary conditions can be set based on the boundary characteristics of the solution region and the physical properties of the deformation field, such as Neumann boundary conditions, Dirichlet boundary conditions, or periodic boundary conditions; no specific limitations are imposed here.

[0096] Specifically, when numerically solving the Poisson equation, considering the large computational cost of directly solving partial differential equations in the spatial domain, the Poisson equation can be transformed from the spatial domain to the frequency domain for solution, including the following steps (a) to (b): (a) Transform the Poisson equation from the spatial domain to the frequency domain, and solve the Poisson equation in the frequency domain based on the preset boundary conditions to obtain the frequency domain deformation coefficients; (b) The frequency domain deformation coefficients are inversely transformed from the frequency domain back to the spatial domain to obtain the deformation rate field of the target study area.

[0097] Here, by transforming the Poisson equation from the spatial domain to the frequency domain, partial differential operations in the spatial domain can be converted into algebraic operations in the frequency domain, thus significantly reducing computational complexity. Furthermore, by solving the Poisson equation in the frequency domain based on preset boundary conditions, the deformation field coefficients in the frequency domain can be obtained efficiently. Further, after obtaining the frequency domain deformation coefficients, they can be inversely transformed back from the frequency domain to the spatial domain to obtain the deformation rate value at each pixel location in the spatial domain, thereby obtaining the complete deformation rate field of the target study area.

[0098] For example, taking the preset boundary conditions as von Neumann boundary conditions and solving the Poisson equation through discrete cosine transform as an example: First, a two-dimensional discrete cosine transform can be performed on the target divergence field to transform it to the frequency domain; then, eigenvalues ​​of the Laplace operator can be constructed in the frequency domain based on the von Neumann boundary conditions, and the frequency domain deformation coefficients can be obtained through frequency domain operations; finally, a two-dimensional inverse discrete cosine transform can be performed on the frequency domain deformation coefficients to inversely transform them from the frequency domain back to the spatial domain, thereby obtaining the deformation rate field of the target study area.

[0099] Here, the Neumann boundary condition, i.e., the normal derivative on the boundary is zero, is: ; Where n represents the outward normal direction of the solution region boundary; Let the target study area be represented as the spatial domain to be solved; This is expressed as the directional derivative of the deformation rate field along the outward normal direction of the boundary. This boundary condition indicates that, at the boundary of the solution domain, the rate of change of the deformation rate field along the outward normal direction is zero, meaning that the deformation rate field does not change in the normal direction at the boundary.

[0100] Here, by solving the Poisson equation using the discrete cosine transform, the eigenvalues ​​of the Laplace operator in the discrete cosine transform domain can be expressed as: ; In the formula, M and N are the grid sizes, corresponding to the total number of rows and columns of the differential interferogram, respectively; k and l are the frequency domain indices, corresponding to the frequency components in the row and column directions, respectively. It is represented as the eigenvalue of the Laplace operator at position (k,l) in the eigenvalue matrix of the discrete cosine transform domain.

[0101] Furthermore, after obtaining the eigenvalues, a two-dimensional discrete cosine transform is performed on the target divergence field to obtain the frequency domain divergence coefficients. Then, the frequency domain deformation coefficients can be obtained through frequency domain operations. ; in, This is represented as the value at position (k,l) in the frequency domain deformation coefficient matrix; This represents the frequency domain divergence coefficient at position (k,l) in the frequency domain after the target divergence field undergoes a two-dimensional discrete cosine transform. When k=l=0, it corresponds to the zero frequency component (i.e., the constant component) in the frequency domain. Setting it to zero indicates the overall translation of the deformation rate field, ensuring a unique solution for the reconstruction result. For other frequency components, the deformation coefficient corresponding to that frequency component can be obtained by dividing the frequency domain divergence coefficient by the eigenvalue at the corresponding position and taking the negative value.

[0102] Furthermore, regarding the frequency domain deformation coefficients By performing a two-dimensional inverse discrete cosine transform, the deformation coefficients in the frequency domain can be restored to the deformation rate values ​​at each pixel location in the spatial domain, thereby obtaining the deformation rate field of the target study area.

[0103] In this way, by fusing and reconstructing multi-directional gradient information into a globally consistent deformation rate field through the Poisson equation, a quantitative transformation from a dimensionless gradient field to a deformation rate field with definite physical units is achieved, avoiding the problems of error accumulation and discontinuity caused by pixel-by-pixel processing.

[0104] In some other embodiments, the Poisson equation can also be solved using fast Fourier transform, sine transform, or other frequency domain solving methods, without being specifically limited here.

[0105] Furthermore, after obtaining the deformation rate field of the target study area, the deformation estimation results of the target study area can be determined based on the deformation rate field, so as to serve as a quantitative basis for the identification and monitoring of landslide geological hazards.

[0106] Specifically, when determining the deformation estimation result based on the deformation rate field, the following steps (1) to (3) may be included: (1) Take the deformation rate value of each pixel in the deformation rate field of the target study area as the deformation rate of the pixel; (2) Based on the deformation rate of all the pixels, identify the deformation region in the target study area, the boundary of the deformation region, and the deformation rate distribution inside the deformation region; (3) The identified deformation region, the boundary of the deformation region, and the deformation rate distribution inside the deformation region are used as the deformation estimation result of the target study area.

[0107] It is understandable that the deformation rate value of each pixel in the deformation rate field can be used as a quantitative indicator of the speed of surface deformation at that location, in order to distinguish between active and stable deformation areas. At the same time, based on the continuous spatial distribution characteristics of the deformation rates of all pixels, the spatial range of the deformation area, the boundary position between the deformation area and the surrounding stable area, and the relative magnitude and spatial variation law of deformation rates between different positions within the deformation area can be determined, thereby comprehensively understanding the overall activity status and internal movement differences of the landslide body.

[0108] In this way, the reconstructed deformation rate field is transformed into a deformation estimation result that includes the deformation region, boundary location, and internal rate distribution. This can fully present the spatial distribution characteristics of surface deformation in the target study area and provide data support for early identification, hazard assessment, and prevention and control decisions of landslide disasters.

[0109] The landslide deformation estimation method, apparatus, storage medium, and computer equipment provided in this disclosure achieve effective characterization of the spatial distribution of surface deformation in the target study area by performing weighted averaging based on the coherence coefficient weight field in the multi-gradient extraction direction and reconstructing the deformation rate field using the Poisson equation, thereby improving the reliability and spatial continuity of the deformation estimation results.

[0110] To verify the effectiveness and superiority of the landslide deformation estimation method proposed in this disclosure, 96 synthetic aperture radar (SAR) images covering a target study area were obtained from a satellite platform for experimental verification. This satellite is a C-band SAR satellite, capable of full coverage and periodic observation of the global land surface. Furthermore, due to its high temporal and spatial resolution and free open access policy, this satellite data can become a commonly used data source in the field of SAR interferometry. Here, after registering all SAR images to the master image, 447 interferometric pairs were generated. Differential interferometry was then performed on the registered SAR images to obtain differential interferograms, and then a coherence diagram corresponding to each differential interferogram was generated.

[0111] Furthermore, in order to compare and verify the technical effect of the method disclosed herein, the deformation estimation results of the target study area were obtained by using the traditional phase gradient rate (PGR) method of the central difference method, the improved phase gradient stacking (IPGS) method, and the method mentioned in this disclosure, respectively.

[0112] Specifically, refer to Figure 6 As shown, Figure 6 (a) is a landslide extent map of the target study area, used to indicate the spatial distribution range of landslide bodies in the area, serving as a reference benchmark for comparing deformation estimation results; Figure 6 (b) is a deformation result obtained by the traditional phase gradient rate method using the center of difference method. It can be seen that this method can only reflect the contour information of the landslide body boundary, but can hardly present the deformation characteristics inside the landslide body. Figure 6 (c) is the deformation result obtained by the improved phase gradient stacking method. This method is better at capturing the overall deformation trend and the boundary is clearer, but it is still sensitive to gradient changes and uniform deformation areas. The continuity of the deformation field and the ability to express internal details are limited. Figure 6 (d) is a deformation result obtained by the method mentioned in this disclosure. This method can quickly and robustly generate a complete deformation rate field of the landslide body without phase unwrapping. The deformation region boundary is clear, and both the internal uniform deformation region and the slow deformation region can be accurately identified. It is almost unaffected by the gradient step size or uniform deformation.

[0113] Here, comparing the three figures, it can be seen that the traditional phase gradient rate method (differential center method) mainly reflects the boundary characteristics of the landslide body, but it can hardly identify the internal uniform or slowly deforming regions, and the boundary becomes blurred when the gradient step size increases, so the result does not correspond to the specific physical deformation. The improved phase gradient stacking method is better at capturing the overall deformation and has clear boundaries, but it is still sensitive to gradient and uniform deformation. The method proposed in this disclosure can quickly and robustly generate the complete deformation rate field of the landslide body without unwrapping, with clear boundaries, and can accurately identify the internal uniform deformation region. It is almost unaffected by the gradient step size or uniform deformation, realizing the innovation of transforming the traditional phase gradient rate gradient stacking result into quantifiable physical deformation, which significantly improves the accuracy and reliability of landslide deformation monitoring.

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

[0115] Based on the same inventive concept, this disclosure also provides a landslide deformation estimation device corresponding to the landslide deformation estimation method. Since the principle of the device in this disclosure is similar to the landslide deformation estimation method described above, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be described again.

[0116] Reference Figure 7 The diagram shown is a schematic representation of a landslide deformation estimation device 700 provided in an embodiment of this disclosure. The device includes: The data acquisition module 701 is used to acquire multiple differential interferograms of the target study area and a coherence coefficient map corresponding to each differential interferogram; The gradient and weight field determination module 702 is used to determine the deformation gradient field of each differential interferogram in multiple preset gradient extraction directions, and to determine the weight field corresponding to each deformation gradient field based on the coherence coefficient map corresponding to the differential interferogram. The average gradient field determination module 703 is used to determine the average deformation gradient field of the target study area in the gradient extraction direction based on the deformation gradient field and the corresponding weight field of all the differential interferograms in the same gradient extraction direction. The deformation estimation module 704 is used to reconstruct the deformation rate field based on the average deformation gradient field in each of the gradient extraction directions using the Poisson equation, thereby obtaining the deformation rate field of the target study area; and to determine the deformation estimation result of the target study area based on the deformation rate field of the target study area.

[0117] In some possible embodiments, the gradient and weight field determination module 702 is specifically used for: For each gradient extraction direction, iterate through all pixel positions of the differential interferogram and perform the following operations for each pixel position: In the complex phase field of the differential interferogram, a first adjacent pixel and a second adjacent pixel are determined along the current gradient extraction direction and located on both sides of the pixel position, respectively at a preset step distance from the pixel position; The complex phase value of the first adjacent pixel and the complex phase value of the second adjacent pixel are multiplied by conjugate, and the phase difference between the first adjacent pixel and the second adjacent pixel is determined based on the result of the conjugate multiplication. Based on the phase difference and the time baseline of the differential interferogram, the deformation rate gradient value between the first adjacent pixel and the second adjacent pixel is determined, and the deformation rate gradient value is used as the gradient value corresponding to the pixel position in the deformation gradient field corresponding to the current gradient extraction direction. After the traversal is completed, the deformation gradient field of the differential interferogram in the current gradient extraction direction is obtained.

[0118] In some possible embodiments, the gradient and weight field determination module 702 is specifically used for: For each gradient extraction direction, iterate through all pixel positions of the differential interferogram and perform the following operations for each pixel position: Determine the first and second adjacent pixels corresponding to the pixel position; The first coherence coefficient corresponding to the first adjacent pixel and the second coherence coefficient corresponding to the second adjacent pixel are obtained from the coherence coefficient map. The minimum value between the first coherence coefficient and the second coherence coefficient is determined as the weight value corresponding to the pixel position in the weight field corresponding to the current gradient extraction direction; After the traversal is completed, the weight field of the differential interferogram in the current gradient extraction direction is obtained.

[0119] In some possible embodiments, the average gradient field determination module 703: For each gradient extraction direction, the deformation gradient fields of all the differential interferograms in the gradient extraction direction are multiplied pixel by pixel with their respective weight fields, and the multiplication results of all the differential interferograms are accumulated pixel by pixel to obtain the gradient accumulation field. The weight fields of all the differential interferograms in the gradient extraction direction are accumulated pixel by pixel to obtain the weight accumulation field; Based on the gradient accumulation field and the weight accumulation field, the average deformation gradient field of the target study area in the gradient extraction direction is determined.

[0120] In some possible embodiments, the deformation estimation module 704 is specifically used for: The average deformation gradient field in the first target direction and the average deformation gradient field in the second target direction are obtained from the average deformation gradient fields in each of the gradient extraction directions; wherein the first target direction and the second target direction are perpendicular to each other; A first partial derivative field is obtained by taking the partial derivative of the average deformation gradient field in the first target direction in the first target direction; and a second partial derivative field is obtained by taking the partial derivative of the average deformation gradient field in the second target direction in the second target direction. Based on the first and second partial derivative fields, the target divergence field is determined; The Poisson equation is constructed by using the target divergence field as the source term of the Poisson equation and the deformation rate field to be reconstructed as the unknown variable of the Poisson equation. The Poisson equation is solved based on preset boundary conditions to obtain the deformation rate field of the target study area.

[0121] In some possible embodiments, the deformation estimation module 704 is specifically used for: The Poisson equation is transformed from the spatial domain to the frequency domain, and the Poisson equation is solved in the frequency domain based on the preset boundary conditions to obtain the frequency domain deformation coefficients. The frequency domain deformation coefficients are inversely transformed from the frequency domain back to the spatial domain to obtain the deformation rate field of the target study area.

[0122] In some possible embodiments, the deformation estimation module 704 is specifically used for: The deformation rate value of each pixel in the deformation rate field of the target study area is taken as the deformation rate of the pixel. Based on the deformation rate of all the pixels, identify the deformation region in the target study area, the boundary of the deformation region, and the deformation rate distribution inside the deformation region; The identified deformation regions, their boundaries, and the deformation rate distribution within the deformation regions are used as the deformation estimation results for the target study area.

[0123] Based on the same technical concept, this disclosure also provides a computer device. (See also...) Figure 8The diagram shows the structure of a computer device 800 provided in this embodiment of the present disclosure, including a processor 801, a memory 802, and a bus 803. The memory 802 stores execution instructions and includes a main memory 8021 and an external memory 8022. The main memory 8021, also called internal memory, is used to temporarily store computational data in the processor 801 and data exchanged with external memory 8022 such as a hard disk. The processor 801 exchanges data with the external memory 8022 through the main memory 8021.

[0124] In this embodiment, the memory 802 is specifically used to store application code that executes the solution of this application, and its 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 foregoing embodiments.

[0125] The memory 802 may be, but is not limited to, random access memory (RAM), read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), etc.

[0126] Processor 801 may be an integrated circuit chip with signal processing capabilities. The aforementioned processor can be a general-purpose processor, including a Central Processing Unit (CPU), a Network Processor (NP), etc.; it can also be a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), a Field Programmable Gate Array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. It can implement or execute the methods, steps, and logic block diagrams disclosed in the embodiments of this invention. The general-purpose processor can be a microprocessor or any conventional processor.

[0127] It is understood that the structures illustrated in the embodiments of this application do not constitute a specific limitation on the computer device 800. In other embodiments of this application, the computer device 800 may include more or fewer components than illustrated, or combine some components, or split some components, or have different component arrangements. The illustrated components may be implemented in hardware, software, or a combination of software and hardware.

[0128] This disclosure also provides a computer-readable storage medium storing a computer program that, when executed by a processor, performs the steps of the landslide deformation estimation method described in the above-described method embodiments. The storage medium can be a volatile or non-volatile computer-readable storage medium.

[0129] This disclosure also provides a computer program product carrying program code. The program code includes instructions that can be used to execute the steps of the landslide deformation estimation method described in the above method embodiments. For details, please refer to the above method embodiments, which will not be repeated here.

[0130] The aforementioned computer program product can be implemented through hardware, software, or a combination thereof. In one optional embodiment, the computer program product is specifically embodied in a computer storage medium; in another optional embodiment, the computer program product is specifically embodied in a software product, such as a software development kit (SDK), etc.

[0131] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems and devices described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here. In the several embodiments provided in this disclosure, it should be understood that the disclosed systems and methods can be implemented in other ways. The device embodiments described above are merely illustrative. For example, the division of units is only a logical functional division; in actual implementation, there may be other division methods. Furthermore, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Another point is that the displayed or discussed mutual coupling or direct coupling or communication connection may be through some communication interfaces; the indirect coupling or communication connection of devices or units may be electrical, mechanical, or other forms.

[0132] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0133] In addition, the functional units in the various embodiments of this disclosure can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0134] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a processor-executable, non-volatile, computer-readable storage medium. Based on this understanding, the technical solution of this disclosure, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this disclosure. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0135] Finally, it should be noted that the above-described embodiments are merely specific implementations of this disclosure, used to illustrate the technical solutions of this disclosure, and not to limit them. The protection scope of this disclosure is not limited thereto. Although this disclosure has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments, or make equivalent substitutions for some of the technical features within the scope of the technology disclosed in this disclosure. Such modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this disclosure, and should all be covered within the protection scope of this disclosure.

Claims

1. A method for estimating landslide deformation, characterized in that, include: Obtain multiple differential interferograms of the target study area and a coherence coefficient map corresponding to each differential interferogram; For each differential interferogram, the deformation gradient field of the differential interferogram in multiple preset gradient extraction directions is determined, and the weight field corresponding to each deformation gradient field is determined based on the coherence coefficient map corresponding to the differential interferogram. Based on the deformation gradient field and corresponding weight field of all the differential interferograms in the same gradient extraction direction, the average deformation gradient field of the target study area in the gradient extraction direction is determined. Using the Poisson equation, the deformation rate field is reconstructed based on the average deformation gradient field in each of the gradient extraction directions to obtain the deformation rate field of the target study area; and the deformation estimation result of the target study area is determined based on the deformation rate field of the target study area.

2. The method according to claim 1, characterized in that, The step of determining the deformation gradient field of the differential interferogram in multiple preset gradient extraction directions includes: For each gradient extraction direction, iterate through all pixel positions of the differential interferogram and perform the following operations for each pixel position: In the complex phase field of the differential interferogram, a first adjacent pixel and a second adjacent pixel are determined along the current gradient extraction direction and located on both sides of the pixel position, respectively at a preset step distance from the pixel position; The complex phase value of the first adjacent pixel and the complex phase value of the second adjacent pixel are multiplied by conjugate, and the phase difference between the first adjacent pixel and the second adjacent pixel is determined based on the result of the conjugate multiplication. Based on the phase difference and the time baseline of the differential interferogram, the deformation rate gradient value between the first adjacent pixel and the second adjacent pixel is determined, and the deformation rate gradient value is used as the gradient value corresponding to the pixel position in the deformation gradient field corresponding to the current gradient extraction direction. After the traversal is completed, the deformation gradient field of the differential interferogram in the current gradient extraction direction is obtained.

3. The method according to claim 2, characterized in that, Determining the weight field corresponding to each deformation gradient field based on the coherence coefficient map corresponding to the differential interferogram includes: For each gradient extraction direction, iterate through all pixel positions of the differential interferogram and perform the following operations for each pixel position: Determine the first and second adjacent pixels corresponding to the pixel position; The first coherence coefficient corresponding to the first adjacent pixel and the second coherence coefficient corresponding to the second adjacent pixel are obtained from the coherence coefficient map. The minimum value between the first coherence coefficient and the second coherence coefficient is determined as the weight value corresponding to the pixel position in the weight field corresponding to the current gradient extraction direction; After the traversal is completed, the weight field of the differential interferogram in the current gradient extraction direction is obtained.

4. The method according to claim 3, characterized in that, Determining the average deformation gradient field of the target study area in the gradient extraction direction includes: For each gradient extraction direction, the deformation gradient fields of all the differential interferograms in the gradient extraction direction are multiplied pixel by pixel with their respective weight fields, and the multiplication results of all the differential interferograms are accumulated pixel by pixel to obtain the gradient accumulation field. The weight fields of all the differential interferograms in the gradient extraction direction are accumulated pixel by pixel to obtain the weight accumulation field; Based on the gradient accumulation field and the weight accumulation field, the average deformation gradient field of the target study area in the gradient extraction direction is determined.

5. The method according to claim 1, characterized in that, The process of reconstructing the deformation rate field using the Poisson equation based on the average deformation gradient field along each gradient extraction direction includes: The average deformation gradient field in the first target direction and the average deformation gradient field in the second target direction are obtained from the average deformation gradient fields in each of the gradient extraction directions; wherein the first target direction and the second target direction are perpendicular to each other; A first partial derivative field is obtained by taking the partial derivative of the average deformation gradient field in the first target direction in the first target direction; and a second partial derivative field is obtained by taking the partial derivative of the average deformation gradient field in the second target direction in the second target direction. Based on the first and second partial derivative fields, the target divergence field is determined; The Poisson equation is constructed by using the target divergence field as the source term of the Poisson equation and the deformation rate field to be reconstructed as the unknown variable to be solved in the Poisson equation. The Poisson equation is solved based on preset boundary conditions to obtain the deformation rate field of the target study area.

6. The method according to claim 5, characterized in that, Solving the Poisson equation based on preset boundary conditions includes: The Poisson equation is transformed from the spatial domain to the frequency domain, and the Poisson equation is solved in the frequency domain based on the preset boundary conditions to obtain the frequency domain deformation coefficients. The frequency domain deformation coefficients are inversely transformed from the frequency domain back to the spatial domain to obtain the deformation rate field of the target study area.

7. The method according to claim 1, characterized in that, The determination of the deformation estimation result of the target study area based on the deformation rate field of the target study area includes: The deformation rate value of each pixel in the deformation rate field of the target study area is taken as the deformation rate of the pixel. Based on the deformation rate of all the pixels, identify the deformation region in the target study area, the boundary of the deformation region, and the deformation rate distribution inside the deformation region; The identified deformation regions, the boundaries of the deformation regions, and the deformation rate distribution within the deformation regions are used as the deformation estimation results for the target study area.

8. A landslide deformation estimation device, characterized in that, include: The data acquisition module is used to acquire multiple differential interferograms of the target study area and a coherence coefficient map corresponding to each differential interferogram; The gradient and weight field determination module is used to determine the deformation gradient field of each differential interferogram in multiple preset gradient extraction directions, and to determine the weight field corresponding to each deformation gradient field based on the coherence coefficient map corresponding to the differential interferogram. The average gradient field determination module is used to determine the average deformation gradient field of the target study area in the gradient extraction direction based on the deformation gradient field and the corresponding weight field of all the differential interferograms in the same gradient extraction direction. The deformation estimation module is used to reconstruct the deformation rate field based on the average deformation gradient field in each of the gradient extraction directions using the Poisson equation, thereby obtaining the deformation rate field of the target study area; and to determine the deformation estimation result of the target study area based on the deformation rate field of the target study area.

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

10. A computer device, comprising a storage medium, a processor, and a computer program stored on the storage medium and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 7.

Images