Space-air-ground integrated cooperative geological disaster real-time monitoring and early warning method

Abstract

The application discloses a space-air-ground integrated cooperative geological disaster real-time monitoring and early warning method and belongs to the technical field of geological disaster early warning, and specifically comprises the following steps: collecting displacement and inclination data of a monitoring area in real time, and marking an abnormal node as a primary early warning node; when an early warning node is detected, scheduling a flight device to collect a ground image and analyze crack geometric features; comparing the node displacement direction with the crack extension direction for consistency, and determining the correlation in combination with the node spatial distribution; dividing the nodes with consistency into an early warning cluster, and acquiring the soil water content distribution in the cluster; comprehensively considering the node displacement, crack features and water content data, generating a geological disaster early warning instruction, and sending the instruction to an execution device; and the application realizes automatic fusion and intelligent research and judgment of space-air-ground multi-source data, and significantly improves the accuracy and timeliness of geological disaster early warning.

Description

Technical Field

[0001] This invention relates to the field of geological disaster early warning technology, specifically to a real-time monitoring and early warning method for geological disasters that integrates air, space, and ground technologies. Background Technology

[0002] Geological hazards, especially sudden geological phenomena such as landslides and collapses, have always been a key focus and challenge in global disaster prevention and mitigation efforts due to their suddenness and destructiveness. Traditional monitoring methods rely on manual inspections and rudimentary observations, which have significant limitations in terms of timeliness and coverage. With the development of information and sensor technologies, geological hazard monitoring has entered a stage characterized by automation, continuity, and three-dimensionality. Constructing a comprehensive monitoring system integrating ground sensing, near-space detection, and remote observation to achieve comprehensive perception and early warning of changes in the state of hazard-prone areas has become an inevitable technological path to improve disaster prevention and mitigation capabilities. Technological advancements in this field are directly related to the accuracy of early warning information and the effectiveness of evacuation responses.

[0003] However, existing technologies still face significant challenges in achieving real-time and accurate early warning. Current geological disaster monitoring and early warning mainly rely on single-type sensors or periodic remote sensing image analysis, resulting in response lag and a high false alarm rate. Specifically, while ground-based sensor networks (such as displacement gauges and inclinometers) can achieve continuous data acquisition at the minute or even second level, capturing deformation and tilt changes at monitoring points in real time, their monitoring range is limited, typically only reflecting the local state of the sensor installation point. When single or a few sensors show data anomalies, existing methods struggle to effectively distinguish whether the anomaly is due to local soil and rock loosening, instrument malfunction, or an early precursor signal of large-scale landslide instability, thus failing to accurately assess the spatial extent and overall risk of the disaster. On the other hand, aerial drones or satellite remote sensing technologies can acquire information on surface deformation and crack development at a macroscopic scale, offering strong spatial coverage; however, their observation behavior is usually constrained by preset mission cycles or fixed revisit cycles, making continuous observation difficult. Especially under adverse weather conditions such as clouds, rain, and fog, optical remote sensing methods often fail. Therefore, when a ground-based sensor network triggers an anomaly alarm, the existing technology system cannot ensure that near-ground observation equipment can be immediately mobilized to the scene to obtain high-definition images and 3D data of the alarm time and subsequent critical periods. As a result, it is impossible to quickly and intuitively verify the authenticity of the ground sensor alarm and the spatial distribution of the anomaly. Summary of the Invention

[0004] The purpose of this invention is to provide a real-time monitoring and early warning method for geological disasters that integrates air, space, and ground technologies, and to solve the following technical problems:

[0005] Existing geological disaster early warning methods suffer from a disconnect between ground-based and space-based monitoring in terms of spatiotemporal response. This makes it impossible to automatically trigger and complete rapid collaborative verification of near-ground observations after a ground sensor alarm, resulting in difficulty in distinguishing between local anomalies and regional instability precursors, leading to delayed early warning response and insufficient reliability.

[0006] The objective of this invention can be achieved through the following technical solutions:

[0007] A real-time monitoring and early warning method for geological disasters that integrates air, space, and ground technologies includes the following steps:

[0008] S1. Real-time acquisition of monitoring data reported by sensor network nodes deployed in the monitoring area, the monitoring data including the displacement and tilt angle of each node position;

[0009] S2. Based on the monitoring data, identify sensor nodes whose displacement change rate or tilt change rate exceeds the primary threshold, and mark the nodes as primary warning nodes.

[0010] S3. When any primary warning node is marked, immediately dispatch flight devices to the airspace of that node, acquire surface images of the area surrounding the node, and analyze the surface images to identify the geometric features of surface cracks.

[0011] S4. Compare the displacement vector direction of the primary warning node with the extension direction of the nearest surface crack, and at the same time determine whether there are other primary warning nodes within a preset range around the node.

[0012] S5. If the directions are consistent and other primary warning nodes exist, then all related primary warning nodes are designated as a warning cluster, and the soil moisture content distribution data of the area covered by the warning cluster is obtained.

[0013] S6. Based on the displacement of nodes within the early warning cluster, the geometric characteristics of surface cracks, and the distribution data of soil moisture content, a geological disaster early warning command is generated and sent to the early warning execution device.

[0014] As a further aspect of the present invention: in S2, the process of marking the primary early warning node is as follows:

[0015] For each sensor node, determine whether the change in its displacement or tilt angle within a single sampling period exceeds the instantaneous threshold. If it does, mark the node as an instantaneous abnormal node and record the abnormality type and amplitude.

[0016] Centered on each instantaneous abnormal node, other sensor nodes are searched within a preset neighborhood radius. The monitoring data of these neighborhood nodes in the most recent time window are checked to determine whether their trend matches the abnormality type of the central abnormal node. When more than a certain number of nodes in the neighborhood show a matching abnormal trend, the initial central instantaneous abnormal node is marked as a primary warning node. If there are not enough nodes in the neighborhood to show a matching trend, the instantaneous abnormality mark of the central node is cleared.

