Non-stationary process early fault diagnosis method based on common trend model

Abstract

The invention discloses a non-stationary industrial process early fault diagnosis method based on a common trend model, and belongs to the field of industrial process monitoring and fault diagnosis. The method comprises the following steps: collecting training data of a non-stationary industrial process under a normal working condition, performing denoising by using singular spectrum analysis, and obtaining a stationary projection matrix and a stationary data matrix by using a common trend model; orthogonal conversion is carried out on the stationary data matrix, a conversion projection matrix is calculated, and an orthogonal conversion data matrix and a statistic data matrix thereof are obtained; calculating a mean value and a covariance matrix of the statistic data matrix to obtain a control limit; the method comprises the following steps: acquiring test data under a real-time working condition of an industrial process, processing to obtain statistic data, calculating a mahalanobis distance index, and comparing with a control limit to judge whether a fault occurs; and if a fault is detected, realizing fault separation according to fault reconstruction and a mapping relation. According to the method, the interference of the non-stationary characteristic on the fault diagnosis model can be reduced, and fault detection and separation in the non-stationary process are realized.

Description

Technical Field

[0001] The present invention belongs to the field of industrial process monitoring and fault diagnosis, and particularly relates to an early fault diagnosis method for non-stationary industrial processes based on a common trend model. Background Art

[0002] With the continuous advancement of industrial automation and intelligent technologies, the complexity of industrial processes has been continuously increasing, and their non-stationary characteristics have gradually emerged and taken the dominant position. This non-stationarity is usually jointly caused by various factors such as equipment aging, production plan adjustment, external environmental changes, and human intervention, resulting in significant changes in the statistical characteristics (such as mean and variance) of process variables over time. However, most current multivariate statistical process monitoring methods are based on the assumption of stationary data distribution and are difficult to effectively handle the complex statistical feature changes in non-stationary processes. This limitation often leads to high false alarm rates and missed alarm rates, significantly reducing the performance of fault diagnosis methods. In actual industrial scenarios, the deficiencies of traditional fault diagnosis methods have gradually become apparent. Therefore, developing a fault diagnosis method that can adapt to non-stationary characteristics is not only a key direction in the current industrial monitoring field but also an important guarantee for the reliable and safe operation of industrial systems.

[0003] In industrial processes, faults usually develop gradually from the early stage, initially showing small amplitudes and unclear characteristics. These early faults have limited direct impact on system operation. However, if not detected and handled in a timely manner, they may gradually evolve into serious faults over time, posing a major threat to system safety and production efficiency. Therefore, industrial monitoring systems need to be highly sensitive to early faults and be able to quickly identify and take effective measures at the budding stage of faults, thereby avoiding fault escalation. In recent years, certain progress has been made in data-driven methods for early fault diagnosis. However, in non-stationary operating conditions full of noise, the weak abnormal signals of early faults are easily masked by non-stationary trends, further increasing the detection difficulty. In addition, research on early fault separation is still relatively scarce. This current situation highlights the importance and urgency of developing an early fault diagnosis method applicable to non-stationary processes.

[0004] Therefore, there is an urgent need for a new fault diagnosis method to achieve the detection and separation of early faults in non-stationary industrial processes. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide an early fault diagnosis method applicable to non-stationary industrial processes. To solve the above technical problem, the present invention provides an early fault diagnosis method based on a common trend model, and the method includes the following steps:

[0006] Step 1: Collect a segment of measurement data under normal conditions of the non-stationary industrial process as the training dataset, perform denoising processing using singular spectrum analysis, and calculate the stationary projection matrix using the common trend model to obtain the stationary data matrix;

[0007] Step 2: Perform orthogonal transformation processing on the stationary data matrix obtained in Step 1, calculate the transformation projection matrix to obtain the orthogonally transformed data matrix. Given an appropriate sliding time window width, extract the statistics of the orthogonally transformed data matrix within each window to obtain the statistic data matrix;

[0008] Step 3: Calculate the mean vector and covariance matrix of the statistic data matrix in Step 2, calculate the Mahalanobis distance index, and use kernel density estimation to calculate the control limit of the Mahalanobis distance index;

