Geological disaster risk identification system and method based on topographic feature extraction

Abstract

The invention relates to the technical field of geological disaster monitoring and risk assessment, in particular to a geological disaster risk identification system and method based on topographic feature extraction, and the method comprises the steps: decomposing digital elevation model data into a reference topographic curved surface and a residual topographic anomaly field; respectively extracting macroscopic terrain factors and microscopic symptom features based on the reference terrain curved surface and the residual terrain abnormal field; and according to a preset geological disaster mode, dynamically configuring a weight, carrying out weighted fusion on the macroscopic and microscopic features, calculating a geological disaster risk index, and realizing risk grade division and visualization. According to the invention, the comprehensiveness, accuracy and reliability of geological disaster risk identification can be improved.

Description

Technical Field

[0001] This invention relates to the field of geological disaster monitoring and risk assessment technology, and in particular to a geological disaster risk identification system based on terrain feature extraction, and a geological disaster risk identification method based on terrain feature extraction. Background Technology

[0002] Geological disasters, as a common type of natural disaster, pose a serious threat to public safety and sustainable economic development due to their suddenness and destructiveness. Against the backdrop of global climate change and increasingly frequent human engineering activities, the types, scale, and frequency of geological disasters are showing a trend of increasing complexity and aggravation. Therefore, establishing an efficient and accurate geological disaster risk identification and early warning system has become a key issue in the field of disaster prevention and mitigation. In this technological context, the technical approach of analyzing and judging based on topographic features has gradually developed into one of the core research directions in this field because it can directly reflect the physical basis of geological body stability and can achieve large-scale, periodic data coverage with the help of modern surveying and mapping technology.

[0003] In the early stages of technological development, related research and applications primarily relied on the deep integration of remote sensing technology and Geographic Information Systems (GIS). Specifically, this approach typically utilizes optical satellite imagery, radar imagery, or aerial photogrammetry to acquire digital elevation model (DEM) data for large areas. Based on this, spatial analysis platforms are used to mathematically interpret and perform derivative calculations on the DEM data, thereby extracting a series of macroscopic topographic factors such as slope, aspect, surface curvature, coefficient of variation of elevation, and topographic relief in batches. By constructing evaluation models and comprehensively analyzing the spatial distribution patterns and combination characteristics of these factors, preliminary zoning and grading assessments of regional geological hazard susceptibility can be conducted.

[0004] However, with the continuous development of related technologies and the increasingly stringent requirements for risk identification accuracy and early warning timeliness in application scenarios, some inherent characteristics of the above-mentioned technological paradigms at the principle level have gradually revealed their inherent limitations in dealing with new challenges. Summary of the Invention

[0005] The purpose of this invention is to provide a geological disaster risk identification system and method based on terrain feature extraction, in order to solve the technical problem that existing technologies rely on macro-topographic factor analysis, which leads to insensitivity to key micro-level disaster precursor features, resulting in missed or incorrect risk assessments.

[0006] According to one aspect of the present invention, a geological hazard risk identification system based on terrain feature extraction is provided for assessing the geological hazard risk of a target area, comprising:

[0007] The data preprocessing unit is configured to receive and process the raw digital elevation model data covering the target area to generate a clean digital elevation model data.

[0008] The terrain decomposition unit, whose input is connected to the output of the data preprocessing unit, is configured to receive the clean digital elevation model data and decompose it into a reference terrain surface data representing the regional geomorphic framework and a residual terrain anomaly field data that reveals local micro-topographic changes.

[0009] The macro-topographic factor calculation unit has its input end connected to the output end of the topographic decomposition unit, and is configured to receive the reference topographic surface data and calculate a set of macro-topographic factors to characterize the macro-topographic stability of the target area based on the data.

[0010] A residual topographic anomaly field precursor feature extraction unit, whose input is connected to the output of the topographic decomposition unit, is configured to receive the residual topographic anomaly field data and extract a set of microscopic precursor geomorphic features related to geological body instability based on the data; and

[0011] The multi-scale feature fusion and risk identification unit has its input end connected to the output end of the macro-topography factor calculation unit and the residual topographic anomaly field sign feature extraction unit, respectively. It is configured to receive the macro-topography factor and the micro-precursor geomorphic features, and perform weighted fusion of the two according to a built-in decision model to calculate and output a risk index that characterizes the degree of geological disaster risk in the target area.

[0012] In some embodiments, the terrain decomposition unit includes:

[0013] A two-dimensional Gaussian low-pass filter is configured to perform spatial domain convolution operations on the input clean digital elevation model data using a preset convolution kernel and standard deviation parameters to filter out high-frequency components in the elevation data, thereby generating the reference topographic surface data. The size of the convolution kernel is set to be larger than an empirical threshold greater than the key scale of the tectonic structures responsible for major geological hazards in the target area.

