CN116736781B - Safety state monitoring method and device for industrial automation control equipment - Google Patents
Safety state monitoring method and device for industrial automation control equipment Download PDFInfo
- Publication number
- CN116736781B CN116736781B CN202311021790.7A CN202311021790A CN116736781B CN 116736781 B CN116736781 B CN 116736781B CN 202311021790 A CN202311021790 A CN 202311021790A CN 116736781 B CN116736781 B CN 116736781B
- Authority
- CN
- China
- Prior art keywords
- data
- distribution
- safety state
- fitting
- control equipment
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 65
- 238000012544 monitoring process Methods 0.000 title claims abstract description 58
- 238000009826 distribution Methods 0.000 claims abstract description 148
- 238000005315 distribution function Methods 0.000 claims abstract description 97
- 230000001186 cumulative effect Effects 0.000 claims abstract description 44
- 230000002159 abnormal effect Effects 0.000 claims abstract description 34
- 238000012163 sequencing technique Methods 0.000 claims abstract description 12
- 230000008447 perception Effects 0.000 claims abstract description 10
- 230000006870 function Effects 0.000 claims description 39
- 238000007476 Maximum Likelihood Methods 0.000 claims description 34
- 238000004364 calculation method Methods 0.000 claims description 28
- 238000001276 Kolmogorov–Smirnov test Methods 0.000 claims description 26
- 238000004458 analytical method Methods 0.000 claims description 15
- 238000012806 monitoring device Methods 0.000 claims 1
- 230000008569 process Effects 0.000 description 11
- 230000006399 behavior Effects 0.000 description 4
- 238000004519 manufacturing process Methods 0.000 description 3
- 206010000117 Abnormal behaviour Diseases 0.000 description 2
- 238000012550 audit Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 230000006854 communication Effects 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 230000005856 abnormality Effects 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000000528 statistical test Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 238000005728 strengthening Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B19/00—Programme-control systems
- G05B19/02—Programme-control systems electric
- G05B19/04—Programme control other than numerical control, i.e. in sequence controllers or logic controllers
- G05B19/042—Programme control other than numerical control, i.e. in sequence controllers or logic controllers using digital processors
- G05B19/0428—Safety, monitoring
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B2219/00—Program-control systems
- G05B2219/20—Pc systems
- G05B2219/24—Pc safety
- G05B2219/24024—Safety, surveillance
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/02—Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Automation & Control Theory (AREA)
- Testing And Monitoring For Control Systems (AREA)
Abstract
The application discloses a safety state monitoring method and device of industrial automation control equipment, comprising the steps of determining a target data distribution model matched with the normal safety state of industrial control equipment before entering a safety state monitoring mode, and determining an abnormal state perception threshold corresponding to the target data distribution model; after entering a safety state monitoring mode, collecting real-time network flow data characteristic values of industrial control equipment, rapidly sequencing, and calculating a corresponding real-time experience cumulative distribution function; calculating the coincidence degree of the corresponding probability distribution function and the real-time experience cumulative distribution function under the normal safety state of the industrial control equipment to obtain a check deviation value; and when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state. The method and the device for monitoring the safety state of the industrial automation control equipment provided by the application have low operation complexity, and are suitable for being deployed on the industrial environment site with limited computing resources to effectively monitor the safety state of the industrial automation control equipment.
Description
Technical Field
The application relates to the technical field of data processing, in particular to a method and a device for monitoring the safety state of industrial automation control equipment.
Background
Industrial automation control devices are an important component in industrial control systems for monitoring and controlling various industrial processes. They are commonly used to control critical infrastructure of plant equipment, energy systems, traffic systems, and the like. The safety threat faced by industrial automation control equipment terminals can lead to production interruption, equipment damage and even casualties, so the safety state monitoring of the industrial automation control terminals is of great importance.
The safety state monitoring of the industrial automation control equipment means that the industrial control terminal is monitored and analyzed in real time to acquire the safety state information of the industrial control terminal. The method comprises the steps of monitoring hardware, an operating system, industrial control software, network connection and the like of the industrial control terminal so as to detect potential security holes, malicious attacks and abnormal behaviors, discover and deal with security threats in time, and ensure reliable operation and production security of the industrial control system.
Therefore, how to ensure the safety state monitoring of the industrial automation control equipment is a technical problem to be solved urgently by those skilled in the art.
Disclosure of Invention
The application provides a method and a device for monitoring the safety state of industrial automation control equipment, which have low operation complexity and are suitable for being deployed on an industrial environment site with limited computing resources to effectively monitor the safety state of industrial control equipment.
In order to solve the above technical problems, an embodiment of the present application provides a method for monitoring a safety state of an industrial automation control device, including:
before entering a safety state monitoring mode, collecting a network flow data characteristic value of industrial control equipment in a normal safety state, and dividing the network flow data characteristic value into a fitting sample and a checking sample;
selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples one by adopting a maximum likelihood estimation method to obtain corresponding fitting probability distribution functions;
rapidly sequencing the check samples by adopting a bucket-tree two-dimensional indexing method to obtain sequenced first data, and calculating an experience cumulative distribution function corresponding to the first data;
calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test;
determining a target data distribution model matched with the normal safety state of industrial control equipment, and determining an abnormal state perception threshold corresponding to the target data distribution model;
after entering a safety state monitoring mode, acquiring a real-time network flow data characteristic value of the industrial control equipment;
inserting the real-time network flow data characteristic values into a barrel-tree two-dimensional index for quick sequencing to obtain sequenced second data, and calculating a real-time experience cumulative distribution function corresponding to the second data;
adopting K-S test to calculate the coincidence degree of the probability distribution function and the real-time experience cumulative distribution function corresponding to the normal safety state of the industrial control equipment, and obtaining a check deviation value;
and when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state.
As one preferable mode, the data distribution model at least comprises a gaussian distribution model, a lognormal distribution model, a poisson distribution model, a weibull distribution model, a gamma distribution model and an exponential distribution model.
As one preferable scheme, the selecting a plurality of data distribution models, and performing fitting analysis on the fitting samples by adopting a maximum likelihood estimation method one by one to obtain corresponding fitting probability distribution functions, specifically including:
selecting a data distribution model from a plurality of data distribution models;
calculating likelihood functions corresponding to the data distribution model, wherein the likelihood functions are as follows:
wherein L (theta) is a likelihood function, X is a network flow data characteristic value set in acquisition time, and theta is a model parameter of a selected data distribution model;
taking the logarithm of the likelihood function to obtain a log likelihood function as follows:
deriving the log-likelihood function to obtain a corresponding maximum likelihood estimation parameter value;
and carrying the maximum likelihood estimation parameter value into an accumulated distribution function corresponding to the selected characteristic distribution model for calculation, wherein the calculation formula is as follows:
wherein ,for accumulating distribution functions +.>Estimating parameter values for the maximum likelihood;
repeating the steps until the maximum likelihood estimation parameter values and the accumulated distribution function corresponding to all the data distribution models are obtained.
As one preferable solution, the quick sorting of the check samples by adopting a bucket-tree two-dimensional indexing method specifically includes:
distinguishing the data in the check sample according to a barrel dividing mode, and storing the data in a corresponding barrel array;
and arranging the data in each bucket array according to a tree array structure.
As one preferable scheme, the calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to the K-S test specifically includes:
according to K-S test, respectively calculating a first parameter, a second parameter and a third parameter of a network flow data distribution characteristic value corresponding to the data stored in each barrel group, wherein the calculation formula is as follows:
wherein ,to reflect the first parameter of maximum forward K-S distance,/>To reflect the second parameter of the maximum inverse K-S distance,/a second parameter of the maximum inverse K-S distance is chosen>To reflect the third parameter of the feature quantity stored in the j-th bucket array, +.>Is->Barrel array corresponding to each interval, < > and the like>For probability distribution function +.>For the empirical cumulative distribution function, the calculation formula is as follows:
wherein ,indicate->Then 1;
and obtaining the K-S check fitting distance reflecting the fitting degree based on the maximum first parameter and the maximum second parameter.
Another embodiment of the present application provides an industrial automation control device safety state monitoring apparatus, comprising a processor configured to:
before entering a safety state monitoring mode, collecting a network flow data characteristic value of industrial control equipment in a normal safety state, and dividing the network flow data characteristic value into a fitting sample and a checking sample;
selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples one by adopting a maximum likelihood estimation method to obtain corresponding fitting probability distribution functions;
rapidly sequencing the check samples by adopting a bucket-tree two-dimensional indexing method to obtain sequenced first data, and calculating an experience cumulative distribution function corresponding to the first data;
calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test;
determining a target data distribution model matched with the normal safety state of industrial control equipment, and determining an abnormal state perception threshold corresponding to the target data distribution model;
after entering a safety state monitoring mode, acquiring a real-time network flow data characteristic value of the industrial control equipment;
inserting the real-time network flow data characteristic values into a barrel-tree two-dimensional index for quick sequencing to obtain sequenced second data, and calculating a real-time experience cumulative distribution function corresponding to the second data;
adopting K-S test to calculate the coincidence degree of the probability distribution function and the real-time experience cumulative distribution function corresponding to the normal safety state of the industrial control equipment, and obtaining a check deviation value;
and when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state.
As one preferable mode, the data distribution model at least comprises a gaussian distribution model, a lognormal distribution model, a poisson distribution model, a weibull distribution model, a gamma distribution model and an exponential distribution model.
As one preferable scheme, the selecting a plurality of data distribution models, and performing fitting analysis on the fitting samples by adopting a maximum likelihood estimation method one by one to obtain corresponding fitting probability distribution functions, specifically including:
selecting a data distribution model from a plurality of data distribution models;
calculating likelihood functions corresponding to the data distribution model, wherein the likelihood functions are as follows:
wherein L (theta) is a likelihood function, X is a network flow data characteristic value set in acquisition time, and theta is a model parameter of a selected data distribution model;
taking the logarithm of the likelihood function to obtain a log likelihood function as follows:
deriving the log-likelihood function to obtain a corresponding maximum likelihood estimation parameter value;
and carrying the maximum likelihood estimation parameter value into an accumulated distribution function corresponding to the selected characteristic distribution model for calculation, wherein the calculation formula is as follows:
wherein ,for accumulating distribution functions +.>Estimating parameter values for the maximum likelihood;
repeating the steps until the maximum likelihood estimation parameter values and the accumulated distribution function corresponding to all the data distribution models are obtained.
As one preferable solution, the quick sorting of the check samples by adopting a bucket-tree two-dimensional indexing method specifically includes:
distinguishing the data in the check sample according to a barrel dividing mode, and storing the data in a corresponding barrel array;
and arranging the data in each bucket array according to a tree array structure.
As one preferable scheme, the calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to the K-S test specifically includes:
according to K-S test, respectively calculating a first parameter, a second parameter and a third parameter of a network flow data distribution characteristic value corresponding to the data stored in each barrel group, wherein the calculation formula is as follows:
wherein ,to reflect the first parameter of maximum forward K-S distance,/>To reflect the second parameter of the maximum inverse K-S distance,/a second parameter of the maximum inverse K-S distance is chosen>To reflect the third parameter of the feature quantity stored in the j-th bucket array, +.>Is->Barrel array corresponding to each interval, < > and the like>For probability distribution function +.>For the empirical cumulative distribution function, the calculation formula is as follows:
wherein ,indicate->Then 1;
and obtaining the K-S check fitting distance reflecting the fitting degree based on the maximum first parameter and the maximum second parameter.
Compared with the prior art, the embodiment of the application has the beneficial effects that at least one of the following points is adopted:
before entering a safety state monitoring mode, fully considering the matching degree of each model in the normal safety state of the industrial control equipment, screening and determining a target data distribution model matched with the normal safety state of the industrial control equipment, determining an abnormal state sensing threshold corresponding to the target data distribution model, then entering a real-time monitoring mode, monitoring the safety state of the industrial automation control equipment, continuously calculating the network flow data characteristic distribution in the running process of the equipment, adopting K-S test to judge the coincidence degree of the network flow data characteristic distribution and the reference characteristic distribution, and if the deviation is larger than the abnormal state sensing threshold, judging that the running state is abnormal. The whole monitoring method does not need to invade the equipment to internally install a specific program to acquire the running state of the system for safety state sensing, has low operation complexity in the whole process, and is suitable for being deployed in the industrial environment site with limited computing resources for effective safety state monitoring of industrial control equipment.
Drawings
FIG. 1 is a flow chart of a method for monitoring a safety state of an industrial automation control device in one embodiment of the present application;
FIG. 2 is a flow chart of a method for monitoring safety status of an industrial automation control device in one embodiment of the application.
Detailed Description
The following description of the embodiments of the present application will be made clearly and fully with reference to the accompanying drawings, in which it is evident that the embodiments described are only some, but not all embodiments of the application, and the purpose of these embodiments is to provide a more thorough and complete disclosure of the present application. All other embodiments, which can be made by those skilled in the art based on the embodiments of the application without making any inventive effort, are intended to be within the scope of the application.
In the description of the present application, the terms "first," "second," "third," and the like are used for descriptive purposes only and are not to be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defining "a first", "a second", "a third", etc. may explicitly or implicitly include one or more such feature. In the description of the present application, unless otherwise indicated, the meaning of "a plurality" is two or more.
In the description of the present application, it should be noted that, unless explicitly specified and limited otherwise, the terms "mounted," "connected," and "connected" are to be construed broadly, and may be either fixedly connected, detachably connected, or integrally connected, for example; can be mechanically or electrically connected; can be directly connected or indirectly connected through an intermediate medium, and can be communication between two elements. The terms "vertical," "horizontal," "left," "right," "upper," "lower," and the like are used herein for descriptive purposes only and not to indicate or imply that the apparatus or element being referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore should not be construed as limiting the application. The term "and/or" as used herein includes any and all combinations of one or more of the associated listed items. The specific meaning of the above terms in the present application will be understood in specific cases by those of ordinary skill in the art.
In the description of the present application, it should be noted that all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs unless defined otherwise. The terminology used in the description of the present application is for the purpose of describing particular embodiments only and is not intended to be limiting of the application, as the particular meaning of the terms described above in the present application will be understood to those of ordinary skill in the art in the detailed description of the application.
The embodiment of the application provides a method for monitoring the safety state of industrial automation control equipment, which needs to be described in advance, wherein the conventional method (or sensing method) for monitoring the safety of an industrial control terminal comprises the following steps: 1) The system operation monitoring system is used for identifying possible security threats by collecting various volumes in the operation process of the industrial control terminal and analyzing the behavior mode, and triggering an alarm or taking corresponding defensive measures; 2) Access control: access to the key system is limited by performing access control on the industrial control terminal, potential security threat is reduced, and external attack threat is prevented by strengthening password policy, limiting privileged access, implementing firewall and access control list and the like; 3) Security audit and log analysis: through the security audit and log analysis of the industrial control terminal, the security events and abnormal behaviors in the system can be tracked and analyzed. This helps to discover potential attack activities and security threats in time and take appropriate action to respond and defend.
However, the above methods belong to the invasive sensing methods, and require the acquisition and analysis of the state by installing a program in the industrial control equipment, which causes a certain interference to the normal operation of the industrial control equipment, and cannot effectively meet the requirements of the industrial control operation site on the system operation continuity, the system stability and the system response instantaneity. In contrast to invasive methods, non-invasive awareness does not require any intervention on the target system, but rather information is obtained by passively listening and observing the behavior, interactions and communications of the target system. The non-invasive sensing method has the advantages that normal operation of the target system is not affected, and safety monitoring can be performed under the condition of high real-time requirements. However, there are some limitations in the non-invasive sensing method, and the existing method relies on passively monitoring and observing the behavior of the target system, which generally requires a large amount of data collection and analysis, and has relatively large consumption of computing resources, relatively low real-time sensing performance, and no timely effective discovery of the abnormality of the system security state.
According to the embodiment of the application, the sensing problem of the safety state is converted into the characteristic matching problem, the network flow data distribution characteristic of the industrial control equipment in the safety operation state is firstly used as the reference distribution characteristic, the deviation between the data distribution characteristic and the reference distribution characteristic of the system in the actual operation process is judged by adopting a statistical test mode to rapidly judge the safety state of the system, and the technical problems that a large amount of observation data needs to be collected, a large amount of calculation resources are needed and the response time is not timely in the traditional non-invasive safety state sensing method can be solved.
Referring to fig. 1 to fig. 2, fig. 1 is a schematic flow chart of a method for monitoring a safety state of an industrial automation control device according to one embodiment of the present application, and fig. 2 is a block flow chart of a method for monitoring a safety state of an industrial automation control device according to one embodiment of the present application, where the monitoring method specifically includes:
s1, before entering a safety state monitoring mode, collecting a network flow data characteristic value of industrial control equipment in a normal safety state, and dividing the network flow data characteristic value into a fitting sample and a verification sample;
s2, selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples by adopting a maximum likelihood estimation method one by one to obtain corresponding fitting probability distribution functions;
s3, rapidly sequencing the check samples by adopting a bucket-tree two-dimensional indexing method to obtain sequenced first data, and calculating an experience cumulative distribution function corresponding to the first data;
s4, calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test;
s5, determining a target data distribution model matched with the normal safety state of the industrial control equipment, and determining an abnormal state perception threshold corresponding to the target data distribution model;
s6, after entering a safety state monitoring mode, acquiring a real-time network flow data characteristic value of the industrial control equipment;
s7, inserting the real-time network flow data characteristic values into a barrel-tree two-dimensional index for quick sorting to obtain sorted second data, and calculating a real-time experience cumulative distribution function corresponding to the second data;
s8, adopting K-S test, and calculating the coincidence degree of the probability distribution function and the real-time experience cumulative distribution function corresponding to the normal safety state of the industrial control equipment to obtain a check deviation value;
and S9, when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state.
In order to monitor the safety state of industrial automation control equipment, the general TCP/IP protocol network data flow of the industrial control equipment is selected as the collected data. The whole implementation process comprises a network flow data distribution characteristic extraction stage of the normal safety state of the industrial control equipment and an abnormal safety state sensing stage of the industrial control equipment. Specifically, the above embodiment includes two steps, namely, a first step (corresponding to the steps S1 to S5) before entering the safety state monitoring mode, and a second step (corresponding to the steps S6 to S9) after entering the safety state monitoring mode.
1. In the network flow data distribution characteristic extraction stage of the normal safety state of the industrial control equipment, the specific steps are as follows:
and step 1, collecting characteristic values of network flow data of industrial control equipment in a normal safety state, and dividing the characteristic values into a fitting sample and a checking sample.
In step 1, the TCP/IP network data of the industrial control device is collected, and various information of the data packet including the source IP address, the destination IP address, the protocol type, the data packet length, etc. are obtained, and the content of the data packet may be further analyzed, and deeper information, such as an application layer protocol, user behavior, etc., may be extracted, where these information constitute the characteristic value of the network traffic data. Preferably, the fitting sample and the verification sample use 1:1, and the acquisition time is preferably 30 minutes.
Step 2, selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples by adopting maximum likelihood estimation one by one to obtain a fitting probability distribution function;
in the step 2, in the first step, a data distribution model is selected from a plurality of data distribution models to be selected, where the data distribution models selectable in the implementation include gaussian distribution, lognormal distribution, poisson distribution, weibull distribution, gamma distribution, exponential distribution, and the like, and the probability density function, cumulative distribution function, and related parameters corresponding to the data distribution models are described in the following table:
in the step 2, a second step of calculating likelihood functions of the data distribution model; the likelihood function is defined as follows:
wherein Representing likelihood functions +.>For collecting characteristic value set of network flow data in time, < +.>Parameters for the feature distribution model are shown in the table above. The likelihood function is given parameter->Down->Is a probability of occurrence of (a).
Data distributed independently and simultaneously->Composition, thus likelihood function->Can be decomposed into->At a given parameter->The probability product is expressed as follows:
in the step 2, in the third step, the log likelihood function takes the log to obtain a log likelihood function, which is expressed as follows:
in the fourth step of the above step 2, the log likelihood function is used for the followingDeriving and zeroing to obtainThe parameter corresponding to the maximum value is called maximum likelihood estimation parameter value +.>。
In the fifth step of the above step 2, the maximum likelihood estimation parameter value is calculatedIs carried into the accumulated distribution function CDF corresponding to the characteristic distribution model to calculate, the accumulated distribution function CDF of each characteristic distribution model is as shown in the table aboveShown.
The function is used for representing when the characteristic distribution model parameters areTime sample variable +.>Less than value->Probability of time.
In the sixth step in the above step 2, the steps 2.1 to 2.5 are repeatedly performed until all feature distribution models are calculated and />。
Step 3, rapidly sorting the check samples by adopting a barrel-tree two-dimensional index, and calculating an experience cumulative distribution function ECDF for the sorted data;
in the step 3, in the first step, the data are distinguished according to the barrel division mode and stored in the barrel array. The process is that the value range of the collected network flow data value is uniformly divided into m intervals, each interval corresponds to one barrel, and then the data value is determined and stored in a specific barrel;
for example, assuming the current 10 data, the data values are 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5.9.5 in order. The array value range is 0-10, the setting is divided into m=5 intervals, and each interval is mapped into a barrel. The ranges of the respective buckets are [0, 2), [2, 4), [4, 6), [6, 8), and [8,10 ], respectively. Each data will then be allocated into a corresponding bucket according to its size. For example, data 0.5 and 1.5 would be placed in the first bucket, 2.5 and 3.5 would be placed in the second bucket, and so on. If there is a new oneData, e.g., 3, will be placed directly into the second bucket. The main advantage of this is that the efficiency of data querying and ordering can be greatly improved by directly accessing the corresponding buckets, rather than traversing all the data. This approach would be able to address one of the needs that would otherwise be requiredThe ordering problem of the temporal complexity is reduced to +.>Is a complex time of the system.
In the step 3, in the second step, the data in the bucket is arranged according to a tree array structure, wherein the tree array size is n+1, and n is the maximum number of the data which can be stored.
The internal structure of a tree array can be seen as a special binary tree, where each node contains a range of prefix sums of elements. This binary tree can be generated by a smart encoding scheme.
In the specific implementation process, the generation mode of the binary tree inside the tree array adopts a low-order generation method: 1) First, given a tree-like array of length n, the index of the array is encoded from 1 to n. Index 0 is not used. 2) For each node in the tree array, its corresponding index value is denoted index, and the range of nodes is from [ index-lowbit (index) +1, index ]. lowbit (index) refers to the value represented by the 1 of the lowest bit of index. lowbit (index) =index & -index. 3) Starting from the root node (index 1), the left and right child nodes of each node are generated in turn. For node index, the index of the left child node is index+lowbit (index), and the index of the right child node is index+1. 4) Repeating the steps until all nodes are generated.
By the generation mode, a binary tree structure with a correct node range can be constructed. Each node stores a prefix sum of elements within a corresponding range.
Step 4, calculating the fitting degree of the fitting probability distribution function and the experience cumulative distribution function by adopting K-S test;
in the step 4, in the first step, the collected network traffic data distribution characteristic value is calculated for the data stored in each barrelThree parameters, the calculation formula of which is defined as follows:
wherein Is->Barrel corresponding to each interval, probability distribution function->Indicating when the distribution model parameter value is +.>The time characteristic value is less than or equal to%>Probability of->Is an abbreviation for empirical cumulative distribution function (Empirical Cumulative Distribution Function) for a non-parametric tool describing the cumulative distribution of random variables. Which is calculated from the sample data. />Representing that the characteristic value is less than or equal to +.>The calculation formula is as follows:
wherein Indicate->Then it is 1.
Namely, the maximum forward K-S distance and the maximum reverse K-S distance obtained by K-S check calculation are +.>For storage at->And (3) rapidly counting the feature quantity in each barrel by using the tree-shaped array in the barrel in the step (3).
In the step 4, in the second step, all buckets in the bucket array are traversed and calculatedMaximum +.>The absolute value of (2) is taken as the K-S check fitting distance of the characteristic distribution model and is recorded as +.>The calculation formula is as follows:
step 5, calculating all characteristic distribution modelsValue, will->The characteristic distribution model with the minimum value is the data distribution model with the highest fitting degree of the normal safety state of the industrial control equipment, and is used as the reference distribution characteristic of the normal safety state of the industrial control equipment. At the same time +.>The value serves as an abnormal state perception threshold.
2. The industrial control equipment abnormal safety state sensing and monitoring stage comprises the following specific steps:
step 1, acquiring a network flow data characteristic value, namely a real-time network flow data characteristic value, in the running process of industrial control equipment;
step 2, inserting the newly acquired numerical values into the barrel-tree two-dimensional indexes to carry out quick sorting, simultaneously counting the number of the inserted numerical values in each barrel, and calculating an experience cumulative distribution function for the sorted data, wherein the specific calculation process is as described above;
step 3, adopting K-S test to calculate the coincidence degree of the probability distribution function corresponding to the standard distribution characteristic of the normal safety state of the industrial control equipment and the experience cumulative distribution function to obtain a check deviation value;
and step 4, judging that the industrial control equipment is in an abnormal safety state when the check deviation value is larger than the abnormal state sensing threshold value.
Another embodiment of the present application provides an industrial automation control device safety state monitoring apparatus, comprising a processor configured to:
before entering a safety state monitoring mode, collecting a network flow data characteristic value of industrial control equipment in a normal safety state, and dividing the network flow data characteristic value into a fitting sample and a checking sample;
selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples one by adopting a maximum likelihood estimation method to obtain corresponding fitting probability distribution functions;
rapidly sequencing the check samples by adopting a bucket-tree two-dimensional indexing method to obtain sequenced first data, and calculating an experience cumulative distribution function corresponding to the first data;
calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test;
determining a target data distribution model matched with the normal safety state of industrial control equipment, and determining an abnormal state perception threshold corresponding to the target data distribution model;
after entering a safety state monitoring mode, acquiring a real-time network flow data characteristic value of the industrial control equipment;
inserting the real-time network flow data characteristic values into a barrel-tree two-dimensional index for quick sequencing to obtain sequenced second data, and calculating a real-time experience cumulative distribution function corresponding to the second data;
adopting K-S test to calculate the coincidence degree of the probability distribution function and the real-time experience cumulative distribution function corresponding to the normal safety state of the industrial control equipment, and obtaining a check deviation value;
and when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state.
Further, in the above embodiment, the data distribution model includes at least a gaussian distribution model, a lognormal distribution model, a poisson distribution model, a weibull distribution model, a gamma distribution model, and an exponential distribution model.
Further, in the above embodiment, the selecting a plurality of data distribution models, and performing fitting analysis on the fitting samples by using a maximum likelihood estimation method one by one to obtain a corresponding fitting probability distribution function, specifically includes:
selecting a data distribution model from a plurality of data distribution models;
calculating likelihood functions corresponding to the data distribution model, wherein the likelihood functions are as follows:
wherein L (theta) is a likelihood function, X is a network flow data characteristic value set in acquisition time, and theta is a model parameter of a selected data distribution model;
taking the logarithm of the likelihood function to obtain a log likelihood function as follows:
deriving the log-likelihood function to obtain a corresponding maximum likelihood estimation parameter value;
and carrying the maximum likelihood estimation parameter value into an accumulated distribution function corresponding to the selected characteristic distribution model for calculation, wherein the calculation formula is as follows:
wherein ,for accumulating distribution functions +.>Estimating parameter values for the maximum likelihood;
repeating the steps until the maximum likelihood estimation parameter values and the accumulated distribution function corresponding to all the data distribution models are obtained.
Further, in the foregoing embodiment, the fast sorting of the check samples by using a bucket-tree two-dimensional indexing method specifically includes:
distinguishing the data in the check sample according to a barrel dividing mode, and storing the data in a corresponding barrel array;
and arranging the data in each bucket array according to a tree array structure.
Further, in the above embodiment, the calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to the K-S test specifically includes:
according to K-S test, respectively calculating a first parameter, a second parameter and a third parameter of a network flow data distribution characteristic value corresponding to the data stored in each barrel group, wherein the calculation formula is as follows:
wherein ,to reflect the first parameter of maximum forward K-S distance,/>To reflect the second parameter of the maximum inverse K-S distance,/a second parameter of the maximum inverse K-S distance is chosen>To reflect the third parameter of the feature quantity stored in the j-th bucket array, +.>Is->Barrel array corresponding to each interval, < > and the like>For probability distribution function +.>For the empirical cumulative distribution function, the calculation formula is as follows:
wherein ,indicate->Then 1;
and obtaining the K-S check fitting distance reflecting the fitting degree based on the maximum first parameter and the maximum second parameter.
The method and the device for monitoring the safety state of the industrial automation control equipment have the beneficial effects that at least one point of the following is:
before entering a safety state monitoring mode, fully considering the matching degree of each model in the normal safety state of the industrial control equipment, screening and determining a target data distribution model matched with the normal safety state of the industrial control equipment, determining an abnormal state sensing threshold corresponding to the target data distribution model, then entering a real-time monitoring mode, monitoring the safety state of the industrial automation control equipment, continuously calculating the network flow data characteristic distribution in the running process of the equipment, adopting K-S test to judge the coincidence degree of the network flow data characteristic distribution and the reference characteristic distribution, and if the deviation is larger than the abnormal state sensing threshold, judging that the running state is abnormal. The whole monitoring method does not need to invade the equipment to internally install a specific program to acquire the running state of the system for safety state sensing, has low operation complexity in the whole process, and is suitable for being deployed in the industrial environment site with limited computing resources for effective safety state monitoring of industrial control equipment.
The foregoing examples illustrate only a few embodiments of the application and are described in detail herein without thereby limiting the scope of the application. It should be noted that it will be apparent to those skilled in the art that several variations and modifications can be made without departing from the spirit of the application, which are all within the scope of the application. Accordingly, the scope of protection of the present application is to be determined by the appended claims.
Claims (8)
1. A method for monitoring the safety state of industrial automation control equipment, comprising the steps of:
before entering a safety state monitoring mode, collecting a network flow data characteristic value of industrial control equipment in a normal safety state, and dividing the network flow data characteristic value into a fitting sample and a checking sample;
selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples one by adopting a maximum likelihood estimation method to obtain corresponding fitting probability distribution functions;
rapidly sequencing the check samples by adopting a bucket-tree two-dimensional indexing method to obtain sequenced first data, and calculating an experience cumulative distribution function corresponding to the first data;
calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test;
determining a target data distribution model matched with the normal safety state of industrial control equipment, and determining an abnormal state perception threshold corresponding to the target data distribution model;
after entering a safety state monitoring mode, acquiring a real-time network flow data characteristic value of the industrial control equipment;
inserting the real-time network flow data characteristic values into a barrel-tree two-dimensional index for quick sequencing to obtain sequenced second data, and calculating a real-time experience cumulative distribution function corresponding to the second data;
adopting K-S test to calculate the coincidence degree of the probability distribution function and the real-time experience cumulative distribution function corresponding to the normal safety state of the industrial control equipment, and obtaining a check deviation value;
when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state;
the method comprises the steps of calculating fitting degrees of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test, and specifically comprises the following steps:
according to K-S test, respectively calculating a first parameter, a second parameter and a third parameter of a network flow data distribution characteristic value corresponding to the data stored in each bucket group, wherein the calculation formula is as follows:
wherein ,to reflect the first parameter of maximum forward K-S distance,/>To reflect the second parameter of the maximum inverse K-S distance,/a second parameter of the maximum inverse K-S distance is chosen>To reflect the third parameter of the feature quantity stored in the j-th bucket array, +.>Is->Barrel array corresponding to each interval, < > and the like>For accumulating distribution functions +.>For the empirical cumulative distribution function, the calculation formula is as follows:
wherein ,indicate->Then 1;
obtaining a K-S check fitting distance reflecting the fitting degree based on the largest first parameter and the largest second parameter;
traversing all buckets in the bucket array, computing themMaximum +.>The K-S check fitting distance of the absolute value of (2) as the characteristic distribution model is recorded as +.>The calculation formula is as follows:
computing all feature distribution modelsValue, will->The characteristic distribution model with the minimum value is the data distribution model with the highest fitting degree of the normal safety state of the industrial control equipment, is used as the standard distribution characteristic of the normal safety state of the industrial control equipment, and simultaneously is the +.>The value serves as an abnormal state perception threshold.
2. The method of claim 1, wherein the data distribution model includes at least a gaussian distribution model, a lognormal distribution model, a poisson distribution model, a weibull distribution model, a gamma distribution model, and an exponential distribution model.
3. The method for monitoring the safety state of industrial automation control equipment according to claim 2, wherein the selecting a plurality of data distribution models, and performing fitting analysis on the fitting samples by adopting a maximum likelihood estimation method one by one to obtain corresponding fitting probability distribution functions, specifically comprises:
selecting a data distribution model from a plurality of data distribution models;
calculating likelihood functions corresponding to the data distribution model, wherein the likelihood functions are as follows:
wherein L (theta) is a likelihood function, X is a network flow data characteristic value set in acquisition time, and theta is a model parameter of a selected data distribution model;
taking the logarithm of the likelihood function to obtain a log likelihood function as follows:
deriving the log-likelihood function to obtain a corresponding maximum likelihood estimation parameter value;
and carrying the maximum likelihood estimation parameter value into an accumulated distribution function corresponding to the selected characteristic distribution model for calculation, wherein the calculation formula is as follows:
wherein ,for accumulating distribution functions +.>Estimating parameter values for the maximum likelihood;
repeating the steps until the maximum likelihood estimation parameter values and the accumulated distribution function corresponding to all the data distribution models are obtained.
4. The method for monitoring the safety state of industrial automation control equipment according to claim 3, wherein the fast sorting of the check samples by adopting a bucket-tree two-dimensional indexing method specifically comprises:
distinguishing the data in the check sample according to a barrel dividing mode, and storing the data in a corresponding barrel array;
and arranging the data in each bucket array according to a tree array structure.
5. An industrial automation control device safety state monitoring apparatus comprising a processor configured to:
before entering a safety state monitoring mode, collecting a network flow data characteristic value of industrial control equipment in a normal safety state, and dividing the network flow data characteristic value into a fitting sample and a checking sample;
selecting a plurality of data distribution models, and carrying out fitting analysis on the fitting samples one by adopting a maximum likelihood estimation method to obtain corresponding fitting probability distribution functions;
rapidly sequencing the check samples by adopting a bucket-tree two-dimensional indexing method to obtain sequenced first data, and calculating an experience cumulative distribution function corresponding to the first data;
calculating the fitting degree of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test;
determining a target data distribution model matched with the normal safety state of industrial control equipment, and determining an abnormal state perception threshold corresponding to the target data distribution model;
after entering a safety state monitoring mode, acquiring a real-time network flow data characteristic value of the industrial control equipment;
inserting the real-time network flow data characteristic values into a barrel-tree two-dimensional index for quick sequencing to obtain sequenced second data, and calculating a real-time experience cumulative distribution function corresponding to the second data;
adopting K-S test to calculate the coincidence degree of the probability distribution function and the real-time experience cumulative distribution function corresponding to the normal safety state of the industrial control equipment, and obtaining a check deviation value;
when the check deviation value is larger than the abnormal state sensing threshold value, judging that the industrial control equipment is in an abnormal safety state;
the method comprises the steps of calculating fitting degrees of each fitting probability distribution function and the corresponding experience cumulative distribution function according to K-S test, and specifically comprises the following steps:
according to K-S test, respectively calculating a first parameter, a second parameter and a third parameter of a network flow data distribution characteristic value corresponding to the data stored in each bucket group, wherein the calculation formula is as follows:
wherein ,to reflect the first parameter of maximum forward K-S distance,/>To reflect the second parameter of the maximum inverse K-S distance,/a second parameter of the maximum inverse K-S distance is chosen>To reflect the third parameter of the feature quantity stored in the j-th bucket array, +.>Is->Barrel array corresponding to each interval, < > and the like>For accumulating distribution functions +.>For the empirical cumulative distribution function, the calculation formula is as follows:
wherein ,indicate->Then 1;
obtaining a K-S check fitting distance reflecting the fitting degree based on the largest first parameter and the largest second parameter;
traversing all buckets in the bucket array, computing themMaximum +.>The K-S check fitting distance of the absolute value of (2) as the characteristic distribution model is recorded as +.>The calculation formula is as follows:
computing all feature distribution modelsValue, will->The characteristic distribution model with the minimum value is the data distribution model with the highest fitting degree of the normal safety state of the industrial control equipment, is used as the standard distribution characteristic of the normal safety state of the industrial control equipment, and simultaneously is the +.>The value serves as an abnormal state perception threshold.
6. The industrial automation control device safety state monitoring device of claim 5, wherein the data distribution model includes at least a gaussian distribution model, a lognormal distribution model, a poisson distribution model, a weibull distribution model, a gamma distribution model, and an exponential distribution model.
7. The apparatus for monitoring the safety state of an industrial automation control device according to claim 6, wherein the selecting the plurality of data distribution models, performing fitting analysis on the fitting samples by using a maximum likelihood estimation method one by one, and obtaining a corresponding fitting probability distribution function, specifically includes:
selecting a data distribution model from a plurality of data distribution models;
calculating likelihood functions corresponding to the data distribution model, wherein the likelihood functions are as follows:
wherein L (theta) is a likelihood function, X is a network flow data characteristic value set in acquisition time, and theta is a model parameter of a selected data distribution model;
taking the logarithm of the likelihood function to obtain a log likelihood function as follows:
deriving the log-likelihood function to obtain a corresponding maximum likelihood estimation parameter value;
and carrying the maximum likelihood estimation parameter value into an accumulated distribution function corresponding to the selected characteristic distribution model for calculation, wherein the calculation formula is as follows:
wherein ,for accumulating distribution functions +.>Estimating parameter values for the maximum likelihood;
repeating the steps until the maximum likelihood estimation parameter values and the accumulated distribution function corresponding to all the data distribution models are obtained.
8. The industrial automation control device safety state monitoring apparatus of claim 7, wherein the fast sorting of the check samples by bucket-tree two-dimensional indexing method comprises:
distinguishing the data in the check sample according to a barrel dividing mode, and storing the data in a corresponding barrel array;
and arranging the data in each bucket array according to a tree array structure.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311021790.7A CN116736781B (en) | 2023-08-15 | 2023-08-15 | Safety state monitoring method and device for industrial automation control equipment |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202311021790.7A CN116736781B (en) | 2023-08-15 | 2023-08-15 | Safety state monitoring method and device for industrial automation control equipment |
Publications (2)
Publication Number | Publication Date |
---|---|
CN116736781A CN116736781A (en) | 2023-09-12 |
CN116736781B true CN116736781B (en) | 2023-11-03 |
Family
ID=87911847
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202311021790.7A Active CN116736781B (en) | 2023-08-15 | 2023-08-15 | Safety state monitoring method and device for industrial automation control equipment |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116736781B (en) |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107167781A (en) * | 2017-03-31 | 2017-09-15 | 西安电子科技大学 | The quantile method of estimation of sea clutter amplitude lognormal distribution parameter |
CN108804806A (en) * | 2018-06-05 | 2018-11-13 | 西南交通大学 | Weibull is distributed the simplification MLE methods of parameter in combined stress CA model |
CN110889554A (en) * | 2019-11-27 | 2020-03-17 | 东南大学 | Power load fluctuation analysis and risk early warning method based on recurrence time interval analysis method |
CN113504515A (en) * | 2021-06-28 | 2021-10-15 | 中国人民解放军海军航空大学航空作战勤务学院 | Method and device for estimating parameters and forming detection threshold of echo extreme value model |
CN113919216A (en) * | 2021-10-08 | 2022-01-11 | 北京航空航天大学 | Parameter uncertainty quantitative measurement method under small sub-sample condition |
CN116167240A (en) * | 2023-03-07 | 2023-05-26 | 四川大学 | Multi-measuring-point combined monitoring method for dam structure damage |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107895058B (en) * | 2017-07-19 | 2018-08-31 | 厦门理工学院 | A kind of method of quick identification wind speed Optimal Distribution rule |
-
2023
- 2023-08-15 CN CN202311021790.7A patent/CN116736781B/en active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107167781A (en) * | 2017-03-31 | 2017-09-15 | 西安电子科技大学 | The quantile method of estimation of sea clutter amplitude lognormal distribution parameter |
CN108804806A (en) * | 2018-06-05 | 2018-11-13 | 西南交通大学 | Weibull is distributed the simplification MLE methods of parameter in combined stress CA model |
CN110889554A (en) * | 2019-11-27 | 2020-03-17 | 东南大学 | Power load fluctuation analysis and risk early warning method based on recurrence time interval analysis method |
CN113504515A (en) * | 2021-06-28 | 2021-10-15 | 中国人民解放军海军航空大学航空作战勤务学院 | Method and device for estimating parameters and forming detection threshold of echo extreme value model |
CN113919216A (en) * | 2021-10-08 | 2022-01-11 | 北京航空航天大学 | Parameter uncertainty quantitative measurement method under small sub-sample condition |
CN116167240A (en) * | 2023-03-07 | 2023-05-26 | 四川大学 | Multi-measuring-point combined monitoring method for dam structure damage |
Non-Patent Citations (2)
Title |
---|
基于Matlab的某电厂给煤机可靠性分析;宋雨佳;徐有宁;;《沈阳工程学院学报》;第16卷(第3期);第39-44页 * |
基于网络流量异常检测的电网工控系统安全监测技术;钟志琛;;电力信息与通信技术(01);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN116736781A (en) | 2023-09-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107241226B (en) | Fuzzy test method based on industrial control private protocol | |
CN111935170B (en) | Network abnormal flow detection method, device and equipment | |
CN111654489B (en) | Network security situation sensing method, device, equipment and storage medium | |
CN102340485B (en) | Network security situation awareness system and method based on information correlation | |
CN107579986B (en) | Network security detection method in complex network | |
CN110896386B (en) | Method, device, storage medium, processor and terminal for identifying security threat | |
CN113408609A (en) | Network attack detection method and system | |
CN116366374B (en) | Security assessment method, system and medium for power grid network management based on big data | |
CN117395076B (en) | Network perception abnormality detection system and method based on big data | |
CN109474603A (en) | Data packet capturing processing method and terminal device | |
CN112749097B (en) | Performance evaluation method and device for fuzzy test tool | |
KR101281456B1 (en) | Apparatus and method for anomaly detection in SCADA network using self-similarity | |
CN115795330A (en) | Medical information anomaly detection method and system based on AI algorithm | |
CN114598506B (en) | Industrial control network security risk tracing method and device, electronic equipment and storage medium | |
CN117857225B (en) | Identity authentication system and method for new energy power station acquisition terminal | |
Sen et al. | Towards an approach to contextual detection of multi-stage cyber attacks in smart grids | |
CN117834311A (en) | Malicious behavior identification system for network security | |
CN116736781B (en) | Safety state monitoring method and device for industrial automation control equipment | |
Yan et al. | Detect and identify DDoS attacks from flash crowd based on self-similarity and Renyi entropy | |
CN115766081A (en) | Abnormal flow detection method and device for power industrial control cloud platform | |
Barsha et al. | Anomaly Detection in SCADA Systems: A State Transition Modeling | |
CN113132301B (en) | Abnormal data collection detection method and device and computer storage medium | |
KR102471618B1 (en) | Netflow based large-scale service network aceess tracking method and device and system therefor | |
CN116527378B (en) | Cloud mobile phone monitoring management method and system | |
CN114448699B (en) | Data detection method, device, electronic equipment and storage medium |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |