Urban road collapse monitoring and early warning method and system based on laser point cloud data

Abstract

The application relates to a kind of city road surface collapse monitoring and early warning method and system based on laser point cloud data, solve the problem that the existing method neglects spatial correlation, prediction lacks physical constraint, and the high false alarm rate caused by rough warning grading, its method includes: obtaining road surface point cloud data by street lamp pole laser scanning equipment;Form time series data after compensation and denoising;Unified coordinate system based on ICP algorithm and calculate monitoring point displacement, use spatial neighborhood relationship to carry out consistency check and eliminate abnormal data;Extract displacement change rate and cumulative displacement amount to input LSTM to predict displacement trend;Combined with displacement rate dynamic adjustment threshold, generate graded warning signal and send to emergency processing center.The application has the following effects: through spatial consistency check and physical constraint prediction, realize the graded warning of high precision and low false alarm, improve the reliability of collapse risk identification.

Description

Technical Field

[0001] This invention relates to the field of urban infrastructure safety monitoring technology, and in particular to a method and system for monitoring and early warning of urban road collapse based on laser point cloud data. Background Technology

[0002] Urban road collapses are characterized by their suddenness and destructiveness. With the aging of underground pipelines and increased disturbance from construction, they have become a significant risk to public safety. Establishing an efficient road monitoring and early warning system to identify potential collapse risks at an early stage is a critical requirement for the safe operation and maintenance of urban infrastructure.

[0003] Currently, road surface deformation monitoring mainly employs fixed-point sensors such as settlement meters or 3D laser scanners for periodic observation. Sensors can monitor the displacement of key points over long periods, while point cloud technology can acquire the regional settlement distribution. After data filtering, noise reduction, and coordinate registration, a threshold method is typically used for early warning; that is, an alarm is triggered when the cumulative displacement or rate of change exceeds an empirical threshold.

[0004] However, existing technologies generally use isolated monitoring points as the basic analysis unit, lacking a systematic consideration of the spatial correlation between monitoring points. As a continuous medium, road surface deformation has inherent spatial synergy, with displacements in adjacent areas influencing each other. Isolated point analysis struggles to distinguish between genuine anomalies and measurement noise, and fails to capture regional coordinated deformation trends, resulting in a high false alarm rate and a tendency to miss early warning opportunities. Summary of the Invention

[0005] In order to achieve high-precision, low-false-alarm graded early warning through spatial consistency verification and physical constraint prediction, and improve the reliability of collapse risk identification, this application provides a method and system for monitoring and early warning of urban road collapse based on laser point cloud data.

[0006] Firstly, this application provides a method for monitoring and early warning of urban road collapse based on laser point cloud data, employing the following technical solution:

[0007] A method for monitoring and early warning of urban road collapses based on laser point cloud data includes:

[0008] Three-dimensional point cloud data is obtained by scanning the road surface with laser point cloud scanning equipment deployed on street light poles according to a preset sampling cycle;

[0009] Environmental compensation and denoising processing are performed on the 3D point cloud data. The processed data is stored according to time index to form a time series data set containing observations of a preset historical period.

[0010] Based on 3D point cloud data from different time periods in the time series dataset, the coordinate system is unified using the iterative nearest point algorithm. Based on the unified 3D point cloud data, each monitoring point on the road surface is determined, and the relative displacement and its horizontal and settlement components of each monitoring point are calculated. The relative displacement is checked for consistency based on the spatial neighborhood relationship of each monitoring point, and the relative displacement that does not meet the preset consistency standard is removed.

[0011] The displacement change rate of each monitoring point is extracted based on the last preset number of continuous sampling periods in the preset historical period, and the cumulative displacement of each monitoring point is extracted based on all sampling periods; the displacement change rate and cumulative displacement are input into a long short-term memory network to predict the displacement evolution trend for a preset time period in the future.

[0012] Based on the preset basic warning threshold, the threshold is dynamically adjusted in combination with the displacement change rate; based on the adjusted threshold, relative displacement and displacement evolution trend, a graded warning signal is generated and sent to the emergency response center and the urban disaster center.

[0013] By adopting the above technical solutions, abnormal data is eliminated through spatial consistency verification, and displacement trends are predicted by LSTM and thresholds are dynamically adjusted to achieve hierarchical early warning, effectively improving the accuracy and reliability of collapse risk identification.

[0014] Secondly, this application provides an urban road collapse monitoring and early warning system based on laser point cloud data, which adopts the following technical solution:

[0015] A city road collapse monitoring and early warning system based on laser point cloud data includes a memory, a processor, and a program stored in the memory and executable on the processor. When the program is loaded and executed by the processor, it implements the city road collapse monitoring and early warning method based on laser point cloud data as described in the first aspect. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating an urban road collapse monitoring and early warning method based on laser point cloud data, according to an embodiment of this application. Detailed Implementation

[0017] The present application will be further described in detail below with reference to the accompanying drawings.

[0018] Reference Figure 1 This application discloses a method for monitoring and early warning of urban road collapse based on laser point cloud data, comprising:

[0019] Step S1: Using a laser point cloud scanning device installed on the street light pole, the road surface is scanned according to a preset sampling period to obtain three-dimensional point cloud data.

[0020] Among them, laser point cloud scanning equipment refers to an active optical measurement device deployed on streetlight poles. It is a solid-state or mechanical rotating radar system that acquires three-dimensional coordinate information of the target surface by emitting pulsed laser beams and receiving echo signals. The preset sampling period refers to the time interval parameter set according to the road traffic load characteristics and soil deformation rate; under normal monitoring mode, it is 1-6 hours, shortened to 10-30 minutes during special periods (flood season, construction period). Three-dimensional point cloud data refers to a digital set of discrete points on the road surface containing their three-dimensional spatial coordinates, reflection intensity, and echo count, stored in PCD or LAS format, with a single frame data size of 10-500MB.

[0021] The necessary procedures are as follows:

[0022] Equipment Deployment and Field of View Configuration: The laser point cloud scanning equipment is fixed to the streetlight pole at a height of 3.5-6 meters above the ground using adjustable angle brackets. The downward angle is set to 15°-35°, ensuring a vertical field of view ≥30° and a horizontal field of view ≥120° covering the target road surface. Deploy the equipment in an array at intervals of 50-200 meters along the road's extension direction, with adjacent equipment overlapping by 10%-20%.

[0023] Periodic scanning and ranging: Scanning is triggered according to a preset sampling period. Time-of-flight ranging or phase ranging is used to scan the road surface with an angular resolution of 0.1°-0.4°. The ranging accuracy is controlled within ±2cm, and the ranging range is 5-150 meters. GPS / BeiDou timestamps are recorded simultaneously to ensure that the time alignment accuracy of multi-temporal data is better than 1 millisecond.

[0024] Point cloud acquisition and preprocessing: The three-dimensional coordinates, reflection intensity value, and echo count of each laser footprint are recorded to generate a raw point cloud containing spatial geometric and physical properties. Preliminary compression is performed using a built-in edge computing unit, and the data is transmitted to the cloud via 4G / 5G or wired network to provide data input for subsequent environmental compensation and noise reduction processing.

[0025] Step S2: Perform environmental compensation and denoising processing on the 3D point cloud data, and store the processed data according to the time index to form a time series data set containing observations of a preset historical period.

[0026] For information on environmental compensation and denoising of 3D point cloud data, please refer to steps S21 to S23.

[0027] The remaining steps are as follows:

[0028] Historical cycle parameter configuration: Based on the characteristics of road settlement evolution and algorithm requirements, set the duration of the historical cycle (e.g., the most recent 7 days) and the number of sampling cycles N to be retained (e.g., the most recent 168 cycles, corresponding to sampling per hour) to form a fixed-capacity circular buffer structure; configure the cycle overlap rate parameter to ensure that adjacent cycle data have temporal continuity.

[0029] Time-series data writing and version control: The single-frame point cloud data processed in steps S21-S23 is used as the latest observation value and written to the storage with timestamp, period number and data quality mark; the time-series chain is formed by append writing mode, each observation value is assigned a unique period ID, and a triple record of "period ID-timestamp-point cloud data" is established to ensure the immutability and traceability of time-series data.

[0030] Sliding window aggregation and set construction: Start a background maintenance thread to manage the observation queue according to the first-in-first-out (FIFO) principle; when a new observation is written, the earliest observation exceeding the preset historical period length is automatically removed to maintain the set containing exactly the data of the most recent N consecutive sampling periods; a doubly linked list index is built for all observations in the set to support ascending order (from early to late) and descending order (from late to early) traversal to form a complete time series observation sequence.

[0031] Multi-dimensional correlation index construction: In the time series dimension, a B+ tree index based on sampling timestamps is established to support range queries of "start time - end time"; in the period dimension, a hash index based on period sequence number is established to support fast truncation of "the most recent K periods"; in the data correlation dimension, a mapping table of monitoring point IDs in point cloud data of each period is maintained to ensure that observations of the same spatial location in different periods can be accurately correlated, providing cross-period data alignment capability for coordinate system one and displacement calculation in step S3.

[0032] Data integrity verification and supplementation: Periodically check the periodic continuity within the time series set, identify missing periods caused by equipment failure or network interruption; trigger the supplementation mechanism to restore the observations of the missing periods from edge device cache or backup storage, ensuring that there are no gaps in the data within the preset historical period, forming a continuous, complete set of observations that can be used for time series analysis.

[0033] Step S3: Based on the three-dimensional point cloud data of different time periods in the time series data set, coordinate system one is established using the iterative nearest point algorithm, and each monitoring point on the road surface is determined based on the unified three-dimensional point cloud data. The relative displacement and its horizontal and settlement components of each monitoring point are calculated. The relative displacement is checked for consistency based on the spatial neighborhood relationship of each monitoring point, and the relative displacement that does not meet the preset consistency standard is removed.

[0034] The iterative nearest-point algorithm refers to a registration method that unifies the coordinate system of multi-period point cloud data by minimizing the sum of squared distances between corresponding points in adjacent sampling periods and solving for the optimal rigid transformation matrix (rotation and translation parameters). Monitoring points are discrete control points extracted from the unified 3D point cloud data and placed at road surface feature locations (such as road marking corners, manhole cover edges, crack endpoints, or intersections of regular grids), used to quantitatively characterize local geometric changes in the road surface. Relative displacement refers to the positional offset vector of the same monitoring point in 3D space relative to the reference period (such as the first historical period or the previous period), which can be decomposed into horizontal displacement components (east-west and north-south offsets) and settlement components (vertical offsets). Spatial neighborhood relationships refer to the topological connections built based on the Euclidean distance between monitoring points, usually represented by a graph structure, where nodes are monitoring points and edges connect pairs of monitoring points whose spatial distance is less than a preset neighborhood radius. Consistency verification refers to the process of using graph convolutional networks to aggregate neighborhood displacement features, calculating the residual between the monitoring point's own displacement and the aggregated neighborhood displacement, and identifying and eliminating abnormal displacement values ​​caused by registration errors or local interference.

[0035] The necessary procedures are as follows:

[0036] Time-series data extraction and multi-station registration: Extract the 3D point cloud data of the current sampling period from the time-series data set, as well as the point cloud data of the first historical period (or deformation stabilization period) as the spatial reference; In view of the fixed viewing angle limitation caused by the layout of street light poles, when the single-station scanning range cannot cover the entire road surface, extract the synchronous point cloud of adjacent pole equipment, and perform multi-station stitching through the ICP algorithm to unify the segmented scanning data into the global reference coordinate system.

[0037] Iterative Nearest Point Coordinate System 1: Construct a kd-tree spatial index for the current period point cloud, and use the ICP algorithm for iterative optimization with the reference period point cloud as the target; set a distance threshold of 5-10cm to establish the relationship between corresponding points, and remove outliers caused by moving objects such as vehicles and pedestrians; solve the rotation matrix and translation vector through SVD decomposition, and terminate the iteration when the root mean square error is less than 1cm to achieve high-precision spatial alignment of the two periods of data.

[0038] Monitoring point deployment and cross-period tracking: In the unified benchmark point cloud road area, monitoring points are deployed using a regular grid method, such as every 1.0m along the longitudinal direction of the road (driving direction) and along the center line of each lane in the transverse direction, and feature monitoring points are added at manhole covers and expansion joints; For each unified point cloud, a neighborhood point set with a radius of 15-20cm is searched with the benchmark coordinates of the monitoring point as the center, and RANSAC plane fitting is used to determine the local road surface elevation, and the accurate three-dimensional coordinates of the monitoring point in each period are obtained by perpendicular projection, and a cross-period coordinate sequence is established.

[0039] Displacement calculation and component decomposition: Calculate the three-dimensional coordinate difference between the current period and the reference period for each monitoring point to obtain the relative displacement vector; establish a local coordinate system for the road (longitudinal along the road design centerline, horizontal perpendicular to the road, and vertical along the gravity direction), and decompose the displacement vector projection into: horizontal displacement component (including longitudinal and lateral offsets, used to identify pavement misalignment and lateral expansion) and settlement component (vertical elevation change, positive values ​​indicate uplift, and negative values ​​indicate subsidence); set a minimum detectable threshold of 3mm for the horizontal component, and distinguish between micro-deformation and significant deformation for the settlement component according to soil deformation characteristics, generating a displacement parameter set for each monitoring point.

[0040] The process of eliminating relative displacements that do not meet the preset consistency criteria can be referred to steps S31 to S33.

