A coal mine industrial control network abnormal flow monitoring method and system

By constructing a two-dimensional feature space and performing clustering in the coal mine industrial control network, and combining density dispersion and Mahalanobis distance variance for weighted sampling, an adaptive isolation tree is constructed and the threshold is dynamically adjusted. This solves the problems of insufficient sparse anomaly sampling and poor adaptability of the isolation tree, and achieves high-precision and stable anomaly flow monitoring.

CN122640206APending Publication Date: 2026-08-25INFORMATION TECH OPERATION & MAINTENANCE BRANCH OF SHAANXI COAL & NORTHERN SHAANXI MINING CO LTD
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202610819317.0
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2026-06-08
Publication Date
2026-08-25

AI Technical Summary

Technical Problem

Existing methods for monitoring abnormal flow in coal mine industrial control networks suffer from insufficient sparse anomaly sampling, poor adaptability of isolation trees, and static thresholds, resulting in inadequate detection accuracy and stability, making them unsuitable for monitoring abnormal flow in complex industrial control environments.

Method used

By constructing a two-dimensional feature space and performing clustering, weighted sampling is performed by combining density dispersion and Mahalanobis distance variance to compensate for sparse clusters, an adaptive isolation tree is constructed, and the judgment threshold is dynamically adjusted to improve detection accuracy and stability.

Benefits of technology

Under the premise of controlled total sampling scale, it effectively protects sparse abnormal traffic, improves the identification accuracy of isolation trees for complex traffic boundaries, reduces false alarm rate and false negative rate, and improves the accuracy and stability of anomaly monitoring.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122640206A_ABST
    Figure CN122640206A_ABST
Patent Text Reader

Abstract

The present application belongs to the technical field of flow monitoring, and particularly relates to a coal mine industrial control network abnormal flow monitoring method and system, which comprises the following steps: constructing a two-dimensional feature space from the load byte entropy value and the time stamp interval, clustering and dividing mutually intersecting flow clusters; determining an initial subsampling number, compensating sparse clusters to a basic share; apportioning and deducting excess samples in inverse proportion to Mahalanobis distance variance, maintaining the sampling scale unchanged; constructing an isolation tree for each subsampling set, setting a non-uniform division probability according to the feature variance, determining the maximum depth of the tree through the trace of the covariance matrix and an S-shaped function; calculating a single-tree score according to the average path length of the sample, and obtaining a comprehensive abnormal score through weighted geometric mean according to the cluster weight; adjusting the quantile threshold value by real-time calculation of the cluster weight information entropy, triggering an abnormal alarm if the score exceeds the limit, and realizing the identification of coal mine industrial control network flow anomalies. The present application can compensate sparse clusters, adaptively isolate trees, dynamically adjust threshold values and improve detection accuracy.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of flow monitoring technology. More specifically, this invention relates to a method and system for monitoring abnormal flow in coal mine industrial control networks. Background Technology

[0002] With the continuous advancement of intelligent coal mine construction, coal mine industrial control networks have become a crucial hub connecting underground production equipment with ground dispatch and control centers, and the integration of IT and OT systems is deepening. Coal mine industrial control networks employ diverse protocol types, and their communication behavior typically exhibits strong periodicity and determinism. Once affected by network attacks, equipment failures, or the injection of abnormal commands, it can lead to production anomalies, interruptions, and even safety risks. Therefore, accurate and timely anomaly monitoring of coal mine industrial control network traffic is a vital link in ensuring safe coal mine production and the stable operation of industrial control systems.

[0003] In real-world coal mine industrial control network environments, normal business traffic typically constitutes the vast majority, while abnormal attack traffic or fault traffic is less numerous and more covertly distributed, exhibiting a clear sample imbalance characteristic. Anomaly detection methods based on a combination of clustering and isolated forest algorithms usually first cluster the traffic feature space and then construct isolation trees through a subsampling mechanism to calculate anomaly scores. However, existing methods often use uniform or simple proportional allocation when allocating subsampling sizes, failing to fully consider the differences in sample density dispersion and multidimensional spatial distribution among different traffic clusters. This can easily lead to sparse clusters containing rare abnormal traffic being weakened during the sampling process, making it difficult to balance the compensation of sparse clusters and the deduction of redundant clusters while keeping the total sampling size under control.

[0004] Furthermore, in the process of constructing the isolation tree, existing methods typically employ completely random partitioning rules and fixed maximum depths, failing to fully utilize spatial structural information such as flow characteristic variance and covariance matrix. This results in insufficient partitioning efficiency of the isolation tree and its inability to adapt to different flow cluster distribution characteristics. In the anomaly detection stage, traditional methods often rely on simple arithmetic average scoring and static fixed thresholds, making it difficult to adapt to the distribution fluctuations of coal mine industrial control network flow caused by changes in production conditions and time windows. They also lack the ability to dynamically adjust for changes in flow cluster weight distribution, easily leading to false alarms or missed alarms in complex industrial control environments, affecting the accuracy and stability of abnormal flow monitoring. Summary of the Invention

[0005] To address the technical problems of insufficient sparse anomaly sampling, poor adaptability of isolation trees, and static thresholds in existing methods, this invention proposes a method and system for monitoring abnormal traffic in coal mine industrial control networks. This system can compensate for sparse clusters, adapt to isolation trees, and dynamically adjust thresholds to improve detection accuracy.

[0006] In a first aspect, the present invention provides a method for monitoring abnormal traffic in a coal mine industrial control network, comprising: S1, collecting coal mine industrial control protocol traffic data, calculating the byte entropy value of the traffic packet payload and the adjacent timestamp interval to construct a two-dimensional feature space and dividing it into several non-overlapping traffic clusters by clustering; S2, using the product of the inverse of the sample density dispersion within each traffic cluster and the variance of the Mahalanobis distance as the initial weight to determine the initial number of subsamples according to the proportion, and compensating sparse clusters with density dispersion higher than the judgment threshold to a preset basic share; S3, distributing and deducting the extra over-allocation quota occupied by sparse clusters according to the inverse of the square root of the variance of the Mahalanobis distance of the remaining traffic clusters, and randomly and uniformly removing redundant samples from the deducted traffic clusters to obtain each subsample set; 4. Construct isolation trees for each subset of samples. Randomly select the splitting dimension with a non-uniform probability set according to the feature variance ratio, and randomly determine the splitting point on the selected splitting dimension. Calculate the trace of the two-dimensional covariance matrix, standardize it, map it with the S-type activation function, multiply it by the base depth value, and round it up to obtain the maximum depth limit of the isolation tree. S5. Substitute the average path length of the traffic sample in each isolation tree into the evaluation function to calculate the single tree score, and perform a weighted geometric average according to the traffic cluster weight to obtain the comprehensive anomaly score. S6. Based on the information entropy of the weight of each traffic cluster in the current time window, adjust the base quantile constant, and use the historical comprehensive anomaly score of the corresponding new quantile as the judgment threshold. When the comprehensive anomaly score exceeds the judgment threshold, output an anomaly alarm.

