Impact force identification on composite panels using fully overlapping group sparsity based on L p -norm regularization

Abstract

Impact force identification is of great importance for composite structural health monitoring due to the poor impact resistance of composite materials. Convex sparse regularization method based on L₁-norm tends to underestimate the amplitude of the impact force. This paper proposes a novel method using fully overlapping group sparsity based on L_p-norm regularization (FOGSL_p) for impact force identification, which can localize the impact force and reconstruct its time history simultaneously with limited measurements in under-determined cases. The FOGSL_p method takes more sparse prior information into account by combining the non-convex L_p-norm ( $0<p<1$ ) and fully overlapping group sparsity to promote sparse solutions, thus improving the accuracy of the reconstructed impact force. An accelerated grouped shrinkage-thresholding algorithm is employed to solve the non-convex optimization problem under the majorization-minimization framework. The Nesterov’s acceleration strategy is modified to accommodate the requirements of non-convex optimization. Simulations and experiments are conducted on composite panels to validate the effectiveness of the FOGSL_p method. Results demonstrate the efficiency and robustness of the FOGSL_p method to localize the impact force and reconstruct its time history simultaneously while 25 potential impact points are monitored using two sensors. Compared with L₁-norm, L_p-norm, and L_2,1-norm regularization methods, the proposed method performs best in both simulations and experiments.

Introduction

Composite materials are widely used in aerospace field to reduce aircraft weight and improve strength. In the practical operation environment, an aircraft composite structure may be subjected to impact forces on the ground and during takeoff. Typical sources of impact include hail, bird strikes, tool drops, runway debris, and ground service equipment.^1,2 However, the impact resistance of composite materials is poor, ³ which may result in barely visible impact damage such as delamination and fiber breakage inside composite materials,^4,5 severely reducing the strength and loading-bearing capacity of composite structures. Therefore, it is of great significance for safe service of aircraft to monitor impact forces acting on composite structures, ⁶ but the impact forces are difficult to be measured directly in practice. In recent years, impact force identification methods have been extensively studied, which localize and reconstruct impact forces based on structural properties like frequency response function (FRF) and measured response.

The problem of impact force identification is typically an ill-conditioned inverse problem, which means that minor errors in response measurements can cause the results of force identification to deviate significantly from the accurate ones. Therefore, various regularization methods are developed to solve the force identification problem. ⁷ For force identification problems with known positions, traditional regularization methods such as Tikhonov regularization are mainly used to reconstruct the time history of the force. Wang et al. ⁸ proposed a sinc-dictionaries-based Tikhonov regularization approach that circumvented the ill-conditions of force identification, which substantially improved the computation efficiency compared with the traditional Tikhonov regularization method. Kalhori et al. ⁹ used the Tikhonov regularization method to solve the inverse problem of impact force identification in the under-determined case, and results indicated that the Tikhonov regularization method performed poorly in solving the under-determined problem, and failed to localize the impact force.

In order to reconstruct the time history of the impact force more accurately, sparse regularization methods based on L₁-norm have gained more attention in the field of impact force identification. Such methods utilize the sparsity of the impact force in the time domain as prior information and enhance the stability and accuracy of the solution of the ill-posed inverse problem. In the framework of L₁-norm sparse regularization, Ginsberg et al. ¹⁰ used three parameters, including the pulse width, pulse magnitude value, and pulse magnitude arrival time, to construct a sparse representation dictionary of impact forces, and an impact force sparse identification method based on basic pursuit denoising was proposed. Huang et al. ¹¹ employed cubic B-spline scaling functions and the TWIST algorithm as pre-processing and post-processing, respectively, to improve impact force identification accuracy with higher computational efficiency. Although the standard sparse regularization method based on the L₁-norm is an easy-to-solve convex optimization problem and improves the accuracy significantly compared with traditional regularization methods like Tikhonov regularization, the problem of insufficient sparsity-promoting ability remains to be solved, as it often results in underestimation of the amplitude of impact forces. ¹²

It is noted that the abovementioned studies only focused on the reconstruction of impact force time history, but the localization of impact force is also a critical part of impact force identification. ¹³ Ciampa et al. ¹⁴ localized the impact sources on composite structures with unknown layup and cross section based on the differences of the stress waves, but the reconstruction of impact force time history was not considered. Li et al. ¹⁵ proposed a deep metric learning method for impact force localization, and the weighted fusion Tikhonov regularized total least squares method was used for impact force reconstruction. Qiu et al. ¹⁶ localized impact forces on a steel panel by using the pattern recognition method and reconstructed the time history of impact force via the Tikhonov regularization method, which are two separate processes. Chen et al. ¹⁷ used two different deep neural networks for localization and time history reconstruction of the impact force separately, while a large amount of data was required to train the networks.

In practical engineering, it makes more sense to achieve accurate localization and time history reconstruction simultaneously for impact forces at unknown locations. For the purpose of taking full advantage of the sparsity of impact force, some non-convex regularizers have been taken into account. Aucejo and De Smet ¹⁸ proposed a multiplicative regularization for force reconstruction and localization, which was able to adjust the parameter p in L_p-norm. Qiao et al. ¹² developed a non-convex optimization model for the under-determined problem of impact force identification based on L_p-norm, using one sensor to monitor nine potential impact locations, and an iteratively reweighed L₁-norm algorithm was introduced to solve the non-convex model. Moreover, considering the impact force is a continuous process over an extremely brief time, it has an inherent group sparsity. Zhang et al. ¹⁹ proposed a new method for moving force identification by combining group sparsity theory and compressed sensing, and signicantly improved the accuracy in practice. Wambacq et al. ²⁰ achieved the localization of impact forces and reconstruction of time histories with limited measurements by taking L_2,1-norm as penalty and proposed an interior point method to improve computational efficiency. From the above research, regularization methods utilizing more prior information can impose stronger constraints on the results of impact force identification, thereby improving the solution accuracy.

Inspired by the successful applications of the non-convex L_p-norm regularization and group sparsity in multiple fields,^12,21–24 this paper develops a novel fully overlapping group sparsity based on non-convex L_p-norm regularization (FOGSL_p) for impact force identification, aiming to localize impact force and reconstruct its time history simultaneously. Firstly, more prior information of sparsity are taken into account by combining the L_p-norm and overlapping group sparsity, which not only enhances the sparsity and improves the accuracy of reconstructed amplitude but also strengthens the ability to solve under-determined problems. Then under the majorization-minimization (MM) framework, an accelerated grouped shrinkage-thresholding algorithm (AGSTA) is proposed to solve the non-convex optimization problem. The performance and effectiveness of the proposed method are verified in simulations and experiments, respectively.

Mathematical description of impact force identification

Inverse problem of impact force identification

For a linear time-invariant system, the relationship between the impact force and the response can be defined as the following Fredholm integral equation:

$$\underset{0}{\overset{t_f}{\int }}k(i,j,\tau ,t)f(j,t)\mathrm{d}t=y(i,\tau )$$

(1)

