# Research on Fault Evolution Feature Extraction and Identification of Sun Gear Cracks in Planetary Gearbox Based on Volterra High-Order Kernel Generalized Frequency Response Graphic Analysis

Authors: Haitao Wang; Zhenya Kang; Lichen Shi; Kun Wang; Kun Li
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In this paper, relying on the Volterra series nonlinear system model and the high-order kernel Hilbert’s reconstructed kernel fast solved algorithm, a fault feature frequency domain identification method based on Volterra high-order kernel generalized frequency response graph analysis is proposed. Firstly, the method uses the system input and output vibration signals to determine the Volterra model. Then, the Volterra high-order kernel function is solved quickly by reproducing kernel Hilbert space method, and the generalized frequency response function is used to identify the model. Finally, multidimensional high-order spectral pattern analysis is used to separate and extract the fault and degree characteristic information implied by frequency and phase coupling in the third-order kernel function. Following the theoretical approach, in the experimental part, this paper uses the planetary gearbox fault loading test rig to complete the data collection and establishes the Volterra experimental model through the measured data. The generalized frequency responses of each order kernel function are compared and analyzed and the capability of distinguishing and the adaptability of different order kernel functions for the degree of crack failure are discussed. The effects of changing the memory length of the Volterra model and the order of the kernel function on the recognition result are verified. The final experimental results show that the use of reproducing kernel Hilbert space can effectively avoid the dimension disaster problem that occurs in the high-order kernel solution process. Moreover, the third-order kernel can describe more intuitively the nonlinear system model under multifactor coupling than the second-order kernel. Finally, Volterra series model the third-order kernel’s generalized frequency response can effectively distinguish between nondefective and faulty gears, and its resolution is enough to distinguish the degree of failure of gear cracks.


---

## Body

## 1. Introduction

Among the advantages of planetary gear transmission system, high transmission ratio, strong load capacity, small size, and light weight contributed to its success in many industrial applications such as wind power equipment, engineering machinery, and aerospace [1, 2]. However, planetary gear transmission systems are subjected to hostile long-lasting working conditions, e.g., variable load, frequent start and stop, and high speed, which inevitably leads to different types of failures of the critical parts of planetary gearboxes [3].Statistical analysis highlighted that the gear crack plays an important role among the most important planetary gearbox failure factors [4]. Moreover, the gear crack easily triggers other types of failures, eventually leading to a chain reaction from small to major faults. Therefore, timely diagnosis and prediction of early failures of gear cracks, advanced control and resolution of crack propagation, and timely planned maintenance will effectively avoid or reduce the probability of a complete machine accident caused by a planetary gearbox crack failure. Therefore, the study and establishment of a set of feature extraction and prediction techniques for early and weak failures of planetary gearbox gears are significant [5].In general, when the internal gear of the planetary gearbox fails, most of the vibration signals of the box body show obvious nonstationary [6–8], nonlinear, and complex features [9]. Therefore, the key to this paper is to choose a nonlinear analysis method that is suitable for the system characteristics to pick up useful information from the vibration signal and complete the crack fault feature identification. At present, among the commonly used analysis technique for fault diagnosis of planetary gearboxes there are new cepstrum method [10], spectral kurtosis (SK) method [11], wavelet transform (WT) [12, 13], singular value decomposition (SVD) [14], empirical mode decomposition (EMD) [15], and adaptive stochastic resonance (ASR) [16, 17]. Some important results have been achieved through research on fault feature identification using the above methods. Barzcz et al. proved the probability of applying spectral kurtosis to detect tooth-like cracks on planetary gear rings of wind turbines [11]. Similarly, Sayena et al. used intricate Morlet wavelets to extract features that can distinguish between healthy carriers and faults in planetary gearboxes [18]. Jiang et al. propose a signal denoising technology based on adaptive Morlet wavelet and singular value decomposition (SVD) and applied it to extract pulse characteristics of second stage reducer planetary gearbox for wind turbine [19]. Based on a set of EMD and an energy segmentation algorithm, Feng et al. proposed a combined amplitude and frequency demodulation technology to diagnose sun gear failures using a planetary gearbox test bench [20]. Further to this, Lei et al. discussed an adaptive stochastic resonance technology to diagnose sun gear faults [15]. From the direction of the above research results, many studies have focused on the qualitative analysis of typical fault types of planetary gearboxes. Of course, it does not involve the light or severity of the fault, nor does it provide a more in-depth study of the mechanism and crack evolution of the fault crack.More recently, many researches focused on investigating the degrees of failure. Ren et al. proposed a method combining the three-dimensional waterfall spectrum and a variable scale wavelet to analyze time-frequency domain characteristics of three different crack depths in a rotor [21]. Moreover, to calculate the time-varying mesh stiffness, Jiateng et al. proposed an analytical-finite element method called the assist-stress intensity factor gear contact model [22]. Haitao Wang et al. attempted to establish a nonlinear system model for rolling bearings using the Volterra series theory and used a multipulse excitation algorithm to quickly estimate the Volterra second-order kernel and their generalized frequency responses [23]. Finally, the high-resolution feature extraction and fault degree identification of the inner ring, outer ring, and rolling element pitting failure of the rolling bearing are finally achieved by analyzing the numerical characteristics of the broad-spectrum frequency response and bispectrum diagonal slice. Although the above studies have achieved specific results in the identification of faults, they still focus on the failure of a single component and do not fully meet the actual engineering needs, that asks for the extraction and prediction of the vibration signal characteristics of the cracks and the degree of evolution of planetary gearbox gears in various mechanical equipment under actual conditions.In response to the above-unsolved problems and considering our previous research results, this paper proposes a novel identification method for crack fault evolutionary characteristics of planetary gearboxes based on the Volterra high-order kernel function theory [24, 25]. The main academic ideas of this method can be summarized as follows: First, a nonlinear Volterra system model of the system is established by providing vibration signals of planetary gearboxes as the input and output; second, through the rational construction of the Volterra series Hilbert reproducing kernels [26–28], the high-order kernels of the Volterra model are calculated. Avoiding the “dimensional disaster” problem that may occur [29], third, through the generalized frequency response characteristics of high-order kernel of the system Volterra model, the crack faults and degree characteristics of planetary gearbox gears are identified efficiently. Finally, several fault evolution feature extraction experiments based on planetary gearbox fault loading experimental platform are simulated and carried out through numerical examples. Further, the correctness, universality, and superiority of the proposed method are verified proving that it further enriches and promotes fundamental research in the field of mechanical and electrical equipment fault diagnosis.


## 2. Theoretical Background

### 2.1. Volterra Series Theory

The Volterra series theory was first proposed from the perspective of mathematics and used to study some integro-differential equations. However, starting from the basic theory of signals and systems, according to the Volterra series, the output of a system can be expressed in terms of an infinite-length algebraic sum of multidimensional convolution integrals of its input values, thereby defining the system model in a specific and unique way.For a linear system, the system input and output relations can be expressed in the time domain as the following convolution operations:(1)yt=∫-∞∞ht-τμτdτwhereμ(t), y(t) are system input and output, respectively, and h(t) is impulse response transfer function (IRF).In the same way, for any continuous nonlinear system, if the inputμ(t) energy of the system is limited under the zero initial condition, the Volterra series represents the response of the system.(2)yt=h0t+∑k=1∞∫0∞⋯∫0∞hkτ1,τ2,…,τk∏i=1kμt-τidτiwherehk(τ1,τ2,…τk) is kernel of the k-order Volterra series and it is the generalization of impulse response transfer function (GIRF) of the linear system in the high-dimensional space. By comparing the linear system with the nonlinear system representation, with regard to linear system, h(t) represents the system property of the linear system, that is, the system transfer function; similarly, for a nonlinear system, the set of h1(τ), h2(τ1,τ2), hn(τ1,…,τn) represents the system properties of nonlinear system. Therefore, one can infer that the system properties of the nonlinear system can be determined and represented by the kernel function of the Volterra nonlinear model of the system.If formula (refer to formula (2)) is discretized in the time domain, a discrete time domain Volterra expression can be obtained:(3)yu=h0+∑n=1∞∑l1=1∞⋯∑ln=1∞hnl1,…,lnμu-l1⋯μnu-lnwhere μu, yu∈R is the system input and output, and hn(l1,l2,…ln) is the nth order kernel function of the discrete Volterra series.To study the kernel function of Volterra series in more depth, the analysis in the time domain is not enough. Therefore, the concept of generalized frequency response function (GFRF) is introduced to realize the research of kernel function in the frequency domain. It is defined as the multidimensional Fourier transform of the kernel function of Volterra series, as follows:(4)Hnω1,ω2,…,ωn=∫-∞∞⋯∫-∞∞hnτ1,τ2,…,τne-jω1τ1+ω2τ2+⋯+ωnτndτ1dτ2⋯dτnTherefore, the higher-order kernel functions of Volterra nonlinear system model represent the basic characteristics of the nonlinear system. At this stage of the analysis, how to quickly and efficiently determine the order Volterra kernel functions of the system to be identified becomes another important and difficult point. Searching in the literature for the solution of the Volterra kernel function shows that most of the calculation methods mainly adopt traditional approximation algorithms such as the least-squares method, multitone excitation, and neural network. This will inevitably expose the deficiencies of this type of method, such as large computational complexity, high complexity, and high-order kernels that are prone to dimension disasters. Therefore, for solving the high-order Volterra kernel function of a nonlinear coupled system under complex conditions, such as the complex planetary gearboxes, this paper proposes a novel Volterra series kernel method based on reproducing kernel Hilbert space (RKHS). This calculation method can more easily determine the Volterra kernel function of nonlinear systems, especially high-order kernel.


