# Identification of the Point Sources in Some Stochastic Wave Equations

Authors: Yaozhong Hu; Guanglin Rang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We introduce and study a type of (one-dimensional) wave equations with noisy point sources. We first study the existence and uniqueness problem of the equations. Then, we assume that the locations of point sources are unknown but we can observe the solution at some other location
continuously in time. We propose an estimator to identify the point source locations and prove the convergence of our estimator.


---

## Body

## 1. Introduction

Assume that there are a certain number of objects in a certain area of ocean or other media. The total number of the objects and the location of each object are unknown. We need to identify the total number and the precise locations of the objects. The objects are also assumed to emit (point source) sound waves and we are able to measure these sound waves received in some known locations. The objective is to use these measurements for our identification problem.This type of problem has been studied mathematically in the framework of inverse problems for partial differential equations (wave equations). The sound travels according to the following second-order wave equation with point sources:(1)∂2u(t,x)∂t2-∑i=1d‍∂2u(t,x)∂xi2=∑k=1N‍σkδ(x-yk)u(0,x)=0,∂u∂t(0,x)=0,
where x∈ℝd, t∈(0,T), N∈ℕ, and y1,…,yN are some given points in ℝd, δ(x) is the Dirac delta function, and σk, k=1,2,…,N, are some known constants. The solution (u(t,x),x∈D,0≤t≤T) is supposed to be known for some space points D and for some interval [0,T]. The total number N and each location yk of point sources are estimated from (u(t,x),x∈D,0≤t≤T). We refer to, for example, [1] and in particular the references therein for some recent study in this area. This theory has found substantial applications in determining the heat sources in heat conduction, the magnetic sources in brain, the earthquake sources of seismic waves, and so on.In practice, the sound wave travels inevitably under some influence of noises. Voluntarily or involuntarily, the point sources themselves may also emit noises to avoid being detected. Thus, we are led to the following stochastic wave equations:(2)∂2u(t,x)∂t2-∑i=1d‍∂2u(t,x)∂xi2=b(t,x,u(t,x))+σ0(t,x,u)B˙0(t,x)+∑k=1N‍σk(t,x,u)δ(x-yk)B˙k(t,x)u(0,x)=v0(x),∂u∂t(0,x)=v-0(x),
where x∈ℝd, t∈(0,T), N∈ℕ, and y1,…,yN are some given points in ℝd, δ(x) is the Dirac delta function, v0 and v1(x) are two given deterministic functions, and B˙k, k=0,1,2,…,N, are independent Gaussian noises which are white in time and correlated in space.Whenσ1=⋯=σk=0, the stochastic wave equation (2) has been studied since long time. Let us mention the first lecture note [2] and a recent lecture note [3]. Many properties such as the sample path Hölder continuity of the solution are obtained (see [4, 5] and the references therein).However, whenσi, i=1,…,N, are not all zero, then (2) is highly singular because of the presence of the Dirac delta functions multiplied by the Gaussian noises. Such equation has not been studied yet. The first objective of this paper is to give the definition of the solution to such an equation and to show the existence of uniqueness of the solution under some appropriate conditions. This will be done in Section 2. Since the case when σ1=⋯=σN=0 has been well studied, we will now assume σ0=0 to simplify our presentation. However, since our objective is the identification of the N and point source positions y1,…,yN, we will not go to spend too much effort here. For this reason, we restrict ourselves to only one space dimension case. We will present higher space dimension case in another project.To well explain our approach of identification, we will further restrict our model. We will consider the special case of (2) where the coefficients σk(t,x,u)=ρk(t) are independent of x and u and v0=v-0=0. Namely, we will concentrate on the stochastic wave equation of the form
(3)∂2u(t,x)∂t2-∑i=1d‍∂2u(t,x)∂xi2=∑k=1N‍ρk(t)δ(x-yk)B˙k(t,x)u(0,x)=0,∂u∂t(0,x)=0.
In this case, we will write down the explicit expression of the solution. It will be done in Section 3. In this section, we also obtain some properties of the solution which will be useful in the later section of the paper.Now, we assume that, in (3), the total number of point sources N and the positions yi, i=1,…,N, are unknown. However, we are able to observe the sound signal received at some given known locations x1,…,xm, continuously in the time interval [0,T]. Namely, we assume that (u(t,x1),…,u(t,xm),0≤t≤T) are known. We would like to use this information to identify N and yi, i=1,…,N. In Section 4, we will develop a new approach to obtain some statistical estimators N^ and y^1,y^2,…,y^N to estimate the total number N and the locations y1,y2,…,yN of the point sources. The approach combines the reciprocity gap functional approach from the theory of partial differential equations with theory from stochastic processes. We show the almost sure convergence of our estimators N^ and y^1,y^2,…,y^N to the true parameters N and y1,y2,…,yN.


## 2. Stochastic Wave Equations with Noisy Point Sources

