Lake water quality change prediction and early warning method based on complex relationship between lake current and water quality

Abstract

The application discloses a lake water quality change forecasting and early warning method based on a complex relationship between lake flow and water quality. The method obtains multi-source monitoring data of a lake area, calculates a lake flow comprehensive characteristic index containing flow velocity, vorticity and hydraulic retention time, constructs a conventional early warning path, identifies the optimal time lag of the lake flow and the water quality, adopts a physical constraint kernel regression with transport and retention constraints to construct a response model, and determines a dynamic early warning threshold that changes with the lake flow condition. A mutation early warning path is constructed, a critical slowing down index of the water quality sequence is calculated, a spatial propagation weight is constructed based on the flow direction of the flow field, and an early warning index with enhanced early warning is synthesized by superimposing the upstream signal on the downstream. The early warning levels of the two paths are fused and decided. The application solves the problems of data model violation of physical laws and single-point early warning lag, and realizes early warning of gradual and sudden water quality risks.

Description

Technical Field

[0001] This invention relates to the field of water environment monitoring and early warning technology, and in particular to a method for forecasting and early warning of water quality changes in lake areas based on the complex relationship between lake flow and water quality. Background Technology

[0002] Lake ecosystems are characterized by nonlinearity and dynamic complexity, with hydrodynamic conditions driving water quality evolution and algal blooms. Analyzing the spatiotemporal coupling mechanism between lake flow fields and water quality indicators, and exploring the nonlinear response of flow field transport to pollutant concentrations, will enable the construction of high-precision forecasting and early warning models. This will help reveal the mechanisms of pollutant migration and transformation, capture early signals of the critical transition from clear to turbid water in the ecosystem, and achieve proactive management of the lake water environment.

[0003] Currently, research and application of lake water quality early warning mainly rely on two technical approaches: one is based on physical mechanism-based hydrodynamic-water quality coupling models, such as the Environmental Fluid Dynamics Model (EFDC) and the Delft 3D hydrodynamic model, which simulate flow fields and material transport processes by numerically solving governing equations; the other is based on data-driven machine learning models, such as using Long Short-Term Memory (LSTM) networks or Support Vector Machines (SVM) to mine statistical patterns in water quality time series. Regarding early warning indicators, existing technologies set concentration thresholds based on national surface water environmental quality standards, or calculate indicators such as variance and autocorrelation coefficients based on the Critical Slowing Theory (CSD) of a single site to identify the risk of sudden changes in water quality.

[0004] Existing technologies mainly suffer from prediction distortion due to the lack of physical mechanism constraints and early warning lag due to insufficient spatial linkage. Therefore, further research and innovation are needed to solve the aforementioned problems of existing technologies. Summary of the Invention

[0005] The purpose of this invention is to propose a method for forecasting and early warning of water quality changes in lake areas based on the complex relationship between lake flow and water quality, in view of the above-mentioned problems existing in the prior art.

[0006] According to one aspect of this application, a method for forecasting and early warning of lake water quality changes based on the complex relationship between lake currents and water quality includes:

[0007] Acquire multi-source monitoring data for the lake area, which should include at least lake flow monitoring data and water quality monitoring data;

[0008] Based on lake flow monitoring data, a comprehensive characteristic index of lake flow reflecting the hydrodynamic conditions of the lake area is calculated;

[0009] Identify the time lag relationship between the comprehensive characteristic index of lake flow and water quality monitoring data, construct a lake flow-water quality response model, determine the dynamic early warning threshold that changes with the comprehensive characteristic index of lake flow based on the lake flow-water quality response model, and determine the regular early warning level based on the dynamic early warning threshold;

[0010] A sliding window analysis was performed on the time series of water quality monitoring data to calculate the critical slowing index. An early warning comprehensive index was synthesized based on the critical slowing index, and the level of sudden change warning was determined based on the early warning comprehensive index.

[0011] The system integrates regular warning levels and sudden change warning levels to generate and output water quality warning information.

[0012] Beneficial effects: This invention solves the problems of data models violating physical laws and single-point early warning lag, and realizes early warning for both gradual and abrupt water quality risks. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0013] Figure 1 A flowchart illustrating a method for forecasting and early warning of lake water quality changes based on the complex relationship between lake currents and water quality, provided for embodiments of this application.

[0014] Figure 2 A flowchart illustrating the time lag relationship between the comprehensive characteristic index of lake flow and water quality monitoring data, provided for embodiments of this application.

[0015] Figure 3 The flowchart for constructing a lake flow-water quality response model is provided for an embodiment of this application.

[0016] Figure 4 This is a flowchart illustrating the dynamic early warning threshold determined by the lake current-water quality response model based on the lake current-water quality response model, as provided in this application embodiment.

[0017] Figure 5 This is a flowchart illustrating the fusion decision-making process for conventional early warning levels and mutation early warning levels, provided as an embodiment of this application. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present invention, 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. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0019] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in sequences other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0020] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0021] Correspondingly, pure data-driven models, under extreme hydrological conditions with insufficient training sample coverage, such as exceptionally low or high water levels, lack the constraints of physical laws and are prone to outputting predictions that violate the laws of mass conservation or transport and diffusion, such as abnormally high concentrations when the flow velocity is high. Mechanistic models, on the other hand, are limited by the uncertainty of parameter calibration.

[0022] Furthermore, existing abrupt change warnings are mostly limited to single-point time series analysis, severing the flow field connection between upstream and downstream areas of lakes. This prevents downstream areas from utilizing upstream precursor signals for early warning, triggering alarms only after pollutants have already migrated to their destination, resulting in poor timeliness. In addition, fixed threshold warnings ignore the dynamic changes in water body environmental capacity under different hydrodynamic conditions, making it difficult to adapt to the complex hydrological rhythms of lakes.

[0023] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0024] On the one hand, this paper provides an exemplary scheme for forecasting and early warning of lake water quality changes based on the complex relationship between lake currents and water quality, solving the technical problem that traditional single water quality early warning methods cannot simultaneously address long-term trend monitoring and the risk of sudden critical transitions. Specifically, this scheme includes:

[0025] Step 101: Obtain multi-source monitoring data for the lake area. The multi-source monitoring data includes at least lake flow monitoring data and water quality monitoring data.

[0026] The multi-source monitoring data can include, but is not limited to, automatic monitoring stations deployed in the lake area, shipborne ADCP (Acoustic Doppler Current Profiler) monitoring, satellite remote sensing data, and meteorological station observations. Lake current monitoring data mainly includes velocity vector fields and water levels. The velocity vector field characterizes the magnitude and direction of the flow velocity at various spatial locations within the lake at corresponding times; water level data is typically collected in real-time by water level gauges. Water quality monitoring data mainly includes key indicators reflecting the eutrophication state of the lake, such as the concentration time series of total phosphorus, total nitrogen, chlorophyll a, and permanganate index.

[0027] Additionally, meteorological monitoring data, such as wind speed, wind direction, temperature, and precipitation, can be acquired simultaneously to assist in the analysis of changes in hydrodynamic conditions. After acquiring the raw data, preprocessing is required, including removing outliers, imputing missing values, and unifying data from different sources to the same time resolution, such as a daily scale.

[0028] Step 102: Based on lake flow monitoring data, calculate the comprehensive characteristic index of lake flow that reflects the hydrodynamic conditions of the lake area.

[0029] Furthermore, the Lake Flow Integrated Characteristic Index (LCI) is a dimensionless index that comprehensively characterizes the local and overall hydrodynamic intensity of a lake area. Since a single flow velocity index is insufficient to fully reflect the transport and retention characteristics of pollutants, this invention introduces three physical components: flow velocity, eddy current, and hydraulic residence time.

[0030] Specifically, based on lake flow monitoring data, the standardized values ​​of flow velocity, flow field vorticity and its standardized value, and hydraulic residence time and its standardized value calculated based on the particle tracking method are determined for each monitoring point or grid point. Next, a weighted synthesis method is used to obtain the comprehensive characteristic index of the lake flow. A higher value indicates stronger hydrodynamic conditions, faster water exchange, and a stronger ability to dilute and disperse pollutants; conversely, a lower value indicates that the water body tends to be stagnant or stagnant, and pollutants are more likely to accumulate.

[0031] Step 103: Identify the time lag relationship between the lake flow comprehensive characteristic index and water quality monitoring data, construct a lake flow-water quality response model, determine the dynamic early warning threshold that changes with the lake flow comprehensive characteristic index based on the lake flow-water quality response model, and determine the regular early warning level based on the dynamic early warning threshold.

