A new iterative time-reversal algorithm capable of identifying and focusing on multiple scatterers in a relatively small number of iterations is developed. It is recognized that the traditional iterated time-reversal method is based on utilizing power iterations to determine the dominant eigenpairs of the time-reversal operator. The convergence properties of these iterations are known to be suboptimal. Motivated by this, a new method based on Lanczos iterations is developed. In several illustrative examples it is demonstrated that for the same number of transmitted and received signals, the Lanczos iterations based approach is substantially more accurate.

## I. Introduction

Time-reversal methods are often used to target inhomogeneities or image scatterers in an acoustic medium.^{1,2} It is being actively pursued in medical ultrasonic imaging, in medical treatment via ultrasound, in underwater communication, in nondestructive evaluation of materials and structures using both linear and nonlinear elastic waves, in imaging the earth’s crust, and imaging its oceans.^{3–7,2}

The basic idea in time-reversal may be explained in the context of an array of transmitters and receivers embedded in an inhomogeneous medium. A source emits a pulse which is multiply scattered, dispersed, and eventually detected by the receivers. The received signal is then reversed in time and transmitted back into the medium. The reversed signal retraces its multiple paths, and medium dispersion recompresses the signal, so that the pulse later focuses at the original source point. In *iterated time-reversal* this process of transmission, time-reversal, and retransmission is repeated and with each iteration the transmitted signal selectively focuses on the strongest scatterer in the medium.

Iterated time-reversal can be interpreted as the power iteration method applied to the time-reversal operator $H$, which is related to the scattering operator $G$ through $H=GGT$, where $GT$ is the transpose of $G$ in the time domain (see below for precise definitions). An element $Gij$ of the discrete scattering operator represents the scattered signal measured at the $ith$ element of an array due to a unit impulse transmitted by the $jth$ element. Each iteration in the time-reversal method corresponds to one iteration of the power method; that is the evaluation of the product $Hv$ for a signal $v$. The power iterations converge to the eigenvector corresponding to the largest eigenvalue of $H$ (which is the same as the singular vector corresponding to the largest singular value of $G$). This vector is the time-reversed scattered field corresponding to the strongest scatterer, restricted to the measurement array. When this field is transmitted, it focuses on the strongest scatterer and this property can be used to selectively transmit energy or information to the scatterer. To locate the scatterer one only needs to search for the maximum of this transmitted field. When the time signal is time-harmonic with frequency $\omega $, the operation of time-reversal of a signal reduces to taking the conjugate of its complex representation and time-reversal operator is related to scattering operator via $H=GG*$, where $G*$ is the complex conjugate of $G$.

The time-reversal iterations described above are effective when there is one dominant scatterer in the medium. Their usefulness is limited when multiple scatterers of comparable strength are present. In the case of multiple scatterers the iterations converge to the strongest scatterer at a rate that is proportional to the ratio of the strength of the second-strongest scatterer to the strongest scatterer. Clearly, in the case of multiple scatterers of similar strength, this convergence can be poor. Further, the recovery of weaker scatterers is cumbersome and requires a large number of iterations, even when their strengths are very different.^{8} The slow convergence of iterated time-reversal can be problematic in applications.^{9} In a temporally (slowly) changing environment, iterated time-reversal is particularly useful to adaptively optimize focusing to the current conditions. Such adaptation is impossible, however, if the convergence time is not significantly shorter than the medium correlation time. Therefore, accelerating the convergence of iterated time-reversal is key to realizing its full potential in applications.

The key to accelerating this process is the recognition that the transducer array can be used as a kind of analog computer. Each time-reversal operation can be interpreted as a matrix-vector product. Thus advanced iterative methods from linear algebra can be implemented in a time-reversal context.

When compared with power iterations the convergence properties of methods based on Lanczos iterations are superior, especially when multiple eigenvalues and eigenvectors are desired.^{10} As we explain below the estimate from Lanczos is always at least as good as that from the power method and is typically far superior. In addition, for multiple eigenvalues the power method requires a complicated multi-step procedure.^{8}

