Multi-scale fusion ecological hydrological interaction quantification method

Abstract

The invention discloses a multi-scale fusion ecological hydrological interaction quantification method, and relates to the technical field of ecological hydrology, and the method comprises the steps: 1, collecting environment driving data and remote sensing data; 2, preprocessing the remote sensing data, and carrying out hydrological process dynamic monitoring and ecological parameter collaborative inversion in combination with environment driving data to respectively obtain hydrological data and ecological indexes; 3, constructing a spatio-temporal data set by using the hydrological data, the ecological indexes and the environment driving data; 4, constructing a bidirectional LSTM neural network model, and performing training optimization by using the spatio-temporal data set to obtain a water volume change predicted value, lag time and corresponding ecological variables; 5, lag effect analysis is carried out according to the water volume change predicted value, the lag time and the corresponding ecological variables, and an interaction quantification network diagram of the ecological hydrological elements is constructed. According to the method, the dynamic coupling relation of the ecological hydrological process of the sand lake basin is quantified, and theoretical support is provided for resource optimization management and ecological restoration engineering.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of ecological hydrology, and more particularly to a multi-scale integrated ecological hydrological interaction quantification method. BACKGROUND

[0002] In the current research on sand lake basins, traditional hydrological observation stations are sparse, ecological parameters are obtained with poor timeliness, and existing coupling models have certain defects, including the use of linear regression and other simplified methods for hydrological-ecological processes, which cannot represent nonlinear interactions, single spatiotemporal scale analysis leading to incomplete expression of system dynamic response mechanisms, and lack of a quantitative framework integrating remote sensing big data and mechanism models.

[0003] Therefore, how to quantitatively represent the hydrological interaction relationship of sand lake basins is a problem that needs to be solved by those skilled in the art. SUMMARY

[0004] Therefore, the present application provides a multi-scale integrated ecological hydrological interaction quantification method, which realizes the quantitative representation of the dynamic coupling relationship of ecological hydrological processes in sand lake basins through multi-source remote sensing data interpretation and machine learning modeling, and provides theoretical support for the optimal management of sand basin water resources and the accurate implementation of ecological restoration projects.

[0005] To achieve the above purpose, the present application adopts the following technical solutions:

[0006] A multi-scale integrated ecological hydrological interaction quantification method, comprising the following steps:

[0007] Step 1: Collect environmental driving data and remote sensing data;

[0008] Step 2: Preprocess the remote sensing data and combine the environmental driving data to dynamically monitor the hydrological process and obtain hydrological data; and cooperatively invert the ecological parameters of the preprocessed remote sensing data to obtain an ecological index;

[0009] Step 3: Construct a spatiotemporal data set using the hydrological data, the ecological index, and the environmental driving data;

[0010] Step 4: Construct a bidirectional LSTM neural network model and perform training and optimization, input the spatiotemporal data set into the trained bidirectional LSTM neural network model, and obtain a water quantity change prediction value, a lag time, and a corresponding ecological variable;

[0011] Step 5: Construct an interaction quantification network graph according to the water quantity change prediction value, the lag time, and the corresponding ecological variable.

[0012] Preferably, the environmental driving data includes precipitation, air temperature, lake water level, etc.; the remote sensing data includes Landsat-8 OLI multispectral image, Sentinel-2 MSI multispectral image, MODIS multispectral image, ERA5-Land reanalysis data and multi-source DEM data, preferably the image in the vegetation growth period; the multi-source DEM data includes spaceborne laser altimetry data and radar interferometric DEM data.

[0013] Preferably, the hydrological data includes monthly lake area and lake water demand; the process of dynamic monitoring of hydrological process is:

[0014] Step 211: extracting the lake boundary according to the multispectral image in the remote sensing data, and calculating the monthly lake area;

[0015] Step 212: constructing the area-storage capacity curve according to the multi-source DEM data of the remote sensing data;

[0016] Step 213: calculating the storage capacity according to the monthly lake area and the area-storage capacity curve.

[0017] Preferably, the specific process of step 211 is:

[0018] Step 2111: pre-processing the multispectral image, including radiation calibration, atmospheric correction, cloud mask processing and image fusion;

[0019] Step 2112: extracting the initial lake boundary grid data from the pre-processed multispectral image by using the pre-trained U-Net++ model and super-resolution reconstruction method in turn;

[0020] Step 2113: post-processing and optimizing the initial lake boundary grid data to obtain accurate lake boundary grid data; the process of post-processing and optimization is: performing open operation on the initial lake boundary grid data to obtain denoising grid data, using the global surface water body distribution dataset to correct the denoising grid data in time and space consistency to obtain accurate lake boundary grid data;

[0021] Step 2114: performing grid statistics according to the accurate lake boundary grid data, calculating the lake area to obtain the monthly lake area; performing grid statistics on the accurate lake boundary grid data based on the WGS84 coordinate system, and converting the actual area according to the statistical result to the pixel resolution as the lake area to obtain the monthly lake area.

[0022] Preferably, the specific process of step 212 is:

[0023] Step 2121: fusing the multi-source DEM data in the remote sensing data to obtain the fused DEM; specifically including:

[0024] The multi-source DEM data is converted to WGS84 coordinate system by affine transformation method and TPS interpolation method, and the converted multi-source DEM data is fused by weighted average based on precision weight to obtain the fused DEM.