[0007] By employing the above technical solution, a two-dimensional feature space is constructed and clustered based on the payload byte entropy value and timestamp interval. Weighted sampling is performed by combining density dispersion and Mahalanobis distance variance. After compensating the basic share for sparse clusters, redundant samples are proportionally deducted. Then, an isolation tree with adaptive feature variance segmentation and dynamic depth adjustment is constructed. Finally, the judgment threshold is dynamically adjusted based on the cluster weight information entropy. This method can effectively protect sparse abnormal traffic under controlled total sampling scale, improve the identification accuracy of the isolation tree for complex traffic boundaries, and adaptively adjust the alarm threshold according to operating condition fluctuations, significantly reducing false alarm and false negative rates, and improving the accuracy and stability of anomaly monitoring.

[0008] Preferably, the step of constructing a two-dimensional feature space by calculating the entropy value of the traffic packet payload byte and the adjacent timestamp interval, and then dividing it into several non-overlapping traffic clusters by clustering, includes: statistically analyzing the frequency of each byte value in the payload of a single data packet, and calculating the payload byte entropy value according to the Shannon information entropy formula; extracting the arrival time of the current traffic packet and the previous adjacent traffic packet, and calculating the absolute difference as the timestamp interval; standardizing the payload byte entropy value and the timestamp interval using Z-standard fraction processing, and combining them to form two-dimensional feature coordinate points; and using a density-based noise application space clustering algorithm to iteratively divide the coordinate points in the two-dimensional feature space to form several non-overlapping traffic clusters.

[0009] By employing the above technical solution, Shannon entropy values ​​are calculated based on the byte frequency of message payloads, and the timestamp intervals of adjacent packets are extracted. After Z-normalization, two-dimensional feature coordinates are formed, and the DBSCAN density clustering algorithm is used to divide non-overlapping traffic clusters. This method can effectively characterize the content randomness and temporal regularity of coal mine industrial control protocol traffic, and automatically identify traffic groups of different densities without pre-setting the number of clusters, providing a reasonable cluster structure basis for subsequent differentiated sampling and anomaly detection.

[0010] Preferably, the step of determining the initial subsample number by using the product of the inverse of the sample density dispersion within each flow cluster and the variance of the Mahalanobis distance as the initial weight and proportionally includes: calculating the two-dimensional covariance matrix and pseudo-inverse matrix of each flow cluster; obtaining the Mahalanobis distance from each sample point within the flow cluster to the center of the flow cluster; and calculating the variance of the Mahalanobis distance within each flow cluster; calculating the average Euclidean distance from each sample point within the flow cluster to its K nearest neighbors; taking the mean of the Euclidean distances as the density dispersion of the flow cluster; multiplying the inverse of the density dispersion of the corresponding flow cluster by the variance of the Mahalanobis distance to obtain the absolute weight; calculating the proportion of the absolute weight of each flow cluster to the total absolute weight of all flow clusters as the initial weight; and multiplying the total sampling scale constant by the initial weight of each flow cluster to obtain the theoretical subsample number of each flow cluster.

[0011] By employing the above technical solution, the Mahalanobis distance variance and sample density dispersion of each traffic cluster are calculated. These two values ​​are multiplied to obtain the initial weight, and the total sampling scale is then allocated according to the weight ratio to obtain the theoretical number of subsamples. This method fully considers the spatial distribution breadth and local density fluctuations of samples within a cluster, giving higher sampling weights to compact and uniformly dense normal business clusters, while giving lower weights to loosely structured and density-dispersed potential anomalous clusters. This initially suppresses redundant normal samples and retains valuable anomalous samples during the sampling process.

[0012] Preferably, the step of determining the initial subsample number by using the product of the inverse of the sample density dispersion within each flow cluster and the variance of the Mahalanobis distance as the initial weight according to the proportion further includes: rounding down the theoretical subsample number to obtain a temporary subsample number, and calculating the difference between the total sampling size constant and the sum of all temporary subsample numbers as the rounding remaining quota; according to the order of the decimal part of the theoretical subsample number of each flow cluster from largest to smallest, sequentially adding the rounding remaining quota to the corresponding flow cluster to obtain the initial subsample number of each flow cluster, so that the sum of all initial subsample numbers is equal to the total sampling size constant.

[0013] By employing the above technical solution, the theoretical subsample number is rounded down, and the remaining quota is supplemented sequentially from largest to smallest fractional part, ensuring that the sum of the initial subsample numbers of each flow cluster is precisely equal to the total sampling scale constant. This method solves the quota loss problem caused by integer rounding, ensures that the total number of samples is strictly conserved, and maintains the relative fairness of the sampling ratio among clusters, providing an accurate numerical basis for subsequent sparse cluster compensation and redundancy deduction.

[0014] Preferably, the step of compensating sparse clusters with density dispersion higher than the judgment threshold to a preset basic share includes: extracting the system's preset density dispersion judgment threshold and preset basic share constant; traversing all traffic clusters and selecting traffic clusters with density dispersion greater than the judgment threshold to form a sparse cluster set; calculating the difference between the current initial subsample number of the sparse cluster and the preset basic share constant; if the initial subsample number is lower than the preset basic share constant, then setting the subsample number of the sparse cluster to the preset basic share, and recording the total amount of additional over-allocation generated by the compensation.

[0015] By employing the above technical solution, sparse clusters are selected, and their initial subsample count is forcibly increased to a preset base share, while additional over-allocation quotas are recorded. This method can proactively protect sparse clusters with few samples and scattered distributions that may contain abnormal traffic, preventing them from being excessively reduced in proportional allocation, and ensuring that key abnormal samples receive sufficient sampling opportunities. This is a crucial step in improving the ability to detect hidden anomalies.

[0016] Preferably, the step of allocating the extra over-allocation quota occupied by sparse clusters according to the reciprocal of the square root of the variance of the Mahalanobis distance of each remaining flow cluster, and randomly and uniformly removing redundant samples from the deducted flow clusters to obtain each sub-sample set, includes: removing the sparse cluster set from all flow clusters to obtain the remaining flow cluster set; calculating the square root of the variance of the Mahalanobis distance of each flow cluster in the remaining flow cluster set, adding the square root to a very small positive number and taking the reciprocal, and using the proportion of each reciprocal to the sum of reciprocals as the allocation coefficient; multiplying the total extra over-allocation quota by the allocation coefficient of each remaining flow cluster, rounding down to obtain the initial number of samples to be deducted for each remaining flow cluster, and allocating the difference between the total extra over-allocation quota and the sum of all initial sample numbers to the corresponding remaining flow clusters in descending order of the deduction quota according to the decimal part of the allocation coefficient to supplement the deduction quota, thus obtaining the deduction quantity for each remaining flow cluster.

