Subgrade settlement and landslide hidden danger identification method and system along traffic corridor

Abstract

The invention discloses a roadbed settlement and landslide hidden danger identification method and system along a traffic corridor, and the method comprises the steps: obtaining a time sequence satellite-borne data set of a target region, and carrying out the format conversion of the data format of the time sequence satellite-borne data set; generating an interference pair network based on the time baseline threshold; generating a filtered differential interference image set; performing differential interference phase unwrapping processing on the filtered differential interference image set to obtain an unwrapped differential interference image set; performing atmospheric delay phase correction on the unwrapped differential interference image set to obtain an unwrapped interference image set after atmospheric correction; calculating to obtain the annual average settlement rate of the linear traffic corridor area and a time sequence accumulated deformation result; based on the annual average settlement rate of the linear traffic corridor area, the spatial settlement gradient field information of the linear traffic corridor area is extracted, the settlement rate and the spatial settlement gradient field information meeting the hidden danger judgment standard are output, the roadbed settlement and landslide hidden danger extraction and recognition result of the target area is obtained, and the road area detection efficiency is improved.

Description

Technical Field

[0001] The present invention relates to the field of image processing, and in particular to a method and system for identifying roadbed settlement and landslide hazards along a traffic corridor. Background Art

[0002] With the rapid development of society, the demand for infrastructure construction in remote areas is growing, especially in plateau regions. The construction of permafrost roadbeds is crucial for ensuring smooth transportation networks. The stability of existing permafrost roadbeds is severely impacted by seasonal freeze-thaw deformation, which not only threatens safe road operations but also has adverse impacts on the ecological environment. Traditional monitoring methods, such as leveling and GNSS, suffer from high costs, low efficiency, and limited spatial coverage, making them inadequate for the full-cycle monitoring of thousands of kilometers of transportation corridors.

[0003] Although existing InSAR technology can achieve large-scale monitoring, it has the following limitations for long linear projects: vegetation cover, summer precipitation, and glacier melting along the line cause radar signal decoherence; long-wavelength atmospheric delay phase masks local millimeter-level deformation signals; and there is a lack of deformation risk anomaly identification and analysis methods based on the spatial characteristics of linear projects. Summary of the Invention

[0004] To achieve the above-mentioned and other related purposes, the present invention discloses a method for identifying roadbed settlement and landslide hazards along a traffic corridor, comprising: Obtain the time-series spaceborne SAR SLC dataset of the target area and convert the format of the time-series spaceborne SAR SLC dataset to obtain the time-series spaceborne SAR SLC image set; According to the time series spaceborne SAR SLC image set, based on the time baseline threshold, an interferometric pair network is generated, which contains multiple groups of interferometric pairs. Perform registration and differential interferometry processing on the interference pairs in the interference pair network and generate filtered differential interferometry atlas; The minimum cost flow method based on irregular grid is used to perform differential interferometry phase unwrapping on the filtered differential interferometry atlas to obtain the unwrapped differential interferometry atlas. The atmospheric delay phase correction method in spatial domain is used to correct the atmospheric delay phase of the unwrapped differential interferogram to obtain the atmospheric-corrected unwrapped interferogram. The atmospherically corrected unwrapped interferograms were used to calculate the time series deformation, and the annual average settlement rate and the time series cumulative deformation results of the linear transportation corridor area were obtained. Based on the average annual settlement rate of the linear traffic corridor area, the Sobel operator is used to extract the spatial settlement gradient field information of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output, and the roadbed settlement and landslide hidden danger extraction and identification results of the target area are obtained.

[0005] Preferably, generating an interference pair network based on a time baseline threshold comprises: The SBAS network method was used to set the time baseline range and spatial baseline threshold to conduct the initial screening of the data set; Based on the principle of the shortest time baseline, the data after the initial screening were screened twice to obtain the interference pairs that met the conditions; According to the principle of minimizing the weighted sum of spatial vertical baseline, time baseline and Doppler centroid frequency difference, a common main image is selected and an interference pair network is generated.

[0006] Preferably, the selecting of the public main image includes: Constructing correlation coefficients , is the correlation coefficient when the mth image is the public main image, ; Among them, Bc, Tc, and fc are the critical values ​​of the spatial vertical baseline, time baseline, and Doppler centroid frequency difference respectively, and 、 、 are the absolute values ​​of the spatial vertical baseline, time baseline, and Doppler centroid frequency difference between the k-th image and the m-th image, respectively. L is the total number of images after screening. Function c is the relationship function, and α, β, and θ are the exponential factors of the three factors. The expression of function c is as follows: ; Where x is the independent variable and a is the parameter; The image corresponding to the minimum correlation coefficient value is selected as the public main image.

[0007] Preferably, the performing of registration and differential interferometry processing and generating a filtered differential interferometry atlas comprises: The interferometer pair is registered using orbit baseline assisted registration, cross-correlation registration and enhanced spectral diversity methods; The phase difference between the registered interference pairs is calculated using the complex conjugate multiplication method, and the flat ground phase and terrain phase signals are removed from the phase difference signal to obtain the original differential interferogram set; Calculate the interference pair correlation coefficient atlas; According to the interference pair correlation coefficient atlas, a filter is used to perform differential interferogram adaptive filtering on the original differential interferogram atlas to generate a filtered differential interferogram atlas.

[0008] Preferably, the calculation of the interference pair correlation coefficient atlas includes: Calculate the coherence of each interferometer pair ; , in, and are the complex values ​​of the common primary image and its secondary image in each interferometer pair, is the number of pixels in the sliding window, Represents the conjugate operation.

[0009] Preferably, the adaptive filtering of the differential interferograms on the original differential interferogram set using a filter includes a circular periodic mean Boxcar filtering method or a spectrum feature-based adaptive filtering Goldstein method.

