Industrial computer network access control method for power monitoring system

Abstract

The application relates to the technical field of industrial network information security, and discloses an industrial computer network access control method for a power monitoring system, which comprises the following steps: obtaining the arrival time and the length of a target node packet, performing zero-order holding alignment on a reference node according to a time sequence to obtain a packet length sequence; calculating the flow intensity according to the packet length and the time difference, adding a zero constant to obtain a joint reference flow; reconstructing the three-dimensional phase space state vector of two nodes according to a delay step; obtaining a geometric centroid and calculating the Euclidean distance to obtain a topological offset; performing the division of the offset and the compensated reference flow to generate an adaptive calibration factor through negative exponential mapping; performing backward difference operation on the phase space vector, extracting the ratio of the outer product and the first-order difference constant cube, and multiplying the calibration factor to obtain a mapping sequence; calculating the mean and variance of the mapping sequence, and constructing a dynamic threshold structure which increases with the increase of the variance; performing the division of the latest sequence item and the threshold structure to obtain a decision index, and executing blocking or release according to the decision index.

Description

Technical Field

[0001] This invention relates to the field of industrial network information security technology, specifically to a method for network access control of industrial control computers for power monitoring systems. Background Technology

[0002] With the expansion of the scale of modern power grids and the advancement of intelligent processes, the network communication security of the internal industrial control computer of the power monitoring system, as the core hub for ensuring the operation of the power grid, is particularly important.

[0003] Currently, network access control for industrial control computers in power monitoring systems typically relies on static access control lists based on IP addresses or MAC addresses, along with various signature-based matching techniques. However, in real-world power monitoring scenarios, attackers often hijack authorized, legitimate industrial control computer nodes to forge identities and send malicious commands. Faced with such malicious injection operations masked by legitimate static identities, traditional static rules only verify surface communication identifiers and cannot identify the deep malicious logic hidden beneath the guise of legitimate communication. To compensate for the shortcomings of static rules, some existing technologies attempt to extract traffic features for dynamic anomaly detection. However, in the actual operation of power grids, normal power load fluctuations cause drastic nonlinear deformations in the communication traffic intensity sequence. It is extremely easy to confuse with real attack traffic, leading to high-frequency false alarms; furthermore, during the detection calculation process, the routing and forwarding of the underlying power switch will generate microsecond-level high-frequency network jitter, and the occasional reconstruction operation of the power network topology will cause the communication baseline of the reference node to drift randomly. These interferences from the underlying hardware and physical topology make conventional traffic detection extremely prone to local feature distortion; in addition, when the power system experiences a global transient failure or faces a network-wide short circuit, the communication behavior of all nodes exhibits sudden and violent fluctuations. Existing methods, which generally use manually set fixed security boundaries or empirical thresholds, will completely lose their adaptive adjustment capabilities under global transient impacts, resulting in the wrong interception of a large area of ​​legitimate sudden emergency traffic.

[0004] Therefore, there is an urgent need for a network access control method for industrial control computers in power monitoring systems to solve the above problems. Summary of the Invention

[0005] This invention provides a network access control method for industrial control computers in power monitoring systems, which helps to solve the problems mentioned in the background art.

[0006] This invention provides the following technical solution: a method for network access control of industrial control computers for power monitoring systems, comprising:

[0007] Obtain the arrival time and length of packets at the target node, and perform zero-order alignment maintenance on the reference nodes according to the time series to obtain the packet length sequence.

[0008] Calculate the flow intensity based on the packet length and the time difference of arrival, use the sum of the time difference and the zero constant as the denominator, and then calculate the baseline flow rate.

[0009] Historical flow intensity is extracted backward according to the set delay step, and the three-dimensional phase space state vectors of the two nodes are spliced ​​and reconstructed.

[0010] The geometric centroid is obtained by taking the mean value of the phase space state vector, and the topological offset is obtained by calculating the Euclidean distance between the centroids of the two nodes.

[0011] The adaptive calibration factor is generated by dividing the topology offset by the compensated baseline flow rate and mapping it using a negative exponential function.

[0012] Perform backward difference operation on the phase space vector, extract the ratio of the outer product to the first-order difference cube of the addition constant, and then multiply by the calibration factor to obtain the mapping sequence;

[0013] Calculate the mean and variance of the mapped sequence, and then add the mean and variance to the quotient of the mean plus a positive bias constant to construct a dynamic threshold structure.

[0014] The decision index is obtained by quoting the latest mapping sequence item with the threshold structure, and the communication port is blocked or allowed based on the index comparison result.

[0015] Optionally, obtaining the arrival time and length of the target node's packets, and performing zero-order alignment maintenance on the reference nodes according to the time series to obtain the packet length sequence, specifically includes:

[0016] Deploy a bypass mirroring agent on the port of the core switch, set the industrial control computer to be monitored as the target node, and specify the industrial control computer running in the same network segment as the reference node;

[0017] Record the network data packets of the target node within the sliding time window, and obtain the arrival timestamp and packet length;

[0018] Set adjacent arrival timestamps to strictly monotonically increasing;

[0019] Set the system's minimum time measurement resolution constant;

[0020] Set the time difference between the arrival timestamp of the current data packet and the previous data packet to be greater than or equal to the system's minimum time measurement resolution constant.

[0021] On the proxy side, zero-order alignment is performed on the data of the reference node. The arrival timestamp of the target node is used as the trigger time. The latest packet length of the reference node before the above time is taken as the alignment packet length to obtain the reference node packet length sequence with equal length alignment.

[0022] Optionally, the step of calculating the flow intensity based on the packet length and the time difference of arrival, using the sum of the time difference and the zero-prevention constant as the denominator, and then calculating the baseline flow rate specifically includes:

[0023] The traffic intensity of the target node and the reference node is calculated based on the packet length and arrival time difference of adjacent data packets, as follows:

[0024] The target node's packet length is used as the numerator, and the difference between the arrival timestamp of the current packet and the previous packet is added to the zero-prevention constant and used as the denominator. The two are then divided to obtain the zero-prevention regularized traffic intensity of the target node.

[0025] Using the same calculation logic, the zero-prevention regularization traffic intensity of the reference node is calculated based on the aligned reference node packet length, the target node time difference, and the aforementioned zero-prevention constant.

[0026] The traffic intensity of the target node and the reference node within the time window is summed and divided by twice the total number of data packets minus one. The joint arithmetic mean is then calculated to obtain the joint network baseline traffic.

[0027] Optionally, the step of extracting historical flow intensity backward according to a set delay step and concatenating and reconstructing the three-dimensional phase space state vectors of the two nodes specifically includes:

[0028] Set a fixed delay step constant, the value of which is the smallest positive integer;

[0029] Historical flow intensity is extracted backward from the time series using the aforementioned delay step constant;

[0030] The current flow intensity, the flow intensity delayed by one step, and the flow intensity delayed by two steps are concatenated to reconstruct the three-dimensional phase space state vector of the target node.

[0031] Using the same reconstruction logic and parallel calculation steps, the three-dimensional phase space state vector of the reference node is reconstructed.

[0032] Optionally, the step of calculating the geometric centroid of the mean value of the phase space state vector and obtaining the topological offset by calculating the Euclidean distance between the centroids of the two nodes specifically includes:

[0033] The three-dimensional phase space state vectors of all target nodes within the effective time window are summed, divided by the total number of effective vectors, and the average value is calculated to obtain the three-dimensional phase space geometric centroid vector of the target node.

