CN112326241B - Nuclear power main pump bearing fault early warning method based on fusion degradation index - Google Patents
Nuclear power main pump bearing fault early warning method based on fusion degradation index Download PDFInfo
- Publication number
- CN112326241B CN112326241B CN202010958398.5A CN202010958398A CN112326241B CN 112326241 B CN112326241 B CN 112326241B CN 202010958398 A CN202010958398 A CN 202010958398A CN 112326241 B CN112326241 B CN 112326241B
- Authority
- CN
- China
- Prior art keywords
- sample
- degradation
- data
- vibration acceleration
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M13/00—Testing of machine parts
- G01M13/04—Bearings
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F04—POSITIVE - DISPLACEMENT MACHINES FOR LIQUIDS; PUMPS FOR LIQUIDS OR ELASTIC FLUIDS
- F04B—POSITIVE-DISPLACEMENT MACHINES FOR LIQUIDS; PUMPS
- F04B49/00—Control, e.g. of pump delivery, or pump pressure of, or safety measures for, machines, pumps, or pumping installations, not otherwise provided for, or of interest apart from, groups F04B1/00 - F04B47/00
- F04B49/10—Other safety measures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M13/00—Testing of machine parts
- G01M13/04—Bearings
- G01M13/045—Acoustic or vibration analysis
Abstract
The invention discloses a nuclear power main pump bearing fault early warning method based on fusion degradation indexes, which comprises the following steps: the vibration acceleration sensor is arranged on the shell of the bearing bushing of the nuclear power main pump so as to acquire vibration acceleration data; cleaning missing points, blank points and outliers in the vibration acceleration data; removing amplitude modulation phenomenon data caused by variable rotation speed from vibration acceleration data to obtain quantized vibration acceleration data; extracting high-dimensional degradation features from the quantized vibration acceleration data; extracting a fusion degradation index for monitoring the fault state of the nuclear power main pump bearing from the high-dimensional degradation characteristic; thresholds for different degradation phases are determined based on the fused degradation indicators.
Description
Technical Field
The invention belongs to the technical field of nuclear power main pump bearing fault detection, and particularly relates to a nuclear power main pump bearing fault early warning method based on fusion degradation indexes.
Background
Nuclear power plant reactor coolant pumps, commonly referred to as primary pumps, are used to drive the circulating flow of coolant within the reactor coolant system, continuously transferring heat generated in the core to the steam generator secondary side feedwater while cooling the core, preventing fuel elements from burning or burning out, are key devices for the reactor coolant system and pressure boundaries. As a main pump of a heart of a nuclear power plant, the running state of the main pump directly relates to the benefit and nuclear safety of the whole nuclear power plant. According to incomplete statistics, since the 80 s of the 20 th century, nuclear power plants around the world have had over a hundred nuclear power plant shutdown events caused by main pump faults, resulting in significant economic losses. Therefore, the main pump is monitored in real time, and the running state of the main pump is predicted in advance, which is very necessary-!
The vibration performance and the running state of the bearing serving as a main pump core component directly reflect the running state of the main pump. However, at present, the monitoring and diagnosis of the nuclear power main pump bearing are very few, and a mode of temperature overrun alarm or vibration root mean square value threshold alarm is generally adopted. These approaches are relatively disfavored and hardly perform the function of early warning, and even if they can be advanced, the advanced time is within several hours, which is insufficient for maintenance of the nuclear power equipment. Maintenance of the main pump is mainly carried out on 'post maintenance' and 'timing maintenance', and the requirement of the intelligent nuclear power industry is difficult to meet.
The problems are solved to a certain extent by the aid of the post maintenance and the timing maintenance, but huge economic loss caused by unexpected shutdown and indirect even unexpected safety accidents cannot be controlled; the temperature alarm is based on indirect monitoring, because after equipment fails, vibration abnormality is firstly shown, and the temperature rise is caused after the abnormality is serious. Based on the vibration root mean square value threshold value alarm belongs to shallow analysis, only if the fault is serious, the fault can be found. The maintenance and early warning modes are difficult to predict faults in advance, avoid accidents and can not serve intelligent operation and maintenance of the nuclear power industry.
The above information disclosed in the background section is only for enhancement of understanding of the background of the invention and therefore may contain information that does not form the prior art that is already known in the country to a person of ordinary skill in the art.
Disclosure of Invention
Aiming at the problems existing in the prior art, the invention provides a nuclear power main pump bearing fault early warning method based on fusion degradation indexes, which comprises the steps of cleaning vibration acceleration data, extracting high-dimensional degradation characteristics, constructing the fusion degradation indexes, obtaining thresholds of different degradation states after exponential fitting, realizing the degradation state monitoring of a nuclear power main pump bearing based on the obtained thresholds, effectively early warning the occurrence and evolution of faults, and improving the reliability of the nuclear power main pump bearing fault monitoring.
The invention aims at realizing the technical scheme that the nuclear power main pump bearing fault early warning method based on the fusion degradation index comprises the following steps:
in the first step, a vibration acceleration sensor is arranged on a nuclear power main pump bearing bushing shell to acquire vibration acceleration data;
in the second step, cleaning missing points, blank points and outliers in the vibration acceleration data;
in the third step, eliminating the amplitude modulation phenomenon data caused by the variable rotation speed from the vibration acceleration data to obtain quantized vibration acceleration data;
In the fourth step, high-dimensional degradation features are extracted from the quantized vibration acceleration data;
in the fifth step, fusion degradation indexes for monitoring the fault state of the nuclear power main pump bearing are extracted from the high-dimensional degradation characteristics;
in a sixth step, thresholds for different degradation phases are determined based on the fused degradation indicators.
In the second step, the missing points in the vibration acceleration data are cleaned through polynomial fitting, and the sample x with the serial numbers of 1,2, … and n is arranged in front of the missing points 1 ，x 2 ，…，x n Fitting on the basis of the data, and solving p (t) i )＝a 0 +a 1 t i +a 2 t i 2 +…+a n t i n And is such thatAnd is composed of p (t) i ) Substitution of deletion points, wherein t i To fit the data samples, p (t i ) At t i Fitting data corresponding to points, a 0 ，a 1 ，…a n Is the fitting coefficient.
In the second step, when the blank points in the vibration acceleration data are cleaned, judging whether 2 or more continuous points with the amplitude of zero exist, and x i ＝x i+1 = … =0, and if present, successive sample points of zero amplitude are removed from the vibration acceleration data.
In the second step, normal data x is obtained when outliers in the vibration acceleration data are cleaned i The range is as follows:mean->The estimation method comprises the following steps: / >The standard deviation sigma estimation method comprises the following steps: />The outlier judging mode is as follows: />
Sample x is determined sample by sample i Whether or not it is in an intervalIn addition, beta is determined i Is less than or equal to 3, or beta i > 3, if sample x i Corresponding beta i Sample x is less than or equal to 3 i Is the normal point, otherwise if sample x i Corresponding beta i > 3, sample x i As outliers, a cleaning process is required.
In the method, the average value of a plurality of samples adjacent to the outlier is taken to replace the value, and the expression is:
wherein N is an outlier sample x i The total number of adjacent sample points is taken, i is the sample number.
In the method, the data of amplitude modulation phenomenon caused by variable rotation speed is removed from vibration acceleration data,
for vibration data x containing variable rotation speed working condition i Taking his Hilbert envelope e i ，e i ＝|x i +jH(x i ) I, wherein H (x i ) Representing de-vibration data x i Is a hilbert transform of (c); x is x i +jH(x i ) Plural, |x i +jH(x i ) The I represents the complex number x i +jH(x i ) Is provided with a die for the mold,
for envelope data e i The magnitudes of the first m magnitudes are ordered in descending order, and the average μ of the first m magnitudes is calculated:wherein max (e i ) m Representing envelope data e i The sum of the first m maximum amplitude values, the quantized vibration acceleration data x 'after eliminating the influence of the rotation speed factor' i The method is characterized by comprising the following steps: / >
In the method, in the fourth step, the high-dimensional degradation feature is extracted from the quantized vibration acceleration data, and the high-dimensional degradation feature includes:
kurtosis characteristics:wherein->
Waveform index:
kurtosis index:wherein->
Margin index:
frequency domain feature 3:wherein->
Frequency domain feature 4:wherein->
Frequency domain features 12:wherein->
Wherein i is a vibration acceleration data sequence number, and n is the total number of data samples; u (u) i Representing a quantized vibration acceleration data sample; y is j Representing quantized vibration acceleration data u i Frequency domain data obtained after FFT conversion; j is quantized vibration acceleration data u i Frequency domain data y of (2) j N represents quantized vibration acceleration data u i Is a sample total number of frequency domain data; f (f) j Representing frequency domain data y j Frequency components in the frequency spectrum.
In the method, in five steps, the fusion degradation index for monitoring the fault state of the nuclear power main pump bearing is extracted from the high-dimensional degradation characteristic,
solving cosine similarity matrix CS and sample cosine similarity mean vector among all samples in high-dimensional degradation feature set X
According to the cosine similarity matrix CS and the sample cosine similarity mean vectorCalculating a similarity matrix S;
calculating to obtain an initial neighborhood parameter K according to the similarity matrix S;
Obtaining the Euclidean distance D between all samples in the high-dimensional degradation feature set X e Distance D from geodesic g ；
Averaging the local density P and global density across all samples in the high dimensional degradation feature set X
Calculating the ratio of the sum of the local Euclidean distances and the sum of the local geodesic distances to obtain the local popular curvature Q and the global popular curvatureMeanwhile, calculating the ratio of the Euclidean distance to the geodesic distance to obtain D';
adjusting the initial neighborhood parameter K based on the neighborhood parameter adjusting method to obtain an adjusted neighborhood parameter K c At the same time, the Euclidean distance matrix D is updated by the ratio D' of Euclidean distance to geodesic distance e Obtaining updated Euclidean distance matrix D e ′；
Based on the adjusted neighborhood parameter K c And updated Euclidean distance matrix D e ' extracting low-dimensional features obtained after high-dimensional degradation feature mapping by a local linear embedding method, and solving to obtain an objective function by calculating a Lagrangian multiplier method: MY (MY) T ＝λ′Y T Solving the formula to obtain a feature vector Y and a feature value lambda 'of the feature matrix M'
Constraint of->
M＝(I i -w i )(I i -w i ) T ，
Wherein W represents a set of all sample weight vectors; x is x i And x j Two samples in the high-dimensional degradation feature X, respectively; w (W) i For sample x i Is a weight vector of (2); w (w) ij Is w i Sample points in (a); i i A D-dimensional vector representing all 1 s; m is a characteristic matrix, wherein M is a characteristic matrix,
y is a matrix formed by feature vectors of the feature matrix M, and the feature value of the matrix M is decomposed to obtain a low-dimensional representation Y of the high-dimensional degradation feature X: X-Y= [ Y ] 2 ，y 3 ，...，y d+1 ]Wherein y is i I=2, 3,..d+1, representing the eigenvector corresponding to the first d non-zero eigenvalues of matrix M, will be eigenvalue y 2 As a fusion degradation index for the fault condition monitoring index of the nuclear power main pump bearing.
In the method, the cosine similarity matrix CS is a symmetric matrix, and the solving formula is as follows:
wherein z is a feature dimension sequence number, D is the dimension of the high-dimensional degradation feature, i and j are sample sequence numbers; CS (i, j) is a sample in the cosine similarity matrix CS; x is x i And x j Two samples in the high-dimensional degradation feature X,and->Respectively sample x i And sample x j In (c) the sample points in (c),
geodesic distance D between all samples of high-dimensional degradation features of bearing g The iterative expression of (2) is:
d g (i，j)＝min{d g (i，j)，d g (i，m)+d g (m，j)}
wherein min is the minimum value; i, j, m are sample sequence numbers, and i, j, m are mutually unequal; n is the total number of samples; d, d g (i, j) is the geodesic distance D g Is a sample of (b).
In the method, in the sixth step, determining the threshold values of the different degradation phases based on the fusion degradation indicator includes exponentially fitting a fusion degradation indicator curve to determine the threshold value intervals of the different degradation states: f (t) =λ+βexp (ηt), wherein: λ, β, η are parameters of an exponential function, which are determined by fitting to a fused degradation indicator.
Advantageous effects
According to the invention, blank data and outlier data in complex vibration acceleration data are cleaned by a data cleaning method, and amplitude modulation phenomenon of the vibration acceleration data generated by a variable speed working condition is cleaned, so that the reliability of the vibration acceleration data to be analyzed is ensured; based on the cleaned vibration acceleration data, a nuclear power main pump bearing fusion degradation index is extracted through the improved local linear embedding method, thresholds of different degradation states are determined on the basis, and the degradation states of the bearing can be intuitively depicted in real time based on different threshold intervals determined by the extracted thresholds, so that the monitoring reliability of the nuclear power main pump bearing is effectively improved, and early warning is timely carried out on the occurrence of faults.
Drawings
Various advantages and benefits of this invention will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments. The drawings are only for purposes of illustrating the preferred embodiments and are not to be construed as limiting the invention. It is evident that the figures described below are only some embodiments of the invention, from which other figures can be obtained without inventive effort for a person skilled in the art. Also, like reference numerals are used to designate like parts throughout the figures.
In the drawings:
FIG. 1 is a flow diagram of a nuclear power main pump bearing fault early warning method based on fusion degradation indexes;
FIG. 2 is an original vibration acceleration diagram of a nuclear power main pump bearing based on a nuclear power main pump bearing fault early warning method based on fusion degradation indexes;
FIG. 3 is a graph of vibration acceleration after cleaning of a nuclear power main pump bearing fault early warning method based on fusion degradation indexes;
FIG. 4 is a square root amplitude chart of 5 rotating speed working conditions of a nuclear power main pump bearing based on a nuclear power main pump bearing fault early warning method based on fusion degradation indexes;
fig. 5 is a graph of a total life trend curve, an exponential fit and a threshold interval of fusion degradation indexes of a nuclear power main pump bearing based on a nuclear power main pump bearing fault early warning method based on fusion degradation indexes.
The invention is further explained below with reference to the drawings and examples.
Detailed Description
Specific embodiments of the present invention will be described in more detail below with reference to fig. 1 to 5. While specific embodiments of the invention are shown in the drawings, it should be understood that the invention may be embodied in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.
It should be noted that certain terms are used throughout the description and claims to refer to particular components. Those of skill in the art will understand that a person may refer to the same component by different names. The specification and claims do not identify differences in terms of components, but rather differences in terms of the functionality of the components. As used throughout the specification and claims, the terms "include" and "comprise" are used in an open-ended fashion, and thus should be interpreted to mean "include, but not limited to. The description hereinafter sets forth a preferred embodiment for practicing the invention, but is not intended to limit the scope of the invention, as the description proceeds with reference to the general principles of the description. The scope of the invention is defined by the appended claims.
For the purpose of facilitating an understanding of the embodiments of the present invention, reference will now be made to the drawings, by way of example, and specific examples of which are illustrated in the accompanying drawings.
For better understanding, fig. 1 is a flowchart of a method for early warning a failure of a nuclear power main pump bearing based on a fusion degradation index, and as shown in fig. 1, the method for early warning the failure of the nuclear power main pump bearing based on the fusion degradation index comprises the following steps:
In the first step S100, a vibration acceleration sensor is provided on the housing of the bearing bush of the nuclear power main pump to collect vibration acceleration data;
in a second step S200, the missing points, blank points and outliers in the vibration acceleration data are cleaned;
in the third step S300, the data of the amplitude modulation phenomenon caused by the variable rotation speed is removed from the vibration acceleration data to obtain quantized vibration acceleration data;
in a fourth step S400, extracting high-dimensional degradation features from the quantized vibration acceleration data;
in a fifth step S500, extracting a fusion degradation index for monitoring the failure state of the nuclear power main pump bearing from the high-dimensional degradation characteristic;
in a sixth step S600, thresholds for different degradation phases are determined based on the fused degradation indicators.
In a preferred embodiment of the method, in the second step S200, missing points in the vibration acceleration data are cleaned by polynomial fitting, and samples x with serial numbers 1,2, …, n before the missing points are cleaned 1 ，x 2 ，…，x n Fitting on the basis of the data, and solving p (t) i )＝a 0 +a 1 t i +a 2 t i 2 +…+a n t i n And make (1)And is composed of p (t) i ) Substitution of deletion points, wherein t i To fit the data samples, p (t i ) At t i Fitting data corresponding to points, a 0 ，a 1 ，…a n Is the fitting coefficient.
In a preferred embodiment of the method, in the second step S200, when the blank points in the vibration acceleration data are cleaned, it is determined whether there are 2 or more continuous points with zero amplitude, x i ＝x i+1 = … =0, and if present, successive sample points of zero amplitude are removed from the vibration acceleration data.
In a preferred embodiment of the method, in a second step S200, normal data x is obtained when cleaning outliers in the vibration acceleration data i The range is as follows:mean->The estimation method comprises the following steps: />The standard deviation sigma estimation method comprises the following steps: />The outlier judging mode is as follows: />
Sample x is determined sample by sample i Whether or not it is in an intervalIn addition, beta is determined i Is less than or equal to 3, or beta i > 3, if sample x i Corresponding beta i Sample x is less than or equal to 3 i Is the normal point, otherwise if sample x i Corresponding beta i > 3, sample x i For outliers, cleaning is requiredAnd (5) washing.
In a preferred embodiment of the method, the average value of a plurality of samples adjacent to the outlier is taken to replace the value, and the expression is:
wherein N is an outlier sample x i The total number of adjacent sample points is taken, i is the sample number.
In a preferred embodiment of the method, the step of eliminating the amplitude modulation phenomenon caused by the variable rotation speed from the vibration acceleration data includes,
for vibration data x containing variable rotation speed working condition i Taking his Hilbert envelope e i ，e i ＝|x i +jH(x i ) I, wherein H (x i ) Representing de-vibration data x i Is a hilbert transform of (c); x is x i +jH(x i ) Plural, |x i +jH(x i ) The I represents the complex number x i +jH(x i ) Is provided with a die for the mold,
for envelope data e i The magnitudes of the first m magnitudes are ordered in descending order, and the average μ of the first m magnitudes is calculated:wherein max (e i ) m Representing envelope data e i The sum of the first m maximum amplitude values, the quantized vibration acceleration data x 'after eliminating the influence of the rotation speed factor' i The method is characterized by comprising the following steps: />
In a preferred embodiment of the method, in the fourth step S400, the high-dimensional degradation feature is extracted from the quantized vibration acceleration data, where the high-dimensional degradation feature includes:
kurtosis characteristics:wherein->
Waveform index:
kurtosis index:wherein->
Margin index:
frequency domain feature 3:wherein->
Frequency domain feature 4:wherein->
Frequency domain features 12:wherein->
Wherein i is a vibration acceleration data sequence number, and n is the total number of data samples; u (u) i Representing a quantized vibration acceleration data sample; y is j Representing quantized vibration acceleration data u i Frequency domain data obtained after FFT conversion; j is the quantization vibrationDynamic acceleration data u i Frequency domain data y of (2) j N represents quantized vibration acceleration data u i Is a sample total number of frequency domain data; f (f) j Representing frequency domain data y j Frequency components in the frequency spectrum.
In a preferred embodiment of the method, in a fifth step S500, extracting the fusion degradation indicator for monitoring the failure state of the main pump bearing of the nuclear power from the high-dimensional degradation characteristic comprises,
Solving cosine similarity matrix CS and sample cosine similarity mean vector among all samples in high-dimensional degradation feature set X
According to the cosine similarity matrix CS and the sample cosine similarity mean vectorCalculating a similarity matrix S;
calculating to obtain an initial neighborhood parameter K according to the similarity matrix S;
obtaining the Euclidean distance D between all samples in the high-dimensional degradation feature set X e Distance D from geodesic g ；
Averaging the local density P and global density across all samples in the high dimensional degradation feature set X
Calculating the ratio of the sum of the local Euclidean distances and the sum of the local geodesic distances to obtain the local popular curvature Q and the global popular curvatureMeanwhile, calculating the ratio of the Euclidean distance to the geodesic distance to obtain D';
adjusting the initial neighborhood parameter K based on the neighborhood parameter adjusting method to obtain an adjusted neighborhood parameter K c At the same time, the Euclidean distance matrix D is updated by the ratio D' of Euclidean distance to geodesic distance e Obtaining updated Euclidean distanceOff matrix D e ′；
Based on the adjusted neighborhood parameter K c And updated Euclidean distance matrix D e Extracting low-dimensional features obtained after high-dimensional degradation feature mapping by a local linear embedding method, and solving to obtain an objective function by calculating a Lagrangian multiplier method: MY (MY) T ＝λ′Y T Solving the formula to obtain a feature vector Y and a feature value lambda 'of the feature matrix M'
Constraint of->
M＝(I i -w i )(I i -w i ) T ，
Wherein W represents a set of all sample weight vectors; x is x i And x j Two samples in the high-dimensional degradation feature X, respectively; w (w) i For sample x i Is a weight vector of (2); w (w) ij Is w i Sample points in (a); i i Expressed as 1 D-dimensional vectors of (c); m is a characteristic matrix, wherein M is a characteristic matrix,
y is a matrix formed by feature vectors of the feature matrix M, and the feature value of the matrix M is decomposed to obtain a low-dimensional representation Y of the high-dimensional degradation feature X: X-Y= [ Y ] 2 ，y 3 ，...，y d+1 ]Wherein y is i I=2, 3,..d+1, representing the eigenvector corresponding to the first d non-zero eigenvalues of matrix M, will be eigenvalue y 2 As a fusion degradation index for the fault condition monitoring index of the nuclear power main pump bearing.
In a preferred embodiment of the method, the cosine similarity matrix CS is a symmetric matrix, and the calculation formula is as follows:
wherein z is a feature dimension number, D isThe dimension of the high-dimension degradation characteristic, i and j are sample serial numbers; CS (i, j) is a sample in the cosine similarity matrix CS; x is x i And x j Two samples in the high-dimensional degradation feature X,and->Respectively sample x i And sample x j In (c) the sample points in (c),
geodesic distance D between all samples of high-dimensional degradation features of bearing g The iterative expression of (2) is:
d g (i，j)＝min{d g (i，j)，d g (i，m)+d g (m，j)}
wherein min is the minimum value; i, j, m are sample sequence numbers, and i, j, m are mutually unequal; n is the total number of samples; dg (i, j) is the geodesic distance D g Is a sample of (b).
In a preferred embodiment of the method, in a sixth step S600, determining the threshold values of the different degradation phases based on the fused degradation indicators comprises exponentially fitting the fused degradation indicator curve to determine the threshold intervals of the different degradation states: f (t) =λ+βexp (ηt), wherein: λ, β, η are parameters of an exponential function, which are determined by fitting to a fused degradation indicator.
In a preferred embodiment, the nuclear power main pump bearing fault early warning method based on the fusion degradation index comprises the following steps:
s100, acquiring vibration acceleration data by a vibration acceleration sensor positioned on a shell of a bearing bushing of a nuclear power main pump;
s200, cleaning missing points, blank points and outliers contained in the collected vibration acceleration data;
s300, eliminating amplitude modulation phenomenon caused by variable rotation speed from vibration acceleration data to obtain quantized vibration acceleration data;
s400, extracting high-dimensional degradation characteristics from quantized vibration acceleration data;
s500, extracting a fusion degradation index from the high-dimensional degradation characteristic, and using the fusion degradation index for monitoring the fault state of the nuclear power main pump bearing;
S600 determines thresholds for different degradation phases based on the fused degradation indicators.
In step S100, vibration acceleration data x of a nuclear power main pump bearing is obtained by using a vibration acceleration sensor i And (3) collecting, wherein i is a vibration acceleration data serial number.
In step S201, blank data generally appears in the running process and in the start-stop stage of the main pump, and is collected by the vibration acceleration sensor, and such data can be classified into two types of single-point missing data and continuous blank invalid data.
1 for single point missing data, fitting is performed on the basis of a plurality of samples before the missing data by a polynomial fitting method, and filling is performed, namely, sample x with the sequence number of 1,2, …, n 1 ，x 2 ，…，x n On the basis of the data, the following polynomial is solved and represented by p (t i ) Instead of the blank point:
p(t i )＝a 0 +a 1 t i +a 2 t i 2 +…+a n t i n and make (1)
Wherein t is i To fit the data sample sequence number, p (t i ) At t i Fitting data corresponding to the points.
2 for continuous blank invalid points, by judging whether there are 2 or more continuous points with zero amplitude, i.e. x i ＝x i+1 = … =0, and if present, successive sample points of zero amplitude are removed from the vibration acceleration data.
In step S202, due to the non-stationarity of the vibration data, the normal vibration data contains a large number of outliers, and the cleaning of the outliers is determined and cleaned by the following improved radon criterion, which comprises the following steps:
1 estimation of the Normal data amplitude Range
The present disclosure determines the range of normal data based on an improved method of the larida criterion. Because noise and outliers exist in the vibration data, the real mean and standard deviation of the vibration data of the bearing are unknown and need to be estimated by the following statistical method:
normal data x i The range is as follows:
mean value ofThe estimation method comprises the following steps: />
The standard deviation sigma estimation method comprises the following steps:
the judgment mode of the outlier is as follows:
judging sample x one by one i Whether or not it is in an intervalIn addition to, i.e. judge beta i Is less than or equal to 3, or beta i > 3, if sample x i Corresponding beta i Sample x is less than or equal to 3 i Is the normal point, otherwise if sample x i Corresponding beta i > 3, sample x i For outliers, a cleaning process is required, and in order to ensure the reliability of the replacement data as much as possible, the average value of a plurality of samples adjacent to the value is taken to replace the value, and the expression is as follows:
wherein N is an outlier sample x i The total number of adjacent sample points is taken, i is the sample number.
In step S300, for the vibration data amplitude modulation phenomenon caused by the variable rotation speed working condition, the cleaning method is as follows:
for vibration data x containing variable rotation speed working condition i Taking his Hilbert envelope e i The formula is shown as follows:
e i ＝|x i +jH(x i )|
Wherein H (x) i ) Representing de-vibration data x i Is a hilbert transform of (c); x is x i +jH(x i ) Plural, |x i +jH(x i ) The I represents the complex number x i +jH(x i ) Is a mold of (a).
For envelope data e i The magnitudes of the first m magnitudes are ordered in descending order, and the average μ of the first m magnitudes is calculated:
wherein max (e i ) m Representing envelope data e i The sum of the first m maximum amplitude values is averaged over a plurality of amplitude values to reduce disturbances caused by non-stationary factors in the vibration data.
Quantized vibration acceleration data x after eliminating influence of rotation speed factor i The method is characterized by comprising the following steps:
in step S400, high-dimensional degradation features are extracted. Vibration acceleration data in three directions can be acquired through a triaxial vibration acceleration sensor, and 15-dimensional high-dimensional degradation characteristics X [ X ] shown below can be extracted after cleaning in the steps S200-S300 1 ，X 2 ，X 3 ，…，X 15 ]：
| Reference numerals | Features (e.g. a character) | Reference numerals | Features (e.g. a character) |
| X 1 | x-axis kurtosis characteristics | X 9 | y-axis waveform index |
| X 2 | x-axis waveform index | X 10 | Y-axis kurtosis index |
| X 3 | x-axis kurtosis index | X 11 | y-axis frequency domain features 3 |
| X 4 | x-axis margin index | X 12 | y-axis frequency domain features 4 |
| X 5 | x-axis frequency domain features 3 | X 13 | z-axis waveform index |
| X 6 | x-axis frequency domain features 4 | X 14 | z-axis frequency domain features 3 |
| X 7 | x-axis frequency domain features 12 | X 15 | z-axis frequency domain features 4 |
| X 8 | y-axis kurtosis characteristics |
The features in the table above are shown below:
kurtosis characteristics: Wherein->
Waveform index:
kurtosis index:wherein->
Margin index:
frequency domain feature 3:wherein->
Frequency domain feature 4:wherein->
Frequency domain features 12:wherein->
In the above calculation expression, i is the vibration acceleration data sequence number, n is the total number of data samples; u (u) i Representing a vibration acceleration data sample; y is j Representing vibration acceleration data u i Frequency domain data obtained after FFT conversion; j is vibration acceleration data u i Frequency domain data y of (2) j N represents vibration acceleration data u i Is a sample total number of frequency domain data; f (f) j Representing frequency domain data y j Frequency components in the frequency spectrum.
In particular, when the vibration acceleration sensor is uniaxial, the degradation characteristic X can be extracted 1 ～X 7 A total of 7 degradation features as high-dimensional degradation features; when the vibration acceleration sensor is in two axial directions, the degradation characteristic X can be extracted 1 ～X 12 A total of 7 degradation features were used as high-dimensional degradation features.
In step S500, a fusion degradation index is extracted from the high-dimensional degradation feature, and the fusion degradation index is used for monitoring the bearing fault state of the wind driven generator transmission system, and the specific process is as follows:
s510, initial neighborhood parameter setting, the general steps of which are as follows:
s511, solving the cosine similarity matrix CS and the sample cosine similarity mean vector among all samples in the high-dimensional degradation feature set X
S512, according to the cosine similarity matrix CS and the sample cosine similarity mean vector in the step S511Calculating a similarity matrix S; />
S513, calculating to obtain an initial neighborhood parameter K according to the similarity matrix S obtained in the step S512.
More specifically, in step S511, the cosine similarity matrix CS is a symmetric matrix, and the calculation formula is as follows:
wherein z is a feature dimension sequence number, D is the dimension of the high-dimensional degradation feature, i and j are sample sequence numbers; CS (i, j) is a sample in the cosine similarity matrix CS; x is x i And x j Two samples in the high-dimensional degradation feature X,and->Respectively sample x i And sample x j Is included in the sample points.
More specifically, in step S512, the initial neighborhood parameter matrix S may be set to:
wherein i, j is the sample number, n is the total number of samples;the mean value of cosine similarity between the ith sample and all other samples; s (i, j) are samples in the initial neighborhood parameter matrix S.
If sample x i And sample x j The cosine similarity between the two samples is larger than the average value of the cosine similarity between all adjacent samplesI.e. consider sample x i And sample x j Greater similarity, sample x j At sample x i Neighbor location, otherwise sample x j Not sample x i Is a neighbor of (a).
More specifically, in step S513, sample x i The neighborhood parameter matrix K of (1) can be expressed as:
wherein i, j is the sample number, n is the total number of samples; s (i, j) is a sample in the initial neighborhood parameter matrix S; k (k) i For samples in the neighborhood parameter matrix K, the sample x is represented i Is used to determine the initial neighborhood parameters of (1).
By the method, the initial neighborhood parameter of all samples is calculated as K 1 ，k 2 ，...，k n ]。
S520, neighborhood parameter adjustment:
the density distribution of the high-dimensional degradation features has a larger influence on the setting of the neighborhood parameters, if the density of the local range of the sample is larger, the Euclidean distance between the sample point and other points in the neighborhood range is smaller, the similarity is higher, and otherwise, the similarity between the sample point and other points in the neighborhood range is lower. At the same time, the magnitude of the intrinsic popular curvature of the high-dimensional degradation feature has a larger influence on the neighborhood parameter, for example, at the larger popular curvature, the neighborhood parameter should be set with a smaller parameter value in order to ensure that the correct low-dimensional feature is extracted. The estimation of the local popular curvature is carried out by adopting the ratio of the Euclidean distance between two samples and the geodesic distance, and the larger the ratio is, the smaller the popular curvature is, and the larger the popular curvature is on the contrary. Thus, it is necessary to comprehensively consider the density distribution of the high-dimensional degradation feature and the magnitude of the intrinsic popular curvature of the high-dimensional degradation feature for the adjustment of the initial neighborhood parameters. The invention estimates the distribution density of the high-dimensional degradation characteristic sample and the inherent popular curvature of the high-dimensional degradation characteristic, and adjusts the neighborhood parameters more accurately.
Based on the related parameters obtained in the step S510, the neighborhood parameter adjustment method provided by the invention comprises the following steps:
s521, solving the Euclidean distance D between all samples in the high-dimensional degradation feature set X e Distance D from geodesic g ；
S522, obtaining local density P and global density average value among all samples in the high-dimensional degradation feature set X
S523, calculating the ratio of the sum of the local Euclidean distances and the sum of the local geodesic distances to obtain the local popular curvature Q and the global popular curvatureMeanwhile, calculating the ratio of the Euclidean distance to the geodesic distance to obtain D'; />
S524, the neighborhood parameter K is adjusted by the neighborhood parameter adjustment method to obtain K c At the same time, the Euclidean distance matrix D is updated by the ratio D' of Euclidean distance to geodesic distance e Obtaining updated Euclidean distance matrix D e 。
In the step S521 of the method,
geodesic distance D between all samples of the bearing high-dimensional degradation feature g The iterative expression of (2) is:
d g (i，j)＝min{d g (i，j)，d g (i，m)+d g (m，j)}
wherein min is the minimum value; i, j, m are sample sequence numbers, and i, j, m are mutually unequal; n is the total number of samples; dg (i, j) is the geodesic distance D g Is a sample of (b).
More specifically, in step S522, sample x i The formula for estimating the local density vector p based on gaussian kernel probability density is as follows:
Wherein n is the total number of samples, i, j is the sample number; p is p i Is a sample in the local density vector P; k (k) i Is a sample in the preliminary neighborhood parameter K; d, d c (i, j) is sample x i And sample x j Euclidean distance between them.
Global density of high-dimensional degradation features, mean value estimated by Gaussian kernel probability densityThe expression is as follows:
wherein n is the total number of samples, i is the sample point number; p is p i Is a sample in the local density vector P.
Local density vector P and global density obtained based on Gaussian kernel probability density estimation methodThe preliminary adaptive neighborhood parameters may be adjusted.
More specifically, in step S523, the local popular curvature matrix Q of the high-dimensional degradation feature is a ratio of the sum of euclidean distances in the sample neighborhood range to the sum of geodesic distances, and the expression is as follows:
wherein i, j is a sample number; d, d e (i, j) and d g (i, j) are respectively Euclidean distance matrices D e Distance matrix D with geodesic g Is a sample of (a); k (k) i Is a sample in the preliminary neighborhood parameter K; q i Is a sample in the local popular curvature matrix Q.
Global popular curvature is as follows:
wherein i is a sample number, and n is the total number of samples; q (i) is a sample in the local popular curvature matrix Q.
The local popular curvature can be estimated by calculating the ratio of the sum of the local Euclidean distances to the sum of the local geodesic distances, and the ratio of the local popular curvature and the global popular curvature obtained by the method can be used for adjusting the preliminary self-adaptive neighborhood parameters.
Meanwhile, D 'is obtained by calculating the ratio of the corresponding Euclidean distance to the geodesic distance between each sample, and the formula is D' (i, j) =d e (i，j)/d g (i, j) and if the ratio d' (i, j) is greater than a constant y, consider sample x i And sample x j Sample d in its Euclidean distance matrix on the same local popular surface with smaller folded surface e The (i, j) value is unchanged; sample x is considered if the ratio d' (i, j) is less than a constant y i And sample x j On two popular sides with larger folding sides, let d e (i, j) = infinity, update euclidean distance matrix D e Can exclude sample x i And sample x j The possibility of becoming a neighborhood can be effectively ensured by the method that the extracted low-dimensional features do not overlap with samples at the places with larger popular curvatures, and the value of gamma is set empiricallyAnd (5) placing.
More specifically, in step S524, based on the above-mentioned parameters, the neighborhood parameter adjustment method is as follows:
Wherein i is a sample number; k (k) i Is a sample in the preliminary neighborhood parameter K;is the adjusted neighborhood parameter K c Is a sample of (a); floor (·) represents rounding down; q i Is a sample in the local popular curvature matrix Q; p is p i Is a sample in the local density vector P.
S530, fusion degradation index extraction:
based on the adjusted neighborhood parameter K obtained in section S520 c And updating the Euclidean distance D e ' extracting the low-dimensional characteristics obtained after the high-dimensional degradation characteristic mapping by a local linear embedding method. Under the following loss function and weight constraint conditions, calculating a Lagrangian multiplier method, and solving to obtain an objective function: MY (MY) T The equation is solved for =λ ' YT to obtain the eigenvector Y and eigenvalue λ ' of the eigenvector M '
Constraint of->
M＝(I i -w i )(I i -w i ) T
Wherein W represents a set of all sample weight vectors; x is x i And x j Two samples in the high-dimensional degradation feature X, respectively; w (w) i For sample x i Is a weight vector of (2); w (w) ij Is w i Sample points in (a); i i A D-dimensional vector representing all 1 s; m is a feature matrix.
Y is a matrix formed by feature vectors of the feature matrix M, in order to reduce high-dimensional data to d dimensions, feature vectors corresponding to the minimum d non-zero feature values of M are taken, and feature value decomposition is carried out on the matrix M to obtain low-dimensional representation Y of the high-dimensional degradation feature X: X-Y= [ Y ] 2 ，y 3 ，...，y d+1 ]Wherein y is i I=2, 3, …, d+1, representing the eigenvector corresponding to the first d nonzero eigenvalues of matrix M, will be eigenvalue y 2 The fusion degradation index extracted by the method is used for monitoring the fault state of the nuclear power main pump bearing.
Step S600, exponentially fitting and fusing degradation index curves, and determining different degradation state threshold intervals according to experience:
f(t)＝λ+βexp(ηt)
wherein: λ, β, η are parameters of an exponential function, which are determined by fitting to a fused degradation indicator.
Empirically determining four degradation phases of a bearing, labeled S 1 、S 2 、S 3 、S 4 Wherein S is 1 The normal wear period is the longest running time for the normal running period; s is S 2 、S 3 The method is a weak fault stage and a fault deep development stage, and maintenance personnel need to provide proper maintenance at the moment so as to prevent the bearing from being damaged more seriously; s is S 4 In the failure stage, the bearing cannot work normally due to the too serious degree of failure, and the dividing threshold is marked as χ 1 、χ 2 、χ 3 . The correspondence is shown in the following table:
the corresponding interval threshold is shown in the following table.
Degradation phase and threshold interval
To facilitate an understanding of embodiments of the present disclosure, embodiment 1
In step S100, vibration acceleration data are collected, and the time duration of the vibration acceleration data of a section of acceleration life test collected on a nuclear power main pump bearing is about 150 hours. As shown in fig. 2, there are a plurality of blank data and a large number of outliers in the segment of data, which need to be cleaned.
In step S200, missing, blank dot, and outlier cleaning:
step S201, missing data and blank invalid data cleaning:
1, for single-point missing data, polynomial fitting is carried out on the basis of a plurality of samples before the missing data by a polynomial fitting method, and the missing data is filled. In this embodiment, the first 30 samples of the missing data are fitted by a polynomial, i.e., the polynomial shown below is solved and calculated from p (t) i ) Substitution of deletion points:
p(t i )＝a 0 +a 1 t i +a 2 t i 2 +…+a n t i n and make (1)
Wherein t is i To fit the data sample sequence number, p (t i ) At t i Fitting data corresponding to the points.
2 for continuous blank invalid points, by judging whether there are 2 or more continuous points with zero amplitude, i.e. x i ＝x i+1 = … =0, and if present, successive invalid sample points of zero amplitude are removed from the vibration acceleration data.
Step S202, outlier cleaning:
1 estimation of the Normal data amplitude Range
Vibration acceleration data x for presence of outliers i Determining a range of normal data based on the radar criterion:
normal data x i The range is as follows:
wherein the average valueStandard deviation->n is vibration acceleration data x i Is a number of samples of (a).
The judgment mode of the outlier is as follows:
judging whether each sample is in a section or not Within, i.e. judge beta i Is less than or equal to 3, or beta i > 3, if sample x i Corresponding beta i Sample x is less than or equal to 3 i Is the normal point, otherwise if sample x i Corresponding beta i > 3, sample x i For outliers, a cleaning process is required, and in order to ensure the reliability of the replacement data as much as possible, the present disclosure takes the average value of 6 samples adjacent to the value to replace the value, and the expression is:
as shown in fig. 3, the vibration acceleration data after cleaning is effectively cleaned by outliers and blank invalid points contained in the data, and the data quality is improved.
Example 2
In step S300, for the vibration acceleration data amplitude modulation phenomenon caused by the variable rotation speed working condition, the present disclosure uses fault simulation data under different rotation speed working conditions extracted by a certain nuclear power main pump bearing fault simulation experiment as analysis data, the rotation speed working conditions are 600RPM, 800RPM, 1000RPM, 1200RPM and 1400RPM, the sample length of each rotation speed data is 100, and the specific analysis is as follows:
for vibration data x containing variable rotation speed working condition i Taking his Hilbert envelope e i As shown in the following formula：
e i ＝|x i +jH(x i )|
Wherein H (x) i ) Representing de-vibration data x i Is a hilbert transform of (c); x is x i +jH(x i ) Plural, |x i +jH(x i ) The I represents the complex number x i +jH(x i ) Is a mold of (a).
For envelope data e i The magnitudes of the first m magnitudes are ordered in descending order, and the average μ of the first m magnitudes is calculated:
wherein max (e i ) m Representing envelope data e i The sum of the first m maximum amplitude values is averaged over a plurality of amplitude values to reduce disturbances caused by non-stationary factors in the vibration data.
Quantized vibration acceleration data u after eliminating influence of rotation speed factor i The method is characterized by comprising the following steps:
as shown in fig. 4, the amplitude of the square root extracted from the processed vibration acceleration data is greatly reduced due to amplitude fluctuation caused by the influence of the rotation speed under different rotation speeds, as shown by the amplitude graph of the square root extracted from the 5 rotation speed working conditions of the fault simulation experimental data of the nuclear power main pump bearing.
In step S400, high-dimensional degradation features are extracted. And a section of bearing life data collected in a nuclear power main pump bearing accelerated life experiment has a data sample length of n=5750 and a data duration of about 150 hours. Vibration acceleration data in three directions can be acquired through a triaxial vibration acceleration sensor, and 15-dimensional high-dimensional degradation characteristics X [ X ] shown below can be extracted after cleaning in the steps S200-S300 1 ，X 2 ，X 3 ，…，X 15 ]：
| Reference numerals | Features (e.g. a character) | Reference numerals | Features (e.g. a character) |
| X 1 | x-axis kurtosis characteristics | X 9 | y-axis waveform index |
| X 2 | x-axis waveform index | X 10 | Y-axis kurtosis index |
| X 3 | x-axis kurtosis index | X 11 | y-axis frequency domain features 3 |
| X 4 | x-axis margin index | X 12 | y-axis frequency domain features 4 |
| X 5 | x-axis frequency domain features 3 | X 13 | z-axis waveform index |
| X 6 | x-axis frequency domain features 4 | X 14 | z-axis frequency domain features 3 |
| X 7 | x-axis frequency domain features 12 | X 15 | z-axis frequency domain features 4 |
| X 8 | y-axis kurtosis characteristics |
The features in the above table are shown below
Kurtosis characteristics:wherein->
Waveform index:/>
kurtosis index:wherein->
Margin index:
frequency domain feature 3:wherein->
Frequency domain feature 4:wherein->
Frequency domain features 12Wherein->
In the above calculation expression, i is the vibration acceleration data sequence number, n is the total number of data samples, and n=5750; u (u) i Representing a vibration acceleration data sample; y is j Representing vibration acceleration data u i Frequency domain data obtained after FFT conversion; j is vibration acceleration data u i Frequency domain data y of (2) j N represents vibration acceleration data u i N=5750; f (f) j Representing frequency domain data y j Frequency components in the frequency spectrum.
In particular, when the vibration acceleration sensor is uniaxial, the degradation characteristic X can be extracted 1 ～X 7 A total of 7 degradation features as high-dimensional degradation features; when the vibration acceleration sensor is in two axial directions, the degradation characteristic X can be extracted 1 ～X 12 A total of 7 degradation features were used as high-dimensional degradation features.
In step S500, a fusion degradation index is extracted from the high-dimensional degradation feature, and the fusion degradation index is used for monitoring the bearing fault state of the wind driven generator transmission system, and the specific process is as follows:
s510, initial neighborhood parameter setting, the general steps of which are as follows:
s511, solving the cosine similarity matrix CS and the sample cosine similarity mean vector among all samples in the high-dimensional degradation feature set X
S512, according to the cosine similarity matrix CS and the sample cosine similarity mean vector in the step S511Calculating a similarity matrix S;
s513, calculating to obtain an initial neighborhood parameter K according to the similarity matrix S obtained in the step S512.
More specifically, in step S511, the cosine similarity matrix CS is a symmetric matrix, and the calculation formula is as follows:
wherein z is a feature dimension sequence number, D is the dimension of the high-dimensional degradation feature, i and j are sample sequence numbers; CS (i, j) is a sample in the cosine similarity matrix CS; x is x i And x j Two samples in the high-dimensional degradation feature X,and->Respectively sample x i And sample x j Is included in the sample points. />
More specifically, in step S512, the initial neighborhood parameter matrix S may be set to:
wherein i, j is the sample number, n is the total number of samples; The mean value of cosine similarity between the ith sample and all other samples; s (i, j) are samples in the initial neighborhood parameter matrix S.
If sample x i And sample x j The cosine similarity between the two samples is larger than the average value of the cosine similarity between all adjacent samplesI.e. consider sample x i And sample x j Greater similarity, sample x j At sample x i Neighbor location, otherwise sample x j Not sample x i Is a neighbor of (a).
More specifically, in step S513, sample x i The neighborhood parameter matrix K of (1) can be expressed as:
wherein i, j is the sample number, n is the total number of samples; s (i, j) is a sample in the initial neighborhood parameter matrix S; k (k) i For samples in the neighborhood parameter matrix K, the sample x is represented i Is used to determine the initial neighborhood parameters of (1).
By the method, the initial neighborhood parameter of all samples is calculated as K 1 ，k 2 ，...，k n ]。
S520, neighborhood parameter adjustment:
based on the related parameters obtained in the step S510, the neighborhood parameter adjustment method provided by the invention comprises the following steps:
s521, solving the Euclidean distance D between all samples in the high-dimensional degradation feature set X e Distance D from geodesic g ；
S522, obtaining local density P and global density average value among all samples in the high-dimensional degradation feature set X
S523, calculating the ratio of the sum of the local Euclidean distances and the sum of the local geodesic distances to obtain the local popular curvature Q and the global popular curvature Meanwhile, calculating the ratio of the Euclidean distance to the geodesic distance to obtain D';
s524, the neighborhood parameter K is adjusted by the neighborhood parameter adjustment method to obtain K c At the same time, the Euclidean distance matrix D is updated by the ratio D' of Euclidean distance to geodesic distance e Obtaining updated Euclidean distance matrix D e ′。
More specifically, in step S521, the geodesic distance matrix D g The method comprises the steps of obtaining by a geodesic distance calculation method in an equidistant feature mapping (ISOMAP) algorithm;
geodesic distance D between all samples of the bearing high-dimensional degradation feature g The iterative expression of (2) is:
d g (i，j)＝min{d g (i，j)，d g (i，m)+d g (m，j)}
wherein min is the minimum value; i, j, m are sample sequence numbers, and i, j, m are mutually unequal; n is the total number of samples; dg (i, j) is the geodesic distance D g Is a sample of (b).
Euclidean distance matrix D e The method is characterized by comprising the following steps:
wherein D is a degradation characteristic dimension, z is a characteristic dimension sequence number, i, j is a sample sequence number;and->Respectively sample x i And x j Sample points in (a); d, d e (i, j) is Euclidean distance matrix D e Is a sample of (b).
More specifically, in step S522, sample x i The formula for estimating the local density vector p based on gaussian kernel probability density is as follows:
/>
wherein n is the total number of samples, i, j is the sample number; p is p i Is a sample in the local density vector P; k (k) i Is a sample in the preliminary neighborhood parameter K; d, d e (i, j) is sample x i And sample x j Euclidean distance between them.
Global density of high-dimensional degradation features, mean value estimated by Gaussian kernel probability densityThe expression is as follows:
wherein n is the total number of samples, i is the sample point number; p is p i Is a sample in the local density vector P.
Local density vector P and global density obtained based on Gaussian kernel probability density estimation methodThe preliminary adaptive neighborhood parameters may be adjusted.
More specifically, in step S523, the local popular curvature matrix Q of the high-dimensional degradation feature is a ratio of the sum of euclidean distances in the sample neighborhood range to the sum of geodesic distances, and the expression is as follows:
wherein i, j is a sample number; d, d e (i, j) and d g (i, j) are respectively Euclidean distance matrices D e Distance matrix D with geodesic g Is a sample of (a); k (k) i Is a sample in the preliminary neighborhood parameter K; q i Is a sample in the local popular curvature matrix Q.
Global popular curvature is as follows:
wherein i is a sample number, and n is the total number of samples; q (i) is a sample in the local popular curvature matrix Q.
The local popular curvature can be estimated by calculating the ratio of the sum of the local Euclidean distances to the sum of the local geodesic distances, and the ratio of the local popular curvature and the global popular curvature obtained by the method can be used for adjusting the preliminary self-adaptive neighborhood parameters.
Meanwhile, D 'is obtained by calculating the ratio of the corresponding Euclidean distance to the geodesic distance between each sample, and the formula is D' (i, j) =d e (i，j)/d g (i, j) and if the ratio d' (i, j) is greater than a constant y, consider sample x i And sample x j Sample d in its Euclidean distance matrix on the same local popular surface with smaller folded surface e The (i, j) value is unchanged; sample x is considered if the ratio d' (i, j) is less than a constant y i And sample x j On two popular sides with larger folding sides, let d e (i, j) = infinity, update euclidean distance matrix D e Can exclude sample x i And sample x j The possibility of being a neighborhood is effectively ensured by the method, and the extracted low-dimensional features are not overlapped at the place with larger popular curvature, wherein the value of gamma is set empirically.
More specifically, in step S524, based on the above-mentioned parameters, the neighborhood parameter adjustment method is as follows:
wherein i is a sample number; k (k) i Is a sample in the preliminary neighborhood parameter K;is the adjusted neighborhood parameter K c Is a sample of (a); floor (·) represents rounding down; q i Is a sample in the local popular curvature matrix Q; p is p i Is a sample in the local density vector P.
S530, fusion degradation index extraction:
based on the adjusted neighborhood parameter K obtained in section S520 c And updating the Euclidean distance D e ' extracting the low-dimensional characteristics obtained after the high-dimensional degradation characteristic mapping by a local linear embedding method. Under the following loss function and weight constraint conditions, calculating a Lagrangian multiplier method, and solving to obtain an objective function: MY (MY) T ＝λ′Y T Solving the formula to obtain a feature vector Y and a feature value lambda 'of the feature matrix M'
Constraint of->
M＝(I i -w i )(I i -w i ) T
Wherein W represents a set of all sample weight vectors; x is x i And x j Two samples in the high-dimensional degradation feature X, respectively; w (w) i For sample x i Is a weight vector of (2); w (w) ij Is w i Sample points in (a); i i A D-dimensional vector representing all 1 s; m is a feature matrix.
Y is the eigenvector of the eigenvector matrix MIn order to reduce the high-dimensional data to d-dimension, the feature vector corresponding to the minimum d non-zero feature values of M is taken, and the low-dimensional representation Y of the high-dimensional degradation feature X can be obtained by carrying out feature value decomposition on the matrix M: X-Y= [ Y ] 2 ，y 3 ，...，y d+1 ]Wherein y is i I=2, 3, …, d+1, representing the eigenvector corresponding to the first d nonzero eigenvalues of matrix M, will be eigenvalue y 2 The fusion degradation index extracted by the method is used for monitoring the fault state of the nuclear power main pump bearing.
In step S600, a degradation trend curve is exponentially fitted, and different degradation state threshold intervals are determined:
f(t)＝λ+βexp(ηt)
wherein: λ, β, η are parameters of an exponential function, which are determined by fitting to a fused degradation indicator.
The four degradation stages of the bearing in different degradation states are marked as S according to experience determination, and the four degradation stages are marked as S according to the combination of degradation indexes and an exponential fit curve as shown in figure 5 1 、S 2 、S 3 、S 4 Wherein S is 1 The normal wear period is the longest running time for the normal running period; s is S 2 、S 3 The method is a weak fault stage and a fault deep development stage, and maintenance personnel need to provide proper maintenance at the moment so as to prevent the bearing from being damaged more seriously; s is S 4 In order to fail, the bearing cannot work normally at this time because the failure degree is too serious, and the corresponding interval threshold is shown in the following table.
TABLE 1 degradation phase and threshold interval
Although embodiments of the present invention have been described above with reference to the accompanying drawings, the present invention is not limited to the specific embodiments and fields of application described above, which are illustrative, instructive, and not restrictive. Those skilled in the art, having the benefit of this disclosure, may effect numerous forms of the invention without departing from the scope of the invention as claimed.
Claims (6)
1. A nuclear power main pump bearing fault early warning method based on fusion degradation indexes comprises the following steps:
in the first step (S100), a vibration acceleration sensor is arranged on a nuclear power main pump bearing bushing shell to acquire vibration acceleration data;
in the second step (S200), the missing points, blank points and outliers in the vibration acceleration data are cleaned, and when the outliers in the vibration acceleration data are cleaned, the normal data x is obtained i The range is as follows:mean->The estimation method comprises the following steps: />The standard deviation sigma estimation method comprises the following steps: />The outlier judging mode is as follows: />
Sample x is determined sample by sample i Whether or not it is in an intervalIn addition, beta is determined i Is less than or equal to 3, or beta i > 3, if sample x i Corresponding beta i Sample x is less than or equal to 3 i Is the normal point, otherwise if sample x i Corresponding beta i > 3, sample x i For outliers, a cleaning process is required, and blank data appears in operationIn the process, the main pump starts and stops the stage, and is collected by a vibration acceleration sensor;
in the third step (S300), the data of the amplitude modulation phenomenon caused by the variable rotation speed is removed from the vibration acceleration data to obtain quantized vibration acceleration data;
in a fourth step (S400), high-dimensional degradation features are extracted from the quantized vibration acceleration data, wherein the high-dimensional degradation features include:
Kurtosis characteristics:wherein->
Waveform index:
kurtosis index:wherein->
Margin index:
frequency domain feature 3:wherein->
Frequency domain feature 4:wherein->
Frequency domain features 12:wherein->
Wherein i is a vibration acceleration data sequence number, and n is the total number of data samples; u (u) i Representing a quantized vibration acceleration data sample; y is j Representing quantized vibration acceleration data u i Frequency domain data obtained after FFT conversion; j is quantized vibration acceleration data u i Frequency domain data y of (2) j N represents quantized vibration acceleration data u i Is a sample total number of frequency domain data; f (f) j Representing frequency domain data y j Frequency components in the frequency spectrum;
the cosine similarity matrix CS is a symmetric matrix, and the calculation formula is as follows:
wherein z is a feature dimension sequence number, D is the dimension of the high-dimensional degradation feature, i and j are sample sequence numbers; CS (i, j) is a sample in the cosine similarity matrix CS; x is x i And x j Two samples in the high-dimensional degradation feature X,and->Respectively sample x i And sample x j In (c) the sample points in (c),
high-dimensional bearing backingGeodesic distance D between all samples of chemical characterization g The iterative expression of (2) is:
d g (i，j)＝min{d g (i，j)，d g (i，m)+d g (m，j)}
wherein min is the minimum value; i, j, m are sample sequence numbers, and i, j, m are mutually unequal; n is the total number of samples; d, d g (i, j) is the geodesic distance D g Is a sample of (a);
in a fifth step (S500), extracting a fusion degradation index for monitoring the failure state of the nuclear power main pump bearing from the high-dimensional degradation characteristic, wherein the extraction of the fusion degradation index for monitoring the failure state of the nuclear power main pump bearing from the high-dimensional degradation characteristic comprises,
solving cosine similarity matrix CS and sample cosine similarity mean vector among all samples in high-dimensional degradation feature set X
According to the cosine similarity matrix CS and the sample cosine similarity mean vectorCalculating a similarity matrix S;
calculating an initial neighborhood parameter K according to the similarity matrix S, wherein the initial neighborhood parameter K is expressed as:
wherein i, j is the sample number, n is the total number of samples; s (i, j) are samples in the similarity matrix s; k (k) i For samples in the neighborhood parameter matrix K, the sample X is represented i Is used for generating initial neighborhood parameters;
obtaining the Euclidean distance D between all samples in the high-dimensional degradation feature set X e Distance D from geodesic g ；
Averaging the local density P and global density across all samples in the high dimensional degradation feature set X
Calculating the ratio of the sum of the local Euclidean distances and the sum of the local geodesic distances to obtain the local popular curvature Q and the global popular curvatureMeanwhile, calculating the ratio of the Euclidean distance to the geodesic distance to obtain D';
Adjusting the initial neighborhood parameter K based on the neighborhood parameter adjusting method to obtain an adjusted neighborhood parameter K c At the same time, the Euclidean distance matrix D is updated by the ratio D' of Euclidean distance to geodesic distance e Obtaining updated Euclidean distance matrix D e ′；
Based on the adjusted neighborhood parameter K c And updated Euclidean distance matrix D e ' extracting low-dimensional features obtained after high-dimensional degradation feature mapping by a local linear embedding method, and solving to obtain an objective function by calculating a Lagrangian multiplier method: MY (MY) T ＝λ′Y T Solving the formula to obtain a feature vector Y and a feature value lambda 'of the feature matrix M'
Constraint of->
M＝(I i -w i )(I i -W i ) T ，
Wherein W represents a set of all sample weight vectors; x is X i And x j Two samples in the high-dimensional degradation feature X, respectively; w (W) i For sample x i Is a weight vector of (2); w (w) ij Is W i Sample points in (a); i i A D-dimensional vector representing all 1 s; m is a characteristic matrix, wherein M is a characteristic matrix,
y is a matrix formed by feature vectors of the feature matrix M, and the matrix M is subjected to feature value decomposition to obtain a low-dimensional table of the high-dimensional degradation feature XShow Y: X-Y= [ Y ] 2 ，y 3 ，...，y d+1 ]Wherein y is i I=2, 3,..d+1, representing the eigenvector corresponding to the first d non-zero eigenvalues of matrix M, will be eigenvalue y 2 As a fusion degradation index for a fault state monitoring index of a nuclear power main pump bearing;
In a sixth step (S600), thresholds for different degradation phases are determined based on the fused degradation indicators.
2. The method according to claim 1, wherein in the second step (S200), missing points in the vibration acceleration data are cleaned by polynomial fitting, and samples x with serial numbers 1,2, …, n preceding the missing points are cleaned up 1 ，x 2 ，…，x n Fitting on the basis of the data, and solving p (t) i )＝a 0 +a 1 t i +a 2 t i 2 +…+a n t i n And make (1)And is composed of p (t) i ) Substitution of deletion points, wherein t i To fit the data samples, p (t i ) At t i Fitting data corresponding to points, a 0 ，a 1 ，…a n Is the fitting coefficient.
3. The method according to claim 1, wherein in the second step (S200), when the blank point in the vibration acceleration data is cleaned, it is determined whether there are 2 or more continuous points with zero amplitude, x i ＝x i+1 = … =0, and if present, successive sample points of zero amplitude are removed from the vibration acceleration data.
4. The method of claim 1, wherein the value is replaced by a mean of a number of samples adjacent to the outlier expressed as:
wherein N is an outlier sample x i The total number of adjacent sample points is taken, i is the sample number.
5. The method of claim 1, wherein removing data of amplitude modulation phenomena caused by varying rotational speed from the vibration acceleration data comprises,
For vibration data x containing variable rotation speed working condition i Taking his Hilbert envelope e i ，e i ＝|x i +jH(x i ) I, wherein H (x i ) Representing de-vibration data x i Is a hilbert transform of (c); x is x i +jH(x i ) Plural, |x i +jH(x i ) The I represents the complex number x i +jH(x i ) Is the modulus, j is the imaginary unit,
for envelope data e i The magnitudes of the first m magnitudes are ordered in descending order, and the average μ of the first m magnitudes is calculated:wherein max (e i ) m Representing envelope data e i The sum of the first m maximum amplitude values, the quantized vibration acceleration data x 'after eliminating the influence of the rotation speed factor' i The method is characterized by comprising the following steps: />
6. The method according to claim 1, wherein in the sixth step (S600) determining thresholds for different degradation phases based on the fused degradation indicators comprises exponentially fitting fused degradation indicator curves to determine different degradation state threshold intervals: f (t) =λ+βexp (ηt), wherein: λ, β, η are parameters of an exponential function, which are determined by fitting to a fused degradation indicator.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010958398.5A CN112326241B (en) | 2020-09-11 | 2020-09-11 | Nuclear power main pump bearing fault early warning method based on fusion degradation index |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010958398.5A CN112326241B (en) | 2020-09-11 | 2020-09-11 | Nuclear power main pump bearing fault early warning method based on fusion degradation index |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112326241A CN112326241A (en) | 2021-02-05 |
| CN112326241B true CN112326241B (en) | 2023-10-13 |
Family
ID=74303157
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010958398.5A Active CN112326241B (en) | 2020-09-11 | 2020-09-11 | Nuclear power main pump bearing fault early warning method based on fusion degradation index |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112326241B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113691398B (en) * | 2021-08-13 | 2023-04-07 | 北京金山云网络技术有限公司 | Node bandwidth prediction method and device, electronic equipment and storage medium |
| CN113537156B (en) * | 2021-09-06 | 2021-12-14 | 航天智控(北京)监测技术有限公司 | Vibration data cleaning method based on interval standard deviation and spectrum analysis |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102944435A (en) * | 2012-10-25 | 2013-02-27 | 北京航空航天大学 | Health assessment and fault diagnosis method for rotating machinery based on fisher discriminant analysis and mahalanobis distance |
| CN109556863A (en) * | 2018-06-13 | 2019-04-02 | 南京工业大学 | A kind of acquisition of large-scale turntable bearing Vibration Signal in Frequency Domain and processing method based on MSPAO-VMD |
| CN109916627A (en) * | 2019-03-27 | 2019-06-21 | 西南石油大学 | Bearing fault detection and diagnosis based on Active Learning |
| CN111175045A (en) * | 2020-01-08 | 2020-05-19 | 西安交通大学 | Method for cleaning vibration acceleration data of locomotive traction motor bearing |
| CN111307453A (en) * | 2020-03-20 | 2020-06-19 | 朗斯顿科技(北京)有限公司 | Transmission system fault diagnosis method based on multi-information fusion |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160245686A1 (en) * | 2015-02-23 | 2016-08-25 | Biplab Pal | Fault detection in rotor driven equipment using rotational invariant transform of sub-sampled 3-axis vibrational data |
| US10495693B2 (en) * | 2017-06-01 | 2019-12-03 | General Electric Company | Wind turbine fault detection using acoustic, vibration, and electrical signals |
-
2020
- 2020-09-11 CN CN202010958398.5A patent/CN112326241B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102944435A (en) * | 2012-10-25 | 2013-02-27 | 北京航空航天大学 | Health assessment and fault diagnosis method for rotating machinery based on fisher discriminant analysis and mahalanobis distance |
| CN109556863A (en) * | 2018-06-13 | 2019-04-02 | 南京工业大学 | A kind of acquisition of large-scale turntable bearing Vibration Signal in Frequency Domain and processing method based on MSPAO-VMD |
| CN109916627A (en) * | 2019-03-27 | 2019-06-21 | 西南石油大学 | Bearing fault detection and diagnosis based on Active Learning |
| CN111175045A (en) * | 2020-01-08 | 2020-05-19 | 西安交通大学 | Method for cleaning vibration acceleration data of locomotive traction motor bearing |
| CN111307453A (en) * | 2020-03-20 | 2020-06-19 | 朗斯顿科技(北京)有限公司 | Transmission system fault diagnosis method based on multi-information fusion |
Non-Patent Citations (1)
| Title |
|---|
| "基于PCA－KNN聚类的通用在线故障诊断算法设计";赵晓君等;《计算机测量与控制》;20150831;第38卷(第8期);第2763页第1栏 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112326241A (en) | 2021-02-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111275288B (en) | XGBoost-based multidimensional data anomaly detection method and device | |
| Huang et al. | Fault diagnosis of bearing in wind turbine gearbox under actual operating conditions driven by limited data with noise labels | |
| CN112326241B (en) | Nuclear power main pump bearing fault early warning method based on fusion degradation index | |
| CN108376298B (en) | Early warning and diagnosing method for temperature faults of engine of wind turbine generator | |
| CN109829402B (en) | GS-SVM-based bearing damage degree diagnosis method under different working conditions | |
| CN105631596B (en) | Equipment fault diagnosis method based on multi-dimensional piecewise fitting | |
| Jin et al. | Anomaly detection of cooling fan and fault classification of induction motor using Mahalanobis–Taguchi system | |
| CN110685868A (en) | Wind turbine generator fault detection method and device based on improved gradient elevator | |
| CN109597315B (en) | Method, equipment and system for identifying health degradation state of mechanical equipment | |
| CN116226646B (en) | Method, system, equipment and medium for predicting health state and residual life of bearing | |
| CN111523081A (en) | Aircraft engine fault diagnosis method based on enhanced gated cyclic neural network | |
| CN111415070A (en) | Wind turbine generator gearbox oil temperature over-temperature fault early warning method based on SCADA data | |
| CN112633098A (en) | Fault diagnosis method and system for rotary machine and storage medium | |
| CN113239534B (en) | Fault and service life prediction method and device of wind generating set | |
| CN111209934A (en) | Fan fault prediction and alarm method and system | |
| CN114841580A (en) | Generator fault detection method based on hybrid attention mechanism | |
| Shi et al. | Study of wind turbine fault diagnosis and early warning based on SCADA data | |
| CN116467592A (en) | Production equipment fault intelligent monitoring method and system based on deep learning | |
| Yuan et al. | Fault detection of rolling bearing based on principal component analysis and empirical mode decomposition | |
| CN113654651B (en) | Method for extracting early degradation features of strong robust signal and monitoring running state of equipment | |
| CN112697268B (en) | Motor anomaly detection integrated algorithm based on T-SNE | |
| CN112699598A (en) | Intelligent diagnosis method and device for abnormal oil temperature of gear box | |
| CN115114770B (en) | Baseline self-adaptive auxiliary power device performance trend analysis method | |
| CN114877925B (en) | Comprehensive energy system sensor fault diagnosis method based on extreme learning machine | |
| CN117370826A (en) | Method for extracting health state characteristics in wind turbine generator operation data |
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 |