Early warning method and system for ground surface settlement along highway based on time sequence InSAR

Abstract

The invention belongs to the technical field of synthetic aperture radar measurement, and particularly relates to a road surface settlement early warning method and system based on a time sequence InSAR, and the method comprises the steps: obtaining a remote sensing image sequence and synchronous meteorological data of a monitoring region, generating a time sequence difference phase sequence and an air temperature difference value sequence of each candidate scattering point, calculating a relative phase temperature sensitivity coefficient between candidate scattering point pairs with a physical adjacent relationship, performing thermodynamic stripping on an original phase difference to obtain a corrected phase difference, constructing a space-time stability index in combination with a space physical distance, and taking the space-time stability index as a topological edge connection weight in an initial space topological graph to obtain a spatial topological graph; and determining an optimal connection path of phase unwrapping by executing a maximum spanning tree optimization algorithm and performing spatial deduction to obtain a real phase sequence after meteorological interference is stripped and calculate a settlement rate so as to trigger settlement early warning based on a judgment result. According to the invention, the capturing precision of the tiny settlement signal along the highway is improved.

Description

Technical Field

[0001] This invention relates to the field of synthetic aperture radar (SAR) measurement technology. More specifically, this invention relates to a method and system for early warning of land subsidence along highways based on time-series InSAR. Background Technology

[0002] With the rapid advancement of remote sensing technology, temporal interferometric synthetic aperture radar (InSAR) technology has been applied to wide-area land subsidence monitoring. Traditional InSAR methods primarily rely on permanent scatterer (PSS) techniques for deformation extraction. This approach acquires multiple InSAR images covering the target area to generate an interferogram sequence. Subsequently, an external digital elevation model is used to remove the flat and topographic phases, resulting in a differential interferogram. Based on this, a deviation index of pixel amplitude is calculated, and a fixed deviation threshold is set to screen for highly coherent PSS candidate points. For these candidate points, traditional techniques construct a spatial topology map, use spatial and temporal baselines to build a phase model, and assume that deformation exhibits a strictly linear evolution trend. Phase unwrapping is performed along the spatial topology map to recover the true phase, and finally, the deformation rate is extracted through spatiotemporal filtering.

[0003] However, in highway scenarios, asphalt or concrete pavements are affected by seasonal changes and exhibit significant thermal expansion and contraction physical characteristics, which manifest as strong nonlinear thermodynamic deformation in microwave interferometry phase. At the same time, high-frequency vehicle loads can cause transient settlement of the roadbed, leading to complex changes in pixel scattering characteristics. Traditional techniques use fixed thresholds for point selection, which easily excludes effective pavement key points with short-term coherence fluctuations caused by heavy vehicle obstruction as noise. More seriously, phase unwrapping models based on strict linear assumptions cannot fit the thermodynamic nonlinear evolution of pavement materials. They often treat the actual thermal expansion and contraction displacement as unwrapping errors and forcefully smooth them, resulting in large-area fractures in the phase integration path in stress concentration areas, which in turn causes serious closure deviations in the extracted local uneven settlement results.

[0004] The cumulative error caused by model mismatch can easily overwhelm the tiny displacement signals generated by early roadbed collapse. Furthermore, it cannot identify and avoid connected edges with loose geological structures or severe traffic fatigue during the spatial topology optimization process. As a result, the deformation monitoring accuracy for highways is insufficient to meet the needs of safe operation and maintenance. The phenomenon of missed or false alarms is serious, and it cannot provide reliable data support and scientific decision-making basis for the preventive maintenance of highway facilities. Summary of the Invention

[0005] To address the technical problems of existing InSAR technology, such as difficulty in removing nonlinear interference from highway thermodynamics and susceptibility to traffic loads leading to the submergence of minute settlement signals by errors and the failure of early warnings, this invention provides solutions in the following aspects.

[0006] In a first aspect, the present invention provides a method for early warning of land subsidence along highways based on time-series InSAR, comprising:

[0007] The remote sensing image sequence and synchronous meteorological data of the monitoring area are acquired, and the temporal differential phase sequence and air temperature difference sequence of candidate scattering points are generated through differential interferometry.

[0008] Based on the time-series differential phase sequence and the air temperature difference sequence, the relative phase temperature sensitivity coefficient between candidate scattering point pairs with physical adjacency is calculated, and the original phase difference between the candidate scattering point pairs is thermodynamically stripped using the relative phase temperature sensitivity coefficient to obtain the corrected phase difference that reflects the permanent deformation characteristics of the roadbed.