[0017] Preferably, obtaining each sub-sample set further includes: in the corresponding remaining flow cluster, using a random uniform sampling method to remove redundant samples of the same number as the number of samples deducted to obtain each sub-sample set, and after removal, verifying that the sum of the sample numbers of all sub-sample sets is equal to a preset total sampling scale constant.

[0018] Preferably, constructing an isolation tree for each sub-sample set includes: calculating the feature variances of the payload byte entropy dimension and the timestamp interval dimension in the sub-sample set respectively; using the ratio of the one-dimensional variance to the sum of the two-dimensional variances as the probability of selecting the dimension for splitting; randomly selecting the splitting dimension according to the probability; and uniformly and randomly selecting splitting points between the maximum and minimum values ​​of the dimension; extracting the sum of the main diagonal elements of the two-dimensional covariance matrix of the sub-sample set as the trace of the matrix; standardizing the trace of the matrix and inputting it into a sigmoid activation function to map the value to the interval 0 to 1; multiplying it by the base depth value; and performing an up-rounding operation on the product result. The resulting integer is the maximum depth limit allowed for the isolation tree to grow during construction.

[0019] Preferably, the step of adjusting the basic quantile constant based on the weighted information entropy of each traffic cluster during the current time window, and using the historical comprehensive anomaly score corresponding to the new quantile as the judgment threshold, includes: Extract the normalized initial weight probability distribution of all traffic clusters within the current time window, and calculate the traffic cluster weight information entropy for the current time window period according to the information entropy formula; calculate the difference rate between the current traffic cluster weight information entropy and the historical traffic cluster weight information entropy of the previous time window; if the difference rate is greater than 0, the basic quantile constant is amplified proportionally to the difference rate and truncated to an upper limit of 1 as the new quantile; if the difference rate is not greater than 0, the basic quantile constant is reduced proportionally to the absolute value of the difference rate as the new quantile; use the score value corresponding to the new quantile in the historical comprehensive anomaly score distribution as the judgment threshold.

[0020] Secondly, the present invention provides a coal mine industrial control network abnormal flow monitoring system, including a processor and a memory, wherein the memory stores computer program instructions, and when the computer program instructions are executed by the processor, the above-mentioned coal mine industrial control network abnormal flow monitoring method is implemented.

[0021] By adopting the above technical solution, a computer program for monitoring abnormal traffic in a coal mine industrial control network is generated and stored in a memory so that it can be loaded and executed by a processor. This allows for the creation of terminal devices based on the memory and processor, making them convenient to use.

[0022] The present invention has the following beneficial effects: This invention constructs a feature space by calculating the payload byte entropy value and timestamp interval and performs clustering to extract coal mine industrial control protocol features to distinguish various traffic behaviors; in the sampling stage, it combines density dispersion and Mahalanobis distance variance to allocate weights, gives basic share compensation to sparse clusters to prevent rare abnormal traffic from being missed, and deducts redundant samples proportionally while maintaining the total sampling scale, effectively controlling computational costs and optimizing the quality of model input features.

[0023] Furthermore, in the construction of the isolation tree, this invention relies on the non-uniform selection of the segmentation dimension based on the feature variance ratio, and uses the trace of the covariance matrix to define the upper limit of the maximum depth, thereby enhancing the recognition accuracy of complex traffic boundaries. At the same time, it dynamically adjusts the judgment threshold according to the fluctuation law of the information entropy of traffic cluster weights within the time window, so that the abnormal alarm mechanism fits the actual network operation status, thereby significantly reducing the false alarm rate and false negative rate in the complex industrial control environment of coal mines, and improving the accuracy and reliability of network security monitoring. Attached Figure Description

[0024] Figure 1 This is a flowchart of a method for monitoring abnormal traffic in a coal mine industrial control network. Figure 2 This is a schematic diagram of two-dimensional feature clustering; Figure 3 This is a schematic diagram of the threshold adjustment mechanism; Figure 4 This is a diagram showing the results of a comparative test. Detailed Implementation

[0025] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0026] In the first embodiment, the present invention proposes a method for monitoring abnormal traffic in a coal mine industrial control network, such as... Figure 1 As shown, it includes the following steps: S1 collects coal mine industrial control protocol traffic data, calculates the byte entropy value of the traffic packet payload and the adjacent timestamp interval to construct a two-dimensional feature space, and then divides it into several non-overlapping traffic clusters through clustering.

[0027] Network probes are deployed via switch port mirroring. A packet sniffing function is used to monitor, identify, and collect Modbus TCP and S7Comm protocol packets transmitted in the coal mine control network in real time. The application layer payload of the packets is extracted, and Shannon entropy is calculated based on byte frequency distribution as the first feature. The arrival timestamp of the current packet is subtracted from the arrival timestamp of the previous packet to obtain the time difference as the second feature. Z-score standardization or minus-max normalization is performed on both the first and second features to eliminate the dimensional and scale differences between information entropy and time intervals. The standardized first and second features are combined to obtain a two-dimensional feature matrix, which is then stored. The KMeans mean clustering algorithm or the DBSCAN density clustering algorithm is called, and iterative spatial clustering is performed on the two-dimensional feature matrix using model fitting and prediction methods. This divides the entire coal mine control network traffic data into multiple traffic clusters that have no overlap or intersection.

[0028] In an optional embodiment, the process of calculating the entropy value of the traffic packet payload byte and the adjacent timestamp interval to construct a two-dimensional feature space and then clustering it into several non-overlapping traffic clusters includes: statistically analyzing the frequency of each byte value within the payload of a single data packet and calculating the payload byte entropy value according to the Shannon information entropy formula; extracting the arrival time of the current traffic packet and the previous adjacent traffic packet and calculating the absolute difference as the timestamp interval; standardizing the payload byte entropy value and the timestamp interval using Z-standard scores and combining them to form two-dimensional feature coordinate points; and using a density-based noise application space clustering algorithm to iteratively divide the coordinate points in the two-dimensional feature space to form several non-overlapping traffic clusters.

[0029] When collecting coal mine industrial control protocol traffic data, for a single data packet with a length of, for example, 64 to 1500 bytes, traverse every byte of the payload portion and count the frequency of each byte value. For i ∈ [0, 255], use Shannon's information entropy formula. The payload byte entropy value is calculated, ranging from 0 to 8 bits. The arrival time of the current traffic packet is extracted. Arrival time of the previous adjacent data packet connected to the same quintuple Calculate the absolute value of the difference between the two. As timestamp intervals, typical values ​​are usually distributed within a continuous range of 1ms to 500ms. After obtaining this two-dimensional feature, a statistical sequence containing N, for example, 10,000 samples is collected within a preset time window, such as 5 minutes.

