Industrial network security situation prediction method and system based on generative large model

Abstract

The application provides an industrial network security situation prediction method and system based on a generative large model, relates to the technical field of industrial network security, and comprises the following steps: acquiring time series network behavior data of a plurality of monitoring nodes in an industrial network, and extracting multi-dimensional features to obtain a security feature vector; the feature vector is decomposed into a periodic baseline component and a transient disturbance component through frequency domain transformation, and a decomposition situation feature is obtained through sparse constraint screening; semantic space mapping is performed by using a generative large model to obtain a semantic enhanced situation representation, the correlation between security anomalies among nodes is calculated based on the decomposition situation feature, and a dynamic correlation matrix is constructed; a propagation operator is constructed based on the dynamic correlation matrix, the semantic enhanced situation representation is subjected to multi-step iterative propagation, a decay factor is introduced to simulate abnormal influence diffusion, and a security situation prediction result in a future time window is obtained.

Description

Technical Field

[0001] This invention relates to the field of industrial network security technology, and in particular to a method and system for predicting industrial network security situation based on generative large models. Background Technology

[0002] As a critical infrastructure supporting industrial production and operation, the security situation prediction of industrial networks is of great significance for ensuring the stable operation of industrial systems. With the rapid development of the Industrial Internet, the industrial network environment is becoming increasingly complex, and the cybersecurity threats it faces are characterized by diversification, concealment, and intelligence. Traditional cybersecurity protection methods mostly adopt a passive defense model, which is insufficient to effectively cope with rapidly evolving security threats. Therefore, conducting research on industrial cybersecurity situation prediction, and realizing the transformation from passive defense to proactive early warning, has become an important research direction in the field of industrial cybersecurity.

[0003] Existing technologies primarily employ single-dimensional feature analysis methods, failing to fully consider the fundamental differences between periodic patterns and sudden anomalies inherent in industrial network security data. Industrial network environments exhibit both normal behavioral fluctuations caused by production cycles and abnormal disturbances triggered by security threats, which are often overlapping and difficult to distinguish in the time domain. Traditional methods lack effective decomposition mechanisms, leading to normal fluctuations being misjudged as security threats, or genuine abnormal signals being masked by periodic noise, severely impacting the accuracy of situational prediction. Existing methods also lack the ability to understand the deep semantic information of industrial network security data. Industrial network security situational awareness is not only reflected in numerical statistical characteristics but also contains complex semantic relationships and contextual dependencies. Traditional prediction models rely mainly on numerical calculations, failing to capture the logical connections, causal relationships, and potential attack intentions between security events, resulting in uninterpretable prediction results and hindering effective support for security decision-making. Summary of the Invention

[0004] This invention provides a method and system for predicting industrial network security situation based on generative large models, which can solve the problems in the prior art.

[0005] A first aspect of this invention provides a method for predicting industrial network security situation based on a generative large model, comprising:

[0006] The network behavior data of multiple monitoring nodes in the industrial network are acquired in a time series, and multi-dimensional feature extraction is performed on the network behavior data to obtain the industrial network security feature vector.

[0007] The industrial network security feature vector is decomposed into periodic baseline components and transient disturbance components through frequency domain transformation, and the transient disturbance components are filtered through sparsity constraints to obtain the decomposed situation features.

[0008] A generative large model is used to perform semantic space mapping between the periodic baseline components and the decomposed situation features to obtain a semantically enhanced situation representation. Based on the decomposed situation features, the correlation of security anomalies among multiple monitoring nodes is calculated, and a dynamic correlation matrix among industrial network nodes is constructed.

[0009] Based on the dynamic correlation matrix, a propagation operator for industrial network security situation is constructed. The semantically enhanced situation representation is propagated through multiple iterations using the propagation operator. In each iteration, the periodic baseline component is used as a propagation stability constraint, and the decomposed situation features are used as a propagation excitation source. In the iteration process, an attenuation factor is introduced to simulate the diffusion and attenuation of abnormal influences over time and space, thereby obtaining the prediction result of industrial network security situation within the future time window.

[0010] Multi-dimensional feature extraction is performed on the network behavior data to obtain the industrial network security feature vector, which includes:

[0011] From the network behavior data, the interaction patterns of industrial network communication protocols in the protocol layer are extracted to obtain protocol layer features, and the statistical periods of network data transmission in the traffic layer are extracted to obtain traffic layer features.

[0012] The protocol layer features are sequence encoded, and the protocol interaction behaviors at different times are arranged in chronological order to form a protocol interaction sequence. The protocol type identifier and interaction state identifier are extracted from each protocol interaction behavior in the protocol interaction sequence and mapped to a vector space to obtain a protocol behavior vector. The protocol behavior vectors are combined in chronological order to form a time-series vector.

[0013] The flow layer features are distributed and modeled. The flow statistics within a preset time window are statistically analyzed. The skewness and kurtosis of the flow statistics are calculated to obtain statistical moments. Based on the statistical moments, the distribution shape parameters of the flow statistics are determined.

[0014] The time-series vector and the distribution pattern parameter are concatenated to form the industrial network security feature vector that integrates the time-series interaction mode and the traffic distribution pattern.

[0015] The industrial network security feature vector is decomposed into periodic baseline components and transient disturbance components through frequency domain transformation. The transient disturbance components are then filtered through sparsity constraints to obtain the decomposed situational features, including:

[0016] The industrial network security feature vector is mapped from the time domain to the frequency domain to obtain frequency domain features. The spectral coefficients in the frequency domain features are arranged in ascending order of frequency and subjected to exponential fitting to obtain an energy decay curve. The fitting residual between the actual energy value at each frequency position and the fitted value of the energy decay curve is calculated. The fitting residual is subjected to second-order difference operation to obtain residual curvature. The frequency position where the sign of the residual curvature changes is identified as the frequency band boundary point. The spectral coefficients before and after the frequency band boundary point are divided into dominant frequency bands and non-dominant frequency bands, respectively.

[0017] The spectral coefficients of the dominant frequency band and the non-dominant frequency band are respectively subjected to inverse frequency domain transformation to obtain the periodic baseline component and the initial transient disturbance component;

[0018] The initial transient disturbance component is decomposed into multiple time scales. The amplitude sequence of the disturbance component at each time scale is calculated. The sparsity is obtained by calculating the ratio of the zero norm of the amplitude sequence to the sequence length. The time scale with the largest sparsity is selected as the dominant scale. The amplitude of the time corresponding to the non-zero position of the amplitude at the dominant scale is extracted in the initial transient disturbance component. After sorting in descending order, the disturbance signal corresponding to the time corresponding to the preset proportion is retained to obtain the filtered transient disturbance component. This component is combined with the periodic baseline component to form the decomposed situation feature.

[0019] Based on the decomposed situation characteristics, the correlation of security anomalies among multiple monitoring nodes is calculated, and a dynamic correlation matrix among industrial network nodes is constructed, including:

[0020] Extract the filtered transient disturbance components from any two monitoring nodes and calculate the time-series cross-response characteristics in the time dimension. Identify the peak amplitude and time offset corresponding to the peak position of the time-series cross-response characteristics. Use the time offset as the anomaly propagation delay. Perform a nonlinear mapping between the peak amplitude and the reciprocal of the anomaly propagation delay to obtain the time-series coupling strength.

[0021] The energy distribution of the periodic baseline components of any two monitoring nodes in the frequency domain is calculated respectively. By constructing the frequency domain coherence metric matrix between the two energy distributions, the frequency domain energy coupling spectrum is obtained. The frequency domain energy coupling spectrum is adaptively weighted and integrated over the entire frequency band and mapped to the zero-to-one interval through hyperbolic tangent transformation to obtain the frequency domain coupling strength.

[0022] The temporal coupling strength and the frequency domain coupling strength are adaptively fused to obtain the safety anomaly correlation between the two monitoring nodes; all node pairs in the monitoring nodes are traversed, and multiple monitoring nodes are used as row and column indices. The safety anomaly correlation between each node pair is used as the matrix element at the corresponding row and column position to construct the dynamic correlation matrix between industrial network nodes.

[0023] The propagation operators for constructing industrial network security posture based on the aforementioned dynamic correlation matrix include:

[0024] Calculate the difference metric between the dynamic correlation matrix and the corresponding transpose matrix. When the difference metric exceeds a preset symmetry threshold, the dynamic correlation matrix and the transpose matrix are averaged element-wise to obtain a symmetric correlation matrix. Otherwise, the dynamic correlation matrix is ​​used as the symmetric correlation matrix.

[0025] By applying an iterative power transformation to the symmetric correlation matrix, the scaling factors of the symmetric correlation matrix in different directions are identified as a set of eigenvalues. For each eigenvalue in the set of eigenvalues, the vector in which the symmetric correlation matrix remains unchanged in the scaling direction corresponding to the eigenvalue is determined as an eigenvector.

