A method for monitoring the stuck drill status of mining drill based on characteristic stationarity test and cointegration analysis

Abstract

The present application relates to a method for monitoring the stuck drill state of a mining drill based on feature stationarity test and cointegration analysis. The method comprises: obtaining stuck drill state monitoring data, including vibration acceleration data, drilling pressure data, drilling displacement data and rotation torque data, extracting time domain features and time-frequency domain features of the vibration acceleration data, taking the drilling pressure data, drilling displacement data and rotation torque data as statistical variables, respectively performing stationarity tests on the time domain features, time-frequency domain features and statistical variables, dividing corresponding stationary signals and non-stationary signals, performing principal component analysis on the corresponding stationary signals, obtaining corresponding principal component analysis residuals, performing cointegration analysis on the non-stationary signals, obtaining corresponding cointegration analysis residuals, performing principal component analysis based on the principal component analysis residuals and the cointegration analysis residuals, obtaining corresponding T 2 Statistics are used to monitor the time when the drill pipe is stuck, which improves the accuracy and reliability of identifying the stuck state of the drill pipe.

Description

Technical Field

[0001] The present application relates to the technical field of mechanical equipment, and in particular to a method for monitoring the stuck drill state of a mining drill based on characteristic stationarity test and cointegration analysis. Background Art

[0002] As one of the important components to ensure the normal operation of the drilling rig, the stability of the drill pipe's operating state has the most direct impact on the safety and production of the equipment. When the drill pipe is affected by external coal seams or excitation, the operating state of the drill pipe will be greatly affected. In severe cases, it may even get stuck, resulting in downtime of the drilling rig, causing great economic losses, and even personal injury. Therefore, in order to further reduce the occurrence of drill pipe sticking and improve production efficiency, it is important to identify the state of the drill pipe sticking phenomenon.

[0003] The identification of the drill rod stuck state is an important prerequisite for ensuring the normal and stable operation of the drilling rig. However, due to the complex environment and noise interference, the recognition accuracy of the drill rod stuck phenomenon is poor. At present, there are relatively few studies on the identification of drill rod stuck in coal mine drilling rigs. The identification method of the drill rod stuck state is usually based on deep learning or image recognition. There is no in-depth consideration of the associated parameters of the drilling tool during drilling, such as torque, thrust pressure, thrust displacement, etc., which affects the accuracy of the identification of the drill rod stuck state and reduces the reliability. Summary of the invention

[0004] Based on this, it is necessary to provide a mining drill rig stuck drill state monitoring method based on feature stationarity test and cointegration analysis, which can improve the recognition accuracy and reliability of the drill rod stuck drill state in order to address the above technical problems.

[0005] A method for monitoring the stuck drill state of a mining drill based on characteristic stationarity test and cointegration analysis, the method comprising:

[0006] Real-time acquisition of drill-stuck status monitoring data of a drill rod of a drilling rig for a preset time period, wherein the drill-stuck status monitoring data includes vibration acceleration data, drilling pressure data, drilling displacement data, and rotational torque data, and the vibration acceleration data includes vibration acceleration signals in the horizontal direction, axial direction, and vertical direction of the drill rod;

[0007] Extracting time domain features and time frequency domain features from the vibration acceleration data to obtain time domain features and time frequency domain features;

[0008] taking the drilling pressure data, the drilling displacement data and the rotation torque data as statistical variables;

[0009] Performing stationarity tests on the time domain features, the time-frequency domain features and the statistical variables respectively, and dividing the stationary signals and non-stationary signals corresponding to the time domain features, the time-frequency domain features and the statistical variables;

[0010] Performing principal component analysis on the time domain features, the time-frequency domain features, and the stationary signals corresponding to the statistical variables respectively to obtain principal component analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables;

[0011] Performing cointegration analysis on the non-stationary signals corresponding to the time domain features, the time-frequency domain features, and the statistical variables respectively, to obtain cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables;

[0012] Perform principal component analysis based on the time domain characteristics, the time-frequency domain characteristics, and the principal component analysis residuals and cointegration analysis residuals corresponding to the statistical variables to obtain the T corresponding to the stuck drill status monitoring data. 2 Statistics;

[0013] According to the T 2 The statistics are compared and analyzed with a preset pipe-stuck state threshold value to monitor the time when the drill rod is stuck.

[0014] In one embodiment, the preset stuck drill state threshold is determined as follows:

[0015] Acquire historical vibration acceleration data, historical drilling pressure data, historical drilling displacement data and historical rotation torque data of a drill pipe of a sample drilling rig as sample data, wherein the historical vibration acceleration data includes vibration acceleration signals of the drill pipe in the horizontal direction, the axial direction and the vertical direction;

[0016] Extract time domain features and time frequency domain features from the historical vibration acceleration data to obtain time domain feature samples and time frequency domain feature samples;

[0017] Taking the historical drilling pressure data, the historical drilling displacement data and the historical rotation torque data as statistical variable samples;

[0018] Performing stationarity tests on the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples respectively, and dividing the stationary signals and non-stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples;

[0019] Performing principal component analysis on the stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples respectively, to obtain principal component analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples;

[0020] Performing cointegration analysis on the non-stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples respectively, to obtain cointegration analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples;

[0021] Perform principal component analysis based on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples to obtain T 2 Statistics sample;

[0022] According to the T 2 The control limits of the statistical sample are used as the stuck pipe state thresholds.

[0023] In one embodiment, performing principal component analysis on the time domain features, the time-frequency domain features, and the stationary signals corresponding to the statistical variables to obtain principal component analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables respectively includes:

[0024] The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the time domain feature to obtain the principal component analysis residual corresponding to the time domain feature. The expression of the principal component analysis residual corresponding to the time domain feature is:

[0025]

[0026] in, is the principal component analysis residual of the time domain characteristic stationary signal sequence, is the standardized time domain characteristic stationary signal sequence variable, P m =(P1,P2,…,P M ) is the variable coefficient of the time domain characteristic stationary signal sequence, M is the length of the time domain characteristic stationary signal sequence;

