# A Comparative Study of Regularization Method in Structure Load Identification

Authors: Bingrong Miao; Feng Zhou; Chuanying Jiang; Xiangyu Chen; Shuwang Yang
Journal: Shock and Vibration
 (2018)

---

## Abstract

This study aims to correlate the vibration data with quantitative indicators of structural health by comparing and validating the feasibility of identifying unknown excitation forces using output vibration responses. First, numerical analysis was performed to investigate the accuracy, convergence, and robustness of the load identification results for different noise levels, sensors numbers, and initial estimates of structural parameters. Then, the laboratory beam structure experiments were conducted. The results show that using the two identification methods Tikhonov (L-curve) and TSVD (GCV-curve) can successfully and accurately identify the different excitation forces of the external hammer. The TSVD based on GCV method has more advantages than the Tikhonov based on L-curve method. The proposed two kinds of load identification procedure based on vibration response can be applied to the safety performance evaluation of the railway track structure in future inverse problems research.


---

## Body

## 1. Introduction

In the past two decades, a large number of researchers have carried out many effective research studies on inverse problems such as structural load identification and damage identification [1–6], and many of the results have been applied to the structural health monitoring technology of structural integrity and safety assessment [7–10]. Most of the load identification techniques use these regularization methods, such as the Tikhonov method and the truncated singular value decomposition (TSVD) method. However, there are few comparative studies on the accuracy, noise resistance, and robustness of these two methods in the identification of engineering structure loads.Thite and Thompson considered the influence of regularizing the matrix inversion, which discussed Tikhonov regularization with the ordinary cross-validation method for selecting the regularization parameter, and they also studied an iterative inversion technique of the method [11]. Jacquelin et al. analyzed the deconvolution problem and utilized general regularization techniques to solve the ill-posed feature [12]. Choi et al. used the Tikhonov regularization method to improve the stability of the load identification results [13]. Wang et al. proposed an improved iteration regularization method for solving linear inverse problems and chose Morozov’s discrepancy principle to verify the posterior regularization parameter [14]. Similarly, Sun et al. combined a new improved regularization operator with L-curve method to overcome the ill condition of load reconstruction and obtained the stable and approximate solutions of inverse problems [15]. With the development of techniques, some new regularization methods have been proposed. Aucejo et al. proposed a method able to proceed without defining any regularization parameter beforehand, and the method is called multiplicative regularization method [16]. Gao and Yu presented a new method based on quadratic programming theory to determine the regularization parameter [17]. Considering the sparsity of force in the time domain, Qiao et al. developed a general sparse regularization method based on minimizing l1-norm of the coefficient vector of basic functions [18]. Different loads on a mechanical system also have different regularization methods. For random dynamic loads, Jia et al. proposed a weighted regularization approach based on the proper orthogonal decomposition (POD), which the regularization parameter was selected by the GCV method [19]. Then, in order to reconstruct the impact force, Gunawan adopted the Wiener filter method and introduced a regularization technique so that the method can provide an optimum solution which balances the aspects of the minimum noise and the minimum residual [20], and Boukria et al. developed an experimental method to identify the impact force on a reinforced concrete slab [21]. In addition, Qiao et al. imposed a novel impact force reconstruction method based on the primal-dual interior point method (PDIPM), which replaced l2-norm by l1-norm [22]. Gunawan extended the Levenberg–Marquardt (LM) method and a variant of the Gauss–Newton methods, in order to solve the ill-posed pulse-type impact-force inverse problem in conjunction with a trust region strategy [23]. A novel method was proposed to identify loads with the combination of the Tikhonov regularization and singular value decomposition (SVD) based on the condition [24].In this paper, the TSVD of regularization parameter selection method based on generalized cross-validation (GCV) method and Tikhonov regularization based on L-curve method are compared and analyzed. The identification accuracy and applicability of the two methods have been summarized through the simulation of the single-input single-output model (SISO) and multiple-input multiple-output model (MIMO), and the identification characteristics of sine wave load and triangle wave load have been demonstrated. Finally, the results of simulation analysis are verified by experiments.


## 2. Construction of Forward Problem

