Landslide deformation monitoring method based on SAR image offset tracking and multi-stage filtering

Abstract

The present application belongs to the field of geological disaster monitoring and early warning, and discloses a landslide deformation monitoring method based on SAR image offset tracking and multi-stage filtering. The registered SAR image data set is divided into image pairs, the original offset field of each image pair is calculated by using the block phase cross-correlation algorithm, and the original offset field is comprehensively processed by multi-stage filtering based on correlation threshold and signal-to-noise ratio threshold filtering, median filtering, statistical aggregation, amplitude threshold filtering, vector consistency filtering and density clustering analysis, so that noise is effectively removed, systematic deviation is corrected, and real and coherent landslide deformation signals are refined, thereby providing a more accurate and reliable landslide deformation monitoring solution. The present application focuses on image pixel feature tracking and post-processing, integrates various technologies to filter noise step by step, can realize the identification and long-term dynamic monitoring of landslides, improves the accuracy and reliability of landslide monitoring, and provides technical support for preventing and controlling landslide disasters.

Description

Technical Field

[0001] This invention belongs to the field of geological disaster monitoring and early warning, specifically involving a landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering. Background Technology

[0002] Landslides are a major concern among geological disasters due to their high hazard and destructive power. Some slowly changing landslides often exhibit progressively accelerating deformation before entering the instability stage, providing a critical window for disaster monitoring and early warning. Continuous and precise capture of this type of deformation is a core element of landslide disaster risk management.

[0003] Traditional contact-based point monitoring technologies, such as Global Navigation Satellite Systems (GNSS) or total stations, while providing high-precision data, are limited by monitoring range, cost, and accessibility in the field, making them unsuitable for meeting the needs of early identification and dynamic monitoring of landslide hazards over wide areas and long time periods. Space remote sensing technology, especially spaceborne synthetic aperture radar (SAR), has become an important technical means for monitoring surface deformation due to its all-weather, all-day imaging capabilities.

[0004] Currently, interferometric measurement technology based on SAR data (InSAR technology) can achieve millimeter-level detection of surface deformation. However, InSAR technology relies on the phase stability of radar signals, and its application in landslide monitoring has inherent technical bottlenecks: 1) In areas with dense vegetation cover or severe deformation gradients, the signal is prone to decoherence, leading to measurement failure; 2) Limited by side-view imaging geometry, it is not sensitive to motions that are nearly orthogonal to the satellite line of sight (such as north-south); 3) For high-rate or nonlinear deformations, phase unwrapping is prone to errors, making it impossible to obtain reliable deformation results.

[0005] To overcome the aforementioned limitations of InSAR technology, Pixel Offset Tracking (POT) technology was proposed. This technology directly obtains a two-dimensional large gradient displacement field in the azimuth and range directions by calculating the cross-correlation of intensity or phase information of consecutive temporal SAR images, and is insensitive to decoherence. However, the application of existing POT technology still faces serious challenges: its measurement accuracy is limited by image resolution, and the original displacement field generally contains significant noise and systematic biases introduced by speckle noise, orbital errors, and matching errors. These interference signals often obscure the true landslide deformation signal, leading to misjudgment of deformation patterns.

[0006] More importantly, current post-processing methods for POT (Protocol on Landslide) technology are still imperfect, lacking a systematic processing flow capable of filtering out multivariate composite noise, correcting systematic biases, and extracting spatially continuous deformation fields. Therefore, how to robustly and accurately extract the true spatiotemporal evolution characteristics of landslides from raw POT observations with low signal-to-noise ratios is a technical challenge that urgently needs to be solved in this field. Summary of the Invention

[0007] The purpose of this invention is to address the problem of extracting spatially continuous landslide deformation signals from noisy, high-dimensional satellite image time series. It provides a landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering. Through hierarchical filtering steps, the landslide deformation signal can be successfully separated and extracted from the original image data. By combining multiple constraints from kinematics, topography, and spatial topology, the false alarm rate is greatly reduced, ensuring the authenticity and reliability of the identified landslide areas.

[0008] The above-mentioned objective of the present invention is achieved through the following technical solution:

[0009] A landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering includes the following steps:

[0010] Step 1: Acquire SAR satellite remote sensing data within a specified time range in the monitoring area, and collect the corresponding high-precision orbit data and geographic elevation data. Preprocess the acquired SAR satellite remote sensing data and register it in the geographic coordinate system to obtain the registered SAR image dataset.

[0011] Step 2: According to the chronological order of the time series, the registered SAR images in the data set are centrally adjacent to form image pairs; based on the phase information of the registered SAR images in the image pairs, the pixel offset in the registered SAR images of the image pairs is calculated using the block phase cross-correlation algorithm, so as to obtain the original offset field of each image pair.

[0012] Step 3: Perform correlation threshold filtering, signal-to-noise ratio threshold filtering, and median filtering on each original offset field in sequence to obtain the offset field after the first-level filtering process;

[0013] Step 4: Statistically aggregate the offset fields after the first-level filtering based on the time series to generate the corresponding median offset fields. Correct each offset field after the first-level filtering based on the median offset fields to obtain the median-corrected offset fields. Accumulate the median-corrected offset fields and remove the corresponding noise points to obtain the statistically aggregated and corrected offset fields. Apply multi-level refining filtering to retain spatially coherent pixels and obtain a smooth offset field.

[0014] Step 5: Convert the azimuth and range offsets of each smoothed offset field to the geographic coordinate system and synthesize them. The synthesized smoothed offset fields are arranged in time series to form a landslide deformation time series dataset.

[0015] As described above, step 1 specifically includes the following steps:

[0016] Define the monitoring area and select the time range and track type;

[0017] Acquire all available SAR satellite remote sensing data and corresponding high-precision orbit data for the monitoring area, and obtain geographic elevation data based on the latitude and longitude of the monitoring area;

[0018] From the acquired SAR satellite remote sensing data, strips covering the monitoring area are selected, and SLC images are generated from the strips using high-precision orbit data. Multi-view processing is then performed on the SLC images to obtain multi-view images.

[0019] Geocoding of multiview images is performed using geoelevation data to obtain geocoded images and geolookup tables;

[0020] Within the selected time range, select one geocoded image as the master image; using the master image as the reference, register all geocoded images using the geolookup table of each geocoded image within the time range to obtain the registered SAR image dataset.

[0021] As described above, step 2 specifically includes the following steps:

[0022] Divide the registered SAR images in the airspace image pair into blocks;

[0023] The spatial domain variation is transformed into the frequency domain for offset tracking using the following formula, resulting in the correlation map between pairs of blocks:

[0024] ,

[0025] In the formula, The correlation maps are for paired blocks. This indicates the size of the block window, which is the size of the block. and These represent the horizontal and vertical pixel coordinates of the pixels corresponding to the blocks, respectively. The Fourier transform result of the reference image blocks is used as an example. The complex conjugate of the Fourier transform result of the same location block in the target image; Indicates modulo, For inverse Fourier transform,

[0026] The reference image is the registered SAR image located earlier in the time series of the image pair, and the target image is the registered SAR image located later in the time series of the image pair.

[0027] The peak value of the correlation spectrum is denoted as the correlation coefficient, and the corresponding pixel is the matching point. The block offset between paired blocks is determined based on the peak value of the correlation spectrum.

[0028] The block offsets corresponding to each pair of blocks are stitched together according to the position of the block relative to the registered SAR image to obtain the original offset field of the corresponding image pair.

[0029] As described above, step 3 specifically includes the following steps:

[0030] Set the correlation threshold and signal-to-noise ratio threshold;

[0031] Perform correlation threshold filtering, signal-to-noise ratio threshold filtering, and median filtering sequentially according to the following requirements:

[0032] Correlation threshold filtering removes matching points whose correlation coefficient is less than the correlation threshold:

[0033] ,

[0034] in, The correlation coefficient of the matching points in the original offset field; The correlation threshold, Represents matching points; This represents the set of matching points in the original offset field; This represents the set of matching points in the offset field after correlation threshold filtering;

[0035] Signal-to-noise ratio (SNR) threshold filtering removes matching points with a SNR lower than the SNR threshold:

[0036] ,

[0037] in, This represents the signal-to-noise ratio at the matching point. , The noise correlation coefficient is equal to the standard deviation of the correlation coefficient of the original migration field; The signal-to-noise ratio threshold; This represents the set of matching points in the offset field after signal-to-noise ratio threshold filtering;

[0038] Median filtering:

[0039] ,

[0040] The result of removing unreliable matching points from the original offset field after applying correlation and signal-to-noise ratio thresholds; This is the offset field after removing the median offset, i.e., the offset field after the first-stage filtering process. This indicates taking the median. and These represent the azmith direction (azimuth) and the range direction (distance), respectively.

[0041] As mentioned above, step 4 involves statistically aggregating and correcting the offset field, specifically including the following steps:

[0042] Step 4.1.1: Based on the following formula, statistically aggregate the real and imaginary parts of the block offset corresponding to each matching point in all the offset fields after the first-level filtering to generate the corresponding median offset field:

[0043] ,

[0044] In the formula, The median displacement vector of the median migration field represents the matching point. Overall movement trend over time and These represent the horizontal and vertical pixel coordinates of the matching point, respectively. For the first Image pairs at matching points The calculated displacement vector, For matching points The total number of valid observations across all image pairs;

[0045] Step 4.1.2: The median offset field corresponding to the real part and the median offset field corresponding to the imaginary part constitute the median offset field;

[0046] Step 4.1.3: Subtract the median offset field from each offset field after the first-stage filtering process to obtain the median-corrected offset field;

[0047] Step 4.1.4: Accumulate all median-corrected offset fields in the time series to obtain the accumulated offset field. Based on the accumulated offset field, remove noise pixels from each median-corrected offset field according to the following rules to obtain the corresponding statistically aggregated corrected offset field:

[0048] Delete pixels in the accumulated offset field whose real part value is equal to 0 and whose real part value sign is opposite to that of most pixels in the accumulated offset field; delete pixels in the accumulated offset field whose imaginary part value is equal to 0 and whose imaginary part value sign is opposite to that of most pixels in the accumulated offset field.

[0049] As described above, step 4 yields the smooth offset field, which specifically includes the following steps:

[0050] Step 4.2.1: Based on the real and imaginary parts of the block offset of each matching point in the offset field after statistical aggregation correction, extract the amplitude value of the corresponding displacement vector and the direction angle of the displacement vector;

[0051] Step 4.2.2: Based on a preset amplitude threshold, filter and remove matching points in the offset field after statistical aggregation correction whose amplitude values ​​are lower than the amplitude threshold;

[0052] Step 4.2.3: Calculate the minimum angular deviation between the displacement vector of each matching point and the overall motion direction of its local neighborhood based on the following formula:

[0053] ,

[0054] It is the direction angle of the displacement vector of a single matching point, denoted as the direction angle. ; It is the direction angle. The median or smoothed average of the displacement vector direction angles of all matching points in the local neighborhood of the corresponding matching point. The modulo operator, The minimum angular deviation between the displacement vector of the matching point and the overall motion direction of its local neighborhood is denoted as the minimum angular deviation. ;

