A landslide risk prediction method based on tectonic stress factors suitable for the Sichuan-Tibet region

Abstract

This invention discloses a landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region, belonging to the field of geological disaster prevention and control technology. It aims to address the insufficient accuracy of traditional landslide prediction methods under strong tectonic stress and complex terrain in the Sichuan-Tibet region. The method includes: collecting data and processing it into DEM raster data; converting the DEM raster data into point-like Shapefile data; extracting the planar coordinates and elevation values ​​of topographic points to form a topographic point matrix with associated global coordinates; calculating the slope length and gradient of the topographic point and the landslide center point, as well as the tectonic stress and topographic stress of the topographic point, to construct a tectonic stress magnitude hazard factor; determining the direction of the regional stress field in molecular regions and the angle between the topographic point and the near-far fault zone, to construct a tectonic stress direction hazard factor; calculating the landslide risk of each topographic point, dividing the area into extremely high to extremely low hazard zones, and outputting a hazard map of the calculated area. This invention is adapted to the geological characteristics of the Sichuan-Tibet region, improves prediction accuracy, and provides a scientific basis for landslide prevention and control in the Sichuan-Tibet region.

Description

Technical Field

[0001] This invention relates to the field of geological disaster prevention and control technology, specifically to a landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region. Background Technology

[0002] Landslides are a common geological hazard in the southeastern edge of the Qinghai-Tibet Plateau. Due to the steep terrain, complex geological structure, concentrated rainfall, and frequent earthquakes in the southeastern edge of the Qinghai-Tibet Plateau, landslides pose a serious threat to the lives and property of residents, infrastructure, and the ecological environment. Therefore, accurately predicting the probability and location of landslides is crucial for disaster prevention and regional planning.

[0003] Currently, landslide prediction mainly employs statistical analysis, physical models, remote sensing and GIS technologies, and integrated modeling methods. While most models incorporate multiple factors, they do not fully consider the actual tectonic stress, topographic slope length, and the influence of topographic stress. For areas with complex geological structures and strong topographic stress on the southeastern edge of the Qinghai-Tibet Plateau, some existing technologies have attempted to couple topography and tectonic stress to improve prediction accuracy. For example, patent CN118568965B proposes a method for estimating the three-dimensional stress state of large-scale slopes in complex terrain areas. By superimposing tectonic stress and topographic stress, a high-precision three-dimensional stress field is generated, providing a stress basis for landslide prediction. However, this method still has significant limitations: First, it fails to solve the problem of accurate slope length calculation under the complex terrain of the Sichuan-Tibet region, while slope length directly affects the quantitative accuracy of topographic stress and the assessment of landslide dynamic conditions; second, it fails to establish a correlation between slope length, tectonic stress, and landslide hazard factors, making it impossible to effectively transform stress calculation results into a basis for landslide hazard classification; third, it lacks adaptability to the coupled effects of multiple factors such as complex topographic stress disturbance, rainfall, and earthquakes in the Sichuan-Tibet region, making it difficult to accurately reflect the triggering mechanism of landslides in this area, leading to deviations between the predicted results and the actual landslide distribution. Therefore, there is an urgent need for a landslide prediction method based on the characteristics of the regional high mountain and canyon terrain and the influence of strong tectonic stress fields to improve the accuracy and practicality of landslide prediction in the Sichuan-Tibet region and provide support for landslide disaster prediction and prevention. Summary of the Invention

[0004] The purpose of this invention is to provide a landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region, so as to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region includes the following steps:

[0007] S1. Collect DEM topographic data, regional tectonic stress field data, GPS motion data and fault distribution data within the calculation area, and process them into DEM raster data using geographic information system tools;

[0008] S2. Using the raster-to-point tool in ArcGIS software, the DEM raster data is converted into point-like Shapefile data. The plane coordinates and elevation values ​​of each terrain point are extracted to form a terrain point array covering the calculation area. The terrain point array consists of several terrain points, and each terrain point is associated with its corresponding global coordinates.

[0009] S3. Calculate the slope length, slope, tectonic stress, and topographic stress of each topographic point in the topographic point matrix, and at the same time calculate the slope length and slope of each landslide center point in the historical landslide database;

[0010] S4. Statistically analyze the distribution patterns of slope length and slope gradient of all landslide center points in the topographic lattice, and construct a structural stress magnitude hazard factor, which includes a slope length hazard factor and a slope gradient hazard factor.

[0011] S5. Divide the calculation area into several sub-regions and determine the regional stress field direction of each sub-region. At the same time, divide each sub-region into near-fault region and far-fault region, and obtain the included angle of each topographic point accordingly.