[0034] The same calculation steps are used to obtain the three-dimensional phase space geometric centroid vector of the reference node;

[0035] Extract the three corresponding dimension components of the centroid vectors of the target node and the reference node in three-dimensional coordinates;

[0036] Calculate the difference between the three corresponding dimensional components, square the difference between each dimension and sum them, then take the square root of the sum to calculate the topological offset across nodes.

[0037] Optionally, the step of dividing the topology offset by the compensated reference flow rate and generating an adaptive calibration factor using a negative exponential function specifically includes:

[0038] Set a minimal positive real constant to prevent division by zero errors;

[0039] Add the calculated network joint baseline traffic to the above minimal positive real constant to obtain the joint baseline traffic after zero-compensation.

[0040] The cross-node topology offset calculated above is used as the numerator, and the joint reference flow after zero compensation is used as the denominator. The ratio is obtained by dividing the two.

[0041] The above ratio results are negativeized and used as the independent variable of the natural exponential function for mapping calculation, thus constructing a dimensionless adaptive calibration factor.

[0042] Optionally, the step of performing backward difference operation on the phase space vector, extracting the ratio of the outer product to the first-order difference cube of the addition constant, and then multiplying by the calibration factor to obtain the mapping sequence specifically includes:

[0043] Calculate the backward first-order difference vector and the backward second-order difference vector of the three-dimensional phase space vector of the target node;

[0044] Perform a three-dimensional vector outer product operation on the first-order difference vector and the second-order difference vector, and take the second norm of the outer product result as the numerator;

[0045] Find the L2 norm of the first-order difference vector, cube the norm, and add the set zero-prevention constant as the denominator.

[0046] Dividing the numerator by the denominator yields a geometric complexity index that reflects the phase trajectory change.

[0047] Cross-node correction is performed by multiplying the calculated adaptive calibration factor with the geometric complexity index to obtain the calibration mapping sequence.

[0048] Optionally, the calculation of the mean and variance of the mapped sequence, by adding the mean and variance to the quotient of the mean plus a positive bias constant, constructs a dynamic threshold structure, specifically including:

[0049] The arithmetic mean of the calibration mapping sequence is obtained by summing the terms of each term in the valid range within the time window and dividing by the total number of valid samples.

[0050] Calculate the sum of squares of the differences between each term in the sequence and the arithmetic mean, divide by the total number of valid samples, and obtain the statistical variance of the calibration mapping sequence;

[0051] Set a positive bias constant to prevent the denominator from being zero;

[0052] Using the above statistical variance as the numerator, add the arithmetic mean to the positive bias constant as the denominator, and divide the two to obtain the quotient.

[0053] By adding the arithmetic mean to the quotient, a dynamic threshold structure that increases with the variance is constructed.

[0054] Optionally, the step of quoting the latest mapping sequence item with the threshold structure to obtain a decision index, and then performing communication port blocking or allowing based on the index comparison result, specifically includes:

[0055] Obtain the latest calibration mapping sequence item corresponding to the end of the current sliding time window;

[0056] Using the latest calibration mapping sequence as the numerator and the dynamic threshold structure constructed above as the denominator, the two are divided to calculate the dimensionless control risk decision index.

[0057] Compare the risk control decision index with the value of one.

[0058] If the control risk decision index is greater than one, then the session connection of the corresponding port of the target industrial control computer will be cut off and a blocking action will be performed.

[0059] If the control risk decision index is less than or equal to one, then the corresponding communication port will be allowed to pass.

[0060] The present invention has the following beneficial effects:

[0061] 1. In a power monitoring system's industrial control computer network access control scenario, a continuous data processing chain is constructed based on the arrival time, packet length, and traffic intensity of network data packets between the target node and the reference node. First, the reference node is aligned to zero order using the target node's time series as a benchmark. Then, dual-node zero-prevention regularized traffic intensity is generated, and further, a three-dimensional phase space state vector is reconstructed. This transforms the originally dispersed underlying communication packet sequence into phase space trajectory data that reflects the correlation between the current state and adjacent historical states. Based on this, cross-node topology offsets are extracted using the spatial distance between the geometric centroids of the target and reference nodes. These topology offsets, along with the zero-prevention compensated network joint benchmark traffic, are used to generate an adaptive calibration factor, enabling subsequent local geometric reconstruction... The complexity can be corrected by combining the cross-node offset state within the current window; at the same time, a calibration mapping sequence is constructed through backward difference, outer product modulus, first-order difference modulus cube, and zero-prevention compensation, and the mean and variance of this sequence are used to form a dynamic threshold structure that changes with the degree of fluctuation within the window. Finally, the ratio of the latest mapping sequence item to the dynamic threshold structure is used as the control risk decision index, realizing closed-loop access control from the underlying packet-level data, phase space feature extraction, cross-node calibration, dynamic threshold construction to port blocking or allowing; this technology chain is beneficial to adjust the decision boundary according to the statistical state within the window when there is overall fluctuation in network traffic, and to form an executable access control output when the phase space trajectory of the target node deviates locally from the reference node and the dynamic threshold;

[0062] 2. By obtaining the arrival time and length of data packets at the target node, and using the target node's time series as a trigger reference, zero-order alignment is performed on reference nodes in the same network segment. This ensures that the packet length sequence of the reference nodes forms an equal-length aligned data sequence according to the arrival time of the target node. This processing method does not resynchronize the physical clocks of the two nodes, but rather extracts the packet length of the most recent data packet of the reference node before the corresponding trigger time under the target node's time reference. This provides a unified sequence position and time reference for subsequent traffic intensity calculations of both nodes. By setting adjacent arrival timestamps to be strictly monotonically increasing and requiring that the arrival time difference between the current data packet and the previous data packet is not less than the system's minimum time measurement resolution constant, clear time constraints are provided for subsequent packet interval calculations. Therefore, the target node packet length, the adjacent arrival time difference of the target node, the aligned packet length of the reference node, and the number of packets in the sliding time window obtained in step one can serve as the basic data source for traffic intensity calculation and subsequent window aggregation processing in step two.

[0063] 3. By using the packet length of the current data packet of the target node as the numerator, adding the arrival time difference between the current data packet and the previous data packet to the system's minimum time measurement resolution constant as the denominator, and performing division, the discrete packet length and packet arrival interval of the target node can be converted into zero-prevention regularized traffic intensity at the corresponding sequence position. Simultaneously, using the same processing logic, the aligned packet length of the reference node is combined with the adjacent arrival time difference under the target node's time reference to generate the zero-prevention regularized traffic intensity of the reference node under the same time reference, providing a comparable data basis for the traffic intensity of the two nodes. Within the time window, the traffic intensity of the target node and the reference node are summarized from the second data packet to the end of the window, and the average denominator is obtained by subtracting one from the total number of data packets and then multiplying by two. This yields the network joint reference traffic. This network joint reference traffic is composed of the traffic intensity of both nodes within the current window, providing a dynamic traffic reference for subsequent scaling of cross-node topology offsets.

