A geological disaster risk assessment data processing method and system

Abstract

The application provides a geological disaster risk assessment data processing method and system, and relates to the technical field of geological disasters, which comprises the following steps: obtaining original observation data of an evaluation area, generating a plurality of error measurement indexes according to the original observation data; mapping the original observation data and the error measurement indexes to a regular grid to obtain a grid field, and constructing a geographic analysis block; constructing a disaster-causing state variable based on the geographic analysis block, performing spatial adjacency consistency verification, and if adjacent geographic analysis blocks are in the same influence area and the KL divergence of the corresponding disaster-causing state variables is less than a probability target value, then the adjacent geographic analysis blocks are merged into a comprehensive risk block; the elevation deviation of a digital elevation model caused by terrain shielding, the spatial interpolation uncertainty caused by sparse rainfall sites, and the distortion of InSAR deformation caused by atmospheric and decorrelation interference are converted into quantifiable error measurement indexes, thereby improving the accuracy of risk assessment.

Description

A method and system for processing geological hazard risk assessment data Technical Field

[0001] This application relates to the field of geological disaster technology, and in particular to a geological disaster risk assessment data processing method and system. Background Technology

[0002] Geological disasters, such as landslides, collapses, debris flows, and ground subsidence, are characterized by their suddenness, complex causes, and high destructiveness, seriously threatening people's lives and property and the stable operation of critical infrastructure. Scientific and reliable risk assessment is a key support for achieving precise disaster prevention and mitigation, optimizing land spatial planning, and improving emergency response capabilities. In recent years, with the development of technologies such as remote sensing (RS), geographic information systems (GIS), global navigation satellite systems (GNSS), the Internet of Things, and big data analytics, geological disaster risk assessment is gradually shifting from traditional experience-based judgment to quantitative models driven by multi-source data. Current mainstream methods typically construct a disaster-causing factor system based on multi-source data such as topography, geology, rainfall, historical disaster sites, and InSAR deformation, and combine this with statistical models or machine learning algorithms to generate regional risk level maps, serving the identification of disaster hazards and risk zoning.

[0003] However, despite significantly enhanced data acquisition capabilities and increasingly complex model structures, existing risk assessment methods still suffer from significant technical deficiencies in data processing and result representation. They generally lack systematic quantification and propagation analysis of multi-level uncertainties during the assessment process, making it difficult to verify the reliability of risk assessment results and limiting their decision-making support capabilities. Specifically, the input data itself contains observable errors. For example, digital elevation models (DEMs) suffer from elevation deviations due to sensor accuracy and terrain occlusion; rainfall data lacks spatial representativeness due to sparse stations or remote sensing inversion errors; and InSAR surface deformation rates are easily distorted by atmospheric delay and decoherence interference. If these errors are directly used as deterministic inputs in modeling without quantification, they will be fully propagated to the final risk results, causing systematic biases. This not only masks the spatial differentiation characteristics of the true risk but may also lead to the underestimation of high-risk areas or the misjudgment of low-risk areas, significantly reducing the accuracy of risk assessments and consequently affecting the scientific nature of disaster prevention resource allocation. Summary of the Invention

[0004] This application aims to at least partially address one of the technical problems in the related art.

[0005] To achieve the above objectives, this application proposes a method for processing geological hazard risk assessment data, including the following steps:

[0006] Step 1: Obtain the raw observation data of the assessment area: digital elevation model data, ground rain gauge data, and surface topography monitoring data;

[0007] Step 2: Generate multiple error metrics based on the original observation data: elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index;

[0008] Step 3: Map the original observation data and error measurement indicators to a regular grid to obtain a grid field, construct multiple influence zones, and divide the grid field into several geographic analysis blocks;

[0009] Step 4: Construct disaster-causing state variables based on the geographic analysis blocks, and perform spatial adjacency consistency test based on the disaster-causing state variables of the geographic analysis blocks. If adjacent geographic analysis blocks are in the same influence area and the KL divergence of the corresponding disaster-causing state variables is less than the probability target value, then the adjacent geographic analysis blocks are merged into a comprehensive risk block.

[0010] Step 5: Generate risk level classification results based on the comprehensive risk blocks.

[0011] Furthermore, generating the elevation error probability distribution includes the following steps:

[0012] Step 2a1: The digital elevation model data is parsed into several grids, each grid representing a cell, and each cell includes an elevation value and the corresponding geographic coordinates.

[0013] Step 2a2: Obtain the nominal vertical accuracy value and mark it as the basic error standard deviation;

[0014] Step 2a3: Traverse all pixels, construct a neighborhood window based on the current pixel, determine the pixel with the largest elevation value within the neighborhood window and mark it as the target pixel, determine the first vector based on the geographic coordinates of the target pixel and the current pixel, and determine the angle between the first vector and the horizontal plane and mark it as the line-of-sight occlusion angle.

[0015] Step 2a4: Determine the local slope based on the neighborhood window, and determine the pixel category based on the local slope and the line-of-sight occlusion angle: shadow occlusion pixel, steep slope distortion pixel, and normal terrain pixel;

[0016] Step 2a5: Determine the magnification factor based on the pixel category, and generate the error correction standard deviation based on the magnification factor and the basic error standard deviation;

[0017] Step 2a6: Using the elevation value of the pixel as the distribution center, construct a normal distribution of elevation error based on the error correction standard deviation.

[0018] Furthermore, if the line-of-sight occlusion angle is greater than the first angle, the pixel is determined to be a shadow occlusion pixel; if the line-of-sight occlusion angle is less than or equal to the first angle and the local slope is greater than the second angle, the pixel is determined to be a steep slope distortion pixel; if the line-of-sight occlusion angle is less than or equal to the first angle and the local slope is less than or equal to the second angle, the pixel is determined to be a normal terrain pixel.

[0019] Furthermore, generating the error covariance matrix includes the following steps:

[0020] Step 2b1: Based on the data from the ground rain gauges, obtain the elevation values ​​corresponding to the ground rain gauges in the assessment area and the effective daily rainfall sequence during the target assessment period;

[0021] Step 2b2: Obtain the variance of the rainfall samples; determine the spatial correlation length based on the effective daily rainfall sequence;

[0022] Step 2b3: Based on the pre-constructed meteorological dataset, extract the average wind direction vector of the assessment area during the target assessment period, and determine the wind direction azimuth based on the average wind direction vector;

[0023] Step 2b4: Construct station pairs based on all surface rain gauges in the assessment area, determine the connection distance between each station pair, determine the elevation difference based on the corresponding elevation values ​​of the station pairs, and determine the azimuth of the station pairs;

[0024] Step 2b5: Construct the angle deviation based on the wind direction and azimuth angle;

[0025] Step 2b6: If the elevation difference is less than the first elevation and the angle deviation meets the target angle, define the first element value based on the sample variance, connection distance, and spatial correlation length; otherwise, define the second element value based on the sample variance, connection distance, and spatial correlation length.

