A method for predicting groundwater rise under heavy rainfall conditions

Abstract

This invention provides a method for predicting groundwater rise under heavy rainfall conditions, relating to the field of hydrogeology. The invention acquires rainfall intensity and groundwater level data, fuses them to generate a spatiotemporal matrix, scans the rainfall intensity data in the matrix to identify heavy rainfall events, and extracts the numerical features of the spatial distribution of rainfall for each event to construct a multidimensional feature vector. This vector is input into a pre-trained feedforward neural network. The model calculates the pressure field state tensor through its internal parameterized decay function and anisotropic diffusion convolution kernel. Trend decomposition of historical water level data yields the background water level digital field, which is superimposed with the pressure field state tensor to obtain the water level response prediction field. Based on the water level change curve of this prediction field, the groundwater level rise rate is output after differentiation, thereby improving the accuracy and timeliness of predicting groundwater dynamics induced by heavy rainfall.

Description

Technical Field

[0001] This invention relates to the field of hydrogeology, specifically to a method for predicting groundwater rise under heavy rainfall conditions. Background Technology

[0002] Heavy rainfall events are usually accompanied by a high-intensity water input in a short period of time. This water input can quickly infiltrate underground aquifers, causing a significant rise in groundwater levels. Accurately predicting the magnitude and rate of groundwater level rise at specific locations under heavy rainfall conditions is of great significance for early warning of geological disasters, safe operation and maintenance of underground engineering projects, and emergency management of water resources. However, the response of groundwater to rainfall is a complex process influenced by multiple coupled factors, including the spatiotemporal heterogeneity of rainfall, the anisotropy of underground media permeability, and the natural dynamic background of groundwater levels. This makes accurate prediction a huge challenge.

[0003] Currently, prediction methods in this field mainly rely on two types of models: one is the hydrogeological numerical model based on physical laws, and the other is the machine learning model based on the statistical relationship of historical data. The hydrogeological numerical model attempts to simulate the response process by solving the partial differential equation of groundwater movement. Its effectiveness is highly dependent on aquifer parameters that are difficult to obtain accurately, such as permeability coefficient, specific yield, and fine boundary conditions. The data-driven machine learning model focuses on learning the statistical relationship between rainfall and water level from historical monitoring data, but usually simplifies the input data to a single point or regional average sequence, failing to fully explore the inherent spatiotemporal structure information of rainfall events. Existing technologies suffer from the following shortcomings. First, at the data utilization level, current methods fail to achieve deep fusion and collaborative analysis of rainfall fields and groundwater level fields within a unified spatiotemporal framework. They often treat these as independent time series, neglecting the crucial impact of rainfall spatial distribution patterns on differentiated pressure excitation at different underground locations. This leads to incomplete input information for the prediction model. Second, at the feature engineering level, existing data-driven models lack effective, physically consistent feature characterization methods for heavy rainfall events. They fail to comprehensively encode the intensity and spatial structure of a rainfall event, such as heterogeneity, center of gravity location, and temporal cumulative effects, into a digital feature that can guide model inference. This results in the model's learning process remaining at a shallow correlation level, lacking comprehensive understanding. The existing methods suffer from poor analytical and physical interpretability. Furthermore, at the model construction level, neither purely physical models nor purely data models can effectively reconcile the relationship between mechanisms and data. Physical models become distorted due to oversimplification or inaccurate parameters, while data models may produce irrational predictions that violate diffusion laws due to a lack of physical constraints. In particular, they cannot simulate the anisotropic spatial diffusion process of pressure in underground media. Finally, at the prediction synthesis level, existing methods usually fail to clearly distinguish and separate the long-term natural trend background of groundwater levels from the short-term rainfall disturbance response, resulting in irrelevant trend noise mixed in the prediction results. It is difficult to clearly separate and accurately quantify the water level rise component directly caused by this heavy rainfall event, which is the key information most needed for early warning decision-making.

[0004] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] The purpose of this invention is to provide a method for predicting groundwater rise under heavy rainfall conditions, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: A method for predicting groundwater rise under heavy rainfall conditions, comprising the following steps: Step 1: Obtain rainfall intensity data and groundwater level data for the target area, and apply a spatial grid interpolation algorithm to convert the groundwater level data into regular grid data with the same spatial resolution as the rainfall intensity data, generating a spatiotemporal matrix containing rainfall intensity channels and groundwater level channels; Step 2: Continuously scan the rainfall intensity data in the spatiotemporal matrix, use the threshold detection method to identify continuous periods exceeding the set heavy rainfall threshold and record them as heavy rainfall events. For each heavy rainfall event, calculate the information entropy of its rainfall intensity distribution, determine the directional offset vector of its rainfall intensity weighted centroid relative to the preset prediction point, and calculate the integral cumulative value of rainfall intensity over time during the duration of the event. Combine the information entropy, directional offset vector and integral cumulative value to construct a multidimensional feature vector. Step 3: Input the multidimensional feature vector into the pre-trained feedforward neural network model. The model performs numerical processing based on the internal parameterized decay function and anisotropic diffusion convolution kernel to calculate the pressure field state tensor of the target region. Step 4: Perform trend decomposition on historical groundwater level data to obtain the background water level digital field. Perform tensor addition operation on the pressure field state tensor and the background water level digital field to obtain the water level response prediction field. Extract the water level change prediction curve of the preset prediction point from the water level response prediction field, and output the groundwater level rise rate after differentiation operation.

[0007] Furthermore, the groundwater level data is converted into regular grid data with the same spatial resolution as the rainfall intensity data. Specifically, the sampling interval between adjacent time sampling points is preset, and multiple groundwater level measurement points with different geographical coordinates are selected within the target area to obtain groundwater level data at each time sampling point during the current monitoring period. Simultaneously, rainfall intensity data covering the same target area within the same monitoring period is acquired. This rainfall intensity data has a unified regular geographical grid and a fixed spatial resolution, and each time sampling point has corresponding gridded rainfall intensity data. Based on the regular geographical grid used for the rainfall intensity data, for each time sampling point, the Kriging spatial interpolation algorithm is used to calculate and generate the corresponding groundwater level estimate for each spatial grid in the regular geographical grid, using the groundwater level data and its geographical coordinates at that time sampling point as input. Through the above operations, the originally discrete groundwater level data is converted into gridded groundwater level data that is aligned with the rainfall intensity data in the spatial grid framework and synchronized with it in time. A spatiotemporal matrix containing rainfall intensity channels and groundwater level channels is generated. Specifically, for the gridded rainfall intensity data and gridded groundwater level data corresponding to each time sampling point within the monitoring period, they are registered according to the same spatial grid coordinates. The two are superimposed in the data channel dimension to form a composite data layer with dual channels. The composite data layers corresponding to all time sampling points are arranged sequentially in chronological order to construct a three-dimensional spatiotemporal matrix. The first and second dimensions of this spatiotemporal matrix correspond to the row index and column index of the spatial grid, respectively, and the third dimension corresponds to the time sampling point index. The rainfall intensity data and the groundwater level data are stored as independent channel data in this spatiotemporal matrix.

