# Load Localization and Reconstruction Using a Variable Separation Method

Authors: Jing Zhang; Fang Zhang; Jinhui Jiang
Journal: Shock and Vibration
 (2019)

---

## Abstract

Load identification is very important in engineering practice. In this paper, a novel method for load reconstruction and localization is proposed. In the traditional load localization method, location information is coupled to the impulse response matrix. The inversion of the impulse response matrix leads the process of load localization to be time-consuming. So we propose a variable separation method to separate the load location information from the impulse response matrix. An error optimization function of load histories in different modes is employed to determine the true load location. After locating the external load, the load time history can be easily reconstructed by the measurement responses and determinate impulse response matrix. This method is verified by simulations of a simply supported beam acted by a sine load and an impact separately. An experiment is also carried out to validate the feasibility and accuracy of the proposed method.


---

## Body

## 1. Introduction

Knowledge of external loads and load locations is crucial in various fields, such as structural dynamic design, noise reduction, and fault diagnosis. However, it is usually difficult to directly measure the structure external loads and the locations due to some physical or economical limitations. This is why indirect methods must be developed to identify structural loads by using measured structural responses, such as acceleration, velocity, displacement, and strain.In recent years, with the deepening of the research on structure dynamics, the techniques of load identification have developed rapidly. Currently, there are two mature methods of load identification: frequency-domain methods and time-domain methods. The frequency-domain method converts the kinetic equation into linear equation in the frequency domain. Bartlett and Flannelly [1] firstly used the frequency-domain methods to identify the hub forces in a helicopter model. Starkey and Merrill [2] used direct inversion of the frequency response function (FRF) to determine the load. In their research, they found the frequency response function (FRF) was ill-posed near the resonance zone, and with the number of the loads increasing, the accuracy of the identification result was reduced. Doyle [3] established a waveguide model to describe the dynamic response and employed the spectral analysis to reconstruct the impacting force acting on a two-dimensional bimaterial beam system. Liu and Shepard [4] used enhanced least squares schemes to reduce the random errors of structural response signals. In their researches, a total least squares scheme was also employed to solve the errors associated with the FRF matrix. The applications of the frequency-domain methods are often restricted since sufficient long data are required to apply the Fourier transformation or other harmonic transformation which is often used in frequency-domain methods. Thus, the time-domain method of load identification was developed. The research of time-domain methods started relatively late and still has many problems need to be solved. Desanghere and Snoeys [5] established the time-domain method of load identification and firstly used the modal coordinate transformation to identify the external load. Chang and Sun [6] employed a deconvolution method to reconstruct the time history of structure loads. The load identification is an inverse problem, where the structural properties and responses are known while the external loads need to be determined. However, this inverse problem is mathematically ill-posed in most case, which means that the uniqueness and the stability of the solution are lost. A little noise in the measured responses of the structure, which is inevitable, would result in the great change of the estimated load. To solve this ill-posed problem, some regularization techniques are introduced to stabilize the identified load. Jacquelin et al. [7] compared the efficiency of two widely used regularization techniques, the Tikhonov method and the truncated singular value decomposition (TSVD), by identifying an impact acting on a circular Kirchhoff plate. In their research, they also compared the effect of two-parameter choice criterion, L-curve method, and generalized cross validation (GCV). Besides these two methods, there are also some commonly used regularization methods, such as the modified TSVD [8], the damped singular value decomposition [9], and the iterative regularization methods [10]. Liu et al. [11, 12] introduced several new regularization filter function and confirmed the effectiveness and the accuracy of proposed methods for solving load identification problems. In many researches including those mentioned above, the main purpose is to reconstruct the time history of the excitation and the load location is assumed as known information, which is often not the case and needs to be identified before the reconstruction of load time history. In fact, the process of load identification consists of two separate parts: localization of the external load and reconstruction of the load time history. In order to identify the load location, wavelets were introduced to locate the external loads. Gaul and Hurlebaus [13] employed the wavelet transform to determine the arrival time of impact wave to different sensors and used an optimization method to identify the impact location. As for impact load, there are several localization methods by analyzing wave propagation, such as time difference of arrival (TDOA) [14, 15], direction of arrival (DOA) [16, 17], and wave’s energy loss with traveled distance [18, 19]. Ginsberg and Fritzen [20] created a sample-force-dictionary as the prior knowledge to transform the impact identification into a sparse recovery task. Li and Lu [21] adopted a complex method to determine the location of the impact and then identify the impact history by using a constrained optimization scheme. They [22] used a two-step iterative approach to both localize and reconstruct a single point force acting on a structure. In their research, they proposed the stabilization diagram of identified locations to determine the appropriate regularization range and the true force location. However, these two studies have the same problem: the location information is coupled to the impulse response matrix, which generally needs to be inverted in the process of load identification. Thus, to identify the load position, a considerable amount of matrix inversion is needed and a lot of time will be consumed.In this paper, a new localization method is proposed, which separates the location variable from the impulse response matrix. In this method, a modal decomposition method is adopted and only once matrix inversion needs to be done to obtain several-order modal force. The load position is identified by optimizing an error function of the modal force. Numerical simulation and identification test on a simply supported beam structure are implemented to demonstrate the accuracy and effectiveness of the proposed method. This paper is organized as follows: in Section2, the inverse model for load localization and reconstruction is discussed, with introducing regularization techniques. The applied method for fast localization is also proposed. In Section 3, both numerical simulation and experiment are carried out to validate the proposed method. Finally, concluding remarks are presented in Section 4.


