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

Authors: Zhongbin Wang; Bin Liang; Lei Si; Kuangwei Tong; Chao Tan
Journal: Shock and Vibration
 (2020)

---

## Abstract

The recognition of shearer cutting state is the key technology to realize the intelligent control of the shearer, which has become a highly difficult subject concerned by the world. This paper takes the sound signal as analytic objects and proposes a novel recognition method based on the combination of variational mode decomposition (VMD), principal component analysis method (PCA), and least square support vector machine (LSSVM). VMD can decompose a signal into various modes by using calculus of variation and effectively avoid the false component and mode mixing problems. On this basis, an improved gravitational search algorithm (IGSA) is designed by using the position update mechanism of Levy flight strategy to find the optimal parameter combination of VMD. Then, the feature extraction is achieved by calculating the envelope entropy and kurtosis of the decomposed intrinsic mode functions (IMFs). To avoid dimensional disasters and reinforce the classification performance, PCA is introduced to choose useful features, and the LSSVM-based classifier is reasonably constructed. Finally, the experimental results indicate that the proposed method is more feasible and superior in the recognition of shearer cutting states.


---

## Body

## 1. 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. The 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. Therefore, 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–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. Therefore, 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–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. Therefore, 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. Then, 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.The rest of this paper can be summarized as follows. Section2 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.


## 2. 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:(1)ukt=Aktcosϕkt,where Akt is the instantaneous amplitude.In fact, the decomposition process of VMD can be essentially divided into two processes: construction and solution of variational problem.


### 2.1. The Construction of Variational Problem

