A method and system for identifying subgrade settlement and landslide hazards along a traffic corridor

Abstract

The application discloses a roadbed settlement and landslide hidden danger identification method and system along a traffic corridor, acquires a time-series satellite-borne data set of a target area, and carries out format conversion on the data format; generates an interference pair network based on a time baseline threshold; generates a filtered differential interferogram set; carries out differential interference phase unwrapping processing on the filtered differential interferogram set, and acquires an unwrapped differential interferogram set; carries out atmospheric delay phase correction on the unwrapped differential interferogram set, and obtains an atmospheric-corrected unwrapped interferogram set; calculates the linear traffic corridor area annual settlement rate and the time series cumulative deformation result; extracts the linear traffic corridor area spatial settlement gradient field information based on the linear traffic corridor area annual settlement rate, outputs the settlement rate and the spatial settlement gradient field information meeting the hidden danger discrimination standard, and obtains the roadbed settlement and landslide hidden danger extraction and identification result of the target area, thereby improving the road area detection efficiency.

Description

TECHNICAL FIELD

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

[0002] With the rapid development of society, the demand for infrastructure construction in remote areas is increasing, especially in plateau areas. The construction of permafrost roadbeds is of great significance to ensuring the smoothness of the transportation network. The stability of existing permafrost roadbeds is severely affected by seasonal freezing and thawing deformation, which not only threatens the safe operation of the road, but also has an adverse impact on the ecological environment. Traditional monitoring methods, such as leveling, GNSS, etc., have problems such as high cost, low efficiency, limited spatial coverage, etc., and are difficult to meet the full-cycle monitoring needs of thousands of kilometers of traffic corridors.

[0003] Although existing InSAR technology can achieve large-scale monitoring, it has the following limitations for long linear projects: vegetation coverage along the line, summer rainfall, and glacial melting leading to loss of radar signal coherence; long-wavelength atmospheric delay phase masks local mm-level deformation signals; lack of deformation hazard anomaly identification and analysis methods for the spatial characteristics of linear projects. SUMMARY

[0004] To achieve the above object and other related objects, the present application discloses a method for identifying roadbed settlement and landslide hazards along a traffic corridor, comprising:

[0005] Obtain a time-series spaceborne SAR SLC data set for a target area, and perform format conversion on the time-series spaceborne SAR SLC data set to obtain a time-series spaceborne SAR SLC image set;

[0006] Based on the time-series spaceborne SAR SLC image set, generate an interference pair network based on a time baseline threshold, the interference pair network containing multiple interference pairs;

[0007] Perform registration and differential interference processing on the interference pairs in the interference pair network, and generate a filtered differential interference image set;

[0008] Using an irregular grid-based minimum cost flow method, perform differential interference phase unwrapping processing on the filtered differential interference image set to obtain an unwrapped differential interference image set;

[0009] Using a spatial domain atmospheric delay phase correction method, perform atmospheric delay phase correction on the unwrapped differential interference image set to obtain an atmospheric-corrected unwrapped interference image set;

[0010] Using the atmospheric-corrected unwrapped interference image set, perform time-series deformation calculation to obtain the annual average settlement rate and the time series cumulative deformation result of the linear traffic corridor area;

[0011] Based on the linear traffic corridor area annual average settlement rate, the Sobel operator is used to extract the linear traffic corridor area spatial settlement gradient field information, and the settlement rate and spatial settlement gradient field information meeting the hidden danger discrimination standard are output after binarization and vectorization, so as to obtain the subgrade settlement and landslide hidden danger extraction and identification result of the target area.

[0012] Preferably, the time baseline threshold is used to generate an interference pair network, including:

[0013] The SBAS network networking method is adopted, the time baseline range and the spatial baseline threshold are set, and the data set data is initially screened;

[0014] The data after the initial screening is secondarily screened according to the principle of the shortest time baseline, and the interference pair meeting the condition is obtained;

[0015] According to the principle of the minimum weighted sum of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference, the common master image is selected, and the interference pair network is generated.

[0016] Preferably, the common master image includes:

[0017] The correlation coefficient is constructed , is the correlation coefficient when the mth image is the common master image,

[0018] ;

[0019] Wherein Bc, Tc and fc are the critical values of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference in turn, and 、 、 are the absolute values of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference of the kth image and the mth image in turn, L is the total number of images after screening, the function c is a relationship function, and alpha, beta and theta are the exponential factors of the three factors, wherein the expression of the function c is as follows:

[0020] ;

[0021] Wherein x is the independent variable, and a is the parameter;

[0022] The image corresponding to the minimum correlation coefficient value is selected as the common master image.

[0023] Preferably, the registration and differential interference processing are carried out, and the filtered differential interference diagram set is generated, including:

[0024] The track baseline auxiliary registration, the cross-correlation registration and the enhanced spectrum diversity method are used to complete the registration of the interference pair;

[0025] The phase difference between the registered interferometric pairs is obtained by using the complex conjugate multiplication method, and the flatland phase and terrain phase signals are removed from the phase difference signal to obtain the original differential interferogram.

[0026] Calculate the correlation coefficient atlas of interferometric pairs;

[0027] Based on the correlation coefficient set of interferometric pairs, an adaptive filtering of the original differential interferogram is performed using a filter to generate a filtered differential interferogram set.

[0028] Preferably, the atlas for calculating the correlation coefficients of interferometric pairs includes:

[0029] Calculate the coherence of each interference pair. ;

[0030] ,

[0031] in, and These are the complex values ​​of the common principal image and its sub-image in each interferometric pair, respectively. The number of pixels in the sliding window. This indicates the conjugate operation.

[0032] Preferably, the adaptive filtering of the original differential interferogram set using filters includes the Boxcar filtering method with circular periodic mean or the Goldstein method based on spectral characteristics.

[0033] Preferably, the step of using a minimum cost flow method based on irregular grids to perform differential interferometric phase unwrapping processing on the filtered differential interferometric atlas to obtain the unwrapped differential interferometric atlas includes:

[0034] A pre-set coherence threshold is used to select highly coherent pixels with coherence greater than the coherence threshold from the filtered differential interferogram set, forming a scatter plot.

[0035] Using the high-coherence pixels as nodes, the high-coherence pixels are connected according to preset connection conditions to form a triangular network diagram;

[0036] The preset connection conditions are as follows:

[0037] ,

