Snow surface rain disaster process simulation method and system fusing hydrothermal coupling model

Abstract

This application provides a method and system for simulating the disaster-causing process of snow surface rain using a hydrothermal coupling model. First, it acquires an initial snow layer state data set before a rainfall event in the area to be simulated. Then, it performs a snow surface rain infiltration process simulation on the initial snow layer state data set, generating data on the distribution of liquid water migration flux and the solid-ice phase transition rate within the snow layer during the simulated rainfall period. Based on these two types of distribution data, it performs a snow layer mechanical stability assessment, generating data on the distribution of pore water pressure field and effective stress field within the snow layer. Based on the pore water pressure field and effective stress field distribution data, it performs a disaster-causing process identification process, generating data on the distribution of disaster-causing risk areas and the evolution path of disaster risks. Finally, it generates a simulation result data set based on the disaster-causing risk area distribution data and the evolution path data. This invention can achieve a systematic and accurate simulation of the disaster-causing process of snow surface rain.

Description

Technical Field

[0001] This application relates to the field of snow disaster risk assessment technology, and more specifically, to a method and system for simulating snow surface rain-induced disaster processes by integrating a hydrothermal coupling model. Background Technology

[0002] When snow-covered areas are affected by rainfall events, rainwater can seep into the snow layer, causing a rapid decrease in the structural strength of the snow layer and triggering instability and landslides. Such disasters occur frequently in high-latitude regions and high-altitude mountainous areas, seriously affecting road traffic, urban facilities, and the safety of people.

[0003] Currently, most simulation studies on rainfall-induced disasters in snow-covered areas adopt a step-by-step approach: some studies focus on analyzing the infiltration behavior of rainwater in the snow layer, usually simplifying the snow layer as a homogeneous porous medium and describing the vertical infiltration of rainwater from the snow surface downwards using only a single water transport equation, without considering the temperature field changes caused by heat conduction inside the snow layer and the resulting ice crystal melting and refreezing processes; other studies focus on determining whether the snow layer has reached a critical state of instability, often directly applying strength criteria from geotechnical mechanics for static calculations under known snow layer physical parameters, without incorporating the impact of continuous changes in moisture and heat inside the snow layer on mechanical parameters during rainfall into the analysis framework; still other studies attempt to jointly calculate multiple physical processes, but usually use the method of solving each equation independently and then transferring the results, with the mutual constraints between the various physical fields only approximated through step-by-step iteration, making it difficult to reflect the synchronous response characteristics of water, heat, and force within the snow layer at the same spatiotemporal scale. Summary of the Invention

[0004] In view of this, the purpose of this application is to provide a method and system for simulating snow surface rain-induced disaster processes by integrating a hydrothermal coupling model.

[0005] In conjunction with the first aspect of this application, a method for simulating snow surface rain-induced disaster processes by incorporating a hydrothermal coupling model is provided, applied to a snow surface rain-induced disaster process simulation system by incorporating a hydrothermal coupling model, the method comprising: Obtain the initial snow layer state data set of the area to be simulated before the occurrence of the rainfall event. The initial snow layer state data set includes spatial distribution data of snow layer thickness, spatial distribution data of snow layer temperature, and spatial distribution data of snow layer liquid water content. The initial snow layer state data set is subjected to snow surface rainfall infiltration process simulation processing to generate data on the distribution of liquid phase water migration flux inside the snow layer and the distribution of solid phase ice phase transformation rate inside the snow layer during the simulated rainfall period. Based on the data on the distribution of liquid water migration flux inside the snow layer and the data on the distribution of solid-ice phase transformation rate inside the snow layer, the mechanical stability assessment of the snow layer is performed to generate data on the distribution of pore water pressure field inside the snow layer and data on the distribution of effective stress field inside the snow layer. Based on the distribution data of pore water pressure field inside the snow layer and the distribution data of effective stress field inside the snow layer, disaster-causing process identification processing is performed to generate disaster-causing risk area distribution data and disaster-causing risk evolution path data for snow layer instability induced by rainfall infiltration. Based on the disaster risk area distribution data and the disaster risk evolution path data, a set of simulation results of the snow surface rain disaster process is generated. The set of simulation results of the snow surface rain disaster process includes snow layer instability trigger time window data and snow layer instability impact range data.

[0006] In conjunction with the second aspect of this application, a simulation system for snow surface rain-induced disaster processes based on a water-thermal coupling model is provided. The simulation system includes a machine-readable storage medium and a processor. The machine-readable storage medium stores machine-executable instructions. When the processor executes the machine-executable instructions, the simulation system implements the aforementioned simulation method for snow surface rain-induced disaster processes based on a water-thermal coupling model.

[0007] In conjunction with a third aspect of this application, a computer-readable storage medium is provided, wherein computer-executable instructions are stored therein, and when the computer-executable instructions are executed, the aforementioned method for simulating the snow surface rain-induced disaster process using a fused hydrothermal coupling model is implemented.

[0008] Combining any of the above aspects, a hydrothermal coupling model simultaneously generates the distribution of liquid water migration flux and solid-state ice phase change conversion rate within the snow layer. This allows the effects of heat conduction driven by temperature gradients during rainwater infiltration, endothermic melting of ice crystals, and exothermic freezing of liquid water to interact with the migration process of water in pores within the same computational system. This reflects the real-time impact of the dynamic adjustment of pore channels within the snow layer due to phase change on water migration paths and fluxes. Using both the liquid water migration flux distribution and the solid-state ice phase change conversion rate distribution as inputs for mechanical stability assessment, the generated pore water pressure field distribution data can reflect the contribution of liquid water accumulation and migration within the snow layer to local water pressure. The effective stress field distribution data can reflect the spatiotemporal changes in the snow layer's bearing capacity under the combined action of skeletal structural changes caused by ice phase change and water pressure. Together, these two data constitute a complete characterization of the mechanical state of the snow layer. Based on the joint distribution data of pore water pressure field and effective stress field, specific regions that gradually lose stability due to continuous increase in water pressure and continuous decrease in effective stress can be identified spatially. Simultaneously, the complete process of these regions from initial weakening to critical instability can be captured temporally. The resulting disaster risk evolution path data can fully present the temporal characteristics of disaster incubation and triggering. The final output data on snow layer instability triggering time window and instability impact range directly transforms the calculation results of the physical process into spatiotemporal information that can be used for disaster prediction, enabling the simulation results to describe both the timing of disaster occurrence and the spatial impact of the disaster. Attached Figure Description

[0009] Figure 1 A flowchart illustrating the snow surface rain-induced disaster simulation method using a fused hydrothermal coupling model provided in this application embodiment. Detailed Implementation

[0010] Figure 1 This document illustrates a flowchart of a method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model, as provided in an embodiment of this application. It should be understood that in other embodiments, the order of some steps in this method may be shared based on actual needs, or some steps may be omitted or maintained. The detailed components of this fusion hydrothermal coupling model-based snow surface rain-induced disaster process simulation method are as follows: Step S110: Obtain the initial snow layer state data set of the area to be simulated before the occurrence of the rainfall event. The initial snow layer state data set includes spatial distribution data of snow layer thickness, spatial distribution data of snow layer temperature, and spatial distribution data of snow layer liquid water content.

[0011] In this embodiment, the simulated area is a snow-covered slope along a mountain road. This slope experiences rainfall events during the snow season, and a snow layer exists on the slope surface before the rainfall events. The initial snow layer state data set is obtained by snow layer detection equipment deployed within the slope area. This snow layer detection equipment includes snow layer thickness detection rods, snow layer temperature chain sensors, and snow layer liquid water content sensors. The snow layer thickness detection rods are deployed at multiple horizontal locations within the slope area, with multiple thickness markers spaced vertically on each rod. The snow layer temperature chain sensors have multiple temperature detection nodes deployed along the snow layer depth direction, with the vertical spacing between adjacent temperature detection nodes being the same. The snow layer liquid water content sensors have multiple water content detection nodes deployed along the snow layer depth direction, with the vertical spacing between adjacent water content detection nodes being the same. The snow layer state data collected by the above snow layer detection equipment does not involve personal privacy data.

[0012] The spatial distribution data for snow layer thickness is a matrix of snow layer thickness values ​​under a spatial grid division of the area to be simulated. The area to be simulated is divided into regular spatial grids in the horizontal direction, and each spatial grid cell has a unique row and column index. The snow layer thickness values ​​are in meters. For example, the area to be simulated is divided into a spatial grid of 20 rows by 15 columns, for a total of 300 spatial grid cells. The snow layer thickness value of each spatial grid cell is read through thickness markers on a snow layer thickness probe. The spatial distribution data for snow layer temperature is the temperature value collected by temperature probe nodes layered along the snow layer depth direction within each spatial grid cell. The temperature value is in degrees Celsius. The spatial distribution data for snow layer liquid water content is the liquid water volume fraction value collected by liquid water content probe nodes layered along the snow layer depth direction within each spatial grid cell. The liquid water volume fraction value is in percentage.

[0013] Step S120: Perform snow surface rainfall infiltration process simulation processing on the initial snow layer state data set to generate data on the distribution of liquid phase water migration flux inside the snow layer and the distribution of solid phase ice phase transformation rate inside the snow layer during the simulated rainfall period.

[0014] In this embodiment, the simulation of snow surface rainfall infiltration uses an initial snow layer state data set as input. A hydrothermal coupled numerical simulation is used to calculate the migration of liquid water and the phase transformation of solid ice within the snow layer during rainfall infiltration. The simulated rainfall period is a continuous time interval from the start to the end of the rainfall event. Liquid water migration flux distribution data describes the amount of liquid water migrating in different directions at various spatial locations within the snow layer. Solid ice phase transformation rate distribution data describes the rates of solid ice to liquid water and vice versa at various spatial locations within the snow layer. The physical governing equations used in the hydrothermal coupled numerical simulation include the unsaturated porous media seepage control equation and the heat conduction control equation. The unsaturated porous media seepage control equation describes the migration process of liquid water in the pore space of the snow layer. The heat conduction control equation describes the evolution of the temperature field within the snow layer over time.

[0015] Step S121: Perform region splitting processing on the spatial distribution data of snow layer thickness in the initial snow layer state data set according to the spatial grid division method to obtain the sequence of spatial units inside the snow layer. Each spatial unit in the sequence of spatial units inside the snow layer has a layered structure attribute.

[0016] In this embodiment, the region splitting process treats each spatial grid cell in the snow layer thickness spatial distribution data as an independent spatial unit. Each spatial unit has horizontal position coordinates and a snow layer thickness value. Each spatial unit is divided into multiple layers along the snow layer depth direction, and the number of layers is determined by dividing the snow layer thickness value of the spatial unit by a preset layer thickness step size. For example, the preset layer thickness step size is 0.05 meters. When the snow layer thickness value of a spatial unit is 0.8 meters, the spatial unit is divided into 16 layers. When the snow layer thickness value of a spatial unit is 0.3 meters, the spatial unit is divided into 6 layers. All layers are numbered sequentially from the snow layer surface to the bottom of the snow layer, with the first layer located at the snow layer surface and the last layer located at the interface between the bottom of the snow layer and the ground. All spatial units are arranged in row and column index order, forming a sequence of spatial units within the snow layer. The layer structure attributes of each spatial unit include the total number of layers, the depth position of each layer, and the spatial coordinates of each layer.

[0017] Step S122: Based on the spatial distribution data of snow layer temperature and the spatial distribution data of snow layer liquid water content, perform initial state mapping processing on each spatial unit in the sequence of spatial units inside the snow layer to obtain the initial temperature state result and the initial liquid water content state result of each spatial unit at each layer position.

[0018] In this embodiment, for each spatial unit, the temperature values ​​of each temperature detection node within that spatial unit are read from the snow layer temperature spatial distribution data. The temperature values ​​of each temperature detection node are linearly interpolated along the snow layer depth direction to each layer location to obtain the initial temperature state result for each layer location. For example, in a certain spatial unit, temperature detection nodes are deployed at depths of 0 m, 0.25 m, 0.5 m, and 0.75 m, with temperature values ​​of -3.5 degrees Celsius, -2.8 degrees Celsius, -1.6 degrees Celsius, and -0.4 degrees Celsius, respectively. This spatial unit is divided into 16 layers, each with a depth interval of 0.05 m. The initial temperature state result for each layer location is obtained through linear interpolation. The liquid water volume fraction values ​​of each liquid water content detection node within the spatial unit are read from the snow layer liquid water content spatial distribution data. The liquid water volume fraction values ​​of each liquid water content detection node are linearly interpolated along the snow layer depth direction to each layer location to obtain the initial liquid water content state result for each layer location. The initial state results for liquid water content are expressed as a percentage.