[0027] The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the time-frequency domain feature to obtain the principal component analysis residual corresponding to the time-frequency domain feature. The expression of the principal component analysis residual corresponding to the time-frequency domain feature is:

[0028]

[0029] in, is the residual of principal component analysis of the characteristic stationary signal sequence in the time-frequency domain, P TP is the variable coefficient of the characteristic stationary signal sequence in the time-frequency domain, is the standardized time-frequency domain characteristic stationary signal sequence variable;

[0030] The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the statistical variable to obtain the principal component analysis residual corresponding to the statistical variable. The expression of the principal component analysis residual corresponding to the statistical variable is:

[0031]

[0032] in, is the residual of principal component analysis of the stationary signal sequence of statistical variables, P S is the coefficient of the stationary signal sequence of statistical variables, is the standardized statistically stable signal sequence variable.

[0033] In one embodiment, the performing cointegration analysis on the non-stationary signals corresponding to the time domain features, the time-frequency domain features, and the statistical variables respectively to obtain the cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables includes:

[0034] The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the time domain feature to obtain the cointegration analysis residual corresponding to the time domain feature. The expression of the cointegration analysis residual corresponding to the time domain feature is:

[0035]

[0036] in, is the residual of cointegration analysis of non-stationary signal sequence with time domain characteristics, is the time series matrix of the non-stationary signal with time domain characteristics, is the cointegration vector of the time domain characteristic non-stationary signal sequence, the superscript T represents the transposition, and N is the length of the time domain characteristic non-stationary signal sequence;

[0037] The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the time-frequency domain features to obtain the cointegration analysis residual corresponding to the time-frequency domain features. The expression of the cointegration analysis residual corresponding to the time-frequency domain features is:

[0038]

[0039] in, is the residual of cointegration analysis of non-stationary signal sequence with time-frequency characteristics, δ TP is the cointegration vector of the non-stationary signal sequence with time-frequency characteristics, is the standardized time-frequency domain characteristic non-stationary signal sequence variable;

[0040] The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the statistical variable to obtain the cointegration analysis residual corresponding to the statistical variable. The expression of the cointegration analysis residual corresponding to the statistical variable is:

[0041]

[0042] in, is the residual of cointegration analysis of non-stationary signal sequence of statistical variables, δ S is the cointegration vector of the non-stationary signal of statistical variables, is the standardized statistical non-stationary signal sequence variable.

[0043] In one embodiment, the principal component analysis is performed based on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables to obtain T 2 Statistics, including:

[0044] The principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics to obtain the comprehensive principal component analysis residuals of the time domain characteristics;

[0045] The principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time-frequency domain features to obtain the comprehensive principal component analysis residuals of the time-frequency domain features;

[0046] Using the principal component analysis method to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the statistical variables, to obtain the comprehensive principal component analysis residuals of the statistical variables;

[0047] The principal component analysis method is used to perform principal component analysis on the comprehensive principal component analysis residuals of the time domain features, the comprehensive principal component analysis residuals of the time-frequency domain features, and the comprehensive principal component analysis residuals of the statistical variables to obtain T 2 Statistics.

[0048] In one embodiment, the T 2 The expression of the statistic is:

[0049]

[0050]

[0051] in, is the T of the stuck drill status monitoring data at the i-th moment 2 Statistics, t i is the score vector, λ is the eigenvalue, y i is the variable vector corresponding to the stuck drill status monitoring data at the i-th moment, q is the weight vector, and Q is the length of the stuck drill status monitoring data.

[0052] In one embodiment, the steps of principal component analysis are:

[0053] Step 1: perform data standardization on the signal sequence for principal component analysis to obtain a standardized signal sequence X = {X1, X2, ...X i …, X I}, the standardized expression is: Among them, X i is the i-th standardized signal, μ is the average value of the signal sequence, x i is the i-th signal to be standardized, σ is the variance of the signal sequence, and I is the length of the signal sequence;

[0054] Step 2: construct the covariance matrix S according to the standardized signal sequence X. The expression of the covariance matrix S is: Among them, the superscript T represents transposition;

[0055] Step 3: According to the covariance matrix S, the standardized signal sequence X is feature decomposed into the principal component subspace and the residual subspace, and the expression is:

[0056]

[0057] in, is the projection of the principal subspace, R is the score moment of the signal sequence in the subspace, P is the load matrix of the signal sequence in the subspace; E is the residual matrix;

[0058] Step 4: Determine the principal component analysis residual corresponding to the signal sequence according to the principal component subspace and the residual subspace.

[0059] In one embodiment, the steps of the cointegration analysis method include:

[0060] Step 1: Determine the homogeneity of the time series: For a time series X t , if after one difference, it becomes a stationary series, then the time series X t It has the characteristics of first-order integration, denoted as I(1);

[0061] Step 2: Construct a cointegration model: For a set of M×N non-stationary time series matrix X=[X1,X2,…,X N ], X t =(x1,x2,…,x M ), X t ~I(1), where N is the number of non-stationary time series and M is the number of sampling points, then there exists a cointegration vector β t =[β1,β2,…,β M ] T The non-stationary variables are expressed as linear combinations, expressed as:

[0062]

[0063] Among them, ε t If is the residual sequence, then a cointegration model is established, and the expression is:

[0064] γ=XB

[0065] Where B = [β 1 ,β 2 ,…,β R ] is the N×R cointegration coefficient matrix, γ is the M×R stationary residual matrix, R is the number of cointegration vectors, R≤N;

[0066] Step 3, build X t The VAR model is expressed as:

[0067] X t =∏ t X t-1 +…+∏ p X t-p +c+μ t

[0068] Among them, t is an N×R coefficient matrix, μ t is an N×1 white noise vector, which follows N(0,Ξ) distribution, c is an N×1 constant vector, and p is the order of the VAR model;