where $f(j,t)$ represents the unknown impact force acting at position j; $k(i,j,\tau ,t)$ is a kernel function; $y(i,\tau )$ is the response such as strain, displacement, and acceleration, etc., at position i and moment τ; [0, $t_f$ ] represents the measurement time interval. For the impact force identification problem, the kernel function of Equation (1) is a convolution kernel, which satisfies $k(i,j,\tau ,t)=h_{\mathrm{ij}}(t-\tau )$ , and $h_{\mathrm{ij}}(t)$ is the impulse response function (IRF) between the impact position j and response measurement position i. Therefore, the convolution relationship between the structural response and the impact force can be expressed by

$$y(t)=h(t)\otimes f(t)=\underset{0}{\overset{t}{\int }}h(t-\tau )f(\tau )\mathrm{d}\tau$$

(2)

where ⊗ is the convolution operation symbol, and τ is the time-shift factor. Equation (2) can be discretized in the time domain for numerical computation as

$$\begin{array}{c}\left[\begin{array}{c}y(\mathrm{d}t)\hfill \\ \multicolumn{1}{c}{y(2\mathrm{d}t)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{y((N-1)\mathrm{d}t)}\\ \multicolumn{1}{c}{y(N\mathrm{d}t)}\end{array}\right]=\hfill \\ \left[\begin{array}{ccccc}h(\mathrm{d}t)\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{h(2\mathrm{d}t)}& h(\mathrm{d}t)\hfill & \cdots \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \cdots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{h((N-1)\mathrm{d}t)}& h((N-2)\mathrm{d}t)\hfill & \cdots \hfill & h(\mathrm{d}t)\hfill & 0\hfill \\ \multicolumn{1}{c}{h(N\mathrm{d}t)}& h((N-1)\mathrm{d}t)\hfill & \cdots \hfill & h(2\mathrm{d}t)\hfill & h(\mathrm{d}t)\hfill \end{array}\right]\hfill \\ \left[\begin{array}{c}f(\mathrm{d}t)\hfill \\ \multicolumn{1}{c}{f(2\mathrm{d}t)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{f((N-1)\mathrm{d}t)}\\ \multicolumn{1}{c}{f(N\mathrm{d}t)}\end{array}\right]\hfill \end{array}$$

(3)

where dt is the sampling time interval, and N is the length of the measured response. Equation (3) can be expressed in a compact matrix-vector form

$${\mathbf{y}}_i={\mathbf{H}}_{\mathrm{ij}}{\mathbf{f}}_j$$

(4)

where ${\mathbf{y}}_i\in {\mathbf{R}}^N$ represents the response, ${\mathbf{f}}_j\in {\mathbf{R}}^N$ is the force vector, and ${\mathbf{H}}_{\mathrm{ij}}\in {\mathbf{R}}^{N\times N}$ is the transfer matrix of a lower triangular Toeplitz matrix between the impact position and response position. For a linear system, the responses of multiple impact forces are the linear superposition of responses when each impact force acts individually, that is,

$${\mathbf{y}}_i=\sum _{j=1}^{N_f}{\mathbf{H}}_{\mathrm{ij}}{\mathbf{f}}_j$$

(5)

where $N_f$ is the number of impact forces. For a more general Multiple-Input Multiple-Output (MIMO) system, Equation (5) can be extended to the following equation:

$$\left[\begin{array}{c}{\mathbf{y}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{y}}_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{y}}_{N_r}}\end{array}\right]=\left[\begin{array}{cccc}{\mathbf{H}}_{11}\hfill & {\mathbf{H}}_{12}\hfill & \cdots \hfill & {\mathbf{H}}_{1N_f}\hfill \\ \multicolumn{1}{c}{{\mathbf{H}}_{21}}& {\mathbf{H}}_{22}\hfill & \cdots \hfill & {\mathbf{H}}_{2N_f}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \cdots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathbf{H}}_{N_{r}1}}& {\mathbf{H}}_{N_{r}2}\hfill & \cdots \hfill & {\mathbf{H}}_{N_r{N}_f}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{f}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{f}}_{N_f}}\end{array}\right]$$

(6)

where ${\mathbf{H}}_{\mathrm{ij}}$ is the transfer matrix between the impact position i and the response position j, and $N_r$ is the number of response measurement positions. Also, Equation (6) can be rewritten in a compact matrix-vector form

$$\mathbf{y}=\mathbf{Hf}+\mathbf{w}$$

(7)

where $\mathbf{H}\in {\mathbf{R}}^{N_{r}N\times N_{f}N}$ represents the transfer matrix of the MIMO system, which is a Toeplitz block matrix, and w represents the inevitable noise or measurement error. When y and H are known and f needs to be solved, it is a typically inverse problem. According to the number of impact forces $N_f$ and responses $N_r$ , this inverse problem can be divided into three types: (1) when $N_r<N_f$ , it is an under-determined problem; (2) when $N_r=N_f$ , it is an even-determined problem; (3) when $N_r>N_f$ , it is an over-determined problem. In practical applications, it is often impractical to install abundant sensors to construct an even-determined or over-determined system, so we focus on under-determined cases in this paper. The least squares method solves Equation (7) by minimizing the residual term

$$\underset{\mathbf{f}}{\mathrm{minimize}}\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2$$

(8)

where $\Vert ·{\Vert }_2$ represents the L₂-norm of a vector. However, the under-determined problem has infinite solutions, and the accurate solutions cannot be selected by least squares method since more prior information is required to obtain stable solutions.

L ₁-norm regularization

Since under-determined problems such as Equation (8) usually have infinite solutions, regularization methods solve this problem by adding constraints to the least squares formula as penalties. The traditional Tikhonov regularization method based on L₂-norm corrects the singular values of the transfer matrix H by adding a penalty to the least squares method, thereby improving the ill-condition of the inverse problem of impact force identification, as shown in the following formula:

$$\underset{\mathbf{f}}{\mathrm{minimize}}(\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda \Vert \mathbf{f}{\Vert }_2^2)$$

(9)

where $\lambda >0$ is the regularization parameter that balances the residual and penalty terms. There exists an analytical solution to Equation (9). However, the solution obtained by the Tikhonov regularization method is usually non-sparse, and large values may also appear in the non-loading stage, ¹³ thereby decreasing the impact force identification accuracy, indicating the need to introduce stronger priors. Fortunately, the impact force is sparse in the time-space domain, that is, the vast majority of elements in the impact force vector f are close to or equal to zero. Theoretically, the optimal sparse solution can be obtained by minimizing L₀-norm

$$\underset{\mathbf{f}}{\mathrm{minimize}}(\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_0)$$

(10)

where $\Vert ·{\Vert }_0$ is the L₀ pseudo norm representing the number of non-zero elements in a vector. To directly solve Equation (10), it is necessary to filter out all possible non-zero elements in the impact force vector f, resulting in a non-deterministic Polynomial (NP) hard problem for large-scale impact force identification problem, because the search space is too large. ²⁵ To solve this puzzle, L₁-norm is used as a relaxation of L₀-norm, ²⁶ which leads to the widely used L₁-norm regularization

$$\underset{\mathbf{f}}{\mathrm{minimize}}(\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_1)$$

(11)