[0025] In step 2122, the adaptive piecewise polynomial fitting is performed by using the fused DEM and the accurate lake boundary grid data to obtain the sub-interval elevation-area relationship and the piecewise polynomial parameters; specifically including:

[0026] The k sub-intervals are divided by Jenks natural breaking according to the accurate lake boundary grid data, and the sub-interval elevation-area relationship is fitted by adaptive piecewise polynomial fitting; the polynomial coefficients of the sub-interval elevation-area relationship are globally optimized by Levenberg-Marquardt algorithm to obtain the piecewise polynomial parameters.

[0027] In step 2123, the numerical integration is performed on the piecewise polynomial parameters by adaptive Simpson integration method to construct the area-storage curve.

[0028] Preferably, in step 213, the storage capacity is calculated according to the monthly lake area by using the area-storage curve; and the storage capacity change trend is obtained by combining the environmental driving data analysis.

[0029] Preferably, the specific process of the ecological parameter collaborative inversion is as follows:

[0030] In step 221, the normalized vegetation index NDVI is calculated according to the remote sensing data.

[0031] In step 222, the gross primary productivity GPP is calculated according to the remote sensing data.

[0032] In step 223, the normalized vegetation index NDVI and the gross primary productivity GPP constitute the ecological index.

[0033] Preferably, the specific process of step 221 is as follows:

[0034] In step 2211, the near-infrared band surface reflectance and the red band surface reflectance of the pre-processed remote sensing data are extracted.

[0035] In step 2212, the normalized vegetation index NDVI is calculated according to the near-infrared band surface reflectance and the red band surface reflectance.

[0036] Preferably, the specific process of step 222 is as follows:

[0037] In step 2221, the photosynthetically active radiation absorption ratio FPAR is obtained according to the remote sensing data, and the photosynthetically active radiation PAR, the air temperature T and the soil moisture W are obtained according to the reanalysis data set.

[0038] Step 2222: calculating a temperature stress factor according to the air temperature T, and calculating a water stress factor according to the soil humidity W;

[0039] Step 2223: calculating the total primary productivity GPP by using a light energy utilization model according to the photosynthetically active radiation absorption ratio FPAR, the photosynthetically active radiation PAR, the temperature stress factor, and the water stress factor.

[0040] Preferably, the monthly lake area, water storage, normalized difference vegetation index NDVI, total primary productivity GPP, and environmental driving data are standardized by using a Z-score normalization method to construct a spatiotemporal dataset.

[0041] Preferably, the bidirectional LSTM neural network model comprises an input layer, a bidirectional LSTM layer, a fully connected layer, and an output layer connected in sequence; the bidirectional LSTM layer comprises a forward LSTM unit, a backward LSTM unit, and a splicing unit; the output layer solves a water loss function and an ecological lag loss function, constructs a comprehensive loss function by weighted summation, and calculates a comprehensive loss to optimize the bidirectional LSTM neural network model by back propagation according to the comprehensive loss function.

[0042] Preferably, the training process of the bidirectional LSTM neural network model comprises:

[0043] Step 41: calculating hydrological features, ecological features, and environmental features according to the spatiotemporal dataset to form an input feature sequence; the hydrological features comprise a moving average water volume and a water level change rate; the ecological features comprise an ecological index variation coefficient and an ecological index spatial gradient; the environmental features comprise a SPEI drought index;

[0044] Step 42: transmitting the input feature sequence to the bidirectional LSTM layer by the input layer;

[0045] Step 43: extracting a historical time series dependency relationship by the forward LSTM unit, extracting a future potential impact feature by the backward LSTM unit, and splicing and merging the historical time series dependency relationship and the future potential impact feature by the splicing unit to obtain a comprehensive feature vector containing historical and future information;

[0046] Step 44: mapping the comprehensive feature vector to a water volume change prediction value, a lag time, and a corresponding ecological variable by the fully connected layer;

[0047] Step 45: calculating the water loss function according to the water volume change prediction value by the output layer, calculating the ecological lag loss function according to the lag time and the corresponding ecological variable, constructing a comprehensive loss function by weighted summation of the water loss function and the ecological lag loss function, calculating a comprehensive loss, and globally searching for a minimum comprehensive loss by using a Bayesian optimization algorithm to search for model hyperparameters and optimize the model.

[0048] Preferably, in step 5, the characteristic contribution is estimated according to the predicted value of water volume change, the lag time and the corresponding ecological variables to obtain the SHAP value; the weight of the eco-hydrological element is calculated according to the SHAP value and the lag time, and an interactive quantitative network diagram is constructed; the eco-hydrological elements include hydrological elements, ecological elements and environmental elements, and the hydrological elements, ecological elements and environmental elements in the interactive quantitative network diagram are used as network nodes, and the weights are marked between the network nodes.

[0049] It can be seen from the above technical solution that compared with the existing technology, the present invention discloses a multi-scale integrated eco-hydrological interaction quantification method, which constructs a spatiotemporal data set through dynamic monitoring of hydrological processes and collaborative inversion of ecological parameters, and uses the spatiotemporal data set to train a machine learning model to learn the eco-hydrological lag effect of the sandy lake basin, analyze the interaction relationship, and obtain an interaction quantification network diagram. According to the network diagram, it can provide data support for resource management optimization and ecological restoration in the sandy lake basin, and achieve more accurate ecological management. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0051] Figure 1 This is a flow chart of a multi-scale integrated eco-hydrological interaction quantification method provided by the present invention. DETAILED DESCRIPTION