[0069] Step 4: Construct a vector error correction model:

[0070] X t Subtract X from both sides of the VAR model expression t-1 The vector error correction model of non-stationary time series is obtained, and the expression is:

[0071]

[0072] Where ΔX t is the corrected non-stationary time series, I N is the N×N identity matrix;

[0073] Decompose Г:

[0074] G = AB T

[0075] Where A is an N×R column-full rank matrix and B is an N×R column-full rank matrix, then:

[0076]

[0077] Therefore, the residual sequence γ is obtained t-1 :

[0078]

[0079] Among them, X t is a first-order integral, so ΔX t and ΔX t-i is also stable, then The elements in are stationary, and accordingly, B T X t-1 The elements in are also stable; therefore, the linear combination coefficient matrix B is the cointegration coefficient matrix, and the column vector in the cointegration coefficient matrix B is the cointegration vector;

[0080] Step 5, obtain the cointegration coefficient matrix B: Calculate the maximum likelihood estimate of the cointegration test and the cointegration coefficient matrix, then ΔX t The probability density function f(ΔX t )for:

[0081]

[0082] Where, Ξ is a variable;

[0083] Then the probability density function of the non-stationary time series matrix is:

[0084]

[0085] in,

[0086] Take the logarithm of the probability density function of the non-stationary time series matrix and establish the maximum likelihood function L(Ω1,…,Ω P-1 ,A,B,Ξ), the expression is:

[0087]

[0088] The cointegration coefficient matrix B is estimated by the maximum likelihood function; the cointegration coefficient matrix B is transformed into a characteristic equation for solving the characteristic equation, and the expression of the characteristic equation is:

[0089]

[0090] in, λ is the eigenvalue, Θ i and Φ i is the weight coefficient;

[0091] Step 6: Construct a non-stationary residual sequence: weight coefficient Θ i and Φ iThe least square method is used for estimation, and the characteristic vector matrix is ​​obtained by solving the characteristic equation. The cointegration vectors in the cointegration coefficient matrix B are included in the characteristic vector matrix; according to the number of cointegration vectors R and the characteristic vector matrix, the cointegration coefficient matrix B is determined to be [β 1 ,P 2 ,…,P R ], the stationary residual sequence is expressed as:

[0092]

[0093] Step 7: Determine the time series X according to the non-stationary residual sequence t The corresponding cointegration analysis residuals.

[0094] The above-mentioned method for monitoring the stuck drill status of a mining drill rig based on feature stationarity test and cointegration analysis obtains the stuck drill status monitoring data of the drill rod of the drilling rig in a preset time period in real time, wherein the stuck drill status monitoring data includes vibration acceleration data, drilling pressure data, drilling displacement data and rotational torque data, and the vibration acceleration data includes vibration acceleration signals in the horizontal direction, axial direction and vertical direction of the drill rod. The vibration acceleration data is subjected to time domain feature and time frequency domain feature extraction to obtain time domain feature and time frequency domain feature, and the drilling pressure data, the drilling displacement data and the rotational torque data are used as statistical variables, and the time domain feature, the time frequency domain feature and the statistical variable are subjected to stationarity test respectively to divide the time domain The stationary signal and the non-stationary signal corresponding to the time-frequency domain feature and the statistical variable are respectively subjected to principal component analysis to obtain the time-domain feature, the time-frequency domain feature and the principal component analysis residuals corresponding to the statistical variable; the time-domain feature, the time-frequency domain feature and the non-stationary signal corresponding to the statistical variable are respectively subjected to cointegration analysis to obtain the time-domain feature, the time-frequency domain feature and the cointegration analysis residuals corresponding to the statistical variable; the principal component analysis is performed based on the principal component analysis residuals and the cointegration analysis residuals corresponding to the time-domain feature, the time-frequency domain feature and the statistical variable to obtain the T corresponding to the stuck drill status monitoring data. 2 Statistics, according to the T 2 The statistics are compared and analyzed with a preset drill-stuck state threshold value to monitor the time when the drill rod is stuck, thereby improving the accuracy and reliability of identifying the drill-stuck state of the drill rod. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] Figure 1 A schematic diagram of a flow chart of a method for monitoring a stuck drill state of a mining drill based on feature stationarity test and cointegration analysis in one embodiment;

[0096] Figure 2It is a processing flow chart of a method for monitoring a stuck drill state of a mining drill based on feature stationarity test and cointegration analysis in one embodiment;

[0097] Figure 3-Figure 6 It is a waveform diagram of drill pipe status monitoring data at a moment in time;

[0098] Figure 7-10 is a waveform diagram of drill pipe status monitoring data at another moment of the drill pipe;

[0099] Fig.11 It is the time domain characteristic residual waveform diagram of the drill pipe before failure;

[0100] Fig.12 It is the residual error waveform of the time-frequency domain characteristics of the drill pipe before failure;

[0101] Fig.13 It is the residual waveform of drill pipe pressure, displacement and torque before failure;

[0102] Fig.14 Schematic diagram of monitoring the stuck drill pipe status. DETAILED DESCRIPTION

[0103] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0104] In one embodiment, Figure 1-2 As shown, a method for monitoring the stuck drill state of a mining drill based on characteristic stationarity test and cointegration analysis is provided. The method is applied to a terminal as an example for explanation, and includes the following steps:

[0105] Step S220, real-time acquisition of drill sticking status monitoring data of the drill rod of the drilling rig for a preset time period, the drill sticking status monitoring data including vibration acceleration data, drilling pressure data, drilling displacement data and rotational torque data, the vibration acceleration data including vibration acceleration signals of the drill rod in the horizontal direction, axial direction and vertical direction.

[0106] Step S240, extracting the time domain features and the time frequency domain features of the vibration acceleration data to obtain the time domain features and the time frequency domain features.

[0107] Step S260: taking the drilling pressure data, drilling displacement data and rotation torque data as statistical variables.