[0041] Step S4: Extract the displacement change rate of each monitoring point based on the last preset number of continuous sampling periods in the preset historical period, and extract the cumulative displacement of each monitoring point based on all sampling periods; input the displacement change rate and cumulative displacement into the long short-term memory network to predict the displacement evolution trend for a preset time period in the future.

[0042] The displacement change rate refers to the displacement change per unit time at a monitoring point within a few recent consecutive sampling periods (e.g., the most recent 12 periods, corresponding to the most recent 12 hours). It is calculated by dividing the difference between the current displacement and the historical displacement by the time interval and is used to characterize the recent development speed of pavement deformation. The cumulative displacement refers to the total cumulative displacement of the monitoring point from the baseline period (the monitoring start date) to the current moment, reflecting the long-term cumulative effect and overall stability of pavement deformation. Long Short-Term Memory (LSTM) networks are a special type of recurrent neural network structure that effectively captures long-term dependencies and short-term fluctuations in displacement time-series data by introducing forget gates, input gates, and output gates. It is suitable for processing pavement settlement time series with nonlinear and non-stationary characteristics. The displacement evolution trend refers to the predicted curve of the displacement of each monitoring point over a predetermined period (e.g., the next 24 hours, 72 hours, or 7 days), including the predicted displacement value and its confidence interval.

[0043] The necessary procedures can be referred to in steps S41 to S44.

[0044] Step S5: Based on the preset basic early warning threshold, the threshold is dynamically adjusted in combination with the displacement change rate; based on the adjusted threshold, relative displacement and displacement evolution trend, a graded early warning signal is generated and sent to the emergency response center and the urban disaster center.

[0045] The necessary procedures can be referred to in steps S51 to S55.

[0046] Before scanning the road surface according to the preset sampling period, it also includes:

[0047] Step S11: Control the laser point cloud scanning device to scan the calibration plate with preset size and reflectivity to obtain the point cloud data of the calibration plate.

[0048] The calibration board refers to a high-reflectivity reference object with a regular geometric shape (such as a planar rectangle or checkerboard pattern) deployed within the monitored road surface area. It provides a spatial reference with known geometric parameters for the laser point cloud scanning equipment. The preset size and reflectivity refer to the physical parameters of the calibration board pre-set according to the laser radar wavelength characteristics and monitoring accuracy requirements. These typically include planar dimensions (e.g., 1.0m × 1.0m or 0.5m × 0.5m), thickness (5-10mm), and surface reflectivity (≥80%, preferably 90%-95%) to ensure stable echo intensity in the laser band. The calibration board point cloud data refers to the discrete set of points containing three-dimensional spatial coordinates and reflection intensity information obtained after the laser point cloud scanning equipment scans the surface of the calibration board. This data is used for subsequent geometric feature extraction and system error analysis.

[0049] The necessary procedures are as follows:

[0050] Calibration plate setup and positioning: The calibration plate is horizontally positioned within the road surface area of ​​the streetlight pole's scanning field of view (usually located 10-30m directly in front of the scanning equipment, avoiding main traffic arteries). An aluminum alloy frame or ceramic substrate is used as the carrier, with a high-reflectivity microprism reflective film pasted on the surface or a diffuse reflection ceramic coating sprayed on it to ensure that the reflectivity is stable at 85%-95% at 905nm or 1550nm laser wavelength, forming a significant contrast with the asphalt pavement (reflectivity 10%-20%). The calibration plate is precisely positioned using a level and a total station to ensure that the plate surface levelness error is less than 0.5°, and the theoretical three-dimensional coordinates of the corner points of the calibration plate are recorded as the reference true value.

[0051] Scanning Equipment Control and Data Acquisition: Control the laser point cloud scanning equipment to perform a specific scan on the calibration board; adjust the pitch and horizontal angles of the scanning equipment to position the calibration board in the center of the scanning field of view (avoiding edge distortion), and set a high angular resolution (e.g., 0.1°-0.2°) to obtain dense point clouds; trigger the equipment to perform repeated scans on the calibration board surface for 3-5 seconds in static scanning mode (the equipment pauses continuous rotation and scans in a fixed posture) or slow scanning mode (reducing the scanning frame rate to 1-2Hz) to obtain multiple frames of raw point cloud data; simultaneously record the ambient temperature, humidity, and internal temperature of the equipment during the scan as metadata for subsequent data filtering.

[0052] Calibration board point cloud data acquisition: Extract point cloud segments with reflection intensity higher than a preset threshold (e.g., the echo intensity of the calibration board surface is more than 3 times higher than the road background) from the raw data stream of the scanning device; Superimpose and average multiple frames of scanning data, remove instantaneous noise points caused by vehicles passing by or dust, and generate calibration board point cloud data with high signal-to-noise ratio; Bind and store the acquired calibration board point cloud data with the corresponding timestamp, device number and environmental parameters to form a raw dataset containing geometric morphology information of the calibration board surface, providing data input for extracting geometric features and calculating system errors in step S12.

[0053] Step S12: Extract the geometric features of the calibration board point cloud data, and calibrate the incident angle error of the laser point cloud scanning device based on the deviation between the geometric features and the theoretical geometric parameters of the calibration board.

[0054] The geometric features of the calibration board point cloud data refer to the set of parameters extracted from the point cloud on the calibration board surface, characterizing its spatial shape and size. These include plane fitting parameters (normal vector, plane equation coefficients), corner point 3D coordinates, edge direction vectors, and plane dimensions (length, width, diagonal length). The theoretical geometric parameters of the calibration board refer to the ideal geometric parameters determined according to the calibration board design drawings, serving as a spatial reference. These include nominal length, nominal width, nominal thickness, and the theoretical spatial coordinates of the corner points. The incident angle error refers to the angular deviation between the actual emission direction and the ideal direction of the laser beam caused by installation attitude deviations (pitch angle, roll angle) or internal optomechanical structural distortions of the lidar. This manifests as a trapezoidal or parallelogram distortion of the rectangular calibration board in the scanned point cloud. Calibration refers to the process of calculating the incident angle error correction based on the deviation between the measured geometric features and theoretical parameters, and eliminating systematic errors through software compensation or physical adjustment.

[0055] The necessary procedures are as follows:

[0056] Geometric feature extraction: For the calibration board point cloud data obtained in step S11, the RANSAC algorithm is used to fit the plane and extract the normal vector, and edge noise is removed; four side lines are extracted in the fitted plane, and the intersection of adjacent side lines is calculated to obtain the three-dimensional coordinates of the four corner points; the Euclidean distance between the corner points is calculated to obtain the measured side length, width and diagonal length.

[0057] Deviation analysis and error calculation: Compare the measured side length with the theoretical nominal size and calculate the proportional distortion rate. When the distortion rate exceeds 1%, a significant incident angle error is identified. Analyze the shape of the quadrilateral formed by the corner points. If it exhibits a trapezoidal distortion that is wider at the top and narrower at the bottom, a pitch angle error is identified. The pitch angle correction is calculated using geometric relationships. Δh is half the difference between the upper and lower side lengths, and d is the scanning distance; similarly, the roll angle error is calculated based on the difference between the left and right side lengths, and the yaw angle error is calculated based on the perpendicularity deviation of the adjacent side.

[0058] Calibration Implementation and Verification: The calculated triaxial angle error is written into the external compensation matrix of the lidar. In subsequent data processing, the original scanning coordinates are rigidly transformed and corrected; or the equipment mounting bracket on the street light pole is physically fine-tuned to eliminate attitude deviation; the calibration plate is re-scanned for verification to ensure that the side length error is less than 5mm and the corner position deviation is less than 1cm. After meeting the road surface monitoring accuracy requirements, the road surface periodic scanning in step S1 is performed.

[0059] The road surface is scanned according to a preset sampling period, including:

[0060] Step S1A involves real-time acquisition of ambient temperature and humidity during the scanning process. Ambient temperature and humidity refer to atmospheric thermodynamic parameters along the laser beam propagation path. Temperature affects air density, and humidity affects refractive index; both together determine the actual laser propagation speed. Real-time acquisition refers to an environmental parameter measurement mechanism triggered synchronously with the laser scanning operation, ensuring that the environmental data used for compensation is consistent with the atmospheric conditions at the actual propagation time of the laser pulse.

[0061] The necessary procedures are as follows:

[0062] Sensor deployment and integration: An integrated temperature and humidity sensor module is deployed inside the main body cavity of the laser point cloud scanning equipment or on the outside of the protective cover adjacent to the laser emission window; the sensor probe is exposed to the outside atmosphere but protected by a waterproof and breathable membrane, and is connected to the radar main control unit through an I2C or RS485 interface to ensure that the environmental conditions along the actual propagation path of the laser beam are collected.

[0063] Synchronous Triggering and Acquisition: A triggering mechanism that is hard synchronized with the laser scanning cycle is adopted. Temperature and humidity sampling is triggered at the start of each scan (such as when each frame of a mechanical rotating radar passes through a 0° azimuth angle). The acquisition frequency is not lower than the scanning frame rate to ensure that the recorded temperature (accuracy ±0.5°C, range -20°C to +50°C) and humidity (accuracy ±3%RH, range 10%-95%RH) correspond to the environmental conditions of that scanning cycle.

[0064] Data encapsulation and quality labeling: The collected environmental parameters are bound to the corresponding 3D point cloud data frames and written into the metadata header of the data file, including a GPS / BeiDou synchronized UTC timestamp (millisecond level), temperature value (degrees Celsius, two decimal places), and humidity value (%RH, one decimal place); anomaly detection is implemented, and when the reading exceeds the physical range or the rate of change is abnormal, the data is marked as invalid and the previous period's valid value or standard atmospheric parameters (20°C, 50%RH) are used to avoid erroneous compensation.

[0065] Step S1B: Calculate the environmentally compensated speed of light value based on temperature, humidity, and a preset speed of light correction formula.

[0066] The preset light speed correction formula refers to a mathematical model, pre-established based on atmospheric optical properties, that describes the quantitative relationship between laser propagation speed and ambient temperature and humidity. It typically employs a simplified engineering approach to avoid complex calculations and meet the computational constraints of real-time compensation. The environmentally compensated light speed value is a physical quantity reflecting the actual laser propagation speed under the current temperature and humidity conditions, calculated by substituting the standard light speed or vacuum light speed into the aforementioned correction formula. It is used for the subsequent precise conversion of laser time-of-flight (ToF) to distance values.

[0067] The necessary procedures are as follows:

[0068] Environmental parameter input and calculation: The real-time ambient temperature T (unit: degrees Celsius, ℃) and relative humidity H (unit: percentage, %) collected in step S1A are substituted into the preset light speed correction formula to calculate the environmentally compensated light speed value c. env : .

[0069] Where c0 is the reference speed of light (taken as the vacuum speed of light, 299,792,458 m / s, or the standard atmospheric speed calibrated by the equipment), T is the current ambient temperature, and H is the current ambient relative humidity. This formula quantifies the comprehensive influence of environmental parameters on the air refractive index through linear coefficients of 0.0037 (temperature influence coefficient) and 0.0005 (humidity influence coefficient), thus achieving rapid correction of the speed of light.

[0070] Calculation Example and Compensation Necessity: For example, when the ambient temperature T=30 ∘ C. When the relative humidity H = 80%, the denominator of the calculation is 1 + 0.0037 × 30 + 0.0005 × 80 = 1.151, then c env The speed of light is approximately 260,463,000 m / s, which is significantly different from the speed of light under standard conditions. If a fixed standard speed of light is used for distance measurement, a systematic error will be introduced. Therefore, this step is necessary to calculate the speed of light after environmental compensation in real time to provide accurate parameters for distance compensation in step S1C.

[0071] Engineering optimization and lookup table method: To reduce edge computing load, pre-build coverage -20 ∘ C to +50 ∘ C. Light speed compensation lookup table (LUT) in the range of 10% to 95% RH, with 1 ∘ C and 5%RH are the step-size discrete storage corresponding to c. env The actual operation involves bilinear interpolation based on the collected temperature and humidity data to quickly obtain the compensated speed of light value, ensuring that the calculation delay is lower than the laser ranging cycle.

[0072] Step S1C: Real-time compensation of the laser ranging value based on the light speed value after environmental compensation.

[0073] The laser ranging value refers to the original distance measurement value between the laser point cloud scanning device and the target point on the road surface, calculated by the laser pulse round-trip flight time based on a preset reference light speed (usually the speed of light in a vacuum or standard atmospheric light). Real-time compensation refers to the process of synchronously and dynamically correcting the original ranging value based on the environmental compensation light speed value calculated in step S1B while the laser scanning data is being output, in order to eliminate systematic ranging errors caused by changes in temperature and humidity and ensure the accuracy of the three-dimensional point cloud coordinates.

[0074] The necessary procedures are as follows:

[0075] Ranging value correction: The original ranging value is obtained from the flight time Δt collected by the laser radar. (c) ref (where c is the reference speed of light); substituting the environmentally compensated speed of light... env Calculate the distance measurement value after compensation: .

[0076] 3D coordinate generation: The compensated distance measurement value d corr With respect to the horizontal angle θ and vertical angle of the laser beam The coordinates (x, y, z) of a 3D point cloud are converted and output to step S2.

[0077] Environmental compensation and denoising processing of 3D point cloud data includes:

[0078] Step S21: Perform system error compensation on the 3D point cloud data based on the environmental parameters and equipment calibration parameters recorded at the scanning time.