[0064] 4. By setting a fixed delay step constant, and in this embodiment taking the delay step size as the smallest positive integer, the traffic intensity at the current sequence position, the historical traffic intensity at the first delay of the previous sequence position, and the historical traffic intensity at the second delay of the sequence position before that can be concatenated in a fixed order, thereby reconstructing the three-dimensional phase space state vectors of the target node and the reference node respectively. Compared with using only a single current traffic intensity, this processing method preserves the sequence correlation between the current state and adjacent historical states, enabling subsequent steps to observe the trajectory changes of node communication traffic in the three-dimensional phase space. Since this phase space reconstruction only calls the historical traffic intensity at the current sequence position and before it, without introducing data from future moments, it can maintain the causal order of data in real-time access control scenarios. At the same time, the target node and the reference node use the same delay extraction and concatenation logic, which is beneficial to maintaining the consistency of the two-node phase space representation in subsequent centroid calculation, topology offset extraction, and cross-node calibration processing.

[0065] 5. By summing the three-dimensional phase space state vectors of the target node from the third sequence position to the end of the window within the effective time window and dividing by the total number of effective vectors, the three-dimensional phase space geometric centroid of the target node within the current sliding time window can be obtained. After obtaining the three-dimensional phase space geometric centroid of the reference node using the same calculation steps, the component differences between the two centroid vectors in the three corresponding dimensions are extracted, and then squared, summed, and square-rooted sequentially to obtain the cross-node topological offset. This cross-node topological offset reflects the overall deviation of the phase space center position of the target node from the phase space center position of the reference node, rather than the physical spatial distance. By aggregating multiple phase space state vectors into a geometric centroid, the influence of instantaneous fluctuations in a single sequence position on the overall phase space center position can be reduced, and the difference between the phase space centers of the two nodes can be converted into a distance-type feature that can be used for subsequent negative exponential mapping.

[0066] 6. By setting a minimal positive flow compensation constant at the same flow intensity scale as the network joint reference flow and adding it to the network joint reference flow, a joint reference flow after zero compensation can be constructed. Based on this, by using the cross-node topology offset as the numerator and the joint reference flow after zero compensation as the denominator, the ratio of the cross-node topology offset to the current window joint flow scale can be obtained. Subsequently, the negative value of this ratio is taken and input into the natural exponential function for mapping, which can generate a dimensionless adaptive calibration factor. This adaptive calibration factor does not change the preceding flow intensity and phase space state vector itself, but participates in the processing as a multiplicative correction when the subsequent geometric complexity index forms a calibration mapping sequence, so that the local phase trajectory geometric complexity can be corrected by combining the cross-node topology offset state and the current joint flow scale, thereby forming a mapping sequence basis that is related to both single-node local changes and double-node overall offsets.

[0067] 7. By performing backward first-order difference processing on the three-dimensional phase space state vector of the target node, the direction and magnitude of state changes between adjacent sequence positions can be obtained. By adding the three-dimensional phase space state vector of the current sequence position to the three-dimensional phase space state vector of the previous sequence position, and then subtracting twice the three-dimensional phase space state vector of the previous sequence position, a backward second-order difference vector can be obtained, which is used to characterize the second-order discrete change of the phase space state change trend. Based on this, a three-dimensional vector outer product is performed on the backward first-order difference vector and the backward second-order difference vector. The second norm of the outer product result is used as the numerator, and the cube of the second norm of the backward first-order difference vector is added to the zero constant of the same scale as the denominator. This yields a geometric complexity index reflecting the degree of local curvature of the phase space trajectory of the target node. Multiplying this geometric complexity index by the adaptive calibration factor forms a calibration mapping sequence term after cross-node correction, so that the output of step six simultaneously reflects the local phase trajectory change of the target node and the cross-node calibration relationship generated in step five.

[0068] 8. By calling the calibration mapping sequence generated in step six, and using the implementation method of taking the smallest positive integer for the delay step, the calibration mapping sequence items from the fifth sequence position to the end of the sliding time window are taken as statistical objects, which can clarify the effective sample range for the calculation of mean and variance. The arithmetic mean of the calibration mapping sequence can be obtained by summing each calibration mapping sequence item in the statistical object and dividing it by the total number of effective samples. The difference between each effective calibration mapping sequence item and the arithmetic mean can be further calculated, the differences can be squared and summarized, and then divided by the total number of effective samples to obtain the statistical variance of the calibration mapping sequence. The statistical variance is used as the numerator, the arithmetic mean is added to the positive bias constant of the same scale as the denominator, and the two are divided to obtain the variance normalization quotient. The arithmetic mean is then added to the variance normalization quotient to construct a dynamic threshold structure that varies with the dispersion of the calibration mapping sequence. This dynamic threshold structure is generated by the statistical results of the mapping sequence in the current window, which enables the access control decision boundary to be adjusted in a data-driven manner according to the fluctuation level in the window, rather than relying solely on a fixed static boundary.

[0069] 9. By obtaining the latest calibration mapping sequence item corresponding to the end of the current sliding time window and using it as the numerator, and the dynamic threshold structure generated in step seven as the denominator, a quotient operation is performed to obtain a dimensionless control risk decision index. This control risk decision index places the local calibration mapping result at the end of the current window and the dynamic threshold structure formed by the current window statistics under the same judgment scale, so that the complex phase space trajectory changes, cross-node calibration relationships and window statistical boundaries are uniformly converted into directly comparable access control indicators. By comparing the control risk decision index with the value one, clear blocking or allowing conditions can be formed. When the control risk decision index is greater than the value one, the blocking action of the session connection of the corresponding port of the target industrial control computer is executed. When the control risk decision index is less than or equal to the value one, the allowing operation is performed on the corresponding communication port. Thus, step eight transforms the results of the preceding data acquisition, feature construction, calibration mapping and dynamic threshold calculation into clear network access control output. Attached Figure Description

[0070] Figure 1 This is a schematic diagram of the overall process of the present invention.

[0071] Figure 2 This is a schematic diagram of the deployment of monitoring nodes and network topology of the system of the present invention.

[0072] Figure 3 This is a schematic diagram of the computational space geometry for 3D phase space reconstruction and cross-node topology offset.

[0073] Figure 4 The graph shows the dynamic threshold and risk decision logic curves that monotonically expand with variance. Detailed Implementation

[0074] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0075] Example 1, refer to Figure 1 A method for network access control of industrial control computers for power monitoring systems, comprising:

[0076] Obtain the arrival time and length of packets at the target node, and perform zero-order alignment maintenance on the reference nodes according to the time series to obtain the packet length sequence.

[0077] Calculate the flow intensity based on the packet length and the time difference of arrival, use the sum of the time difference and the zero constant as the denominator, and then calculate the baseline flow rate.

[0078] Historical flow intensity is extracted backward according to the set delay step, and the three-dimensional phase space state vectors of the two nodes are spliced ​​and reconstructed.

[0079] The geometric centroid is obtained by taking the mean value of the phase space state vector, and the topological offset is obtained by calculating the Euclidean distance between the centroids of the two nodes.

[0080] The adaptive calibration factor is generated by dividing the topology offset by the compensated baseline flow rate and mapping it using a negative exponential function.

[0081] Perform backward difference operation on the phase space vector, extract the ratio of the outer product to the first-order difference cube of the addition constant, and then multiply by the calibration factor to obtain the mapping sequence;

[0082] Calculate the mean and variance of the mapped sequence, and then add the mean and variance to the quotient of the mean plus a positive bias constant to construct a dynamic threshold structure.

[0083] The decision index is obtained by quoting the latest mapping sequence item with the threshold structure, and the communication port is blocked or allowed based on the index comparison result.

[0084] The process of obtaining the arrival time and length of packets from the target node, and performing zero-order alignment maintenance on the reference nodes according to the time series to obtain the packet length sequence, specifically includes:

[0085] Deploy a bypass mirroring agent on the port of the core switch, set the industrial control computer to be monitored as the target node, and specify the industrial control computer running in the same network segment as the reference node;

[0086] Record the network data packets of the target node within the sliding time window, and obtain the arrival timestamp and packet length;

[0087] Set adjacent arrival timestamps to strictly monotonically increasing;

[0088] Set the system's minimum time measurement resolution constant;

[0089] Set the time difference between the arrival timestamp of the current data packet and the previous data packet to be greater than or equal to the system's minimum time measurement resolution constant.

[0090] On the proxy side, zero-order alignment is performed on the data of the reference node. The arrival timestamp of the target node is used as the trigger time. The latest packet length of the reference node before the above time is taken as the alignment packet length to obtain the reference node packet length sequence with equal length alignment.

[0091] The calculation of flow intensity based on packet length and arrival time difference, using the sum of time difference and zero-resistance constant as the denominator, and then determining the baseline flow rate, specifically includes:

[0092] The traffic intensity of the target node and the reference node is calculated based on the packet length and arrival time difference of adjacent data packets, as follows:

[0093] The target node's packet length is used as the numerator, and the difference between the arrival timestamp of the current packet and the previous packet is added to the zero-prevention constant and used as the denominator. The two are then divided to obtain the zero-prevention regularized traffic intensity of the target node.

[0094] Using the same calculation logic, the zero-prevention regularization traffic intensity of the reference node is calculated based on the aligned reference node packet length, the target node time difference, and the aforementioned zero-prevention constant.

[0095] The traffic intensity of the target node and the reference node within the time window is summed and divided by twice the total number of data packets minus one. The joint arithmetic mean is then calculated to obtain the joint network baseline traffic.

[0096] The step of extracting historical flow intensity backward according to a set delay step and concatenating and reconstructing the three-dimensional phase space state vectors of the two nodes specifically includes:

[0097] Set a fixed delay step constant, the value of which is the smallest positive integer;

[0098] Historical flow intensity is extracted backward from the time series using the aforementioned delay step constant;

[0099] The current flow intensity, the flow intensity delayed by one step, and the flow intensity delayed by two steps are concatenated to reconstruct the three-dimensional phase space state vector of the target node.

[0100] Using the same reconstruction logic and parallel calculation steps, the three-dimensional phase space state vector of the reference node is reconstructed.

[0101] The geometric centroid of the mean value of the phase space state vector is calculated, and the topological offset is obtained by calculating the Euclidean distance between the centroids of the two nodes. Specifically, this includes:

[0102] The three-dimensional phase space state vectors of all target nodes within the effective time window are summed, divided by the total number of effective vectors, and the average value is calculated to obtain the three-dimensional phase space geometric centroid vector of the target node.

[0103] The same calculation steps are used to obtain the three-dimensional phase space geometric centroid vector of the reference node;

[0104] Extract the three corresponding dimension components of the centroid vectors of the target node and the reference node in three-dimensional coordinates;

[0105] Calculate the difference between the three corresponding dimensional components, square the difference between each dimension and sum them, then take the square root of the sum to calculate the topological offset across nodes.

[0106] The step of dividing the topology offset by the compensated baseline flow rate and generating an adaptive calibration factor using a negative exponential function specifically includes:

[0107] Set a minimal positive real constant to prevent division by zero errors;

[0108] Add the calculated network joint baseline traffic to the above minimal positive real constant to obtain the joint baseline traffic after zero-compensation.

[0109] The cross-node topology offset calculated above is used as the numerator, and the joint reference flow after zero compensation is used as the denominator. The ratio is obtained by dividing the two.

[0110] The above ratio results are negativeized and used as the independent variable of the natural exponential function for mapping calculation, thus constructing a dimensionless adaptive calibration factor.

[0111] The step of performing backward difference operations on the phase space vectors, extracting the ratio of the outer product to the first-order difference cube of the addition constant, and then multiplying by the calibration factor to obtain the mapping sequence, specifically includes:

[0112] Calculate the backward first-order difference vector and the backward second-order difference vector of the three-dimensional phase space vector of the target node;

[0113] Perform a three-dimensional vector outer product operation on the first-order difference vector and the second-order difference vector, and take the second norm of the outer product result as the numerator;

[0114] Find the L2 norm of the first-order difference vector, cube the norm, and add the set zero-prevention constant as the denominator.

[0115] Dividing the numerator by the denominator yields a geometric complexity index that reflects the phase trajectory change.

[0116] Cross-node correction is performed by multiplying the calculated adaptive calibration factor with the geometric complexity index to obtain the calibration mapping sequence.

[0117] The calculation of the mean and variance of the mapped sequence involves adding the mean and variance to the quotient of the mean plus a positive bias constant to construct a dynamic threshold structure, specifically including:

[0118] The arithmetic mean of the calibration mapping sequence is obtained by summing the terms of each term in the valid range within the time window and dividing by the total number of valid samples.

[0119] Calculate the sum of squares of the differences between each term in the sequence and the arithmetic mean, divide by the total number of valid samples, and obtain the statistical variance of the calibration mapping sequence;

[0120] Set a positive bias constant to prevent the denominator from being zero;

[0121] Using the above statistical variance as the numerator, add the arithmetic mean to the positive bias constant as the denominator, and divide the two to obtain the quotient.

[0122] By adding the arithmetic mean to the quotient, a dynamic threshold structure that increases with the variance is constructed.

[0123] The step of quoting the latest mapping sequence item with the threshold structure to obtain a decision index, and then executing communication port blocking or allowing based on the index comparison result, specifically includes:

[0124] Obtain the latest calibration mapping sequence item corresponding to the end of the current sliding time window;

[0125] Using the latest calibration mapping sequence as the numerator and the dynamic threshold structure constructed above as the denominator, the two are divided to calculate the dimensionless control risk decision index.

[0126] Compare the risk control decision index with the value of one.

[0127] If the control risk decision index is greater than one, then the session connection of the corresponding port of the target industrial control computer will be cut off and a blocking action will be performed.

[0128] If the control risk decision index is less than or equal to one, then the corresponding communication port will be allowed to pass.

[0129] This paper addresses the network access control scenario of industrial control computers in power monitoring systems. A continuous data processing chain is constructed based on the arrival time, packet length, and traffic intensity of network data packets between target and reference nodes. First, zero-order alignment is performed on the reference node using the target node's time series as a benchmark. Then, dual-node zero-prevention regularized traffic intensity is generated, and a three-dimensional phase space state vector is further reconstructed. This transforms the originally dispersed underlying communication packet sequence into phase space trajectory data that reflects the correlation between the current state and adjacent historical states. Based on this, cross-node topology offsets are extracted using the spatial distance between the geometric centroids of the target and reference nodes. These topology offsets, along with the zero-prevention compensated network joint benchmark traffic, generate an adaptive calibration factor to address subsequent local geometrical complexities. The degree can be corrected by combining the cross-node offset state within the current window; at the same time, a calibration mapping sequence is constructed through backward difference, outer product modulus, first-order difference modulus cube, and zero-prevention compensation, and the mean and variance of this sequence are used to form a dynamic threshold structure that changes with the degree of fluctuation within the window. Finally, the ratio of the latest mapping sequence item to the dynamic threshold structure is used as the control risk decision index, realizing closed-loop access control from the underlying packet-level data, phase space feature extraction, cross-node calibration, dynamic threshold construction to port blocking or allowing; this technology chain is beneficial to adjust the judgment boundary according to the statistical state within the window when there is overall fluctuation in network traffic, and to form an executable access control output when the phase space trajectory of the target node deviates locally from the reference node and the dynamic threshold;

