Spatial resolution of the PBP algorithm in acoustical inverse scattering

F. Natterer Institut für Numerische und instrumentelle Mathematik

August 1997

1 Introduction

The PBP algorithm [1] solves the inverse scattering problem

Δu + k^{2}(1 + f)u = 0 in IR^{2} ,

(1.1)

u = e^{ikx·θ}(1 + v)

where e^{ikx·θ}v satisfies the radiation condition, θ ∈ S^{1}, and k > 0 fixed. The sought - for function f is supported in |x| ≤ ρ, and the data is g_{+} = v on Γ_{+}(θ), g_{−} = v on Γ_{−}(θ) where

Γ_{±}(θ) = {x = ±ρθ + sθ^{⊥} : s ∈ IR^{1}} .

Hence g_{+} is the transmitted signal and g_{−} the backscatter for incidence direction θ. The spatial resolution we expect is λ/2 with λ = 2π/k the wavelength of the irradiating wave. This means that functions f of (essential) bandwidth 2k should be recovered reliably. It is well known [2] that this is true within the Born approximation, where the frequency components < \sqrt{2}k are determined by the transmitted signal g_{+}(θ, s), and those between \sqrt{2}k and 2k by the backscatter g_{−}(θ, s). We note that the condition for the Born approximation to hold is

|Rf| ≪ \frac{2π}{k}

where R is the Radon transform. This is a slight extension of the condition in [2].

In this paper we analyse PBP within the Born approximation. We show that the algorithm as it stands can’t possibly achieve resolution λ/2 or 2k. Based on our analysis we suggest a modificantion of PBP and demonstrate by numerical examples taken from [9] that the modified PBP algorithm does achieve resolution 2k even in the fully nonlinear case, i.e. when the Born approximation is not valid. We do not aim at superresolution [7].

2 The PBP algorithm

The PBP algorithm is essentially an extension of the ART algorithm of CT (see e.g. [3]) to the nonlinear System

R_{j}(f) = g_{+}(θ_{j}, ·) ,  j = 1, …, p  (2.1)

where R_{j}(f) = v^{j} on Γ_{+}(θ_{j}) with v^{j} a low-pass filtered version of the solution v of (1.1) with θ = θ_{j} and initial values

v = g_{−}(θ_{j}, ·) ,  \frac{∂v}{∂ν} = \frac{∂g_{−}}{∂ν} (θ_{j}, ·) on Γ_{−}(θ_{j}) .  (2.2)

Here, ν denotes the interior normal on Γ_{−}(θ_{j}). It has been shown in [4] that v^{j} can be determined in a stable way from these initial values provided the cut-off frequency of the low-pass filter is put at some value slightly smaller than k. We remark that ∂g_{−}/∂ν can be computed easily from g_{−} by solving the exterior Helmholtz problem where f = 0. We also remark that (2.1) involves implicitly the backscatter g_{−} via the initial conditions (2.2) even though g_{−} does not show up explicitly.

PBP solves (2.1) iteratively. We describe the step f^{0} → f^{1} where f^{0} is an initial approximation:

f_{0} = f^{0}

For j = 1, …, p :

f_{j} = f_{j−1} − ωR_{j}^{′^{∗}}(f_{j−1})S_{j}^{−1}(R_{j}(f_{j−1}) − g_{j}) ,  (2.3)
f^{1} = f_{p} .

Here we have put g_{j} = g_{+}(θ_{j}, ·). S_{j} is a certain operator on L_{2}(Γ_{+}(θ_{j})) which is supposed to mimic the action of R′_{j}(f_{j−1})R_{j}^{′^{∗}}(f_{j−1}). In [5] we have shown that for linear operators R_{j} the process converges for 0 < ω < 2 provided that S_{j} ≥ R_{j}R_{j}^{∗}. In our implementation of PBP we use either a multiple of the unit operator or the operator R′_{j}(0)R_{j}^{′^{∗}}(0). The latter choice leads to convergence in a few steps provided the initial approximation is chosen reasonably.

We now give an analysis of PBP in the Born appximation. The algorithm is then essentially identical to the one given in [6], except that R_{j}(f_{j−1}) is evaluated in a different way. The differential equation for v is

