An irrigation district data anomaly analysis method and server

Abstract

The application discloses an irrigation area data anomaly analysis method and a server, relates to the field of data anomaly analysis, and comprises the following steps: acquiring irrigation area canal system topological data and sensor deployment information, and constructing a topological causal diagram; acquiring time series data for a sensor pair in the topological causal diagram, and determining a space-time lag parameter; constructing a physical conduction prediction model; acquiring actual measurement time series data of a downstream sensor and real-time time series data of an upstream sensor, and generating physical expected time series data of the downstream sensor by using the physical conduction prediction model according to the real-time time series data and the space-time lag parameter; and calculating a residual sequence based on the actual measurement time series data and the physical expected time series data of the downstream sensor, and determining whether the data of the downstream sensor is abnormal based on the residual sequence. The application can effectively solve the problem of high false alarm rate in irrigation area application.

Description

Technical Field

[0001] This application relates to the field of data anomaly analysis, and in particular to a method and server for analyzing data anomalies in irrigation districts. Background Technology

[0002] With the rapid development of information technology and intelligentization in modern agriculture, automated monitoring systems for irrigation districts have become a key technological infrastructure for improving water resource utilization efficiency and ensuring food security. Modern irrigation districts typically cover thousands to tens of thousands of hectares of farmland, with intricate canal networks including main canals, branch canals, and distribution canals, as well as numerous gates, pumping stations, and water distribution outlets. To achieve real-time monitoring and precise scheduling of water conditions throughout the irrigation district, the system deploys a large number of different types of sensor devices, including water level gauges, flow meters, gate opening sensors, and water quality monitors. These sensors generate massive amounts of time-series monitoring data daily. This time-series data is the core basis for intelligent management of the irrigation district, and its accuracy and reliability directly affect the scientific nature of irrigation decisions. However, due to the complexity and harshness of the irrigation district environment, sensor devices frequently face various failure risks: zero-point drift due to equipment aging, signal interference caused by extreme weather, measurement deviations caused by weeds or silt blockage, and data loss due to unstable communication links. If these abnormal data are not detected and processed in a timely manner, they will seriously mislead irrigation scheduling decisions, potentially leading to localized droughts or floods, resulting in agricultural losses and water waste. Therefore, building an efficient and accurate data anomaly detection system has become a core technical problem that urgently needs to be solved in the field of automated irrigation monitoring.

[0003] Traditional anomaly detection techniques are primarily based on statistical principles, treating each sensor as an independent monitoring object for analysis. In terms of statistical methods, the most common is the 3σ criterion based on the normal distribution assumption. This method calculates the mean and standard deviation of historical data, marking data points that deviate from the mean by more than three times the standard deviation as anomalies.

[0004] However, the aforementioned traditional technologies have a fundamental flaw: they treat each sensor within the irrigation area as an independent entity, completely ignoring the inherent physical connections and spatiotemporal transmission laws inherent in the irrigation area as a complete water conservancy system. Irrigation systems follow strict physical conservation laws, and there are clear causal relationships and spatiotemporal lag effects between upstream and downstream sensors. For example, when an upstream gate opens, the downstream water level and flow rate will inevitably increase accordingly after a certain period. While this change may statistically appear abnormal, deviating from historical averages, it is logically perfectly normal. Traditional independent detection methods cannot recognize this logical normality, leading to a large number of false alarms. Therefore, traditional anomaly detection technologies suffer from a high false alarm rate in irrigation area applications and cannot meet the stringent data quality requirements of intelligent irrigation management. Summary of the Invention

[0005] This application provides a method and server for analyzing anomalies in irrigation district data, which effectively solves the problem of high false alarm rates in irrigation district applications.

[0006] To achieve the above objectives, the embodiments of this application adopt the following technical solutions:

[0007] Firstly, a method for anomaly analysis of irrigation district data is provided, applied to a server. This method includes:

[0008] Obtain irrigation district canal system topology data and sensor deployment information, and construct a topological causal graph based on the irrigation district canal system topology data and sensor deployment information to characterize the upstream and downstream water flow relationship between multiple sensors in the irrigation district. In the topological causal graph, the nodes are sensors and the directed edges are direct water flow paths between sensors.

[0009] For sensor pairs connected by directed edges in a topological causal graph and having direct upstream and downstream relationships, time-series data of the sensor pairs under historical normal operating conditions are obtained, and spatiotemporal lag parameters are determined based on the time-series data. The spatiotemporal lag parameters are used to characterize the time required for upstream data changes to be transmitted to the downstream.

[0010] For at least one downstream sensor, a physical conduction prediction model is constructed based on the historical normal time-series data of one or more upstream sensors and the corresponding spatiotemporal lag parameters between the downstream sensor and one or more upstream sensors. The physical conduction prediction model is used to predict the expected physical time-series data of the downstream sensor based on the time-series data of the upstream sensor.

[0011] During the real-time monitoring phase, the actual measurement time series data of downstream sensors and the real-time time series data of one or more upstream sensors are acquired. Then, using the physical conduction prediction model, the expected physical time series data of downstream sensors are generated based on the real-time time series data and the spatiotemporal lag parameters.

[0012] Based on the actual measurement time series data and the physical expected time series data of the downstream sensor, a residual sequence is calculated, and the presence of anomalies in the data of the downstream sensor is determined based on the residual sequence.

[0013] In one possible implementation of the first aspect, determining the spatiotemporal lag parameter based on time-series data includes:

[0014] Extract the timing data of the upstream sensor and the timing data of the downstream sensor of the sensor pair;

[0015] The correlation coefficient values ​​at different time offsets are obtained by calculating the time series data of upstream and downstream sensors based on the preset cross-correlation function.

[0016] Determine the target time offset that maximizes the correlation coefficient value, and use the target time offset as the spatiotemporal lag parameter between sensor pairs.

[0017] In another possible implementation of the first aspect, before using the target time offset as a spatiotemporal hysteresis parameter between sensor pairs, the following is also included:

[0018] Calculate the spatiotemporal lag parameter differences between adjacent sensor pairs within the same channel level to identify spatiotemporal lag parameters of local anomalies;

[0019] The theoretical propagation time range is calculated based on the preset channel length and design flow rate, and the target time offset is limited to the theoretical propagation time range.

[0020] When the target time offset exceeds a reasonable range, the short-term spatiotemporal lag parameter is determined based on the historical time series data within the recent time window, and the long-term spatiotemporal lag parameter is determined based on the historical time series data within the long-term time window, wherein the historical time series data within the long-term time window covers at least one complete irrigation cycle.

[0021] The short-term and long-term lag parameters are weighted and fused to obtain the corrected spatiotemporal lag parameters.

[0022] In another possible implementation of the first aspect, before calculating the spatiotemporal hysteresis parameter differences between adjacent sensor pairs within the same channel level, the following is included:

[0023] The target time offset is verified in a hierarchical manner, which includes lateral verification between sensors of the same level and longitudinal verification between sensors of different levels.

[0024] The target time offset is subjected to hierarchical verification, including:

[0025] The channel-level hierarchy of sensors is determined based on the topological cause-effect graph, and the sensors within the same channel level and the upstream and downstream sensors across channel levels are determined based on the channel-level hierarchy.

[0026] Construct a cross-validation set based on sensors within the same channel level, and calculate the statistical distribution of spatiotemporal lag parameters of the sensors in the cross-validation set;

[0027] A vertical verification set is constructed based on upstream and downstream sensors across channels, and the changing trend of the spatiotemporal lag parameters of the sensors in the vertical verification set is calculated.

[0028] The target time offset is compared with the statistical distribution of the horizontal validation set and the changing trend of the vertical validation set to determine the rationality of the target time offset.

[0029] Record the verification results and generate a tiered verification report.

[0030] In another possible implementation of the first aspect, the physical conduction prediction model is a regression model; for at least one downstream sensor, a physical conduction prediction model is constructed based on historical normal time-series data from one or more upstream sensors and the corresponding spatiotemporal lag parameters between the downstream sensor and one or more upstream sensors, including:

[0031] The historical normal time series data of downstream sensors are used as the target variable, and the historical normal time series data of one or more upstream sensors after time alignment with spatiotemporal lag parameters are used as the feature variable.

[0032] The regression model is trained based on the target variable and feature variables to determine the model parameters and obtain the trained physical conduction prediction model.

[0033] In another possible implementation of the first aspect, a residual sequence is calculated based on the actual measurement time-series data and the physically expected time-series data from the downstream sensor, including:

[0034] The actual measurement time series data and the physical expected time series data of the downstream sensor are time-aligned, and the difference between the actual measurement time series data and the physical expected time series data is calculated at each time point to obtain the original residual sequence;

[0035] The original residual sequence is standardized to obtain the standardized residual sequence;

[0036] A moving average filter is used to smooth the standardized residual sequence in order to remove high-frequency noise interference;

[0037] The smoothed data is used as a residual sequence.

[0038] In another possible implementation of the first aspect, determining whether there are anomalies in the data from downstream sensors based on the residual sequence includes:

[0039] Based on the residual sequence, an anomaly score sequence is obtained by using a sliding window calculation.

[0040] The dynamic anomaly judgment threshold is obtained, which is based on the historical residual sequence generated by the downstream sensor under historical normal operation conditions.

[0041] The current abnormal score in the abnormal score sequence is compared with the dynamic abnormal judgment threshold. If the current abnormal score is greater than the dynamic abnormal judgment threshold, the data is determined to be abnormal.

[0042] In another possible implementation of the first aspect, obtaining the dynamic anomaly judgment threshold includes:

[0043] Based on the topological causal graph, the spatial proximity relationship of the sensors is obtained, and a neighborhood influence weight matrix is ​​constructed based on the spatial proximity relationship. The neighborhood influence weight is calculated according to the hydraulic connectivity and distance between the sensors.

[0044] Obtain historical residual data of neighboring sensors of downstream sensors within the same time period, and calculate the weighted neighborhood residual benchmark based on the neighborhood influence weight matrix;

[0045] Calculate the historical residual statistical characteristics of downstream sensors, determine the historical residual distribution based on the historical residual statistical characteristics, and fit the historical residual distribution using a probability density function;

[0046] The weighted neighborhood residual benchmark is used as a spatial reference value, and the deviation between the spatial reference value and the center of the historical residual distribution is calculated.

[0047] Spatially correct the historical residual distribution based on the deviation to obtain the spatially corrected residual distribution.

[0048] Calculate the regional residual mean and standard deviation based on the spatially corrected residual distribution;

[0049] The regional residual mean and standard deviation are combined to form a regional anomaly judgment benchmark;

[0050] Based on the regional anomaly judgment benchmark and the preset confidence interval, the dynamic anomaly judgment threshold is calculated.

[0051] In another possible implementation of the first aspect, the method also includes:

[0052] When there are anomalies in the data from downstream sensors, the anomaly type is determined. The anomaly types include sensor malfunction anomalies, data transmission anomalies, and environmental change anomalies.

[0053] Determine the scope of the anomaly's impact based on the anomaly type and the topological cause-effect graph;

[0054] An anomaly report is generated based on the scope of the anomaly's impact, and the anomaly report is sent to the monitoring terminal connected to the server.

[0055] Secondly, this application provides a server, comprising:

[0056] The memory is configured to store instructions; and

[0057] The processor is configured to retrieve the instructions from the memory and, when executing the instructions, to implement the above-described method for analyzing irrigation district data anomalies.