where $\Vert \mathbf{f}{\Vert }_1=\sum _{i=1}^N|f_i|$ represents L₁-norm, which can induce the small values in the impact force vector f to become zero, thereby obtaining a sparse solution. Different from L₂-norm regularization, L₁-norm regularization has no analytical solution. However, Equation (11) is convex and thus convex optimization algorithms such as interior-point method, ²⁷ proximal gradient descent, ²⁸ and fast iterative shrinkage-thresholding algorithm (FISTA) ²⁹ can be used to obtain the global optimal solution. FISTA can be expressed as Algorithm 1, where $\mathrm{maxeig}\left({\mathbf{H}}^{\mathrm{T}}\mathbf{H}\right)$ represents the maximum eigenvalue of ${\mathbf{H}}^{\mathrm{T}}\mathbf{H}$ , and the soft-thresholding/shrinkage operator soft(·) for each element $x_i$ of the input vector x is defined as

$$\mathrm{soft}(x_i,T)=\mathrm{sign}(x_i)·\max (|x_i|-T,0)$$

(12)

where T is the threshold. Compared with the traditional L₂-norm regularization method, the L₁-norm regularization method can obtain sparse solutions. However, as the convex relaxation of the L₀-norm, L₁-norm imposes a weaker sparse constraint on the result than L₀-norm, resulting in defectively sparse force identification results with underestimated impact amplitudes.

Algorithm 1.

Fast iterative shrinkage-thresholding algorithm (FISTA) for L₁-norm regularization.

Set the step size $\alpha =1/L$ , and $L=\mathrm{maxeig}\left({\mathbf{H}}^{\mathrm{T}}\mathbf{H}\right)$ is the Lipschitz constant. input: $\lambda ,\mathbf{H},\mathbf{y}$ Step 0: ${\mathbf{f}}^{(0)}=0,{\mathbf{x}}^{(0)}=0,t_0=1,k=0.$ Step k: 1. ${\mathbf{g}}^{(k)}={\mathbf{x}}^{(k)}-\alpha ({\mathbf{H}}^{\mathrm{T}}(\mathbf{H}{\mathbf{x}}^{(k)}-\mathbf{y}))$ , 2. ${\mathbf{f}}^{(k+1)}=\mathrm{soft}({\mathbf{g}}^{(k)},\lambda \alpha )$ , 3. $t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}$ , 4. ${\mathbf{x}}^{(k+1)}={\mathbf{f}}^{(k+1)}+\frac{t_k-1}{t_{k+1}}({\mathbf{f}}^{(k+1)}-{\mathbf{f}}^{(k)})$ , 5. $k=k+1,$ 6. Repeat until convergence. output: ${\mathbf{f}}^{(k+1)}$

Fully overlapping group sparsity based on L_p-norm regularization

In this section, we describe the proposed method for impact force identification. Firstly, the non-convex L_p-norm is introduced for more accurate estimation of the amplitudes of impact forces. Then, the L_p-norm is combined with group sparsity, and thus more prior information is obtained to reduce the degrees of freedom of the solution in under-determined cases. Finally, AGSTA algorithm is proposed to solved the non-convex optimization problem. The overall scheme framework of the proposed impact force identification method is depicted in Figure 1.

Figure 1.

The framework of the proposed impact force identification method.

L_p-norm regularization

To address the amplitude-underestimation problem of L₁-norm regularization, the non-convex L_p-norm $(0<p<1)$ , which is more approximate to L₀-norm, is used to replace L₁-norm. Hence, Equation (11) can be rewritten as

$$\underset{\mathbf{f}}{\mathrm{minimize}}(\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda \Vert \mathbf{f}{\Vert }_p^p)$$

(13)

where ${\Vert \mathbf{f}\Vert }_p={\left(\sum _{i=1}^N|f_i{|}^p\right)}^{1/p}$ represents the L_p-norm $(0<p<1)$ .

In order to describe the difference between L_p-norm and L₁-norm more intuitively, we show the visual comparison results of different norm penalties in Figure 2, where $\mathrm{grad}(·)$ represents the gradient of different penalty terms. From Figure 2(a), we can see clearly that the L_p-norm penalty term matches L₀-norm better than L₁-norm. From the perspective of gradient descent, the L_p-norm penalty term has a larger gradient at the small value element, which is easier to cause iterating the small element to zero than the L₁-norm with a constant gradient value, resulting in a more sparse solution, as shown in Figure 2(b).

Figure 2.

Comparison of L₁-norm and L_p-norm penalty terms: (a) one-dimension example of L_p-norm corresponding to differentp values, and (b) gradient of L_p-norm and L₁-norm penalty terms.

L_p-norm regularization with fully overlapping sparsity

Since the impact force usually occurs in a brief time and it is a continuous process in the time domain, the impact force vector f naturally has the feature of group sparsity. That is, after grouping the elements in f, the elements in the group may all be zeros or non-zero values, as shown in Figure 3. Therefore, more satisfactory impact force identification results can be obtained at a group level instead of individual elements of the force vector. A general regularization method to exploit group sparsity is using the L _q ,1-norm $(q>1)$ as the regularization term, which can be expressed as

$$\underset{\mathbf{f}}{\mathrm{minimize}}(\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_{q,1})$$

(14)

where $\Vert ·{\Vert }_{q,1}$ is the L_q,1-norm mixed by the L₁-norm and the L_q-norm $(q>1)$ . For elements in each group, the L_q-norm is calculated first, and then all groups are combined by L₁-norm. When $q=2$ , Equation (14) is the famous L_2,1-norm regularization, ³⁰ that is,

$$\underset{\mathbf{f}}{\mathrm{minimize}}(\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_{2,1})$$

(15)

where $\Vert ·{\Vert }_{2,1}$ represents the L_2,1-norm in which the L₁-norm assembles the L₁-norms of different groups. For the impact force identification problem, the time of occurrence and duration of the impact force are unknown, so it is difficult for us to set the boundaries of the groups in advance. Hence, non-overlapping group sparsity, as shown in Figure 4(a), which requires all the groups to be individual, may lead to considerable errors in the results of impact force identification. In this paper, fully overlapping group sparsity is used to fix the defect of non-overlapping group sparsity, as shown in Figure 4(b). Then the L_2,1-norm in Equation (15) can be written as

$$\Vert \mathbf{f}{\Vert }_{2,1}=\sum _{i=1}^N\left({\Vert {\mathbf{f}}_{i,K}\Vert }_2\right)$$

(16)

where $\Vert {\mathbf{f}}_{i,K}\Vert =[f_i{f}_{i+1}\cdots f_{i+K-1}]$ is the $i\mathrm{th}$ group, N is the length of the force vector, and K is the fixed group size.

Figure 3.

The group sparse structure of the impact force in the time-space domain.

Figure 4.

Examples of group structure for a force vector of six elements: (a) the non-overlapping group structure with $K=3$ , and (b) the fully overlapping group structure with $K=3$ .

Compared with L₁-norm, L_p-norm is closer to the morphological features of L₀-norm, so it is expected to enhance the sparsity of the solution and improve the accuracy of the impact force identification results. Therefore, in this paper, the L₁-norm as the outer layer in L_2,1-norm regularizer in Equation (15) is replaced by the L_p-norm. Then we can obtain the FOGSL_p model for impact force identification problem

$$\underset{\mathbf{f}}{\mathrm{minimize}}\left(\underset{G\left(\mathbf{f}\right)}{\underbrace{\frac{1}{2}\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2+\lambda \Vert \mathbf{f}{\Vert }_{2,p}^p}}\right)$$