Δv + 2ikθ · ∇v + k^{2}vf = −k^{2}f .  (2.4)

The Born approximation is obtained by dropping k^{2}vf, yielding

Δv + 2ikθ · ∇v = −k^{2}f .  (2.5)

The initial value problem for this equation with initial values on Γ_{−}(θ) is easily solved by taking the Fourier transform of v(tθ + sθ^{⊥}) with respect to s. Putting κ = \sqrt{k^{2} − σ^{2}} we obtain for |σ| < k

v̂(tθ + σθ^{⊥}) = −\frac{k^{2}}{κ} ∫_{−ρ}^{t} e^{−ik(t−t′)} sin(κ(t − t′))f̂(t′θ + σθ^{⊥})dt′

(2.6)

+e^{−ik(t+ρ)} (ĝ_{−}(θ, σ) cos(κ(t + ρ)) + \frac{1}{κ}(\frac{∂ĝ_{−}}{∂ν} + ikĝ_{−})(θ, σ) sin(κ(t + ρ))) .

For t = ρ and θ = θ_{j} the left hand side is just v̂^{j} on Γ_{+}(θ_{j}). Hence

(R_{j}f)^{∧}(σ) = −\frac{k^{2}}{κ} ∫_{−ρ}^{ρ} e^{−ik(ρ−t)} sin(κ(ρ − t))f̂(tθ + σθ^{⊥})dt

+ Ĝ_{−}(θ_{j}, σ)

where G_{−} depends only on g_{−}, ∂g_{−}/∂ν on Γ_{−}(θ_{j}) but not on f. Expressing f̂ by the 2D Fourier transform of f which we denote by f̂, too, yields

(R_{j}f)^{∧}(σ) = −\frac{k^{2}}{(2π)^{1/2}κ} ∫_{IR^{1}} ∫_{IR^{1}} e^{−ik(ρ−t)} sin(κ(ρ − t))f̂(τθ_{j} + σθ_{j}^{⊥})e^{itτ}dτdt
+ Ĝ_{−}(θ_{j}, σ) .

The integral over [−ρ, ρ] has been extended to all of IR^{1} because f(x) vanishes outside |x| < ρ. The t integral can be evaluated, yielding

(R_{j}f)^{∧}(σ) = −\frac{e^{−ikρ}(2π)^{1/2}}{2iκ} ∫_{IR^{1}} (e^{iκρ}δ(k − κ + τ) − e^{−iκρ}δ(k + κ + τ))f̂(τθ_{j} + σθ_{j}^{⊥})dτ

+Ĝ_{−}(θ_{j}, σ)

(2.7)

= −\frac{e^{−ikρ}(2π)^{1/2}}{2iκ} (e^{iκρ}f̂((κ − k)θ_{j} + σθ_{j}^{⊥}) − e^{−iκρ}f̂((−κ − k)θ_{j} + σθ_{j}^{⊥}))

+Ĝ_{−}(θ_{j}, σ) .

The values of f̂ which appear on the right hand side are just the values on the Ewald Spheres (or better Ewald circles in 2D), see Fig. 1.

θ −kθ σ ξ_{−}  ξ_{+}

Fig. 1: The points ξ_{±}(θ, σ) = (±κ − k)θ + σθ^{⊥} in Fourier space where f̂ is taken in (2.7).

From (2.7) we draw an important conclusion: Our operator R_{j} provides information only about a linear combination of f̂((±κ − k)θ_{j} + σθ_{j}^{⊥}). This is in contrast to diffraction tomography where each of these two numbers is determined separately by the transmitted and backscattered fields, resp.. It is clear that f̂ can’t be determined uniquely from (2.7). For instance, a function f̂ with

e^{iκρ}f̂(ξ_{+}(θ_{j}, σ)) = e^{−iκρ}f̂(ξ_{−}(θ_{j}, σ))

does not contribute to R_{j}(f). Such a function is easily constructed as a radially symmetric function.

The conclusion is that the operator R_{j} is not a good one for the solution of the inverse problem. The only reason for this choice is that R_{j}(f) can be evaluated very efficiently by initial value techniques. In the next section we show how to extend R_{j} to an operator which still can be evaluated by initial value techniques but which, at least in the Born approximation, is stably invertible.

