An anomaly detection method, system, device and medium based on global reachability density and adaptive threshold mechanism

Abstract

The application discloses an abnormality detection method based on global reachability density and an adaptive threshold mechanism, and comprises the following steps: collecting single-index historical water quality monitoring data of a target monitoring site, extracting a subsequence with a fixed length through a sliding window mode, constructing a near-neighbor index structure and extracting neighborhood information, obtaining the k near-neighbor set of each data point, respectively calculating the local reachability density of each data point in its local neighborhood, respectively calculating the global reachability density of each data point, statistically summarizing the GRD values of all data points in the sample set, calculating the mean value and the standard deviation, defining a dynamically adjusted distance threshold based on the coefficient of variation, further constructing an abnormality score function fusing the local and global density features, and realizing abnormal point determination according to the distance threshold. The application is suitable for complex time-series water quality monitoring data, can effectively improve the identification accuracy of weak abnormalities, local drifts and extreme mutations, and has strong universality and robustness.

Description

Technical Field

[0001] This invention relates to the field of water quality monitoring data anomaly analysis technology, and in particular to an anomaly detection method, system, device and medium based on global reachability density and adaptive threshold mechanism. Background Technology

[0002] The Local Outlier Factor (LOF) is a density-based anomaly detection algorithm specifically designed to identify data points whose local density is significantly lower than their neighbors (i.e., "local outliers"). While LOF has some applicability in anomaly detection tasks across various fields, it still suffers from the following significant shortcomings in practical applications such as automatic water quality monitoring in environmental engineering, severely limiting its practical value and widespread application:

[0003] (1) High false alarm rate in high-density areas, interfering with normal alarms:

[0004] Traditional LOF algorithms rely too heavily on local neighborhood density, which can easily lead to misjudgments in high-density data areas with large natural fluctuations. This can cause data points that should be considered normal to be identified as abnormal, resulting in a large number of invalid alarms and affecting the reliability and practicality of the monitoring system.

[0005] (2) Lack of a global density perspective makes it difficult to identify relative anomalies:

[0006] The LOF algorithm and its mainstream variants only rely on local reachability density for anomaly detection, ignoring the distribution background of data in the overall space. They cannot effectively characterize the relative position of samples in the global structure, resulting in poor stability of detection results in multi-scale data scenarios.

[0007] (3) The anomaly detection threshold is fixed, making it difficult to adapt to data heterogeneity:

[0008] Most existing methods use manually set or empirical parameters as anomaly judgment thresholds, lacking dynamic adjustment mechanisms. This leads to instability when facing heterogeneity in the water environment (such as dry season versus flood season, and differences between different watersheds), resulting in problems of being overly sensitive or overly insensitive.

[0009] (4) Slow response to systemic trend changes, high risk of missed detection:

[0010] Environmental data often contains non-sudden but trend-based abnormal patterns, such as a slow decline or long-term fluctuation in water quality indicators. Traditional LOF algorithms are weak in responding to these types of changes, and are prone to missing detections or mistaking them for normal evolution processes, affecting the timeliness of event warnings.

[0011] (5) High computational complexity limits real-time monitoring applications:

[0012] In high-frequency sampling and multi-parameter parallel water quality automatic monitoring systems, existing algorithms have a heavy computational burden in steps such as nearest neighbor search and density estimation. Especially when processing millions of monitoring records, the response speed cannot meet the engineering requirements for real-time performance and scalability.

[0013] (6) The ability to identify anomalies decreases when anomalies are clustered:

[0014] When anomalies form a dense distribution within a specific spatiotemporal range (such as pollution clumps or sewage discharge events), their internal density is relatively high. Existing algorithms are prone to misclassifying them as "normal clusters," causing the overall anomaly clusters to be ignored, which seriously affects the accurate identification of pollution events.

[0015] Therefore, it is urgent to solve the above problems. Summary of the Invention

[0016] Purpose of the invention: The first objective of this invention is to provide an anomaly detection method based on global reachability density and adaptive threshold mechanism, which is applicable to complex time-series water quality monitoring data and can effectively improve the accuracy of identifying weak anomalies, local drift and extreme mutations. It has strong versatility and robustness.

[0017] The second objective of this invention is to provide an anomaly detection system based on global reachability density and an adaptive threshold mechanism.

[0018] A third objective of this invention is to provide an electronic device.

[0019] A fourth objective of this invention is to provide a computer-readable storage medium.

[0020] Technical Solution: To achieve the above objectives, this invention discloses an anomaly detection method based on global reachability density and an adaptive threshold mechanism, comprising the following steps:

[0021] (1) Collect historical water quality monitoring data of a single indicator from the target monitoring stations, extract fixed-length subsequence samples using the sliding window method to form a sample set, construct a nearest neighbor index structure, and calculate the value of each data point. p of k Nearest neighbor set This allows for the extraction of neighborhood information;

[0022] (2) Calculate each data point p Local reachability density in its local neighborhood ;

[0023] (3) Within the sample set, based on k Calculate the reachability distance of the farthest point among the nearest neighbors for each data point. p Global reachability density ;

[0024] (4) For all data points in the sample set GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated. The standard deviation is then used as the criterion for further analysis. with the mean The ratio is defined as the coefficient of variation. The dynamic distance threshold is determined by combining the adjustment coefficient α. ;

[0025] ,

[0026] (5) Based on local reachability density and global reachability density, construct a nearest neighbor fusion anomaly score. GRD-LOF and with dynamic distance threshold Comparisons are made to enable automatic identification and judgment of anomalies;

[0027] (6) Introduce an incremental update mechanism. After each round of new data points arrive, only the newly added points and their associated information are partially updated.

[0028] Optionally, step (1) specifically includes the following steps:

[0029] (1.1) Obtain continuous historical monitoring data of a certain water quality indicator from the target monitoring station, construct a one-dimensional time series, denoted as . Where T represents the total number of observation time steps, X For this water quality index to be in continuous T Historical water quality monitoring sequences at each time step , indicating a point in time t The monitored values;

[0030] (1.2) For the sequence X Use a length of w The sliding window extracts subsequences to form a sample set: ,in ,in, w For window length, The total number of subsequences, and each sample subsequence Indicates continuity w The indicator monitoring value at each moment;