[0009] Based on the degree of discrete fluctuation of the corrected phase difference on the time axis, and combined with the spatial physical distance between the candidate scattering point pairs, a spatiotemporal stability index that takes into account both temporal immunity and spatial coherence is constructed.

[0010] Using the spatiotemporal stability index as the topological edge connection weight in the initial spatial topology graph, the optimal connection path for phase unwrapping is determined by executing the maximum spanning tree optimization algorithm. Spatial deduction is then performed according to the connection order of the optimal connection path to obtain the true phase sequence after removing meteorological interference. The settlement rate of candidate scattering points is determined based on the true phase sequence, and a settlement warning is triggered based on the determination result of the settlement rate.

[0011] Preferably, the step of generating the temporal differential phase sequence of candidate scattering points and the air temperature difference sequence through differential interferometry includes:

[0012] The remote sensing image sequence was registered and differentially interferometrically processed using an external digital elevation model to obtain a differential interferogram sequence. The amplitude deviation features of the pixels in the differential interferogram sequence were extracted, and each pixel was divided into multiple level categories using the natural breakpoint classification method. The pixel points in the level category with the smallest amplitude deviation feature were selected as candidate scattering points. The phase values ​​of the candidate scattering points at each observation time were extracted from the differential interferogram sequence and arranged in chronological order to obtain the temporal differential phase sequence of the corresponding candidate scattering points.

[0013] The air temperature difference at each observation time is obtained by subtracting the regional meteorological station temperature data from the pre-set reference temperature. The air temperature difference at each observation time is then arranged in chronological order to form an air temperature difference sequence.

[0014] Preferably, the relative phase temperature sensitivity coefficient satisfies the expression:

[0015] ;

[0016] In the formula, Indicates the first The candidate scattering point and the first The relative phase temperature sensitivity coefficient between candidate scattering points, the th The candidate scattering point and the first Each candidate scattering point has a physical adjacency relationship; Indicates the total number of observation times; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Represents the sequence of air temperature difference values. The air temperature difference at each observation time; It represents a very small positive number and is used to prevent the denominator from being zero in calculations.

[0017] Preferably, the corrected phase difference satisfies the expression:

[0018] ;

[0019] In the formula, Indicates the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering point and the first The relative phase temperature sensitivity coefficient between candidate scattering points; Represents the sequence of air temperature difference values. The air temperature difference at each observation time.

[0020] Preferably, the spatiotemporal stability index satisfies the expression:

[0021] ;

[0022] In the formula, Indicates the first The candidate scattering point and the first The spatiotemporal stability index between candidate scattering points, the th The candidate scattering point and the first Each candidate scattering point has a physical adjacency relationship; Indicates the first The candidate scattering point and the first The spatial physical distance between candidate scattering points; Indicates the distance to zero constant; This represents the function for extracting the maximum value. Represents an exponential function with the natural constant as its base; This represents the phase normalization constant; Indicates the total number of observation times; Indicates the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time; Indicates the first The candidate scattering point and the first The average corrected phase difference between the candidate scattering points at all observation times; Represents the absolute value symbol.

[0023] Preferably, before executing the maximum spanning tree optimization algorithm, the method further includes:

[0024] The basic skeleton network of the initial spatial topology graph is obtained based on the minimum spanning tree algorithm; the average and standard deviation of the spatial physical distance of all connected edges in the basic skeleton network are calculated; the sum of the average and three times the standard deviation is determined as the connected distance threshold; redundant topological edges with spatial physical distances greater than the connected distance threshold are removed from the initial spatial topology graph.

[0025] Preferably, the step of spatial extrapolating according to the connectivity order of the optimal connection path to obtain the true phase sequence after removing meteorological interference includes:

[0026] From all candidate scattering points, a candidate scattering point located on the periphery of the monitoring area with a stable geological structure is selected as a stable reference point, and its absolute phase reference is set to 0. Using the stable reference point as the starting root node, along the connection direction of the optimal connection path, two adjacent candidate scattering points on the optimal connection path are respectively taken as the starting candidate scattering point and the target candidate scattering point. Integer ambiguity is calculated for the temporal differential phase difference between the starting candidate scattering point and the target candidate scattering point to recover the relative phase difference. The relative phase difference is accumulated to the absolute phase value of the starting candidate scattering point to obtain the absolute phase value of the target candidate scattering point. The spatial extrapolation process is executed synchronously at all observation times to reconstruct the unwrapped phase sequence of each candidate scattering point.