[0038] Let be the distance between any two highly coherent pixels. For the coherence between any two highly coherent pixels, This represents the maximum distance between all highly coherent pixels. It represents the minimum coherence among all highly coherent pixels;

[0039] For each edge in the triangle network, calculate its phase difference of unwrapping:

[0040] ,

[0041] Wherein and are the phases of the two high-coherence pixels of the edge;

[0042] Find the minimum integer number of periods , so that the unwrapping phase difference is minimized, including:

[0043] ;

[0044] Wherein, is the unwrapping phase difference;

[0045] Based on the preset weight of each edge in the triangle network, set the unwrapping order, and use the minimum production tree to construct a path graph;

[0046] Select the point with the most stable imaging signal in the filtered differential interferogram set as the reference point, propagate the phase value point by point according to the path graph, and obtain the unwrapped differential interferogram set.

[0047] Preferably, the setting principle of the weight includes: preferentially unwrapping high-coherence pixels with high coherence and short distance.

[0048] Preferably, the atmospheric delay phase correction method in the spatial domain is used to correct the unwrapped differential interferogram set to obtain an atmospheric-corrected unwrapped interferogram set, including:

[0049] Select the coherent points with coherence greater than a preset threshold in the unwrapped differential interferogram set, establish a linear model of atmospheric delay phase and height by least squares fitting, construct an interferometric atmospheric delay phase map based on the linear model, and remove the interferometric atmospheric delay phase map from the unwrapped differential interferogram set to obtain a height-corrected unwrapped phase map set;

[0050] The height-corrected unwrapped phase map set 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 height-corrected unwrapped phase map set to obtain an atmospheric-corrected unwrapped interferogram set.

[0051] Preferably, the atmospheric-corrected unwrapped interferogram set is used to calculate the time series deformation to obtain the annual average subsidence rate and the time series cumulative deformation result of the linear traffic corridor region, including:

[0052] The atmospheric corrected unwrapped interferogram set is inversely calculated by using a minimum norm criterion and a singular value decomposition method to obtain real ground surface deformation rate and time sequence cumulative deformation result of the target region.

[0053] Preferably, the extraction of the spatial settlement gradient field information by using the Sobel operator comprises:

[0054] The east-west direction gradient is calculated by using a Sobel X convolution kernel.

[0055] The south-north direction gradient is calculated by using a Sobel Y convolution kernel.

[0056] The gradient module is solved and converted into mm / m units.

[0057] Preferably, the hidden danger discrimination criterion adopts a double-index threshold method, including a pixel region with an absolute value of the settlement rate greater than a settlement rate threshold value and a pixel region with a spatial gradient vector field module length greater than a module length threshold value.

[0058] Preferably, the extraction and identification result of the subgrade settlement and landslide hidden danger of the target region comprises:

[0059] The pixel meeting any one or both of the double-index threshold values is marked as a hidden danger feature point.

[0060] The hidden danger feature point is subjected to vectorization processing to generate a polygon vector hidden danger area patch, and the extraction and identification result of the subgrade settlement and landslide hidden danger of the target region is obtained.

[0061] In a second aspect, the present application provides a subgrade settlement and landslide hidden danger identification system along a traffic corridor, comprising:

[0062] A satellite-borne satellite SAR data reading and data conversion module acquires a time sequence satellite-borne SAR SLC data set of a target region, and converts the time sequence satellite-borne SAR SLC data set to obtain a time sequence satellite-borne SAR SLC image set.

[0063] An interference network construction module generates an interference pair network based on a time baseline threshold value according to the time sequence satellite-borne SAR SLC image set, and the interference pair network contains multiple groups of interference pairs.

[0064] A phase optimization processing module performs registration and differential interference processing on the interference pairs in the interference pair network, and generates a filtered differential interferogram set.

[0065] A phase unwrapping module adopts an irregular grid-based minimum cost flow method to perform differential interference phase unwrapping processing on the filtered differential interferogram set to obtain an unwrapped differential interferogram set.

[0066] The atmospheric correction module corrects the unwrapped differential interferogram set by using a spatial domain atmospheric delay phase correction method to obtain an atmospheric corrected unwrapped interferogram set.

[0067] The time series deformation inversion module calculates time series deformation by using the atmospheric corrected unwrapped interferogram set to obtain linear traffic corridor area annual average settlement rate and time series cumulative deformation results.

[0068] The hidden danger identification module extracts spatial settlement gradient field information of the linear traffic corridor area by using a Sobel operator based on the linear traffic corridor area annual average settlement rate, and outputs settlement rate and spatial settlement gradient field information meeting hidden danger identification standards through binarization and vectorization to obtain subgrade settlement and landslide hidden danger extraction and identification results of the target area.

[0069] In a third aspect, the present application provides a computer readable storage medium having a computer program stored thereon, the program being executed by a processor to implement the above method.

[0070] The above technical solution has at least the following beneficial effects: based on space-borne SAR SLC data, the traditional long-time series InSAR ground settlement and landslide monitoring technology is improved, in the development of long linear engineering, traffic corridor along the subgrade settlement and landslide hidden danger identification system, the constraints of deformation extraction such as seasonal decorrelation and macro atmospheric disturbance can be significantly overcome, in addition, through the double-index early warning mechanism based on annual average settlement rate and annual settlement spatial gradient field, the rapid identification and extraction of geological settlement disasters and landslide stress of long linear traffic corridor with a span of thousands of kilometers in remote road area can be realized, the road detection efficiency is improved, in addition, the present application has the advantages of low cost, wide monitoring range, flexible monitoring area and high monitoring applicability. BRIEF DESCRIPTION OF DRAWINGS

[0071] The above and other features, advantages, and aspects of the present disclosure will become more apparent with reference to the following detailed description in conjunction with the accompanying drawings that illustrate the present disclosure by way of example. The drawings provided herein are for illustrative purposes only and constitute part of the detailed description. In the drawings, the same or similar reference signs refer to the same or similar elements, wherein:

[0072] Figure 1 The flowchart of the embodiment of the present application is shown;

[0073] Figure 2 The principle diagram of the embodiment of the present application is shown;

[0074] Figure 3 The optical image of the target area of the embodiment of the present application is shown;

[0075] Figure 4 The annual average settlement rate result diagram of the target area of the embodiment of the present application is shown;

[0076] Figure 5 A schematic diagram of the spatial gradient field of the annual average settlement of the target area in the embodiment of the present application is shown in the figure.