[0058] Compared to traditional independent detection methods based on statistical principles, this technical solution significantly improves the accuracy and reliability of anomaly detection in irrigation districts by constructing a topological causal graph and a physical transmission prediction model, effectively solving the core problem of high false alarm rates in traditional methods. By fully considering the upstream and downstream causal relationships of the irrigation district's water conservancy system, it can accurately distinguish between statistical anomalies and physical anomalies, avoiding misjudging normal data that conforms to the water flow transmission law as anomalies, thereby greatly reducing the false alarm rate and improving the practicality of the anomaly detection system. The introduction of spatiotemporal lag parameters makes anomaly detection more accurate, enabling the establishment of a physical prediction model based on the actual propagation time and attenuation law of water flow in the canal system, making anomaly judgments more consistent with the actual operating characteristics of the irrigation district's water conservancy system. The multi-sensor collaborative analysis mechanism based on the topological causal graph overcomes the limitations of isolated information in traditional single-point detection methods. By integrating time-series data from multiple upstream sensors to predict downstream conditions, it improves the robustness and fault tolerance of anomaly detection, maintaining high detection accuracy even if some sensors malfunction. The residual sequence analysis method identifies anomalies by comparing the deviation between actual measured values ​​and expected values ​​derived from physical laws, making the physical meaning of anomaly detection clearer and the detection results more interpretable and reliable. Furthermore, by automatically extracting spatiotemporal lag parameters and building predictive models through machine learning methods, the workload of manual parameter tuning is reduced, and the adaptive capability and intelligence level are improved. This not only significantly enhances the accuracy and efficiency of anomaly detection but also provides more reliable data quality assurance for intelligent irrigation management, strongly supporting the scientific and precise nature of irrigation decisions.

[0059] Other features and advantages of the embodiments of this application will be described in detail in the following detailed description section. Attached Figure Description

[0060] Figure 1 A flowchart illustrating a method for analyzing anomalies in irrigation district data provided in this application embodiment;

[0061] Figure 2 A schematic diagram of the topology and cause-effect graph of an irrigation canal system provided in this application embodiment;

[0062] Figure 3 This is a schematic diagram illustrating the construction of a physical conduction prediction model provided in an embodiment of this application;

[0063] Figure 4 This is a schematic diagram of an adjacent sensor pair provided in an embodiment of this application. Detailed Implementation

[0064] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only for illustration and explanation of the embodiments of this application and are not intended to limit the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0065] It should be noted that if the embodiments of this application involve directional indicators (such as up, down, left, right, front, back, etc.), the directional indicators are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicators will also change accordingly.

[0066] Furthermore, if the embodiments of this application involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, features defined with "first" or "second" may explicitly or implicitly include at least one of those features. Additionally, the technical solutions of various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. If the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed in this application.

[0067] Figure 1 The illustration shows a flowchart of a method for analyzing anomalies in irrigation district data according to an embodiment of this application. Figure 1 As shown in the figure, this application provides a method for analyzing anomalies in irrigation district data, applied to a server. The method may include the following steps.

[0068] S110. Obtain irrigation district canal system topology data and sensor deployment information, and construct a topological causal graph based on the irrigation district canal system topology data and sensor deployment information to characterize the upstream and downstream water flow relationship between multiple sensors in the irrigation district. The nodes of the topological causal graph are sensors, and the directed edges are direct water flow paths between sensors.

[0069] S120. For sensor pairs connected by directed edges in the topological causal graph and having direct upstream and downstream relationships, obtain the time series data of the sensor pairs under historical normal operation, and determine the spatiotemporal lag parameter based on the time series data. The spatiotemporal lag parameter is used to characterize the time required for upstream data changes to be transmitted to the downstream.

[0070] S130. For at least one downstream sensor, a physical conduction prediction model is constructed based on the historical normal time series data of one or more upstream sensors and the corresponding spatiotemporal lag parameters between the downstream sensor and one or more upstream sensors. The physical conduction prediction model is used to predict the expected physical time series data of the downstream sensor based on the time series data of the upstream sensor.

[0071] S140. In the real-time monitoring stage, the actual measurement time series data of the downstream sensor and the real-time time series data of one or more upstream sensors are acquired, and the physical transmission prediction model is used to generate the physical expected time series data of the downstream sensor based on the real-time time series data and the spatiotemporal lag parameters.

[0072] S150. Based on the actual measurement time series data and the physical expected time series data of the downstream sensor, the residual sequence is calculated, and the abnormality of the data of the downstream sensor is determined based on the residual sequence.

[0073] In practical implementation, the server first obtains complete topological structure information of the irrigation district from multiple data sources. The irrigation district canal system topology data mainly comes from engineering design drawings, Geographic Information System (GIS) databases, and field survey records, including detailed parameters such as spatial location, length, elevation, and cross-sectional shape of main canals, branch canals, and distribution canals. Sensor deployment information records the specific installation location, equipment model, measurement parameter type, and correspondence with the canals for each sensor. Based on this fundamental data, the server uses graph theory algorithms to construct a topological causal graph, treating each sensor as a node in the graph, and establishing directed edges between nodes according to the natural flow direction of water, such as... Figure 2 As shown in the figure, the irrigation district's canal system structure includes main canals and branch canals. Circular nodes represent sensors (such as water level gauges and flow meters) at different locations, and directed edges between nodes represent direct water flow paths and upstream-downstream causal relationships between sensors. For example, water flowing from a water level gauge directly affects a downstream flow meter.

[0074] In the specific construction process, the water flow direction is first determined based on the channel's elevation information and design slope. Then, the connection relationships of the sensors are traced along the water flow path. For example, if sensor A is located upstream of a branch canal and sensor B is located 500 meters downstream of the same branch canal, and there are no other branching or confluence nodes between them, a directed edge is established between A and B, indicating that the water flows directly from A to B. For complex branching and confluence nodes, the server establishes many-to-many connections based on hydraulic principles. The topological causal graph constructed in this way not only accurately reflects the physical structure of the irrigation area but, more importantly, clarifies the causal relationships between the sensors, providing a physical basis for subsequent anomaly detection.

[0075] For each pair of sensors with a direct upstream-downstream relationship identified in the topological cause-effect graph, the server needs to accurately calculate the spatiotemporal lag parameters between them. The spatiotemporal lag parameters reflect the time required for a water flow change signal to propagate from the upstream sensor location to the downstream sensor location. This time is affected by various factors such as channel length, flow velocity, water depth, and channel roughness.

[0076] The calculation process employs cross-correlation function analysis. In practice, the server first extracts time-series data of the sensor pairs during normal operation from the historical database. The time span is typically chosen to include data from multiple complete irrigation cycles to ensure statistical reliability. Then, the correlation coefficient between the two time series is calculated at different time offsets τ.

[0077] ;

[0078] in and Let represent the measurements taken by the upstream and downstream sensors at time t, respectively, and N be the total number of data points. The server iterates through different values ​​of τ to find the time offset that maximizes the correlation coefficient R(τ). This offset is the desired spatiotemporal lag parameter. For example, if the opening of the upstream gate causes an increase in flow, and the downstream sensor observes a corresponding rise in water level 2 hours later, then the spatiotemporal lag parameter for that sensor pair is 2 hours. This parameter determination method based on actual observation data can accurately capture the comprehensive influence of various complex factors on water flow propagation, significantly improving the accuracy of subsequent prediction models.

[0079] Based on historical normal operation data and established spatiotemporal lag parameters, the server constructs a physical conduction prediction model for each downstream sensor. The physical conduction prediction model employs multiple linear regression or nonlinear regression methods. Its basic principle is to use real-time data from one or more upstream sensors, combined with spatiotemporal lag effects, to predict the expected physical values ​​that the downstream sensors should observe at the current moment. Figure 3 This illustration shows a schematic diagram of the construction of a physical conduction prediction model provided in an embodiment of this application, such as... Figure 3 As shown, during model construction, the historical time-series data of each upstream sensor are first time-aligned according to the corresponding spatiotemporal lag parameters to ensure that the upstream change signal and the downstream response signal correctly correspond on the time axis. Then, the time-aligned upstream data is used as the feature variable, and the actual observation data of the downstream sensor is used as the target variable. The regression model is trained using the least squares method or gradient descent algorithm. For example, for a certain downstream water level sensor, its prediction model might be expressed as:

[0080] ;

[0081] in Downstream expected water level, , , These are the upstream flow rate, water level, and gate opening, respectively. , , For the corresponding spatiotemporal lag parameters, , , , These are the model parameters. After the model is trained, the server will conduct validation tests to ensure that the prediction accuracy meets the requirements. This prediction model based on physical causality can accurately reflect the inherent laws of the irrigation area's water conservancy system, providing a physical benchmark for anomaly detection.

[0082] During the real-time monitoring phase, the server continuously receives real-time data streams from various sensors and uses the established physical conduction prediction model to generate the expected physical time-series data for downstream sensors. Specifically, the server establishes a data buffering mechanism to collect the latest measurement data from each upstream sensor in real time and performs time correction on this data according to pre-determined spatiotemporal lag parameters. For example, if an upstream sensor measures a flow rate of 5 cubic meters per second at 10:00, and its corresponding spatiotemporal lag parameter is 1.5 hours, this data will be used to predict the expected physical value of the downstream sensor at 11:30. The server processes the continuous data stream through a sliding window mechanism, immediately triggering the prediction calculation process whenever new upstream data arrives.

[0083] The prediction calculation employs a pre-trained physical transmission prediction model, inputting current upstream real-time data into the model and outputting the expected physical values ​​of downstream sensors. To improve prediction stability, the server also considers the combined influence of multiple upstream sensors, generating the final expected value through weighted averaging or multiple regression. This real-time prediction mechanism ensures that the anomaly detection system can obtain accurate physical benchmarks in a timely manner, providing a reliable basis for subsequent anomaly judgment.

[0084] Based on a comparison between actual measured data and expected physical data, the server calculates a residual sequence and uses this to determine if there are any anomalies in the sensor data. The residual calculation process first aligns the actual measured time series data and the expected physical time series data to ensure accurate correspondence at each time point. Then, it calculates the difference between the two point by point to obtain the original residual sequence. Where e(t) is the residual at time t, These are actual measured values. This is the expected physical value.

[0085] To eliminate the influence of units and numerical ranges, the server performs standardization on the original residual sequence: ,in and These represent the mean and standard deviation of the residual sequence, respectively. Next, a moving average filter is used to smooth the standardized residuals, removing the effects of high-frequency noise and random disturbances. Anomaly detection employs a dynamic threshold mechanism. The server dynamically adjusts the anomaly detection threshold based on the statistical characteristics of historical residual sequences and the residual performance of neighboring sensors. When the residual value at a certain moment exceeds the dynamic threshold, the data at that moment is considered anomaly. For example, if the actual measurement value of a water level sensor is 30 cm higher than the expected physical value, and this deviation persists for more than 10 minutes, while also exceeding the anomaly threshold determined based on historical data, then the sensor data is considered anomaly. This anomaly detection method based on physical expectations can effectively distinguish between statistical anomalies and physical anomalies, significantly reducing the false alarm rate and improving the accuracy and practicality of detection.

