CN115758094A

CN115758094A - Method for extracting degradation index of speed reducer and predicting service life of speed reducer

Info

Publication number: CN115758094A
Application number: CN202310023178.7A
Authority: CN
Inventors: 韩旭; 王嘉; 陶友瑞; 张露予; 赵赢; 李本旺; 尹慧强
Original assignee: Hebei University of Technology
Current assignee: Hebei University of Technology
Priority date: 2023-01-09
Filing date: 2023-01-09
Publication date: 2023-03-07
Anticipated expiration: 2043-01-09
Also published as: CN115758094B

Abstract

The application provides a speed reducer degradation index extraction and service life prediction method, wherein the degradation index extraction method comprises the following steps: simultaneously extracting a current signal and a vibration signal in the harmonic reducer, respectively preprocessing the current signal and the vibration signal, and respectively performing modal decomposition and feature extraction; obtaining a plurality of current components and vibration components, and calculating the correlation among the components; then, correlation screening is carried out, and feature level data fusion dimensionality reduction is carried out to obtain a plurality of principal elements; and obtaining a degradation index curve by calculating the relationship between the pivot and the contribution rate of the pivot. The current and vibration information are combined and preprocessed, so that the influence of noise signals on fault characteristics can be reduced, and the fault missing detection rate and the fault false detection rate are effectively reduced; and modal decomposition feature extraction is carried out, and coefficient fusion dimensionality reduction is carried out by utilizing the correlation among all components, so that fault information contained in signals in the early stage of equipment operation can be effectively extracted, and the accuracy of an analysis result is improved.

Description

Method for extracting degradation index and predicting service life of speed reducer

Technical Field

The disclosure generally relates to the field of harmonic reducer degradation models, and particularly relates to a reducer degradation index extraction and life prediction method.

Background

According to the data it is shown that 40-70% of industrial robot failures are caused by the speed reducer. Because the working environment of the industrial robot is severe, the working mode relates to frequent starting and stopping, overload operation and the like, and the speed reducer, particularly the harmonic speed reducer with the tail end of the mechanical arm directly participating in operation, is easy to have the faults of flexible wheel abrasion, fatigue fracture, flexible bearing damage and the like, so that the transmission precision of the robot is reduced, the mechanical arm does not stably operate and the like, and even the whole machine fails.

Current signals acquired in an actual industrial environment are often subjected to strong noise interference to enable fault information to be submerged, observed current waveforms are mostly non-stationary signals, and a traditional signal processing method based on Fourier transform is not suitable any more. In order to better analyze and process such non-stationary signals, some scholars propose time-frequency analysis methods, which mainly include short-time fourier transform, wavelet transform, empirical mode decomposition, and the like.

In the prior art, a vibration signal or a current signal is independently adopted to carry out motor fault analysis or service life prediction, so that the fault undetected rate and the false detected rate are increased due to signal loss and signal interference. Meanwhile, the signals contain weak fault information in the early stage of the operation of the equipment, the fault characteristics contained in the signals are difficult to extract, and part of fault information can be separately expressed in one signal. It is difficult to preprocess against noise interference in the current signal, and therefore the accuracy of the subsequent analysis result is affected.

Disclosure of Invention

In view of the above-mentioned drawbacks or deficiencies in the prior art, it is desirable to provide a harmonic reducer degradation indicator extraction method based on current and vibration signals.

One aspect of the present application provides a method for extracting degradation indicators of a speed reducer, including:

acquiring an initial current signal and an initial vibration signal;

respectively carrying out signal preprocessing on the initial current signal and the initial vibration signal to obtain a first current signal and a first vibration signal;

performing modal decomposition and feature extraction on the first current signal and the first vibration signal respectively; obtaining a plurality of current signal components and a plurality of vibration signal components;

optionally selecting one of said current signal components as a reference current component; optionally selecting one of said vibration signal components as a reference vibration component;

respectively calculating the correlation of each current signal component which is not taken as the reference current component to obtain a plurality of current correlation coefficients;

respectively calculating the correlation of each vibration signal component which is not taken as the reference vibration component to obtain a plurality of vibration correlation coefficients;

eliminating current signal components with current correlation coefficients larger than a first set value, and taking the residual current signal components as current characteristic parameters; eliminating vibration signal components with vibration correlation coefficients larger than a second set value, and taking the residual vibration signal components as vibration characteristic parameters;

performing characteristic level data fusion on the current characteristic parameters and the vibration characteristic parameters to obtain a characteristic fusion parameter matrix;

reducing the dimension of the characteristic fusion parameter matrix to obtain a plurality of principal elements;

respectively calculating the contribution rates of the pivot elements to the fault information to obtain a plurality of contribution rates;

and multiplying the plurality of principal elements by the contribution rates of the principal elements, and summing to obtain a degradation index curve.

According to the technical scheme provided by the embodiment of the application, the current signal component with the highest kurtosis is selected as a reference current component; selecting the vibration signal component with the highest kurtosis as a reference vibration component.

According to the technical scheme provided by the embodiment of the application, the signal preprocessing comprises the following steps: filtering and denoising; respectively carrying out noise reduction processing on the initial current signal and the initial vibration signal by utilizing a correlation kurtosis deconvolution method;

substituting the initial current signal or the initial vibration signal as a first input signal into the correlation kurtosis deconvolution method, specifically comprising:

s3.1: acquiring a first input signal, a pulse period, a filter length and a time shift period;

s3.2: calculating to obtain a first filter function according to the first input signal, the pulse period, the filter length and the time shifting period;

s3.3: filtering the first input signal according to the first filter function to obtain a filtered signal;

s3.4: calculating the correlation kurtosis of the first input signal according to the first input signal; calculating according to the filtering signal to obtain the correlation kurtosis of the filtering signal;

s3.5: calculating the difference value of the correlation kurtosis of the first input signal and the correlation kurtosis of the filtering signal to obtain a correlation kurtosis difference value;

s3.6: and (3) judging: when the correlation kurtosis difference is smaller than a third set value, directly performing S3.9; otherwise, carrying out the next step;

s3.7: calculating according to the filtering signal, the pulse period, the filter length and the time shifting period to obtain a second filter function;

s3.8: taking the second filter function as the first filter function and returning to S3.3;

s3.9: taking the filtered signal as a first output signal;

if the initial current signal is used as a first input signal, the first output signal is used as the first current signal; and if the initial vibration signal is taken as a first input signal, taking the first output signal as the first vibration signal.

