# Kalman-based estimation of loading conditions from ultrasonic guided wave measurements

Andre Dalmora<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Alexandre Imperiale<fnref n="1" />, Sébastien Imperiale<fnref n="2" />,<fnref n="3" />, Philippe Moireau<fnref n="2" />,<fnref n="3" />

E-mail: philippe.moireau@inria.fr

March 2024

**Abstract.** Ultrasonic guided wave-based Structural Health Monitoring (SHM) of structures can be perturbed by Environmental and Operations Conditions (EOCs) that alter wave propagation. In this work, we present an estimation procedure to reconstruct an EOC-free baseline of the structure from the only available Ultrasonic guided wave measurements. This procedure could typically be used as a prior step to increase the robustness of a more general ultrasonic imaging algorithm or SHM process dedicated to flaw detection. Our approach is model-based, i.e. we use a precise modeling of the wave propagation altered by structure loading conditions. This model is coupled with the acquired data through a data assimilation procedure to estimate the deformation caused by the unknown loading conditions. From a methodological point of view, our approach is original since we have proposed an iterated Reduced-Order Unscented Kalman strategy, which we justify as an alternative to a Levenberg-Marquardt strategy for minimizing the non quadratic least-squares estimation criteria. Therefore, from a data assimilation perspective, we provide a quasi-sequential strategy that can valuably replace more classical variational approaches. Indeed, our resulting algorithm proves to be computationally very effective, allowing us to successfully apply our strategy to realistic 3D industrial SHM configurations.

## 1. Introduction

In various cutting-edge industrial fields, e.g. nuclear power generation, transportation, or aeronautics, the safe and reliable use of critical parts of structures is of paramount importance. To meet safety regulations in these areas, one must often be able to assess the integrity of the materials or equipment that make up these critical parts. To this end, numerous Non-Destructive Techniques (NDT) have been developed over the years. They are means of examining the material in question and obtaining quantitative information about its integrity without damaging it. Among them, Structural Health Monitoring (SHM) is an approach that – compared to other NDT techniques – generates a continuous stream of field data by incorporating actuators and sensors in situ. In other

<footnote n="1">¹Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France</footnote>

<footnote n="2">²Inria, Team-MΞDISIM, Inria-Saclay Ile de France 91120 Palaiseau, France</footnote>

<footnote n="3">³LMS, Ecole Polytechnique, CNRS – Institut Polytechnique de Paris, France</footnote>
words, SHM systems monitor structures as they are used. One way to implement such
systems is to rely on ultrasonic Guided Waves (GWs) [28, 38] because of their attractive
properties, such as propagation over long distances to study a large volume of material,
or sensitivity to local thickness variations due to dispersion phenomena. Nevertheless,
there are a number of challenges in the actual implementation of GW-based SHM
systems in realistic structures. To name just two of them that drive the goals of our work,
let us mention that the monitoring system should be minimal to avoid overloading the
structure of interest, and that the Environmental and Operations Conditions (EOCs)
have a non-negligible impact on GW propagation [21]. Typical EOCs that can affect
wave propagation include temperature variations and the mechanical loading sustained
by the structure during its use. We focus on the latter because it is more general in
terms of mathematical formulation and can indeed lead to a change in wave velocities or
even induce anisotropy. These effects are often referred to by the term “acoustoelastic”
propagation [40, 1]. These acoustoelastic effects can alter the data registered by
ultrasonic sensors, potentially affecting the precision of GW-based monitoring systems.
In this context, the main objective of our work is the following: using the available data
– the ultrasonic measurements – we aim to remove the bias caused by the mechanical
loading conditions in order to reconstruct an EOC-free baseline. This objective solves
the above two problems since the influence of mechanical loading is captured without
additional sensors. This EOC-free baseline can then be used to find evidence of potential
defects or damage within the structure in the ultrasonic signals – assuming here that
the contributions of the defects and the EOCs are separable in the time or frequency
domain.

In essence, our goal is to reconstruct a pre-deformation of the structure using only
the GW measurements. From the wave propagation point of view, this is an inverse
problem that we solve as a minimization problem of a fidelity-to-data functional under
the condition that wave propagation dynamics is satisfied. An important feature arising
from our EOC context is that the propagation model is in fact the one obtained by
linearizing the nonlinear elastodynamics model around the sought pre-deformation, as
presented in details in previous works [12] and recalled later in this article.

Solving this non-linear optimization problem can be carried out through different
approaches. A first one is to resort to gradient descent iterations or quasi-Newton
processes [26, 11]. At each iteration, the gradient of the cost function can be obtained
by solving the so-called forward and backward adjoint problems. This approach is
arguably the most common one to address this type of problem, it is referred to as
the variational approach (4D-Var) [14, 8] in the data assimilation community with now
extension to mechanical systems [22], while in the geophysics community, it has been
labeled as Full Waveform Inversion (FWI) [44] with adaptation to tomography [7]. A
strong advantage of FWI is its robustness w.r.t. the size of the parametric space, i.e.
the space in which lies the (discrete) solution of the minimization problem. In fact,
numerous successful applications of this method have led to the reconstruction of wave
velocity maps over large propagation domains. However, one significant difficulty of
this approach is managing the adjoint dynamics, which contains the tangent of the
propagation model around the forward trajectory. In the context of our work, this
tangent model is intricate. In particular, it entails the third derivative of the hyperelastic
potential ruling the constitutive behavior of the material. Also, storing the forward
trajectory to evaluate this tangent is prohibitive, since the state space is very large
in the context of high-frequency time-domain wave propagation. Note that however,
in the case of inviscid wave propagation, one can save storage space by simply back-
propagating the forward trajectory – a technique exploited in other wave propagation
inverse problems [37, 13]. However, this is done at the cost of yet another call to the
wave propagation solver.

A second approach is to resort to sequential methods and, in particular, to Kalman
filtering approaches that can be developed in a stochastic or a deterministic context
[5, 41]. The term sequential stems from the fact that the main building blocks of
this approach are (exclusively) forward problems, where the dynamics are modified by
the addition of a feedback loop. This feedback loop is proportional to the discrepancy
between the synthetic data generated from the model's current trajectory and the actual
measurements. The resulting modified dynamics is often referred to as a sequential
estimator (in the stochastic context) or an observer (in the deterministic context) of
the target trajectory. In Kalman-based filtering, the gain operator in the feedback
loop is computed from a covariance operator satisfying a Riccati equation [3]. For
Linear-Quadratic (LQ) problems, i.e. linear state dynamics with linear parameter-state
coupling and a quadratic misfit functional, the observer at the final time corresponds
exactly to the solution of the minimization problem [30]. Also, in the special case
where the uncertainty is limited exclusively to the parameter space, this method
leads to a specific application of the Reduced-Order Kalman Filter [32, 30]. This
method can be generalized to cases with non-linear parameter-state coupling, leading
to the so-called Reduced Order Extended Kalman Filter (EKF) or its gradient-free
version, the Reduced-Order Unscented Kalman Filter (UKF) [32, 31]. One of the main
advantages of these sequential estimators lies in their ability to provide a solution to
the minimization problem in one pass, consisting of embarrassingly parallel forward
problems. Moreover, thanks to their tangent-free alternative, e.g. UKF, they are
easy to interface with legacy code used in parallel as a black box, propagating each
UKF sigma-point, also referred to as particle. Nevertheless, the usual computational
bottleneck of these sequential methods is to store and invert the (dense) covariance
matrix whose dimension corresponds to the size of the parameter space. This limits
the application of Kalman filters to relatively small parameter spaces compared to the
typical configurations handled by variational or FWI methods. Moreover, although
Kalman filters are exactly equivalent to the minimization problem LQ problems their
extension to nonlinear minimization problems is either prohibitively expensive [29, 30]
or approximate when relying on Extended or Unscented Kalman filters [32, 31, 30].

In our work, we develop a new estimation procedure that combines elements of the
two approaches. Namely, we first apply a Levenberg-Marquardt (LM) algorithm [23, 4]
to derive from the initial minimization problem a set of LQ subproblems satisfied by parameter increments. Each of these LQ subproblems is then solved using a Kalman filter approach, a sequential strategy already studied for wave or elasticity problems [32, 27, 19]. We then revisit the UKF approach, to obtain a tangent-free algorithm with increasing convergence at each iteration of the outer loop of LM descent. Thus, we avoid the differentiation of the acoustoelastic wave propagation model and provide an estimation algorithm that can be easily interoperated with black box industrial codes. To be compatible with the dimensionality constraints of Kalman filtering, the parameter space is built from a modal decomposition of the pre-deformation, which in practice leads to the estimation of tens to hundreds of components on a modal basis. Moreover, the number of LM iterations to achieve convergence is rather small, leading overall to an almost sequential estimation approach. To illustrate the power of our approach, both in computation efficiency and parameter estimation, we provide extensive 3D results with synthetic noisy data in configurations associated with realistic SHM applications.

The structure of this article is as follows. In Section 2, we provide details on the direct problem, i.e. the time-domain acoustoelasticity propagation model. This leads us to the definition of the inverse problem we consider. In particular, we give a precise definition of the observation operator that generates the GW measurements and the modal basis that forms the parameter space. In Section 3, we develop our new method for identifying parameters using a combination of the LM algorithm and Kalman-based filtering. In a first step, we provide meaningful insights into this approach by considering continuous-time dynamics. In a second step, we give extensive details leading to the final fully discrete observer. Finally, in Section 4, we give relevant numerical illustrations of the estimation of the 3D pre-deformation with our approach and noisy synthetic data in realistic industrial configurations.

## 2. Problem setting

The starting point of our study is a problem in non-linear elastodnymics written in Lagrangian coordinates. Defining by $\Omega$ the reference domain with Lipschitz boundaries and the material position $\mathbf{x}$ in the reference configuration, we have

