A Gauss–Newton full-waveform inversion for material profile reconstruction in viscoelastic semi-infinite solid media

Abstract

An inversion framework employing a Gauss–Newton method is developed to reconstruct material profiles in heterogeneous, viscoelastic, semi-infinite domains. In particular, a full-waveform inversion approach is investigated to image the elastic and attenuating parameters of a layered media. To account for the viscoelasticity of the medium, a Generalized Maxwell Body with one spring and two Maxwell elements in parallel (GMB2) is adopted in the forward and inverse wave propagation problems. Perfectly-matched-layers were introduced as wave absorbing buffers to simulate the semi-infinite extent of the domain. Using transient wave equations endowed with the GMB2 constitutive relation and the PML, a partial-differential-equations-constrained optimization scheme was implemented that lead to classic KKT (Karush–Kuhn–Tucker) conditions including time-dependent state, adjoint, and time-invariant control problems. An optimal solution of the viscoelastic parameters was obtained using a reduced-space approach based on a line search algorithm where the search direction was computed by the Gauss–Newton method. Considerable improvements on the accuracy and convergence rate of solutions were made by the developed Gauss–Newton inversion procedure compared to previous research using the Fletcher–Reeves method.

Keywords:

• full-waveform inversion
• Gauss–Newton method
• semi-infinite domain
• viscoelasticity
• Generalized Maxwell Body
• perfectly-matched-layers (PMLs)
• Lagrangian
• PDE-constrained optimization
• KKT system
• reduced-space approach

AMS Subject Classifications:

• 65M32
• 65M60
• 49M25
• 74D99

1. Introduction

Recently, perfectly-matched-layers (PMLs) have been widely used as wave-absorbing boundaries for elastodynamic wave propagation problems in unbounded domains. Compared to the remarkable progress made for the elastodynamic PML, there has been substantially less development of the PML in viscoelastic media. Among the developments for viscoelastic media is the work by Bécache et al. [1] who adapted the PML to the Zener viscoelastic model and extended a mixed finite element method proposed in [2] for numerical solution. Ma and Liu associated the PML with the intrinsic attenuation parameter $Q$ and implemented it with an explicit finite element method.[3] They used the coarse-grain approach for modeling $Q$ in a structured finite-element mesh.[4] Later, Ely et al. [5] extended the work of [3] by using a generalized finite difference method for viscoelastic wave modeling in 3D heterogeneous media.

An advancing, unsplit-field, convolutional PML was also proposed for 3D viscoelastic seismic wave modeling in the context of a fourth-order staggered finite-difference method.[6] This formulation was particularly useful when two or more standard linear solid mechanisms were used. Dhemaied et al. [7] introduced modeling of seismic wave propagation in viscoelastic media using the auxiliary differential equation method. The formulation was based on a convolutional PML in a homogeneous viscoelastic medium. Based on Carcione’s theory,[8] Yan and Liu [9] presented a two-dimensional first order velocity-stress formulation for viscoelastic, tilted, and transversely isotropic media. They introduced the PML in a viscoelastic media, and used a finite difference scheme with a rotated staggered grid to solve the equations. Most recently, Assimaki et al. [10] demonstrated time-domain forward and inverse modeling of lossy soils in viscoelastic heterogeneous PML-truncated domains. They employed a Generalized Maxwell Body with one spring and two Maxwell elements in parallel (GMB2) as the constitutive model of a viscoelastic medium. Their approach yielded the reconstruction of the GMB2 material properties in a PML-truncated domain via a full-waveform inversion method.

This paper extends the inversion framework of [10] and discusses improvements on the accuracy and convergence rate of solutions in the full-waveform inversion. In particular, it is sought to recover the spatial distribution of the GMB2 material parameters $E_1$ and $E_2$ from surficial measurements of the domain’s response to interrogating waves via a Gauss–Newton full-waveform inversion method in the time-domain. The reconstructed profiles of $E_1$ and $E_2$ can be utilized to construct the profile of the quality factor, $Q$, a commonly used metric of intrinsic and apparent attenuation.

To reconstruct viscoelastic profiles, a partial-differential-equations (PDE) constrained optimization approach is introduced, where an objective functional representing the misfit between measured and computed surficial displacement responses is minimized subject to the satisfaction of a set of viscoelastic wave equations. The PML-endowed wave equations developed in [10] based on the GMB2 constitutive model are used as the constraint in the optimization framework. Specifically, a Lagrangian functional is constructed by augmenting the misfit through the side imposition of the governing PDEs via Lagrange multipliers. Enforcing the stationarity of the Lagrangian leads to KKT (Karush–Kuhn–Tucker) conditions comprising time-dependent state, adjoint, and time-independent control problems. Here, the state and adjoint problems are initial- and final-value problems, respectively, and both are resolved using a mixed finite element method. The optimal solution of material parameters is obtained via a reduced-space approach, where the updates of $E_1$ and $E_2$ are determined by a Gauss–Newton–Krylov optimization method at each inversion iteration.

In the mixed finite element method, special pairs of basic functions were used to avoid numerical instabilities that may occur when dealing with forward state as well as backward adjoint problems.[11, 12] In the implementation of the developed inversion procedure, storage problems due to reverse time backward fields were not prohibitive because of the dimensional simplification to 1D heterogeneous media.

The significance of this research lies on the fact that the use of PML in conjunction with the Gauss–Newton–Krylov method can be one of the best solution strategies possible for the inversion of viscoelastic properties in semi-infinite domains in terms of the accuracy and convergence of solutions. Through a series of numerical examples, it is reported that substantial improvements on the accuracy and convergence rate of solutions were attained by the developed inversion procedure, compared to the results of [10] where the Fletcher–Reeves method was used.

2. Forward problem

The forward problem is posed in a horizontally layered, vertically heterogeneous viscoelastic solid medium of semi-infinite extent as depicted in Figure 1(a). The forward problem of interest herein is to determine compressional wave motion in such a layered medium due to a stress load $p(t)$ applied on the surface. The problem can be reduced to a one-dimensional problem in the direction of the wave propagation. To simulate the wave motion in the semi-infinite extent, the semi-infinite domain is truncated and a PML wave-absorbing boundary is introduced beyond the truncation surface, as shown in Figure 1(b). Figure 1(c) shows the schematic of a PML-truncated finite computational domain comprising the regular domain ($0\le x\le L$) and the PML ($L\le x\le L_t$).

Figure 1. (a) Original domain: horizontally-layered heterogeneous semi-infinite solid medium; (b) PML-truncated domain; (c) One-dimensional schematic of the PML-truncated domain (PML location: $L\le x\le L_t$).

Figure 2. Schematic representation of a Generalized Maxwell Body.

In recent work, Assimaki et al. [10] discussed the constitutive model of a reduced Generalized Maxwell Body (GMB2) that can effectively capture the frequency-independence of $Q$. The model consists of three elastic springs ($E_1,E_2,E_3$) and two dashpots (${\mathit{\eta }}_2,{\mathit{\eta }}_3$) in such a way that one spring and two Maxwell elements are connected in parallel, as shown in Figure 2. For the two Maxwell elements, relaxation time parameters ${\mathit{\tau }}_1$ and ${\mathit{\tau }}_2$ were defined as:(1) $$\begin{array}{c}\hfill {\mathit{\tau }}_1=\frac{{\mathit{\eta }}_2}{E_2},{\mathit{\tau }}_2=\frac{{\mathit{\eta }}_3}{E_3},\end{array}$$(1)

where ${\mathit{\tau }}_1$ and ${\mathit{\tau }}_2$ can be shown to be constants within a certain frequency range.[4, 10] For the three moduli of the GMB2, it was shown that:(2) $$\begin{array}{cc}\hfill E_2& =E_3,\hfill \end{array}$$(2) (3) $$\begin{array}{cc}\hfill \frac{E_2}{E_1+2E_2}& =\mathit{\alpha }{Q}^{\mathit{\beta }},\hfill \end{array}$$(3)

where $\mathit{\alpha }=1.768$ and $\mathit{\beta }=-0.98$ are constants. Equations (2(2) $$\begin{array}{cc}\hfill E_2& =E_3,\hfill \end{array}$$(2) ) and (3(3) $$\begin{array}{cc}\hfill \frac{E_2}{E_1+2E_2}& =\mathit{\alpha }{Q}^{\mathit{\beta }},\hfill \end{array}$$(3) ) show that the Generalized Maxwell Body depicted in Figure 2 can be described by any of the pairs ($E_1,E_2$), ($E_1$, $Q$), or ($E_2$, $Q$). If $E_1$ and $E_2$ are taken as two independent material parameters, the damping coefficients ${\mathit{\eta }}_2$ and ${\mathit{\eta }}_3$ can be obtained from (1(1) $$\begin{array}{c}\hfill {\mathit{\tau }}_1=\frac{{\mathit{\eta }}_2}{E_2},{\mathit{\tau }}_2=\frac{{\mathit{\eta }}_3}{E_3},\end{array}$$(1) ) as ${\mathit{\eta }}_2={\mathit{\tau }}_1{E}_2$ and ${\mathit{\eta }}_3={\mathit{\tau }}_2{E}_2$, respectively. Based on the GMB2 constitutive model, Assimaki et al. showed that viscoelastic wave motion in PML-truncated solid domains can be described by a system of coupled PDEs in the time-domain.[10] The governing PDEs are adopted herein to establish the forward problem in one-dimensional viscoelastic solid domains truncated by the PML. The problem can be stated as: find $u=u(x,t)$ and $\mathit{\sigma }=\mathit{\sigma }(x,t)$, such that:(4a) $$\begin{array}{cc}& \frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}=0,\hfill \end{array}$$(4a) (4b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}=0,\hfill \end{array}$$(4b) (4c) $$\begin{array}{cc}& \mathrm{f}or0<x<L_t,0<t\le T,\hfill \end{array}$$(4c) and(5a) $$\begin{array}{cc}& u(L_t,t)=0,\hfill \end{array}$$(5a) (5b) $$\begin{array}{cc}& \mathit{\sigma }(0,t)=p(t),\hfill \end{array}$$(5b) (5c) $$\begin{array}{cc}& u(x,0)=\frac{\mathit{\partial }u}{\mathit{\partial }t}(x,0)=0,\hfill \end{array}$$(5c) (5d) $$\begin{array}{cc}& \tilde{\mathit{\sigma }}(x,0)=\mathit{\sigma }(x,0)=\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}(x,0)=0,\hfill \end{array}$$(5d) where $x$ denotes location and $t$ denotes time. In the above, $u(x,t)$, $\mathit{\sigma }(x,t)$, and $\mathit{\rho }$ denote displacement, stress, and mass density, respectively, and $\tilde{\mathit{\sigma }}\equiv {\int }_0^t\mathit{\sigma }(x,z)dz$ is an auxiliary stress-memory term. Additionally, $h$ is defined as $h(x)\equiv c_{\mathrm{ref}}g(x)$, where $c_{\mathrm{ref}}$ is a reference wave velocity of the medium and $g(x)$ is a function that plays a central role in artificial wave attenuation within the PML. The attenuation function $g(x)$ was defined in [11] as:(6) $$\begin{array}{c}\hfill g(x)=\left\{\begin{array}{cc}0,\hfill & 0\le x\le L,\hfill \\ {\displaystyle \frac{3}{2L_{\mathrm{P}ML}}\mathrm{l}og\left(\frac{1}{|R|}\right){\left(\frac{x-L}{L_{\mathrm{P}ML}}\right)}^2,}\hfill & L\le x\le L_t,\hfill \end{array}\right.\end{array}$$(6)

in which $L_{\mathrm{P}ML}=L_t-L$ is the length of the PML and $|R|$ is a user-tunable reflection coefficient that controls the amount of reflection from the fixed PML end $(x=L_t)$. In Equation (4b(4b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}=0,\hfill \end{array}$$(4b) ), the coefficients ${\mathit{\kappa }}_i(i=1,\dots ,4)$ are auxiliary variables in terms of $h$, relaxation times, and the two material parameters $E_1$ and $E_2$ of the GMB2 constitutive model. They are expressed as:(7) $$\begin{array}{cc}\hfill {\mathit{\kappa }}_1& ={\mathit{\tau }}_1+{\mathit{\tau }}_2+h{\mathit{\tau }}_1{\mathit{\tau }}_2,\hfill \end{array}$$(7) (8) $$\begin{array}{cc}\hfill {\mathit{\kappa }}_2& =1+h\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right),\hfill \end{array}$$(8) (9) $$\begin{array}{cc}\hfill {\mathit{\kappa }}_3& =\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\left(E_1+E_2\right),\hfill \end{array}$$(9) (10) $$\begin{array}{cc}\hfill {\mathit{\kappa }}_4& ={\mathit{\tau }}_1{\mathit{\tau }}_2\left(E_1+2E_2\right).\hfill \end{array}$$(10)

Equations (4a(4a) $$\begin{array}{cc}& \frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}=0,\hfill \end{array}$$(4a) ) and (4b(4b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}=0,\hfill \end{array}$$(4b) ) are the displacement ($u$)-stress ($\mathit{\sigma }$) mixed equations governing the propagation of waves in one-dimensional viscoelastic PML-truncated domains. Equations (5a(5a) $$\begin{array}{cc}& u(L_t,t)=0,\hfill \end{array}$$(5a) ) and (5b(5b) $$\begin{array}{cc}& \mathit{\sigma }(0,t)=p(t),\hfill \end{array}$$(5b) ) represent displacement and traction boundary conditions, respectively, while Equations (5c(5c) $$\begin{array}{cc}& u(x,0)=\frac{\mathit{\partial }u}{\mathit{\partial }t}(x,0)=0,\hfill \end{array}$$(5c) ) and (5d(5d) $$\begin{array}{cc}& \tilde{\mathit{\sigma }}(x,0)=\mathit{\sigma }(x,0)=\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}(x,0)=0,\hfill \end{array}$$(5d) ) indicate silent initial conditions for the system. This initial-boundary-value-problem (IBVP) (Equations (4) and (5)) can be solved using a mixed finite element method.