[0086] This embodiment achieves a fundamental shift from traditional statistical methods to physically driven methods by constructing a topological causal graph and a physical transmission prediction model based on physical causal relationships, effectively solving the key problem of high false alarm rates in irrigation area anomaly detection. Precise calculation of spatiotemporal lag parameters ensures accurate modeling of causal relationships between upstream and downstream sensors, making anomaly detection more consistent with the true physical characteristics of the water conservancy system. A multi-sensor collaborative prediction mechanism enhances the robustness and fault tolerance of the detection system; even if some sensors fail, the entire system can still maintain high detection accuracy. Dynamic threshold adjustment and residual sequence analysis methods improve the scientific rigor and reliability of anomaly judgment, making the detection results more physically meaningful and interpretable. Real-time monitoring and prediction mechanisms guarantee the timeliness of anomaly detection, meeting the real-time requirements of intelligent irrigation area management. The entire technical solution not only significantly improves the accuracy and efficiency of anomaly detection but also provides more reliable data quality assurance for irrigation area automated monitoring systems, strongly supporting scientific and precise irrigation decisions and promoting the efficient use of water resources and the sustainable development of agricultural production.

[0087] In one embodiment of this invention, determining the spatiotemporal lag parameter based on time-series data includes the following steps:

[0088] S210. Extract the timing data of the upstream sensor and the timing data of the downstream sensor of the sensor pair;

[0089] S220. Calculate the time series data of the upstream sensor and the downstream sensor based on the preset cross-correlation function to obtain the correlation coefficient values ​​under different time offsets;

[0090] S230. Determine the target time offset that maximizes the correlation coefficient value, and use the target time offset as the spatiotemporal lag parameter between sensor pairs.

[0091] In determining the spatiotemporal lag parameters, the server needs to extract time-series data for specific sensor pairs from the historical database. The data extraction process first identifies sensor combinations with direct upstream-downstream causal relationships based on the sensor pair relationships determined in the topological causal graph. For example, on the main canal of an irrigation district, sensor S001 is located at the upstream gate, and sensor S002 is located at the water diversion point 2 kilometers downstream; the two constitute a typical sensor pair.

[0092] When extracting data, the server sets a time window and data quality standards, typically selecting historical data covering at least 3-6 complete irrigation cycles to ensure coverage of various typical operating conditions. During data extraction, the server automatically filters out obvious outliers, missing values, and data from equipment maintenance periods, retaining only high-quality time-series data from normal operating conditions. For water level sensors, the extracted data includes water level readings every 5 minutes; for flow sensors, it extracts flow measurement values ​​at corresponding time intervals.

[0093] To ensure data time synchronization, the server aligns the data according to the timestamps of each sensor, eliminating time errors caused by clock skew or data transmission delays. The extracted data undergoes preprocessing, including outlier detection, interpolation imputation, and smoothing filtering, to ensure the accuracy of subsequent cross-correlation calculations.

[0094] The correlation coefficient is calculated based on a preset cross-correlation function. This process reveals the time dependency between upstream and downstream sensor data through mathematical statistical methods. The implementation principle of the cross-correlation function is based on signal processing theory. By calculating the correlation strength of two time-series signals at different time offsets, the time delay characteristics of signal propagation are identified.

[0095] In the specific calculation process, the server uses a standardized cross-correlation function: ,in and Let τ and N represent the standardized data values ​​of the upstream and downstream sensors at time t, respectively, where τ is the time offset and N is the data sequence length. , Let be the mean of their respective sequences. , The standard deviation is denoted as .

[0096] The calculation process employs a sliding window method, with the time offset τ starting from 0 and increasing in fixed steps to a preset maximum value. Typically, the maximum offset is set to 1.5 times the theoretical maximum propagation time estimated based on the channel length and minimum flow velocity. For example, for a 3-kilometer-long channel section, assuming a minimum flow velocity of 0.1 meters per second and a theoretical maximum propagation time of approximately 8.3 hours, the maximum time offset is set to 12 hours. At each offset τ, the server calculates the corresponding correlation coefficient value. This value reflects the degree of linear correlation between upstream data at time t and downstream data at time t+τ.

[0097] To improve computational accuracy and noise resistance, the server also employs windowing functions to preprocess the data. Commonly used windowing functions include the Hanning window or the Hamming window to reduce the impact of boundary effects and random noise. The entire calculation process generates a curve showing the correlation coefficient changing with the time offset; the peak position of this curve corresponds to the optimal spatiotemporal lag parameter.

[0098] After the correlation coefficient is calculated, the server needs to identify and determine the target time offset that will maximize the correlation coefficient. The determination of the target time offset uses a peak detection algorithm, which can automatically identify the global maximum point in the correlation coefficient curve.

[0099] In practice, the server first smooths the calculated correlation coefficient sequence, using a moving average filter to eliminate high-frequency noise and local fluctuations, ensuring the identification of true peaks. Then, a peak detection algorithm is applied, which considers not only numerical magnitude but also the significance and stability of the peaks. For example, if the correlation coefficient curve reaches its maximum value of 0.85 at a time offset of 2.5 hours, and this peak has a significant advantage over surrounding local maximum values ​​(at least 0.1 higher), then 2.5 hours is determined as the target time offset.

[0100] To avoid erroneous peak identification, the server sets a minimum peak threshold, typically requiring the maximum correlation coefficient to be no less than 0.6, ensuring that the detected lag relationship is statistically significant. In some complex cases, the correlation coefficient curve may show multiple high peaks. In such cases, the server will combine prior physical knowledge to make a judgment and prioritize the peak that falls within a reasonable time range.

[0101] Once the target time offset is determined, this value is set as the spatiotemporal lag parameter of the sensor pair and stored in the parameter database for use in subsequent physical conduction prediction models. Simultaneously, the server records the maximum value of the correlation coefficient as a reliability indicator for this parameter, used for weight allocation and uncertainty assessment during subsequent model training.

[0102] This embodiment achieves accurate quantification of the time delay relationship between upstream and downstream sensors through data extraction, cross-correlation function calculation, and peak identification. Quality control in the data extraction stage ensures the accuracy of the computational basis and effectively avoids the adverse effects of abnormal data on parameter estimation. The application of standardized cross-correlation functions accurately captures the propagation characteristics of water flow signals in the canal system, reflecting the true physical conduction laws. The intelligent peak detection algorithm, combining statistical principles and physical constraints, ensures the accuracy and stability of target time offset identification. Through precise spatiotemporal lag parameters, the subsequent physical conduction prediction model can more accurately reflect the dynamic characteristics of the irrigation district's water conservancy system, thereby significantly improving the accuracy of anomaly detection.

[0103] In one embodiment of this invention, before using the target time offset as the spatiotemporal lag parameter between sensor pairs, the following steps are further included:

[0104] S310. Calculate the difference in spatiotemporal lag parameters between adjacent sensor pairs within the same channel level to identify spatiotemporal lag parameters of local anomalies.

[0105] S320. Calculate the theoretical propagation time range based on the preset channel length and design flow rate, and limit the target time offset within the theoretical propagation time range;

[0106] S330. When the target time offset exceeds a reasonable range, determine the short-term spatiotemporal lag parameter based on the historical time series data within the recent time window, and determine the long-term spatiotemporal lag parameter based on the historical time series data within the long-term time window, wherein the historical time series data within the long-term time window covers at least one complete irrigation cycle.

[0107] S340. Weighted fusion of short-term and long-term lag parameters yields the corrected spatiotemporal lag parameters.

[0108] The short-term and long-term lag parameters are weighted and fused to obtain the corrected spatiotemporal lag parameters, including:

[0109] S341. Calculate the short-term weight coefficient based on the data quality and coverage of the short-term time window, and calculate the long-term weight coefficient based on the data stability and integrity of the long-term time window.

[0110] S342. Based on the stage characteristics of the current irrigation cycle, dynamically adjust the ratio of short-term weight coefficient and long-term weight coefficient, wherein the short-term weight coefficient is increased during the peak irrigation period and the long-term weight coefficient is increased during the stable irrigation period.

[0111] S343. Using a weighted average algorithm, the short-term and long-term spatiotemporal lag parameters are weighted and fused according to the adjusted short-term and long-term weight coefficients.

[0112] S344. Smooth the weighted fusion results to eliminate the impact of sudden changes during weight switching.

[0113] S345. Use the smoothed result as the corrected spatiotemporal lag parameter.

[0114] During the verification of spatiotemporal lag parameters, the server identifies potential local anomalies by calculating the differences in spatiotemporal lag parameters between adjacent sensor pairs within the same channel level. In specific implementation, the server first classifies the sensors according to different levels such as main canals, branch canals, and distribution canals based on the channel level information in the topological causal graph, and then identifies spatially adjacent sensor pairs within each level.

[0115] like Figure 4 As shown, for example, on a main canal, sensors S001, S002, and S003 are arranged sequentially in the direction of water flow, forming two adjacent sensor pairs: (S001, S002) and (S002, S003). For these adjacent sensor pairs, since they are located in the same canal environment and have similar hydraulic conditions and physical characteristics, their spatiotemporal lag parameters should exhibit a certain continuity and similarity.

[0116] The server calculates the differences in spatiotemporal lag parameters between adjacent sensors: ,in Let represent the spatiotemporal lag parameter of the i-th sensor pair. An anomaly detection threshold is established by statistically analyzing the distribution characteristics of differences between all adjacent parameters within the same channel level. When a parameter difference significantly exceeds the normal range, for example, if the spatiotemporal lag parameters of two adjacent sensor pairs are 2.1 hours and 4.8 hours respectively, with a difference of 2.7 hours, which far exceeds the average difference of parameters within that channel level plus twice the standard deviation, it is marked as a locally anomalous parameter.

[0117] Calculating the theoretical propagation time range based on hydraulic theory verifies the physical rationality of the spatiotemporal lag parameters. The server extracts detailed channel parameters for each sensor pair from the irrigation district's basic database, including channel length L, channel cross-sectional shape, design slope i, and roughness coefficient n. Based on these parameters, the Manning formula is used to calculate the channel's design flow velocity range. , where R is the hydraulic radius, which is calculated based on the channel geometry and the design water depth.

[0118] Considering the variations in flow rate during actual operation, the server will set a reasonable flow rate range, typically 0.3 to 3 times the designed flow rate, to cover various operating conditions from low traffic to full load. The theoretical propagation time range is calculated by dividing the length by the flow rate range. , For example, for a channel section with a length of 1500 meters, the design flow velocity is 0.8 meters per second. Considering that the flow velocity varies from 0.24 to 2.4 meters per second, the theoretical propagation time ranges from 10.4 minutes to 104 minutes.

[0119] The server compares the target time offset calculated by the cross-correlation function with this theoretical range. If the offset falls within a reasonable range, the parameter is considered to have physical reliability; if it exceeds the range, the subsequent parameter correction process is triggered.

[0120] When the target time offset exceeds the theoretically reasonable range, the server initiates a multi-timescale analysis mechanism to recalculate the spatiotemporal lag parameters based on both short-term and long-term historical data. The calculation of short-term spatiotemporal lag parameters is based on historical data within the recent time window, typically selecting data from the last 1 to 2 months. This data reflects the current operating status and environmental conditions and is highly sensitive to seasonal changes and equipment status changes.

[0121] Short-term parameters are calculated using the same cross-correlation function method as described above, but with a smaller data window and a higher parameter update frequency, enabling rapid response to changes in system status. Long-term spatiotemporal lag parameters are calculated based on a larger time window, with data covering at least one complete irrigation cycle, typically 6 to 12 months. This data range can encompass various typical operating conditions, including different stages such as spring irrigation, summer irrigation, and autumn irrigation, as well as different hydrological conditions such as dry season and wet season.

[0122] Long-term parameters exhibit better statistical stability, reflecting the inherent physical relationship between sensor pairs and are less susceptible to short-term fluctuations. For example, the short-term spatiotemporal lag parameter of a sensor pair is 1.8 hours, reflecting the propagation characteristics under current high flow velocities; while the long-term parameter is 2.3 hours, representing the propagation time under historical average conditions.