[0032] Alternatively, the time lag relationship between the comprehensive characteristic index of the lake flow and the water quality monitoring data can be identified to obtain the optimal time lag. A lake flow-water quality response model can be constructed based on this model to determine the dynamic early warning threshold that changes with the comprehensive characteristic index of the lake flow. The real-time water quality monitoring data can be compared with the dynamic early warning threshold, and the conventional early warning level can be determined based on the comparison results.

[0033] Since the impact of lake current changes on water quality is not instantaneous but exhibits a certain lag effect, it is necessary to identify the optimal time lag between the two. Using time-lag cross-correlation analysis, the number of days in which lake current changes lead water quality changes can be determined. Based on this time-lag relationship, a response model describing the change in water quality concentration with the lake current comprehensive characteristic index is constructed.

[0034] Based on this, instead of using traditional fixed concentration thresholds, such as those based solely on national standards, the warning thresholds are dynamically adjusted according to current lake flow conditions. For example, when the comprehensive characteristic index of the lake flow is high, the water body has a strong self-purification capacity, and the warning threshold can be appropriately relaxed; when the comprehensive characteristic index of the lake flow is low, the risk of water retention is high, and the warning threshold is correspondingly tightened. Furthermore, real-time water quality monitoring data is compared with the currently calculated dynamic warning thresholds to determine the conventional warning level, such as normal, blue warning, yellow warning, or red warning.

[0035] Step 104: Perform sliding window analysis on the time series of water quality monitoring data, calculate the critical slowing index, synthesize an early warning comprehensive index based on the critical slowing index, and determine the sudden change warning level based on the early warning comprehensive index.

[0036] In this step, the approach is based on the critical slowing principle in complex systems theory, used to capture early signals before abrupt changes in ecosystem states, such as algal blooms or rapid water quality deterioration. Specifically, a recent segment of water quality monitoring data is extracted using a sliding window, and the changing trend of its statistical characteristics, i.e., the critical slowing index, is calculated. This index mainly includes the variance growth rate, the autocorrelation enhancement rate, and the skewness change rate.

[0037] As the system approaches a critical point, these indicators typically show an increasing trend. Furthermore, individual indicators are combined into an Early Warning System (EWS) that reflects the overall risk of mutation. Based on the magnitude of this indicator or its degree of anomalousness relative to historical benchmarks, the mutation warning level is determined, such as attention, alert, or mutation warning.

[0038] Step 105: Perform a fusion decision on the regular warning level and the sudden change warning level to generate and output water quality warning information.

[0039] Based on this, the two warning results are merged. The fusion logic is used to balance gradual and abrupt risks. For example, when both are normal, the final warning is normal; when only one triggers the warning, the trigger level is output; when both trigger the warning simultaneously, it indicates overlapping risks, requiring an upgrade of the warning level to alert for compound risks. The output water quality warning information can include the comprehensive warning level, warning type (gradual / abrupt / compound), spatial distribution of the warning area, and possible warning lead time, providing management departments with information for taking emergency measures.

[0040] On the other hand, this describes the optional implementation methods for constructing and calculating the comprehensive characteristic index of lake currents, particularly the technical details of calculating hydraulic residence time using the Lagrange particle tracking method. In this embodiment, multi-dimensional hydrodynamic characteristics are introduced, solving the problem that traditional methods relying solely on flow velocity cannot accurately assess the risk of water retention. Accordingly, this embodiment can be carried out using the following steps:

[0041] Step 201: The lake current comprehensive characteristic index is obtained by weighted summation of the standardized values ​​of flow velocity, eddy current, and residence time.

[0042] Among them, the lake current comprehensive characteristic index is an indicator that comprehensively reflects the local flow field intensity, rotational retention characteristics, and overall exchange capacity. Its calculation formula can be expressed as:

[0043] LCI(x, y, t) = α _1 ×V _norm (x, y, t) + α _2 ×RT _norm (x, y, t) + α _3 ×ω _norm (x, y, t);

[0044] Where LCI(x, y, t) is the lake current comprehensive characteristic index at position (x, y) at time t, V _norm RT is the normalized value of the flow rate. _norm ω is the standardized value of the dwell time. _norm Let α be the vorticity normalized value. _1 α _2 α _3 These are the corresponding weight coefficients, and they satisfy α. _1 +α _2 +α _3 =1.

[0045] Regarding the weighting coefficient α _1 α _2 α _3 The value can be determined using the entropy weight method, and the specific process is as follows:

[0046] A historical period of monitoring data was selected as the sample set;

[0047] Normalize each indicator;

[0048] Calculate the information entropy E of the m-th indicator. _m =-(1 / ln(n))×∑p _mi ×ln(p _mi );

[0049] Among them, E _mLet p be the information entropy of the m-th indicator, n be the number of samples, ∑ represent the summation over all samples, and p _mi Let be the proportion of the i-th sample value of the m-th indicator to the total sum of that indicator, and ln be the natural logarithm.

[0050] Furthermore, the weight α of each indicator is calculated based on information entropy. _m =(1-E _m ) / ∑(1-E _k ); where α _m Let ∑(1-E) be the weight of the m-th indicator. _k ) represents the sum of the information entropy differences of all indicators, and k represents the traversal index of the indicator.

[0051] In some implementations, when there is a lack of sufficient historical data for entropy weight calculation, an equal-weight method can be used, i.e., setting α... _1 =α _2 =α _3 =1 / 3, or the weight can be determined based on expert experience using the Analytic Hierarchy Process (AHP).

[0052] Step 202: The velocity standardization value is determined based on the ratio of the velocity magnitude in the lake flow monitoring data to the multi-year average velocity; the eddy normalization value is determined based on the ratio of the absolute value of the flow field eddy calculated from the lake flow monitoring data to the multi-year average eddy; and the residence time standardization value is determined based on the ratio of the average residence time of the entire lake to the actual residence time of the grid.

[0053] Specifically, the calculation of the normalized velocity value is used to eliminate the influence of velocity dimensions and reflect the intensity of the current velocity relative to the background velocity. The calculation formula is as follows:

[0054] V _norm (x, y, t) = |V(x, y, t)| / V _ref ;

[0055] Where |V(x, y, t)| represents the flow velocity at that location and at that moment, typically expressed in meters per second; V _ref For reference flow velocity, the multi-year average flow velocity can be used. In other words, V(x, y, t) characterizes the magnitude and direction of the flow velocity in the lake area at spatial location (x, y) and time t.

[0056] The vorticity normalization value is used to characterize the rotational properties of a fluid. High vorticity usually corresponds to a vortex structure, which can easily lead to the local retention of pollutants. The vorticity ω of the flow field is calculated as: ω(x, y, t) = (Яv / Яx) - (Яu / Яy); where u and v are the eastward and northward components of the flow velocity, respectively, and Я represents the partial derivative operation. Next, normalization is performed:

[0057] ω _norm (x,y,t)=|ω(x,y,t)| / ω_ref ;

[0058] Where |ω(x, y, t)| is the absolute value of vorticity, ω _ref For reference vorticity, the average value of the absolute vorticity over many years can be taken.

[0059] It should be understood that vorticity is a physical quantity in fluid mechanics used to characterize the intensity of the rotational motion of a fluid particle. Its value is equal to twice the angular velocity of the fluid particle rotating about its own axis. When the local flow field in a lake area exhibits high vorticity characteristics, it means that there is a vortex structure or circulation system in that area. The vortex structure causes the water to form a relatively closed circulating flow path in and around the vortex center, making it difficult for water masses entering this area to migrate smoothly downstream, resulting in the retention and accumulation of pollutants within the vortex area.

[0060] Correspondingly, the presence of eddies also weakens the efficiency of material exchange between pollutants and surrounding clean water, reducing dilution and diffusion. Therefore, high eddy areas are high-risk areas for lake water quality deterioration. Incorporating eddy into the construction of a comprehensive lake flow characteristic index can identify hydrodynamically sensitive areas with the risk of pollutant retention.

[0061] Furthermore, the standardized residence time value was defined using inverse standardization, making it positively correlated with hydrodynamic intensity; that is, the shorter the residence time, the stronger the hydrodynamic force. _norm The larger the value, the more it can be calculated using the following formula:

[0062] RT _norm (j, k) = RT _ref / RT(j,k);

[0063] Where RT(j,k) is the actual dwell time of the (j,k)th grid, RT _ref For reference, the average stay time across the entire lake can be used.

[0064] Step 203: Deploy virtual tracer particles within the computational domain of the two-dimensional hydrodynamic model.

