Method for identifying key driving factors of evolution of ecosystem structure and function in estuary and adjacent sea area

Abstract

The application provides a method for identifying key driving factors of estuary and adjacent sea ecosystem structure and function evolution, and belongs to the technical field of ecological systems. The application preprocesses multi-source long time series data and constructs an ecological system evolution index system, realizes interpretable original physical variable space contribution degree decomposition under strong collinearity conditions by using a two-stage collinearity decomposition framework, extracts time-varying time lag by using wavelet coherence analysis and adaptive dynamic time warping, models multi-state time lag transition probability by using a hidden Markov model to form a time-varying causal structure, inputs a reservoir pool causal reinforcement artificial intelligence model for comprehensive inference, and outputs contribution degree scores, stability scores and confidence scores of each candidate driving factor, thereby solving the technical problem that key driving factors of an estuary ecological system cannot be stably and quantitatively attributed under the coexistence conditions of strong collinearity, nonlinear driving response relationship and high-dimensional sparse observation.

Description

Technical Field

[0001] This invention belongs to the field of ecosystem technology, and more specifically, relates to a method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems. Background Technology

[0002] Estuarine and adjacent marine ecosystems are regulated by a variety of physical, chemical, and biological processes, making the identification of driving factors a core scientific task in ecosystem management. Currently, researchers typically employ multiple linear regression, generalized additive models, or traditional attribution methods based on variance decomposition to quantify the associations between candidate driving factors such as tides, runoff, and nutrients and ecological response variables. These methods are applicable to scenarios where variables are independent and the response relationship is approximately linear, and have been widely used in attribution studies of nearshore eutrophication, dissolved oxygen anomalies, and algal blooms.

[0003] However, strong collinearity often exists among candidate driving factors in estuarine environments. Under these conditions, traditional regularization methods systematically compress regression coefficients, leading to the underestimation or confounding of the independent contributions of each factor. Simultaneously, the impact of driving factors on ecosystems generally exhibits seasonally dependent, time-varying, and time-lag characteristics, which traditional static correlation analysis cannot capture. Furthermore, estuarine ecological observation data are often sparsely sampled with high missing rates, making traditional methods prone to overfitting in high-dimensional sparse feature spaces and failing to guarantee generalization ability during cross-regional migration. In current estuarine ecological attribution research, due to the simultaneous presence of the triple challenges of strong collinearity, time-varying causal structures, and high-dimensional sparse observations, existing methods cannot achieve interpretable independent contribution decomposition under collinearity conditions, nor can they robustly infer time-varying driving mechanisms under sparse data conditions. In other words, existing technologies suffer from the technical problem of being unable to achieve robust and interpretable quantitative attribution of key driving factors in estuarine ecosystems under the coexistence of strong collinearity, nonlinear driving-response relationships, and high-dimensional sparse observations. Summary of the Invention

[0004] In view of this, the present invention provides a method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems, which can solve the technical problem in the prior art that it is impossible to make robust and interpretable quantitative attribution of key driving factors of estuarine ecosystems under the conditions of strong collinearity, nonlinear driving response relationship and high-dimensional sparse observation.

[0005] This invention is implemented as follows: This invention provides a method for identifying key driving factors in the structural and functional evolution of estuarine and adjacent marine ecosystems, comprising the following steps: We acquired multi-source long-term series data of the study area, preprocessed the multi-source long-term series data, and constructed an indicator system for the evolution of ecosystem structure and function based on the preprocessed multi-source long-term series data. Calculate the variance inflation factor among all candidate driving factors in the preprocessed multi-source long-term series data, compare the variance inflation factor with the variance inflation factor threshold, and when the variance inflation factor exceeds the variance inflation factor threshold, initiate a two-stage collinearity decomposition framework to obtain the spatial contribution of the original physical variables. Wavelet coherence analysis is used to extract the phase lag spectrum from the preprocessed multi-source long-time series data. The optimal time lag is estimated in different seasonal windows by combining adaptive dynamic time warping. Hidden Markov model is used to integrate the multi-state time lag transition probabilities to form a time-varying causal structure. Based on the knowledge graph of ecological mechanism, the first-order physical driving factors are forcibly retained. The selection frequency of each candidate driving factor is statistically selected by the stability selection of elastic network. The candidate driving factors with selection frequency exceeding the frequency threshold are retained as stable feature subsets. The stable feature subsets are evaluated by leave-one-out cross-validation to obtain the generalization ability evaluation results. The candidate driving factors in the stable feature subset are ranked by the driving factor importance ranking algorithm based on information geometry geodesic distance. The geodesic curvature integral value is calculated as the driving importance score. The ratio of the highest driving importance score to the second highest driving importance score is compared with the importance ratio threshold to obtain the dominant driving factor or composite driving factor group. The time-varying causal structure, driving importance score, and stable feature subset are input into the reserve pool causal enhancement model for comprehensive inference. The contribution score, stability score, and confidence score of each candidate driving factor are output to form a list of key driving factors and their period of action, scale of action, and driving type.

[0006] The multi-source long-term series data includes marine meteorological data, hydrodynamic data, water quality and environmental data, ecological observation data, land-based input data, and human activity intensity data. The preprocessing includes three steps: spatiotemporal alignment, quality control, and missing time step completion.

[0007] The missing time step completion refers to inputting multi-source long-term series data, which has undergone spatiotemporal alignment and quality control, into the causal enhancement model of the reserve pool to complete the missing time steps and obtain complete multi-source long-term series data.

[0008] Specifically, the two-stage collinearity decomposition framework first uses principal component analysis to orthogonalize candidate driving factors to eliminate linear collinearity, then calculates the Shapley additive explanatory value of each candidate driving factor in a low-dimensional orthogonal space using an ensemble tree model, and finally maps the Shapley additive explanatory value back to the original physical variable space through inverse transformation.

[0009] The phase lag spectrum is the output of wavelet coherence analysis, which represents the change matrix of the phase difference between candidate driving factors and ecological response variables at each frequency component over time, and serves as the input for adaptive dynamic time warping.

[0010] Specifically, the adaptive dynamic time warping involves using dynamic programming to search for the optimal alignment path between two time series that minimizes the cumulative distance within different seasonal windows, thereby independently estimating the optimal time delay within each seasonal window. The optimal time delay serves as the input of the observation sequence to the hidden Markov model.

[0011] Specifically, the selection of elastic network stability involves repeatedly running elastic network regularized regression on multiple bootstrap subsamples and calculating the selection frequency of each candidate driving factor into the model. The selection frequency is the ratio of the number of times each candidate driving factor is selected into the model to the total number of bootstraps.

[0012] The driving factor importance ranking algorithm based on information geometry geodesic distance uses the Fisher information matrix as the Riemann metric tensor, treats each candidate driving factor as a coordinate in a family of parameterized probability distributions, calculates the geodesic in the coordinate direction of each candidate driving factor using exponential mapping, numerically integrates the curvature of each geodesic, and uses the integral value of the geodesic curvature as the driving importance score.

[0013] The reservoir causal enhancement model consists of four sequentially connected parts: a missing perception input layer, an echo state network reservoir, a sparse causal transformer, and an output readout layer.

[0014] The missing information input layer receives a mixed input of high-frequency continuous physical-driven time series and monthly-scale sparse ecological observations. It replaces the missing time steps with learnable missing embedding vectors and applies a zero-weight mask to the missing positions during loss calculation.