*Motivated by these observations we describe a new iterated time-reversal method based on Lanczos iterations that inherits their advantages over power iterations*. Through numerical examples we demonstrate here that when compared with the power iterations based method it requires fewer iterations to converge to an accurate answer. This observation directly translates to fewer transmit and receive cycles in an experimental setting. We also note that the experimental realization of the proposed method requires only a slight modification of the current protocol for time-reversal and no modification to the actual hardware. It is worth reemphasizing that the method proposed here and described below offers rapid convergence to multiple eigenvectors with few transmission cycles. The method *does not* depend upon measuring the entire scattering operator.

## II. Iterated time-reversal

Time-reversal concepts are conveniently described in terms of scattering and time-reversal operators. Consider a surface in $R3$ denoted by $\Gamma $ that is continuously embedded with transmitters and receivers (discrete transmitters and receivers may be treated identically). In the following treatment, we consider the standard development in time-reversal using single-frequency harmonic signals. Thus the transmitted and received signals are functions of $x\u220a\Gamma $ and frequency $\omega $. The dependence of signals and operators on the frequency is implicit, so $\omega $ is dropped in the notation. We require that signals be square-integrable, that is they are functions $v\u220aV\u2261L2(\Omega )$. We denote the $L2$ inner product on $\Gamma $ by $(u,v)=\u222b\Gamma u*vdx$ and the corresponding norm by ‖⋅‖. The kernel of the scattering operator of the problem is denoted by $g(x;x\u2032)$, which is the Green’s function corresponding to the appropriate wave propagation problem. It represents the response at a location $x$ due to a unit impulse at $x\u2032$. The received signal, corresponding to a transmitted signal $v(x)$, is given by

The relation above may be written succinctly as $r(x)=G[v(x)]$, by defining the scattering operator $G:V\u2192V$. The kernel of the time-reversal operator is

where the asterisk denotes the operation of conjugation. The time-reversed signal corresponding to a transmitted signal $v$ is given by Eq. (1) with $g$ replaced by $h$. This can be written succinctly as $r(x)=H[v(x)]$, by defining the time-reversal operator $H:V\u2192V$.

Let ${s(i),\varphi (i)(x)}$ be the $ith$ singular value and (right) singular vector of the scattering operator, $G$, arranged and ordered so that $(\varphi (i),\varphi (j))=\delta ij$ and $\u2223s(1)\u2223\u2265\u2223s(2)\u2223\u2265\cdots $. Then from Eq. (2) and the time-reversibility of the wave equation, we conclude that the eigenpairs of $H$ are ${\lambda (i),\varphi (i)(x)}$, where $\lambda (i)=\u2223s(i)\u22232$. We note that the singular vectors of $G$ are also eigenvectors of $H$. Under certain conditions it can be shown that if an eigenvector of the time-reversal operator, $H$, is used as the transmitted signal, the wave field selectively focuses on the scatterer associated with the corresponding eigenvalue. Thus determining the eigenvectors of $H$, or, equivalently, the singular vectors of $G$, allows targeting specific scatterers in the propagation domain. Next we describe how this is accomplished by the traditional iterative time-reversal method via power iterations. Thereafter we describe the use of Lanczos iterations.

**Power iterations:** We define the following sequence as a single time-reversal iteration: given $v(0)(x)$; transmit $v(0)*(x)$; receive $v(1\u22152)(x)=G[v(0)*(x)]$; transmit $v(1\u22152)*(x)$: and receive $v(1)(x)=G[v(1\u22152)*(x)]$. Using the definition of the scattering and the time-reversal operators, this sequence can be written as $v(1)(x)=H[v(0)(x)]$, where $H=GG*$. In the time-reversal method based on power iterations the sequence described above is repeated $n$ times to yield $v(n)=Hn[v(0)]$. Using the spectral decomposition of $H$ it is easy to show that as $n$ increases $v(n)$ converges to the eigenvector corresponding to the largest eigenvalue of $H$. That is $v(n)\u2215\Vert v(n)\Vert \u2248\varphi (1)$, for $n$ large.

The eigenvector corresponding to the second-largest eigenvalue can be determined once the first eigenvector has been determined by performing power iterations with $v(0)\u2212(v(0),\varphi (1))\varphi (1)$, where $\varphi (1)$ is the estimate of the first eigenvector. That is, for every iterate, the component in the direction of $\varphi (1)$ is annihilated. This or a similar process can be repeated to locate weaker scatterers.^{8}

