A method and system for evaluating rainfall erosivity in karst mountainous areas

Abstract

The application discloses a kind of karst mountain area-oriented rainfall erosion force evaluation method and system, it is related to ecological environment monitoring technical field, the method is based on unmanned aerial vehicle orthographic image to construct ground point cloud, the lowest point is selected according to grid to generate ground candidate point, and according to elevation difference and slope determination ground point generates digital elevation model, and the digital elevation model is divided into evaluation unit and exports static ground data;According to the real-time acquisition of hydrological time series data of monitoring device arranged in each evaluation unit, and the erosion rainfall day data is counted;According to erosion rainfall day data, construct rain cluster concentration structure index Gjz, and carry out rain cluster distribution determination;When rain cluster distribution determination is rainfall in time distribution concentration, respectively construct rock-soil discontinuous coupling index Kdo and sink-sink switching index Dsq, and with rain cluster concentration structure index Gjz fitting calculation, construct comprehensive erosion driving value Eqd to carry out comprehensive erosion evaluation.

Description

Technical Field

[0001] This invention relates to the field of ecological environment monitoring technology, specifically to a method and system for assessing rainfall erosivity in karst mountainous areas. Background Technology

[0002] The hydrological environment plays a crucial role in the evolution of landforms and soil migration. Rainfall, as an active factor in the hydrological environment, plays a direct role in the erosion process of surface materials through its functions in slope runoff generation, interface confluence, fissure infiltration, and scour transport. Therefore, research on rainfall erosivity assessment is of practical significance. Especially in karst mountainous areas, where surface soil and rock are interwoven, bare rock is frequently exposed, and fissures and sinkholes are well-developed, rainfall often does not move along a single slope runoff path after reaching the surface. Instead, it is repeatedly distributed among bare rock boundaries, soil layers, fissure openings, and water-retaining depressions. This results in significant differences in erosion responses in different areas under the same rainfall conditions. Therefore, it is necessary to develop rainfall erosivity assessment methods that correspond to the surface structure, water flow switching paths, and rainfall timing characteristics of karst mountainous areas.

[0003] Existing methods for analyzing rainfall erosion in karst mountainous areas primarily rely on conventional indicators such as rainfall amount, intensity, slope, or length. These assessment approaches tend to focus on general slope erosion processes, neglecting the interfacial runoff phenomena under conditions of alternating bare rock and soil layers, the diversion effect of fissures, and the retention processes at the edges of depressions and sinkholes. This makes it difficult to reflect the coexisting hydrodynamic characteristics of surface runoff, infiltration, and transition in karst landscapes. Furthermore, some methods emphasize static topographic parameter calculations, lacking specific characterization of intraday peak surges, duration variations, and short-duration concentrated rainfall events, resulting in inaccurate identification of temporally concentrated rain clusters. Other methods, while incorporating monitoring data, fail to link surface structure parameters with real-time hydrological processes, making it difficult to make targeted determinations of the erosion-driving state in different assessment units during a single rainfall event. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method and system for assessing rainfall erosivity in karst mountainous areas, thus solving the problems mentioned in the background section.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for assessing rainfall erosivity in karst mountainous areas, comprising the following steps: S1. Construct a surface point cloud based on UAV orthophotos, select the lowest point according to the grid to generate ground candidate points, and generate a digital elevation model based on the elevation difference and slope of the ground points. Divide the digital elevation model into evaluation units and output static surface data. S2. Collect hydrological time-series data in real time based on the monitoring devices deployed in each assessment unit, and compile daily data on erosive rainfall. S3. Construct the rain cluster concentration structure index Gjz based on daily data of erosive rainfall, and determine the distribution of rain clusters; S4. When the distribution of rain clusters is determined to be concentrated in time, the discontinuous coupling index of soil and rock and the infiltration switching index Dsq are constructed respectively, and the comprehensive erosion driving value Eqd is constructed by fitting and calculating with the concentrated structure index of rain clusters Gjz to conduct comprehensive erosion assessment.

[0006] Preferably, S1 includes S11; S11. Based on the orthophoto data of the karst mountain area collected by UAV, perform aerial triangulation on the orthophoto data to obtain photogrammetric exterior orientation elements, and perform multi-view dense matching based on the photogrammetric exterior orientation elements to generate surface point cloud data. Divide the surface point cloud data into grids according to the set grid size, and select the lowest elevation point in each grid as ground candidate points. Use the ground candidate points as seed points to calculate the elevation difference and local slope between the remaining elevation points and the adjacent seed points. When the elevation difference is not greater than the preset elevation difference threshold and the local slope is not greater than the preset slope threshold, the current elevation point is determined to be a ground point; otherwise, it is determined to be a non-ground point. Construct an irregular ground triangular network based on all ground points, and perform raster interpolation on the irregular ground triangular network to generate a digital elevation model. The elevation difference threshold and slope threshold are set by using the lowest elevation point in each grid as the ground candidate point, establishing the candidate point neighborhood based on the relationship between adjacent grids, calculating the elevation difference sequence and slope sequence between adjacent ground candidate points, and setting the elevation difference threshold by calculating the sum of the median value and interquartile range of the elevation difference sequence using statistical methods, and setting the slope threshold by calculating the sum of the median value and interquartile range of the slope sequence.

[0007] Preferably, S1 further includes S12; S12. Extract the slope aspect, slope length and flow direction of each grid in the digital elevation model, trace the confluence path according to the flow direction, determine the confluence path with the largest flow volume as the dominant scour line, and divide the target area of ​​the same dominant scour line to obtain several evaluation units. Based on orthophoto data, the boundaries of bare rock, soil layers, and fissures in karst mountain areas are identified. The number of fissures (Ls) within each assessment unit is counted, and the total length (Lr) of the boundary between the bare rock and soil layers and the area (Mj) of each assessment unit are calculated. Then, the bare rock cutting density (Ym) and fissure density (Lm) are calculated separately. , Among them, Ym u Lr u Mj u 、Ls u and Lm uThese represent the bare rock cutting density, total boundary length, area, number of fracture openings, and fracture opening density of the u-th evaluation unit, respectively. Sampling lines were laid out along the main slope, and the lengths of the soil strips intersecting with the soil layer area of ​​each sampling line were counted to obtain a soil strip length sequence. The median value of the soil layer length Lz was obtained by extracting the median of the soil strip length sequence according to the statistical method. The shortest distance Zd from each fracture opening to the dominant scour line was calculated using statistical methods to obtain a distance sequence. The average distance Qj between the fracture opening and the dominant scour line in the distance sequence was then calculated, specifically as follows: , among which, Qj u Zd represents the average distance from the fracture opening to the dominant scour line in the u-th evaluation unit. u，i This represents the shortest distance from the i-th fracture opening in the u-th evaluation unit to the dominant scour line; Based on the bare rock cutting density Ym, the median value of soil layer continuity length Lz, the fracture opening density Lm, and the average distance from the fracture opening to the dominant scour line Qj, output the static surface data of each evaluation unit.

[0008] Preferably, S2 includes S21; S21. Based on the monitoring devices set up in the corresponding karst mountain area of ​​each assessment unit, the hydrological time series data are collected in real time and written into the data table of the corresponding assessment unit according to the location of the monitoring devices. The monitoring devices include a rainfall monitoring device, an interface runoff monitoring device, a fissure runoff monitoring device, and a water level monitoring device at the front edge of a depression or sinkhole. Among them, the rainfall monitoring device is used to collect rainfall sequences, the interface runoff monitoring device is set at the junction of the bare rock area and the soil layer area to collect the runoff volume entering the soil layer area along the junction of the bare rock area and the soil layer area, the fissure mouth runoff monitoring device is set at the upstream confluence of the fissure mouth to collect the runoff volume entering the fissure mouth, and the water level monitoring device is set at the water collection point at the front edge of the depression or sinkhole to collect the water level.

