The angular spectrum approach (ASA)—a frequency domain method to calculate the acoustic field—enables highly efficient passive source localization and modeling forward propagation in homogeneous media. If the medium is continuously stratified, a first-order analytical solution may be obtained for the field at arbitrary depth. Simulations show that the proposed stratified ASA solution enables accurate source localization as compared to the uncorrected ASA (error from 1.2 ± 0.3 to 0.49 ± 0.3 wavelengths) at scalings relevant to biomedical, underwater, and atmospheric acoustic applications, and requiring milliseconds on nonspecialized hardware. The results suggest the proposed correction enables efficient and accurate localization in stratified environments.

## 1. Introduction

The angular spectrum approach (ASA) is a method of solving the wave equation in the spatial frequency domain.^{1} It may be considered a decomposition of the harmonic field into plane waves that travel with a continuous spectrum of directions, with appropriate phases and amplitudes such that their summation yields fields of arbitrary spatial distribution.^{2} Plane waves are mathematically convenient, as the propagation through space of each component may be modeled with a simple phase delay. Coupled with fast Fourier transform (FFT) algorithms, the method enables exceptionally efficient computations and is naturally suited to real-time applications with relatively narrow-band sources. Of particular interest for the ASA here is the recovery of the volumetric wave field from the measured pressure on a planar surface (i.e., boundary condition), a process known as acoustical holography.^{3} From the ASA reconstruction, a passive acoustic map (PAM) may be formed at specific frequencies, and thus sources at these frequencies to be localized as intensity peaks in the PAM.

In practical applications, this surface measurement is made by a series of sensors (e.g., an ultrasonic imaging transducer or linear array of hydrophones) that discretely sample the pressure field due to sources in the region of interest. Such source localization from PAMs is a problem of significant interest in biomedical,^{4–6} underwater,^{7–9} and aeroacoustic^{10,11} applications. However, in such applications, the propagation environment typically varies as a function of the vertical direction. While recent work^{12} has shown that heterogeneity can be accounted for with the ASA, an analytical solution might ensure high computational efficiency for large imaging domains to reduce errors in the source localization.

Herein, we derive a first order analytical correction for the ASA and the case of a stratified medium. Then, through simulated acoustic propagation, we demonstrate corrected sub-wavelength localization in stratified environments at scales relevant to biological, atmospheric, and underwater acoustics. Finally, we consider the point source localization errors, and the computational cost of the different algorithm permutations.

## 2. Methods

### 2.1 Angular spectrum approach for stratified media

The angular spectrum $P(kx,ky,z)$ of a monochromatic ($\u221d\u2212i\omega t$) field $p\u0303(x,y,z)$ defined by

Taking the two-dimensional (2D) Fourier transform to the homogeneous Helmholtz equation $(\u22072+\omega 2/c02)p\u0303=0$ with a constant speed of sound *c*_{0} gives an ordinary differential equation for *P*,

where $kz2=(\omega /c0)2\u2212kx2\u2212ky2$. With knowledge of the boundary condition *P*_{0} at *z* = 0 and the assumption that there are no backward-travelling waves, then Eq. (2) has the solution

The acoustic field in any plane may then be reconstructed by evaluation of the inverse transform of *P* given by Eq. (3).

Suppose now that the sound speed has weak spatial dependence (i.e., that the sound speed changes over scales that are large compared with the wavelength), such that $c0\u2192c(r)$ in the Helmholtz equation. Then, with the definitions $\mu =c02/c2(r)$ and $\lambda =(1\u2212\mu )\omega 2/c02$ for a constant reference sound speed *c*_{0}, it can be shown that the governing equation is

Here, $\Lambda =Fk[\u2009\lambda \u2009]$, and $*$ represents the 2D convolution over the wavenumber components *k _{x}* and

*k*. If

_{y}*c*varies only with the axial coordinate

*z*, the convolution may be evaluated to give

Assume a Wentzel-Kramers-Brillouin (WKB)-type solution^{13,14} of the form