[0077] Figure 6 A schematic diagram of the roadbed settlement and landslide hazard identification result of the target area in the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0078] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0079] Referring to Figure 1 and Figure 2 , the embodiment of the present application provides a roadbed settlement and landslide hazard identification method along a traffic corridor, comprising:

[0080] Referring to Figure 3 , a time-series spaceborne SAR SLC data set of a target area is obtained, and the time-series spaceborne SAR SLC data set is format-converted to obtain a time-series spaceborne SAR SLC image set;

[0081] According to the time-series spaceborne SAR SLC image set, an interference pair network is generated based on a time baseline threshold, and the interference pair network contains multiple groups of interference pairs;

[0082] The interference pairs in the interference pair network are subjected to registration and differential interference processing, and a filtered differential interference image set is generated;

[0083] An irregular grid-based minimum cost flow method is used to perform differential interference phase unwrapping processing on the filtered differential interference image set to obtain an unwrapped differential interference image set;

[0084] A spatial domain atmospheric delay phase correction method is used to perform atmospheric delay phase correction on the unwrapped differential interference image set to obtain an atmospheric-corrected unwrapped interference image set;

[0085] The atmospheric-corrected unwrapped interference image set is used to calculate a time-series deformation to obtain a linear traffic corridor area annual average settlement rate and a time sequence cumulative deformation result;

[0086] Based on the linear traffic corridor area annual average settlement rate, the Sobel operator is used to extract the linear traffic corridor area spatial settlement gradient field information, and the settlement rate and spatial settlement gradient field information meeting the hidden danger discrimination standard are output after binarization and vectorization, to obtain the subgrade settlement and landslide hidden danger extraction and identification result of the target area.

[0087] Preferably, a time-series spaceborne SAR SLC dataset of the target area is acquired, and the time-series spaceborne SAR SLC dataset is format-converted to obtain a time-series spaceborne SAR SLC image set, which includes:

[0088] The dedicated SAR image data storage format is converted into a general binary encoding format, and SAR parameter conversion, data radiation scaling, thermal noise removal and the like are sequentially performed according to the metafile information, to obtain general format SLC data and imaging parameter files necessary for carrying out time-series InSAR processing. The converted SLC data and imaging parameter files are stored in the processing platform.

[0089] The spaceborne SAR SLC image has the technical advantages of high resolution and ultra-large observation range, and can complete the deformation observation of the ground surface with a resolution of 10 m within a range of tens of thousands of square kilometers in one satellite pass, thereby improving the identification efficiency of geological subsidence disasters and landslide stress. The spaceborne SAR satellite usually adopts a fixed orbit for shooting, and can repeatedly observe the earth with a fixed period, and can complete one global land surface observation in 10-28 days. The observation plan is not affected by factors such as rugged terrain and harsh weather environment of the observation area, and has high monitoring environmental applicability.

[0090] SLC data, i.e. single-view complex data, is one of the most original and core image formats in synthetic aperture radar, each pixel of which is a complex number without multi-view (Multi-looking) processing, and the original spatial resolution is maintained. The data stores the imaging amplitude information and phase information of the ground surface corresponding to each pixel at any imaging moment, wherein the phase stores information related to the ground surface deformation, and is the basis of InSAR ground surface deformation monitoring data.

[0091] The header file (Header File) of the SAR image contains key files of metadata such as geometric information, orbit parameters, radar parameters, time information, processing information and the like. These parameters play a crucial role in InSAR processing, and determine whether the SAR image can be correctly interpreted, whether SAR image registration, InSAR interference, unwrapping and deformation conversion and the like can be performed.

[0092] Preferably, according to the time-series spaceborne SAR SLC image set, an interference pair network is generated based on a time baseline threshold, and the interference pair network contains multiple groups of interference pairs, including:

[0093] The SBAS network method is adopted, the time baseline range and the spatial baseline threshold are set, and the data set data is initially screened;

[0094] The data after the initial screening is secondarily screened according to the principle that the time baseline is the shortest, and the interference pairs meeting the conditions are obtained.

[0095] According to the principle that the weighted sum of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference is the minimum, the common main image is selected, and the interference pair network is generated.

[0096] In the InSAR technology, the core of the SBAS (Small Baseline Subset) network method is to balance the deformation monitoring accuracy and the data processing efficiency by screening the high-coherent interference pair combination through the constraints of the time baseline, the spatial baseline and the Doppler baseline threshold.

[0097] The purpose of the time baseline constraint is to reduce the time de-coherence caused by the ground surface changes such as vegetation growth and freeze-thaw cycle, and to ensure that the interference phase mainly reflects the ground surface deformation rather than the environmental noise. The longer the time baseline ΔT is, the greater the change of the ground surface scatterer characteristics is, which leads to the phase de-coherence.

[0098] The purpose of the spatial baseline constraint is to reduce the geometric de-coherence and to reduce the terrain phase error caused by the difference in the satellite viewing angle. The larger the spatial vertical baseline B⊥ is, the more significant the geometric distortion of the same ground object point in the two images is.

[0099] The purpose of the Doppler centroid frequency constraint is to avoid the registration error caused by the azimuth spectrum shift. The Doppler centroid frequency difference is not greater than one resolution unit.

[0100] In the demonstrative embodiment, the water along the line of a certain area is mainly affected by the seasonal freeze-thaw, which leads to the InSAR measurement de-coherence. Considering the characteristics of the seasonal freeze-thaw activities in the certain area, the SBAS time baseline is set to 0-100 days, and the maximum threshold of the spatial baseline is 150 m. In the 16 potential interference pairs meeting the baseline, the interference pairs are sorted in ascending order according to the time baseline threshold, and the 8 interference pairs with the shortest time baseline are selected as the preferred interference pairs.

[0101] The time series InSAR analysis technology is developed for multiple time series SAR images covering the same area. In the processing process, one of them is selected as the common main image, and the remaining images are registered and sampled to the main image space. In the demonstrative embodiment, the optimal common main image is selected according to the principle that the sum of the time baseline, the spatial baseline and the Doppler centroid frequency baseline of the time series interference pair is the minimum. The specific selection method can be converted into an optimal solution problem of a correlation coefficient. The mathematical model is constructed for the above factors:

[0102] The correlation coefficient , is the correlation coefficient value when the m-th image is assumed to be the common master image:

[0103] ;

