A load identification method based on wavelet multi-resolution analysis

Journal of Sound and Vibration 333 (2014) 381–391

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

A load identification method based on wavelet multi-resolution analysis Zong Li a, Zhipeng Feng b, Fulei Chu a,n a b

Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China Institute of Vehicular Engineering, University of Science and Technology Beijing, Beijing 100083, China

a r t i c l e i n f o

abstract

Article history: Received 6 January 2013 Received in revised form 27 May 2013 Accepted 19 September 2013 Handling Editor: L.G. Tham Available online 16 October 2013

Load identification, as a kind of indirect identification method, uses system characteristic and responses to calculate loads. A method based on wavelet multi-resolution analysis is proposed in this paper. By wavelet decomposition and transform at certain resolutions, the proposed method transforms the convolution relation between responses and loads in time domain into the linear multiplicative relation between system responses and wavelet responses in the wavelet domain. Loads can be identified as long as the linear multiplicative relation is solved. Qualitative and quantitative rules are proposed for selecting parameters that affect the accuracy of the proposed method, and are illustrated via numerical investigations. The method is illustrated by a multi-input-multi-output numerical simulation. A multi-input-multi-output laboratory experiment is performed to compare the proposed method with the frequency method on the identification ability. & 2013 Elsevier Ltd. All rights reserved.

1. Introduction With the increase in engineering requirements and the improvement of engineering technology, one kind of problem regarding the determination of the exerted loads on a dynamic system urgently needs to be solved. However, in many practical problems, dynamic loads, such as aerodynamic loads exerted on aircrafts, wind loads applied on ocean platforms, the interaction between road surface and running tires and so on, are very difficult to be directly measured owing to technical or economic limitations. Load identification, also known as load reconstruction or load deconvolution in some literatures, uses system characteristics and responses to calculate loads. Compared with direct measurement, it is an indirect identification method. It can be used to calculate loads in severe environments and loads that cannot be acquired by transducers because the laying will change system characteristic or block the motion of system components. Frequency method is one of the most commonly used identification methods. It uses a system's Frequency Response Function (FRF) and response spectrums to calculate load spectrums, and then derives loads in the time domain via the inverse Fourier transform. Some studies [1–3] determined forces in high-speed compression tests and in an experimental helicopter model with this method. However, references [4–6] found that it was very sensitive to the ill-condition of FRF and data noise. Using mathematical regularization technology, the identification ability of the method was improved [7–9]. Liu [10] analyzed two commonly used regularization methods: the truncated singular value decomposition method (TSVD) and the Tikhonov filter method in detail. Some advantages and disadvantages of the frequency method can be concluded from many literatures: owing to the linear multiplicative relation between responses and loads in the frequency domain, the

n

Corresponding author. Tel.: þ86 10 6279 2842; fax: þ 86 10 6278 8308. E-mail address: [email protected] (F. Chu).

0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.09.026

382

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

