A Survey on Inverse Problems for Applied Sciences

Abstract

The aim of this paper is to introduce inversion-based engineering applications and to investigate some of the important ones from mathematical point of view. To do this we employ acoustic, electromagnetic, and elastic waves for presenting different types of inverse problems. More specifically, we first study location, shape, and boundary parameter reconstruction algorithms for the inaccessible targets in acoustics. The inverse problems for the time-dependent differential equations of isotropic and anisotropic elasticity are reviewed in the following section of the paper. These problems were the objects of the study by many authors in the last several decades. The physical interpretations for almost all of these problems are given, and the geophysical applications for some of them are described. In our last section, an introduction with many links into the literature is given for modern algorithms which combine techniques from classical inverse problems with stochastic tools into ensemble methods both for data assimilation as well as for forecasting.

1. Introduction

Finding causes from the knowledge of their effects, which actually constitutes the idea of solving inverse problems, is necessary for the identifications of practical applications in different critical areas such as mine detection, medical imaging, remote sensing, nondestructive testing, and geophysical explorations. It is also a crucial ingredient of forecasting in basic areas such as weather prediction or projections for climate change.

However, small changes in the effects might result in large differences in the causes, or the same effect might be obtained from more than one cause. Therefore, it is difficult or sometimes impossible to find actual reason uniquely by observing only effects. This ill-posed fundamental characteristic of inverse problems influences the complexity and features of their solutions dramatically [1–6].

It is a generally accepted view that the first mathematical investigation of inverse problems is the study of Abel’s on a mechanical problem for finding the curve of an unknown path in 1826. On the other hand, the invention of radar and sonar during the Second World War inspired researchers to focus especially on inverse scattering problems whose aims are not only to determine locations of the targets from the transmitter/receiver antennas but also to construct their detailed images. This motivation induced the progress in developing new reconstruction methods and their extensions to other research branches such as nondestructive testing, biomedical imaging, seismology, and atmospheric profile inversions. This agenda has led to a large number of methods and tools since 1980, with a continuing flow of new ideas into the field.

Several inversion methods in scattering and tomography have been suggested in the 1980s and 1990s, with qualitative methods arising since around 1996. Since 2000, an increased interest in the treatment of sparsity has led to new regularization tools and new approaches. Also, the need to carry out inversion as a part of dynamical systems and forecasting, also called data assimilation, is leading to an increasing synergy of filtering methods from stochastics and inversion tools from regularization theory. Since around 2000, multiphysics-based methods have become very popular, leading to increased insight for example into the human body or the atmosphere. Nowadays the advent of powerful computers and high technologies made it possible to evaluate and process large volume of data for finding sufficiently accurate solutions of practical inverse problems. For a more detailed information into the different areas of inverse problems we refer to [7–20].

The plan of the paper is as follows. A short summary to the inverse problems in the historical perspective has just been reported in the introduction section. The following two sections are devoted for the presentation of the mathematical basics used in the solution of well-known inverse problems for acoustic, electromagnetic, and elastic waves. Clearly, we cannot aim to cover their entire literature with this survey but describe some interesting and important lines of development. Afterwards, some of important inversion-based practical applications are introduced, and as a selected topic, inverse problems for the neural field equation is described in details. In the final section, conclusions and concluding remarks are given.

2. Inverse Problems for Time-Harmonic Acoustic and Electromagnetic Waves

In this section, we focus on the mathematical investigation of inverse problems whose aims are to reconstruct geometrical and/or physical properties of penetrable/impenetrable objects from the knowledge of scattered acoustic or electromagnetic waves at certain measurement points. Physically speaking, identifications of the unknown and inaccessible objects with certain waves are the main problems of the nondestructive testing, radar/sonar applications, tumor detection, and so forth.

Let us assume a time-harmonic plane wave 𝑢^𝑖=𝑒^{𝑖(𝑘𝑥⋅𝑑−𝜔𝑡)}, which was insonificated by a sufficiently far acoustic source operating at a constant frequency 𝜔. The acoustic wave propagates in the direction 𝑑, with a speed of sound 𝑐, in a homogeneous medium having a wave number 𝑘=𝜔/𝑐. In this background medium, we consider a bounded penetrable/impenetrable scatterer with support given by a domain 𝐷∈ℝ^𝑚(𝑚=2,3) and define 𝜈 as the unit outward normal to the boundary 𝜕𝐷.

The scattered field 𝑢^𝑠, which appears from the interaction of the incident wave with the obstacle, satisfies the Sommerfeld radiation condition at infinity lim_𝑟→∞𝑟^(𝑚−1)/2𝜕𝑢^𝑠𝜕𝑟−𝑖𝑘𝑢^𝑠=0,𝑟=|𝑥|, (1) and ensures that the scattered wave has the form of an outgoing wave. Furthermore, the so-called far-field pattern 𝑢_∞ can be obtained from the asymptotic behavior of the scattered wave such that 𝑢^𝑠(𝑥)=𝑒^𝑖𝑘𝑟𝑟^(𝑚−1)/2𝑢_∞(𝑥)1+𝑂𝑟,𝑥𝑥=|𝑥|,𝑟=|𝑥|→∞. (2) The total field 𝑢 is simply the superposition of incident and scattered fields 𝑢=𝑢^𝑖+𝑢^𝑠, which satisfies the homogeneous Helmholtz equation in the exterior domain of the obstacle Δ𝑢+𝑘²𝑢=0inℝ^𝑚⧵𝐷for(𝑚=2,3) (3) and the boundary condition depending on scatterer type. Commonly used boundary conditions in the literature are given in the following.

The Dirichlet boundary condition: 𝑢=0on𝜕𝐷, (4) the Neumann boundary condition: 𝜕𝑢𝜕𝜈=0on𝜕𝐷, (5) the impedance boundary condition: 𝜂𝑢+𝑖𝑘𝜕𝑢𝜕𝜈=0on𝜕𝐷, (6) the conductive boundary condition: 𝑢₀=𝑢,𝜕𝑢−𝜕𝜈𝜕𝑢₀𝜕𝜈=𝜆𝑢on𝜕𝐷. (7)

Parameters 𝜂 and 𝜆 are the impedance and conductivity functions, respectively. In latest equation, 𝑢₀ is the total field inside of the scatterer, and this boundary condition can be considered as a more general form of the transmission conditions.

On the other hand, in order to model acoustic scattering by an inhomogeneous medium, it is sufficient to replace (3)–(7) with Δ𝑢+𝑘²𝑛(𝑥)𝑢=0inℝ^𝑚for(𝑚=2,3), (8) where 𝑛(𝑥)=𝑐²/𝑐²_𝑖(𝑥)+𝑖𝜎_𝑖(𝑥) is the refractive index, 𝑐 is the sound speed in the homogeneous medium, 𝑐_𝑖(𝑥) is the speed of sound in the inhomogeneous medium, and 𝜎_𝑖(𝑥)≥0 models the absorption.