[0079] Systematic error compensation refers to the process of eliminating systematic deviations in 3D point clouds through mathematical models or coordinate transformations based on prior parameters of known error sources. This includes compensation for ranging errors caused by the environment and compensation for equipment geometric attitude errors. Equipment calibration parameters refer to the correction parameters characterizing the installation attitude and internal geometric characteristics of the lidar, obtained through the calibration process in steps S11-S12. These parameters include incident angle errors (pitch angle, roll angle, yaw angle deviations), ranging zero-point errors, and geometric distortion correction matrices.

[0080] The necessary procedures are as follows:

[0081] Fine-grained compensation for environmental ranging errors: The ambient temperature T and humidity H recorded in step S1A at the scanning time are read, and the environmental compensation light speed cenv calculated in step S1B is called to verify and compensate the original ranging values ​​of the 3D point cloud; through... Correct the distance values ​​at each point, based on the scanning angle. Recalculate the three-dimensional coordinates (x, y, z) to eliminate the ranging system error caused by temperature and humidity fluctuations.

[0082] Equipment geometric error compensation: Load the incident angle error parameters (pitch angle correction Δα, roll angle correction Δβ, yaw angle correction Δγ) obtained from calibration in step S12 and the ranging zero point offset Δd; construct the rigid transformation matrix R. calib With translation vector t calib Perform coordinate system correction on the point cloud coordinates: , where n is the laser beam direction vector; eliminate systematic geometric errors caused by the installation posture deviation of the street light pole and the internal optomechanical distortion of the equipment, and output the three-dimensional point cloud data after environment-equipment joint compensation.

[0083] Step S22: Divide the three-dimensional point cloud data after system error compensation into spatial voxels of a preset size, and retain the point with the highest density in each voxel.

[0084] In this context, a spatial voxel refers to a regular three-dimensional cubic mesh unit that divides the three-dimensional point cloud space according to a preset size (such as a side length of 0.1m or a volume of 0.1m³). It is used to discretize continuous space and realize the block processing of point cloud data. The point with the highest density refers to the geometric representative point at the location where the point cloud distribution is most dense within a single spatial voxel. It is usually represented by the three-dimensional centroid (arithmetic mean) of all points within the voxel, or by selecting the measured point closest to the centroid as the representative.

[0085] The necessary procedures are as follows:

[0086] 3D voxel mesh generation: Set the voxel side length L (usually 0.1m, balancing data compression rate and road detail preservation requirements), and calculate the extreme values ​​of the point cloud data in the X, Y, and Z directions (x, y, z). min ,y min ,z min ) is used as the origin of the grid; for each point P n (x n ,y n ,z n ), calculate the integer grid index of its corresponding voxel:

[0087] ,in, For floor operations, (i,j,k) represents the spatial location encoding of the point in the 3D mesh, dividing the discrete point cloud data into corresponding voxel units.

[0088] Voxel density peak extraction and representative point preservation: For each non-empty voxel unit, extract the set of three-dimensional coordinates of the N points contained within it. Calculate the geometric centroid (density center) within this voxel:

[0089] The center of mass This refers to the spatial location of the "highest density point" within the voxel; the centroid coordinates are retained as the unique representative point of the voxel, and the remaining N-1 redundant points within the voxel are eliminated (or the original measurement point with the closest Euclidean distance from the centroid within the voxel is selected as the representative to maintain the measurement accuracy), thereby realizing the downsampling conversion of point cloud data from the original high-density sampling to a regular voxel grid, compressing the data volume while retaining the key geometric features of the road surface.

[0090] Step S23: Calculate the average distance between each point and all points in its preset radius neighborhood, and remove outliers whose average distance deviates from the overall mean by three times the standard deviation.

[0091] The preset radius neighborhood refers to a three-dimensional spherical spatial range with the point to be determined as its center and a preset distance r (e.g., 0.2m) as its radius, used to define the set of locally geometrically related points. The average distance is the arithmetic mean of the Euclidean distances between the point to be determined and all its neighboring points within its preset radius neighborhood, used to quantify the spatial distribution relationship between the point and its local neighborhood. Outliers are abnormal points whose average distance deviates significantly from the overall statistical distribution, usually caused by moving objects (vehicles, pedestrians), suspended particles, or measurement noise, and lack local geometric continuity.

[0092] The necessary procedures are as follows:

[0093] Neighborhood search and average distance calculation: For each point P after voxel downsampling in step S22 i A spherical neighborhood is established with a preset radius r = 0.2m, and all neighboring points within the search area are then searched. Calculate the Euclidean distances from point Pi to each neighboring point. And calculate the average neighborhood distance: .

[0094] Overall statistics and anomaly detection: Traverse all N points and calculate the overall mean of the average neighborhood distance. with standard deviation Set a statistical threshold T = μ + 3σ (three standard deviations criterion). When a certain point satisfies... When this happens, the point is determined to be an outlier.

[0095] Outlier removal and data cleaning: will satisfy Outliers are removed from the point cloud set; points with fewer than 3 points in their neighborhood are marked as edge points or removed, and finally clean point cloud data after statistical denoising is output.

[0096] The relative displacement is checked for consistency based on the spatial neighborhood relationship of each monitoring point. Relative displacements that do not meet the preset consistency criteria are removed, including:

[0097] Step S31: Taking each monitoring point as a node and pairs of monitoring points with a spatial distance less than the preset neighborhood radius as edges, construct a graph structure representing the spatial proximity relationship of the monitoring points.

[0098] Among them, the graph structure refers to the mathematical topological structure G=(V,E) composed of a node set and an edge set, which is used to represent the spatial proximity relationship between monitoring points and the spatial propagation path of pavement deformation.

[0099] A node refers to a vertex in the graph, corresponding to each monitoring point determined in step S3. Each node carries the unique identifier, three-dimensional spatial coordinates, and displacement characteristic attributes of the monitoring point. An edge refers to the topological relationship connecting two nodes, which is established if and only if the spatial distance between two monitoring points is less than the preset neighborhood radius, representing the geometric continuity and deformation correlation of a local area of the pavement. The preset neighborhood radius refers to the distance threshold for constructing the spatial adjacency relationship (usually set to 2m to 5m), which is comprehensively determined based on the thickness of the pavement structure layer, the influence range of soil body deformation, and the layout spacing of monitoring points.

[0100] The necessary processes are as follows:

[0101] Node initialization and attribute assignment: From the set of monitoring points obtained in step S3, establish a graph node for each monitoring point , and assign node attributes including: the unique ID of the monitoring point, three-dimensional spatial coordinates (x i , y i , z i ), relative displacement Δd i ]> and its horizontal components (Δx i , Δy i ) and settlement component Δz i .

[0102] Spatial neighborhood determination and edge construction: Set the preset neighborhood radius R (such as 3m), and traverse all pairs of monitoring points; for any two monitoring points P i and P j , calculate their Euclidean distance in the horizontal plane (ignoring the elevation difference to reflect the lateral influence range of the pavement); if d ij < R, then establish an undirected edge i between nodes v j and v , and assign an edge weight (Gaussian kernel function, σ is the bandwidth parameter) or , indicating that the closer the spatial distance between monitoring points, the stronger the deformation correlation.