[0010] Preferably, the step of performing differential interference phase unwrapping on the filtered differential interference atlas using a minimum cost flow method based on an irregular grid, and obtaining the unwrapped differential interference atlas comprises: A coherence threshold is preset, and high-coherence pixels with a coherence greater than the coherence threshold are selected from the filtered differential interference pattern set to form a scattered point combination; Taking the high coherence pixels as nodes, the high coherence pixels are connected according to preset connection conditions to form a triangular network diagram; The preset connection conditions are: , is the distance between any two highly coherent pixels, is the coherence between any two highly coherent pixels, is the maximum distance between all highly coherent pixels, is the minimum value of coherence among all high coherence pixels; For each edge in the triangular network, calculate its winding phase difference: , in and are the phases of the two highly coherent pixels on this edge; Finding the minimum integer number of cycles , so that the unwrapping phase difference is minimized, including: ; in, is the unwrapping phase difference; The untangling order is set based on the weight of each edge in the preset triangular network, and the path graph is constructed using the minimum production tree; The point with the most stable imaging signal in the filtered differential interferogram is selected as the reference point, and the phase value is propagated point by point according to the path diagram to obtain the unwrapped differential interferogram.

[0011] Preferably, the weight setting principle includes: giving priority to disentanglement of highly coherent pixels with high coherence and close distance.

[0012] Preferably, the atmospheric delay phase correction method in the spatial domain is used to perform atmospheric delay phase correction on the unwrapped differential interferogram to obtain the atmospherically corrected unwrapped interferogram, which includes: Selecting coherent points with coherence greater than a preset threshold in the unwrapped differential interferogram set, establishing a linear model of atmospheric delay phase and elevation through least squares fitting, constructing an interferometric atmospheric delay phase map based on the linear model, and removing the interferometric atmospheric delay phase map from the unwrapped differential interferogram set to obtain an unwrapped phase map set after elevation correction; The unwrapped phase atlas after elevation correction is divided into blocks according to a preset size and subjected to FFT transformation. A Gaussian low-pass filter is used to estimate the long-wavelength atmospheric delay phase, which is then removed from the unwrapped phase atlas after elevation correction to obtain the atmospherically corrected unwrapped interferometry atlas.

[0013] Preferably, the time series deformation calculation is performed using the atmospherically corrected unwrapped interferogram to obtain the average annual settlement rate and the time series cumulative deformation results of the linear transportation corridor area, including: The minimum norm criterion and singular value decomposition method are used to invert the time series deformation matrix of the atmospherically corrected unwrapped interferogram to obtain the true surface deformation rate and time series cumulative deformation results of the target area.

[0014] Preferably, the extracting of spatial sedimentation gradient field information using the Sobel operator includes: Calculate the east-west gradient through the Sobel X convolution kernel; Calculate the north-south gradient using the Sobel Y convolution kernel; Solve for the gradient modulus and convert it to mm / m units.

[0015] Preferably, the hidden danger identification standard adopts a dual index threshold method, including pixel areas where the absolute value of the sedimentation rate is greater than the sedimentation rate threshold, and pixel areas where the modulus of the sedimentation spatial gradient vector field is greater than the modulus threshold.

[0016] Preferably, obtaining the extraction and identification results of roadbed settlement and landslide hazards in the target area includes: Pixels that meet any or both of the dual indicator thresholds are marked as hidden danger feature points; Vectorization is performed on the hidden danger feature points to generate polygonal vector hidden danger area patches, and the roadbed settlement and landslide hidden danger extraction and identification results of the target area are obtained.

[0017] In a second aspect, the present invention provides a system for identifying roadbed settlement and landslide hazards along a traffic corridor, comprising: The spaceborne satellite SAR data reading and data conversion module obtains the time-series spaceborne SAR SLC dataset of the target area and converts the format of the time-series spaceborne SAR SLC dataset to obtain the time-series spaceborne SAR SLC image set; The interferometric network construction module generates an interferometric pair network based on the time baseline threshold according to the time series spaceborne SAR SLC image set. The interferometric pair network contains multiple groups of interferometric pairs. Phase optimization processing module, which performs registration and differential interferometry processing on the interference pairs in the interference pair network and generates filtered differential interferometry atlas; The phase unwrapping module uses the minimum cost flow method based on irregular grids to perform differential interference phase unwrapping on the filtered differential interference atlas to obtain the unwrapped differential interference atlas; The atmospheric correction module uses the spatial domain atmospheric delay phase correction method to perform atmospheric delay phase correction on the unwrapped differential interferogram set to obtain the atmospheric corrected unwrapped interferogram set; The time series deformation inversion module uses the atmospherically corrected unwrapped interferogram to calculate the time series deformation, obtaining the annual average settlement rate and time series cumulative deformation results of the linear transportation corridor area; The hidden danger identification module uses the Sobel operator to extract the spatial settlement gradient field information of the linear traffic corridor area based on the average annual settlement rate of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output to obtain the roadbed settlement and landslide hidden danger extraction and identification results of the target area.

[0018] In a third aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon, which implements the above method when executed by a processor.

[0019] The above technical solution includes at least the following beneficial effects: based on spaceborne SAR SLC data, the traditional long-time series InSAR surface settlement and landslide monitoring technology is improved. In the implementation of long linear projects, roadbed settlement and landslide hazard identification systems along transportation corridors, it can significantly overcome the constraints of deformation extraction such as seasonal incoherence and macro-atmospheric disturbances. In addition, a dual-indicator early warning mechanism based on the annual average settlement rate and the annual settlement spatial gradient field is proposed, which can realize the rapid identification and extraction of geological settlement hazards and landslide threats in long linear transportation corridors with a span of thousands of kilometers in remote areas, thereby improving the efficiency of road area detection. In addition, the present application has the advantages of low cost, wide monitoring range, flexible monitoring area and high monitoring applicability. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] The above and other features, advantages and aspects of the embodiments of the present disclosure will become more apparent with reference to the following detailed description in conjunction with the accompanying drawings. The accompanying drawings are provided for a better understanding of the present disclosure and do not constitute a limitation of the present disclosure. In the accompanying drawings, the same or similar reference numerals represent the same or similar elements, among which: Figure 1 Schematic diagram of a flow chart of an embodiment of the present invention; Figure 2 Schematic diagram of the principle of an embodiment of the present invention; Figure 3 An optical image of a target area according to an embodiment of the present invention; Figure 4 This is a schematic diagram of the annual average sedimentation rate results of the target area in an embodiment of the present invention; Figure 5 This is a schematic diagram of the spatial gradient field results of the average annual settlement in the target area according to an embodiment of the present invention; Figure 6 This is a schematic diagram of the identification results of roadbed settlement and landslide hazards in the target area of ​​an embodiment of the present invention. DETAILED DESCRIPTION