[0108] Step S280, performing stationarity tests on the time domain features, the time-frequency domain features and the statistical variables respectively, and dividing the stationary signals and non-stationary signals corresponding to the time domain features, the time-frequency domain features and the statistical variables.

[0109] Step S300, principal component analysis is performed on the time domain features, time-frequency domain features and stationary signals corresponding to the statistical variables to obtain principal component analysis residuals corresponding to the time domain features, time-frequency domain features and statistical variables.

[0110] In one embodiment, principal component analysis is performed on the time domain features, the time-frequency domain features, and the stationary signals corresponding to the statistical variables to obtain principal component analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables, including:

[0111] The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the time domain characteristics to obtain the principal component analysis residual corresponding to the time domain characteristics. The expression of the principal component analysis residual corresponding to the time domain characteristics is:

[0112]

[0113] in, is the principal component analysis residual of the time domain characteristic stationary signal sequence, is the standardized time domain characteristic stationary signal sequence variable, P m =(P1,P2,…,P M ) is the variable coefficient of the time domain characteristic stationary signal sequence, M is the length of the time domain characteristic stationary signal sequence;

[0114] The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the time-frequency domain characteristics to obtain the principal component analysis residual corresponding to the time-frequency domain characteristics. The expression of the principal component analysis residual corresponding to the time-frequency domain characteristics is:

[0115]

[0116] in, is the residual of principal component analysis of the characteristic stationary signal sequence in the time-frequency domain, P TP is the variable coefficient of the characteristic stationary signal sequence in the time-frequency domain, is the standardized time-frequency domain characteristic stationary signal sequence variable;

[0117] The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the statistical variable to obtain the principal component analysis residual corresponding to the statistical variable. The expression of the principal component analysis residual corresponding to the statistical variable is:

[0118]

[0119] in, is the residual of principal component analysis of the stationary signal sequence of statistical variables, P t is the coefficient of the stationary signal sequence of statistical variables, is the standardized statistically stable signal sequence variable.

[0120] Among them, principal component analysis (PCA) is an effective method to achieve data dimensionality reduction and construct correlations between stationary sequences.

[0121] In one embodiment, the steps of principal component analysis are:

[0122] Step 1: perform data standardization on the signal sequence for principal component analysis to obtain a standardized signal sequence X = {X1, X2, ...X i …, X I}, the standardized expression is: Among them, X i is the i-th standardized signal, μ is the average value of the signal sequence, x i is the i-th signal to be standardized, σ is the variance of the signal sequence, and I is the length of the signal sequence;

[0123] Step 2: construct the covariance matrix S according to the standardized signal sequence X. The expression of the covariance matrix S is: Among them, the superscript T represents transposition;

[0124] Step 3: According to the covariance matrix S, the standardized signal sequence X is feature decomposed into the principal component subspace and the residual subspace, and the expression is:

[0125]

[0126] in, is the projection of the principal subspace, R is the score moment of the signal sequence in the subspace, P is the load matrix of the signal sequence in the subspace; E is the residual matrix;

[0127] Step 4: Determine the principal component analysis residual corresponding to the signal sequence according to the principal component subspace and the residual subspace.

[0128] Step S320, performing cointegration analysis on the non-stationary signals corresponding to the time domain features, the time-frequency domain features and the statistical variables, respectively, to obtain cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features and the statistical variables.

[0129] In one embodiment, cointegration analysis is performed on the non-stationary signals corresponding to the time domain features, the time-frequency domain features, and the statistical variables, respectively, to obtain the cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables, including:

[0130] The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the time domain characteristics to obtain the cointegration analysis residual corresponding to the time domain characteristics. The expression of the cointegration analysis residual corresponding to the time domain characteristics is:

[0131]

[0132] in, is the residual of cointegration analysis of non-stationary signal sequence with time domain characteristics, is the time series matrix of the non-stationary signal with time domain characteristics, is the cointegration vector of the time domain characteristic non-stationary signal sequence, the superscript T represents the transposition, and N is the length of the time domain characteristic non-stationary signal sequence;

[0133] The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the time-frequency domain characteristics to obtain the cointegration analysis residual corresponding to the time-frequency domain characteristics. The expression of the cointegration analysis residual corresponding to the time-frequency domain characteristics is:

[0134]

[0135] in, is the residual of cointegration analysis of non-stationary signal sequence with time-frequency characteristics, δ TP is the cointegration vector of the non-stationary signal sequence with time-frequency characteristics, is the standardized time-frequency domain characteristic non-stationary signal sequence variable;

[0136] The cointegration analysis method is used to perform cointegration analysis on the non-stationary signals corresponding to the statistical variables to obtain the cointegration analysis residuals corresponding to the statistical variables. The expression of the cointegration analysis residuals corresponding to the statistical variables is:

[0137]

[0138] in, is the residual of cointegration analysis of non-stationary signal sequence of statistical variables, δ S is the cointegration vector of the non-stationary signal of statistical variables, is the standardized statistical non-stationary signal sequence variable.

[0139] Among them, cointegration analysis is an effective method to describe the relationship between non-stationary variables.

[0140] In one embodiment, the steps of cointegration analysis include:

[0141] Step 1: Determine the homogeneity of the time series: For a time series X t , if after one difference, it becomes a stationary series, then the time series X t It has the characteristics of first-order integration, denoted as I(1);

[0142] Step 2: Construct a cointegration model: For a set of M×N non-stationary time series matrix X=[X1,X2,…,X N ], X t =(x1,x2,…,xM ), X t ~I(1), where N is the number of non-stationary time series and M is the number of sampling points, then there exists a cointegration vector β t =[β1,β2,…,β M ] T The non-stationary variables are expressed as linear combinations, expressed as:

[0143]

[0144] Among them, ε t If is the residual sequence, then a cointegration model is established, and the expression is:

[0145] γ=XB