$$ \left\{ \begin{array}{ll} \varrho_0 \partial_{tt}^2 \mathbf{u}_{\text{tot}}(\mathbf{x}, t) - \nabla \cdot \mathbf{T}(\mathbf{x}, t) = \varrho_0 \mathbf{f}_{\text{tot}}(\mathbf{x}, t) & (\mathbf{x}, t) \in \Omega \times (0, T), \\ \mathbf{u}_{\text{tot}}(\mathbf{x}, t) = 0 & (\mathbf{x}, t) \in \Gamma_D \times (0, T), \\ \mathbf{T} \cdot \mathbf{n}(\mathbf{x}, t) = 0 & (\mathbf{x}, t) \in \Gamma_N \times (0, T), \\ \mathbf{u}(\mathbf{x}, 0) = 0 & \mathbf{x} \in \Omega, \end{array} \right. \quad (1) $$

where, for the sake of simplicity, we consider homogeneous boundary conditions defined in the reference configuration. In (1) the first Piola-Kirchhoff stress tensor $\mathbf{T}$ is here defined from the derivative of an hyperelastic potential $\mathcal{W}$, with respect to the deformation gradient $\mathbf{F}(\mathbf{x}, t) = \mathbf{Id} + \nabla \mathbf{u}_{\text{tot}}(\mathbf{x}, t)$, namely $\mathbf{T} = \mathrm{D}_{\mathbf{F}}\mathcal{W}(\mathbf{x}, \mathbf{F})$.

Let $\mathcal{K} \subset \mathcal{V}$ where $\mathcal{V}$ is defined by $\mathcal{V} = \{\mathbf{w} \in \mathrm{H}^1(\Omega)^3 \mid \mathbf{w}|_{\Gamma_D} = 0\}$ be the space of admissible displacements. The weak form of the equations above, also known as the
principle of virtual work, reads

$$ \forall \mathbf{w} \in \mathcal{V}, \quad \int_{\Omega} \varrho_0 \partial_{tt}^2 \mathbf{u}_{\text{tot}} \cdot \mathbf{w} \, \mathrm{d}\mathbf{x} + \int_{\Omega} \mathbf{T} : \nabla \mathbf{w} \, \mathrm{d}\mathbf{x} = \int_{\Omega} \varrho_0 \mathbf{f}_{\text{tot}} \cdot \mathbf{w} \, \mathrm{d}\mathbf{x}. \quad (2) $$

Regarding NDT applications, we assume that our system is inspected while being loaded by external unknown forces. These structural loading forces $f_0$ are considered to be volume distributed, for the sake of simplicity, and quasi-static – namely not depending on time – albeit of possibly strong amplitude, as opposed to the ultrasonic excitation $f$ which is of high-frequency and low amplitude. In fact, the quasi-static assumption is considered with respect to the time scale of ultrasonic excitation so the total loading decomposes into

$$ \mathbf{f}_{\text{tot}}(\mathbf{x}) = \mathbf{f}_0(\mathbf{x}) + \delta\mathbf{f}(\mathbf{x}, t), \quad (\mathbf{x}, t) \in \Omega \times (0, T), \quad (3) $$

where $\delta$ is a small parameter representing the fact that the amplitude of the ultrasonic excitation is small compared to the external loading. This allows us to separate scales and to consider that the resulting displacement is decomposed into

$$ \mathbf{u}_{\text{tot}}(\mathbf{x}, t) = \mathbf{u}_0(\mathbf{x}) + \delta\mathbf{u}(\mathbf{x}, t) + O(\delta^2). \quad (4) $$

Injecting such ansatz into the principle of virtual work (2) and identifying the zero order terms w.r.t. $\delta$, we formally showed in [12] that $\mathbf{u}_0 \in \mathcal{K}$ is a solution of a large displacement static problem

$$ \forall \mathbf{w} \in \mathcal{V}, \quad \int_{\Omega} \mathbf{T}(\mathbf{x}, \mathbf{F}_0) : \nabla \mathbf{w} \, \mathrm{d}\mathbf{x} = \int_{\Omega} \varrho_0 \mathbf{f}_0 \cdot \mathbf{w} \, \mathrm{d}\mathbf{x}, \quad (5) $$

where $\mathbf{F}_0 = \mathbf{Id} + \nabla \mathbf{u}_0$ is the deformation gradient associated with $\mathbf{u}_0$. Then, identifying the first order term in $\delta$, we find that $\mathbf{u}$ should be a solution to the wave-propagation problem given by, $\forall \mathbf{w} \in \mathcal{V}$,

$$ \int_{\Omega} \varrho_0 \partial_{tt}^2 \mathbf{u} \cdot \mathbf{w} \, \mathrm{d}\mathbf{x} + \int_{\Omega} \nabla \mathbf{u} : D_F^2 \mathcal{W}(\mathbf{x}, \mathbf{F}_0) : \nabla \mathbf{w} \, \mathrm{d}\mathbf{x} = \int_{\Omega} \varrho_0 \mathbf{f} \cdot \mathbf{w} \, \mathrm{d}\mathbf{x}. \quad (6) $$

From a mathematical viewpoint, the existence of a solution for the general problem formulation (1) is still an open problem. However, under suitable conditions, we can consider that problem (5) and (6) admit one and only one solution. For (5), this is typically the case when defining a hyperelastic law using a polyconvex potential [9], giving the existence of a displacement $\mathbf{u}_0$ in the space of admissible displacements.

Then, moving to (6), this formulation can be recast into a general, second-order in time, weak formulation. For conciseness, we define $\mathcal{H} = L^2(\Omega)^3$ as well as the scalar products

$$ \forall (\mathbf{u}, \mathbf{w}) \in \mathcal{H} \times \mathcal{H}, \quad (\mathbf{u}, \mathbf{w})_{\mathcal{H}} = \int_{\Omega} \varrho_0 \mathbf{u} \cdot \mathbf{w} \, \mathrm{d}\mathbf{x}, $$

$$ \forall (\mathbf{u}, \mathbf{w}) \in \mathcal{V} \times \mathcal{V}, \quad (\mathbf{u}, \mathbf{w})_{\mathcal{V}} = \int_{\Omega} \nabla \mathbf{u} : D_F^2 \mathcal{W}(\mathbf{x}, \mathbf{Id}) : \nabla \mathbf{w} \, \mathrm{d}\mathbf{x} + (\mathbf{u}, \mathbf{w})_{\mathcal{H}}. $$

Note that $D_F^2 \mathcal{W}(\mathbf{x}, \mathbf{Id})$ actually corresponds to the standard Hooke's law in linear elasticity. Then, we introduce the linear operator $A_0(\mathbf{u}_0) \in \mathcal{L}(\mathcal{D}(A_0), \mathcal{H})$ such that

$$ \forall (\mathbf{u}, \mathbf{w}) \in \mathcal{D}(A_0) \times \mathcal{V}, \quad (A_0(\mathbf{u}_0)\mathbf{u}, \mathbf{w})_{\mathcal{H}} = \int_{\Omega} \nabla \mathbf{u} : D_F^2 \mathcal{W}(\mathbf{x}, \mathbf{F}_0) : \nabla \mathbf{w} \, \mathrm{d}\mathbf{x} $$
The definition of $\mathcal{D}(A_0)$ is omitted for simplicity, but should be seen as a space of functions with square integrable second-order derivatives. The weak formulation (6) can be written in the following form

$$ \forall \mathbf{w} \in \mathcal{V}, \quad \frac{d^2}{dt^2}(\mathbf{u}, \mathbf{w})_{\mathcal{H}} + (A_0(\mathbf{u}_0)\mathbf{u}, \mathbf{w})_{\mathcal{H}} = (\mathbf{f}, \mathbf{w})_{\mathcal{H}}. \qquad (7) $$

We now make the following assumption: there exists $\lambda \ge 0$ and $\alpha > 0$ such that

$$ \forall \mathbf{w} \in \mathcal{D}(A_0), \quad (A_0(\mathbf{u}_0)\mathbf{w}, \mathbf{w})_{\mathcal{H}} + \lambda \|\mathbf{w}\|_{\mathcal{H}}^2 \ge \alpha \|\mathbf{w}\|_{\mathcal{V}}^2. \qquad (8) $$

Such assumption is satisfied for $\mathbf{u}_0$ small enough whereas for large displacement $\mathbf{u}_0$, it may be violated. This is however very dependent on the potential $\mathcal{W}$ that describes the elastic behavior of the medium. Some choices tend to make the assumption above more restrictive than others – see e.g. the examples in [12]. Thanks to this assumption – see for instance [16] – there exists one, and only one, variational solution of (7). Then, for a given $\mathbf{u}_0 \in \mathcal{K}$, by introducing the operator

$$ A(\mathbf{u}_0) = \begin{pmatrix} 0 & \mathrm{Id} \\ -A_0(\mathbf{u}_0) & 0 \end{pmatrix} \quad \text{with} \quad \mathcal{D}(A) = \mathcal{D}(A_0) \times \mathcal{V} \subset \mathcal{Z} := \mathcal{V} \times \mathcal{H}, $$

we can rewrite the wave dynamics in the state-space form

$$ \left\{ \begin{array}{ll} \frac{d}{dt}z = A(\mathbf{u}_0)z + r, & \text{in } (0, T) \\ z(0) = 0, & \text{with } z = \begin{pmatrix} \mathbf{u} \\ \mathbf{v} \end{pmatrix}, \text{ and } r = \begin{pmatrix} 0 \\ \mathbf{f} \end{pmatrix}, & (9) \end{array} \right. $$

and $\mathbf{v}$ denotes the velocity unknown. In the sequel, we will denote $\dot{z} = \frac{d}{dt}z$. This state-space form allows to define solutions using semi-group theory [6]. Namely for $r \in L^2((0, T), \mathcal{Z})$ there exists one and only one mild solution of (9) in $C^0((0, T), \mathcal{Z})$ which is also a variational solution of (7). We here underline that the variational formalism is necessary to further justify finite element discretization of (7), whereas the mild solution formulation in the sense of semi-group will simplify the presentation of the estimation problem.

In this modeling context, we consider that we have at our disposal some recorded measurements of a target wave field $\tilde{\mathbf{u}}$ from $d$ sensors occupying the subregions $\{\omega_i\}_{i=1}^d$. Typically we consider recordings $y = (y_i)_{1 \le i \le d}$ given by

$$ [0, T] \ni t \mapsto \tilde{y}_i(t) = \frac{1}{|\omega_i|} \int_{\omega_i} \tilde{\mathbf{u}}(t, \mathbf{x}) \cdot \mathbf{d}_i \, \mathrm{d}\mathbf{x} \in \mathbb{R}, \quad 1 \le i \le d, $$

where the fields $\mathbf{d}_i$ represent the sensitivity of the sensors to a displacement field. In fact, we do not have at our disposal $[0, T] \ni t \mapsto \tilde{y}(t) \in \mathbb{R}$ but rather the perturbed measurement $y_\gamma \in L^2((0, T); \mathcal{Y})$ with $\mathcal{Y} = \mathbb{R}^d$ such that for a given noise level $\gamma > 0$,

$$ \int_0^T \|\tilde{y}_i - y_{\gamma,i}\|_{\mathcal{Y}}^2 \, \mathrm{d}t \lesssim \gamma^2 T. $$

From a state-space point of view, this allows to define an observation operator $C \in \mathcal{L}(\mathcal{Z}, \mathcal{Y})$ by

$$ C: z = \begin{pmatrix} \mathbf{u} \\ \mathbf{v} \end{pmatrix} \mapsto \left( \frac{1}{|\omega_i|} \int_{\omega_i} \mathbf{u}(\mathbf{x}) \cdot \mathbf{d}_i \, \mathrm{d}\mathbf{x} \right)_{1 \le i \le d}. \qquad (10) $$
**Remark 1** We would like to point out that in an alternative modeling approach, we could have considered that the measurements are recorded from the deformed configuration as the sensors are operating in this configuration. When rewriting the integrals with respect to the reference configuration, this should lead to the definition of an observation operator $C(\mathbf{u}_0)$ that depends on the deformation.

We can now introduce the inverse problem addressed in our work. We assume that $\mathbf{u}_0$ is unknown and we want to reconstruct it from the available measurements. More precisely, we would like to specify the true $\check{\mathbf{u}}_0$ from an a priori displacement $\hat{\mathbf{u}}_0$ assuming that the true displacement is a regular perturbation with respect to the a priori, typically there exists a constant $M$ such that

$$ \|\check{\mathbf{u}}_0 - \hat{\mathbf{u}}_0\|_{\mathcal{V}}^2 \leq M^2. $$

Note that $M^{-1}$ actually represents the level of confidence in the given parameter prior $\hat{\mathbf{u}}_0$. In practice, we propose to decompose any $\mathbf{u}_0 \in \mathcal{H}$ on the basis made of the eigenvectors $(\varphi_j)_{j \geq 0}$ of the compact operator $\Lambda_0 = A_0(\mathbf{0})^{-1} \in \mathcal{L}(\mathcal{V}', \mathcal{V})$. We recall that there exists $(\mu_j)_{j \geq 0} \in (\mathbb{R}_+^*)^N$ such that

$$ A_0(\mathbf{0})\varphi_j = \mu_j \varphi_j \text{ with } \|\varphi_j\|_{\mathcal{H}}^2 = 1 \quad (11) $$

and $\lambda_j = \sqrt{\mu_j^{-1}}$ tends to 0 as $j$ tends to infinity. Therefore, we can decompose

$$ \mathbf{u}_0(\mathbf{x}) = \sum_{j \geq 0} \theta_j \varphi_j(\mathbf{x}) \quad (12) $$

and we can enforce

$$ \mathbf{u}_0 \in \mathcal{V} \Leftrightarrow \sum_{j \geq 0} \frac{\theta_j^2}{\lambda_j^2} < +\infty. $$

In practice, we may even want to consider that $\mathbf{u}_0$ belongs to a more regular space $\mathcal{V}^m$ leading to increased convergence rate of the sequence $(\theta_j)_{k \geq 0}$ to 0, typically

$$ \mathbf{u}_0 \in \mathcal{V}^m \Rightarrow \sum_{j \geq 0} \frac{\theta_j^2}{\lambda_j^{2m}} < +\infty, \quad (13) $$

which is equivalent to estimating $\theta \in \mathcal{P} \subset \ell_2(\mathbb{R})$, equipped with a norm

$$ \|\theta\|_{\mathcal{P}}^2 = (\theta, \Lambda_0^{-m}\theta)_{\ell_2}, \text{ with } (\Lambda_0)_{ij} = \lambda_j^{2m}\delta_{ij} \text{ and } \delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{otherwise.} \end{cases} \quad (14) $$

When discretizing $\mathbf{u}_0$, we will consider a finite-dimensional space corresponding to a finite number $N_p$ of modes, hence $\theta \in \mathbb{R}^{N_p}$ and $\Lambda_0$ becomes typically a diagonal matrix.

Therefore, and to conclude our problem setting, our objective is to estimate $\check{\theta} \in \mathcal{P}$ from the available measurements, namely to invert the following operator

$$ \Psi_T: \begin{cases} \mathcal{P} \rightarrow L^2((0, T); \mathcal{Y}) \\ \theta \mapsto \left[ t \mapsto \int_0^t Ce^{A(\theta)(t-s)}r(s) \, ds \right] \end{cases} $$

where we replace the $\mathbf{u}_0$ dependency by the $\theta$ dependency in the definition of the operator $A_0$. In the following, this inversion is based on a least-squares minimization using Levenberg-Marquardt strategy.
### 3. Identification method

#### 3.1. From Levenberg-Marquardt minimization scheme to an iterated Extended Kalman strategy

Following [4], we approximate pseudo-inverse of $\Psi_T$ with a Levenberg-Marquardt minimization scheme [23] which consists of the following iterative procedure.

$$ \theta^{k+1} = \theta^k + \alpha^k \\ = \theta^k + \left[ D\Psi_T(\theta^k)^* D\Psi_T(\theta^k) + \varepsilon \frac{\gamma^2}{M^2} \Lambda_0^{-m} \right]^{-1} D\Psi_T(\theta^k)^* (y_\gamma - \Psi_T(\theta^k)), $$

where the number of iterations of this scheme is finite, thus acting as a regularization, and typically controlled by classical Morozov-like criteria, as advised in [4]. The parameter $\frac{\gamma^2}{M^2}$ is a scaling parameter taking into account the prior over noise ratio, while $\varepsilon$ will give us an additional degree of freedom for weighting the regularization in the LM algorithm, see below and in particular Remark 2. Then, we can see the increment as

$$ \alpha^k = \underset{\alpha \in \mathcal{P}}{\operatorname{argmin}} \left\{ \frac{\varepsilon}{2M^2} \|\alpha\|_{\mathcal{P}}^2 + \frac{1}{2\gamma^2} \|y_\gamma - \Psi_T(\theta^k) - D\Psi_T(\theta^k)\alpha\|_{L^2((0,T);\mathcal{Y})}^2 \right\}, \quad (15) $$

**Remark 2** Note that by penalizing $\|\alpha\|_{\mathcal{P}}^2$ we penalize an incremental displacement $\|\tilde{u}_0(\alpha)\|_{V^m}^2$. Therefore, we can force the first increments of the seeking displacement to belong to $\mathcal{K}$ since a sufficiently smooth displacement with a sufficiently small amplitude will necessarily be admissible. However, this may come at the price of too much regularization, which can only be compensated by a higher number of iterations of the LM algorithm. This is illustrated in the numerical section.

Let us now specify the tangent operator

$$ \forall \theta \in \mathcal{P}, \quad D\Psi_T(\theta) : \left. \begin{array}{l} \mathcal{H} \to L^2((0,T); \mathcal{Y}) \\ \alpha \mapsto y = C\zeta_{|\theta,\alpha} \end{array} \right\} $$

and $\zeta_{|\theta,\alpha}$ is a mild solution of

$$ \left\{ \begin{array}{l} \dot{\zeta}_{|\theta,\alpha}(t) = A(\theta)\zeta_{|\theta,\alpha}(t) + B(\theta, z_{|\theta}(t))\alpha, \quad t \in [0, T] \\ \zeta_{|\theta,\alpha}(0) = 0, \end{array} \right. $$

where we have introduced the linear operator representing the tangent of $A(\theta)$ w.r.t. the parameter $\theta$

$$ \forall \theta \in \mathcal{P}, \quad \forall z \in \mathcal{D}(A), \quad \mathcal{L}(\mathcal{P}, \mathcal{Z}) \ni B(\theta, z) : \mathcal{P} \ni \alpha \mapsto (D_{\theta}A(\theta)\alpha)z \in \mathcal{Z}. $$

Therefore, we face a linear-quadratic optimal control problem, that can be minimized using the following Reduced-Order Kalman Filter (ROKF) sequential estimator reformulated from the initial design found in [32] – see also the more recent review [30]. To this end, let us first introduce the set of sensitivity operators $(L_{|\theta}(t))_{t \ge 0}$ in $\mathcal{L}(\mathcal{P}, \mathcal{Z})$ defined for all time $t \ge 0$ by

$$ L_{|\theta}(t) : \alpha \mapsto \zeta(t) \text{ the mild solution of} \quad (16) $$

$$ \left\{ \begin{array}{l} \dot{\zeta}(s) = A(\theta)\zeta(s) + B(\theta, z_{|\theta}(s))\alpha, \quad s \in [0, t] \\ \zeta(0) = 0. \end{array} \right. \quad (17) $$
We easily verify that $L_{|\theta} \in C^0([0, T]; \mathcal{L}(\mathcal{P}, \mathcal{Z}))$ the space of continuous mapping from $[0, T]$ to $\mathcal{L}(\mathcal{P}, \mathcal{Z})$ endowed with the uniform convergence topology [6].

Then, we introduce the time-dependent Riccati operator $(\Lambda_{|\theta})_{t \ge 0}$ which can be interpreted as a parameter-covariance operator [6] – solution in $C^0([0, T]; S_+^*(\mathcal{P}))$, with $S_+^*(\mathcal{P})$ the space of symmetric positive definite bounded operators – of

$$ \left\{ \begin{array}{l} \dot{\Lambda}_{|\theta}(t) = -\frac{1}{\gamma^2} \Lambda_{|\theta}(t) L_{|\theta}(t)^* C^* C L_{|\theta}(t) \Lambda_{|\theta}(t), \quad t \in [0, T] \\ \Lambda_{|\theta}(0) = \frac{M^2}{\epsilon} \Lambda_0^m. \end{array} \right. \qquad (18) $$

We would like to underline the fact that here the adjoint operator $L_{|\theta}(t)^*$ is defined with respect to the $\ell^2$-norm as $\theta \in \mathcal{P} \subset \ell^2(\mathbb{R})$.

We then define the sequential estimator

$$ \left\{ \begin{array}{l} \dot{\hat{\zeta}}^k(t) = A(\theta^k)\hat{\zeta}^k(t) + B(\theta^k, z_{|\theta^k}(t))\hat{\alpha}^k(t) + L_{|\theta^k}(t)\dot{\hat{\alpha}}^k(t), \quad t \in [0, T] \\ \dot{\hat{\alpha}}^k(t) = \frac{1}{\gamma^2} \Lambda_{|\theta^k}(t) L_{|\theta^k}(t)^* C^* (y_\gamma(t) - Cz_{|\theta^k}(t) - C\hat{\zeta}^k(t)), \quad t \in [0, T] \\ \hat{\zeta}^k(0) = 0 \\ \hat{\alpha}^k(0) = 0 \end{array} \right. \qquad (19) $$

which, ultimately, sequentially solves the minimization problem (15) as recall in the next theorem justified in [30].

**Theorem 1** At every iteration $k$, the mild solution $\hat{\alpha}^k \in C^0([0, T]; \mathcal{P})$ and $\hat{\zeta}^k \in C^0([0, T]; \mathcal{Z})$ of (19) satisfies

$$ \hat{\alpha}^k(T) = \alpha^k \text{ and } \hat{\zeta}^k(T) = \zeta_{|\theta^k, \alpha^k}(T). \qquad (20) $$

In the previous theorem, we understand that we can sequentially compute $t \mapsto \hat{\alpha}^k$ but also $t \mapsto \hat{\zeta}^k$ which is interpreted as a sequential estimator of the associated trajectory $t \mapsto \zeta_{|\theta^k, \alpha^k}$. Note that this estimator can be computed together with $\hat{\alpha}^k$ in only one coupled forward dynamics.

In fact, instead of a parameter increment, we are more interested in reconstructing the parameter itself, namely $t \mapsto \bar{\theta}^k(t) = \theta^k + \hat{\alpha}^k(t)$, and its associated trajectory $t \mapsto \bar{z}^k(t) = z_{|\bar{\theta}^k}(t)$. As $t \mapsto \hat{\alpha}^k(t)$ is a time-dependent trajectory, application of a chain rule yields

$$ \dot{\bar{z}}^k = \frac{d}{dt} [z_{|\bar{\theta}^k}(\cdot)] = \dot{z}_{|\bar{\theta}^k}(\cdot) + D_{\theta^k}(z_{|\theta^k}(\cdot)) \Big|_{\theta^k=\bar{\theta}^k(\cdot)} \dot{\hat{\alpha}}^k, $$

and the sensitivity operator $L_{|\bar{\theta}^k}(\cdot) = D_{\theta^k}(z_{|\theta^k}(\cdot))|_{\theta^k=\bar{\theta}^k(\cdot)}$ is nothing else than, for all time $t \ge 0$,

$$ L_{|\bar{\theta}^k}(\cdot)(t) : \alpha \mapsto \zeta(t) \text{ is a mild solution of} \\ \left\{ \begin{array}{l} \dot{\zeta}(s) = A(\bar{\theta}^k(s))\zeta(s) + B(\bar{\theta}^k(s), z_{|\bar{\theta}^k}(s))\alpha, \quad s \in [0, t] \\ \zeta(0) = 0 \end{array} \right. \qquad (21) $$

Here note that in the previous dynamics, $s \mapsto \bar{\theta}^k(s)$ is a function of time, implying a time-dependent operator $s \mapsto A(\bar{\theta}^k(s))$. The existence of a solution for such dynamics is fundamentally based on the study of evolution equation operator in general [35, Chapter 5].
As $\dot{\hat{\alpha}}^k = \dot{\bar{\theta}}^k$, and by approximating $z_{|\theta^k} + \hat{\zeta}^k \simeq z_{|\theta^k + \hat{\alpha}^k} = \bar{z}^k$ we could compute the sequential estimator

$$ 
\left\{ 
\begin{array}{ll}
\dot{\bar{z}}^k = A(\bar{\theta}^k)\bar{z}^k + r + L_{|\bar{\theta}^k} \dot{\bar{\theta}}^k, & t \in [0, T], \\
\dot{\bar{\theta}}^k = \frac{1}{\gamma^2} \Lambda_{|\theta^k}(t) L_{|\theta^k}(t)^* C^* (y_\gamma - C\bar{z}^k), & t \in [0, T], \\
\bar{z}^k(0) = 0 \\
\bar{\theta}^k(0) = \bar{\theta}^k,
\end{array}
\right. 
\quad (22) $$

such that at final time $\bar{\theta}^k(T) \simeq \bar{\theta}^{k+1}$.

One drawback in solving (22) is that it necessitates to propagate two sensitivity operators $L_{|\theta^k}$ – which implies storing the complete trajectory $z_{|\theta^k}$ – and $L_{|\bar{\theta}^k}$. To circumvent this drawback we propose to solve instead

$$ 
\left\{ 
\begin{array}{ll}
\dot{\hat{z}}^k = A(\hat{\theta}^k)\hat{z}^k + r + L^k \dot{\hat{\theta}}^k, & t \in [0, T], \\
\dot{\hat{\theta}}^k = \frac{1}{\gamma^2} \Lambda^k(t) L^k(t)^* C^* (y_\gamma - C\hat{z}^k), & t \in [0, T], \\
\hat{z}^k(0) = 0 \\
\hat{\theta}^k(0) = \hat{\theta}^{k-1}(T)
\end{array}
\right. 
\quad (23) $$

where, from now on, we simply use the notation $L^k = L_{|\hat{\theta}^k}$ and $\Lambda^k = \Lambda_{|\hat{\theta}^k}$ and we see that the sensitivity and the covariance operator are updated through time. Indeed, $\Lambda^k$ is an operator that should be seen as a mild solution of the dynamics

$$ 
\left\{ 
\begin{array}{l}
\dot{\Lambda}^k(t) = -\frac{1}{\gamma^2} \Lambda^k(t) L^k(t)^* C^* L^k(t) \Lambda^k(t), \quad t \in [0, T] \\
\Lambda^k(0) = \frac{M^2}{\varepsilon} \Lambda_0^m,
\end{array}
\right. 
\quad (24) $$

while $L^k$ can be seen from (21) as a mild solution in $C^0([0, T], \mathcal{L}(\mathcal{P}, \mathcal{Z}))$ of

$$ 
\left\{ 
\begin{array}{l}
\dot{L}^k(t) = A(\hat{\theta}^k)L^k(t) + B(\hat{\theta}^k(t), \hat{z}^k(t)) \quad t \in [0, T] \\
L^k(0) = 0.
\end{array}
\right. 
\quad (25) $$

We here recognize in (23) – combined with (24)-(25) – an iterated version of the Reduced-Order-Extended-Kalman-Filter (ROEKF) estimator proposed in [32] for joint state and parameter estimation for wave-like equations, here formally generalized to infinite dimensional systems. One very strong advantage of such an estimator is that it is fully sequential in the sense that it does not require storing any trajectory in the iteration procedure. Note that the existence of a solution of (23) is much more intricate and is based on justifying that $L^k$ is well-defined even for time-dependent parameter $\hat{\theta}^k$. This question is typically covered by the definition of stable families of generators in evolution equation [35], see also [2] for similar questions when defining the Extended Kalman Filter for infinite dimensional systems.

In this work, we call iROEKF the proposed iterated ROEKF. However, we must warn the reader that in our case the iteration $k$ is to be considered as an outer loop of the LM descent. This iROEKF is therefore different from the classical iterated Extended Kalman Filter presented in the literature [42].

**Remark 3** The separation of amplitudes (4) comes from the fact that most of the ultrasonic experiments performed for NDT (and in particular guided-wave based SHM)
are based on “low-amplitude” elastic wave propagation. Hence, it seems reasonable to assume infinitesimal deformations for the ultrasonic beam. If this assumption is, in fact, not verified (for instance in the context of high intensity focused ultrasound) then a fully nonlinear model should be considered. This model would couple the quasi-static deformation with the ultrasonic wave propagation, at the cost of a significant increase in the computational load for solving the forward problem. Note that the presented filtering approach could take this into account using a nonlinear state equation in the LM iterations, namely the first equation in (23) would become

$$ \dot{\hat{z}}^k(t) = A(\hat{\theta}^k(t), \hat{z}^k(t)) + r(t) + L^k(t)\dot{\hat{\theta}}^k(t), \quad t \in [0, T], $$

and using the tangent of the operator w.r.t. the state variable in the sensitivity dynamics, namely (25) would read

$$ \dot{L}^k(t) = D_z A(\hat{\theta}^k(t), \hat{z}^k(t)) L^k(t) + B(\hat{\theta}^k(t), \hat{z}^k(t)), \quad t \in [0, T]. $$

## 3.2. Space-time-discretized version of the iterated reduced-order Extended Kalman Filter

Since we ultimately solve a discretized version of (23), we now present a discretization strategy based on a stable discretization of the Levenberg-Marquardt increment estimator (19). Our discretization strategy is based on the fundamental principle *discretize-then-optimize*, which means that we first discretize the direct problem and then reapply the equivalence of Levenberg-Marquardt and dynamic programming at the discretized level, as recommended in [30]. For (7), we consider a finite-element discretization based on the Spectral Finite Element Method (SFEM) [10] leading to the formulation

$$ \forall \mathbf{w}_h \in \mathcal{V}_h, \quad (\partial_{tt}^2 \mathbf{u}_h, \mathbf{w}_h)_{\mathcal{V}_h} + (A_{0h}(\mathbf{u}_{h0})\mathbf{u}_h, \mathbf{w}_h)_{\mathcal{V}_h} = (\mathbf{f}, \mathbf{w}_h)_{\mathcal{V}_h}. \quad (26) $$

To construct the finite dimensional space $\mathcal{V}_h$ we assume given a partition $\mathcal{T}_h$ of quadrangles (in 2D) or hexahedra (in 3D) of the domain $\Omega$, namely

$$ \bar{\Omega} = \bigcup_{K \in \mathcal{T}_h} K, \quad \forall (K, L) \in \mathcal{T}_h \times \mathcal{T}_h, \quad \mathring{K} \cap \mathring{L} = \emptyset, $$

without hanging nodes and with maximum diameter given by $h$. Then $\mathcal{V}_h \subset \mathcal{V}$ is obtained using a $\mathbb{Q}_k$-Lagrangian basis on a set of nodes $\{\xi_i\}_{i=1}^{N_h}$, such that $\mathcal{V}_h = (\text{span}\{\varphi_i\}_{i=1}^{N_h})^3$, with

$$ \forall 1 \le i, j \le N_h, \quad \varphi_i(\xi_j) = \delta_{ij}, \quad \forall K \in \mathcal{T}_h, \quad \varphi_{i|K} \circ F_K \in \mathbb{Q}_k. $$

where the mapping from the reference element to any element $K \in \mathcal{T}_h$ is denoted $F_K$ – see [10] for more details. The nodes $\{\xi_i\}$ are obtained using Gauss-Lobatto integration points on a reference square (in 2D) or cube (in 3D). One fundamental aspect of the efficiency of the SFEM is that $\mathcal{V}_h$ is equipped with a scalar product that – using a quadrature formula on the nodes $\{\xi_i\}$ – leads at the algebraic level to a diagonal mass matrix. This is called a mass-lumping strategy and preserves a near-optimal accuracy
[15]. The mentioned quadrature formula is also employed to compute the operator $A_h(\mathbf{u}_{h0})$, thus given by, for any pairs of test functions $(\mathbf{u}_h, \mathbf{w}_h) \in \mathcal{V}_h^2$,

$$ (A_{0h}(\mathbf{u}_{h0})\mathbf{u}_h, \mathbf{w}_h)_{\mathcal{V}_h} = \int_{\Omega} \nabla \mathbf{u}_h : D_{\mathbf{F}}^2 \mathcal{W}(\mathbf{x}, \mathbf{F}_{0h}) : \nabla \mathbf{w}_h \, d\mathbf{x}, $$

where $\int_{\Omega}$ stands for an integral computed from the quadrature formula, and $\mathbf{F}_{0h} = \mathbf{Id} + \nabla \mathbf{u}_{h0}$ with $\mathbf{u}_{h0}$ an interpolation in $\mathcal{V}_h$ of $\mathbf{u}_0$. Note that the operator inherits some properties of the operator $A(\mathbf{u}_{h0})$, in particular it is self-adjoint. However it is not clear that it satisfies a positivity property of the form (8) for at least two reasons, the interpolation of $\mathbf{u}_{h0}$ and the use of the quadrature formula. Therefore we are led to the following assumption that is sufficient to have a well-posed discrete problem,

$$ \forall \mathbf{w}_h \in \mathcal{V}_h, \quad (A_{0h}(\mathbf{u}_{h0})\mathbf{w}_h, \mathbf{w}_h)_{\mathcal{V}_h} \geq 0. $$

The spatial discretization is then followed by an explicit time-discretization. Being given a time step $\Delta t$, the solution $\mathbf{u}_h(t)$ is approximated at time $t^n = n\Delta t$ by solving for any test function $\mathbf{w}_h \in \mathcal{V}_h$

$$ \left( \frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{w}_h \right)_{\mathcal{V}_h} + (A_{0h}(\mathbf{u}_{h0})\mathbf{u}_h^n, \mathbf{w}_h)_{\mathcal{V}_h} = (\Pi_h \mathbf{f}(t^n), \mathbf{w}_h)_{\mathcal{V}_h}, \quad (27) $$

where $\Pi_h \mathbf{f}$ is a projection of $\mathbf{f}$ in $\mathcal{V}_h$. Thanks to the mass-lumping strategy the computation of $\mathbf{u}_h^{n+1}$ involves – at the algebraic level – only the inversion of a diagonal mass matrix. This discretization is stable for small enough time step $\Delta t$, namely the time step must satisfy the CFL condition

$$ \Delta t \le 2 \left( \sup_{\mathbf{w}_h \in \mathcal{V}_h} \frac{(A_{0h}(\mathbf{u}_{h0})\mathbf{w}_h, \mathbf{w}_h)_{\mathcal{V}_h}}{\|\mathbf{w}_h\|_{\mathcal{V}_h}^2} \right)^{-1/2}. \quad (28) $$

The CFL condition depends on the mesh size $h$ (in particular $\Delta t$ behaves like $O(h)$) and on the mesh quality, but also on the gradient of the displacement field $\mathbf{u}_{h0}$. We can finally rewrite the time-discretization (27) as

$$ \left\{ \begin{array}{l} \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t} = \mathbf{v}_h^{n+1}, \\ \\ \frac{\mathbf{v}_h^{n+1} - \mathbf{v}_h^n}{\Delta t} + A_{0h}(\theta_h)\mathbf{u}_h^n = \Pi_h \mathbf{f}(t^n), \end{array} \right. \quad (29) $$

where we replace the $\mathbf{u}_{h0}$ dependency by the $\theta_h \in \mathcal{P}_h \simeq \mathbb{R}^{N_p}$ dependency such that

$$ \mathbf{u}_{h0} = \sum_{j=1}^{N_p} \theta_{h,j} \varphi_{h,j} \quad \text{with} \quad A_{0h}(0) \varphi_{h,j} = \mu_{h,j} \varphi_{h,j}, $$

namely the $(\mu_{h,j}, \varphi_{h,j})$ are the eigenelements of the discrete operators $A_{0h}(0)$ (ordered increasingly with respect to the eigenvalues). Then, by defining $z_h^n = (\mathbf{u}_h^n, \mathbf{v}_h^n)^\top \in \mathcal{Z}_h = \mathcal{V}_h \times \mathcal{V}_h$, we rewrite (29) in the following abstract state-space form

$$ \left\{ \begin{array}{l} z_h^{n+1} = \Phi_{h,\Delta t}(\theta_h) z_h^n + r_{h,\Delta t}^{n+1}, \\ z_h^0 = 0, \end{array} \right. \quad (30) $$

where

$$ \Phi_{h,\Delta t}(\theta_h) : z_h \mapsto \left( \begin{array}{cc} \mathrm{Id}_{\mathcal{V}_h} - \Delta t^2 A_{0h}(\theta_h) & \Delta t \mathrm{Id}_{\mathcal{V}_h} \\ \Delta t A_{0h}(\theta_h) & \mathrm{Id}_{\mathcal{V}_h} \end{array} \right) z_h, $$
and

$$r_{h,\Delta t}^{n+1} = \begin{pmatrix} \Delta t^2 \Pi_h \boldsymbol{f}(t^n) \\ \Delta t \Pi_h \boldsymbol{f}(t^n) \end{pmatrix}.$$

Using the time-discretized system (30), we propose to solve the Levenberg-Marquardt procedure

$$\begin{aligned} \theta_h^{k+1} &= \theta_h^k + \alpha_h^k \\ &= \theta_h^k + \left[ D\Psi_{h,N}(\theta^k)^* D\Psi_{h,N}(\theta^k) + \varepsilon \frac{\gamma^2}{M^2} \Lambda_{0h}^{-m} \right]^{-1} D\Psi_{h,N}(\theta_h^k)^* (y_\gamma - \Psi_{h,N}(\theta_h^k)), \end{aligned}$$

where

$$\Psi_{h,N} : \left\{ \begin{array}{ll} \mathcal{P}_h & \rightarrow \ell^2((1:N); \mathcal{Y}) \\ \theta_h & \mapsto \left\{ \sum_{j=1}^{n} C_h \Phi_{h,\Delta t}^{n-j}(\theta_h) r_{h,\Delta t}^{j} \right\}_{n=1}^{N}, \end{array} \right.$$

which can be equivalently rewritten in the following form

$$\alpha_h^k = \underset{\alpha_h \in \mathcal{P}}{\operatorname{argmin}} \left\{ \frac{\varepsilon}{2M^2} \|\alpha_h\|_{\mathcal{P}}^2 + \frac{\Delta t}{2\gamma^2} \|y_\gamma - \Psi_{h,N}(\theta^k) - D\Psi_{h,N}(\theta_h^k)\alpha_h\|_{\ell^2((1:N);\mathcal{Y})}^2 \right\}.$$

Then, we define the fully discrete sensitivity operator, for all $n \in [0:N]$, as

$$L_{h|\theta_h}^{k,n} : \alpha_h \mapsto \zeta_h^n \quad \text{the solution of} \\ \left\{ \begin{array}{l} \zeta_h^{j+1} = \Phi_{h,\Delta t}(\theta_h)\zeta_h^j + B_{h,\Delta t}(\theta_h, z_h^j)\alpha_h, \quad j \in [0:n-1] \\ \zeta_h^0 = 0, \end{array} \right.$$

with

$$B_{h,\Delta t}(\theta_h, z_h)\alpha_h = \begin{pmatrix} -\Delta t^2 D_{\theta} A_{0h}(\theta_h)\alpha_h & 0 \\ \Delta t D_{\theta} A_{0h}(\theta_h)\alpha_h & 0 \end{pmatrix} z_h.$$

Moreover, we define the fully discrete covariance operator $\Lambda_{h|\theta_h}^{k,n}$ satisfying the dynamics

$$\left\{ \begin{array}{l} (\Lambda_{h|\theta_h}^{k,n+1})^{-1} = (\Lambda_{h|\theta_h}^{k,n})^{-1} + \frac{\Delta t}{\gamma^2} L_{h|\theta_h}^{k,n+1*} C_h^* C_h L_{h|\theta_h}^{k,n+1}, \quad n \in [0:N] \\ \Lambda_{h|\theta_h}^{k,0} = \frac{M^2}{\varepsilon} \Lambda_{0h}^m. \end{array} \right.$$

If we now choose to solve from $\hat{\zeta}_h^{k,0} = 0$, and $\hat{\alpha}_h^{k,0} = 0$, for all $k \le 0$ and $n \in [0:N-1]$

$$\left\{ \begin{aligned} \hat{\zeta}_h^{k,n+1-} &= \Phi_{h,\Delta t}(\theta_h^k)\hat{\zeta}_h^{k,n} + B_{h,\Delta t}(\theta_h^k, z_h^n)\hat{\alpha}_h^{k,n} \\ \hat{\alpha}_h^{k,n+1} &= \hat{\alpha}_h^{k,n} + \frac{\Delta t}{\gamma^2} \Lambda_{h|\theta_h}^{k,n+1} L_{h|\theta_h}^{k,n+1*} C_h^*(y_\gamma^{n+1} - C_h z_h^{n+1} - C_h \hat{\zeta}_h^{k,n+1-}) \\ \hat{\zeta}_h^{k,n+1} &= \hat{\zeta}_h^{k,n+1-} + L_{h|\theta_h}^{k,n+1}(\hat{\alpha}_h^{k,n+1} - \hat{\alpha}_h^{k,n}) \end{aligned} \right. \quad (31)$$

where the “n-” and “n” exponents correspond to the two steps of a splitting time-scheme for $\hat{\zeta}_h^{k,n}$ – where the “n-” step is often referred as the prediction (or forecast) step and the “n” step is referred as the correction (or analysis) step – we have at the discrete level the following equivalence proved in [30].

**Theorem 2** At every iteration $k$ of the Levenberg-Marquardt minimization, we have

$$\hat{\alpha}_h^{k,N} = \alpha_h^k \quad \text{and} \quad \hat{\zeta}_h^{k,N} = \zeta_{h|\theta^k, \alpha^k}^N \quad (32)$$
The last theorem helps to understand that we keep at the discrete-time level the same equivalence between the least square minimization of a discretized criterion and the sequential approach. Therefore pursuing the same strategy than for the continuous-time level leads us to defining the Discrete-Time iterated Reduced Order Extended Kalman Filter (DT-iROEKF) as an approximation of the Levenberg-Marquardt minimization strategy for discrete-time dynamics. The resulting joint state-parameter sequential estimator based on this DT-iROEKF is initialized from $\hat{z}_h^k = 0, \hat{\theta}_h^{k,0} = \hat{\theta}_h^{k-1,N}$ if $k > 1$ and $\hat{\theta}_h^{1,0} = \hat{\theta}_{0h}$, while $L^{k,0} = 0, \Lambda^{k,0} = \varepsilon^{-1}M^2\Lambda_{h0}^m$. Then the recursive dynamics reads for all $k \ge 1$ and for all $n \in [0:N-1]$

— Prediction / Forecast:

$$ 
\begin{cases}
\hat{z}_h^{k,n+1-} = \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n})\hat{z}_h^{k,n} + r_{h,\Delta t}^{n+1}, \\
L_h^{k,n+1} = \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n})L_h^{k,n} + B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n}).
\end{cases}
$$

— Correction / Analysis:

$$ 
\begin{cases}
(\Lambda_h^{k,n+1})^{-1} = (\Lambda_h^{k,n})^{-1} + \frac{\Delta t}{\gamma^2} L_h^{k,n+1*} C_h^* C_h L_h^{k,n+1}, \\
\hat{\theta}_h^{k,n+1} = \hat{\theta}_h^{k,n} + \frac{\Delta t}{\gamma^2} \Lambda_h^{k,n+1} L_h^{k,n+1*} C_h^* (y_\gamma^{n+1} - C_h \hat{z}_h^{k,n+1-}), \\
\hat{z}_h^{k,n+1} = \hat{z}_h^{k,n+1-} + L_h^{k,n+1} (\hat{\theta}_h^{k,n+1} - \hat{\theta}_h^{k,n}).
\end{cases}
$$

Here again, the approximation relies on the fact that the sensitivity and the covariances are computed sequentially from the current estimation $\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n}$.

**Remark 4** Note that, in our work, we omit the effect of numerical errors and, in a larger scope, modelling errors. Taking these into account could typically be done by adding, in addition to the parameter-reduced optimal Kalman filter, a state observer. Due to the large dimension of the (discrete) state space – inherent to any finite element discretizations of problem (6) – Kalman filtering is prohibitive due to the dense Ricatti operator. This is often referred to as the “curse of dimensionality” in the literature. To overcome this difficulty, in the specific case of wavelike problems, a potential strategy proposed in [32, 24] is to resort to a Luenberger observer. This filtering strategy bears significantly lower computational cost compared to Kalman filters, and could be incorporated in our filtering procedure. Note however that, by doing so, a fine balance between parameter sensitivity and corrections of numerical or modelling errors should then be sought.

### 3.3. From a square root to a tangent-free formulation

We are now ready to formulate a practical implementation of our DT-iROEKF. In this respect, we recast the DT-iROEKF formulation into a so-called square-root form as commonly done in data assimilation for the robust implementation of the Kalman approaches [41]. Then we show how we can replace, up to second-order terms, the tangent computations by finite difference scheme over wave solution, hence leading to an original tangent-free approach.
The SEEK algorithm for parameter identification. In this section, we start by formulating a square-root approach of the reduced-order Extended Kalman filter that we are going to apply to our state-parameter decomposition. This strategy is an adaption to parameter estimation of the square root algorithm initially presented in [39]. Following their strategy, we choose to store not $L_h^{k,n}$ but $S_h^{k,n} := L_h^{k,n} D_h^{k,n}$ and $S_h^{k,n+1-} := L_h^{k,n+1} D_h^{k,n}$ with $D_h^{k,n}$ a Cholesky factorization of $\Lambda_h^{k,n}$. Therefore, by equivalently representing $D_h^{k,n}$ as a collection of $N_p$ members

$$ \{D_{h,j}^{k,n}\}_{j=1}^{N_p} \in \mathcal{P}_h, $$

that can be seen as the “columns” of the matrix representations of the operator, each member is then comparable to a parameter increment around the observer trajectory $\hat{\theta}_h^{k,n}$. Hence, each member is associated with a reconstructed displacement

$$ \hat{\mathbf{u}}_{0h}^{k,n} = \sum_{1 \le i \le N_p} (\hat{\theta}_h^{k,n} + D_{h,i}^{k,n}) \varphi_{h,i}, $$

with $\varphi_{h,i}$ the normalized eigenvectors of the operator $A_{0h}$. In the same manner $S_h^{k,n}$ and $S_h^{k,n+1-}$ can be equivalently represented as a collection of $N_p$ members

$$ \{S_{h,j}^{k,n}\}_{j=1}^{N_p} \in \mathcal{Z}_h. $$

Each member is comparable to the variation of a wave solution for the specific parameter increment $D_{h,j}^{k,n}$. Let us now show that a $S_h^{k,n}$ and $D_h^{k,n}$ can be sequentially computed instead of $L_h^{k,n}$ and $\Lambda_h^{k,n}$

**Theorem 3** For all $k \ge 1$ From $S_h^{k,n} = 0$ and $D_h^{k,0} := \text{chol}(\Lambda_h^{k,0})$, we compute sequentially for all $n \in [0:N-1]$,

— Prediction / Forecast:

$$ S_{h,j}^{k,n+1-} := \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) S_{h,j}^{k,n} + B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n}) D_{h,j}^{k,n}, \quad \forall j \in [1:N_p]. \quad (33) $$

— Correction / Analysis: introducing the Gramian operator $G^{k,n+1}$ and $Q^{k,n+1}$ the Cholesky decomposition of its inverse by

$$ G^{k,n+1} := \text{Id}_{\mathcal{P}_h} + \frac{\Delta t}{\gamma^2} S_h^{k,n+1-*} C_h^* C_h S_h^{k,n+1-}, \quad (G^{k,n+1})^{-1} = Q^{k,n+1} Q^{k,n+1*}, $$

and defining

$$ D_h^{k,n+1} := D_h^{k,n} Q^{k,n+1}, \quad \text{and} \quad S_h^{k,n+1} := S_h^{k,n+1-} Q^{k,n+1}. \quad (34) $$

Then, we have the following identities :

$$ \forall n \in [0:N]. \quad \Lambda_h^{k,n} = D_h^{k,n} D_h^{k,n*}, \quad \text{and} \quad S_h^{k,n} = L_h^{k,n} D_h^{k,n}. \quad (35) $$

Moreover, we have

$$ \forall n \in [0:N-1], \quad S_h^{k,n+1-} = L_h^{k,n+1} D_h^{k,n}. \quad (36) $$
**Proof:** We proceed by induction. At $n = 0$ $D_h^{k,0} := \text{chol}(\Lambda_h^{k,0})$ and $L_h^{k,0} = 0$ implies (35). Assuming that (35) is valid at iteration $n$, we have for all $j \in [1:N_p]$,

$$ 
\begin{aligned}
S_{h,j}^{k,n+1-} &= \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) S_{h,j}^{k,n} + B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n+1-}) D_{h,j}^{k,n} \\
&= \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) L_{h,j}^{k,n} D_{h,j}^{k,n} + B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n+1-}) D_{h,j}^{k,n} \\
&= L_{h,j}^{k,n+1} D_{h,j}^{k,n},
\end{aligned}
$$