[0012] S6. Statistically analyze the distribution pattern of the included angles of all topographic points in the topographic point matrix and the landslide direction recorded in the historical landslide database, and construct the tectonic stress direction hazard factor.

[0013] S7. Based on the hazard factors of tectonic stress magnitude and tectonic stress direction, calculate the landslide risk of each topographic point, and divide the calculation area into extremely high, high, medium, low and extremely low hazard zones according to the landslide risk value, and output the hazard map of the calculation area.

[0014] Preferably, the global coordinates are defined based on the OXYZ right-handed coordinate system, where the X-axis points east, the Y-axis points north, and the Z-axis is vertically upward with the average level of the calculation area as the zero point of the Z-axis; the global coordinates are represented as (x, y, z), where x and y are the planar coordinates of the terrain point, and z is the elevation value of the terrain point.

[0015] Preferably, the slope calculation process is as follows: The DEM raster data is processed using ArcGIS Pro's slope analysis tool to generate a slope raster map. After converting the slope raster map into vector point data, it is matched with the terrain point matrix using a spatial connection tool to obtain the slope 'a' of each terrain point. Ai Topographic points are represented by Ai.

[0016] Preferably, the specific steps of the method for calculating the slope length are as follows:

[0017] S31. Starting from a topographic point, search for topographic points in the topographic point array within a 10km radius around the starting point in the east, south, west and north directions with a step size of 100m. Search for topographic points in the topographic point array step by step in each direction, record the global coordinates of the topographic points searched in each direction, and form a sequence of elevation points in each direction.

[0018] S32. For each elevation point sequence in each direction, the priority order from high to low is compression effect judgment, gully isolation judgment, flat terrain judgment and terrain trend judgment. The terrain points in the elevation point sequence are dynamically terminated one by one. If any judgment condition is met, the search for terrain points in the elevation point sequence in that direction is stopped, and the terrain point searched at this time is determined as the direction endpoint of that direction.

[0019] The rule for determining the compression effect is as follows: calculate the slope 'a' of the line connecting the searched terrain point within the current elevation point sequence to the starting point. jk The terrain points found within the current elevation point sequence are denoted as Ak, and the starting point is denoted as Aj. The calculation method is as follows: Among them, z k and z Aj These are the elevation values ​​of the terrain points found within the current elevation point sequence and the starting point, respectively, x. k and y k It is the planar coordinate of the terrain point searched within the current elevation point sequence, x Aj and y Aj These are the planar coordinates of the starting point;

[0020] If a jk If the value is greater than or equal to the preset compression threshold S, then it is determined that Ak exerts a compression effect along the slope on the starting point Aj, indicating that it is still within the effective extension segment of the slope. It is necessary to continue searching for the next topographic point in the elevation point sequence and re-determine the compression effect; if a jk If the preset compression threshold S is less than the threshold value, then trench isolation is determined.

[0021] The rule for determining gully isolation is as follows: compare the elevation value of Ak with the elevation value of the terrain point A(k-1) searched in the previous elevation point sequence. If z k <z k-1 - If a pre-set gully threshold X is set, it is determined that there is a gully terrain isolation between Ak and the terrain point A(k-1) searched in the previous elevation point sequence. The search must be stopped immediately, and A(k-1) is determined as the endpoint of the direction; if z k ≥z k-1 - If a gully threshold X is preset, then the gully isolation situation is excluded, and then the flat terrain judgment and terrain trend judgment are entered;

[0022] The rule for determining flat terrain is: if Ak satisfies if a jk <Preset compression threshold S and z k ≥z k-1 - If a pre-set gully threshold X is defined, i.e., non-compression and non-gully terrain, then the area is considered to be in a flat terrain segment. Ak is taken as the starting point of the flat segment, and the search continues from this starting point to the next terrain point within the elevation point sequence, accumulating the horizontal distance. The calculation method for the accumulated horizontal distance is as follows: Where D is the cumulative horizontal distance, D prev It is the cumulative horizontal distance of the previous search, and its initial value is 0, x k-1 and y k-1 These are the planar coordinates of the terrain points searched within the previous elevation point sequence;

[0023] If D > 700m, the flat terrain segment is determined to have exceeded the effective extension range of the slope, the search is stopped, and the starting point of the flat segment is determined as the endpoint of the direction; if D ≤ 700m, but a terrain point Am (m > k) is found midway through the search for the next terrain point in the elevation point sequence starting from the starting point of the flat segment, satisfying a km If the value is greater than or equal to the preset compression threshold S, then the user re-enters the slope section with compression effect, stops the search, and determines Am as the directional endpoint in that direction, where a km Let a be the slope of the line connecting Am and Ak. jk ;