[0103] Graph structure data generation: Construct an adjacency matrix (N is the total number of monitoring points), where A ij = w ij If there exists an edge e ijOtherwise A ij =0; Construct the degree matrix D(diagonal matrix, D ii =∑ j A ij The generated graph structure is G=(V,E,A,X), where the node feature matrix is... This provides data input for the graph convolutional network neighborhood aggregation in step S32.

[0104] Step S32: Use a graph convolutional network to perform neighborhood aggregation on the graph structure. Neighborhood aggregation takes the relative displacement of each monitoring point and its horizontal and settlement components as input features. Determine the aggregation weight based on the spatial distance and displacement similarity between the monitoring point and its neighboring nodes. Then, based on the aggregation weight, perform a weighted summation of the relative displacement of each monitoring point in the neighborhood to obtain the aggregated displacement feature corresponding to each monitoring point.

[0105] Graph convolutional networks (GCNNs) are neural network mechanisms that perform neighborhood feature aggregation on a graph structure. They generate new node feature representations by fusing feature information from the central node and its spatial neighbors. Neighborhood aggregation refers to the process of weighted fusion of the displacement features of the central monitoring point and its neighbors within the graph structure, aiming to generate robust displacement estimates by utilizing the spatial continuity constraints of road deformation. Aggregation weights are weighting coefficients assigned to the central node and its neighbors, determined by a spatial distance factor (negatively correlated with geometric distance) and a displacement similarity factor (negatively correlated with displacement feature differences), characterizing the influence of neighboring nodes on the deformation state of the central node. Aggregated displacement features are new feature vectors generated after weighted fusion, comprehensively reflecting the collaborative estimation of the central monitoring point's own displacement and the consistent deformation trend of its spatial neighborhood.

[0106] The necessary procedures are as follows:

[0107] Input feature retrieval: Based on the graph structure G=(V,E,A,X) constructed in step S31, the node feature matrix X is retrieved, where each row xi=[Δd i ,Δx i ,Δy i ,Δz i The relative displacement of monitoring point i and its horizontal and settlement components.

[0108] Composite aggregation weight determination: For each neighboring node of the central node vi (including its own node v) i ), calculate the composite aggregation weights:

[0109] Spatial distance weight: The edge weights constructed directly in step S31 (Based on spatial distance d) ijThe Gaussian kernel function or its reciprocal function indicates that the closer the distance, the stronger the correlation.

[0110] Displacement similarity weight: Euclidean distance calculation based on node feature vectors ,in This is a displacement tolerance parameter, characterizing the sensitivity to displacement differences; the more similar the displacement characteristics, the greater the weight.

[0111] Weight normalization: Calculating composite weights And perform normalization processing. This ensures that the sum of weights within the neighborhood is 1.

[0112] Weighted summation and feature generation: The displacement features of each monitoring point in the neighborhood are weighted and summed based on normalized weights to generate the central node v. i The aggregate displacement characteristic h i :

[0113] The formula weights and fuses the relative displacement and the horizontal and settlement components separately to generate a vector. This means that the aggregated displacement features that integrate spatial neighborhood consistency constraints are passed to step S33 as a benchmark reference value for consistency verification.

[0114] The graph convolutional network used in this application is a neural network architecture based on spatial domain information transfer. Its structure includes an input layer, at least one graph convolutional layer, and a feature aggregation layer. The feature propagation function of the graph convolutional layer is as follows:

[0115] In the formula, σ is the output feature matrix of the (l+1)th graph convolutional layer; σ is the non-linear activation function (such as ReLU); D is the degree matrix with self-loops, D = D + I, where D is the degree matrix and I is the identity matrix; Ã is the adjacency matrix with self-loops, Ã = A + I, where A is the adjacency matrix constructed in step S31; Let l be the input feature matrix of the l-th layer; Let be the learnable weight matrix of the l-th layer.

[0116] This architecture uses a symmetric normalization transformation of the normalized adjacency matrix à and the degree matrix D to weightedly fuse the displacement features of the monitoring point and its spatial neighboring nodes, thereby achieving feature propagation under spatial consistency constraints.

[0117] The training data comes from the historical relative displacement dataset of each monitoring point calculated in step S3. The network training adopts a semi-supervised or unsupervised learning method, and by minimizing the aggregation feature reconstruction error or node classification loss, the network weights W^(l) can optimally aggregate neighborhood information. During operation, the model takes the graph structure constructed in step S31 and the relative displacement of each node as input, and outputs the aggregated displacement features of each monitoring point after spatial consistency constraints.

[0118] This network is entirely dedicated to the technical purpose of eliminating spatial anomalies. Its inputs, outputs, and internal structure are all based on graph topology data mapping and do not involve any non-technical rules. Those skilled in the art can build and use this network based on the above description.

[0119] Step S33: Calculate the absolute value of the residual between the relative displacement of each monitoring point and the aggregate displacement characteristic. When the absolute value of the residual exceeds the preset residual threshold, the relative displacement is judged as abnormal and removed.

[0120] The absolute value of the residual refers to the difference between the original relative displacement of the monitoring point and the displacement feature generated by the graph convolutional network, used to quantify the degree of inconsistency between the monitoring point and the overall deformation trend of the spatial neighborhood. The preset residual threshold is a threshold parameter for judging displacement anomalies, denoted as δ. th (e.g., 5mm to 8mm), determined comprehensively based on the statistical characteristics of registration error and the physical constraints of the continuity of local road deformation in step S3. Anomaly detection and removal refers to marking monitoring points with residuals exceeding a preset threshold as outliers and removing them from the input dataset for subsequent displacement analysis and trend prediction, in order to avoid false deformation signals caused by registration errors, local noise, or interference from moving objects.

[0121] The necessary process is described below:

[0122] Residual calculation: Receive the original relative displacement x of each monitoring point output in step S3. i =[Δd i ,Δx i ,Δy i ,Δz i ] and the aggregate displacement feature h generated in step S32 i =[h d,i ,h x,i ,h y,i ,h z,i ]; Calculate the Euclidean norm (absolute value of the residual vector): Alternatively, a simplified form can be used to calculate only the settlement component: r i =∣Δz i -h z,i |

[0123] Threshold determination and anomaly labeling: The absolute value of the residual r... i Compared with the preset residual threshold δ th Comparison; when r i >δ th When the displacement of a monitoring point deviates significantly from the consistency trend of its spatial neighborhood, it is marked as an anomaly. Threshold δ th Typically, the registration error in step S3 is taken as 2 to 3 times (e.g., 6 mm).

[0124] Anomaly Removal and Data Output: Monitoring points marked as anomalies are removed from the dataset and do not participate in the displacement rate calculation and trend prediction in step S4; alternatively, neighborhood interpolation is used to estimate their reasonable displacement values ​​to ensure the spatiotemporal continuity of the monitoring network; the clean displacement dataset after removal is output to step S4.

[0125] After constructing a graph structure representing the spatial proximity relationships of monitoring points, and before performing neighborhood aggregation on the graph structure using a graph convolutional network, the following steps are also included:

[0126] Step S3A: Obtain geological attribute data and historical settlement data for each monitoring point. The geological attribute data includes at least two of the following: soil layer type, soil compression coefficient, and shear strength parameters. The historical settlement data includes the cumulative settlement time series within a preset historical period.

[0127] Among them, geological attribute data refers to parameters characterizing the engineering properties of the soil at the monitoring point, including at least two of the following: soil layer type, soil compression coefficient, and shear strength parameters. Historical settlement data refers to the time series data of cumulative settlement within a preset historical period.

[0128] The necessary process is described as follows: extract the soil layer type, compression coefficient, and shear strength parameters of each monitoring point from the geological survey database; extract the cumulative settlement time series of the most recent N sampling periods of each monitoring point from the time series data set in step S2; establish a mapping table between monitoring point ID and geological attributes and historical settlement data, normalize the continuous variables, and output to step S3B.

[0129] Step S3B: Based on geological attribute data, a clustering algorithm is used to divide the monitoring points into different geological zones, and monitoring points with geological attribute similarity exceeding a first preset threshold are divided into the same geological zone.

[0130] Geological zoning refers to the category into which monitoring points with similar geological and engineering characteristics are grouped together. The first preset threshold is the similarity threshold (e.g., Euclidean distance of 0.5) for determining whether monitoring points belong to the same geological zoning.

[0131] The necessary steps are as follows: Based on the normalized geological attribute feature vectors, clustering is performed using the K-means algorithm (with a preset number of clusters K) or a hierarchical clustering algorithm; the Euclidean distance between monitoring points of geological attributes is calculated. Monitoring points with a distance less than the first preset threshold θ1 are divided into the same geological zone and assigned a zone label g. i Generate a geological zoning label mapping table and output it to step S3E.

[0132] Step S3C: Real-time collection of environmental load data at each monitoring point location. The environmental load data includes at least one of the following: surface traffic vibration intensity, rainfall, and energy of nearby construction disturbance.

[0133] The necessary process is as follows: Traffic vibration (unit: m / s²) is collected in real time using accelerometers, rain gauges, or construction blasting monitoring instruments deployed around the monitoring points. 2 At least one of the following: rainfall (in mm / h), construction disturbance energy (in J); data are collected to a processing center via wired or wireless network, time-aligned and spatially matched, and an environmental load vector e is generated for each monitoring point. i =[e vib ,e rain ,e cons Output to step S3D.

[0134] Step S3D involves inputting environmental load data into a pre-trained load influence prediction model and outputting the expected disturbance coefficient of each monitoring point under the current environmental load. The expected disturbance coefficient characterizes the degree of disturbance of the monitoring point displacement by the environmental load and has a value range of [0,1].

[0135] The necessary process is as follows: The environmental load data e collected in step S3C... i The load impact prediction model is input into a pre-trained model; the model learns the environment-displacement sensitivity pattern based on historical training data and outputs the expected disturbance coefficients for each monitoring point. The closer the value is to 1, the more significant the environmental load disturbance effect; the expected disturbance coefficient sequence {u i} Bind to the monitoring point ID and output to step S3E.

[0136] Step S3E: Based on the expected disturbance coefficient, combined with the geological zoning attributes and spatial distance of each monitoring point, calculate the dynamic correlation weight between each monitoring point. The dynamic correlation weight is negatively correlated with spatial distance, positively correlated with the membership degree of belonging to the same geological zoning, and positively correlated with the similarity of the expected disturbance coefficient.

[0137] The dynamic association weight refers to the composite weight of spatial distance, geological zoning membership, and environmental disturbance similarity. Membership refers to the degree to which a monitoring point belongs to a specific geological zoning (1 for the same zoning, 0 for different zonings).

[0138] The necessary procedures are as follows:

[0139] Geological zoning labels based on step S3B i The expected perturbation coefficient u in step S3D i and the spatial distance d in step S31 ij Calculate the dynamic association weights:

[0140] .

[0141] Where α, β, and γ are weighting coefficients (satisfying α + β + γ = 1), and σ is the spatial distance attenuation scale. This is a geological zoning indicator function; the weights are locally normalized. Generate a dynamic weight matrix and output it to step S3F.

[0142] Step S3F: Based on the dynamic association weight, the edge weights of the graph structure are dynamically updated. Edges with dynamic association weights lower than the second preset threshold are removed. Virtual edges are added between non-adjacent monitoring points with dynamic association weights higher than the third preset threshold and spatial distances greater than the preset neighborhood radius, thereby generating a dynamic adaptive graph structure.

[0143] Edge pruning refers to deleting edges whose dynamic association weight is below a threshold. Virtual edges are logical edges connecting non-adjacent monitoring point pairs whose spatial distance is greater than a preset neighborhood radius but whose dynamic association weight is extremely high. The dynamic adaptive graph structure refers to the graph updated after pruning and addition operations. .

[0144] Set a second preset threshold θ2 (e.g., 0.3), traverse the initial graph edge set, and set w... ij Edges less than θ2 are pruned; a third preset threshold θ3 (e.g., 0.8) is set, and all non-adjacent monitoring point pairs (spatial distance d) are traversed. ij >R neighbor ), will w ij Add virtual edges between point pairs of θ3 Based on the remaining edges after pruning and the newly added virtual edges, construct the updated adjacency matrix A. dynamic Generate a dynamic adaptive graph structure and output it to step S3G.

[0145] In step S3G, the dynamic adaptive graph structure and the updated edge weights are used as input to the graph convolutional network for subsequent neighborhood aggregation. During the neighborhood aggregation process, the relative displacement of the monitoring points is weighted and summed according to the dynamic association weights to generate aggregated displacement features.

[0146] In this context, "dynamic adaptive graph structure as input" refers to using the updated graph topology and dynamic weights as the computation graph for the graph convolutional network. "Weighted summation" refers to using dynamically associated weights to weight and fuse the features of neighboring nodes.

[0147] The necessary procedures are as follows:

[0148] The dynamic adaptive graph structure generated in step S3F (including the adjacency matrix A) dynamic With weight matrix W dynamic Input to a graph convolutional network; for each monitoring point node v i Aggregate the features of its neighborhood (including virtual neighborhood): .

[0149] Where, x j =[Δd j ,Δx j ,Δy j ,Δz j ] represents the relative displacement and its components of the monitoring point determined in step S3, w ij The dynamic correlation weights are calculated in step S3E; the aggregated displacement feature h, which integrates geological, environmental, and spatial correlations, is generated. i Output to step S33.

[0150] By inputting the displacement change rate and cumulative displacement into a long short-term memory network, the displacement evolution trend over a predetermined time period is predicted, including:

[0151] Step S41: The displacement change rate and cumulative displacement are used as time-series input features and input into the long short-term memory network to generate the original displacement prediction value.

[0152] The temporal input features refer to a two-dimensional time series matrix composed of displacement change rates and cumulative displacement quantities organized in the order of sampling time. Rows correspond to historical time steps, and columns correspond to feature dimensions, serving as the input tensor of the Long Short-Term Memory (LSTM) network. The original displacement prediction values ​​refer to the future displacement estimation sequence directly output by the LTM network based on historical temporal features through forward propagation calculation. This reflects the data-driven trend extrapolation results and has not yet been corrected by soil physical constraints.

[0153] The necessary procedures are as follows:

[0154] Construction of the temporal feature matrix: From the clean displacement dataset output in step S33, extract the displacement change rate and cumulative displacement of each monitoring point for the most recent T consecutive sampling periods (time steps, such as the most recent 24 periods); combine the two-dimensional features in chronological order to construct a temporal input feature matrix of shape (T,2), with each row corresponding to a feature vector of one time step [v t D cum (t)]; Normalize the matrix (e.g., Z-score standardization) to eliminate dimensional differences.

[0155] Forward propagation and generation of original predicted values ​​in the LSTM network: The temporal input feature matrix is ​​fed into a pre-trained Long Short-Term Memory (LSTM) network; the network processes temporal dependencies through its memory units, capturing the long-term trend and short-term fluctuation features of displacement changes; a linear transformation is performed through the fully connected output layer to generate a sequence of original displacement prediction values ​​for the next τ sampling periods (prediction step size, such as the next 24 hours or 7 days). : .in, The original displacement prediction vector has a dimension of τ×1, and its elements are... The predicted displacement values ​​corresponding to the next τ sampling periods; W out h is the weight matrix of the fully connected output layer, with dimension τ×H (H is the dimension of the LSTM hidden layer state), used to map the high-dimensional hidden state to the prediction output space; T b is the hidden state vector of the LSTM network at the last input time step T, with dimension H×1, which encodes the complete displacement evolution pattern and temporal dependencies over the past T time steps; out τ×1 is the bias vector of the fully connected output layer, used to adjust the baseline offset of the predicted output.

[0156] The original displacement prediction value sequence The model is generated solely based on data-driven modeling, without considering physical constraints such as soil compression limit and shear strength. The output is then sent to step S42 for physical residual calculation and constraint projection correction.

[0157] The Long Short-Term Memory (LSTM) network constructed in this application is a recurrent neural network suitable for time series prediction. Its structure includes an input layer, two stacked LSTM hidden state layers, and a fully connected output layer. Each LSTM hidden state layer includes forget gate, input gate, and output gate mechanisms, and the hidden state vector has a dimension of H (such as 128 or 256), which is used to capture long-term dependencies and short-term fluctuation characteristics in displacement time series data.

[0158] The training data comes from the clean displacement time series data (displacement change rate and cumulative displacement) of each monitoring point in the monitoring area during historical periods, after consistency verification in step S3. The training objective is to minimize the mean square error between the network-predicted displacement and the actual observed displacement. During operation, the model takes the time series features of the most recent T periods extracted in step S4 as input and outputs a sequence of original displacement prediction values ​​for the next τ periods.

[0159] This network is designed solely for the technical purpose of predicting road surface displacement trends. Its inputs, outputs, and internal structure are all based on time-series data mapping and do not involve any non-technical rules. Those skilled in the art can build and use this network based on the above description.

[0160] Step S42: Calculate the physical residual based on the original displacement prediction value and the physical constraints of soil deformation. The physical constraints of soil deformation include the settlement limit constraint based on the soil compressibility parameter and the shear failure limit constraint based on the soil shear strength parameter.

[0161] Among them, physical residual refers to the degree of deviation of the original displacement prediction value from the feasible boundary of the physical constraint domain of soil deformation, and is used to quantify the degree of anomaly in the prediction results that violate the laws of soil mechanics. Settlement limit constraint refers to the maximum allowable vertical settlement of the pavement determined based on soil compressibility parameters (compression modulus, compression coefficient), reflecting the ultimate compressive deformation capacity of the soil under load. Shear failure limit constraint refers to the maximum allowable horizontal displacement of the pavement determined based on soil shear strength parameters (cohesion, internal friction angle), reflecting the ultimate ability of the soil to resist shear slip failure.

[0162] The necessary procedures are as follows:

[0163] Constraint boundary setting: Based on the subgrade soil type and geotechnical test parameters of the monitoring area, set the physical constraint threshold for soil deformation; determine the settlement limit value z based on the soil compressibility parameter. max (For example, 300mm for cohesive soil and 150mm for sandy soil); Determine the ultimate horizontal displacement u based on soil shear strength parameters and the Mohr-Coulomb criterion. max (e.g., take 50mm or according to c, calculate).

[0164] Displacement component decomposition and residual calculation: Receive the original displacement prediction value from step S41 and decompose it into settlement components. With horizontal displacement components Calculate the deviation from the constraint boundary respectively:

[0165] Settlement residual: .

[0166] Horizontal shear residual: .

[0167] Combined physical residuals (Euclidean distance or weighted combination): .

[0168] Residual output: Generate the physical residual sequence for each prediction time, mark the predicted values ​​that exceed the physical feasible region of the soil and their degree of exceeding the limit, and output them to step S43 for constrained projection correction.

[0169] Step S43 involves performing constrained projection correction on the original displacement prediction values ​​based on the physical residuals to generate the final displacement prediction values ​​that satisfy the physical constraints of soil deformation. For details, please refer to steps S431 to S433, which will not be elaborated here.

[0170] Step S44: Generate the displacement evolution trend based on the final displacement prediction value.

[0171] Among them, the displacement evolution trend refers to the time series curve constructed based on the final displacement prediction value, which characterizes the development and change of road surface displacement over time within a preset period of time, including displacement value, development rate and relative relationship with the warning threshold.

[0172] The necessary procedures are as follows:

[0173] The final displacement prediction value sequence output in step S43 Arrange the data in chronological order to construct a time-series forecast curve; calculate trend characteristic parameters, including the maximum predicted displacement within a preset future time period. Average displacement change rate The timeline curves and characteristic parameters are encapsulated into structured data and output to step S5 for graded early warning determination and risk assessment.

[0174] Constrained projection correction of the original displacement prediction values ​​based on physical residuals includes:

[0175] Step S431: Based on the settlement limit constraint and the shear failure limit constraint, construct the soil deformation feasible region; according to the magnitude of the physical residual, dynamically adjust the boundary relaxation of the soil deformation feasible region to generate adaptive feasible region boundary conditions.

[0176] The feasible region for soil deformation refers to the closed area enclosed by the settlement limit constraint and the shear failure limit constraint in the three-dimensional displacement space, representing the allowable range of soil deformation without violating physical laws. Boundary relaxation is a constraint boundary adjustment coefficient dynamically determined based on the magnitude of the physical residuals, used to quantify the looseness or tightness of the feasible region boundary, balancing physical strictness and predictive flexibility. Adaptive feasible region boundary conditions refer to the time-varying constraint boundaries dynamically adjusted by the boundary relaxation. When the physical residuals are large, they automatically tighten to enforce physical consistency; when the residuals are small, they are moderately relaxed to retain data-driven characteristics.

[0177] The necessary procedures are as follows:

[0178] Feasible domain space construction: Based on the settlement limit value z determined in step S42 max Based on the horizontal displacement limit value umax, construct the feasible deformation domain Ω of the soil in the three-dimensional displacement space (x,y,z): .

[0179] Dynamic calculation of boundary relaxation: Based on the comprehensive physical residual r calculated in step S42 phy Determine the boundary relaxation coefficient γ (range [γ]). min ,1], where γmin (This is the minimum allowable relaxation coefficient, usually taken as 0.8). Where k is the residual sensitivity coefficient (e.g., 0.1 mm). -1 When the physical residual r phy The larger the value (the more serious the deviation of the prediction from physical laws), the smaller the relaxation coefficient γ, and the tighter the boundary.

