Regularized Cubic B-Spline Collocation Method With Modified L-Curve Criterion for Impact Force Identification

The time history of the impact force, especially for the peak force, is vital to monitor the performance of mechanical products over the life-time. However, considering the limitation of sensing technology and inaccessibility of installing, it is always difficult or even impossible to measure the force directly in engineering practice. Therefore, a regularized cubic B-spline collocation (RCBSC) method combined with the modified L-curve criterion (RCBSC-ML) is presented to identify the impact force by easily measurable dynamic response. Because the profile of cubic B-spline is close to that of impact force, a linear combination of cubic B-spline functions is used to approximate the unknown impact force. Then the force identification equation is reformed, wherein the unknown variable changes from impact force vector to weight coefficient vector. Furthermore, the modified L-curve criterion is developed to select the optimal number of collocation points which is the regularization parameter of RCBSC method. Rich simulation studies on a five degree-of-freedoms structure and experimental studies on a simplified artillery test-bed are performed to validate the performance of RCBSC-ML method. The results illustrate that RCBSC-ML method can be used to accurately identify the unknown impact force, and compared with modified generalized cross validation criterion, the presented modified L-curve criterion is more suitable to determine the optimal regularization parameter of RCBSC method. Furthermore, compared with classical Tikhonov regularization method, RCBSC-ML method is more efficient and accurate in impact force identification.