(17)

where $\Vert \mathbf{f}{\Vert }_{2,p}$ is a composite norm penalty that fuses the non-convex property of L_p-norm and the group sparsity of the force vector f, which is defined as

$$\Vert \mathbf{f}{\Vert }_{2,p}={\left(\sum _{i=1}^N(\Vert {\mathbf{f}}_{i,K}{\Vert }_2^p)\right)}^{1/p}$$

(18)

Equation (17) is a non-convex optimization problem. Fortunately, under the MM optimization framework, Equation (17) can be solved by constructing easy-to-optimize surrogate functions of non-convex terms. ³¹

Solution strategy for the FOGSL_p method

The key to MM algorithms is constructing an appropriate majorizer $M(\mathbf{f},{\mathbf{f}}_{\mathrm{k}})$ , which is the upper bounds of the optimization objective function $G({\mathbf{f}}_{\mathrm{k}})$ corresponding to the $k\mathrm{th}$ iterate ${\mathbf{f}}_{\mathrm{k}}$ . The majorizer is larger than the objective function $G(\mathbf{f})$ except for the tangent at ${\mathbf{f}}_{\mathrm{k}}$ , that is, $M({\mathbf{f}}_{\mathrm{k}},{\mathbf{f}}_{\mathrm{k}})=G({\mathbf{f}}_{\mathrm{k}})$ . When the majorizer is continuously minimized, the descending search of the objective function is completed. In this paper, majorizers of the residual term and the regularizer are both constructed to solve the problem.

Firstly, the majorizer of the residual term $G_1(\mathbf{f})=(1/2)\Vert \mathbf{y}-\mathbf{Hf}{\Vert }_2^2$ can be obtained by the second-order Taylor expansion, ³² that is,

$$\begin{array}{cc}\multicolumn{1}{c}{M_1^{(k)}\left(\mathbf{f}\right)=}& \frac{1}{2}\Vert \mathbf{y}-\mathbf{H}{\mathbf{f}}^{(k)}{\Vert }_2^2\hfill \\ \multicolumn{1}{c}{+{(\mathbf{f}-{\mathbf{f}}^{(k)})}^{\mathrm{T}}{\mathbf{H}}^{\mathrm{T}}(\mathbf{H}{\mathbf{f}}^{(k)}-\mathbf{y})}\\ \multicolumn{1}{c}{+\frac{L}{2}{(\mathbf{f}-{\mathbf{f}}^{(k)})}^{\mathrm{T}}\mathbf{I}(\mathbf{f}-{\mathbf{f}}^{(k)})}\\ \multicolumn{1}{c}{=\frac{L}{2}\Vert \mathbf{f}-({\mathbf{f}}^{(k)}-\frac{1}{L}{\mathbf{H}}^{\mathrm{T}}(\mathbf{y}-\mathbf{H}{\mathbf{f}}^{(k)})){\Vert }_2^2+C_1}\end{array}$$

(19)

where L is the Lipschitz constant that is greater or equal to the largest eigenvalue of ${\mathbf{H}}^{\mathrm{T}}\mathbf{H}$ , and $C_1$ is a constant independent of the vector f. $M_1\left(\mathbf{f}\right)$ meets the requirements of the majorizer, that is,

$$\begin{array}{c}{M}_1^{(k)}\left(\mathbf{f}\right)\ge G_1\left(\mathbf{f}\right)\forall \mathbf{f}\hfill \\ M_1^{(k)}\left({\mathbf{f}}^{(k)}\right)=G_1\left({\mathbf{f}}^{(k)}\right)\hfill \end{array}$$

(20)

For the regularizer $G_2(\mathbf{f})=\lambda \Vert \mathbf{f}{\Vert }_{2,p}^p$ , different from the residual term $G_1(\mathbf{f})$ , a nonquadratic majorizer based on L₁-norm is a natural choice. ³³ For a scalar variable f, the L_2,p-norm $(0<p<1)$ regularization term is $|f{|}^p$ , which satisfies the inequality when $f^{(k)}\ne 0$

$$|f{|}^p\le |f|p|f^{(k)}{|}^{p-1}+(1-p)|f^{(k)}{|}^p$$

(21)

which is an equality when $f=f^{(k)}$ . Notice that the majorizer in Equation (21) is undefined at $f^{(k)}=0$ . For any $f^{(k)}\in R$ , the majorizer can be defined as ³³