[0180] Adaptive boundary condition generation: The original constraint boundaries are multiplied by relaxation coefficients to generate adaptive feasible region boundary conditions.

[0181] Adaptive Settlement Limit: .

[0182] Adaptive horizontal displacement limit: .

[0183] Forming a feasible region that dynamically adjusts with the magnitude of the physical residual. The output is sent to step S432 to constrain the projection.

[0184] Step S432: The original displacement prediction value is projected into the adaptive feasible domain boundary condition using a constraint projection algorithm based on graph Laplace regularization. Graph Laplace regularization maintains the spatial continuity of the displacement field during the projection process based on the spatial adjacency relationship between each monitoring point, and generates a preliminary corrected displacement prediction value.

[0185] Among them, graph Laplacian regularization refers to a regularization method that, based on the Laplacian matrix L (defined as the difference between the degree matrix and the adjacency matrix) of the graph structure in step S31, penalizes the displacement differences between adjacent monitoring points in the projection optimization, thus maintaining the spatial continuity of the displacement field. The constrained projection algorithm refers to mapping the original displacement prediction values ​​to the adaptive feasible region Ω generated in step S431. adaptive The mathematical optimization process within the system ensures that the predicted values ​​meet the physical constraints. The preliminary corrected displacement prediction value refers to the displacement estimate that, after constraint projection, is both within the physical feasible region and maintains spatial continuity, and serves as the initial value for the iterative correction in step S433.

[0186] The necessary procedures are as follows:

[0187] Graph Laplacian Matrix Construction: Using the adjacency matrix A and degree matrix D constructed in step S31, calculate the graph Laplacian matrix. (N is the number of monitoring points), representing the spatial correlation strength between monitoring points.

[0188] Constrained Projection Optimization Solution: Establish a constrained projection optimization problem with graphical Laplacian regularization, and transform the original displacement prediction values... Project onto the feasible region:

[0189] .

[0190] .

[0191] in, Let λ be the initial corrected displacement prediction value to be determined, and λ be the regularization coefficient (usually taken as 0.01 to 0.1) used to balance projection accuracy and spatial continuity; the first term of the objective function ensures that the projected value is close to the original prediction, and the second term... The displacement difference between adjacent monitoring points is penalized, and the displacement field is forced to be smooth.

[0192] Feasible region projection and spatial continuity preservation: The above optimization problem is solved using quadratic programming or the alternating direction multiplier method (ADMM); when the original predicted value exceeds the adaptive feasible region, it is projected to the nearest point on the boundary; at the same time, the Laplace regularization term ensures that the corrected displacement of adjacent monitoring points has spatial consistency; and preliminary corrected displacement prediction values ​​that satisfy physical constraints and are spatially continuous are generated. Output to step S433.

[0193] Step S433: Calculate the correction residual between the preliminary corrected displacement prediction value and the original displacement prediction value. When the correction residual is greater than the preset convergence threshold, backpropagate the correction residual to the hidden state layer of the long short-term memory network, update the internal parameters of the network and regenerate the original displacement prediction value until the convergence condition is met, and obtain the final displacement prediction value.

[0194] The corrected residual refers to the deviation between the preliminary corrected displacement prediction and the original displacement prediction, used to quantify the adjustment magnitude of the constraint projection on the original prediction. The preset convergence threshold is the error threshold parameter used to determine iterative convergence (e.g., 1 mm or 5% relative error); iteration terminates when the corrected residual is less than this threshold. Backpropagation refers to the optimization process of transferring the gradient information of the corrected residual from the output layer to the hidden state layer of the Long Short-Term Memory (LSTM) network to update the network's internal parameters. The hidden state layer refers to the network layer in the LSM network containing the hidden state vector hT that encodes historical temporal information, carrying the model's memory and inference capabilities. The final displacement prediction value refers to the displacement prediction sequence obtained after iterative optimization, which simultaneously satisfies physical constraints and data distribution.

[0195] The necessary procedures are as follows:

[0196] Corrected residual calculation: Receive the preliminary corrected displacement prediction value generated in step S432. Compared with the original displacement prediction value in step S41 Calculate the corrected residuals: Or use the maximum absolute error , where N is the number of monitoring points and τ is the prediction step size.

[0197] Convergence judgment and backpropagation: The corrected residual rcorr With the preset convergence threshold δ conv (e.g., 2mm) comparison; if r corr >δ conv Calculate the residual gradient The gradient is propagated to the hidden state layer h of the Long Short-Term Memory network through backpropagation. T And the network weight parameter W, update the network parameters using the Adam or SGD optimizer: , where η is the learning rate (e.g., 0.001).

[0198] Iterative Update and Re-prediction: Using the updated network parameters, re-execute the forward propagation in step S41 to generate new original displacement prediction values. The new original predicted values ​​are input again into step S432 for constrained projection correction, generating new preliminary corrected values. .

[0199] Convergence determination and final output: Repeat the above process until the corrected residual satisfies r. corr ≤δ conv Or, when the maximum number of iterations is reached (e.g., 10), the preliminary corrected displacement prediction value that meets the convergence condition at this point is taken as the final displacement prediction value x. final Output to step S44.

[0200] Based on the adjusted thresholds, relative displacement, and displacement evolution trends, a tiered early warning signal is generated and sent to the emergency response center and the urban disaster center, including:

[0201] Step S51: Based on the predicted displacement value and the predicted displacement rate of change in the displacement evolution trend, calculate the displacement acceleration factor. The displacement acceleration factor is determined by the time-series difference of the predicted displacement rate of change.

[0202] The displacement acceleration factor refers to the rate of change of displacement over a predetermined time period. It is calculated by time-series difference of the predicted displacement rate and is used to characterize the acceleration or deceleration trend of pavement deformation. A positive value indicates accelerated settlement, while a negative value indicates slowed deformation. The predicted displacement rate refers to the future displacement change per unit time calculated by displacement difference between adjacent time periods based on the final displacement prediction value generated in step S44, reflecting the rate of deformation development in each future time period.

[0203] Temporal differencing is a mathematical method that performs successive differences on adjacent sampling points of a time series, and obtains an acceleration measure by calculating the difference quotient of the rates of change between adjacent time points.

[0204] The necessary procedures are as follows:

[0205] Calculation of displacement rate of change prediction: Obtain the sequence of final displacement prediction values ​​within a preset time period from step S44. (where τ is the prediction step size, d) t+k (where t is the current time, and t is the predicted displacement value for the k-th future period). Calculate the predicted displacement rate of change between adjacent periods:

[0206] Where ΔT is the sampling period duration (in hours), v t+k It represents the average deformation rate within the k-th future period.

[0207] Displacement acceleration factor calculation: Perform time-series differencing on the predicted displacement rate of change sequence to calculate the displacement acceleration factor. A positive value and an increasing absolute value of this factor indicate that pavement settlement has entered an accelerated development stage, while a negative value indicates that the deformation rate has slowed down and is tending to stabilize. The calculated displacement acceleration factor sequence... The output is sent to step S52 as a key component of the multi-dimensional risk assessment feature vector, used to identify the accelerating deterioration trend of collapse risk.

[0208] Step S52: Combine the deviation of displacement acceleration factor, relative displacement amount from the adjusted threshold, and compliance assessment results of soil deformation physical constraints to construct a multi-dimensional risk assessment feature vector; the compliance assessment results are determined based on the historical distribution statistical characteristics of physical residuals.

[0209] Here, the deviation degree refers to the relative difference between the current relative displacement of the monitoring point and the threshold adjusted in step S5, used to quantify the urgency of the current deformation state exceeding the safety boundary. The conformity assessment result refers to a quantitative index evaluating the consistency between the current prediction result and the physical laws of the soil, based on the historical distribution statistical characteristics (such as mean, standard deviation, or quantiles) of the physical residuals from step S42. A higher value indicates a better conformity between the prediction and the laws of soil mechanics. The multi-dimensional risk assessment feature vector refers to a multi-dimensional feature set that integrates displacement dynamics, threshold deviation degree, and physical law conformity, used as input for Bayesian network inference to comprehensively determine the collapse risk.

[0210] The necessary procedures are as follows:

[0211] Feature calculation for each dimension:

[0212] Displacement acceleration factor: directly using a calculated in step S51 t+k This characterizes the trend of accelerated deformation;

[0213] Threshold deviation: Calculate the relative displacement d current With the adjusted threshold S adj deviation ratio r dev =dcurrent / S adj Or calculate the safety margin m=(S adj -d current ) / S adj ;

[0214] Physical constraint compliance: statistical historical physical residuals The distribution of the physical residual r is used to calculate the current physical residual r. phy The cumulative distribution function value or standard score z=(r phy -μ hist ) / σ hist Converted into compliance score (Φ is the standard normal CDF), with a range of [0,1]. The closer to 1, the higher the degree of conformity.

[0215] Feature vector construction: The above features are concatenated according to their dimensions to construct a multi-dimensional risk assessment feature vector. or extended form , where v current The current displacement change rate is used; the features are normalized (e.g., Min-Max normalization to [0,1]) to eliminate dimensional differences, and the output is sent to step S53 as the input variable of the Bayesian network.

[0216] Step S53: Input the multi-dimensional risk assessment feature vector into the Bayesian network inference model, calculate the collapse risk probability and its posterior probability distribution of each monitoring point through probabilistic inference, and extract the uncertainty confidence interval based on the posterior probability distribution.

[0217] The Bayesian network inference model refers to a probabilistic graphical model based on a directed acyclic graph (DAG) representing the dependency between risk features and collapse states, quantifying uncertainty propagation through node conditional probabilities. The collapse risk probability refers to the posterior probability value of a monitoring point collapsing given observed features. The posterior probability distribution refers to the conditional probability distribution of collapse risk obtained by fusing prior knowledge and observational data. The uncertainty confidence interval refers to the range of risk probability estimates extracted based on the posterior distribution (e.g., 95% confidence level).

[0218] The necessary procedures are as follows:

[0219] Evidence Input and Model Initialization: Construct a Bayesian network containing input layer nodes (corresponding to the dimensions of the multi-dimensional risk assessment feature vector in step S52: displacement acceleration factor, threshold deviation, and physical constraint compliance) and output layer nodes (collapse risk state); discretize the risk state into three levels: {low risk, medium risk, high risk} or define it as a continuous probability variable; use the feature vector generated in step S52 as the input node of the observation evidence input network to activate evidence propagation.

[0220] Probabilistic reasoning and posterior calculation: Probabilistic reasoning is performed using the belief propagation algorithm or Markov chain Monte Carlo (MCMC) sampling; the posterior probability of collapse risk under given features is calculated based on Bayes' theorem.

[0221] Where P(R) is the prior probability of risk and P(f|R) is the likelihood function; calculate and output the posterior probability of each monitoring point belonging to each risk level, and extract the posterior probability value p=P(R=high risk|f) of the high risk level as the collapse risk probability, and obtain the complete posterior probability distribution P(R|f).

[0222] Uncertainty confidence interval extraction: Calculate the uncertainty confidence interval based on the posterior probability distribution P(R|f); for discrete risk levels, calculate the confidence range of the probability for each level; for continuous risk probability values, extract the upper and lower quantiles [q] at the 95% confidence level based on the posterior cumulative distribution function. 0.025 ,q 0.975 ] Calculate the confidence interval width Δ=q 0.975 -q 0.025 Output the collapse risk probability p, posterior distribution, and confidence interval width Δ of each monitoring point to step S54.