[0014] A matrix difference operator is configured to receive the clean digital elevation model data and the reference terrain surface data generated by the two-dimensional Gaussian low-pass filter, and perform pixel-by-pixel matrix subtraction to generate the residual terrain anomaly field data, wherein the value of each pixel in the residual terrain anomaly field data represents the degree of deviation of the true elevation of that point relative to its local macroscopic geomorphic background.

[0015] In some embodiments, the macro-topographic factor calculation unit integrates multiple topographic parameter calculation modules, and the macro-topographic factors include at least: slope and aspect, plane curvature and profile curvature, topographic relief and surface roughness.

[0016] In some embodiments, the residual terrain anomaly field feature extraction unit includes three parallel feature extraction subunits: a linear structure feature extraction subunit, a morphological feature extraction subunit, and a texture complexity feature extraction subunit.

[0017] In some embodiments, the linear structural symptom extraction subunit is internally configured with a Hessian matrix-based structural enhancement filter.

[0018] In some embodiments, the morphological sign extraction subunit is internally configured with a multi-scale top-hat and bottom-hat transformation operator based on mathematical morphology theory.

[0019] In some embodiments, the texture complexity symptom extraction subunit is internally configured with a texture parameter calculator based on the gray-level co-occurrence matrix.

[0020] In some embodiments, the multi-scale feature fusion and risk identification unit is internally configured with a geological disaster risk index calculation engine and a geological disaster model and disaster-causing mechanism knowledge base.

[0021] According to another aspect of the present invention, a method for identifying geological hazard risks based on terrain feature extraction is provided, comprising:

[0022] Acquire and preprocess digital elevation model data covering the target area;

[0023] Using two-dimensional Gaussian low-pass filtering and matrix difference operations, the preprocessed digital elevation model data is decomposed into benchmark topographic surface data and residual topographic anomaly field data.

[0024] Based on the aforementioned benchmark terrain surface data, a set of macroscopic terrain factors, including slope, aspect, and curvature, is calculated.

[0025] Based on the residual terrain anomaly field data, linear structure features were extracted using a Hessian matrix analysis method, morphological features were extracted using a mathematical morphology multi-scale top and bottom cap transformation method, and texture complexity features were extracted using a gray-level co-occurrence matrix method.

[0026] Based on a pre-defined geological hazard model that matches the regional geological background, weights are dynamically configured to weight and fuse the extracted macroscopic topographic factors and microscopic symptom features to calculate the geological hazard risk index. Based on the calculated risk index, the target area is classified into risk levels, and a visual risk distribution map is generated.

[0027] Compared with existing technologies, this invention has the following advantages: By constructing an innovative technical framework of decomposition-parallel analysis-fusion, it successfully resolves the inherent contradiction between macroscopic analysis and microscopic identification in traditional methods; this invention can deeply mine key disaster precursor information that was previously ignored from conventional high-precision DEMs without relying on ultra-high-cost data sources, and places the macroscopic disaster-inducing environment and microscopic instability signs in a unified analytical framework based on geomechanical mechanisms for collaborative evaluation, which greatly improves the comprehensiveness, accuracy and reliability of geological disaster risk identification, and has significant technological progress and important application value. Attached Figure Description

[0028] To more clearly illustrate the technical solutions 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.

[0029] Figure 1 This is a structural block diagram of the geological disaster risk identification system provided in an embodiment of the present invention;

[0030] Figure 2 This is a flowchart illustrating the geological disaster risk identification method provided in this embodiment of the invention.

[0031] In the figure: 10 - Data preprocessing unit; 20 - Terrain decomposition unit; 30 - Macro-terrain factor calculation unit; 40 - Residual terrain anomaly field feature extraction unit; 50 - Multi-scale feature fusion and risk identification unit; 60 - Risk information visualization and output unit. Detailed Implementation

[0032] The following will refer to the appendices in the embodiments of the present invention. Figure 1 -Appendix Figure 2 The technical solutions in the embodiments of the present invention will be clearly and completely described together. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0033] Example 1

[0034] The geological hazard risk identification system of the present invention can be physically carried by a server or workstation equipped with a high-performance graphics processing unit (GPU), with Linux as the operating system, and with Geographic Information System (GIS) software (such as GDAL, SAGA GIS library) and scientific computing and image processing libraries (such as NumPy, SciPy, OpenCV) installed.