[0065] In other words, virtual tracer particles are deployed within the computational domain of a two-dimensional hydrodynamic model based on the lake region.

[0066] In this embodiment, the Lagrange particle tracking method is used to calculate the aforementioned hydraulic residence time. Accordingly, a two-dimensional hydrodynamic model of the lake area is used. This model can be constructed based on shallow water equations and calibrated using measured data. A large number of virtual tracer particles are uniformly deployed throughout the computational domain. For example, the initial positions of the particles can be determined at a density of no less than 100 particles per square kilometer. The particles are treated as point masses moving with the water flow, used to simulate the transport trajectory of pollutants.

[0067] Step 204: Based on lake current monitoring data, construct the Lagrange particle motion equation containing advection transport and random diffusion terms, and use the Lagrange particle motion equation to track the motion trajectory of virtual tracer particles.

[0068] Alternatively, based on lake current monitoring data, a Lagrange particle motion equation can be established, and the motion trajectory of virtual tracer particles can be simulated and tracked through this equation.

[0069] The motion of particles is influenced not only by the advection of the average flow velocity but also by the random effects of turbulent diffusion. Therefore, the Lagrange equations of motion for particles are constructed as follows:

[0070] dx _i / dt=u(x _i y _i ,t)+sqrt(2×D _h )×ξ _x (t);

[0071] dy _i / dt=v(x _i y _i ,t)+sqrt(2×D _h )×ξ _y (t);

[0072] Among them, dx _i / dt and dy _i / dt represent the velocities of the i-th particle in the east and north directions, respectively, u(x) _i y _i ,t) and v(x) _i y _i ,t) represents the velocity component at the particle's location, D _h ξ is the horizontal diffusion coefficient, sqrt(...) denotes the square root operation, and ξ... _x (t) and ξ _y (t) is a random number that follows a standard normal distribution and is used to simulate turbulent diffusion effects.

[0073] In other words, within each time step Δt, the particle position is updated as follows:

[0074] x _i (t+Δt)=x _i (t)+u(x _i ,y _i ,t)×Δt+sqrt(2×D _h ×Δt)×ξ _x ;

[0075] y _i (t+Δt)=y _i (t)+v(x_i ,y _i ,t)×Δt+sqrt(2×D _h ×Δt)×ξ _y ;

[0076] Where, x _i y _i These correspond to the spatial coordinates of the i-th particle in the X and Y directions, respectively, u(x _i y _i ,t),v(x _i y _i , t) respectively correspond to the particle at position (x) _i y _i The velocity components in the X and Y directions at time t, where Δt corresponds to the time step, D _h For the corresponding horizontal diffusion coefficient, sqrt(...) represents the square root operation, ξ _x and ξ _y This corresponds to Gaussian white noise with a mean of zero.

[0077] Optionally, a fourth-order Runge-Kutta method can be used for numerical integration. Compared to the first-order Euler method, the fourth-order Runge-Kutta method has higher computational accuracy and stability, and can more accurately simulate the motion path of particles in complex flow fields. By updating the integration at each time step, the position coordinates of each particle can be continuously recorded.

[0078] Step 205: Calculate the time taken for the virtual tracer particle to move from the release point to leave the boundary of the computational domain, and calculate the actual dwell time of each computational grid.

[0079] During the tracking process, the total time elapsed for each particle from the moment of release until its trajectory crosses the lake's exit boundary or reaches the preset maximum simulation duration is recorded and denoted as T. _i For each grid cell (j, k) within the computational domain, count all particles that pass through or are released from that grid, and calculate RT(j, k) = (1 / N). _jk )×∑T _i Where RT(j,k) is the actual dwell time of the grid, and N _jk Let ∑T be the total number of relevant particles. _i This represents the sum of the total residence times of the particles. Using this statistical method, a spatial distribution map of the hydraulic residence time across the entire lake area can be obtained, identifying the main flow areas with rapid water exchange and the stagnant water areas with slow water exchange.

[0080] Furthermore, hydraulic retention time, defined as the average time required for a lake to complete one water renewal cycle, characterizes the speed at which external water replaces the existing water in the lake.

[0081] In summary, incorporating the reciprocal form of the hydraulic residence time into the lake flow comprehensive characteristic index can make the index positively correlated with the water body's self-purification capacity. That is, the shorter the hydraulic residence time, the larger the standardized value, and the higher the lake flow comprehensive characteristic index, reflecting the physical law that better hydrodynamic conditions and lower pollution risk.

[0082] On the other hand, it describes the specific implementation process of time-varying time-delay identification and early warning prediction, solving the problem of early warning failure in different hydrological seasons due to the assumption of a constant time delay, and realizing the leap from lagging monitoring to advanced early warning. Furthermore,

[0083] Step 301: Standardize the time series of the lake flow comprehensive characteristic index and the time series of water quality monitoring data respectively.

[0084] In other words, standardization is performed to eliminate dimensional differences, specifically involving time series of lake flow comprehensive characteristic indices and time series of water quality monitoring data.

[0085] Since the lake current comprehensive characteristic index (LCI(t)) is a dimensionless comprehensive indicator, while water quality monitoring data C(t) (such as total phosphorus concentration) has specific physical units (such as milligrams per liter), their dimensions and orders of magnitude may differ. Therefore, Z-score standardization is required for both. Specifically, for the lake current comprehensive characteristic index series, the standardization formula is:

[0086] V _tilde (t)=(LCI(t)-μ _V ) / σ _V ;

[0087] Among them, V _tilde (t) represents the standardized lake current characteristic sequence, LCI(t) is the comprehensive lake current characteristic index at the original time t, and μ _V σ is the arithmetic mean of the sequence. _V denoted as the standard deviation of the sequence.

[0088] Similarly, for water quality monitoring data sequences, the standardization formula is:

[0089] C _tilde (t)=(C(t)-μ _C ) / σ _C ;

[0090] Among them, C _tilde (t) represents the standardized water quality sequence, C(t) represents the water quality concentration at the original time t, and μ _C σ is the arithmetic mean of the water quality series. _CThe standard deviation of the water quality series is represented by . This process converts both series into standard series with a mean of 0 and a standard deviation of 1, facilitating comparison of fluctuation trends.

[0091] Step 302: Within the time delay search range, calculate the time delay cross-correlation spectrum between the two standardized time series.

[0092] Cross-correlation spectrum (CCF) with time delay is used to quantify the correlation between two time series at different time shifts. Accordingly, a physically reasonable time delay search range needs to be set. Considering the material transport scale of lakes, especially large lakes like Poyang Lake, a maximum search time delay τ is typically set. _max The search period is 30 days. The search range is defined as τ belonging to the interval [-τ]. _max ,+τ _max For each integer time delay τ within this range, calculate the cross-correlation coefficient:

[0093] ρ _VC (τ)=(1 / (N-|τ|))×∑V _tilde (t)×C _tilde (t+τ);

[0094] Where, ρ _VC (τ) represents the cross-correlation coefficient with a time lag of τ, ranging from [-1, 1]; N is the total length of the time series; |τ| is the absolute value of the time lag; ∑ represents the summation operation on all samples within the effective overlapping time period. By iterating through all the set τ values, the spectrum reflecting the change in correlation with lag time can be obtained, namely the time-lag cross-correlation spectrum.

[0095] Step 303: Identify the time delay corresponding to the largest absolute value of the cross-correlation coefficient from the time delay cross-correlation spectrum, and use it as the optimal time delay.

[0096] The optimal time delay represents the time scale at which changes in lake currents have the greatest impact on water quality changes. Specifically, it is determined by finding the point of maximum absolute value in the calculated time delay cross-correlation spectrum.

[0097] τ _star =argmax _τ (|ρ _VC (τ)|);

[0098] Where, argmax _τ This represents the value of the independent variable τ that maximizes the objective function. It should be understood that if τ... _star A value greater than 0 indicates that changes in lake current precede changes in water quality, meaning that lake current is the cause and water quality is the effect; this time lag is the material transport time. If τ _star A value less than 0 may reflect the mechanism by which algal blooms alter water density and influence circulation. This invention primarily focuses on τ. _starWhen the value is greater than 0, water quality early warning based on lake flow can be achieved.

[0099] Step 304: Use the sliding window method to traverse the time series and identify the optimal time delay sequence that changes over time.

[0100] Considering the hydrological conditions of lakes, such as rapid flow velocity during the high-water season and slow flow velocity during the low-water season, exhibiting significant seasonal variations, the impact of lake currents on water quality is not constant over time. Therefore, this embodiment abandons the assumption of a globally fixed time delay and employs a sliding window technique to identify time-varying time delays. The specific operation is as follows:

[0101] A fixed-length time window W is defined, for example, W can be 90 days, and a sliding step size Δ, for example, Δ can be 1 day. For each starting position of the window, a subsequence of data within that window is extracted, and the above calculation process is repeated to identify the local optimal time delay within that specific time window. As the window slides across the entire time series, an optimal time delay sequence that varies with time is generated. This sequence can dynamically reflect the changes in material transport rates under different seasons and hydrological conditions.

[0102] Step 305: Construct a multiple regression relationship between the optimal time delay and the lake's concurrent water level and total inflow.

[0103] Furthermore, a quantitative regression model between time lag and hydrological driving factors was established to predict future time lags, i.e., the lead time for early warning, based on current hydrological conditions. Hydrodynamic principles indicate that the water level H(t) determines the lake's volume and flow field morphology, while the total inflow Q... _in The velocity (t) determines the flow rate, and both together determine the transport time of matter in the lake. Based on this, the following multiple linear regression model is constructed:

[0104] τ _star (t)=β _0 +β _1 ×H(t)+β _2 ×Q _in (t)+ε(t);

[0105] Where, τ _star H(t) represents the identified time-varying time-delay sequence, H(t) represents the lake water level at the corresponding time, in meters, and Q(t) represents the lake level at the corresponding time. _in (t) represents the total inflow into the lake at the corresponding time, in cubic meters per second, β _0 β is the constant intercept term. _1 β is the water level influence coefficient. _2 Let β be the flow rate influence coefficient, and ε(t) be the residual term. By fitting historical data using the least squares method, the model parameter β can be determined. _0 β _1 and β _2 .

[0106] In some alternative implementations, the regression model can be extended to a nonlinear form to take into account the nonlinear response relationship, for example by introducing logarithmic or interaction terms for water level and flow, or by using machine learning methods such as support vector regression (SVR) to construct more complex mapping relationships and improve fitting accuracy.

[0107] Step 306: Based on the real-time acquired current water level and current total inflow into the lake, the effective time lag is estimated using a multiple regression relationship. This effective time lag is then used as the lead time for issuing early warning information, thus enabling the forecasting function.

[0108] In actual operation, the system receives the current water level H from the hydrological telemetry station in real time. _now The inflow data Q provided by the hydrological station _now Substituting the two real-time variables into the regression model, the current estimated time lag is calculated, i.e.:

[0109] τ _eff =β _0 +β _1 ×H _now +β _2 ×Q _now ;

[0110] Where, τ _eff This is the current effective time delay, representing approximately τ if the lake current undergoes an abnormal change at the current moment. _eff After that, the water quality will deteriorate accordingly. Therefore, the system can... _eff It is released directly as the advance warning amount for early warning information. For example, if the calculated τ _eff If the timeframe is 5.5 days, the system can issue early warnings indicating that the water quality is at risk of exceeding standards in 5 to 6 days, allowing management departments sufficient time for emergency response, such as advance gate scheduling and ecological water replenishment.

[0111] In some embodiments, possible implementation methods for physical constraint kernel regression and dynamic threshold determination are provided, particularly how to introduce physical mechanism constraints to overcome the problem of pure data-driven models violating physical common sense in sparse data regions. This scheme also solves the technical problem that traditional fixed threshold early warning methods cannot adapt to the dynamic environmental capacity of lakes. Specifically, this method can be:

[0112] Step 401: Based on the optimal time delay identified in the time delay relationship, the time series of the lake flow comprehensive characteristic index is shifted and aligned with the time series of water quality monitoring data to construct an input-output paired dataset; the Nadaraya-Watson kernel regression method is used to perform nonparametric fitting of the paired dataset using the Gaussian kernel function to obtain the estimation function as the lake flow-water quality response model.

[0113] In this embodiment, a sample set needs to be constructed for model training. Because the impact of lake currents on water quality has a lag, directly pairing data from the same moment cannot reflect the true causal relationship. Therefore, the optimal time delay τ is utilized... _star The time series of the lake current composite characteristic index LCI(t) is shifted backward by τ. _star This aligns the time frame of the water quality monitoring data C(t).

[0114] Furthermore, the constructed paired dataset D can be represented as a set {(LCI(t-τ)} _star ), C(t))}, where t takes over all valid observation times. In some implementations, when more physical variables need to be considered, the paired dataset can also be extended to include the normalized velocity value V. _norm and standardized value of dwell time RT _norm That is, D={(LCI(t-τ)} _star V _norm (t-τ _star ), RT _norm (t-τ _star ), C(t))}.

[0115] For this dataset, the Nadaraya-Watson kernel regression method can be used for nonparametric estimation. Kernel regression does not presuppose the specific form of the function; instead, its shape is determined by the data, allowing it to fit the complex nonlinear relationships of the lake ecosystem. Its core estimation formula is:

[0116] g _hat (v)=∑(K _h (LCI _i -v)×C _i ) / ∑K _h (LCI _i -v);

[0117] Among them, g _hat (v) represents the estimated water quality response when the lake current characteristic value is v, LCI _i Let C be the comprehensive characteristic index of the lake currents for the i-th sample. _i Let K be the water concentration of the i-th sample, ∑ represent the summation over all samples, and K is the water concentration of the ith sample. _h Let h be the kernel function with bandwidth h.

[0118] Alternatively, a Gaussian kernel function can be used, whose expression can be described as:

[0119] K _h (z)=(1 / (sqrt(2π)×h))×exp(-z 2 / (2×h 2 ));

[0120] Among them, K _h (z) represents the value of the kernel function at z, π is the mathematical constant pi, h is the bandwidth parameter, and exp is the exponential function. The bandwidth parameter h controls the smoothness of the model: too small a value of h will lead to overfitting, while too large a value of h will lead to underfitting. It should be understood that leave-one-out cross-validation (LOOCV) can also be used, that is, selecting the value of h that minimizes the cross-validation error to determine the optimal bandwidth.

[0121] Step 402: Invoke the optimization objective function that includes data fitting terms and physical constraint penalty terms;

[0122] Define physical constraint penalty terms, which include at least: a transport dilution constraint based on the transport dilution principle, used to penalize the case where the partial derivative of the response model with respect to the normalized value of velocity is positive; and a retention accumulation constraint based on the retention accumulation principle, used to penalize the case where the partial derivative of the response model with respect to the normalized value of residence time is positive.

[0123] By minimizing the objective function, a lake current-water quality response model that satisfies physical constraints is obtained.

[0124] Traditional kernel regression relies on data, but under extreme hydrological conditions, such as catastrophic floods or extremely low water levels, the scarcity of historical observation data may lead the model to predict results that violate physical laws, such as increased concentrations at higher flow velocities. Therefore, this invention proposes a physics-guided kernel regression method. Accordingly, the following optimization objective function is constructed:

[0125] J(g)=L _data (g)+λ _1 ×P _transport (g)+λ _2 ×P _retention (g)+λ _3 ×P _bound (g);

[0126] Where J(g) is the overall objective function, L _data (g) is a loss term that measures the goodness of fit to the data, such as weighted least squares error, P. _transport P _retention P _boundλ represents the penalties for transport dilution, retention accumulation, and boundary constraints, respectively. _1 , λ _2 , λ _3 This represents the corresponding penalty coefficient.

[0127] Specifically, the transport-dilution constraint is based on the fundamental principle of mass transport: the higher the flow velocity, the stronger the dilution and diffusion effect, and the pollutant concentration should decrease. Therefore, the response function g is normalized to the flow velocity V. _norm The partial derivatives should be less than or equal to zero. If the model predicts a positive partial derivative, violating physical laws, a penalty is applied using the following integral formula:

[0128] P _transport (g)=∫(max(0,Яg / ЯV _norm ) 2 )dV _norm ;

[0129] Where ∫ represents integration, Яg / ЯV _norm The max(0,...) function is used to accumulate the penalty only when the partial derivative of the response function with respect to the normalized value of the velocity.

[0130] Similarly, the retention accumulation constraint is based on the principle of water exchange, with the standardized value of residence time RT. _norm The larger the value, the shorter the actual residence time, the faster the exchange, the less likely pollutants are to accumulate, and the lower the concentration should be. Therefore, the response function affects RT. _norm The partial derivatives should also be less than or equal to zero. Its penalty term is defined as:

[0131] P _retention (g)=∫(max(0,Яg / ЯRT _norm ) 2 )dRT _norm ;