[0019] Step S123: Perform rainfall infiltration boundary condition setting processing on the overall boundary region of the snow layer internal spatial unit sequence to obtain the top rainfall infiltration boundary parameters and the bottom water discharge boundary parameters of the snow layer internal spatial unit sequence.

[0020] In this embodiment, the top rainfall infiltration boundary parameter is the infiltration rate of rainwater into the snow surface during a rainfall event. This parameter is determined based on rainfall intensity data, and the infiltration rate equals the rainfall intensity multiplied by the snow surface infiltration coefficient. The snow surface infiltration coefficient is a dimensionless parameter, ranging from 0 to 1. The rainfall intensity data is obtained from hourly rainfall data recorded by meteorological observation stations along the mountainous road in the area to be simulated. For example, if the rainfall intensity is 5 mm / h and the snow surface infiltration coefficient is 0.7, the infiltration rate is 3.5 mm / h. The bottom water discharge boundary parameter is the water discharge rate at the interface between the bottom of the snow layer and the ground. This parameter is determined based on the ground permeability characteristics. When the ground is a highly permeable gravel layer, a larger value is used for the bottom water discharge rate. When the ground is a less permeable clay layer, a smaller value is used for the bottom water discharge rate.

[0021] Step S124: Perform liquid phase water migration driven calculation processing between layers in each spatial unit based on the top rainfall infiltration boundary parameters and the bottom water discharge boundary parameters to obtain liquid phase water migration direction data and liquid phase water migration amount data between layers in each spatial unit at each time node during the simulated rainfall period.

[0022] In this embodiment, the calculation of liquid phase water migration is performed using the unsaturated porous media seepage control equation. For each spatial unit, the driving force of liquid phase water migration between adjacent layers is calculated. This driving force is determined by both gravitational potential difference and matrix potential difference. The gravitational potential difference equals the depth difference between adjacent layers multiplied by the liquid phase water density multiplied by the gravitational acceleration. The matrix potential difference is determined based on the difference in liquid water content between adjacent layers and the snow layer matrix characteristic parameters. The direction of liquid phase water migration is determined by the sign of the driving force; a positive sign indicates downward migration, and a negative sign indicates upward migration. The amount of liquid phase water migration equals the driving force multiplied by the snow layer hydraulic conductivity coefficient and then multiplied by the time step. The snow layer hydraulic conductivity coefficient is an inherent hydraulic parameter of the snow layer, related to the snow layer density and porosity. The time step is set according to the simulation accuracy requirements; for example, a time step of 60 seconds is used. The entire time step covers the entire simulated rainfall period.

[0023] Step S125: Perform spatial aggregation processing on each of the liquid water migration data according to the spatial unit arrangement order to obtain the liquid water migration flux distribution data inside the snow layer. The liquid water migration flux distribution data inside the snow layer reflects the overall migration total characteristics of liquid water in the vertical direction and the migration distribution characteristics in the horizontal direction.

[0024] In this embodiment, spatial aggregation processing sums the liquid water migration data between all layers within the same spatial cell along the depth direction to obtain the total vertical migration amount for that spatial cell. The total vertical migration amounts of all spatial cells are then arranged by their row and column indices to obtain the vertical migration amount matrix in the liquid water migration flux distribution data. The number of rows and columns in the vertical migration amount matrix is ​​consistent with the number of rows and columns in the spatial grid of the region to be simulated. Simultaneously, the differences in liquid water migration amounts at the same layer position between adjacent spatial cells are arranged by their row and column indices to obtain the horizontal migration distribution matrix in the liquid water migration flux distribution data. The horizontal migration distribution matrix reflects the distribution characteristics of liquid water migrating horizontally from spatial cells with high liquid water content to those with low liquid water content.

[0025] Step S126: Based on the liquid water migration direction data and the liquid water migration amount data, combined with the initial temperature state results of each spatial unit's stratification location, perform solid-phase ice phase transition determination processing to obtain solid-phase ice phase transition type data and solid-phase ice phase transition degree data at each stratification location of each spatial unit.

[0026] In this embodiment, the solid-state ice phase transition determination process compares the initial temperature state at each layer location of each spatial unit with the freezing point temperature. The freezing point temperature is set to 0 degrees Celsius. When the initial temperature state is higher than the freezing point temperature and solid ice exists at the layer location, the solid-state ice phase transition type is determined to be melting. When the initial temperature state is lower than the freezing point temperature and liquid water exists at the layer location, the solid-state ice phase transition type is determined to be freezing. The solid-state ice phase transition degree data is determined based on the magnitude of temperature deviation from the freezing point temperature and the liquid water migration data. The solid-state ice phase transition degree corresponding to the melting type is equal to the magnitude of temperature deviation from the freezing point temperature multiplied by the melting rate coefficient. The solid-state ice phase transition degree corresponding to the freezing type is equal to the magnitude of temperature deviation from the freezing point temperature multiplied by the freezing rate coefficient. Both the melting rate coefficient and the freezing rate coefficient are physical constants, with the melting rate coefficient being greater than the freezing rate coefficient.

[0027] Step S127: Spatial correlation and aggregation processing is performed on the solid ice phase change type data and the solid ice phase change degree data to obtain the solid ice phase change conversion rate distribution data inside the snow layer. The solid ice phase change conversion rate distribution data inside the snow layer reflects the rate attribute of solid ice to liquid water conversion and the rate attribute of liquid water to solid ice conversion at each layer position of each spatial unit.

[0028] In this embodiment, the spatial association aggregation process combines the solid-ice phase transition type data and solid-ice phase transition degree data for each spatial unit and each layer position into a binary data structure. The first element of the binary data structure is the solid-ice phase transition type identifier, and the second element is the solid-ice phase transition degree scalar. All binary data structures are arranged according to the row and column index of the spatial unit and the layer position, forming the solid-ice phase transition rate distribution data.

[0029] Step S128: Perform layered temperature state update processing based on the liquid water migration flux distribution data in the snow layer and the solid ice phase change degree data in the solid ice phase change rate distribution data in the snow layer to obtain the updated temperature state data of each layer position in each spatial unit.

[0030] In this embodiment, the layered temperature state update process performs addition and subtraction operations on the product of the initial temperature state result of each layer location in each spatial unit and the product of the liquid water migration amount data multiplied by the specific heat capacity of liquid water and the temperature difference. Simultaneously, it performs addition and subtraction operations on the product of the solid ice phase change degree data multiplied by the latent heat of solid ice phase change, to obtain the updated temperature state data. Both the specific heat capacity of liquid water and the latent heat of solid ice phase change are physical constants. The temperature difference is determined based on the heat exchange caused by the liquid water migration. When liquid water migrates downwards and the upper layer temperature is higher than the lower layer temperature, the temperature difference is negative, indicating that the downward migration of upper liquid water causes the lower layer temperature to rise. When liquid water migrates upwards and the lower layer temperature is lower than the upper layer temperature, the temperature difference is positive, indicating that the upward migration of lower liquid water causes the upper layer temperature to decrease.

[0031] Step S129: Use the updated temperature state data as the initial temperature state input data for the next time node, return to execute the rainfall infiltration boundary condition setting process and the liquid phase water migration driving calculation process until the preset simulation termination time is reached, and obtain the liquid phase water migration flux distribution data and the solid phase ice phase change conversion rate distribution data of the snow layer corresponding to each time node during the simulated rainfall period.

[0032] In this embodiment, the simulation termination time is the end time of the rainfall event. The continuous time interval from the start time to the end time of the rainfall event is divided into multiple time steps. Within each time step, a complete process of setting rainfall infiltration boundary conditions, liquid phase water migration-driven calculation, solid phase ice phase change determination, and stratification temperature state update is performed. After all time steps are completed, the liquid phase water migration flux distribution data and solid phase ice phase change conversion rate distribution data corresponding to each time node are stored in chronological order. The storage structure is a key-value pair data structure with time nodes as keys and liquid phase water migration flux distribution data and solid phase ice phase change conversion rate distribution data as values.

[0033] Step S130: Perform a snow layer mechanical stability assessment based on the distribution data of liquid water migration flux inside the snow layer and the distribution data of solid-ice phase transformation rate inside the snow layer to generate pore water pressure field distribution data and effective stress field distribution data inside the snow layer.

[0034] In this embodiment, the snow layer mechanical stability assessment process transforms the hydrothermal coupling simulation results into a description of the internal mechanical state of the snow layer. Pore water pressure field distribution data describes the pressure state of liquid water in the pore space within the snow layer. Effective stress field distribution data describes the actual stress state borne by the solid skeleton of the snow layer. The snow layer mechanical stability assessment process is performed based on the effective stress principle and pore water pressure theory.

[0035] Step S131: Extract the liquid water migration direction data and liquid water migration amount data of each layer position in each spatial unit from the liquid water migration flux distribution data inside the snow layer, and determine the flow thrust distribution data of liquid water acting on the solid skeleton at each layer position based on the liquid water migration direction data and the liquid water migration amount data.

[0036] In this embodiment, the flow thrust distribution data is calculated by multiplying the liquid water migration amount data by the liquid water density, and then multiplying by the unit direction vector corresponding to the liquid water migration direction data to obtain the flow thrust vector. The magnitude of the flow thrust vector is the liquid water migration amount multiplied by the liquid water density, and the direction of the flow thrust vector is consistent with the liquid water migration direction. The liquid water density is a physical constant with a value of 1000 kg / m³. All flow thrust vectors constitute the flow thrust distribution data. Each flow thrust vector in the flow thrust distribution data corresponds to a layer position in a spatial unit.

[0037] Step S132: Based on the solid-ice phase change type data and solid-ice phase change degree data in the solid-ice phase change rate distribution data inside the snow layer, perform solid volume ratio update processing at each layer location to obtain solid volume ratio change sequence data at each layer location. The solid volume ratio change sequence data reflects the solid volume reduction caused by the conversion of solid ice to liquid water and the solid volume increase caused by the conversion of liquid water to solid ice.

[0038] In this embodiment, for each stratification location, the initial solid volume percentage is determined by the initial density of the snow layer and the ice density. The initial solid volume percentage equals the initial density of the snow layer divided by the ice density. The initial density of the snow layer is estimated using an empirical formula based on the spatial distribution data of snow layer temperature and thickness. The ice density is a physical constant, with a value of 917 kg / m³. When the solid-ice phase transition type is melting, the decrease in the solid volume percentage at that stratification location equals the solid-ice phase transition degree data multiplied by the time step. When the solid-ice phase transition type is freezing, the increase in the solid volume percentage at that stratification location equals the solid-ice phase transition degree data multiplied by the time step. The changes in the solid volume percentage over all time steps constitute the solid volume percentage change sequence data.

[0039] Step S133: Based on the solid volume ratio change sequence data of each layer location, perform pore space ratio update processing for each layer location to obtain pore space ratio update data for each layer location. The pore space ratio update data reflects the dynamic change distribution characteristics of the pore space inside the snow layer as a result of the solid-ice phase transformation process.

[0040] In this embodiment, the pore space percentage at each layer location is equal to 1 minus the solid phase volume percentage. The solid phase volume percentage at each time step in the solid phase volume percentage change sequence data is substituted into the pore space percentage calculation formula to obtain the pore space percentage for the corresponding time step. The pore space percentages for all time steps constitute the updated pore space percentage data. The pore space percentage is a dimensionless parameter, with a value between 0 and 1.

[0041] Step S134: Perform a one-to-one spatial unit position correlation processing on the distribution data of the flow thrust of liquid water acting on the solid skeleton at each layer position and the updated data of the pore space ratio at each layer position to obtain the pore water pressure field distribution data of each layer position of each spatial unit. The pore water pressure field distribution data reflects the pressure magnitude distribution and pressure direction distribution of liquid water in the pore space.