[0030] Extract the sample mean of the payload byte entropy set and the timestamp interval set within the time window, respectively. and standard deviation Execute the Z-standard score formula Perform Gaussian normalization to normalize the entropy features. With time interval characteristics Combined into two-dimensional feature coordinate points A two-dimensional feature space with Cartesian coordinates is constructed. The DBSCAN clustering algorithm is then used for unsupervised scanning learning of this feature space. The neighborhood search radius is set between 0.10 and 0.25, and the minimum number of samples required to trigger the formation of core objects within the neighborhood is set between 5 and 20. By continuously calculating the Euclidean distance between coordinate points in the space, connected sample points with reachable density are merged into clusters. Density iterative optimization generates, for example, 3 to 8 independent flow clusters with no overlapping boundaries, as well as scattered noise outliers. This achieves the decoupling and classification of normal polling traffic, abnormal sudden large commands, or congestion delay responses onto the manifold space. The two-dimensional feature clustering is as follows: Figure 2 As shown.

[0031] S2 uses the product of the inverse of the sample density dispersion within each flow cluster and the variance of the Mahalanobis distance as the initial weight to determine the initial number of subsamples according to the proportion, and compensates the sparse clusters with density dispersion higher than the judgment threshold to the preset basic share.

[0032] For each defined traffic cluster, the KDTree spatial indexing algorithm is used to calculate the Euclidean distance from each traffic sample within the cluster to its K nearest neighbors. The standard deviation or mean of these distances is then used as the sample density dispersion of the traffic cluster. The covariance matrix of the two-dimensional features within the cluster is calculated and combined with the pseudo-inverse function of the linear algebra submodule to obtain the pseudo-inverse covariance matrix. The Mahalanobis distance from the samples within the cluster to the cluster center is calculated, and the variance is calculated using the variance calculation function. The inverse of the sample density dispersion is multiplied by the Mahalanobis distance variance to obtain the absolute weight of each traffic cluster. The absolute weight of a single cluster is divided by the sum of the absolute weights of all clusters to obtain the normalized initial weight. This initial weight is multiplied by the system's preset total sampling size constant to obtain the theoretical subsample number for each traffic cluster. The theoretical subsample number is rounded down to obtain the temporary subsample number. The remaining quota is then allocated in descending order of the fractional part of the theoretical subsample number, ensuring that the sum of the initial subsample numbers for all traffic clusters equals the total sampling size constant. Set a threshold for determining the system density dispersion and iterate through each flow cluster. When the sample density dispersion value of a flow cluster is greater than the threshold, the flow cluster is marked as a sparse cluster, and the initial subsample value of the sparse cluster is modified and increased to the base share value preset by the system to complete the compensation operation.

[0033] In an optional embodiment, the initial number of subsamples is determined proportionally using the product of the inverse of the sample density dispersion within each flow cluster and the variance of the Mahalanobis distance as the initial weight. This includes: calculating the two-dimensional covariance matrix and pseudo-inverse matrix of each flow cluster; obtaining the Mahalanobis distance from each sample point within the flow cluster to the center of the flow cluster; and calculating the variance of the Mahalanobis distance within each flow cluster; calculating the average Euclidean distance from each sample point within the flow cluster to its K nearest neighbors; taking the mean of the Euclidean distances as the density dispersion of the flow cluster; and multiplying the inverse of the density dispersion of the corresponding flow cluster by the variance of the Mahalanobis distance to obtain the absolute weight. The initial weight is calculated as the proportion of the absolute weight of each traffic cluster to the total absolute weight of all traffic clusters. The total sampling scale constant is multiplied by the initial weight of each traffic cluster to obtain the theoretical subsample number of each traffic cluster. The theoretical subsample number is rounded down to obtain the temporary subsample number. The difference between the total sampling scale constant and the sum of all temporary subsample numbers is calculated as the rounding remaining quota. The rounding remaining quota is sequentially added to the corresponding traffic clusters in descending order of the decimal part of the theoretical subsample number of each traffic cluster to obtain the initial subsample number of each traffic cluster, so that the sum of all initial subsample numbers is equal to the total sampling scale constant.

[0034] For each independent flow cluster identified by clustering, extract the feature coordinates of all samples within the cluster to construct a 2×2 two-dimensional covariance matrix. And the pseudo-inverse matrix is ​​obtained using the singular value decomposition method. To avoid computational divergence when highly repetitive coal mine control commands lead to non-full-rank and non-invertible matrices, the arithmetic mean centroid coordinates of the flow cluster in the feature space are calculated as the center point. The Mahalanobis distance from each sample point within the cluster to this center point is calculated using a pseudo-inverse matrix, and the variance of the Mahalanobis distance within each flow cluster is statistically analyzed. This method quantitatively characterizes the distribution breadth of the cluster in the feature space. The K-nearest neighbor Euclidean distance method is used, with K configured between 5 and 15, preferably 10. The mean spatial Euclidean distance from each sample point within the cluster to its K nearest neighbors is calculated, and the arithmetic mean of these mean distances for all samples within the cluster is defined as the density dispersion of that specific flow cluster. The density dispersion of normal business clusters is low, while the dispersion of abnormal clusters with covert detection characteristics is high.

[0035] After obtaining the two evaluation variables mentioned above, a risk-aware weight fusion strategy is executed, which combines the reciprocal of the dispersion of the corresponding traffic cluster density with the variance of the Mahalanobis distance. The total base is obtained by summing the absolute weights of all traffic clusters within the current time window, and then the absolute weights of each independent cluster are calculated. Divide by the total base to obtain the normalized initial weights. Read the preset total sampling size constant, and combine the constant with the initial weights of each independent cluster. Multiply to obtain the theoretical subsample number, round down the theoretical subsample number, and then allocate the remaining quota in descending order of the fractional part of the theoretical subsample number to obtain the initial subsample number for each flow cluster. and make all The sum of these equals the total sampling size constant.

[0036] In an optional embodiment, sparse clusters with density dispersion higher than a judgment threshold are compensated to a preset base share, including: extracting the system's preset density dispersion judgment threshold and a preset base share constant; traversing all traffic clusters and selecting traffic clusters with density dispersion greater than the judgment threshold to form a sparse cluster set; calculating the difference between the current initial subsample number of the sparse cluster and the preset base share constant; if the initial subsample number is lower than the preset base share constant, then setting the subsample number of the sparse cluster to the preset base share, and recording the total amount of additional over-allocation generated by the compensation.