[0146] Where B = [β 1 ,β 2 ,…,β R ] is the N×R cointegration coefficient matrix, γ is the M×R stationary residual matrix, R is the number of cointegration vectors, R≤N;

[0147] Step 3, build X t The VAR model is expressed as:

[0148] X t =∏ t X t-1 +…+∏ p X t-p +c+μ t

[0149] Among them, t is an N×R coefficient matrix, μ t is an N×1 white noise vector, which follows N(0,Ξ) distribution, c is an N×1 constant vector, and p is the order of the VAR model;

[0150] Step 4: Construct a vector error correction model:

[0151] X t Subtract X from both sides of the VAR model expression t-1 The vector error correction model of non-stationary time series is obtained, and the expression is:

[0152]

[0153] Where ΔX t is the corrected non-stationary time series, I N is the N×N identity matrix;

[0154] Decompose Г:

[0155] G = AB T

[0156] Where A is an N×R column-full rank matrix and B is an N×R column-full rank matrix, then:

[0157]

[0158] Therefore, the residual sequence γ is obtained t-1 :

[0159]

[0160] Among them, X t is a first-order integral, so ΔX t and ΔX t-i is also stable, then The elements in are stationary, and accordingly, B T X t-1 The elements in are also stable; therefore, the linear combination coefficient matrix B is the cointegration coefficient matrix, and the column vector in the cointegration coefficient matrix B is the cointegration vector;

[0161] Step 5, obtain the cointegration coefficient matrix B: Calculate the maximum likelihood estimate of the cointegration test and the cointegration coefficient matrix, then ΔX t The probability density function f(ΔX t )for:

[0162]

[0163] Where, Ξ is a variable;

[0164] Then the probability density function of the non-stationary time series matrix is:

[0165]

[0166] in,

[0167] Take the logarithm of the probability density function of the non-stationary time series matrix and establish the maximum likelihood function L(Ω1,…,Ω P-1 ,A,B,Ξ), the expression is:

[0168]

[0169] The cointegration coefficient matrix B is estimated by the maximum likelihood function; the cointegration coefficient matrix B is transformed into the characteristic equation, and the expression of the characteristic equation is:

[0170]

[0171] in, λ is the eigenvalue, Θ i and Φ i is the weight coefficient;

[0172] Step 6: Construct a non-stationary residual sequence: weight coefficient Θ i and Φ i The least squares method is used for estimation, and the characteristic vector matrix is ​​obtained by solving the characteristic equation. The cointegration vector in the cointegration coefficient matrix B is included in the characteristic vector matrix; according to the number of cointegration vectors R and the characteristic vector matrix, the cointegration coefficient matrix B is determined to be [β 1 ,β 2 ,…,β R ], the stationary residual sequence is expressed as:

[0173]

[0174] Step 7: Determine the time series X based on the non-stationary residual sequence t The corresponding cointegration analysis residuals.

[0175] Step S340, performing principal component analysis based on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, and obtaining the T corresponding to the stuck drill status monitoring data. 2 Statistics.

[0176] In one embodiment, principal component analysis is performed based on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics, time-frequency domain characteristics, and statistical variables to obtain T 2 Statistics, including:

[0177] The principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics to obtain the comprehensive principal component analysis residuals of the time domain characteristics; the principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time-frequency domain characteristics to obtain the comprehensive principal component analysis residuals of the time-frequency domain characteristics; the principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the statistical variables to obtain the comprehensive principal component analysis residuals of the statistical variables; the principal component analysis method is used to perform principal component analysis on the comprehensive principal component analysis residuals of the time domain characteristics, the comprehensive principal component analysis residuals of the time-frequency domain characteristics and the comprehensive principal component analysis residuals of the statistical variables to obtain T 2 Statistics.

[0178] In one embodiment, T 2 The expression of the statistic is:

[0179]

[0180]

[0181] in, is the T of the stuck drill status monitoring data at the i-th moment 2 Statistics, t i is the score vector, λ is the eigenvalue, y i is the variable vector corresponding to the stuck drill status monitoring data at the i-th moment, q is the weight vector, and Q is the length of the stuck drill status monitoring data.

[0182] Step S360, according to T 2 The statistics are compared and analyzed with the preset pipe-stuck state threshold to monitor the time when the drill pipe is stuck.

[0183] It should be understood that T 2 Statistics can be used to evaluate the main changes in the observed quantity. If the change in the latent variable is greater than the normal change, it means that the equipment status is abnormal.

[0184] In one embodiment, the preset stuck pipe state threshold is determined as follows:

[0185] The historical vibration acceleration data, historical drilling pressure data, historical drilling displacement data and historical rotation torque data of the drill pipe of the sample drilling rig are obtained as sample data, and the historical vibration acceleration data include the vibration acceleration signals of the drill pipe in the horizontal direction, the axial direction and the vertical direction; the time domain features and the time frequency domain features of the historical vibration acceleration data are extracted to obtain the time domain feature samples and the time frequency domain feature samples; the historical drilling pressure data, the historical drilling displacement data and the historical rotation torque data are used as statistical variable samples; the time domain feature samples, the time frequency domain feature samples and the statistical variable samples are respectively tested for stationarity, and the time domain feature samples, the time frequency domain feature samples and the statistical variable samples are divided into and the stationary signals and non-stationary signals corresponding to the statistical variable samples; perform principal component analysis on the stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples, respectively, and obtain the principal component analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples; perform cointegration analysis on the non-stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples, respectively, and obtain the cointegration analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples; perform principal component analysis on the principal component analysis residuals and the cointegration analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples, and obtain T 2 Statistics sample; according to T 2 The control limits of the statistical sample are used as the stuck pipe state thresholds.

[0186] Among them, the T corresponding to the stuck drill status monitoring data detected by PCA 2 When the statistic exceeds the preset stuck drill status threshold, the operating status of the drill pipe can be determined. It can be expressed as:

[0187]

[0188] Where F(k,mk,α) is an F distribution with mk degrees of freedom and confidence level α; as control limits at confidence level α.