[0015] The echo state network reservoir is composed of a large random sparse cyclic reservoir. Its internal weights are fixed after initialization and do not participate in training; only the linear weights of the output readout layer are trained. The attention matrix of the sparse causal transformer is forcibly subjected to a lower triangular causal mask, and additional causal masks are applied to the attention weights. Regularization is used to generate sparse attention maps.

[0016] Specifically, the training of the reservoir causal enhancement model is performed using a model-independent meta-learning variant framework on the source domain dataset for meta-learning pre-training. After the meta-learning pre-training is completed, the output readout layer and sparse causal transformer are fine-tuned using labeled samples from the target domain dataset, and the weights inside the echo state network reservoir are frozen.

[0017] The reservoir causal enhancement model calculates the memory pressure index based on the sequence length of the current batch, the number of nodes in the echo state network reservoir, and the number of layers in the sparse causal transformer. It then dynamically adjusts the number of CUDA streams, batch size, gradient checkpoints, and floating-point arithmetic precision based on the memory pressure index.

[0018] The Fisher information matrix is ​​calculated using a low-rank approximation method for low-rank decomposition, and the low-rank number is determined by Pareto front analysis of reconstruction error and computation time.

[0019] Wherein, the variance inflation factor threshold is 10, the frequency threshold is 80%, the importance ratio threshold is 1.5, the low-rank number ranges from 5 to 30, the number of nodes in the echo state network reserve pool ranges from 4000 to 6000, the connectivity ranges from 0.5% to 2%, and the number of sparse causal transformer layers ranges from 2 to 4. The range of regularization weights is to Time series with missing rates exceeding 40% were preprocessed using linear interpolation and mean padding. 10 to 20 labeled samples were retained from the target domain dataset for rapid adaptation.

[0020] This invention proposes a method for identifying key driving factors in estuarine ecosystems that integrates two-stage collinear decomposition, time-varying causal structure modeling, and reservoir causal enhancement model. It provides a systematic solution to the attribution problem of strong collinearity, nonlinear driving response relationships, and high-dimensional sparse observations.

[0021] This invention employs a two-stage collinearity decomposition framework. First, candidate driving factors are orthogonalized to eliminate linear collinearity. Then, the contribution values ​​in the low-dimensional space are mapped back to the original physical variable space via inverse transformation, thus avoiding the underestimation of contribution caused by coefficient compression in traditional regularization methods. Time-varying delays are extracted through wavelet coherence analysis and adaptive dynamic time warping, and multi-state delay transition probabilities are modeled using a Hidden Markov Model. This enables the method to capture seasonally dependent dynamic causal structures, rather than relying on static correlation assumptions. Through the meta-learning pre-training mechanism of the reservoir causal enhancement model, the model can quickly adapt with a small number of target domain samples. Combined with the specialized processing of sparse observations by the missing perception input layer, the risk of overfitting under high-dimensional sparse conditions is effectively suppressed.

[0022] In summary, this invention solves the technical problem mentioned in the background art of being unable to perform robust and interpretable quantitative attribution of key driving factors of estuarine ecosystems under the conditions of strong collinearity, nonlinear driving response relationships and high-dimensional sparse observations. Attached Figure Description

[0023] Figure 1 This is a flowchart of the method of the present invention.

[0024] Figure 2 This is a spatial distribution diagram of the original physical variables of the candidate driving factors.

[0025] Figure 3 This is a time-varying causal structure probability diagram of runoff and chlorophyll a concentration.

[0026] Figure 4 A ranking plot of importance scores driven by candidate driving factors.

[0027] Figure 5 This is a sparse attention graph for a sparse causal transformer.

[0028] Figure 6 A graph showing the dynamic changes in the seasonal contribution scores of each candidate driving factor. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0030] like Figure 1 The diagram shown is a flowchart of a method for identifying key driving factors in the structural and functional evolution of estuary and adjacent marine ecosystems, provided by this invention. This method includes the following steps: S01. Obtain multi-source long-term series data of the study area from monthly to annual scales, preprocess the multi-source long-term series data, and construct an indicator system for the evolution of ecosystem structure and function based on the preprocessed multi-source long-term series data. S02. Calculate the variance inflation factor among all candidate driving factors in the preprocessed multi-source long-term series data, compare the variance inflation factor with the variance inflation factor threshold, and when the variance inflation factor exceeds the variance inflation factor threshold, start the two-stage collinearity decomposition framework to obtain the spatial contribution of the original physical variables. S03. Wavelet coherence analysis is used to extract the phase lag spectrum from the preprocessed multi-source long-time series data. The optimal time lag is estimated in different seasonal windows by combining adaptive dynamic time warping. The multi-state time lag transition probability is integrated by a hidden Markov model to form a time-varying causal structure. S04. Based on the ecological mechanism knowledge graph, the first-order physical driving factors are forcibly retained. The selection frequency of each candidate driving factor is statistically analyzed on the bootstrap subsample using the elastic network stability selection. Candidate driving factors with selection frequencies exceeding the frequency threshold are retained as stable feature subsets. The stable feature subsets are evaluated using leave-one-out cross-validation to obtain the generalization ability evaluation results. S05. Sort each candidate driving factor in the stable feature subset by the driving factor importance ranking algorithm based on information geometric geodesic distance, calculate the geodesic curvature integral value as the driving importance score, and compare the ratio of the highest driving importance score to the second highest driving importance score with the importance ratio threshold to obtain the dominant driving factor or composite driving factor group. S06. The time-varying causal structure, the driving importance score, and the stable feature subset are input into the reserve pool causal enhancement model for comprehensive inference. The contribution score, stability score, and confidence score of each candidate driving factor are output to form a list of key driving factors and their period of action, scale of action, and driving type.

[0031] The multi-source long-term series data includes marine meteorological data, hydrodynamic data, water quality and environmental data, ecological observation data, terrestrial input data, and human activity intensity data. The preprocessing includes three stages: spatiotemporal alignment, quality control, and missing time step completion. Spatiotemporal alignment refers to interpolating data from different sources with different spatial and temporal resolutions to the same spatial grid and time step. Missing time step completion refers to inputting the spatiotemporally aligned and quality-controlled multi-source long-term series data into the causal enhancement model of the reservoir to complete the missing time steps, thus obtaining complete multi-source long-term series data. The ecosystem structure and function evolution index system refers to a set of quantitative indicators used to characterize changes in the state of the ecosystem, covering biodiversity index, community stability index, primary productivity proxy quantity, and food web transmission efficiency.

[0032] The variance inflation factor (VIF) is a statistic that measures the degree to which a given independent variable in a multiple regression is linearly explained by other independent variables; a larger value indicates more severe collinearity. The VIF threshold is set to 10. This threshold is derived by simulating attribution bias under different degrees of collinearity using Monte Carlo methods on multiple estuary datasets. The VIF value corresponding to the first occurrence of an attribution bias rate exceeding 30% is used as a candidate threshold, and the average value is determined after 10 iterations. The two-stage collinearity decomposition framework refers to first orthogonalizing the candidate driving factors using principal component analysis to eliminate linear collinearity, then calculating the Shapley additive explanatory value of each candidate driving factor using an ensemble tree model in a low-dimensional orthogonal space, and finally mapping the Shapley additive explanatory value to the desired value through inverse transformation. The two-step attribution process projects back to the original physical variable space. The Shapley additive explanatory value is an application of Shapley values ​​in game theory to the interpretation of machine learning models. It is used to decompose the model output into the sum of the marginal contributions of each input feature, satisfying the axioms of efficiency, symmetry, and virtuality. The contribution of the original physical variable space is the quantified contribution share obtained by mapping the Shapley additive explanatory values ​​of each candidate driving factor to the original physical variable space after inverse transformation. It serves as one of the input sources for the contribution score in step S06. The two-stage collinearity decomposition framework solves the problem that the contribution of the original physical variable space is systematically underestimated due to coefficient compression in traditional regularization methods, and realizes interpretable quantitative decomposition of the independent contributions of each candidate driving factor under strong collinearity conditions.