[0042] In this embodiment, for each spatial unit and each layer location, the flow thrust vector at that layer location is divided by the proportion of pore space at that layer location to obtain the pore water pressure vector. The magnitude of the pore water pressure vector is the magnitude of the pore water pressure, and the direction of the pore water pressure vector is the direction of the pore water pressure. The unit of measurement for the magnitude of the pore water pressure is Pascal. All pore water pressure vectors constitute the pore water pressure field distribution data. The pore water pressure field distribution data is organized and stored using the spatial unit row and column indices and the layer location as indices.

[0043] Step S135: Based on the solid volume ratio change sequence data of each layer position and the initial gravity load data of each layer position of each spatial unit, perform the total stress field generation process inside the snow layer to obtain the total stress distribution data of each layer position of each spatial unit. The total stress distribution data includes the vertical stress component caused by the gravity of the snow layer itself and the boundary stress component applied from the outside.

[0044] In this embodiment, for each spatial unit and each layer location, the vertical stress component is equal to the sum of the density of all snow layers above that layer location multiplied by the gravitational acceleration and the layer thickness. The snow layer density is calculated based on the solid volume ratio and ice density at that layer location. The gravitational acceleration is a physical constant, with a value of 9.8 m / s². The layer thickness is a preset layer thickness step. The boundary stress component is the external stress applied at the boundary of the region to be simulated. The external stress originates from the lateral compression of the surrounding snow mass of the region to be simulated. The total stress is equal to the sum of the vertical stress component and the boundary stress component. The unit of measurement for the total stress is Pascal.

[0045] Step S136: Perform effective stress solution processing on the total stress distribution data and the pore water pressure field distribution data of each layer position of each spatial unit to obtain the effective stress distribution data of each layer position of each spatial unit. The effective stress distribution data reflects the actual stress state and stress transmission path characteristics of the solid skeleton.

[0046] In this embodiment, the effective stress equals the total stress minus the pore water pressure. For each spatial unit and each layer location, the effective stress value is obtained by subtracting the pore water pressure value from the pore water pressure field distribution data in the total stress distribution data. All effective stress values ​​constitute the effective stress distribution data. The unit of measurement for effective stress is Pascal. The effective stress distribution data is organized and stored using the spatial unit row and column indices and the layer location as indices.

[0047] Step S140: Based on the distribution data of pore water pressure field inside the snow layer and the distribution data of effective stress field inside the snow layer, perform disaster-causing process identification processing to generate disaster-causing risk area distribution data and disaster-causing risk evolution path data for snow layer instability induced by rainfall infiltration.

[0048] In this embodiment, the disaster-causing process identification process identifies spatial regions that may lead to snow layer instability from pore water pressure field distribution data and effective stress field distribution data, and tracks the evolution path of these spatial regions during simulated rainfall periods. The disaster-causing risk region distribution data consists of a set of spatial units identified as having snow layer instability risk. The disaster-causing risk evolution path data is the spatial trajectory of the disaster-causing risk regions expanding within the simulated area as the simulated rainfall period progresses.

[0049] Step S141: Extract pore water pressure numerical data of each layer position of each spatial unit from the pore water pressure field distribution data inside the snow layer, and perform pore water pressure exceeding limit area determination processing based on the pore water pressure numerical data to obtain the set data of pore water pressure exceeding limit spatial units and the exceeding limit degree data of each exceeding limit spatial unit.

[0050] In this embodiment, for each spatial unit and each layer location, the pore water pressure value is compared with a preset pore water pressure exceeding threshold. The preset pore water pressure exceeding threshold is determined based on the shear strength parameters of the snow layer solid skeleton. When the pore water pressure value is not lower than the preset pore water pressure exceeding threshold, the spatial unit and that layer location are determined to be a pore water pressure exceeding threshold location. Spatial units containing pore water pressure exceeding threshold locations are marked as pore water pressure exceeding threshold spatial units. The degree of exceeding the threshold is equal to the pore water pressure value minus the preset pore water pressure exceeding threshold. All pore water pressure exceeding threshold spatial units constitute the pore water pressure exceeding threshold spatial unit set data.

[0051] Step S142: Extract the effective stress numerical data of each layer position of each spatial unit from the effective stress field distribution data inside the snow layer, and perform effective stress mutation region determination processing based on the effective stress numerical data to obtain the set data of effective stress mutation spatial units and the mutation degree data of each mutation spatial unit.

[0052] In this embodiment, for each spatial unit at each layer location, the effective stress value is calculated to be different from the effective stress value at the same layer location in adjacent spatial units. When the absolute value of the difference is not lower than a preset effective stress mutation threshold, the layer location of the spatial unit is determined to be an effective stress mutation location. The spatial unit containing the effective stress mutation location is marked as an effective stress mutation spatial unit. The mutation degree is equal to the absolute value of the difference minus the preset effective stress mutation threshold. All effective stress mutation spatial units constitute the effective stress mutation spatial unit set data.

[0053] Step S143: Perform spatial position cross-matching processing on the set data of spatial units with pore water pressure exceeding the limit and the set data of spatial units with effective stress change to obtain the data of spatial units with pore water pressure exceeding the limit and effective stress change. The data of spatial units with overlapping regions is used as the distribution data of disaster-causing risk areas for snow layer instability induced by rainfall infiltration.

[0054] In this embodiment, the spatial location cross-matching process checks whether each spatial cell exists simultaneously in both the pore water pressure exceeding limit spatial cell set data and the effective stress mutation spatial cell set data. If both exist simultaneously, the spatial cell is marked as an overlapping spatial cell. All overlapping spatial cells constitute the disaster risk area distribution data. The disaster risk area distribution data is a list structure, where each element in the list is the row and column index of an overlapping spatial cell.

[0055] Step S144: Perform disaster risk spatial connectivity analysis processing based on the spatial position arrangement order of each overlapping spatial unit in the spatial unit overlapping area data to obtain the spatial connectivity path set data between disaster risk areas. The spatial connectivity path set data includes spatial connection direction data and spatial connection order data between adjacent disaster risk areas.

[0056] In this embodiment, for all overlapping spatial units, it is checked whether there is a shared boundary between each pair of adjacent overlapping spatial units. When a shared boundary exists, it is determined that there is a spatial connectivity relationship between the two adjacent overlapping spatial units. The spatial connection direction is determined based on the difference between the row and column indices of the adjacent overlapping spatial units. The sign of the row and column index difference corresponds to the spatial connection direction. The spatial connection order is determined based on the order in which adjacent overlapping spatial units are marked as overlapping spatial units during the simulated rainfall period. All spatial connectivity relationships constitute a set of spatial connectivity path data.

[0057] Step S145: Perform disaster risk evolution trend tracking processing based on the spatial connection direction data and spatial connection sequence data in the spatial connectivity path set data to obtain the spatial diffusion direction sequence data and spatial diffusion rate sequence data of the disaster risk in the area to be simulated.

[0058] In this embodiment, the spatial diffusion direction sequence data is the direction sequence of the expansion of the disaster-causing risk area along the spatial connection direction within the area to be simulated, step by step. The spatial diffusion rate sequence data is the number of newly added overlapping spatial units in each time step divided by the time step. The spatial diffusion rate is measured in units per second.

[0059] Step S146: The disaster risk spatial connectivity path set data, the spatial diffusion direction sequence data, and the spatial diffusion rate sequence data are comprehensively merged to generate the disaster risk evolution path data. The disaster risk evolution path data reflects the spatial trajectory characteristics and expansion speed change characteristics of the disaster risk as it gradually expands from the initial triggering area to the surrounding area.

[0060] In this embodiment, the disaster risk evolution path data is a structured data object containing a set of spatial connectivity paths, a sequence of spatial diffusion directions, and a sequence of spatial diffusion rates. The disaster risk evolution path data is organized and stored using time nodes as indexes.

[0061] Step S150: Generate a set of simulation results for the snow surface rain disaster process based on the disaster risk area distribution data and the disaster risk evolution path data. The set of simulation results for the snow surface rain disaster process includes snow layer instability trigger time window data and snow layer instability impact range data.

[0062] In this embodiment, the simulation results dataset for the snow surface rain-induced disaster process is a summary of the results of the snow layer instability process during the simulated rainfall period. The snow layer instability trigger time window data represents the time range from the start to the end of snow layer instability. The snow layer instability impact range data represents the spatial coverage of the instability area within the area to be simulated.

[0063] Step S151: Extract the spatial unit identifier data corresponding to each disaster risk area and the disaster risk start time point data corresponding to the spatial unit identifier data from the disaster risk area distribution data. Perform trigger time series sorting processing on each disaster risk area according to the disaster risk start time point data to obtain the time dimension distribution sequence data of each disaster risk area.

[0064] In this embodiment, the disaster risk initiation time point data is the time node when the spatial unit is first marked as an overlapping spatial unit. All disaster risk areas are arranged in ascending order according to their disaster risk initiation time points to obtain the time-dimensional distribution sequence data. The time-dimensional distribution sequence data is a list structure, where each element is a data pair containing a spatial unit identifier and the disaster risk initiation time point.

[0065] Step S152: Perform snow layer instability triggering time window range delineation processing based on the time dimension distribution sequence data of each disaster risk area to obtain snow layer instability triggering time window data, which includes the earliest triggering time point data and the latest triggering time point data of snow layer instability.

[0066] In this embodiment, the earliest trigger time for snow instability is the disaster risk initiation time of the first disaster-causing risk area in the time-dimensional distribution sequence data. The latest trigger time for snow instability is the disaster risk initiation time of the last disaster-causing risk area in the time-dimensional distribution sequence data. The snow instability trigger time window data is a binary data structure containing the earliest and latest trigger times.

[0067] Step S153: Extract the spatial diffusion direction sequence data and spatial diffusion rate sequence data of disaster risk from the disaster risk evolution path data, and perform snow layer instability influence range boundary estimation processing based on the disaster risk spatial diffusion direction sequence data and the disaster risk spatial diffusion rate sequence data to obtain the snow layer instability influence range boundary coordinate sequence data of each time node during the simulated rainfall period.

[0068] In this embodiment, for each time node, the spatial coordinates of all overlapping spatial units at that time node are collected, and the spatial coordinates of overlapping spatial units located at the edge are connected sequentially to form the boundary coordinate sequence data of the snow layer instability influence range. An overlapping spatial unit located at the edge is an overlapping spatial unit whose at least one adjacent spatial unit is not an overlapping spatial unit. The boundary coordinate sequence data of the snow layer instability influence range is an ordered coordinate list.

[0069] Step S154: Based on the earliest and latest trigger time points in the snow instability trigger time window data, and combined with the boundary coordinate sequence data of the snow instability influence range, perform dynamic mapping processing of the snow instability influence range over time to obtain the snow instability influence range data. The snow instability influence range data reflects the expansion coverage area change characteristics of the instability area in the time dimension and the expansion direction change characteristics in the spatial dimension.

[0070] In this embodiment, the dynamic mapping process binds the boundary coordinate sequence data of the snow layer instability influence range at each time node to that time node, thus obtaining the snow layer instability influence range data. The snow layer instability influence range data is a key-value pair data structure with the time node as the key and the boundary coordinate sequence as the value.

[0071] Step S210: After obtaining the initial snow layer state data set of the area to be simulated before the occurrence of the rainfall event, and before performing the snow surface rainfall infiltration process simulation processing on the initial snow layer state data set, perform the validity determination processing of the initial snow layer state data set to obtain a valid initial snow layer state data set.

[0072] In this embodiment, the validity determination process performs quality review and correction on three types of data in the initial snow layer state dataset. These three types of data are spatial distribution data of snow layer thickness, spatial distribution data of snow layer temperature, and spatial distribution data of snow layer liquid water content. The quality review includes spatial continuity checks, temperature gradient rationality checks, and physical boundary checks of liquid water content. The quality correction includes missing data completion, outlier data correction, temperature interpolation correction, and adjustment for out-of-limit values.

[0073] Step S211: Perform a spatial continuity check on the spatial distribution data of snow layer thickness in the initial snow layer state data set to obtain the continuity index data of the spatial distribution data of snow layer thickness at each spatial location.

[0074] In this embodiment, the spatial continuity check process calculates the absolute value of the difference between the snow layer thickness value of each spatial grid cell and the snow layer thickness values ​​of its four adjacent spatial grid cells (upper, lower, left, and right). The arithmetic mean of the four absolute differences is used as the continuity index data for that spatial grid cell. The continuity index data is in meters. The smaller the continuity index data, the more continuous the snow layer thickness transition between the spatial grid cell and its adjacent spatial grid cells. For example, if the snow layer thickness value of a certain spatial grid cell is 0.8 meters, and the snow layer thickness values ​​of its four adjacent spatial grid cells (upper, lower, left, and right) are 0.75 meters, 0.82 meters, 0.78 meters, and 0.85 meters, respectively, and the four absolute differences are 0.05 meters, 0.02 meters, 0.02 meters, and 0.05 meters, respectively, the continuity index data is 0.035 meters.