### 2.2. Improved Calculation Method for Solving High-Order Kernel

#### 2.2.1. Volterra Series in Hilbert Space

The discrete form of the bounded Volterra series is defined as(5)yu=h0+∑n=1∞∑l1=0m-1∑l2=0m-1∑lj=0m-1hnl1,l2,…,ln∏i=1nμt-liwherem is length of the system memory andt is time. Let the nonlinear system have an input vectorU and ask to compute the Volterra series kernel to thep-thorder.(6)U=u1,…,ui,…,um=ut,…,ut-i-1,…ut-m+1Define the spatial projection as(7)ηiu=s∣s=u1i1u2i2u3i3…umim,i=i1+i2+⋯+im,0≤i1,i2,…,im≤IIf we let η=η1,η2,…,ηi,…,ηp, different nonlinear systems can only correspond to a unique set of Volterra kernel vectors, such as h0,h1,h2,h3,…,hn. If H=h1,h2,h3,…,hn is given, the Volterra series can be represented as a compact linear form h0+η,H.Let all the nonlinear time-invariant dynamic systems (LTIs) that can be modeled by the Volterra series be the spaceH, then the Volterra series can be expressed as(8)yu=h0+η,H=h0+∑i=1qhiηiuwhere hi∈R, η is spatial projection vector input by the system, q→∞, i.e., y(u) is an infinite series.(9)y1u,y2u=h0+∑i=1qhiηi,h0+∑i=1qhi′ηi=h02+∑i=1qhihi′Since a linear space H can be called a Hilbert space after establishing an inner product, the space H at this time is Hilbert space.If the input function in Hilbert space isk(u,v), then there is(10)ku,v=1+∑i=1qηiuηivThe input function k(u,v) is constructed from the spatial projection η(u) of the original system input value u. It can be shown that function k(u,v) also has the following two important properties.( 1 ) Fixv, then k(·,v) is an element in the Hilbert spaceH, where (·) represents any one of input vector.Proof:(11)k·,v=1+∑i=1qηi·ηi·=1+∑i=1qηivηi·Let (ηi(v))=hi, h0=1, then (11) can be transformed to (12)k·,v=h0+∑i=1qhiηi·The control can be obtained k·,v∈H.( 2 ) For any element in the Hilbert spaceH, y·=h0+∑i=1qhiηi·, k·,v and its inner product are values y(v) where y is at v.Proof:(13)k·,v,y·=1+∑i=1qηivηi·,h0+∑i=1qhiηi·=h0+∑i=1qhiηiv=yvAccording to the kernel Hilbert space theory, the linear functionalk(·,v) of the Hilbert space with both properties mentioned above is called a “reproducing kernel Hilbert space” (RKHS). According to the RKHS theory, the nonlinear system can be described as an element in Hilbert space; assuming some input values and corresponding output values of the nonlinear system were known, then the system can be reconstructed. It is a coincidence that this is the process of modeling nonlinear systems using the Volterra series theory, proving that there is interchangeability between the two. Moreover, this important theoretical discovery lays on solid mathematical basis for solving the kernel function of the nonlinear system Volterra model.Assume that the system input value is(u1,u2,…,uN) and the observation value after passing through the nonlinear system is y=(y1,y2,…,yN). Then, according to the second property of the RKHS, it is known that yv=k·,v,y·. The solution for y(·) in (13) requires the value of y(t) is known for any input, as in(14)
k
·
,
u
1
,
y
·
=
y
u
1
k
·
,
u
2
,
y
·
=
y
u
2
⋮
⋮
⋮
k
·
,
u
N
,
y
·
=
y
u
NAccording to linear space theory, a nonlinear model to be identified corresponds to the output valuey(·) of the input (u1,u2,…,uN), which is the projection of the response of any input signal of the system in k(·,u1), k(·,u2),…,k(·,uN), the linear subspace into which the system responds. According to the Hilbert space theory of reproducing kernels, the better approximation of the Volterra series representation of a nonlinear model can be represented by a linear combination of k(·,u1), k(·,u2),…,k(·,uN), i.e.,(15)y·=∑j=1Nqjk·,uj


#### 2.2.2. Explicit Expressions of Each Intermediate Parameter and Volterra Series Kernel

This paper replaces the (·) in (15) with U∈u1,u2,…,ui,…,uN for a more clear and visual representation of the Volterra series kernel h.(16)yui=∑j=1Nqjkui,uj,i∈1,…,Nif(17)K=k11…k1N⋮⋱⋮kN1⋯kNNAmong them there is kij=kui,uj. By rewriting (16) into a matrix, we get(18)y=KQFrom the solution(19)Q=K-1ywhere K-1 represents the generalized inverse of the matrix K. According to (16), (20)yu=∑j=1Nqjku,uj=∑j=1Nqj1+∑i=1qηiuηiuj=∑j=1Nqj+∑i=1q∑j=1NqjηiujηiuAfter comparing the coefficients of (20) and (16), we can get the Volterra series kernel as(21)h0=∑j=1Nqj,hi=∑j=1NqjηiujThis solution highlights that improved method for solving the kernel of the Volterra series model based on RKHS avoids the exponentially explosive dimensionality of the iterations of conventional calculation methods and also solves Volterra high-order kernel problems, providing a new way of thinking.


#### 2.2.3. Examples and Simulations