so that (36) is satisfied. Then from Woodbury inversion formula, we have that

$$ 
\begin{aligned}
(G^{k,n+1})^{-1} &= (\text{Id}_{P_h} + \frac{\Delta t}{\gamma^2} S_h^{k,n+1-*} C_h^* C_h S_h^{k,n+1-})^{-1} \\
&= \text{Id}_{P_h} - \Delta t S_h^{k,n+1-*} C_h^* [\gamma^2 \text{Id}_y + \Delta t C_h S_h^{k,n+1-} S_h^{k,n+1-*} C_h^*]^{-1} C_h S_h^{k,n+1-} \\
&= Q^{k,n+1} Q^{k,n+1*}.
\end{aligned}
$$

Therefore, definition of $D_h^{k,n+1}$ yields

$$ 
\begin{aligned}
D_h^{k,n+1} D_h^{k,n+1*} &= D_h^{k,n} Q^{k,n+1} Q^{k,n+1*} D_h^{k,n} \\
&= \Lambda_h^{k,n} \\
&\quad - \Delta t \Lambda_h^{k,n} L_h^{k,n+1*} C_h^* [\gamma^2 \text{Id}_y + \Delta t C_h L_h^{k,n+1} \Lambda_h^{k,n} L_h^{k,n+1*} C_h^*]^{-1} C_h L_h^{k,n+1} \Lambda_h^{k,n},
\end{aligned}
$$

