A geological disaster risk dynamic assessment method based on multi-source spatio-temporal big data

Abstract

This application relates to the field of geological disaster monitoring and early warning technology, and in particular to a dynamic geological disaster risk assessment method based on multi-source spatiotemporal big data. The method includes: acquiring multi-source spatiotemporal data of a target area and performing spatiotemporal registration; extracting characteristics such as deformation acceleration, rainfall accumulation, water saturation, and vegetation degradation; and calculating an initial geological disaster risk value. For high-risk areas exceeding a first risk threshold, spatial correlation analysis is performed using geological attribute information to obtain risk correction factors; the target geological disaster risk value is calculated; and early warning information is generated. This application solves the problems of insufficient multi-source data fusion and lack of dynamic updating capability in existing technologies by fusing multi-source time-series monitoring data with geological attribute information for dynamic risk assessment, thereby improving the accuracy and timeliness of geological disaster risk assessment.

Description

Technical Field

[0001] This application relates to the field of geological disaster monitoring and early warning technology, and in particular to a dynamic assessment method for geological disaster risk based on multi-source spatiotemporal big data. Background Technology

[0002] Geological disasters such as landslides, collapses, and debris flows are frequent natural disasters worldwide, seriously threatening people's lives and property and the stable operation of major engineering facilities. Existing geological disaster risk assessment methods mainly rely on single-source ground patrol data or human experience judgment, making it difficult to comprehensively analyze multi-dimensional time-series monitoring information such as surface deformation, rainfall, soil moisture, and vegetation cover. This results in insufficient comprehensiveness and accuracy of risk assessments. Furthermore, traditional methods lack real-time response mechanisms for dynamic risk changes after the disaster body enters the accelerated deformation stage. After the early warning information is issued, it cannot be dynamically updated based on the latest meteorological forecast data and monitoring data, resulting in a problem of delayed early warning. Summary of the Invention

[0003] To help solve the technical problems of insufficient fusion of multi-source monitoring data and lack of dynamic updating capability of early warning information in the existing technology, this application provides a dynamic assessment method for geological disaster risk based on multi-source spatiotemporal big data.

[0004] This application provides a dynamic assessment method for geological hazard risk based on multi-source spatiotemporal big data, which adopts the following technical solution: A dynamic assessment method for geological hazard risk based on multi-source spatiotemporal big data includes: Acquire multi-source spatiotemporal data of the target area, including time-series surface deformation data, time-series rainfall data, time-series soil moisture data, and time-series vegetation cover data; Spatiotemporal registration is performed on the time-series surface deformation data, the time-series rainfall data, the time-series soil moisture data, and the time-series vegetation cover data to obtain multiple spatiotemporal matching sequences; Based on each spatiotemporal matching sequence, features of deformation acceleration, rainfall accumulation, water saturation, and vegetation degradation are extracted respectively. Based on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics, the initial geological hazard risk value is calculated; If the initial geological hazard risk value exceeds the first risk threshold, the target area will be marked as a high-risk area; Obtain the geological attribute information corresponding to the high-risk area, including stratigraphic lithology information, geological structure information, and historical disaster point distribution information; Based on the geological attribute information, spatial correlation analysis is performed on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics in the high-risk area to obtain risk correction factors; Calculate the target geological hazard risk value based on the initial geological hazard risk value and the risk correction factor; If the target geological hazard risk value exceeds the second risk threshold, a geological hazard early warning message will be generated.

[0005] Optionally, the calculation of the initial geological hazard risk value based on the deformation acceleration characteristic, the rainfall accumulation characteristic, the water saturation characteristic, and the vegetation degradation characteristic includes: Obtain the sequence length corresponding to each spatiotemporal matching sequence; Based on the sequence length, the deformation acceleration feature is piecewise linearly fitted to obtain the deformation rate within each segment, and the difference in deformation rate between adjacent segments is obtained as the rate change value. If the rate change value exceeds the preset rate change threshold, it is determined that there is an accelerated deformation stage, and the ratio of the duration of the accelerated deformation stage to the total observation time is obtained as the acceleration percentage. Based on the rainfall accumulation characteristics, a sequence of consecutive rainless days is obtained, and it is determined whether the length of the sequence of consecutive rainless days exceeds a first length threshold. If the length of the consecutive rainless day sequence exceeds the first length threshold, then within a preset time window after rainfall recovery, the surface deformation response amplitude corresponding to a unit rainfall amount is extracted. Based on the water saturation characteristics, a soil moisture content time series is obtained, and it is determined whether the soil moisture content time series continuously exceeds a preset moisture content threshold and the duration exceeds a second length threshold. If the soil moisture content time series continuously exceeds the preset moisture content threshold and the duration exceeds the second length threshold, then the slope of the soil moisture content time series curve is obtained as the saturation slope. Based on the vegetation degradation characteristics, a vegetation coverage decline sequence is obtained, and the maximum decline magnitude and duration of the vegetation coverage decline sequence are extracted. The initial geological hazard risk value is calculated based on the acceleration ratio, the surface deformation response amplitude, the saturation slope, the maximum decrease amplitude, and the decrease duration.

[0006] Optionally, the step of performing spatial correlation analysis on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics within the high-risk area based on the geological attribute information to obtain risk correction factors includes: Based on the stratigraphic lithology information, the distribution boundary of soft rock layers in the high-risk area is obtained; Determine whether the spatial distribution of the deformation acceleration characteristics overlaps with the distribution boundary of the soft rock layer; If the spatial distribution of the deformation acceleration feature overlaps with the boundary of the soft rock layer distribution, the average deformation rate of the deformation acceleration feature points in the overlapping area is obtained as the first average rate, and the overall average rate of all deformation acceleration feature points in the high-risk area is obtained as the second average rate. The ratio of the first average rate to the second average rate is calculated as the lithological amplification factor; Based on the geological structural information, the centerline of the fault zone is extracted; Calculate the vertical distance from each region with significant rainfall accumulation characteristics to the center line of the fault zone to obtain multiple vertical distance values, wherein the region with significant rainfall accumulation characteristics is a spatially continuous region where the accumulated rainfall exceeds a preset accumulated threshold. The average of the multiple vertical distance values ​​is used to construct the association distance; Risk correction factors are obtained based on the lithological amplification factor, the structural correlation distance, the water saturation characteristics, and the vegetation degradation characteristics.

[0007] Optionally, obtaining the risk correction factor based on the lithological amplification factor, the structural correlation distance, the water saturation characteristics, and the vegetation degradation characteristics includes: Based on the distribution information of the historical disaster points, the water saturation characteristic value at the time of occurrence of each historical disaster point is obtained; Determine whether the difference between the water saturation characteristic value in the current high-risk area and the water saturation characteristic value at the time of the occurrence of the historical disaster point is within a preset difference range; If there are historical disaster points whose difference is within the preset difference range, then the number of such historical disaster points is obtained as the number of associated disaster points, and the total number of such historical disaster points is obtained. The ratio of the number of associated disaster points to the total number is calculated as the water content similarity. Based on the aforementioned vegetation degradation characteristics, degraded patches with a vegetation coverage decline rate exceeding a preset decline rate threshold are extracted. Determine whether there is a spatial superposition relationship between the degraded plaque and the deformation acceleration feature; If the aforementioned spatial superposition relationship exists, the area of ​​the superposition region is obtained as the superposition area, and the total area of ​​the degraded patches is obtained as the total degradation area; The ratio of the superimposed area to the total degraded area is calculated as the vegetation deformation correlation degree; The risk correction factor is determined based on the lithological amplification factor, the structural correlation distance, the water content similarity, and the vegetation deformation correlation.

[0008] Optionally, after marking the target area as a high-risk area, the method further includes: Obtain meteorological forecast data corresponding to the high-risk area, including hourly rainfall forecasts and hourly wind forecasts for a preset future period; Based on the hourly rainfall forecast values, the cumulative forecast rainfall is calculated, and it is determined whether the cumulative forecast rainfall exceeds a preset cumulative rainfall threshold. If the cumulative forecast rainfall exceeds the preset cumulative rainfall threshold, then the maximum hourly rainfall intensity is extracted from the hourly rainfall forecast values; Based on the hourly wind forecast values, calculate the average wind speed and the maximum instantaneous wind speed; Determine whether the current surface deformation rate of the high-risk area exceeds a preset deformation rate threshold; If the surface deformation rate exceeds the preset deformation rate threshold, then a dynamic risk increment value is obtained based on the maximum hourly rainfall intensity and the maximum instantaneous wind speed. The target geological hazard risk value is updated based on the dynamic risk increment value.

[0009] Optionally, obtaining the dynamic risk increment value includes: Determine whether the maximum hourly rainfall intensity exceeds a first intensity threshold; If the maximum hourly rainfall intensity exceeds the first intensity threshold, then the slope and aspect information of the high-risk area are obtained; Based on the slope information and the aspect information, calculate the surface runoff convergence coefficient of the high-risk area under the maximum hourly rainfall intensity; The rainwater infiltration rate is calculated based on the surface runoff convergence coefficient and the maximum hourly rainfall intensity. Determine whether the maximum instantaneous wind speed exceeds the first wind speed threshold; If the maximum instantaneous wind speed exceeds the first wind speed threshold, then obtain the vegetation cover type and average vegetation height of the high-risk area; Based on the vegetation cover type and the average vegetation height, calculate the wind load transfer coefficient of the high-risk area under the maximum instantaneous wind speed; The dynamic risk increment value is calculated based on the rainwater infiltration rate and the wind load transfer coefficient.

[0010] Optionally, after generating the geological disaster early warning information, the method further includes: Obtain real-time surface deformation data corresponding to multiple monitoring points within the high-risk area; Based on the real-time surface deformation data of each monitoring point, the real-time deformation rate of each monitoring point is calculated as the first rate; Determine whether there is a monitoring point where the first rate continuously increases and exceeds the second rate for three consecutive observation cycles as the target monitoring point, where the second rate is the deformation rate within each segment. If the target monitoring point exists, obtain the spatial coordinates of the target monitoring point, and calculate the spatial distribution density of all target monitoring points based on the spatial coordinates; Determine whether the spatial distribution density exceeds a preset density threshold; If the spatial distribution density exceeds the preset density threshold, then the deformation acceleration direction corresponding to all target monitoring points is extracted; Calculate the angle between every two deformation acceleration directions, and obtain the average of all angles as the direction consistency index; If the directional consistency index is less than the preset consistency threshold, it is determined that there is an overall slippage trend, and the warning level of the geological disaster early warning information is upgraded by one level.

[0011] Optionally, after determining that there is an overall slippage trend, the method further includes: Obtain digital elevation model data for the high-risk area; Based on the digital elevation model data, the slope runoff network of the high-risk area is extracted; Determine whether the spatial distribution of the target monitoring points coincides with the gully lines in the slope runoff network; If the spatial distribution of the target monitoring points coincides with the valley line, then the longitudinal slope of the valley in the overlapping section is obtained; Based on the longitudinal slope of the valley, the potential slip velocity is calculated; Determine whether the potential slip speed exceeds a preset speed threshold; If the potential slip velocity exceeds the preset velocity threshold, then the protection target information downstream of the high-risk area is obtained, and the protection target information includes the protection target type and the distance to the protection target; Calculate the impact range level based on the type of protected target and the distance to the protected target; Based on the aforementioned impact range level, the warning level of the geological disaster early warning information is further upgraded.

[0012] Optionally, after generating the geological disaster early warning information, the method further includes: Obtain emergency resource data for the administrative region where the high-risk area is located, including the location of emergency teams, the location of emergency material warehouses, and the location of emergency shelters; Based on the spatial range of the high-risk area, calculate the first arrival time of each emergency team from its location to the nearest edge of the high-risk area; Based on the location of each emergency supplies warehouse, calculate the second arrival time of supplies to the high-risk area; Based on the location of each emergency shelter, calculate the evacuation time from the residential area within the high-risk zone to the nearest emergency shelter; If the first arrival time exceeds the first time threshold, or the second arrival time exceeds the second time threshold, or the evacuation time exceeds the third time threshold, then resource gap information is generated; The resource shortage information and the geological disaster early warning information are combined and sent to the superior emergency management platform.

