A geological disaster early warning method based on remote sensing technology

Abstract

This invention discloses a geological disaster early warning method based on remote sensing technology, belonging to the field of geological disaster early warning technology. The method includes: acquiring remote sensing image data with at least two different imaging principles within the same geographic unit and the same observation time window; performing geometric correction on the remote sensing image data under a unified coordinate system; and converting the remote sensing image data into corresponding physical quantity data and completing scale consistency processing through sensor response conversion and radiometric calibration; using the remote sensing image data and digital elevation data under the unified coordinate system, calculating the slope and aspect of each pixel, and constructing the light propagation path based on the solar altitude angle and solar azimuth angle. This invention models the influence range of environmental factors from the perspective of propagation structure by constructing light propagation paths and heat diffusion propagation paths under a unified coordinate system and calculating the direction vector field formed by the gravity direction of the slope. This allows anomaly identification to be based not only on numerical changes but also on the analysis of propagation mechanisms.

Description

Technical Field

[0001] This invention relates to the field of geological disaster early warning technology, specifically a geological disaster early warning method based on remote sensing technology. Background Technology

[0002] With the widespread application of remote sensing technology in geological disaster monitoring, landslide identification and early warning based on optical and thermal infrared images has become an important technical means. Existing methods usually identify potential landslide areas by extracting surface deformation features, temperature anomaly features, or image texture change features, and combine them with time series analysis to achieve dynamic monitoring. In practical applications, influenced by factors such as changes in solar altitude angle, topographic shading, surface heat diffusion, and lithological differences, changes in sunlight and ground temperature anomalies may produce spatially similar anomalous features to landslide activity. Existing technologies typically identify anomalous areas through threshold determination or single feature analysis. However, under complex terrain and multiple environmental factors coupled together, some environmentally driven anomalies may have similar spatial morphology or numerical changes to actual landslide precursors. Therefore, under the condition of multi-source remote sensing data, how to construct the environmental propagation structure and analyze the propagation mechanism of anomalies in order to distinguish between environmentally driven anomalies and landslide-driven anomalies has become an important research direction for improving the accuracy of remote sensing geological disaster early warning; therefore, this invention proposes a geological disaster early warning method based on remote sensing technology. Summary of the Invention

[0003] The purpose of this invention is to provide a geological disaster early warning method based on remote sensing technology to solve the problems mentioned in the background art.

[0004] This invention can be achieved through the following technical solution: a geological disaster early warning method based on remote sensing technology, comprising: Step 1: Acquire remote sensing image data with at least two different imaging principles within the same geographic unit and the same observation time window, perform geometric correction on the remote sensing image data under a unified coordinate system, and convert the remote sensing image data into corresponding physical quantity data and complete the scale uniformity processing through sensor response conversion and radiometric calibration. Step 2: Using remote sensing image data and digital elevation data under a unified coordinate system, calculate the slope and aspect of each pixel, construct the light propagation path based on the solar altitude angle and solar azimuth angle, construct the heat diffusion propagation path based on the main direction of heat diffusion, and calculate the gravity direction of the slope to form a direction vector field under a unified coordinate system. Step 3: Calculate the environmentally explainable coverage area based on the light propagation path and the heat diffusion propagation path, and identify anomalous areas that have spatial differences from the environmentally explainable coverage area; Step 4: Calculate the change based on remote sensing image data, and perform residual extraction on the change under the spatial constraints of the environmentally interpretable coverage area to obtain the abnormal residual region. Calculate the minimum distance from the abnormal residual region to the nearest boundary of the light propagation path or the heat diffusion propagation path, as well as the angle between the abnormal expansion direction and the light propagation direction and the main heat diffusion direction, to obtain the abnormality detachment parameter. Step 5: When the following conditions are met simultaneously, a landslide driving anomaly is identified and an early warning is triggered: The anomaly in the abnormal region is explained as a percentage lower than a preset percentage threshold; the anomaly detachment parameter is greater than a preset distance threshold; the number of continuously connected pixels in the abnormal residual region is greater than a preset connected number threshold; and the projection component of the anomaly in the slope gravity direction is greater than the projection components in the light propagation direction and the main heat diffusion direction, respectively.

[0005] A further technical improvement of the present invention is that constructing the light propagation path includes the following steps: In the unified coordinate system after geometric correction in step one, the center point of the pixel corresponding to the digital elevation data is taken as the ray starting point, and the direction of solar incidence is determined according to the solar altitude angle and solar azimuth angle as the ray direction. The occlusion relationship is determined by comparing the height of the ray at the current pixel position with the ground elevation of the current pixel along the ray direction with a preset step size. In the absence of obstruction, calculate the angle of incidence between the sun's incident direction and the slope normal vector of the current pixel, determine the light intensity attenuation coefficient based on the angle of incidence, and record the light intensity attenuation coefficient as the light intensity value of the current pixel. During the ray calculation process, the propagation distance is accumulated and compared with the maximum propagation distance determined based on the solar altitude angle and the current pixel slope. The ray calculation is terminated when the propagation distance exceeds the maximum propagation distance. A preset angle perturbation is applied to the solar incident direction to form multiple adjacent rays. The above-mentioned occlusion judgment and light intensity attenuation coefficient calculation process is repeated for each adjacent ray. Only when each adjacent ray meets the unobstructed condition in the corresponding pixel is the pixel determined to be a light propagation path.

[0006] A further technical improvement of the present invention is that constructing the heat diffusion propagation path includes the following steps: Under the unified coordinate system after geometric correction in step one, the slope and aspect of each pixel are calculated based on the digital elevation data, and the main direction of heat diffusion is determined according to the direction of the maximum gradient on the slope. The thermal diffusion distance is calculated pixel by pixel in the main direction of thermal diffusion. For each pixel, the thermal conductivity is determined according to its lithology and surface cover type. The slope distance between the current pixel and the adjacent pixel is weighted according to the thermal conductivity to form the corrected propagation distance. The heat input energy per unit area is calculated based on the angle between the solar incident direction and the current pixel slope normal vector, and the heat input energy per unit area is accumulated by gradually attenuating along the corrected propagation distance to obtain the cumulative heat diffusion energy value. During the recursive process, the cumulative heat diffusion energy value is compared with a preset energy threshold. The recursive process is terminated when the cumulative heat diffusion energy value is less than the preset energy threshold. Simultaneously, the direction of the main diffusion axis of the heat diffusion path is calculated, and the angle between the main diffusion axis and the light propagation direction determined in step two is compared with a preset direction threshold. When the angle is greater than the preset direction threshold, the recursion is terminated, and only the pixels that simultaneously satisfy the modified propagation distance constraint, energy threshold constraint and direction consistency constraint are determined as the heat diffusion propagation path.