This paper establishes the forward problem of load identification based on the green kernel function. According to the characteristics of the pulse function, the integral of the product of the pulse function with any other function is equal to the function value of the pulse moment. Therefore, any dynamic load of a structure can be expressed as a superposition of the unit pulse signal in the time domain:(1)ft=∫0tftδt−τdτ,where ft is the external load and δ is the Dirac function.It is assumed that the Green function from the loading point to the responding measurement point isz(t), which is the response of the system under the unit pulse load δt. According to the principle of superposition, the relation between excitation and response of system can be written as convolution form (Duhamel integral):(2)zt=∫0tfτHt−τdτ,where Ht−τ denotes the Green kernel function of the structure impulse response, z(t) is the measured responses which can be displacement, velocity, acceleration, strain, and stress. For systems under zero initial conditions, the integral form of Equation (2) in time domain is discretized into linear equations as follows:(3)z1z2⋮zm=h1h2h1⋮⋮⋱hmhm−1⋯h1f0f1⋮fm−1Δt,where Δt is the sampling interval and m is the sampling points.Then, it is easily to written as(4)Z=HF.Equation (4) is represented as a linear discrete equation for a SISO system. When a certain number of impact loads are applied to the different locations simultaneously, the total responses of the system can be expressed as the superposition of the dynamic responses at each known measurement point:(5)Zj=∑i=1MHjiFi,where Hji is the Green function matrix between the points of application of force and sensor location. Therefore, Equation (5) is transformed into a matrix form as follows:(6)Z1Z2⋮ZN=H11H12⋯H1MH21H22⋯H2M⋮⋮⋮⋮HN1HN2⋯HNMF1F2⋮FM,where M and N are the number of points of application and the points of measurement, respectively, and it is needed to ensure that N ≥ M. Equation (6) expresses the linear discrete equation model for MIMO system. In fact, the solution F of the linear system is very sensitive to data errors. And an attempt to solve directly the equation will likely lead to an erroneous solution. Then regularization is used to acquire a meaningful solution.


## 3. Regularized Method

### 3.1. Tikhonov Regularization

The Tikhonov regularization method can be expressed as an optimization problem as follows [25]:(7)minf∈RnHf−z22+α2f22,where ⋅22 is the square of the 2-norm for vector αα>0 is the regularization parameter.The more exact representation of the objective function in Equation (7) can be expressed as(8)fTHTH+α2If−2zTHf+zTz.According to the concept that the gradient of the objective function is equal to zero, the least squares solution of Equation (8) is written as(9)HTH+α2If=HTz.Deriving from Equation (9), the identified load can be obtained by(10)F=HHH+λI−1HHZ,where I is the unit matrix and H indicates Hermitian transpose.


### 3.2. The TSVD Method

The method for dealing with ill-conditioned matrices is to derive new problems by using the rank deficient coefficient matrix of rank deficit. The rank deficient matrixH (which can be deduced by SVD) is approximated by k, the nearest matrix of Hk, which is extended by the decomposition of the truncated singular value of the rank K in Equation (6). Hk can be written as(11)Hk=∑i=1kuiσiviT,k≤n,where σi is the singular values of the matrix H, and the vectors ui,vi are the left and right singular vectors of the matrix H. The truncated singular value decomposition (TSVD) method is used to solve the following equation of the problem [26]:(12)minHkf−z2.The solution to Equation (12) can be written as(13)fk=∑i=1nuiTzσivi.The solution offk is the only solution to the regularization of the numerical zero space in the matrix H without components. In order to avoid excessive amplification of the identified results' disturbance error, the right end of the resulting generalized solution formula is truncated by using the method. In addition, the parts of the previous k corresponding to greater singular values are reserved, and the smaller singular values will be filtered.


## 3.1. Tikhonov Regularization

The Tikhonov regularization method can be expressed as an optimization problem as follows [25]:(7)minf∈RnHf−z22+α2f22,where ⋅22 is the square of the 2-norm for vector αα>0 is the regularization parameter.The more exact representation of the objective function in Equation (7) can be expressed as(8)fTHTH+α2If−2zTHf+zTz.According to the concept that the gradient of the objective function is equal to zero, the least squares solution of Equation (8) is written as(9)HTH+α2If=HTz.Deriving from Equation (9), the identified load can be obtained by(10)F=HHH+λI−1HHZ,where I is the unit matrix and H indicates Hermitian transpose.


