A Novel Shearer Cutting State Recognition Method Based on Improved Variational Mode Decomposition and LSSVM with Acoustic Signals

,


Introduction
In the coal mining face, the shearer plays an indispensable role, and the continuous, safe, and reliable operation of the shearer is an important guarantee to achieve high production and efficiency of coal mining. In order to improve the intelligent level of the shearer, it is necessary to realize the accurate identification of the shearer cutting state in the working process and to adjust the height of the cutting drum and the traction speed of the body adaptively [1,2]. However, how to accurately and quickly identify the current geological conditions of coal seam is still a technical problem in the field of fully mechanized mining technology and also a major technical bottleneck restricting the intelligent level of shearers.
In the actual coal production process, experienced shearer operators often judge the coal-rock cutting state according to the sound of the shearer cutting coal and manually adjust the traction speed and drum height. e research shows that the sound produced by the interaction between the drum and coal-rock in the cutting process of the shearer contains abundant information about the cutting state. erefore, the recognition of the shearer cutting state based on sound provides a new idea for the research in this field.
In recent years, some acoustic-based fault diagnosis techniques have been developed [3][4][5]. Yao et al. [6] established a convolutional neural network for gear fault diagnosis by using the time and frequency domain sound signals. Van Hecke et al. [7] used the heterodyne-based frequency reduction technique and spectral averaging to process acoustic emission signals and extract condition indicators for bearing fault diagnosis. In [8], an intelligent fault diagnosis architecture for the fused magnesia furnace based on the analysis of the acoustic spectrum submanifolds was proposed. Kumar et al. [9] proposed a bearing fault diagnosis method based on the statistical features of sound signal and Bayes classifier. Some researchers also used array measurement technology in acoustic-based diagnosis. Hou et al. [10,11] developed one near-field acoustical holography-(NAH-) based fault diagnosis method through the combination of NAH technology and blocking feature extraction. Lu et al. [12,13] focused on applying distribution information of sound field to gearbox fault diagnosis and presented a new diagnosis method for gearbox based on NAH and spatial distribution features of sound field. Although the acoustic-based fault diagnosis techniques have been successfully applied to many rotating parts, such as gears and bearings, there is still little research in the field of coal-rock cutting state recognition.
In the actual coal mining process, the sound signal collected by the sound sensor belongs to a coupling signal, including the sound emitted by other mechanical and electrical equipment and personnel, in addition to the sound signal of cutting coal-rock. However, the coupling relationship between different sound sources is discrepant, and the decoupling process involves many parameters that are hard to determine. It is extremely difficult to describe the sound signal by using the mathematical modeling method. erefore, starting with the acquisition of original sound signal, this paper extracts the characteristic information that can represent the cutting state of the shearer on the basis of signal decomposition. Common signal processing methods include empirical mode decomposition (EMD), ensemble empirical mode decomposition (EEMD), local mean decomposition (LMD), and variational mode decomposition (VMD). Among these methods, EMD, firstly proposed by Huang et al. in 1998 [14], has been widely used and achieved significant results [15][16][17]. It is a self-adaptive time-frequency decomposition method and can self-adaptively decompose any complex signal into a series of intrinsic mode functions (IMFs). However, in the process of EMD decomposition, there is often a phenomenon of mode mixing, which leads to the incomplete separation of high-frequency and low-frequency signals in IMFs [18,19]. In order to solve the above problems, Huang et al. [20] proposed EEMD on the basis of EMD algorithm in 2003. By adding random noise into the original signal, EEMD can eliminate the influence of noise on the decomposition process to a certain extent and achieve complete mode separation. However, the decomposition precision of EEMD is affected by the added noise, which will lead to overdecomposition or insufficient decomposition, thus generating false IMF components and endpoint effect [21,22]. LMD was proposed by Smith in 2005 and has been firstly applied in the process of electroencephalogram [23], and it then attracted many scholars' attention in the fault diagnosis of rotating machinery [24,25]. However, mode mixing problem is still occurring in the strong noise environment. Recently, a novel self-adaptive signal analysis method of VMD was proposed by Dragomiretskiy and Zosso in 2014 [26]. Compared with EMD and LMD, VMD can effectively avoid the false component and mode mixing problems [27]. In [28], the kernel-based mutual information (KEMI) fitness function was presented to find the optimizing parameters of VMD for the purpose of easy identification of single and multiple defects of bearing. In [29], a new fault diagnosis method is proposed by using particle swarm optimization variational mode decomposition (PSO-VMD) to distinguish different work states of hydraulic pumps. However, the decomposition accuracy of VMD is affected by the penalty factor α and the number of modes K. Too large or too small K will still lead to overdecomposition or insufficient decomposition. erefore, an appropriate intelligent algorithm is needed to optimize these parameters.
Based on the foregoing descriptions, a novel swarm intelligence algorithm of gravitational search algorithm (GSA) with Levy flight strategy is proposed to optimize the parameters of VMD. en, the envelope entropy and kurtosis of each IMF component are extracted, and the principal component analysis method (PCA) is introduced to achieve useful feature selection. Finally, a multistate classifier is designed based on least square support vector machine (LSSVM) to identify the cutting state types of the shearer. e rest of this paper can be summarized as follows. Section 2 describes the basic theory of VMD. In Section 3, the improved gravitational search algorithm is proposed and the flowchart is designed. Section 4 presents the parameter optimization of VMD by using improved GSA. Section 5 gives a detailed description for the proposed shearer cutting state recognition method. Some experiments are provided to verify the effectiveness and superiority of the proposed method in Section 6. Finally, the conclusions are summarized in Section 7.