[0130] Example 2: A network access control method for industrial control computers in a power monitoring system, further comprising:

[0131] Step 1: Obtain the arrival time and length of the packets at the target node. Perform zero-order alignment on the reference nodes according to the time series to obtain the packet length sequence. This includes the following data acquisition and alignment operations.

[0132] like Figure 2 As shown, a bypass mirroring agent is first deployed on the port of the core switch. The industrial control computer to be monitored is set as the target node, and the industrial control computer running in the same network segment is specified as the reference node. The bypass mirroring agent records the network data packets of the target node within the sliding time window and obtains the arrival timestamp and packet length of each data packet.

[0133] Strictly monotonically increasing constraints are set for the adjacent arrival timestamps of the target node, and a system minimum time measurement resolution constant is set. The arrival timestamp of the current data packet is subtracted from the arrival timestamp of the previous data packet to obtain the arrival time difference between adjacent data packets, and this time difference is required to be no less than the system minimum time measurement resolution constant.

[0134] The minimum time measurement resolution constant of the system is a manually set basic time parameter, and its value range is recommended to be from one millionth of a second to one thousandth of a second. The value is based on the measurement accuracy limit of the underlying hardware clock. If the value is too large, it will cause excessive interception of normal high-frequency communication packets and result in feature loss. If the value is too small, it will not be able to effectively handle the risk of division by zero in the underlying calculation under concurrent burst traffic.

[0135] On the bypass mirror agent, zero-order alignment is performed on the data of the reference node. The arrival timestamp of each data packet of the target node is used as the trigger time. The most recent data packet of the reference node before the trigger time is found, and the packet length of the reference node data packet is truncated to the reference node aligned packet length corresponding to the current trigger time. This forms a reference node packet length sequence that is consistent with the number of data packets of the target node. At the same time, the total number of data packets in the sliding time window is recorded.

[0136] The total number of data packets within the sliding time window is a manually set data volume parameter, and its value range is recommended to be 100 to 1000. The value is based on the balance between statistical significance and real-time delay. If the value is too large, it will cause the system to lag in responding to sudden attacks. If the value is too small, it will lead to insufficient samples and cause baseline jitter and false alarms.

[0137] By acquiring the arrival time and length of data packets at the target node, and using the target node's time series as a trigger reference, zero-order alignment is performed on reference nodes in the same network segment. This ensures that the reference node's packet length sequence forms an equal-length aligned data sequence according to the target node's arrival time. This processing method does not resynchronize the physical clocks of the two nodes, but rather extracts the packet length of the most recent data packet from the reference node before the corresponding trigger time, based on the target node's time reference. This provides a unified sequence position and time reference for subsequent traffic intensity calculations of both nodes. By setting adjacent arrival timestamps to be strictly monotonically increasing and requiring that the arrival time difference between the current data packet and the previous data packet is not less than the system's minimum time measurement resolution constant, clear time constraints are provided for subsequent packet interval calculations. Thus, the target node packet length, the adjacent arrival time difference of the target node, the aligned packet length of the reference node, and the number of packets in the sliding time window formed in step one can serve as the basic data source for traffic intensity calculation in step two and subsequent window aggregation processing.

[0138] Step two: Calculate the flow intensity based on the packet length and the time difference of arrival. Use the sum of the time difference and the zero-prevention constant as the denominator, and calculate the baseline flow rate. This specifically includes:

[0139] The target node packet length, the time difference between adjacent arrivals of the target node, the aligned packet length of the reference node, and the system minimum time measurement resolution constant obtained in step one are used to calculate the zero-preservation traffic intensity of the target node and the reference node respectively.

[0140] When calculating the zero-prevention regularized traffic intensity of the target node, the packet length of the current data packet of the target node is used as the numerator, and the arrival time difference between the current data packet and the previous data packet is added to the system minimum time measurement resolution constant and used as the denominator. The numerator is divided by the denominator to obtain the zero-prevention regularized traffic intensity corresponding to the current data packet of the target node. This traffic intensity is used to characterize the data transmission intensity formed by the target node in the adjacent packet interval.

[0141] When calculating the zero-prevention regularized traffic strength of the reference node, the same calculation logic as that of the target node is used. The reference node aligned packet length obtained in step one is used as the numerator, and the arrival time difference between the current data packet and the previous data packet of the target node is added to the system minimum time measurement resolution constant and used as the denominator. The numerator is divided by the denominator to obtain the zero-prevention regularized traffic strength of the reference node under the time reference of the target node.

[0142] After completing the target node traffic intensity and reference node traffic intensity corresponding to the second data packet to the end of the sliding time window within the time window, the target node traffic intensity and reference node traffic intensity at each valid sequence position are added together, and then all the sums are summed. The total number of data packets in the sliding time window is reduced by 1 and then multiplied by 2 to obtain the average denominator. The total summation result is divided by the average denominator to obtain the network joint reference traffic in the current sliding time window. This network joint reference traffic serves as the basic traffic scale for subsequent cross-node calibration processing.

[0143] By using the packet length of the current data packet at the target node as the numerator, adding the arrival time difference between the current data packet and the previous data packet to the system's minimum time measurement resolution constant as the denominator, and performing division, the discrete packet length and packet arrival interval of the target node can be converted into zero-prevention regularized traffic intensity at the corresponding sequence position. Simultaneously, using the same processing logic, the reference node's aligned packet length is combined with the adjacent arrival time difference under the target node's time reference to generate the reference node's zero-prevention regularized traffic intensity under the same time reference, providing a comparable data basis for the dual-node traffic intensity. Within the time window, the traffic intensity of the target node and the reference node are summarized from the second data packet to the end of the window, and the average denominator is obtained by subtracting one from the total number of data packets and then multiplying by two. This yields the network joint reference traffic. This network joint reference traffic is composed of the dual-node traffic intensity within the current window, providing a dynamic traffic reference for subsequent scaling of cross-node topology offsets.

[0144] Step 3: Extract historical flow intensity backward according to the set delay step, and reconstruct the three-dimensional phase space state vectors of the two nodes. This specifically includes:

[0145] A fixed delay step constant is set. In this embodiment, the delay step constant is taken as the smallest positive integer. Based on the delay step, the current flow intensity is extracted at the current sequence position, the delayed historical flow intensity is extracted once at the previous sequence position, and the delayed historical flow intensity is extracted twice at the previous sequence position.

[0146] Based on the target node zero-prevention regularized flow intensity sequence generated in step two, the current flow intensity is extracted at the current sequence position, and then the historical flow intensity at a position delayed by one step is extracted according to the set delay step, and the historical flow intensity at a position delayed by two steps is extracted. The above three flow intensities are spliced ​​together in the order of current term, first delay term and second delay term to form the three-dimensional phase space state vector of the target node at the current sequence position.

[0147] Using the same delay extraction and splicing logic, based on the reference node zero-prevention regularized flow intensity sequence generated in step two, a three-dimensional phase space state vector of the reference node at the corresponding sequence position is formed. The three-dimensional phase space state vectors of the target node and the reference node are used to characterize the phase trajectory state formed by the corresponding node under the combined action of the current flow state and its adjacent historical states.