[0189] In one embodiment, the operation state of the drill rod can be divided into two operation states: normal operation and stuck drill. The stuck drill state refers to the state where the drill rod is stuck in the coal rock and can neither continue drilling nor be successfully pulled out. The normal state refers to the state where the drill rod drills normally.

[0190] The preset stuck drill status threshold can be combined with T 2 The setting is based on a comprehensive analysis of statistical samples, process regulations and experience. For example, the drill stuck pressure threshold can be set to 50N.

[0191] It should be understood that the preset stuck drill state threshold corresponds to the stuck drill state judgment. For example, the T corresponding to the stuck drill state monitoring data 2 If the maximum value of the statistic is greater than the preset stuck drill state threshold, the drill pipe is in a stuck drill state; the T corresponding to the stuck drill state monitoring data 2 If the maximum value of the statistic is less than the preset drill-stuck state threshold, the drill pipe is in a normal operating state.

[0192] The above-mentioned method for monitoring the stuck drill status of a mining drill rig based on feature stationarity test and cointegration analysis obtains the stuck drill status monitoring data of the drill rod of the drill rig in a preset time period in real time. The stuck drill status monitoring data includes vibration acceleration data, drilling pressure data, drilling displacement data and rotation torque data. The vibration acceleration data includes vibration acceleration signals in the horizontal direction, axial direction and vertical direction of the drill rod. The time domain features and time frequency domain features of the vibration acceleration data are extracted to obtain the time domain features and time frequency domain features. The drilling pressure data, drilling displacement data and rotation torque data are used as statistical variables, and the time domain features, time frequency domain features and statistical variables are respectively extracted. Stationarity test, divide the stationary signals and non-stationary signals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, perform principal component analysis on the stationary signals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, respectively, obtain the principal component analysis residuals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, perform cointegration analysis on the non-stationary signals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, respectively, obtain the cointegration analysis residuals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics, time-frequency domain characteristics and statistical variables, and obtain the T corresponding to the stuck drill status monitoring data 2 Statistics, according to T 2The statistics are compared and analyzed with a preset drill-stuck state threshold to monitor the time when the drill rod is stuck, thereby improving the accuracy and reliability of identifying the drill-stuck state of the drill rod.

[0193] In order to verify the mining drill rig stuck drill state monitoring method based on characteristic stationarity test and cointegration analysis in this application, the following verification is carried out:

[0194] Use a properly assembled and normally operating drill pipe, apply a certain speed and thrust pressure to the drill pipe through the drilling rig experimental platform, make it start running, drill into the test block, and stop running when the drill pipe gets stuck. Collect its vibration acceleration signal, and select the data of the two moments when the drill pipe is stuck and not stuck as input. The waveforms of the vibration acceleration of the drill pipe at two different moments are shown in the figure below. Figure 3-Figure 6 as well as Figure 7-10 As shown in Figure 2. The signal acquisition frequency is 1kHz.

[0195] The T was established by the residuals after PCA analysis. 2 The statistic is compared with the preset drill-stuck threshold to determine whether the drill rod is stuck. If the maximum value of the statistic is greater than the preset drill-stuck threshold, considering the operation status of drilling into the coal body, when the statistic gradually increases, the corresponding drill rod generates more coal dust, and the risk of drill-stuck increases. The drill rod is in the drill-stuck state. Here, T 2 The control limit threshold for the statistic is 6.8286.

[0196] For drills that are not stuck Figure 7-10 Data is extracted and its time domain and time frequency domain features are used for subsequent calculations. In order to reduce the impact of outliers and facilitate calculations, the signal is downsampled before feature extraction to reduce the sampling frequency and the amount of data, but it does not affect the evaluation results of the operating status of the drill pipe. After interval sampling, there are 100 values ​​in each sample point. The extracted time domain and time frequency domain features are respectively obtained based on PCA analysis and cointegration analysis. The residuals are as follows: Fig.11 , Fig.12 The drilling pressure, drilling displacement and rotation torque of the drill pipe are established based on PCA analysis and cointegration analysis, as shown in Fig.13 As shown. PCA analysis is performed on the residuals of the generated population to obtain T 2 Statistics and preset stuck drill status thresholds are as follows Fig.14 shown.

[0197] Set the preset stuck drill status threshold by T 2 The control limit threshold of the statistic indicates that T 2 The control limit threshold of the statistic is 6.8286. 2 The change curve is judged by the drill stuck threshold value. The judgment result is as follows: Fig.14 shown.

[0198] pass Figure 3-Figure 6 It is found that the length of the collected signal is 577002, every 10000 values ​​represent a sampling point, and the interval between each sampling point is 1s, so the total sampling time is 57.7s. The preset drill stuck state threshold is determined by T 2 The control limit threshold of the statistic indicates that T 2 The control limit threshold of the statistic is 6.8286, and the intersection with the signal is the 5570th sample point, corresponding to a usage time of 55.7s.

[0199] Fig.11 It is the residual graph, which shows the residuals obtained after applying PCA and cointegration analysis to the extracted time domain features of the data; Fig.12 This is the residual waveform diagram, from which we can see the residuals obtained after using PCA and cointegration analysis on the time-frequency domain features extracted from the data. Fig.14 PCA analysis is used to obtain T for the overall residual 2 Statistics and control limit threshold diagram. The preset stuck drill status threshold is determined by T 2 The control limit threshold of the statistic indicates that setting T 2 The control limit threshold of the statistic is 6.8286. Fig.14 As shown, a cross position can be obtained, which is the 5570th sample point, and the corresponding drill rod running time is 55.7s. It can be found that the mining drill rig stuck drill state monitoring method based on characteristic stationarity test and cointegration analysis of the present application can well judge the running state of the drill rod.

[0200] In summary, the method for monitoring the stuck drill status of a mining drill based on feature stationarity test and cointegration analysis provided in the embodiment of the present application can effectively judge the operating status of the drill rod, can timely and accurately judge the time when the drill rod is stuck, and can achieve good monitoring of the operating status of the drill rod.