[0009] Preferably, S2 further includes S22 and S23; S22. Divide the rainfall sequence by natural days, calculate the cumulative rainfall for each natural day, obtain the daily rainfall P for each natural day, and within each natural day, perform cumulative calculations on the rainfall sequence within each natural day using a calculation duration of 60 minutes and a sliding step size of 10 minutes. Determine the maximum value among all 60-minute cumulative rainfall values ​​as the maximum continuous 60-minute rainfall H for the natural day. Within each natural day, determine the period with 10-minute rainfall of not less than 0.5 mm as the effective rainfall period, and sum the durations of all effective rainfall periods. The erosion duration T is determined as a natural day; within each natural day, the change in rainfall intensity between two adjacent 10-minute intervals with a difference of not less than 6 mm / h is recorded as a rainfall peak jump, and the total number of rainfall peak jumps N within each natural day is counted; the flow volume recorded by the interface flow monitoring device and the flow volume recorded by the fissure flow monitoring device are combined and divided by the flow area in the corresponding assessment unit to obtain the interface ingress water depth U for the current natural day, and the maximum water level value among the accumulated water levels recorded by the water level monitoring device is determined as the depression retention water depth S for the current natural day; S23. Rainfall days that meet the condition of daily rainfall P≥30mm are defined as erosive rainfall days. The daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, total number of rain peak jumps N, depression water depth S, and interface transition water depth U of the erosive rainfall days are written into the corresponding evaluation unit to form erosive rainfall day data.

[0010] Preferably, S3 includes S31; S31. Construct a rain cluster concentration structure index Gjz based on erosive rainfall day data to reflect the aggregation state of rainfall over time on erosive rainfall days, specifically: Among them, Gjz j N j T j and H j P represents the rain cluster concentration structure index, total number of rain peak jumps, erosion duration, and maximum continuous 60-minute rainfall on the j-th erosive rainfall day. j T0 represents the daily rainfall on the j-th day of erosive rainfall, and T0 represents the time base, which is the standard time scale used when calculating the erosion duration and maximum continuous rainfall of the rainfall process data.

[0011] Preferably, S3 further includes S32; S32. Collect a sample set of rain cluster concentration structure index Gjz under the same historical karst mountain landform type, monitoring area and assessment scale, and calculate the mean and standard deviation by statistical method. Set the sum of the mean and standard deviation as the rain cluster concentration judgment threshold Yg, and then use it to judge the rain cluster distribution with the real-time acquired rain cluster concentration structure index Gjz, as follows. When the rain cluster concentration structure index Gjz ≤ the rain cluster concentration judgment threshold Yg, it indicates that the rainfall on the current rainy day is distributed evenly over time. When the rain cluster concentration structure index Gjz > the rain cluster concentration judgment threshold Yg, it indicates that the rainfall on the current rainy day is concentrated in time.

[0012] Preferably, S4 includes S41; S41. When the distribution of rain clusters is determined to be concentrated in the time of rainfall, construct the soil-rock discontinuity coupling index Kdo and the sinking-infiltration switching index Dsq based on static surface data and daily data of erosive rainfall, respectively. The soil-rock discontinuous coupling index Kdo is used to reflect the degree of participation of karst surface structure in rainfall erosion pathways, specifically: Among them, Kdo u Ym u Lz u 、Lm u and Qj u These represent the soil-rock discontinuity coupling index, bare rock cutting density, median soil layer continuity length, fracture density, and average distance from the fracture to the dominant scour line for the u-th evaluation unit, respectively. The infiltration-sinking switching index Dsq is used to reflect the distribution of water along the soil-rock interface, fissures, and depressions during the current rainfall day, specifically: , among which, Dsq j U j and S j , representing the infiltration-drainage switching index, interface ingress depth, and depression retention depth, respectively, on the j-th day of erosive rainfall.

[0013] Preferably, S4 further includes S42 and S43; S42. The rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo, and the sink-infiltration switching index Dsq are fitted and calculated to construct the comprehensive erosion driving value Eqd, which represents the daily comprehensive erosion driving state of the assessment unit on each erosive rainfall day. Specifically: Where k represents the total number of monitoring days, Eqd u，j This represents the composite erosion driving value of the u-th assessment unit on the j-th day of erosive rainfall; S43. Collect a sample set of comprehensive erosion driving values ​​Eqd under the same historical karst landform type, monitoring area and assessment scale, and calculate the mean and standard deviation using statistical methods. Set the sum of the mean and standard deviation as the erosion judgment threshold Ye, and then perform comprehensive erosion assessment with the real-time acquired comprehensive erosion driving values ​​Eqd, as follows; When the comprehensive erosion driving value Eqd ≤ erosion judgment threshold Ye, it means that the surface structure of the current assessment unit is stable under the influence of rainfall erosion force, and normal monitoring should be maintained at this time. When the comprehensive erosion driving value Eqd > the erosion judgment threshold Ye, it indicates that the rainfall erosion force has an impact on the surface structure of the current assessment unit. At this time, static surface data when the rain cluster distribution is determined to be that the rainfall is distributed evenly in time are collected to perform a second iteration assessment. If the second iteration assessment shows that there is stability, it indicates that the rainfall erosion force is in the erosion accumulation state. If the second iteration assessment still shows an impact, it indicates that the rainfall erosion force is in the erosion destruction state.

[0014] A rainfall erosivity assessment system for karst mountainous areas includes a surface modeling module, a hydrological data acquisition module, a structure determination module, and an erosion assessment module. The surface modeling module is used to construct a surface point cloud based on UAV orthophotos, select the lowest point according to the grid to generate ground candidate points, and determine the ground points based on the elevation difference and slope to generate a digital elevation model. The digital elevation model is divided into evaluation units and static surface data is output. The hydrological acquisition module is used to collect hydrological time-series data in real time based on the monitoring devices deployed in each assessment unit, and to statistically analyze daily data on erosive rainfall. The structure determination module is used to construct the rain cluster concentration structure index Gjz based on the daily data of erosive rainfall, and to determine the distribution of rain clusters; The erosion assessment module is used to construct the soil-rock discontinuity coupling index Kdo and the sink-infiltration switching index Dsq when the rain cluster distribution is determined to be concentrated in time. These are then fitted with the rain cluster concentration structure index Gjz to calculate and construct the comprehensive erosion driving value Eqd for comprehensive erosion assessment.

[0015] This invention provides a method and system for assessing rainfall erosivity in karst mountainous areas. It offers the following advantages: (1) In this method, S1 first performs aerial triangulation and multi-view dense matching based on the orthophoto data collected by the UAV to generate surface point cloud data. Then, the lowest elevation point is selected as the ground candidate point according to the grid division, and the ground point is determined by the elevation difference threshold and the slope threshold to generate a digital elevation model. Subsequently, the slope aspect, slope length and flow direction of each grid of the digital elevation model are extracted, and the confluence path is extracted according to the flow direction. The confluence path with the largest flow volume is determined as the dominant scour line, and the target area is divided into evaluation units. At the same time, the bare rock cutting density Ym, the median value of the continuous soil length Lz, the fissure density Lm and the average distance Qj from the fissure to the dominant scour line of each evaluation unit are statistically analyzed, and static surface data are output. The task completed in this step is to put the spatial relationship between the bare rock area, soil area, fissure and dominant scour line in the karst mountain area into the evaluation unit scale, which is used to describe the basic control of the rock-soil phase, interface cutting and fissure distribution on the erosion path. Compared to current techniques that rely solely on slope, slope length, or conventional topographic relief for analysis, this scheme incorporates the bare rock cutting density Ym, the median length of continuous soil layers Lz, the fissure density Lm, and the average distance Q from the fissure to the dominant scour line into the evaluation basis. This allows the discontinuous features, interface features, and fissure features of the karst surface structure to be expressed within the same evaluation framework.