[0052] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0053] The embodiment of the present invention discloses a multi-scale integrated eco-hydrological interaction quantification method, such as Figure 1 As shown, the following steps are included:

[0054] S1: Collect environmental driving data and remote sensing data;

[0055] S2: Preprocess the remote sensing data and combine it with environmental driving data to dynamically monitor the hydrological process and obtain hydrological data; perform collaborative inversion of ecological parameters on the preprocessed remote sensing data to obtain ecological indices;

[0056] S3: Constructing a spatio-temporal dataset using hydrological data, ecological index and environmental driving data;

[0057] S4: Constructing a bidirectional LSTM neural network model and training and optimizing it using the spatio-temporal dataset to obtain water quantity change prediction values, lag time and corresponding ecological variables;

[0058] S5: Constructing an interaction quantification network diagram according to the water quantity change prediction values, lag time and corresponding ecological variables.

[0059] Further, the environmental driving data includes precipitation, air temperature, lake water level and the like; the remote sensing data includes Landsat-8 OLI multispectral images, Sentinel-2 MSI multispectral images, MODIS multispectral images, ERA5-Land reanalysis data and multi-source DEM data, and the images of the vegetation growth period are preferably selected; the multi-source DEM data includes spaceborne laser altimetry data and radar interferometric DEM data.

[0060] Further, the hydrological data includes monthly lake area and storage capacity; the ecological index includes normalized vegetation index and gross primary productivity.

[0061] In one specific embodiment, the environmental driving data such as precipitation, air temperature and lake water level in the study area is collected; multi-temporal remote sensing images (Landsat-8, Sentinel-2) are obtained using the GEE platform, MODIS data remote sensing images are collected, and multi-source DEM data is also collected to form remote sensing data; the hydrological process is dynamically monitored by interpreting the lake boundary through a deep learning model (U-Net architecture) to calculate the monthly lake area; a high-precision fused DEM is constructed, and an adaptive segmented polynomial fitting is used to generate an area-storage capacity curve; based on the lake area estimation results and the area-storage capacity curve, the dynamic change of lake storage capacity is calculated; the specific process is as follows:

[0062] S211: Lake boundary interpretation and area calculation;

[0063] S2111: Remote sensing data acquisition and preprocessing;

[0064] Landsat-8 OLI (30-meter resolution) and Sentinel-2 MSI (10-meter resolution) multispectral images, and MODIS data (for cloud mask completion) are collected, and the images of the vegetation growth period (May-October) are preferably selected to ensure that at least one scene (one image covering a specific geographic area taken by a satellite during one pass) of effective data with low cloud cover (i.e. cloud interference is removed through cloud mask processing (such as QA band analysis, FMask algorithm) to ensure that the image can be used for lake boundary extraction) and reasonable preprocessing (i.e. completing radiation calibration, atmospheric correction, image fusion and the like to meet the analysis accuracy requirements) per month;

[0065] Data preprocessing is performed, including radiation calibration (converting DN values to surface reflectance), atmospheric correction (eliminating aerosol effects using FLAASH (Landsat) or Sen2Cor (Sentinel-2)), cloud mask completion processing (generating a mask based on the QA band of MODIS data and a cloud detection algorithm (FMask), if a certain area of Landsat-8 / Sentinel-2 is contaminated by clouds, then use the cloud-free image of MODIS on the same day or adjacent day (≤3 days) to replace the area), image fusion (using Gram-Schmidt Pan Sharpening method to improve the spatial resolution of multi-source data), etc.

[0066] S2112: Lake boundary extraction;

[0067] (a) Initial lake boundary raster data extraction;

[0068] The U-Net++ model (nested skip connection structure) is used, the backbone network of the U-Net++ model is ResNet-50, and the input channel number is 6 (blue, green, red, near-infrared, short-wave infrared 1, and short-wave infrared 2 bands); The training strategy of the U-Net++ model is: 200 labeled samples (boundary accuracy ± 0.5 pixels) of typical lakes in sandy land obtained through field investigation are used as the training set, synthetic data with sandy land specific noise (such as wind-sand shielding and salinization interference) are introduced for data enhancement, and the loss function (combination of DiceLoss and Boundary-aware Loss) is used to strengthen the boundary pixel weight; The pre-trained U-Net++ model is used to extract the preliminary lake boundary raster data from the preprocessed remote sensing data;

[0069] Then, for the mixed pixels (such as water-sand transition zone) in the preliminary lake boundary raster data, the super-resolution reconstruction technology (ESRGAN model) is used to realize sub-pixel level (1-3 meters) boundary positioning, and the initial lake boundary raster data is obtained;

[0070] (b) Post-processing optimization, obtaining high-precision, denoised accurate lake boundary raster data, the accurate lake boundary raster data is binary raster data (1 for water body, 0 for non-water body);

[0071] Perform open operation (erode first, then dilate) on the preliminary extracted initial lake boundary raster data, the kernel size is 3x3, and the iteration number is 1 time; Eliminate isolated noise regions with an area less than 10 pixels; Combine the JRC Global Surface Water dataset to modify the extraction results for spatial and temporal consistency, the specific steps are as follows:

[0072] First, based on the historical mean, the area outlier detection is performed, and the current extracted lake area At The deviation of A ht from the historical same period mean A t exceeds 20%, i.e., |A ht | / A ht >0.2, then the current extracted lake area is considered abnormal, and the abnormal value A t is replaced by the corrected value A c , and the corrected value A c is represented as:

[0073]

[0074] In the formula, A t-1 and A t+1 represent the effective areas of the adjacent previous and next months t-1 and t+1 respectively, and A c is the corrected lake area; the historical same period mean A ht is obtained through the JRC Global Surface Water dataset, and only the 10-year same month mean value is taken;

[0075] Secondly, small water bodies are removed based on NDWI verification; for lakes with an area <0.1 km 2 , the normalized water index (NDWI) is used for auxiliary verification to remove false detection areas with an NDWI mean value <0.2; the calculation formula of NDWI is as follows:

[0076] NDWI=(pGreen-pNIR) / (pGreen+pNIR)

[0077] In the formula, pGreen represents the surface reflectivity of the green band (Band 3 of Landsat-8 or Band 3 of Sentinel-2); and pNIR represents the surface reflectivity of the near-infrared band (Band 5 of Landsat-8 or Band 8 of Sentinel-2);

[0078] S2113: Lake area calculation and verification;

[0079] Based on the grid statistics under the WGS84 coordinate system, the actual lake area is converted into the actual lake area according to the pixel resolution, and the actual lake area is taken as the monthly lake area. The actual lake area calculation formula is as follows:

[0080] A=N×(p×cosθ) 2

[0081] Where A is the actual lake area (i.e., monthly lake area); N is the number of water body pixels, based on the accurate lake boundary raster data generated by post-processing optimization in S2112(b), the total number of pixels with a raster median value of 1 in the raster data is counted; p is the pixel size, based on the pre-processed Landsat-8 and Sentinel-2 multispectral images in S2111, determined by the spatial resolution of the multispectral images (Landsat-8 is 30 meters / pixel, and Sentinel-2 is 10 meters / pixel); and θ is the latitude band projection correction factor, extracted from the WGS84 image geographic coordinates converted from the MODIS data output in S2111, i.e., the latitude of the center point of the image;

[0082] Thirty GPS verification points (Trimble R10, planar accuracy 2 cm) were deployed in a typical lake for ground measurement verification, and the root mean square error (RMSE) was used as an index for model accuracy verification;

[0083]

[0084] In the formula, (x model,i ,y model,i ) is the coordinates of the i-th lake boundary point extracted from the remote sensing image, (x GPS,i ,y GPS,i ) is the coordinates of the i-th lake boundary point measured by GPS on the ground, and n is the sample size; the overall RMSE of the lake boundary is required to be ≤5 meters, and the local RMSE of the transition zone is required to be ≤3 meters, and the final output is a binary boundary raster;

[0085] S212: Constructing area-storage capacity curve;

[0086] S2121: Multi-source DEM data fusion;

[0087] Obtain multi-source DEM data such as spaceborne laser altimetry data (ICESat-2 ATL08, elevation accuracy ±0.1 m), radar interferometric DEM (TanDEM-X 90m, relative accuracy ±2 m); use affine transformation + TPS interpolation to unify all DEMs to the WGS84 coordinate system, with a planar error <0.5 pixels;

[0088] Generate a fused DEM based on the accuracy weight of the data source:

[0089] H fused =(w1H ICESat-2 +w2H TanDEM-X ) / (w1+w2)

[0090] In the formula, H fused is the fused DEM, i.e., the fused high-precision elevation data obtained by weighted averaging of multi-source data; and H ICESat-2The data is the onboard laser altimeter data of ICESat-2 satellite, with an elevation accuracy of ±0.1m; H TanDEM-X is the radar interferometry DEM data from the TanDEM-X satellite, with a relative accuracy of ±2 meters. w1 and w2 are the weights of different data sources, which are assigned in inverse proportion to the accuracy of each data source. For example, the accuracy of ICESat-2 is ±0.1 meters, so its weight w1 is 10. The accuracy of TanDEM-X is ±2 meters, so its weight w2 is 0.5. The essence of fused DEM is to integrate multiple data sources such as ICESat-2 (laser altimetry) and TanDEM-X (radar interferometry), leveraging their respective advantages (ICESat-2's high accuracy and TanDEM-X's wide coverage) to reduce the systematic errors of a single data source through weight assignment, making the generated fused DEM more reliable.

[0091] S2122: Adaptive piecewise polynomial fitting;

[0092] The lake area was defined as the geographic boundary using the precise lake boundary grid data extracted from S2112. The Jenks natural break method was used to divide the lake area into k subintervals (k = 3-5), ensuring that the TRI coefficient of variation within each subinterval was <15%. For each subinterval j, a cubic polynomial was fitted to establish the subinterval elevation-area relationship:

[0093] A j (z) = p j1 z 3 +p j2 z 2 +p j3 z+p j4