[0009] Step 4: Collect a segment of measurement data under real-time conditions of the industrial process as the test data, process the test data using the denoising parameters in Step 1, and perform stationary processing on the denoised test data using the stationary projection matrix to obtain the stationary test data. Then, perform orthogonal transformation processing on the stationary test data using the transformation projection matrix in Step 3 to obtain the orthogonally transformed test data, and further calculate its statistics;

[0010] Step 5: Calculate the Mahalanobis distance index of the statistics of the orthogonally transformed test data in Step 4 using the mean and covariance parameters in Step 3, and compare it with the control limit in Step 3 to determine whether a fault has occurred, realizing fault detection;

[0011] Step 6: If a fault is detected in Step 5, calculate the reconstruction contribution value of the corresponding statistic of the orthogonally transformed test variable to the fault detection index through the fault reconstruction method, so as to determine the abnormal orthogonally transformed test variable, and determine the fault variable according to the mapping relationship between it and the original measurement variable, realizing fault isolation.

[0012] Preferably, Step 1 is specifically:

[0013] Record the collected fault-free training dataset as matrix where represents a sample containing m measurements, and there are N independent samples in total; perform denoising processing on the training dataset using singular spectrum analysis. For each measurement variable {x i (k); k = 1, 2,..., N, i = 1, 2,..., m}, it is all transformed into its trajectory matrix G:

[0014]

[0015] where L is the selected embedding window width, 1 < L < N, K = N - L + 1; perform singular value decomposition on the trajectory matrix G:

[0016] G = UΣV T

[0017] Where U and V are the left singular vector matrix and the right singular vector matrix respectively, and Σ is the singular value matrix; define the matrix:

[0018]

[0019] Perform eigenvalue decomposition on it to obtain L non-zero eigenvalues λ 1 ≥ λ 2 ≥ … ≥ λ L > 0, and the corresponding eigenvectors are v 1 , v 2 ,..., v L ; Denote d = rank(G), and the left singular vector of the trajectory matrix is denoted as Then G is re-expressed as:

[0020]

[0021] Select the first r components with large contributions according to the magnitudes of the eigenvalues to reconstruct the original measured variable to complete the denoising process, where r is calculated by the cumulative variance contribution rate; 1 ≤ r ≤ d; Using the diagonal averaging method, let the elements of G be g(k - i + 1, i), L * = min(L, K), K * = max(L, K), and restore the matrix G to the form of the original measured variable:

[0022]

[0023] Each denoised measured variable is denoted as They together form a new denoised dataset which is the denoised sample at time k;

[0024] For the denoised dataset The common trend model is expressed as:

[0025]

[0026] Where X ns is non-stationary data, and X s is stationary data; To identify the stationary data of the common trend model, Kasa decomposition is adopted:

[0027]

[0028] Where the stationary projection matrix is a linear combination of the cointegration matrix B, denoted as The cointegration matrix B is obtained by solving the following eigenvalue equation:

[0029]

[0030] Among them, e 0 and e 1 are obtained by the following method:

[0031]

[0032] Among them, the coefficient matrix Θ i and Φ i are obtained by least squares estimation, where p is the model lag order; finally, the stationary data matrix

[0033] Preferably, step 2 is specifically as follows:

[0034] Perform orthogonal decomposition on the stationary data matrix to obtain the eigenvector matrix P s , and the orthogonally transformed data matrix is obtained by projection:

[0035] T = X s W t = X s P s

[0036]

[0037] Among them, the transformation projection matrix is denoted as Finally, the orthogonally transformed data matrix Given an appropriate sliding time window width w, extract statistics from the orthogonally transformed data matrix within each window. The sample statistics of the measurement variables within the current window have the following vector form:

[0038]

[0039] Among them, μ 1 , μ 2 represents the sample mean of the measurement variables x 1 , x 2 , represents the sample variance of the measurement variables x 1 , x 2 , max 1 , max 2 represents the sample maximum of the measurement variables x 1 , x 2 , min 1 , min 2 represents the sample minimum of the measurement variables x 1 , x 2The sample minimum value, and so on; k represents the statistic vector calculated using the data within a window of width w with the current moment as the base point; for the training data set, by continuously moving the sliding time window, each window obtains a vector, and the vectors obtained from all windows are arranged row by row to construct a statistic data matrix.

[0040] Preferably, step 3 is specifically as follows:

[0041] For the statistic data matrix S constructed in step 2, its mean vector is denoted as with m s elements, and each element is obtained by calculating the mean of the corresponding column of the statistic data matrix; its covariance matrix is denoted as C s , which is calculated from the statistic data matrix with the mean zeroed; the control limit of the Mahalanobis distance index is estimated by kernel density estimation and is denoted as σ s ; α is the significance level, taking the value of 0.01;

[0042] The Mahalanobis distance index is calculated according to the following formula:

[0043]

[0044] The control limit σ s is calculated according to the following formula:

[0045]

[0046] where P(·) is its probability distribution function, and its probability density function is estimated by kernel density estimation:

[0047]

[0048] where h is the bandwidth, and the kernel function K(·) selects the Gaussian kernel function, that is, K(·) = φ(·), and φ(·) is the standard Gaussian density function.

[0049] Preferably, step 4 is specifically as follows:

[0050] Denote the sample of the test data at the current moment in the test stage as x(k'), and use the denoising parameter in step 1 to process the test data to obtain the denoised test data and use the stationary projection matrix W s to perform stationary processing on the denoised test data to obtain the stationary test data:

[0051]

[0052] Use the transformation projection matrix W in step 2 t to perform orthogonal transformation processing on the stationary test data, thereby obtaining the orthogonally transformed test data:

[0053]

[0054] Using the same window parameter w as in step 2, calculate the statistic vector of the orthogonal transformation test data at the current moment as s(k′), which has a similar form to s(k) shown in step 2.

[0055] Preferably, step 5 is specifically as follows:

[0056] Subtract the mean s in step 3 from s(k′) in step 4, and calculate the Mahalanobis distance index according to the following formula:

[0057]

[0058] where C s is the statistic covariance matrix in step 3; compare the Mahalanobis distance index with the control limit σ s in step 3, that is, when D s (k′) > σ s , it is considered that a fault has occurred in the non-stationary industrial process.

[0059] Preferably, step 6 is specifically as follows:

[0060] When a fault is detected in step 5, calculate the statistic reconstruction contribution value of the i-th orthogonal transformation test variable according to the following formula:

[0061]

[0062] where Ξ i represents the fault direction corresponding to the i-th orthogonal transformation test variable in the statistic space when a fault occurs, and obtain the abnormal orthogonal transformation test variable corresponding to the maximum GC i ; according to obtain the mapping relationship between the orthogonal transformation test variable and the original measurement variable Based on the principle that the variable with a larger weight coefficient occupies the dominant position, determine the corresponding fault variable; if the i-th orthogonal transformation test variable is an abnormal orthogonal transformation test variable, the variable corresponding to the number with the largest absolute value in the i-th row of the mapping matrix W is the fault variable, thereby realizing fault separation.

[0063] The beneficial technical effects brought by the present invention:

[0064] The present invention provides an early fault diagnosis method for non-stationary industrial processes. By using the denoised data for modeling, the accuracy and robustness of the detection of non-stationary industrial processes are greatly improved; the common trend model is adopted to effectively separate the long-term trends in non-stationary variables, convert the denoised data into stationary data, significantly reduce the interference of non-stationary characteristics on the detection model, and thus effectively reduce the false alarm rate; by combining the method of extracting statistics after orthogonal transformation of the data, the data characteristics in the industrial process can be comprehensively captured with fewer statistical features, realizing the early and accurate identification of various types of minor faults, and thus effectively improving the detection rate. In addition, after fault detection, fault separation can be carried out in a timely manner, further enhancing the practicability and reliability of the method.

[0065] Other advantages, objectives, and features of the present invention will be described to some extent in the subsequent specification, and, in part, will become apparent from the specification, or will be understood by implementing the present invention. The objectives and other advantages of the present invention can be achieved and obtained through the structures specifically pointed out in the specification, claims, and drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] The drawings are used to provide a further understanding of the present invention, and constitute a part of the specification, and are used together with the embodiments of the present invention to explain the present invention, and do not constitute a limitation to the present invention. In the drawings:

[0067] Figure 1 is a schematic flowchart of an early fault diagnosis method for a non-stationary process based on a common trend model according to an embodiment of the present invention;

[0068] Figure 2 is a schematic flowchart of an offline modeling process according to an embodiment of the present invention;

