A method for evaluating lake ecological network resilience based on adaptive food web model

Abstract

The application discloses a kind of lake ecological network resilience evaluation methods based on adaptive food web model, method includes: obtaining the time series data of biological of target lake in recent years, obtain preprocessed time series data;Submerged plant area time series data in target lake is extracted;Based on the time series data after preprocessing, the adaptive food web model of area time series data composite disturbance of water level amplitude-human fishing is constructed;The model parameters of the adaptive food web model are estimated by Markov chain Monte Carlo method, and the optimal value of model parameter is obtained;Based on the adaptive food web model under the optimal value of model parameter, time-varying jacobian matrix is constructed, and the maximum lyapunov index under the interaction of different water level amplitude coefficients and human fishing intensity is calculated, to identify the water level amplitude threshold and human fishing intensity threshold leading to steady-state mutation, provide scientific basis for water level regulation and fishery management optimization.

Description

Technical Field

[0001] This invention belongs to the field of ecological protection technology and relates to a method for assessing the resilience of lake ecological networks based on an adaptive food web model. Background Technology

[0002] Shallow lakes, as an important component of Earth's ecosystem, play a vital role in maintaining biodiversity, purifying water quality, providing freshwater resources, and supporting fisheries. However, in recent years, shallow lakes worldwide have faced a dual challenge: firstly, intensified fishing pressure, typically manifested in overfishing leading to a sharp decline in fish biomass and subsequently triggering a series of extinctions; secondly, frequent abnormal water level fluctuations, mainly caused by seasonal or sudden rises and falls due to climate change and water conservancy project scheduling. These fluctuations typically alter the physicochemical properties of the water, directly impacting the photosynthesis, respiration, and nutrient cycling of aquatic plants. These two factors, as the two core drivers threatening the resilience of lake ecosystems, have attracted significant attention. Resilience refers to the ability of an ecosystem to recover its function or adapt to a new stable state through internal regulatory mechanisms after being disturbed. In other words, the system possesses self-repair capabilities, but current scientific understanding of it is largely at the qualitative descriptive level. Different systems face complex types of disturbances (natural factors such as climate fluctuations and biological invasions; anthropogenic factors such as pollution and overfishing), making the resilience of lake ecosystems complex and difficult to predict.

[0003] Traditional research often focuses on simple food web models, which are relatively simplistic and fail to incorporate biological behavioral adaptation strategies (such as dynamic adjustments in foraging and tolerance / adaptation to environmental stresses). This leads to an oversimplified characterization of ecosystem dynamics. While current research has made some progress, such as revealing the unidirectional effects of water level fluctuations on submerged plant distribution and fishing pressure on fish biomass, or simulating ecological responses to single disturbances using static food web models, significant shortcomings remain. These include a lack of quantification of synergistic mechanisms, static model simplification, reliance on data, and limitations in validation. This makes governance prone to idealized extrapolations, failing to adapt to complex real-world scenarios, and ultimately deviating from expected results.

[0004] Therefore, it is necessary to propose an assessment method for system resilience under the background of complex disturbances, which uses dynamic feedback modeling and collaborative mechanism quantification to scientifically reveal the resilience response law of the ecosystem under the complex disturbances of water level fluctuations and human fishing, and provide a basis for the optimization of water level regulation and fisheries management. Summary of the Invention

[0005] This invention aims to develop a method for assessing the resilience of lake ecological networks based on an adaptive food web model. This method, combining remote sensing imagery and field survey data, can accurately predict the response of lake resilience to water level fluctuations and human fishing, providing a quantitative basis for lake ecological management.

[0006] The technical solution adopted in this invention is as follows: A method for assessing the resilience of lake ecological networks based on an adaptive food web model, comprising the following steps: Step S1: Obtain time-series data on phytoplankton, zooplankton, benthic animals, and anchovy catch in the target lake over the years, and perform preprocessing to obtain preprocessed time-series data; Step S2: Acquire multispectral remote sensing images of the target lake, and extract the temporal data of submerged plants in the target lake using multispectral remote sensing image interpretation technology; Step S3: Based on the preprocessed time-series data and the area time-series data, construct an adaptive food web model subject to the combined disturbance of water level amplitude and human fishing; the adaptive food web model includes the water level amplitude coefficient and the intensity of human fishing. The adaptive food web model consists of differential equations with seven state variables, each corresponding to a dynamic model, which together constitute a dynamic system describing the dynamics of the lake ecological network. The seven state variables are: phytoplankton, submerged plants, zooplankton, benthic animals, anchovy catch, optimal water level for submerged plants, and anchovy feeding preference for benthic animals. The dynamic models are: phytoplankton population dynamic model, submerged plant population dynamic model, zooplankton population dynamic model, benthic animal population dynamic model, lake anchovy catch dynamic model, submerged plant optimal water level dynamic model, and lake anchovy adaptive foraging dynamic model under balanced catch effort conditions. Step S4: Use the Markov chain Monte Carlo method to estimate the model parameters of the adaptive food web model and obtain the optimal values ​​of the model parameters; Step S5: Construct a time-varying Jacobian matrix based on the adaptive food web model with optimal model parameters, and calculate the maximum Lyapunov exponent under the interaction of different water level amplitude coefficients and human fishing intensity, so as to identify the water level amplitude threshold and human fishing intensity threshold that lead to steady-state abrupt change.