[0223] The Bayesian network inference model constructed in this application is a probabilistic graphical model based on a directed acyclic graph. Its structure includes an input node layer (corresponding to the dimensions of the multi-dimensional risk assessment feature vector in step S52: displacement acceleration factor, threshold deviation, and physical constraint compliance) and an output node layer (collapse risk state, discretized into {high risk, medium risk, low risk}). Nodes are connected by directed edges, representing the causal dependency between risk features and collapse states. The initial values ​​of the conditional probability table (CPT) for each node are set based on historical collapse case data statistics and prior knowledge of geotechnical engineering.

[0224] Model training is equivalent to its parameter learning process. It utilizes a complete historical case dataset containing features and known risk level labels, learning and updating CPT parameters through maximum likelihood estimation or maximum a posteriori probability estimation. During operation, the model takes the multi-dimensional risk assessment feature vector generated in step S52 as input, performs probabilistic inference through a belief propagation algorithm or Markov chain Monte Carlo sampling, and outputs the posterior probability distribution of collapse risk at each monitoring point and its uncertainty confidence interval.

[0225] This model serves the technical purpose of quantifying the uncertainty of collapse risk. Its inputs, outputs, and internal structure are all based on probabilistic reasoning and do not involve any non-technical rules. Those in the relevant technical field can build and use the model based on the above description.

[0226] Step S54: Based on the combination relationship between the collapse risk probability and the uncertainty confidence interval, dynamically classify the early warning level: when the collapse risk probability exceeds the high risk threshold and the uncertainty confidence interval width is lower than the preset reliability threshold, a level one early warning signal is triggered; when the collapse risk probability is in the medium risk interval or the uncertainty confidence interval width exceeds the preset reliability threshold, a level two early warning signal is triggered.

[0227] The high-risk threshold refers to the probability threshold for determining the collapse risk in a Level 1 warning (e.g., 0.7 or 0.8). When the risk probability exceeds this value, it indicates an extremely high probability of collapse. The preset reliability threshold refers to the confidence interval width threshold for determining the credibility of the prediction results (e.g., 0.1 or 0.15). When the interval width is lower than this value, the uncertainty is considered low and the prediction is reliable. The medium-risk interval refers to the interval where the risk probability is between low and high risk (e.g., [0.3, 0.7]), indicating that there is a potential collapse risk but it has not yet reached an emergency level. The combination relationship refers to the logical combination (AND / OR relationship) of the collapse risk probability and the uncertainty confidence interval width, used to distinguish between high-certainty high-risk states and low-certainty suspicious states.

[0228] The necessary procedures are as follows:

[0229] Threshold parameter setting: Based on historical collapse case statistics and engineering experience, a high-risk threshold p is set. high (e.g., 0.8), medium-risk range [p] low ,p high (e.g., [0.3, 0.8)) and preset reliability threshold Δ reliable (e.g., 0.15); where the reliability threshold represents the maximum acceptable level of uncertainty, and predictions below this value are considered reliable.

[0230] Combinatorial condition determination and dynamic level classification: Receive the collapse risk probability p and uncertainty confidence interval width Δ of each monitoring point output from step S53, and perform combinational logic determination:

[0231] Level 1 warning signal: When p > p high And Δ < Δ reliable The timing of the event indicates a high risk of collapse and a reliable prediction result, requiring immediate emergency response.

[0232] Level II warning signal: When (p low ≤p≤p high ) or (Δ≥Δ reliable When triggered, the former indicates a medium risk that requires enhanced monitoring, while the latter indicates that although the risk may not be high, the prediction uncertainty is large and requires vigilance and verification.

[0233] No warning: When p <p low And Δ < Δ reliable At that time, the risk is low and certain, and no warning is triggered.

[0234] Warning status marking and output: Mark the warning level (Level 1 / Level 2 / None) of each monitoring point in the monitoring point data structure, and output it along with the risk probability value, confidence interval width and judgment basis to step S55 to generate a structured warning message of the corresponding level.

[0235] Step S55: Based on the warning level and multi-dimensional risk assessment feature vector, generate a structured warning message corresponding to the warning level. The structured warning message includes risk location, impact range assessment and on-site handling suggestions.

[0236] Structured early warning messages refer to formatted messages (such as JSON or XML) containing standardized data fields, used for lossless transmission of risk information between monitoring platforms and emergency command systems. Risk location refers to the geographical coordinates (such as WGS84 latitude and longitude or CGCS2000 coordinates) and road marker information of the monitoring point or risk center that triggers the early warning. Impact range assessment refers to the estimation of the length of road sections, number of lanes, and surrounding facilities potentially affected by the collapse, based on risk propagation models or spatial neighborhood relationships. On-site handling recommendations refer to standardized emergency response measures (such as road closures, speed limit patrols, and personnel evacuation) automatically matched according to the early warning level.

[0237] The necessary procedures are as follows:

[0238] Risk location information extraction and encoding: Based on the unique ID of the early warning monitoring point marked in step S54, extract its precise three-dimensional coordinates (longitude, latitude, elevation) and the name of the road and station number from the monitoring point attribute database; encode the location information into a standardized geographic data format (such as GeoJSON point features or WKT text strings) to ensure that the GIS system can directly parse and locate it.

[0239] Impact range assessment calculation: Based on the spatial neighborhood graph structure constructed in step S31, the monitoring point that triggers the warning is used as the seed node. A breadth-first search (BFS) is performed along the graph edges to include the spatial neighborhood monitoring points within the preset distance threshold (e.g., 50 meters for level 1 warning and 30 meters for level 2 warning) into the impact range. The length of the affected road segment, the number of covered lanes, and key surrounding facilities (e.g., underground pipelines, adjacent buildings) are counted to generate a polygonal description or road segment interval code of the spatial impact range.

[0240] Matching and Assembling Disposal Recommendations: Establish a mapping rule between warning levels and disposal recommendation templates (Level 1 warning corresponds to "immediately close roads, prohibit passage, and activate emergency plans," while Level 2 warning corresponds to "speed limit passage, enhanced monitoring, and preparation of emergency resources"); based on the warning level determined in step S54, call the corresponding template, fill in the specific risk location coordinates, impact range parameters, and real-time risk assessment values ​​(such as the predicted maximum settlement); assemble and generate a structured warning message containing the fields {warning level, risk location, impact range, disposal recommendation, generation timestamp, and monitoring point list}, and output it to step S56 for push.

[0241] Step S56: The structured early warning messages are sent in real time to the emergency response center and the urban disaster center through a parallel push link according to the priority order of the early warning level, and the risk status of the monitoring points that trigger the early warning and the related monitoring points in their spatial neighborhood are marked.

[0242] Parallel push links refer to a message distribution mechanism that simultaneously establishes independent data transmission channels to both the emergency response center and the urban disaster center, ensuring that the two command centers receive early warning information synchronously. Priority order refers to the order in which messages are sent according to the early warning level (Level 1 is higher than Level 2), with higher-level early warnings occupying network bandwidth and computing resources first. Risk status labeling refers to marking the current risk level (Level 1 / Level 2) and the early warning timestamp in the monitoring point attribute data structure for GIS visualization and historical tracing. Associated monitoring points refer to neighboring monitoring points that are directly connected to the triggering early warning monitoring point or are within a preset distance threshold, based on the spatial neighborhood graph structure of step S31.

[0243] The necessary procedures are as follows:

[0244] Prioritization and Link Preparation: Receive the structured early warning message queue generated in step S55 and prioritize them according to the early warning level (Level 1 early warning messages are at the front of the queue, followed by Level 2 early warning messages); establish a dual-channel parallel push link, with Channel A connecting to the emergency response center via a 4G / 5G private network or MQTT protocol (focusing on on-site handling), and Channel B connecting to the city disaster center via the government extranet or RESTful API (focusing on macro-level coordination). The two channels are initialized independently and maintain a long connection.

[0245] Parallel push and acknowledgment mechanism: Structured warning messages are sent synchronously through channel A and channel B in priority order, and asynchronous non-blocking I / O is used to ensure real-time performance under high concurrency; a message delivery acknowledgment mechanism is implemented, requiring the receiver to return an ACK acknowledgment message within 30 seconds. If no acknowledgment is received within the time limit, a retransmission strategy is triggered (up to 3 retries). If it still fails, the system switches to the SMS / telephone backup link; message sending timestamp, receiver acknowledgment status and transmission delay indicators are recorded.

[0246] Risk status labeling and neighborhood diffusion: For monitoring points that receive the early warning triggered in step S54, the risk status field is marked with the corresponding level (level 1 or level 2) in their data structure and the current timestamp is written; based on the spatial neighborhood map constructed in step S31, the monitoring point that triggered the early warning is used as the seed node, and its direct neighboring nodes (related monitoring points whose spatial distance is less than the preset neighborhood radius) are traversed. The risk status of these related monitoring points is marked as "affected" or inherits the original early warning level (determined according to the propagation intensity), so as to realize the spatial diffusion labeling of risk status in the monitoring network and provide data support for GIS linkage display and range early warning.

[0247] Risk status labeling for monitoring points that trigger warnings and their associated monitoring points in their spatial neighborhood includes:

[0248] Step S551: Using each monitoring point that triggers the early warning as a seed node, a breadth-first search is performed based on spatial neighborhood relationships to include monitoring points whose uncertainty confidence interval width exceeds the preset confidence threshold and are connected to the seed node into the same clustering area, thereby generating a high uncertainty clustering area.

[0249] In this context, the seed node refers to the monitoring point that triggers the early warning (level one or level two) in step S54, serving as the starting node for identifying risk clustering areas. Breadth-first search is a search algorithm that traverses the graph structure layer by layer along spatial neighborhood relationships, starting from the seed node. The preset confidence threshold is a threshold parameter (e.g., 0.2) for determining that the uncertainty of a monitoring point is too high; when the uncertainty confidence interval width exceeds this value, the prediction reliability is considered insufficient. A high-uncertainty clustering area refers to a spatially continuous region formed by monitoring points connected by the seed node through breadth-first search, where the uncertainty confidence interval width of each point exceeds the preset confidence threshold.

[0250] The necessary procedures are as follows:

[0251] Extract all monitoring points that triggered the early warning marked in step S54 as a set of seed nodes; for each seed node, based on the spatial neighborhood graph constructed in step S31, initialize a breadth-first search queue and include the seed node itself in the current cluster area.

[0252] Take a node from the queue and traverse the unvisited associated monitoring points in its spatial neighborhood; check the uncertainty confidence interval width of each neighborhood point (from step S53). If it exceeds the preset confidence threshold (e.g., Δ>0.2) and is connected to the current seed node, then include the neighborhood point in the same cluster area and add it to the search queue; continue traversing until the queue is empty, forming a high uncertainty cluster area obtained by expanding the seed node.

[0253] Merge and remove duplicates from clustered regions that may overlap between different seed nodes to generate a complete set of high-uncertainty clustered regions; record the list of monitoring points contained in each region, the geometric boundary of the region, and the average uncertainty level, and output them to step S552 for risk propagation weight calculation.

[0254] Step S552: Calculate the risk propagation weight between adjacent monitoring points in the high uncertainty cluster area; the risk propagation weight is positively correlated with the mean of the collapse risk probability and negatively correlated with the mean of the uncertainty confidence interval width.

[0255] Among them, the risk propagation weight is a quantitative indicator of the intensity of risk transmission between adjacent monitoring points within a high uncertainty cluster area, representing the ability of risk to spread from one point to the surrounding area. The mean of the collapse risk probability is the arithmetic mean of the collapse risk probabilities of two adjacent monitoring points, reflecting the overall risk level of the local area. The mean of the uncertainty confidence interval width is the arithmetic mean of the uncertainty confidence interval widths of two adjacent monitoring points, reflecting the reliability of the prediction results.

[0256] The necessary procedures are as follows:

[0257] Adjacent monitoring point pair extraction: Extract all monitoring point pairs (i,j) that have spatial adjacency relationships from the high uncertainty clustering area generated in step S551.

[0258] Risk probability and uncertainty mean calculation: For each pair of adjacent monitoring points, calculate the mean of the collapse risk probability. and the mean of the confidence interval width of uncertainty , where p i ,p j The collapse risk probabilities from step S53, Δi, Δj, are derived from the confidence interval widths from step S53.

[0259] Risk propagation weight calculation: Constructing a risk propagation weight function based on positive and negative correlations: .

[0260] in, To avoid small constants (such as 0.01) dividing by zero, this formula reflects the following: the higher the mean risk probability of adjacent points (the larger the numerator), the greater the risk propagation weight; the higher the mean uncertainty (the larger the denominator), the smaller the propagation weight. Alternatively, an exponential decay form can be used. .

[0261] Where λ is the uncertainty sensitivity coefficient (e.g., 5.0), it also satisfies the constraints of positive correlation between risk probability and negative correlation between uncertainty.

[0262] Weight normalization and output: Locally normalize the risk propagation weights for each pair of adjacent points. Ensure that the risk propagation weights from a single node to all neighboring regions sum to 1; output the normalized risk propagation weight matrix to step S553 for calculating the risk diffusion coefficient.

[0263] Step S553: ​​Based on the risk propagation weight and the spatial distance between adjacent monitoring points, calculate the risk diffusion coefficient of each monitoring point using the graph diffusion kernel function.

[0264] The graph diffusion kernel function is a mathematical function that simulates the attenuation and propagation of risk between monitoring points based on a graph structure. By combining risk propagation weights and spatial distance parameters, it quantifies the attenuation law of risk diffusion from the source point to the neighborhood. The risk diffusion coefficient is a quantitative index that comprehensively characterizes the cumulative effect of the risk of each monitoring point itself and the risk propagated in the neighborhood. The larger the value, the more significant the impact of risk diffusion on that point. The distance attenuation scale parameter is a characteristic length that controls the rate of risk attenuation with spatial distance, denoted as σ (usually taken as 1m to 2m), which is empirically set based on the influence range of road surface soil deformation.

[0265] The necessary procedures are as follows:

[0266] Risk propagation weight calculated based on step S552 Spatial distance between adjacent monitoring points obtained in step S31 Construct a Gaussian graph diffusion kernel function: Where σ is the distance decay scale parameter, controlling the rate at which risk decreases with distance; K ij The kernel function represents the intensity of risk diffusion from monitoring point j to monitoring point i. The closer the distance and the higher the propagation weight, the larger the kernel function value.

[0267] The probability of collapse at each monitoring point calculated in step S53 is p. j As the initial risk source strength, the risk diffusion coefficient c of each monitoring point is calculated by weighted aggregation using the graph diffusion kernel function. i :

[0268] .

[0269] Alternatively, a normalization method can be used to ensure numerical stability: .

[0270] This coefficient integrates the risk inherent to the monitoring point itself with the risk contribution from spatial neighborhood propagation, generating a risk diffusion coefficient sequence {c} for each monitoring point. i The output is sent to step S554 for the generation of continuous risk impact hot zones.

[0271] Step S554: Based on the risk diffusion coefficient of each monitoring point, a continuous risk impact hot zone is generated using a spatial interpolation method. The spatial range of the risk impact hot zone extends outward from the boundary of the high uncertainty cluster area to the location where the risk diffusion coefficient decays to a preset background threshold.

[0272] The necessary procedures are as follows:

[0273] Discrete point data preparation: Extract the risk diffusion coefficient {c} of each monitoring point calculated in step S553. i} and its corresponding spatial coordinates (x) i ,y i ,z i Meanwhile, the boundary polygon of the high uncertainty clustering region determined in step S551 is obtained as the interpolation constraint boundary.