[0026] Step 2b7: Construct a symmetric matrix. Define the diagonal elements of the symmetric matrix as the sample variances and the off-diagonal elements as the first or second element values. All the sample variances, first element values, and second element values ​​constitute the error covariance matrix.

[0027] Furthermore, the deformation rate uncertainty index is generated, including the following steps:

[0028] Step 2c1: Obtain SAR images and corresponding imaging timestamps based on surface topography monitoring data;

[0029] Step 2c2: Determine the interferometric pair based on the SAR images. The interferometric pair consists of two SAR images, corresponding imaging stamps, complex interferograms, and coherence coefficient diagrams.

[0030] Step 2c3: Perform phase unwrapping processing based on the complex interferogram to obtain a phase difference map;

[0031] Step 2c4: Based on the phase difference map, determine the surface deformation rate and coherence coefficient sequence of each pixel during the observation period, and determine the evaluation coherence based on the coherence coefficient sequence;

[0032] Step 2c5: Determine the delayed phase map based on the SAR image;

[0033] Step 2c6: Based on the phase difference map, the delayed phase map, and the surface deformation rate, generate the atmospheric corrected phase residual at the pixel for each interferometric pair; determine the standard deviation of the phase residual based on the atmospheric corrected phase residual.

[0034] Step 2c7: Determine the stable reference region based on the evaluation region and generate the deformation rate uncertainty index.

[0035] Furthermore, based on the assessment region, a stable reference region is determined and a deformation rate uncertainty index is generated, including the following steps:

[0036] Step 2c71: Select a sub-region within the assessment area with a local slope of less than 5°, land use type of exposed bedrock or desert, and no historical geological disaster record as a stability reference area;

[0037] Step 2c72: Based on the stable reference area, extract the surface deformation rate of all the pixels, and generate the standard deviation of the surface deformation rate according to the surface deformation rate.

[0038] Step 2c73: Based on the evaluation of coherence, the standard deviation of phase residuals, and the standard deviation of surface deformation rate, the deformation rate uncertainty index is obtained.

[0039] Furthermore, a spatial adjacency consistency test is performed on the disaster-causing state variables of the geographic analysis blocks. If the KL divergence of the disaster-causing state variables of adjacent geographic analysis blocks is less than the probability target value and they are located in the same influence zone, then the adjacent geographic analysis blocks are merged into a comprehensive risk block, including the following steps:

[0040] Step 41: Extract adjacent block pairs based on geographic analysis blocks;

[0041] Step 42: Traverse all adjacent geographic analysis block pairs, obtain the first factor based on the disaster-causing state variables of the adjacent block pairs, and generate KL divergence and the second factor based on the first factor; if both KL divergence and the second factor are less than the probability target value, it means that the disaster-causing degree of the adjacent block pairs is consistent, and proceed to the next step; otherwise, end.

[0042] Step 43: If adjacent block pairs are all in the same influence zone, then the consistency condition is met and the adjacent block pairs are merged into a comprehensive risk block.

[0043] Step 44: Based on the comprehensive risk block, re-execute step 4 to obtain the disaster state update variables;

[0044] Step 45: Repeat steps 41 to 44 until no adjacent geospatial block pairs are merged in any round of traversal.

[0045] Furthermore, the disaster-causing state variables include continuous variables and discrete variables.

[0046] Furthermore, the first factor is a continuous variable, and the second factor is the deviation value of the continuous variable.

[0047] This application also provides a geological hazard risk assessment data processing system, including the following modules:

[0048] Data acquisition module: used to acquire raw observation data of the assessment area: digital elevation model data, ground rain gauge data, and surface topography monitoring data;

[0049] Data processing module: used to generate multiple error metrics based on the raw observation data: elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index;

[0050] The partitioning module is used to map the original observation data and error metrics to a regular grid to obtain a grid field, construct multiple influence zones, and divide the grid field into several geographic analysis blocks.

[0051] Merging module: used to construct disaster-causing state variables based on the geographic analysis blocks, and perform spatial adjacency consistency checks based on the disaster-causing state variables of the geographic analysis blocks. If adjacent geographic analysis blocks are in the same influence area and the KL divergence of the corresponding disaster-causing state variables is less than the probability target value, then the adjacent geographic analysis blocks are merged into a comprehensive risk block.

[0052] Risk classification module: Generates risk classification results based on the comprehensive risk blocks.

[0053] Compared with existing technologies, the geological disaster risk assessment data processing method and system provided in this application abandons the existing approach of relying on statistical models or machine learning algorithms for risk assessment. Based on digital elevation model data, ground rain gauge data, and surface topography monitoring data, it generates elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index, respectively. This transforms the elevation deviation caused by topographic shading in the digital elevation model, the spatial interpolation uncertainty caused by sparse rainfall stations, and the distortion caused by atmospheric and decoherent interference in InSAR deformation into quantifiable error metrics, thereby improving the accuracy of risk assessment and ensuring the scientific nature of disaster prevention resource allocation.

[0054] This application maps the original observation data and error measurement indicators together to a regular grid, and combines the influence area to divide the geographic analysis blocks. During the regional growth process, it introduces the similarity judgment of the standard deviation of elevation error correction and the uncertainty index of deformation rate to ensure that the same geographic analysis block is not only located in the same tectonic control area, but also has similar data quality levels.

[0055] This application transforms the original parameters into hazard-causing state variables that carry uncertainty information. This fundamentally avoids the systematic bias caused by directly using erroneous observations for risk assessment.

[0056] This application introduces a dual criterion of KL divergence and a second factor during the merging stage of adjacent block pairs. It ensures that adjacent block pairs, belonging to the same influence zone, are merged not only based on statistical similarity but also controlled by the consistency of the constructed control background. This effectively prevents false similarities caused by local data noise or interpolation bias from being misjudged as real risks, thereby suppressing the disorderly diffusion of uncertainty during the spatial aggregation process. Attached Figure Description

[0057] The above and / or additional aspects and advantages of this application will become apparent and readily understood from the following description of the embodiments taken in conjunction with the accompanying drawings, wherein:

[0058] Figure 1 is a flowchart of a geological disaster risk assessment data processing method provided in an embodiment of this application;

[0059] Figure 2 is a structural diagram of a geological disaster risk assessment data processing system provided in an embodiment of this application;

[0060] Figure 3 is a block diagram of an electronic device provided in an embodiment of this application. Detailed Implementation

[0061] The embodiments of this application are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain this application, and should not be construed as limiting this application.

[0062] The following description, with reference to the accompanying drawings, illustrates a geological disaster risk assessment data processing method and system according to an embodiment of this application.

[0063] As shown in Figure 1, a method for processing geological hazard risk assessment data includes the following steps:

[0064] Step 1: Obtain the raw observation data of the assessment area: digital elevation model data, ground rain gauge data, and surface topography monitoring data.

[0065] The digital elevation model (DEM) data is derived from ASTER GDEM V3 data with a resolution of 30 meters, covering the entire assessment area. Each cell contains an elevation value (unit: meters) and its corresponding WGS84 geographic coordinates.

