Geological disaster detection method and system based on image processing

Abstract

The present application relates to the technical field of disaster detection, in particular to a geological disaster detection method and system based on image processing, comprising the following steps: collecting multi-temporal topographic images and dividing them into blocks, constructing a main boundary angle difference value array, calculating the included angle change to screen image blocks and combining fault regions, extracting crack contours to calculate the opening angle and construct the offset sequence, generating an enhanced feature map based on gradient projection texture field interpolation, constructing a residual data field to divide grid points, and determining abnormal regions. In the present application, the main boundary angle difference value array of the blocks is constructed and affine aggregation is performed in combination with the consistency of the edge gray gradient, so that the precise reconstruction and continuity recovery of the cross-scale nonlinear fault structure in the topographic image are realized, the gradient derivative quantity is generated by using the crack endpoint opening angle offset sequence and the texture vector field is expanded, which can sensitively capture the weak morphological changes of the local topography in the time evolution process, and solve the problem that the small disaster precursor features are difficult to identify in a complex background.

Description

Technical Field

[0001] This invention relates to the field of disaster detection technology, and in particular to a geological disaster detection method and system based on image processing. Background Technology

[0002] The field of geological hazard detection technology involves monitoring, identifying, and assessing areas in the geological environment that may trigger disasters such as landslides, collapses, and debris flows. Core aspects include the collection and analysis of information on surface deformation, geomorphological features, changes in soil structure, and the propagation of rock fissures, as well as the determination of disaster risk areas and the generation of early warning information. This technological field relies on remote sensing monitoring, 3D modeling, image processing, and geological modeling, combining multi-source image information such as satellite imagery, UAV aerial images, and ground cameras. By comparing and analyzing abnormal surface changes, it assists in disaster prediction. In this field, image processing technology has become an important tool for identifying early signs of disasters. By extracting and analyzing image texture, edges, shape, color, and other feature parameters, a preliminary assessment of dangerous areas can be achieved.

[0003] Traditional geological hazard detection methods rely on image processing to identify geological hazards. The main technical challenge is extracting and analyzing topographic change features from images to identify and monitor signs of hazards such as landslides and collapses. Traditional methods typically employ edge detection operators based on grayscale changes to extract areas of mountain contour variation, utilize morphological processing to identify crack orientation, analyze the consistency of rock mass distribution in the image using grayscale co-occurrence matrices, and combine threshold segmentation to determine potential hazard areas. Furthermore, after image registration, traditional methods also use template matching to compare multi-temporal images to identify areas of significant surface change, which are then used for geological hazard risk assessment.

[0004] However, existing geological hazard detection technologies still have significant limitations in practical applications. First, traditional image processing methods based on grayscale thresholding or simple edge detection often struggle to effectively distinguish minute cracks (typically on the millimeter or centimeter scale) in the early stages of a disaster when faced with complex backgrounds (such as vegetation cover or gravel accumulation). These minute cracks have extremely weak grayscale features and are easily masked by image noise, resulting in a high rate of missed detections.

[0005] Secondly, geological fault structures often exhibit nonlinear broken lines or curves across scales, rather than ideal straight lines. Existing block-based processing methods are usually based on the assumption of continuity within the blocks, ignoring the abrupt angular changes at the block boundaries. This leads to "broken chains" or spatial misalignments when reconstructing faults, making it impossible to accurately restore the true fault zone orientation.

[0006] Furthermore, existing monitoring technologies primarily focus on geomorphological features at a single moment, lacking dynamic analysis of the evolution of texture fields over time. For example, the opening of cracks is often accompanied by a slight deflection of the texture direction of the surrounding rock mass; this change in the "field" often occurs earlier than a simple change in the "line." However, existing technologies rarely utilize this gradient-derived quantity to drive texture enhancement, making it difficult to capture early precursory information of disasters. Therefore, there is an urgent need for a geological disaster detection method that can sensitively capture subtle local geomorphological anomalies and accurately reconstruct cross-scale nonlinear fault structures. Summary of the Invention

[0007] To address the technical problems existing in the prior art, embodiments of the present invention provide a geological disaster detection method based on image processing, comprising the following steps: S1: Collect multi-temporal geomorphic images of the target monitoring area, divide the multi-temporal geomorphic images into image blocks according to the spatial scale sequence, and extract the structural fold vectors at the boundaries of the image blocks to construct the block main boundary angle difference array; S2: Based on the block main boundary angle difference array, calculate the angle change value between adjacent image blocks, filter image blocks that meet the angle and consistency thresholds and combine them affinely to generate tomographic aggregation region data; S3: Extract the fracture edge contour line based on the fault aggregation area data, calculate the dual-temporal endpoint subtraction angle of the fracture edge contour line endpoints and map it to the direction field to construct the fracture endpoint subtraction angle offset sequence. S4: For the crack endpoint angular offset sequence, calculate the gradient derivation amount, project the gradient derivation amount to the texture field of the multi-temporal geomorphic image, interpolate the texture vectors that meet the projection screening threshold, and generate enhanced texture feature map data. S5: Based on the enhanced texture feature map data, construct a registration residual data field and divide it into grid windows. Calculate the number of morphological contours, fluctuation frequency, and span ratio within the grid windows. Determine the evolutionary anomalous region based on the difference response classification and superposition results of the number of morphological contours, fluctuation frequency, and span ratio.

[0008] As a further aspect of the present invention, the block main boundary angle difference array includes a block index identifier, boundary normal vector angle difference elements, and structural fold direction deviation coefficient; the fault aggregation region data includes the spatial location coordinates of the aggregation block, region boundary fitting parameters, and internal texture continuity index; the crack endpoint angular offset sequence includes the time-varying difference value of the angular offset, the two-dimensional projection coordinates of the endpoint, and the offset trend vector set; the enhanced texture feature map data includes the enhanced gradient magnitude matrix, the corrected texture direction vector map, and the interpolated region location mask; the evolutionary anomalous region includes the spatial boundary range of the anomalous region, the evolutionary trend classification label, and the local difference cumulative intensity index.

[0009] As a further aspect of the present invention, the specific steps of S1 are as follows: S101: Collect multi-temporal geomorphic images of the target monitoring area, call the preset spatial scale sequence parameters, map the multi-temporal geomorphic images to the corresponding resolution level, perform image cutting processing according to the grid division rules corresponding to each level, traverse each local region unit generated after cutting, extract the pixel matrix data inside it, index and mark the region units of multiple scales under the same temporal phase, and generate a multi-scale geomorphic image block set. S102: Call the multi-scale geomorphic image block set, identify the geometric abrupt change position of the edge pixel of each block unit, establish a local coordinate system centered on the abrupt change position, calculate the tilt angle of the edge tangent at the abrupt change point relative to the local coordinate axis, convert the tilt angle into a two-dimensional unit vector with direction attribute, select vectors with a modulus value greater than the noise suppression benchmark value as feature descriptors, and generate a block boundary structure folding direction vector set. S103: Based on the set of folding direction vectors of the block boundary structure, determine the adjacent objects of each block unit in the spatial topology, extract the direction vectors between the adjacent object pairs and perform dot product operation and inverse cosine transformation to obtain angle values, calculate the absolute value of the difference between the corresponding angle values ​​of adjacent blocks, arrange the calculated difference data in a two-dimensional matrix according to the spatial position index of the block, and construct the block main boundary angle difference array.

[0010] As a further aspect of the present invention, the specific steps of S2 are as follows: S201: Based on the block main boundary angle difference array, analyze the structural fold angle deviation values ​​between adjacent blocks, call the multi-scale landform image block set, obtain the edge gray matrix at the junction of adjacent blocks, use the gradient operator to calculate the gradient direction vector of the edge pixels, calculate the mean cosine similarity of adjacent gradient vectors, quantify the degree of matching of texture direction, and generate a block edge geometric and texture feature parameter set. S202: Call the block edge geometry and texture feature parameter set, call the preset angle filtering threshold and consistency filtering threshold, perform dual logic verification on the structural angle deviation value and texture direction matching degree in the parameter set, filter the block adjacent records with deviation values ​​less than the angle filtering threshold and matching degree greater than the consistency filtering threshold, and generate the target fracture area block index list. S203: Based on the target fault region block index list, extract the corresponding block image data from the multi-scale geomorphic image block set, calculate the coordinate offset of the feature control points on the shared boundary of adjacent blocks, construct a six-parameter affine transformation matrix, perform translation and rotation correction on the blocks, and stitch and fuse the corrected blocks in the spatial coordinate system to generate fault aggregation region data.

[0011] As a further aspect of the present invention, the specific steps of S3 are as follows: S301: Based on the fault aggregation area data, perform morphological refinement, extract the fracture skeleton, use the edge detection operator to track the pixel abrupt boundary of the fracture area, use the linear regression algorithm to fit the geometric extension trajectory of the fracture main axis segment, traverse the topological nodes of the trajectory, locate the pixel positions at the beginning and end, and generate the coordinate set of the fracture edge contour endpoints. S302: Call the coordinate set of the crack edge contour endpoints, construct a local analysis window centered on the endpoints in the topographic images of the first and second time phases, identify the tangent directions of the crack edges on both sides within the window, calculate the angle between the tangent directions, quantify the opening degree at the endpoints, extract and temporally correlate the opening degree values ​​in the two time phases, and generate a dual-time phase endpoint opening angle feature parameter set. S303: Based on the dual-phase endpoint sub-angle characteristic parameter set, establish a two-dimensional direction mapping field, convert the sub-angle values ​​of multiple phases into polar coordinate vectors in the direction field, calculate the rotation angle deviation and modulus scaling ratio of the vector in the time dimension, arrange and aggregate the deviation data in an orderly manner according to the extension direction of the crack principal axis, and generate the crack endpoint sub-angle offset sequence.

[0012] As a further aspect of the present invention, the specific steps of S4 are as follows: S401: For the rotation angle deviation data recorded in the crack endpoint angle offset sequence, a fixed step size sliding window is set to perform interval difference operation along the sequence index, calculate the change slope of the data within the window, quantify the evolution rate of the angle in the time dimension, and combine the rate value with the spatial coordinate information of the original sequence to construct the angle offset gradient vector set. S402: Based on the angular offset gradient vector set, call the multi-temporal landform image and use the structure tensor operator to calculate the local texture principal direction of each pixel in the image, construct a dense vector field representing the direction of landform texture, map the angular offset gradient vector set to the corresponding coordinate position of the vector field, calculate the projection component between the gradient vector and the texture principal direction vector point by point, construct a multi-dimensional data structure including spatial position, texture direction and projection weight, and generate a texture gradient projection mapping matrix; S403: Based on the texture gradient projection mapping matrix, analyze the geometric angle between the gradient vector and the main direction vector of the texture at each point, compare the angle value with the preset projection screening threshold, filter the target texture vector whose angle value meets the threshold restriction condition, perform weighted extended interpolation operation on the pixel area where the vector is located, fill the non-continuous gaps in texture information, and generate enhanced texture feature map data.

[0013] As a further aspect of the present invention, the projection filtering threshold is set by traversing the texture vector field and calculating the absolute value of the directional angle between the texture vector and the corresponding gradient derivative at each pixel position, constructing a set of directional angle numerical distributions across the entire field, performing statistical analysis on the set of directional angle numerical distributions across the entire field, calculating its arithmetic mean and standard deviation, and using the value obtained by subtracting the standard deviation by a preset multiple from the arithmetic mean as the projection filtering threshold.

[0014] As a further aspect of the present invention, the specific steps of S5 are as follows: S501: Call the enhanced texture feature map data, use the cross-correlation matching algorithm to calculate the pixel displacement vector between multiple temporal feature maps, perform geometric correction and resampling on the image data to eliminate spatial misalignment, calculate the pixel intensity difference matrix at the corresponding position after correction, and construct the registration residual data field. According to the preset physical size parameters, the data field is spatially discretized and cut to generate a registration residual grid cell set. S502: Traverse each independent cell in the registered residual grid cell set, extract the morphological contour of the response region using the edge operator, count the number of closed contours in the cell as the number of residual morphological response contour lines, calculate the number of alternations of curvature signs on the contour lines, determine the fluctuation frequency, measure the maximum projection length of the contour and calculate its ratio with the cell side length to obtain the maximum contour span ratio, combine the three indicators to generate a grid morphological response feature vector table; S503: Based on the grid morphological response feature vector table, calculate the numerical difference of the feature index in the continuous time series, map each grid cell to the preset evolution intensity category according to the difference amplitude, perform spatial superposition operation on the classification results under multiple time phases, count the cumulative occurrence of high intensity categories, delineate the connected grid regions with the cumulative occurrence exceeding the anomaly discrimination threshold as risk ranges, and generate evolution anomaly regions.

[0015] As a further aspect of the present invention, the method for setting the anomaly discrimination threshold is as follows: a set of monitoring images that have been manually verified and confirmed to be in a geologically stable state within a historical period is selected as a reference sample set. Registration residual calculation, grid division and feature extraction are performed on the reference sample set. The cumulative occurrence frequency of high-intensity categories in each grid unit within the corresponding time span is counted. A probability density distribution curve of the cumulative occurrence frequency is constructed. The expected value and standard deviation of the distribution curve are calculated. The sum of the expected value and the standard deviation of a preset multiple is used as the anomaly discrimination threshold.