[0013] In summary, this application includes the following beneficial technical effects: By acquiring and spatiotemporally registering time-series surface deformation, rainfall, soil moisture, and vegetation cover data, this application integrates four types of multi-source monitoring data—deformation, rainfall, soil moisture, and vegetation cover—into a unified spatiotemporal reference system. This addresses the technical shortcomings of existing technologies, such as insufficient integration of multi-source monitoring data and difficulty in simultaneously utilizing multi-dimensional information for comprehensive assessment. Furthermore, by extracting features related to deformation acceleration, rainfall accumulation, water saturation, and vegetation degradation, initial geological hazard risk values ​​are calculated. Spatial correlation analysis, combined with stratigraphic lithology, geological structure, and historical hazard point distribution information, yields risk correction factors, further improving the accuracy of risk assessment. Simultaneously, after generating early warning information, this application can dynamically identify target monitoring points, determine overall slippage trends, and upgrade the early warning level based on real-time monitoring data. This achieves dynamic updating capabilities for early warning information, resolving the problem of existing technologies where early warning information cannot be adjusted in real-time based on the latest monitoring data, resulting in early warning lag. Attached Figure Description

[0014] Figure 1 This is a flowchart of the main process of a dynamic assessment method for geological disaster risk based on multi-source spatiotemporal big data according to an embodiment of this application. Detailed Implementation

[0015] Firstly, this application discloses a dynamic assessment method for geological disaster risks based on multi-source spatiotemporal big data.

[0016] Reference Figure 1 A dynamic assessment method for geological hazard risk based on multi-source spatiotemporal big data, comprising steps S101 to S109: Step S101: Obtain multi-source spatiotemporal data of the target area. The multi-source spatiotemporal data includes time-series surface deformation data, time-series rainfall data, time-series soil moisture data, and time-series vegetation cover data.

[0017] Specifically, the target area refers to a continuous geographic spatial range requiring geological hazard risk assessment, which can be an independent landslide body, a complete watershed unit, or an administrative region. Time-series surface deformation data refers to time-series surface displacement data obtained using synthetic aperture radar interferometry (SAR). This technique calculates the cumulative deformation of each observation point relative to the initial time at each observation moment by processing the phase difference information of multi-temporal satellite radar imagery. The deformation measurement accuracy can reach the millimeter level, and the data format is the cumulative deformation value corresponding to each observation moment. Time-series precipitation data refers to time-series precipitation depth data obtained from meteorological satellite inversion products or ground rain gauge networks. Each data point records the rainfall within a unit time window, measured in millimeters. The time window is typically daily or hourly. Time-series soil moisture content data refers to the time series data of surface soil volumetric moisture content obtained through microwave remote sensing satellite inversion. Microwave signals are very sensitive to the soil dielectric constant, which is directly related to soil moisture content. Therefore, soil moisture content can be obtained through microwave brightness temperature data inversion. Each data point represents the volume percentage of water in a unit volume of soil, with a value ranging from 0% to 100%. Time-series vegetation cover data refers to the time series data of the proportion of the vertical projection area of ​​green vegetation canopy to a unit land area, calculated based on optical satellite imagery. It is usually converted into vegetation cover by calculating the normalized vegetation index, with a value ranging from 0 to 1, where 0 represents bare soil or water without vegetation cover, and 1 represents complete vegetation cover.

[0018] Step S102: Perform spatiotemporal registration on time-series surface deformation data, time-series rainfall data, time-series soil moisture data, and time-series vegetation cover data to obtain multiple spatiotemporal matching sequences.

[0019] Specifically, spatiotemporal registration refers to the data processing process of unifying data from different sources and with different spatiotemporal resolutions under the same spatiotemporal reference system. It includes two stages: spatial registration and temporal registration. In spatial registration, all data are unified to the same coordinate system, typically using the WGS84 geographic coordinate system or the universal transverse Mercator projection coordinate system. Data with different spatial resolutions are resampled to the same grid size using bilinear interpolation. The grid size is determined based on the target area and monitoring accuracy requirements, with a typical value of 30 meters by 30 meters. Bilinear interpolation is an interpolation method that calculates the value of a point by weighted averaging the values ​​of its four nearest neighbors. The weights depend on the distance from the interpolated point to each neighbor. In temporal registration, the data with the lowest temporal resolution is used as the baseline, and other data are interpolated to the same time node using linear interpolation. Linear interpolation assumes that the data values ​​between two known time points change linearly with time, calculating the value at intermediate times through proportional relationships. After spatiotemporal registration, each grid cell simultaneously possesses four attributes at each time node: deformation value, rainfall, soil moisture content, and vegetation cover, forming a quintuple. The set of quintuples of all grid cells at all time nodes constitutes a spatiotemporal matching sequence, and each spatiotemporal matching sequence corresponds to a set of data records of a grid cell in the time dimension.

[0020] Step S103: Based on each spatiotemporal matching sequence, extract features of deformation acceleration, rainfall accumulation, water saturation, and vegetation degradation respectively.

[0021] Specifically, deformation acceleration features refer to the stage-specific characteristic of a continuously increasing deformation rate identified from deformation time series. The extraction method involves performing second-order difference calculations on the deformation time series of each grid cell. The first-order difference of the deformation time series represents the deformation rate, and the second-order difference represents the acceleration of the deformation rate. When the second-order difference value is positive for three consecutive time nodes and exceeds a preset threshold of 0.05 mm / s² day, the location is determined to have deformation acceleration features. Simultaneously, the acceleration start time, acceleration duration, and average acceleration value of the acceleration phase are recorded. Rainfall accumulation features refer to the characteristic of accumulated rainfall at different time scales exceeding the historical normal level for the same period, calculated from rainfall time series. The extraction method involves using the sliding window method to calculate accumulated rainfall at three scales: 15 days, 30 days, and 60 days. The sliding window method involves sliding a fixed-length window along the time axis and calculating the cumulative value of the data within the window. When the accumulated rainfall exceeds the 85th percentile of the historical average for the same period in that region, it is marked as having significant rainfall accumulation features.

[0022] Soil water saturation refers to a state where soil moisture content consistently exceeds a certain proportion of field capacity. Field capacity refers to the maximum moisture content that soil can retain after gravity drainage, and it is an important characteristic value of soil moisture movement. Field capacity varies among different soil types: approximately 10% to 15% for sandy soil, 20% to 30% for loam, and 30% to 40% for clay. When soil moisture content consistently exceeds 90% of field capacity for more than 72 hours, it is considered water saturation, and the onset time and duration of saturation are recorded. Vegetation degradation refers to a trend of continuous decline in vegetation cover. It is extracted by comparing vegetation cover in different years during the same season and calculating the interannual rate of change of vegetation cover. The interannual rate of change is the difference between the current year's vegetation cover and the vegetation cover of the same season in the previous year, divided by the previous year's vegetation cover. When the interannual rate of change is negative and its absolute value exceeds 0.15, vegetation degradation is considered to exist, and the onset time, magnitude, and duration of degradation are recorded.

[0023] Step S104: Calculate the initial geological hazard risk value based on deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics.

[0024] Specifically, the initial geological hazard risk value refers to a quantitative indicator reflecting the current probability of a disaster occurring, calculated directly based solely on multi-source monitoring data without considering background geological conditions. The value ranges from 0 to 100. The calculation method consists of three steps: First, each characteristic is classified and quantitatively scored. The deformation acceleration characteristic is divided into four levels based on the average acceleration during the acceleration phase: acceleration less than 0.01 mm / sqm (0 points), acceleration between 0.01 and 0.05 mm / sqm (0.3 points), acceleration between 0.05 and 0.10 mm / sqm (0.6 points), and acceleration greater than 0.10 mm / sqm (1.0 point). The rainfall accumulation characteristic is divided into four levels based on the multiple by which the accumulated rainfall exceeds the historical threshold: not exceeding the threshold (0 points), exceeding the threshold by 1 to 1.5 times (0.2 points), exceeding by 1.5 to 2 times (0.5 points), and exceeding by more than 2 times (1.0 point).

[0025] Water saturation characteristics are categorized into four levels based on the duration of saturation: unsaturated (0 points), saturation lasting 24-48 hours (0.2 points), saturation lasting 48-96 hours (0.6 points), and saturation lasting over 96 hours (1.0 point). Vegetation degradation characteristics are categorized into four levels based on the degree of degradation: degradation less than 5% (0 points), degradation between 5% and 15% (0.2 points), degradation between 15% and 30% (0.5 points), and degradation greater than 30% (1.0 point). The second step involves weighting and summing the four scores according to preset weighting coefficients: accelerated deformation characteristics have a weight of 0.4, water saturation characteristics have a weight of 0.3, cumulative rainfall characteristics have a weight of 0.2, and vegetation degradation characteristics have a weight of 0.1. The weighting coefficients are determined based on the contribution of each factor to the occurrence of landslide disasters. The third step is to multiply the weighted sum by 100 to obtain the initial geological hazard risk value on a percentage basis.

[0026] Step S105: If the initial geological hazard risk value exceeds the first risk threshold, the target area is marked as a high-risk area.

[0027] Specifically, the first risk threshold is a critical score used to distinguish between general risk areas and high-risk areas. It is determined based on historical disaster statistics and operational experience; in this embodiment, it is set to 60 points. A high-risk area refers to a continuous spatial region where the initial geological hazard risk value exceeds the first risk threshold. The marking of high-risk areas is implemented using a region growing algorithm: First, all grid cells with an initial geological hazard risk value exceeding 60 points are selected as seed points. Then, starting from each seed point, adjacent grid cells are searched in four connected directions (up, down, left, and right). If the risk value of an adjacent grid cell also exceeds 60 points, it is classified into the same high-risk area. This process is repeated until no new grid cells are added. Isolated high-risk areas with an area less than 5000 square meters are considered noise data and are removed from subsequent processing.

[0028] Step S106: Obtain the geological attribute information corresponding to the high-risk area. The geological attribute information includes stratigraphic lithology information, geological structure information, and historical disaster point distribution information.

[0029] Specifically, geological attribute information refers to a spatial data set describing the geological background and disaster history of a target area, used to assess the controlling effect of background geological conditions on the development of geological hazards. Stratigraphic and lithological information refers to data obtained from regional geological maps describing the types of strata, formation ages, weathering degrees, and engineering mechanical parameters. Geological map scales are typically 1:50,000 to 1:200,000. Lithologies closely related to geological hazards include soft rocks such as shale, mudstone, phyllite, and schist, as well as special soils such as weathered granite and loess. The uniaxial compressive strength of soft rocks is generally less than 30 MPa, and they are easily softened and disintegrated when exposed to water. Geological structural information refers to data obtained from active fault distribution maps describing the spatial location, strike, dip, angle, activity age, and activity rate of fault zones. Active faults are those with evidence of activity since the Late Pleistocene; areas within 5 kilometers of active fault zones require close attention. Historical disaster point distribution information refers to a dataset obtained from the geological disaster survey database that records the location, type, scale, occurrence time, triggering factors, and disaster situation of geological disasters that have occurred in history. Each disaster point record includes attribute fields such as longitude coordinates, latitude coordinates, disaster type, volume, occurrence date, and rainfall.

[0030] Step S107: Based on geological attribute information, conduct spatial correlation analysis on the characteristics of accelerated deformation, accumulated rainfall, water saturation, and vegetation degradation in high-risk areas to obtain risk correction factors.

[0031] Specifically, spatial correlation analysis refers to a data processing method that spatially overlays and statistically analyzes monitoring feature layers and geological attribute layers to assess the control effect of geological conditions on disaster development. Risk correction factor refers to an adjustment coefficient used to amplify or reduce the initial geological disaster risk value, with a value ranging from 0.5 to 2.0. A value greater than 1 indicates that geological conditions have a promoting effect, while a value less than 1 indicates an inhibiting effect.

[0032] The spatial correlation analysis comprises four parallel analytical tasks: First, overlaying the deformation acceleration feature point layer with the soft rock layer distribution boundary layer to calculate the proportion of deformation acceleration points falling within the soft rock layer; a higher proportion indicates a stronger control effect of lithology on deformation. Second, establishing multi-level buffer zones with radii of 1 km, 2 km, and 5 km along the fault zone centerline, and calculating the spatial overlap area between areas with significant rainfall accumulation characteristics and each level of buffer zone; a larger overlap area indicates a stronger guiding effect of tectonic activity on rainfall accumulation. Third, comparing the current water saturation feature layer values ​​with the water saturation feature values ​​at the time of each historical disaster point, and statistically analyzing the proportion of disaster points with similar values; a higher proportion indicates a more similarity between current hydrological conditions and historical disaster conditions. Fourth, overlaying the vegetation degradation feature layer with the deformation acceleration feature layer to calculate the spatial overlap rate between vegetation degradation patches and deformation acceleration areas; a higher overlap rate indicates a stronger spatial coupling between vegetation degradation and surface deformation. Then, the weights of the four analysis results are determined by the Analytic Hierarchy Process (AHP), a multi-objective decision analysis method that combines qualitative and quantitative analysis. The weight coefficients of each factor are calculated by constructing a judgment matrix. In this embodiment, the typical values ​​of the four weights are 0.35, 0.30, 0.20, and 0.15. The risk correction factor is obtained after weighted summation.

