A Sampling Framework for Solving Physics-driven Inverse Source Problems

TL;DR

This paper introduces a novel sampling framework leveraging frequency estimation and Green's functions to solve inverse source problems in physics-driven fields, enabling accurate source recovery from sensor data.

Contribution

The paper develops a new sampling framework that uses multidimensional frequency estimation and Green's functions to solve inverse source problems for PDE-governed fields.

Findings

Framework effectively recovers sources from sensor data.

Explicit formulas for sensor-to-measurement mapping are derived.

Numerical simulations confirm robustness and accuracy.

Abstract

Partial differential equations are central to describing many physical phenomena. In many applications these phenomena are observed through a sensor network, with the aim of inferring their underlying properties. Leveraging from certain results in sampling and approximation theory, we present a new framework for solving a class of inverse source problems for physical fields governed by linear partial differential equations. Specifically, we demonstrate that the unknown field sources can be recovered from a sequence of, so called, generalised measurements by using multidimensional frequency estimation techniques. Next we show that---for physics-driven fields---this sequence of generalised measurements can be estimated by computing a linear weighted-sum of the sensor measurements; whereby the exact weights (of the sums) correspond to those that reproduce multidimensional exponentials, when used to linearly combine translates of a particular prototype function related to the Green's function of our underlying field. Explicit formulae are then derived for the sequence of weights, that map sensor samples to the exact sequence of generalised measurements when the Green's function satisfies the generalised Strang-Fix condition. Otherwise, the same mapping yields a close approximation of the generalised measurements. Based on this new framework we develop practical, noise robust, sensor network strategies for solving the inverse source problem, and then present numerical simulation results to verify their performance.

Equations

$$\mathcal{D}u(\mathbf{x},t)=f(\mathbf{x},t),$$

$$u(\mathbf{x},t)=(g\ast f)(\mathbf{x},t),$$

$${\nabla }^{2}u(\mathbf{x})=f(\mathbf{x}).$$

$$g(\mathbf{x})=\frac{1}{2\pi }\log \left(\Vert \mathbf{x}\Vert \right),$$

$$g(\mathbf{x})=-\frac{1}{4\pi \Vert \mathbf{x}\Vert }.$$

$$\frac{\partial }{\partial t}u(\mathbf{x},t)=\mu {\nabla }^{2}u(\mathbf{x},t)+f(\mathbf{x},t),$$

$$g(\mathbf{x},t)=\frac{1}{(4\pi \mu t{)}^{d/2}}{e}^{-\frac{{\Vert \mathbf{x}\Vert }^2}{4\mu t}}H(t),$$

$${\nabla }^{2}u(\mathbf{x},t)-\frac{1}{c^2}\frac{{\partial }^2}{\partial t^2}u(\mathbf{x},t)=f(\mathbf{x},t),$$

$$g(\mathbf{x},t)=\frac{c}{2\pi \sqrt{c^2{t}^2-\Vert \mathbf{x}{\Vert }^2}}H(ct-\Vert \mathbf{x}\Vert ).$$

$$g(\mathbf{x},t)=\frac{1}{4\pi \Vert \mathbf{x}\Vert }\delta (t-\Vert \mathbf{x}\Vert /c).$$

$${\phi }_n(t_l)=u({\mathbf{x}}_n,t_l).$$

$${\phi }_{n,l}^{ϵ}={\phi }_n(t_l)+{ϵ}_{n,l},$$

$$\mathrm{SNR}\stackrel{\text{def}}{=}10{\log }_{10}\left(\frac{{\sum }_{n=1}^N{\sum }_{l=0}^L|{\phi }_n(t_l){|}^2}{N(L+1){\sigma }^2}\right).$$

$$f(\mathbf{x},t)=\sum _{m=1}^M{c}_m\delta (\mathbf{x}-{\mathbf{\xi }}_m,t-{\tau }_m),$$

$$\mathcal{Q}(\mathbf{k},r)={⟨f(\mathbf{x},t),{\Psi }_{\mathbf{k}}(\mathbf{x}){\Gamma }_r(t)⟩}_{\mathbf{x},t},$$

$$\mathcal{Q}(\mathbf{k},r)=\sum _{m=1}^M{c}_m{\Psi }_{\mathbf{k}}({\mathbf{\xi }}_m){\Gamma }_r({\tau }_m).$$

$$\mathcal{Q}(\mathbf{k},r)=\sum _{m=1}^M{c}_m{e}^{\mathrm{j}r{\tau }_m/T}{e}^{\mathrm{j}{k}_1{\xi }_{1,m}+\mathrm{j}{k}_2{\xi }_{2,m}+\cdots +\mathrm{j}{k}_d{\xi }_{d,m}}.$$

$$\mathcal{Q}(\mathbf{k},1)\stackrel{\text{def}}{=}\mathcal{Q}(k_1,k_2,\dots ,k_d,1)=\sum _{m=1}^M{b}_m\prod _{i=1}^d{v}_{i,m}^{k_i},$$

$$\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r){\phi }_n(t_l)=\hat{\mathcal{Q}}(\mathbf{k},r),$$

$$\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r)g({\mathbf{x}}_n-\mathbf{x},t_l-t)\equiv {\Psi }_{\mathbf{k}}(\mathbf{x}){\Gamma }_r(t).$$

$$u(\mathbf{x},t)$$

$$={\int }_{{\mathbf{x}}^{\prime }\in {\mathbb{R}}^d}{\int }_{t^{\prime }\in \mathbb{R}}g({\mathbf{x}}^{\prime },t^{\prime })f(\mathbf{x}-{\mathbf{x}}^{\prime },t-t^{\prime })\mathrm{d}{t}^{\prime }\mathrm{d}{\mathbf{x}}^{\prime }$$

$$=⟨f({\mathbf{x}}^{\prime },t^{\prime }),g(\mathbf{x}-{\mathbf{x}}^{\prime },t-t^{\prime }){⟩}_{{\mathbf{x}}^{\prime },t^{\prime }},$$

$${\phi }_n(t_l)=u({\mathbf{x}}_n,t_l)=⟨f(\mathbf{x},t),g({\mathbf{x}}_n-\mathbf{x},t_l-t){⟩}_{\mathbf{x},t}.$$

$$\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r){\phi }_n(t_l)=\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r)u({\mathbf{x}}_n,t_l)$$

$$=\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r){⟨f(\mathbf{x},t),g({\mathbf{x}}_n-\mathbf{x},t_l-t)⟩}_{\mathbf{x},t}$$

$$={⟨f(\mathbf{x},t),\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r)g({\mathbf{x}}_n-\mathbf{x},t_l-t)⟩}_{\mathbf{x},t},$$

$$\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r)g({\mathbf{x}}_n-\mathbf{x},t_l-t)\equiv {\Psi }_{\mathbf{k}}(\mathbf{x}){\Gamma }_r(t)$$

$$\sum _{n=1}^N\sum _{l=0}^L{w}_{n,l}(\mathbf{k},r)g({\mathbf{x}}_n-\mathbf{x},t_l-t)=e^{\mathrm{j}\mathbf{k}\cdot \mathbf{x}}{e}^{\mathrm{j}rt/T}$$

$$\sum _{n,l}{w}_{n,l}(\mathbf{k},r)g({\mathbf{x}}_n-\mathbf{x},t_l-t)=e^{\mathrm{j}(\mathbf{k},r/T)\cdot (\mathbf{x},t)}=e^{\mathbf{\kappa }\cdot (\mathbf{x},t)},$$

$$V_{{\mathbf{\Delta }}_{\mathbf{x}}}(p)={\mathrm{span}}_{\mathbf{n}\in {\mathbb{Z}}^d}\left\{p\left(\mathbf{x}/{\mathbf{\Delta }}_{\mathbf{x}}-\mathbf{n}\right)\right\}\text{,}$$

$$\tilde{h}(\mathbf{x})=\sum _{\mathbf{n}\in {\mathbb{Z}}^d}{a}_{\mathbf{n}}p\left(\mathbf{x}/{\mathbf{\Delta }}_{\mathbf{x}}-\mathbf{n}\right),$$

Full text

I Introduction

Sensor networks, and the use thereof, for sensing and monitoring physical fields is receiving significant research attention due, in part, to the significant advances made over the last few decades in the fields of (wireless) networking, communications and in the fabrication of microprocessors [3, 4]. During this period many interesting applications in localisation, tracking and parameter estimation have been considered [5, 6]. The sensor nodes are deployed over a region of interest to obtain spatiotemporal samples of some physical phenomena. Often, these phenomena are driven by natural mechanisms that typically involve the transportation of matter/particles or the transportation of energy, from one point to another and thus can be described by partial differential equations (PDEs). For example, the mode of transport governing the dispersion of plumes in environmental monitoring [7], spreading of fungal diseases in precision agriculture [8], biochemical and nuclear wastes [9] is well-known to be diffusion and the corresponding diffusion field is the variation in concentration of the released substance over space and time.

Besides diffusion, there exist numerous modes of transport and corresponding fields, such as wave and potential fields that have also received considerable research attention from the signal processing community. Such efforts have focussed on developing robust sensor data fusion schemes that either infer the sources inducing the measured field or directly reconstruct the field. Often these physical fields of interest are spatially non-bandlimited and so require an extremely dense set of samples in order to achieve a faithful recovery, using the classical linear bandlimited (BL) reconstruction framework.

