A method for quantifying the spatiotemporal pattern of landslide disasters at a wide range of scales

Abstract

The present application relates to the technical field of geological disaster monitoring and prevention, and discloses a kind of wide-area scale landslide disaster spatio-temporal pattern quantification method, comprising: setting virtual space base map division physical calculation grid, access deformation and rainfall time series;Perform resampling, calculate maximum short-time distance cross-correlation coefficient and deformation rate absolute value average generation nonlinear coupling index;The index is mapped to generate spatio-temporal third-order tensor, and the spatial energy projection matrix is restored by high-order singular value decomposition;According to energy threshold, generate grid splitting or merging instruction to update grid structure;Extract matrix to calculate the coordinate bias of energy geometric center, generate spatio-temporal pattern migration vector, output rate and direction parameter and dispatch external equipment.The present application eliminates the difference of data acquisition frequency and local noise, realizes dynamic distribution of calculation grid, and converts energy distribution into kinematic parameters to directly guide equipment scheduling.

Description

Technical Field

[0001] This invention relates to the field of geological disaster monitoring and prevention technology, specifically a method for quantifying the spatiotemporal patterns of landslide disasters on a wide-area scale. Background Technology

[0002] Wide-scale landslide disaster monitoring typically requires comprehensive analysis of surface deformation and rainfall data. Existing monitoring methods, when processing multi-source sensor data, suffer from difficulties in aligning time series data from different modes across the same time dimension due to differences in data acquisition frequencies between rainfall and surface deformation monitoring equipment. This directly impacts the accuracy of calculating the coupling relationship between rainfall and deformation.

[0003] When conducting quantitative analysis of geological hazards over a wide spatial area, the base map grid is typically divided statically. This fixed grid resolution is difficult to adapt to the differences in energy distribution of geological activity in different regions, easily leading to excessive consumption of computational resources in stable and large background areas, while failing to provide the required computational accuracy in areas of actual geological activity undergoing deformation. In addition, multi-source wide-area monitoring data are often mixed with local random noise, and due to the inconsistent spatial resolution of each sensor, there are data misalignment problems, which further interfere with the extraction of real geological activity signals.

[0004] Existing quantitative analysis results are mostly presented in the form of static spatial probability distributions, lacking analysis of the dynamic migration characteristics of surface deformation during spatiotemporal evolution. Because the high-dimensional spatial distribution state cannot be transformed into specific kinematic parameters, existing methods are unable to calculate the absolute migration rate and evolution direction between continuous time slices, thus failing to provide direct quantitative control basis for the dynamic scheduling and deployment of external physical monitoring equipment. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for quantifying the spatiotemporal patterns of landslide disasters at a wide-scale, solving the problems of inconsistent data acquisition frequencies from multiple sources of sensors leading to difficulties in data alignment and calculation, difficulty in unifying the resolution of wide-area spatial grids which are susceptible to local random noise interference, and the waste of computational resources caused by static grid partitioning and the lack of clear kinematic parameters to directly guide the scheduling of external physical monitoring equipment in wide-scale landslide monitoring.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale, comprising:

[0007] Set a global virtual space base map, divide the physical calculation grid according to the initial resolution, and connect the surface deformation time series and rainfall time series within the current time sliding window; Within the physical computing grid, time-dimensional resampling operations are performed on the land surface deformation time series and rainfall time series to calculate the maximum short-term cross-correlation coefficient. Based on the maximum short-term cross-correlation coefficient and the mean absolute value of deformation rate, the nonlinear coupling index between rainfall and deformation is calculated. The nonlinear coupling index of rainfall and deformation is numerically mapped to the global virtual space base map to generate a third-order spatiotemporal tensor. The third-order spatiotemporal tensor is then subjected to high-order singular value decomposition to reconstruct and calculate the spatial energy projection matrix. The energy distribution within the spatial energy projection matrix is ​​determined based on the set energy threshold, and a mesh splitting command or mesh merging command is generated to update the physical calculation mesh structure corresponding to the next time sliding window. Extract the spatial energy projection matrix from continuous time slices, calculate the spatial coordinate offset of the energy geometric centroid, generate the spatiotemporal pattern migration vector, output the velocity and direction parameters corresponding to the spatiotemporal pattern migration vector, and schedule external physical monitoring equipment based on the velocity and direction parameters.

[0008] The process of setting up a global virtual spatial base map and dividing it into physical computation grids specifically includes acquiring the geographic boundary data of the wide-area geological region to be monitored, establishing a two-dimensional spatial matrix under a defined geographic coordinate system, defining the first spatial dimension of the two-dimensional spatial matrix as the longitude direction and the second spatial dimension as the latitude direction. The number of longitude grids corresponding to the highest spatial computational precision is used as the maximum computational resolution index of the first spatial dimension, and the number of latitude grids corresponding to the highest spatial computational precision is used as the maximum computational resolution index of the second spatial dimension, generating a global virtual spatial base map represented by a set of two-dimensional grid coordinates.

[0009] During the system initialization phase, the global virtual spatial base map is aggregated according to a preset initial spatial resolution. This maps the coordinates of multiple adjacent grids in the global virtual spatial base map to the same physical computation grid, generating an initial physical computation grid with spatial coordinate indices. The principle behind this setting is to establish a unified and fixed underlying coordinate reference system to adapt to dynamically changing computation grid structures and to extract the physical state of geological bodies and hydrological elements in time segments.

[0010] Within the physical computation grid, a time-dimensional resampling operation is performed on the land surface deformation time series and rainfall time series. Specifically, this involves obtaining the average soil permeability and porosity, and calculating the effective physical seepage period required for rainfall moisture to infiltrate to the estimated slip surface based on the average soil permeability and porosity. The effective physical seepage period is then used as the cumulative window length and continuously slid across the rainfall time series to calculate the cumulative rainfall value within the window, generating a cumulative rainfall series.

[0011] For land surface deformation time series, a conformal piecewise cubic interpolation algorithm is used for continuous resampling in the time domain. A piecewise cubic polynomial is constructed between two adjacent discrete land surface deformation data points, resampling the land surface deformation time series into a daily deformation rate series with the same data length as the cumulative rainfall series. This operation is used to eliminate the differences in data acquisition frequency between different types of physical sensors.

[0012] The process of calculating the maximum short-term cross-correlation coefficient and the nonlinear coupling index between rainfall and deformation includes introducing a time-shifting variable to perform relative shifting on the cumulative rainfall series and the daily deformation rate series. The maximum cross-correlation coefficient between the cumulative rainfall series and the daily deformation rate series under different time-shifting variables is extracted as the maximum short-term cross-correlation coefficient. The absolute values ​​of all values ​​in the daily deformation rate series within the current time sliding window are extracted, and the arithmetic mean of all absolute values ​​is calculated to generate the mean absolute value of the deformation rate.