[0033] Step S108: Calculate the target geological hazard risk value based on the initial geological hazard risk value and the risk correction factor.

[0034] Specifically, the target geological hazard risk value refers to the final quantitative risk index obtained after correction based on geological attribute information, reflecting the probability of hazard occurrence indicated by current monitoring data under a specific geological background. The calculation formula is: target geological hazard risk value equals initial geological hazard risk value multiplied by a risk correction factor. When the risk correction factor is greater than 1, it indicates that geological conditions promote hazard occurrence, and the initial risk value is amplified; when the risk correction factor is less than 1, it indicates that geological conditions inhibit hazard occurrence, and the initial risk value is reduced; when the risk correction factor equals 1, it indicates that geological conditions are moderate, and no adjustment is made to the initial risk. If the calculated result exceeds 100 points, 100 points is taken as the upper limit; if it is less than 0 points, 0 points is taken as the lower limit.

[0035] Step S109: If the target geological hazard risk value exceeds the second risk threshold, a geological hazard early warning message is generated.

[0036] Specifically, the second risk threshold refers to the critical score of the risk value used to trigger the issuance of a geological disaster early warning. It is determined according to the early warning issuance standards and emergency response levels; in this embodiment, it is set to 75 points. Geological disaster early warning information refers to a notification message containing the geographical location, area, risk level, expected impact range, and issuance time of the high-risk area. The risk level is divided into three levels based on the target geological disaster risk value range: a yellow warning is issued between 75 and 85 points, indicating a relatively high probability of disaster, and it is recommended to pay attention and prepare for emergencies; an orange warning is issued between 85 and 95 points, indicating a high probability of disaster, and it is recommended to organize the evacuation of threatened personnel; and a red warning is issued above 95 points, indicating a very high probability of disaster, and it is recommended to take immediate emergency avoidance measures. The generated early warning information includes: the coordinates of the center point and boundary coordinates of the high-risk area, the total area of ​​the high-risk area, the current risk level, the expected risk change trend within the next 72 hours, a list of protected targets that may be affected downstream of the high-risk area, and recommended emergency measures. The early warning information is sent to the mobile terminals of relevant responsible persons via SMS gateways, to village-level monitors via dedicated early warning terminals, and broadcast to the public in the threatened area via an emergency broadcast system.

[0037] In one embodiment of this example, step S104, based on deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics, calculates the initial geological hazard risk value, including steps S201 to S209: Step S201: Obtain the sequence length corresponding to each spatiotemporal matching sequence.

[0038] Specifically, sequence length refers to the number of valid observation data points in each grid cell along the time dimension, i.e., the number of time nodes from the first valid data point to the last valid data point. Valid data points are time nodes where all four types of data are complete and without missing data. This is because synthetic aperture radar interferometry data may have missing data due to phase incoherence caused by vegetation cover or steep terrain, and optical remote sensing data may have missing data due to cloud cover. The sequence length is obtained by iterating through each node on the time axis and checking whether the node simultaneously contains deformation, rainfall, soil moisture content, and vegetation cover data. If all four are present, the count is incremented by one, and the final count is the sequence length. The sequence length is used to determine the minimum window size for subsequent piecewise linear fitting and to assess the sufficiency of the sample size for statistical analysis.

[0039] Step S202: Based on the sequence length, perform piecewise linear fitting on the deformation acceleration feature, obtain the deformation rate within each segment, and obtain the difference in deformation rate between adjacent segments as the rate change value.

[0040] Specifically, piecewise linear fitting refers to a data processing method that divides a deformation time series into several consecutive time periods according to time sequence, and uses the least squares method to fit a straight line within each time period to approximate the deformation change over time. The least squares method is a mathematical optimization technique that solves for the optimal line parameters by minimizing the sum of squared vertical distances between the fitted line and the actual data points. The segmentation is based on the points of change in the deformation trend. The optimal number of segment points is determined using the Bayesian information criterion, which balances model fit and complexity based on the likelihood function and penalty term. Segmentation stops when increasing the number of segment points does not significantly improve the fitting effect. The deformation rate refers to the change in surface deformation per unit time, i.e., the slope of the fitted straight line, measured in millimeters per day or millimeters per month. A positive slope indicates that the surface is moving downhill, and a negative slope indicates that the surface is moving uphill. The rate change value is the absolute value of the difference between the deformation rates fitted from two adjacent time periods, reflecting whether there is a sudden change in the deformation rate. A larger rate change value indicates a more drastic change in the deformation rate, which may be a precursor signal that the disaster body is entering an accelerated destruction phase.

[0041] Step S203: If the rate change value exceeds the preset rate change threshold, it is determined that there is an accelerated deformation stage, and the ratio of the duration of the accelerated deformation stage to the total observation time is obtained as the acceleration percentage.

[0042] Specifically, the preset rate change threshold is an empirical critical value used to determine whether the deformation rate change is significant. It is determined based on the monitoring accuracy of the synthetic aperture radar interferometry data and the regional deformation background value; in this embodiment, it is set to 0.1 mm per day. When the rate change value exceeds this threshold for two consecutive time periods, and the sign of the rate change value is positive (i.e., the rate in the later time period is greater than the rate in the earlier time period), the deformation is determined to have entered the accelerated deformation stage. The accelerated deformation stage refers to the time interval during which the deformation rate continuously increases, from the moment the rate first increases until the moment the rate stops increasing or the data sequence ends. The acceleration proportion is the ratio of the duration of the accelerated deformation stage to the total observation time from the first observation moment to the last observation moment, ranging from 0 to 1. The physical meaning of the acceleration proportion is the proportion of time the target area is in an unstable state. The larger the acceleration proportion, the longer the target area is in an accelerated deformation state, and the higher the probability of overall instability.

[0043] Step S204: Based on the rainfall accumulation characteristics, obtain a sequence of consecutive rainless days and determine whether the length of the sequence of consecutive rainless days exceeds the first length threshold.

[0044] Specifically, a consecutive dry day sequence refers to a time interval consisting of consecutive days with daily rainfall less than 0.1 mm. 0.1 mm is a commonly used threshold in meteorology to distinguish between rainy and dry periods. The method for extracting all consecutive dry day intervals from a rainfall time series is as follows: traverse the time axis, identify all dates with rainfall less than 0.1 mm, connect these dates in chronological order to form consecutive intervals, and record the start date and duration of each interval. The first length threshold is a critical value for distinguishing between short-term and long-term droughts, determined based on regional climate characteristics; 15 days is used in humid regions such as southern China, and 30 days in arid and semi-arid regions such as northwest China. A consecutive dry day sequence exceeding the first length threshold indicates a drought event. After a drought event, the slope soil shrinks and cracks, making it easier for rainwater to infiltrate rapidly along the cracks during rainfall, leading to a sharp increase in pore water pressure and a significantly increased probability of landslides.

[0045] Step S205: If the length of the consecutive rainless day sequence exceeds the first length threshold, then within the preset time window after rainfall recovery, extract the surface deformation response amplitude corresponding to the unit rainfall.

[0046] Specifically, the preset time window after rainfall recovery refers to a specific number of days after the drought ends; in this embodiment, it is 30 days. Within this time window, the slope soil is in a sensitive state after drought, with rapid rainfall infiltration and the most significant response of surface deformation to rainfall. The surface deformation response amplitude corresponding to a unit of rainfall is the ratio of the cumulative deformation caused by each rainfall event to the total rainfall of that event within the time window, expressed in millimeters per millimeter. The specific extraction method is as follows: identify each rainfall event within the preset time window, defined as a time interval consisting of consecutive rainy days, with at least one dry day between two rainfall events; for each rainfall event, record the total rainfall during the event and the cumulative deformation after the event ends; calculate the ratio of deformation to rainfall, and take the maximum value of all event ratios as the surface deformation response amplitude. The physical meaning of this index is the degree of slope deformation response to a unit of rainfall; a larger response amplitude indicates a more unstable slope.

[0047] Step S206: Based on the water saturation characteristics, obtain the soil moisture content time series, and determine whether the soil moisture content time series continuously exceeds the preset moisture content threshold and the duration exceeds the second length threshold.

[0048] Specifically, the soil moisture content time series refers to the data set of soil moisture content values ​​extracted from the spatiotemporal matching sequence, arranged in chronological order. The preset moisture content threshold is a critical value used to determine whether the soil is in a supersaturated state; in this embodiment, it is set to 90% of field capacity. Field capacity refers to the maximum moisture content that soil can maintain after gravity drainage, and it is an important characteristic point on the soil moisture characteristic curve. Field capacity varies for different soil types and needs to be determined through soil texture analysis. The second length threshold is a critical value for the duration at which water saturation is effective; in this embodiment, it is set to 72 hours, or three days. All time intervals in the soil moisture content time series that continuously satisfy a moisture content greater than or equal to 90% of field capacity are found. The duration of each interval is checked; if an interval exists with a duration greater than or equal to 72 hours, it is determined to be a valid water saturation event.

[0049] Step S207: If the soil moisture content time series continuously exceeds the preset moisture content threshold and the duration exceeds the second length threshold, then the slope of the obtained soil moisture content time series curve is the saturation slope.

[0050] Specifically, the saturation slope refers to the slope of a straight line obtained by linearly fitting a curve showing the change in soil moisture content over time for a period of time when the soil moisture content exceeds a preset threshold. The linear fitting uses the least squares method, with time as the independent variable and moisture content as the dependent variable, to obtain a straight line. The slope of this line is the saturation slope. The physical meaning of the saturation slope is as follows: When the saturation slope is positive, it indicates that the moisture content is still slowly increasing, meaning that the rainwater infiltration rate is greater than the drainage rate within the soil, and the pore water pressure continues to accumulate and increase, which is extremely detrimental to slope stability. When the saturation slope is zero, it indicates that the moisture content remains stable, and infiltration and drainage reach a dynamic balance. When the saturation slope is negative, it indicates that the moisture content begins to decrease, the drainage rate exceeds the infiltration rate, and the pore water pressure gradually dissipates, which is beneficial to slope stability. The absolute value of the saturation slope reflects the drastic change in moisture content; the larger the absolute value, the faster the change.

[0051] Step S208: Based on vegetation degradation characteristics, obtain the vegetation coverage decline sequence, and extract the maximum decline magnitude and duration of the vegetation coverage decline sequence.

[0052] Specifically, a vegetation cover decline sequence refers to a data set of continuous periods in which vegetation cover shows a downward trend over time. The method for identifying decline trend segments from a vegetation cover time series is as follows: A sliding window method is used to calculate the local trend, with a window length of 5 time nodes. Within each window, the slope of a linear regression is calculated. If the slope is continuously negative and passes a significance test, it is determined to be a decline trend segment. For each identified decline trend segment, the difference between the vegetation cover at the start and end times is calculated to obtain the decline magnitude. The start and end dates of the decline segment are recorded, and the difference in the number of days between them is calculated to obtain the decline duration. Then, the maximum decline magnitude and its corresponding decline duration are extracted from all decline segments. The combination of these two indicators is used to assess the most severe degree of vegetation degradation. Vegetation degradation leads to a decrease in root system soil-fixing capacity, a reduction in soil shear strength, and increased surface evaporation, thus exacerbating drying shrinkage and cracking, all of which increase the risk of landslides.

[0053] Step S209: Calculate the initial geological hazard risk value based on the acceleration ratio, surface deformation response amplitude, saturation slope, maximum decrease amplitude, and decrease duration.

[0054] Specifically, the initial geological hazard risk value is calculated using a multi-factor comprehensive scoring method. The formula is as follows: the initial risk value equals the acceleration percentage multiplied by the first weighting coefficient of 0.35, plus the surface deformation response amplitude multiplied by the second weighting coefficient of 0.25, plus the absolute value of the saturation slope multiplied by the third weighting coefficient of 0.20, plus the maximum decrease amplitude multiplied by the fourth weighting coefficient of 0.12, plus the decrease duration divided by the standard duration multiplied by the fifth weighting coefficient of 0.08, and then multiplied by 100 to convert to a percentage score. Each weighting coefficient is determined using the analytic hierarchy process (AHP). The AHP judgment matrix is ​​constructed based on expert experience; a consistency ratio of less than 0.1 indicates reasonable weight allocation. The standard duration is 365 days, used to normalize the decrease duration to the range of 0 to 1. The saturation slope is taken as an absolute value because the drastic change in water content, whether rising or falling, reflects the instability of the system, but the sign information of the rise and fall has already been used in other steps; therefore, the absolute value is used here. After calculation, the initial risk value is limited to the range of 0 to 100 points; any value exceeding this range is treated as a boundary value.