In the localisation of neuronal source activities from electroencephalographic (EEG) signals [10, 11]—for bioengineering applications—the use of Poisson’s equation to model the brain activity is important, since this PDE accurately describes the relationship between the measured electrical potentials (the field) and the current dipoles (the sources). To alleviate some of the limitations of a BL reconstruction when solving this EEG-related inverse source problem (ISP), many approaches utilising least-squares [12], sparsity-based [13, 14] and, more recently, cosparsity-based [15, 16] regularisation have been proposed. Moreover, Bayesian modelling [17], beamforming (see [18, 19] and references therein), as well as subspace techniques [20] have also been explored.

Wave fields are also prevalent in applications such as acoustic tomography [21], speech and sound enhancement [22], sound/wave source localisation [22] and many more. In such situations the wave equation provides a physical law that describes the propagation of such fields. To solve the associated ISPs, several interesting techniques have been put forward; the classical and most commonly used techniques are based on maximum likelihood estimation and beamforming [23]. Recently however, Dokmanic̀ et al [24, 25] proposed an approach that exploits certain salient properties of euclidean distance matrices to solve the simultaneous localisation and mapping problem for acoustic sources. Kitić et al [16] also formulated an optimisation problem by using the cosparse regularisation framework, whilst a finite rate of innovation (FRI) based method is introduced in [26] to solve this acoustic ISP using boundary-only measurements.

For the diffusion field reconstruction problem, Reise et al [27] propose the use of hybrid shift-invariant subspaces, whilst recovery algorithms based on the use of finite element methods (FEMs) [28, 29] and compressive sensing (CS) [30, 31] have also been proposed. Moreover statistical estimation techniques, based on Bayesian estimation and Kalman filtering, have also been studied, see [32, 33, 34] for instance. Meanwhile, Dokmanic et al [35] retrieve the single source parameters by approximating the resulting field using a truncated Fourier series, and Lu et al demonstrate that by solving a set of linear equations the single diffusion source parameters can be estimated [36]. Furthermore, in [37] we showed that given a proper sequence of generalised measurements it is possible to recover the unknown parameters for a specific class of diffusion source distributions. We will leverage from that concept, in this current work, to devise a framework for solving a more general class of physics-driven ISPs. Specifically, the new approach will apply to higher dimensional ISPs governed by a larger class of PDE models.

The rest of this paper is organised as follows. In Section II, we give mathematical descriptions of some common linear, constant coefficient PDEs encountered in many applications and then state the related ISP and also discuss the sensor network model. In Section III we outline how to solve the ISP of interest assuming we have access to a set of generalised measurements. Specifically, we discuss and state explicitly how to select properly the spatiotemporal sensing functions in order to be able to solve the $d$-dimensional ISP. Then we explore a new approach for computing the desired generalised measurements in Section IV, based on taking proper linear combinations of the sensor data. In particular we realise that this leads to the well-known exponential reproduction problem—using translates of a prototype function—encountered in approximation theory and in the FRI framework [38]. Here however, the prototype function coincides with the space- and time-reversed Green’s function of the physical field. Then we derive conditions on the Green’s function, for which this exponential reproduction problem can be exactly or approximately solved. Section V, discusses how to adapt this new framework to solving ISPs using sensor networks. We develop explicit centralised and distributed estimation strategies whilst considering both uniform and nonuniform sensor placements. Then in Section VI, we provide numerical simulation results to validate the new framework, comparing it against a sparsity-based recovery method; and finally, we conclude the paper in Section VII.

II Physics-driven inverse problems: Problem formulation

The term physics-driven is used in this paper to describe physical phenomena, specifically physical fields, that propagate through space and time according to some linear partial differential equation (PDE). In a more general form, such phenomena can be written as,

[TABLE]

where $\mathcal{D}$ denotes a linear differential operator, whilst $f(\mathbf{x},t)$ is the source of the field $u(\mathbf{x},t)$ propagating through space $\mathbf{x}\in\Omega\subset\mathbb{R}^{d}$ and time $t\in\mathbb{R}_{+}$. According to the method of Green’s functions, under certain boundary conditions, the system (1) admits the solution

[TABLE]

where $g(\mathbf{x},t)$ is the so called Green’s function of the underlying field. The following fields will be considered in this paper:

1. (a)

Potential fields: these are encountered frequently in many situations arising in electrostatics. Mathematically, the potential field satisfies

[TABLE]

The Green’s function for this PDE in 2D (i.e. $d=2$) is,

[TABLE]

whilst for $d=3$ the Green’s function becomes

[TABLE] 2. (b)

Diffusion fields: refer to physical phenomena such as the propagation of heat, plumes and leakages that can be described mathematically by,

[TABLE]

with the following Green’s function, for $d\geq 1$:

[TABLE]

where $H(t)$ is the unit step and $\mu$ is the diffusivity. 3. (c)

Wave fields: describe many situations arising for example in acoustics and electromagnetism. The wave equation is given by:

[TABLE]

where the wave field $u(\mathbf{x},t)$ induced by the source distribution $f(\mathbf{x},t)$ propagates through the medium at a speed $c$. The Green’s function for the $2$-D wave equation (i.e. $d=2$) is given by:

[TABLE]

Moreover, in $3$-D (i.e. $d=3$), it can be shown that:

[TABLE]

The Green’s functions above assume a Sommerfeld radiation condition – i.e. a quiescent condition at an initial time, such that $\left.u(\mathbf{x},t)\right|_{t=0}=\left.\frac{\partial}{\partial t}u(\mathbf{x},t)\right|_{t=0}=0$ and a convergence condition at infinity, meaning $\left.u(\mathbf{x},t)\right|_{\|\mathbf{x}\|\rightarrow\infty}=\left.\frac{\partial}{\partial x_{1}}u(\mathbf{x},t)\right|_{\|\mathbf{x}\|\rightarrow\infty}=\left.\frac{\partial}{\partial x_{2}}u(\mathbf{x},t)\right|_{\|\mathbf{x}\|\rightarrow\infty}=0$. See for example, [39] for a derivation of these expressions.

Given the PDE models of such fields, the aim of this paper is to develop a framework for solving the associated ISP, from spatiotemporal sensor network measurements of the induced field111For instance using a suitable microphone array for audio fields, or an array of thermal sensors to monitor the temperature of a room.. Precisely, the ISP considered here is the following:

Problem 1

Let $\mathcal{S}=\{\mathbf{x}_{n}\}_{n=1}^{N}$ denote a network of $N$ sensors, so that the $n$-th sensor situated at $\mathbf{x}_{n}$ collects samples $\varphi_{n}(t_{l})=u(\mathbf{x}_{n},t_{l})$ of the field $u$, at times $t_{l}$ for $l=0,1,\ldots,L$. Given these spatiotemporal samples and knowledge of the Green’s function of the field, we intend to estimate the unknown source distribution $f(\mathbf{x},t)$.

II-A Sensor Network Model and Assumptions

The sensor networks used to monitor our physics-driven fields are such that:

1. (a)

They comprise $N$ sensor nodes deployed (uniformly or randomly) over the region of interest ($\Omega\subset\mathbb{R}^{d}$). For example in $2$-D the sensors all lie in the same plane. 2. (b)

The sensor locations $\mathbf{x}_{n}\in\Omega$ are known and each sensor samples the field locally at time instants $t_{l}$ for $l=0,1,\ldots,L$. Hence the noiseless field samples are simply the field (2), evaluated at $\mathbf{x}=\mathbf{x}_{n}$ and $t=t_{l}$ as follows:

[TABLE] 3. (c)

The sensor noise can be modelled by a zero mean additive white Gaussian noise (AWGN) process, so that the noisy measurements are:

[TABLE]

where $\epsilon_{n,l}\sim\mathcal{N}(0,\sigma^{2})$ and the (average) signal-to-noise ratio ($\mathrm{SNR}$) is

[TABLE] 4. (d)

The sensor nodes are synchronised. Hence the sensors sample the field at the same instants.

Within this setting, we consider two scenarios: a centralised scenario where we assume that all sensors’ readings are available at a fusion centre at which all the processing is performed and a distributed scenario where sensors can perform processing locally but can only communicate with neighbouring sensors. In this distributed setup, the sensors must recover the unknown distribution $f(\mathbf{x},t)$ through localised/in-network data processing and communications alone.

To model the network topology, in the distributed setting, we assume a connected random geometric graph (RGG), denoted by $\mathcal{G}(N,r_{\mathrm{con}})$, with $N$ sensor nodes and connectivity radius $r_{\mathrm{con}}$. To realise this, we place $N$ nodes uniformly at random over a unit square/cube and then put an edge between a pair of nodes if their Euclidean distance is at most $r_{\mathrm{con}}$. An example, with $N=10$ nodes is shown in Figure 1, the shaded circular region is the communication radius of the red sensor.

III Choosing the sensing functions: Source recovery from generalised measurements

The focus of this paper will be predominantly on fields induced by multiple localised and instantaneous sources. This source distribution $f(\mathbf{x},t)$ admits the parametrisation:

[TABLE]

where $M$ is the total number of sources, $c_{m},\tau_{m}\in\mathbb{R}$ and $\bm{\xi}_{m}=(\xi_{i,m})_{i=1}^{d}\in\mathbb{R}^{d}$ are the intensity, activation time and location of the $m$-th source respectively. As a result of this parametrisation, the ISP becomes one of estimating $M$ triples $\{(c_{m},\tau_{m},\bm{\xi}_{m})\}_{m=1}^{M}$.

The proposed framework will be based on estimating the unknown source parameters from the following multidimensional sequence of generalised measurements,

[TABLE]

where $\{\Psi_{\bm{k}}(\mathbf{x})\}_{\bm{k}\in\mathbb{N}^{d}}$ and $\{\Gamma_{r}(t)\}_{r\in\mathbb{N}}$, are families of properly chosen spatial and temporal sensing functions respectively. We will discuss how to choose these in the sequel.