[0016] A geological hazard detection system based on image processing, the system comprising: The boundary difference analysis module is used to acquire multi-temporal geomorphic images of the target monitoring area, divide the multi-temporal geomorphic images into image blocks according to the spatial scale sequence, and extract the structural fold vectors at the boundaries of the image blocks to construct the block main boundary angle difference array. The fault region identification module calculates the angle change value between adjacent image blocks based on the block main boundary angle difference array, filters image blocks that meet the angle and consistency thresholds and combines them affinely to generate fault aggregation region data. The fracture migration analysis module extracts the fracture edge contour line based on the fault aggregation area data, calculates the dual-temporal endpoint subtraction angle of the fracture edge contour line endpoints and maps it to the direction field, and constructs the fracture endpoint subtraction angle migration sequence. The texture enhancement processing module calculates the gradient derivation amount for the crack endpoint angular offset sequence, projects the gradient derivation amount to the texture field of the multi-temporal geomorphic image, and interpolates the texture vectors that meet the projection screening threshold to generate enhanced texture feature map data. The abnormal region identification module constructs a registration residual data field and divides it into grid windows based on the enhanced texture feature map data. It calculates the number of morphological contours, fluctuation frequency, and span ratio within the grid windows. Based on the difference response classification and superposition results of the number of morphological contours, fluctuation frequency, and span ratio, it determines the evolutionary abnormal region.

[0017] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by constructing a block-based main boundary angle difference array and combining it with edge gray-scale gradient consistency for affine aggregation, the accurate reconstruction and continuity restoration of cross-scale nonlinear fault structures in geomorphological images are achieved. The gradient derivative is generated by using the crack endpoint angular offset sequence and driving the texture vector field expansion, thereby directionally enhancing the extension characteristics of the micro-crack ends. Combined with the spatial morphological response raster analysis of the registration residual data field, it can keenly capture the subtle morphological changes of local landforms during the time evolution process, solving the problem of the difficulty in identifying the precursory features of minor disasters in complex backgrounds. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying 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.

[0019] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a detailed schematic diagram of S1 of the present invention; Figure 3 This is a detailed schematic diagram of S2 of the present invention; Figure 4 This is a detailed schematic diagram of S3 of the present invention; Figure 5 This is a detailed schematic diagram of S4 of the present invention; Figure 6 This is a detailed schematic diagram of S5 of the present invention; Figure 7 This is a system module diagram of the present invention. Detailed Implementation

[0020] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0021] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0022] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, their intended meanings are consistent. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, their intended meanings are consistent.

[0023] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0024] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0025] Please see Figure 1 This invention provides a geological hazard detection method based on image processing, comprising the following steps: S1: Collect multi-temporal geomorphic images of the target monitoring area, divide the multi-temporal geomorphic images into regions according to the spatial scale sequence to generate image blocks, extract the structural bending angle direction vectors located at the boundary of the image blocks, and construct the block main boundary angle difference array based on the structural bending angle direction vectors. S2: Based on the block principal boundary angle difference array, calculate the change value of the boundary principal vector angle between adjacent image blocks, filter image blocks whose boundary principal vector angle change value is lower than the angle filtering threshold and whose edge gray level continuous gradient direction consistency is higher than the consistency filtering threshold, and perform affine combination operation on image blocks that meet the filtering conditions to generate tomographic aggregation region data. S3: Extract the fracture edge contour line based on the fault aggregation area data, select the endpoint of the fracture edge contour line as the sampling reference to calculate the endpoint subtraction angle of the fracture main axis segment in the first and second time phases, map the endpoint subtraction angle to the two-dimensional direction field, and construct the fracture endpoint subtraction angle offset sequence. S4: Perform interval difference calculation on the crack endpoint angular offset sequence to generate gradient derivation quantity, project the gradient derivation quantity onto the texture vector field of the multi-temporal geomorphic image, and perform extended interpolation operation on the texture vector in the texture vector field where the angle between the gradient direction and the gradient derivation quantity is less than the projection screening threshold to generate enhanced texture feature map data. S5: Construct a registration residual data field based on enhanced texture feature map data, divide the registration residual data field into grid windows of fixed size, calculate the number of residual morphological response contour lines, fluctuation frequency, and maximum contour span ratio within the grid windows, and determine the evolutionary anomalous region based on the differential response classification and superposition results of the number of residual morphological response contour lines, fluctuation frequency, and maximum contour span ratio in multiple time slices.

[0026] The block-based main boundary angle difference array includes the block index identifier, boundary normal vector angle difference elements, and structural fold direction deviation coefficient; the fault aggregation region data includes the spatial location coordinates of the aggregation block, region boundary fitting parameters, and internal texture continuity index; the crack endpoint angular offset sequence includes the time-varying difference value of the angular offset, the two-dimensional projection coordinates of the endpoint, and the offset trend vector set; the enhanced texture feature map data includes the enhanced gradient magnitude matrix, the corrected texture direction vector map, and the interpolated region location mask; the evolutionary anomalous region includes the spatial boundary range of the anomalous region, the evolutionary trend classification label, and the local difference cumulative intensity index.

[0027] Preferably, the various thresholds involved in the embodiments of the present invention (such as angle screening threshold, consistency screening threshold, projection screening threshold, and anomaly detection threshold) support adaptive adjustment based on statistics. Specifically, during the initialization phase, the system analyzes historical background images of the target monitoring area (i.e., images from geologically stable periods), calculates the statistical distribution histogram of relevant feature parameters (such as the angle between the texture vector and the gradient direction), and uses the boundary values ​​of a preset confidence interval (e.g., a 95% confidence interval) as the baseline threshold. Furthermore, as the surface vegetation cover changes due to seasonal variations, the system can periodically update the aforementioned statistical distribution parameters, thereby dynamically adjusting the threshold to reduce the false alarm rate caused by environmental changes.

[0028] Please see Figure 2 The specific steps of S1 are as follows: S101: Collect multi-temporal geomorphic images of the target monitoring area, call the preset spatial scale sequence parameters, map the multi-temporal geomorphic images to the corresponding resolution levels, perform image cutting processing according to the grid division rules corresponding to each level, traverse each local region unit generated after cutting, extract the pixel matrix data inside it, index and mark the region units of multiple scales under the same temporal phase, and generate a multi-scale geomorphic image block set.

[0029] Image acquisition equipment captures high-resolution optical remote sensing images in key monitoring areas prone to landslides or collapses during different months before and after the rainy season, forming a geomorphic image data set that includes a temporal dimension. The image processing server reads pre-set spatial scale sequence parameters, for example, setting the first-level scale to 1:500, the second-level scale to 1:2000, and the third-level scale to 1:5000, and resamples the original multi-temporal geomorphic images to the aforementioned three resolution levels. For the first-level high-resolution images, a non-overlapping grid size of 50 pixels by 50 pixels is used for segmentation; for the second-level medium-resolution images, a grid size of 100 pixels by 100 pixels is used; and for the third-level medium-resolution images, a grid size of 250 pixels by 250 pixels is used. The image processing server scans each rectangular local region unit generated after segmentation, reading the grayscale value and color channel value of each pixel within the region unit to construct a pixel matrix. Subsequently, the image processing server assigns a unique spatial location code to each regional unit. This code is composed of a time identifier, a scale level identifier, and a row and column number. For example, the slice located in the tenth row and fifteenth column of the first-level image collected in June is marked as a specific index code, and the corresponding slice of the same geographical location in the August image is marked as an associated index code. This completes the data structure organization of all slices and generates a multi-scale geomorphic image block set.