[0007] Further, in step S2, the temporal data of submerged plants in the target lake are extracted using multispectral remote sensing image interpretation technology; the specific steps are as follows: Step S21: Acquire multispectral remote sensing images of the target lake during a period of less algal bloom, and perform radiometric calibration and geometric correction preprocessing to obtain preprocessed multispectral remote sensing images. Step S22: Calculate the aquatic vegetation index of each pixel in the preprocessed multispectral remote sensing image, set a discrimination threshold, and determine the pixels with the aquatic vegetation index greater than the discrimination threshold as aquatic vegetation areas. Step S23: Within the aquatic vegetation area, calculate the normalized vegetation index of each pixel, set a classification threshold, and determine the pixels with the normalized vegetation index less than the classification threshold as submerged plant areas. Step S24: Based on the spatial resolution of the preprocessed multispectral remote sensing image, sum up the areas of all pixels in the submerged plant area to obtain the temporal data of the submerged plant area in the target lake.

[0008] Furthermore, in step S3, the specific differential equations of the adaptive food web model specifically include: The seven state variables are: phytoplankton, submerged plants, zooplankton, benthic animals, anchovy catch, optimal water level for submerged plants, and anchovy feeding preference for benthic animals. The dynamic changes of these seven state variables are described by the aforementioned phytoplankton population dynamic model, submerged plant population dynamic model, zooplankton population dynamic model, benthic animal population dynamic model, anchovy catch dynamic model, submerged plant optimal water level dynamic model, and an anchovy adaptive foraging dynamic model under the condition of balanced catch effort. The specific equations are as follows: (1) Phytoplankton population dynamics model The formula is as follows: ; In the formula, Indicates phytoplankton, Indicates the effective photosynthetic rate of phytoplankton. Indicates the natural mortality rate of phytoplankton. Z represents the predation rate of zooplankton on phytoplankton, where Z represents zooplankton. This indicates the processing time of zooplankton for their prey, where P represents benthic animals. , 1- - These represent the feeding preferences of benthic animals for phytoplankton, submerged plants, and zooplankton, respectively. , , These represent the predation rates of benthic animals on phytoplankton, submerged plants, and zooplankton, respectively. This indicates the processing time of benthic animals for their prey; (2) Dynamic model of submerged plant populations The formula is as follows: ; In the formula, Indicates submerged plants. Indicates the effective photosynthetic rate of submerged plants. Indicates the mortality rate of submerged plants. Indicates the suitability of submerged plants; The mortality rate of submerged plants is a function of the deviation between the target lake level and the optimal water level for submerged plants; the formula is as follows: ; In the formula, Indicates the baseline mortality rate of submerged plants. Represented by natural base As the base, The coefficient representing the sensitivity of submerged plant mortality to changes in the target lake's water level is L, where L represents the target lake's water level and x represents the optimal water level for submerged plants. This indicates the niche width by which submerged plants adapt to changes in water level; (3) Zooplankton population dynamics model The formula is as follows: ; In the formula, This represents the assimilation coefficient of zooplankton. Indicates the natural mortality rate of zooplankton. , These represent the predation rates of the lake anchovy on zooplankton and benthic animals, respectively. This indicates the time it takes for the lake anchovy to process its prey. Indicates the fitness of zooplankton. 1- These represent the feeding preferences of lake anchovies for benthic and zooplankton, respectively. (4) Benthic animal population dynamics model The formula is as follows: ; In the formula, This represents the assimilation coefficient of benthic animals. Indicates the natural mortality rate of benthic animals. Indicates the fitness of benthic animals; (5) Dynamic model of lake anchovy catch under the condition of fishing effort equilibrium It is derived from the following formula: ; ; In the formula, This represents a population dynamic model of the lake anchovy. This represents the assimilation coefficient of the lake anchovy. The natural mortality rate of the lake anchovy is represented by q, and the human overfishing coefficient is represented by q. E indicates the suitability of the lake anchovy; E represents the fishing effort. This represents a dynamic model of fishing effort. Indicates the price of lake anchovies. Represents the sensitivity coefficient. This indicates the highest price for lake anchovies. This represents the investment cost per unit of effort. This represents the sensitivity coefficient of lake anchovy prices to lake anchovy catch volume. Indicates the price elasticity index; The fishing effort E is in equilibrium, i.e., the fishing effort dynamic model. At that time, the fishing effort under equilibrium conditions can be obtained. According to the relationship between the catch Q of anchovies and the catch F of anchovies: Q = The lake anchovy under equilibrium conditions can be obtained as follows: , can make , Indicating the intensity of human fishing, the lake anchovy under equilibrium conditions Substitute it back into the zooplankton population dynamics model and benthic animal population dynamics model The final zooplankton population dynamics model can be obtained by organizing the data. and benthic animal population dynamics model ; Based on the above-mentioned equilibrium condition of fishing effort E, the dynamic model of lake anchovy catch under the equilibrium condition of fishing effort can be obtained, and the formula is as follows: ; In the formula, Indicates the suitability of the lake anchovy. Indicates the suitability of lake anchovy catch; (6) Dynamic model of optimal water level for submerged plants The formula is as follows: ; In the formula, This represents the adaptation coefficient of submerged plants to changes in the water level of the target lake; (7) Dynamic model of adaptive foraging of lake anchovies The formula is as follows: ; In the formula, This indicates the foraging adaptability of the lake anchovy. This indicates the fitness selection gradient for anchovies. This represents the suitability selection gradient for lake anchovy catch.

[0009] Furthermore, in the phytoplankton population dynamics model, the effective photosynthetic rate of phytoplankton... The formula is as follows: ; In the formula, Indicates the maximum growth rate of phytoplankton. Indicates the turbidity of water. This represents the light attenuation coefficient of phytoplankton. Indicates the light compensation coefficient for phytoplankton. Indicates the intensity of incident light. This indicates the remaining light intensity reaching below the phytoplankton layer; Among them, the remaining light intensity reaching below the phytoplankton layer The formula is as follows: .