Moreover, observe that when we substitute (14) into the inner product (15), we obtain:

[TABLE]

Hence, our first task will be to select $\Psi_{\bm{k}}(\mathbf{x})$ and $\Gamma_{r}(t)$ so that we can recover $\{c_{m},\tau_{m},\bm{\xi}_{m}\}_{m=1}^{M}$ from $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k},r}$. We then discuss in Section V how to obtain the measurements (15) from the sensors’ readings. In the same spirit of [37], we propose the use of exponentials with purely imaginary exponents as our sensing functions for two main reasons. The first is that, the sum (16) becomes a superposition of multidimensional exponentials, i.e. a multidimensional system of superimposed sinusoids, which can be efficiently solved by using multidimensional extensions of Prony’s frequency estimation methods [40, 41, 42]. The second reason is because the use of imaginary exponentials improves the stability of the estimation problem, since in this case the magnitude of the terms in the generalised sequence remain bounded. We now examine explicitly the $d$-dimensional source recovery problem. Note however that, in most natural applications, one would mostly be interested in two- and three-dimensional fields.

III-A Sensing sources in space and time

The temporal sensing function is chosen to be the exponential $\Gamma_{r}(t)=e^{\mathrm{j}rt/T}$, where $T=t_{L}$ is the instant at which the sensors measure the last sample of the field. In addition, we choose $\Psi_{\bm{k}}(\mathbf{x})=e^{\mathrm{j}\bm{k}\cdot\mathbf{x}}$ with $\bm{k}\stackrel{{\scriptstyle\text{\tiny{def}}}}{{=}}(k_{1},k_{2},\ldots,k_{d})\in\mathbb{N}^{d}$. Under this particular selection of spatial and temporal sensing functions, the expression (16) becomes:

[TABLE]

Notice now that, for some fixed $r\neq 0$ (take, for instance, the particular case $r=1$) then expression (17) is of the form:

[TABLE]

with $b_{m}{=}c_{m}e^{\mathrm{j}\tau_{m}/T}$ and $v_{i,m}{=}e^{\mathrm{j}\xi_{i,m}}$. This is a multidimensional Prony-like system which, as described in Appendix A, can be solved to obtain $\{(b_{m},v_{1,m},v_{2,m},\ldots,v_{d,m})\}_{m=1}^{M}$ simultaneously from $\{\mathcal{Q}(\bm{k},1)\}_{\bm{k}}$ provided $k_{i}{=}0,1,\ldots,K_{i}$ and $K_{i}{\geq}2M{-}1$ for any $i\in{1,2,\ldots,d}$. See also [41, 40], for more on the topic of multidimensional frequency estimation. Then it follows immediately that the unknown source parameters are given by: $c_{m}=|b_{m}|$, $\tau_{m}=T\arg(b_{m})$ and $\bm{\xi}_{m}=-\mathrm{j}\left(\log(v_{1,m}),\log(v_{2,m}),\ldots,\log(v_{d,m})\right)$.

Having outlined how to recover the unknown point source parameters from $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k}}$, we must now focus on retrieving these generalised measurements from the spatiotemporal samples of the field.

In [37, 43] it was shown that for the case of the diffusion field, the generalised measurements can be found by imposing that $\Psi_{k}(\mathbf{x})$ be analytic. A similar strategy has been used in [26, 44] for the wave and Poisson equation. The disadvantage of these approaches is that

(i) they cannot be easily extended to $d>2$ and (ii) the constraint that $\Psi_{k}$ is analytic leads in some cases to less stable reconstruction algorithms. Consequently, in the following section we outline a new and more versatile approach to compute $\mathcal{Q}(\bm{k},r)$ from the sensor data.

Notice that if we allow $k_{1}$ to be imaginary and impose, $k_{1}=-\mathrm{j}k_{2}=\mathrm{j}k$ we obtain the same $\mathcal{Q}(k,r)$ of [37], so this new approach is a generalisation.

IV Multidimensional generalised measurements from sensor samples

To compute the desired set of multidimensional generalised measurements we consider taking weighted linear combinations of the spatiotemporal samples:

[TABLE]

hence the goal is to find the weights $\{w_{n,l}(\bm{k},r)\}_{n,l}$, such that the left hand side of (18), that is $\widehat{\mathcal{Q}}(\bm{k},r)$, coincides with the right hand side of (15) for any $\bm{k}$ and $r$. The result below follows from this consideration.

Proposition 1

Computing the multidimensional sequence of generalised measurements $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k}}$, in (15) for any $r\in\mathbb{N}$, by taking weighted linear combinations (18) of the sensor data $\varphi_{n}(t_{l})$ is equivalent to reproducing the function $\Psi_{\bm{k}}(\mathbf{x})\Gamma_{r}(t)$ from space- and time-reversed translates of the Green’s function $g(\mathbf{x},t)$ of the underlying field. Specifically, imposing $\widehat{\mathcal{Q}}(\bm{k},r)={\mathcal{Q}}(\bm{k},r)$ implies that:

[TABLE]

Furthermore when $\Psi_{\bm{k}}(\mathbf{x})$ and $\Gamma_{r}(t)$ are chosen to be exponentials, this results in a multidimensional exponential reproduction problem.

Proof:

We commence this proof by noting that (2) can be written as:

[TABLE]

where $\mathbf{x}^{\prime}=(x^{\prime}_{1},x^{\prime}_{2},\ldots,x^{\prime}_{d})\in\mathbb{R}^{d}$ and $\mathrm{d}\mathbf{x}^{\prime}=\prod_{i=1}^{d}\mathrm{d}x_{i}^{\prime}$. Consequently, the discrete measurement obtained by the $n$-th sensor (located at $\mathbf{x}_{n}$) at some time instant $t_{l}\geq 0$ is

[TABLE]

Replacing (19) into the left hand side (lhs) of (18) yields:

[TABLE]

where $\{w_{n,l}(\bm{k},r)\}_{n,l}\in\mathbb{C}$ denote the specific sequence of weights we wish to compute222In the last equality, we are able to pass the summation inside the inner product because it is finite. If it were infinite then we would require that the sum converges absolutely; which is ensured if $g$ and its translates form a Riesz basis.. In particular if we require this weighted sum of the sensor data to yield the exact multidimensional measurements, i.e. $\widehat{\mathcal{Q}}(\bm{k},r)=\mathcal{Q}(\bm{k},r)=\left\langle f(\mathbf{x},t),\Psi_{\bm{k}}(\mathbf{x})\Gamma_{r}(t)\right\rangle_{\mathbf{x},t}$, then by comparing the inner products in (15) and (20), we realise that we must choose the sequence of weights $\{w_{n,l}(\bm{k},r)\}_{n,l}$ such that, for each $\bm{k}$ and $r$, the identity

[TABLE]

is satisfied. This proves the first claim of the proposition.

The second statement follows immediately. If we impose the choice of sensing functions of Section III, wherein $\Psi_{\bm{k}}(\mathbf{x})=e^{\mathrm{j}\bm{k}\cdot\mathbf{x}}$ and $\Gamma_{r}(t)=e^{\mathrm{j}rt/T}$, then (21) becomes

[TABLE]

which can be written in a more compact form as follows:

[TABLE]

*where $\bm{\kappa}=\mathrm{j}(\bm{k},r/T)\in\mathbb{R}^{d+1}$. Hence we can immediately observe that the required coefficients are those that reproduce the $(d+1)$-dimensional (space and time varying) exponentials by summing translates of $g(-\mathbf{x},-t)$, which is the space- and time-reversed Green’s function of the underlying field $u(\mathbf{x},t)$.

∎*

Hence with access to the desired coefficients, $\{w_{n,l}(\bm{k},r)\}_{n,l}$, that are capable of reproducing exponentials from the translates of $g(\mathbf{x},t)$, all we would need to do is: evaluate the sequence $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k}}$ for a fixed $r\neq 0$ using (18) and then from $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k}}$, extract the unknown source parameters as described in Section III. Therefore the only missing piece in our framework is how to obtain the exponential reproducing coefficients. For this, one needs to understand when the exponential reproduction problem is exactly or approximately feasible and then establish appropriate schemes to find the desired coefficients.

We address these questions by leveraging from results in generalised sampling and approximation theory. Then for both (exact and approximate reproduction) cases, we derive closed-form formulae to compute the “best” weights $w_{n,l}(\bm{k},r)$, when the translates of the approximant are assumed to be regular. In the sensor network setup this is equivalent to having uniform spatiotemporal samples with sampling intervals $\bm{\Delta}_{\mathbf{x}}=\left(\Delta_{x_{1}},\Delta_{x_{2}},\ldots,\Delta_{x_{d}}\right)$ and $\Delta_{t}$. Finally in the nonuniform sampling case, where it is generally not possible to obtain simple closed-form expressions for the desired exponential reproducing coefficients $w_{n,l}(\bm{k},r)$, we propose two approaches:

1. (a)

Formulating and solving the linear system that comes from discretizing (23). 2. (b)

Interpolating and resampling the sensor data uniformly.

IV-A Function spaces, generalised sampling and function approximation

In the generalised (uniform) sampling paradigm—see [45, 46, 47] and references therein—the primary goal is to reconstruct some functions of a continuous variable from a discrete set of measurements collected on a uniform grid. Often, this reconstruction will be an approximation of the original signal in some function spaces with a further property that the approximation error decays to zero, in the limit as the “density” of the sampling grid increases.