I. INTRODUCTION
Impact events are one of the main factors that aggravate the vibration of mechanical systems, and the time history of impact force (especially for the peak force) is the vital index for structural health monitoring, reliability analysis and optimal design [1]- [4]. There are two main ways to obtain the force acted on the mechanical structure, namely direct measurement and indirect identification. The former means that a suitable sensor is installed at the acting position of the force to directly measure the force value. However, considering the limitation of sensing technology and inaccessibility of installing, it is always difficult or even impossible to measure The associate editor coordinating the review of this manuscript and approving it for publication was Baoping Cai . dynamic force directly in engineering practice. Therefore, in the past few decades, many scholars are committed to develop the force identification technique to indirectly obtain the time history or peak of impact force [5]- [7].
Based on the dynamic response (such as acceleration, velocity, displacement, or strain) of mechanical systems which can be easily measured, the appropriate force identification method is applied to solve force identification equation to obtain the unknown force. Unfortunately, the transfer matrix in the force identification equation is usually an illconditioned matrix [8]. The inverse of ill-conditioned transfer matrix can amplify the small error (unavoidable noise interference) existing in the measured response, resulting in the solution seriously deviating from the true value of the force. It means that the force identification problem of mechanical VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ structures does not satisfy the Hadmard stability condition, and it is a typical ill-posed inverse problem [9], [10]. Regularization method is the most common technique for solving such ill-posed inverse problem [11]- [13]. A set of neighboring well-posed questions is used to approximate the original ill-posed problem, and the solution of the wellposed question can be viewed as the approximate solution of the original problem. Leclere et al. [14] used the truncated singular value decomposition (TSVD) method to identify the bearing capacity acting on the main shaft of the diesel engine. Chen and Chan [15] developed TSVD method to identify the moving force, and considered the effect of sensor combinations and sensor numbers on the identification accuracy. Based on a time-domain deconvolution model, Jacquelin et al. [16] compared the performance of TSVD method and Tikhonov regularization method in identifying the impact force of aluminium plate, and the result shows that global identification accuracy by Tikhonov regularization method is better than that by TSVD method. Choi et al [17] studied the performance of OCV (ordinary cross validation) criterion, GCV (generalized cross validation) criterion, and L-curve criterion in determining Tikhonov regularization parameters through force identification of a four-sided simply supported plate. It shows that, when the noise level is high, the performance of L-curve criterion is better than that of OCV criterion and GCV criterion. When the noise level is low, the above conclusion is reversed. Trabelsi et al. [18] combined Tikhonov regularization method and L-curve criterion to identify the aerodynamic force acting on the axial fan blade model. For the four possible excitation positions on a composite board, Kalhori et al. [19] used Tikhonov regularization method to identify the force and its position based on an underdetermined equation. The above methods have excellent performance and have been widely used in force identification in engineering practice. However, these methods include the complex operation such as matrix factorization, and it would reduce computational efficiency. Furthermore, for the impact force, the identification performance (especially for peak accuracy) of these above method needs to be improved [20].
The force identification equation has the form of the first-order Fredholm integral equation [21]. Therefore, basis function collocation method (also called basis function expansion method) that has been widely used for solving firstorder Fredholm integral equation can also be applied to solve the force identification equation, and it can be viewed as a special kind of regularization method. For the basis function collocation method, various orthogonal basis functions are utilized to approximate the force identification model (force vector and response vector). And then, the unknown variable in the force identification model changes from unknown force vector to weight coefficient vector. It could help to reduce the time consuming and improve the ill-posedness of inverse matrix. The commonly used basis functions are polynomials, power series, wavelet functions, and spline functions [22]- [24]. Hu et al. [25] proposed a Chebyshev polynomial collocation method to identify the impact force acting on CFRP (Carbon Fibre Reinforced Plastics) laminated plate, and they also pointed out that, in order to approximate the impact force, the number of Chebyshev polynomial can be set from 20 to 40. Liu et al. [26] combined Gegenbauer polynomial and Tikhonov regularization method to identify dynamic force acting on the stochastic structure. Li and Deng [27] proposed the second-order Taylor-series expansion method to identify the dynamic force in time domain, and the simulation study on a cantilever model shows that the proposed method has high identification accuracy. Li et al. [28] developed a comprehensive algorithm combining Taylor polynomial iteration and cubic Catmull-Rom spline interpolation to identify the distributed dynamic excitation in the domain. Li et al. [29] proposed a novel function collocation method based on wavelet multi-resolution analysis. Through wavelet decomposition, the convolution relation between the exciting force and dynamic response is transformed into a linear multiplicative relation in the wavelet domain. Based on representing the unknown force by the expanded coefficients at different scales and shifts, Tran and Inoue [30] proposed a wavelet deconvolution technique to identify the impact force acting on the thin-walled column.
When the basis function collocation method is applied to solve the force identification problem, the identification accuracy depends on the approximation degree of the basis function to the force identification model and the improvement degree of the ill-posedness of inverse matrix. Type and number of basis function are two key elements of the basis function collocation method. It means that selected basis function can well match the profile of the force to be identified, and we also should determine the optimal number of basis functions (like determining optimal Tikhonov regularization parameter) to improve the ill-posedness of force identification problem. Cubic B-spline function has the characteristic of local smoothness, and its profile can well match impact force studied in this paper. Based on an alternative cubic B-spline formula and QR decomposition, Gunawan et al. [31] presented a two-step B-splines regularization method to identify impact force. In order to improve the computational efficiency of the identification method presented in [31], Qiao et al. [20] combined cubic B-spline collocation (CBSC) method and the modified GCV criterion to conduct impact force identification. The identification method does not include time-consuming operations such as QR decomposition and SVD, and the modified GCV criterion is utilized to determine a suitable number of cubic B-spline functions. However, as shows in [20], the identification accuracy (especially peak identification accuracy) still needs to be improved when the frequency range of impact force is large, and the value of modified GCV function tends to be stable when the number of cubic B-spline functions is greater than a certain value. It is difficult to determine a reliable minimum of modified GCV function. Furthermore, transfer matrix in the force identification equation is a Toeplitz matrix composed of unit impulse response function, which is obtained by inverse Fourier transform (IFT) of frequency response function (FRF). On the one hand, IFT would introduce numerical truncating errors. On the other hand, once the dimension of the force vector to be identified increases, a new vibration testing needs to be conducted to obtain enough data to form the corresponding transfer matrix.
In this paper, a novel method combined regularized cubic B-spline collocation (RCBSC) method and a modified L-curve criterion (RCBSC-ML) is proposed to solve ill-posed inverse problem to obtain the unknown impact force. For the proposed identification method, force identification equation is built on the base of the state-space model of mechanical structures, in which the transfer matrix is formed by coefficient matrices of the state-space model. Then, cubic B-spline function is utilized to approximate the profile of impact force to form a new transfer matrix and weight coefficient vector. By the above processing, solving the unknown force vector is transformed into solving the unknown weight coefficient vector, which can greatly improve the ill-posedness of the inverse problem. Afterwards, the modified L-curve criterion is applied to determine the optimal number of cubic B-spline functions or number of collocation points (these two parameters can be converted to each other), and the weight coefficient vector can be obtained by the simple least squares algorithm. Finally, the unknown impact force can be identified by the weight coefficient vector and cubic B-spline function.
The remaining part of this paper is organized as follows. In section II, force identification equation on the base of state-space model is built. In section III, the regularized cubic B-spline collocation method is introduced to identify the impact force, and a modified L-curve criterion is developed to determine the optimal number of cubic B-spline functions. In section IV, numerical study on a five degree-of-freedom chain type structure is applied to verify the feasibility of RCBSC-ML method, and some influence factors are also discussed to show the robustness of the identification method. In section V, experimental data from a simplified artillery test-bed are utilized to verify the effectiveness of RCBSC-ML method. Finally, several conclusions are given in section VI.