Analogically one can also consider a scattering problem for electromagnetic waves assuming that the electric and magnetic field components of the incident plane wave are given by 𝐸^𝑖(𝑥,𝑡)=𝑖𝑘𝑑×𝑝×𝑑𝑒^{𝑖(𝑘𝑥⋅𝑑−𝜔𝑡)},𝐻^𝑖(𝑥,𝑡)=𝑖𝑘𝑑×𝑝𝑒^{𝑖(𝑘𝑥⋅𝑑−𝜔𝑡)}, (9) where 𝑘=𝜔√𝜖𝜇 is the wave number and 𝜖 and 𝜇 are the electric permittivity and magnetic permeability of the host medium, respectively. Here, 𝜔 represents the frequency of the wave source, 𝑑 is the direction, and 𝑝 is the polarization of the electromagnetic wave. Then, we consider that the wave interacts with a 3D obstacle 𝐷∈ℝ³ in the homogeneous medium and that scattered electromagnetic fields (𝐸^𝑠,𝐻^𝑠) occur similar to the previous case. However, in electromagnetics, the scattered field has to satisfy the Silver-Müller radiation condition lim_𝑟→∞𝐻^𝑠×𝑥−𝑟𝐸^𝑠=0. (10) The far field pattern of the corresponding scattered field can be obtained by 𝐸^𝑠(𝑥)=𝑒^𝑖𝑘𝑟𝑟𝐸_∞(𝑥)1+𝑂𝑟,𝑥𝑥=|𝑥|,𝑟=|𝑥|⟶∞. (11) Furthermore, the total electric field 𝐸=𝐸^𝑖+𝐸^𝑠 and the total magnetic field 𝐻=𝐻^𝑖+𝐻^𝑠 in the medium satisfy the Maxwell equations curl𝐸−𝑖𝑘𝐻=0,curl𝐻+𝑖𝑘𝐸=0inℝ³⧵𝐷 (12) and the chosen boundary conditions on the surface of the scatterer. The boundary conditions which were presented for acoustic waves are now given for electromagnetic waves in (13)–(16).

The Dirichlet-type perfect conductor boundary condition: 𝜈×𝐸=0on𝜕𝐷, (13) the Neumann-type boundary condition: 𝜈×𝐻=0on𝜕𝐷, (14) the impedance boundary condition: 𝜂𝜈×𝐸+((𝑖𝑘𝜈×𝜈×𝐻))=0on𝜕𝐷, (15) the conductive boundary condition: 𝜈×𝐸₀()=𝜈×𝐸,𝜈×𝜈×𝐻−𝜉𝜈×𝜈×𝐻₀=𝜆(𝜈×𝐸)on𝜕𝐷. (16) The parameter 𝜉 in the last equation is the ratio of the wave numbers 𝜉=𝑘₀/𝑘, where 𝑘₀ is the wave number of the object's interior domain.

For the case of modeling scattering of electromagnetic waves by an inhomogeneous medium, (12)–(16) are replaced by curl𝐸−𝑖𝑘𝐻=0,curl𝐻+𝑖𝑘𝑛(𝑥)𝐸=0inℝ³, (17) where 𝑛(𝑥)=1/𝜖(𝜖_𝑖+𝑖𝜎_𝑖(𝑥)/𝜔) is the refractive index, 𝜖 is the permittivity of the homogeneous medium, 𝜖_𝑖 is the permittivity of the inhomogeneous medium, 𝜎_𝑖(𝑥) is the conductivity, and 𝜔 is the frequency of the wave.

We recommend to the readers the books [9, 21–23] for the detailed investigations of the electromagnetic and acoustic waves with the same notation.

2.1. Boundary Reconstruction Problems of Acoustic Waves

In this section, our main concern is to present mathematical treatments of some available methods for reconstructing shapes of obstacles from the knowledge of scattered field. Solution of this type of problems with acoustic waves are important from practical point of view since the real life experiments with acoustic waves are easier, safer, and lower in cost.

To provide a simpler mathematical presentation, we consider a two dimensional simply connected object 𝐷∈ℝ² having the Dirichlet condition on its boundary Γ. It is further assumed that the incident acoustic wave 𝑢^𝑖, which interacts with the object, and scattered near/far fields are given. This inverse problem is nonlinear due to the mathematical relation between the scattered wave and the shape of the cylinder, and it is ill-posed since the determination of Γ does not depend continuously on the scattered field [18]. For the solution of this type of inverse problems, there are a variety of methods of which we present a short list in the following without claiming to cover the entire literature.

Iterative Methods(i)Landweber iterations [24, 25].(ii)Regularized Newton [22, 26–50].(iii)Newton-Kantorovich [51, 52].

Decomposition Methods(i)Colton-Monk [53–55].(ii)Kirsch-Kress [56–58].(iii)Angell-Kleinmann-Roach [59–61].(iv)Hybrid method [62–67].(v)Potthast's point source method [23, 68, 69].

Probe and Sampling Methods(i)Linear sampling [70–74].(ii)Factorization [75, 76].(iii)Singular sources [19, 23, 77].(iv)Probe method [78, 79].(v)Enclosure method [80, 81].(vi)No-response test [82–84].

In the first group of methods [24, 52], the inverse obstacle problem is considered as an ill-posed nonlinear operator equation, 𝐴Γ₀=𝑢_∞, (18) for an initial boundary Γ₀. Often, it is assumed to be a star-like parametrization Γ₀={𝑧(𝑡)∶=𝑟(𝑡)(cos𝑡,sin𝑡)∶𝑡∈[0,2𝜋)}. Here, the operator 𝐴∶Γ→𝑢_∞ is actually defined to map the boundary Γ of the scatterer onto the far field 𝑢_∞. Newton's method to solve a nonlinear equation employs successive linearization; that is, we replace (18) by an equation for the update ℎ_𝑛𝐴𝑧_𝑛+𝐴^′𝑧_𝑛ℎ_𝑛=𝑢_∞, (19) where 𝑧_𝑛+1=𝑧_𝑛+ℎ_𝑛,𝑛=0,1,2,…, (20) leading to a sequence of linear equations (19) and updates (20) with starting value 𝑧₀∶=Γ₀. Equation (19) is ill-posed and needs to be regularized. As a simple option, Tikhonov’s regularization with regularization parameter 𝛼>0 can be applied, leading to 𝛼ℎ_𝑛+𝐴𝑧_𝑛^′∗𝐴^′𝑧_𝑛ℎ_𝑛=𝐴𝑧_𝑛^′∗𝑢_∞−𝐴𝑧_𝑛 (21) in each iteration step 𝑛=0,1,2,…. Iterative methods for ill-posed equations need to be stopped after a finite number of iterations to keep the ill-posedness under control; compare, for example, [23, 85].

The main idea of the decomposition methods is to split full nonlinear shape reconstruction problem given by (18) into a linear ill-posed equation which is solved first and a nonlinear well-posed equation to be solved in a second step.

In order to discuss the Kirsch-Kress method [56–58] via potential approach let us employ an initial boundary Γ₀ as an approximation of the actual boundary Γ. We assume that Γ₀ is in the interior of the true scatterer.

In this case, the approximate total field 𝑢(𝑥) can be approximated by the sum of the incident wave 𝑢^𝑖 and a single-layer potential, 𝑢(𝑥)=𝑢^𝑖(𝑥)+_Γ0Φ𝑥,𝑦𝜑𝑦𝑑𝑠𝑦,𝑥∈ℝ²⧵𝐷, (22) where Φ(𝑥,𝑦) is the fundamental solution to the Helmholtz equation in two dimensions and 𝜑 is a continuous density (source) function defined over Γ₀. The far field pattern of the single-layer potential should be measured far field pattern 𝑢^(∞), which leads to the linear equation 𝑢^(∞)(𝑥)=𝑐_Γ0𝑒^−𝑖𝜅𝑥_⋅𝑦𝜑𝑦𝑑𝑠𝑦,𝑥∈𝕊, (23) with some constant 𝑐 depending on the dimension of the space under consideration. The first step consists of solving the ill-posed linear equation (23) to calculate 𝜑.

Then, when we consider a Dirichlet boundary condition, the shape is found in a second nonlinear step as the zero curve of the total field 𝑢(𝑥). To find this zero curve, we introduce an operator 𝐺, which maps Γ to the values of the approximate total field 𝑢, on Γ such that 𝐺∶Γ↦𝑢∣_Γ. (24) Then the problem is reduced to the solution of the following optimization problem: 𝐺(Γ)=0, (25) which can be done by minimizing the defect ‖𝐺(Γ)‖_𝐿²_(Γ) in a least square sense.