[0037] The system reads a preset density dispersion threshold constant, which, based on historical flow distribution characteristics, typically ranges from 0.15 to 0.35, preferably 0.25, along with a preset basic share constant. After obtaining the threshold, it iterates through all the divided flow clusters, compares the actual density dispersion index of the current flow cluster, and selects flow clusters with dispersion greater than the threshold to construct a sparse cluster set.

[0038] Upon entering the compensation phase, for each sparse cluster in the sparse cluster set, the initial subsample number allocated to the sparse cluster according to the aforementioned ratio is obtained, and the difference between this initial value and the preset basic share constant is calculated. If the initial subsample number is lower than the preset basic share constant, the subsample number of the sparse cluster is set and increased to the preset basic share. The additional allocation amount generated by each compensation operation is accumulated and recorded. After completing the traversal, the total amount of additional over-allocation is calculated. This integer result will serve as the baseline parameter for subsequent global traffic rebalancing and amortization operations.

[0039] S3, the extra over-domination quota occupied by the sparse cluster is distributed and deducted according to the inverse of the square root of the variance of the Mahalanobis distance of the remaining flow clusters, and redundant samples are randomly and uniformly removed from the deducted flow clusters to obtain each subset.

[0040] The additional over-dispatch amount is obtained by subtracting the sum of the initial subsamples originally allocated to the sparse clusters from the sum of the base share values ​​after compensation for all sparse clusters. For the remaining ordinary flow clusters not marked as sparse clusters, the square root of the variance of the Mahalanobis distance obtained in the previous calculation stage is extracted, and the reciprocal is taken. The reciprocal of the individual remaining ordinary flow cluster is divided by the sum of the reciprocals of all remaining ordinary flow clusters to obtain the corresponding deduction ratio. The aforementioned additional over-dispatch amount is multiplied by the deduction ratio of each remaining ordinary flow cluster to obtain the theoretical deduction amount for each remaining ordinary flow cluster. The theoretical deduction amount is rounded down, and the remaining deduction quota is allocated in descending order according to the decimal part of the theoretical deduction amount, so that the sum of all deduction amounts equals the additional over-dispatch amount. The sampling selection function without replacement is called to randomly select redundant samples equal to the deduction amount from the initial subsample candidate sample set corresponding to each remaining ordinary flow cluster, and the redundant samples are removed from the set using the deletion function. After this adjustment of scale, the sum of the final subsamples of all flow clusters can be made equal to the initial preset total sampling scale constant, thus obtaining the subsample sets corresponding to each coal mine industrial control flow cluster.

[0041] In an optional embodiment, the excess over-dispatch quota occupied by sparse clusters is allocated and deducted in proportion to the reciprocal of the square root of the variance of the Mahalanobis distance of the remaining flow clusters. Redundant samples are then randomly and uniformly removed from the deducted flow clusters to obtain each subset, including: removing the sparse cluster set from all flow clusters to obtain the remaining flow cluster set; calculating the square root of the variance of the Mahalanobis distance of each flow cluster in the remaining flow cluster set, adding the square root to a minimum positive number, taking the reciprocal, and using the proportion of each reciprocal to the sum of reciprocals as the allocation coefficient; multiplying the total excess over-dispatch quota by... The allocation coefficient of each remaining flow cluster is rounded down to obtain the initial sample quantity to be deducted for each remaining flow cluster. The difference between the total amount of extra over-allocation quota and the sum of all initial sample quantities is allocated to the corresponding remaining flow clusters in descending order of the deduction quota according to the fractional part of the allocation coefficient to supplement the deduction quota, thus obtaining the deduction quantity for each remaining flow cluster. In the corresponding remaining flow cluster, redundant samples with the same deduction quantity are removed by random uniform sampling to obtain each sub-sample set. After removal, it is verified that the sum of the sample quantities of all sub-sample sets is equal to the preset total sampling scale constant.

[0042] After the compensation operation on sparse clusters causes the global sampling sum to exceed the preset total size, in order to maintain the conservation of the total sampling size, the excess portion needs to be amortized and deducted to the non-sparse remaining traffic clusters. By removing the above sparse cluster set from all traffic clusters, a set of remaining traffic clusters consisting of normal service traffic is obtained. The pre-calculated Mahalanobis distance variance of each cluster in this set of remaining traffic clusters is obtained, and the square root of the variance is calculated to shrink the numerical scale and avoid individual values ​​being too large. A very small positive number is added, such as... To prevent overflow errors caused by division by zero, the reciprocals are calculated. The proportion of each reciprocal to the sum of the reciprocals of the remaining flow clusters is then used as the allocation factor for that flow cluster. This mechanism allows flow clusters with smaller Mahalanobis variance and more stable internal structures to bear a higher deduction ratio, thereby minimizing the sampling attenuation of flow clusters with large distribution variations.

[0043] After calculating the allocation factor, the allocation factor is combined with the total amount of additional overspending generated above. Multiply to obtain the predicted deduction value in real numbers. Round down to obtain the initial deduction amount for each cluster; sum all the initial deduction amounts and add them to the total outstanding amount. The difference is calculated to determine the remaining gaps to be allocated. The fractional parts of the first-round product results are sorted in descending order, and one deduction quota is allocated sequentially to the corresponding ranked remaining flow cluster until the gaps are allocated, generating integer deduction instructions for each remaining flow cluster. Within each remaining flow cluster, redundant samples equal to the deduction quantity are uniformly removed using pseudo-random sampling. The samples remaining after this removal constitute a subset equal to the preset total sampling size constant, used for subsequent training of the isolation tree model.

[0044] S4. Construct an isolation tree for each subset of samples. Randomly select the splitting dimension with a non-uniform probability set according to the feature variance ratio, and randomly determine the splitting point on the selected splitting dimension. Calculate the trace of the two-dimensional covariance matrix, standardize it, map it with the S-type activation function, multiply it by the base depth value, and round it up to get the maximum depth limit of the isolation tree.

[0045] Each generated sub-sample set is traversed and an isolation tree is constructed recursively from top to bottom. When selecting the splitting feature dimension for each generated internal node, the feature variance of the current sub-sample set in both the payload byte entropy dimension and the adjacent timestamp interval dimension is calculated. The variance of a single dimension is divided by the sum of the variances of the two dimensions to obtain the feature variance ratio. This feature variance ratio is used as the parameter vector of a multinomial probability distribution, and unequal probability sampling is performed to determine the feature dimension used for this split. A uniformly distributed random number generation function is used to select the specific splitting point between the maximum and minimum values ​​contained in this dimension. Before each split, the covariance matrix of the two-dimensional feature matrix of the sample set corresponding to the current node is calculated, and the trace summation function is called to obtain the sum of the main diagonal elements of the matrix, i.e., the trace of the matrix. Obtain the mean and standard deviation of the traces of multiple subset matrices within the training period or the current time window. Subtract the mean from the trace of the current matrix and divide by the sum of the standard deviation and the smallest positive number to complete the Z-score standardization. Input the standardized result into a Sigmoid activation function with the natural constant as the base and nonlinearly map it to a floating-point value between 0 and 1. Multiply the floating-point value by the base depth parameter value configured in the system initialization and use the floor function to convert the multiplication result into an integer. Set this integer result as the maximum depth limit of the isolation tree currently being built. When the recursive growth depth of the tree reaches this limit, stop splitting downwards and convert the current node into a leaf node.