[0075] Step S212: Identify the missing data region and abnormal data region in the spatial distribution data of snow layer thickness based on the continuity index data, and obtain the missing region identification data and abnormal region identification data.

[0076] In this embodiment, the missing data region is the area where spatial grid cells with snow layer thickness values ​​of null or zero are located. The abnormal data region is the area where spatial grid cells with continuity index data exceeding a preset continuity anomaly threshold are located. The preset continuity anomaly threshold is determined based on the terrain undulation characteristics of the area to be simulated and the acquisition accuracy of the spatial distribution data of snow layer thickness. The row and column indices of the spatial grid cells in the missing data region are collected as missing region identifier data. The row and column indices of the spatial grid cells in the abnormal data region are collected as abnormal region identifier data.

[0077] Step S213: Perform missing data region completion processing and abnormal data region correction processing based on the missing region identifier data and the abnormal region identifier data to obtain the corrected spatial distribution data of snow layer thickness.

[0078] In this embodiment, the missing data region completion process takes the arithmetic mean of the snow thickness values ​​of the four adjacent spatial grid cells (up, down, left, and right) of each missing data region as the completed snow thickness value for that region. The abnormal data region correction process takes the arithmetic mean of the snow thickness values ​​of the four adjacent spatial grid cells (up, down, left, and right) of each abnormal data region and replaces the original snow thickness value for that region. After all completion and correction are completed, the corrected spatial distribution data of snow thickness is obtained.

[0079] Step S214: Perform a temperature gradient rationality check on the spatial distribution data of snow layer temperature in the initial snow layer state data set to obtain temperature gradient rationality index data at each spatial location.

[0080] In this embodiment, the temperature gradient rationality check process calculates the temperature difference between adjacent layers along the snow depth direction within each spatial grid cell. The temperature gradient rationality index is equal to the temperature difference divided by the layer thickness step size. The temperature gradient rationality index is measured in degrees Celsius per meter. When the absolute value of the temperature gradient rationality index exceeds the preset upper limit of the temperature gradient rationality, the layer location is marked as an unreasonable temperature gradient location. The preset upper limit of the temperature gradient rationality is determined based on the thermal conductivity physical properties of the snow layer.

[0081] Step S215: Based on the temperature gradient rationality index data, identify the temperature abrupt change location data and the temperature reverse change location data in the snow layer temperature spatial distribution data to obtain the temperature anomaly location set data.

[0082] In this embodiment, the location of a temperature abrupt change is the stratum where the absolute value of the temperature gradient rationality index data exceeds the preset upper limit of the temperature gradient rationality. The location of a temperature reversal is the stratum where the direction of the temperature difference between adjacent strata is opposite to the overall temperature change direction. The overall temperature change direction is the direction of the temperature change trend from the snow surface to the bottom of the snow layer. All temperature abrupt change locations and temperature reversal change locations constitute the temperature anomaly location set data.

[0083] Step S216: Perform interpolation correction processing on the spatial distribution data of snow layer temperature based on the set of temperature anomaly locations to obtain the corrected spatial distribution data of snow layer temperature.

[0084] In this embodiment, for each anomalous location in the temperature anomaly location dataset, the arithmetic mean of the temperature values ​​of the two adjacent non-anomaly layers above and below the anomalous location is taken as the corrected temperature value for that anomalous location. When the anomalous location is the first layer of the snow surface, the arithmetic mean of the temperature values ​​of the two adjacent non-anomaly layers below the anomalous location is taken as the corrected temperature value. When the anomalous location is the last layer at the bottom of the snow layer, the arithmetic mean of the temperature values ​​of the two adjacent non-anomaly layers above the anomalous location is taken as the corrected temperature value. After all anomalous locations are corrected, the corrected spatial distribution data of snow layer temperature is obtained.

[0085] Step S217: Perform physical boundary review processing on the spatial distribution data of liquid water content in the snow layer in the initial snow layer state data set to obtain the physical boundary compliance index data of liquid water content at each spatial location.

[0086] In this embodiment, the physical boundary check process for liquid water content performs boundary checks on the liquid water volume fraction value at each layer location within each spatial grid cell. The physical boundary for the liquid water volume fraction value is from 0 to 100. When the liquid water volume fraction value is below 0, the physical boundary compliance index data is marked as below the lower limit. When the liquid water volume fraction value is above 100, the physical boundary compliance index data is marked as above the upper limit. When the liquid water volume fraction value is between 0 and 100, the physical boundary compliance index data is marked as normal.

[0087] Step S218: Adjust the out-of-limit data in the spatial distribution data of liquid water content in the snow layer according to the physical boundary compliance index data of liquid water content, and obtain the adjusted spatial distribution data of liquid water content in the snow layer.

[0088] In this embodiment, liquid water volume fraction values ​​below the lower limit are adjusted to 0. Liquid water volume fraction values ​​above the upper limit are adjusted to 100. After all excess values ​​are adjusted, the adjusted spatial distribution data of snow liquid water content is obtained.

[0089] Step S219: The corrected snow layer thickness spatial distribution data, the corrected snow layer temperature spatial distribution data, and the adjusted snow layer liquid water content spatial distribution data are aggregated and merged to generate the effective initial snow layer state data set.

[0090] In this embodiment, the effective initial snow layer state data set is a data object containing corrected spatial distribution data of snow layer thickness, corrected spatial distribution data of snow layer temperature, and adjusted spatial distribution data of snow layer liquid water content. This effective initial snow layer state data set replaces the original initial snow layer state data set as input data for subsequent simulation processing of snow surface precipitation infiltration.

[0091] Step S310: Before performing the snow layer mechanical stability assessment processing based on the distribution data of liquid water migration flux inside the snow layer and the distribution data of solid-ice phase change conversion rate inside the snow layer, perform the coupling coordination verification processing between liquid water migration flux and solid-ice phase change conversion rate to obtain the coupling coordination verification result.

[0092] In this embodiment, the coupling consistency verification process checks the temporal and spatial consistency between the liquid phase water migration flux distribution data and the solid phase ice phase change conversion rate distribution data. Temporal consistency is checked using a temporal synchronization index, and spatial consistency is checked using a spatial consistency score.

[0093] Step S311: Extract the total liquid water migration data for each time node during the simulated rainfall period from the liquid water migration flux distribution data inside the snow layer, and extract the total solid-ice phase change conversion data for the corresponding time node from the solid-ice phase change conversion rate distribution data inside the snow layer.

[0094] In this embodiment, the total liquid water migration data is the sum of liquid water migration data at all depths across all spatial units at each time point. The total liquid water migration data is measured in cubic meters. The total solid-state ice phase change data is the sum of the solid-state ice phase change degree data at all depths across all spatial units at each time point. The total solid-state ice phase change data is also measured in cubic meters.

[0095] Step S312: Calculate the time synchronization index data between the total liquid phase water migration data and the total solid phase ice phase transformation data based on the total liquid phase water migration data and the total solid phase ice phase transformation data. The time synchronization index data reflects the relative offset between the total liquid phase water migration data and the total solid phase ice phase transformation data in the time dimension.

[0096] In this embodiment, the time synchronization index data is calculated by taking the Pearson correlation coefficient between the total liquid water migration data and the total solid-phase ice phase transformation data for each time point. The closer the Pearson correlation coefficient is to 1, the better the time synchronization. The time synchronization index data is equal to the Pearson correlation coefficient.

[0097] Step S313: Based on the liquid water migration direction data of each spatial unit in the liquid water migration flux distribution data and the solid ice phase change type data of each spatial unit in the solid ice phase change rate distribution data, perform spatial correspondence analysis to obtain spatial consistency score data between the liquid water migration direction and the solid ice phase change type of each spatial unit.

[0098] In this embodiment, for each spatial unit, when the liquid water migration direction data is downward and the solid ice phase change type data is melting, the spatial consistency score is 1. When the liquid water migration direction data is upward and the solid ice phase change type data is freezing, the spatial consistency score is 1. In all other combinations, the spatial consistency score is 0. A spatial consistency score of 1 indicates that the liquid water migration direction and the solid ice phase change type are spatially consistent. A spatial consistency score of 0 indicates that the liquid water migration direction and the solid ice phase change type are spatially inconsistent.

[0099] Step S314: Based on the spatial consistency score data, identify the set of spatial units where spatial mismatch occurs between the liquid phase water migration direction and the solid phase ice phase change type, and obtain spatial mismatch unit identification data and mismatch type data corresponding to each mismatch unit.

[0100] In this embodiment, spatial units with a spatial consistency score of 0 are marked as spatial mismatch units. The mismatch type data can be either a downward migration and freezing mismatch or an upward migration and melting mismatch. A downward migration and freezing mismatch indicates that the spatial unit is freezing while the liquid water is migrating downwards. An upward migration and melting mismatch indicates that the spatial unit is melting while the liquid water is migrating upwards.

[0101] Step S315: Based on the time synchronization index data and the spatial mismatch unit identification data, perform a comprehensive judgment process on the coupling relationship between liquid water migration flux and solid ice phase change conversion rate to obtain the coupling coordination verification result. The coupling coordination verification result includes coupling coordination qualified region data and coupling coordination needing adjustment region data.

[0102] In this embodiment, when the time synchronization index data is not lower than a preset synchronization threshold and the number of spatial mismatch units does not exceed a preset mismatch number threshold, all regions are determined to be regions with qualified coupling coordination. When the time synchronization index data is lower than the preset synchronization threshold or the number of spatial mismatch units exceeds the preset mismatch number threshold, the region where the spatial mismatch unit is located is determined to be a region where coupling coordination needs to be adjusted. The coupling coordination verification results are used to guide the subsequent snow layer mechanical stability assessment process in performing correction processing on the regions where coupling coordination needs to be adjusted.

[0103] For example, in step S410: after performing the disaster-causing process identification processing based on the distribution data of the pore water pressure field inside the snow layer and the distribution data of the effective stress field inside the snow layer, and before generating the snow surface rain disaster-causing process simulation result data set based on the distribution data of the disaster-causing risk area and the disaster-causing risk evolution path data, the spatial correlation consistency verification processing between the disaster-causing risk area distribution data and the disaster-causing risk evolution path data is performed to obtain the spatial correlation consistency verification result.

[0104] In this embodiment, the spatial correlation consistency verification process examines the spatial location matching degree and temporal consistency between the disaster risk area distribution data and the disaster risk evolution path data. The spatial location matching degree is measured by spatial overlap data. Temporal consistency is measured by temporal consistency data.

[0105] Step S411: Extract the spatial coordinate data of each disaster-causing risk area and the trigger time point data of each disaster-causing risk area from the disaster-causing risk area distribution data; extract the spatial coordinate data of each path node in the disaster-causing risk evolution path and the time point data corresponding to each path node from the disaster-causing risk evolution path data.

[0106] In this embodiment, the spatial coordinates of the disaster-causing risk area are the row and column indices of the overlapping spatial units. The spatial coordinates of the path nodes are the row and column indices of each connected node in the spatial connectivity path. The trigger time point data is the time node when the overlapping spatial unit is first marked as an overlapping spatial unit. The time point data corresponding to the path node is the time node when the path node is included in the spatial connectivity path.

[0107] Step S412: Perform spatial coordinate matching processing on the spatial location coordinate data of each disaster-causing risk area and the spatial location coordinate data of each path node to obtain the spatial overlap data between the disaster-causing risk area and the evolution path node.

[0108] In this embodiment, the spatial overlap data is calculated by dividing the number of elements in the intersection of the spatial coordinate data set of the disaster-causing risk area and the spatial coordinate data set of the path nodes by the total number of elements in the spatial coordinate data set of the disaster-causing risk area. The spatial overlap data is a dimensionless parameter, ranging from 0 to 1. The closer the spatial overlap data is to 1, the higher the degree of spatial matching between the disaster-causing risk area and the evolution path nodes.