## 2. Load Localization and Reconstruction Scheme

### 2.1. Inverse Model

First, we construct the relationship between a concentrated force and responses of a structure. Without loss of generality, a multiple degree-of-freedom system is considered as the structure model, and the degree of freedom of the system is assumed asN. Generally, the dynamic equilibrium equation of the MDOF system is expressed as(1)Mu¨t+Cu˙t+Kut=Pt,where M, C, and K represent the mass, damping, and stiffness matrixes, respectively. ut is the response vector of the system, and Pt is the external load vector.Under the zero initial condition, the relation between the responses and the load can be written as the form of a convolution:(2)ut=ht∗Pt,where ht is the unit impulse response matrix, and it can be expressed through the following modal superposition:(3)ht=∑r=1NφrφrTmrωr1−ξr2e−ξrωrtsin1−ξr2ωrt,where mr, ξr, and ωr are the rth-order modal mass, damping ratio, and natural circular frequency, respectively. φr is the rth-order mode of vibration.For a concentrated force,sf is assumed as the load location, and ft is the load time history. The response of the ith degree of freedom is obtained:(4)uit=∫0t∑r=1Nhrt−τφriφrsffτdτ,where hrt is the rth-order impulse response function:(5)hr=1mrωr1−ξr2e−ξrωrtsin1−ξr2ωrt.The summation part in equation (4) can be represented by(6)hi,sft=∑r=1Nhrtφriφrsf.The subscript ofhi,sft denotes the ith DOF unit impulse response with load acting on the location sf. The convolution in equation (4) can be taken into a discrete form:(7)uit1uit2⋮uitQ=hi,sft10⋯0hi,sft2hi,sft1⋯0⋮⋮⋱⋮hi,sftQhi,sftQ−1⋯hi,sft1ft0ft1⋮ftQ−1Δt,where Δt denotes the sampling interval, Q denotes the number of sampling points, and ti=iΔti=0,1,2,…,Q. Equation (7) can be simplified as(8)ui=Hi,sf⋅f.For a certain system, the only unknown variables in equation (8) are the load location sf and load time history f. To determine the location sf, two or more measurement points are chosen to obtain corresponding load time history f. A minimum optimization problem is established to find the suitable value of location sf, which makes the error function value of several load vectors f minimum. In this process, there are lots of matrix inversion because the location variable sf is included in the matrix H. Generally, the size of the matrix H is large, due to the number of sampling points Q. Thus, the process of load localization will cost plenty of time which may cause unknown problems during researches and practical engineering situations.


### 2.2. New Load Localization Method

To reduce the time of load localization, a new method is considered which separates the load locationsf from the matrix H. In this method, the number of matrix inversions is reduced significantly so that the load location can be identified faster.From equation (6), hi,sft is expressed as the form of the modal superposition. By substituting equations (6) and (8), the following equation can be derived:(9)ui=∑r=1NφriφrsfHr⋅f,where(10)Hr=hrt10⋯0hrt2hrt1⋯0⋮⋮⋱⋮hrtQhrtQ−1⋯hrt1.Combining the componentφrsf and load time history f, equation (9) is turned into(11)ui=∑r=1NφriHr⋅Sr,where Sr=φrsf⋅fr=0,1,2,…,N. Then, summation in equation (11) can be transformed to a matrix form:(12)ui=Wi⋅S,where(13)Wi=φ1iH1φ2iH2⋯φNiHNQ×NQ,S=S1S2⋯SN1×NQT.From equation (12), we can see that the load location sf is included in the vector S, and the matrix Wi only contains the information of measurement point and properties of the system. Thus, the matrix Wi is defined when the measurement point is settled. However, the vector S cannot be calculated by directly inversing equation (12) because equation (12) is an undetermined system of equation.Thus, to obtain the vectorS uniquely, the modal truncation method is employed and the m modes which contribute the most to the response are chosen. Equation (12) is transformed to(14)ui=Witrun⋅Strun,where(15)Witrun=φ1iH1φ2iH2⋯φmiHmQ×mQ,Strun=S1S2⋯Sm1×mQT.In addition, the responses of other measurement points are necessary. Assuming that responses ofm points are known, then we have(16)Y=Wtrun⋅Strun,where(17)Y=u1Tu2T⋯umTT,Wtrun=φ11H1φ21H2⋯φm1Hmφ12H1φ22H2⋯φm2Hm⋮⋮⋱⋮φ1mH1φ2mH2⋯φmmHm.Equation (16) can be inversed directly because Wtrun is a mQ×mQ matrix. Then Sii=1,2,…,m are calculated, and m load vectors can be obtained:(18)fi=Siφisf,i=1,2,…,m.For the true load locationsf, the m load vectors fii=1,2,…,m are equal. Thus, an optimization function is introduced:(19)ηsf=∑i=1mfi−1m∑r=1mfr2.From equation (19), the optimization function ηsf gets the minimum value when the variable sf is the true load location. Thus, the load localization is transformed to a minimum optimization problem.In this method of load localization, only once matrix inversion is computed, which reduces the operation time greatly and increases the efficiency of load localization. For a certain system, the matrixWtrun is determined when the measurement points, the sampling interval, and the number of sampling points are all defined. Therefore, the load location can be quickly identified, as the responses are measured by the sensors. In addition, the selection of modes is a very important part in this method. The identified location may not be accurate, if a specific mode is not kept. Thus, the selected modes should reflect the response as far as possible. As the number of selected modes is limited by the number of measurement points, the identified result of a load with a wide range of frequencies would not be effective when the number of measuring points is not large enough. Furthermore, the normalized modal shape values must be measured, which is difficult to realize in the real work. This problem needs to be solved with the help of finite element techniques.Due to the modal truncation, the load time history which is obtained in this method may not be accurate enough. Thus, after identifying the load location, the load time history can be reconstructed by equation (8).