[0008] Furthermore, a threshold detection method is used to identify continuous periods exceeding a set heavy rainfall threshold and record them as heavy rainfall events. Specifically, the rainfall intensity data in the spatiotemporal matrix is ​​scanned sequentially along the time dimension. For each time sampling point, the average rainfall intensity value of all spatial grids in the target area corresponding to that time sampling point is calculated. The average rainfall intensity value of each time sampling point is continuously compared with a preset heavy rainfall threshold to identify all time sampling points whose average rainfall intensity value exceeds the heavy rainfall threshold. Among the identified time sampling points, a series of time sampling points that are temporally continuous and whose time interval between two adjacent time sampling points does not exceed the preset maximum allowable interval are merged into an independent heavy rainfall event. At the same time, the start time and end time of each heavy rainfall event are recorded. For any heavy rainfall event, the information entropy of its rainfall intensity distribution is calculated as follows: the time sampling point when the rainfall intensity reaches the global maximum value during the duration of the heavy rainfall event is taken as the peak time; the rainfall intensity values ​​of all spatial grids in the target area at the peak time are extracted, and these rainfall intensity values ​​are normalized so that their sum is 1, thus forming a spatial probability distribution representing the spatial distribution of rainfall intensity; based on the Shannon entropy calculation formula in information theory, the information entropy of this spatial probability distribution is calculated, and the resulting entropy value is the information entropy index representing the heavy rainfall event. For any heavy rainfall event, the directional offset vector of its weighted centroid of rainfall intensity relative to the preset prediction point is determined. Specifically, the peak time of the heavy rainfall event is extracted, and the average rainfall intensity value of each spatial grid at the peak time is used as the weight. The weighted average values ​​of all spatial grids in the east-west and north-south directions are calculated respectively. The coordinate point formed by these two weighted average values ​​is the coordinate of the weighted centroid of rainfall intensity at the peak time of the heavy rainfall event. The plane vector pointing from the coordinate of the weighted centroid of rainfall intensity to the coordinate of the preset prediction point is calculated. The plane vector is normalized to make its magnitude 1, resulting in a unit vector that only represents directional information. The resulting unit vector is the directional offset vector.

[0009] Furthermore, for any heavy rainfall event, the cumulative integral of rainfall intensity over time over its duration is calculated. Specifically, for each time sampling point within the entire duration, the sum of rainfall intensities of all spatial grids within the target area at that time sampling point is calculated to obtain the total instantaneous rainfall intensity of the region at that time sampling point. The total instantaneous rainfall intensity of the region at all time sampling points throughout the entire duration of the heavy rainfall event is summarized to construct a time series of total instantaneous rainfall intensity of the region. The time series of total instantaneous rainfall intensity of the region is then integrated over the entire duration, and the scalar result obtained from the integration is the cumulative integral value. Information entropy, direction offset vector, and integral cumulative value are combined to construct a multidimensional feature vector. Specifically, the information entropy, two component values ​​of the direction offset vector, and integral cumulative value calculated for the same heavy rainfall event, a total of four numerical features, are arranged and combined into a one-dimensional numerical array in a preset fixed order to form a multidimensional feature vector.

[0010] Furthermore, the feedforward neural network model performs numerical processing based on its internal parameterized decay function and anisotropic diffusion convolution kernel. Specifically, the feedforward neural network model designs two computational components in its network structure, including a parameterized decay function component and an anisotropic diffusion convolution kernel component. The parameterized decay function component is parameterized based on the integral accumulation value in the multidimensional feature vector. Its function is to define a decay function that follows an exponential law to simulate the characteristic that the intensity of the virtual pressure field induced by rainfall infiltration decays over time. The anisotropic diffusion convolutional kernel component is a two-dimensional weight matrix whose internal parameters are constrained during the training of the feedforward neural network model. This constraint enables the convolutional kernel to possess anisotropic spatial diffusion coefficients. Specifically, it has a dominant diffusion coefficient in the direction parallel to the direction offset vector in the multidimensional feature vector, and a different diffusion coefficient in another dominant direction perpendicular to this direction. This simulates the non-uniform propagation phenomenon of the virtual pressure field in space along the direction offset vector. The multidimensional feature vector is input into the pre-trained feedforward neural network model. Through the collaborative computation of the two components, the feedforward neural network model simulates the dynamic process of the virtual pressure field propagating from the end of the current monitoring period to the future, and finally outputs the pressure field state tensor arranged in the order of future predicted time points.

[0011] Furthermore, trend decomposition is performed on historical groundwater level data to obtain a background water level digital field. Specifically, this involves: extracting long-term time-series water level data corresponding to each spatial grid within the target area from historical gridded groundwater level data; applying a time series decomposition algorithm to the long-term time-series water level data of each spatial grid to decompose it into a long-term trend term and a seasonal periodic term; superimposing the two to obtain a background water level change benchmark sequence for each spatial grid; extrapolating the background water level change benchmark sequence of each spatial grid to these future prediction time points to generate the background water level value of each spatial grid at each future prediction time point; finally, reorganizing the background water level values ​​of all spatial grids at the same future prediction time point according to the spatial grid structure to obtain a background water level digital field corresponding to the pressure field state tensor in the time dimension. The water level response prediction field is obtained by performing tensor addition operations on the pressure field state tensor and the background water level digital field. Specifically, the pressure-water level conversion coefficient is determined based on the physical characteristics of the aquifer. For each future prediction time point, the pressure value of each spatial grid of the pressure field state tensor corresponding to that future prediction time point is multiplied by the pressure-water level conversion coefficient to obtain the water level change field. The water level change field is then added grid by grid to the same spatial grid value of the background water level digital field corresponding to the same future prediction time point to obtain the water level response prediction field for that future prediction time point. This operation is repeated for all future prediction time points to obtain a sequence of water level response prediction fields arranged in chronological order.

[0012] Furthermore, the water level change prediction curves of preset prediction points are extracted from the water level response prediction field, and the groundwater level rise rate is output after differential operation. Specifically, from the water level response prediction field sequence arranged in the order of future prediction time points, the water level prediction values ​​of the pre-selected groundwater prediction points at each future prediction time point are extracted in time order. These water level prediction values ​​are connected in time order to form a prediction curve of water level changing with time. The average slope of the prediction curve is calculated using the numerical differential method. The average slope value is the groundwater level rise rate and is output as the final result.

[0013] Compared with the prior art, the beneficial effects of the present invention are: This invention, by creating a spatiotemporal matrix containing both rainfall and water level channels, achieves for the first time deep fusion and collaborative expression of two types of heterogeneous monitoring data within a unified spatial grid and time series framework, overcoming the information fragmentation problem caused by the isolated processing of the two data in existing technologies. Furthermore, by performing spatiotemporal scanning and feature extraction on heavy rainfall events, a multidimensional feature vector is constructed that integrates the spatial heterogeneity of the event (information entropy), the directional relationship of the rainfall subject relative to the prediction point (direction offset vector), and the total energy input of the event (integral cumulative value). This digitizes and vectorizes the key physical attributes of a complex rainfall event, providing subsequent models with driving signals that are both rich in information and conform to physical intuition. The core innovation of the pre-trained feedforward neural network model used in this invention lies in the integration of a parameterized decay function and an anisotropic diffusion convolution kernel. These two components are not ordinary mathematical operation units, but rather embed the physical laws of the decay and diffusion of groundwater pressure with time and space into the black box calculation process of the neural network in the form of learnable parameters. This makes the information propagation and transformation process inside the model spontaneously subject to the soft constraints of physical laws when the model makes predictions based on data features, thereby simulating the complex process of non-uniform diffusion of pressure in underground media. This invention employs time-series decomposition technology to pre-extract the long-term background trend and periodic fluctuations in historical water levels, forming a background water level digital field. This field is then superimposed and synthesized with the pressure field state tensor, representing the instantaneous response to rainfall. This decomposition-superposition prediction framework clearly distinguishes and quantifies the water level benchmark determined by long-term natural factors and the water level change increment directly caused by this specific heavy rainfall event in the final prediction result. This not only significantly improves the accuracy and relevance of the prediction of the water level rise component, making the prediction results more instructive for disaster early warning, but also overcomes the shortcomings of existing methods where trend noise and event signals are mixed in the prediction results, making it difficult to directly use them to assess the immediate risk of the event. Finally, by extracting the prediction curve and performing differential operations, the key decision-making indicator of the water level rise rate is directly output, enabling the technical effects of this invention to effectively serve disaster prevention and mitigation practices and provide high-value prediction information. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a curve showing the fitting of the integral cumulative value and decay rate in this invention. Figure 3 This is a bar graph showing the convolution kernel parameter 1, convolution kernel parameter 2, and decay rate of the present invention. Detailed Implementation