[0123] The calculation of short-term and long-term weighting coefficients comprehensively considers multiple dimensions of indicators, including data quality, coverage, stability, and completeness. The calculation of short-term weighting coefficients is primarily based on the quality and coverage of data within a short time window. In practice, the server first assesses the completeness of the short-term data. ,in For the number of valid data points, This represents the total number of theoretical data points.

[0124] Then, data quality metrics are evaluated, including data consistency, continuity, and the proportion of outliers. The formula for calculating the short-term weighting coefficient is: Where α is the basic weighting factor, Score the data quality. This represents the proportion of outliers. The calculation of long-term weighting coefficients focuses on data stability and historical representativeness, employing a similar calculation framework: Where β is the long-term fundamental weighting factor, As an indicator of long-term data stability, Assess historical representativeness. For long-term data integrity.

[0125] The weight ratios are dynamically adjusted based on the stage characteristics of the current irrigation cycle to optimize the parameter fusion effect. The server first establishes a stage identification model for the irrigation cycle, dividing a complete irrigation year into different operational stages, including the spring preparation period, the spring irrigation peak period, the summer management period, the summer irrigation peak period, the autumn closing period, and the winter maintenance period. Each stage has different hydrological characteristics, flow patterns, and operational rules.

[0126] The dynamic adjustment algorithm automatically identifies the current irrigation stage based on the current time and historical operating patterns, and adjusts the weight allocation strategy accordingly. During peak irrigation periods, due to frequent changes in water flow conditions and rapid updates in system operating status, a short-term weight coefficient is added to allow parameters to quickly respond to the current actual operating state. The adjustment formula is: , ,in This is a peak period adjustment factor. This is an indicator of peak intensity.

[0127] During periods of stable irrigation, the system operates relatively steadily. At this time, the long-term weighting coefficient is increased to improve the stability of the parameters by utilizing historical experience. For example, during the peak spring irrigation season, the short-term weighting is adjusted from 0.82 to 0.95, and the long-term weighting from 0.76 to 0.63; while during the winter maintenance period, the short-term weighting is adjusted to 0.68, and the long-term weighting is increased to 0.89.

[0128] The weighted fusion calculation employs an improved weighted average algorithm, comprehensively considering adjusted short-term and long-term weight coefficients to fuse the lag parameters of the two time scales. The standard weighted average formula is: ,in The spatiotemporal lag parameters after fusion. and These are the short-term and long-term spatiotemporal lag parameters, respectively.

[0129] To further improve the fusion effect, the server also introduced a confidence-weighted mechanism, using the maximum correlation coefficient in the calculation process of short-term and long-term parameters as the confidence index: , The revised fusion formula is as follows: .

[0130] To eliminate the potential abrupt changes during weight switching, the server smooths the weighted fusion results. The smoothing process employs an exponential smoothing algorithm, which effectively suppresses sudden parameter jumps and maintains the stationarity of the time series. The exponential smoothing formula is: ,in This is a smoothing coefficient, typically set between 0.1 and 0.3.

[0131] The server dynamically adjusts the smoothing coefficient based on the magnitude and frequency of parameter changes. When a large parameter jump is detected, the smoothing coefficient is automatically reduced to enhance the smoothing effect. Furthermore, an outlier suppression mechanism is employed; when the change in the fused parameters relative to historical smoothed values ​​exceeds a preset threshold, additional smoothing processing is initiated. The smoothed parameters are stored in the parameter database as the final corrected spatiotemporal lag parameters for use by the physical transmission prediction model.

[0132] This embodiment significantly improves the accuracy, stability, and adaptability of spatiotemporal lag parameters through multi-level quality control and parameter fusion techniques. Local consistency checks effectively identify anomalous results in statistical calculations, preventing individual erroneous parameters from negatively impacting overall system performance. Theoretical verification based on physical constraints ensures the engineering rationality of the parameters, preventing unrealistic statistical fitting results. The multi-timescale analysis mechanism fully utilizes the advantages of data from different time windows, maintaining sensitivity to the current state while possessing the stability of historical experience. The dynamic weight adjustment algorithm optimizes the parameter fusion strategy based on the stage characteristics of the irrigation cycle, making the parameters more consistent with the actual operating patterns of the irrigation district.

[0133] In one embodiment of this example, before calculating the spatiotemporal hysteresis parameter differences between adjacent sensor pairs within the same channel level, the following steps are included:

[0134] S410. Perform hierarchical verification of the target time offset, wherein hierarchical verification includes lateral verification between sensors of the same level and longitudinal verification between sensors of different levels.

[0135] The target time offset is subjected to hierarchical verification, including:

[0136] S411. Determine the channel-level hierarchy of sensors based on the topological cause-effect graph, and determine the sensors within the same channel level and the upstream and downstream sensors across channel levels based on the channel-level hierarchy.

[0137] S412. Construct a horizontal verification set based on sensors within the same channel level, and calculate the statistical distribution of the spatiotemporal lag parameters of the sensors in the horizontal verification set.

[0138] S413. Construct a vertical verification set based on upstream and downstream sensors at the cross-channel level, and calculate the changing trend of the spatiotemporal lag parameters of the sensors in the vertical verification set.

[0139] S414. Compare the target time offset with the statistical distribution of the horizontal validation set and the changing trend of the vertical validation set to determine the rationality of the target time offset.

[0140] S415. Record the verification results and generate a graded verification report.

[0141] In the hierarchical verification process of spatiotemporal lag parameters, the server needs to establish a comprehensive verification framework that systematically examines the target time offset through both horizontal and vertical dimensions. The core idea of ​​hierarchical verification is based on the hierarchical structural characteristics of irrigation district water conservancy systems. Channels at different levels have different hydraulic characteristics and propagation laws, while channels at the same level have similar physical properties.

[0142] Horizontal verification primarily examines the consistency and reasonableness of the target time offset within the same channel level. This is achieved by statistically comparing parameters with those of other sensors at the same level to identify outliers that significantly deviate from the normal range. Vertical verification, on the other hand, focuses on the parameter variation patterns between different channel levels, verifying whether the target time offset conforms to the physical transmission characteristics from higher to lower channel levels. For example, the propagation of water flow from the main canal to branch canals and then to distribution canals typically shows a trend of shorter propagation times at lower channel levels, because lower-level channels usually have smaller cross-sections and higher flow velocities.

[0143] Determining the channel hierarchy of sensors based on topological causal graphs is a fundamental technical step in graded verification. The server first extracts channel hierarchy information from the topological causal graph. This information is typically defined during the irrigation district design phase, including different levels such as primary channels (main canals), secondary channels (branch canals), tertiary channels (distributary canals), and quaternary channels (farm canals). The determination of the hierarchy is primarily based on the channel's design flow rate, service area, and importance within the overall irrigation system.

[0144] For example, main canals typically bear the primary water supply for the entire irrigation district, with the largest designed flow rate and the widest service area; branch canals divert water from the main canals to serve specific areas; and distribution canals and farm canals are further subdivided, directly serving the fields. After determining the hierarchical structure, the server establishes a mapping relationship between sensors and canal levels, with each sensor categorized into the corresponding canal level based on its installation location. For sensors installed at canal intersections, they are categorized into higher-level canals using the principle of "highest priority."

[0145] The construction and statistical distribution calculation of lateral verification sets are core technical steps in verifying the consistency of sensor parameters at the same level. Based on a defined channel-level hierarchy, the server groups sensors belonging to the same level according to their spatial location and water flow path, constructing lateral verification sets. For example, the verification set at the main canal level includes all sensor pairs installed on the main canal, while the verification set at the branch canal level includes sensor pairs on each branch canal.

[0146] For each validation set, the server collects the spatiotemporal lag parameters of all sensor pairs within it, forming a parameter sample set. The statistical distribution is calculated using various descriptive statistical indicators, including mean, standard deviation, median, quartiles, skewness, and kurtosis. Specifically, the parameter samples are first preprocessed to remove obvious outliers and missing values, and then the basic statistics are calculated: mean. Standard deviation .

[0147] To more accurately describe the distribution characteristics, quantile information was also calculated, including the 25th quantile Q1, the median Q2, and the 75th quantile Q3, as well as the interquartile range (IQR) based on these quantiles: IQR = Q3 - Q1. Furthermore, the Shapiro-Wilk test was used to determine whether the parametric distribution conformed to a normal distribution. If it did not conform to a normal distribution, nonparametric statistical methods were employed for subsequent analysis.

[0148] The construction and trend analysis of the longitudinal validation set focuses on revealing the regular changes in spatiotemporal lag parameters across channel levels. Based on the water flow path in the topological causal graph, the server identifies sensor pairs with direct upstream-downstream relationships but belonging to different channel levels, and constructs the longitudinal validation set. These sensor pairs typically appear at channel branching points, confluence points, or cascade connections, such as sensor pairs at the end of a main canal and the beginning of a branch canal, or sensor pairs at the end of a branch canal and the beginning of a branch canal.

[0149] After collecting sufficient parameter data from sensor pairs across different levels, the server uses trend analysis to calculate the patterns of change. First, the sensor pairs are sorted according to channel level hierarchy (from highest to lowest), and then the average rate of change of parameters between adjacent levels is calculated. A model of the relationship between levels and parameters can be established using linear regression or nonlinear fitting methods. Commonly used models include linear models. or exponential model .

[0150] Comparative analysis of the target time offset and the validation baseline is a crucial step in determining the reasonableness of the parameters. For cross-validation, the server compares the target time offset with the statistical distribution of the corresponding level of validation sets, primarily using the standardized score method. ,in Let |Z| be the target time offset, and μ and σ be the mean and standard deviation of the parameters at the same level, respectively. According to the 3σ rule of normal distribution, the parameter is considered to be within the normal range when |Z| < 2, marked as suspicious when 2 ≤ |Z| < 3, and judged as abnormal when |Z| ≥ 3.

[0151] For non-normally distributed cases, an outlier detection method based on quartiles is used, with an upper boundary set at Q3 + 1.5 × IQR and a lower boundary at Q1 - 1.5 × IQR. Parameters exceeding this range are marked as suspicious. For longitudinal validation, the server calculates the theoretical expected value based on the sensor pair's level combination using an established trend model. , where f is the fitted trend function, and These are the levels of upstream and downstream sensors, respectively.

[0152] Then calculate the relative deviation between the target value and the expected value: A deviation of less than 20% is considered reasonable, a deviation between 20% and 40% is marked as suspicious, and a deviation exceeding 40% is judged as abnormal. Taking into account the results of cross-sectional and longitudinal verification, a weighted scoring mechanism is used to determine the final level of reasonableness: completely reasonable, basically reasonable, suspicious, or abnormal.

[0153] Recording verification results and generating tiered verification reports are the summary and output stages of the entire verification process. The server establishes a comprehensive verification result recording mechanism, meticulously documenting the verification process and results for each target time offset. Recorded information includes verification time, sensor pair identification, target time offset value, channel level, horizontal verification statistics (mean, standard deviation, Z-score, etc.), vertical verification results (expected value, actual deviation, trend fitting parameters, etc.), overall reasonableness level, anomaly identification, and handling recommendations.

[0154] The tiered validation report uses a structured format and includes an execution summary, validation method description, statistical results summary, list of abnormal parameters, trend analysis charts, and handling recommendations. The statistical summary includes the distribution characteristics of parameters at each level, the number and proportion of abnormal parameters, and the validation pass rate. The list of abnormal parameters details all parameters marked as suspicious or abnormal, including the specific degree of deviation and possible cause analysis. The handling recommendations section proposes specific handling solutions for different types of abnormal parameters, such as recalculation, manual inspection, and equipment maintenance.