[0024] The rule for determining the terrain trend is: if Ak and the terrain point A(k+1) searched in the next elevation point sequence satisfy z twice consecutively... k <z (k-1) And z k+1 <z k If the elevation value decreases twice in a row, it is determined that the terrain in that direction has entered a downhill section or is close to the top of the slope. The search should be stopped and the terrain point found in the previous elevation point sequence should be determined as the endpoint of that direction.

[0025] S33. For the endpoints of each direction (x) dir y dir , z dir ), calculate the elevation value z of the endpoint in each direction respectively. dir Elevation value z of the starting point Aj The difference is used to obtain the elevation difference Δz in each direction. dir , where Δz dir A positive value indicates that the endpoint of the direction is higher than the starting point, suggesting that the terrain along that direction shows an upward trend. Δz dir A negative value indicates that the endpoint of a direction is lower than the starting point, suggesting a downward trend in the terrain along that direction. Simultaneously, based on the global coordinates of the endpoints and starting points in each direction, the slope vector V corresponding to the elevation difference in each direction is calculated.dir =(Δx) dir Δy dir Δz dir ), where Δx dir =x dir -x Aj Δy dir =y dir -y Aj ;

[0026] S34. Filtering for elevation differences in each direction that satisfy Δz dir The slope vectors corresponding to slopes greater than 0 are denoted as the set of filtered slope vectors {V1, V2, ..., Vm} (m≤4, where m is a slope vector with a positive elevation difference). The set of slope vectors is then combined according to the slope vector composition rules to form a slope length vector, and the magnitude of this slope length vector is the slope length L. Ai ;

[0027] The slope vector synthesis rule is as follows:

[0028] When m = 1, the slope vector in that direction is directly used as the slope length vector V. slope That is, the slope length vector is V1; at this time, it indicates that the topographic point extends to higher terrain in only one direction, and this direction vector can fully reflect the slope extension characteristics.

[0029] When m = 2, if the slope vectors in two directions are adjacent slope vectors (i.e., east and north, south and west, east and south, and north and west), then the vector components corresponding to the slope vectors in the two directions are added together or their maximum values ​​are taken, i.e., Δx slope =Δx1 + Δx2, Δy slope =Δy1+Δy2,Δz slope =max(Δz1, Δz2), to obtain the slope length vector V. slope =(Δx) slope Δy slope Δz slope If the slope vectors of the two directions are slope vectors of opposite directions, i.e. east and west and south and north, choose the larger slope vector between Δz1 and Δz2 as the slope length vector. That is, if Δz1>Δz2, the slope length vector is V1, and otherwise the slope length vector is V2.

[0030] When m = 3, the vector components corresponding to the slope vectors in the three directions are added together, i.e., Δx slope =Δx1 + Δx2 + Δx3, Δy slope =Δy1+Δy2+Δy3,Δz slope =Δz1+Δz2+Δz3, thus obtaining the slope length vector V slope =(Δx) slope Δy slopeΔz slope );

[0031] When m = 4, select Δz from the slope vectors in the four directions. dir The slope vector with the largest value is taken as the slope length vector.

[0032] Preferably, the calculation process for the terrain stress is as follows:

[0033] Obtain the slope length vector V of the terrain point slope The elevation difference term Δz slope The elevation difference term Δz slope and the slope a of the terrain point Ai Substitute into the topographic stress calculation formula to calculate the topographic stress S topo-Ai ;

[0034] The formula for calculating terrain stress is:

[0035] S topo-Ai =–ρgh L sina Ai cosa Ai ;

[0036] Among them, h L ρ is the elevation difference term in the slope length vector of the topographic point, ρ is the average density of rock in the calculation area, and g is the gravitational acceleration.

[0037] Preferably, the process for constructing slope length hazard factors and slope gradient hazard factors is as follows:

[0038] Records of landslide events that have occurred in the past 30 years within the calculation area were selected from the historical landslide database. Each landslide event required the coordinates of the landslide center point, verified by on-site investigation or high-resolution remote sensing imagery. For each landslide center point, the slope length L was obtained using a calculation method that was completely consistent with the slope length and gradient of the topographic points in the topographic matrix. i (i = 1, 2, ..., n, where n is the total number of landslide events) and slope a i This leads to the formation of the slope length point set L. landslide ={L1,L2,…,L} and the slope point set a landslide ={a1,a2,…,a};