[0104] where Bc, Tc, fc are the critical values of the spatial vertical baseline, the temporal baseline and the Doppler centroid frequency difference respectively, and , , are the absolute values of the spatial vertical baseline, the temporal baseline and the Doppler centroid frequency difference of the k-th image and the m-th image respectively, L is the total number of images after screening, the function c is a relationship function, and a, b, q are the index factors of the three factors, wherein the expression of the function c is as follows:

[0105] ;

[0106] where x is the independent variable of the function c, and a is a parameter;

[0107] The image corresponding to the minimum correlation coefficient value is selected as the common master image.

[0108] Preferably, the registration and differential interference processing are carried out on the interference pairs in the interference pair network, and a filtered differential interference image set is generated, including:

[0109] The interference pairs are registered by using the orbit baseline assisted registration, the cross-correlation registration and the enhanced spectrum diversity method;

[0110] The phase difference between the registered interference pairs is calculated by using the complex conjugate multiplication method, and the flat phase and the terrain phase signals are removed from the phase difference signal to obtain an original differential interference image set;

[0111] A correlation coefficient image set of the interference pairs is calculated;

[0112] According to the correlation coefficient image set of the interference pairs, a filter is used to carry out adaptive filtering on the original differential interference image set to generate a filtered differential interference image set.

[0113] Preferably, the registration method of the time sequence SAR image set includes orbit baseline assisted registration, intensity cross-correlation registration and enhanced spectrum diversity registration.

[0114] The orbit baseline assisted registration is an important first step of the registration processing, and the purpose is to carry out coarse registration through orbit imaging parameters before sub-pixel accuracy registration, reduce the initial geometric deviation between the master image and the slave image, and improve the efficiency and stability of the subsequent cross-correlation or spectrum method. In the exemplary embodiment, first, the orbit information of the master image and the slave image is read, including the position and velocity vectors, the satellite position difference corresponding to the observation time difference of the master image and the slave image is calculated, and the spatial baseline vector is obtained and decomposed into vertical baseline and parallel baseline . Then, combined with DEM or assuming a flat height, the sub-image pixels are projected into the master image coordinate system using the imaging geometry, resulting in a rough registration coordinate transformation field. Then, the sub-image is resampled into the master image geometric framework using the bilinear interpolation method. The relationship between the imaging geometry and the geodetic projection is

[0115]

[0116]

[0117] is the master image satellite position vector, is the sub-image satellite position vector, is the ground point coordinate vector, and is the distance from the master and sub-images to the ground point. Then, the corresponding position of the sub-image pixel in the master image is calculated by geometric projection.

[0118] Intensity Cross-Correlation (ICC) is used for the coarse registration between the SLC master image and the sub-image. A small window, i.e., the template window, of the master image is used to slide in the sub-image, and the cross-correlation coefficient between the two windows is calculated. The maximum correlation position is the best matching point, and the sub-pixel level offset is obtained using sub-pixel interpolation, such as quadratic surface fitting. The calculation formula is as follows:

[0119] Let the master image window be and the sub-image be The cross-correlation coefficient is:

[0120]

[0121] where and are the mean values of the pixels in the window.

[0122] Enhanced Spectral Diversity (ESD) registration is used for the fine registration of SLC images. It uses the overlapping spectrum between images acquired by the same satellite on the same orbit, such as the TOPS mode of Sentinel-1, to estimate the azimuth direction, i.e., the along-track direction, by calculating the interference phase difference between adjacent bursts. Specifically, for the master image and the sub-image, the frequency spectrum in the azimuth direction of their overlapping regions is extracted, and the cross-correlation of the two frequency spectrums is calculated:

[0123] Let the interference phase between the two sub-aperture images be:

[0124]

[0125] In the formula, and are the spectra of the primary and secondary images in the overlapping spectral region, respectively;

[0126] The fine drift is obtained by phase gradient analysis:

[0127]

[0128] where is the wavelength, is the slant range, is the satellite velocity, is the incidence angle, is the azimuth distance.

[0129] The phase difference between the registered interferometric pairs is obtained by complex conjugate multiplication, and the flat phase and terrain phase signals are removed from the phase difference signal.

[0130] The two registered images are conjugate multiplied to generate an interferogram, and the interferometric phase value of the resolution cell is:

[0131]

[0132] Since the SAR system is a time ranging system, the observed value is related to the slant range of the radar wave to the ground cell, and is also affected by the scattering characteristics of the ground resolution cell, including the surface humidity, roughness, complex dielectric characteristics, etc. Assuming that the surface scattering characteristics remain unchanged during the two imaging times, the interferometric phase can be represented as:

[0133]

[0134] The difference in slant range between the radar wave to the ground target during the two imaging times is mainly contributed by the following aspects: the distance difference caused by the reference ellipsoid, i.e. the flat ground effect; the distance difference caused by the terrain undulation; the distance difference caused by the surface deformation during imaging; the signal propagation delay caused by atmospheric disturbance; the phase ; system thermal noise, speckle noise of the interferogram, etc. ; system thermal noise, speckle noise of the interferogram, etc. ; system thermal noise, speckle noise of the interferogram, etc. ; system thermal noise, speckle noise of the interferogram, etc. ; system thermal noise, speckle noise of the interferogram, etc.

[0135] + +

[0136] where is the interferometric phase, flat earth phase, topographic phase, phase due to deformation of the LOS vector during the two acquisitions, delay phase due to the non-uniformity of the atmosphere during the two SAR acquisitions, random noise phase.

[0137] flat earth phase is the phase contribution due to the curvature of the Earth and the satellite orbit geometry, and is independent of the topography. After the interferogram is generated, the flat earth phase is first removed in order to facilitate subsequent processing. For each pixel in the interferogram, its flat earth phase can be expressed as:

[0138]

[0139] where is the radar wavelength, is the vertical baseline, is the slant range from the satellite to the target, is the incidence angle.

[0140] The topographic phase is the phase contribution due to the surface elevation. In differential interferometry, the topographic phase is simulated and removed using an external digital elevation model (DEM), and the resulting phase is the deformation phase. The topographic phase is calculated as:

[0141]

[0142] where is the actual elevation of each pixel provided by the DEM, is the ellipsoidal elevation based on the WGS84 ellipsoid.

[0143] differential interferometric phase after topographic phase differencing At this point, only the deformation phase, the atmospheric phase and the noise phase remain. The use of filters to carry out adaptive filtering of the differential interferogram according to the interferometric coherence coefficient map set includes a circular period mean Boxcar filtering method and a Goldstein adaptive filtering method based on spectral characteristics. The embodiment of the present application uses a Goldstein frequency domain adaptive filter to suppress noise by adjusting the amplitude of the interferometric phase spectrum. The core idea is to perform strong filtering in areas with low coherence and weak filtering in areas with high coherence.

