A highway pavement uneven deformation identification method based on time-series remote sensing data

Abstract

The application discloses a highway pavement uneven deformation identification method based on time sequence remote sensing data and relates to the technical field of engineering measurement and deformation monitoring. The method comprises the following steps: S1, a virtual measurement control point sequence is arranged along a target highway design centerline at a preset measurement density to form a deformation measurement reference; multi-temporal synthetic aperture radar image data and synchronous environment data are acquired, and registration and resampling are performed according to the deformation measurement reference to generate a standardized time sequence observation data sequence; and S2, the time sequence observation data sequence is subjected to interference processing, and a dynamic stability index is set by combining the environment data to analyze the physical law of phase change. The method realizes spatial continuous, non-contact and periodic deformation measurement and evaluation of the whole route of the target highway by using the wide coverage advantage of satellite remote sensing and performing unified calculation based on the virtual measurement reference network arranged continuously along the highway line shape.

Description

Technical Field

[0001] This invention relates to the field of engineering measurement and deformation monitoring technology, specifically a method for identifying uneven deformation of highway pavement based on time-series remote sensing data. Background Technology

[0002] Uneven deformation of highway pavement is one of the main manifestations of highway distress, affecting driving comfort and safety, and in severe cases, leading to pavement structural damage and major traffic accidents. Traditional deformation measurement methods (such as leveling) are limited by efficiency and cost, making it difficult to achieve efficient and continuous deformation surveys along the entire highway. Time-series synthetic aperture radar interferometry provides a new data source for wide-area deformation monitoring, but its existing general methods have significant shortcomings when applied to highly continuous linear engineering projects like highways: they fail to integrate highway network topology and pavement mechanical continuity constraints, resulting in insufficient spatial continuity and engineering reliability of the final calculated deformation, making it difficult to use as accurate engineering measurement results to define and quantify uneven deformation segments of the pavement.

[0003] The invention patent with publication number CN120031392A discloses a method, system, and storage medium for road surface deformation monitoring and control. It collects road surface data by deploying a sensor network, uses an LSTM network to predict deformation risks and formulate graded strategies, and then automatically optimizes traffic control schemes through reinforcement learning algorithms to achieve intelligent monitoring and adaptive maintenance decisions for road surface conditions.

[0004] However, the above and similar technical solutions have the following shortcomings: they focus on real-time or short-term responses to deformations that have already occurred, making it difficult to effectively identify and preventively warn of slow-accumulating deformations; their monitoring capabilities are limited by preset fixed sensor locations and mobile inspection routes, making it difficult to achieve synchronous assessment of the overall deformation status of the road network in a wide-area, full-coverage, and spatially continuous manner. Summary of the Invention

[0005] The purpose of this invention is to provide a method for identifying uneven deformation of highway pavement based on time-series remote sensing data, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for identifying uneven deformation of highway pavement based on time-series remote sensing data, comprising:

[0007] S1. Along the design centerline of the target highway, a sequence of virtual measurement control points is set up according to a preset measurement density to form a deformation measurement reference; multi-temporal synthetic aperture radar image data and synchronous environmental data are acquired, and registration and resampling are performed according to the deformation measurement reference to generate a standardized time-series observation data sequence.

[0008] S2. Perform interferometric processing on the time-series observation data sequence, analyze the physical laws of phase change in conjunction with environmental data, and set dynamic stability index; based on the deformation measurement benchmark and highway network topology constraints, select road surface steady-state deformation measurement points from the time-series observation data sequence according to the dynamic stability index;

[0009] S3. Based on the topological constraints of the highway network and the preset pavement mechanical structure parameters, the temporal phase of the measurement points is constrained, unwrapped, and filtered, and the temporal deformation sequence of each measurement point along the imaging geometric line of sight is calculated.

[0010] S4. Based on the time-series deformation sequence of at least two different observation geometries, a unified measurement adjustment model is constructed with the three-dimensional displacement of each measurement point as the parameter to be estimated. Vertical dominant constraints and spatial continuous smoothness constraints are introduced, and the vertical time-series deformation sequence of each measurement point is adaptively solved by partitioned variance component estimation.

[0011] S5. Perform multi-scale sliding analysis and abrupt change detection on the vertical time-series deformation sequence along the highway design centerline, divide the uniform deformation segments and quantify the deformation differences between segments, and output the measurement results of non-uniform deformation according to the predefined risk assessment rules.

[0012] Furthermore, the step S1 of arranging a sequence of virtual measurement control points according to a preset measurement density specifically includes:

[0013] Acquire historical time-series deformation measurement data covering the target highway, extract historical deformation anomaly areas and deformation gradient areas, and extract historical disease association areas based on the highway's historical maintenance records;

[0014] The extracted feature regions are spatially fused to generate a comprehensive risk indication map along the highway alignment, and a benchmark deployment density model related to the risk level of the map is established.

[0015] Identify key sections of highway design and use them as optimization factors to optimize the density of corresponding sections in the benchmark layout density model;

[0016] Based on the optimized baseline layout density model, the sequence of virtual measurement control point locations is adaptively determined, and each virtual measurement control point is assigned an engineering attribute label; the engineering attribute label includes at least the road structure type and the design station number.

[0017] Furthermore, the method for setting the dynamic stability index in S2 specifically includes:

[0018] Differential interferometry is performed on the time-series observation sequence to generate the original time-series differential phase sequence for each spatial measurement point;

[0019] By combining synchronously acquired environmental data, a statistical response model of phase change and environmental factors is established; using this model, the periodic phase components caused by environmental factors are removed from the original time-series differential phase sequence to obtain the environmentally compensated time-series phase sequence.

[0020] Based on the time-series phase sequence after environmental compensation, a composite dynamic stability index is constructed that combines temporal stability with spatial consistency. The evaluation threshold of the dynamic stability index is set differently according to the road structure type of the measurement point.

[0021] Furthermore, the method for selecting the road surface steady-state deformation measurement points in S2 includes:

[0022] Measurement points with time-dimensional stability better than the first threshold are selected to form a candidate point set; within the candidate point set, isolated outliers that do not conform to the deformation behavior of their direct topological neighbors are removed based on the spatial dimension consistency index.

[0023] The candidate points after topology optimization are mapped to the nearest virtual measurement control points based on their coordinates in the three-dimensional spatial network nodes.