[0027] The unwound phase sequence is filtered in the time dimension using a time high-pass filter window to separate the high-frequency time sequence; the high-frequency time sequence is filtered in the spatial dimension using a spatial low-pass filter window to extract the atmospheric delay component sequence; the atmospheric delay component sequence is subtracted from the unwound phase sequence to obtain the true phase sequence.

[0028] Preferably, determining the sedimentation rate of candidate scattering points based on the true phase sequence includes:

[0029] The true phase sequence is converted into a line-of-sight deformation sequence, and the line-of-sight deformation sequence is converted into a vertical time-series settlement sequence by using the geometric projection transformation relationship of the radar incident angle.

[0030] The sedimentation rate of candidate scattering points was obtained by fitting the time-series sedimentation sequence using a linear regression algorithm.

[0031] Preferably, triggering a settlement early warning based on the settlement rate determination includes:

[0032] The maximum absolute value of the instantaneous settlement rate of the stable reference point in all observation periods is calculated as the extreme value of the background settlement rate caused by the background disturbance of the natural environment at the stable reference point. The extreme value of the background settlement rate is numerically superimposed with the allowable settlement rate critical limit in the highway subgrade engineering structure specification to obtain the settlement rate warning threshold.

[0033] A settlement warning command is generated in response to a candidate scattering point having a settlement rate less than a negative value of the settlement rate warning threshold.

[0034] Secondly, the present invention provides a roadside surface subsidence early warning system based on time-series InSAR, including a processor and a memory. The memory stores computer program instructions, and when the computer program instructions are executed by the processor, the above-mentioned roadside surface subsidence early warning method based on time-series InSAR is implemented.

[0035] By adopting the above technical solution, the above-mentioned method for early warning of land subsidence along highways based on time-series InSAR is generated into a computer program and stored in a memory so that it can be loaded and executed by a processor. A terminal device can then be made based on the memory and processor for convenient use.

[0036] The beneficial effects of this invention are as follows: By deeply coupling synchronous meteorological data with remote sensing image interferometric phase, this invention accurately extracts the thermal response properties of road surface materials using the relative phase temperature sensitivity coefficient, achieving intrinsic stripping of nonlinear thermodynamic expansion disturbances and solving the deformation extraction deviation caused by the mismatch between traditional linear deformation models and the elastic physical properties of the road surface; This invention constructs a spatiotemporal stability index by comprehensively considering the phase time fluctuation and spatial distance attenuation law, measuring the structural stability under traffic load, and using this to drive the maximum spanning tree algorithm to autonomously plan the untangling path in the spatial topology network, ensuring that the phase integration derivation process actively avoids fragile areas with severe shear fatigue and loose geological structure, blocking the global error diffusion caused by local phase jumps, improving the acquisition accuracy and reliability of small settlement signals along the highway, and solving the technical problem of settlement signal submersion and early warning failure caused by environmental interference and model errors. Attached Figure Description

[0037] Figure 1 This is a flowchart illustrating the method for early warning of land subsidence along highways based on time-series InSAR in this invention.

[0038] Figure 2 This is a schematic diagram of the air temperature difference sequence;

[0039] Figure 3 A schematic diagram showing the phase difference before and after correction for candidate scattering point pairs with physical adjacency.

[0040] Figure 4 This is a diagram showing the final spatial topology and the optimal connection path.

[0041] Figure 5 This is a schematic diagram showing the distribution of sedimentation rates at candidate scattering points and the triggering of early warnings. Detailed Implementation

[0042] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0043] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0044] This invention discloses a method for early warning of land subsidence along highways based on time-series InSAR, referring to... Figure 1 This includes steps S1-S4:

[0045] S1. Acquire remote sensing image sequences and synchronous meteorological data of the monitoring area, and generate temporal differential phase sequences and air temperature difference sequences of candidate scattering points through differential interferometry.

[0046] It should be noted that although synthetic aperture radar images contain rich information on surface microwave scattering, the original echo signals are severely affected by topographic undulations and atmospheric composition, making it impossible to directly extract small subsidence data of the road surface. Therefore, this invention acquires remote sensing image sequences and synchronous meteorological data of the monitoring area, and generates temporal differential phase sequences and air temperature difference sequences of candidate scattering points at each observation time through differential interferometry, providing a data basis for subsequent removal of environmental thermodynamic interference.

[0047] Specifically, it acquires multi-phase synthetic aperture radar image sequences covering the monitoring area, as well as regional meteorological station temperature data at the corresponding satellite transit times.