[0010] Furthermore, in the submerged plant population dynamics model, the effective photosynthetic rate of submerged plants... The formula is as follows: ; In the formula, This represents the maximum growth rate coefficient of submerged plants. Indicates the light attenuation coefficient of submerged plants. This represents the light compensation coefficient for submerged plants. This indicates the remaining light intensity reaching below the submerged plant layer; Among them, the remaining light intensity reaching below the submerged plant layer The formula is as follows: .

[0011] Furthermore, in the submerged plant population dynamics model, the target lake water level The formula is as follows: ; In the formula, Indicates the target lake in time water level, This represents the average of the lowest water levels of the target lake over the years. Indicates the water level amplitude coefficient. Represents the sine function. The time adjustment coefficient represents the periodic change.

[0012] Further, in step S4, the model parameters of the adaptive food web model are estimated using the Markov chain Monte Carlo method to obtain the optimal values ​​of the model parameters; specifically, this includes: Step S41: The preprocessed time series data and the area time series data of submerged plants in the target lake are used together as observations. The Markov chain Monte Carlo method is used to estimate the model parameters in the adaptive food web model to obtain the optimal values ​​of the model parameters. Step S42: Based on the optimal values ​​of the model parameters estimated by the Markov chain Monte Carlo method, the adaptive food web model is simulated to predict the long-term coexistence state of each state variable.

[0013] Further, in step S5, a time-varying Jacobian matrix is ​​constructed based on the adaptive food web model with optimal model parameters. The maximum Lyapunov exponent is calculated under the interaction of different water level amplitude coefficients and human fishing intensity, thereby identifying the water level amplitude threshold and human fishing intensity threshold that lead to steady-state abrupt changes; specifically including: Step S51: Based on the adaptive food web model with optimal model parameters, construct a time-varying Jacobian matrix J(t); wherein, the element in the i-th row and j-th column of the time-varying Jacobian matrix J(t) is the partial derivative of the differential equation of the i-th state variable with respect to the j-th state variable, and the values ​​of i and j range from 1 to 7. Step S52: Discretize the time evolution process of the dynamic system into N time steps, each time step having a length of... Substituting the value of the state variable at the k-th time step into the time-varying Jacobian matrix, the Jacobian matrix at the k-th time step is calculated. k=1,2,…,N; select an initial unit perturbation vector. Calculate the Lyapunov index The formula is as follows: ; In the formula, The Lyapunov index represents the Lyapunov index. Denotes the Euclidean vector norm; Step S53: By iterating through different water level amplitude coefficients and human fishing intensities, the dynamic system is run and the corresponding Lyapunov index is calculated. This Lyapunov index is the maximum Lyapunov index under the corresponding water level amplitude coefficient and human fishing intensity. A two-dimensional matrix diagram of the maximum Lyapunov index is constructed with the water level amplitude coefficient as the abscissa and the human fishing intensity as the ordinate. When the maximum Lyapunov index LLE is greater than 0, the dynamic system is determined to be in a low-toughness state. When the maximum Lyapunov index LLE is less than or equal to 0, the dynamic system is determined to be in a high-toughness state. In the two-dimensional matrix diagram, the parameter boundary value where the maximum Lyapunov index LLE changes from less than or equal to 0 to greater than 0 is determined as the water level amplitude threshold and the human fishing intensity threshold.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows: (1) Breaking through the limitations of traditional models that isolate and analyze single factors, an adaptive food web model driven by multiple factors is established. This invention constructs an adaptive food web model that integrates water level changes, human fishing activities, food webs (including plankton, submerged plants, benthic animals, and anchovies), and biological behavioral strategies. This model innovatively introduces two types of organisms' active adaptation capabilities to the environment: one is the adaptive dynamic adjustment of submerged plants to water level, and the other is the dynamic adjustment of anchovies' foraging preferences. In the dynamic process of lake organisms, the adaptive food web model comprehensively simulates the "behavioral adaptation-system recovery" response process under the combined disturbance of water level amplitude and human fishing.

[0015] (2) Breaking through the bottleneck of traditional methods that rely on direct data from Anchovy, this paper indirectly characterizes the biomass data of Anchovy through the historical catch volume and the equilibrium catch effort. Multi-source data fusion and Markov chain Monte Carlo (MCMC) methods are used for parameter estimation to improve model accuracy. The adaptive food web model based on optimal parameters enhances its ability to simulate and predict system dynamics under complex perturbations.

[0016] (3) The maximum Lyapunov index (LLE) is used to quantify the system resilience. By comparing the dynamic changes of the maximum Lyapunov index (LLE) of the system under the combined effects of multiple factors and the system with only a single factor, this invention reveals that the synergistic effect of the water level adaptability of submerged plants and the foraging adaptability of lake anchovies under the combined disturbance of water level amplitude and human fishing can enhance the system resilience. Based on this analysis, the adjustable water level range and the fishing intensity threshold are further determined, providing a scientific basis for the optimization of water level regulation and fishery management, and supporting the ecological restoration and ecological security management of lakes. Attached Figure Description

[0017] Figure 1 A flowchart of the method provided for this invention; Figure 2 A diagram of a lake food web structure provided by the present invention; Figure 3 The fitting effect diagram of parameter estimation using the Markov chain Monte Carlo method (MCMC) provided by the present invention; Figure 4 The optimal parameter distribution diagram provided by the present invention for parameter estimation using the Markov chain Monte Carlo method (MCMC); Figure 5 Simulation diagram of the optimal parameterized adaptive food web model provided for this invention; Figure 6 This invention provides a lake ecological network resilience assessment diagram under the combined disturbance of water level amplitude and human fishing. Detailed Implementation

[0018] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. The technical solutions adopted by the present invention are as follows: Figure 1 As shown.