[0039] Based on the actual distribution of slope extension lengths in the high mountain and canyon terrain of the Sichuan-Tibet region, the slope lengths of the terrain lattice are divided into six sub-intervals: [0, 500) m, [500, 1000) m, [1000, 1500) m, and [1500, +∞) m, each denoted as L. o (o=1,2,…,4), count the slope length point set L within each slope length subinterval. landslide Number of slope lengths included That is, belonging to Lo L i The number of slope lengths is counted, and the total number of slope lengths in the terrain lattice within the slope length sub-interval is also counted. That is, the slope length in the terrain lattice belongs to L o The number of terrain points is determined by the formula. Calculate the percentage of slope length in each slope length sub-interval. This is used to reflect the slope length probability density of landslides occurring within a slope length interval. Based on the high-incidence slope range of landslides in the Sichuan-Tibet region, the slope of the topographic lattice is divided into 5 slope sub-intervals: [0, 15°), [15°, 30°), [30°, 45°), [45°, 60°), and [60°, 90°]. Each slope sub-interval is denoted as a. m (m=1,2,…,5), the slope percentage of each slope sub-interval is calculated using the same logic as the calculation of the slope length percentage of each slope length sub-interval.

[0040] Using the minimum-maximum normalization method, the slope length ratio is... and slope ratio Mapping to the (0,1) interval yields the normalized probability distributions of slope length and slope gradient, and based on these, a slope length risk factor f is constructed. L and slope risk factor f a Both the slope length hazard factor and the slope gradient hazard factor are expressed in piecewise function form;

[0041] The slope length hazard factor is:

[0042]

[0043] Among them, val1, val2, val3 and val4 are the normalized slope length probability distributions;

[0044] The slope hazard factor is:

[0045]

[0046] Among them, vaa1, vaa2, vaa3, vaa4, and vaa5 are normalized slope probability distributions.

[0047] Preferably, the process of determining the regional stress field direction of each sub-region is as follows: statistically analyze the average value of the focal mechanism solution direction, the average value of the measured stress direction, and the average value of the GPS displacement direction for each sub-region. If the directional deviation between the average value of the focal mechanism solution direction, the average value of the measured stress direction, and the average value of the GPS displacement direction is less than 30°, then the arithmetic mean of the average value of the focal mechanism solution direction, the average value of the measured stress direction, and the average value of the GPS displacement direction is taken as the regional stress field direction of the sub-region; otherwise, the average value of the GPS displacement direction is taken as the regional stress field direction of the sub-region.

[0048] Preferably, the historical landslide database is used to record the landslide center point and the landslide direction at that center point for each landslide event.

[0049] Preferably, the near-fault region is a sub-region with a horizontal distance of less than or equal to 10 km from the fault, and the far-fault region is a sub-region with a horizontal distance of more than 10 km from the fault, and the included angle θ of each topographic point in the near-fault region is defined as follows: Ai This includes the angle θ1 between the slope aspect of the topographic point and the direction perpendicular to the fault strike, and the angle θ2 between the slope aspect of the topographic point and the direction of the regional stress field, as well as the angle θ for each topographic point within the far-fault region. Ai θ2 is the angle between the slope aspect of the topographic point and the direction of the regional stress field; the slope aspect of the topographic point is obtained by processing DEM raster data using ArcGIS Pro software, and the fault strike is obtained based on fault distribution data.

[0050] Preferably, the process for constructing slope length hazard factors and slope gradient hazard factors is as follows:

[0051] Landslide events within the computational area are extracted from the historical landslide database, and the landslide direction S at the landslide center point corresponding to each landslide event is obtained. i (i = 1, 2, ..., n, where n is the total number of landslide events), the landslide slides towards S i Including θ1 and θ2, thus forming the landslide sliding set S landslide ={S1,S2,…,S}, referring to the statistical patterns of near-fault landslide angles in the Sichuan-Tibet region, θ1 is divided into three first-angle sub-intervals: A1[0,30°), A2[30°,60°), and A3[30°,90°), and θ2 is divided into three second-angle sub-intervals: B1[0,30°), B2[30°,60°), and B3[30°,90°), forming nine combined angle sub-intervals, namely A1B1, A1B2, A1B3, A2B1, A2B2, A2B3, A3B2, A3B2, and A3B3. Each combined angle sub-interval is denoted as L. AiBj(i = 1, 2, 3; j = 1, 2, 3), and referring to the statistical law of the included angle of the distant fault landslide in the Sichuan-Tibet region, θ1 is re-divided into two third included angle sub-intervals: C1[0, 45°) and C2[45°, 90°), and each third included angle sub-interval is denoted as L. Ci (i = 1, 2, 3);