[0048] The multi-phase synthetic aperture radar image sequence was registered and differentially interferometrically processed using an external digital elevation model to obtain a differential interferogram sequence. The amplitude deviation features of the pixels in the differential interferogram sequence were extracted. Each pixel was divided into multiple levels and categories using the natural breakpoint classification method. The pixel points in the level category with the smallest amplitude deviation feature were selected as candidate scattering points. The three-dimensional spatial coordinates of each candidate scattering point were extracted using the geocoding information of the external digital elevation model.

[0049] The phase values ​​of each candidate scattering point along the highway at each observation time are extracted from the differential interferogram sequence and used as the temporal differential phase of each candidate scattering point at each observation time. The temporal differential phases of the same candidate scattering point at all observation times are arranged in chronological order to obtain the temporal differential phase sequence of that candidate scattering point. Simultaneously, the regional meteorological station temperature data at each observation time are subtracted from a pre-set reference temperature to obtain the air temperature difference at each observation time. The air temperature difference values ​​at all observation times are arranged in chronological order to form an air temperature difference sequence. The reference temperature is the average air temperature at all observation times.

[0050] For example, Figure 2 This is a schematic diagram of the air temperature difference sequence. Figure 2 It demonstrates the objective fluctuation pattern of the air temperature difference in the monitoring area as a function of the observation time.

[0051] S2. Based on the time-series differential phase sequence and the air temperature difference sequence, calculate the relative phase temperature sensitivity coefficient between candidate scattering point pairs that have a physical adjacency relationship, and use the relative phase temperature sensitivity coefficient to perform thermodynamic stripping on the original phase difference between the candidate scattering point pairs to obtain the corrected phase difference that reflects the permanent deformation characteristics of the roadbed.

[0052] It should be noted that, due to the thermal expansion and contraction elastic physical properties of highway pavement materials and the microscopic differences in material distribution between adjacent candidate scattering points, the physical response of each measuring point to air temperature fluctuations is different. If the original temporal differential phase sequence is directly used for spatial path unwrapping, it is very easy to cause the elastic deformation error to be maliciously propagated in the spatial topology map. Therefore, this invention combines the air temperature difference sequence to extract the relative phase temperature sensitivity coefficient, performs thermodynamic stripping on the temporal differential phase sequence, and restores the corrected phase difference that reflects the stability of the geological structure.

[0053] Specifically, an initial spatial topology map is constructed based on the three-dimensional spatial coordinates of each candidate scattering point, and candidate scattering point pairs with physical adjacency are determined. To avoid calculation errors caused by spatial phase jumps, the absolute value of the time-series differential phase difference of the candidate scattering point pairs with physical adjacency at the same observation time must be less than pi.

[0054] Based on the temporal differential phase sequence and the air temperature difference sequence, the relative phase temperature sensitivity coefficient between candidate scattering point pairs with physical adjacency is calculated:

[0055] ;

[0056] In the formula, Indicates the first The candidate scattering point and the first The relative phase temperature sensitivity coefficient between candidate scattering points, the th The candidate scattering point and the first Each candidate scattering point has a physical adjacency relationship; Indicates the total number of observation times; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Represents the sequence of air temperature difference values. The air temperature difference at each observation time; It represents a very small positive number and is used to prevent the denominator from being zero in calculations.

[0057] In the formula, when the temporal differential phase difference between candidate scattering points With air temperature difference The more drastic the change, the higher the calculated relative phase temperature sensitivity coefficient. The larger the value, the lower the relative phase temperature sensitivity coefficient. The smaller the value, the better the intrinsic thermal response properties of the material hidden in the microwave interference phase, thus achieving the effect of removing periodic elastic expansion disturbances.

[0058] In the formula, the smallest positive number The experience range is usually to In this embodiment, the smallest positive number The value takes The reason for this value is that its magnitude is much smaller than the normal fluctuation variance of the sum of squares of the ambient temperature differences. While ensuring the non-zero characteristic of the denominator, it avoids interference with the accuracy of the true relative phase temperature sensitivity coefficient. In other embodiments, implementers can use extremely small positive numbers based on the floating-point precision of the calculation system. The settings.

[0059] Further, calculate the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time:

[0060] ;

[0061] In the formula, Indicates the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering point and the first The relative phase temperature sensitivity coefficient between candidate scattering points; Represents the sequence of air temperature difference values. The air temperature difference at each observation time. When the air temperature difference... relative phase temperature sensitivity coefficient The larger the product, the larger the bias component subtracted from the original temporal differential phase difference, thus removing the periodic contamination of the spatial gradient by the environmental thermal field, and improving the phase difference correction. It presents the characteristics of permanent roadbed settlement more purely.