[0144] coherence This is an index that measures the pixel similarity at the same location between two SAR images, used to assess the quality of interferometric phase. Its value ranges from 0 (completely incoherent) to 1 (completely coherent). This exemplary example uses... For pixel windows, coherence is calculated as follows:

[0145]

[0146] In this formula, and These are the complex values ​​of the main image and the secondary image, respectively. The number of pixels in the sliding window. This indicates the conjugate operation.

[0147] After calculating coherence Next, Goldstein frequency domain adaptive filtering is performed. First, the interferogram is divided into blocks, using 32x32 pixels in this exemplary embodiment; then, a two-dimensional fast Fourier transform (FFT) is performed on each block to obtain the spectrum. :

[0148]

[0149]

[0150] For the fast Fourier transform of the interferogram, These are the filter parameters. This is a user-defined constant, set to 1 in this exemplary instance.

[0151] Then the filtering function is constructed. :

[0152]

[0153] Applying a Goldstein frequency adaptive filter to the interferogram :

[0154]

[0155] The method employing a minimum cost flow approach based on irregular grids performs differential interferometric phase unwrapping processing on the filtered differential interferometric atlas to obtain the unwrapped differential interferometric atlas, including:

[0156] A coherence threshold is preset. In this embodiment of the invention, the coherence threshold is set to 0.3. High coherence pixels with coherence greater than the coherence threshold are selected from the filtered differential interferometric dataset to form a scatter plot.

[0157] Taking the high-coherent pixels as nodes, connecting the high-coherent pixels according to a preset connection condition to form a triangular network diagram;

[0158] The preset connection condition is:

[0159]

[0160] is a distance between any two high-coherent pixels, is a coherence between any two high-coherent pixels, is a maximum value of distances between all high-coherent pixels, is a minimum value of coherences between all high-coherent pixels;

[0161] Calculating a phase difference on an edge: for each edge in the graph , calculating a wrapped phase difference

[0162]

[0163] wherein and are phases of the two high-coherent pixels of the edge;

[0164] Finding a minimum integer such that the unwrapped phase difference is minimum: defining the integer number of periods , the unwrapped phase difference is

[0165]

[0166] Seeking the smoothest phase difference, i.e., seeking the minimum edge in .

[0167] Constructing a minimum spanning tree: based on the weight of the edge Setting an unwrapping order, in the present exemplary embodiment, using

[0168]

[0169] Further, preferentially unwrapping points with high coherence and short distance, using the minimum spanning tree to construct a path diagram.

[0170] Phase unwrapping: selecting a point with the most stable integrated signal of the time-series image as a reference point, and propagating the phase value point by point along the path:

[0171]

[0172] Preferably, the method for correcting atmospheric delay phase in a spatial domain is used to correct the atmospheric delay phase of the unwrapped differential interferogram set, to obtain an atmospheric-corrected unwrapped interferogram set, comprising:

[0173] Select the coherent points in the unwrapped differential interferogram set whose coherence is greater than a preset threshold, establish a linear model of atmospheric delay phase and elevation by least square fitting, construct an interferometric atmospheric delay phase map based on the linear model, and remove the interferometric atmospheric delay phase map from the unwrapped differential interferogram set to obtain an elevation-corrected unwrapped phase map set;

[0174] The elevation-corrected unwrapped phase map set is divided into blocks according to a preset size and subjected to FFT transformation, a Gaussian low-pass filter is used to estimate a long-wavelength atmospheric delay phase, the long-wavelength atmospheric delay phase is removed from the elevation-corrected unwrapped phase map set to obtain an atmospheric-corrected unwrapped interferogram set. In the embodiment of the application, the spatial domain atmospheric delay phase correction of the unwrapped differential interferogram set is mainly implemented by using the atmospheric delay phase elevation linear fitting and the spatial domain long-wavelength filtering method.

[0175] The atmospheric delay phase elevation linear fitting mainly considers that the atmospheric tropospheric delay is mainly affected by the elevation. It is assumed that the atmospheric delay is in a certain linear relationship with the elevation in space, a function model between the atmospheric delay and the elevation is established by regression fitting, and then the component is removed from the interferogram. The specific execution 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 used as the reference:

[0176] 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;

[0177] Construct a linear model: for all high-quality points satisfying the coherence greater than 0.3 Extract the unwrapped phase values of the same-named pixels to perform least square fitting to obtain the linear coefficient between the phase and the elevation . Let represent the atmospheric delay phase of the pixel position, represent the terrain elevation, represent the high linear fitting coefficient, and the basic linear model is

[0178]

[0179] wherein is a residual term, reflecting the non-elevation-related atmospheric error and other noise.

[0180] The regression coefficient is used to construct the interferometric atmospheric delay phase map of the target area, and is subtracted from the interferogram to obtain the elevation-linearly calibrated differential interferogram set and the unwrapped phase map set. The specific calculation formula is as follows: ​​

[0181]

[0182] The spatial domain long-wavelength filtering method suppresses atmospheric delay phase. Mainly considering that atmospheric delay generally has spatial long-wavelength and low-frequency characteristics, and the ground deformation such as fault slip has high-frequency variation characteristics, the interferogram is high-pass filtered or low-pass fitted to subtract, retain high-frequency deformation, and eliminate low-frequency atmosphere. The specific implementation path is as follows, wherein the formula used in this step and the parameters used in the formula are not universal with the parameters in the foregoing formula. Please refer to the explanation in this step:

[0183] Full consideration is given to the fact that long-wavelength atmospheric characteristics generally have spatial low-frequency characteristics, and the wavelength is often greater than 20 km. The target area deformation signal generally has a spatial size of less than 5 km. Therefore, the linear elevation phase-corrected unwrapped phase set is divided into blocks according to the size of 20 km*20 km, and is converted to the frequency domain by using fast Fourier transform. The calculation formula is:

[0184]

[0185] A Gaussian low-pass filter is used, and the cut-off wavelength is set to 20 km. The calculation formula is:

[0186]

[0187] wherein , is the cut-off wavelength.

[0188] The Gaussian low-pass filter is applied to the spectrum diagram to estimate the long-wavelength atmospheric delay phase

[0189]