[0021] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative efforts shall fall within the scope of protection of the present invention.

[0022] Reference Figure 1 and Figure 2 The embodiment of the present invention provides a method for identifying roadbed settlement and landslide hazards along a traffic corridor, comprising: Reference Figure 3, obtain the time-series spaceborne SAR SLC dataset of the target area, and convert the format of the time-series spaceborne SAR SLC dataset to obtain the time-series spaceborne SAR SLC image set; According to the time series spaceborne SAR SLC image set, based on the time baseline threshold, an interferometric pair network is generated, which contains multiple groups of interferometric pairs. Perform registration and differential interferometry processing on the interference pairs in the interference pair network and generate filtered differential interferometry atlas; The minimum cost flow method based on irregular grid is used to perform differential interferometry phase unwrapping on the filtered differential interferometry atlas to obtain the unwrapped differential interferometry atlas. The atmospheric delay phase correction method in spatial domain is used to correct the atmospheric delay phase of the unwrapped differential interferogram to obtain the atmospheric-corrected unwrapped interferogram. The atmospherically corrected unwrapped interferograms were used to calculate the time series deformation, and the annual average settlement rate and the time series cumulative deformation results of the linear transportation corridor area were obtained. Based on the average annual settlement rate of the linear traffic corridor area, the Sobel operator is used to extract the spatial settlement gradient field information of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output, and the roadbed settlement and landslide hidden danger extraction and identification results of the target area are obtained.

[0023] Preferably, obtaining a time-series spaceborne SAR SLC dataset of the target area and performing format conversion on the time-series spaceborne SAR SLC dataset to obtain a time-series spaceborne SAR SLC image set includes: The dedicated SAR image data storage format is converted to a universal binary encoding format. Based on the metafile information, SAR parameter conversion, data radiometric calibration, and thermal noise removal are performed sequentially to obtain the universal format SLC data and imaging parameter files required for time-series InSAR processing. The converted SLC data and imaging parameter files are stored in the processing platform.

[0024] Spaceborne SAR SLC imagery boasts high resolution and a vast observation range. A single satellite pass can monitor deformation over tens of thousands of square kilometers with a resolution of 10 meters, thereby improving the identification of geological subsidence hazards and landslide threats. Spaceborne SAR satellites typically use fixed orbits for repeated Earth observations at fixed intervals, completing a global land surface observation cycle every 10 to 28 days. This observation program is unaffected by factors such as rugged terrain and adverse climate conditions in the observation area, ensuring high adaptability to various monitoring environments.

[0025] SLC data, or single-look complex data, is one of the most primitive and core image formats used in synthetic aperture radar (SAR). Each pixel is complex, without multi-looking processing, preserving its original spatial resolution. This data stores both the amplitude and phase information of the ground surface corresponding to each pixel at any given moment. The phase information contains information related to surface deformation, forming the foundation of InSAR surface deformation monitoring data.

[0026] The SAR image header file contains key metadata, including geometry, orbital parameters, radar parameters, time information, and processing information. These parameters play a crucial role in InSAR processing, determining whether the SAR image can be correctly interpreted and whether operations such as SAR image registration, InSAR interferometry, unwrapping, and deformation transformation can be performed.

[0027] Preferably, an interference pair network is generated based on a time baseline threshold according to a time-series spaceborne SAR SLC image set, wherein the interference pair network includes multiple groups of interference pairs including: The SBAS network method was used to set the time baseline range and spatial baseline threshold to conduct the initial screening of the data set; Based on the principle of the shortest time baseline, the data after the initial screening were screened twice to obtain the interference pairs that met the conditions; According to the principle of minimizing the weighted sum of spatial vertical baseline, time baseline and Doppler centroid frequency difference, a public main image is selected and an interference pair network is generated.

[0028] In InSAR technology, the core of the SBAS (Small Baseline Subset) network method is to screen out high-coherence interference pairs by constraining the time baseline, spatial baseline, and Doppler baseline thresholds to balance deformation monitoring accuracy and data processing efficiency.

[0029] The purpose of the time baseline constraint is to reduce temporal decoherence caused by surface changes, such as vegetation growth and freeze-thaw cycles, to ensure that the interferometric phase primarily reflects surface deformation rather than environmental noise. The longer the time baseline ΔT, the greater the change in the surface scatterer characteristics, leading to phase decoherence.

[0030] The purpose of the spatial baseline constraint is to reduce geometric decorrelation and reduce the terrain phase error caused by the difference in satellite viewing angles. The larger the spatial vertical baseline B⊥ is, the more significant the geometric distortion of the same ground object point in the two images will be.

[0031] The purpose of the Doppler center frequency constraint is to avoid the registration error caused by the azimuth spectrum offset. The Doppler center frequency difference should not be greater than one resolution unit.

[0032] In this exemplary embodiment, seasonal freeze-thaw activity along a certain route in a certain region causes waterlogging along the road, leading to loss of coherence in InSAR measurements. To fully account for the region's seasonal freeze-thaw activity, the SBAS temporal baseline was set to 0–100 days, with a maximum spatial baseline threshold of 150 meters. Among the 16 potential interferometer pairs that met these baselines, the pairs were sorted in ascending order by the temporal baseline threshold, and the eight pairs with the shortest temporal baselines were selected as the preferred interferometer pairs.