[0109] Step S413: Perform time series alignment processing on the trigger time point data of each disaster-causing risk area and the time point data corresponding to each path node to obtain time consistency data between the trigger time point of the disaster-causing risk area and the time point of the evolution path node.

[0110] In this embodiment, the time consistency data is the number of path nodes whose absolute value of the difference between the trigger time data and the path node time data does not exceed the preset time alignment tolerance, divided by the total number of path nodes. The time consistency data is a dimensionless parameter, with a value between 0 and 1. The closer the time consistency data is to 1, the higher the consistency between the trigger time of the disaster-causing risk area and the time of the evolution path node.

[0111] Step S414: Based on the spatial overlap data and the temporal consistency data, perform a comprehensive spatial correlation scoring process between the disaster risk area distribution data and the disaster risk evolution path data to obtain the spatial correlation consistency verification result. The spatial correlation consistency verification result includes spatial correlation fully matched area data and spatial correlation deviation area data.

[0112] In this embodiment, the spatial correlation comprehensive score is equal to the spatial overlap data multiplied by the temporal consistency data. The spatial correlation comprehensive score is a dimensionless parameter, with a value between 0 and 1. When the spatial correlation comprehensive score is not lower than the preset comprehensive matching threshold, the region is marked as a spatially correlated fully matched region. When the spatial correlation comprehensive score is lower than the preset comprehensive matching threshold, the region is marked as a region with spatial correlation deviation.

[0113] Step S415: Based on the spatial correlation deviation area data, perform spatial location correction processing between the disaster risk area distribution data and the disaster risk evolution path data to obtain the corrected disaster risk area distribution data and the corrected disaster risk evolution path data.

[0114] In this embodiment, for each deviation location in a spatially correlated deviation area, the spatial coordinates of that deviation location are replaced with the spatial coordinates of the corresponding path node in the disaster risk evolution path data. After all deviation locations are corrected, corrected disaster risk area distribution data and corrected disaster risk evolution path data are obtained. The corrected data is used for the subsequent generation and processing of the snow surface rain disaster process simulation result dataset.

[0115] Step S510: After generating the simulation result data set of snow surface rain disaster-causing process based on the disaster risk area distribution data and the disaster risk evolution path data, perform reliability evaluation processing on the simulation result data set of snow surface rain disaster-causing process to obtain reliability evaluation results.

[0116] In this embodiment, the reliability evaluation process performs a rationality review on the trigger time window data and impact range data in the simulation results dataset of the snow surface rain-induced disaster process. The rationality review includes a temporal rationality evaluation and a boundary shape complexity analysis.

[0117] Step S511: Extract the earliest and latest snow instability trigger time points from the snow instability trigger time window data in the snow surface rain-induced disaster simulation result data set, and extract the boundary coordinate sequence data of the instability influence range at each time node from the snow instability influence range data.

[0118] Step S512: Based on the earliest and latest trigger time data of snow layer instability, perform the trigger time window length calculation to obtain the total length data of the snow layer instability trigger time window. Based on the boundary coordinate sequence data of the instability influence range at each time node, perform the instability influence range total area evolution trend calculation to obtain the cumulative curve data of the instability influence range expansion over time.

[0119] In this embodiment, the total length of the snow layer instability triggering time window is equal to the latest triggering time of snow layer instability minus the earliest triggering time of snow layer instability. The total length of the snow layer instability triggering time window is measured in seconds. The calculation and processing of the total area evolution trend of the instability impact range involves calculating the area value enclosed by the boundary coordinate sequence data of the instability impact range at each time node. All area values ​​are arranged in chronological order to form the cumulative curve data of the expansion of the instability impact range over time. The area values ​​are measured in square meters.

[0120] Step S513: Based on the total length of the snow layer instability triggering time window and the cumulative curve data of the instability influence range over time, perform time rationality evaluation processing on the simulation result data set of the snow surface rain disaster process to obtain time rationality evaluation index data.

[0121] In this embodiment, the time rationality evaluation index data is determined based on the degree of matching between the total length of the snow instability triggering time window and the preset reasonable time window length range. When the total length of the snow instability triggering time window falls within the preset reasonable time window length range, the time rationality evaluation index data is marked as reasonable. When the total length of the snow instability triggering time window exceeds the preset reasonable time window length range, the time rationality evaluation index data is marked as unreasonable. The preset reasonable time window length range is determined based on the duration of the rainfall event and the physical evolution timescale of snow instability.

[0122] Step S514: Perform instability influence range boundary shape complexity analysis based on the instability influence range boundary coordinate sequence data at each time node to obtain instability influence range boundary shape complexity index data.

[0123] In this embodiment, the boundary shape complexity index of the instability impact range is calculated by dividing the perimeter of the area enclosed by the boundary coordinate sequence data of the instability impact range at each time point by the area value. The larger the boundary shape complexity index, the more complex the boundary shape. The boundary shape complexity index of the instability impact range is measured in units of 1 meter.

[0124] Step S515: Perform a comprehensive evaluation and merging process on the time rationality evaluation index data and the instability impact range boundary shape complexity index data to obtain the reliability evaluation result. The reliability evaluation result includes simulation result reliability level identification data and simulation result uncertainty area identification data.

[0125] In this embodiment, when the time rationality evaluation index data is reasonably marked and the instability impact range boundary shape complexity index data does not exceed the preset complexity upper limit, the simulation result reliability level identifier data is set to a high reliability level. Otherwise, the simulation result reliability level identifier data is set to a pending review level. The simulation result uncertainty area identifier data is the boundary area where the instability impact range boundary shape complexity index data exceeds the preset complexity upper limit.

[0126] Step S610: After generating the simulation result data set of snow surface rain disaster-causing process based on the disaster risk area distribution data and the disaster risk evolution path data, perform output presentation processing of the simulation result data set of snow surface rain disaster-causing process to obtain visualized output data.

[0127] In this embodiment, the output presentation processing converts the simulation results data set of the snow surface rain-induced disaster process into a visual data format for display on the display terminal. The display terminal is the display device of the snow disaster monitoring and management center. The visual output data includes three parts: snow layer instability trigger time window annotation data, instability impact range boundary superimposed display data, and instability impact range expansion trend curve display data.

[0128] Step S611: Extract the earliest and latest snow instability trigger time points from the snow instability trigger time window data in the snow surface rain-induced disaster simulation result data set, and extract the boundary coordinate sequence data of the instability influence range at each time node from the snow instability influence range data.

[0129] In this embodiment, both the earliest and latest trigger times for snow layer instability are timestamp data. The boundary coordinate sequence data of the instability impact range is an ordered list of boundary spatial coordinates at each time point. All of the above data were read by field name from the simulation results dataset of the snow surface rain-induced disaster process.

[0130] Step S612: Based on the earliest and latest snow instability trigger time data, perform the annotation and display processing of the snow instability trigger time window to obtain the snow instability trigger time window annotation data.

[0131] In this embodiment, the snow instability trigger time window annotation data is a text annotation data object containing the earliest and latest snow instability trigger times. This text annotation data object is displayed on the terminal as a highlighted segment on a timeline, with the segment's starting point corresponding to the earliest snow instability trigger time and the segment's ending point corresponding to the latest snow instability trigger time.

[0132] Step S613: Based on the coordinate sequence data of the boundary of the instability influence range at each time node, perform the superposition display processing of the boundary of the instability influence range on the spatial base map of the area to be simulated to obtain the superposition display data of the boundary of the instability influence range.

[0133] In this embodiment, the spatial base map is a digital elevation model (DEM) rendering of the area to be simulated. The DEM rendering uses color coding to represent the terrain undulations of the area to be simulated. The overlay display processing converts the coordinate sequence data of the boundary of the instability influence range at each time point into polygon overlay layer data on the base map. Each polygon in the polygon overlay layer data corresponds to the boundary of the instability influence range at a time point. Polygons at different time points are distinguished by different colors, with the colors ranging from light to dark representing the expansion of the instability influence range over time.

[0134] Step S614: Based on the cumulative curve data of the expansion of the instability impact range over time, perform curve plotting processing of the expansion trend of the instability impact range to obtain the display data of the expansion trend curve of the instability impact range.

[0135] In this embodiment, the trend curve display data of the expansion of the instability impact range is a curve data object with time as the horizontal axis and the impact range as the vertical axis. The curve in the curve data object is formed by connecting all data points in the cumulative curve data of the expansion of the instability impact range over time. The curve data object marks the data point positions corresponding to the earliest and latest trigger times of snow layer instability.

[0136] Step S615: Perform comprehensive view layout processing on the snow layer instability trigger time window annotation data, the instability influence range boundary superimposed display data, and the instability influence range expansion trend curve display data to generate the visualization output data. The visualization output data is used to present the full picture of the simulation results of the snow surface rain disaster process on the display terminal.

[0137] In this embodiment, the integrated view layout processing combines the snow layer instability trigger time window annotation data, the instability impact range boundary overlay display data, and the instability impact range expansion trend curve display data into a multi-view linked layout data object. The three views in the multi-view linked layout data object share a timeline. When the user selects a time node in any view, the other two views synchronously update the display content corresponding to that time node. This multi-view linked layout data object is the visualization output data.

[0138] Step S710: While performing the snow layer mechanical stability assessment based on the snow layer internal liquid water migration flux distribution data and the snow layer internal solid-ice phase change conversion rate distribution data, perform the correlation state tracking processing between the snow layer internal liquid water migration flux distribution data and the snow layer internal solid-ice phase change conversion rate distribution data to obtain the correlation state tracking result.

[0139] In this embodiment, the correlation state tracking process performs a layer-by-layer analysis of the spatial correlation state between liquid water migration flux distribution data and solid ice phase change conversion rate distribution data, identifies the directional cooperative region and directional antagonistic region between liquid water migration and solid ice phase change, and tracks the dominant state transfer path in the vertical direction of each spatial unit.

[0140] Step S711: Based on the liquid water migration direction data and the liquid water migration amount data in the liquid water migration flux distribution data inside the snow layer, determine the active liquid water migration region and the passive liquid water migration region at each layer position of each spatial unit. The active liquid water migration region is the region where liquid water moves downward along the direction of gravity and is dominant. The passive liquid water migration region is the region where liquid water is forced to move upward due to the volume expansion caused by the ice phase change and is dominant.

[0141] In this embodiment, for each spatial unit and each layer location, when the liquid water migration direction data is downward and the liquid water migration amount data exceeds a preset active migration threshold, the layer location is marked as an active liquid water migration region. When the liquid water migration direction data is upward and the liquid water migration amount data exceeds a preset passive migration threshold, the layer location is marked as a passive liquid water migration region. Both the preset active migration threshold and the preset passive migration threshold are empirical parameters.

[0142] Step S712: Based on the solid ice phase change type data and the solid ice phase change degree data in the solid ice phase change rate distribution data inside the snow layer, determine the solid ice melting-dominant region and the solid ice freezing-dominant region at each layer location of each spatial unit. The solid ice melting-dominant region is the region where the solid ice to liquid water conversion rate exceeds the reverse conversion rate, and the solid ice freezing-dominant region is the region where the liquid water to solid ice conversion rate exceeds the reverse conversion rate.

[0143] In this embodiment, for each spatial unit and each layer location, when the solid-state ice phase change type data is of the melting type and the solid-state ice phase change degree data exceeds a preset melting-dominant threshold, the layer location is marked as a solid-state ice melting-dominant region. When the solid-state ice phase change type data is of the freezing type and the solid-state ice phase change degree data exceeds a preset freezing-dominant threshold, the layer location is marked as a solid-state ice freezing-dominant region. Both the preset melting-dominant threshold and the preset freezing-dominant threshold are empirical parameters.

[0144] Step S713: Perform a layer-by-layer superposition judgment process on the active migration region of liquid water, the passive migration region of liquid water, the melting-dominant region of solid ice, and the freezing-dominant region of solid ice according to the spatial unit position to obtain directional cooperative region data and directional antagonistic region data between the migration direction of liquid water and the phase change direction of solid ice.