[0030] S102: Call the multi-scale geomorphic image block set, identify the geometric abrupt change location of the edge pixels of each block unit, establish a local coordinate system centered on the abrupt change location, calculate the tilt angle of the edge tangent at the abrupt change point relative to the local coordinate axis, convert the tilt angle into a two-dimensional unit vector with direction attribute, select vectors with a magnitude greater than the noise suppression benchmark value as feature descriptors, and generate a block boundary structure folding direction vector set.

[0031] For each terrain image block, the image processing server uses an edge detection operator to scan the rock texture and terrain contours within the block, locating edge paths where pixel grayscale values ​​change drastically. Along these edge paths, curvature changes are detected point-by-point, and locations with curvature radii less than a preset inflection point radius threshold are identified as geometric abrupt change locations, i.e., the bends in the rock structure. A local Cartesian coordinate system is constructed with this abrupt change location as the origin, the horizontal edge direction of the image block as the x-axis, and the vertical edge direction as the y-axis. The image processing server selects edge pixels in the neighborhood of the abrupt change location and uses the least squares method to fit a tangent line passing through that point, calculating the angle between this tangent line and the horizontal axis of the local coordinate system. For example, if the angle between the fitted tangent line and the positive horizontal axis is 45 degrees, it is converted into a two-dimensional unit vector containing horizontal and vertical components. At this point, minor pseudo-features caused by image noise or vegetation cover need to be removed. The image processing server calculates the magnitude of the resulting vector and sets a noise suppression baseline value. The noise suppression benchmark is set by selecting images with flat, featureless regions for edge extraction, calculating the average value of the random vector magnitudes generated by the background noise, and using three times this average value as the noise suppression benchmark. If the currently calculated vector magnitude is less than the benchmark, it is discarded; if it is greater than the benchmark, the vector is retained and its corresponding spatial coordinates and direction attributes are recorded to generate a set of folding direction vectors for the block boundary structure.

[0032] S103: Based on the set of folding direction vectors of the block boundary structure, determine the adjacent objects of each block unit in the spatial topology, extract the direction vectors between adjacent object pairs and perform dot product operation and inverse cosine transformation to obtain angle values, calculate the absolute value of the difference between the corresponding angle values ​​of adjacent blocks, arrange the calculated difference data in a two-dimensional matrix according to the spatial position index of the block, and construct the block main boundary angle difference array.

[0033] The image processing server identifies adjacent block units in the four directions (top, bottom, left, and right) based on the spatial index encoding of the block units, establishing adjacency pairs. For each pair of adjacent blocks, it retrieves the most significant structural angle direction vector at the boundary from the block boundary structural angle direction vector set. For example, the principal vector of the current block boundary is vector A, and the principal vector of the adjacent boundary of the right-side adjacent block is vector B. The image processing server performs a dot product operation on vectors A and B, divides the result of the dot product by the product of the magnitudes of vectors A and B to obtain a cosine value, and then performs an inverse cosine transform on this cosine value to obtain the angle between the two vectors. Assuming the angle corresponding to vector A is 30 degrees and the angle corresponding to vector B is 35 degrees, the absolute value of the difference between the two angle values ​​is 5 degrees. The image processing server traverses all adjacent block pairs within the entire monitoring area, sequentially performing the above angle difference calculation, and fills the calculated difference data into the matrix unit corresponding to the original image grid distribution. For example, if the monitoring area is divided into a grid of 100 rows and 100 columns, an array of 100 rows and 100 columns is constructed, and the angle difference between the block in the i-th row and j-th column and the block in the i-th row and j-th column plus one is stored in the corresponding position of the array, thereby forming an array of block main boundary angle differences that reflects the continuous changes in the geomorphic structure of the entire area.

[0034] Please see Figure 3 The specific steps of S2 are as follows: S201: Based on the block main boundary angle difference array, the structural fold angle deviation values ​​between adjacent blocks are parsed, the multi-scale geomorphic image block set is called to obtain the edge gray matrix at the junction of adjacent blocks, the gradient direction vector of the edge pixels is calculated using the gradient operator, and the cosine similarity mean of adjacent gradient vectors is calculated to quantify the degree of matching of texture direction, and generate the block edge geometric and texture feature parameter set.

[0035] The image processing server reads each element from the array of angle differences between the main boundaries of the blocks, which represents the structural fold angle deviation between adjacent blocks. Simultaneously, the server extracts image data of the corresponding junctions between adjacent blocks from the multi-scale geomorphic image block set based on the index, constructing an edge grayscale matrix. For each pixel in this edge grayscale matrix, the Sobel or Prewitt operator is applied to calculate its grayscale gradient components in the horizontal and vertical directions, respectively, thus synthesizing a gradient direction vector. The image processing server selects pairs of pixels within a preset width range on both sides of the junction boundary line and calculates the cosine similarity of the gradient direction vectors of each pair of pixels. For example, if the gradient direction of the left pixel is 60 degrees and the gradient direction of the corresponding pixel on the right is 62 degrees, then their cosine similarity is close to one. The server accumulates the cosine similarities of all pairs of pixels on the boundary line and calculates the arithmetic mean. This average value represents the degree of matching of the texture direction of adjacent blocks at the junction. If the texture is continuous, the mean value approaches one; if there is a break or misalignment, the mean value decreases significantly. The image processing server packages and stores the parsed structural folding deviation values ​​and the calculated texture direction matching values ​​to generate a set of block edge geometry and texture feature parameters.

[0036] S202: Call the block edge geometry and texture feature parameter set, call the preset angle filtering threshold and consistency filtering threshold, perform dual logical verification on the structural fold angle deviation value and texture direction matching degree in the parameter set, filter the block adjacent records with deviation value less than the angle filtering threshold and matching degree greater than the consistency filtering threshold, and generate the target fracture area block index list.

[0037] The specific process for setting the angle screening threshold is as follows: select the edge segments of the continuous fault structure that have been confirmed in the historical geological disaster samples as reference objects, calculate the statistical distribution characteristics of the structural bending angle change rate at adjacent small segments along the extension direction of the reference object, and select the upper limit boundary value covering the main peak interval of the statistical distribution characteristics as the angle screening threshold.

[0038] The process of setting the consistency screening threshold is as follows: by sampling the texture gradient of the non-fractured intact rock surface area of ​​the geomorphological image in the monitoring area, calculating the average level of the gradient direction consistency coefficient between local pixel neighborhoods in the intact area, and subtracting the tolerance margin based on the rock surface roughness estimation from the average level, the value is used as the consistency screening threshold.