Linear sampling [70–74] is based on solving the far field equation _Ω𝑢_∞()𝑔(𝑑)(𝑑)=𝑒𝑥;𝑑𝑑𝑠^𝑖𝜋/4√𝑒8𝜋𝑘^−𝑖𝑘𝑥^⋅𝑦𝑥∈𝕊. (26) Then, the density 𝑔 satisfies 𝑔⋅,𝑦_𝐿2(Ω)⟶∞,𝑦⟶𝜕𝐷, (27) that is, when the source point 𝑦 approaches the boundary Γ. This behaviour can be used to visualize the shape of the scatterer from the knowledge of the far field pattern 𝑢_∞(⋅,𝑑) for all directions 𝑑∈𝕊.

We also refer to the orthogonality sampling method, which has recently been suggested in [86]. It is particularly suited to deal with multifrequency data as is naturally obtained when acoustic pulses are used to probe an object or region in space. The method has been independently suggested by Ito et al. [87, 88] and successfully applied as a first step in a larger inversion procedure. A convergence analysis of the method in the limit of small scatterers has recently been achieved by Griesmaier [89].

2.2. Parameter on the Boundary Reconstruction Problems of Acoustic Waves

The inverse problem considered in this section is finding a continuous function which is defined on the boundary of the obstacle from the knowledge of the scattered acoustic waves for a given shape in two dimensions. We discuss the methods whose aims are to reconstruct 𝜂, in the impedance boundary condition (6), and 𝜆, in the conductive boundary condition (7). Note that conductivity function is employed to model the inhomogeneity which might exist on the boundary of an object in a more realistic way.

In [90–92], 3D obstacles with impedance boundary condition are studied for acoustic case, where in [91] electromagnetic case is included. In the same field, the papers [64, 93–96] are focused on 2D geometries for the less complexity of governing numerical experiments.

In [90], the impedance reconstruction problem is solved theoretically with a technique depending on Backus and Gilbert’s method which is applicable to linear moment problems. To this aim, approximate Green's functions are used to reduce the nonlinear problem to two linear moment problems. On the other hand, the study [92] is devoted for the reconstruction of impedance functions via the Kirsch-Kress and Colton-Monk decomposition methods. Furthermore, some interesting papers appeared on the impedance reconstructions, recently [97, 98].

The paper [93] introduced a new method for impedance reconstructions in the spirit of the Kirsch-Kress decomposition method. That is, the scattered field is represented via single-layer potential over the known boundary of the impedance cylinder Γ, instead of defining an auxiliary initially guessed curve. Then, the density function 𝜑 is solved from the following ill-posed integral equation 𝑢_∞(𝑥)=𝑒^𝑖𝜋/4√8𝜋𝑘^Γ𝑒^−𝑖𝑘𝑥^⋅𝑦𝜑𝑦𝑑𝑠𝑦,𝑥∈Ω, (28) through Tikhonov’s regularization for known far field, where Ω is a unit circle. From the knowledge of the density function now the total field and the normal derivative of the total field can be computed on the boundary of the obstacle via jump relations [22]. Finally, 𝜂 is obtained from (6) in the least squares sense. This method is also extended for the reconstructions of the conductivity functions of the obstacles in free space [99], for the obstacles buried in penetrable cylinders [100] and for a combination of a shape and conductivity function reconstruction problem [101], firstly by Yaman [6].

Moreover, [64, 95, 96] are devoted for the shape and impedance reconstructions of 2D obstacles in acoustics. To do this, the hybrid method is employed by Serranho [64, 66]. In [95], a level set algorithm is combined with boundary integral equations in acoustic case to reconstruct the shape and impedance of 2D obstacles from multi-illuminations, and in [96], it has been shown that the knowledge of the scattered fields corresponding to three incident waves can be used for the determination of the shape and the impedance via integral equation methods and conformal mapping techniques.

3. Inverse Problems for Differential Equations of Elastodynamics

3.1. Differential Equations of Elastodynamics

Let 𝑥=(𝑥₁,𝑥₂,𝑥₃)∈ℝ³ be a 3D space variable, and let 𝑡∈ℝ be a 1D time variable; let 𝐮(𝑥,𝑡)=(𝑢₁(𝑥,𝑡),𝑢₂(𝑥,𝑡),𝑢₃(𝑥,𝑡)) be the displacement vector function of an inhomogeneous anisotropic elastic material characterizing by density 𝜌 and the elastic moduli 𝐶_{𝑖𝑗𝑘𝑙}. The density 𝜌 and elastic moduli 𝐶_{𝑖𝑗𝑘𝑙} are varying functions of position 𝑥=(𝑥₁,𝑥₂,𝑥₃). Combining the properties of the strain-energy function with Hooke's law we find [102] that 𝐶_{𝑖𝑗𝑘𝑙} satisfy the following property and strong convexity 𝐶_{𝑖𝑗𝑘𝑙}=𝐶_{𝑘𝑙𝑖𝑗}=𝐶_{𝑙𝑘𝑖𝑗}=𝐶_{𝑖𝑗𝑙𝑘},³³_𝑖=1³_𝑗=1³_𝑘=1𝑙=1𝐶_{𝑖𝑗𝑘𝑙}𝜀_𝑖𝑗𝜀_𝑘𝑙>0 (29) for any nonzero 3×3 real symmetric matrix (𝜀_𝑖𝑗)_3×3. Equations for motion in inhomogeneous anisotropic elastic materials are, in our notation (see, e.g., [102]), 𝜌𝜕²𝑢_𝑖𝜕𝑡²=³³_𝑗=1³_𝑘=1𝑙=1𝜕𝜕𝑥_𝑗𝐶_{𝑖𝑗𝑘𝑙}𝜕𝑢_𝑘𝜕𝑥_𝑙+𝑓_𝑖,𝑖=1,2,3, (30) where 𝑓_𝑖 are components of the body forces 𝐟(𝑥,𝑡)=(𝑓₁(𝑥,𝑡),𝑓₂(𝑥,𝑡),𝑓₃(𝑥,𝑡)) acting per unit volume on the particle originally at position 𝑥 at some reference time 𝑡.

Let us examine 𝐮(𝑥,𝑡)=𝐔(𝑡−𝑇(𝑥))𝐴(𝑥) (31) as an approximate solution of (30) for 𝐟=0. Near a wavefront 𝑡=𝑇(𝑥), we assume that components of 𝐔=(𝑈₁,𝑈₂,𝑈₃) are fluctuating much more rapidly than 𝐴(𝑥) or 𝐶_{𝑖𝑗𝑘𝑙}, and the successive derivatives 𝜕𝐔/𝜕𝑡 and 𝜕²𝐔/𝜕𝑡² are fluctuating still more rapidly. Substituting (31) into (30) we find (see, e.g., [102]) 𝜌𝛿_𝑖𝑘−³³_𝑗=1𝑙=1𝜕𝑇𝜕𝑥_𝑗𝜕𝑇𝜕𝑥_𝑙𝜕²𝑈_𝑘𝜕𝑡²𝐴=𝐸_𝑖(𝐔𝐴), (32) where 𝐸_𝑖(𝐔𝐴) includes merely first-order derivatives of 𝐔, 𝐔 itself, the elastic moduli, amplitude function 𝐴(𝑥), and gradients of these.