[0033] Time-series InSAR analysis technology analyzes multiple time-series SAR images covering the same area. During the processing, one image must be selected as the common master image, and the remaining images are registered and sampled into the master image space. In an exemplary embodiment, the optimal common master image is selected based on the principle of minimizing the sum of the temporal baseline, spatial baseline, and Doppler centroid frequency baseline of the time-series interferometer pair. The specific selection method can be transformed into a problem of finding the optimal solution to the correlation coefficient, constructing a mathematical model for these multiple factors: Constructing correlation coefficients , is the correlation coefficient value when assuming that the mth image is the public master image: ; Among them, Bc, Tc, and fc are the critical values ​​of the spatial vertical baseline, time baseline, and Doppler centroid frequency difference, respectively. 、 、 are the absolute values ​​of the spatial vertical baseline, time baseline, and Doppler centroid frequency difference between the k-th image and the m-th image, respectively. L is the total number of images after screening. Function c is the relationship function, and α, β, and θ are the exponential factors of the three factors. The expression of function c is as follows: ; Where x is the independent variable of function c and a is the parameter; The image corresponding to the smallest correlation coefficient value is selected as the public main image.

[0034] Preferably, performing registration and differential interferometry processing on the interference pairs in the interference pair network and generating a filtered differential interferometry atlas comprises: The interferometer pair is registered using orbit baseline assisted registration, cross-correlation registration and enhanced spectral diversity methods; The phase difference between the registered interference pairs is calculated using the complex conjugate multiplication method, and the flat ground phase and terrain phase signals are removed from the phase difference signal to obtain the original differential interferogram set; Calculate the interference pair correlation coefficient atlas; According to the interference pair correlation coefficient atlas, a filter is used to perform differential interferogram adaptive filtering on the original differential interferogram atlas to generate a filtered differential interferogram atlas.

[0035] Preferably, the time-series SAR image set registration method includes track baseline auxiliary registration, intensity cross-correlation registration and enhanced spectral diversity registration.

[0036] Orbital baseline assisted registration is an important first step in the registration process. Its purpose is to perform coarse registration using orbital imaging parameters before sub-pixel precision registration, thereby reducing the initial geometric deviation between the primary and secondary images and improving the efficiency and stability of subsequent cross-correlation or spectral methods. In an exemplary embodiment, the orbital information of the primary and secondary images, including position and velocity vectors, is first read, and the satellite position difference corresponding to the observation time difference between the primary and secondary images is calculated to obtain the spatial baseline vector , and decomposed into vertical baselines and parallel baselines Then, combined with DEM or assumed flat ground height, the secondary image pixels are projected into the primary image coordinate system using the imaging geometry relationship to obtain a rough registration coordinate transformation field. The secondary image is then resampled into the primary image geometric frame using the bilinear interpolation method. The imaging geometry is related to the geodesic projection. is the main image satellite position vector, is the secondary image satellite position vector, is the ground point coordinate vector, and is the distance from the primary and secondary images to the ground point. Then the corresponding position of the secondary image pixel in the primary image is calculated through geometric projection.

[0037] Intensity Cross-Correlation (ICC) is used for coarse registration between the primary and secondary SLC images. It mainly uses a small area of ​​the primary image, that is, the template window, to slide in the secondary image and calculate the cross-correlation coefficient between the two windows. The position with the maximum correlation is the best matching point. Sub-pixel interpolation, such as quadratic surface fitting, is used to obtain the sub-pixel offset. The calculation formula is as follows: Set the main image window to , the secondary image is , the mutual correlation coefficient is: in and is the mean value of the pixels in the window.

[0038] Enhanced Spectral Diversity (ESD) registration is used for precise registration of SLC images. It uses the overlapping spectra between images acquired continuously by the same satellite in the same orbit, such as the TOPS mode of Sentinel-1, to estimate the azimuth direction, that is, the offset along the orbital direction, by calculating the interferometric phase difference between adjacent sub-strips (bursts). Specifically, for the primary and secondary images, the azimuth spectra of their overlapping areas are extracted, and the cross-correlation of the two spectra is calculated: Assume that the interference phase between the two sub-aperture images is: In this formula, and are the spectra of the primary image and the secondary image in the overlapping spectrum region; The fine drift is obtained by phase gradient analysis: in is the wavelength, is the slope distance, is the satellite speed, is the angle of incidence, is the azimuth distance.

[0039] The phase difference between the registered interference pairs is calculated using the complex conjugate multiplication method, and the flat bottom phase and terrain phase signals are removed from the phase difference signal.

[0040] After the two images are registered, they are conjugate multiplied to generate an interference pattern. The interference phase value of the resolution unit is: Since the SAR system is a time-range measurement, the observed value It is related to the slant range from the radar wave to the ground unit. In addition, it is also affected by the scattering characteristics of the ground resolution unit, including surface humidity, roughness, complex dielectric properties, etc. Assuming that the surface scattering characteristics remain unchanged during the two imaging, the interference phase can be expressed as: The difference in slant range between the radar wave and the ground target during the two imaging The main contributions are as follows: distance difference caused by the reference ellipsoid , that is, the flat land effect; the distance difference caused by the terrain undulation ;Distance difference caused by surface deformation during imaging ;Signal propagation delay caused by atmospheric disturbances ; Phase of system thermal noise, interferogram speckle noise, etc. ; Substituting the above components into, we get + + in is the interference phase, is the flat-Earth phase caused by the reference ellipsoid, is the terrain phase caused by ground undulation, is the LOS deformation phase caused by surface displacement during the two imaging periods, is the delayed phase caused by atmospheric inhomogeneity when acquiring two SAR images, is the random noise phase.

[0041] Flat Earth Phase It is the phase contribution caused by the curvature of the earth and the geometry of the satellite orbit, and has nothing to do with the terrain. After generating the interferogram, the flat earth phase must be removed first for subsequent processing. For each pixel in the interferogram, its flat earth phase It can be expressed as: in is the radar wavelength, is the vertical baseline, is the slant range from the satellite to the target, is the angle of incidence.

[0042] Terrain phase is the phase contribution caused by surface elevation. In differential interferometry, an external digital elevation model (DEM) is used to simulate and remove the terrain phase to obtain the deformation phase. The calculation formula is: in The actual elevation of each pixel provided by the DEM, It is the ellipsoidal elevation based on the WGS84 ellipsoid.

[0043] Differential interferometry phase after topographic phase difference processing At this point, only the deformation phase, atmospheric phase, and noise phase are retained. The adaptive filtering of differential interferograms using filters based on the interferometric coherence coefficient atlas includes the circular periodic mean Boxcar filtering method and the spectral-based adaptive filtering Goldstein method. This embodiment of the present invention utilizes the Goldstein frequency-domain adaptive filter to suppress noise by adjusting the amplitude of the interferometric phase spectrum. Its core concept is to perform strong filtering in areas of low coherence and weak filtering in areas of high coherence.

[0044] Coherence It is an indicator that measures the similarity of pixels at the same position in two SAR images and is used to evaluate the interferometric phase quality. Its value range is between 0 (completely incoherent) and 1 (completely coherent). Pixel window, coherence is calculated as follows: In this formula, and are the complex values ​​of the primary and secondary images, is the number of pixels in the sliding window, Represents the conjugate operation.

[0045] After calculating the coherence Then, Goldstein frequency domain adaptive filtering is performed. First, the interference pattern is divided into blocks. In this exemplary embodiment, 32x32 pixels are used for block division. Then, a two-dimensional fast Fourier transform (FFT) is performed on each small block to obtain the spectrum. : is the fast Fourier transform of the interferogram, is the filtering parameter, It is a custom constant and is set to 1 in this example.

[0046] Then build the filter function : Apply a Goldstein frequency adaptive filter to the interferogram : The method of using the minimum cost flow method based on irregular grid to perform differential interference phase unwrapping processing on the filtered differential interference atlas, and obtaining the unwrapped differential interference atlas includes: A coherence threshold is preset. In the embodiment of the present invention, the coherence threshold is set to 0.3. Highly coherent pixels with a coherence greater than the coherence threshold are selected from the filtered differential interferogram to form a scattered point combination. Taking the high coherence pixels as nodes, the high coherence pixels are connected according to preset connection conditions to form a triangular network diagram; The preset connection conditions are: is the distance between any two highly coherent pixels, is the coherence between any two highly coherent pixels, is the maximum distance between all highly coherent pixels, is the minimum value of coherence among all high coherence pixels; Calculate the phase difference on the edge: for each edge in the graph , calculate the winding phase difference in and are the phases of the two highly coherent pixels on this edge; Finding the smallest integer , so that the unwrapping phase difference is minimized: define the integer number of cycles , the phase difference after unwrapping is Seek the smoothest phase difference, that is Seeking The smallest side.

[0047] Building a minimum spanning tree: based on edge weights Set the untangle order. In this example, Then, points with high coherence and close distance are preferentially untangled, and a path graph is constructed using the minimum production tree.

[0048] Phase unwrapping: Select the most stable point of the time-series image signal as the reference point, and propagate the phase value point by point along the path: Preferably, the atmospheric delay phase correction method in the spatial domain is used to perform atmospheric delay phase correction on the unwrapped differential interferogram to obtain the atmospherically corrected unwrapped interferogram, which includes: Selecting coherent points with coherence greater than a preset threshold in the unwrapped differential interferogram set, establishing a linear model of atmospheric delay phase and elevation through least squares fitting, constructing an interferometric atmospheric delay phase map based on the linear model, and removing the interferometric atmospheric delay phase map from the unwrapped differential interferogram set to obtain an unwrapped phase map set after elevation correction; The unwrapped phase atlas after elevation correction is divided into blocks according to a preset size and subjected to FFT transformation, a Gaussian low-pass filter is used to estimate the long-wavelength atmospheric delay phase, and the long-wavelength atmospheric delay phase is removed from the unwrapped phase atlas after elevation correction to obtain an unwrapped interference atlas after atmospheric correction. In an embodiment of the present invention, the atmospheric delay phase correction in the spatial domain is performed on the unwrapped differential interference atlas, mainly utilizing the atmospheric delay phase elevation linear fitting and the spatial domain long-wavelength filtering method.

[0049] The atmospheric delay phase height linear fitting mainly considers that the atmospheric troposphere delay is mainly affected by the height. Assuming that the atmospheric delay has a certain linear relationship with the height in space, a function model between the atmospheric delay and the height is established by regression fitting, and then the component is removed from the interferogram. The specific implementation steps are as follows, wherein the formula used in the step and the parameters used in the formula are not universal with the parameters in the foregoing formula, and the explanation in the step is for reference: Select high-coherence backbone points: select high-quality points with a coherence greater than 0.3 from the coherence coefficient map of the interference pair; Construct a linear model: for all high-quality points satisfying the coherence greater than 0.3 Extract the same-name pixel unwrapping phase value to perform least square fitting to obtain the linear coefficient between the phase and the height . Let represent the atmospheric delay phase of the pixel position, represent the terrain height, represent the high linear fitting coefficient, and the basic linear model is wherein is a residual term, reflecting the non-height-related atmospheric error and other noise.

[0050] Use the regression coefficient to construct the interferometric atmospheric delay phase map of the target area, and subtract from the interferogram to obtain the differential interferogram set and the unwrapping phase set after the height linear calibration. The specific calculation formula is as follows: The spatial domain long wavelength filtering method suppresses the atmospheric delay phase, mainly considering that the atmospheric delay generally has spatial long wavelength and low frequency characteristics, while the surface deformation such as fault slip has high frequency variation characteristics. The interferogram is high-pass filtered or low-pass fitted to remove, retain the high frequency deformation, and remove the low frequency atmosphere. The specific implementation path is as follows, wherein the formula used in the step and the parameters used in the formula are not universal with the parameters in the foregoing formula, and the explanation in the step is for reference: Fully consider that the long wavelength atmospheric characteristics generally have spatial low frequency characteristics, and the wavelength is generally greater than 20 km. The deformation signal of the target area generally has a spatial size of less than 5 km. Therefore, the unwrapping phase set after the linear height phase correction is divided into blocks according to the size of 20 km*20 km, and the fast Fourier transform is used to convert to the frequency domain, and the calculation formula is as follows: A Gaussian low-pass filter is used, and the cut-off wavelength is set to 20 km. The calculation formula is as follows: in , is the cut-off wavelength.

[0051] Applying Gaussian low-pass filter to spectrum to estimate long-wavelength atmospheric delay phase By inverse fast Fourier transform Transform back to the image domain and get: Removing the long-wavelength atmospheric delay term from the interferometric phase pattern yields a clean unwrapped differential interferogram: Preferably, refer to Figure 4 The above-mentioned atmospheric-corrected unwrapped interferogram is used to calculate the time series deformation, and the annual average settlement rate and time series cumulative deformation results of the linear transportation corridor area are obtained, including: The minimum norm criterion and singular value decomposition method are used to invert the time series deformation matrix of the unwrapped differential interferogram after atmospheric correction to obtain the true annual average deposition rate and time series cumulative deformation results of the target area.

[0052] In the above embodiment of the present invention, the formula and parameters used in this step are not universal with the parameters in the aforementioned formula, and the explanation in this step shall prevail: Assume that the time series of the study area is obtained The time series of SAR images with coverage is: Select one of them as the super main image, and all images are composed according to the small baseline principle The scene interference pattern is: by is the reference time, any time Relative to Differential phase at time is an unknown number, and the differential interferometry phase is obtained during data processing. is the observed quantity. When only the phase change phase is retained in the interference pattern, the moment 、 The corresponding pixel in the differential interference pattern The phase value can be expressed as: is the radar wavelength; is the deformation along the radar line of sight between time A and B, 、 are the phase changes caused by the deformation.

[0053] Since SBAS-InSAR technology calculates the deformation of each pixel in the differential interferogram pixel by pixel in the time series, we take a certain pixel as an example to introduce the solution model of SBAS-InSAR technology. Suppose the phase of a certain pixel in the unwrapped differential interferogram group of the small baseline set matrix forms a vector : middle, is the phase value relative to the reference point. The time series corresponding to the main and auxiliary images are: If the main and auxiliary images are arranged in time sequence, the phase in the differential interferogram is expressed as: The equation shown is the one containing An unknown number equation, can be simplified to for The coefficient matrix of each row has a non-zero value, which is an interference pair. and , the other elements in the matrix are zero, then If all interference pairs belong to the same subbaseline set, then the matrix Rank , at this time the rank is full, and the least squares method is used to solve it: When SBAS consists of multiple subsets, the matrix Rank deficiency, For a singular matrix, the singular value decomposition method is used to solve it: Pair Matrix The singular value decomposition is: Where, for orthogonal matrix of order; The diagonal elements of are singular values ; for The orthogonal matrix of order. At this time, the least squares norm solution is solved by the generalized inverse matrix: Where, .

[0054] In order to obtain a solution that conforms to physical meaning, the solution of phase is converted into the solution of deformation rate, and the parameter vector to be solved is: Can be transformed and simplified to: for Matrix element , and the other elements are 0. Perform singular value decomposition to solve the deformation rate in each time period , and then calculate the time series shape variable based on this.

[0055] Preferably, refer to Figure 5 and Figure 6 In an embodiment of the present invention, based on the average annual settlement rate of the linear traffic corridor area, the Sobel operator is used to extract the spatial settlement gradient field information of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output to obtain the roadbed settlement and landslide hidden danger extraction and identification results of the target area. The formula and parameters used in this step are not universal with the parameters in the aforementioned formula, and the explanation in this step shall prevail.

[0056] The calculation and analysis of InSAR deformation phase gradients is an important method for identifying key geological elements such as deformation boundaries, fault activity, and landslide slip zones. In an InSAR deformation rate map, the pixel value of each pixel represents the deformation along the radar line of sight, and its spatial gradient represents the intensity of the spatial variation of the deformation rate. The phase gradient calculation method is: Further solving the gradient modulus is: The specific implementation steps are as follows: The east-west sedimentation gradient field of the annual average sedimentation rate is calculated using a 3*3 convolution kernel. The specific calculation kernel is: The north-south direction sedimentation gradient field of the annual average sedimentation rate is calculated using a 3*3 convolution kernel. The specific calculation kernel is: Use the grid calculator to calculate the modulus length of the east-west and north-south gradients of the annual sedimentation rate results. The specific calculation formula is as follows: ; in 、 These are the east-west gradient and north-south gradient of the annual average sedimentation rate results, respectively.

[0057] The Sobel operator output is in mm / pixel, which needs to be converted to mm / m by dividing by the pixel size to obtain the normalized annual sedimentation rate gradient norm. For this example, the mapping resolution for this region is 30 m. The specific conversion formula in the raster calculator is: In order to accurately identify the significant abnormal areas in the region, this scheme uses a strict double-threshold criterion to perform binary element determination on the above deformation field: Sedimentation rate criterion: Extract all pixel areas where the absolute value of the sedimentation rate exceeds a preset threshold of 15 mm / year. This threshold specifically screens out areas with significantly abnormal sedimentation rates and is a key indicator for identifying severe sedimentation or potential instability.

[0058] Spatial gradient criterion: All pixel regions where the modulus of the subsidence spatial gradient vector field exceeds a preset threshold of 50 mm / 100 m are simultaneously extracted. This high gradient threshold is intended to capture areas where the subsidence rate changes dramatically over a short distance. This is often a typical surface response characteristic of instability phenomena such as local geological structural variations (such as faults and weak interlayers), slope shear deformation, or potential landslide back-edge extension.