## 3.2. The TSVD Method

The method for dealing with ill-conditioned matrices is to derive new problems by using the rank deficient coefficient matrix of rank deficit. The rank deficient matrixH (which can be deduced by SVD) is approximated by k, the nearest matrix of Hk, which is extended by the decomposition of the truncated singular value of the rank K in Equation (6). Hk can be written as(11)Hk=∑i=1kuiσiviT,k≤n,where σi is the singular values of the matrix H, and the vectors ui,vi are the left and right singular vectors of the matrix H. The truncated singular value decomposition (TSVD) method is used to solve the following equation of the problem [26]:(12)minHkf−z2.The solution to Equation (12) can be written as(13)fk=∑i=1nuiTzσivi.The solution offk is the only solution to the regularization of the numerical zero space in the matrix H without components. In order to avoid excessive amplification of the identified results' disturbance error, the right end of the resulting generalized solution formula is truncated by using the method. In addition, the parts of the previous k corresponding to greater singular values are reserved, and the smaller singular values will be filtered.


## 4. Methods for Selecting the Tikhonov Regularization Parameter

### 4.1. L-Curve Criterion

H
F
−
Z and F are both functions of regularization parameter λ. The minimum between the regularization solution and the residual norm has been weighed through the L-curve method. When the regularization parameter changes, it refers to the change of the demonstration number and the norm of the residual norm. The solution norm and the residual norm are well balanced at the inflection point of the L-curve (the maximum curvature of the L-curve), and the regularization parameter is the optimal regularization parameter. Therefore, the reasonable value of the point obtained by the L-curve criterion is located at the inflection point of the curve. If the relationships ρλ=HF−Z,ηλ=F are assumed, the curvature of the L-curve is calculated as follows:(14)Kλ=ρ′η″−ρ″η′ρ′2+η′23/2.Unfortunately, since optimal regularization parameter value search range is too large, the computational cost will be high or the L-curve graph is not satisfied.


### 4.2. Generalized Cross-Validation (GCV) Method

The basic idea of generalized cross-validation (GCV) means if the arbitrary elementzi of the z is ignored, the corresponding regular resolution will be predicted better, and the regularization parameter is chosen as the orthogonal transformation of z. Therefore, choosing the regularization parameter is required to minimize the GCV function.(15)G=Hfreg−z22traceIm−HHI2.Equation (14) can also be expressed in discrete form.(16)Gk=∑i=1NuiTz/σi2+k22∑i=1N1/σi2+k22,where HI denotes the regularized solution produced by freg multiplied with z (i.e., freg=HIz) and k is the regularization parameter. A balance between the disturbance error and the regularization error is actually searched by using the GCV method, and the inflection point of the curve is produced.


## 4.1. L-Curve Criterion

H
F
−
Z and F are both functions of regularization parameter λ. The minimum between the regularization solution and the residual norm has been weighed through the L-curve method. When the regularization parameter changes, it refers to the change of the demonstration number and the norm of the residual norm. The solution norm and the residual norm are well balanced at the inflection point of the L-curve (the maximum curvature of the L-curve), and the regularization parameter is the optimal regularization parameter. Therefore, the reasonable value of the point obtained by the L-curve criterion is located at the inflection point of the curve. If the relationships ρλ=HF−Z,ηλ=F are assumed, the curvature of the L-curve is calculated as follows:(14)Kλ=ρ′η″−ρ″η′ρ′2+η′23/2.Unfortunately, since optimal regularization parameter value search range is too large, the computational cost will be high or the L-curve graph is not satisfied.


## 4.2. Generalized Cross-Validation (GCV) Method