[0132] Among them, Яg / ЯRT _norm This is the partial derivative of the response function with respect to the standardized value of the dwell time.

[0133] In this embodiment, a boundary bounded constraint P is also set. _bound This constraint penalizes any value that exceeds physically reasonable limits, such as negative concentrations or predicted values ​​that exceed historical high concentrations by a certain percentage. Furthermore, by minimizing the objective function J(g) using numerical optimization algorithms such as Sequential Quadratic Programming (SQP), a response model that both fits the observed data and obeys physical laws can be obtained.

[0134] Optionally, the obtained discrete function values ​​can be post-processed by order-preserving regression to eliminate any minor fluctuations that may remain in the numerical calculation.

[0135] Step 403: Calculate the local data density of the lake current comprehensive characteristic index at different value points; dynamically adjust the adaptive penalty coefficient according to the local data density so that the value of the adaptive penalty coefficient increases in areas with low local data density, thereby enhancing the dominant role of physical constraints.

[0136] Accordingly, an adaptive penalty mechanism is introduced. In data-intensive regions, the observed data should be trusted, and the weight of physical constraints should be reduced; in data-sparse regions, such as extreme conditions that have never occurred in history, the data reliability is low, and deductions should be made based on physical laws, thus increasing the weight of physical constraints.

[0137] Accordingly, the local data density ρ(v) at the value v of the independent variable is calculated. This can be achieved using kernel density estimation or by counting the number of samples within a radius of v and bandwidth h. Next, the adaptive penalty coefficient λ(v) is calculated using the following formula:

[0138] λ _m (v)=λ _m_0 ×exp(-ρ(v) / ρ _ref );

[0139] Where, λ _m (v) is the penalty coefficient for the m-th constraint at v, λ _m_0 The basic penalty coefficient can be determined by generalized cross-validation (GCV), ρ _ref This serves as a reference density, such as the median of the density of all samples.

[0140] As the formula shows, when ρ(v) approaches 0, i.e., when the data is sparse, the exponential term approaches 1, the penalty coefficient is large, and the model is subject to strong physical constraints, forcing its prediction results to conform to the dilution and retention laws; when ρ(v) is large, i.e., when the data is dense, the exponential term approaches 0, the penalty coefficient is extremely small, and the model is mainly composed of L _data The dominant element can capture the unique nonlinear details in the data.

[0141] In other embodiments, the local data density of the lake current composite characteristic index at different value points is calculated;

[0142] Based on this, the adaptive penalty coefficient is dynamically adjusted so that in areas where the local data density is lower than the density threshold, the value of the adaptive penalty coefficient increases, thereby enhancing the dominant role of physical constraints.

[0143] Step 404: Divide the historical water quality monitoring data into several data groups according to the magnitude of the corresponding lake flow comprehensive characteristic index;

[0144] Calculate the conditional quantiles of water quality concentrations for each data group at the preset risk level;

[0145] Based on this (conditional quantiles of each data group), the quantile regression method is used to fit the functional relationship between the dynamic warning threshold and the comprehensive characteristic index of the lake flow, which serves as the basis for determining the dynamic warning threshold.

[0146] Traditional fixed thresholds, such as the constant Class III water standard, ignore the dynamic changes in lake environmental capacity with hydrological conditions. In this embodiment, when the lake current dynamics are strong (high LCI), the water body's self-purification capacity is enhanced, and it can tolerate slightly higher pollutant concentrations without algal blooms; conversely, strict control is required.

[0147] Accordingly, the range of LCI values ​​is determined and divided into K intervals, for example, K=10. Historical monitoring data is then compared with (LCI...) _t C _t The data is projected into the corresponding intervals, forming K data sets. Next, for each data set, a specific quantile of its water quality concentration is calculated. For example, if the risk level for the warning is set to α=0.1, that is, allowing a 10% risk of exceeding the standard, then the 90th quantile of the data set is calculated, resulting in a series of discrete point pairs.

[0148] Furthermore, quantile regression is used to fit the above point pairs to obtain a continuous threshold function. Typically, the threshold function C is assumed to be... _th (LCI) follows a quadratic polynomial form:

[0149] C _th (LCI)=β _0 +β _1 ×LCI+β _2 ×LCI^2;

[0150] Where, β _0 β _1 β _2 These are the regression coefficients. The above coefficients can be estimated by minimizing the weighted absolute deviation objective function. The resulting curve C... _th (LCI) is the dynamic early warning threshold.

[0151] In practical applications, multiple curves can be fitted for different risk levels, such as α=0.3, 0.1, and 0.05, corresponding to blue, yellow, and red warning thresholds, respectively. When the real-time monitored water quality concentration C... _now Exceeding the corresponding LCI _now C below _th When the value is reached, the corresponding level of warning is triggered. This method makes the warning standards more scientific and accurate, reducing missed reports due to decreased environmental capacity during the dry season and false reports due to increased background values ​​during the wet season.

[0152] In other embodiments, an alternative implementation of spatially propagated enhanced mutation early warning is described, which upgrades traditional single-point time-series early warning to spatially linked early warning, improving the ability and timeliness of capturing sudden aquatic environmental events such as algal blooms.

[0153] From a dynamic perspective, stable ecosystems possess strong negative feedback regulation capabilities. When disturbed by external forces and deviating from equilibrium, the internal negative feedback mechanism can quickly pull the state variables back to their equilibrium position. However, as the system gradually approaches a critical point, the strength of the negative feedback regulation gradually weakens, the system's response to disturbances becomes sluggish, and its resilience approaches zero. Near the critical point, even small random disturbances can cause significant fluctuations in the system's state, and these fluctuations are difficult to quickly smooth out.

[0154] Critical slowing down exhibits detectable early warning signals in the statistical characteristics of time series, primarily manifested in the following three aspects: First, increased variance. Due to decreased resilience, the system's ability to suppress random disturbances weakens, leading to increased fluctuations in state variables and a growing trend in the variance of the time series. Second, enhanced autocorrelation. Slower recovery speed implies a stronger dependence of the system's current state on past states, resulting in a stronger correlation between observations at adjacent times, with the first-order autocorrelation coefficient approaching one. Third, changes in skewness. As the system approaches an alternative steady state in a certain direction, the probability distribution of state variables gradually shifts in that direction, exhibiting a systematic change in distribution skewness.

[0155] Lake ecosystems are typical complex systems with multiple stable states, exhibiting two alternative stable states: clear water and turbid water. When a lake transitions from a clear water state to a turbid water state, the water quality time series often displays the aforementioned critical slowdown characteristics before reaching the critical point. Therefore, monitoring the changing trends of variance, autocorrelation coefficient, and skewness of water quality time series to capture early signals of lake ecosystems approaching critical transitions can help achieve early warning before abrupt changes in water quality occur.

[0156] This embodiment specifically includes:

[0157] Step 501, the critical slowing indicators include the variance growth rate, the autocorrelation enhancement rate, and the skewness change rate;

[0158] The critical slowdown index is calculated by: extracting time series segments of water quality monitoring data using a sliding window, calculating the variance, lag 1 autocorrelation coefficient and skewness within each window, and calculating the growth trend of statistical characteristics over time based on the Mann-Kendall trend test method, thereby obtaining the variance growth rate, autocorrelation enhancement rate and skewness change rate.

[0159] Furthermore, a sliding window processing method needs to be applied to the water quality time series C(t) obtained for each monitoring point. The window length is set to W. _s For example, 60 days, step size Δ _s For example, 1 day. For each window position, calculate the following three key statistical characteristics of the data within the window:

[0160] Firstly, variance, which reflects the fluctuation range of the data, is calculated using the following formula:

[0161] Var=(1 / (W _s -1))×∑((C(t)-μ) 2 );

[0162] Where Var is the variance, μ is the mean of the data within the window, and ∑ represents the summation of all data within the window.

[0163] Secondly, the lag-1 autocorrelation coefficient reflects the strength of the system's memory, and its calculation formula is as follows:

[0164] AR1=∑((C(t)-μ)×(C(t+1)-μ)) / ∑((C(t)-μ) 2 );

[0165] Where AR1 is the autocorrelation coefficient, the numerator represents the sum of the products of the deviations between adjacent time points, and the denominator is the variance term.

[0166] Third, skewness, which reflects the asymmetry of the data distribution. Its calculation formula is:

[0167] Skew=(1 / W _s )×∑(((C(t)-μ) / σ) 3 );

[0168] Where Skew is the skewness and σ is the standard deviation of the data within the window.