[0007] A further technical improvement of the present invention is that the direction of gravity of the slope is determined in the following way under a unified coordinate system: The negative direction of the maximum gradient direction on the slope is determined based on the slope and aspect of each pixel and is taken as the initial gravity direction. Calculate the aspect difference between adjacent pixels and correct the initial gravity direction based on the aspect difference; The corrected direction vector is projected onto the light propagation direction and the main heat diffusion direction respectively, and the projection components with an angle less than a preset direction threshold with the light propagation direction or the main heat diffusion direction are removed. The direction vector after removing the projection components is normalized in a unified coordinate system to form the gravity direction vector field of the slope.

[0008] A further technical improvement of the present invention is that the determination of the explanatory ratio of anomalies includes the following steps: Under a unified coordinate system, each pixel in the abnormal area is taken as the starting pixel for reverse tracing, and accessibility is determined pixel by pixel in reverse along the direction of light propagation or the main direction of heat diffusion. During the reverse determination process, adjacent pixels meet the directional consistency condition, and the light intensity value or cumulative heat diffusion energy value changes monotonically along the reverse tracing path and does not cross the boundary of the light propagation path or the boundary of the heat diffusion propagation path. When the reverse tracing path can be continuously connected to the pixel where the ray origin of the light propagation path is located or the initial pixel recursively derived from the main direction of heat diffusion, the corresponding pixel will be counted as an interpretable pixel. The ratio of the number of explainable pixels to the total number of pixels in the anomalous area is calculated, and this ratio is determined as the anomalous explainable proportion.

[0009] A further technical improvement of the present invention is that: within multiple consecutive observation time windows, the anomaly explainable ratio and the anomaly deviation parameter obtained in step four are determined respectively for each observation time window; When the explanatory proportion of anomalies decreases successively within multiple consecutive observation time windows, and the anomaly deviation parameter increases successively within multiple consecutive observation time windows, the minimum explanatory proportion of anomalies and the maximum anomaly deviation parameter within multiple consecutive observation time windows are used as the explanatory proportion of anomalies and the anomaly deviation parameter for determining landslide-driven anomalies.

[0010] A further technical improvement of the present invention is that the minimum distance from the abnormal residual region to the boundary of the light propagation path or the boundary of the heat diffusion propagation path in step four is an effective distance calculated along the slope direction and satisfying the propagation consistency constraint. The effective distance is calculated pixel by pixel in reverse along the direction of light propagation or the main direction of heat diffusion in a unified coordinate system. The recursive path must not cross the boundary of the light propagation path or the boundary of the heat diffusion path. During the recursive process, the angle between the connection direction between adjacent pixels and the light propagation direction or the main direction of heat diffusion is less than the preset direction threshold, and the light intensity value or the cumulative heat diffusion energy value changes monotonically along the recursive path. Among the paths that simultaneously satisfy the following conditions: the recursive path does not cross the boundary of the light propagation path or the boundary of the heat diffusion propagation path, the angle between the connecting directions is less than a preset direction threshold, and the light intensity value or the cumulative heat diffusion energy value changes monotonically along the recursive path, the path length with the smallest distance is selected as the minimum distance.

[0011] A further technical improvement of the present invention is that the counting of the number of continuously connected pixels in the abnormal residual region in step five includes the following limitations: Under a unified coordinate system, for a connected region formed by continuous connected pixels, calculate the maximum projected distance of the connected region in the direction of slope gravity and the maximum projected distance in the direction perpendicular to the direction of slope gravity, and calculate the ratio between the two. When the ratio is greater than the preset aspect ratio threshold, the connected region will be counted in the number of continuous connected pixels; wherein, the preset aspect ratio threshold and the preset connected number threshold are determined according to the slope, lithology type and land cover type calculated from digital elevation data.

[0012] Compared with the prior art, the present invention has the following beneficial effects: This invention constructs light propagation paths and heat diffusion propagation paths under a unified coordinate system and calculates the direction vector field formed by the gravity direction of the slope. It models the influence range of environmental factors from the perspective of propagation structure, so that anomaly identification is no longer based solely on numerical changes, but on the analysis of propagation mechanisms, thereby improving the reliability of anomaly source determination. Furthermore, this invention achieves a quantitative determination of whether an anomaly has deviated from the environmental propagation structure by calculating the interpretable coverage area of ​​the environment, the anomaly detachment degree parameter, and the projection relationship of the anomaly in the multi-directional vector field. It also combines the anomaly interpretability ratio and the number of continuously connected pixels to carry out multi-condition joint control, so that the determination of landslide-driven anomalies has a structural basis. On the other hand, this invention further enhances the ability to identify real landslide-driving anomalies by analyzing the changing trends of abnormal parameters within a time window and constraining the morphological structure of connected regions. Under conditions of complex terrain and coupling of multiple environmental factors, it helps to reduce the risk of false alarms and improve the stability and accuracy of geological disaster early warning. Attached Figure Description

[0013] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings.

[0014] Figure 1 This is a schematic diagram of the method flow of the present invention. Detailed Implementation

[0015] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided.

[0016] Please see Figure 1 As shown, the present invention provides a geological disaster early warning method based on remote sensing technology, comprising: Step 1: Acquire remote sensing image data of at least two different imaging principles within the same geographic unit and the same observation time window. Perform geometric correction on the remote sensing image data under a unified coordinate system. Convert the remote sensing image data into corresponding physical quantity data and complete scale consistency processing through sensor response conversion and radiometric calibration. By using the constraints of the same geographic unit and the same observation time window, ensure that the different remote sensing image data are comparable in time. Through geometric correction under a unified coordinate system, achieve one-to-one spatial correspondence between remote sensing image data of different imaging principles at the pixel level, and avoid false anomalies introduced by registration errors. Specifically, at least two remote sensing image data based on different imaging principles are acquired within the same geographic unit and the same observation time window. The preferred remote sensing image data includes one optical remote sensing image and one thermal infrared remote sensing image. Geometric correction is performed on the remote sensing image data under a unified coordinate system. Specifically, control points of ground features under the same coordinate datum as the digital elevation data are selected for quadratic polynomial registration and bicubic interpolation resampling to ensure that the pixel grid of the two types of remote sensing image data is consistent with that of the digital elevation data. This satisfies the requirement that "geometric correction under a unified coordinate system enables a one-to-one spatial correspondence between remote sensing image data based on different imaging principles at the pixel level." Furthermore, sensor response conversion and radiometric calibration are used to convert the remote sensing image data into corresponding physical quantity representations. Furthermore, scale uniformity processing is completed. Specifically, for optical remote sensing images, digital quantization values ​​are converted into radiance and then into surface reflectance. For thermal infrared remote sensing images, digital quantization values ​​are converted into radiance and then into brightness temperature, and surface emissivity is corrected to obtain surface temperature. This ensures that the observation values ​​output by different sensors meet the comparability requirements of differential operations and threshold determination under a clearly defined physical representation. For example, optical remote sensing images and thermal infrared remote sensing images at 10:30 can be selected within the same observation time window, and the spatial resolution is uniformly resampled to 10 meters. The residual of corresponding points after geometric correction is controlled within 0.3 pixels to ensure the stability of the correspondence between the ray origin pixel and digital elevation data when "constructing the light propagation path" in the subsequent process.