[0055] If the minimum angle deviation If the set threshold is exceeded, the matching point corresponding to this displacement vector is removed;

[0056] Step 4.2.4: Treat the high-density area of ​​matching points as the landslide body, and remove the isolated low-density matching points.

[0057] As described above, step 5 specifically includes the following steps:

[0058] Interpolate the smoothed offset field;

[0059] The smoothed offset field is transformed to the geographic coordinate system, and the eastward displacement component and the northward velocity component are generated and synthesized to obtain the synthesized smoothed offset field.

[0060] The synthesized smooth offset fields within the selected time range are arranged in chronological order to form a landslide deformation time series dataset.

[0061] A computer device includes a memory and a processor, the memory storing a computer program, characterized in that the processor executes the computer program to implement steps 1 to 5 of the landslide deformation monitoring method based on SAR image offset tracking and multi-level filtering as described in any one of the claims.

[0062] A computer-readable storage medium having a computer program stored thereon, characterized in that, when the computer program is executed by a processor, it implements steps 1 to 5 of any of the landslide deformation monitoring methods based on SAR image offset tracking and multi-level filtering as described in the present invention.

[0063] A computer program product, comprising a computer program, characterized in that, when the computer program is executed by a processor, it implements steps 1 to 5 of any of the landslide deformation monitoring methods based on SAR image offset tracking and multi-level filtering described in the present invention.

[0064] Compared with the prior art, the present invention has the following advantages:

[0065] 1. A more accurate and reliable method for landslide deformation monitoring based on SAR imagery is provided: This method integrates a multi-level filtering comprehensive processing flow. By constructing a signal processing flow that removes specific types of interference step by step, it first removes large outliers and systematic biases, then removes isolated noise points, and then suppresses spatiotemporal random noise and retains continuous landslide signals. Therefore, it can effectively remove noise, correct systematic biases, and refine real and coherent landslide deformation signals, ensuring the accuracy and reliability of the results.

[0066] 2. Enhance the spatial continuity and physical realism of landslide deformation: Based on the characteristics of landslide deformation, vector consistency filtering and clustering methods can effectively eliminate abnormal isolated noise points, retain spatially coherent deformation patterns, and enhance the spatial continuity and physical realism of landslide deformation.

[0067] 3. Identification and monitoring of small and slow landslides: Small and slow landslides have small amplitude and slow speed, and many general tools cannot effectively capture their deformation. This method can retain continuous small deformation signals in time and space, and realize the identification and monitoring of small and slow landslides.

[0068] Compared with traditional early warning systems, this invention is optimized specifically for the characteristics of SAR pixel tracking, enabling more accurate assessment of its motion patterns. Furthermore, it improves the accuracy and reliability of monitoring through automated processing and multi-level filtering, achieving more comprehensive and efficient risk management and disaster reduction measures for landslide disasters. Attached Figure Description

[0069] Figure 1 This is a flowchart of the present invention.

[0070] Figure 2 This is a schematic diagram showing the results of Embodiment 2 of the present invention. Figure 2 (a) ~ Figure 2 (j) is a deformation diagram of the landslide area corresponding to the time period in 2024. Figure 2(k) is a schematic diagram of the landslide area boundary. Figure 2 (l) is a diagram illustrating the changes in shape and deformation. Detailed Implementation

[0071] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to embodiments. It should be understood that the embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0072] Example 1:

[0073] A landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering, such as Figure 1 As shown, it includes the following steps:

[0074] Step 1: Acquire SAR satellite remote sensing data within a specified time range in the monitoring area, and collect the corresponding high-precision orbital data and geographic elevation data. Preprocess the acquired SAR satellite remote sensing data and register it in the geographic coordinate system to obtain the registered SAR image dataset.

[0075] In this embodiment, the acquired SAR satellite remote sensing data is Sentinel-1 SAR satellite remote sensing data. High-precision orbit and geographic elevation data corresponding to the Sentinel-1 SAR satellite remote sensing data are also collected. Specifically, the following steps are included:

[0076] Step 1.1: Define the monitoring area and select the time range and orbit type (ascending or descending).

[0077] Step 1.2: Obtain all available SAR satellite remote sensing data for the monitoring area through an automated interface, and acquire the high-precision orbit data corresponding to these SAR satellite remote sensing data (usually released about 20 days after the satellite acquires the SAR satellite remote sensing data, with accuracy down to the centimeter level). Obtain the geographic elevation data (Digital Elevation Model, DEM) based on the latitude and longitude of the monitoring area. In this embodiment, all available SAR satellite remote sensing data for the monitoring area are retrieved and acquired through the European Space Agency's data automation interface. The open-source tool sentineleof is used to automatically acquire the high-precision orbit data corresponding to all SAR satellite remote sensing data along the path. Based on the latitude and longitude of the corners recorded in the SLC parameter file, the tool elevation is used to automatically download the geographic elevation data covering the monitoring area based on the latitude and longitude. The dem_import command is used to convert the format and import it into the GAMMA software.

[0078] Step 1.3: Select strips covering the monitoring area from the acquired SAR satellite remote sensing data, generate SLC images from the strips using high-precision orbit data, and then perform multi-view processing on the SLC images to remove some noise effects, thus obtaining multi-view images.

[0079] Step 1.4: Using the obtained geographic elevation data (DEM data), geocode the generated multi-view image to obtain a geocoded image and a geographic lookup table; the geographic lookup table is used to find the correspondence between the pixel in the image coordinate system and the coordinate position in the geographic coordinate system in the geocoded image, thereby reducing geometric distortion caused by changes in orbit, terrain and satellite attitude.