(22) y t = 1 + sin ⁡ 2.1 u t + 1.3 u t - 1 + 0.6 u t - 2 e - 3 2.1 u t + 1.3 u t - 1 + 0.6 u t - 2 + 2 3Let the excitation signal input by the system be a random signal uniformly distributed over the range(-0.15,0.15). Observe 500 sets of output data. Then use the Volterra series expansion formula identified by the proposed method as(23)y^t=1+0.2339ut+0.1439ut-1+0.0662ut-2+0.4353ut2+0.5398utut-1+0.2492utut-2+0.1671ut-12+0.1539ut-1ut-2+0.0355ut-22-0.8931u3+1.6595u2ut-1+0.7638u2ut-2+1.0275uut-12+0.9479uut-1ut-2+0.2181uut-22+0.2124ut-13+0.2938ut-12ut-2+0.1342ut-1ut-22+0.0219ut-23-0.7509uut-1ut-22+1.2365u4+3.0615u3ut-1+2.6321u2ut-1ut-2+1.6258uut-12ut-2+1.4142u3ut-2+2.8427u2ut-12+0.6075u2ut-22+1.1746uut-13+0.1159uut-23+0.3359ut-13ut-2+0.2319ut-12ut-22+0.0718ut-1ut-23+0.1814ut-14+0.0081ut-24+O·5(  1) Comparison of Time Domain Output.For nonlinear systems, this paper applies u(t) excitation signals(24)ut=0124131133tsin⁡tThe number of sampling points was 120, and the time domain response of the system output based on the simulation is shown in Figure1; except for some points in the middle and the tail of the signal, the system well approximates the output at most sampling points.Figure 1
System time domain response comparison.(  2) Comparison of Generalized Frequency Response of First-Order Kernel Functions. Because the effect of the DC component on the amplitude-frequency response is the same as that of the phase-frequency response, only the first-order kernel exists in the Volterra series of the nonlinear system when the zero-order or high-order Volterra kernels are ignored. Therefore, in the case of neglecting these Volterra kernels, only the first-order kernel exists in the Volterra series of nonlinear systems. The comparison results of the frequency response of the first-order kernel functions of an ideal system with those of the identified model are shown in Figure 2 and Figure 3, respectively. Figures show that the amplitude error at any point never exceeds 1.2 dB, while the phase angle error never exceeds 4 degrees.Figure 2
Comparison of amplitude-frequency response of the first-order kernel function of the system.Figure 3
Comparison of phase-frequency response of the first-order kernel function of the system.(  3) Comparison of Generalized Frequency Response of Second-Order Kernel Function.The comparison between the first-order kernel function amplitude-frequency response and the first-order kernel function phase-frequency response of the above system can be intuitively obtained. The system identification result has a good identification accuracy from the perspective of frequency response analysis. However, because the identified system is a nonlinear system, the frequency response of the first-order kernel function can only reflect the linear characteristics of the system. Therefore, the similarity of the frequency response of the first-order kernel function does not mean that the identified model can completely represent the original system. For this reason, the input signal u(t) and the output signal y(t) in this paper are used to solve the second-order kernel function amplitude-frequency response of the system (Figure 4(a)) and the second-order kernel function phase-frequency response (Figure 4(b)).Figure 4
Second-order kernel responses of the real system.
(a)
3D diagram of amplitude-frequency response curve of Volterra second-order (b)
3D diagram of phase-frequency response curve of Volterra second-orderThe frequency response of the second-order kernel function of the system obtained by reproducing kernel Hilbert space method is shown in Figure5(a) and its phase-frequency response in Figure 5(b). The amplitude frequency and phase-frequency errors between the identification model and the real system model are shown in Figures 6(a) and 6(b).Figure 5
Second-order kernel responses of the identified model.
(a)
3D diagram of amplitude-frequency response curve of Volterra second-order (b)
3D diagram of phase-frequency response curve of Volterra second-orderFigure 6
Error analysis of real and identification models for nonlinear systems.
(a)
Magnitude response error (b)
Phase-frequency response errorTo sum up, the correctness, generalization, and universality of the Volterra kernel function algorithm based on RKHS are proved by comparison and analysis of the first-order kernel, second-order kernel amplitude frequency, and phase-frequency response of the nonlinear simulation example model. By analogy, according to the calculation method described in this paper, the third-order kernel, fourth-order kernel, and even higher orders of the nonlinear system Volterra model can be solved, thus providing a reliable and complete mathematical theoretical basis for the subsequent acquisition, extraction, and tracking algorithm for the evolution characteristics of the early weak crack faults in the planetary gearbox gear.Section3 presents two basic studies on the mechanism of crack failure of the planetary gearbox and the theoretical model of planetary gear system and also discusses the influence of the crack evolution on the frequency response of the dynamic model.


## 2.1. Volterra Series Theory

The Volterra series theory was first proposed from the perspective of mathematics and used to study some integro-differential equations. However, starting from the basic theory of signals and systems, according to the Volterra series, the output of a system can be expressed in terms of an infinite-length algebraic sum of multidimensional convolution integrals of its input values, thereby defining the system model in a specific and unique way.For a linear system, the system input and output relations can be expressed in the time domain as the following convolution operations:(1)yt=∫-∞∞ht-τμτdτwhereμ(t), y(t) are system input and output, respectively, and h(t) is impulse response transfer function (IRF).In the same way, for any continuous nonlinear system, if the inputμ(t) energy of the system is limited under the zero initial condition, the Volterra series represents the response of the system.(2)yt=h0t+∑k=1∞∫0∞⋯∫0∞hkτ1,τ2,…,τk∏i=1kμt-τidτiwherehk(τ1,τ2,…τk) is kernel of the k-order Volterra series and it is the generalization of impulse response transfer function (GIRF) of the linear system in the high-dimensional space. By comparing the linear system with the nonlinear system representation, with regard to linear system, h(t) represents the system property of the linear system, that is, the system transfer function; similarly, for a nonlinear system, the set of h1(τ), h2(τ1,τ2), hn(τ1,…,τn) represents the system properties of nonlinear system. Therefore, one can infer that the system properties of the nonlinear system can be determined and represented by the kernel function of the Volterra nonlinear model of the system.If formula (refer to formula (2)) is discretized in the time domain, a discrete time domain Volterra expression can be obtained:(3)yu=h0+∑n=1∞∑l1=1∞⋯∑ln=1∞hnl1,…,lnμu-l1⋯μnu-lnwhere μu, yu∈R is the system input and output, and hn(l1,l2,…ln) is the nth order kernel function of the discrete Volterra series.To study the kernel function of Volterra series in more depth, the analysis in the time domain is not enough. Therefore, the concept of generalized frequency response function (GFRF) is introduced to realize the research of kernel function in the frequency domain. It is defined as the multidimensional Fourier transform of the kernel function of Volterra series, as follows:(4)Hnω1,ω2,…,ωn=∫-∞∞⋯∫-∞∞hnτ1,τ2,…,τne-jω1τ1+ω2τ2+⋯+ωnτndτ1dτ2⋯dτnTherefore, the higher-order kernel functions of Volterra nonlinear system model represent the basic characteristics of the nonlinear system. At this stage of the analysis, how to quickly and efficiently determine the order Volterra kernel functions of the system to be identified becomes another important and difficult point. Searching in the literature for the solution of the Volterra kernel function shows that most of the calculation methods mainly adopt traditional approximation algorithms such as the least-squares method, multitone excitation, and neural network. This will inevitably expose the deficiencies of this type of method, such as large computational complexity, high complexity, and high-order kernels that are prone to dimension disasters. Therefore, for solving the high-order Volterra kernel function of a nonlinear coupled system under complex conditions, such as the complex planetary gearboxes, this paper proposes a novel Volterra series kernel method based on reproducing kernel Hilbert space (RKHS). This calculation method can more easily determine the Volterra kernel function of nonlinear systems, especially high-order kernel.


## 2.2. Improved Calculation Method for Solving High-Order Kernel

### 2.2.1. Volterra Series in Hilbert Space

The discrete form of the bounded Volterra series is defined as(5)yu=h0+∑n=1∞∑l1=0m-1∑l2=0m-1∑lj=0m-1hnl1,l2,…,ln∏i=1nμt-liwherem is length of the system memory andt is time. Let the nonlinear system have an input vectorU and ask to compute the Volterra series kernel to thep-thorder.(6)U=u1,…,ui,…,um=ut,…,ut-i-1,…ut-m+1Define the spatial projection as(7)ηiu=s∣s=u1i1u2i2u3i3…umim,i=i1+i2+⋯+im,0≤i1,i2,…,im≤IIf we let η=η1,η2,…,ηi,…,ηp, different nonlinear systems can only correspond to a unique set of Volterra kernel vectors, such as h0,h1,h2,h3,…,hn. If H=h1,h2,h3,…,hn is given, the Volterra series can be represented as a compact linear form h0+η,H.Let all the nonlinear time-invariant dynamic systems (LTIs) that can be modeled by the Volterra series be the spaceH, then the Volterra series can be expressed as(8)yu=h0+η,H=h0+∑i=1qhiηiuwhere hi∈R, η is spatial projection vector input by the system, q→∞, i.e., y(u) is an infinite series.(9)y1u,y2u=h0+∑i=1qhiηi,h0+∑i=1qhi′ηi=h02+∑i=1qhihi′Since a linear space H can be called a Hilbert space after establishing an inner product, the space H at this time is Hilbert space.If the input function in Hilbert space isk(u,v), then there is(10)ku,v=1+∑i=1qηiuηivThe input function k(u,v) is constructed from the spatial projection η(u) of the original system input value u. It can be shown that function k(u,v) also has the following two important properties.( 1 ) Fixv, then k(·,v) is an element in the Hilbert spaceH, where (·) represents any one of input vector.Proof:(11)k·,v=1+∑i=1qηi·ηi·=1+∑i=1qηivηi·Let (ηi(v))=hi, h0=1, then (11) can be transformed to (12)k·,v=h0+∑i=1qhiηi·The control can be obtained k·,v∈H.( 2 ) For any element in the Hilbert spaceH, y·=h0+∑i=1qhiηi·, k·,v and its inner product are values y(v) where y is at v.Proof:(13)k·,v,y·=1+∑i=1qηivηi·,h0+∑i=1qhiηi·=h0+∑i=1qhiηiv=yvAccording to the kernel Hilbert space theory, the linear functionalk(·,v) of the Hilbert space with both properties mentioned above is called a “reproducing kernel Hilbert space” (RKHS). According to the RKHS theory, the nonlinear system can be described as an element in Hilbert space; assuming some input values and corresponding output values of the nonlinear system were known, then the system can be reconstructed. It is a coincidence that this is the process of modeling nonlinear systems using the Volterra series theory, proving that there is interchangeability between the two. Moreover, this important theoretical discovery lays on solid mathematical basis for solving the kernel function of the nonlinear system Volterra model.Assume that the system input value is(u1,u2,…,uN) and the observation value after passing through the nonlinear system is y=(y1,y2,…,yN). Then, according to the second property of the RKHS, it is known that yv=k·,v,y·. The solution for y(·) in (13) requires the value of y(t) is known for any input, as in(14)
k
·
,
u
1
,
y
·
=
y
u
1
k
·
,
u
2
,
y
·
=
y
u
2
⋮
⋮
⋮
k
·
,
u
N
,
y
·
=
y
u
NAccording to linear space theory, a nonlinear model to be identified corresponds to the output valuey(·) of the input (u1,u2,…,uN), which is the projection of the response of any input signal of the system in k(·,u1), k(·,u2),…,k(·,uN), the linear subspace into which the system responds. According to the Hilbert space theory of reproducing kernels, the better approximation of the Volterra series representation of a nonlinear model can be represented by a linear combination of k(·,u1), k(·,u2),…,k(·,uN), i.e.,(15)y·=∑j=1Nqjk·,uj