[0017] As a further aspect of the present invention: in step S3, the process of scheduling the flight device to acquire images is as follows:

[0018] A scanning track containing multiple waypoints at different altitudes and inclination angles is generated. The flight device is equipped with a multispectral imager and a lidar and flies along the scanning track. The multispectral imager acquires visible light and near-infrared band images, and the lidar emits laser pulses and receives echoes to generate a three-dimensional point cloud. The images of different bands are weighted and fused by pixel values ​​to enhance the contrast of linear low reflectivity areas in the images.

[0019] The fused image is aligned with the 3D point cloud, so that each 3D point is associated with a spectral reflectance value. In the aligned data, a continuous spatial point sequence with a reflectance value lower than the reflectance threshold and an elevation difference between adjacent points higher than the elevation change threshold is searched. The continuous spatial point sequence that meets the conditions is connected and fitted into a 3D spatial polyline segment, which is the geometric center line of the surface crack. The distance between the two boundary points is calculated along the normal section of the geometric center line to obtain the width of the crack at each position. The total length and overall direction of the geometric center line are then calculated.

[0020] As a further aspect of the present invention: the process of connecting and fitting a sequence of continuous spatial points that meet the conditions into a three-dimensional spatial polyline segment is as follows:

[0021] For each spatial point that satisfies the conditions of reflectivity and elevation change, its three-dimensional coordinates are recorded. Starting from any point that meets the conditions, the next point that meets the conditions is searched within its neighboring spatial range, which is defined by a three-dimensional search radius. When the next point is found, a connection is established between the two points. The above neighbor search and connection process is repeated to form an initial path composed of multiple spatial line segments. The initial path is smoothed by fitting the endpoints of all line segments with a piecewise cubic polynomial curve to generate a smooth three-dimensional spatial curve. The three-dimensional spatial curve is discretized into a series of equally spaced polyline segments to obtain the final three-dimensional spatial polyline segment.

[0022] As a further aspect of the present invention: in step S4, the process of performing consistency comparison and association judgment is as follows:

[0023] For each primary early warning node, calculate the projection direction of its three-dimensional displacement vector on the horizontal plane. For each surface crack identified around the node, calculate the average projection direction of its geometric center line on the horizontal plane, and calculate the angle between the displacement projection direction and the average projection direction of the crack.

[0024] When the included angle is less than the first angle threshold, it is determined that the node is consistent with the direction of the crack. For all nodes and cracks with consistent directions, the association between the node and the crack is established. Based on the geographical location of all primary warning nodes, the Euclidean distance between any two nodes is calculated. When the distance between two nodes is less than the distance threshold, and the included angle between the average projection directions of the cracks associated with the two nodes on the horizontal plane is less than the second angle threshold, it is determined that there is a strong spatial association between the two nodes.

[0025] As a further aspect of the present invention: the process of calculating the angle between the average projection directions of the cracks associated with each of the two nodes on the horizontal plane is as follows:

[0026] For two primary early warning nodes that need to be compared, all surface cracks associated with each node are obtained. For each associated crack, the projection line segment of its geometric center line on the horizontal plane is extracted, and the direction of the line connecting the start and end points of each projection line segment is calculated as the projection direction vector of the crack. For each node, the projection direction vectors of all associated cracks of the node are vector-added to obtain a composite vector. The azimuth angles of the composite vectors corresponding to the two nodes are calculated, and the absolute value of the difference between the two azimuth angles is calculated to obtain the angle between the average projection directions of the cracks associated with the two nodes on the horizontal plane.

[0027] As a further aspect of the present invention: in step S5, the process of defining the early warning cluster is as follows:

[0028] All primary early warning nodes with strong spatial correlations are grouped into the same set of correlated nodes. The geographical center of the set of correlated nodes is calculated, and the maximum distance from all nodes in the set to the center is calculated. A circular area is delineated as the physical coverage of the early warning cluster, with the geographical center as the center and a fixed multiple of the maximum distance as the radius.

[0029] Within the physical coverage area, a brightness temperature data raster map generated by the microwave radiometer of the flight device is acquired. Based on the pre-stored mapping table between brightness temperature values ​​and soil volumetric moisture content, each pixel value in the brightness temperature data raster map is converted into a soil volumetric moisture content value, thereby generating a two-dimensional soil moisture content distribution map within the physical coverage area.

[0030] As a further aspect of the present invention: in step S6, the process of generating a geological disaster early warning command is as follows:

[0031] Obtain a two-dimensional soil moisture content distribution map within the physical coverage area of ​​the early warning cluster, calculate its average soil moisture content, obtain the latest displacement rate of all primary early warning nodes within the early warning cluster, calculate its average displacement rate, and obtain the latest average width expansion rate of all relevant surface cracks within the early warning cluster.

[0032] An early warning decision table is established, wherein the row index of the early warning decision table is a discrete range of average soil moisture content, the column index is a discrete range of average displacement rate, and the table cell value is the early warning level corresponding to the range of average width expansion rate. Based on the calculated average soil moisture content, average displacement rate, and average width expansion rate, the early warning decision table is queried to determine the corresponding early warning level, and a geological disaster early warning instruction containing the coordinates of the physical coverage area of ​​the early warning cluster, the early warning level, and the three average parameter values ​​is generated.

[0033] The beneficial effects of this invention are:

[0034] This invention achieves real-time on-site verification of abnormal areas by automatically scheduling flight devices for rapid response and image acquisition after real-time identification of primary early warning nodes in a ground sensor network. This effectively overcomes the limitations of space-air monitoring cycles and ground-based monitoring space, significantly shortening the time lag from anomaly detection to on-site verification. By fusing multispectral and lidar data to accurately identify the geometric features of surface cracks and comparing the consistency between node displacement direction and crack extension direction, combined with spatial correlation analysis between nodes, it can effectively distinguish between local interference and regional instability precursors, thereby significantly reducing the false alarm rate. Furthermore, by acquiring soil moisture distribution data within the early warning cluster and integrating it with multi-dimensional parameters such as displacement and cracks for comprehensive decision-making, a dynamic early warning process with real-time correlation analysis of multi-source data is constructed. Finally, based on an early warning decision-making mechanism that integrates multiple disaster-causing factors, early warning instructions are generated, comprehensively improving the accuracy, reliability, and timeliness of early warning information, and providing reliable technical support for the precise prevention and control of geological disasters. Attached Figure Description

[0035] The invention will now be further described with reference to the accompanying drawings.

[0036] Figure 1 This is a flowchart illustrating the present invention. Detailed Implementation

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

[0038] Please see Figure 1 As shown, this invention is a real-time monitoring and early warning method for geological disasters that integrates air, space, and ground technologies, and includes the following steps:

[0039] S1. Real-time acquisition of monitoring data reported by sensor network nodes deployed in the monitoring area, the monitoring data including the displacement and tilt angle of each node position;

[0040] S2. Based on the monitoring data, identify sensor nodes whose displacement change rate or tilt change rate exceeds the primary threshold, and mark the nodes as primary warning nodes.

[0041] S3. When any primary warning node is marked, immediately dispatch flight devices to the airspace of that node, acquire surface images of the area surrounding the node, and analyze the surface images to identify the geometric features of surface cracks.

[0042] S4. Compare the displacement vector direction of the primary warning node with the extension direction of the nearest surface crack, and at the same time determine whether there are other primary warning nodes within a preset range around the node.

[0043] S5. If the directions are consistent and other primary warning nodes exist, then all related primary warning nodes are designated as a warning cluster, and the soil moisture content distribution data of the area covered by the warning cluster is obtained.

[0044] S6. Based on the displacement of nodes within the early warning cluster, the geometric characteristics of surface cracks, and the distribution data of soil moisture content, a geological disaster early warning command is generated and sent to the early warning execution device.

[0045] In a preferred embodiment of the present invention, the process of marking the primary early warning node in step S2 is as follows:

[0046] It continuously receives real-time monitoring data streams from a sensor network deployed throughout the structure. Each sensor node uploads a data packet every 30 seconds, containing the displacement in millimeters and the tilt angle measured by that node within the current sampling period.

[0047] For each sensor node, the processing logic first calculates the instantaneous change of its key monitoring indicators between two adjacent sampling periods. Specifically, the program reads the displacement data D of that node at the current time T. t And read the displacement data D from the previous time step T-1. t-1 Calculate the displacement difference ΔD=D t -D t-1 Simultaneously, read the tilt angle value A between the current and previous moments. t and A t-1Calculate the difference in angle change ΔA = A t -A t-1 .

[0048] The program then compares the calculated ΔD and ΔA with preset instantaneous thresholds. These instantaneous thresholds are set based on the structure's safety specifications and historical baseline data. For example, the instantaneous threshold for displacement change might be set to 2.0 mm, and the instantaneous threshold for tilt angle change might be set to 0.5 degrees.

[0049] The judgment logic is as follows: If the absolute value of ΔD is greater than 2.0 mm, or the absolute value of ΔA is greater than 0.5 degrees, then the node is determined to have experienced a transient anomaly within the current sampling period. The program immediately marks the node as a transient anomaly node. The marking operation includes setting an anomaly flag in the node's status attributes and creating an anomaly record entry. This entry records the precise timestamp T of the anomaly occurrence, the anomaly type is determined based on the threshold exceeded (e.g., displacement exceeding limit or tilt exceeding limit), and the specific anomaly amplitude, i.e., the value of ΔD or ΔA. For example, if node ID45 calculates ΔA to be 0.8 degrees at time 14:25:30, exceeding the 0.5-degree threshold, it is marked as a transient anomaly node, the anomaly type is recorded as tilt exceeding limit, and the amplitude is 0.8 degrees.

[0050] After detecting and marking transient anomalies, the program proceeds to the second stage: neighborhood correlation verification. The purpose of this stage is to eliminate isolated anomalies caused by sensor malfunctions or transient interference, and to identify potential risk points with spatial correlation.

[0051] The program iterates through all nodes marked as transient anomalies. For each such central node, it uses its coordinates in the structural space as the center and a preset neighborhood radius, such as 5 meters, as the range. It then searches the sensor network topology to find all other sensor nodes located within this circular area; these nodes are called neighboring nodes. For example, using the transient anomaly node ID45 as the center, the search finds eight neighboring nodes within a 5-meter radius, with IDs 46, 47, 48, 49, 50, 51, 52, and 53.

[0052] Next, the program examines the monitoring data of these eight neighboring nodes within the most recent time window. The length of the time window typically covers multiple sampling periods, such as the most recent 10 minutes, corresponding to 20 data points. The goal of the examination is to determine whether the data change trends of these neighboring nodes match the anomaly type of the central anomalous node.

[0053] The rules for determining trend consistency depend on the anomaly type. If the anomaly type of the central node ID45 is tilt exceeding the limit with an amplitude of positive 0.8 degrees (i.e., tilt increasing), the program analyzes the tilt time series data of each neighboring node within the most recent 10 minutes. It calculates the linear regression slope of the tilt value of each neighboring node within that time window, or simply calculates the change in its first and last values. If the tilt of a neighboring node also shows a significant increasing trend within the window period, for example, its cumulative tilt increase exceeding 0.2 degrees, and the direction of change is the same as that of the central node, then that neighboring node is determined to exhibit an anomaly trend consistent with that of the central node. For displacement anomalies, the program determines whether the direction and magnitude of the displacement change have spatial consistency.

[0054] The program counts the number of nodes exhibiting an abnormal trend across all neighboring nodes. This number needs to be compared with a preset threshold. The threshold can be set as a percentage of the total number of neighboring nodes, such as 60%. For example, suppose that out of 8 neighboring nodes, 6 nodes show a tilting increase trend consistent with ID45 within the last 10 minutes.