[0017] Based on the principle of "under the unified coordinate system after geometric correction in step one, using the pixel center point corresponding to the digital elevation data as the ray starting point, and determining the solar incidence direction as the ray direction based on the solar altitude angle and solar azimuth angle," the three-dimensional coordinates of each pixel center point are set as the ray starting point. The solar incidence direction is generated as a unit direction vector from the solar altitude angle and solar azimuth angle and used for pixel-by-pixel advancement. The preset step size for "calculating pixel-by-pixel along the ray direction at a preset step size" is explicitly set to one pixel center spacing, meaning that each advancement crosses one pixel grid along the ray direction. If the pixel resolution is 10 meters, the step size is 10 meters. At the k-th step of each advancement, the height of the ray at the current pixel position is calculated as follows: (Solar altitude angle), where The ray originates from the surface elevation data corresponding to the digital elevation data, where s is the preset step size. The occlusion relationship is determined by comparing the ray's height at the current pixel location with the current pixel's surface elevation. ≥ The system determines when there is obstruction and terminates the contribution of the ray to the light propagation of subsequent pixels. This represents the surface elevation value of the pixel corresponding to the k-th advance position in the digital elevation data. When there is no obstruction, the following steps are executed: "Calculate the angle of incidence between the sun's incident direction and the slope normal vector of the current pixel, determine the light intensity attenuation coefficient based on the angle of incidence, and record the light intensity attenuation coefficient as the light intensity value of the current pixel;" where the slope normal vector is obtained from the slope and aspect calculated from digital elevation data, the angle of incidence is the angle between the sun's incident direction and the slope normal vector, and the light intensity attenuation coefficient is taken as max(0, cos(angle of incidence)) and recorded as the light intensity value; Furthermore, the feasible calculation of the maximum propagation distance is performed by "accumulating the propagation distance during the ray calculation process and comparing the propagation distance with the maximum propagation distance determined based on the solar altitude angle and the current pixel slope": The slope propagation distance is the path length along the slope direction; with the ray starting point as the center, the digital elevation data is scanned pixel by pixel along the sun incident direction within the preset search steps N, and the maximum rise of the ground elevation of each pixel on the scanning path relative to the ground elevation of the ray starting point is recorded as the maximum elevation difference ΔH. The maximum horizontal propagation distance is calculated based on the solar altitude angle. This maximum horizontal propagation distance equals ΔH divided by the tangent of the solar altitude angle. Then, the maximum horizontal propagation distance is converted to the maximum slope propagation distance based on the current pixel slope. This maximum slope propagation distance equals the maximum horizontal propagation distance divided by the cosine of the current pixel slope. Therefore, we obtain: Maximum propagation distance on slope Dmax = ΔH / (tan (solar altitude angle) × cos (current pixel slope)); Subsequently, the maximum propagation distance Dmax on the slope is converted into the upper limit of the number of steps Kmax, which is equal to Dmax divided by the ray step size s and rounded down; the ray calculation is terminated when the cumulative slope propagation distance exceeds Dmax.

[0018] For example: The solar altitude angle is 45°, the current pixel slope is 30°, the step size is 10m, the scanning step size for pixel-by-pixel scanning is taken as the ray step size s, the preset search number N is 200 steps, and the corresponding maximum search distance is 2000m. If the maximum elevation difference ΔH = 150m is obtained from the scan, then: first step: The maximum horizontal propagation distance D_horizontal = 150 / tan45° = 150m; Step Two: cos30°≈0.866; The maximum propagation distance on the slope, Dmax, is approximately 173.2 m (150 / 0.866). Step 3: Kmax = ⌊173.2 / 10⌋ = 17 steps; The ray calculation is terminated when the cumulative slope propagation distance exceeds 173.2m.

[0019] Four adjacent rays are formed by applying angular perturbations of ±0.5° and ±1.0° to the solar azimuth angle. The process of occlusion judgment and light intensity attenuation coefficient calculation is repeated for each adjacent ray. Only when each adjacent ray meets the unobstructed condition in the corresponding pixel is the pixel determined to be a light propagation path, thereby reducing the influence of local elevation noise on the single ray judgment result and stabilizing the light propagation path boundary.

[0020] After completing the ray tracing and multi-ray consistency screening, a set of pixels representing the light propagation path is formed. Under a unified coordinate system, the light intensity value is synchronously output for each pixel belonging to the light propagation path. This allows for sharing the same spatial grid and directional representation when constructing the heat diffusion propagation path in step two. In step three, the light propagation path can be included as part of the environmental factor propagation structure in the calculation of the environmentally interpretable coverage area. Simultaneously, the boundary of the light propagation path is defined as the adjacency boundary between the set of pixels representing the light propagation path and the set of pixels representing non-light propagation paths. This boundary is used in step four to calculate the minimum distance from the abnormal residual region to the boundary of the light propagation path, ensuring that subsequent differential operations and the calculation of the detachment degree parameter can directly reference the light propagation path results on the same pixel grid without secondary conversion.

[0021] Step Two: Using remote sensing imagery and digital elevation data under a unified coordinate system, calculate the slope and aspect of each pixel. Construct light propagation paths based on solar altitude and azimuth angles, and heat diffusion propagation paths based on the main direction of heat diffusion. Calculate the gravity direction of the slope to form a directional vector field under the unified coordinate system. Calculating the slope and aspect of each pixel using digital elevation data provides basic parameters for characterizing the directional constraints of terrain on light shading, heat diffusion, and material sliding. Constructing light propagation paths using solar altitude and azimuth angles clearly expresses the spatial range of light propagation due to shading on the terrain surface. Constructing heat diffusion propagation paths using the main direction of heat diffusion clearly expresses the spatial extension characteristics of surface temperature affected by heat propagation. Calculating the gravity direction of the slope introduces the dominant directional reference for potential slope movement. Finally, a directional vector field under the unified coordinate system is formed, enabling the light propagation direction, the main direction of heat diffusion, and the gravity direction of the slope to participate in subsequent anomaly source determination under the same coordinate and pixel system.

[0022] Specifically, after geometric correction in step one, digital elevation data is read under the unified coordinate system. For each pixel, the east-west and north-south elevation gradients are calculated using a 3x3 neighborhood elevation difference method. The east-west gradient is calculated by dividing the elevation difference between adjacent east-west pixels by twice the pixel resolution, and the north-south gradient is calculated by dividing the elevation difference between adjacent north-south pixels by twice the pixel resolution. Then, the gradient magnitude and gradient direction are calculated based on the east-west and north-south elevation gradients. The gradient magnitude is converted into a slope value, and the gradient direction is converted into an aspect value. Based on this, the direction vector is determined according to the direction of the maximum gradient on the slope, and the slope... The direction of the maximum gradient on the slope serves as the common directional reference for subsequent steps such as "determining the main direction of heat diffusion based on the direction of the maximum gradient on the slope" and "determining the negative direction of the direction of the maximum gradient on the slope based on the slope and aspect of each pixel as the initial gravity direction." For example, if the pixel resolution is 10m and a certain pixel is calculated to have a slope of 28° and an aspect of 135° through a three-by-three neighborhood, then the direction of the maximum gradient on the slope can be determined to be the 135° direction, and the corresponding negative direction is the 315° direction. This ensures that the subsequent heat diffusion propagation path and the slope gravity direction are defined under the same pixel coordinates and the same direction definition system.