[0055] In one embodiment of this example, step S107, based on geological attribute information, performs spatial correlation analysis on the characteristics of accelerated deformation, accumulated rainfall, water saturation, and vegetation degradation in the high-risk area to obtain risk correction factors, including steps S301 to S308: Step S301: Based on stratigraphic lithology information, obtain the distribution boundary of soft rock layers in the high-risk area.

[0056] Specifically, soft rock layers refer to rock types with a uniaxial compressive strength of less than 30 MPa that are easily softened and disintegrated when exposed to water. Uniaxial compressive strength is the ability of a rock to resist axial pressure under unconfined conditions and is a core indicator for evaluating the hardness of a rock, determined through rock mechanics tests. Soft rock layers mainly include mudstone, shale, phyllite, slate, and schist. These rocks contain a large amount of clay minerals, which absorb water and swell upon contact with water, leading to structural damage and a significant decrease in strength. The distribution boundary of soft rock layers refers to the spatial distribution polygonal boundary of soft rock layer units extracted from the stratigraphic lithology information layer after filtering out soft rock layer units based on stratigraphic codes and rock names. In this embodiment, the attribute table of the stratigraphic lithology vector map layer is read, stratigraphic units belonging to soft rock layers are filtered out based on the stratigraphic code field, and the polygonal geometric objects of these units are exported to form a soft rock layer distribution boundary layer, stored in vector format, containing the outer closed outline of each soft rock layer block.

[0057] Step S302: Determine whether the spatial distribution of deformation acceleration characteristics overlaps with the distribution boundary of soft rock layers.

[0058] Specifically, spatial overlap determination refers to using spatial overlay analysis methods in Geographic Information Systems (GIS) to perform spatial intersection queries between the deformation acceleration feature point layer and the soft rock layer distribution boundary layer. Each point in the deformation acceleration feature point layer represents the center location of a grid cell determined to have deformation acceleration features. The calculation method for spatial intersection queries is as follows: traverse each deformation acceleration feature point, and determine whether the point's coordinates fall inside the soft rock layer distribution boundary polygon. The algorithm uses a ray method, i.e., radiate a ray from the point in any direction, and calculate the number of intersection points between the ray and the polygon boundary. If the number of intersection points is odd, the point is inside the polygon; if it is even, it is outside. Count the number of deformation acceleration feature points that fall inside the soft rock layer and compare it with the total number of deformation acceleration feature points to calculate the overlap ratio.

[0059] Step S303: If the spatial distribution of deformation acceleration features overlaps with the boundary of soft rock layer distribution, the average deformation rate of deformation acceleration feature points in the overlapping area is obtained as the first average rate, and the overall average rate of all deformation acceleration feature points in the high-risk area is obtained as the second average rate.

[0060] Specifically, the overlapping region refers to the part where the deformation acceleration feature points and the soft rock layer distribution boundary coincide spatially, that is, the set of deformation acceleration feature points located inside the soft rock layer distribution boundary. The first average rate is calculated as follows: traverse all deformation acceleration feature points located inside the soft rock layer, extract the deformation rate value corresponding to each point, add these rate values ​​together, and divide by the total number of points to obtain the arithmetic mean. The second average rate is calculated as follows: traverse all deformation acceleration feature points in the high-risk area, regardless of whether they are located inside the soft rock layer, extract the deformation rate value of each point, and calculate the arithmetic mean. Comparing the two average rates can reveal the amplification effect of the soft rock layer on the deformation rate.

[0061] Step S304: Calculate the ratio of the first average rate to the second average rate as the lithology amplification factor.

[0062] Specifically, the lithological amplification factor refers to the ratio of the deformation rate of a soft rock layer area to the overall average deformation rate. The calculation formula is the first average rate divided by the second average rate. The physical meaning of the lithological amplification factor is as follows: When the ratio is greater than 1, it indicates that the deformation rate of the soft rock layer area is higher than the overall average level, meaning that the soft rock layer has an amplifying effect on surface deformation. The more concentrated the distribution of soft rock and the higher its softening degree, the larger the amplification factor. When the ratio is equal to 1, it indicates that the deformation rate of the soft rock layer area is comparable to the overall average level, and the soft rock layer has no significant amplification or inhibition effect on deformation. When the ratio is less than 1, it indicates that the deformation rate of the soft rock layer area is lower than the overall average level, meaning that the soft rock layer may be discontinuous or completely weathered in this area, and its control effect on deformation is weak. The typical range of the lithological amplification factor is 0.8 to 2.5, with 2.5 taken when it exceeds 2.5 and 0.8 taken when it is below 0.8.

[0063] Step S305: Extract the center line of the fault zone based on geological structural information.

[0064] Specifically, the centerline of a fault zone refers to the linear geometric element representing the main extension trajectory of a fault structure on a plane. Fault structures are geological formations where crustal rock layers fracture under stress and undergo relative displacement along the fracture surface; they are a crucial controlling factor in geological hazards. Active faults and seismogenic faults are identified from geological structural information. Active faults are those with evidence of activity since the Late Pleistocene, while seismogenic faults are those that have experienced moderate to strong earthquakes. The method for extracting the spatial linear geometric elements of these faults is as follows: Linear elements are read from the geological structural vector layer, and filtered according to the fault type and activity age in the attribute fields to extract line elements that meet the criteria. For fault breccia zones with a certain width, i.e., the band-shaped area formed by the fracturing of rocks on both sides of the fault, its geometric central axis is taken as the centerline. Each centerline contains attribute information such as fault name, strike, dip, dip angle, activity age, and activity rate.

[0065] Step S306: Calculate the vertical distance from each region with significant rainfall accumulation characteristics to the center line of the fault zone, and obtain multiple vertical distance values. The region with significant rainfall accumulation characteristics is a spatially continuous region where the accumulated rainfall exceeds a preset accumulated threshold.

[0066] Specifically, salient areas of accumulated rainfall characteristics refer to continuous spatial regions extracted from the accumulated rainfall characteristic layer using spatial clustering methods where the accumulated rainfall exceeds a preset accumulated threshold. The preset accumulated threshold is 1.5 times the historical average accumulated rainfall for the same period in history, where the historical period refers to the same season or month. The spatial clustering method employs a four-connected region growing algorithm, grouping adjacent grid cells with accumulated rainfall exceeding the threshold into the same salient area. Vertical distance refers to the length of the perpendicular segment drawn from the geometric center of the salient area of ​​accumulated rainfall characteristics to the fault zone centerline. The geometric center refers to the location of the centroid of the polygon in the region; for simple polygons, the centroid coordinates are obtained by calculating the arithmetic mean of the polygon's vertex coordinates. The perpendicular line is the line drawn from the geometric center point to the line containing the fault zone centerline, with the foot of the perpendicular being the intersection of the perpendicular line and the line. Each salient area corresponds to one vertical distance value; therefore, multiple salient areas result in multiple vertical distance values.

[0067] Step S307: Calculate the average of multiple vertical distance values ​​to construct the association distance.

[0068] Specifically, the tectonic correlation distance is the arithmetic mean obtained by summing the vertical distances corresponding to all areas with significant rainfall accumulation characteristics and dividing by the total number of significant areas, expressed in kilometers. The physical meaning of the tectonic correlation distance is to characterize the overall spatial coupling relationship between rainfall accumulation characteristics and fault structures. A smaller tectonic correlation distance indicates that the rainfall accumulation area is closer to the fault zone, and the stronger the guiding effect of the fault zone on rainfall convergence. This is because fault zones often form low-lying or fractured zones, which are conducive to the convergence of surface runoff and groundwater, resulting in higher soil moisture content near the fault zone and a greater risk of geological hazards. Conversely, a larger tectonic correlation distance indicates that the rainfall accumulation area is farther from the fault zone, the fault zone has a weaker control over hydrological processes, and the degree to which geological hazard risk is controlled by tectonic structures is lower.

[0069] Step S308: Obtain the risk correction factor based on the lithological amplification factor, structural correlation distance, water saturation characteristics, and vegetation degradation characteristics.

[0070] Specifically, the preliminary formula for obtaining the risk correction factor is: the risk correction factor equals the lithological amplification factor multiplied by the baseline distance divided by the structural correlation distance, then multiplied by the spatial proportion of water saturation characteristics, and finally multiplied by the spatial proportion of vegetation degradation characteristics. The baseline distance is a reference distance value used to convert the structural correlation distance into a dimensionless adjustment factor. In this embodiment, it is taken as 5 kilometers, meaning that no risk correction is applied when the structural correlation distance is 5 kilometers. The spatial proportion of water saturation characteristics refers to the proportion of the area of ​​areas with significant water saturation characteristics within the high-risk area to the total area of ​​the high-risk area. Areas with significant water saturation characteristics refer to continuous areas where soil moisture content consistently exceeds 90% of field capacity for more than 72 hours. The spatial proportion of vegetation degradation characteristics refers to the proportion of the area of ​​areas with significant vegetation degradation characteristics within the high-risk area to the total area of ​​the high-risk area. Areas with significant vegetation degradation characteristics refer to continuous areas where the annual vegetation cover decline rate exceeds 15%. This formula takes into account the amplification effect of lithology on deformation, the convergence effect of tectonics on rainfall, the spatial distribution of hydrological conditions, and the spatial range of vegetation degradation. The product of these four factors is the risk correction factor.

[0071] In one embodiment of this example, step S308, based on lithological amplification factor, structural correlation distance, water saturation characteristics, and vegetation degradation characteristics, obtains risk correction factors, including steps S401 to S409: Step S401: Based on the distribution information of historical disaster points, obtain the water saturation characteristic value at the time of occurrence of each historical disaster point.

[0072] Specifically, the water saturation characteristic value at the time of a historical disaster refers to the average soil moisture content extracted from historical meteorological data and remote sensing data within three days before and after the date of each disaster. The specific method is as follows: For each historical disaster point, based on its recorded date of occurrence, soil moisture data for the three days before and after that date is retrieved. If multi-source data is available, satellite remote sensing inversion data is used first; if remote sensing data is missing, regional hydrological model inversion data or meteorological station interpolation data is used as a substitute. The soil moisture values ​​within a maximum of seven days are summed and divided by the number of days with valid data to obtain the average soil moisture content. This average value serves as the water saturation characteristic value at the time of the disaster, reflecting the soil moisture conditions at the time of the disaster.

[0073] Step S402: Determine whether the difference between the water saturation characteristic value in the current high-risk area and the water saturation characteristic value at the time of the historical disaster is within the preset difference range.

[0074] Specifically, the preset difference range refers to the allowable deviation interval used to judge the similarity between current and historical hydrological conditions. In this embodiment, it is set to ±5 percentage points, meaning that the current water saturation characteristic value is considered similar when the difference between the current value and the historical value does not exceed 5%. The judgment method is as follows: take the latest soil moisture content data in the current high-risk area as the current water saturation characteristic value. For each historical disaster point, calculate the difference between the current value and the historical value, i.e., the current value minus the historical value, and then take the absolute value of the difference to determine whether the absolute value is less than or equal to 5 percentage points. This judgment is used to screen out historical disaster events similar to the current hydrological conditions.

[0075] Step S403: If there are historical disaster points with differences within the preset difference range, the number of historical disaster points is obtained as the number of associated disaster points, and the total number of historical disaster points is obtained.

[0076] Specifically, the number of associated disaster points refers to the number of disaster points whose absolute difference is less than or equal to 5 percentage points after traversing all historical disaster points. The specific statistical method is as follows: initialize the counter to 0, perform the judgment in step S402 for each historical disaster point, and if the judgment result is true (i.e., the difference is within the preset range), increment the counter by 1. After traversal, the value of the counter is the number of associated disaster points. The total number of historical disaster points refers to the total number of historical disaster points within the high-risk area; this value is fixed and does not change with current conditions. Comparing the number of associated disaster points with the total number of historical disaster points reflects the frequency of current hydrological conditions in historical records.

[0077] Step S404: Calculate the ratio of the number of associated disaster points to the total number as the water content similarity.

[0078] Specifically, water content similarity refers to the ratio of the number of associated disaster sites to the total number of historical disaster sites, ranging from 0 to 1. The physical meaning of water content similarity is: under current soil moisture conditions, what percentage of historical disaster events occurred under similar moisture conditions? A higher water content similarity indicates that the current hydrological conditions are closer to the hydrological conditions that historically triggered disasters, suggesting that the current conditions are within a hydrological range conducive to disaster occurrence, and the likelihood of disaster recurrence is greater. A lower water content similarity indicates that the current hydrological conditions have historically triggered fewer disasters, or that the historical disaster record for the region is incomplete.