where *A* is a complex amplitude. Substitution of Eq. (6) into Eq. (5) and evaluation of the derivatives yields

To first order, the first term in Eq. (7) can be neglected to obtain a first-degree ordinary differential equation for *A*, which may be integrated directly to give

Neglecting the higher order terms requires that changes in sound speed occur on scales that are very large compared to the wavelength ($|(dc/dz)\Delta z/c|$), and that the relative magnitude of these changes is not too large ($|c/c0|\u223c1$, see the supplementary material^{23}). Application of the boundary condition (i.e., the Fourier transform of the measured field at *z* = 0) gives

The angular spectrum at arbitrary *z* is then

### 2.2 Simulations and numerical implementation

To determine the improvement in source localization, point sources were simulated in a half-space above a virtual receiving array in k-Wave;^{15} see Fig. 1(a). For computational efficiency, simulations were performed in 2D, and thus in the simulations and reconstructions, $y=ky=0$. The PAM reconstruction routines were written in matlab and computation times are reported for a standard desktop computer (Intel Core i7, four cores at 2.8 GHz and 16 GB memory); no parallel or graphical processing techniques were employed. Sources were simulated as Gaussian pulses with 5% bandwidth. The resulting pressure was measured by a virtual linear array; the acoustic data were then beamformed to compute $P(kx,z)$ with Eq. (3) or Eq. (10) for the uncorrected and corrected cases, respectively. Finally, the intensity field was then calculated from

Maps were formed over a frequency range *ω* of three frequency bins centered at the source frequency. The reconstructed source location $(xr,zr)$ was taken to be the position of peak intensity of *I*, and the error was defined relative to the known source position $(xtrue,ztrue)$ as

To prevent spatial aliasing, all initial spectra were zero padded such that their computational spatial extent was four times larger than their physical extent. As the medium is stratified, the material properties are known for all *x*, and thus the computational domain may be extended arbitrarily from the simulation domain. The measured spectra $p\u03030$ were windowed with a Tukey window with cosine fraction *R* = 0.25 to taper the high spatial frequency components that would be introduced by the discontinuity caused by the zero-padding.^{3} All reconstructions were performed with data from the virtual sensors spaced at $\Delta x$ and at depth steps $\Delta z=f0/6c0$ (i.e., 1/6 of the reference wavelength); see Fig. 1(g). The discrete integration in Eq. (10) was computed as a Riemann sum with the same $\Delta z$. Finally, the time step was chosen to satisfy the Courant–Friedrichs–Lewy condition for the frequencies of interest, and ensure stability of the simulations.^{15} In all cases, a perfectly matched layer of at least two wavelengths was used at all boundaries (with the number of points chosen to improve the simulation efficiency). As the ASA is relatively insensitive to noise,^{6} reconstructions did not generally include added noise. To confirm that this effect persisted with the proposed correction, select reconstructions were performed with added white noise with amplitudes up to five times the peak signal level.

In the biomedical-scale applications, a mean sound speed of 1540m/s was augmented by a 25% Gaussian profile with variance of 30 mm. A mean density of 1043 kg/m^{3} and attenuation of 0.54 dB/cm/MHz were defined for the entire medium^{16} and 99 narrowband sources of 1 MHz were distributed uniformly in the medium as in Fig. 2(a); the simulation grid had resolution $\Delta x=\Delta z=0.2\u2009mm$, and $\Delta t=50\u2009ns$. For underwater simulations, 70 narrowband sources of 1 kHz were distributed uniformly underwater as shown in Fig. 3(a), and the length of the array located on the water surface was 190 m. The sound speed profile was taken from Ref. 17 at 0.5°N and 100.5°W, with attenuation defined from an empirical frequency-dependent model;^{18} the simulation grid for the underwater simulations had resolution $\Delta x=\Delta z=60\u2009cm$, and $\Delta t=5\u2009\mu s$. Last, for the atmospheric simulations, 209 narrowband sources of 1 Hz were distributed uniformly above ground as shown in Fig. 4(a), and the length of the array located on ground was 6 km. The medium density and sound speeds were defined from the 1972 US Standard Atmosphere,^{19} and the altitude- and frequency-dependent attenuation was taken from Ref. 20. In the atmospheric simulations, the simulation grid had resolution $\Delta x=\Delta z=10\u2009m$, and $\Delta t=20\u2009ms$.