### 2.2.2. Explicit Expressions of Each Intermediate Parameter and Volterra Series Kernel

This paper replaces the (·) in (15) with U∈u1,u2,…,ui,…,uN for a more clear and visual representation of the Volterra series kernel h.(16)yui=∑j=1Nqjkui,uj,i∈1,…,Nif(17)K=k11…k1N⋮⋱⋮kN1⋯kNNAmong them there is kij=kui,uj. By rewriting (16) into a matrix, we get(18)y=KQFrom the solution(19)Q=K-1ywhere K-1 represents the generalized inverse of the matrix K. According to (16), (20)yu=∑j=1Nqjku,uj=∑j=1Nqj1+∑i=1qηiuηiuj=∑j=1Nqj+∑i=1q∑j=1NqjηiujηiuAfter comparing the coefficients of (20) and (16), we can get the Volterra series kernel as(21)h0=∑j=1Nqj,hi=∑j=1NqjηiujThis solution highlights that improved method for solving the kernel of the Volterra series model based on RKHS avoids the exponentially explosive dimensionality of the iterations of conventional calculation methods and also solves Volterra high-order kernel problems, providing a new way of thinking.


### 2.2.3. Examples and Simulations

(22) y t = 1 + sin ⁡ 2.1 u t + 1.3 u t - 1 + 0.6 u t - 2 e - 3 2.1 u t + 1.3 u t - 1 + 0.6 u t - 2 + 2 3Let the excitation signal input by the system be a random signal uniformly distributed over the range(-0.15,0.15). Observe 500 sets of output data. Then use the Volterra series expansion formula identified by the proposed method as(23)y^t=1+0.2339ut+0.1439ut-1+0.0662ut-2+0.4353ut2+0.5398utut-1+0.2492utut-2+0.1671ut-12+0.1539ut-1ut-2+0.0355ut-22-0.8931u3+1.6595u2ut-1+0.7638u2ut-2+1.0275uut-12+0.9479uut-1ut-2+0.2181uut-22+0.2124ut-13+0.2938ut-12ut-2+0.1342ut-1ut-22+0.0219ut-23-0.7509uut-1ut-22+1.2365u4+3.0615u3ut-1+2.6321u2ut-1ut-2+1.6258uut-12ut-2+1.4142u3ut-2+2.8427u2ut-12+0.6075u2ut-22+1.1746uut-13+0.1159uut-23+0.3359ut-13ut-2+0.2319ut-12ut-22+0.0718ut-1ut-23+0.1814ut-14+0.0081ut-24+O·5(  1) Comparison of Time Domain Output.For nonlinear systems, this paper applies u(t) excitation signals(24)ut=0124131133tsin⁡tThe number of sampling points was 120, and the time domain response of the system output based on the simulation is shown in Figure1; except for some points in the middle and the tail of the signal, the system well approximates the output at most sampling points.Figure 1
System time domain response comparison.(  2) Comparison of Generalized Frequency Response of First-Order Kernel Functions. Because the effect of the DC component on the amplitude-frequency response is the same as that of the phase-frequency response, only the first-order kernel exists in the Volterra series of the nonlinear system when the zero-order or high-order Volterra kernels are ignored. Therefore, in the case of neglecting these Volterra kernels, only the first-order kernel exists in the Volterra series of nonlinear systems. The comparison results of the frequency response of the first-order kernel functions of an ideal system with those of the identified model are shown in Figure 2 and Figure 3, respectively. Figures show that the amplitude error at any point never exceeds 1.2 dB, while the phase angle error never exceeds 4 degrees.Figure 2
Comparison of amplitude-frequency response of the first-order kernel function of the system.Figure 3
Comparison of phase-frequency response of the first-order kernel function of the system.(  3) Comparison of Generalized Frequency Response of Second-Order Kernel Function.The comparison between the first-order kernel function amplitude-frequency response and the first-order kernel function phase-frequency response of the above system can be intuitively obtained. The system identification result has a good identification accuracy from the perspective of frequency response analysis. However, because the identified system is a nonlinear system, the frequency response of the first-order kernel function can only reflect the linear characteristics of the system. Therefore, the similarity of the frequency response of the first-order kernel function does not mean that the identified model can completely represent the original system. For this reason, the input signal u(t) and the output signal y(t) in this paper are used to solve the second-order kernel function amplitude-frequency response of the system (Figure 4(a)) and the second-order kernel function phase-frequency response (Figure 4(b)).Figure 4
Second-order kernel responses of the real system.
(a)
3D diagram of amplitude-frequency response curve of Volterra second-order (b)
3D diagram of phase-frequency response curve of Volterra second-orderThe frequency response of the second-order kernel function of the system obtained by reproducing kernel Hilbert space method is shown in Figure5(a) and its phase-frequency response in Figure 5(b). The amplitude frequency and phase-frequency errors between the identification model and the real system model are shown in Figures 6(a) and 6(b).Figure 5
Second-order kernel responses of the identified model.
(a)
3D diagram of amplitude-frequency response curve of Volterra second-order (b)
3D diagram of phase-frequency response curve of Volterra second-orderFigure 6
Error analysis of real and identification models for nonlinear systems.
(a)
Magnitude response error (b)
Phase-frequency response errorTo sum up, the correctness, generalization, and universality of the Volterra kernel function algorithm based on RKHS are proved by comparison and analysis of the first-order kernel, second-order kernel amplitude frequency, and phase-frequency response of the nonlinear simulation example model. By analogy, according to the calculation method described in this paper, the third-order kernel, fourth-order kernel, and even higher orders of the nonlinear system Volterra model can be solved, thus providing a reliable and complete mathematical theoretical basis for the subsequent acquisition, extraction, and tracking algorithm for the evolution characteristics of the early weak crack faults in the planetary gearbox gear.Section3 presents two basic studies on the mechanism of crack failure of the planetary gearbox and the theoretical model of planetary gear system and also discusses the influence of the crack evolution on the frequency response of the dynamic model.


## 2.2.1. Volterra Series in Hilbert Space