[0079] Step S405: Based on vegetation degradation characteristics, extract degraded patches where the vegetation coverage decline rate exceeds a preset decline rate threshold.

[0080] Specifically, the preset degradation rate threshold refers to the critical degradation rate value used to determine whether vegetation degradation is significant. In this embodiment, an annual degradation rate of 15% is used, meaning that areas where vegetation cover decreases by more than 15% within a year are considered degraded patches. The method for extracting degraded patches is as follows: pixels with an annual vegetation cover degradation rate greater than or equal to 15% are selected from the vegetation degradation feature layer. The annual degradation rate is calculated by dividing the difference between the current year's vegetation cover and the vegetation cover of the same season in the previous year by the vegetation cover of the previous year. Then, connectivity analysis is performed on these pixels, and an eight-connected region growing algorithm is used to group adjacent pixels into the same patch. Eight connectivity refers to eight directions: up, down, left, right, upper left, upper right, lower left, and lower right. Each connected region is treated as a degraded patch, and its boundary polygon and area are recorded.

[0081] Step S406: Determine whether there is a spatial superposition relationship between the spatial locations of the degraded plaques and the deformation acceleration features.

[0082] Specifically, spatial overlay relationship refers to whether there is an overlap between the degraded patch layer and the deformation acceleration feature point layer in space. The determination method is as follows: using the spatial overlay analysis function of the geographic information system, an intersection query is performed between the degraded patch polygon layer and the deformation acceleration feature point layer. If at least one deformation acceleration feature point falls inside the degraded patch polygon, a spatial overlay relationship is determined to exist; if all deformation acceleration feature points fall outside all degraded patch polygons, a spatial overlay relationship is determined not to exist.

[0083] Step S407: If there is a spatial superposition relationship, the area of ​​the superposition region is obtained as the superposition area, and the total area of ​​the degraded patches is obtained as the total degraded area.

[0084] Specifically, the superimposed area refers to the spatial overlap between degraded patches and deformation acceleration feature points, that is, the area covered by deformation acceleration feature points located within the degraded patches. Since deformation acceleration feature points are discrete points, the superimposed area is calculated as follows: count the number of deformation acceleration feature points located within the degraded patches, multiply by the area of ​​the grid cell represented by each point. The grid cell area equals the spatial resolution, i.e., 30 meters multiplied by 30 meters equals 900 square meters. The total degradation area refers to the sum of the areas of all degraded patches, that is, the sum of the polygonal areas of each degraded patch.

[0085] Step S408: Calculate the vegetation deformation correlation degree as the ratio of the superimposed area to the total degraded area.

[0086] Specifically, the vegetation deformation correlation degree is the ratio of the superimposed area to the total degraded area, ranging from 0 to 1. The physical meaning of the vegetation deformation correlation degree is: in areas where vegetation degradation has occurred, what proportion of the area also experiences accelerated surface deformation? A higher vegetation deformation correlation degree indicates a stronger spatial coupling between vegetation degradation and surface deformation, and a greater likelihood that vegetation degradation induces or promotes surface deformation; conversely, a lower degree indicates that the two may be driven by different factors. This indicator is used to quantify the contribution of vegetation degradation to geological hazards.

[0087] Step S409: Determine the risk correction factor based on lithological amplification factor, structural correlation distance, water content similarity, and vegetation deformation correlation.

[0088] Specifically, the formula for determining the risk correction factor is: the risk correction factor equals the lithological amplification factor multiplied by the baseline distance divided by the structural correlation distance, then multiplied by the water-bearing similarity, and finally multiplied by the vegetation deformation correlation. The baseline distance is 5 kilometers, consistent with the baseline distance in step S308. The comprehensive physical meaning of this formula is: the risk correction factor is obtained by multiplying four independent factors. The lithological amplification factor reflects the amplification effect of geological lithology on deformation; the reciprocal of the structural correlation distance reflects the control effect of fault zones on disasters; the water-bearing similarity reflects the similarity between hydrological conditions and historical disaster-causing conditions; and the vegetation deformation correlation reflects the spatial coupling relationship between vegetation degradation and deformation. These four factors together determine the direction and magnitude of the correction of the initial risk by the geological background conditions.

[0089] In one embodiment of this example, after marking the target area as a high-risk area in step S105, steps S501 to S507 are further included: Step S501: Obtain the meteorological forecast data corresponding to the high-risk area. The meteorological forecast data includes the hourly rainfall forecast and hourly wind forecast for the future preset time period.

[0090] Specifically, meteorological forecast data refers to predicted weather changes obtained from meteorological departments, typically acquired in real-time through meteorological data interfaces. The future preset time period refers to the time interval extending forward from the current moment; in this example, it is 72 hours, i.e., the next three days. The hourly rainfall forecast value refers to the predicted rainfall depth for each hour within the next 72 hours, measured in millimeters per hour. The data format is a time series, with each time point corresponding to a rainfall value. The hourly wind force forecast value refers to the predicted wind speed for each hour within the next 72 hours, including both average wind speed and maximum gust speed, measured in meters per second. The average wind speed is used to assess sustained wind load, while the maximum gust speed is used to assess the impact of instantaneous wind load on surface stability.

[0091] Step S502: Based on the hourly rainfall forecast, calculate the cumulative forecast rainfall and determine whether the cumulative forecast rainfall exceeds the preset cumulative rainfall threshold.

[0092] Specifically, the cumulative forecast rainfall refers to the total rainfall measured in millimeters, obtained by summing the hourly rainfall forecasts for the next preset time period. The calculation method is as follows: iterate through each hour of the next 72 hours, summing the hourly rainfall forecasts to obtain the cumulative forecast rainfall. The preset cumulative rainfall threshold is the critical value of cumulative rainfall used to trigger dynamic risk increment calculations. It is determined based on historical rainfall intensity and geological disaster-triggered rainfall statistics for the region; in this embodiment, it is set to 100 millimeters. The physical meaning of this threshold is: when the cumulative rainfall over the next three days exceeds 100 millimeters, it is considered a heavy rainfall event, sufficient to significantly impact slope stability, requiring the initiation of a dynamic risk assessment.

[0093] Step S503: If the cumulative forecast rainfall exceeds the preset cumulative rainfall threshold, extract the maximum hourly rainfall intensity from the hourly rainfall forecast values.

[0094] Specifically, the maximum hourly rainfall intensity refers to the highest value among all hourly rainfall forecasts for a predetermined future time period, measured in millimeters per hour. The extraction method involves iterating through each hour of the next 72 hours, recording the hourly rainfall value, and identifying the highest value, which is the maximum hourly rainfall intensity. Maximum hourly rainfall intensity characterizes the extreme nature of rainfall. Short-duration heavy rainfall can quickly form surface runoff and infiltrate, making it more likely to trigger geological disasters than long-duration but low-intensity rainfall.

[0095] Step S504: Calculate the average wind speed and the maximum instantaneous wind speed based on the hourly wind forecast values.

[0096] Specifically, the average wind speed refers to the arithmetic mean of all hourly average wind speed forecasts over a predetermined time period, measured in meters per second. The calculation method is as follows: Iterate through each hour of the next 72 hours, sum the hourly average wind speed values, and divide by 72 to obtain the average wind speed. The maximum instantaneous wind speed refers to the maximum value among all hourly maximum gust wind speed forecasts over a predetermined time period, also measured in meters per second. The extraction method is as follows: Iterate through each hour of the next 72 hours, record the maximum gust wind speed value for each hour, and find the maximum value. The average wind speed is used to assess the long-term impact of continuous wind loads on surface stability, while the maximum instantaneous wind speed is used to assess the instantaneous impact of gust wind loads on surface stability.

[0097] Step S505: Determine whether the current surface deformation rate in the high-risk area exceeds the preset deformation rate threshold.

[0098] Specifically, the current surface deformation rate refers to the deformation rate value obtained from the latest observation in the high-risk area, expressed in millimeters per day. For multiple grid cells within a high-risk area, the average deformation rate of all grid cells is taken as the representative deformation rate of that high-risk area. The preset deformation rate threshold is a critical rate value used to determine whether the deformation has reached a level requiring attention. It is determined based on the regional background deformation value and monitoring accuracy; in this embodiment, it is set to 2 millimeters per month, which is approximately 0.067 millimeters per day. In this embodiment, when the surface deformation rate exceeds this value, it indicates that the slope is already in an unstable state and is more prone to instability and failure under extreme weather conditions.

[0099] Step S506: If the surface deformation rate exceeds the preset deformation rate threshold, obtain the dynamic risk increment value based on the maximum hourly rainfall intensity and the maximum instantaneous wind speed.

[0100] Specifically, the dynamic risk increment value refers to the quantitative indicator of additional risk brought about by extreme weather conditions superimposed on existing deformation, with a value range of 0 to 30. The calculation of the dynamic risk increment value comprehensively considers the weakening effect of rainfall infiltration on soil shear strength and the disturbance effect of wind on surface stability. The specific calculation method is as follows: divide the maximum hourly rainfall intensity by the reference rainfall intensity and multiply by the first coefficient, add the maximum instantaneous wind speed divided by the reference wind speed and multiply by the second coefficient, and then multiply by the normalized value of the deformation rate to obtain the dynamic risk increment value. The reference rainfall intensity is taken as 20 mm / h, the reference wind speed as 15 m / s, the first coefficient as 0.6, the second coefficient as 0.4, and the normalized value of the deformation rate as the ratio of the current deformation rate to the deformation rate threshold, but not greater than 2.

[0101] Step S507: Update the target geological hazard risk value based on the dynamic risk increment value.

[0102] Specifically, updating the target geological hazard risk value means adding the target geological hazard risk value calculated in step S108 to the dynamic risk increment value obtained in step S506 to obtain the final risk value after considering meteorological forecast factors. The calculation formula for the update operation is: the updated risk value equals the original target geological hazard risk value plus the dynamic risk increment value. If the sum exceeds 100 points, 100 points is taken as the upper limit; if the sum is less than 0 points, 0 points is taken as the lower limit. This update operation is used to realize the dynamic rolling assessment of geological hazard risk, enabling early warning information to reflect the impact of future changes in meteorological conditions on the probability of disaster occurrence, and improving the timeliness and accuracy of early warning.

[0103] In one embodiment of this example, step S506, obtaining the dynamic risk increment value, includes steps S601 to S608: Step S601: Determine whether the maximum hourly rainfall intensity exceeds the first intensity threshold.

[0104] Specifically, the first intensity threshold refers to the critical value of hourly rainfall intensity used to distinguish between ordinary rainfall and heavy rainfall. It is determined based on regional climate characteristics and statistical analysis of rainfall intensity that triggers geological disasters; in this embodiment, it is set at 15 mm per hour. This threshold corresponds to the meteorological standard for short-duration heavy rainfall. When the hourly rainfall intensity exceeds 15 mm, the rainfall intensity is relatively high, capable of forming surface runoff and rapidly infiltrating into the soil within a short period. If the maximum hourly rainfall intensity does not exceed this threshold, the impact of rainfall on geological disasters is relatively small, and the contribution of the rainfall component to the dynamic risk increment value can be ignored.

[0105] Step S602: If the maximum hourly rainfall intensity exceeds the first intensity threshold, obtain the slope and aspect information of the high-risk area.

[0106] Specifically, slope information refers to spatial data on the degree of surface inclination, representing the angle between the surface tangent plane and the horizontal plane at each grid cell, measured in degrees, ranging from 0 to 90 degrees. A steeper slope indicates a steeper surface and faster rainwater runoff, but may reduce infiltration. Aspect information refers to the direction of the surface orientation, representing the projection direction of the normal to the surface tangent plane at each grid cell onto the horizontal plane, with true north as 0 degrees and increasing clockwise to 360 degrees. Slope aspect determines the amount of solar radiation received by the surface, affecting vegetation growth and soil moisture evaporation, and also influencing the dominant direction of rainwater runoff. Both slope and aspect information are obtained from digital elevation model data through spatial analysis calculations. Slope is obtained by calculating the ratio of the elevation difference between adjacent grid cells to the horizontal distance, while aspect is obtained by calculating the direction with the greatest rate of elevation change.

[0107] Step S603: Based on slope and aspect information, calculate the surface runoff convergence coefficient of the high-risk area under the maximum hourly rainfall intensity.