**Lanczos iterations:** Lanczos iterations may be used to determine the eigenvectors of the time-reversal operator.^{10}

Select arbitrary $v(1)$ such that $(v(1),v(1))=1$.

Set $v(0)=0$ and $\beta 1=0$.

For $j=1,\u2026,n$

$w(j)=H[v(j)]\u2212\beta jv(j\u22121)$

$\alpha j=(w(j),v(j))$

$w(j)=w(j)\u2212\alpha jv(j)$

$\beta j+1=(w(j),w(j))$

$v(j+1)=w(j)\u2215\beta j+1$

After $n$ iterations we will have created the $n\xd7n$ tridiagonal matrix, $T=tridiag(\beta i+1,\alpha i,\beta i+1)$, $i=1,\u2026,n$. Let ${\gamma (i),\psi (i)}$ be the $ith$ eigenpair for $T$. *Then the eigenvalues of* $T$ *approximate the eigenvalues of* $H$, *that is* $\lambda (i)\u2248\gamma (i)$ *and the eigenvectors of* $H$ *can be determined from the eigenvectors of* $T$ *and the Lanczos vectors* $v(i)$, *that is* $\varphi (i)\u2248\u2211j=1n\psi j(i)v(j)$.

There are several comments to be made here. First, in implementing this algorithm in an experiment the Lanczos vectors $v(i)(x)$ would have to be stored. Though this is an added memory cost over power iterations, it is not a big penalty as they may be stored on the hard disk of a computer. Second, like the power iterations this algorithm also involves a time-reversal step (step 3a). Thus its implementation in an experiment is similar to that of power iterations, except that every time-reversal step is followed by some additional signal processing. Finally, at the completion of $n$ iterations, approximations to $n$ eigenpairs are available. Out of these a small fraction $(\u224820%)$ can be expected to be accurate. The accuracy may be increased by increasing the number of iterations

**Theoretical advantages of Lanczos:** The Lanczos method finds the best approximation of the $n$ largest and smallest eigenvectors with an $n+1$ dimensional Krylov space spanned by the vectors ${v(0),Hv(0),H2v(0),\u2026,Hnv(0)}$. The power method’s estimate for the largest eigenvector is $Hnv(0)$, which is contained within the Krylov space. Clearly then, *for the same number of iterations, the Lanczos iterated time-reversal method always provides a more accurate estimate of the eigenvector than does standard iterated-time reversal*. Furthermore, from those same iterations, we have very accurate estimates of eigenvectors for smaller eigenvalues; no further iterations are required.

## III. Numerical examples

**A simplified problem:** We consider a three-dimensional domain with $N$ point scatterers located at $xj$, $j=1,\u2026,N$, clustered about the origin within a sphere of radius $\delta $. The transmitters∕receivers are continuously embedded in the $y\u2212z$ plane at $x=D$ (denoted by $\Gamma $), with $D\u2aa2\delta $. We consider the time-harmonic case with frequency $\omega $ and wavenumber $k=\omega \u2215c$, where $c$ is the sound speed in the medium. Further, we neglect multiple scattering and assume that scatterers are ideally separated.^{11} That is, it is possible to focus on one scatterer without directing any energy to another scatterer. For this case ${s(j),\varphi (j)(y,z)}$ is the $jth$ eigenpair of the scattering operator, where $s(j)$ is proportional to the strength of the $jth$ scatterer and $\varphi (j)(y,z)=g(j)(D,y,z)\u2215\Vert g(j)(D,y,z)\Vert $. In this expression ‖⋅‖ denotes the $L2$ norm on $\Gamma $ and $g(j)(x)=eik\u2223x\u2212xj\u2223\u22154\pi \u2223x\u2212xj\u2223$. The time-reversal operator is given by $H=GG*$, where ∗ denotes the complex conjugate. The eigenpairs of $H$ are given by ${\u2223s(j)\u22232,\varphi (j)(y,z)}$. Our goal is to use time-reversal methods to determine the $\varphi (j)$.

First we consider two scatterers $(N=2)$. We start with an initial arbitrary signal $v(0)(y,z)$ and let $\rho j=(\varphi (j),v(0))$ be the components of this signal along $\varphi (j)$. Then after $n$ steps of time-reversal iterations using the power method, the signal $v(n)$, which approximates $\varphi (1)$, is given by