The discrete form of the bounded Volterra series is defined as(5)yu=h0+∑n=1∞∑l1=0m-1∑l2=0m-1∑lj=0m-1hnl1,l2,…,ln∏i=1nμt-liwherem is length of the system memory andt is time. Let the nonlinear system have an input vectorU and ask to compute the Volterra series kernel to thep-thorder.(6)U=u1,…,ui,…,um=ut,…,ut-i-1,…ut-m+1Define the spatial projection as(7)ηiu=s∣s=u1i1u2i2u3i3…umim,i=i1+i2+⋯+im,0≤i1,i2,…,im≤IIf we let η=η1,η2,…,ηi,…,ηp, different nonlinear systems can only correspond to a unique set of Volterra kernel vectors, such as h0,h1,h2,h3,…,hn. If H=h1,h2,h3,…,hn is given, the Volterra series can be represented as a compact linear form h0+η,H.Let all the nonlinear time-invariant dynamic systems (LTIs) that can be modeled by the Volterra series be the spaceH, then the Volterra series can be expressed as(8)yu=h0+η,H=h0+∑i=1qhiηiuwhere hi∈R, η is spatial projection vector input by the system, q→∞, i.e., y(u) is an infinite series.(9)y1u,y2u=h0+∑i=1qhiηi,h0+∑i=1qhi′ηi=h02+∑i=1qhihi′Since a linear space H can be called a Hilbert space after establishing an inner product, the space H at this time is Hilbert space.If the input function in Hilbert space isk(u,v), then there is(10)ku,v=1+∑i=1qηiuηivThe input function k(u,v) is constructed from the spatial projection η(u) of the original system input value u. It can be shown that function k(u,v) also has the following two important properties.( 1 ) Fixv, then k(·,v) is an element in the Hilbert spaceH, where (·) represents any one of input vector.Proof:(11)k·,v=1+∑i=1qηi·ηi·=1+∑i=1qηivηi·Let (ηi(v))=hi, h0=1, then (11) can be transformed to (12)k·,v=h0+∑i=1qhiηi·The control can be obtained k·,v∈H.( 2 ) For any element in the Hilbert spaceH, y·=h0+∑i=1qhiηi·, k·,v and its inner product are values y(v) where y is at v.Proof:(13)k·,v,y·=1+∑i=1qηivηi·,h0+∑i=1qhiηi·=h0+∑i=1qhiηiv=yvAccording to the kernel Hilbert space theory, the linear functionalk(·,v) of the Hilbert space with both properties mentioned above is called a “reproducing kernel Hilbert space” (RKHS). According to the RKHS theory, the nonlinear system can be described as an element in Hilbert space; assuming some input values and corresponding output values of the nonlinear system were known, then the system can be reconstructed. It is a coincidence that this is the process of modeling nonlinear systems using the Volterra series theory, proving that there is interchangeability between the two. Moreover, this important theoretical discovery lays on solid mathematical basis for solving the kernel function of the nonlinear system Volterra model.Assume that the system input value is(u1,u2,…,uN) and the observation value after passing through the nonlinear system is y=(y1,y2,…,yN). Then, according to the second property of the RKHS, it is known that yv=k·,v,y·. The solution for y(·) in (13) requires the value of y(t) is known for any input, as in(14)
k
·
,
u
1
,
y
·
=
y
u
1
k
·
,
u
2
,
y
·
=
y
u
2
⋮
⋮
⋮
k
·
,
u
N
,
y
·
=
y
u
NAccording to linear space theory, a nonlinear model to be identified corresponds to the output valuey(·) of the input (u1,u2,…,uN), which is the projection of the response of any input signal of the system in k(·,u1), k(·,u2),…,k(·,uN), the linear subspace into which the system responds. According to the Hilbert space theory of reproducing kernels, the better approximation of the Volterra series representation of a nonlinear model can be represented by a linear combination of k(·,u1), k(·,u2),…,k(·,uN), i.e.,(15)y·=∑j=1Nqjk·,uj


## 2.2.2. Explicit Expressions of Each Intermediate Parameter and Volterra Series Kernel

This paper replaces the (·) in (15) with U∈u1,u2,…,ui,…,uN for a more clear and visual representation of the Volterra series kernel h.(16)yui=∑j=1Nqjkui,uj,i∈1,…,Nif(17)K=k11…k1N⋮⋱⋮kN1⋯kNNAmong them there is kij=kui,uj. By rewriting (16) into a matrix, we get(18)y=KQFrom the solution(19)Q=K-1ywhere K-1 represents the generalized inverse of the matrix K. According to (16), (20)yu=∑j=1Nqjku,uj=∑j=1Nqj1+∑i=1qηiuηiuj=∑j=1Nqj+∑i=1q∑j=1NqjηiujηiuAfter comparing the coefficients of (20) and (16), we can get the Volterra series kernel as(21)h0=∑j=1Nqj,hi=∑j=1NqjηiujThis solution highlights that improved method for solving the kernel of the Volterra series model based on RKHS avoids the exponentially explosive dimensionality of the iterations of conventional calculation methods and also solves Volterra high-order kernel problems, providing a new way of thinking.


## 2.2.3. Examples and Simulations

(22) y t = 1 + sin ⁡ 2.1 u t + 1.3 u t - 1 + 0.6 u t - 2 e - 3 2.1 u t + 1.3 u t - 1 + 0.6 u t - 2 + 2 3Let the excitation signal input by the system be a random signal uniformly distributed over the range(-0.15,0.15). Observe 500 sets of output data. Then use the Volterra series expansion formula identified by the proposed method as(23)y^t=1+0.2339ut+0.1439ut-1+0.0662ut-2+0.4353ut2+0.5398utut-1+0.2492utut-2+0.1671ut-12+0.1539ut-1ut-2+0.0355ut-22-0.8931u3+1.6595u2ut-1+0.7638u2ut-2+1.0275uut-12+0.9479uut-1ut-2+0.2181uut-22+0.2124ut-13+0.2938ut-12ut-2+0.1342ut-1ut-22+0.0219ut-23-0.7509uut-1ut-22+1.2365u4+3.0615u3ut-1+2.6321u2ut-1ut-2+1.6258uut-12ut-2+1.4142u3ut-2+2.8427u2ut-12+0.6075u2ut-22+1.1746uut-13+0.1159uut-23+0.3359ut-13ut-2+0.2319ut-12ut-22+0.0718ut-1ut-23+0.1814ut-14+0.0081ut-24+O·5(  1) Comparison of Time Domain Output.For nonlinear systems, this paper applies u(t) excitation signals(24)ut=0124131133tsin⁡tThe number of sampling points was 120, and the time domain response of the system output based on the simulation is shown in Figure1; except for some points in the middle and the tail of the signal, the system well approximates the output at most sampling points.Figure 1
System time domain response comparison.(  2) Comparison of Generalized Frequency Response of First-Order Kernel Functions. Because the effect of the DC component on the amplitude-frequency response is the same as that of the phase-frequency response, only the first-order kernel exists in the Volterra series of the nonlinear system when the zero-order or high-order Volterra kernels are ignored. Therefore, in the case of neglecting these Volterra kernels, only the first-order kernel exists in the Volterra series of nonlinear systems. The comparison results of the frequency response of the first-order kernel functions of an ideal system with those of the identified model are shown in Figure 2 and Figure 3, respectively. Figures show that the amplitude error at any point never exceeds 1.2 dB, while the phase angle error never exceeds 4 degrees.Figure 2
Comparison of amplitude-frequency response of the first-order kernel function of the system.Figure 3
Comparison of phase-frequency response of the first-order kernel function of the system.(  3) Comparison of Generalized Frequency Response of Second-Order Kernel Function.The comparison between the first-order kernel function amplitude-frequency response and the first-order kernel function phase-frequency response of the above system can be intuitively obtained. The system identification result has a good identification accuracy from the perspective of frequency response analysis. However, because the identified system is a nonlinear system, the frequency response of the first-order kernel function can only reflect the linear characteristics of the system. Therefore, the similarity of the frequency response of the first-order kernel function does not mean that the identified model can completely represent the original system. For this reason, the input signal u(t) and the output signal y(t) in this paper are used to solve the second-order kernel function amplitude-frequency response of the system (Figure 4(a)) and the second-order kernel function phase-frequency response (Figure 4(b)).Figure 4
Second-order kernel responses of the real system.
(a)
3D diagram of amplitude-frequency response curve of Volterra second-order (b)
3D diagram of phase-frequency response curve of Volterra second-orderThe frequency response of the second-order kernel function of the system obtained by reproducing kernel Hilbert space method is shown in Figure5(a) and its phase-frequency response in Figure 5(b). The amplitude frequency and phase-frequency errors between the identification model and the real system model are shown in Figures 6(a) and 6(b).Figure 5
Second-order kernel responses of the identified model.
(a)
3D diagram of amplitude-frequency response curve of Volterra second-order (b)
3D diagram of phase-frequency response curve of Volterra second-orderFigure 6
Error analysis of real and identification models for nonlinear systems.
(a)
Magnitude response error (b)
Phase-frequency response errorTo sum up, the correctness, generalization, and universality of the Volterra kernel function algorithm based on RKHS are proved by comparison and analysis of the first-order kernel, second-order kernel amplitude frequency, and phase-frequency response of the nonlinear simulation example model. By analogy, according to the calculation method described in this paper, the third-order kernel, fourth-order kernel, and even higher orders of the nonlinear system Volterra model can be solved, thus providing a reliable and complete mathematical theoretical basis for the subsequent acquisition, extraction, and tracking algorithm for the evolution characteristics of the early weak crack faults in the planetary gearbox gear.Section3 presents two basic studies on the mechanism of crack failure of the planetary gearbox and the theoretical model of planetary gear system and also discusses the influence of the crack evolution on the frequency response of the dynamic model.