[0013] The absolute value of the maximum short-term cross-correlation coefficient is extracted, and then multiplied by the mean absolute value of the deformation rate to generate the nonlinear coupling index between rainfall and deformation for the current physical computing grid. This calculation process uses the mean absolute value of the deformation rate as a weight to filter out background static areas with high correlation coefficients but weak actual deformation magnitudes, while retaining geologically active areas with significant deformation and strong correlation with rainfall.

[0014] The nonlinear coupling index between rainfall and deformation is numerically mapped to a global virtual space base map to generate a spatiotemporal third-order tensor. This involves establishing a standard spatiotemporal third-order tensor storage space with the maximum computational resolution index of the first and second spatial dimensions as the first two spatial dimensions, and the set historical sequence buffer window depth as the third dimension. The latitude and longitude coordinates of the global virtual space base map covered by the physical computation grid are extracted. The nonlinear coupling index between rainfall and deformation from the physical computation grid is copied and filled into the coordinates of all corresponding virtual space base map grids covered by the physical computation grid, generating two-dimensional data slices with regular internal data structures. These two-dimensional data slices are then pushed into the final time index layer of the standard spatiotemporal third-order tensor storage space.

[0015] A higher-order singular value decomposition (SVD) is performed on the third-order spatiotemporal tensor to reconstruct the spatial energy projection matrix. Specifically, the generated third-order spatiotemporal tensor is transformed into a core tensor and a modular product of three orthogonal factor matrices corresponding to the longitude, latitude, and time dimensions, respectively. The squares of the Frobenius norms of each principal component term in the core tensor are calculated as the energy values ​​of the corresponding principal component terms. The principal component terms of the core tensor are then sorted in descending order according to their energy values. The top few principal component terms whose cumulative energy percentage reaches the energy retention threshold are extracted to form a truncated core tensor.

[0016] By performing inverse modular multiplication on the truncated core tensor and the orthogonal factor matrices corresponding to the three dimensions, a noise-removed spatiotemporal third-order tensor is reconstructed. The two-dimensional data surface corresponding to the current time slice is extracted as the spatial energy projection matrix. All elements in the spatial energy projection matrix are iterated and compared; elements less than zero are forced to zero, while elements greater than or equal to zero undergo non-negative thresholding. The principle behind this numerical mapping and higher-order decomposition reconstruction lies in addressing the data misalignment problem caused by inconsistent grid resolution, removing local random noise signals in the wide-area space through tensor dimensionality reduction.

[0017] The steps for generating the energy thresholds described above are as follows: Extract the latitude and longitude coordinates of the global virtual space map contained in each physical computing grid; sum the energy values ​​corresponding to these coordinates in the spatial energy projection matrix; divide by the number of virtual space map grids contained in the physical computing grid to calculate the energy density value of each physical computing grid. Calculate the expected value and standard deviation of the energy density values ​​of all current physical computing grids. Calculate the upper energy threshold using the expected value, standard deviation, and upper limit control coefficient; and calculate the lower energy threshold using the expected value, standard deviation, and lower limit control coefficient.

[0018] The process of generating grid splitting or merging commands involves several steps. First, when the energy density of a physical computing grid is determined to be greater than the upper energy threshold, a grid splitting command is triggered. This command divides the current physical computing grid into four equal sub-grids along the longitude and latitude centerlines, based on quadtree spatial indexing rules. Then, four adjacent physical computing grids of the same level belonging to the same parent node are searched. If the energy density of these four adjacent grids is determined to be less than the lower energy threshold, a grid merging command is triggered. This command merges the four sub-grids back into a single parent grid at the next higher level. This operation establishes a dynamic allocation mechanism for computing resources, allocating high-precision computing grids to active geological deformation zones with high-energy characteristics and reducing the grid resolution in low-energy, stable background zones.

[0019] The spatial coordinate offset of the energy geometric centroid is calculated by using the actual geographic longitude and latitude coordinates of each grid element in the spatial energy projection matrix as the base position, and the denoised energy intensity value of the corresponding grid element as the weight. The longitude coordinates of each grid element are multiplied by the denoised energy intensity value, summed, and then divided by the total energy value of the current spatial energy projection matrix to obtain the longitude coordinates of the weighted energy geometric centroid. The latitude coordinates of each grid element are multiplied by the denoised energy intensity value, summed, and then divided by the total energy value of the current spatial energy projection matrix to obtain the latitude coordinates of the weighted energy geometric centroid. A spatial vector subtraction operation is performed on the two-dimensional coordinates of the weighted energy geometric centroid corresponding to the previous time slice and the current time slice to obtain the spatial coordinate offset.

[0020] The process of generating the spatiotemporal pattern migration vector and its output parameters involves using the weighted energy geometric centroid of the current time slice as the vector endpoint and the weighted energy geometric centroid of the previous time slice as the vector starting point to generate the spatiotemporal pattern migration vector. The starting and ending latitude and longitude coordinates of the spatiotemporal pattern migration vector are extracted. Reference ellipsoid parameters are introduced, and the great circle distance along the Earth's surface between the starting and ending points is calculated. This great circle distance is used as the absolute physical movement distance between consecutive time slices. Dividing the absolute physical movement distance by the time step span between consecutive time slices yields the absolute migration rate parameter.

[0021] The clockwise deflection angle of the spatiotemporal pattern migration vector relative to true north is extracted as the evolution direction parameter. The absolute migration rate parameter and the evolution direction parameter are output for the scheduling of external physical monitoring equipment. The principle of this step is to transform the high-dimensional spatial energy distribution state into low-dimensional kinematic geometric parameters, quantify the migration path of deformation energy in space, and quantitatively describe the macroscopic evolution trend of the disaster body.

[0022] This invention provides a method for quantifying the spatiotemporal patterns of landslide disasters on a wide-area scale. It has the following beneficial effects: 1. This invention calculates the effective physical seepage cycle by acquiring soil permeability and porosity, and then accumulates this data over a rainfall time series. It combines this with a deformation sequence obtained through segmented three-stage interpolation resampling to calculate the product of the maximum short-term cross-correlation coefficient and the mean absolute value of deformation. This operation eliminates the differences in data acquisition frequencies between different types of physical sensors and uses the mean absolute value of deformation rate to filter out background static areas with weak deformation, thus preserving signals from geologically active areas that exhibit genuine deformation and are strongly correlated with rainfall.

[0023] 2. This invention maps the nonlinear coupling index of rainfall and deformation to a global virtual space base map to generate a spatiotemporal third-order tensor, and performs high-order singular value decomposition to extract the truncated core tensor for inverse reconstruction. It uses a unified virtual space coordinate system to solve the data misalignment problem caused by inconsistent grid resolution in wide-area computing. At the same time, it uses tensor dimensionality reduction operation to separate and remove local random noise in the wide-area space and extract the spatial energy projection matrix that reflects the main geological deformation characteristics.

[0024] 3. This invention triggers the splitting and merging of quadtree grids based on the comparison results of the energy density value of the computational grid with a set threshold, and calculates the spatial coordinate offset of the weighted energy geometric centroid to generate a spatiotemporal pattern migration vector. It dynamically allocates the high-precision computational grid to high-energy geological deformation areas to control the overall computational load, and converts the high-dimensional spatial energy distribution state into specific migration rate and direction parameters, directly providing a quantitative control basis for the scheduling of external physical monitoring equipment. Attached Figure Description

