Impact Force Localization and Reconstruction via ADMM-based Sparse Regularization Method

In practice, simultaneous impact localization and time history reconstruction can hardly be achieved, due to the ill-posed and under-determined problems induced by the constrained and harsh measuring conditions. Although ℓ 1 regularization can be used to obtain sparse solutions, it tends to underestimate solution amplitudes as a biased estimator. To address this issue, a novel impact force identification method with ℓ p regularization is proposed in this paper, using the alternating direction method of multipliers (ADMM). By decomposing the complex primal problem into sub-problems solvable in parallel via proximal operators, ADMM can address the challenge effectively. To mitigate the sensitivity to regularization parameters, an adaptive regularization parameter is derived based on the K -sparsity strategy. Then, an ADMM-based sparse regularization method is developed, which is capable of handling ℓ p regularization with arbitrary p values using adaptively-updated parameters. The effectiveness and performance of the proposed method are validated on an aircraft skin-like composite structure. Additionally, an investigation into the optimal p value for achieving high-accuracy solutions via ℓ p regularization is conducted. It turns out that ℓ 0.6 regularization consistently yields sparser and more accurate solutions for impact force identification compared to the classic ℓ 1 regularization method. The impact force identification method proposed in this paper can simultaneously reconstruct impact time history with high accuracy and accurately localize the impact using an under-determined sensor configuration.