[0031] (1.3) For each sample subsequence Perform z-score standardization to obtain normalized sample subsequences. ,in, For sample subsequence The mean, For sample subsequence Standard deviation of the normalized sample This constitutes the normalized sample set. ;

[0032] (1.4) For the sample set Establish a spatial index structure;

[0033] (1.5) For the sample set Each normalized subsequence Using the constructed spatial index structure, a k-nearest neighbor query is performed based on Euclidean distance to obtain its k most similar neighbor subsequences: .

[0034] Optionally, step (2) specifically includes the following steps:

[0035] (2.1) For the normalized sample set any subsequence in , and its first j A neighbor The reachability distance between them is defined as:

[0036] ,

[0037] Sample points Compared with its own k The Euclidean distance between the nearest neighbors, representing the lower bound of its local density, is defined as:

[0038] ,

[0039] in Represents sample points of k A set of nearest neighbors; For sample points its neighbors q The Euclidean distance between them is defined as: ,in Represents sample points In the t The value of each component For the nearest neighbor q in the th... t The value of each component;

[0040] The Euclidean distance between two sample points is defined as:

[0041] ,

[0042] in, Represents sample points The t One portion, Represents sample points Thet One component;

[0043] (2.2) Calculate the local reachability density LRD , after defining arbitrary sample points After calculating the reachability distance between each of its neighbors, compute the local reachability density of any subsequence. It is the reciprocal of the average reachable distance within its neighborhood, calculated using the following formula:

[0044] ,

[0045] The denominator represents the sample points. To its k The average reachable distance of the nearest neighbor.

[0046] Optionally, step (3) specifically includes the following steps:

[0047] For each normalized sample subsequence Its global reachability density Defined as the reciprocal of the average reachability distance between itself and all other sample points except itself, the formula is:

[0048] ,

[0049] Where n is the total number of samples. To remove Other samples besides; For point With point The reachable distance between them is defined as follows:

[0050] ,

[0051] For sample points To its own k The Euclidean distance between the nearest neighbors is calculated using the following formula:

[0052] ,

[0053] in s It is the farthest nearest neighbor sample point;

[0054] For sample points With sample points The Euclidean distance between them is calculated using the following formula:

[0055] ,

[0056] in Represents sample points Thet One portion, Represents sample points The t One portion, w This is the length of the sliding window.

[0057] Optionally, step (5) specifically includes the following steps:

[0058] (5.1) Calculate the nearest neighbor fusion anomaly score for each data point. GRD-LOF, The calculation formula is:

[0059] ,

[0060] ,

[0061] in, For point of k Nearest neighbor set; Point for neighbors Global reachability density; For point Local reachability density;

[0062] (5.2) When GRD-LOF The scores are not normally distributed; the quantile method is used to analyze all points. GRD-LOF Scores are sorted in ascending order, and a high percentile value is set as an outlier threshold. If a certain data point... GRD-LOF If the score exceeds this percentile value, it is judged as an outlier;

[0063] when GRD-LOF When the scores exhibit a "centralized distribution with small fluctuations but extreme high values," the standard deviation discriminant method is used to statistically analyze all points. GRD-LOF The mean μ and standard deviation σ of the scores are used to set a threshold for discrimination. If a certain data point... GRD-LOF If the value exceeds the discrimination threshold, it is determined to be an outlier. The discrimination threshold is calculated as follows:

[0064] ,

[0065] in β For dynamic adjustment coefficients;

[0066] when GRD-LOF If the score does not meet the above two characteristics, the empirical multiplier method is used. If a certain data point GRD-LOF Score exceeded λ Distance threshold times If it is an outlier, then it is considered an anomaly. λ This is the experience multiplier.

[0067] Optionally, step (6) specifically includes the following steps:

[0068] (6.1) Sample sliding window update

[0069] A fixed-length sliding window mechanism is used to maintain the sample set; at each new time step t+1, a normalized subsequence is extracted from the historical monitoring sequence. Add to the sample set, and simultaneously remove the earliest one. To maintain a constant total sample size of n, a local update is performed on the index structure used for nearest neighbor search, i.e., new sample points are inserted. Delete the excluded points This avoids rebuilding the entire index structure;

[0070] (6.2) Incremental calculation of LRD and GRD

[0071] For the newly added subsequence Re-execute k Nearest neighbor search to determine its nearest neighbor set And calculate its local reachability density based on reachability distance. and global reachability density To maintain the statistical characteristics of the GRD distribution in real time, the following weighted formula is used to calculate the mean of the GRD. and standard deviation Make incremental adjustments:

[0072] ,

[0073] ,

[0074] ,

[0075] in, For newly added sample points GRD value, Points to be removed GRD value, n This represents the total number of samples within the sliding window; This is the variation in standard deviation;

[0076] (6.3) Linked update of anomaly score and judgment threshold

[0077] For new sample points According to its GRD and LRD Calculate its anomaly score GRD-LOF The abnormal score sequence is then incorporated into the abnormal score sequence. Based on the updated abnormal score sequence, its quantiles, mean, and standard deviation are dynamically maintained, and the judgment threshold is adjusted in real time.

[0078] Based on the same inventive concept, the present invention provides an anomaly detection system based on global reachability density and an adaptive threshold mechanism, comprising:

[0079] The neighborhood information extraction module is used to collect historical water quality monitoring data for a single indicator from target monitoring stations. It extracts fixed-length sub-sequence samples using a sliding window method to form a sample set, constructs a nearest neighbor index structure, and calculates the nearest neighbor index for each data point. p of k Nearest neighbor set This allows for the extraction of neighborhood information;

[0080] The local reachability density calculation module is used to calculate the density of each data point. p Local reachability density in its local neighborhood ;

[0081] The global reachability density calculation module is used to calculate the reachability density within the sample set based on... k Calculate the reachability distance of the farthest point among the nearest neighbors for each data point. p Global reachability density ;

[0082] The dynamic distance threshold determination module is used to determine the distance threshold for all data points in the sample set. GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated. The standard deviation is then used as the criterion for further analysis. with the mean The ratio is defined as the coefficient of variation. The dynamic distance threshold is determined by combining the adjustment coefficient α. ;

[0083] ,