From Eq. (3) it is clear that the error is $(\rho 2\u2215\rho 1)\u2223s(2)\u2215s(1)\u22232n$. Assuming that the initial guess is arbitrary and thus not aligned with either $\varphi (1)$ or $\varphi (2)$ we expect $\rho 2\u2215\rho 1=O(1)$ and conclude that the error is determined by the ratio of the strength of the scatterers. For scatterers of similar strength, say $s(2)=0.98s(1)$, after ten time-reversal iterations the error is $\u22480.9820=0.66$. For scatterers with disparate strengths, say $s(2)=0.70s(1)$, this error is much smaller $(\u22480.08%)$. Because the Lanczos method recovers in $m$ iterations the eigenpairs of an operator of rank $m$,^{10} we recover $\varphi (1)$ and $\varphi (2)$ in exactly two iterations.

Next we consider scatterers with $s(j)$ selected from a normal distribution of mean 0.5 and standard deviation 1. We use time-reversal methods to determine the eigenfunctions corresponding to the two strongest scatterers. For power iterations we use the first ten iterations to determine a guess for $\varphi (1)$, then use another ten iterations to determine a guess for $\varphi (2)$. This corresponds to a total of 20 time-reversal steps. In the Lanczos iterations based method also we use a total of 20 iterations. In order to obtain reliable indicators of performance we use 10 000 realizations. For each realization we determine the normalized error in estimating the first two eigenfunctions and plot a histogram of the log of this error (see Fig. 1).

Figure 1(a) shows histograms of the error in estimating the first eigenfunction by both the power method and the Lanczos method. We observe that the typical error in the Lanczos estimate is roughly six orders of magnitude smaller than that for the power method. A similar trend is observed in the histograms shown in Fig. 1(b) for the second eigenfunction error, though there the power method estimate is frequently unusable.

**Two-dimensional problems:** Here we consider multiple scattering in two dimensions with 40 transmitters and receivers that produce waves with a wavelength of $0.25units$ and are equally spaced in the interval (0, 6) along the $x$-axis. The medium contains $28(weak)$ point scatterers located within a square of edge length 6 that is centered at $(x,y)=(4.5,4.5)$. The strength of these scatterers is randomly distributed in the interval ]0, 0.2[. There are two other scatterers that are significantly stronger (as described below). Eight time-reversal iterations based on power and Lanczos iterations are used to locate the position of the two strongest scatterers. For this purpose, the computed value of the eigenvector of the time-reversal operator is back-propagated from the transmitters, once the iterations are complete. It is expected that the back-propagated field will focus on the scatterer we wish to target.

In the first case the strength of strongest scatterer is 2.1 and it is placed at $(x,y)=(4,3)$, while the strength of the second-strongest scatterer is 2 and it is placed at $(x,y)=(2,3)$. The back-propagated fields are plotted in Fig. 2. We observe that for both eigenvalues the Lanczos iterations focus tightly on the corresponding scatterer, while the power iterations based method is unable to distinguish between the two scatterers.

We next place the strongest scatterer further from the array [$strength=3$; located at $(x,y)=(2,4.5)$], and the weaker scatterer closer to the array [$strength=2$; located at $(x,y)=(4,3)$], so that the weaker scatterer produces a stronger scattered field. In Fig. 3 we plot the back-propagated field for the two largest eigenvalues. We observe that the largest eigenvalue corresponds to the second-strongest scatterer and that in this case also Lanczos iterations are better at selectively focusing on scatterers.

## IV. Conclusions

We have presented a new iterated time-reversal method for efficiently identifying the eigenvectors of a scattering operator, which can be used (e.g., by DORT,^{12} MUSIC^{13}) to focus on or locate strong scatterers in a propagating medium. Our method makes use of Lanczos iterations in order to converge to the eigenpairs of the time-reversal operator. In numerical tests this method is shown to be remarkably superior over traditional iterated time-reversal. The new approach requires relatively little additional hardware and processing; it requires the storage of several (as many as time-reversal iterations) time signals for each transmitter on a hard drive, and some extra numerical processing. Its rapid convergence makes it well suited to applications.