### 2.3. Regularization Method

Since the matrixWtrun is generally ill-posed, the direct inversion of equation (16) would not give a stable result. The identified result is sensitive to the errors including signal sampling errors, modal truncation errors, and rounding errors. A little error will cause the direct inverse solution change a great deal.Considering the errors, the measurement response can be expressed as(20)Yδ=Wtrun⋅Strun+δ,where δ is the error data of responses and Yδ is the measurement response including errors data. The singular value decomposition of Wtrun is(21)Wtrun=UΩVT=∑r=1mQurσrvrT,where U=u1,u2,…,umQ and V=v1,v2,…,vmQ are matrixes with orthonormal columns and Ω=Diagσii=1,2,…,mQ has nonnegative diagonal singular values appearing in the nonincreasing order. Thus, equation (20) gives a formulation of the identified result(22)S∗=Wtrun−1Yδ=VDiagσi−1UTYδ=Strun+∑i=1mQσi−1uiTδvi.As shown in equation (22), the ill-posed factor of the load identification model is controlled by using the smaller singular value σi rather than the maximal singular value σ1. With the singular values going down to zero, the error of measurement responses would strongly influence the stability of identified results. To reduce this effect, the regularization method is introduced to filter the small singular values. σi−1 in equation (22) is coupled with a regularization operator pα,σi, and the parameter α can make pα,σi decay to zero when σi approaches zero. Thus, the stable identified result can be obtained:(23)S∗=VDiagpα,σiσi−1UTYδ=∑i=1mQpα,σiσi−1uiTYδvi.The regularization operator is selected as follows:(24)pα,σi=σi2σi2+α.This regularization method has turned into the famous regularization method called the Tikhonov regularization method. The regularization parameterα is generally determined with a certain criterion, such as L-curve criterion and GCV criterion.


## 2.1. Inverse Model

First, we construct the relationship between a concentrated force and responses of a structure. Without loss of generality, a multiple degree-of-freedom system is considered as the structure model, and the degree of freedom of the system is assumed asN. Generally, the dynamic equilibrium equation of the MDOF system is expressed as(1)Mu¨t+Cu˙t+Kut=Pt,where M, C, and K represent the mass, damping, and stiffness matrixes, respectively. ut is the response vector of the system, and Pt is the external load vector.Under the zero initial condition, the relation between the responses and the load can be written as the form of a convolution:(2)ut=ht∗Pt,where ht is the unit impulse response matrix, and it can be expressed through the following modal superposition:(3)ht=∑r=1NφrφrTmrωr1−ξr2e−ξrωrtsin1−ξr2ωrt,where mr, ξr, and ωr are the rth-order modal mass, damping ratio, and natural circular frequency, respectively. φr is the rth-order mode of vibration.For a concentrated force,sf is assumed as the load location, and ft is the load time history. The response of the ith degree of freedom is obtained:(4)uit=∫0t∑r=1Nhrt−τφriφrsffτdτ,where hrt is the rth-order impulse response function:(5)hr=1mrωr1−ξr2e−ξrωrtsin1−ξr2ωrt.The summation part in equation (4) can be represented by(6)hi,sft=∑r=1Nhrtφriφrsf.The subscript ofhi,sft denotes the ith DOF unit impulse response with load acting on the location sf. The convolution in equation (4) can be taken into a discrete form:(7)uit1uit2⋮uitQ=hi,sft10⋯0hi,sft2hi,sft1⋯0⋮⋮⋱⋮hi,sftQhi,sftQ−1⋯hi,sft1ft0ft1⋮ftQ−1Δt,where Δt denotes the sampling interval, Q denotes the number of sampling points, and ti=iΔti=0,1,2,…,Q. Equation (7) can be simplified as(8)ui=Hi,sf⋅f.For a certain system, the only unknown variables in equation (8) are the load location sf and load time history f. To determine the location sf, two or more measurement points are chosen to obtain corresponding load time history f. A minimum optimization problem is established to find the suitable value of location sf, which makes the error function value of several load vectors f minimum. In this process, there are lots of matrix inversion because the location variable sf is included in the matrix H. Generally, the size of the matrix H is large, due to the number of sampling points Q. Thus, the process of load localization will cost plenty of time which may cause unknown problems during researches and practical engineering situations.