[0055] Since 6 nodes constitute more than 60% of 8 (4.8 nodes), meeting the quantity threshold, the program determines that the anomaly has spatial correlation. The program upgrades the status of the initial central transient anomaly node ID45, marking it as a primary warning node. The upgrade operation includes setting a higher-level warning flag in the node's status attributes and may associate its anomaly record with the list of neighboring node IDs participating in the verification.

[0056] Conversely, if only two out of eight neighboring nodes exhibit a similar trend, failing to reach the quantity threshold, the program determines that this transient anomaly may be an isolated event. The program will clear the transient anomaly marker for node ID45, restore its state to normal, and archive this anomaly record for subsequent analysis, but will not escalate it to an alert.

[0057] In another preferred embodiment of the present invention, in step S3, the process of scheduling the flight device to acquire images is as follows:

[0058] First, flight mission planning is required. Based on the terrain boundary coordinates of the area to be monitored and the preset imaging resolution requirements, the program automatically generates a scanning track containing multiple waypoints at different altitudes and tilt angles. Track planning ensures full coverage and overlapping data acquisition of the target area. For example, for a rectangular area approximately 500 meters long and 300 meters wide, the program might generate a zigzag track containing 50 waypoints, with flight altitudes alternating between 80 and 150 meters, and camera tilt angles switching between vertically downward and tilted at 30 degrees to achieve multi-angle observation.

[0059] The flight vehicle, equipped with a multispectral imager and lidar, autonomously flies along a planned scanning path. The multispectral imager synchronously acquires digital images in the visible and near-infrared bands at a fixed frequency. The visible band typically includes three channels: red, green, and blue, while the near-infrared band is a separate channel. The lidar device scans at a rate of hundreds of thousands of laser pulses per second, precisely recording the time difference between the emission and return of each pulse. Combined with the flight vehicle's positioning and attitude data, the three-dimensional spatial coordinates of each laser footprint are calculated, thereby generating high-density three-dimensional point cloud data. After a flight mission is completed, the obtained dataset includes a set of multispectral images with geographic location information and a set of three-dimensional point clouds covering the same area.

[0060] Next, the multispectral imagery is processed. The program reads visible and near-infrared images taken at the same location, first performing radiometric and geometric corrections to eliminate sensor errors and terrain distortion. Then, pixel-level weighted fusion is performed on the corrected images of different bands. The weighting coefficients are set according to the target enhancement requirements; for example, higher positive weights are assigned to the near-infrared band, and lower negative weights to the red band. This weighting operation highlights the differences between vegetation and soil, and particularly enhances the contrast between linear low-reflectivity areas in the image, such as potential cracks, and the surrounding background. The fused image is then a composite image with enhanced features.

[0061] Next, coordinate alignment, i.e., registration between the point cloud and the image, is performed. The program utilizes high-precision timestamps recorded by the flight device and synchronized geographic location and attitude data to project the 3D coordinates of each LiDAR point cloud onto the fused multispectral image plane. By solving the projection transformation matrix, the corresponding pixel position on the 2D image is found for each 3D point cloud point, and the fused spectral reflectance value is read from that pixel position. This establishes an enhanced point cloud dataset where each 3D point contains both spatial coordinates and spectral reflectance attributes.

[0062] Crack feature extraction is performed on this augmented point cloud data. The program first searches for potential crack points throughout the entire dataset. The search criteria are based on reflectance values ​​and local terrain features. Cracks typically exhibit differences in opening and filling material, have low reflectance, and may have a small elevation difference between their two sides. The program sets a reflectance threshold, for example, 0.15, and an elevation change threshold, for example, 0.05 meters. The algorithm iterates through each point, checking if its reflectance value is below 0.15, and simultaneously calculates the absolute value of the elevation difference between that point and its nearest neighbor. If the difference is greater than 0.05 meters, the point is marked as a candidate crack point.

[0063] After the search is complete, a series of discrete candidate points satisfying the initial conditions are obtained. The program needs to connect these points to form a continuous crack line. The specific process is as follows: For each marked candidate point, the program records its precise three-dimensional coordinates. The connection process starts from an arbitrary candidate point, which is used as the starting point of the current path. The program defines a three-dimensional search radius, for example, 0.2 meters, and searches for other candidate points within this spherical space centered on the current point. If other candidate points exist within the range, the program calculates the direction vector between the current point and each point within the range, and selects the point whose direction is most consistent with the extension direction of the current path or is the closest to it as the next connection point. A virtual connecting line segment is established between the two points.

[0064] After the connection is completed, the newly added point becomes the current point. The program repeats the above neighbor search and connection process to find the next point that meets the conditions. This process continues, like stringing beads, gradually forming an initial 3D path composed of multiple short line segments. To prevent the path from extending indefinitely or forming loops, the program sets stopping conditions, such as stopping the growth of the current path when there are no new candidate points within the search range or when the path length exceeds a preset maximum value.

[0065] Once a path has grown, the program selects a new starting point from the remaining unconnected candidate points and repeats the process until all candidate points have been attempted to be connected. This may ultimately generate multiple independent initial paths, each representing a potential crack segment.

[0066] Since the paths obtained by direct connection may be jagged, the program smooths each initial path. The smoothing process uses a piecewise cubic polynomial curve fitting algorithm, such as cubic spline interpolation. The program uses the endpoints of all line segments on the path as control points to fit a smooth, continuous three-dimensional curve. This curve more realistically reflects the natural meandering shape of the crack.

[0067] To facilitate subsequent geometric measurements and analysis, the program re-discretizes the fitted smooth 3D spatial curve. It samples a series of new points along the curve at fixed step sizes, such as 0.1 meters, and connects these new points sequentially with straight line segments. The result is a 3D spatial polygonal line segment composed of a series of equally spaced short polygonal line segments; this polygonal line segment represents the geometric centerline of the identified surface crack.