identification problem can be solved easily. However, the use of the Fourier transform requires an adequate data length. Inverse Fourier transformation must be performed at every discrete frequency and problems with numerical instability may arise when the frequency approaches the resonance frequency. Besides, ill-condition matrix and noise affect the identification accuracy severely. In contrast to the Fourier transform, which expresses a function by a linear combination of sinusoidal functions having infinite duration, the wavelet transform expresses a function as a linear combination of wavelets that have finite durations. This is one of the most important features of wavelet transform [9]. Doyle [11] transformed a time domain convolution into the wavelet domain product to solve the identification problem, which started the application of the wavelet transform in load identification. However, the functions of wavelets were not fully exploited because the scale factor was fixed; thus, the proposed method was a single-resolution one-dimension wavelet transform in essence. Mao [12] used a similar method to identify single-input impulse load. In the loading phase, the result matched the original load well, whereas in the unloading phase there was a big oscillation. Nordberg [13] parameterized road profile into coefficients by the wavelet transform method to reduce the number of unknowns. Practical data was used to validate the proposed method and the results were satisfactory. However, due to the use of iteration, the method needed initial values that were difficult to be set. Being similar to the wavelet basis functions, Gunawan [22–24] used harmonic functions, quadratic splines functions and B-splines functions as the basis functions to approximate forces, especially impact forces. Besides, the ill-posed problem and regularization were discussed in these papers. Liu [25] used basis functions that spanned only the forcing space to reconstruct distributed forces. The method required fewer basis functions. A novel load identification method based on wavelet multi-resolution analysis using the Impulse Response Function (IRF) of the system is proposed in this paper. The proposed method makes full use of the excellent features of the wavelet, such as adjustable time-frequency resolution, and several other advantages, such as that the stationary white noise transformed by orthogonal wavelet is still stationary white noise. The scale factor is adjustable in this method and this without using the iteration algorithm. The initial values do not need to be confirmed in advance. This method decomposes and transforms signals by two-dimension wavelets at a certain resolution and the convolution relation between responses and loads in time domain is transformed into the linear multiplicative relation between system responses and wavelet responses in the wavelet domain. The physical significance of the proposed method is explained during derivation. Through numerical analyses, qualitative and quantitative rules are proposed for selecting parameters that affect the accuracy of this method. The performance of the method is illustrated by a multi-input-multi-output (MIMO) numerical simulation. A MIMO laboratory experiment is also performed to compare the proposed method with the frequency method regarding identification ability. The results show that the proposed method has better identification ability and lower noise sensibility. 2. Load identification based on wavelet multi-resolution analysis 2.1. Daubechies wavelet This method adopts Daubechies (db) wavelets as the basic wavelets. db Wavelets, which have some excellent features, such as orthogonality and compact support, were created and used by the Belgian scholar Daubechies in the 1990s. db Wavelets are named dbN, where N represents the order. In all the interesting examples, the orthonormal wavelet bases can be associated with a multi-resolution analysis framework [20]. The concept of the multi-resolution analysis was introduced by Mallat in [21]. db Wavelets are orthogonal wavelets, which mean at a level (or scale) j of the wavelet basis function ϕ(t), there is Z 1 ϕð2j t  k1 Þϕð2j t  k2 Þ dt ¼ δðk1  k2 Þ; k1 ; k2 A Z: (1) 1

db Wavelets are compactly supported in [0, 2N  1]. Furthermore, for level j and shift factor k, the support interval is [2  jk, 2  j(kþ2N  1)]. 2.2. Method deduction In time domain, system responses are determined by the convolution of inputs and system characteristic yðtÞ ¼ hðtÞf ðtÞ;

(2)

where y(t) is the system response, h(t) is the unit impulse response function and f(t) is the load. Decomposing f(t) with the wavelet yields J

J

k¼I

k¼I

f ðtÞ ¼ ∑ αðkÞ2j=2 ϕð2j t  kÞ ¼ ∑ PðkÞϕð2j t  kÞ;

(3)

PðkÞ ¼ αðkÞ2j=2 ;

(4)

where α(k) are the scale coefficients, P(k) are the weight coefficients, I and J are lower and upper limits of the integral, respectively (the difference is the number of weight coefficients), j is the decomposition level, k is the shifting factor and

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

383

ϕ(2jt k) is the wavelet basis function. Strictly speaking, f(t) in Eq.(3) is an approximation of the real load. But actually, with the increase in j, f(t) will gradually approximate the real load. This will be discussed in detail in Section 2.3. Substituting Eq.(3) into Eq.(2), we have yðtÞ ¼ ∑ PðkÞφk ðtÞ;

(5)

k

φk ðtÞ ¼ hðtÞϕð2j t  kÞ:

(6)

where φk(t)is the convolution of the unit impulse response function and the wavelet basis function at level j. Similar to Eq. (2), φk(t) is defined as a wavelet response function, which is the output generated by wavelet basis function passing through the same system. φk(t) comprehensively represents the system characteristic and the wavelet function. Eq. (5) expresses the linear relation between the system response y(t) and the wavelet response φk(t). Using Eq. (3), the convolution relation in Eq. (2) is first transformed into the linear relation in Eq. (5). Applying the wavelet transform to both sides of Eq. (5) simultaneously yields Z 1 Z 1 ϕð2j t  iÞyðtÞ dt ¼ ∑ PðkÞ ϕð2j t iÞφk ðtÞ dt; i ¼ I  J: (7) 1

1

k

Let

Z uðiÞ ¼

1

ϕð2j t  iÞyðtÞ dt;

(8)

ϕð2j t  iÞφk ðtÞ dt;

(9)

1

Z g ik ¼

1 1

where u(i) is the system response in wavelet domain and gik is the wavelet response in wavelet domain. Then Eq. (7) can be simplified to J

uðiÞ ¼ ∑ PðkÞg ik ; i ¼ I  J

(10)

uðJ  I þ 1Þ1 ¼ GðJ  I þ 1ÞðJ  I þ 1Þ PðJ  I þ 1Þ1 ;

