# An Adaptive Observer-Based Algorithm for Solving Inverse Source Problem for the Wave Equation

Authors: Sharefa Asiri; Chadia Zayane-Aissa; Taous-Meriem Laleg-Kirati
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Observers are well known in control theory. Originally designed to estimate the hidden states of dynamical systems given some measurements, the observers scope has been recently extended to the estimation of some unknowns, for systems governed by partial differential equations. In this paper, observers are used to solve inverse source problem for a one-dimensional wave equation. An adaptive observer is designed to estimate the state and source components for a fully discretized system. The effectiveness of the algorithm is emphasized in noise-free and noisy cases and an insight on the impact of measurements’ size and location is provided.


---

## Body

## 1. Introduction

In this paper we are interested in an inverse source problem for the wave equation. This problem appears frequently in many fields, especially, in modern seismology [1]. One important application of this problem is to distinguish between different types of seismic events (e.g., earthquake, implosion, or explosion) [2]. It is also important in monitoring hydraulic fracturing by which fractures are created in rocks such that entrapped hydrocarbons can be released and extracted [3].Inverse problems are usually solved using optimization techniques, where an appropriate cost function is minimized. However the ill-posedness of such problems generates instability. Regularization techniques are then used to restore the stability. Among the regularization techniques, Tikhonov regularization [4] is probably the most used one. For instance, it has been applied to the wave equation in [5, 6]. Other techniques have been also proposed to solve inverse problems; for example, in [7], a new minimization algorithm has been proposed to solve an inverse problem for the wave equation with unknown wave speed. Most of the proposed methods end up with an optimization step which generally turns to be computationally heavy, especially in the case of large number of unknowns, and may require an extensive storage.The objective of this paper is to present an alternative algorithm, based on observers, to solve the inverse source problem for the wave equation. Observers are well known in control theory for state estimation in finite dimensional dynamical systems. Presenting the distinctive feature and main advantage of operating recursively on direct problems, observers are gaining more and more interest in a wide variety of problems, including partial differential equations (PDEs) systems. For instance, in [8] states and parameters are estimated using an observer depending on a discretized space for a mechanical system. In [9], the initial state of a distributed parameter system has been estimated using two observers: one for the forward time and the other for the backward time. A similar approach has been used in [10], using the forward-backward approach to solve inverse source problem for the wave equation. An adaptive observer was applied in [11] for parameter estimation and stabilization of one-dimensional wave equation where the boundary observation suffers from an unknown constant disturbance. A similar work was proposed in [12] with the state as unknown and the boundary observation suffers from an arbitrary long time delay.Dealing with PDEs, either with observers or classical inverse problems methods, poses the challenge of approximating infinite dimensional systems. As regards observers, we can distinguish three approaches for studying such systems. The first approach considers the design of the observer in the continuous domain which requires mathematical analysis [13], and the application of the observer to real application will require some adaptation. The second approach consists in the semidiscretization of the equation in space. The result of this semidiscretization can be usually written in the standard state-space representation in the continuous domain (in time) which makes the extension of the known methods in control theory easier. The third approach is the full-discretization of the PDE in space and time. In this case we can write the system in a discrete state-space representation. We have chosen this latter approach since it is more suitable for real implementation. We show that it can give good results provided that some conditions, aimed at minimizing the effect of numerical errors resulting from discretization, are met.Another challenge, related to solving inverse problems in general, arises when it comes to measurement constraints. Indeed, from a practical point of view, we usually do not have enough measurements to estimate all the unknowns. Dealing with this source of ill-posedness, means, in observers theory framework, satisfying the equivalent property of observability. Indeed, given the PDE system together with the measurements, we can test in a prior step whether the unknown variables can be estimated fully or partially, regardless of the kind of observer to be used. For instance, in [8, 9, 11, 12], the measurements were taken as the time derivative of the solution of the wave equation. This kind of measurements gives a typical observability condition which has a positive effect on the stabilization, but it is less readily available than field measurements. Hence, some authors sought to solve inverse problems for wave equation using observers based on partial filed measurements, that is, measurements taken from the solution of the wave equation, as in [14–16].In this paper, we consider a fully discretized version of a one-dimensional wave equation and we propose a new algorithm for inverse source problem based on adaptive observer for the joint estimation of the states and the source term from partial measurements of the field. Adaptive observers are widely used in control theory for parameter estimation in adaptive control or fault estimation in fault detection and isolation [17, 18]. In Section 2, problem statement is detailed. Then, the observer design is presented in Section 3. Finally, numerical results are presented and discussed.