[0274] Spatial Interpolation and Continuous Field Generation: Inverse distance weighted (IDW) or Kriging interpolation methods are used to interpolate the risk diffusion coefficients of discrete monitoring points to a continuous spatial grid (grid resolution such as 0.5m × 0.5m). Taking inverse distance weighted interpolation as an example, the interpolated risk diffusion coefficient at any spatial location P(x,y) is:

[0275] , where d(P,Pi) is the distance from the interpolation point to the i-th monitoring point, β is the distance decay exponent (usually taken as 2), and N is the number of monitoring points participating in the interpolation; interpolation is only performed in the local range extending outward from the boundary of the high uncertainty clustering region in step S551 to reduce the computational load.

[0276] Hot zone boundary determination and range clipping: based on a preset background threshold c bg (e.g., 0.05), after interpolation, the contour lines of the risk diffusion coefficient are c(x,y)=c bg The boundary of the hot zone is used as the outer boundary; the boundary of the high uncertainty aggregation area is used as the inner boundary of the hot zone, forming a continuous transition region from the inside to the outside where the risk diffusion coefficient decreases from the peak value to the background value; the interpolation results are cropped, retaining only c(P)≥c bg Within the spatial range, a closed polygonal thermal zone boundary is generated.

[0277] Risk impact hotspot output: Output continuous risk impact hotspots in the form of vector polygons (such as GeoJSON format) or raster heat maps (such as TIFF format), record the risk diffusion coefficient value at each location, and pass it to step S555 for overlay display with the risk map and linkage display with displacement evolution trend.

[0278] Step S555: Overlay and mark the risk impact hotspots on the risk map, and display them in conjunction with the temporal changes in displacement evolution trends and uncertainty confidence intervals.

[0279] The risk map refers to a base map of urban geographic information that includes road networks, topography, underground pipelines, and buildings, providing spatial location references for risk hotspots. Overlay annotation refers to overlaying the generated risk impact hotspots onto the risk map as a thematic layer, visually displaying the spatial distribution of risks through color gradients or transparency. Linked display refers to an interactive and correlated display mechanism between spatial data (risk hotspots) and temporal data (displacement evolution trends, confidence intervals), supporting the synchronous visualization of multi-dimensional risk information in both spatial and temporal dimensions.

[0280] The necessary procedures are as follows:

[0281] Data layer preparation and loading: Obtain the risk impact hotspot data (vector polygon or raster format) generated in step S554, as well as the displacement evolution trend sequence from step S44 and the uncertainty confidence interval time series data from step S53; In the Geographic Information System (GIS) platform, the risk impact hotspot is overlaid as a dynamic thematic layer on the basic risk map, and the continuous spatial distribution of the risk diffusion coefficient is rendered using color coding (e.g., red - high risk diffusion, orange - medium risk, green - background threshold).

[0282] Spatial-temporal data binding: Establish an index association between the spatial location of risk hotspots and the temporal data of corresponding monitoring points to ensure that any location on the map can be mapped to its displacement evolution trend curve and confidence interval zone; bind the prediction time step as a temporal attribute to the hotspot layer to support dynamic updates by time dimension.

[0283] Interactive and Linked Visualization: Enables two-way linked display function: When a user clicks or hovers over a specific location in a hot zone on the risk map, a pop-up window or sidebar simultaneously displays the displacement evolution trend curve (including future predicted values) and the time-series change band of the uncertainty confidence interval (such as the upper and lower bounds of the 95% confidence interval); or provides a time axis slider control, which dynamically updates the range of the hot zone with the predicted time when dragged, intuitively showing the diffusion process of risk evolution over time, and providing spatiotemporal integrated visualization decision support for emergency command.

[0284] The steps for dynamically adjusting this threshold based on the rate of displacement change include:

[0285] Step S5A: Extract the relative displacement and its horizontal and settlement components of each monitoring point. Combine the physical residuals calculated based on the physical constraints of soil deformation, calculate the approximation degree of the current state of each monitoring point to the settlement limit constraint and the shear failure limit constraint, and use it as the local progressive failure degree of each monitoring point.

[0286] Among them, the local progressive failure degree refers to the degree of deviation of the current displacement state and physical laws of the comprehensive monitoring point. It is a dimensionless index that quantitatively characterizes the relative proximity of the soil from the current state to the ultimate failure state. The value ranges from 0 to 1, and the larger the value, the higher the risk of failure. The approximation degree refers to the ratio of the relative distance between the current displacement and physical residual and the physical constraint boundary of the soil. The relative displacement and its horizontal and settlement components are derived from the calculation results of step S3. The physical residual is derived from the calculation results of step S42 based on the physical constraint conditions of soil deformation, and its magnitude reflects the degree of deviation of the current state from the physical laws. The settlement limit constraint and shear failure limit constraint are derived from the physical boundaries set in step S42 based on the soil compressibility parameters and shear strength parameters.

[0287] The necessary procedures are as follows:

[0288] Displacement state and physical residual extraction: Extract the relative displacement of each monitoring point from step S3. Horizontal displacement components and sedimentation components Extract the comprehensive physical residual corresponding to the monitoring point from step S42. .

[0289] Approximation Calculation: Based on the physical residual, an equivalent penalty is applied to the current displacement to calculate the vertical approximation of the current state of each monitoring point from the settlement limit constraint. Horizontal approximation of the distance from the shear failure limit constraint :

[0290] .

[0291] In the formula, For monitoring points The vertical approximation; For monitoring points The settlement component; This is the vertical physical residual penalty coefficient, with a value greater than 0, used to convert the degree of deviation from physical laws into an equivalent vertical approximation quantity; For monitoring points The comprehensive physical residual; This represents the settlement limit value.

[0292] .

[0293] In the formula, For monitoring points Horizontal approximation; For monitoring points The horizontal displacement component; This is the penalty coefficient for horizontal physical residuals, and its value is greater than 0. This represents the limit value of horizontal displacement.

[0294] Determination of local progressive destructive degree: The maximum value of the vertical approximation degree and the horizontal approximation degree is taken as the local progressive destructive degree. :

[0295] .

[0296] Alternatively, a weighted combination approach can be used:

[0297] .

[0298] In the formula, The weighting coefficient has a range of values. .

[0299] right Perform truncation if Then take ,like Then take Ensure that its value falls within Within the interval; generate a local progressive damage sequence for each monitoring point and output it to step S5B.

[0300] Step S5B: Construct a graph structure based on the spatial neighborhood relationship of the monitoring points, use graph Laplacian regularization to calculate the spatial transmission weight of local asymptotic damage between adjacent monitoring points, and integrate the displacement change rate of each monitoring point to generate the field instability diffusion index of each monitoring point affected by the spatial neighborhood.

[0301] The field instability diffusion index is a composite index that comprehensively characterizes the influence of the transmission of damage state within its spatial neighborhood and the driving effect of local deformation rate on a monitoring point, reflecting the degree of spread of overall regional instability risk. The spatial transmission weight is a coefficient determined based on the principle of graph Laplace regularization, reflecting the intensity of the diffusion and transmission of local progressive damage between adjacent monitoring points; the smaller the difference in damage degree and the tighter the spatial connection, the greater the weight. The graph structure is derived from the topological network representing spatial proximity relationships constructed in step S31. The displacement change rate is derived from the recent displacement development rate index extracted in step S4.

[0302] The necessary procedures are as follows:

[0303] Graph Laplacian Matrix Call: Calls the graph adjacency matrix constructed in step S31. degree matrix Calculate the graph Laplacian matrix :

[0304] .

[0305] In the formula, The graph is a Laplace matrix; For a degree matrix, its diagonal elements ; It is an adjacency matrix.

[0306] Transmission weight calculation: Based on the principle of graph Laplace regularization, calculate the weight of adjacent monitoring points. and The original transmission weight of local progressive damage between :

[0307] .

[0308] In the formula, For monitoring points and The original transmission weights between them; For monitoring points The degree of localized, progressive damage; For monitoring points The degree of localized, progressive damage; The bandwidth parameter for the degree of damage is positive; For adjacency matrix elements, when monitoring points and When the spatial distance is less than the preset neighborhood radius ,otherwise This formula penalizes edges with excessively large differences in destructive power, ensuring that destructive states are smoothly propagated between similar and adjacent nodes.

[0309] Normalization: The original transmission weights are normalized to obtain the monitoring points. Neighboring nodes Spatial transmission weight :

[0310] .

[0311] In the formula, For monitoring points Neighboring nodes Spatial transmission weight; For monitoring points The set of neighboring nodes; Indicates monitoring point The sum of the original propagation weights between the monitoring point and all its neighboring nodes; after normalization, the monitoring point The sum of the neighborhood weights is 1, that is... .

[0312] Calculation of the field instability diffusion index: based on the local progressive degree of damage at each monitoring point. As the source strength, the spatial transmission weight and displacement change rate driving factor are combined to calculate the strength of each monitoring point. Field instability diffusion index :

[0313] .

[0314] In the formula, For monitoring points The field instability diffusion index; For monitoring points The set of neighboring nodes; For monitoring points Neighboring nodes Spatial transmission weight; For monitoring points The degree of localized, progressive damage; The rate excitation coefficient is a positive number. For monitoring points The rate of displacement change; the greater the rate of displacement change, the stronger the driving effect of the region on risk diffusion; generate the field instability diffusion index sequence for each monitoring point and output it to step S5C.

[0315] In step S5C, based on the evolution mechanism of soil and rock disasters, the local progressive damage degree, field instability diffusion index and displacement acceleration factor of each monitoring point are integrated to calculate the dynamic adjustment coefficient for each monitoring point; the displacement acceleration factor is determined based on the time difference of the displacement change rate.

[0316] The dynamic adjustment coefficient is a dimensionless factor used to adaptively reduce the preset basic early warning threshold. Its value is positively correlated with the degree of risk deterioration, and its range is [value range missing]. The evolution mechanism of soil and rock hazards refers to the physical process by which road collapse spreads from local mechanical damage to the spatial field, ultimately accelerating instability over time. This corresponds to three dimensions: mechanical state, spatial transmission, and temporal acceleration. The local progressive damage degree characterizes the degree of deterioration of the mechanical state; the field instability diffusion index characterizes the breadth and intensity of spatial transmission; and the displacement acceleration factor characterizes the accelerating trend of temporal evolution. Its calculation method is consistent with step S51, determined by the temporal difference of the displacement change rate, reflecting whether the deformation has entered the acceleration stage.

[0317] The necessary procedures are as follows:

[0318] Displacement acceleration factor calculation: Based on the displacement change rate time series data extracted in step S4, the monitoring point is calculated by the difference in rate between adjacent time moments. displacement acceleration factor The calculation method is the same as step S51:

[0319] .

[0320] In the formula, For monitoring points Displacement acceleration factor; For monitoring points In the The rate of displacement change per sampling period; For monitoring points In the The rate of displacement change per sampling period; This refers to the sampling period duration; A positive value indicates accelerated deformation, while a negative value indicates deceleration.

[0321] Multi-dimensional feature normalization: reducing local progressive destruction Field instability diffusion index and displacement acceleration factor Min-Max normalization was applied to each of the following methods. Intervals eliminate dimensional differences.

[0322] Local gradual destruction normalization:

[0323] .

[0324] In the formula, For monitoring points Normalized local progressive destruction degree; It represents the minimum value among all the local progressive damage levels at all monitoring points; This represents the maximum value among all monitoring points for localized progressive damage.