Variational Mode Decomposition
Variational mode decomposition (VMD) uses the alternating direction multiplier method to realize the adaptive separation of the signal frequency domain and each component and has better noise robustness. In VMD, the intrinsic mode function (IMF) is defined as an AM-FM signal, with the expression as follows: where A k (t) is the instantaneous amplitude. In fact, the decomposition process of VMD can be essentially divided into two processes: construction and solution of variational problem.

e Construction of Variational Problem.
If the bandwidth of each "mode" is finite and the central frequency is different, the variational problem can be described as K mode functions. u k (t) can be calculated on the premise that the sum of the modal functions equals to the input signal f(t) so that the sum of the estimated bandwidths of the mode functions is minimized. Specific steps are as follows: Hilbert transformation was performed on each mode function u k (t) to obtain the transformed analytic signal, and the unilateral spectrum can be solved: e center frequency of each mode analytical signal is estimated, and the frequency spectrum of each mode is modulated and distributed to the corresponding fundamental frequency band: e bandwidth of each mode function is estimated; then, the corresponding constraint variational problem model is as follows: where ω k (t) � ω 1 , ω 2 , . . . is the center frequency of u k .

e Solution of Variational Problem
By introducing the penalty factor α and Lagrangian multiplier operator λ(t) into the solution process, the extended Lagrangian expression is as follows: u n+1 k , ω n+1 k , and n+1 λ are updated alternately by alternating directions of multiplication operators. u n+1 k can be calculated as follows: Formula (6) is transformed into frequency domain: us, the updating method of u n+1 k is as follows: Similarly, the updating methods of ω n+1 k and n+1 λ can be obtained: e specific steps of VMD algorithm are as follows: Shock and Vibration 3 Step 1: initialize u 1 k , ω 1 k ,{λ 1 }, and K.
Step 5: determine if the stop condition is satisfied: If so, stop the iteration and obtain K mode components. Otherwise, repeat Step 2-Step 4 until the maximum number of iterations is reached.

Standard Gravitational Search Algorithm.
Gravitational search algorithm (GAS) was proposed by Rashedi in 2009 on the basis of the law of universal gravitation, as shown in Figure 1. D is the dimension of search space, N is the number of particles, and X i is the position of particle i, which can be expressed as follows: rough many experiments, it is found that when solving the gravity in GSA, the square R 2 of the distance between two objects in the law of universal gravitation is changed to R, which is faster and more effective in calculation. erefore, the gravitational force F d ij (t) of particle j acting on particle i in the dth dimension at a certain time t can be calculated by the following formula: where M ai (t) and M pi (t) are the inertial mass of particle i and j, respectively, ε is an infinitesimal constant, R ij (t) is the Euclidean distance between particle i and j, and G(t) is the gravitational constant at time t, which can be calculated as follows: where G 0 is the gravitational constant at the initial time t 0 and is usually set to 100 and T is the number of iterations. In order to increase the randomness and diversity of GSA, the total force acting on particle i in the dth dimension is the sum of forces acting on other N − 1 particles, and the size of F d i (t) is calculated as follows: where rand j is a random constant between 0 and 1. e acceleration a d i (t) of particle i in the dth dimension at time t is calculated as follows: e velocity and position update mechanism of particle i can be expressed as follows: e inertial mass of a particle is related to its current fitness value. e greater the inertial mass of a particle is, the closer the position of the particle is to the optimal value. Meanwhile, the greater the gravitational force on the particle, the lower the moving speed. e inertial mass of the particle can be calculated as follows: where M ii is the inertial mass of the i-th particle, fit(t) is the fitness value of the ith particle at time t, and worst(t) and best(t) are the optimal fitness value and the worst fitness value, respectively. According to foregoing description, the specific implementation steps of GSA are as follows: Step 1: determine the search space, set the population size N, the maximum number of iterations T, and other relevant parameters and initialize the particle velocity v d i (t) and position x d i (t) Step 2: calculate particle fitness value fit(t) and update the gravitational constant G(t), the worst fitness value worst(t), the optimal fitness value best(t), and the particle inertial mass M i (t) Step 3: calculate the total force F i d acting on the ith particle and the acceleration a d i (t) Step 4: update the velocity v d i (t + 1) and position x d i (t + 1) by using equation (16) Step 5: repeat Steps 2-4, until the optimization precision is satisfied or the maximum number of iterations is reached, and then output the optimization result

Levy Flight Strategy.
Levy flight is a non-Gaussian random walk strategy, and the step size obeys Levy distribution, At present, some scholars have proved that the foraging trajectories of many animals and insects in nature conform to Levy distribution, such as the foraging flight trajectories of albatrosses, the intermittent foraging flight trajectories of fruit flies, and the foraging behaviors of bees.
e Levy distribution can be expressed as follows: where s is the random step length, λ 0 is an exponential parameter, and 1 ≤ λ 0 ≤ 3. Because the variance of Levy distribution varies exponentially with time, the search performance is more effective than Brownian random motion in large-scale search. e random flight stride length adopts an alternating mode of frequent short-distance flight and occasional long-distance flight.
e characteristics of Levy flight can increase the population diversity and expand the search range of the algorithm. erefore, using Levy flight path in intelligent optimization algorithm can make it easier to jump out of local optimum. e common formula for calculating the random step length s of Levy flight is as follows: where u and v obey normal distribution,

Improved Gravitational Search Algorithm Based on Levy
Flight. In GSA, the optimal solution of the algorithm is the position of the particle with the largest inertial mass, which can be found by the gravitational force between the particles. With the continuous evolution of the population, individuals move to the optimal solution region. In the later stage of evolution, the clustering phenomenon is more serious, and the pace of individual movement is becoming smaller, so it is easy to fall into the local extreme value, and premature convergence will occur. e updating strategy of Levy flight adopts an alternating mode of frequent short-distance flight and occasional longdistance flight, which can increase the diversity of the population, expand the search range, and help the algorithm to jump out of local optimum. erefore, Levy flight is introduced to update the position and velocity of the particles in GSA, and then an improved GSA (IGSA) is presented in this paper. e specific improvement ideas are as follows: when updating the position and velocity of particles, the particles are sorted according to the fitness values. en, a proportional coefficient θ is determined to divide the sorted particles into two parts of A 1 and A 2 . For the A 1 part with better fitness values, the original updating mechanism of GSA is still adopted and the updating mechanism based on random step size of Levy flight is used to update the position of particles in the A 2 part with poor fitness values. e position updating mechanism based on Levy flight is as follows: where β � 1/t is the step size factor, which will decrease with the increase of iterations. Levy_step is random step size of Levy flight and Levy_step � s. e flowchart of the improved GSA based on Levy flight is shown in Figure 2, and the specific implementation steps are as follows: Step 1: set the population size N, the maximum number of iterations T and θ, and other related parameters; initialize the velocity v d i (t) and position x d i (t) of the particle.
Step 2: calculate the fitness value fit i (t), update the universal gravitation constant G(t), the worst fitness value worst(t), the optimal fitness value best(t), and the inertial mass M i (t).
Step 3: calculate the total force F i d on the particle and its acceleration a d i (t).
Step 4: the particles are sorted according to their fitness values and divided into the better part A 1 and the worse part A 2 according to the proportional coefficient θ.
Step 5: the particles in A 1 are updated by using equation (16) and the particles in A 2 part are updated by using equation (20).
Step 6: calculate the fitness value of the updated particle. If it is better than the previous particle, update the position of the particle; otherwise, continue to use the position of the previous particle.
Step 7: Repeat Step 2-Step 6 until the optimization precision is satisfied or the maximum number of iterations is reached, and output the optimization result.