Consider the $d$-dimensional generating function $p$, whose uniform translates generates the space

[TABLE]

then any function $\hat{h}(\mathbf{x})\in V_{\bm{\Delta}_{\mathbf{x}}}(p)\subset\mathcal{L}^{2}$ is characterised by the sequence of coefficients $a_{\bm{n}}$, such that:

[TABLE]

where $\mathcal{L}^{2}$ denotes the space of square-integrable functions. In fact for this series to be well posed we require that:

1. (a)

For convergence, $\left\{a_{\bm{n}}\right\}_{\bm{n}}$ must be square-summable. 2. (b)

For uniqueness and stability of this discrete representation, $\left\{p(\mathbf{x}-\bm{n})\right\}_{\bm{n}}$ must form a Riesz basis of $V_{\bm{1}}(p)$. 3. (c)

Finally, $p(\mathbf{x})$ must satisfy the partition of unity condition

[TABLE]

for all $\mathbf{x}\in\mathbb{R}^{d}$ in order to guarantee that by choosing $\bm{\Delta}_{\mathbf{x}}$ in (24) sufficiently small, we can approximate any function $\tilde{h}(\mathbf{x})$ as closely as we want (a detailed proof of this fact can be found in [46, Appendix B]).

If we instead want to reconstruct some signal $h(\mathbf{x})\in\mathcal{L}^{2}\setminus\!V_{\bm{\Delta}_{\mathbf{x}}}$, then the reconstruction $\tilde{h}(\mathbf{x})$ in (24) should produce the best approximation of $h(\mathbf{x})$ in the space $V_{\bm{\Delta}_{\mathbf{x}}}$, and hence minimise the approximation error in the least-squares sense. This is achieved by computing the orthogonal projections of $h(\mathbf{x})$ onto $V_{\bm{\Delta}_{\mathbf{x}}}(p)$ which is obtained by choosing $a_{\bm{n}}=\left\langle h(\mathbf{x}),p_{\text{dual}}\left({\mathbf{x}}/{\bm{\Delta}_{\mathbf{x}}}-\bm{n}\right)\right\rangle$, where $p_{\text{dual}}$ is the dual of $p$ and is given by [48],

[TABLE]

Here $\hat{p}(\bm{\omega})=\mathcal{F}_{\mathbf{x}}\{p\}=\int_{\mathbf{x}}p(\mathbf{x})e^{-\mathrm{j}\bm{\omega}\cdot\mathbf{x}}\mathrm{d}\mathbf{x}$ denotes the multidimensional Fourier transform of $p$. Please note also that in the rest of this paper, it is assumed that all transforms are taken in the sense of distributions.

Within our proposed framework we are seeking the specific coefficients $\{w_{n,l}(\bm{k},r)\}_{n,l}$ that reproduce the exponential function using shifted versions of the Green’s function of the underlying physical field. Hence this is a special case of the above, where:

1. (a)

the signal we want to reconstruct is a specific $(d+1)$-dimensional exponential (i.e. $e^{\bm{\kappa}\cdot(\mathbf{x},t)}$), and 2. (b)

the generating function $p$ is precisely the Green’s function of the underlying PDE.

Under these conditions we want to find the best representation of the exponentials, in the space spanned by the translates of the Green’s functions.

IV-B Exact and approximate Strang-Fix theory for exponential reproduction from uniform translates

Building on the discussion of the previous section, we now focus on the approximation of exponentials from uniform translates of a single prototype function. We begin by proving the following lemma which is an extension of [49] using a multidimensional generalisation of the proof given in [38].

Lemma 1** (Generalised multidimensional Strang-Fix conditions [50])**

Let $p(\mathbf{x})$ be of compact support and its multidimensional bilateral Laplace transform be $P(\mathbf{s})=\int_{\mathbf{x}\in\mathbb{R}^{d}}g(\mathbf{x},t)e^{-\mathbf{x}\cdot\mathbf{s}}\mathrm{d}\mathbf{x}$, then the following conditions are equivalent:

1. (a)

For any $\bm{n}\in\mathbb{Z}^{d}\setminus\!\{\bm{0}\}$, where $\bm{0}$ is the zero $d$-vector,

[TABLE] 2. (b)

For some coefficients $a_{\bm{n}}\in\mathbb{C}$,

[TABLE]

is an exponential in $\mathbf{x}$, i.e. $\sum_{\bm{n}\in\mathbb{Z}^{d}}{a_{\bm{n}}p(\mathbf{x}-\bm{n})}=Ce^{\bm{\kappa}\cdot\mathbf{x}}$, for some $C\neq 0$.

Proof:

See Appendix C. ∎

The assumption of compactness here ensures first that the bilateral Laplace transform exists, so that the conditions (26) are well-defined and second, that (27) converges. This assumption is sufficient, but not necessary, since for example a suitable polynomial decay [51], or even milder restriction on $p(\mathbf{x})$ [52, 53] would still guarantee convergence. For the class of functions satisfying the generalised Strang-Fix conditions above, the desired coefficients can be computed exactly. To treat our multidimensional problem, we now extend formally the one-dimensional formulae obtained in [38] to the multidimensional case. This new formulae are still valid even when the prototype function is not separable with respect to its variables – i.e we do not require that $p(\mathbf{x})=\prod_{i=1}^{d}{p_{i}(x_{i})}$. The absence of this separability property is of paramount importance for us, especially because the spatial and temporal dimensions for most non-static fields encountered, in reality, are neither separable nor homogeneous. We begin our derivation by recalling that we are after the coefficients $\{a_{\bm{n}}\}_{\bm{n}}$ such that

[TABLE]

According to generalised sampling theory, the sequence of weights that minimises the approximation error in the least-squares sense is given by

[TABLE]

where the second line follows from the change of variable $\mathbf{x}^{\prime}=\mathbf{x}-\bm{n}$. Thus finding $a_{\bm{0}}$ allows us to compute $a_{\bm{n}}$ for all $\bm{n}\in\mathbb{Z}^{d}$ using (29). To find $a_{\bm{0}}$ substitute (29) into (28) to get

[TABLE]

We then apply Poisson summation formula on the lattice to the l.h.s. of this expression, which if $p(\mathbf{x})$ is well-behaved reduces to

[TABLE]

Finally from the condition (26) we get $a_{\bm{0}}=\frac{1}{P(\bm{\kappa})}$, since all terms in the summation vanish for $\bm{n}\neq\bm{0}$. Hence for any $\bm{n}\in\mathbb{Z}^{d}$ it follows that,

[TABLE]

IV-B1 Approximate Strang-Fix in multidimensions and the approximation error

In the derivation of (30) we imposed some regularity conditions on $p(\mathbf{x})$, specifically for the l.h.s of the Poisson summation formula to converge, the function must decay sufficiently quickly. The strongest constraint on $p(\mathbf{x})$ however is due to (26), where $P(\bm{\kappa}+\mathrm{j}2\pi\bm{n})=0$ for $\bm{n}\in\mathbb{Z}^{d}\setminus\!\{\bm{0}\}$.

For general physical fields of interest to us, the approximant $p(\mathbf{x})$ will be replaced by the corresponding Green’s function $g$ of the field. Whilst these will generally not satisfy the Strang-Fix condition (26), we still wish to reproduce exponentials approximately with them. To this end, we will extend the approximate Strang-Fix method introduced in [38], which relaxes the assumptions on the generators (for the $1$-D exponential case), so that we are now after the coefficients that gives the best approximate exponential reproduction, for any $p(\mathbf{x})$. Mathematically this means that we desire

[TABLE]

where $p(\mathbf{x})$ does not necessarily satisfy the generalised Strang-Fix conditions (26). There are a few possible choices one may make for the “best” approximation coefficients. For any choice, the associated approximation error is, $\varepsilon(\mathbf{x})=e^{\bm{\kappa}\cdot\mathbf{x}}\left(1-a_{\bm{0}}\sum_{\bm{n}}P(\bm{\kappa}+\mathrm{j}2\pi\bm{n})e^{\mathrm{j}2\pi\bm{n}\cdot\mathbf{x}}\right)$ which can be minimised in the least-squares sense, by computing the orthogonal projection of $e^{\bm{\kappa}\cdot\mathbf{x}}$ on the the subspace $V_{\bm{1}}(p)$. This yields $a_{\bm{0}}=\frac{P(-\bm{\kappa})}{R_{p}(e^{\bm{\kappa}})}$, where $R_{p}(e^{\bm{\kappa}})=\sum_{\bm{\ell}\in\mathbb{Z}^{d}}r_{p}[\bm{\ell}]e^{-\bm{\kappa}\cdot\bm{\ell}}$ is the multidimensional $z$-transform of the autocorrelation sequence $r_{p}[\bm{\ell}]=\left\langle p(\mathbf{x}-\bm{\ell}),p(\mathbf{x})\right\rangle_{\mathbf{x}}$, evaluated at $\mathbf{z}=e^{\bm{\kappa}}=(e^{\kappa_{1}},e^{\kappa_{2}},\ldots,e^{\kappa_{d}})$ [38, 53].

Moreover, observe that the square error $\varepsilon^{2}(\mathbf{x})$ is minimised when $1-a_{\bm{0}}\sum_{\bm{n}}P(\bm{\kappa}+\mathrm{j}2\pi\bm{n})e^{\mathrm{j}2\pi\bm{n}\cdot\mathbf{x}}=0$; in addition if the Laplace transform $P$ of the generator decays quickly, i.e. assuming $P(\bm{\kappa}+\mathrm{j}2\pi\bm{n})\approx 0$ for any $\bm{n}\in\mathbb{Z}^{d}\setminus\{\bm{0}\}$, then $a_{\bm{0}}=1/P(\bm{\kappa})$ is a good proxy for the minimiser of $\varepsilon(\mathbf{x})$ and therefore:

[TABLE]

These are the coefficients we shall utilise in this work, for their simplicity and accuracy. Moreover, the approximation error using (32) is given by:

[TABLE]

which will be small if $P(\bm{\kappa}+\mathrm{j}2\pi\bm{n})$ decays quickly to zero as $|\bm{n}|$ increases.

IV-C Computing the analysis coefficients for space-time fields

Equipped with the formulae in Section IV-B, we can now compute explicitly the weights $w_{n,l}(\bm{k},r)$ in (18). We distinguish two cases:

(i) the case where the sensors are uniformly spaced and (ii) the case where the location of the sensors are arbitrary.

IV-C1 Uniform sensor placement

Denote the sensor measurements by $\varphi_{\bm{n}}(t_{l})=u(\bm{n}\bm{\Delta}_{\mathbf{x}},l\Delta_{t})$, where $\bm{n}\bm{\Delta}_{\mathbf{x}}=(n_{1}\Delta_{x_{1}},n_{2}\Delta_{x_{2}},\ldots,n_{d}\Delta_{x_{d}})$ and $n_{i}=0,1,\ldots,N_{i}-1$ for $i=1,\ldots,d$. Note that we can reconcile the vector index sensor measurement $\varphi_{\bm{n}}(t_{l})$ with the scalar indexed one $\varphi_{n}(t_{l})$, by simply taking the lexicographic ordering of the elements of $\{\bm{n}\Delta_{\mathbf{x}}\}_{\bm{n}\in\mathbb{N}^{d}}$ to give $\{\mathbf{x}_{n}\}_{n=1}^{N}$, where $N=\prod_{i=1}^{d}N_{i}$.

Consequently for physics-driven fields, it is clear from (23) that the prototype function is the space- and time-reversed Green’s function: $p=g(-\mathbf{x},-t)$, with bilateral Laplace transform $G(-\mathbf{s}_{\mathbf{x}},-s_{t})$. Therefore for translates $\bm{\Delta}_{\mathbf{x}}\in\mathbb{R}_{+}^{d}$ and $\Delta_{t}\in\mathbb{R}_{+}$ the corresponding exponential reproducing coefficients are:

[TABLE]

where $\bar{g}(\mathbf{x},t)=g(\bm{\Delta}_{\mathbf{x}}\mathbf{x},\Delta_{t}t)$ and $\bar{\bm{\kappa}}=\mathrm{j}(\bm{\Delta}_{\mathbf{x}}\bm{k},\Delta_{t}r/T)$. Since, $\bar{g}(\mathbf{x},t)\Leftrightarrow\bar{G}(\mathbf{s}_{\mathbf{x}},s_{t})=\frac{1}{\Delta_{t}\prod_{i=1}^{d}\Delta_{x_{i}}}G(\frac{\mathbf{s}_{\mathbf{x}}}{\bm{\Delta}_{\mathbf{x}}},\frac{s_{t}}{\Delta_{t}})$, we can find $\bar{G}(-\mathrm{j}\bm{\Delta}_{\mathbf{x}}\bm{k},-\Delta_{t}\mathrm{j}r/T)$ and substitute it into (34), to conclude that for any $\bm{\Delta}_{\mathbf{x}}\in\mathbb{R}^{d}_{+}$ and $\Delta_{t}\in\mathbb{R}_{+}$,

[TABLE]

As an example, we show in Figure 2, the approximation of the 2-D spatial exponentials, $e^{\mathrm{j}k_{1}x_{1}+\mathrm{j}k_{2}x_{2}}$, using the 2-D Green’s function of Poisson’s equation (4) and the coefficients (35), for $r=0$, $k_{1}=2$, and $k_{2}=2,3,4$.

IV-C2 Least-squares scheme for nonuniformly placed sensors

For the case of non-uniformly placed sensors, it is generally not possible to find similar closed-form expressions for the desired coefficients $\{w_{n,l}(\bm{k},r)\}_{n,l}$. It is possible however to formulate a linear system of equations to find $\{w_{n,l}(\bm{k},r)\}_{n,l}$ since the approximating function $g$ and the exponentials $\Gamma_{r}(t)\Psi_{\bm{k}}(\mathbf{x})$ we want to approximate are known. One approach is to discretise (21) as follows. First, for each $l$, formulate the following linear system at some fixed time snapshot $t_{j}>0$:

[TABLE]

[TABLE]

Solving this system gives the coefficients $\{w_{n,l}(\bm{k},r)\}_{n}$ for a fixed $l$. Second, we stack (36) for each $l$ and take several $t_{j}$’s, $j=1,\ldots,J$ to get

[TABLE]

where $\mathbf{G}\in\mathbb{R}^{IJ\times N(L+1)}$ is a discretisation of $g(\mathbf{x},t)$, $\mathbf{p}(\bm{k},r)\in\mathbb{R}^{IJ}$ are discretisations of the spatiotemporal sensing functions, whilst $\mathbf{w}(\bm{k},r)\in\mathbb{R}^{N(L+1)}$ are the desired weights for each $\bm{k}$ and $r$. Consequently, in order to recover the desired field analysis coefficients, we would need to solve the system (50). In general, this system admits a least-squares solution if $IJ\geq N(L+1)$, where the observation matrix $\mathbf{G}$ can be constructed from the Green’s function of the problem at hand (i.e. (4), (5), (7), (9) and so on).

Although straightforward to formulate, the conditioning of such a system can be poor in some instances. Specifically, the condition number of $\mathbf{G}$ depends directly on the sensor locations $\mathbf{x}_{n}$, the sampling instants $t_{l}$ and the Green’s function $g(\mathbf{x},t)$ of the underlying phenomena.

IV-C3 Interpolation scheme for nonuniform sampling

Another simple, yet effective scheme for handling non-uniform sensor placement is to interpolate the field and resample it on a uniform grid. In so doing, we can obtain an approach that can still exploit the closed form expression (35) even when the spatial sampling is irregular.

Essentially we want to return to the situation where translates of the Green’s function are on a uniform lattice. To do this, we assume that the spatial field samples are interpolated, on a uniform grid, using the interpolator family $\{\gamma_{n}(\mathbf{x})\}$, such that $\hat{u}(\mathbf{x},t_{l})=\sum_{n=1}^{N}\varphi_{n}(t_{l})\gamma_{n}(\mathbf{x}-\mathbf{x}_{n})$. Then this new approximation of the underlying field is resampled uniformly at the new locations $\{\mathbf{x}_{\bar{\bm{n}}}=(\bar{n}_{1}\Delta_{x_{1}},\bar{n}_{2}\Delta_{x_{2}},\ldots,\bar{n}_{d}\Delta_{x_{d}})\}_{\bar{\bm{n}}\in\mathbb{N}^{d}}$ to give the corresponding data samples $\hat{\varphi}_{\bar{\bm{n}}}(t_{l})=\hat{u}(\mathbf{x}_{\bar{\bm{n}}},t_{l})$. Finally, since we are back to the uniform case, the corresponding exponential reproducing coefficients can be recovered using (35). Then weighting the interpolated measurements $\hat{\varphi}_{\bar{\bm{n}}}(t_{l})$ by the obtained coefficients will produce an estimate for the desired sequence of generalised measurements. Besides avoiding the inversion of poorly-conditioned matrices, the interpolation approach is less intensive computationally compared to matrix inversion, particularly when the matrix $\mathbf{G}$ (in (50)) is large. In our experimental results, we have used linear splines as the interpolator family $\{\gamma_{n}(\mathbf{x})\}$.

V The Sensor Network Algorithms: Source estimation from field samples

Based on the framework outlined in Sections III and IV, we now develop practical sensor network algorithms for estimating the sources of a physical field from its sensor measurements, and therefore solve the class of ISPs driven by linear PDEs with constant coefficients. First, we outline a centralised algorithm, which follows straightforwardly from Equation 23, then we consider the distributed case which relies on average consensus algorithms.

V-A Centralised source estimation

Assuming the source measurements have been made available at the fusion centre and that the PDE model—and Green’s function—of the monitored phenomena are known, then the point and instantaneous source estimation scheme can be summarised as in Algorithm 1, when the number of field sources $M$ is known.

Remark 1

Notice that Algorithm 1 requires $M$ as input, however when $M$ is not known, the following interesting phenomenon persists. Let ${M^{\prime}}$ be our initial guess for the unknown $M$, if it is incorrect (i.e. when ${M^{\prime}}\neq M$), we observe even in noisy settings, that:

1. (a)

If $M^{\prime}<M$, the most dominant $M^{\prime}$ sources are recovered from the field measurements. 2. (b)

Otherwise if $M^{\prime}>M$, all $M$ sources are recovered. In addition to them however, $M^{\prime}-M$ spurious sources—which may be attributed to noise and other sources of model mismatch—are also estimated. These extra $M^{\prime}-M$ spurious sources will either fall outside of the monitored region $\Omega$ or their estimated intensities will be very small in comparison to the true sources.

As we will see in Section V-C, this observation can be conveniently exploited to develop a suitable estimation scheme even when $M$ is not known a priori.

V-B Distributed source estimation

In the distributed set up, we want each node in the network to first estimate $\mathcal{Q}(\bm{k},r)$ through localised interactions with its neighbouring nodes. These localised interactions in our field estimation setting are based on the use of consensus algorithms.