[0023] Under the unified coordinate system after geometric correction in step one, the slope and aspect of each pixel are calculated based on digital elevation data, and the main direction of heat diffusion is determined according to the direction of the maximum gradient on the slope. After determining the main direction of heat diffusion for each pixel, the calculation is performed pixel by pixel along the main direction of heat diffusion. For each pixel, the thermal conductivity is determined according to its lithology and land cover type. The thermal conductivity is determined through two levels: a preset lithology-thermal conductivity correspondence table and a land cover type-correction coefficient correspondence table. The lithology type is mapped from the rasterized geological map vector data to the pixels under the unified coordinate system, and the land cover type is mapped from the optical remote sensing image classification results to the corresponding pixels. Example values ​​are: granite 2.6 W / (m·K), shale 1.8 W / (m·K), sandstone 2.2 W / (m·K), bare land correction coefficient 1.0, grassland 0.85, and forest 0.75. The slope distance between the current pixel and its adjacent pixels is weighted according to the thermal conductivity to form a corrected propagation distance. The slope distance is calculated by dividing the length of the line connecting the centers of adjacent pixels by the cosine of the slope. If the pixel resolution is 10m and the slope is 30°, the slope distance is approximately 10 / cos(30°) ≈ 11.55m. The corrected propagation distance is calculated by dividing the slope distance by the normalized value of the thermal conductivity, thus increasing the corrected propagation distance in areas with lower thermal conductivity. Subsequently, the heat input energy per unit area is calculated based on the angle between the solar incidence direction and the slope normal vector of the current pixel. The solar incidence direction is calculated by generating a unit vector from the solar altitude angle and the solar azimuth angle, and the slope normal vector is calculated from the slope and aspect. The heat input energy per unit area is calculated as incident irradiance × incident angle cosine × observation time window duration. For example, with an incident irradiance of 800W / m², an incident angle cosine of 0.7, and an observation time window of 6... If the heat input energy per unit area is 336000 J / m², then the cumulative heat diffusion energy value is calculated by progressively attenuating the heat input energy value per unit area along the corrected propagation distance. The attenuation coefficient is set as exp(-β×corrected propagation distance), where β is the attenuation parameter after normalization of thermal conductivity. For example, β=0.02m⁻¹. Finally, the cumulative heat diffusion energy value is compared with the preset energy threshold during the recursive process. The recursive process is terminated when the cumulative heat diffusion energy value is less than the preset energy threshold. The preset energy threshold is obtained by selecting a historically stable area within the same geographic unit as the calibration area, statistically analyzing the mean and standard deviation of the cumulative heat diffusion energy value within multiple observation time windows in the calibration area, and taking the mean minus twice the standard deviation as the preset energy threshold. For example, if the mean is 120000 J / m² and the standard deviation is 15000 J / m², then the preset energy threshold is 90000 J / m².

[0024] After completing the recursive calculation, the "simultaneous calculation of the diffusion principal axis direction of the heat diffusion path" is performed on the set of pixels along the heat diffusion path. The diffusion principal axis direction is obtained by weighting the direction of the line connecting the centers of adjacent pixels in the heat diffusion path according to the cumulative heat diffusion energy value. Then, the angle between the diffusion principal axis direction and the light propagation direction determined in step two is compared with a preset direction threshold. When the angle is greater than the preset direction threshold, the recursion is terminated. The angle is calculated by multiplying the unit vector by the inverse cosine. The preset direction threshold is determined by the 90th percentile of the distribution of the angle between the diffusion principal axis direction and the light propagation direction in the calibration area, for example, 25°. Finally, only pixels that simultaneously satisfy the modified propagation distance constraint, energy threshold constraint, and direction consistency constraint are identified as heat diffusion propagation paths. The heat diffusion propagation path is then output.

[0025] The negative direction of the maximum gradient direction on the slope is determined based on the slope and aspect of each pixel and taken as the initial gravity direction. After obtaining the initial gravity direction, the aspect difference between adjacent pixels is calculated, and the initial gravity direction is corrected based on the aspect difference. The direction correction is achieved by weighting and normalizing the absolute values ​​of the aspect differences of eight neighboring pixels. Subsequently, the corrected direction vector is projected onto the light propagation direction and the principal direction of heat diffusion, respectively, and the projection components are obtained by vector dot product. Then, projection components with an angle less than a preset direction threshold with the light propagation direction or the principal direction of heat diffusion are removed. The angle calculation is the same as above, and the preset direction threshold remains 25°. Finally, the direction vector after removing the projection components is normalized in a unified coordinate system to form a slope gravity direction vector field. Normalization is achieved by dividing the direction vector by its magnitude. If the magnitude is less than 1×10⁻ 6 The initial gravity direction unit vector is then used as a substitute, thereby forming a slope gravity direction vector field under a unified coordinate system, which together with the light propagation direction and the main heat diffusion direction constitutes a direction vector field under a unified coordinate system.

[0026] Step 3: Calculate the interpretable environmental coverage area based on the light propagation path and the heat diffusion propagation path, and identify anomalous areas that spatially differ from the interpretable environmental coverage area. Utilizing the propagation range and direction constraints described by the light propagation path and the heat diffusion propagation path, calculate the interpretable environmental coverage area to clearly define the spatial impact range that may be caused by changes in light shading or heat diffusion. Based on this, spatially correlate the observed changes with the interpretable environmental coverage area to identify anomalous areas that spatially differ from the interpretable environmental coverage area. This initially highlights anomalies that may be caused by geological processes such as landslides from the interpretable environmental propagation range, providing a clear spatial object for subsequent differential operations to extract anomalous residual areas and calculate anomalous deviation parameters.