[0080] Step 1.5: Within the time range selected in Step 1.1, select a geocoded image as the master image; using the master image as the reference, use the geolookup table of each geocoded image within the time range to register all geocoded images, so that the same point (such as the same rock or the same building) between each geocoded image is aligned in the pixel grid, thereby obtaining the registered SAR image dataset.

[0081] This embodiment uses the compressed package of raw SAR satellite remote sensing data in GAMMA software, employing the SLC_mosaic_S1_TOPS command to select strips and directly extract SLC images and SLC parameter files using high-precision orbital data. Imported geo-elevation data is used to geocode the selected master image and other multi-view images within the time range, generating a geo-lookup table and refining it to map pixels in the image coordinate system to their geographic coordinates. The ScanSAR_coreg command is used for registration, initial offset estimation, and resampling to obtain a unified, registered SAR image dataset based on the master image. In this embodiment, the multi_look command in GAMMA software is used to process the SLC images to obtain multi-view images, and a parameter file is created to record command calls.

[0082] Step 2: Following the chronological order of the time series, adjacent registered SAR images in the registered SAR image dataset obtained in Step 1 are grouped into image pairs. The registered SAR image earlier in the time series within the image pair is the reference image, and the registered SAR image later in the time series is the target image. Based on the phase information of the registered SAR images in the image pair, the pixel offset in the registered SAR images of the image pair is calculated using the Blocked-based Phase CrossCorrelation (BPCC) algorithm, thus obtaining the original offset field of each image pair. The original offset field of each image pair includes the offset of each image in the azimuth and range directions. The Blocked-based Phase CrossCorrelation (BPCC) algorithm is based on the Fourier transform displacement theorem and specifically includes the following steps:

[0083] Step 2.1: Divide the registered SAR images in the spatial image pair into blocks. Specifically, this can be done by sliding the block window in each registered SAR image of the image pair according to a preset step size, thereby extracting the corresponding blocks one by one. The block window size is the block size. The blocks can be arranged closely or overlapped. Overlapping blocks can improve the smoothness and robustness of the results, but will increase the computational load. In addition, in order to reduce the spectral leakage problem caused by manually cutting the image blocks, a window function (Hanning window or Gaussian window) needs to be applied to each small block to make the pixel values ​​at the block edge smoothly transition to 0, effectively suppressing high-frequency noise after Fourier transform and improving matching accuracy.

[0084] Step 2.2: For each pair of blocks, the spatial domain variation is transformed into the frequency domain according to the following formula to perform offset tracking and obtain the correlation map between each pair of blocks:

[0085] ,

[0086] In the formula, The correlation maps are for paired blocks. This indicates the size of the block window, which is the size of the block. and These represent the horizontal and vertical pixel coordinates of the pixels corresponding to the blocks, respectively. The Fourier transform result of the reference image blocks is used as an example. The complex conjugate of the Fourier transform result of the same location block in the target image; Calculate the phase difference between the two; Indicates modulo, Used to normalize the phase difference results and eliminate the influence of amplitude information; The cross-power spectrum of the paired blocks was calculated; The cross-power spectrum is converted into a correlation spectrum by inverse Fourier transform.

[0087] Step 2.3: The values ​​in the correlation map range from 0 to 1, and are close to 0 in most positions, with a significant peak only at a certain pixel. The peak value of the correlation map is recorded as the correlation coefficient, and the corresponding pixel is the matching point. The offset between paired blocks is determined based on the peak value of the correlation map. , This represents the offset in the horizontal pixel coordinate direction. This represents the offset in the vertical pixel coordinate direction; this is the block offset between different pairs of blocks in the reference image and the target image.

[0088] Step 2.4: Set the block offsets corresponding to each pair of blocks. The original offset fields of the corresponding image pairs are obtained by stitching together the blocks according to their positions relative to the registered SAR images.

[0089] In this embodiment, the offset of each block in the original offset field is stored in complex form, where the real part of the block offset represents the range offset and the imaginary part represents the azimuth offset.