Thus, the left-hand side of (32) must be much smaller than 𝜕²𝐔/𝜕𝑡². We conclude that the matrix of coefficients of 𝜕²𝐔(𝑡−𝑇(𝑥))/𝜕𝑡²𝐴(𝑥) must be singular: det𝜌𝛿_𝑖𝑘−³³_𝑗=1𝑙=1𝜕𝑇𝜕𝑥_𝑗𝜕𝑇𝜕𝑥_𝑙=0. (33) This equation determines the possible wavefronts in an elastic medium, since it gives a constraint on the function 𝑇(𝑥). In an inhomogeneous isotropic medium, 𝐶_{𝑖𝑗𝑘𝑙}=𝜆𝛿_𝑖𝑗𝛿_𝑘𝑙+𝜇𝛿_𝑖𝑘𝛿_𝑗𝑙+𝛿_𝑖𝑙𝛿_𝑗𝑘, (34) where 𝛿_𝑖𝑗 is the Kronecker symbol; that is, 𝛿_𝑖𝑗=1 if 𝑖=𝑗 and 𝛿_𝑖𝑗=0 if 𝑖≠𝑗; moreover, 𝜆=𝜆(𝑥), 𝜇=𝜇(𝑥) are known as the Lame functions.

In an inhomogeneous isotropic medium, the special form (34) of 𝐶_{𝑖𝑗𝑘𝑙} makes it possible to get (33) as follows: 𝜌∇𝑇⋅∇𝑇−𝜌𝜆+2𝜇∇𝑇⋅∇𝑇−𝜇²=0. (35) This is, 𝑇(𝑥) satisfies the eikonal equation |∇𝑇|²=1𝐶²_𝑃(𝑥) (36) or eikonal equation |∇𝑇|²=1𝐶²_𝑆(𝑥), (37) where 𝐶²_𝑃=(𝜆+2𝜇)/𝜌 is the local 𝑃-wave speed and 𝐶²_𝑆=𝜇/𝜌 is the local 𝑆-wave speed.

3.2. Inverse Kinematic Problem of Seismic

Suppose that a point source at position 𝑥⁰∈ℝ³ becomes active at a time chosen to be the origin, 𝑡=0. In homogeneous isotropic medium, wavefronts emanate from the source as ever-expanding spheres, with radius 𝐶_𝑃𝑡 (for 𝑃-waves) and 𝐶_𝑆𝑡 for 𝑆-waves, arriving at the general position 𝑥 at time 𝑡=|𝑥−𝑥⁰|/𝐶_𝑃 and 𝑡=|𝑥−𝑥⁰|/𝐶_𝑆. We introduce the function 𝑇(𝑥,𝑥⁰) as the travel time required for the wavefront to reach 𝑥 from 𝑥⁰. The function 𝑇(𝑥,𝑥⁰) satisfies (36) for 𝑃-waves and (37) for 𝑆-waves.

One of the first inverse problem, stated and studied in geophysics, was the inverse kinematic problem. The physical interpretation of this problem is the following. Let us assume that Earth is an isotropic inhomogeneous elastic medium and the measurements of the seismic waves, arising from a point source 𝑥⁰ and propagating in Earth, are given for points on its surface. These measurements contain data of the travel time 𝑇(𝑥,𝑥⁰) of seismic waves between the point of the source 𝑥⁰ and any point of the Earth's surface. The inverse kinematic problem is to find the speed of the seismic waves inside of Earth using the measurement data. Mathematically we can state the inverse kinematic problem as follows. Let 𝐷 be a domain bounded by the surface 𝑆, and let 𝑇(𝑥,𝑥⁰) be the function of the travel time required by a signal with unknown speed 𝐶(𝑥)>0 to reach 𝑥 from 𝑥⁰. Find 𝐶(𝑥) for all 𝑥 from 𝐷 if the function 𝑇(𝑥,𝑥⁰) is given for all points 𝑥⁰∈𝑆₁ and 𝑥∈𝑆₂, where 𝑆₁⊆𝑆 and 𝑆₂⊆𝑆 are subsets of 𝑆.

Herglotz [103] and Wiechert and Zoeppritz [104] were the first who studied the inverse kinematic problem in assumptions 𝐶(𝑥)1𝑛(𝑟)|𝑥|,𝑑,𝑟=((𝑟𝑥𝑑𝑟𝑟𝑛))>0,₀isafixedpointfrom𝑆, (38) if 𝑇(𝑥,𝑥⁰) is known for any 𝑥 from 𝑆. Gerver and Markushevich [105] have showed that the condition (𝑑/𝑑𝑟)(𝑟𝑛(𝑟))>0 can be illiminated, but in this case the inverse kinematic problems have many solutions and the set of these solutions has been described. The first theoretical study of the inverse kinematic problem for a horizontal inhomogeneous medium has been made by Lavrenti'ev and Romanov [106]. The first result of the study of the multidimentional inverse kinematic problem in a linear approximation, when the function 𝐶(𝑥) depends on 3D space variable 𝑥=(𝑥₁,𝑥₂,𝑥₃)∈ℝ³, was obtained by Romanov [107]. The existence, uniqueness, and stability estimate theorems for the multidimensional inverse kinematic problems have been established in some classes of analytic and smooth functions by Anikonov [108], Bukhgeĭm [109], and Mukhometov [110].

In a recent time, the Earth is modeled as an anisotropic elastic medium which is located in the given 3D domain. The wave speed is given by a symmetric positive definite matrix 𝐺=(𝑔_𝑖𝑗)_3×3, that is, a Riemannian metric in mathematical terms (see, e.g., [111]). The problem is to determine the metric in a given domain from the lengths of geodesics joining points on the boundary of the domain. The linearization of this problem leads to a problem of the integral geometry [111–113].

The regular study of the problems of finding the isotropic and anisotropic Riemannian metrics and integral geometry problems has been made in the works [111, 114–117]. The modern numerical algorithms for the computation of the inverse kinematic problems of seismic have been developed in the works [118, 119].

3.2.1. Inverse Problems in Ray Statements

Let us note that isotropic inhomogeneous elastic medium is completely characterized by three functions: the density 𝜌(𝑥) and speeds 𝐶_𝑃(𝑥), 𝐶_𝑆(𝑥) of 𝑃- and 𝑆-waves. Using the point source at position 𝑥⁰ of the boundary 𝑆 of the given domain 𝐷, which becomes active at the time 𝑡=0, we measure the function 𝑇(𝑥,𝑥⁰) of the travel time required for the fronts of 𝑃- and 𝑆-waves to reach 𝑥 from 𝑥⁰. We use these information for solving two inverse kinematic problems for 𝑃- and 𝑆-waves. The solutions of these problems are speeds 𝐶_𝑃(𝑥), 𝐶_𝑆(𝑥). To complete the identification of unknown isotropic inhomogeneous medium we need to determine the last unknown function 𝜌(𝑥) after finding 𝐶_𝑃(𝑥), 𝐶_𝑆(𝑥). An inverse problem to recover 𝜌(𝑥) in a given bounded domain 𝐷, containing an isotropic inhomogeneous elastic medium, has been solved by Yakhno [120–123]. In these works, the displacement fields 𝐮(𝑥,𝑥⁰,𝑡) have been measured for all points 𝑥 and 𝑥⁰ running the boundary 𝑆 of 𝐷 for all times from a time interval containing the time 𝑡=𝑇(𝑥,𝑥⁰) of arriving of the 𝑃-waves.

3.3. One-Dimensional Inverse Dynamic Problems