[0016] (2) In S2 and S3, this method collects rainfall sequences, runoff volumes and water levels in real time based on rainfall monitoring devices, interface runoff monitoring devices, fissure runoff monitoring devices and water level monitoring devices at the front edge of depressions or sinkholes. It also collects daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, total number of rain peak jumps N, interface transition water depth U and depression retention water depth S according to natural days. Then, the rainfall days that meet the condition of daily rainfall P≥30mm are identified as erosive rainfall days, forming erosive rainfall day data. On this basis, the rain cluster concentration structure index Gjz is constructed based on the erosive rainfall day data, and the rain cluster concentration structure index Gjz is compared with the rain cluster concentration judgment threshold Yg to determine whether the rainfall of the current rainfall day is distributed gently or concentrated in time. The task accomplished in this section is to incorporate the total amount, duration, local concentration periods, and abrupt changes of rainfall into the same set of erosive rainfall daily data, and to synchronously record the distribution of water bodies at the interface transition and in depression storage, in order to reflect the temporal organization of erosive rainfall days. Compared with current techniques that mainly use daily rainfall or short-duration rainfall intensity for erosion identification, this scheme, in addition to daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, and total number of rain peak abrupt changes N, also incorporates depression storage depth S and interface transition depth U into the data chain, and uses the rain cluster concentration structure index Gjz to determine the aggregation state of rainfall in the time dimension, thus establishing a complete correspondence between the rainfall process and the karst surface response.

[0017] (3) In S4, when the rain cluster distribution is determined to be concentrated in time, the method constructs the soil-rock discontinuity coupling index Kdo and the infiltration-drainage switching index Dsq based on static surface data and data of days with erosive rainfall, respectively. Then, the rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo, and the infiltration-drainage switching index Dsq are fitted and calculated to construct the comprehensive erosion driving value Eqd. The comprehensive erosion driving value Eqd is compared with the erosion judgment threshold Ye to carry out comprehensive erosion assessment. When the comprehensive erosion driving value Eqd exceeds the erosion judgment threshold Ye, static surface data when the rain cluster distribution is determined to be gradual in time are collected to perform a second iteration assessment to distinguish between the erosion accumulation state and the erosion damage state. The task completed in this step is to integrate the rain cluster concentration structure, the soil-rock discontinuity coupling relationship, and the infiltration-drainage switching relationship into the comprehensive erosion driving value Eqd, which is used to determine the daily comprehensive erosion driving state of the assessment unit on each day of erosive rainfall. Compared to current techniques that rely on a single rainfall or topographic index, this scheme connects the rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo, and the sink-infiltration switching index Dsq. This not only determines whether rainfall erosion has an impact, but also distinguishes between the state of erosion accumulation and the state of erosion destruction. The aim is to perform unit-based determination and state-based identification of rainfall erosion in karst mountainous areas. Attached Figure Description

[0018] Figure 1 This is a schematic diagram illustrating the steps of a rainfall erosivity assessment method for karst mountainous areas according to the present invention. Figure 2 This is a schematic diagram of the process of a rainfall erosivity assessment system for karst mountainous areas according to the present invention; Figure 3 This is a logic diagram of a rainfall erosivity assessment method for karst mountainous areas according to the present invention. Detailed Implementation

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

[0020] Please see Figure 1 This invention provides a method for assessing rainfall erosivity in karst mountainous areas. To achieve the above objectives, this invention employs the following technical solution, comprising the following steps: S1. Construct a surface point cloud based on UAV orthophotos, select the lowest point according to the grid to generate ground candidate points, and generate a digital elevation model based on the elevation difference and slope of the ground points. Divide the digital elevation model into evaluation units and output static surface data. S2. Collect hydrological time-series data in real time based on the monitoring devices deployed in each assessment unit, and compile daily data on erosive rainfall. S3. Construct the rain cluster concentration structure index Gjz based on daily data of erosive rainfall, and determine the distribution of rain clusters; S4. When the distribution of rain clusters is determined to be concentrated in time, the discontinuous coupling index of soil and rock and the infiltration switching index Dsq are constructed respectively, and the comprehensive erosion driving value Eqd is constructed by fitting and calculating with the concentrated structure index of rain clusters Gjz to conduct comprehensive erosion assessment.

[0021] In this embodiment, S1 constructs a surface point cloud based on UAV orthophotos, selects the lowest point according to the grid to generate ground candidate points, and determines the ground points based on elevation difference and slope to generate a digital elevation model. Then, the digital elevation model is divided into evaluation units and static surface data is output. The key here is not just to create a topographic map, but to first separate the spatial differences of surface undulations, exposed bare rock, shallow soil cover, and slope gullies in karst mountainous areas, and then put these differences into individual evaluation units. For example, on the edge slope of a peak-cluster depression, in the orthophotos obtained after a UAV flight, bare rock ridges, thin soil zones, and small runoff channels are often interspersed. If conventional topographic maps or low-resolution elevation data are used directly, many small undulations and local depressions will be merged. After processing according to this method, the lowest point is first selected from the grid as a ground candidate point, and then the ground points are screened according to elevation difference and slope. The resulting digital elevation model can distinguish the main slope, secondary slope breaks, and local depressions, and then evaluation units are delineated according to terrain features. This approach involves identifying surface structure, spatial zoning, and organizing static surface data. The goal is to ensure that each subsequent set of rainfall data corresponds to a specific surface unit, rather than relying on averages across the entire hillside or small watershed. Compared to the common technique of using a single set of regional topographic parameters to represent an entire slope, this step integrates imagery, point clouds, digital elevation models, and assessment unit division into a continuous process, better suited to the spatially diverse characteristics of karst mountain areas. In S2 and S3, the scheme uses monitoring devices deployed in each assessment unit to collect real-time hydrological time-series data and compiles daily data on erosive rainfall. Based on this data, a rain cluster concentration structure index (Gjz) is constructed, and rain cluster distribution is determined. This process can be understood as: first, recording "how much rain fell," and then distinguishing "how the rain fell." For example, in the same daily rainfall event, one assessment unit experiences rainfall in both the morning and afternoon, scattered throughout the day. This differs from another assessment unit where dense rainfall occurs within an hour in the evening, with less rainfall throughout the rest of the day. Even with similar total rainfall, the resulting erosion of the slope will differ. This scheme doesn't just consider the daily total amount; instead, it organizes the real-time collected hydrological time-series data into erosive rainfall daily data, and then further constructs a rain cluster concentration structure index (Gjz) to determine whether the rainfall is concentrated or gradual over time. Taking the upper edge of a sinkhole slope as an example, if the rainfall on a certain day is mainly concentrated within a short period, the rain cluster concentration structure index (Gjz) is high, and the scheme will identify this day as having "concentrated rainfall over time." If the rainfall is dispersed over a longer period, the rain cluster concentration structure index (Gjz) is low, and it will be classified as "gradually distributed rainfall over time." This part undertakes the tasks of rainfall day screening, rainfall time-series identification, and rain cluster distribution determination, aiming to transform simple rainfall observations into process data that can be used for erosion analysis.Compared to current methods that rely solely on cumulative rainfall, single rainfall intensity, or empirical thresholds to assess erosion risk, the Rain Cluster Concentration Structure Index (Gjz) isolates the temporal organization of rainfall events, moving beyond simply assessing "how heavy the rainfall" to focusing on "whether the rainfall is concentrated." In S4, when the rain cluster distribution is determined to be concentrated in time, the scheme further constructs the soil-rock discontinuity coupling index (Kdo) and the sink-infiltration switching index (Dsq), and fits them with the Rain Cluster Concentration Structure Index (Gjz) to calculate a comprehensive erosion driving value (Eqd) for comprehensive erosion assessment. In other words, the scheme does not directly provide a final judgment under all rainfall conditions; instead, it first screens the rainfall cluster concentration structure index (Gjz) and then, under the condition of "concentrated temporal distribution," incorporates surface structure factors and hydrological transformation factors into the analysis. To illustrate with a specific scenario, consider an assessment unit located in the middle of a slope where bare rock and soil layers intersect. The static surface data from the previous stage shows that the soil-rock interface is fractured and the runoff paths are relatively clear. If a concentrated rainfall event with a high rain cluster structure index (Gjz) occurs, the discontinuous coupling relationship between the soil and rock will be reflected through the discontinuous coupling index (Kdo), and the switching state of water between surface runoff and infiltration will be reflected through the runoff-infiltration switching index (Dsq). Finally, these factors, along with the rain cluster structure index (Gjz), contribute to the calculation of the comprehensive erosion driving value (Eqd). The resulting comprehensive erosion driving value (Eqd) is no longer simply a result of "heavy rainfall" or "steep slope," but rather incorporates surface structure, rainfall concentration, and runoff-infiltration conversion into a single decision chain. Compared to current common single-index erosion assessment methods, this scheme has a more complete process: S1 is responsible for static surface data modeling; S2 is responsible for hydrological time-series data collection and daily statistics of erosive rainfall; S3 is responsible for constructing the rain cluster concentration structure index Gjz and determining rain cluster distribution; and S4 is responsible for fitting the soil-rock discontinuity coupling index Kdo and the sink-infiltration switching index Dsq with Gjz to form a comprehensive erosion driving value Eqd for comprehensive erosion assessment. Ultimately, this addresses the problems of "insufficient spatial unit detail, unclear process identification, and coarse comprehensive judgment" in the assessment of rainfall erosivity in karst mountain areas, making the assessment results more suitable for monitoring area identification, key unit screening, and subsequent treatment deployment. Example