The basic idea of generalized cross-validation (GCV) means if the arbitrary elementzi of the z is ignored, the corresponding regular resolution will be predicted better, and the regularization parameter is chosen as the orthogonal transformation of z. Therefore, choosing the regularization parameter is required to minimize the GCV function.(15)G=Hfreg−z22traceIm−HHI2.Equation (14) can also be expressed in discrete form.(16)Gk=∑i=1NuiTz/σi2+k22∑i=1N1/σi2+k22,where HI denotes the regularized solution produced by freg multiplied with z (i.e., freg=HIz) and k is the regularization parameter. A balance between the disturbance error and the regularization error is actually searched by using the GCV method, and the inflection point of the curve is produced.


## 5. Numerical Example and Test Validation

As shown in Figure1, the example model is a beam with one end fixed to simulate a steel cantilever beam. The material properties of the beam are as follows: Young’s modulus is 2.1 × 1011 Pa, density is 7.8 × 103 Kg/m3, and Poisson’s ratio is 0.3. The size of the beam is 1 × 0.05 × 0.005 m. The beam is divided into ten elements, marked elements 1 to 10. Then, the number of finite element model nodes is expressed from the mark number of one to eleven. Numerical simulation is used to verify the efficiency of two proposed force identification method as well as compared to the identified results. Furthermore, the identified results of different loads applied on the same structure are analyzed.Figure 1
The finite element model of cantilever beam.For the purpose of studying the effect of noise on identified results, the noise is added to the response data.(17)z=zcal+lp⋅stdzcal⋅nnoise,where z is the response of load identification and zcal represents the response obtained by simulation. stdzcal denotes the standard deviation of the computed response. lp and nnoise are the noise level and white noise sequence with the same length as zcal, which is a mean value of 0 and a variance of 1. And the relative error between the real load and the identified load is defined as follows:(18)RE=fidentify−frealfreal×100%.


### 5.1. Load Identification of Single-Input Single-Output System (SISO System)

A sinusoidal loadft=10sin100πt0≤t≤0.12 is applied at node 11. The direction of the load is perpendicular to the Z axis of the beam, and a sampling period is Δt=0.0005s. The identified load is obtained by using the response of node 10. At the noise levels 1%, 2%, and 5%, the two identification methods mentioned above are used to identify the loads. Figure 2 shows the results of the two methods when the noise level is 1%. The overall relative error of the two methods are shown in Table 1 at 0%,1%, 2%, and 5% noise levels. Among them, the noise level is 0%, which means no noise.Figure 2
Comparison of results using (a) Tikhonov regularization based on L-curve and (b) TSVD regularization based on GCV.
(a)
(b)Table 1
Results of load identification of SISO system.
Noise level (%)
TSVD + GCV (RE/%)
Tikhonov + L-curve (RE/%)
0
7.14
7.36
1
8.20
7.35
2
11.35
8.86
5
24.44
16.78As shown in Figure2, the two methods can accurately identify the load applied to the system. Obviously, it can be found clearly that the identification accuracy of Figure 2(a) is less than that of Figure 2(b) in the initial period of time. It refers that the two methods of the identified results show the analogous situation in the absence of noise in Table 2. It can be found that the identification error of Tikhonov based on the L-curve method is greater than the TSVD method based on GCV in the case of noise. In addition, the latter method also has weak sensitivity to noise and has better adaptability. In fact, periodic synthetic signal gives naturally better results than the ones obtained from nonperiodic signal for inverse analysis, because SVD (TSVD) is based on Fourier series.Table 2
Comparison of experimental and simulated natural frequencies of test specimens.
Modal orders
Natural frequency value of test (Hz)
Natural frequency value of simulation (Hz)
Error (%)
1
20.42
20.271
0.7
2
55.66
56.933
2.2
3
184.40
186.49
1.12
4
268.93
280.64
4.17


### 5.2. Load Identification of Multiple-Input Multiple-Output System (MIMO System)