[0039] The image processing server first determines the screening criteria. For the angle screening threshold, historical geological images of known fault slip are selected, and multiple small segments of the continuous fault edges are extracted. The rate of change of the angle between each adjacent segment is calculated. Statistical data shows that 95% of the rate of change of the angle of continuous fault edges is concentrated between 0 and 15 degrees, exhibiting a normal distribution. Therefore, the upper limit of this main peak range, 15 degrees, is selected as the angle screening threshold. For the consistency screening threshold, images of intact granite rock surfaces that have not undergone fracturing are selected for analysis, and the average value of the consistency of their local texture gradient direction is calculated to be 0.9. Considering the roughness effect caused by natural weathering of the rock surface, the estimated tolerance margin is 0.1. Therefore, 0.8, obtained by subtracting 0.1 from 0.9, is set as the consistency screening threshold. During the screening process, the image processing server reads records from the parameter set one by one. For example, if the structural angle deviation of an adjacent segment is 5 degrees, the texture matching degree is 0.85. Since the deviation is less than 15 degrees and the degree of agreement is greater than 0.85, the record passes the verification; conversely, if the deviation is 20 degrees or the degree of agreement is 0.5, it is determined that the continuity condition is not met. The image processing server extracts the index numbers of all block pairs that pass the double verification and summarizes them to generate a block index list of the target fracture region.

[0040] S203: Based on the target fault region block index list, extract the corresponding block image data from the multi-scale geomorphic image block set, calculate the coordinate offset of the feature control points on the shared boundary of adjacent blocks, construct a six-parameter affine transformation matrix, perform translation and rotation correction on the blocks, and stitch and fuse the corrected blocks in the spatial coordinate system to generate fault aggregation region data.

[0041] The image processing server uses an index list to locate image blocks belonging to potential fault zones and identifies salient feature points on shared boundaries, such as rock protrusions or vegetation roots, as control points. The positional differences of these control points in the coordinate systems of adjacent blocks are calculated. For example, if a control point's coordinates in the left block are (100, 100), and its corresponding coordinates in the right block are (98, 102), then there is a horizontal offset of -2 and a vertical offset of +2. Based on at least three sets of control point coordinate pairs, a six-parameter affine transformation matrix, including horizontal translation, vertical translation, rotation angle, horizontal scaling, vertical scaling, and shear parameters, is calculated using the least squares method. This matrix is ​​applied to transform the coordinates of each pixel in the right block, ensuring precise spatial alignment with the left block. After geometric correction, the image processing server uses a weighted average fusion algorithm to process the pixel grayscale values ​​in overlapping areas, eliminating stitching gaps and seamlessly merging the scattered image blocks into a complete local image of the fault zone, generating fault aggregation region data.

[0042] Please see Figure 4 The specific steps of S3 are as follows: S301: Based on the fault aggregation region data, perform morphological refinement, extract the fracture skeleton, use the edge detection operator to track the pixel abrupt boundary of the fracture region, use the linear regression algorithm to fit the geometric extension trajectory of the fracture main axis segment, traverse the topological nodes of the trajectory, locate the pixel positions at the beginning and end, and generate the coordinate set of the fracture edge contour endpoints.

[0043] The image processing server binarizes the fault aggregation area data to distinguish the fracture area from the background rock mass. Then, a morphological thinning algorithm is applied to continuously peel away the edge pixels of the fracture area until only a single-pixel-width central skeleton line remains. The Canny edge detection operator is used to search for abrupt gray-level changes along both sides of the skeleton line to establish the solid boundary of the fracture. A linear regression algorithm is then used to fit the thinned skeleton line pixel set. For example, the horizontal and vertical coordinates of all pixels on the skeleton line are substituted into the linear equation, and the slope and intercept of the best-fit line are obtained by minimizing the sum of squared errors. This line represents the extension trajectory of the fracture's main axis. The image processing server scans along this trajectory, identifying endpoints with only one adjacent pixel, i.e., the start and end points of the fracture. The precise row and column coordinates of these endpoints in the global coordinate system are recorded; for example, the start coordinates are (500, 600), and the end coordinates are (550, 680). All identified endpoint coordinates are compiled to generate a set of fracture edge contour endpoint coordinates.

[0044] S302: Call the coordinate set of the crack edge contour endpoints, construct a local analysis window centered on the endpoints in the topographic images of the first and second time phases, identify the tangent directions on both sides of the crack within the window, calculate the angle between the tangent directions, quantify the opening degree at the endpoints, extract and temporally correlate the opening degree values ​​in the two time phases, and generate a dual-time phase endpoint opening angle feature parameter set.

[0045] In the first time-phase image acquired in June and the second time-phase image acquired in August, local analysis windows of 20 pixels by 20 pixels were extracted. Within each window, the rock wall edges on both sides of the fracture terminus were identified, and tangents were fitted to the edges. The angle formed by the two tangents at the endpoints was calculated. For example, in the first time-phase image, the opening angle of a fracture endpoint was measured to be 10 degrees; in the second time-phase image, the opening angle of the same endpoint was measured to be 12 degrees. The image processing server paired and recorded these two values, adding time labels to clarify the evolution of the opening angle over time. The same operation was performed on all fracture endpoints within the monitoring area to obtain a massive amount of opening angle pair data, generating a set of dual-time-phase endpoint opening angle feature parameters.

[0046] S303: Based on the dual-phase endpoint subtraction angle characteristic parameter set, a two-dimensional direction mapping field is established, the subtraction angle values ​​of multiple phases are converted into polar coordinate vectors in the direction field, the rotation angle deviation and magnitude scaling ratio of the vector in the time dimension are calculated, and the deviation data are arranged and assembled in an orderly manner according to the extension direction of the crack principal axis to generate the crack endpoint subtraction angle offset sequence.

[0047] The image processing server constructs a two-dimensional polar coordinate system as a direction mapping field. The angular displacement value of the first time phase is mapped to the magnitude of a vector, and the direction of the crack principal axis is mapped to the polar angle of the vector; similarly, the data of the second time phase is mapped. For example, the vector magnitude of the first time phase is 10, and the polar angle is 45 degrees; the vector magnitude of the second time phase is 12, and the polar angle is 47 degrees. The image processing server calculates the rotation angle deviation between the two vectors, i.e., 47 degrees minus 45 degrees equals 2 degrees; it also calculates the magnitude scaling ratio, i.e., 12 divided by 10 equals 1.2. This means that a slight rotational distortion occurs at the crack endpoints while they are opening. Based on the order of the crack principal axis from start to finish, the server sequentially arranges the deviation data of each endpoint and key node along the path, forming a data sequence reflecting the overall dynamic deformation trend of the crack, generating a crack endpoint angular displacement sequence.

[0048] Please see Figure 5 The specific steps of S4 are as follows: S401: For the rotation angle deviation data recorded in the crack endpoint angle offset sequence, a fixed step size sliding window is set to perform interval difference operation along the sequence index, calculate the slope of the data change within the window, quantify the evolution rate of the angle in the time dimension, and combine the rate value with the spatial coordinate information of the original sequence to construct the angle offset gradient vector set.