## 3. Failure Mechanism Analysis

### 3.1. Planetary Gear System Dynamic Modeling

According to the dynamic model of planetary gear transmission system shown in Figure7, the system consists of a sun gear, a planet carrier, a ring gear, and three planetary gears. θc, θs, and θpi are rotation angle of the planet carrier, the sun gear, and the planet wheels, respectively; Kspi, Krpi, Cspi, Crpi are the meshing stiffness and meshing damping coefficient of the sun gear, the ring gear, and the planet gears, respectively; 2bspi, 2brpi are the meshing gear pair side clearance; TD and TL are input torque and load torque, respectively.Figure 7
Dynamic model of planetary transmission system.From the Lagrange equation, the differential equation for the dynamic model shown in Figure7 can be obtained:(25)Jsθ¨s+∑i=13Cspiθ˙srbs-θ˙pirbpi-θ˙crccos⁡α0-e˙spit+Kspitfθsrbs-θpirbpi-θcrccos⁡α0-espit,bspirbs=TDJpiθ¨pi-Cspiθ˙srbs-θ˙pirbpi-θ˙crccos⁡α0-e˙spit-Crpiθ˙pirbpi-θ˙crccos⁡α0-e˙rpit+Kspitfθsrbs-θpirbpi-θcrccos⁡α0-espit,bspi-Krpitfθpirbpi-θcrccos⁡α0-erpit,brpirbpi=0Jc+∑i=13mpirc2θ¨c-∑i=13Cspiθ˙srbs-θ˙pirbpi-θ˙crccos⁡α0-e˙spit+Kspitfθsrbs-θpirbpi-θcrccos⁡α0-espit,bspi+Krpitfθpirbpi-θcrccos⁡α0-erpit,brpirccos⁡α0=-TLwhere Jpi is dynamic moment of inertia of the planet gear i; Jc is moment of inertia of the planet carrier; rbs is base circle radius of the sun wheel; rbpi is radius of the planet’s base circle; rc is planet carrier radius; mpi is quality of the planet’s wheels; α0 is pressure angle; espi(t), erpi(t) are integrated meshing error of each gear pair; f(x,b) represents a gap nonlinear function and can be expressed as(26)fx,b=x-bx<b,0-b≤x≤b,x+bx>b.To facilitate the solution, introduce relative displacement coordinates:(27)xspi=θsrbs-θpirbpi-θcrccos⁡α0-espitxsc=θsrbs-2θcrccos⁡α0Defining the time nominal scale ωn and the displacement nominal scale bc, the dimensional-dynamic differential equation is (28)x-¨spi+Me1Ms+1Mc∑i=13Kspiτfx-spi,b-spi+2Me1Ms+1Mc∑i=13ζspix-˙spi+2MeMpiζspix-˙spi+MeMpiKspiτfx-spi,b-spi-MeMpiKrpiτfx-sc-x-spi-e-spiτ-e-rpiτ,b-rpi-2MeMpiζrpix-˙sc-x-˙spi-e-˙spiτ-e-˙rpiτ+MeMc∑i=13Krpiτfx-sc-x-spi-e-spiτ-e-rpiτ,b-rpi+2MeMc∑i=13ζrpix-˙sc-x-˙spi-e-˙spiτ-e-˙rpiτ=TDMsrbsωn2bc+TLMcrcωn2bccos⁡α0-e-¨spiτx-¨sc+Me1Ms+2Mc∑i=13Kspiτfx-spi,b-spi+2Me1Ms+2Mc∑i=13ζspix-˙spi+2MeMc∑i=13Krpiτfx-sc-x-spi-e-spiτ-e-rpiτ,b-rpi-4MeMc∑i=13ζrpix-˙sc-x-˙spi-e-˙spiτ-e-˙rpiτ=TDMsrbsωn2bc+2TLMcrcωn2bccos⁡α0whereMc=Jc/rbc2+∑i=1Nmpi/cos⁡α02; rbc is planet carrier equivalent radius; ζspi is meshing damping ratio of planeti and sun gear; ζrpi is meshing damping ratio of the planetary geari and the ring gear; Ms and Mp are equivalent masses of the sun gear and the planet gear, respectively, on their respective base circle radii.Research on the meshing stiffness of the sun gear failure shows that when the gear is slightly broken, it will cause a change in the meshing stiffness. Therefore, the gear failure factorλ∈0,1 is introduced. In the normal state λ=0, the severe fault condition of the broken tooth λ=1, the slight fault state of the gear crack 0<λ<1, and the time-varying meshing stiffness of the meshing gear pair is uniformly expressed as(29)Kt=kav+2kanπ∑n=1∞sin⁡nωmt,n=1,3,5-λkbz+∑n=1∞kbnπsin⁡2nπzcos⁡nωft+kbnπ1-cos⁡2nπzsin⁡nωftBecause the crack fault of the gear is mainly presented on the equation of the meshing stiffness, the relevant research on the meshing stiffness of the fault of a gear broken tooth by Shi Lichen et al. is used for reference [30]. In this paper, take λ=0.4 and λ=0.6, respectively, to establish the sun gear wheel crack failure meshing stiffness equation, then its expression is as follows:(30)Ksp1t=Kmt-λkfsp1tKsp2t=Kmt-λkfsp2tKsp3t=Kmt-λkfsp3tKrpit=Kmti=1,2,3


### 3.2. Dynamic Model Response Analysis

Based on the established planetary gearbox system dynamics model described in Section3.1 and on the system parameters listed in Tables 1 and 2, by solving and simulating the dynamic equations, the model response spectrum of the sun gear broken tooth can be obtained as shown in Figures 8 and 9.Table 1
Planetary gearbox configuration parameters I.
Parameter Module Number of solar gear teeth Number of planetary gear teeth Number of annular gear teeth Input torque m/mm f i n/Hz Digital quantity 2 18 27 72 15Table 2
Planetary gearbox configuration parameters II.
Parameter Meshing frequency Failure frequency of sun gearfs/Hz Transfer frequency of planet carrierfc/Hz Pressure angle f m/Hz α 0/(°) Digital quantity 216 12 3 20Figure 8
Spectrogram of model vibration response (sun gear crack of 0.2 mm).Figure 9
Spectrogram of model vibration response (sun gear crack of 0.4 mm).When a 0.2 mm gear crack fault occurs in the sun gear (Figure8), due to the amplitude modulation effect of the planetary gear revolution, the frequency spectrum in the spectrum diagram appears as the side band of the planet carrier transponder f-c; i.e., the amplitude at the kf-m±nf-c site is more obvious (being k, n an integer). For example, consider the normalized frequency f-m≈0.009041. The peak at f-m-f-c ≈ 0.008916, f-m+f-c≈0.009166 is large, and the side frequency with the sun gear failure frequency f-s as the modulation frequency also appears in the spectrum. If the normalized frequency f-m+f-s-f-c≈0.009416, f-m+f-s≈ 0.009541, f-m+f-s+2f-c≈0.009791, and f-m+2f-s-f-c≈0.009916 at the peak, side band is obvious. That is, the amplitude in the spectrogram appears at the normalized frequency kf-m±mf-s±nf-c (being k, m, n an integer). Comparing the results with those of Figure 9, when a 0.4 mm gear crack fault occurs in the sun gear, the crack fault frequency calculated by the theoretical model is basically very close to the characteristic frequency at 0.2 mm fault. The maximum normalized frequency discrimination does not exceed 0.0003, so a high-resolution and reasonable discrimination of the difference in the degree of minor failure of the crack cannot be made.At this stage, by summarizing the mechanism analysis of the cracks and the results of the frequency characteristics of the dynamic response of different degrees of failure, one can note that different degrees of early micro-crack fault eigen frequency have very high similarity and low frequency feature resolution, which is not enough to effectively distinguish faults.Therefore, in Section4, this paper carries out an experimental study of the fault feature extraction and identification method for planetary gearbox gears based on Volterra high-order kernel generalized frequency response graphical analysis. The experimental results are compared and analyzed and they prove that the proposed fault diagnosis method has certain innovation and advantages in the identification of crack fault feature resolution and degree.