## 2. Problem Statement

Consider the one-dimensional wave equation with Dirichlet boundary conditions defining in the domainΩ
∶
=
(
x
,
t
)
∈
(
0
,
l
)
×
(
0
,
T
]:(1)
u
t
t
x
,
t
-
c
2
u
x
x
x
,
t
=
f
x
,
u
0
,
t
=
0
,
u
l
,
t
=
0
,
u
x
,
0
=
r
1
x
,
u
t
x
,
0
=
r
2
x
,where x is the space coordinate, t is the time coordinate, r
1
(
x
) and r
2
(
x
) are the initial conditions in L
2
[
0
,
l
], f
(
x
)
∈
L
2
[
0
,
l
] is the source function which is assumed, for simplicity, to be independent of time, and c is the velocity which is known. The notations u
a and u
a
a refer to the first and second derivatives of u with respect to a, respectively.Our inverse problem falls in the estimation of the sourcef
(
x
) in (1) using an adaptive observer with partial measurements of the field u available. We first propose to rewrite (1) in a system of first order PDEs by introducing two auxiliary variables v
(
x
,
t
)
=
u
(
x
,
t
) and w
(
x
,
t
)
=
u
t
(
x
,
t
) and let(2)
ξ
x
,
t
=
v
x
,
t
,
w
x
,
t
t
r
,where t
r refers to transpose. Then (1) can be written as follows:(3)
∂
ξ
x
,
t
∂
t
=
A
ξ
x
,
t
+
F
,
v
0
,
t
=
0
,
v
l
,
t
=
0
,
v
x
,
0
=
r
1
x
,
v
t
x
,
0
=
r
2
x
,
z
=
H
ξ
x
,
t
,where the operator A is given by A
=
0
I
c
2
(
∂
2
/
∂
x
2
)
0, F
=
0
f, z is the output, and H is the observation operator such that H
=
[
H
0
0
], where H
0 is a restriction operator on the measured domain.Discretizing system (3) using implicit Euler scheme in time and central finite difference discretization for the space gives the following discrete state-space representation:(4)
ξ
j
+
1
=
G
ξ
j
+
B
f
j
+
b
,
z
j
=
H
ξ
j
,
f
j
+
1
=
f
j
,
j
=
1,2
,
…
,
N
t
,where(5)
G
=
Δ
t
E
+
I
Δ
t
I
E
I
;
E
=
c
2
Δ
t
Δ
x
2
-
2
1
1
-
2
⋱
⋱
⋱
1
1
-
2
;
B
=
Δ
t
2
I
Δ
t
I
;
H
=
H
m
0
;and H
m
=
0
⋯
0
⋮
I
m
⋮
0
⋯
0
;
I
m is the identity matrix of dimension m, where m refers to the number of measurements, and b is a term that includes the boundary conditions such that(6)
b
=
c
2
Δ
t
2
Δ
x
2
v
1
j
0
1
×
N
x
-
2
c
2
Δ
t
2
Δ
x
2
v
N
x
j
c
2
Δ
t
Δ
x
2
v
1
j
0
1
×
N
x
-
2
c
2
Δ
t
Δ
x
2
v
N
x
j
t
r
.This system is linear multiple-input multiple-output discrete time invariant. If N
x refers to the space grid size and m refers to the number of measurements, then the state matrix G is of dimension 2
N
x
×
2
N
x, the observer matrix H is of dimension m
×
2
N
x, and the input matrix B is of dimension 2
N
x
×
N
x.The numerical scheme (4) is consistent. In addition, it is stable if and only if c
Δ
t
≤
Δ
x (the CFL condition). Thus, if c
Δ
t
≤
Δ
x, scheme (4) converges as Δ
t
→
0, Δ
x
→
0 to (3) and therefore to (1).