II. FORCE IDENTIFICATION EQUATION
For mechanical structures, the output response y(t) can be described by the convolution of the excitation force f (t) and the unit impulse response function h(t), and the expression is given by where the symbol ⊗ denotes the convolution operation. τ is the time-shift factor satisfying t ≥ τ . The output response y(t) is not only the displacement response, but also the acceleration response, velocity response, strain, and other measured physical quantities. It can be noted from (1) that, when the output response y(t) and the unit impulse response function h(t) are known, the unknown excitation force f (t) can be obtained by deconvolution operation of the above convolution formula. It means that a force identification equation can be formed by the deconvolution of (1) [20].For this force identification equation, the transfer matrix is composed of the unit impulse response function h(t). In engineering practice, h(t) is usually obtained by inverse Fourier transform (IFT) of frequency response function (FRF). The measurement error of FRF at high-frequency part is large and discrete Fourier transform (DFT) can introduce numerical truncating errors. Additionally, in the force identification equation, the dimension of unknown force vector should correspond to that of the transfer matrix, and this dimension is also related to analysis time. It means that, when the dimension of the force vector to be identified changes, a new vibration testing needs to be conducted to obtain the corresponding transfer matrix.
In addition to the deconvolution operation of (1), the force identification equation can also be obtained from the statespace model of mechanical structures. Based on defining a state variable z(t) = x(t)ẋ(t) T , the differential equation of motion can be converted to the continuous-time state-space model [32] where subscript · c corresponds to the continuous-time model. A c , B c , C c , and D c are the coefficient matrices of the continuous-time state-space model. n denotes the order of the model. C α , K, and M are the damping, stiffness, and mass matrices, respectively. C ac , C ve , and C di are output influence matrices related to acceleration, velocity, and displacement, respectively. It means that, for different response types, the values of these output influence matrices are different. For instance, when the output response is the acceleration response, these influence matrices are C ac = I (I denotes unit matrix) and C ve = C di = 0, respectively. Additionally, the corresponding coefficient matrices in Considering that the measured input and output data are discrete, the continuous-time model given in (2) and (3) should be converted to the corresponding discrete-time model where subscript · d corresponds to the discrete-time model. t denotes sampling interval, and k t is the non-negative integer. The Zero-Order-Hold (ZOH) discretization method [33] is adopted in this paper. It means that the input signal of discrete-time model is regarded as constant within the sample interval. Based on ZOH discretization, these coefficient matrices of the above two models follow the relationship By iterating (4) and (5), a linear regression equation can be obtained as The mechanical system is the causal system. The system generates the output response only when the input signal acts on the mechanical system. It means that there is the zero initial condition for causal system (x(0) = 0 andẋ(0) = 0). Based on zero initial condition and (7), the output response at different sampling times can be expressed as where N denotes the number of sampling points. Equation (8) can also rewritten as a compact matrix-vector form It can be noted from (9) that input vector (excitation force vector) f can be directly related to output vector (dynamic response vector) Y by a Toeplitz matrix (transfer matrix) H. When Y and H are known, unknown force f can be obtained by solving (9). Therefore, equation (9) can be viewed as the force identification equation. For the output response Y, it can be measured by the suitable sensor. For the transfer matrix H, it is composed of coefficient matrices of discrete-time model, and these coefficient matrices can be calculated by subspace identification (SI) method [34]. SI method is developed from the cybernetics and has been widely used in mechanical structure analysis, and it is a classic parameter identification method of linear structure. Once input and output vectors of the state space model are formed by the measured dynamic data of structural vibration test, SI method can be applied to calculate coefficient matrices of the state space model, wherein it contains matrix projection technique. SI method has the characteristics of high computational efficiency and strong robustness, and the whole calculation process is carried out in the time domain. It means that DFT is not performed in the calculation process. Furthermore, it can be observed from (2) and (3) that these coefficient matrices are formed by system characteristic matrices (M, K, and C α ), and these characteristic matrices are invariant for the given structure. It means that these coefficient matrices are invariant for different excitation cases. When the dimension of f changes, we just use original coefficient matrices to rebuild a corresponding transfer matrix as the same law shown in (8), rather than conducting a new vibration test. Therefore, force identification equation shown in (8) has superior performance and is adopted in this paper to conduct force identification.

III. APPLICATION OF CUBIC B-SPLINE FUNCTION FOR IDENTIFYING IMPACT FORCE
It can be seen from (9) that inverse operation of the transfer matrix is one key step for identifying unknown force. Considering that the condition number of transfer matrix is always large in mechanical structures, the small error existing in the measured response can be greatly amplified and result in the worthless solution. It means that force identification is a typical ill-posed inverse problem. Tikhonov method and TSVD method which are well-known regularization methods can be used to obtain the stable solution of such ill-posed inverse problem, and have been widely utilized to conduct force identification in engineering practice. However, SVD technique contained in the calculation process of these methods would reduce computational efficiency, and the identification performance for impact force needs to be improved [6]. Different from the above two regularization methods, basis function collocation method which is a special kind of regularization method does not involve matrix factorization. By adjusting the number of basis functions, the ill-posed force identification problem in the time domain or frequency domain is transformed into a well-posed problem in basis function space. The selected basis function should well match the profile of the force to be identified. More importantly, the optimal number of basis functions should be determined because it acts as a regularization parameter. Cubic B-spline function has the characteristic of local smoothness, and its profile can well match impact force. Therefore, in this paper, cubic B-spline collocation method is utilized to address the impact force identification, and modified L-curve criterion is proposed to select optimal number of cubic B-spline functions.