[0090] The block phase cross-correlation algorithm BPCC in this invention calculates the cross-power spectrum of paired blocks (i.e. By normalizing the spectrum, the image amplitude information is eliminated, which provides strong robustness to complex radiation differences (nonlinear changes in pixel brightness) commonly caused by imaging geometric changes in SAR images. At the same time, by using only phase information to perform inverse Fourier transform, a very sharp correlation peak is generated at the actual offset position between the two images, making the offset positioning extremely accurate, thereby improving the reliability and accuracy of monitoring.

[0091] Block-based phase cross-correlation exhibits strong robustness against speckle noise in SAR images. Operating in the frequency domain, it focuses on the overall structural information of the blocks (reflected in the phase at frequency), and normalization effectively suppresses this noise. It is also less sensitive to changes in backscattering characteristics. SAR image intensity is highly sensitive to the physical properties of ground features; changes in the surface (soil moisture, plant growth, or snowmelt, etc.) can significantly alter backscattering intensity, while phase-based calculations are almost unaffected. The block-based approach and Fast Fourier Transform (FFT) result in higher computational efficiency, and the phase matching at frequency makes the correlation peaks very prominent, leading to a high signal-to-noise ratio and thus high-precision migration measurement results.

[0092] Compared to the commonly used pixel offset tracking algorithm POT, the biggest advantage of this step is its robustness to radiation changes and the high accuracy and reliability brought by sharp peaks, ensuring the potential to detect minute displacements. However, the signal-to-noise ratio of the offset result is lower, so we perform subsequent processing to suppress noise and enhance the signal.

[0093] Step 3: Based on the original offset fields of each image pair obtained in Step 2, perform the first-level median filtering to remove systematic bias: sequentially perform correlation threshold filtering and signal-to-noise ratio threshold filtering to filter out unreliable large values, and apply median filtering to correct the systematic bias of the entire monitoring area, obtaining the offset field after the first-level filtering; specifically including the following steps:

[0094] The original offset fields of all image pairs are quality assessed and corrected. Predetermined correlation and signal-to-noise ratio (SNR) thresholds are set, and correlation threshold filtering and SNR threshold filtering are performed sequentially to remove unreliable matching points. Furthermore, to correct for systematic biases caused by inter-image pair translation, median filtering is applied to calculate and remove the median offset across the entire monitoring area. The correlation filtering rules are as follows:

[0095] Correlation threshold filtering removes matching points whose correlation coefficient is less than the correlation threshold:

[0096] ,

[0097] in, The correlation coefficient of the matching points in the original offset field; The correlation threshold, Represents matching points; This represents the set of matching points in the original offset field; This represents the set of matching points in the offset field after correlation threshold filtering;

[0098] Signal-to-noise ratio (SNR) threshold filtering removes matching points with a SNR lower than the SNR threshold:

[0099] ,

[0100] in, This represents the signal-to-noise ratio at the matching point. , The noise correlation coefficient (equal to the standard deviation of the correlation coefficient of the original offset field); The signal-to-noise ratio threshold; This represents the set of matching points in the offset field after signal-to-noise ratio threshold filtering;

[0101] Median filtering:

[0102] ,

[0103] The result is the original offset field after removing unreliable matching points by applying correlation and signal-to-noise ratio thresholds. This is the offset field after removing the median offset, i.e., the offset field after the first-stage filtering process. This indicates taking the median. , These represent the azmith direction (azimuth) and the range direction (distance), respectively.

[0104] In this embodiment, the original offset fields obtained in step 2 are stored in complex form in GAMMA software. Each complex form of the original offset field is separated and converted into two real form files. The median_filter command is used to sequentially perform correlation threshold filtering and signal-to-noise ratio threshold filtering on the original offset fields to filter out unreliable large values. Then, median filtering is used to remove isolated outliers from the original offset fields, resulting in the offset fields after the first-level filtering for each image pair.

[0105] This step removes systematic errors such as those related to the track and atmosphere that may be greater than the landslide signal by setting correlation thresholds, signal-to-noise ratio thresholds, and median filtering, thereby highlighting the local relative motion of the landslide.

[0106] Step 4: Perform a second-level spatiotemporal refinement filter on the offset field after the first-level filtering process in Step 3: Statistically aggregate the offset field after the first-level filtering process to obtain the median offset field. Based on the median offset field, filter is performed from the time dimension to eliminate temporal random noise in the offset field after the first-level filtering process and enhance the continuous deformation signal. Apply multi-level refinement filtering to retain spatially coherent matching points, thereby preserving spatially coherent deformation patterns and removing isolated outliers to obtain a smooth offset field. Specifically, this includes the following steps:

[0107] Step 4.1: Statistically aggregate the offset fields after the first-level filtering based on the time series to generate the corresponding median offset fields. Correct each offset field after the first-level filtering based on the median offset fields to obtain the median-corrected offset fields. Accumulate the median-corrected offset fields and remove the corresponding noise points to obtain the statistically aggregated and corrected offset fields.

[0108] Step 4.1.1: Based on the following formula, statistically aggregate the real and imaginary parts of the block offset corresponding to each matching point in all the offset fields after the first-level filtering to generate the corresponding median offset field:

[0109] ,

[0110] In the formula, The median displacement vector of the median migration field represents the matching point. Overall movement trend over time and These represent the horizontal and vertical pixel coordinates of the matching point, respectively. For the first Image pairs at matching points The calculated displacement vector, For matching points The total number of valid observations across all image pairs;

[0111] Step 4.1.2: The median offset field corresponding to the real part and the median offset field corresponding to the imaginary part constitute the median offset field;

[0112] Step 4.1.3: Subtract the median offset field from each offset field after the first-stage filtering process to obtain the median-corrected offset field;

[0113] Step 4.1.4: Accumulate all median-corrected offset fields in the time series to obtain the accumulated offset field. Based on the accumulated offset field, remove noise pixels from each median-corrected offset field according to the following rules to obtain the corresponding statistically aggregated corrected offset field:

[0114] Delete pixels in the accumulated offset field whose real part value is equal to 0 and whose real part value sign is opposite to that of most pixels in the accumulated offset field; delete pixels in the accumulated offset field whose imaginary part value is equal to 0 and whose imaginary part value sign is opposite to that of most pixels in the accumulated offset field.

[0115] The above operations preserve long-term deformation signals while removing noise. Statistically enhanced signals are considered to indicate the presence of deformation with the same trend; otherwise, they are considered noise. Landslide movements are generally long-term movements, meaning that as long as they are not considered noise, landslide movements can be considered to exist.

[0116] Step 4.2: Apply the following multi-stage refined filtering to preserve spatially coherent deformation patterns and eliminate isolated outliers. The multi-stage refined filtering process includes, but is not limited to, amplitude threshold filtering, vector consistency filtering, and density-based clustering analysis for denoising.

[0117] Step 4.2.1: Based on the real and imaginary parts of the block offset of each matching point in the offset field after statistical aggregation correction, extract the amplitude value of the corresponding displacement vector and the direction angle of the displacement vector (i.e., the amplitude value and direction angle of the resultant vector of the real and imaginary parts).

[0118] Step 4.2.2: Based on the preset amplitude threshold, filter and remove matching points in the offset field after statistical aggregation correction whose amplitude values ​​are lower than the amplitude threshold, and exclude minor noise and inactive areas;

[0119] Step 4.2.3: Based on angular consistency, remove matching points corresponding to displacement vectors whose orientation angles are inconsistent with those of their local neighbors, thereby enhancing the coherence of motion patterns. The specific process is as follows:

[0120] The minimum angular deviation between the displacement vector of each matching point and the overall motion direction of its local neighborhood is calculated based on the following formula:

[0121] ,

[0122] It is the direction angle of the displacement vector of a single matching point, denoted as the direction angle. ; It is the direction angle. The median or smoothed average of the displacement vector direction angles of all matching points in the local neighborhood of the corresponding matching point. The local neighborhood refers to a set of all matching points within a neighborhood of size m×m, centered on the matching point of the displacement vector to be evaluated in the displacement vector field of the entire monitoring area. m is the side length of the local neighborhood. The spatial range covered by this set is the local neighborhood of the vector. The local neighborhood is set here to detect whether the motion patterns of the main deformation area (landslide) and its surroundings are consistent. Therefore, the size of the local neighborhood is set according to the specific landslide area. The modulo operator, That is Divide by The remainder is used to normalize the direction angle to In the interval, This is the minimum angular deviation between the displacement vector of this matching point and the overall motion direction of its local neighborhood, denoted as the minimum angular deviation. The range of values ​​is The plus or minus sign represents clockwise or counterclockwise. This indicates that the direction of the displacement vector at the matching point is consistent with the overall motion direction of its local neighborhood. This indicates that the direction of the displacement vector at the matching point is completely opposite to the overall motion direction of its local neighborhood.

[0123] If the minimum angle deviation If the set threshold is exceeded, the matching point corresponding to this displacement vector is judged as an abnormal matching point and removed.

[0124] Step 4.2.4: Density-Based Spatial Clustering of Application with Noise (DBSCAN) to remove noisy pixels: By identifying the spatial distribution density of matching points corresponding to displacement vectors, high-density areas of matching point distribution are identified as landslide bodies, while isolated low-density matching points are identified as noise and removed.

[0125] This embodiment statistically aggregates all the offset fields after the first-level filtering to obtain the median offset field, and then obtains the corresponding statistically aggregated corrected offset field based on the median offset field. Subsequently, multi-level refined filtering is applied, setting an amplitude threshold to remove matching points corresponding to abnormal displacement vectors in the statistically aggregated corrected offset field, and utilizing the vector consistency formula... By filtering the matching points corresponding to spatially discontinuous displacement vectors, a coherent deformation pattern is preserved. Density-based clustering analysis is applied to identify the displacement vector density, determine the landslide body, and exclude low-density matching points, resulting in a multi-level filtered and refined displacement field.

[0126] This step utilizes the essential difference between landslide signals and noise. First, statistical aggregation superimposes multiple measurement results to reduce random noise and capture persistent landslide signals. Then, through spatial filtering and cluster analysis, it further eliminates residual noise by utilizing the prior knowledge that landslide deformation must be spatially concentrated.

[0127] This invention obtains the offset results between consecutive time points through step 2 and obtains a series of independent displacement observations within discrete time intervals after pairwise filtering in step 4. In step 4, "statistical aggregation" accumulates the discrete displacements after quality control, and the continuous signal will continuously increase to form a clear trend, while random noise will be weakened, thereby reconstructing the cumulative deformation of each time point in the time series.

[0128] Step 5: Convert the smoothed offset field (including the offsets in the azimuth and range directions) obtained in Step 4 to the geographic coordinate system and synthesize it (in Step 1.4, the registration is to convert the pixels in the image coordinate system to the corresponding geographic coordinate system for grid registration, while the subsequent calculations are based on the image coordinate system, not the geographic coordinate system, so the coordinate system needs to be converted here). The synthesized smoothed offset field is arranged in chronological order to form a landslide deformation time series dataset.

[0129] This embodiment interpolates the smoothed offset field before converting it to the geographic coordinate system to maintain its spatial integrity.

[0130] Then, the interpolated smoothed offset field is converted to the geographic coordinate system to generate eastward and northward displacement components, which are then synthesized to obtain the corresponding synthesized smoothed offset field.

[0131] The synthesized smooth offset fields within the selected time range are arranged in chronological order to form a landslide deformation time series dataset, thus completing the long-term dynamic landslide deformation monitoring results. The remaining results are then mapped and output.

[0132] This invention provides an innovative processing flow that can robustly and accurately extract real deformation signals from noisy raw data, resulting in a time series dataset that can be used for in-depth analysis of various subsequent models.

[0133] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above methods.

[0134] Example 2:

[0135] Taking a landslide in a selected area from January 11 to December 11, 2024 as an example, its steep slope is adjacent to a highway. Long-term rainfall and freeze-thaw cycles have weakened the slope stability, threatening not only road safety but also severely impacting rescue and subsequent maintenance work. To prevent similar disasters and ensure subsequent repair work, the landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering, as described in this invention, was implemented for landslide deformation monitoring. Sentinel-1 SAR satellite remote sensing data covering the area in 2024 was selected and paired monthly. Corresponding high-precision orbital data was used to eliminate orbital influences, and geographic elevation data covering the monitoring area was selected for unified registration and deskewing of the SAR images. The block phase cross-correlation algorithm (BPCC) was used to calculate the migration of the paired image pairs monthly. For each image pair, pairwise filtering was used to set correlation coefficients and signal-to-noise ratio thresholds, and reliable migration values ​​were selected using the formula... Subtract the system mean deviation. Statistically aggregate the offsets of all image pairs to obtain the final offset field after the first-stage filtering process, and apply vector consistency. To enhance motion pattern coherence, unreliable displacement vector matching points are removed. Cluster analysis is used to identify the spatial distribution density of displacement vectors to eliminate noise, generating spatially continuous and high-precision landslide deformation results. The Block Phase Cross-Correlation (BPCC) algorithm uses a 32×32 block window, a 7×7 window for median filtering, and correlation and signal-to-noise ratio thresholds of 0.1 and 3, respectively. The vector consistency threshold is set to 20°. Deformation is calculated based on monthly pairings from 2024, as shown below. Figure 2 The landslide area shown corresponds to the time series variation diagram of landslide deformation, providing guidance for landslide identification, monitoring, and early warning. Figure 2 (a) The corresponding time period is from January 11, 2024 to February 15, 2024. Figure 2 (b) The corresponding time period is from February 15, 2024 to March 10, 2024. Figure 2(c) The corresponding time period is from March 10, 2024 to April 15, 2024. Figure 2 (d) corresponds to the period from April 15, 2024 to May 9, 2024. Figure 2 (e) The corresponding time period is from May 9, 2024 to June 14, 2024. Figure 2 (f) corresponds to the period from June 14, 2024 to August 13, 2024. Figure 2 (g) The corresponding time period is from August 13, 2024 to September 18, 2024. Figure 2 (h) corresponds to the period from September 18, 2024 to October 12, 2024. Figure 2 (i) The corresponding time period is from October 12, 2024 to November 17, 2024. ​ (j) corresponds to the period from November 17, 2024 to December 11, 2024.

[0136] This invention achieves more comprehensive landslide deformation monitoring by combining pixel feature tracking with various data processing techniques such as system bias correction, multi-level filtering, vector consistency, and cluster analysis. The landslide monitoring process built upon this approach can adjust the identifiable degree of landslide deformation according to the set landslide risk, effectively improving the accuracy of landslide displacement calculation and the consistency of spatial patterns.

[0137] Example 3:

[0138] A landslide deformation monitoring device based on SAR image migration tracking and multi-level filtering, used to implement the landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering described in Example 1, includes:

[0139] Data acquisition and preprocessing module: used to implement step 1 of embodiment 1;

[0140] Original offset field acquisition module: used to implement step 2 of embodiment 1;

[0141] First-level filtering module: used to implement step 3 of embodiment 1;

[0142] Secondary filtering module: used to implement step 4 of embodiment 1;

[0143] Landslide deformation time series dataset acquisition module: used to implement step 5 of embodiment 1.

[0144] Example 4:

[0145] This embodiment provides a computer device, including a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements steps 1 to 5 in the above embodiment 1.

[0146] Example 5:

[0147] This embodiment provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements steps 1 to 5 in Embodiment 1 above.

[0148] Example 6:

[0149] This embodiment provides a computer program product, including a computer program that, when executed by a processor, implements steps 1 to 5 in embodiment 1 above.

[0150] The above are merely preferred embodiments of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering, characterized in that, Includes the following steps: Step 1: Acquire SAR satellite remote sensing data within a specified time range in the monitoring area, and collect the corresponding high-precision orbit data and geographic elevation data. Preprocess the acquired SAR satellite remote sensing data and register it in the geographic coordinate system to obtain the registered SAR image dataset. Step 2: According to the chronological order of the time series, the registered SAR images in the data set are centrally adjacent to form image pairs; based on the phase information of the registered SAR images in the image pairs, the pixel offset in the registered SAR images of the image pairs is calculated using the block phase cross-correlation algorithm, so as to obtain the original offset field of each image pair. Step 3: Perform correlation threshold filtering, signal-to-noise ratio threshold filtering, and median filtering on each original offset field in sequence to obtain the offset field after the first-level filtering process; Step 4: Statistically aggregate the offset fields after the first-level filtering based on the time series to generate the corresponding median offset fields. Correct each offset field after the first-level filtering based on the median offset fields to obtain the median-corrected offset fields. Accumulate the median-corrected offset fields and remove the corresponding noise points to obtain the statistically aggregated and corrected offset fields. Multi-level refining filtering is applied to preserve spatially coherent pixels and obtain a smooth offset field; Step 5: Convert the azimuth and range offsets of each smoothed migration field to the geographic coordinate system and synthesize them. The synthesized smoothed migration fields are arranged in time series to form a landslide deformation time series dataset. Step 2 specifically includes the following steps: Divide the registered SAR images in the airspace image pair into blocks; The spatial domain variation is transformed into the frequency domain for offset tracking using the following formula, resulting in the correlation map between pairs of blocks: ， In the formula, The correlation maps are for paired blocks. This indicates the size of the block window, which is the size of the block. and These represent the horizontal and vertical pixel coordinates of the pixels corresponding to the blocks, respectively. The Fourier transform result of the reference image blocks is used as an example. The complex conjugate of the Fourier transform result of the same location block in the target image; Indicates modulo, For inverse Fourier transform, The reference image is the registered SAR image located earlier in the time series of the image pair, and the target image is the registered SAR image located later in the time series of the image pair. The peak value of the correlation spectrum is denoted as the correlation coefficient, and the corresponding pixel is the matching point. The block offset between paired blocks is determined based on the peak value of the correlation spectrum. The block offsets corresponding to each pair of blocks are stitched together according to the position of the block relative to the registered SAR image to obtain the original offset field of the corresponding image pair. Step 4, which involves obtaining the smooth offset field, specifically includes the following steps: Step 4.2.1: Based on the real and imaginary parts of the block offset of each matching point in the offset field after statistical aggregation correction, extract the amplitude value of the corresponding displacement vector and the direction angle of the displacement vector; Step 4.2.2: Based on a preset amplitude threshold, filter and remove matching points in the offset field after statistical aggregation correction whose amplitude values ​​are lower than the amplitude threshold; Step 4.2.3: Calculate the minimum angular deviation between the displacement vector of each matching point and the overall motion direction of its local neighborhood based on the following formula: ， It is the direction angle of the displacement vector of a single matching point, denoted as the direction angle. ; It is the direction angle. The median or smoothed average of the displacement vector direction angles of all matching points in the local neighborhood of the corresponding matching point. The modulo operator, The minimum angular deviation between the displacement vector of the matching point and the overall motion direction of its local neighborhood is denoted as the minimum angular deviation. ; If the minimum angle deviation If the set threshold is exceeded, the matching point corresponding to this displacement vector is removed; Step 4.2.4: Treat the high-density area of ​​matching points as the landslide body, and remove the isolated low-density matching points.

2. The landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering according to claim 1, characterized in that, Step 1 specifically includes the following steps: Define the monitoring area and select the time range and track type; Acquire all available SAR satellite remote sensing data and corresponding high-precision orbit data for the monitoring area, and obtain geographic elevation data based on the latitude and longitude of the monitoring area; From the acquired SAR satellite remote sensing data, strips covering the monitoring area are selected, and SLC images are generated from the strips using high-precision orbit data. Multi-view processing is then performed on the SLC images to obtain multi-view images. Geocoding of multiview images is performed using geoelevation data to obtain geocoded images and geolookup tables; Within the selected time range, select one geocoded image as the master image; using the master image as the reference, register all geocoded images using the geolookup table of each geocoded image within the time range to obtain the registered SAR image dataset.

3. The landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering according to claim 1, characterized in that, Step 3 specifically includes the following steps: Set the correlation threshold and signal-to-noise ratio threshold; Perform correlation threshold filtering, signal-to-noise ratio threshold filtering, and median filtering sequentially according to the following requirements: Correlation threshold filtering removes matching points whose correlation coefficient is less than the correlation threshold: ， in, The correlation coefficient of the matching points in the original offset field; This is the correlation threshold. Represents matching points; This represents the set of matching points in the original offset field; This represents the set of matching points in the offset field after correlation threshold filtering; Signal-to-noise ratio (SNR) threshold filtering removes matching points with a SNR lower than the SNR threshold: ， in, This represents the signal-to-noise ratio at the matching point. , The noise correlation coefficient is equal to the standard deviation of the correlation coefficient of the original migration field. The signal-to-noise ratio threshold; This represents the set of matching points in the offset field after signal-to-noise ratio threshold filtering; Median filtering: ， The result of removing unreliable matching points from the original offset field after applying correlation and signal-to-noise ratio thresholds; This is the offset field after removing the median offset, i.e., the offset field after the first-stage filtering process. This indicates taking the median. and These represent the azmith direction (azimuth) and the range direction (distance), respectively.

4. The landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering according to claim 1, characterized in that, Step 4, which obtains the statistically aggregated and corrected offset field, specifically includes the following steps: Step 4.1.1: Based on the following formula, statistically aggregate the real and imaginary parts of the block offset corresponding to each matching point in all the offset fields after the first-level filtering to generate the corresponding median offset field: ， In the formula, The median displacement vector of the median migration field represents the matching point. Overall movement trend over time and These represent the horizontal and vertical pixel coordinates of the matching point, respectively. For the first Image pairs at matching points The calculated displacement vector, For matching points The total number of valid observations across all image pairs; Step 4.1.2: The median offset field corresponding to the real part and the median offset field corresponding to the imaginary part constitute the median offset field; Step 4.1.3: Subtract the median offset field from each offset field after the first-stage filtering process to obtain the median-corrected offset field; Step 4.1.4: Accumulate all median-corrected offset fields in the time series to obtain the accumulated offset field. Based on the accumulated offset field, remove noise pixels from each median-corrected offset field according to the following rules to obtain the corresponding statistically aggregated corrected offset field: Delete pixels in the accumulated offset field whose real part value is equal to 0 and whose real part value sign is opposite to that of most pixels in the accumulated offset field. Delete pixels whose imaginary part value is equal to 0 in the accumulated offset field, and pixels whose imaginary part value has the opposite sign to that of most pixels in the accumulated offset field.

5. The landslide deformation monitoring method based on SAR image migration tracking and multi-level filtering according to claim 1, characterized in that, Step 5 specifically includes the following steps: Interpolate the smoothed offset field; The smoothed offset field is transformed to the geographic coordinate system, and the eastward displacement component and the northward velocity component are generated and synthesized to obtain the synthesized smoothed offset field. The synthesized smooth offset fields within the selected time range are arranged in chronological order to form a landslide deformation time series dataset.

6. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements steps 1 to 5 of the landslide deformation monitoring method based on SAR image offset tracking and multi-level filtering as described in any one of claims 1 to 5.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements steps 1 to 5 of the landslide deformation monitoring method based on SAR image offset tracking and multi-level filtering as described in any one of claims 1 to 5.

8. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by the processor, it implements steps 1 to 5 of the landslide deformation monitoring method based on SAR image offset tracking and multi-level filtering as described in any one of claims 1 to 5.

Images