Let(Ω,ℱ,P) be a basic probability space with right continuous filtration of σ-algebras {ℱt,t⩾0} satisfying the usual conditions. Let B(t,x)=(B1(t,x),B2(t,x),…,BN(t,x)),t⩾0,x∈ℝd be N dimensional Gaussian random fields. The formal derivatives B˙k(t,x)∶=(∂d+1/∂t∂x1⋯∂xd)Bk(t,x), k=1,2,…,N, are called the Gaussian noises. We assume that these Gaussian noises are white in time and correlated in space with covariance r(x,y) (see [5]). Namely, we assume that
(4)𝔼[B˙k(t,x)B˙k(s,y)]=δ(t-s)rk(x,y),
where 𝔼 denotes the expectation on (Ω,ℱ,P), rk:ℝd×ℝd→ℝ, and (k=1,2,…,N) are some symmetric positive definite functions of x and y. This is interpreted as
(5)𝔼[Bk(t,x)Bk(s,y)]=(t∧s)∫0x1‍⋯∫0xd‍∫0y1‍⋯∫0yd‍rk(u,v)dudv.
For any (deterministic) smooth function φ in 𝒞0∞([0,T]×ℝd), the stochastic integral Bk(φ)=∫0T‍∫ℝd‍φ(s,x)Bk(ds,dx) is well-defined in the sense of Walsh ([2]). The following fact is well-known: for any two smooth functions φ and ψ in 𝒞0∞([0,T]×ℝd), we have
(6)𝔼(Bk(φ)Bk(ψ))=∫0T‍ds∬ℝdφ(s,x)rk(x,y)ψ(s,y)dxdy.
We also call 𝔼(Bk(φ)Bk(ψ)) the covariance functional of B˙k, denoted by Jk(φ,ψ).Sometimes, we also use the Fourier transform theory to study the stochastic integral and the stochastic equations. We useℱφ(s,·)(ξ)=∫ℝd‍e2πixξφ(s,x)dx to denote the Fourier transformation of ϕ(s,·) and φ~(s,x)=φ(s,-x).If we assume that the noise is spatially homogeneous, that is,rk(x,y)=rk(x-y), then there exists a nonnegative tempered measure μk which is the Fourier transform of rk(x)dx. With this notation, we can also write
(7)Jk(φ,ψ)=∫0∞‍ds∫ℝd‍μk(dξ)ℱφ(s,·)(ξ)ℱψ~(s,·)(ξ).
From the general theory of stochastic integral (see, e.g., [3]), we see that if (g(t,x),0≤t≤T,x∈ℝd) is a real valued ℱt-adapted process such that
(8)∫0Tds∬ℝdrk(x,y)𝔼[|g(s,x)||g(s,y)|]dxdy<∞,
then the stochastic integral ∫0T‍∫ℝd‍g(t,x)Bk(dt,dx) is well-defined and
(9)𝔼[∫0T‍∫ℝd‍g(s,x)Bk(ds,dx)]2=∫0T‍ds∬ℝdrk(x,y)𝔼[g(s,x)g(s,y)]dxdy.
Now, let y be a given fixed point in ℝd and let (λ(t,x),0≤t≤T,x∈ℝd) be a real valued ℱt-adapted process. We want to define the stochastic integral ∫0t‍∫ℝd‍λ(s,x)δ(x-y)F(dsdx), where δ(x-y) is the Dirac delta function. To this end, we will use smooth functions to approximate the Dirac delta function.Letϕ be a smooth function with compact support K=[-1,1]d. Set ϕn(x-y)=ndϕ(n(x-y)). It is clear that ϕn converges to δ(x-y) in the sense of distribution as n→∞. For each n, the stochastic integral
(10)In∶=∫0T‍∫ℝd‍λ(s,x)ϕn(x-y)Bk(ds,dx)
is well-defined. We want to know whether In has a limit in L2(Ω,ℱ,P) or not. First, we have the following computations:
(11)𝔼[In-Im]2=𝔼(∫0T‍∫ℝd‍λ(s,x)[ϕn(x-y)kkkkkkkjjjjjjjjjjjjjjjkklllk-ϕm(x-y)]Bk(ds,dx)∫222)2=∫0T‍∬ℝd𝔼[λ(s,x)λ(s,z)]kkkkkkkkkk×(ϕn(x-y)-ϕm(x-y))kkkkkkkkkk×(ϕn(z-y)-ϕm(z-y))kkkkkkkkkk×rk(x,z)dxdzds.
Passing the difference ϕn(x-y)-ϕm(x-y) to that of λ and r, we have
(12)𝔼[In-Im]2=∫0T‍∬ℝd𝔼[λ(s,y+xn)λ(s,y+zn)kkkkkkkkkkkkk×r(y+xn,y+zn)-λ(s,y+xn)kkkkkkkkkkkkk×λ(s,y+zm)r(y+xn,y+zm)]llllllllllllllllllll×ϕ(x)ϕ(z)dxdzds+∫0T‍∬ℝd𝔼[λ(s,y+xm)λ(s,y+zm)kkkkkkkkkkkkkkk×r(y+xm,y+zm)-λ(s,y+xm)kkkkkkkkkkkkkkk×λ(s,y+zn)r(y+xm,y+zn)]kkkkkkkkkkkk×ϕ(x)ϕ(z)dxdzds≕I1+I2.I1 can be written as
(13)I1=∫0T‍∬ℝd𝔼{λ(s,y+xn)dddddddddkk×(λ(s,y+zn)kkkkkkkkkkkkkkk×[r(y+xn,y+zn)kkkkkkkkkkkkkkkkkk-r(y+xn,y+zm)]kkkkkkkkkkk+r(y+xn,y+zm)kkkkkkkkkkk×[λ(s,y+zn)kkkkkkkkkkkkkk-λ(s,y+zm)])}kkkkkkk×ϕ(x)ϕ(z)dxdzds.
Now, we let r(x,y) be a continuous function of x and y. We also assume that there is a ε0>0 such that
(14)sup⁡|y-x|≤ε0𝔼∫0T‍[λ(s,y)2]ds<∞,lim⁡ε→0sup⁡|y-x|≤ε𝔼∫0T‍[[λ(s,y)-λ(s,x)]2]ds=0.
We write I1=I11+I12, where I11 and I12 are defined and estimated as follows. For I11, we have
(15)I11∶=∫0T‍∬ℝd𝔼{λ(s,y+xn)λ(s,y+zn)}kkkkkkkkk×[r(y+xn,y+zn)-r(y+xn,y+zm)]kkkkkkkkkkk×ϕ(x)ϕ(z)dxdzds≤Csup⁡|x|≤1,|z|≤1|r(y+xn,y+zn)-r(y+xn,y+zm)|×∫0T‍∬ℝd𝔼[λ(s,y+xn)]2ϕ(x)ϕ(z)dsdxdz≤Csup⁡|x|≤1,|z|≤1|r(y+xn,y+zn)-r(y+xn,y+zm)|×∫ℝd‍∫0T‍𝔼[λ(s,y+xn)]2dsϕ(x)dx≤Csup⁡|x|≤1,|z|≤1|r(y+xn,y+zn)-r(y+xn,y+zm)|×sup⁡|x|≤1∫0T‍𝔼[λ(s,y+xn)]2ds⟶0.
For I12, we have (recall K=[-1,1]d)
(16)I12∶=∫0T‍∬ℝd𝔼{λ(s,y+xn)kkkkkkkkkkk×r(y+xn,y+zm)kkkkkkkkkkkkkkkkkkkkkk×[λ(s,y+zn)kkkkkkkkkkkkkkkk-λ(s,y+zm)]}kkkkkkkkkkkk×ϕ(x)ϕ(z)dxdzds≤C∫0T‍∬K𝔼{λ(s,y+xn)kkkkkkkkkkkk×[λ(s,y+zn)kkkkkkkkkkkkkkk-λ(s,y+zm)]}dxdzds≤C{∫0T‍∫K‍𝔼λ2(s,y+xn)dxds}1/2×{∫0T‍∫K‍𝔼[λ(s,y+zn)kkkkkkkkkkk-λ(s,y+zm)]2dzds}1/2≤C{sup⁡|x|≤1∫0T‍𝔼λ2(s,y+xn)dx}1/2ds×{sup⁡|x|≤1∫0T‍𝔼[λ(s,y+zn)kkkkkkkkkkkkk-λ(s,y+zm)]2ds}1/2⟶0.
Thus, we have I1 converging to 0.In the same way, we can show thatI2 converges to 0 under the same conditions.Theorem 1.
LetB˙k(t,x) be a Gaussian noise which is white in time and correlated in space with covariance rk(x,y). Assume that r(x,y) is a continuous function of x and y. Let (λ(t,x),0≤t≤T,x∈ℝd) be an ℱt-adapted processes such that conditions (14) hold. Then, the stochastic integral ∫0T‍∫ℝd‍λ(t,x)δ(x-y)Bk(ds,dx) exists and
(17)𝔼[∫0T‍∫ℝd‍λ(s,x)δ(x-y)Bk(ds,dx)]2=rk(y,y)∫0T‍𝔼[λ(s,y)]2ds.Proof.
The above argument shows the existence of the stochastic integral∫0T‍∫ℝd‍λ(t,x)δ(x-y)Bk(ds,dx). We have
(18)𝔼[In]2=∫0T‍∬ℝd𝔼[λ(s,x)λ(s,z)]kkkkk×ϕn(x-y)ϕn(z-y)kkkkkkkkkk×rk(x,z)dxdzds=∫0T‍∬ℝd𝔼[λ(s,y+xn)λ(s,y+zn)kkkkkkkkkkkk×rk(y+xn,y+zn)]kkkkkkkkk×ϕ(x)ϕ(z)dxdzds.
Now, the Fatou lemma yields (17).We also need to bound general moments of the stochastic integral∫0T‍∫ℝd‍λ(t,x)δ(x-y)Bk(ds,dx). We have the following Burkholder type inequality.Theorem 2.
LetB˙k(t,x) be a Gaussian noise which is white in time and correlated in space with covariance rk(x,y). Assume that r(x,y) is a continuous function of x and y. Let (λ(t,x),0≤t≤T,x∈ℝd) be an ℱt-adapted processes such that
(19)sup⁡|y-x|≤ε0𝔼(∫0T‍[λ(s,y)2]ds)p/2<∞,lim⁡ε→0sup⁡|y-x|≤ε𝔼(∫0T‍[[λ(s,y)-λ(s,x)]2]ds)p/2=0.
Then, the stochastic integral ∫0T‍∫ℝd‍λ(t,x)δ(x-y)Bk(ds,dx) exists and
(20)𝔼[∫0T‍∫ℝd‍λ(s,x)δ(x-y)Bk(ds,dx)]p≤Cprkp/2(y,y)𝔼(∫0T‍[λ(s,y)]2ds)p/2.Now, we turn to consider the existence and uniqueness of the solution to the stochastic wave equation (2). We will follow the idea of mild solution. Since the Green function is more sophisticated to study in the dimension higher than 1, we will only study the one space dimensional wave equation in this paper. Higher dimension case needs much more care. Consider
(21)∂2u(t,x)∂t2-∂2u(t,x)∂x2=b(t,x,u(t,x))+σ0(t,x,u(t,x))B˙0(t,x)+∑k=1N‍σk(t,x,u(t,x))×δ(x-yk)B˙k(t,x)u(0,x)=v0(x),∂u∂t(0,x)=v-0(x).
In the one-dimensional case, the associated fundamental solution (Green’s function) of the wave operator ∂2/∂t2-∂2/∂x is
(22)G(t,x,y):=12I{|x-y|≤t}.
First, we give the following definition about the solution.Definition 3.
A random field{u(t,x),t>0,x∈ℝ} is called a solution to the stochastic wave equation if it satisfies the following identity:
(23)u(t,x)=∑k=1N‍∫0t‍∫ℝ‍G(t-s,x-z)kkkkkkk×σk(s,z,u(s,z))δ(z-yk)Bk(ds,dz)+∫0t‍∫ℝ‍G(t-s,x-z)b(s,z,u(s,z))dsdz+ddt∫ℝ‍G(t,x-z)v0(z)dz+∫ℝ‍G(t,x-y)v-0(z)dz.Theorem 4.
Assume thatb(t,x,u) and σk(t,x,u), k=0,1,…,N, satisfy the global Lipschitz condition and the linear growth condition in u uniformly in t∈[0,T] and x∈ℝ. Assume that v0 is bounded continuous functions in ℝ and v-0 is bounded continuously differentiable functions in ℝ with bounded derivative. Then, there is a unique solution to (21).Proof.
Since the solution of wave equation has the past-light cone property (see [6], p. 63), we can study the solution on a bounded domain. To simplify the presentation, we assume b=0.
Let us define𝔹p as the set of all mappings u:[0,T]×ℝ×Ω→ℝ such that u(t,x,ω) is continuous in (t,x)∈[0,T]×ℝ for almost all ω∈Ω and
(24)sup⁡0≤t≤T,x∈ℝ𝔼|u(t,x)|p<∞.
It is clear that 𝔹p is a Banach space with the norm
(25)∥u∥𝔹p∶=sup⁡0≤t≤T,x∈ℝ(𝔼|u(t,x)|p)1/p.
We will prove the existence and uniqueness of the solutions in 𝔹=∩p≥2𝔹p to (21) by the Picard iteration. We define {ui}i=0∞ recursively by the following
(26)ui(t,x)=T(ui-1)(t,x)∶=∑k=1N‍∫0t‍∫ℝ‍G(t-s,x-z)kkkkkkk×σk(s,z,ui-1(s,z))kkkkkkkkk×δ(z-yk)Bk(ds,dz)+ddt∫ℝ‍G(t,x-z)v0(z)dz+∫ℝ‍G(t,x-y)v-0(z)dz
for any (t,x)∈[0,T]×ℝ. We also define u0=0.
First, we show the well-posedness of the above stochastic integral for everyi=1,2,….
For anyp>1, let ui-1∈𝔹2p. By the Burkholder inequality, we have
(27)𝔼|ui(t,x+h)-ui(t,x)|2p⩽C∑k=1N‍rp(yk,yk)𝔼kkk×{∫0t‍[G(t-s,x+h-yk)kkkkkkkkkkkkkk-G(t-s,x-yk)]2kkkkkkkkk×|σk(s,yk,ui-1(s,yk))|2ds∫0t‍}p≕C∑k=1N‍I1k.
By the Hölder inequality, we have
(28)I1k⩽C[∫0t‍|G(t-s,x+h-yk)kkkkkkkk-G(t-s,x-yk)|qds∫0t‍]p/q×∫0t‍|G(t-s,x+h-yk)-G(t-s,x-yk)|pkkkkkk×𝔼|σk(s,yk,ui-1(s,yk))|2pds⩽C[∫0t‍|G(t-s,x+h-yk)kkkkkkk-G(t-s,x-yk)|qds∫0t‍]p/q×∫0t‍|G(t-s,x+h-yk)-G(t-s,x-yk)|pdskkkkkk×[1+sup⁡s,yk𝔼|ui-1(s,yk)|2p]⩽C|h|1+p/q,
where q=p/(p-1) is the conjugate number of p.
In the similar way, we can show𝔼|u(t,x)-u(t+s,x)|2p⩽C|s|1+p/q for any (t,x)∈[0,T]×ℝ. Thus by induction we have
(29)𝔼|ui(t+s,x+h)-ui(t,x)|2p⩽C(|s|p+|h|p),
for all i∈ℕ and (t,x)∈[0,T]×ℝ. Let 𝒞α denote the set of all functions of t∈[0,T] and x∈ℝ, which are Hölder continuous on both t and x on any compact subinterval of ℝ. Hence, by Kolmogorov’s theorem, we have, for any (fixed) 0<α<1/2, ui∈𝒞α for every i∈ℕ. It is easy to check 𝔼|u(t,x)|2p⩽C. Now, we verify easily
(30)∫0t‍𝔼[G(t-s,x-y)σ(s,y,ui(s,y))]2ds⩽C.
Moreover,
(31)∫0t‍𝔼[G(t-s,x-y)σ(s,y,ui(s,y))kkkk-G(t-s,x-z)σ(s,z(ui(s,z)))]2ds⩽∫0t‍𝔼[(s,z(ui(s,z)))2|G(t-s,x-y)kkkkkkkkkk-G(t-s,x-z)|2|σ(s,y,ui(s,y))|2kkkkkkkkkk+G2(t-s,x-z)kkkkkkkkkk×(σ(s,y,ui(s,y))kkkkkkkkkkσii×-σ(s,z,ui(s,z)))2]ds⩽C|y-z|
for any x∈ℝ. Thus, G(t-s,x-y)σ(s,y,ui(s,x)),i=1,2,…,N, satisfy conditions (14). Therefore, (26) is well-defined for all i∈ℕ.
Next, we show that{ui,i=1,2,…} is a Cauchy sequence in 𝔹p for any p. For n,m∈ℕ, by Burkholder’s inequality and Jensen’s inequality, we have
(32)𝔼|un(t,x)-um(t,x)|2p⩽C∑k=1N‍rp(yk,yk)𝔼{∫0t‍G2(t-s,x-yk)kkkkkkkkkkkkkkkkkkkk×(σk(s,yk,un-1(s,yk))kkkkkkkkkkkkkkkkkkkkkkk-σk(s,yk,um-1(s,yk)))2ds∫0t‍}p⩽C∑k=1N‍rp(yk,yk)×∫0t‍𝔼|un-1(s,yk)-um-1(s,yk)|2pds.
Let us denote by 𝔹p,t the set of all mappings u:[0,t]×ℝ×Ω→ℝ such that u(s,x,ω) is continuous in (s,x)∈[0,t]×ℝ for almost all ω∈Ω and
(33)sup⁡0≤t≤T,x∈ℝ𝔼|u(t,x)|p<∞.
We also denote
(34)∥u∥𝔹p,t∶=sup⁡0≤s≤t,x∈ℝ(𝔼|u(s,x)|p)1/p.
Thus, (32) can be written as
(35)∥un-um∥𝔹p,t2p2p≤C∫0t‍∥un-um∥𝔹p,s2p2pds.
Now, a routine argument shows that un is a Cauchy sequence in 𝔹2p. The limit of this sequence is denoted by u.
Lettingi→∞ in (29), we have
(36)𝔼|u(t+s,x+h)-u(t,x)|2p⩽C(|s|p+|h|p)
and sup⁡x∈ℝ𝔼|u(t,x)|2p⩽C. This implies that the stochastic integral ∫0t‍G(t-s,x-y)σk(s,y,(u(s,y)))δ(y-yk)Bk(ds,dy) is well-defined. It is easy to see
(37)sup⁡0≤t≤T,x∈ℝ𝔼|∫0t‍G(t-s,x-y)δ(y-yk)kkkkkkkkk×(σk(s,y,(ui(s,y)))-σk(s,y,(u(s,y))))kkkkkkkkk×Bk(ds,dy)∫0t‍|2p⟶0
as i→∞. Now, letting i tend to infinity on both sides of (26), we see that u satisfies (23). The uniqueness can be proved in similar way. Thus, we complete the proof of this theorem.