The vertical inhomogeneous model of Earth is one of the popular models of geophysics [102], and the inverse problems of recovering the density 𝜌(𝑥₃) and Lame functions 𝜇(𝑥₃) and 𝜆(𝑥₃), depending on one variable 𝑥₃ and appearing in equations of elastodynamics (30) for the case of inhomogeneous isotropic elastic medium, have been studied by many authors [1, 2, 7–15, 21, 22, 24, 26–28, 36–46, 53–55, 59–61, 70–73, 78–81, 85, 87–91, 97, 98, 103, 105, 108, 109, 113, 118, 119, 124–214]. Because the unknown functions depend on one variable, the inverse problems of their recovering are called one-dimensional inverse problems although all differential equations of elasticity contain 3D space variable 𝑥=(𝑥₁,𝑥₂,𝑥₃) and 1D time variable 𝑡. Alekseev and Dobrinsky [124, 125] were the first who described the importance of one-dimensional inverse problems of elastodynamics in geophysics and studied them as problems of the recovery of smooth functions 𝜌(𝑥₃), 𝜇(𝑥₃), and 𝜆(𝑥₃) of one variable 𝑥₃. The uniqueness of the solutions of these inverse problems has been studied firstly by Blagoveschenskii [147] and Romanov [215] for the isotropic elastic media and then by Volkova and Romanov [216] for anisotropic elastic media. The regular study of the theory, methods, and applications of one-dimensional inverse problems for dynamical differential equations of isotropic and anisotropic elastic media has been made in works [122, 159, 170, 190, 217–221] and others. The recent development of theory, methods, and applications of one-dimensional inverse problems of dynamic elasticity can be found in the works [150, 211, 222, 223].

We note that a model of Earth as a composite medium consisting of a finite number of different elastic layers is very popular in geophysics. In this case, the one-dimensional inverse problem consists of finding 𝜌(𝑥₃) and the Lame functions 𝜇(𝑥₃) and 𝜆(𝑥₃) as functions of one variable 𝑥₃ with piecewise constant values. The computation of solutions of this type of one-dimensional inverse problems has been studied in [186, 187, 220]. The modern theory and methods of the construction of solutions of one-dimensional inverse problems for the equations of elastodynamics in elastic composite media can be found in the works [142, 214, 224–228].

3.4. Multidimensional Inverse Dynamic Problems in Linear (Born) Approximation

Linearized multidimensional inverse dynamic problems (or inverse problems in the Born approximation) take an important place through all statements of multidimensional inverse problems for equations of elastodynamics. The statements of these problems have natural physical and mathematical sense. From the physical point of view, an isotropic inhomogeneous elastic body, which is characterized by the Lame functions 𝜇₁(𝑥) and 𝜆₁(𝑥) and density 𝜌₁(𝑥), is included in a vertical inhomogeneous (or homogeneous) elastic medium. Let, for example, the half space 𝑥₃>0 contain this medium, and let the characteristics 𝜇₁(𝑥), 𝜆₁(𝑥), and 𝜌₁(𝑥) of the elastic body be unknown functions. The linearized inverse problem is to find these unknown functions if we measure the first act of scattering the displacement field on the surface 𝑥₃=0 arising from the forces located on the same surface 𝑥₃=0. From the mathematical point of view, we consider the differential equations of elastodynamics in a half space 𝑥₃>0 with boundary conditions on 𝑥₃=0. We assume that the Lame functions 𝜇(𝑥) and 𝜆(𝑥) and density 𝜌(𝑥) appearing in differential equations and boundary conditions can be presented in the form 𝜇(𝑥)=𝜇₀𝑥₃+𝜇₁(𝑥),𝜆(𝑥)=𝜆₀𝑥₃+𝜆₁(𝑥),𝜌(𝑥)=𝜌₀𝑥₃+𝜌₁(𝑥), (39) where 𝜇₀(𝑥₃), 𝜆₀(𝑥₃), and 𝜌₀(𝑥₃) are functions depending on 𝑥₃ and characterizing vertical inhomogeneous medium and 𝜇₁(𝑥), 𝜆₁(𝑥), and 𝜌₁(𝑥) are functions of 3D space variable 𝑥=(𝑥₁,𝑥₂,𝑥₃) characterizing the elastic body which is included in the vertical inhomogeneous medium.

We assume that the displacement field 𝐮(𝑥,𝑡) is presented in the form 𝐮(𝑥,𝑡)=𝐮₀(𝑥,𝑡)+𝐮₁(𝑥,𝑡), where 𝐮₀(𝑥,𝑡) is the displacement field of the vertical inhomogeneous medium arising from the given forces, and 𝐮₁(𝑥,𝑡) is the first act of scattering 𝐮₀(𝑥,𝑡) on the inhomogeneous inclusion with characteristics 𝜇₁(𝑥), 𝜆₁(𝑥), and 𝜌₁(𝑥). The equations of elastodynamics with boundary conditions are linearized around 𝜇₀(𝑥₃), 𝜆₀(𝑥₃), 𝜌₀(𝑥₃), and 𝐮₀(𝑥,𝑡). The unknown functions 𝜇₁(𝑥), 𝜆₁(𝑥), and 𝜌₁(𝑥) appear in inhomogeneous terms of linearized equations for 𝐮₁(𝑥,𝑡). We need to recover 𝜇₁(𝑥), 𝜆₁(𝑥), and 𝜌₁(𝑥) if we know 𝐮₁(𝑥,𝑡) for 𝑥₃=0 (see, e.g., [122]). The uniqueness of the solution of a linearized multidimensional inverse problem has been studied by Romanov [215]. The existence theorem and computation of a solutions of a linearized multidimensional inverse problems of elastodynamics have been obtained in the works [122, 229]. The recovery of the function characteristics of an elastic body in linear approximation was a subject of the works [146, 154, 202, 230]. The linearized inverse problems of determining the function characteristics of transversally isotropic elastic media from the measurements of reflected waves have been developed by Sharafutdinov [231, 232]. The linearized inverse problems for nonhomogeneous bodies have been stated and developed by Steinberg [233].

3.5. Multidimensional Inverse Dynamic Problems in the Statements of the Dirichlet-To-Neumann Map

The inverse problems of determining the elastic moduli and density as functions of the space variables in a bounded domain from observed data of the solution on the boundary (or a part of the boundary) of this domain are geophysical motivated. One important class of these problems is inverse problems for equations of elastodynamics in terms of the Dirichlet-to-Neumann map. The Dirichlet-to-Neumann map models surface measurements by giving the correspondence between a displacement at the boundary 𝑆 of the given bounded domain 𝐷 and surface traction ³³_𝑗=1³_𝑘=1𝑙=1𝜈_𝑗𝐶_{𝑖𝑗𝑘𝑙}𝜕𝑢_𝑘𝜕𝑥_{𝑙𝑆×(0,𝑇)},𝑖=1,2,3, (40) where 𝜈=(𝜈₁,𝜈₂,𝜈₃) is the unit outer normal to 𝑆, (0,𝑇) is the observed time interval, and 𝑢_𝑘, 𝑘=1,2,3 are components of the displacement vector function 𝐮(𝑥,𝑡) satisfying (30). The details of the use of the Dirichlet-to-Neumann map in modeling surface measurements in inverse problems can be found in [234]. The inverse problems in the Dirichlet-to-Neumann map statements are successfully applied to study the unique determination of the solutions of the inverse problems of elasticity as in the static isotropic and anisotropic cases [235–239] as well as in dynamic case [240].

3.6. Carleman Estimates and Uniqueness of Solutions for Multidimensional Inverse Dynamic Problems

For the study of the inverse problems for the scalar partial differential equations with a finite number of observation, Bukhgeim and Klibanov [160] proposed a remarkable method based on a Carleman estimate. Later, the Carleman estimate method has been generalized to study the uniqueness and stability estimate of the solutions of the inverse problems for equations of elastodynamics by Isakov [208], Ikehata et al. [203], Eller et al. [182], Imanuvilov et al. [204], Imanuvilov and Yamamoto [205, 206], and Lin and Nakamura [241].