[0069] Figure 3 is a schematic flowchart of an online diagnosis process according to an embodiment of the present invention;

[0070] Figure 4 is a schematic diagram of the detection result of a deviation fault based on the method proposed by the present invention according to an example of the present invention;

[0071] Figure 5 is a schematic diagram of the detection result of a precision degradation fault based on the method proposed by the present invention according to an example of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0072] The present invention will be further described in detail below in conjunction with the drawings and specific embodiments:

[0073] Figure 1It is a schematic diagram of an early fault diagnosis process for non-stationary processes based on a common trend model according to an embodiment of the present invention. This method first uses singular spectrum analysis to denoise the measurement data. Secondly, it uses the common trend model to obtain a stationary projection matrix and then obtains a stationary data matrix. On this basis, the stationary data matrix is subjected to orthogonal transformation processing, the transformation projection matrix is calculated, and then the orthogonally transformed data and its statistical data are obtained. Then, the Mahalanobis distance index is used for fault detection, and finally, the methods of fault reconstruction and linear mapping are used for fault separation. This method includes the following steps:

[0074] Step S110: Collect a section of measurement data under normal working conditions of the non-stationary industrial process as the training data set, use singular spectrum analysis for denoising processing, and use the common trend model to calculate the stationary projection matrix and then obtain the stationary data.

[0075] Step S120: Perform orthogonal transformation processing on the stationary data obtained in step S110, calculate the transformation projection matrix to obtain the orthogonally transformed data matrix, and given an appropriate sliding time window width, extract the statistics of the orthogonally transformed data matrix within each window to obtain the statistical data matrix.

[0076] Step S130: Calculate the mean vector and covariance matrix of the statistical data matrix in step S120, calculate the Mahalanobis distance index, and use kernel density estimation to calculate the control limit of the Mahalanobis distance index.

[0077] Step S140: Collect a section of measurement data under the real-time working conditions of the industrial process as the test data, use the denoising parameters in step S110 to process the test data, and use the stationary projection matrix to perform stationary processing on the denoised test data to obtain the stationary test data. Then, use the transformation projection matrix in step S120 to perform orthogonal transformation processing on the stationary test data to obtain the orthogonally transformed test data, and further calculate its statistics.

[0078] Step S150: Use the mean and covariance parameters in step S130 to calculate the Mahalanobis distance index of the statistics of the orthogonally transformed test data in step S140, and compare it with the control limit in step S130 to determine whether a fault has occurred, and achieve fault detection.

[0079] Step S160: If a fault is detected in step S150, calculate the reconstruction contribution value of the statistics corresponding to the orthogonally transformed test variables to the fault detection index through the fault reconstruction method, so as to determine the abnormal orthogonally transformed test variables, and determine the fault variables according to the mapping relationship between them and the original measurement variables, and achieve fault separation. Among them, steps S110, S120, and S130 are the offline modeling processes, and steps S140, S150, and S160 are the online diagnosis processes.

[0080] Figure 2 It is a schematic diagram of the detailed step - by - step process of the offline modeling in the embodiment of the present invention, and the specific implementation steps are as described below.

[0081] (1) Step S210: Collect a segment of measurement data under normal operating conditions of the non - stationary industrial process as the training data set, store it as a two - dimensional matrix, denoise the training data set, and then calculate the stationary data matrix. Specifically, denote the collected training data set as the matrix where represents a sample containing m measurements, and there are N independent samples in total. Each measurement variable is denoted as {x i (k); k = 1, 2, …, N, i = 1, 2,..., m}. Use singular spectrum analysis to denoise each measurement variable, and each denoised measurement variable is denoted as They together constitute a new denoised data set is the denoised sample at time k. Obtain the stationary projection matrix and then obtain the stationary data matrix

[0082] (2) Step S220: Calculate the orthogonal transformation data matrix and its statistic data matrix. Specifically, perform orthogonal transformation on the stationary data obtained in step S210, calculate the transformation projection matrix and then obtain the orthogonal transformation data matrix Given an appropriate sliding time - window width w, extract statistics within each window for the orthogonal transformation data matrix to obtain the statistic data matrix

[0083] (3) Step S230: Calculate the mean vector and covariance matrix of the statistic data matrix constructed in step S220, and calculate the Mahalanobis distance index. Specifically, denote the mean vector of the constructed statistic matrix as with m s elements, and each element is obtained by taking the mean of the corresponding column of the statistic data matrix; its covariance matrix is denoted as which can be calculated from the statistic data matrix with the mean zeroed.