[0094] Among them, z is the elevation independent variable, according to the pixel elevation value H in the subinterval j of the fused DEM i Get; A j (z) is the cumulative lake area corresponding to the elevation variable z in subinterval j, obtained by fitting the monthly lake area; p j1 、p j2 、p j3 is a polynomial coefficient, which is determined by fitting the elevation-area data (Hi, A) in the subinterval j. It controls the slope and curvature of the fitting curve and is also called the nonlinear terrain morphological coefficient. It reflects the steepness and gentleness of the terrain and is used to characterize the terrain morphological characteristics of the interval. j4 is a constant term, i.e., the elevation datum offset;

[0095] In order to ensure that the area-elevation relationship is continuous and smooth, can truly reflect the lake topography, and guarantee the convergence of the area-reservoir capacity curve constructed later when integrating the reservoir capacity, constraints are designed for the piecewise polynomial parameter optimization process, that is, adjacent intervals are connected at the point z. j Satisfy Aj (z j )=A j+1 (z j ) and the first-order derivative is continuous, which directly determines the accuracy and reliability of the area-capacity curve;

[0096] The Levenberg-Marquardt algorithm is used to globally optimize the polynomial coefficients of the elevation-area relationship of all subintervals. The objective function is:

[0097]

[0098] Among them, A obs (z i ) is at a specific elevation z i The actual observed lake area is derived from the grid statistics of the precise lake boundary grid data; model (z i ) is the lake area predicted by the subinterval elevation-area relationship; λ is the smoothing penalty coefficient, which is used to balance the data fitting accuracy and the smoothness of the piecewise polynomial (the default value is λ = 0.1); is the elevation-area relationship polynomial of the j-th subinterval at the intersection point z j The first-order derivative (i.e., slope) at reflects the rate of change of elevation-area in this interval; The j+1th subinterval polynomial is at the same intersection point z j The first-order derivative at , ensures smooth connection between adjacent intervals;

[0099] Solve the objective function according to the constraints and obtain the global optimization parameter matrix of the optimization output, that is, a set of piecewise polynomial parameters that can simultaneously meet high-precision fitting and global smooth connection, which can be directly used for area-capacity curve integration;

[0100] S2123: Construct area-capacity curve;

[0101] The area-capacity curve is generated by numerical integration based on the piecewise polynomial parameters and is expressed as:

[0102]

[0103] The specific implementation uses the adaptive Simpson integration method, and the absolute error limit is set to 10 -4 m 3 ; Integrate the elevation-area polynomial of each subinterval, accumulate the integral results of all subintervals, and obtain the total reservoir capacity V(z). The integral result is expressed in the form of a piecewise function as the mathematical model of the area-reservoir capacity curve;

[0104] The ±1σ (σ = 0.5 m) Gaussian noise is added to the fused DEM elevation value, and the 95% confidence interval of the storage capacity is generated by repeating the integral calculation 1000 times. The influence of the uncertainty of the fused DEM on the area-storage curve is quantified by the noise and the confidence interval.

[0105] S213: calculating the dynamic lake water demand;

[0106] According to the sequence of monthly lake area, the monthly lake area A is converted into the corresponding storage capacity V through the area-storage curve; through the storage capacity time series curve, the Mann-Kendall trend test method is used to analyze the seasonal fluctuations and long-term trends of the storage capacity changes.

[0107] In one specific embodiment, the pre-processed (radiometric calibration, atmospheric correction) Landsat-8 multispectral image is obtained, and the ecological index is obtained through ecological parameter collaborative inversion. The specific process of ecological parameter collaborative inversion is as follows:

[0108] S221: calculating the normalized difference vegetation index NDVI using Landsat series multispectral images;

[0109] NDVI is used to quantify vegetation coverage and growth status, and its calculation depends on the specific band reflectance data of multispectral images;

[0110] The near-infrared band (NIR, Band 5 (0.85-0.88 μm) of Landsat-8) and the red band (Red, Band 4 (0.64-0.67 μm) of Landsat-8) of the pre-processed (radiometric calibration, atmospheric correction) Landsat-8 multispectral image are selected, and the NDVI calculation formula is used for calculation;

[0111] The NDVI calculation formula is:

[0112] NDVI = (ρ NIR -ρ Red ) / (ρ NIR +ρ Red )

[0113] ρ NIR represents the near-infrared surface reflectance; ρ Red represents the red band surface reflectance; the NDVI value range is [-1, 1], and the vegetation coverage area is usually > 0.2, and the water body or bare land is < 0;

[0114] S222: Inverting GPP raster data based on the light use efficiency model (PML_V2); GPP (gross primary productivity) in GPP raster data represents the total amount of carbon fixed by vegetation through photosynthesis, with the unit of gC / m 2 / day; the specific inversion process is as follows:

[0115] S2221: Obtain FPAR (fraction of photosynthetically active radiation) by MODIS multi-spectral image; obtain meteorological data including PAR (photosynthetically active radiation), air temperature (T) and soil moisture (W) data by ERA5-Land reanalysis data;

[0116] S2222: Calculate GPP by light use efficiency model (PML_V2) to obtain GPP grid data, which is expressed as:

[0117] GPP = FPAR · PAR · ∈max · f(T) · f(W)

[0118] Where ∈max represents the maximum light use efficiency corresponding to the vegetation type; f(T) represents the temperature stress factor (0-1), which is calculated based on air temperature T; f(W) represents the water stress factor (0-1), which is calculated based on soil moisture W; FPAR, PAR and stress factors are calculated pixel by pixel, ∈max is obtained by looking up the table combined with the vegetation type, and finally GPP grid data is output.