In the papers [203–206, 208] the authors assume some geometric constraints on the surface under observation for proving the uniqueness and stability of the solutions of the multidimensional inverse problems using the Carleman estimates. Later, the stability estimate theorem for solutions of a multidimensional inverse problem for equations of elastodynamics has been proven for an arbitrary subboundary by Bellassoued et al. [143].

4. Inversion Based Applications

In the following, we present a general overview on some inversion-based application areas. We then give more details about particular selected applications in the subsequent sections. Further, different applications whose details are skipped in this section such as remote sensing, nuclear science, and geophysics can be found in [164, 184, 242–247].

Underwater Acoustics and Traveltime Tomography. Inverse problems related to underwater acoustics are a critical research area due to their wide range of important practical applications. Solutions of such problems can provide a heat distribution of seawater, an ocean depth at a point, sediment properties of a seabed, locations and/or shapes of submerged objects (mines, submarines, sunken wrecked ships buried to the seafloor, etc.), positions of underwater cables (communication, gas, petroleum, waste, etc.), profiles of seamounts, and a layered seabed for seismic applications or paths for current flows and migrations of sea animals [21]. In principle, all these applications can roughly be classified into two groups of problems such as remote sensing of passive marine environments and localization/characterization of acoustic sources which are in the sea or buried in the seafloor [171]. For the solution of these problems, one can use the speed of acoustic wave as a tool which is obtained from traveling time of an acoustic wave between two certain points or the amplitude and/or phase information of acoustic waves reflected/scattered from a target(s). In water-type medium, propagating acoustic waves are collected at a selected number of separate hydrophones to obtain measured field data.

Generally, reconstructions of desired parameters from the knowledge of scattered acoustic waves lead to nonlinear and ill-posed inverse problems. Therefore, it is always a complicated issue to find unique and stable solutions, and one has to apply some additional techniques, that is, regularizations, for the proper treatment of the ill-posedness. More specifically, researchers have been applying the methods of Ray’s theory, modal travel time/phase inversions, boundary integral equation approaches, Rayleigh’s hypothesis, linear sampling, complete family method, convex scattering support theory, and so forth, in underwater acoustics [21, 71, 153, 165, 171, 189, 248–252].

Nondestructive Testing. The inspection of an object without touching or without changing its characteristic properties is a general definition of nondestructive testing (NDT) in the literature. In practice, NDT is commonly performed by following visual, penetrant, magnetic particle, radiographic, thermal infrared, Eddy’s current, and ultrasonic type testing procedures for the identifications of cracks, shapes of surface discontinuities, or corrosion damages in power plants, rails, pipelines, bridges, buildings, airplanes, railroad tank cars, and so forth, [196, 253]. The main idea of these methods, which actually shows the conceptual relation between inverse problems and NDT applications, is the evaluation of the responses obtained from the interaction of a test object with a particular effect, that is, electromagnetic/acoustic field, gamma or X-radiation, fluorescent dye, and so forth. Therefore, approaches which are used for inverse problems, that is, synthetic aperture focussing technique (SAFT, [180, 254]), diffraction tomography (DT, [177–179, 255–257]), multiple signal classification (MUSIC, [129, 155, 167, 176, 198]), linear sampling method (LSM, [70–74]), factorization [75, 76, 258], point source [23, 68, 69], no response test [82–84], and so forth, are also employed for NDT problems [259].

More specifically, SAFT is an algorithm which uses the collection of echo signals over a specific aperture to obtain a reconstruction by performing time shifting and superposition of adjacent signals. DT is based on a linear solution of the wave equation which can be obtained via the Born or Rytov approximations. Here, the linearization approaches define mainly the success and the solution space of the inverse problem. MUSIC was initially employed as a direct imaging algorithm to obtain locations of point scatterers [176] and extended to find also the geometry of targets [198]. The method employs the eigenvalue structure of time-reversal matrix which is obtained from measured data at different receiver antennas. The main idea of linear sampling method is to find an indicator function such that its value provides whether an arbitrarily tested space coordinate lies inside or outside of the object.

Biomedical Imaging Techniques, Tomography, and EIT. Nowadays medical imaging, which can be considered as one of the most developed area of inverse problems in practice, provides high-resolution reconstructions in the order of millimeters. In the last decays, the main effort is given for the implementation of harmless, fast, cheap, robust, and reliable techniques to use in practice for obtaining high-resolution images in real time. In the similar direction, early studies of bioimaging started with the reconstructions of 2D images of human body parts via the inverse Radon transform [145, 162, 260] of measured X-rays which were attenuated inside the body [168, 261]. Afterwards, 3D images were assembled via computed tomography (CT-scan) from a series of X-ray data measured on different planes (sinograms) in 2D [144, 199, 212, 262, 263]. Even though satisfactory results were obtained with the X-ray radiology especially for the bone structures [131, 264] the method was found not sufficiently efficient due to attenuation characteristics of X-rays and harm risks of using ionized radiation on humans. On the other hand, electrical impedance tomography (EIT), after its main idea and formulation were introduced by Calderon in 1980 [163] and D. Isaacson and E. L. Isaacson [207], has gained high interest both from theoretical and physical point of view. In principle, in EIT low-frequency electrical currents are applied to the body part under investigation. Then electrical properties of body tissues are computed from the measurements of electric currents and voltage at the boundary [126, 136, 151, 152, 156, 166, 197, 265–268]. EIT is successfully applied for diagnosis of breast cancer, monitoring brain and gastrointestinal functions, detection of blood clots in the lungs, and so forth [138]. Furthermore, electroencephalography (EEG) and magnetoencephalography (MEG) are used for passive monitoring of neuron activities in the brain from the weak electric or magnetic fields, respectively, [130, 137, 169, 191, 269, 270].

MRI, PET, and SPECT. A different approach which is based on using properties of subatomic particles with the connection to electromagnetism opened a new area for obtaining high spatial resolutions in bioimaging, for example, magnetic resonance imaging (MRI), positron emission tomography (PET), and single-photon emission computed tomography (SPECT). In MRI, the patient stays in a tunnel under a strong magnetic field typically 0.2–1.5 T. This large static magnetic field aligns protons of many atoms either parallel or antiparallel existing in the body. In the meantime, weak radio frequency fields are applied systematically to the patient for altering the alignment of the magnetization. As a result of this procedure, rotating magnetic fields induce a voltage at the receiver coils of the magnet which is used to reconstruct the image of the scanning area [131, 192, 271–274]. In PET and in SPECT, chosen molecules are labeled with radioactive atoms having short half-life and injected to the patients’ bloodstream at very low concentrations in order not to violate radiation exposure limits. For PET, labeled atoms are chosen to emit positrons, and for SPECT, they are chosen to emit photons when they are decaying. Gamma rays, which occur when an electron and an emitted positron annihilate in PET, and photons which are released in SPECT can be visualized out of body by using scintigraphic detectors for clinical applications of oncology, cardiology, pharmacology, and so forth [127, 175, 275–282].

MW/Ultrasound Tomography, Optical Imaging, and Cognition. Another group of techniques such as microwave tomography, ultrasound, and optical imaging, which are commonly used for solutions of inverse problems in different areas, are also applied to biomedical applications especially for investigating soft body tissues [135, 148, 149, 185, 210, 283–289]. For instance, optical tomography is used for the detection of cancerous cells in breast and brain. Acoustic waves are employed for the imaging of liver, kidney, fetus in pregnant women, and so forth, and microwaves are used in mammography and diagnosis of leukemia [131, 144].