which, with a second use of Woodbury inversion formula gives $D_h^{k,n+1} D_h^{k,n+1*} = \Lambda_h^{k,n+1}$. Finally, we have that

$$ S_h^{k,n+1} = S_h^{k,n+1-} Q^{k,n+1} = L_h^{k,n+1} D_h^{k,n} Q^{k,n+1} = L_h^{k,n+1} D_h^{k,n+1}. $$

which concludes the proof. ■

Note that in practice, one can use the component-wise expression of the Gramian matrix, namely

$$ G_{ij}^{k,n+1} = \delta_{ij} + \frac{\Delta t}{\gamma^2} (C_h S_{h,i}^{k,n+1-}, C_h S_{h,j}^{k,n+1-})_y, \quad \forall 1 \le i, j \le N_p. $$

In the same fashion, we have

$$ S_{h,j}^{k,n+1} = \sum_{1 \le i \le N_p} Q_{ij}^{k,n+1} S_{h,i}^{k,n+1-}, \quad \text{and} \quad D_{h,i}^{k,n+1} = \sum_{1 \le i \le N_p} Q_{ij}^{k,n+1} D_{h,i}^{k,n}. $$

From the previous theorem, we then easily verify that the DT-iROEKF estimator can be computed using only $S_h^{k,n-}$, $S_h^{k,n}$ and $D_h^{k,n}$ with

— Prediction / Forecast:

$$ \hat{z}_h^{k,n+1-} = \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) \hat{z}_h^{k,n} + r_{h,\Delta t}^{n+1}, \quad (37) $$

— Correction / Analysis:

$$ 
\left\{ 
\begin{aligned}
\hat{\theta}_h^{k,n+1} &= \hat{\theta}_h^{k,n} + \frac{\Delta t}{\gamma^2} D_h^{k,n+1} S_h^{k,n+1*} C_h^* (y_\gamma^{n+1} - C_h \hat{z}_h^{k,n+1-}), \\
\hat{z}_h^{k,n+1} &= \hat{z}_h^{k,n+1-} + \frac{\Delta t}{\gamma^2} S_h^{k,n+1} S_h^{k,n+1*} C_h^* (y_\gamma^{n+1} - C_h \hat{z}_h^{k,n+1-}).
\end{aligned}
\right. 
\quad (38) $$

And again in practice, evaluating (38) is performed using the component-wise expression

$$ 
\left\{ 
\begin{aligned}
\hat{\theta}_h^{k,n+1} &= \hat{\theta}_h^{k,n} + \frac{\Delta t}{\gamma^2} \sum_{1 \le j \le N_p} (C_h S_{h,j}^{k,n+1}, y_\gamma^{n+1} - C_h \hat{z}_h^{k,n+1-})_y D_{h,j}^{k,n+1}, \\
\hat{z}_h^{k,n+1} &= \hat{z}_h^{k,n+1-} + \frac{\Delta t}{\gamma^2} \sum_{1 \le j \le N_p} (C_h S_{h,j}^{k,n+1}, y_\gamma^{n+1} - C_h \hat{z}_h^{k,n+1-})_y S_{h,j}^{k,n+1}.
\end{aligned}
\right. 
$$
A tangent-free alternative. Ultimately, (33) can be seen as the tangent equation of

$$
\begin{align*}
S_{h,j}^{k,n+1-} &= \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) S_{h,j}^{k,n} + B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n}) D_{h,j}^{k,n} \\
&= \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n})(\hat{z}_h^{k,n} + S_{h,j}^{k,n}) + B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n}) D_{h,j}^{k,n} - \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) \hat{z}_h^{k,n} \\
&= \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n} + D_{h,j}^{k,n})(\hat{z}_h^{k,n} + S_{h,j}^{k,n}) - \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n}) \hat{z}_h^{k,n} + O(\|D_{h,j}^{k,n}\|^2) \\
&= \Phi_{h,\Delta t}(\hat{\theta}_h^{k,n} + D_{h,j}^{k,n})(\hat{z}_h^{k,n} + S_{h,j}^{k,n}) + r_{h,\Delta t}^{n+1} - \hat{z}_h^{k,n+1-} + O(\|D_{h,j}^{k,n}\|^2).
\end{align*}
$$

Therefore, by computing the forecast of each member $S_{h,j}^{k,n+1-}$ by a finite difference between two wave equations, we do not need to compute the operator $B_{h,\Delta t}(\hat{\theta}_h^{k,n}, \hat{z}_h^{k,n+1})$. This means that the forecast step can be performed by simply forecasting $\hat{z}_h^{k,n}$ and each member $S_{h,j}^{k,n}$.

### 3.4. The Reduced-Order UKF alternative

In the DT-iROEKF approach, we have seen that we can propagate the estimator $\hat{\theta}_h^{k,n}$, $\hat{z}_h^{k,n}$ and the increment set $\{S_{h,j}^{k,n+1-}\}_{j=1}^{N_p}$ and $\{D_{h,j}^{k,n+1-}\}_{j=1}^{N_p}$. These increments give the sensitivity direction around the estimator $\hat{\theta}_h^{k,n}$, $\hat{z}_h^{k,n}$. Revisiting [36, 31], the Reduced-Order Unscented Kalman alternative (ROUKF) to ROEKF, we are going to replace the computation of the increment sets by sets of so-called sigma-points around $\hat{\theta}_h^{k,n}$ and $\hat{z}_h^{k,n}$ so that the estimator is computed by an averaging formula like in a finite difference stencil. In essence, while EKF propagates a point-estimator and $N_p$ directions around this point-estimator, UKF computes a stencil of sigma-points such that the point-estimator is a resulting average of this point. In practice we will need at least $N_s = N_p + 1$ sigma-points when using a simplex stencil, but more complex stencils with $N_s \ge N_p + 1$ can also be proposed with additional trajectories to be computed. Note that our presentation clarifies the generalization of ROUKF to general complex stencil initially proposed in [31, 33]. This may be of particular interest since increasing the number of sigma points results in a potentially more precise sampling and rendering of the sensitivity directions around the estimator.

*Empirical mean and empirical covariance.* Let us introduce a set of strictly positive weights $\beta = \{\beta_j\}_{j=1}^{N_s} \in (\mathbb{R}_+^*)^{N_s}$ such that

$$ \sum_{1 \le j \le N_s} \beta_j = 1. $$

From this set of weights, we introduce the strictly positive diagonal matrix $D_\beta$ of dimension $N_s \times N_s$ such that

$$ \forall i, j \in [0 : N_s], \quad (D_\beta)_{ij} = \beta_i \delta_{ij}, $$

and the following weight averaging operator defined, for all sample set $\{s_j\}_{j=1}^{N_s} \in \mathcal{X}^{N_s}$ of elements of a Hilbert space $\mathcal{X}$ (typically $\mathcal{Z}_h$, $\mathcal{P}_h$ or $\mathcal{Y}$), by

$$ \mathbb{E}_{\beta}(\{s_j\}_{j=1}^{N_s}) = \sum_{1 \le j \le N_s} \beta_j s_j \in \mathcal{Z}. $$
This can be understood as an empirical mean of the {$s_j$}$_{j=1}^{N_s}$. In the same fashion, we define the empirical covariance as the following linear operator

$$
\mathbb{P}_{\beta}(\{s_j\}_{j=1}^{N_s}) : \mathcal{X} \to \mathcal{X} \\
s \mapsto \sum_{1 \le j \le N_s} \beta_j (s, s_j - \mathbb{E}_{\beta}(\{s_j\}_{j=1}^{N_s}))_{\mathcal{X}} (s_j - \mathbb{E}_{\beta}(\{s_j\}_{j=1}^{N_s})),
$$