## 2.2. New Load Localization Method

To reduce the time of load localization, a new method is considered which separates the load locationsf from the matrix H. In this method, the number of matrix inversions is reduced significantly so that the load location can be identified faster.From equation (6), hi,sft is expressed as the form of the modal superposition. By substituting equations (6) and (8), the following equation can be derived:(9)ui=∑r=1NφriφrsfHr⋅f,where(10)Hr=hrt10⋯0hrt2hrt1⋯0⋮⋮⋱⋮hrtQhrtQ−1⋯hrt1.Combining the componentφrsf and load time history f, equation (9) is turned into(11)ui=∑r=1NφriHr⋅Sr,where Sr=φrsf⋅fr=0,1,2,…,N. Then, summation in equation (11) can be transformed to a matrix form:(12)ui=Wi⋅S,where(13)Wi=φ1iH1φ2iH2⋯φNiHNQ×NQ,S=S1S2⋯SN1×NQT.From equation (12), we can see that the load location sf is included in the vector S, and the matrix Wi only contains the information of measurement point and properties of the system. Thus, the matrix Wi is defined when the measurement point is settled. However, the vector S cannot be calculated by directly inversing equation (12) because equation (12) is an undetermined system of equation.Thus, to obtain the vectorS uniquely, the modal truncation method is employed and the m modes which contribute the most to the response are chosen. Equation (12) is transformed to(14)ui=Witrun⋅Strun,where(15)Witrun=φ1iH1φ2iH2⋯φmiHmQ×mQ,Strun=S1S2⋯Sm1×mQT.In addition, the responses of other measurement points are necessary. Assuming that responses ofm points are known, then we have(16)Y=Wtrun⋅Strun,where(17)Y=u1Tu2T⋯umTT,Wtrun=φ11H1φ21H2⋯φm1Hmφ12H1φ22H2⋯φm2Hm⋮⋮⋱⋮φ1mH1φ2mH2⋯φmmHm.Equation (16) can be inversed directly because Wtrun is a mQ×mQ matrix. Then Sii=1,2,…,m are calculated, and m load vectors can be obtained:(18)fi=Siφisf,i=1,2,…,m.For the true load locationsf, the m load vectors fii=1,2,…,m are equal. Thus, an optimization function is introduced:(19)ηsf=∑i=1mfi−1m∑r=1mfr2.From equation (19), the optimization function ηsf gets the minimum value when the variable sf is the true load location. Thus, the load localization is transformed to a minimum optimization problem.In this method of load localization, only once matrix inversion is computed, which reduces the operation time greatly and increases the efficiency of load localization. For a certain system, the matrixWtrun is determined when the measurement points, the sampling interval, and the number of sampling points are all defined. Therefore, the load location can be quickly identified, as the responses are measured by the sensors. In addition, the selection of modes is a very important part in this method. The identified location may not be accurate, if a specific mode is not kept. Thus, the selected modes should reflect the response as far as possible. As the number of selected modes is limited by the number of measurement points, the identified result of a load with a wide range of frequencies would not be effective when the number of measuring points is not large enough. Furthermore, the normalized modal shape values must be measured, which is difficult to realize in the real work. This problem needs to be solved with the help of finite element techniques.Due to the modal truncation, the load time history which is obtained in this method may not be accurate enough. Thus, after identifying the load location, the load time history can be reconstructed by equation (8).


## 2.3. Regularization Method

Since the matrixWtrun is generally ill-posed, the direct inversion of equation (16) would not give a stable result. The identified result is sensitive to the errors including signal sampling errors, modal truncation errors, and rounding errors. A little error will cause the direct inverse solution change a great deal.Considering the errors, the measurement response can be expressed as(20)Yδ=Wtrun⋅Strun+δ,where δ is the error data of responses and Yδ is the measurement response including errors data. The singular value decomposition of Wtrun is(21)Wtrun=UΩVT=∑r=1mQurσrvrT,where U=u1,u2,…,umQ and V=v1,v2,…,vmQ are matrixes with orthonormal columns and Ω=Diagσii=1,2,…,mQ has nonnegative diagonal singular values appearing in the nonincreasing order. Thus, equation (20) gives a formulation of the identified result(22)S∗=Wtrun−1Yδ=VDiagσi−1UTYδ=Strun+∑i=1mQσi−1uiTδvi.As shown in equation (22), the ill-posed factor of the load identification model is controlled by using the smaller singular value σi rather than the maximal singular value σ1. With the singular values going down to zero, the error of measurement responses would strongly influence the stability of identified results. To reduce this effect, the regularization method is introduced to filter the small singular values. σi−1 in equation (22) is coupled with a regularization operator pα,σi, and the parameter α can make pα,σi decay to zero when σi approaches zero. Thus, the stable identified result can be obtained:(23)S∗=VDiagpα,σiσi−1UTYδ=∑i=1mQpα,σiσi−1uiTYδvi.The regularization operator is selected as follows:(24)pα,σi=σi2σi2+α.This regularization method has turned into the famous regularization method called the Tikhonov regularization method. The regularization parameterα is generally determined with a certain criterion, such as L-curve criterion and GCV criterion.


## 3. Tests and Results