[0026] The eigenvalue set is sorted in descending order of eigenvalues. The cumulative energy contribution rate sequence of the eigenvalue set is calculated. The eigenvalues ​​that reach a preset energy percentage and their corresponding eigenvectors in the cumulative energy contribution rate sequence are selected. The corresponding eigenvectors are arranged in columns and multiplied with their own transpose to construct a spectral projection operator. The spectral projection operator is used to perform a linear transformation on the symmetric correlation matrix. The linearly transformed symmetric correlation matrix is ​​used as the propagation operator for industrial network security situation.

[0027] The semantically enhanced situation representation is propagated through a multi-step iterative process using the propagation operator. In each iteration, the periodic baseline component is used as a propagation stability constraint, and the decomposed situation features are used as the propagation excitation source.

[0028] The distribution of the semantically enhanced situational representation across multiple monitoring nodes is used as an iterative state vector.

[0029] The amplitude of the periodic baseline component corresponding to each monitoring node at the current moment is extracted and mapped to a preset constraint interval through nonlinear saturation transformation to obtain the stability constraint coefficient of the monitoring node. The stability constraint coefficients of all monitoring nodes are arranged into a diagonal matrix to form a propagation stability constraint. The energy intensity of the transient disturbance component in the decomposed situation feature corresponding to each monitoring node at the current moment is extracted. The energy intensity is integrated in the time domain to obtain the cumulative excitation intensity of the monitoring node. The cumulative excitation intensities of all monitoring nodes are arranged into a column vector to form a propagation excitation source.

[0030] The propagated state vector is obtained by convolving the propagation operator with the iterative state vector. An adaptive threshold pruning transformation based on the propagation stability constraint is applied to the propagated state vector to perform nonlinear compression mapping on the state components outside the dynamic boundary range of the propagation stability constraint. The vector is then combined with the propagation excitation source to obtain the excitation propagation state.

[0031] The iteration process introduces a decay factor to simulate the diffusion and decay of anomalous effects over time and space, including:

[0032] Set the initial value of the iteration step counter to zero, and record the node that first shows an abnormal response as the set of abnormal source nodes;

[0033] The current value of the iteration step counter is used as the exponent to calculate the time decay component. The topological distance between each monitoring node and the set of abnormal source nodes is calculated and a Gaussian radial basis transformation is applied to obtain the spatial decay component. The time decay component and the spatial decay component are subjected to a Hadamard product to form an adaptive decay weight distribution. The local variance of the propagation state after excitation is calculated, and the monitoring nodes whose local variance exceeds a preset fluctuation threshold are identified as high fluctuation region nodes. The decay weights corresponding to the high fluctuation region nodes are attenuated and compensated. The adaptive decay weight distribution after attenuation compensation is fused with the propagation state after excitation node by node to obtain the propagation state after decay. The iteration step counter is then incremented by one.

[0034] Calculate the relative rate of change between the decayed propagation state and the decayed propagation state of the previous iteration step. If the relative rate of change is lower than a preset convergence threshold or the iteration step counter reaches a preset maximum iteration step, the decayed propagation state is used as the industrial network security situation prediction result.

[0035] A second aspect of this invention provides an industrial network security situation prediction system based on a generative large model, comprising:

[0036] The first unit is used to acquire network behavior data of multiple monitoring nodes in the industrial network in a time series, and to extract multi-dimensional features from the network behavior data to obtain an industrial network security feature vector.

[0037] The second unit is used to decompose the industrial network security feature vector into periodic baseline components and transient disturbance components through frequency domain transformation, and to filter the transient disturbance components through sparsity constraints to obtain the decomposed situation features.

[0038] The third unit is used to perform semantic space mapping between the periodic baseline components and the decomposed situation features using a generative large model to obtain a semantically enhanced situation representation. Based on the decomposed situation features, it calculates the security anomaly correlation between multiple monitoring nodes and constructs a dynamic correlation matrix between industrial network nodes.

[0039] The fourth unit is used to construct a propagation operator for industrial network security situation based on the dynamic correlation matrix. The propagation operator is used to propagate the semantically enhanced situation representation in multiple steps. In each step of the iteration, the periodic baseline component is used as a propagation stability constraint, the decomposed situation features are used as a propagation excitation source, and an attenuation factor is introduced in the iteration process to simulate the diffusion and attenuation of the abnormal influence over time and space, so as to obtain the prediction result of industrial network security situation within the future time window.

[0040] A third aspect of the embodiments of the present invention,

[0041] An electronic device is provided, comprising:

[0042] processor;

[0043] Memory used to store processor-executable instructions;

[0044] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0045] Fourth aspect of the present invention,

[0046] A computer-readable storage medium is provided, having stored thereon computer program instructions, which, when executed by a processor, implement the aforementioned method.

[0047] The beneficial effects of this application are as follows:

[0048] The industrial network security situation prediction method based on generative large models provided by this invention decomposes industrial network behavior data into periodic baseline components and transient disturbance components through frequency domain transformation, and uses sparsity constraints to screen key disturbance features. This method can effectively separate normal business patterns from abnormal behavior patterns, improve the accuracy of identifying industrial network security threats, avoid the problem of confusing normal fluctuations with real threats in traditional methods, and enhance the accuracy of situational awareness.

[0049] This invention uses a generative large model to perform semantic space mapping on the decomposed situation features, transforming numerical features into a representation with rich semantic information. This makes the network security situation more interpretable and expressive. At the same time, it constructs a dynamic correlation matrix based on the correlation of security anomalies, which can capture the correlation between different monitoring nodes in the industrial network and the threat propagation path. This realizes the transformation from single-point defense to global collaborative protection, and improves the integrity and systematicness of industrial network security protection.

[0050] This invention constructs a propagation operator and combines it with periodic baseline components as stability constraints and decomposed situational features as propagation excitation sources for multi-step iterative propagation. It also introduces an attenuation factor to simulate the spatiotemporal diffusion characteristics of anomaly effects. This allows for accurate prediction of the evolution trend of the security situation within future time windows, realizing a shift from passive response to proactive prediction in security protection. This provides valuable response time for industrial network security protection and significantly improves the forward-looking and proactive defense capabilities of industrial control systems against network security threats. Attached Figure Description

[0051] Figure 1 This is a flowchart illustrating the industrial network security situation prediction method based on a generative large model, as described in an embodiment of the present invention.

[0052] Figure 2 This is a schematic diagram of the frequency domain decomposition process for industrial network security feature vectors. Detailed Implementation

[0053] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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.

[0054] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0055] Figure 1 This is a flowchart illustrating the industrial network security situation prediction method based on a generative large model, as described in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0056] The network behavior data of multiple monitoring nodes in the industrial network are acquired in a time series, and multi-dimensional feature extraction is performed on the network behavior data to obtain the industrial network security feature vector.

[0057] The industrial network security feature vector is decomposed into periodic baseline components and transient disturbance components through frequency domain transformation, and the transient disturbance components are filtered through sparsity constraints to obtain the decomposed situation features.

[0058] A generative large model is used to perform semantic space mapping between the periodic baseline components and the decomposed situation features to obtain a semantically enhanced situation representation. Based on the decomposed situation features, the correlation of security anomalies among multiple monitoring nodes is calculated, and a dynamic correlation matrix among industrial network nodes is constructed.

[0059] Based on the dynamic correlation matrix, a propagation operator for industrial network security situation is constructed. The semantically enhanced situation representation is propagated through multiple iterations using the propagation operator. In each iteration, the periodic baseline component is used as a propagation stability constraint, and the decomposed situation features are used as a propagation excitation source. In the iteration process, an attenuation factor is introduced to simulate the diffusion and attenuation of abnormal influences over time and space, thereby obtaining the prediction result of industrial network security situation within the future time window.

[0060] In one optional implementation, multi-dimensional feature extraction is performed on the network behavior data to obtain an industrial network security feature vector, including:

[0061] From the network behavior data, the interaction patterns of industrial network communication protocols in the protocol layer are extracted to obtain protocol layer features, and the statistical periods of network data transmission in the traffic layer are extracted to obtain traffic layer features.

[0062] The protocol layer features are sequence encoded, and the protocol interaction behaviors at different times are arranged in chronological order to form a protocol interaction sequence. The protocol type identifier and interaction state identifier are extracted from each protocol interaction behavior in the protocol interaction sequence and mapped to a vector space to obtain a protocol behavior vector. The protocol behavior vectors are combined in chronological order to form a time-series vector.

[0063] The flow layer features are distributed and modeled. The flow statistics within a preset time window are statistically analyzed. The skewness and kurtosis of the flow statistics are calculated to obtain statistical moments. Based on the statistical moments, the distribution shape parameters of the flow statistics are determined.

[0064] The time-series vector and the distribution pattern parameter are concatenated to form the industrial network security feature vector that integrates the time-series interaction mode and the traffic distribution pattern.

[0065] In the acquired network behavior data, at the protocol level, the specific type of industrial network communication protocol is identified from the raw network data packets. This includes commonly used industrial control system protocols such as Modbus, OPC, and DNP3. By parsing the packet header information and payload content, key fields such as function code, register address, and data object are extracted from each protocol data unit. For the Modbus protocol, the function code is recorded as either reading coil status or writing to a holding register, along with the start address and data length. For the OPC protocol, the subscription item, data item identifier, and read / write operation type are extracted. These extracted protocol fields collectively constitute the interaction mode characteristics of the protocol layer, reflecting the communication intent and control logic between devices.