[0035] The logical structure of the system is shown in the appendix. Figure 1 As shown, its overall architecture consists of a series of functionally interconnected units, including a data preprocessing unit 10, a terrain decomposition unit 20, a macro-terrain factor calculation unit 30, a residual terrain anomaly field feature extraction unit 40, a multi-scale feature fusion and risk identification unit 50, and a risk information visualization and output unit 60. These units work together through internally defined data interfaces and processing flows to complete the entire transformation process from raw terrain data to the final risk distribution map.

[0036] Specifically, the data preprocessing unit 10, acting as the system's input portal, is responsible for receiving and standardizing the raw digital elevation model (DEM) data covering the target area. DEM data is the foundation of this system's analysis, and its accuracy directly affects the system's ability to resolve micro-topographic features. Therefore, this embodiment limits the use of data with a spatial resolution of at least 5 meters, preferably using 1-meter resolution DEM data acquired through LiDAR or UAV oblique photogrammetry. The internal processing flow of the data preprocessing unit 10 begins with the data format conversion module. This module has built-in parsers for various mainstream raster data formats, capable of uniformly converting input DEM files in formats such as GeoTIFF, Erdas Imagine (.img), and ASCII Grid into a binary floating-point raster format for efficient storage and computation within the system. Next, the coordinate system module performs georegistration checks on the transformed data. Equipped with the PROJ coordinate transformation library, it can accurately reproject data from different geodetic datums (such as WGS-84, Beijing 54, Xi'an 80) and projection methods (such as UTM, Gauss-Kruger) to the user-specified target coordinate system, such as the corresponding projection zone under the CGCS2000 National Geodetic Coordinate System, thus ensuring the consistency of the datum for all subsequent spatial analyses. The final step in the processing flow is completed by the outlier removal module. To prevent isolated elevation peaks or depressions introduced by sensor noise or data stitching errors from interfering with the analysis, this module employs a filter based on local statistics. Specifically, it traverses the entire DEM with a 5x5 pixel window, calculating the mean μ and standard deviation σ of the elevation values ​​of all pixels within the window except the center pixel. If the elevation value of the center pixel deviates from the mean μ by more than 3σ, it is identified as an outlier, and its elevation value is replaced with the mean μ of the window, thereby generating a clean DEM data with continuous elevation and logical consistency. The clean DEM data is passed to the topographic decomposition unit 20 as standardized input.

[0037] The terrain decomposition unit 20 is the core embodiment of the technical concept of this invention. Its function is to decompose the raw DEM data containing complex information into two components with clear physical meaning. This unit receives the clean DEM output by the data preprocessing unit 10, and its core is a two-dimensional Gaussian low-pass filter. The function of this filter is to smooth the terrain, simulate the large-scale geomorphic evolution process in nature, and thus extract the reference terrain surface representing the geomorphic skeleton of the region. The key parameters of the filter, namely the size of the convolution kernel and the standard deviation (σ) of the Gaussian function, need to be set according to the type of geological hazard in the area to be evaluated. For example, for large bedrock landslides, the scale of the hazard-generating structure is usually on the order of hundreds of meters, and the size of the convolution kernel can be set to 51x51 pixels (corresponding to a range of 51 meters by 51 meters); for smaller loess landslides or collapses, the scale of the controlling terrain is smaller, and the size of the convolution kernel can be adjusted accordingly to 31x31 pixels. In this embodiment, the standard deviation σ of the Gaussian function is strictly set to one-sixth of the kernel size N, i.e., σ = N / 6. This relationship ensures that the Gaussian function decays to near zero at the kernel boundary, avoiding truncation effects. The specific process of the Gaussian low-pass filtering operation is to perform a two-dimensional convolution of this Gaussian convolution kernel with the DEM data, generating a smooth reference topographic surface raster data that retains only the macroscopic topographic undulation trend. Subsequently, the matrix difference operator inside the topographic decomposition unit 20 starts working. This operator simultaneously reads the original clean DEM data and the newly generated reference topographic surface data, both of which have identical dimensions and registration information. It performs a subtraction operation pixel by pixel, i.e., for each pixel (i,j), it calculates its residual value Residual(i,j) = DEM_original(i,j) - Surface_base(i,j). The set of residual values ​​for all pixels constitutes a new raster image, i.e., residual topographic anomaly field data. Positive values ​​in the data field indicate a macro-topographic trend where the surface elevation is higher than the surrounding area, possibly corresponding to bulging deformation or compression ridges in the front of the landslide body; negative values ​​indicate a macro-topographic trend where the surface elevation is lower than the surrounding area, possibly corresponding to tensional fissures, differential settlement zones, or gullies at the rear edge of the landslide; and areas with values ​​close to zero indicate a smooth surface, consistent with the macro-topographic background. Through this decomposition process, the topographic decomposition unit 20 divides the original topographic information into two parts, and transmits the baseline topographic surface data to the macro-topographic factor calculation unit 30, while simultaneously transmitting the residual topographic anomaly field data to the residual topographic anomaly field feature extraction unit 40, laying the foundation for subsequent parallel analysis.