[0052] Statistical analysis of the landslide direction set S within each combined angle sub-interval landslide The number of landslides included That is, belonging to L AiBj S i The number of points is counted, and the number of angles of topographic points belonging to the near-fault region in the topographic lattice of the combined angle sub-interval is also counted. That is, the included angle in the terrain lattice belongs to L AiBj And this represents the number of topographic points in the near-fault region, calculated using the formula... Calculate the sliding proportion of each combined angle sub-interval. This is used to reflect the probability density of landslide directions occurring within the combined angled sub-intervals; simultaneously, it statistically analyzes the set of landslide directions S within the third angled sub-interval. landslide The number of landslides included That is, belonging to L Ci S i The number of points is counted, and the number of angles of topographic points belonging to the distant fault region in the topographic lattice within the third angle sub-interval is also counted. That is, the included angle in the terrain lattice belongs to L Ci And this represents the number of topographic points in the far-fault region, calculated using the formula... Calculate the sliding proportion of each third included angle sub-interval The probability density of landslide direction used to reflect the occurrence of landslides within the third included sub-interval;

[0053] Using the min-max normalization method, the sliding proportion is determined. and sliding proportion Mapping to the (0,1) interval yields the normalized probability distribution of landslide direction, and based on this, a tectonic stress direction hazard factor f is constructed. θ The structural stress direction hazard factor is expressed in piecewise function form;

[0054] The structural stress direction hazard factor is:

[0055]

[0056] in,

[0057]

[0058] Furthermore, sal1, sal2, Saa3, Saa1, Saa2, Saa3, Saa4, Saa5, Saa6, Saa7, Saa8, and Saa9 represent the normalized probability distributions of the landslide direction.

[0059] Preferably, the landslide risk of each topographic point is calculated by substituting the slope, slope length, tectonic stress, topographic stress, and included angle of the topographic point into the landslide risk formula constructed from the slope length risk factor, slope risk factor, tectonic stress direction risk factor, tectonic stress, and topographic stress. The tectonic stress of the topographic point is calculated with reference to a method for estimating the three-dimensional stress state of large-scale slopes in complex terrain areas proposed in patent CN118568965B.

[0060] The landslide risk formula is as follows:

[0061]

[0062] Among them, Risk Ai It is a landslide risk, S topo-Ai It is topographic stress, H Ai It is the structural stress, ω1 and ω2 are weighting coefficients, and ω1+ω2=1;

[0063] Using ArcGIS Pro software, the landslide risk of each topographic point is assigned to the corresponding global coordinates to generate a landslide risk raster map. Based on the classification criteria for dangerous areas, the calculation area is divided into extremely high, high, medium, low and extremely low dangerous areas for color rendering. Extremely high dangerous areas are red, high dangerous areas are orange, medium dangerous areas are yellow, low dangerous areas are blue and extremely low dangerous areas are green, resulting in a dangerous map of the calculation area.

[0064] The criteria for classifying dangerous zones are as follows: landslide risk ≥ 0.8 is an extremely high-risk zone, 0.6 ≤ landslide risk < 0.8 is a high-risk zone, 0.4 ≤ landslide risk < 0.6 is a medium-risk zone, 0.2 ≤ landslide risk < 0.4 is a low-risk zone, and landslide risk < 0.2 is an extremely low-risk zone.

[0065] Compared with the prior art, the present invention has the following beneficial effects:

[0066] This invention considers the relationship between topographic stress and tectonic stress related to slope length under the strong tectonic stress field unique to the Sichuan-Tibet region, and its impact on landslide distribution. It provides a basis for better identifying landslide hazard areas and more accurately predicting the probability of landslide disaster risk assessment. Thus, it provides a strong and quantitative scientific basis for the risk assessment, early warning and scientific prevention and control measures for landslide disasters in this region, effectively serving the actual needs of disaster prevention and mitigation in the Sichuan-Tibet region. Attached Figure Description

[0067] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0068] Figure 1 This is a flowchart of the method steps of the present invention. Detailed Implementation

[0069] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be described in detail below. Obviously, the described embodiments are merely some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0070] Examples, such as Figure 1 As shown, a landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region includes the following steps:

[0071] S1. Collect DEM topographic data, regional tectonic stress field data, GPS motion data and fault distribution data within the calculation area, and process them into DEM raster data using geographic information system tools;

[0072] S2. Using the raster-to-point tool in ArcGIS software, the DEM raster data is converted into point-like Shapefile data. The plane coordinates and elevation values ​​of each terrain point are extracted to form a terrain point array covering the calculation area. The terrain point array consists of several terrain points, and each terrain point is associated with its corresponding global coordinates.