[0015] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0016] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0017] Example: Please see Figures 1-3 The present invention provides a technical solution: A method for predicting groundwater rise under heavy rainfall conditions, comprising the following steps: Step 1: Obtain rainfall intensity data and groundwater level data for the target area, and apply a spatial grid interpolation algorithm to convert the groundwater level data into regular grid data with the same spatial resolution as the rainfall intensity data, generating a spatiotemporal matrix containing rainfall intensity channels and groundwater level channels; The groundwater level data is converted into regular grid data with the same spatial resolution as the rainfall intensity data. Specifically, the sampling interval between adjacent time sampling points is preset. Multiple groundwater level measurement points with different geographical coordinates are selected within the target area to obtain groundwater level data at each time sampling point during the current monitoring period. Simultaneously, rainfall intensity data covering the same target area within the same monitoring period is acquired. This rainfall intensity data has a unified regular geographical grid and a fixed spatial resolution, and each time sampling point has corresponding gridded rainfall intensity data. Based on the regular geographical grid used for the rainfall intensity data, for each time sampling point, the Kriging spatial interpolation algorithm is used to calculate and generate the corresponding groundwater level estimate for each spatial grid in the regular geographical grid, using the groundwater level data and its geographical coordinates at that time sampling point as input. Through the above operations, the originally discrete groundwater level data is converted into gridded groundwater level data that is aligned with the rainfall intensity data in the spatial grid framework and synchronized with it in time. In this invention, the current monitoring period refers to the core time window used to construct the spatiotemporal matrix, identify heavy rainfall events, and extract event features. This period is located after the historical period and before the prediction period, and is a short, continuous period immediately adjacent to the prediction start time. The determination of the current monitoring period follows the following principles: taking the peak time of the heavy rainfall event to be predicted as the center, extending forward 72 hours as the starting point to ensure that the background state of groundwater level before the event occurs can be fully captured; extending backward to the end time of the event as the ending point to fully cover the entire process of the event from occurrence to dissipation. The total length of the entire current monitoring period depends on the duration of the event, generally several days to a week, and only includes the current heavy rainfall event and its direct impact period, excluding the delayed response period after the event ends. It is important to emphasize that the current monitoring period and the historical period used for trend decomposition are two independent periods that do not overlap in time. The historical period is located before the current monitoring period and is a longer historical time period used to extract long-term patterns of background water levels (such as interannual variations and seasonal fluctuations). The sampling interval between adjacent time sampling points is preset to 1 hour. This setting is determined by the data source: rainfall intensity data provided by operational meteorological radar or automatic rain gauge networks, as well as the standard reporting interval of most automated groundwater monitoring systems, is 1 hour. Choosing a 1-hour interval can effectively capture the evolution of heavy rainfall processes (such as peak value and duration) on a time scale, while avoiding the huge storage and computational burden caused by minute-level data. It is a common standard for balancing the accuracy of hydrological process characterization and data processing efficiency. In the spatial resolution of the regular geographic grid, the grid size directly adopts the original resolution of the input gridded rainfall intensity data. For example, if the spatial resolution of the radar quantitative precipitation estimation product used is 1 km × 1 km, then the spatial grid size in this invention is 1 square kilometer. In a specific embodiment of the present invention, considering the fundamental differences in the spatial representation of rainfall intensity data and groundwater level data, direct correlation analysis would lead to spatial information mismatch, thereby affecting the accuracy of the prediction model. Rainfall intensity data generally comes from regional observation methods such as weather radar, and is naturally presented in the form of a regular grid, with each spatial grid representing the rainfall intensity of a fixed geographical area at a certain moment. Groundwater level data, on the other hand, comes from discrete monitoring wells deployed at specific locations. Each well can only provide information on water level changes at its location and cannot directly reflect the water level field of the entire area. Therefore, converting discrete groundwater level monitoring data into regular grid data with the same spatial resolution as rainfall intensity data through spatial interpolation algorithms is the primary prerequisite for conducting collaborative analysis of the two within a unified spatial framework. This step ensures that all subsequent calculations are based on the same geographic spatial grid, enabling accurate correlation and comparison between the spatial distribution characteristics of rainfall and the spatial response of groundwater level. It should be clarified that the regular geographic grid framework mentioned in this invention refers to a predefined digital spatial division system that covers the entire target area. The framework consists of a fixed number of rows and columns, and each spatial grid has a uniform size and a clear geographic coordinate range, thereby discretizing the continuous geographic space into a series of regular spatial grids. This framework directly adopts the coordinate system and resolution used in the self-acquired gridded rainfall intensity data to ensure the consistency of the spatial reference. The process of converting discrete groundwater level data into regular grid data is specifically implemented through a spatial gridding interpolation algorithm. This invention preferably uses the Kriging spatial interpolation algorithm to complete this task. The reason for choosing the Kriging algorithm is that it not only considers the distance relationship between the monitoring point and the grid point to be estimated, but also incorporates the structure and correlation of the spatial distribution of groundwater level through the variogram model, which can provide the statistically optimal and unbiased estimate, thereby more realistically reflecting the spatial continuous change trend of the water level field. In practice, the first step is to unify the time reference for rainfall intensity data and groundwater level data. This invention sets a unified time sampling interval to ensure that each identical time sampling point is processed. For a given time sampling point, the algorithm takes the measured water level values ​​and their precise geographic coordinates of all groundwater level monitoring points at that moment as input, and the center coordinates of each spatial grid to be calculated in the regular geographic grid used for rainfall intensity data as the estimated position. Based on the optimal linear unbiased estimation principle of the Kriging algorithm, the algorithm calculates the estimated groundwater level of that spatial grid at that moment. This process is repeated for all time sampling points, and finally a gridded groundwater level dataset that is completely synchronized with the rainfall intensity data in terms of time series and completely aligned with the rainfall intensity data in terms of spatial grid is generated. A spatiotemporal matrix containing rainfall intensity channels and groundwater level channels is generated. Specifically, for the gridded rainfall intensity data and gridded groundwater level data corresponding to each time sampling point within the monitoring period, they are registered according to the same spatial grid coordinates. The two are superimposed in the data channel dimension to form a composite data layer with dual channels. The composite data layers corresponding to all time sampling points are arranged sequentially in chronological order to construct a three-dimensional spatiotemporal matrix. The first and second dimensions of the spatiotemporal matrix correspond to the row index and column index of the spatial grid, respectively, and the third dimension corresponds to the time sampling point index. The rainfall intensity data and the groundwater level data are stored as independent channel data in the spatiotemporal matrix. After completing the spatial gridding alignment of rainfall intensity data and groundwater level data, a unified data structure needs to be constructed to integrate the two and perform spatiotemporal analysis. This invention defines and generates a three-dimensional data container called a spatiotemporal matrix. The specific operation for generating this matrix is ​​as follows: for each time sampling point, the gridded rainfall intensity data and gridded groundwater level data that have been aligned at this moment are paired according to a completely consistent spatial grid row and column order, and in terms of data organization, they are regarded as two independent data layers, namely channels. The rainfall intensity value constitutes the rainfall intensity channel, and the groundwater level estimate constitutes the groundwater level channel. Then, all these dual-channel data layers corresponding to all time sampling points are stacked in strict chronological order on the third dimension. The resulting spatiotemporal matrix has its first and second dimensions defining the spatial location, and the third dimension defining the temporal evolution. Through this structure, the originally independent multi-source spatiotemporal data are integrated into a unified computing framework, enabling subsequent steps to efficiently perform continuous spatiotemporal scanning and feature extraction of the rainfall process and water level status of the entire region.