[0155] This embodiment significantly improves the quality control level of spatiotemporal lag parameters by establishing a comprehensive verification system with both horizontal and vertical dimensions. The accurate establishment of the channel-level hierarchical structure provides a classification basis for the verification process, ensuring its relevance and systematic nature. Horizontal verification, through statistical distribution analysis of parameters at the same level, effectively identifies anomalous parameters deviating from group characteristics, ensuring the consistency and rationality of parameters from similar sensors. Vertical verification, through trend analysis of parameters across different levels, reveals the inherent physical laws of the irrigation district's water conservancy system, ensuring the logical rationality of parameters across different levels.

[0156] In one embodiment of this invention, the physical conduction prediction model is a regression model. For at least one downstream sensor, a physical conduction prediction model is constructed based on historical normal time-series data from one or more upstream sensors and the corresponding spatiotemporal lag parameters between the downstream sensor and one or more upstream sensors, including the following steps:

[0157] S510. Use the historical normal time series data of downstream sensors as the target variable, and use the historical normal time series data of one or more upstream sensors after time alignment with the spatiotemporal lag parameter as the feature variable.

[0158] S520. Train the regression model based on the target variable and feature variables to determine the model parameters and obtain the trained physical conduction prediction model.

[0159] In the data preparation phase of building a physical conduction prediction model, the server needs to perform time alignment and variable organization on historical time-series data. The target variable is set by selecting actual measurement data from downstream sensors during their historical normal operation as the object that the model needs to learn and predict. This data typically includes indicators such as water level, flow rate, and water quality parameters. The time resolution is set according to monitoring requirements, with common sampling intervals being 5 minutes, 15 minutes, or 1 hour.

[0160] The selection of target variable data requires quality control, including removing outlier data during equipment failures, missing data during sensor maintenance, and outliers that significantly deviate from the normal range. The construction of feature variables requires aligning historical data from one or more upstream sensors with corresponding spatiotemporal lag parameters according to time offsets.

[0161] In practice, if the spatiotemporal lag parameter between upstream sensor A and downstream sensor is 2.5 hours, then the data from upstream sensor A at time t will be used to predict the state of downstream sensor at time t+2.5 hours. This time alignment ensures the correct establishment of causal relationships, enabling the model to learn the true physical conduction laws.

[0162] In the case of multiple upstream sensors, the server will perform time alignment according to their respective spatiotemporal lag parameters to form a multidimensional feature vector. For example, if a downstream water level sensor is affected by the upstream flow sensor (lag of 1.8 hours), the upstream water level sensor (lag of 2.2 hours), and the upstream gate opening sensor (lag of 2.0 hours), then the feature vector corresponding to the target water level at time t+2.5 hours is [Q(t+0.7), H(t+0.3), G(t+0.5)].

[0163] The training process of the regression model is a core technical step in building the physical conduction prediction model. This process automatically extracts the quantitative relationship between upstream and downstream sensors from historical data using machine learning algorithms. The server uses multiple linear regression as the basic model framework, and the model is as follows: , where y is the predicted value from the downstream sensor. Let i be the i-th upstream feature variable. ε represents the corresponding regression coefficient, and ε is the error term.

[0164] The model parameters are determined using the least squares method, by minimizing the sum of squared errors between the predicted and actual values ​​to find the optimal parameters. , where m is the number of training samples. The specific solution process uses the normal equation or gradient descent algorithm. The analytical solution of the normal equation is: , where X is the feature matrix and y is the target vector.

[0165] To improve the model's generalization ability and prevent overfitting, the server also incorporates regularization techniques, including L1 regularization (Lasso regression) and L2 regularization (Ridge regression). The objective function after regularization is: , where λ is the regularization parameter, p=1 corresponds to L1 regularization, and p=2 corresponds to L2 regularization.

[0166] The training process employs cross-validation to select the optimal regularization parameter, typically using 5-fold or 10-fold cross-validation. The dataset is divided into training and validation sets, and a grid search is used to find the λ value that minimizes the validation error. For example, the prediction model of a downstream water level sensor, after training, yields: The coefficients represent the degree of influence of upstream flow, upstream water level, and gate opening on downstream water level, respectively.

[0167] After the model is trained, the server will calculate a number of performance metrics, including the coefficient of determination R², mean squared error MSE, and mean absolute error MAE, to evaluate the model’s fit and prediction accuracy, and ensure that the model meets the accuracy requirements of practical applications.

[0168] During model training, the server also implemented several quality control and optimization measures to ensure the stability of the constructed physical transmission prediction model. The first step was feature engineering optimization, which used methods such as correlation analysis and principal component analysis to identify and remove redundant features, retaining only those that significantly impact the target variable.

[0169] Correlation analysis uses the Pearson correlation coefficient to measure the linear relationship between features: When the correlation coefficient between features is too high (usually >0.8), select the feature with the stronger influence.

[0170] Secondly, outlier handling involves using statistical distribution-based methods (such as the 3σ criterion) or regression residual-based methods (such as studentized residuals) to identify outliers in the training data. Outliers are then appropriately processed or removed to avoid adverse effects on model parameter estimation.

[0171] Third is model diagnosis, which involves testing the reasonableness of the model assumptions through residual analysis, including tests for the normality of residuals (Shapiro-Wilk test), homoscedasticity (Breusch-Pagan test), and independence (Durbin-Watson test). If the test results indicate that the model assumptions are not met, appropriate improvement measures are adopted, such as logarithmic transformation, weighted regression, or time series models.

[0172] Fourth is model stability verification. This involves training the model on data from different periods using a sliding time window approach, comparing the changes in model parameters, and ensuring consistent performance under different operating conditions. For example, the model can be trained using data from spring, summer, and autumn, and the differences in model parameters between seasons can be compared. If the differences are too large, it is necessary to consider introducing seasonal factors or segmented modeling.

[0173] Finally, the predictive performance is evaluated using time-series cross-validation to simulate real-world prediction scenarios, assess the model's real-time predictive capabilities and long-term stability, and ensure that the model not only performs well on historical data but also maintains its predictive accuracy on future data.

[0174] This embodiment successfully established a quantitative prediction model based on physical causality by employing time alignment, regression algorithms, and quality control measures. The time alignment mechanism ensures the correct modeling of the true physical relationships between upstream and downstream sensors, avoiding modeling biases caused by neglecting time delay effects in traditional methods, and significantly improving the model's physical meaning and prediction accuracy. The multivariate regression framework can simultaneously consider the combined influence of multiple upstream factors, more comprehensively reflecting the complexity and variability of irrigation system water conservancy, and improving the model's adaptability and robustness. Regularization techniques and cross-validation methods effectively prevent model overfitting, ensuring good generalization ability and enabling the model to maintain stable performance even when faced with unseen data.

[0175] In one embodiment of this invention, the residual sequence is calculated based on the actual measurement time series data and the physical expected time series data of the downstream sensor, including the following steps:

[0176] S610. Time-align the actual measurement time series data and the physical expected time series data of the downstream sensor, and calculate the difference between the actual measurement time series data and the physical expected time series data at each time point to obtain the original residual sequence.

[0177] S620. Standardize the original residual sequence to obtain the standardized residual sequence;

[0178] S630. A moving average filter is used to smooth the standardized residual sequence in order to remove high-frequency noise interference.

[0179] S640. Use the smoothed data as a residual sequence.

[0180] In the initial stage of residual sequence computation, the server needs to time-align the actual measurement data from downstream sensors with the expected physical data generated by the physical transmission prediction model. The main challenges in the time alignment process include differences in data sampling frequencies, timestamp accuracy, and potential time offsets that may occur during data transmission.

[0181] The server first performs a unified formatting process on the timestamps of the two time series data sets, converting all time information into a standard time format with precision typically maintained at the second or millisecond level. For cases with inconsistent sampling frequencies, interpolation methods are used to unify the data onto the same time grid. Commonly used interpolation methods include linear interpolation, spline interpolation, or nearest neighbor interpolation; the selection criteria are primarily based on the data's variation characteristics and accuracy requirements.

[0182] For example, if actual measurement data is sampled every 5 minutes, while expected data is generated every 10 minutes, then linear interpolation is performed on the expected data to generate the expected value at the same time point as the actual data. After time alignment, the server calculates the residual value for each time point: Where e(t) is the original residual at time t, These are actual measured values. This is the expected physical value.

[0183] This calculation process generates a residual sequence of the same length as the original time-series data. Each residual value represents the degree of deviation between the actual observation and the physical expectation at that moment. A positive residual indicates that the actual value is higher than the expected value, and a negative residual indicates that the actual value is lower than the expected value. The absolute value of the residual reflects the severity of the deviation. For example, if a water level sensor's actual measurement at 10:30 is 2.35 meters and the physical expectation is 2.28 meters, then the original residual at that moment is 0.07 meters.

[0184] Standardization is used to eliminate the influence of dimensions and differences in numerical range. This process transforms the original residual sequence into dimensionless standardized values, facilitating subsequent statistical analysis and threshold comparisons. The core principle of standardization is based on the Z-score transformation in statistics. By calculating the mean and standard deviation of the residual sequence, each residual value is transformed into a standardized position relative to the overall distribution.

[0185] In practice, the server first calculates the statistical characteristics of the original residual sequence, including the mean and standard deviation. Then, it performs a standardization transformation on each residual value to obtain a standardized residual sequence.

[0186] The standardized residuals have a standard normal distribution with a mean of 0 and a standard deviation of 1, which allows the residuals of different sensors and different measurement parameters to be compared and analyzed under the same standard.

[0187] Moving average filtering is a signal processing technique that eliminates high-frequency noise interference and smooths residual sequences. This process suppresses random noise through local averaging while preserving the truly meaningful trends in the residual sequence. The basic principle of the moving average filter is to replace the original value of the center point with the average value of the data points within a sliding window. The filtering formula is: ,in is the filtered value, k is the half-width of the window, and 2k+1 is the window length.

[0188] The choice of window length needs to strike a balance between noise suppression and signal fidelity. A window that is too small will result in ineffective filtering, while a window that is too large may smooth out real anomalous signals. The server employs an adaptive window selection strategy, dynamically adjusting the window size based on the statistical characteristics of the residual sequence. Specifically, the autocorrelation function of the residual sequence is calculated, and the lag time at which the autocorrelation function first drops below 0.1 is selected as a reference value for the window length. For typical 5-minute sampling data, the window length is usually set to 3-7 data points.

[0189] In addition to the simple moving average, the server also implements a weighted moving average filter, giving higher weight to the center point: ,in The weights are typically distributed using a triangular or Gaussian weighting system. To handle boundary effects, padding strategies are used at both ends of the sequence; common methods include zero-padding, boundary value extension, or mirror padding.

[0190] The data processed by moving average filtering is output as the final residual sequence. This sequence integrates the original bias information, standardization, and noise suppression effects, becoming the direct input for anomaly detection algorithms. The final residual sequence has several characteristics: First, data smoothness, as filtering eliminates measurement noise and random disturbances, making true anomaly patterns more prominent and easier to identify; second, standardized numerical range, as all residual values ​​are converted into dimensionless standardized fractions, facilitating the setting of a uniform anomaly judgment threshold; and third, temporal continuity, maintaining the same time resolution and sampling frequency as the original time-series data, ensuring the timeliness of anomaly detection.

[0191] The server performs a quality assessment on the final residual sequence, calculating performance parameters such as the change in signal-to-noise ratio (SNR) before and after filtering, smoothness index, and fidelity of outliers. The SNR is calculated using power spectral analysis, comparing the energy changes of high-frequency noise components before and after filtering. Smoothness is measured by calculating the standard deviation of the differences between adjacent data points; a smaller value indicates a smoother sequence.

[0192] This embodiment constructs high-quality anomaly detection foundation data by organically combining multiple technical aspects such as time alignment, standardization, and noise filtering. Time alignment ensures an accurate correspondence between actual measured values ​​and expected physical values, avoiding false anomalies caused by time deviations and improving the accuracy and physical meaning of residual calculations. Standardization eliminates dimensional differences between different sensors and measurement parameters, establishes a unified anomaly evaluation standard, makes the setting of anomaly detection thresholds more reasonable, and enhances the system's adaptability to different types of sensors. Moving average filtering effectively suppresses measurement noise and environmental interference, significantly improves the signal-to-noise ratio of the residual sequence, makes true anomaly signals more prominent, and reduces the false alarm rate caused by noise.

[0193] In one embodiment of this invention, determining whether there are anomalies in the data from downstream sensors based on the residual sequence includes the following steps:

[0194] S710. Based on the residual sequence, the anomaly score sequence is calculated using a sliding window.

[0195] S720. Obtain the dynamic anomaly judgment threshold, wherein the dynamic anomaly judgment threshold is obtained based on the historical residual sequence generated by the downstream sensor under historical normal operating conditions.

[0196] S730. Compare the current abnormal score in the abnormal score sequence with the dynamic abnormal judgment threshold. If the current abnormal score is greater than the dynamic abnormal judgment threshold, then the data is determined to be abnormal.

[0197] Anomaly score sequences are calculated using a sliding window method based on residual sequences. This process transforms point-like residual values ​​into a comprehensive anomaly score through local statistical analysis, improving the stability and accuracy of anomaly detection. The basic principle of the sliding window method is to set a fixed-length time window on the residual sequence, which moves along the time axis. At each position, the statistical characteristics of the residual data within the window are calculated as the anomaly score for the current time.

[0198] In practice, the server first determines the length of the sliding window. The choice of window length needs to comprehensively consider the sensitivity and stability requirements of anomaly detection. A window that is too short may lead to excessive sensitivity to noise and generate a large number of false alarms; a window that is too long may delay anomaly detection and reduce the system's real-time response capability. It is usually determined based on the sensor's sampling frequency and the expected detection delay. For data with a 5-minute sampling interval, the window length is generally set to 6-12 data points, corresponding to a time span of 30 minutes to 1 hour.

[0199] The outlier score is calculated using a combination of statistical indicators, including the mean, standard deviation, maximum absolute value, and coefficient of variation of the residuals within the window. The formula for calculating the composite outlier score is: Where AS(t) is the anomaly score at time t. The mean of the residuals within the window. The standard deviation of the residuals within the window. The maximum absolute value of the residual within the window. The coefficient of variation is 1. arrive These are the weighting coefficients.

[0200] The weighting coefficients are determined through training with historical data based on the contribution of different statistical indicators to anomaly detection. Typically, the residual mean and maximum absolute value have higher weights, set to 0.3 and 0.4 respectively, while the standard deviation and coefficient of variation have weights of 0.2 and 0.1 respectively. For example, if the residual sequence within a sliding window at a certain time is [0.2, -0.5, 0.8, 1.2, -0.3, 0.6], the calculated mean is 0.17, the standard deviation is 0.53, the maximum absolute value is 1.2, and the coefficient of variation is 3.12. Therefore, the anomaly score is 0.3 × 0.17 + 0.4 × 0.53 + 0.3 × 1.2 + 0.1 × 3.12 = 0.975.

[0201] The dynamic anomaly detection threshold needs to be dynamically adjusted based on the historical normal operating characteristics of downstream sensors to adapt to the impact of different operating conditions and seasonal variations. The calculation of the dynamic threshold is based on the statistical characteristics of the residual sequence under historical normal operating conditions. This historical data typically covers at least one complete irrigation cycle and includes various typical operating conditions and environmental conditions.

[0202] The server first extracts the residual sequence of specified downstream sensors during normal operation from the historical database. The definition of normal operation is based on multiple criteria: normal sensor equipment status, no maintenance work, no extreme weather impact, and stable operation of the irrigation area. The extracted historical residual data is processed using the same method as the current anomaly score calculation, employing the same sliding window length and statistical indicators to generate a historical anomaly score sequence.

[0203] The dynamic threshold is calculated using the quantile method, which determines an appropriate threshold level by analyzing the statistical distribution of historical outlier scores. A commonly used formula is: ,in The threshold is the 95th percentile of historical outlier scores, IQR is the interquartile range, and k is an adjustment factor (typically between 0.5 and 1.5). This method ensures that under normal operating conditions, approximately 95% of outlier scores are below the threshold, thus keeping the false alarm rate low.

[0204] To further improve the adaptability of the thresholds, the server also introduces a time-dependent adjustment mechanism. Considering the seasonal and periodic characteristics of irrigation district operation, the thresholds are set differently for different periods. In practice, historical data is grouped by month or season, and the dynamic threshold for each group is calculated separately. The threshold at the current moment is obtained through time-weighted interpolation. ,in Let t be the weight of historical period i, with the weight being higher for historical periods closer to the current time.

[0205] The comparison between the anomaly score and the dynamic threshold quantitatively determines whether there are anomalies in the sensor data and outputs a clear detection result. The basic logic of the comparison is that if the anomaly score at the current moment is greater than the corresponding dynamic anomaly judgment threshold, the sensor data at that moment is determined to be abnormal; otherwise, the data is considered normal.

[0206] To improve the accuracy of judgments and reduce false alarms, the server implemented several optimization strategies. The first is a continuity check. Considering that occasional numerical fluctuations may lead to single-point misjudgments, abnormal states are required to have a certain duration before being confirmed as genuine anomalies. In practice, a minimum duration threshold is set, typically 15-30 minutes. An anomaly is only officially reported when the anomaly score continuously exceeds the threshold for the specified time.

[0207] Secondly, confidence assessment involves calculating the confidence level for anomalies, in addition to simple binary judgments. The higher the confidence level, the more obvious the anomaly. Thirdly, severity is graded. Based on the degree to which the anomaly score exceeds the threshold, anomalies are divided into four levels: mild, moderate, severe, and extremely severe, providing a reference for subsequent response processing. The grading criteria are: mild anomaly (1-1.5 times the threshold), moderate anomaly (1.5-2 times the threshold), severe anomaly (2-3 times the threshold), and extremely severe anomaly (>3 times the threshold).

[0208] Fourth is contextual analysis, which combines the status of other related sensors within the same time period to verify the rationality of the anomaly. If multiple related sensors show anomalies at the same time, the credibility of the anomaly judgment is enhanced; if only a single sensor is abnormal while other sensors are normal, further verification is required.

[0209] For example, a downstream sensor had an anomaly score of 2.35 at 14:25, while the dynamic threshold was 1.73, exceeding the threshold by 35.8%, thus being classified as a moderate anomaly with a confidence level of 35.8%. This anomaly lasted for 25 minutes, meeting the continuity requirement, and two related upstream sensors also showed slight anomalies during the same period. Based on these combined findings, the anomaly was determined to be a genuine anomaly.

[0210] This embodiment constructs a high-precision and highly adaptable anomaly detection system by organically combining sliding window statistical analysis, dynamic threshold adjustment, and a multi-level judgment mechanism. The sliding window anomaly score calculation, through comprehensive evaluation of multi-dimensional statistical indicators, transforms point-like residual information into continuous anomaly scores, which not only improves the stability and noise resistance of detection but also accurately reflects the development trend and severity of anomalies. The dynamic threshold mechanism, based on the statistical characteristics and temporal correlation of historical data, achieves adaptive threshold adjustment, effectively addressing seasonal variations and diverse operating conditions in irrigation districts, and significantly reducing the false alarm and false negative rates of traditional fixed threshold methods.

[0211] In one embodiment of this invention, obtaining the dynamic anomaly judgment threshold includes the following steps:

[0212] S810. Based on the topological causal graph, obtain the spatial proximity relationship of the sensors, and construct the neighborhood influence weight matrix based on the spatial proximity relationship. The neighborhood influence weight is calculated according to the hydraulic connectivity and distance between the sensors.

[0213] S820: Obtain historical residual data of neighboring sensors of downstream sensors within the same time period, and calculate the weighted neighborhood residual benchmark based on the neighborhood influence weight matrix.

[0214] S830. Calculate the historical residual statistical characteristics of the downstream sensor, determine the historical residual distribution based on the historical residual statistical characteristics, and fit the historical residual distribution using a probability density function.

[0215] S840. Use the weighted neighborhood residual benchmark as the spatial reference value, and calculate the deviation between the spatial reference value and the center of the historical residual distribution.

[0216] S850. Spatial correction is performed on the historical residual distribution based on the deviation to obtain the spatially corrected residual distribution.

[0217] S860. Calculate the regional residual mean and standard deviation based on the spatially corrected residual distribution;

[0218] S870. Combine the regional residual mean and standard deviation to form a regional anomaly judgment benchmark;

[0219] S880. Based on the regional anomaly judgment benchmark and the preset confidence interval, the dynamic anomaly judgment threshold is calculated.

[0220] Among them, the dynamic anomaly judgment threshold is calculated based on the regional anomaly judgment benchmark and the preset confidence interval, including:

[0221] S881. Determine the central location and degree of variation of the residual distribution based on the regional anomaly judgment criteria;

[0222] S882. Determine the upper and lower boundary values ​​of the confidence interval, and calculate the critical value within the confidence interval. Use the critical value as the initial anomaly judgment threshold.

[0223] S883. The initial anomaly judgment threshold is dynamically adjusted based on the center location and the degree of variation to obtain the adjusted threshold.

[0224] S884. Use the adjusted threshold as the dynamic anomaly detection threshold.

[0225] The initial anomaly detection threshold is dynamically adjusted based on the center location and degree of variation to obtain the adjusted threshold, including:

[0226] S8831. Calculate the position adjustment factor based on the degree of offset of the center position, and calculate the variation adjustment factor based on the fluctuation range of the degree of variation.

[0227] S8832. Based on the data characteristics of the current monitoring period, determine the weight ratio of the location adjustment factor and the variation adjustment factor;

[0228] S8833. A linear combination method is used to weight and fuse the position adjustment factor and the variation adjustment factor with the initial anomaly judgment threshold.

[0229] S8834. Perform boundary constraint processing on the weighted fusion result to ensure that the adjusted threshold is within a reasonable range.

[0230] S8835. Use the result of boundary constraint processing as the adjusted threshold.

[0231] The spatial proximity relationships of sensors are obtained based on a topological causal graph, and a neighborhood influence weight matrix is ​​constructed. This process aims to quantify the degree of spatial influence between different sensors, providing a basis for subsequent regional anomaly detection. The server first extracts the spatial location information of the sensors from the topological causal graph, including attributes such as geographic coordinates, channel affiliation, and water flow path relationships.

[0232] The definition of spatial proximity considers not only geometric distance but, more importantly, hydraulic connectivity, i.e., whether there is a direct or indirect water flow connection between sensors. The server employs a multi-level proximity recognition algorithm. The first level is direct connectivity, i.e., sensor pairs directly connected by a single channel; the second level is indirect connectivity, i.e., sensor pairs connected by multiple channels but accessible by water flow; and the third level is regional proximity, i.e., sensor pairs that, although not directly hydraulically connected, are located within the same irrigation subsystem.