[0025] Figure 1 This is a schematic diagram of the overall system architecture of an embodiment of the present invention; Figure 2 This is a flowchart of the macroscopic spatiotemporal quantization closed-loop processing according to an embodiment of the present invention; Figure 3 This is a comparison chart of experimental results from an embodiment of the present invention. Detailed Implementation

[0026] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0027] See attached document Figure 1 This invention provides a system for quantifying the spatiotemporal patterns of landslide disasters on a wide-area scale. The system operates in a computer equipment and network environment and includes: a spatiotemporal mapping and alignment module, a feature coupling and tensor purification module, a grid adaptive control module, and a parameter calculation module.

[0028] The system is connected to physical sensing devices, including a synthetic aperture radar (SAR) satellite system and a meteorological monitoring network. The SAR satellite system acquires surface physical echo phase data. The meteorological monitoring network includes regional weather radars and ground-based rain gauge nodes, acquiring hydrological and meteorological data. The physical sensing devices transmit data to the system via communication interfaces. The computer equipment on which the system operates includes a graphics processing unit (GPU) and a central processing unit (CPU), and is configured with a spatial database and a time-series database. The spatial database stores spatial metadata of the global virtual base map and the index relationships of the physical computing grid. The time-series database stores the time-series data collected by the physical sensing devices.

[0029] The spatiotemporal mapping and alignment module is configured to establish a global virtual space base map and generate an initial physical calculation grid at a set resolution. The module acquires the land surface deformation time series and precipitation time series within a time sliding window.

[0030] The feature coupling and tensor purification module is communicatively connected to the spatiotemporal mapping and alignment module. The feature coupling and tensor purification module is configured to execute a temporal resampling algorithm to unify the sampling frequency of different data sources, calculate the short-term cross-correlation coefficient between rainfall time series and surface deformation time series within the physical computing grid, and generate a nonlinear coupling index between rainfall and deformation. In the spatial dimension, the feature coupling and tensor purification module performs numerical mapping operations, constructing a third-order spatiotemporal tensor based on the nonlinear coupling index between rainfall and deformation. Finally, the module performs high-order singular value decomposition on the third-order spatiotemporal tensor, extracting a set number of principal components and inversely generating a spatial energy projection matrix.

[0031] The mesh adaptive control module is communicatively connected to the feature coupling and tensor purification module. The mesh adaptive control module is configured to receive the spatial energy projection matrix and compare the energy values ​​in the matrix with a set energy threshold. Based on the comparison result, the mesh adaptive control module generates mesh splitting or merging instructions. The module then sends these instructions to the spatiotemporal mapping and alignment module to update the physical computation mesh resolution for the next time sliding window.

[0032] The parametric calculation module communicates with the mesh adaptive control module. The parametric calculation module is configured to obtain the spatial energy projection matrix corresponding to two consecutive time slices. Based on spatial coordinates and energy values, the module calculates the energy geometric barycentric coordinates and generates a spatiotemporal pattern migration vector based on the difference in energy geometric barycentric coordinates between adjacent time slices. The module outputs the magnitude and argument of the spatiotemporal pattern migration vector.

[0033] See attached document Figure 2 This invention provides a method for quantifying the spatiotemporal patterns of landslide disasters on a wide-area scale, comprising the following steps: S10: Set the global virtual space base map, divide the physical calculation grid according to the initial resolution, and connect the surface deformation time series and rainfall time series within the current time sliding window; S20, within each physical computing grid, perform time dimension resampling operations on the surface deformation time series and rainfall time series, calculate the maximum short-term cross-correlation coefficient of the resampled surface deformation time series and rainfall time series, and calculate the nonlinear coupling index between rainfall and deformation based on the maximum short-term cross-correlation coefficient and the mean absolute value of deformation rate. S30 numerically maps the nonlinear coupling index of rainfall and deformation of the physical computing grid to the global virtual space base map, generates a regular spatiotemporal third-order tensor, performs high-order singular value decomposition on the spatiotemporal third-order tensor, extracts principal component data, and reconstructs the spatial energy projection matrix corresponding to the current time slice. S40, determine the energy distribution within the spatial energy projection matrix according to the set energy threshold, generate a mesh splitting command for regions with energy values ​​greater than the upper threshold, generate a mesh merging command for regions with energy values ​​less than the lower threshold, feed back the mesh splitting command and mesh merging command to step S10, and update the physical calculation mesh structure corresponding to the next time sliding window. S50 extracts the spatial energy projection matrix under continuous time slices, calculates the spatial coordinate offset of the energy geometric centroid corresponding to each spatial energy projection matrix, generates the spatiotemporal pattern migration vector based on the spatial coordinate offset result, and outputs the velocity and direction parameters corresponding to the spatiotemporal pattern migration vector.

[0034] See attached document Figure 2 Step S10 specifically includes setting a global virtual space base map, dividing the physical computation grid, and integrating multi-source time series data. The principle of this step is to establish a unified and fixed underlying coordinate reference system to accommodate the subsequently dynamically changing computation grid structure and to extract the physical state of geological bodies and hydrological elements in time segments. Step S10 specifically includes the following sub-steps: S101, Construct a global virtual spatial base map. Acquire geographic boundary data of the wide-area geological region to be monitored. Based on the geographic boundary data of the wide-area geological region to be monitored, establish a two-dimensional spatial matrix in the set geographic coordinate system. Define the first spatial dimension of the two-dimensional spatial matrix as the longitude direction and the second spatial dimension as the latitude direction. Set the highest spatial calculation accuracy supported by the system, which is determined by the highest spatial resolution of the raw imagery provided by the Synthetic Aperture Radar satellite system. Use the number of longitude grids corresponding to the highest spatial calculation accuracy as the maximum calculation resolution index of the first spatial dimension, denoted as . The number of dimensional grids corresponding to the highest spatial computational precision is used as the index of the maximum computational resolution of the second spatial dimension, denoted as . The resulting global virtual space base map can be represented by a two-dimensional grid coordinate set, defined by the following formula: ; In the formula, A set of two-dimensional grid coordinates representing the global virtual space base map; The coordinate index variable representing the first spatial dimension; The coordinate index variable representing the second spatial dimension; Indicates the maximum computational resolution index of the first spatial dimension; This represents the maximum computational resolution index for the second spatial dimension. The total number of grids contained in the global virtual space base map is... and The product of the coordinates. For the acquisition of geographic boundary data and the projection transformation of geographic coordinate systems, those skilled in the art can use conventional geographic information system software tools for processing. The coordinate mapping and projection calculations are well-known techniques in this field and will not be elaborated upon here.