[0027] Specifically, under a unified coordinate system, the obtained light propagation paths and heat diffusion propagation paths are read. The set of pixels covered by the light propagation paths and the set of pixels covered by the heat diffusion propagation paths are combined to obtain a candidate coverage set. Within the candidate coverage set, the light intensity value corresponding to the light propagation path and the cumulative heat diffusion energy value corresponding to the heat diffusion propagation path are used as environmental impact intensity constraints. Threshold-based screening is applied to form environmentally interpretable coverage areas. The light intensity value threshold is obtained by selecting stable areas without geological disaster records within the same geographic unit, statistically analyzing the distribution of light intensity values ​​within the stable area, and taking the 10th percentile as the lower limit threshold. The cumulative heat diffusion energy value threshold is obtained by statistically analyzing the distribution of cumulative heat diffusion energy values ​​within the stable area and taking the 10th percentile as the lower limit threshold. For example, if the 10th percentile of the light intensity value in the stable area is 0.18, and the cumulative heat diffusion energy value is 10... If the percentile is 90,000 J / m², only candidate coverage pixels with an illuminance value of not less than 0.18 or a cumulative thermal diffusion energy value of not less than 90,000 J / m² are retained as environmentally interpretable coverage areas. Based on this, the change in remote sensing image data is calculated to characterize the "observed change". For optical remote sensing images, the change is the difference in surface reflectance between two adjacent periods within the same observation time window. For thermal infrared remote sensing images, the change is the difference in surface temperature between two adjacent periods. For each change map, the mean plus three times the standard deviation of the stable area within the same geographic unit is used as the change threshold to extract areas of significant change. The areas of significant change are spatially correlated with the environmentally interpretable coverage areas. The set of pixels in the areas of significant change that fall outside the environmentally interpretable coverage areas is identified as anomaly areas with spatial differences from the environmentally interpretable coverage areas, thus providing a clear object for the selection of starting pixels for subsequent reverse tracing.

[0028] Under a unified coordinate system, each pixel within the anomaly region is used as the starting pixel for reverse tracing. Reachability is determined pixel by pixel along either the light propagation direction or the main heat diffusion direction. The step size for reverse pixel-by-pixel tracing is one pixel center-to-center distance; for a pixel resolution of 10m, the step size is 10m. During the reverse tracing process, adjacent pixels must meet a directional consistency condition. This condition is defined as the angle between the line connecting the centers of adjacent pixels and the light propagation direction or the main heat diffusion direction being less than a preset directional threshold. This preset directional threshold is determined by the 90th percentile of the angle distribution between the line connecting the centers of adjacent pixels and the light propagation direction within the stable region, for example, 25°. Simultaneously, immutable constraints are imposed on monotonic changes: when tracing back along the light propagation direction, the light intensity value along the reverse tracing path must not decrease progressively, meaning the light intensity value of the next pixel is not less than the light intensity value of the previous pixel; when tracing back along the main heat diffusion direction, the reverse... The cumulative thermal diffusion energy value along the tracing path does not decrease progressively, meaning the cumulative thermal diffusion energy value of the subsequent pixel is not less than the cumulative thermal diffusion energy value of the preceding pixel; and it does not cross the boundary of the light propagation path or the boundary of the thermal diffusion propagation path. "Crossing" is defined as the reverse tracing path jumping from a set of pixels belonging to the light propagation path to a set of pixels not belonging to the light propagation path, or from a set of pixels belonging to the thermal diffusion propagation path to a set of pixels not belonging to the thermal diffusion propagation path. When the reverse tracing path can be continuously connected to the pixel where the ray origin of the light propagation path is located or the initial pixel recursively in the main direction of thermal diffusion, the corresponding pixel is counted as an interpretable pixel. Finally, the ratio of the number of interpretable pixels to the total number of pixels in the anomaly area is calculated, and this ratio is determined as the anomaly interpretability ratio. For example, if the anomaly area contains a total of 500 pixels, of which 120 satisfy the reverse tracing accessibility criteria and are counted as interpretable pixels, then the anomaly interpretability ratio is 0.24.

[0029] Within multiple consecutive observation time windows, maintaining the same geographical unit and processing flow, the corresponding anomalous area and anomaly explainable proportion are repeatedly obtained for each observation time window, and the anomaly deviation parameter obtained in step four is read. The observation time window can be set as a 10-minute window centered on a fixed time within the same day and rolled over multiple days. For example, three consecutive observation time windows are 10:30 on day 1, 10:30 on day 2, and 10:30 on day 3. For each observation time window, the anomaly explainable proportion sequence and the anomaly deviation parameter sequence are calculated. When the anomaly explainable proportion successively increases within multiple consecutive observation time windows... When the anomaly deviation parameter decreases and increases successively within multiple consecutive observation time windows, the minimum anomaly explainable ratio and the maximum anomaly deviation parameter within these multiple observation time windows are used as the anomaly explainable ratio and anomaly deviation parameter for determining landslide-driven anomalies. For example, if the anomaly explainable ratios within three observation time windows are 0.31, 0.27, and 0.22, and the anomaly deviation parameters are 1.6, 1.9, and 2.3, then the minimum anomaly explainable ratio of 0.22 and the maximum anomaly deviation parameter of 2.3 are taken as the subsequent judgment inputs, thereby strengthening the suppression of environmentally driven fluctuations in the time dimension.

[0030] Finally, after obtaining the explanatory proportion and anomaly deviation parameters used for judgment, step five of weight 1 is followed to simultaneously test whether the explanatory proportion of the anomaly region is lower than a preset proportion threshold, the anomaly deviation parameter is greater than a preset distance threshold, the number of continuously connected pixels in the anomaly residual region is greater than a preset connected number threshold, and the projection component of the anomaly in the slope gravity direction is greater than the projection components in the light propagation direction and the main heat diffusion direction, respectively. The preset proportion threshold is determined by taking the 10th percentile of the historical distribution of the explanatory proportion of the anomaly within the stable region, for example, 0.30; the preset distance threshold is determined by taking the 90th percentile of the historical distribution of the anomaly deviation parameter within the stable region, for example, 1.8. The preset connectivity threshold is determined by taking the 95th percentile of the number of connected pixels in the connected noise patch within the stable region, for example, 60 pixels. When the anomaly explainable ratio is 0.22 and lower than 0.30, the anomaly deviation parameter is 2.3 and greater than 1.8, the largest connected patch in the anomaly residual region contains 85 pixels and is greater than 60, and the projection component of the anomaly expansion direction in the slope gravity direction is greater than the projection components in the light propagation direction and the main heat diffusion direction, respectively, it is determined to be a landslide-driven anomaly and an early warning is triggered. Thus, "steps three through four" and the two subordinate constraints form a continuous, reproducible, and calculable complete chain in the same embodiment.