[0019] Example As a typical shallow lake, Taihu Lake has faced the dual challenges of abnormal water level fluctuations and increased pressure from fishing in recent years. These two factors work together to affect the resilience of the ecosystem, making it an ideal area to verify the method of this invention.

[0020] The specific process is as follows: A method for assessing the resilience of lake ecological networks based on an adaptive food web model, comprising the following steps: Step S1: Obtain time-series data on phytoplankton, zooplankton, benthic animals, and anchovy catches in the shallow waters of Taihu Lake over the years, and perform preprocessing to obtain preprocessed time-series data, as follows: Step S11: Collect time-series data on phytoplankton, zooplankton, benthic animals, and anchovy catch in shallow lakes of Taihu Lake over the past 20 years; Step S12: Classify and organize the collected time series data, remove outliers, and complete the data preprocessing (such as taking the square root of the original data) to obtain the preprocessed time series data. Step S2: Acquire multispectral remote sensing images of shallow lakes in Taihu Lake, and extract the temporal data of submerged vegetation area in the shallow lakes of Taihu Lake using multispectral remote sensing image interpretation technology, as detailed below: Step S21: Acquire preprocessed surface reflectance images from Landsat 5, 7, and 8 on the Earth Engine (GEE) platform, extract cloud / cloud shadow markers using the Quality Assessment (QA) band and calculate cloud coverage, retain images with cloud content below 15%, and then extract the study area through regional cropping to obtain preprocessed multispectral remote sensing images. Step S22: Based on the preprocessed multispectral remote sensing image, calculate the aquatic vegetation index (AVI) value for each pixel within the shallow lake area of ​​Taihu Lake according to the formula for calculating the aquatic vegetation index (AVI). ; in, The AVI value represents the aquatic vegetation index. This represents the third component of the preprocessed multispectral remote sensing image after tassel transformation; Set the discrimination threshold a=-0.026, classify the aquatic vegetation index AVI value, and classify the pixels with an aquatic vegetation index AVI value greater than the discrimination threshold a as aquatic vegetation area, and the rest as non-aquatic vegetation areas, and extract all pixels of aquatic vegetation areas. Step S23: Calculate the Normalized Difference Vegetation Index (NDVI) value based on each pixel in the multispectral remote sensing image of the aquatic vegetation area. The calculation formula is as follows: ; in, This represents the Normalized Difference Vegetation Index (NDVI) value. Red band, It is in the near-infrared band; The classification threshold b for the Normalized Difference Vegetation Index (NDVI) is set to 0.2. Pixels with an NDVI value greater than the classification threshold b are classified as floating-leaved plants / emergent plants, and pixels with an NDVI value less than the classification threshold b are classified as submerged plants, thus completing the hierarchical classification and identification of aquatic vegetation types. Step S24: Based on the spatial resolution of the preprocessed multispectral remote sensing image, the area of ​​the pixels identified as submerged plants is statistically analyzed to generate temporal data of the area of ​​submerged plants in shallow lakes of Taihu Lake.

[0021] In step S3, based on the preprocessed time-series data and the area time-series data, an adaptive food web model is constructed under the combined disturbance of water level amplitude and human fishing; the adaptive food web model includes the water level amplitude coefficient and the intensity of human fishing. In step S3, the specific differential equations of the adaptive food web model include: The seven state variables are: phytoplankton, submerged plants, zooplankton, benthic animals, anchovy catch, optimal water level for submerged plants, and anchovy feeding preference for benthic animals. The dynamic changes of these seven state variables are described by the aforementioned phytoplankton population dynamic model, submerged plant population dynamic model, zooplankton population dynamic model, benthic animal population dynamic model, anchovy catch dynamic model, submerged plant optimal water level dynamic model, and an anchovy adaptive foraging dynamic model under the condition of balanced catch effort. The specific equations are as follows: (1) Phytoplankton population dynamics model The formula is as follows: ; In the formula, Indicates phytoplankton, Indicates the effective photosynthetic rate of phytoplankton. Indicates the natural mortality rate of phytoplankton. Z represents the predation rate of zooplankton on phytoplankton, where Z represents zooplankton. This indicates the processing time of zooplankton for their prey, where P represents benthic animals. , 1- - These represent the feeding preferences of benthic animals for phytoplankton, submerged plants, and zooplankton, respectively. , , These represent the predation rates of benthic animals on phytoplankton, submerged plants, and zooplankton, respectively. This indicates the processing time of benthic animals for their prey; (2) Dynamic model of submerged plant populations The formula is as follows: ; In the formula, Indicates submerged plants. Indicates the effective photosynthetic rate of submerged plants. Indicates the mortality rate of submerged plants. Indicates the suitability of submerged plants; The mortality rate of submerged plants is a function of the deviation between the target lake level and the optimal water level for submerged plants; the formula is as follows: ; In the formula, Indicates the baseline mortality rate of submerged plants. Represented by natural base As the base, The coefficient representing the sensitivity of submerged plant mortality to changes in the target lake's water level is L, where L represents the target lake's water level and x represents the optimal water level for submerged plants. This indicates the niche width by which submerged plants adapt to changes in water level; (3) Zooplankton population dynamics model The formula is as follows: ; In the formula, This represents the assimilation coefficient of zooplankton. Indicates the natural mortality rate of zooplankton. , These represent the predation rates of the lake anchovy on zooplankton and benthic animals, respectively. This indicates the time it takes for the lake anchovy to process its prey. Indicates the fitness of zooplankton. 1- These represent the feeding preferences of lake anchovies for benthic and zooplankton, respectively. (4) Benthic animal population dynamics model The formula is as follows: ; In the formula, This represents the assimilation coefficient of benthic animals. Indicates the natural mortality rate of benthic animals. Indicates the fitness of benthic animals; (5) Dynamic model of lake anchovy catch under the condition of fishing effort equilibrium It is derived from the following formula: ; ; In the formula, This represents a population dynamic model of the lake anchovy. This represents the assimilation coefficient of the lake anchovy. The natural mortality rate of the lake anchovy is represented by q, and the human overfishing coefficient is represented by q. E indicates the suitability of the lake anchovy; E represents the fishing effort. This represents a dynamic model of fishing effort. Indicates the price of lake anchovies. Represents the sensitivity coefficient. This indicates the highest price for lake anchovies. This represents the investment cost per unit of effort. This represents the sensitivity coefficient of lake anchovy prices to lake anchovy catch volume. Indicates the price elasticity index; The fishing effort E is in equilibrium, i.e., the fishing effort dynamic model. At that time, the fishing effort under equilibrium conditions can be obtained. According to the relationship between the catch Q of anchovies and the catch F of anchovies: Q = The lake anchovy under equilibrium conditions can be obtained as follows: , can make , Indicating the intensity of human fishing, the lake anchovy under equilibrium conditions Substitute it back into the zooplankton population dynamics model and benthic animal population dynamics model The final zooplankton population dynamics model can be obtained by organizing the data. and benthic animal population dynamics model ; Based on the above-mentioned equilibrium condition of fishing effort E, the dynamic model of lake anchovy catch under the equilibrium condition of fishing effort can be obtained, and the formula is as follows: ; In the formula, Indicates the suitability of the lake anchovy. Indicates the suitability of lake anchovy catch; (6) Dynamic model of optimal water level for submerged plants The formula is as follows: ; In the formula, This represents the adaptation coefficient of submerged plants to changes in the water level of the target lake; (7) Dynamic model of adaptive foraging of lake anchovies The formula is as follows: ; In the formula, This indicates the foraging adaptability of the lake anchovy. This indicates the fitness selection gradient for anchovies. This represents the suitability selection gradient for lake anchovy catch.