[0201] It should be understood that although Figure 1 The steps in the flowchart are shown in sequence as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified in this document, there is no strict order restriction for the execution of these steps, and these steps can be executed in other orders. Moreover, Figure 1 At least part of the steps may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least part of the sub-steps or stages of other steps.

[0202] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0203] The above-mentioned embodiments only express several implementation methods of the present application, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for a person of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. Therefore, the protection scope of the patent of the present application shall be subject to the attached claims.

Claims

1. A method for monitoring the stuck drill state of a mining drill based on characteristic stationarity test and cointegration analysis, characterized in that: The method comprises: Real-time acquisition of drill-stuck status monitoring data of a drill rod of a drilling rig for a preset time period, wherein the drill-stuck status monitoring data includes vibration acceleration data, drilling pressure data, drilling displacement data, and rotational torque data, and the vibration acceleration data includes vibration acceleration signals in the horizontal direction, axial direction, and vertical direction of the drill rod; Extracting time domain features and time frequency domain features from the vibration acceleration data to obtain time domain features and time frequency domain features; taking the drilling pressure data, the drilling displacement data and the rotation torque data as statistical variables; Performing stationarity tests on the time domain features, the time-frequency domain features and the statistical variables respectively, and dividing the stationary signals and non-stationary signals corresponding to the time domain features, the time-frequency domain features and the statistical variables; Performing principal component analysis on the time domain features, the time-frequency domain features, and the stationary signals corresponding to the statistical variables respectively to obtain principal component analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables; Performing cointegration analysis on the non-stationary signals corresponding to the time domain features, the time-frequency domain features, and the statistical variables respectively, to obtain cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features, and the statistical variables; Perform principal component analysis based on the time domain characteristics, the time-frequency domain characteristics, and the principal component analysis residuals and cointegration analysis residuals corresponding to the statistical variables to obtain the T corresponding to the stuck drill status monitoring data. 2 Statistics; According to the T 2 Comparing and analyzing the statistic with a preset drill-stuck state threshold value to monitor the time when the drill rod is stuck; The principal component analysis is performed based on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics, the time-frequency domain characteristics and the statistical variables to obtain T 2 Statistics, including: The principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain characteristics to obtain the comprehensive principal component analysis residuals of the time domain characteristics; The principal component analysis method is used to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the time-frequency domain features to obtain the comprehensive principal component analysis residuals of the time-frequency domain features; Using the principal component analysis method to perform principal component analysis on the principal component analysis residuals and cointegration analysis residuals corresponding to the statistical variables, to obtain the comprehensive principal component analysis residuals of the statistical variables; The principal component analysis method is used to perform principal component analysis on the comprehensive principal component analysis residuals of the time domain features, the comprehensive principal component analysis residuals of the time-frequency domain features, and the comprehensive principal component analysis residuals of the statistical variables to obtain T 2 Statistics; The T 2 The expression of the statistic is: in, is the T of the stuck drill status monitoring data at the i-th moment 2 Statistics, t i is the score vector, λ is the eigenvalue, y i is the variable vector corresponding to the stuck drill status monitoring data at the i-th moment, q is the weight vector, and Q is the length of the stuck drill status monitoring data.

2. The method according to claim 1, characterized in that The preset stuck drill state threshold is determined as follows: Acquire historical vibration acceleration data, historical drilling pressure data, historical drilling displacement data and historical rotation torque data of a drill pipe of a sample drilling rig as sample data, wherein the historical vibration acceleration data includes vibration acceleration signals of the drill pipe in the horizontal direction, the axial direction and the vertical direction; Extract time domain features and time frequency domain features from the historical vibration acceleration data to obtain time domain feature samples and time frequency domain feature samples; Taking the historical drilling pressure data, the historical drilling displacement data and the historical rotation torque data as statistical variable samples; Performing stationarity tests on the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples respectively, and dividing the stationary signals and non-stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples; Performing principal component analysis on the stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples respectively, to obtain principal component analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples; Performing cointegration analysis on the non-stationary signals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples respectively, to obtain cointegration analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples; Perform principal component analysis based on the principal component analysis residuals and cointegration analysis residuals corresponding to the time domain feature samples, the time-frequency domain feature samples and the statistical variable samples to obtain T 2 Statistics sample; According to the T 2 The control limits of the statistical sample are used as the stuck pipe state thresholds.

3. The method according to claim 1, characterized in that The performing principal component analysis on the time domain features, the time-frequency domain features and the stationary signals corresponding to the statistical variables respectively to obtain the principal component analysis residuals corresponding to the time domain features, the time-frequency domain features and the statistical variables comprises: The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the time domain feature to obtain the principal component analysis residual corresponding to the time domain feature. The expression of the principal component analysis residual corresponding to the time domain feature is: in, is the principal component analysis residual of the time domain characteristic stationary signal sequence, is the standardized time domain characteristic stationary signal sequence variable, P m =(P1,P2,…,P M ) is the variable coefficient of the time domain characteristic stationary signal sequence, M is the length of the time domain characteristic stationary signal sequence; The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the time-frequency domain feature to obtain the principal component analysis residual corresponding to the time-frequency domain feature. The expression of the principal component analysis residual corresponding to the time-frequency domain feature is: in, is the residual of principal component analysis of the characteristic stationary signal sequence in the time-frequency domain, P TP is the variable coefficient of the characteristic stationary signal sequence in the time-frequency domain, is the standardized time-frequency domain characteristic stationary signal sequence variable; The principal component analysis method is used to perform principal component analysis on the stationary signal corresponding to the statistical variable to obtain the principal component analysis residual corresponding to the statistical variable. The expression of the principal component analysis residual corresponding to the statistical variable is: in, is the residual of principal component analysis of the stationary signal sequence of statistical variables, P S is the coefficient of the stationary signal sequence of statistical variables, is the standardized statistically stable signal sequence variable.