We assume the same sensor network as described in Section II-A, comprising of “smart” sensor nodes. In addition to knowing the Green’s function of the monitored phenomena333This comes for free since the SNs are designed to sense a particular phenomena: i.e. if we are sensing acoustic fields then we use the wave equation, for temperature and leakages we use the diffusion equation, and so on., these nodes are able to perform mathematical computations and can also learn the network topology upon deployment. Knowledge of the network topology and the Green’s function of the underlying phenomena, means that all sensors can compute independently their sets of exponential reproducing coefficients $\{w_{n,l}(\bm{k},r)\}_{l}$ by using, either (35) if they are on a uniform grid or the least-squares scheme described in Section IV-C2 when they are nonuniformly placed.

They may start to sense the field in time to measure $\varphi_{n}(t_{l})$. To estimate the unknown sources the sensors must exchange and aggregate their sensor data, using average consensus as described in Appendix B. Specifically, gossiping is initiated when the $n$-th sensor contacts and exchanges its local measure,

[TABLE]

with a neighbour. After several rounds of gossip, each sensor will converge to the generalised measurements,

[TABLE]

Consequently $N$-D ESPRIT, or a similar Prony-like method, is then applied by each sensor—independently—on $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k}}$, to compute locally the unknown source parameters.

V-C Estimation in presence of noise and model mismatch

Although the sensors actually acquire noisy measurements, i.e. $\{\varphi_{n,l}^{\epsilon}\}_{n,l}$ in (12), of the underlying field the same recovery schemes outlined in Sections V-A and V-B are still effective. Specifically mapping the noisy sensor samples, $\{\varphi_{n,l}^{\epsilon}\}_{n,l}$, to generalised measurements using (18) gives the noisy sequence,

[TABLE]

Thus to recover the field sources from $\{\mathcal{Q}_{\epsilon}(\bm{k},r)\}_{\bm{k}=(0,\ldots,0)}^{(K_{1},\ldots,K_{d})}$, using $N$-D ESPRIT, we construct a (noisy) multilevel Hankel matrix $\mathbf{H}_{\epsilon}$ according to (58). If $M$ is known, then choosing $K_{i}\geq 2M-1$ for all $i=1,\ldots,d$, and retaining only those singular vectors due to the $M$ largest singular values of $\mathbf{H}_{\epsilon}$ is implicitly denoising. Generally, choosing large $K_{i}$ for all $i$ promotes robustness even in low SNR regimes.

Often however, the multilevel Hankel structure of $\mathbf{H}_{\epsilon}$ is lost when we retain only its $M$ largest singular values (and zero the rest). However we may restore its structure by averaging the appropriate elements. These steps can be repeated until convergence. This is the fundamental idea behind Cadzow-like [54] algorithms common in the FRI literature. Furthermore, these approaches are most effective when the noise in $\mathbf{H}_{\epsilon}$ is white, so we first need to apply a noise prewhitening transform similar to the approach in [43].

V-C1 Recovering an unknown number of sources

By combining the observations in Remark 1 with the fact that the singular value decomposition of the multilevel Hankel matrix $\mathbf{H}$ (or $\mathbf{H}_{\epsilon}$) also encodes information about $M$—i.e. most dominant singular values are due to the sources—an iterative estimation scheme similar to those in [37, 43] can be devised. The scheme relies on finding a time interval over which only a fixed number of field sources are active. Finding such a window allows us to reliably estimate these active sources and then adjust the spatiotemporal sensor measurements by removing their contribution to the sensor measurements. Given $\{\varphi_{n}(t_{l})\}_{l=0}^{L}$, the strategy (to find a window with one active source) is as follows:

1. (a)

Assume that there are $M^{\prime}\geq 2$ sources and approximate $\{\mathcal{Q}(\bm{k},0)\}_{\bm{k}}$ where $k_{i}=0,1,\ldots,2M^{\prime}-1$, using only the samples $\{\varphi_{n}(t_{l})\}_{l=0}^{L^{\prime}}$ with $L^{\prime}<L$ from the time window $[0,L^{\prime}\Delta_{t}]$. The centralised or distributed approaches (in Section V) can be used. 2. (b)

Proceed to estimate the $M^{\prime}$ source intensities $\{c^{\prime}_{m^{\prime}}\}_{m^{\prime}=1}^{M^{\prime}}$ and locations $\{\bm{\xi}^{\prime}_{m^{\prime}}\}_{m^{\prime}=1}^{M^{\prime}}$ using the multidimensional ESPRIT method. 3. (c)

Check that the estimated sources are valid. In particular, a source $(c^{\prime}_{m^{\prime}},\bm{\xi}^{\prime}_{m^{\prime}})$ is valid if the conditions below are simultaneously satisfied:

1. (a)

$c^{\prime}_{m^{\prime}}$ is greater than some predetermined threshold, and; 2. (b)

$\bm{\xi}^{\prime}_{m^{\prime}}$ is within the monitored region.

Check all $M^{\prime}$ sources and let $M_{\text{vs}}$ be the number of valid sources found. 4. (d)

There are three cases, if

1. i)

$M_{\text{vs}}>1$: Reduce the time widow by reducing $L^{\prime}$ and return to step (a). 2. ii)

$M_{\text{vs}}<1$: Increase the time widow by increasing $L^{\prime}$ and return to step (a). 3. iii)

$M_{\text{vs}}=1$: Estimate the source parameters from the measurements $\{\varphi_{n}(t_{l})\}_{l=0}^{L^{\prime}}$ using Algorithm 1. 5. (e)

Synthesise the field due to this source using equation (2) and adjust the sensor measurements by removing the contribution of this source. Increment $L^{\prime}$ and return to step (a). 6. (f)

Stop when the $L^{\prime}=L$ or when the field measurements are below some predefined threshold.

V-D Filtering in the time-domain

Using the framework summarised by Equation 23, we have been able to devise practical sensor network algorithms to solve the ISP of interest. During the sensing phase, although spatial prefiltering is generally not realisable, we are still able to perform prefiltering in time. The prefiltered samples obtained by the $n$-th sensor, using the filter $h(t)$ are:

[TABLE]

where $\star$ is the time-convolution operator. In light of this new formulation, the generator that will be used to reproduce exponentials from its space-time translates is

[TABLE]

with $G_{f}(\mathbf{s}_{\mathbf{x}},s_{t})={G}(-\mathbf{s}_{\mathbf{x}},-s_{t})H(s_{t}).$

We now have the freedom to design $H(s_{t})$ in such a way that $G_{f}(\mathbf{s}_{\mathbf{x}},s_{t})$ has some desirable properties. For our framework it is favourable to choose $H(s_{t})$ such that $G_{f}(\mathbf{s}_{\mathbf{x}},s_{t})$ decays quickly, at least in the $s_{t}$-domain. This will reduce the approximation error (33) as discussed in Section IV-B.

VI Numerical Simulations

In this section, we present some numerical results to validate the proposed framework for solving ISPs driven by linear PDEs. To investigate several scenarios, we present results for cases where the measured phenomena have been generated by:

1. (a)

the diffusion equation – we simulate a $2$-D multiple source field assuming a nonuniform sensor placement and investigate the proposed interpolation approach outlined in Section IV-C3. The results are shown in Figure 3. We also devote Section VI-D to comparing the performance and computational complexity of the proposed framework to a sparse synthesis recovery scheme [16]. 2. (b)

Poisson’s equation – we consider a single source field in $3$-D with nonuniform spatial samples, and utilise the linear system approach to find the desired analysis coefficients. Corresponding results are presented in Figure 4. 3. (c)

the wave equation – we simulate a single source wave field in $3$-D, sampled using a distributed SN of uniformly placed sensors. We assume that the sensors filter the field in time using a cubic spline before sampling. The results obtained are summarised in Figure 5.

The sensor measurements are simulated numerically and artificially corrupted by AWGN as defined in (13), using Matlab. We then apply our estimation algorithms on the measurements. For statistical significance, multiple independent trials are performed within each setup, with both a new noise process and random sensor placement (for the nonuniform sampling setups). In Figures 3, 4 and 5, the green ‘$\bullet$’, red ‘$\times$’ and blue ‘$+$’ denote the sensor locations, the estimated and true source locations, respectively.

VI-A Inverse source problem for the diffusion equation

We present numerical results in the nonuniform sampling scenario, where the diffusion field is induced by three localised and instantaneous sources. We use the $2$-D test function family $\{\Psi_{\bm{k}}(\mathbf{x})\Gamma_{r}(t)=e^{\mathrm{j}k_{1}x_{1}+\mathrm{j}k_{2}x_{2}}e^{\mathrm{j}rt/T}\}_{\bm{k},r}$ with $r=1,k_{1}=k_{2}=1,2,\ldots,15$ and present in Figure 3 the source estimation results obtained by using Algorithm 1 with a linear interpolator and a resampling grid with $\Delta_{x_{1}}=\Delta_{x_{2}}=1/30$. For statistical significance, we perform 20 independent trials of the experiment.

Moreover for the $2$-D field, with Green’s function (7), we obtain its Laplace transform as (see Appendix D-A):

[TABLE]

provided $\Re{\left(s_{t}-{\mu}\|\mathbf{s}_{\mathbf{x}}\|^{2}\right)}>0$, where $\Re{\left(z\right)}$ is used to denote the real part of a complex number $z$. Hence, by substituting (53) into (35) with $\bm{\kappa}=\mathrm{j}(\bm{k},r/T)$, the desired exponential reproducing coefficients are of the form:

[TABLE]