[0068] After obtaining the geometric centerline, the program calculates the crack width along this centerline. For each sampling point on the centerline, the program calculates the normal plane of the centerline at that point. Within the normal plane, it searches for boundary points in the augmented point cloud from that point outwards. The boundary point is determined by the transition point where the reflectivity changes from below a threshold to above a threshold. The program scans along the normal direction on both sides within the normal plane to find the nearest boundary point and calculates the three-dimensional Euclidean distance between these two points; this distance represents the crack width at that location. By traversing all sampling points on the centerline, the crack width distribution along its length can be obtained.

[0069] Finally, the program performs overall geometric statistics on the cracks. The total length of the geometric centerline is obtained by summing the lengths of all polyline segments. The overall orientation is determined by calculating the principal direction of the projection of the centerline onto the horizontal plane; for example, the orientation angle can be calculated by performing linear regression on the projection line and determining its angle with true north. The program outputs the three-dimensional coordinate sequence of the geometric centerline of each crack, the width at each point, the total length, and the overall orientation, completing the entire extraction process from remote sensing data to quantitative crack parameters.

[0070] In another preferred embodiment of the present invention, the process of consistency comparison and association judgment in step S4 is as follows:

[0071] The input data for this module comes from two upstream processing results: one is all the primary early warning nodes and their attribute information output by S2; the other is all the identified surface cracks and their geometric parameters output by S3.

[0072] First, the module performs directional consistency analysis on each primary early warning node. Each primary early warning node contains a key attribute: a three-dimensional displacement vector. This vector describes the spatial displacement change of the node over a monitoring period and typically consists of three components: east-west displacement DX, north-south displacement DY, and vertical displacement DZ. For example, a node might record displacement vector components where DX equals 0.05 meters, DY equals 0.12 meters, and DZ equals -0.02 meters.

[0073] Calculate the projection direction of the 3D displacement vector onto the horizontal plane. The program ignores the vertical component DZ and uses only the horizontal components DX and DY. The direction angle of the projected vector is obtained through mathematical calculation. The direction angle is usually expressed as a clockwise rotation with true north as 0 degrees. Specifically, the arctangent of the angle between the projected vector and the true north axis is calculated first. The formula involves the ratio of DX to DY. For example, for a DX of 0.05 meters and a DY of 0.12 meters, the calculated direction angle is approximately 157.38 degrees. This angle represents the main displacement trend direction of the node on the horizontal plane.

[0074] Simultaneously, the program needs to process all surface cracks identified in the area surrounding the nodes. For each crack, S3 has calculated and output the average projection direction of its geometric centerline on the horizontal plane. This direction is also an azimuth angle based on true north. For example, the average projection azimuth angle of a crack F1 is 160 degrees.

[0075] Next, the program calculates the absolute difference between the displacement projection direction angle of the primary warning node and the average projection direction angle of each crack. This difference is the angle between the two directions. For example, the difference between the node direction angle of 157.38 degrees and the crack F1 direction angle of 160 degrees is 2.62 degrees.

[0076] The program presets a first angle threshold to determine if the directions are consistent. This threshold is usually set to 30 degrees. The judgment logic is as follows: if the calculated angle is less than 30 degrees, the primary warning node is determined to be consistent with the horizontal movement trend of the surface crack. For example, 2.62 degrees is less than 30 degrees, therefore the node is determined to be consistent with the direction of crack F1.

[0077] For a primary warning node, there may be multiple cracks around it. The program will compare its direction with all cracks as described above. All cracks that pass the consistency check (i.e., the included angle is less than 30 degrees) will be recorded as associated cracks of that node. The program will create a list of associated cracks for each node. For example, node N45 may be aligned with cracks F1 and F3, so its list of associated cracks will include F1 and F3. At this point, the module has established a preliminary association between structural anomaly nodes and surface cracks.

[0078] After completing the independent association of all nodes, the module proceeds to the second stage of spatial association analysis between nodes. The goal is to identify groups of nodes that are spatially close and have similar deformation patterns, which may indicate a larger, coherent unstable region.

[0079] The program first calculates the three-dimensional Euclidean distance between any two nodes based on their geographical coordinates. The node coordinates are determined by measurements taken during deployment. For example, node N45 has coordinates of X100.5 meters, Y200.3 meters, and Z50.1 meters; node N67 has coordinates of X95.2 meters, Y205.8 meters, and Z50.5 meters. The calculated straight-line distance between them is approximately 7.2 meters.

[0080] The program presets a distance threshold, such as 10 meters. The judgment logic is: if the distance between two nodes is less than 10 meters, they are considered to be spatially adjacent. For example, the distance between N45 and N67 is 7.2 meters, which is less than 10 meters, thus satisfying the proximity condition.

[0081] Spatial proximity alone is insufficient to determine a strong correlation; it is also necessary to examine whether the deformation trends reflected by the two nodes are consistent. This is achieved by comparing the overall average direction of their respective associated crack groups. The calculation process is as follows.

[0082] For two nodes that need to be compared, such as N45 and N67, the program retrieves their respective lists of associated cracks. Assume N45 is associated with cracks F1 and F3, and N67 is associated with cracks F2 and F3.

[0083] For each associated crack, the program extracts the projection segment of its geometric centerline onto the horizontal plane. This projection segment consists of a series of two-dimensional coordinate points. The program calculates the direction of the line connecting the start and end points of this projection segment, which serves as a vector representing the projection direction of the crack. For example, the projection segment of crack F1 has a start point coordinate of X1, Y1 and an end point coordinate of X2, Y2. The calculated direction vector angle is 160 degrees.