## 3. Results

Figure 2 demonstrates the improvement in source localization due to the phase correction for 1 MHz sources. Without the phase correction, the error was 2.05 ± 1.0 mm, whereas with the phase correction it was 0.97 ± 0.2 mm. The error in both the corrected and uncorrected case was principally in the vertical (axial) direction ($|\u03f5z,avg/\u03f5x,avg|=34.6$ and 35.6 for the corrected and uncorrected cases, respectively) and is plotted in Fig. 2(b). To understand the effect of the beamforming aperture, we beamformed and localized using three different aperture sizes. For the 50 mm aperture, some localizations in the corrected case have larger errors; however, for larger apertures (75 and 100 mm) that cover the horizontal extent of the sources, the corrected localization accuracy was within approximately one half wavelength [Fig. 2(c), orange]. In the uncorrected case, (purple) the error was largest for depths near 25 and 75 mm. These are the extrema of $dc/dz$, and thus where discounting the medium variation is most egregious.

Finally, the absolute localization accuracy did not depend strongly on the wavelength [Fig. 2(d)]. The absolute localization error was similar for all frequencies (mean 0.91, 0.97, and 0.88 mm for 0.5, 1.0, and 1.5 MHz, respectively), as was the improvement relative to the uncorrected case at that frequency (mean improvement 55%, 57%, and 53%). The short wavelength criterion [Eq. (S-1) in the supplementary material^{23}] was roughly met in all cases; for the longest wavelength (500 kHz), $|d2A/dz2|\u223c0.4$.

For the simulations with data from measured underwater sound speed profile, the magnitude of the $|c\u2212c0|/c0$ was small—approximately 2% [Fig. 3(a)]. Accordingly, the source localization even in the uncorrected case was reasonable, with errors on the order of a wavelength (1.38 ± 0.12 m compared with the wavelength of 1.5 m) over all depths [Fig. 3(b)]. The correction successfully reduced the error by approximately fivefold, to 0.26 ± 0.08 m. The axial component comprised most of the total error in the uncorrected case ($|\u03f5z,avg/\u03f5x,avg|=9.2$); however, in the corrected case, the axial and transverse errors were comparable ($|\u03f5z,avg/\u03f5x,avg|=1.5$).

In the atmospheric case, the sound speed decreased approximately linearly over the altitudes considered [up to 9 km, Fig. 4(a)]. Unlike the other cases considered here, the simulation medium also had significant density variations [$|\rho \u2212\rho 0|/\rho 0\u223c0.5$, Fig. 4(b)]. Over the range of positions considered, the mean source localization error was reduced from 474.4 ± 395.6 m in the uncorrected case to 219.4 ± 226.3 m with the correction, compared with a wavelength of 343 m [Fig. 4(c)].

### 3.1 Effect of noise

To demonstrate the robustness of both the homogeneous and stratified ASA solutions to noise, the data from the simulations whose results are reported in Fig. 2 (100 mm aperture, 1 MHz sources) were beamformed after the addition of uniform white noise with amplitude from one to five times that of the maximum value of the recorded RF data. The resulting PAMs show clear peaks with comparable error to the noiseless case [Fig. 5(a)] despite apparent total corruption of the time series data [Fig. 5(b)]. Averaged over the entire grid [Fig. 2(a)], the error without noise was 0.97 ± 0.20 and 2.0 ± 0.97 mm in the corrected and uncorrected cases, respectively. At the noise level of 5, the average errors were 0.93 ± 0.25 mm with the correction and 2.0 ± 1.0 mm uncorrected. The noise manifests in the PAMs as irregular interference patterns [arrow in Fig. 5(a)], which can cause aberrant high intensities and thus spurious localizations at still larger noise values. However, Fig. 5 confirms that both the conventional and stratified medium ASA techniques are suitable for noisy conditions.