(11)

k¼I

where Eq. (11) is the matrix form of Eq. (10) and subscripts are dimensions of matrixes. Then P can be obtained as P ¼ G  1 u:

(12)

Substituting Eq. (12) into Eq. (3), the solution of f(t) can then be obtained. Using the Fourier transform, the traditional frequency method transforms the convolution relation in the time domain into the linear multiplicative relation in the frequency domain, and then solves the linear problem in the frequency domain. The proposed method first transforms the convolution relation in the time domain into the multiplicative relation in the time domain (left-hand side of Eq. (5)) and the wavelet domain (right-hand side of Eq. (5)) by wavelet decomposition (from Eqs. (2)–(6)), and then solves the linear problem in the wavelet domain by wavelet transform (from Eqs. (7)–(12)). These two methods are similar in form. However, through the analyses carried out below, the proposed method obtains better identification ability and lower noise sensibility than the frequency method. For a system with L inputs and M outputs, formulations can be derived as follows: J

uðmÞ ðiÞ ¼ ∑ P ðlÞ ðkÞg ðmlÞ ; m ¼ 1  M; l ¼ 1  L; ik

(13)

k¼I

where uðmÞ ðiÞ ¼

Z

1

ϕð2j t  iÞyðmÞ ðtÞ dt;

1

(14)

J

ðlÞ

f ðtÞ ¼ ∑ P ðlÞ ðkÞϕð2j t  kÞ;

(15)

k¼I

¼ g ðmlÞ ik

Z

1 1

ϕð2j t iÞφðmlÞ ðtÞ dt; k ðmlÞ

φðmlÞ ¼h k (m)

(16)

ϕð2j t  kÞ; (l)

is the system response in wavelet domain of the output m, P is the weight coefficient of the input l, u response in wavelet domain of the input l and the output m. Eq. (13) can be written in matrix form as uMðJ  I þ 1Þ1 ¼ GMðJ  I þ 1ÞLðJ  I þ 1Þ PLðJ  I þ 1Þ1 ;

(17) g ðmlÞ ik

is the wavelet (18)

384

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

when L is not equal to M and G is not a square matrix. The computation of Eq. (18) involves the pseudo-inverse problem. If the matrix is ill-conditioned, errors in u will actually be amplified and the overall accuracy of the identified loads will be reduced dramatically. One possible way to address this problem of error amplification is to implement a regularization method. Advantages and disadvantages of some regularization methods, like the TSVD method [14–16] and the Tikhonov method [14,16,17], had been discussed by many researchers [10,18], which would not be further analyzed here. This paper adopts the commonly used method, the Tikhonov method. It targets to find a suitable P so as to minimize the following equation: ‖GP u‖2 þ λ‖P‖2 ;

(19)

where ‖‖ denotes the Euclidean norm and the parameter λ is a positive constant called the regularization parameter. λ can adjust the relative size of the residual error ‖GP  u‖ and solution ‖P‖. For the purpose of better approximation, λ should be smaller, whereas to obtain better numerical stability, λ should be larger. This paper uses the L-curve method [19] to find the optimal regularization parameter. With the help of the Tikhonov method and the L-curve method, P in Eq. (18) can be solved. Then f(t) in Eq. (15) can be obtained. In addition, from Eqs. (2)–(18), we see that four parameters need to be determined: lower and upper limits I and J of the integral, level j and order N. They will be quantitatively and qualitatively discussed in Sections 2.3, 3 and 4. 2.3. Method analysis Using the wavelet multi-resolution analysis, a signal x(t) can be reconstructed at level j as follows: D0 xðtÞ þ D1 xðtÞ þ ⋯Dj xðtÞ þ P j xðtÞ; Di xðtÞ ¼

P i xðtÞ ¼

1

∑

k ¼ 1 1

∑

k ¼ 1

(20)

di ðkÞ2  2=i φð2  i t  kÞ;

(21)

αi ðkÞ2  2=i ϕð2  i t kÞ;

(22)

where Dix(t) and Pix(t) are detailed and approximation components of the original signal x(t), respectively. The two spaces associated with Dix(t) and Pix(t) are defined, respectively, to be Wi and Vi. Then V i ¼ Wi þ 1  Wi þ 2  ⋯  Wj  V j: In the definition given by Mallat [9], the multi-resolution analysis has such a property   1 [ V j ¼ L2 ðRÞ: lim V j ¼ Closure j-  1

j-  1

(23)

(24)

This property, in other words, expresses that when j - 1, namely 2  j-1, the approximation component will converge to itself lim jP j xðtÞ  xðtÞj ¼ 0:

j-  1

(25)

In this paper, for convenient calculation, j in Eq. (25) is replaced by –j, which means that with the increase of j, approximation will converge to its original signal. Actually, according to the analyses carried out below, level j does not need to be so large that the approximation component can well converge to the original. Some good properties of wavelets are reflected in this method because of the use of wavelets. db Wavelets are compactly supported, which means that this method has better localization ability. Owing to the use of orthogonal wavelets, similar to the frequency method, orthogonal bases are used. Another good property is multi-resolution analysis. With increasing level j, wavelets are more compactly supported. This means that wavelets have a shorter duration in the time domain, which is good for the analysis of high-frequency transient signals. The number of weight coefficients, namely the difference between I and J, needs to be confirmed in Eq. (3). It needs to be pointed out that in dyadic wavelets, the unit of k is second and db orthogonal wavelet bases have orthogonality only when the shifting of wavelets is an integral multiple of k. For J, we should ensure the starting time of the last basis function is not larger than the ending time of the signals. For I, the ending time of first basis should not be smaller than the starting time of the signals. Under such conditions, basis functions can cover the whole time of signals. Hence, at level j, we have 2  j t l rt end ;

(26)

2  j ðt f þ 2N 1Þ Z t begin ;

(27)

where tl is the starting time of the last basis function, tf is the starting time of the first basis function, tbegin is the starting time of the signals and tend is the ending time of the signals. Eqs. (26) and (27) can be solved as follows: t l r 2j t end ;

(28)

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

385

t f Z2j t begin ð2N 1Þ;

(29)

J ¼ maxðt l Þ

J A Z;

(30)

I ¼ minðt f Þ

I A Z:

(31)

3. Numerical studies As shown in Fig. 1, a 4dofs spring-mass-damping system is modeled with the software of ANSYS. Mass M1 ¼M4 ¼1, M2 ¼M3 ¼ 2. Stiffness coefficients K1 ¼K5 ¼1, K2 ¼K3 ¼K4 ¼2. Damping coefficients C1 ¼C2 ¼C3 ¼C4 ¼C5 ¼1. A sinusoidal excitation f1 and an impulse excitation f4, described by Eqs. (32) and (33), are applied simultaneously on M1 and M4. Responses are acquired on M2 and M3. The length of each data is 3000 and the time step is 0.0078 s. db3 Wavelet is used, and not the Tikhonov method. f 1 ¼ sin ð2pi0:2tÞ;

(32)

2

f 4 ¼ e  ðt  5=0:05Þ :

(33)

The wavelet method is used to identify these two loads. To analyze multi-resolution feature, j varies from 1 to 7. Corresponding I and J are illustrated in Table 1. Two evaluation methods, the correlation coefficient method and the relative error method, as described in Eqs. (34) and (35), are adopted to show the difference between identified loads and original loads. Results are shown in Fig. 2. corr ¼

covðf o ; f i Þ ; sf o sf i

relative error ¼

(34)

norm2 ðf i f o Þ ; norm2 ðf o Þ

(35)

where fo is the original loads, fi is the identified loads, cov is the covariance, s is the standard derivation and norm2 is the 2-norm. To identify fast varying loads, such as steep front edge, back edge, sharp pulse and so on, a high resolution is required, so as to ‘capture' the short duration of the fast varying loads. With the increase of j, time duration of the basis functions decreases and the resolution is improved. From Fig. 2 it can be seen that the identification accuracy of f4 is improved significantly with the increase of j. A good identification result of f4 is achieved at level 5 (j¼5), as shown in Fig. 3. The identified load matches the original one well. The correlation coefficient and the relative error are 0.9989 and 4.7%, respectively. However, the resolution is not very important for identifying slowly varying signals. As seen in Fig. 2, the sinusoidal load f1has excellent identification results at any level. When j is smaller than 5, its identification accuracy is much better than f4, which requires a high resolution. With the increase of j, the relative error of f1 decreases first and then increases (it should be pointed out that the change is small, only 4%). The identification of slowly varying signals does not require a high resolution, but if the resolution is too bad, the basis functions cannot ‘capture' enough signal details. This leads to a poor identification result. This is why the relative error of level 1 is bigger than those of levels 2–4. However, when the level keeps increasing from 5 to 7, the relative error increases instead. This is because the number of weight coefficients is so large that it leads to over-reconstructions of the signals. Over-reconstructions make signals oscillate, as marked with the red circle in Fig. 3. For single input, it is not difficult to choose the decomposition level by the characteristic of the signal. However, for the case of multiple inputs, like in this example where fast and slowly varying loads are applied simultaneously, it is not easy to

Fig. 1. A 4dofs spring-mass-damp system model.

Table 1 The size of lower and upper limits of the integral I and J changing with j. Level j

1

2

3

4

5

6

7

Lower limit of the integral I Upper limits of the integral J The number of weight coefficients J I þ1

4 46 51

4 93 98

4 187 192

4 374 379

4 748 753

4 1497 1502

4 2995 3000

386

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

Fig. 2. Identification results changing with j for two loads: solid circle: f1, empty circle: f4.

Fig. 3. Comparison of original and identified loads when j¼5: (a) original load f1, (b) original load f4, (c)identified load f1 and (d) identified load f4 (I=  4, J=748).

choose a level. As a matter of experience, fast varying signals are more sensitive to resolution than slowly varying signals. Therefore, when dealing with such a problem, a higher level is recommended. In this example, level 5 is chosen. The correlation coefficient and the relative error of f1 and f4 are 0.9992, 4% and 0.9989, 4.7%, respectively. From Eqs. (28) and (29), it can be seen that level j directly influences the number of weight coefficients. If j is low, a small number of weight coefficients is used. Conversely, if j is high, more weight coefficients are used. Owing to different characteristic of the system and loads in experiments, it is not easy to determine an accurate j. Here the range of j is given. The number of weight coefficients has a relation with the length of data (Len) shown as follows: J  I þ1 ¼ kLen

(36)

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

387

where k A[0,1]. When in use, too many weight coefficients will lead to an over-reconstruction whereas a small number is incomplete and will not provide a good identification result. As a matter of experience, k A[0.2,0.7] is a more appropriate range. For example, in Section 3, k is chosen as 0.25 (j¼5). In Section 4.2, k is chosen as 0.5 (j¼11). In other experiments not illustrated in this paper, k belongs to [0.2,0.7]. Therefore, from Eqs. (28)–(31) and (36) and the range of k, the range of j can be calculated. Furthermore, if a priori information of the loads is known, follow the rule: fast varying signals need a higher level whereas slowly varying signals do not. 4. Simulation and experiment 4.1. Numerical simulation A numerical simulation model of cantilever beam is built considering the real beam experimental system to be introduced in Section 4.2, as shown in Fig. 8. From right to left is node 1 to 4, as shown in Fig. 4. The impulse excitation in Eq. (33) and the sinusoidal excitation in Eq. (32) are exerted on node 1 and node 4 simultaneously. Responses are acquired at nodes 2 and 3. The length of each data is 3000 and the time step is 0.0024 s. The wavelet method is used to identify the two loads. db6 Wavelet and the Tikhonov method are used. As analyzed in Section 3, a higher decomposition level should be adopted when identifying fast varying signals. From Fig. 5 it can be clearly seen that with the increase of j, the relative error of f1 decreases, namely the identification accuracy is improved gradually. f4 is a slowly varying signal and its identification results become good quickly with the increase of j (except for level 7; a mutational problem may arise when dealing with inverse problem and it is not a common case). In level 9, relative errors of the two loads are only 6.89% and 2.57%, respectively. When j keeps increasing from 9 to 10, the accuracy of f4 decreases due to the oscillation. The relative error at level 10 is 3.7 times than at level 9. The decomposition level 9 is chosen to analyze the sensitivity to noise of this method. The orthogonal wavelet transform has an important characteristic: that the stationary white noise transformed by orthogonal wavelet is still stationary white noise. Therefore, for additive white noise model yðtÞ ¼ yo ðtÞ þ nðtÞ

(37)

where yo(t) is a signal without pollution, n(t) is noise and y(t) is the polluted signal. Applying the wavelet transform to Eq. (37), the correlation of y(t) will be reduced dramatically and its energy will concentrate on a few of the wavelet weight coefficients with large amplitude whereas n(t) is still white noise, and its energy distributes in the whole wavelet domain with small amplitude. Further verification is performed by numerical simulation. All other conditions being the same, white noise is added with SNR starting from 0 dB to 40 dB to original responses. Identification results are shown in Fig. 6. It can be seen that this method is not sensitive to white noise. Excellent identifications can be achieved when SNR is higher than 15 dB. Relative errors are only 6.92% and 8.88%, respectively, when SNR is 15 dB. For laboratory experimental data, the calculated IRF or Frequency Response Function (FRF) may be polluted due to input or output noise. In this case, the identification accuracy will drop dramatically. Therefore, it is important to verify the sensitivity of a method to a polluted system. From Eqs. (6)–(11), it can be seen that G is transformed twice by wavelet from h(t). Therefore, its energy highly concentrates on a few of the coefficients with big amplitude. Further verification is performed by numerical simulation. All other conditions being the same, original IRF is added with white noise with SNR starting from 30 dB to 60 dB.

Fig. 4. Node location of the beam model.

350.00%

relative error

300.00% 250.00% 200.00% 150.00% 100.00% 50.00% 0.00%

j

6

7

8

f1 112.96%

5

87.39%

69.93%

44.61%

6.89%

9

2.84%

10

f4 137.88%

17.12%

325.94%

9.13%

2.57%

9.41%

Fig. 5. Relative errors of identified loads and original loads when j increases from 5 to 10.

388

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

60.00%

relative error

50.00% 40.00% 30.00% 20.00% 10.00% 0.00% SNR (dB)

40

35

30

25

20

15

10

5

0

f1 6.89% 6.89% 6.89% 6.90% 6.92% 6.99% 7.20% 7.94% 9.74% f4 2.61% 2.67% 3.17% 3.87% 5.50% 8.88% 14.36% 29.23% 49.67% Fig. 6. Relative errors of identification results with SNR (dB) of responses starting from 0 dB to 40 dB.

250.00%

relative error

200.00% 150.00% 100.00% 50.00% 0.00%

55

50

40

35

30

f1 6.89%

6.89%

6.90%

6.91%

7.04%

7.48%

30.00%

f4 2.60%

2.79%

3.71%

5.23%

11.34% 20.88% 216.63%

SNR

60

45

Fig. 7. Relative errors of identification results with SNR of system starting from 30 dB to 60 dB.

Table 2 The first two smallest relative errors and corresponding orders N at levels 8–10. The first two smallest relative errors

j 8

Corresponding N

N ¼4 N¼ 8

9 8.83% 11.37%

N¼ 7 N ¼ 11

10 3.75% 4.59%

N¼ 6 N ¼3

6.12% 7.31%

Identification results are shown in Fig. 7. It can be seen that the identification is acceptable when SNR is higher than 35 dB. In these conditions, accuracy is not affected too much by the noisy system. Accuracy drops dramatically and loads cannot be identified when SNR is lower than 35 dB. According to Eqs. (2)–(18), four parameters need to be determined: order N, level j and lower and upper limits I and J of the integral. I and J can be determined by Eqs. (31) and (30). j can be chosen following the rules proposed in Section 3. Now there only remains N to be determined. N represents the length of the wavelet basis function. The length can be calculated by ð2N  1Þ2  j 2ITER ;

(38)

where ITER is the iterations of the wavelet basis function and is used to generate wavelet function. All other conditions being the same and N ranging from 3 to 15, relative errors are calculated with j starting from 8 to 10. The orders N corresponding to the first two smallest relative errors are chosen at the every level. Results are shown in Table 2. It can be seen that with the increase of j, smaller N can be chosen. For example, N¼8 can be chosen at level 8, N¼7 at level 9 and N¼6 at level 10, as marked with bold in Table 2. This is because a larger j represents more wavelet basis functions, so the length can decrease. It should be remembered that although N affects identification accuracy, it is not a key factor. For example at level 9, the maximum difference between relative errors is smaller than 5%. It can be concluded that order N affects accuracy partly, but can be decreased by increasing j.

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

389

4.2. Laboratory experiment A 2-input-2-output cantilever beam system is set up as shown in Fig. 8. The size of the rectangular cantilever beam is 45.5 cm  6 cm  0.5 cm. SINOCERA vibration exciters are used to impose the loads. DYTRAN 3056b1 piezoelectric acceleration transducers, PAK Mobile MKII:MF03-M128data acquisition equipment and IBM T61 laptop are used to acquire the data. KISTLER force-hammer is used to acquire the FRF of the beam by the force-hammer excitation method. Then IRF can be calculated. As illustrated in Fig. 4, nodes from the right to left are Nos. 1–4. An impulse load is exerted on node 1 by a force-hammer. A sinusoidal load with a frequency of 15 Hz and a peak of 15 N is exerted on node 4 by a vibration exciter. Responses are collected on nodes 2 and 3. The proposed method and the frequency method are used to identify loads. Sample frequency is 4096 Hz. The length of each data is 3000. db6 Wavelet is used. J is 1499 and I is  10. The Tikhonov filter and L-curve are used in both methods. Relative errors of the original loads and the identified loads are shown in Table 3, where it can be seen that the wavelet method has better identification accuracy. It can be seen from Fig. 9 that the identified loads by the wavelet method (red thick line) match the originals (blue thin line) well. Fig. 10 is the comparison of the zoomed-in impulse load. It can be seen that the original and the identified by the wavelet method match well. In particular, the relative error of the impact duration is smaller than 0.2%. Peak is another key parameter. For the original peak of 120.42 N, the one identified using the wavelet method is 120.04 N; however, the one identified by the frequency method is 82.29 N. Relative error of the peak by the wavelet method with respect to the original peak is only 0.32%, whereas the latter is up to 32.67%. In the unloading phase, owing to the inverse of ill-conditioned matrix and the effect of noise, the load identified by the wavelet method has slight oscillation, e.g. from 0.28 s to 0.36 s, as shown in Fig. 10. Besides, at the end (from 0.7 s) of the impulse load, there is an up-tail as shown in Fig. 9. These two factors increase the overall relative error of the wavelet method, which is 13.64% as listed in Table 3. Compared with the wavelet method, impulse load identified by the frequency method has a severely low-frequency oscillation during the whole duration, as shown in Fig. 9. Its oscillation period approaches 0.1 s, which is exactly the period of sinusoidal load. It means that the impulse load identified by the frequency method is affected by the sinusoidal load. Its relative error is up to 47.98%. For the sinusoidal load, both methods identify the period accurately. For the real peak-to-peak value of 30 N, the one identified by the wavelet method is 30.22 N; however, the other identified by the frequency method is 21.42 N. Relative error of the peak-to-peak value obtained by the wavelet method with respect to the original is 0.73% whereas for the latter it is up to 28.6%. Moreover, the sinusoidal load identified by the frequency method has an obvious impact at 0.27 s, as shown in Fig. 9. It corresponds exactly to the time when the impulse load is exerted. It means that the sinusoidal load identified by the frequency method is affected by impulse load. Its relative error is up to 21.35%. Compared with the frequency method, the sinusoidal load identified by the wavelet method does not impact at 0.27 s. Owing to the excellent adjustable resolution, the two loads have little effect on each other. In this analysis, a relatively higher level, 11, is chosen. Higher level means a shorter wavelet duration, which correctly distinguishes sinusoidal loads from the impulse at 0.27 s. The resolution of the frequency method is not adjustable. Therefore, it cannot distinguish fast varying loads from slowly varying loads, which leads to the interaction. It should be pointed out that, as aforementioned, for the decomposition level in the wavelet method, it is not true that the higher the better. More time and more computer resources will be needed for calculation; in the meantime, the oscillation problem may arise. As shown in Fig. 9, the impulse load identified by the wavelet method has an up-tail at the end of it.

Fig. 8. Cantilever beam experiment.

Table 3 Relative error of original and identified loads using wavelet method and frequency method. Relative error

f1 Impulse load (%)

f4 Sinusoidal load (%)

Frequency method Wavelet method

47.98 13.64

21.35 4.73

390

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

Fig. 9. Time domain waveforms of original loads (blue), identified by the wavelet method (red) and the frequency method (black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Zoomed-in impulse load f1 in time domain.

5. Conclusions A novel load identification method based on wavelet multi-resolution analysis is proposed. The method transforms the convolution relation between responses and loads in time domain into the linear multiplicative relation between system responses and wavelet responses in the wavelet domain. Loads can be identified as long as the linear relation is solved without knowing the initial values. Compared with the frequency method, which has an analogous form, the proposed method has a better identification ability and excellent robustness to response noise and system noise. Moreover, this paper proposes qualitative and quantitative selection rules about parameters that affect the accuracy of this method, such as level j, order N and lower and upper limits I and J of integral. In laboratory experiments, the proposed method precisely identified

Z. Li et al. / Journal of Sound and Vibration 333 (2014) 381–391

391

the sinusoidal and the impulse loads exerted on a cantilever beam simultaneously. The results are much better than those obtained by the frequency method. Acknowledgments This research was supported by Natural Science Foundation of China (Grant no. 11272170) and the Tsinghua University Initiative Scientific Research Program (2011Z08137). References [1] A.J. Holzer, A technique for obtaining compressive strength at high strain rates using short load cells, International Journal of Mechanical Sciences 20 (1978) 553–560. [2] F.D. Bartlett, W.G. Flannelly, Modal verification of force determination for measuring vibratory loads, Journal of American Helicopter Society 24 (2) (1979) 10–18. [3] B. Hillary, D.J. Ewnis, The use of strain gauges in force determination and frequency response function measurement, Proceedings of the 2nd IMAC, Orlando, FL, USA, 1984, pp. 627–634. [4] J.F. Doyle, Further developments in determining the dynamic contact law, Experimental Mechanics 24 (4) (1984) 265–270. [5] J.F. Doyle, Experimentally determining the contact force during the transverse impact of an orthotropic plate, Journal of Sound and Vibration 118 (3) (1987) 441–448. [6] K.K. Stevens, Force identification problems—an overview, Proceedings of the 1987 Society of Experimental Mechanics Spring Conference on Experimental Mechanics, Houston, USA, June 1987, pp. 838–844. [7] A.N. Thite, D.J. Thompson, The quantification of structure-borne transmission paths by inverse methods—Part 1: improved singular value rejection methods, Journal of Sound and Vibration 264 (2) (2003) 411–431. [8] A.N. Thite, D.J. Thompson, The quantification of structure-borne transmission paths by inverse methods—Part 2: use of regularization techniques, Journal of Sound and Vibration 264 (2) (2003) 433–451. [9] H. Inoue, J.J. Harrigan, S.R. Reid, Review of inverse analysis for indirect measurement of impact force, Applied Mechanics Reviews 54 (2001) 503–524. [10] Y. Liu, W.S. Shepard Jr., Dynamic force identification based on enhanced least squares and total least-squares schemes in the frequency domain, Journal of Sound and Vibration 282 (2005) 37–60. [11] J.F. Doyle, A wavelet deconvolution method for impact force identification, Experimental Mechanics 37 (4) (1997) 403–408. [12] Y.M. Mao, X.L. Guo, Y. Zhao, Experimental study of hammer impact identification on a steel cantilever beam, Experimental Techniques 34 (3) (2010) 82–85. [13] T.P. Nordberg, An iterative approach to road/profile identification utilizing wavelet parameterization, Vehicle System Dynamics 42 (6) (2004) 413–432. [14] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer Academic, Netherlands, 1996. [15] H. Tanaka, Y. Ohkami, Estimation of impact force on a space vehicle based on an inverse analysis technique, Transactions of the Japan Society of Mechanical Engineers Series C (1997) 1172–1178. (63C). [16] C.W. Groetsch, Inverse Problems in the Mathematical Sciences, Vieweg, Braunschweig, 1993. [17] C.R. Vogel, Computational Methods for Inverse Problems, SIAM, Philadelphia, 2002. [18] M.T. Martin, J.F. Doyle, Impact force identification from wave propagation responses, International Journal of Impact Engineering 18 (1) (1996) 65–77. [19] P.C. Hansen, D.P.O.' Leary, The use of the L-curve in the regularization of discrete ill-posed problems, SIAM Journal on Scientific and Statistical Computing 14 (1993) 1487–1503. [20] A. Cohen, I. Daubechies, J. Feauveau, Biorthogonal bases of compactly supported wavelets, Communications on Pure and Applied Mathematics 45 (5) (1992) 485–560. [21] S. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2(R), Transactions of the American Mathematical Society 315 (1989) 69–88. [22] F.E. Gunawan, H. Homma, A solution of the ill-posed impact-force inverse problems by the weighted least squares method, Journal of Solid Mechanics and Materials Engineering (Scopus Indexed) 2 (2) (2008) 188–198. [23] F.E. Gunawan, H. Homma, Impact-force estimation by quadratic spline approximation, Journal of Solid Mechanics and Materials Engineering 2 (2008) 1092–1103. [24] F.E. Gunawan, H. Homma, Y. Kanto, Two-step b-splines regularization method for solving an ill-posed problem of impact-force reconstruction, Journal of Sound and Vibration (Scopus Indexed) 297 (2006) 200–214. [25] Y. Liu, W.S. Shepard Jr., An improved method for the reconstruction of a distributed force acting on a vibrating structure, Journal of Sound and Vibration 291 (2006) 369–387.