3. Inverse problem in viscoelastic PML-truncated domains

3.1. PDE-constrained optimization problem

The problem of interest in this inverse problem is to reconstruct the spatial variation of viscoelastic parameters $E_1$ and $E_2$ of the GMB2 in PML-truncated domains given surficial measurements of the domain’s response to a known excitation. To accomplish this reconstruction, the following objective functional is to be minimized:(11) $$\begin{array}{c}\hfill \mathcal{J}=\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2),\end{array}$$(11)

subject to Equations (4)–(5).

The above problem is an optimization problem constrained by the PML-endowed viscoelastic wave equations. In Equation (11(11) $$\begin{array}{c}\hfill \mathcal{J}=\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2),\end{array}$$(11) ), the objective functional $\mathcal{J}$ comprises a misfit functional representing the difference between computed ($u(0,t)$) and measured ($u_{\mathrm{m}}(0,t)$) surficial responses in a least-squares sense, augmented by regularization terms ${\mathcal{R}}_{E_1}(E_1)$ and ${\mathcal{R}}_{E_2}(E_2)$. The computed response refers to the one obtained by solving the IBVP (Equations (4) and (5)) under assumed trial profiles of $E_1(x)$ and $E_2(x)$. To relieve the inherent ill-posedness of the inverse problem, $\mathcal{J}$ includes regularization terms with respect to $E_1$ and $E_2$, respectively. When Tikhonov (TN) regularization is employed,[13] ${\mathcal{R}}_{E_1}(E_1)$ and ${\mathcal{R}}_{E_2}(E_2)$ become:(12a) $$\begin{array}{cc}\hfill {\mathcal{R}}_{E_1}(E_1)& =\frac{R_{E_1}}{2}{\int }_0^{L_t}{\left(\frac{\mathrm{d}{E}_1}{\mathrm{d}x}\right)}^2\mathrm{d}x,\hfill \end{array}$$(12a) (12b) $$\begin{array}{cc}\hfill {\mathcal{R}}_{E_2}(E_2)& =\frac{R_{E_2}}{2}{\int }_0^{L_t}{\left(\frac{\mathrm{d}{E}_2}{\mathrm{d}x}\right)}^2\mathrm{d}x.\hfill \end{array}$$(12b) When Total Variation (TV) regularization is used,[14] the corresponding terms become:(13a) $$\begin{array}{cc}\hfill {\mathcal{R}}_{E_1}(E_1)& =R_{E_1}{\int }_0^{L_t}{\left[{\left(\frac{\mathrm{d}{E}_1}{\mathrm{d}x}\right)}^2+e\right]}^{\frac{1}{2}}\mathrm{d}x,\hfill \end{array}$$(13a) (13b) $$\begin{array}{cc}\hfill {\mathcal{R}}_{E_2}(E_2)& =R_{E_2}{\int }_0^{L_t}{\left[{\left(\frac{\mathrm{d}{E}_2}{\mathrm{d}x}\right)}^2+e\right]}^{\frac{1}{2}}\mathrm{d}x,\hfill \end{array}$$(13b) where $R_{E_1}$ and $R_{E_2}$ are regularization factors for the target inversion variables $E_1$ and $E_2$. In Equation (13a(13a) $$\begin{array}{cc}\hfill {\mathcal{R}}_{E_1}(E_1)& =R_{E_1}{\int }_0^{L_t}{\left[{\left(\frac{\mathrm{d}{E}_1}{\mathrm{d}x}\right)}^2+e\right]}^{\frac{1}{2}}\mathrm{d}x,\hfill \end{array}$$(13a) ), a small auxiliary parameter $e$ was added to make it possible to differentiate the TV regularization. The objective functional $\mathcal{J}$ is sought to be minimized, in an attempt to reconstruct the viscoelastic material properties $E_1(x)$ and $E_2(x)$.

3.2. First-order optimality conditions

To solve the PDE-constrained optimization problem defined in Equation (11(11) $$\begin{array}{c}\hfill \mathcal{J}=\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2),\end{array}$$(11) ) for $E_1$ and $E_2$, a Lagrangian functional is introduced, where the objective functional $\mathcal{J}$ is augmented by the weak imposition of the governing PDEs and a boundary condition:(14) $$\begin{array}{cc}& \mathcal{L}(u,\mathit{\sigma },{\mathit{\lambda }}_u,{\mathit{\lambda }}_{\mathit{\sigma }},E_1,E_2),\hfill \\ \hfill & =\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2)\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_u\left[\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}\right]\mathrm{d}t\mathrm{d}x\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_{\mathit{\sigma }}\left[{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}\right]\mathrm{d}t\mathrm{d}x.\hfill \\ \hfill \end{array}$$(14)

In Equation (14(14) $$\begin{array}{cc}& \mathcal{L}(u,\mathit{\sigma },{\mathit{\lambda }}_u,{\mathit{\lambda }}_{\mathit{\sigma }},E_1,E_2),\hfill \\ \hfill & =\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2)\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_u\left[\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}\right]\mathrm{d}t\mathrm{d}x\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_{\mathit{\sigma }}\left[{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}\right]\mathrm{d}t\mathrm{d}x.\hfill \\ \hfill \end{array}$$(14) ), ${\mathit{\lambda }}_u$ and ${\mathit{\lambda }}_{\mathit{\sigma }}$ are Lagrange multipliers used to side-impose the governing PDEs (4a(4a) $$\begin{array}{cc}& \frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}=0,\hfill \end{array}$$(4a) ) and (4b(4b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}=0,\hfill \end{array}$$(4b) ). The boundary and initial conditions (Equation (5)) can be explicitly imposed. It is then sought to satisfy the stationarity of $\mathcal{L}$, by requiring that the first variations of $\mathcal{L}$ with respect to state, adjoint, and control variables vanish. This results in the following optimality conditions:

3.2.1. The first optimality condition

First, the first variations of $\mathcal{L}$ with respect to the Lagrange multipliers ${\mathit{\lambda }}_u$ and ${\mathit{\lambda }}_{\mathit{\sigma }}$ must vanish, i.e.:(15a) $$\begin{array}{cc}\hfill {\mathit{\delta }}_{{\mathit{\lambda }}_u}\mathcal{L}& =0,\hfill \end{array}$$(15a) (15b) $$\begin{array}{cc}\hfill {\mathit{\delta }}_{{\mathit{\lambda }}_{\mathit{\sigma }}}\mathcal{L}& =0.\hfill \end{array}$$(15b) Equation (15a(15a) $$\begin{array}{cc}\hfill {\mathit{\delta }}_{{\mathit{\lambda }}_u}\mathcal{L}& =0,\hfill \end{array}$$(15a) ) results in the mixed state (or forward) problem, identical to the IBVP given by Equations (4) and (5). Using operator notation, the state problem can be written as:(16) $$\begin{array}{cc}\hfill \mathcal{A}(\mathit{E})\mathit{u}& =\mathbf{0},\mathrm{f}or0<x<L_t,0<t\le T,\hfill \end{array}$$(16) (17) $$\begin{array}{cc}\hfill {\mathcal{B}}_s\mathit{u}& =\mathbf{0},\hfill \end{array}$$(17) (18) $$\begin{array}{cc}\hfill {\mathcal{I}}_s\mathit{u}& =\mathbf{0},\hfill \end{array}$$(18)

where:(19) $$\begin{array}{cc}& \mathcal{A}(\mathit{E})\mathit{u}=\left[\begin{array}{c}{\displaystyle \frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}}\hfill \\ {\displaystyle {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \end{array}\right],\hfill \\ \hfill & {\mathcal{B}}_s\mathit{u}=\left[\begin{array}{c}u(L_t,t)\hfill \\ \mathit{\sigma }(0,t)-p(t)\hfill \end{array}\right],{\mathcal{I}}_s\mathit{u}=\left[\begin{array}{c}u(x,0)\\ {\displaystyle \frac{\mathit{\partial }u}{\mathit{\partial }t}(x,0)}\\ \tilde{\mathit{\sigma }}(x,0)\\ \mathit{\sigma }(x,0)\\ {\displaystyle \frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}(x,0)}\end{array}\right],\hfill \\ \hfill & \mathit{u}=\left[\begin{array}{c}u\\ \mathit{\sigma }\end{array}\right],\mathit{E}=\left[\begin{array}{c}{E}_1\\ E_2\end{array}\right].\hfill \end{array}$$(19)

In Equation (19(19) $$\begin{array}{cc}& \mathcal{A}(\mathit{E})\mathit{u}=\left[\begin{array}{c}{\displaystyle \frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}}\hfill \\ {\displaystyle {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \end{array}\right],\hfill \\ \hfill & {\mathcal{B}}_s\mathit{u}=\left[\begin{array}{c}u(L_t,t)\hfill \\ \mathit{\sigma }(0,t)-p(t)\hfill \end{array}\right],{\mathcal{I}}_s\mathit{u}=\left[\begin{array}{c}u(x,0)\\ {\displaystyle \frac{\mathit{\partial }u}{\mathit{\partial }t}(x,0)}\\ \tilde{\mathit{\sigma }}(x,0)\\ \mathit{\sigma }(x,0)\\ {\displaystyle \frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}(x,0)}\end{array}\right],\hfill \\ \hfill & \mathit{u}=\left[\begin{array}{c}u\\ \mathit{\sigma }\end{array}\right],\mathit{E}=\left[\begin{array}{c}{E}_1\\ E_2\end{array}\right].\hfill \end{array}$$(19) ), $\mathcal{A}(\mathit{E})$ is defined as a state operator on the vector of state variables $\mathit{u}={\left[u\mathit{\sigma }\right]}^{\mathrm{T}}$, while ${\mathcal{B}}_s$ and ${\mathcal{I}}_s$ are defined as the operators indicating boundary and initial conditions, respectively.

3.2.2. The second optimality condition

Similarly, the vanishing of the first variations of $\mathcal{L}$ with respect to the state variables $u$ and $\mathit{\sigma }$ is enforced, i.e.:(20a) $$\begin{array}{cc}\hfill {\mathit{\delta }}_u\mathcal{L}& =0,\hfill \end{array}$$(20a) (20b) $$\begin{array}{cc}\hfill {\mathit{\delta }}_{\mathit{\sigma }}\mathcal{L}& =0,\hfill \end{array}$$(20b) which yields the following mixed adjoint problem: Find ${\mathit{\lambda }}_u={\mathit{\lambda }}_u(x,t)$ and ${\mathit{\lambda }}_{\mathit{\sigma }}={\mathit{\lambda }}_{\mathit{\sigma }}(x,t)$, such that:(21a) $$\begin{array}{cc}& \mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathrm{\Theta }}{\mathit{\partial }x}=0,\hfill \end{array}$$(21a) (21b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}-{\mathit{\kappa }}_1\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2{\mathit{\lambda }}_{\mathit{\sigma }}-h{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}+\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}=0,\hfill \end{array}$$(21b) (21c) $$\begin{array}{cc}& \mathrm{f}or0<x<L_t,0\le t<T,\hfill \end{array}$$(21c) and(22a) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_u(L_t,t)& =0,\hfill \end{array}$$(22a) (22b) $$\begin{array}{cc}\hfill \mathrm{\Theta }(0,t)& =u(0,t)-u_{\mathrm{m}}(0,t),\hfill \end{array}$$(22b) (22c) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_u(x,T)& =\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}(x,T)=0,\hfill \end{array}$$(22c) (22d) $$\begin{array}{cc}\hfill {\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}(x,T)& ={\mathit{\lambda }}_{\mathit{\sigma }}(x,T)=\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}(x,T)=0,\hfill \end{array}$$(22d) where:(23) $$\begin{array}{cc}\hfill \mathrm{\Theta }(x,t)& =E_1{\mathit{\lambda }}_{\mathit{\sigma }}-{\mathit{\kappa }}_3\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2},\hfill \end{array}$$(23) (24) $$\begin{array}{cc}\hfill {\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}& ={\int }_0^t{\mathit{\lambda }}_{\mathit{\sigma }}(x,z)\mathrm{d}z.\hfill \end{array}$$(24)