### 3.2 Computational efficiency

The time required for evaluation of Eq. (11) depends on the number of frequency bins, signal duration, and the size of the computational grid. For the cases described here, reconstruction times were on the order of tens of microseconds for the uncorrected case (three frequency bins) and on the order of 100 ms when the correction was used. In all cases, this reduced from tens to a few hundred nanoseconds per pixel in the reconstructed image (see Table 1).

Environment . | Uncorrected (ms) . | Corrected (ms) . | Grid size ($Nx\xd7Nz\xd7Nt$) . |
---|---|---|---|

Biomedical | 43.4 ± 43.5 | 124.7 ± 125.0 | $588\xd7600\xd72000$ |

Underwater | 150.5 ± 2.2 | 437.5 ± 12.9 | $314\xd7320\xd77000$ |

Atmospheric | 40.4 ± 1.3 | 112.3 ± 5.0 | $588\xd7600\xd72500$ |

Environment . | Uncorrected (ms) . | Corrected (ms) . | Grid size ($Nx\xd7Nz\xd7Nt$) . |
---|---|---|---|

Biomedical | 43.4 ± 43.5 | 124.7 ± 125.0 | $588\xd7600\xd72000$ |

Underwater | 150.5 ± 2.2 | 437.5 ± 12.9 | $314\xd7320\xd77000$ |

Atmospheric | 40.4 ± 1.3 | 112.3 ± 5.0 | $588\xd7600\xd72500$ |

## 4. Discussion

In this paper, we have presented an analytical phase correction for ASA beamforming in a stratified medium. Under the assumption that the sound speed changes slowly compared to a wavelength and that the magnitude of the change is small, an analytical solution the governing equation accurate to first order may be obtained. Beamforming of passively collected data with the correction mitigates aberration caused by the constant sound speed assumption of the canonical ASA (Fig. 1). Through simulations of point source radiation at biomedical (Fig. 2), underwater (Fig. 3), and atmospheric (Fig. 4) scalings, the error was reduced by 62.6%, averaged across all cases. Collectively, the results indicate that the method has value for passive source localization at several realistic environments where such information is of considerable interest.^{21}

While the computational cost was approximately three-fold higher in computational time (mean 147 ms), this is still orders of magnitude more efficient than time-domain methods that can handle such corrections natively. The computation time could be reduced further by reducing the number of receivers (provided the spacing between remains smaller than a half wavelength i.e., the spatial Nyquist frequency),^{3} reducing the signal duration or sampling rate *f _{s}* (provided that the full source signal is captured and $fs\u22652fmax$, i.e., the Nyquist limit), or by use of parallel computations to reconstruct the maps at each frequency. Thus the localization method extends the real-time capabilities of the ASA.

The method has a few limitations. First, the derivations require that the relative changes in sound speed occur over long scales (i.e., $|(dc/dz)\Delta z/c|$ is very small over a wavelength), and secondarily that these changes have relatively small magnitudes (i.e., $|c/c0|\u223c1$). Thus, discrete impedance changes (e.g., layers) violate these assumptions (however, such cases may be treated with a layer-by-layer approach).^{22} Second, only forward propagation is included, such that reflections (e.g., from the ground or water surface) are not accounted for by the method.

## 5. Conclusion

Here, we derive an analytical, first-order correction to the angular spectrum approach for media with stratified material properties. This method extends the ability of the inherently efficient method to localize individual sources at subwavelength accuracy with little additional computational burden.

## Acknowledgments

This work was supported by NIH Grant No. R00EB016971 (NIBIB) and NSF Grant No. 1830577 (LEAP HI).