[0059] Binarization and vectorization: Pixel locations that meet any or both of the above criteria are marked as "hazard feature points" and assigned a value of 1. All other areas are marked as background and assigned a value of 0, thereby generating a preliminary binary mask map of the hazard areas. Furthermore, this binary raster data is vectorized to generate potential hazard patches (Polygon Features) with clear geographic boundaries. This vectorization process not only clearly defines the hazard area but also facilitates subsequent spatial overlay analysis, risk grading, and precise positioning of engineering remediation measures.

[0060] The output of this step, i.e., the hidden danger area map extracted and vectorized based on the dual physical quantity threshold criterion, constitutes the core input data layer for subsequent detailed risk analysis, cause inference, and early warning decision-making.

[0061] In a preferred embodiment, the present invention provides a system for identifying roadbed settlement and landslide hazards along a traffic corridor, comprising: The spaceborne satellite SAR data reading and data conversion module obtains the time-series spaceborne SAR SLC dataset of the target area and converts the format of the time-series spaceborne SAR SLC dataset to obtain the time-series spaceborne SAR SLC image set; The interferometric network construction module generates an interferometric pair network based on the time baseline threshold according to the time series spaceborne SAR SLC image set. The interferometric pair network contains multiple groups of interferometric pairs. Phase optimization processing module, which performs registration and differential interferometry processing on the interference pairs in the interference pair network and generates filtered differential interferometry atlas; The phase unwrapping module uses the minimum cost flow method based on irregular grids to perform differential interference phase unwrapping on the filtered differential interference atlas to obtain the unwrapped differential interference atlas; The atmospheric correction module uses the spatial domain atmospheric delay phase correction method to perform atmospheric delay phase correction on the unwrapped differential interferogram set to obtain the atmospheric corrected unwrapped interferogram set; The time series deformation inversion module uses the atmospherically corrected unwrapped interferogram to calculate the time series deformation, obtaining the annual average settlement rate and time series cumulative deformation results of the linear transportation corridor area; The hidden danger identification module uses the Sobel operator to extract the spatial settlement gradient field information of the linear traffic corridor area based on the average annual settlement rate of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output to obtain the roadbed settlement and landslide hidden danger extraction and identification results of the target area.

[0062] In a preferred embodiment, the present invention discloses a computer-readable storage medium having a computer program stored thereon, which implements the above method when executed by a processor.

[0063] Those skilled in the art will understand that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by those skilled in the art in the art to which the present invention pertains. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with those in the context of the prior art and, unless specifically defined, will not be interpreted in an idealized or overly formal sense.

[0064] For the method embodiments, the description is made in a series of action combinations for simplicity and clarity, but those skilled in the art should know that the present application is not limited to the action sequence described, because according to the present application, certain steps can be performed in other sequences or at the same time. In addition, those skilled in the art should know that the embodiments described in the specification are all preferred embodiments, and the actions involved are not necessarily essential to the present application.

[0065] From the above description of the embodiments, those skilled in the art can clearly understand that the present application can be implemented by means of software and the necessary general hardware platform. Based on such understanding, the technical solutions of the present application can be embodied in the form of a software product, which can be stored in a storage medium, such as a ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments of the present application.

[0066] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for identifying roadbed settlement and landslide hazards along a traffic corridor, characterized in that: include: Obtain the time-series spaceborne SAR SLC dataset of the target area and convert the format of the time-series spaceborne SAR SLC dataset to obtain the time-series spaceborne SAR SLC image set; According to the time series spaceborne SAR SLC image set, based on the time baseline threshold, an interferometric pair network is generated, which contains multiple groups of interferometric pairs. Perform registration and differential interferometry processing on the interference pairs in the interference pair network and generate filtered differential interferometry atlas; The minimum cost flow method based on irregular grid is used to perform differential interferometry phase unwrapping on the filtered differential interferometry atlas to obtain the unwrapped differential interferometry atlas. The atmospheric delay phase correction method in spatial domain is used to correct the atmospheric delay phase of the unwrapped differential interferogram to obtain the atmospheric-corrected unwrapped interferogram. The atmospherically corrected unwrapped interferograms were used to calculate the time series deformation, and the annual average settlement rate and the time series cumulative deformation results of the linear transportation corridor area were obtained. Based on the average annual settlement rate of the linear traffic corridor area, the Sobel operator is used to extract the spatial settlement gradient field information of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output, and the roadbed settlement and landslide hidden danger extraction and identification results of the target area are obtained.

2. The method according to claim 1, characterized in that The generating of the interference pair network based on the time baseline threshold comprises: The SBAS network method was used to set the time baseline range and spatial baseline threshold to conduct the initial screening of the data set; Based on the principle of the shortest time baseline, the data after the initial screening were screened twice to obtain the interference pairs that met the conditions; According to the principle of minimizing the weighted sum of spatial vertical baseline, time baseline and Doppler centroid frequency difference, a public main image is selected and an interference pair network is generated.

3. The method according to claim 2, characterized in that The selecting of the public main image includes: Constructing correlation coefficients , is the correlation coefficient when the mth image is the public main image, ； Among them, Bc, Tc, and fc are the critical values ​​of the spatial vertical baseline, time baseline, and Doppler centroid frequency difference, respectively. 、 、 are the absolute values ​​of the spatial vertical baseline, time baseline, and Doppler centroid frequency difference between the k-th image and the m-th image, respectively. L is the total number of images after screening. Function c is the relationship function, and α, β, and θ are the exponential factors of the three factors. The expression of function c is as follows: ； Where x is the independent variable and a is the parameter; The image corresponding to the minimum correlation coefficient value is selected as the public main image.

4. The method according to claim 1, wherein The performing of registration and differential interferometry processing and generating a filtered differential interferometry atlas comprises: The interferometer pair is registered using orbit baseline assisted registration, cross-correlation registration and enhanced spectral diversity methods; The phase difference between the registered interference pairs is calculated using the complex conjugate multiplication method, and the flat ground phase and terrain phase signals are removed from the phase difference signal to obtain the original differential interferogram set; Calculate the interference pair correlation coefficient atlas; According to the interference pair correlation coefficient atlas, a filter is used to perform differential interferogram adaptive filtering on the original differential interferogram atlas to generate a filtered differential interferogram atlas.