3 The extended PBP algorithm

Our aim is to separate the two contributions of f̂ in (2.7). An obvious way to do this is to include in R_{j}(f) the derivative of (2.7) with respect to ρ. It will turn out that

R_{j}(f) = ( \begin{matrix}v^{j} \\ \frac{∂}{∂ν} v^{j} − ikv^{j} \end{matrix} )_{Γ_{+}(θ_{j})}  (3.1)

is a good choice. Here, v^{j} is the same function as in (2.1). The system we have to solve is now

R_{j}(f) = ( \begin{matrix}g_{+}(θ_{j}, ·) \\ (\frac{∂g_{+}}{∂ν} − ikg_{+}) (θ_{j}, ·) \end{matrix} ) ,  j = 1, …, p .  (3.2)

Again ∂g_{+}/∂ν can be computed easily by solving the exterior Helmholtz equation.

We compute the operator R′_{j}(f) : L_{2}(|x| < ρ) → L_{2}(Γ_{+}(θ_{j}))^{2} and its adjoint. Replacing v by v + w, f by f + h in (2.4) and neglecting wh we obtain

Δw + 2ikθ · w + k^{2}fw = −k^{2}(1 + v)h in IR^{2} ,

(3.3)

w = 0 , \frac{∂w}{∂ν} = 0 on Γ_{−}(θ) .

Let w^{j} be a low-pass filtered version of w for θ = θ_{j}. Then,

R′_{j}(f)h = ( \begin{matrix}w^{j} \\ \frac{∂}{∂ν} w^{j} − ikw^{j} \end{matrix} )_{Γ_{+}(θ_{j})} .  (3.4)

In order to determine R′_{j}(f)^{∗} : L_{2}(Γ_{+}(θ_{j}))^{2} → L_{2}(|x| < ρ) we take

g = ( \begin{matrix}g_{0} \\ g_{1} \end{matrix} ) ∈ L_{2}(Γ_{+}(θ_{j}))^{2} and compute

(R′_{j}(f)h, g) = (w^{j}, g_{0})_{L_{2}(Γ_{+}(θ_{j}))} + (\frac{∂w^{j}}{∂ν} − ikw^{j}, g_{1})_{L_{2}(Γ_{+}(θ_{j}))} .

Let z^{j} be a low-pass filtered version of the solution of the initial (better: final) value problem

Δz + 2ikθ · ∇z + k^{2}\overline{f}z = 0 in IR^{2} ,

(3.5)

z = −g_{1} , \frac{∂z}{∂ν} − ikz = g_{0} on Γ_{+}(θ)

for θ = θ_{j}. Let Q(θ) be the strip between Γ_{−}(θ), Γ_{+}(θ). Green’s second formula and an integration by parts yields for any w, z

∫_{Q(θ)} \{(Δw + 2ikθ · ∇w − k^{2}fw)\overline{z} − w\overline{(Δz + 2ikθ · ∇z − k^{2}\overline{f}z)}\} dx

= ∫_{Γ_{+}(θ)∪Γ_{−}(θ)} (w\frac{∂\overline{z}}{∂ν} − \frac{∂w}{∂ν}\overline{z}) ds + 2ik ( ∫_{Γ_{+}(θ)} w\overline{z}ds − ∫_{Γ_{−}(θ)} w\overline{z}ds) .

For w = w^{j}, z = z^{j}, θ = θ_{j} this reads

−k^{2} ∫_{Q(θ)} (1 + v)h\overline{z}dx = ∫_{Γ_{+}(θ)} \{w\overline{(\frac{∂z}{∂ν} − ikz)} − (\frac{∂w}{∂ν} − ikw) \overline{z}\} ds ,

hence

(h, −k^{2}(1 + \overline{v})z)_{L_{2}(|x|<ρ)} = (w, g_{0})_{L_{2}(Γ_{+}(θ))} + (\frac{∂w}{∂ν} − ikw, g_{1})_{L_{2}(Γ_{+}(θ))} .

This means that