[0018] Step 2: Continuously scan the rainfall intensity data in the spatiotemporal matrix, use the threshold detection method to identify continuous periods exceeding the set heavy rainfall threshold and record them as heavy rainfall events. For each heavy rainfall event, calculate the information entropy of its rainfall intensity distribution, determine the directional offset vector of its rainfall intensity weighted centroid relative to the preset prediction point, and calculate the integral cumulative value of rainfall intensity over time during the duration of the event. Combine the information entropy, directional offset vector and integral cumulative value to construct a multidimensional feature vector. In a specific embodiment of the present invention, the core objective of conducting in-depth analysis of the spatiotemporal matrix generated in step one is to accurately extract those heavy rainfall segments that have a significant impact on groundwater dynamics from the continuous rainfall process, and to construct a digital feature expression that can comprehensively characterize its spatiotemporal structure. Directly using the original, continuous grid rainfall data for prediction will introduce a large amount of noise information that is unrelated to the rapid response of groundwater and lead to a dramatic increase in computational complexity. Therefore, identifying and extracting independent heavy rainfall events and calculating a set of feature indicators with clear physical or statistical significance for them is a key step in transforming complex natural environmental processes into standard inputs that can be efficiently processed by machine learning models. A threshold detection method is used to identify continuous periods exceeding a set heavy rainfall threshold and record them as heavy rainfall events. Specifically, the rainfall intensity data in the spatiotemporal matrix is ​​scanned sequentially along the time dimension. For each time sampling point, the average rainfall intensity value of all spatial grids in the target area corresponding to that time sampling point is calculated. The average rainfall intensity value of each time sampling point is continuously compared with a preset heavy rainfall threshold to identify all time sampling points whose average rainfall intensity value exceeds the heavy rainfall threshold. Among the identified time sampling points, a series of time sampling points that are temporally continuous and whose time interval between two adjacent time sampling points does not exceed the preset maximum allowable interval are merged into an independent heavy rainfall event. At the same time, the start time and end time of each heavy rainfall event are recorded. The first step is the automatic identification of heavy rainfall events. This invention uses a threshold-based detection method to achieve this goal. The preset heavy rainfall threshold is a critical value determined in advance based on the statistical characteristics of the historical hourly rainfall intensity sequence of the target area. Its physical meaning is to distinguish between ordinary rainfall and heavy rainfall that causes significant groundwater level fluctuations. An objective method for determining this threshold is to statistically analyze the long-term historical hourly rainfall intensity sequence of the target area and take its 95th percentile as the threshold. This value means that only 5% of the hourly rainfall events in history have exceeded this value, which can effectively screen out extreme rainfall events. The purpose of setting this threshold is to filter out ordinary rainfall that has a weak disturbance to the groundwater level from continuous rainfall and focus on heavy rainfall processes that trigger significant responses. In practice, the rainfall intensity channel data in the spatiotemporal matrix is ​​read sequentially along the time dimension. For each time sampling point, the average rainfall intensity value of all spatial grids in the target area under that time sampling point is calculated. This value represents the overall intensity level of rainfall at that time sampling point. Then, the average rainfall intensity value of each time sampling point is compared with the preset heavy rainfall threshold one by one. All time sampling points with average rainfall intensity values ​​exceeding the threshold are marked as candidate points. Then, time continuity analysis is performed on these candidate points. This invention sets a maximum allowable time interval to tolerate short-term fluctuations or minor discontinuities in the data. A series of candidate points that are continuous in time and whose actual time interval between any two adjacent candidate points does not exceed this maximum allowable interval are grouped into an independent heavy rainfall event. At the same time, the start time and end time of the heavy rainfall event are recorded, thereby clearly defining its range on the time axis. The maximum allowable interval is set to 2 to 6 hours. For example, setting it to 3 hours is acceptable. In meteorology, a continuous rainfall process may be interrupted below the threshold due to short-term fluctuations in intensity. Setting a reasonable maximum allowable interval allows rainfall processes caused by the same weather system with short-term weakening to be grouped into the same event, thereby more completely depicting the total impact of a weather process. If the interval is too long, different events will be incorrectly merged. For any heavy rainfall event, the information entropy of its rainfall intensity distribution is calculated as follows: the time sampling point when the rainfall intensity reaches the global maximum value during the duration of the heavy rainfall event is taken as the peak time; the rainfall intensity values ​​of all spatial grids in the target area at the peak time are extracted, and these rainfall intensity values ​​are normalized so that their sum is 1, thus forming a spatial probability distribution characterizing the spatial distribution of rainfall intensity, where the value of each spatial grid represents the relative probability of the rainfall intensity at that location relative to the total intensity of the area; based on the Shannon entropy calculation formula in information theory, the information entropy of this spatial probability distribution is calculated, and the resulting entropy value is the information entropy index characterizing the heavy rainfall event. This information entropy index can sensitively reflect the degree of concentration or dispersion of rainfall in space: the lower the entropy value, the more concentrated the rainfall is in a few spatial grids; the higher the entropy value, the more uniformly dispersed the rainfall in space, providing a quantitative indicator of the spatial structure of rainfall for the model; For any heavy rainfall event, the directional offset vector of its rainfall intensity weighted centroid relative to the preset prediction point is determined. Specifically, the peak time of the heavy rainfall event is extracted, and the average rainfall intensity value of each spatial grid at the peak time is used as the weight. The weighted average values ​​of all spatial grids in the east-west and north-south directions are calculated respectively. The coordinate point formed by these two weighted average values ​​is the coordinate of the rainfall intensity weighted centroid at the peak time of the heavy rainfall event. The plane vector pointing from the coordinate of the rainfall intensity weighted centroid to the coordinate of the preset prediction point is calculated. The plane vector is normalized to make its magnitude 1, resulting in a unit vector that only represents directional information. The resulting unit vector is the directional offset vector. The weighted centroid coordinates of rainfall intensity are the center position coordinates obtained by spatial weighting calculation based on rainfall intensity. They represent the location of the mass center or main energy accumulation area of ​​the heavy rainfall event in space. The calculation method is to use the rainfall intensity values ​​of each grid at the aforementioned peak time as weights, and calculate the weighted average of the east-west and north-south coordinates of the grid to obtain the latitude and longitude coordinates of the weighted centroid. The preset prediction point is the specific geographical coordinate of the designated location that needs to be focused on and used for groundwater level prediction. The direction offset vector is defined as the directed line segment from the weighted centroid of this event to the preset prediction point. In order to eliminate the influence of distance and purify the direction information, this invention normalizes the plane vector by scaling its length to 1, resulting in a unit vector containing only direction information. This direction offset vector reveals the orientation relationship of the main area of ​​heavy rainfall relative to the point to be predicted, which is crucial for determining the main direction of pressure propagation and water level response. In this invention, the preset prediction point refers to a specific geographical location within the target area that is pre-selected and requires key monitoring and prediction of groundwater level changes. This point is essentially a spatial anchor point with precise geographical coordinates, corresponding to a physical observation well or key engineering point with long-term continuous monitoring data. Its selection is not arbitrary and must simultaneously meet three core conditions: located in the same hydrogeological unit as the rainfall area to ensure hydraulic connection; possessing sufficient historical water level data to support model training; and located in a key area with practical disaster prevention or management significance (such as an urban flooding risk area or an engineering sensitive point). The coordinates of this preset prediction point are deeply integrated into the entire prediction process, playing a pivotal role. In the feature extraction stage, it is used together with the centroid coordinates of the rainfall event to calculate the direction offset vector, quantifying the relative orientation of the rainfall subject and the prediction target. In the neural network prediction stage, this vector guides the anisotropic diffusion convolution kernel inside the model to adjust its principal axis direction, thereby simulating the transmission of pressure along this advantageous path. Finally, in the result output stage, the grid value corresponding to the coordinates of this point is directly extracted to generate a unique water level change curve and rise rate prediction value. The directional offset vector is a unit vector that only indicates direction. Its magnitude is always 1. The calculation steps are to first find the original vector (with the centroid pointing to the prediction point), and then divide it by its own magnitude. The physical meaning of the weighted centroid is the balance point of rainfall energy on the spatial grid. Normalization is to purify the directional features and eliminate the influence of distance factors in the subsequent neural network model. Only the relationship between direction and pressure diffusion anisotropy needs to be learned, while distance information is implicit or expressed by other features, such as centroid position and decay function. For any heavy rainfall event, the cumulative integral of rainfall intensity over time is calculated over its duration. Specifically, for each time sampling point within the entire duration, the sum of rainfall intensities across all spatial grids within the target area at that time sampling point is calculated to obtain the total instantaneous rainfall intensity of the region at that time sampling point. The total instantaneous rainfall intensity of the region at all time sampling points throughout the entire duration of the heavy rainfall event is summarized to construct a time series of the total instantaneous rainfall intensity of the region. Numerical integration is performed on the time series of the total instantaneous rainfall intensity of the region over the entire duration, for example, using the trapezoidal rule. The scalar result obtained from the integration is the cumulative integral value. This cumulative integral value represents the total energy or total water volume input into the groundwater of the target area during this heavy rainfall event (ignoring evaporation, runoff, etc.). It comprehensively reflects the intensity and duration of rainfall and is a key energy indicator for measuring the total rainfall of the event. It is also the core input for determining the initial pressure intensity in the subsequent decay function. Information entropy, direction offset vector and integral cumulative value are combined to construct a multidimensional feature vector. Specifically, the information entropy, two component values ​​of the direction offset vector and integral cumulative value calculated for the same heavy rainfall event, a total of four numerical features, are arranged and combined into a one-dimensional numerical array in a preset fixed order to form a multidimensional feature vector. Finally, the three features mentioned above are standardized and combined to construct the multidimensional feature vector defined in this invention. Specifically, the information entropy scalar value, the two component values ​​of the direction offset vector, and the integral accumulation scalar value, all calculated for the same heavy rainfall event, are arranged and concatenated in a pre-set, fixed order to form a one-dimensional numerical array. This multidimensional feature vector serves as the unique digital signature of the heavy rainfall event and will be used as the unified input for subsequent machine learning models. It concisely expresses the spatiotemporal structure of the event and is a bridge connecting physical observation and prediction.