[0062] For example, Figure 3This diagram illustrates the phase difference correction before and after correction for candidate scattering point pairs with physical adjacency, reflecting the effect of eliminating differences in road surface temperature sensitivity on improving radar interferometric phase quality.

[0063] S3. Based on the degree of discrete fluctuation of the corrected phase difference on the time axis, and combined with the spatial physical distance between the candidate scattering point pairs, a spatiotemporal stability index that takes into account both temporal immunity and spatial coherence is constructed.

[0064] It should be noted that, since road surfaces are prone to microscopic shear fatigue damage when subjected to high-frequency dynamic loads, the physical fatigue caused by traffic rolling manifests as highly unpredictable and violent oscillations in the temporal phase gradient. Therefore, this invention constructs a spatiotemporal stability index that takes into account both temporal immunity and spatial coherence by evaluating the degree of discrete fluctuation of the corrected phase difference on the time axis and combining the physical propagation attenuation properties between microwave scatterers. This index is used to identify and block fragile connected edges in the subsequent initial spatial topology graph.

[0065] Specifically, the spatial physical distance between candidate scattering point pairs with physical adjacency is calculated using the three-dimensional spatial coordinates of each candidate scattering point. Based on the corrected phase difference and the spatial physical distance, the spatiotemporal stability index between candidate scattering point pairs with physical adjacency is calculated.

[0066] ;

[0067] In the formula, Indicates the first The candidate scattering point and the first The spatiotemporal stability index between candidate scattering points, the th The candidate scattering point and the first Each candidate scattering point has a physical adjacency relationship; Indicates the first The candidate scattering point and the first The spatial physical distance between candidate scattering points; Indicates the distance to zero constant; This represents the maximum value extraction function, used to guarantee... Not less than the distance zero constant To avoid the risk of the denominator being zero; Represents an exponential function with the natural constant as its base; This represents the phase normalization constant; Indicates the total number of observation times; Indicates the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time; Indicates the first The candidate scattering point and the first The average corrected phase difference between the candidate scattering points at all observation times; Represents the absolute value symbol.

[0068] Among them, the phase difference correction Deviation from the average The absolute deviation essentially reflects the severity of nonlinear fatigue oscillations generated by the road surface under high-frequency heavy traffic loads, while the spatial physical distance... This reflects the natural attenuation of the geological correlation between microwave scatterers. The greater the absolute deviation and the greater the spatial physical distance... A larger value indicates a severe internal structural loosening or hidden cracks between adjacent candidate scattering point pairs, and weak spatial energy transfer correlation. In this case, the calculated spatiotemporal stability index... The smaller the value, the smoother the deformation transition and the more stable the geological structure in the region. This invention constructs a spatiotemporal stability index, transforming the actual physical damage characteristics and distance attenuation laws into penalty terms for network weights, thereby achieving the effect of identifying and blocking fragile connected edges while preserving high-confidence phase integral channels.

[0069] In the formula, the distance to zero constant is... The empirical range is typically 0.1 meters to 1 meter; in this embodiment, the distance is zero constant. The value of 0.5 meters is chosen because the physical limits of candidate scattering points along the highway in the highest resolution image are usually within this range. This avoids range collapse caused by multiple scatterers in the same pixel. In other embodiments, the implementer can determine the range protection zero constant based on the spatial resolution of the actual radar image. The settings.

[0070] Phase normalization constant This is used to convert purely numerical radians into a strictly dimensionless proportion, preventing extreme numerical divergence caused by directly applying the exponent to large phase errors. Its value should cover the limit error range of conventional interferometric phase unwrapping. Based on this mathematical and physical boundary mapping criterion, the phase normalization constant... The empirical range is typically from 1 radian to 3.14 radians. In this embodiment, the phase normalization constant is... Take pi In other embodiments, the implementer may set the phase normalization constant according to the numerical stability requirements of the system. The settings.

[0071] S4. Using the spatiotemporal stability index as the topological edge connection weight in the initial spatial topology graph, the optimal connection path for phase unwrapping is determined by executing the maximum spanning tree optimization algorithm, and spatial deduction is performed according to the connection order of the optimal connection path to obtain the true phase sequence after removing meteorological interference. The settlement rate of candidate scattering points is determined based on the true phase sequence, and a settlement warning is triggered based on the determination result of the settlement rate.