## 3. Observer Design

We propose to use an adaptive observer for the joint estimation of the statesv and w and the source f. This observer has been proposed in [17] and it has been developed for joint estimation of the state and the parameters. However, we propose to generalize the idea behind this observer to estimate the input considering each spatial sample of the input as an independent parameter. The adaptive observer is given by the following system of equations:(7)
z
^
j
=
H
ξ
^
j
,
Υ
j
+
1
=
G
-
L
H
Υ
j
+
B
,
f
^
j
+
1
=
f
^
j
+
σ
Υ
j
T
H
T
z
j
-
z
^
j
,
ξ
^
j
+
1
=
G
ξ
^
j
+
B
f
^
j
+
b
+
L
z
j
-
z
^
j
+
Υ
j
+
1
f
^
j
+
1
-
f
^
j
,where L is the observer gain matrix of dimension 2
N
x
×
m, ξ
^
j and f
^
j are the state and source estimates, respectively, Υ
j is a matrix sequence obtained by linearly filtering B, and σ is a scalar gain satisfying the following assumption as in [17].Assumption 1.
The scalar gainσ satisfies the following: (1)
σ
H
Υ
j
2
≤
1;
(2)
(
1
/
κ
)
∑
i
=
j
j
+
κ
-
1
σ
Υ
i
T
H
T
H
Υ
i
≥
β
I for some constant β
>
0, integer κ
>
0, and all j.Remark 2.
In the proposed method, no particular form for the matricesG, B, and H is required. However, all these matrices are assumed bounded. In our problem, wave equation with constant velocity, these matrices are actually constant and therefore always bounded.Remark 3.
From Remark2, any consistent and stable numerical method can be used to discretize system (3) provided that it ends up with bounded matrices (G, B, and H).Under Assumption1, Algorithm (7) converges exponentially fast when j tends to infinity in noise-free case, and the estimation errors remain bounded in the noisy case as long as the noises are bounded. Moreover, the estimation errors converge in the mean to zero if the noises have zero means; see Theorems 1 and 2 with their proofs in [17].


## 4. Numerical Simulations

To test the performance of the observer, we generated a set of synthetic data using the following parameters:Δ
x
=
0.01, l
=
2, Δ
t
=
0.01, and T
=
100. Thus, N
x
=
201 and N
t
=
10001. The velocity is chosen to be c
2
=
0.9, and the source is equal to f
(
x
)
=
3
sin
⁡
(
5
x
). The matrix sequence Υ
j and the scalar gain σ are chosen such that Assumption 1 is satisfied. The algorithm was implemented in Matlab and the tests were run for two main cases: noise-free and noisy datasets. In the noise-corrupted case, zero mean white Gaussian random noises were added to the states and to the measurements with standard deviations σ
ξ
=
0.007816 and σ
z
=
0.01044, respectively. The gain matrix is selected to have fast and accurate convergence of the observer. We took advantage of the particular structure of G to design the gain L. Indeed the matrix G is sparse, so we selected L to be also a sparse matrix. The number of unknown entries is then reduced and we identified them such that the eigenvalues of (
G
-
L
H
) are inside the unit circle. In general, standard pole placement can be used to select the gain matrix L.Figure1 shows the error in the estimated state and Figure 2 presents the exact and the estimated source; both figures exhibit noise-free and noisy cases with respect to full and partial measurements. For the partial measurements, we supposed that the field is available on half of the space domain only. Tables 1 and 2 show the minimum square error (MSE) in the estimated source in noise-free and noisy cases, respectively. Both tables show the error in case of full measurements, partial measurements taken from the middle, and partial measurements taken from the end.Table 1
Source estimation errors in the noise-free case.
Measurements
M
S
E
=
1
N
x
∑
i
=
1
N
x
f
i
-
f
^
i
Full
1.8168
×
1
0
-
14
Partial (middle)
0.3354
Partial (end)
0.2096Table 2
Source estimation errors in the noisy case.
Measurements
M
S
E
=
1
N
x
∑
i
=
1
N
x
f
i
-
f
^
i
Full
0.2865
Partial (middle)
0.4014
Partial (end)
0.3213Figure 1
The state error (ξ
-
ξ
^): (a) and (b) present the noise-free case with respect to full measurements and partial measurements, respectively. (c) and (d) show the noise-corrupted case with respect to full measurements and partial measurements, respectively. In the partial measurements cases, 50
% of the state components are taken from the end.Figure 2
The exact sourcef (blue) and the estimated source f
^ (black): (a) and (b) present the noise-free case with respect to full measurements and partial measurements, respectively. (c) and (d) show the noise-corrupted case with respect to full measurements and partial measurements, respectively. In the partial measurements cases, 50
% of the state components are taken from the end.In the noise-free case, the adaptive observer used in this paper provides a good estimate of the unknown source for the wave equation both when the field is available on the whole space domain and when it is available only on half the domain. However, we noticed in the second case a small error at the end of the interval. In the noisy case also, the reconstruction is good but can be improved by a good choice of the gainL.