If the bandwidth of each “mode” is finite and the central frequency is different, the variational problem can be described asK mode functions. uk(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 functionuk(t) to obtain the transformed analytic signal, and the unilateral spectrum can be solved:(2)δt+jπt×ukt.The 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:(3)δt+jπt×ukte−jωkt.The bandwidth of each mode function is estimated; then, the corresponding constraint variational problem model is as follows:(4)min∑k∂tδt+jπt×ukte−jωkt2,s.t.∑kuk=ft,where ωkt=ω1,ω2,… is the center frequency of uk.


### 2.2. The 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:(5)Luk,ωk,λ=α∑k∂tδt+jπt∗ukte−jωkt2+ft−∑kukt22+λt,ft−∑kukt.ukn+1, ωkn+1, and λn+1 are updated alternately by alternating directions of multiplication operators. ukn+1 can be calculated as follows:(6)ukn+1=argminα∂tδt+jπt∗ukte−jωkt22+ft−∑iuit+λt222uk∈X.Formula (6) is transformed into frequency domain:(7)u^kn+1=argminu^k,uk∈Xαjω1+sgnω+ωk⋅u^kω+ωk22+f^ω−∑iu^iω+λ^ω222.Thus, the updating method ofu^kn+1 is as follows:(8)u^kn+1=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2.Similarly, the updating methods ofωkn+1 and λn+1 can be obtained:(9)ωkn+1=∫0∞ωu^knω2dω∫0∞u^knω2dω,(10)λ^n+1ω=λ^nω+τf^ω−∑ku^ωn+1ω.The specific steps of VMD algorithm are as follows:Step 1: initializeu^k1,ωk1,{λ^1}, and K.Step 2: updateu^kand ωk according to Formulas (8) and (9).Step 3:k = k + 1. Repeat step 2 until k = K.Step 4: updateλ according to equation (10).Step 5: determine if the stop condition is satisfied:∑ku^kn+1−u^kn22/u^kn22<ε. If so, stop the iteration and obtain K mode components. Otherwise, repeat Step 2–Step 4 until the maximum number of iterations is reached.


## 2.1. The Construction of Variational Problem

If the bandwidth of each “mode” is finite and the central frequency is different, the variational problem can be described asK mode functions. uk(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 functionuk(t) to obtain the transformed analytic signal, and the unilateral spectrum can be solved:(2)δt+jπt×ukt.The 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:(3)δt+jπt×ukte−jωkt.The bandwidth of each mode function is estimated; then, the corresponding constraint variational problem model is as follows:(4)min∑k∂tδt+jπt×ukte−jωkt2,s.t.∑kuk=ft,where ωkt=ω1,ω2,… is the center frequency of uk.


## 2.2. The 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:(5)Luk,ωk,λ=α∑k∂tδt+jπt∗ukte−jωkt2+ft−∑kukt22+λt,ft−∑kukt.ukn+1, ωkn+1, and λn+1 are updated alternately by alternating directions of multiplication operators. ukn+1 can be calculated as follows:(6)ukn+1=argminα∂tδt+jπt∗ukte−jωkt22+ft−∑iuit+λt222uk∈X.Formula (6) is transformed into frequency domain:(7)u^kn+1=argminu^k,uk∈Xαjω1+sgnω+ωk⋅u^kω+ωk22+f^ω−∑iu^iω+λ^ω222.Thus, the updating method ofu^kn+1 is as follows:(8)u^kn+1=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2.Similarly, the updating methods ofωkn+1 and λn+1 can be obtained:(9)ωkn+1=∫0∞ωu^knω2dω∫0∞u^knω2dω,(10)λ^n+1ω=λ^nω+τf^ω−∑ku^ωn+1ω.The specific steps of VMD algorithm are as follows:Step 1: initializeu^k1,ωk1,{λ^1}, and K.Step 2: updateu^kand ωk according to Formulas (8) and (9).Step 3:k = k + 1. Repeat step 2 until k = K.Step 4: updateλ according to equation (10).Step 5: determine if the stop condition is satisfied:∑ku^kn+1−u^kn22/u^kn22<ε. If so, stop the iteration and obtain K mode components. Otherwise, repeat Step 2–Step 4 until the maximum number of iterations is reached.


## 3. Improved Gravitational Search Algorithm

### 3.1. 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 Figure1. D is the dimension of search space, N is the number of particles, and Xi is the position of particle i, which can be expressed as follows:(11)Xi=xi1,xi2,…,xid,…,xiD,i=1,2,…,N.Figure 1
Law of universal gravitation.Through many experiments, it is found that when solving the gravity in GSA, the squareR2 of the distance between two objects in the law of universal gravitation is changed to R, which is faster and more effective in calculation. Therefore, the gravitational force Fijdt of particle j acting on particle i in the dth dimension at a certain time t can be calculated by the following formula:(12)Fijdt=GtMpit⋅MaitRijt+ε,where Mai(t) and Mpi(t) are the inertial mass of particle i and j, respectively, ε is an infinitesimal constant, Rij(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:(13)Gt=G0×e−αt/T,where G0 is the gravitational constant at the initial time t0 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 particlei in the dth dimension is the sum of forces acting on other N − 1 particles, and the size of Fidt is calculated as follows:(14)Fid=∑j=1,j≠iNrandjFijdt,where randj is a random constant between 0 and 1.The accelerationaidt of particle i in the dth dimension at time t is calculated as follows:(15)aidt=FidtMit.The velocity and position update mechanism of particlei can be expressed as follows:(16)vidt+1=randi×vidt+aidt,xidt+1=xidt+vidt+1,The inertial mass of a particle is related to its current fitness value. The 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. The inertial mass of the particle can be calculated as follows:(17)Mai=Mpi=Mii=Mi,Mit=mit∑j=1Nmjt,mit=fitt−worsttbestt−worstt,where Mii 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 sizeN, the maximum number of iterations T, and other relevant parameters and initialize the particle velocity vidt and position xidtStep 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 Mi(t)Step 3: calculate the total forceFdi acting on the ith particle and the acceleration aidtStep 4: update the velocityvidt+1 and position xidt+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


### 3.2. 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.The Levy distribution can be expressed as follows:(18)Levys∼s−λ0,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. The random flight stride length adopts an alternating mode of frequent short-distance flight and occasional long-distance flight. The characteristics of Levy flight can increase the population diversity and expand the search range of the algorithm. Therefore, using Levy flight path in intelligent optimization algorithm can make it easier to jump out of local optimum.The common formula for calculating the random step lengths of Levy flight is as follows:(19)s=uv1/β,where u and v obey normal distribution, u∼N0,σu2σu=Γ1+βsinπβ/2/Γ1+β/2β2β−1/21/β, β is usually set to 1.5, v∼N0,σu2, and σv=1.


### 3.3. 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.The updating strategy of Levy flight adopts an alternating mode of frequent short-distance flight and occasional long-distance flight, which can increase the diversity of the population, expand the search range, and help the algorithm to jump out of local optimum. Therefore, 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. The specific improvement ideas are as follows: when updating the position and velocity of particles, the particles are sorted according to the fitness values. Then, a proportional coefficientθ is determined to divide the sorted particles into two parts of A1 and A2. For the A1 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 A2 part with poor fitness values.The position updating mechanism based on Levy flight is as follows:(20)xidt+1=xidt+β×Levy_step,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.The flowchart of the improved GSA based on Levy flight is shown in Figure2, and the specific implementation steps are as follows:Step 1: set the population sizeN, the maximum number of iterations T and θ, and other related parameters; initialize the velocity vidt and position xidt of the particle.Step 2: calculate the fitness value fiti(t), update the universal gravitation constant G(t), the worst fitness value worst(t), the optimal fitness value best(t), and the inertial mass Mi(t).Step 3: calculate the total forceFdi on the particle and its acceleration aidt.Step 4: the particles are sorted according to their fitness values and divided into the better partA1 and the worse part A2 according to the proportional coefficient θ.Step 5: the particles inA1 are updated by using equation (16) and the particles in A2 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.Figure 2
Flowchart of IGSA.


## 3.1. 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 Figure1. D is the dimension of search space, N is the number of particles, and Xi is the position of particle i, which can be expressed as follows:(11)Xi=xi1,xi2,…,xid,…,xiD,i=1,2,…,N.Figure 1
Law of universal gravitation.Through many experiments, it is found that when solving the gravity in GSA, the squareR2 of the distance between two objects in the law of universal gravitation is changed to R, which is faster and more effective in calculation. Therefore, the gravitational force Fijdt of particle j acting on particle i in the dth dimension at a certain time t can be calculated by the following formula:(12)Fijdt=GtMpit⋅MaitRijt+ε,where Mai(t) and Mpi(t) are the inertial mass of particle i and j, respectively, ε is an infinitesimal constant, Rij(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:(13)Gt=G0×e−αt/T,where G0 is the gravitational constant at the initial time t0 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 particlei in the dth dimension is the sum of forces acting on other N − 1 particles, and the size of Fidt is calculated as follows:(14)Fid=∑j=1,j≠iNrandjFijdt,where randj is a random constant between 0 and 1.The accelerationaidt of particle i in the dth dimension at time t is calculated as follows:(15)aidt=FidtMit.The velocity and position update mechanism of particlei can be expressed as follows:(16)vidt+1=randi×vidt+aidt,xidt+1=xidt+vidt+1,The inertial mass of a particle is related to its current fitness value. The 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. The inertial mass of the particle can be calculated as follows:(17)Mai=Mpi=Mii=Mi,Mit=mit∑j=1Nmjt,mit=fitt−worsttbestt−worstt,where Mii 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 sizeN, the maximum number of iterations T, and other relevant parameters and initialize the particle velocity vidt and position xidtStep 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 Mi(t)Step 3: calculate the total forceFdi acting on the ith particle and the acceleration aidtStep 4: update the velocityvidt+1 and position xidt+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


## 3.2. 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.The Levy distribution can be expressed as follows:(18)Levys∼s−λ0,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. The random flight stride length adopts an alternating mode of frequent short-distance flight and occasional long-distance flight. The characteristics of Levy flight can increase the population diversity and expand the search range of the algorithm. Therefore, using Levy flight path in intelligent optimization algorithm can make it easier to jump out of local optimum.The common formula for calculating the random step lengths of Levy flight is as follows:(19)s=uv1/β,where u and v obey normal distribution, u∼N0,σu2σu=Γ1+βsinπβ/2/Γ1+β/2β2β−1/21/β, β is usually set to 1.5, v∼N0,σu2, and σv=1.


## 3.3. 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.The updating strategy of Levy flight adopts an alternating mode of frequent short-distance flight and occasional long-distance flight, which can increase the diversity of the population, expand the search range, and help the algorithm to jump out of local optimum. Therefore, 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. The specific improvement ideas are as follows: when updating the position and velocity of particles, the particles are sorted according to the fitness values. Then, a proportional coefficientθ is determined to divide the sorted particles into two parts of A1 and A2. For the A1 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 A2 part with poor fitness values.The position updating mechanism based on Levy flight is as follows:(20)xidt+1=xidt+β×Levy_step,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.The flowchart of the improved GSA based on Levy flight is shown in Figure2, and the specific implementation steps are as follows:Step 1: set the population sizeN, the maximum number of iterations T and θ, and other related parameters; initialize the velocity vidt and position xidt of the particle.Step 2: calculate the fitness value fiti(t), update the universal gravitation constant G(t), the worst fitness value worst(t), the optimal fitness value best(t), and the inertial mass Mi(t).Step 3: calculate the total forceFdi on the particle and its acceleration aidt.Step 4: the particles are sorted according to their fitness values and divided into the better partA1 and the worse part A2 according to the proportional coefficient θ.Step 5: the particles inA1 are updated by using equation (16) and the particles in A2 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.Figure 2
Flowchart of IGSA.


## 4. The Combination of VMD and IGSA

In the decomposition process of VMD, the artificial trial method is mostly adopted to determine the number of modesK, 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. Therefore, a variational modal decomposition method based on improved gravitational search algorithm (IGSA-VMD) is proposed in this paper. The 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. The better the orthogonality between the mode components, the smaller the correlation between them and the better the decomposition effect. When the number of modesK and the penalty factor α are not appropriate, the orthogonality among the mode components will be affected. Therefore, the index of orthogonality (IO) can quantitatively analyze the accuracy of signal decomposition. The smaller the IO value is, the higher the accuracy of signal decomposition and the lower the confusion degree among the modes. The formula for calculating IO value is as follows:(21)IO=∑j=1K∑h=1h≠jK∫0Lujtuktdt∫0Lx2tdt=∑j=1K∑h=1h≠jK∑i=1Hujiuki∑i=1Hx2i,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 parametersK and α of VMD by using IGSA algorithm, [K, α] is taken as the position of particle i. The signal is decomposed with this parameter combination and the IO value of decomposed components is calculated. The optimal value of the whole optimization process is the minimum IO value and the optimal particle position is the best parameter combination of [Kbest, αbest]. The flowchart of IGSA-VMD is shown in Figure 3.Figure 3
Flowchart of IGSA-VMD.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. The mathematical expression and waveform of the signal are shown in Figure4:(22)Xt=cos2π⋅20t+13cos2π⋅50t+14cos2π⋅60t+116cos2π⋅100t.Figure 4
Waveform of the simulated signal.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. The simulated signal is decomposed and the IO values are calculated. The 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. Therefore, the optimal number of modes K can be accurately identified according to the evaluation index IO.Figure 5
IO values of differentK.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. The number of particles is set to 30 and the number of maximum iterations is set to 100. The iterative curve of optimal fitness value in the optimization process is shown in Figure6, and the optimization results are listed in Table 1.Figure 6
Optimal fitness value iteration curve of the simulated signal.Table 1
Optimization results of the simulated signal.
MethodsPosition of optimal particleOptimal IO valueGSA-VMD[4, 4936]2.1714 × 10−4IGSA-VMD[4, 4999]1.3710 × 10−4It 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 α. The 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. The simulation results verify the effectiveness of the proposed IGSA-VMD method.Figure 7
Decomposition results of the simulated signal.


## 5. The Proposed Shearer Cutting State Recognition Method

Figure8 illustrates the main point of the proposed shearer cutting state recognition method, which can be summarized as follows.Figure 8
The main point of the proposed shearer cutting state recognition method.


### 5.1. 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.


### 5.2. 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 signalx with length N, the envelope entropy Ep can be calculated as follows:(23)Ep=−∑i=1Npilgpi,pi=ai∑i=1Nai,where a is the envelope of demodulation signal obtained by Hilbert transform of signal x.The kurtosisKu of each IMF component can be calculated as follows:(24)Ku=Ex−u4σ4,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 setX = {x1, x2, …, xn}, there are n sample points and each sample point contains p indicators, namely, xi ∈ Rp, i = 1, 2, …, n; then,(25)X=x11x12⋯x1px21x22⋯x2p⋮⋮⋮⋮xn1xn2⋯xnp=x1x2⋯xp,where xi = (x1j, x2j, …, xnj)T, j = 1, 2, …, p.Standardize original data: the original data are standardized by using the following equation:(26)yij=xij−μjstdj,whereμj and stdj are the mean and standard deviation of the jth dimensional vector in X.Construct the covariance matrixS: the covariance matrix S can be calculated as follows:(27)S=1p−1∑i=1pyi,yiT.The eigenvalues (λ1, λ2, …, λp) of matrix S and corresponding eigenvectors ui=u1i,u2i,…,upi,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. The calculation formula of the contribution rate is as follows:(28)δi=λi∑i=1pλi.The cumulative contribution rate of the firstk principal components can be calculated as follows:(29)Zk=∑i=1kλi∑i=1pλi.If the cumulative contribution rate is higher than the preset threshold, thek principal components can be used to evaluate the original data.


### 5.3. 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. The RBF kernel can be expressed as follows:(30)Kx,xc=exp−x−xc22γ2,where γ 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 parameterC and the width parameter γ, are optimized by particle swarm optimization algorithm (PSO) according to the literature [32]. The parameters of LSSVM are set as C (0.1, 1000) and γ (0.001, 10), according to multiple attempts.


## 5.1. 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.


## 5.2. 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 signalx with length N, the envelope entropy Ep can be calculated as follows:(23)Ep=−∑i=1Npilgpi,pi=ai∑i=1Nai,where a is the envelope of demodulation signal obtained by Hilbert transform of signal x.The kurtosisKu of each IMF component can be calculated as follows:(24)Ku=Ex−u4σ4,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 setX = {x1, x2, …, xn}, there are n sample points and each sample point contains p indicators, namely, xi ∈ Rp, i = 1, 2, …, n; then,(25)X=x11x12⋯x1px21x22⋯x2p⋮⋮⋮⋮xn1xn2⋯xnp=x1x2⋯xp,where xi = (x1j, x2j, …, xnj)T, j = 1, 2, …, p.Standardize original data: the original data are standardized by using the following equation:(26)yij=xij−μjstdj,whereμj and stdj are the mean and standard deviation of the jth dimensional vector in X.Construct the covariance matrixS: the covariance matrix S can be calculated as follows:(27)S=1p−1∑i=1pyi,yiT.The eigenvalues (λ1, λ2, …, λp) of matrix S and corresponding eigenvectors ui=u1i,u2i,…,upi,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. The calculation formula of the contribution rate is as follows:(28)δi=λi∑i=1pλi.The cumulative contribution rate of the firstk principal components can be calculated as follows:(29)Zk=∑i=1kλi∑i=1pλi.If the cumulative contribution rate is higher than the preset threshold, thek principal components can be used to evaluate the original data.


## 5.3. 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. The RBF kernel can be expressed as follows:(30)Kx,xc=exp−x−xc22γ2,where γ 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 parameterC and the width parameter γ, are optimized by particle swarm optimization algorithm (PSO) according to the literature [32]. The parameters of LSSVM are set as C (0.1, 1000) and γ (0.001, 10), according to multiple attempts.


## 6. Experimental Validation

In this section, an experimental shearer cutting coal-rock system is established to collect cutting sound signal. The 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 Figure9.Figure 9
The experiment system of shearer cutting coal-rock.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 Figure10. The 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. Thus, five cutting states of F1, F2, F3, F4, and F5 are obtained and listed in Table 2.Figure 10
Four self-made typical coal-rock specimens.Table 2
The class labels of five cutting states.
Cutting stateCoal-rock specimenClass labelF1None1F2Coal seam with hardnessf = 22F3Coal seam with hardnessf = 33F4Coal seam with gangue4F5Rock5The 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. The 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. The time domain waveforms of sound signal with five different cutting states are shown in Figure11. 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.Figure 11
The sound signals of five cutting states.In this experiment, the collected data are divided into several nonoverlapping samples. The data length of each sample is 0.5 s and each cutting state contains 60 samples. Thus, 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. These details are given in Table3.Table 3
The detailed arrangements of the experimental datasets for classification.
Cutting stateClass labelNumber of training samplesNumber of testing samplesF113030F223030F333030F443030F553030According to the procedure of the proposed cutting state recognition method, the collected signals are firstly decomposed by using IGSA-VMD. The 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. The results are shown in Figure12.Figure 12
The decomposed IMF components of sound signal based on IGSA-VMD.Then, the envelope entropy and kurtosis are used to extract sensitive feature information of each IMF component to describe different cutting states. Thus, each sample data can be represented by a feature vector of 16 dimensions, as shown in Table4, and an feature matrix of 300 × 16 can be obtained.Table 4
The feature vector of cutting sound signal.
Sample numberFeature vector1[0.7918, 0.8348, 0.6840, 0.6662, 0.6728, 0.2799, 0.5645, 0.3014, 0.0778, 0.0571, 0.1491, 0.1436, 0.1539, 0.3018, 0.2394, 0.3293]2[0.7580, 0.8097, 0.7496, 0.6550, 0.6429, 0.2441, 0.5348, 0.2728, 0.1093, 0.0638, 0.1180, 0.1777, 0.1924, 0.3794, 0.2982, 0.4088]3[0.7569, 0.7500, 0.6837, 0.6211, 0.6323, 0.2123, 0.5254, 0.2497, 0.0766, 0.0786, 0.1650, 0.2138, 0.1962, 0.4333, 0.2939, 0.4577]4[0.7979, 0.6867, 0.7044, 0.6364, 0.6128, 0.2586, 0.5207, 0.2799, 0.0858, 0.1743, 0.1943, 0.2460, 0.2906, 0.4333, 0.3933, 0.4679]5[0.7069, 0.7678, 0.6523, 0.6789, 0.6609, 0.2893, 0.5471, 0.3269, 0.0940, 0.1060, 0.1646, 0.1755, 0.2290, 0.3997, 0.3537, 0.4320]6[0.7064, 0.7961, 0.5993, 0.6139, 0.2582, 0.3263, 0.5198, 0.3258, 0.1006, 0.0940, 0.2121, 0.1974, 0.3416, 0.3454, 0.2613, 0.3433]7[0.7707, 0.6345, 0.6857, 0.6141, 0.6036, 0.2438, 0.5117, 0.2614, 0.1335, 0.1619, 0.1396, 0.2814, 0.2743, 0.4821, 0.3710, 0.5139]……298[0.5633, 0.6539, 0.6756, 0.6330, 0.3202, 0.4556, 0.5804, 0.4202, 0.3001, 0.1814, 0.1643, 0.2107, 0.3230, 0.3072, 0.2490, 0.3292]299[0.5550, 0.7038, 0.5930, 0.6688, 0.3250, 0.4413, 0.5975, 0.4143, 0.3332, 0.1067, 0.1729, 0.1559, 0.2629, 0.2559, 0.1962, 0.2670]300[0.5804, 0.7032, 0.6115, 0.6708, 0.2965, 0.4490, 0.5986, 0.4080, 0.4209, 0.1266, 0.1668, 0.1492, 0.2701, 0.2550, 0.1930, 0.2702]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. Therefore, PCA is introduced to achieve the dimension reduction processing. The histogram of principal component contribution rate is shown in Figure13. When k is equal to 5, the cumulative contribution rate Zk exceeds 90%. Then, the new matrix composed of the first 5 principal components is selected to replace the original feature vector, as shown in Table 5.Figure 13
Histogram of the principal component contribution rate.Table 5
The new feature vectors after dimension reduction.
Sample numberNew feature vector1[0.5014, 0.4801, 0.4248, 0.5117, 0.5391]2[0.4255, 0.2261, 0.3565, 0.5393, 0.5422]3[0.2646, 0.0710, 0.1204, 0.1617, 0.1227]4[0.2550, 0.1696, 0.1234, 0.1359, 0.1282]5[0.1544, 0.1578, 0.1280, 0.1221, 0.1007]6[0.1191, 0.1568, 0.1555, 0.1401, 0.1010]7[0.0741, 0.0696, 0.0689, 0.0566, 0.0628]……298[0.0624, 0.0673, 0.0654, 0.0570, 0.0542]299[0.0269, 0.0306, 0.0317, 0.0265, 0.0243]300[0.0253, 0.0299, 0.0268, 0.0251, 0.0240]The 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 Figure14 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.Figure 14
Recognition accuracy of the four methods.Table 6
The cutting state recognition results of the four methods.
MethodsAccuracy (%)Standard deviationMax.Min.MeanThe proposed method97.5092.5095.540.01669GSA − VMD + PCA95.0090.8392.920.01684VMD + PCA94.1788.3391.250.02165EEMD + PCA95.0090.0092.750.01675Furthermore, 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. The obtained state identification accuracies are described in Table7. 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. The 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. The comparison results further illustrate the necessity of selecting useful features by using PCA. The comparison results also prove the advantage and superiority of IGSA-VMD in terms of signal decomposition and state feature extraction.Table 7
The cutting state recognition results of the four methods without using PCA.
MethodsAccuracy (%)Standard deviationMax.Min.MeanIGSA-VMD88.3378.3384.340.01824GSA-VMD86.6774.1782.450.01956VMD81.6770.8378.870.2248EEMD84.1773.3380.960.01758In 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. The 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. The foregoing comparison results indicate the effectiveness and superiority of proposed method in identifying the shearer cutting state in terms of recognition accuracy.Table 8
The recognition results based on the proposed method and MRVM-based method.
MethodsAccuracy (%)Standard deviationMax.Min.MeanThe proposed method97.5092.5095.540.01669MRVM-based method [33]93.3388.3390.490.01735


## 7. 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. The 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. The effectiveness of IGSA-VMD is validated by using the simulated signals. Then, 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.



---

*Source: 8835462-2020-12-16.xml*