[0072] It should be noted that, since traditional unwrapping paths often blindly cross areas with hidden cracks in the road surface, causing local phase jumps that lead to global unwrapping collapse errors, this invention introduces the spatiotemporal stability index into the topology optimization algorithm to drive the integral path to actively avoid low-stability areas with severe fatigue and loose geological structures, separate the atmospheric delay component and extract the settlement rate to accurately trigger settlement early warning.

[0073] Specifically, the spatiotemporal stability index is used as the connection weight of the topological edges in the initial spatial topology graph. Redundant topological edges with a spatial physical distance greater than the connectivity distance threshold are removed from the initial spatial topology graph to generate the final spatial topology graph. The connectivity distance threshold is obtained by: obtaining the basic skeleton network of the initial spatial topology graph based on the minimum spanning tree algorithm; calculating the average spatial physical distance of all connected edges in the basic skeleton network and three times the standard deviation; and determining the connectivity distance threshold by the sum of the two values.

[0074] The maximum spanning tree algorithm is applied to the final spatial topology graph to find the optimal connection path with the maximum cumulative edge weight. For example, Figure 4 To obtain the final spatial topology map and optimal connection path diagram of a 500-meter section of a monitored area, this paper shows the spatial connectivity pattern of the candidate scattering point network in this local area after screening by connectivity distance threshold and weighting based on spatiotemporal stability index.

[0075] Furthermore, from all candidate scattering points, a candidate scattering point located on the periphery of the monitoring area with a stable geological structure is selected as a stable reference point. Since this stable reference point is in a geologically stable state and has no actual settlement or deformation, its absolute phase reference at each observation time is set to 0. Using this stable reference point as the starting root node for spatial integration derivation, the adjacent candidate scattering points on the optimal connection path are solved step by step according to the connectivity order of the optimal connection path. The specific derivation logic is as follows: along the connectivity direction of the optimal connection path, two adjacent candidate scattering points on the optimal connection path are respectively taken as the starting candidate scattering point and the target candidate scattering point, where the starting candidate scattering point in the first derivation is the stable reference point. First, the integer ambiguity of the temporal differential phase difference between the starting candidate scattering point and the target candidate scattering point is calculated to recover the true relative phase difference. Then, this relative phase difference is accumulated to the solved absolute phase value of the starting candidate scattering point, thereby calculating the absolute phase value of the target candidate scattering point. Since the spatial incremental transmission process is carried out synchronously and independently at all observation times, it is possible to reconstruct the original relative and entangled phase differences in the network into the absolute phase sequence of each candidate scattering point over time, thus obtaining the untangled phase sequence.

[0076] Since atmospheric delay phase exhibits random high-frequency abrupt changes in time characteristics and large-scale low-frequency smoothness in space characteristics, a high-pass filter is applied to the unwrapped phase sequence in the time dimension using a preset time high-pass filter window to separate the high-frequency time sequence containing atmospheric delay and random noise. Then, a low-pass filter is applied to this high-frequency time sequence in the spatial dimension using a preset spatial low-pass filter window to smooth out local random noise and purify the atmospheric delay component sequence. In a preferred embodiment of the invention, the low-pass filtering process uses a Gaussian filtering algorithm to effectively suppress spatial high-frequency random noise. Furthermore, the corresponding atmospheric delay component sequence is subtracted one by one from the unwrapped phase sequence to obtain the true phase sequence stripped of meteorological interference.

[0077] The empirical range of the time high-pass filter window is usually 200 to 500 days. In this embodiment, 365 days is used because this window can effectively filter out long-period deformation trends to extract short-period atmospheric fluctuations. The empirical range of the space low-pass filter window is usually 800 to 2000 meters. In this embodiment, 1200 meters is used because this window can smooth local high-frequency noise to obtain a large-scale atmospheric delay surface. In other embodiments, the implementer can set the time high-pass filter window and the space low-pass filter window according to the atmospheric flow characteristics of the monitoring area.

[0078] Furthermore, for each candidate scattering point along the highway, the true phase sequence of the corresponding candidate scattering point is converted into a line-of-sight deformation sequence using the inherent radar microwave wavelength of the radar satellite hardware. Then, combined with the line-of-sight geometric projection transformation relationship of the radar incident angle, the line-of-sight deformation sequence is converted into a vertical temporal settlement sequence. A linear regression algorithm is used to fit and calculate the temporal settlement sequence, yielding the settlement rate of each candidate scattering point along the highway.