[0024] For each successfully mapped candidate point, the engineering attribute label corresponding to the virtual measurement control point it is mapped to is bound, and it is finally determined as the road surface steady-state deformation measurement point.

[0025] Furthermore, S3 specifically includes:

[0026] Based on the highway network topology, adjacent measurement points are connected, and each connection is assigned a composite weight of mechanical transmission constraint weight and observation quality weight to obtain a composite weight network;

[0027] Based on the composite weighted network, adaptive phase unwrapping of structural constraints is performed;

[0028] A road surface mechanical response-driven state model is used to perform physical constraint filtering on the unwrapped phase;

[0029] The filtered phase estimate is converted into a temporal deformation sequence along the imaging geometry line of sight, and the unwrapping reliability index is output simultaneously.

[0030] Furthermore, the adaptive phase unwrapping of the structural constraints includes: selecting the main deformation transmission path in the composite weight network based on the spatial distribution of the mechanical transmission constraint weights; when unwrapping along the main path is blocked due to poor local observation quality, dynamic path reorganization is triggered, and based on the similarity of the mechanical transmission constraint weights, a parallel alternative transmission path is automatically searched and switched in the topology network to maintain the continuity of unwrapping.

[0031] Furthermore, the method for constructing the road surface mechanical response driving state model includes:

[0032] Multiple sub-mechanical response models, each corresponding to a typical pavement structure type, are initialized in parallel to form a competitive model set;

[0033] During the filtering process, the matching likelihood of each sub-mechanical response model to the current time phase is calculated, and weights are dynamically allocated or a dominant model is selected for state estimation.

[0034] Furthermore, the three-dimensional displacement includes a vertical displacement component and a horizontal displacement component; the vertical dominant constraint is achieved by assigning a higher prior weight to the vertical displacement component in the parameter to be estimated than that to the horizontal displacement component; the spatial continuous smooth constraint applies a smoothing penalty term to the difference in the three-dimensional displacement of all adjacent point pairs based on the connection relationship between adjacent measurement points in the highway network topology.

[0035] Furthermore, the partitioning of the variance component estimation is based on the source track of the observation and the road structure type to which the measurement point belongs; in each iteration, the variance component of each partition is independently estimated based on the adjustment residual, and the weight of the corresponding partition observation in the measurement adjustment model is dynamically adjusted accordingly until the weights converge.

[0036] Furthermore, S5 specifically includes:

[0037] Multi-scale sliding analysis was performed on the vertical temporal deformation sequence along the highway design centerline to extract the deformation trend, fluctuation intensity and frequency domain characteristics in each window, and then fused to generate a multi-dimensional feature sequence.

[0038] For the multi-dimensional feature sequence, a probabilistic mutation detection method that does not rely on a fixed threshold is used to calculate the mutation confidence of the multi-dimensional feature sequence at each mileage point along the route. The continuous mileage section where the mutation confidence is consistently higher than the preset judgment threshold is identified as the deformation transition zone, and its center mileage position is extracted as the boundary point of uneven deformation.

[0039] Based on the boundary points, uniform deformation segments are divided, the dominant deformation modes of each segment are identified, and the differences in average rate, trend and fluctuation between adjacent segments are quantified to form a difference vector.

[0040] The difference vector and the engineering attribute labels bound to the measurement points constituting each uniform deformation segment are input into a predefined risk assessment rule base. Through rule matching, the comprehensive risk level, dominant cause inference and treatment priority of each segment are output, and measurement results including design station number and assessment conclusion are generated.

[0041] Compared with the prior art, the beneficial effects of the present invention are:

[0042] A method for identifying uneven deformation of highway pavement based on time-series remote sensing data is proposed. This method deeply integrates deformation measurement results with road design attributes and performs risk assessment based on a predefined rule base. It directly outputs an engineering report containing risk level, cause inference, and treatment priority, thereby enhancing the decision support value of monitoring results.

[0043] Meanwhile, by leveraging the wide-area coverage advantage of satellite remote sensing and performing unified calculations based on a virtual measurement benchmark network continuously deployed along the highway alignment, continuous, non-contact, and periodic deformation measurement and assessment of the entire target highway route was achieved. By processing historical archived remote sensing data, a long-term road surface deformation sequence was constructed, and multi-scale sliding analysis and probabilistic abrupt change detection were combined with multi-scale trend analysis and abrupt change detection to perform multi-scale sliding analysis and probabilistic abrupt change detection on the time-series deformation sequence, enabling early identification and warning of slowly accumulating uneven deformation. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of the road surface uneven deformation identification method of the present invention;

[0045] Figure 2 This is a schematic diagram illustrating the effect of the virtual measurement control point sequence layout according to the present invention;

[0046] Figure 3 This is a schematic diagram of the adaptive phase unwrapping method for structural constraints according to the present invention;

[0047] Figure 4 This is a schematic diagram of the orbit observation fusion and adaptive adjustment method for solving vertical deformation according to the present invention;

[0048] Figure 5 This is a schematic diagram of the intelligent division and evaluation method for road deformation sections according to the present invention. Detailed Implementation

[0049] 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 embodiments of the present invention, and not all embodiments. 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.

[0050] like Figure 1 As shown, the present invention provides a technical solution: a method for identifying uneven deformation of highway pavement based on time-series remote sensing data, comprising:

[0051] S1. Along the design centerline of the target highway, a sequence of virtual measurement control points is set up according to a preset measurement density to form a deformation measurement reference; multi-temporal synthetic aperture radar image data and synchronous environmental data are acquired, and registration and resampling are performed according to the deformation measurement reference to generate a standardized time-series observation data sequence.

[0052] The specific steps of setting up a sequence of virtual measurement control points according to a preset measurement density include:

[0053] Historical time-series deformation measurement data covering the target highway is obtained, historical deformation anomaly areas and deformation high gradient areas are extracted, and historical disease association areas are extracted based on the highway's historical maintenance records.