[0046] In an optional embodiment, an isolation tree is constructed for each sub-sample set. A splitting dimension is randomly selected with a non-uniform probability based on the feature variance ratio, and splitting points are randomly determined on the selected splitting dimension. The trace of the two-dimensional covariance matrix is ​​calculated, standardized, mapped by a sigmoid activation function, multiplied by the base depth value, and rounded up to obtain the maximum depth limit of the isolation tree. This includes: calculating the feature variances of the payload byte entropy dimension and the timestamp interval dimension in the sub-sample set respectively; using the ratio of the one-dimensional variance to the sum of the two-dimensional variances as the probability of splitting the selected dimension; randomly selecting a splitting dimension based on the probability; uniformly and randomly selecting splitting points between the maximum and minimum values ​​of the dimension; extracting the sum of the main diagonal elements of the two-dimensional covariance matrix of the sub-sample set as the trace of the matrix; standardizing the trace of the matrix and inputting it into a sigmoid activation function to map the value to the interval 0 to 1; multiplying it by the base depth value; and performing a rounding up operation on the product result. The resulting integer is the maximum depth limit allowed for the isolation tree during construction.

[0047] When constructing isolation tree nodes for each subset of samples, to optimize the random feature selection strategy, features of data samples within the node are extracted, and the feature variances of the payload byte entropy value dimension and the timestamp interval dimension are calculated separately. The sum of these two one-dimensional variances is used as the probability of selecting the corresponding dimension for segmentation, with the ratio of each one-dimensional variance to the sum of the two-dimensional variances as the probability. A probability wheel is constructed in this way, so that dimensions with larger variances have a higher selection probability. After randomly selecting the segmentation dimension based on this probability, the maximum and minimum values ​​of the samples within the current node in that dimension are obtained. Within the numerical interval formed by these two values, a value is randomly selected using a uniform distribution as the segmentation point for that node.

[0048] When setting the maximum depth limit for the isolation tree, the two-dimensional covariance matrix of the subsampled sets is extracted. The trace of the matrix is ​​obtained by summing the elements on the main diagonal, and this trace is used as a measure of the discreteness and fluctuation intensity of the feature space. Z-normalization based on the statistical parameters of the trace values ​​of multiple subsampled sets is performed on the trace of the matrix, and a sigmoid activation function is input for nonlinear mapping. This mapped value is multiplied by a preset base depth constant, and the product is rounded up. This generates a maximum growth depth limit for the isolation tree that is strongly correlated with the discreteness of the current sample data, avoiding ineffective tree depth extension and wasted computational power.

[0049] S5, substitute the average path length of the traffic sample in each isolation tree into the evaluation function to calculate the single tree score, and perform a weighted geometric average according to the traffic cluster weight to obtain the comprehensive anomaly score.

[0050] For each communication traffic sample to be tested that is connected to the coal mine industrial control network in real time, the sample is traversed from top to bottom through all the completed isolation trees starting from the root node according to the feature dimensions and split points set by each node. The number of decision edges traversed by the sample to reach the leaf node is recorded in a loop as the basic path length. Considering the limitation of the maximum depth of the isolation tree in the previous steps, the harmonic series formula of the average path length when the standard binary search tree fails is called. The compensation value of the path where the leaf node cannot be fully expanded is accumulated according to the number of remaining samples in the leaf node, so as to obtain the actual average path length of the sample to be tested on each single isolation tree. The actual average path length is divided by the theoretical average path length calculated using the normalized formula based on the total number of samples in the corresponding sub-sample set. The negative form of the division result is then used as the exponent variable of a base-2 exponential function to calculate the single-tree score of the corresponding isolation tree. The initial weights of each traffic cluster, which have been calculated and normalized in the second stage, are extracted as exponential weights. The single-tree scores calculated for each isolation tree are then subjected to idempotent operations of the corresponding exponential weight powers, and all the exponentiation results are multiplied together. Alternatively, a very small positive number is added to the single-tree scores of each isolation tree, the natural logarithm is taken, and the scores are weighted and summed according to the initial weights of the corresponding traffic clusters. Then, exponential mapping is performed to obtain a weighted geometric mean comprehensive anomaly score. This comprehensive anomaly score is used to represent the overall deviation of the measured communication traffic from the historical normal traffic cluster structure.

[0051] S6 adjusts the basic quantile constant based on the entropy of the weight information of each traffic cluster in the current time window, and uses the historical comprehensive anomaly score corresponding to the new quantile as the judgment threshold. When the comprehensive anomaly score exceeds the judgment threshold, an anomaly alarm is output.

[0052] In the operation of the coal mine control anomaly monitoring system, a fixed-length sliding time window is set. The discrete probability distribution sequence of the normalized initial weights of each flow cluster generated in all data update cycles within this time window is extracted. The Shannon information entropy formula is used to sum the products of all probability distribution values ​​and their corresponding base-2 logarithmic values, and the negative sign is taken to obtain the weight information entropy within the current sliding time window. The weight information entropy values ​​stored in the current sliding time window are compared with those stored in the previous sliding time window in the system memory. When the current weight information entropy is detected to be greater than the weight information entropy of the previous period, indicating an increase in system flow distribution complexity or operating condition volatility, a preset basic quantile constant is added to a pre-configured fixed step increment. The built-in minimum value function is then called to find the minimum value between the added updated quantile value and the floating-point constant 1 to achieve a hard upper limit truncation constraint. Conversely, when the current weight information entropy is detected to be less than or equal to the weight information entropy of the previous period, the basic quantile constant is subtracted from the fixed step increment. A percentile statistical function is used to transform the basic quantile constant, after the above-mentioned rule update and correction, into a specific percentage statistical parameter. This percentage is then used to sort and extract the comprehensive anomaly score sequence of all historical normal industrial control communications persistently stored in the time-series database. The specific historical score value corresponding to this percentile is then updated and set as the anomaly judgment threshold for the current monitoring period. The comprehensive anomaly score of the current real-time input traffic sample is continuously judged. When the comprehensive anomaly score exceeds the judgment threshold, a system interrupt handling routine is triggered, and a message queue middleware is invoked to push a red anomaly alarm message log to the host computer software interface of the coal mine industrial control safety control center to complete the abnormal behavior alarm output. This threshold adjustment method is used to suppress false alarms caused by operating condition switching during traffic distribution fluctuations and to improve the detection sensitivity of subtle anomaly deviations when the traffic structure tends to stabilize.