[0066] At the traffic level, network data transmission is statistically analyzed at fixed time intervals, with a time window of 10 seconds. Within each time window, basic indicators such as the number of data packets, total bytes, average packet length, maximum packet length, and minimum packet length are recorded. For actual monitoring data from an industrial production line, 320 data packets were captured within a 10-second window, totaling 84,000 bytes, with an average packet length of 262 bytes, a maximum packet length of 1,500 bytes, and a minimum packet length of 64 bytes. These statistical indicators constitute the basic statistical periodic characteristics of the traffic layer, reflecting the network load status and data transmission patterns.

[0067] Sequence encoding is performed on the extracted protocol layer features, arranging the monitored protocol interactions in chronological order to form a continuous protocol interaction sequence. Within a specific monitoring period, interactions such as the master station sending a read request to the slave station, the slave station returning a data response, the master station sending a write command, and the slave station returning an acknowledgment are captured sequentially. Each protocol interaction in the sequence is analyzed in fine granularity, extracting a protocol type identifier and an interaction status identifier. The protocol type identifier records which industrial protocol the interaction belongs to (e.g., Modbus protocol is labeled as type 1, OPC protocol as type 2). The interaction status identifier records the specific operation of the interaction (e.g., read operation is labeled as status 01, write operation as status 02, successful response as status 10, and failed response as status 11).

[0068] A 64-dimensional vector space is established, mapping protocol type identifiers to the first 32 dimensions and interaction state identifiers to the last 32 dimensions. For Modbus read requests, the protocol type is assigned a value of 1.0 in the 5th dimension of the first 32 dimensions, with the remaining positions being 0; the read operation state is assigned a value of 1.0 in the 2nd dimension of the last 32 dimensions, with the remaining positions being 0. This forms a 64-dimensional protocol behavior vector, fully describing the characteristics of this protocol interaction. Successive protocol behavior vectors are concatenated in chronological order. If 20 protocol interaction behaviors are captured within a monitoring period, the 20 64-dimensional protocol behavior vectors are connected in chronological order to form a 1280-dimensional time-series vector, which fully preserves the temporal dependencies and evolution patterns of the protocol interactions.

[0069] Distribution modeling was performed on the traffic layer characteristics, and in-depth analysis was conducted on the obtained traffic statistics within a preset 10-second time window. The skewness index of the traffic statistics was calculated to measure the symmetry of the data distribution relative to the mean. Specifically, the difference between the length of each data packet and the average packet length was cubed, the sum of the cubed differences of all data packets was divided by the total number of data packets, and then divided by the cube of the standard deviation to obtain the skewness value. Within a certain time window, a skewness value of 0.75 was obtained, indicating a right-skewed distribution of packet lengths and the presence of a small number of oversized data packets. The kurtosis index of the traffic statistics was calculated to measure the sharpness of the data distribution at its peaks. The difference between the length of each data packet and the average packet length was raised to the fourth power, the sum of the fourth power differences of all data packets was divided by the total number of data packets, then divided by the fourth power of the standard deviation, and finally subtracted by 3 to obtain the kurtosis value. Within this time window, a kurtosis value of 2.13 was obtained, indicating that the packet length distribution is more angular than a normal distribution, and the data concentration is relatively high.

[0070] The calculated skewness and kurtosis are used as statistical moments to characterize the distribution pattern of the flow statistics. Based on these two statistical moments, specific morphological parameters of the flow distribution are determined: skewness of 0.75 characterizes the degree of skewness, and kurtosis of 2.13 characterizes the degree of steepness. The scale parameter of the distribution is further calculated, and the standard deviation of 285.6 characterizes the dispersion of the data. Simultaneously, the location parameter of the distribution is determined, and the average packet length of 262.0 characterizes the central location of the data. These four parameters together constitute a complete set of parameters describing the flow distribution pattern, organized into a four-dimensional vector with values ​​of 0.75, 2.13, 285.6, and 262.0, respectively.

[0071] A vector concatenation operation is performed to fuse the temporal vector and the distribution morphological parameter vector. The 1280-dimensional temporal vector and the 4-dimensional distribution morphological parameter vector are connected along the feature dimension to form a comprehensive feature vector of dimension 1284. The first 1280 dimensions of this vector carry the temporal pattern information of protocol interactions, reflecting the communication behavior sequence and operational logic between devices; the last 4 dimensions carry the morphological information of traffic distribution, reflecting the statistical characteristics and load patterns of network transmission. This fusion method preserves both the fine-grained interaction details at the protocol level and covers the macroscopic statistical laws at the traffic level, forming a complete industrial network security feature vector. This feature vector can simultaneously characterize the temporal evolution characteristics and traffic distribution characteristics of industrial network behavior, providing a comprehensive feature representation foundation for subsequent security threat detection and abnormal behavior identification.

[0072] In one optional implementation, the industrial network security feature vector is decomposed into periodic baseline components and transient disturbance components through frequency domain transformation, and the transient disturbance components are filtered through sparsity constraints to obtain the decomposed situational features, including:

[0073] The industrial network security feature vector is mapped from the time domain to the frequency domain to obtain frequency domain features. The spectral coefficients in the frequency domain features are arranged in ascending order of frequency and subjected to exponential fitting to obtain an energy decay curve. The fitting residual between the actual energy value at each frequency position and the fitted value of the energy decay curve is calculated. The fitting residual is subjected to second-order difference operation to obtain residual curvature. The frequency position where the sign of the residual curvature changes is identified as the frequency band boundary point. The spectral coefficients before and after the frequency band boundary point are divided into dominant frequency bands and non-dominant frequency bands, respectively.

[0074] The spectral coefficients of the dominant frequency band and the non-dominant frequency band are respectively subjected to inverse frequency domain transformation to obtain the periodic baseline component and the initial transient disturbance component;

[0075] The initial transient disturbance component is decomposed into multiple time scales. The amplitude sequence of the disturbance component at each time scale is calculated. The sparsity is obtained by calculating the ratio of the zero norm of the amplitude sequence to the sequence length. The time scale with the largest sparsity is selected as the dominant scale. The amplitude of the time corresponding to the non-zero position of the amplitude at the dominant scale is extracted in the initial transient disturbance component. After sorting in descending order, the disturbance signal corresponding to the time corresponding to the preset proportion is retained to obtain the filtered transient disturbance component. This component is combined with the periodic baseline component to form the decomposed situation feature.

[0076] like Figure 2 As shown, the method includes:

[0077] After receiving real-time network traffic data, the industrial network security monitoring system performs frequency domain transformation on the industrial network security feature vector composed of this data. The Fast Fourier Transform (FFT) algorithm is used to map the feature vector from the time domain to the frequency domain, generating a corresponding frequency domain feature representation. Assuming the original feature vector contains 1024 sampling points in the time domain, the transformation yields 512 complex spectral coefficients, each containing amplitude and phase information.

[0078] The obtained spectral coefficients are preprocessed. The energy value at each frequency position is obtained by calculating the square of the modulus of the complex spectral coefficients. The energy values ​​corresponding to these 512 frequency positions are arranged in ascending order of frequency from low to high to form an energy distribution sequence. An exponential function fitting operation is performed on this energy distribution sequence. The least squares method is used to determine the coefficient parameters of the exponential function during the fitting process. It is assumed that the fitted energy decay curve represents an exponential decrease in the baseline value with increasing frequency. The fitted value is 8500 units at the 10 Hz position, 320 units at the 100 Hz position, and 45 units at the 200 Hz position.

[0079] Calculate the difference between the actual energy value and the fitted value of the energy decay curve at each frequency position. This difference is defined as the fitting residual. At a frequency of 10 Hz, the actual energy value is 8620 units, the fitted value is 8500 units, and the fitting residual is 120 units. At a frequency of 55 Hz, the actual energy value is 1850 units, the fitted value is 1200 units, and the fitting residual is 650 units. Construct a residual sequence with a length of 512 for all frequency positions. Perform a second-order difference operation on the residual sequence. Specifically, calculate the rate of change of the rate of change of three adjacent residual values. Assuming the residual value is r_k at the k-th sampling point, r_k1 at the (k+1)-th position, and r_k2 at the (k+2)-th position, the first-order difference is r_k1-r_k and r_k2-r_k1. The second-order difference is the sum of the previous and subsequent first-order difference values, and this value reflects the curvature characteristics of the residual curve.

[0080] The residual curvature values ​​are iterated across all frequency positions, and the locations where the curvature sign changes are detected. Before the 78 Hz position, the residual curvature remains positive, indicating that the residual curve bulges upwards. At the 78 Hz position, the residual curvature changes from positive to negative, indicating that the direction of curvature has reversed. This frequency position is identified as the frequency band boundary. Based on this boundary point, the entire spectral coefficient sequence is divided into two frequency band regions: the interval from 0 Hz to 78 Hz is defined as the dominant band, which contains 79 spectral coefficients; the interval from 78 Hz to 256 Hz is defined as the non-dominant band, which contains 433 spectral coefficients.