As seen in Figure 3, all source parameters have been recovered reliably; in particular, the source locations (Figure 3(a)) and activation times (Figure 3(b)), which are usually the two main parameters of interest, vary only marginally around their true values despite the low measurement SNR.

Remark 2

The requirement that $\Re{\left(s_{t}-{\mu}\|\mathbf{s}_{\mathbf{x}}\|^{2}\right)}>0$ is necessary for the transform integral to converge. By substituting $(\mathbf{s}_{\mathbf{x}},s_{t})=-\mathrm{j}(\bm{k},r/T)$ into (53), it is easy to see that this is satisfied when $\bm{k}\neq\{\bm{0}\}$.

VI-B Inverse source problem for Poisson’s Equation

Since the Poisson field is static, we only focus on recovering the source location and intensity. Specifically, we estimate the unknown source from a single time-snapshot of the $3$-D potential field measurements by formulating the linear system discussed in Section IV-C, with $r=0$ and $K_{1}=K_{2}=K_{3}=2$. Therefore according to (36), we get:

[TABLE]

where $I=1000$ and $\mathbf{x}^{\prime}_{i}$ is obtained from a lexicographic ordering of $\{(i_{1}\delta_{x_{1}}+\varepsilon,i_{2}\delta_{x_{2}}+\varepsilon,i_{3}\delta_{x_{3}}+\varepsilon)\}_{i_{1},i_{2},i_{3}=0}^{9}$ with $\delta_{x_{1}}=\delta_{x_{1}}=\delta_{x_{1}}=0.03$,444The slight shift $\varepsilon\geq 0$ is used here to avoid the singularity of $g(\mathbf{x})$ at $\|\mathbf{x}\|=0$. and $g(\mathbf{x})=-\frac{1}{4\pi\|\mathbf{x}\|}$.

The simulation results are shown in Figure 4, with the source location and intensity estimates obtained over the $20$ independent trials. All estimates show only small variations around the ground truth parameters. In particular, the variation in location estimates is much smaller than the average sensor spacing. Furthermore the condition number for the matrix $\mathbf{G}$ is reasonably small falling in the range $15$ – $38$.

VI-C Distributed acoustic source localisation: inverse source problem for the wave equation

Here the simulated physical phenomena is a $3$-D wave field, obeying (8) and, induced by a single source. We consider the distributed estimation setup, wherein the SN model is assumed to be the RGG $\mathcal{G}(N,0.4)$. Each of the $N=27$ regularly placed sensors acquires the field by filtering it (in time) using a third order B-spline before sampling. We apply the approach described in Section V-D by first computing the desired coefficients $w_{\bm{n},l}(\bm{k},r)$ from the filtered Green’s function $g_{f}(\mathbf{x},t)$, which has the bilateral Laplace transform, $G_{f}(\mathbf{s}_{\mathbf{x}},s_{t})=G(\mathbf{s}_{\mathbf{x}},s_{t})H(s_{t})$, where $H(s_{t})=\left(\frac{1-e^{-s_{t}}}{s_{t}}\right)^{4},$ is the Laplace transform for the third order B-spline555The zero order B-spline is taken as the indicator function on [0,1], from this we define the $n$-th order B-spline as the spline generated by convolving $(n+1)$ zero order B-splines., whereas $G(\mathbf{s}_{\mathbf{x}},s_{t})=\frac{1}{\|\mathbf{s}_{\mathbf{x}}\|^{2}-({s_{t}/c})^{2}}$. A derivation of this expression is provided in Appendix D-B. Combining these gives us the desired coefficients,

[TABLE]

where $\Delta_{x_{1}}=\Delta_{x_{2}}=\Delta_{x_{3}}=0.1$ and $\Delta_{t}=1$. Moreover, we choose $K_{1}=K_{2}=K_{3}=5$ and $r=1$ for the sensing function family $\{\Psi_{\bm{k}}(\mathbf{x})\Gamma_{r}(t)=e^{k_{1}x_{1}+\mathrm{j}k_{1}x_{2}+\mathrm{j}k_{2}x_{3}}e^{\mathrm{j}rt/T}\}_{\bm{k},r}$. Given these coefficients, the sensor nodes can perform a distributed estimation of the unknowns via gossiping, as described in Section V-B. The estimation results for 20 independent trials in this uniform sampling case is presented in Figure 5, wherein the recovered source locations, in (a), and activation times, in (b), are plotted over the true values. From these plots, we can conclude that the estimates are reliable, even when we have noisy measurements. In addition by considering a single experiment only, Figure 5(c)–(f) shows the evolution of the location and activation time estimates for three randomly chosen sensors (with each pairwise gossip round). Specifically, Figure 5(c), (d) and (e) displays the location estimates (in each spatial dimension), i.e. $x_{1},x_{2}$ and $x_{3}$ respectively. Similarly, Figure 5(f) shows the evolution of the activation time estimate. In both cases we notice that the estimates converge to the true value, despite the low measurement SNR.

VI-D Estimation performance comparison

We now numerically compare the proposed scheme against the sparse synthesis (SS) formulation described in [16]. For the SS method, we formulate the LASSO problem:

[TABLE]

where the dictionary $\mathbf{D}\in\mathbb{R}^{N(L+1)\times N_{x}N_{t}}$ is formed by discretising the Green’s function appropriately assuming uniform spatial and temporal grids with $N_{x}=51$ and $N_{t}=101$ divisions, $\mathbf{f}\in\mathbb{R}^{N(L+1)}$ is the source discretisation, whilst $\bm{\varphi}\in\mathbb{R}^{N_{x}N_{t}}$ is the vectorised sensor samples. The problem is then solved using the alternating direction method of multipliers (ADMM) [55]. The sparse synthesis formulation exhibits the same estimation performance as the sparse analysis formulation, proposed in [15, 16], so we only focus on one of them; although SS is much more computation intensive. The proposed sampling-based method uses $r=1$ and $K=50$. Note that for both approaches their activation time estimates are further refined by performing a line search in a small neighbourhood of the initial estimate.

The field used is $1$-D diffusion induced by a single source with parameters $c_{1}=5$, $\xi_{1}=0.1207m$ and $\tau_{1}=1.2175s$ over $x\in\Omega=[0,0.3]m$ and $t\in[0,20]s$. All reported statistics were computed from $200$ independent trials, each having a new noise realisation.

Figure 6 shows that for very low SNR combined with low samples $N=6$ and $L+1=6$, the sparsity-based recovery method achieves a lower mean absolute error (MAE) than the proposed method. However at higher temporal sampling frequency, for e.g. $L+1=21$, the current method consistently outperforms the SS method even at an extremely low SNR$=0$dB. Moreover, for fixed SNR $=20$dB we observe, in Figure 7, that the MAE remains constant even as the number of temporal samples increases – largely due to the discretisation; where as the current approach improves drastically with the number of temporal samples, outperforming the SS method for most values of $N$ and $L$. Finally by comparing the runtimes of both schemes in Figure 8, we notice that the complexity of the present method increases only marginally, whilst the SS method measures around one to two orders of magnitude slower as the number of temporal and spatial samples increase.

VII Conclusion

In this work we established a new general framework to solve ISPs for fields described by linear PDEs with constant coefficients. This formulation reduced the ISPs, for such physics-driven fields, to the problem of reproducing exponentials using shifted versions of the space- and time-reversed Green’s function of the corresponding field, under the assumption that the sources are localised. Consequently, we proposed practical sensor network algorithms for both uniform and nonuniform sampling setup. Finally, we validated our framework on some popular PDE models encountered in various applications, and also demonstrated that our method compares favourably against sparse recovery based methods.

Appendix A Estimating Parameters of Multidimensional Superimposed Exponentials

Consider the sequence (of generalised measurements)

[TABLE]

for any $r$, it is a superposition of $M$ $d$-dimensional damped complex sinusoids. Moreover, $k_{i}=0,1,\ldots,K_{i}-1$ for each $i=1,2,\ldots,d$. This multidimensional Prony-like system is prevalent in spectral estimation and array processing applications for example. Many approaches have been put forward for recovering the unknown frequencies and amplitudes, particularly for $d=2$ [40, 56, 41].

In our framework, discussed in Section III, solving this system for $\{b_{m},c_{m},(v_{i,m})_{i=1}^{d}\}_{m=1}^{M}$ allows us to recover the unknown source parameters $\{c_{m},\tau_{m},\bm{\xi}_{m}\}_{m=1}^{M}$. We utilise the $N$-D ESPRIT algorithm by Sahnoun et al [42] (based on the $2$-D ESPRIT algorithm of [56]) to solve this system. The $N$-D ESPRIT algorithm is as follows:

1. (a)

Choose $L_{1},L_{2},\ldots,L_{d}\in\mathbb{N}$ such that $1\leq L_{i}\leq K_{i}$ and set $J_{i}=K_{i}-L_{i}+1$. 2. (b)

Construct the multilevel Hankel matrix $\mathbf{H}$

[TABLE]

where the vectors $\mathbf{q}_{\bm{j}}{\stackrel{{\scriptstyle\text{\tiny{def}}}}{{=}}}\mathrm{vec}\mathopen{}\left(\mathcal{Q}(\bm{j}-\bm{1}:\bm{j}+\bm{L}-2\cdot\bm{1},r)\right)$, $\bm{j}{=}(j_{1},j_{2},\ldots,j_{d})\in\mathbb{N}^{d}_{+}$, and $\bm{L}=(L_{1},L_{2},\ldots,L_{d})$. Here $\bm{j}:\bm{j}+\bm{L}-\bm{1}\stackrel{{\scriptstyle\text{\tiny{def}}}}{{=}}(j_{1}:j_{1}+L_{1}-1,j_{2}:j_{2}+L_{2}-1,\ldots,j_{d}:j_{d}+L_{d}-1)$ is the (MATLAB-like) notation used for extracting subarrays, and $\mathrm{vec}\mathopen{}(\cdot\mathclose{)}$ is the vectorisation operator. 3. (c)