[0084] The anomaly identification and judgment module is used to construct a nearest-neighbor fusion-type anomaly score based on local and global reachability densities. GRD-LOF and with dynamic distance threshold Comparisons are made to enable automatic identification and judgment of anomalies;

[0085] The incremental update module is used to introduce an incremental update mechanism, which only updates the newly added points and their associated information after each round of new data points arrive.

[0086] Optionally, the neighborhood information extraction module obtains continuous historical monitoring data of a certain water quality indicator from the target monitoring station, constructs a one-dimensional time series, denoted as... Where T represents the total number of observation time steps, X For this water quality index to be in continuous T Historical water quality monitoring sequences at each time step , indicating a point in timet The monitored values;

[0087] For sequence X Use a length of w The sliding window extracts subsequences to form a sample set: ,in ,in, w For window length, The total number of subsequences, and each sample subsequence Indicates continuity w The indicator monitoring value at each moment;

[0088] For each sample subsequence Perform z-score standardization to obtain the normalized sample subsequences: ,in, For sample subsequence The mean, For sample subsequence Standard deviation of the normalized sample This constitutes the normalized sample set. ;

[0089] For sample set Establish a spatial index structure;

[0090] For sample set Each normalized subsequence Using the constructed spatial index structure, a k-nearest neighbor query is performed based on Euclidean distance to obtain its k most similar neighbor subsequences: .

[0091] Based on the same inventive concept, the present invention provides an electronic device comprising one or more processors, one or more memories, and one or more programs, wherein the programs are stored in the memory and configured to be executed by the processor, and when the programs are loaded onto the processor, they implement the steps of an anomaly detection method based on a global reachability density and adaptive threshold mechanism as described above.

[0092] Based on the same inventive concept, the present invention provides a computer-readable storage medium storing a computer program, the computer program including program instructions, which, when executed by a processor, cause the processor to perform steps according to an anomaly detection method based on a global reachability density and adaptive threshold mechanism as described above.

[0093] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages:

[0094] (1) This invention introduces the Global Reachability Density (GRD) index as a global reference benchmark for anomaly detection, thereby enhancing the robustness of density assessment in the traditional LOF algorithm and reducing the false alarm rate in high-density areas.

[0095] (2) By fusing local density and global density information, this invention enhances the perception of global data structure, enabling it to not only focus on the local differences around data points, but also identify their relative degree of anomaly in the entire data space, thereby improving the detection stability in multi-scale density environments.

[0096] (3) The present invention utilizes the coefficient of variation (CV) to dynamically adjust the distance threshold required for anomaly determination, thereby achieving adaptive adjustment of the anomaly detection threshold, avoiding the limitations of manually setting fixed parameters, and enhancing the generalization ability of the model under different water bodies, time periods and environmental conditions.

[0097] (4) By introducing global quantile constraints on reachability, this invention enhances the model’s ability to perceive slow-changing or trend-type anomalies (such as continuous rise / fall of water quality indicators), improves the sensitivity to the identification of systemic trend changes, reduces misjudgment of normal evolution processes, and improves the response to pollution evolution processes.

[0098] (5) This invention adopts efficient nearest neighbor query structures such as KD-Tree or Ball-Tree, combined with quantile statistics and sliding window mechanism, to optimize the algorithm calculation efficiency, adapt to the needs of large-scale high-frequency data processing, realize low latency and high throughput online anomaly detection, and meet the requirements of real-time water quality monitoring system for response speed.

[0099] (6) This invention, through z-score normalization and density-based unsupervised logic, is applicable to water quality time series analysis under different indicators and scales, and has excellent scale adaptability and transferability.

[0100] (7) This invention addresses the problem of density reversal caused by local abnormal clusters. By using global density constraints and abnormal score distribution correction mechanisms, it effectively identifies abnormal clusters, improves the ability to identify dense abnormal points, avoids ignoring abnormal clusters, and enhances the algorithm's ability to accurately identify contaminated clumps or sudden events. Attached Figure Description

[0101] Figure 1 This is a schematic diagram of the process of the present invention;

[0102] Figure 2 This is a comparative diagram showing the CHLA outlier identification performed by the present invention and the algorithm before the improvement;

[0103] Figure 3 This is a comparative diagram showing the COND outlier identification performed by the present invention and the algorithm before the improvement.

[0104] Figure 4 This is a comparative diagram showing the CODCR anomaly identification performed by the present invention and the algorithm before the improvement.

[0105] Figure 5 This is a system framework diagram of the present invention. Detailed Implementation

[0106] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0107] It should be understood that the invention can be embodied in various forms and should not be construed as being limited to the embodiments set forth herein. Rather, providing these embodiments will make the disclosure thorough and complete, and will fully convey the scope of the invention to those skilled in the art. In the drawings, for clarity, the dimensions and relative dimensions of the parts may be exaggerated. The same reference numerals denote the same parts throughout.

[0108] Example 1: As Figure 1 As shown, this invention discloses an anomaly detection method based on global reachability density and an adaptive threshold mechanism, comprising the following steps:

[0109] (1) Collect historical water quality monitoring data of a single indicator from the target monitoring stations, extract fixed-length subsequence samples using the sliding window method to form a sample set, construct a nearest neighbor index structure, and calculate the value of each data point. p of k Nearest neighbor set This allows for the extraction of neighborhood information;

[0110] (1.1) Obtain continuous historical water quality monitoring data for a specific water quality indicator (such as dissolved oxygen, total phosphorus, ammonia nitrogen, etc.) from the target monitoring station, and construct a one-dimensional time series, denoted as: Where T represents the total number of observation time steps, X For this water quality index to be in continuous T Historical water quality monitoring sequences at each time step , indicating a point in time t The monitored values;

[0111] (1.2) For the sequence X Use a length of w The sliding window extracts subsequences to form a sample set: ,in , w For window length, For the number of samples, each sample subsequence Indicates continuityw The indicator monitoring value at each moment;

[0112] (1.3) To eliminate amplitude bias and scaling effects between different subsequences, for each sample subsequence z-score standardization is performed to obtain normalized sample subsequences. ,in For sample subsequence The mean, For sample subsequence Standard deviation of the normalized sample This constitutes the normalized sample set. ;