Introduction
Characterized by exceptional properties of high strengthto-weight ratio, composites have been widely used in such industries as the aerospace to make transportation lighter [1].Unfortunately, composite structures are quite sensitive to impact forces from foreign object debris (FOD) like stones, hail and birds during service, and are then susceptible to the barely visible impact damage (BVID), including delamination and debonding [2,3].Such internal damage is tricky to be detected by the routine procedure like preflight visual inspection, and can wherefore wreak havoc on structural integrity, jeopardizing service security and even cause catastrophic failure [4].To detect the damage at an early stage, impact force monitoring has received considerable attention for its progress in reducing human interference and decreasing maintenance cost [5].It is of great help for monitoring the safety and integrity of composite structures to know and record the impact event such as its time history and location [2].Measuring impact forces directly with transducers tends to be impractical, considering that the excitation sources may be unknown or the excitation location may not be available for the transducer installation.In consequence, an indirect alternative that has come into the spotlight these days is impact force identification (IFI) technique, which combines together the accessible responses collected by easy-to-use sensors like strain gages and the dynamic model of the monitored structure.However, participate to solve the impact forces in this inverse problem accurately, tailored algorithms are required.
Normally, IFI incorporates two main tasks including force history reconstruction and impact localization [6,7].Jacquelin et al. [8] employed truncated singular value decomposition (TSVD) and Tikhonov regularization (TR) to reconstruct impact forces in time domain respectively for a comparative study.Thite et al. [9] conducted the force reconstruction through a plate simulation study in frequency domain using TSVD, and proposed a singular value rejection threshold criterion.Tran et al. [10] developed a deconvolution technique for impact force reconstruction in wavelet domain based on TSVD.Zhang et al. [11] presented a Bayesian regularization (BR) approach for force reconstruction to cope with an uncertain structural model, jointly using Monte Carlo Markov chain to determine the unknown forces.Li et al. [12] dealt with the force identification problem via BR in conjunction with adaptive ℓ q (q = 1 or 2) norm, treat- ing force history, precision parameters and q as random unknowns.When it comes to IFI, TSVD, TR and BR often cannot produce the force history with high accuracy generally, since they penalize solutions smoothly and false forces appear at unloading stage.
In composite structural health monitoring, localization is often considered as the first step in quickly measuring impact forces to avoid fatal damage to the structure.Qiu et al. [13] employed pattern recognition to localize the impact firstly in time domain and then reconstructed its history using the TR method.Yan et al. [14] presented a two-step method to accomplish the IFI of the composite plate, localizing the impact using a nonlinear unscented Kalman filter first and then reconstructing impact force time history by BR.When the number of monitored locations rises, accessible measurements are limited and under-determined cases might be faced.Traditional methods are usually incapable of dealing with these problems, especially for identifying time history and localization simultaneously [15].As an eye-catching regularization method, sparse regularization (SR) explores the inherent sparsity of impact forces in the joint timespace domain [16][17][18], making it very popular in fields such as fault diagnosis [19,20] and image processing [21].Under such strong constraints brought by sparse priors, SR can achieve the force reconstruction and impact localization at the same time, even in under-determined circumstances [22].Ginsberg et al. [17] addressed IFI via SR, and achieves simultaneous localization and reconstruction of the impact force with a basis pursuit denoising algorithm.Also other types of loads, say, moving forces [23] and distributed forces [24] are able to be identified by SR.To further promote the identification accuracy, enhanced SR [25] and group SR [16,26] were proposed respectively.
These aforementioned methods are ℓ 1 -norm based SR, as a relaxation of ℓ 0 -norm, which unavoidably produce biased estimations and underestimate the true solution [27][28][29].Then non-convex penalties are getting more and more attention in place of ℓ 1 -norm to improve the spar- sity and accuracy of IFI [22,[30][31][32].As a typical non-convex regularizer, the quasi-norm ℓ p ( 0 ≤ p < 1 ) has gained popularity in a wide range of fields [27,28] because the non-convex ℓ p regularization performs much better than the convex ℓ 1 regularization in promoting sparsity and obtaining more accurate solutions when p < 1 .Sev- eral typical algorithms have been employed to resolve ℓ p -norm regularized IFI problems.Aucejo [31] introduced an iteratively reweighted least-squares (IRLS) algorithm to solve the multi-parameter multiplicative ℓ p regulariza- tion model of IFI in frequency domain.Qiao [32] adopted an iteratively reweighted ℓ 1 -norm (IRL1) algorithm to tackle an additive ℓ p regularization regularized IFI model in time domain.Both IRLS and IRL1 suffer from several limitations when implemented in IFI: (i) The two algorithms optimize the surrogate functions respectively rather than the original function, so the solution accuracy is limited [33]; (ii) These two algorithms converge to an undesired local minima when an inappropriate initial solution is utilized, so pre-computing the initial solution is generally required, which is time-consuming and unstable [31,32]; (iii) Both algorithms usually admit slow convergence when the transfer matrix is ill-conditioned [33].Thus, the non-convex ℓ p optimization problem still remains difficult for IFI.
In this contribution, in order to address the aforesaid issues existing in ℓ p regularization, a novel ℓ p sparse regularization approach for IFI is proposed based on the principle of ADMM [34], which can simultaneously reconstruct force time-history and localize the impact with an under-determined sensor configuration.The proposed ADMM-based algorithm for IFI shows the following merits: (i) By virtue of ADMM, the primal complicated problem is split into three easier sub-problems, so the large-scale problems caused by multi-point impact monitoring can be solved; (ii) Each subproblem is resolved directly in virtue of proximal operators, which can promote high-accuracy IFI results; (iii) This algorithm ensures that it converges to a stationary point, providing feasible solutions for IFI; (iv) The non-convex optimization requires no well-chosen initial solution in advance; (v) The algorithm solves non-convex sparse regularization problems with arbitrary p values.Nonetheless, the solutions produced by regularization technique are considerably sensitive to the regularization parameter, and inappropriate choices result in poor estimations.Unfortunately, there is still no versatile criteria for regularization parameter selection [35].In this article, an adaptive parameter strategy is proposed for IFI: (i) The IFI problem is translated into a K-sparsity problem by considering the joint space-time sparsity priori; (ii) An adaptive regularization parameter is derived according to the relationship between the sparsity value K and the regularization parameter λ; (iii) The optimal sparsity value K is determined by minimizing the identification mean squared error (MSE) which is feasibly estimated by Monte Carlo Generalized Stein Unbiased Risk Estimate (MC-GSURE) technique [36]; (iv) Once K is decided, an adaptive λ is automatically updated at each iteration and fine-tuning is not required even if the noise level of measurements changes.Finally, an ADMM-based Sparse Regularization method named ADMM-SpaRe is proposed for IFI.
This paper is organized in the following way.In Section 2, inverse analysis for impact force reconstruction and localization is presented as the theoretical background.In Section 3, the ℓ p sparse regularization model for IFI is established, an efficient algorithm termed ADMM-SpaRe is proposed to resolve this model, and parameter setting strategies are discussed in detail.Section 4 presents a series of experimental studies to prove the effectiveness and performance of the proposed method and compares the results of this method with those of the ℓ 0 and ℓ 1 regularization methods.The con- clusions are summarized in Section 5.

Inverse Analysis for Impact Force reconstruction and Localization
For a linear time-invariant (LTI) multi-input/multioutput (MIMO) mechanical system, when subjected to impact forces, its structural responses can be obtained via the convolution integral as where x i (t) represents the response at point i , and the impulse response function (IRF) h ij (t) denotes the lin- ear transmission between the external force at the input position j (j ∈ 1, 2, ..., n f ) and the response at the output position i (i ∈ {1, 2, ..., n r }) in continuous (1) time domain.Zero initial conditions are assumed, i.e., f j (t) = h ij (t) = x i (t) = 0 when t < 0 are assumed.Dis- cretizing the convolution integral with a fixed sampling interval t leads to the following algebraic equation [17,37], in which H ij ∈ ℜ N ×N is a Toeplitz-like matrix and N is the analyzed data length.Eq. ( 2) can be rewritten in matrix-vector form as, For a MIMO dynamic system, a more generalized form of Eq. ( 3) can be expressed as, Then, a compact output form for impact responses can be written as, with the transfer matrix H ∈ ℜ n r •N ×n f •N , the impact force vector f ∈ ℜ n f •N and the measurement vector x ∈ ℜ n r •N .Once f is obtained, both the time history and localization of impact forces are known.
As for the forward problem, once the transfer matrix H and the external force vector f are known, responses x can be readily calculated using Eq.(5).When the transfer matrix H and responses of the mechanical system are known to monitor the external impact forces acting on it, namely solving f in Eq. ( 5), the inverse problem is termed impact force identification (IFI).Generally, according to the relationship between the quantity of inputs and outputs, IFI problems are divided into three categories: 1) the over-determined case when n r > n f ; 2) the even- determined case when n r = n f ; 3) the under-determined case when n r < n f .Whatever the case, IFI may be intrin- sically ill-conditioned, for even minute variations can cause large fluctuations in force identification.Furthermore, this inverse problem might be ill-posed, as the use of noisy and limited measurements cannot consequentially produce a unique and stable solution, especially in the under-determined case, and tackling the underdetermined case is the focus of this work.
Inevitably, noise is introduced into the measurements causing the solution to "runs away" uncontrollably, and it (2) (5) x = Hf is hardly possible to obtain the desired result through classical methods, like QR factorization or least squares.But it is the stable and unique solution that allows for a meaningful interpretation in relation to the basic inverse problem model.Thus, the impact force identification problem is reformulated by numerical treatment.The well-known regularization technique is resorted to, containing additional assumptions of the solution like sparsity.In this work, an ADMM-based ℓ p sparse regularization method with adap- tive regularization parameters is proposed for IFI in the under-determined case.