[0084] (4) Step S240: Calculate the control limit of the Mahalanobis distance index. Specifically, the control limit can be calculated by kernel density estimation, denoted as σ s . α is the significance level, usually taking the value of 0.01.

[0085] Figure 3 It is a schematic diagram of the detailed step - by - step process of the online detection in the embodiment of the present invention, and the specific implementation steps are as described below.

[0086] (1) Step S310: Collect a segment of measurement data under the real-time working conditions of the industrial process as test data, process the test data using the same denoising parameters as in step S210 of the offline modeling process, and perform a stationary processing on the denoised test data using the stationary projection matrix in the offline modeling process to obtain stationary test data. Specifically, similar to the offline modeling process, the denoised data matrix and sample vector are and Use the stationary projection matrix W in step S210 of the offline modeling process s Calculate the stationary test data:

[0087]

[0088] (2) Step S320: Perform an orthogonal transformation on the stationary test data to obtain the orthogonally transformed test data and its statistical data. Specifically, use the transformation projection matrix W in step S220 of the offline modeling process t Calculate the orthogonally transformed test data:

[0089]

[0090] Using the same window parameter w as in step S220 of the offline modeling process, calculate the statistical vector s(k′) of the orthogonally transformed test data at the current moment.

[0091] (3) Step S330: Calculate the Mahalanobis distance index of the statistics of the orthogonally transformed test data in step S320 using the mean and covariance parameters in the offline modeling process, and compare it with the corresponding control limit in step S240 of the offline modeling process to determine whether a fault has occurred. Specifically, subtract the mean s in step S230 of the offline modeling process from s(k′) in step S320 and calculate the Mahalanobis distance index according to the following formula:

[0092]

[0093] where Cs is the covariance matrix in step S230 of the offline modeling process. And compare D s (k′) with the control limit σ s in step S230 of the offline modeling process. If D s (k′) > σ s , it is considered that a fault has occurred in the non-stationary industrial process.

[0094] (6) Step S340: If a fault is analyzed in step S330, perform fault isolation. Specifically, when a fault is detected in step S330, calculate the statistical reconstruction contribution value of the i-th orthogonally transformed test variable according to the following formula:

[0095]

[0096] where Ξ i represents the fault direction corresponding to the occurrence of a fault in the statistic space when the i-th orthogonal transformation test variable fails. Obtain the maximum GC i corresponding abnormal orthogonal transformation test variable. According to obtain the mapping relationship between the orthogonal transformation test variable and the original measurement variable Based on the principle that the variable with a larger weight coefficient occupies a dominant position, determine the corresponding fault variable. That is, if the i-th orthogonal transformation test variable is an abnormal orthogonal transformation test variable, the variable corresponding to the number with the largest absolute value in the i-th row of the mapping matrix W is the fault variable, thereby realizing fault separation.

[0097] The early fault diagnosis method for non-stationary processes based on the common trend model in the embodiments of the present invention uses a large amount of normal data easily obtained from industrial processes for modeling and applies the model to online diagnosis. It does not require a complex mechanism model of the industrial process and does not require fault data that is difficult to obtain or label, and is easy to implement; it can simultaneously detect and separate sensor faults, providing valuable reference information for component repair and replacement; it is sensitive to incipient faults and can further improve the diagnostic performance for minor faults by selecting an appropriate window width.

[0098] The diagnostic method in the embodiments of the present invention uses singular spectrum analysis technology to denoise the data, then applies the common trend model to smooth the denoised data to obtain stationary data, and then performs orthogonal transformation on the stationary data to obtain orthogonal transformation data, further calculates its statistics, and finally constructs a Mahalanobis distance index for fault detection, and applies the idea of fault reconstruction and the mapping relationship to achieve fault location. Compared with traditional multivariate statistical analysis methods, it can detect early faults in non-stationary industrial processes and achieve fault separation.

[0099] Example

[0100] To help understand the present invention and intuitively show its effect on non-stationary industrial processes, an example is described below. This example is based on the Matlab tool and uses a numerical example in the existing literature (Yuanlin Lin, et al., Industrial & Engineering Chemistry Research, 2017, 56(31): 8895 - 8905) to describe the present invention and shows the effect of the present invention in combination with the accompanying drawings.