[0054] It is important to note that, using the historical time-series deformation monitoring data, the standard deviation of the deformation rate time series for each measurement point is calculated. Points with a standard deviation greater than 1.5 times the average standard deviation of all points along the entire road section are spatially clustered and marked as deformation anomaly areas. The annual average deformation rate difference between adjacent monitoring points (within 100 meters) along the route direction is calculated, and road sections with an absolute difference greater than 5 mm / year are marked as high deformation gradient areas. The station intervals corresponding to typical defects described in the maintenance records such as "uneven settlement," "wavy bulging," and "severe rutting" are directly mapped onto the target highway electronic map and marked as defect-related areas.

[0055] The extracted feature regions are spatially fused to generate a comprehensive risk indication map along the highway alignment, and a benchmark deployment density model related to the risk level of the map is established.

[0056] It is important to note that the vector data of the three types of feature areas mentioned above are spatially overlaid on the same highway alignment map to generate a comprehensive risk indication map. In this map, each road segment is assigned a risk level (e.g., low, high, medium). The initial density rules of the preset baseline layout density model are as follows: for high-risk road segments, the preset measurement density is one virtual measurement control point every 30m; for medium-risk road segments, the density is one point every 60m; and for low-risk road segments, the density is one point every 100m.

[0057] Key sections of highway design are identified and used as optimization factors to optimize the density of corresponding sections in the benchmark layout density model.

[0058] It is important to note that key sections such as bridge approach transition sections, culvert top overburden sections, and high embankment sections are identified from the target highway design documents, and their chainage ranges are listed as optimization factors. These key sections are spatially correlated with the benchmark layout density model. When a road section in the model falls within a key section, its layout density is doubled based on the density value corresponding to the original risk level.

[0059] Based on the optimized baseline layout density model, the sequence of virtual measurement control point locations is adaptively determined, and each virtual measurement control point is assigned an engineering attribute label.

[0060] It is important to note that, starting from the highway's initial chainage, the plane coordinates (obtained by interpolation from the design alignment coordinate data) and elevation (obtained by interpolation from the design longitudinal profile) of each virtual measurement control point are calculated and determined sequentially along the highway design centerline, strictly following the intervals specified by the optimized benchmark layout density model, forming an ordered position sequence. An attribute record is created for each point in the sequence, and the engineering attribute label must at least include the road structure type (obtained from the design documents, such as flexible base asphalt pavement) and the design chainage (corresponding highway mileage). Figure 2 This is a sequence of virtual measurement control points (top view) adaptively deployed along the design centerline of the highway. Hollow small circles represent points deployed at a lower density (e.g., 100m) in general areas, while solid large circles represent points deployed at a higher density (e.g., 30m) in critical sections. The black center solid line represents the highway design direction and station direction. "k100+000", "k100+200", etc. are design station numbers. The shaded area represents the risk area / critical section.

[0061] Centered on the aforementioned virtual measurement control points, extend a lane width (e.g., 3.75m) to the left and right along the cross-section of the highway on the horizontal plane to generate a regular grid with a resolution of 1m×1m. Merge all grids generated from the control points to form a continuous three-dimensional spatial grid covering the target highway surface and adjacent areas, where each node has geodetic coordinates. For each SAR image, using its precise orbital data and radar parameters (center incident angle, heading angle), and based on the range-Doppler geolocation model, batch calculate the precise pixel coordinates (row and column numbers, usually floating-point numbers) of each node in the three-dimensional spatial grid corresponding to that image. Based on the calculated floating-point pixel coordinates, use a bilinear interpolation algorithm to extract the complex observation values ​​(including amplitude and phase information) corresponding to each node in the three-dimensional spatial grid from the complex matrix of the SAR image. The observation values ​​of all nodes for a single image constitute an observation set. Repeat the above steps to generate the corresponding observation sets for all real-phase SAR images with all orbits (e.g., ascending, descending, etc.). A structured database is constructed using the 3D spatial grid node ID as the first dimension, time as the second dimension, and track identifier as the third dimension to store all complex observation values, resulting in a standardized four-dimensional observation data volume. Based on the highway design cross-section diagram, all nodes whose horizontal positions fall within the design pavement boundary line are selected from the 3D spatial network and defined as valid pavement observation points. The standardized time-series observation data sequence is the sequence of complex observation values, sorted by timestamp, extracted from the standardized four-dimensional observation data volume and corresponding to each valid pavement observation point.

[0062] S2. The time-series observation data sequence is subjected to interferometric processing, and the physical laws of phase change are analyzed in combination with environmental data to set dynamic stability index. Based on the deformation measurement benchmark and highway network topology constraints, the road surface steady-state deformation measurement points are selected from the time-series observation data sequence according to the dynamic stability index.

[0063] The method for setting the dynamic stability index specifically includes:

[0064] Differential interferometry is performed on the time-series observation sequence to generate the original time-series differential phase sequence for each spatial measurement point.

[0065] It is important to note that the small baseline set (SBAS) strategy is used to generate differential interferograms from the standardized continuous observation data sequence generated in step S1. Specifically, a temporal baseline threshold of ≤48 days and a spatial vertical baseline threshold of ≤150m are set, and interferometry is performed on all pairable image combinations. For each pair of images that satisfy the baseline (master image and slave image), the differential interferometric phase of each pair of effective road observation points is calculated. :

[0066]

[0067] in, The phase extracted from complex observations, To obtain the image phase, Main image phase, For orbital error phase, This refers to the topographic phase. The differential interferometric phases of all interferometric pairs are arranged in chronological order (based on the time of the main image) to form the original temporal differential phase sequence for each effective observation point on the road surface. The length of this sequence is equal to the number of interferometric pairs, and each value represents the relative phase change within a time interval.

[0068] By combining synchronously acquired environmental data, a statistical response model of phase change and environmental factors is established. Using this model, the periodic phase components caused by environmental factors are removed from the original time-series differential phase sequence to obtain the environmentally compensated time-series phase sequence.

[0069] It is important to note that the environmental data acquired covering the target area and at the observation time should have a temporal resolution synchronized with the SAR image acquisition time, and a spatial resolution of 0.25° × 0.25°. The environmental data should include at least: surface temperature, total columnar water vapor content, and surface air pressure. For each effective observation point on the road surface, a multiple linear regression model should be established between its original temporal differential phase sequence and the synchronized environmental data, for example:

[0070]