[0148] By setting a fixed delay step constant, and in this embodiment taking the delay step size as the smallest positive integer, the traffic intensity at the current sequence position, the historical traffic intensity at the first delay of the previous sequence position, and the historical traffic intensity at the second delay of the sequence position before that can be concatenated in a fixed order, thus reconstructing the three-dimensional phase space state vectors of the target node and the reference node respectively. Compared with using only a single current traffic intensity, this processing method preserves the sequence correlation between the current state and adjacent historical states, allowing subsequent steps to observe the trajectory changes of node communication traffic in the three-dimensional phase space. Since this phase space reconstruction only calls the historical traffic intensity at the current sequence position and before it, without introducing data from future moments, it can maintain the causal order of data in real-time access control scenarios. At the same time, the target node and the reference node use the same delay extraction and concatenation logic, which is beneficial to maintaining the consistency of the two-node phase space representation in subsequent centroid calculation, topology offset extraction, and cross-node calibration processing.

[0149] Step four involves finding the geometric centroid of the mean state vector in phase space and calculating the Euclidean distance between the centroids of the two nodes to obtain the topological offset. This includes:

[0150] The target node three-dimensional phase space state vector sequence generated in step three is called. Under the implementation method of taking the minimum positive integer for the delay step size, the target node three-dimensional phase space state vector from the third sequence position to the end of the sliding time window all participate in the centroid calculation.

[0151] The three-dimensional phase space state vectors of all target nodes involved in the centroid calculation are summed, and the summation result is divided by the total number of effective vectors to obtain the three-dimensional phase space geometric centroid vector of the target node. The total number of effective vectors is the total number of data packets in the sliding time window minus 2. The three-dimensional phase space geometric centroid vector of the target node is used to characterize the phase space center position of the target node in the current time window.

[0152] Using the same calculation steps, the three-dimensional phase space state vector sequence of the reference node generated in step three is called. Starting from the third sequence position to the end of the sliding time window, the three-dimensional phase space state vector of the reference node is extracted. All the extracted three-dimensional phase space state vectors of the reference node are summed and divided by the same total number of effective vectors to obtain the three-dimensional phase space geometric centroid vector of the reference node.

[0153] Reference Figure 3 Extract the three corresponding dimensional components of the centroid vector of the target node and the centroid vector of the reference node in three-dimensional coordinates, calculate the difference between each corresponding dimensional component, square the difference of each dimension and sum them, and then take the square root of the sum to obtain the cross-node topological offset. This cross-node topological offset is used to characterize the overall deviation of the phase space center position of the target node from the phase space center position of the reference node.

[0154] By summing the three-dimensional phase space state vectors of the target node from the third sequence position to the end of the window within the effective time window and dividing by the total number of effective vectors, the three-dimensional phase space geometric centroid of the target node within the current sliding time window can be obtained. After obtaining the three-dimensional phase space geometric centroid of the reference node using the same calculation steps, the component differences between the two centroid vectors in the three corresponding dimensions are extracted and then squared, summed, and square-rooted in sequence to obtain the cross-node topological offset. This cross-node topological offset reflects the overall deviation of the phase space center position of the target node from the phase space center position of the reference node, rather than the physical spatial distance. By aggregating multiple phase space state vectors into a geometric centroid, the influence of instantaneous fluctuations in a single sequence position on the overall phase space center position can be reduced, and the difference between the phase space centers of the two nodes can be converted into a distance-type feature that can be used for subsequent negative exponential mapping.

[0155] Step 5: Divide the topology offset by the compensated baseline flow rate, and use a negative exponential function to generate an adaptive calibration factor. This includes:

[0156] Set a minimal positive flow compensation constant to prevent the network joint baseline flow from being zero, and call the network joint baseline flow generated in step two and the cross-node topology offset generated in step four.

[0157] Add the network joint reference flow to the minimum positive flow compensation constant to obtain the joint reference flow after zero compensation. Use the cross-node topology offset as the numerator and the joint reference flow after zero compensation as the denominator. Divide the numerator by the denominator to obtain the ratio of the cross-node topology offset to the joint reference flow scale.

[0158] The ratio result is negative, and the negative result is used as the input of the natural exponential function for mapping to obtain a dimensionless adaptive calibration factor. The adaptive calibration factor is used to participate in the cross-node correction of the subsequent geometric complexity index.

[0159] The minimum positive flow compensation constant is a manually set basic flow compensation parameter. It is on the same flow intensity scale as the network joint reference flow and participates in the calculation as the zero-prevention compensation amount of the joint reference flow. The minimum positive flow compensation constant is used to maintain the computable state of exponential mapping processing when the network joint reference flow is extremely weak.

[0160] By setting a minimal positive flow compensation constant at the same flow intensity scale as the network joint reference flow and adding it to the network joint reference flow, a joint reference flow with zero compensation can be constructed. Based on this, by using the cross-node topology offset as the numerator and the joint reference flow with zero compensation as the denominator, the ratio of the cross-node topology offset to the current window joint flow scale can be obtained. Subsequently, the negative value of this ratio is taken and input into the natural exponential function for mapping, which can generate a dimensionless adaptive calibration factor. This adaptive calibration factor does not change the preceding flow intensity and phase space state vector itself, but participates in the processing as a multiplicative correction when the subsequent geometric complexity index forms a calibration mapping sequence. This allows the local phase trajectory geometric complexity to be corrected by combining the cross-node topology offset state and the current joint flow scale, thereby forming a mapping sequence basis that is related to both single-node local changes and double-node overall offsets.

[0161] Step six involves performing backward difference operations on the phase space vectors, extracting the ratio of the outer product to the first-order difference cube of the addition constant, and then multiplying by the calibration factor to obtain the mapping sequence. Specifically, this includes:

[0162] The target node's three-dimensional phase space state vector sequence generated in step three and the adaptive calibration factor generated in step five are used to perform backward differential processing on the three-dimensional phase space state vector at the current sequence position of the target node.

[0163] In the backward first-order difference processing, the three-dimensional phase space state vector of the target node at the current sequence position is subtracted from the three-dimensional phase space state vector of the target node at the previous sequence position to obtain the backward first-order difference vector. This backward first-order difference vector is used to characterize the direction and magnitude of the change in phase space state between adjacent sequence positions.

[0164] In the backward second-order difference processing, the three-dimensional phase space state vector of the target node at the current sequence position is added to the three-dimensional phase space state vector of the target node at the previous sequence position to obtain the first and last state superposition result. Then, the three-dimensional phase space state vector of the target node at the previous sequence position is doubled to obtain the intermediate state multiplication result. The first and last state superposition result is subtracted from the intermediate state multiplication result to obtain the backward second-order difference vector. This backward second-order difference vector is used to characterize the second-order discrete change of the phase space state change trend.

[0165] Perform a three-dimensional vector outer product operation on the backward first-order difference vector and the backward second-order difference vector, and calculate the second-norm modulus of the outer product result. Use this modulus as the numerator of the geometric complexity calculation. At the same time, calculate the second-norm modulus of the backward first-order difference vector, cube the first-order difference modulus, and add the cubed result to the set zero-prevention constant. Use the sum as the denominator of the geometric complexity calculation.

[0166] Dividing the numerator by the denominator yields a geometric complexity index that reflects the phase trajectory change. This geometric complexity index is used to characterize the degree of local bending change of the phase space trajectory of the target node at the current sequence position.