[0019] Step 3: Input the multidimensional feature vector into the pre-trained feedforward neural network model. The model performs numerical processing based on the internal parameterized decay function and anisotropic diffusion convolution kernel to calculate the pressure field state tensor of the target region. In a specific embodiment of the present invention, the core task of step 3 is to map the multidimensional feature vector constructed in step 2, which characterizes the spatiotemporal structural features of a single heavy rainfall event, into a pressure field state tensor that can describe the initial spatial distribution of the virtual pressure disturbance induced by the event on the aquifer. The key to this step is to design a computational model that can combine event features with hydrogeological and physical mechanisms, rather than establishing a black-box end-to-end mapping. To this end, the present invention constructs and adopts a pre-trained feedforward neural network model. This feedforward neural network model innovatively integrates two learning computational components with clear physical meaning in the standard spatial network structure: a parameterized decay function component and an anisotropic diffusion convolution kernel component. This design enables the model not only to learn complex patterns in the data, but also to comply with the basic physical laws of groundwater pressure propagation to a certain extent, thereby improving the rationality of the prediction results. The feedforward neural network model performs numerical processing based on an internal parameterized decay function and anisotropic diffusion convolution kernel. Specifically, the feedforward neural network model designs two computational components in its network structure, including a parameterized decay function component and an anisotropic diffusion convolution kernel component. The parameterized decay function component is parameterized based on the integral accumulation value in the multidimensional feature vector. Its function is to define a decay function that follows an exponential law to simulate the characteristic that the intensity of the virtual pressure field induced by rainfall infiltration decays over time. First, it is necessary to clarify the meaning of the pressure field state tensor. In this invention, the pressure field state tensor is a three-dimensional data structure. Its first two dimensions correspond to the rows and columns of the regular geographic grid of the target area. In some implementations, the third dimension represents different time steps or different pressure components, but here the main output is the pressure spatial distribution field, which is regarded as a tensor with the third dimension of 1. It characterizes the spatial distribution of a virtual pressure field generated by a heavy rainfall event through infiltration at the top of a hypothetical homogeneous unconfined aquifer at that moment. This pressure is not directly measured, but is an intermediate physical quantity used to connect rainfall input and water level response. Its unit is related to the water head height. The virtual pressure field is a distribution of pressure increments caused by the additional load exerted by rainwater infiltration on underground aquifers, derived from the model. It can be imagined as the pressure changes at different locations below the water surface after a stone (heavy rainfall event) is thrown into a pond (aquifer), except that this occurs in a porous medium. This pressure itself is a model variable, its core purpose being as an intermediate physical quantity. It can be effectively learned and inferred from rainfall characteristics through a neural network, and can also be converted into the actual observed groundwater level changes through a clear physical formula (i.e., multiplied by a conversion coefficient). Therefore, the virtual pressure field is the bridge between data-driven and physical mechanisms in this invention. To calculate this virtual pressure field from rainfall characteristics, this invention designs a pre-trained feedforward neural network model. The innovation of this model lies not in its depth, but in the key physical mechanisms for simulating groundwater dynamics embedded within it, specifically implemented through two computational components: a parameterized decay function component and an anisotropic diffusion convolution kernel component. The parameterized decay function component in the feedforward neural network model simulates the physical process of the intensity of the aforementioned virtual pressure pulse decaying over time. Decay is an inherent characteristic of groundwater pressure disturbances propagating in the medium, primarily caused by the water storage and release effects of the aquifer and hydraulic diffusion. This component does not use a fixed decay formula but defines an exponential decay function framework controlled by the model parameters, in the form of: The decay relationship, in which the key parameter is the decay rate. Instead of being a preset constant, it is dynamically generated through an auxiliary sub-neural network. This sub-network is primarily based on the cumulative integral value in the input feature vector, because the cumulative integral value characterizes the total rainfall energy of the event; the greater the energy, the stronger the initial pressure. The larger the value, the different the attenuation process will be due to the nonlinear response of the medium. It is the initial pressure intensity to be determined. It is a time variable, and this relationship will be incorporated into the forward propagation calculation of the network, affecting the change of the pressure field with the simulation time step; The anisotropic diffusion convolutional kernel component is a two-dimensional weight matrix whose internal parameters are constrained during the training of the feedforward neural network model. This constraint enables the convolutional kernel to possess anisotropic spatial diffusion coefficients. Specifically, it has a dominant diffusion coefficient in the direction parallel to the direction offset vector in the multidimensional feature vector, and a different diffusion coefficient in another dominant direction perpendicular to this direction. This simulates the non-uniform propagation phenomenon of the virtual pressure field in space along the direction offset vector. The multidimensional feature vector is input into the pre-trained feedforward neural network model. Through the collaborative computation of the two components, the feedforward neural network model simulates the dynamic process of the virtual pressure field propagating from the end of the current monitoring period to the future. Finally, it outputs pressure field state tensors arranged in the order of future prediction time points. Each pressure field state tensor corresponds to a future prediction time point after the end of the current monitoring period. These tensors together cover the entire future prediction period, providing a time-matching basis for the superposition operation with the background water level digital field in subsequent steps. The anisotropic diffusion convolutional kernel component in the feedforward neural network model functions to simulate the physical process of non-uniform diffusion of virtual pressure fields in space. The permeability of underground media often has directionality, such as along the strike of strata, the direction of fracture development, or the direction of ancient river channel extension, which has stronger water conductivity and leads to faster pressure propagation. This invention utilizes the spatial structural characteristics of rainfall events themselves to indicate this dominant direction. This component is essentially a two-dimensional convolutional kernel, i.e., a weight matrix, but its weight distribution is constrained to have anisotropic characteristics. Anisotropy here specifically refers to the different smoothing (or diffusion) effect intensities of the convolutional kernel in different spatial directions. This constraint is achieved by associating the internal weight matrix of the convolutional kernel with the direction offset vector in the input multidimensional feature vector during the model training phase. The direction offset vector indicates the direction from the main rainfall area to the prediction point. During training, by designing a specific loss function or network structure, the weight distribution of the convolutional kernel can generate stronger spatial correlation transmission, i.e., a higher equivalent diffusion coefficient, in the direction parallel to the offset vector, while its transmission effect is weaker in the direction perpendicular to the vector. In this way, when the convolutional kernel acts on an initial pressure distribution in the neural network, it will cause the pressure value to be smoothed or diffused over a larger range in the dominant direction, thus mathematically simulating the physical phenomenon that pressure is transmitted faster along a specific direction. For the Secondary diffusion-decay iteration ( ): in, Let be the pressure field state tensor of the k-th iteration, with its initial value... Based on the multidimensional feature vectors of the base fully connected layer Generation, final pressure field , This is a two-dimensional convolution operation, typically with zero padding to maintain spatial dimensions, used to simulate the spatial diffusion of pressure within an aquifer. For anisotropic diffusion convolution kernels, For convolution kernel parameters, It is determined by parameters A parameterized neural network mapping function whose input is a direction offset vector. (Unit vector), its magnitude is always equal to 1, and The X component of the direction offset vector. The Y component of the direction offset vector, and the range of values ​​for both are... The specific values ​​are determined based on the relative orientation of the rainfall weighted centroid and the preset prediction point; The constraint is that it is parallel to The direction in which the diffusion coefficient is dominant, i.e. This is used to simulate the faster transmission of groundwater pressure along the dominant direction (from the main rainfall event to the preset prediction point). Parallel to the direction offset vector Equivalent diffusion coefficient in the direction, The offset vector perpendicular to the direction Equivalent diffusion coefficient in the direction; This is element-wise multiplication (Hadamard product). The decay rate parameter controls how quickly pressure disturbances decay over time, reflecting the aquifer's water storage and release capacity. It is determined by parameters Parameterized sub-neural networks This is the accumulated points. To simulate the time step, it is always consistent with the preset time interval. The total number of iterations is determined by the total simulation duration. and step length Decide, This is used to simulate the time evolution from the start of rainfall to the predicted time. Sub-neural network The input is the direction offset vector, and the output is a set of feature mapping parameters. ,in Defined as the principal axis diffusion gain coefficient, named convolution kernel parameter 1, it is used to adjust the diffusion intensity along the direction of the offset vector. Defined as the orthogonal axis scaling factor, named kernel parameter 2, it is used to adjust the diffusion distribution range perpendicular to this direction. These two parameters are then mapped to a two-dimensional convolution kernel. The coefficients of the weight distribution matrix; Values ​​ranging from 0.8 to 2.0 are significantly greater than The value ranges from 0.2 to 1.1, and varies with the components of the direction offset vector. The value changes, The dynamic fluctuations demonstrate that the convolution kernel is dynamically generated based on the spatial structure. hypothesis function It is a simple nonlinear function (simulating the nonlinear mapping of a neural network): ,in , , For constant parameters, along with (Rainfall energy) increases, This also increases accordingly, reflecting the characteristic that the pressure disturbance brought about by heavy rainfall disappears faster or the nonlinear response is more intense; Table 1 shows the specific data for some heavy rainfall events and their parameters.