The adjoint Equations (21a(21a) $$\begin{array}{cc}& \mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathrm{\Theta }}{\mathit{\partial }x}=0,\hfill \end{array}$$(21a) ) and (21b(21b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}-{\mathit{\kappa }}_1\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2{\mathit{\lambda }}_{\mathit{\sigma }}-h{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}+\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}=0,\hfill \end{array}$$(21b) ) are also mixed with ${\mathit{\lambda }}_u$ and ${\mathit{\lambda }}_{\mathit{\sigma }}$ and PML-endowed. The adjoint solutions depend on state solutions since the adjoint problem is driven by the misfit between computed and observed responses per the boundary condition (22b(22b) $$\begin{array}{cc}\hfill \mathrm{\Theta }(0,t)& =u(0,t)-u_{\mathrm{m}}(0,t),\hfill \end{array}$$(22b) ). As indicated by Equations (22c(22c) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_u(x,T)& =\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}(x,T)=0,\hfill \end{array}$$(22c) ) and (22d(22d) $$\begin{array}{cc}\hfill {\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}(x,T)& ={\mathit{\lambda }}_{\mathit{\sigma }}(x,T)=\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}(x,T)=0,\hfill \end{array}$$(22d) ), the adjoint problem is a final-value problem, as opposed to the initial-value state problem. Using operator notation, the adjoint problem can be rewritten as:(25) $$\begin{array}{cc}\hfill {\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }& =\mathbf{0},\mathrm{f}or0<x<L_t,0\le t<T,\hfill \end{array}$$(25) (26) $$\begin{array}{cc}\hfill {\mathcal{B}}_a\mathit{\lambda }& ={\left[0\left\{u(0,t)-u_{\mathrm{m}}(0,t)\right\}\right]}^{\mathrm{T}},\hfill \end{array}$$(26) (27) $$\begin{array}{cc}\hfill {\mathcal{I}}_a\mathit{\lambda }& =\mathbf{0},\hfill \end{array}$$(27)

where:(28) $$\begin{array}{cc}& {\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }=\left[\begin{array}{c}{\displaystyle \mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathrm{\Theta }}{\mathit{\partial }x}}\hfill \\ {\displaystyle {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}-{\mathit{\kappa }}_1\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2{\mathit{\lambda }}_{\mathit{\sigma }}-h{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}+\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}}\hfill \end{array}\right],\hfill \\ \hfill & {\mathcal{B}}_a\mathit{\lambda }=\left[\begin{array}{c}{\mathit{\lambda }}_u(L_t,t)\hfill \\ \mathrm{\Theta }(0,t)\hfill \end{array}\right],{\mathcal{I}}_a\mathit{\lambda }=\left[\begin{array}{c}{\mathit{\lambda }}_u(x,T)\\ {\displaystyle \frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}(x,T)}\\ {\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}(x,T)\\ {\mathit{\lambda }}_{\mathit{\sigma }}(x,T)\\ {\displaystyle \frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}(x,T)}\end{array}\right],\mathit{\lambda }=\left[\begin{array}{c}{\mathit{\lambda }}_u\\ {\mathit{\lambda }}_{\mathit{\sigma }}\end{array}\right].\hfill \end{array}$$(28)

In Equation (28(28) $$\begin{array}{cc}& {\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }=\left[\begin{array}{c}{\displaystyle \mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathrm{\Theta }}{\mathit{\partial }x}}\hfill \\ {\displaystyle {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}-{\mathit{\kappa }}_1\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2{\mathit{\lambda }}_{\mathit{\sigma }}-h{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}+\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}}\hfill \end{array}\right],\hfill \\ \hfill & {\mathcal{B}}_a\mathit{\lambda }=\left[\begin{array}{c}{\mathit{\lambda }}_u(L_t,t)\hfill \\ \mathrm{\Theta }(0,t)\hfill \end{array}\right],{\mathcal{I}}_a\mathit{\lambda }=\left[\begin{array}{c}{\mathit{\lambda }}_u(x,T)\\ {\displaystyle \frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}(x,T)}\\ {\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}(x,T)\\ {\mathit{\lambda }}_{\mathit{\sigma }}(x,T)\\ {\displaystyle \frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}(x,T)}\end{array}\right],\mathit{\lambda }=\left[\begin{array}{c}{\mathit{\lambda }}_u\\ {\mathit{\lambda }}_{\mathit{\sigma }}\end{array}\right].\hfill \end{array}$$(28) ), ${\mathcal{A}}^{\ast }(\mathit{E})$ represents an adjoint operator on the vector of adjoint variables $\mathit{\lambda }={\left[{\mathit{\lambda }}_u{\mathit{\lambda }}_{\mathit{\sigma }}\right]}^{\mathrm{T}}$, while ${\mathcal{B}}_a$ and ${\mathcal{I}}_a$ are defined as the operators indicating the boundary and final conditions of the adjoint problem, respectively.

3.2.3. The third optimality condition

Lastly, the first variations of $\mathcal{L}$ with respect to the material parameters $E_1$ and $E_2$ must vanish, i.e.:(29a) $$\begin{array}{cc}\hfill {\mathit{\delta }}_{E_1}\mathcal{L}& =0,\hfill \end{array}$$(29a) (29b) $$\begin{array}{cc}\hfill {\mathit{\delta }}_{E_2}\mathcal{L}& =0,\hfill \end{array}$$(29b) which entails the following time-independent control problem:

Find $E_1=E_1(x)$ and $E_2=E_2(x)$, such that:(30a) $$\begin{array}{cc}\hfill -R_{E_1}\frac{{\mathrm{d}}^2{E}_1}{\mathrm{d}{x}^2}+{\int }_0^T\left[-{\mathit{\lambda }}_{\mathit{\sigma }}\frac{\mathit{\partial }u}{\mathit{\partial }x}+\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}+{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\right]\mathrm{d}t& =0,\mathrm{f}or0\le x\le L_t,\hfill \end{array}$$(30a) (30b) $$\begin{array}{cc}\hfill \frac{\mathrm{d}{E}_1}{\mathrm{d}x}(0)=\frac{\mathrm{d}{E}_1}{\mathrm{d}x}(L_t)& =0,\hfill \end{array}$$(30b) (31a) $$\begin{array}{cc}& -R_{E_2}\frac{{\mathrm{d}}^2{E}_2}{\mathrm{d}{x}^2}+{\int }_0^T\left[\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}+2{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\right]\mathrm{d}t=0,\mathrm{f}or0\le x\le L_t,\hfill \end{array}$$(31a) (31b) $$\begin{array}{cc}& \frac{\mathrm{d}{E}_2}{\mathrm{d}x}(0)=\frac{\mathrm{d}{E}_2}{\mathrm{d}x}(L_t)=0.\hfill \end{array}$$(31b) To arrive at Equations (30) and (31a(31a) $$\begin{array}{cc}& -R_{E_2}\frac{{\mathrm{d}}^2{E}_2}{\mathrm{d}{x}^2}+{\int }_0^T\left[\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}+2{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\right]\mathrm{d}t=0,\mathrm{f}or0\le x\le L_t,\hfill \end{array}$$(31a) ), a TN regularization scheme was used in Equation (11(11) $$\begin{array}{c}\hfill \mathcal{J}=\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2),\end{array}$$(11) ). Using operator notation, the control problem (Equations (30) and (31a(31a) $$\begin{array}{cc}& -R_{E_2}\frac{{\mathrm{d}}^2{E}_2}{\mathrm{d}{x}^2}+{\int }_0^T\left[\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}+2{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\right]\mathrm{d}t=0,\mathrm{f}or0\le x\le L_t,\hfill \end{array}$$(31a) )) can be written as:(32) $$\begin{array}{cc}\hfill \mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda })& =\mathbf{0},\mathrm{f}or0\le x\le L_t,\hfill \end{array}$$(32) (33) $$\begin{array}{cc}\hfill {\mathcal{B}}_c\mathit{E}& =\mathbf{0}.\hfill \end{array}$$(33)

where:(34) $$\begin{array}{cc}& \mathcal{G}\mathit{E}=\left[\begin{array}{c}{\displaystyle -R_{E_1}\frac{{\mathrm{d}}^2{E}_1}{\mathrm{d}{x}^2}}\hfill \\ {\displaystyle -R_{E_2}\frac{{\mathrm{d}}^2{E}_2}{\mathrm{d}{x}^2}}\hfill \end{array}\right],\hfill \end{array}$$(34) (35) $$\begin{array}{cc}& \mathcal{S}(\mathit{u},\mathit{\lambda })=\left[\begin{array}{c}{\displaystyle {\int }_0^T\left[-{\mathit{\lambda }}_{\mathit{\sigma }}\frac{\mathit{\partial }u}{\mathit{\partial }x}+\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}+{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\right]\mathrm{d}t}\hfill \\ {\displaystyle {\int }_0^T\left[\left({\mathit{\tau }}_1+{\mathit{\tau }}_2\right)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}+2{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\right]\mathrm{d}t}\hfill \end{array}\right],\hfill \end{array}$$(35) (36) $$\begin{array}{cc}& {\mathcal{B}}_c\mathit{E}={\left[\begin{array}{cccc}{\displaystyle \frac{\mathrm{d}{E}_1}{\mathrm{d}x}(0)}& {\displaystyle \frac{\mathrm{d}{E}_1}{\mathrm{d}x}(L_t)}& {\displaystyle \frac{\mathrm{d}{E}_2}{\mathrm{d}x}(0)}& {\displaystyle \frac{\mathrm{d}{E}_2}{\mathrm{d}x}(L_t)}\end{array}\right]}^{\mathrm{T}}.\hfill \end{array}$$(36)

3.3. KKT conditions

The time-dependent state and adjoint problems and the time-independent control problem construct KKT conditions whose solutions ($\mathit{u},\mathit{\lambda },\mathit{E}$) are the optimal solutions of the inverse problem defined in Equation (11(11) $$\begin{array}{c}\hfill \mathcal{J}=\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2),\end{array}$$(11) ). The KKT conditions can be written as: In Equation (37(37) $$\begin{array}{c}\hfill {\mathfrak{g}}_{\mathit{E}}=\left[\begin{array}{c}{\mathrm{\nabla }}_{E_1}\mathcal{L}\\ {\mathrm{\nabla }}_{E_2}\mathcal{L}\end{array}\right]=\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda }).\end{array}$$(37) ), ${\mathrm{\nabla }}_{\mathit{E}}\mathcal{L}$ indicates ${\left[{\mathrm{\nabla }}_{E_1}\mathcal{L}{\mathrm{\nabla }}_{E_2}\mathcal{L}\right]}^{\mathrm{T}}$, and ${\mathrm{\nabla }}_x\mathcal{L}$ denotes the gradient of the Lagrangian in the form of a Fréchet derivative with respect to a variable $x$. One can use a full-space method to obtain $\mathit{u},\mathit{\lambda }$, and $\mathit{E}$ simultaneously from Equation (37(37) $$\begin{array}{c}\hfill {\mathfrak{g}}_{\mathit{E}}=\left[\begin{array}{c}{\mathrm{\nabla }}_{E_1}\mathcal{L}\\ {\mathrm{\nabla }}_{E_2}\mathcal{L}\end{array}\right]=\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda }).\end{array}$$(37) ). However, the required computational cost is substantial. Therefore, a reduced-space method is introduced to seek optimal solutions of viscoelastic parameters $E_1$ and $E_2$. In particular, the parameters $E_1$ and $E_2$ are iteratively updated in the reduced-space of the control variables $E_1$ and $E_2$ via Gauss–Newton line search algorithms. The reduced gradient required for updating $E_1$ and $E_2$ can be taken as the middle part of the left-hand-side of the vector Equation (37(37) $$\begin{array}{c}\hfill {\mathfrak{g}}_{\mathit{E}}=\left[\begin{array}{c}{\mathrm{\nabla }}_{E_1}\mathcal{L}\\ {\mathrm{\nabla }}_{E_2}\mathcal{L}\end{array}\right]=\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda }).\end{array}$$(37) ) as:(37) $$\begin{array}{c}\hfill {\mathfrak{g}}_{\mathit{E}}=\left[\begin{array}{c}{\mathrm{\nabla }}_{E_1}\mathcal{L}\\ {\mathrm{\nabla }}_{E_2}\mathcal{L}\end{array}\right]=\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda }).\end{array}$$(37)

The reduced gradient can be computed as the following: with the estimated viscoelastic parameters $E_1(x)$ and $E_2(x)$, the state problem (Equations (16(16) $$\begin{array}{cc}\hfill \mathcal{A}(\mathit{E})\mathit{u}& =\mathbf{0},\mathrm{f}or0<x<L_t,0<t\le T,\hfill \end{array}$$(16) )–(18(18) $$\begin{array}{cc}\hfill {\mathcal{I}}_s\mathit{u}& =\mathbf{0},\hfill \end{array}$$(18) )) is solved for $\mathit{u}$; with the state solutions $\mathit{u}$, the adjoint problem (Equations (25(25) $$\begin{array}{cc}\hfill {\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }& =\mathbf{0},\mathrm{f}or0<x<L_t,0\le t<T,\hfill \end{array}$$(25) )–(27(27) $$\begin{array}{cc}\hfill {\mathcal{I}}_a\mathit{\lambda }& =\mathbf{0},\hfill \end{array}$$(27) )) is solved for $\mathit{\lambda }$; and using $\mathit{u}$ and $\mathit{\lambda }$, the reduced gradient is computed per Equation (38(38) $$\begin{array}{c}\hfill \left[\begin{array}{ccccc}\mathcal{B}& & {\mathcal{C}}^{\ast }(\mathit{\lambda })& & {\mathcal{A}}^{\ast }(\mathit{E})\\ \mathcal{C}(\mathit{\lambda })& & \mathcal{F}& & {\mathcal{D}}^{\ast }(\mathit{u})\\ \mathcal{A}(\mathit{E})& & \mathcal{D}(\mathit{u})& & 0\end{array}\right]\left[\begin{array}{c}\overline{\mathit{u}}\\ \overline{\mathit{E}}\\ \overline{\mathit{\lambda }}\end{array}\right]=-\left[\begin{array}{c}{\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }\\ \mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda })\\ \mathcal{A}(\mathit{E})\mathit{u}\end{array}\right].\end{array}$$(38) ). In this procedure, it is noted that the first and last KKT conditions in Equation (37(37) $$\begin{array}{c}\hfill {\mathfrak{g}}_{\mathit{E}}=\left[\begin{array}{c}{\mathrm{\nabla }}_{E_1}\mathcal{L}\\ {\mathrm{\nabla }}_{E_2}\mathcal{L}\end{array}\right]=\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda }).\end{array}$$(37) ) are satisfied at each inversion iteration, while the second KKT condition becomes gradually satisfied as the solution converges to optimal values.