A. CUBIC B-SPLINE COLLOCATION METHOD
B-spline is one kind of spline functions originally proposed by American mathematician Schoenberg in the mid-1940s. It has been widely used in numerical analysis, difference problem and optimal control [35], [36]. Cubic B-spline is the piecewise cubic polynomial with second order continuous at piecewise points.
For the cubic B-spline collocation method, cubic B-spline is as the basis function and its linear combination is utilized to approximate unknown impact force vector f. And then, force identification problem is transformed into solving the unknown weight coefficient vector of the linear combination, which can greatly improve the ill-posedness of the inverse problem.
In order to apply the cubic B-spline collocation method for identifying impact force, the time range . These knots which are also called collocation points of basis function is given by where m s is the number of collocation points. Based on these collocation points, the impact force can be expressed in term of cubic B-spline function where α i is the weight coefficient corresponding to cubic B-spline function S ci (t), and S ci (t) is defined as It can be seen from (12) that a complete cubic B-spline function S ci (t) contains five collocation points  (11). It also means that m b (number of cubic B-spline functions) and m s (number of collocation points) can be converted to each other. Furthermore, S ci (t) is the complete cubic B-spline function, and its function expression is shown in (12). However, for i ≤ 2 and i ≥ m s − 1, some piecewise polynomials in (12) are replaced by 0, and the corresponding functional expressions are given by From (12) to (18), incomplete cubic B-spline S ci (t) (i = 0, 1, 2, m s − 1, m s , m s + 1) curve is the part of the complete cubic B-spline S ci (t) (3 ≤ i ≤ m s − 2) curve, and the profiles of these curves are identical. It means that, as shown in (19), these m s + 2 cubic B-spline functions can be obtained by translating a complete cubic B-spline function over the collocation points.
where i 2 ∈ [3, m s −2] means that S ci 2 (t) denotes the complete cubic B-spline calculated by (12). For instance, S c2 (t) can be obtained  (9) can be rewritten as Once the number of collocation points m s is determined, matrix can be calculated by (12) to (18). Therefore, for the cubic B-spline collocation method, original transfer matrix H in force identification equation is replaced by a new transfer matrix = H , and the weight coefficient vector α is the new unknown variable. It means that force identification problem is transformed into solving the unknown α. As discussed above, the condition number of original transfer matrix H is always large in mechanical structures. However, matrix composed of cubic B-spline functions is a banded sparse matrix, and the new transfer matrix is a diagonally dominant matrix. Therefore, the condition number of is far less than that of H. It can result in the improvement of illposeness of solving (12), and the unknown α can be obtained directly by least squares algorithm where symbol · † denotes pseudo inverse operation of matrix. After obtaining weight coefficient vector α, impact force to be identified can be estimated by

B. REGULARIZATION AND MODIFIED L-CURVE CRITERION
Based on (12), the cubic B-spline curves with different m s (number of collocation points) are shown in Fig. 1. The figure shows that the profile of the cubic B-spline curve is similar to that of impact force. Thus, it is reasonable to use cubic B-spline to approximate impact force. It can also be seen from Fig. 1 that the shape of cubic B-spline is related to m s . As m s increases, range of the impact loading area decreases, and the corresponding frequency band of impact force would increases. It means that the approximate degree between cubic spline function and impact force can be adjusted by changing m s . Furthermore, m s can also adjust the reduction degree of cond(H) (condition number of transfer matrix H) by matrix ∈ R N ×(m s +2) . Therefore, the optimal m s should be selected by an appropriate criterion to obtain cubic B-spline matrix . And then, unknown impact force can be calculated by (21) and (22). The above identification process is called the regularized cubic B-spline collocation (RCBSC) method. Determining optimal parameter m s plays a key role of successfully solving such ill-posed problem. It is called regularization parameter of the cubic B-spline collocation method like Tikhonov regularization parameter and truncation threshold of TSVD method.
Currently, some posterior criterions are usually applied to select the regularization parameter. The optimal parameter is determined by extracting useful information from the measured output response as much as possible without knowing the noise level. GCV criterion and L-curve criterion which are common posterior criterions have been widely used to determine the regularization parameters of Tikhonov regularization and TSVD methods [12], [17], [37]. Additionally, scholars have carried out a lot of research to improve the applicability and robustness of the above two selection criterions. Based on a so-called U function developed form classical L-curve criterion, Krawczyk-StańDo and Rudnicki proposed a U-curve criterion to determine the regularization parameter of Tikhonov regularization method [38]. Bazán and Francisco introduced an improved fixed-point algorithm to determine the regularization parameter even in the case where the L-curve has several local corners [39]. Pazos and Bhaya proposed an adaptive selection method via Liapunov optimizing control to determine the Tikhonov regularization parameter [40]. However, RCBSC method discussed in this paper is a special regularization method. The corresponding regularization parameter is the number of cubic B-spline basis functions, rather than an additional variable. Therefore, the previous selection methods cannot be used directly to determine it. Qiao et al. developed a modified GCV criterion to determine the number of collocation points for the cubic B-spline collocation method [20], and the modified GCV function is expressed by For the modified GCV criterion, parameter m s corresponding to the minimum value of the modified GCV function is the optimal regularization parameter. This criterion can solve the problem of determining m s to a certain extent, so as to realize the impact force identification by RCBSC method. However, as shown in [20], for the case that the impact loading area is small (the corresponding frequency band of impact force is large), the peak error by RCBSC method combined with modified GCV criterion is large. As discussed above, when range of the impact loading area decreases, optimal m s required for RCBSC method would increase. From (23), in order to achieve minimum value of GCV cubic (m s ), residual norm ||Y − Hf|| 2 and m s need to be as small as possible. It means that, for RCBSC method, the modified GCV criterion directly limits the number of collocation points. Therefore, the performance of modified GCV criterion needs to be improved to obtain accurate identification result of impact force by RCBSC method. Furthermore, Choi et al. illustrates that the performance of L-curve criterion is better than that of GCV criterion when noise level is high [17]. Through the above discussion, L-curve criterion is developed to select the optimal m s for RCBSC method in this paper.
For original L-curve criterion, when it is used to select the optimal Tikhonov regularization parameter or TSVD truncation threshold, a series of possible regularization parameter are drawn in a logarithmic coordinate system. Solution norm ||f|| 2 and residual norm ||Y − Hf|| 2 are respectively vertical coordinate and horizontal coordinate. However, for RCBSC method studied in this paper, solution norm is the weight coefficient norm ||α|| 2 , and residual norm is ||Y − α|| 2 . Therefore, a series of possible m s can be drawn in a new logarithmic coordinate system, wherein ||α|| 2 and ||Y − α|| 2 are the vertical coordinate and horizontal coordinate, respectively. The shape of the above curve used to describe parameter m s is similar to the letter 'L', so this curve is called L curve. The position of the largest curvature at the corner of L curve means the corresponding solution can make a compromise between the solution norm and residual norm, and the corresponding m s is viewed as the optimal regularization parameter of RCBSC method. For the modified L-curve criterion, the curvature can be expressed as in which It is worth noting that, for Tikhonov regularization method, 'L curve' is a continuous curve with L shape. However, for RCBSC method, 'L curve' consists of a series of discrete points. Furthermore, for classic TSVD method, when determining the transfer matrix H of different truncation thresholds k T , it is only necessary to increase the column vectors of the transfer matrix corresponding to k T . It means h k T1 i = h k T2 i , and h k T1 i denotes i-th column of H when truncation threshold is k T1 . It also means that the variety in the truncation threshold only changes the number of columns (dimensions) of matrix H and does not change the elements of matrix H. However, for RCBSC method, the cubic B-spline functions calculated by (12) to (18)  when the number of collocation points is m s1 . It also means that, once m s changes, all elements in matrix should be recalculated by (12) to (18), rather than simply increasing or decreasing the column vectors of matrix .
In section I, the number of basis functions is regarded as regularization parameter so as to expediently introduce the regularization effect of basis function collocation method. However, in section III, based on the relationship between cubic B-spline function and the number of collocation points m s , m s is viewed as regularization parameter. In fact, as shown in (11), the number of cubic B-spline functions m b is equal to m s +2. It means that m b and m s can be converted to each other. Therefore, the number of cubic B-spline functions can also be viewed as the regularization parameter of RCBSC method, and the modified L-curve criterion can also be used to select the optimal m b . In the following simulation and experimental studies, m s is uniformly set as the regularization parameter so as to be consistent with the theoretical analysis in section III.