[0190] By inverse fast Fourier transform, the is converted back to the image domain to obtain:

[0191]

[0192] The long-wavelength atmospheric delay term is removed from the interferometric phase diagram to obtain a clean unwrapped differential interferogram:

[0193]

[0194] Preferably, with reference to Figure 4 , the linear traffic corridor area annual settlement rate and the time sequence cumulative deformation result are obtained by using the atmospheric-corrected unwrapped interferogram set to calculate the time sequence deformation.

[0195] The minimum norm criterion and singular value decomposition method are used for time sequence deformation matrix inversion of the atmospheric correction and unwrapped differential interferogram set, so that the real annual average subsidence rate and time sequence cumulative deformation variable result of the target region are obtained.

[0196] 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 used as a reference.

[0197] It is assumed that the time sequence of the study area is obtained The SAR image coverage is selected, and the image acquisition time sequence is:

[0198]

[0199] One of them is selected as a super master image, and all images are composed into a small baseline set interferogram, so that:

[0200]

[0201] Take as the reference time, and the difference phase of the target region at any time relative to the time is an unknown number, and the difference interference phase is obtained in the data processing process. When only the phase change phase is retained in the interferogram, the time , The phase value of the pixel in the corresponding difference interferogram can be expressed as:

[0202]

[0203] is the radar wavelength; is the deformation variable along the radar line of sight between times A and B, , respectively, are the change phase values caused by the deformation variable.

[0204] Since the SBAS-InSAR technology is to calculate the deformation of each pixel in the time sequence in the difference interferogram, the calculation model of the SBAS-InSAR technology is introduced by taking a pixel as an example. Let the phase of a pixel in the unwrapped difference interferogram group of the small baseline set matrix form a vector :

[0205]

[0206] In the formula, the phase value relative to the reference point is . The time sequences corresponding to the main and auxiliary images are:

[0207]

[0208] If the primary and secondary images are arranged in time sequence, the phase in the differential interferogram is represented as:

[0209]

[0210] The equation shown is a system of equations with unknowns and equations, which can be simplified as

[0211]

[0212] is a coefficient matrix of , and the non-zero value of each row is an interference pair, with and , and other elements in the matrix are zero, then

[0213]

[0214] If all interference pairs belong to the same sub-baseline set, the matrix has a rank of , and the rank is full at this time. The least squares method is used for solving:

[0215]

[0216] When the SBAS consists of multiple subsets, the matrix is rank-deficient, and is a singular matrix. At this time, the singular value decomposition method is used for solving:

[0217]

[0218] The singular value decomposition of the matrix has:

[0219]

[0220] In the formula, is an orthogonal matrix; the diagonal elements of are singular values; is an orthogonal matrix. At this time, the least squares norm solution through the generalized inverse matrix is:

[0221]

[0222] In the formula, .

[0223] In order to obtain a physically meaningful solution, the phase is solved into the deformation rate, and the parameter vector to be solved is:

[0224]

[0225] which can be transformed and simplified as:

[0226]

[0227] = 0 Matrix; matrix element , and other element values are 0. For singular value decomposition is performed, and the deformation rate in each time period is solved, and then the time sequence deformation variable is calculated according to the deformation rate.

[0228] Preferably, referring to Figure 5 and Figure 6 , in the embodiment of the present application, based on the linear traffic corridor area annual settlement rate, the Sobel operator is used to extract the linear traffic corridor area spatial settlement gradient field information, and the settlement rate and the spatial settlement gradient field information meeting the hidden danger discrimination standard are output through binarization and vectorization, so as to obtain the subgrade settlement and landslide hidden danger extraction and identification result of the target area. The formula used in this step and the parameters used in the formula are not universal with the parameters in the foregoing formula, and the explanation in this step is used as the standard.

[0229] The calculation and analysis of the InSAR deformation phase gradient is an important method for identifying key geological elements such as deformation boundary, fault activity and landslide slip zone. In the InSAR deformation rate graph, the pixel value of each pixel represents the deformation variable along the radar line-of-sight direction, and the spatial gradient represents the intensity of the change of the deformation rate in space. The phase gradient calculation method is:

[0230]

[0231] Further solving the gradient module length has:

[0232]

[0233] The specific implementation steps are as follows:

[0234] The 3*3 convolution kernel is used to calculate the east-west direction settlement gradient field of the annual settlement rate result, and the specific calculation kernel is:

[0235]

[0236] The 3*3 convolution kernel is used to calculate the north-south direction settlement gradient field of the annual settlement rate result, and the specific calculation kernel is:

[0237]

[0238] The east-west and north-south gradients of the annual settlement rate results are calculated using a grid calculator, and the specific calculation formula is as follows:

[0239]

[0240] The east-west and north-south gradients of the annual settlement rate results are calculated using a grid calculator, and the specific calculation formula is as follows:

[0241] The unit of the Sobel operator output is mm / pixel, which needs to be divided by the pixel size to convert to mm / m, to obtain the normalized annual settlement rate result gradient module length. The area example uses a mapping resolution of 30 m, so the specific conversion formula in the grid calculator is:

[0242]

[0243] To accurately identify the significant abnormal areas in the area, the present scheme uses a strict double-threshold criterion to perform binary element judgment on the above deformation field:

[0244] Settlement rate criterion: Extract all pixel regions with an absolute value of settlement rate greater than a preset threshold of 15 mm / year. This threshold specifically filters out areas with significantly abnormal settlement rates, which is a key indicator for identifying severe settlement or potential instability.

[0245] Spatial gradient criterion: Extract all pixel regions with a module length of the settlement spatial gradient vector field greater than a preset threshold of 50 mm / 100m. This high gradient threshold aims to capture areas where the settlement rate changes sharply over short distances, which is typically a typical surface response characteristic of local geological structure variation (such as faults, weak interlayers), slope shear deformation, or potential landslide rear edge tension, etc.

[0246] Binary and vectorization processing: Mark the pixel positions that meet the above criteria or both criteria as "hidden danger feature points", i.e. assign a value of 1, and mark the remaining areas as background, i.e. assign a value of 0, to generate a preliminary hidden danger area binary mask map. Further, the binary grid data is subjected to vectorization processing to generate potential hidden danger area polygons (Polygon Features) with clear geographical boundaries. This vectorization process not only clearly defines the scope of hidden dangers, but also facilitates subsequent spatial overlay analysis, risk level classification, and precise positioning of engineering treatment measures.

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