3.4. Newton updates – KKT system

The inversion algorithm requires calculating material property updates ${\overline{E}}_1$ and ${\overline{E}}_2$. In a full-space framework, Newton updates ($\overline{\mathit{u}}$, $\overline{\mathit{\lambda }}$, $\overline{\mathit{E}}$) of the state, adjoint, and control variables should satisfy the following KKT system: where the matrix on the left-hand-side is the full-space Hessian of $\mathcal{L}$, and the vector on the right-hand-side is the full-space gradient defined in Equation (37(37) $$\begin{array}{c}\hfill {\mathfrak{g}}_{\mathit{E}}=\left[\begin{array}{c}{\mathrm{\nabla }}_{E_1}\mathcal{L}\\ {\mathrm{\nabla }}_{E_2}\mathcal{L}\end{array}\right]=\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda }).\end{array}$$(37) ). The Hessian consists of the second Fréchet derivatives of the Lagrangian with respect to $\mathit{u}$, $\mathit{E}$, and $\mathit{\lambda }$ (see Appendix 1 for the expression of the derivatives). Partitioning the Hessian as in Equation (39(39) $$\begin{array}{cc}\hfill \mathcal{B}\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{uu}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{u\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{\mathit{\sigma }u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{\mathit{\sigma }\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}\overline{u}(x,t)\mathit{\delta }(x)\\ 0\end{array}\right],\hfill \end{array}$$(39) ), the KKT system can be written as:(38) $$\begin{array}{c}\hfill \left[\begin{array}{ccccc}\mathcal{B}& & {\mathcal{C}}^{\ast }(\mathit{\lambda })& & {\mathcal{A}}^{\ast }(\mathit{E})\\ \mathcal{C}(\mathit{\lambda })& & \mathcal{F}& & {\mathcal{D}}^{\ast }(\mathit{u})\\ \mathcal{A}(\mathit{E})& & \mathcal{D}(\mathit{u})& & 0\end{array}\right]\left[\begin{array}{c}\overline{\mathit{u}}\\ \overline{\mathit{E}}\\ \overline{\mathit{\lambda }}\end{array}\right]=-\left[\begin{array}{c}{\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }\\ \mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda })\\ \mathcal{A}(\mathit{E})\mathit{u}\end{array}\right].\end{array}$$(38)

In Equation (40(40) $$\begin{array}{cc}\hfill \mathcal{C}(\mathit{\lambda })\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{E_{1}u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_1\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{E_{2}u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_2\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}{\displaystyle -{\mathit{\lambda }}_{\mathit{\sigma }}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \\ {\displaystyle -({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }{x}^2}-2{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \end{array}\right],\hfill \end{array}$$(40) ), the operators in the Hessian are defined as:(39) $$\begin{array}{cc}\hfill \mathcal{B}\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{uu}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{u\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{\mathit{\sigma }u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{\mathit{\sigma }\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}\overline{u}(x,t)\mathit{\delta }(x)\\ 0\end{array}\right],\hfill \end{array}$$(39) (40) $$\begin{array}{cc}\hfill \mathcal{C}(\mathit{\lambda })\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{E_{1}u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_1\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{E_{2}u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_2\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}{\displaystyle -{\mathit{\lambda }}_{\mathit{\sigma }}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \\ {\displaystyle -({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }{x}^2}-2{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \end{array}\right],\hfill \end{array}$$(40) (41) $$\begin{array}{cc}\hfill {\mathcal{C}}^{\ast }(\mathit{\lambda })\overline{\mathit{E}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{u{E}_1}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{u{E}_2}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{\mathit{\sigma }{E}_1}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{\mathit{\sigma }{E}_2}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}{\overline{E}}_1\hfill \\ {\overline{E}}_2\hfill \end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}{\displaystyle \frac{\mathit{\partial }}{\mathit{\partial }x}\left({\mathit{\lambda }}_{\mathit{\sigma }}{\overline{E}}_1\right)-({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }}{\mathit{\partial }x}\left(({\overline{E}}_1+{\overline{E}}_2)\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\right)+{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }}{\mathit{\partial }x}\left(({\overline{E}}_1+2{\overline{E}}_2)\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\right)}\\ {\displaystyle 0}\end{array}\right],\hfill \end{array}$$(41) (42) $$\begin{array}{cc}\hfill \mathcal{D}(\mathit{u})\overline{\mathit{E}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{{\mathit{\lambda }}_u{E}_1}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{{\mathit{\lambda }}_u{E}_2}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}{E}_1}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}{E}_2}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}{\overline{E}}_1\hfill \\ {\overline{E}}_2\hfill \end{array}\right]\hfill \\ & =\left[\begin{array}{c}0\\ {\displaystyle {\overline{E}}_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-({\mathit{\tau }}_1+{\mathit{\tau }}_2)({\overline{E}}_1+{\overline{E}}_2)\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\tau }}_1{\mathit{\tau }}_2({\overline{E}}_1+2{\overline{E}}_2)\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}}\end{array}\right],\hfill \end{array}$$(42) (43) $$\begin{array}{cc}\hfill {\mathcal{D}}^{\ast }(\mathit{u})\overline{\mathit{\lambda }}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{E_1{\mathit{\lambda }}_u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_1{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{E_2{\mathit{\lambda }}_u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_2{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}{\overline{\mathit{\lambda }}}_u\hfill \\ {\overline{\mathit{\lambda }}}_{\mathit{\sigma }}\hfill \end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}{\displaystyle -{\int }_0^T{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}t+{\int }_0^T({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}t+{\int }_0^T{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\mathrm{d}t}\\ {\displaystyle {\int }_0^T({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}t+2{\int }_0^T{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\mathrm{d}t}\end{array}\right],\hfill \end{array}$$(43) (44) $$\begin{array}{cc}\hfill \mathcal{F}\overline{\mathit{E}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{E_1{E}_1}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_1{E}_2}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{E_2{E}_1}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_2{E}_2}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}{\overline{E}}_1\hfill \\ {\overline{E}}_2\hfill \end{array}\right]\hfill \\ & =\left\{\begin{array}{cc}{\displaystyle \left[\begin{array}{c}{\displaystyle -R_{E_1}\frac{{\mathrm{d}}^2{\overline{E}}_1}{\mathrm{d}{x}^2}}\\ {\displaystyle -R_{E_2}\frac{{\mathrm{d}}^2{\overline{E}}_2}{\mathrm{d}{x}^2}}\end{array}\right],}\hfill & \text{for TN regularization,}\hfill \\ \left[\begin{array}{c}{\displaystyle -e{R}_{E_1}\frac{\mathrm{d}}{\mathrm{d}x}\left\{{\left({\left(\frac{\mathrm{d}{E}_1}{\mathrm{d}x}\right)}^2+e\right)}^{-\frac{3}{2}}\frac{\mathrm{d}{\overline{E}}_1}{\mathrm{d}x}\right\}}\\ {\displaystyle -e{R}_{E_2}\frac{\mathrm{d}}{\mathrm{d}x}\left\{{\left({\left(\frac{\mathrm{d}{E}_2}{\mathrm{d}x}\right)}^2+e\right)}^{-\frac{3}{2}}\frac{\mathrm{d}{\overline{E}}_2}{\mathrm{d}x}\right\}}\end{array}\right],\hfill & \text{for TV regularization.}\hfill \end{array}\right.\hfill \end{array}$$(44) (45) $$\begin{array}{cc}\hfill \overline{\mathit{u}}& =\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right],\overline{\mathit{\lambda }}=\left[\begin{array}{c}{\overline{\mathit{\lambda }}}_u\hfill \\ {\overline{\mathit{\lambda }}}_{\mathit{\sigma }}\hfill \end{array}\right],\overline{\mathit{E}}=\left[\begin{array}{c}{\overline{E}}_1\hfill \\ {\overline{E}}_2\hfill \end{array}\right].\hfill \end{array}$$(45)

At each inversion iteration, the updates $\overline{\mathit{u}}$, $\overline{\mathit{E}}$, and $\overline{\mathit{\lambda }}$ can be computed from Equation (39(39) $$\begin{array}{cc}\hfill \mathcal{B}\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{uu}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{u\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{\mathit{\sigma }u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{\mathit{\sigma }\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}\overline{u}(x,t)\mathit{\delta }(x)\\ 0\end{array}\right],\hfill \end{array}$$(39) ) once the state and adjoint solutions $\mathit{u}$ and $\mathit{\lambda }$ are obtained. In general, the Newton Hessian is known to be indefinite if the solution is not close to the optimum. Therefore, finding a suitable search direction (or update) may become onerous, especially at early inversion iterations. To mitigate this problem, it is common to use a Gauss–Newton Hessian in the KKT system that is positive definite even when the solution is far from the optimum.[15] The Gauss–Newton approximation is equivalent to ignoring terms that depend on adjoint variables $\mathit{\lambda }$ in the Newton Hessian.[15] Introducing the Gauss–Newton Hessian, Equation (40(40) $$\begin{array}{cc}\hfill \mathcal{C}(\mathit{\lambda })\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{E_{1}u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_1\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{E_{2}u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{E_2\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}{\displaystyle -{\mathit{\lambda }}_{\mathit{\sigma }}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \\ {\displaystyle -({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }{x}^2}-2{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}\hfill \end{array}\right],\hfill \end{array}$$(40) ) can be modified as follows:(46) $$\begin{array}{c}\hfill \left[\begin{array}{ccccc}\mathcal{B}& & 0& & {\mathcal{A}}^{\ast }(\mathit{E})\\ 0& & \mathcal{F}& & {\mathcal{D}}^{\ast }(\mathit{u})\\ \mathcal{A}(\mathit{E})& & \mathcal{D}(\mathit{u})& & 0\end{array}\right]\left[\begin{array}{c}\overline{\mathit{u}}\\ \overline{\mathit{E}}\\ \overline{\mathit{\lambda }}\end{array}\right]=-\left[\begin{array}{c}{\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }\\ \mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda })\\ \mathcal{A}(\mathit{E})\mathit{u}\end{array}\right].\end{array}$$(46)

To obtain the material updates $\overline{\mathit{E}}$ from Equation (48(48) $$\begin{array}{cc}\hfill {\mathcal{B}}_s\overline{\mathit{u}}& =\mathbf{0},\hfill \end{array}$$(48) ), the reduced-space method is again used. Notice that the first and last terms of the right-hand-side vector in Equation (48(48) $$\begin{array}{cc}\hfill {\mathcal{B}}_s\overline{\mathit{u}}& =\mathbf{0},\hfill \end{array}$$(48) ) vanish, since the adjoint and state equations are satisfied at each inversion iteration. Therefore, the third block Equation in (48(48) $$\begin{array}{cc}\hfill {\mathcal{B}}_s\overline{\mathit{u}}& =\mathbf{0},\hfill \end{array}$$(48) ) can be solved for $\overline{\mathit{u}}$ given the state solution $\mathit{u}$ and the estimate of the material updates $\overline{\mathit{E}}$. Using the computed $\overline{\mathit{u}}$, the first block equation can be solved for $\overline{\mathit{\lambda }}$. Then, the second block equation in (48(48) $$\begin{array}{cc}\hfill {\mathcal{B}}_s\overline{\mathit{u}}& =\mathbf{0},\hfill \end{array}$$(48) ) can be solved for $\overline{\mathit{E}}$ using the solution $\overline{\mathit{\lambda }}$. This block elimination process is equivalent to solving the following incremental state, adjoint, and control problems:

Incremental state problem

Find $\overline{\mathit{u}}(x,t)\equiv {\left[\overline{u}(x,t)\overline{\mathit{\sigma }}(x,t)\right]}^{\mathrm{T}}$ such that:(47) $$\begin{array}{cc}\hfill \mathcal{A}\overline{\mathit{u}}& =-\mathcal{D}\overline{\mathit{E}},\mathrm{f}or0<x<L_t,0<t\le T,\hfill \end{array}$$(47) (48) $$\begin{array}{cc}\hfill {\mathcal{B}}_s\overline{\mathit{u}}& =\mathbf{0},\hfill \end{array}$$(48) (49) $$\begin{array}{cc}\hfill {\mathcal{I}}_s\overline{\mathit{u}}& =\mathbf{0}.\hfill \end{array}$$(49) Incremental adjoint problem

Find $\overline{\mathit{\lambda }}(x,t)\equiv {\left[{\overline{\mathit{\lambda }}}_u(x,t){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}(x,t)\right]}^{\mathrm{T}}$ such that:(50) $$\begin{array}{cc}\hfill {\mathcal{A}}^{\ast }\overline{\mathit{\lambda }}& =-\mathcal{B}\overline{\mathit{u}},\mathrm{f}or0<x<L_t,0\le t<T,\hfill \end{array}$$(50) (51) $$\begin{array}{cc}\hfill {\mathcal{B}}_a\overline{\mathit{\lambda }}& =\mathbf{0},\hfill \end{array}$$(51) (52) $$\begin{array}{cc}\hfill {\mathcal{I}}_a\overline{\mathit{\lambda }}& =\mathbf{0}.\hfill \end{array}$$(52) Incremental control problem

Find $\overline{\mathit{E}}(x,t)\equiv {\left[{\overline{E}}_1(x,t){\overline{E}}_2(x,t)\right]}^{\mathrm{T}}$ such that:(53) $$\begin{array}{cc}\hfill \mathcal{F}\overline{\mathit{E}}+{\mathcal{D}}^{\ast }\overline{\mathit{\lambda }}& =-\left[\mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda })\right],\mathrm{f}or0<x<L_t,\hfill \end{array}$$(53) (54) $$\begin{array}{cc}\hfill {\mathcal{B}}_c\overline{\mathit{E}}& =\mathbf{0}.\hfill \end{array}$$(54)

In Equations (49(49) $$\begin{array}{cc}\hfill {\mathcal{I}}_s\overline{\mathit{u}}& =\mathbf{0}.\hfill \end{array}$$(49) )–(56(56) $$\begin{array}{cc}\hfill \overline{\mathit{\lambda }}& ={\left({\mathcal{A}}^{\ast }\right)}^{-1}\mathcal{B}{\mathcal{A}}^{-1}\mathcal{D}\overline{\mathit{E}},\hfill \end{array}$$(56) ), differential operators $\mathcal{A}$, ${\mathcal{B}}_s$, ${\mathcal{I}}_s$, ${\mathcal{A}}^{\ast }$, ${\mathcal{B}}_a$, ${\mathcal{I}}_a$, $\mathcal{G}$, and ${\mathcal{B}}_c$ are the same as the state, adjoint, and control operators in the KKT conditions. The incremental state and adjoint problems can also be solved by using a mixed finite element method. Using operator notation, solutions of the incremental state and adjoint problems can be expressed as:(55) $$\begin{array}{cc}\hfill \overline{\mathit{u}}& =-{\mathcal{A}}^{-1}\mathcal{D}\overline{\mathit{E}},\hfill \end{array}$$(55) (56) $$\begin{array}{cc}\hfill \overline{\mathit{\lambda }}& ={\left({\mathcal{A}}^{\ast }\right)}^{-1}\mathcal{B}{\mathcal{A}}^{-1}\mathcal{D}\overline{\mathit{E}},\hfill \end{array}$$(56)

in which $\mathcal{A}$ and ${\left({\mathcal{A}}^{\ast }\right)}^{-1}$ denote inverse operators of the state and adjoint problems, respectively. Substituting Equation (58(58) $$\begin{array}{c}\hfill W_{\mathit{E}}=\mathcal{F}+{\mathcal{D}}^{\ast }{\left({\mathcal{A}}^{\ast }\right)}^{-1}\mathcal{B}{\mathcal{A}}^{-1}\mathcal{D}.\end{array}$$(58) ) into Equation (55(55) $$\begin{array}{cc}\hfill \overline{\mathit{u}}& =-{\mathcal{A}}^{-1}\mathcal{D}\overline{\mathit{E}},\hfill \end{array}$$(55) ) yields the following reduced KKT system:(57) $$\begin{array}{c}\hfill W_{\mathit{E}}\overline{\mathit{E}}=-{\mathfrak{g}}_{\mathit{E}},\end{array}$$(57)

where $W_{\mathit{E}}$ denotes the differential operator of the reduced Hessian in the space of control variables $\mathit{E}$, which can be expressed as:(58) $$\begin{array}{c}\hfill W_{\mathit{E}}=\mathcal{F}+{\mathcal{D}}^{\ast }{\left({\mathcal{A}}^{\ast }\right)}^{-1}\mathcal{B}{\mathcal{A}}^{-1}\mathcal{D}.\end{array}$$(58)

Equation (59(59) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}\right)\mathrm{d}x+{\int }_0^{L_t}\frac{\mathrm{d}w}{\mathrm{d}x}\mathit{\sigma }\mathrm{d}x=-w(0)p(t),\hfill \end{array}$$(59) ) indicates a continuous form of the reduced KKT system. To represent the system in a discrete form, and thus to obtain the material updates at discretized grid points, it is necessary to construct the reduced Hessian in matrix form. However, such a matrix representation of the reduced Hessian is computationally impractical, since its size becomes excessively large once discretization is introduced in both space and time. In this work, a Krylov subspace method is used to iteratively determine material updates $\overline{\mathit{E}}$ from Equation (59(59) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}\right)\mathrm{d}x+{\int }_0^{L_t}\frac{\mathrm{d}w}{\mathrm{d}x}\mathit{\sigma }\mathrm{d}x=-w(0)p(t),\hfill \end{array}$$(59) ) without requiring the construction of a discrete Hessian.