[0071] Where i is the index of the measurement point; These represent the changes in environmental factors between the time points when the master and slave images were acquired; a, b, and c are regression coefficients, which can be solved using the least squares method. The residuals contain deformation signals and noise. Using the regression model, the phase components predicted by changes in environmental factors are subtracted from the original time-series differential phase sequence at each point to obtain the residual sequence. This refers to the time-phase sequence after environmental compensation.

[0072] Based on the environmentally compensated temporal phase sequence, a composite dynamic stability index is constructed that combines temporal dimension stability with spatial dimension consistency.

[0073] It is important to note that the time series standard deviation of the environmentally compensated time-series phase sequence at each point should be calculated. And calculate the standard deviation of the sequence on the first difference. Then the time dimension stability index Defined as:

[0074]

[0075] The larger the value, the more stable the phase at that point in time. In the three-dimensional spatial network, for each effective observation point on the road surface, find its two closest neighboring points (one on the left and one on the right) along the road cross-section, forming its direct topological neighborhood. Calculate the average Pearson correlation coefficient between this point and the environmentally compensated temporal phase sequences of its two neighboring points. And calculate the root mean square value of the difference between the phase sequence at this point and the phase sequences of two neighboring points. Then the spatial dimension consistency index Defined as:

[0076]

[0077] The larger this value, the more consistent the deformation behavior of this point is with that of its neighboring points. Ultimately, this results in a composite dynamic stability index. Defined as:

[0078]

[0079] in, These are weighting coefficients, which can be preset based on experience. This emphasizes the dominant role of time stability. The evaluation threshold for the composite dynamic stability index is differentiated based on the virtual measurement control point associated with the point and the road structure type in the engineering attribute label: for bridge approach pavements or cement concrete pavements, which are more sensitive to environmental factors such as temperature and prone to abrupt changes, a stricter threshold is set. ≥1.5; For conventional asphalt concrete pavements, a general threshold is set. ≥1.0; For pavement sections of embankment subgrades, a wider threshold is set because they may have more complex deformation patterns. ≥0.8.

[0080] The method for selecting the road surface steady-state deformation measurement points includes:

[0081] Measurement points with time-dimensional stability better than the first threshold are selected to form a candidate point set; within the candidate point set, isolated outliers that do not conform to the deformation behavior of their direct topological neighbors are removed based on the spatial dimension consistency index.

[0082] It should be noted that, for all valid road surface observation points, their calculations are performed. .Will The value is greater than the value corresponding to its road structure type. In the threshold Points with a reference value (e.g., 0.6 for conventional asphalt pavement) are initially selected as the candidate point set, which is the first threshold screening. Within the candidate point set, their spatial dimension consistency index is checked. If the index value of a point is lower than the empirical threshold corresponding to its road structure type (e.g., the threshold for conventional asphalt pavement is set to 0.4), and its index is significantly lower than (e.g., lower than 50%) the average value of other candidate points in the direct topological neighborhood, then the point is determined to be an isolated outlier and is removed from the candidate point set.

[0083] The candidate points, after topology optimization, are mapped to the nearest virtual measurement control points based on their coordinates in the three-dimensional spatial network nodes.

[0084] It is important to note that the candidate points after topology optimization are themselves nodes in a three-dimensional spatial network. Based on their three-dimensional coordinates, the Euclidean distance between them and all the virtual measurement control points set up in step S1 is calculated, and they are mapped to the nearest virtual measurement control point. If the minimum distance is greater than a preset tolerance (e.g., 10m), the mapping is considered to have failed, and the candidate point is discarded.

[0085] For each successfully mapped candidate point, the engineering attribute label corresponding to the virtual measurement control point it is mapped to is bound, and it is finally determined as the road surface steady-state deformation measurement point.

[0086] S3. Based on the topological constraints of the highway network and the preset pavement mechanical structure parameters, the temporal phase of the measurement points is constrained, unwrapped, and filtered, and the temporal deformation sequence of each measurement point along the imaging geometric line of sight is calculated.

[0087] Based on the highway network topology, adjacent measurement points are connected, and each connection is assigned a composite weight of mechanical transmission constraint weight and observation quality weight to obtain a composite weight network.

[0088] It is important to note that all pavement steady-state deformation measurement points are linearly connected along the route direction according to their actual station number sequence on the highway design drawing, forming a highway network topology. Each measurement point is only bidirectionally connected to its preceding and following points (with a spacing equal to the optimized layout density in S1). The preset pavement mechanical structure parameters for each measurement point are obtained from the highway design documents, including at least: base modulus (MPa), surface layer thickness (cm), and subgrade type (e.g., cohesive soil, sandy soil). For two adjacent measurement points i and j, if their road structure types are the same and the difference in design parameters (base modulus, surface layer thickness) is less than 10%, they are considered to have strong mechanical conductivity and are assigned a higher foundation weight. If the road structure types are different, or the key design parameters differ by more than 30%, a mechanical discontinuity interface is considered to exist, and a lower isolation weight is assigned. For cases where the parameter difference is between 10% and 30%, linear interpolation is used, and the result is calculated based on the percentage of parameter difference. The value ranges from 0.1 to 1.0. Then the mechanical transmission constraint weight of link ij is... Defined as:

[0089]

[0090] Using the calculations for each measurement point in the preceding steps To measure the phase observation quality. For connection ij, define its observation quality weight. This is the average of the quality indicators at both ends, after normalization:

[0091]

[0092] in, Let i be the composite dynamic stability index value. Let j be the composite dynamic stability index value. This represents the maximum value of the composite dynamic stability index across all measurement points. Observation quality weight. The value range is (0,1). For each topological connection ij, its final composite weight is... for:

[0093]

[0094] All measurement points and their pathways Weighted connections constitute a composite weighted network, which is a weighted undirected graph.

[0095] Based on the aforementioned composite weighted network, adaptive phase unwrapping of structural constraints is performed, such as... Figure 3 As shown.

[0096] It is important to note that when choosing a composite weight network... The continuous connection path with the highest value is selected as the main path for deformation propagation. An anchor point on the main path, where the known deformation is negligible, is chosen, such as at the starting point of the route. Points with extremely high values ​​or points known to be stable through field measurements, and their unwrapping phase is assumed. The value is 0. Starting from the anchor point, along the main path in station order, for each subsequent point j and its preceding point i on the path, perform the following operations (performed independently for each interference time t):

[0097] First, obtain the entanglement phase difference between the two points at time t. :