5. The method according to claim 4, characterized in that The calculation of the interference pair correlation coefficient atlas includes: Calculate the coherence of each interferometer pair ; ， in, is the complex value of the i-th pixel of the common main image in the interferometer, is the conjugate operation value of the complex value of the i-th pixel of the secondary image in the interferometer, is the number of pixels of the sliding window.

6. The method according to claim 4, characterized in that The method of using a filter to perform adaptive differential interferogram filtering on the original differential interferogram set includes a circular periodic mean Boxcar filtering method or a spectrum feature-based adaptive filtering Goldstein method.

7. The method according to claim 1, characterized in that The method of using the minimum cost flow method based on irregular grid to perform differential interference phase unwrapping processing on the filtered differential interference atlas, and obtaining the unwrapped differential interference atlas includes: A coherence threshold is preset, and high-coherence pixels with a coherence greater than the coherence threshold are selected from the filtered differential interference pattern set to form a scattered point combination; Taking the high coherence pixels as nodes, the high coherence pixels are connected according to preset connection conditions to form a triangular network diagram; The preset connection conditions are: ， is the distance between any two highly coherent pixels, is the coherence between any two highly coherent pixels, is the maximum distance between all highly coherent pixels, is the minimum value of coherence among all high coherence pixels; Compute the winding phase difference for each edge in a triangulated network: ， in and are the phases of two highly coherent pixels on one edge respectively; Finding the minimum integer number of cycles , so that the unwrapping phase difference is minimized, including: ； in, is the unwrapping phase difference; The untangling order is set based on the weight of each edge in the preset triangular network, and the path graph is constructed using the minimum production tree; The point with the most stable imaging signal in the filtered differential interferogram is selected as the reference point, and the phase value is propagated point by point according to the path diagram to obtain the unwrapped differential interferogram.

8. The method according to claim 7, characterized in that The weight setting principle includes: giving priority to disentanglement of highly coherent pixels with high coherence and close distance.

9. The method according to claim 1, characterized in that The atmospheric delay phase correction method in the spatial domain is used to perform atmospheric delay phase correction on the unwrapped differential interferogram to obtain the atmospherically corrected unwrapped interferogram. Selecting coherent points with coherence greater than a preset threshold in the unwrapped differential interferogram set, establishing a linear model of atmospheric delay phase and elevation through least squares fitting, constructing an interferometric atmospheric delay phase map based on the linear model, and removing the interferometric atmospheric delay phase map from the unwrapped differential interferogram set to obtain an unwrapped phase map set after elevation correction; The unwrapped phase atlas after elevation correction is divided into blocks according to a preset size and subjected to FFT transformation. A Gaussian low-pass filter is used to estimate the long-wavelength atmospheric delay phase, which is then removed from the unwrapped phase atlas after elevation correction to obtain the atmospherically corrected unwrapped interferometry atlas.

10. The method according to claim 1, characterized in that The above method uses the atmospheric-corrected unwrapped interferogram to calculate the time series deformation, and obtains the annual average settlement rate and time series cumulative deformation results of the linear transportation corridor area: The minimum norm criterion and singular value decomposition method are used to invert the time series deformation matrix of the atmospherically corrected unwrapped interferogram to obtain the true surface deformation rate and time series cumulative deformation results of the target area.

11. The method according to claim 1, characterized in that The method of extracting spatial sedimentation gradient field information using the Sobel operator includes: Calculate the east-west gradient through the Sobel X convolution kernel; Calculate the north-south gradient using the Sobel Y convolution kernel; Solve for the gradient modulus and convert it to mm / m units.

12. The method according to claim 1, characterized in that The hidden danger identification standard adopts a dual indicator threshold method, including: The pixel area where the absolute value of the sedimentation rate is greater than the sedimentation rate threshold, and the pixel area where the modulus of the sedimentation spatial gradient vector field is greater than the modulus threshold.

13. The method according to claim 1, wherein The extraction and identification results of roadbed settlement and landslide hazards in the target area include: Pixels that meet any or both of the dual indicator thresholds are marked as hidden danger feature points; Vectorization is performed on the hidden danger feature points to generate polygonal vector hidden danger area patches, and the roadbed settlement and landslide hidden danger extraction and identification results of the target area are obtained.

14. A roadbed settlement and landslide hazard identification system along a traffic corridor, characterized in that: include: The spaceborne satellite SAR data reading and data conversion module obtains the time-series spaceborne SAR SLC dataset of the target area and converts the format of the time-series spaceborne SAR SLC dataset to obtain the time-series spaceborne SAR SLC image set; The interferometric network construction module generates an interferometric pair network based on the time baseline threshold according to the time series spaceborne SAR SLC image set. The interferometric pair network contains multiple groups of interferometric pairs. Phase optimization processing module, which performs registration and differential interferometry processing on the interference pairs in the interference pair network and generates filtered differential interferometry atlas; The phase unwrapping module uses the minimum cost flow method based on irregular grids to perform differential interference phase unwrapping on the filtered differential interference atlas to obtain the unwrapped differential interference atlas; The atmospheric correction module uses the spatial domain atmospheric delay phase correction method to perform atmospheric delay phase correction on the unwrapped differential interferogram set to obtain the atmospheric corrected unwrapped interferogram set; The time series deformation inversion module uses the atmospherically corrected unwrapped interferogram to calculate the time series deformation, obtaining the annual average settlement rate and the time series cumulative deformation results of the linear transportation corridor area; The hidden danger identification module uses the Sobel operator to extract the spatial settlement gradient field information of the linear traffic corridor area based on the average annual settlement rate of the linear traffic corridor area. After binarization and vectorization, the settlement rate and spatial settlement gradient field information that meet the hidden danger identification criteria are output to obtain the roadbed settlement and landslide hidden danger extraction and identification results of the target area.

15. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the method according to any one of claims 1 to 13 is implemented.

Images