4. The method according to claim 3, characterized in that The performing cointegration analysis on the non-stationary signals corresponding to the time domain features, the time-frequency domain features and the statistical variables respectively to obtain the cointegration analysis residuals corresponding to the time domain features, the time-frequency domain features and the statistical variables comprises: The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the time domain feature to obtain the cointegration analysis residual corresponding to the time domain feature. The expression of the cointegration analysis residual corresponding to the time domain feature is: in, is the residual of cointegration analysis of non-stationary signal sequence with time domain characteristics, is the time series matrix of the non-stationary signal with time domain characteristics, is the cointegration vector of the time domain characteristic non-stationary signal sequence, the superscript T represents the transposition, and N is the length of the time domain characteristic non-stationary signal sequence; The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the time-frequency domain features to obtain the cointegration analysis residual corresponding to the time-frequency domain features. The expression of the cointegration analysis residual corresponding to the time-frequency domain features is: in, is the residual of cointegration analysis of non-stationary signal sequence with time-frequency characteristics, δ TP is the cointegration vector of the non-stationary signal sequence with time-frequency characteristics, is the standardized time-frequency domain characteristic non-stationary signal sequence variable; The cointegration analysis method is used to perform cointegration analysis on the non-stationary signal corresponding to the statistical variable to obtain the cointegration analysis residual corresponding to the statistical variable. The expression of the cointegration analysis residual corresponding to the statistical variable is: in, is the residual of cointegration analysis of non-stationary signal sequence of statistical variables, δ S is the cointegration vector of the non-stationary signal of statistical variables, is the standardized statistical non-stationary signal sequence variable.

5. The method according to claim 3, characterized in that: The steps of principal component analysis are: Step 1: perform data standardization on the signal sequence for principal component analysis to obtain a standardized signal sequence X = {X1, X2, ...X i …, X I }, the standardized expression is: Among them, X i is the i-th standardized signal, μ is the average value of the signal sequence, x i is the i-th signal to be standardized, σ is the variance of the signal sequence, and I is the length of the signal sequence; Step 2: construct the covariance matrix S according to the standardized signal sequence X. The expression of the covariance matrix S is: Among them, the superscript T represents transposition; Step 3: According to the covariance matrix S, the standardized signal sequence X is feature decomposed into the principal component subspace and the residual subspace, and the expression is: in, is the projection of the principal subspace, R is the score moment of the signal sequence in the subspace, P is the load matrix of the signal sequence in the subspace; E is the residual matrix; Step 4: Determine the principal component analysis residual corresponding to the signal sequence according to the principal component subspace and the residual subspace.

6. The method according to claim 4, characterized in that The steps of the cointegration analysis method include: Step 1: Determine the homogeneity of the time series: For a time series X t , if after one difference, it becomes a stationary series, then the time series X t It has the characteristics of first-order integration, denoted as I(1); Step 2: Construct a cointegration model: For a set of M×N non-stationary time series with cointegration relationship Matrix X = [X1,X2,…,X N ], X t =(x1,x2,…,x M ), X t ~I(1), where N is the number of non-stationary time series and M is the number of sampling points, then there exists a cointegration vector β t =[β1,β2,…,β M ] T The non-stationary variables are expressed as linear combinations, expressed as: Among them, ε t If is the residual sequence, then a cointegration model is established, and the expression is: γ=XB Where B = [β 1 ,β 2 ,…,β R ] is the N×R cointegration coefficient matrix, γ is the M×R stationary residual matrix, R is the number of cointegration vectors, R≤N; Step 3, build X t The VAR model is expressed as: X t =∏ t X t-1 +...+∏ p X t-p +c+μ t Among them, t is an N×R coefficient matrix, μ t is an N×1 white noise vector, which follows N(0,Ξ) distribution, c is an N×1 constant vector, and p is the order of the VAR model; Step 4: Construct a vector error correction model: X t Subtract X from both sides of the VAR model expression t-1 The vector error correction model of non-stationary time series is obtained, and the expression is: Where ΔX t is the corrected non-stationary time series, I N is the N×N identity matrix; Decompose Г: G=AB T Where A is an N×R column-full rank matrix and B is an N×R column-full rank matrix, then: Therefore, the residual sequence γ is obtained t-1 : Among them, X t is a first-order integral, so ΔX t and ΔX t-i is also stable, then The elements in are stationary, and accordingly, B T X t-1 The elements in are also stable; therefore, the linear combination coefficient matrix B is the cointegration coefficient matrix, and the column vector in the cointegration coefficient matrix B is the cointegration vector; Step 5, obtain the cointegration coefficient matrix B: Calculate the maximum likelihood estimate of the cointegration test and the cointegration coefficient matrix, then ΔX t The probability density function f(ΔX t )for: Where, Ξ is a variable; Then the probability density function of the non-stationary time series matrix is: in, Take the logarithm of the probability density function of the non-stationary time series matrix and establish the maximum likelihood function L(Ω1,…,Ω P-1 ,A,B,Ξ), the expression is: The cointegration coefficient matrix B is estimated by the maximum likelihood function; the cointegration coefficient matrix B is transformed into a characteristic equation for solving the characteristic equation, and the expression of the characteristic equation is: in, λ is the eigenvalue, Θ i and Φ i is the weight coefficient; Step 6: Construct a non-stationary residual sequence: weight coefficient Θ i and Φ i The least square method is used for estimation, and the characteristic vector matrix is ​​obtained by solving the characteristic equation. The cointegration vectors in the cointegration coefficient matrix B are included in the characteristic vector matrix; according to the number of cointegration vectors R and the characteristic vector matrix, the cointegration coefficient matrix B is determined to be [β 1 ,β 2 ,…,β R ], the stationary residual sequence is expressed as: Step 7: Determine the time series X according to the non-stationary residual sequence t The corresponding cointegration analysis residuals.

Images