[0033] The wavelet coherence analysis is a time-frequency analysis method used to extract the phase relationship and coherence degree between two sets of time series in a two-dimensional time and frequency plane. The phase lag spectrum is the output of the wavelet coherence analysis, representing the change matrix of the phase difference between candidate driving factors and ecological response variables at each frequency component over time, and serves as the input to adaptive dynamic time warping. The adaptive dynamic time warping is a sliding window variant of the dynamic time warping algorithm, which uses dynamic programming to search for the optimal alignment path between the two time series that minimizes the cumulative distance within different seasonal windows, thereby independently estimating the optimal alignment path within each seasonal window. The time lag; the optimal time lag is the optimal sequence alignment offset output by adaptive dynamic time warping within each seasonal window, which serves as the input of the observed sequence to the Hidden Markov Model; the Hidden Markov Model is a statistical model that describes the probabilistic relationship between the hidden state sequence and the observable sequence. Here, the optimal time lag within each seasonal window is used as the observed sequence to model the transition probability of multi-state time lags; the time-varying causal structure is the output of the Hidden Markov Model, representing the probability graph of the evolution of the causal relationship between each candidate driving factor and the ecological response variable over time, which serves as the input for the comprehensive inference of the reservoir causal enhancement model in step S06.

[0034] The ecological mechanism knowledge graph refers to a structured knowledge base built based on the generally accepted driver-response mechanism relationship in the field of ecology. It is used to forcibly retain first-order physical driving factors with clear physical and ecological significance before feature selection. The first-order physical driving factors are candidate driving factors marked in the ecological mechanism knowledge graph as directly affecting ecological response variables. The elastic network stability selection refers to repeatedly running elastic network regularized regression on multiple bootstrap subsamples and calculating the selection frequency of each candidate driving factor into the model. The selection frequency is the ratio of the number of times each candidate driving factor is selected into the model on the bootstrap subsamples to the total number of bootstrap iterations. The frequency threshold is set to 80%, and the source of the frequency threshold is: based on a simulated dataset... Based on the known true driving factors, the selection frequency range of 60% to 95% is traversed. The selection frequency value corresponding to the maximum harmonic mean of feature selection accuracy and recall is used as the candidate frequency threshold, which is determined through 5 independent experiments. The stable feature subset is the set of candidate driving factors whose selection frequency exceeds the frequency threshold from the stability selection output of the elastic network, and is used as the input for steps S05 and S06. The leave-one-region cross-validation refers to a validation method in which one region is reserved as an independent test set and the remaining regions are used as the training set for model training and evaluation in a dataset of multiple geographical regions. The generalization ability evaluation result is the prediction performance index on the independent test sets of each region output by leave-one-region cross-validation, which is used as one of the correction bases for the confidence score in step S06.

[0035] The driving factor importance ranking algorithm based on information geometry geodesic distance uses the Fisher information matrix of the statistical manifold in information geometry theory as the Riemannian metric tensor. It treats each candidate driving factor as a coordinate in a family of parameterized probability distributions, calculates the geodesic along the coordinate direction of each candidate driving factor using an exponential mapping, numerically integrates the curvature of each geodesic, and uses the integral value of the geodesic curvature as the driving factor importance score. A direction with greater geodesic curvature indicates a greater degree of bending in the ecological response distribution caused by a small change in the candidate driving factor. To reduce computational complexity, the algorithm employs... The Fisher information matrix is ​​decomposed using a low-rank approximation method. The low-rank number ranges from 5 to 30. The low-rank number is determined by Pareto front analysis of reconstruction error and computation time on simulated datasets of different dimensions, and the low-rank number corresponding to the inflection point of the Pareto front is selected through three iterative experiments. The geodesic curvature integral value is the output of the driving factor importance ranking algorithm based on information geometry geodesic distance, serving as the driving importance score for each candidate driving factor, and is used for the determination of the dominant driving factor and the grouping of composite driving factors in step S05. The importance ratio threshold is set to 1.5. This threshold is derived by iterating through a synthetic dataset of known dominant driving factors, traversing the ratio range from 1.1 to 3.0, and using the ratio corresponding to the first time the correct judgment rate of the dominant factor reaches 90% as a candidate threshold. This threshold is determined by averaging the results of five independent experiments. The dominant driving factor is the candidate driving factor corresponding to the highest driving importance score, serving as the primary component of the key driving factor list. The composite driving factor group is a set of candidate driving factors with similar driving importance scores that do not exceed the importance ratio threshold distinction standard, serving as one component of the key driving factor list. The driving factor importance ranking algorithm based on information geometry geodesic distance elevates the importance assessment of candidate driving factors from a Euclidean geometric framework of linear correlation or variance contribution to a Riemannian geometric framework. This framework can capture the influence of candidate driving factors on the nonlinear bending of the ecological response distribution, enabling robust ranking of key driving factors even in scenarios with strong nonlinearity and distribution shift in the driving-response relationship. This provides a theoretical basis with statistical manifold geometric significance for attribution conclusions.

[0036] The specific structure of the reservoir causal enhancement model is as follows: the reservoir causal enhancement model consists of four sequentially connected parts: a missing information input layer, an echo state network reservoir, a sparse causal transformer, and an output readout layer. The missing information input layer receives a mixed input of high-frequency continuous physical-driven time series and monthly-scale sparse ecological observations. Missing time steps are replaced with learnable missing embedding vectors, and a zero-weight mask is applied to the missing locations during loss calculation, preventing gradient updates to the predictions of missing locations during the training phase. The high-frequency continuous physical-driven time series refers to hourly sampled continuous observation sequences of tide level, salinity, and dissolved oxygen, serving as one of the inputs to the missing information input layer. The monthly-scale sparse ecological observations refer to monthly-scale sampled ecological observation records with missing data, also serving as one of the inputs to the missing information input layer. The echo-state network (ESN) reservoir consists of a large random sparse cyclic reservoir with 4000 to 6000 nodes, and a connectivity rate set between 0.5% and 2%. The node count and connectivity rate are determined by performing a grid search on multiple estuary datasets with the goal of minimizing the validation set prediction error. The search is conducted within a range of 3000 to 8000 nodes and a connectivity rate of 0.1% to 5%, and the parameter values ​​corresponding to the minimum error are determined after five independent experiments. The internal weights of the ESN reservoir are fixed after initialization and do not participate in training; only the linear weights of the output readout layer are trained. The ESN reservoir maps the output of the missing perceptual input layer to high-dimensional nonlinear temporal features, utilizing the memory properties of chaotic dynamics to preserve long-range temporal dependencies. These high-dimensional nonlinear temporal features serve as the input to a sparse causal transformer. The sparse causal transformer is a small transformer structure with 2 to 4 layers. A lower triangular causal mask is forcibly applied to the attention matrix to ensure the consistency of the temporal causal direction, and additional weights are applied to the attention weights. Regularization is used to generate a sparse attention graph. The non-zero edges of each attention head in the sparse attention graph are directly interpreted as the probability of the existence of edges between candidate driving factors and ecological response variables in the time-series causal graph. The layer number is determined by iterating through multiple estuarine datasets with the causal inference accuracy of the sparse attention graph as the target, within a range of 2 to 6 layers, and selecting the layer number corresponding to the highest accuracy after three independent experiments. The sparse attention graph, as the edge existence probability output of the time-series causal graph, is input together with the time-varying causal structure into the comprehensive inference stage of the reserve pool causal enhancement model in step S06. The output readout layer is a single-layer linear mapping that maps the output of the sparse causal transformer to the contribution score, stability score, and confidence score vector of each candidate driving factor.

[0037] Regarding resource allocation, a memory pressure index adjustment function is designed. This function calculates the memory pressure index based on the sequence length of the current batch, the number of nodes in the echo state network reserve pool, and the number of layers in the sparse causal transformer. The formula is as follows: ;in The current batch sequence length (unit: time steps). The length of the baseline sequence (unit: time step). The current number of nodes in the network reserve pool in the echo state. The number of reference nodes, This represents the current sparse causal transformer layer number. Baseline number of layers , , All parameters were determined by preliminary experiments on the target hardware when the video memory utilization rate reached exactly 70%; when At that time, one independent CUDA stream is allocated to both the echo state network reservoir and the sparse causal transformer, with the batch size remaining at the default value; when At that time, the computation of the echo state network reservoir and the computation of the sparse causal transformer are distributed to two parallel CUDA streams, and the batch size is reduced to 50% of the default value; when When the gradient checkpoint mechanism is activated, the number of nodes in the echo state network reserve pool is dynamically reduced to [a specific value]. And assign the three CUDA streams to matrix multiplication kernels at different levels; when When switching to half-precision floating-point operations, the batch size is reduced to 25% of the default value, and the CPU offload mechanism is enabled to temporarily store the output readout layer weights in memory.

[0038] The steps for establishing the training dataset for the causal enhancement model of the reservoir specifically include: extracting the high-frequency continuous physical driving time series and the monthly sparse ecological observation records from the historical observation databases of estuaries in multiple geographical regions, processing them uniformly through spatiotemporal alignment and quality control processes, and dividing them into source domain datasets and target domain datasets according to regions; the source domain dataset is used for meta-learning pre-training; only 10 to 20 labeled samples are retained in the target domain dataset for rapid adaptation, and the remaining samples are used as the test set; time series with a missing rate exceeding 40% are preprocessed with linear interpolation and mean filling before being processed by the missing rate-aware input layer; the 40% is derived by: traversing multiple estuary datasets with the goal of minimizing interpolation error, within a missing rate range of 30% to 60%, and determining the missing rate value corresponding to the first significant increase in interpolation error after 5 independent experiments.

[0039] The specific steps of training the reservoir causal enhancement model include: performing meta-learning pre-training on the source domain dataset using a model-independent meta-learning variant framework; calculating task gradients using a single region dataset as the task in the inner loop; and updating meta-parameters using the average gradients from multiple tasks in the outer loop; after meta-learning pre-training, fine-tuning the output readout layer and the sparse causal transformer using 10 to 20 labeled samples from the target domain dataset; freezing the weights within the echo state network reservoir; and the training loss consisting of the prediction error term and... The regularization terms are formed by weighted summation. Regularization weights are applied via grid search. to Within the range, the The regularization weights are derived as follows: with the goal of achieving the highest causal inference accuracy of sparse attention maps on multiple estuary datasets, a grid search is performed within the specified range, and the weight value corresponding to the highest accuracy is determined after three independent experiments.

[0040] The reservoir causal enhancement model captures the nonlinear long-range memory features of the high-frequency continuous physical-driven time series by fixing a large echo state network reservoir. Simultaneously, it infers a sparse temporal causal graph in the low-dimensional monthly-scale ecological observation space through a sparse causal mask attention mechanism, thus fusing the high-frequency continuous physical-driven time series with the cross-scale information of the monthly-scale sparse ecological observation within a unified end-to-end framework. The missing value-aware input layer avoids the interference of missing values ​​on gradient propagation, and the meta-learning pre-training endows the model with the ability to quickly adapt to new regions with a small number of target domain dataset samples. This allows the model to output a robust time-varying causal structure even with only 100 to 300 effective data points, fundamentally alleviating the overfitting problem in high-dimensional sparse feature spaces and providing a parameter-efficient adaptation path for cross-regional transfer and reuse.

[0041] Optionally, the present invention also provides a computer-based method for forming a key driver identification system for the structural and functional evolution of estuarine and adjacent marine ecosystems. The computer is equipped with a readable storage medium that stores program instructions, which execute the above-described method when the computer is run.

[0042] The specific implementation of step S01 is as follows: First, multi-source long-term series data from monthly to annual scales are collected for the study area. The data types include marine meteorological data, hydrodynamic data, water quality and environmental data, ecological observation data, terrestrial input data, and human activity intensity data. Due to differences in spatial and temporal resolution among data from different sources, bilinear interpolation or nearest-neighbor interpolation methods are used in the spatiotemporal alignment stage to uniformly interpolate all data to the same spatial grid and time step, ensuring that subsequent analysis is conducted in a unified spatiotemporal coordinate system. In the quality control stage, outlier removal and consistency checks are performed on each series to remove abnormal records introduced by instrument malfunctions or transmission errors. For time series with a missing rate of no more than 40%, the missing time step completion stage inputs the data into the causal enhancement model in the reserve pool, utilizing the model's sequence modeling capabilities to complete the missing positions. For time series with a missing rate exceeding 40%, preliminary preprocessing is performed using linear interpolation and mean imputation before being processed by the causal enhancement model in the reserve pool to reduce the interference of high missing rates on model inference. After preprocessing, an indicator system for the evolution of ecosystem structure and function was constructed based on complete multi-source long-term series data. This system includes four quantitative indicators: biodiversity index, community stability index, primary productivity proxy quantity, and food web transmission efficiency. These indicators serve as the source of ecological response variables in subsequent steps.

[0043] The specific implementation of step S02 is as follows: For each candidate driving factor in the preprocessed multi-source long-term series data, the variance inflation factor is calculated pairwise. The larger the variance inflation factor, the higher the degree to which the corresponding candidate driving factor is linearly explained by other candidate driving factors, i.e., the more severe the collinearity. When the variance inflation factor of any candidate driving factor exceeds the threshold of 10, a two-stage collinearity decomposition framework is initiated. In the first stage, principal component analysis is performed on all candidate driving factors to orthogonalize the original variable space, eliminating linear collinearity among candidate driving factors and obtaining uncorrelated low-dimensional orthogonal principal components. In the second stage, an ensemble tree model is constructed in the low-dimensional orthogonal space, and the marginal contribution decomposition of the model output is performed using Shapley additive explanatory values. Shapley additive explanatory values ​​satisfy the axioms of efficiency, symmetry, and virtuality, and can unbiasedly allocate the contribution of each orthogonal component in a game-theoretic sense. Finally, through the inverse transformation of principal component analysis, the Shapley additive explanatory values ​​in the orthogonal space are mapped back to the original physical variable space to obtain the original physical variable space contribution of each candidate driving factor, which serves as one of the input sources for the contribution score in step S06.