[0035] S102, Initial Physical Calculation Grid Delineation. During system initialization, the global virtual spatial base map is aggregated using a preset initial spatial resolution to generate an initial physical calculation grid with spatial coordinate indices. To meet the symmetric splitting and merging requirements of the subsequent quadtree algorithm, the number of global virtual spatial base map grids included in the initial physical calculation grid in both the first and second spatial dimensions is set to an integer power of 2. The initial spatial resolution is determined based on the total area of ​​the wide-area geological region to be monitored and the memory capacity of the computer equipment. The physical size corresponding to the initial spatial resolution is larger than the physical size corresponding to the highest spatial calculation precision. The grid aggregation operation is completed by mapping the coordinates of multiple adjacent grids in the global virtual spatial base map to the same physical calculation grid. The system records and stores the longitude and latitude index sets corresponding to each physical calculation grid. The physical calculation grid constitutes the smallest spatial unit for subsequent time series data extraction and algebraic operations.

[0036] S103, Integrate heterogeneous time series data. Establish a time-sliding window mechanism. The time-sliding window has a fixed time span and slides continuously on the time axis according to a set time step. The fixed time span is determined based on the typical meteorological cycle of the wide-area geological region to be monitored, for example, it is set to include a complete rainfall cycle. The set time step is determined based on the revisit period of the synthetic aperture radar satellite. Within the current time-sliding window, the system retrieves and retrieves the surface deformation time series and rainfall time series from the external spatial database based on the spatial coordinate index of each physical computing grid.

[0037] The surface deformation time series is a set of discrete data points on the line-of-sight deformation rate obtained based on temporal interferometric synthetic aperture radar (TIAR) technology, with data units of mm / year. The precipitation time series is a set of precipitation observations collected by weather radar or rain gauges, with data units of mm. Due to the physical revisit period characteristics of TAR satellites, the sampling frequency of the surface deformation time series is lower than that of the precipitation time series. Both the surface deformation and precipitation time series acquired by the system carry absolute timestamps and are bound and stored with the corresponding physical computation grid spatial index. For phase unwrapping and deformation calculation of the temporal interferometric synthetic aperture radar data, those skilled in the art can use conventional interferometric measurement algorithms. The radar echo signal processing flow is a well-known technology in this field and will not be described in detail here.

[0038] See attached document Figure 2 Step S20 specifically includes performing time-dimension resampling operations within each physical computing grid and calculating the nonlinear coupling index between rainfall and deformation. These operations eliminate differences in data acquisition frequencies among different types of physical sensors and quantify the characteristics of rainfall-induced surface deformation. Step S20 specifically includes the following sub-steps: S201, Generate the cumulative rainfall sequence. Set the discrete time step, which is determined based on the fundamental time unit for effective physical seepage in the geological body, specifically set to 1 day. Apply a sliding window accumulation algorithm to the rainfall time series. Specifically, obtain empirical parameters of the soil and rock mass in the wide-area geological region to be monitored, including average soil permeability and porosity. Calculate the time required for rainfall to infiltrate to the estimated sliding surface based on the average soil permeability and porosity, and set this time as the effective physical seepage period. In practical engineering, this effective physical seepage period is typically taken as 7 to 15 discrete time steps. Use the effective physical seepage period as the accumulation window length, continuously slide it across the rainfall time series with discrete time steps, calculate the accumulated rainfall value within the window, and generate a cumulative rainfall sequence matching the discrete time step.

[0039] S202, Generate Daily Deformation Rate Sequence. For the surface deformation time series, a conformal piecewise cubic interpolation algorithm is used for continuous resampling in the time domain. Specifically, a piecewise cubic polynomial is constructed between two adjacent low-frequency discrete surface deformation data points, using the discrete time step as the resampling interval. The polynomial nodes are constrained to maintain monotonicity and continuity of the first derivative. This resampling process resamples the surface deformation time series into a daily deformation rate sequence with the same data length as the cumulative rainfall sequence. This resampling operation eliminates phase artifacts caused by sampling discontinuities in the interferometric synthetic aperture radar data. The specific polynomial coefficients of the conformal piecewise cubic interpolation algorithm can be solved using conventional numerical analysis methods by those skilled in the art; the mathematical solution process is well-known in the field and will not be elaborated here.

[0040] S203, Calculate the maximum short-term cross-correlation coefficient. Perform short-term cross-correlation analysis on the resampled cumulative rainfall sequence and daily deformation rate sequence. Within the current time sliding window, introduce a time shift variable to relatively shift the cumulative rainfall sequence and daily deformation rate sequence. Set the value range of the time shift variable to zero to the total duration corresponding to the effective physical infiltration period, and calculate the cross-correlation coefficient between the two sequences under different time shift variables. Extract the maximum value among all calculation results as the maximum short-term cross-correlation coefficient. The system synchronously records the optimal lag time corresponding to the generation of this maximum short-term cross-correlation coefficient. The optimal lag time characterizes the physical delay response period of rainfall infiltration into the ground and the resulting surface deformation.

[0041] S204, Calculate the nonlinear coupling index between rainfall and deformation. The nonlinear coupling index is calculated based on the maximum short-term cross-correlation coefficient and the mean absolute value of the deformation rate. Specifically, the absolute values ​​of all values ​​in the daily deformation rate sequence within the current time sliding window are extracted, and the arithmetic mean of all absolute values ​​is calculated to generate the mean absolute value of the deformation rate. To ensure the non-negative physical property of the subsequent energy tensor features, the absolute value of the maximum short-term cross-correlation coefficient is extracted, and the absolute value of the maximum short-term cross-correlation coefficient is multiplied by the mean absolute value of the deformation rate. The calculation formula is as follows: ; In the formula, This represents the nonlinear coupling index between rainfall and deformation in the current physical computing grid; This represents the absolute value of the maximum short-term cross-correlation coefficient; This represents the mean absolute value of the deformation rate. The product operation uses the mean absolute value of the deformation rate as a weight to filter out background static areas with high correlation coefficients but weak actual deformation levels, retaining geologically active areas with significant deformation and strong correlation to rainfall. The calculated nonlinear coupling index between rainfall and deformation is bound and stored with the spatial coordinate index of the corresponding physical computation grid.

[0042] See attached document Figure 2 Step S30 specifically includes numerically mapping the nonlinear coupling index of rainfall and deformation from the physical computation grid to the global virtual space base map, and performing high-order tensor mode purification. The principle of this step is to solve the data misalignment problem caused by inconsistent grid resolution, and to remove local random noise signals in the wide-area space through tensor dimensionality reduction operations, retaining surface deformation data with co-evolutionary characteristics. Step S30 specifically includes the following sub-steps: S301, Construct a standard spatiotemporal third-order tensor storage space. Establish a standard spatiotemporal third-order tensor storage space with the maximum computational resolution index of the first and second spatial dimensions of the global virtual space base map as the first two spatial dimensions, and the set historical sequence buffer window depth as the third dimension. Set the historical sequence buffer window depth as... , The value of this parameter is determined jointly by the computer memory threshold and the preset backtracking time span. During system operation, this value is maintained at a constant value. A sliding tensor queue. This standard spacetime third-order tensor has a strictly multidimensional orthogonal grid structure, used to carry subsequent higher-order algebraic operations.