Data Assimilation. Over the past two decades, it has became feasible to simulate atmospheric and geophysical processes from large-scale atmospheric flow down to convective processes on a kilometer scale. This led to the need to determine initial conditions for simulations and forecasts from a collection of diverse direct and indirect measurements, and the field of data assimilation arose.

Inverse Problems in Biological and Environmental Applications. Inverse problems are of growing importance in many parts of medicine or biology as well as in environmental applications. Here, we will provide a brief introduction into two areas, first into recent results of neural field theory and second into the basic setup of data assimilation as it is, for example, used in numerical weather prediction or for climate projections, which usually incorporates various inverse problems.

4.1. Inverse Neural Field Theory

Neural activity is often modelled by the activity potential 𝑢(𝑥) in some domain 𝐷. The activity potential 𝑢 satisfies some integro-differential equations, which in its simplest form has been suggested by Wilson and Cowan [290, 291] and Amari [128]: 𝜏𝜕𝑢(𝑥,𝑡)()+𝜕𝑡=−𝑢𝑥,𝑡_𝐷𝑤𝑥,𝑦𝑓𝑢𝑦,𝑡𝑑𝑦,𝑥∈𝐷,𝑡≥0, (41) where 𝑓 is some nonlinear function which usually is used to model the firing rate of neurons in dependence of their current action potential and 𝜏>0 is a constant. The kernel 𝑤(𝑥,𝑦) models the strength of the influence of an excitation at point 𝑦 to the neural field at point 𝑥. The first term of the right-hand side of (41) generates the decay of the activity potential 𝑢 in the absence of excitation or inhibition. Neural fields have been widely studied in recent years, with applications to a wide range of medical phenomena starting from electroencephalogram (EEG) and magnetoencephalogram (MEG) rhythms to robotic behaviour—for an extensive literature list we refer to [172].

Neural network, that is, the discrete version of neural fields, has attracted strong interest over many years and is a standard tool today in many applied parts of science. Neural networks have also been used to solve inverse problems, compare, for example, [292, 293] for some recent papers and further citations. However, here we want to look into inverse problems which arise in the modelling of neural activities themselves. Training of neural networks or neural fields is, in general, an ill-posed inverse problem, as we will see in due course.

Inverse neural field theory investigates the construction of connectivity kernels 𝑤 given some dynamics 𝑢(𝑥,𝑡) for 𝑥∈𝐷 and some time intervals 𝑡∈[0,𝑇]. This so-called full-field neural inverse problem is linear and ill-posed in the sense of Hadamard, as can be readily seen by the following transform. We define 𝜓()()𝑥,𝑡∶=𝜏𝜕𝑢𝑥,𝑡()[],𝜑𝑢[].𝜕𝑡+𝑢𝑥,𝑡,𝑥∈𝐷,𝑡∈0,𝑇𝑦,𝑡∶=𝑓𝑦,𝑡,𝑦∈𝐷,𝑡∈0,𝑇 (42) If we further define the operator 𝑊𝜑(𝑥)∶=_𝐷𝑤𝑥,𝑦𝜑𝑦𝑑𝑦,𝑥∈𝐷, (43) equation (41) obtains the form 𝜓(⋅,𝑡)=𝑊𝜑(⋅,𝑡),𝑡∈[0,𝑇]. (44) The task is to find the operator 𝑊 given the family of states 𝜑(⋅,𝑡) and corresponding images 𝜓(⋅,𝑡) for 𝑡∈[0,𝑇]. Changing our perspective slightly, introducing the operator 𝐾𝑔(𝑡)∶=_𝐷𝜑𝑦,𝑡𝑔𝑦𝑑𝑦,𝑡∈[0,𝑇], (45) we transform (41) or (44) into 𝜓(𝑥,𝑡)=𝐾𝑤(𝑥,⋅)(𝑡),𝑡∈[0,𝑇] (46) for each fixed 𝑥∈𝐷. For each 𝑥∈𝐷, (46) is an integral equation of the first kind for the function 𝑤(𝑥,⋅). Usually we need some smoothness of the potential 𝑢 with respect to its arguments 𝑦∈𝐷 and 𝑡∈[0,𝑇]. In this case, the operator 𝐾 is an integral operator with continuous kernel, which is known to be compact in either 𝐿²(𝐷) or 𝐶(𝐷) (cf. [22]) and cannot have a bounded inverse. Thus, the inverse neural field problem is ill-posed.

The particular form (46) provides a basis for kernel construction, that is, for the solution of the neural field problem. We may apply spectral regularization schemes as described in [1, 22, 85], for example, Tikhonov regularization in a similar sense of (21) 𝑅_𝛼∶=𝛼𝐼+𝐾^∗𝐾⁻¹𝐾^∗,𝛼>0, (47) with regularization parameter 𝛼>0 as regularized inverse to calculate 𝑤(𝑥,⋅) from the knowledge of 𝜓 and 𝜑. For smooth dynamics, the problem is exponentially ill-posed. We refer to [139–141, 294] for the analysis and many examples for the inverse neural field problem. The problem of the ill-posedness of the neural inverse problem is addressed in [295], where a dimensional reduction approach is suggested to decompose the large and strongly ill-posed full problem into more stable individual tasks.

4.2. Data Assimilation with an Application to Numerical Weather Prediction (NWP)

In inverse problems, we are usually given measured data, and the task is to gain insight into some unknown parameter functions insight of an inaccessible body or area of space or to reconstruct the shape of scatterers or inclusions, but often, the quantities under reconstruction are not static but dynamic and change over time. Then, the inverse task is not only carried out once but the reconstruction is repeated over time with some cycling frequency given by a time interval 𝛿𝑇.

Let 𝑓_𝑘∈𝑌 for 𝑘∈ℕ be an element of a Hilbert space 𝑌 representing our measurement data. The task is to reconstruct some state 𝑥_𝑘∈𝑋 in a Hilbert space 𝑋, where the measurement is described by an observation operator 𝐻_𝑘∶𝑋→𝑌. An underlying dynamical system is given by some operator 𝑀_𝑘∶𝑋→𝑋, mapping the state 𝑥_𝑘 at time 𝑡_𝑘 onto the state 𝑥_𝑘+1 at time 𝑡_𝑘+1. With the reconstruction 𝑥^(𝑎)_𝑘−1 from the previous time step, we can calculate a state 𝑥_𝑘^(𝑏) at time 𝑡_𝑘, which serves as a priori knowledge for the current reconstruction when data 𝑓_𝑘 are given.

The goal of data assimilation in its simplest form is to employ measurement data 𝑓_𝑘, the operator 𝐻_𝑘, and the background state 𝑥_𝑘 to calculate a so-called analysis 𝑥_𝑘^(𝑎), which is the best possible estimate of the true state 𝑥^(true)_𝑘 of the dynamical system at time 𝑡_𝑘 under the given assumptions.

Data assimilation is needed as soon as dynamic situations are studied. For example, the reconstruction of a current density from magnetic fields leads to static magnetic tomography, as, for example, solved in [193–195, 296], but when the underlying currents show a dynamic behaviour, we need to carry out dynamic magnetic tomography [297, 298], which is a data assimilation problem.