[0073] S3. Calculate the slope length, slope, tectonic stress, and topographic stress of each topographic point in the topographic point matrix, and at the same time calculate the slope length and slope of each landslide center point in the historical landslide database;

[0074] S4. Statistically analyze the distribution patterns of slope length and slope gradient of all landslide center points in the topographic lattice, and construct a structural stress magnitude hazard factor, which includes a slope length hazard factor and a slope gradient hazard factor.

[0075] S5. Divide the calculation area into several sub-regions and determine the regional stress field direction of each sub-region. At the same time, divide each sub-region into near-fault region and far-fault region, and obtain the included angle of each topographic point accordingly.

[0076] S6. Statistically analyze the distribution pattern of the included angles of all topographic points in the topographic point matrix and the landslide direction recorded in the historical landslide database, and construct the tectonic stress direction hazard factor.

[0077] S7. Based on the hazard factors of tectonic stress magnitude and tectonic stress direction, calculate the landslide risk of each topographic point, and divide the calculation area into extremely high, high, medium, low and extremely low hazard zones according to the landslide risk value, and output the hazard map of the calculation area.

[0078] Furthermore, the working principle of the present invention will be illustrated below through embodiments:

[0079] The calculation area is the central segment of the Longmenshan fault zone in the Sichuan-Tibet region. This area is a typical zone of strong tectonic stress in the Sichuan-Tibet region, with a tectonic stress value of 0.8 × 10⁻⁶. 7 -1.2×10 7 Pa, with its complex terrain including high mountains, canyons, ravines, and flat sections of wide valleys, revealed 217 recorded landslides, which fully verifies the effectiveness of the method of this invention.

[0080] DEM topographic data at 30m resolution, regional tectonic stress field data from 2010 to 2023, GPS motion data from 2018 to 2023, and 1:50,000 fault distribution data (including the strike attributes of the Longmenshan main central fault and main boundary fault) were collected within the computational area. ArcGIS Pro 3.0 software was used to convert the DEM topographic data, regional tectonic stress field data, GPS motion data, and fault distribution data into shapefile point feature format. The projection conversion tool was used to unify them to the CGCS2000 coordinate system, which is consistent with the OXYZ global coordinate system. ArcGIS Pro's inverse distance weighting tool was used to interpolate the DEM point feature data, setting the output cell size to 30m to match the original DEM resolution, generating DEM raster data. Then, the raster-to-point tool was used to convert the DEM raster into point shapefile data, extracting the global coordinates (x, y, z) of each topographic point to form a topographic point matrix covering the area. The topographic point matrix contains 12,860 topographic points with a spacing of 30m.

[0081] A typical terrain point Ai was selected from the terrain point matrix. The global coordinates of terrain point Ai are (3320500m, -790200m, 3520m), located 5km west of the main boundary fault of Longmenshan. The terrain is a canyon slope. Starting from Ai, terrain points within a 10km range were searched in four directions (east, south, west, and north) with a step size of 100m. The elevation point sequence in each direction was recorded. For example, the elevation point sequence in the east direction is Ai→A1(3320600m, -790200m, 3535m)→A2(3320700m, -790200m, 3550m)…). The preset compression threshold S is 8° and the preset gully… The ravine threshold X is 50m. Dynamic termination is determined for the elevation point sequence in each direction. For example, in the eastward elevation point sequence, the global coordinates of A8 (3321300m, -790200m, 3640m) and the slope of the line connecting A8 and Ai is 7.5°, less than S, and the elevation difference between A8 and A7 is 45m < X. Therefore, it enters the flatness determination stage. The search continues until A8, where the cumulative horizontal distance is 1200m, greater than 700m, thus determining A8 as the eastward endpoint. Similarly, the endpoints for the south, west, and north directions are determined. Based on this, the slope vectors for each direction are calculated, and slope vectors with positive elevation differences in the four directions are selected. The slope length vector V of the terrain point Ai is synthesized according to the slope vector synthesis rules. slope Given (800m, 600m, 120m), calculate the modulus to obtain the slope length L of terrain point Ai. Ai The slope is 986m; ArcGIS Pro's slope analysis tool is used to process the DEM raster to generate a slope raster map, which is then converted into vector point data and matched with the terrain point matrix. The slope of the resulting terrain point Ai is 31.5°.