[0108] Specifically, the surface runoff convergence coefficient (SFC) is a dimensionless index characterizing the ability of rainfall to form runoff on the surface and converge in a specific area, with a value ranging from 0 to 1. The SFC is calculated as follows: First, the runoff velocity coefficient is calculated based on slope information; the steeper the slope, the faster the runoff velocity. The runoff velocity coefficient is taken as the sine of the slope angle. Second, the convergence direction is calculated based on slope aspect information, dividing the slope aspect into eight main directions, and the connectivity of each grid cell towards the downstream direction is statistically analyzed. Then, the D8 algorithm (eight-direction unidirectional flow algorithm) is used to simulate the surface runoff path. The D8 algorithm assumes that water flows from each grid cell to the cell with the steepest slope among its eight neighboring areas. Finally, the upstream catchment area of ​​each grid cell is calculated, and the ratio of the catchment area to the total area is taken as the SFC for that grid cell. This coefficient reflects the spatial concentration of runoff under the maximum hourly rainfall intensity.

[0109] Step S604: Calculate the rainwater infiltration rate based on the surface runoff convergence coefficient and the maximum hourly rainfall intensity.

[0110] Specifically, rainwater infiltration rate refers to the speed at which rainwater seeps from the surface into the soil per unit time, measured in millimeters per hour. The formula for calculating rainwater infiltration rate is: Rainwater infiltration rate equals the maximum hourly rainfall intensity multiplied by the infiltration coefficient. The infiltration coefficient depends on the surface runoff convergence coefficient and soil type. The specific calculation method for the infiltration coefficient is: Infiltration coefficient equals 1 minus the surface runoff convergence coefficient, then multiplied by the soil infiltration capacity coefficient. The soil infiltration capacity coefficient is determined based on soil texture: 0.8 for sandy soil, 0.5 for loam, and 0.2 for clay. The physical meaning of this formula is: The larger the surface runoff convergence coefficient, the more concentrated the runoff, the deeper the surface water accumulation, and the higher the infiltration rate; the larger the soil infiltration capacity coefficient, the better the soil permeability, and the easier it is for rainwater to infiltrate.

[0111] Step S605: Determine whether the maximum instantaneous wind speed exceeds the first wind speed threshold.

[0112] Specifically, the first wind speed threshold refers to the critical wind speed value used to distinguish between ordinary winds and strong winds. In this embodiment, it is set to 15 meters per second, corresponding to a level 7 wind. When the maximum instantaneous wind speed exceeds 15 meters per second, the force of the wind on surface vegetation and loose deposits cannot be ignored. Strong winds can affect slope stability by blowing down trees, disturbing loose soil, and increasing the dynamic load on slopes. If the maximum instantaneous wind speed does not exceed this threshold, the impact of wind on geological hazards is relatively small, and the contribution of wind load to the dynamic risk increment value can be ignored.

[0113] Step S606: If the maximum instantaneous wind speed exceeds the first wind speed threshold, obtain the vegetation cover type and average vegetation height of the high-risk area.

[0114] Specifically, vegetation cover type refers to the classification of vegetation types within high-risk areas, mainly including forests, shrublands, grasslands, and farmland. Vegetation cover type is obtained from land use cover classification maps, which are typically generated based on multi-temporal optical remote sensing image classification. Different vegetation types exhibit significant differences in root depth, canopy structure, and wind resistance. Average vegetation height refers to the average vertical height of the vegetation canopy within high-risk areas, measured in meters, and is obtained from lidar data or forest resource survey data. The average height of forests is generally 5 to 30 meters, shrublands 1 to 5 meters, and grasslands 0.1 to 1 meter. Vegetation height determines the lever arm of wind loads; the greater the height, the greater the overturning moment generated by the wind load on the ground surface.

[0115] Step S607: Based on the vegetation cover type and average vegetation height, calculate the wind load transfer coefficient of the high-risk area under the maximum instantaneous wind speed.

[0116] Specifically, the wind load transfer coefficient is a dimensionless index used to characterize the proportion of wind load transferred to the ground surface through vegetation roots, with a value ranging from 0 to 1. The calculation method for the wind load transfer coefficient is as follows: First, determine the wind resistance coefficient based on the vegetation cover type: 0.8 for arbor forests, 0.5 for shrub forests, and 0.2 for grasslands. Second, calculate the wind load amplification factor based on the average vegetation height: the amplification factor equals 1 plus the average vegetation height divided by a reference height of 10 meters. Then, calculate the wind load transfer coefficient as the wind resistance coefficient multiplied by the amplification factor, and then multiplied by the wind speed ratio, which is the ratio of the maximum instantaneous wind speed to the reference wind speed, but not greater than 2. In this embodiment, the higher the vegetation, the greater the wind resistance, the higher the proportion of wind load transferred to the ground surface through the root system, and the stronger the disturbance effect on slope stability.

[0117] Step S608: Calculate the dynamic risk increment value based on the rainwater infiltration rate and wind load transfer coefficient.

[0118] Specifically, the formula for calculating the dynamic risk increment is as follows: The dynamic risk increment equals the rainwater infiltration rate divided by the reference infiltration rate, multiplied by the first weight, plus the wind load transfer coefficient multiplied by the second weight, then multiplied by the deformation rate amplification factor, and finally multiplied by the maximum risk increment of 30 points. The reference infiltration rate is taken as 5 mm / h, the first weight as 0.6, the second weight as 0.4, and the deformation rate amplification factor as the ratio of the current deformation rate to the deformation rate threshold, but not greater than 2. The physical meaning of this formula is: the dynamic risk increment is obtained by weighted summation of the contributions of rainfall infiltration and wind load, multiplied by the deformation rate amplification factor to reflect the sensitivity of existing deformation to meteorological factors, and finally scales the result to the range of 0 to 30 points.

[0119] In one embodiment of this example, after generating geological disaster early warning information in step S109, steps S701 to S708 are further included: Step S701: Obtain real-time surface deformation data corresponding to multiple monitoring points in the high-risk area.

[0120] Specifically, monitoring points refer to locations within high-risk areas used for continuous observation of surface deformation. These monitoring points can be permanent scattering body points obtained through synthetic aperture radar interferometry (SAR) processing, or they can be deployed Global Navigation Satellite System (GNSS) monitoring stations or crack gauges. Permanent scattering bodies refer to ground features that maintain stable scattering characteristics in multi-temporal radar imagery, such as buildings and exposed rocks. Real-time surface deformation data refers to the latest deformation information obtained from the monitoring points, including cumulative deformation, deformation rate, and deformation direction. This data is acquired in real-time from the monitoring database via a data interface and is used to identify the presence of abnormal deformation monitoring points.

[0121] Step S702: Based on the real-time surface deformation data of each monitoring point, calculate the real-time deformation rate of each monitoring point as the first rate.

[0122] Specifically, the first rate refers to the deformation rate value at the current moment for each monitoring point extracted from real-time surface deformation data, measured in millimeters per day. For synthetic aperture radar interferometry monitoring points, the deformation rate is obtained through linear fitting of multi-period observation data; for Global Navigation Satellite System monitoring stations, the deformation rate is calculated through differential calculation of the location time series. The first rate reflects the current deformation activity level of the monitoring point and is a core indicator for determining whether the warning level needs to be raised.

[0123] Step S703: Determine whether there is a monitoring point where the first rate continuously increases and exceeds the second rate for three consecutive observation cycles, and use this as the target monitoring point. The second rate is the deformation rate within each segment.

[0124] Specifically, the second rate refers to the deformation rate within each segment obtained in step S202, representing the historical background deformation rate of that monitoring point. "Continuous increase" means that the deformation rate in each of the three consecutive observation periods is greater than the deformation rate in the previous period. The observation period depends on the temporal resolution of the data source; a typical observation period for synthetic aperture radar interferometry data is 12 days. The judgment criteria are: for each monitoring point, extract the first rate values ​​from the most recent three observation periods, and check whether the rate value in each period is greater than the rate value in the previous period, and whether the rate value in each period is greater than the second rate of that monitoring point, i.e., the historical background rate. If both conditions are met simultaneously, the monitoring point is marked as a target monitoring point. A target monitoring point represents an anomaly where deformation is accelerating and has exceeded historical normal levels.

[0125] Step S704: If there are target monitoring points, obtain the spatial coordinates of the target monitoring points and calculate the spatial distribution density of all target monitoring points based on the spatial coordinates.

[0126] Specifically, spatial coordinates refer to the geographic coordinates of the target monitoring points, including longitude and latitude, used to determine the spatial location of the monitoring points. Spatial distribution density refers to the number of target monitoring points per unit area, used to assess the spatial concentration of target monitoring points. The specific calculation method is as follows: First, project the spatial coordinates of all target monitoring points onto a plane coordinate system; second, use the kernel density estimation method to calculate the point density at each location. Kernel density estimation is a non-parametric statistical method that obtains a continuous density surface by placing a kernel function at each point location and superimposing them; then, extract the density values ​​on all grid cells within the high-risk area, and calculate the average density as the spatial distribution density, with units of points per square kilometer.

[0127] Step S705: Determine whether the spatial distribution density exceeds the preset density threshold.

[0128] Specifically, the preset density threshold refers to a critical density value used to determine whether target monitoring points are significantly concentrated in space; in this embodiment, it is set to 5 per square kilometer. The physical meaning of this threshold is: when there are more than 5 acceleration anomaly monitoring points per square kilometer, it indicates that there is a widespread and spatially continuous unstable phenomenon in the area. If the spatial distribution density exceeds the preset density threshold, it indicates that the target monitoring points are not isolated anomalies, but rather form an anomaly region of a certain scale, requiring further analysis of its spatial distribution pattern.

[0129] Step S706: If the spatial distribution density exceeds the preset density threshold, extract the deformation acceleration direction corresponding to all target monitoring points.

[0130] Specifically, the preset density threshold refers to a critical density value used to determine whether target monitoring points are significantly concentrated in space; in this embodiment, it is set to 5 per square kilometer. The physical meaning of this threshold is: when there are more than 5 acceleration anomaly monitoring points per square kilometer, it indicates that there is a widespread and spatially continuous unstable phenomenon in the area. If the spatial distribution density exceeds the preset density threshold, it indicates that the target monitoring points are not isolated anomalies, but rather form an anomaly region of a certain scale, requiring further analysis of its spatial distribution pattern.

[0131] Step S707: Calculate the angle between every two deformation acceleration directions, and obtain the average value of all angles as the direction consistency index.

[0132] Specifically, the angle between any two deformation acceleration directions refers to the minimum angle between the deformation directions of any two points randomly selected from all target monitoring points, expressed in degrees. The formula for calculating the angle is: the angle equals the absolute value of the difference between the two direction angles; if the difference is greater than 180 degrees, then 360 degrees is subtracted from the difference. This process is repeated for all monitoring point pairs, calculating the angle between each pair to obtain a set of angle values. The direction consistency index is the arithmetic mean of this set of angle values, reflecting the degree of consistency in the deformation directions of all target monitoring points. A smaller direction consistency index indicates a more consistent deformation direction across all monitoring points, and a more pronounced overall slip trend; a larger direction consistency index indicates a more dispersed deformation direction, potentially manifesting as localized collapse or discrete settlement.

[0133] Step S708: If the directional consistency index is less than the preset consistency threshold, it is determined that there is an overall slippage trend, and the warning level of the geological disaster warning information is upgraded by one level.

[0134] Specifically, the preset consistency threshold refers to the critical angle value used to determine whether the deformation direction has statistical consistency; in this embodiment, it is set to 30 degrees. When the direction consistency index is less than 30 degrees, it indicates that the deformation direction of all target monitoring points is highly consistent, conforming to the kinematic characteristics of overall slippage. Overall slippage trend refers to the potential for overall sliding failure in high-risk areas. This failure mode is large-scale and severely harmful, requiring an increased emergency response level. Raising the warning level by one level means: if the current warning is yellow, it will be upgraded to orange; if the current warning is orange, it will be upgraded to red; if the current warning is red, it will remain red. After the warning level is upgraded, the corresponding emergency response measures need to be upgraded simultaneously.

[0135] In one embodiment of this example, after determining the existence of an overall slippage trend in step S708, steps S801 to S809 are further included: Step S801: Obtain digital elevation model data for high-risk areas.