[0066] The data from the ground rain gauge stations, provided by the local meteorological bureau, includes several days of rainfall records. Each record contains the station's latitude and longitude, altitude, and daily effective rainfall (unit: mm).

[0067] The surface topography monitoring data includes SAR imagery in wide interferometric (IW) mode with a spatial resolution of 5 m × 20 m, and also includes the corresponding imaging timestamps.

[0068] Step 2: Generate multiple error metrics based on the original observation data: elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index.

[0069] Generating the elevation error probability distribution includes the following steps:

[0070] Step 2a1: The digital elevation model data is parsed into several grids, each grid representing a pixel, and each pixel includes an elevation value and the corresponding geographic coordinates.

[0071] The digital elevation model data is parsed into a regular raster structure, with each raster corresponding to a cell, storing its elevation value and the geographic coordinates of its center point.

[0072] Step 2a2: Obtain the nominal vertical accuracy value and mark it as the standard deviation of the basic error.

[0073] Read the nominal vertical accuracy value from the metadata file accompanying the digital elevation model.

[0074] Step 2a3: Traverse all pixels, construct a neighborhood window based on the current pixel, determine the pixel with the largest elevation value within the neighborhood window and mark it as the target pixel, determine the first vector based on the geographic coordinates of the target pixel and the current pixel, and determine the angle between the first vector and the horizontal plane and mark it as the line-of-sight occlusion angle.

[0075] For each pixel in the digital elevation model (DEM) data, a 7×7 square neighborhood window is constructed centered on that pixel. This window covers a local terrain region extending 3 pixels horizontally and 3 pixels vertically from the current pixel. Within this window, elevation values ​​are compared pixel by pixel, and the pixel with the highest elevation is identified and designated as the target pixel. Geographic coordinates are latitude and longitude, using UTM (Universal Transverse Mercator) projection to convert the latitude and longitude of the current and target pixels to planar coordinates within the same projection zone (e.g., UTMZone 48N). The first vector is obtained by subtracting the planar coordinates of the current and target pixels. The projected length of the first vector on the horizontal plane is calculated using Euclidean distance. Based on the planar coordinates and projected length of the current and target pixels, the angle between the first vector and the horizontal plane is calculated using the arctangent function. This angle represents the viewing angle from the current pixel to the target pixel in the neighborhood, used to determine whether the current pixel is in a terrain shadow or a blocked viewing area. It is a key basis for distinguishing pixel error amplification.

[0076] Step 2a4: Determine the local slope based on the neighborhood window, and determine the pixel category based on the local slope and the line-of-sight occlusion angle: shadow occlusion pixel, steep slope distortion pixel, and normal terrain pixel.

[0077] Based on the neighborhood window where the current pixel is located, the surface tilt of the local area is estimated by the least squares plane fitting method to obtain the local slope: the optimal plane equation is fitted with the plane coordinates and elevation values ​​of the neighborhood window, and the local slope is obtained by minimizing the sum of squared residuals. The above-mentioned least squares plane fitting process is a conventional technique, and will not be described in detail in this embodiment.

[0078] If the line-of-sight occlusion angle is greater than the first angle, the pixel is identified as a shadow occlusion pixel; if the line-of-sight occlusion angle is less than or equal to the first angle and the local slope is greater than the second angle, the pixel is identified as a steep slope distortion pixel; if the line-of-sight occlusion angle is less than or equal to the first angle and the local slope is less than or equal to the second angle, the pixel is identified as a normal terrain pixel.

[0079] In this embodiment, each pixel of the digital elevation model data is classified based on its line-of-sight occlusion angle and local slope to distinguish terrain regions with different error sources. The first angle is set to 30°, and the second angle to 45°. If a pixel's line-of-sight occlusion angle is greater than 30°, it is identified as a shadowed pixel; if the line-of-sight occlusion angle is less than or equal to 30° and the local slope is greater than 45°, it is identified as a steep slope distortion pixel; if the line-of-sight occlusion angle is less than or equal to 30° and the local slope is less than or equal to 45°, it is identified as a normal terrain pixel. When the line-of-sight obstruction angle exceeds 30°, it indicates that there is significant terrain obstruction of the sensor signal round-trip path within the neighborhood window. This easily forms shadows or incoherent areas in optical or radar measurements, resulting in a systematic underestimation or absence of elevation. 30° is a significant threshold for the obstruction effect widely used in InSAR and stereo photogrammetry. Although there is no obstruction in areas with slopes exceeding 45°, the accuracy of digital elevation model inversion decreases significantly due to pixel projection distortion (such as radar overlay and optical stretching) and increased surface complexity. 45° is close to the natural angle of repose of most soil and rock materials and is also a recognized critical value for steep slopes in data quality assessment.

[0080] Step 2a5: Determine the magnification factor based on the pixel category, and generate the error correction standard deviation based on the magnification factor and the basic error standard deviation.

[0081] Based on the aforementioned pixel classification results, differentiated magnification factors are set to reflect their actual impact on elevation error. The magnification factor is set to 1.8 for shadowed pixels, 1.3 for steep slope distortion pixels, and 1.0 for regular terrain pixels. These magnification factors are used to modulate the baseline error standard deviation. Multiplying the baseline error standard deviation by the corresponding magnification factor generates an error correction standard deviation adapted to the terrain conditions.

[0082] Step 2a6: Using the elevation value of the pixel as the distribution center, construct a normal distribution of elevation error based on the error correction standard deviation.

[0083] The elevation value of this pixel is taken as the mean (i.e., the center of the distribution) of the probability distribution; the standard deviation of the error correction of this pixel is taken as the standard deviation of the normal distribution; thus, the elevation error of this pixel is defined to follow a one-dimensional normal distribution:

[0084] ; P H (h) is the probability density of the elevation value h, which, under a given normal distribution, represents the likelihood of the elevation value h occurring. This represents the mean of a normal distribution of elevation error. This represents the standard deviation of the normal distribution of elevation error. This represents the corrected standard deviation of the error.

[0085] This embodiment effectively addresses the systematic elevation bias in Digital Elevation Models (DEMs) caused by sensor accuracy limitations and terrain shading effects by constructing an elevation error probability distribution as a key error metric. Instead of treating DEM elevations as deterministic true values, this embodiment establishes a normal distribution of elevation errors using pixel elevations as the mean and the error correction standard deviation as the dispersion parameter. This explicitly characterizes the uncertainty level of elevation data under different terrain conditions, avoiding the risk misjudgment caused by directly applying distorted terrain features to risk assessment due to ignoring elevation errors in traditional methods.

[0086] Generating the error covariance matrix includes the following steps:

[0087] Step 2b1: Based on the data from the ground rain gauges, obtain the elevation values ​​corresponding to the ground rain gauges in the assessment area and the effective daily rainfall sequence during the target assessment period.

[0088] Location data refers to the latitude and longitude of surface rain gauges. The main flood season from June to September each year is selected as the target evaluation period. Stations with a missing rate of more than 20% or consecutive missing days are removed to obtain the final valid stations. Each station contains several valid daily rainfall values.