[0022] In the phytoplankton population dynamics model, the effective photosynthetic rate of phytoplankton The formula is as follows: ; In the formula, Indicates the maximum growth rate of phytoplankton. Indicates the turbidity of water. This represents the light attenuation coefficient of phytoplankton. Indicates the light compensation coefficient for phytoplankton. Indicates the intensity of incident light. This indicates the remaining light intensity reaching below the phytoplankton layer; Among them, the remaining light intensity reaching below the phytoplankton layer The formula is as follows: .

[0023] In the submerged plant population dynamics model, the effective photosynthetic rate of submerged plants The formula is as follows: ; In the formula, This represents the maximum growth rate coefficient of submerged plants. Indicates the light attenuation coefficient of submerged plants. This represents the light compensation coefficient for submerged plants. This indicates the remaining light intensity reaching below the submerged plant layer; Among them, the remaining light intensity reaching below the submerged plant layer The formula is as follows: .

[0024] In the submerged plant population dynamics model, the target lake water level The formula is as follows: ; In the formula, Indicates the target lake in time water level, This represents the average of the lowest water levels of the target lake over the years. Indicates the water level amplitude coefficient. Represents the sine function. The time adjustment coefficient represents the periodic change.

[0025] Step S4 involves estimating the model parameters of the adaptive food web model using the Markov chain Monte Carlo method to obtain the optimal values ​​of the model parameters; specifically, this includes: Step S41: Combining historical data on submerged plant area, phytoplankton concentration, zooplankton concentration, benthic abundance, and anchovy catch in shallow Taihu Lake from 2000 to 2020, the model parameters of the adaptive food web model are estimated using the Markov Chain Monte Carlo (MCMC) method. This method employs a delayed rejection adaptive Monte Carlo (DRAM) algorithm within a Bayesian framework, using modeling software for parameter sampling. The error between the model prediction and the observed data follows an independent and identically distributed normal distribution, thus constructing the likelihood function using the sum of squared errors. The delayed rejection adaptive Monte Carlo (DRAM) algorithm uses a scaling factor of 0.001 on the identity matrix as the initial value for the proposed covariance, adaptively adjusting it every 500 steps, running a total of 10,000 iterations. The posterior mean, standard deviation, Monte Carlo standard error (MC_err), autocorrelation time (tau), and convergence diagnostic index (Geweke) are calculated using Markov chain statistical diagnostics (chainstats). The system dynamics are reconstructed by using the posterior mean of the parameters to drive the differential equation model, and the observation data is superimposed to visualize the fitting effect. Finally, the predicted distribution is generated to evaluate the model's fitting ability.

[0026] The results of the model parameter estimation analysis are shown in Table 1. The core statistics show that the Monte Carlo error (MC_err) is much lower than the standard deviation (std), the sampling accuracy meets the requirements, and the convergence diagnostic index (Geweke) is around 95%, indicating that the Markov chain has converged.

[0027] Results of Markov Chain Monte Carlo Method (MCMC): Fitting results of parameter estimation using Markov Chain Monte Carlo Method (MCMC) (e.g.) Figure 3 ) and the distribution of optimal parameter values ​​using the Markov chain Monte Carlo method (MCMC) for parameter estimation (e.g. Figure 4 ),like Figure 3 As shown, Figure 3 It contains five sub-graphs: top, top-middle, middle, bottom-middle, and bottom. Figure 3The horizontal axis of the five subplots represents the year, and the vertical axis represents the biomass of each species. From top to bottom, they show the observed values ​​of five species—phytoplankton, submerged plants, zooplankton, benthic animals, and anchovies—between 2000 and 2020, along with the dynamic model fitting curves for each species. The dynamic model effectively captures the dynamic trends of each variable, and combined with... Figure 4 The results show that the posterior distribution of the 44 parameters is approximately normal. Combined with the analysis of the above parameter estimates, this indicates that the estimates of the 44 parameters are stable and reliable.