R′_{j}(f)^{∗}g = −k^{2}(1 + \overline{v})z^{j} .  (3.6)

We may convince ourselves that R_{j}(f) from (3.1) contains, at least within the Born approximation, enough information about f to determine it up to spatial frequencies 2k. For this we look at the operator R′_{j}(0). Solving (3.3) for f = 0 yields according to (2.6)

ŵ(tθ + σθ^{⊥}) = −\frac{k^{2}}{κ} ∫_{−ρ}^{t} e^{−ik(t−t′)} sin(κ(t − t′))ĥ(t′θ + σθ^{⊥})dt′ .

As in (2.6), ŵ, ĥ denote the 1D Fourier transforms of w(tθ+sθ^{⊥}), h(t′θ+sθ^{⊥}) with respect to s. Observing that \frac{∂w}{∂ν} = −\frac{∂w}{∂t} on Γ_{+}(θ) yields

(R′_{j}(0)h)^{∧}(σ) = k^{2} ∫_{−ρ}^{ρ} e^{−ik(ρ−t)} ( \begin{matrix}\frac{−1}{κ} sin κ(ρ − t) \\ cos κ(ρ − t) \end{matrix} ) ĥ(tθ + σθ^{⊥})dt . (3.7)

Inserting here the 2D Fourier transform of h, again denoted by ĥ, we obtain in a way completely analogous to (2.7)

(R′_{j}(0)h)^{∧}(σ) = k^{2} \frac{1}{2}(2π)^{1/2}e^{−ikρ} ·

( \begin{matrix}\frac{i}{κ}(e^{iκρ}ĥ((κ − k)θ_{j} + σθ_{j}^{⊥}) & − & e^{−iκρ}ĥ((−κ − k)θ_{j} + σθ_{j}^{⊥})) \\ e^{iκρ}ĥ((κ − k)θ_{j} + σθ_{j}^{⊥}) & + & e^{−iκρ}ĥ((−κ − k)θ_{j} + σθ_{j}^{⊥}) \end{matrix} ) .

Let K : L_{2}(Γ_{+}(θ)) → L_{2}(Γ_{+}(θ)) be defined by (Kg)^{∧} = κĝ. Then,

((±iK, I)R′_{j}(0)h)^{∧}(σ) = (2π)^{1/2}k^{2}e^{−i(k±κ)ρ}ĥ((∓κ − k)θ_{j} + σθ_{j}^{⊥}) . (3.8)

This shows that R′_{j}(0)h determines ĥ(ξ) for |ξ| ≤ 2k, with + corresponding to the backscatter and − to the transmitted field.

4 Implementation of the extended PBP method

There are two ways of implementing ART for (3.2). In the first one we use exactly (2.3) with R_{j} as in (3.1). For S_{j} we take the operator R′_{j}(0)R′_{j}(0)^{∗}. For R′_{j}(0) we already found in (3.7) a handy expression. R′_{j}(0)^{∗} is determined essentially in the same way. Solving (3.5) for f = 0 by taking Fourier transforms in Γ_{+}(θ) yields

ẑ(tθ + σθ^{⊥}) = (ĝ_{0}(σ)\frac{sin(κ(ρ − t))}{κ} − ĝ_{1}(σ) cos(κ(ρ − t))) e^{ik(ρ−t)} ,

hence, according to (3.6),

(R′_{j}(0)^{∗}g)^{∧}(tθ + σθ^{⊥}) =

−k^{2} (ĝ_{0}(σ)\frac{sin(κ(ρ − t))}{κ} − ĝ_{1}(σ) cos(κ(ρ − t))) e^{ik(ρ−t)} .

From this and (3.7) we obtain for the operator S_{j} in Fourier space

(S_{j}g)^{∧}(σ) = (R′_{j}(0)R_{j}^{′^{∗}}(0)g)^{∧}(σ)

= −k^{4} ∫_{−ρ}^{ρ} ( \begin{matrix}−\frac{1}{κ} sin(κ(ρ − t)) \\ cos(κ(ρ − t)) \end{matrix} ) (ĝ_{0}(σ)\frac{sin(κ(ρ − t))}{κ} − ĝ_{1}(σ) cos(κ(ρ − t))) dt