The Combination of VMD and IGSA
In the decomposition process of VMD, the artificial trial method is mostly adopted to determine the number of modes K, which has a great subjective impact on the decomposition effect. For the penalty factor α, due to its large value range, it is generally selected according to the empirical value. In order to achieve the best decomposition effect, multiple repeated tests are often needed to determine the best parameter combination, with a large amount of calculation and low reliability. erefore, a variational modal decomposition method based on improved gravitational search algorithm (IGSA-VMD) is proposed in this paper. e number of modes K and penalty factor α of VMD are accordingly optimized by using IGSA algorithm to achieve the optimal decomposition effect of signals.
In order to describe the decomposition effect of VMD, the index of orthogonality (IO) is adopted in this paper. e better the orthogonality between the mode components, the smaller the correlation between them and the better the decomposition effect. When the number of modes K and the penalty factor α are not appropriate, the orthogonality Shock and Vibration 5 among the mode components will be affected. erefore, the index of orthogonality (IO) can quantitatively analyze the accuracy of signal decomposition. e smaller the IO value is, the higher the accuracy of signal decomposition and the lower the confusion degree among the modes. e formula for calculating IO value is as follows: where x(t) is the original signal, u(t) is the decomposed mode component, L is the signal length, and H is the total number of sampling points.
In the process of optimizing parameters K and α of VMD by using IGSA algorithm, [K, α] is taken as the position of particle i. e signal is decomposed with this parameter combination and the IO value of decomposed components is Update using equation (16) Update using equation (20) Update particle position Retain current particle e optimal value of the whole optimization process is the minimum IO value and the optimal particle position is the best parameter combination of [K best , α best ]. e flowchart of IGSA-VMD is shown in Figure 3.

Yes No
In order to verify the decomposition performance of proposed IGSA-VMD, a simulated signal is used in this paper, which is a mixed signal composed of four cosine signals with different frequencies.
e mathematical expression and waveform of the signal are shown in Figure 4: To verify the validity of the fitness function, the penalty factor α is selected as 2000 according to the empirical value and the number of modes K is set from 2 to 15, respectively. e simulated signal is decomposed and the IO values are calculated.
e results are shown in Figure 5. From this figure, the IO can reach to the smallest value of 8.1282 × 10 −4 for K � 4, indicating that the simulated signal is composed of four signals with different frequencies, which is consistent with the actual situation. erefore, the optimal number of modes K can be accurately identified according to the evaluation index IO.
In order to further improve the decomposition effect and testify the optimization performance of IGSA, GSA, and IGSA both are used to optimize the parameters of VMD to display the comparisons. e number of particles is set to 30 and the number of maximum iterations is set to 100. e iterative curve of optimal fitness value in the optimization process is shown in Figure 6, and the optimization results are listed in Table 1.
It can be observed that IGSA algorithm obviously has better searching performance than GSA algorithm, with faster convergence speed and higher convergence accuracy. IGSA can obtain excellent parameters combination to achieve better decomposition effect of VMD. In addition, the optimal IO value based on IGSA-VMD can reach to 1.3710 × 10 −4 , which is superior to the empirical approach (8.1282 × 10 −4 ), signifying the necessity of simultaneously optimizing parameters K and α. e decomposition results of simulated signal based on IGSA-VMD are shown in Figure 7. Four IMF components can be accordingly obtained and the central frequencies of the four IMF components are 19.9997 Hz, 49.9992 Hz, 60.0012 Hz, and 100.0082 Hz, respectively, which are basically consistent with the frequency components of the original signal. e simulation results verify the effectiveness of the proposed IGSA-VMD method. Figure 8 illustrates the main point of the proposed shearer cutting state recognition method, which can be summarized as follows.