[0084] For each node, the program performs vector addition on the projected direction vectors of all its associated cracks. Vector addition considers not only direction but also the length or weight of each vector. A simple implementation is to assign the same unit length to each crack vector. For example, node N45 has two associated cracks with direction angles of 160 degrees and 165 degrees respectively. Converting these two direction angles to unit vectors and combining them might yield a composite vector with an direction angle of approximately 162.5 degrees. This composite vector's direction angle represents the average projected direction of all associated cracks at that node on the horizontal plane.

[0085] After calculating the azimuth angles of the composite vectors corresponding to the two nodes, calculate the absolute value of the difference between these two azimuth angles. This absolute value is the angle between the average projection directions of the cracks associated with each node on the horizontal plane. For example, the composite direction angle of node N45 is 162.5 degrees, and the composite direction angle of node N67 is 158.0 degrees, with an absolute difference of 4.5 degrees.

[0086] The program presets a second angle threshold, which is usually more stringent than the first angle threshold, for example, set to 20 degrees. The judgment logic is: if the angle between the average directions of the associated crack groups of two adjacent nodes is less than 20 degrees, then their deformation trends are considered to be consistent.

[0087] Ultimately, determining a strong spatial correlation requires both conditions to be met simultaneously: first, the Euclidean distance between the two primary warning nodes is less than a distance threshold of 10 meters; second, the angle between the average projection directions of the crack groups associated with each of the two nodes is less than a second angle threshold of 20 degrees. Only when both conditions are met can the program determine that a strong spatial correlation exists between the two nodes.

[0088] For example, the distance between nodes N45 and N67 is 7.2 meters, which is less than 10 meters, and the average directional angle between the associated crack groups of the two is 4.5 degrees, which is less than 20 degrees. Therefore, it is determined that there is a strong spatial correlation between N45 and N67.

[0089] The module outputs a relational network containing a list of cracks associated with each primary warning node, as well as a list of all node pairs identified as having strong spatial correlations. This network provides crucial spatial and causal evidence for subsequent comprehensive risk assessment, helping to identify potential patterns that may represent overall structural instability from discrete anomalies.

[0090] In another preferred embodiment of the present invention, the process of defining the early warning cluster in step S5 is as follows:

[0091] The process of defining early warning clusters begins with processing the list of strongly spatially correlated node pairs output by S4. The module reads this list, where each record contains identifiers for two strongly correlated primary early warning nodes. The program first needs to convert these pairwise relationships into one or more interconnected sets of nodes.

[0092] The program initializes an empty list of sets to store different sets of associated nodes. It takes the first pair of nodes from the strongly associated list, such as nodes A and B, creates a new set of associated nodes, and adds nodes A and B to this set. The program then iterates through the remaining pairs of associated nodes. For each new pair of nodes, such as nodes C and D, the program checks if these two nodes already exist in any of the existing sets of associated nodes. If node C already exists in set 1, and node D does not exist in any set, then node D is also added to set 1. If node C is in set 1 and node D is in set 2, then set 1 and set 2 are merged into a new set containing all nodes. If neither node appears in any existing set, a new set is created using them. Through this iterative merge algorithm, all nodes connected by direct or indirect strong associations are eventually aggregated into the same set of associated nodes, while nodes that are not associated with each other belong to different sets. For example, this might result in two sets: set S1 containing nodes A, B, and C; and set S2 containing nodes D and E.

[0093] For each set of associated nodes, the program needs to define the physical coverage area of ​​its corresponding early warning cluster. First, the geographical center of the set is calculated. The program reads the three-dimensional geographic coordinates of each node in the set and calculates the average of the X-coordinate, the average of the Y-coordinate, and the average of the Z-coordinate of all nodes. These three averages constitute the coordinates of the geographical center point O of the set. For example, if set S1 contains three nodes with X-coordinates of 100.0 meters, 102.0 meters, and 98.0 meters, then the X-coordinate of the center point O is 100.0 meters.

[0094] Next, the Euclidean distance from all nodes in the set to the center point O is calculated. For each node, the program calculates the difference between its coordinates and the coordinates of the center point O, sums the squares, and then takes the square root to obtain the distance from the node to the center point. After traversing all nodes in the set, the maximum distance value Dmax is found. For example, after calculation, the maximum distance from the node in set S1 to the center O is 15.3 meters.

[0095] The program then defines a circular area centered at the geographic center point O, with a radius equal to Dmax multiplied by a fixed expansion factor K. The expansion factor K ensures the warning range covers the potential impact area of ​​the node cluster and is typically set empirically, for example, K equals 1.2. The warning radius R is calculated as Dmax * 1.2. In the example above, R = 15.3 * 1.2 = 18.36 meters. Therefore, a circular area centered at O ​​with a radius of 18.36 meters is determined as the physical coverage area of ​​the warning cluster. The program records the coordinates of the circle's center and its radius.

[0096] After delineating the physical coverage area, the program needs to acquire soil moisture information within that area. This is achieved by processing scanning data from the microwave radiometer onboard the flight device. During flight operations, the microwave radiometer generates a brightness temperature raster map covering the entire monitoring area. Each raster pixel value represents the radiation brightness temperature value of a ground resolution unit in a specific microwave frequency band, measured in Kelvin.

[0097] The program crops the corresponding sub-region data from the complete brightness temperature raster map based on the circular boundary coordinates of the warning cluster. Next, it converts the brightness temperature values ​​to soil volumetric moisture content. The module pre-stores a mapping table, built using historical experimental data, to correlate brightness temperature values ​​with soil volumetric moisture content. For example, the table records a brightness temperature value of 270K as potentially corresponding to a moisture content of 15%, and a brightness temperature value of 280K as potentially corresponding to a moisture content of 10%. The conversion process iterates through each pixel in the cropped brightness temperature raster map, looks up its brightness temperature value in the mapping table, and calculates the corresponding soil volumetric moisture content value as a percentage through linear interpolation. For example, a pixel with a brightness temperature of 275K, after looking up the mapping table and interpolating, yields a moisture content of approximately 12.5%.