[0169] When an ecosystem approaches a critical transition point, due to decreased resilience, three indicators will show an increasing trend: increased variance, enhanced autocorrelation, and increased skewness. The Mann-Kendall (MK) trend test can be used to quantify this increasing trend. For any indicator's time series X={x _1 x _2 , ..., x _n}, where n represents the length of the time series, and its MK statistic S _MK The calculation is as follows:

[0170] S _MK =∑ _i (∑ _j (sign(x _j -x _i )));

[0171] In this loop, the outer loop i ranges from 1 to n-1, the inner loop j ranges from i+1 to n, and sign(...) is the sign function, which is used when x... _j >x _i x is 1 if x is less than x, -1 if x is equal to x. _j x _i This represents the index values ​​at different times in the sequence. Based on S _MK Standardized growth rate indicators can be calculated, such as variance growth rate (VGR), autocorrelation enhancement rate (AGR), and skewness change rate (SGR). A larger positive value indicates a stronger growth trend and that the system is closer to its critical point.

[0172] Step 502: Calculate local early warning indicators for each monitoring point in the lake area; construct a spatial propagation weight matrix between monitoring points based on the lake flow direction relationship; for any target monitoring point, identify its upstream monitoring point according to the spatial propagation weight matrix, and weight and superimpose the local early warning indicators of the upstream monitoring point onto the local early warning indicators of the target monitoring point to obtain an early warning comprehensive indicator enhanced by spatial propagation.

[0173] The direction of lake flow is determined based on the velocity vector in the lake flow monitoring data.

[0174] In this embodiment, a spatial propagation mechanism is introduced to use upstream information to provide early warnings downstream. Accordingly, for each monitoring point i, the three growth rate indicators are weighted and synthesized to obtain the local early warning indicator EWS. _local (i):

[0175] EWS _local (i)=w _1 ×VGR(i)+w _2 ×AGR(i)+w _3 ×SGR(i);

[0176] Among them, w _1 w _2 w _3 The weighting coefficients can be determined through historical backtracking analysis or by using equal weights (1 / 3 each).

[0177] Next, we compute the spatial propagation-enhanced early warning composite index EWS. _spatial (i). This indicator not only includes local signals but also incorporates signals from upstream monitoring point j (EWS). _local (j). Its synthesis formula is:

[0178] EWS _spatial (i)=α _local ×EWS _local (i)+α _upstream ×∑_j (W _ji_norm ×EWS _local (j)×T _delay (j, i));

[0179] Where, α _local and α _upstream The weights are assigned to the local and upstream signals respectively, for example, 0.6 and 0.4, W. _ji_norm T represents the normalized spatial propagation weights from upstream point j to target point i. _delay (j, i) is the attenuation factor considering propagation time, ∑ _j This represents the summation over all upstream points that affect point i. Through this summation, even if the local index of downstream point i has not yet increased, as long as its upstream point j issues a warning signal, the EWS of point i will be affected. _spatial The value will also rise in advance, achieving an early warning that is not yet known.

[0180] For the target monitoring point i, the original weights W_ji of all upstream points j that affect it are normalized to obtain W. _ji_norm =W _ji / ∑ _k W _ki , where ∑ _k W _ki It is the sum of the original weights of all upstream points to point i.

[0181] Step 503: The spatial propagation weight is determined based on the angle between the direction vector from the upstream monitoring point to the target monitoring point and the direction of the lake flow at the upstream monitoring point, as well as the distance between the two monitoring points. The spatial propagation weight is positively correlated with the cosine of the angle and decreases exponentially with the increase of distance.

[0182] The calculation of the weights requires determining whether point j is truly located within the effective influence zone upstream of point i, which depends on the distance d between the two points. _ji And the flow direction relationship. The specific calculation formula is:

[0183] W _ji =exp(-d _ji / L _c )×max(0,cos(θ _ji ));

[0184] Among them, W _ji d represents the original propagation weight from point j to point i; _ji L is the Euclidean distance between two points; _c The characteristic attenuation length, for example, can be 1 / 3 of the average scale of a lake, used to control the influence range of distance; the greater the distance, the smaller the weight. θ _jiLet be the angle between the vector pointing from point j to point i and the actual lake current velocity vector at point j.

[0185] When the included angle θ _ji When the angle is less than 90 degrees, cos(θ) _ji A positive angle (θ) indicates that point i is downstream of point j, allowing matter and signals to propagate downstream. Therefore, the weight is positive and increases as the angle decreases, or in other words, the weight is greatest when directly downstream. When the angle θ is positive... _ji When the angle is greater than or equal to 90 degrees, cos(θ) _ji If the value is negative or zero, it means that point i is located upstream or laterally in an irrelevant region of point j, and the material cannot flow downstream to reach it. Therefore, the weight is truncated to 0 by the max function.

[0186] Considering that signal propagation takes time, the timeliness during propagation is determined by the time decay factor T. _delay To reflect this, it can be calculated using the following method:

[0187] T _delay (j, i) = exp(-d _ji / (V _mean ×τ _c ));

[0188] Among them, V _mean Let d be the average flow velocity between two points. _ji / V _mean The time required for material transport, τ, was estimated. _c This is the characteristic time constant. This factor indicates that if the propagation time is too long, for example, if there is a large time beyond the characteristic time, the reference value of the upstream signal will decrease.

[0189] Step 504: Calculate the spatial synchronicity index of multiple monitoring points in the lake area that simultaneously exhibit critical slowdown characteristics;

[0190] A correction coefficient is generated based on the spatial synchronicity index. This correction coefficient is then used to amplify and correct the comprehensive early warning index, so that when the spatial synchronicity index increases, the level of sudden change warning is correspondingly increased.

[0191] When lakes face the risk of systemic ecological collapse, they exhibit synchronous critical slowdown characteristics over a large area. Optionally, a spatial synchronicity index Sync=(1 / (M×(M-1)))×∑ _i (∑ _j (I _i ×I _j ×exp(-d _ij / L _sync )));

[0192] Where M is the total number of monitoring points; I _i and I _jThe indicator function takes the value 1 when the local EWS of point i or j exceeds a certain threshold, such as the historical 75th percentile, and 0 otherwise; d _ij L represents distance; _sync This is the length of the synchronicity feature. This indicator quantifies the degree of spatial clustering and synchronicity of high-risk points.

[0193] Based on this indicator, a correction coefficient γ is constructed. _sync =1+β _sync ×max(0, Sync-Sync) _th ); where β _sync This is the amplification factor, for example, 2.0, Sync _th This serves as the baseline threshold for synchronicity. The final comprehensive early warning index, EWS, is... _final (i)=γ _sync ×EWS _spatial (i).

[0194] Based on this, when multiple points across the entire lake are detected to be synchronously abnormal, gamma is used. _sync By amplifying the warning indicator values ​​of all points, it is easier to trigger higher-level warnings and prevent signal fluctuations at a single site from being misjudged as local noise.

[0195] Step 505: Calculate the baseline mean and standard deviation of the early warning comprehensive index using historical data from periods without sudden changes; calculate the standardized anomaly of the early warning comprehensive index at the current moment relative to the baseline mean; determine the sudden change warning level based on the range of standardized anomalies, where a higher standardized anomaly corresponds to a higher sudden change warning level. This is used to map continuously changing EWS values ​​to discrete warning levels.

[0196] Accordingly, a period of historically stable water quality without abrupt changes was selected as the baseline period, and the EWS during this period was calculated. _final The mean μ of the indicator _EWS and standard deviation σ _EWS For the current indicator value EWS obtained through real-time calculation _now Calculate its standardized outlier, i.e.:

[0197] Z _EWS =(EWS _now -μ _EWS ) / σ _EWS ;

[0198] In the formula, Z _EWS It reflects the degree to which the current state deviates from the normal baseline, expressed in standard deviations. According to Z... _EWS The magnitude of the mutation warning level is determined according to the preset range.

[0199] For example, if Z _EWSA value less than or equal to 1.5 is considered normal; if Z... _EWS Values ​​greater than 1.5 and less than or equal to 2.0 are considered "attention-worthy"; if Z... _EWS A value greater than 2.0 and less than or equal to 2.5 is considered a warning sign; if Z... _EWS A value greater than 2.5 is considered a potential mutation warning. The grading standard based on statistical significance has good universality and can adapt to the dimensional differences of different lakes and indicators.

[0200] According to one aspect of this application, suppose that monitoring stations A, B, and C are sequentially deployed along the main current of a lake area, with station A located upstream, station C downstream, and station B midstream. The distances between the three stations are as follows: 5 kilometers from station A to station B, and 3 kilometers from station B to station C. The current average flow velocity of the lake is 0.05 m / s, and the flow direction is from station A to station C.