[0098]

[0099] Where t is the time-series index. For the time interval corresponding to the t-th interference pair, the differential interferometric phase at point j after environmental compensation is given. Let be the differential interference phase of point i after environmental compensation within the time interval corresponding to the t-th interference pair.

[0100] Secondly, calculate the possible untangling phase difference between the two points. :

[0101]

[0102] Where k is an integer (usually -1, 0, or 1). Based on the physical constraints (smooth and continuous deformation at adjacent points), an optimal value of k is selected. Specifically, a smoothness cost function is calculated:

[0103]

[0104] in, It is the composite weight of ij. It is a predicted value extrapolated from the phase difference of historically untangled points. The k value that minimizes the smoothness cost is selected.

[0105] Finally, the untangling phase of point j at time t for:

[0106]

[0107] Based on the above operations, the phase of all points on the main path at time t is unwound sequentially.

[0108] When integrating along the main path to a connection (m,n), if its If the entanglement phase difference of the connection is extremely low (e.g., <0.2), or if the entanglement phase difference is abnormally high for most times t (e.g., >π), making it impossible to reliably determine the value of k, then the untangling is considered blocked. In the composite weighted network, starting from the blocked point m and ending at the blocked point n, K (e.g., K = 3) parallel alternative propagation paths are searched, avoiding the blocked connection (m, n). Parallelism means that the alternative path is spatially parallel to the main path through the mn segment, and the intermediate nodes on the path should be roughly located within the station interval of m and n, without large station backtracking or jumps. The search algorithm is Dijkstra's algorithm, and the connection cost is defined as 1 / W. c To prioritize W c For each alternative path, starting from the existing unwrapped phase value at point m, integrate to point n using the one-dimensional path integral method to obtain R possible candidate unwrapped phase values ​​for point n. Calculate the standard deviation of these candidate values. If the standard deviation is less than a preset threshold (e.g., π / 2), take the average value as the unwrapped phase value for point n and update the network connection state. If the standard deviation is too large, temporarily shelve point n until more points are unwrapped before processing it.

[0109] A road surface mechanical response-driven state model is used to perform physical constraint filtering on the unwrapped phase.

[0110] It is important to note that, based on common highway structures, multiple sub-mechanical response models corresponding to typical pavement structure types are initialized in parallel as a competing model set. For example, Model A (elastic layered system model) is applicable to most asphalt or cement concrete pavements, and its state equation describes the deformation increment. It is linearly related to the current cumulative deformation d(t-1) and the ambient temperature T(t):

[0111]

[0112] Among them, a A b is the springback / creep coefficient. A W is the temperature sensitivity coefficient. A(t) For process noise, This represents the effective temperature change within the pavement structural layers (usually the surface layer or base layer). Model B (viscoelastic creep model) is suitable for pavements with high asphalt content or soft soil subgrade sections, and its equation of state includes a memory term:

[0113]

[0114] Among them, a B c is the instantaneous response coefficient.B Let d(tf) be the creep memory coefficient, and w be the deformation displacement d at the f-th historical moment before the current time t. f W is the weight (which can be set based on f). B(t) This represents process noise. Model C (rigid body displacement model) is suitable for bridge expansion joints or severely cracked areas, assuming that the deformation is mainly discontinuous jumps, and its state equation is:

[0115]

[0116] W C(t) This represents process noise, and abrupt changes are allowed when its variance is large. The state vector for each of the above models is [d(t), v(t)]. T , where v(t) is the deformation rate, using an extended Kalman filter framework.

[0117] The input to the model is For each time step t, each model independently performs prediction and update steps. After the update at step t, for model M (M = A, B, C), the difference between its observed predicted value and the actual value is calculated. covariance matrix S M(t) Then the matching likelihood L of the model at the current time step is... M(t) The calculation is as follows:

[0118]

[0119] Weights of model M at time t :

[0120]

[0121] Finally, the filtered state estimate (displacement and velocity) of the measurement point at time t is a weighted sum of the state estimates from the three models:

[0122]

[0123] Among them, X M(t) Let be the optimal estimation vector of the state of the current measurement point obtained by the Mth sub-mechanical response model after independently performing state estimation through its own extended Kalman filter at time t. If the weights of a certain model... If the value is greater than 0.7 for more than 5 consecutive time steps, it is determined to be the dominant model. In this case, only this model can be used for filtering to improve the rate until its weight continues to decrease.

[0124] The filtered phase estimate is converted into a temporal deformation sequence along the imaging geometry line of sight, and the unwrapping reliability index is output simultaneously.

[0125] It is important to note that arranging the displacement components in the filtered state vector of each measurement point in chronological order yields the temporal deformation sequence of that point along the imaging geometric line of sight, in mm. The unwrapping reliability indices include: a path reliability index, which is the proportion of dynamic path reorganization occurring during unwrapping of that point out of the total number of integration steps; a lower proportion indicates higher reliability. A model consistency index is calculated by determining the weight quotient H of the competing model sets during the filtering process.

[0126]

[0127] The lower the entropy value, the more dominant a certain model is, the clearer the mechanical behavior, and the more reliable the filtering results. A weighted average of these two indices is used to generate a detangling reliability index between 0 and 1 for each measurement point; the higher the value, the more reliable the deformation sequence at that point.

[0128] S4. Based on the time-series deformation sequence of at least two different observation geometries, a unified measurement adjustment model is constructed with the three-dimensional displacement of each measurement point as the parameter to be estimated. Vertical dominant constraints and spatial continuous smoothness constraints are introduced, and the vertical time-series deformation sequence of each measurement point is adaptively solved by partitioned variance component estimation.

[0129] like Figure 4 As shown, this invention provides a method for solving vertical deformation through multi-track observation fusion and adaptive adjustment.

[0130] It is important to note that the input data consists of time-series deformation sequences from at least two different satellite orbits (e.g., orbit A is an ascending orbit, and orbit B is a descending orbit) output from step S3. For each measurement point i, each orbit g, and each time t, a deformation d along the line-of-sight (LOS) direction of that orbit is obtained. LOS(i,g,t) The unit is mm. For each measurement point i and each time t, the three-dimensional displacement to be estimated is defined as:

[0131] X i(t) =(U i(t) V i(t) W i(t) ) T

[0132] Among them, U i(t) V represents the east-west displacement (east is positive), in mm. i(t) The displacement is in the north-south direction (north is positive), and the unit is mm; W i(t) The vertical displacement component (positive upwards) is represented in mm.

[0133] For each observation d LOS(i,g,t) Establish its three-dimensional displacement X at time t. i(t) Geometric relationships:

[0134] d LOS(i,g,t) =(s e(o) ,s n(o) ,s u(o) )×X i(t) +v(i,g,t)

[0135] Among them, (s e(o) ,s n(o) ,s u(o) ) represents the projected components of the unit vector of orbit g in the LOS direction along the east-west, north-south, and vertical directions; v(i,g,t) represents the observation noise. The observation equations for all measurement points, all orbits, and all times are combined and written in matrix form:

[0136] L = B × X + V

[0137] Where L represents all d LOS The column vector consists of the observation values, X is the column vector consisting of all the three-dimensional displacement parameters to be estimated, B is the design matrix (composed of LOS projection vectors), and V is the observation noise vector.

[0138] The introduction of vertical dominant constraints involves attaching a virtual prior observation equation to each parameter to be estimated, constraining it to be near 0. Specifically: for the vertical displacement component W at all points... i(t) Add the prior constraint equation: 0=W i(t) +w w And assign a high prior weight P W (e.g., P) W =1.0, corresponding to a priori standard deviation of 1.0 mm). For the horizontal displacement component U at all points i(t) and V i(t) Similarly, add the prior constraint equation: 0=U i(t) +w U 0=V i(t) +w V However, it is given a lower prior weight P. h (e.g., P) h =0.1, corresponding to a prior standard deviation of approximately 3.16 mm. Where, w w w u w v This represents the random error term in the prior constraint equation. By setting P... W >P h This allows for the allocation of higher prior weights to vertical displacement components than to horizontal displacement components during adjustment, guiding the solution results to better align with the actual physical laws of road surfaces, where vertical deformation is the dominant feature.

[0139] The introduced spatially continuous smoothing constraint, based on the highway network topology, establishes three spatial smoothing constraint equations for each pair of adjacent measurement points (i,j) at each time t:

[0140] 0=U i(t) -U j(t) +w u_s

[0141] 0=V i(t) -V j(t) +w v_s

[0142] 0=W i(t) -W j(t) +w w_s

[0143] Among them, w u_s w v_s w w_s Let P be the random error term of the spatial smoothing constraint equations. Assign a smoothing penalty weight P to these smoothing constraint equations. s This weight value is set based on the expected deformation gradient between adjacent points. For example, assuming that the three-dimensional displacement difference between adjacent points (30m apart) within a day does not exceed 0.5 mm, then P can be set. s =1 / 0.5 2 =4.0, which constitutes a smoothing penalty term for the difference in displacement between adjacent points.

[0144] The adaptive solution of the vertical temporal deformation sequence of each measurement point through partitioned variance component estimation includes: all d LOS The observations are partitioned based on the following criteria: the source track (e.g., ascending track observation set and descending track observation set); and the road structure type of the measurement point (e.g., asphalt pavement observation set, concrete pavement observation set, and bridge approach section observation set). These two criteria are then combined to form multiple partitions, such as ascending track-asphalt pavement partitions and descending track-bridge approach section partitions. Initial weights are assigned to the observations in each partition, typically set to be equal. The observation equations, vertical dominant constraint equations, and spatial smoothing constraint equations are combined to form a unified measurement adjustment model. The parameters are solved using weighted least squares, and the residual vector of this adjustment is calculated. For the z-th observation partition, the variance components of that partition are independently estimated using the Helmholtz variance component estimation formula based on its residuals. Based on the estimated variance components, the weights of the observations in each partition are dynamically adjusted in the next iteration; the new weights are inversely proportional to the variance components.

[0145] Check the weights of all partition observations relative to the changes in the previous iteration. If the largest relative change is less than a preset threshold (e.g., 1%), the weights are considered to have converged, and the iteration stops. Otherwise, the iteration continues. After convergence, the vertical displacement component of each measurement point i at each time t, calculated by the final adjustment, is the desired vertical temporal deformation sequence.

[0146] S5. Perform multi-scale sliding analysis and abrupt change detection on the vertical time-series deformation sequence along the highway design centerline, divide the uniform deformation segments and quantify the deformation differences between segments, and output the measurement results of non-uniform deformation according to the predefined risk assessment rules.

[0147] like Figure 5 As shown, the present invention provides a method for intelligent division and evaluation of road deformation segments;

[0148] Specifically:

[0149] Multi-scale sliding analysis was performed on the vertical temporal deformation sequence along the highway design centerline to extract the deformation trend, fluctuation intensity and frequency domain characteristics within each window, and then fused to generate a multi-dimensional feature sequence.

[0150] It is important to note that the vertical temporal deformation sequence W at each measurement point should be... i(t) The points are arranged in ascending order according to their associated design station numbers, forming a continuous sequence along the route. Multiple windows of different lengths are set to capture deformation patterns in different spatial ranges, such as short windows (100m), medium windows (300m), and long windows (800m). The sliding step size is uniformly set to 20m (approximately the average point spacing). For each window position, calculation is only performed if the window contains at least three measurement points. Feature extraction is performed independently for each window and each window length: deformation trend, and W values ​​for all measurement points within the window. i(t) The sequence is calculated, and its average deformation rate (mm / year) over the most recent year is taken as the median rate of all points within the window as the deformation trend of that window. Fluctuation intensity is calculated for W at each measurement point within the window. i(t) For the sequence, calculate the standard deviation of the remaining sequence after detrending (i.e., subtracting the linear fit value), and take the average of the standard deviations of all points as the fluctuation intensity of the window, in mm. For frequency domain characteristics, calculate the W value for all measurement points within the window. i(t)The sequence is averaged to obtain the average deformation sequence for that window. A Fast Fourier Transform is then performed on this average sequence to extract its dominant frequency component (i.e., the frequency with the highest energy, measured in units of 1 / year). This value reflects the periodicity of the deformation; if there is no obvious periodicity, it is recorded as 0. For each sliding step position along the route (i.e., every 20m), there are three feature values ​​at three window scales, resulting in a total of nine feature values. Arranging these nine feature values ​​in a fixed order forms a nine-dimensional feature vector for that mileage position. Arranging the feature vectors of all mileage positions in order of their station numbers constitutes a multi-dimensional feature sequence along the entire target highway.