[0043] S302, Perform heterogeneous grid numerical mapping. Obtain the nonlinear coupling index of rainfall and deformation for each physical computing grid within the current time sliding window. Since the physical computing grid will undergo adaptive reconstruction in subsequent steps, its spatial resolution differs from the highest spatial computation accuracy of the global virtual space map. To eliminate the heterogeneous data structure problem caused by this scale difference, a numerical mapping operation is performed. Specifically, the latitude and longitude coordinates of the global virtual space map covered by the physical computing grid are extracted. The nonlinear coupling index of rainfall and deformation of this physical computing grid is copied and filled into the coordinates of all corresponding virtual space map grids it covers. This fill operation ensures a rigorous mapping of the nonlinear coupling characteristics of spatial rainfall and deformation in a regular high-dimensional matrix. After traversing all physical computing grids and completing the filling, a two-dimensional data slice with a regular internal data structure is generated. This two-dimensional data slice is then pushed into the final time index layer of the standard spatiotemporal third-order tensor storage space to construct the updated spatiotemporal third-order tensor.

[0044] S303, Higher-Order Singular Value Decomposition and Core Tensor Generation. Higher-order singular value decomposition is performed on the generated spatiotemporal third-order tensor. Higher-order singular value decomposition transforms the original spatiotemporal third-order tensor into a core tensor and the modular product of three orthogonal factor matrices corresponding to the longitude, latitude, and time dimensions, respectively. The core tensor reflects the energy distribution of the original multidimensional data along each eigenvector direction, and the orthogonal factor matrices contain the spatial and temporal basis features. In geological disaster monitoring, wide-area environmental noise typically exhibits a low-energy discrete distribution, while actual landslide deformation exhibits a high-energy clustered distribution. For the specific matrix expansion and singular value solving process of higher-order singular value decomposition, those skilled in the art can use the conventional alternating least squares method, and its multilinear algebraic calculation steps are well-known techniques in the field, and will not be elaborated here.

[0045] S304, Energy Truncation and Spatial Energy Projection Matrix Reconstruction. Calculate the Frobenius norm of each principal component term in the core tensor, and use the square of this Frobenius norm as the energy value of the corresponding principal component term. Sort the principal component terms of the core tensor in descending order according to their energy values. Before calculation... The cumulative energy percentage of each principal component. An energy retention threshold is set; in practical engineering, this threshold is typically set to 85% to 95%. When the cumulative energy percentage reaches this energy retention threshold, the previous energy is extracted. The principal component terms constitute the truncated core tensor. The truncation operation removes low-energy random noise from the wide-area background. Using the truncated core tensor and the orthogonal factor matrices corresponding to the three dimensions, an inverse modular multiplication operation is performed to reconstruct the noise-removed spatiotemporal third-order tensor, and the two-dimensional data surface corresponding to the current time slice is extracted as the spatial energy projection matrix. The calculation formula is as follows: ; In the formula, This represents the spatial energy projection matrix corresponding to the current time slice; Indicates before extraction The truncated core tensor after each principal component term; This represents the orthogonal factor matrix corresponding to the first spatial dimension; This represents the orthogonal factor matrix corresponding to the second spatial dimension; This represents the orthogonal factor matrix corresponding to the time dimension; , and These represent the modular multiplication operations of the tensor along the first, second, and third dimensions, respectively; subscripts This represents the latest slice index value of the current time sliding window in the third dimension, i.e. .

[0046] Because tensor dimensionality reduction and reconstruction operations are prone to Gibbs oscillations at numerical margins, the system projects the reconstructed spatial energy matrix... All elements in the matrix undergo non-negative thresholding, which involves iterating through all element values ​​in the matrix, setting elements less than zero to zero, and retaining elements greater than or equal to zero. After non-negative thresholding, the element values ​​in the generated spatial energy projection matrix represent the denoised geological deformation energy intensity in the corresponding spatial coordinate system, providing a quantitative benchmark for subsequent grid feedback control.

[0047] See attached document Figure 2Step S40 specifically includes determining the energy distribution within the spatial energy projection matrix based on a set energy threshold and generating grid closed-loop feedback control commands. The principle behind this step is to establish a dynamic allocation mechanism for computing resources. In wide-area geological disaster monitoring, most areas belong to stable background zones, with only a small portion exhibiting deformation evolution. Focusing a high-precision computational grid on active geological deformation zones with high-energy characteristics to capture deformation details, while simultaneously reducing the grid resolution of low-energy stable background zones to release computing power, achieves a balanced configuration of system performance and spatial resolution. Step S40 specifically includes the following sub-steps: S401: Calculate the energy density value of the physical computing grid and dynamically set the energy threshold. Obtain the spatial energy projection matrix corresponding to the current time slice. Extract the latitude and longitude coordinates of the global virtual spatial base map contained in each physical computing grid. Sum the energy values ​​corresponding to these coordinates in the spatial energy projection matrix and divide by the number of virtual spatial base map grids contained in that physical computing grid to calculate the energy density value of each physical computing grid. Area normalization can eliminate the energy accumulation error caused by the area difference between different grid levels. Calculate the expected value and standard deviation of the energy density values ​​of all current physical computing grids. Construct a dynamic upper and lower energy threshold using the expected value and standard deviation, calculated as follows: ; ; In the formula, Indicates the upper limit threshold of energy; Indicates the lower energy threshold; This represents the mathematical expectation of the energy density values ​​of all physical computation grids; This represents the standard deviation of the energy density values ​​for all physical computation grids. This represents the upper limit control coefficient, which is set to 1.5 to 2.0 in the actual engineering range; This represents the lower limit control coefficient, which is set to 0.5 to 1.0 in actual engineering applications. By introducing statistical parameters, the system can adaptively adjust the threshold boundary according to the overall energy distribution characteristics of different geological regions.

[0048] S402 generates a mesh splitting command based on a quadtree structure. It iterates through all physical computation meshes within the current time sliding window, comparing their energy density values ​​with the energy upper limit threshold. When the energy density value of a physical computation mesh is determined to be greater than the energy upper limit threshold, it is determined that the mesh region exhibits significant rainfall-induced deformation characteristics. The system triggers a mesh splitting command for this physical computation mesh. The mesh splitting command, according to the quadtree spatial indexing rules, divides the current physical computation mesh into four equal-area sub-meshes along the longitude and latitude centerlines. Before performing the division operation, the system compares the resolution of the current physical computation mesh with the highest spatial computation accuracy of the global virtual space base map. If the current mesh resolution has reached the limit physical size set by the highest spatial computation accuracy, the splitting is terminated to prevent coordinate index overflow.

