A multi-scale fused ecological hydrological interaction quantification method

Abstract

The application discloses a kind of multi-scale fusion ecological hydrology interaction quantification method, it is related to the technical field of ecological hydrology, comprising: step 1: collection environment driving data and remote sensing data;Step 2: remote sensing data is preprocessed, and water hydrological process dynamic monitoring and ecological parameter collaborative inversion are carried out in combination with environment driving data, respectively obtain hydrological data and ecological index;Step 3: utilize hydrological data, ecological index and environment driving data to build space-time data set;Step 4: build bidirectional LSTM neural network model, and utilize space-time data set to train optimization, obtain water quantity change predicted value, lag time and corresponding ecological variable;Step 5: according to water quantity change predicted value, lag time and corresponding ecological variable carry out lag effect analysis, and construct the interaction quantification network diagram of ecological hydrological element.The application realizes sand lake basin ecological hydrological process dynamic coupling relationship quantification, provides theoretical support 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, and the image in the vegetation growth period is preferably selected; 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 processes 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: sequentially using the pre-trained U-Net++ model and super-resolution reconstruction method to extract the initial lake boundary grid data from the pre-processed multispectral image;

[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 an opening operation on the initial lake boundary grid data to obtain denoised grid data, and using the global surface water body distribution dataset to correct the denoised grid data for spatiotemporal consistency to obtain accurate lake boundary grid data;

[0021] Step 2114: performing grid statistics according to the accurate lake boundary grid data to calculate the lake area and 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 a 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 degree of each ecological variable is estimated according to the water quantity change prediction value, the lag time and the corresponding ecological variable, and the SHAP value is obtained; the weight of the eco-hydrological element is calculated according to the SHAP value and the lag time, and the interaction quantization network diagram is constructed; the eco-hydrological element includes a hydrological element, an ecological element and an environmental element, and the hydrological element, the ecological element and the environmental element in the interaction quantization network diagram are taken as network nodes, and the weights between the network nodes are marked.

[0049] According to the technical solution, compared with the prior art, the application provides a multi-scale fusion ecological hydrological interaction quantization method, a time-space data set is constructed through hydrological process dynamic monitoring and ecological parameter collaborative inversion, a machine learning model is trained using the time-space data set, and the ecological hydrological lag effect of the sandy lake basin is learned, the interaction relationship is analyzed, and an interaction quantization network diagram is obtained, which can provide data support for resource management optimization and ecological restoration of the sandy lake basin, and realize more accurate ecological management. BRIEF DESCRIPTION OF DRAWINGS

[0050] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or the prior art description will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and those skilled in the art can obtain other drawings according to the provided drawings without creative labor.

[0051] Figure 1 A multi-scale fusion ecological hydrological interaction quantization method flowchart is provided. DETAILED DESCRIPTION

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

[0053] The embodiments of the present application disclose a multi-scale fusion ecological hydrological interaction quantization method, as shown in Figure 1 The method comprises the following steps:

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

[0055] S2: pre-processing the remote sensing data, and combining the environmental driving data to perform hydrological process dynamic monitoring to obtain hydrological data; and performing ecological parameter collaborative inversion on the pre-processed remote sensing data to obtain an ecological index;

[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 in 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 by 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 result and the area-storage capacity curve, the dynamic change of the 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 in 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 in one pass) of effective data with low cloud cover (i.e. cloud interference is removed by 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) is obtained 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 elevation data is from the onboard laser altimetry of the ICESat-2 satellite, with an accuracy of ±0.1 meters; H TanDEM-X The data is derived from radar interferometry DEM data from the TanDEM-X satellite, with a relative accuracy of ±2 meters. w1 and w2 represent the weights of different data sources, allocated inversely proportional to their respective accuracy. For example, ICESat-2 has an accuracy of ±0.1 meters, so its weight w1 is 10; TanDEM-X has an accuracy of ±2 meters, so its weight w2 is 0.5. The essence of DEM fusion is to integrate multi-source data from ICESat-2 (laser altimeter) and TanDEM-X (radar interferometry), leveraging their respective advantages (ICESat-2's high accuracy and TanDEM-X's wide coverage) and reducing systematic errors from single data sources through weight allocation, thus making the generated DEM more reliable.

[0091] S2122: Adaptive piecewise polynomial fitting;

[0092] The precise lake boundary raster data extracted from S2112 was used as the geographic boundary to define the lake region. The Jenks natural faulting method was used to divide the lake region into k sub-intervals (k = 3–5), ensuring that the TRI coefficient of variation within each sub-interval was <15%. For each sub-interval j, a cubic polynomial was fitted to establish the sub-interval elevation-area relationship.

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

[0094] Where z is the elevation independent variable, which is determined by the pixel elevation value H in sub-interval j of the fused DEM. i A j (z) represents the cumulative lake area corresponding to the elevation variable z within sub-interval j, obtained by fitting the lake area monthly; p j1 p j2 p j3 These are polynomial coefficients, determined by fitting the elevation-area data (Hi, A) within the sub-interval j. They control the slope and curvature of the fitted curve and are also known as nonlinear terrain morphology coefficients. They reflect the steepness and gentleness of the terrain and are used to characterize the terrain morphology of that interval; p j4 This is a constant term, namely, the offset of the elevation datum.

[0095] To ensure a continuous and smooth area-elevation relationship that accurately reflects the lake's topography, and to guarantee the convergence of the subsequently constructed area-storage capacity curve during storage capacity integration, constraints are designed for the piecewise polynomial parameter optimization process. Specifically, the constraints apply to adjacent intervals at the connection point z. j A satisfiesj (z j ) = A j+1 (z j Furthermore, the first derivative is continuous, and its effect directly determines the accuracy and reliability of the area-storage capacity curve.

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

[0097]

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

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

[0100] S2123: Constructing the area-storage capacity curve;

[0101] Numerical integration based on piecewise polynomial parameters generates an area-storage capacity curve, represented as:

[0102]

[0103] The specific implementation uses the adaptive Simpson integral method, with the absolute error limit set to 10. -4 m 3 Integrate the elevation-area polynomial for each subinterval, sum the integral results of all subintervals, and obtain the total reservoir capacity V(z). The integral result is expressed in piecewise function form 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 fusion method for quantifying eco-hydrological interactions, characterized in that, Includes the following steps: Step 1: Collect environmental driving data and remote sensing data; Step 2: Preprocess the remote sensing data and combine it with environmentally driven data for dynamic monitoring of hydrological processes to obtain hydrological data; perform ecological parameter co-inversion on the preprocessed remote sensing data to obtain ecological indices; the hydrological data includes monthly lake area and water storage; the process of dynamic monitoring of hydrological processes is as follows: Step 211: Extract lake boundaries based on remote sensing data and calculate the monthly lake area; Step 212: Construct an area-storage 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; Step 3: Construct a spatiotemporal dataset using hydrological data, ecological indices, and environmentally driven data; Step 4: Construct a bidirectional LSTM neural network model and train and optimize it. Input the spatiotemporal dataset into the trained bidirectional LSTM neural network model to obtain the predicted water volume change, lag time, and corresponding ecological variables. Step 5: Construct an interaction quantification network diagram based on the predicted water volume change, lag time, and corresponding ecological variables.

2. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 1, characterized in that, The specific process of step 211 is as follows: Step 2111: Preprocess the remote sensing data; Step 2112: The pre-trained U-Net++ model and the super-resolution reconstruction method are used sequentially to extract the initial lake boundary raster data from the pre-processed remote sensing data; Step 2113: Post-process the initial lake boundary raster data to obtain accurate lake boundary raster data; post-processing optimization includes opening operation and spatiotemporal consistency correction; Step 2114: Perform raster statistics based on the accurate lake boundary raster data, calculate the lake area, and obtain the monthly lake area.

3. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 1, characterized in that, The specific process of step 212 is as follows: Step 2121: Perform weighted average fusion of the remote sensing data to obtain the fused DEM; Step 2122: Adaptive piecewise polynomial fitting is performed using fused DEM and precise lake boundary raster data to obtain the sub-interval elevation-area relationship and piecewise polynomial parameters; specifically including: Based on accurate lake boundary raster data, k sub-intervals were divided using Jenks natural faults. The elevation-area relationship of the sub-intervals was fitted by an adaptive piecewise polynomial. The polynomial coefficients of the elevation-area relationship of the sub-intervals were globally optimized using the Levenberg-Marquardt algorithm to obtain the piecewise polynomial parameters. Step 2123: Numerical integration of the piecewise polynomial parameters is performed using the adaptive Simpson integral method to construct the area-storage capacity curve.

4. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 1, characterized in that, The specific process of ecological parameter collaborative inversion is as follows: Step 221: Calculate the Normalized Difference Vegetation Index (NDVI) based on remote sensing data; Step 222: Calculate the total primary productivity (GPP) based on remote sensing data; Step 223: The ecological index is composed of the Normalized Difference Vegetation Index (NDVI) and the Total Primary Productivity (GPP).

5. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 4, characterized in that, The specific process of step 221 is as follows: Step 2211: Extract the near-infrared and red-band surface reflectance from the preprocessed remote sensing data; Step 2212: Calculate the Normalized Difference Vegetation Index (NDVI) based on the near-infrared and red-band surface reflectance.

6. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 4, characterized in that, The specific process of step 222 is as follows: Step 2221: Obtain the 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 based on the air temperature; Calculate water stress factors based on soil moisture content; Step 2223: Calculate the total primary productivity (GPP) using the light energy utilization efficiency model based on the photosynthetically active radiation absorption ratio, photosynthetically active radiation, temperature stress factor, and water stress factor.

7. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 1, characterized in that, Z-score normalization was used to standardize hydrological data, ecological indices, and environmental driving data to construct a spatiotemporal dataset.

8. The multi-scale fusion method for quantifying eco-hydrological interactions according to claim 1, characterized in that, The bidirectional LSTM neural network model consists of an input layer, a bidirectional LSTM layer, a fully connected layer, and an output layer connected sequentially. The bidirectional LSTM layer includes forward LSTM units, backward LSTM units, and a concatenation unit. The output layer solves for the water loss function and the ecological lag loss function, constructs a comprehensive loss function through weighted summation, calculates the comprehensive loss based on the comprehensive loss function, and optimizes the bidirectional LSTM neural network model through backpropagation. The bidirectional LSTM neural network model is trained using historical spatiotemporal datasets. The training process includes: Step 41: Calculate hydrological, ecological, and environmental features based on the spatiotemporal dataset to form the 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 potential future impact features, and the concatenation unit concatenates and merges the historical temporal dependencies and potential future impact features to obtain a comprehensive feature vector containing historical and future information. Step 44: The fully connected layer maps the integrated feature vector to the predicted water change value, lag time, and corresponding ecological variables; Step 45: The output layer calculates the water loss function based on the predicted water change value, and calculates the ecological lag loss function based on the lag time and the corresponding ecological variables. The water loss function and the ecological lag loss function are weighted and summed to construct the comprehensive loss function. The comprehensive loss is calculated, and the Bayesian optimization algorithm is used to perform a global search on the model hyperparameters to minimize the comprehensive loss and optimize the model.

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

Images