[0022] Please refer to Figure 3 Specifically: S1 includes S11; S11. Based on the orthophoto data of the karst mountain area collected by UAV, perform aerial triangulation on the orthophoto data to obtain photogrammetric exterior orientation elements, and perform multi-view dense matching based on the photogrammetric exterior orientation elements to generate surface point cloud data. Divide the surface point cloud data into grids according to the set grid size, and select the lowest elevation point in each grid as ground candidate points. Use the ground candidate points as seed points to calculate the elevation difference and local slope between the remaining elevation points and the adjacent seed points. When the elevation difference is not greater than the preset elevation difference threshold and the local slope is not greater than the preset slope threshold, the current elevation point is determined to be a ground point; otherwise, it is determined to be a non-ground point. Construct an irregular ground triangular network based on all ground points, and perform raster interpolation on the irregular ground triangular network to generate a digital elevation model. The elevation difference threshold and slope threshold are set by using the lowest elevation point in each grid as the ground candidate point, establishing the candidate point neighborhood based on the relationship between adjacent grids, calculating the elevation difference sequence and slope sequence between adjacent ground candidate points, and setting the elevation difference threshold by calculating the sum of the median value and interquartile range of the elevation difference sequence using statistical methods, and setting the slope threshold by calculating the sum of the median value and interquartile range of the slope sequence.

[0023] S1 further includes S12; S12. Extract the slope aspect, slope length and flow direction of each grid in the digital elevation model, trace the confluence path according to the flow direction, determine the confluence path with the largest flow volume as the dominant scour line, and divide the target area of ​​the same dominant scour line to obtain several evaluation units. Based on orthophoto data, the boundaries of bare rock, soil layers, and fissures in karst mountain areas are identified. The number of fissures (Ls) within each assessment unit is counted, and the total length (Lr) of the boundary between the bare rock and soil layers and the area (Mj) of each assessment unit are calculated. Then, the bare rock cutting density (Ym) and fissure density (Lm) are calculated separately. , Among them, Ym u Lr u Mj u 、Ls u and Lm u These represent the bare rock cutting density, total boundary length, area, number of fracture openings, and fracture opening density of the u-th evaluation unit, respectively. Sampling lines were laid out along the main slope, and the lengths of the soil strips intersecting with the soil layer area of ​​each sampling line were counted to obtain a soil strip length sequence. The median value of the soil layer length Lz was obtained by extracting the median of the soil strip length sequence according to the statistical method. The shortest distance Zd from each fracture opening to the dominant scour line was calculated using statistical methods to obtain a distance sequence. The average distance Qj between the fracture opening and the dominant scour line in the distance sequence was then calculated, specifically as follows: , among which, Qj uZd represents the average distance from the fracture opening to the dominant scour line in the u-th evaluation unit. u，i This represents the shortest distance from the i-th fracture opening in the u-th evaluation unit to the dominant scour line; Based on the bare rock cutting density Ym, the median value of soil layer continuity length Lz, the fracture opening density Lm, and the average distance Q from the fracture opening to the dominant scour line, output the static surface data of each evaluation unit. In this embodiment, aerial triangulation is performed based on orthophoto data of karst mountains collected by UAV to obtain photogrammetric exterior orientation elements. Then, multi-view dense matching is performed based on these photogrammetric exterior orientation elements to generate surface point cloud data. Subsequently, the surface point cloud data is divided into grids according to a set grid size, and the lowest elevation point within each grid is selected as a candidate ground point. Each candidate ground point is used as a seed point. The elevation difference and local slope between the remaining elevation points and adjacent seed points are calculated. Points with elevation differences not exceeding an elevation difference threshold and local slopes not exceeding a slope threshold are classified as ground points; otherwise, they are classified as non-ground points. The elevation difference threshold and slope threshold are not directly given. Instead, after establishing the neighborhood of candidate points based on the relationship between adjacent grids, the median and interquartile range of the elevation difference and slope sequences between adjacent ground candidate points are calculated using statistical methods. Then, an irregular triangular mesh is constructed based on all ground points, and rasterization interpolation is performed on the irregular triangular mesh to generate a digital elevation model. On this basis, the slope aspect, slope length, and flow direction of each grid in the digital elevation model are extracted, and the confluence path is extracted according to the flow direction. The confluence path with the largest flow volume is determined as the dominant scour line. The target area of ​​the same dominant scour line is then segmented to obtain several evaluation units. At the same time, bare rock boundaries, soil layers, and fracture openings are identified based on orthophoto data. The number of fissures (Ls), the total length of the boundary between the bare rock area and the soil layer area (Lr), and the area (Mj) of each assessment unit are statistically analyzed to form the bare rock cutting density (Ym) and fissure density (Lm). Sampling lines are laid out along the main slope direction to obtain the soil strip length sequence. The median value of the continuous soil layer length (Lz) is extracted using statistical methods. Then, a distance sequence is formed from the shortest distance (Zd) from each fissure to the dominant scour line, and the average distance (Qj) between the fissure and the dominant scour line is calculated. Finally, based on the bare rock cutting density (Ym), the median value of the continuous soil layer length (Lz), the fissure density (Lm), and the average distance (Qj) between the fissure and the dominant scour line, the static surface data of each assessment unit is output. This step refers to... The aim is to unify the distribution of bare rock and soil layers, the distribution of fissures, the differentiation of slope runoff paths, and the control relationships of dominant scour lines in karst mountainous areas within the scale of assessment units. This allows the construction of digital elevation models, the identification of surface structures, and the division of assessment units to form a coherent data chain. Compared to techniques that rely solely on conventional topographic maps, manual interpretation, or analysis of simple slope parameters, this step combines the ground point screening process with the representation of karst surface structures. This makes the surface elevation reconstruction more closely resemble the actual undulation state, the assessment unit boundaries more closely resemble the real runoff patterns, and the static surface data more completely expresses the cutting of rock-soil interfaces, the discontinuous distribution of soil layers, and the spatial location of fissures. Example

[0024] Please refer to Figure 3 Specifically: S2 includes S21; S21. Based on the monitoring devices set up in the corresponding karst mountain area of ​​each assessment unit, the hydrological time series data are collected in real time and written into the data table of the corresponding assessment unit according to the location of the monitoring devices. The monitoring devices include a rainfall monitoring device, an interface runoff monitoring device, a fissure runoff monitoring device, and a water level monitoring device at the front edge of a depression or sinkhole. Among them, the rainfall monitoring device is used to collect rainfall sequences, the interface runoff monitoring device is set at the junction of the bare rock area and the soil layer area to collect the runoff volume entering the soil layer area along the junction of the bare rock area and the soil layer area, the fissure mouth runoff monitoring device is set at the upstream confluence of the fissure mouth to collect the runoff volume entering the fissure mouth, and the water level monitoring device is set at the water collection position at the front edge of the depression or sinkhole to collect the water level.