### 3.1. Numerical Simulation

In this section, two numerical simulations of a simply supported beam model are implemented to verify the proposed method. In this beam model, the length is1m long, and the cross section is 6cm×1cm. The material of the beam has the density 7800kg/m3, Young’s module 210GPa and Poisson’s ratio 0.3. The beam is evenly divided into 1000 elements, and there are 1001 nodes in total. The first two-order natural frequencies of the beam are f1=23.5Hz and f2=94.1Hz. Two displacement sensors are installed at x1=0.23m and x2=0.45m from one simply supported end, and the beam model is excited by a sine load and an impact separately, as shown in Figure 1.Figure 1
Beam simulation model.


#### 3.1.1. Identification of a Sine Load

The beam is excited by a loadF1, which is shown in Figure 2. F1 is a sine load with frequency f=60Hz acting at the point a1=0.35m from one simply supported end. The measurement responses of the two points x1 and x2 are computed through the analytical response expression of simply supported beams, and 5% (5% of the maximum of the measured signals) Gaussian white noise has been added to the responses. The responses of the two points are shown in Figures 3 and 4.Figure 2
Real load history and identified loads ofF1.Figure 3
Displacement response of point 1 byF1.Figure 4
Displacement response of point 2 byF1.To localize the load, the first two-order modes are adopted. There is only one matrix inversion needing to be done in the identification procedure. Furthermore, the Tikhonov regularization method is used to overcome the ill-condition of the matrix inversion.The regularization parameterα is determined by the L-curve method, which is shown in Figure 5. The red dot indicates the regularization parameter α=0.0028. The two vectors S1 and S2 are calculated through equation (23). To identify the true load location, the 999 nodes are assumed as the load location, except the two ends. At each assumed load location, the two load vectors f1 and f2 can be obtained by using equation (18) and corresponding optimization function value η is computed. Figure 6 shows the curve of η and the load location with minimum value of η is 0.356m from one simply supported end, which is very close to the true load location a1=0.35m. f1 and f2 in the identified load location are shown in Figure 2. At the same time, the relative error and the correlation coefficient of f1 and f2 are listed in Table 1.Figure 5
The L-curve plot ofF1.Figure 6
The optimization function value ofF1.Table 1
The relative error and the correlation coefficient of the identifiedF1.
f
1
f
2
Pointx1
Pointx2
Relative error
0.1599
0.1845
0.231
0.1159
Correlation coefficient
0.9876
0.9876
0.9815
0.9945In Table1, the identified f1 and f2 can be used to describe the external load, but they still have some noticeable oscillatory components. To reconstruct the load history more accurately, the responses of measure points x1 and x2 are, respectively, used to identify the load history by regularizing equation (8) as the load location has already been determined. The reconstruction of the load history is shown in Figure 7. The relative error and the correlation coefficient of the reconstruction are listed in Table 1. It can be found that the reconstruction by measurement response of point 2 is closest to the real load. This is because the response value and the SNR (signal-to-noise ratio) of point x2 are larger to these of point x1. The final result is considered satisfactory.Figure 7
The reconstruction of the loadF1.


#### 3.1.2. Identification of an Impact

In this section, a simulation of the beam excited by using an impact loadF2 is carried out to verify the proposed method. The impact load F2 is shown in Figure 8, acting at the point a2=0.26m from one simply supported end. The measurement responses of the two points x1 and x2 are computed through the analytical response expression of simply supported beams, and 5% (5% of the maximum of the measured signals) Gaussian white noise has been added to responses. The responses of the two points are shown in Figures 9 and 10.Figure 8
Real load history and identified loads ofF2.Figure 9
Displacement response of point 1 byF2.Figure 10
Displacement response of point 2 byF2.In the same way, the first two-order modes are used to identify the impact location, and the Tikhonov regularization method is employed to reduce the influence of the errors on the identified result. The regularization parameterα is determined by the L-curve method, which is shown in Figure 11. The red dot indicates the regularization parameter α=0.0054.Figure 11
The L-curve plot ofF2.To localize the impact load, the optimization function valueη at each assumed load location is computed, and the curve of η is shown in Figure 12. The identified load location is 0.259m from one simply supported end. This identified result is very close to the actual impact location a2=0.26m. The corresponding f1 and f2 in the identified load location are compared in Figure 8. We also use the response of measure point x2 to reconstruct the impact history, which is shown in Figure 13. The relative error and the correlation coefficient of the reconstruction are listed in Table 2. The reconstruction can describe the external load well, so the identified result is accepted.Figure 12
The optimization function value ofF2Figure 13
The reconstruction of the impactF2.Table 2
The relative error and the correlation coefficient of the identifiedF2.
f
1
f
2
Pointx1
Pointx2
Relative error
0.2051
0.1808
0.1339
0.1138
Correlation coefficient
0.9713
0.9792
0.9867
0.9911


### 3.2. Experiment