[0167] The zero-prevention constant is a basic calculation parameter set manually. It is used as a compensation quantity of the same scale as the result of the first-order difference vector second norm modulus long cube processing to participate in the denominator construction. It is recommended that its value range be between one part per billion and one part per million. Its value is based on the machine calculation precision limit of the system's floating-point operation. If the value is too large, it will excessively suppress the curvature characteristics of the normal trajectory and cause the algorithm to lose sensitivity. If the value is too small, it will be unable to handle the risk of calculation division by zero overflow when the network traffic is completely still and the phase space states overlap.

[0168] After obtaining the geometric complexity index, the adaptive calibration factor generated in step five is multiplied by the geometric complexity index to obtain the calibration mapping sequence item after cross-node correction, and the calibration mapping sequence is continuously generated according to the effective sequence position from the fifth sequence position to the end of the sliding time window.

[0169] By performing backward first-order difference processing on the three-dimensional phase space state vector of the target node, the direction and magnitude of state changes between adjacent sequence positions can be obtained. By adding the three-dimensional phase space state vector of the current sequence position to the three-dimensional phase space state vector of the previous sequence position, and then subtracting twice the three-dimensional phase space state vector of the previous sequence position, a backward second-order difference vector can be obtained, which is used to characterize the second-order discrete change of the phase space state change trend. Based on this, a three-dimensional vector outer product is performed on the backward first-order difference vector and the backward second-order difference vector. The second norm of the outer product result is used as the numerator, and the cube result of the second norm of the backward first-order difference vector is added to the zero constant of the same scale as the denominator. This yields a geometric complexity index reflecting the degree of local curvature of the phase space trajectory of the target node. Multiplying this geometric complexity index by the adaptive calibration factor forms a calibration mapping sequence term after cross-node correction, so that the output of step six simultaneously reflects the local phase trajectory change of the target node and the cross-node calibration relationship generated in step five.

[0170] Step 7: Calculate the mean and variance of the mapped sequence. Add the mean and variance to the quotient of the mean plus a positive bias constant to construct the dynamic threshold structure. This specifically includes:

[0171] The calibration mapping sequence generated in step six is ​​called, and under the implementation method of taking the smallest positive integer for the delay step, the calibration mapping sequence items from the fifth sequence position to the end of the sliding time window are used as statistical objects.

[0172] The arithmetic mean of the calibration mapping sequence is obtained by summing the values ​​of each calibration mapping sequence in the above statistical objects and dividing the sum by the total number of valid samples. The total number of valid samples is the total number of data packets in the sliding time window minus 4.

[0173] After obtaining the arithmetic mean of the calibration mapping sequence, the difference between each valid calibration mapping sequence item and the arithmetic mean is calculated. Each difference is squared, and all the squared results are summed. The summation result is then divided by the total number of valid samples to obtain the statistical variance of the calibration mapping sequence.

[0174] A positive bias constant is set to prevent the denominator from being zero and to be on the same computational scale as the arithmetic mean of the calibration mapping sequence. The statistical variance is used as the numerator, and the arithmetic mean is added to the positive bias constant as the denominator. The numerator is divided by the denominator to obtain the variance normalization quotient. The arithmetic mean is then added to the variance normalization quotient to construct a dynamic threshold structure that varies with the degree of dispersion of the calibration mapping sequence.

[0175] The positive bias constant is a manually set basic threshold parameter. It is used as a bias of the same scale as the arithmetic mean of the calibration mapping sequence to participate in the denominator construction of the dynamic threshold structure. Its value range is recommended to be one ten-thousandth to one hundredth. Its value is based on the basic noise curvature level under the resting steady state of the normal power network. If the value is too large, the dynamic threshold structure will be over-inflated and real and hidden malicious attacks will be missed. If the value is too small, the small fluctuation variance will cause the threshold to jump and generate false alarm blocking when the network is extremely stable.

[0176] By calling the calibration mapping sequence generated in step six, and using the implementation method of taking the smallest positive integer for the delay step, the calibration mapping sequence items from the fifth sequence position to the end of the sliding time window are taken as statistical objects, which can clarify the effective sample range for the calculation of mean and variance. Summing each calibration mapping sequence item in the statistical object and dividing by the total number of effective samples can yield the arithmetic mean of the calibration mapping sequence. Further calculating the difference between each effective calibration mapping sequence item and the arithmetic mean, squaring each difference and summing them, and then dividing by the total number of effective samples can yield the statistical variance of the calibration mapping sequence. Using the statistical variance as the numerator, adding the arithmetic mean to the positive bias constant of the same scale as the denominator, and dividing the two to obtain the variance normalization quotient, and then adding the arithmetic mean to the variance normalization quotient, can construct a dynamic threshold structure that varies with the dispersion of the calibration mapping sequence. This dynamic threshold structure is generated by the statistical results of the mapping sequence in the current window, which can make the access control decision boundary adjust data-driven according to the fluctuation level in the window, rather than relying solely on a fixed static boundary.

[0177] Step 8: Calculate the decision index by quoting the latest mapping sequence item with the threshold structure, and then block or allow communication ports based on the index comparison result. Specifically, this includes:

[0178] Call the calibration mapping sequence generated in step six and the dynamic threshold structure generated in step seven to obtain the latest calibration mapping sequence item corresponding to the end of the current sliding time window;

[0179] Using the latest calibration mapping sequence as the numerator and the dynamic threshold structure as the denominator, the numerator is divided by the denominator to calculate the dimensionless control risk decision index, which is used to characterize the degree of geometric distortion of the current communication behavior phase trajectory relative to the dynamic threshold structure.

[0180] like Figure 4 As shown, the control risk decision index is compared with the value 1. If the control risk decision index is greater than the value 1, it is determined that the phase trajectory geometric distortion of the current communication behavior exceeds the tolerance boundary, and the session connection of the corresponding port of the target industrial control computer is cut off, and the blocking action is executed.

[0181] If the control risk decision index is less than or equal to 1, it is determined that the current communication is within the nonlinear safety margin, and a release operation is performed on the communication port corresponding to the target industrial control computer.

[0182] By obtaining the latest calibration mapping sequence item corresponding to the end of the current sliding time window and using it as the numerator, and the dynamic threshold structure generated in step seven as the denominator, a quotient operation is performed to obtain a dimensionless control risk decision index. This control risk decision index places the local calibration mapping result at the end of the current window and the dynamic threshold structure formed by the current window statistics under the same judgment scale, so that the complex phase space trajectory changes, cross-node calibration relationships and window statistical boundaries are uniformly converted into directly comparable access control indicators. By comparing the control risk decision index with the value one, clear blocking or allowing conditions can be formed. When the control risk decision index is greater than the value one, the blocking action of the session connection of the corresponding port of the target industrial control computer is executed; when the control risk decision index is less than or equal to the value one, the allowing operation is performed on the corresponding communication port. Thus, step eight transforms the results of the preceding data acquisition, feature construction, calibration mapping and dynamic threshold calculation into clear network access control output.