[0325] Normalization of the field instability diffusion index:

[0326] .

[0327] In the formula, For monitoring points Normalized field instability diffusion index; It is the minimum value among all the field instability diffusion indices at all monitoring points; This is the maximum value among all monitoring points' field instability diffusion indices.

[0328] Displacement acceleration factor normalization:

[0329] .

[0330] In the formula, For monitoring points Normalized displacement acceleration factor; It is the minimum value among all displacement acceleration factors at all monitoring points; This is the maximum value among all displacement acceleration factors at all monitoring points.

[0331] Dynamic adjustment coefficient calculation: Based on the evolution mechanism of soil and rock hazards, an increase in any of the above three dimensions indicates a sharp rise in risk, and the early warning threshold should be more sensitive; a fusion calculation model based on an exponential function is constructed to calculate the monitoring points. Dynamic adjustment coefficient :

[0332] .

[0333] In the formula, For monitoring points The dynamic adjustment coefficient; , , These are the weighting coefficients for each dimension of the features in the index model, determined based on regression analysis of historical disaster cases and taking positive values. It is a natural exponential function.

[0334] Alternatively, a fusion computation model based on power functions can be used:

[0335] .

[0336] In the formula, , , The weight coefficients for each dimension of the features in the power function model are determined based on regression analysis of historical disaster cases and are set to positive values. , , These are non-linear exponents, all taking positive values; the exponential model and the power function model described above should be used interchangeably. Dynamic adjustment coefficients. The larger the value, the more drastic the evolution of comprehensive risk, and the output is sent to step S5D.

[0337] Step S5D involves using a dynamic adjustment coefficient to nonlinearly reduce the preset basic warning threshold, thereby generating an adaptive dynamic warning threshold corresponding to the spatial location of each monitoring point.

[0338] Nonlinear reduction refers to the process of lowering the basic early warning threshold using a nonlinear decreasing function based on the magnitude of a dynamically adjusted coefficient. This ensures that the threshold drops rapidly when the risk evolves dramatically to improve early warning sensitivity, and remains stable when the risk is stable to avoid false alarms. The preset basic early warning threshold refers to a fixed displacement alarm threshold value set based on engineering experience or specifications. The adaptive dynamic early warning threshold refers to a personalized alarm threshold that is dynamically updated according to the spatial location of each monitoring point and the real-time risk evolution status, replacing the traditional one-size-fits-all fixed threshold.

[0339] The necessary procedures are as follows:

[0340] Basic threshold retrieval: Retrieves the system's preset basic warning threshold. This threshold typically includes a cumulative displacement threshold and a velocity threshold.

[0341] Nonlinear reduction calculation: using the dynamic adjustment coefficient calculated in step S5C An exponential decay function is used to nonlinearly reduce the basic early warning threshold, and the monitoring points are calculated. Adaptive dynamic early warning threshold :

[0342] .

[0343] In the formula, For monitoring points Adaptive dynamic early warning threshold; Basic warning threshold; To reduce the sensitivity coefficient, a positive value is used; It is a natural exponential function.

[0344] Alternatively, the reciprocal function form can be used:

[0345] .

[0346] In the formula, either the exponential decay function or the reciprocal function is used. When the dynamic adjustment coefficient... hour, The threshold remains unchanged from the base value; when hour, It decreases at a nonlinear rate, thereby automatically lowering the alarm threshold when local damage intensifies, field spreads, or deformation accelerates, enabling early detection of high-risk areas.

[0347] Minimum Threshold Constraint: To prevent the adaptive dynamic warning threshold from approaching zero due to excessively large dynamic adjustment coefficients, thus causing false triggering of measurement noise, a preset minimum threshold is set. (Usually, the basic early warning threshold is used) (10% to 20%), and impose a lower limit truncation constraint on the reduced adaptive dynamic early warning threshold:

[0348] .

[0349] In the formula, This is the preset minimum threshold. The function aims to maximize the value of the warning threshold, ensuring that even under extreme risk evolution, the warning threshold remains within a safe boundary that is not lower than the lower limit of the measurement noise identification limit, thus balancing warning sensitivity and noise resistance stability.

[0350] Adaptive threshold output: The set of adaptive dynamic early warning thresholds corresponding to the spatial location of each monitoring point after being subject to lower limit truncation constraints. The static threshold in substitution weight 1 is output for use in the hierarchical early warning judgment logic of subsequent step S5. When the relative displacement or displacement evolution trend exceeds the adaptive threshold, an early warning is triggered.

[0351] Based on the same inventive concept, embodiments of the present invention provide an urban road collapse monitoring and early warning system based on laser point cloud data, including a memory and a processor, wherein the memory stores data that can run on the processor to implement the following... Figure 1 The procedure for the method shown.

[0352] The embodiments described in this specific implementation are preferred embodiments of this application and are not intended to limit the scope of protection of this application. Therefore, all equivalent changes made in accordance with the structure, shape and principle of this application should be covered within the scope of protection of this application.

Claims

1. A method for monitoring and early warning of urban road collapse based on laser point cloud data, characterized in that, include: Three-dimensional point cloud data is obtained by scanning the road surface with laser point cloud scanning equipment deployed on street light poles according to a preset sampling cycle; Environmental compensation and denoising processing are performed on the 3D point cloud data. The processed data is stored according to time index to form a time series data set containing observations of a preset historical period. Based on the three-dimensional point cloud data of different time periods in the time series data set, the coordinate system is unified by the iterative nearest point algorithm, and the monitoring points of the road surface are determined based on the unified three-dimensional point cloud data. The relative displacement and its horizontal and settlement components of each monitoring point are calculated. Based on the spatial neighborhood relationship of each monitoring point, the relative displacement is checked for consistency, and the relative displacement that does not meet the preset consistency standard is eliminated. The displacement change rate of each monitoring point is extracted based on the last preset number of continuous sampling periods in the preset historical period, and the cumulative displacement of each monitoring point is extracted based on all sampling periods. The displacement change rate and cumulative displacement are input into a long short-term memory network to predict the displacement evolution trend over a preset time period. Based on a preset basic warning threshold, the threshold is dynamically adjusted in conjunction with the displacement change rate. Based on the adjusted threshold, relative displacement, and displacement evolution trend, a graded early warning signal is generated and sent to the emergency response center and the urban disaster center.

2. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 1, characterized in that, Before scanning the road surface according to the preset sampling period, it also includes: The laser point cloud scanning device is controlled to scan a calibration plate of preset size and reflectivity to obtain point cloud data of the calibration plate; Geometric features of the calibration board point cloud data are extracted, and the incident angle error of the laser point cloud scanning equipment is calibrated based on the deviation between the geometric features and the theoretical geometric parameters of the calibration board.

3. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 1, characterized in that, The road surface is scanned according to a preset sampling period, including: Real-time acquisition of ambient temperature and humidity at the time of scanning; Based on temperature and humidity, as well as a preset light speed correction formula, calculate the light speed value after environmental compensation. The laser ranging value is compensated in real time based on the speed of light value after environmental compensation.

4. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 1, characterized in that, Environmental compensation and denoising processing of 3D point cloud data includes: Based on the environmental parameters and equipment calibration parameters recorded at the scanning time, systematic error compensation is performed on the 3D point cloud data; The 3D point cloud data after system error compensation is divided into spatial voxels of a preset size, and the point with the highest density is retained in each voxel. Calculate the average distance between each point and all points in its preset radius neighborhood, and remove outliers whose average distance deviates from the overall mean by three times the standard deviation.

5. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 1, characterized in that, The relative displacement is checked for consistency based on the spatial neighborhood relationship of each monitoring point. Relative displacements that do not meet the preset consistency criteria are removed, including: Using each monitoring point as a node and monitoring point pairs whose spatial distance is less than the preset neighborhood radius as edges, a graph structure representing the spatial proximity relationship of monitoring points is constructed. Graph convolutional networks are used to perform neighborhood aggregation on graph structures. Neighborhood aggregation takes the relative displacement of each monitoring point and its horizontal and settlement components as input features. Aggregation weights are determined based on the spatial distance and displacement similarity between the monitoring point and its neighboring nodes. The relative displacements of each monitoring point in the neighborhood are weighted and summed based on the aggregation weights to obtain the aggregated displacement features corresponding to each monitoring point. Calculate the absolute value of the residual between the relative displacement of each monitoring point and the aggregate displacement characteristic. When the absolute value of the residual exceeds the preset residual threshold, the relative displacement is judged as abnormal and removed.

6. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 1, characterized in that, By inputting the displacement change rate and cumulative displacement into a long short-term memory network, the displacement evolution trend over a predetermined time period is predicted, including: The displacement change rate and cumulative displacement are used as time-series input features and fed into a long short-term memory network to generate the original displacement prediction value. The physical residual is calculated based on the original displacement prediction value and the physical constraints of soil deformation. The physical constraints of soil deformation include the settlement limit constraint based on the soil compressibility parameter and the shear failure limit constraint based on the soil shear strength parameter. Based on the physical residual, the original displacement prediction value is corrected by constraint projection to generate the final displacement prediction value that satisfies the physical constraint conditions of soil deformation. The displacement evolution trend is generated based on the final displacement prediction value.

7. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 6, characterized in that, Constrained projection correction of the original displacement prediction values ​​based on physical residuals includes: Based on settlement limit constraints and shear failure limit constraints, a feasible domain for soil deformation is constructed; according to the magnitude of physical residuals, the boundary relaxation of the feasible domain for soil deformation is dynamically adjusted to generate adaptive feasible domain boundary conditions. A constraint projection algorithm based on graph Laplace regularization is adopted to project the original displacement prediction values ​​into the adaptive feasible region boundary conditions. Graph Laplace regularization maintains the spatial continuity of the displacement field during the projection process based on the spatial adjacency relationship between each monitoring point, and generates preliminary corrected displacement prediction values. The correction residual between the preliminary corrected displacement prediction value and the original displacement prediction value is calculated. When the correction residual is greater than the preset convergence threshold, the correction residual is backpropagated to the hidden state layer of the long short-term memory network to update the internal parameters of the network and regenerate the original displacement prediction value until the convergence condition is met, and the final displacement prediction value is obtained.

8. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 5, characterized in that, Based on the adjusted thresholds, relative displacement, and displacement evolution trends, a tiered early warning signal is generated and sent to the emergency response center and the urban disaster center, including: Based on the predicted displacement value and the predicted displacement rate of change in the displacement evolution trend, the displacement acceleration factor is calculated. The displacement acceleration factor is determined by the time difference of the predicted displacement rate of change. By integrating the deviation of displacement acceleration factor, relative displacement amount from the adjusted threshold, and the conformity assessment results of soil deformation physical constraints, a multi-dimensional risk assessment feature vector is constructed; the conformity assessment results are determined based on the historical distribution statistical characteristics of physical residuals. The multi-dimensional risk assessment feature vector is input into the Bayesian network inference model. The collapse risk probability and its posterior probability distribution of each monitoring point are calculated through probabilistic inference. The uncertainty confidence interval is extracted based on the posterior probability distribution. Based on the combination relationship between the collapse risk probability and the uncertainty confidence interval, the early warning level is dynamically divided: when the collapse risk probability exceeds the high-risk threshold and the uncertainty confidence interval width is lower than the preset reliability threshold, a level one early warning signal is triggered; when the collapse risk probability is in the medium-risk range or the uncertainty confidence interval width exceeds the preset reliability threshold, a level two early warning signal is triggered. Based on the warning level and multi-dimensional risk assessment feature vector, a structured warning message corresponding to the warning level is generated. The structured warning message includes risk location, impact range assessment and on-site handling suggestions. Structured early warning messages are sent in real time to the emergency response center and the urban disaster center in parallel push links according to the priority order of the early warning level, and the risk status of the monitoring points that trigger the early warning and the related monitoring points in their spatial neighborhood are marked.

9. The urban road collapse monitoring and early warning method based on laser point cloud data according to claim 8, characterized in that, Risk status is marked for the monitoring points that trigger the early warning and their associated monitoring points in the spatial neighborhood, including: Using each monitoring point that triggers an early warning as a seed node, a breadth-first search is performed based on spatial neighborhood relationships. Monitoring points whose uncertainty confidence interval width exceeds a preset confidence threshold and are connected to the seed node are included in the same cluster area, generating a high uncertainty cluster area. Calculate the risk propagation weight between adjacent monitoring points within a high uncertainty cluster area; the risk propagation weight is positively correlated with the mean of the collapse risk probability and negatively correlated with the mean of the uncertainty confidence interval width; Based on the risk propagation weight and the spatial distance between adjacent monitoring points, the risk diffusion coefficient of each monitoring point is calculated using the graph diffusion kernel function. Based on the risk diffusion coefficient of each monitoring point, a continuous risk impact hot zone is generated using a spatial interpolation method. The spatial range of the risk impact hot zone extends outward from the boundary of the high uncertainty cluster area to the location where the risk diffusion coefficient decays to a preset background threshold. The risk impact hotspots are overlaid and marked on the risk map, and displayed in conjunction with the displacement evolution trend and the time-series changes of the uncertainty confidence interval.

10. A monitoring and early warning system for urban road collapse based on laser point cloud data, characterized in that, The system includes a memory, a processor, and a program stored in the memory and executable on the processor, which, when loaded and executed by the processor, implements a method for monitoring and early warning of urban road collapse based on laser point cloud data as described in any one of claims 1 to 9.

Images