$$m(f;f^{(k)})=\{\begin{array}{cc}R\left(f^{(k)}\right)|f|+\beta ,\hfill & f^{(k)}\ne 0\hfill \\ \multicolumn{1}{c}{+\infty ,}& f^{(k)}=0\wedge f\ne 0\hfill \\ \multicolumn{1}{c}{0,}& f^{(k)}=0\wedge f=0\hfill \end{array}$$

(22)

where $R\left(f^{(k)}\right)=\lambda p|f^{(k)}{|}^{p-1}$ , and $\beta =(1-p)|f^{(k)}{|}^p$ is a constant independent of f. For force vector f, we can get the majorizer for the L_2,p-norm regularization term by setting $f=\Vert {\mathbf{f}}_{i,K}{\Vert }_2$ and $f^{(k)}=\Vert {\mathbf{f}}_{i,K}^{(k)}{\Vert }_2$ , that is,

$$M_2^{(k)}\left(\mathbf{f}\right)=\sum _{i=1}^{N}m({\Vert {\mathbf{f}}_{i,K}\Vert }_2;{\Vert {\mathbf{f}}_{i,K}^{(k)}\Vert }_2)$$

(23)

Similarly, the $M_2\left(\mathbf{f}\right)$ also satisfies the requirements of the majorizer. However, it is noticed that the individual elements in $\Vert {\mathbf{f}}_{i,K}{\Vert }_2$ are still inseparable, making subsequent iterations difficult. Fortunately, due to the sparse prior of the impact force vector f, we can employ the following inequality to decouple the elements in $\Vert {\mathbf{f}}_{i,K}{\Vert }_2$ , that is,

$$\Vert {\mathbf{f}}_{i,K}{\Vert }_2\le \Vert {\mathbf{f}}_{i,K}{\Vert }_1=\left|f_i\right|+\left|f_{i+1}\right|+\cdots +\left|f_{i+K-1}\right|$$

(24)

In order to further illustrate the rationality of choosing Equation (24) for decoupling, we show the $\Vert {\mathbf{f}}_{i,K}{\Vert }_2$ and $\Vert {\mathbf{f}}_{i,K}{\Vert }_1$ , respectively, in Figure 5 when f is a binary vector. One can see that the more sparse ${\mathbf{f}}_{i,K}$ is, the closer $\Vert {\mathbf{f}}_{i,K}{\Vert }_2$ and $\Vert {\mathbf{f}}_{i,K}{\Vert }_1$ are. Hence, taking $\Vert {\mathbf{f}}_{i,K}{\Vert }_1$ as the upper bound of $\Vert {\mathbf{f}}_{i,K}{\Vert }_2$ fits closely with the sparse priors of f. Finally, the surrogate function for the entire objective function Equation (17) can be obtained after combining the $M_1\left(\mathbf{f}\right)$ , $M_2\left(\mathbf{f}\right)$ , and Equation (24), that is,

$$M^{(k)}\left(\mathbf{f}\right)=\frac{L}{2}\Vert \mathbf{f}-{\mathbf{D}}^{(k)}{\Vert }_2^2+\sum _{i=1}^{N}m({\Vert {\mathbf{f}}_{i,K}\Vert }_1;{\Vert {\mathbf{f}}_{i,K}^{(k)}\Vert }_2)+C$$

(25)

where C is a constant from the superposition of the constants in $M_1^{(k)}\left(\mathbf{f}\right)$ and $M_2^{(k)}\left(\mathbf{f}\right)$ , which is still independent of f and irrelevant to the optimization process, and

$${\mathbf{D}}^{(k)}=({\mathbf{f}}^{(k)}-\frac{1}{L}{\mathbf{H}}^{\mathrm{T}}(\mathbf{y}-\mathbf{H}{\mathbf{f}}^{(k)}))$$

(26)

Figure 5.

The upper bound of the $\Vert {\mathbf{f}}_{i,k}{\Vert }_2$ of a two-dimensional vector $\mathbf{f}=[f_1{f}_2]$ .

The optimization problem of Equation (25), that is,

$${\mathbf{f}}^{(k+1)}=\underset{\mathbf{f}}{\mathrm{argmin}}{M}^{(k)}\left(\mathbf{f}\right)$$

(27)

is a weighted L₁-norm optimization problem, and the update rule with respect to each element $f_i$ is ³³

$$f_i^{(k+1)}=\mathrm{soft}(D_i^{(k)},\overline{R}\left(f_i^{(k)}\right))$$

(28)

where

$$\overline{R}\left(f_i^{(k)}\right)=\{\begin{array}{cc}+\infty ,\hfill & f_i^{(k)}=0\hfill \\ \multicolumn{1}{c}{R\left(f_i^{(k)}\right),}& f_i^{(k)}\ne 0\hfill \end{array}$$

(29)

Since $\mathrm{soft}(f,+\infty )=0$ for any $f\in R$ , if an element becomes zero in the optimization process, it will be kept as zero forever, leading to the desired sparse result. Considering that the impact force identification problem is usually ill-posed, and ${\mathbf{H}}^{\mathrm{T}}\mathbf{H}$ always has a large condition number, the convergence rate of Equation (27) may be slow. Therefore, we utilize the Nesterov’s acceleration technique here to improve the convergence rate, ³⁴ which is also used in FISTA. This technique uses ${\mathbf{f}}^{(k)}$ and ${\mathbf{f}}^{(k+1)}$ to obtain a new iteration starting point ${\mathbf{x}}^{(k+1)}$ to update, that is,

$${\mathbf{x}}^{(k+1)}={\mathbf{f}}^{(k+1)}+\frac{{\theta }_k-1}{{\theta }_{k+1}}({\mathbf{f}}^{(k+1)}-{\mathbf{f}}^{(k)})$$

(30)

where $\theta$ is a coefficient that changes with iteration. Since there are multiple local optimal solutions to Equation (17), the Nesterov’s acceleration process is improved by adding a judgment step to achieve a continuous minimization of the objective function. From the above derivation, an AGSTA can be summarized as Algorithm 2.

Algorithm 2

Accelerated grouped shrinkage-thresholding algorithm (AGSTA) algorithm.

Set the step size $\alpha =1/L$ , and $L=\mathrm{maxeig}({\mathbf{H}}^{\mathrm{T}}\mathbf{H})$ is the Lipschitz constant. input: $\lambda ,\mathbf{H},\mathbf{y},\mathrm{and}K$ Step 0: ${\mathbf{f}}^{(0)}=0,{\mathbf{x}}^{(0)}=0,{\theta }_0=1,k=0.$ Step k: 1. ${\mathbf{D}}^{(k)}={\mathbf{x}}^{(k)}-\alpha \left({\mathbf{H}}^{\mathrm{T}}\left(\mathbf{H}{\mathbf{x}}^{(k)}-\mathbf{y}\right)\right)$ , 2. $f_i^{(k+1)}=\mathrm{soft}(D_i^{(k)},\overline{R}\left(f_i^{(k)}\right))$ , 3. ${\theta }_{k+1}=\frac{1+\sqrt{1+4{\theta }_k^2}}{2}$ , 4. ${\mathbf{x}}^{(k+1)}=\min ({\mathbf{f}}^{(k+1)},{\mathbf{f}}^{(k+1)}+\frac{{\theta }_k-1}{{\theta }_{k+1}}({\mathbf{f}}^{(k+1)}-{\mathbf{f}}^{(k)}))$ , 5. $k=k+1,$ 6. Repeat until convergence. output: ${\mathbf{f}}^{(k+1)}$

One can see that when $K=1$ and $p=1$ , L_2,p-norm is an equivalence of L₁-norm, and Algorithm 2 degrades into FISTA for L₁-norm regularization. When $K=1$ and $0<p<1$ , L_2,p-norm is equivalent to the non-convex L_p-norm. Furthermore, with $p=1$ and $K>1$ , L_2,p-norm is the L_2,1-norm commonly used in group sparsity regularization. Thus, the L₁-norm ( $K=1$ , $p=1$ ) solved by FISTA, L_p-norm $(K=1)$ and L_2,1-norm $(p=1)$ regularization solved by AGSTA are employed in the following simulations and experiments sections for comparison with the FOGSL_p method.

Numerical study

In this section, numerical simulations are carried out to verify the performance and effectiveness of the proposed method on a composite panel with both ends fixed as depicted in Figure 6. The simulation results under different noise levels are compared with existing sparse regularization methods.

Figure 6.

The composite panel with both ends fixed.

In the area of structural health monitoring, the peak value of impact force is closely related to the damage of the structure. The time history of the impact force contains information about the process of impact occurrence. To quantitatively verify the performance of the proposed method, relative error (RE) and peak relative error (PRE) between the exact force ${\mathbf{f}}_i$ and the identified force ${\tilde{\mathbf{f}}}_i$ at the impact position i are defined as

$$\mathrm{RE}(\%)=\frac{{\Vert {\mathbf{f}}_i-{\tilde{\mathbf{f}}}_i\Vert }_2}{{\Vert {\mathbf{f}}_i\Vert }_2}\times 100$$

(31)

$$\mathrm{PRE}(\%)=\frac{{\Vert f_{\underset{i}{\mathrm{p}}}-{\tilde{f}}_{\underset{i}{\mathrm{p}}}\Vert }_2}{{\Vert f_{\underset{i}{\mathrm{p}}}\Vert }_2}\times 100$$

(32)

where $f_{\underset{i}{\mathrm{p}}}$ and ${\tilde{f}}_{\underset{i}{\mathrm{p}}}$ are the peak values of the exact impact force ${\mathbf{f}}_i$ and the identified one ${\tilde{\mathbf{f}}}_i$ at the impact position i. In addition, the index of localization accuracy (LA) for quantitative analysis of the localization of impact force is defined as

$$\mathrm{LA}(\%)=\frac{{\Vert {\tilde{\mathbf{f}}}_i\Vert }_2}{\sum _{j=1}^{n_f}{\Vert {\tilde{\mathbf{f}}}_j\Vert }_2}\times 100$$

(33)

To control the computation time within an acceptable range, Algorithms 1 and 2 need to set the termination criterion, which can be expressed as

$$\frac{{\Vert {\tilde{\mathbf{f}}}^{(k+1)}-{\tilde{\mathbf{f}}}^{(k)}\Vert }_1}{{\Vert {\tilde{\mathbf{f}}}^{(k)}\Vert }_1}\le \epsilon$$

(34)

where the value of ε is set to ${10}^{-5}$ .

Problem description

The size of the composite panel is 400 mm × 400 mm × 6 mm, and its left and right ends are both fixed. The material properties of the composite panel are listed in Table 1. As illustrated in Figure 6, the surface of the panel is divided into 64 regular grids with a size of 50 mm × 50 mm, and 25 potential impact points are located at the grid nodes, while r1–r24 are 24 candidate strain gauges. In the numerical simulations, the impact force acts randomly on one of the potential impact points. Responses of two strain gauges at r1 and r13 are selected in this section for force time history reconstruction and localization. The problem of impact force identification is severely under-determined in this condition.

Table 1.

Material properties of the composite panel.

Properties	Values
$E_1(\mathrm{GPa})$	138
$E_2(\mathrm{GPa})$	8.8
$G_{12}(\mathrm{GPa})$	4.5
${\upsilon }_{12}$	0.28
$\rho (\mathrm{kg}/{\mathrm{m}}^3)$	1560
Layups	[45/ −45/0/ −45/45/0/ −45/45/90/45/−45/45/ $\overline{0}$ ] s

The impact force is assumed to be a Gaussian distribution force, which has an analytical formula

$$f=A{e}^{-\frac{{(t-\mu )}^2}{2b^2}}$$

(35)

where A decides the amplitude of the impact force, $\mu$ determines the peak time of the impact force, and b controls the duration of the impact force. In this case, the parameters are set to $A=120\mathrm{N}$ , $\mu =3\mathrm{ms}$ , and $b=2.9\times {10}^{-4}$ to match the impact force profile in the experiments. Figure 7 illustrates the shape profile of the impact force in the time domain and the frequency domain. The spectrum of the impact force is concentrated below 5000 Hz, so the sampling frequency $f_s$ is set to 10,240 Hz according to the Nyquist-Shannon sampling theorem. The data length N is 256. The strain response without noise is calculated by transfer matrix H and impact force vector f. The transfer matrix H is obtained by the Green function method. ³⁵ To simulate the practical measurement environment, different levels of Gaussian white noise are added to the calculated strain response with signal-to-noise ratio (SNR) ranging from 5 dB to 30 dB with a step size of 5 dB:

$${\mathbf{y}}_n=\mathbf{y}+\sigma \mathbf{n}$$

(36)

where $\sigma$ is the standard deviation of noise, and n is a standard normally distributed vector generated by the Mersenne Twister algorithm. ³⁶ All programs are run on a processor Intel^® Core⁽™⁾ i5-10400F, 32 GB RAM.

Figure 7.

Synthesized impact excitation signal in (a) the time domain and (b) the frequency domain.

Number of the sensors

As stated above, this paper focuses on the situation where the number of sensors is less than the number of impact sources, that is, the under-determined case, as using as few sensors as possible for impact force localization and reconstruction is a challenging problem. And we try to realize impact force identification using a sparse and under-determined sensor configuration. Therefore, it is worth exploring how many sensors can ensure the accuracy of impact force identification. Monte Carlo simulations are performed with the number of sensors varying from 1 to 6. In this case, the impact force defined by Equation (35) is applied to point 6. The measurements are randomly chosen from 24 candidate sensors, and 100 sets of signals are selected for each number of sensors. Moreover, all the strain signals are corrupted by Gaussian white noise with $\mathrm{SNR}=15\mathrm{dB}$ . The sampling length is set to 256, so the size of the transfer matrix H varies from 256 × 6400 to 1536 × 6400 as the number of sensers changes.

Figure 8 depicts the average REs and PREs corresponding to different numbers of sensors. Notice that although the average PRE is below 5% when using only one sensor, the average RE is over 10%, which is more than twice that using two sensors. Therefore, in the subsequent cases, two sensors at r1 and r13 are used for impact force identification. The size of the corresponding transfer matrix is 512 × 6400, whose condition number tends to $+\infty$ , leading to the inverse problem of impact force identification being a severe ill-posed problem.

Figure 8.

The results of the Monte Carlo simulations while the impact force is applied to point 6: (a) the average REs of FOGSL_p using different numbers of sensors ranging from 1 to 6, and (b) the average PREs of FOGSL_p using different numbers of sensors ranging from 1 to 6.

Regularization parameter selection principle

In the proposed FOGSL_p method, there are three parameters to be selected, namely p, K, and λ. It has been revealed that the L_1/2-norm regularization is representative among the L_p-norm $(0<p<1)$ regularization methods in previous studies.^12,37–39 When $0.5\le p<1$ , the smaller the p is, the more sparse the solution obtained by the L_p-norm $(0<p<1)$ regularization; when $0<p\le 0.5$ , the performance of the L_p-norm regularization method is not significantly different. Therefore, the parameter p is set to 0.5 here. The K of the FOGSL_p and L_2,1 norm method is set to 15 to accommodate the duration of the impact force. As we know, the sparse regularizer can induce the sparse solution of the impact force identification problem, and the regularization parameter λ establishes a balance between the sparsity of the solution and the residual value of the solution and the exact value. With the λ getting smaller, the sparsity of the solution decreases, and a zero solution is obtained as the λ approaches infinity. Hence, the λ should be appropriately selected to obtain satisfactory impact force identification results. Various principles such as the L-curve method, ⁴⁰ Bayesian inversion method,^41,42 and hand-turning according to the experience ¹³ have been used to select the parameter λ. Here, the λ is selected over a wide range to obtain the optimal ratio of λ and ${\lambda }_{\max }$ corresponding to the minimum RE. The transfer matrix H changes when impact forces are applied to different locations, which alters the optimization objective function. Hence, we choose optimal values of λ for different impact force identifications in this section by minimizing RE. In this case, the impact force is applied at point 16, and the noise is added to the responses with $\mathrm{SNR}=15\mathrm{dB}$ . The same selection process for the optimal λ is cycled 50 times at the same noise level, and the results are averaged to eliminate the random error caused by noise. The average REs of FOGSL_p corresponding to different values of λ are shown in Figure 9(a), and the optimum $\lambda /{\lambda }_{\max }$ is 0.004 here. The values of λ of the other three methods, that is, L₁-norm regularization, L_p-norm regularization, and L_2,1-norm regularization methods, are selected by the same principle, and the results are shown in Figure 9(b), (c), and (d), respectively.

Figure 9.

The average REs corresponding to different values of $\lambda$ using (a) the FOGSL_p method, (b) the L₁-norm method, (c) the L_p-norm method, and (d) the L_2,1-norm method.

Impact force identification using the proposed method

In this section, impact force identification results of the proposed FOGSL_p method are compared with the other three regularization methods to validate the effectiveness of the FOGSL_p method. In this case, the impact force defined by Equation (35) is applied at point 18, chosen randomly from all 25 potential locations. The strain responses of r1 and r13 are collected simultaneously and mixed by noise with $\mathrm{SNR}=15\mathrm{dB}$ . Figure 10 shows the localization results of the FOGSL_p method, the L₁-norm method, the L_p-norm method, and the L_2,1-norm method. The proposed FOGSL_p method identifies the location of impact force from 25 potential impact points precisely with $\mathrm{LA}=100\%$ , and the identification results of other locations are zero. The L_p-norm method comes in second with $\mathrm{LA}=98.17\%$ and identifies the wrong impact force at point 5. The L₁-norm and L_2,1-norm methods perform unsatisfactorily, with $\mathrm{LA}=82.34\%\mathrm{and}85.10\%$ , respectively, and minor error impacts occur at most non-loading points. Figure 11 shows the reconstructed time histories of impact forces via the four methods. The reconstructed impact force via the FOGSL_p method is almost identical to the actual value, with $\mathrm{RE}=1.64\%$ and $\mathrm{PRE}=-0.79\%$ . Notice that although the L_p-norm method can decrease the PRE compared with the L₁-norm method, the reconstructed time history has a large deviation with RE of larger than 10% (15.19%), while the PRE and RE of the FOGSL_p method are both below 5%, improving identification accuracy significantly compared with the L₁-norm sparse regularization method.

Figure 10.

The localization results of the impact force at point 18: (a) the FOGSL_p method solution; (b) the L₁-norm method solution; (c) the L_p-norm method solution; and (d) the L_2,1-norm method solution.

Figure 11.

Impact force time history reconstruction results at point 18: (a) the FOGSL_p method solution; (b) the L₁-norm method solution; (c) the L_p-norm method solution; and (d) the L_2,1-norm method solution.

To analyze the computational efficiency of the AGSTA, the cost curves with different values of p and K are shown in Figure 12. In this case, the FOGSL_p method solved by AGSTA ( $p=0.5$ and $K=15$ ) converges to the set condition after 322 iterations, costing 1.12 s in total. The L₁-norm method solved by FISTA and L_2,1-norm method solved by AGSTA ( $p=1$ and $K=15$ ) spend more time (3.25 and 3.51 s) with 1644 and 1147 iterations, respectively, because of the constant gradient of the L₁-norm. The sudden drop in Figure 12(b) is due to the Nesterov’s acceleration strategy employed by FISTA to speed up the iterations by interpolation. It is noted that though the L_p-norm method converges as fast as FOGSL_p, the final identification result has much lower accuracy. The FOGSL_p method significantly improves the identification accuracy by taking similar calculation time to L_p-norm method.

Figure 12.

The cost function values with corresponding iterations using (a) the FOGSL_p method, (b) the L₁-norm method, (c) the L_p-norm method, and (d) the L_2,1-norm method.

In addition, to further verify the performance of the proposed method, the accuracy of the identification results at 12 randomly selected potential impact points is analyzed, as shown in Figure 13. One can see that the FOGSL_p method can accurately localize impact forces and reconstruct their time histories of impact forces acting on each point, and all of its REs and PREs remain below 5%. Moreover, the LAs of the FOGSL_p method are 100% at all positions, indicating its stability in localizing impact force. Table 2 indicates the respective iteration numbers and computation time of the FOGSL_p method and the other three methods. The L_p-norm method sovled by AGSTA ( $K=1$ and $p=0.5$ ) requires the least iterations in most cases, but has a lower reconstruction accuracy for impact force time history, which results from the fast gradient descent rate of the L_p-norm penalty, but due to its non-convexity can only converge to a local optimal solution. As a comparison, the L_2,1-norm method solved by AGSTA ( $K=15$ and $p=1$ ) takes more iterations and more computation time, but it is more capable of reconstructing the impact force time history. The FOGSL_p method retains the ability of L_p-norm to promote sparse solutions, and at the same time incorporates more prior information through group sparsity, so it can significantly improve the impact force reconstruction accuracy while improving computational efficiency.

Figure 13.

Impact force identification results at different positions with a 15 dB of noise: (a) REs of the FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods; (b) PREs of the FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods, and (c) LAs of the FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods.

Table 2.

Computation time and iterations for identifing impact forces at different positions with a 15 dB noise level.

Impact location		21	23	16	3	7	24	20	17	1	22	10	5
Iterations	FOGSLp	740	233	321	379	437	291	355	1097	637	434	376	366
	L 1-norm	1777	1017	1657	1040	1388	1339	1805	1235	1494	1379	1265	1551
	Lp-norm	252	238	775	176	367	317	831	321	256	1061	632	552
	L 2,1-norm	1242	993	712	777	413	782	1476	886	916	974	889	1293
Computation time (s)	FOGSLp	2.35	0.73	1.03	1.22	1.39	0.95	1.14	3.47	2.01	1.38	1.23	1.18
	L 1-norm	3.49	2.10	3.32	2.11	2.71	2.59	3.50	2.41	2.91	2.68	2.46	3.02
	Lp-norm	0.78	0.70	2.34	0.52	1.12	0.95	2.46	0.98	0.75	3.17	1.87	1.66
	L 2,1-norm	3.62	2.98	2.14	2.31	1.24	2.30	4.36	2.63	2.70	2.86	2.64	3.82

The effect of the noise level

To inspect the robustness to noise levels of the proposed FOGSL_p method, different noise levels, that is, 30, 25, 20, 15, 10, and 5 dB, are considered here. The impact force is applied to point nine randomly here. Considering that the parameter $\lambda$ is bound up with the optimization objective function, values of $\lambda$ are selected according to the minimum RE principle at each noise level. The same program at each noise level is run one hundred times and the results are averaged to eliminate random errors. One can see that the FOGSL_p method accurately localizes the impact force at different noise levels and maintains its LA at $100\%$ for all cases, demonstrating more consistent performance than the other three methods. As shown in Figure 14, the REs and PREs of all the four methods increase with the decrease of SNR. One can note that the RE and PRE of the FOGSL_p method remain below $5\%$ even at 5 dB noise level, standing in sharp contrast to the notable increase of the other methods. These results suggest that the proposed FOGSL_p method possesses excellent robustness and great practical application potential.

Figure 14.

Impact force identification results under different noise levels: (a) REs of FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods, (b) PREs of FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods, and (c) LAs of FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods.

Experimental validation

In this section, impact force identification experiments are conducted on a composite panel with both ends fixed with similar properties in the simulations. The experiment setup is shown in Figure 15. Here, two strain gauges are also used to monitor 25 potential impact points, keeping the same transfer matrix dimension as in the simulation study. Once more, the performance and feasibility of the FOGSL_p method are compared with the L₁-norm method, L_p-norm method, and L_2,1-norm method.

Problem description

The layups of the composite panel employed in the impact force identification experiments are consistent with those in the simulation study, but the dimension slightly differs from that in the simulations, which is 400 × 300 × 6 mm as shown in Figure 15. The dynamic behavior of the compiste panel is considered linear, as in the simulations. The surface of the panel is divided into 48 regular grids of 50 × 50 mm, of which 25 nodes are set as potential impact points, and nine strain gauges are placed at r1–r9, respectively, to measure the responses of the impact. The impact force is applied by an impact hammer (PCB® 086C03; PCB Piezotronics, Inc., Buffalo, NY, USA, with sensitivity 2.25 mV/N), and a PZT force sensor integrated into the hammer head records the time history of the applied impact force for reference. Referring to the simulation study, here we select the strain responses of r4 and r5 (with sensitivity 51.1 and 47.5 mV/µε, respectively) to identify the impact force. The impact force and strain signals are acquired simultaneously with a sampling frequency of 10,240 Hz using the Leuven Measurement System (LMS) SCADAS III data acquisition system. The transfer matrix H is obtained in the following ways. First, the modal test module in the LMS system is employed to obtain the FRF $h_{i,j}(\omega )$ between the impact point j and the response point i. Then the IRF $h_{i,j}(t)$ in the time domain is obtained by inverse fast fourier transform. And finally, the $h_{i,j}(t)$ is arranged into a Toeplitz form matrix H after discretization. Consistent with the simulation study, the data length is 256 sampling points, and thus the dimension of the transfer matrix H is 512 × 6400 (Figure 15).

Figure 15.

The set-up of the impact force identification experiment on the composite panel.

Results and discussion

In this case, the impact force is imposed at point 23, and the strain responses are measured by strain gauges at r4 and r5. Impact modal tests are conducted to obtain the transfer matrices, and the impact forces and corresponding strain signals are shown in Figure 16. Then different impact forces are imposed at point 23 for validation purposes, as depicted in Figures 17 and 20. Since the actual value of the impact force is unknown in practical applications, the selection principle for parameter $\lambda$ by minimizing RE adopted in the simulations is no longer suitable here. Hence, here we use the L-curve criterion to select the $\lambda$ . The parameter $\lambda$ corresponding to the maximum curvature of the L-curve is chosen as the optimal $\lambda$ .

Figure 16.

The impact forces and corresponding strain signals used to obtain the transfer matrices between point 23 and r4, and between point 23 and r5.

Figure 17.

The shape of the impact force excited by the hammer at point 23 in (a) the time domain and (b) the frequency domain.

Figure 18 shows the localization results of the impact force on the composite panel using the FOGSL_p method and the other three sparse regularization methods from two strain responses. One can see that the FOGSL_p method accurately localizes the impact force acting at point 23 with $\mathrm{LA}=100\%$ and $\mathrm{PRE}=-1.05\%$ . In contrast, the L₁-norm method and the L_2,1-norm method perform relatively poorly, identifying small erroneous force values at other points where no impact force is applied (with $\mathrm{LA}=82.29\%\mathrm{and}83.53\%$ , respectively), while underestimating the amplitude of the impact force, with PRE of $-2.65\%$ and $-3.83\%$ , respectively. The L_p-norm method performs the worst, identifying relatively large impact forces at three wrong points, with PRE exceeding $20\%$ , which can be regarded as a failure of identification.

Figure 18.

The localization results of impact force excited by the hammer at point 23: (a) the FOGSL_p method solution; (b) the L₁-norm method solution; (c) the L_p-norm method solution; and (d) the L_2,1-norm method solution.

Figure 19 shows the reconstructed time history of the impact force. It can be seen that two impact forces are applied in an extremely brief time in the experiment. As demonstrated in the magnified time history figure, the FOGSL_p method accurately reconstructs the time histories of two consecutively applied impact forces with the RE of $7.3\%$ over the entire analysis period. The time history of the impact force reconstructed by the L_p-norm method has a significant deviation from the actual value with the second impact being unidentified, and the RE reaches $50.19\%$ , which is unacceptable. The reconstruction results of another four impact forces excited at point 23 with different durations and magnitudes are shown in Figure 20. It is shown that FOGSL_p can accurately reconstruct the time histories of different impact forces, and the REs are all below $10\%$ , while the PREs are controlled to be below $5\%$ .

Figure 19.

The time history reconstruction results of the impact force excited by the hammer at point 23: (a) the FOGSL_p method solution; (b) the L₁-norm method solution; (c) the L_p-norm method solution; and (d) the L_2,1-norm method solution.

Figure 20.

The reconstruction results of different impacts excited at point 23 using the FOGSL_p method: (a)–(d) four impact forces with different durations and magnitudes.

To further verify the computational efficiency of the proposed method, the convergence rates of different methods are analyzed, as shown in Figure 21. In this case, the FOGSL_p method converges to the set accuracy after 834 iterations, taking a total of 2.64 s. In contrast, the L₁-norm and L_2,1-norm methods converge after more iterations, and the computation time is more than twice that of the FOGSL_p method, costing 5.64 and 5.06 s, respectively. Although the L_p-norm method also converges fast, the final identification result is not satisfying, which is due to the fact that it only obtains the local optimal solution and cannot guarantee the global optimal solution. Figure 22 shows the corresponding REs and PREs of the FOGSL_p method and the other three methods at 12 randomly selected locations. One can see that the FOGSL_p method can keep its PREs below $10\%$ in most cases, while its most REs does not exceed $20\%$ . In contrast, the other three methods have lower identification accuracy and fail to identify the impact forces at several points. The PRE and RE of L_p-norm method solution at point 15 even exceed 60%. Moreover, the FOGSL_p method accurately localizes the impact forces in most cases, maintaining its most LA at 100%. The above impact force identification results indicate that the proposed FOGSL_p method can accurately localize the impact force and reconstruct its time history with high accuracy under highly under-determined conditions.

Figure 21.

The computation time and iterations of FOGSL_p and the other three methods while the impact force is excited at point 23.

Figure 22.

Impact force identification results at different positions: (a) REs of FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods; (b) PREs of FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods; and (c) LAs of FOGSL_p, L₁-norm, L_p-norm, and L_2,1-norm methods.

Conclusions

To sum up, a novel FOGSL_p method is proposed for the inverse problem to localize and reconstruct impact forces simultaneously. More critical piors of impact forces are applied to penalty function by combining the promotion of sparsity of the non-convex L_p-norm and the fully overlapping group sparsity, and therefore stronger constraints are introduced to stabilize the solution in under-determined case. Then an AGSTA algorithm is proposed to solve the non-convex optimization problem under the MM framework, which constructs a surrogate function to optimize the non-convex and element-coupled objective function. Moreover, modified Nesterov’s acceleration strategy is adopted to increase the convergence rate. The performance of the proposed method is verified from the perspective of simulations and experiments in under-determined case by comparing with L₁-norm, L_p-norm, and L_2,1-norm regularization methods. The RE-minimization principle and the L-curve criterion are applied for the regularization parameter selection in simulations and experiments, respectively. Both simulations and experiments indicate that the FOGSL_p method can obtain satisfactory results by using two sensors while monitoring 25 potential positions, which is a severely under-determined condition, and the FOGSL_p method always outperforms the three other methods with lower REs and PREs and higher LAs. Furthermore, in simulations, the FOGSL_p method consistently performs best under different noise levels, whose PREs and REs remain below 5% even at 5 dB noise level, which demonstrates the excellent robustness of the FOGSL_p method. Finally, the experimental results illustrate the effectiveness of the proposed FOGSL_p method in practical applications.