4. Finite element implementation of the state and adjoint problems

To obtain solutions $\mathit{u}={[u\mathit{\sigma }]}^{\mathrm{T}}$ of the state problem (Equations (4)–(5)), a mixed finite element method can be employed, where both displacement ($u$) and stress ($\mathit{\sigma }$) are treated as independent unknowns.[12] The variational form can be obtained by multiplying Equations (4a(4a) $$\begin{array}{cc}& \frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}=0,\hfill \end{array}$$(4a) ) and (4b(4b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}=0,\hfill \end{array}$$(4b) ) by appropriate test functions $w(x)$ and $q(x)$, and then integrating over the entire domain, $\mathrm{\Omega }=\left\{x:0<x<L_t\right\}$. Then, using the definition of function spaces $H^1(\mathrm{\Omega }):=\left\{h(x):{\int }_{\mathrm{\Omega }}{h}^2+{(h^{\prime })}^2\mathrm{d}x<\infty \right\}$ and $L^2(\mathrm{\Omega }):=\left\{h(x):{\int }_{\mathrm{\Omega }}{h}^2\mathrm{d}x<\infty \right\}$, the weak form of Equation (4) can be stated as: find $u\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))\times (0,T]$ and $\mathit{\sigma }\in L^2(\mathrm{\Omega })\times (0,T]$ such that the following equations are satisfied for all $w\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))$ and $q\in L^2(\mathrm{\Omega })$:(59) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}\right)\mathrm{d}x+{\int }_0^{L_t}\frac{\mathrm{d}w}{\mathrm{d}x}\mathit{\sigma }\mathrm{d}x=-w(0)p(t),\hfill \end{array}$$(59) (60) $$\begin{array}{cc}& {\int }_0^{L_t}{\mathit{\tau }}_1{\mathit{\tau }}_{2}q\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}\mathrm{d}x+{\int }_0^{L_t}{\mathit{\kappa }}_{1}q\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}\mathrm{d}x+{\int }_0^{L_t}{\mathit{\kappa }}_{2}q\mathit{\sigma }\mathrm{d}x+{\int }_0^{L_t}hq\tilde{\mathit{\sigma }}\mathrm{d}x\hfill \\ \hfill & -{\int }_0^{L_t}{E}_{1}q\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}x-{\int }_0^{L_t}{\mathit{\kappa }}_3{\mathit{\tau }}_1{\mathit{\tau }}_{2}q\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\mathrm{d}x-{\int }_0^{L_t}{\mathit{\kappa }}_{4}q\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}\mathrm{d}x=0.\hfill \end{array}$$(60)

Based on the Galerkin method, the trial functions $u_h\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))\times (0,T]$ and ${\mathit{\sigma }}_h\in L^2(\mathrm{\Omega })\times (0,T]$ can be spatially discretized as:(61) $$\begin{array}{cc}\hfill u(x,t)& \cong u_h(x,t)=\sum _{j=1}^N{\mathit{\varphi }}_j(x)u_j(t),\hfill \end{array}$$(61) (62) $$\begin{array}{cc}\hfill \mathit{\sigma }(x,t)& \cong {\mathit{\sigma }}_h(x,t)=\sum _{j=1}^N{\mathit{\psi }}_j(x){\mathit{\sigma }}_j(t),\hfill \end{array}$$(62)

where ${\mathit{\varphi }}_j$ and ${\mathit{\psi }}_j$ are basis functions associated with nodal displacement $u_j$ and nodal stress ${\mathit{\sigma }}_j$, respectively. Similarly, the test functions $w_h\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))$ and $q_h\in L^2(\mathrm{\Omega })$ can be expressed as:(63) $$\begin{array}{cc}\hfill w(x)& \cong w_h(x)=\sum _{i=1}^N{w}_i{\mathit{\varphi }}_i(x),\hfill \end{array}$$(63) (64) $$\begin{array}{cc}\hfill q(x)& \cong q_h(x)=\sum _{i=1}^N{q}_i{\mathit{\psi }}_i(x).\hfill \end{array}$$(64)

To ensure the stability of solutions from the mixed finite element method, the basic functions ${\mathit{\varphi }}_i$ and ${\mathit{\psi }}_i$ were chosen to be piecewise quadratic. Introducing Equations (63(63) $$\begin{array}{cc}\hfill w(x)& \cong w_h(x)=\sum _{i=1}^N{w}_i{\mathit{\varphi }}_i(x),\hfill \end{array}$$(63) )–(66(66) $$\begin{array}{cc}& {\int }_0^{L_t}{\mathit{\tau }}_1{\mathit{\tau }}_{2}q\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\mathrm{d}x-{\int }_0^{L_t}{\mathit{\kappa }}_{1}q\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\mathrm{d}x+{\int }_0^{L_t}{\mathit{\kappa }}_{2}q{\mathit{\lambda }}_{\mathit{\sigma }}\mathrm{d}x-{\int }_0^{L_t}hq{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}\mathrm{d}x+{\int }_0^{L_t}q\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}\mathrm{d}x=0.\hfill \\ \hfill \end{array}$$(66) ) to Equations (61(61) $$\begin{array}{cc}\hfill u(x,t)& \cong u_h(x,t)=\sum _{j=1}^N{\mathit{\varphi }}_j(x)u_j(t),\hfill \end{array}$$(61) ) and (62(62) $$\begin{array}{cc}\hfill \mathit{\sigma }(x,t)& \cong {\mathit{\sigma }}_h(x,t)=\sum _{j=1}^N{\mathit{\psi }}_j(x){\mathit{\sigma }}_j(t),\hfill \end{array}$$(62) ) yields semi-discrete equations of motion that can be solved using a Newmark-$\mathit{\beta }$ method.

For the adjoint problem defined in Equations (21)–(22), the mixed finite element method can also be used to obtain approximate solutions for ${\mathit{\lambda }}_v$ and ${\mathit{\lambda }}_{\mathit{\sigma }}$. The variational form can be obtained by multiplying (21a(21a) $$\begin{array}{cc}& \mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathrm{\Theta }}{\mathit{\partial }x}=0,\hfill \end{array}$$(21a) ) and (21b(21b) $$\begin{array}{cc}& {\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}-{\mathit{\kappa }}_1\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2{\mathit{\lambda }}_{\mathit{\sigma }}-h{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}+\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}=0,\hfill \end{array}$$(21b) ) by individual test functions $w(x)$ and $q(x)$, and then integrating over the entire domain $\mathrm{\Omega }$: find ${\mathit{\lambda }}_u\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))\times [0,T)$ and ${\mathit{\lambda }}_{\mathit{\sigma }}\in L^2(\mathrm{\Omega })\times [0,T)$ such that the following equations are satisfied for all $w\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))$ and $q\in L^2(\mathrm{\Omega })$:(65) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}\right)\mathrm{d}x-{\int }_0^{L_t}\frac{\mathit{\partial }w}{\mathit{\partial }x}\mathrm{\Theta }\mathrm{d}x=-w(0)\left[u(0,t)-u_m(0,t)\right],\hfill \end{array}$$(65) (66) $$\begin{array}{cc}& {\int }_0^{L_t}{\mathit{\tau }}_1{\mathit{\tau }}_{2}q\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\mathrm{d}x-{\int }_0^{L_t}{\mathit{\kappa }}_{1}q\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\mathrm{d}x+{\int }_0^{L_t}{\mathit{\kappa }}_{2}q{\mathit{\lambda }}_{\mathit{\sigma }}\mathrm{d}x-{\int }_0^{L_t}hq{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}\mathrm{d}x+{\int }_0^{L_t}q\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}\mathrm{d}x=0.\hfill \\ \hfill \end{array}$$(66)

In deriving Equation (67(67) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_u(x,t)& \cong {({\mathit{\lambda }}_u)}_h(x,t)=\sum _{j=1}^N{\mathit{\varphi }}_j(x){\mathit{\lambda }}_{uj}(t),\hfill \end{array}$$(67) ), a boundary condition (Equation (22b(22b) $$\begin{array}{cc}\hfill \mathrm{\Theta }(0,t)& =u(0,t)-u_{\mathrm{m}}(0,t),\hfill \end{array}$$(22b) )) was used as well as the condition on the test function, $w(L_t)=0$. Using the same basis functions used in the state problem, the trial functions ${({\mathit{\lambda }}_u)}_h\in L^2(\mathrm{\Omega };H^1(\mathrm{\Omega }))\times [0,T)$ and ${({\mathit{\lambda }}_{\mathit{\sigma }})}_h\in L^2(\mathrm{\Omega })\times [0,T)$ are spatially discretized as:(67) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_u(x,t)& \cong {({\mathit{\lambda }}_u)}_h(x,t)=\sum _{j=1}^N{\mathit{\varphi }}_j(x){\mathit{\lambda }}_{uj}(t),\hfill \end{array}$$(67) (68) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_{\mathit{\sigma }}(x,t)& \cong {({\mathit{\lambda }}_{\mathit{\sigma }})}_h(x,t)=\sum _{j=1}^N{\mathit{\psi }}_j(x){\mathit{\lambda }}_{\mathit{\sigma }j}(t),\hfill \end{array}$$(68)