For experimental validation of the proposed localization method, an experimental test of a steel simply supported beam was set up. As shown in Figure14, the steel rectangular beam is simply supported at the two ends and the geometric dimension of this beam is measured. The length, width, and thickness are l=0.695m, w=0.04m, and h=0.007m, respectively. The sensors used in this experiment are two PCB piezoelectric accelerometers and one force transducer of Model Number 208C02. An impact hammer with a rubber head is used to excite the beam, and the sampling frequency is 4096 Hz.Figure 14
Experimental model of the simply supported beam.The modal parameters are obtained by the system identification techniques. The experiment values of the first three-order natural frequencies are recorded in Table3. The experiment values of the first three-order modal damping ratio are ξ1=0.0026, ξ2=0.0013, and ξ3=0.0009. According to the natural frequencies obtained by the experiment, the finite element model of the beam is established and the modal shape values φ can be calculated from this finite element model. The beam is evenly divided into 695 elements, and there are 696 nodes in total. The simulation values of the first three-order natural frequencies are also recorded in Table 3.Table 3
Natural frequency of the beam (unit: Hz).
1st
2nd
3rd
Test
39.37
153.41
346.74
Simulation
38.48
153.91
346.30Eleven mark points are distributed on the beam, and two accelerometers are installed at points 3 and 6, shown in Figure14. The beam is impacted at points A, B, and C separately. The real location of these impacts is shown in Table 4. The measurement of the accelerometers is integrated twice to obtain the displacement responses of the measure points. As shown in Figure 15, an FRF analysis shows that the displacement responses are mainly involved in the first two modes, which are used to identify the impact location. After computing the optimization function value η at each node, the impact locations are identified and shown in Figures 16–18 separately.Table 4
Identified locations of the experiment.
A
B
C
True locations (m)
0.15
0.25
0.43
Identified locations (m)
0.151
0.252
0.431
Absolute error (m)
0.001
0.002
0.001Figure 15
FRF of the displacement responses by impact A.Figure 16
The optimization function value of impact A.Figure 17
The optimization function value of impact B.Figure 18
The optimization function value of impact C.The details of the identified locations are shown in Table4. It can be seen that the identified locations coincide with the true impact locations well for all three impact points A, B, and C. After determining the impact location, the force histories can be reconstructed. By using the responses of sensor 2 installed at point 6, the comparison of the measured and identified time histories of the three impacts are shown in Figures 19–21, respectively. The relative error and the correlation coefficient of the reconstructions are listed in Table 5. For all three impact test, the identified results have good consistency. Thus, the experiment result is satisfactory and verifies the credibility of the proposed method.Figure 19
The reconstruction of the impact A.Figure 20
The reconstruction of the impact B.Figure 21
The reconstruction of the impact C.Table 5
The relative error and the correlation coefficient of the identified impacts.
A
B
C
Relative error
0.1440
0.1195
0.1207
Correlation coefficient
0.9943
0.9957
0.9965


## 3.1. Numerical Simulation

In this section, two numerical simulations of a simply supported beam model are implemented to verify the proposed method. In this beam model, the length is1m long, and the cross section is 6cm×1cm. The material of the beam has the density 7800kg/m3, Young’s module 210GPa and Poisson’s ratio 0.3. The beam is evenly divided into 1000 elements, and there are 1001 nodes in total. The first two-order natural frequencies of the beam are f1=23.5Hz and f2=94.1Hz. Two displacement sensors are installed at x1=0.23m and x2=0.45m from one simply supported end, and the beam model is excited by a sine load and an impact separately, as shown in Figure 1.Figure 1
Beam simulation model.


### 3.1.1. Identification of a Sine Load

The beam is excited by a loadF1, which is shown in Figure 2. F1 is a sine load with frequency f=60Hz acting at the point a1=0.35m from one simply supported end. The measurement responses of the two points x1 and x2 are computed through the analytical response expression of simply supported beams, and 5% (5% of the maximum of the measured signals) Gaussian white noise has been added to the responses. The responses of the two points are shown in Figures 3 and 4.Figure 2
Real load history and identified loads ofF1.Figure 3
Displacement response of point 1 byF1.Figure 4
Displacement response of point 2 byF1.To localize the load, the first two-order modes are adopted. There is only one matrix inversion needing to be done in the identification procedure. Furthermore, the Tikhonov regularization method is used to overcome the ill-condition of the matrix inversion.The regularization parameterα is determined by the L-curve method, which is shown in Figure 5. The red dot indicates the regularization parameter α=0.0028. The two vectors S1 and S2 are calculated through equation (23). To identify the true load location, the 999 nodes are assumed as the load location, except the two ends. At each assumed load location, the two load vectors f1 and f2 can be obtained by using equation (18) and corresponding optimization function value η is computed. Figure 6 shows the curve of η and the load location with minimum value of η is 0.356m from one simply supported end, which is very close to the true load location a1=0.35m. f1 and f2 in the identified load location are shown in Figure 2. At the same time, the relative error and the correlation coefficient of f1 and f2 are listed in Table 1.Figure 5
The L-curve plot ofF1.Figure 6
The optimization function value ofF1.Table 1
The relative error and the correlation coefficient of the identifiedF1.
f
1
f
2
Pointx1
Pointx2
Relative error
0.1599
0.1845
0.231
0.1159
Correlation coefficient
0.9876
0.9876
0.9815
0.9945In Table1, the identified f1 and f2 can be used to describe the external load, but they still have some noticeable oscillatory components. To reconstruct the load history more accurately, the responses of measure points x1 and x2 are, respectively, used to identify the load history by regularizing equation (8) as the load location has already been determined. The reconstruction of the load history is shown in Figure 7. The relative error and the correlation coefficient of the reconstruction are listed in Table 1. It can be found that the reconstruction by measurement response of point 2 is closest to the real load. This is because the response value and the SNR (signal-to-noise ratio) of point x2 are larger to these of point x1. The final result is considered satisfactory.Figure 7
The reconstruction of the loadF1.