[0113] (1.4) For the sample set Establish a spatial index structure: to accelerate subsequent... k The nearest neighbor search process, for the sample set Establish a spatial index structure: if the window length w ≤ 20, KD-Tree is preferred; if the window length w > 20, Ball-Tree is preferred.

[0114] The index is constructed as follows:

[0115] KD-Tree: Based on dimensional rotation and median partitioning, samples are recursively divided into left and right subtrees to construct a balanced binary tree structure;

[0116] Ball-Tree: Using the two data points that are farthest apart in space as the center of the sphere, the data is divided into multiple spherical subspaces to construct a hierarchical spherical structure;

[0117] (1.5) Calculation k Nearest neighbor set : For the sample set Each normalized subsequence Determined based on spatial index structure k Objects, based on Euclidean distance k Nearest neighbor search to find the most similar k Neighbor subsequences: .

[0118] (2) For each data point p Calculate the data point p achievable density in its local neighborhood ;

[0119] (2.1) Define Reachability Distance: For a normalized sample set any subsequence in , and its first jA neighbor The reachability distance between them is defined as:

[0120] ,

[0121] Sample points Compared with its own k The Euclidean distance between the nearest neighbors, representing the lower bound of its local density, is defined as:

[0122] ,

[0123] in Represents sample points of k A set of nearest neighbors; For sample points its neighbors q The Euclidean distance between them is defined as: ,in Represents sample points In the t The value of each component For the nearest neighbor q in the th... t The value of each component;

[0124] The Euclidean distance between two sample points is defined as:

[0125] ,

[0126] in, Represents sample points The t One portion, Represents sample points The t The smaller the distance between the components, the more similar the patterns.

[0127] (2.2) Calculate the local reachability density LRD (Local Reachability Density), when any sample point is defined After calculating the reachability distance between each of its neighbors, compute the local reachability density of any subsequence. It is the reciprocal of the average reachable distance within its neighborhood, calculated using the following formula:

[0128] ,

[0129] The denominator represents the sample points. Average reachable distance to its k nearest neighbors The larger the value, the denser the area where the point is located; conversely, the smaller the value, the sparser the area.

[0130] (3) Introduce global reachability density GRD Global Reachability Density, as a complement to traditional local density, provides a stable global density benchmark for anomaly detection by characterizing the reachability level of samples throughout the entire data space.

[0131] Within the sample set, based on k Calculate the reachability distance of the farthest point among the nearest neighbors for each data point. p Global reachability density ;

[0132] For each normalized sample subsequence Its global reachability density Defined as the reciprocal of the average reachability distance between it and all other sample points except itself, the formula is as follows:

[0133] ,

[0134] in n This represents the total number of samples. To remove Other samples besides; For point With point The reachable distance between them is defined as follows:

[0135] ,

[0136] For sample points To its own k The Euclidean distance between the nearest neighbors is calculated using the following formula:

[0137]

[0138] in s It is the farthest nearest neighbor sample point;

[0139] For sample points With sample points The Euclidean distance is calculated using the following formula:

[0140] ,

[0141] Where represents the sample point The t One portion, Represents sample points The t One portion,w This is the length of the sliding window.

[0142] (4) For all data points in the sample set GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated. The standard deviation is then used as the criterion for further analysis. with the mean The ratio is defined as the coefficient of variation. The dynamic distance threshold is determined by combining the adjustment coefficient α. To avoid false alarms or missed alarms under fixed thresholds in different density regions;

[0143] (4.1) For all data points in the sample set GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated:

[0144] ,

[0145]

[0146] in n This represents the total number of samples. For the first i sample points Global reachability density;

[0147] (4.2) Calculate the coefficient of variation:

[0148] ,

[0149] (4.3) Define dynamic distance threshold :

[0150]

[0151] in For adjustment coefficients, The value ranges from 1 to 2, and is used to control the model sensitivity; when When the size is large, the model becomes more conservative, identifying only strong outliers. When the size is smaller, the model is more sensitive and can identify weak boundary anomalies.

[0152] (5) Based on local reachability density and global reachability density, construct a nearest neighbor fusion anomaly score. GRD-LOF and with dynamic distance threshold Comparisons are made to enable automatic identification and judgment of anomalies;

[0153] (5.1) Calculate the nearest neighbor fusion anomaly score GRD-LOF for each data point. The calculation formula is as follows:

[0154] ,

[0155] ,

[0156] in, For point of k Nearest neighbor set; Point for neighbors Global reachability density, For point Local reachability density;

[0157] (5.2) To adapt to the distribution characteristics of different types of data, GRD-LOF The anomaly detection mechanism introduces three strategies, each applicable to the following three types of scenarios;

[0158] (5.2.1) When GRD-LOF The scores are not normally distributed; the quantile method is used to analyze all points. GRD-LOF Scores are sorted in ascending order, and a high percentile value is set as an outlier threshold. If a certain data point... GRD-LOF If the score exceeds this percentile value, it is judged as an outlier;

[0159] Specifically, for GRD-LOF A normality test was performed on the scores. If the p-value for the Shapiro–Wilk normality test was <0.01, the null hypothesis was rejected, indicating that... GRD-LOF When the scores are not normally distributed, the quantile method is preferred; specifically, all points... GRD-LOF The scores are sorted, and a certain high percentile value (such as 95%, 97.5%, or 99%) is set as the outlier threshold. Scores exceeding this threshold are considered outliers, forming an outlier set.

[0160] ,

[0161] (5.2.2) When GRD-LOF When the scores exhibit a "centralized distribution with small fluctuations but extreme high values," the standard deviation discriminant method is used to statistically analyze all points. GRD-LOF The mean μ and standard deviation σ of the scores are used to set a threshold for discrimination. If a certain data point... GRD-LOF If the value exceeds the discrimination threshold, it is determined to be an outlier. The discrimination threshold is calculated as follows:

[0162] ,

[0163] in β For dynamic adjustment of the coefficient, it is recommended to take a value of 2 to 3, which is suitable for scenarios where abnormal scores are approximately normally distributed;

[0164] When the GRD-LOF score shows a "centralized distribution with small fluctuations but extreme high values", it is a typical low variance and heavy-tailed distribution or a prominent outlier. Specifically, the mean and median are very close, indicating that most of the data is concentrated; the standard deviation is significantly small, indicating that the overall data fluctuation is low; the maximum value (max) is significantly higher than the upper quartile (Q3 + 1.5 × IQR). When there are extreme high values ​​on the surface, the standard deviation method is preferred.