To study the differences between the two methods, a sinusoidal load in theZ direction and a triangular load in the Y direction are applied at the nodes 11 and 7, respectively, and the responses of the nodes 10 and 6 are measured. Among them, the sinusoidal load f1t and the triangular load f2t are described as follows:(19)f1t=10sin2πt/td,0≤t≤4td,0,t>4td,f2t=5t/tf,0≤t≤tf,−5t/tf+10,tf≤t≤3tf,5t/tf−20,3tf≤t≤5tf,−5t/tf+30,5tf≤t≤7tf,5t/tf−40,7tf≤t≤8tf,0,t>8tf,where td is a period of sinusoidal load and tf is a quarter period of triangle load. In this case, the sampling time interval Δt is 0.001 s, 0.0005 s, and 0.00025 s, respectively. Clearly, Figures 3 and 4 show the identification result of the above two kinds of loads in the different sampling time. Due to the great influence of high levels of noise on load identification results of MIMO system, the main features of the load to the maximum amplitude are concerned. And at the noise level 1%, the identified results of the two methods are shown in Figures 5 and 6.Figure 3
The sinusoidal load identification results with no noise and different sampling time.Figure 4
The triangle load identification results with no noise and different sampling time.Figure 5
The sinusoidal load identification results with 1 percent.Figure 6
The triangular load identification results with 1 percent.According to Figures3 and 4, it can be found that the sampling time interval has a great influence on the identification of triangle wave load, while the sine wave load is not. It means that the sampling time required for the triangular wave load is shorter than the sinusoidal load. In addition, the TSVD method has a more accurate effect than the Tikhonov method, and the former method is not sensitive to noise shown in Figures 5 and 6. Simultaneously, it can be found that the actual noise has a great influence on the identification results.


## 5.1. Load Identification of Single-Input Single-Output System (SISO System)

A sinusoidal loadft=10sin100πt0≤t≤0.12 is applied at node 11. The direction of the load is perpendicular to the Z axis of the beam, and a sampling period is Δt=0.0005s. The identified load is obtained by using the response of node 10. At the noise levels 1%, 2%, and 5%, the two identification methods mentioned above are used to identify the loads. Figure 2 shows the results of the two methods when the noise level is 1%. The overall relative error of the two methods are shown in Table 1 at 0%,1%, 2%, and 5% noise levels. Among them, the noise level is 0%, which means no noise.Figure 2
Comparison of results using (a) Tikhonov regularization based on L-curve and (b) TSVD regularization based on GCV.
(a)
(b)Table 1
Results of load identification of SISO system.
Noise level (%)
TSVD + GCV (RE/%)
Tikhonov + L-curve (RE/%)
0
7.14
7.36
1
8.20
7.35
2
11.35
8.86
5
24.44
16.78As shown in Figure2, the two methods can accurately identify the load applied to the system. Obviously, it can be found clearly that the identification accuracy of Figure 2(a) is less than that of Figure 2(b) in the initial period of time. It refers that the two methods of the identified results show the analogous situation in the absence of noise in Table 2. It can be found that the identification error of Tikhonov based on the L-curve method is greater than the TSVD method based on GCV in the case of noise. In addition, the latter method also has weak sensitivity to noise and has better adaptability. In fact, periodic synthetic signal gives naturally better results than the ones obtained from nonperiodic signal for inverse analysis, because SVD (TSVD) is based on Fourier series.Table 2
Comparison of experimental and simulated natural frequencies of test specimens.
Modal orders
Natural frequency value of test (Hz)
Natural frequency value of simulation (Hz)
Error (%)
1
20.42
20.271
0.7
2
55.66
56.933
2.2
3
184.40
186.49
1.12
4
268.93
280.64
4.17


## 5.2. Load Identification of Multiple-Input Multiple-Output System (MIMO System)