Introducing Equations (69(69) $$\begin{array}{c}\hfill {\mathbf{W}}_{\mathit{E}}\overline{\mathbf{E}}=-{\mathbf{g}}_{\mathit{E}},\end{array}$$(69) ) and (70(70) $$\begin{array}{c}\hfill ||{\left({\mathbf{W}}_{\mathit{E}}\right)}_k{\overline{\mathbf{E}}}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_k||\le {\mathit{\eta }}_k||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k||,\end{array}$$(70) ) into the variational form, while also using Equations (65(65) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}\right)\mathrm{d}x-{\int }_0^{L_t}\frac{\mathit{\partial }w}{\mathit{\partial }x}\mathrm{\Theta }\mathrm{d}x=-w(0)\left[u(0,t)-u_m(0,t)\right],\hfill \end{array}$$(65) ) and (66(66) $$\begin{array}{cc}& {\int }_0^{L_t}{\mathit{\tau }}_1{\mathit{\tau }}_{2}q\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\mathrm{d}x-{\int }_0^{L_t}{\mathit{\kappa }}_{1}q\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\mathrm{d}x+{\int }_0^{L_t}{\mathit{\kappa }}_{2}q{\mathit{\lambda }}_{\mathit{\sigma }}\mathrm{d}x-{\int }_0^{L_t}hq{\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}\mathrm{d}x+{\int }_0^{L_t}q\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }x}\mathrm{d}x=0.\hfill \\ \hfill \end{array}$$(66) ), one can obtain semi-discrete equations of motion for the adjoint problem. Recall that the adjoint problem is a final-value problem per Equations (22c(22c) $$\begin{array}{cc}\hfill {\mathit{\lambda }}_u(x,T)& =\frac{\mathit{\partial }{\mathit{\lambda }}_u}{\mathit{\partial }t}(x,T)=0,\hfill \end{array}$$(22c) ) and (22d(22d) $$\begin{array}{cc}\hfill {\tilde{\mathit{\lambda }}}_{\mathit{\sigma }}(x,T)& ={\mathit{\lambda }}_{\mathit{\sigma }}(x,T)=\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}(x,T)=0,\hfill \end{array}$$(22d) ). Consequently, a reverse time integration is needed when solving the adjoint semi-discrete form.

5. Inversion process

Once the state and adjoint solutions are obtained at each inversion iteration, the material updates $\overline{\mathit{E}}$ can be calculated by resolving Equation (59(59) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}\right)\mathrm{d}x+{\int }_0^{L_t}\frac{\mathrm{d}w}{\mathrm{d}x}\mathit{\sigma }\mathrm{d}x=-w(0)p(t),\hfill \end{array}$$(59) ). In a discrete form, Equation (59(59) $$\begin{array}{cc}& {\int }_0^{L_t}w\left(\mathit{\rho }\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}\right)\mathrm{d}x+{\int }_0^{L_t}\frac{\mathrm{d}w}{\mathrm{d}x}\mathit{\sigma }\mathrm{d}x=-w(0)p(t),\hfill \end{array}$$(59) ) can be written as:(69) $$\begin{array}{c}\hfill {\mathbf{W}}_{\mathit{E}}\overline{\mathbf{E}}=-{\mathbf{g}}_{\mathit{E}},\end{array}$$(69)

where ${\mathbf{W}}_{\mathit{E}}$ and ${\mathbf{g}}_{\mathit{E}}$ denote the discrete form of the reduced Hessian and the reduced gradient, respectively. $\overline{\mathbf{E}}\equiv {\left[{\overline{\mathbf{E}}}_1^{\mathrm{T}}{\overline{\mathbf{E}}}_2^{\mathrm{T}}\right]}^{\mathrm{T}}$ is the vector of nodal viscoelastic parameter updates. The number of elements in ${\mathbf{g}}_{\mathit{E}}$ and $\overline{\mathbf{E}}$ equals twice the number of nodes. The gradient vector ${\mathbf{g}}_{\mathit{E}}$ can be constructed by evaluating Equation (38(38) $$\begin{array}{c}\hfill \left[\begin{array}{ccccc}\mathcal{B}& & {\mathcal{C}}^{\ast }(\mathit{\lambda })& & {\mathcal{A}}^{\ast }(\mathit{E})\\ \mathcal{C}(\mathit{\lambda })& & \mathcal{F}& & {\mathcal{D}}^{\ast }(\mathit{u})\\ \mathcal{A}(\mathit{E})& & \mathcal{D}(\mathit{u})& & 0\end{array}\right]\left[\begin{array}{c}\overline{\mathit{u}}\\ \overline{\mathit{E}}\\ \overline{\mathit{\lambda }}\end{array}\right]=-\left[\begin{array}{c}{\mathcal{A}}^{\ast }(\mathit{E})\mathit{\lambda }\\ \mathcal{G}\mathit{E}+\mathcal{S}(\mathit{u},\mathit{\lambda })\\ \mathcal{A}(\mathit{E})\mathit{u}\end{array}\right].\end{array}$$(38) ) at each nodal point using the discrete solutions of state and adjoint problems under the current estimate of viscoelastic parameters. To resolve Equation (71(71) $$\begin{array}{c}\hfill {\mathit{\eta }}_k=\frac{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}-{\left({\mathbf{W}}_{\mathit{E}}\right)}_{k-1}{\overline{\mathbf{E}}}_{k-1}||}{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}||},\end{array}$$(71) ) for $\overline{\mathbf{E}}$, a Krylov subspace method is utilized instead of directly solving it by factorization. The Krylov subspace method avoids factorization of system matrix by using matrix-vector multiplication in the iterative process to yield solutions.[16–21] In this work, the Conjugate Gradient (CG) method is adopted as the Krylov subspace method.

Outlined in Table 1 is the procedure to compute ${\mathbf{W}}_{\mathit{E}}$ multiplied by an arbitrary vector $\widehat{\mathbf{E}}$. This procedure will be referred to as Algorithm 1 in the remaining discussion. The CG method is outlined in Table 2. The CG method (Algorithm 2) iteratively computes viscoelastic parameter updates $\overline{\mathbf{E}}$ within each inversion iteration. Step 7 in the algorithm requires the calculation of the reduced Gauss–Newton Hessian multiplied by $\widehat{\mathbf{E}}$, which is computed with Algorithm 1.

Table 1. Algorithm 1: Calculation of a reduced Gauss-Newton Hessian multiplied by a vector $\widehat{\mathbf{E}}$.

Step 6 in Algorithm 2 states a termination condition of the CG iteration. Specifically, the update $\overline{\mathbf{E}}$ of nodal viscoelastic parameters at the $k$-th inversion iteration should satisfy the following condition:(70) $$\begin{array}{c}\hfill ||{\left({\mathbf{W}}_{\mathit{E}}\right)}_k{\overline{\mathbf{E}}}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_k||\le {\mathit{\eta }}_k||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k||,\end{array}$$(70)

In (72(72) $$\begin{array}{c}\hfill {\mathbf{E}}_{k+1}={\mathbf{E}}_k+\mathit{\alpha }{\overline{\mathbf{E}}}_k,\end{array}$$(72) ), ${\mathit{\eta }}_k$ is a forcing term that can be determined by:(71) $$\begin{array}{c}\hfill {\mathit{\eta }}_k=\frac{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}-{\left({\mathbf{W}}_{\mathit{E}}\right)}_{k-1}{\overline{\mathbf{E}}}_{k-1}||}{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}||},\end{array}$$(71)

Alternative forms are also possible.[22] ${\mathit{\eta }}_k$ is typically taken as ${\mathit{\eta }}_k<1$, and for ${\mathit{\eta }}_k\to 0$, the convergence rate has been reported to be super-linear.[15] Using the Gauss–Newton update ${\overline{\mathbf{E}}}_k$ at the $k$-th inversion iteration, the nodal viscoelastic parameter vector ${\mathbf{E}}_k$ can be updated as:(72) $$\begin{array}{c}\hfill {\mathbf{E}}_{k+1}={\mathbf{E}}_k+\mathit{\alpha }{\overline{\mathbf{E}}}_k,\end{array}$$(72)

where $\mathit{\alpha }$ is the step length in the direction of ${\overline{\mathbf{E}}}_k$. In this work, an inexact line search method was used to determine a step length that achieves sufficient reduction in the objective functional $\mathcal{J}$. A backtracking line search algorithm to perform the inexact line search is outlined in Table 3.

Table 2. Algorithm 2: Conjugate Gradient method.

Table 3. Algorithm 3: backtracking line search procedure.

Table 4. Algorithm 4: Gauss-Newton inversion process.

In Algorithm 3, $\overline{\mathit{\alpha }}$ denotes the initial step length and is typically chosen to be $\overline{\mathit{\alpha }}=1$ for the Newton method. Step 4 specifies the Armijo condition enforcing the sufficient decrease in the objective functional $\mathcal{J}$. If the Armijo condition is violated, an appropriate step length will be determined after a finite number of trials by setting $\mathit{\alpha }\leftarrow \mathit{\rho }\mathit{\alpha }$. $\mathit{\rho }=0.5$ is a typical value and $\mathit{\mu }$ is usually chosen to be small.[23] $\mathit{\mu }={10}^{-10}$ was used herein. The entire Gauss–Newton inversion process discussed so far is summarized in Table 4.

6. Numerical examples

6.1. Smoothly-varying profiles

To provide numerical examples, a heterogeneous semi-infinite solid medium with smoothly-varying material properties was defined. The medium was modeled as a one-dimensional PML-truncated domain, with the regular domain occupying $0m\le x<100m$, and the PML placed at $100m\le x\le 110m$, as shown in Figure 3. The medium’s viscoelastic properties were described using the Generalized Maxwell Body. Figure 4(a) and (b) depict the target $E_1$ and $E_2$ profiles within the PML-truncated domain. Both profiles show bell-shaped variations along the depth direction with peak values of $E_1=15\mathrm{MPa}$ and $E_2=4\mathrm{MPa}$ at $x=50\mathrm{m}$, respectively, as represented by the following equations:(73) $$\begin{array}{cc}\hfill E_1(x)& ={10}^7+\left(0.5\times {10}^7\right)e^{-\frac{{(x-50)}^2}{200}}{\mathrm{N}/\mathrm{m}}^2,\hfill \end{array}$$(73) (74) $$\begin{array}{cc}\hfill E_2(x)& =\left(3\times {10}^6\right)+{10}^6{e}^{-\frac{{(x-50)}^2}{450}}{\mathrm{N}/\mathrm{m}}^2.\hfill \end{array}$$(74)

With the $E_1$ and $E_2$ profiles, the target damping parameter $Q$ can be determined by Equation (3(3) $$\begin{array}{cc}\hfill \frac{E_2}{E_1+2E_2}& =\mathit{\alpha }{Q}^{\mathit{\beta }},\hfill \end{array}$$(3) ). Figure 4(c) presents the target $Q$ distribution in the entire domain. The density $\mathit{\rho }$ is $2000{\mathrm{kg}/\mathrm{m}}^3$, and a reflection coefficient of $|R|={10}^{-4}$ was imposed on the PML.

Figure 3. (a) Heterogeneous semi-infinite solid medium with smoothly-varying material properties; (b) Schematic of a one-dimensional PML-truncated domain.

Figure 4. Smooth target profiles for $E_1$, $E_2$, and $Q$ factor.

Excitation of the medium was applied on the surface ($x=0$) by two Gaussian pulse-type loads with maximum frequencies of about 9 and 30 Hz, respectively. Figure 5 shows the time history and the frequency spectrum of the excitations. Within this frequency range, it was shown that the two relaxation time parameters ${\mathit{\tau }}_1$ and ${\mathit{\tau }}_2$ are approximately 0.45 and 0.02, respectively.[10] Figure 6 shows the measured displacement responses $u(0,t)$ on the surface, which was synthesized by solving the forward problem using the two Gaussian pulse loads. For the solution of all inverse problems in this section, quadratic elements with lengths of $L_e=0.25$ m were used.

Figure 5. Time history of surface stress load $p(t)$ ($\mathrm{\Delta }t=0.002$ s) and its frequency spectrum.

Figure 6. Measured displacement responses $u(0,t)$ for the surface loads $p(t)$ with (a) $f_{\max }=9$ Hz and (b) $f_{\max }=30$ Hz.

Figure 7. Initial guess, target, and inverted viscoelastic material profiles with smooth variation.

Figure 7(a) and (b) show the reconstructed profiles of $E_1$ and $E_2$ of the PML-truncated viscoelastic medium using Tikhonov regularization. The inversion process started with initial guesses $E_1=1.2\times {10}^7\mathrm{Pa}$ and $E_2=3\times {10}^6\mathrm{Pa}$. Both profiles were updated simultaneously by applying the stress load $p(t)$ with $f_{\max }=9$ Hz. Then, using the updated profiles as new initial guesses, subsequent inversion was performed only for $E_2(x)$ while $E_1(x)$ was fixed and the higher frequency load ($f_{\max }=30$ Hz) was applied. By the described procedure, the smoothly-varying target profiles could be recovered fairly well as shown in Figure 7(a) and (b). Figure 7(c) represents the target, initial guess, and inverted $Q$ profiles of the medium with the smoothly-varying $E_1(x)$ and $E_2(x)$.

6.2. Layered profiles

To further demonstrate the developed inversion algorithm, a heterogeneous semi-infinite solid medium with 4 horizontal layers was considered, as depicted in Figure 8(a). The medium can also be modeled as a one-dimensional PML-truncated domain as shown in Figure 8(b). The Generalized Maxwell Body was used to describe the medium’s viscoelastic properties.