### 3.1.2. Identification of an Impact

In this section, a simulation of the beam excited by using an impact loadF2 is carried out to verify the proposed method. The impact load F2 is shown in Figure 8, acting at the point a2=0.26m from one simply supported end. The measurement responses of the two points x1 and x2 are computed through the analytical response expression of simply supported beams, and 5% (5% of the maximum of the measured signals) Gaussian white noise has been added to responses. The responses of the two points are shown in Figures 9 and 10.Figure 8
Real load history and identified loads ofF2.Figure 9
Displacement response of point 1 byF2.Figure 10
Displacement response of point 2 byF2.In the same way, the first two-order modes are used to identify the impact location, and the Tikhonov regularization method is employed to reduce the influence of the errors on the identified result. The regularization parameterα is determined by the L-curve method, which is shown in Figure 11. The red dot indicates the regularization parameter α=0.0054.Figure 11
The L-curve plot ofF2.To localize the impact load, the optimization function valueη at each assumed load location is computed, and the curve of η is shown in Figure 12. The identified load location is 0.259m from one simply supported end. This identified result is very close to the actual impact location a2=0.26m. The corresponding f1 and f2 in the identified load location are compared in Figure 8. We also use the response of measure point x2 to reconstruct the impact history, which is shown in Figure 13. The relative error and the correlation coefficient of the reconstruction are listed in Table 2. The reconstruction can describe the external load well, so the identified result is accepted.Figure 12
The optimization function value ofF2Figure 13
The reconstruction of the impactF2.Table 2
The relative error and the correlation coefficient of the identifiedF2.
f
1
f
2
Pointx1
Pointx2
Relative error
0.2051
0.1808
0.1339
0.1138
Correlation coefficient
0.9713
0.9792
0.9867
0.9911


## 3.1.1. Identification of a Sine Load

The beam is excited by a loadF1, which is shown in Figure 2. F1 is a sine load with frequency f=60Hz acting at the point a1=0.35m from one simply supported end. The measurement responses of the two points x1 and x2 are computed through the analytical response expression of simply supported beams, and 5% (5% of the maximum of the measured signals) Gaussian white noise has been added to the responses. The responses of the two points are shown in Figures 3 and 4.Figure 2
Real load history and identified loads ofF1.Figure 3
Displacement response of point 1 byF1.Figure 4
Displacement response of point 2 byF1.To localize the load, the first two-order modes are adopted. There is only one matrix inversion needing to be done in the identification procedure. Furthermore, the Tikhonov regularization method is used to overcome the ill-condition of the matrix inversion.The regularization parameterα is determined by the L-curve method, which is shown in Figure 5. The red dot indicates the regularization parameter α=0.0028. The two vectors S1 and S2 are calculated through equation (23). To identify the true load location, the 999 nodes are assumed as the load location, except the two ends. At each assumed load location, the two load vectors f1 and f2 can be obtained by using equation (18) and corresponding optimization function value η is computed. Figure 6 shows the curve of η and the load location with minimum value of η is 0.356m from one simply supported end, which is very close to the true load location a1=0.35m. f1 and f2 in the identified load location are shown in Figure 2. At the same time, the relative error and the correlation coefficient of f1 and f2 are listed in Table 1.Figure 5
The L-curve plot ofF1.Figure 6
The optimization function value ofF1.Table 1
The relative error and the correlation coefficient of the identifiedF1.
f
1
f
2
Pointx1
Pointx2
Relative error
0.1599
0.1845
0.231
0.1159
Correlation coefficient
0.9876
0.9876
0.9815
0.9945In Table1, the identified f1 and f2 can be used to describe the external load, but they still have some noticeable oscillatory components. To reconstruct the load history more accurately, the responses of measure points x1 and x2 are, respectively, used to identify the load history by regularizing equation (8) as the load location has already been determined. The reconstruction of the load history is shown in Figure 7. The relative error and the correlation coefficient of the reconstruction are listed in Table 1. It can be found that the reconstruction by measurement response of point 2 is closest to the real load. This is because the response value and the SNR (signal-to-noise ratio) of point x2 are larger to these of point x1. The final result is considered satisfactory.Figure 7
The reconstruction of the loadF1.


## 3.1.2. Identification of an Impact

