CN110084185B - Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train - Google Patents

Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train Download PDF

Info

Publication number
CN110084185B
CN110084185B CN201910336833.8A CN201910336833A CN110084185B CN 110084185 B CN110084185 B CN 110084185B CN 201910336833 A CN201910336833 A CN 201910336833A CN 110084185 B CN110084185 B CN 110084185B
Authority
CN
China
Prior art keywords
matrix
vibration signal
small
snaking
energy
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201910336833.8A
Other languages
Chinese (zh)
Other versions
CN110084185A (en
Inventor
宁静
冉伟
陈春俊
李艳萍
张兵
方明宽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN201910336833.8A priority Critical patent/CN110084185B/en
Publication of CN110084185A publication Critical patent/CN110084185A/en
Application granted granted Critical
Publication of CN110084185B publication Critical patent/CN110084185B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Abstract

The invention discloses a method for quickly extracting small snaking running characteristics of a high-speed train. Aiming at the problem that the existing method cannot quickly identify the small snaking motion state (normal operation, small convergence, small divergence and large snaking) of the high-speed train, the invention adopts an average empirical mode decomposition (EEMD) method to decompose the preprocessed signals, converts the result into an energy matrix, then performs combined approximate diagonalization (JAD) processing on the energy matrix under the non-stationary condition, fuses the energy matrices of a plurality of sensors, and obtains the fused small snaking motion state characteristic of the high-speed train. The fusion characteristics are put into a least square method support vector machine for training and recognition, and the method is verified to be capable of quickly and accurately separating four running states of normal running, small-amplitude convergence, small-amplitude divergence and large-amplitude snaking of a high-speed train, so that the running safety of the train is guaranteed.

Description

Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train
Technical Field
The invention belongs to the technical field of train running state feature extraction, and particularly relates to a rapid extraction method for small-amplitude snaking running features of a high-speed train.
Background
The high-speed train rapidly develops in China, and by 2018, China is the country with the highest mileage of high-speed rail passing vehicles and mileage under construction in the world, and high-speed rails bring convenience to people's trip, but with the increase of the speed of the high-speed rails, the operation safety and stability become the focus of attention.
During the operation of the high-speed train, research shows that the high-speed train can generate unstable snaking movement after the running speed of the train exceeds a certain critical value. Hunting instability can produce large transverse wheel rail force, and can cause train derailment in severe cases. Therefore, safe operation of the hunting unstable high-speed train is significantly hindered. At present, no unified assessment standard exists for vehicle hunting instability at home and abroad, and most of the foreign relevant standards judge whether the hunting instability exists through axle transverse force, wheel rail transverse force, vehicle transverse acceleration and the like in the process of train running. The running safety monitoring standard of the railway passenger train in China also uses the index of the transverse acceleration of a framework to evaluate the transverse stability of the locomotive, and the transverse acceleration of a bogie of the train is regulated by the whole train test specification of the high-speed motor train unit to continuously reach or even exceed a limit value for 6 times after being filtered by 10Hz, so that the condition that the train loses the snaking stability can be judged. The existing peak monitoring method in China monitors the high-speed train based on the specifications so as to determine whether the hunting instability phenomenon exists in the running process of the train.
When experimental data are analyzed and processed, the train is found to be easy to have small snaking movement when running at high speed, namely small snaking. The small snaking refers to the part of the transverse acceleration signal of the framework, the amplitude of which does not reach the safety limit value in the running process of the train. Small snaking has many hazards, such as affecting passenger riding comfort, causing wheel track wear to be aggravated, reducing train service time, and bringing severe hazards to train safe operation, and small snaking abnormity is a sign of severe snaking instability. When the invention researches the small snaking state of the high-speed train, the invention discovers that the small snaking can be divided into two types according to the development of the small snaking motion of the train: small amplitude convergence and small amplitude divergence. The small-amplitude convergent state of the high-speed train can be slowly converged to a normal running state, and the transverse vibration amplitude of the small-amplitude divergent state of the high-speed train can be continuously enlarged until the high-speed train becomes snake instability, so that the high-speed train is greatly damaged in running, and therefore the two stages for distinguishing the small-amplitude snake in time are very important for the running stability of the train.
Referring to fig. 1, fig. 1(a) is a small snaking convergence diagram of a high-speed train, and fig. 1(b) is a small snaking divergence diagram of the high-speed train, it can be seen that the lateral amplitude of the bogie in fig. 1(a) is gradually reduced to a normal operation state, and the lateral amplitude of the bogie in fig. 1(b) is gradually increased to a snaking instability state. Therefore, it is necessary to distinguish the two states of small snaking in time.
Researchers have conducted relevant studies on the small snaking evolution of high-speed trains. Ning J provides a method for extracting small-amplitude snake-shaped evolution characteristics of a high-speed train based on multi-scale arrangement entropy and local tangent space arrangement, and the method can accurately extract the operation characteristics of a plurality of operation states, but needs longer consumed operation time. The running speed of the high-speed train is extremely high, and the characteristic extraction result cannot be quickly reflected to the train due to the overlong operation time, so that the train misses the best adjustment window. In order to solve the problems, the invention provides a method for quickly extracting small snaking running characteristics of a high-speed train.
Meanwhile, in the test, when data collection is carried out on the running state of the high-speed train, a multi-sensor network system mode is generally adopted, when the hunting instability phenomenon of the high-speed train and the small hunting phenomenon of the high-speed train are analyzed, a certain section of data information in a certain position sensor is selected by many researchers for research and analysis, and at the moment, whether the research conclusion can reflect all fault information of the train is uncertain. Therefore, aiming at the problem to be solved urgently in the research of the snaking motion of the high-speed train in China, the invention combines theoretical analysis and engineering practice, takes a small snaking evolution characteristic extraction method as a research core, integrates the vibration data of a plurality of position sensors of the train, extracts the small snaking evolution characteristic of the high-speed train according to the measurement data of the acceleration signal of the bogie frame when the high-speed train operates, and provides a reference basis for the state monitoring and safety early warning of the high-speed train in China.
Disclosure of Invention
The method for rapidly extracting the small snake-shaped running characteristic of the high-speed train solves the problems that the running state characteristic of the high-speed train in the prior art is not completely extracted and the small snake-shaped evolution characteristic of the train cannot be rapidly extracted.
In order to achieve the purpose of the invention, the invention adopts the technical scheme that: a method for quickly extracting small-amplitude snaking running characteristics of a high-speed train comprises the following steps:
s1, extracting an original vibration signal when the high-speed train runs, and preprocessing the original vibration signal;
s2, decomposing the preprocessed vibration signal into a plurality of IMF components through EEMD;
s3, constructing an energy matrix corresponding to the vibration signal according to the IMF component;
and S4, processing the energy matrix through a non-stationary state JADE method to obtain a fusion characteristic matrix, namely a small snaking running characteristic of the high-speed train. .
Further, the original vibration signals in step S1 include a 1-bit axle box transverse vibration signal, a 1-bit frame transverse vibration signal and a 4-bit frame transverse vibration signal;
further, the step S2 is specifically:
s21, adding a Gaussian white noise sequence omega (t) with the same amplitude to the preprocessed vibration signal a (t) to obtain an overall signal a' (t);
wherein the overall signal a' (t) is:
a'(t)＝a(t)+ω(t)
s22, decomposing the overall signal a' (t) added with white noise each time according to an EMD method, and determining each order component ci：
The decomposition formula of the overall signal a' (t) is as follows:
Figure BDA0002039410560000031
wherein: c. CiIs IMF component of each order;
i represents the ith IMF component;
r is a residual term;
n is the number of decomposed IMFs.
S23, adding different white noise sequences omega with the same amplitude each timej(t), repeating steps S21-S22 to obtain:
Figure BDA0002039410560000041
in the formula: a isj' (t) is the overall signal after the j-th added Gaussian white noise sequence;
n is the number of times of adding Gaussian white noise sequence;
ωj(t) is the j-th added Gaussian white noise sequence;
cijthe ith IMF component is decomposed when the white Gaussian noise sequence is added for the jth time;
rjand residual values are decomposed when the white Gaussian noise sequence is added for the jth time.
S24, according to the zero-mean principle of Gaussian white noise frequency, eliminating the influence of the Gaussian white noise on the IMF component, and obtaining the IMF component corresponding to the preprocessed original vibration signal as follows:
Figure BDA0002039410560000042
in the formula: c. Ci(t) is an IMF component corresponding to the ith original vibration signal;
cij(t) is the i-th IMF component obtained by EEMD decomposition of the original vibration signal.
Further, the step S3 is specifically:
s31, forming an IMF energy moment according to the IMF components as follows:
Figure BDA0002039410560000043
in the formula: p is the total number of sampling points;
m is a sampling point;
cithe IMF component corresponding to the selected vibration signal;
Δ t is the sampling period.
S32, according to the IMF energy moment, constructing a feature vector M of the normalized IMF energy moment as follows:
Figure BDA0002039410560000051
in the formula: q1.. Q is the number of samples of the original vibration signal for one position;
n is the number of decomposed IMFs;
s33, according to the steps S31-S32, the corresponding eigenvectors M of IMF components of all vibration signals of the same vibration signal form an energy matrix E0；
S34, obtaining energy matrixes E corresponding to the 1-bit axle box transverse vibration signal, the 1-bit framework transverse vibration signal and the 4-bit upper framework transverse vibration signal according to the steps S31-S331Energy matrix E2And energy matrix E3；
S35, converting the energy matrix E1Energy matrix E2Energy matrix E3And combining to obtain an energy matrix E corresponding to the vibration signals.
Further, the step S4 is specifically:
s41, decomposing the time sequence T into K time intervals TKFurther generate corresponding K covariance matrixes
Figure BDA0002039410560000052
The time sequence T is a corresponding time node when the original vibration signal is extracted;
s42, determining a diagonalized matrix corresponding to the K covariance matrices
Figure BDA0002039410560000053
Therein, diagonalizing the matrix
Figure BDA0002039410560000054
Comprises the following steps:
Figure BDA0002039410560000055
wherein K is 2.., K;
Figure BDA0002039410560000056
is a first covariance matrix
Figure BDA0002039410560000057
A diagonalized matrix of
Figure BDA0002039410560000058
V is the first covariance matrix
Figure BDA0002039410560000059
A first covariance matrix of
Figure BDA00020394105600000510
The superscript H of the feature vector of (1) is a conjugate transpose operator;
s43, according to the diagonalized matrix
Figure BDA0002039410560000061
Determining a unitary matrix U such that the following value is maximized;
Figure BDA0002039410560000062
wherein: the function diag represents a diagonal element;
p is the dimensionality after dimensionality reduction;
b and d denote the b-th row and d-th column of the matrix.
S44, setting a threshold epsilon, and aligning the unitary matrix U and the covariance matrix according to the rotation matrix
Figure BDA0002039410560000063
Performing iterative update until the updated unitary matrix U and covariance matrix
Figure BDA0002039410560000064
When the values of all the non-diagonal elements in (1) are less than the threshold epsilon, the final unitary matrix is obtained
Figure BDA0002039410560000065
Complete the K covariance matrices
Figure BDA0002039410560000066
Joint diagonalization of (a);
k ═ 2.. K;
the rotation matrix G (i, j, θ) is:
Figure BDA0002039410560000067
in the formula: i and j are respectively the ith row and the jth column of the rotation matrix;
theta is an intermediate calculation parameter;
s45 according to the final unitary matrix
Figure BDA0002039410560000068
Determining a conversion matrix A;
wherein, the transformation matrix A is:
Figure BDA0002039410560000069
in the formula: superscript # is the pseudo-inverse transform operator;
s46, obtaining a decomposed fusion characteristic matrix Z according to the conversion matrix A;
Z＝A·E。
further, in step S44, the iterative formula for the unitary matrix U is:
U←UG(1,2,θ)
in the formula: g (1,2, theta) is the element value of the 1 st row and 2 nd column of the rotation matrix G (i, j, theta);
for covariance matrix
Figure BDA0002039410560000071
The iterative formula of (a) is:
Figure BDA0002039410560000072
further, the small snaking operation characteristic of the high-speed train in the step S4 includes a normal operation state characteristic, a small convergent state characteristic, a small divergent state characteristic, and a large snaking state characteristic.
The invention has the beneficial effects that: the rapid extraction method for the small snaking running characteristic of the high-speed train analyzes the running signal of the high-speed train, successfully separates the normal running state, the small convergence state and the large snaking running state of the high-speed train, rapidly and accurately detects the running state of the high-speed train, and extracts the running state characteristic of the high-speed train by using the LSSVM learned by a machine, thereby predicting whether the high-speed train has the large snaking instability state or not and being beneficial to improving the accuracy and timeliness of train monitoring.
Drawings
FIG. 1 is a schematic diagram of a small snake-shaped operation characteristic of a high-speed train in the background art of the invention.
FIG. 2 is a flow chart of a rapid extraction method of small snaking running characteristics of a high-speed train in the invention.
FIG. 3 is a time domain diagram of simulation signals x (t), y (t), z (t) according to an embodiment of the present invention.
Fig. 4 is a simulation signal JADE method clustering effect graph.
Fig. 5 is a graph of the train speed profile and the lateral acceleration of the bogie frame in the embodiment of the present invention.
FIG. 6 is a decomposition result diagram of the normal train operation signal EEMD according to the embodiment of the present invention.
FIG. 7 is a decomposition result diagram of the train small amplitude convergence signal EEMD in the embodiment of the present invention.
FIG. 8 is a decomposition result diagram of the EEMD of the small divergent train signal in the embodiment of the present invention.
FIG. 9 is a decomposition result diagram of a large train serpentine signal EEMD according to an embodiment of the present invention.
Fig. 10 is a JADE feature extraction diagram of position 1 in the embodiment provided by the invention.
Fig. 11 is a JADE feature extraction diagram of position 2 in the embodiment of the present invention.
Fig. 12 is a JADE feature extraction diagram of position 3 in an embodiment provided by the present invention.
Fig. 13 is a multi-position fused JADE feature extraction diagram in an embodiment provided by the invention.
Detailed Description
The following description of the embodiments of the present invention is provided to facilitate the understanding of the present invention by those skilled in the art, but it should be understood that the present invention is not limited to the scope of the embodiments, and it will be apparent to those skilled in the art that various changes may be made without departing from the spirit and scope of the invention as defined and defined in the appended claims, and all matters produced by the invention using the inventive concept are protected.
As shown in fig. 2, a method for rapidly extracting small snaking operation characteristics of a high-speed train includes the following steps:
s1, extracting an original vibration signal when the high-speed train runs, and preprocessing the original vibration signal;
the original vibration signals comprise 1-bit axle box transverse vibration signal, 1-bit framework transverse vibration signal and 4-bit upper framework transverse vibration signal;
the method for preprocessing the original vibration signal specifically comprises the following steps: and sequentially carrying out resampling processing, band-pass filtering processing, zero-averaging processing, trend item removal and other processing on the extracted original vibration signal.
S2, decomposing the preprocessed vibration signal into a plurality of IMF components through EEMD;
s3, constructing an energy matrix corresponding to the vibration signal according to the IMF component;
and S4, processing the energy matrix through a non-stationary state JADE method to obtain a fusion characteristic data matrix, namely the small snaking running characteristic of the high-speed train.
The small snaking running characteristic of the high-speed train comprises a normal running state characteristic, a small convergence state characteristic, a small divergence state characteristic and a large snaking state characteristic.
Wherein, step S2 specifically includes:
s21, adding a Gaussian white noise sequence omega (t) with the same amplitude to the preprocessed vibration signal a (t) to obtain an overall signal a' (t);
wherein the overall signal a' (t) is:
a'(t)＝a(t)+ω(t)
s22, decomposing the overall signal a' (t) added with white noise each time according to an EMD method, and determining each order component ci：
The decomposition formula of the overall signal a' (t) is as follows:
Figure BDA0002039410560000091
in the formula: c. CiIs IMF component of each order;
i represents the ith IMF component;
r is a residual term;
n is the number of decomposed IMFs.
S23, adding different white noise sequences omega with the same amplitude each timej(t), repeating steps S21-S22 to obtain:
Figure BDA0002039410560000092
in the formula: a'j(t) is the overall signal after the j-th added Gaussian white noise sequence;
n is the number of times of adding Gaussian white noise sequence;
ωj(t) is the j-th added Gaussian white noise sequence;
cijthe ith IMF component is decomposed when the white Gaussian noise sequence is added for the jth time;
rjand residual values are decomposed when the white Gaussian noise sequence is added for the jth time.
S24, according to the zero-mean principle of Gaussian white noise frequency, eliminating the influence of the Gaussian white noise on the IMF component, and obtaining the IMF component corresponding to the preprocessed original vibration signal as follows:
Figure BDA0002039410560000101
in the formula: c. Ci(t) is an IMF component corresponding to the ith original vibration signal;
cij(t) is the i-th IMF component obtained by EEMD decomposition of the original vibration signal.
When the high-speed train runs, the signal is generally a non-stationary signal, and the result is not ideal when the conventional IMF energy-based method embodies the hidden features in the signal, but the present invention uses the IMF energy moment method to extract the energy features of the high-speed train, so the step S3 is specifically:
s31, forming an IMF energy moment according to the IMF components as follows:
Figure BDA0002039410560000102
in the formula: p is the total number of sampling points;
m is a sampling point;
cithe IMF component corresponding to the selected vibration signal;
Δ t is the sampling period.
S32, according to the IMF energy moment, constructing a feature vector M of the normalized IMF energy moment as follows:
Figure BDA0002039410560000103
in the formula: q is 1,. Q is the number of samples of the original vibration signal for one position;
n is the number of the decomposed IMF components;
s33, according to the steps S31-S32, corresponding eigenvectors M of IMF components of all vibration signals of the same vibration signal form an energy matrix Ee0；
S34, obtaining energy matrixes E corresponding to the 1-bit axle box transverse vibration signal, the 1-bit framework transverse vibration signal and the 4-bit upper framework transverse vibration signal according to the steps S31-S331Energy matrix E2And energy matrix E3；
S35, converting the energy matrix E1Energy matrix E2And energy matrix E3And combining to obtain an energy matrix E corresponding to the vibration signals.
The feature vector of the original vibration signal sample at three positions is [ M1；M2；...；MQ]I.e. an energy matrix E can be formed1、E2And E3The dimension of the energy matrix is high and the features are not obvious. In order to better extract the vibration signal characteristics of the high-speed train in consideration of the non-stationarity of the high-speed train in operation, the method uses a JADE algorithm in a non-stationary state to extract the fusion characteristics of the energy matrixes at three positions.
Considering the non-stationarity of the high-speed train signal, in order to use the JADE method in the non-stationary environment, the invention decomposes the whole time sequence T into K time intervals T1,...,TKAnd then K covariance matrices can be generated
Figure BDA0002039410560000111
The method of the invention jointly diagonalizes the K covariance matrices to find aA unitary matrix U which can diagonalize K covariance matrixes, so as to reduce the dimension of an energy matrix E and extract dimension reduction characteristic data; therefore, step S4 specifically includes:
s41, decomposing the time sequence T into K time intervals TKFurther generate corresponding K covariance matrixes
Figure BDA0002039410560000112
The time sequence T is a corresponding time node when the original vibration signal is extracted;
for a time interval T1,...,TKThe covariance matrix is:
Figure BDA0002039410560000113
wherein: e (-) is the desired operator;
Etfor a time T ∈ TkAn inner energy matrix E;
s42, determining K covariance matrixes
Figure BDA0002039410560000114
Corresponding diagonalized matrix
Figure BDA0002039410560000115
It is desirable to diagonalize multiple matrices simultaneously
Figure BDA0002039410560000116
First, the first matrix needs to be diagonalized, and then the remaining K-1 matrices need to be diagonalized, wherein the matrix W is a covariance matrix
Figure BDA0002039410560000117
A diagonalized matrix of;
wherein the diagonalized matrix of the first covariance matrix is:
Figure BDA0002039410560000121
where V is the first covariance matrix
Figure BDA0002039410560000122
A first covariance matrix of
Figure BDA0002039410560000123
The superscript H of the feature vector of (1) is a conjugate transpose operator;
thus, diagonalizing the matrix
Figure BDA0002039410560000124
Comprises the following steps:
Figure BDA0002039410560000125
wherein K is 2.., K;
Figure BDA0002039410560000126
is a first covariance matrix
Figure BDA0002039410560000127
A diagonalized matrix of;
s43, according to the diagonalized matrix
Figure BDA0002039410560000128
Determining a unitary matrix U such that the following value is maximized;
the problem of jointly diagonalizing the matrix in this step can be transformed to find a unitary matrix U such that the following equation is minimal:
Figure BDA0002039410560000129
in the formula: the function off represents the sum of the squares of the off-diagonal elements;
b. d denotes the b-th row and d-th column of the data matrix, respectively.
Since the sum of squares of the above equation remains the same when multiplied by the orthogonal matrix, we can equivalently maximize the sum of squares of the diagonal elements, i.e. the problem turns to find a unitary matrix U that maximizes the following equation, so there is:
Figure BDA00020394105600001210
wherein: the function diag represents a diagonal element;
p is the dimensionality after dimensionality reduction;
b and d are the b-th row and d-th column of the matrix.
S44, setting a threshold epsilon, and aligning the unitary matrix U and the covariance matrix according to the rotation matrix
Figure BDA00020394105600001211
Performing iterative update until the updated unitary matrix U and covariance matrix
Figure BDA00020394105600001212
When the values of all the non-diagonal elements in (1) are less than the threshold epsilon, the final unitary matrix is obtained
Figure BDA0002039410560000131
Complete the K covariance matrices
Figure BDA0002039410560000132
Joint diagonalization of (a);
the method converts a plurality of covariance matrixes into diagonal matrixes by using Jacobi rotation, and the basic idea of the Jacobi method is to zero a pair of non-zero non-diagonal elements in a target matrix and reduce the square sum of the non-diagonal elements through one-time orthogonal transformation. The above process is repeated to make the square sum of the non-diagonal elements of the transformed matrix tend to zero, so that the matrix is approximated to a diagonal matrix to obtain all eigenvalues and eigenvectors.
The method comprises the following steps: substituting the unit matrix I as the initial value of the matrix U, and the value of theta is derived from the matrix
Figure BDA0002039410560000133
Of (2) element(s)
Figure BDA0002039410560000134
The calculation results in that,
Figure BDA0002039410560000135
converting the matrix set to diagonal form using the following rotation matrix G (i, j, θ);
Figure BDA0002039410560000136
in the formula: i and j are respectively the ith row and the jth column of the rotation matrix;
theta is an intermediate calculation parameter.
The iterative formula of the unitary matrix U is as follows:
U←UG(1,2,θ)
wherein G (1,2, theta) is the element value of the 1 st row and 2 nd column of the rotation matrix G (i, j, theta);
for covariance matrix
Figure BDA0002039410560000137
The iterative formula of (a) is:
Figure BDA0002039410560000138
s45 according to the final unitary matrix
Figure BDA0002039410560000139
Determining a conversion matrix A;
wherein, the transformation matrix A is:
Figure BDA0002039410560000141
in the formula, the superscript # is a pseudo-inverse transformation operator;
s46, obtaining a decomposed fusion characteristic matrix Z according to the conversion matrix A;
Z＝A·E。
in the invention, a least square method support vector machine (LSSVM) is used for verifying the extracted fusion characteristic matrix, and the fusion characteristic matrix is input into the LSSVM for training and testing, so that the accuracy of characteristic data is tested, and the accuracy of prediction is enhanced; the LSSVM method has the core that the quadratic programming problem is replaced by solving a linear equation set, so that an insensitive loss function is avoided, the operation complexity is reduced, and the operation efficiency is improved.
In the process, the dimension reduction characteristic data matrix Z is taken as a sample, and because the first-dimension data in Z contains the most information, the first-dimension vector Z of the dimension reduction matrix Z is taken as: { zq,yq},q＝1,...,Q，y∈[-1,1]Z is dimension reduction characteristic data, y is a classification label, Q is the number of samples, and the function of the hyperplane classified by the least square support vector machine is as follows:
Figure BDA0002039410560000142
wherein:
Figure BDA0002039410560000143
is an adjustable weight vector;
g is an offset;
phi (z) is a non-linear function, and can deal with the problem of sample linear inseparability, and the basic idea is to map linearly inseparable data in an input space to a feature vector of a high-dimensional feature space, so that the data becomes separable.
The optimization problem of the least squares support vector machine can be expressed as:
Figure BDA0002039410560000144
in the formula: z is an input vector, namely dimension reduction characteristic data;
ξl>0 is a relaxation variable used for measuring the deviation degree;
gamma is a penalty factor;
superscript H is the conjugate transpose.
Lagrange multiplication is used on the above equation to obtain:
Figure BDA0002039410560000151
wherein: beta is Lagrangian factor;
and then respectively pair
Figure BDA0002039410560000152
The derivation can be found as follows:
Figure BDA0002039410560000153
obtained by the above formula:
Figure BDA0002039410560000154
wherein I is an identity matrix;
Ωql＝yqylφ(zl)Hφ(zq)＝yqylD(zq,zl) Q, l ═ 1.., η is the kernel matrix;
Figure BDA0002039410560000155
to satisfy the kernel function of Mercer's theorem, q, o denote the elements of sample z, and σ denotes the parameter that controls the kernel function width.
With the above equation, the LSSVM optimization problem can be converted to a solution problem of a linear system of equations, so the hyperplane function of the LSSVM classifier can be written as:
Figure BDA0002039410560000161
in the formula: (z) is a least squares support vector machine objective function;
sgn (·) is a step function;
βlis lagrange factor, l is 1lThe number of (2);
ylto classify labels, y e [ -1,1]；
D(z,zl) Is a kernel function that satisfies the Mercer theorem;
g is the bias.
Through the process, different motion states of the high-speed train are predicted and extracted in time, the high-speed train is prevented from snake-shaped instability, and the running safety of the train is guaranteed.
In a first embodiment of the invention, the proposed method was tested in simulation in order to verify its superiority. The analog signals x (t), y (t), z (t) are selected to analyze the feature extraction method, as shown in the following formula. The simulation signals x (t), y (t), z (t) are obtained by amplitude modulation addition of three original signals with frequencies of 10Hz, 20Hz and 30Hz, and Gaussian white noise with the signal-to-noise ratio of 10 is added to the three simulation signals on the basis of the three original signals. The sampling frequency was 500Hz and the resulting time domain plot of the three signals is shown in figure 3 below.
Processing the simulation signal according to the framework and performing dimensionality reduction processing to obtain low-dimensional features, and projecting the features in a two-dimensional space, as shown in fig. 4 below;
the respective feature heights of the simulation signals x (t), y (t), z (t) are gathered, only one point can be seen in a two-dimensional space, the feature arrangement positions of different simulation signals are greatly different, the features of different signals can be well distinguished, the recognition effect is excellent, and the correctness of the method is further verified.
In the second embodiment of the present invention, an example of extracting the small snaking operation state features of the high-speed train by using the method of the present invention is provided:
(1) acquiring experimental data: the experimental data are from scientific experiments of a certain type of motor train unit, the signal sampling frequency is 2500Hz, the train speed information is provided by a vehicle-mounted wireless GPS, and the acceleration information is provided by a plurality of sensors such as a 2-train 1-bit axle box transverse acceleration sensor, a 2-train 1-bit framework transverse acceleration sensor, a 2-train 4-bit axle box upper framework transverse acceleration sensor and the like. The line adopts the German Bo grid ballastless track technology, and the steel rail is a domestic steel rail with the fixed length of 100 m. The time-speed curve of the train and the time-acceleration curve of the bogie frame are shown in fig. 5(a), and the total time length is 2491 seconds. Through observation, the experimental train has repeated hunting instability when the speed is 330-350 km/h.
Raw data from the truck frame sensors is relatively coarse and requires some pre-processing. The frequency of the snaking motion of the high-speed train is 2-12.07 Hz, band-pass filtering of 2-12.07 Hz is carried out on the original data, so that a large number of interference signals in the original signals are effectively eliminated, accuracy and efficiency of feature identification are improved, resampling with the sampling frequency of 250Hz is carried out on the original data according to the Shannon sampling theorem, and then zero averaging, trend item elimination and smoothing processing are carried out on the data, so that smooth and accurate original data are obtained.
(2) EEMD decomposition: when the running speed of the high-speed train is 330-350 km/h, 10 samples in four states of normal running, small convergence, small divergence and large snaking are selected from the original data of one position, the total number of the samples is 40, and the total number of the samples in the original data of the three positions is 120. If the sample data is too long, the calculation amount and the operation time can be increased; and if the sample data is short, the running state of the train cannot be completely represented, and in combination with the fact that the sample length is 500 and the sampling time is 2 seconds, the method is selected.
The samples obtained from each original data are processed separately, EEMD decomposition is carried out on signals in four states of normal operation, small convergence, small divergence and large snaking of the high-speed train, IMF components obtained by EEMD decomposition in four operation states of normal operation, small convergence, small divergence and large snaking of the high-speed train are shown in figures 6 to 9, and the energy of the original signals is mainly concentrated in the first IMF components.
(3) IMF energy diagnosis is extracted from IMF decomposed by EEMD, converted into IMF energy moment, IMF energy characteristics are extracted, and energy matrix E of three positions can be obtained1、E2、E3Are combined to obtain an energy matrixE, since the energy feature dimension is 9, too high a dimension results in poor recognition of the energy feature.
(4) And fusing the vibration matrix signals of the high-speed train at different positions by using a non-stationary state JADE method, combining a plurality of diagonalized energy matrixes, and fusing the characteristics to extract the characteristics of the high-speed train in different running states. The data is reduced to three dimensions and the features are projected into a three-dimensional interface as shown in fig. 13. For comparison, fig. 10 to 12 are JADE feature extraction diagrams of the journal box position (position 1), the 1-position frame position (position 2), and the 4-position frame position (position 3), respectively. In addition, in order to highlight the rapidity of the method provided by the text, the running time comparison of a plurality of methods is shown in table 1, a CPU of a PC platform is an Intel Core i5-4460, a 12GB memory, and a video card is NVIDIA GeForce GT 720;
table 1: run-time comparison of multiple methods
Figure BDA0002039410560000181
(5) In order to verify the accuracy of feature extraction, the dimension reduction features are used as the input of the LSSVM, and the accuracy of feature extraction is verified. The data are divided into two types of training and testing, wherein each type has 20 groups of data (5 groups of normal operation, 5 groups of small convergence, 5 groups of small divergence and 5 groups of large snake movement). For comparison, features generated by other methods are also put into the LLSVM together, wherein one half of the features are used for training and the other half of the features are used for testing. Because the vibration data of the high-speed train at a plurality of positions are utilized, and the JADE method in a non-steady state is utilized to fuse and reduce the data at the plurality of positions, the recognition rate of the method provided by the invention is the highest and reaches 100 percent, and the accuracy and the superiority of the feature extraction method provided by the invention are also proved.
The invention has the beneficial effects that: the rapid extraction method for the small snaking running characteristic of the high-speed train analyzes the running signal of the high-speed train, successfully separates the normal running state, the small convergence state and the large snaking running state of the high-speed train, timely and accurately detects the running state of the high-speed train, and extracts the small snaking running characteristic of the high-speed train by using the LSSVM learned by a machine, thereby predicting whether the high-speed train has the large snaking instability state or not and being beneficial to improving the accuracy and the timeliness of train monitoring.

Claims (3)

1. A method for rapidly extracting small-amplitude snaking running characteristics of a high-speed train is characterized by comprising the following steps:
s1, extracting an original vibration signal when the high-speed train runs, and preprocessing the original vibration signal;
s2, decomposing the preprocessed vibration signal into a plurality of IMF components through EEMD;
s3, constructing an energy matrix corresponding to the vibration signal according to the IMF component;
s4, processing the energy matrix through a non-stationary state JADE method to obtain a fusion characteristic matrix, namely a small snaking running characteristic of the high-speed train;
the original vibration signals in the step S1 include a 1-bit axle box transverse vibration signal, a 1-bit framework transverse vibration signal and a 4-bit upper framework transverse vibration signal;
the method for preprocessing the original vibration signal specifically comprises the following steps: the extracted original vibration signal is subjected to resampling processing, band-pass filtering processing, zero-averaging processing and trend item removing processing in sequence;
the step S2 specifically includes:
s21, adding a Gaussian white noise sequence omega (t) with the same amplitude to the preprocessed vibration signal a (t) to obtain an overall signal a' (t);
wherein the overall signal a' (t) is:
a'(t)＝a(t)+ω(t)
s22, decomposing the overall signal a' (t) added with white noise each time according to an EMD method, and determining each order component ci：
The decomposition formula of the overall signal a' (t) is as follows:
Figure FDA0002790213230000011
in the formula: c. CiIs IMF component of each order;
i represents the ith IMF component;
r is a residual term;
n is the number of decomposed IMFs;
s23, adding different white noise sequences omega with the same amplitude each timej(t), repeating steps S21-S22 to obtain:
Figure FDA0002790213230000021
in the formula: a isj' (t) is the overall signal after the j-th added Gaussian white noise sequence;
n is the number of times of adding Gaussian white noise sequence;
ωj(t) is the j-th added Gaussian white noise sequence;
cijthe ith IMF component is decomposed when the white Gaussian noise sequence is added for the jth time;
rjresidual error values obtained by time division for j-th addition of the Gaussian white noise sequence;
s24, according to the zero-mean principle of Gaussian white noise frequency, eliminating the influence of the Gaussian white noise on the IMF component, and obtaining the IMF component corresponding to the preprocessed original vibration signal as follows:
Figure FDA0002790213230000022
in the formula (I); c. Ci(t) is an IMF component corresponding to the ith original vibration signal;
cij(t) an ith IMF component obtained by EEMD decomposition of the original vibration signal;
the step S3 specifically includes:
s31, forming an IMF energy moment according to the IMF components as follows:
Figure FDA0002790213230000023
in the formula: p is the total number of sampling points;
m is a sampling point;
cithe IMF component corresponding to the selected vibration signal;
Δ t is the sampling period;
s32, according to the IMF energy moment, constructing a feature vector M of the normalized IMF energy moment as follows:
Figure FDA0002790213230000031
in the formula: q1.. Q is the number of samples of the original vibration signal for one position;
n is the number of decomposed IMFs;
s33, according to the steps S31-S32, the corresponding eigenvectors M of IMF components of all vibration signals of the same vibration signal form an energy matrix E0；
S34, obtaining energy matrixes E corresponding to the 1-bit axle box transverse vibration signal, the 1-bit framework transverse vibration signal and the 4-bit upper framework transverse vibration signal according to the steps S31-S331Energy matrix E2Energy matrix E3；
S35, converting the energy matrix E1Energy matrix E2And energy matrix E3Combining to obtain an energy matrix E corresponding to the vibration signal;
the step S4 specifically includes:
s41, decomposing the time sequence T into K time intervals TKFurther generate corresponding K covariance matrixes
Figure FDA0002790213230000032
The time sequence T is a corresponding time node when the original vibration signal is extracted;
s42, determining a diagonalized matrix corresponding to the K covariance matrices
Figure FDA0002790213230000033
Therein, diagonalizing the matrix
Figure FDA0002790213230000034
Comprises the following steps:
Figure FDA0002790213230000035
in the formula: k2.., K;
Figure FDA0002790213230000036
is a first covariance matrix
Figure FDA0002790213230000037
A diagonalized matrix of
Figure FDA0002790213230000038
V is the first covariance matrix
Figure FDA0002790213230000039
A first covariance matrix of
Figure FDA00027902132300000310
The superscript H of the feature vector of (1) is a conjugate transpose operator;
s43, according to the diagonalized matrix
Figure FDA00027902132300000311
Determining a unitary matrix U such that the following value is maximized;
Figure FDA0002790213230000041
wherein: the function diag represents a diagonal element;
p is the dimensionality after dimensionality reduction;
b. d-row b and column d of the matrix;
s44, setting a threshold epsilon, and aligning the unitary matrix U and the covariance matrix according to the rotation matrix
Figure FDA0002790213230000042
Performing iterative update until the updated unitary matrix U and covariance matrix
Figure FDA0002790213230000043
When the values of all the non-diagonal elements in (1) are less than the threshold epsilon, the final unitary matrix is obtained
Figure FDA0002790213230000044
Complete the K covariance matrices
Figure FDA0002790213230000045
Joint diagonalization of (a);
k ═ 2.. K;
the rotation matrix G (i, j, θ) is:
Figure FDA0002790213230000046
in the formula: i and j are respectively the ith row and the jth column of the rotation matrix;
theta is an intermediate calculation parameter;
s45 according to the final unitary matrix
Figure FDA0002790213230000047
Determining a conversion matrix A;
wherein, the transformation matrix A is:
Figure FDA0002790213230000048
in the formula: superscript # is the pseudo-inverse transform operator;
s46, obtaining a decomposed fusion characteristic matrix Z according to the conversion matrix A;
Z＝A·E。
2. the method for rapidly extracting small snaking features of high-speed train according to claim 1, wherein in step S44, the iterative formula for the unitary matrix U is:
U←UG(1,2,θ)
in the formula: g (1,2, theta) is the element value of the 1 st row and 2 nd column of the rotation matrix G (i, j, theta);
for covariance matrix
Figure FDA0002790213230000051
The iterative formula of (a) is:
Figure FDA0002790213230000052
3. the method for rapidly extracting small snaking characteristics of high-speed trains according to claim 1, wherein the small snaking characteristics of high-speed trains in step S4 include normal operation status characteristics, small convergent status characteristics, small divergent status characteristics and large snaking status characteristics.
CN201910336833.8A 2019-04-25 2019-04-25 Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train Expired - Fee Related CN110084185B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910336833.8A CN110084185B (en) 2019-04-25 2019-04-25 Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910336833.8A CN110084185B (en) 2019-04-25 2019-04-25 Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train

Publications (2)

Publication Number Publication Date
CN110084185A CN110084185A (en) 2019-08-02
CN110084185B true CN110084185B (en) 2021-03-16

Family

ID=67416623

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910336833.8A Expired - Fee Related CN110084185B (en) 2019-04-25 2019-04-25 Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train

Country Status (1)

Country Link
CN (1) CN110084185B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112948981B (en) * 2021-04-08 2022-06-21 西南交通大学 Interval prediction method for small-amplitude snaking evolution trend of high-speed train
CN113624521B (en) * 2021-08-20 2024-04-02 厦门物之联智能科技有限公司 Train serpentine instability monitoring method and system based on axle box vibration
CN114997252B (en) * 2022-08-05 2022-10-25 西南交通大学 Vehicle-mounted detection method for wheel polygon based on inertia principle
CN115017965B (en) * 2022-08-09 2022-11-01 西南交通大学 Snake classification method based on HHT energy and maximum Lyapunov exponent

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101900708A (en) * 2010-08-18 2010-12-01 哈尔滨工业大学 Vibration and audio signal-based high-speed train track defect detecting method
CN102175768A (en) * 2011-02-22 2011-09-07 哈尔滨工业大学 Method and device for detecting defects and failures of high-speed rail based on vibration signals
CN103163505A (en) * 2013-01-31 2013-06-19 西安电子科技大学 Time-varying narrow-band interference suppression method based on joint approximate diagonalization of eigen-matrices (JADE)
CN104502126A (en) * 2014-12-28 2015-04-08 华东交通大学 Modal intervals-based high-speed train bogie fault diagnosis method
CN105699080A (en) * 2015-12-18 2016-06-22 华北电力大学（保定） Wind turbine generator set bearing fault feature extraction method based on vibration data
CN105890744A (en) * 2016-03-30 2016-08-24 中车青岛四方机车车辆股份有限公司 Method for determining snakelike instability of rail vehicle
CN107192553A (en) * 2017-06-28 2017-09-22 石家庄铁道大学 Gear-box combined failure diagnostic method based on blind source separating
CN107563403A (en) * 2017-07-17 2018-01-09 西南交通大学 A kind of recognition methods of bullet train operating condition
CN108254179A (en) * 2017-08-08 2018-07-06 常州路航轨道交通科技有限公司 A kind of bullet train wheel set bearing method for diagnosing faults based on MEEMD arrangement entropys

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4481190B2 (en) * 2005-02-17 2010-06-16 株式会社沖データ Image processing method, image processing apparatus, and image forming apparatus
US8044629B2 (en) * 2008-08-29 2011-10-25 Northern Illinois University Self-tuning vibration absorber
JP5858149B2 (en) * 2012-04-25 2016-02-10 トヨタ自動車株式会社 Meander determination device

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101900708A (en) * 2010-08-18 2010-12-01 哈尔滨工业大学 Vibration and audio signal-based high-speed train track defect detecting method
CN102175768A (en) * 2011-02-22 2011-09-07 哈尔滨工业大学 Method and device for detecting defects and failures of high-speed rail based on vibration signals
CN103163505A (en) * 2013-01-31 2013-06-19 西安电子科技大学 Time-varying narrow-band interference suppression method based on joint approximate diagonalization of eigen-matrices (JADE)
CN104502126A (en) * 2014-12-28 2015-04-08 华东交通大学 Modal intervals-based high-speed train bogie fault diagnosis method
CN105699080A (en) * 2015-12-18 2016-06-22 华北电力大学（保定） Wind turbine generator set bearing fault feature extraction method based on vibration data
CN105890744A (en) * 2016-03-30 2016-08-24 中车青岛四方机车车辆股份有限公司 Method for determining snakelike instability of rail vehicle
CN107192553A (en) * 2017-06-28 2017-09-22 石家庄铁道大学 Gear-box combined failure diagnostic method based on blind source separating
CN107563403A (en) * 2017-07-17 2018-01-09 西南交通大学 A kind of recognition methods of bullet train operating condition
CN108254179A (en) * 2017-08-08 2018-07-06 常州路航轨道交通科技有限公司 A kind of bullet train wheel set bearing method for diagnosing faults based on MEEMD arrangement entropys

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
A multi-sensor fusion framework for detecting small amplitude hunting of high-speed trains;Jing Ning 等;《Journal of Vibration and Control》;20180717;第24卷(第17期);摘要，图5 *
Blind Source Separation Based on Joint Diagonalization in R: The Packages JADE and BSSasymp;Jari Miettinen 等;《Journal of Statistical Software》;20170111;第76卷(第2期);第1-31页 *
EMD-ISOMAP高速列车小幅蛇行异常特征提取;崔万里 等;《中国测试》;20161231;第42卷(第12期);第105-110页 *
Feature fusion using kernel joint approximate diagonalization of eigen-matrices for rolling bearing fault identification;Yongbin Liu 等;《Journal of Sound and Vibration》;20160921;第389-401页 *
JADE 盲源分离算法在变压器振动信号监测中的应用;徐智超 等;《传感器与微系统》;20170920;第36卷(第9期);第157-160页 *
基于EEMD-SVD-LTSA的高速列车蛇行演变特征提取框架;冉伟 等;《电子测量技术》;20190308;第42卷(第5期);摘要，第1-2节，图1 *
高速列车转向架蛇行失稳的MEEMD-LSSVM预测模型;叶运广 等;《铁道学报》;20180115;第40卷(第1期);第38-43页 *

Also Published As

Publication number Publication date
CN110084185A (en) 2019-08-02

Similar Documents

Publication Publication Date Title
CN110084185B (en) Method for rapidly extracting small-amplitude snaking operation characteristics of high-speed train
CN102629298B (en) A kind of Rail Transit System security of operation appraisal procedure
CN105346561A (en) Rail turnout disease detection system based on operating vehicle and rail turnout disease detection method based on operating vehicle
CN103196682B (en) Based on the train suspension system fault separating method of the information fusion of D-S evidence theory
CN106021789B (en) Railway vehicle suspension system Fault Classification and system based on fuzzy intelligence
CN109532940A (en) The fastener of high speed non-fragment orbit loosens degree detecting method
CN110411766A (en) The snakelike unstability detection method of train bogie, device, system and storage medium
CN102622519A (en) Method for estimating safety domain of track irregularity amplitude
CN106096096A (en) Train suspension system failure analysis methods based on MPCA and system
CN114881076B (en) Rail corrugation identification method, device, equipment and medium based on support vector machine
Li et al. Driver fatigue detection using approximate entropic of steering wheel angle from real driving data
CN112381027B (en) Wheel polygon wave depth estimation method based on train axle box vertical acceleration signal
Hong et al. High-speed rail suspension system health monitoring using multi-location vibration data
CN114997252B (en) Vehicle-mounted detection method for wheel polygon based on inertia principle
CN114861741B (en) Snake state identification method based on wheel set transverse displacement
CN114997218A (en) Recognition and detection method for polygonal abrasion of wheels of railway vehicle
CN107563403B (en) Working condition identification method for high-speed train operation
Zhang et al. Fault diagnosis on train brake system based on multi-dimensional feature fusion and GBDT enhanced classification
CN107688820A (en) A kind of Elevator Fault Diagnosis method based on BCSA Support Vector Machines Optimizeds
CN116467570B (en) Heavy-duty locomotive coupler swing angle quantitative identification method based on data driving
Zhang et al. Multi-sensor graph transfer network for health assessment of high-speed rail suspension systems
CN108509999A (en) It is a kind of indignation drive detection and safe early warning method
CN113627358A (en) Multi-feature fusion turnout intelligent fault diagnosis method, system and equipment
CN112308824B (en) Curve radius classification identification method and device based on track geometric detection data
Qin et al. Bogie fault identification based on EEMD information entropy and manifold learning

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210316