Formulation of IFI via ℓ p Regularization
To stabilize the resolving of Eq. ( 5) for f , since the inverse problem may be ill-conditioned or ill-posed (e.g., the condition number of H is quite large or H is singular), the ℓ p regularization technique is brought in, where the first term 1  2 �x − Hf � 2 2 measures the differences between the linear model Hf and the output x , the second term f p p measures the sparsity of f , and the positive constant is the so-called regularization parameter.When p = 2 , this model is known as the Tikhonov regularization and is one of the most widely employed models for regularization of linear inverse problems.When p = 1 , the model is called the least absolute shrinkage and selection operator(LASSO).Here, the ℓ p - norm ( 0 < p ≤ 1 ) is defined loosely, Additionally, when p = 0 , the ℓ 0 -norm defines that f 0 is the quantity of non-zero entries of f .These quasi- norms have been extensively applied in many fields thanks to their excellent properties [27,28,38].
The regularization method based on ℓ p -norm is hard to be optimized.Especially, the resulting optimization problem is not convex while 0 ≤ p < 1 .The non-convex prob- lem is usually intractable as undesired local minima tend to be acquired reluctantly.To deal with the ℓ p regulari- zation model for IFI, an efficient algorithm is developed under ADMM framework.

Reformulating the IFI Model under ADMM Principle
As one of splitting methods, ADMM is a simple but powerful optimization algorithm that works for both convex and non-convex problems [34].By means of ADMM, the problem ( 6) can be reformulated as, where the data fidelity term , the penalty term G(z)= �z� p p , and an additional vector of variable z is introduced for variable splitting.Then, the augmented Lagrangian can be derived as [34], where u is a scaled vector of Lagrangian multipliers, and ρ denotes the positive penalty parameter which can be adjusted to promote the rate of convergence.Minimizing L ρ (f , z, u) with respect to f and then with respect to z , and conducting a multiplier update, yields three steps, where f-update and z-update can be accomplished by computing the proximal operator of the convex quadratic term in Eq. ( 10) and the ℓ p -norm ( 0 ≤ p ≤ 1 ) in Eq. ( 11), respectively, which will be studied in the following section.Then the u-update shall be carried out through Eq. ( 7) 12).A critical point of the augmented Lagrangian can be reached in the end [39].