## 3.1. Planetary Gear System Dynamic Modeling

According to the dynamic model of planetary gear transmission system shown in Figure7, the system consists of a sun gear, a planet carrier, a ring gear, and three planetary gears. θc, θs, and θpi are rotation angle of the planet carrier, the sun gear, and the planet wheels, respectively; Kspi, Krpi, Cspi, Crpi are the meshing stiffness and meshing damping coefficient of the sun gear, the ring gear, and the planet gears, respectively; 2bspi, 2brpi are the meshing gear pair side clearance; TD and TL are input torque and load torque, respectively.Figure 7
Dynamic model of planetary transmission system.From the Lagrange equation, the differential equation for the dynamic model shown in Figure7 can be obtained:(25)Jsθ¨s+∑i=13Cspiθ˙srbs-θ˙pirbpi-θ˙crccos⁡α0-e˙spit+Kspitfθsrbs-θpirbpi-θcrccos⁡α0-espit,bspirbs=TDJpiθ¨pi-Cspiθ˙srbs-θ˙pirbpi-θ˙crccos⁡α0-e˙spit-Crpiθ˙pirbpi-θ˙crccos⁡α0-e˙rpit+Kspitfθsrbs-θpirbpi-θcrccos⁡α0-espit,bspi-Krpitfθpirbpi-θcrccos⁡α0-erpit,brpirbpi=0Jc+∑i=13mpirc2θ¨c-∑i=13Cspiθ˙srbs-θ˙pirbpi-θ˙crccos⁡α0-e˙spit+Kspitfθsrbs-θpirbpi-θcrccos⁡α0-espit,bspi+Krpitfθpirbpi-θcrccos⁡α0-erpit,brpirccos⁡α0=-TLwhere Jpi is dynamic moment of inertia of the planet gear i; Jc is moment of inertia of the planet carrier; rbs is base circle radius of the sun wheel; rbpi is radius of the planet’s base circle; rc is planet carrier radius; mpi is quality of the planet’s wheels; α0 is pressure angle; espi(t), erpi(t) are integrated meshing error of each gear pair; f(x,b) represents a gap nonlinear function and can be expressed as(26)fx,b=x-bx<b,0-b≤x≤b,x+bx>b.To facilitate the solution, introduce relative displacement coordinates:(27)xspi=θsrbs-θpirbpi-θcrccos⁡α0-espitxsc=θsrbs-2θcrccos⁡α0Defining the time nominal scale ωn and the displacement nominal scale bc, the dimensional-dynamic differential equation is (28)x-¨spi+Me1Ms+1Mc∑i=13Kspiτfx-spi,b-spi+2Me1Ms+1Mc∑i=13ζspix-˙spi+2MeMpiζspix-˙spi+MeMpiKspiτfx-spi,b-spi-MeMpiKrpiτfx-sc-x-spi-e-spiτ-e-rpiτ,b-rpi-2MeMpiζrpix-˙sc-x-˙spi-e-˙spiτ-e-˙rpiτ+MeMc∑i=13Krpiτfx-sc-x-spi-e-spiτ-e-rpiτ,b-rpi+2MeMc∑i=13ζrpix-˙sc-x-˙spi-e-˙spiτ-e-˙rpiτ=TDMsrbsωn2bc+TLMcrcωn2bccos⁡α0-e-¨spiτx-¨sc+Me1Ms+2Mc∑i=13Kspiτfx-spi,b-spi+2Me1Ms+2Mc∑i=13ζspix-˙spi+2MeMc∑i=13Krpiτfx-sc-x-spi-e-spiτ-e-rpiτ,b-rpi-4MeMc∑i=13ζrpix-˙sc-x-˙spi-e-˙spiτ-e-˙rpiτ=TDMsrbsωn2bc+2TLMcrcωn2bccos⁡α0whereMc=Jc/rbc2+∑i=1Nmpi/cos⁡α02; rbc is planet carrier equivalent radius; ζspi is meshing damping ratio of planeti and sun gear; ζrpi is meshing damping ratio of the planetary geari and the ring gear; Ms and Mp are equivalent masses of the sun gear and the planet gear, respectively, on their respective base circle radii.Research on the meshing stiffness of the sun gear failure shows that when the gear is slightly broken, it will cause a change in the meshing stiffness. Therefore, the gear failure factorλ∈0,1 is introduced. In the normal state λ=0, the severe fault condition of the broken tooth λ=1, the slight fault state of the gear crack 0<λ<1, and the time-varying meshing stiffness of the meshing gear pair is uniformly expressed as(29)Kt=kav+2kanπ∑n=1∞sin⁡nωmt,n=1,3,5-λkbz+∑n=1∞kbnπsin⁡2nπzcos⁡nωft+kbnπ1-cos⁡2nπzsin⁡nωftBecause the crack fault of the gear is mainly presented on the equation of the meshing stiffness, the relevant research on the meshing stiffness of the fault of a gear broken tooth by Shi Lichen et al. is used for reference [30]. In this paper, take λ=0.4 and λ=0.6, respectively, to establish the sun gear wheel crack failure meshing stiffness equation, then its expression is as follows:(30)Ksp1t=Kmt-λkfsp1tKsp2t=Kmt-λkfsp2tKsp3t=Kmt-λkfsp3tKrpit=Kmti=1,2,3


## 3.2. Dynamic Model Response Analysis

Based on the established planetary gearbox system dynamics model described in Section3.1 and on the system parameters listed in Tables 1 and 2, by solving and simulating the dynamic equations, the model response spectrum of the sun gear broken tooth can be obtained as shown in Figures 8 and 9.Table 1
Planetary gearbox configuration parameters I.
Parameter Module Number of solar gear teeth Number of planetary gear teeth Number of annular gear teeth Input torque m/mm f i n/Hz Digital quantity 2 18 27 72 15Table 2
Planetary gearbox configuration parameters II.
Parameter Meshing frequency Failure frequency of sun gearfs/Hz Transfer frequency of planet carrierfc/Hz Pressure angle f m/Hz α 0/(°) Digital quantity 216 12 3 20Figure 8
Spectrogram of model vibration response (sun gear crack of 0.2 mm).Figure 9
Spectrogram of model vibration response (sun gear crack of 0.4 mm).When a 0.2 mm gear crack fault occurs in the sun gear (Figure8), due to the amplitude modulation effect of the planetary gear revolution, the frequency spectrum in the spectrum diagram appears as the side band of the planet carrier transponder f-c; i.e., the amplitude at the kf-m±nf-c site is more obvious (being k, n an integer). For example, consider the normalized frequency f-m≈0.009041. The peak at f-m-f-c ≈ 0.008916, f-m+f-c≈0.009166 is large, and the side frequency with the sun gear failure frequency f-s as the modulation frequency also appears in the spectrum. If the normalized frequency f-m+f-s-f-c≈0.009416, f-m+f-s≈ 0.009541, f-m+f-s+2f-c≈0.009791, and f-m+2f-s-f-c≈0.009916 at the peak, side band is obvious. That is, the amplitude in the spectrogram appears at the normalized frequency kf-m±mf-s±nf-c (being k, m, n an integer). Comparing the results with those of Figure 9, when a 0.4 mm gear crack fault occurs in the sun gear, the crack fault frequency calculated by the theoretical model is basically very close to the characteristic frequency at 0.2 mm fault. The maximum normalized frequency discrimination does not exceed 0.0003, so a high-resolution and reasonable discrimination of the difference in the degree of minor failure of the crack cannot be made.At this stage, by summarizing the mechanism analysis of the cracks and the results of the frequency characteristics of the dynamic response of different degrees of failure, one can note that different degrees of early micro-crack fault eigen frequency have very high similarity and low frequency feature resolution, which is not enough to effectively distinguish faults.Therefore, in Section4, this paper carries out an experimental study of the fault feature extraction and identification method for planetary gearbox gears based on Volterra high-order kernel generalized frequency response graphical analysis. The experimental results are compared and analyzed and they prove that the proposed fault diagnosis method has certain innovation and advantages in the identification of crack fault feature resolution and degree.


## 4. Experimental Study