[0031] Step 4: Calculate the variation based on remote sensing image data, and perform residual extraction on the variation under the spatial constraints of the environmentally interpretable coverage area to obtain the anomalous residual region. Calculate the minimum distance from the anomalous residual region to the nearest boundary of the light propagation path or heat diffusion propagation path, as well as the angle between the anomalous expansion direction and the main directions of light propagation and heat diffusion, to obtain the anomalous deviation parameter. By calculating the variation based on remote sensing image data and performing residual extraction on the variation under the spatial constraints of the environmentally interpretable coverage area, the anomalous residual region is obtained, making the residual parts in the anomalous region that cannot be explained by the environmentally interpretable coverage area explicitly separated. The process involves several steps: first, determining the distance between the anomaly residual region and the nearest boundary of either the light propagation path or the heat diffusion propagation path; second, calculating the minimum distance between the anomaly residual region and the nearest boundary of either the light propagation path or the heat diffusion propagation path, to measure the degree to which the anomaly residual region deviates from the range of environmental propagation in space; and third, calculating the angle between the anomaly expansion direction and the light propagation direction and the main heat diffusion direction, to measure whether the directional characteristics of the anomaly expansion are consistent with the direction of environmental propagation. Finally, a parameter representing the degree of anomaly detachment is formed, providing a quantifiable basis for determining whether the anomaly originates from environmental propagation or whether it is more consistent with the evolution characteristics of the slope itself, and providing direct input for the threshold determination in step five.

[0032] Specifically, remote sensing image data and environmentally interpretable coverage areas are read under a unified coordinate system, maintaining a one-to-one correspondence between pixels. The process of calculating changes and extracting residuals is performed according to the principle of "calculating changes based on remote sensing image data and performing residual extraction on changes under the spatial constraints of environmentally interpretable coverage areas to obtain abnormal residual areas". For optical remote sensing images, the change map is calculated using the difference in surface reflectance between two adjacent periods, and for thermal infrared remote sensing images, the change map is calculated using the difference in surface temperature between two adjacent periods. Subsequently, using the environmentally explainable coverage area as a spatial mask in the change map, pixels falling within the environmentally explainable coverage area are marked as environmentally explainable change pixels and removed, while only pixels located outside the environmentally explainable coverage area are retained as candidate anomalous pixels. Based on this, a threshold is determined for the change in candidate abnormal pixels. The threshold is determined by selecting historically stable areas within the same geographic unit as calibration areas, statistically analyzing the mean and standard deviation of the change in the calibration areas, and taking the mean plus three times the standard deviation as the threshold. For example, if the mean value of the optical change calibration area is 0.012 and the standard deviation is 0.006, then the threshold value is 0.030; if the mean value of the thermal infrared change calibration area is 0.35K and the standard deviation is 0.25K, then the threshold value is 1.10K. When the change in a candidate anomalous pixel exceeds the corresponding threshold, the pixel is identified as an anomalous residual region pixel, thus obtaining a set of pixels in the anomalous residual region, providing a clear spatial object for subsequent minimum distance and angle calculations.

[0033] Under a unified coordinate system, the boundary of the light propagation path is determined by the adjacency boundary between the set of pixels along the light propagation path and the set of pixels along the non-light propagation path, and the boundary of the heat diffusion propagation path is determined by the adjacency boundary between the set of pixels along the heat diffusion propagation path and the set of pixels along the non-heat diffusion propagation path. Then, based on the minimum distance from the abnormal residual region to the boundary of the light propagation path or the boundary of the heat diffusion propagation path in step four, the effective distance calculated along the slope direction and satisfying the propagation consistency constraint is taken. Each pixel in the abnormal residual region is taken as the starting pixel, and the path is calculated pixel by pixel in reverse along the light propagation direction or the main direction of heat diffusion. The path length is defined as the cumulative sum of the slope distance between the centers of adjacent pixels. The slope distance is obtained by dividing the horizontal distance between the centers of pixels by the cosine of the slope. For example, when the pixel resolution is 10m and the slope is 30°, the slope distance between adjacent pixels is 10 / cos30°≈11.55m, which is used for the cumulative path length.

[0034] The effective distance is calculated pixel-by-pixel along the reverse recursive path in a unified coordinate system. Each pixel in the recursive path should belong to the corresponding set of pixels for light propagation or heat diffusion propagation. During the recursion, the angle between the direction of the line connecting the centers of adjacent pixels and the direction of light propagation or the main direction of heat diffusion is calculated by multiplying the unit vector by the inverse cosine and compared with a preset direction threshold. The preset direction threshold is determined by the 90th percentile of the distribution of the angle between the connecting direction and the propagation direction in the calibration area, for example, 25°. At the same time, it is required that the light intensity value or the cumulative heat diffusion energy value along the recursive path changes monotonically, that is, the light intensity value gradually does not decrease when recursing in reverse along the direction of light propagation, and the cumulative heat diffusion energy value gradually does not decrease when recursing in reverse along the main direction of heat diffusion.

[0035] For each pixel in the abnormal residual region, calculate the path length that satisfies the above propagation consistency constraint, and select the minimum path length as the minimum distance from the pixel to the boundary of the light propagation path or the boundary of the heat diffusion propagation path.

[0036] In this embodiment, the abnormal detachment degree parameter is defined as: Abnormal detachment degree parameter = minimum distance × sin (angle); The included angle is the minimum angle between the direction of abnormal expansion, the direction of light propagation, and the main direction of heat diffusion.

[0037] Among the paths that simultaneously satisfy the following conditions: the recursive path does not cross the boundary of the light propagation path or the boundary of the heat diffusion propagation path, the angle between the connecting directions is less than a preset direction threshold, and the light intensity value or the cumulative heat diffusion energy value changes monotonically along the recursive path, the path length with the smallest distance is selected as the minimum distance. When multiple reachable paths exist for a pixel in an abnormal residual region, the length of each path is calculated, and the minimum value is taken as the minimum distance from the pixel to the nearest boundary of the light propagation path or the heat diffusion propagation path. Then, the minimum value of the minimum distance is taken within the abnormal residual region, or the minimum distance distribution is obtained by statistical analysis of pixels. For example, if a pixel in an abnormal residual region is recursively derived along the light propagation direction and two candidate paths are obtained, with path lengths of 92.4m and 115.5m respectively, then the minimum distance of the pixel is taken as 92.4m. At the same time, the angle is calculated according to another part of step four, "and the angle between the abnormal expansion direction and the light propagation direction and the main direction of heat diffusion". The abnormal expansion direction is obtained by calculating the main axis direction of the connected patches in the abnormal residual region. The main axis direction is determined by the main direction feature vector of the pixel coordinates within the connected patches. Finally, the minimum distance and the angle are used together as the calculation input for the abnormality detachment degree parameter, and provide the parameter basis for the threshold determination in step five.