[0025] S2 also includes S22 and S23; S22. Divide the rainfall sequence by natural days, calculate the cumulative rainfall for each natural day, obtain the daily rainfall P for each natural day, and within each natural day, perform cumulative calculations on the rainfall sequence within each natural day using a calculation duration of 60 minutes and a sliding step size of 10 minutes. Determine the maximum value among all 60-minute cumulative rainfall values ​​as the maximum continuous 60-minute rainfall H for the natural day. Within each natural day, determine the period with 10-minute rainfall of not less than 0.5 mm as the effective rainfall period, and sum the durations of all effective rainfall periods. The erosion duration T is determined as a natural day; within each natural day, the change in rainfall intensity between two adjacent 10-minute intervals with a difference of not less than 6 mm / h is recorded as a rainfall peak jump, and the total number of rainfall peak jumps N within each natural day is counted; the flow volume recorded by the interface flow monitoring device and the flow volume recorded by the fissure flow monitoring device are combined and divided by the flow area in the corresponding assessment unit to obtain the interface ingress water depth U for the current natural day, and the maximum water level value among the accumulated water levels recorded by the water level monitoring device is determined as the depression retention water depth S for the current natural day; S23. Rainfall days that meet the condition of daily rainfall P≥30mm are defined as erosive rainfall days. The daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, total number of rain peak jumps N, depression water depth S, and interface transition water depth U of the erosive rainfall days are written into the corresponding evaluation unit to form erosive rainfall day data.

[0026] In this embodiment, based on the assessment units obtained from S1, rainfall monitoring devices, interface runoff monitoring devices, fissure runoff monitoring devices, and water level monitoring devices at the front edge of depressions or sinkholes are respectively deployed within the corresponding karst mountain areas. The installation positions of each monitoring device are matched one-to-one with their respective assessment units. Subsequently, the rainfall monitoring devices continuously collect rainfall sequences, the interface runoff monitoring devices collect the runoff volume entering the soil layer area at the boundary between the bare rock area and the soil layer area, the fissure runoff monitoring devices collect the runoff volume entering the fissure, and the water level monitoring devices collect the water level at the catchment location of depressions or sinkholes. The above hydrological time series data are then processed according to... The location of the monitoring device is entered into the data table of the corresponding assessment unit. For example, in an assessment unit located at the intersection of bare rock and soil layers, with a sinkhole's front water collection point below, after continuous monitoring, the rainfall sequence can be divided into natural days, and the daily rainfall P for each natural day can be calculated. Within a single natural day, the maximum continuous 60-minute rainfall H is calculated with a calculation duration of 60 minutes and a slip step size of 10 minutes. Then, the erosion duration T is obtained by accumulating the time periods with 10-minute rainfall of not less than 0.5 mm. The change in rainfall intensity between two adjacent 10-minute time periods with a difference of not less than 6 mm / h is counted as the total number of rain peak jumps N. At the same time, the interface runoff monitoring... The combined flow volume recorded by the flow monitoring device and the flow monitoring device at the fissure opening are divided by the catchment area within the corresponding assessment unit to obtain the interface transition water depth U. The maximum water level value among the accumulated water levels recorded by the water level monitoring device is determined as the depression retention water depth S. Finally, rainfall days that meet the condition of daily rainfall P ≥ 30 mm are determined as erosive rainfall days. The corresponding daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, total number of rain peak jumps N, depression retention water depth S, and interface transition water depth U are written into the corresponding assessment unit to form erosive rainfall day data. The purpose of this step is to integrate rainfall input, interface flow, and fissure flow. By unifying the scattered monitoring information on flow and depression retention into the same assessment unit and the same calendar day scale, subsequent analysis not only focuses on the rainfall record itself, but also simultaneously reflects the transfer status of rainwater between the soil-rock interface, fissure openings, and the front edge of depressions or sinkholes. Compared with existing techniques that rely solely on data from a single rain gauge station and only count cumulative rainfall or conventional rainfall intensity indicators, this step incorporates daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, total number of rain peak jumps N, depression retention depth S, and interface transition depth U into the daily data of erosive rainfall, placing the characteristics of the rainfall process and karst surface runoff within the same data chain. Example

[0027] Please refer to Figure 3 Specifically: S3 includes S31; S31. Construct a rain cluster concentration structure index Gjz based on erosive rainfall day data to reflect the aggregation state of rainfall over time on erosive rainfall days, specifically: Among them, Gjz j N j T j and H j P represents the rain cluster concentration structure index, total number of rain peak jumps, erosion duration, and maximum continuous 60-minute rainfall on the j-th erosive rainfall day. j T0 represents the daily rainfall on the j-th day of erosive rainfall, and T0 represents the time base, which is the standard time scale used when calculating the erosion duration and maximum continuous rainfall of the rainfall process data.

[0028] S3 further includes S32; S32. Collect a sample set of rain cluster concentration structure index Gjz under the same historical karst mountain landform type, monitoring area and assessment scale, and calculate the mean and standard deviation by statistical method. Set the sum of the mean and standard deviation as the rain cluster concentration judgment threshold Yg, and then use it to judge the rain cluster distribution with the real-time acquired rain cluster concentration structure index Gjz, as follows. When the rain cluster concentration structure index Gjz ≤ the rain cluster concentration judgment threshold Yg, it indicates that the rainfall on the current rainy day is distributed evenly over time. When the rain cluster concentration structure index Gjz > the rain cluster concentration judgment threshold Yg, it indicates that the rainfall on the current rainy day is concentrated in time.