[0165] Q3 refers to the third quartile (75th quartile), which is the larger 25% dividing point in the data. IQR represents the interquartile range, which indicates the distribution range of the middle 50% of the data. IQR = Q3 - Q1, where Q1 refers to the first quartile (25th quartile), which is the smaller 25% dividing point in the data. IQR represents the interquartile range.

[0166] (5.2.3) When GRD-LOF If the score does not meet the above two characteristics, the empirical multiplier method is used. If a certain data point GRD-LOF Score exceeded λ Distance threshold times If it is an outlier, then it is considered an anomaly. λ For experience multipliers, λ A value of 1.5 to 2.0 can be used to control the screening sensitivity; λ =1.5: Slightly sensitive, suitable for scenarios where anomalies are obvious or where the requirements for anomaly identification are high (such as water quality early warning). λ =2.0: Medium sensitivity, suitable for most routine monitoring or environments with relatively stable data distribution.

[0167] (6) To meet the real-time processing requirements of high-frequency flow cytometry data in the automatic water quality monitoring system, an incremental update mechanism is introduced to realize density calculation indicators (such as...). GRD , LRD , GRD-LOF Dynamic maintenance of scores and distance thresholds; specifically, after each round of new data points arrive, only the newly added points and their associated information are locally updated to avoid reconstructing the entire dataset;

[0168] (6.1) Sample sliding window update

[0169] A fixed-length (e.g., 24-hour) sliding window mechanism is used to maintain the sample set; at each new time step t+1, a normalized subsequence is extracted from the historical monitoring sequence. Add to the sample set, and simultaneously remove the earliest one. To maintain a constant total sample size of n, a local update (such as KD-Tree / Ball-Tree) is performed on the index structure used for nearest neighbor search, i.e., new sample points are inserted. Delete the excluded points This avoids rebuilding the entire index structure;

[0170] (6.2) Incremental updates of LRD and GRD

[0171] For the newly added subsequence Re-execute k Nearest neighbor search to determine its nearest neighbor set And calculate its local reachability density based on reachability distance. and global reachability density To maintain the statistical characteristics of the GRD distribution in real time, the following weighted formula is used to calculate the mean of the GRD. and standard deviation Make incremental adjustments:

[0172] ,

[0173] ,

[0174]

[0175] in For new members GRD value, Points to be removed GRD value, n This represents the number of samples within the window. This is the standard deviation variation term, which can be quickly updated based on approximate derivation;

[0176] (6.3) Anomaly scores are updated in conjunction with thresholds.

[0177] For new sample points According to its GRD and LRD Instant score calculation GRD-LOF The updated anomaly score sequence is then incorporated into the anomaly score sequence. Based on this updated sequence, the quantiles, mean, and standard deviation are dynamically maintained. The judgment threshold is then dynamically adjusted in real time using preset rules (such as quantile method, standard deviation discrimination method, or empirical multiple method) to achieve highly responsive anomaly detection decisions. The results can be saved as .csv, .json, or database tables for later use.

[0178] like Figure 2 , Figure 3 and Figure 4 As shown, the gray lines represent the original data, the black dots represent non-anomalies, the yellow dots represent anomalies identified by the algorithm before the improvement, the red dots represent anomalies identified by the present invention, and the blue dots represent anomalies that can be identified by both the algorithm before the improvement and the present invention.

[0179] like Figure 2 As shown, the CHLA index sequence from station 71553 in October 2019 was selected as a sample. The original data generally exhibited high-frequency oscillations, with localized short-term spikes. Using the anomaly detection method proposed in this invention, multiple anomalies with instantaneous spikes (marked in red) were successfully identified, while maintaining high accuracy in identifying key extreme points. In contrast, the traditional LOF method misidentified several high-frequency fluctuations (marked in yellow), leading to a decrease in anomaly identification accuracy. This example verifies the robustness and confidence level of this invention in handling abrupt change detection tasks in highly volatile, non-stationary time series.

[0180] like Figure 3 As shown, an analysis of the COND index sequence for site 71553 in September 2018 reveals an overall decreasing trend in the raw data, accompanied by several short-term sharp declines. The anomaly detection method proposed in this invention effectively identifies anomalous points (marked in red) during the trend evolution process and accurately detects abnormally high points (marked in blue) at key trend turning points, demonstrating adaptability to trend-type anomalies and multi-source fusion identification capabilities. In contrast, the traditional LOF method exhibits false identifications in multiple locally stable intervals (marked in yellow), demonstrating insufficient adaptability to trend structures. Therefore, this invention is more suitable for anomaly detection scenarios in dynamic sequences containing trend characteristics.

[0181] like Figure 4 As shown, using the CODCR data from station eh202402 in January 2025 as the research object, the original sequence showed relatively small overall fluctuations, but multiple isolated mutation points appeared in the mid-to-late stages. The anomaly detection method proposed in this invention can accurately identify local deviation points (marked in red) in the mid-term slight increase phase and has good sensitivity and recognition ability for late-term surge peaks (marked in blue). In contrast, the traditional LOF method misjudges anomalies in many normal fluctuations (marked in yellow), affecting detection accuracy. The results show that this invention is suitable for detecting nested outliers in a low-fluctuation background, possessing superior anomaly boundary determination ability and generalization adaptability.

[0182] Traditional LOF algorithms exhibit high false alarm rates for indicators such as TP and CHLA, easily misclassifying high or low values ​​within normal fluctuations as abnormal. This false alarm phenomenon leads to unnecessary frequent alarms in practical water quality monitoring applications, reducing system reliability. This invention, by introducing global reachability density constraints and dynamic threshold adjustment, significantly suppresses misjudgments of normal fluctuations. Particularly in the CHLA indicator, GRD-LOF can distinguish between normal concentration fluctuations and abrupt abnormal changes, reducing the probability of incorrect detection of data within normal ranges and effectively improving the practicality and stability of the detection system.

[0183] Example 2:Figure 5 As shown, an anomaly detection system based on global reachability density and adaptive threshold mechanism in this invention includes:

[0184] The neighborhood information extraction module is used to collect historical water quality monitoring data for a single indicator from target monitoring stations. It extracts fixed-length sub-sequence samples using a sliding window method to form a sample set, constructs a nearest neighbor index structure, and calculates the nearest neighbor index for each data point. p of k Nearest neighbor set This allows for the extraction of neighborhood information;

[0185] The neighborhood information extraction module obtains continuous historical monitoring data of a certain water quality indicator from the target monitoring station, constructs a one-dimensional time series, and denoted as... Where T represents the total number of observation time steps, X For this water quality index to be in continuous T Historical water quality monitoring sequences at each time step , indicating a point in time t The monitored values;

[0186] For sequence X Use a length of w The sliding window extracts subsequences to form a sample set: ,in ,in, w For window length, The total number of subsequences, and each sample subsequence Indicates continuity w The indicator monitoring value at each moment;

[0187] For each sample subsequence Perform z-score standardization to obtain the normalized sample subsequences: ,in, For sample subsequence The mean, For sample subsequence Standard deviation of the normalized sample This constitutes the normalized sample set. ;

[0188] For sample set Establish a spatial index structure;

[0189] For sample set Each normalized subsequence Using the constructed spatial index structure, a k-nearest neighbor query is performed based on Euclidean distance to obtain its k most similar neighbor subsequences: .

[0190] The local reachability density calculation module is used to calculate the density of each data point. p Local reachability density in its local neighborhood ;

[0191] In the local reachability density calculation module, for the normalized sample set any subsequence in , and its first j A neighbor The reachability distance between them is defined as:

[0192] ,

[0193] Sample points Compared with its own k The Euclidean distance between the nearest neighbors, representing the lower bound of its local density, is defined as:

[0194] ,

[0195] in Represents sample points of k A set of nearest neighbors; For sample points its neighbors q The Euclidean distance between them is defined as: ,in Represents sample points In the t The value of each component For the nearest neighbor q in the th... t The value of each component;

[0196] The Euclidean distance between two sample points is defined as:

[0197] ,

[0198] in, Represents sample points The t One portion, Represents sample points The t One component;

[0199] Calculate local reachability density LRD , after defining arbitrary sample points After calculating the reachability distance between each of its neighbors, compute the local reachability density of any subsequence. It is the reciprocal of the average reachable distance within its neighborhood, calculated using the following formula:

[0200] ,

[0201] The denominator represents the sample points. To its k The average reachable distance of the nearest neighbor.

[0202] The global reachability density calculation module is used to calculate the reachability density within the sample set based on... k Calculate the reachability distance of the farthest point among the nearest neighbors for each data point. p Global reachability density ;

[0203] The global reachability density calculation module calculates each normalized sample subsequence. Its global reachability density Defined as the reciprocal of the average reachability distance between itself and all other sample points except itself, the formula is:

[0204] ,

[0205] Where n is the total number of samples. To remove Other samples besides; For point With point The reachable distance between them is defined as follows:

[0206] ,

[0207] For sample points To its own k The Euclidean distance between the nearest neighbors is calculated using the following formula:

[0208] ,

[0209] in s It is the farthest nearest neighbor sample point;

[0210] For sample points With sample points The Euclidean distance between them is calculated using the following formula:

[0211] ,

[0212] in Represents sample points The t One portion, Represents sample points The t One portion, w This is the length of the sliding window.

[0213] The dynamic distance threshold determination module is used to determine the distance threshold for all data points in the sample set. GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated. The standard deviation is then used as the criterion for further analysis. with the mean The ratio is defined as the coefficient of variation. The dynamic distance threshold is determined by combining the adjustment coefficient α. ;

[0214] .

[0215] The anomaly identification and judgment module is used to construct a nearest-neighbor fusion-type anomaly score based on local and global reachability densities. GRD-LOF and with dynamic distance threshold Comparisons are made to enable automatic identification and judgment of anomalies;

[0216] The anomaly detection and judgment module calculates the nearest neighbor fusion anomaly score for each data point. GRD-LOF, The calculation formula is:

[0217] ,

[0218] ,

[0219] in, For point of k Nearest neighbor set; Point for neighbors Global reachability density; For point Local reachability density;

[0220] when GRD-LOF The scores are not normally distributed; the quantile method is used to analyze all points. GRD-LOF Scores are sorted in ascending order, and a high percentile value is set as an outlier threshold. If a certain data point... GRD-LOF If the score exceeds this percentile value, it is judged as an outlier;

[0221] when GRD-LOF When the scores exhibit a "centralized distribution with small fluctuations but extreme high values," the standard deviation discriminant method is used to statistically analyze all points. GRD-LOF The mean μ and standard deviation σ of the scores are used to set a threshold for discrimination. If a certain data point... GRD-LOF If the value exceeds the discrimination threshold, it is determined to be an outlier. The discrimination threshold is calculated as follows:

[0222] ,

[0223] in βFor dynamic adjustment coefficients;

[0224] when GRD-LOF If the score does not meet the above two characteristics, the empirical multiplier method is used. If a certain data point GRD-LOF Score exceeded λ Distance threshold times If it is an outlier, then it is considered an anomaly. λ This is the experience multiplier.

[0225] The incremental update module is used to introduce an incremental update mechanism, which only updates the newly added points and their associated information after each round of new data points arrive.

[0226] Incremental update module for sample sliding window updates:

[0227] A fixed-length sliding window mechanism is used to maintain the sample set; at each new time step t+1, a normalized subsequence is extracted from the historical monitoring sequence. Add to the sample set, and simultaneously remove the earliest one. To maintain a constant total sample size of n, a local update is performed on the index structure used for nearest neighbor search, i.e., new sample points are inserted. Delete the excluded points This avoids rebuilding the entire index structure;

[0228] Incremental calculation of LRD and GRD in the incremental update module:

[0229] For the newly added subsequence Re-execute k Nearest neighbor search to determine its nearest neighbor set And calculate its local reachability density based on reachability distance. and global reachability density To maintain the statistical characteristics of the GRD distribution in real time, the following weighted formula is used to calculate the mean of the GRD. and standard deviation Make incremental adjustments:

[0230] ,

[0231] ,

[0232] ,