[0233] The calculation of neighborhood influence weights takes into account both hydraulic connectivity and geometric distance. ,in Let i be the influence weight between sensors i and j. For hydraulic connectivity weights, The distance weight is used for α, which is the balance coefficient (usually taken as 0.7). The hydraulic connectivity weight is calculated based on the complexity of the connectivity path: 1.0 for direct connectivity, 0.8 for first-order indirect connectivity, 0.6 for second-order indirect connectivity, and so on.

[0234] The distance weights use a Gaussian decay function: ,in Let σ represent the actual distance between the sensors, and σ be the attenuation parameter. For example, if a downstream sensor S005 is directly connected to an upstream sensor S003 at a distance of 800 meters, then the hydraulic connectivity weight is 1.0, the distance weight is 0.82, and the overall influence weight is 0.934. By calculating the influence weights of all sensor pairs, a complete neighborhood influence weight matrix is ​​constructed.

[0235] Acquiring historical residual data from neighboring sensors and calculating a weighted neighborhood residual benchmark is a crucial step in establishing a spatial reference standard. The server first identifies neighboring sensors that significantly influence the target downstream sensor based on a neighborhood influence weight matrix. Sensors with influence weights greater than 0.1 are typically selected as effective neighbors, forming a neighborhood sensor set. Then, residual data from these neighboring sensors within the same timeframe as the target sensor are extracted from the historical database, ensuring consistency of the time window and comparability of the data.

[0236] The selection of the historical time period needs to cover a sufficient sample size and representative operating conditions, typically including at least 3-6 months of normal operation data. For each neighboring sensor, the server calculates its historical residual sequence using the same method as the target sensor, including steps such as data preprocessing, standardization, and filtering.

[0237] The weighted neighborhood residual benchmark is calculated using a weighted average method: ,in As the weighted neighborhood residual benchmark, N is the set of neighborhood sensors. To influence the weight, Let be the historical residual mean of sensor j. To improve the stability of the benchmark, the server also calculates the variance and confidence interval of the neighborhood residuals: .

[0238] The calculation of historical residual statistical characteristics and probability density function fitting of downstream sensors are the core steps in establishing a quantitative anomaly judgment model. The server extracts the complete residual sequence of the target downstream sensor under normal operating conditions from the historical database. The data time span usually covers 12-24 months to ensure that it includes various seasonal changes and operating conditions.

[0239] The calculation of statistical characteristics involves multiple dimensions: descriptive statistics (mean, median, standard deviation, skewness, kurtosis), distribution shape parameters (quartiles, percentiles), and time-series characteristics (autocorrelation, periodicity), etc. Specifically, the mean reflects the central tendency of the residuals, the standard deviation reflects the dispersion of the data, and skewness... Reflects the symmetry of the distribution, kurtosis It reflects the sharpness of the distribution.

[0240] The probability density function is fitted using multiple candidate distribution models, including normal distribution, log-normal distribution, gamma distribution, Weibull distribution, etc. The parameters of each distribution are determined by the maximum likelihood estimation method, and then the goodness of fit is evaluated by methods such as the Anderson-Darling test to select the best-fit distribution.

[0241] Calculating the deviation between the spatial reference value and the historical residual distribution center is a step in achieving spatial correction. The deviation is calculated based on the weighted neighborhood residual benchmark and the historical residual distribution center obtained in the preceding steps. The distribution center is usually represented by the mean or median, and the selection criterion is based on the distribution characteristics of the data: the mean is used for symmetrical distributions, and the median is used for skewed distributions.

[0242] The formula for calculating the deviation is: Where Bias is the deviation. As a weighted neighborhood residual benchmark, This represents the center of the historical residual distribution. To assess the significance of the bias, the server also calculates the standardized value of the bias: ,in This represents the standard deviation of the historical residuals.

[0243] Server establishment deviation classification standard: minor deviation ( ), medium deviation ( Significant deviation ( ), serious deviation ( Different levels of deviation require different levels of correction strategies.

[0244] Spatial correction of historical residual distribution based on deviation is a technical step to eliminate systematic bias and improve regional consistency. The basic principle of spatial correction is to eliminate systematic bias through translation transformation while maintaining the data distribution shape and relative relationships. The correction transformation uses a linear translation model: ,in The corrected residual value, α is the original residual value, α is the correction strength coefficient, and Bias is the deviation calculated above.

[0245] Determining the correction strength coefficient requires striking a balance between bias elimination and data fidelity. The server employs an adaptive correction strength strategy: , where t is the correction threshold parameter (usually set to 1.5). Partial correction is used when the standardization deviation is small, and full correction is used when the deviation is large.

[0246] Calculating the regional residual mean and standard deviation based on the spatially corrected residual distribution is a step in establishing a unified anomaly detection benchmark. The server re-analyzes the spatially corrected residual data to calculate new distribution parameters. The formula for calculating the regional residual mean is: , where r'_i is the spatially corrected residual value.

[0247] The formula for calculating the regional standard deviation is: This parameter reflects the true variability of the residuals after eliminating systematic bias. To further improve the representativeness of regional parameters, the server also calculates weighted statistics that take into account the statistical characteristics of neighboring sensors.

[0248] The formation of regional anomaly identification criteria comprehensively considers the central tendency and variability characteristics of residuals, establishing a multi-level anomaly identification system. The regional anomaly identification criteria adopt a mean-standard deviation model: , where k is the standard deviation factor, determined based on the required detection sensitivity and false alarm tolerance.

[0249] The server establishes multi-level judgment benchmarks: early warning benchmark (k=1.5), anomaly benchmark (k=2.0), severe anomaly benchmark (k=2.5), and extreme anomaly benchmark (k=3.0). In addition to the standard deviation multiple method based on the normal distribution assumption, the server also uses a non-parametric method based on quantiles: slight anomaly benchmark (90th quantile), moderate anomaly benchmark (95th quantile), and severe anomaly benchmark (99th quantile).

[0250] The dynamic anomaly detection threshold is calculated based on regional anomaly detection criteria and pre-set confidence intervals. First, the central location and variability of the residual distribution are determined according to the regional anomaly detection criteria; the central location is determined using the regional mean. The degree of variation is expressed using regional standard deviation. Quantification. Then determine the boundaries of the confidence interval. Commonly used confidence levels include 90%, 95%, and 99%, with corresponding critical values ​​of 1.645σ, 1.96σ, and 2.576σ, respectively.

[0251] The initial anomaly detection threshold is calculated using a one-sided test. ,in This represents the critical value of the standard normal distribution corresponding to the confidence level. Dynamic adjustment based on the central location and degree of variability is a technique to improve the adaptability of the threshold.

[0252] First, calculate the position adjustment factor: ,in This represents the mean residual value for the current monitoring window. The variation adjustment factor reflects changes in data volatility. ,in This represents the standard deviation of the residuals for the current monitoring window.

[0253] The weighting ratios are determined based on the data characteristics of the current monitoring period. When the data center's location shifts significantly, the weight of the location adjustment factor is increased; when data volatility is significant, the weight of the variation adjustment factor is increased. The weighting formula is as follows: , .

[0254] Weighted fusion is performed using a linear combination method: The adjusted threshold is: .

[0255] Finally, boundary constraints are applied to set upper and lower boundaries for the threshold to prevent abnormal adjustments: , The final dynamic anomaly detection threshold is: .

[0256] This embodiment constructs a high-precision, highly adaptable threshold generation system by organically combining spatial correction, regional modeling, and adaptive adjustment. The quantification of spatial proximity relationships and the construction of the neighborhood influence weight matrix accurately reflect the physical connectivity characteristics of the sensor network. The calculation of the weighted neighborhood residual benchmark fully utilizes spatial correlation information, providing a regional reference standard. Probability density function fitting and statistical characteristic analysis establish a mathematical model, ensuring the statistical rationality of the threshold calculation. The spatial correction mechanism effectively eliminates the systematic bias of individual sensors, improving the consistency of detection standards within the region. The multi-level dynamic adjustment strategy enables the threshold to adapt to the current operating state and environmental conditions, significantly improving the accuracy and timeliness of detection.

[0257] In one embodiment of this invention, the method further includes the following steps:

[0258] S910. When there is an anomaly in the data of the downstream sensor, determine the anomaly type. The anomaly types include sensor malfunction anomaly, data transmission anomaly, and environmental change anomaly.

[0259] S920. Determine the scope of the anomaly's impact based on the anomaly type and the topological cause-effect graph;

[0260] S930: Generate an anomaly report based on the scope of the anomaly's impact and send the anomaly report to the monitoring terminal connected to the server.

[0261] After detecting anomalies in downstream sensor data, the server needs to further determine the specific type of anomaly. This is a crucial technical step in implementing anomaly response and fault handling. Anomaly type determination is based on multidimensional feature analysis and pattern recognition technology. The server has established a comprehensive anomaly classification system, primarily including three categories: sensor malfunction anomalies, data transmission anomalies, and environmental mutation anomalies.

[0262] Sensor malfunctions primarily refer to physical failures or performance degradation of the sensor hardware itself. Typical characteristics include data deviating from the normal range for extended periods, readings exhibiting a monotonous trend, abrupt changes in values, or complete cessation of updates. Data transmission anomalies refer to situations where the sensor is functioning normally, but problems occur during data transmission, manifesting as data packet loss, incorrect timestamps, repeated values, or transmission delays. Environmental mutation anomalies refer to situations where both the sensor and transmission system are functioning normally, but the monitoring environment undergoes changes beyond the normal range, such as extreme weather, human error, or equipment adjustments leading to drastic changes in measurement parameters.

[0263] The server uses a multi-feature fusion classification algorithm to identify anomaly types. The main features include: time features (duration of anomaly, frequency of occurrence, periodicity), numerical features (deviation amplitude, rate of change, distribution pattern), correlation features (relationship with adjacent sensors, similarity with historical patterns), and physical features (whether it conforms to physical laws, and the rationality of transmission delay).

[0264] For example, if a water level sensor displays the same value continuously for two hours, and this value significantly deviates from the expected physical value, while data from adjacent sensors shows normal changes, it is initially determined that the sensor is faulty. The server further classifies the data using machine learning algorithms such as decision trees or support vector machines. The decision rules are based on statistical analysis of historical failure cases and quantitative analysis of expert experience. The classification results include not only the anomaly type but also a confidence score and a ranked list of possible causes.

[0265] Determining the scope of anomaly impact based on anomaly type and topological cause-effect graph is a technical step in assessing the severity of a fault and developing response strategies. This process involves a systematic analysis based on the physical connectivity and operational relevance of the irrigation district's water conservancy system. The server first identifies downstream sensors and irrigation areas potentially affected by the anomaly based on the water flow path relationships in the topological cause-effect graph.

[0266] The scope of impact is determined using a multi-level diffusion model. The first level is the direct impact scope, which includes the direct downstream nodes of the abnormal sensors and the control equipment that relies on their data for decision-making. The second level is the indirect impact scope, which includes remote sensors and irrigation units that may be affected through multi-level transmission. The third level is the system impact scope, which includes related systems that may affect the overall scheduling strategy and operational efficiency.

[0267] Different types of anomalies exhibit different impact patterns and diffusion characteristics. Sensor malfunctions primarily affect monitoring and control systems that directly rely on the sensor's data; the degree of impact is directly proportional to the sensor's importance in the system. A single sensor failure can lead to scheduling anomalies across the entire irrigation zone. Data transmission anomalies have relatively minor impacts, mainly affecting data integrity and real-time performance. They do not directly affect physical water flow but may result in monitoring blind spots and decision-making delays.

[0268] The impact of sudden environmental changes is the most complex, requiring analysis based on the specific type of environmental change and physical transmission laws. For example, an upstream gate malfunction could affect the water level distribution of the entire downstream canal system within hours. The server uses graph theory algorithms to calculate the impact range, starting from the anomalous node and performing a breadth-first search along the directed edges of the topological causal graph, while incorporating spatiotemporal lag parameters to estimate the propagation time of the impact.

[0269] The degree of influence is quantified using a decay function model: ,in Let d be the initial impact intensity, t be the propagation distance, λ be the propagation time, and μ be the attenuation coefficients. For example, if a water level sensor malfunctions upstream of a branch canal, analysis reveals that it may affect three downstream water distribution points, 12 field gates, and irrigation scheduling covering 1200 hectares of farmland. The impact is expected to propagate to the furthest point within 1.5 hours, with the impact intensity decreasing to 15% of the initial value at the terminal.

[0270] Generating an anomaly report based on the scope of its impact and sending it to the monitoring terminal is the final output of the anomaly handling process. This process transforms complex technical analysis results into intuitive and actionable management information. The anomaly report adopts a structured format and contains multiple standardized information modules: the basic information module records basic data such as the time, location, sensor identification, and anomaly type; the technical analysis module describes in detail the anomaly's characteristic parameters, detection process, confidence level assessment, and classification results; the impact assessment module describes the anomaly's scope, severity, expected development trend, and potential risks; and the handling recommendation module provides specific response measures and priority ranking based on the anomaly type and scope of its impact.

[0271] The report generation employs a combination of template-based and intelligent methods. The server pre-sets report templates for different anomaly types, including standardized descriptive language and format specifications. Then, it automatically fills in relevant content based on the characteristic parameters of the specific anomaly. For example, the report template for sensor malfunction anomalies includes fixed modules such as a description of the malfunction phenomenon, possible causes of the malfunction, equipment maintenance suggestions, and temporary alternative solutions. The system will automatically generate personalized handling suggestions based on information such as the specific sensor model, service life, and maintenance records.

[0272] The report also includes visualization elements, such as map markers of anomaly locations, graphical representations of the affected area, and time-series charts of anomaly trends, enabling managers to intuitively understand the full picture of the anomaly. The server has established a multi-tiered report delivery mechanism, determining the notification recipients and priorities based on the severity and scope of the anomaly.

[0273] General anomalies are sent to on-duty personnel via standard email or system message. Important anomalies are notified to relevant personnel simultaneously via SMS and telephone. Emergency anomalies activate a multi-level linkage mechanism, simultaneously notifying all relevant personnel, including on-site operators, technical supervisors, and management leaders.

[0274] The server records the sending status and receipt confirmation of reports to ensure that abnormal information is delivered in a timely manner. Emergency reports that are not confirmed in a timely manner will be automatically escalated.

[0275] This embodiment constructs a complete anomaly response system through anomaly classification, impact analysis, and reporting mechanisms. The multi-dimensional feature fusion-based anomaly classification algorithm accurately identifies different types of anomalies, providing a basis for targeted handling and avoiding resource waste and time delays caused by blind processing. The impact range analysis based on topological causal graphs fully utilizes the physical connectivity of the system, accurately assessing the degree of impact and diffusion trend of anomalies, providing a reference for risk control and resource allocation. The structured anomaly reporting standardizes the content and format of information transmission, improving the accuracy and completeness of information delivery and reducing handling delays caused by information asymmetry.

[0276] This application embodiment also provides a server, including:

[0277] The memory is configured to store instructions; and

[0278] The processor is configured to retrieve instructions from memory and, when executing instructions, to implement the aforementioned method for analyzing irrigation district data anomalies.

[0279] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0280] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0281] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0282] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0283] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0284] Memory may include non-persistent memory in computer-readable media, such as random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM. Memory is an example of computer-readable media.