According to the technical scheme provided by the embodiment of the application, modal decomposition and feature extraction are respectively carried out on the first current signal and the first vibration signal by utilizing a complementary set empirical mode decomposition algorithm to obtain a plurality of current signal components and a plurality of vibration signal components;

taking the first current signal component or the first vibration signal component as a second input signal; the specific steps of performing modal decomposition and feature extraction by using the complementary set empirical mode decomposition algorithm comprise:

s4.1: acquiring a second input signal;

s4.2: adding multiple groups of white noise to the second input signal respectively to obtain multiple first set signals;

subtracting a plurality of groups of white noise from the second input signal respectively to obtain a plurality of second set signals;

s4.3: respectively calculating all maximum value points and all minimum value points of each first set signal to obtain a plurality of first maximum value points and a plurality of first minimum value points of each first set signal;

respectively calculating all maximum value points and all minimum value points of each second set signal to obtain a plurality of second maximum value points and a plurality of second minimum value points of each second set signal;

s4.4: respectively connecting the plurality of first maximum points and the plurality of first minimum points of each first set signal by using a spline curve to obtain a first maximum envelope curve and a first minimum envelope curve of each first set signal;

respectively connecting the second maximum points and the second minimum points of each second set signal by spline curves to obtain a second maximum envelope curve and a second minimum envelope curve of each second set signal;

s4.5: respectively calculating the average value of the first maximum envelope line and the first minimum envelope line of each first set signal to obtain a first average line of each first set signal;

respectively calculating the average value of the second maximum envelope line and the second minimum envelope line of each second set signal to obtain a second average line of each second set signal;

s4.6: calculating the difference between each first set signal and the corresponding first average line respectively to obtain a plurality of first difference functions;

respectively calculating the difference value of each second set signal and the corresponding second average line to obtain a plurality of second difference value functions;

s4.7: and (3) judging: when the plurality of first difference functions and the plurality of second difference functions all meet the condition of intrinsic mode functions, S4.9 is carried out; otherwise, carrying out the next step;

s4.8: taking each first difference function as a corresponding first set signal; taking each second difference function as a corresponding second set signal, and returning to S4.3;

s4.9: taking a plurality of the first difference functions as a plurality of first intrinsic mode functions and a plurality of the second difference functions as a plurality of second intrinsic mode functions;

s4.10: calculating the component of each first intrinsic mode function to obtain a plurality of first component parameters; calculating the component of each second eigenmode function to obtain a plurality of second component parameters; taking the first component parameter and the second component parameter as a set of component parameters;

s4.11: and (3) judging: when the number of the component parameters reaches a fourth set value, the next step is carried out; otherwise, returning to S4.3;

s4.12; respectively carrying out set average processing on each group of component parameters to obtain a plurality of signal component mean values;

if the first current signal component is taken as the second input signal, taking a plurality of signal component mean values as a plurality of current signal components; and if the first vibration signal component is taken as the second input signal, taking the average value of the plurality of signal components as a plurality of vibration signal components.

According to the technical scheme provided by the embodiment of the application, the Copula function is utilized to calculate the correlation of other current signal components to the reference current component and the correlation of other vibration signal components to the reference vibration component;

a plurality of the current signal components not serving as the reference current component or a plurality of the vibration signal components not serving as the reference vibration component are taken as a plurality of third input values; taking the reference current component or the reference vibration component as a fourth input value;

the specific steps of calculating the correlation by using the Copula function comprise:

s5.1: acquiring a plurality of third input values and fourth input values;

s5.2: selecting one of the third input values, which is not subjected to the calculation of the kernel distribution estimation value, as a first sample; taking the fourth input value as a second sample;

s5.3: calculating a kernel distribution estimated value of the first sample and the second sample; obtaining a first variable and a second variable;

s5.4: calculating linear correlation parameters of the first variable and the second variable; obtaining a correlation coefficient between the first sample and the second sample;

s5.5: and (3) judging: when the third input value which is not subjected to the calculation of the kernel distribution estimated value does not exist, the next step is carried out; otherwise, returning to S5.2;

s5.6: using a plurality of correlation coefficients as a plurality of second output values;

wherein if the third input value is the current signal component, the plurality of second output values are used as a plurality of current correlation coefficients; and if the third input value is the vibration signal component, taking a plurality of second output values as the vibration correlation coefficient.

According to the technical scheme provided by the embodiment of the application, a kernel principal component analysis dimensionality reduction algorithm is used for carrying out feature level data fusion on the current feature parameters and the vibration feature parameters; obtaining a characteristic fusion parameter matrix; reducing the dimension of the characteristic fusion parameter matrix to obtain a plurality of principal elements;

taking the current characteristic parameter or the vibration characteristic parameter as a plurality of fifth input values;

the specific steps of performing data fusion and dimensionality reduction by using the kernel principal component analysis dimensionality reduction algorithm comprise:

acquiring a plurality of fifth input values;

performing matrix combination operation on the fifth input values to obtain a feature matrix;

carrying out standardized operation on the feature matrix to obtain a standardized feature matrix;

calculating a kernel matrix according to the standardized feature matrix to obtain a feature fusion parameter matrix;

centralizing the feature fusion parameter matrix to obtain a centralized feature fusion parameter matrix;

calculating to obtain a plurality of eigenvalues of the centralized characteristic fusion parameter matrix and a plurality of eigenvectors corresponding to the eigenvalues, and arranging the eigenvalues in a descending order;

sequentially calculating the accumulated contribution rates of the characteristic values according to a descending order to obtain a plurality of accumulated contribution rates;

and selecting a plurality of eigenvectors corresponding to a plurality of eigenvalues which are larger than the fifth set value and have the least required eigenvalue and the cumulative contribution rate as a plurality of principal elements.

According to the technical scheme provided by the embodiment of the application, the current signal component and the vibration signal component both comprise:

time domain characteristic parameters, including:

mean, root mean square value, variance, square root amplitude, peak, skewness, kurtosis, form factor, peak factor, pulse factor, and margin factor;

frequency domain characteristic parameters, including:

mean, root mean square, and standard deviation frequencies;

entropy domain feature parameters, including:

time domain entropy and frequency domain entropy.

Another aspect of the present application provides a method for predicting a lifetime of a speed reducer, including: the method for extracting a degradation index of a speed reducer according to claim 1, further comprising: the method for predicting the service life of the harmonic reducer by using the longicorn whisker algorithm model and the neural network model comprises the following steps:

acquiring a neural network model, and setting an initial weight and an initial threshold;