[0029] In this embodiment, the total number of rain peak jumps N, erosion duration T, maximum continuous 60-minute rainfall H, and daily rainfall P corresponding to the erosive rainfall day data are first collected. Using time reference T0 as a unified scale, a rain cluster concentration structure index Gjz is constructed to represent the aggregation state of the current erosive rainfall day in the time dimension. This formula reads the total number of rain peak jumps N, erosion duration T, maximum continuous 60-minute rainfall H, and daily rainfall P from the erosive rainfall day data, and performs dimensionless processing of the erosion duration T using time reference T0 to form T0. j / T0 is used to measure the duration of a rainfall event over several standard time units, obtaining a converted value of the duration relative to a uniform time scale, which is then combined with Nj to form... The reasoning is as follows: within the same erosive rainfall day, the greater the total number of rain peak jumps (Nj), the more frequent the fluctuations in rainfall intensity and the more uneven the rainfall distribution over time, resulting in a larger rain cluster concentration structure index (Gjz); the longer the erosion duration (Tj), the more spread out the rainfall distribution. j The larger / T0)+1 is, the lower the rain cluster concentration structure index Gjz will be. Therefore, this part essentially describes the relationship between "jump frequency and duration of expansion"; then calculate Hj / P j This indicates the proportion of the maximum continuous 60-minute rainfall to the total daily rainfall, and is expressed through... Compress the amplification of this percentage to prevent H j / P j An excessively large value exerts an overly strong pull on the rain cluster concentration structure index Gjz. If the proportion of the maximum continuous 60-minute rainfall H in the daily rainfall P is too high, it indicates that the rainfall is more easily concentrated and compressed into a short duration window, and this term will increase accordingly. The square root treatment is used to preserve H. j / P j It differentiates the degree of concentration while avoiding [the situation in H]. j / P j When the value approaches 1, the amplification is too rapid, making it easier to compare the rain cluster concentration structure index Gjz across different erosive rainfall days. Finally, the two parts are multiplied to obtain the rain cluster concentration structure index Gjz. For example, in a certain karst mountainous area, on the j-th erosive rainfall day, monitoring yields N... j =5、T j =180min, T0=60min, H j =42mm, P j =70mm, then T j / T0=3,(N j +1) / ((T j / T0)+1)=6 / 4=1.5,H j / P j =0.6, ≈1.265, therefore Gjz j ≈1.5×1.265≈1.90, indicating that the current day's rainfall process involves both numerous sudden spikes in rainfall peaks and a high proportion of 60-minute rainfall, demonstrating a prominent concentration of rainfall over time; if another day with similar erosive rainfall also has P j =70mm, but N j =1、T j =420min, H j =18mm, then ≈0.28, which shows that under the same daily rainfall conditions, Gjzj is significantly smaller due to fewer sudden jumps, longer duration, and lower peak proportion. Therefore, the reasoning principle of this formula is to combine the "time fluctuations reflected by the total number of rain peak jumps Nj", the "degree of process unfolding reflected by the erosion duration Tj", and the "degree of local compression reflected by the maximum continuous 60-minute rainfall Hj relative to the daily rainfall Pj" into the same expression. The rain cluster concentration structure index Gjz is used to characterize the aggregation state of the rainfall on the erosive rainfall day in the time dimension, which is then compared with the rain cluster concentration judgment threshold Yg to distinguish whether the rainfall on the current rainfall day is distributed gently or concentrated in time. Subsequently, a sample set of historical rain cluster concentration structure index Gjz under the same karst landform type, the same monitoring area, and the same assessment scale was collected. The mean and standard deviation were calculated using statistical methods, and the sum of the mean and standard deviation was set as the rain cluster concentration judgment threshold Yg. The real-time rain cluster concentration structure index Gjz was then compared with the rain cluster concentration judgment threshold Yg. When the rain cluster concentration structure index Gjz ≤ the rain cluster concentration judgment threshold Yg, it was determined that the rainfall of the current rainy day was distributed evenly in time. When the rain cluster concentration structure index Gjz > the rain cluster concentration judgment threshold Yg, it was determined that the rainfall of the current rainy day was concentrated in time. The key point of this step is to consolidate the previously scattered records of the total number of rain peak jumps N, erosion duration T, maximum continuous 60-minute rainfall H, and daily rainfall P into a rain cluster concentration structure index Gjz that can be directly used in subsequent analysis. A unified judgment caliber is formed by using the rain cluster concentration judgment threshold Yg, making it comparable between different erosive rainfall days. Compared with the current technical means that mostly judge based on cumulative rainfall or single short-duration rainfall intensity, this step incorporates rainfall jumps, duration, and local high-value segments into the same judgment chain, distinguishing rainfall processes with "similar total amount but different temporal organization", so that the rain cluster distribution judgment is closer to the actual temporal characteristics of erosive rainfall in karst mountainous areas. Example

[0030] Please refer to Figure 3 Specifically: S4 includes S41; S41. When the distribution of rain clusters is determined to be concentrated in the time of rainfall, construct the soil-rock discontinuity coupling index Kdo and the sinking-infiltration switching index Dsq based on static surface data and daily data of erosive rainfall, respectively. The soil-rock discontinuous coupling index Kdo is used to reflect the degree of participation of karst surface structure in rainfall erosion pathways, specifically: Among them, Kdo u Ym u Lz u 、Lm u and Qj uThese represent the soil-rock discontinuity coupling index, bare rock cutting density, median soil layer continuity length, fracture density, and average distance from the fracture to the dominant scour line for the u-th evaluation unit, respectively. The infiltration-sinking switching index Dsq is used to reflect the distribution of water along the soil-rock interface, fissures, and depressions during the current rainfall day, specifically: , among which, Dsq j U j and S j , representing the infiltration-drainage switching index, interface ingress depth, and depression retention depth, respectively, on the j-th day of erosive rainfall.

[0031] S4 further includes S42 and S43; S42. The rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo, and the sink-infiltration switching index Dsq are fitted and calculated to construct the comprehensive erosion driving value Eqd, which represents the daily comprehensive erosion driving state of the assessment unit on each erosive rainfall day. Specifically: Where k represents the total number of monitoring days, Eqd u，j This represents the composite erosion driving value of the u-th assessment unit on the j-th day of erosive rainfall; S43. Collect a sample set of comprehensive erosion driving values ​​Eqd under the same historical karst landform type, monitoring area and assessment scale, and calculate the mean and standard deviation using statistical methods. Set the sum of the mean and standard deviation as the erosion judgment threshold Ye, and then perform comprehensive erosion assessment with the real-time acquired comprehensive erosion driving values ​​Eqd, as follows; When the comprehensive erosion driving value Eqd ≤ erosion judgment threshold Ye, it means that the surface structure of the current assessment unit is stable under the influence of rainfall erosion force, and normal monitoring should be maintained at this time. When the comprehensive erosion driving value Eqd > the erosion judgment threshold Ye, it indicates that the rainfall erosion force has an impact on the surface structure of the current assessment unit. At this time, static surface data when the rain cluster distribution is determined to be that the rainfall is distributed evenly in time are collected to perform a second iteration assessment. If the second iteration assessment shows that there is stability, it indicates that the rainfall erosion force is in the erosion accumulation state. If the second iteration assessment still shows an impact, it indicates that the rainfall erosion force is in the erosion destruction state.

[0032] In this embodiment, the rain cluster distribution determination result is used as the entry point. When the rain cluster distribution is determined to be concentrated in terms of rainfall over time, static surface data and data of erosive rainfall days are retrieved. Within the scale of the assessment unit, the soil-rock discontinuity coupling index Kdo and the infiltration-sinking switching index Dsq are constructed respectively. The soil-rock discontinuity coupling index Kdo is used to characterize the participation of bare rock cutting, soil layer continuity, fissure distribution and their relationship with the dominant scour line in the u-th assessment unit on the rainfall erosion path. The infiltration-sinking switching index Dsq is used to represent the distribution state of water body in the j-th erosive rainfall day between the soil-rock interface transfer, fissure inflow and depression retention. Then, the rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo and the infiltration-sinking switching index Dsq are fitted and calculated to obtain the comprehensive erosion driving value Eqd.