[0020] Table 1 Data Statistics Table Analysis of the data revealed a clear correlation between different characteristic parameters. The data showed a significant positive correlation between the cumulative integral value and the decay rate: as the cumulative integral value increased, the decay rate also increased. For example, the cumulative integral value of heavy rainfall event 14 was 3.210, corresponding to a decay rate of 0.018; while for heavy rainfall event 2, the cumulative integral value was 35.678, and the decay rate increased to 0.051. This indicates that heavy rainfall events with greater total rainfall energy exhibit a faster dissipation process of the pressure disturbance induced in the aquifer, consistent with the physical intuition that high energy input leads to a stronger but also shorter hydrological response. When analyzing the relationship between directional offset and convolution kernel parameters, it can be found that the magnitude of the convolution kernel parameters is closely related to the specific direction of the directional offset vector. For example, the directional offset X component of heavy rainfall event 3 is 0.985, and the corresponding convolution kernel parameter 1 is 1.543; while the directional offset X component of heavy rainfall event 15 is -0.985, and the corresponding convolution kernel parameter 1 is 1.398. The convolution kernel parameter 2 also shows a similar change pattern. This correlation indicates that the neural network can adaptively adjust the shape of the pressure diffusion operator according to the specific orientation of the rainfall spatial distribution relative to the prediction point, thereby simulating the physical mechanism of the anisotropic transmission process of groundwater pressure within the model. Further observation of the data reveals that even events with similar cumulative integral values ​​exhibit subtle differences in their decay rates. For example, the cumulative integral values ​​of heavy rainfall event 1 and heavy rainfall event 12 are 8.214 and 10.123, respectively, with corresponding decay rates of 0.027 and 0.030. Simultaneously, samples with different directional offset vectors, such as heavy rainfall event 3 and heavy rainfall event 11, also show different numerical characteristics in their convolution kernel parameter combinations. These complex correspondences collectively demonstrate that the feedforward neural network model does not perform a simple linear mapping, but rather comprehensively analyzes and represents the multidimensional characteristics of rainfall events through its internal network parameters, ultimately forming a digital twin capable of reflecting complex hydrophysical processes.

[0021] The entire pre-trained feedforward neural network model is built around the two aforementioned physical guidance components. The model receives a multi-dimensional feature vector as input, consisting of information entropy, orientation offset vector, and integral accumulation. Internally, the input features are first encoded and transformed through a series of fully connected layers. Then, the parameterized decay function component determines the decay behavior based on the integral accumulation branch. Simultaneously, the anisotropic diffusing convolution kernel component determines the anisotropic principal axis direction of its weight matrix based on the orientation offset vector branch. The core computational process of the feedforward neural network model includes the following steps: an initial, coarse-grained feature vector is generated by a basic fully connected layer based on the feature vector. The initial pressure distribution is roughened; then, anisotropic diffusion convolution kernels are used to perform multiple convolution operations on this initial distribution to simulate the spatial diffusion process of pressure. During this process, a parameterized decay function is applied to the pressure field after each diffusion to simulate its decay over time. Finally, after a series of such diffusion-decay cycle calculations, the network outputs a stable and refined two-dimensional grid data, namely the pressure field state tensor covering the entire target area. This tensor is the virtual groundwater pressure disturbance spatial pattern predicted by the model and uniquely determined by the characteristics of this heavy rainfall event. It will serve as the basis for the next step of calculating the actual changes in water level.