To study the differences between the two methods, a sinusoidal load in theZ direction and a triangular load in the Y direction are applied at the nodes 11 and 7, respectively, and the responses of the nodes 10 and 6 are measured. Among them, the sinusoidal load f1t and the triangular load f2t are described as follows:(19)f1t=10sin2πt/td,0≤t≤4td,0,t>4td,f2t=5t/tf,0≤t≤tf,−5t/tf+10,tf≤t≤3tf,5t/tf−20,3tf≤t≤5tf,−5t/tf+30,5tf≤t≤7tf,5t/tf−40,7tf≤t≤8tf,0,t>8tf,where td is a period of sinusoidal load and tf is a quarter period of triangle load. In this case, the sampling time interval Δt is 0.001 s, 0.0005 s, and 0.00025 s, respectively. Clearly, Figures 3 and 4 show the identification result of the above two kinds of loads in the different sampling time. Due to the great influence of high levels of noise on load identification results of MIMO system, the main features of the load to the maximum amplitude are concerned. And at the noise level 1%, the identified results of the two methods are shown in Figures 5 and 6.Figure 3
The sinusoidal load identification results with no noise and different sampling time.Figure 4
The triangle load identification results with no noise and different sampling time.Figure 5
The sinusoidal load identification results with 1 percent.Figure 6
The triangular load identification results with 1 percent.According to Figures3 and 4, it can be found that the sampling time interval has a great influence on the identification of triangle wave load, while the sine wave load is not. It means that the sampling time required for the triangular wave load is shorter than the sinusoidal load. In addition, the TSVD method has a more accurate effect than the Tikhonov method, and the former method is not sensitive to noise shown in Figures 5 and 6. Simultaneously, it can be found that the actual noise has a great influence on the identification results.


## 6. Experimental Validation

The cantilever beam structure applied for force identification is set up as shown in Figure7. The material is structural steel, the length is 1.196 m, the width is 0.06 m, and the thickness is 0.0055 m. In this experimental study, the beam is clamped at one end and the excitation point is selected to be near the end of the beam. For the material properties of the steel beam, modulus of elasticity is 2.1 × 1011 N/m2, density is 7.8 × 103 Kg/m3, and Poisson ratio is 0.3. Four points of the beam are attached by four accelerometers, which are IEPE voltage output piezoelectric acceleration sensor. A shaker is used to generate the sinusoidal force. The measured signals containing the force and acceleration are synchronously recorded by DH9522 data acquisition system. For the purpose of modal frequency, the finite element model of the structure is established and calculated according to the actual information of the geometry, material, and boundary conditions. A modal impact hammer is used to impose impact forces acting on the cantilever beam. The natural frequencies of simulation and experiment are compared as shown in Table 2, and the finite element model is modified according to the experimental model. And the mentioned methods of identification are verified by experiments.Figure 7
Experimental setup: (a) site map; (b) sketch map.
(a)
(b)From the natural frequency of Table2, the relative error of the natural frequency of the test and simulation is relatively small and is less than 5%. It means that the established finite element model is more reliable and can be used for the identification of sinusoidal loads directly.The dynamic response of the structure is measured by experiment, and the Green kernel function between the loading point and the measuring point is calculated by the finite element method. Therefore, the load identification equation based on the Green kernel function has been established. As shown in Figure8, the identified results of the TSVD based on the GCV method and the Tikhonov based the on L-curve method are compared and analyzed. The conclusion of the identified result shows that the identification load reconstructed by the TSVD based on the GCV method has higher accuracy than the Tikhonov based on the L-curve method, and it is in suitable agreement with the actual load. And using the Tikhonov regularization method leads to obtain worse correlation error of identification near the initial moment.Figure 8
The identified results in the experimental.


## 7. Conclusions

This study proposes a comparative research method for the inverse problem of load identification including the TSVD (GCV) method and Tikhonov based on the L-curve method. At the same time, numerical methods and experiments are used to verify the recognition effectiveness and robustness of these two approaches based on vibration response to identify the external excitation force. For these two kinds of regularization methods under different noise levels, the equations based on the established input-output relations for unknown external excitation forces are iteratively estimated. Some detailed conclusions are briefly listed below:(1)
In the absence of noise, both the two mentioned methods of the above identification methods can accurately identify the load imposed on the system.(2)
By virtue of the proposed method, the noise should be considered to the simulation and experiment. The results of numerical simulations and experimental are indicated that the identification accuracy of the TSVD based on the GCV method is much higher than that of the Tikhonov based on the L-curve method. In other words, the former method is less sensitive to noise.(3)
The difference between the two methods is mainly due to the fact that the L-curve generated by the experiment is unobservable, and the optimal point on the L-curve is difficult to localize with accuracy.(4)
The results of different loads in different sampling time can get the conclusion that the identified effect of triangle wave load changes obviously with the change of sampling time. Therefore, the sampling time of linear load such as triangle wave is short; in other words, the sampling frequency could be large.



---

*Source: 9204865-2018-10-15.xml*