[0028] Table 1. Parameter estimation results of the adaptive food web model based on the Markov chain Monte Carlo method.

[0029] Step S42: Based on the optimal values ​​of the model parameters estimated using the Markov chain Monte Carlo method, the adaptive food web model is simulated to predict the long-term coexistence state of each state variable. The results are as follows: Figure 5 As shown: Figure 5 It contains four sub-images: top left, top right, bottom left, and bottom right. Figure 5 The graphs show the changes in submerged plant biomass over time under four combinations of two different adaptations. The horizontal axis of each subgraph represents the simulated time, and the vertical axis represents the simulated submerged plant biomass. In the top left subgraph (A: g=0.1, V=0.1), submerged plants exhibit adaptation to water level, and anchovies exhibit adaptive foraging. Under these conditions, phytoplankton, submerged plants, zooplankton, benthic animals, and anchovies can coexist for a long period, with a richer biomass. This difference stems from the ecological benefit effect of the synergistic regulation of two factors: the anchovy's foraging preferences adjust to alleviate the pressure of a single prey, while the submerged plant's water level adaptation improves vegetation photosynthesis and rhizosphere carbon fixation capacity to enhance habitat heterogeneity; in the upper right subgraph (B: g=0.1, V=0), the anchovy possesses foraging adaptation, and various species can coexist for a long time; in the lower left subgraph (C: g=0, V=0.1), even with the introduction of submerged plant water level adaptation, the anchovy's foraging adaptation cannot respond to vegetation changes, still causing food chain disruption and water quality deterioration, and all species failed to coexist for a long time; in the lower right subgraph (D: g=0, When V=0, submerged plants are the first to become extinct because they lack the adaptability to water level and the adaptability of anchovies to forage. This is followed by benthic animals (due to lack of food and habitat), phytoplankton (due to water quality deterioration caused by the reduction of submerged plants), and zooplankton (due to the chain reaction of food chains). Ultimately, anchovies approach extinction due to food scarcity and environmental degradation, and the ecosystem collapses.

[0030] Step S5: Construct a time-varying Jacobian matrix based on the adaptive food web model with optimal model parameters, and calculate the maximum Lyapunov exponent under the interaction (compound perturbation) of different water level amplitude coefficients and human fishing intensity, thereby identifying the water level amplitude threshold and human fishing intensity threshold that lead to steady-state abrupt changes; specifically including: Step S51: Based on the adaptive food web model with optimal model parameters, construct a time-varying Jacobian matrix J(t); wherein, the element in the i-th row and j-th column of the time-varying Jacobian matrix J(t) is the partial derivative of the differential equation of the i-th state variable with respect to the j-th state variable, and the values ​​of i and j range from 1 to 7. Step S52: Discretize the time evolution process of the dynamic system into N time steps, each time step having a length of... Substituting the value of the state variable at the k-th time step into the time-varying Jacobian matrix, the Jacobian matrix at the k-th time step is calculated. k=1,2,…,N; select an initial unit perturbation vector. Calculate the Lyapunov index The formula is as follows: ; In the formula, The Lyapunov index represents the Lyapunov index. Denotes the Euclidean vector norm; Step S53: By iterating through different water level amplitude coefficients and human fishing intensities, the dynamic system is run and the corresponding Lyapunov index is calculated. This Lyapunov index is the maximum Lyapunov index under the corresponding water level amplitude coefficient and human fishing intensity. A two-dimensional matrix diagram of the maximum Lyapunov index is constructed with the water level amplitude coefficient as the abscissa and the human fishing intensity as the ordinate. When the maximum Lyapunov index LLE is greater than 0, the dynamic system is determined to be in a low-toughness state. When the maximum Lyapunov index LLE is less than or equal to 0, the dynamic system is determined to be in a high-toughness state. In the two-dimensional matrix diagram, the parameter boundary value where the maximum Lyapunov index LLE changes from less than or equal to 0 to greater than 0 is determined as the water level amplitude threshold and the human fishing intensity threshold.

[0031] The lake ecological network resilience assessment map provided by this invention under the combined disturbance of water level amplitude and human fishing is as follows: Figure 6 As shown: Based on the maximum Lyapunov index (LLE) Figure 6 Classified into three levels of resilience zones: Zone I with higher resilience ( ), moderate toughness zone II (0 < LLE < 0.15), and low toughness zone III (LLE ≥ 0.15). Figure 6The diagram includes left and right subplots. In the left subplot (A: g=0.1, V=0.1), submerged plants adapt to water levels and anchovies adapt their foraging, with highly resilient areas accounting for 6.13% of the total area and less resilient areas accounting for 30.99%. In the right subplot (B: g=0.1, V=0), anchovies adapt their foraging, with highly resilient areas accounting for only 3.93% and less resilient areas accounting for 38.85%. Under the combined disturbance of water level amplitude and human fishing, the introduction of both submerged plant adaptability to water levels and anchovy foraging adaptability further enhances the resilience of the dynamic system to high-intensity external disturbances.