IV. NUMERICAL SIMULATION VERIFICATION
In this section, a five degree-of-freedoms (DOFs) chain type structure is established in MATLAB so as to verify the effectiveness of RCBSC method for impact force identification, and the advantage of modified L-curve criterion on determining optimal RCBSC regularization parameter m s is discussed by comparing with modified GCV criterion. The numerical studies include single-time impact identification and consecutive impact identification. Furthermore, some influence factors, including the noise level, impact position, and response measurement position, are fully discussed to further illustrate the performance of RCBSC method.
In order to quantitatively evaluate the performance of RCBSC method combined with modified L-curve criterion (RCBSC-ML), total relative error ε t and peak relative error ε p are defined as where f me denotes true force, it is set in the simulation or measured in the experiment. f id denotes identified force. max(·) represents the maximum of a vector, and | · | represents the absolute value of an element. In engineering practice, the measured response would inevitably be polluted by the interference information. Therefore, in order to simulate the actual situation, simulated output response y is added by an additive noise vector e n , and the contaminated response vector y e is defined as y e = y + e n = y + L noise · std(y) · η n (28) where L noise is the noise level, it indicates the intensity of noise interference. std(y) represents the standard deviation of simulated output response y. η n ∈ R N ×1 is a vector consist of the pseudorandom number, and it obeys uniform distribution of interval (−1, 1).

A. NUMERICAL SIMULATION MODEL
The numerical simulation model which is a five DOFs chain type structure is shown in Fig. 2. In the figure, m 1 , k 1 , c 1 , and x 1 (t) are the mass, stiffness, damping, and displacement response, respectively. The other structural parameters are defined as a similar way. In the simulation, these structural VOLUME 8, 2020 parameters of the numerical model are set as In the simulation, excitation force f (t) acting on the mass m 3 is set as impact force, and it is simulated by the following formula where F max denotes the peak force or the large amplitude of impact force. f 0 is applied to adjust the frequency band of impact force. t m indicates the time corresponding to reach the peak force. In the simulation, sampling frequency is 2000Hz, and acquisition duration is 2s. Based on the equation of motion, fourth-order Runge-Kutta Method is utilized to calculate the dynamic output response of the simulation model in MATLAB.