[0038] Step 5: When the following conditions are met simultaneously, a landslide driving anomaly is identified and an early warning is triggered: The anomaly region exhibits the following characteristics: the anomaly explainable proportion is below a preset threshold; the anomaly detachment parameter is above a preset distance threshold; the number of continuously connected pixels in the anomaly residual region is above a preset connectivity threshold; and the anomaly's projection component in the slope gravity direction is greater than its projection components in the light propagation direction and the main heat diffusion direction, respectively. By utilizing the anomaly explainable proportion being below the preset threshold, the overall explainability of the anomaly within the environmentally explainable coverage area is constrained, making the anomaly more likely to be classified as non-environmentally propagated. Furthermore, by utilizing the anomaly detachment parameter being above the preset distance threshold, the spatial location and directional characteristics of the anomaly residual region are constrained to influence the environmental propagation structure. The degree of separation is reduced, decreasing the probability of false anomalies caused by environmental propagation being misjudged; the number of continuously connected pixels in the anomaly residual region is greater than the preset connected number threshold, constraining the anomaly to have sufficient spatial continuity to eliminate scattered noise or local unstable disturbances; the projection component of the anomaly in the slope gravity direction is greater than the projection components in the light propagation direction and the main heat diffusion direction, respectively, constraining the anomaly expansion to be more consistent with the evolution direction of the slope under gravity rather than the propagation direction of light or heat diffusion; when the above conditions are met simultaneously, an early warning is triggered, so that the final output of the landslide-driven anomaly judgment is supported in terms of spatial structure, propagation mechanism and directional consistency.

[0039] Specifically, the set of pixels in the abnormal residual region is read under a unified coordinate system, and the abnormal residual region is marked with connected components using the eight-neighbor connectivity rule to obtain a connected region formed by continuously connected pixels. The eight-neighbor connectivity rule is defined as follows: if two pixels are adjacent in the row and column direction or diagonal direction and both belong to the abnormal residual region, they are determined to be connected. For each connected region, the number of pixels it contains, the set of pixel center coordinates, and the slope gravity direction to which it belongs are recorded. The slope gravity direction is calculated in step two and has formed a direction vector field, thereby ensuring that the subsequent calculation of the "maximum projection distance" and the comparison of direction projection are completed under the same coordinate and the same pixel scale.

[0040] The counting of continuous connected pixels in the abnormal residual region in step five includes the following limitations: Under a unified coordinate system, for the connected region formed by continuous connected pixels, calculate the maximum projection distance of the connected region in the slope gravity direction and the maximum projection distance in the direction perpendicular to the slope gravity direction, and calculate the ratio between the two; when the ratio is greater than a preset aspect ratio threshold, the connected region is counted in the number of continuous connected pixels; the maximum projection distance is calculated as follows: project the line vector connecting the center points of any two pixels in the connected region in the slope gravity direction and the direction perpendicular to the slope gravity direction, and take the projection... The maximum length is taken as the maximum projection distance in the corresponding direction. To facilitate reproduction, the maximum projection distance of a connected region in the direction of slope gravity can be obtained by traversing the boundary pixel pairs of the connected region and calculating the projection length. The maximum projection distance in the direction perpendicular to the direction of slope gravity is obtained in the same way. For example, if the maximum projection distance of a connected region in the direction of slope gravity is 210m and the maximum projection distance in the direction perpendicular to the direction of slope gravity is 55m, then the ratio of the two is 3.82. If the preset aspect ratio threshold is 3.00, then the connected region meets the inclusion condition and is included in the number of continuous connected pixels.

[0041] The preset aspect ratio threshold and preset connectivity threshold are determined based on the slope, lithology, and land cover type calculated from digital elevation data. Specifically, within the same geographic unit, a historically stable area is selected as the calibration area. Pixels within the calibration area are grouped according to slope, lithology, and land cover type. The slope is calculated from digital elevation data and divided into three groups according to fixed grouping rules: less than 15°, 15° to 30°, and greater than 30°. The lithology type is mapped to pixels from a geological map using rasterization and existing classification coding. The land cover type is mapped to pixels from optical remote sensing imagery and existing classification coding is used. For each group, the empirical distribution of aspect ratio and the number of connected pixels is calculated, and the preset aspect ratio threshold is set to the 90th percentile of the aspect ratio distribution for that group. The preset connectivity threshold is set to the 90th percentile of the connectivity distribution for that group. The threshold values ​​are determined by the 95th percentile of the number of connected pixels. For example, in a group with a slope of 15° to 30°, lithology of shale, and land cover of grassland, the 90th percentile of the aspect ratio is 2.80 and the 95th percentile of the number of connected pixels is 60. Therefore, the preset aspect ratio threshold is 2.80 and the preset connected pixel threshold is 60. In a group with a slope greater than 30°, lithology of granite, and land cover of bare land, the 90th percentile of the aspect ratio is 3.20 and the 95th percentile of the number of connected pixels is 80. Therefore, the preset aspect ratio threshold is 3.20 and the preset connected pixel threshold is 80. This ensures that the threshold values ​​have a consistent statistical basis for determination as the slope, lithology, and land cover change.

[0042] A landslide-driven anomaly is identified and an early warning is triggered when the following conditions are met simultaneously: the anomaly's explainable proportion in the anomaly area is lower than a preset proportion threshold; the anomaly's detachment degree parameter is greater than a preset distance threshold; the number of continuously connected pixels in the anomaly residual area is greater than a preset connected number threshold; and the anomaly's projection component in the slope's gravity direction is greater than its projection components in the light propagation direction and the main heat diffusion direction, respectively. The number of continuously connected pixels in the abnormal residual region is based on the statistical results after the aforementioned limitation. That is, only connected regions that meet the condition of "ratio greater than the preset aspect ratio threshold" are included before being summarized. For example, in a group with a slope of 15° to 30°, lithology of shale, and surface cover of grassland, the preset aspect ratio threshold is 2.80 and the preset connected number threshold is 60. If there are two connected regions in the abnormal residual region, connected region one with an aspect ratio of 3.10 and 72 pixels, and connected region two with an aspect ratio of 2.10 and 95 pixels, then only connected region one is included in the number of continuously connected pixels and meets the judgment requirement of "greater than the preset connected number threshold". On this basis, the statistical results are judged simultaneously with the anomaly explainable ratio, anomaly deviation degree parameter, and directional projection component comparison conditions. Thus, without changing the judgment logic of step five, the ability to distinguish between noise clusters and strip-shaped anomalies is improved by limiting the statistical method of the number of continuously connected pixels, and the implementation description of the landslide-driven anomaly triggering early warning is completed.

[0043] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0044] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A geological disaster early warning method based on remote sensing technology, characterized in that, include: Step 1: Acquire remote sensing image data with at least two different imaging principles within the same geographic unit and the same observation time window, perform geometric correction on the remote sensing image data under a unified coordinate system, and convert the remote sensing image data into corresponding physical quantity data and complete the scale uniformity processing through sensor response conversion and radiometric calibration. Step 2: Using remote sensing image data and digital elevation data under a unified coordinate system, calculate the slope and aspect of each pixel, construct the light propagation path based on the solar altitude angle and solar azimuth angle, construct the heat diffusion propagation path based on the main direction of heat diffusion, and calculate the gravity direction of the slope to form a direction vector field under a unified coordinate system. Step 3: Calculate the environmentally explainable coverage area based on the light propagation path and the heat diffusion propagation path, and identify anomalous areas that have spatial differences from the environmentally explainable coverage area; Step 4: Calculate the change based on remote sensing image data, and perform residual extraction on the change under the spatial constraints of the environmentally interpretable coverage area to obtain the abnormal residual region. Calculate the minimum distance from the abnormal residual region to the nearest boundary of the light propagation path or the heat diffusion propagation path, as well as the angle between the abnormal expansion direction and the light propagation direction and the main heat diffusion direction, to obtain the abnormality detachment parameter. Step 5: When the following conditions are met simultaneously, a landslide driving anomaly is identified and an early warning is triggered: The anomaly in the abnormal region is explained as a percentage lower than a preset percentage threshold; the anomaly detachment parameter is greater than a preset distance threshold; the number of continuously connected pixels in the abnormal residual region is greater than a preset connected number threshold; and the projection component of the anomaly in the slope gravity direction is greater than the projection components in the light propagation direction and the main heat diffusion direction, respectively.