[0049] The image processing server sets a sliding window of length five, which slides progressively along the angular offset sequence at the crack endpoints. At each window position, the rotation angle deviation values ​​of the five consecutive nodes contained within the window are extracted. The first-order difference of these five values ​​is calculated, and the average value is obtained to obtain the slope of change in that local interval. For example, if the deviation values ​​of the five nodes are 2 degrees, 2.1 degrees, 2.2 degrees, 2.3 degrees, and 2.4 degrees respectively, the average difference value is 0.1 degrees, which represents the cumulative rate of angular offset in spatial extension. The image processing server assigns a direction attribute (along the crack tangential) to this slope value and, combined with the spatial coordinates of the center node of the window in the original image, constructs a gradient vector with position, magnitude, and direction. After traversing the entire sequence, a set of angular offset gradient vectors is generated.

[0050] S402: Based on the angular offset gradient vector set, multi-temporal landform images are called and the local texture principal direction of each pixel in the image is calculated using the structure tensor operator. A dense vector field representing the direction of landform texture is constructed. The angular offset gradient vector set is mapped to the corresponding coordinate position of the vector field. The projection component between the gradient vector and the texture principal direction vector is calculated point by point. A multi-dimensional data structure including spatial position, texture direction and projection weight is constructed to generate the texture gradient projection mapping matrix.

[0051] The image processing server first calculates the structure tensor for the multi-temporal geomorphic image, and solves for the eigenvalues ​​and eigenvectors of the structure tensor matrix for each pixel. The eigenvector corresponding to the largest eigenvalue is defined as the principal direction of the texture at that point, thus establishing a dense vector field covering the entire image. Subsequently, each vector in the angular offset gradient vector set generated in S401 is projected onto this dense vector field according to its coordinates. For each projection point, the dot product between the angular offset gradient vector and the principal direction vector of the base image texture at that point is calculated. This dot product result is the projection component, which physically represents the projection intensity of the crack opening deformation trend along the rock texture direction. For example, if the gradient vector direction is parallel to the texture direction, the projection component is the largest, indicating a very high probability that the crack will extend along the texture. The image processing server combines and stores the pixel coordinates, the principal direction vector of the texture, and the calculated projection component values ​​to generate a texture gradient projection mapping matrix.

[0052] S403: Based on the texture gradient projection mapping matrix, the geometric angle between the gradient vector and the main direction vector of the texture at each point is analyzed. The angle value is compared with the preset projection screening threshold. Target texture vectors whose angle values ​​meet the threshold restriction conditions are selected. Weighted extended interpolation operation is performed on the pixel area where the vector is located to fill the discontinuous gaps in texture information and generate enhanced texture feature map data.

[0053] The projection filtering threshold is set by traversing the texture vector field and calculating the absolute value of the directional angle between the texture vector and the corresponding gradient derivative at each pixel position, constructing a set of directional angle values ​​for the entire field, performing statistical analysis on the set of directional angle values ​​for the entire field, calculating its arithmetic mean and standard deviation, and using the value obtained by subtracting the standard deviation by a preset multiple from the arithmetic mean as the projection filtering threshold.

[0054] The image processing server first determines the projection screening threshold. By calculating the absolute value of the angle between all texture vectors and gradient derivatives in the entire image, an angle distribution set is obtained. The arithmetic mean of this set is set to 60 degrees, the standard deviation to 15 degrees, and the preset multiplier to 2; therefore, the projection screening threshold is calculated as 60 degrees minus 30 degrees, which is 30 degrees. Next, the server traverses the mapping matrix, filtering out texture vectors with angles less than 30 degrees. For the regions containing these vectors that meet the criteria, it indicates that the texture direction is highly consistent with the crack propagation trend, representing potential disaster hazard points. Using these points as centers, the image processing server employs a Gaussian weighted interpolation algorithm, referencing neighborhood texture features, to fill and enhance areas with interrupted or blurred textures, making the originally inconspicuous micro-crack textures stand out in the feature map, generating enhanced texture feature map data.

[0055] Please see Figure 6 The specific steps of S5 are as follows: S501: Call the enhanced texture feature map data, use the cross-correlation matching algorithm to calculate the pixel displacement vector between multiple temporal feature maps, perform geometric correction and resampling on the image data to eliminate spatial misalignment, calculate the pixel intensity difference matrix at the corresponding position after correction, and construct the registration residual data field. According to the preset physical size parameters, the data field is spatially discretized and cut to generate the registration residual grid cell set.

[0056] The image processing server uses the enhanced first-phase texture feature map as a reference and the second-phase feature map as a floating image. A normalized cross-correlation matching algorithm is employed to search for the patch in the floating image most similar to the reference image, calculating the displacement vector for each pixel. Based on the displacement vector, the second-phase image undergoes reverse distortion correction to ensure precise overlap with the reference image. Subsequently, the pixel values ​​of the reference image are subtracted from the corresponding pixel values ​​of the corrected image to obtain a difference matrix. The non-zero values ​​in this difference matrix represent the local non-rigid deformation residuals, excluding the overall translation. The image processing server constructs a registration residual data field and divides it into numerous small grid cells according to, for example, actual ground physical dimensions of 10 meters by 10 meters, generating a registration residual grid cell set.

[0057] S502: Traverse each independent cell in the registered residual grid cell set, extract the morphological contour of the response region using the edge operator, count the number of closed contours in the cell as the number of residual morphological response contour lines, calculate the number of alternations of curvature signs on the contour lines to determine the fluctuation frequency, measure the maximum projection length of the contour and calculate its ratio with the cell side length to obtain the maximum contour span ratio, combine the three indicators to generate a grid morphological response feature vector table.

[0058] The image processing server extracts contour lines from the residual data within each grid cell to generate a morphological profile. It counts the number of closed contour loops within a cell; for example, if three independent closed loops are detected, the count is three. The curvature is calculated along each contour line, and the number of times the curvature changes from positive to negative or vice versa is counted; for example, if a contour line changes its curvature direction five times, the frequency of change is five. The maximum horizontal span of the contour line is measured and set to eight pixels, while the grid cell side length is ten pixels, resulting in a maximum contour span ratio of 0.8. The image processing server uses this set of data—"count three, frequency five, ratio 0.8"—as a feature vector, associates it with the index of that grid cell, and after traversing all cells, generates a grid point morphological response feature vector table.

[0059] S503: Based on the grid morphological response feature vector table, calculate the numerical difference of feature indicators in the continuous time series, map each grid cell to the preset evolution intensity category according to the difference amplitude, perform spatial superposition operation on the classification results under multiple time phases to count the cumulative occurrence of high intensity categories, delineate the connected grid regions with the cumulative occurrence exceeding the anomaly discrimination threshold as risk ranges, and generate evolution anomaly regions.