Retrieve the singular value decomposition (SVD) of $\mathbf{H}$, and form $\mathbf{U}\in\mathbb{C}^{(L_{1}L_{2}\cdots L_{d})\times M}$, which is a matrix of the $M$ most dominant left singular vectors of $\mathbf{H}$. 4. (d)

Compute the matrices $\mathbf{F}_{i}$, using

[TABLE]

where ${}_{i}\underline{\mathbf{U}}=\left.{}_{i}\underline{\mathbf{I}}\right.\mathbf{U}$ and ${}^{i}\overline{\mathbf{U}}=\left.{}^{i}\overline{\mathbf{I}}\right.\mathbf{U}$. Moreover, ${}_{i}\underline{\mathbf{I}}=\mathbf{I}_{\prod_{j=1}^{i-1}L_{j}}\otimes\underline{\mathbf{I}}_{L_{i}}\otimes\mathbf{I}_{\prod_{j=i+1}^{d}L_{j}}$ and ${}^{i}\overline{\mathbf{I}}=\mathbf{I}_{\prod_{j=1}^{i-1}L_{j}}\otimes\overline{\mathbf{I}}_{L_{i}}\otimes\mathbf{I}_{\prod_{j=i+1}^{d}L_{j}}$. The symbol $\otimes$ denotes the Kronecker product, $\mathbf{I}_{n}$ is the $n\times n$ identity matrix and the overbar (respectively underbar) represents the operation of deleting the first (respectively last) row of a matrix. 5. (e)

For some random choice of $\beta_{1},\beta_{2},\ldots,\beta_{d}$, compute the linear combination of matrices:

[TABLE] 6. (f)

Diagonalise the matrix $\mathbf{K}$ to find $\mathbf{T}$, such that

[TABLE] 7. (g)

Then for each $i=1,2,\ldots,d$, transform $\mathbf{F}_{i}$ using the matrix $\mathbf{T}$,

[TABLE] 8. (h)

Compute the unknown frequencies $\{v_{i,1},v_{i,2},\ldots,v_{i,M}\}$ as the diagonal of $\mathbf{D}_{i}$, i.e. $\mathrm{diag}\mathopen{}\left(\mathbf{D}_{i}\right)$ for each $i=1,2,\ldots,d$. 9. (i)

Using $\{v_{i,1},v_{i,2},\ldots,v_{i,M}\}_{i=1}^{d}$, solve the least-squares problem associated with (57) to find $\{c_{m}^{\prime}=c_{m}b_{m}^{r}\}_{m=1}^{M}$.

Appendix B Gossiping for Distributed Average Consensus

The problem of achieving consensus or agreement amongst agents of a network in a distributed fashion, is well-studied. Gossip algorithms, based on the early works of Tsitsiklis et al [57], have been applied to the distributed averaging problem, as they possess the attractive property of not requiring a specialised routing strategy. Consider for example pairwise randomised gossip [58], which involves at each time step two random but connected nodes updating their values with a weighted average of their current values.

Let us denote the value at node $n$ after the $i$-th pairwise gossip by $y_{n,i}$, then the initial value is $y_{n,0}$. At each iteration, a random node $n$ wakes up and contacts a randomly chosen neighbor $n^{\prime}$, they both update their estimates with $y_{n,i{+}1}=y_{n^{\prime},i{+}1}=(y_{n,i}{+}y_{n^{\prime},i})/2$. Furthermore, let $\mathbf{y}(i)=[y_{1,i},y_{2,i},\ldots,y_{N,i}]^{\!\mathsf{T}}$ then this pairwise gossip algorithm can be summarised mathematically as:

[TABLE]

where $\mathbf{P}(i)$’s are doubly stochastic matrices drawn randomly at the $i$-th iteration. Moreover, $\mathbf{P}(i)$ is a diagonal identity matrix everywhere else apart from the elements $(n,n),(n,n^{\prime}),(n^{\prime},n)$ and $(n^{\prime},n^{\prime})$ which are all equal to $1/2$.

If the network is connected and the nodes communicate using this scheme, then it can be shown that each node is guaranteed to converge to the global network average $\bar{y}=\frac{1}{N}\sum_{n=1}^{N}y_{n,0}$, after enough iterations, i.e. $\lim_{i\rightarrow\infty}\mathbf{y}(i)=\mathbf{1}\bar{y}$. Performance guarantees and convergence results have also been studied (see [58] and the references therein).

The localised interactions in our inverse source problem set up will be based on the use of gossip algorithms to compute the generalised measurements $\{\mathcal{Q}(\bm{k},r)\}_{\bm{k}}=\sum_{n,l}w_{n,l}(\bm{k},r)\varphi_{n}(t_{l})$, in a distributed manner. Specifically, we assume each node $n$ knows only its sequence of coefficients $\{w_{n,l}(k,r)\}_{l}$ and senses the field $\varphi_{n}(t_{l})$ in time, then the sensors communicate their local measures $y_{n}(\bm{k},r)=N\sum_{l}w_{n,l}(\bm{k},r)\varphi_{n}(t_{l})$. On convergence, the nodes will have $\mathcal{Q}(\bm{k},r)=\frac{1}{N}\sum_{n=1}^{N}y_{n}(\bm{k},r),$ which is the desired generalised measurements.

Appendix C The Generalised Strang-Fix Conditions

We first state the classical Strang-Fix condition [49], for multidimensional polynomial reproduction.

Lemma 2** (Strang-Fix condition [49, 38])**

Any compactly supported kernel $\psi(\mathbf{x})$ whose derivatives up to and including order $p$ are in $\mathcal{L}^{2}$, is able to reproduce polynomials, i.e.:

[TABLE]

where $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})$ with $\sum_{i}\alpha_{i}\leq p$, if and only if

[TABLE]

where $\hat{\psi}(\mathrm{j}\bm{\omega})$ is the Fourier transform of $\psi$, and $\bm{\ell}\in\mathbb{Z}^{d}\setminus\{\bm{0}\}$.

Proof:

The proof of this lemma depends on the Poisson summation formula – which connects the values of a function $\psi$ on the lattice $\mathbb{Z}^{d}$ with its Fourier transform $\hat{\psi}$ on the lattice $2\pi\mathbb{Z}^{d}$. Note that both sides of the summation converge (absolutely) when $\psi$ and $\hat{\psi}$ decay sufficiently quickly. The proof is based on taking the Fourier transform of $\psi(\bm{n})=\bm{n}^{\bm{\alpha}}\varphi(\mathbf{x}-\bm{n})$ and applying the frequency differentiation property to it. The reverse implication follows immediately from this expression, whereas the forward follows by induction. A complete proof appears in **[49]**. ∎

To demonstrate the multidimensional generalised Strang-Fix condition, we adapt the proof of [38]. First, we require that the function $\psi(\mathbf{x})=e^{-\bm{\kappa}\cdot\mathbf{x}}s(\mathbf{x})$ is able to reproduce the polynomial $\mathbf{x}^{\bm{\alpha}}$. Consequently, we obtain

[TABLE]

such that $w_{\bm{n}}=\left.c_{\bm{\alpha},\bm{n}}e^{\bm{\kappa}\cdot\bm{n}}\right\rvert_{\bm{\alpha}=\bm{0}}$. Moreover, it follows that (65) holds true provided $\psi(\mathbf{x})$ satisfies the classical Strang-Fix condition with $\bm{\alpha}=\bm{0}$. This is the case when $\hat{\psi}(\bm{0})\neq\bm{0}$ and $\hat{\psi}(2\pi)=\bm{0}$. Consequently, since the Fourier transform of $\psi$ is by construction related to the Laplace transform of $s(\mathbf{x})$, by $\hat{\psi}(\bm{\omega})=\left.S(\bm{\kappa}+\mathbf{s}_{\mathbf{x}})\right\rvert_{\mathbf{s}_{\mathbf{x}}=\mathrm{j}\bm{\omega}}$, we obtain the desired conditions:

[TABLE]

Appendix D Bilateral Laplace transforms

D-A Green’s function of the two-dimensional diffusion equation

Herein, we derive expression (53), i.e. the multidimensional bilateral Laplace transform of the diffusion Green’s function, $g(\mathbf{x},t)=\frac{1}{(4\pi\mu t)^{\frac{d}{2}}}e^{-\frac{\|\mathbf{x}\|^{2}}{4\mu t}}H(t).$ Thus,

[TABLE]

We first consider the spatial integral:

[TABLE]

where the second equality follows by completing the square and (68) uses the fact that $\int_{x=-\infty}^{\infty}e^{-\frac{x^{2}}{a}}\mathrm{d}x=\sqrt{a\pi},$ thus,

[TABLE]

with $a=4\mu t$. We can now substitute (68) back into (67) and proceed as follows:

[TABLE]

provided $\Re{\left(s_{t}-{\mu}\|\mathbf{s}_{\mathbf{x}}\|^{2}\right)}>0$, as required.

D-B Green’s function of the three-dimensional wave equation

We now obtain the multidimensional bilateral Laplace transform of the Green’s function for the wave equation (8), which by definition must satisfy,

[TABLE]

We begin by taking the Laplace transform of the PDE above, as follows:

[TABLE]

[TABLE]

We can now rearrange to obtain the Laplace transform of $g$,

[TABLE]