The data assimilation scheme can also use measurements from an interval of time steps, for example, a window consisting of 𝐿 time steps with data 𝑓_𝑘−𝐿,𝑓_{𝑘−𝐿+1},…,𝑓_𝑘. In this case, the task is to fit the trajectory given by 𝐻𝑥_𝑘−𝐿,𝐻𝑥_{𝑘−𝐿+1}=𝐻𝑀_𝑘−𝐿𝑥_𝑘−𝐿,𝐻𝑥_{𝑘−𝐿+2}=𝐻𝑀_{𝑘−𝐿+1}𝑀_𝑘−𝐿𝑥_𝑘−𝐿,…,𝐻𝑥_𝑘=𝐻𝑀_𝑘−1𝑀_𝑘−2⋯𝑀_𝑘−𝐿𝑥_𝑘−𝐿 (48) to the measured data 𝑓_𝑘−𝐿,…,𝑓_𝑘 given the a priori information 𝑥^(𝑏)_𝑘−𝐿.

Let us denote the analysis or reconstruction, respectively, for this window [𝑡_𝑘−𝐿,𝑡_𝑘] at time 𝑡_𝑘−𝐿 by 𝑥^(𝑎,𝐿)_𝑘−𝐿. Then, it leads to a state estimate at time 𝑡_𝑘 by 𝑥_𝑘^(𝑎)∶=𝑀_𝑘−1𝑀_𝑘−2⋯𝑀_𝑘−𝐿𝑥^(𝑎,𝐿)_𝑘−𝐿,𝑘∈ℕ. (49) An algorithm, which calculates 𝑥^(𝑎)_𝑘−𝐿 at the beginning of the window [𝑡_𝑘−𝐿,𝑡_𝑘], is usually denoted as smoother. The calculation of 𝑥_𝑘^(𝑎) at the end of the interval is called a filter. There is a wide range of literature about data assimilation, in particular from the perspective of filtering. We refer to [133, 134, 173, 209, 213, 299–304] for introduction and further details, in particular in atmospherical sciences.

We can reformulate the data assimilation task as a minimization, here written in its usual finite-dimensional version. If we carry out minimization independently at every time 𝑡_𝑘, the scheme is known as three-dimensional variational assimilation (3dVar) and minimizes 𝐽𝑥_𝑘∶=𝑥_𝑘−𝑥_𝑘^(𝑏)2_𝐵−1+𝑓_𝑘−𝐻𝑥_𝑘²_𝑅−1, (50) where 𝐵 and 𝑅 are weight matrices in 𝑋 and 𝑌, respectively. In a stochastic framework, 𝐵 is the covariance matrix of background distribution of states in 𝑋 and 𝑅 is the covariance matrix of the data error in 𝑌.

The 3dVar algorithm is basically a version of the well-known Tikhonov regularization in a Hilbert space where weighted norms are used. Also, it can be seen as a version of the Bayes approach for a Gaussian densities; compare [209] for an introduction and [188] for a recent survey.

When a window of measurements is used, the corresponding minimization problem is given by 𝐽𝑥_𝑘−𝐿=𝑥_𝑘−𝐿−𝑥^(𝑏)_𝑘−𝐿²_𝐵−1+^𝐿_ℓ=1𝑓_{𝑘−𝐿+ℓ}−𝐻𝑥_{𝑘−𝐿+ℓ}²_𝑅−1, (51) where 𝑥_{𝑘−𝐿+ℓ} is given by (48). The minimization of (51) is denoted as four-dimensional variational assimilation (4dVar). The 4dVar algorithm is known to be expensive both in terms of computing time and programming effort, since usually the adjoint tangent linear model of the dynamics needs to be calculated and implemented, but 4dVar has turned out to be the algorithm which provides best scores in the area of numerical weather forecasting (NWP); compare the scores which are available from the European Centre for Medium Range Weather Forecast (ECMWF) [174].

For linear systems the minimization of the 4dVar functional (51) can be carried out iteratively with individual steps as in (50). The key difference to 3dVar is a dynamic update of the covariance matrix 𝐵 in every time step. The method is the well-known Kalman Filter (KF) [213]. The Kalman filter is often studied in a stochastic framework, for example, [209, 304], where it calculates the maximum probability estimator (MAP) for the Bayes formula in a Gaussian framework and when the prior density has a covariance matrix 𝐵 and mean 𝑥_𝑘^(𝑏). A simple algebraic equivalence proof of the Kalman smoother/filter to 4dVar is presented, for example, in [188]. We need to remark that for large-scale systems, the update of the 𝐵 matrix in the Kalman filter is not feasible, which has led to the development of ensemble and particle methods.

Ensemble or particle methods which are well-known in stochastic estimation since around 1980 (Markov Chain Monte Carlo methods, MCMC) have become very popular in the area of data assimilation. The basic idea here is to estimate the covariance matrix 𝐵 by a stochastic covariance estimator, based on a set of ensembles or particles. Then, an assimilation step as in (50) is carried out. The method is known as ensemble Kalman filter (EnKF). Here, a particle or ensemble member, respectively, is basically a state 𝑥_𝑘−1, and its trajectory 𝑥_𝑘=𝑀_𝑘−1𝑥_𝑘−1. When a distribution of states is given, we can calculate mean and covariance by standard tools. The propagation of these states through time using the dynamical system 𝑀 leads to an estimate of mean and covariance of the propagated distribution at a later point in time.

The estimation of the covariance matrix 𝐵 by an ensemble with 40, 100, or 200 members in a high-dimensional space leads to spurious covariances, which make the ensemble approach useless for such systems. As a consequence, ensemble data assimilation systems have employed localization [132, 181, 200, 305, 306]. Here, calculations are carried out in a local region only, not the global space [307, 308]. The local ensemble Kalman filter (LEKF) is studied in [306], and its computationally more efficient version has been developed by Hunt et al. [201], the local ensemble transformed Kalman filter (LETKF).

More recently, there has been increasing interest in the study of stability of data assimilation algorithms over time, in particular when ill-posed operators are involved in the observation process; compare [157, 158, 309, 310]. The key object here is to study the mapping of the data sequence onto the analysis sequence 𝑓_{𝑘𝑘∈ℕ}↦𝑥_𝑘^(𝑎)_𝑘∈ℕ. (52) Stability—estimating ‖𝑥_𝑘^(𝑎)‖ for 𝑘→∞—has been derived under particular conditions, leading to tools to control the stability over time by scaling the background covariance matrix 𝐵 (or alternatively the data covariance matrix 𝑅) that have been developed [311].

5. Conclusions

In this paper, inverse problems are reviewed starting with a historical perspective. Mathematical backgrounds for the fundamental inverse problems in acoustics, electromagnetics, and elastics are given. The main ideas of many algorithms whose aims to find locations, shapes, and boundary-parameters of obstacles are introduced by studying two dimensional acoustic cases for the sake of simplicity of the theoretical investigations.

Furthermore, the results of the study of the inverse problems for the time-dependent differential equations of isotropic and anisotropic elasticity are presented chronologically from one-dimensional kinematic and dynamic inverse problems with essentially over-determined data up to multidimensional inverse problems with a finite number of observations. The physical and geophysical interpretations are supplied for almost all reviewed works.

We note that the study of one-dimensional and linearized multidimensional inverse problems for equations of elastodynamics contains the theory (existence, uniqueness, and stability theorems) and computational methods of finding solutions. The study of the multidimensional dynamic problems for the differential equations of elasticity contains mainly the investigations on the uniqueness of the solutions as well as stability to the variation of data. At the same time, computational methods of solving the multidimensional inverse dynamic problems are not as widely developed. The theory and methods of solving the direct and inverse problems for classical equations of elastodynamics have multiple applications in geophysics and engineering.

In the last part of the paper, practical applications of inverse problems are summarized. Moreover, a brief introduction into recent results of neural field theory is provided since inverse problems arising in cognitive sciences have become popular in an interdisciplinary community of cognitive neuroscience. As a final application basic setup of data assimilation as it is, for example, used in numerical weather prediction or for climate projections is investigated in details.