[0022] Step 4: Perform trend decomposition on historical groundwater level data to obtain background water level digital field. Perform tensor addition operation on pressure field state tensor and background water level digital field to obtain water level response prediction field. Extract water level change prediction curves of preset prediction points from water level response prediction field, and output the groundwater level rise rate after differential operation. In a specific embodiment of the present invention, the core task of step four is to fuse the pressure field state tensor generated in step three, which characterizes the virtual pressure disturbance induced by heavy rainfall, with the background water level digital field that reflects the natural fluctuation law of groundwater, and finally calculate the predicted value of the groundwater level rise rate at a preset specific location. The technical necessity lies in the fact that the groundwater level at any given time is the result of the combined effects of multiple factors, including long-term regional replenishment and discharge trends, periodic seasonal changes, and the instantaneous response triggered by this heavy rainfall event. If these components are not separated and synthesized, the prediction results will not be able to accurately depict the impact of this event alone. First, it's important to clarify that the pressure field state tensor output by the feedforward neural network model in step 3 implicitly defines a sequence of future predicted time points. Specifically, based on the characteristics of heavy rainfall events identified during the current monitoring period, the model outputs a series of pressure field distributions uniformly arranged at preset time steps through iterative calculations using a built-in parameterized decay function and anisotropic diffusion convolution kernels. The times corresponding to these pressure field distributions constitute the sequence of future predicted time points. The starting point of this sequence is after the end of the current monitoring period, and its length and time step are determined by the model's design parameters. Therefore, the pressure field state tensor naturally defines a set of future predictions. Each pressure field state tensor corresponds to a unique future prediction time point. Trend decomposition is performed on historical groundwater level data to obtain the background water level digital field. Specifically, long-term time series water level data corresponding to each spatial grid in the target area is extracted from the gridded groundwater level data of historical periods before the current monitoring period. The historical period refers to the independent time interval before the current monitoring period used to extract the long-term evolution law of the groundwater system. It is strictly distinguished from the current monitoring period on the time axis and does not overlap with it. Generally, complete hydrological year data of three to five consecutive years are selected to ensure that the long-term trend and seasonal cycle can be accurately separated. A time series decomposition algorithm is applied to the long-term water level data of each spatial grid to decompose it into a long-term trend term and a seasonal periodic term. The two are superimposed to obtain the background water level change benchmark sequence for each spatial grid. This background water level change benchmark sequence describes the natural law of water level change over time at that point without the interference of the current heavy rainfall event. The pressure field state tensor output by the feedforward neural network model is a three-dimensional tensor sequence arranged in the order of future prediction time points. Each pressure field state tensor in this sequence corresponds to a future prediction time point after the end of the current monitoring period. Using the future prediction time point sequence defined by the pressure field state tensor as a benchmark, the background water level change benchmark sequence of each spatial grid is extrapolated to these future prediction time points, thereby generating the background water level values ​​of each spatial grid at each future prediction time point. This extrapolation operation is the key step in establishing a correspondence between the background water level digital field and the pressure field state tensor. It is precisely because the two are generated based on the exact same future prediction time point sequence that the subsequent superposition operation has temporal consistency. Finally, the background water level values ​​of all spatial grids at the same future prediction time point are reorganized according to the original spatial grid structure to obtain the two-dimensional background water level digital field corresponding to that future prediction time point. For each future prediction time point, a corresponding background water level digital field can be generated. All these background water level digital fields are arranged in chronological order to form a background water level digital field sequence. Each background water level digital field in this sequence has a one-to-one correspondence with the pressure field state tensor of the pressure field state tensor sequence corresponding to the same future prediction time point. The background water level digital field is defined as: a water level field that excludes short-term random fluctuations and interference from specific heavy rainfall events, and only contains the long-term evolution law and periodic changes of groundwater itself. The specific operation to obtain this field is as follows: for each spatial grid in the regular geographic grid in the target area, extract its corresponding long-term water level data; apply a time series decomposition algorithm (such as seasonal decomposition method) to this time series to decompose it into a long-term trend term, a seasonal periodic term, and a residual term; superimpose the long-term trend term and the seasonal periodic term to obtain the background water level change benchmark sequence of the spatial grid. This sequence describes how the water level at this point should have changed if there were no sudden events. It should be noted that the final water level response prediction field is formed by superimposing the pressure disturbance caused by this heavy rainfall event (converted into water level change through pressure-water level conversion coefficient) with the natural background water level at the same time. Only by superimposing the pressure field state tensor and the background water level digital field at the same future prediction time point can the true state of the groundwater level affected by the two factors at that time be accurately reflected. Therefore, ensuring that the two are generated based on the same future prediction time point sequence is the logical premise for achieving accurate fusion in a physical sense. In this invention, the historical period refers to an independent time interval preceding the current monitoring period, used to extract the long-term evolution patterns of the groundwater system. This historical period is strictly distinguished from the current monitoring period on the time axis and does not overlap with it. Its data is specifically used for trend decomposition in step four to obtain the background water level digital field. Specifically, the gridded groundwater level data of the historical period comes from the regular grid data covering the target area generated in step one, but its time range is strictly limited to a long time series preceding the start of the current monitoring period. Generally, complete hydrological annual data of 3 to 5 consecutive years are selected. This time length is chosen based on the following considerations: it is sufficient to cover multiple complete seasonal cycles to accurately separate the periodic term, while containing sufficient interannual variation to fit a reliable long-term trend term. Since the historical period data is strictly preceding the current monitoring period, it does not contain any information about the heavy rainfall event that needs to be predicted, nor does it contain any impact of the event on the groundwater system. To align with future forecast periods, the background water level change benchmark sequence needs to be extrapolated. The future forecast time point sequence is a series of moments that are located after the end of the current monitoring period and are evenly arranged according to a preset time step. The length and step size of this sequence must be strictly consistent with the length and time resolution of the pressure field state tensor sequence output by the feedforward neural network model in step 3. Extrapolating the background water level change benchmark sequence of each spatial grid to this series of future forecast time points yields the background water level prediction value of each spatial grid at each future forecast time point. Finally, the prediction values ​​of all spatial grids at the same future forecast time point are reassembled according to the original spatial grid row and column structure to obtain the two-dimensional background water level digital field corresponding to that future forecast time point. Arranging the background water level digital fields corresponding to all future forecast time points in chronological order constitutes a background water level digital field sequence that strictly corresponds to the pressure field state tensor sequence in the time dimension. Tensor addition operations are performed between the pressure field state tensor and the background water level digital field to synthesize the final water level response prediction field. The core of this approach is the introduction of a key parameter: the pressure-water level conversion coefficient. This coefficient is a scalar determined based on the physical characteristics of the aquifer, and its physical meaning is equivalent to the reciprocal of the aquifer's specific yield or storage coefficient. It establishes a linear conversion relationship between the virtual pressure increment and the actual water level change. The specific calculation process unfolds for each common future prediction time point: the pressure-water level conversion coefficient is determined based on the aquifer's physical characteristics; for each future prediction time point, the pressure value of each spatial grid corresponding to that future prediction time point in the pressure field state tensor is multiplied by the pressure-water level conversion coefficient to obtain the water level change field; this water level change field is then added grid by grid to the same spatial grid value corresponding to the same future prediction time point in the background water level digital field to obtain the water level response prediction field for that future prediction time point. This operation is repeated for all future prediction time points to obtain a sequence of water level response prediction fields arranged in chronological order. Each water level response prediction field in this sequence integrates the natural background changes and the abnormal response triggered by the current event, jointly reflecting the complete process of the dynamic evolution of groundwater level in the future period. The pressure-water level conversion coefficient is physically approximately equal to the reciprocal of the aquifer specific yield or related to the storage coefficient. For example, for a typical unconfined aquifer with a specific yield of about 0.1 to 0.2, the coefficient is about 5 to 10 (meter head per meter water level). The specific value needs to be determined through hydrogeological tests (such as pumping tests) or historical data inversion in the study area. This coefficient is the core bridge connecting the virtual pressure field and the actual water level changes. The method extracts the water level change prediction curves of preset prediction points from the water level response prediction field sequence, and outputs the groundwater level rise rate after differential operation. Specifically, it extracts the water level prediction values ​​of pre-selected groundwater prediction points at each future prediction time point from the water level response prediction field sequence arranged in chronological order, connects these water level prediction values ​​in chronological order to form a prediction curve of water level changing with time, and calculates the instantaneous rate of change of the prediction curve at each future prediction time point using a numerical differential method (such as the central difference method). This instantaneous rate of change is the groundwater level rise rate and is output as the final result. A positive value indicates a rise in water level and a negative value indicates a fall in water level. Specifically, for a discrete water level prediction sequence, at the target future prediction time point, the instantaneous slope of that point is calculated using the water level prediction values ​​of adjacent time points. This slope directly reflects the dynamic change rate of the groundwater level at that moment and is a core indicator for emergency response decisions, such as the assessment of urban flooding risk.

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

[0024] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0025] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0026] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for predicting groundwater rise under heavy rainfall conditions, characterized in that, The specific steps include: Step 1: Obtain rainfall intensity data and groundwater level data for the target area, and apply a spatial grid interpolation algorithm to convert the groundwater level data into regular grid data with the same spatial resolution as the rainfall intensity data, generating a spatiotemporal matrix containing rainfall intensity channels and groundwater level channels; Step 2: Continuously scan the rainfall intensity data in the spatiotemporal matrix, use the threshold detection method to identify continuous periods exceeding the set heavy rainfall threshold and record them as heavy rainfall events. For each heavy rainfall event, calculate the information entropy of its rainfall intensity distribution, determine the directional offset vector of its rainfall intensity weighted centroid relative to the preset prediction point, and calculate the integral cumulative value of rainfall intensity over time during the duration of the event. Combine the information entropy, directional offset vector and integral cumulative value to construct a multidimensional feature vector. Step 3: Input the multidimensional feature vector into the pre-trained feedforward neural network model. The model performs numerical processing based on the internal parameterized decay function and anisotropic diffusion convolution kernel to calculate the pressure field state tensor of the target region. Step 4: Perform trend decomposition on historical groundwater level data to obtain the background water level digital field. Perform tensor addition operation on the pressure field state tensor and the background water level digital field to obtain the water level response prediction field. Extract the water level change prediction curve of the preset prediction point from the water level response prediction field, and output the groundwater level rise rate after differentiation operation.