[0233] in, For newly added sample points GRD value, Points to be removed GRD value, n This represents the total number of samples within the sliding window; This is the variation in standard deviation;

[0234] The incremental update module links the anomaly score with the judgment threshold for continuous updating.

[0235] For new sample points According to its GRD and LRD Instant score calculation GRD-LOF The updated anomaly score sequence is then incorporated into the anomaly score sequence. Based on this updated sequence, the quantiles, mean, and standard deviation are dynamically maintained. The judgment threshold is then dynamically adjusted in real time using preset rules (such as quantile method, standard deviation discrimination method, or empirical multiple method) to achieve highly responsive anomaly detection decisions. The results can be saved as .csv, .json, or database tables for later use.

[0236] Example 3: An electronic device of the present invention includes one or more processors, one or more memories, and one or more programs. The programs are stored in the memory and configured to be executed by the processor. When the programs are loaded into the processor, they implement the steps of an anomaly detection method based on a global reachability density and adaptive threshold mechanism as described above.

[0237] Example 4: A computer-readable storage medium of the present invention stores a computer program, the computer program including program instructions, which, when executed by a processor, cause the processor to perform steps according to an anomaly detection method based on a global reachability density and adaptive threshold mechanism as described above.

Claims

1. A method for detecting anomalies in water quality monitoring data based on global reachability density and adaptive threshold mechanism, characterized in that, Includes the following steps: (1) Collect historical water quality monitoring data of a single indicator from the target monitoring stations, extract fixed-length subsequence samples using the sliding window method to form a sample set, construct a nearest neighbor index structure, and calculate the value of each data point. p of k Nearest neighbor set This allows for the extraction of neighborhood information; (2) Calculate each data point p Local reachability density in its local neighborhood ; (3) Within the sample set, based on k Calculate the reachability distance of the farthest point among the nearest neighbors for each data point. p Global reachability density Step (3) specifically includes the following steps: For each normalized sample subsequence Its global reachability density Defined as the reciprocal of the average reachability distance between itself and all other sample points except itself, the formula is: ， Where n is the total number of samples. To remove Other samples besides; For point With point The reachable distance between them is defined as follows: ， For sample points To its own k The Euclidean distance between the nearest neighbors is calculated using the following formula: ， in s It is the farthest nearest neighbor sample point; For sample points With sample points The Euclidean distance between them is calculated using the following formula: ， in Represents sample points The t One portion, Represents sample points The t One portion, w The length of the sliding window; (4) For all data points in the sample set GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated. The standard deviation is then used as the criterion for further analysis. with the mean The ratio is defined as the coefficient of variation. The dynamic distance threshold is determined by combining the adjustment coefficient α. ; ， (5) Based on local reachability density and global reachability density, construct a nearest neighbor fusion anomaly score. GRD-LOF and with dynamic distance threshold The comparison is performed to achieve automatic identification and judgment of anomalies; step (5) specifically includes the following steps: (5.1) Calculate the nearest neighbor fusion anomaly score for each data point. GRD-LOF, The calculation formula is: ， ， in, For point of k Nearest neighbor set; Point for neighbors Global reachability density; For point Local reachability density; (5.2) When GRD-LOF The scores are not normally distributed; the quantile method is used to analyze all points. GRD-LOF Scores are sorted in ascending order, and a high percentile value is set as an outlier threshold. If a certain data point... GRD-LOF If the score exceeds this percentile value, it is judged as an outlier; when GRD-LOF When the scores exhibit a "centralized distribution with small fluctuations but extreme high values," the standard deviation discriminant method is used to statistically analyze all points. GRD-LOF The mean μ and standard deviation σ of the scores are used to set a threshold for discrimination. If a certain data point... GRD-LOF If the value exceeds the discrimination threshold, it is determined to be an outlier. The discrimination threshold is calculated as follows: ， in β For dynamic adjustment coefficients; when GRD-LOF If the score does not meet the above two conditions, the empirical multiplier method is used. If a certain data point GRD-LOF Score exceeded λ Distance threshold times If it is an outlier, then it is considered an anomaly. λ Experience multiplier; (6) Introduce an incremental update mechanism. After each round of new data points arrive, only the newly added points and their associated information are partially updated.

2. The method for detecting anomalies in water quality monitoring data based on global reachability density and adaptive threshold mechanism according to claim 1, characterized in that, Step (1) specifically includes the following steps: (1.1) Obtain continuous historical monitoring data of a certain water quality indicator from the target monitoring station, construct a one-dimensional time series, denoted as . Where T represents the total number of observation time steps, X For this water quality index in continuous T Historical water quality monitoring sequences at each time step , indicating a point in time t The monitored values; (1.2) For the sequence X Use a length of w The sliding window extracts subsequences to form a sample set: ,in ,in, w For window length, The total number of subsequences, and each sample subsequence Indicates continuity w The indicator monitoring value at each moment; (1.3) For each sample subsequence Perform z-score standardization to obtain normalized sample subsequences. ,in, For sample subsequence The mean, For sample subsequence Standard deviation of the normalized sample This constitutes the normalized sample set. ; (1.4) For the sample set Establish a spatial index structure; (1.5) For the sample set Each normalized subsequence Using the constructed spatial index structure, a k-nearest neighbor query is performed based on Euclidean distance to obtain its k most similar neighbor subsequences: .

3. The method for detecting anomalies in water quality monitoring data based on global reachability density and adaptive threshold mechanism according to claim 2, characterized in that, Step (2) specifically includes the following steps: (2.1) For the normalized sample set any subsequence in , and its first j A neighbor The reachability distance between them is defined as: ， Sample points Compared with its own k The Euclidean distance between the nearest neighbors, representing the lower bound of its local density, is defined as: ， in Represents sample points of k A set of nearest neighbors; For sample points Its neighbors q The Euclidean distance between them is defined as: ,in Represents sample points In the t The value of each component For the nearest neighbor q in the th... t The value of each component; The Euclidean distance between two sample points is defined as: ， in, Represents sample points The t One portion, Represents sample points The t One component; (2.2) Calculate the local reachability density LRD , after defining arbitrary sample points After calculating the reachability distance between each of its neighbors, compute the local reachability density of any subsequence. It is the reciprocal of the average reachable distance within its neighborhood, calculated using the following formula: ， The denominator represents the sample points. To its k The average reachable distance of the nearest neighbor.