[0151] For the multi-dimensional feature sequence, a probabilistic mutation detection method that does not rely on a fixed threshold is used to calculate the mutation confidence of the multi-dimensional feature sequence at each mileage point along the route. The continuous mileage segment where the mutation confidence is consistently higher than the preset judgment threshold is identified as the deformation transition zone, and its center mileage position is extracted as the boundary point of uneven deformation.

[0152] It is important to note that this step is based on the Bayesian online change point detection algorithm. It assumes that the feature vectors of adjacent uniformly deformed segments follow the same multidimensional normal distribution; a significant change in the distribution parameters is considered a sudden change.

[0153] Starting from the route's origin, the algorithm processes data sequentially along a multi-dimensional feature sequence. For the currently processed mileage point e, the algorithm maintains a running length r (i.e., how many mileage points have been traversed since the last hypothetical change point). Given all the data conditions, the posterior probability of a change point occurring at the current point c (i.e., r returning to 0) is P(change point|data). This posterior probability is calculated recursively by comparing the likelihoods of two hypotheses: the current data point belongs to the current segment; the current data point belongs to a new segment. In the calculation, the distribution of the feature vectors is described by the estimated values ​​of the mean and covariance matrix of the current segment's data. The calculated posterior probability is used as the change confidence of the mileage point, with a value between 0 and 1. The change confidence sequence of the entire route is scanned, and all continuous mileage segments with a change confidence greater than a preset threshold (e.g., 0.7) are selected.

[0154] Based on the boundary points, uniform deformation segments are divided, the dominant deformation modes of each segment are identified, and the differences in average rate, trend and fluctuation between adjacent segments are quantified to form a difference vector.

[0155] It is important to note that each uniform deformation segment does not contain boundary points, and its deformation characteristics are assumed to be relatively uniform. For each uniform deformation segment, analyze the W values ​​of all measurement points it contains. i(t)The rules for determining the dominant deformation mode include: if the absolute value of the average deformation rate of more than 80% of the points within a segment is <2 mm / year and the average value of the fluctuation intensity is <3 mm, it is determined to be a stable segment; if the average deformation rate of more than 80% of the points has the same sign and the absolute value of the rate is >5 mm / year, while the fluctuation intensity is relatively small, it is determined to be a uniform subsidence / uplift segment; if the average value of the fluctuation intensity is >5 mm, regardless of the trend, it is determined to be a fluctuating deformation segment; if the eigenvalues ​​of the frequency domain features within the segment show a significant and consistent non-zero period, it is determined to be a periodic deformation segment.

[0156] For adjacent uniformly deformed segments, a difference vector is calculated. Specifically, for uniformly deformed segment 1 and uniformly deformed segment 2, a 5-dimensional difference vector is calculated, including: △Average rate = |Average rate of all points in uniformly deformed segment 1 - Average rate of all points in uniformly deformed segment 2|, in mm / year; △Trend type, which is 1 if the dominant modes of the two segments are different, otherwise 0; △Fluctuation intensity = |Average fluctuation intensity value of uniformly deformed segment 1 - Average fluctuation intensity value of uniformly deformed segment 2|, in mm; △Frequency domain characteristics = |Average frequency domain characteristics of uniformly deformed segment 1 - Average frequency domain characteristics of uniformly deformed segment 2|, in 1 / year; Boundary strength = Spatial gradient of deformation rate within 50m of the boundary point, in (mm / year) / m.

[0157] The difference vector and the engineering attribute labels bound to the measurement points constituting each uniform deformation segment are input into a predefined risk assessment rule base. Through rule matching, the comprehensive risk level, dominant cause inference and treatment priority of each segment are output, and measurement results including design station number and assessment conclusion are generated.

[0158] It should be noted that the predefined risk assessment rule base includes the following example rules: Example Rule 1: If the average rate (Δ) is greater than 10 mm / year and the boundary gradient is greater than 0.2 mm / year / m, and the boundary point is located on the bridge approach road surface, then the comprehensive risk level of uniformly deformed road segment 1 is high risk, the dominant cause is the potential area of ​​bridge approach slab settlement, and the treatment priority is P1 (highest); if the trend type (Δ) is 1 and the fluctuation intensity (Δ) is greater than 4 mm, and one of the two segments is a soft soil subgrade, then the comprehensive risk level of uniformly deformed road segment 1 is medium risk, the dominant cause is the difference in subgrade soil quality, and the treatment priority is P2; if the average rate (Δ) is less than 3 mm / year and the fluctuation intensity (Δ) is less than 2 mm, then the comprehensive risk level of uniformly deformed road segment 1 is low risk, the dominant cause is normal deformation, and the treatment priority is P4.

[0159] The difference vectors between each uniformly deformed segment and its preceding and following adjacent segments, along with the segment's own engineering attribute labels, are input into a rule base for matching. If multiple rules are triggered, the conclusion with the highest overall risk level is adopted. Finally, a record is generated for each uniformly deformed segment, including: design station interval, overall risk level, inferred dominant cause, and treatment priority. All evaluation records for uniformly deformed segments are integrated to form a structured report on the measurement results of uneven deformation of highway pavement.

[0160] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended embodiments and their equivalents.

Claims