[0098] After converting all pixels, the program generates a new raster map with the same resolution and range as the original brightness temperature map, but the value of each pixel is the soil volumetric moisture content. This map is the two-dimensional soil moisture content distribution map within the physical coverage area of ​​the warning cluster. The program binds and stores this distribution map with the boundary information of the warning cluster for use by the next module.

[0099] In another preferred embodiment of the present invention, the process of generating a geological disaster early warning command in step S6 is as follows:

[0100] For each warning cluster to be processed, the module first extracts data from its associated two-dimensional soil moisture content distribution map. The program reads the soil moisture content values ​​of all cells in the distribution map, calculates the arithmetic mean of these values, and obtains the average soil moisture content value θ within the range of that warning cluster. avg For example, θ is calculated. avg It is 18.3%.

[0101] Secondly, the module obtains the latest displacement rates of all primary early warning nodes within the early warning cluster. Displacement rate refers to the rate of change of a node's displacement within the most recent calculation period. The program reads the latest reported or calculated three-dimensional displacement rate vector from the real-time database or node status cache for each node in the cluster and calculates the magnitude of this vector, i.e., the total displacement rate. For example, if an early warning cluster contains three nodes with latest displacement rates of 2.1 mm / day, 1.8 mm / day, and 2.5 mm / day, respectively, the program calculates the arithmetic mean of these three values ​​to obtain the average displacement rate V of the early warning cluster. avg For example, 2.13 millimeters per day.

[0102] Third, the module obtains the latest average width expansion rate of all relevant surface cracks within the warning cluster. The width expansion rate of each crack is obtained through time-series analysis of crack width data identified at different times. The program reads all crack identifiers associated with the warning cluster and retrieves the average width expansion rate value of each crack in the most recent monitoring period from the crack status database. For example, the expansion rates of three associated cracks are 0.5 mm / day, 0.7 mm / day, and 0.6 mm / day, respectively. The program calculates the arithmetic mean of these values ​​to obtain the average crack width expansion rate W of the warning cluster. avg For example, 0.6 millimeters per day.

[0103] After obtaining three key average parameters, the module determines the warning level by querying the warning decision table. The warning decision table is a predefined three-dimensional lookup table with the following structure: The row index is a discretized interval of average soil moisture content, for example, interval 1 is less than 15%, interval 2 is 15% to 25%, and interval 3 is greater than 25%. The column index is a discretized interval of average displacement rate, for example, interval A is less than 1 mm / day, interval B is 1 to 3 mm / day, and interval C is greater than 3 mm / day. Each table cell at the intersection of row and column is not a single value, but a sub-rule whose content is associated with the average crack width expansion rate W. avg The interval. For example, for the intersection of row index interval 2 and column index interval B, the rule might be: if W avg If the daily rainfall (Wavg) is less than 0.5 mm, the warning level is blue; if W... avgIf the rainfall is between 0.5 and 1.0 mm per day, the warning level is yellow; if W avg If the rainfall exceeds 1.0 mm per day, the warning level is orange.

[0104] The module's query process consists of three steps. The first step involves calculating the average soil moisture content θ. avg This determines the row index range to which it belongs. For example, θ avg The value is 18.3%, falling within interval 2. The second step is to determine the average displacement rate V. avg This determines the column index range to which it belongs. For example, V avg The rate is 2.13 mm per day, belonging to interval B. The third step is to locate the table cell where row interval 2 intersects column interval B, and based on the rules defined within that cell, combined with the average crack width propagation rate W... avg The value of W determines the final warning level. For example, W avg The value is 0.6 mm per day. According to the example rule above, its value falls between 0.5 and 1.0 mm per day, so the warning level is determined to be yellow.

[0105] Finally, the module generates a complete geological disaster early warning instruction. This instruction is a structured data message containing the following core fields: a unique identifier for the early warning cluster, the center coordinates and radius of the physical coverage area of ​​the cluster, the determined early warning level code, and three average parameter values ​​used as the basis for decision-making. For example, an instruction might contain the following: Early warning cluster CL001, center point coordinates X100.0, Y200.0, radius 18.36 meters, early warning level yellow, average moisture content 18.3%, average displacement rate 2.13 mm / day, average crack propagation rate 0.6 mm / day. This instruction is immediately sent to the early warning information dissemination and emergency response platform, thus completing a fully automated process from multi-source data fusion analysis to quantitative early warning decision-making.

[0106] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A real-time monitoring and early warning method for geological disasters integrating air-space-ground coordination, characterized in that, Includes the following steps: S1. Real-time acquisition of monitoring data reported by sensor network nodes deployed in the monitoring area, the monitoring data including the displacement and tilt angle of each node position; S2. Based on the monitoring data, identify sensor nodes whose displacement change rate or tilt change rate exceeds the primary threshold, and mark the nodes as primary warning nodes. S3. When any primary warning node is marked, immediately dispatch flight devices to arrive at the airspace of that node, acquire surface images of the preset area around the node, and analyze the surface images to identify the geometric features of surface cracks. S4. Compare the displacement vector direction of the primary warning node with the extension direction of the nearest surface crack, and at the same time determine whether there are other primary warning nodes within a preset range around the node. S5. If the directions are consistent and other primary warning nodes exist, then all related primary warning nodes are designated as a warning cluster, and the soil moisture content distribution data of the area covered by the warning cluster is obtained. S6. Based on the displacement of nodes in the early warning cluster, the geometric characteristics of surface cracks, and the distribution data of soil moisture content, a geological disaster early warning command is generated and the early warning command is sent to the early warning execution device. In step S4, the process of consistency comparison and association determination is as follows: For each primary warning node, calculate the projection direction of its three-dimensional displacement vector on the horizontal plane. For each surface crack identified around the node, calculate the average projection direction of its geometric center line on the horizontal plane, and calculate the angle between the displacement projection direction and the average projection direction of the crack. When the included angle is less than the first angle threshold, it is determined that the node is consistent with the direction of the crack. For all nodes and cracks with consistent directions, the association between the node and the crack is established. Based on the geographical location of all primary warning nodes, the Euclidean distance between any two nodes is calculated. When the distance between two nodes is less than the distance threshold, and the included angle between the average projection directions of the cracks associated with the two nodes on the horizontal plane is less than the second angle threshold, it is determined that there is a strong spatial association between the two nodes. The process of calculating the angle between the average projection directions of the cracks associated with each of the two nodes on the horizontal plane is as follows: For two primary early warning nodes that need to be compared, all surface cracks associated with each node are obtained. For each associated crack, the projection line segment of its geometric center line on the horizontal plane is extracted, and the direction of the line connecting the start and end points of each projection line segment is calculated as the projection direction vector of the crack. For each node, the projection direction vectors of all associated cracks of the node are converted into unit vectors and then vector superposition is performed to obtain a composite vector. The azimuth angles of the composite vectors corresponding to the two nodes are calculated, and the absolute value of the difference between the two azimuth angles is calculated to obtain the angle between the average projection directions of the cracks associated with the two nodes on the horizontal plane.