[0089] Step 2b2: Obtain the variance of the rainfall samples; determine the spatial correlation length based on the effective daily rainfall sequence.

[0090] The mean of effective daily rainfall is calculated based on the effective daily rainfall values, and then substituted into the variance formula to calculate the sample variance of rainfall. The preferred spatial correlation length is 48.6 km, as the spatial correlation length of rainfall in similar geomorphic units is generally concentrated in the range of 40–55 km. In this embodiment, L=48.6 km, obtained by fitting measured data, falls precisely within this typical range and is consistent with the average influence scale of major weather systems in the region. That is, when the station spacing exceeds approximately 50 km, the correlation of daily rainfall significantly decreases to a stable state. Therefore, using 48.6 km as the spatial correlation length not only conforms to the local meteorological and topographical coupling pattern but also ensures that the subsequent error covariance matrix accurately reflects the uncertain spatial structure of rainfall interpolation.

[0091] Step 2b3: Based on the pre-constructed meteorological dataset, extract the average wind direction vector of the assessment area during the target assessment period, and determine the wind direction azimuth based on the average wind direction vector.

[0092] The pre-constructed meteorological dataset has a spatial resolution of 0.25°×0.25° and a temporal resolution of hourly, including zonal and meridional wind speed components at a height of 10 meters. For each time step, the regional zonal and meridional average wind speed components are calculated; subsequently, the zonal and meridional average wind speed components for all time steps are time-averaged to obtain the regional average wind direction vector for the entire evaluation period.