## 3. Some Preliminaries for Estimation of the Point Sources

To simplify the presentation of estimation method, we assumeσ0(t,x,u)=0 in (21). We also assume that the noises are space independent. Without loss of generality, we also assume v0(x)=v-0(x)=0. Moreover, we assume that σk(t,x,u)=ρk. This means we will consider the following stochastic wave equation:
(38)∂2u(t,x)∂t2-∂2u(t,x)∂x2=∑k=1N‍ρkδ(x-yk)B˙k(t)u(0,x)=0,∂u∂t(0,x)=0.From the Duhamel principle, we know that the solution of one-dimensional stochastic wave equations (38) is given by
(39)u(t,x)=∫0t‍∫R‍1{|x-y|≤t-s}(y,s)∑k=1N‍ρkδ(y-yk)dydBsk=∑k=1N‍∫0t‍∫R‍ρk1{|x-y|≤t-s}(y,s)δ(y-yk)dydBsi.In the remainder of this paper, we assume that the parametersN, y1,…,yNρ1,…,ρN are unknown. However, at some fixed location x∈ℝ, we can observe the sound wave signal {u(t,x), 0<t⩽T} continuously over the time interval [0,T]. We would like to use {u(t,x),0<t⩽T} to identify N, yk, and ρk,k=1,2,…,N. If one can observe the sound wave signals at some other locations x1,…,xm, we can use the similar approach to (better) estimate N, y1,…,yNρ1,…,ρN from all the observations {u(t,x1),…,u(t,xm),0<t⩽T}.Putαk=|yk-x|,k=1,2,…,N. We arrange the real positive numbers αk in increasing order. For example, we can assume α1<α2<⋯<αN.Ift<α1, then, for any k and s≤t, α1=|yk-x|≤t-s cannot hold, and hence u(t,x)=0.Ifα1≤t<α2 and 0<s<t-|y1-x|=t-α1, then y1∈{y,|y-x|≤t-s}; hence
(40)u(t,x)=∫0t-α1‍ρ1dBs1=ρ1Bt-α11.
Similarly, we have, for α2≤t<α3,
(41)u(t,x)=ρ1Bt-α11+ρ2Bt-α22;⋮
for αN-1≤t<αN,
(42)u(t,x)=ρ1Bt-α11+ρ2Bt-α22+⋯+ρN-1Bt-αN-1N-1;
and, for αN≤t,
(43)u(t,x)=ρ1Bt-α11+ρ2Bt-α22+⋯+ρNBt-αNN.Now, for any processZt, we define its quadratic variation process (if it exists) as
(44)[Z,Z]t=[Z]t=lim⁡|π|→0∑j=0n-1‍(Ztj+1-Ztj)2,
where π:0=t0<t1<⋯<tn=t is a partition of the interval [0,t] and |π|=max⁡0≤j≤n-1(tj+1-tj). For d dimensional Brownian motion B={B1,…,Bd}, we know that [Bi,Bj]t=tδij for i,j=1,2,…,d and any t>0, where δij=1, if i=j, 0, otherwise.The quadratic variation process ofu(t,x) is given by
(45)h(t)=Δ[u(·,x)]t={0,t<α1;ρ1(t-α1),α1≤t<α2;⋮⋮∑k=1N-1‍ρk(t-αk)αN-1≤t<αN;∑k=1N‍ρk(t-αk),αN≤t.
Since we can observe {u(t,x)} at the space location x continuously from the time interval [0,T], we know that the quadratic variation process h(t), 0≤t≤T is also an observable.LetT be large enough; for example, T>αN. We denote α0=0,αN+1=T. On [0,T], h(t) has its second-order distributional derivative by
(46)h′′(t)=∑k=1N‍ρkδ(t-αk).
For g∈C2[0,T] with g(0)=0, we can define linear operator ℛ(g), which is called thereciprocity gap functionalas follows:
(47)ℛ(g)=Δ∫0T‍h(t)g′′(t)dt-h(T)g′(T)+h′(T)g(T).
By the integration by parts formula, we have
(48)ℛ(g)=∑k=1N‍ρkg(αk).
For any function independent of the unknown parameters, we know from the definition that ℛ(g) is also independent of the unknown parameters. Namely, ℛ(g) is observable.To obtain our estimators for the parameters, we take(49)gn(t)=(sinπt2T)n,n=1,2,….
Then,
(50)ℛ(gn)=∑k=1N‍ρk(sinπαk2T)n=∑k=1N‍ρkSkn,
where Sk=sin(παk/2T).Clearly, since{ℛ(gn),n=1,2,…} are constructed from u(t,x), we see that they are observable. Furthermore, once we know Sk, then we can get αk from the identity αk=(2T/π)arcsin(Sk).TheSk can be obtained from ℛ(gn) by the following way. We let Hμ(μ=0,1) be the N×N matrices of the following forms:
(51)Hμ=(ℛ(gμ)ℛ(gμ+1)⋯ℛ(gμ+N-1)ℛ(gμ+1)ℛ(gμ+2)⋯ℛ(gμ+N)⋮⋮⋮⋮ℛ(gμ+N-1)ℛ(gμ+N)⋯ℛ(gμ+2N-2)).The following result is from [7, 8]. For more details, see [7, 8] and references therein.Theorem 5.
LetHμ(μ=0,1) be defined above; then S1,S2,…,SN are the eigenvalues of Hermite matrix H0-1H1.Proof.
First, we will introduce some intermediate quantities as follows. Forj∈N, we define
(52)Aj=((S1)j(S2)j⋯(SN)j(S1)j+1(S2)j+1⋯(SN)j+1⋮⋮⋯⋮(S1)N+j-1(S2)N+j-1⋯(SN)N+j-1)N×N.
We denote the diagonal matrix
(53)D=diag⁡(S1,S2,…,SN).
We define the vectors
(54)μj=(ℛ(gj)ℛ(gj+1)⋮ℛ(gN+j-1))N×1
and vector ρ=(ρ1,ρ2,…,ρN)T.
It is obvious that, for allj⩾0, we have
(55)μj=Ajρ,Aj+1=AjD.
Thus, for all j, we have
(56)μj+1=A0Dj+1ρ=A0DA0-1A0Djρ=A0DA0-1μj.
Denote Γ=A0DA0-1. Then, μj+1=Γμj and the matrix Γ has the Sk as its eigenvalue corresponding to the eigenvector βk=(1,Sk1,…,SkN-1)T for k=1,2,…,N.
On the other hand, we have(57)H1=(μ1,μ2,…,μN)=Γ(μ0,μ1,…,μN-1)=ΓH0.
Because A0 is a Vandermonde matrix and S1,S2,…,SN are assumed to be distinct, we also can show that vectors μ0,μ1,…,μN-1 are independent, which means that H0 is invertible. Hence, Γ=H1H0-1. The conclusion follows.Remark 6.
LetM>N be integers. Defining μj(M)=(ℛ(gj),ℛ(gj+1),…,ℛ(gM+j-1))T (T denotes the transpose), one can see that the vectors μ0(M),μ1(M),…,μM-1(M) are images of ρ,Dρ,…,DM-1ρ∈RN by the action of the matrix
(58)A0(M)=(11⋯1(S1)1(S2)1⋯(SN)1⋮⋮⋯⋮(S1)M-1(S2)M-1⋯(SN)M-1)M×N,
from which and from the independence of the vectors μ0(M),μ1(M),…,μM-1(M) one can deduce that the determinant of matrix (μ0,μ1,…,μM-1) is equal to zero. For more details, see [7].Remark 7.
One can find, from the proof, that vectorρ is the solution of linear equations μ0=A0ρ.


## 4. Estimations for Point Sources from Discrete Time Observations

By Theorem5, we know that the parameters we want to estimate are contained in the eigenvalues of Hermite matrix H0-1H1.We assume in this section that the wave signals are observed at the locationx∈ℝ but at discrete time instants 0=t0<t1<⋯<tm=T. We denote Δ=tj+1-tj(j=0,1,…,m-1).We define the approximationhΔ of the quadratic variations h of solution process u(t,x) for any Δ>0. Assume that there are integers l1,l2,…,lN such that tlj⩽αj<tlj+1 for all j=1,2,…,N. For any 0<t<T, since there is an integer jt such that tjt⩽t<tjt+1, one defines
(59)hΔ(t)=Δ∑k=1jt‍|u(tk,x)-u(tk-1,x)|2.
This is also an approximation of the quadratic variations process h(t) of solution process u(t,x) based on the observation time instants 0=t0<t1<⋯<tm=T. We know that for t fixed hΔ(t)→L2jt-α1-α2-⋯-αj  as  Δ→0 with αj⩽t<αj+1. In fact, we can show that convergence holds in the almost sure sense (see Lemma 10).To usehΔ(t) to obtain estimators of the parameters, we compute the following. For any twice differentiable function g, by Abel sum formula, we have
(60)∫0T‍hΔ(t)g′′(t)dt=hΔ(T)g′(T)-∫0T‍g′(t)dhΔ(t)=hΔ(T)g′(T)-∑i=1m‍g′(ti)|u(ti,x)-u(ti-1,x)|2=hΔ(T)g′(T)-|u(T,x)-u(tMΔ-1,x)|2×∑i=1m‍g′(ti)+∑i=1m-1‍(Zi+1-Zi)∑k=1i‍g′(tk),
where Zi=|u(ti,x)-u(ti-1,x)|2, i=1,2,…,m. Thus, we can define
(61)ℛΔ(g)=Δ∫0T‍hΔ(t)Δg(t)dt-hΔ(T)g′(T)kkk+|u(T,x)-u(tMΔ-1,x)|2∑i=1m‍g′(ti).
For any nonnegative integer μ, we introduce
(62)HΔ,μ,m=(ℛΔ(gμ)ℛΔ(gμ+1)⋯ℛΔ(gμ+m-1)ℛΔ(gμ+1)ℛΔ(gμ+2)⋯ℛΔ(gμ+m)⋮⋮⋮⋮ℛΔ(gμ+m-1)ℛΔ(gμ+m)⋯ℛΔ(gμ+2m-2)).First, we estimateN. It is known that when m>N, H0,m will be degenerate. However, there is a difference between HΔ,0,m and H0,m and there may be also error from the computer computation. To make the method robust, we introduce a small number δ. We propose to estimate N by the following estimator:
(63)N^=max⁡{m,|det⁡(HΔ,0,m)|>δ}.As in the case of continuous time observation we letSΔ,1↑,SΔ,2↑,…,SΔ,N↑ be the n eigenvalues (in increasing order) of the Hermite matrix ΓΔ=HΔ,0-1HΔ,1. Then, we obtain the estimations of the locations of the point sources by
(64)α^Δ,k=|y^k-x|=gk-1(SΔ,k↑),
for k=1,2,…,N. The estimation of yk follows from (64). Note that gk,k=1,2,…, are given by (49).Similar toA0,μ0 in Theorem 5, we define
(65)A0(Δ)=(11⋯1SΔ,1↑SΔ,2↑⋯SΔ,N↑⋮⋮⋯⋮(SΔ,1↑)N-1(SΔ,2↑)N-1⋯(SΔ,N↑)N-1)N×N,μ0(Δ)=(ℛΔ(g0)ℛΔ(g1)⋮ℛΔ(gN-1))N×1.
As in Remark 7, the estimation ρ^=(ρ^1,ρ^2,…,ρ^N) of the strength ρ of the point sources can be obtained by solving the following linear equation:
(66)μ0(Δ)=A0(Δ)ρ^.Now, we can summarize the estimation procedure as follows.Step 1.
ComputehΔ(t) on [0,T] and ℛΔ(gi) according to transformation (61).Step 2.
Identify the value ofN as the maximum M such that
(67)det⁡(μ0(Δ,M),μ1(Δ,M),…,μM-1(Δ,M))≠0.Step 3.
Compute the eigenvaluesSΔ,k↑ of matrix HΔ,0-1HΔ,1.Step 4.
Compute the intensities by formulaeρ(Δ)=A0-1(Δ)μ0(Δ).


## 5. Convergence of Estimations

In this section, we will show that the estimations obtained in the previous section converge to the true values a.s. as time spaceΔ tends to zero.Theorem 8.
For any functiong∈C2([0,T]) with g(0)=0, we have ℛΔ(g) converging to ℛ(g) a.s. as Δ→0.This will be fulfilled by the following two lemmas.Lemma 9.
For any functiong∈C([0,T]), ∫0T‍hΔ(t)g(t)dt converges to ∫0T‍h(t)g(t)dt in the sense of L2 as Δ→0.Proof.
It is sufficient to show that, for any functiong∈C([0,T]), ∫0T‍hΔ(t)g(t)dt converges to ∫0T‍h(t)g(t)dt in the sense of L2 as Δ→0. By Hölder inequality, one has, for any function g∈C([0,T]),
(68)𝔼(∫0T‍hΔ(t)g(t)dt-∫0T‍h(t)g(t)dt)2⩽∫0T‍𝔼(hΔ(t)-h(t))2dt×∫0T‍g2(t)dt=∑j=1N‍∫αjαj+1‍𝔼[hΔ(t)-h(t)]2dt×∫0T‍g2(t)dt=∫0T‍g2(t)dt∑j=1N‍Ij,
where
(69)Ij=∫αjtlj+1‍𝔼(hΔ(t)-h(t))2dt+∑k=lj+1lj+1-1‍∫tktk+1‍𝔼(hΔ(t)-h(t))2dt+∫tlj+1αj+1‍𝔼(hΔ(t)-h(t))2dt.
For lj+1⩽k⩽lj+1-1, we have
(70)∫tktk+1‍𝔼(hΔ(t)-h(t))2dt=∫tktk+1‍𝔼(∑i=1k‍|u(ti,x)-u(ti-1,x)|2kkkkkkkkkkkk-(∑i=1j‍ρi(t-αi)))2dt=∫tktk+1‍𝔼(∑s=1j-1‍∑i=ls+1ls+1‍|u(ti,x)-u(ti-1,x)|2kkkkkkkkkkk+∑i=lj+1k‍|u(ti,x)-u(ti-1,x)|2kkkkkkkkkkkk-(∑i=1j‍ρi(t-αi))∑s=1j-1‍)2dt=∫tktk+1‍𝔼(∑s=1j-1‍(∑i=ls+2ls+1‍|(ρ1Bti-α11+⋯+ρsBti-αss)kkkkkkkkkkkkkkkkk-(ρ1Bti-1-α11+⋯+ρsBti-1-αss)|2kkkkkkkkkkkkkkkkkk+|(ρ1Btls+1-α11+⋯+ρsBtls+1-αss)kkkkkkkkkkkkkkkkkk-(ρ1Btls-α11+⋯+ρs-1Btls-αs-1s-1)|2∑i=ls+2ls+1‍)kkkkkkkkkkkk+(∑i=lj+2k‍|(ρ1Bti-α11+⋯+ρjBti-αjj)kkkkkkkkkkkkkkkkkkkk-(ρ1Bti-1-α11+⋯+ρjBti-1-αjj)|2kkkkkkkkkkkkkkkkkkk+|(ρ1Btlj+1-α11+⋯+ρjBtlj+1-αjj)kkkkkkkkkkkkkkkkkkkk-(ρ1Btlj-α11+⋯kkkkkkkkkkkkkkkkkkkkkkkkk+ρj-1Btlj-αj-1j-1)|2∑i=lj+2k‍)kkkkkkkkkkk-(∑i=1j‍ρi(t-αi))∑s=1j-1‍)2dt.
Let X be a standard normal random variable. Denote the constant 𝔼(X2-1)2 by A (actually A=2). By the independent increments property of Brownian motion and the independence between Bk and Bl for k≠l, we have
(71)∫tktk+1‍𝔼(hΔ(t)-h(t))2dt=C∫tktk+1‍(A∑s=1j-1‍(((s-1)Δ+tls+1-αs)2(ls+1-ls-1)s2Δ2kkkkkkkkkkkkkkk+((s-1)Δ+tls+1-αs)2)kkkkkkkkkkkkkk+Aj2Δ2(k-lj-1)kkkkkkkkkkkkkk+A((j-1)Δ+tlj+1-αj)2A∑s=1j-1‍)kkkkkk+[(α1+⋯+tlj+1-αj)(jt-α1-⋯-αj)kkkkkkkkkk-((α1+⋯+tlj+1-αj)jkΔ-l1Δ-⋯-ljΔ+tl1+1kkkkkkkkkkk-α1+⋯+tlj+1-αj(α1+⋯+tlj+1-αj))]2dt⩽C∫tktk+1‍(∑s=1j-1‍(ls+1-ls)s2Δ2kkkkkkkkkkkkk+j2Δ2(k-lj)+9j2Δ2∑s=1j-1‍)dt=CΔ2+o(Δ2),
where the constant C depends only on T, A, αj, ρj  (j=1,2,…,N) and o(x)/x→0 as x→0. In the similar way, we can get the same estimates for the first and last terms in Ij. Since the total number of the terms in Ij is less than ((αj+1-αj)/Δ)+2, Ij→0 as Δ→0. Hence, we complete the proof of this lemma.Lemma 10.
For any functiong∈C([0,T]), ∫0T‍hΔ(t)g(t)dt converges to ∫0T‍h(t)g(t)dt a.s. as Δ→0.Proof.
By integration by parts, we know(72)∫0T‍hΔ(t)g(t)dt=G(T)hΔ(T)-∑k=1N‍G(tk)|u(tk,x)-u(tk-1,x)|2,
where G(t)=∫0t‍g(s)ds. Therefore, for any Δ>0, the square integrable functional ∫0T‍hΔ(t)g(t)dt of N dimensional Brownian motion B=(B1,B2,…,BN) belongs to the direct sum of R and the second chaos.
To this end, we define the second quantization operatorΓ(α) (for 0≤α≤1) from L2 to L2 by
(73)Γ(α)F=∑n=0∞‍αnJnF,
where F∈L2(Ω).
Settingm=α-2+1≥2, then, by the hypercontractivity of OU semigroup (see [9, 10]), one has
(74)∥Γ(α)F∥m≤∥F∥2.
Letting F=∫0T‍hΔ(t)g(t)dt-∫0T‍h(t)g(t)dt∈Range(J0⊕J2), we have
(75)Γ(α)F=𝔼F+α2(F-𝔼F)=(1-α2)𝔼F+α2F.
Thus, the hypercontractivity justifies the following form:
(76)∥α2F∥m⩽∥F∥2+∥(1-α2)𝔼F∥m⩽∥F∥2+∥(1-α2)𝔼F∥2⩽CΔ1/2.
In the last inequality, we use the facts appearing in the end of the proof of Lemma 9. Therefore,
(77)𝔼|F|m≤C(m-1)mΔm/2.
Now, we identify the Δ→0 with 1/n→0 as n→0. Denote the set {|F|>n-1/8} by An,n=1,2,…. Then, P(An)≤C/nm/2. Choose m large enough such that ∑0∞‍P(An)<∞ holds. Thus, the Borel-Cantelli lemma can be applied to show F=∫0T‍hΔ(t)g(t)dt-∫0T‍h(t)g(t)dt→0 a.s. as Δ tends to zero.Proof of Theorem8.
As the proof of Lemma10, we just need to show that the limit holds in the sense of L2. Actually, as Δ→0, one has
(78)|u(T,x)-u(tMΔ-1,x)|2∑i=1MΔ‍g′(ti)=(|u(T,x)-u(tMΔ-1,x)|2-NΔ)×∑i=1MΔ‍g′(ti)+NΔ∑i=1MΔ‍g′(ti)⟶N∫0T‍g′(t)dt=Ng(T)=h′(T)g(T).
Hence, the result of Theorem 8 is concluded.SinceWeyl’s perturbation theorem (see [11, 12]) says that
(79)max⁡k|λk↑(Δ)-λk↑|⩽∥HΔ,0-1HΔ,1-H0-1H1∥,
by Theorem 8, we can get the convergence of estimations of the locations by an obvious way, which is stated as follows.Theorem 11.
For anyk=1,2,…,N, one has SΔ,k↑→Sk a.s. as Δ→0.Next, we will give the convergence of estimators of the intensities of the point sources.Theorem 12.
Letρ(Δ)=(ρ1(Δ),ρ2(Δ),…,ρN(Δ)) be the solution of linear equations μ0(Δ)=A0(Δ)ρ(Δ); that is, ρ(Δ)=A0-1(Δ)μ0(Δ). Then, one has ρ(Δ)→ρ a.s. as Δ→0.Proof.
Notice thatρ solves the equation μ0=A0ρ. By Theorems 8 and 11, one can get the conclusion of this theorem.



---

*Source: 219876-2014-04-22.xml*