Resolving Subproblems Using Proximal Operators
The subproblems in Eq. (10) and Eq. ( 11) can be handled efficiently by proximal operators.Commonly, given a function h(v) , its proximal operator is given by [40]: where w denotes the input vector.To begin with, for f -update, the proximal operator of the convex quadratic function h which is equivalent to the minimization in Eq. (10).The optimal solution of Eq. ( 14) is calculated by taking the derivative with respect to f and setting it equal to 0 , that is, and hence the f-update can be gained in a closed-form as In the end, a compact expression of f-update step leads to, (13)  The z-update can be accomplished by computing the proximal operator of ℓ p -norm(0 ≤ p ≤ 1 ).To make the description easy to follow, a general proximal operator form is studied with h(v)=�v� p p scaled by µ , which is, (17) Different p s make for different solutions of Eq. ( 18), and detailed studies come as follows.
I) When p = 0 , the corresponding proximal operator known as hard-thresholding [38] is defined componentwise as: where w i represents the ith element in vector w , and T (µ) is the threshold, where T (µ) = √ 2µ.II) When p = 1 , the proximal operator is the familiar soft-thresholding operator [41], where T (µ) = µ.
III) When 0 < p < 1 , the corresponding proximal operator has no analytic solution, but can be resolved via an inexact proximal operator called generalized shrinkage / thresholding (GST) operator [28], where T p (µ)=(2µ(1 − p)) is the thresholding value derived in Ref. [28], and .., M .Empirically M is chosen as 2, and v (0) i is initialized as |w i | .Notably, the hard-thresholding and soft-threholding functions are special cases of GST with p = 0 and p = 1 [28], and when p = 1/2 or 2/3 , a closed-form solution can be deduced [35,42].Hence, considering the generalized p values, the GST thresholding is utilized for ℓ p reg- ularization herein.Accordingly, combining Eq. ( 11) with Eq. ( 18), and after some variable substitution, especially v = z , w = f (k+1) + u (k) and µ = ρ , the proximal operator for the z-update is given by, and then the update step for z can be compactly written as, By virtue of the ADMM, the optimization for problem (8) can be handled succinctly by three main steps, briefly indicated in Eq. ( 17), Eq. ( 23) and Eq. ( 12). (

Deriving the Adaptive Regularization Parameter Using K-sparsity
As we all know, the results of regularization techniques are quite sensitive and vulnerable to the selection of the regularization parameter .And it is generally difficult to choose the optimal which mainly depends on the noise level and the transfer matrix [43].Once noise levels or transfer matrices change, careful and energydraining operations to find optimal are inevitable.To overcome these obstacles, a robust and adaptive strategy for the selection is developed according to K-sparsity [35], because IFI is actually a K-sparsity problem in time domain, thanks to its joint space-time sparsity priori.Specifically speaking, the z-update in Eq. ( 23) is involved, and z expects the identical sparsity as f according to the constraints in Eq. ( 8).Instead of a fixed , the can be set adaptively in light of the K-sparsity priori, that is, the adaptive threshold T is set to the K th largest element in the input vector w (k) = f (k+1) + u (k) (in absolute value) at each iteration [35], which is expressed as where w (k)  [K ] denotes the K th largest component (in absolute value) in the vector w (k) at the kth iteration.Resolving Eq. ( 24) can lead to a reliable choice of at each iteration, and here an analytic solution of can be gained as, Finally, the derived ℓ p regularization algorithm with adaptive regularization parameters named ADMM-SpaRe for IFI is stated in Algorithm 1 (shown in Table 1).Besides, the Algorithm 1 stops when the relative maximum change of adjacent iterations becomes less than a small value ε , say ε=10 −5 here, i.e., ∞ .Note that early termination in the f-update and z-update steps can even give (24) an overall enhancement in algorithm efficiency [34,39,44].Presetting K makes an adaptive threshold instead of a fixed one.Note that in applications the parameter K can be decided exactly based on the physical prior, or can be set as an upper-bound estimate of the sparsity under discussion via simulation/experimentation studies [35].In this contribution, K is selected by minimiz- ing MC-GSURE MSE.

Parameter Setting
There are three parameters including 0 ≤ p ≤ 1 , ρ > 0 , K > 0 to be determined in Algorithm 1. Firstly, taking sparsity as a prior can bring strong constraints and then reduce the degrees of freedom of the IFI model, which makes it feasible to tackle the under-determined case.Generally speaking, setting 0 ≤ p ≤ 1 makes sparse solu- tions accessible for IFI, but ℓ 1 regularization often pro- duces biased estimates.Numerous studies point out that when p < 1 , the non-convex ℓ p regularization can pro- mote sparser and more accurate results than the convex ℓ 1 regularization, but the optimal p is usually different in different applications [28].For example, ℓ 2/3 performs better than ℓ 1/2 in image deconvolution [42].There- fore, the optimal p is studied by investigating different p s in a fine division between 0 and 1 in the following experiments.
Secondly, for the sake of improving convergence, the penalty parameter ρ is adjusted using residual balancing strategy [45], i.e., where r (k) = f (k) − z (k) and s (k) = − ρ(z (k+1) − z (k) ) are the primal and dual residuals in the kth iteration of ADMM-SpaRe, respectively.( 26) For a better choice of K in practice, a practical and reli- able thought is based on minimizing MSE.However, the identification MSE can not be known in practice.MC-GSURE is used to estimate MSE availably.Due to a rankdeficient model dealt with in the under-determined IFI case, generalized SURE to get an unbiased estimate of MSE is defined as [36], where P is a projection matrix P := H T (HH T ) −1 H , h K (u) returns the identification results of Eq. ( 6) using K-sparsity, u = 1 σ 2 H T x ( σ is the standard deviation of noise) represents the sufficient statistic for the model, fML = H T (HH T ) −1 x is the maximum likelihood estima- tion, and the divergence div u (Ph K (u)) is approximated by the MC method as [46]: where b is a Gaussian vector with zero mean and unit standard deviation, and δ is a small positive parameter.Thus far, an overall flowchart of the proposed ADMM-SpaRe algorithm for IFI is given in Figure 1.

Experimental Validation
In this section, a series of laboratory experiments are implemented on a composite panel to verify the effectiveness and investigate the performance of the proposed ADMM-SpaRe method for IFI tasks.Impact forces are identified utilizing only structural responses, and impact events can be localized and reconstructed concurrently based on an under-determined sensor configuration.The optimal sparsity value K is investi- gated by calculating MC-GSURE MSE.A representative p is chosen for IFI through experimental study.The robustness of ADMM-SpaRe to noise interference is also examined.

Experimental Setup
An aircraft skin-like composite laminated plate, which is fabricated by T700/QY8911 carbon fiber reinforced pre-impregnated unidirectional layers and clamped at both ends, is limned in Figure 2, with the size of 400 × 300 × 6 mm and the lay-up: [+45/0/ − 45/90] 6s . ( Firstly, as shown in Figure 2(a), impact forces are produced by an instrumented hammer with a force transducer (PCB 086C03, sensitivity 2.25 mV/N), and the vibration data (i.e., strain responses) are measured by ten stain gages (PCB 740B02).Only two strain gages (strain gage 1 and strain gage 9, with sensitivity 52.2 mV/g and 48.6 mV/g, separately) are utilized here for IFI. Figure 2(b) visualizes the spatial arrangement of the key points on the laminate: nine uniformly distributed points are monitored as potential impact locations; ten strain gages are equally mounted on both sides of the monitoring area.All of these points fall on the grid with size of 50 × 50 mm.Secondly, the data including impact forces and strain responses are recorded by LMS Test.Lab system synchronously with a sampling frequency of 10240 Hz.Thirdly, the transfer matrix H is constructed experi- mentally, which is briefly divided into two steps: (i) The IRF h ij (t) between the output location i and input loca- tion j can be gained by means of inverse fast Fourier transform (IFFT) of the frequency response function (FRF) h ij (w) which can be acquired by impact testing via the LMS modal testing module; (ii) The block Toeplitz-like matrix H in Eq. ( 5) is assembled according to Eq. ( 4).The analyzed data length N of each discrete IRF is 512, making each submatrix in Eq. ( 4) have a dimension of 512.Herein, only two strain gages are employed to monitor the nine potential impact positions, resulting in a quite under-determined system with the transfer matrix H ∈ ℜ 1024×4608 .Finally, all experimental analyses are conducted under MATLAB 2018b and on a personal computer equipped with Intel (R) Core (TM) i5-8400, 24 GB RAM.
To evaluate the performance of ADMM-SpaRe quantitatively, several accuracy assessment indicators are defined.Note that the impact force measured by the force transducer attached on the hammer is taken as reference in the experiments.Firstly, the global relative error (RE) between the identified force vector f and the referenced one f is calculated as: and then, the RE at each location i can be computed as: Next, the peak relative error (PRE) at the impact location i , assessing the reconstruction quality of the impact amplitude which matters much for IFI, is expressed as: ( 29) Finally, an indicator LA i at the impact location i , evalu- ating the localization accuracy (LA), is represented as: In this work, it is assumed that one impact force is exerted on one of the potential impact locations during the observation period.Theoretically, LA i should be one for the actual impact position, and be zero for nonimpact locations.

Sparsity Estimation via MC-GSURE
As aforementioned, once K is set, adaptive regularization parameters are updated automatically at each iteration according to Eq. ( 25).K is actually determined by the sparse physical priori of the impact forces, and requires no adjustment regardless of changes in the noise level of the measurement.Next, how to choose the optimal K and why K-sparsity is used instead of a fixed are respec- tively stated through experimental studies.Without loss of generality, p is equal to 1/2 and the true impact force acts on #1. 31 different values are allocated to K and in a range of 0.01 − 1 (in increments of 0.033) and 2 − 60 (in increments of 2) respectively.For MC-GSURE MSE computing, the noise level of the measurement can be obtained by the well-known mean absolute derivative (MAD) rule [47].For the MC process in Eq. ( 28), δ is set to 10 −5 and an average of five independent realizations (31 are taken.Then the optimal and K are both found by minimizing MSE, as shown in Figure 3. From Figure 3(a1) and (a2), firstly, pretty accurate estimates of MSE can be acquired by MC-GSURE.Secondly, for K-sparsity, the optimal K = 10 is obtained, exactly matching with the sparsity of the true impact force.From Figure 3(b1), when K is larger than the real sparsity, a wide range of overestimations of K can lead to stable and high-accuracy solutions, and when K = 60 , six times the real K , the accurate identification accuracy can be real- ized with RE less than 10%.Looking into Figure 3(b2), the fixed strategy is not comparable to the K-sparsity strategy with lower accuracy.According to the current division of , there is a small adjustment range and consequent results are unstable.What is noteworthy is that once K is set, regularization parameters are updated automatically without considering noise levels.However, when the fixed strategy is used, the optimal needs to be adjusted continually according to the corresponding noise level [43].As a result, the K-sparsity strategy based on MC-GSURE is employed in the following experimental study.

The Optimal Selection for p
Regrading ℓ p regularization, the optimal p varies in dif- ferent application areas [28,48].Thus, the optimal p for IFI is investigated.The impact forces applied at nine different locations in Figure 2(b) are individually identified by ADMM-SpaRe with p ranging from 0 to 1 at the inter- val of 0.1.Then, the box-plots in Figure 4 give the identified results indicated by REs, PREs, and LA i s, and each box visually shows the identification accuracy of all nine monitored positions under each p value.The box-plots say that ADMM-SpaRe possesses good identification accuracy and stability when p = 0.4 − 0.6 .To be specific, when setting p = 0.4, 0.5 and 0.6, induced REs and PREs are less than 20% and 10% respectively, and LA i s are more than 0.9, which are hard to be guaranteed by other p settings, especially when p = 0 or 1.Moreover, when p = 0.6 , ADMM-SpaRe achieves the highest identification accuracy.An exemplary reconstruction accuracy on position #7 are, respectively, REi ℓ 0.4 = 13.7% , REi ℓ 0.5 = 12.9% , REi ℓ 0.6 = 12.4% , PREi ℓ 0.4 = 4.5% , PREi ℓ 0.5 = 5.3% , and PREi ℓ 0.6 = 4.4% , and ℓ 0.6 regularization achieves the highest identifica- tion accuracy.Therefore, p = 0.6 can be taken as the pre- ferred presentative of ℓ p regularization for IFI because of its high and stable identification accuracy.

Robustness to Noise Interference
Low noise levels are involved in laboratory experiments indeed.In order to reveal the robustness of ADMM-SpaRe-ℓ 0.6 regularization to noise interference, a noise model is taken into consideration as follows: where x and x ∈ ℜ N denote the original measured sig- nals and signals with additive noise, β is the noise level, respectively, σ signifies the standard deviation of the measured data, and randn(N , 1) generates values from a normal distribution with zero mean and unit standard deviation.In this study case, the real impact force is exerted on location #9 and the strain responses are still collected by sensors #1 and #9 to ensure generality.Subsequently, β accounts for different noise levels, vary- ing from 0 to 50% in 2% increments.IFI is realized using ADMM-SpaRe-ℓ 0.6 , and the outcomes of ℓ 0 and ℓ 1 regu- larization are also considered for comparison, solved by the well-known hard and soft thresholding formulas, respectively.To eliminate the randomness, 100 independent realizations are performed for different methods under each noise level.Figure 6 shows the average values of indicators including RE, PRE i , and LA i .As described in Figure 6, ℓ 0.6 obtains the best robust- ness behavior over ℓ 0 and ℓ 1 .In particular, despite the increase of noise, ℓ 0.6 regularization can achieve high performances under 50% noise level, keeping the REs less than 10%, the PREs less than 5%, and the LA i s more than (33) x = x + βσ randn(N , 1), 0.95.Looking into the results of ℓ 1 regularization, large identification errors appear under each noise level, inferior to ℓ 0.6 regularization.Besides, reasonable identifi- cation results can be acquired through ℓ 0 regularization under low noise level, i.e., under less than 15%, as displayed in Figure 6, while identification precision rapidly declines as the noise level rises up.
To better observe the robustness to noise, the indicators representing the identification accuracy of ℓ p (0 ≤ p ≤ 1) regularization are displayed in Figure 7 under the 50% noise level after 100 independent realizations.It can be discovered that ℓ 0.6 regularization gains the highest iden- tification accuracy, while the accuracy decreases with ps away from 0.6.The noisy signals under the 50% noise level are shown in Figure 8. Correspondingly, detailed IFI results of ℓ p ( p = 0, 0.6, 1 ), covering force time history and impact localization, are illustrated in Figure 9.It can be found visually that spurious forces appear on the nonimpact positions, especially in the results produced by ℓ 0 and ℓ 1 regularization, and by contrast, ℓ 0.6 regularization can perform better to alleviate this problem.In short, it is justified that ℓ 0.6 regularization is more powerful than ℓ 0 and ℓ 1 regularization, and ℓ 0 regularization is competitive to ℓ 1 regularization under relatively low noise level.

IFI Results and Discussion
In this case, for the sake of generality, impact forces generated by the impact hammer with a steel hammer cap are applied to a randomly selected location #4 from nine potential locations as displayed in Figure 2(b).Also, the responses are collected by strain gage #1 and #9, as shown in Figure 10.As previously mentioned, the noise levels are low in experimental scenarios.IFI is accomplished by the proposed ADMM-SpaRe-ℓ p method with p = 0.6 , and additionally the results produced by ℓ 0 and ℓ 1 regularization are taken as comparisons.

Table 1 ADMM-SpaRe algorithm pseudo code
These are not hard to be found from the local analysis of identification results as exhibited in Figure 11  In experimental tests, an interesting phenomenon that cannot be ignored is that double-impact sometimes occurs.In this case,the double-impact produced by the impact hammer on ♯4 is identified using the same three regularization methods.Figure 12 illustrates the strain responses gathered by strain gages ♯1 and ♯9 separately, and it is hard to tell whether a double-impact occurs.On this occasion, ℓ 0.6 regularization still retains the highest identification accuracy with RE ℓ 0.6 = 16.9% , compared to RE ℓ 0 = 18.7% for ℓ 0 regularization.Notably, ℓ 1 regu- larization performs unsatisfyingly with RE ℓ 1 = 60.1% .Then, from Figure 13(b1), ℓ 0.6 regularization is able to fully realize the localization with LAi ℓ 0.6 = 1 .How- ever, as can be seen from Figure 13(b2) and (b3), ℓ 0 and ℓ 1 regularization fail to do so due to spurious forces appearing at the unloading positions, with LAi ℓ 0 = 0.89 and LAi ℓ 1 = 0.87 .Through local analysis of the identi- fied results, it can be discovered that the time history of impact forces can be well-recovered by ℓ 0.6 regularization with REi ℓ 0.6 = 16.9% , and next comes ℓ 0 regularization with REi ℓ 0 = 17.1% .While ℓ 1 regularization fails with REi ℓ 1 = 59.5% .Overall, the ℓ 0.6 regularization method solved by ADMMSpaRe is more powerful than the other two methods, and ℓ 0 regularization is comparable to ℓ 0.6 regularization in this case, but ℓ 1 regularization cannot manage to realize the identification under this doubleimpact circumstance.
In order to further verify the performance of the proposed method in identifying impact forces with different shapes, a rubber hammer cap and a plastic hammer cap are used on the impact hammer to impact the #4 position respectively (Figure 13).
In the case of rubber hammer cap, the impact force localization and time-history reconstruction results of ℓ 0.6 , ℓ 0 and ℓ 1 regularization methods are shown in Fig- ure 14.It can be seen from Figure 14(a1-a3) that all three methods can accurately localize the impact force, but the global identification result of ℓ 0.6 regularization is most consistent with the spatial sparsity of the impact force, with RE ℓ 0.6 = 4.0% and LAi ℓ 0.6 = 1.00 .The global relative errors and localization accuracy of the other two methods are RE ℓ 0 = 6.8% , RE ℓ 1 = 10.7% , and LAi ℓ 0 = 0.94 , LAi ℓ 1 = 0.98 , respectively.It can be observed from the reconstruction results Figure 14(b1-b3) that the reconstruction accuracy of ℓ 0.6 regularization is the highest, and the identified impact amplitude is also closest to the real force, with REi ℓ 0.6 = 4.0% and PREi ℓ 0.6 = 0.2% .In contrast, the time-history reconstruction errors of ℓ 0 and ℓ 1 regularization methods are REi ℓ 0 = 6.3% , REi ℓ 1 = 10.6% , PREi ℓ 0 = 1.6% , and PREi ℓ 1 = 5.4%.
In the case of plastic hammer cap, the impact force identification results of ℓ 0.6 , ℓ 0 and ℓ 1 regularization methods are shown in Figure 15.The same conclusion can be drawn, that is, the impact force localization and reconstruction results of ℓ 0.6 regularization match the real impact situation best.By quantifying the quality of the identification results of the three methods, the identification errors of these three methods can be obtained as listed in Table 2. Therefore, it can be concluded that ℓ 0.6 regularization can reconstruct and localize impact forces of different shapes with higher accuracy.

Conclusions
This paper develops an ADMM-based ℓ p sparse regu- larization method termed ADMM-SpaRe to cope with the ill-posed IFI inverse problem.Unlike the previous algorithms which approximate the minimum of the nonconvex ℓ p -minimization by optimizing the proxy functions.ADMM-SpaRe splits the complicated and variable-coupled IFI model into three easier subproblems, and each subproblem can be resolved directly by its proximal operator with high accuracy.Then force reconstruction and impact localization are realized simultaneously with the under-determined sensor configuration where merely two strain gages monitoring nine potential impact locations.To tackle the regularization parameter selection problem, the K -sparsity strategy on account of the joint space-time sparsity priori of IFI is adopted.The optimal K can be accurately found by minimizing the MC-GSURE MSE.A wider range of K s can be chosen to achieve accurate IFI, and once K is set, regularization parameters are updated automatically at each iteration and require no constant fine-tuning, making ADMM-SpaRe an adaptive and robust approach.Moreover, the optimal p = 0.6 is rec- ommended for the IFI problem through experimental comparative studies, which can promote sparser and more accurate solutions.Then, to verify the robustness of ADMM-SpaRe to noise, different levels of white Gaussian noise are added into the measured signals.ℓ 0.6 regularization can accomplish relatively accurate IFI even under 50% noise level.Finally, additional impact scenarios including single impact and double impact, rubber hammer cap impact, and plastic hammer cap impact are investigated.It is found that ℓ 0.6 regulariza- tion adopting GST thresholding performs better than the ℓ 0 and ℓ 1 regularization conducted using classic soft thresholding and hard thresholding.In addition, the proposed method is applicable to the impact force identification of different impact times and impact morphologies.

Figure 1 AFigure 2
Figure 1 A brief flowchart of ADMM-SpaRe for IFI

Figure 4
Figure 4 Identification results of nine monitored locations under varying p values: (a) REs, (b) PREs, (c) LA i s

Figure 10 Figure 11
Figure 10 Strain responses in the single-impact case from (a) strain gage #1 and (b) strain gage #9

Figure 12
Figure12 Strain responses in the double-impact case from (a) strain gage #1 and (b) strain gage #9

Table 2
IFI result errors of plastic hammer cap impact case (%)