[0248] In one preferred embodiment, the present application provides a roadbed settlement and landslide hidden danger identification system along a traffic corridor, comprising:

[0249] The satellite-borne SAR data reading and data conversion module acquires a time-series satellite-borne SAR SLC dataset and performs format conversion on the time-series satellite-borne SAR SLC dataset to obtain a time-series satellite-borne SAR SLC image set.

[0250] The interference network construction module generates an interference pair network based on a time baseline threshold according to the time-series satellite-borne SAR SLC image set, and the interference pair network contains multiple groups of interference pairs.

[0251] The phase optimization processing module performs registration and differential interference processing on the interference pairs in the interference pair network and generates a filtered differential interference graph set.

[0252] The phase unwrapping module uses an irregular grid-based minimum cost flow method to perform differential interference phase unwrapping processing on the filtered differential interference graph set to obtain an unwrapped differential interference graph set.

[0253] The atmospheric correction module uses a spatial domain atmospheric delay phase correction method to perform atmospheric delay phase correction on the unwrapped differential interference graph set to obtain an atmospheric-corrected unwrapped interference graph set.

[0254] The time-series deformation inversion module uses the atmospheric-corrected unwrapped interference graph set to perform time-series deformation calculation to obtain the linear traffic corridor region annual settlement rate and the time sequence cumulative deformation result.

[0255] The hidden danger identification module extracts the linear traffic corridor region spatial settlement gradient field information using a Sobel operator based on the linear traffic corridor region annual settlement rate, performs binarization and vectorization, and outputs the settlement rate and spatial settlement gradient field information that meet the hidden danger discrimination standard to obtain the roadbed settlement and landslide hidden danger extraction and identification result of the target region.

[0256] In one preferred embodiment, the present application discloses a computer readable storage medium having a computer program stored thereon, which is executed by a processor to implement the above method.

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

[0258] For the sake of simplicity, the method embodiments are described as a series of actions. However, those skilled in the art should understand that the embodiments of the present invention are not limited to the described order of actions, because according to the embodiments of the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to the embodiments of the present invention.

[0259] As can be seen from the above description of the embodiments, those skilled in the art can clearly understand that this application can be implemented by means of software plus necessary general-purpose hardware platforms. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in various embodiments or some parts of the embodiments of this application.

[0260] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for identifying roadbed settlement and landslide hazards along a transportation corridor, characterized in that, The method comprises the following steps: acquire a time-series spaceborne SAR SLC dataset, and perform format conversion on the time-series spaceborne SAR SLC dataset to obtain a time-series spaceborne SAR SLC image set; generate an interferometric pair network based on a time baseline threshold according to the time-series spaceborne SAR SLC image set, and the interferometric pair network comprises a plurality of interferometric pairs; perform registration and differential interference processing on the interferometric pairs in the interferometric pair network, and generate a filtered differential interferogram set; perform differential interference phase unwrapping processing on the filtered differential interferogram set by using an irregular grid-based minimum cost flow method to obtain an unwrapped differential interferogram set; perform atmospheric delay phase correction on the unwrapped differential interferogram set by using a spatial domain atmospheric delay phase correction method to obtain an atmospheric-corrected unwrapped interferogram set; perform time-series deformation calculation by using the atmospheric-corrected unwrapped interferogram set to obtain an annual average settlement rate and a time series cumulative deformation result of the linear traffic corridor region; extract spatial settlement gradient field information of the linear traffic corridor region by using a Sobel operator based on the annual average settlement rate of the linear traffic corridor region, binarize and vectorize the settlement rate and the spatial settlement gradient field information to output settlement rate and spatial settlement gradient field information meeting a hidden danger discrimination standard, and obtain subgrade settlement and landslide hidden danger extraction and identification results of the target region; wherein the differential interference phase unwrapping processing on the filtered differential interferogram set by using the irregular grid-based minimum cost flow method to obtain the unwrapped differential interferogram set comprises: preset a coherence threshold, and select high-coherence pixels with coherence greater than the coherence threshold from the filtered differential interferogram set to form a scatter point combination; connect the high-coherence pixels according to a preset connection condition to form a triangular network graph, with the high-coherence pixels as nodes; wherein the preset connection condition is: ， is the distance between any two high-coherent pixels, is the coherence between any two high-coherent pixels, is the maximum distance between all high-coherent pixels, is the minimum coherence between all high-coherent pixels; calculate the wrapped phase difference of each edge in the triangular network graph; ， wherein and are the phases of the two high-coherence pixels of one edge, respectively; Finding the minimum number of integer periods minimizing the unwrapped phase difference, comprising: ； wherein is the unwrapped phase difference; set an unwrapping sequence based on the weight of each edge in the preset triangular network graph, and use a minimum production tree to construct a path graph; select a point with the most stable imaging signal in the filtered differential interferogram set as a reference point, and propagate the phase value point by point according to the path graph to obtain the unwrapped differential interferogram set; the atmospheric delay phase correction on the unwrapped differential interferogram set by using the spatial domain atmospheric delay phase correction method to obtain the atmospheric-corrected unwrapped interferogram set comprises: select coherent points with coherence greater than a preset threshold from the unwrapped differential interferogram set, establish a linear model of atmospheric delay phase and elevation by least square fitting, construct an interferometric atmospheric delay phase graph based on the linear model, and remove the interferometric atmospheric delay phase graph from the unwrapped differential interferogram set to obtain an elevation-corrected unwrapped phase graph set; block the elevation-corrected unwrapped phase graph set according to a preset size and perform FFT transformation, estimate long-wavelength atmospheric delay phase by using a Gaussian low-pass filter, remove the long-wavelength atmospheric delay phase from the elevation-corrected unwrapped phase graph set to obtain the atmospheric-corrected unwrapped interferogram set; The linear traffic corridor area annual settlement rate and the time series cumulative deformation result are obtained by using the atmospheric correction disentangled interferogram set for timing deformation calculation: The real ground deformation rate and time series cumulative deformation variable result of the target area are obtained by using the minimum norm criterion and singular value decomposition method for time series deformation matrix inversion of the atmospheric correction disentangled interferogram set.

2. The method of claim 1, wherein, The interference pair network is generated based on the time baseline threshold, including: The time baseline range and the spatial baseline threshold are set by using the SBAS network method, and the data set data is initially screened; The shortest time baseline is taken as the principle, and the data after the initial screening is secondarily screened to obtain the interference pairs meeting the conditions; The common master image is selected according to the minimum principle of the weighted sum of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference, and the interference pair network is generated.