[0053] In an optional embodiment, based on the information entropy of the weights of each traffic cluster during the current time window, the basic quantile constant is adjusted, and the historical comprehensive anomaly score corresponding to the new quantile is used as the judgment threshold. This includes: extracting the normalized initial weight probability distribution of all traffic clusters within the current time window; calculating the information entropy of the traffic cluster weights during the current time window according to the information entropy formula; calculating the difference rate between the current traffic cluster weight information entropy and the historical traffic cluster weight information entropy of the previous time window; if the difference rate is greater than 0, then the basic quantile constant is increased proportionally to the difference rate and truncated to an upper limit of 1 as the new quantile; if the difference rate is not greater than 0, then the basic quantile constant is reduced proportionally to the absolute value of the difference rate as the new quantile; and using the score value corresponding to the new quantile in the historical comprehensive anomaly score distribution as the judgment threshold.

[0054] Extract the normalized initial weight probability distribution of all traffic clusters within the current time window, and then apply the Shannon information entropy formula. ,in Let c be the normalized initial weights for the c-th traffic cluster; when When =0, this item is treated as 0. The entropy of traffic cluster weight information for the current time window is calculated. This metric can comprehensively represent the balance and dispersion of traffic cluster weight distribution within the current time window, and is used to reflect whether the weights of different service clusters are redistributed due to operating condition switching, equipment start-up and shutdown, or abnormal traffic disturbances.

[0055] The difference rate between the weight information entropy of the current time window and the weight information entropy of the previous historical time window is calculated. The difference rate is the difference between the weight information entropy of the current traffic cluster and the weight information entropy of the traffic cluster in the previous time window, divided by the sum of the absolute value of the weight information entropy of the traffic cluster in the previous time window and the smallest positive number. This is used to extract the phase transition coefficient representing network state fluctuations. If the difference rate is determined to be greater than zero, it indicates that the complexity of the current network traffic distribution or the volatility of the operating conditions has increased. In this case, a false alarm suppression mechanism is adopted. The preset basic quantile constant is amplified proportionally to the difference rate, and then truncated to a value no greater than 1 using the supremum function as the new quantile. This allows the basic quantile constant to be appropriately increased to raise the judgment threshold when the weight of the normal business cluster is redistributed due to changes in coal mine production scheduling, equipment start-up and shutdown, or communication cycle changes, thereby reducing false alarms caused by normal operating condition switching. If the difference rate is not greater than zero, it indicates that the current network state is stable. In this case, the basic quantile constant is reduced proportionally to the absolute value of the difference rate to enhance the response capability to slight abnormal deviations. Based on the calculated new quantiles, the corresponding quantile score value is searched in the historically accumulated comprehensive anomaly score distribution set. This extracted true score value serves as the alarm threshold for judging abnormal network behavior within the current time window, thereby adjusting the detection sensitivity. The process of adjusting the anomaly judgment threshold according to the complexity of network traffic distribution or operational fluctuations is as follows: Figure 3 As shown.

[0056] The test dataset used was industrial control protocol traffic captured from a real coal mine network. A total of 1000 packet samples were extracted, including 950 normal business polling packets and 50 simulated slow probe and abnormal burst packets. The experiment was conducted on a hardware platform with 64GB of RAM and a 16-core CPU. To verify the effectiveness of each improved module in this scheme, three control models were set up: the baseline model used the standard isolated forest algorithm; the ablation model removed the long-tail sparse cluster compensation strategy and the decision threshold adjustment module; and the complete model adopted a complete scheme including two-dimensional feature clustering, weighted sampling, deep adaptive isolated tree construction, and a threshold adjustment mechanism.

[0057] The same test sample set was input into the three models mentioned above for anomaly detection, and the precision, recall, F1 score, and false positive rate were statistically analyzed. Experimental results show that the baseline model achieved a precision of 83.2%, a recall of 79.4%, an F1 score of 81.2%, and a false positive rate of 6.5%; the ablation model achieved a precision of 90.5%, a recall of 85.6%, an F1 score of 88%, and a false positive rate of 3.2%; and the complete model achieved a precision of 96.8%, a recall of 97.1%, an F1 score of 96.9%, and a false positive rate of 0.8%. Figure 4 As shown, the results indicate that the modules in this scheme, such as sparse cluster compensation, weighted sampling, isolation tree depth adjustment, and dynamic threshold adjustment, can improve the abnormal traffic detection effect.

[0058] Furthermore, the basic share compensation mechanism of long-tail sparse micro-clusters helps reduce the risk of excessive dilution of mutated and abnormal data during weight allocation, thereby improving the sampling coverage of concealed long-cycle detection messages; the non-uniform segmentation strategy based on feature variance and the tree depth adjustment method based on covariance matrix trace help improve the isolation tree's ability to distinguish between normal homogeneous flow and abnormal deviation samples; by adjusting the judgment threshold in combination with the change of weight information entropy over time periods, the risk of false alarms can be reduced when coal mine production conditions fluctuate, and the sensitivity of anomaly detection can be improved when the flow structure is relatively stable, thereby improving the overall stability and reliability of the system monitoring.

[0059] In a second embodiment, the present invention also provides a coal mine industrial control network abnormal flow monitoring system, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement a coal mine industrial control network abnormal flow monitoring method according to the present invention.

[0060] The system also includes other components well known to those skilled in the art, such as communication buses and communication interfaces, the settings and functions of which are known in the art and will not be described in detail here.

[0061] In the description of this specification, "multiple" or "several" means at least two, such as two, three or more, unless otherwise expressly and specifically defined.

Claims

1. A method for monitoring abnormal traffic in a coal mine industrial control network, characterized in that, include: S1. Collect coal mine industrial control protocol flow data, calculate the flow packet payload byte entropy value and adjacent timestamp interval to construct a two-dimensional feature space and divide it into several non-intersecting flow clusters through clustering; S2. Use the product of the inverse of the sample density dispersion within each flow cluster and the Mahalanobis distance variance as the initial weight to determine the initial number of subsamples according to the proportion, and compensate sparse clusters with density dispersion higher than the judgment threshold to the preset basic share; S3. Distribute the extra over-allocation quota occupied by sparse clusters according to the proportion of the inverse of the square root of the Mahalanobis distance variance of the remaining flow clusters, and randomly and uniformly remove redundant samples from the deducted flow clusters to obtain each subsample set; S4. Construct an isolation tree for each subsample set, randomly select the split dimension according to the non-uniform probability set by the feature variance ratio, randomly determine the split point on the selected split dimension, calculate the trace of the two-dimensional covariance matrix, process it through standardization, map it through the S-type activation function, multiply it by the basic depth value, and round it up to obtain the maximum depth upper limit of the isolation tree; S5: Substitute the average path length of the traffic sample in each isolation tree into the evaluation function to calculate the single tree score and perform a weighted geometric average according to the traffic cluster weight to obtain the comprehensive anomaly score; S6: Based on the information entropy of the weight of each traffic cluster in the current time window, adjust the basic quantile constant, and use the historical comprehensive anomaly score of the corresponding new quantile as the judgment threshold. When the comprehensive anomaly score exceeds the judgment threshold, an anomaly alarm is output.