establishing a longicorn stigma algorithm model;

optimizing the initial weight and the initial threshold by using a longicorn algorithm model; obtaining an optimized weight and an optimized threshold;

the optimized weight and the optimized threshold are used as the weight and the threshold of the neural network model again; obtaining an optimized neural network model;

obtaining a sample set of harmonic reducer degradation indicators, wherein the sample set comprises a plurality of sample points, and each sample point comprises a plurality of principal elements;

dividing the sample set into a training set and a test set;

training the optimized neural network model by taking a plurality of principal elements of each sample point in the training set as input and taking the residual service life of a training sample as output; obtaining a trained optimized neural network model;

taking a plurality of pivot elements of each sample point in the test set as input, and testing the trained optimized neural network model; and obtaining the residual service life of the test sample.

According to the technical scheme provided by the embodiment of the application, the specific process of optimizing the initial weight and the initial threshold by using the longicorn whisker algorithm comprises the following steps:

establishing a longicorn whisker algorithm model, which comprises the following steps:

setting a decreasing factor and a total iteration number; determining a search direction; setting a step factor;

creating a left-right whisker space coordinate simulating a search behavior of a longicorn left-right whisker;

setting a fitness function and a position updating formula;

inputting the initial weight and the initial threshold value into the longicorn whisker algorithm model, and performing iterative operation, wherein the specific iterative process comprises the following steps:

inputting a weight value and a threshold value into the longicorn whisker algorithm model as an initial position;

performing a plurality of iterations comprising:

judging the odor intensity of the left beard and the right beard;

move to the side where the odor is stronger: updating the individual position according to the position updating formula;

judging whether an ending condition is reached; the end conditions are as follows: the total iteration number has been reached or the fitness function is less than a sixth set value;

and outputting the optimized weight and the optimized threshold after the end condition is reached.

According to the technical scheme provided by the embodiment of the application, the residual service life of the test sample comprises a residual service life value of each sample point in the test set;

obtaining a degradation index of a sample point according to the residual service life value;

taking the stage where the sampling point with the degradation index smaller than the first degradation value is as a healthy stage;

taking the stage where the sampling point with the degradation index larger than the first degradation value and smaller than the second degradation value is as a slight degradation stage;

and taking the stage where the sampling point with the degradation index larger than the second degradation value is positioned as a serious degradation stage.

The beneficial effect of this application lies in:

the current signal and the vibration signal in the harmonic reducer are extracted simultaneously, and are preprocessed respectively, and modal decomposition and characteristic extraction are carried out respectively; obtaining a plurality of current components and vibration components, and calculating the correlation between the current components and the correlation between the vibration components; then, correlation screening is carried out, and feature level data fusion and dimension reduction are carried out; obtaining a plurality of principal elements; and obtaining a degradation index curve by calculating the relationship between the principal element and the contribution rate of the principal element. Meanwhile, the current information and the vibration information are combined and the signals are preprocessed, so that the influence of noise signals on fault characteristics can be reduced, and the fault missing rate and the fault false rate are effectively reduced; modal decomposition feature extraction is carried out, coefficient fusion dimensionality reduction is carried out by utilizing the correlation among all components, more than 85% of fault information contained in signals in the early stage of operation of the equipment can be effectively extracted, and the accuracy of analysis results is improved.

Drawings

Other features, objects and advantages of the present application will become more apparent upon reading of the following detailed description of non-limiting embodiments thereof, made with reference to the accompanying drawings in which:

fig. 1 is a schematic diagram of a method for extracting a degradation indicator of a speed reducer according to this embodiment;

FIG. 2 is a schematic diagram illustrating a life prediction method for a speed reducer according to an embodiment of the present invention;

FIG. 3 is a graph of the variation trend of the degradation index of the harmonic reducer;

FIG. 4 is a waveform diagram of an initial current signal;

FIG. 5 is a waveform diagram of an initial vibration signal;

FIG. 6 is a waveform diagram of a first current signal;

FIG. 7 is a waveform diagram of a first vibration signal;

FIG. 8 is a graph of the variation of the signal components after filtering of the current data;

FIG. 9 is a graph of the variation of the signal components after vibration data filtering;

fig. 10 is a time domain waveform at 100 sample points after the current signal MCKD is processed;

fig. 11 is a time domain waveform at 600 sample points after the current signal MCKD is processed;

FIG. 12 is a time domain waveform at 1060 sample points after processing of current signal MCKD;

FIG. 13 is a 16-dimensional characteristic time domain waveform diagram for the # 3 reducer;

FIG. 14 is a schematic diagram of the correlation between 16 features and signal components;

FIG. 15 is a comparison graph of the prediction result and the true value of the training sample of the BAS-BP prediction model;

FIG. 16 is a graph comparing the predicted results and the actual values of the test samples of the BAS-BP prediction model;

FIG. 17 is a diagram of the prediction error of the BAS-BP prediction model.

Detailed Description

The present application will be described in further detail with reference to the following drawings and examples. It is to be understood that the specific embodiments described herein are merely illustrative of the relevant invention and not restrictive of the invention. It should be noted that, for convenience of description, only the portions related to the present invention are shown in the drawings.

It should be noted that the embodiments and features of the embodiments in the present application may be combined with each other without conflict. The present application will be described in detail below with reference to the accompanying drawings in conjunction with embodiments.

Example 1

Referring to fig. 1, a schematic diagram of a method for extracting a degradation indicator of a speed reducer according to this embodiment includes:

s1: acquiring an initial current signal and an initial vibration signal;

s2: respectively carrying out signal preprocessing on the initial current signal and the initial vibration signal to obtain a first current signal and a first vibration signal;

s3: respectively carrying out modal decomposition and feature extraction on the first current signal and the first vibration signal; obtaining a plurality of current signal components and a plurality of vibration signal components;

s4: optionally selecting one of said current signal components as a reference current component; optionally one of said vibration signal components as a reference vibration component;

s5: respectively calculating the correlation of each current signal component which is not taken as the reference current component to obtain a plurality of current correlation coefficients;

s6: eliminating current signal components with current correlation coefficients larger than a first set value, and taking the residual current signal components as current characteristic parameters; eliminating vibration signal components with vibration correlation coefficients larger than a second set value, and taking the residual vibration signal components as vibration characteristic parameters;

s7: performing characteristic level data fusion on the current characteristic parameters and the vibration characteristic parameters to obtain a characteristic fusion parameter matrix;

