M Aucejo

O De Smet

J.-F Deü

On a space-time regularization for force reconstruction problems

Keywords: Inverse problem, Force reconstruction, Space-time approach

Tikhonov and LASSO regularizations are commonly used to solve force reconstruction problems in time domain. Unfortunately, these particular forms of additive regularization are not well adapted to tackle both localization and time reconstruction problems simultaneously, since they are generally restricted to the reconstruction of sources sharing the same space and time characteristics. To alleviate this limitation, a multiplicative space-time regularization is introduced. The proposed regularization strategy takes advantage of one's prior knowledge of the space-time characteristics of excitation sources. It also introduces a novel reconstruction model based on the generalized-α method, which is unconditionally stable and second-order accurate. The validity of the proposed method is assessed numerically and experimentally. In particular, comparisons with standard regularization terms point out the practical benefit in exploiting both spatial and temporal prior information simultaneously in terms of quality and robustness of reconstructed solutions.

Multiplicative regularization, Generalized-α method

Introduction

The resolution of force reconstruction problems in time domain still remains an active topic in the mechanical engineering community as suggested by the intensive dedicated literature. Except particular strategies [START_REF] Nordström | A dynamic programming algorithm for input estimation on linear time-variant systems[END_REF][START_REF] Wu | Statistical moving load identification including uncertainty[END_REF][START_REF] Yan | Impact load identification of composite structure using genetic algorithms[END_REF][START_REF] Jones | Detection of impact location and magnitude for isotropic plates using neural networks[END_REF][START_REF] Doyle | A wavelet deconvolution method for impact force identification[END_REF][START_REF] Qiao | A force identification method using cubic b-spline scaling functions[END_REF][START_REF] Qiao | The application of cubic b-spline collocation method in impact force identification[END_REF][START_REF] O'donoughue | Time-resolved identification of mechanical loadings on plates using the virtual fields method and deflectometry measurements[END_REF], the most widespread approaches to solve this problem nowadays are certainly Kalman filtering [START_REF] Gillijns | Unbiased minimum-variance input and state estimation for linear discrete-time systems[END_REF][START_REF] Lourens | An augmented Kalman filter for force identification in structural dynamics[END_REF][START_REF] Azam | A dual Kalman filter approach for state estimation via output-only acceleration measurements[END_REF][START_REF] Petersen | Estimation of the full-field dynamic response of a floating bridge using Kalman-type filtering algorithms[END_REF] and regularization techniques [START_REF] Jacquelin | Force reconstruction: analysis and regularization of a deconvolution[END_REF][START_REF] Mao | Experimental study of hammer impact identification on a steel cantilever beam[END_REF][START_REF] Li | A load identification method based on wavelet multi-resolution analysis[END_REF][START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF]. However, both methods operate in different context. Kalman filtering is an online strategy making assumptions on the evolution of the systems [START_REF] Haug | Bayesian Estimation and Tracking -A practical guide[END_REF][START_REF] Sarkka | Bayesian filtering and smoothing[END_REF]. On the contrary, regularization is an offline approach that allows exploiting one's prior knowledge on the sources to identify through the definition of the regularization term [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Chambolle | Total variation minimization and a class of binary mrf models[END_REF][START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF][START_REF] Li | A hierarchical bayesian method for vibration-based time domain force reconstruction problems[END_REF]. The latter characteristic is at core of the present paper and explains why the next of this introduction is only focused on the analysis of existing regularization strategies.

In general, two categories of reconstruction problem can arise in practical situations. The first one is related to the localization of excitation sources, while the second one consists in reconstructing the time signal of prelocalized sources. Regarding the localization problem, the regularization term reflects the spatial prior information available on the sources to reconstruct.

It is often expressed as the q -norm of the solution sought. Such a norm is flexible enough to express one's prior knowledge on the forces to identify, since smooth solutions are obtained for q = 2 [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF][START_REF] Boyd | Convex optimization[END_REF], while localized ex-citation fields are promoted for q ≤ 1 [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Grasmair | Non-convex sparse regularization[END_REF][START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF]. On the other hand, for reconstructing the excitation signal of prelocalized sources, the regularization term has to reflect the associated prior information. It should be noted that the choice of a particular value of q is strongly dependent on the length on the time window on which the reconstruction is performed. Indeed, a signal can be considered continuous over a short duration and localized on the long term. Generally, when the force signal exhibits a certain continuity, the regularization term is defined from the 2 -norm of the solution vector to identify and leads to the Tikhonov regularization [START_REF] Jacquelin | Force reconstruction: analysis and regularization of a deconvolution[END_REF][START_REF] Mao | Experimental study of hammer impact identification on a steel cantilever beam[END_REF][START_REF] Turco | Tools for the numerical solution of inverse problems in structural mechanics: review and research perspectives[END_REF]. On the contrary, when the excitation signal is rather impulsive, the regularization term is constructed from the 1 -norm and gives rise to the LASSO regularization [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Qiao | Sparse deconvolution for the large-scale ill-posed inverse problem of impact force reconstruction[END_REF]. Consequently, it appears that most of the methods proposed in the literature are theoretically not always well adapted to tackle both the localization and time reconstruction problems at the same time, except for configurations where the force vector to identify has the desired structure or the spatial distribution of the sources and the nature of the excitation signals share the same space-time characteristics 1 such as the sparsity as shown in Refs. [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF]. To the best of our knowledge, only a few methods have been developed to address these issues [START_REF] Sun | A new improved regularization method for dynamic load identification[END_REF][START_REF] Liu | A novel computational inverse technique for load identification using the shape function method of moving least square fitting[END_REF][START_REF] Liu | Dynamic load identification for stochastic structures based on gegenbauer polynomial approximation and regularization method[END_REF][START_REF] Liu | An efficient method to reduce ill-posedness for structural dynamic load identification[END_REF]. However, these methods generally addressed the space-time reconstruction problem in a separated manner [START_REF] Turco | A strategy to identify excitating forces acting on structures[END_REF][START_REF] Li | A novel approach for distributed dynamic load reconstruction by space-time decoupling[END_REF][START_REF] Qiao | A novel method for force identification based on the discrete cosine transform[END_REF][START_REF] Kalhori | Inverse estimation of impact on a composite panel using a single piezoelectric sensor[END_REF]]. 1 This is typically the case of a hammer impact excitation reconstructed over a long duration, for which the corresponding spatial distribution is sparse (point force) as well as the time excitation signal (impulsive excitation). Another example is a distributed harmonic excitation.

It is thus of primary interest to simultaneously exploit both the spatial and the temporal features of excitation sources to constrain the space of admissible solutions to aid the reconstruction process in finding an optimal solution. These requirements are actually met by regularization terms derived from mixed p,q -norms. Indeed, mixed 2,q -norm for q ≤ 1 has revealed all its potential in signal and image recovery applications [START_REF] Kowalski | Sparse regression using mixed norms[END_REF][START_REF] Zheng | Subspace weighted 2,1 minimization for sparse signal recovery[END_REF][START_REF] Wang | On recovery of block-sparse signals via mixed 2 / q (0 < q ≤ 1) norm minimization[END_REF]. In the context of force reconstruction, Rezayat et al. first introduced an additive regularization using a regularization term based on a mixed 2,1 -norm to identify broadband point forces in the frequency domain [START_REF] Rezayat | Identification of dynamic forces using group-sparsity in frequency domain[END_REF]. Recently, Aucejo and De Smet have extended this idea by developing, for frequency domain applications, a multiplicative regularization based on the definition of a mixed p,q -norm to accurately reflect experimenters knowledge on the type (localized or distributed) of the excitation forces, as well as on the nature of the force spectrum [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF].

In the present paper, a space-time regularization is proposed for time domain applications in order to solve both localization and time reconstruction problems within a unique framework. Actually, the proposed method relies on three pillars. The first one is related to the definition of the reconstruction model. Here, the reconstruction model is obtained from a state-space model of the structure built from a generalized-α integration scheme [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method[END_REF]. In doing so, the reconstruction model is unconditionally stable and second-order accurate. The second key feature is the space-time regularization term defined, as in [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF], from a mixed p,q -norm to properly exploit the space-time infor-mation available on the sources to identify. The last pillar of the proposed approach concerns the formulation of the inverse problem, since it is based on the multiplicative regularization recently introduced by the authors in the context of mechanical source identification [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF][START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF][START_REF] Aucejo | Further investigation on "A multiplicative regularization for for reconstruction[END_REF]. This particular strategy is generally computationally more efficient than the corresponding additive regularization, because it is free from the preliminary definition of any regularization parameter [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]. This explains its use in the present paper. To clearly introduce the main features of the proposed regularization strategy, this article is divided into four parts. Before considering the core of the paper, the need for another regularization strategy for dealing with time domain applications is explained in section 2. Section 3 is devoted to the introduction of the three pillars of the space-time regularization, as well as that of the related resolution algorithm. Numerical and experimental validations of the space-time regularization are proposed in section 4 and 5. Obtained results point out the practical interest in exploiting both spatial and temporal prior information simultaneously in terms of quality and robustness of reconstructed solutions.

The need for another regularization method in time domain

As highlighted in the introduction, the most widespread regularization strategies to solve the source identification problem in time domain are the Tikhonov regularization [START_REF] Jacquelin | Force reconstruction: analysis and regularization of a deconvolution[END_REF][START_REF] Mao | Experimental study of hammer impact identification on a steel cantilever beam[END_REF][START_REF] Li | A load identification method based on wavelet multi-resolution analysis[END_REF][START_REF] Kazemi Amiri | Derivation of a new parametric impulse response matrix utilized for nodal wind load identification by response measurement[END_REF][START_REF] Zemcik | Reconstruction of impact on composite airfoil segment using piezoelectric sensors[END_REF] and the LASSO regularization [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Qiao | Impact-force sparse reconstruction from highly incomplete and inaccurate measurements[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF][START_REF] Qiao | Sparse deconvolution for the large-scale ill-posed inverse problem of impact force reconstruction[END_REF], which belong to the class of additive regularization methods.

The related multiplicative regularizations can be formally written under the following generic form:

F = argmin F\{0} Y -HF 2 2 • F q q , ( 1
)
where H is the convolution matrix, Y is the measured output vector, F is the force vector and q is the norm parameter. This formulation gives rise to the multiplicative Tikhonov regularization (mTIK) when q = 2 and to the multiplicative LASSO regularization (mLASSO) when q = 1.

Unfortunately, these formulations can lead to inaccurate reconstructions when one wants to identify both the location and the time signal of the excitation field. Indeed, the norm parameter q helps to exploit one's prior knowledge on the sources to identify, since choosing q = 2 promotes distributed solutions [START_REF] Boyd | Convex optimization[END_REF][START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF], while setting q ≤ 1 enforces the sparsity of the solution vector [START_REF] Chartrand | Nonconvex regularization for image segmentation[END_REF][START_REF] Nikolova | Efficient reconstruction of piecewise constant images using nonsmooth nonconvex minimization[END_REF]. In other words, mTIK and mLASSO regularizations are generally limited to the reconstruction of sources sharing the same space and time characteristics such as the sparsity (e.g. reconstruction of a hammer impact).

To illustrate this, let us consider the reconstruction of a harmonic point force imposed by a shaker on the simply-supported 1m-long beam described in section 4 at x 0 = 0.3 m from its left end [see section 4 for details]. To implement the overall reconstruction process, it is first assumed that the bending displacement field Y, measured in 9 points with a Digital Image Correlation system for instance, is corrupted by an additive Gaussian white noise such that the corresponding signal-to-noise ratio (SNR) is equal to 25 dB. Then, the convolution matrix H is computed using the procedure described in section 3.1. Finally, the identification of the excitation field is performed by supposing that all the reconstruction points (i.e. the points where the excitation field is identified) are collocated with the measurement points, including the excitation point [see the application of mTIK and mLASSO regularizations. A thorough analysis of these results show that none of the considered regularization strategy is able to provide a satisfying reconstruction in the present application. Indeed, mTIK regularization allows identifying the shape of the excitation signal, but not its amplitude, which is underestimated. Furthermore, in absence of information on the number of actual sources and their location, the inspection of the waterfall representation can lead to erroneous conclusions, since the excitation field is more distributed than expected. On the other hand, mLASSO regularization correctly retrieves the point source location, but the reconstructed excitation signal is clearly different from the reference one. A better reconstruction can be obtained by choosing the norm parameter q ∈ [START_REF] Nordström | A dynamic programming algorithm for input estimation on linear time-variant systems[END_REF][START_REF] Wu | Statistical moving load identification including uncertainty[END_REF] to achieve a compromise between the spatial sparsity and the time continuity of the excitation field [see Fig. 4]. However, this choice is very challenging when one has only a rough idea of the actual excitation field. Consequently, the proposed example shows the importance of exploiting both spatial and temporal prior information on the sources to identify for improving the quality of the reconstructed solutions.

Theoretical background

Regularization consists in including in the formulation of the inverse problem some prior information on the excitation field to identify in order to constrain the space of solutions. In the present paper, the regularization problem is expressed as:

F = argmin F F(Y -HF) • R(F), (2)
where

• F(Y -HF) is the data-fidelity term which controls the a priori on the noise corrupting the data [START_REF] Green | Statistics of Images, the TV Algorithm of Rudin-Osher-Fatemi for Image Denoising and an Improved Denoising Algorithm[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF][START_REF] Le | A variational approach to reconstructing images corrupted by poisson noise[END_REF]. Here, it is assumed that the reconstruction model is linear and defined such that:

Y = HF; (3)
• R(F) is the regularization term that encodes prior information on the excitation field F, i.e. what is known about the excitation field before making any measurement [START_REF] Chambolle | Total variation minimization and a class of binary mrf models[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Hintermuller | Nonconvex TV q -models in image restoration: Analysis and a trust-region regularization based superlineraly convergent solver[END_REF].

At this stage, it is clear that the quality of the reconstructed excitation field is not only conditioned to the adequacy of the data-fidelity and regularization terms with the actual noise and the actual space-time source characteristics, but also to the quality of the reconstruction model describing the dynamic behavior of the considered structure. For all these reasons, a particular attention is paid to the definition of each element of the regularization problem. That is also why, the present paper introduces an original reconstruction model, as well as an adapted regularization term for encoding available space-time prior information.

Derivation of the reconstruction model

The reconstruction model given by Eq. ( 3) can be classically established either from the discretization of the Duhamel's integral [START_REF] Jacquelin | Force reconstruction: analysis and regularization of a deconvolution[END_REF][START_REF] Kazemi Amiri | Derivation of a new parametric impulse response matrix utilized for nodal wind load identification by response measurement[END_REF][START_REF] Qiao | Sparse deconvolution for the large-scale ill-posed inverse problem of impact force reconstruction[END_REF][START_REF] Qiao | Sparse regularization for force identification using dictionaries[END_REF][START_REF] Law | Moving loads -Dynamic analysis and Identification techniques[END_REF] or from the discretization of the state-space representation of the mechanical system [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Law | Moving loads -Dynamic analysis and Identification techniques[END_REF][START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF][START_REF] Bernal | Sequential deconvolution input reconstruction[END_REF]. In the present paper, the reconstruction model is built from a discretized state-space representation, because we are prone to think that such an approach can deal with measurement data of various types (displacement, velocity, acceleration, strain, etc.) in a simple, efficient and convenient manner.

Discretized state-space representation of the dynamical system

The most widespread approach to discretize the state-space representation of a dynamical system is the zero-order-hold (ZOH) sampling technique, which consists in assuming that the input of the system f(t) is constant within a time step. Formally, the corresponding discretized state-space representation, composed of the state and output equations, is written:

     x k+1 = A x k + B f k y k = O x k + D f k , (4)
where x k , f k and y k are the state, input and output vectors at sample k, while A, B, O and D are, respectively, the discretized state, input, output and feedthrough matrices.

Although widely used, this model is conditionally stable and its performances strongly depends on the sampling rate. That is why, Liu et al. proposed a state-space model established from the Newmark integration scheme [START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF]. However, despite Newmark approach is unconditionally stable for a well-chosen set of parameters (namely γ = 1/2 and β = 1/4) [START_REF] Géradin | Mechanical Vibrations -Theory and Application to Structural Dynamics[END_REF], it does not allow introducing numerical damping on the response of the high frequency modes present in the model. Consequently, if those modes are not properly discretized in space and time, numerical errors can affect the computed mechanical fields [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method[END_REF]. A practical way of introducing numerical damping in the Newmark method consists in choosing the parameters γ and β such that [START_REF] Géradin | Mechanical Vibrations -Theory and Application to Structural Dynamics[END_REF]:

γ = 1 2 + α and β = 1 4 1 + α 2 with 0 ≤ α < 1. (5)
As appealing it seems, the resulting Newmark scheme is first-order accurate only. To preserve the second-order accuracy, some other numerical integration methods have been introduced. One of the most general approaches is the generalized-α method [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method[END_REF], which includes the Newmark method, the HHT (Hilber-Hughes-Taylor) scheme [START_REF] Hughes | Computational methods for transient analysis, chapter Analysis of transient algorithms with particular reference to stability behavior[END_REF] and the WBZ (Wood-Bossak-Zienkiewicz) approach [START_REF] Wood | An Alpha Modification of Newmark's Method[END_REF]. For a proper choice of parameters, the generalized-α method is unconditionally stable and second-order accurate and minimizes the numerical damping at low frequencies for a given highfrequency damping level. That is why, we propose to extend the approach introduced by Liu et al. in Ref. [START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF] to derive a state-space model based on the generalized-α method.

Derivation of the state equation

Let us consider a regular partition of the time domain such that

t 0 < • • • < t k < • • • < t f (
t 0 and t f : initial and final instants) and let h = t k+1 -t k denote the time step size. To derive the state equation from the generalized-α method for a mechanical system described by its mass, stiffness and damping matrices (M, K, C), one has to start from the following equations [START_REF] Erlicher | The analysis of the generalized-α method for non-linear dynamic problems[END_REF]:

M .. d k+1-αm + C . d k+1-α f + K d k+1-α f = S f f k+1-α f , (6a)
.

d k+1 = . d k + (1 -γ) h .. d k + γ h .. d k+1 , (6b)
d k+1 = d k + h . d k + 1 2 -β h 2 .. d k + β h 2 .. d k+1 , (6c)
where γ and β are the parameters of the Newmark method, α m and α f are two averaging parameters associated to the inertia and internal/external forces, while d k = d(t k ) and f k = f(t k ) are, respectively, the displacement and excitation vectors at time t k , S f being the selection matrix of the excitation degrees of freedom. In the previous equation, the quantity d k+1-α must be read as:

d k+1-α = (1 -α) d k+1 + α d k . (7)
For linear systems, it is a common practice to use a modally reduced order model to limit the computational costs. In this regard, when the displacement vector is expanded on its modal basis, the previous equation of motion becomes:

..

q k+1-αm + Z n . q k+1-α f + Ω 2 n q k+1-α f = Φ T n S f f k+1-α f , (8a)
.

q k+1 = . q k + (1 -γ) h .. q k + γ h .. q k+1 , (8b)
q k+1 = q k + h . q k + 1 2 -β h 2 .. q k + β h 2 .. q k+1 , (8c)
where q k is the generalized displacement, Φ n is the matrix of the massnormalized mode shapes,

Ω n = diag(ω 1 , . . . , ω n ) and Z n = diag(2ξ 1 ω 1 , . . . , 2ξ n ω n ),
where ω n and ξ n are respectively the modal angular frequency and the modal damping ratio of the mode n.

The first step to derive the state equation consists in expressing the acceleration and velocity vectors at instant t k+1 from the generalized acceleration, velocity and displacement vectors at instant t k and the generalized displacement vector at instant t k+1 [START_REF] Liu | Explicit form of an implicit method for inverse force identification[END_REF]. After some manipulations of Eqs. (8b) et (8c), it comes:

.

q k+1 = γ β h q k+1 -q k + 1 - γ β . q k + 1 - γ 2β h .. q k , (9a)
..

q k+1 = 1 β h 2 q k+1 -q k - 1 β h . q k - 1 2β -1 .. q k . ( 9b
)
By introducing the last two relations into Eq. (8a), one gets:

q k+1 = L d Φ T n S f [(1 -α f ) f k+1 + α f f k ] + K d q k + C d . q k + M d .. q k , (10)
where

L d = 1 -α m β h 2 I + (1 -α f ) γ β h Z n + (1 -α f ) Ω 2 n -1
,

K d = L d 1 -α m β h 2 I + (1 -α f ) γ β h Z n -α f Ω 2 n , C d = L d 1 -α m β h I + (1 -α f ) γ β -1 Z n -α f Z n , M d = L d (1 -α m ) 1 2β -1 I + (1 -α f ) γ 2β -1 h Z n -α m I .
Replacing now Eq.( 10) in Eq. (9a), one finds:

.

q k+1 = L v Φ T n S f [(1 -α f ) f k+1 + α f f k ] + K v q k + C v . q k + M v .. q k , (12)
where

L v = γ β h L d , K v = γ β h [K d -I] , C v = γ β h C d + 1 -γ β I and M v = γ β h M d + 1 -γ 2β h I.
Similarly, after introducing Eq.( 10) in Eq. (9b), it readily comes:

..

q k+1 = L a Φ T n S f [(1 -α f ) f k+1 + α f f k ] + K a q k + C a . q k + M a .. q k , (13)
where

L a = 1 β h 2 L d , K a = 1 β h 2 [K d -I] , C a = 1 β h 2 [C d -h I] and M a = 1 β h 2 M d - 1 2β -1 I.
All things considered, if the state vector x k is defined such that 10), ( 12) and ( 13) can be gathered to form the desired state equation:

x k = [q T k . q T k .. q T k ] T , then Eqs. (
x k+1 = A x k + B + f k+1 + B -f k , (14)
where

A =      K d C d M d K v C v M v K a C a M a      , B + = (1 -α f )      L d L v L a      S f , B -= α f      L d L v L a      S f .
Finally, it should be noted that the proposed state equation based on the generalized-α method is unconditionally stable and second-order accurate provided that [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method[END_REF][START_REF] Géradin | Mechanical Vibrations -Theory and Application to Structural Dynamics[END_REF]:

α f = ρ ∞ 1 + ρ ∞ , α m = 3 α f -1, γ = 1 2 + α f -α m and β = 1 4 (1 + α f -α m ) 2 , ( 15
)
where ρ ∞ is the spectral radius belonging to the interval [0, 1]. In the present paper, ρ ∞ is set to 1.

Derivation of the output equation

Considering the previously defined state vector, the output equation can be expressed as:

y k = S d d k + S v . d k + S a .. d k = S d Φ n q k + S v Φ n . q k + S a Φ n .. q k = O x k , (16)
where S d , S v and S a are, respectively, the selection matrices of the displacement, velocity and acceleration data measured on the structure and

O = [S d , S v , S a ]Φ n .

State-space model

From the state and output equations presented above, one can form the state-space representation of the considered dynamical system derived from the generalized-α method, namely:

     x k+1 = A x k + B + f k+1 + B -f k y k = O x k . (17)

Reconstruction model

The implementation of a regularization strategy requires the definition of an ad-hoc reconstruction model expressed as a linear system of equations relating the excitation field F to identify to the measured data Y through the convolution matrix H [see Eq. ( 3)]. To derive such a relation from the state-space model given by Eq. ( 17), it is first necessary to express the state vector at instant t k from the state and input vectors at instants t j ≤ t k . In doing so, it can be shown that the state vector x k writes:

x k = A k x 0 + A k-1 B -f 0 + k-1 j=1 A k-j-1 B f j + B + f k , for t k > t 0 , (18)
where B = A B + + B -and x 0 et f 0 are the known initial state and input vectors.

As a result, the output vector y k is given, for t 0 < t k ≤ t n (n: number of time samples), by the relation:

y k = OA k x 0 + OA k-1 B -f 0 + k-1 j=1 OA k-j-1 B f j + OB + f k . (19)
In the end, the reconstruction model used in this paper simply writes:

Y = H F, (20)
where

Y =            y 1 y 2 y 3 . . . y n            -            OA OA 2 OA 3 . . . OA n            x 0 -            OB - OAB - OA 2 B - . . . OA n-1 B -            f 0 , H =         OB + 0 • • • 0 OB OB + . . . . . . . . . . . . . . . 0 OA n-2 B OA n-3 B • • • OB +         , F =            f 1 f 2 f 3 . . . f n            .
Here, the initial state and input vectors are chosen so that x 0 = 0 and f 0 = 0. Consequently, when measuring output data, it is compulsory to start the recordings before exciting the structure. A final word must be added about the dimensions of the proposed reconstruction model, because it can be thought that the size of the reconstruction problem should be necessarily large. Actually, in practical applications, the number of measured responses is always limited as well as the number of reconstruction points. In other words, this means that the output vector Y and the excitation vector F are not evaluated on all the degrees of freedom of the structure, but only on a (possibly small) subset2 . To fix ideas, if the number of measurement points is N m and the number of reconstruction points is N r (N m not necessarily equal to N r ), then the dimensions of Y and F are respectively equal to

(N m • n) × 1 and (N r • n) × 1.

Space-time regularization

The core idea behind the space-time regularization is to exploit available prior information on the type of the sources (localized or distributed) and the nature of the excitation signal simultaneously, while properly exploiting prior information on the noise corrupting the data. As already evoked in the preamble of this section, all this information is encoded in the data-fidelity and regularization terms. That is why, a particular attention is paid to their definition, which is detailed in the following subsections.

Definition of the regularization term

An appropriate definition of the regularization term is all the more crucial since it conditions the quality of the reconstructed solution [START_REF] Ginsberg | New approach for impact detection by finding sparse solution[END_REF][START_REF] Boyd | Convex optimization[END_REF][START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF].

In this paper, the regularization term is defined, as for the space-frequency regularization developed by the authors, from the general mixed p,q -norm in order to exploit one's prior knowledge on the spatial and temporal characteristics of the forces to reconstruct [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF].

Formally, the proposed space-time regularization term is expressed as:

R(F) = F q p,q . (22)
To better understand the influence of the space-time regularization term during the resolution process, it is interesting to represent the unknown force vector F as a matrix, where the rows correspond to the time signal at a particular location and the columns to the excitation field at a specific instant.

In doing so, it comes:

F = [f 1 . . . f j . . . f n ] =            f 11 • • • f 1j • • • f 1n . . . . . . . . . f i1 • • • f ij • • • f in . . . . . . . . . f m1 • • • f mj • • • f mn            , ( 23
)
where m is the number of reconstruction points of the excitation field.

Consequently, by recalling that the mixed p,q -norm is defined by the relation:

F p,q =   m i=1 n j=1 |f ij | p q p   1 q , ∀ (p, q) ∈ R * + , (24)
it is clear that the proposed regularization term introduces an explicit coupling between the coefficients of F and allows promoting some structures observed in real signals [START_REF] Gramfort | Mixed-norm estimates for the m/eeg inverse problem using accelerated gradient methods[END_REF]. To illustrate this particular property of the mixed norms, let us consider the case for which (p, q) = (2, 1). In this situation, the matrix F is supposed to be sparse along the lines (space) and full along the rows (time). In other words, one promotes the spatial sparsity of the excitation field and the continuity of the time signals. Consequently, the proposed space-time regularization term is highly flexible, since it allows dealing with various force distributions and various excitation signals within a unique framework. However, it can be argued that in real-life applications finding proper values for p and q is far from an easy task, because the force distribution is unknown or there is no idea about the force to recover. We are prone to think that even in these situations rough information about the force to recover is available. This information is provided by a careful analysis of the mechanical system and recorded time signals. Once useful information is obtained, one has to translate it into mathematical terms through the choice of the norm parameters p and q. To this end, some guidelines, based on our own experience [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF][START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Aucejo | Bayesian source identification using local priors[END_REF][START_REF] Aucejo | Bayesian formulations for force reconstruction problems[END_REF] and the existing literature on convex and non-convex optimization [START_REF] Boyd | Convex optimization[END_REF][START_REF] Chartrand | Nonconvex regularization for image segmentation[END_REF][START_REF] Jain | Non-convex optimization for machine learning[END_REF], can be used. Regarding the spatial distribution, if the excitation field is supposed to be distributed, then it is reasonable to set q = 2. On the contrary, if the excitation field is rather localized, then choosing q ≤ 1 allows promoting the sparsity of the solution.

Regarding now the mathematical transcription of one's prior knowledge of excitation signal, it is reasonable to set p = 2 if the time signal is supposed to be continuous over the considered duration, while p ≤ 1 if the signal is rather impulsive.

Definition of the data-fidelity term

The data-fidelity term is a measure of the difference between the measured output data Y and the reconstructed output data HF. Accordingly, it reflects prior information on the bias between measured and reconstructed output data. In other words, it defines one's prior knowledge on the noise corrupting the data. A common assumption consists in considering that the measured data are corrupted by an additive Gaussian white noise, leading to a datafidelity term of the form:

F(Y -HF) = Y -HF 2 2 . ( 25
)

Generic form of the space-time regularization

The generic form of the proposed space-time (ST) regularization is simply obtained by introducing Eqs. ( 22) and (25) into Eq. ( 2). In doing so, one obtains:

F = argmin F\{0} Y -HF 2 2 • F q p,q . (26)

Resolution algorithm

The solution of the proposed space-time regularization can only be found in an iterative manner, insofar as the minimization problem has no explicit solutions. That is why, the regularized solution is computed here from an adapted version of the Iteratively Reweighted Least Squares (IRLS) algorithm [START_REF] Rodriguez | An iteratively weighted norm algorithm for total variation regularization[END_REF].

The basic idea behind the IRLS algorithm is to replace the direct resolution of the minimization problem by an equivalent iterative process having a unique solution at each iteration. To this end, the mixed p,q -norm is replaced by a weighted 2 -norm, that is:

∀(p, q), F q p,q = W(F) 1/2 F 2 2 , (27)
where W(F) is a weighting matrix depending explicitly on the solution vector F.

As part of an iterative scheme, the aim is to find the solution F

), in order to satisfy Eq. ( 27) when the iterative process has converged. Thereby, F

at iteration r + 1 is solution of the following minimization problem:

F (r+1) = argmin F\{0} Y -H F 2 2 • W (r) 1/2 F 2 2 = H T H + α (r+1) W (r) -1 H T Y, (28)
where α (r+1) is the adaptive regularization parameter, defined such that [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF][START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF]:

α (r+1) = Y -H F (r) 2 2 W (r) 1/2 F (r) 2 2 . ( 29
)
Practically, the coefficients of the weighting matrix W (r) are defined such that [START_REF] Aucejo | A space-frequency multiplicative regularization for force reconstruction problems[END_REF]:

W (r) I = W s(r) i • W t(r) i,j , (30)
where I = j + n(i -1) is a global index and W s(r) i and W t(r) i,j are the weighting coefficients associated to the space (s) and the time (t) domains. By construction, the weighting coefficients are written:

W t(r) i,j = max p-2 , fij p-2 (31)
and

W s(r) i = max 2 q p -1 , W t(r) i 1/2 F[i] (r) 2 2 q p -1 , ( 32
)
where

F[i] = [ fi1 . . . fij . . . fin ] and W t(r) i = diag W t(r) i,1 , . . . , W t(r) i,n
.

In the previous relations, is a damping parameter that allows avoiding infinite weights. The damping parameter is automatically selected once for all at the beginning of the iterative process from the cumulative histogram of the initial solution | F | are less than or equal to .

Finally, to completely implement the IRLS algorithm, it remains to define an initial solution and a reliable stopping criterion. In the present paper, the initial solution is chosen as the solution of the standard Tikhonov regularization, namely:

F (0) = argmin F Y -H F 2 2 + α (0) F 2 2 = H T H + α (0) I -1 H T Y, ( 33
)
where I is the identity matrix of dimension m • n and α (0) is the initial adaptive regularization parameter chosen as in Ref [START_REF] Aucejo | A multiplicative regularization for force reconstruction[END_REF].

Regarding the choice of stopping criterion, the iterative process is stopped when the relative variation δ of the adaptive regularization between two successive iterations is less than or equal to some tolerance. Formally, the relative variation δ is defined by:

δ = α (k) -α (k-1) α (k-1) . ( 34
)
Experimentally, it has been found that setting the tolerance to 10 -8 allows obtaining consistent reconstructions.

Numerical validation

This numerical validation intends to investigate the practical interest of applying the proposed approach for solving force reconstruction problems.

Description of the numerical test case

In this section, we consider the space-time reconstruction of a harmonic point force of amplitude 15 N oscillating at 100 Hz imposed by a shaker from the measurement of the displacement field only. The harmonic excitation signal f (t) imposed on the structure is defined such that:

f (t) = A sin(2π f 0 t), (35)
with A = 15 N and f 0 = 100 Hz.

The studied structure is a simply-supported IPN steel beam of length 1 m, cross-sectional area 1060 mm 2 and second moment of area 171 mm 4 . The coordinate of the point force, measured from the left end of the beam, is

x 0 = 0.3 m.
In the present numerical validation, it is supposed that only displacement data are available. To synthesize these experimental data, the equation of motion of the beam is first expanded on the corresponding modal basis, which is computed analytically. Then, the resulting modally reduced system, containing the first 18 bending modes (resonance frequencies less than 1 kHz), is solved using the Newmark method for (γ, β) = (1/2, 1/4) and h = 0.1 ms (sampling frequency: 10 kHz). Here, a modal damping ratio equal to 1% for all the modes has been imposed. In a last step, the obtained output fields is corrupted by an additive Gaussian white noise with a SNR equal to 25 dB.

Regarding the reconstruction model, the output matrix O should be defined so that the output vector y k corresponds to the measured displacement field. As a result, the output matrix is expressed as:

O = S d Φ n 0 0 , (36)
where S d is the selection matrix of the displacement data measured on the structure.

Finally, it should be noted that the displacement field is computed on nine measurement positions equally distributed along the beam between 0.1 m and 0.9 m. Then, the reconstruction of the excitation fields is performed by supposing that all the measurement points are collocated with the identification points, including the excitation points, as presented in Fig. 1.

Application

In section 2, it has been shown that standard regularization terms don't provide satisfying results, because they don't fully exploit space-time prior information on the excitation field to reconstruct. This is the reason why the space-time regularization is developed in this paper.

To apply the space-time regularization, the mixed norm parameters must be properly chosen. For this purpose, the test case must be carefully analyzed. Here, the excitation field is necessarily spatially sparse, insofar as a point force is to the structure by a shaker. On the other hand, mechanical signals are continuous by nature. It is thus reasonable to promote this feature in the formulation of the reconstruction problem. This information can also be inferred after the analysis of the output signals. From the previous considerations, the norm parameters defining the regularization term are chosen such that (p, q) = (2, 0.5).

As shown in Figure 5, the proposed ST regularization allows reconstructing more consistently the target excitation field than mTIK and mLASSO regularization approaches. By extension, it indicates that the standard regularizations classically used in the literature dedicated to time domain source identification are not always the best ones. Actually, the choice of a particular approach must be done after a thorough study of the experimental set-up.

Nevertheless, ST regularization is a safe choice when one wants to localize and quantify the excitation sources during a certain period of time. As a side note, it should be noted that displacement data has been used in the proposed numerical validation to reconstruct the unknown excitation field.

For the sake of completness, reconstruction results obtain from acceleration data are presented in Appendix A.

Experimental validation

To confirm the conclusions drawn in the previous numerical experiment, the space-time is applied to a real structure. The main objective of this experimental validation is to demonstrate the advantages of the proposed strategy in terms of identification quality and robustness in operating conditions. In this section, the reconstruction of a hammer impact is considered.

Description of the experimental set-up

The structure under test is a thin aluminum plate of 60 cm in length, 40 cm in width and 6 mm in thickness, clamped along its length in a wooden support. The effective width of the plate resulting from the mounting conditions is 39.1 cm. To perform all the subsequent measurements the system is suspended to a rigid structure through a set of elastic bungee cords [see Fig. 6] and all the measurements (FRF and time signals) are recorded and processed using an M+P data acquisition system.

O = 0 0 S a Φ n , (37)
where S a is the selection matrix of the acceleration data measured on the structure.

To establish the reconstruction model from a modal reduction, the modes of the system must be known beforehand. For this purpose, an experimental modal analysis (EMA) has been conducted from a roving hammer test performed on a grid of 17× 17 points using three reference accelerometers and a impact hammer equipped with a steel tip to properly excite all the modes below 6500 Hz. The goal of the EMA is to extract the modes from the real structure in order to limit the influence of modeling errors when establishing the reconstruction model. The locations and the associated identification numbers of the references accelerometers with respect to the measurement grid are presented in Fig. 7.

Once all the measurements collected, the modal parameters have been extracted from a MDOF curve fitting algorithm. In total, the first 83 flexible modes and 2 suspension modes (at 1 Hz and 2.5 Hz) have been extracted.

As an illustration, the modal parameters (resonance frequency and modal damping) of the first 10 flexible modes extracted by EMA are given in Table 1. Furthermore, a global overview of these modal parameters is presented in Figure 8.

Measurement of acceleration data and excitation signal

The output data resulting from a hammer impact have been collected using four accelerometers mounted on the structure. The input and output measurement devices has been located at nodes of the grid defined for EMA to ensure the consistency of the reconstruction process. The locations and the identification numbers (ID) of the measured input and output data are defined in Fig. 10. In particular, it can be seen that no accelerometer is located at the excitation point (ID: 236). This means that the reconstruction have been performed in non-collocated configuration.

From a practical point of view, the hammer is equipped with a soft rubber tip, so as to fix the cut-off frequency of the excitation around 500 Hz. In this way, the convergence of the reconstruction model in terms of modal series is ensured. Regarding the signal processing parameters, the sampling frequency has been set to 16384 Hz for 32768 lines. In other words, the data acquisition time is fixed to 2 s, while the length of the time step is equal to 61 µs. To As indicated previously, because no accelerometer is mounted at the excitation point, the reconstruction of the excitation signal is necessarily performed in non-collocated configuration. For the sake of clarity, the identification numbers of the measurement and reconstruction points are given in Table 2.

Application

To compare ST regularization with mTIK and mLASSO regularizations, the reconstruction is first performed on a sequence of 24 ms defined between 96 ms and 120 ms from the acceleration data measured on this interval. As shown in the zoomed portion of the excitation signal presented in Fig. 11, it is reasonable to promote the continuity of the reconstructed excitation signals. On the other hand, the excitation field is necessarily sparse, since the structure is excited by a hammer impact. Consequently, the norm parameters defining the space-time regularization term are set to (p, q) = (2, 0.5).

When applying mTIK (q = 2), mLASSO (q = 1) and ST regularizations on the selected sequence, we observe that all the three regularization approaches give very similar reconstructed excitation fields, which are in very good agreement with the target excitation field [see Fig. 13]. To quantify the accuracy of the reconstructed solution with respect to the regularization strategy, the correlation and the relative error (RE) are evaluated and given in Table 3. Formally, the correlation is given by:

C = F T meas • F id F meas 2 • F id 2 , (38)
where F meas is the force measured by the force sensor and F id is the force identified by mTIK, mLASSO or ST regularizations. Regarding now, the relative error, it is defined such that: However, the results significantly differ when considering a shorter reconstruction sequence. Indeed, when the reconstructions are performed on a sequence of 5.4 ms defined between 96 ms and 101.4 ms, mTIK, mLASSO and ST regularizations behave rather differently, as shown in Fig. 14

E = F meas -F id 2 2 F meas 2 2 . ( 39
)

Conclusion

The present paper has introduced a space-time multiplicative regularization to reconstruct mechanical sources in time domain. The proposed regularization method relies on three particular features. First, a novel reconstruction model based on a discretized state-space representation of the mechanical system, derived from a generalized-α integration scheme, has been developed. This allows obtaining an unconditionally stable and second-order

Fig. 1 Figure 1 :

 11 Figure 1: Locations of the measurement and reconstruction points along the beam -(•) Ends of the beam, ( ) Measurement points and (•) Reconstruction points

Figure 2 :

 2 Figure 2: Reconstruction of the excitation field from mTIK regularization -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

Figure 3 :Figure 4 :

 34 Figure 3: Reconstruction of the excitation field from mLASSO regularization -(a) Waterfall representation and (b) Identified time signal at excitation point -(-) Reference signal and (--) Reconstructed signal

  at iteration r by setting W (r) = W( F

  Actually, its value is calculated so that 5% of the values of | F (0)

Figure 5 :

 5 Figure 5: Reconstruction of the excitation field from the space-time regularization for (p, q) = (2, 0.5) -(a) Waterfall representation and (b) Identified time signal at excitation -(-) Reference signal and (--) Reconstructed signal

Figure 6 :

 6 Figure 6: Experimental set-up

Finally, the mode shapes of flexible modes 5 Figure 7 :Figure 8 :

 578 Figure 7: Locations of the reference accelerometers with respect to the grid used for EMA -(•) Roving hammer excitation points and ( ) Measurement points (reference accelerometers)

Figure 9 :

 9 Figure 9: Mode shapes of three identified flexible modes -Mode 5: 295.13 Hz, (b) Mode 30: 1824.20 Hz and (c) Mode 60: 3986.80 Hz

Figure 10 :Figure 11 :Figure 12 : 5 . 1 . 3 .

 101112513 Figure 10: Location of the accelerometers and the hammer impact with respect to grid used for EMA -(•) Nodes of the grid used for EMA, ( ) Accelerometers et (♦) Hammer impact

Figure 13 :

 13 Figure 13: Reconstruction of the excitation field corresponding to a hammer impact over 24 ms -(a) Waterfall representation -mTIK regularization, (b) Waterfall representation -mLASSO regularization, (c) Waterfall representation -ST regularization for (p, q) = (2, 0.5) and (d) Identified time signal at excitation point -(-) Reference signal and (--) Tikhonov regularization, (-• -) LASSO regularization and (• • • ) ST regularization

Figure 14 :

 14 Figure 14: Reconstruction of the excitation field corresponding to a hammer impact over 5.4 ms -(a) Waterfall representation -Tikhonov regularization, (b) Waterfall representation -LASSO regularization, (c) Waterfall representation -ST regularization for (p, q) = (2, 0.5) and (d) Identified time signal at excitation point -(-) Reference signal and (--) Tikhonov regularization, (-• -) LASSO regularization and (• • • ) ST regularization

Figure A. 15 :

 15 Figure A.15: Reconstruction of the excitation field corresponding to harmonic point force -(a) Waterfall representation -mTIK regularization, (b) Waterfall representation -mLASSO regularization, (c) Waterfall representation -ST regularization for (p, q) = (2, 0.5) and (d) Identified time signal at excitation point -(-) Reference signal and (--) Tikhonov regularization, (-• -) LASSO regularization and (• • • ) ST regularization

Table 1 :

 1 Resonance frequency and modal damping of the first 10 flexible modes

	ID Frequency (Hz) Modal damping (%)
	1	65.26	0.74
	2	122.61	0.55
	3	141.86	0.43
	4	221.02	0.36
	5	295.13	0.31
	6	315.81	0.30
	7	489.74	0.37
	8	501.56	0.18
	9	537.30	0.20
	10	540.96	0.19

Table 2 :

 2 Identification numbers of measurement and reconstruction points with respect to grid used for EMA

	Measurement Reconstruction
	19	19
	66	66
	228	228
	153	236

Table 3 :

 3 Correlation and RE for mTIK, mLASSO and ST regularizations

	Reg.strategy Correlation RE (%)
	mTIK	0.9963	0.74
	mLASSO	0.9997	0.06
	ST	0.9993	0.15
	This surprising result could potentially be explained by the fact that the
	forces are directly related to the accelerations through the second Newton's
	law, making the reconstruction easier.

Table 4 .

 4 and Obtained results show that, in this particular situation, only ST

	regularization allows obtaining a consistent space-time reconstruction of the
	target excitation field. Consequently, this experimental validation indicates
	that properly exploiting simultaneously the space-time characteristics of exci-
	tation sources is beneficial in terms of quality and robustness of reconstructed
	solutions.

Table 4 :

 4 Correlation and RE for mTIK, mLASSO and ST regularizations

	Reg. strategy Correlation RE (%)
	mTIK	0.6267	67.4
	mLASSO	0.7703	52.5
	ST	0.9999	0.01

If the actual excitation source lies outside the chosen subset, the proposed approach, like any other regularization methods, provides a set of pseudo-forces (minimizing the considered objective functional) and not the actual excitation field.

Acknowledgements

The work presented in this paper has been performed within the research collaboration framework agreement CT-413-1601-CNES "Experimental reconstruction of transient mechanical sources" with the Centre National d'Étude Spatiales (Cnes). The authors are grateful to the Cnes for funding this research.

Appendix A. Numerical validation based on acceleration data

In section 4, displacement data has been used to reconstruct the unknown excitation field. In the engineering practice, acceleration data are preferred because they can be measured accurately at relatively low cost. For this reason, acceleration-based reconstructions are presented in this appendix. This