[0081] An inverse Fourier transform was performed on the 79 spectral coefficients of the dominant frequency band, converting the frequency domain representation back to the time domain representation, yielding a periodic baseline component. This component exhibits smooth periodic fluctuations in the time domain, with amplitude variations ranging from -15 to +18 and a period length of approximately 26 sampling points. This periodic baseline component reflects the normal operation mode of industrial network traffic. An inverse Fourier transform was performed on the 433 spectral coefficients of the non-dominant frequency band, yielding an initial transient disturbance component. This component exhibits irregular abrupt changes in the time domain, with its amplitude close to zero most of the time, and significant spike fluctuations only appearing at a few moments.

[0082] The initial transient disturbance component is decomposed into sub-components at multiple time scales using wavelet packet decomposition. The first scale corresponds to the finest time granularity, with a time window width of 2 sampling points; the eighth scale corresponds to the coarsest time granularity, with a time window width of 256 sampling points. At each time scale, the amplitude sequence of the disturbance component is extracted. For the third scale, the time window width is 8 sampling points, dividing the 1024 sampling points into 128 time windows. Within each window, the maximum absolute value of the amplitude is calculated, forming an amplitude sequence of length 128.

[0083] For the amplitude sequence at each time scale, a sparsity index is calculated by counting the number of non-zero elements in the amplitude sequence. An amplitude threshold of 0.01 is set, and elements with absolute values ​​less than this threshold are considered zero. In the amplitude sequence at the third scale, 19 out of 128 elements have amplitudes exceeding the threshold and are defined as non-zero elements. The zero norm is defined as the number of non-zero elements; at this scale, the zero norm is 19. Sparsity is calculated by dividing the zero norm by the sequence length; the sparsity at the third scale is 19 / 128 = 0.148. The sparsity is calculated sequentially for all eight scales, resulting in a sparsity sequence of 0.652, 0.485, 0.148, 0.095, 0.063, 0.042, 0.031, and 0.028.

[0084] The time scale with the highest sparsity is identified as the dominant scale. In this case, the sparsity of the first scale is 0.652, which is the highest among all scales. Therefore, the first scale is selected as the dominant scale. The amplitude distribution is re-analyzed under this dominant scale to locate the time points where the amplitude is non-zero. The first scale divides the 1024 sampling points into 512 time windows, of which the amplitude of 337 windows exceeds the threshold. The center time of these 337 time windows is recorded, and these times are mapped to the initial transient perturbation components to extract the amplitude of the perturbation signal at the corresponding time.

[0085] The 337 extracted amplitude values ​​were sorted in descending order of absolute value. With a retention rate of 30%, the number of time points to be retained was calculated to be 337 × 0.3 = 101. The top 101 time points with the largest amplitudes were selected, and the corresponding disturbance signals were retained. The amplitudes of the disturbance signals at the remaining time points were set to zero. The filtered transient disturbance components retained non-zero amplitudes at only 101 locations out of 1024 sampling points, while the remaining 923 locations had zero amplitudes, achieving a sparse representation of transient disturbances. The filtered transient disturbance components were concatenated with the periodic baseline components to form a decomposed situational feature containing two parts of information. This feature vector is 2048 dimensions long; the first 1024 dimensions store the periodic baseline components, and the last 1024 dimensions store the filtered transient disturbance components, fully characterizing the normal patterns and abnormal disturbance characteristics of the industrial network security situation.

[0086] In one optional implementation, based on the decomposed situation characteristics, calculating the correlation of security anomalies among multiple monitoring nodes and constructing a dynamic correlation matrix among industrial network nodes includes:

[0087] Extract the filtered transient disturbance components from any two monitoring nodes and calculate the time-series cross-response characteristics in the time dimension. Identify the peak amplitude and time offset corresponding to the peak position of the time-series cross-response characteristics. Use the time offset as the anomaly propagation delay. Perform a nonlinear mapping between the peak amplitude and the reciprocal of the anomaly propagation delay to obtain the time-series coupling strength.

[0088] The energy distribution of the periodic baseline components of any two monitoring nodes in the frequency domain is calculated respectively. By constructing the frequency domain coherence metric matrix between the two energy distributions, the frequency domain energy coupling spectrum is obtained. The frequency domain energy coupling spectrum is adaptively weighted and integrated over the entire frequency band and mapped to the zero-to-one interval through hyperbolic tangent transformation to obtain the frequency domain coupling strength.

[0089] The temporal coupling strength and the frequency domain coupling strength are adaptively fused to obtain the safety anomaly correlation between the two monitoring nodes; all node pairs in the monitoring nodes are traversed, and multiple monitoring nodes are used as row and column indices. The safety anomaly correlation between each node pair is used as the matrix element at the corresponding row and column position to construct the dynamic correlation matrix between industrial network nodes.

[0090] Transient disturbance components from any two nodes are extracted from the data stream of the monitoring nodes. The first node is designated as node A, and the second as node B. The transient disturbance component of node A is a numerical sequence sampled from the time series, containing 256 sampling points at a sampling interval of 10 milliseconds. The transient disturbance component of node B also contains 256 sampling points at the same sampling interval. Using the transient disturbance component of node A as the reference sequence and the transient disturbance component of node B as the comparison sequence, the comparison sequence is slid along the time axis, with each slide involving one sampling point. At each sliding position, the reference sequence and the slid comparison sequence are multiplied point-by-point. All product values ​​are summed to obtain the response value at that sliding position. This process is repeated for all sliding positions, with the sliding range set from -128 to +128 sampling points, resulting in a time-series cross-response characteristic curve containing 257 response values.

[0091] The maximum response value is searched in the time-series cross-response characteristic curve; this maximum response value is the peak amplitude. Assume that in a certain calculation, the peak amplitude occurs at a sliding position of +15 sampling points, with a value of 8.73. The physical time corresponding to this sliding position is 15 × 10^-15 milliseconds, or 150 milliseconds. This time difference is identified as the anomalous propagation delay. The reciprocal of the anomalous propagation delay is calculated; the reciprocal of 150 milliseconds is 6.67, in seconds. Multiplying the peak amplitude 8.73 by the reciprocal 6.67 yields 58.23. This intermediate result is nonlinearly mapped using an exponential function. Using 58.23 as the input parameter of the exponential function, with the base set to the base of the natural logarithm, the calculated exponent value is divided by this exponent value plus 1. The mapping result is close to 1.0, and this mapping result is the temporal coupling strength.

[0092] Extract the periodic baseline components of nodes A and B, both of which are time series with 512 sampling points. Perform a Discrete Fourier Transform (DFT) on the periodic baseline component of node A to convert the time-domain signal into a frequency-domain representation, resulting in 257 frequency components. Each frequency component contains a real part and an imaginary part. Calculate the sum of the squares of the real and imaginary parts of each frequency component to obtain the energy value at that frequency. Iterate through all 257 frequency components to obtain the energy distribution sequence of node A in the frequency domain. Use the same method to perform a DFT on the periodic baseline component of node B and calculate the energy distribution sequence.

[0093] For each frequency point, the complex frequency components of nodes A and B at that frequency point are extracted. Assuming that at frequency index 42, the real part of the complex component of node A is 3.2 and the imaginary part is 1.8, and the real part of the complex component of node B is 2.7 and the imaginary part is 2.1, the conjugate product of the complex components of node A and node B is calculated to obtain the cross-spectral density. The square of the modulus of the cross-spectral density is divided by the product of the energy values ​​of node A and node B to obtain the coherence metric value of 0.82 for that frequency point. Traversing all 257 frequency points yields the frequency domain coherence metric matrix, which is actually a vector containing 257 coherence metric values, forming the frequency domain energy coupling spectrum.

[0094] A weighted integral is applied to the frequency domain energy coupling spectrum. The weighting function is defined to increase linearly with frequency, with a weighting coefficient of 0.5 for low frequencies and 1.5 for high frequencies. The coherence metric at each frequency point is multiplied by its corresponding weighting coefficient. All products are summed and divided by the sum of the weighting coefficients to obtain the weighted integral value. Assuming the calculated weighted integral value is 0.68, this value is input into the hyperbolic tangent function. The hyperbolic tangent function maps the input value to the interval -1 to 1. Since the input value is positive, the output value falls between 0 and 1. A linear transformation is then performed on the output value of the hyperbolic tangent function, mapping the interval from -1 to +1 to the interval 0 to 1. This transformation is achieved by adding 1 to the hyperbolic tangent output value and dividing by 2, resulting in a frequency domain coupling strength of 0.74.