Signal Decomposition.
IGSA-VMD method is firstly employed to decompose the sound signal to eliminate the effects of distractions, and several IMF components are acquired.

Feature Extraction and Dimensionality Reduction.
In this paper, the envelope entropy and kurtosis are used to extracted features from the sound signals. For a zero-mean signal x with length N, the envelope entropy E p can be calculated as follows: where a is the envelope of demodulation signal obtained by Hilbert transform of signal x. e kurtosis K u of each IMF component can be calculated as follows: where μ and σ are the mean and standard deviation of signal x, respectively. After feature extraction, sixteen eigenvalues can be acquired to describe the sound signal. If too many extracted features are used, it will inevitably cause dimensional disaster and decrease the network performance. Principal component analysis method (PCA) is used to extract sensitive features in this paper.
For a given sample set X � {x 1 , x 2 , . . ., x n }, there are n sample points and each sample point contains p indicators, namely, x i ∈ R p , i � 1, 2, . . ., n; then, Standardize original data: the original data are standardized by using the following equation: where μ j and std j are the mean and standard deviation of the jth dimensional vector in X. Construct the covariance matrix S: the covariance matrix S can be calculated as follows: e eigenvalues (λ 1 , λ 2 , . . ., λ p ) of matrix S and corresponding eigenvectors u i � [u 1i , u 2i , . . . , u pi ], i � 1, 2, . . . , p can be solved. Simultaneously, the eigenvalues are arranged in the descending order. Determine principal component: the larger the contribution rate of the variable index is, the more original data information can be presented by the variable index. e calculation formula of the contribution rate is as follows: e cumulative contribution rate of the first k principal components can be calculated as follows: If the cumulative contribution rate is higher than the preset threshold, the k principal components can be used to evaluate the original data.

Multistate Classifier Based on Least Square Support
Vector Machine. Naturally, after feature extraction, a multistate classifier is designed to identify cutting state types of the shearer. LSSVM has been widely used in many fields due to its high calculation efficiency and generalization ability [30,31]. In LSSVM, radial basis function (RBF) kernel function has universal application and good performance and is used in this paper. e RBF kernel can be expressed as follows: where c is the width of the basis function. For the multiclass identification, several techniques have been developed to solve the multiclass problem on the basis of LSSVM, such as one versus one (OVO), one versus rest (OVR), and hierarchical LSSVM. Because the OVR and hierarchical LSSVM require more time and datasets for training procedures, the OVO is chosen to identify different shearer cutting states. In addition, the parameters of LSSVM, penalty parameter C and the width parameter c, are optimized by particle swarm optimization algorithm (PSO) according to the literature [32]. e parameters of LSSVM are set as C (0.1, 1000) and c (0.001, 10), according to multiple attempts.