1. A method for identifying uneven deformation of highway pavement based on time-series remote sensing data, characterized in that, include: S1. Along the design centerline of the target highway, a sequence of virtual measurement control points is set up according to the preset measurement density to form a deformation measurement benchmark; Acquire multi-temporal synthetic aperture radar image data and synchronous environmental data, and perform registration and resampling based on the deformation measurement benchmark to generate a standardized time-series observation data sequence; S2. Perform interferometric processing on the time-series observation data sequence, analyze the physical laws of phase change in conjunction with environmental data, and set dynamic stability index; based on the deformation measurement benchmark and highway network topology constraints, select road surface steady-state deformation measurement points from the time-series observation data sequence according to the dynamic stability index; S3. Based on the topological constraints of the highway network and the preset pavement mechanical structure parameters, the temporal phase of the measurement points is constrained, unwrapped, and filtered, and the temporal deformation sequence of each measurement point along the imaging geometric line of sight is calculated. S4. Based on the time-series deformation sequence of at least two different observation geometries, a unified measurement adjustment model is constructed with the three-dimensional displacement of each measurement point as the parameter to be estimated. Vertical dominant constraints and spatial continuous smoothness constraints are introduced, and the vertical time-series deformation sequence of each measurement point is adaptively solved by partitioned variance component estimation. S5. Perform multi-scale sliding analysis and abrupt change detection on the vertical time-series deformation sequence along the highway design centerline, divide the uniform deformation segment and quantify the deformation difference between segments, and output the measurement results of non-uniform deformation according to the predefined risk assessment rules. S3 specifically includes: Based on the highway network topology, adjacent measurement points are connected, and each connection is assigned a composite weight of mechanical transmission constraint weight and observation quality weight to obtain a composite weight network; Based on the composite weighted network, adaptive phase unwrapping of structural constraints is performed; A road surface mechanical response-driven state model is used to perform physical constraint filtering on the unwrapped phase; The filtered phase estimate is converted into a temporal deformation sequence along the imaging geometry line of sight, and the unwrapping reliability index is output simultaneously. The adaptive phase unwrapping of the structural constraints includes: selecting the main deformation transmission path in the composite weight network based on the spatial distribution of the mechanical transmission constraint weights; when unwrapping along the main path is blocked due to poor local observation quality, dynamic path reorganization is triggered, and based on the similarity of the mechanical transmission constraint weights, the network is automatically searched for and switched to a parallel alternative transmission path to maintain the continuity of unwrapping. The construction methods of the road surface mechanical response driving state model include: Multiple sub-mechanical response models, each corresponding to a typical pavement structure type, are initialized in parallel to form a competitive model set; During the filtering process, the matching likelihood of each sub-mechanical response model to the current time phase is calculated, and weights are dynamically allocated or a dominant model is selected for state estimation.

2. The method for identifying uneven deformation of highway pavement based on time-series remote sensing data according to claim 1, characterized in that: The sequence of virtual measurement control points arranged according to a preset measurement density in S1 specifically includes: Acquire historical time-series deformation measurement data covering the target highway, extract historical deformation anomaly areas and deformation gradient areas, and extract historical disease association areas based on the highway's historical maintenance records; The extracted feature regions are spatially fused to generate a comprehensive risk indication map along the highway alignment, and a benchmark deployment density model related to the risk level of the map is established. Identify key sections of highway design and use them as optimization factors to optimize the density of corresponding sections in the benchmark layout density model; Based on the optimized baseline layout density model, the sequence of virtual measurement control point locations is adaptively determined, and each virtual measurement control point is assigned an engineering attribute label; the engineering attribute label includes at least the road structure type and the design station number.

3. The method for identifying uneven deformation of highway pavement based on time-series remote sensing data according to claim 1, characterized in that: The method for setting the dynamic stability index in S2 specifically includes: Differential interferometry is performed on the time-series observation data sequence to generate the original time-series differential phase sequence for each spatial measurement point; By combining synchronously acquired environmental data, a statistical response model of phase change and environmental factors is established; using this model, the periodic phase components caused by environmental factors are removed from the original time-series differential phase sequence to obtain the environmentally compensated time-series phase sequence. Based on the time-series phase sequence after environmental compensation, a composite dynamic stability index is constructed that combines temporal stability with spatial consistency. The evaluation threshold of the dynamic stability index is set differently according to the road structure type of the measurement point.

4. The method for identifying uneven deformation of highway pavement based on time-series remote sensing data according to claim 1, characterized in that: The method for selecting the road surface steady-state deformation measurement points in S2 includes: Measurement points with time-dimensional stability better than the first threshold are selected to form a candidate point set; within the candidate point set, isolated outliers that do not conform to the deformation behavior of their direct topological neighbors are removed based on the spatial dimension consistency index. The candidate points after topology optimization are mapped to the nearest virtual measurement control points based on their coordinates in the three-dimensional spatial network nodes. For each successfully mapped candidate point, the engineering attribute label corresponding to the virtual measurement control point it is mapped to is bound, and it is finally determined as the road surface steady-state deformation measurement point.

5. The method for identifying uneven deformation of highway pavement based on time-series remote sensing data according to claim 1, characterized in that: The three-dimensional displacement includes a vertical displacement component and a horizontal displacement component; the vertical dominant constraint is achieved by assigning a higher prior weight to the vertical displacement component in the parameter to be estimated than that to the horizontal displacement component; the spatial continuous smooth constraint is based on the connection relationship between adjacent measurement points in the highway network topology, and applies a smoothing penalty term to the difference in the three-dimensional displacement of all adjacent point pairs.

6. The method for identifying uneven deformation of highway pavement based on time-series remote sensing data according to claim 1, characterized in that: The partitioning basis for the estimated partitioned variance components is the source track of the observations and the road structure type to which the measurement points belong; In each iteration, the variance components of each partition are independently estimated based on the adjustment residuals, and the weights of the corresponding partition observations in the measurement adjustment model are dynamically adjusted accordingly until the weights converge.

7. The method for identifying uneven deformation of highway pavement based on time-series remote sensing data according to claim 1, characterized in that: S5 specifically includes: Multi-scale sliding analysis was performed on the vertical temporal deformation sequence along the highway design centerline to extract the deformation trend, fluctuation intensity and frequency domain characteristics in each window, and then fused to generate a multi-dimensional feature sequence. For the multi-dimensional feature sequence, a probabilistic mutation detection method that does not rely on a fixed threshold is used to calculate the mutation confidence of the multi-dimensional feature sequence at each mileage point along the route. The continuous mileage section where the mutation confidence is consistently higher than the preset judgment threshold is identified as the deformation transition zone, and its center mileage position is extracted as the boundary point of uneven deformation. Based on the boundary points, uniform deformation segments are divided, the dominant deformation modes of each segment are identified, and the differences in average rate, trend and fluctuation between adjacent segments are quantified to form a difference vector. The difference vector and the engineering attribute labels bound to the measurement points constituting each uniform deformation segment are input into a predefined risk assessment rule base. Through rule matching, the comprehensive risk level, dominant cause inference and treatment priority of each segment are output, and measurement results including design station number and assessment conclusion are generated.

Images