[0033] For each assessment unit, the soil-rock discontinuity coupling index Kdo is calculated based on static surface data. This formula decomposes "bare rock cutting - soil layer continuity" and "fracture density - fracture location" into two product terms and then takes the square root. The former corresponds to the degree to which slope runoff is reorganized at the soil-rock interface, and the latter corresponds to the degree to which water flow is pulled by fractures and deviates from the dominant scour line. The product form indicates that these two types of structural conditions must be considered together. The square root form is used to compress the amplification effect of extreme values ​​in the comprehensive calculation, and the "1+" before each term is used to maintain the dimensionless expression and avoid the loss of continuity of the integral formula when a term is 0. Subsequently, for erosive rainfall days, the infiltration-drainage switching index Dsq is calculated based on the data of erosive rainfall days. This formula uses the proportion of the depression retention water depth S in the interface transition water depth U to characterize "drainage dominance". The pattern remains "retention-driven." When the proportion of the interface transition water depth U is relatively large, it indicates that the water body is more inclined to be transported downstream along the interface and fissures, and the infiltration-capacity switching index Dsq is correspondingly higher. When the proportion of the depression retention water depth S is relatively large, it indicates that the water body stays more in the depression retention area, and the infiltration-capacity switching index Dsq is closer to 1. Finally, the rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo, and the infiltration-capacity switching index Dsq of the erosive rainfall days are fitted and calculated to obtain the comprehensive erosion driving value Eqd. The concentrated rainfall process within the total monitoring days is accumulated to obtain the comprehensive erosion driving expression of the assessment unit within the monitoring period. For example, in the u-th assessment unit between the toe of a peak cluster slope and a sinkhole, if Ym=2.4, Lz=1.8, Lz=0.5, Lm=1.1, and Qj=2.0, then... ≈7.17 indicates that the cutting of bare rock, the distribution of soil layers, and the location of fractures within the current assessment unit have a strong influence on the erosion path; if the measured values ​​for a certain day of corrosive rainfall are U=12 mm and S=8 mm... If the rainfall is mm, then Dsq = 1 + 12 / (12 + 8) = 1.60, indicating that water is shifting more along the interface and fissures on the current day of erosive rainfall. If Gjz = 1.6 on the same day, then the comprehensive erosion driving value Eqd corresponding to the same day's item on the current day of erosive rainfall is 1.6 × 7.17 × 1.60 ≈ 18.36. The purpose of this process is to put the surface structure conditions, the concentration of rainfall time, and the distribution of water between the interface transition and the depression retention state into the same calculation link, so that the soil-rock discontinuity coupling index Kdo is responsible for describing "how the surface guides the erosion path", the infiltration-confluence switching index Dsq is responsible for describing "how the water body switches between confluence and retention on the same day", the rain cluster concentration structure index Gjz is responsible for describing "how rainfall is concentrated in time", and the comprehensive erosion driving value Eqd combines the three into a single interpretable quantity that can be directly used for comprehensive erosion assessment. A sample set of comprehensive erosion driving values ​​Eqd under the same historical karst landform type, monitoring area, and assessment scale was collected. The mean and standard deviation were calculated using statistical methods, and the sum of the mean and standard deviation was set as the erosion judgment threshold Ye. The real-time acquired comprehensive erosion driving values ​​Eqd were then compared with the erosion judgment threshold Ye: when the comprehensive erosion driving value Eqd ≤ the erosion judgment threshold Ye, the surface structure of the current assessment unit was determined to be stable under the influence of rainfall erosion and maintained normal monitoring; when the comprehensive erosion driving value Eqd > the erosion judgment threshold Ye, the surface structure of the current assessment unit was determined to have been affected by rainfall erosion. Furthermore, static surface data when the rain cluster distribution was determined to be a time-slow distribution of rainfall were collected. A second-stage iterative assessment is conducted. If the second-stage assessment indicates stability, the state is determined to be one of erosion accumulation. If the second-stage assessment still indicates an impact, the state is determined to be one of erosion destruction. The purpose of this step is to integrate the rainfall concentration index Gjz (reflecting the degree of rainfall time concentration), the soil-rock discontinuity coupling index Kdo (reflecting the soil-rock discontinuity coupling relationship), and the runoff-infiltration switching index Dsq (reflecting the runoff-infiltration switching state) into the comprehensive erosion driving value Eqd's comprehensive judgment chain. Compared with the current technical means of making judgments based solely on rainfall amount, rainfall intensity, or single slope topographic parameters, this step incorporates the characteristics of soil-rock interweaving, fissure development, depression retention, and interface transfer in karst mountainous areas into the comprehensive erosion assessment. Example

[0034] Please refer to Figure 2 A rainfall erosivity assessment system for karst mountainous areas includes a surface modeling module, a hydrological data acquisition module, a structure determination module, and an erosion assessment module. The surface modeling module is used to construct a surface point cloud based on UAV orthophotos, select the lowest point according to the grid to generate ground candidate points, and determine the ground points based on the elevation difference and slope to generate a digital elevation model. The digital elevation model is divided into evaluation units and static surface data is output. The hydrological acquisition module is used to collect hydrological time-series data in real time based on the monitoring devices deployed in each assessment unit, and to statistically analyze daily data on erosive rainfall. The structure determination module is used to construct the rain cluster concentration structure index Gjz based on the daily data of erosive rainfall, and to determine the distribution of rain clusters; The erosion assessment module is used to construct the soil-rock discontinuity coupling index Kdo and the sink-infiltration switching index Dsq when the rain cluster distribution is determined to be concentrated in time. These are then fitted with the rain cluster concentration structure index Gjz to calculate and construct the comprehensive erosion driving value Eqd for comprehensive erosion assessment.

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

Claims

1. A method for assessing rainfall erosivity in karst mountainous areas, characterized in that: Includes the following steps: S1. Construct a surface point cloud based on UAV orthophotos, select the lowest point according to the grid to generate ground candidate points, and generate a digital elevation model based on the elevation difference and slope of the ground points. Divide the digital elevation model into evaluation units and output static surface data. S2. Collect hydrological time-series data in real time based on the monitoring devices deployed in each assessment unit, and compile daily data on erosive rainfall. S3. Construct the rain cluster concentration structure index Gjz based on daily data of erosive rainfall, and determine the distribution of rain clusters; Among them, Gjz j N j T j and H j P represents the rain cluster concentration structure index, total number of rain peak jumps, erosion duration, and maximum continuous 60-minute rainfall on the j-th erosive rainfall day. j T0 represents the daily rainfall on the j-th day of erosive rainfall, and T0 represents the time base, which is the standard time scale used when calculating the erosion duration and maximum continuous rainfall of the rainfall process data. S4. When the rain cluster distribution is determined to be concentrated in the time distribution of rainfall, the soil-rock discontinuous coupling index Kdo and the sinking-infiltration switching index Dsq are constructed respectively, and fitted with the rain cluster concentrated structure index Gjz to calculate the comprehensive erosion driving value Eqd for comprehensive erosion assessment. Among them, Kdo u Ym u Lz u 、Lm u and Qj u These represent the soil-rock discontinuity coupling index, bare rock cutting density, median soil layer continuity length, fracture density, and average distance from the fracture to the dominant scour line for the u-th evaluation unit, respectively. Among them, Dsq j U j and S j These represent the infiltration-drainage switching index, the interface ingress depth, and the depression retention depth, respectively, on the j-th day of erosive rainfall. Where k represents the total number of monitoring days, Eqd u，j This represents the composite erosion driving value of the u-th assessment unit on the j-th day of erosive rainfall.

2. The method for assessing rainfall erosivity in karst mountainous areas according to claim 1, characterized in that: S1 includes S11; S11. Based on the orthophoto data of the karst mountain area collected by UAV, perform aerial triangulation on the orthophoto data to obtain photogrammetric exterior orientation elements, and perform multi-view dense matching based on the photogrammetric exterior orientation elements to generate surface point cloud data. Divide the surface point cloud data into grids according to the set grid size, and select the lowest elevation point in each grid as ground candidate points. Use the ground candidate points as seed points to calculate the elevation difference and local slope between the remaining elevation points and the adjacent seed points. When the elevation difference is not greater than the preset elevation difference threshold and the local slope is not greater than the preset slope threshold, the current elevation point is determined to be a ground point; otherwise, it is determined to be a non-ground point. Construct an irregular ground triangular network based on all ground points, and perform raster interpolation on the irregular ground triangular network to generate a digital elevation model. The elevation difference threshold and slope threshold are set by using the lowest elevation point in each grid as the ground candidate point, establishing the candidate point neighborhood based on the relationship between adjacent grids, calculating the elevation difference sequence and slope sequence between adjacent ground candidate points, and setting the elevation difference threshold by calculating the sum of the median value and interquartile range of the elevation difference sequence using statistical methods, and setting the slope threshold by calculating the sum of the median value and interquartile range of the slope sequence.