B. SIMULATION RESULTS AND DISCUSSIONS 1) SINGLE-TIME IMPACT CASE
Single-time impact means that the impact force acted on the structure contains only one peak. Symbol (Fi1, Ri2) is defined in this paper, and it denotes that impact position and response measurement position are point i1 and point i2, respectively. In this section, (F3, R2) case is used to illustrate the performance of RCBSC method combined with modified L-curve criterion on impact force identification. It means that the acceleration response measured at mass m 2 is applied to identify the impact force acted on mass m 3 . Noise level is set as L noise = 20%. Parameters of simulated impact force shown in (30) are set as F max = 50N, f 0 = 500Hz, and t m = 0.5s. As discussed in section II, in order to identify the impact force at mass m 3 , a vibration test should be conducted at first, and SI identification is utilized to obtain coefficient matrices of the state space model to form the transfer matrix H 23 between mass m 2 and mass m 3 . The condition number of transfer matrix H 23 is cond(H 23 ) = 1.13 × 10 21 . It illustrates that transfer matrix H 23 is an ill-conditioned matrix, and  force identification problem of the simulation model shown in Fig. 2 is ill-posed. Therefore, the regularization method should be applied to improve the ill-posedness of force identification to obtain the stable solution.
In order to identify the impact force by RCBSC method, determining optimal number of collocation points m s (or number of cubic B-spline functions m b ) is one key step. The total relative error ε t and peak relative error ε p of identified impact force for different m s are shown in Fig. 3. It shows that, as discussed in section III, RCBSC method is semiconvergent for force identification. Excessive basis functions may increase the total relative error and peak relative error. Therefore, the number of cubic B-spline functions (or number of collocation points) is regarded as the regularization parameter, and its optimal value should be determined by a reasonable criterion. Furthermore, Fig. 3 also illustrates that, for impact force identification, parameter m s corresponding to the minimum total relative error ε t is not the same as that corresponding to the minimum peak relative error ε p . It means that the identification result of impact force usually cannot be optimal for both ε t and ε p . In vibration analysis of mechanical structure, compared with variety of impact force over entire lifetime, engineers pay more attention to its peak force only in impact loading area. The peak of impact force is a key index of structural health monitoring (SHM). Therefore, compared with total relative error, scholars prefer to improve the identification performance in the peak relative error.
Regularization parameters of RCBSC method determined by modified GCV criterion and modified L-curve criterion are shown in Fig. 4. The figure shows that, based on modified GCV criterion and modified L-curve criterion, optimal regularization parameters m s are 83 and 115, respectively. The corresponding identification results by RCBSC method are 36344 VOLUME 8, 2020 shown in Fig. 5. In the figure, the measured impact force (red line) means its true value set in the simulation calculation. RCBSC-MGCV (black dash-dotted line) means the result obtained by RCBSC method combined with modified GCV criterion. RCBSC-ML (blue dotted line) means the result obtained by RCBSC method combined with modified L-curve criterion. For the identification result corresponding to m s = 83 determined by modified GCV criterion, the corresponding total relative error is 14.47%, the identified peak force is 45.77N, and peak relative error is 8.46%. For the identification result corresponding to m s = 115 determined by modified L-curve criterion, the corresponding total relative error is 20.64%, the identified peak force is 50.46N, and peak relative error is 0.92%. Through the above information, the peak relative error by RCBSC-ML method is smaller than that by RCBSC-MGCV method.
As shown in Fig. 3, parameter m s corresponding to minimum ε p is larger than that corresponding to minimum ε t . It means that enough collocation points are needed to make ε p small. As discussed in section III, for the modified GCV criterion, parameter m s is determined by the minimum value of modified GCV function shown in (23), and this function would directly limit m s so that m s cannot increase to its optimal value corresponding to minimum peak relative error. However, for the modified L-curve criterion, parameter m s is determined by the relationship between residual norm and solution norm. This criterion cannot directly limit m s . Furthermore, it can be seen from Figs. 3 and 4 that regularization parameter (83) determined by modified GCV criterion is close to that (77) determined by minimum total relative error, and regularization parameter (115) determined by modified L-curve criterion is close to that (122) determined by minimum peak relative error. Therefore, compared with modified GCV criterion, modified L-curve criterion presented in this paper is more suitable to determine the optimal regularization parameter of RCBSC method.
For RCBSC method, after determining the number of collocation points, force identification equation shown in (20) is solved by the least squares algorithm (pseudo inverse operation). The least squares algorithm is essentially an l 2 -norm regularization method. The solution is calculated by energy minimization principle, and it is always nonzero throughout time history [6]. Therefore, as shown in Fig. 5, the identification result in the non-loading area fluctuates around zero  rather than being close to zero. Furthermore, as the number of collocation number m s increases, the dimension of weight coefficient vector α to be solved increases, and the number of non-zero elements in vector α would also increase so as to aggravate the fluctuations in the non-loading area. This is why the total relative error by RBCSC-ML method is greater than that by RBCSC-MGCV method. However, as shown in Fig. 5, for the identification result by RBCSC-ML method, it is in good agreement with the measured result in loading area, and the amplitude of the fluctuation in non-loading area is not large. It illustrates that RBCSC-ML method presented in this paper is suitable for identifying the impact force.
To show the robustness of RBCSC-ML method, identification results for different noise levels are shown in Table 1. The identification operation is still for case (F3, R2), and simulation parameters are unchanged except noise level. It can be observed from Table 1 that total relative errors increase with the noise levels. This is because the non-zero characteristic of the l 2 -norm solution can amplify the noise interference existing in the measured response. However, these errors mainly exist in the non-loading area, and it has less influence on the peak accuracy of impact force. As shown in the table, for different noise levels, the range of peak relative error is 0.66%-1.05%. Even if the noise level is 40%, the peak relative error is just 1.05%. It illustrates that RCBSC-ML method has strong robustness in impact force identification.
In order to further illustrate performance of RCBSC-ML method, identification results for different excitation positions and measurement positions are shown in Table 2. Nose level is L noise = 25%, and parameter f 0 of impact force is set as 400Hz. It can be seen from Table 2 that, for different excitation positions and measurement positions, the range of total relative error is 12.19%-18.42%, and the range of peak relative error is 0.84%-1.63%. Therefore, RCBSC-ML method presented in this paper is still effective for different excitation positions and measurement positions.