2. The method according to claim 1, characterized in that, The calculated traffic packet payload byte entropy value and adjacent timestamp intervals construct a two-dimensional feature space, which is then clustered into several non-overlapping traffic clusters, including: The frequency of each byte value in the payload of a single data packet is counted, and the payload byte entropy value is calculated according to the Shannon information entropy formula; the arrival time of the current traffic packet and the previous adjacent traffic packet are extracted, and the absolute value of the difference is calculated as the timestamp interval. The payload byte entropy value and the timestamp interval are standardized using Z-standard fraction processing and combined to form two-dimensional feature coordinate points; A density-based noise-based spatial clustering algorithm is used to iteratively divide the coordinate points in the two-dimensional feature space, forming several non-overlapping flow clusters.

3. The method according to claim 1, characterized in that, The method of determining the initial subsample number by using the product of the inverse of the sample density dispersion within each flow cluster and the variance of the Mahalanobis distance as the initial weight and according to the proportion includes: Calculate the two-dimensional covariance matrix and pseudo-inverse matrix of each flow cluster, obtain the Mahalanobis distance from each sample point in the flow cluster to the center of the flow cluster, and calculate the variance of the Mahalanobis distance in each flow cluster. Calculate the average Euclidean distance from each sample point within the flow cluster to its K nearest neighbor sample, and take the mean of the Euclidean distances as the density dispersion of the flow cluster. The absolute weight is obtained by multiplying the inverse of the density dispersion of the corresponding flow cluster by the variance of the Mahalanobis distance. The proportion of the absolute weight of each flow cluster to the total absolute weight of all flow clusters is calculated as the initial weight. The theoretical subsample number of each flow cluster is obtained by multiplying the total sampling scale constant by the initial weight of each flow cluster.

4. The method according to claim 3, characterized in that, The method of determining the initial subsample number by using the product of the inverse of the sample density dispersion within each flow cluster and the variance of the Mahalanobis distance as the initial weight and according to the proportion also includes: The theoretical subsample number is rounded down to obtain the temporary subsample number, and the difference between the total sampling size constant and the sum of all temporary subsample numbers is calculated as the remaining quota after rounding. The remaining quota after rounding is sequentially added to the corresponding flow clusters in descending order of the fractional part of the theoretical subsample number of each flow cluster to obtain the initial subsample number of each flow cluster, so that the sum of all the initial subsample numbers is equal to the total sampling scale constant.

5. The method according to claim 1, characterized in that, The compensation of sparse clusters with density dispersion higher than the judgment threshold to a preset base share includes: Extract the system's preset density dispersion judgment threshold and preset basic share constant; Traverse all traffic clusters and select those with a density dispersion greater than the judgment threshold to form a sparse cluster set; Calculate the difference between the current initial subsample number of the sparse cluster and the preset basic share constant. If the initial subsample number is lower than the preset basic share constant, set the subsample number of the sparse cluster to the preset basic share and record the total amount of additional over-allocation generated by the compensation.

6. The method according to claim 5, characterized in that, The process of allocating the excess over-allocation quota for sparse clusters according to the reciprocal of the square root of the variance of the Mahalanobis distance of each remaining flow cluster, and randomly and uniformly removing redundant samples from the deducted flow clusters to obtain each subset, includes: removing the sparse cluster set from all flow clusters to obtain the remaining flow cluster set; calculating the square root of the variance of the Mahalanobis distance of each flow cluster in the remaining flow cluster set, adding the square root to a minimum positive number and taking the reciprocal, and using the proportion of each reciprocal to the sum of reciprocals as the allocation coefficient; multiplying the total excess over-allocation quota by the allocation coefficient of each remaining flow cluster, rounding down to obtain the initial number of samples to be deducted for each remaining flow cluster, and allocating the difference between the total excess over-allocation quota and the sum of all initial sample numbers to the corresponding remaining flow clusters in descending order of the deduction quota to supplement the deduction quota, thus obtaining the deduction quantity for each remaining flow cluster.

7. The method according to claim 6, characterized in that, The process of obtaining each sub-sample set also includes: In the corresponding remaining traffic cluster, redundant samples of the same number are removed using a random uniform sampling method to obtain each sub-sample set. After removal, the sum of the sample counts of all sub-sample sets is verified to be equal to the preset total sampling scale constant.

8. The method according to claim 1, characterized in that, The construction of an isolation tree for each subset of sampling sets includes: Calculate the feature variance of the payload byte entropy value dimension and the timestamp interval dimension in the subsample set respectively. Use the ratio of the one-dimensional variance to the sum of the two-dimensional variances as the probability of selecting the dimension for splitting. Randomly select the splitting dimension according to the probability, and uniformly and randomly select the splitting point between the maximum and minimum values ​​of the dimension. The sum of the main diagonal elements of the two-dimensional covariance matrix of the subsample set is extracted as the trace of the matrix. After standardization, the trace of the matrix is ​​input into the S-type activation function to map the value to the interval 0 to 1. The result is multiplied by the base depth value, and the product is rounded up. The resulting integer is the maximum depth limit that the isolation tree is allowed to grow during construction.

9. The method according to claim 1, characterized in that, The process of adjusting the basic quantile constant based on the weighted information entropy of each traffic cluster within the current time window, and using the historical comprehensive anomaly score corresponding to the new quantile as the judgment threshold, includes: Extract the normalized initial weight probability distribution of all traffic clusters within the current time window, and calculate the traffic cluster weight information entropy for the current time window period according to the information entropy formula; calculate the difference rate between the current traffic cluster weight information entropy and the historical traffic cluster weight information entropy of the previous time window; if the difference rate is greater than 0, the basic quantile constant is amplified proportionally to the difference rate and truncated to an upper limit of 1 as the new quantile; if the difference rate is not greater than 0, the basic quantile constant is reduced proportionally to the absolute value of the difference rate as the new quantile; use the score value corresponding to the new quantile in the historical comprehensive anomaly score distribution as the judgment threshold.

10. A coal mine industrial control network abnormal flow monitoring system, characterized in that, include: A processor and a memory, wherein the memory stores computer program instructions that, when executed by the processor, implement the abnormal flow monitoring method for coal mine industrial control networks according to any one of claims 1-9.