High-precision insar phase gradient rate calculation method, device and medium

Abstract

The present disclosure relates to the technical field of ground large gradient deformation monitoring, and provides a high-precision InSAR phase gradient rate calculation method, device and medium, comprising: acquiring SAR data set and DEM data of a region to be monitored and generating differential interferograms; mapping the wrapped phase values corresponding to each pixel point in each differential interferogram into complex exponential signals, determining the real part and imaginary part data matrixes of each differential interferogram based on the real part and imaginary part data of each wrapped phase value extracted therefrom; determining the spatial domain directional gradient field in multiple directions of each graph based on the real part and imaginary part data matrixes corresponding to each differential interferogram and the Riesz-Gauss transfer function family; and performing stack average processing on the spatial domain directional gradient field of all differential interferograms in the same direction to determine the phase gradient rate graph of the region to be monitored. The present disclosure enhances the robustness of phase gradient rate calculation of the ground large gradient deformation region and improves the ground large gradient deformation monitoring accuracy.

Description

Technical Field

[0001] This disclosure relates to the field of surface large gradient deformation monitoring technology, and more specifically, to a high-precision InSAR phase gradient rate calculation method, device and medium. Background Technology

[0002] Synthetic Aperture Radar Interferometry (InSAR) technology has been widely used in the field of surface deformation monitoring, and its derived time-series analysis methods (such as PS-InSAR and SBAS-InSAR) can achieve millimeter-level deformation detection. However, traditional deformation inversion mostly relies on phase unwrapping, which is prone to unwrapping errors in low-coherence regions and propagates along the integration path, affecting the final accuracy. Phase Gradient Rate (PGR) provides an alternative technique that does not require phase unwrapping: by directly calculating the spatial phase gradient amplitude of the wrapped interferogram and performing time-dimensional stacked averaging, it can quickly identify deformation abrupt locations such as landslide edges, fault traces, and subsidence basin boundaries, showing unique advantages in large-scale geological hazard surveys. The computational quality of PGR is highly dependent on the noise propagation characteristics of the gradient operator, which directly determines the signal-to-noise ratio and reliability of the final detection results.

[0003] In related technologies, spatial domain gradient operators are commonly used for PGR calculation, mainly including the central difference method and the Sobel operator. The central difference method approximates the gradient by using the entangled phase difference between adjacent pixels. Its frequency domain transfer function exhibits high-pass characteristics, inevitably amplifying high-frequency noise components while extracting the gradient signal, which is particularly detrimental to low-coherence regions. Although the Sobel operator introduces a [1,2,1] weighted average perpendicular to the gradient direction to suppress some noise, its smoothness is determined by a fixed step size parameter, which cannot be continuously adjusted according to the data noise level, and the boundary diffusion effect becomes significant as the step size increases. More importantly, both are finite difference templates with fixed directions, resulting in inconsistent responses to deformation gradients with different orientations, leading to directional bias. These limitations cause related PGR methods to face insufficient signal-to-noise ratio in low-coherence regions such as vegetation cover and complex terrain, limiting their deformation boundary detection capabilities. Summary of the Invention

[0004] This disclosure provides at least one high-precision InSAR phase gradient rate calculation method, apparatus, and medium, which enhances the robustness of phase gradient rate calculation in large surface gradient deformation regions and improves the accuracy of large surface gradient deformation monitoring.

[0005] This disclosure provides a high-precision InSAR phase gradient rate calculation method, including: Acquire SAR datasets and DEM data for the area to be monitored; and generate differential interferograms based on the SAR datasets and DEM data; For each differential interferogram, the winding phase value corresponding to each pixel in the differential interferogram is mapped to a complex exponential signal; and based on the real part data and imaginary part data of the winding phase value corresponding to each pixel extracted from the complex exponential signal, the real part data matrix and imaginary part data matrix of the differential interferogram are determined. Based on the real part data matrix and the imaginary part data matrix corresponding to each difference interferogram, and the pre-constructed family of Riesz-Gauss transfer functions, the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions is determined; wherein, the pre-constructed family of Riesz-Gauss transfer functions includes Riesz-Gauss transfer functions corresponding to the multiple gradient extraction directions respectively. The spatial domain directional gradient fields of all differential interferograms in the same gradient extraction direction are stacked and averaged in the time dimension, and the phase gradient rate map of the region to be monitored is determined based on the stacked and averaged results of multiple gradient extraction directions.

[0006] This disclosure provides a high-precision InSAR phase gradient rate calculation device, comprising: The data acquisition module is used to acquire SAR datasets and DEM data of the area to be monitored; and to generate differential interferograms based on the SAR datasets and DEM data. The data mapping module is used to map the entangled phase value corresponding to each pixel in the differential interferogram to a complex exponential signal for each differential interferogram, and to determine the real part data matrix and the imaginary part data matrix of the differential interferogram based on the real part data and imaginary part data of the entangled phase value corresponding to each pixel extracted from the complex exponential signal. The gradient field determination module is used to determine the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions based on the real part data matrix and the imaginary part data matrix corresponding to each difference interferogram, as well as the pre-constructed Riesz-Gauss transfer function family; wherein, the pre-constructed Riesz-Gauss transfer function family includes Riesz-Gauss transfer functions corresponding to the multiple gradient extraction directions respectively. The data stacking module is used to perform time-dimensional stacked averaging of the spatial domain directional gradient fields of all differential interferograms in the same gradient extraction direction, and to determine the phase gradient rate map of the region to be monitored based on the stacked averaging results of multiple gradient extraction directions.