Figure 9(a) and (b) depict the target $E_1$ and $E_2$ profiles within the PML-truncated domain, respectively. Both profiles have four distinct layers with the same material interface locations. In this setting, $Q$ is different for each of the layers. Figure 9(c) presents the target $Q$ distribution along the entire domain. The density, reflection coefficient of PML, element length, surface loads, and relaxation time parameters are the same as in the case of smoothly-varying material profiles.

Figure 8. (a) Heterogeneous semi-infinite solid medium with 4 horizontal layers; (b) Schematic of a one-dimensional PML-truncated domain.

Figure 9. Target viscoelastic material profiles with 4 layers.

Figure 10. Initial guess, target, and inverted $E_1$ profiles when $E_2$ is fixed to the target.

Three cases were considered for the inversion of layered viscoelastic parameters. First, it was sought to recover the $E_1$ profile, while assuming that the $E_2$ profile was already known. Secondly, the $E_2$ profile was sought, while the distribution of $E_1$ was assumed to be known. Finally, it was sought to reconstruct both profiles simultaneously.

6.2.1. Inverting for $E_1$ when $E_2$ is known

This section presents the attempt to reconstruct $E_1(x)$ when $E_2(x)$ was already known. Equation (71(71) $$\begin{array}{c}\hfill {\mathit{\eta }}_k=\frac{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}-{\left({\mathbf{W}}_{\mathit{E}}\right)}_{k-1}{\overline{\mathbf{E}}}_{k-1}||}{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}||},\end{array}$$(71) ) was used to update $E_1$. Figure 10 shows the reconstructed $E_1$ profile using TV regularization, when $E_2$ remained the same as the target. The true profile is recovered fairly well. The inversion started with a homogeneous initial guess of $E_1={10}^7\mathrm{N}/{\mathrm{m}}^2$ and the regularization factor continuation scheme was utilized to assist in convergence. From this observation, it is noted that the elastic part ($E_1(x)$) of the Generalized Maxwell Body can be reconstructed quite well when $E_2(x)$ is known.

6.2.2. Inverting for $E_2$ when $E_1$ is fixed

In this section, the inversion process was implemented for recovering $E_2(x)$ when $E_1(x)$ was fixed to the target. Equation (71(71) $$\begin{array}{c}\hfill {\mathit{\eta }}_k=\frac{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}-{\left({\mathbf{W}}_{\mathit{E}}\right)}_{k-1}{\overline{\mathbf{E}}}_{k-1}||}{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}||},\end{array}$$(71) ) was used to update $E_2$ while the $E_1$ profile was held constant. A homogeneous initial guess of $E_2=5\times {10}^6\mathrm{N}/{\mathrm{m}}^2$ was utilized and the target profile was recovered using TN regularization. Figure 11(a) shows that the $E_2$ profile was reconstructed fairly well when $E_1$ was known. Figure 11(b) shows the $Q$ profile recovered using Equation (3(3) $$\begin{array}{cc}\hfill \frac{E_2}{E_1+2E_2}& =\mathit{\alpha }{Q}^{\mathit{\beta }},\hfill \end{array}$$(3) ) as well as the reconstructed $E_2$ profile. From these results, it is apparent that the damping-related modulus ($E_2(x)$) of the viscoelastic medium can be recovered quite well provided that the elastic property ($E_1(x)$) is known in advance.

Figure 11. (a) Initial guess, target, and inverted $E_2$ profiles when $E_1$ is fixed to the target; (b) Reconstructed $Q$ factor profile when inverting for $E_2$ alone.

6.2.3. Simultaneous inversion for $E_1$ and $E_2$

Finally, simultaneous recovery of both $E_1(x)$ and $E_2(x)$ was attempted. This problem entailed reconstruction of both elastic and damping parameters of the Generalized Maxwell Body, given the same surficial measurements used in the case of the single variable inversion. To this end, Equation (71(71) $$\begin{array}{c}\hfill {\mathit{\eta }}_k=\frac{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_k-{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}-{\left({\mathbf{W}}_{\mathit{E}}\right)}_{k-1}{\overline{\mathbf{E}}}_{k-1}||}{||{\left({\mathbf{g}}_{\mathit{E}}\right)}_{k-1}||},\end{array}$$(71) ) was used to update $E_1$ and $E_2$ profiles at the same time. One difficulty in this problem is that solution multiplicity becomes exacerbated due to the presence of dual inversion variables. To avoid falling to undesirable local minima, initial guesses and regularization factors should be carefully chosen for both inversion parameters, but there is, in general, no guarantee of success.

First, two-parameter inversion was considered for a constant target profile $E_1(x)={10}^7\mathrm{N}/{\mathrm{m}}^2$ and the 4-layer target profile for $E_2(x)$ depicted in Figure 4(b). Figure 12 shows the target, initial guess, and inverted profiles for $E_1(x)$ and $E_2(x)$, respectively. As shown in the figure, the reconstruction of $E_1(x)$ was reasonably accurate, while the reconstruction of $E_2(x)$ provided moderately reasonable estimates of the thickness, location, and material property values of each layer. Homogeneous initial guesses of $E_1=7.5\times {10}^6\mathrm{N}/{\mathrm{m}}^2$ and $E_2=5\times {10}^6\mathrm{N}/{\mathrm{m}}^2$ were used and then $E_1$ and $E_2$ were updated simultaneously at every inversion iteration with the aid of the regularization factor continuation scheme.[24, 25] Figure 13 shows the recovered $Q$ profile obtained from the inverted $E_1$ and $E_2$ results. The agreement between the reconstructed properties and the target properties was marginally successful.

Figure 12. Initial guess, target, and inverted profiles for $E_1$ and $E_2$, respectively.

Figure 13. Target and inverted $Q$ factor profiles; $E_1$ is constant and $E_2$ has 4 layers.

Figure 14. Initial guess, target, and inverted profiles for $E_1$ and $E_2$; both profiles have 4 layers.

Figure 15. Target and inverted $Q$ factor profiles; both $E_1$ and $E_2$ have 4 layers.

Next, two-parameter inversion was attempted when both profiles were layered. In particular, it was attempted to recover the layered profiles $E_1(x)$ and $E_2(x)$ depicted in Figure 4. In this case, simultaneously inverting for both properties at every inversion iteration did not lead to satisfactory reconstruction of the target profiles, regardless of the choice of initial guesses, step lengths, and regularization factors. To attempt to overcome this difficulty, a staggered inversion method was used, where single variable inversion for $E_1$ was followed by simultaneous inversion with a source frequency continuation scheme.[25] In this example, the details are as follows:

(1)	Constant and linear initial guesses were used for $E_1(x)$ and $E_2(x)$, respectively. The initial profiles were:(75) $$\begin{array}{cc}\hfill E_1(x)& ={10}^7\mathrm{N}/{\mathrm{m}}^2\hfill \end{array}$$(75) (76) $$\begin{array}{cc}\hfill E_2(x)& =(-0.042x+6)\times {10}^7\mathrm{N}/{\mathrm{m}}^2\hfill \end{array}$$(76)
(2)	Then, the initial guess of $E_2(x)$ was fixed, and inversion was accomplished only for $E_1(x)$ using the Gaussian pulse load with $f_{\max }=9$ Hz shown in Figure 5. With the updated $E_1(x)$ and the initial guess of $E_2(x)$, simultaneous inversion was attempted using a higher frequency load. In this example, a Gaussian load with $f_{\max }=30$ Hz was used for the simultaneous inversion.

Figure 14 shows the results of the two-parameter inversion implemented by the staggered inversion procedure. As can be seen in the figure, the inverted $E_1(x)$ captures the target pretty well, while the accuracy of $E_2(x)$ was limited. The inaccurate recovery pattern near the interfaces in $E_2(x)$ does not match observations from the case of the single variable inversion for $E_2(x)$ (Figure 11(a)). However, the variation of the layers can be roughly observed, with the value of $E_2$ closely matching the target around the middle of each layer.

Figure 16. Reconstructed profiles obtained from Gauss–Newton and Fletcher–Reeves optimization algorithms; $E_1$ is constant and $E_2$ has 4 layers.

Figure 17. Response misfit versus iteration numbers in the simultaneous inversion; $E_1$ is constant and $E_2$ has 4 layers.

Figure 15 shows the $Q$ profile constructed from the inverted $E_1(x)$ and $E_2(x)$ profiles. Regarding the $Q$ profile, there are some errors, but the recovered Q profile characterizes the distribution of the domain’s attenuation reasonably well.

6.2.4. Gauss–Newton versus Fletcher–Reeves optimization method

To investigate the efficiency of the developed Gauss–Newton (GN) approach for full-waveform inversion, results from the GN method were compared with those from the Fletcher–Reeves (FR) conjugate gradient optimization method. Figure 16 compares inverted $E_1(x)$ and $E_2(x)$ profiles obtained from the GN and FR methods, respectively, when $E_1$ was held constant and $E_2$ had 4 layers. Both methods resulted in profiles close to the target. However, it is noted that the accuracy of solutions was improved for the GN profile, especially when $x$ was greater than 50 m. Additionally, it took considerably fewer iterations in the case of GN than in the case of FR. Figure 17 exhibits response misfit versus iteration numbers in the simultaneous inversion. The GN method required 130 iterations, whereas it took 820 iterations when the FR method was used for the inversion.

Figure 18. Reconstructed profiles obtained from Gauss–Newton and Fletcher–Reeves optimization algorithms; both $E_1$ and $E_2$ have 4 layers.

Figure 19. Response misfit versus iteration numbers in the simultaneous inversion; both $E_1$ and $E_2$ have 4 layers.

In Figure 18, GN and FR profiles are compared when both $E_1$ and $E_2$ have 4 layers. The accuracy of solutions was improved by the developed GN approach compared to the FR, with the GN profile sharper than the FR profile at interfaces. The developed approach also attained substantial improvement on the convergence rate of solutions for both load frequencies $f_{\max }=9\text{and}30\text{Hz}$, as seen in Figure 19. Using $f_{\max }=9$ Hz, the GN and FR methods required 110 and 7780 iterations, respectively, to approximate the target. These iteration numbers were 350 and 2860 for GN and FR methods, respectively, when $f_{\max }=30$ Hz.

For the case of the constant $E_1$ profile, the CPU times to perform GN and FR inversions were about 430 and 3100 s, respectively. For the case where both $E_1$ and $E_2$ were layered, the CPU times for GN and FR inversions were about 770 and 5300 s, respectively. In both cases, the CPU time of GN inversion was less than 15% that of FR inversion.

7. Conclusions

This study presents the development of a Gauss–Newton full-waveform inversion in the time-domain for the reconstruction of elastic and attenuating properties of semi-infinite solid media truncated by PML, using solely the surficial response of the domain to dynamic excitation. The nonlinear optimization problem was resolved using an iterative procedure in the reduced space of the control variables, which were viscoelastic parameters of the PML-truncated domains. The search direction of $E_1(x)$ and $E_2(x)$ were determined by a Gauss–Newton–Krylov method on a reduced KKT system obtained by second-order Fréchet derivatives of the Lagrangian functional.

The numerical results indicate that, while the elastic properties are well captured, reconstruction of the attenuation profiles is more challenging. Still, satisfactory results were obtained when a staggered scheme was employed. Approximate solutions of the Gauss–Newton inversion exhibited faster convergence behavior and better robustness, while requiring fewer iterations and less CPU time to reconstruct the profiles compared to the Fletcher-Reeves optimization method.

Additional information

Funding

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea [2014R1A1A1006419].

Notes

No potential conflict of interest was reported by the authors.

Appendix 1

Second-order Fréchet derivatives of the Lagrangian functional

The second-order Fréchet derivatives of $\mathcal{L}$ in the KKT system (Equation (39(39) $$\begin{array}{cc}\hfill \mathcal{B}\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{uu}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{u\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{\mathit{\sigma }u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{\mathit{\sigma }\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}\overline{u}(x,t)\mathit{\delta }(x)\\ 0\end{array}\right],\hfill \end{array}$$(39) )) can be obtained by taking the second variation of the Lagrangian with respect to the state, adjoint, and control variables. For example, ${\mathrm{\nabla }}_{uu}^2\mathcal{L}$ can be obtained by the following variational calculus; given arbitrary functions $\widehat{u}(x,t)$ and $\overline{u}(x,t)$ with $\mathit{ϵ}$, $\mathit{\tau }\in \mathbb{R}$, the perturbation of the Lagrangian (14(14) $$\begin{array}{cc}& \mathcal{L}(u,\mathit{\sigma },{\mathit{\lambda }}_u,{\mathit{\lambda }}_{\mathit{\sigma }},E_1,E_2),\hfill \\ \hfill & =\frac{1}{2}{\int }_0^T{\left[u(0,t)-u_{\mathrm{m}}(0,t)\right]}^2\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2)\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_u\left[\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }u}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}\right]\mathrm{d}t\mathrm{d}x\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_{\mathit{\sigma }}\left[{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+h\tilde{\mathit{\sigma }}-E_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}\right]\mathrm{d}t\mathrm{d}x.\hfill \\ \hfill \end{array}$$(14) ) can be written as:(A1) $$\begin{array}{cc}& \mathcal{L}(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u},\mathit{\sigma },{\mathit{\lambda }}_u,{\mathit{\lambda }}_{\mathit{\sigma }},E_1,E_2)\hfill \\ \hfill & =\frac{1}{2}{\int }_0^T{\int }_0^{L_t}{\left[u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}-u_m\right]}^2\mathit{\delta }(x)\mathrm{d}x\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2)\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_u\left[\frac{{\mathit{\partial }}^2\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}\right]\mathrm{d}t\mathrm{d}x\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_{\mathit{\sigma }}[{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+\overline{\mathit{\sigma }}-E_1\frac{\mathit{\partial }\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }x}\hfill \\ \hfill & -{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^2\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^3\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }x\mathit{\partial }{t}^2}]\mathrm{d}t\mathrm{d}x.\hfill \end{array}$$(A1)