= k^{4} ( \begin{matrix}c_{1}/κ^{2} & , & c_{2} \\ c_{2} & , & c_{3} \end{matrix} ) ( \begin{matrix}ĝ_{0}(σ) \\ ĝ_{1}(σ) \end{matrix} ) ,  (4.1)

c_{1} = ∫_{−ρ}^{ρ} sin^{2}(κ(ρ − t))dt , c_{3} = ∫_{−ρ}^{ρ} cos^{2}(κ(ρ − t))dt

c_{2} = \frac{−1}{κ} ∫_{−ρ}^{ρ} cos(κ(ρ − t)) sin(κ(ρ − t))dt

Introducing α = 2ρκ be we obtain

c_{1} = ρ (1 − \frac{cos α sin α}{α}) , c_{3} = ρ (1 + \frac{cos α sin α}{α}) ,

(4.2)

c_{2} = −2ρ^{2} (\frac{sin α}{α})^{2} .

The operator S_{j}^{−1} in (2.3) in Fourier space is the inverse of (4.1), i.e.

S_{j}^{−1} = \frac{1}{Δk^{4}} ( \begin{matrix}c_{3} & −c_{2} \\ −c_{2} & c_{1}/κ^{2} \end{matrix} ) ,  (4.3)
Δ = \frac{c_{1}c_{3}}{κ^{2}} − c_{2}^{2} = \frac{ρ^{2}}{κ^{2}} (1 − (\frac{sin α}{κ})^{2}) .

With (4.3), (3.6) all the operations in (2.3) are known. Alternatively we may work with the operator

K R_{j} , K = ( \begin{matrix}iK & I \\ −iK & I \end{matrix} )

instead of R_{j}. Then, (2.3) reads

f_{j} = f_{j−1} − ωR_{j}^{′^{∗}}(f_{j−1})K^{∗}(KS_{j}K^{∗})^{−1}K(R_{j}f_{j−1} − g_{j})

with S_{j} as above. From (4.1) we obtain

(KS_{j}K^{∗}g)^{∧} = \binom{iκ , 1}{−iκ 1} k^{4} ( \begin{matrix}c_{1}/κ^{2} & , & c_{2} \\ c_{2} & , & c_{3} \end{matrix} ) ( \begin{matrix}−iκ & iκ \\ 1 & 1 \end{matrix} ) ĝ
= 2ρk^{4} ( \begin{matrix}1 & , & c \\ \overline{c} & , & 1 \end{matrix} ) ĝ , c = \frac{sin α}{α}e^{−iα} .

Hence the operator (KS_{j}K^{∗})^{−1} in Fourier space is

\frac{1}{2ρk^{4}} ( \begin{matrix}1 & , & c \\ \overline{c} & , & 1 \end{matrix} )^{−1} = \frac{1}{2ρk^{4}(1 − c\overline{c})} ( \begin{matrix}1 & , & −c \\ −\overline{c} & , & 1 \end{matrix} ) .

In the second implementation we consider the high and low frequency parts R_{j,±}(f) = (±iK, I)R_{j}(f) separately, i.e. we apply ART to the systems of 2p equations

R_{j,±}(f) = g_{j,±} ,  j = 1, …, p  (4.4)

where

g_{j,±} = ±iKg_{+}(θ_{j}, ·) + (\frac{∂g_{+}}{∂ν} − ikg_{+}) (θ_{j}, ·) .

The adjoint of R′_{j,±}(f) is

R′_{j,±}(f)^{∗} = R′_{j}(f)^{∗} ( \begin{matrix}∓iK \\ I \end{matrix} ) .

Hence the ART update for (4.4) is

f_{j} = f_{j−1} − ωR_{j}^{′^{∗}}(f_{j−1}) ( \begin{matrix}∓iK \\ I \end{matrix} ) S_{j,±}^{−1}(R_{j,±}(f_{j−1}) − g_{j,±}) .

A good choice for S_{j,±} is

S_{j,±} = R′_{j,±}(0)R_{j,±}^{′∗}(0)
= (±iK, I)S_{j} ( \begin{matrix}∓iK \\ I \end{matrix} )