In this section, a simulation of the beam excited by using an impact loadF2 is carried out to verify the proposed method. The impact load F2 is shown in Figure 8, acting at the point a2=0.26m from one simply supported end. The measurement responses of the two points x1 and x2 are computed through the analytical response expression of simply supported beams, and 5% (5% of the maximum of the measured signals) Gaussian white noise has been added to responses. The responses of the two points are shown in Figures 9 and 10.Figure 8
Real load history and identified loads ofF2.Figure 9
Displacement response of point 1 byF2.Figure 10
Displacement response of point 2 byF2.In the same way, the first two-order modes are used to identify the impact location, and the Tikhonov regularization method is employed to reduce the influence of the errors on the identified result. The regularization parameterα is determined by the L-curve method, which is shown in Figure 11. The red dot indicates the regularization parameter α=0.0054.Figure 11
The L-curve plot ofF2.To localize the impact load, the optimization function valueη at each assumed load location is computed, and the curve of η is shown in Figure 12. The identified load location is 0.259m from one simply supported end. This identified result is very close to the actual impact location a2=0.26m. The corresponding f1 and f2 in the identified load location are compared in Figure 8. We also use the response of measure point x2 to reconstruct the impact history, which is shown in Figure 13. The relative error and the correlation coefficient of the reconstruction are listed in Table 2. The reconstruction can describe the external load well, so the identified result is accepted.Figure 12
The optimization function value ofF2Figure 13
The reconstruction of the impactF2.Table 2
The relative error and the correlation coefficient of the identifiedF2.
f
1
f
2
Pointx1
Pointx2
Relative error
0.2051
0.1808
0.1339
0.1138
Correlation coefficient
0.9713
0.9792
0.9867
0.9911


## 3.2. Experiment

For experimental validation of the proposed localization method, an experimental test of a steel simply supported beam was set up. As shown in Figure14, the steel rectangular beam is simply supported at the two ends and the geometric dimension of this beam is measured. The length, width, and thickness are l=0.695m, w=0.04m, and h=0.007m, respectively. The sensors used in this experiment are two PCB piezoelectric accelerometers and one force transducer of Model Number 208C02. An impact hammer with a rubber head is used to excite the beam, and the sampling frequency is 4096 Hz.Figure 14
Experimental model of the simply supported beam.The modal parameters are obtained by the system identification techniques. The experiment values of the first three-order natural frequencies are recorded in Table3. The experiment values of the first three-order modal damping ratio are ξ1=0.0026, ξ2=0.0013, and ξ3=0.0009. According to the natural frequencies obtained by the experiment, the finite element model of the beam is established and the modal shape values φ can be calculated from this finite element model. The beam is evenly divided into 695 elements, and there are 696 nodes in total. The simulation values of the first three-order natural frequencies are also recorded in Table 3.Table 3
Natural frequency of the beam (unit: Hz).
1st
2nd
3rd
Test
39.37
153.41
346.74
Simulation
38.48
153.91
346.30Eleven mark points are distributed on the beam, and two accelerometers are installed at points 3 and 6, shown in Figure14. The beam is impacted at points A, B, and C separately. The real location of these impacts is shown in Table 4. The measurement of the accelerometers is integrated twice to obtain the displacement responses of the measure points. As shown in Figure 15, an FRF analysis shows that the displacement responses are mainly involved in the first two modes, which are used to identify the impact location. After computing the optimization function value η at each node, the impact locations are identified and shown in Figures 16–18 separately.Table 4
Identified locations of the experiment.
A
B
C
True locations (m)
0.15
0.25
0.43
Identified locations (m)
0.151
0.252
0.431
Absolute error (m)
0.001
0.002
0.001Figure 15
FRF of the displacement responses by impact A.Figure 16
The optimization function value of impact A.Figure 17
The optimization function value of impact B.Figure 18
The optimization function value of impact C.The details of the identified locations are shown in Table4. It can be seen that the identified locations coincide with the true impact locations well for all three impact points A, B, and C. After determining the impact location, the force histories can be reconstructed. By using the responses of sensor 2 installed at point 6, the comparison of the measured and identified time histories of the three impacts are shown in Figures 19–21, respectively. The relative error and the correlation coefficient of the reconstructions are listed in Table 5. For all three impact test, the identified results have good consistency. Thus, the experiment result is satisfactory and verifies the credibility of the proposed method.Figure 19
The reconstruction of the impact A.Figure 20
The reconstruction of the impact B.Figure 21
The reconstruction of the impact C.Table 5
The relative error and the correlation coefficient of the identified impacts.
A
B
C
Relative error
0.1440
0.1195
0.1207
Correlation coefficient
0.9943
0.9957
0.9965


## 4. Conclusion

In this paper, the problem of load localization and reconstruction from structural responses is addressed. In order to improve the problem that a lot of matrix inversions in the identification process consume much operation time, a variable separation method is proposed. This method separates the load location variable from the impulse response matrix which needs to be inversed in the localization process. We take several order modes in which the response is mainly included and identify the stable time history vectors of the corresponding modes. In this step, there is only once matrix inversion and it can be replaced by singular value decomposition method. The regularization method is employed to overcome the ill-posed problem and obtain the stable time history vector. Then we compute the force history through dividing the time history vector by corresponding mode shape value. The load location is determined with error function using the minimum optimization method. After identifying the location, this problem is transformed into the classic reconstruction of force history. By using the regularization method and L-curve criterion, the force history is inversely computed. The proposed method is fully demonstrated and verified with simulations of a simply supported beam separately acted by a sine load and an impact. An experiment is also implemented to prove the validity of the method.



---

*Source: 4207473-2019-04-01.xml*