2) CONSECUTIVE IMPACT CASE
The effectiveness of RCBSC method combined with modified L-curve criterion (RCBSC-ML method) for single-time impact case has been verified in the previous section. However, in engineering practice, the mechanical structure is usually impacted continuously by the excitation source. For instance, during the service of a wind turbine, its blades are subject to consecutive impact from wind loads. Additionally, during the continuous firing of the artillery, the barrel is also subjected to the consecutive impact from bore pressure. These peak forces of consecutive impact case should be accurately obtained to monitor the operation state of mechanical structures. Therefore, the identification performance of RCBSC-ML method for consecutive impact case is studied in this section.
Case (F4-2, R2) means the acceleration response measured at mass m 2 is applied to identify the consecutive two-time impact acted on mass m 4 . Based on (30), the consecutive twotime impact f 42 (t) is simulated by The acquisition duration is 2s, and noise level is L noise = 20%. The structural parameters of the simulation model and other simulation calculation parameters are consistent with the single-time impact case. The identification process of consecutive impact case is same as that of single impact case. For case (F4-2, R2), the transfer matrix H 24-2 between mass m 2 and m 4 should be obtained by (8) and SI identification method at first. The condition number of transfer matrix H 24-2 is cond(H 24-2 ) = 1.91 × 10 21 . It shows that consecutive impact identification problem is also ill-posed. When RCBSC-ML method is utilized to conduct the consecutive impact identification, modified L-curve criterion for choosing number of collocation points is shown in Fig. 6(a). From the figure, determined optimal parameter m s is 192, and the corresponding condition number of new transfer matrix = H 24-2 is cond( ) = 235.67 which is much smaller than cond(H 24-2 ). It means that cubic B-spline function matrix can reduce the condition number of the original transfer matrix H 24-2 , and the ill-posedness of the new transfer matrix in the matrix inversion process can be greatly improved. The identification result of case (F4-2, R2) is shown in Fig. 6(b). It shows  that, similar to single-time impact case, the identification result in the non-loading area fluctuates around zero, and the identification result is in good agreement with the measured result in the loading area. For case (F4-2, R2), the total relative error of identification result by RCBSC-ML method is 16.68%. The identified two peak forces are 39.55N and 58.99N, respectively. The corresponding peak relative errors are 1.13% and 1.68%.
Case (F3-3, R1) is conducted to further illustrate the effectiveness of RCBSC-ML method for consecutive impact case. For (F3-3, R1) case, the acceleration response measured at mass m 1 is applied to identify the consecutive three-time impact acted on mass m 3 , and the impact force f 3-3 (t) is simulated by (32). The identification result by RCBSC-ML method is shown in Fig. 7. The optimal parameter m s determined by modified L-curve criterion is 222. The identified three peak forces are 50.72N, 29.37N and 71.90N, respectively. The corresponding peak relative errors are 1.44%, 2.10% and 2.71%. In conclusion, RCBSC method combined with modified L-curve criterion (RCBSC-ML method) is effective for consecutive impact case. + · · · · · · + 70e (−2π×400(t−2.5) 2 ) (32)

V. EXPERIMENTAL VERIFICATION A. EXPERIMENTAL SETUP
In order to validate the performance of the presented RCBSC-ML method in actual experiment environment, a simplified artillery test-bed is established to conduct the experimental testing, and the experimental test system is shown in Fig. 8. The artillery test-bed which is simplified from one kind of real 155mm artillery mainly includes the barrel, cradle and base, wherein the barrel and cradle are fixed on the base through a locking mechanism. When the artillery fires, the projectile or bore pressure would impose the impact force on the barrel. The impact event can cause intense vibration of the barrel. Once this impact event is abnormal, the vibration of artillery would change, which would reduce the firing accuracy. Therefore, the impact force acted on the artillery should be identified so as to monitor the firing performance of artillery. The experimental test system shown in Fig. 8 is mainly composed of the impact hammer, accelerometer, data acquisition instrument, computer and simplified artillery test-bed. In the experiment, the impact hammer is applied to impose impact at the excitation position (point 4 and point 5). The accelerometer is installed at the measurement position (point 1, point 2, and point 3), and the measured dynamic response is utilized to conduct impact force identification. The impact force can also be measured by a force transducer contained in the impact hammer, and the measured impact force is used to illustrate the identification accuracy of RCBSC-ML method. During recording signals by the data acquisition instrument, the sampling frequency is set as 2560 Hz, and acquisition duration is 1s.