[0136] Specifically, digital elevation model (DEM) data refers to datasets that digitally simulate ground topography using limited terrain elevation data, serving as the foundational data for constructing digital terrain models. DEMs are stored in a regular grid format, with each grid point recording the ground elevation value in meters. DEM data can be obtained from Space Shuttle Radar Topographic Mapping Mission (SRTM) data, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER GDEM) global digital elevation model data, or domestically produced high-resolution satellite stereo images. SRTM data is global terrain data acquired by the Synthetic Aperture Radar (SAR) aboard the US Space Shuttle using interferometry, covering the land surface from 60°N to 56°S latitude, with spatial resolutions of 30 meters and 90 meters. ASTER GDEM data is a global digital elevation model jointly released by the Japanese Ministry of Economy, Trade and Industry and NASA, generated based on optical stereo image technology, with a spatial resolution of 30 meters. This embodiment uses 30-meter resolution DEM data to ensure consistency with the spatial resolution of multi-source spatiotemporal data. The elevation value recorded in each grid cell of the digital elevation model data is used for subsequent spatial analysis operations such as slope calculation, aspect calculation, runoff network extraction, and valley gradient calculation.

[0137] Step S802: Based on the digital elevation model data, extract the slope runoff network in the high-risk area.

[0138] Specifically, a slope runoff network refers to the runoff path system formed by surface runoff flowing along a slope under the influence of gravity, including linear topographic features such as gullies, ravines, and streams. The slope runoff network is extracted using the D8 algorithm, an eight-direction unidirectional algorithm that assumes that water flows from each grid cell to the cell with the steepest slope among its eight neighboring cells. The specific extraction steps are as follows: First, the digital elevation model data is filled with depressions. This means filling the concave areas in the digital elevation model to be level with the surrounding terrain to prevent water flow from stagnating in the depressions and causing the runoff network to be interrupted. Second, the flow direction of each grid cell is calculated. The flow direction refers to the direction in which water flows from that cell to its downstream cells, and eight possible directions are represented by power values ​​from 1 to 128. Then, the runoff accumulation of each grid cell is calculated. The runoff accumulation refers to the total number of upstream grid cells flowing to that cell, reflecting the size of the catchment area at that location in the runoff network. Finally, a runoff accumulation threshold is set. In this embodiment, 1000 grid cells correspond to a catchment area of ​​0.9 square kilometers. Grid cells with runoff accumulation exceeding the threshold are extracted as gully lines. Connecting all the gully lines into a network yields the slope runoff network.

[0139] Step S803: Determine whether the spatial distribution of the target monitoring points coincides with the gully lines in the slope runoff network.

[0140] Specifically, a gully line refers to a linear topographic feature where runoff concentrates in a slope runoff network; it is a topographic concave line formed by surface water erosion. The spatial overlap determination of target monitoring points with the gully line employs a buffer zone analysis method. The specific steps are as follows: First, a buffer zone is established for the extracted gully line, with a width of 30 meters, i.e., a range of 15 meters on each side of the gully line. Then, it is determined whether the spatial coordinates of each target monitoring point fall within the gully line buffer zone. If at least one target monitoring point falls within the gully line buffer zone, it is determined that the spatial distribution overlaps with the gully line; if all target monitoring points fall outside the gully line buffer zone, it is determined that there is no overlap. The physical significance of this determination is that landslides often develop along both sides of gullies, and the distribution of monitoring points along the gully is an important spatial characteristic of the overall landslide movement.

[0141] Step S804: If the spatial distribution of the target monitoring points coincides with the valley line, then obtain the longitudinal gradient of the valley in the overlapping section.

[0142] Specifically, the overlapping segment refers to the valley line segment where the target monitoring point is located, that is, the valley interval between the projection points of the upstream and downstream target monitoring points. The valley gradient refers to the elevation drop per unit length along the direction of water flow in the valley; it is a dimensionless topographic parameter, usually expressed as a percentage or per thousand. The calculation method for the valley gradient is as follows: First, extract the start and end points of the overlapping segment valley line. The start point is the valley location corresponding to the upstream target monitoring point, and the end point is the valley location corresponding to the downstream target monitoring point. Second, read the elevation values ​​of the start and end points from the digital elevation model data. Then, calculate the horizontal distance between the start and end points, in meters. Finally, the valley gradient equals the difference between the elevation of the start point and the elevation of the end point, divided by the horizontal distance. A larger valley gradient indicates a steeper valley, greater potential energy of the landslide body, and a faster potential sliding velocity.

[0143] Step S805: Calculate the potential slip velocity based on the longitudinal gradient of the valley.

[0144] Specifically, potential slip velocity refers to the velocity that a landslide body might reach if a landslide were to occur, measured in meters per second (m / s). The potential slip velocity is calculated using the energy method, based on the principle of converting landslide potential energy into kinetic energy. The specific formula is: potential slip velocity equals the acceleration due to gravity multiplied by the longitudinal slope of the valley, multiplied by the square root of the slip distance, and then multiplied by an empirical coefficient. The acceleration due to gravity is taken as 9.8 m / s², the slip distance is the length of the overlapping section of the valley, and the empirical coefficient is between 0.3 and 0.6; in this embodiment, it is taken as 0.45 to account for frictional energy loss and energy dissipation during the sliding process. The physical meaning of this formula is: the greater the longitudinal slope of the valley and the longer the slip distance, the greater the kinetic energy gained by the landslide body, and the faster the potential slip velocity.

[0145] Step S806: Determine whether the potential slip speed exceeds the preset speed threshold.

[0146] Specifically, the preset speed threshold refers to the critical speed value used to determine whether the potential landslide speed has reached a dangerous level. In this embodiment, it is set to 5 meters per second, equivalent to 18 kilometers per hour. The physical meaning of this threshold is: when the landslide speed exceeds 5 meters per second, the landslide body has a large impact kinetic energy, capable of destroying buildings and infrastructure along its path, posing a serious threat to the safety of people downstream. If the potential landslide speed does not exceed 5 meters per second, it means that even if a landslide occurs, its movement speed is relatively slow, and the harm is relatively small; if it exceeds this threshold, it is considered a high-speed landslide, and the warning level needs to be further increased.

[0147] Step S807: If the potential slip velocity exceeds the preset velocity threshold, obtain the protection target information downstream of the high-risk area. The protection target information includes the protection target type and the distance to the protection target.

[0148] Specifically, the downstream area of ​​a high-risk zone refers to the region located below the high-risk zone along the direction of water flow in the valley, typically extending to the valley outlet or major tributary. Protected targets are objects located downstream of the high-risk zone that may be directly impacted by landslides or debris flows. Protected targets are categorized into three types based on importance: the first type includes densely populated residential areas and schools where evacuation is difficult; the second type includes major transportation routes such as highways and railways, and important infrastructure such as power stations and water plants; the third type includes general land resources such as farmland and forest land. The distance to a protected target refers to the straight-line distance or valley length from the edge of the high-risk zone to the location of the protected target, measured in kilometers. Protected target information is extracted from land use planning maps, transportation route maps, and administrative division maps.

[0149] Step S808: Calculate the impact range level based on the type of protected target and the distance to the protected target.

[0150] Specifically, the impact range level is a quantitative indicator used to assess the potential geographical area and severity of a disaster. It is divided into three levels: Level 1, Level 2, and Level 3, with Level 1 having the most severe impact and Level 3 the least. The calculation method for the impact range level is as follows: First, a base level is determined based on the type of protected target. The base level for Category I protected targets is Level 1, for Category II targets it is Level 2, and for Category III targets it is Level 3. Second, adjustments are made based on the distance to the protected targets. If the distance is less than 1 kilometer, the level is adjusted one level towards a more severe impact; if the distance is greater than 5 kilometers, the level is adjusted one level towards a less severe impact; if the distance is between 1 and 5 kilometers, the level remains unchanged. Finally, the adjusted level is the final impact range level. This level reflects the degree of socio-economic impact that a disaster may cause.

[0151] Step S809: Based on the impact range level, further upgrade the warning level of the geological disaster early warning information.

[0152] Specifically, raising the warning level again refers to raising the warning level again based on the impact range level, building upon the level raised by one level in step S708. The specific rules are as follows: if the impact range level is Level 1, the warning level is raised by one level (yellow to orange, orange to red, red remains red); if the impact range level is Level 2, the warning level remains unchanged; if the impact range level is Level 3 and the current warning level is yellow, it remains yellow; if it is orange or red, it remains unchanged. The combined effect of the two warning level upgrades is that, in cases of overall landslide trend, important downstream protected targets, and a relatively fast potential landslide speed, the warning level can be raised directly from yellow to red. The final determination of the warning level is released to the public as an important component of geological disaster early warning information.

[0153] In one embodiment of this example, after generating geological disaster early warning information in step S109, steps S901 to S905 are further included: Step S901: Obtain emergency resource data for the administrative region where the high-risk area is located. The emergency resource data includes the location of emergency teams, the location of emergency material warehouses, and the location of emergency shelters.

[0154] Specifically, administrative regions refer to the county-level or township-level administrative divisions where high-risk areas are located. Emergency resource data is typically statistically analyzed and managed according to administrative regions. Emergency team locations refer to the coordinates of the deployment sites of professional emergency rescue teams, including fire rescue teams, geological disaster investigation teams, and mine rescue teams. Each team's information includes its type, number of personnel, and equipment configuration. Emergency material warehouse locations refer to the coordinates of warehouses storing emergency materials, including rescue tools, life-saving equipment, communication equipment, and living supplies. Each warehouse's information includes the type, quantity, and expiration date of the materials. Emergency shelter locations refer to the coordinates of sites planned for accommodating disaster-stricken populations during disasters, including open spaces such as school playgrounds, gymnasiums, squares, and parks. Each shelter's capacity, facilities, and access routes are recorded. Emergency resource data is obtained from the emergency management department's resource database and needs to be updated regularly to maintain its timeliness.

[0155] Step S902: Based on the spatial range of the high-risk area, calculate the first arrival time of each emergency team from its location to the nearest edge of the high-risk area.

[0156] Specifically, the spatial extent of the high-risk zone refers to the boundary polygon of the high-risk zone marked in step S105, encompassing all grid cells covered by the high-risk zone. The nearest edge refers to the endpoint of the shortest path from the emergency team's location to the boundary polygon of the high-risk zone, i.e., the point on the boundary of the high-risk zone closest to the emergency team. The first arrival time is calculated as follows: First, a network analysis algorithm is used to calculate the shortest driving path from the emergency team's location to the nearest edge of the high-risk zone. The network analysis is based on road network data, considering factors such as road grade, traffic speed, and traffic conditions. Then, the shortest path length is divided by the expected driving speed to obtain the travel time. The expected driving speed is determined according to the road grade: 80 km / h for expressways, 60 km / h for national highways, 40 km / h for provincial highways, and 20 km / h for county and township roads. Finally, the travel time is added to the team assembly time of 15 minutes to obtain the first arrival time, in minutes. The first arrival time reflects the fastest speed at which emergency rescue forces can reach the disaster site.

[0157] Step S903: Based on the location of each emergency supplies warehouse, calculate the second arrival time of supplies to be allocated to high-risk areas.

[0158] Specifically, material allocation refers to the process of transporting emergency supplies from warehouses to the edge of high-risk areas. The calculation method for the second arrival time is similar to that of the first arrival time, but it needs to consider the loading time. The specific calculation steps are as follows: First, for each emergency supply warehouse, calculate the shortest driving route and driving time from the warehouse location to the nearest edge of the high-risk area; then, estimate the loading time based on the type and quantity of supplies, generally taking 30 minutes for general supplies and 60 minutes for large equipment; finally, add the driving time to the loading time to obtain the second arrival time, in minutes. If multiple emergency supply warehouses exist, the minimum second arrival time among all warehouses is taken as the material allocation time to the high-risk area. The second arrival time reflects the fastest speed at which emergency supplies can reach the disaster site.

[0159] Step S904: Based on the location of each emergency shelter, calculate the evacuation time from the residential area in the high-risk zone to the nearest emergency shelter.

[0160] Specifically, settlements refer to human dwellings within or around high-risk areas that may be affected by disasters, extracted from settlement distribution maps or census data. Evacuation time refers to the time required for threatened residents to move from their residences to safe shelters, measured in minutes. The calculation method for evacuation time is as follows: First, for each settlement within the high-risk area, calculate the walking distance from that settlement to each emergency shelter. The walking path is based on road network and walkability analysis, considering road accessibility and terrain undulation. Then, divide the walking distance by the walking speed to obtain the walking time. The walking speed is determined based on the population composition, generally taken as 1.2 meters per second (approximately 4.3 kilometers per hour). Finally, add the resident response time of 5 minutes to the walking time to obtain the evacuation time. For each settlement, the minimum evacuation time corresponding to all emergency shelters is taken as the evacuation time for that settlement. Then, the maximum evacuation time of all settlements within the high-risk area is taken as the overall evacuation time index for the high-risk area.