s8: reducing the dimension of the characteristic fusion parameter matrix to obtain a plurality of principal elements;

s9: respectively calculating the contribution rates of the multiple principal elements to fault information to obtain multiple contribution rates;

s10: and multiplying the plurality of principal elements by the contribution rates of the principal elements, and summing to obtain a degradation index curve.

In some embodiments, the retarder is a harmonic retarder. The application takes a harmonic reducer as an example. Meanwhile, the method can be further applied to other speed reducers through non-creative modification.

In some embodiments, a median filtering method is used to filter the acquired original signal, so as to obtain a filtered current signal and a filtered vibration signal that do not contain high-frequency sampling noise after filtering. After the original motor stator current signal and the original vibration signal are collected, the high-frequency noise component in the signal is eliminated by using median filtering. The method can ensure that the edges of the signals are not blurred while noise is filtered. And a high-quality source signal is provided for subsequent analysis, and the effectiveness of the extraction of the degradation characteristics is ensured.

In some embodiments, referring to FIG. 3, to facilitate observing the degradation index as a function of the degree of degradation of the harmonic reducer, a normalization (Min-max normalization) process is performed on the degradation index data over the full life cycle. Different degradation stages of the harmonic reducer under the full life cycle are divided, and the harmonic reducer totally undergoes three stages from a normal state to failure: healthy stage, mild degenerative and severe degenerative stage. The degradation indexes of the harmonic reducer in different time periods can be displayed most visually, and observation is facilitated. The number of the sample points represents the number of sampling points from the beginning of the harmonic reducer to the whole time line; and after sampling, taking the sampling point as a sample point, wherein the number of the sampling points is equal to that of the sample point. The degradation indicator is the ratio of the current elapsed time to the total service life.

Specifically, the degradation indicator curve is normalized according to formula (one);

(one);

wherein the content of the first and second substances,

maximum boundary for normalized data;

a minimum boundary for normalized data;

in order to normalize the values of the data before they are normalized,

as a value after normalization.

Further, selecting the current signal component with the highest kurtosis as a reference current component; selecting the vibration signal component with the highest kurtosis as a reference vibration component.

In some embodiments, the most kurtosis signal component is selected as the reference component; the calculation complexity can be reduced to the maximum extent, and the contained fault characteristics are reserved; the calculation result is closer to the real situation.

Further, the signal preprocessing comprises: filtering and denoising; respectively performing noise reduction processing on the initial current signal and the initial vibration signal by using a correlation kurtosis deconvolution Method (MCKD);

substituting the initial current signal or the initial vibration signal as a first input signal into the correlation kurtosis deconvolution method, the method specifically comprising:

s3.8: taking the second filter function as the first filter function, and returning to S3.3;

s3.9: taking the filtered signal as a first output signal;

In some embodiments, the initial current signal and the initial vibration signal are separately denoised using a correlated kurtosis deconvolution Method (MCKD).

The related kurtosis has the characteristic that the larger the impact is, the larger the kurtosis value is, so that the intensity of the fault impact signal in the current signal can be accurately reflected. Meanwhile, the correlation kurtosis can also keep the characteristics of a correlation function, and can extract the impact signal of a specific period.

Therefore, the related kurtosis deconvolution method can fully consider the periodicity of impact components in the signals, and highlight fault pulses in the signals by taking the related kurtosis as an optimization object to obtain noise reduction current signals and noise reduction vibration signals.

In some embodiments, the third set point is 0.01; the noise influence can be effectively reduced and the fault pulse in the signal can be highlighted.

In some embodiments, the kurtosis is calculated according to equation (two);

(II);

wherein a (n) is a first input signal,

the first input signal is delayed by M pulse periods, N is the number of data points, L is the length of a filter, T is the pulse period, M is the number of time shift periods, and M =1,2, \ 8230;, M.

Specifically, assuming that the first input signal is a (n), the first input signal may be represented by a fault impact signal x (n), a response h (n) of the system transmission path, and a noise component e (n) as shown in formula (three);

and (III).

Deriving a formula (IV) by which the first input signal is restored to the fault impact signal;

and (IV).

The MCKD algorithm takes an output signal r (n) with the maximum correlation kurtosis after the first input signal a (n) is filtered as an objective function, and the filter also needs to be selected to meet the maximum correlation kurtosis and is obtained according to the technology of the formula (five);

(V);

wherein L is the length of the filter,

the filter is searched for optimally.

Order to

An electronic filter

Can be expressed by the formula (six);

and (VI).

Wherein, the first and the second end of the pipe are connected with each other,

an open root number representing the sum of squares of the components;

mixing k =1,2, \ 8230l; the bands are sequentially entered into formula (five) to obtain formula (six), wherein B,

，

，

both are matrix representations of the coefficients.

B is obtained by calculation of a formula (VII);

(VII);

calculated by formula (eight);

is composed of

The transposed matrix of (2);

(eight);

calculated by formula (nine);

(nine);

calculated by formula (ten);

and (ten).

Among others, the details of the MCKD algorithm are described in the www.elsevier.com/locate/ymssp website under the name "Maximum corrected Kurtosis reduction and application on gear

the publication of the tooth chip fault detection "belongs to the prior art, and is not described herein.

Further, performing modal decomposition and feature extraction on the first current signal and the first vibration signal respectively by using a complementary set empirical mode decomposition algorithm (CEEMD) to obtain a plurality of current signal components and a plurality of vibration signal components;

s4.1: acquiring a second input signal;

s4.2: adding a plurality of groups of white noises to the second input signal respectively to obtain a plurality of first set signals;

respectively subtracting a plurality of groups of white noise from the second input signal to obtain a plurality of second set signals;

s4.5: respectively calculating the average value of the first maximum envelope curve and the first minimum envelope curve of each first set signal to obtain a first average line of each first set signal;

s4.6: respectively calculating the difference value of each first set signal and the corresponding first average line to obtain a plurality of first difference value functions;

s4.7: and (3) judging: s4.9 is carried out when the plurality of first difference functions and the plurality of second difference functions all meet the condition of intrinsic mode functions; otherwise, carrying out the next step;

s4.9: taking a plurality of said first difference functions as a plurality of first eigenmode functions and a plurality of said second difference functions as a plurality of second eigenmode functions;

In particular, the spline curve refers to a curve given a set of control points, the approximate shape of which is controlled by these points; the method can be generally divided into two types, namely an interpolation spline and an approximation spline, wherein the interpolation spline is generally used for designing digital drawing or animation, and the approximation spline is generally used for constructing the surface of an object.