[0038] The macro-topographic factor calculation unit 30 is tasked with calculating a set of traditional topographic parameters characterizing regional geomorphic stability based on smoothed benchmark topographic surface data. Because its input data has been filtered to remove noise and disturbances from the microscopic surface, its calculation results better reflect controlling, large-scale geomorphic features, improving the stability and accuracy of these traditional factors. This unit integrates a series of mature topographic analysis algorithm modules. The slope and aspect calculation module employs the third-order inverse distance squared weighted finite difference method. Compared to the simple third-order finite difference method, this algorithm assigns different weights to the eight neighboring points around the central pixel, with closer points receiving higher weights, thus more accurately calculating the slope (maximum rate of elevation change) and aspect (azimuth angle of the direction of maximum slope) for each pixel. The curvature calculation module uses the nine-parameter quadratic polynomial fitting method proposed by Zevenbergen & Thorne in 1987. For each pixel, this module uses its own elevation and the nine elevation values ​​within its 3x3 neighborhood to fit a local quadratic surface equation z = ax². 2 +by 2 +cxy+dx+ey+f. By solving the coefficients of the equation, the section curvature at that point can be directly calculated as (-(2b) / ((d)). 2 +e 2 )+1)) and plane curvature ((2a) / ((d) 2 +e 2 The former reflects the vertical slope undulation, while the latter reflects the horizontal confluence and divergence characteristics. Furthermore, the terrain relief calculation module quantifies the local terrain cutting depth by setting a statistical window (e.g., a circular window with a radius of 100 meters) and calculating the difference between the maximum and minimum elevations of the reference terrain surface within the window. The surface roughness calculation module characterizes the ruggedness and complexity of the surface by calculating the ratio of the surface area (based on slope calculation) of each pixel unit to its horizontal projected area. All these calculated macroscopic terrain factors, including slope map, aspect map, plan curvature map, profile curvature map, terrain relief map, and surface roughness map, are organized into a multi-layered macroscopic feature dataset and then fed as a whole into the multi-scale feature fusion and risk identification unit 50.

[0039] Meanwhile, the residual topographic anomaly feature extraction unit 40 performs in-depth analysis on another key data set output by the topography decomposition unit 20: the residual topographic anomaly field. The design goal of this unit is to identify and quantify micro-topographic anomalies directly related to the early instability of geological bodies, which are often overlooked in macro-topographic analysis. To this end, this unit integrates three feature extraction sub-units targeting different types of micro-signs: a linear structure sign extraction sub-unit, a morphological sign extraction sub-unit, and a texture complexity sign extraction sub-unit.

[0040] Furthermore, the linear structure feature extraction sub-unit focuses on detecting linear or banded structures in the residual field, which geologically correspond to the initiation of surface tensional fractures, shear fracture surfaces, or erosion grooves. The core algorithm of this sub-unit is a structure enhancement filter based on Hessian matrix eigenvalue analysis.

[0041] The workflow is as follows: First, the residual topographic anomaly field is treated as a two-dimensional grayscale image, and the second-order partial derivative of each pixel is calculated to form a 2x2 Hessian matrix H = [[I_xx,I_xy],[I_yx,I_yy]], where I_xx, etc., represent the second-order derivative of the image in the corresponding direction. This is approximated by performing finite difference on the Gaussian smoothed residual field. Second, the two eigenvalues ​​λ1 and λ2 of this Hessian matrix are calculated (assuming |λ1|≥|λ2|). Based on the geometric characteristics of tubular structures (such as fissures), the ideal Hessian response is: a large second-order derivative in the direction across the structure (λ1 has a large absolute value and is negative, corresponding to a valley shape), while the second-order derivative is close to zero in the direction extending along the structure (the absolute value of λ2 is close to zero). Based on this principle, this sub-unit constructs a linear structural response function, such as the tubular enhancement function defined in the Frangi filter: V(λ1,λ2)=exp(-R_B 2 / (2*β2))*(1-exp(-S 2 / (2*c 2 ))), where R_B=|λ2 / λ1| is a blob metric, S=sqrt(λ1 2 +λ2 2 ) is a structural strength metric, and β and c are parameters controlling the sensitivity of the filter. This function provides extremely high response values ​​for pixels that satisfy |λ2|≈0 and λ1 is negative. By applying this function to the entire residual field, a linear structural saliency map can be generated, and the highlighted linear regions in the map represent potential geological fissures or early gullies identified by the system.

[0042] Furthermore, the morphological feature extraction subunit aims to identify and quantify local anomalies in the residual field that exist in clumps or patches. These anomalies typically indicate deformations with more three-dimensional volumetric effects, such as creeping bulges at the leading edge of a landslide (positive anomalous clumps) or saucer-shaped depressions at the trailing edge of a slope due to traction (negative anomalous clumps). This subunit achieves this goal using multi-scale top-hat and bottom-hat transformations from mathematical morphology. The top-hat transformation is defined as subtracting the original image from its "opening" operation (erosion followed by dilation), which effectively extracts all bright patches smaller than the structuring element in the image. The bottom-hat transformation is defined as subtracting the original image from its "closing" operation (dilation followed by erosion), which extracts all dark holes or depressions. To capture deformable bodies of different sizes, this sub-unit does not use a single-sized structuring element. Instead, it employs a sequence of circular structuring elements with radii progressively increasing from 3 pixels to 15 pixels (with a step size of 2 pixels), resulting in circular kernels with radii of 3, 5, 7, 9, 11, 13, and 15. For each structuring element, a top-hat transformation and a bottom-hat transformation are performed on the residual topographic anomaly field. This generates a series (in this example, 7) of bulge feature maps and a series (7) of subsidence feature maps at different scales. These maps collectively constitute a comprehensive characterization of local bulges and depressions.

[0043] Furthermore, the texture complexity symptom extraction sub-unit identifies subtle surface disturbances or fragmentation caused by stress adjustments or gradual destruction within the geological body by quantifying the local texture features of the residual field. These phenomena manifest as chaotic and irregular variations in grayscale values ​​in the residual field. The core of this sub-unit is a texture analyzer based on the Gray-Level Co-occurrence Matrix (GLCM).

[0044] The workflow is as follows: First, to reduce computational complexity, the floating-point residual field data is linearly quantized into 64 gray levels. Then, a sliding window size (e.g., 11x11 pixels) and a displacement vector (e.g., d = (1,1), representing the lower right diagonal direction) are set. The analyzer traverses the entire quantized residual field with this window, and within each window, it counts the gray value co-occurrence frequency of all pixel pairs that satisfy the displacement vector relationship, thereby constructing a 64x64 gray-level co-occurrence matrix. Based on this matrix, a set of predefined statistics that reflect texture characteristics are calculated. In this embodiment, these statistics are defined as four core metrics: Angular Second Moment (ASM), calculated as ΣΣp(i,j)², which measures the uniformity of the texture; Contrast (CON), calculated as ΣΣ(ij)²p(i,j), which measures the sharpness and groove depth of the texture; Correlation (COR), calculated as ΣΣ[(i-μ_i)(j-μ_j)p(i,j)] / (σ_iσ_j), which measures the directionality of the texture; and Entropy (ENT), calculated as -ΣΣp(i,j)log(p(i,j)), which measures the complexity and randomness of the texture. By performing a sliding window analysis on the entire residual field, four texture feature maps are finally generated, representing the spatial distribution of Angular Second Moment, Contrast, Correlation, and Entropy, respectively.

[0045] All the microscopic symptom feature maps extracted from the above three sub-units, including linear structure saliency maps, multi-scale bulging / settling feature map sets, and four texture feature maps, are integrated into a multi-dimensional microscopic symptom feature dataset, and transmitted together with the macroscopic feature dataset to the multi-scale feature fusion and risk identification unit 50.

[0046] The multi-scale feature fusion and risk identification unit 50 is the decision-making center of the entire system. It receives two major categories of feature datasets, macroscopic and microscopic, from upstream units. Its core task is to intelligently fuse this multi-source information based on a decision model with built-in geomechanical prior knowledge and calculate the comprehensive geological hazard risk index (GHRI) for each pixel. Its internal GHRI calculation engine first performs normalization processing, using the min-max normalization method to uniformly linearly map the numerical range of all input feature layers (whether it is slope, curvature, linear structural response value, or texture entropy) to the interval [0, 1] to eliminate dimensional differences and ensure the comparability of different features in the subsequent weighted model. After normalization, GHRI is calculated through a hierarchical weighted summation model, the general form of which is: GHRI=α*Σ(W_macro_i*F_macro_i)+β*Σ(W_micro_j*F_micro_j). Wherein, F_macro_i and F_micro_j are the normalized macroscopic and microscopic feature values, W_macro_i and W_micro_j are the weights of each specific feature, and α and β are the overall weights of the two categories of macroscopic and microscopic features, satisfying α+β=1. The key innovation of this invention lies in the fact that these weight coefficients are not static constants, but are dynamically adjusted by a built-in "Geological Hazard Model and Disaster-Causing Mechanism Knowledge Base." This knowledge base is a structured database that pre-stores typical combinations of disaster-causing factors and weight configuration schemes for different types of geological hazards (such as bedding bedrock landslides, loess collapsible landslides, high-altitude long-distance debris flows, and unstable rockfalls) under specific geological environments (such as lithology, structure, and rainfall conditions). Before performing a risk assessment, the operator needs to select the most suitable geological hazard model based on the actual conditions of the study area. For example, if the assessment area is a granite residual slope, prone to rainfall-induced shallow landslides, the knowledge base will automatically apply a weighting scheme. This scheme may assign higher weights to macroscopic features such as slope and runoff capacity (reflected by planar curvature), while emphasizing the weights of morphological features such as dish-shaped depressions (indicating potential catchment and saturated areas) and textural complexity such as entropy and contrast (indicating soil disturbance). When the assessment object becomes a hard rock slope in a tectonic fracture zone, the knowledge base switches to a different weighting scheme, significantly increasing the weight of macroscopic features such as slope aspect and structural surface dip consistency, and placing great emphasis on microscopic features such as linear structural features (indicating the exposure of continuous structural surfaces). This dynamic weighting mechanism based on prior knowledge allows the risk assessment model to closely align with the physical mechanisms of specific geological hazards, thereby achieving highly adaptive and accurate identification.

[0047] Finally, the risk information visualization and output unit 60 is responsible for organizing and presenting the final results calculated by the risk identification unit 50. This unit receives a continuous GHRI index raster map covering the entire study area. Its internal risk level classification module uses a statistical clustering algorithm called "Jenks Natural Breaks," which can identify inherent groupings in the data and automatically divide continuous GHRI index values ​​into five distinct discrete risk levels: very low risk, low risk, medium risk, high risk, and very high risk. Subsequently, the color mapping module assigns standardized color bands from green (very low risk) to yellow (medium risk) to red (very high risk) to these five levels, generating an intuitive and easy-to-read geological hazard risk distribution map. This unit also provides an interactive query tool, allowing users to click on any point on the generated risk map. The system will pop up an information box listing the geographic coordinates, total GHRI index, specific risk level, and the original values, normalized values, and weights of all macroscopic and microscopic features that constitute the index, helping users understand the specific reasons for the risk rating of that point. The final risk distribution map, risk zoning statistics table, and related raster data layers can be exported to common GIS file formats (such as GeoTIFF with layer styles, Esri Shapefile), or generated into printable high-resolution maps, providing direct, quantitative, and insightful scientific basis for land spatial planning, site selection of major projects, and deployment of disaster monitoring and early warning networks.

[0048] To verify the effectiveness and superiority of the method proposed in this invention, a section of a riverside highway corridor in a mountainous area of ​​southwestern China was selected as the study area. This region is located in a tectonic zone, with bedrock consisting of interbedded layered sandstone and mudstone. The dip and aspect of the rock strata are nearly consistent with the slope in some sections, and several bedding-parallel landslides have occurred historically. The study area covers 5 km x 8 km, and a 1-meter resolution digital elevation model (DEM) obtained from UAV aerial surveys was used as the basic data.

[0049] First, the data is processed according to the technical solution of this invention. The data preprocessing unit 10 performs format unification, coordinate registration (CGCS2000 projected coordinate system), and outlier removal on the original DEM. In the terrain decomposition unit 20, considering the landslide-prone scale of the study area, a 41x41 pixel Gaussian low-pass filter (standard deviation σ = 41 / 6 ≈ 6.83) is selected to generate a reference terrain surface, and then the residual terrain anomaly field is obtained through matrix differencing. The macroscopic terrain factor calculation unit 30 calculates six factors based on the reference terrain surface, including slope, aspect, profile curvature, and planar curvature. The residual terrain anomaly field feature extraction unit 40 operates in parallel: the linear structure feature extraction subunit uses a Frangi filter (parameters β=0.5, c=10) to extract linear negative anomaly structures; the morphological feature extraction subunit uses a sequence of circular structuring elements with radii of 3, 7, 11, and 15 pixels to perform multi-scale top-hat and bottom-hat transformations, and synthesizes the maximum values ​​of the results at each scale to obtain a comprehensive bulging and subsidence feature map; the texture complexity feature extraction subunit uses a 15x15 pixel window to calculate the angular second moment, contrast, correlation, and entropy of the GLCM.

[0050] In the multi-scale feature fusion and risk identification unit 50, based on the geological background of the study area as "bedslide along bedding planes," a corresponding weighting scheme was retrieved from the knowledge base. This scheme sets the overall weights for macroscopic and microscopic features at α = 0.4 and β = 0.6, reflecting a high degree of attention to microscopic instability precursors. Within macroscopic features, the weights are: slope (0.3), consistency between slope aspect and rock strata dip (calculated from externally input rock strata attitude and slope aspect data) (0.4), profile curvature (0.2), and others (0.1). Within microscopic features, the weights are: linear structure significance (0.5), morphological bulging features (0.2), morphological subsidence features (0.1), texture entropy (0.1), and texture contrast (0.1). Based on this dynamic weighted model, the GHRI index is calculated, and the risk information visualization and output unit 60 uses the natural breakpoint method to divide the area into five risk levels and generate maps.

[0051] The results showed that the method of this invention identified a total of 7 high-risk areas and 2 extremely high-risk areas. By comparing with historical disaster data and high-resolution remote sensing images, 8 of these 9 areas corresponded perfectly to known old landslides or potentially unstable slopes with obvious deformation signs (such as back wall tensile cracking, lateral shearing, and leading edge bulging). In particular, one area assessed as high-risk only showed a slight tonal anomaly and vegetation disturbance in the central part of the slope in the remote sensing image, which was difficult to identify using traditional methods. However, the residual field analysis of this invention clearly revealed a discontinuous and significant linear negative anomaly band (corresponding to tensile cracks) at its rear edge, and the central part of the slope had high morphological bulging and texture entropy values, thus successfully capturing the area.

[0052] For comparison, a representative traditional geological hazard risk assessment method was used to analyze the same 1-meter resolution DEM data of the same study area. This method does not perform topographic decomposition and directly calculates macroscopic topographic factors (slope, aspect, curvature, topographic relief, and elevation) on the original DEM. The risk assessment model uses expert scoring combined with the Analytic Hierarchy Process (AHP) to determine the fixed weights of each factor, set as follows: slope (0.35), elevation (0.15), aspect (0.10), curvature (0.20), and topographic relief (0.20). This method does not include the extraction and analysis of microscopic topographic anomalies. After weighted overlay of the factor layers, the natural breakpoint method was also used to divide the area into five risk levels.

[0053] The comparative results showed that a total of 12 high-risk areas and 4 extremely high-risk areas were identified. The identified high-risk areas were generally large, encompassing numerous steep but geologically stable slopes showing no signs of deformation. In comparison with known disaster sites, the method successfully identified 6 of the 9 known unstable slopes, but missed 3. The 3 missed sites were all small-scale landslides or landslides in the early stages of deformation; their macroscopic topographic features (such as slope) were not extremely prominent, but they showed significant signs of local deformation. Furthermore, the weakly deformed areas successfully identified by the method of this invention were classified as medium-risk areas in the comparative results, failing to raise any alerts.

[0054] The following table compares the embodiments of the present invention with the comparative examples:

[0055]

[0056]

[0057] The quantitative comparison of the data in the table above clearly shows that the technical solution provided by this invention has significantly improved the recognition accuracy compared with the traditional evaluation method that relies solely on macroscopic terrain factors, while greatly reducing the false negative rate and the false positive area ratio.

[0058] Example 2

[0059] Based on the same inventive concept as the geological hazard risk identification system based on terrain feature extraction in Embodiment 1 above, such as Figure 2 As shown, the present invention also provides a geological hazard risk identification method based on terrain feature extraction, including:

[0060] The process involves acquiring and preprocessing digital elevation model (DEM) data covering the target area; decomposing the preprocessed DEM data into baseline topographic surface data and residual topographic anomaly field data using two-dimensional Gaussian low-pass filtering and matrix difference operations; calculating a set of macroscopic topographic factors, including slope, aspect, and curvature, based on the baseline topographic surface data; extracting linear structure features using Hessian matrix analysis, extracting morphological features using mathematical morphology multi-scale top and bottom cap transformation, and extracting texture complexity features using gray-level co-occurrence matrix; dynamically configuring weights based on a preset geological hazard model matching the regional geological background, and weighted fusing the extracted macroscopic topographic factors and microscopic feature characteristics to calculate a geological hazard risk index; and classifying the target area into risk levels based on the calculated risk index and generating a visualized risk distribution map.

[0061] The specific example of the geological hazard risk identification system based on terrain feature extraction in the aforementioned embodiment 1 is also applicable to the geological hazard risk identification method based on terrain feature extraction in this embodiment. Through the foregoing detailed description of the geological hazard risk identification system based on terrain feature extraction, those skilled in the art can clearly understand the geological hazard risk identification method based on terrain feature extraction in this embodiment. Therefore, for the sake of brevity, it will not be described in detail here.

[0062] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered illustrative and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the scope of the invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0063] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A geological hazard risk identification system based on terrain feature extraction, used for geological hazard risk assessment of a target area, characterized in that, include: The data preprocessing unit (10) is configured to receive and process the raw digital elevation model data covering the target area to generate a clean digital elevation model data. The terrain decomposition unit (20) is connected to the output of the data preprocessing unit (10) and is configured to receive the clean digital elevation model data and decompose it into a reference terrain surface data representing the regional geomorphic skeleton and a residual terrain anomaly field data that reveals local micro-topographic changes. The macro-topographic factor calculation unit (30) is connected to the output of the topographic decomposition unit (20) and is configured to receive the reference topographic surface data and calculate a set of macro-topographic factors to characterize the macro-topographic stability of the target area based on the data. The residual topographic anomaly field sign feature extraction unit (40) has its input end connected to the output end of the topographic decomposition unit (20) and is configured to receive the residual topographic anomaly field data and extract a set of microscopic precursor geomorphic features related to the instability of the geological body based on the data. as well as The multi-scale feature fusion and risk identification unit (50) has its input end connected to the output end of the macro-topography factor calculation unit (30) and the residual topographic anomaly field symptom feature extraction unit (40), respectively. It is configured to receive the macro-topography factor and the micro-precursor geomorphic features, and perform weighted fusion of the two according to a built-in decision model to calculate and output a risk index that characterizes the degree of geological disaster risk in the target area.

2. The method according to claim 1, characterized in that, The terrain decomposition unit (20) includes: A two-dimensional Gaussian low-pass filter is configured to perform spatial domain convolution operations on the input clean digital elevation model data using a preset convolution kernel and standard deviation parameters to filter out high-frequency components in the elevation data, thereby generating the reference topographic surface data. The size of the convolution kernel is set to be larger than an empirical threshold greater than the key scale of the tectonic structures responsible for major geological hazards in the target area. A matrix difference operator is configured to receive the clean digital elevation model data and the reference terrain surface data generated by the two-dimensional Gaussian low-pass filter, and perform pixel-by-pixel matrix subtraction to generate the residual terrain anomaly field data, wherein the value of each pixel in the residual terrain anomaly field data represents the degree of deviation of the true elevation of that point relative to its local macroscopic geomorphic background.

3. The method according to claim 1, characterized in that, The macro-topographic factor calculation unit (30) integrates multiple topographic parameter calculation modules. The macro-topographic factors include at least: slope and aspect, plane curvature and profile curvature, topographic relief and surface roughness.

4. The method according to claim 1, characterized in that, The residual terrain anomaly field feature extraction unit (40) includes three parallel feature extraction subunits: linear structure feature extraction subunit, morphological feature extraction subunit, and texture complexity feature extraction subunit.

5. The method according to claim 4, characterized in that, The linear structure symptom extraction subunit is internally configured with a structure enhancement filter based on the Hessian matrix.

6. The method according to claim 4, characterized in that, The morphological sign extraction subunit is internally configured with a multi-scale top-hat and bottom-hat transformation operator based on mathematical morphology theory.

7. The method according to claim 4, characterized in that, The texture complexity symptom extraction subunit is internally configured with a texture parameter calculator based on the gray-level co-occurrence matrix.

8. The method according to claim 1, characterized in that, The multi-scale feature fusion and risk identification unit (50) is internally configured with a geological disaster risk index calculation engine and a geological disaster model and disaster-causing mechanism knowledge base.

9. A method for identifying geological hazard risks based on terrain feature extraction, characterized in that, include: Acquire and preprocess digital elevation model data covering the target area; Using two-dimensional Gaussian low-pass filtering and matrix difference operations, the preprocessed digital elevation model data is decomposed into benchmark topographic surface data and residual topographic anomaly field data. Based on the aforementioned benchmark terrain surface data, a set of macroscopic terrain factors, including slope, aspect, and curvature, is calculated. Based on the residual terrain anomaly field data, linear structure features were extracted using a Hessian matrix analysis method, morphological features were extracted using a mathematical morphology multi-scale top and bottom cap transformation method, and texture complexity features were extracted using a gray-level co-occurrence matrix method. Based on a pre-defined geological hazard model that matches the regional geological background, weights are dynamically configured, and the extracted macroscopic topographic factors and microscopic symptom features are weighted and fused to calculate the geological hazard risk index. Based on the calculated risk index, the target area is classified into risk levels, and a visual risk distribution map is generated.

Images