Experimental Validation
In this section, an experimental shearer cutting coal-rock system is established to collect cutting sound signal. e experimental setup consists of the shearer, scraper conveyor, hydraulic support, coal-rock specimen and fixed rack, sound sensor, signal conditioner, data acquisition card, and industrial computer, as shown in Figure 9.
In order to better simulate the real coal-rock properties of underground coal mining face, the coal powder, sand, and cement are mixed together with water according to a certain proportion, and four typical coal-rock specimens are developed, as shown in Figure 10. e four specimens have different hardness and gangue mixing ratios. Each specimen has a length of 1800 mm, a height of 1500 mm, and a   thickness of 800 mm. us, five cutting states of F 1 , F 2 , F 3 , F 4 , and F 5 are obtained and listed in Table 2. e sound sensor is installed on the right rocker arm facing the coal-rock specimen, and the sensor belongs to the external bias capacitive transducer with large vibration film lateral entry. e frequency response range is 20-20 kHz and can withstand the maximum sound pressure level of 158 decibels. In order to avoid the impact of rock or coal caving during the experiment, the sensor is fixed inside the metal shell, which transmits sound signals to industrial computer through signal conditioner and data acquisition card for subsequent processing. e time domain waveforms of sound signal with five different cutting states are shown in Figure 11. Because the collision between cutting teeth and coal-rock with different properties can cause complicated dynamic characteristics, resulting in indistinguishable complexity of the measured sound signals, the cutting state of the shearer cannot be directly distinguished through observing the sound signals alone.
In this experiment, the collected data are divided into several nonoverlapping samples. e data length of each sample is 0.5 s and each cutting state contains 60 samples.  us, totally, 300 samples can be generated. Half of the samples are used to create the training set, and the other half are used to generate the testing set. ese details are given in Table 3.
According to the procedure of the proposed cutting state recognition method, the collected signals are firstly decomposed by using IGSA-VMD. e optimal particle position is located at [8,3945] and the corresponding IO value is 0.1791, which indicates that the sound signal can be decomposed into eight IMF components. e results are shown in Figure 12.
en, the envelope entropy and kurtosis are used to extract sensitive feature information of each IMF component to describe different cutting states. us, each sample data can be represented by a feature vector of 16 dimensions, as shown in Table 4, and an feature matrix of 300 × 16 can be obtained.
Obviously, the dimension of the feature vector is too large, and there may be redundant relations, which will cause the dimension disaster, and then affect the recognition accuracy. erefore, PCA is introduced to achieve the dimension reduction processing. e histogram of principal  Coal seam with f = 2 Coal seam with f = 3 Coal seam with gangue Rock Figure 10: Four self-made typical coal-rock specimens.    component contribution rate is shown in Figure 13. When k is equal to 5, the cumulative contribution rate Z k exceeds 90%. en, the new matrix composed of the first 5 principal components is selected to replace the original feature vector, as shown in Table 5.
e new feature vectors are used to train the multiclassifier LSSVM for cutting state recognition of the shearer. In order to verify the recognition performance of the proposed method, LSSVM combined with GSA-VMD and PCA (simplified into GSA − VMD + PCA), LSSVM combined with VMD and PCA (simplified into VMD + PCA), and LSSVM combined with EEMD and PCA are also employed to conduct the pattern identifications. In order to improve the fairness of the comparison results, this paper runs 20 times for each algorithm. Finally, the obtained classification results are shown in Figure 14 and Table 6. Obviously, the following conclusions can be summarized. Firstly, the highest classification accuracy (97.50%-92.50%) can be achieved by the proposed method, which verifies the feasibility and superiority of the proposed method in identifying shearer cutting states. Secondly, the GSA − VMD + PCA and EEMD + PCA have the similar recognition accuracy (about 90%-95%) and the VMD + PCA has the least recognition accuracy (88.33-94.17%). In addition, the proposed method has the smallest standard deviation of the 20 runs, which indicates the stability of the algorithm.
Furthermore, in order to investigate the effectiveness and necessity of the PCA method, five elements are randomly selected from the initial feature vectors and then used to train and test the LSSVM-based classifier. e obtained state identification accuracies are described in Table 7. It can be found that the classification accuracies of four methods are obviously lower than that of the methods combined with PCA in Table 6. e reason is that the randomly selected features do not contain enough abundant characteristic information, resulting in weak generalization ability and   identification performance of the classifier. e comparison results further illustrate the necessity of selecting useful features by using PCA. e comparison results also prove the advantage and superiority of IGSA-VMD in terms of signal decomposition and state feature extraction. In order to verify the effectiveness of the proposed method, the multiclass relevance vector machine-(MRVM-) based method developed in [33] is compared to achieve the cutting state recognition. e obtained recognition accuracies are listed in Table 8. It is obvious that the classification accuracy (90.49%) of MRVM-based method is much lower than that of the proposed method (95.54%), and the stability of proposed method is also superior according the standard deviation. e foregoing comparison results indicate the effectiveness and superiority of proposed method in identifying the shearer cutting state in terms of recognition accuracy.

Conclusions
In this paper, a novel recognition method for the shearer cutting state is proposed based on the combination of IGSA, VMD, PCA, and LSSVM. e position update mechanism based on Levy flight strategy is introduced to enhance the searching performance of GSA, and the IGSA algorithm is presented to optimize the optimal parameter combination of VMD to ensure the decomposition effect. e effectiveness of IGSA-VMD is validated by using the simulated signals. en, the feature to represent the cutting state information is extracted by using envelope entropy and kurtosis of the decomposed IMFs. Furthermore, PCA is introduced to select the sensitive features to enhance the recognition performance of LSSVM. Finally, some experiments are provided, and results indicate that the proposed method has a better state identification ability with higher accuracy.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare no conflict of interest.