Differentiating Equation (A1(A1) $$\begin{array}{cc}& \mathcal{L}(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u},\mathit{\sigma },{\mathit{\lambda }}_u,{\mathit{\lambda }}_{\mathit{\sigma }},E_1,E_2)\hfill \\ \hfill & =\frac{1}{2}{\int }_0^T{\int }_0^{L_t}{\left[u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}-u_m\right]}^2\mathit{\delta }(x)\mathrm{d}x\mathrm{d}t+{\mathcal{R}}_{E_1}(E_1)+{\mathcal{R}}_{E_2}(E_2)\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_u\left[\frac{{\mathit{\partial }}^2\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }t}-\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }x}\right]\mathrm{d}t\mathrm{d}x\hfill \\ \hfill & +{\int }_0^{L_t}{\int }_0^T{\mathit{\lambda }}_{\mathit{\sigma }}[{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\mathit{\sigma }}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\mathit{\sigma }}{\mathit{\partial }t}+{\mathit{\kappa }}_2\mathit{\sigma }+\overline{\mathit{\sigma }}-E_1\frac{\mathit{\partial }\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }x}\hfill \\ \hfill & -{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^2\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^3\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u}\right)}{\mathit{\partial }x\mathit{\partial }{t}^2}]\mathrm{d}t\mathrm{d}x.\hfill \end{array}$$(A1) ) with respect to $\mathit{ϵ}$ and $\mathit{\tau }$, respectively, and then evaluating the resulting equation at $\mathit{ϵ}=\mathit{\tau }=0$ yields:(A2) $$\begin{array}{c}\hfill \frac{{\mathit{\partial }}^2}{\mathit{\partial }\mathit{ϵ}\mathit{\partial }\mathit{\tau }}\mathcal{L}\left(u+\mathit{ϵ}\widehat{u}+\mathit{\tau }\overline{u},\mathit{\sigma },{\mathit{\lambda }}_v,{\mathit{\lambda }}_{\mathit{\sigma }},E_1,E_2\right){|}_{\mathit{ϵ}=\mathit{\tau }=0}={\int }_0^T{\int }_0^{L_t}\widehat{u}\overline{u}\mathit{\delta }(x)\mathrm{d}x\mathrm{d}t.\end{array}$$(A2)

Therefore,(A3) $$\begin{array}{c}\hfill \left({\mathrm{\nabla }}_{uu}^2\mathcal{L}\right)\overline{u}=\overline{u}(x,t)\mathit{\delta }(x),\mathrm{f}or0<x<L_t,0\le t\le T.\end{array}$$(A3)

The rest of the Fréchet derivatives in Equation (39(39) $$\begin{array}{cc}\hfill \mathcal{B}\overline{\mathit{u}}& :=\left[\begin{array}{ccc}{\mathrm{\nabla }}_{uu}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{u\mathit{\sigma }}^2\mathcal{L}\hfill \\ {\mathrm{\nabla }}_{\mathit{\sigma }u}^2\mathcal{L}\hfill & \hfill & {\mathrm{\nabla }}_{\mathit{\sigma }\mathit{\sigma }}^2\mathcal{L}\hfill \end{array}\right]\left[\begin{array}{c}\overline{u}\hfill \\ \overline{\mathit{\sigma }}\hfill \end{array}\right]=\left[\begin{array}{c}\overline{u}(x,t)\mathit{\delta }(x)\\ 0\end{array}\right],\hfill \end{array}$$(39) ) can be derived similarly. There result the following expressions:

Table A1. Second-order Fréchet derivatives of $\mathcal{L}$.

$\left({\mathrm{\nabla }}_{uu}^2\mathcal{L}\right)\overline{u}=\overline{u}(x,t)\mathit{\delta }(x)$,	$\left({\mathrm{\nabla }}_{u\mathit{\sigma }}^2\mathcal{L}\right)\overline{\mathit{\sigma }}=0$,
${\displaystyle \left({\mathrm{\nabla }}_{u{E}_1}^2\mathcal{L}\right){\overline{E}}_1=\frac{\mathit{\partial }}{\mathit{\partial }x}\left({\mathit{\lambda }}_{\mathit{\sigma }}{\overline{E}}_1\right)-({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }}{\mathit{\partial }x}\left({\overline{E}}_1\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\right)+{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }}{\mathit{\partial }x}\left({\overline{E}}_1\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\right)}$,
${\displaystyle \left({\mathrm{\nabla }}_{u{E}_2}^2\mathcal{L}\right){\overline{E}}_2=-({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }}{\mathit{\partial }x}\left({\overline{E}}_2\frac{\mathit{\partial }{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }t}\right)+2{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }}{\mathit{\partial }x}\left({\overline{E}}_2\frac{{\mathit{\partial }}^2{\mathit{\lambda }}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\right)}$,
${\displaystyle \left({\mathrm{\nabla }}_{u{\mathit{\lambda }}_u}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_u=\mathit{\rho }\frac{{\mathit{\partial }}^2{\overline{\mathit{\lambda }}}_u}{\mathit{\partial }{t}^2}-\mathit{\rho }h\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_u}{\mathit{\partial }t}}$,	${\displaystyle \left({\mathrm{\nabla }}_{u{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}=\frac{\mathit{\partial }}{\mathit{\partial }x}\left(E_1{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}-{\mathit{\kappa }}_3\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^2{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}\right)}$,
$\left({\mathrm{\nabla }}_{\mathit{\sigma }u}^2\mathcal{L}\right)\overline{u}=0$,	$\left({\mathrm{\nabla }}_{\mathit{\sigma }\mathit{\sigma }}^2\mathcal{L}\right)\overline{\mathit{\sigma }}=0$,
$\left({\mathrm{\nabla }}_{\mathit{\sigma }{E}_1}^2\mathcal{L}\right){\overline{E}}_1=0$,	$\left({\mathrm{\nabla }}_{\mathit{\sigma }{E}_2}^2\mathcal{L}\right){\overline{E}}_2=0$,
${\displaystyle \left({\mathrm{\nabla }}_{\mathit{\sigma }{\mathit{\lambda }}_u}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_u=\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_u}{\mathit{\partial }x}}$,	${\displaystyle \left({\mathrm{\nabla }}_{\mathit{\sigma }{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}={\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }{t}^2}-{\mathit{\kappa }}_1\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}-h{\overline{\tilde{\mathit{\lambda }}}}_{\mathit{\sigma }}}$,
${\displaystyle \left({\mathrm{\nabla }}_{E_{1}u}^2\mathcal{L}\right)\overline{u}=-{\mathit{\lambda }}_{\mathit{\sigma }}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2},\left({\mathrm{\nabla }}_{E_1\mathit{\sigma }}^2\mathcal{L}\right)\overline{\mathit{\sigma }}=0}$,
$\left({\mathrm{\nabla }}_{E_1{E}_1}^2\mathcal{L}\right){\overline{E}}_1=\left\{\begin{array}{cc}{\displaystyle -R_{E_1}\frac{{\mathrm{d}}^2{\overline{E}}_1}{\mathrm{d}{x}^2},}\hfill & \mathrm{f}orTNregularization\hfill \\ {\displaystyle -e{R}_{E_1}\frac{\mathrm{d}}{\mathrm{d}x}\left\{{\left({\left(\frac{\mathrm{d}{E}_1}{\mathrm{d}x}\right)}^2+e\right)}^{-\frac{3}{2}}\frac{\mathrm{d}{\overline{E}}_1}{\mathrm{d}x}\right\},}\hfill & \mathrm{f}orTVregularization\hfill \end{array}\right.$
$\left({\mathrm{\nabla }}_{E_1{E}_2}^2\mathcal{L}\right){\overline{E}}_2=0$,
$\left({\mathrm{\nabla }}_{E_1{\mathit{\lambda }}_u}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_u=0$,
${\displaystyle \left({\mathrm{\nabla }}_{E_1{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}=-{\int }_0^T{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}t+{\int }_0^T({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}t+{\int }_0^T{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\mathrm{d}t}$,
(Continued)

Table A2. (Continued)

${\displaystyle \left({\mathrm{\nabla }}_{E_{2}u}^2\mathcal{L}\right)\overline{u}=-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }{x}^2}-2{\mathit{\tau }}_1{\mathit{\tau }}_2{\mathit{\lambda }}_{\mathit{\sigma }}\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2},\left({\mathrm{\nabla }}_{E_2\mathit{\sigma }}^2\mathcal{L}\right)\overline{\mathit{\sigma }}=0}$,
$\left({\mathrm{\nabla }}_{E_2{E}_1}^2\mathcal{L}\right){\overline{E}}_1=0$,
$\left({\mathrm{\nabla }}_{E_2{E}_2}^2\mathcal{L}\right){\overline{E}}_2=\left\{\begin{array}{cc}{\displaystyle -R_{E_2}\frac{{\mathrm{d}}^2{\overline{E}}_2}{\mathrm{d}{x}^2},}\hfill & \mathrm{f}orTNregularization\hfill \\ {\displaystyle -e{R}_{E_2}\frac{\mathrm{d}}{\mathrm{d}x}\left\{{\left({\left(\frac{\mathrm{d}{E}_2}{\mathrm{d}x}\right)}^2+e\right)}^{-\frac{3}{2}}\frac{\mathrm{d}{\overline{E}}_2}{\mathrm{d}x}\right\},}\hfill & \mathrm{f}orTVregularization\hfill \end{array}\right.$
$\left({\mathrm{\nabla }}_{E_2{\mathit{\lambda }}_u}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_u=0$,
${\displaystyle \left({\mathrm{\nabla }}_{E_2{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}={\int }_0^T({\mathit{\tau }}_1+{\mathit{\tau }}_2)\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{\mathit{\partial }u}{\mathit{\partial }x}\mathrm{d}t+2{\int }_0^T{\mathit{\tau }}_1{\mathit{\tau }}_2\frac{\mathit{\partial }{\overline{\mathit{\lambda }}}_{\mathit{\sigma }}}{\mathit{\partial }t}\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}\mathrm{d}t}$,
${\displaystyle \left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{u}u}^2\mathcal{L}\right)\overline{u}=\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }{t}^2}+\mathit{\rho }h\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }t}}$,	${\displaystyle \left({\mathrm{\nabla }}_{{\mathit{\lambda }}_u\mathit{\sigma }}^2\mathcal{L}\right)\overline{\mathit{\sigma }}=-\frac{\mathit{\partial }\overline{\mathit{\sigma }}}{\mathit{\partial }x}}$,
$\left({\mathrm{\nabla }}_{{\mathit{\lambda }}_u{E}_1}^2\mathcal{L}\right){\overline{E}}_1=0$,	$\left({\mathrm{\nabla }}_{{\mathit{\lambda }}_u{E}_2}^2\mathcal{L}\right){\overline{E}}_2=0$,
$\left({\mathrm{\nabla }}_{{\mathit{\lambda }}_u{\mathit{\lambda }}_u}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_u=0$,	$\left({\mathrm{\nabla }}_{{\mathit{\lambda }}_v{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}=0$,
${\displaystyle \left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}u}^2\mathcal{L}\right)\overline{u}=-E_1\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}-{\mathit{\kappa }}_3\frac{{\mathit{\partial }}^2\overline{u}}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\kappa }}_4\frac{{\mathit{\partial }}^3\overline{u}}{\mathit{\partial }x\mathit{\partial }{t}^2}}$,
${\displaystyle \left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}\mathit{\sigma }}^2\mathcal{L}\right)\overline{\mathit{\sigma }}={\mathit{\tau }}_1{\mathit{\tau }}_2\frac{{\mathit{\partial }}^2\overline{\mathit{\sigma }}}{\mathit{\partial }{t}^2}+{\mathit{\kappa }}_1\frac{\mathit{\partial }\overline{\mathit{\sigma }}}{\mathit{\partial }t}+{\mathit{\kappa }}_2\overline{\mathit{\sigma }}+h\overline{\tilde{\mathit{\sigma }}}}$,
${\displaystyle \left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}{E}_1}^2\mathcal{L}\right){\overline{E}}_1={\overline{E}}_1\frac{\mathit{\partial }u}{\mathit{\partial }x}-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\overline{E}}_1\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-{\mathit{\tau }}_1{\mathit{\tau }}_2{\overline{E}}_1\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}}$,
${\displaystyle \left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}{E}_2}^2\mathcal{L}\right){\overline{E}}_2=-({\mathit{\tau }}_1+{\mathit{\tau }}_2){\overline{E}}_2\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }x\mathit{\partial }t}-2{\mathit{\tau }}_1{\mathit{\tau }}_2{\overline{E}}_2\frac{{\mathit{\partial }}^{3}u}{\mathit{\partial }x\mathit{\partial }{t}^2}}$,
$\left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}{\mathit{\lambda }}_u}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_u=0$,	$\left({\mathrm{\nabla }}_{{\mathit{\lambda }}_{\mathit{\sigma }}{\mathit{\lambda }}_{\mathit{\sigma }}}^2\mathcal{L}\right){\overline{\mathit{\lambda }}}_{\mathit{\sigma }}=0$.