[0201] At a certain moment, due to changes in the water quality of the upstream water, the local early warning indicator (EWS) at station A is triggered. _local (A) was the first to rise to an abnormal level, with a calculated value of 0.8, exceeding the warning threshold of 0.5. Meanwhile, the local indicators for sites B and C remained within the normal range. _local (B) is 0.3, EWS _local (C) is 0.2.

[0202] If the traditional single-point early warning method is used, only station A will trigger an early warning, while stations B and C will not issue any warning information because their local indicators have not exceeded the threshold. However, from a physical mechanism perspective, the abnormal signal from station A will inevitably propagate downstream with the lake flow, and stations B and C will face the risk of water quality deterioration in the near future.

[0203] The spatial propagation enhancement method of this invention is calculated as follows. Accordingly, the spatial propagation weights are calculated. Since station A is located directly upstream of station B, the angle between the line connecting the two points and the lake flow direction is close to 0 degrees, and the cosine value is close to 1. Assuming the characteristic attenuation length is 10 kilometers, the propagation weight of station A to station B is W(A, B) = exp(-5 / 10) × cos(0°) ≈ 0.61. Similarly, the propagation weight of station A to station C is approximately 0.45, and the propagation weight of station B to station C is approximately 0.74.

[0204] Furthermore, we calculate the early warning indicator after spatial propagation enhancement. Assuming a local weight of 0.6 and an upstream weight of 0.4, the spatial enhancement indicator for site B is: EWS _spatial (B) = 0.6 × 0.3 + 0.4 × 0.61 × 0.8 × T _delay Assume a time decay factor T. _delayIf the value is 0.9, the calculated result is approximately 0.36. The spatial enhancement index for station C, taking into account the contributions from both station A and station B, is approximately 0.34.

[0205] After spatial propagation enhancement, although the index values ​​of stations B and C are still below the warning threshold, their values ​​are significantly higher than those considering only local signals. If the anomaly at station A continues or worsens, the spatial enhancement indexes of stations B and C will further increase and may trigger a warning, achieving the effect of early warning.

[0206] The estimation of the warning lead time can be calculated based on the distance between stations and the lake flow velocity. The time required for the signal from station A to travel to station B is approximately 5000m / 0.05m / s = 100000s, or about 1.2 days. The time required to travel to station C is approximately 8000m / 0.05m / s = 160000s, or about 1.9 days. That is, through the spatial propagation enhancement mechanism, station B can receive the warning information about 1 day in advance, and station C can receive the warning information about 2 days in advance, thus providing downstream management departments with an emergency response window.

[0207] The above examples demonstrate that the space propagation-enhanced early warning mechanism can overcome the lag problem of traditional single-point early warning methods. By integrating the leading signals from upstream monitoring stations, downstream areas can receive risk warnings before pollutants actually arrive, thus improving the timeliness and practical value of the early warning system.

[0208] In other embodiments, exemplary schemes for integrating decision-making and early warning information output are provided, specifically including:

[0209] Step 601: When both the normal warning level and the sudden change warning level are normal, the final warning level is determined to be normal.

[0210] When only one of the regular warning level and the sudden warning level triggers a warning, the final warning level is determined to be the level that triggered the warning.

[0211] When both the regular warning level and the sudden change warning level are triggered simultaneously, it is determined to be a compound risk. The higher of the two warning levels is selected, and the final warning level is raised by one level based on the higher warning level.

[0212] In this embodiment, the warning level is quantified into numerical values, for example, Normal = 0, Blue / Attention = 1, Yellow / Alert = 2, and Red / Sudden Change Warning = 3. The fusion decision formula can be expressed as:

[0213] Integrated Early Warning Level _fusion =min(3, max(Level) _routine Level _sudden )+I_both );

[0214] Among them, Level _routine This is the standard warning level value, Level _sudden This represents the numerical value for the mutation warning level. I _both For dual-trigger identifier: when Level _routine >0 and Level _sudden When >0, I _both =1, otherwise I _both =0. min(3,...) ensures that the final level does not exceed the highest level (red).

[0215] For example, if the regular warning is yellow (2) and the mutation warning is normal (0), the final warning is yellow, indicating a gradual risk. If the regular warning is blue (1) and the mutation warning is alert (2), both are triggered, indicating a compound risk. The larger value (2) is taken, and the warning level is increased by one level (+1), ultimately resulting in a red warning (3). This upgrade mechanism reflects the high importance attached to compound risks, because when water quality has exceeded the standard (regular trigger) and ecosystem stability is being lost (mutation trigger), the probability of catastrophic consequences is high.

[0216] Step 602: Generate and output water quality early warning information.

[0217] Furthermore, the output information includes not only the fused warning level but also multi-dimensional decision support information. Specifically, water quality warning information may include:

[0218] Comprehensive warning level, such as red warning.

[0219] Risk type diagnosis indicates whether it is a gradual water quality deterioration, a potential sudden risk, or a complex high risk.

[0220] Early warning lead time is determined by the effective time lag τ. _eff Provide, for example, an expected appearance in 5 days.

[0221] Spatial distribution hotspot map, a lake area risk map drawn based on the warning level of each grid point, shows high-risk areas and their expansion trend with the flow field.

[0222] A snapshot of key indicators, including current LCI values, real-time concentrations, dynamic thresholds, and EWS values, is provided for expert review.

[0223] The above information can be sent to environmental protection departments or water management agencies through visual dashboards, mobile push notifications, or data interfaces.

[0224] According to another aspect of this application, after introducing a physical constraint penalty term, the obtained discrete response function value may still exhibit slight non-monotonic fluctuations during the numerical solution process, especially when using approximate optimization algorithms. To satisfy the monotonically decreasing physical constraints of transport dilution and retention accumulation, this embodiment may employ an order-preserving regression algorithm to correct the preliminary results.

[0225] Assume that the response function at discrete grid points v _1 <v _2 <... <v _J The estimated value on is g _hat ={g _1 g _2 , ..., g _J}, where j is the total number of discrete grid points. Based on physical constraints, the desired corrected value g is... _mono Satisfy g _mono_1 ≥g _mono_2 ≥...≥g _mono_J That is, it is monotonically decreasing.

[0226] Optionally, the Pooled Adjacent Violation Algorithm (PAVA) can be used for correction. The steps are as follows:

[0227] Check adjacent pairs in the sequence (g _j g _j+1 If g is satisfied. _j ≥g _j+1 If a violation of monotonicity is found, i.e., g remains unchanged. _j <g _j+1 If two points or intervals containing them are to be merged, their weighted average will be used instead, as shown in the following formula:

[0228] The corrected result g after interval merging _new =(w _j ×g _j +w _j+1 ×g _j+1 ) / (w _j +w _j+1 );

[0229] Among them, w _j and w _j+1 This represents the weight of the corresponding point, typically set to 1 or the local sample size. The checking and merging process is repeated until the entire sequence satisfies the monotonically decreasing condition. The corrected response function will conform to physical laws, eliminating numerical noise and improving the stability and interpretability of the dynamic threshold.

[0230] According to one aspect of this application, the method of the present invention can be encoded as a computer program and run on a computer system. The computer system includes:

[0231] A processor is used to execute program instructions stored in memory. The processor can be a central processing unit (CPU), a graphics processing unit (GPU), or a dedicated digital signal processor (DSP). In this invention, the processor is specifically used to perform high-computational tasks such as LCI exponent calculation, Lagrange particle tracking integral, kernel regression optimization solution, and time-delay cross-correlation analysis.

[0232] The memory is used to store program code, lake current / water quality monitoring data, pre-built response models, and historical databases. Memory includes high-speed random access memory (RAM) and non-volatile memory such as hard drives and SSDs.

[0233] The input / output interface is used to receive multi-source monitoring data from monitoring stations, satellites, or external databases, and to send early warning information to display terminals or remote servers.

[0234] The communication module is used to enable real-time data transmission via wired or wireless networks, such as 4G / 5G, long-range LoRa radio, and satellite communication.

[0235] In other embodiments, specific numerical examples are provided to illustrate the calculation process of the lake flow composite characteristic index. Assume that the raw data at a certain monitoring point at a certain moment is: flow velocity magnitude |V| = 0.15 meters per second, and the multi-year average flow velocity V _ref =0.10 m / s; absolute vorticity |ω| = 0.002 s / s, multi-year average vorticity ω _ref =0.0025 seconds; the actual residence time of this grid is RT=25 days, and the average residence time of the entire lake is RT. _ref =30 days. The standardized values ​​are calculated as follows:

[0236] V _norm =0.15 / 0.10=1.50; ω _norm =0.002 / 0.0025=0.80; RT _norm =30 / 25=1.20.

[0237] Assuming an equal-weighted method is used (α1=α2=α3=1 / 3), the comprehensive characteristic index of the lake flow at this monitoring point is:

[0238] LCI=(1 / 3)×1.50+(1 / 3)×1.20+(1 / 3)×0.80=1.17.

[0239] A value greater than 1 indicates that the current hydrodynamic conditions at this monitoring point are stronger than the multi-year average, water exchange is more active, and the ability to dilute pollutants is relatively strong.

[0240] In some embodiments, during the data preprocessing stage, this embodiment also includes strategies for handling boundary conditions and anomalies, as follows:

[0241] Optionally, when the duration of missing data at a monitoring point is less than a preset threshold of the sliding window length, a time series interpolation method, such as linear interpolation or spline interpolation, is used to fill the missing data; when the duration of missing data exceeds the preset threshold, the calculation result of that window is marked as unreliable and its weight is reduced during fusion decision-making.

[0242] Optionally, outliers in the raw data can be identified using the 3σ criterion or the interquartile range (IQR) method. For data points identified as outliers, check for sensor malfunctions or recording errors; if confirmed as genuine extreme events, the data is retained and labeled in subsequent analyses.

[0243] Optionally, when the real-time monitored comprehensive characteristic index of the lake flow or the water level exceeds the coverage of historical training data, the system automatically triggers an extrapolation warning, indicating that the model prediction may have entered a region with high uncertainty. At this time, the adaptive penalty coefficient in the physical constraint kernel regression will automatically increase to ensure that the prediction results do not violate basic physical laws.

[0244] Accordingly, this scheme adopts a physics-guided kernel regression method. This method incorporates physical constraint penalty terms such as transport dilution and retention accumulation into the objective function of nonparametric estimation, and adaptively adjusts the penalty weights through local data density. This dual-driven mechanism of mechanism and data forces the model to follow the basic physical laws of material transport even under extreme conditions where historical data is scarce. This eliminates logical errors such as anomalous concentration increases at high flow rates, and improves the model's generalization ability and physical interpretability.

[0245] Furthermore, this scheme also constructs a spatial propagation-enhanced early warning mechanism based on flow field relationships. It calculates the flow direction angle and distance attenuation factor to construct a spatial propagation weight matrix, weighting and superimposing the critical slowing signals (variance, autocorrelation, etc.) from upstream monitoring points to the downstream, and then correcting them using spatial synchronicity indicators. This allows the downstream region to utilize upstream precursor signals for advanced sensing, solving the lag problem of single-point early warning and improving emergency response time for sudden events such as algal blooms.

[0246] In addition, this scheme also establishes a dynamic threshold system based on the comprehensive characteristic index of lake flow, realizing the dynamic adjustment of early warning standards with hydrodynamic conditions.

[0247] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately. Technical contents not described in detail in this invention can be implemented using conventional methods in the art.

Claims

1. A method for forecasting and early warning of water quality changes in lake areas based on the complex relationship between lake currents and water quality, characterized in that, include: Acquire multi-source monitoring data for the lake area, which should include at least lake flow monitoring data and water quality monitoring data; Based on lake flow monitoring data, a comprehensive characteristic index of lake flow reflecting the hydrodynamic conditions of the lake area is calculated; Identify the time lag relationship between the comprehensive characteristic index of lake flow and water quality monitoring data, construct a lake flow-water quality response model, determine the dynamic early warning threshold that changes with the comprehensive characteristic index of lake flow based on the lake flow-water quality response model, and determine the regular early warning level based on the dynamic early warning threshold; A sliding window analysis was performed on the time series of water quality monitoring data to calculate the critical slowing index. An early warning comprehensive index was synthesized based on the critical slowing index, and the level of sudden change warning was determined based on the early warning comprehensive index. The system integrates regular and sudden warning levels to generate and output water quality warning information. The comprehensive characteristic index of the lake current is obtained by weighted summation of the standardized values ​​of velocity, eddy current, and residence time. The standardized value of velocity is determined by the ratio of the velocity magnitude in the lake current monitoring data to the multi-year average velocity; the standardized value of eddy current is determined by the ratio of the absolute value of the eddy current in the flow field calculated from the lake current monitoring data to the multi-year average eddy current; and the standardized value of residence time is determined by the ratio of the average residence time of the entire lake to the actual residence time of the grid. The process of identifying the time lag relationship between the comprehensive characteristic index of lake currents and water quality monitoring data includes: standardizing the time series of the comprehensive characteristic index of lake currents and the time series of water quality monitoring data respectively; calculating the time lag cross-correlation spectrum between the two standardized time series within the time lag search range; identifying the time lag corresponding to the largest absolute value of the cross-correlation coefficient from the time lag cross-correlation spectrum as the optimal time lag; and using the sliding window method to traverse the time series and identify the optimal time lag sequence that changes over time. The construction of the lake current-water quality response model includes: based on the optimal time delay identified in the time delay relationship, shifting the time series of the lake current comprehensive characteristic index and aligning it with the time series of water quality monitoring data to construct an input-output paired dataset; using the Nadaraya-Watson kernel regression method, and using the Gaussian kernel function to perform nonparametric fitting on the paired dataset to obtain the estimation function of the lake current-water quality response model. The synthesis of an early warning comprehensive index based on the critical slowing index includes: calculating local warning indices for each monitoring point in the lake area; constructing a spatial propagation weight matrix between monitoring points based on the lake flow direction relationship; for any target monitoring point, identifying its upstream monitoring points according to the spatial propagation weight matrix, and weighting and superimposing the local warning indices of the upstream monitoring points onto the local warning indices of the target monitoring point to obtain an early warning comprehensive index enhanced by spatial propagation; wherein, the lake flow direction relationship is determined based on the velocity vector in the lake flow monitoring data; The spatial propagation weight is determined based on the angle between the direction vector from the upstream monitoring point to the target monitoring point and the direction of the lake flow at the upstream monitoring point, as well as the distance between the two monitoring points. The spatial propagation weight is positively correlated with the cosine of the angle and decreases exponentially with increasing distance.

2. The method according to claim 1, characterized in that, The lake current-water quality response model was constructed based on the physically constrained kernel regression method, specifically including: Call the optimization objective function that includes data fitting terms and physical constraint penalty terms; Define physical constraint penalty terms, which include at least: a transport dilution constraint based on the transport dilution principle, used to penalize the case where the partial derivative of the response model with respect to the standardized value of velocity is positive; and a retention accumulation constraint based on the retention accumulation principle, used to penalize the case where the partial derivative of the response model with respect to the standardized value of residence time is positive. By minimizing the objective function, a lake current-water quality response model that satisfies physical constraints is obtained.

3. The method according to claim 1, characterized in that, Based on the lake current-water quality response model, dynamic early warning thresholds that vary with the comprehensive characteristic index of the lake current are determined, including: Historical water quality monitoring data are divided into several data groups according to the magnitude of the corresponding lake flow comprehensive characteristic index. Calculate the conditional quantiles of water quality concentrations for each data group at the preset risk level; Based on this, the quantile regression method was used to fit the functional relationship between the dynamic early warning threshold and the comprehensive characteristic index of the lake flow, which served as the basis for determining the dynamic early warning threshold.

4. The method according to claim 1, characterized in that, Critical slowdown indicators include variance growth rate, autocorrelation enhancement rate, and skewness change rate; The critical slowdown index is calculated by: extracting time series segments of water quality monitoring data using a sliding window, calculating the variance, lag 1 autocorrelation coefficient and skewness within each window, calculating the growth trend of statistical characteristics over time based on the Mann-Kendall trend test method, and obtaining the variance growth rate, autocorrelation enhancement rate and skewness change rate, respectively.

5. The method according to claim 1, characterized in that, The decision-making process integrates routine early warning levels and sudden change early warning levels, specifically including: When both the regular warning level and the sudden change warning level are normal, the final warning level is determined to be normal. When only one of the regular warning level and the sudden warning level triggers a warning, the final warning level is determined to be the level that triggered the warning. When both the regular warning level and the sudden change warning level are triggered simultaneously, it is determined to be a compound risk. The higher of the two warning levels is selected, and the warning level is raised by one level based on the higher warning level.

Images