In particular, the noise-reduced current signal and the noise-reduced vibration generally exhibit non-stationary, non-linear characteristics, and the fault features are difficult to extract. The CEEMD algorithm can be used for carrying out smooth decomposition processing on the non-smooth and non-linear signals, and characteristic information of data is kept in the decomposition process. Meanwhile, the method has the self-adaptive characteristic, so that the essential characteristics of faults in the signals can be reflected better. By adding paired positive and negative white noises to the original signal, the residual noise is eliminated, and the decomposition precision is improved.

Specifically, the eigenmode function conditions include: the difference between the extreme points and the zero points is not more than 1; and the mean value of the maximum envelope curve and the minimum envelope curve is 0.

In the scheme provided by the application, the mean value of the first maximum envelope curve and the first minimum envelope curve is 0; the average value of the second maximum envelope curve and the second minimum envelope curve is 0.

In some embodiments, 10 sets of white noise are added to the second input signal, and 10 sets of white noise are subtracted from the second input signal. Wherein the first component parameter comprises 10 components and the second component parameter comprises 10 components; a set of component parameters has 20 components. The fourth setting is 8; therefore, the CEEMD algorithm has 8 component parameters obtained by performing multiple cycles.

In some embodiments, the first and second sets of signals are calculated according to equation (eleven);

(eleven);

wherein the content of the first and second substances,

is a first set of signals that is to be transmitted,

is a second set of signals;

is a second input signal;

is white noise.

In some embodiments, the first average and the second average are calculated according to equation (twelve);

(twelve);

wherein the content of the first and second substances,

as the first maximum packetThe envelope or the second maxima envelope,

is a first minimum value envelope line or a second minimum value envelope line;

is the first average line or the second average line.

In some embodiments, modal decomposition is performed to obtain a first component parameter and a second component parameter, represented by formula (thirteen);

(thirteen);

wherein m represents the number of the first component parameters and the number of the second component parameters obtained after modal decomposition;

a component parameter representing the jth added white noise obtained by the CEEMD algorithm modal decomposition;

represents the j-th white noise-subtracted component parameter obtained by the modal decomposition of the CEEMD algorithm.

In some embodiments, each set of component parameters is subjected to a set averaging process according to formula (fourteen);

(fourteen);

where M represents the number of times the algorithm is repeated.

Further, the Copula function is used for calculating the correlation of other current signal components to the reference current component and the correlation of other vibration signal components to the reference vibration component;

a plurality of the current signal components that are not the reference current component or a plurality of the vibration signal components that are not the reference vibration component are taken as a plurality of third input values; taking the reference current component or the reference vibration component as a fourth input value;

s5.1: acquiring a plurality of third input values and fourth input values;

s5.2: selecting one of the third input values which is not subjected to the kernel distribution estimation value calculation as a first sample; taking the fourth input value as a second sample;

s5.4: calculating a linear correlation parameter of the first variable and the second variable; obtaining a correlation coefficient between the first sample and the second sample;

s5.6: taking a plurality of correlation coefficients as a plurality of second output values;

wherein if the third input value is the current signal component, the plurality of second output values are used as a plurality of current correlation coefficients; and if the third input value is the vibration signal component, taking a plurality of second output values as the vibration correlation coefficients.

Specifically, when the signal components are too many, information redundancy is generated and the calculation complexity is increased, and the Copula function is used for calculating the correlation between the signal components to remove the redundant signal components. And calculating to obtain a signal component which can show the degradation trend of the harmonic speed reducer clearly, and analyzing the signal component as a reference component. The computational complexity can be effectively reduced, and the efficiency of the degradation feature extraction is improved.

In some embodiments, the first set point and the second set point are both set to 0.3; redundant signal components can be eliminated, and signal components containing the degradation trend of the harmonic reducer are reserved; the computational complexity can be effectively reduced, and the efficiency of the degradation feature extraction is improved.

In some embodiments, the kernel distribution estimate

Calculated according to the formula (fifteen);

(fifteen);

wherein, N represents the number of samples; b represents the window width (or bandwidth), which affects the smoothness of the probability density estimation; k (-) is called a kernel function; x is a random variable;

sample data points for random variables (i =1,2, \8230; N).

In some embodiments, the linear correlation parameter is

Calculating according to a formula (sixteen);

(sixteen);

wherein rho is the correlation parameter for measuring the correlation between u and v (-1, 1); u and v represent two variables;

is the inverse of the standard normal distribution function; exp denotes an exponential function with a natural constant as the base.

Further, performing characteristic level data fusion on the current characteristic parameters and the vibration characteristic parameters by using a Kernel Principal Component Analysis (KPCA) algorithm; obtaining a characteristic fusion parameter matrix; reducing the dimension of the characteristic fusion parameter matrix to obtain a plurality of principal elements;

acquiring a plurality of fifth input values;

carrying out standardization operation on the feature matrix to obtain a standardized feature matrix;

and selecting a plurality of eigenvectors corresponding to a plurality of eigenvalues corresponding to the accumulated contribution rate which is greater than the fifth set value and has the least required eigenvalue as a plurality of principal elements.

Specifically, a certain correlation still exists between the signal components after the redundancy is removed, so that the accurate construction of the degradation index is influenced. Therefore, redundant current and vibration signal components are removed by fusing two noise reduction signals through a kernel principal component analysis and dimensionality reduction algorithm (KPCA), and low-dimensional principal elements, namely the number of the principal elements required by representing that the cumulative contribution rate reaches 85%, are extracted from high-dimensional signal components, so that the problem of high calculation complexity caused by signal component redundancy is solved, and the calculation efficiency is improved.

In some embodiments, selecting principal elements with an aggregate contribution rate of 85% can greatly reduce computational complexity while preserving most of the fault signature contained in the current and vibration signals during early operation of the harmonic reducer. The extracted signals can truly reflect the actual degradation condition of the harmonic reducer.

In some embodiments, according to formula (seventeen)Calculating to obtain a feature fusion parameter matrix

；

(seventeen);

wherein the content of the first and second substances,

all columns of the ith sample;

exp represents an exponential function with a natural constant as a base;

represents the open square of the sum of the squares of all elements;

the bandwidth of the Gaussian kernel is shown and is used for controlling the local action range of the Gaussian kernel function.

In some embodiments, the centralized feature fusion parameter matrix

Calculating according to a formula (eighteen);

(eighteen);

is an N × N matrix, each element is

。