[0145] In this embodiment, the layer-by-layer stacking judgment process checks the marking of each layer position in each spatial unit within either the active or passive migration region of liquid water, and then checks its marking within either the melting or freezing region of solid ice. When a layer position is simultaneously located in both the active liquid water migration region and the melting region of solid ice, it is marked as a directionally cooperative region. Similarly, when it is simultaneously located in both the passive liquid water migration region and the freezing region of solid ice, it is marked as a directionally cooperative region. When it is simultaneously located in both the active liquid water migration region and the freezing region of solid ice, it is marked as a directionally antagonistic region. Finally, when it is simultaneously located in both the passive liquid water migration region and the melting region of solid ice, it is marked as a directionally antagonistic region. All directionally cooperative regions constitute directionally cooperative region data, and all directionally antagonistic regions constitute directionally antagonistic region data.

[0146] Step S714: Based on the directional coordination region data, perform positive promotion relationship extraction processing between liquid water migration flux and solid ice phase change conversion rate to obtain positive promotion relationship distribution data. The positive promotion relationship distribution data includes mutually reinforcing region identifiers of liquid water migration promoting ice melting and ice melting promoting liquid water migration.

[0147] In this embodiment, for each directional cooperation region in the directional cooperation region data, the positive promotion relationship type is determined based on the specific combination of markers for that directional cooperation region. When the directional cooperation region is simultaneously located in both the liquid phase water active migration region and the solid phase ice melting dominant region, the positive promotion relationship type is liquid phase water migration promoting ice melting. When the directional cooperation region is simultaneously located in both the liquid phase water passive migration region and the solid phase ice freezing dominant region, the positive promotion relationship type is ice freezing promoting liquid phase water migration. All positive promotion relationship types and their corresponding region identifiers constitute the positive promotion relationship distribution data.

[0148] Step S715: Based on the directional resistance region data, perform reverse inhibition relationship extraction processing between liquid phase water migration flux and solid phase ice phase change conversion rate to obtain reverse inhibition relationship distribution data. The reverse inhibition relationship distribution data includes mutually weakening region identifiers of liquid phase water migration inhibiting freezing and freezing inhibiting liquid phase water migration.

[0149] In this embodiment, for each directional opposition region in the directional opposition region data, the type of reverse inhibition relationship is determined based on the specific combination of markers for that directional opposition region. When the directional opposition region is simultaneously located in both the liquid phase water active migration region and the solid phase ice freezing dominant region, the reverse inhibition relationship type is liquid phase water migration inhibiting ice freezing. When the directional opposition region is simultaneously located in both the liquid phase water passive migration region and the solid phase ice melting dominant region, the reverse inhibition relationship type is ice melting inhibiting liquid phase water migration. All reverse inhibition relationship types and their corresponding region identifiers constitute the reverse inhibition relationship distribution data.

[0150] Step S716: Based on the positive promotion relationship distribution data and the negative inhibition relationship distribution data, perform relative dominance ranking processing between liquid water migration flux and solid ice phase change conversion rate at each layer position of each spatial unit to obtain the relative dominance ranking result for each layer position of each spatial unit. The relative dominance ranking result is used to indicate which of the liquid water migration flux and the solid ice phase change conversion rate is dominant at the current layer position.

[0151] In this embodiment, for each spatial unit and each layer location, the relative magnitudes of the liquid water migration data and the solid-ice phase transition degree data at that layer location are compared. When the liquid water migration data is greater than the solid-ice phase transition degree data, the relative dominance ranking result is liquid water migration dominant. When the solid-ice phase transition degree data is greater than the liquid water migration data, the relative dominance ranking result is solid-ice phase transition dominant. When the two are equal, the relative dominance ranking result is a balanced state.

[0152] Step S717: Based on the relative dominance ranking results, perform dominant state classification processing for each layer position of each spatial unit to obtain dominant state classification set data, which includes liquid phase water migration dominant state classification results and solid phase ice phase change dominant state classification results.

[0153] In this embodiment, for each spatial unit and each layer location, when the relative dominance ranking result is liquid water migration dominant, the dominant state classification result for that layer location is marked as liquid water migration dominant. When the relative dominance ranking result is solid-ice phase change dominant, the dominant state classification result for that layer location is marked as solid-ice phase change dominant. All dominant state classification results constitute the dominant state classification set data.

[0154] Step S718: Based on the dominant state classification results of each spatial unit and each layer position in the dominant state classification set data, perform the dominant state transfer path generation process in the vertical direction within each spatial unit to obtain the dominant state transfer path data in the vertical direction of each spatial unit. The dominant state transfer path data in the vertical direction reflects the alternating transformation law of liquid phase water migration dominant state and solid phase ice phase change dominant state in the vertical direction of the snow layer with depth.

[0155] In this embodiment, for each spatial unit, starting from the first layer on the snow surface, the dominant state classification results of each layer position in the spatial unit are read layer by layer along the depth direction. The dominant state classification results are arranged sequentially along the depth direction to form the vertical dominant state transition path data of the spatial unit. The vertical dominant state transition path data is a sequence data structure, where each element in the sequence is the dominant state classification result of a layer position.

[0156] Step S719: Perform spatial stitching processing on the vertical dominant state transition path data of each spatial unit according to the horizontal arrangement order of the spatial units to obtain the spatial dimension dominant state distribution map data of the entire snow layer in the simulated area. The spatial dimension dominant state distribution map data is used to indicate the overall spatial pattern of the liquid phase water migration dominant region and the solid phase ice phase change dominant region in the snow layer thickness direction.

[0157] In this embodiment, the spatial stitching process arranges the vertical dominant state transition path data of all spatial units according to the row and column indices of the spatial units, forming a spatial dimension dominant state distribution map data. The spatial dimension dominant state distribution map data is a three-dimensional matrix data structure, with the first dimension being the row index of the spatial units, the second dimension being the column index of the spatial units, and the third dimension being the hierarchical position index. Each element in the three-dimensional matrix is ​​the dominant state classification result.

[0158] Step S7110: Based on the spatial dimension dominant state distribution map data, perform the interaction evolution stage division processing between the liquid phase water migration flux and the solid phase ice phase change conversion rate inside the snow layer, and generate the associated state tracking result. The associated state tracking result includes the initial cooperative stage identifier, the alternating dominant stage identifier, and the final equilibrium stage identifier between liquid phase water migration and solid phase ice phase change.

[0159] In this embodiment, the interaction evolution stage segmentation process performs stage identification along the time dimension on the spatial dominant state distribution map data. When the spatial distribution ratio of the liquid phase water migration-dominant state classification result and the solid phase ice phase change-dominant state classification result in the spatial dominant state distribution map data remains relatively stable, this time period is marked as the initial synergistic stage. When the spatial distribution ratio of the two dominant states in the spatial dominant state distribution map data undergoes rapid alternation, this time period is marked as the alternating dominance stage. When the spatial distribution ratio of the two dominant states in the spatial dominant state distribution map data tends to stabilize again, this time period is marked as the final equilibrium stage. All stage identifiers constitute the associated state tracking result.

[0160] Step S810: After performing the snow layer mechanical stability assessment processing based on the snow layer internal liquid water migration flux distribution data and the snow layer internal solid-ice phase transformation rate distribution data, and before performing the disaster-causing process identification processing based on the snow layer internal pore water pressure field distribution data and the snow layer internal effective stress field distribution data, the transmission path identification processing between the snow layer internal pore water pressure field distribution data and the snow layer internal effective stress field distribution data is performed to obtain the transmission path identification result.

[0161] In this embodiment, the transmission path identification process identifies the spatial transformation path between pore water pressure field distribution data and effective stress field distribution data, and determines the priority transmission channel and secondary transmission channel for pore water pressure to be converted into effective stress inside the snow layer.

[0162] Step S811: Extract the pore water pressure gradient direction data of each layer position of each spatial unit from the pore water pressure field distribution data inside the snow layer. The pore water pressure gradient direction data reflects the direction of the most drastic change in pore water pressure inside the snow layer.

[0163] In this embodiment, for each spatial unit and each layer location, the pore water pressure difference between that layer location and its six adjacent layer locations (up, down, left, right, front, and back) is calculated. The adjacent direction with the largest pore water pressure difference is taken as the pore water pressure gradient direction data.

[0164] Step S812: Extract the effective stress gradient direction data of each layer position of each spatial unit from the effective stress field distribution data inside the snow layer. The effective stress gradient direction data reflects the direction of the most drastic change in effective stress in the space inside the snow layer.

[0165] In this embodiment, for each spatial unit and each layer location, the effective stress difference between that layer location and its six adjacent layer locations (up, down, left, right, front, and back) is calculated. The adjacent direction with the largest effective stress difference is taken as the effective stress gradient direction data.

[0166] Step S813: Perform directional comparison processing on the pore water pressure gradient direction data and the effective stress gradient direction data according to the spatial unit position to obtain the directional angle data between the pore water pressure gradient direction and the effective stress gradient direction at each layer position of each spatial unit.

[0167] In this embodiment, the direction angle data is calculated by performing a dot product operation between the unit direction vector corresponding to the pore water pressure gradient direction data and the unit direction vector corresponding to the effective stress gradient direction data, and then taking the inverse cosine of the dot product result to obtain the direction angle data. The direction angle data is measured in degrees and ranges from 0 to 180 degrees.

[0168] Step S814: Determine the direct transmission direction region and the indirect transmission direction region between pore water pressure and effective stress based on the directional angle data. The direct transmission direction region is the region where the directional angle meets the directional consistency condition, and the indirect transmission direction region is the region where the directional angle does not meet the directional consistency condition.

[0169] In this embodiment, the directional consistency condition is that the directional angle data does not exceed a preset directional consistency angle threshold. The preset directional consistency angle threshold is 45 degrees. When the directional angle data does not exceed 45 degrees, the layered location is marked as a directly transmitted directional region. When the directional angle data exceeds 45 degrees, the layered location is marked as an indirectly transmitted directional region.

[0170] Step S815: Extract the conversion path node data of pore water pressure directly converting into effective stress from the directly transmitted direction region to obtain the direct conversion path node set data.

[0171] In this embodiment, for each layer location in the direct transfer direction region, the spatial coordinates of that layer location, along with the pore water pressure value and the effective stress value, are combined to form conversion path node data. All conversion path node data constitute the direct conversion path node set data.

[0172] Step S816: Extract the indirect conversion path node data from the indirect transmission direction region, which shows that the pore water pressure is converted into effective stress after being transmitted through the intermediate medium, and obtain the indirect conversion path node set data.

[0173] In this embodiment, for each layer location in the indirect transmission direction region, the spatial coordinates of that layer location, the pore water pressure value, the intermediate medium transmission parameters, and the effective stress value are combined to form indirect conversion path node data. All indirect conversion path node data constitute the indirect conversion path node set data.

[0174] Step S817: Based on the direct conversion path node set data and the indirect conversion path node set data, perform overall topology construction processing of the pore water pressure transmission path to obtain pore water pressure transmission path topology map data, which includes direct transmission branches and indirect transmission branches.

[0175] In this embodiment, a direct transmission branch is a path formed by sequentially connecting spatially adjacent conversion path nodes in the direct conversion path node set data. An indirect transmission branch is a path formed by sequentially connecting spatially adjacent indirect conversion path nodes in the indirect conversion path node set data. All direct and indirect transmission branches constitute the pore water pressure transmission path topology data.

[0176] Step S818: Based on the path length data of each transmission branch and the transmission efficiency data of the nodes on each transmission branch in the pore water pressure transmission path topology map data, perform the transmission contribution calculation of each transmission branch to obtain the transmission contribution ranking result of each transmission branch.

[0177] In this embodiment, the transfer contribution is calculated by dividing the average ratio of the effective stress value to the pore water pressure value of all nodes on the transfer branch by the path length data. The transfer contribution ranking results are arranged from largest to smallest.

[0178] Step S819: Based on the ranking results of the transmission contribution, select the transmission branch with the highest transmission contribution as the main transmission path data for the transmission of pore water pressure to effective stress.

[0179] In this embodiment, the transmission branch ranked first in the transmission contribution ranking results is marked as the main transmission path data.

[0180] Step S8110: The main transmission path data and the key indirect conversion node data in the indirect conversion path node set data are associated and stored to generate the transmission path identification result. The transmission path identification result is used to indicate the priority transmission channel and secondary transmission channel of pore water pressure being converted into effective stress inside the snow layer.

[0181] In this embodiment, the key indirect conversion node data refers to the indirect conversion path nodes in the indirect conversion path node set data whose transmission efficiency exceeds a preset transmission efficiency threshold. The transmission path identification result is a data object containing both the main transmission path data and the key indirect conversion node data.