[0101] (1) Generate training data.

[0102] This example uses the following equation to generate N = 7000 normal samples:

[0103]

[0104] where

[0105]

[0106] There are 5 measurement variables in the simulation, \(x = [x 1 , x 2 , …, x 5 \) T , that is, \(m = 5\); \(e 1 \) and \(e 2 \) terms represent zero-mean Gaussian white noise, and the variances of the 5 components are all 0.25. All samples are stored as a two-dimensional data matrix in the way that each row represents a sample and each column represents a variable.

[0107] (2) Given an appropriate denoising parameter \(L\), calculate the denoised data.

[0108] In this example, the denoising parameter \(L\) is selected as \(L = 200\), and denoising processing is performed on each measurement variable \(\{x i (k); k = 1, 2, …, 7000, i = 1, 2, …, 5\}\) to obtain all denoised measurement variables They together constitute a new denoised data set

[0109] (3) Perform stationarization processing on the denoised data to obtain a stationary projection matrix and a stationary data matrix.

[0110] Perform stationarization processing on the denoised data, and calculate the cointegration vector matrix Stationary projection matrix and stationary data matrix

[0111] (4) Perform orthogonal transformation processing on the stationary data matrix to obtain an orthogonal transformation data matrix and its statistical quantity data matrix.

[0112] Perform orthogonal decomposition on the stationary data matrix to obtain a transformation projection matrix Furthermore, obtain the orthogonal transformation data matrix Given a sliding time window width \(w = 200\), extract statistical quantities of the orthogonal transformation data within each window. The statistical quantities of the measurement variables within the window at the current moment have the following vector form:

[0113]

[0114] A vector such as \(s(k)\) can be obtained for each window. The vectors obtained from all windows are arranged by rows to construct a statistical quantity data matrix

[0115] (5) Calculate the Mahalanobis distance index and give the control limits.

[0116] Calculate the mean vector and covariance matrix of the statistical data matrix, and give the control limits of the Mahalanobis distance index. Here, the kernel density estimation is used to calculate the control limit σ s , specify the significance level as α = 0.01, and the control limit of the Mahalanobis distance index is 293.56.

[0117] (6) Construct test data.

[0118] In this example, two types of faults are considered, namely: deviation fault, expressed as y = y * + f; precision degradation fault, expressed as y = y * + e. Based on the above expressions, two sets of test data are generated respectively, and each set of data contains 3000 samples. The fault is applied starting from the 601st sample and continues until the 3000th sample. The first type of fault is applied to the first measurement variable with an amplitude of 0.25; the second type of fault is applied to the second measurement variable, where the standard deviation of the zero-mean white noise e is 0.6.

[0119] (7) Perform fault detection on the test data.

[0120] Figure 4 and Figure 5 respectively show the fault detection results of the method proposed in the present invention for the above two types of faults. It can be seen that after the fault occurs, the statistical Mahalanobis distance index D s both increase significantly and exceed the control limit. For the two types of faults, the false alarm rate (FAR) of this detection index is 0.00% for both, and both fluctuate near the significance level of 1%, which belongs to a reasonable range. The fault detection rate (FDR) of this detection index is 98.12% and 99.00% respectively, and the results are very ideal.

[0121] In addition, a comparison is also made with the traditional method. Here, the schematic diagram of the detection results is not given, and only the FAR and FDR evaluation indexes of the method are listed. For the traditional principal component analysis PCA method, the FAR of its T 2 index is 0.00% for both, and the FDR is 0.00% for both. However, the FAR of the Q index is 90.18% and 99.83% respectively, and the FDR is 84.12% and 70.09% respectively. It can be seen from the Q index that its FAR is very high, indicating that it fails for non-stationary processes. For the traditional canonical variate analysis CVA method, its T 2The FARs of the indicators are 66.06% and 95.84% respectively, the FDRs are 75.91% and 1.00% respectively, the FARs of the Q indicators are 0.50% and 72.11% respectively, the FDRs are 20.55% and 79.37% respectively. Through the T 2 indicator, it can be seen that its FAR is still very high, indicating that it also fails for non-stationary processes. For the cointegration analysis CA method, its FAR is within a reasonable range (not listed here), and the FDRs are 16.88% and 51.98% respectively. From the FDR, it can be seen that it cannot effectively or well detect faults. For the time series model VAR / VARI method, its T s 2 and the T n 2 s indicators' FARs are all within a reasonable range (not listed here), and its T s 2 indicator's FDRs are 6.79% and 58.82% respectively, and its T n 2 s indicator's FDRs are 1.63% and 30.39% respectively. From the FDR, it can be seen that it cannot effectively or well detect faults. Through comparison, it can be found that the method proposed in the present invention can well detect two types of non-stationary process faults.

[0122] (6) If a fault occurs, perform fault isolation.

[0123] In terms of fault isolation, it is evaluated by the correct isolation rate (CIR) indicator. For the above two types of faults, the CIRs of the method proposed in the present invention are 93.33% and 90.08% respectively, and the effect is very satisfactory. The method proposed in the present invention can well isolate early faults.

[0124] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the substantial scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for early fault diagnosis of non-stationary processes based on a common trend model, characterized in that: The steps include: Step 1: Collect a section of measurement data under normal working conditions of a non-stationary industrial process as a training data set, use singular spectrum analysis to perform denoising, and use the common trend model to calculate the stationary projection matrix to obtain the stationary data matrix; Step 2: Perform an orthogonal transformation on the stationary data matrix obtained in step 1, calculate the transformation projection matrix and then obtain the orthogonal transformation data matrix. Given an appropriate sliding time window width, extract the statistics of the orthogonal transformation data matrix in each window to obtain the statistical data matrix; Step 3: Calculate the mean vector and covariance matrix of the statistical data matrix in step 2, calculate the Mahalanobis distance index, and use kernel density estimation to calculate the control limit of the Mahalanobis distance index; Step 4: Collect a section of measurement data under real-time working conditions of the industrial process as test data, process the test data using the denoising parameters in step 1, and use the stationary projection matrix to stabilize the denoised test data to obtain stationary test data, and then use the conversion projection matrix in step 3 to perform orthogonal transformation on the stationary test data to obtain orthogonal transformation test data, and further calculate its statistics; Step 5: Calculate the Mahalanobis distance index of the statistic of the orthogonal transformation test data in step 4 using the mean and covariance parameters in step 3, and compare it with the control limit in step 3 to determine whether a fault occurs, thereby realizing fault detection; Step 6: If a fault is detected in step 5, the reconstruction contribution value of the corresponding statistic of the orthogonal transformation test variable to the fault detection index is calculated through the fault reconstruction method, so as to determine the abnormal orthogonal transformation test variable, and determine the fault variable according to the mapping relationship between it and the original measurement variable to achieve fault separation.

2. The non-stationary process early fault diagnosis method based on the common trend model according to claim 1 is characterized in that: Step 1 is as follows: The collected fault-free training data set is recorded as a matrix in represents a sample containing m measurements, with a total of N independent samples; the training data set is denoised using singular spectrum analysis, and for each measurement variable {x i (k); k=1,2,…,N,i=1,2,…,m}, are all transformed into their trajectory matrix G: Where L is the selected embedding window width, 1<L<N, K=N-L+1; perform singular value decomposition on the trajectory matrix G: G=UΣV T Where U and V are the left singular vector matrix and the right singular vector matrix respectively, Σ is the singular value matrix; define the matrix: Perform eigenvalue decomposition and obtain L non-zero eigenvalues ​​λ1≥λ2≥…≥λ L > 0, the corresponding eigenvectors are v1,v2,…,v L ; Let d = rank(G), and the left singular vector of the trajectory matrix is Then G can be reformulated as: According to the size of the eigenvalue, the first r components with the largest contribution are selected to reconstruct the original measured variable to complete the denoising process. r is calculated by the cumulative variance contribution rate; 1≤r≤d; using the diagonal average method, let the element of G be g(k-i+1,i), L * =min(L,K),K * =max(L,K), and restore the matrix G to the form of the original measurement variables: Each denoised measurement variable is recorded as Together they form a new denoised dataset is the denoised sample at time k; For the denoised data set The common trend model is expressed as: Among them, X ns is non-stationary data, X s is stationary data; to identify the stationary data of the common trend model, Kasa decomposition is used: Among them, the stationary projection matrix is ​​a linear combination of the cointegration matrix B, denoted as The cointegration matrix B is obtained by solving the following characteristic equation: in, e0 and e1 are obtained by: in, Coefficient matrix Θ i and Φ i It is obtained through least squares estimation, where p is the lag order of the model; finally, the stable data matrix is ​​obtained 3. The non-stationary process early fault diagnosis method based on the common trend model according to claim 2 is characterized in that: Step 2 is as follows: Perform orthogonal decomposition on the stationary data matrix to obtain the eigenvector matrix P s , the orthogonal transformation data matrix is ​​obtained by projection: T=X s W t =X s P s Among them, the transformation projection matrix is ​​recorded as Finally, the orthogonal transformation data matrix is ​​obtained Given a suitable sliding time window width w, statistics are extracted from the orthogonal transformation data matrix in each window. The sample statistics of the measured variables in the window at the current moment have the following vector form: Among them, μ1, μ2 represent the sample means of the measured variables x1, x2, represents the sample variance of the measured variables x1 and x2, max1 and max2 represent the sample maximum values ​​of the measured variables x1 and x2, min1 and min2 represent the sample minimum values ​​of the measured variables x1 and x2, and so on; k represents the statistical vector calculated using the window data with a width of w based on the current moment; for the training data set, by continuously moving the sliding time window, each window obtains a vector, and the vectors obtained from all windows are arranged in rows to construct a statistical data matrix 4. The non-stationary process early fault diagnosis method based on the common trend model according to claim 3 is characterized in that: Step 3 is as follows: The statistical data matrix S constructed in step 2 has its mean vector recorded as With m s elements, each element is obtained by averaging the corresponding columns of the statistic data matrix; Its covariance matrix is ​​denoted as C s , calculated from the statistical data matrix after the mean is reset to zero; the control limit of the Mahalanobis distance index is estimated by kernel density estimation, denoted by σ s ; α is the significance level, which is 0.01; The Mahalanobis distance index is calculated as follows: Control limit σ s Calculate according to the following formula: Among them, P(·) is its probability distribution function, and its probability density function is estimated using kernel density estimation: Wherein, h is the bandwidth, and the kernel function K(·) selects the Gaussian kernel function, that is, K(·) = φ(·), where φ(·) is the standard Gaussian density function.

5. The non-stationary process early fault diagnosis method based on common trend model according to claim 4 is characterized in that: Step 4 is as follows: The sample of the test data at the current moment in the test phase is recorded as x(k'). The denoising parameters in step 1 are used to process the test data to obtain the denoised test data. And using the stationary projection matrix W s The denoised test data is stabilized to obtain stable test data: Using the transformation projection matrix W in step 2 t Perform orthogonal transformation on the stationary test data to obtain orthogonal transformation test data: Using the same window parameter w as in step 2, the statistic vector of the orthogonal transformation test data at the current moment is calculated as s(k′), which has a similar form to s(k) shown in step 2.

6. The non-stationary process early fault diagnosis method based on the common trend model according to claim 5 is characterized in that: Step 5 is as follows: Subtract the mean value from step 3 from s(k′) in step 4 The Mahalanobis distance index is calculated as follows: Among them, C s This is the statistical covariance matrix in step 3; the Mahalanobis distance index and the control limit σ in step 3 s In contrast, when D s (k′)>σ s , then it is considered that a failure has occurred in the non-stationary industrial process.

7. The non-stationary process early fault diagnosis method based on common trend model according to claim 6 is characterized in that: Step 6 is as follows: When a fault is detected in step 5, the statistical reconstruction contribution value of the i-th orthogonal transformation test variable is calculated as follows: Among them, i represents the fault direction corresponding to the fault of the i-th orthogonal transformation test variable in the statistical space when the fault occurs, and the maximum GC is obtained. i The corresponding abnormal orthogonal transformation test variable; according to Get the mapping relationship between the orthogonal transformation test variable and the original measurement variable Based on the principle that variables with larger weight coefficients occupy the dominant position, the corresponding fault variable is determined; if the i-th orthogonal transformation test variable is an abnormal orthogonal transformation test variable, the variable corresponding to the number with the largest absolute value in the i-th row of the mapping matrix W is the fault variable, thereby realizing fault separation.

Images