[0119] In one specific embodiment, based on the above hydrological process dynamic monitoring and ecological parameter collaborative inversion results, hydrological dimension (reservoir capacity of the lake, lake area), ecological index (NDVI, GPP time series data) and environmental driving data (precipitation, temperature, lake water level) are directly collected, Z-score normalization method is used for standardization processing of the data, dimension difference is eliminated, and a spatio-temporal data set is constructed; the normalized data is integrated into a spatio-temporal cube database, supporting dynamic analysis of multi-dimensional ecological hydrological interaction relationship.

[0120] In one specific embodiment, a bidirectional long short-term memory (LSTM) neural network model is constructed for analyzing the ecological hydrological lag effect of the sandy lake basin; the specific process of training the model is as follows:

[0121] S41: Prepare the input feature sequence of the model according to the spatio-temporal data set: contains hydrological features (moving average water volume, water level change rate), ecological features (ecological index variation coefficient, spatial gradient), and environmental features (SPEI drought index);

[0122] ① Moving average water volume

[0123] Based on the time series data of monthly lake water storage, the moving average is applied to obtain the moving average water volume V' t , which is expressed as:

[0124]

[0125] In the formula, V t-iHt-i represents the lake storage in the t-i th month, which is obtained according to the area-storage curve; n represents the size of the sliding window (such as 3 months, 6 months), which is used to smooth short-term fluctuations and reflect long-term trends;

[0126] ②Water level change rate

[0127] Based on the collected environmental driving data, monthly lake water level data are obtained, and the water level sequence is subjected to first-order difference to obtain the water level change rate ΔH t , the formula is

[0128] ΔH t = H t -H t-1

[0129] In the formula, H t represents the lake water level value in the t th month, and k represents the time interval (1 month);

[0130] ③Coefficient of variation of ecological index (Coefficient of Variation, CV)

[0131] Based on the average NDVI and GPP time series data of the basin, the coefficient of variation of the time series of the ecological index is calculated. Taking NDVI as an example, the calculation formula is:

[0132]

[0133] In the formula, CV NDVI represents the coefficient of variation of the NDVI time series data; σ NDVI represents the standard deviation of the NDVI time series data, reflecting the volatility of the time series; μ NDVI represents the mean of the NDVI time series data, reflecting the average ecological state; i represents the i th time point; and n represents the number of time series data points.

[0134] ④Ecological index spatial gradient

[0135] Based on the spatial distribution grid data of NDVI and GPP of the basin, the spatial change rate is calculated by using the spatial gradient operator (Sobel operator) to obtain the spatial gradient Gradient(x,y) of the ecological index. Taking NDVI as an example, the calculation formula is:

[0136]

[0137] In the formula, x is the longitude, and y is the dimension; represents the spatial gradient operator;

[0138] ⑤SPEI drought index

[0139] Based on the collected lake basin precipitation and temperature data in the environmental driving data, the calculation method of SPEI (Standardized Precipitation Evapotranspiration Index) is as follows: first, based on the temperature, the potential evapotranspiration is calculated by using the Penman formula, then the difference between the precipitation and the potential evapotranspiration is calculated, then the difference sequence is accumulated in a 3-month scale, and finally the SPEI value is obtained by standardizing the Log-Logistic probability distribution, so as to obtain the SPEI drought index;

[0140] S42: The constructed bidirectional LSTM neural network model comprises an input layer, a bidirectional LSTM layer, a full connection layer and an output layer;

[0141] ①Input layer: receiving the time series data of the above hydrological features, ecological features and environmental features, the dimension is T*N, wherein the time step T=12, and the feature number N=6;

[0142] ②Bidirectional LSTM layer: containing forward LSTM unit and reverse LSTM unit, respectively used for extracting historical time sequence dependence and future potential influence features; the output is merged by concatenation to form a comprehensive feature vector containing historical and future information;

[0143] ③Full connection layer: mapping the 128-dimensional comprehensive feature vector to two output targets (water quantity change prediction value, lag time and corresponding ecological variable, which can be used to judge the ecological lag effect through the lag time and the corresponding ecological variable).

[0144] ④Output layer: setting a comprehensive loss function to optimize the model, including:

[0145] (a) Water loss function: based on the water quantity change prediction error between the water quantity change prediction value and the measured value, the water loss function Loss1 is calculated by mean square error (MSE);

[0146]

[0147] (b) Ecological lag loss function: based on the cross-correlation analysis of the hydrological process and the lagged ecological process output by the model, the ecological variable is determined through the cross-correlation coefficient at different lag times, and the ecological lag loss function Loss2 is calculated;

[0148]

[0149] Wherein, R(τ) is the cross-correlation coefficient; ΔH t is the water level change rate, E t+τ is the ecological variable of the lag time τ; and are the mean values; T is the total time length; τ' is all the considered lag times; τ is the best lag time corresponding to the maximum cross-correlation coefficient;

[0150] (c) The comprehensive loss function is a weighted sum of Loss1 and Loss2:

[0151] Loss = Loss1 + λLoss2, λ = 1;