[0095] An adaptive fusion weight parameter is defined, which dynamically adjusts based on the relative magnitudes of temporal and frequency domain coupling strengths. The absolute value of the difference between the temporal and frequency domain coupling strengths is calculated; a larger absolute value indicates a greater divergence between the two strengths, requiring more reliance on the larger strength index during fusion. The absolute value of the difference is mapped to a weight range of 0.3 to 0.7 using an S-curve function, with the center point of the S-curve function set at a difference of 0.2. Assuming a temporal coupling strength of 1.0, a frequency domain coupling strength of 0.74, and an absolute difference of 0.26, the fusion weight for the temporal coupling strength is 0.62, and the fusion weight for the frequency domain coupling strength is 0.38, obtained through S-curve mapping. Multiplying the temporal coupling strength by its fusion weight yields 0.62, and multiplying the frequency domain coupling strength by its fusion weight yields 0.28. The sum of these two weighted results gives a security anomaly correlation of 0.90.

[0096] Iterate through all monitoring nodes. Assuming there are 16 monitoring nodes in the system, calculate the safety anomaly correlation of each node with the other 15 nodes. The safety anomaly correlation between node 1 and node 2 is 0.90, and the safety anomaly correlation between node 1 and node 3 is 0.45. Repeat this process for all node pairs. Create a 16x16 matrix structure, where both row and column indices correspond to the monitoring node numbers. Fill the matrix with the correlation value of 0.90 between node 1 and node 2 in the first row and second column, and simultaneously fill it in the second row and first column to ensure matrix symmetry. The diagonal positions of the matrix correspond to the correlation between the same node and itself, set to 1.0. Complete the calculation and matrix filling for all 240 node pairs to obtain a complete dynamic correlation matrix. This matrix reflects the safety anomaly correlations between nodes in the industrial network in real time.

[0097] In one optional implementation, the propagation operator for constructing the industrial network security posture based on the dynamic correlation matrix includes:

[0098] Calculate the difference metric between the dynamic correlation matrix and the corresponding transpose matrix. When the difference metric exceeds a preset symmetry threshold, the dynamic correlation matrix and the transpose matrix are averaged element-wise to obtain a symmetric correlation matrix. Otherwise, the dynamic correlation matrix is ​​used as the symmetric correlation matrix.

[0099] By applying an iterative power transformation to the symmetric correlation matrix, the scaling factors of the symmetric correlation matrix in different directions are identified as a set of eigenvalues. For each eigenvalue in the set of eigenvalues, the vector in which the symmetric correlation matrix remains unchanged in the scaling direction corresponding to the eigenvalue is determined as an eigenvector.

[0100] The eigenvalue set is sorted in descending order of eigenvalues. The cumulative energy contribution rate sequence of the eigenvalue set is calculated. The eigenvalues ​​that reach a preset energy percentage and their corresponding eigenvectors in the cumulative energy contribution rate sequence are selected. The corresponding eigenvectors are arranged in columns and multiplied with their own transpose to construct a spectral projection operator. The spectral projection operator is used to perform a linear transformation on the symmetric correlation matrix. The linearly transformed symmetric correlation matrix is ​​used as the propagation operator for industrial network security situation.

[0101] After obtaining the dynamic correlation matrix, it needs to undergo symmetry detection and processing. This process involves performing a matrix transpose operation and calculating the difference metric. The dynamic correlation matrix is ​​denoted as matrix A, with dimensions N×N, where N represents the total number of monitoring nodes in the industrial network. Performing a transpose operation on matrix A yields the transpose matrix A_T, where the element in the i-th row and j-th column of the transpose matrix is ​​equal to the element in the j-th row and i-th column of the original matrix A. When calculating the difference metric, the difference between matrix A and matrix A_T is calculated element-wise, and the sum of the absolute values ​​of all differences is taken as the difference metric. In a specific case, assuming an industrial network contains 4 monitoring nodes, the first row of the dynamic correlation matrix A has the following four elements: 0.8, 0.5, 0.3, 0.1; the second row has: 0.4, 0.9, 0.6, 0.2; the third row has: 0.2, 0.5, 0.7, 0.4; and the fourth row has: 0.1, 0.3, 0.5, 0.6. After transposing the matrix, we obtain matrix A_T. The four elements in the first row are 0.8, 0.4, 0.2, and 0.1, and the four elements in the second row are 0.5, 0.9, 0.5, and 0.3. When calculating the difference metric, the absolute value of the difference between the elements in the first row and second column is the absolute value of 0.5 - 0.4, which is 0.1. The absolute value of the difference between the elements in the first row and third column is the absolute value of 0.3 - 0.2, which is also 0.1. The absolute values ​​of the differences between all off-diagonal elements are calculated sequentially and summed to obtain a difference metric of 0.8. The preset symmetry threshold is set to 0.5. Since the difference metric of 0.8 exceeds the threshold of 0.5, symmetry processing is required. Symmetry processing is achieved by averaging element-wise. Specifically, the elements at corresponding positions are added and then divided by 2 to obtain a symmetric correlation matrix S. The element in the first row and second column is (0.5 + 0.4) / 2 = 0.45, and the element in the second row and first column is also 0.45, ensuring the symmetry of the matrix.

[0102] After symmetricizing the incidence matrix, iterative power transformation is needed to identify its inherent structural features. This process involves repeatedly multiplying the symmetric incidence matrix by the initial vector and observing the vector's changing trend. An initial vector *v* is selected, with dimensions N×1, and all elements can be initialized with equal positive values. The symmetric incidence matrix *S* is multiplied by vector *v* to obtain a new vector *w*. The length of *w* is calculated by taking the square root of the sum of the squares of its elements. Each element of the new vector *w* is divided by its length for normalization, resulting in a normalized vector *v_new*. This matrix-vector multiplication and normalization operation is repeated until the difference between the normalized vectors obtained from two consecutive iterations is less than a preset convergence threshold. At this point, the vector converges to a specific direction of the symmetric incidence matrix; the scaling factor corresponding to this direction is the eigenvalue, and the converged vector is the corresponding eigenvector. In the above 4-node case, the initial vector is selected with all elements being 1. After 12 iterations, the vector converges. At this time, the scaling factor is calculated to be 1.75, and the corresponding four elements of the feature vector are 0.52, 0.58, 0.47, and 0.38, respectively.

[0103] To obtain the complete set of eigenvalues ​​and eigenvectors, all N eigenvalues ​​need to be identified. After obtaining the first eigenvalue and eigenvector, a projection elimination operator is constructed. This operator is used to eliminate the influence of the direction corresponding to the identified eigenvector. The projection elimination operator is constructed by multiplying the eigenvector with its transpose. The residual matrix is ​​obtained by subtracting the projection elimination operator from the symmetric incidence matrix. The iterative power transformation process is repeated on the residual matrix to identify the second eigenvalue and eigenvector. This process is repeated until all N eigenvalues ​​and their corresponding eigenvectors are identified. In the 4-node case, four eigenvalues ​​were identified: 1.75, 0.92, 0.63, and 0.28, corresponding to four eigenvectors that form the eigenvector set.

[0104] After obtaining the eigenvalue set, it needs to be sorted in descending order of numerical value. In the 4-node case, the sorted eigenvalue sequence is 1.75, 0.92, 0.63, and 0.28. When calculating the cumulative energy contribution rate sequence, the squares of all eigenvalues ​​are summed to obtain the total energy value. The squares of the four eigenvalues ​​are 3.0625, 0.8464, 0.3969, and 0.0784, respectively, with a total energy value of 4.3842. The energy contribution rate of the first eigenvalue is 3.0625 / 4.3842, which equals 0.6988, and the cumulative energy contribution rate is 0.6988. The energy contribution rate of the second eigenvalue is 0.8464 / 4.3842, which equals 0.1930, and the cumulative energy contribution rate is 0.6988 + 0.1930, which equals 0.8918. The preset energy percentage threshold is set to 0.90, therefore, the first two eigenvalues ​​and their corresponding eigenvectors need to be selected.

[0105] The two selected eigenvectors are arranged column-wise to form an eigenma matrix V, which has a dimension of 4×2. The transpose of eigenma matrix V, V_T, is calculated, with a dimension of 2×4. Multiplication of eigenma matrix V with its transpose V_T yields the spectral projection operator P, which has a dimension of 4×4. During the calculation of the spectral projection operator P, the element in the i-th row and j-th column of the resulting matrix is ​​equal to the sum of the products of all elements in the i-th row of eigenma matrix V and the corresponding elements in the j-th column of the transpose V_T. In the 4-node case, the element in the first row and first column of the spectral projection operator P is calculated by multiplying the first element of the first eigenvector by itself and then multiplying the first element of the second eigenvector by itself.

[0106] The spectral projection operator P is linearly transformed with the symmetric incidence matrix S, specifically through matrix multiplication. First, the product of P and S is calculated to obtain an intermediate matrix M. Then, the product of M and P yields the final propagation operator L. Propagation operator L preserves the structural information of the main energy contribution directions in the symmetric incidence matrix while filtering out noise interference from secondary directions, thus achieving an accurate characterization of the propagation characteristics of industrial network security situations. This propagation operator can be used for subsequent situation evolution prediction and risk propagation path analysis.

[0107] In one optional implementation, the semantically enhanced situation representation is propagated through multiple iterations using the propagation operator. In each iteration, the periodic baseline component is used as a propagation stability constraint, and the decomposed situation features are used as a propagation stimulus source, including:

[0108] The distribution of the semantically enhanced situational representation across multiple monitoring nodes is used as an iterative state vector.

[0109] The amplitude of the periodic baseline component corresponding to each monitoring node at the current moment is extracted and mapped to a preset constraint interval through nonlinear saturation transformation to obtain the stability constraint coefficient of the monitoring node. The stability constraint coefficients of all monitoring nodes are arranged into a diagonal matrix to form a propagation stability constraint. The energy intensity of the transient disturbance component in the decomposed situation feature corresponding to each monitoring node at the current moment is extracted. The energy intensity is integrated in the time domain to obtain the cumulative excitation intensity of the monitoring node. The cumulative excitation intensities of all monitoring nodes are arranged into a column vector to form a propagation excitation source.

[0110] The propagated state vector is obtained by convolving the propagation operator with the iterative state vector. An adaptive threshold pruning transformation based on the propagation stability constraint is applied to the propagated state vector to perform nonlinear compression mapping on the state components outside the dynamic boundary range of the propagation stability constraint. The vector is then combined with the propagation excitation source to obtain the excitation propagation state.

[0111] The distribution of semantically enhanced situational awareness across multiple monitoring nodes is transformed into a state vector form required for iterative computation. Assuming the system contains 256 monitoring nodes, each node's semantically enhanced situational awareness is a 512-dimensional feature vector. These feature values ​​are arranged in spatial topological order to construct a 131072-dimensional iterative state vector. Each element of this vector corresponds to the situational strength value of a specific monitoring node in a specific semantic dimension. The i-th element to the i+511-th element corresponds to the complete feature expression of the first monitoring node, the 512th element to the 1023rd element corresponds to the feature expression of the second monitoring node, and so on, completing the state vectorization encoding of all nodes.

[0112] The amplitude data of the periodic baseline component corresponding to each monitoring node at the current calculation time is extracted. Taking monitoring node No. 37 as an example, its periodic baseline component amplitude at the current time is 0.682, which reflects the basic stability level of the node's situational characteristics. A nonlinear saturation transformation is applied to this amplitude. The amplitude input value is processed through a hyperbolic tangent mapping function. This mapping function first linearly amplifies the input value, setting the amplification factor to 2.5, obtaining an intermediate value of 1.705. Then, a saturation function compresses the intermediate value to a preset constraint range of -1 to +1. After the transformation, the stability constraint coefficient of this monitoring node is obtained as 0.936, which characterizes the acceptable range of state changes during propagation. The same nonlinear saturation transformation operation is performed on all 256 monitoring nodes, obtaining 256 stability constraint coefficients. These coefficients range from 0.821 to 0.978. Arrange these 256 constraint coefficients on the main diagonal of the diagonal matrix to form a 256x256 diagonal structure matrix. All elements in the matrix except those on the main diagonal are set to zero, forming a matrix representation of the propagation stability constraints.

[0113] Transient disturbance component data were extracted from the decomposed situational features corresponding to each monitoring node. The transient disturbance component of monitoring node No. 37 at the current moment is represented by a feature vector containing 512 elements. The energy intensity of this vector was calculated using the square root method. The 512 elements were squared and summed to obtain a total value of 146.73. The square root of the sum was then calculated to obtain an energy intensity value of 12.11. Time-domain integration was performed on this energy intensity value. In the 20 consecutive time steps prior to the current moment, the historical energy intensity sequence of this node was 11.84, 11.97, 12.05, and 12.11. The trapezoidal integral rule was used to accumulate the energy intensity in these 20 time steps, with a time step interval of 0.1 seconds. The integration result yielded a cumulative excitation intensity of 23.87 for this monitoring node. The same energy intensity extraction and time-domain integration were performed on all 256 monitoring nodes to obtain 256 cumulative excitation intensity values, which reflect the cumulative effect of transient disturbances on each node. The 256 cumulative excitation intensity values ​​were arranged into a 256-dimensional column vector according to the node number. This column vector constitutes the vector representation of the propagating excitation source, and the value of the 37th element in the vector is 23.87.

[0114] The propagation operator is constructed as a 256x256 sparse matrix structure. Each row of the matrix contains several non-zero elements representing the propagation connections between the corresponding monitoring node and its neighboring nodes. Taking the matrix row corresponding to monitoring node 37 as an example, the elements in columns 36, 37, 38, and 52 are 0.15, 0.50, 0.20, and 0.15, respectively, while the elements in the remaining columns are zero, indicating that this node mainly has propagation coupling relationships with its four neighboring nodes. The propagation operator matrix is ​​multiplied by a 131072-dimensional iterative state vector. For the 512 consecutive elements in the iterative state vector corresponding to node 37, a weighted sum is performed with the corresponding weights in row 37 of the propagation operator matrix. The state vector contributions from nodes 36, 38, and 52 are also considered to complete the propagation calculation for the 512 feature dimensions of this node. The same convolution operation is performed on the feature vectors of all 256 nodes to generate a 131072-dimensional post-propagation state vector.

[0115] An adaptive threshold pruning transformation based on propagation stability constraints is applied to the propagated state vector. The dynamic boundary range is determined according to the stability constraint coefficients of each node in the diagonal structure matrix. The constraint coefficient of monitoring node No. 37 is 0.936. The baseline state mean of this node before propagation is set to 0.428. Therefore, the dynamic boundary range of this node is the baseline state mean minus the constraint coefficient, resulting in a lower bound of -0.508, and the baseline state mean plus the constraint coefficient, resulting in an upper bound of 1.364. The 512 eigenvalues ​​corresponding to node No. 37 in the propagated state vector are then examined. The 89th dimension eigenvalue was found to be 1.576, exceeding the upper bound, and the 203rd dimension eigenvalue was -0.691, exceeding the lower bound. Nonlinear compression mapping was applied to the eigenvalues ​​exceeding the upper bound, calculating the excess magnitude as 0.212. Multiplying this excess magnitude by a decay factor of 0.3 yielded a compressed excess of 0.064, and the compressed eigenvalue was updated to 1.428. For the eigenvalues ​​exceeding the lower bound, the excess magnitude was calculated to be 0.183, and similarly multiplied by a decay factor of 0.3, yielded a compressed excess of 0.055, and the compressed eigenvalue was updated to -0.563. Boundary detection and nonlinear compression mapping operations were performed on the feature vectors of all 256 nodes to complete the adaptive threshold pruning transformation.

[0116] The propagated state vector after pruning is combined with the propagation excitation source column vector. For monitoring node 37, its pruned 512-dimensional feature vector is fused with the corresponding cumulative excitation intensity of 23.87. Specifically, the cumulative excitation intensity value is distributed across the 512 feature dimensions according to the excitation allocation weight vector. The excitation allocation weight vector is determined based on the sensitivity of each feature dimension to transient disturbances. In the weight vector of node 37, the weight of the 89th dimension is 0.032, and the excitation increment obtained for this dimension is 0.764. This increment is superimposed on the pruned feature value of 1.428 to obtain the excitation value of 2.192. The excitation increment calculation and superposition operation are performed on all 512 feature dimensions of this node to complete the generation of the propagated state after excitation for a single node. The same excitation fusion process is repeated for all 256 monitoring nodes to form a complete 131072-dimensional propagated state vector after excitation. This vector serves as the input state vector for the next iteration step, and the multi-step iterative propagation calculation is performed repeatedly until the preset iteration termination condition is reached.

[0117] In one alternative implementation, and by introducing a decay factor during the iteration process to simulate the diffusion and decay of the anomalous effect over time and space, the following is included:

[0118] Set the initial value of the iteration step counter to zero, and record the node that first shows an abnormal response as the set of abnormal source nodes;

[0119] The current value of the iteration step counter is used as the exponent to calculate the time decay component. The topological distance between each monitoring node and the set of abnormal source nodes is calculated and a Gaussian radial basis transformation is applied to obtain the spatial decay component. The time decay component and the spatial decay component are subjected to a Hadamard product to form an adaptive decay weight distribution. The local variance of the propagation state after excitation is calculated, and the monitoring nodes whose local variance exceeds a preset fluctuation threshold are identified as high fluctuation region nodes. The decay weights corresponding to the high fluctuation region nodes are attenuated and compensated. The adaptive decay weight distribution after attenuation compensation is fused with the propagation state after excitation node by node to obtain the propagation state after decay. The iteration step counter is then incremented by one.

[0120] Calculate the relative rate of change between the decayed propagation state and the decayed propagation state of the previous iteration step. If the relative rate of change is lower than a preset convergence threshold or the iteration step counter reaches a preset maximum iteration step, the decayed propagation state is used as the industrial network security situation prediction result.

[0121] In the iterative prediction of industrial network security posture, the initialization phase requires establishing an iteration step counter and setting it to zero. This counter tracks the time steps experienced during the diffusion and evolution process, providing a time-dimensional baseline parameter for subsequent decay calculations. Simultaneously, the real-time response status of all monitored nodes is scanned to identify the node location where abnormal characteristics are first detected. Anomaly identification is based on indicators such as sudden traffic changes, protocol deviations, or abnormal behavior. The identifiers of these nodes are stored in an anomaly source node set, which may contain a single node or multiple nodes, depending on whether the attack mode is a single-point intrusion or coordinated penetration. For example, in an industrial control system, if two SCADA servers, numbered N15 and N23, simultaneously exhibit abnormal login attempt frequencies at timestamp T0, the anomaly source node set is initialized as {N15, N23}.