[0032] Further analysis of the threshold response: Left subplot (A: g=0.1, V=0.1). Under these conditions, the dynamic system simultaneously exhibits adaptation to submerged plant water levels and foraging adaptation to anchovies. When the maximum Lyapunov exponent (LLE) of the dynamic system in the left subplot is ≤0, the water level amplitude coefficient (γ) ranges from 0 to γ ​​to 0.35, indicating an increase in human fishing intensity ( The range is 0 < <0.665; Right subplot (B: g=0.1, V=0), under these conditions, the dynamic system only exhibits the foraging adaptability of anchovies. Correspondingly, when the maximum Lyapunov exponent LLE ≤ 0 in the dynamic system of the right subplot, the range of water level amplitude coefficient narrows to 0 < γ < 0.2, and the range of human fishing intensity narrows to 0 < <0.655. When either the water level amplitude coefficient or the intensity of human fishing exceeds its corresponding threshold, the resilience of the dynamic system will be significantly reduced.

[0033] Based on the lake ecological network resilience assessment map under the combined disturbance of water level amplitude and human fishing, the resilience changes of the lake under the interaction of water level amplitude coefficient and human fishing intensity can be clarified, providing a reference for water level regulation and fisheries management optimization.

Claims

1. A method for assessing the resilience of lake ecological networks based on an adaptive food web model, characterized in that, Includes the following steps: Step S1: Obtain time-series data on phytoplankton, zooplankton, benthic animals, and anchovy catch in the target lake over the years, and preprocess the data to obtain preprocessed time-series data. Step S2: Acquire multispectral remote sensing images of the target lake, and extract the temporal data of submerged plants in the target lake using multispectral remote sensing image interpretation technology; Step S3: Based on the preprocessed time-series data and the area time-series data, construct an adaptive food web model subject to the combined disturbance of water level amplitude and human fishing; the adaptive food web model includes the water level amplitude coefficient and the intensity of human fishing. The adaptive food web model consists of differential equations with seven state variables, each corresponding to a dynamic model, which together constitute a dynamic system describing the dynamics of the lake ecological network. The seven state variables are: phytoplankton, submerged plants, zooplankton, benthic animals, anchovy catch, optimal water level for submerged plants, and anchovy feeding preference for benthic animals. The dynamic models are: phytoplankton population dynamic model, submerged plant population dynamic model, zooplankton population dynamic model, benthic animal population dynamic model, lake anchovy catch dynamic model, submerged plant optimal water level dynamic model, and lake anchovy adaptive foraging dynamic model under balanced catch effort conditions. Among them, the submerged plant population dynamic model The formula is as follows: ； In the formula, Indicates submerged plants. Indicates the effective photosynthetic rate of submerged plants. Indicates the mortality rate of submerged plants. Indicates the suitability of submerged plants; , 1- - These represent the feeding preferences of benthic animals for phytoplankton, submerged plants, and zooplankton, respectively. , , These represent the predation rates of benthic animals on phytoplankton, submerged plants, and zooplankton, respectively. Indicates benthic animals, This indicates the processing time of benthic animals for their prey. Indicates phytoplankton, Indicates zooplankton; The mortality rate of submerged plants is a function of the deviation between the target lake level and the optimal water level for submerged plants; the formula is as follows: ； In the formula, Indicates the baseline mortality rate of submerged plants. Represented by natural base As the base, The coefficient representing the sensitivity of submerged plant mortality to changes in the target lake's water level is L, where L represents the target lake's water level and x represents the optimal water level for submerged plants. This indicates the niche width by which submerged plants adapt to changes in water level; target lake water level The formula is as follows: ； In the formula, Indicates the target lake in time water level, This represents the average of the lowest water levels of the target lake over the years. Indicates the water level amplitude coefficient. Represents the sine function. A time adjustment coefficient representing periodic changes; The optimal water level dynamic model for submerged plants The formula is as follows: ； In the formula, This represents the adaptation coefficient of submerged plants to changes in the water level of the target lake; The adaptive foraging dynamic model of the lake anchovy The formula is as follows: ； In the formula, This indicates the foraging adaptability of the lake anchovy. 1- These represent the feeding preferences of the lake anchovy towards benthic and zooplankton, respectively. This represents the assimilation coefficient of the lake anchovy. , These represent the predation rates of the lake anchovy on zooplankton and benthic animals, respectively. This indicates the time it takes for the lake anchovy to process its prey; Step S4: Use the Markov chain Monte Carlo method to estimate the model parameters of the adaptive food web model and obtain the optimal values ​​of the model parameters; Step S5: Construct a time-varying Jacobian matrix based on the adaptive food web model with optimal model parameters, and calculate the maximum Lyapunov exponent under the interaction of different water level amplitude coefficients and human fishing intensity, so as to identify the water level amplitude threshold and human fishing intensity threshold that lead to steady-state abrupt change.

2. The method for assessing the resilience of lake ecological networks based on an adaptive food web model according to claim 1, characterized in that, In step S2, the temporal data of submerged plants in the target lake are extracted using multispectral remote sensing image interpretation technology; the specific steps are as follows: Step S21: Acquire multispectral remote sensing images of the target lake and perform radiometric calibration and geometric correction preprocessing to obtain preprocessed multispectral remote sensing images; Step S22: Calculate the aquatic vegetation index of each pixel in the preprocessed multispectral remote sensing image, set a discrimination threshold, and determine the pixels with the aquatic vegetation index greater than the discrimination threshold as aquatic vegetation areas. Step S23: Within the aquatic vegetation area, calculate the normalized vegetation index of each pixel, set a classification threshold, and determine the pixels with the normalized vegetation index less than the classification threshold as submerged plant areas. Step S24: Based on the spatial resolution of the preprocessed multispectral remote sensing image, sum up the areas of all pixels in the submerged plant area to obtain the temporal data of the submerged plant area in the target lake.

3. The method for assessing the resilience of lake ecological networks based on an adaptive food web model according to claim 2, characterized in that, In step S3, the specific differential equations of the adaptive food web model include: The phytoplankton population dynamics model The formula is as follows: ； In the formula, Indicates the effective photosynthetic rate of phytoplankton. Indicates the natural mortality rate of phytoplankton. This indicates the predation rate of zooplankton on phytoplankton. Indicates the time it takes for zooplankton to process their food; The zooplankton population dynamics model The formula is as follows: ； In the formula, This represents the assimilation coefficient of zooplankton. Indicates the natural mortality rate of zooplankton. Indicates the fitness of zooplankton; The benthic animal population dynamics model The formula is as follows: ； In the formula, This represents the assimilation coefficient of benthic animals. Indicates the natural mortality rate of benthic animals. Indicates the fitness of benthic animals; Dynamic model of anchovy catch under equilibrium fishing effort It is derived from the following formula: ； ； In the formula, This represents a population dynamic model of the lake anchovy. This represents the assimilation coefficient of the lake anchovy. The natural mortality rate of the lake anchovy is represented by q, and the human overfishing coefficient is represented by q. E indicates the suitability of the lake anchovy; E represents the fishing effort. This represents a dynamic model of fishing effort. Indicates the price of lake anchovies. Represents the sensitivity coefficient. This indicates the highest price for lake anchovies. This represents the investment cost per unit of effort. This represents the sensitivity coefficient of lake anchovy prices to lake anchovy catch volume. Indicates the price elasticity index; The fishing effort E is in equilibrium, i.e., the fishing effort dynamic model. At that time, the fishing effort under equilibrium conditions is obtained. According to the relationship between the catch Q of anchovies and the catch F of anchovies: Q = The lake anchovy under equilibrium conditions is obtained as follows: ,make , Indicating the intensity of human fishing, the lake anchovy under equilibrium conditions Substitute it back into the zooplankton population dynamics model and benthic animal population dynamics model The final zooplankton population dynamics model was obtained. and benthic animal population dynamics model ; Based on the above-mentioned equilibrium condition of fishing effort E, the dynamic model of lake anchovy catch under the equilibrium condition of fishing effort is obtained, and the formula is as follows: ； In the formula, This indicates the suitability of the lake anchovy catch.

4. The method for assessing the resilience of lake ecological networks based on an adaptive food web model according to claim 3, characterized in that, In the phytoplankton population dynamics model, the effective photosynthetic rate of phytoplankton The formula is as follows: ； In the formula, Indicates the maximum growth rate of phytoplankton. Indicates the turbidity of water. This represents the light attenuation coefficient of phytoplankton. Indicates the light compensation coefficient for phytoplankton. Indicates the intensity of incident light. This indicates the remaining light intensity reaching below the phytoplankton layer; Among them, the remaining light intensity reaching below the phytoplankton layer The formula is as follows: 。 5. The method for assessing the resilience of lake ecological networks based on an adaptive food web model according to claim 3, characterized in that, In the submerged plant population dynamics model, the effective photosynthetic rate of submerged plants The formula is as follows: ； In the formula, This represents the maximum growth rate coefficient of submerged plants. Indicates the light attenuation coefficient of submerged plants. This represents the light compensation coefficient for submerged plants. This indicates the remaining light intensity reaching below the submerged plant layer; Among them, the remaining light intensity reaching below the submerged plant layer The formula is as follows: 。 6. The method for assessing the resilience of lake ecological networks based on an adaptive food web model according to claim 1, characterized in that, Step S4 involves estimating the model parameters of the adaptive food web model using the Markov chain Monte Carlo method to obtain the optimal values ​​of the model parameters; specifically, this includes: Step S41: The preprocessed time series data and the area time series data of submerged plants in the target lake are used together as observations. The Markov chain Monte Carlo method is used to estimate the model parameters in the adaptive food web model to obtain the optimal values ​​of the model parameters. Step S42: Based on the optimal values ​​of the model parameters estimated by the Markov chain Monte Carlo method, the adaptive food web model is simulated to predict the long-term coexistence state of each state variable.

7. The method for assessing the resilience of lake ecological networks based on an adaptive food web model according to claim 1, characterized in that, Step S5: Construct a time-varying Jacobian matrix based on the adaptive food web model with optimal model parameters, and calculate the maximum Lyapunov exponent under the interaction of different water level amplitude coefficients and human fishing intensity, thereby identifying the water level amplitude threshold and human fishing intensity threshold that lead to steady-state abrupt changes; specifically including: Step S51: Based on the adaptive food web model with optimal model parameters, construct a time-varying Jacobian matrix J(t); wherein, the element in the i-th row and j-th column of the time-varying Jacobian matrix J(t) is the partial derivative of the differential equation of the i-th state variable with respect to the j-th state variable, and the values ​​of i and j range from 1 to 7. Step S52: Discretize the time evolution process of the dynamic system into N time steps, each time step having a length of... Substituting the value of the state variable at the k-th time step into the time-varying Jacobian matrix, the Jacobian matrix at the k-th time step is calculated. k=1,2,…,N; select an initial unit perturbation vector. Calculate the Lyapunov index The formula is as follows: ； In the formula, The Lyapunov index represents the index of Lyapunov. Denotes the Euclidean vector norm; Step S53: By iterating through different water level amplitude coefficients and human fishing intensities, the dynamic system is run and the corresponding Lyapunov index is calculated. This Lyapunov index is the maximum Lyapunov index under the corresponding water level amplitude coefficient and human fishing intensity. A two-dimensional matrix diagram of the maximum Lyapunov index is constructed with the water level amplitude coefficient as the abscissa and the human fishing intensity as the ordinate. When the maximum Lyapunov index LLE is greater than 0, the dynamic system is determined to be in a low-toughness state. When the maximum Lyapunov index LLE is less than or equal to 0, the dynamic system is determined to be in a high-toughness state. In the two-dimensional matrix diagram, the parameter boundary value where the maximum Lyapunov index LLE changes from less than or equal to 0 to greater than 0 is determined as the water level amplitude threshold and the human fishing intensity threshold.

Images