## 5. Comparison between Observer-Based and Tikhonov-Based Approaches

To assess the observer performance in the source estimation, a comparison to optimization-based method has been performed. The initial problem has been first written in a suitable input/output formulation allowing us to have comparable problems. Then, due to the ill-posedness of the problem, a Tikhonov regularization has been used when solving the optimization problem. The two methods, the observer-based method and the optimization-based method, have been used for the estimation of the source for the same set of measurements and the same level of noise.To formulate the initial problem, we propose to derive the state and output at timek
+
p, p
∈
N, from the state at time k and the input sequence, using the state-space matrices G, B, and H. Thus, by repeating substitution from (4) and for some p
∈
N, we obtain a new state-space representation, where the transmission matrix is given by a Hankel matrix, as follows:(8)
ξ
j
+
p
=
G
p
ξ
j
+
C
p
f
p
j
+
b
,
z
p
j
=
O
ξ
j
+
τ
f
p
j
,where(9)
f
p
j
=
f
j
f
j
+
1
⋯
f
j
+
p
-
1
t
r
;
z
p
j
=
z
j
z
j
+
1
⋯
z
j
+
p
-
1
t
r
;
C
p
=
G
p
-
1
B
⋯
G
B
B
;
O
=
H
H
G
⋯
H
G
p
-
1
t
r
;
τ
=
0
m
×
N
x
0
0
⋯
0
H
B
0
m
×
N
x
⋱
⋱
⋮
H
G
B
H
B
0
m
×
N
x
⋱
0
⋮
⋱
⋱
⋱
0
H
G
p
-
2
B
⋯
H
G
B
H
B
0
m
×
N
x
; and b
=
1
p
⊗
b (⊗ is the Kronecker product).Thus, from the second equation in (8), new set of measurements can be defined as(10)
z
~
p
j
=
τ
f
p
j
,where z
~
p
j
=
z
p
j
-
O
ξ
j.The aim is to estimate the sourcef at time step j by minimizing the following cost function where Tikhonov regularization is used:(11)
J
α
f
p
j
=
1
2
τ
f
p
j
-
z
~
j
2
2
+
α
2
f
p
j
2
2
,where α is the regularization parameter. There are well-studied approaches for selecting this parameter such as L-curve, GCV, and NCP [19].For the numerical simulations, it is important to note that the size of the Hankel matrix depends on the space stepΔ
x, time step Δ
t, and the final time T. To have reasonable size which allows computation using Matlab, the values of these parameters have been chosen as follows: Δ
x
=
0.1, Δ
t
=
0.05, and T
=
2, respectively. This decrease in the final time T will affect the estimation errors convergence as discussed in Section 3, especially in the noisy case. Consequently we consider a small modification on the observer’s structure (7) in order to increase the robustness of the algorithm. This modification has been inspired by sliding mode observers and consists in adding t
a
n
h to the correction term as described in the following:(12)
z
^
j
=
H
ξ
^
j
,
Υ
j
+
1
=
G
-
L
H
Υ
j
+
B
,
f
^
j
+
1
=
f
^
j
+
Σ
Υ
j
T
H
T
z
j
-
z
^
j
+
γ
1
tanh
⁡
γ
2
z
j
-
z
^
j
,
ξ
^
j
+
1
=
G
ξ
^
j
+
B
f
^
j
+
b
+
L
z
j
-
z
^
j
+
γ
3
tanh
⁡
γ
4
z
j
-
z
^
j
+
Υ
j
+
1
f
^
j
+
1
-
f
^
j
,where γ
1, γ
2, γ
3, and γ
4 are scalers.The results for source estimation using observer-based and Tikhonov-based methods with full and partial measurements in noise-free and noisy cases are depicted in Figure3. The corresponding MSE are presented in Tables 3 and 4 for noise-free and noise-corrupted cases, respectively.Table 3
MSE for the source estimation using observer and Tikhonov methods in the noise-free case.
Measurements
Observer
Tikhonov
Full
1.2212
×
1
0
-
5
9.9112
×
1
0
-
13
Partial
0.5033
1.5539Table 4
MSE for the source estimation using observer and Tikhonov methods in the noisy case.
Measurements
Observer
Tikhonov
Full
0.2060
0.2355
Partial
0.6805
1.5466Figure 3
The exact sourcef (solid blue line) and the estimated one f
^ using observer solution (red) and modified-Tikhonov (green): (a) and (b) present the noise-free case with respect to full measurements and partial measurements, respectively. (c) and (d) show the noise-corrupted case with respect to full measurements and partial measurements, respectively. In the partial measurements cases, 50
% of the state components are taken from the end.Under the described conditions, the observer approach gives comparable results, in some cases better results than the optimization-based methods.