In some embodiments, the calculation of the cumulative contribution rate is: calculating the contribution rate of the maximum characteristic value to achieve a first accumulative high efficiency; calculating the contribution rate of the sum of the first two characteristic values to obtain a second accumulated contribution rate; and analogizing in turn, judging that the calculation is stopped when the cumulative contribution rate reaches the fifth set value, and taking the feature vector corresponding to the feature value contained in the latest cumulative contribution rate as a principal element.

In some embodiments, the fifth set point is set to 85%, and the proportion of the fault information included in the screening feature value can be adjusted by setting the fifth set point; the method can be adaptively adjusted according to actual computing power and fault feature extraction requirements, and the application range of the algorithm is enlarged.

Further, the current signal component and the vibration signal component each include:

time domain feature parameters, including:

frequency domain characteristic parameters, including:

average frequency, root mean square frequency and standard deviation frequency;

entropy domain feature parameters, including:

time domain entropy and frequency domain entropy.

In some embodiments, the time domain entropy and the frequency domain entropy are time domain shannon entropy and frequency domain shannon entropy, respectively.

Specifically, the mean value is calculated according to the formula (nineteen);

(nineteen);

the average value of all the collected current signals or vibration signals is represented and is the average intensity of the whole body; the collected current signal is

M =1,2, \8230, n is the number of sampling points, the same applies below.

The root mean square value

Calculating according to a formula (twenty);

(twenty).

The variance is

Calculating according to a formula (twenty-one);

(twenty one).

The square root amplitude

Calculating according to a formula (twenty two);

(twenty-two).

The peak value

Calculating according to a formula (twenty three);

(twenty three);

the skewness

Calculating according to a formula (twenty-four);

(twenty-four).

The kurtosis

Calculating according to a formula (twenty five);

(twenty-five).

The form factor

Calculating according to a formula (twenty-six);

(twenty-six).

The crest factor

Calculating according to a formula (twenty-seven);

(twenty-seven).

The pulse factor

Calculated according to a formula (twenty eight);

(twenty eight).

The margin factor

Calculating according to a formula (twenty-nine);

(twenty-nine).

The average frequency

Calculating according to a formula (thirty);

(thirty);

representing the frequency domain signal (k =1,2, \8230; m, m = n/2 for the signal length).

The root mean square frequency

Calculated according to the formula (thirty-one);

(thirty-one).

The standard deviation frequency

Calculating according to a formula (thirty-two);

(thirty-two).

The time domain Shannon entropy

Calculated according to the formula (thirty-three)；

(thirty-three);

wherein the content of the first and second substances,

representing time domain signals

The value probability of (2).

The frequency domain Shannon entropy

Calculating according to a formula (thirty-four);

(thirty-four);

wherein the content of the first and second substances,

representing frequency domain signals

The value probability of (2).

Example 2

Referring to fig. 2, a schematic view of a life prediction method of a speed reducer provided in this embodiment includes: the method for extracting a degradation indicator of a speed reducer according to claim 1, further comprising: the service life prediction of the harmonic reducer is carried out by utilizing a longicorn whisker algorithm model and a neural network model (BAS-BP prediction model), and the method comprises the following steps:

s8.1: acquiring a neural network model, and setting an initial weight and an initial threshold;

s8.2: establishing a longicorn stigma algorithm model;

s8.3: optimizing the initial weight and the initial threshold by using a longicorn whisker algorithm model; obtaining an optimized weight and an optimized threshold;

s8.4: the optimized weight and the optimized threshold are used as the weight and the threshold of the neural network model again; obtaining an optimized neural network model;

s8.5: obtaining a sample set of harmonic reducer degradation indicators, wherein the sample set comprises a plurality of sample points, and each sample point comprises a plurality of principal elements;

s8.6: dividing the sample set into a training set and a test set;

s8.7: training the optimized neural network model by taking a plurality of principal elements of each sample point in the training set as input and taking the residual service life of a training sample as output; obtaining a trained optimized neural network model;

s8.8: taking a plurality of pivot elements of each sample point in the test set as input, and testing the trained optimized neural network model; and obtaining the residual service life of the test sample.

In some embodiments, the proportion of the sample set divided into the training set and the test set can be set according to experimental requirements, and is set to be 7. The neural network model is trained through a large number of sample points, so that the prediction result of the model is closer to the actual result, and the accuracy of data is improved.

Further, the specific process of optimizing the initial weight and the initial threshold by using the longicorn whisker algorithm includes:

setting a fitness function and a position updating formula;

performing a plurality of iterations comprising:

judging the strength of left and right fibrous odor;

In some embodiments, the weight of the neural network model and the number k of the threshold values are calculated according to a formula (thirty-five);

(thirty-five);

wherein the content of the first and second substances,

the number of nodes of an input layer;

number of hidden layer nodes;

is the number of output layer nodes.

In some embodiments, the decrement factor eta in the longicorn algorithm model is set to 0.95 and the total number of iterations n is set to 30.

In some embodiments, the process of determining the search direction in the longicorn whisker algorithm model comprises:

creating a random orientation vector of the longicorn whiskers and carrying out normalization processing, wherein the random orientation vector

Calculating according to a formula (thirty-six);

(thirty-six);

where rands (·) is a random function and k is the spatial dimension.

In some embodiments, the step-size factor in the longicorn whisker algorithm model is calculated according to formula (thirty-seven);

(thirty-seven);

wherein the factor is decreased

In the range of [0,1]，

Is the current iteration number, and n is the total iteration number.

The step length of the iteration is the current step length;

the last iteration step size.

In some embodiments, the left and right whisker spatial coordinates that mimic the longicorn left and right whisker search behavior are created according to the formula (thirty-eight):

(thirty-eight);

wherein in the formula

Representing the position coordinates of the longicorn right-hand hair at the t-th iteration,

representing the position coordinates of the longicorn left hair at the t-th iteration,

representing longicorn in the t-th iterationThe coordinates of the center of mass,

indicating the distance between the two whiskers.

In some embodiments, the left and right whisker odor intensity is judged according to a fitness function formula (thirty-nine);

(thirty-nine);

wherein the content of the first and second substances,

the odor intensity of the left palpus is shown,

the strength of the smell of the right beard is higher;

outputting a left whisker output value of the ith sample model;

outputting a value for the ith sample model right whisker;

the actual value of the ith sample at the left whisker;

is the actual value of the ith sample at the right whisker.

In some embodiments, the individual locations are updated according to formula (forty);

(forty);

wherein the content of the first and second substances,