B. EXPERIMENTAL RESULTS AND DISCUSSIONS
The impact force identification is conducted in the simplified artillery test-bed shown in Fig. 8. Case (F4, R2) is taken as an instance to validate the effectiveness of RCBSC-ML method in the actual experimental environment. Case (F4, R2) means that the measured acceleration response at point 2 is used to identify the impact force at point 4. The transfer matrix between excitation position (point 4) and response measurement position (point 2) can be obtained by (8) and SI method, and the condition number of the transfer matrix is 2.50×10 18 . It shows that the impact force identification problem on the test-bed is a severely ill-posed inverse problem.
Regularization parameters of RCBSC method determined by modified GCV criterion and modified L-curve criterion are shown in Fig. 9. From the figure, optimal regularization parameters obtained by the above two criterions are  541 and 638, respectively. The corresponding identification results by RCBSC method are shown in Fig. 10. It can be noted from the figure that the experimental identification result is similar to the simulation result. The identification result in the non-loading area fluctuates around zero, and the identification result basically coincides with the measured result in the loading area. For RCBSC-MGCV method, the total relative error of identification result is 17.81%, identified peak force is 77.39N, and the peak relative error is 5.32% (measured peak force is 81.74N). For RCBSC-ML method, the total relative error of identification result is 21.47%, identified peak force is 80.62N, and the peak relative error is 1.37%. It shows that the peak relative error by RCBSC-ML method is smaller than that by RCBSC-MGCV method in the experiment.
In order to further illustrate the advantage of RCBSC-ML method presented in this paper, RCBSC-ML method and classical Tikhonov regularization method are used to conduct impact force identification of case (F5, R3), and the identification result is shown in Fig. 11. The number of collocation points determined by modified L-curve criterion is 550, and Tikhonov regularization parameter determined by original L-curve criterion is 0.0026. The total relative errors by Tikhonov regularization method and RCBSC-ML method are respectively 71.33% and 23.97%. The corresponding peak relative errors are 9.35% and 1.56%. For Tikhonov regularization method, matrix factorization process would reduce the computational efficiency. Additionally, Tikhonov regularization method is also an l 2 -norm regularization method [6]. Therefore, as shown in Fig. 11, the identification result by Tikhonov regularization method also fluctuates around zero in the non-loading area, and this fluctuation amplitude is greater than that for RCBSC-ML method. For RCBSC-ML method, as shown in (20), cubic B-spline matrix is applied to form a new transfer matrix so as to improve the ill-posedness of the matrix inverse operation. The new transfer matrix is a banded sparse matrix, and its dimension is much smaller than that of original transfer matrix. Additionally, RCBSC-ML method does not involve matrix factorization. Therefore, the computational efficiency of RCBSC-ML method is higher than that of Tikhonov regularization method. Furthermore, the profile of the cubic B-spline curve is similar to that of impact force. It means that RCBSC-ML method already has a priori information of approximating impact force before conducting force identification. Therefore, both total relative error and peak relative error of RCBSC-ML method are smaller than those of Tikhonov regularization method.

VI. CONCLUSION
This paper presents a regularized cubic B-spline collocation method combined with the modified L-curve criterion (RCBSC-ML) to conduct impact force identification. Considering that the profile of cubic B-spline is close to that of impact force, a linear combination of cubic B-spline functions can be used to approximate the impact force to be identified. Through the above processing, the unknown variable in the reformed force identification equation changes from impact force vector to weight coefficient vector. The new transfer matrix consisting of original transfer matrix and cubic B-spline matrix is a banded sparse matrix, and its condition number is much small than that of original transfer matrix. Therefore, the ill-posedness of transfer matrix inverse operation is improved. The number of cubic B-spline functions or collocation points (these two parameters can be converted to each other) controls the improvement degree of the ill-posedness, and it is called as regularization parameter of RCBSC method, like Tikhonov regularization parameter and truncation threshold of TSVD method. Modified L-curve criterion is developed to select the optimal regularization parameter of RCBSC method, and it is determined by the largest curvature of L curve which is drawn by the solution norm and residual norm.
Simulation studies including single-time impact case and consecutive impact case on a five DOFs chain type structure are conducted to verify the feasibility of RCBSC-ML method. It also illustrates that the regularization parameter determined by modified GCV criterion is close to that determined by minimum total relative error, and regularization parameter determined by modified L-curve criterion is close to that determined by minimum peak relative error. The peak of impact force is a key index of structural health monitoring, and scholars prefer to improve the identification performance in the peak relative error. Therefore, compared with modified GCV criterion, the presented modified L-curve criterion is more suitable to determine the optimal regularization parameter of RCBSC method. Additionally, some influence factors, including noise level, excitation position and measurement position are fully discussed to show the robustness of RCBSC-ML method. Experimental studies on a simplified artillery test-bed are conducted to verify the effectiveness of RCBSC-ML method. The experiment results also show that, compared with classical Tikhonov regularization method, RCBSC-ML method is more efficient and accurate in impact force identification. In further work, we will develop the modified L-curve criterion to select the optimal regularization parameter of other basis function collocation methods.