[0152] (d) Hyperparameter optimization: the hyperparameters of the above bidirectional LSTM model are globally searched using a Bayesian optimization algorithm, wherein the hyperparameters include the number of LSTM layers, the number of hidden units, the learning rate, and the sliding window length; the Bayesian optimization algorithm is based on a Gaussian process regression (GPR) model, and the search space is defined as: the number of LSTM layers: 1 to 3 layers, the number of hidden units: 32 to 128, the initial learning rate: 0.0001 to 0.01, and the sliding window length: 6 to 24 months; the optimization goal is to minimize the joint loss (Loss) of the validation set; the number of iterations is set to 100 times, and the Top 5 parameter combinations are selected for cross-validation;

[0153] (e) Model training;

[0154] The data set is divided into a training set (70%), a validation set (15%), and a test set (15%), and the time series data is input in batches (such as each batch containing 12 consecutive months of data), the weights are updated by back propagation, and the validation set loss is monitored; the total loss of the validation set is terminated when it does not decrease for 5 consecutive times.

[0155] In one specific embodiment, the process of interaction relationship analysis is:

[0156] S51: Feature contribution estimation: based on the trained model, the SHAP (Shapley Additional Explanation) value is used to quantify the feature contribution of each input feature to the model output, so as to analyze the driving factors of the ecological hysteresis effect;

[0157] The calculation method of the SHAP value is to perturb each input feature, and calculate the marginal contribution of the hysteresis effect explanation through the change of the model output; taking the water level change rate ΔH t in the hydrological feature as an example, the SHAP value calculation process is as follows:

[0158] ① Fix other features, and only perturb the water level change rate ΔH t ;

[0159] ② Calculate the difference in model output before and after adding ΔH t , f(S) is the ecological hysteresis response value without adding H, and f(S∪{H}) is the response value after adding ΔH t .

[0160] ③ Traverse all possible feature subsets S, and obtain φ H by weighted sum according to the Shapley value formula.

[0161] The calculation formula is:

[0162]

[0163] wherein φ Hi is the SHAP value of feature i; F is the set of all features; S is a subset that does not contain feature i, f(S) is the predicted value of the model when using the subset S, and the predicted value includes the water volume change prediction value; f(S∪{i}) is the predicted value after adding feature i;

[0164] S52: Constructing an interaction network topology;

[0165] According to the SHAP value and the lag time, an interaction network topology diagram of the ecological hydrological elements is generated, and a key feedback path is identified; the hydrological elements (H, including the water storage volume and the water level change rate), the ecological elements (E, including the NDVI and the GPP), and the environmental elements (Env, including the air temperature and the precipitation) are defined as network nodes; based on the SHAP value and the lag time τ, the weight of the edge is defined, for example, the edge weight ω H→E The assignment calculation formula is:

[0166]

[0167] wherein φ H is the SHAP value (contribution degree) of the hydrological element H to the ecological element E, which is obtained through a bidirectional LSTM model and SHAP analysis; τ is the optimal lag time, that is, the delay time of the hydrological change to the ecological response, which is determined through cross-correlation analysis;

[0168] Finally, the node colors are distinguished according to the element types, the edge widths are designed according to the weight sizes, the lag time τ and the weight value are labeled by adding edge labels, and finally a visual network diagram reflecting the complex feedback of the ecological hydrology is formed, which provides a reference for the basin management.

[0169] The various embodiments in the specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be mutually referred to. For the device disclosed by the embodiments, since it corresponds to the method disclosed by the embodiments, the description is relatively simple, and the related parts can be referred to the method part.

[0170] The foregoing description of the disclosed embodiments enables a person skilled in the art to make or use the application. Modifications of these embodiments will occur to persons of skill in the art, and that the appended claims are intended to cover all such modifications that do not depart from the true spirit and scope of the application. Therefore, the application is not limited to the embodiments shown but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A multi-scale integrated eco-hydrological interaction quantification method, characterized by: The following steps are involved: Step 1: Collect environmental driving data and remote sensing data; Step 2: Preprocess the remote sensing data and combine it with environmental driving data to dynamically monitor the hydrological process and obtain hydrological data; The ecological parameters are collaboratively inverted on the pre-processed remote sensing data to obtain the ecological index; Step 3: Construct a spatiotemporal dataset using hydrological data, ecological indices, and environmental driver data; Step 4: Build a bidirectional LSTM neural network model and perform training and optimization. Input the spatiotemporal dataset into the trained bidirectional LSTM neural network model to obtain the predicted value of water volume change, lag time, and corresponding ecological variables. Step 5: Construct an interaction quantitative network diagram based on the predicted water volume change value, lag time, and corresponding ecological variables.

2. A multi-scale fusion eco-hydrological interaction quantification method according to claim 1, characterized in that: Hydrological data include monthly lake area and water storage; the process of dynamic monitoring of hydrological processes is as follows: Step 211: extracting lake boundaries based on remote sensing data and calculating monthly lake areas; Step 212: constructing an area-reservoir capacity curve based on remote sensing data; Step 213: Calculate the water storage capacity based on the monthly lake area and the area-storage capacity curve.