[0082] Referring to the method of patent CN118568965B, the horizontal tectonic stress tensor H0 calculated using regional tectonic stress field data is (1.0 × 10⁻⁶). 7 Pa, 0.6 × 10 7 Given that Pa, 0), and taking the highest point within a 10km range as 5000m, the tectonic stress H of Ai after superposition is... Ai Approximately 1.2 × 10 7 Pa; taking the slope length vector elevation difference as 120m, and the average rock density as 2600kg / m³. 3 The acceleration due to gravity is 9.8 m / s². 2 Converting the slope of 31.5° to 0.55 radians and substituting it into the topographic stress calculation formula, the topographic stress S at topographic point Ai is calculated. topo-Ai -1.3×10 6 Pa, the negative sign indicates that the stress direction is the same as the sliding direction (opposite to the uphill direction), and the actual value should focus on the absolute value.

[0083] 217 landslide center points were extracted from the historical landslide database. The slope length of each center point was calculated using the aforementioned slope length calculation method, forming a slope length point set. Statistics showed that 72% of the slope length points in the 500-1500m range were within this range, and 35% of the topographic points in the same range. The normalized slope length hazard factor f for this 500-1500m range was then calculated. L The value is 0.85. The slope length of the topographic point Ai, 986m, falls within this interval, therefore f L The value is 0.85; statistically, the slope at the center point of the landslide accounts for 68% of the range of 15-45°, and this range accounts for 40% of the topographic data points. The normalized slope hazard factor f for this 15-45° range is... a The value is 0.82, and the slope of topographic point Ai is 31.5°, which falls within this range. Therefore, f a The value is 0.82; the computational domain is divided into 100km sections. 2 In the sub-region, the focal mechanism solution direction of the sub-region where the topographic point Ai is located is approximately 102-130° (0 for north, positive for clockwise), the measured stress direction is approximately 129°, and the GPS displacement direction is approximately 113°. The deviation of these three values ​​is less than 30°. Taking the arithmetic mean of these three values, approximately 119°, as the stress field direction of this sub-region, the topographic point Ai is located within 5km to 10km of the fault, which is considered a near-fault region. The calculated angle θ1 between the slope aspect of topographic point Ai (105°) and the perpendicular direction of the fault strike (125°) is 20°, and the angle θ2 between the slope aspect of topographic point Ai and the regional stress field direction (83°) is 22°. Statistical analysis of historical near-fault landslide data shows that the highest proportion of landslides occurred in the interval where θ1 < 30° and θ2 < 30°, corresponding to the tectonic stress direction hazard factor f. θ It is 0.88.

[0084] Using the landslide risk formula with ω1 = 0.7 and ω2 = 0.3, and substituting the slope, slope length, tectonic stress, topographic stress, and included angle of topographic point Ai, the landslide risk of topographic point Ai was finally obtained as 0.688. According to the criteria for classifying dangerous areas, 0.6 ≤ landslide risk < 0.8 is considered a high-risk area, therefore topographic point Ai is determined to belong to the high-risk area.

[0085] The landslide risk of 12,860 topographic points in the topographic point matrix was calculated using the same steps as those used to calculate the landslide risk of topographic point Ai. The landslide risk was then assigned to the corresponding global coordinates using ArcGIS Pro to generate a landslide risk raster map. Color rendering was performed based on the classification standards to generate a hazard map of the middle section of the Longmenshan fault zone in the Sichuan-Tibet region.

[0086] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them; modifications to the technical solutions described in the foregoing embodiments, or equivalent substitutions of some of the technical features, do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region, characterized in that, Includes the following steps: S1. Collect DEM topographic data, regional tectonic stress field data, GPS motion data and fault distribution data within the calculation area, and process them into DEM raster data using geographic information system tools; S2. Using the raster-to-point tool in ArcGIS software, the DEM raster data is converted into point-like Shapefile data. The plane coordinates and elevation values ​​of each terrain point are extracted to form a terrain point array covering the calculation area. The terrain point array consists of several terrain points, and each terrain point is associated with its corresponding global coordinates. S3. Calculate the slope length, slope, tectonic stress, and topographic stress of each topographic point in the topographic point matrix, and at the same time calculate the slope length and slope of each landslide center point in the historical landslide database; S4. Statistically analyze the distribution patterns of slope length and slope gradient of all landslide center points in the topographic lattice, and construct a structural stress magnitude hazard factor, which includes a slope length hazard factor and a slope gradient hazard factor. S5. Divide the calculation area into several sub-regions and determine the regional stress field direction of each sub-region. At the same time, divide each sub-region into near-fault region and far-fault region, and obtain the included angle of each topographic point accordingly. S6. Statistically analyze the distribution pattern of the included angles of all topographic points in the topographic point matrix and the landslide direction recorded in the historical landslide database, and construct the tectonic stress direction hazard factor. S7. Based on the hazard factors of tectonic stress magnitude and tectonic stress direction, calculate the landslide risk of each topographic point, and divide the calculation area into extremely high, high, medium, low and extremely low hazard zones according to the landslide risk value, and output the hazard map of the calculation area. The landslide risk of each topographic point is calculated by substituting the slope, slope length, tectonic stress, topographic stress and included angle of the topographic point into the landslide risk formula constructed by the hazard factors of tectonic stress magnitude and tectonic stress direction.

2. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 1, characterized in that, The global coordinates are defined based on the OXYZ right-handed coordinate system, where the X-axis points east, the Y-axis points north, and the Z-axis is vertically upward with the average level of the calculation area as the zero point of the Z-axis; the global coordinates are represented as (x, y, z), where x and y are the planar coordinates of the terrain point, and z is the elevation value of the terrain point.

3. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 1, characterized in that, The slope calculation process is as follows: the slope analysis tool in ArcGIS Pro is used to process the DEM raster data to generate a slope raster map. After converting the slope raster map into vector point data, the slope of each terrain point is obtained by matching it with the terrain point matrix using the spatial connection tool.

4. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 3, characterized in that, The steps for calculating the slope length are as follows: S31. Starting from a topographic point, search for topographic points in the topographic point matrix within a 10 km radius around the starting point in the east, south, west and north directions with a step size of 100 m. Record the topographic points searched in each direction to form a sequence of elevation points in each direction. S32. The endpoint of each direction is determined by dynamic termination determination for the sequence of elevation points in each direction. The dynamic termination determination includes compression determination, gully isolation determination, flat terrain determination, and terrain trend determination. S33. Calculate the difference between the elevation value of the endpoint and the elevation value of the starting point in each direction to obtain the elevation difference in each direction. At the same time, based on the global coordinates of the endpoint and the starting point in each direction, calculate the slope vector corresponding to the elevation difference in each direction. S34. Combine the slope vectors with positive elevation differences in each direction into a slope length vector according to the slope vector synthesis rules. The magnitude of the slope length vector is the slope length.

5. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 4, characterized in that, The calculation process for the terrain stress is as follows: Obtain the elevation difference term from the slope length vector of the terrain point, and substitute the elevation difference term and the slope of the terrain point into the terrain stress calculation formula to calculate the terrain stress.

6. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 5, characterized in that, The process of constructing the structural stress magnitude risk factor is as follows: Extract the slope length and slope of all landslide center points to form a slope length point set and a slope slope point set. Calculate the proportion of the slope length point set in the slope length of the terrain point array and the proportion of the slope slope point set in the slope of the terrain point array, respectively, to obtain the slope length probability distribution and the slope probability distribution, and normalize them. Based on the normalized slope length probability distribution and slope probability distribution, construct slope length risk factors and slope risk factors with a value range of (0,1), respectively.

7. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 1, characterized in that, The process of determining the regional stress field direction of each sub-region is as follows: Calculate the average value of the focal mechanism solution direction, the average value of the measured stress direction, and the average value of the GPS displacement direction for each sub-region. If the directional deviation between the average value of the focal mechanism solution direction, the average value of the measured stress direction, and the average value of the GPS displacement direction is less than 30°, then the arithmetic mean of the average value of the focal mechanism solution direction, the average value of the measured stress direction, and the average value of the GPS displacement direction is taken as the regional stress field direction of that sub-region. Otherwise, the average value of the GPS displacement direction is taken as the regional stress field direction of that sub-region.

8. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 7, characterized in that, The near-fault region is a sub-region with a horizontal distance of less than or equal to 10 km from the fault, and the far-fault region is a sub-region with a horizontal distance of more than 10 km from the fault. The included angle of each topographic point in the near-fault region includes the angle between the slope direction of the topographic point and the direction perpendicular to the fault strike, and the angle between the slope direction of the topographic point and the direction of the regional stress field. The included angle of each topographic point in the far-fault region is the angle between the slope direction of the topographic point and the direction of the regional stress field.

9. The landslide hazard prediction method based on tectonic stress factors applicable to the Sichuan-Tibet region according to claim 8, characterized in that, The criteria for classifying dangerous zones are as follows: landslide risk ≥ 0.8 is an extremely high-risk zone, 0.6 ≤ landslide risk < 0.8 is a high-risk zone, 0.4 ≤ landslide risk < 0.6 is a medium-risk zone, 0.2 ≤ landslide risk < 0.4 is a low-risk zone, and landslide risk < 0.2 is an extremely low-risk zone.

Images