[0093] The wind direction azimuth angle is calculated based on the average wind vector (defined as the angle rotated clockwise from due north to the direction of the wind's origin). The formula is as follows:

[0094] arctan is the arctangent function in the four quadrants, α is the wind direction azimuth, and π is pi. This represents the zonal average wind speed component. This represents the meridional average wind speed component.

[0095] Step 2b4: Construct station pairs based on all ground rain gauges in the assessment area, determine the connection distance between each station pair, determine the elevation difference based on the corresponding elevation values ​​of the station pairs, and determine the azimuth of the station pairs.

[0096] The Euclidean distance between each pair of stations is calculated as the connection distance. The azimuth angle is defined as the clockwise angle between the direction from one surface rain gauge station to another relative to true north. It can be calculated using the four-quadrant arctangent function, which will not be described in detail in this embodiment.

[0097] Step 2b5: Construct the angle deviation based on the wind direction and azimuth.

[0098] Step 2b6: If the elevation difference is less than the first elevation and the angle deviation meets the target angle, define the first element value based on the sample variance, connection distance and spatial correlation length; otherwise, define the second element value based on the sample variance, connection distance and spatial correlation length.

[0099] The first elevation is 200 meters, and the target angle is ±45°. If the elevation difference is less than the first elevation and the angle deviation is less than 45° or greater than 315°, this embodiment assumes that the stations are located in similar terrain and arranged along the prevailing wind direction, and their rainfall errors have a strong spatial correlation. Therefore, the first element value is constructed as follows:

[0100] , For sample variance, Let L be the connection distance between the i-th and j-th surface rain gauges, and L be the spatial correlation length. Otherwise, construct the second element value:

[0101] .

[0102] Only when the elevation difference between two stations is less than 200 meters and the angle between the line connecting them and the prevailing wind is within ±45° is the rainfall error considered to have a strong correlation, and the covariance is amplified by 1.2 times after exponential decay; in other cases, it is considered to be controlled by different microclimates or topographic effects, and the correlation is suppressed to 0.5 times.

[0103] Step 2b7: Construct a symmetric matrix. Define the diagonal elements of the symmetric matrix as the sample variances and the off-diagonal elements as the first or second element values. All the sample variances, first element values, and second element values ​​constitute the error covariance matrix.

[0104] Construct an N×N symmetric matrix, where N≥2. Let the diagonal elements of the symmetric matrix represent the sample variances, and define the remaining elements as either the first or second element value. All sample variances, first element values, and second element values ​​constitute the error covariance matrix.

[0105] Generating the deformation rate uncertainty index includes the following steps:

[0106] Step 2c1: Obtain SAR images and corresponding imaging timestamps based on surface topography monitoring data.

[0107] Step 2c2: Determine the interferometric pair based on the SAR images. The interferometric pair consists of two SAR images, corresponding imaging stamps, complex interferograms, and coherence coefficient diagrams.

[0108] From all single-look complex-format SAR images acquired within the assessment area coverage period, pairwise screening is performed based on preset time baseline thresholds (≤270 days) and vertical baseline thresholds (≤300 meters) to form several candidate interferometric pairs. Each interferometric pair consists of a primary image and a secondary image, and its corresponding precise imaging time (i.e., imaging stamp) is recorded. Precise registration is performed on each image pair: coarse registration is performed using orbital parameters and a digital elevation model (DEM), followed by sub-pixel-level fine registration in the azimuth and range directions using frequency domain cross-correlation or phase correlation methods, allowing the secondary image to be resampled to the geometric grid of the primary image. After registration, for each pixel location, the complex conjugate of the complex values ​​of the primary and secondary images is multiplied to generate a complex interferogram. Within a local window (usually 5×5 or 9×9 pixels), the coherence coefficient is calculated based on the amplitude data of the complex interferogram and the primary and secondary images to obtain a correlation coefficient map.

[0109] Step 2c3: Perform phase unwrapping processing based on the complex interferogram to obtain a phase difference map;

[0110] According to the selected threshold method (the threshold is preferably 0.5), each pixel on the coherence coefficient map is traversed. If the coherence coefficient of a pixel is greater than or equal to the set threshold, a high value (such as 1) is assigned to the pixel on the quality map, indicating that it is a high-quality, high-reliability pixel; conversely, if the coherence coefficient is lower than the threshold, a low value (such as 0) is assigned, indicating that the pixel may be noise or the data is unreliable due to other reasons, thus obtaining the quality map.

[0111] The coherence coefficient is normalized and directly assigned as a weight value to the corresponding pixel to characterize the reliability of the phase data at that location. Based on this, a weighted least squares phase unwrapping method is employed: a large sparse linear equation system covering the entire interferometric region is constructed, with the unknown being the true (unwrapped) phase value of each effective pixel. The core constraint of this equation system is that the phase gradient between adjacent pixels should be as close as possible to the observation gradient derived from the wrapped phase. The quality map is embedded in the equation system as a weight matrix, ensuring that the phase relationship in high-coherence regions dominates the solution process, while the influence of low-coherence regions is suppressed. This linear system is efficiently solved using iterative algorithms (such as the preconditional conjugate gradient method), ultimately outputting a continuous phase field without 2π jumps. This phase field is the unwrapped phase difference distribution; the phase difference values ​​of all pixels together constitute the phase difference map, accurately reflecting the cumulative phase difference changes caused by terrain undulations or surface deformation during the imaging of two SAR images. This can be further combined with radar wavelength, incident angle, and baseline parameters to convert into elevation error or deformation, providing crucial input for geological hazard risk assessment. Untangling is a conventional technique, and will not be described in detail in this embodiment.

[0112] Step 2c4: Based on the phase difference map, determine the surface deformation rate and coherence coefficient sequence of each pixel during the observation period, and determine the evaluation coherence based on the coherence coefficient sequence.

[0113] The operating wavelength M was acquired using a synthetic aperture radar (SAR) system, and the cumulative deformation was calculated as (M / 4Π)*X using the phase difference value X from the phase difference map. The actual acquisition time interval between the main and auxiliary images was obtained, and the cumulative deformation was divided by the time interval to obtain the surface deformation rate, which reflects the relative surface deformation of each pixel's coverage area during the observation period.

[0114] Based on the coherence coefficient map, the corresponding coherence coefficient value is extracted for each pixel from all the interference pairs it participates in, forming a coherence coefficient time series. In this embodiment, the average coherence coefficient of the coherence coefficient time series is used as the evaluation coherence.

[0115] Step 2c5: Determine the delayed phase map based on the SAR image.

[0116] Three-dimensional atmospheric profile data corresponding to the SAR image imaging time are extracted from the reanalysis meteorological dataset, including parameters such as air pressure, temperature, and water vapor content. This meteorological data is used to drive an empirical model to calculate the dry and wet delay components at each SAR pixel location, and these are converted into equivalent radar line-of-sight (LOS) phase delays. The dry and wet delay phases are then superimposed to generate a delay phase map covering the entire interferometric region. This is a conventional technique, and will not be described in detail in this embodiment.

[0117] Step 2c6: Based on the phase difference map, the delayed phase map, and the surface deformation rate, generate the atmospheric corrected phase residual at the pixel for each interferometric pair; determine the standard deviation of the phase residual based on the atmospheric corrected phase residual.

[0118] The atmospheric correction phase residual and the standard deviation of the phase residual are expressed as follows:

[0119] ;

[0120] ;

[0121] Where, r k (x, y) represents the atmospheric-corrected phase residual of the k-th interferometer pair at pixel (x, y). Let be the phase difference of the k-th interference pair after unwrapping at pixel (x,y). Let be the atmospheric delayed phase map value of pixel (x,y), and v(x,y) be the surface deformation rate at pixel (x,y). The actual acquisition time interval between the primary and secondary images of the k-th interferometric pair is given. Let K be the standard deviation of the phase residual at pixel (x,y) for the k-th interferometer pair, where K is the total number of interferometer pairs. Let be the mean atmospheric correction phase residual of pixel (x,y) across all interference pairs.

[0122] Since the phase values ​​themselves are periodic ([-π,π]), if there are obvious phase jumps in the residual sequence, normalization should be performed before calculation (such as mapping all residuals to the interval [0,2π)) to avoid underestimation of the standard deviation due to periodicity.

[0123] The atmospheric correction phase residual is obtained by subtracting the atmospheric delay phase from the unwrapped phase difference of each pixel in the interferometric pair, and further subtracting the phase contribution corresponding to the theoretical deformation calculated based on the long-term deformation rate. This effectively removes two main types of systematic errors: one is the line-of-sight path delay caused by water vapor and dry air, and the other is the linear cumulative deformation caused by stable subsidence or uplift. After correction, the residual phase mainly contains information such as nonlinear deformation, land cover changes, or incompletely modeled minor atmospheric disturbances. This not only significantly improves the purity of the deformation signal but also provides a more sensitive criterion for identifying potential disaster precursors.

[0124] Building upon this, the standard deviation of the atmospheric-corrected phase residual for each pixel across all interferometric pairs is further calculated. This quantifies the dispersion of phase fluctuations at that location, i.e., the level of uncertainty. Regions with smaller standard deviations indicate highly consistent phase behavior, minimal influence from random noise, and reliable deformation rate results. Conversely, regions with larger standard deviations may be in an active deformation state or severely affected by decoherence, residual atmospheric disturbances, etc. In summary, introducing the atmospheric-corrected phase residual and its standard deviation not only achieves refined separation and quantification of InSAR observation errors but also enhances the spatial reliability of the deformation field. This processing strategy, which balances physical mechanisms and statistical robustness, significantly improves the applicability and early warning effectiveness of surface deformation monitoring in complex mountainous or urban environments, laying a solid data foundation for the fusion analysis of multi-source disaster-causing factors.

[0125] Step 2c7: Determine the stable reference region based on the evaluation region and generate the deformation rate uncertainty index:

[0126] Step 2c71: Select sub-regions within the assessment area with local slopes less than 5°, land use types of exposed bedrock or desert, and no historical geological disaster records as stability reference areas.

[0127] Based on the assessment area, flat or gentle slope areas with a slope of less than 5° were selected; combined with the global land use / cover dataset, the land use type was further limited to exposed bedrock or desert to exclude non-tectonic deformation disturbances caused by vegetation cover, cultivation disturbance, or urban activities; the historical geological disaster database was overlaid to eliminate any candidate areas within the buffer zone (set to a radius of 500 meters) of historical disaster points; within the intersection of the above three conditions, a spatially continuous sub-region with an area of ​​not less than 1 square kilometer and good coherence was selected as the stable reference area.

[0128] Step 2c72: Based on the stable reference area, extract the surface deformation rate of all the pixels, and generate the standard deviation of the surface deformation rate according to the surface deformation rate.

[0129] The surface deformation rate value corresponding to each pixel in the stable reference area is extracted. Statistical analysis is performed on the sample set composed of all extracted rate values ​​to calculate its arithmetic mean and standard deviation.

[0130] Step 2c73: Based on the evaluation of coherence, the standard deviation of phase residuals, and the standard deviation of surface deformation rate, the deformation rate uncertainty index is obtained.

[0131] U(x,y) is the uncertainty index of the deformation rate of pixel (x,y). For evaluating the coherence of pixel (x, y), Let (x, y) be the standard deviation of the phase residuals of pixel (x, y). This represents the maximum standard deviation of the phase residual. The standard deviation of the surface deformation rate. This represents the maximum standard deviation of the Earth's surface deformation rate.

[0132] Step 3: Map the original observation data and error metrics to a regular grid to obtain a grid field, construct multiple influence zones, and divide the grid field into several geographic analysis blocks.

[0133] A regular grid structure covering the evaluation area is constructed. This structure is formed by dividing the area at equal intervals of longitude and latitude, with each grid cell having a side length of no more than 50 meters. The original observation data and error metrics are resampled to the regular grid structure to obtain the following grid field with the same resolution.

[0134] The affected areas include fault-affected areas, landform abrupt change lines, and areas of artificial disturbance.

[0135] In the assessment area, line features of active fault zones are extracted based on a pre-constructed regional geological map, and a 500-meter bidirectional buffer operation is performed to generate a surface layer representing the fault-affected area. Topographic analysis is performed using a digital elevation model covering the assessment area to obtain the local slope and curvature of each raster. Raster locations satisfying the condition of an absolute curvature greater than 0.02 and a local slope variation exceeding 15 degrees within an adjacent 100-meter range are extracted as geomorphic abrupt change lines, and a binary geomorphic category layer is constructed based on this (1 represents a geomorphic abrupt change area, 0 represents a non-geomorphic abrupt change area). A 100-meter buffer is applied to the railway centerline, the reservoir's normal water level inundation boundary, and the open-pit mine boundary to create an artificial disturbance zone layer, used to identify areas of surface instability that may be caused by human engineering activities. The fault-affected area, geomorphic abrupt change lines, and artificial disturbance zone are overlaid to form an initial segmentation guide line. When applying the region growing algorithm, starting from any unassigned raster, it is checked whether it shares the same set of attributes with adjacent rasters. Specifically, these attributes include:

[0136] Fault influence zone attribution: Determine whether the grid is located within the fault influence zone;

[0137] Landform category: Based on the results of landform abrupt change line analysis, determine the landform type to which the raster belongs;

[0138] Artificial disturbance status: Confirm whether the grid is within an artificial disturbance zone, such as a railway buffer zone, a reservoir flooding boundary buffer zone, or an open-pit mine mining boundary buffer zone.

[0139] If two adjacent rasters are consistent in all three aspects, they are considered to belong to the same geographic feature block and are merged into the same block. This process is repeated until all rasters have been assigned to their respective blocks, ultimately forming a set of non-overlapping, spatially contiguous geographic analysis blocks.

[0140] Step 4: Construct disaster-causing state variables based on the geographic analysis blocks, and perform spatial adjacency consistency test based on the disaster-causing state variables of the geographic analysis blocks. If adjacent geographic analysis blocks are in the same influence area and the KL divergence of the corresponding disaster-causing state variables is less than the probability target value, then the adjacent geographic analysis blocks are merged into a comprehensive risk block.

[0141] This embodiment performs statistical analysis on the pixels of each geographic analysis block using multivariate parameters (raw observations and error metrics) to form a multidimensional probability distribution structure that comprehensively characterizes its potential disaster risk, i.e., the disaster-causing state variable. The multivariate parameters specifically include:

[0142] Local slope: The local slope of each pixel is obtained based on the local slope. The mean and standard deviation of the first sample are calculated based on the local slope. The variance of the first sample is obtained based on the standard deviation of the first sample. A truncated normal distribution is constructed based on the mean and variance of the first sample, with a lower limit of 0° and an upper limit of 90°, to reflect the overall distribution characteristics of the terrain slope within the block.

[0143] Historical rainfall: Obtain the historical annual rainfall for each geographic analysis block during the target evaluation period, calculate the second sample mean and second sample standard deviation based on the historical annual rainfall, calculate the second sample variance based on the second sample standard deviation, and construct a normal distribution based on the second sample mean and second sample variance.

[0144] Surface deformation rate: The deformation rate uncertainty index of each geographic analysis block is statistically analyzed and the mean value of the deformation rate uncertainty index is obtained. If the mean value of the deformation rate uncertainty index is greater than 0.3, it indicates that the deformation rate has low reliability. In this case, the upper 90th percentile of the deformation rate in the current geographic analysis block is used as the representative value, and the deformation rate not exceeding the representative value is taken as a landslide susceptibility condition. Otherwise, the deformation rate standard deviation is calculated based on the deformation rate, the deformation rate variance is constructed based on the deformation rate standard deviation, and a normal distribution of the deformation rate is constructed based on the deformation rate mean and deformation rate variance.

[0145] Rock type category: Statistical analysis of the percentage of pixels corresponding to each type of rock within the block, and construction of discrete probability distribution based on the percentage of pixels and the total number of pixels.

[0146] Landslide susceptibility status indicators: Perform the following checks for each geographic analysis block:

[0147] The percentage of pixels within a statistical geographic analysis block that simultaneously meet both of the following conditions:

[0148] Local slope >30°;

[0149] The lithology belongs to the low-permeability group (such as mudstone, shale, and other rocks that are easily softened and have low shear strength).

[0150] If the percentage exceeds the preset threshold (30%), the area is determined to have the geological and topographical combination conditions that make it prone to landslides, and a landslide-prone sign is set; otherwise, a non-landslide-prone sign is set.

[0151] The probability distributions of the aforementioned local slope, historical rainfall, deformation rate, lithology, and landslide susceptibility indicators are combined with state variables to form a multidimensional joint probability representation structure, which serves as the disaster-causing state variable for this geographic analysis block.

[0152] Step 41: Extract adjacent block pairs based on geographic analysis blocks.

[0153] Each geospatial block is treated as an independent spatial polygon object, and the topology analysis function in the geographic information system is used to identify blocks that share a boundary (i.e., the boundary length is greater than zero, excluding cases of only point contact). To avoid duplication, each pair of adjacent blocks is sorted by its unique identifier (such as block ID) to ensure that each pair of adjacent relationships is recorded only once.

[0154] Step 42: Traverse all adjacent geographic analysis block pairs, obtain the first factor based on the disaster-causing state variables of the adjacent block pairs, and generate KL divergence and the second factor based on the first factor; if both KL divergence and the second factor are less than the probability target value, it means that the disaster-causing degree of the adjacent block pairs is consistent, and proceed to the next step; otherwise, end.

[0155] The disaster-causing state variables include continuous variables (local slope, historical rainfall, deformation rate) and discrete variables (lithological type and landslide susceptibility indicators). The first factor is a continuous variable, which is the upper 90th percentile of all continuous variables. The second factor is the deviation value of the continuous variable (after normalization).

[0156] Three continuous hazard-causing variables—local slope, historical annual precipitation, and surface deformation rate—were extracted from all pixels within adjacent geographic analysis blocks. For each variable, an empirical probability distribution was constructed within the adjacent geographic analysis block pair. This was achieved by dividing the variable values ​​into several equally wide intervals, statistically analyzing the frequency of pixel occurrences within each interval, and normalizing these intervals to a probability quality function, thus forming two comparable probability distributions. Based on these two probability distributions, the KL divergence between them was calculated. KL divergence measures the information difference in the distribution of hazard-causing variables of one block relative to another: if the distributions of slope, precipitation, or deformation rate are highly similar between adjacent geographic analysis blocks, the KL divergence value is small; if the distribution patterns differ significantly, the KL divergence value is large. To avoid directional bias, a symmetrical form of KL divergence was used, i.e., the KL divergence was calculated once for each geographic analysis block and then averaged. Finally, the symmetrical KL divergences of the three continuous variables were arithmetically averaged to obtain the comprehensive KL divergence value for the adjacent block pair.

[0157] The target probabilities are set as follows: the KL divergence target value is 0.15, and the second factor target value is 0.2.

[0158] If the KL divergence of a pair of adjacent geographic analysis blocks is <0.15 and the second factor is <0.2, then the pair of geographic analysis blocks is considered to have a high degree of consistency in disaster-causing status, that is, their topographic, hydrological and deformation characteristics are similar in statistical distribution and high-risk tail. The pair of geographic analysis blocks will be retained and proceed to the next step.

[0159] If any indicator fails to meet the above conditions, it is determined that the disaster-causing characteristics of the two geographic analysis blocks are significantly different and inconsistent. Therefore, the subsequent processing of the geographic analysis block pair will be terminated, and it will no longer participate in the next step of analysis.

[0160] Step 43: If adjacent block pairs are all in the same influence zone, then the consistency condition is met and the adjacent block pairs are merged into a comprehensive risk block.

[0161] For any pair of adjacent geographic analysis blocks filtered in step 42, check their affiliation status in the three types of influence zones in turn:

[0162] Fault influence zone attribution: Determine whether the two blocks in the adjacent geographic analysis block pair are both completely located within the fault influence zone, or both completely located outside the fault influence zone.

[0163] Landform abrupt change status: Determine whether the two belong to the same landform abrupt change area or the same non-landform abrupt change area;

[0164] Artificial disturbance status: Determine whether both fall into at least one artificial disturbance zone, or whether neither falls into any artificial disturbance zone.

[0165] Only when the above three attribution statuses are completely consistent is the adjacent geographic analysis block pair considered to meet the construction consistency condition.

[0166] Once the construction consistency condition is met, the adjacent geographic analysis blocks are merged into a single comprehensive risk block, and the newly generated block inherits the geometric boundaries of the original two blocks.

[0167] Step 44: Based on the comprehensive risk block, re-execute step 4 to obtain the disaster state update variables;

[0168] Step 45: Repeat steps 41 to 44 until no adjacent geospatial block pairs are merged in any round of traversal.

[0169] Step 5: Generate risk level classification results based on the comprehensive risk blocks.

[0170] If a comprehensive risk area is marked as prone to landslides and meets any of the following conditions, it is classified as extremely high risk:

[0171] The upper 90th percentile of the surface deformation rate is greater than 5.0 mm per year;

[0172] The 90th percentile of the terrain slope is greater than 45 degrees.

[0173] The upper 90th percentile of the historical maximum 24-hour rainfall is greater than 200 mm.

[0174] If a comprehensive risk block is set as a landslide-prone state marker, but does not reach the extremely high risk condition, and the mean value of the deformation rate uncertainty index is less than or equal to 0.3, it is judged as high risk.

[0175] If a comprehensive risk block is set as a landslide-prone state marker and the average deformation rate uncertainty index is greater than 0.3, it is judged as medium risk.

[0176] If a comprehensive risk area is set to a non-landslide-prone status, it is considered low risk.

[0177] Each comprehensive risk zone is spatially represented according to its corresponding risk level, forming a comprehensive geological hazard risk zoning map covering the entire assessment area. In this process, each comprehensive risk zone is assigned a corresponding risk level and visually coded using different colors or symbols. For example: low-risk areas are represented by green; medium-risk areas by yellow; high-risk areas by orange; and extremely high-risk areas by red.

[0178] This application determines risk levels using multi-dimensional parameters such as landslide susceptibility indicators, the first factor, and the mean of the deformation rate uncertainty index. Compared to existing methods that classify risks solely based on average deformation rate or simple slope thresholds, this significantly improves the accuracy and precision of the assessment results. This embodiment uses the upper 90th percentile of each disaster-causing variable within a block as the risk discrimination criterion, highlighting the concentration trend of potentially high-risk pixels in local areas. Simultaneously, it quantifies data reliability using the deformation rate uncertainty index, distinguishing between areas with high deformation but low reliability and truly high-risk areas in risk assessment, effectively reducing false alarms. Furthermore, the risk level classification is based on comprehensive risk blocks formed after dual screening for consistency in disaster severity and impact area attribution, ensuring consistency in spatial attributes such as fault influence, geomorphic abrupt changes, and human disturbance within each risk unit, thereby avoiding the conflation of areas with different causal mechanisms. The final risk zoning results more closely reflect actual disaster development patterns, providing more targeted and operable support for geological disaster monitoring and early warning, key hazard identification, and optimal allocation of disaster prevention resources.

[0179] As shown in Figure 2, this embodiment also discloses a geological hazard risk assessment data processing system, including the following modules:

[0180] Data acquisition module: used to acquire raw observation data of the assessment area: digital elevation model data, ground rain gauge data, and surface topography monitoring data;

[0181] Data processing module: used to generate multiple error metrics based on the raw observation data: elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index;

[0182] The partitioning module is used to map the original observation data and error metrics to a regular grid to obtain a grid field, construct multiple influence zones, and divide the grid field into several geographic analysis blocks.

[0183] Merging module: used to construct disaster-causing state variables based on the geographic analysis blocks, and perform spatial adjacency consistency checks based on the disaster-causing state variables of the geographic analysis blocks. If adjacent geographic analysis blocks are in the same influence area and the KL divergence of the corresponding disaster-causing state variables is less than the probability target value, then the adjacent geographic analysis blocks are merged into a comprehensive risk block.

[0184] Risk classification module: Generates risk classification results based on the comprehensive risk blocks.

[0185] To implement the above embodiments, this application also proposes an electronic device. Please refer to FIG2 and FIG3, which are schematic diagrams of the structure of the electronic device provided in the embodiments of this application. As shown in FIG3, the electronic device 500 includes: a processor 501 and a memory 502 communicatively connected to the processor 501; the memory 502 stores computer execution instructions; the processor 501 executes the computer execution instructions stored in the memory to implement the method provided in the foregoing embodiments.

[0186] To implement the above embodiments, this application also proposes a computer-readable storage medium storing computer-executable instructions, which, when executed by a processor, are used to implement the methods provided in the foregoing embodiments.

[0187] To implement the above embodiments, this application also proposes a computer program product, including a computer program that, when executed by a processor, implements the methods provided in the foregoing embodiments.

[0188] The storage medium mentioned above can be a read-only memory, a disk, or an optical disk, etc. Although embodiments of this application have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting this application. Those skilled in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of this application.

Claims

1. A method for processing geological hazard risk assessment data, characterized in that, Includes the following steps: Step 1: Obtain the original observation data of the assessment area: digital elevation model data, ground rain gauge data, and surface topography monitoring data; Step 2: Generate multiple error metrics based on the original observation data: elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index; The generation of the deformation rate uncertainty index includes the following steps: Step 2c1, acquiring SAR images and corresponding imaging timestamps based on surface topography monitoring data; Step 2c2, determining interferometric pairs based on the SAR images, wherein the interferometric pair consists of two SAR images, corresponding imaging timestamps, complex interferograms, and coherence coefficient maps; Step 2c3, performing phase unwrapping processing based on the complex interferograms to obtain a phase difference map; Step 2c4, determining the surface deformation rate and coherence coefficient sequence for each pixel during the observation period based on the phase difference map, and determining the coherence assessment based on the coherence coefficient sequence; Step 2c5, determining a delayed phase map based on the SAR images; Step 2c6, generating the interferometric pair for each pixel based on the phase difference map, delayed phase map, and surface deformation rate. Step 2c7: Determine the atmospheric corrected phase residual; determine the standard deviation of the phase residual based on the atmospheric corrected phase residual; Step 2c7: Determine the stable reference area based on the assessment area and generate the deformation rate uncertainty index; Step 3: Map the original observation data and error metric to a regular grid to obtain a grid field, construct multiple influence areas and divide the grid field into several geographic analysis blocks; Step 4: Construct disaster-causing state variables based on the geographic analysis blocks, perform spatial adjacency consistency test based on the disaster-causing state variables of the geographic analysis blocks, if adjacent geographic analysis blocks are in the same influence area and the KL divergence of the corresponding disaster-causing state variables is less than the probability target value, then merge the adjacent geographic analysis blocks into a comprehensive risk block; Step 5: Generate risk level classification results based on the comprehensive risk blocks.

2. The geological hazard risk assessment data processing method according to claim 1, characterized in that, The process of generating an elevation error probability distribution includes the following steps: Step 2a1, parsing the digital elevation model data into several grids, each grid representing a pixel, and each pixel including an elevation value and corresponding geographic coordinates; Step 2a2, obtaining the nominal vertical accuracy value and marking it as the basic error standard deviation; Step 2a3, traversing all pixels, constructing a neighborhood window based on the current pixel, determining the pixel with the largest elevation value within the neighborhood window and marking it as the target pixel, determining a first vector based on the geographic coordinates of the target pixel and the current pixel, and determining the angle between the first vector and the horizontal plane and marking it as the line-of-sight occlusion angle; Step 2a4, determining the local slope based on the neighborhood window, and determining the pixel category based on the local slope and the line-of-sight occlusion angle: shadow occlusion pixel, steep slope distortion pixel, and normal terrain pixel; Step 2a5, determining the magnification factor based on the pixel category, and generating an error correction standard deviation based on the magnification factor and the basic error standard deviation; Step 2a6, using the elevation value of the pixel as the distribution center, constructing a normal distribution of elevation error based on the error correction standard deviation.

3. The geological hazard risk assessment data processing method according to claim 2, characterized in that, If the line-of-sight occlusion angle is greater than the first angle, the pixel is identified as a shadow occlusion pixel; if the line-of-sight occlusion angle is less than or equal to the first angle and the local slope is greater than the second angle, the pixel is identified as a steep slope distortion pixel; if the line-of-sight occlusion angle is less than or equal to the first angle and the local slope is less than or equal to the second angle, the pixel is identified as a normal terrain pixel.

4. The geological hazard risk assessment data processing method according to claim 1, characterized in that, Generating the error covariance matrix includes the following steps: Step 2b1, obtaining the elevation values ​​corresponding to the ground rain gauges in the assessment area and the effective daily rainfall sequence during the target assessment period based on the ground rain gauge data; Step 2b2, obtaining the rainfall sample variance; determining the spatial correlation length based on the effective daily rainfall sequence; Step 2b3, extracting the average wind direction vector of the assessment area during the target assessment period based on the pre-constructed meteorological dataset, and determining the wind direction azimuth based on the average wind direction vector; Step 2b4, constructing station pairs based on all ground rain gauges in the assessment area, determining the connection distance of each station pair, and determining the connection distance of each station pair based on the corresponding... Determine the elevation difference based on the elevation value; determine the azimuth of the station pair; step 2b5, construct the angle deviation based on the wind direction and azimuth; step 2b6, if the elevation difference is less than the first elevation and the angle deviation meets the target angle, define the first element value based on the sample variance, connection distance, and spatial correlation length; otherwise, define the second element value based on the sample variance, connection distance, and spatial correlation length; step 2b7, construct a symmetric matrix, define the diagonal elements of the symmetric matrix as the sample variance, and define the off-diagonal elements as the first or second element value, and all the sample variances, first element values, and second element values ​​constitute the error covariance matrix.

5. The geological hazard risk assessment data processing method according to claim 1, characterized in that, The determination of a stable reference area and the generation of a deformation rate uncertainty index based on the assessment area includes the following steps: Step 2c71, selecting a sub-area with a local slope of less than 5°, land use type of exposed bedrock or desert, and no historical geological disaster records as a stable reference area; Step 2c72, extracting the surface deformation rate of all pixels based on the stable reference area, and generating the surface deformation rate standard deviation based on the surface deformation rate; Step 2c73, obtaining the deformation rate uncertainty index based on the assessment coherence, phase residual standard deviation, and surface deformation rate standard deviation.

6. The geological hazard risk assessment data processing method according to claim 1, characterized in that, Based on the disaster-causing state variables of geographic analysis blocks, a spatial adjacency consistency test is performed. If the KL divergence of the disaster-causing state variables of adjacent geographic analysis blocks is less than the probability target value and they are in the same influence zone, then the adjacent geographic analysis blocks are merged into a comprehensive risk block. This includes the following steps: Step 41, extracting adjacent block pairs based on geographic analysis blocks; Step 42, traversing all adjacent geographic analysis block pairs, obtaining a first factor based on the disaster-causing state variables of the adjacent block pairs, and generating a KL divergence and a second factor based on the first factor; if both the KL divergence and the second factor are less than the probability target value, it means that the disaster-causing degree of the adjacent block pairs is consistent, and proceeding to the next step; otherwise, the process ends; Step 43, if the adjacent block pairs are all in the same influence zone, it is determined that the construction consistency condition is met and the adjacent block pairs are merged into a comprehensive risk block; Step 44, based on the comprehensive risk block, step 4 is re-executed to obtain the disaster-causing state update variable; Step 45, steps 41 to 44 are repeated until no adjacent geographic analysis block pairs are merged in any round of traversal.

7. The geological hazard risk assessment data processing method according to claim 6, characterized in that, The disaster-causing state variables include continuous variables and discrete variables.

8. The geological hazard risk assessment data processing method according to claim 7, characterized in that, The first factor is a continuous variable, and the second factor is the deviation value of the continuous variable.

9. A geological hazard risk assessment data processing system, used to execute the geological hazard risk assessment data processing method according to any one of claims 1-8, characterized in that, It includes the following modules: Data acquisition module: used to acquire raw observation data of the assessment area: digital elevation model data, ground rain gauge data, and surface topography monitoring data; Data processing module: used to generate multiple error metrics based on the raw observation data: elevation error probability distribution, error covariance matrix, and deformation rate uncertainty index; The partitioning module maps the original observation data and error metrics to a regular grid to obtain a grid field, constructs multiple influence zones, and divides the grid field into several geographic analysis blocks. The merging module constructs disaster-causing state variables based on the geographic analysis blocks, performs spatial adjacency consistency checks on the disaster-causing state variables of the geographic analysis blocks, and merges adjacent geographic analysis blocks into a comprehensive risk block if adjacent geographic analysis blocks are in the same influence zone and the KL divergence of the corresponding disaster-causing state variables is less than the probability target value. Risk classification module: Generates risk classification results based on the comprehensive risk blocks.

Images