where $(\cdot, \cdot)_X$ is the scalar product in $\mathcal{X}$. As soon as we are considering a finite-
dimensional Euclidian space $\mathcal{X} = \mathbb{R}^{N_x}$, then

$$
\mathbb{P}_{\beta}(\{s_j\}_{j=1}^{N_s}) = \sum_{1 \le j \le N_s} \beta_j (s_j - \mathbb{E}_{\beta}(\{s_j\}_{j=1}^{N_s}))(s_j - \mathbb{E}_{\beta}(\{s_j\}_{j=1}^{N_s}))^{\top}.
$$

Sampling around the observer trajectory. Let us now introduce a scaling factor $\rho \in \mathbb{R}_+^*$ and define a particular finite sequence of sigma-points {$e_j$}$_{j=1}^{N_s}$ in the parameter space $\mathcal{P}_h$ such that

$$
\mathbb{E}_{\beta}(\{e_j\}_{j=1}^{N_s}) = 0 \quad \text{and} \quad \mathbb{P}_{\beta}(\{e_j\}_{j=1}^{N_s}) = \text{Id}_{\mathcal{P}_h}. \tag{39}
$$

Note that in this case, if all the sigma-points are roughly at a distance $\delta$ from 0, then

$$
N_p = \mathrm{Tr}(\mathbb{P}_{\beta}(\{s_j\}_{j=1}^{N_s})) \approx \delta \sum_{1 \le j \le N_s} \beta_i = \delta.
$$

Therefore, $\|e_j\|_P = O(N_p^{-1/2})$ which implies that the sigma-points can be relatively far
from 0 as the parameter dimension increases.

Again, defining $D_h^{k,n}$ such as $\Lambda_h^{k,n} = D_h^{k,n} D_h^{k,n*}$ and using $L_h^{k,n}$ at iteration $n$, we define the *sampling* operation as constructing the following set of trajectories