[0049] S403, generate a mesh merging command based on a quadtree structure. During the traversal and comparison process, the energy density value of the physical computation mesh is compared with the lower energy threshold. Based on the node hierarchy of the quadtree, find four adjacent physical computation meshes of the same level belonging to the same parent node. When it is determined that the energy density values ​​of these four adjacent physical computation meshes are all less than the lower energy threshold, the local area is determined to be in a relatively stable geological state. The system triggers a mesh merging command for these four physical computation meshes. The mesh merging command merges the above four sub-meshes back into a parent mesh of the next higher level. Before executing the merging operation, the system verifies the size of the merged mesh to ensure that it is not greater than the physical size corresponding to the initial spatial resolution set in step S10. For the node splitting, merging, and hierarchical traversal of the quadtree data structure, those skilled in the art can use conventional spatial database indexing algorithms for processing. Its basic data structure operations are well-known technologies in the field and will not be described in detail here.

[0050] S404 executes a closed-loop update of the physical grid structure. It collects all grid splitting and merging instructions generated in the current time slice. These instructions are fed back to the system's spatial database, overwriting and updating the stored longitude and latitude index sets of the physical computation grid. The system then proceeds to the next time sliding window's computation flow according to the updated spatial coordinate index relationships. This operation implements a closed-loop feedback mechanism where the front-end data sampling grid is inversely controlled by the back-end algebraic solution energy characteristics, ensuring that the deformation and precipitation data extraction in the next time slice matches the spatial scale of geological hazard evolution.

[0051] See attached document Figure 2Step S50 specifically includes extracting the spatial energy projection matrix under continuous time slices, calculating the spatial coordinate offset of the energy geometric centroid corresponding to each spatial energy projection matrix, and generating spatiotemporal pattern migration vector output parameters. The principle of this step is to transform the high-dimensional abstract spatial energy distribution state into low-dimensional kinematic geometric parameters with intuitive physical meaning. In a wide-area geological environment, the deformation of local landslide bodies often exhibits spatial correlation and expansion. By introducing the concept of the center of mass from physics, and replacing mass with deformation energy, the transfer trajectory of the core deformation energy in space can be tracked, which can quantitatively describe the macroscopic spread and evolution trend of wide-area landslide disaster groups, providing clear geometric indicators for the scheduling and deployment of on-site physical monitoring equipment. Step S50 specifically includes the following sub-steps: S501, Extract continuous time slice data. Obtain the spatial energy projection matrix after purification and reconstruction in step S30 and under grid closed-loop control in step S40. Extract two adjacent time slices from the system's spatial database, denoted as the previous time slice and the current time slice, respectively. Read the spatial energy projection matrices corresponding to the previous time slice and the current time slice, as well as the latitude and longitude coordinates set of the global virtual spatial base map corresponding to each grid element in the matrix.