[0044] The specific implementation of step S03 is as follows: Wavelet coherence analysis is a time-frequency analysis method that can simultaneously characterize the phase relationship and coherence degree between two sets of time series in a two-dimensional plane of time and frequency. In this step, wavelet coherence analysis is performed pairwise on each candidate driving factor and ecological response variable in the preprocessed multi-source long-term series data, outputting the phase lag spectrum, which is the matrix showing the phase difference over time at each frequency component, reflecting the degree to which the candidate driving factor leads or lags the ecological response variable at different time scales. Using the phase lag spectrum as input, adaptive dynamic time warping searches for the optimal alignment path between the two time series that minimizes the cumulative distance within different seasonal windows using dynamic programming, thereby independently estimating the optimal time lag within each seasonal window and capturing the seasonal dependence of the driving factor's effect time lag. The optimal time-delay sequence output within each seasonal window is used as the observation sequence and input into the Hidden Markov Model. The Hidden Markov Model establishes a probabilistic graphical model for the implicit multi-state time-delay transition process, outputs the transition probability matrix between each state, forms a time-varying causal structure, and represents the probability graph of the evolution of the causal relationship between each candidate driving factor and the ecological response variable over time, which is used in step S06.

[0045] The specific implementation of step S04 is as follows: This step adopts a feature selection strategy that combines mechanistic constraints and data-driven approaches. First, based on the ecological mechanism knowledge graph, candidate driving factors labeled as first-order physical driving factors in the known driving response mechanism relationships are forcibly retained to ensure that variables with clear physical and ecological significance are not mistakenly eliminated by the data-driven selection process. Subsequently, elastic network stability selection is performed on all candidate driving factors. Elastic network regularized regression is repeatedly run on multiple bootstrap subsamples, and the selection frequency of each candidate driving factor into the model is counted. Candidate driving factors with a selection frequency exceeding a threshold of 80% are retained as a stable feature subset. This threshold is determined by maximizing the harmonic mean of feature selection accuracy and recall on the simulated dataset. After the stable feature subset is determined, leave-one-region cross-validation is used to evaluate its generalization ability. Each time, one geographical region is left as an independent test set, and the remaining regions are used as the training set. The prediction performance index on the independent test set of each region is output as the generalization ability evaluation result. This result is used as the basis for correcting the confidence score in step S06.

[0046] The specific implementation of step S05 is as follows: The driving factor importance ranking algorithm based on information geometric geodesic distance uses the Fisher information matrix of the statistical manifold as the Riemannian metric tensor, and regards each candidate driving factor as the coordinate axis direction in a family of parameterized probability distributions. For each candidate driving factor, a geodesic is generated on the statistical manifold through exponential mapping along its coordinate direction, and the curvature along the geodesic is numerically integrated to obtain the integral value of the geodesic curvature as the driving importance score. The larger the geodesic curvature, the greater the bending of the ecological response distribution caused by a small change in the candidate driving factor, that is, the more significant its driving effect. To reduce computational complexity, the Fisher information matrix is ​​decomposed into low-rank components using a low-rank approximation method. The low-rank number is determined by Pareto front analysis of reconstruction error and computation time, and its value ranges from 5 to 30. After ranking the driving importance scores of each candidate driving factor, the ratio of the highest score to the second highest score is calculated. If the ratio exceeds the importance ratio threshold of 1.5, the candidate driving factor corresponding to the highest score is determined as the dominant driving factor; otherwise, candidate driving factors with similar scores are grouped into the composite driving factor group.

[0047] The specific implementation of step S06 is as follows: The time-varying causal structure, driving importance score, and stable feature subset are used as inputs and fed into the reservoir causal enhancement model for comprehensive inference. The reservoir causal enhancement model consists of four sequentially connected parts: a missing-aware input layer, an echo-state network reservoir, a sparse causal transformer, and an output readout layer. The missing-aware input layer replaces missing time steps with learnable missing embedding vectors and applies zero-weight masks to missing positions during training to prevent missing values ​​from interfering with gradient propagation. The echo-state network reservoir has 4000 to 6000 nodes, a connectivity rate of 0.5% to 2%, fixed internal weights that do not participate in training, and maps the input to high-dimensional nonlinear temporal features, preserving long-range temporal dependencies. The sparse causal transformer has 2 to 4 layers, applies a lower triangular causal mask to the attention matrix, and uses... Regularization generates a sparse attention graph, where each non-zero edge directly corresponds to the probability of edge existence in the temporal causal graph. The output readout layer is a single-layer linear mapping that outputs the contribution score, stability score, and confidence score of each candidate driving factor, ultimately forming a list of key driving factors along with their duration, scale, and driving type.

[0048] It should be noted that the key technologies of this invention include: a two-stage collinearity decomposition framework that combines principal component analysis orthogonalization with the inverse transformation of Shapley additive explanatory values ​​to achieve unbiased decomposition of the independent contributions of each candidate driving factor under collinearity conditions without compressing the coefficients of the original variables, overcoming the inherent defect of traditional regularization methods where the contribution is systematically underestimated due to coefficient compression; a driving factor importance ranking algorithm based on information geometry geodesic distance that elevates the importance measurement from a linearly correlated Euclidean framework to a Riemannian geometric framework, using the geodesic curvature integral value on the statistical manifold to capture the nonlinear bending effect of candidate driving factors on the ecological response distribution, giving the ranking results in strongly nonlinear scenarios theoretical support of statistical manifold geometric meaning; and a reservoir causal enhancement model that captures the nonlinear long-range memory of high-frequency physical drivers by fixing a large echo state network reservoir, and endowing the model with the ability to quickly adapt to a small number of target domain samples through meta-learning pre-training. The synergistic effect of these three key technologies makes this invention complementary in three dimensions: collinearity resolution, nonlinear ranking, and sparse data inference, thereby achieving robust and interpretable quantitative attribution in complex scenarios with three coexisting challenges.

[0049] It should be noted that this invention also solves the following technical problem: In estuarine ecological attribution studies, the time lag of the effect of driving factors on ecological responses often changes significantly with the seasons. Traditional static correlation analysis assumes that the causal relationship is constant and cannot capture this seasonally dependent time-varying time lag characteristic, leading to systematic biases in attribution conclusions under different seasonal scenarios. This invention extracts the phase lag spectrum in the time-frequency two-dimensional plane through wavelet coherence analysis, then independently estimates the optimal time lag within each seasonal window using adaptive dynamic time warping, and finally models the multi-state time lag transition probability using a hidden Markov model to form a probability graph that reflects the time-varying causal structure. Thus, without assuming a constant causal relationship, it achieves dynamic modeling of seasonally dependent time-varying time lags, solving the technical problem that traditional static correlation analysis cannot capture the dynamic evolution of the time lag of the driving factors with the seasons.

[0050] Specifically, the principle of this invention is: The fundamental reason why this invention can solve the above-mentioned technical problems is that its solution designs theoretically supported processing mechanisms for the three challenges, and the mechanisms are logically connected.

[0051] First, the two-stage collinearity decomposition framework utilizes the orthogonality of principal component analysis to map collinear variables to an uncorrelated low-dimensional space. The Shapley additive explanatory values ​​calculated in this space satisfy the axioms of efficiency and symmetry, enabling unbiased decomposition of the marginal contributions of each orthogonal component. An inverse transformation then restores the physical variable space, theoretically guaranteeing the interpretability of independent contributions under collinearity conditions. Second, the information geometry geodesic distance ranking algorithm elevates the importance assessment of candidate factors from an Euclidean geometric framework to a Riemannian geometric framework. It uses the geodesic curvature integral value on the statistical manifold to measure the nonlinear bending degree of each factor's response to the ecological response distribution, providing a geometrically meaningful theoretical basis for importance ranking in strongly nonlinear scenarios. Finally, the reservoir causal enhancement model captures high-frequency, physically driven, nonlinear long-range memory by fixing a large echo state network reservoir. Simultaneously, it infers a sparse temporal causal graph in a low-dimensional monthly-scale observation space using a sparse causal mask attention mechanism. Meta-learning pre-training endows the model with the ability to adapt quickly to a small number of target domain samples, alleviating the overfitting problem under high-dimensional sparse conditions from a parameter efficiency perspective. The three elements work together to form a logically consistent and theoretically verifiable complete attribution framework.