$$
\forall j \in [1:N_s], \quad \left\{
\begin{array}{l}
\hat{z}_{h,j}^{k,n} = \hat{z}_h^{k,n} + \rho L_h^{k,n} D_h^{k,n} e_j, \\
\hat{\theta}_{h,j}^{k,n} = \hat{\theta}_h^{k,n} + \rho D_h^{k,n} e_j,
\end{array}
\right.
\qquad (40)
$$

such that, by construction,

$$
\mathbb{E}_{\beta}(\{\hat{z}_{h,j}^{k,n}\}_{j=1}^{N_s}) = \hat{z}_h^{k,n} \quad \text{and} \quad \mathbb{E}_{\beta}(\{\hat{\theta}_{h,j}^{k,n}\}_{j=1}^{N_s}) = \hat{\theta}_h^{k,n},
$$

while

$$
\frac{1}{\rho^2} \mathbb{P}_{\beta}(\{\hat{z}_{h,j}^{k,n}\}_{j=1}^{N_s}) = L_h^{k,n} \Lambda_h^{k,n} L_h^{k,n*} \quad \text{and} \quad \frac{1}{\rho^2} \mathbb{P}_{\beta}(\{\hat{\theta}_{h,j}^{k,n}\}_{j=1}^{N_s}) = \Lambda_h^{k,n}.
$$

Different kinds of sigma-points can be used to perform this sampling operation – see for
instance [31]. In Figure 1 we illustrate three different (albeit standard) sets of sigma-
points, and we compare it with the sampling strategy derived from the tangent-free
EKF approach presented in the previous section.

*Prediction step of the reduced-order UKF.* The prediction step of the reduced-order UKF is simply obtained by propagating $\hat{z}_{h,j}^{k,n}$ using the discretized dynamics, namely

$$
\hat{z}_{h,j}^{k,n+1-} = \Phi_{h,\Delta t}(\hat{\theta}_{h,j}^{k,n}) \hat{z}_{h,j}^{k,n} + r_{h,\Delta t}^{n+1}
$$

Using the empirical average operator, we can derive the state and parameter sensitivity
around the estimator trajectory, namely we define the predictions

$$
\forall j \in [1:N_s], \quad \left\{
\begin{array}{l}
\sum_{h,j}^{k,n+1-} = \hat{z}_{h,j}^{k,n+1-} - \mathbb{E}_{\beta}(\{\hat{z}_{h,j}^{k,n+1-}\}_{j=1}^{N_s}), \\
\Delta_{h,j}^{k,n} = \hat{\theta}_{h,j}^{k,n} - \mathbb{E}_{\beta}(\{\hat{\theta}_{h,j}^{k,n}\}_{j=1}^{N_s}).
\end{array}
\right.
\qquad (41)
$$
<figure><img src="image_1.png" /><figcaption>Figure 1: Definition of standard sigma-points for the simple case of $\mathcal{P} = \mathbb{R}^2$ and comparison with the state and sensitivity trajectories associated with the tangent-free EKF approach.</figcaption></figure>

The next lemma shows that the quantities defined in (41) are actually consistent with $S_h^{k,n+1-}$ and $D_h^{k,n}$ except that $N_s > N_p$.

**Lemma 4** The predicted state sensitivity defined in (41) satisfies

$$ \forall j \in [1:N_s], \quad \sum_{h,j}^{k,n+1-} = \rho L_h^{k,n+1} D_h^{k,n} e_j + O(\rho^2), \quad (42) $$

while the parameter sensitivity satisfies

$$ \forall j \in [1:N_s], \quad \Delta_{h,j}^{k,n} = \rho D_h^{k,n} e_j, \quad (43) $$

**Proof:** From its definition, and using the definition of the sampling operation (40), one can easily verify the expression satisfied by the parameter sensitivity. Concerning the state sensitivity, on the one hand, the state prediction of each sigma-point satisfies for any $j \in [1:N_s]$,

$$ \begin{aligned} \hat{z}_{h,j}^{k,n+1-} &= \Phi_{h,\Delta t}(\hat{\theta}_{h,j}^{k,n}) \hat{z}_{h,j}^{k,n} + r_{h,\Delta t}^{n+1} \\ &= \Phi_{h,\Delta t}(\hat{\theta}_{h}^{k,n}) \hat{z}_{h}^{k,n} + r_{h,\Delta t}^{n+1} + \Phi_{h,\Delta t}(\hat{\theta}_{h}^{k,n})(\hat{z}_{h,j}^{k,n} - \hat{z}_{h}^{k,n}) \\ &\quad + B_{h,\Delta t}(\hat{\theta}_{h}^{k,n}, \hat{z}_{h}^{k,n})(\hat{\theta}_{h,j}^{k,n} - \hat{\theta}_{h}^{k,n}) + O(\rho^2). \end{aligned} $$

On the other hand, we can introduce

$$ \hat{z}_{h}^{k,n+1-} = \mathbb{E}_{\beta}(\hat{z}_{h,j}^{k,n+1-}) = \mathbb{E}_{\beta}(\{\Phi_{h,\Delta t}(\hat{\theta}_{h,j}^{k,n})\hat{z}_{h,j}^{k,n} + r_{h,\Delta t}^{n+1}\}_{j=1}^{N_s}), $$

and a Taylor expansion gives

$$ \hat{z}_{h}^{k,n+1-} = \Phi_{h,\Delta t}(\hat{\theta}_{h}^{k,n}) \hat{z}_{h}^{k,n} + r_{h,\Delta t}^{n+1} + O(\rho^2). $$

Using the sampling operations defined in (40) we obtain

$$ \begin{aligned} \sum_{h,j}^{k,n+1-} &= \hat{z}_{h,j}^{k,n+1-} - \hat{z}_{h}^{k,n+1-} \\ &= \Phi_{h,\Delta t}(\hat{\theta}_{h}^{k,n})(\hat{z}_{h,j}^{k,n} - \hat{z}_{h}^{k,n}) + B_{h,\Delta t}(\hat{\theta}_{h}^{k,n}, \hat{z}_{h}^{k,n})(\hat{\theta}_{h,j}^{k,n} - \hat{\theta}_{h}^{k,n}) + O(\rho^2) \\ &= \rho(\Phi_{h,\Delta t}(\hat{\theta}_{h}^{k,n})L_{h}^{k,n} + B_{h,\Delta t}(\hat{\theta}_{h}^{k,n}, \hat{z}_{h}^{k,n}))D_{h}^{k,n}e_j + O(\rho^2) \\ &= \rho L_{h}^{k,n+1}D_{h}^{k,n}e_j + O(\rho^2). \end{aligned} $$
<figure><img src="image_2.png" /></figure>

In the following, we denote by $\Sigma_h^{k,n+1-} : \mathbb{R}^{N_s} \to \mathcal{Z}$ the operator defined by

$$ \Sigma_h^{k,n+1-} : \mathbb{R}^{N_s} \ni \alpha \mapsto \sum_{1 \le i \le N_s} \alpha_i \Sigma_{h,i}^{k,n+1-} \in \mathcal{Z} $$

and similarly for $\Delta_h^{k,n}$. Namely a matrix representation of $\Sigma_h^{k,n+1-}$ and $\Delta_h^{k,n}$ corresponds to the column concatenation of the state and parameter sensitivity as

$$ \Sigma_h^{k,n+1-} = (\Sigma_{h,1}^{k,n+1-} | \cdots | \Sigma_{h,N_s}^{k,n+1-}), \quad \Delta_h^{k,n} = (\Delta_{h,1}^{k,n} | \cdots | \Delta_{h,N_s}^{k,n}), \quad (44) $$

which are matrices with dimensions $\dim(\mathcal{Z}_h) \times N_s$ and $N_p \times N_s$ respectively. These matrices satisfy the following properties, which will enable us to derive the correction step of the reduced-order UKF.

**Lemma 5** The state sensitivity matrix defined in (44) satisfies

$$ \frac{1}{\varrho^2} \Sigma_h^{k,n+1-} D_\beta \Sigma_h^{k,n+1-*} = L_h^{k,n+1} \Lambda_h^{k,n} L_h^{k,n+1*} + O(\varrho), $$

while the parameter sensitivity matrix satisfies

$$ \frac{1}{\varrho^2} \Delta_h^{k,n} D_\beta \Delta_h^{k,n*} = \Lambda_h^{k,n}. $$

Both matrices are linked through the identity

$$ \frac{1}{\varrho^2} \Sigma_h^{k,n+1-} D_\beta \Delta_h^{k,n*} = L_h^{k,n+1} \Lambda_h^{k,n} + O(\varrho). $$

**Proof:** Starting with the parameter sensitivity matrix, we have

$$ \Delta_h^{k,n} D_\beta \Delta_h^{k,n*} = \sum_{j=1}^{N_s} \beta_j \Delta_{h,j}^{k,n} \Delta_{h,j}^{k,n*} = \varrho^2 \sum_{j=1}^{N_s} \beta_j (D_h^{k,n} e_j)(D_h^{k,n} e_j)^* = D_h^{k,n} (\varrho^2 \sum_{j=1}^{N_s} \beta_j e_j e_j^*) D_h^{k,n*}. $$

Furthermore, using the properties (39) satisfied by the sigma-points combined with the identity $D_h^{k,n} D_h^{k,n*} = \Lambda_h^{k,n}$, yields the desired results. Following the same line of arguments, one has

$$ \Sigma_h^{k,n+1-} D_\beta \Sigma_h^{k,n+1-*} = L_h^{k,n+1} D_h^{k,n} (\varrho^2 \sum_{j=1}^{N_s} \beta_j e_j e_j^*) D_h^{k,n*} L_h^{k,n+1*} + O(\varrho^3), \\ = \varrho^2 L_h^{k,n+1} \Lambda_h^{k,n} L_h^{k,n+1*} + O(\varrho^3). $$

For the last relation, we combine (42) and (43) to obtain $\Sigma_h^{k,n+1-} = L_h^{k,n+1} \Delta_h^{k,n} + O(\varrho^2)$, and conclude using $\frac{1}{\varrho^2} \Delta_h^{k,n} D_\beta \Delta_h^{k,n*} = \Lambda_h^{k,n}$. ■
Correction step of the reduced-order UKF. We can now provide the main result for
deriving the correction step of the reduced-order UKF observer.

**Theorem 6** For all $k \ge 1$ and $n \in [0:N-1]$, defining the Gramian operator by

$$
\tilde{G}^{k,n+1} = D_{\beta}^{-1} + \frac{\Delta t}{\gamma^2} \Sigma_h^{k,n+1-*} C_h^* C_h \Sigma_h^{k,n+1-},
$$

and an alternative Cholesky decomposition of its inverse by

$$
(\tilde{G}^{k,n+1})^{-1} = \tilde{Q}^{k,n+1} D_{\beta} \tilde{Q}^{k,n+1*},
$$

along with $\Delta_h^{k,n+1} = \Delta_h^{k,n} \tilde{Q}^{k,n+1}$, yields

$$
\Lambda_h^{k,n+1} = \frac{1}{\varrho^2} \Delta_h^{k,n+1} D_\beta \Delta_h^{k,n+1*} + O(\varrho). \quad (45)
$$

**Proof:** To start with, let us remark that applying Woodbury inversion formula yields

$$
\begin{align*}
(G^{k,n+1})^{-1} &= (D_{\beta}^{-1} + \frac{\Delta t}{\gamma^2} \Sigma_h^{k,n+1-*} C_h^* C_h \Sigma_h^{k,n+1-})^{-1} \\
&= D_{\beta} - \Delta t D_{\beta} \Sigma_h^{k,n+1-*} C_h^* (\gamma^2 \mathrm{Id}_y + \Delta t C_h \Sigma_h^{k,n+1-} D_{\beta} \Sigma_h^{k,n+1-*} C_h^*)^{-1} C_h \Sigma_h^{k,n+1-} D_{\beta} \\
&= \tilde{Q}^{k,n+1} D_{\beta} \tilde{Q}^{k,n+1*}.
\end{align*}
$$

By definition, we have

$$
\Delta_h^{k,n+1} D_\beta \Delta_h^{k,n+1*} = \Delta_h^{k,n} \tilde{Q}^{k,n+1} D_\beta \tilde{Q}^{k,n+1*} \Delta_h^{k,n*},
$$

which, using the relations in Lemma 5, entails

$$
\frac{1}{\varrho^2} \Delta_h^{k,n+1} D_\beta \Delta_h^{k,n+1*} = \Lambda_h^{k,n} \\
- \Delta t \Lambda_h^{k,n} L_h^{k,n+1*} C_h^* [\gamma^2 \mathrm{Id}_y + \Delta t C_h L_h^{k,n+1} \Lambda_h^{k,n} L_h^{k,n+1*} C_h^*]^{-1} C_h L_h^{k,n+1} \Lambda_h^{k,n} + O(\varrho),
$$

$$
\text{which gives us the desired result } \frac{1}{\varrho^2} \Delta_h^{k,n+1} D_\beta \Delta_h^{k,n+1*} = \Lambda_h^{k,n+1} + O(\varrho).
$$

<figure><img src="image_3.png" /></figure>

Then, following the definition of the correction (or analysis) step of the SEEK algorithm,
namely (34), we define the corrected stat sensitivity operator as $\Sigma_h^{k,n+1} = \Sigma_h^{k,n+1-} \tilde{Q}^{k,n+1}$.
Component-wise, the state and parameter sensitivity reads

$$
\Sigma_{h,j}^{k,n+1} = \sum_{1 \le i \le N_s} \tilde{Q}_{ij}^{k,n+1} \Sigma_{h,i}^{k,n+1-} \quad \text{and} \quad \Delta_{h,j}^{k,n+1} = \sum_{1 \le i \le N_s} \tilde{Q}_{ij}^{k,n+1} \Delta_{h,i}^{k,n},
$$

where for all $n$, $\tilde{Q}^{k,n} \in \mathbb{M}_{N_s}(\mathbb{R})$ is given by the alternative Cholesky decomposition of
$(\tilde{G}^{k,n})^{-1} = \tilde{Q}^{k,n} D_{\beta} \tilde{Q}^{k,n*}$, and $G^{k,n}$ reads

$$
\tilde{G}_{ij}^{k,n} = \beta_i^{-1} \delta_{ij} + \frac{\Delta t}{\varrho^2 \gamma^2} (C_h \Sigma_{h,i}^{k,n-} , C_h \Sigma_{h,j}^{k,n-}) y, \quad \forall i,j \in [1:N_s].
$$

The corrected estimator then takes the following form

$$
\left\{
\begin{aligned}
\hat{\theta}_h^{k,n+1} &= \hat{\theta}_h^{k,n} + \frac{\Delta t}{\varrho^2 \gamma^2} \sum_{1 \le j \le N_s} (C \Sigma_{h,j}^{k,n+1}, y_\gamma^{n+1} - C \hat{z}_h^{k,n+1-}) y \Delta_{h,j}^{k,n+1}, \\
\hat{z}_h^{k,n+1} &= \hat{z}_h^{k,n+1-} + \frac{\Delta t}{\varrho^2 \gamma^2} \sum_{1 \le j \le N_s} (C \Sigma_{h,j}^{k,n+1}, y_\gamma^{n+1} - C \hat{z}_h^{k,n+1-}) y \Sigma_{h,j}^{k,n+1}.
\end{aligned}
\right.
$$
Resampling around the observer trajectory As a final step of our iterated UKF algorithm, we derive relations that enable us to rewrite the sampling operation (40) exclusively in terms of the computed state and parameter sensitivities for the subsequent time step n + 1.

**Lemma 7** The corrected state and parameter sensitivity matrices satisfy

$$ \forall j \in [1:N_s], \quad \left\{ \begin{array}{l} \sum_{h,j}^{k,n+1} = \varrho L_h^{k,n+1} D_h^{k,n+1} e_j + O(\varrho^2), \\ \Delta_{h,j}^{k,n+1} = \varrho D_h^{k,n+1} e_j + O(\varrho^2). \end{array} \right. \qquad (46) $$

**Proof:** Using the definition of $\sum_{h}^{k,n+1}$ and the relation (42) we obtain

$$ \begin{aligned} \sum_{h,j}^{k,n+1} &= \sum_{1 \le i \le N_s} \tilde{Q}_{ij}^{k,n+1} \sum_{h,i}^{k,n+1-} = \varrho \sum_{1 \le i \le N_s} \tilde{Q}_{ij}^{k,n+1} L_h^{k,n+1} D_h^{k,n} e_i + O(\varrho^2) \\ &= \varrho L_h^{k,n+1} \left[ \sum_{1 \le i \le N_s} \tilde{Q}_{ij}^{k,n+1} D_h^{k,n} e_i \right] + O(\varrho^2). \\ &= L_h^{k,n+1} \left[ \sum_{1 \le i \le N_s} \tilde{Q}_{ij}^{k,n+1} \Delta_{h,i}^{k,n} \right] + O(\varrho^2). \\ &= L_h^{k,n+1} \Delta_{h,j}^{k,n+1} + O(\varrho^2) = \varrho L_h^{k,n+1} D_h^{k,n+1} e_j + O(\varrho^2). \end{aligned} $$

In the same manner, using the definition of $\Delta_{h}^{k,n+1}$ and the relation (43) yields
$$ \Delta_{h,j}^{k,n+1} = \varrho D_h^{k,n+1} e_j + O(\varrho^2). \quad \blacksquare $$

Hence, $\sum_{h,j}^{k,n+1}$ and $\Delta_{h,j}^{k,n+1}$ are consistent with $S_h^{k,n+1}$ and $D_h^{k,n+1}$, up to second order terms, as for the prediction step. Once again, the only – but fundamental – difference compared to the EKF version is the fact that the Gramian is defined in $M_{N_s}(\mathbb{R})$ instead of $M_{N_p}(\mathbb{R})$. Note that, one important aspect of the relations (46) is that it enables us to recreate the sample at n + 1 by discarding second-order terms in the following relations

$$ \hat{z}_{h,j}^{k,n+1} = \hat{z}_h^{k,n+1} + \sum_{h,j}^{k,n+1} \pm O(\varrho^2), \quad \text{and} \quad \hat{\theta}_{h,j}^{k,n+1} = \hat{\theta}_h^{k,n} + \Delta_{h,j}^{k,n+1} \pm O(\varrho^2). $$

We provide in Figure 2 an illustration of its main steps, in order to fully apprehend the various unknowns appearing in the algorithm.

**Final algorithm.** We can now summarize the complete algorithm used in this work, here named the Discrete-Time iterated Reduced-Order Unscented Kalman Filter (DT-IROUKF):

* Choice of the UKF sigma-points $\{e_j\}_{j=1}^{N_s}$ and the associated weights $\{\beta_j\}_{j=1}^{N_s}$.

* Initialization of the parameter prior from the modal decomposition

$$ \hat{u}_0(\mathbf{x}) = \sum_{1 \le j \le N_p} \hat{\theta}_{0h,j} \varphi_j(\mathbf{x}). $$

* For all $k \ge 1$ until convergence:

* Initialization at $n = 0$:
<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of the main steps of the reduced-order UKF algorithm in the specific case of $\mathcal{P} = \mathbb{R}^2$ and focusing specifically on the state estimation part.</figcaption></figure>

Parameter initialization:

$$ \hat{\theta}_{h}^{k,0} = \hat{\theta}_{h}^{k-1,N} \quad \text{if } k > 1, \quad \hat{\theta}_{h}^{1,0} = \hat{\theta}_{0h} \quad \text{otherwise.} $$

Parameter and state particles initialization:

$$ \forall j \in [1:N_s], \quad \left\{ \begin{array}{l} \hat{\theta}_{h,j}^{k,0} = \hat{\theta}_{h}^{k,0} + \varrho D^{k,0} e_j, \\ \hat{z}_{h,j}^{0} = 0, \end{array} \right. $$

with $D^{k,0} = \text{chol}(\Lambda_0^m)$.

For $0 \le n \le N$:

Prediction step:

$$ \forall j \in [1:N_s], \quad \left\{ \begin{array}{l} \hat{z}_{h,j}^{k,n+1-} = \Phi_{h,\Delta t}(\hat{\theta}_{h,j}^{k,n}) \hat{z}_{h,j}^{k,n} + r_{h,\Delta t}^{n+1}, \\ \sum_{h,j}^{k,n+1-} = \hat{z}_{h,j}^{k,n+1-} - \mathbb{E}_{\beta}(\{\hat{z}_{h,j}^{k,n+1-}\}_{j=1}^{N_s}), \\ \Delta_{h,j}^{k,n} = \hat{\theta}_{h,j}^{k,n} - \mathbb{E}_{\beta}(\{\hat{\theta}_{h,j}^{k,n}\}_{j=1}^{N_s}). \end{array} \right. $$

Correction step:

- Compute Gramian matrix and alternative Cholesky decomposition w.r.t. $D_\beta$ of its inverse:

$$ \left\{ \begin{array}{l} G^{k,n+1} = (\beta_i^{-1}\delta_{ij} + \frac{\Delta t}{\gamma^2}(C_h\sum_{h,i}^{k,n+1-}, C_h\sum_{h,j}^{k,n+1-}))y)_{1\le i,j\le N_s}, \\ Q^{k,n+1} = \text{chol}_{D_\beta}((G^{k,n+1})^{-1}). \end{array} \right. $$

- Compute corrected sensitivity matrices:
$$ \forall j \in [1:N_s], \quad \left\{ \begin{array}{l} \Sigma_{h,j}^{k,n+1} = \sum_{1 \le i \le N_p} Q_{ij}^{k,n+1} \Sigma_{h,i}^{k,n+1-}, \\ \Delta_{h,j}^{k,n+1} = \sum_{1 \le i \le N_p} Q_{ij}^{k,n+1} \Delta_{h,i}^{k,n}. \end{array} \right. $$

– Compute corrected state and parameter:

$$ \left\{ \begin{aligned} \hat{\theta}_h^{k,n+1} &= \hat{\theta}_h^{k,n} + \frac{\Delta t}{\gamma^2} \sum_{1 \le j \le N_s} (C \Sigma_{h,j}^{k,n+1}, y_\gamma^{n+1} - C \hat{z}_h^{k,n+1-}) y \Delta_{h,j}^{k,n+1} \\ \hat{z}_h^{k,n+1} &= \hat{z}_h^{k,n+1-} + \frac{\Delta t}{\gamma^2} \sum_{1 \le j \le N_s} (C \Sigma_{h,j}^{k,n+1}, y_\gamma^{n+1} - C \hat{z}_h^{k,n+1-}) y \Delta_{h,j}^{k,n+1} \end{aligned} \right. $$

→ Sampling step:

$$ \forall j \in [1:N_s], \quad \left\{ \begin{array}{l} \hat{\theta}_{h,j}^{k,n+1} = \hat{\theta}_h^{k,n+1} + \Delta_{h,j}^{k,n+1}, \\ \hat{z}_{h,j}^{k,n+1} = \hat{z}^{k,n+1} + \Sigma_{h,j}^{k,n+1}. \end{array} \right. $$

Additionally, for output purposes, we can re-define up to second-order terms the parameter covariance operator from the stored {$\Delta_{h,j}^{k,n}$}$_{j=1}^{N_s}$ and {$\hat{z}_{h,j}^{n+1-}$}$_{j=1}^{N_s}$:

$$ \Lambda_h^{k,n} : \mathcal{P}_h \ni \theta \mapsto \frac{1}{\varrho^2} \sum_{1 \le j \le N_s} \beta_j(\Delta_{h,j}^{k,n}, \theta) \mathcal{P}_h \Delta_{h,j}^{k,n} \in \mathcal{P}_h. $$

In the same fashion, the parameter sensitivity operator reads

$$ L_h^{k,n+1} : \mathcal{P}_h \ni \theta \mapsto \frac{1}{\varrho} \sum_{1 \le j \le N_s} \beta_j((\Lambda_h^{k,n})^{-1}\theta, \Delta_{h,j}^{k,n}) \mathcal{P}_h (\hat{z}_{h,j}^{k,n+1-} - \hat{z}_h^{k,n+1-}) \in \mathcal{Z}_h. $$

**Remark 5 (Scaling factor)** The scaling factor $\varrho$ is here defined to relate the reduced-order EKF with the reduced-order UKF showing that EKF is an asymptotic version of UKF with $\varrho$ tending to 0. But in practice, the UKF algorithm can be used with any $\varrho$ as it was initially designed to avoid relying on tangent operator when finite difference can produced better statistical estimates [25, 31]

Finally let us give some implementation details, see also Figure 3. We implemented the DT-iROUKF algorithm in a specific Python library called AKILLES‡ based on the message passing library $\oslash$MQ§. The objective of this library is to couple the estimator library with any model, written in any language that can receive and send vector quantities through $\oslash$MQ. In AKILLES, the Gramian, Cholesky decomposition, corrected states and parameters are computed by the main process of AKILLES called Agamemnon who communicates through $\oslash$MQ with Achilles' army of Myrmidons – here they represent wave solvers built upon the SFEM kernel of the commercial software CIVA$\|$. The physical outputs are finally processed after the correction step by the process Achilles, a specific Myrmidon wave solver in CIVA in charge of extra output computations. Note that

<footnote>‡ https://gitlab.inria.fr/AKILLES/AKILLES</footnote>

<footnote>§ https://zeromq.org</footnote>

<footnote>|| https://www.extende.com</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Parallel implementation of the DT-iROUKF in AKILLES.</figcaption></figure>

the non-linear mechanical deformation and the corresponding modal decomposition are initially computed using the Finite Element Library MoReFEM¶ for solving large deformation mechanical problems.

## 4. Numerical results

In this section, we illustrate the use of the presented algorithm for estimating the deformation caused by the mechanical loading in structures using GW measurements. Three different cases are presented in increasing order of complexity. The first case is a small-scale problem with the objective of testing the robustness and analyzing the estimation results of the proposed strategies. The second case is an aluminum plate under traction forces inspired from [20], and done in the context of acoustoelasticity to retrieve the material third-order elastic constants. The third case is a pipe subjected to a 4-point bending test inspired from [43], and related to the detection of welding defects under operational conditions employing ultrasound. All the illustrated cases are presented with the structure described below.

*Noisy synthetic data generation.* First, let us describe the numerical parameters for generating noisy synthetic data, and comprising mechanical loading, ultrasonic excitation and signal acquisition. A quasi-static structural problem is defined to compute the structural deformation caused by mechanical loading. For all the cases, the constitutive behavior of the material is modeled by the Compressible Neo Hookean
(CNH) hyperelastic law. With $\lambda$ and $\mu$ being the Lamé parameters, its hyperelastic potential reads

$$ \mathcal{W}^{\text{CNH}} = \frac{\lambda}{2}(\sqrt{I_3} - 1)^2 + \frac{\mu}{2}(I_1 - 3 - \log(I_3)), $$

with $I_1 = \text{tr} \mathbf{C}$ and $I_3 = \det \mathbf{C}$ being invariants of the right Cauchy-Green deformation tensor $\mathbf{C} = \mathbf{F}^\mathrm{T}\mathbf{F}$. The mesh, type of forces and boundary conditions configuring the quasi-static problem are given for every configurations considered in the following. After defining the configuration for the structural deformation problem, we define the configuration for the excitation, propagation and measurement of ultrasonic waves. A force acting upon the outer surface generates a wavefield radially with respect to the structure surface. It represents a force $\mathbf{f}$ defined in a sufficiently thin layer on the outer surface. Its ring geometry has an inner radius of 10mm and an outer radius of 20mm. The excitation signal is a 5-cycle cosine Hanning windowed at a specific frequency for each case. We consider zero initial conditions for the wavefield. The observation data are generated by using the operator in (10) with the domain $\omega_i$ defined upon the outer surface as a sufficiently thin layer in which the wavefield does no vary in the thickness direction. This amounts to modeling a point or surface probe with specific sensitivity. One may note that during the estimation procedure, as the estimated deformation changes while the wave propagation problem is running, the CFL condition (28) also changes. However, the time step is computed once at the initialization and this changes may cause numerical instability. To avoid such an issue, the time step is computed to satisfy the CFL condition for the expected deformation extrema.

To emulate signals obtained from an acquisition system in real scenarios, we generate synthetic observation data using the target deformation and the wave propagation solver. Gaussian noise is added to the simulated data to represent potential electronic and environmental noise. The Gaussian noise is added to the simulated measurements as

$$ y_{\gamma,i}^{n} = \tilde{y}_{i}^{n} + \frac{\gamma}{\sqrt{\Delta t}} \|\tilde{y}_{i}\|_{\ell_{2}} \chi_{i}^{n}, $$

where $\chi_i^n$ is a realization of a random variable with distribution $\mathcal{N}(0, 1)$. The noise level is set by adjusting $\gamma$.

*Reducing the parametric space using mode decomposition.* To represent the parametric space on a reduced basis, we use the eigenmodes of the quasi-static problem as defined in (11). First, we compute a finite set of eigenmodes $\mathcal{I}$ ordered in ascending eigenvalue order, i.e. from the lower to highest spatial frequency. A reduced parametric space for reconstructing the deformation is selected as a subset $\mathcal{I}^* \subset \mathcal{I}$. This selection is made by having a guess $\mathbf{u}_0^*$ of the deformation decomposed in the eigenbasis of $\mathcal{I}$ and selecting the smallest subset $\mathcal{I}^*$ that satisfy a representation error criterion. More precisely, by decomposing $\mathbf{u}_0^*$ in the eigenmodes of $\mathcal{I}$ we select a set of modes associated with the
components $\theta_i^*$ that have a minimum relative $\ell_2$-error $\tau^*$ defined as

$$ \tau^* = \frac{\|\mathbf{u}_0^* - \sum_{j \in \mathcal{I}^*} \theta_j^* \varphi_j\|}{\|\mathbf{u}_0^*\|} $$

To better visualize the relevance of the selected modes we plot the representation error for the first $m$ modes in $\mathcal{I}$, namely $1 - \tau_m$, where

$$ \tau_m = \frac{\|\mathbf{u}_0^* - \sum_{j \in \mathcal{I}_m} \theta_j^* \varphi_j\|}{\|\mathbf{u}_0^*\|} $$

and $\mathcal{I}_m$ being the set of the $m$ lowest frequency modes in $\mathcal{I}$.

**Initializing the estimator.** To initialize the estimator, we must give an a priori parameter and its covariance. Even though one could use $\theta^*$ as an a priori parameter value, we set $\hat{\theta}_h^{1,0}$ to zero in order to test the robustness of our estimator with respect to the initial guess. Regarding the initial covariance, we consider that $\max_{i \in \mathcal{I}^*} \theta_i^*$ is a relevant estimation of the standard deviation of the lowest frequencies, assuming that the lowest frequency is the most relevant for reconstruction. Then, the initial covariance $\Lambda_0$ is computed as, based on (14),

$$ (\Lambda_0)_{ij} = \frac{\lambda_i^2}{\lambda_{\min}^2} (\max_{i \in \mathcal{I}^*} \theta_i^*)^2 \delta_{ij}, \quad \forall i, j \in \mathcal{I}^*, $$

where the subscript min is associated with the lowest-frequency component in $\mathcal{I}^*$. The constant $M$ used to normalize the regularization term is computed as

$$ M = \sum_{i \in \mathcal{I}^*} \frac{\lambda_{\min}^2}{\lambda_i^2}. $$

Finally, the constant $\varepsilon$, which weights the regularization and misfit terms, is adjusted for each case. The parameter $m$ appearing in the parametric space norm is set as $m = 1$.

**Presentation of the estimation results.** At each outer iteration $k$ of our estimator, the last estimated parameter $\hat{\theta}_h^{k-1,N}$ is used as an initial parameter for the next one $\hat{\theta}_h^{k,0}$ while the initial covariance restarts at $\Lambda_0$. We plot for each case the evolution in $n \in [1:N]$ of the estimated parameters $\hat{\theta}_h^{k,n}$ and a region representing an estimation of the associated standard deviation $\sqrt{(\Lambda^{k,n+1})_{ii}}$. Regarding the later, the initial covariance is propagated through the estimation process while retrieving information from the observations. It provides qualitative information about how the uncertainty evolves, although a quantitative interpretation of this uncertainty must be considered carefully.

**Sensitivity analysis.** Instead of using the subset $\mathcal{I}^*$, the estimation could be done using a larger set of eigenmodes, $\mathcal{I}$ for instance, eliminating the need for $\mathbf{u}_0^*$. However, the ill-posedness of the problem combined with the potential instability of the forward model, as the estimated pre-deformation varies in time, makes it unviable. To circumvent this difficulty, we assume to know a set $\mathcal{I}^*$ of eigenmodes that best represents the target deformation. Nonetheless, one may like to gather $\mathcal{I}^*$ directly from the measurements,
instead of using $\mathbf{u}_0^*$. As an attempt to do so, we perform a sensitivity analysis using one Levenberg-Marquardt iteration of the proposed estimator and the set $\mathcal{I}$. As we run the estimator in the whole set $\mathcal{I}$, we must avoid instability by constraining the trajectory $\hat{\theta}_h^{1,n}$. This can be done by setting the covariance as

$$ (\Lambda_0)_{ii} = \gamma^2 = 10^{-8} $$

with a parameter $\varepsilon = 1$, initially. Although at the end of the iteration the estimated $\hat{\theta}_h^{1,N}$ have no quantitative significance, it still retains qualitative information relating the parametric space to the observations. To visualize and assess this relation, we plot the normalized estimated parameters at the end of the estimation,

$$ \frac{\hat{\theta}_{h,i}^{1,N}}{\max_{i \in \mathcal{I}} \hat{\theta}_{h,i}^{1,N}}, \quad \forall i \in \mathcal{I}, $$

as well as the estimated parameters $\hat{\theta}_{h,i}^{1,n}$ over time. Additionally to the estimated values, the covariance matrix is also an output of the estimation process and can be used for this purpose. We plot its normalized inverse,

$$ (\Lambda_0^{-\frac{1}{2}}) ((\Lambda^{0,N})^{-1} - \Lambda_0^{-1}) (\Lambda_0^{-\frac{1}{2}}), $$

representing the normalized Gramian matrix, related to the system's observability [34].

## 4.1. Results and discussion

### 4.1.1. A test case on robustness to noise
To illustrate and test the estimator robustness to noise, we model an aluminum cube of dimensions 60 × 60 × 60 mm³. The configuration for the quasi-static problem is shown in Figure 4a. The target deformation shown in Figure 4b is computed by solving the quasi-static problem for a body force of 50N/mm³ in the Y direction.

<figure><img src="image_6.png" /><figcaption>(a) Configuration of forces the quasi-static (b) Target deformation (visualization scaled problem. 10x).</figcaption></figure>

<caption>Figure 4: Quasi-static problem configuration for the illustration on the cube.</caption>
From zero initial conditions, the wave field is excited at one of its faces at 100kHz in a ring region as described previously. The configuration for the wave propagation model is depicted in Figure 5a, where the characteristic Gauss-Lobatto points of the SFEM discretization are represented. Examples of extracted signals are plotted in Figure 5b. In this case, the observations are point measurements at ten positions, where the three components of the displacement field are obtained. Considering the definition of the observation operator (10), $\omega_i$ is a sufficiently small volume for which the displacement field is considered constant. The three components are taken as sensitivity $d_i$, resulting in $d = 30$, the number of measured signals. Each face, except the excited one, has two measurement points on the surface at a distance of 15 mm from the edges. We run the wave propagation solver up to 100µs to generate the synthetic observation data. In Figure 5b we show the extracted field components for the point located at (15,15,0).

<figure><img src="image_7.png" /><figcaption>(a) Configuration of the wave propagation problem. An example of the wave field is (b) Example of measured signals (three components) at (15,15,0).</figcaption></figure>

Figure 5: Configuration and examples of measured signals for the wave propagation problem. The measurement regions are represented in blue and the excited region in red.

To test the effect of noise in the estimation, different sets of observed data are constructed with different levels of noise $\gamma = [10^{-4}, 0.01, 0.1, 0.2]$. Examples of the different levels of added noise are shown in Figure 6.

The eigenmodes associated with the hundred lowest eigenvalues are computed, embodying $\mathcal{I}$. We decompose a guess of the deformation generated by an axial body force with low amplitude, $u_0^*$, in the basis associated with $\mathcal{I}$. Using this decomposition we plot $1 - \tau_m$ and the associated relevant modes in Figure 7. Then, the most
<figure><img src="image_8.png" /><figcaption>Figure 6: Example of different noise levels added to the synthetic data for the cube case.</figcaption></figure>

relevant modes can be selected and are therefore highlighted, making, for $\tau^* = 0.01$, $\mathcal{I}^* = \{3, 15, 19, 36, 39, 47, 90, 98\}$, meaning that with 8 selected modes we can reconstruct 99% of its deformation.

<figure><img src="image_9.png" /><figcaption>(a) Reconstruction error per number of (b) Deformation associated with the relevant modes. eigenmodes.</figcaption></figure>

<caption>Figure 7: Reconstruction of $\boldsymbol{u}^*$ using the lowest frequency modes from which the selection of the most relevant modes is done.</caption>

We run the estimation algorithm for the different noise levels $\gamma$. We first use $\varepsilon = 1$
as the regularization and misfit terms in function (15) are normalized. The iterative
process is shown in Figure 8 with the time evolution of the estimated parameters $\hat{\theta}_h^k$.
To illustrate the evolution of the estimated deformation, we reconstruct it at given
iterations and time and compare it with the target deformation in Figure 9, for the case
with $\gamma = 0.2$.

*Discussion.* At relatively low noise levels ($\gamma \le 0.01$), the estimate converges in one iteration, and at each subsequent iteration, the estimate deviates slightly from the converged value and converges again. At higher noise levels ($\gamma \ge 0.1$), a few iterations may be required to converge with significant errors w.r.t. the target, although it presents a good estimation of the reconstructed deformation, as shown in Figure 9. In addition, the estimator updates the estimated parameter at a lower rate as $\gamma$ increases, as we can see in the case of $\gamma = 0.2$. This is due to the fact that $\gamma$ weights the mismatch term, as can be seen in (15), and gives less credibility to the observed data when noise increases. In such a case, the parameter $\varepsilon$ can be adjusted to increase the rate at which the estimator updates the parameter based on the observed data. For example, we set $\varepsilon = 0.025$ and show the corresponding estimate in Figure 10. We can see that when the update rate is increased, the noise causes larger instabilities in the estimated parameter during an iteration. The variable $\varepsilon$ can be adjusted on a case-by-case basis. As discussed earlier, the tangent stiffness operator in (6) can be ill-posed, i.e. the coercivity assumption (8) for the forward problem is not guaranteed. Therefore, $\varepsilon$ must be adjusted taking into account the direct-problem stability issues due to strong oscillations during the estimation.

As mentioned earlier, we can restrict the estimation to a reduced parametric space
without resorting to an initial deformation prior, but by retrieving information from the
observations. For this case, we evaluate the possibility of doing this through a sensitivity
analysis using a low noise level, $\gamma = 10^{-4}$, and $\varepsilon = 1$. The normalized estimated
parameters are shown in Figure 11a and the quantity $\mathcal{I}^*$ is highlighted. Based on these
results, the most important modes cannot be completely distinguished, but one could
empirically set a threshold to select them. At the risk of overlooking some important
modes, an iterative selection process can be proposed but it is not discussed here. The
evolution of these parameters during estimation is shown in Figure 11b, with the $\mathcal{I}^*$ set
highlighted. The main modes do not appear to exhibit any particular dynamic behavior
that could be used to distinguish them. Finally, Figure 11c plots the Gramian for the
50 first modes in $\mathcal{I}$. The Gramian quantifies how observable and distinguishable the
different modes are from each other and information about the overall observability of
the inverse problem can be retrieved. High observability can lead to overestimation,
as in the case of the higher frequency modes, for example, mode 47. Inversely, low
observability can lead to underestimation, as in the case of mode 3. The diagonal values
quantifies the observability of the associated modes while off-diagonal terms quantifies
how indistinguishable the associated modes are from each other.
4.1.2. Aluminum plate under traction We present here a case inspired from [20] to model guided wave propagation in axially loaded structures. The specimen is an aluminum plate of dimensions 610 × 305 × 6.35 mm³ under traction forces. The configuration for the quasi-static problem is shown in Figure 12a. We compute the target deformation with an axial stress of 57.5 MPa, resulting in the deformation illustrated in Figure 12b.

The ultrasonic excitation is done as previously described at the center of the plate (upper surface) at 100kHz. The observation data are obtained by defining the observation operator in twelve surface regions $\omega_i$ distributed in an ellipse as depicted in Figure 13a. The mesh used in the wave propagation problem is also depicted. The synthetic observed data are generated by running the simulation up to 200µs. An example of the obtained signal is shown in Figure 13b.

For the larger set of eigenmodes $\mathcal{I}$ we compute those associated with the 120 lowest eigenvalues. As the configuration has no essential boundary conditions, rigid body movements are penalized and the first 6 eigenmodes, associated to them, are removed. From a guess of the deformation $\mathbf{u}_0^*$ generated by the traction surface forces with lower amplitude, we can select a smaller set of eigenmodes by analyzing its decomposition on a parametric basis. In Figure 14 we plot $1 - \tau_m$ from which the most relevant modes can be selected. For $\tau^* = 0.05$, we have $\mathcal{I}^* = \{39, 67, 75, 83, 107, 118\}$.

After acquiring the observed data from the simulated measurements on the target deformation, we add different levels of noise, $\gamma = [0.01, 0.1]$, as plotted in Figure 15.

*Discussion.* As in the previous case of the cube, a few specific modes are important to the reconstruction of the deformation although the first relevant mode is of relatively high frequency. This must be taken into account when computing $\mathcal{I}$ as we risk not including the relevant modes. The estimation shown in Figure 16 required more iterations to achieve apparent convergence, even with the lowest noise. As noise increases, the differentiation between modes becomes more difficult although the most important mode, 39, remains pronounced. For visualization purposes, the estimated deformation is reconstructed at different times of the estimation and depicted in Figure 17. As previously done, by adjusting $\varepsilon$ we can achieve a faster response of the estimate by reducing the number of iterations needed until apparent convergence, as shown in the estimation for this adjusted case, Figure 18.

Regarding performance, the model for the quasi-static problem has 7750 degrees of freedom (DoFs) and the computation of the first 120 eigenmodes takes 9 minutes. Finite elements of order one is used for the quasi-static problem as numerical locking is negligible in pure extensional deformation. The wave propagation problem has 420k DoFs and a total of 2223 time steps. One run of the wave propagation problem to generate the observations takes less than a minute and requires 180MB of RAM. During an estimation of 6 mode components, 7 wave propagation problems are run in parallel and each iteration of the estimation procedure takes 12 minutes. The computational cost of solving a wave propagation problem is higher during estimation due to parameter
change at each time step, requiring extra computations related to the update of the
tangent stiffness operator. We use a laptop computer equipped with a *Intel i9-9880H*
CPU and 32GB of RAM for computing the presented results.

Finally, the sensitivity analysis is done for a negligible noise $\gamma = 10^{-4}$ and $\varepsilon = 1$. The final estimated values for the sensitivity analysis are plotted in Figure 19a. In this case, the relevant modes are well distinguishable. The observability can be analyzed in the Gramian matrix, Figure 19c. Due to its sparsity, we conclude that there are only a few observable modes and the fact that the diagonal is not pronounced means that the modes are less indistinguishable from each other. This can be seen in the results with the reduced set of eigenmodes, where the estimation of less relevant modes are highly affected by noise. The high observability shown in Figure 19a is due to the fact that the deformation of the structure is of relatively high frequency, hence more observable according to the Gramian. Also, the relevant modes are not clustered around a certain frequency, which helps their spatial differentiation. The evolution of the parameters shown in Figure 19b can be analyzed and provide some insights into the process of estimation. Modes 67 and 87 present non-monotonic behavior since the deformation related to these modes is concentrated farther from the center where the waves were excited (see Figure 14), hence the wave takes longer to travel through these regions and its carried information to be considered.

Regarding performance of the sensitivity analysis, 115 wave propagation problems had to be run in parallel, taking 5 hours to complete and requiring 271GB of RAM, in total. The memory usage for each problem was significantly higher than when estimating only 6 modes as each problem had to store the deformation gradient of all 114 modes. Other strategies for accessing such gradients can be proposed, leaving room for improvement. A desktop workstation equipped with 2\*Intel Xeon Platinum 8276 CPU and 512GB of RAM was used to run these problems. Note that, although we use only one workstation, our implementation with $\oplus$MQ allows parallelization through workstations and clusters.
<figure><img src="image_10.png" /><figcaption>Figure 8: Evolution in time for the estimation $\hat{\theta}_h^k$ for the cube case. The estimated modes components are plotted with a highlighted region corresponding to the standard deviation. The target values are represented as dashed lines. Different $\gamma$ noise levels are tested with a regularization $\varepsilon = 1$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Snapshots at given iterations *k* and time *t* of the estimated deformation compared with the target deformation (transparent) for a high level of noise, *γ* = 0.2. From left to right, top to bottom we have (*k*=0, *t*=0µs) (*k*=0, *t*=25µs) (*k*=0, *t*=35µs) (*k*=1, *t*=25µs) (*k*=2, *t*=75µs) and (*k*=3, *t*=100µs).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Evolution in time for the estimation $\hat{\theta}_h^k$ for the cube case. The estimated modes components are plotted with a highlighted region corresponding to the standard deviation. The target values are represented as dashed lines. A weak regularization, $\varepsilon = 0.0025$, and noise $\gamma = 0.2$ is used.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>(a) Normalized estimated value for each parameter at the end of the sensitivity analysis.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Evolution of $\hat{\theta}_h$ for every parameter in $\mathcal{I}$. analysis.
(c) Normalized Gramian at the end of the</figcaption></figure>

Figure 11: The estimation is done for the larger set $\mathcal{I}$ to perform a sensitivity analysis in the cube illustration. Different outputs of the estimation procedure are plotted and the set $\mathcal{I}^*$ is highlighted.

<figure><img src="image_15.png" /></figure>

Figure 12: Quasi static problem configuration and target deformation for the plate illustration.
<figure><img src="image_16.png" /><figcaption>(a) Wave propagation problem configuration, mesh, excitation and measurement regions.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Example of measured signal.</figcaption></figure>

Figure 13: Configuration and examples of measured signals for the wave propagation problem. The measurement regions are represented in blue and the excited region in red.

<figure><img src="image_18.png" /><figcaption>(a) Reconstruction error per number of modes. eigenmodes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Deformation associated with the relevant</figcaption></figure>

Figure 14: Reconstruction of $u^*$ using the lowest frequency modes from which the selection of the most relevant modes is done.
<figure><img src="image_20.png" /><figcaption>Figure 15: Example of different noise levels for the aluminum plate illustration.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 16: Evolution in time for the estimation $\hat{\theta}_h^k$ for the plate case. The estimated modes components are plotted with a highlighted region corresponding to the standard deviation. The target values are represented as dashed lines. Different $\gamma$ noise levels are tested with a neutral regularization $\varepsilon = 1$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 17: Snapshots at given iterations *k* and time *t* of the estimated deformation compared with the target deformation (transparent). Estimation for $\gamma = 0.1$. From left to right, top to bottom we have ($k=0, t=60\mu s$) ($k=0, t=130\mu s$) ($k=0, t=180\mu s$) and ($k=1, t=200\mu s$). Visualization is scaled 300x.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 18: Evolution in time for the estimation $\hat{\theta}_h^k$ for the plate case. The estimated modes components are plotted with a highlighted region corresponding to the standard deviation. The target values are represented as dashed lines. A weaker regularization, $\varepsilon = 0.1$, is used.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>(a) Normalized estimated value for each parameter at the end of the sensitivity analysis.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b) Evolution of $\hat{\theta}_h$ for every parameter in $\mathcal{I}$. analysis.
(c) Normalized Gramian at the end of the</figcaption></figure>

Figure 19: The estimation is done for the larger set $\mathcal{I}$ to perform a sensitivity analysis in the plate illustration. Different outputs of the estimation procedure are plotted and the set $\mathcal{I}^*$ is highlighted..
4.1.3. Case of a 4-point bending on a steel pipe To increase complexity, we consider a configuration inspired from [43], where the authors use ultrasonic guided waves to detect weld damage in a steel pipe while it was subjected to a 4-point bending test. They assess the performance of traditional imaging techniques to detect the damage while the specimen is subjected to mechanical load. To overcome the bias introduced by the mechanical load to the ultrasonic signal, they perform signal acquisitions in pristine conditions for different loads and use them as baselines to, through comparison, detect changes in the signal due to potential damage. The baseline should be acquired for every expected load condition, which requires manipulations and reduces the robustness of a detection system. An alternative to acquiring several baselines is to estimate the deformation of the structure using the already available signals. For that purpose, we illustrate the use of our method in this more complex and application-related configuration. The modeled pipe has a length of 2.94m, a diameter of 0.1973m and is 8mm thick. The quasi-static problem is depicted in Figure 20a. For a total applied force of 220kN, the target deformation is depicted in Figure 20b.

<figure><img src="image_26.png" /><figcaption>Figure 20: Quasi static problem configuration and target deformation for the pipe
illustration.</figcaption></figure>

The excitation is done at the top outer-upper surface of the steel pipe at 30kHz, the observations are obtained from applying the observation operator, as previously described, with radial sensitivity in 48 regions distributed in 4 circular evenly distributed sections as depicted in Figure 21a. Each section of transducers is rotated 15 degrees with respect to each other. The synthetic observed data are obtained from simulating the wave propagation problem up to 1000µs. An example of the obtained signal is shown in Figure 21b.

The set $\mathcal{I}$ is embodied by the eigenmodes associated with the 60 lowest eigenvalues. From a guess of the deformation $\mathbf{u}_0^*$ generated by the bending surface forces, we can analyze its decomposition in the eigenbasis. In Figure 22 we plot $1 - \tau_m$ from which the most relevant modes can be selected and are therefore highlighted, making, for $\tau^* = 0.05$, $\mathcal{I}^* = \{0, 1, 2, 4, 6, 7, 8, 9\}$. Differently from the previous cases, the relevant modes are concentrated in the lowest frequencies and are poorly distributed.

Due to the stability issues associated with the potential non-coercivity of the wave propagation problem, the estimation has to be adjusted with $\epsilon = 10^4$, setting a relatively low update rate for the estimation, hence avoiding instabilities during estimation. The
<figure><img src="image_27.png" /><figcaption>(a) Wave propagation problem configuration, mesh, excitation and measurement regions.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) Example of measured signal.</figcaption></figure>

Figure 21: Configuration and examples of measured signals for the wave propagation problem. The measurement regions are represented in blue and the excited region in red.

<figure><img src="image_29.png" /><figcaption>Figure 22: Reconstruction error of $u^*$ using the lowest frequency modes from which the selection of the most relevant modes is done.</figcaption></figure>

evolution of the parameter can be seen in Figure 23 for two different levels of noise, $\gamma = 10^{-4}$ and $\gamma = 0.1$. For both levels of noise, the estimation succeeds in differentiating the modes. The estimated deformation are reconstructed at different steps of the estimation procedure and is depicted in Figure 24.

*Discussion.* Using high regularization and a large amount of observed data, the estimate is stable and shows good convergence to the target in a few iterations, even when noise increases. By comparison with the plate case, we believe that this is because sensors are well positioned and the acquisition time is long enough so that this case benefits from rather strong observability, allowing it to well distinguish the deformation modes. This illustrates the importance of a robust and well-designed acquisition configuration for a correct estimation.
<figure><img src="image_30.png" /><figcaption>Figure 23: Evolution in time for the estimation $\hat{\theta}_h^k$ for the pipe case. The estimated modes components are plotted with a highlighted region corresponding to the standard deviation. The target values are represented as dashed lines. Different $\gamma$ noise levels are tested with a strong regularization $\varepsilon = 10^4$.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 24: Snapshots at given iterations *k* and time *t* of the estimated deformation compared with the target deformation (transparent). Estimation for $\gamma = 0.01$. From left to right, top to bottom we have ($k$=0, $t$=60µs) ($k$=0, $t$=130µs) ($k$=0, $t$=180µs) and ($k$=1, $t$=200µs). Visualization is scaled 20x.</figcaption></figure>

Regarding performance of the estimation, The model for the quasi-static problem had 148k DoFs and the computation of the first 60 eigenmodes takes 8 hours to complete using a *Intel(R) Xeon(R) W-3245* CPU. Finite elements of order two is used for the quasi-static problem. The wave propagation problem has 2M DoFs and a total of 4662 time steps. One run of the wave propagation problem to generate the observations takes 5 minutes and requires 1GB of RAM. During an estimation of 8 mode components, 9 wave propagation problems are computed in parallel and each iteration of the estimation procedure takes 2 hours. A workstation equipped with a *Intel i9-9880H* CPU and 32GB of RAM is used for the estimation.

Note, however, that while we obtain satisfactory results when the estimation is performed in a reduced set of modes, the sensitivity analysis performs poorly when it comes to selecting the relevant modes from a larger set, as can be seen in Figure 25a. The high-frequency modes are overestimated because they are more observable, as seen in the Gramian (Figure 25c). This overestimation of the less relevant modes due to observability is a cause of instability when trying to estimate the deformation using the
larger set of modes. The evolution of the sensitivity analysis in Figure 25b shows that some irrelevant modes compensate for each other when trying to minimize the misfit because they are less distinguishable. Without the potential instability issue of the wave propagation problem, the Gramian diagonal indicates that the modes are readily observable. The Gramian diagonal can also be used to group the modes that conflict with each other in minimizing the mismatch, allowing selection of distinguishable modes and adaptive estimation. This is typical of ill-posed inverse problems, where estimation is efficient once we accept a regularization that here comes from a reconstruction in a parameter space with a small dimension.

Regarding performance of the sensitivity analysis, 61 wave propagation problems are computed in parallel, taking 5 hours to complete and using 400GB of RAM, in total. A desktop workstation equipped with 2× Intel Xeon Platinum 8276 and 512GB of RAM was used. The same remarks done for the plate case regarding performance are valid here.

## 5. Conclusion and perspectives

In this paper, we show how to reconstruct in nondestructive tests a baseline free of environmental loading conditions using available SHM measurements. Our model-based approach is inspired by Full-Wave Inversion (FWI) strategies, but here we take advantage of dynamic programming principles to avoid multiple iterations of adjoint minimization. In fact, we still rely on an iterative strategy, the Levenberg-Marquardt algorithm, to transform the initial nonlinear inverse problem into successive linear-quadratic estimation problems. Each of these linear-quadratic estimation problems is then solved using a Kalman-based approach. By relying on the Unscented-Kalman filter, we also avoid the computation of the tangent operator required by the Levenberg-Marquardt algorithm. The final algorithm converges in a limited number of iterations compared to standard FWI and is therefore mostly sequential and inherently parallel with respect to the sigma-points covering the reconstruction space. As a result, our algorithm efficiency remains comparable to solving the direct guided wave propagation problem. The main limitation of our approach is that, due to followed dynamic programming point of view, we are limited to a low-dimensional representation of the estimated deformation, i.e., only about a hundred modes can be estimated. Moreover, and this will be the goal of future developments, our approach first considers the reconstruction of the displacement field, while from the identification perspective, it might be better suited to reconstruct a strain tensor, for instance, the Cauchy-Green strain or its invariants. In particular, we believe that this research direction will help us to better control some deformation constraints related to a potentially ill-posed wave propagation problem that limits the trust region neighborhood of our current method. Moreover, our method now needs to be extended to more general observational operators, some of which may depend on the estimated deformation, making the presented approach more complex. Moreover, our approach suffers from the
<figure><img src="image_32.png" /><figcaption>(a) Normalized estimated value for each parameter at the end of the sensitivity analysis.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(b) Evolution of $\hat{\theta}_h$ for every parameter in $\mathcal{I}$. analysis.
(c) Normalized Gramian at the end of the</figcaption></figure>

Figure 25: The estimation is done for the larger set $\mathcal{I}$ to perform a sensitivity analysis in the pipe illustration. Different outputs of the estimation procedure are plotted and the set $\mathcal{I}^*$ is highlighted.

same limitation as the classical FWI method in terms of the choice of the discrepancy measure between the model and the data. Recently, [17, 18] has been shown that the use of Wasserstein-based distances can overcome cycle-skipping effects in wave field inversions.

## Acknowledgment and fundings

The authors would like to deeply thank Laurent Steff for his implementation of the AKILLES Library and Jérôme Diaz for his help with the MoReFEM Library. The authors would also like to thank the Isaac Newton Institute for Mathematical Sciences for support (EPSRC grant number EP/R014604/1) and hospitality during the programme "The mathematical and statistical foundation of future data-driven engineering" when
part of this work on this paper was undertaken. This research was funded by the
following project: “GW4SHM”(gw4shm.eu) project from the European Union’s Horizon
2020 Research and Innovation program under the Marie Skłodowska-Curie, grant
number 860104.

References

[1] Z. Abiza, M. Destrade, and R. W. Ogden. Large acoustoelastic effect. *Wave Motion*, 49(2):364–374, 2012.

[2] S Afshar and K. Germ and F. Morris. Well-posedness of extended Kalman filter equations for semilinear infinite-dimensional systems. In *59th IEEE Conference on Decision and Control*, 1210–1215, 2020.

[3] M. Aussal and P. Moireau. Kernel representation of kalman observer and associated h-matrix based discretization. *ESAIM: Control, Optimisation and Calculus of Variations*, 28:78, 2022.

[4] G. Bal, W. Naetar, Scherzer. O., and J. Schotland. The levenberg–marquardt iteration for numerical inversion of the power density operator. *Journal of Inverse and Ill-Posed Problems*, 21(2):265–280, 2013.

[5] A. Bensoussan. *Filtrage optimal des systèmes linéaires*, 3. Dunod, 1971.

[6] A. Bensoussan, G. Da Prato, M. C. Delfour, and S. K. Mitter. *Representation and Control of Infinite Dimensional Systems*. Systems & Control: Foundations & Applications. Birkhäuser Boston, Boston, MA, 2007.

[7] S. Bernard, V. Monteiller, D. Komatitsch, and P. Lasaygues. Ultrasonic computed tomography based on full-waveform inversion for bone quantitative imaging. *Physics in Medicine & Biology*, 62(17):7011, August 2017.

[8] J. Blum, F.-X. Le Dimet, and I. M. Navon. Data assimilation for geophysical fluids. In *Handbook of numerical analysis*, 14:385–441. Elsevier, 2009.

[9] P. G. Ciarlet. *Mathematical elasticity. Volume I, Three-dimensional elasticity*. North-Holland ; Sole distributors for the U.S.A. and Canada, Elsevier Science Pub. Co., Amsterdam, New York, 1988.

[10] G. C. Cohen. *Higher-order numerical methods for transient wave equations*. Scientific computation. Springer-Verlag, Berlin Heidelberg, 1st edition, 2002.

[11] P. Courtier, J.-N. Thépaut, and A. Hollingsworth. A strategy for operational implementation of 4D-Var, using an incremental approach. *Quarterly Journal of the Royal Meteorological Society*, 120(519):1367–1387, July 1994.

[12] A. Dalmora, A. Imperiale, S. Imperiale, and P. Moireau. A time-domain spectral finite element method for acoustoelasticity: modeling the effect of mechanical loading on guided wave propagation. *Wave Motion*, 129, 103328.

[13] M. De Buhan and M. Kray. A new approach to solve the inverse scattering problem for waves: combining the TRAC and the adaptive inversion methods. *Inverse Problems*, 29(8):085009, 2013.

[14] F.-X. Le Dimet and O. Talagrand. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. *Tellus A*, 38A(2):97–110, 1986.

[15] M. Duruflé, P. Grob, and P. Joly. Influence of Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain. *Numerical Methods for Partial Differential Equations*, 25(3):526–551, 2009.

[16] G. Duvaut and J. L. Lions. *Inequalities in mechanics and physics*. Springer-Verlag, 1976.

[17] B. Engquist, B. D. Froese, and Y. Yang. Optimal transport for seismic full waveform inversion. *Communications in Mathematical Sciences*, 14(8):2309–2330, 2016.

[18] B. Engquist, K. Ren, and Y. Yang. The quadratic Wasserstein metric for inverse data matching. *Inverse Problems*, 36(5):055001, 2020.
[19] T. Furuya and R. Potthast. Inverse medium scattering problems with Kalman filter techniques. *Inverse Problems*, 38(9):095003, 2022.

[20] N. Gandhi, J. E. Michaels, and S. J. Lee. Acoustoelastic Lamb wave propagation in biaxially stressed plates. *The Journal of the Acoustical Society of America*, 132(3):1284–1293, 2012.

[21] R. Gorgin, Y. Luo, and Z. Wu. Environmental and operational conditions effects on Lamb wave based structural health monitoring systems: A review. *Ultrasonics*, 106114, 2020.

[22] W. Haik, Y Maday, and L. Chamoin. A real-time variational data assimilation method with data-driven model enrichment for time-dependent problems. *Computer Methods in Applied Mechanics and Engineering*, 405:115868, 2023.

[23] M. Hanke. A regularizing levenberg - marquardt scheme, with applications to inverse groundwater filtration problems. *Inverse Problems*, 13(1):79, 1997.

[24] A. Imperiale, D. Chapelle, and P. Moireau. Sequential data assimilation for mechanical systems with complex image data: application to tagged-MRI in cardiac mechanics. *Advanced Modeling and Simulation in Engineering Sciences*, 8(1):2, 2021.

[25] Julier, Simon J. and Uhlmann, J. K. and Durrant-Whyte, H. F. A new method for the nonlinear transformation of means and covariances in filters and estimators *IEEE Transactions on Automatic Control* 45(3):477–482, 2000.

[26] D. G. Luenberger and Y. Ye. *Linear and nonlinear programming*, volume 2 of *International series in operations research and management science*. Springer, New York, NY, 3rd edition, 2008.

[27] B. Marchand, L. Chamoin, and C. Rey. Real-time updating of structural mechanics models using Kalman filtering, modified constitutive relation error, and proper generalized decomposition. *International Journal for Numerical Methods in Engineering*, 107(9):786–810, 2016.

[28] M. Mitra and S. Gopalakrishnan. Guided wave based structural health monitoring: A review. *Smart Materials and Structures*, 25(5):053001, May 2016.

[29] P. Moireau. A Discrete-time Optimal Filtering Approach for Non-linear Systems as a Stable Discretization of the Mortensen Observer. *ESAIM: Control, Optimisation and Calculus of Variations*, 24(4):1815–1847, 2019.

[30] P. Moireau. Discrete-time formulations as time discretization strategies in data assimilation. In *Handbook of Numerical Analysis*, Handbook of Numerical Analysis. Elsevier, Chap. 9, Volume 24, 297–339, 2023.

[31] P. Moireau and D. Chapelle. Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems. *ESAIM: Control, Optimisation and Calculus of Variations*, 17(2):380–405, 2011.

[32] P. Moireau, D. Chapelle, and P. Le Tallec. Joint state and parameter estimation for distributed mechanical systems. *Computer Methods in Applied Mechanics and Engineering*, 197(6-8):659–677, 2008.

[33] P. Moireau and D. Chapelle. Erratum of article "reduced-order Unscented Kalman filtering with application to parameter identification in large-dimensional systems". *ESAIM: COCV*, 17(2):406–409, 2011.

[34] B. Moore. Principal component analysis in linear systems: Controllability, observability, and model reduction. *IEEE Transactions on Automatic Control*, 26(1):17–32, 1981.

[35] A. Pazy. *Semigroups of Linear Operators and Applications to Partial Differential Equations*, volume 44 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1983.

[36] D. T. Pham, J. Verron, and L. Gourdeau. Filtres de Kalman singuliers évolutifs pour l'assimilation de données en océanographie. *Comptes Rendus de l'Académie des Sciences - Series IIA - Earth and Planetary Science*, 326(4):255–260, 1998.

[37] K. Ramdani, M. Tucsnak, and G. Weiss. Recovering the initial state of an infinite-dimensional system using observers. *Automatica*, 46(10):1616–1625, 2010.

[38] F. Ricci, E. Monaco, N. D. Boffa, L. Maio, and V. Memmolo. Guided waves for structural health monitoring in composites: A review and implementation strategies. *Progress in Aerospace Sciences*, 129:100790, 2022.
[39] D. Rozier, F. Birol, E. Cosme, P. Brasseur, J. M. Brankart, and J. Verron. A Reduced-Order Kalman Filter for Data Assimilation in Physical Oceanography. *SIAM Review*, 49(3):449–465, 2007.

[40] M. Shams, M. Destrade, and R. W. Ogden. Initial stresses in elastic solids: Constitutive laws and acoustoelasticity. *Wave Motion*, 48(7):552–567, 2011.

[41] D. Simon. *Optimal state estimation: Kalman, H [infinity] and nonlinear approaches*. Wiley-Interscience, Hoboken, N.J, 2006.

[42] S. Särkkä. *Bayesian Filtering and Smoothing*. Cambridge University Press, 1st edition, 2013.

[43] K. Tschöke, B. Weihnacht, E. Schulze, T. Gaul, L. Schubert, and R. Neubeck. Determination of Defect Sizes with the help of Structural-Health-Monitoring Methods based on Guided Waves. In *7th European-American Workshop on Reliability of NDE*, 2017.

[44] J. Virieux, A. Asnaashari, R. Brossier, L. Métivier, Alessandra Ribodetti, and Wei Zhou. An introduction to full waveform inversion. In *Encyclopedia of exploration geophysics*, pages R1–1. Society of Exploration Geophysicists, 2017.