[0079] To obtain the settlement rate warning threshold: Extract the time-series settlement sequence of the stable reference point, calculate the ratio of the settlement value difference between adjacent observation times in the time-series settlement sequence to the corresponding time interval, and obtain the instantaneous settlement rate of the stable reference point in each observation period. Extract the maximum absolute value of the instantaneous settlement rate in all observation periods as the extreme value of the background settlement rate of the stable reference point caused by the background disturbance of the natural environment. Since the absolute value has been extracted, the extreme value of the background settlement rate is a positive value representing the rate of disturbance error. The extreme value of the background settlement rate is numerically superimposed with the allowable settlement rate critical limit specified in the highway subgrade engineering structure specification. The result of the numerical superposition of the two is used as the settlement rate warning threshold, where the allowable settlement rate critical limit is a positive value representing the maximum allowable settlement rate.

[0080] Since surface subsidence manifests as negative displacement in vertical space, a positive subsidence rate at a candidate scattering point indicates upward uplift of the corresponding surface, while a negative rate indicates downward subsidence. Therefore, when the subsidence rate at a candidate scattering point is less than the negative value of the subsidence rate warning threshold, the absolute rate of downward subsidence at that point exceeds the safety limit. In this case, a subsidence warning command is generated for the three-dimensional spatial coordinates of the candidate scattering point.

[0081] For example, Figure 5 This diagram illustrates the distribution of settlement rates at candidate scattering points and the triggering of early warnings. When the columnar settlement rate of some candidate scattering points breaks through the negative value of the settlement rate early warning threshold, it indicates that the absolute downward settlement rate at the physical coordinates of the candidate scattering point has exceeded the allowable range for engineering safety. Based on this, the system generates and outputs targeted settlement early warning commands, thereby achieving quantitative identification of fatigue settlement risks on highway surfaces.

[0082] This invention also discloses a highway surface subsidence early warning system based on temporal InSAR, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement the highway surface subsidence early warning method based on temporal InSAR according to this invention.

[0083] The system also includes other components well known to those skilled in the art, such as communication buses and communication interfaces, the settings and functions of which are known in the art and will not be described in detail here.

[0084] In the description of this specification, "multiple" or "several" means at least two, such as two, three or more, unless otherwise expressly and specifically defined.

[0085] While this specification has shown and described numerous embodiments of the invention, it will be apparent to those skilled in the art that such embodiments are provided by way of example only. Many modifications, alterations, and alternatives will occur to those skilled in the art without departing from the spirit and essence of the invention. It should be understood that various alternatives to the embodiments of the invention described herein may be employed in the practice of this invention.

Claims

1. A method for early warning of surface subsidence along highways based on time-series InSAR, characterized in that, include: The remote sensing image sequence and synchronous meteorological data of the monitoring area are acquired, and the temporal differential phase sequence and air temperature difference sequence of candidate scattering points are generated through differential interferometry. Based on the time-series differential phase sequence and the air temperature difference sequence, the relative phase temperature sensitivity coefficient between candidate scattering point pairs with physical adjacency is calculated, and the original phase difference between the candidate scattering point pairs is thermodynamically stripped using the relative phase temperature sensitivity coefficient to obtain the corrected phase difference that reflects the permanent deformation characteristics of the roadbed. Based on the degree of discrete fluctuation of the corrected phase difference on the time axis, and combined with the spatial physical distance between the candidate scattering point pairs, a spatiotemporal stability index that takes into account both temporal immunity and spatial coherence is constructed. Using the spatiotemporal stability index as the topological edge connection weight in the initial spatial topology graph, the optimal connection path for phase unwrapping is determined by executing the maximum spanning tree optimization algorithm. Spatial deduction is then performed according to the connection order of the optimal connection path to obtain the true phase sequence after removing meteorological interference. The settlement rate of candidate scattering points is determined based on the true phase sequence, and a settlement warning is triggered based on the determination result of the settlement rate.

2. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, The generation of the temporal differential phase sequence of candidate scattering points and the air temperature difference sequence through differential interferometry includes: The remote sensing image sequence was registered and differentially interferometrically processed using an external digital elevation model to obtain a differential interferogram sequence. The amplitude deviation features of the pixels in the differential interferogram sequence were extracted, and each pixel was divided into multiple level categories using the natural breakpoint classification method. The pixel points in the level category with the smallest amplitude deviation feature were selected as candidate scattering points. The phase values ​​of the candidate scattering points at each observation time were extracted from the differential interferogram sequence and arranged in chronological order to obtain the temporal differential phase sequence of the corresponding candidate scattering points. The air temperature difference at each observation time is obtained by subtracting the regional meteorological station temperature data from the pre-set reference temperature. The air temperature difference at each observation time is then arranged in chronological order to form an air temperature difference sequence.

3. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, The relative phase temperature sensitivity coefficient satisfies the expression: ； In the formula, Indicates the first The candidate scattering point and the first The relative phase temperature sensitivity coefficient between candidate scattering points, the th The candidate scattering point and the first Each candidate scattering point has a physical adjacency relationship; Indicates the total number of observation times; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Represents the sequence of air temperature difference values. The air temperature difference at each observation time; It represents a very small positive number and is used to prevent the denominator from being zero in calculations.

4. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, The corrected phase difference satisfies the expression: ； In the formula, Indicates the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering points at the th in the The temporal differential phase at each observation time; Indicates the first The candidate scattering point and the first The relative phase temperature sensitivity coefficient between candidate scattering points; Represents the sequence of air temperature difference values. The air temperature difference at each observation time.

5. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, The spatiotemporal stability index satisfies the expression: ； In the formula, Indicates the first The candidate scattering point and the first The spatiotemporal stability index between candidate scattering points, the th The candidate scattering point and the first Each candidate scattering point has a physical adjacency relationship; Indicates the first The candidate scattering point and the first The spatial physical distance between candidate scattering points; Indicates the distance to zero constant; This represents the function for extracting the maximum value. Represents an exponential function with the natural constant as its base; This represents the phase normalization constant; Indicates the total number of observation times; Indicates the first The candidate scattering point and the first The candidate scattering points at the th in the Corrected phase difference at each observation time; Indicates the first The candidate scattering point and the first The average corrected phase difference between the candidate scattering points at all observation times; Represents the absolute value symbol.

6. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, Before executing the maximum spanning tree optimization algorithm, the following steps are also included: The basic skeleton network of the initial spatial topology graph is obtained based on the minimum spanning tree algorithm; the average and standard deviation of the spatial physical distance of all connected edges in the basic skeleton network are calculated; the sum of the average and three times the standard deviation is determined as the connected distance threshold; redundant topological edges with spatial physical distances greater than the connected distance threshold are removed from the initial spatial topology graph.

7. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, The spatial extrapolation, performed according to the connectivity order of the optimal connection path, yields the true phase sequence after removing meteorological interference, including: From all candidate scattering points, a candidate scattering point located on the periphery of the monitoring area with a stable geological structure is selected as a stable reference point, and its absolute phase reference is set to 0. Using the stable reference point as the starting root node, along the connection direction of the optimal connection path, two adjacent candidate scattering points on the optimal connection path are respectively taken as the starting candidate scattering point and the target candidate scattering point. Integer ambiguity is calculated for the temporal differential phase difference between the starting candidate scattering point and the target candidate scattering point to recover the relative phase difference. The relative phase difference is accumulated to the absolute phase value of the starting candidate scattering point to obtain the absolute phase value of the target candidate scattering point. The spatial extrapolation process is executed synchronously at all observation times to reconstruct the unwrapped phase sequence of each candidate scattering point. The unwound phase sequence is filtered in the time dimension using a time high-pass filter window to separate the high-frequency time sequence; the high-frequency time sequence is filtered in the spatial dimension using a spatial low-pass filter window to extract the atmospheric delay component sequence; the atmospheric delay component sequence is subtracted from the unwound phase sequence to obtain the true phase sequence.

8. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 1, characterized in that, The step of determining the sedimentation rate of candidate scattering points based on the true phase sequence includes: The true phase sequence is converted into a line-of-sight deformation sequence, and the line-of-sight deformation sequence is converted into a vertical time-series settlement sequence by using the geometric projection transformation relationship of the radar incident angle. The sedimentation rate of candidate scattering points was obtained by fitting the time-series sedimentation sequence using a linear regression algorithm.

9. The method for early warning of surface subsidence along highways based on time-series InSAR according to claim 7, characterized in that, Triggering a settlement warning based on the settlement rate determination result includes: The maximum absolute value of the instantaneous settlement rate of the stable reference point in all observation periods is calculated as the extreme value of the background settlement rate caused by the background disturbance of the natural environment at the stable reference point. The extreme value of the background settlement rate is numerically superimposed with the allowable settlement rate critical limit in the highway subgrade engineering structure specification to obtain the settlement rate warning threshold. A settlement warning command is generated in response to a candidate scattering point having a settlement rate less than a negative value of the settlement rate warning threshold.

10. A highway surface subsidence early warning system based on time-series InSAR, characterized in that, include: A processor and a memory, wherein the memory stores computer program instructions that, when executed by the processor, implement the method for early warning of road surface subsidence based on time-series InSAR according to any one of claims 1-9.

Images