2. The geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, Constructing a light propagation path involves the following steps: In the unified coordinate system after geometric correction in step one, the center point of the pixel corresponding to the digital elevation data is taken as the ray starting point, and the direction of solar incidence is determined according to the solar altitude angle and solar azimuth angle as the ray direction. The occlusion relationship is determined by comparing the height of the ray at the current pixel position with the ground elevation of the current pixel along the ray direction with a preset step size. In the absence of obstruction, calculate the angle of incidence between the sun's incident direction and the slope normal vector of the current pixel, determine the light intensity attenuation coefficient based on the angle of incidence, and record the light intensity attenuation coefficient as the light intensity value of the current pixel. During the ray calculation process, the propagation distance is accumulated and compared with the maximum propagation distance determined based on the solar altitude angle and the current pixel slope. The ray calculation is terminated when the propagation distance exceeds the maximum propagation distance. A preset angle perturbation is applied to the solar incident direction to form multiple adjacent rays. The above-mentioned occlusion judgment and light intensity attenuation coefficient calculation process is repeated for each adjacent ray. Only when each adjacent ray meets the unobstructed condition in the corresponding pixel is the pixel determined to be a light propagation path.

3. The geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, Constructing a heat diffusion propagation path involves the following steps: Under the unified coordinate system after geometric correction in step one, the slope and aspect of each pixel are calculated based on the digital elevation data, and the main direction of heat diffusion is determined according to the direction of the maximum gradient on the slope. The thermal diffusion distance is calculated pixel by pixel in the main direction of thermal diffusion. For each pixel, the thermal conductivity is determined according to its lithology and surface cover type. The slope distance between the current pixel and the adjacent pixel is weighted according to the thermal conductivity to form the corrected propagation distance. The heat input energy per unit area is calculated based on the angle between the solar incident direction and the current pixel slope normal vector, and the heat input energy per unit area is accumulated by gradually attenuating along the corrected propagation distance to obtain the cumulative heat diffusion energy value. During the recursive process, the cumulative heat diffusion energy value is compared with a preset energy threshold. The recursive process is terminated when the cumulative heat diffusion energy value is less than the preset energy threshold. Simultaneously, the direction of the main diffusion axis of the heat diffusion path is calculated, and the angle between the main diffusion axis and the light propagation direction determined in step two is compared with a preset direction threshold. When the angle is greater than the preset direction threshold, the recursion is terminated, and only the pixels that simultaneously satisfy the modified propagation distance constraint, energy threshold constraint and direction consistency constraint are determined as the heat diffusion propagation path.

4. The geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, The direction of gravity on a slope is determined in a unified coordinate system as follows: The negative direction of the maximum gradient direction on the slope is determined based on the slope and aspect of each pixel and is taken as the initial gravity direction. Calculate the aspect difference between adjacent pixels and correct the initial gravity direction based on the aspect difference; The corrected direction vector is projected onto the light propagation direction and the main heat diffusion direction respectively, and the projection components with an angle less than a preset direction threshold with the light propagation direction or the main heat diffusion direction are removed. The direction vector after removing the projection components is normalized in a unified coordinate system to form the gravity direction vector field of the slope.

5. A geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, Determining the explanatory percentage of anomalies includes the following steps: Under a unified coordinate system, each pixel in the abnormal area is taken as the starting pixel for reverse tracing, and accessibility is determined pixel by pixel in reverse along the direction of light propagation or the main direction of heat diffusion. During the reverse determination process, adjacent pixels meet the directional consistency condition, and the light intensity value or cumulative heat diffusion energy value changes monotonically along the reverse tracing path and does not cross the boundary of the light propagation path or the boundary of the heat diffusion propagation path. When the reverse tracing path can be continuously connected to the pixel where the ray origin of the light propagation path is located or the initial pixel recursively derived from the main direction of heat diffusion, the corresponding pixel will be counted as an interpretable pixel. The ratio of the number of explainable pixels to the total number of pixels in the anomalous area is calculated, and this ratio is determined as the anomalous explainable proportion.

6. A geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, Within multiple consecutive observation time windows, the anomaly explainable proportion and the anomaly deviance parameter obtained in step four are determined for each observation time window. When the explanatory proportion of anomalies decreases successively within multiple consecutive observation time windows, and the anomaly deviation parameter increases successively within multiple consecutive observation time windows, the minimum explanatory proportion of anomalies and the maximum anomaly deviation parameter within multiple consecutive observation time windows are used as the explanatory proportion of anomalies and the anomaly deviation parameter for determining landslide-driven anomalies.

7. A geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, The minimum distance from the abnormal residual region to the boundary of the light propagation path or the boundary of the heat diffusion propagation path in step four is the effective distance calculated along the slope direction and satisfying the propagation consistency constraint. The effective distance is calculated pixel by pixel in reverse along the direction of light propagation or the main direction of heat diffusion in a unified coordinate system. The recursive path must not cross the boundary of the light propagation path or the boundary of the heat diffusion path. During the recursive process, the angle between the connection direction between adjacent pixels and the light propagation direction or the main direction of heat diffusion is less than the preset direction threshold, and the light intensity value or the cumulative heat diffusion energy value changes monotonically along the recursive path. Among the paths that simultaneously satisfy the following conditions: the recursive path does not cross the boundary of the light propagation path or the boundary of the heat diffusion propagation path, the angle between the connecting directions is less than a preset direction threshold, and the light intensity value or the cumulative heat diffusion energy value changes monotonically along the recursive path, the path length with the smallest distance is selected as the minimum distance.

8. A geological disaster early warning method based on remote sensing technology according to claim 1, characterized in that, The counting of the number of continuously connected cells in the abnormal residual region in step five includes the following limitations: Under a unified coordinate system, for a connected region formed by continuous connected pixels, calculate the maximum projected distance of the connected region in the direction of slope gravity and the maximum projected distance in the direction perpendicular to the direction of slope gravity, and calculate the ratio between the two. When the ratio is greater than the preset aspect ratio threshold, the connected region will be counted in the number of continuous connected pixels; wherein, the preset aspect ratio threshold and the preset connected number threshold are determined according to the slope, lithology type and land cover type calculated from digital elevation data.

Images