with S_{j} from (4.1). In Fourier space this operator is

(S_{j,±}g)^{∧}(σ) = (∓iκ, 1)k^{4} ( \begin{matrix}c_{1}/κ^{2} & , & c_{2} \\ c_{2} & , & c_{3} \end{matrix} ) ( \begin{matrix}∓iκ \\ 1 \end{matrix} ) ĝ(σ)
= 2ρk^{4}ĝ(σ)

as is easily seen from (4.2). Thus the update simply reads

f_{j} = f_{j−1} − \frac{ω}{2ρk^{4}} R′_{j}(f_{j−1})^{∗} ( \begin{matrix}∓iK \\ I \end{matrix} ) (R_{j,±}(f_{j−1}) − g_{j,±}) .

Note that this iteration is more time consuming than the first one since R′_{j}(f_{j−1})^{∗} has to be evaluated twice as often for one complete sweep.

5 Numerical example

We consider the function

f(x) = 0.2 (1 + cos(11x_{1}) + sin(11x_{2})) , |x_{i}| ≤ π

which has been reconstructed in [9] succesfully from the scattering operator of that paper at frequencies k = 1, 2, …, 11. According to the theoretical resolution (see §1) f should be computable from the scattered field for some k with 2k ≤ 11. In fact the reconstruction in [9] for k ≤ 6 is almost perfect. We note that the condition for the Born approximation to be valid is |Rf| ≪ 2π/k. This means that approximately 0.2 · 2π ≪ 2π/k or k ≪ 5. Thus the Born approximation can only be used for k = 1, 2, 3, say, but not for k ≥ 4.

We did the reconstruction using the first version of the extended PBP algorithm with k = 6.1. This is only slightly larger than the smallest value 5.5 for which f can be recovered. The result is displayed in Fig. 1. The reconstruction after 6 sweeps with ω = 0.2 is almost perfect (with the exception of a Gibbs phenomenon close to the boundary of the object, and some disturbances in the vanishing imaginary part). As initial approximation we used the function f without the cosine and sine terms. The computing time on a SUN ultra 1 was 1 minute per sweep. We used 128 directions and we worked on a 129 × 129 grid corresponding to a square of sidelength 10.

Acknowledgement: We want to thank V. Palamodov for an analysis of the null space of (2.7) and N. Grinberg for the idea to include the normal derivative in the operator R_{j} from (3.2).

Figure Captions

Top left  : Original on a 129 × 129 grid.
Middle left : Initial approximation.
Top right  : Reconstruction after 6 sweeps.
Middle right : Reconstruction after 18 sweeps.
Bottom  : Horizontal cross sections.

Real part left, imaginary part right.

Literatur

[1] Natterer, F. - Wübbeling, F.: A propagation - backpropagation method for ultrasound tomography, Inverse Problems 11, 1225-1232 (1995).

[2] Kak, A.C. - Slaney, M.: Principle of Computerized Tomography, IEEE Press, New York 1987.

[3] Natterer, F.: The Mathematics of Computerized Tomography. Wiley-Teubner 1986.

[4] Natterer, F.: An initial value approach to the inverse Helmholtz problem at fixed frequency, in: Engl, H.W. - Louis, A.K. - Rundell, W. (eds.): Inverse Problems in Medical Imaging and Nondestructive Testing, Springer 1997.

[5] Natterer, F.: Numerical solution of bilinear inverse problems, Preprint 19/96-N, Fachbereich Mathematik, University of Muenster, 1996.

[6] Ladas, K.T. - Devaney, A.J.: Generalized ART algorithm for diffraction tomography, Inverse Problems 7, 109-125 (1991).

[7] Schatzberg, A. - Devaney, A.J.: Superresolution in diffraction tomography, Inverse Problems 8, 149-164 (1992).

[8] Natterer, F.: Finite difference methods for inverse problems, in: Ang, D.D. et al. (eds.): Inverse Problems and Applications to Geophysics, Industry, Medicine and Technology, Publications of the HoChiMinh City Mathematical Society 2, 1995.

[9] Chen, Y.: Inverse scattering via Heisenberg’s uncertainty principle, Inverse Problems 13, 253-282 (1997).