2. The integrated air-space-ground geological disaster real-time monitoring and early warning method according to claim 1, characterized in that, In S2, the process of marking the primary early warning node is as follows: For each sensor node, determine whether the change in its displacement or tilt angle within a single sampling period exceeds the instantaneous threshold. If it does, mark the node as an instantaneous abnormal node and record the abnormality type and amplitude. Centered on each instantaneous abnormal node, other sensor nodes are searched within a preset neighborhood radius. The monitoring data of these neighborhood nodes in the most recent time window are checked to determine whether their trend matches the abnormality type of the central instantaneous abnormal node. When more than a certain number of nodes in the neighborhood show a matching abnormal trend, the initial central instantaneous abnormal node is marked as a primary warning node. If there are not enough nodes in the neighborhood to show a matching trend, the instantaneous abnormality mark of the central instantaneous abnormal node is cleared.

3. The integrated air-space-ground geological disaster real-time monitoring and early warning method according to claim 1, characterized in that, In S3, the process of scheduling the flight device to acquire images is as follows: A scanning track containing multiple waypoints at different altitudes and inclination angles is generated. The flight device is equipped with a multispectral imager and a lidar and flies along the scanning track. The multispectral imager acquires visible light and near-infrared band images, and the lidar emits laser pulses and receives echoes to generate a three-dimensional point cloud. The images of different bands are weighted and fused by pixel values ​​to enhance the contrast of linear low reflectivity areas in the images. The fused image is aligned with the 3D point cloud, so that each 3D point is associated with a spectral reflectance value. In the aligned data, a continuous spatial point sequence with a reflectance value lower than the reflectance threshold and an elevation difference between adjacent points higher than the elevation change threshold is searched. The continuous spatial point sequence that meets the conditions is connected and fitted into a 3D spatial polyline segment, which is the geometric center line of the surface crack. The distance between the two boundary points is calculated along the normal section of the geometric center line to obtain the width of the crack at each position. The total length and overall direction of the geometric center line are then calculated.

4. The integrated air-space-ground geological disaster real-time monitoring and early warning method according to claim 3, characterized in that, The process of connecting and fitting a continuous sequence of spatial points that meet the conditions into a three-dimensional spatial polyline segment is as follows: For each spatial point that satisfies the conditions of reflectivity and elevation change, its three-dimensional coordinates are recorded. Starting from any point that meets the conditions, the next point that meets the conditions is searched within its neighboring spatial range, which is defined by a three-dimensional search radius. When the next point is found, a connection is established between the two points. The neighbor search and connection process is repeated to form an initial path composed of multiple spatial line segments. The initial path is smoothed by fitting the endpoints of all line segments with a piecewise cubic polynomial curve to generate a smooth three-dimensional spatial curve. The three-dimensional spatial curve is discretized into a series of equally spaced polyline segments to obtain the final three-dimensional spatial polyline segment.

5. The integrated air-space-ground geological disaster real-time monitoring and early warning method according to claim 3, characterized in that, In S5, the process of defining the early warning cluster is as follows: All primary early warning nodes with strong spatial correlations are grouped into the same set of correlated nodes. The geographical center of the set of correlated nodes is calculated, and the maximum distance from all nodes in the set to the center is calculated. A circular area is delineated as the physical coverage of the early warning cluster, with the geographical center as the center and a fixed multiple of the maximum distance as the radius. Within the physical coverage area, a brightness temperature data raster map generated by the microwave radiometer of the flight device is acquired. Based on the pre-stored mapping table between brightness temperature values ​​and soil volumetric moisture content, each pixel value in the brightness temperature data raster map is converted into a soil volumetric moisture content value, thereby generating a two-dimensional soil moisture content distribution map within the physical coverage area.

6. The integrated air-space-ground geological disaster real-time monitoring and early warning method according to claim 5, characterized in that, In S6, the process of generating geological disaster early warning instructions is as follows: Obtain a two-dimensional soil moisture content distribution map within the physical coverage area of ​​the early warning cluster, calculate its average soil moisture content, obtain the latest displacement rate of all primary early warning nodes within the early warning cluster, calculate its average displacement rate, and obtain the latest average width expansion rate of all relevant surface cracks within the early warning cluster. An early warning decision table is established, wherein the row index of the early warning decision table is a discrete range of average soil moisture content, the column index is a discrete range of average displacement rate, and the table cell value is the early warning level corresponding to the range of average width expansion rate. Based on the calculated average soil moisture content, average displacement rate, and average width expansion rate, the early warning decision table is queried to determine the corresponding early warning level, and a geological disaster early warning instruction containing the coordinates of the physical coverage area of ​​the early warning cluster, the early warning level, and the three average parameter values ​​is generated.

Images