2. The groundwater rise prediction method based on heavy rainfall conditions according to claim 1, characterized in that: The groundwater level data is converted into regular grid data with the same spatial resolution as the rainfall intensity data. Specifically, the sampling interval between adjacent time sampling points is preset, and multiple groundwater level measurement points with different geographical coordinates are selected in the target area to obtain the groundwater level data of each groundwater level measurement point at each time sampling point within the current monitoring period. At the same time, the rainfall intensity data covering the same target area within the same monitoring period is obtained. This rainfall intensity data has a unified regular geographical grid and a fixed spatial resolution, and each time sampling point has corresponding gridded rainfall intensity data. Based on the regular geographic grid used for rainfall intensity data, for each time sampling point, the Kriging spatial interpolation algorithm is used to calculate and generate the corresponding groundwater level estimate for each spatial grid in the regular geographic grid, taking the groundwater level data and its geographic coordinates at that time sampling point as input. Through the above operation, the originally discrete groundwater level data is transformed into gridded groundwater level data that is aligned with the rainfall intensity data on the spatial grid framework and synchronized with it in time. A spatiotemporal matrix containing rainfall intensity channels and groundwater level channels is generated. Specifically, for the gridded rainfall intensity data and gridded groundwater level data corresponding to each time sampling point within the monitoring period, they are registered according to the same spatial grid coordinates. The two are superimposed in the data channel dimension to form a composite data layer with dual channels. The composite data layers corresponding to all time sampling points are arranged sequentially in chronological order to construct a three-dimensional spatiotemporal matrix. The first and second dimensions of this spatiotemporal matrix correspond to the row index and column index of the spatial grid, respectively, and the third dimension corresponds to the time sampling point index. The rainfall intensity data and the groundwater level data are stored as independent channel data in this spatiotemporal matrix.

3. The groundwater rise prediction method based on heavy rainfall conditions according to claim 2, characterized in that: A threshold detection method is used to identify continuous periods exceeding a set heavy rainfall threshold and record them as heavy rainfall events. Specifically, the rainfall intensity data in the spatiotemporal matrix is ​​scanned sequentially along the time dimension. For each time sampling point, the average rainfall intensity value of all spatial grids in the target area corresponding to that time sampling point is calculated. The average rainfall intensity value of each time sampling point is continuously compared with a preset heavy rainfall threshold to identify all time sampling points whose average rainfall intensity value exceeds the heavy rainfall threshold. Among the identified time sampling points, a series of time sampling points that are temporally continuous and whose time interval between two adjacent time sampling points does not exceed the preset maximum allowable interval are merged into an independent heavy rainfall event. At the same time, the start time and end time of each heavy rainfall event are recorded. For any heavy rainfall event, the information entropy of its rainfall intensity distribution is calculated as follows: the time sampling point when the rainfall intensity reaches the global maximum value during the duration of the heavy rainfall event is taken as the peak time; the rainfall intensity values ​​of all spatial grids in the target area at the peak time are extracted, and these rainfall intensity values ​​are normalized so that their sum is 1, thus forming a spatial probability distribution representing the spatial distribution of rainfall intensity; based on the Shannon entropy calculation formula in information theory, the information entropy of this spatial probability distribution is calculated, and the resulting entropy value is the information entropy index representing the heavy rainfall event. For any heavy rainfall event, the directional offset vector of its weighted centroid of rainfall intensity relative to the preset prediction point is determined. Specifically, the peak time of the heavy rainfall event is extracted, and the average rainfall intensity value of each spatial grid at the peak time is used as the weight. The weighted average values ​​of all spatial grids in the east-west and north-south directions are calculated respectively. The coordinate point formed by these two weighted average values ​​is the coordinate of the weighted centroid of rainfall intensity at the peak time of the heavy rainfall event. The plane vector pointing from the coordinate of the weighted centroid of rainfall intensity to the coordinate of the preset prediction point is calculated. The plane vector is normalized to make its magnitude 1, resulting in a unit vector that only represents directional information. The resulting unit vector is the directional offset vector.

4. The groundwater rise prediction method based on heavy rainfall conditions according to claim 3, characterized in that: For any heavy rainfall event, the cumulative integral of rainfall intensity over time is calculated over its duration. Specifically, for each time sampling point within the entire duration, the sum of rainfall intensities of all spatial grids within the target area at that time sampling point is calculated to obtain the total instantaneous rainfall intensity of the region at that time sampling point. The total instantaneous rainfall intensity of the region at all time sampling points throughout the entire duration of the heavy rainfall event is summarized to construct a time series of the total instantaneous rainfall intensity of the region. The time series of the total instantaneous rainfall intensity of the region is then integrated over the entire duration, and the scalar result obtained from the integration is the cumulative integral value. Information entropy, direction offset vector, and integral cumulative value are combined to construct a multidimensional feature vector. Specifically, the information entropy, two component values ​​of the direction offset vector, and integral cumulative value calculated for the same heavy rainfall event, a total of four numerical features, are arranged and combined into a one-dimensional numerical array in a preset fixed order to form a multidimensional feature vector.

5. The groundwater rise prediction method based on heavy rainfall conditions according to claim 4, characterized in that: The feedforward neural network model performs numerical processing based on an internal parameterized decay function and anisotropic diffusion convolution kernel. Specifically, the feedforward neural network model designs two computational components in its network structure, including a parameterized decay function component and an anisotropic diffusion convolution kernel component. The parameterized decay function component is parameterized based on the integral accumulation value in the multidimensional feature vector. Its function is to define a decay function that follows an exponential law to simulate the characteristic that the intensity of the virtual pressure field induced by rainfall infiltration decays over time. The anisotropic diffusion convolution kernel component is a two-dimensional weight matrix whose internal parameters are constrained during the training of the feedforward neural network model. This constraint enables the convolution kernel to have anisotropic spatial diffusion coefficients. Specifically, it has a dominant diffusion coefficient in the direction parallel to the direction offset vector in the multidimensional feature vector, and a different diffusion coefficient in another dominant direction perpendicular to this direction, thereby simulating the non-uniform transmission phenomenon of the virtual pressure field in space along the direction offset vector. The multidimensional feature vector is input into the pre-trained feedforward neural network model. The feedforward neural network model simulates the dynamic process of the virtual pressure field propagating from the end of the current monitoring period to the future through the collaborative computation of two components. Finally, it outputs the pressure field state tensor arranged in the order of the future predicted time points.

6. The groundwater rise prediction method based on heavy rainfall conditions according to claim 5, characterized in that: To obtain a background water level digital field, historical groundwater level data is decomposed into trends. Specifically, long-term time-series water level data corresponding to each spatial grid in the target area is extracted from historical gridded groundwater level data. A time series decomposition algorithm is applied to the long-term time-series water level data of each spatial grid to decompose it into a long-term trend term and a seasonal periodic term. The two are superimposed to obtain a background water level change benchmark sequence for each spatial grid. The background water level change benchmark sequence of each spatial grid is extrapolated to these future prediction time points to generate the background water level value of each spatial grid at each future prediction time point. Finally, the background water level values ​​of all spatial grids at the same future prediction time point are reorganized according to the spatial grid structure to obtain a background water level digital field corresponding to the pressure field state tensor in the time dimension. The water level response prediction field is obtained by performing tensor addition operations on the pressure field state tensor and the background water level digital field. Specifically, the pressure-water level conversion coefficient is determined based on the physical characteristics of the aquifer. For each future prediction time point, the pressure value of each spatial grid of the pressure field state tensor corresponding to that future prediction time point is multiplied by the pressure-water level conversion coefficient to obtain the water level change field. The water level change field is then added grid by grid to the same spatial grid value of the background water level digital field corresponding to the same future prediction time point to obtain the water level response prediction field for that future prediction time point. This operation is repeated for all future prediction time points to obtain a sequence of water level response prediction fields arranged in chronological order.

7. The groundwater rise prediction method based on heavy rainfall conditions according to claim 6, characterized in that: The water level change prediction curves of preset prediction points are extracted from the water level response prediction field, and the groundwater level rise rate is output after differential operation. Specifically, from the water level response prediction field sequence arranged in the order of future prediction time points, the water level prediction values ​​of the pre-selected groundwater prediction points at each future prediction time point are extracted in time order. These water level prediction values ​​are connected in time order to form a prediction curve of water level changing with time. The average slope of the prediction curve is calculated using the numerical differential method. The average slope value is the groundwater level rise rate and is output as the final result.

Images