[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 high-precision InSAR phase gradient rate calculation 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 high-precision InSAR phase gradient rate calculation method as described in any of the possible embodiments above.

[0009] The high-precision InSAR phase gradient rate calculation method, apparatus, and medium provided in this disclosure transform the originally discontinuous periodic phase angle into a continuous complex domain signal by mapping the wrapped phase value to a complex exponential signal and extracting the real and imaginary parts. This allows the gradient operator to be executed directly in the complex domain without triggering the ambiguity caused by the 2π jump. Based on this, a pre-constructed family of Riesz-Gauss transfer functions is used to simultaneously extract gradients in multiple directions and attenuate high-frequency noise in the frequency domain. The continuity of the complex domain ensures that the gradient operation can obtain stable results without unwrapping, thereby fundamentally avoiding error propagation caused by unwrapping errors. Furthermore, the spatial domain directional gradient fields in the same direction of all differential interferograms are subjected to time-dimensional stacking and averaging. Based on the redundant observation of the long-term series, the weight of random noise and local outliers in a single image is reduced, while consistent deformation gradient information is preserved. This achieves direct estimation of the phase gradient rate without phase unwrapping, improving the robustness and accuracy of monitoring results in large-gradient deformation areas of the Earth's surface.

[0010] 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

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

[0012] Figure 1A flowchart of a high-precision InSAR phase gradient rate calculation method provided by an embodiment of this disclosure is shown; Figure 2 A flowchart of a differential interferogram generation method provided in an embodiment of this disclosure is shown; Figure 3 A flowchart of a method for determining a spatial domain directional gradient field provided by an embodiment of this disclosure is shown; Figure 4 A flowchart of a spatial domain oriented gradient field stacking method provided in an embodiment of this disclosure is shown; Figure 5 A comparative schematic diagram of phase gradient rate results obtained by different gradient operators provided in the embodiments of this disclosure is shown, wherein... Figure 5 (a) is a schematic diagram of the location of the study area provided in an embodiment of this disclosure. Figure 5 (b) is a schematic diagram of the phase gradient results of the central difference method provided in the embodiments of this disclosure. Figure 5 (c) is a schematic diagram of the phase gradient result of the Sobel operator provided in the embodiments of this disclosure. Figure 5 (d) is a schematic diagram of the phase gradient result of the Riesz-Gauss operator provided in the embodiments of this disclosure; Figure 6 A comparative schematic diagram of the PGR calculation results of three gradient operators in a first typical landslide area provided in the embodiments of this disclosure is shown, wherein... Figure 6 (a) is an optical remote sensing image of the first typical landslide area and a schematic diagram of the landslide boundary provided in the embodiments of this disclosure. Figure 6 (b) is a schematic diagram of the central difference method normalized PGR results for the first typical landslide area provided in the embodiments of this disclosure. Figure 6 (c) is a schematic diagram of the Sobel operator-normalized PGR results for the first typical landslide area provided in the embodiments of this disclosure. Figure 6 (d) is a schematic diagram of the Riesz-Gauss operator normalized PGR results for the first typical landslide area provided in the embodiments of this disclosure; Figure 7 A comparative schematic diagram of the PGR calculation results of three gradient operators in a second typical landslide area provided in the embodiments of this disclosure is shown, wherein... Figure 7 (a) is an optical remote sensing image of a second typical landslide area and a schematic diagram of the landslide boundary provided in an embodiment of this disclosure. Figure 7 (b) is a schematic diagram of the central difference method normalized PGR results for a second typical landslide area provided in the embodiments of this disclosure. Figure 7 (c) is a schematic diagram of the Sobel operator-normalized PGR results for a second typical landslide area provided in the embodiments of this disclosure. Figure 7(d) is a schematic diagram of the Riesz-Gauss operator normalized PGR results for a second typical landslide area provided in the embodiments of this disclosure; Figure 8 This illustration shows a comparison diagram of the PGR calculation results of three gradient operators in a third typical landslide area provided in this embodiment of the present disclosure, wherein... Figure 8 (a) is an optical remote sensing image of a third typical landslide area and a schematic diagram of the landslide boundary provided in an embodiment of this disclosure. Figure 8 (b) is a schematic diagram of the central difference method normalized PGR results for a third typical landslide area provided in the embodiments of this disclosure. Figure 8 (c) is a schematic diagram of the Sobel operator-normalized PGR results for a third typical landslide area provided in this embodiment of the present disclosure. Figure 8 (d) is a schematic diagram of the Riesz-Gauss operator normalized PGR results for a third typical landslide area provided in the embodiments of this disclosure; Figure 9 A schematic diagram of the structure of a high-precision InSAR phase gradient rate calculation device provided in an embodiment of this disclosure is shown. Figure 10 A schematic diagram of the structure of a computer device provided in an embodiment of this disclosure is shown. Detailed Implementation

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

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

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

[0016] Interferometric Synthetic Aperture Radar (InSAR) technology, with its all-weather, wide-area, and high-precision surface deformation monitoring capabilities, has become an important technical means for studying geological hazards such as landslides, ground subsidence, and coseismic deformation caused by earthquakes. Temporal InSAR methods (such as PS-InSAR and SBAS-InSAR) can obtain millimeter-level precision time series of surface deformation by jointly processing multiple SAR images, providing crucial data support for the early identification and dynamic monitoring of wide-area geological hazards.

[0017] However, extracting deformation information from the interference phase usually requires a crucial step of phase unwrapping, i.e., the affected... The fuzzy-constrained entangled phase is restored to a continuous absolute phase. This process is prone to unwrapping errors in low-coherence regions, and these errors accumulate and propagate along the integration path, thus affecting the deformation inversion accuracy of the entire interferogram. In addition, phase unwrapping has a large computational cost, limiting its processing efficiency in applications involving large scales and long time series.

[0018] Phase Gradient Rate (PGR) provides a technique for deformation detection without phase unwrapping. PGR is defined as the time average of the spatial phase gradient amplitude across interferograms, calculated directly from the wrapped phase, thus avoiding the phase unwrapping step. Unlike the deformation phase, which characterizes cumulative displacement, PGR depicts locations where significant spatial abrupt changes in deformation rate occur, such as landslide edges, fault traces, and subsidence basin boundaries. This rapid response characteristic gives PGR a unique advantage in large-scale geological hazard surveys.

[0019] The calculation of PGR consists of two stages: first, the spatial phase gradient field is extracted from each wound interferogram using a gradient operator; then, the gradient magnitudes of each image are stacked and averaged over time. This time stacking allows for... Proportional suppression of random noise (where The noise level of the gradient field of a single interferogram depends on the gradient operator used. Therefore, the noise transmission characteristics of the gradient operator directly determine the signal-to-noise ratio and detection quality of the PGR.

[0020] Research has found that current PGR calculations commonly employ spatial domain gradient operators, primarily the central difference method and the Sobel operator. The central difference method approximates the gradient by using the entangled phase difference between adjacent pixels; the Sobel operator, building upon this, introduces a weighted smoothing kernel and uses a 3×3 convolution window to achieve a certain degree of noise suppression. These methods are computationally efficient and simple to implement, but they have inherent limitations: the numerical differentiation of noisy data is a classical ill-posed problem, and the finite difference method inevitably amplifies high-frequency noise components while extracting the gradient signal. Specifically, the frequency domain transfer function of the central difference method is... The phase gradient calculation method exhibits high-pass characteristics, increasing approximately linearly with frequency in the range of 0 to 0.25 cycles / pixel and peaking at f = 0.25, but its ability to suppress high-frequency noise is insufficient. While the Sobel operator introduces a [1,2,1] weighted average perpendicular to the gradient extraction direction to achieve partial noise suppression, its smoothness is fixed by the step size parameter and cannot be continuously adjusted; furthermore, the larger the step size, the more severe the boundary diffusion. In addition, both methods use finite difference templates with fixed orientations, resulting in inconsistent responses to deformation gradients with different orientations and exhibiting directional bias. Therefore, existing phase gradient calculation methods face insufficient signal-to-noise ratio in low-coherence regions, limiting the detection performance of PGR in complex surface environments.

[0021] Based on the above research, this disclosure provides a high-precision InSAR phase gradient rate calculation method, device, and medium. By mapping the wrapped phase value to a complex exponential signal and extracting the real and imaginary parts, the originally discontinuous periodic phase angle is transformed into a continuous complex domain signal, thus allowing the gradient operator to be executed directly in the complex domain without triggering the ambiguity caused by the 2π jump. On this basis, a pre-constructed family of Riesz-Gauss transfer functions is used to simultaneously achieve gradient extraction and high-frequency noise attenuation in multiple directions in the frequency domain. The continuity of the complex domain ensures that the gradient operation can obtain stable results without unwrapping, thereby fundamentally avoiding error propagation caused by unwrapping errors. Furthermore, the spatial domain directional gradient fields in the same direction of all differential interferograms are subjected to time-dimensional stacking and averaging. Based on the redundant observation of the long-term series, the weight of random noise and local outliers in a single image is reduced, while retaining consistent deformation gradient information. This achieves direct estimation of the phase gradient rate without phase unwrapping, improving the robustness and accuracy of monitoring results in large-gradient deformation areas of the Earth's surface.

[0022] To facilitate understanding of this embodiment, the executing entity of the high-precision InSAR phase gradient rate calculation method provided in this disclosure will first be described in detail. The executing entity of the high-precision InSAR phase gradient rate calculation method provided in this disclosure 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.

[0023] The high-precision InSAR phase gradient rate calculation method provided in this application embodiment will be described in detail below with reference to the accompanying drawings. See also Figure 1 The diagram shows a flowchart of a high-precision InSAR phase gradient rate calculation method provided in this embodiment of the present disclosure. The method includes the following steps S101 to S104: S101, acquire the SAR dataset and DEM data of the area to be monitored; and generate a differential interferogram based on the SAR dataset and DEM data.

[0024] Understandably, SAR (Synthetic Aperture Radar) is a technology that images the ground using electromagnetic waves. It can acquire surface information around the clock and in all weather conditions, and is widely used in areas such as surface change monitoring and disaster early warning. SAR datasets are image data collected by radar sensors that contain electromagnetic wave reflection information of the monitored area, providing accurate information on surface topography and deformation. DEM (Digital Elevation Model) is a digital model describing the distribution of ground elevation obtained through remote sensing or ground surveying. DEM data can be used to aid in understanding terrain undulations and to remove phase components caused by terrain in subsequent interferometric processing, thereby highlighting surface deformation signals. It can also be used to correct geometric distortions caused by terrain during registration, improving the accuracy of image alignment.

[0025] In SAR interferometry, differential interferograms are images that reflect surface deformation obtained by calculating the phase difference between two SAR images. They can show the changes that occur on the surface between two points in time and are often used to monitor phenomena such as landslide creep, ground subsidence, and coseismic deformation.

[0026] Specifically, refer to Figure 2As shown, determining the differential interferogram based on the SAR dataset and DEM data may include the following steps S201~S204: S201, Determine the initial SAR image set based on the SAR dataset.

[0027] The initial SAR image set includes multiple SAR images of the area to be monitored. The area to be monitored usually refers to a specific area that needs to be monitored for ground changes, such as high-risk areas like landslide-prone areas, mining subsidence areas, or urban land subsidence areas.

[0028] S202, determine the SAR main image in the initial SAR image set, and based on the SAR main image, register other SAR images in the initial SAR image set except for the SAR main image according to the DEM data, and determine the target SAR image set based on the registration results.

[0029] Specifically, the SAR master image refers to the reference image selected from the entire image set, typically the one closest in time or of the highest quality. Registration involves aligning other images with the master image, ensuring precise spatial matching of SAR images from different times and angles. During registration, DEM data is used to correct geometric errors in the images, thereby ensuring spatial consistency. Finally, these registered SAR images and the SAR master image are used to generate a target SAR image set for subsequent analysis. The target SAR image set includes multiple registered target SAR images.

[0030] S203, determine the baseline set based on the target SAR image set according to the preset spatiotemporal baseline threshold.

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

[0032] S204, Generate a differential interferogram corresponding to each interference pair based on the DEM data.

[0033] Specifically, an interferometric pair is formed from two radar observations (or images), each with slightly different times and locations. By comparing these two images, a differential interferogram reflecting surface changes (such as subsidence and displacement) can be obtained. The differential interferogram is obtained by subtracting topographic errors (usually corrected using DEM data) and can display fine information about relative surface changes. The differential interferogram of each interferometric pair reveals the phase difference between the two observations, reflecting information about ground displacement or deformation.

[0034] Here, after the differential interferometry processing in steps S201~S204 above, the entanglement phase in the differential interferogram corresponding to each interference pair can be expressed as: ; in, Represented as the winding phase value at spatial location x, this value is obtained through the winding operator. Limit the true phase principal value to Within the range; It is represented as the true phase component caused by surface deformation, which is proportional to the displacement of the surface along the radar line of sight. It is represented as the residual systematic error phase, which includes the non-deformation component composed of the residual terrain phase, atmospheric delay phase, and orbital error phase that were not completely eliminated after external DEM correction; The phase is represented as the decoherent noise phase, which originates from factors such as thermal noise, registration error, and time decoherence. Its statistical characteristics follow a complex circular Gaussian distribution, and its variance is determined by the Cramér-Rao lower bound. Represented as the entanglement operator, its function is to map the true phase value modulo 2π within the parentheses to the principal value interval. within, that is , where k is an integer that makes the result fall within the principal value interval.

[0035] In this disclosure, to reduce the risk of phase gradients exceeding π radians due to deformation accumulation in long-term baseline interferograms, interferometric pairs with the shortest time baselines are preferentially selected for subsequent phase gradient rate calculations. For example, among interferometric pairs with time baseline entries of 12 days, 24 days, and 36 days, the 12-day time baseline interferometric pair is preferentially selected. This ensures that the true phase difference between adjacent pixels generally does not exceed π, thus satisfying the phase continuity assumption and improving the accuracy of subsequent phase gradient estimation.

[0036] S102, for each of the differential interferograms, the winding phase value corresponding to each pixel in the differential interferogram is mapped to a complex exponential signal; and based on the real part data and imaginary part data of the winding phase value corresponding to each pixel extracted from the complex exponential signal, the real part data matrix and imaginary part data matrix of the differential interferogram are determined.

[0037] Here, after obtaining each differential interferogram through differential interferometry processing, each pixel in the image corresponds to a wrapped phase value. This wrapped phase value is obtained by performing an arctangent operation on the interferometric phase and taking the principal value. Its numerical range is limited to between negative and positive pi. The wrapped phase value refers to the remaining part after the modulus of the interferometric phase is 2π. It cannot directly reflect integer phase changes beyond this range, but it contains comprehensive information on surface deformation, residual errors, and noise.

[0038] Understandably, the wrapped phase value in the differential interferogram is confined to the principal value range of negative π to positive π, while the actual surface deformation phase may exceed this range, potentially causing a 2π jump between adjacent pixels. This phase jump manifests as a spatial discontinuity. If the wrapped phase is directly numerically differentiated or frequency-domain filtered, spurious high-frequency components will be generated at the jump point, interfering with the extraction of the true gradient signal. To circumvent the discontinuity of the wrapped phase, the wrapped phase value of each pixel can be mapped to a complex exponential signal. The complex exponential signal is defined by Euler's formula as a function with the natural constant as the base and the imaginary unit multiplied by the wrapped phase as the exponent, resulting in a complex number with a constant modulus of 1. The argument of this complex signal is equal to the original wrapped phase value, and it itself is a point on a continuously varying unit circle in the complex plane. When the wrapped phase jumps from π to -π, the complex exponential signal does not jump because π and -π correspond to the same complex point, -1. Therefore, complex exponential signals are spatially continuous and differentiable, making them suitable as inputs for subsequent frequency domain gradient estimation.

[0039] For example, the complex exponential signal at each pixel can be represented by the following formula: ; in, Represented as a complex exponential signal located at spatial location x, it is a complex number; Let represent the imaginary unit, satisfying ; Represented as a base of the natural constant e and an exponent of 1 / e. The exponential function, according to Euler's formula, is equal to ; Thus, by introducing a complex exponential signal, the originally discontinuous winding phase gradient can be transformed into a combination of derivatives of a continuous function, specifically expressed as: ; in, Represented as entangled phase The sine value, i.e., the imaginary part of the complex exponential signal; Represented as entangled phase The cosine value is the real part of the complex exponential signal. This formula transforms the calculation of the phase gradient into a linear combination of the derivatives of the continuously differentiable cosine and sine fields, thus avoiding the problem of directly dealing with the entangled phase jump.

[0040] Furthermore, after obtaining the complex exponential signal for each pixel, the real and imaginary data can be extracted from it. The real part of the complex exponential signal is equal to the cosine of the wrapped phase, and the imaginary part is equal to the sine of the wrapped phase; both components are real numbers and change continuously throughout the image space. Based on the extracted real data of each pixel, arranged according to the row and column positions of the pixels, the real data matrix of the differential interferogram can be determined. The real data matrix is ​​a two-dimensional array whose row and column numbers correspond to the row and column numbers of the differential interferogram, respectively, and the value of each element represents the cosine of the wrapped phase value at the corresponding pixel position. Similarly, based on the imaginary data of each pixel, arranged in the same row and column order, the imaginary data matrix can be obtained. The imaginary data matrix represents a two-dimensional array with the same size as the real data matrix, and the value of each element represents the sine of the wrapped phase value at the corresponding pixel position.

[0041] S103, based on the real part data matrix and the imaginary part data matrix corresponding to each difference interferogram, and the pre-constructed Riesz-Gauss transfer function family, determine the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions.

[0042] Here, the pre-constructed family of Riesz-Gauss transfer functions is a collection of direction-specific filters, including Riesz-Gauss transfer functions corresponding to multiple gradient extraction directions. A Riesz-Gauss transfer function is a complex filter defined in the frequency domain, consisting of the Riesz transform component multiplied by a Gaussian low-pass filter component. It is primarily used to simultaneously achieve gradient extraction and high-frequency noise attenuation in the frequency domain. Specifically, gradient extraction is the process of calculating the rate of change of phase with space at each pixel location in an image, quantifying the rate of change of surface deformation in space. High-frequency noise attenuation refers to the process of suppressing high-frequency noise components in the image, reducing the interference of random noise on the gradient estimation results.

[0043] Understandably, since the real and imaginary data matrices corresponding to each differential interferogram are continuously differentiable two-dimensional signals, they can be efficiently processed in the frequency domain. Furthermore, the pre-constructed Riesz-Gauss transfer function possesses isotropic gradient extraction and frequency-adaptive noise suppression capabilities. Therefore, based on the real and imaginary data matrices corresponding to each differential interferogram, and the pre-constructed Riesz-Gauss transfer function, the phase gradient in each direction can be calculated through frequency domain multiplication and inverse transformation to determine the spatial domain directional gradient field of each differential interferogram in different directions. Here, the spatial domain directional gradient field is represented as a two-dimensional array with the same size as the original differential interferogram, where the value of each element represents the phase change rate along a specific direction at that pixel location, reflecting the intensity of the surface deformation boundary in that direction.

[0044] In some possible embodiments, when determining the spatial domain orientation gradient field of each differential interferogram, reference is made to... Figure 3 As shown, for each differential interferogram, the following steps S301~S304 can be performed: S301, perform a frequency domain transformation on the real part data matrix of the differential interferogram to obtain the frequency domain real part data matrix corresponding to the differential interferogram; and perform a frequency domain transformation on the imaginary part data matrix of the differential interferogram to obtain the frequency domain imaginary part data matrix corresponding to the differential interferogram.

[0045] It is understandable that convolution operations in the spatial domain are transformed into element-wise multiplication operations in the frequency domain, which is more computationally efficient. Therefore, we can first perform frequency domain transformation on the real and imaginary data matrices to convert the spatial domain signal to the frequency domain. The Fast Fourier Transform algorithm can be used to significantly reduce the computational complexity.

[0046] Here, the specific method for performing frequency domain transformation on the real and imaginary data matrices can be selected according to the actual data size and processing platform, such as two-dimensional fast Fourier transform, two-dimensional discrete cosine transform, or two-dimensional wavelet transform. The embodiments of this disclosure use two-dimensional fast Fourier transform to obtain the optimal computational efficiency, and no other feasible methods are specifically limited here.

[0047] S302, for each of the preset multiple gradient extraction directions, the frequency domain real part data matrix and frequency domain imaginary part data matrix corresponding to the differential interferogram are multiplied element-wise in the frequency domain with the Riesz-Gauss transfer function corresponding to the direction to obtain the frequency domain real part gradient data matrix and frequency domain imaginary part gradient data matrix in the direction.

[0048] Furthermore, after obtaining the frequency domain real part data matrix and the frequency domain imaginary part data matrix, for each preset gradient extraction direction, the frequency domain real part data matrix can be multiplied element-wise with the Riesz-Gauss transfer function corresponding to the direction, and the frequency domain imaginary part data matrix can be multiplied element-wise with the same transfer function in the frequency domain to determine the frequency domain real part gradient data matrix and the frequency domain imaginary part gradient data matrix in that direction.

[0049] Among them, the preset multiple gradient extraction directions refer to the set of spatial directions selected in advance for calculating the phase gradient. These directions can be preset according to the anisotropic characteristics of the deformation field, such as 0-degree direction, 45-degree direction, 90-degree direction and 135-degree direction, or three directions of 0-degree, 60-degree and 120-degree. No specific limitation is made here.

[0050] For example, the above-mentioned family of Riesz-Gauss transfer functions can be pre-constructed through the following steps (1) to (5): (1) Construct the Riesz transform transfer function. The amplitude-frequency response of the Riesz transform transfer function along any gradient extraction direction is constant, which is used to realize isotropic gradient extraction. (2) Construct a Gaussian low-pass transfer function. The amplitude-frequency response of the Gaussian low-pass transfer function decays exponentially with increasing frequency, which is used to suppress high-frequency noise components. (3) Multiply the Riesz transform transfer function and the Gaussian low-pass transfer function in the frequency domain to obtain the basic Riesz-Gauss transfer function; (4) For each preset gradient extraction direction, project the basic Riesz-Gauss transfer function along the gradient extraction direction to obtain the Riesz-Gauss transfer function corresponding to the gradient extraction direction. (5) Construct a family of Riesz-Gauss transfer functions based on the Riesz-Gauss transfer functions corresponding to all gradient extraction directions.

[0051] Here, the Riesz transform transfer function is a linear operator in the frequency domain, a generalization of the two-dimensional Hilbert transform. In the frequency domain, the Riesz transform transfer function component along direction θ is negative j multiplied by the direction cosine and divided by the radial frequency; its amplitude is always 1 and does not change with frequency. The Riesz transform transfer function can extract the partial derivatives of an image along any direction while maintaining the signal's energy distribution, thus achieving an unbiased response to the gradient direction. The Gaussian low-pass transfer function in the frequency domain is a function with the natural constant e as the base and a negative exponent, divided by twice the square of the frequency domain Gaussian width. This function can be regarded as an exponentially decaying function with the square of the radial frequency as the independent variable; as the frequency increases, the transfer coefficient decreases exponentially, thus effectively suppressing high-frequency noise components.

[0052] Furthermore, after constructing the basic Riesz-Gauss transfer function, the Riesz transform transfer function and the Gaussian low-pass transfer function can be multiplied point-by-point in the frequency domain to obtain a gradient extraction function with low-pass characteristics. The basic Riesz-Gauss transfer function can preserve gradient information in the low-frequency band and attenuate noise in the high-frequency band, achieving a balance between gradient extraction and noise suppression.

[0053] It is understandable that transfer functions corresponding to different gradient extraction directions have different directional projection factors. Therefore, to obtain gradient fields in multiple directions, the basic Riesz-Gauss transfer function can be directionally projected along different gradient extraction directions. That is, the direction variable in the original transfer function is replaced with the direction cosine of that direction, resulting in a Riesz-Gauss transfer function corresponding to each gradient extraction direction. Then, by combining the Riesz-Gauss transfer functions corresponding to all gradient extraction directions, a family of Riesz-Gauss transfer functions can be obtained. This family of functions covers the required multiple gradient extraction directions, with each member function responsible for gradient extraction and noise suppression in its corresponding direction.

[0054] For example, the Riesz-Gauss transfer functions for each direction can be expressed as: ; in, It is represented as the gradient extraction direction angle, in radians or degrees, and is used to specify the direction along which the phase gradient is calculated. Represented as a spatial domain Gaussian smoothing parameter, in pixels, it controls the cutoff frequency of the Gaussian low-pass filter. The larger the value, the stronger the high-frequency attenuation; and These represent the frequency coordinates in the horizontal and vertical directions in the frequency domain, respectively. This is expressed as given smoothing parameters and direction The Riesz-Gauss transfer function in frequency coordinates Complex values ​​at; Represented as radial frequency; Represented as the Gaussian width in the frequency domain, and the smoothing parameter in the spatial domain. Inversely proportional, The units are consistent with the frequency coordinates; It is expressed as the square of the Gaussian width in the frequency domain, used to normalize the square of the frequency in the exponential term; Represented as the Riesz transform along the direction The term multiplied by the projection and the radial frequency is used to extract the directional gradient; Represented as a Gaussian low-pass term, where This determines the decay rate. After multiplying the two, the amplitude-frequency response of the entire transfer function is dominated by Gaussian terms in the low-frequency range, thus achieving the extraction of gradients while suppressing high-frequency noise.

[0055] Here, when At that time, the directional projection can degenerate into That is, the standard x-direction Riesz transform.

[0056] Thus, and Substitute into the above equation, along The phase gradient Riesz-Gauss estimate in the direction can be expressed as: ; in, It is represented as the corresponding pixel value in the imaginary part data matrix, that is, the sine value of the wrapped phase of that pixel; This is represented as applying the Riesz-Gauss operator to the real part of the data matrix. The following obtained along The value of the spatial real part gradient field of the orientation at that pixel point; It is represented as the corresponding pixel value in the real part data matrix, that is, the cosine value of the wrapped phase of that pixel; This is represented as applying the Riesz-Gauss operator to the imaginary part of the data matrix. The following obtained along The value of the spatial imaginary part gradient field of the direction at this pixel.

[0057] Specifically, the process of performing element-wise multiplication of the real part data matrix and the imaginary part data matrix in the frequency domain with the Riesz-Gauss transfer function corresponding to the direction can include the following steps (a)~(b): (a) For each gradient extraction direction, obtain the pre-constructed Riesz-Gauss transfer function corresponding to the gradient extraction direction; wherein the Riesz-Gauss transfer function is represented in the frequency domain as a numerical matrix with the same size as the real part data matrix in the frequency domain; (b) Multiply the real part data matrix in the frequency domain with the elements at the same position in the numerical matrix to obtain the real part gradient data matrix in the frequency domain; and multiply the imaginary part data matrix in the frequency domain with the elements at the same position in the numerical matrix to obtain the imaginary part gradient data matrix in the frequency domain.

[0058] Understandably, after pre-constructing the Riesz-Gauss transfer function corresponding to each gradient extraction direction based on the above steps, when processing each direction of each differential interferogram, the Riesz-Gauss transfer function corresponding to that direction can be obtained first. This transfer function is represented in the frequency domain as a numerical matrix with the same size as the real part data matrix in the frequency domain. Each element of this numerical matrix is ​​a complex number, its magnitude determined by the Gaussian low-pass transfer function, and its phase determined by the direction projection term of the Riesz transform. This allows for the application of different weights to each frequency component based on its frequency level and the gradient extraction direction. By element-wise multiplying this numerical matrix with the real part data matrix in the frequency domain, frequency domain filtering of the real part data can be performed to obtain the real part gradient data matrix in the frequency domain; similarly, element-wise multiplying the same numerical matrix with the imaginary part data matrix in the frequency domain yields the imaginary part gradient data matrix in the frequency domain.

[0059] In this way, by using the element-wise multiplication operation in the frequency domain and the pre-constructed direction-dependent Riesz-Gauss transfer function, gradient extraction and noise suppression can be completed simultaneously in a single frequency domain multiplication, achieving efficient gradient estimation of the continuous components of the entangled phase in different directions.

[0060] S303, perform spatial domain transformation on the frequency domain real part gradient data matrix in each direction to obtain the spatial domain real part gradient field in the direction; perform spatial domain transformation on the frequency domain imaginary part gradient data matrix in each direction to obtain the spatial domain imaginary part gradient field in the direction.

[0061] Specifically, after calculating the real part gradient data matrix and the imaginary part gradient data matrix in the frequency domain for each direction, they can be transformed back from the frequency domain to the spatial domain to obtain a spatial domain gradient field that can be used for subsequent synthesis.

[0062] Here, the specific method for spatial domain transformation of the real part gradient data matrix and the imaginary part gradient data matrix in the frequency domain can be selected according to computational efficiency and hardware platform, such as two-dimensional fast Fourier transform, two-dimensional discrete cosine transform or two-dimensional wavelet transform, etc. The embodiments of this disclosure adopt two-dimensional fast Fourier transform to ensure consistency with frequency domain transformation.

[0063] S304. Based on the spatial domain real part gradient field and spatial domain imaginary part gradient field in each direction, the spatial domain directional gradient field in the stated direction is synthesized.

[0064] Furthermore, after performing the inverse transformation of the frequency domain results, the spatial domain real part gradient field and spatial domain imaginary part gradient field in the same direction can be synthesized to obtain the spatial domain directional gradient field in that direction. Here, the synthesis operation can be performed based on the relationship between the phase gradient and the cosine and sine partial derivatives. That is, at each pixel position, the negative imaginary part value is multiplied by the real part gradient value, and then the real part value is multiplied by the imaginary part gradient value to obtain the phase gradient value of that pixel. Repeating the above synthesis operation for all directions yields the spatial domain directional gradient fields in multiple directions for each differential interferogram. These directional gradient fields reflect the spatial rate of change of surface deformation in different directions.

[0065] S104, perform time-dimensional stacked averaging on the spatial domain gradient fields of all differential interferograms in the same gradient extraction direction, and determine the phase gradient rate map of the region to be monitored based on the stacked averaging results of multiple gradient extraction directions.

[0066] Understandably, after obtaining the spatial domain gradient fields of all differential interferograms in each gradient extraction direction, a time-dimensional stacked averaging process can be performed to compress random noise and enhance stable deformation boundary signals. Here, the time-dimensional stacked averaging process refers to accumulating the gradient values ​​of all differential interferograms at the same spatial location and in the same gradient extraction direction, and then dividing by the total number of differential interferograms. This can reduce noise levels while preserving consistent gradient information across different time phases.

[0067] Specifically, in order to obtain robust stacked averaging results, when performing stacked averaging in the time dimension, refer to Figure 4 As shown, the steps S401 to S402 may be included: S401, group all differential interferograms according to the same gradient extraction direction, with each group corresponding to a gradient extraction direction.

[0068] Here, we can first obtain the spatial domain oriented gradient fields of all calculated differential interferograms in multiple directions. These oriented gradient fields are then divided according to their corresponding gradient extraction directions, and all gradient fields in the same direction are grouped into the same group to obtain multiple groups. For example, taking four preset gradient extraction directions as an example, the four groups correspond to the 0° direction, the 45° direction, the 90° direction, and the 135° direction, respectively.

[0069] S402, for each gradient extraction direction group, the gradient values ​​of all directional gradient fields in the group at the same pixel position are accumulated to obtain the gradient sum of each pixel position; the gradient sum of each pixel position is divided by the total number of directional gradient fields in the group to obtain the stacked average gradient value of each pixel position; the stacked average gradient values ​​of all pixel positions constitute the stacked average gradient field of the gradient extraction direction.

[0070] Specifically, when processing a grouping in a certain direction, the gradient fields of each direction within that group can be traversed. For each pixel coordinate, the values ​​of all gradient fields at that coordinate are summed to obtain the gradient summation of that pixel location across all interferograms. Repeating the above operation for all pixel locations yields a gradient summation matrix with the same size as the image.

[0071] Furthermore, after obtaining the gradient summation at each pixel location, the gradient summation can be used as the numerator, and the total number of directional gradient fields within the group can be used as the denominator for division. The result is the stacked average gradient value of that pixel.

[0072] Here, the stacked average gradient values ​​calculated for each pixel position are rearranged into a two-dimensional array according to the row and column coordinates of the pixels to obtain the stacked average gradient field for that gradient extraction direction. The stacked average gradient field can reflect the average intensity of the phase gradient in that direction over a long time series. Noise components are suppressed due to the cancellation of positive and negative signals, while the gradient signal at the deformation boundary is preserved.

[0073] Furthermore, after obtaining the stacked average gradient field of all gradient extraction directions, the phase gradient rate map of the area to be monitored can be determined based on these stacked average gradient fields. This phase gradient rate map can intuitively show the areas where the surface deformation rate changes significantly in space. The high value areas on it correspond to the abnormal deformation gradient locations such as landslide boundaries, fault traces, or the edges of subsidence basins, which can be used for early identification and survey of geological hazards.

[0074] Specifically, fusing stacked average gradient fields from multiple directions into a final phase gradient rate map may include the following steps (I) to (III): (I) Take the absolute value of the stacked average gradient field in each gradient extraction direction to obtain the absolute gradient field in each direction.

[0075] Here, in order to eliminate the influence of the positive or negative sign of the gradient direction on the fusion result, we can first perform a pixel-by-pixel absolute value operation on the stacked average gradient field of each direction, converting negative gradient values ​​into positive values, and obtain the absolute gradient field of each direction.

[0076] (II) The absolute gradient fields of all gradient extraction directions are added together at the pixel level to obtain the initial map of the fused phase gradient rate.

[0077] Furthermore, after completing the absolute value operation, the absolute value gradient fields in each direction can be numerically added at the same pixel position to achieve the fusion of gradient information in multiple directions, resulting in an initial image that integrates the deformation boundary intensity of each orientation.

[0078] (III) The initial phase gradient rate map is filtered to obtain the phase gradient rate map of the region to be monitored.

[0079] Specifically, since there may be residual isolated noise points in the initial phase gradient rate map, which appear as random spikes or salt-and-pepper outliers, the initial phase gradient rate map can be filtered in order to suppress noise interference and enhance the clarity of deformation boundaries. For example, median filtering can be used, which replaces the original value of a pixel with the median value of all pixels in the neighborhood of that pixel. This eliminates isolated noise while preserving edge clarity, and finally obtains a high-quality phase gradient rate map of the area to be monitored.

[0080] For example, in the PGR calculation process, taking the stacking and averaging of the four gradient extraction directions (i.e., 0°, 45°, 90°, and 135°) of N differential interferograms as an example, its expression in the phase gradient rate calculation can be expressed as: ; ; in, It is represented as the average gradient value obtained by stacking and averaging the spatial domain directional gradient fields along the gradient extraction direction d in N differential interferograms at pixel position x in the time dimension; It is represented as the final phase gradient rate value at pixel position x, which is the result of taking the absolute value of the stacked average gradient field in four directions and then calculating the arithmetic mean.

[0081] It is important to note that the order of stacking followed by taking the absolute value is crucial. During the stacking stage, since the noise is a zero-mean random variable with random positive and negative values, the positive and negative components partially cancel each other out during summation. Meanwhile, the gradient signal of the deformation boundary has a consistent direction (the sign of the deformation rate does not change with time), so the signal amplitude is preserved after stacking, and the noise is compressed. If the absolute value is taken before stacking, the noise will follow a semi-normal distribution (with a non-zero mean), and stacking cannot eliminate the bias, resulting in an overall increase in the PGR background level, severely weakening the discernibility of weak deformation boundaries.

[0082] The high-precision InSAR phase gradient rate calculation method, apparatus, and medium provided in this disclosure map the wrapped phase value into a complex exponential signal and extract the real and imaginary parts. Combined with a pre-constructed Riesz-Gauss transfer function, gradient extraction and high-frequency noise attenuation are completed simultaneously in the frequency domain. While suppressing interferogram noise and residual phase interference, the gradient information of the deformation edge is preserved. After time-dimensional stacking and averaging, the phase gradient rate information of the large gradient deformation region on the surface can be stably obtained, thereby improving the accuracy and robustness of phase gradient rate estimation.

[0083] To verify the performance of the proposed method in calculating the phase gradient rate, this disclosure obtained real data from a reservoir area and collected 30 Sentinel-1 down-orbit images covering the study area, spanning from January 2, 2024 to December 27, 2024. After registering all SAR images to the master image, a maximum time baseline of 36 days and a maximum spatial baseline of 200 meters were set, generating 74 interferometric pairs. (Refer to...) Figure 5 As shown, the reservoir area is located on the southeastern boundary of a plateau (i.e., Figure 5 (a) Located in a river basin, the area is frequently affected by geological disasters due to both tectonic activity and river erosion. The study area has drastic topography and dense vegetation cover, resulting in generally low SAR interferometric coherence, making it an ideal test area for verifying the noise resistance performance of gradient operators.

[0084] Furthermore, after obtaining differential interferograms from the registered SAR images through differential interferometry, the phase gradients of each differential interferogram were calculated using the central difference method, the Sobel operator, and the Riesz-Gauss operator, respectively. These were then stacked and averaged over time and fused in multiple directions to obtain the final normalized phase gradient rate result, i.e., the result obtained using the central difference method. Figure 5 (b) The result of the Sobel operator is Figure 5 (c) and the results of the Riesz-Gauss operator are Figure 5 (d). Here, from the perspective of the overall background noise level, the central difference method ( Figure 5 (b) The background PGR value is generally high, the yellow-green mottled texture is widely distributed in the figure, a large amount of high-value noise is scattered in the mountainous and valley areas, and the background uniformity is poor. Sobel operator ( Figure 5 (c) shows a significant improvement in background noise compared to the central difference method, with a more uniform green hue and a significant reduction in the number of high-value noise points. However, sporadic medium-intensity noise remnants can still be identified in some ridge areas. The Riesz-Gauss operator ( Figure 5 (d) has the cleanest and most uniform background, with a consistent low PGR value (dark green) across almost the entire area. Only in areas with significant deformation gradients, such as the sides of the valley and near the hydropower station, does it show a higher PGR value.

[0085] Here, in order to quantitatively evaluate the noise suppression performance of the three operators, we also... Figure 5 (b)- Figure 5Statistical analysis was performed on the normalized PGR results in (d). Background pixels were defined as those with PGR values ​​below the median plus the interquartile range. Contrast ratio was defined as the ratio of the mean of the top 5% highest PGR pixels to the background mean. Signal-to-noise ratio (SNR) was defined as the difference between the mean of the top 5% highest PGR pixels and the background mean divided by the background standard deviation, with local standard deviations calculated within a 5×5 sliding window. In terms of background noise levels, the mean background value of the central difference method was 0.16, and the standard deviation was 0.05; the mean background value of the Sobel operator was 0.13, and the standard deviation was 0.04; while the mean background value of the Riesz-Gauss operator was only 0.06, and the standard deviation was only 0.02. Its mean background value was only 37.5% of that of the central difference method, and its standard deviation was only 40% of that of the central difference method. In terms of signal discrimination capability, the central difference method has a contrast ratio of 2.07 and a signal-to-noise ratio (SNR) of 3.24; the Sobel operator has a contrast ratio of 2.22 and an SNR of 3.45; and the Riesz-Gauss operator achieves a contrast ratio of 3.43 and an SNR of 5.52, representing a 65.7% improvement in contrast and a 70.4% improvement in SNR compared to the central difference method. Regarding spatial uniformity, the local standard deviation of the 5×5 window for the central difference method is 0.05, for the Sobel operator it is 0.04, and for the Riesz-Gauss operator it is 0.02, a reduction of 60%. Thus, through the quantitative comparative analysis of background mean, background standard deviation, contrast ratio, SNR, and local standard deviation, it can be shown that the Riesz-Gauss operator produces a smoother and more uniform background, and its noise suppression performance is significantly better than that of the central difference method and the Sobel operator.

[0086] For example, refer to respectively Figure 6-8 The image shows a comparison of PGR calculation results for three typical landslide areas using three different gradient operators. In each image, the top left corner... Figure 6 (a) Figure 7 (a) Figure 8 (a) All images are optical remote sensing images of the study area, with white lines marking the boundaries of the landslides; upper right sub-image. Figure 6 (b) Figure 7 (b) Figure 8 (b) are all normalized PGR results calculated using the central difference method; lower left. Figure 6 (c) Figure 7 (c) Figure 8 (c) All are normalized PGR results calculated using the Sobel operator; bottom right Figure 6 (d) Figure 7 (d) Figure 8 (d) All are normalized PGR results calculated using the Riesz-Gauss operator proposed in this invention. The color scale ranges from 0 to 1, with green indicating a stable region with low PGR values ​​and red indicating anomaly regions with high PGR values.

[0087] Specifically, refer to Figure 6 As shown, Figure 6 (a) is an optical remote sensing image of the first typical landslide area, with white lines outlining the boundary of the landslide body. Figure 6 (b) is the normalized PGR result calculated using the central difference method. The background noise in the figure is significant, and the whole image shows a yellow-green noise pattern. The contrast between the landslide body and the surrounding stable area is low, and the deformation boundary features are not clear enough. Figure 6 (c) shows the normalized PGR result calculated using the Sobel operator. The background noise in the figure is improved compared to the central difference method, but there are still obvious noise textures in the stable region. The distinction between the high PGR signal and the background noise in the landslide area is limited. Figure 6 (d) shows the normalized PGR results calculated using the Riesz-Gauss operator. The background in the figure is the cleanest, the stable area shows uniform low green values, and the high PGR anomalies (red areas) within the landslide body are concentrated and have clear boundaries, which match well with the landslide boundaries marked in the optical image.

[0088] Similarly, refer to Figure 7 , Figure 7 (a) is an optical remote sensing image of the second typical landslide area, with the white lines marking the landslide boundary. Figure 7 (b) is the result calculated by the central difference method. The background noise is prominent, the landslide boundary is blurred, and a large amount of high-value noise interference makes it difficult to identify the deformation edge. Figure 7 (c) shows the results of the Sobel operator calculation. The background noise has been reduced, but the landslide boundary is still not sharp enough, and the gradient response broadens in some areas. Figure 7 (d) shows the results calculated using the Riesz-Gauss operator, with a uniformly low background, clear and continuous landslide boundaries, and a height consistent with the optical image boundaries. (Refer to...) Figure 8 , Figure 8 (a) is an optical remote sensing image of the third typical landslide area, with white lines indicating the landslide boundaries. Figure 8 (b) shows the result of the central difference method, with a high noise level, and the landslide body is almost submerged in the background noise. Figure 8 (c) shows the Sobel operator result, which has limited noise suppression and the landslide boundary is faintly visible but lacks contrast. Figure 8 (d) shows the results of the Riesz-Gauss operator. The background is clean, and the red high-value areas of the landslide body are sharp and have complete boundaries, which is significantly better than the previous two methods. The comparison of the three methods shows that the Riesz-Gauss operator, through the combination of the full-pass differential of the Riesz transform and the Gaussian low-pass filter, effectively suppresses the transmission of high-frequency noise to the PGR results, and significantly improves the signal-to-noise ratio and spatial identification of the deformation boundary.

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

[0090] Based on the same inventive concept, this disclosure also provides a high-precision InSAR phase gradient rate calculation device corresponding to the high-precision InSAR phase gradient rate calculation method. Since the principle of the device in this disclosure is similar to the high-precision InSAR phase gradient rate calculation 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.

[0091] Reference Figure 9 The diagram shown is a schematic of a high-precision InSAR phase gradient rate calculation device 900 provided in an embodiment of this disclosure. The device includes: Data acquisition module 901 is used to acquire SAR dataset and DEM data of the area to be monitored; and generate differential interferograms based on the SAR dataset and DEM data; The data mapping module 902 is used to map the entangled phase value corresponding to each pixel in the differential interferogram to a complex exponential signal for each differential interferogram, and to determine the real part data matrix and the imaginary part data matrix of the differential interferogram based on the real part data and imaginary part data of the entangled phase value corresponding to each pixel extracted from the complex exponential signal. The gradient field determination module 903 is used to determine the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions based on the real part data matrix and the imaginary part data matrix corresponding to each difference interferogram, as well as the pre-constructed Riesz-Gauss transfer function family; wherein, the pre-constructed Riesz-Gauss transfer function family includes Riesz-Gauss transfer functions corresponding to the multiple gradient extraction directions respectively; The data stacking module 904 is used to perform time-dimensional stacked averaging on the spatial domain directional gradient fields of all differential interferograms in the same gradient extraction direction, and to determine the phase gradient rate map of the region to be monitored based on the stacked averaging results of multiple gradient extraction directions.

[0092] In some possible embodiments, the data acquisition module 901 is specifically used for: An initial SAR image set is determined based on the SAR dataset; wherein, the initial SAR image set includes multiple SAR images of the area to be monitored; The SAR master image in the initial SAR image set is determined. Based on the SAR master image, other SAR images in the initial SAR image set except for the SAR master image are registered according to the DEM data. The target SAR image set is determined based on the registration results. The target SAR image set includes multiple registered target SAR images. A baseline set is determined based on a preset spatiotemporal baseline threshold and the target SAR image set; wherein, the baseline set includes multiple interferometric pairs, and each interferometric pair includes two target SAR images; Based on the DEM data, a differential interferogram is generated corresponding to each interference pair.

[0093] In some possible embodiments, the gradient field determination module 903 is specifically used for: For each difference interferogram, perform the following operations: A frequency domain transformation is performed on the real part data matrix of the differential interferogram to obtain the frequency domain real part data matrix corresponding to the differential interferogram; and a frequency domain transformation is performed on the imaginary part data matrix of the differential interferogram to obtain the frequency domain imaginary part data matrix corresponding to the differential interferogram. For each of the preset multiple gradient extraction directions, the real part data matrix and the imaginary part data matrix in the frequency domain corresponding to the differential interferogram are multiplied element-wise in the frequency domain with the Riesz-Gauss transfer function corresponding to the direction to obtain the real part gradient data matrix and the imaginary part gradient data matrix in the frequency domain in the direction. The real part gradient data matrix in the frequency domain is transformed in the spatial domain in each direction to obtain the real part gradient field in the spatial domain in that direction; the imaginary part gradient data matrix in the frequency domain is transformed in the spatial domain in each direction to obtain the imaginary part gradient field in the spatial domain in that direction. The spatial domain directional gradient field is synthesized based on the spatial domain real part gradient field and spatial domain imaginary part gradient field in each direction.

[0094] In some possible embodiments, the Riesz-Gauss transfer function is pre-constructed in the following manner: Construct a Riesz transform transfer function, wherein the amplitude-frequency response of the Riesz transform transfer function along any gradient extraction direction is constant, and use it to achieve isotropic gradient extraction; A Gaussian low-pass transfer function is constructed, wherein the amplitude-frequency response of the Gaussian low-pass transfer function decays exponentially with increasing frequency, in order to suppress high-frequency noise components. Multiplying the Riesz transform transfer function and the Gaussian low-pass transfer function in the frequency domain yields the basic Riesz-Gauss transfer function; For each preset gradient extraction direction, the basic Riesz-Gauss transfer function is projected along the gradient extraction direction to obtain the Riesz-Gauss transfer function corresponding to the gradient extraction direction. Based on the Riesz-Gauss transfer functions corresponding to all gradient extraction directions, the family of Riesz-Gauss transfer functions is constructed.

[0095] In some possible embodiments, the gradient field determination module 903 is specifically used for: For each gradient extraction direction, a pre-constructed Riesz-Gauss transfer function corresponding to the gradient extraction direction is obtained; wherein, the Riesz-Gauss transfer function is represented in the frequency domain as a numerical matrix with the same size as the real part data matrix in the frequency domain; Multiply the real part data matrix in the frequency domain with the elements at the same positions in the numerical matrix to obtain the real part gradient data matrix in the frequency domain; and multiply the imaginary part data matrix in the frequency domain with the elements at the same positions in the numerical matrix to obtain the imaginary part gradient data matrix in the frequency domain.

[0096] In some possible embodiments, the data stacking module 904 is specifically used for: All differential interferograms are grouped according to the same gradient extraction direction, with each group corresponding to a gradient extraction direction; For each gradient extraction direction group, the gradient values ​​of all directional gradient fields within the group at the same pixel position are accumulated to obtain the gradient sum of each pixel position; the gradient sum of each pixel position is divided by the total number of directional gradient fields within the group to obtain the stacked average gradient value of each pixel position; the stacked average gradient values ​​of all pixel positions constitute the stacked average gradient field of the gradient extraction direction.

[0097] In some possible embodiments, the data stacking module 904 is specifically used for: The absolute value of the stacked average gradient field in each gradient extraction direction is taken to obtain the absolute value gradient field in each direction. The absolute gradient fields of all gradient extraction directions are added at the pixel level to obtain the initial phase gradient rate map after fusion. The initial phase gradient rate map is filtered to obtain the phase gradient rate map of the region to be monitored.

[0098] Based on the same technical concept, this disclosure also provides a computer device. (See also...) Figure 10The diagram shown is a structural schematic of a computer device 1000 provided in an embodiment of this disclosure, including a processor 1001, a memory 1002, and a bus 1003. The memory 1002 stores execution instructions and includes a main memory 10021 and an external memory 10022. The main memory 10021, also called internal memory, is used to temporarily store computational data in the processor 1001 and data exchanged with external memory 10022 such as a hard disk. The processor 1001 exchanges data with the external memory 10022 through the main memory 10021.

[0099] In this embodiment, the memory 1002 is specifically used to store application code that executes the solution of this application, and its execution is controlled by the processor 1001. That is, when the computer device 1000 is running, the processor 1001 communicates with the memory 1002 through the bus 1003, so that the processor 1001 executes the application code stored in the memory 1002, and then executes the method described in any of the foregoing embodiments.

[0100] The memory 1002 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.

[0101] Processor 1001 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.

[0102] It is understood that the structures illustrated in the embodiments of this application do not constitute a specific limitation on the computer device 1000. In other embodiments of this application, the computer device 1000 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.

[0103] This disclosure also provides a computer-readable storage medium storing a computer program that, when executed by a processor, performs the steps of the high-precision InSAR phase gradient rate calculation method described in the above-described method embodiments. The storage medium can be either volatile or non-volatile computer-readable storage.

[0104] 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 high-precision InSAR phase gradient rate calculation method described in the above method embodiments. For details, please refer to the above method embodiments, which will not be repeated here.

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

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

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

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

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

[0110] 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 high-precision InSAR phase gradient rate calculation method, characterized in that, include: Acquire SAR datasets and DEM data for the area to be monitored; Differential interferograms are generated based on the SAR dataset and DEM data; For each differential interferogram, the winding phase value corresponding to each pixel in the differential interferogram is mapped to a complex exponential signal; and based on the real part data and imaginary part data of the winding phase value corresponding to each pixel extracted from the complex exponential signal, the real part data matrix and imaginary part data matrix of the differential interferogram are determined. Based on the real part data matrix and the imaginary part data matrix corresponding to each difference interferogram, and the pre-constructed family of Riesz-Gauss transfer functions, the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions is determined; wherein, the pre-constructed family of Riesz-Gauss transfer functions includes Riesz-Gauss transfer functions corresponding to the multiple gradient extraction directions respectively. The spatial domain directional gradient fields of all differential interferograms in the same gradient extraction direction are stacked and averaged in the time dimension, and the phase gradient rate map of the region to be monitored is determined based on the stacked and averaged results of multiple gradient extraction directions.

2. The method according to claim 1, characterized in that, The generation of differential interferograms based on the SAR dataset and DEM data includes: An initial SAR image set is determined based on the SAR dataset; wherein, the initial SAR image set includes multiple SAR images of the area to be monitored; The SAR master image in the initial SAR image set is determined. Based on the SAR master image, other SAR images in the initial SAR image set except for the SAR master image are registered according to the DEM data. The target SAR image set is determined based on the registration results. The target SAR image set includes multiple registered target SAR images. A baseline set is determined based on a preset spatiotemporal baseline threshold and the target SAR image set; wherein, the baseline set includes multiple interferometric pairs, and each interferometric pair includes two target SAR images; Based on the DEM data, a differential interferogram is generated corresponding to each interference pair.

3. The method according to claim 1, characterized in that, The determination of the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions, based on the real and imaginary data matrices corresponding to each difference interferogram and a pre-constructed family of Riesz-Gauss transfer functions, includes: For each difference interferogram, perform the following operations: A frequency domain transformation is performed on the real part data matrix of the differential interferogram to obtain the frequency domain real part data matrix corresponding to the differential interferogram; and a frequency domain transformation is performed on the imaginary part data matrix of the differential interferogram to obtain the frequency domain imaginary part data matrix corresponding to the differential interferogram. For each of the preset multiple gradient extraction directions, the real part data matrix and the imaginary part data matrix in the frequency domain corresponding to the differential interferogram are multiplied element-wise in the frequency domain with the Riesz-Gauss transfer function corresponding to the direction to obtain the real part gradient data matrix and the imaginary part gradient data matrix in the frequency domain in the direction. The real part gradient data matrix in the frequency domain is transformed in the spatial domain in each direction to obtain the real part gradient field in the spatial domain in that direction; the imaginary part gradient data matrix in the frequency domain is transformed in the spatial domain in each direction to obtain the imaginary part gradient field in the spatial domain in that direction. The spatial domain directional gradient field is synthesized based on the spatial domain real part gradient field and spatial domain imaginary part gradient field in each direction.

4. The method according to claim 3, characterized in that, The Riesz-Gauss transfer function family is pre-constructed in the following manner: Construct a Riesz transform transfer function, wherein the amplitude-frequency response of the Riesz transform transfer function along any gradient extraction direction is constant, and use it to achieve isotropic gradient extraction; A Gaussian low-pass transfer function is constructed, wherein the amplitude-frequency response of the Gaussian low-pass transfer function decays exponentially with increasing frequency, in order to suppress high-frequency noise components. Multiplying the Riesz transform transfer function and the Gaussian low-pass transfer function in the frequency domain yields the basic Riesz-Gauss transfer function; For each preset gradient extraction direction, the basic Riesz-Gauss transfer function is projected along the gradient extraction direction to obtain the Riesz-Gauss transfer function corresponding to the gradient extraction direction. Based on the Riesz-Gauss transfer functions corresponding to all gradient extraction directions, the family of Riesz-Gauss transfer functions is constructed.

5. The method according to claim 4, characterized in that, The step of performing frequency domain element-wise multiplication of the real and imaginary data matrices corresponding to the differential interferogram with the Riesz-Gauss transfer function corresponding to the direction includes: For each gradient extraction direction, a pre-constructed Riesz-Gauss transfer function corresponding to the gradient extraction direction is obtained; wherein, the Riesz-Gauss transfer function is represented in the frequency domain as a numerical matrix with the same size as the real part data matrix in the frequency domain; Multiply the real part data matrix in the frequency domain with the elements at the same positions in the numerical matrix to obtain the real part gradient data matrix in the frequency domain; and multiply the imaginary part data matrix in the frequency domain with the elements at the same positions in the numerical matrix to obtain the imaginary part gradient data matrix in the frequency domain.

6. The method according to claim 1, characterized in that, The step of performing time-dimension stacked averaging on the spatial domain directional gradient fields of all differential interferograms along the same gradient extraction direction includes: All differential interferograms are grouped according to the same gradient extraction direction, with each group corresponding to a gradient extraction direction; For each gradient extraction direction group, the gradient values ​​of all directional gradient fields within the group at the same pixel position are accumulated to obtain the gradient sum of each pixel position; the gradient sum of each pixel position is divided by the total number of directional gradient fields within the group to obtain the stacked average gradient value of each pixel position; the stacked average gradient values ​​of all pixel positions constitute the stacked average gradient field of the gradient extraction direction.

7. The method according to claim 6, characterized in that, The determination of the phase gradient rate map of the region to be monitored based on the stacked averaging results of multiple gradient extraction directions includes: The absolute value of the stacked average gradient field in each gradient extraction direction is taken to obtain the absolute value gradient field in each direction. The absolute gradient fields of all gradient extraction directions are added at the pixel level to obtain the initial phase gradient rate map after fusion. The initial phase gradient rate map is filtered to obtain the phase gradient rate map of the region to be monitored.

8. A high-precision InSAR phase gradient rate calculation device, characterized in that, include: The data acquisition module is used to acquire SAR datasets and DEM data of the area to be monitored. Differential interferograms are generated based on the SAR dataset and DEM data; The data mapping module is used to map the entangled phase value corresponding to each pixel in the differential interferogram to a complex exponential signal for each differential interferogram, and to determine the real part data matrix and the imaginary part data matrix of the differential interferogram based on the real part data and imaginary part data of the entangled phase value corresponding to each pixel extracted from the complex exponential signal. The gradient field determination module is used to determine the spatial domain directional gradient field of each difference interferogram in multiple preset gradient extraction directions based on the real part data matrix and the imaginary part data matrix corresponding to each difference interferogram, as well as the pre-constructed Riesz-Gauss transfer function family; wherein, the pre-constructed Riesz-Gauss transfer function family includes Riesz-Gauss transfer functions corresponding to the multiple gradient extraction directions respectively. The data stacking module is used to perform time-dimensional stacked averaging of the spatial domain directional gradient fields of all differential interferograms in the same gradient extraction direction, and to determine the phase gradient rate map of the region to be monitored based on the stacked averaging results of multiple gradient extraction directions.

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