The experimental study was carried out using the planetary gearbox fault diagnosis test bench shown in Figure10, and the detailed parameters of the planetary gearbox are listed in Table 1. In order to simulate the gear crack fault, two cracks of different depths (0.2 and 0.4 mm) were machined on one of the gears of the first stage sun gear. In the experiment, two accelerometers installed on the motor shell by using a magnetic base acquired the vibration signal. The acceleration sensors are installed at the driving end of the motor shell and the 12 o’clock position of the fan end, respectively. The load of the driving end is 0, the drive end speed is 1979 r/min, and the vibration signal enters the 32-bit data acquisition card through the signal conditioning unit to collect the high-speed data. The experimental observation data is sampled and 4000 sets of data are collected to carry on the experimental analysis and processing.Figure 10
Planetary gearbox experimental platform for fault identification.The experimental technology program continues the path developed in the previous research [23]. That is, Hilbert’s reproducing kernel space method is first used to calculate and solve the input and output data. The first-order and second-order kernel functions of the Volterra model of planetary gearbox nondefective gear and that of the two faulted gears are obtained, respectively. Moreover, use the kernel function to draw the bispectrum three-dimensional map of the two gears shown in Figure 11.Figure 11
Bispectral three-dimensional map with a memory length of 4.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthThe comparison of the results in Figure11 shows that there is a small difference in the shape of the peaks and the energy distribution among the bispectrum 3D maps of the nondefective gear and that of defective gears. Moreover, because of the too low resolution, the curves do not give enough information about fault feature information and details of its fault evolution. Therefore, to observe the features shown in Figure 11, diagonal slices are reduced in dimension obtaining a two-dimensional map (see Figure 12). The map highlights that the frequencies of fault features corresponding to nondefective gears and defective gears with the 0.2 mm and 0.4 mm cracks are all around 70 Hz, and effective distinctions among defects cannot be made. It can be seen that the spectra of second-order kernel and diagonal slice diagrams are difficult to provide a certain resolution for distinguishing crack failures. The relevant conclusions obtained in the failure mechanism analysis of the third part of this paper are verified. The above experiments can initially prove that in the process of modeling and analyzing Volterra of a more complex nonlinear system such as a planetary gearbox, because of the uncertainty of the coupling factors between the components and the transmission path in the gearbox, the Volterra model’s second-order kernel function and its corresponding generalized frequency response eventually lose the ability to identify faults.Figure 12
Bispectrum slice with memory length 4.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthThe comparison of results in Figure12 shows that although the characteristic frequency points can be captured, the correlation calculation result graph cannot extract fault features efficiently. Combining the characteristics of the Volterra series itself, there are two variable parameters for the nonlinear coupling resolution of the Volterra model built by the system to adjust its characteristics: one is the model memory length m and the other is the model order n. In this research, the change of the length of the memory, from m=4 to m=8 and m=16, was tried first and the calculation results are given in Figures 13–16.Figure 13
Bispectral 3D plot,m=8.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthFigure 14
Bispectrum slice,m=8.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthFigure 15
Bispectral 3D plot,m=16.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthFigure 16
Bispectrum slice,m=16.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthThe plots show that changing the memory length improves the gear fault identification of planetary gearboxes slightly because the change highlights the fault frequency point, which corresponds to the highest peak. However, the amplitude of other noise frequency points also increases; because of its irregular and random behavior, the equipment failure still cannot be identified. Therefore, one can conclude that changing the system memory length of the Volterra model is also not ideal for feature extraction and identification of gear cracks in planetary gearboxes.In summary, the study readjusted the experimental scheme and tried to increase the Volterra ordern from 2 to 3. However, the problems faced are very serious; that is, the computational complexity of high-order Volterra series kernels becomes more difficult with the increase of order. The dimensions of its parameters and computational data increase exponentially and fall into the vortex of dimension disasters. In order to solve the problem of dimensional disasters, this paper used the reproducing kernel Hilbert space method to calculate the Volterra series higher-order kernel. Data acquisition for input and output of the system is carried out before the experiment is extended. The third-order kernel of the Volterra nonlinear model of the planetary gearbox is obtained by commercial software MATLAB, and the generalized frequency response graph of the third-order kernel is drawn in a 4D plot, using the dimensionality reduction technique to visualize the results (see Figure 17). This technique allows to extract and compare the characteristics of gear crack faults through the generalized frequency response graph of the third-order kernel so that the failure coupling characteristics of Volterra high-order kernel functions can be expressed more clearly and intuitively.Figure 17
Higher-order analysis chart.
(a)
Nondefective gear (b)
0.2 mm crack depth (c)
0.4 mm crack depthBy observing the evolution of the Volterra third-order kernel generalized frequency response of the epicyclic gearbox as shown in Figure17, it can be seen that the high-order spectra of nondefective gear and the high-order spectra with gear cracks of 0.2 mm and 0.4 mm are very different. Comparing the plots in Figures 17(a), 17(b), and 17(c), some fault features emerge. For example, as can be seen from the whole, the shape of the nondefective gear is similar to a wavy, reversed cone, but that of the defective gear is like a column. Moreover, the colors highlight that the high-order spectral set in nondefective gear tends to be light red and yellow, whereas those of gears with crack failure are mainly blue and there are a few dark red areas.


## 5. Discussion of Experimental Results

The experimental part is divided into three stages.In the first stage, the main research results and ideas of the research group were applied. The experimental analysis of the cracks in the planetary gearbox is performed through the Volterra second-order kernel generalized frequency response characteristics. Although the algorithm can obtain the frequency response graph of the system, the characteristics of the crack failure and evolution cannot be extracted clearly from the feature analysis of diagonal slice graphics. The second-order kernel fault characteristics of the cracks and evolution of gears in planetary gearboxes do not fully comply with the frequency response characteristics of pitting failures in the inner and outer rings of rolling bearings because of the ambiguity of the frequency of the fault. Moreover, there is a broad, random interference which causes significant difficulties in fault feature extraction. These issues show that the second-order kernel generalized frequency response of the Volterra model is not suitable for fault diagnosis analysis of nonlinear coupled system like a planetary gearbox.The second stage mainly analyzes the reasons why the first phase of the experiment did not meet the expectations of the study and its solutions. Through a new analysis of the Volterra series theory, this paper found an improved method for improving the identification resolution of the Volterra model, that is, increase the memory lengthm of Volterra model and the order n of the kernel function. In subsequent experiments, the memory length m was first scaled up. The experimental results show that the increase of memory length contributes to the nonlinear coupling characteristics of Volterra second-order kernel function and its generalized frequency response, but some shortcomings such as large frequency-coupled crosstalk randomness still exist.In the third stage, mainly based on the summarization of relevant experiences in the first and second stages, a tentative study was conducted by increasing the order of kernel functionn. The experimental scheme has been improved from the original second-order kernel to the third-order kernel, but subject to the constraints of traditional algorithms, the third-order kernel calculation is more difficult. Therefore, in this experimental scheme, the third-order kernel solution of the planetary gearbox crack failure Volterra model was solved using the related results based on the improved calculation method of the Hilbert regenerative kernel studied in Section 2. By drawing the dimensionality reduction curve of higher order kernel generalized frequency response graph, the fault characteristics are more clearly defined, and it is found that the graph has higher resolution.Combining all the above experimental contents, one can conclude that the Volterra high-order kernel generalized frequency response method can effectively extract fault features of planetary gearbox cracks and it is also possible to perform high-resolution graphical analysis on the evolution of early weak faults, thus introducing a new method of fault diagnosis.


## 6. Conclusion

This paper addresses a series of nonlinear problems related to the diagnosis of gear cracks in planetary gearboxes. By introducing Volterra nonlinear system modeling theory and proposing a high-order kernel improved algorithm based on reproducing kernel Hilbert space, the problem of dimension disasters is easily avoided in the process of solving high-order kernel functions of the Volterra model. Based on the failure mechanism and experimental comparison of the measured data of the planetary gearbox experimental loading platform, the research highlights that the Volterra second-order kernel function does not have enough capacity to represent a nonlinear frequency coupling system. Through repeated attempts in the experiment, the experimental scheme is adjusted for increasing the memory length of the Volterra model and increasing the order of the kernel function. The experiments show that the method of improving the model order is more suitable for the acquisition of crack fault characteristics of planetary gearbox gears under actual conditions than increasing the length of Volterra memory. A comparative test between a nondefective gear and two defective gears with crack depths of 0.2 mm and 0.4 mm, respectively, proved the conclusions. Finally, a new method for nonlinear frequency domain coupled feature identification of planetary gearboxes based on the generalized frequency response of higher-order kernel functions of Volterra series theory is proposed. This new method has very important theoretical and practical research significance because of the very high frequency response resolution of fault signatures in nonlinear systems.



---

*Source: 6251051-2018-10-11.xml*