[0122] The current value of the iteration step counter is used as the control parameter for the decay rate. During the calculation, the base of the natural logarithm is used as the base, and the iteration step count is multiplied by a negative decay coefficient and then used as the exponent for power operation. The decay coefficient typically ranges from 0.05 to 0.2; a larger coefficient value indicates that the impact of the anomaly weakens rapidly over time. Assuming the decay coefficient is set to 0.1, when the iteration step counter value is 3, the time decay component is calculated as the base of the natural logarithm raised to the power of -0.3, with a value of approximately 0.74, indicating that the intensity of the anomaly impact retains 74% of its original value. This mechanism simulates the physical characteristic of the impact of a security incident gradually weakening over time in reality.

[0123] The construction of the spatial attenuation component relies on network topology analysis. Each monitoring node is traversed, and the shortest path hop count to all members of the anomaly source node set is calculated. If a node has different distances to multiple anomaly sources, the minimum value is selected as the node's characteristic distance. Taking node N08 as an example, if its topological distance to N15 is 2 hops and its topological distance to N23 is 4 hops, then 2 is selected as the characteristic distance of N08. Subsequently, a Gaussian radial basis transformation is applied to the characteristic distance. This transformation requires setting an influence radius parameter, typically taken as 1 / 3 to 1 / 2 of the network's average diameter. Specifically, the square of the characteristic distance is divided by negative twice the square of the influence radius, and then the quotient is subjected to natural exponentiation. If the influence radius is set to 3 hops, the spatial attenuation component of node N08 is calculated as the base of the natural logarithm of -4 divided by the 18th power, approximately 0.80, indicating that the spatial intensity of the anomaly source influence at this location retains 80% of the original value.

[0124] An adaptive decay weight distribution is formed through the Hadamard product operation. This involves multiplying the temporal decay component and spatial decay component element-wise for each monitoring node. This operation ensures that both temporal evolution and spatial propagation decay effects are considered simultaneously. Continuing the previous example, the decay weight of node N08 at iteration step 3 is the product of 0.74 and 0.80, which is 0.592. After all nodes in the network complete this calculation, a weight vector of equal length to the number of nodes is formed. This vector describes the distribution intensity map of the anomaly's impact throughout the network. Nodes closer to the anomaly source and with earlier occurrences receive higher weight values, while nodes farther from the anomaly source or with later occurrences have weight values ​​close to zero.

[0125] A sliding window mechanism is used to statistically analyze the local variance of the state values ​​of each node and its neighboring nodes. The window size is typically set to the node's degree plus one, ensuring that it includes the node itself and all its directly connected neighbors. The arithmetic mean of the propagated state values ​​of each node within the window is calculated. Then, the squared difference between each state value and the mean is calculated. The sum of all squared differences is divided by the number of nodes within the window to obtain the local variance. Taking node N12 as an example, assuming its own state value is 0.85, and the state values ​​of its three neighboring nodes are 0.78, 0.92, and 0.81 respectively, the mean within the window is 0.84, and the four squared differences are 0.0001, 0.0036, 0.0064, and 0.0009 respectively. The local variance is the sum of these four values ​​divided by 4, resulting in 0.00275.

[0126] The identification of nodes in high-fluctuation areas is based on the comparison between local variance and a preset fluctuation threshold. The fluctuation threshold is determined according to the fluctuation statistical characteristics of historical network operation data, and is usually set to three to five times the standard deviation of the mean local variance during normal operation. If the mean local variance of an industrial network during normal operation is 0.0015 and the standard deviation is 0.0005, then the fluctuation threshold can be set to 0.003. Node N12 is not marked when its local variance is below this threshold (0.00275), while node N19 is identified as a node in a high-fluctuation area if its local variance reaches 0.0042. Such nodes are located at the forefront of abnormal propagation or are critical hubs in the network topology, requiring enhanced monitoring accuracy.

[0127] The attenuation compensation mechanism adjusts the weights of nodes in high-fluctuation regions by introducing a compensation factor to offset part of the attenuation effect. The compensation factor is calculated by taking the square root of the ratio of the node's local variance to the fluctuation threshold to avoid over-amplification. For example, node N19 has a local variance of 0.0042 and a threshold of 0.003, with a ratio of 1.4. Taking the square root yields a compensation factor of approximately 1.18. Multiplying the node's original adaptive attenuation weight by the compensation factor, if the original weight was 0.45, the compensated weight is adjusted to 0.531. This mechanism ensures that regions with abnormally active evolution receive sufficient attention, preventing the attenuation model from excessively suppressing the propagation of potential threats.

[0128] The node-by-node fusion operation multiplies the adaptive attenuation weight distribution after attenuation compensation with the corresponding element-wise propagation state after excitation. The post-excitation state value of each monitoring node is multiplied by its corresponding adjusted weight value to generate the attenuated state value for that node. Assuming node N19 has a post-excitation state value of 0.68 and a compensated weight of 0.531, its attenuated state value is calculated to be 0.361. After fusion of all nodes in the network, a new state vector is formed, namely the attenuated propagation state. This state reflects the distribution of the abnormal impact after spatiotemporal attenuation and fluctuation compensation. After completing this round of iteration calculation, the iteration step counter is incremented by one to prepare for the next iteration.

[0129] The relative rate of change is calculated to determine whether the iterative process has reached a stable state. It extracts the decayed propagation state vector of the current iteration step and the decayed propagation state vector of the previous iteration step, calculates the absolute value of the difference between corresponding elements of the two vectors, divides it by the absolute value of the corresponding element of the previous iteration step, and takes the arithmetic mean of the quotients over all nodes to obtain the relative rate of change. If a network contains fifty nodes, and the average relative rate of change after the seventh iteration is 0.018, while the preset convergence threshold is set to 0.02, then it is determined to have reached stability. In this case, the decayed propagation state of the seventh iteration is output as the final industrial network security situation prediction result. If the relative rate of change remains above the threshold but the iteration step counter reaches the preset maximum iteration step count, such as twenty steps, the iteration is terminated and the current state is output to avoid unlimited consumption of computational resources. This prediction result provides the security response team with the probability distribution of attacks affecting each monitored node within a future time window, supporting the optimized allocation of defense resources.

[0130] This invention provides an industrial network security situation prediction system based on a generative large model, the system comprising:

[0131] The first unit is used to acquire network behavior data of multiple monitoring nodes in the industrial network in a time series, and to extract multi-dimensional features from the network behavior data to obtain an industrial network security feature vector.

[0132] The second unit is used to decompose the industrial network security feature vector into periodic baseline components and transient disturbance components through frequency domain transformation, and to filter the transient disturbance components through sparsity constraints to obtain the decomposed situation features.

[0133] The third unit is used to perform semantic space mapping between the periodic baseline components and the decomposed situation features using a generative large model to obtain a semantically enhanced situation representation. Based on the decomposed situation features, it calculates the security anomaly correlation between multiple monitoring nodes and constructs a dynamic correlation matrix between industrial network nodes.

[0134] The fourth unit is used to construct a propagation operator for industrial network security situation based on the dynamic correlation matrix. The propagation operator is used to propagate the semantically enhanced situation representation in multiple steps. In each step of the iteration, the periodic baseline component is used as a propagation stability constraint, the decomposed situation features are used as a propagation excitation source, and an attenuation factor is introduced in the iteration process to simulate the diffusion and attenuation of the abnormal influence over time and space, so as to obtain the prediction result of industrial network security situation within the future time window.

[0135] A third aspect of the present invention provides an electronic device, comprising:

[0136] processor;

[0137] Memory used to store processor-executable instructions;