the spatial position of the longicorn at the time of the t +1 time of searching,

for the spatial position of the longicorn at the time of the t-th search,

the step size of the move is detected for the t-th longicorn,

is a function of the sign.

In some embodiments, the sixth set point is 0.001, and the sixth set point can be set small enough to approach the global optimal solution at the end of the iteration; the calculation efficiency of the model is greatly accelerated.

Further, the remaining service life of the test sample comprises a remaining service life value of each sample point in the test set;

taking the stage where the sampling point with the degradation index larger than the first degradation value and smaller than the second degradation value is positioned as a slight degradation stage;

In some embodiments, referring to fig. 3, the first degradation value is 0.45 and the second degradation value is 0.75. The number of the sample points represents the number of sampling points on the whole time line from the beginning of the operation of the harmonic reducer to the current state. The degradation indicator is the ratio of the current elapsed time to the total service life.

Specifically, a small amount of noise signals are generated between parts in a running-in stage of the harmonic reducer in the initial stage of operation, so that the initial degradation index is large; at the end of the severe degradation stage, the wear of parts severely generates a fault signal which is difficult to control, so that the degradation index is smaller. The harmonic reducer is run in before and after the sampling point 200, the degradation index extracted at the moment of normal operation is close to 0, and the current state of the harmonic reducer is good; when the sampling point exceeds 1000, the degradation index is close to 1, which indicates that the harmonic reducer has serious degradation at the moment and is close to failure and rejection.

Specifically, the operating condition of the harmonic reducer is judged according to the three stages, the working state of the harmonic reducer can be predicted, the harmonic reducer can be replaced before the fault occurs, and the fault is avoided.

The specific working process is as follows:

the output torque of the load motor is set to be different multiples of the rated torque of the speed reducer, the rotating speed of the driving motor is set to be 2000 r/min, and the harmonic speed reducer is enabled to run in an accelerated degradation state until the harmonic speed reducer is completely failed. The experimental setting of the input speed-load torque test working conditions of the 3 harmonic reducers is as follows: under the condition that the indoor environment temperature is 25 ℃ and the standard air pressure is adopted, the output end loads corresponding to experimental samples with the output rotating speeds of 2000rpm, no. 1, no. 2 and No. 3 of the servo motors are 12.4 N.m, 15.5 N.m and 18.6 N.m respectively, and are 2 times, 2.5 times and 3 times of the rated torque of the speed reducer respectively. Current signals in the experiment are acquired through a data acquisition card and a current probe, and vibration signals are acquired through an acceleration sensor and a data acquisition card.

The detailed information of all harmonic reducer samples after accelerated life test is shown in table 1, and comprises working conditions, sampling intervals, total number of data samples, duration of accelerated life and failure modes corresponding to the harmonic reducer samples.

TABLE 1 harmonic reducer testing information table

In which, taking the No. 3 harmonic reducer as an example, the degradation index extraction and the life prediction are performed.

Acquiring an initial current signal and an initial vibration signal of a No. 3 harmonic reducer; as shown in fig. 4 and 5.

Respectively carrying out signal preprocessing on the initial current signal and the initial vibration signal to obtain a first current signal and a first vibration signal of the No. 3 harmonic reducer; as shown in fig. 6 and 7.

The signal component (RMS) change curves of the 3# harmonic reducer after filtering processing of the current data and the vibration data in the full life cycle of 0 to 621min (1242 sample points in total) are shown in fig. 8 and fig. 9.

Taking a current signal as an example, the vibration signal is processed in the same way, and the MCKD processing is carried out on 1242 sample points in the whole life cycle of the No. 3 harmonic reducer, so that the periodic pulse component in the current signal is extracted, and the fault impact component in the signal is effectively highlighted. Fig. 10, 11, and 12 are time domain waveforms at the 100 th, 600 th, and 1060 th sample points of three different degradation states after the current signal is filtered and MCKD processed, respectively. As is apparent from the graph, the fault pulse component contained in the current signal after the MCKD processing gradually increases as the degradation degree of the harmonic reducer increases.

and performing modal decomposition and feature extraction on the first current signal and the first vibration signal by using a Complementary Ensemble Empirical Mode Decomposition (CEEMD) algorithm on the MCKD to obtain 16-dimensional features (P1 to P16 respectively represent 16-dimensional features and respectively correspond to 16 signal components) in a time domain, a frequency domain and an entropy domain. The 16-dimensional characteristic time-domain waveform diagram of the No. 3 reducer is shown in FIG. 13.

eliminating the current correlation coefficient larger than a first set value; and rejecting the vibration correlation coefficient larger than the second set value.

Referring to fig. 14, copula functions are applied to compare the correlation of 16-dimensional features with reference signal components (RMS), 16-dimensional features corresponding to 16 signal components, respectively.

If the correlation is too high, the redundancy of the feature dimension is larger, and the feature dimension is rejected. The application selects 0.3 as the first set value and the second set value for analysis. The analysis results are shown in fig. 14, and P1 to P16 represent 16-dimensional features, respectively. The correlation between the characteristic dimensions P3, P4, P8 and the characteristic dimension P2 is relatively large (note that P2 is no longer compared with itself and therefore the correlation is set to 0; in order to prevent rejection according to the first setting or the second setting).

Performing characteristic level data fusion on the current characteristic parameters and the vibration characteristic parameters to obtain characteristic fusion parameters;

reducing the dimension of the feature fusion parameters; obtaining a plurality of principal elements;

performing feature level data fusion and dimensionality reduction by adopting a KPCA (kernel principal component analysis) method, and performing nuclear parameter analysis; the optimal state is reached when the kernel parameter c =17, and the number of principal elements at this time is 4, which is the minimum number of principal elements required for the cumulative contribution rate to reach 85%.

And multiplying the first four principal elements by the corresponding contribution rates respectively, and then summing the four principal elements to obtain the degradation index. In order to facilitate observation of the change of the degradation index with the degradation degree of the harmonic reducer, the degradation index under the full life cycle is subjected to mean value standardization, and the result is shown in fig. 3.

Inputting the obtained four degeneration principal elements into a BAS-BP prediction model (the optimized neural network model optimized by the Tianniu algorithm), outputting the residual service life, and completing the life prediction.

Referring to table 2, principal component (KPCA) values and actual lifetime (RUL) values of a portion of the samples during the degradation process are obtained.

Table 2 pivot values and actual life values of partial samples

Wherein, the network model training parameters are shown in table 3.

TABLE 3 network model training parameter Table