4. The anomaly detection method based on global reachability density and adaptive threshold mechanism according to claim 3, characterized in that, Step (6) specifically includes the following steps: (6.1) Sample sliding window update A fixed-length sliding window mechanism is used to maintain the sample set; at each new time step t+1, a normalized subsequence is extracted from the historical monitoring sequence. Add to the sample set, and simultaneously remove the earliest one. To maintain a constant total sample size of n, a local update is performed on the index structure used for nearest neighbor search, i.e., new sample points are inserted. Delete the excluded points This avoids rebuilding the entire index structure; (6.2) Incremental calculation of LRD and GRD For the newly added subsequence Re-execute k Nearest neighbor search to determine its nearest neighbor set And calculate its local reachability density based on reachability distance. and global reachability density To maintain the statistical characteristics of the GRD distribution in real time, the following weighted formula is used to calculate the mean of the GRD. and standard deviation Perform incremental adjustments: ， ， ， in, For newly added sample points GRD value, Points to be removed GRD value, n This represents the total number of samples within the sliding window; This is the variation in standard deviation; (6.3) Linked update of anomaly score and judgment threshold For new sample points According to its GRD and LRD Calculate its anomaly score GRD-LOF The abnormal score sequence is then incorporated into the abnormal score sequence. Based on the updated abnormal score sequence, its quantiles, mean, and standard deviation are dynamically maintained, and the judgment threshold is adjusted in real time.

5. An anomaly detection system based on global reachability density and adaptive threshold mechanism, characterized in that, include: The neighborhood information extraction module is used to collect historical water quality monitoring data for a single indicator from target monitoring stations. It extracts fixed-length sub-sequence samples using a sliding window method to form a sample set, constructs a nearest neighbor index structure, and calculates the nearest neighbor index for each data point. p of k Nearest neighbor set This allows for the extraction of neighborhood information; The local reachability density calculation module is used to calculate the density of each data point. p Local reachability density in its local neighborhood ; Global The reachability density calculation module is used to calculate the reachability density within the sample set based on... k Calculate the reachability distance of the farthest point among the nearest neighbors for each data point. p Global reachability density The global reachability density calculation module calculates each normalized sample subsequence. Its global reachability density Defined as the reciprocal of the average reachability distance between itself and all other sample points except itself, the formula is: ， Where n is the total number of samples. To remove Other samples besides; For point With point The reachable distance between them is defined as follows: ， For sample points To its own k The Euclidean distance between the nearest neighbors is calculated using the following formula: ， in s It is the farthest nearest neighbor sample point; For sample points With sample points The Euclidean distance between them is calculated using the following formula: ， in Represents sample points The t One portion, Represents sample points The t One portion, w The length of the sliding window; The dynamic distance threshold determination module is used to determine the distance threshold for all data points in the sample set. GRD The values ​​are statistically summarized, and their mean and standard deviation are calculated. The standard deviation is then used as the criterion for further analysis. with the mean The ratio is defined as the coefficient of variation. The dynamic distance threshold is determined by combining the adjustment coefficient α. ; ， The anomaly identification and judgment module is used to construct a nearest-neighbor fusion-type anomaly score based on local and global reachability densities. GRD-LOF and with dynamic distance threshold The comparison is performed to achieve automatic identification and judgment of anomalies; the anomaly identification and judgment module calculates the nearest neighbor fusion anomaly score for each data point. GRD-LOF, The calculation formula is: ， ， in, For point of k Nearest neighbor set; Point for neighbors Global reachability density; For point Local reachability density; when GRD-LOF The scores are not normally distributed; the quantile method is used to analyze all points. GRD-LOF Scores are sorted in ascending order, and a high percentile value is set as an outlier threshold. If a certain data point... GRD-LOF If the score exceeds this percentile value, it is judged as an outlier; when GRD-LOF When the scores exhibit a "centralized distribution with small fluctuations but extreme high values," the standard deviation discriminant method is used to statistically analyze all points. GRD-LOF The mean μ and standard deviation σ of the scores are used to set a threshold for discrimination. If a certain data point... GRD-LOF If the value exceeds the discrimination threshold, it is determined to be an outlier. The discrimination threshold is calculated as follows: ， in β For dynamic adjustment coefficients; when GRD-LOF If the score does not meet the above two conditions, the empirical multiplier method is used. If a certain data point GRD-LOF Score exceeded λ Distance threshold times If it is an outlier, then it is considered an anomaly. λ Experience multiplier; The incremental update module is used to introduce an incremental update mechanism, which only updates the newly added points and their associated information after each round of new data points arrive.

6. An anomaly detection system based on global reachability density and adaptive threshold mechanism according to claim 5, characterized in that, The neighborhood information extraction module obtains continuous historical monitoring data of a certain water quality indicator from the target monitoring station, constructs a one-dimensional time series, and denoted as... Where T represents the total number of observation time steps, X For this water quality index in continuous T Historical water quality monitoring sequences at each time step , indicating a point in time t The monitored values; For sequence X Use a length of w The sliding window extracts subsequences to form a sample set: ,in ,in, w For window length, The total number of subsequences, and each sample subsequence Indicates continuity w The indicator monitoring value at each moment; For each sample subsequence Perform z-score standardization to obtain the normalized sample subsequences: ,in, For sample subsequence The mean, For sample subsequence Standard deviation of the normalized sample This constitutes the normalized sample set. ; For sample set Establish a spatial index structure; For sample set Each normalized subsequence Using the constructed spatial index structure, a k-nearest neighbor query is performed based on Euclidean distance to obtain its k most similar neighbor subsequences: .

7. An electronic device, characterized in that: It includes one or more processors, one or more memories, and one or more programs, said programs being stored in the memory and configured to be executed by the processor, said programs being loaded onto the processor to implement the steps of a water quality monitoring data anomaly detection method based on a global reachability density and adaptive threshold mechanism according to any one of claims 1 to 4.

8. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, the computer program including program instructions, which, when executed by a processor, cause the processor to perform the steps of a water quality monitoring data anomaly detection method based on a global reachability density and adaptive threshold mechanism according to any one of claims 1 to 4.

Images