[0060] The specific method for setting the anomaly discrimination threshold is as follows: Select monitoring images that have been manually verified and confirmed to be in a geologically stable state within a historical period as the benchmark sample set. Perform registration residual calculation, grid division and feature extraction on the benchmark sample set. Count the cumulative occurrence of high intensity categories in each grid unit within the corresponding time span, construct the probability density distribution curve of the cumulative occurrence, calculate the expected value and standard deviation of the distribution curve, and use the sum of the expected value and the standard deviation of the preset multiple as the anomaly discrimination threshold.

[0061] The image processing server first sets an anomaly detection threshold. Image data from stable areas that have not experienced disasters in the past five years are selected, and the cumulative occurrence count is calculated following the steps described above. Statistical results show that the cumulative occurrence count of high-intensity categories in stable area grids follows a normal distribution with a mathematical expectation of 2 and a standard deviation of 1. A preset multiplier of 3 is set, so the anomaly detection threshold is 2 + 3, or five times. During real-time monitoring, the server compares the changes in the characteristic indicators of the current grid unit over time. If the change is large, it is classified as "high-intensity evolution." The number of times the grid has been classified as "high-intensity evolution" in the last ten monitoring periods is counted. If the cumulative occurrence count of a grid reaches six times, exceeding the threshold of five times, the grid is determined to be an anomaly. The image processing server merges all spatially connected grid areas that are determined to be anomalies, delineates the specific geographical range, and finally generates an evolution anomaly area.

[0062] Please see Figure 7 A geological hazard detection system based on image processing includes: The boundary difference analysis module is used to acquire multi-temporal geomorphic images of the target monitoring area, divide the multi-temporal geomorphic images into image blocks according to the spatial scale sequence, and extract the structural fold vectors at the boundaries of the image blocks to construct the block main boundary angle difference array. The fault region identification module calculates the angle change value between adjacent image blocks based on the block main boundary angle difference array, filters image blocks that meet the angle and consistency thresholds and combines them affinely to generate fault aggregation region data. The fracture migration analysis module extracts the fracture edge contour line based on the fault aggregation area data, calculates the dual-temporal endpoint subtraction angle of the fracture edge contour line endpoints and maps it to the direction field, and constructs the fracture endpoint subtraction angle migration sequence. The texture enhancement processing module calculates the gradient derivation amount for the crack endpoint angular offset sequence, projects the gradient derivation amount to the texture field of the multi-temporal geomorphic image, and interpolates the texture vectors that meet the projection screening threshold to generate enhanced texture feature map data. The abnormal region identification module constructs a registration residual data field based on enhanced texture feature map data and divides it into grid windows. It calculates the number of morphological contours, fluctuation frequency, and span ratio within the grid windows. Based on the difference response classification and superposition results of the number of morphological contours, fluctuation frequency, and span ratio, it determines the evolutionary abnormal region.

[0063] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A geological hazard detection method based on image processing, characterized in that, Includes the following steps: S1: Collect multi-temporal geomorphic images of the target monitoring area, divide the multi-temporal geomorphic images into images according to spatial scale sequence and generate image blocks, extract the structural folding vectors at the boundaries of the image blocks to construct the block main boundary angle difference array; S2: Based on the block main boundary angle difference array, calculate the angle change value between adjacent image blocks, filter image blocks that meet the angle and consistency thresholds and combine them affinely to generate tomographic aggregation region data; S3: Extract the fracture edge contour line based on the fault aggregation area data, calculate the dual-temporal endpoint subtraction angle of the fracture edge contour line endpoints and map it to the direction field to construct the fracture endpoint subtraction angle offset sequence. S4: For the crack endpoint angular offset sequence, calculate the gradient derivation amount, project the gradient derivation amount to the texture field of the multi-temporal geomorphic image, interpolate the texture vectors that meet the projection screening threshold, and generate enhanced texture feature map data. S5: Based on the enhanced texture feature map data, construct a registration residual data field and divide it into grid windows. Calculate the number of morphological contours, fluctuation frequency, and span ratio within the grid windows. Determine the evolutionary anomalous region based on the difference response classification and superposition results of the number of morphological contours, fluctuation frequency, and span ratio.

2. The geological disaster detection method based on image processing according to claim 1, characterized in that, The block-based main boundary angle difference array includes block index identifiers, boundary normal vector angle difference elements, and structural fold direction deviation coefficients; the fault aggregation region data includes aggregation block spatial location coordinates, region boundary fitting parameters, and internal texture continuity index; the crack endpoint angular offset sequence includes angular time-varying difference values, endpoint two-dimensional projection coordinates, and offset trend vector set; the enhanced texture feature map data includes enhanced gradient magnitude matrix, corrected texture direction vector map, and interpolated region location mask; the evolutionary anomalous region includes anomalous region spatial boundary range, evolutionary trend classification label, and local difference cumulative intensity index.

3. The geological disaster detection method based on image processing according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Collect multi-temporal geomorphic images of the target monitoring area, call the preset spatial scale sequence parameters, map the multi-temporal geomorphic images to the corresponding resolution level, perform image cutting processing according to the grid division rules corresponding to each level, traverse each local region unit generated after cutting, extract the pixel matrix data inside it, index and mark the region units of multiple scales under the same temporal phase, and generate a multi-scale geomorphic image block set. S102: Call the multi-scale geomorphic image block set, identify the geometric abrupt change position of the edge pixel of each block unit, establish a local coordinate system centered on the abrupt change position, calculate the tilt angle of the edge tangent at the abrupt change point relative to the local coordinate axis, convert the tilt angle into a two-dimensional unit vector with direction attribute, select vectors with a modulus value greater than the noise suppression benchmark value as feature descriptors, and generate a block boundary structure folding direction vector set. S103: Based on the set of folding direction vectors of the block boundary structure, determine the adjacent objects of each block unit in the spatial topology, extract the direction vectors between the adjacent object pairs and perform dot product operation and inverse cosine transformation to obtain angle values, calculate the absolute value of the difference between the corresponding angle values ​​of adjacent blocks, arrange the calculated difference data in a two-dimensional matrix according to the spatial position index of the block, and construct the block main boundary angle difference array.

4. The geological disaster detection method based on image processing according to claim 1, characterized in that, The specific steps of S2 are as follows: S201: Based on the block main boundary angle difference array, analyze the structural fold angle deviation values ​​between adjacent blocks, call the multi-scale landform image block set, obtain the edge gray matrix at the junction of adjacent blocks, use the gradient operator to calculate the gradient direction vector of the edge pixels, calculate the mean cosine similarity of adjacent gradient vectors, quantify the degree of matching of texture direction, and generate a block edge geometric and texture feature parameter set. S202: Call the block edge geometry and texture feature parameter set, call the preset angle filtering threshold and consistency filtering threshold, perform dual logic verification on the structural angle deviation value and texture direction matching degree in the parameter set, filter the block adjacent records with deviation values ​​less than the angle filtering threshold and matching degree greater than the consistency filtering threshold, and generate the target fracture area block index list. S203: Based on the target fault region block index list, extract the corresponding block image data from the multi-scale geomorphic image block set, calculate the coordinate offset of the feature control points on the shared boundary of adjacent blocks, construct a six-parameter affine transformation matrix, perform translation and rotation correction on the blocks, and stitch and fuse the corrected blocks in the spatial coordinate system to generate fault aggregation region data.