[0052] S502, Calculate the two-dimensional coordinates of the weighted energy geometric centroid. Using the actual geographic latitude and longitude coordinates of each grid element in the spatial energy projection matrix as the base position, and the denoised energy intensity value of the corresponding grid as the weight, calculate the two-dimensional coordinates of the weighted energy geometric centroid for the previous and current time slices respectively. Specifically, multiply the actual geographic longitude coordinates of each grid element by the denoised energy intensity value, sum the results, and then divide by the total energy value of the current spatial energy projection matrix to obtain the longitude coordinates of the centroid; similarly, calculate the latitude coordinates of the centroid. This centroid is used to characterize the core location of geological deformation energy concentration within the region, and its calculation formula is as follows: ; ; In the formula, The longitude coordinates representing the weighted energy geometric centroid; The latitudinal coordinates representing the weighted energy geometric centroid; This represents the total number of grid elements contained in the spatial energy projection matrix. The value of is equal to the product of the maximum number of grids in the longitude direction and the maximum number of grids in the latitude direction of the global virtual space base map; Indicates the first The actual geographic longitude coordinates corresponding to the center point of each grid element; Indicates the first The actual geographic latitude coordinates corresponding to the center point of each grid element; Indicates the first The denoised energy intensity value of each grid element. Before performing the above calculation, the system will determine the total energy value of the current spatial energy projection matrix (i.e., If the value is zero or less than the set machine minimum, it indicates that the current wide-area geological region is in a stable dormant state. The system will automatically terminate the centroid calculation of the current time slice and directly assign the absolute migration rate parameter of the current time slice to zero. If the value is false, the system will use the above formula to obtain the two-dimensional coordinates of the weighted energy geometric centroid of the previous time slice and the two-dimensional coordinates of the weighted energy geometric centroid of the current time slice.

[0053] S503, Generate Spatiotemporal Pattern Migration Vector. The system determines whether the weighted energy geometric centroid coordinates of both the previous and current time slices have been successfully calculated. If either time slice fails to calculate its centroid coordinates due to being in a stable dormant state, the system directly sets the current spatiotemporal pattern migration vector to zero. If both time slices have valid centroid coordinates, the system performs a spatial vector subtraction operation on the weighted energy geometric centroid coordinates of the previous and current time slices. Using the weighted energy geometric centroid coordinates of the current time slice as the vector endpoint and the weighted energy geometric centroid coordinates of the previous time slice as the vector starting point, a spatiotemporal pattern migration vector representing the macroscopic spread trajectory is generated. This spatiotemporal pattern migration vector indicates the specific spatial offset path of the geological disaster group deformation energy transfer in the geographic coordinate system.

[0054] S504 calculates parameters and outputs them to the hardware scheduling interface. Based on the generated spatiotemporal pattern migration vector, it calculates its physical movement distance and geometric argument. The system first determines whether the spatiotemporal pattern migration vector is a zero vector; if it is, the absolute physical movement distance, absolute migration rate parameter, and evolution direction parameter are all set to zero, skipping subsequent geometric calculations; if it is a non-zero vector, the system extracts the latitude and longitude coordinates of the starting and ending points of the spatiotemporal pattern migration vector, introduces the WGS84 reference ellipsoid parameters, calculates the great circle distance along the Earth's surface between the starting and ending points, and uses this great circle distance as the absolute physical movement distance of the energy center of gravity between continuous time slices, with the unit being meters. Dividing this absolute physical movement distance by the time step span between continuous time slices yields the absolute migration rate parameter.

[0055] The geometric argument represents the clockwise deflection angle of the spatiotemporal pattern migration vector relative to true north, and is used as an evolution direction parameter. The system outputs the absolute migration rate parameter and the evolution direction parameter to the external physical monitoring equipment scheduling interface. The scheduling interface adjusts the preset flight path of the UAV inspection equipment according to the evolution direction parameter, and densifies the deployment of ground-based global navigation satellite system monitoring base stations in the target area according to the absolute migration rate parameter. For the specific calculation of the great circle distance and azimuth between two points in the geographic coordinate system, those skilled in the art can use conventional geodetic formulas for processing, and the geometric calculation is a well-known technology in the field, which will not be elaborated here.

[0056] Specific application examples: Application Scenario Setting: A mountainous area in western China prone to heavy rainfall (approximately 20km × 20km) was selected as the monitoring area. Five consecutive time-sliding windows were monitored, each spanning 30 days with a step size of 12 days (consistent with the revisit period of Synthetic Aperture Radar satellites), for a total of 60 days. The method provided in this invention (mesh adaptation and tensor purification mechanism) was compared and verified with the traditional fixed-mesh static threshold method.

[0057] Specific implementation and chart data mapping: Step S10: Global Base Map Construction and Physical Calculation Mesh Initialization Acquire SAR images of the region with a maximum spatial resolution of 10m, and construct an image with a size of [size missing]. , (Total 4×10) 6 A global virtual space base map (number of grids). In the first time sliding window of system initialization, the physical computation grid is divided with an initial spatial resolution of 160m.

[0058] Chart mapping (resource consumption characteristics): as attached Figure 3 As shown in (a), in the initial stage (time sliding window 1), both the traditional method and the method of the present invention need to perform equal-resolution operations on the global matrix, so the initial memory consumption of both is about 850MB.

[0059] Step S20: Multi-source time series access and coupling index calculation Within a specific physical computation grid, empirical permeability parameters are extracted for the region. An effective physical seepage cycle of 7 days is set, and a cumulative rainfall sequence is generated. For the line-of-sight deformation discrete points calculated by InSAR (assuming the average annualized creep rate of this local grid is...),... (mm / year), and daily deformation rate sequences were generated through conformal piecewise cubic interpolation. After translation calculation, the absolute value of the maximum short-term cross-correlation coefficient between rainfall and deformation for this grid was... The nonlinear coupling index between rainfall and deformation, calculated using the formula, is 12.75.

[0060] Step S30: Construction of the third-order spacetime tensor and higher-order singular value decomposition The nonlinear coupling exponents of each grid are mapped to a global virtual space base map, and a spatiotemporal third-order tensor is constructed by combining a historical buffer window (depth set to 5). Higher-order singular value decomposition is performed, the energy retention threshold is set to 90%, low-frequency environmental random noise (such as vegetation weathering and atmospheric phase delay) is truncated, and the spatial energy projection matrix is ​​reconstructed.

[0061] Chart mapping (false alarm rate characteristics): as attached Figure 3 As shown in (b), due to the lack of tensor purification mechanism in traditional methods, it is impossible to remove background false positive signals with high correlation but weak deformation level. Its false alarm rate is consistently distributed between 15.2% and 17.1% in 5 windows. The method of the present invention removes environmental noise by truncating the low-energy principal components in the core tensor, and stably controls the false alarm rate in the low range of 3.7% to 4.1%.

[0062] Step S40: Adaptive Energy Threshold Determination and Mesh Reconstruction Based on the current spatial energy projection matrix, calculate the mathematical expectation of the energy density values ​​of all physical computation grids. The standard deviation is The system dynamically generates an upper energy threshold of 4.25 and a lower energy threshold of 1.25. For the aforementioned deformation-active region with a coupling index of 12.75, the system triggers a mesh splitting command to refine its physical calculation mesh to 80m (and can continue to subdivide it in subsequent slices until the set maximum spatial calculation accuracy of 10m is reached); for a large number of stable dormant regions with energy densities below 1.25, a mesh merging command is triggered.

[0063] Chart mapping (computing power release characteristics): as attached Figure 3 As shown in time sliding windows 2 to 5 in (a), due to the activation of mesh merging and splitting instructions, the method of this invention reduces the number of computation nodes in the stable background region. Compared with the traditional method, which always maintains a high memory consumption of over 850MB, the memory usage of the method of this invention decreases stepwise as the time sliding window progresses, and eventually stabilizes at about 320MB, releasing about 62% of the computing resources.

[0064] Step S50: Spatiotemporal pattern migration vector calculation and terminal scheduling The spatial energy projection matrix of continuous time slices is extracted, and the two-dimensional coordinates of the weighted energy geometric centroid are calculated. The calculation shows that the energy centroid of the current time slice has moved 315 meters compared to the previous time slice, corresponding to an absolute migration rate of 26.25 m / day, and a geometric argument of [missing information]. (Southeast direction).

[0065] Here, the absolute physical movement distance of 315 meters represents the spatial macroscopic spread trajectory of the area where the deformation energy of the geological disaster group is concentrated, rather than the actual sliding displacement of the strata themselves. The actual surface deformation monitored by InSAR remains at the millimeter level, but this parameter indicates that the active area of ​​deformation has expanded 315 meters southeastward within 12 days. Based on the migration vector parameter of this spatiotemporal pattern (velocity 26.25 m / day, direction 135°), the system guides the UAV inspection equipment to prioritize coverage of the southeastern area of ​​the target through the scheduling interface, and densifies the deployment of ground GNSS monitoring base stations.

Claims

1. A method for quantifying the spatiotemporal patterns of landslide disasters on a wide-area scale, characterized in that, include: Set a global virtual space base map, divide the physical calculation grid according to the initial resolution, and connect the surface deformation time series and rainfall time series within the current time sliding window; Within the physical computing grid, a time dimension resampling operation is performed on the surface deformation time series and the rainfall time series to calculate the maximum short-term cross-correlation coefficient. Based on the maximum short-term cross-correlation coefficient and the mean absolute value of the deformation rate, the nonlinear coupling index between rainfall and deformation is calculated. The nonlinear coupling index of rainfall and deformation is numerically mapped to the global virtual space base map to generate a third-order spatiotemporal tensor. The third-order spatiotemporal tensor is then subjected to high-order singular value decomposition to reconstruct and calculate the spatial energy projection matrix. The energy distribution within the spatial energy projection matrix is ​​determined based on the set energy threshold, and a mesh splitting command or mesh merging command is generated to update the physical calculation mesh structure corresponding to the next time sliding window. Extract the spatial energy projection matrix under continuous time slices, calculate the spatial coordinate offset of the energy geometric centroid, generate a spatiotemporal pattern migration vector, output the velocity and direction parameters corresponding to the spatiotemporal pattern migration vector, and schedule external physical monitoring equipment based on the velocity and direction parameters.

2. The method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale according to claim 1, characterized in that, The process of setting a global virtual space base map and dividing the physical computing grid according to the initial resolution includes: Obtain the geographic boundary data of the wide geological area to be monitored, establish a two-dimensional spatial matrix under the set geographic coordinate system, and define the first spatial dimension of the two-dimensional spatial matrix as the longitude direction and the second spatial dimension as the latitude direction; The number of longitude grids corresponding to the highest spatial calculation precision is used as the maximum calculation resolution index of the first spatial dimension, and the number of latitude grids corresponding to the highest spatial calculation precision is used as the maximum calculation resolution index of the second spatial dimension, thereby generating a global virtual space base map represented by a two-dimensional grid coordinate set; During the system initialization phase, the global virtual space base map is meshed according to the preset initial spatial resolution, and the coordinates of multiple adjacent grids in the global virtual space base map are mapped to the same physical computing grid to generate an initial physical computing grid with spatial coordinate index.

3. The method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale according to claim 1, characterized in that, Performing time-dimensional resampling operations on the surface deformation time series and rainfall time series within the physical computing grid includes: The average soil permeability and porosity are obtained. Based on the average soil permeability and porosity, the effective physical seepage period required for rainfall to infiltrate to the estimated sliding surface is calculated. The effective physical seepage period is used as the cumulative window length and continuously slides on the rainfall time series. The cumulative rainfall value within the window is calculated to generate a cumulative rainfall series. For the aforementioned surface deformation time series, a shape-preserving piecewise cubic interpolation algorithm is used for continuous resampling in the time domain. A piecewise cubic polynomial is constructed between two adjacent discrete surface deformation data points, and the surface deformation time series is resampled into a daily deformation rate series with the same data length as the cumulative rainfall series.

4. The method for quantifying the spatiotemporal pattern of landslide disasters at a wide-area scale according to claim 3, characterized in that, The calculation of the maximum short-term cross-correlation coefficient, and the calculation of the nonlinear coupling index between rainfall and deformation based on the maximum short-term cross-correlation coefficient and the mean absolute value of deformation rate, include: A time-shifting variable is introduced to perform relative shifting on the cumulative rainfall sequence and the daily deformation rate sequence. The maximum value of the cross-correlation coefficient between the cumulative rainfall sequence and the daily deformation rate sequence under different time-shifting variables is extracted as the maximum short-term cross-correlation coefficient. Extract the absolute values ​​of all values ​​in the daily deformation rate sequence within the current time sliding window, calculate the arithmetic mean of all absolute values, and generate the mean absolute value of deformation rate. Extract the absolute value of the maximum short-term cross-correlation coefficient, and multiply the absolute value of the maximum short-term cross-correlation coefficient by the mean of the absolute values ​​of the deformation rates to generate the nonlinear coupling index of rainfall and deformation for the current physical computing grid.

5. The method for quantifying the spatiotemporal pattern of landslide disasters at a wide-area scale according to claim 1, characterized in that, The nonlinear coupling index between rainfall and deformation is numerically mapped to a global virtual space base map to generate a spatiotemporal third-order tensor, including: Establish a standard spatiotemporal third-order tensor storage space with the maximum computational resolution index of the first spatial dimension and the maximum computational resolution index of the second spatial dimension as the first two spatial dimensions, and the set historical sequence buffer window depth as the third dimension. Extract the latitude and longitude coordinate set of the global virtual space base map covered by the physical computing grid, copy the nonlinear coupling index of rainfall and deformation of the physical computing grid equally and fill it into the coordinates of all corresponding virtual space base map grids covered by the physical computing grid, and generate two-dimensional data slices with regular internal data structure. The two-dimensional data slices are pushed into the last time index layer of the standard spatiotemporal third-order tensor storage space to construct a spatiotemporal third-order tensor.

6. The method for quantifying the spatiotemporal pattern of landslide disasters at a wide-area scale according to claim 5, characterized in that, High-order singular value decomposition is performed on the third-order spacetime tensor to reconstruct and calculate the spatial energy projection matrix, including: A high-order singular value decomposition operation is performed on the generated spatiotemporal third-order tensor to convert the spatiotemporal third-order tensor into a core tensor and a modular product of three orthogonal factor matrices corresponding to the longitude direction, latitude direction and time dimension, respectively. The square of the Frobenius norm of each principal component term in the core tensor is calculated as the energy value of the corresponding principal component term. The principal component terms of the core tensor are sorted in descending order according to the energy value. The top few principal component terms whose cumulative energy percentage reaches the energy retention threshold are extracted to form the truncated core tensor. Using the truncated core tensor and the orthogonal factor matrices corresponding to the three dimensions, an inverse modular multiplication operation is performed to reconstruct the noise-removed spatiotemporal third-order tensor, and the two-dimensional data surface corresponding to the current time slice is extracted as the spatial energy projection matrix. Iterate through and compare all element values ​​in the spatial energy projection matrix, forcibly set elements less than zero to zero, and perform non-negative thresholding on elements greater than or equal to zero.

7. The method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale according to claim 1, characterized in that, The set energy threshold is generated through the following steps: Extract the set of latitude and longitude coordinates of the global virtual space base map contained in each physical computing grid, sum the energy values ​​corresponding to these coordinates in the spatial energy projection matrix, and divide by the number of virtual space base map grids contained in the physical computing grid to calculate the energy density value of each physical computing grid. Calculate the expected value and standard deviation of the energy density values ​​for all current physical computing grids; The upper energy threshold is calculated using the expected value, standard deviation, and upper limit control coefficient. The lower energy threshold is calculated using the expected value, standard deviation, and lower limit control coefficient.

8. A method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale, as described in claim 7, is characterized in that... The generated mesh splitting or mesh merging instructions include: When it is determined that the energy density value of a certain physical computing grid is greater than the energy upper limit threshold, a grid splitting command is triggered for the physical computing grid. The grid splitting command divides the current physical computing grid into four subgrids of equal area along the longitude and latitude center lines according to the quadtree spatial indexing rules. Find four adjacent physical computing grids of the same level that belong to the same parent node. When it is determined that the energy density values ​​of these four adjacent physical computing grids are all less than the lower energy threshold, trigger a grid merging command for these four physical computing grids. The grid merging command merges these four sub-grids back into a parent grid of the next higher level.

9. A method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale, as described in claim 1, is characterized in that... The spatial coordinate offset of the calculated energy geometric centroid includes: Using the actual geographic longitude and actual geographic latitude coordinates of each grid element in the spatial energy projection matrix as the base position, and using the denoised energy intensity value of the corresponding grid element as the weight, the actual geographic longitude coordinates of each grid element are multiplied and summed, and then divided by the total energy value of the current spatial energy projection matrix to obtain the longitude coordinates of the weighted energy geometric centroid. Multiply the actual geographic latitude coordinates of each grid element by the denoised energy intensity value and sum them, then divide by the total energy value of the current spatial energy projection matrix to obtain the latitude coordinates of the weighted energy geometric centroid. Perform a spatial vector subtraction operation on the two-dimensional coordinates of the weighted energy geometric centroid corresponding to the previous time slice and the current time slice to obtain the spatial coordinate offset.

10. A method for quantifying the spatiotemporal patterns of landslide disasters at a wide-area scale, as described in claim 9, is characterized in that... The generation of the spatiotemporal pattern migration vector outputs the velocity and direction parameters corresponding to the spatiotemporal pattern migration vector, including: The two-dimensional coordinates of the weighted energy geometric centroid of the current time slice are used as the vector endpoint, and the two-dimensional coordinates of the weighted energy geometric centroid of the previous time slice are used as the vector starting point to generate a spatiotemporal pattern migration vector. Extract the latitude and longitude coordinates of the starting point and the ending point of the spatiotemporal pattern migration vector, introduce the reference ellipsoid parameters, calculate the great circle distance along the Earth's surface between the starting point and the ending point, and use the great circle distance as the absolute physical movement distance between continuous time slices. Divide the absolute physical distance traveled by the time step span between consecutive time slices to obtain the absolute migration rate parameter. The clockwise deflection angle of the spatiotemporal pattern migration vector relative to the due north direction is extracted as the evolution direction parameter. The absolute migration rate parameter and the evolution direction parameter are output for scheduling of external physical monitoring equipment.

Images