## 6. Discussion

We have studied the effect of number of measurements on the convergence of the proposed observer. Obviously, increasing number of measurements means increasing information on the state, thus insuring the observability condition for all the states. However, for some applications, only few measurements can be available and the idea is to study the effect of this number on the convergence of the observer.The analysis of the error of estimation of the source with respect to the number of measurements shows that numerical issues may happen when we reduce the number of measurements below a threshold. These numerical problems come in fact from the ill-conditioning of the observability matrixW:(13)
W
=
H
H
G
H
G
2
⋯
H
G
n
-
1
t
r
. The decay of the condition number of the observability matrix W as a function of the number of measurements is illustrated in Figure 4. It is well known in control theory that the rank of W gives the number of observable states. It is known also that a high condition number for the observability matrix leads to nearly unobservable states [20].Figure 4
Number of measurements versus the condition number of the observability matrixW.It is also important to study the effect of discretization on the performance of the method, which was not included in this paper as the objective was to assess the possibility of using this method and to discuss its performance in presence of noise. The scheme that we used in the paper works well but the step discretization may affect the performance and especially in case of few measurements. More investigations on this question are required.The objective of this paper was to propose a new method, as an alternative to the standard optimization methods, in order to solve inverse source problems for the wave equation. Of course this problem has several important applications in different fields ranging from geophysics to medical field and especially when few measurements are available. The first results obtained on simulations are promising and the observer approach seems to be suitable for real online estimation problems thanks to its recursive structure. However, we still have to investigate more the approach before we can claim its performance on real application. One of the points to assess carefully is the number of available measurements. Through this work, we studied the effect of measurements on the performance and from the comparison to the optimization-based methods in Section5 it is clear that the adaptive observer gives interesting results which is promising for real applications. The second important point is the effect of noise. Even if we succeed to obtain good results in noisy cases, we believe that some improvement can be suggested leading to some modifications of the observer structure aiming at improving the robustness properties.


## 7. Conclusion

In this paper, an adaptive observer for the joint estimation of the source and the states in the wave equation has been designed. Numerical simulations for the source and states estimation using observer have been presented, and they have proven the capability of observer to estimate both the source and the states in noise-free and noisy cases. A comparison between observer algorithm and an optimization-based method has been performed. This comparison considered also the different cases of noise (noise-free and noise-corrupted) with full and partial measurements. The results show the outperformance of the observer-based approach.



---

*Source: 796539-2015-10-19.xml*