[0285] Computer-readable media includes both permanent and non-permanent, removable and non-removable media that can store information using any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transferable medium that can be used to store information accessible by a computing device. As defined herein, computer-readable media does not include transient computer-readable media, such as modulated data signals and carrier waves.

[0286] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0287] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A method for anomaly analysis of irrigation district data, characterized in that, Applied to a server, the method includes: Obtain irrigation district canal system topology data and sensor deployment information, and construct a topological causal graph based on the irrigation district canal system topology data and sensor deployment information to characterize the upstream and downstream water flow relationship between multiple sensors in the irrigation district. In the topological causal graph, the nodes are sensors and the directed edges are direct water flow paths between sensors. For sensor pairs connected by directed edges in a topological causal graph and having direct upstream and downstream relationships, time-series data of the sensor pairs under historical normal operating conditions are obtained, and spatiotemporal lag parameters are determined based on the time-series data. The spatiotemporal lag parameters are used to characterize the time required for upstream data changes to be transmitted to the downstream. For at least one downstream sensor, a physical conduction prediction model is constructed based on the historical normal time-series data of one or more upstream sensors and the corresponding spatiotemporal lag parameters between the downstream sensor and one or more upstream sensors. The physical conduction prediction model is used to predict the expected physical time-series data of the downstream sensor based on the time-series data of the upstream sensor. During the real-time monitoring phase, the actual measurement time series data of downstream sensors and the real-time time series data of one or more upstream sensors are acquired. Then, using the physical conduction prediction model, the expected physical time series data of downstream sensors are generated based on the real-time time series data and the spatiotemporal lag parameters. Based on the actual measurement time series data and the physical expected time series data of the downstream sensor, the residual sequence is calculated, and the presence of anomalies in the data of the downstream sensor is determined based on the residual sequence. Among them, determining the spatiotemporal lag parameters based on time series data includes: Extract the timing data of the upstream sensor and the timing data of the downstream sensor of the sensor pair; The correlation coefficient values ​​at different time offsets are obtained by calculating the time series data of upstream and downstream sensors based on the preset cross-correlation function. Determine the target time offset that maximizes the correlation coefficient value, and use the target time offset as the spatiotemporal lag parameter between sensor pairs.

2. The method according to claim 1, characterized in that, Before using the target time offset as a spatiotemporal hysteresis parameter between sensor pairs, the following is also included: Calculate the spatiotemporal lag parameter differences between adjacent sensor pairs within the same channel level to identify spatiotemporal lag parameters of local anomalies; The theoretical propagation time range is calculated based on the preset channel length and design flow rate, and the target time offset is limited to the theoretical propagation time range. When the target time offset exceeds a reasonable range, the short-term spatiotemporal lag parameter is determined based on the historical time series data within the recent time window, and the long-term spatiotemporal lag parameter is determined based on the historical time series data within the long-term time window, wherein the historical time series data within the long-term time window covers at least one complete irrigation cycle. The short-term and long-term lag parameters are weighted and fused to obtain the corrected spatiotemporal lag parameters.

3. The method according to claim 2, characterized in that, Before calculating the spatiotemporal hysteresis parameter differences between adjacent sensor pairs within the same channel level, the following is included: The target time offset is verified in a hierarchical manner, which includes lateral verification between sensors of the same level and longitudinal verification between sensors of different levels. The target time offset is subjected to hierarchical verification, including: The channel-level hierarchy of sensors is determined based on the topological cause-effect graph, and the sensors within the same channel level and the upstream and downstream sensors across channel levels are determined based on the channel-level hierarchy. Construct a cross-validation set based on sensors within the same channel level, and calculate the statistical distribution of spatiotemporal lag parameters of the sensors in the cross-validation set; A vertical verification set is constructed based on upstream and downstream sensors across channels, and the changing trend of the spatiotemporal lag parameters of the sensors in the vertical verification set is calculated. The target time offset is compared with the statistical distribution of the horizontal validation set and the changing trend of the vertical validation set to determine the rationality of the target time offset. Record the verification results and generate a tiered verification report.

4. The method according to claim 1, characterized in that, The physical conduction prediction model is a regression model. For at least one downstream sensor, based on historical normal time-series data from one or more upstream sensors and the corresponding spatiotemporal lag parameters between the downstream sensor and one or more upstream sensors, a physical conduction prediction model is constructed, including: The historical normal time series data of downstream sensors are used as the target variable, and the historical normal time series data of one or more upstream sensors after time alignment with spatiotemporal lag parameters are used as the feature variable. The regression model is trained based on the target variable and feature variables to determine the model parameters and obtain the trained physical conduction prediction model.

5. The method according to claim 1, characterized in that, Based on the actual measurement time-series data and the physical expected time-series data from the downstream sensors, the residual sequence is calculated, including: The actual measurement time series data and the physical expected time series data of the downstream sensor are time-aligned, and the difference between the actual measurement time series data and the physical expected time series data is calculated at each time point to obtain the original residual sequence; The original residual sequence is standardized to obtain the standardized residual sequence; A moving average filter is used to smooth the standardized residual sequence in order to remove high-frequency noise interference; The smoothed data is used as a residual sequence.

6. The method according to claim 1, characterized in that, Determining whether there are anomalies in the data from downstream sensors based on residual sequences includes: Based on the residual sequence, an anomaly score sequence is obtained by using a sliding window calculation. The dynamic anomaly judgment threshold is obtained, which is based on the historical residual sequence generated by the downstream sensor under historical normal operation conditions. The current abnormal score in the abnormal score sequence is compared with the dynamic abnormal judgment threshold. If the current abnormal score is greater than the dynamic abnormal judgment threshold, the data is determined to be abnormal.

7. The method according to claim 6, characterized in that, Obtain the dynamic anomaly detection threshold, including: Based on the topological causal graph, the spatial proximity relationship of the sensors is obtained, and a neighborhood influence weight matrix is ​​constructed based on the spatial proximity relationship. The neighborhood influence weight is calculated according to the hydraulic connectivity and distance between the sensors. Obtain historical residual data of neighboring sensors of downstream sensors within the same time period, and calculate the weighted neighborhood residual benchmark based on the neighborhood influence weight matrix; Calculate the historical residual statistical characteristics of downstream sensors, determine the historical residual distribution based on the historical residual statistical characteristics, and fit the historical residual distribution using a probability density function; The weighted neighborhood residual benchmark is used as a spatial reference value, and the deviation between the spatial reference value and the center of the historical residual distribution is calculated. Spatially correct the historical residual distribution based on the deviation to obtain the spatially corrected residual distribution. Calculate the regional residual mean and standard deviation based on the spatially corrected residual distribution; The regional residual mean and standard deviation are combined to form a regional anomaly judgment benchmark; Based on the regional anomaly judgment benchmark and the preset confidence interval, the dynamic anomaly judgment threshold is calculated.

8. The method according to claim 1, characterized in that, The method also includes: When there are anomalies in the data from downstream sensors, the anomaly type is determined. The anomaly types include sensor malfunction anomalies, data transmission anomalies, and environmental change anomalies. Determine the scope of the anomaly's impact based on the anomaly type and the topological cause-effect graph; An anomaly report is generated based on the scope of the anomaly's impact, and the anomaly report is sent to the monitoring terminal connected to the server.

9. A server, characterized in that, include: The memory is configured to store instructions; as well as The processor is configured to retrieve the instructions from the memory and, when executing the instructions, to implement the irrigation district data anomaly analysis method according to any one of claims 1 to 8.

Images