The comparison of the results of the BAS-BP prediction model with the actual life and the error are shown in fig. 15, 16, and 17.

And calculating the goodness of fit respectively input by the result obtained by the BAS-BP model and the actual result. The detailed calculation process and principle of goodness of fit belong to the prior art, and are not described again.

Fig. 15 is a comparison of prediction results of the BAS-BP model training set, the horizontal axis is a training sample, the vertical axis is an output value, and the goodness of fit R =0.99032; fig. 16 is a comparison of prediction results of the test set of the BAS-BP model, where the horizontal axis represents test samples, the vertical axis represents output values, and the goodness of fit R =0.99373; FIG. 17 shows the prediction error of the BAS-BP model, with the horizontal axis showing the sample and the vertical axis showing the error value. The actual life value and the predicted value have a high degree of agreement with each other except for 20 sample points at which operation is started and failure is imminent. There are also only a few errors between the two, and in most cases very small errors. The BAS-BP prediction model can accurately predict the service life of the harmonic reducer in most of time.

The above description is only a preferred embodiment of the application and is illustrative of the principles of the technology employed. It will be appreciated by a person skilled in the art that the scope of the invention according to the present application is not limited to the specific combination of the above-mentioned features, but also covers other embodiments where any combination of the above-mentioned features or their equivalents is made without departing from the inventive concept. For example, the above features may be replaced with (but not limited to) features having similar functions disclosed in the present application.

Claims

1. A method for extracting degradation indexes of a speed reducer is characterized by comprising the following steps:

acquiring an initial current signal and an initial vibration signal;

respectively carrying out modal decomposition and feature extraction on the first current signal and the first vibration signal; obtaining a plurality of current signal components and a plurality of vibration signal components;

optionally one of said current signal components as a reference current component; optionally selecting one of said vibration signal components as a reference vibration component;

eliminating current signal components with current correlation coefficients larger than a first set value, and taking the residual current signal components as current characteristic parameters; eliminating vibration signal components with the vibration correlation coefficient larger than a second set value, and taking the residual vibration signal components as vibration characteristic parameters;

2. The method according to claim 1, wherein the current signal component with the highest kurtosis is selected as a reference current component; selecting the vibration signal component with the highest kurtosis as a reference vibration component.

3. The method for extracting a degradation indicator of a speed reducer according to claim 1, wherein the signal preprocessing comprises: filtering and denoising; respectively carrying out noise reduction processing on the initial current signal and the initial vibration signal by utilizing a correlation kurtosis deconvolution method;

s3.5: calculating a difference value of the correlation kurtosis of the first input signal and the correlation kurtosis of the filtering signal to obtain a correlation kurtosis difference value;

s3.9: taking the filtered signal as a first output signal;

if the initial current signal is used as a first input signal, the first output signal is used as the first current signal; and if the initial vibration signal is used as a first input signal, taking the first output signal as the first vibration signal.

4. The method for extracting the degradation index of the speed reducer according to claim 1, wherein a complementary set empirical mode decomposition algorithm is used for carrying out mode decomposition and feature extraction on the first current signal and the first vibration signal respectively to obtain a plurality of current signal components and a plurality of vibration signal components;

s4.1: acquiring a second input signal;

s4.4: respectively connecting the plurality of first maximum points and the plurality of first minimum points of each first set signal by spline curves to obtain a first maximum envelope line and a first minimum envelope line of each first set signal;

respectively calculating the average value of the second maximum envelope curve and the second minimum envelope curve of each second set signal to obtain a second average line of each second set signal;

s4.10: calculating the component of each first intrinsic mode function to obtain a plurality of first component parameters; calculating the component of each second intrinsic mode function to obtain a plurality of second component parameters; taking the first component parameter and the second component parameter as a set of component parameters;

if the first current signal component is taken as the second input signal, taking the average value of a plurality of signal components as a plurality of current signal components; and if the first vibration signal component is taken as the second input signal, taking the average value of the plurality of signal components as a plurality of vibration signal components.

5. The method for extracting the degradation index of the speed reducer according to claim 1, wherein the Copula function is used for calculating the correlation of other current signal components to the reference current component and the correlation of other vibration signal components to the reference vibration component;

s5.1: acquiring a plurality of third input values and fourth input values;

6. The extraction method of the degradation index of the speed reducer according to claim 1, wherein a kernel principal component analysis dimensionality reduction algorithm is used for performing feature level data fusion on the current feature parameter and the vibration feature parameter; obtaining a characteristic fusion parameter matrix; reducing the dimension of the characteristic fusion parameter matrix to obtain a plurality of principal elements;

acquiring a plurality of fifth input values;

7. The method according to claim 1, wherein the current signal component and the vibration signal component each include:

time domain characteristic parameters, including: mean, root mean square value, variance, square root amplitude, peak, skewness, kurtosis, form factor, peak factor, pulse factor, and margin factor;

frequency domain characteristic parameters, including: average frequency, root mean square frequency and standard deviation frequency;

entropy domain feature parameters, including: time domain entropy and frequency domain entropy.

8. A method for predicting a life of a speed reducer, comprising: the method for extracting a degradation indicator of a speed reducer according to claim 1, further comprising: the method for predicting the service life of the harmonic reducer by using the longicorn whisker algorithm model and the neural network model comprises the following steps:

establishing a longicorn stigma algorithm model;

optimizing the initial weight and the initial threshold by using a longicorn whisker algorithm model; obtaining an optimized weight and an optimized threshold;

dividing the sample set into a training set and a test set;

taking a plurality of pivot elements of each sample point in the training set as input, and taking the residual service life of a training sample as output to train the optimized neural network model; obtaining a trained optimized neural network model;

9. The method for predicting the life of the speed reducer according to claim 8, wherein the specific process of optimizing the initial weight and the initial threshold by using the longicorn whisker algorithm comprises the following steps:

setting a decreasing factor and a total iteration number; determining a searching direction; setting a step factor;

creating a left and right whisker space coordinate simulating the searching behavior of the left and right whiskers of the longicorn;

setting a fitness function and a position updating formula;

performing a plurality of iterations comprising:

judging the strength of left and right fibrous odor;

judging whether an ending condition is reached; the end conditions are as follows: the total iteration number is reached or the fitness function is smaller than a sixth set value;

10. The method of claim 8, wherein the test sample remaining useful life comprises a remaining useful life value for each sample point in the test set;

taking the stage where the sampling point with the degradation index smaller than the first degradation value is located as a health stage;