3. The method for assessing rainfall erosivity in karst mountainous areas according to claim 2, characterized in that: S1 further includes S12; S12. Extract the slope aspect, slope length and flow direction of each grid in the digital elevation model, trace the confluence path according to the flow direction, determine the confluence path with the largest flow volume as the dominant scour line, and divide the target area of ​​the same dominant scour line to obtain several evaluation units. Based on orthophoto data, the boundaries of bare rock, soil layers, and fissures in karst mountain areas are identified. The number of fissures (Ls) within each assessment unit is counted, and the total length (Lr) of the boundary between the bare rock and soil layers and the area (Mj) of each assessment unit are calculated. Then, the bare rock cutting density (Ym) and fissure density (Lm) are calculated separately. , , among which, Ym u Lr u Mj u 、Ls u and Lm u These represent the bare rock cutting density, total boundary length, area, number of fracture openings, and fracture opening density of the u-th evaluation unit, respectively. Sampling lines were laid out along the main slope, and the lengths of the soil strips intersecting with the soil layer area of ​​each sampling line were counted to obtain a soil strip length sequence. The median value of the soil layer length Lz was obtained by extracting the median of the soil strip length sequence according to the statistical method. The shortest distance Zd from each fracture opening to the dominant scour line was calculated using statistical methods to obtain a distance sequence. The average distance Qj between the fracture opening and the dominant scour line in the distance sequence was then calculated, specifically as follows: , among which, Qj u Zd represents the average distance from the fracture opening to the dominant scour line in the u-th evaluation unit. u，i This represents the shortest distance from the i-th fracture opening in the u-th evaluation unit to the dominant scour line; Based on the bare rock cutting density Ym, the median value of soil layer continuity length Lz, the fracture opening density Lm, and the average distance from the fracture opening to the dominant scour line Qj, output the static surface data of each evaluation unit.

4. The method for assessing rainfall erosivity in karst mountainous areas according to claim 3, characterized in that: S2 includes S21; S21. Based on the monitoring devices set up in the corresponding karst mountain area of ​​each assessment unit, the hydrological time series data are collected in real time and written into the data table of the corresponding assessment unit according to the location of the monitoring devices. The monitoring devices include a rainfall monitoring device, an interface runoff monitoring device, a fissure runoff monitoring device, and a water level monitoring device at the front edge of a depression or sinkhole. Among them, the rainfall monitoring device is used to collect rainfall sequences, the interface runoff monitoring device is set at the junction of the bare rock area and the soil layer area to collect the runoff volume entering the soil layer area along the junction of the bare rock area and the soil layer area, the fissure mouth runoff monitoring device is set at the upstream confluence of the fissure mouth to collect the runoff volume entering the fissure mouth, and the water level monitoring device is set at the water collection position at the front edge of the depression or sinkhole to collect the water level.

5. The method for assessing rainfall erosivity in karst mountainous areas according to claim 4, characterized in that: S2 also includes S22 and S23; S22. Divide the rainfall sequence into natural days, calculate the cumulative rainfall for each natural day, obtain the daily rainfall P for each natural day, and perform cumulative calculation on the rainfall sequence within each natural day with a calculation duration of 60 minutes and a sliding step size of 10 minutes. Determine the maximum value among all 60-minute cumulative rainfall as the maximum continuous 60-minute rainfall H for the natural day. Within each natural day, the period with a rainfall of not less than 0.5 mm in 10 minutes is defined as the effective rainfall period, and the sum of the durations of all effective rainfall periods is defined as the erosion duration T of the natural day. Within each natural day, a change in rainfall intensity of not less than 6 mm / h between two adjacent 10-minute intervals is recorded as a rainfall peak jump, and the total number of rainfall peak jumps N within each natural day is counted; the flow volume recorded by the interface flow monitoring device and the flow volume recorded by the fissure flow monitoring device are combined and divided by the flow area within the corresponding assessment unit to obtain the interface transition water depth U for the current natural day; the maximum water level value among the accumulated water levels recorded by the water level monitoring device is determined as the depression retention water depth S for the current natural day; S23. Rainfall days that meet the condition of daily rainfall P≥30mm are defined as erosive rainfall days. The daily rainfall P, maximum continuous 60-minute rainfall H, erosion duration T, total number of rain peak jumps N, depression water depth S, and interface transition water depth U of the erosive rainfall days are written into the corresponding evaluation unit to form erosive rainfall day data.

6. The method for assessing rainfall erosivity in karst mountainous areas according to claim 5, characterized in that: S3 includes S31; S31. Construct a rain cluster concentration structure index Gjz based on erosive rainfall day data to reflect the aggregation state of rainfall over time on erosive rainfall days.

7. The method for assessing rainfall erosivity in karst mountainous areas according to claim 6, characterized in that: S3 further includes S32; S32. Collect a sample set of rain cluster concentration structure index Gjz under the same historical karst mountain landform type, monitoring area and assessment scale, and calculate the mean and standard deviation by statistical method. Set the sum of the mean and standard deviation as the rain cluster concentration judgment threshold Yg, and then use it to judge the rain cluster distribution with the real-time acquired rain cluster concentration structure index Gjz, as follows. When the rain cluster concentration structure index Gjz ≤ the rain cluster concentration judgment threshold Yg, it indicates that the rainfall on the current rainy day is distributed evenly over time. When the rain cluster concentration structure index Gjz > the rain cluster concentration judgment threshold Yg, it indicates that the rainfall on the current rainy day is concentrated in time.

8. The method for assessing rainfall erosivity in karst mountainous areas according to claim 7, characterized in that: S4 includes S41; S41. When the distribution of rain clusters is determined to be concentrated in the time of rainfall, construct the soil-rock discontinuity coupling index Kdo and the sinking-infiltration switching index Dsq based on static surface data and daily data of erosive rainfall, respectively. The soil-rock discontinuous coupling index Kdo is used to reflect the degree of participation of karst surface structure in rainfall erosion paths. The infiltration-sinking switching index Dsq is used to reflect the distribution of water along the soil-rock interface, fissures, and depressions during the current rainfall day.

9. A method for assessing rainfall erosivity in karst mountainous areas according to claim 8, characterized in that: S4 further includes S42 and S43; S42. The rain cluster concentration structure index Gjz, the soil-rock discontinuity coupling index Kdo, and the sink-infiltration switching index Dsq are fitted and calculated to construct the comprehensive erosion driving value Eqd, which represents the single-day comprehensive erosion driving state of the evaluation unit on each erosive rainfall day. S43. Collect a sample set of comprehensive erosion driving values ​​Eqd under the same historical karst landform type, monitoring area and assessment scale, and calculate the mean and standard deviation using statistical methods. Set the sum of the mean and standard deviation as the erosion judgment threshold Ye, and then perform comprehensive erosion assessment with the real-time acquired comprehensive erosion driving values ​​Eqd, as follows; When the comprehensive erosion driving value Eqd ≤ erosion judgment threshold Ye, it means that the surface structure of the current assessment unit is stable under the influence of rainfall erosion force, and normal monitoring should be maintained at this time. When the comprehensive erosion driving value Eqd > the erosion judgment threshold Ye, it indicates that the rainfall erosion force has an impact on the surface structure of the current assessment unit. At this time, static surface data when the rain cluster distribution is determined to be that the rainfall is distributed evenly in time are collected to perform a second iteration assessment. If the second iteration assessment shows that there is stability, it indicates that the rainfall erosion force is in the erosion accumulation state. If the second iteration assessment still shows an impact, it indicates that the rainfall erosion force is in the erosion destruction state.

10. A rainfall erosivity assessment system for karst mountainous areas, applied to the rainfall erosivity assessment method for karst mountainous areas as described in any one of claims 1-9, characterized in that: It includes a surface modeling module, a hydrological data acquisition module, a structure determination module, and an erosion assessment module; The surface modeling module is used to construct a surface point cloud based on UAV orthophotos, select the lowest point according to the grid to generate ground candidate points, and determine the ground points based on the elevation difference and slope to generate a digital elevation model. The digital elevation model is divided into evaluation units and static surface data is output. The hydrological acquisition module is used to collect hydrological time-series data in real time based on the monitoring devices deployed in each assessment unit, and to statistically analyze daily data on erosive rainfall. The structure determination module is used to construct the rain cluster concentration structure index Gjz based on the daily data of erosive rainfall, and to determine the distribution of rain clusters; The erosion assessment module is used to construct the soil-rock discontinuity coupling index Kdo and the sink-infiltration switching index Dsq when the rain cluster distribution is determined to be concentrated in time. These are then fitted with the rain cluster concentration structure index Gjz to calculate and construct the comprehensive erosion driving value Eqd for comprehensive erosion assessment.

Images