3. The multi-scale fusion eco-hydrological interaction quantification method according to claim 2 is characterized in that: The specific process of step 211 is: Step 2111: pre-processing remote sensing data; Step 2112: extracting initial lake boundary raster data from the preprocessed remote sensing data using the pre-trained U-Net++ model and super-resolution reconstruction method; Step 2113: performing post-processing optimization on the initial lake boundary raster data to obtain accurate lake boundary raster data; the post-processing optimization includes opening operation and spatiotemporal consistency correction; Step 2114: Perform grid statistics based on the precise lake boundary grid data, calculate the lake area, and obtain the monthly lake area.

4. A multi-scale fusion eco-hydrological interaction quantification method according to claim 3, characterized in that: The specific process of step 212 is: Step 2121: Perform weighted averaging fusion on the remote sensing data to obtain a fused DEM; Step 2122: Using the fused DEM and the precise lake boundary raster data, perform adaptive piecewise polynomial fitting to obtain the sub-interval elevation-area relationship and piecewise polynomial parameters; specifically, the following steps are involved: Based on the precise lake boundary grid data, Jenks natural fractures were used to divide the lake into k sub-intervals. The elevation-area relationship of the sub-intervals was fitted using an adaptive piecewise polynomial. The polynomial coefficients of the sub-interval elevation-area relationship were globally optimized using the Levenberg-Marquardt algorithm to obtain the piecewise polynomial parameters. Step 2123: Adopt the adaptive Simpson integration method to perform numerical integration on the piecewise polynomial parameters and construct the area-capacity curve.

5. The multi-scale integrated eco-hydrological interaction quantification method according to claim 2 is characterized in that: The specific process of collaborative inversion of ecological parameters is as follows: Step 221: Calculating the Normalized Difference Vegetation Index (NDVI) based on remote sensing data; Step 222: Calculate gross primary productivity (GPP) based on remote sensing data; Step 223: The normalized difference vegetation index NDVI and the gross primary productivity GPP are combined to form an ecological index.

6. The multi-scale fusion eco-hydrological interaction quantification method according to claim 5 is characterized in that: The specific process of step 221 is: Step 2211: extracting the surface reflectance of the near infrared band and the surface reflectance of the red light band of the pre-processed remote sensing data; Step 2212: Calculate the Normalized Difference Vegetation Index (NDVI) based on the surface reflectance in the near-infrared band and the surface reflectance in the red band.

7. The multi-scale integrated eco-hydrological interaction quantification method according to claim 5 is characterized in that: The specific process of step 222 is: Step 2221: Obtain photosynthetically active radiation absorption ratio based on remote sensing data, and obtain photosynthetically active radiation, air temperature, and soil moisture based on the reanalysis dataset; Step 2222: Calculate the temperature stress factor according to the air temperature; Calculate water stress factor based on soil moisture; Step 2223: Calculate the gross primary productivity (GPP) using a light energy utilization rate model based on the photosynthetically active radiation absorption ratio, photosynthetically active radiation, temperature stress factor, and water stress factor.

8. The multi-scale integrated eco-hydrological interaction quantification method according to claim 1 is characterized in that: The Z-score normalization method was used to standardize the hydrological data, ecological index and environmental driving data to construct a spatiotemporal dataset.

9. The multi-scale integrated eco-hydrological interaction quantification method according to claim 1 is characterized in that: The bidirectional LSTM neural network model consists of a sequentially connected input layer, a bidirectional LSTM layer, a fully connected layer, and an output layer. The bidirectional LSTM layer includes a forward LSTM unit, a backward LSTM unit, and a splicing unit. The output layer solves the water loss function and the ecological hysteresis loss function, constructs a comprehensive loss function through weighted summation, and calculates the comprehensive loss based on the comprehensive loss function. Back propagation optimizes the bidirectional LSTM neural network model. The bidirectional LSTM neural network model is trained using historical spatiotemporal datasets. The training process includes: Step 41: Calculate hydrological characteristics, ecological characteristics, and environmental characteristics based on the spatiotemporal dataset to form an input feature sequence; Step 42: The input layer transmits the input feature sequence to the bidirectional LSTM layer; Step 43: The forward LSTM unit extracts historical temporal dependencies, the backward LSTM unit extracts future potential impact features, and the concatenation unit concatenates and merges the historical temporal dependencies and future potential impact features to obtain a comprehensive feature vector containing both historical and future information. Step 44: The fully connected layer maps the comprehensive feature vector to the predicted value of water volume change, lag time and corresponding ecological variables; Step 45: The output layer calculates and solves the water loss function based on the predicted value of water volume change, solves the ecological hysteresis loss function based on the lag time and the corresponding ecological variables, constructs a comprehensive loss function by weighted summation of the water loss function and the ecological hysteresis loss function, calculates the comprehensive loss, and uses the Bayesian optimization algorithm to perform a global search on the model hyperparameters to minimize the comprehensive loss and optimize the model.

10. The multi-scale integrated eco-hydrological interaction quantification method according to claim 2, characterized in that: In step 5, the characteristic contribution is estimated based on the predicted value of water volume change, lag time and corresponding ecological variables to obtain the SHAP value; the weights of the eco-hydrological elements are calculated based on the SHAP value and lag time, and an interaction quantitative network diagram is constructed; the eco-hydrological elements include hydrological elements, ecological elements and environmental elements. In the interaction quantitative network diagram, the hydrological elements, ecological elements and environmental elements are used as network nodes, and the weights are marked between the network nodes.

Images