[0182] Step S910: While performing the disaster-causing process identification processing based on the pore water pressure field distribution data and the effective stress field distribution data inside the snow layer, spatial matching processing is performed between the disaster-causing risk area distribution data and the location data of existing weak structures inside the area to be simulated to obtain the spatial matching result.

[0183] In this embodiment, spatial matching processing performs spatial overlay analysis on the distribution data of disaster-risk areas and the location data of existing weak structures within the area to be simulated. The location data of existing weak structures within the area to be simulated consists of location information of weak structures such as bedrock fissures, faults, and weak interlayers at the bottom of the snow layer, obtained in advance through geological exploration. Spatial matching processing is used to identify the spatial correspondence between disaster-risk areas and the locations of existing weak structures, and to determine the evolutionary tendency of disaster-risk areas under the influence of the locations of existing weak structures.

[0184] Step S911: Extract the spatial boundary coordinate data of each disaster risk area from the disaster risk area distribution data. The spatial boundary coordinate data is used to define the external contour range of each disaster risk area within the area to be simulated.

[0185] In this embodiment, the disaster risk area distribution data is a set of overlapping spatial units. For each spatially continuous set of overlapping spatial units, the spatial coordinates of the overlapping spatial units located at the edge of the set are extracted. The edge spatial coordinates are then connected sequentially to form the spatial boundary coordinate data of the disaster risk area. The spatial boundary coordinate data is an ordered coordinate list.

[0186] Step S912: Obtain the location data of existing weak structures within the area to be simulated. The location data of existing weak structures includes the location data of existing cracks, faults, and weak interlayers within the area to be simulated.

[0187] In this embodiment, the location data of existing weak structures are read from the geological exploration database of the area to be simulated. The fracture location data is linear feature data, with each fracture represented by its starting and ending coordinates. The fault location data is planar feature data, with each fault represented by a sequence of boundary coordinates. The weak interlayer location data is planar feature data, with each weak interlayer represented by a sequence of boundary coordinates. All of the above-mentioned existing weak structure location data are in vector data format.

[0188] Step S913: Spatial overlay processing is performed on the spatial boundary coordinate data of the disaster risk area and the location data of the existing weak structure to obtain the spatial overlap area data between the disaster risk area and the location of the existing weak structure.

[0189] In this embodiment, spatial overlay processing performs spatial intersection analysis on the polygons enclosed by the spatial boundary coordinate data of each disaster-risk area and the linear or planar features represented by the location data of each existing weak structure. When the polygons of the disaster-risk area and the location features of the existing weak structure spatially intersect, the spatial extent of the intersection is extracted as spatially overlapping area data. The spatially overlapping area data is a set of polygon features.

[0190] Step S914: Based on the spatial overlap area data, determine the identification data of the disaster-causing risk sub-regions containing existing weak structure locations within the disaster-causing risk area, and obtain the risk sub-region set data containing weak structure locations.

[0191] In this embodiment, for each disaster-risk area, it is checked whether the polygon enclosed by the spatial boundary coordinates of the disaster-risk area intersects with any spatially overlapping area data. When an intersection exists, the disaster-risk area is marked as a disaster-risk sub-region containing the location of weak structures. The identifiers of all disaster-risk sub-regions constitute the risk sub-region set data containing the location of weak structures.

[0192] Step S915: Based on the spatial overlap area data, determine the weak structure location identification data located within the disaster risk area in the existing weak structure locations, and obtain the set data of weak structure locations located within the risk area.

[0193] In this embodiment, for each existing weak structure location feature, it is checked whether the feature intersects with any spatially overlapping area data. When an intersection exists, the existing weak structure location feature is marked as a weak structure location located within the risk area. All weak structure location identifiers located within the risk area constitute the set data of weak structure locations located within the risk area.

[0194] Step S916: Based on the risk sub-region set data containing the location of weak structures and the set data of weak structure locations located within the risk region, perform spatial correlation quantification processing between the disaster-causing risk region and the existing weak structure to obtain the spatial correlation quantification result. The spatial correlation quantification result includes the spatial overlap area ratio data and spatial overlap pattern data between each disaster-causing risk region and the existing weak structure location.

[0195] In this embodiment, for each disaster-risk sub-region containing weak structural locations, the ratio of the sum of the areas covered by all spatially overlapping areas within that sub-region to the total area of ​​that sub-region is calculated to obtain the spatial overlap area ratio data. Simultaneously, the boundary shape features of all spatially overlapping areas within that sub-region are extracted, and these features are compared with the boundary shape features of the sub-region itself to obtain spatial overlap morphology data. The spatial overlap morphology data can be either a strip-like overlap morphology or a cluster-like overlap morphology.

[0196] Step S917: Based on the spatial overlap area ratio data in the spatial correlation degree quantification result, perform internal priority sorting processing on all disaster risk sub-regions containing weak structure locations to obtain internal disaster priority sorting results. The internal disaster priority sorting results indicate the disaster risk sub-region identifier with the highest disaster priority.

[0197] In this embodiment, all disaster-risk sub-regions containing weak structural locations are sorted from largest to smallest by the percentage of their spatial overlap area. The disaster-risk sub-region ranked first in the sorting results is marked as the disaster-risk sub-region with the highest disaster priority. The internal disaster priority sorting result is a list containing all disaster-risk sub-regions arranged from highest to lowest priority.

[0198] Step S918: Based on the spatial overlap pattern data in the spatial correlation degree quantification result, perform disaster risk area boundary expansion possibility determination processing on all disaster risk sub-regions containing weak structure locations to obtain boundary expansion possibility determination results. The boundary expansion possibility determination results indicate whether the disaster risk area has a tendency to expand outward along the existing weak structure location direction.

[0199] In this embodiment, for each disaster-risk sub-region containing the location of a weak structure, when the spatial overlap pattern data is a strip-shaped overlap pattern, the disaster-risk sub-region has a tendency to expand outward along the direction of the existing weak structure location, and the boundary expansion probability determination result is high expansion probability. When the spatial overlap pattern data is a cluster-shaped overlap pattern, the disaster-risk sub-region does not have a tendency to expand outward along the direction of the existing weak structure location, and the boundary expansion probability determination result is low expansion probability.

[0200] Step S919: The internal disaster priority ranking results and the boundary expansion possibility determination results are integrated and merged to generate the spatial matching results. The spatial matching results are used to indicate the spatial correspondence between the disaster risk area and the location of the existing weak structure, as well as the evolution tendency of the disaster risk area under the influence of the location of the existing weak structure.

[0201] In this embodiment, the spatial matching result is a data object containing the internal disaster priority ranking result and the boundary expansion probability determination result. The spatial matching result is used in the subsequent generation and processing of the snow surface rain disaster process simulation result dataset to perform priority weighting and boundary correction on the disaster risk area distribution data.

[0202] In some embodiments, the snow surface rain-induced disaster process simulation system using a fused hydrothermal coupling model for performing the above methods can be any electronic device with data computing, processing, and storage capabilities. This fused hydrothermal coupling model snow surface rain-induced disaster process simulation system can be used to implement the snow surface rain-induced disaster process simulation method using a fused hydrothermal coupling model provided in the above embodiments.

[0203] Typically, a snow surface rain-induced disaster simulation system incorporating a hydrothermal coupling model includes a processor and a memory. The memory may include one or more computer-readable storage media, which may be non-transitory. The memory may also include high-speed random access memory and non-volatile memory, such as one or more disk storage devices or flash memory devices. In some embodiments, the non-transitory computer-readable storage media in the memory is used to store a computer program configured to be executed by one or more processors to implement the above-described snow surface rain-induced disaster simulation method incorporating a hydrothermal coupling model.

[0204] Finally, it should be noted that the above-disclosed embodiments are merely preferred embodiments of the present invention and are only used to illustrate the technical solutions of the present invention, not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for simulating snow surface rain-induced disaster processes using a hydrothermal coupling model, characterized in that, The method includes: Obtain the initial snow layer state data set of the area to be simulated before the occurrence of the rainfall event. The initial snow layer state data set includes spatial distribution data of snow layer thickness, spatial distribution data of snow layer temperature, and spatial distribution data of snow layer liquid water content. The initial snow layer state data set is subjected to snow surface rainfall infiltration process simulation processing to generate data on the distribution of liquid phase water migration flux inside the snow layer and the distribution of solid phase ice phase transformation rate inside the snow layer during the simulated rainfall period. Based on the data on the distribution of liquid water migration flux inside the snow layer and the data on the distribution of solid-ice phase transformation rate inside the snow layer, the mechanical stability assessment of the snow layer is performed to generate data on the distribution of pore water pressure field inside the snow layer and data on the distribution of effective stress field inside the snow layer. Based on the distribution data of pore water pressure field inside the snow layer and the distribution data of effective stress field inside the snow layer, disaster-causing process identification processing is performed to generate disaster-causing risk area distribution data and disaster-causing risk evolution path data for snow layer instability induced by rainfall infiltration. Based on the disaster risk area distribution data and the disaster risk evolution path data, a set of simulation results of the snow surface rain disaster process is generated. The set of simulation results of the snow surface rain disaster process includes snow layer instability trigger time window data and snow layer instability impact range data.

2. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The initial snow layer state data set is subjected to a snow surface rainfall infiltration process simulation process to generate data on the distribution of liquid phase water migration flux inside the snow layer and the distribution of solid-ice phase transformation rate inside the snow layer during the simulated rainfall period, including: The spatial distribution data of snow layer thickness in the initial snow layer state data set is processed by performing regional splitting according to the spatial grid division method to obtain the sequence of spatial units inside the snow layer. Each spatial unit in the sequence of spatial units inside the snow layer has a hierarchical structure attribute. Based on the spatial distribution data of snow layer temperature and the spatial distribution data of snow layer liquid water content, an initial state mapping process is performed on each spatial unit in the sequence of spatial units inside the snow layer to obtain the initial temperature state result and the initial liquid water content state result of each spatial unit at each layer position. Rainfall infiltration boundary condition setting processing is performed on the overall boundary region of the snow layer internal spatial unit sequence to obtain the top rainfall infiltration boundary parameters and the bottom water discharge boundary parameters of the snow layer internal spatial unit sequence. Based on the top rainfall infiltration boundary parameters and the bottom water discharge boundary parameters, liquid phase water migration-driven calculations are performed between the layers in each spatial unit to obtain liquid phase water migration direction data and liquid phase water migration amount data at each time node during the simulated rainfall period. Spatial aggregation processing is performed on each of the liquid water migration data according to the spatial unit arrangement order to obtain the liquid water migration flux distribution data inside the snow layer. The liquid water migration flux distribution data inside the snow layer reflects the overall migration total characteristics of liquid water in the vertical direction and the migration distribution characteristics in the horizontal direction. Based on the liquid water migration direction data and the liquid water migration amount data, combined with the initial temperature state results of each spatial unit's stratification location, a solid-ice phase transition determination process is performed to obtain solid-ice phase transition type data and solid-ice phase transition degree data at each stratification location of each spatial unit. Spatial correlation and aggregation processing is performed on the solid ice phase change type data and the solid ice phase change degree data to obtain the solid ice phase change conversion rate distribution data inside the snow layer. The solid ice phase change conversion rate distribution data inside the snow layer reflects the rate attribute of solid ice to liquid water conversion and the rate attribute of liquid water to solid ice conversion at each layer position of each spatial unit. Based on the liquid water migration flux distribution data inside the snow layer and the solid ice phase change degree data in the solid ice phase change rate distribution data inside the snow layer, the layered temperature state update processing is performed to obtain the updated temperature state data of each layer position in each spatial unit. The updated temperature state data is used as the initial temperature state input data for the next time node. The process of setting the rainfall infiltration boundary conditions and the liquid water migration driving calculation is then performed until the preset simulation termination time is reached. This process yields the distribution data of liquid water migration flux inside the snow layer and the distribution data of solid-ice phase transformation rate inside the snow layer for each time node during the simulated rainfall period.

3. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The process of performing a snow layer mechanical stability assessment based on the distribution data of liquid water migration flux within the snow layer and the distribution data of solid-ice phase transformation rate within the snow layer generates data on the distribution of pore water pressure field and effective stress field within the snow layer, including: Extract the liquid water migration direction data and liquid water migration amount data of each layer position in each spatial unit from the liquid water migration flux distribution data inside the snow layer, and determine the flow thrust distribution data of liquid water acting on the solid skeleton at each layer position based on the liquid water migration direction data and liquid water migration amount data. Based on the solid-ice phase transition type data and solid-ice phase transition degree data in the solid-ice phase transition rate distribution data inside the snow layer, the solid volume ratio update processing of each layer position is performed to obtain the solid volume ratio change sequence data of each layer position. The solid volume ratio change sequence data reflects the solid volume reduction caused by the conversion of solid ice to liquid water and the solid volume increase caused by the conversion of liquid water to solid ice. Based on the solid volume ratio change sequence data of each layer location, the pore space ratio update processing of each layer location is performed to obtain the pore space ratio update data of each layer location. The pore space ratio update data reflects the dynamic change distribution characteristics of the pore space inside the snow layer as the solid-ice phase transformation process occurs. The flow thrust distribution data of liquid water acting on the solid skeleton at each layer location is correlated with the pore space ratio update data at each layer location to obtain the pore water pressure field distribution data at each layer location of each spatial unit. The pore water pressure field distribution data reflects the pressure magnitude distribution and pressure direction distribution of liquid water in the pore space. Based on the solid volume ratio change sequence data of each layer position and the initial gravity load data of each layer position of each spatial unit, the total stress field generation process inside the snow layer is performed to obtain the total stress distribution data of each layer position of each spatial unit. The total stress distribution data includes the vertical stress component caused by the snow layer's own gravity and the boundary stress component applied from the outside. Effective stress distribution data of each layer location in each spatial unit and pore water pressure field distribution data of each layer location in each spatial unit are processed by effective stress solution to obtain effective stress distribution data of each layer location in each spatial unit. The effective stress distribution data reflects the actual stress state and stress transmission path characteristics of the solid skeleton.

4. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The process of identifying the disaster-causing process is performed based on the distribution data of the pore water pressure field inside the snow layer and the distribution data of the effective stress field inside the snow layer, generating disaster-causing risk area distribution data and disaster-causing risk evolution path data for snow layer instability induced by rainfall infiltration, including: The pore water pressure data of each layer in each spatial unit is extracted from the pore water pressure field distribution data inside the snow layer. Based on the pore water pressure data, the pore water pressure exceeding the limit area is determined, and the set data of the pore water pressure exceeding the limit spatial unit and the degree of exceeding the limit of each exceeding spatial unit are obtained. Effective stress numerical data of each layer position of each spatial unit are extracted from the effective stress field distribution data inside the snow layer. Effective stress abrupt change region determination processing is performed based on the effective stress numerical data to obtain effective stress abrupt change spatial unit set data and change degree data of each abrupt change spatial unit. The spatial location cross-matching process is performed on the set of spatial units with pore water pressure exceeding the limit and the set of spatial units with effective stress change to obtain the data of the overlapping area of ​​spatial units with pore water pressure exceeding the limit and effective stress change. The data of the overlapping area of ​​spatial units is used as the disaster risk area distribution data of snow layer instability induced by rainfall infiltration. Based on the spatial positional arrangement of each overlapping spatial unit in the spatial unit overlapping area data, a disaster risk spatial connectivity analysis is performed to obtain a set of spatial connectivity paths between disaster risk areas. The set of spatial connectivity paths includes spatial connection direction data and spatial connection order data between adjacent disaster risk areas. Based on the spatial connection direction data and spatial connection sequence data in the spatial connectivity path set data, disaster risk evolution trend tracking processing is performed to obtain the spatial diffusion direction sequence data and spatial diffusion rate sequence data of disaster risk in the area to be simulated. The disaster risk spatial connectivity path set data, the spatial diffusion direction sequence data, and the spatial diffusion rate sequence data are comprehensively merged and processed to generate the disaster risk evolution path data. The disaster risk evolution path data reflects the spatial trajectory characteristics and expansion speed change characteristics of the disaster risk as it gradually expands from the initial triggering area to the surrounding area.

5. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The data set of simulation results for the snow surface rain-induced disaster process is generated based on the disaster risk area distribution data and the disaster risk evolution path data. This data set includes snow layer instability triggering time window data and snow layer instability impact range data, including: Extract the spatial unit identifier data corresponding to each disaster risk area and the disaster risk start time point data corresponding to the spatial unit identifier data from the disaster risk area distribution data. Perform trigger time sequence sorting processing on each disaster risk area according to the disaster risk start time point data to obtain the time dimension distribution sequence data of each disaster risk area. Based on the time dimension distribution sequence data of each disaster-causing risk area, the snow layer instability triggering time window range is defined to obtain the snow layer instability triggering time window data, which includes the earliest triggering time point data and the latest triggering time point data of snow layer instability. The spatial diffusion direction sequence data and spatial diffusion rate sequence data of disaster risk are extracted from the disaster risk evolution path data. Based on the spatial diffusion direction sequence data and spatial diffusion rate sequence data of disaster risk, the boundary estimation processing of the snow layer instability influence range is performed to obtain the snow layer instability influence range boundary coordinate sequence data of each time node during the simulated rainfall period. Based on the earliest and latest trigger time points in the snow instability trigger time window data, and combined with the boundary coordinate sequence data of the snow instability influence range, dynamic mapping processing of the snow instability influence range over time is performed to obtain the snow instability influence range data. The snow instability influence range data reflects the characteristics of the expansion coverage area of ​​the instability region in the time dimension and the characteristics of the expansion direction in the spatial dimension.

6. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The method further includes: After obtaining the initial snow layer state data set of the area to be simulated before the occurrence of a rainfall event, and before performing snow surface rainfall infiltration process simulation processing on the initial snow layer state data set, a validity determination process for the initial snow layer state data set is performed to obtain a valid initial snow layer state data set, including: A spatial continuity check is performed on the spatial distribution data of snow layer thickness in the initial snow layer state data set to obtain the continuity index data of the spatial distribution data of snow layer thickness at each spatial location; Based on the continuity index data, identify the missing data regions and abnormal data regions in the spatial distribution data of snow layer thickness, and obtain the missing region identification data and abnormal region identification data; Based on the missing region identifier data and the abnormal region identifier data, perform missing data region completion processing and abnormal data region correction processing to obtain corrected snow layer thickness spatial distribution data. Perform a temperature gradient rationality check on the spatial distribution data of snow layer temperature in the initial snow layer state data set to obtain temperature gradient rationality index data at each spatial location; Based on the temperature gradient rationality index data, identify the temperature abrupt change location data and temperature reverse change location data in the snow layer temperature spatial distribution data to obtain the temperature anomaly location set data; Based on the set of temperature anomaly locations, interpolation correction processing is performed on the spatial distribution data of snow layer temperature to obtain the corrected spatial distribution data of snow layer temperature. Perform physical boundary review processing on the spatial distribution data of liquid water content in the snow layer in the initial snow layer state data set to obtain the physical boundary compliance index data of liquid water content at each spatial location; The excessive value data in the spatial distribution data of liquid water content in the snow layer are adjusted according to the physical boundary compliance index data of liquid water content to obtain the adjusted spatial distribution data of liquid water content in the snow layer. The corrected spatial distribution data of snow layer thickness, the corrected spatial distribution data of snow layer temperature, and the adjusted spatial distribution data of snow layer liquid water content are aggregated and merged to generate the effective initial snow layer state data set.

7. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The method further includes: Before performing the snow layer mechanical stability assessment based on the distribution data of liquid water migration flux and solid-ice phase transition rate within the snow layer, a coupling compatibility verification process is performed between the liquid water migration flux and the solid-ice phase transition rate to obtain the coupling compatibility verification results, including: Extract the total liquid water migration data at each time point during the simulated rainfall period from the liquid water migration flux distribution data inside the snow layer, and extract the total solid-ice phase change conversion data at the corresponding time point from the solid-ice phase change conversion rate distribution data inside the snow layer. Based on the total liquid phase water migration data and the total solid phase ice phase transformation data, the time synchronization index data between the total liquid phase water migration data and the total solid phase ice phase transformation data is calculated. The time synchronization index data reflects the relative offset between the total liquid phase water migration data and the total solid phase ice phase transformation data in the time dimension. Based on the liquid water migration direction data of each spatial unit in the liquid water migration flux distribution data and the solid ice phase change type data of each spatial unit in the solid ice phase change rate distribution data, spatial correspondence analysis is performed to obtain spatial consistency score data between the liquid water migration direction and the solid ice phase change type of each spatial unit. Based on the spatial consistency score data, identify the set of spatial units where spatial mismatch occurs between the liquid phase water migration direction and the solid phase ice phase change type, and obtain spatial mismatch unit identification data and mismatch type data corresponding to each mismatch unit; Based on the time synchronization index data and the spatial mismatch unit identification data, a comprehensive judgment process is performed on the coupling relationship between liquid water migration flux and solid ice phase change conversion rate to obtain the coupling coordination verification result. The coupling coordination verification result includes coupling coordination qualified region data and coupling coordination needing adjustment region data.

8. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The method further includes: After performing the disaster-causing process identification processing based on the pore water pressure field distribution data and the effective stress field distribution data within the snow layer, and before generating the snow surface rain disaster-causing process simulation result dataset based on the disaster-causing risk area distribution data and the disaster-causing risk evolution path data, spatial correlation consistency verification processing is performed between the disaster-causing risk area distribution data and the disaster-causing risk evolution path data to obtain spatial correlation consistency verification results, including: The spatial coordinate data of each disaster-causing risk area and the trigger time data of each disaster-causing risk area are extracted from the disaster-causing risk area distribution data. The spatial coordinate data of each path node in the disaster-causing risk evolution path and the time point data corresponding to each path node are extracted from the disaster-causing risk evolution path data. Spatial coordinate matching processing is performed between the spatial location coordinate data of each disaster-causing risk area and the spatial location coordinate data of each path node to obtain spatial overlap data between the disaster-causing risk area and the evolution path node. The trigger time data of each disaster-causing risk area and the time data corresponding to each path node are time-series aligned to obtain time consistency data between the trigger time of the disaster-causing risk area and the time of the evolution path node. Based on the spatial overlap data and the temporal consistency data, a comprehensive spatial correlation scoring process is performed between the disaster risk area distribution data and the disaster risk evolution path data to obtain the spatial correlation consistency verification result. The spatial correlation consistency verification result includes spatial correlation completely matched area data and spatial correlation deviation area data. Based on the spatial correlation of the biased area data, spatial location correction processing is performed between the disaster risk area distribution data and the disaster risk evolution path data to obtain the corrected disaster risk area distribution data and the corrected disaster risk evolution path data.

9. The method for simulating snow surface rain-induced disaster processes using a fusion hydrothermal coupling model according to claim 1, characterized in that, The method further includes: After generating the snow surface rain disaster-causing process simulation result dataset based on the disaster-causing risk area distribution data and the disaster-causing risk evolution path data, a reliability evaluation process is performed on the snow surface rain disaster-causing process simulation result dataset to obtain the reliability evaluation results, including: Extract the earliest and latest snow instability trigger time points from the snow instability trigger time window data in the snow surface rain disaster simulation result data set, and extract the boundary coordinate sequence data of the instability influence range at each time node from the snow instability influence range data. Based on the earliest and latest trigger time data of snow layer instability, the trigger time window length is calculated to obtain the total length of the snow layer instability trigger time window. Based on the boundary coordinate sequence data of the instability influence range at each time node, the evolution trend of the total area of ​​the instability influence range is calculated to obtain the cumulative curve data of the expansion of the instability influence range over time. Based on the total length of the snow layer instability triggering time window and the cumulative curve data of the instability influence range over time, the time rationality evaluation data set of the snow surface rain disaster process simulation results is performed to obtain time rationality evaluation index data. Based on the boundary coordinate sequence data of the instability impact range at each time point, perform instability impact range boundary shape complexity analysis to obtain instability impact range boundary shape complexity index data; The time rationality evaluation index data and the instability impact range boundary shape complexity index data are comprehensively evaluated and merged to obtain the reliability evaluation result. The reliability evaluation result includes simulation result reliability level identification data and simulation result uncertainty area identification data.

10. A simulation system for snow surface rain-induced disaster processes integrating a hydrothermal coupling model, characterized in that, The method includes a processor and a computer-readable storage medium storing machine-executable instructions, which, when executed by a computer, implement the snow surface rain-induced disaster simulation method according to any one of claims 1-9, which incorporates a water-thermal coupling model.

Images