[0138] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0139] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0140] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0141] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for predicting industrial network security situation based on generative large models, characterized in that, include: The network behavior data of multiple monitoring nodes in the industrial network are acquired in a time series, and multi-dimensional feature extraction is performed on the network behavior data to obtain the industrial network security feature vector. The industrial network security feature vector is decomposed into periodic baseline components and transient disturbance components through frequency domain transformation. The transient disturbance components are then filtered using sparsity constraints to obtain the decomposed situational features, including: The industrial network security feature vector is mapped from the time domain to the frequency domain to obtain the frequency domain feature. The spectral coefficients in the frequency domain feature are arranged in ascending order of frequency and subjected to exponential fitting to obtain the energy attenuation curve. Calculate the fitting residual between the actual energy value at each frequency position and the fitted value of the energy attenuation curve, perform a second-order difference operation on the fitting residual to obtain the residual curvature, identify the frequency position where the sign of the residual curvature changes as the frequency band boundary point, and divide the spectral coefficients before and after the frequency band boundary point into dominant frequency bands and non-dominant frequency bands respectively. The spectral coefficients of the dominant frequency band and the non-dominant frequency band are respectively subjected to inverse frequency domain transformation to obtain the periodic baseline component and the initial transient disturbance component; The initial transient disturbance component is decomposed into multiple time scales. The amplitude sequence of the disturbance component at each time scale is calculated. The sparsity is obtained by calculating the ratio of the zero norm of the amplitude sequence to the sequence length. The time scale with the largest sparsity is selected as the dominant scale. The amplitude of the time corresponding to the non-zero position of the amplitude at the dominant scale is extracted in the initial transient disturbance component. After sorting in descending order, the disturbance signal corresponding to the time corresponding to the preset proportion is retained to obtain the filtered transient disturbance component. This component is combined with the periodic baseline component to form the decomposed situation feature. A generative large model is used to perform semantic space mapping between the periodic baseline components and the decomposed situation features to obtain a semantically enhanced situation representation. Based on the decomposed situation features, the correlation of security anomalies among multiple monitoring nodes is calculated, and a dynamic correlation matrix among industrial network nodes is constructed. Based on the dynamic correlation matrix, a propagation operator for industrial network security situation is constructed. The semantically enhanced situation representation is propagated through multiple iterations using the propagation operator. In each iteration, the periodic baseline component is used as a propagation stability constraint, and the decomposed situation features are used as a propagation excitation source. In the iteration process, an attenuation factor is introduced to simulate the diffusion and attenuation of abnormal influences over time and space, thereby obtaining the prediction result of industrial network security situation within the future time window.

2. The method according to claim 1, characterized in that, Multi-dimensional feature extraction is performed on the network behavior data to obtain the industrial network security feature vector, which includes: From the network behavior data, the interaction patterns of industrial network communication protocols in the protocol layer are extracted to obtain protocol layer features, and the statistical periods of network data transmission in the traffic layer are extracted to obtain traffic layer features. The protocol layer features are sequence encoded, and the protocol interaction behaviors at different times are arranged in chronological order to form a protocol interaction sequence. The protocol type identifier and interaction state identifier are extracted from each protocol interaction behavior in the protocol interaction sequence and mapped to a vector space to obtain a protocol behavior vector. The protocol behavior vectors are combined in chronological order to form a time-series vector. The flow layer features are distributed and modeled. The flow statistics within a preset time window are statistically analyzed. The skewness and kurtosis of the flow statistics are calculated to obtain statistical moments. Based on the statistical moments, the distribution shape parameters of the flow statistics are determined. The time-series vector and the distribution pattern parameter are concatenated to form the industrial network security feature vector that integrates the time-series interaction mode and the traffic distribution pattern.

3. The method according to claim 1, characterized in that, Based on the decomposed situation characteristics, the correlation of security anomalies among multiple monitoring nodes is calculated, and a dynamic correlation matrix among industrial network nodes is constructed, including: Extract the filtered transient disturbance components from any two monitoring nodes and calculate the time-series cross-response characteristics in the time dimension. Identify the peak amplitude and time offset corresponding to the peak position of the time-series cross-response characteristics. Use the time offset as the anomaly propagation delay. Perform a nonlinear mapping between the peak amplitude and the reciprocal of the anomaly propagation delay to obtain the time-series coupling strength. The energy distribution of the periodic baseline components of any two monitoring nodes in the frequency domain is calculated respectively. By constructing the frequency domain coherence metric matrix between the two energy distributions, the frequency domain energy coupling spectrum is obtained. The frequency domain energy coupling spectrum is adaptively weighted and integrated over the entire frequency band and mapped to the zero-to-one interval through hyperbolic tangent transformation to obtain the frequency domain coupling strength. The temporal coupling strength and the frequency domain coupling strength are adaptively fused to obtain the security anomaly correlation between the two monitoring nodes. By traversing all node pairs in the monitoring nodes, using multiple monitoring nodes as row and column indices, and using the security anomaly correlation between each node pair as matrix elements at the corresponding row and column positions, the dynamic correlation matrix between industrial network nodes is constructed.

4. The method according to claim 1, characterized in that, The propagation operators for constructing industrial network security posture based on the aforementioned dynamic correlation matrix include: Calculate the difference metric between the dynamic correlation matrix and the corresponding transpose matrix. When the difference metric exceeds a preset symmetry threshold, the dynamic correlation matrix and the transpose matrix are averaged element-wise to obtain a symmetric correlation matrix. Otherwise, the dynamic correlation matrix is ​​used as the symmetric correlation matrix. By applying an iterative power transformation to the symmetric correlation matrix, the scaling factors of the symmetric correlation matrix in different directions are identified as a set of eigenvalues. For each eigenvalue in the set of eigenvalues, the vector in which the symmetric correlation matrix remains unchanged in the scaling direction corresponding to the eigenvalue is determined as an eigenvector. The eigenvalue set is sorted in descending order of eigenvalues. The cumulative energy contribution rate sequence of the eigenvalue set is calculated. The eigenvalues ​​that reach a preset energy percentage and their corresponding eigenvectors in the cumulative energy contribution rate sequence are selected. The corresponding eigenvectors are arranged in columns and multiplied with their own transpose to construct a spectral projection operator. The spectral projection operator is used to perform a linear transformation on the symmetric correlation matrix. The linearly transformed symmetric correlation matrix is ​​used as the propagation operator for industrial network security situation.

5. The method according to claim 1, characterized in that, The semantically enhanced situation representation is propagated through a multi-step iterative process using the propagation operator. In each iteration, the periodic baseline component is used as a propagation stability constraint, and the decomposed situation features are used as the propagation excitation source. The distribution of the semantically enhanced situational representation across multiple monitoring nodes is used as an iterative state vector. The amplitude of the periodic baseline component corresponding to each monitoring node at the current time is extracted and mapped to a preset constraint interval through nonlinear saturation transformation to obtain the stability constraint coefficient of the monitoring node. The stability constraint coefficients of all monitoring nodes are arranged into a diagonal matrix to form a propagation stability constraint. Extract the energy intensity of the transient disturbance component in the decomposed situation features corresponding to each monitoring node at the current moment, perform time-domain integration on the energy intensity to obtain the cumulative excitation intensity of the monitoring node, and form a column vector of the cumulative excitation intensities of all monitoring nodes to form a propagation excitation source; The propagated state vector is obtained by convolving the propagation operator with the iterative state vector. An adaptive threshold pruning transformation based on the propagation stability constraint is applied to the propagated state vector to perform nonlinear compression mapping on the state components outside the dynamic boundary range of the propagation stability constraint. The vector is then combined with the propagation excitation source to obtain the excitation propagation state.

6. The method according to claim 5, characterized in that, The iteration process introduces a decay factor to simulate the diffusion and decay of anomalous effects over time and space, including: Set the initial value of the iteration step counter to zero, and record the node that first shows an abnormal response as the set of abnormal source nodes; The current value of the iteration step counter is used as the exponent to calculate the time decay component. The topological distance between each monitoring node and the set of abnormal source nodes is calculated and Gaussian radial basis transformation is applied to obtain the spatial decay component. The time decay component and the spatial decay component are subjected to Hadamard product operation to form an adaptive decay weight distribution. Calculate the local variance of the propagation state after excitation, identify monitoring nodes whose local variance exceeds a preset fluctuation threshold as high fluctuation region nodes, perform attenuation compensation on the attenuation weights corresponding to the high fluctuation region nodes, fuse the adaptive attenuation weight distribution after attenuation compensation with the propagation state after excitation node by node to obtain the propagation state after attenuation, and increment the iteration step counter by one. Calculate the relative rate of change between the decayed propagation state and the decayed propagation state of the previous iteration step. If the relative rate of change is lower than a preset convergence threshold or the iteration step counter reaches a preset maximum iteration step, the decayed propagation state is used as the industrial network security situation prediction result.

7. An industrial network security situation prediction system based on a generative large model, used to implement the method as described in any one of claims 1-6, characterized in that, include: The first unit is used to acquire network behavior data of multiple monitoring nodes in the industrial network in a time series, and to extract multi-dimensional features from the network behavior data to obtain an industrial network security feature vector. The second unit is used to decompose the industrial network security feature vector into periodic baseline components and transient disturbance components through frequency domain transformation, and to filter the transient disturbance components through sparsity constraints to obtain the decomposed situation features. The third unit is used to perform semantic space mapping between the periodic baseline components and the decomposed situation features using a generative large model to obtain a semantically enhanced situation representation. Based on the decomposed situation features, it calculates the security anomaly correlation between multiple monitoring nodes and constructs a dynamic correlation matrix between industrial network nodes. The fourth unit is used to construct a propagation operator for industrial network security situation based on the dynamic correlation matrix. The propagation operator is used to propagate the semantically enhanced situation representation in multiple steps. In each step of the iteration, the periodic baseline component is used as a propagation stability constraint, the decomposed situation features are used as a propagation excitation source, and an attenuation factor is introduced in the iteration process to simulate the diffusion and attenuation of the abnormal influence over time and space, so as to obtain the prediction result of industrial network security situation within the future time window.

8. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images