[0052] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0053] The specific implementation of step S01 is as follows: Acquire multi-source long-term series data of the study area at monthly to annual scales, covering marine meteorological data, hydrodynamic data, water quality and environmental data, ecological observation data, terrestrial input data, and human activity intensity data. Preprocessing includes three stages: spatiotemporal alignment, quality control, and missing time step completion. Spatiotemporal alignment unifies the interpolation of data from different sources and at different resolutions to the same spatial grid and time step. Quality control removes outlier observations. Missing time step completion inputs the aligned data into a causal enhancement model in the buffer pool to fill in the missing data, obtaining complete multi-source long-term series data. Based on this, an indicator system for the evolution of ecosystem structure and function is constructed, covering biodiversity index, community stability index, primary productivity proxy quantity, and food web transmission efficiency.

[0054] The specific implementation of step S02 is as follows: Calculate the variance inflation factor among all candidate driving factors, setting the variance inflation factor threshold to 10. When the variance inflation factor of a candidate driving factor exceeds the threshold, initiate the two-stage collinearity decomposition framework. In the first stage, for... A data matrix consisting of candidate driving factors ( For the sample size, Principal component analysis was performed to orthogonalize the covariance matrix (to the number of candidate driving factors), and the covariance matrix was decomposed into:

[0055] In the formula, Let covariance matrix be the variance matrix. The eigenvector matrix, It is a diagonal eigenvalue matrix. ( ) is the first 1 eigenvalue, for The transpose of the matrix. Take the first part. One principal component, yielding a low-dimensional orthogonal projection matrix. ,in for The former The truncated matrix formed by the column eigenvectors, This is the orthogonalized low-dimensional representation. The second phase, in Spatially, the Shapley additive interpretation values ​​for each principal component direction are calculated using an ensemble tree model. ( Then, through an inverse transformation, it is mapped back to the original physical variable space to obtain the first... Spatial contribution of each candidate driving factor to the original physical variables:

[0056] In the formula, for No. Line number Column elements, For the first Spatial contribution of original physical variables of each candidate driving factor ( ), dimensions and Both are dimensionless contribution shares, which serve as one of the inputs for contribution scoring in step S06.

[0057] The specific implementation of step S03 is as follows: wavelet coherence analysis is used to analyze the time series of candidate driving factors. Time series of ecological response variables Extract the phase lag spectrum in the time-frequency plane. The continuous wavelet transform is defined as:

[0058] In the formula, For the mother wavelet function, It is a scale parameter (inversely proportional to frequency, dimensionless). This is the translation parameter (with the same dimensions as time, and the unit is time step). express The complex conjugate, for In scale Translation parameters Wavelet transform coefficients at position, dimensions and same. Definition and Completely similar Replace with time series of ecological response variables That is, dimensions and Same. Cross-wavelet spectrum is defined as... The upper horizontal line indicates complex conjugation. Dimensions are and Product of dimensions. The wavelet coherence coefficient is defined as:

[0059] In the formula, For time-frequency smoothing operators, Represents the modulus of a complex number. These are dimensionless wavelet coherence coefficients, with dimensions of both the numerator and denominator. and The square of the product of dimensions, divided by the product, yields a dimensionless quantity. A value closer to 1 indicates a higher degree of coherence between the two sequences at that time-frequency point. Phase lag spectrum matrix. The The elements are:

[0060] In the formula, For the first The scale corresponding to each frequency component ( ), Take a complex argument, and output the unit in radians (rad). The total number of frequency components. The total number of time steps. For the first The frequency component is at the ... The phase difference at each time step is expressed in radians. As input for adaptive dynamic time warping, in each seasonal window Within this framework, dynamic programming is used to search for the alignment path that minimizes the cumulative distance between two time series, and the optimal time delay is estimated. :

[0061] In the formula, To align the path, For path Index pairs on, This represents the Euclidean distance (in radians) of the phase lag spectrum at the corresponding index. For seasonal windows The optimal time delay estimate (in time steps) within the time window. The optimal time delay sequence for each seasonal window. The observation sequence, as a hidden Markov model, models the multi-state time-delay transition probability and outputs a time-varying causal structure probability map, which serves as the input for step S06.

[0062] The specific implementation of step S04 is as follows: First-order physical driving factors are forcibly retained based on the ecological mechanism knowledge graph, and then the stability of the elastic network is selected... The elastic net regularized regression is repeatedly run on the bootstrap subsamples. Solving on the second-order bootstrap sample:

[0063] In the formula, For the first Number of samples in the second bootstrap subsample For the response variable vector, For the first Candidate driving factor matrix of the second bootstrap subsample For the regression coefficient vector, for Regularization coefficient (dimensionless). for Regularization coefficient (dimensionless). for Norm, For the first The elastic network estimation coefficient vector on the second bootstrap subsample. Statistical analysis of the first... Frequency of selection of each candidate driving factor :

[0064] In the formula, This represents the total number of bootstrapping attempts. For the first The feature subset selected from the elastic net on the secondary bootstrap subsamples. For indicator functions, The selected frequencies are dimensionless. The frequency threshold was experimentally determined to be 80%, and these were retained. The candidate driving factors constitute a stable feature subset. The generalization ability of the stable feature subset is evaluated by leave-one-out cross-validation, and the prediction performance index on the independent test set of each region is output as one of the bases for correcting the confidence score in step S06.

[0065] The specific implementation of step S05 is as follows: Ranking the candidate driving factors in the stable feature subset using a driving factor importance ranking algorithm based on information geometric geodesic distance. This is done using the Fisher information matrix. ( (The parameter vector of the parameterized probability distribution family) is the Riemann metric tensor, for Perform low-rank approximate decomposition:

[0066] In the formula, For the front A matrix composed of singular vectors. For the corresponding singular value diagonal matrix, ( ) is the first A singular value, a low-rank number. The value ranges from 5 to 30, and is determined at the inflection point of the Pareto front. for The transpose of . The Christofel notation derived from the Fisher information matrix is ​​defined as:

[0067] In the formula, for The inverse matrix of the first element, for The element, For parameter vectors The One portion, for right The partial derivative, Christofel notation It is a dimensionless geometric quantity. The first is calculated using an exponential mapping. Coordinates of candidate driving factors geodesic lines on It satisfies the geodesic equation:

[0068] In the formula, geodesic The Each coordinate component For dimensionless curve parameters, the dimensions of all terms on the left-hand side of the equation are the dimensions of the second derivative of the curve parameters, ensuring dimensional consistency. Numerical integration of the geodesic curvature yields the... The driving importance score of each candidate driving factor :

[0069] In the formula, For The Riemann norm is a metric. The dimensionless integral value of curvature. A higher score indicates a greater degree of bending in the ecological response distribution caused by a small change in the candidate driving factor. The highest driving factor importance score is assigned accordingly. With the second highest driver importance score Comparison of the ratio with the importance ratio threshold of 1.5: If If the condition is met, the dominant driving factor will be output; otherwise, a group of composite driving factors will be output.

[0070] The specific implementation of step S06 is as follows: The reservoir causal enhancement model consists of four sequentially connected parts: a missing-aware input layer, an echo-state network reservoir, a sparse causal transformer, and an output readout layer. The missing-aware input layer replaces missing time steps with learnable missing embedding vectors. During loss calculation, a zero-weight mask is applied to the missing positions to ensure that no gradient updates occur at the missing positions. The internal state update equation of the echo-state network reservoir is:

[0071] In the formula, for The state vector of the reserve pool at any given time. For the number of nodes in the reserve pool, Leakage rate (dimensionless). The hyperbolic tangent activation function is used. For the input weight matrix, The input feature dimension, This is the internal weight matrix of the reserve pool (fixed after initialization and not trained). for Input vector at time step, for The time-order noise vector is both sides of the equation. 3D vector. The sparse causal transformer applies a lower triangular causal mask to the attention matrix. ( (where the sequence length is in time steps), and applies it to the attention weights. Regularization, the total training loss is:

[0072] In the formula, This is the prediction error term (dimensionless normalization loss). for Regularization weights (dimensionless, grid search range) to ), For the number of converter layers, For the first Layer attention matrix, for The zero norm is the number of non-zero elements (dimensionless). For dimensionless total loss, non-zero edges in the sparse attention graph are directly interpreted as the probability of edge existence in the temporal causal graph. Each regional task updates its task parameters using gradient descent:

[0073] In the formula, For the meta-parameter vector, The learning rate (dimensionless step size) is the inner loop learning rate. For the first The loss function value (dimensionless) for each task. The gradient vector of the loss with respect to the meta-parameters. For the first The parameter vector after task adaptation, the dimensions of both sides of the equation are all consistent with... Same. The outer loop updates the meta-parameters using the mean of gradients from multiple tasks:

[0074] In the formula, The outer loop learning rate (dimensionless step size). The total number of tasks in the source domain. In order to be in For meta-parameters The obtained gradient vector has dimensions on both sides of the equation that are the same as... Same. The formula for the memory pressure index adjustment function is:

[0075] In the formula, The current batch sequence length (unit: time steps). The length of the baseline sequence (unit: time step). The number of nodes in the baseline reserve pool, The three baseline values ​​are determined by pre-experimentation on the target hardware when the video memory utilization rate reaches exactly 70%, representing the number of base converter layers. This is a dimensionless memory pressure index, where all three factors are ratios of the same type of quantity. When At that time, each is assigned an independent computation stream, and the batch size remains at the default value; when When two parallel computing streams are allocated, the batch size is reduced to 50% of the default value; when When gradient checkpoints are enabled, the number of nodes in the reserve pool is dynamically reduced to [a specific value]. Allocate 3 computational streams; when At this time, half-precision floating-point operations are switched, the batch size is reduced to 25% of the default value, and the CPU offload mechanism is enabled. The output readout layer is a single-layer linear mapping that maps the output of the sparse causal transformer to the contribution score, stability score, and confidence score vector of each candidate driving factor. Combined with the time-varying causal structure and stable feature subset, a list of key driving factors and their period of action, scale of action, and driving type are finally formed.

[0076] To better understand and implement this invention, the following is a specific application scenario example 2: Technicians set up a test environment, selected a typical tidal estuary as the research object, and used multi-source historical observation data from 2005 to 2023 for a total of 19 years as input. The key driving factors of the evolution of the estuary ecosystem structure and function were identified using the method of this invention.

[0077] In step S01, technicians collected monthly marine meteorological data (including wind speed, air temperature, and precipitation), hydrodynamic data (including tide level, runoff, and salinity), water quality data (including dissolved oxygen, inorganic nitrogen, phosphate, and chlorophyll a), ecological observation data (including phytoplankton abundance and benthic biodiversity), land-based input data (including total nitrogen and total phosphorus fluxes carried by runoff into the sea), and human activity intensity data (including the rate of change in reclaimed area and aquaculture density index) covering 19 years. Spatiotemporal alignment was performed using bilinear interpolation, uniformly interpolating all data to a 0.05°×0.05° spatial grid and a monthly time step. After quality control, the overall missing rate of the ecological observation data was approximately 28%, lower than the 40% preprocessing threshold, and the missing time steps were directly filled in by the reserve pool causal enhancement model. The constructed ecosystem evolution index system is shown in Table 1.

[0078] Table 1. Indicator System for Ecosystem Structure and Function Evolution

[0079] In step S02, technicians calculated the variance inflation factor for each of the 12 candidate driving factors. The results showed that the variance inflation factor between runoff and total nitrogen flux into the sea reached 14.3, exceeding the threshold of 10, indicating strong collinearity between the two. A two-stage collinearity decomposition framework was then initiated. Principal component analysis orthogonalized the 12 candidate driving factors into 8 principal components. An ensemble tree model was constructed in a low-dimensional orthogonal space, and Shapley additive explanatory values ​​were calculated. Then, an inverse transformation was performed to restore the original physical variable space, obtaining the original physical variable space contribution of each candidate driving factor, such as... Figure 2 As shown (the spatial distribution of the original physical variables of each candidate driving factor).

[0080] In step S03, wavelet coherence analysis is used to perform time-frequency analysis on runoff and chlorophyll a concentration. Significant coherence is observed between the two in the 8-16 month period. The phase difference indicates that runoff leads chlorophyll a concentration by approximately 1 to 3 months, with this lag being approximately 1 month in the summer window and approximately 3 months in the winter window, demonstrating a clear seasonal dependence. Adaptive dynamic time warping estimates the optimal lag within each seasonal window. After modeling the multi-state time-lag transition probability using a hidden Markov model, a time-varying causal structure probability diagram is output, as shown below. Figure 3 As shown. Figure 3 As shown, the driving effect of runoff on chlorophyll a is in a short time lag state in summer and switches to a long time lag state in winter. The state transition probability is given by the hidden Markov model, which is consistent with the understanding of the ecological mechanism of seasonal hydrological rhythm.

[0081] In step S04, based on the ecological mechanism knowledge graph, runoff, tidal mixing intensity, and solar radiation are forcibly labeled as first-order physical driving factors and retained. The stability of the elastic network is determined by statistically analyzing the selection frequency of each candidate driving factor across 500 bootstrapping subsamples; the results are shown in Table 2.

[0082] Table 2. Selection Frequency of Candidate Driving Factors in Elastic Network Stability Selection

[0083] Five candidate driving factors with a selection frequency exceeding 80% were included in the stable feature subset. Leave-one-out cross-validation had a mean root mean square error of 0.18 across four independent test regions, indicating good generalization ability.

[0084] In step S05, the driving factor importance ranking algorithm based on information geometry geodesic distance calculates the geodesic curvature integral value for the five candidate driving factors. The Fisher information matrix uses a low-rank approximation with a rank of 12, which is determined by Pareto front analysis after balancing reconstruction error and computation time. The driving importance scores of each candidate driving factor are as follows: Figure 4 As shown in the ranking chart of candidate driving factors by importance score, runoff has the highest driving importance score, followed by total nitrogen flux into the sea. The ratio of the two is 1.38, which is lower than the importance ratio threshold of 1.5. Therefore, runoff and total nitrogen flux into the sea together constitute a composite driving factor group, indicating that in this estuary, runoff and nutrient input synergistically drive changes in primary productivity, and there is no single dominant factor.

[0085] In step S06, the time-varying causal structure, driving importance score, and stable feature subset are input into the reservoir causal enhancement model for comprehensive inference. The number of reservoir nodes in the echo state network of the reservoir causal enhancement model is set to 5000, the connectivity rate is set to 1%, and the number of sparse causal transformer layers is set to 3. The regularization weights were determined by grid search. The model was pre-trained using historical data from four source regions, and then fine-tuned using 15 labeled samples from the target region to optimize the output readout layer and the sparse causal transformer. The overall inference results are shown in Table 3.

[0086] Table 3. Comprehensive scoring results of key driving factors

[0087] The final list of key driving factors indicates that runoff and total nitrogen flux into the sea are the key composite driving factors for the evolution of primary productivity in this estuary, with an effect spanning the entire year, ranging from seasonal to interannual scales, and a continuous driving type. Tidal mixing intensity and solar radiation are seasonal auxiliary driving factors, with their effects mainly concentrated in spring and summer. Phosphate concentration is an intermittent driving factor, with its effect associated with extreme precipitation events. The sparse attention map output by the sparse causal transformer is shown below. Figure 5 As shown in the sparse attention graph of the sparse causal transformer, each non-zero edge corresponds to a temporal causal relationship, and the edge weight reflects the probability of the existence of the causal relationship. This visually illustrates the temporal causal network structure between each candidate driving factor and the ecological response variable. The dynamic changes in the contribution scores of each candidate driving factor in different seasons are shown below. Figure 6 As shown in the figure (dynamic change of seasonal contribution scores of each candidate driving factor), the relative contribution weight of the composite driving factor group fluctuates significantly in different seasons, further verifying the necessity of time-varying causal structure modeling.

[0088] Compared to traditional multiple linear regression and static correlation attribution methods, this invention achieves three advancements in technical principles: First, the two-stage collinearity decomposition framework, through the combined use of orthogonalization and inverse transformation, avoids the confusion of independent contributions of each candidate driving factor by coefficient compression, enabling runoff and total nitrogen flux into the sea to still obtain interpretable independent contributions under strong collinearity conditions; Second, the time-varying causal structure modeling link enables the method to distinguish between two driving states: short time lag in summer and long time lag in winter, while traditional methods, due to the assumption of constant time lag, cannot detect this seasonal difference; Third, the meta-learning pre-training mechanism allows the model to complete rapid adaptation with only 15 labeled samples in the target domain, effectively suppressing the risk of overfitting in high-dimensional sparse feature spaces from the perspective of parameter efficiency, and ensuring the generalization ability during cross-regional migration.

[0089] It should be noted that the variables involved in this invention are explained in detail in Tables 4 and 5.

[0090] Table 4. Variable Explanation Table (Part 1)

[0091] Table 5. Variable Explanation Table (Part Two)

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

Claims

1. A method for identifying key driving factors in the structural and functional evolution of estuarine and adjacent marine ecosystems, characterized in that, Includes the following steps: We acquired multi-source long-term series data of the study area, preprocessed the multi-source long-term series data, and constructed an indicator system for the evolution of ecosystem structure and function based on the preprocessed multi-source long-term series data. Calculate the variance inflation factor among all candidate driving factors in the preprocessed multi-source long-term series data, compare the variance inflation factor with the variance inflation factor threshold, and when the variance inflation factor exceeds the variance inflation factor threshold, initiate a two-stage collinearity decomposition framework to obtain the spatial contribution of the original physical variables. Wavelet coherence analysis is used to extract the phase lag spectrum from the preprocessed multi-source long-time series data. The optimal time lag is estimated in different seasonal windows by combining adaptive dynamic time warping. Hidden Markov model is used to integrate the multi-state time lag transition probabilities to form a time-varying causal structure. Based on the knowledge graph of ecological mechanism, the first-order physical driving factors are forcibly retained. The selection frequency of each candidate driving factor is statistically selected by the stability selection of elastic network. The candidate driving factors with selection frequency exceeding the frequency threshold are retained as stable feature subsets. The stable feature subsets are evaluated by leave-one-out cross-validation to obtain the generalization ability evaluation results. The candidate driving factors in the stable feature subset are ranked by the driving factor importance ranking algorithm based on information geometry geodesic distance. The geodesic curvature integral value is calculated as the driving importance score. The ratio of the highest driving importance score to the second highest driving importance score is compared with the importance ratio threshold to obtain the dominant driving factor or composite driving factor group. The time-varying causal structure, driving importance score, and stable feature subset are input into the reserve pool causal enhancement model for comprehensive inference, and the contribution score, stability score and confidence score of each candidate driving factor are output to form a list of key driving factors and their period of action, scale of action and driving type. Specifically, the two-stage collinearity decomposition framework first uses principal component analysis to orthogonalize candidate driving factors to eliminate linear collinearity, then uses an ensemble tree model to calculate the Shapley additive explanatory value of each candidate driving factor in a low-dimensional orthogonal space, and finally maps the Shapley additive explanatory value back to the original physical variable space through inverse transformation. The reservoir causal enhancement model consists of four sequentially connected parts: a missing perception input layer, an echo state network reservoir, a sparse causal transformer, and an output readout layer. The missing perception input layer receives a mixed input of high-frequency continuous physical driving time series and monthly-scale sparse ecological observations, replaces missing time steps with learnable missing embedding vectors, and applies a zero-weight mask to the missing positions during loss calculation. The echo state network reservoir is composed of a large random sparse cyclic reservoir. Its internal weights are fixed after initialization and do not participate in training; only the linear weights of the output readout layer are trained. The attention matrix of the sparse causal transformer is forcibly subjected to a lower triangular causal mask, and additional causal masks are applied to the attention weights. Regularization is used to generate sparse attention maps.

2. The method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems according to claim 1, characterized in that, The multi-source long-term series data includes marine meteorological data, hydrodynamic data, water quality and environmental data, ecological observation data, land-based input data, and human activity intensity data. The preprocessing includes three steps: spatiotemporal alignment, quality control, and missing time step completion.

3. The method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems according to claim 2, characterized in that, The missing time step completion refers to inputting multi-source long-term series data, which has undergone spatiotemporal alignment and quality control, into the causal enhancement model of the reserve pool to complete the missing time steps and obtain complete multi-source long-term series data.

4. The method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems according to claim 3, characterized in that, The phase lag spectrum is the output of wavelet coherence analysis, which represents the change matrix of the phase difference between candidate driving factors and ecological response variables at each frequency component over time, and serves as the input for adaptive dynamic time warping.

5. The method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems according to claim 4, characterized in that, The adaptive dynamic time warping specifically involves using dynamic programming to search for the optimal alignment path between two time series that minimizes the cumulative distance within different seasonal windows, thereby independently estimating the optimal time delay within each seasonal window. The optimal time delay serves as the input of the observation sequence to the hidden Markov model.

6. The method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems according to claim 5, characterized in that, The stability selection of the elastic network specifically involves repeatedly running the elastic network regularized regression on multiple bootstrap subsamples and calculating the selection frequency of each candidate driving factor into the model. The selection frequency is the ratio of the number of times each candidate driving factor is selected into the model to the total number of bootstraps.

7. The method for identifying key driving factors of the evolution of the structure and function of estuary and adjacent marine ecosystems according to claim 6, characterized in that, The driving factor importance ranking algorithm based on information geometry geodesic distance uses the Fisher information matrix as the Riemann metric tensor, treats each candidate driving factor as a coordinate in a family of parameterized probability distributions, calculates the geodesic along the coordinate direction of each candidate driving factor using exponential mapping, numerically integrates the curvature of each geodesic, and uses the integral value of the geodesic curvature as the driving factor importance score.

Images