5. The geological hazard detection method based on image processing according to claim 1, characterized in that, The specific steps for S3 are as follows: S301: Based on the fault aggregation area data, perform morphological refinement, extract the fracture skeleton, use the edge detection operator to track the pixel abrupt boundary of the fracture area, use the linear regression algorithm to fit the geometric extension trajectory of the fracture main axis segment, traverse the topological nodes of the trajectory, locate the pixel positions at the beginning and end, and generate the coordinate set of the fracture edge contour endpoints. S302: Call the coordinate set of the crack edge contour endpoints, construct a local analysis window centered on the endpoints in the topographic images of the first and second time phases, identify the tangent directions of the crack edges on both sides within the window, calculate the angle between the tangent directions, quantify the opening degree at the endpoints, extract and temporally correlate the opening degree values ​​in the two time phases, and generate a dual-time phase endpoint opening angle feature parameter set. S303: Based on the dual-phase endpoint sub-angle characteristic parameter set, establish a two-dimensional direction mapping field, convert the sub-angle values ​​of multiple phases into polar coordinate vectors in the direction field, calculate the rotation angle deviation and modulus scaling ratio of the vector in the time dimension, arrange and aggregate the deviation data in an orderly manner according to the extension direction of the crack principal axis, and generate the crack endpoint sub-angle offset sequence.

6. The geological disaster detection method based on image processing according to claim 1, characterized in that, The specific steps of S4 are as follows: S401: For the rotation angle deviation data recorded in the crack endpoint angle offset sequence, set a fixed step size sliding window to perform interval difference operation along the sequence index, calculate the change slope of the data in the window, quantify the evolution rate of the angle in the time dimension, and combine the rate value with the spatial coordinate information of the original sequence to construct the angle offset gradient vector set. S402: Based on the angular offset gradient vector set, call the multi-temporal landform image and use the structure tensor operator to calculate the local texture principal direction of each pixel in the image, construct a dense vector field representing the direction of landform texture, map the angular offset gradient vector set to the corresponding coordinate position of the vector field, calculate the projection component between the gradient vector and the texture principal direction vector point by point, construct a multi-dimensional data structure including spatial position, texture direction and projection weight, and generate a texture gradient projection mapping matrix; S403: Based on the texture gradient projection mapping matrix, analyze the geometric angle between the gradient vector and the main direction vector of the texture at each point, compare the angle value with the preset projection screening threshold, filter the target texture vector whose angle value meets the threshold restriction condition, perform weighted extended interpolation operation on the pixel area where the vector is located, fill the non-continuous gaps in texture information, and generate enhanced texture feature map data.

7. The geological disaster detection method based on image processing according to claim 6, characterized in that, The projection filtering threshold is set by traversing the texture vector field and calculating the absolute value of the directional angle between the texture vector and the corresponding gradient derivative at each pixel position, constructing a set of directional angle values ​​for the entire field, performing statistical analysis on the set of directional angle values ​​for the entire field, calculating its arithmetic mean and standard deviation, and using the value obtained by subtracting the standard deviation by a preset multiple from the arithmetic mean as the projection filtering threshold.

8. The geological disaster detection method based on image processing according to claim 1, characterized in that, The specific steps of S5 are as follows: S501: Call the enhanced texture feature map data, use the cross-correlation matching algorithm to calculate the pixel displacement vector between multiple temporal feature maps, perform geometric correction and resampling on the image data to eliminate spatial misalignment, calculate the pixel intensity difference matrix at the corresponding position after correction, and construct the registration residual data field. According to the preset physical size parameters, the data field is spatially discretized and cut to generate a registration residual grid cell set. S502: Traverse each independent cell in the registered residual grid cell set, extract the morphological contour of the response region using the edge operator, count the number of closed contours in the cell as the number of residual morphological response contour lines, calculate the number of alternations of curvature signs on the contour lines, determine the fluctuation frequency, measure the maximum projection length of the contour and calculate its ratio with the cell side length to obtain the maximum contour span ratio, combine the three indicators to generate a grid morphological response feature vector table; S503: Based on the grid morphological response feature vector table, calculate the numerical difference of the feature index in the continuous time series, map each grid cell to the preset evolution intensity category according to the difference amplitude, perform spatial superposition operation on the classification results under multiple time phases, count the cumulative occurrence of high intensity categories, delineate the connected grid regions with the cumulative occurrence exceeding the anomaly discrimination threshold as risk ranges, and generate evolution anomaly regions.

9. The geological disaster detection method based on image processing according to claim 8, characterized in that, The specific method for setting the anomaly discrimination threshold is as follows: Select monitoring images that have been manually verified and confirmed to be in a geologically stable state within a historical period as a reference sample set. Perform registration residual calculation, grid division and feature extraction on the reference sample set. Count the cumulative occurrence of high-intensity categories in each grid unit within the corresponding time span. Construct a probability density distribution curve of the cumulative occurrence. Calculate the expected value and standard deviation of the distribution curve. Use the sum of the expected value and the standard deviation of a preset multiple as the anomaly discrimination threshold.

10. A geological hazard detection system based on image processing, characterized in that, The system is used to implement the image processing-based geological disaster detection method according to any one of claims 1-9, the system comprising: The boundary difference analysis module is used to acquire multi-temporal geomorphic images of the target monitoring area, divide the multi-temporal geomorphic images into image blocks according to the spatial scale sequence, and extract the structural fold vectors at the boundaries of the image blocks to construct the block main boundary angle difference array. The fault region identification module calculates the angle change value between adjacent image blocks based on the block main boundary angle difference array, filters image blocks that meet the angle and consistency thresholds and combines them affinely to generate fault aggregation region data. The fracture migration analysis module extracts the fracture edge contour line based on the fault aggregation area data, calculates the dual-temporal endpoint subtraction angle of the fracture edge contour line endpoints and maps it to the direction field, and constructs the fracture endpoint subtraction angle migration sequence. The texture enhancement processing module calculates the gradient derivation amount for the crack endpoint angular offset sequence, projects the gradient derivation amount to the texture field of the multi-temporal geomorphic image, and interpolates the texture vectors that meet the projection screening threshold to generate enhanced texture feature map data. The abnormal region identification module constructs a registration residual data field and divides it into grid windows based on the enhanced texture feature map data. It calculates the number of morphological contours, fluctuation frequency, and span ratio within the grid windows. Based on the difference response classification and superposition results of the number of morphological contours, fluctuation frequency, and span ratio, it determines the evolutionary abnormal region.

Images