[0183] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0184] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for network access control of industrial control computers in a power monitoring system, characterized in that, include: Obtain the arrival time and length of packets at the target node, and perform zero-order alignment maintenance on the reference nodes according to the time series to obtain the packet length sequence. Calculate the flow intensity based on the packet length and the time difference of arrival, use the sum of the time difference and the zero constant as the denominator, and then calculate the baseline flow rate. Historical flow intensity is extracted backward according to the set delay step, and the three-dimensional phase space state vectors of the two nodes are spliced ​​and reconstructed. The geometric centroid is obtained by taking the mean value of the phase space state vector, and the topological offset is obtained by calculating the Euclidean distance between the centroids of the two nodes. The adaptive calibration factor is generated by dividing the topology offset by the compensated baseline flow rate and mapping it using a negative exponential function. Perform backward difference operation on the phase space vector, extract the ratio of the outer product to the first-order difference cube of the addition constant, and then multiply by the calibration factor to obtain the mapping sequence; Calculate the mean and variance of the mapped sequence, and then add the mean and variance to the quotient of the mean plus a positive bias constant to construct a dynamic threshold structure. The decision index is obtained by quoting the latest mapping sequence item with the threshold structure, and the communication port is blocked or allowed based on the index comparison result.

2. The industrial control computer network access control method for power monitoring systems according to claim 1, characterized in that, The process of obtaining the arrival time and length of packets from the target node, and performing zero-order alignment maintenance on the reference nodes according to the time series to obtain the packet length sequence, specifically includes: Deploy a bypass mirroring agent on the port of the core switch, set the industrial control computer to be monitored as the target node, and specify the industrial control computer running in the same network segment as the reference node; Record the network data packets of the target node within the sliding time window, and obtain the arrival timestamp and packet length; Set adjacent arrival timestamps to strictly monotonically increasing; Set the system's minimum time measurement resolution constant; Set the time difference between the arrival timestamp of the current data packet and the previous data packet to be greater than or equal to the system's minimum time measurement resolution constant. On the proxy side, zero-order alignment is performed on the data of the reference node. The arrival timestamp of the target node is used as the trigger time. The latest packet length of the reference node before the above time is taken as the alignment packet length to obtain the reference node packet length sequence with equal length alignment.

3. The industrial control computer network access control method for power monitoring systems according to claim 1, characterized in that, The calculation of flow intensity based on packet length and arrival time difference, using the sum of time difference and zero-resistance constant as the denominator, and then determining the baseline flow rate, specifically includes: The traffic intensity of the target node and the reference node is calculated based on the packet length and arrival time difference of adjacent data packets, as follows: The target node's packet length is used as the numerator, and the difference between the arrival timestamp of the current packet and the previous packet is added to the zero-prevention constant and used as the denominator. The two are then divided to obtain the zero-prevention regularized traffic intensity of the target node. Using the same calculation logic, the zero-prevention regularization traffic intensity of the reference node is calculated based on the aligned reference node packet length, the target node time difference, and the aforementioned zero-prevention constant. The traffic intensity of the target node and the reference node within the time window is summed and divided by twice the total number of data packets minus one. The joint arithmetic mean is then calculated to obtain the joint network baseline traffic.

4. The industrial control computer network access control method for power monitoring systems according to claim 1, characterized in that, The step of extracting historical flow intensity backward according to a set delay step and concatenating and reconstructing the three-dimensional phase space state vectors of the two nodes specifically includes: Set a fixed delay step constant, the value of which is the smallest positive integer; Historical flow intensity is extracted backward from the time series using the aforementioned delay step constant; The current flow intensity, the flow intensity delayed by one step, and the flow intensity delayed by two steps are concatenated to reconstruct the three-dimensional phase space state vector of the target node. Using the same reconstruction logic and parallel calculation steps, the three-dimensional phase space state vector of the reference node is reconstructed.

5. A method for network access control of industrial control computers for power monitoring systems according to claim 1, characterized in that, The geometric centroid of the mean value of the phase space state vector is calculated, and the topological offset is obtained by calculating the Euclidean distance between the centroids of the two nodes. Specifically, this includes: The three-dimensional phase space state vectors of all target nodes within the effective time window are summed, divided by the total number of effective vectors, and the average value is calculated to obtain the three-dimensional phase space geometric centroid vector of the target node. The same calculation steps are used to obtain the three-dimensional phase space geometric centroid vector of the reference node; Extract the three corresponding dimension components of the centroid vectors of the target node and the reference node in three-dimensional coordinates; Calculate the difference between the three corresponding dimensional components, square the difference between each dimension and sum them, then take the square root of the sum to calculate the topological offset across nodes.

6. A method for network access control of an industrial control computer for a power monitoring system according to claim 1, characterized in that, The step of dividing the topology offset by the compensated baseline flow rate and generating an adaptive calibration factor using a negative exponential function specifically includes: Set a minimal positive real constant to prevent division by zero errors; Add the calculated network joint baseline traffic to the above minimal positive real constant to obtain the joint baseline traffic after zero-compensation. The calculated cross-node topology offset is used as the numerator, and the joint baseline flow after zero compensation is used as the denominator. The two are divided to obtain the ratio result. The above ratio results are negativeized and used as the independent variable of the natural exponential function for mapping calculation, thus constructing a dimensionless adaptive calibration factor.

7. A method for network access control of industrial control computers for power monitoring systems according to claim 1, characterized in that, The step of performing backward difference operations on the phase space vectors, extracting the ratio of the outer product to the first-order difference cube of the addition constant, and then multiplying by the calibration factor to obtain the mapping sequence, specifically includes: Calculate the backward first-order difference vector and the backward second-order difference vector of the three-dimensional phase space vector of the target node; Perform a three-dimensional vector outer product operation on the first-order difference vector and the second-order difference vector, and take the second norm of the outer product result as the numerator; Find the L2 norm of the first-order difference vector, cube the norm, and add the set zero-prevention constant as the denominator. Dividing the numerator by the denominator yields a geometric complexity index that reflects the phase trajectory change. Cross-node correction is performed by multiplying the calculated adaptive calibration factor with the geometric complexity index to obtain the calibration mapping sequence.

8. A method for network access control of industrial control computers for power monitoring systems according to claim 1, characterized in that, The calculation of the mean and variance of the mapped sequence involves adding the mean and variance to the quotient of the mean plus a positive bias constant to construct a dynamic threshold structure, specifically including: The arithmetic mean of the calibration mapping sequence is obtained by summing the terms of each term in the valid range within the time window and dividing by the total number of valid samples. Calculate the sum of squares of the differences between each term in the sequence and the arithmetic mean, divide by the total number of valid samples, and obtain the statistical variance of the calibration mapping sequence; Set a positive bias constant to prevent the denominator from being zero; Using the above statistical variance as the numerator, add the arithmetic mean to the positive bias constant as the denominator, and divide the two to obtain the quotient. By adding the arithmetic mean to the quotient, a dynamic threshold structure that increases with the variance is constructed.

9. A method for network access control of an industrial control computer for a power monitoring system according to claim 1, characterized in that, The step of quoting the latest mapping sequence item with the threshold structure to obtain a decision index, and then executing communication port blocking or allowing based on the index comparison result, specifically includes: Obtain the latest calibration mapping sequence item corresponding to the end of the current sliding time window; Using the latest calibration mapping sequence as the numerator and the dynamic threshold structure constructed above as the denominator, the two are divided to calculate the dimensionless control risk decision index. Compare the risk control decision index with the value of one. If the control risk decision index is greater than one, then the session connection of the corresponding port of the target industrial control computer will be cut off and a blocking action will be performed. If the control risk decision index is less than or equal to one, then the corresponding communication port will be allowed to pass.

Images