3. The method of claim 2, wherein, The common master image includes: correlation coefficient , is the correlation coefficient when the mth image is the common master image, ； wherein Bc, Tc, fcare the critical values of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference, respectively, , , are the absolute values of the spatial vertical baseline, the time baseline and the Doppler centroid frequency difference of the kth image and the mth image, respectively, L is the total number of images after screening, c is a function, and a, β, θ are the index factors of the three factors, wherein the expression of the function c is as follows: ； Wherein x is the independent variable, and a is the parameter; The image corresponding to the minimum correlation coefficient value is selected as the common master image.

4. The method of claim 1, wherein, The registration and differential interference processing are carried out, and the filtered differential interferogram set is generated, including: The orbital baseline auxiliary registration, cross-correlation registration and enhanced spectral diversity method are used to complete the registration of the interference pairs; The phase difference between the registered interference pairs is calculated by using the complex conjugate multiplication method, and the flat phase and the terrain phase signals are removed from the phase difference signal to obtain the original differential interferogram set; The interference pair correlation coefficient set is calculated; According to the interference pair correlation coefficient set, the filter is used to carry out adaptive filtering on the original differential interferogram set to generate the filtered differential interferogram set.

5. The method of claim 4, wherein, The interference pair correlation coefficient set includes: calculating the coherence of each interference pair ; ， wherein, is a complex value of the i-th pixel of the primary image in the interference pair, is a conjugate operation value of a complex value of the i-th pixel of the secondary image in the interference pair, is a number of pixels of the sliding window.

6. The method of claim 4, wherein, The adaptive filtering of the original differential interferogram set using the filter includes the circular period mean Boxcar filtering method or the adaptive filtering Goldstein method based on the frequency spectrum characteristics.

7. The method of claim 1, wherein, The setting principle of the weight includes: preferentially disentangling the high-coherent pixels with high coherence and short distance.

8. The method of claim 1, wherein, The spatial settlement gradient field information is extracted by using the Sobel operator, including: The east-west direction gradient is calculated by using the Sobel X convolution kernel; The south-north direction gradient is calculated by using the Sobel Y convolution kernel; The gradient module length is solved and converted into mm / m unit.

9. The method of claim 1, wherein, The hidden danger discrimination standard adopts a double index threshold method, including: The pixel region with a settlement rate absolute value greater than the settlement rate threshold, and the pixel region with a settlement space gradient vector field module length greater than the module length threshold.

10. The method of claim 9, wherein, The roadbed settlement and landslide hidden danger extraction and identification result of the target area includes: The pixel meeting any one or both of the double index threshold values is marked as a hidden danger feature point; The hidden danger feature point is vectorized to generate a polygon vector hidden danger area patch, and the roadbed settlement and landslide hidden danger extraction and identification result of the target area is obtained.

11. A system for identifying subgrade settlement and landslide hazards along a transportation corridor, comprising: It includes: The satellite-borne SAR data reading and data conversion module obtains the time series satellite-borne SAR SLC data set of the target area, and converts the time series satellite-borne SAR SLC data set to obtain the time series satellite-borne SAR SLC image set; The interference network construction module generates an interference pair network based on a time baseline threshold according to a time sequence spaceborne SAR SLC image set, and the interference pair network contains multiple groups of interference pairs; The phase optimization processing module performs registration and differential interference processing on the interference pairs in the interference pair network, and generates a filtered differential interference image set; The phase unwrapping module uses an irregular grid-based minimum cost flow method to perform differential interference phase unwrapping processing on the filtered differential interference image set to obtain an unwrapped differential interference image set, including: presetting a coherence threshold, selecting high-coherence pixels with coherence greater than the coherence threshold from the filtered differential interference image set to form a scatter point combination; The high-coherence pixels are connected as nodes according to a preset connection condition to form a triangular network diagram; The preset connection condition is: ， is the distance between any two high-coherent pixels, is the coherence between any two high-coherent pixels, is the maximum distance between all high-coherent pixels, is the minimum coherence between all high-coherent pixels; Calculate the wrapped phase difference of each edge in the triangular network: ， wherein and are the phases of the two highly coherent pixels of one edge, respectively; Finding the minimum number of integer periods minimizing the unwrapped phase difference, comprising: ； wherein is the unwrapped phase difference; Determine the unwrapping order based on the weight of each edge in the preset triangular network, and use the minimum production tree to construct a path diagram; Select the most stable point in the imaging signal in the filtered differential interference image set as a reference point, and propagate the phase value point by point according to the path diagram to obtain the unwrapped differential interference image set; The atmospheric correction module uses a spatial domain atmospheric delay phase correction method to perform atmospheric delay phase correction on the unwrapped differential interference image set to obtain an atmospheric corrected unwrapped interference image set, including: selecting coherent points with coherence greater than a preset threshold from the unwrapped differential interference image set, establishing a linear model of atmospheric delay phase and elevation by least squares fitting, constructing an interference atmospheric delay phase diagram based on the linear model, and removing the interference atmospheric delay phase diagram from the unwrapped differential interference image set to obtain an elevation corrected unwrapped phase image set; Block the elevation corrected unwrapped phase image set according to the preset size and perform FFT transformation, estimate the long-wavelength atmospheric delay phase using a Gaussian low-pass filter, remove the long-wavelength atmospheric delay phase from the elevation corrected unwrapped phase image set, and obtain the atmospheric corrected unwrapped interference image set; The time sequence deformation inversion module uses the atmospheric corrected unwrapped interference image set to calculate the time sequence deformation to obtain the annual average settlement rate and the time sequence cumulative deformation result of the linear traffic corridor region, including: using the minimum norm criterion and singular value decomposition method to perform time sequence deformation matrix inversion on the atmospheric corrected unwrapped interference image set to obtain the true ground deformation rate and time sequence cumulative deformation result of the target region; The hidden danger identification module uses the Sobel operator to extract the spatial settlement gradient field information of the linear traffic corridor region based on the annual average settlement rate of the linear traffic corridor region, and outputs the settlement rate and spatial settlement gradient field information that meet the hidden danger discrimination standard after binarization and vectorization to obtain the subgrade settlement and landslide hidden danger extraction and identification result of the target region.

12. A computer readable storage medium having stored thereon a computer program, characterized in that The program is executed by the processor to implement the method of any one of claims 1-10.

Images