[0161] Step S905: If the first arrival time exceeds the first time threshold, or the second arrival time exceeds the second time threshold, or the evacuation time exceeds the third time threshold, then resource gap information is generated.

[0162] Specifically, the first time threshold refers to the acceptable upper limit of the arrival time of emergency teams, determined according to the geological disaster emergency response specifications. In this embodiment, it is set to 60 minutes, meaning that emergency teams should arrive at the disaster site within one hour of the warning being issued. The second time threshold refers to the acceptable upper limit of the emergency material allocation time, set to 120 minutes in this embodiment, meaning that materials should arrive at the disaster site within two hours. The third time threshold refers to the acceptable upper limit of the evacuation time for residents, set to 30 minutes in this embodiment, meaning that threatened people should be transferred to safe shelters within 30 minutes. Resource gap information refers to notification information that identifies which specific indicator has a gap and the degree of the gap when the emergency response capacity is insufficient. The generation logic of resource gap information is as follows: if the first arrival time exceeds 60 minutes, it is recorded that there is a gap in the emergency team's response capacity, and the degree of the gap is the ratio of the time exceeding the threshold; if the second arrival time exceeds 120 minutes, it is recorded that there is a gap in the emergency material allocation capacity; if the evacuation time exceeds 30 minutes, it is recorded that there is a gap in the evacuation capacity. Resource gap information also includes the specific location of the gap and the number of people affected.

[0163] Step S906: Combine the resource shortage information with the geological disaster early warning information and send it to the superior emergency management platform.

[0164] Specifically, the superior emergency management platform refers to the emergency management agency at the next higher level in the administrative region where the high-risk area is located. For example, the superior of the county-level platform is the city-level platform, and the superior of the city-level platform is the provincial-level platform. Merged transmission refers to integrating geological disaster early warning information and resource gap information into a single, complete emergency information message, which is then sent to the superior emergency management platform via a dedicated emergency communication network. The merged information includes: basic information about the high-risk area, such as its location, area, and risk level; the early warning level, such as yellow, orange, or red; a detailed description of the resource gap, such as the number of minutes needed for emergency teams, supplies, and evacuation; and specific requests for superior coordination and support, such as requesting additional emergency teams, allocation of emergency supplies, and assistance in organizing the evacuation of residents. Upon receiving the merged information, the superior emergency management platform initiates a cross-regional emergency resource dispatch procedure based on the resource gap situation, dispatching support forces and allocating scarce resources to the gap area. This step achieves the linkage between early warning information and emergency response resources, improving the collaborative efficiency of geological disaster emergency management.

[0165] The above are all 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 dynamic assessment method for geological hazard risk based on multi-source spatiotemporal big data, characterized in that, include: Acquire multi-source spatiotemporal data of the target area, including time-series surface deformation data, time-series rainfall data, time-series soil moisture data, and time-series vegetation cover data; Spatiotemporal registration is performed on the time-series surface deformation data, the time-series rainfall data, the time-series soil moisture data, and the time-series vegetation cover data to obtain multiple spatiotemporal matching sequences; Based on each spatiotemporal matching sequence, features of deformation acceleration, rainfall accumulation, water saturation, and vegetation degradation are extracted respectively. Based on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics, the initial geological hazard risk value is calculated; If the initial geological hazard risk value exceeds the first risk threshold, the target area will be marked as a high-risk area; Obtain the geological attribute information corresponding to the high-risk area, including stratigraphic lithology information, geological structure information, and historical disaster point distribution information; Based on the geological attribute information, spatial correlation analysis is performed on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics in the high-risk area to obtain risk correction factors; Calculate the target geological hazard risk value based on the initial geological hazard risk value and the risk correction factor; If the target geological hazard risk value exceeds the second risk threshold, a geological hazard early warning message will be generated.

2. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 1, characterized in that, The calculation of the initial geological hazard risk value based on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics includes: Obtain the sequence length corresponding to each spatiotemporal matching sequence; Based on the sequence length, the deformation acceleration feature is piecewise linearly fitted to obtain the deformation rate within each segment, and the difference in deformation rate between adjacent segments is obtained as the rate change value. If the rate change value exceeds the preset rate change threshold, it is determined that there is an accelerated deformation stage, and the ratio of the duration of the accelerated deformation stage to the total observation time is obtained as the acceleration percentage. Based on the rainfall accumulation characteristics, a sequence of consecutive rainless days is obtained, and it is determined whether the length of the sequence of consecutive rainless days exceeds a first length threshold. If the length of the consecutive rainless day sequence exceeds the first length threshold, then within a preset time window after rainfall recovery, the surface deformation response amplitude corresponding to a unit rainfall amount is extracted. Based on the water saturation characteristics, a soil moisture content time series is obtained, and it is determined whether the soil moisture content time series continuously exceeds a preset moisture content threshold and the duration exceeds a second length threshold. If the soil moisture content time series continuously exceeds the preset moisture content threshold and the duration exceeds the second length threshold, then the slope of the soil moisture content time series curve is obtained as the saturation slope. Based on the vegetation degradation characteristics, a vegetation coverage decline sequence is obtained, and the maximum decline magnitude and duration of the vegetation coverage decline sequence are extracted. The initial geological hazard risk value is calculated based on the acceleration ratio, the surface deformation response amplitude, the saturation slope, the maximum decrease amplitude, and the decrease duration.

3. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 1, characterized in that, Based on the geological attribute information, spatial correlation analysis is performed on the deformation acceleration characteristics, rainfall accumulation characteristics, water saturation characteristics, and vegetation degradation characteristics within the high-risk area to obtain risk correction factors, including: Based on the stratigraphic lithology information, the distribution boundary of soft rock layers in the high-risk area is obtained; Determine whether the spatial distribution of the deformation acceleration characteristics overlaps with the distribution boundary of the soft rock layer; If the spatial distribution of the deformation acceleration feature overlaps with the boundary of the soft rock layer distribution, the average deformation rate of the deformation acceleration feature points in the overlapping area is obtained as the first average rate, and the overall average rate of all deformation acceleration feature points in the high-risk area is obtained as the second average rate. The ratio of the first average rate to the second average rate is calculated as the lithological amplification factor; Based on the geological structural information, the centerline of the fault zone is extracted; Calculate the vertical distance from each region with significant rainfall accumulation characteristics to the center line of the fault zone to obtain multiple vertical distance values, wherein the region with significant rainfall accumulation characteristics is a spatially continuous region where the accumulated rainfall exceeds a preset accumulated threshold. The average of the multiple vertical distance values ​​is used to construct the association distance; Risk correction factors are obtained based on the lithological amplification factor, the structural correlation distance, the water saturation characteristics, and the vegetation degradation characteristics.

4. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 3, characterized in that, The risk correction factor obtained based on the lithological amplification factor, the structural correlation distance, the water saturation characteristics, and the vegetation degradation characteristics includes: Based on the distribution information of the historical disaster points, the water saturation characteristic value at the time of occurrence of each historical disaster point is obtained; Determine whether the difference between the water saturation characteristic value in the current high-risk area and the water saturation characteristic value at the time of the occurrence of the historical disaster point is within a preset difference range; If there are historical disaster points whose difference is within the preset difference range, then the number of such historical disaster points is obtained as the number of associated disaster points, and the total number of such historical disaster points is obtained. The ratio of the number of associated disaster points to the total number is calculated as the water content similarity. Based on the aforementioned vegetation degradation characteristics, degraded patches with a vegetation coverage decline rate exceeding a preset decline rate threshold are extracted. Determine whether there is a spatial superposition relationship between the degraded plaque and the deformation acceleration feature; If the aforementioned spatial superposition relationship exists, the area of ​​the superposition region is obtained as the superposition area, and the total area of ​​the degraded patches is obtained as the total degradation area; The ratio of the superimposed area to the total degraded area is calculated as the vegetation deformation correlation degree; The risk correction factor is determined based on the lithological amplification factor, the structural correlation distance, the water content similarity, and the vegetation deformation correlation.

5. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 1, characterized in that, After marking the target area as a high-risk area, the method further includes: Obtain meteorological forecast data corresponding to the high-risk area, including hourly rainfall forecasts and hourly wind forecasts for a preset future period; Based on the hourly rainfall forecast values, the cumulative forecast rainfall is calculated, and it is determined whether the cumulative forecast rainfall exceeds a preset cumulative rainfall threshold. If the cumulative forecast rainfall exceeds the preset cumulative rainfall threshold, then the maximum hourly rainfall intensity is extracted from the hourly rainfall forecast values; Based on the hourly wind forecast values, calculate the average wind speed and the maximum instantaneous wind speed; Determine whether the current surface deformation rate of the high-risk area exceeds a preset deformation rate threshold; If the surface deformation rate exceeds the preset deformation rate threshold, then a dynamic risk increment value is obtained based on the maximum hourly rainfall intensity and the maximum instantaneous wind speed. The target geological hazard risk value is updated based on the dynamic risk increment value.

6. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 5, characterized in that, The acquisition of dynamic risk increment values ​​includes: Determine whether the maximum hourly rainfall intensity exceeds a first intensity threshold; If the maximum hourly rainfall intensity exceeds the first intensity threshold, then the slope and aspect information of the high-risk area are obtained; Based on the slope information and the aspect information, calculate the surface runoff convergence coefficient of the high-risk area under the maximum hourly rainfall intensity; The rainwater infiltration rate is calculated based on the surface runoff convergence coefficient and the maximum hourly rainfall intensity. Determine whether the maximum instantaneous wind speed exceeds the first wind speed threshold; If the maximum instantaneous wind speed exceeds the first wind speed threshold, then obtain the vegetation cover type and average vegetation height of the high-risk area; Based on the vegetation cover type and the average vegetation height, calculate the wind load transfer coefficient of the high-risk area under the maximum instantaneous wind speed; The dynamic risk increment value is calculated based on the rainwater infiltration rate and the wind load transfer coefficient.

7. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 1, characterized in that, After generating the geological disaster early warning information, the following is also included: Obtain real-time surface deformation data corresponding to multiple monitoring points within the high-risk area; Based on the real-time surface deformation data of each monitoring point, the real-time deformation rate of each monitoring point is calculated as the first rate; Determine whether there is a monitoring point where the first rate continuously increases and exceeds the second rate for three consecutive observation cycles as the target monitoring point, where the second rate is the deformation rate within each segment. If the target monitoring point exists, obtain the spatial coordinates of the target monitoring point, and calculate the spatial distribution density of all target monitoring points based on the spatial coordinates; Determine whether the spatial distribution density exceeds a preset density threshold; If the spatial distribution density exceeds the preset density threshold, then the deformation acceleration direction corresponding to all target monitoring points is extracted; Calculate the angle between every two deformation acceleration directions, and obtain the average of all angles as the direction consistency index; If the directional consistency index is less than the preset consistency threshold, it is determined that there is an overall slippage trend, and the warning level of the geological disaster early warning information is upgraded by one level.

8. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 7, characterized in that, After determining that there is an overall slippage trend, the method further includes: Obtain digital elevation model data for the high-risk area; Based on the digital elevation model data, the slope runoff network of the high-risk area is extracted; Determine whether the spatial distribution of the target monitoring points coincides with the gully lines in the slope runoff network; If the spatial distribution of the target monitoring points coincides with the valley line, then the longitudinal slope of the valley in the overlapping section is obtained; Based on the longitudinal slope of the valley, the potential slip velocity is calculated; Determine whether the potential slip speed exceeds a preset speed threshold; If the potential slip velocity exceeds the preset velocity threshold, then the protection target information downstream of the high-risk area is obtained, and the protection target information includes the protection target type and the distance to the protection target; Calculate the impact range level based on the type of protected target and the distance to the protected target; Based on the aforementioned impact range level, the warning level of the geological disaster early warning information is further upgraded.

9. The method for dynamic assessment of geological disaster risk based on multi-source spatiotemporal big data according to claim 1, characterized in that, After generating the geological disaster early warning information, the following is also included: Obtain emergency resource data for the administrative region where the high-risk area is located, including the location of emergency teams, the location of emergency material warehouses, and the location of emergency shelters; Based on the spatial range of the high-risk area, calculate the first arrival time of each emergency team from its location to the nearest edge of the high-risk area; Based on the location of each emergency supplies warehouse, calculate the second arrival time of supplies to the high-risk area; Based on the location of each emergency shelter, calculate the evacuation time from the residential area within the high-risk zone to the nearest emergency shelter; If the first arrival time exceeds the first time threshold, or the second arrival time exceeds the second time threshold, or the evacuation time exceeds the third time threshold, then resource gap information is generated; The resource shortage information and the geological disaster early warning information are combined and sent to the superior emergency management platform.

Images