Traditional matched-field processing (MFP) refers to array processing algorithms, which fully exploit the physics of wave propagation to localize underwater acoustic sources. As a generalization of plane wave beamforming, the “steering vectors,” or replicas, are solutions of the wave equation descriptive of the ocean environment. Thus, model-based MFP is inherently sensitive to environmental mismatch, motivating the development of robust methods. One such method is the array invariant (AI), which instead exploits the dispersion characteristics of broadband signals in acoustic waveguides, summarized by a single parameter known as the waveguide invariant *β*. AI employs conventional plane wave beamforming and utilizes coherent multipath arrivals (eigenrays) separated into beam angle and travel time for source-range estimation. Although originating from the ideal waveguide, it is applicable to many realistic shallow-water environments wherein the dispersion characteristics are similar to those in ideal waveguides. First introduced in 2006 and denoted by *χ*, the dispersion-based AI has been fully integrated with *β*. The remarkable performance and robustness of AI were demonstrated using various experimental data collected in shallow water, including sources of opportunity. Further, it was extended successfully to a range-dependent coastal environment with a sloping bottom, using an iterative approach and a small-aperture array. This paper provides an overview of AI, covering its basic physics and connection with *β*, comparison between MFP and AI, self-calibration of the array tilt, and recent developments such as adaptive AI, which can handle the dependence of *β* on the propagation angle, including steep-angle arrivals.

## I. INTRODUCTION

Over the last five decades matched-field processing (MFP) has been extensively studied in the literature for localization of underwater acoustic sources.^{1–6} MFP refers to array processing algorithms, which exploit the full field structure of signals propagating in an ocean waveguide. It is a generalization of conventional plane wave beamforming in which the “steering vectors,” or replicas, are not plane waves but solutions of the wave equation that is descriptive of the ocean environment (i.e., Green's function). The generalized matched-field beamformer matches the measured field at the array (typically vertical) with replicas of the expected field for all possible source locations in range and depth, referred to as the ambiguity surface shown in Fig. 1(a), where the main peak (best match) corresponds to the source location. Exploiting the unique spatial structure of the field (complexity) requires accurate knowledge of the environment, which includes the sound speed in the water column and layering within the sediment even in the most simple cases. As a result, model-based MFP is inherently sensitive to the environmental mismatch that is inevitable in real-world situations, motivating the development of robust methods such as focalization^{7} and multiple-constraint^{8,9} and frequency-difference^{10} MFP. For a good overview and history of MFP, we refer the reader to Refs. 1 and 2, and the references therein.

As a paradigm shift in source ranging, the so-called array invariant (AI) emerged recently; it requires no modeling of the acoustic field and little to no environmental information in shallow water.^{11–15} Denoted by *χ*, AI exploits the dispersion characteristics of broadband signals that can be summarized in a single scalar parameter known as the waveguide invariant,^{16–18} denoted by *β*. In contrast to model-based MFP, the dispersion-based AI employs conventional plane wave beamforming and utilizes coherent multipath arrivals (eigenrays) separated into beam angle (*θ*) and travel time (*t*), referred to as the “beam-time migration” wherein the source-range information is embedded [see Fig. 1(b)]. On the other hand, the well-known waveguide invariant is defined by the functional relationship between phase and group speed (or slowness) and provides a constant value for a group of normal modes that behave similarly in terms of dispersion.^{19,20} Since beam-time migration corresponds to the dispersion curve in phase versus group slowness, a simple formula can be derived for the source range (*r*_{0}) (refer to Sec. II),

where *c* is the nominal sound speed of 1500 m/s. The AI parameter (*χ*) defined in Eq. (3) is estimated directly from the beam-time migration data shown in Fig. 1(b), thus requiring knowledge of *β* for source-range estimation.

The original version of AI developed in 2006 from the ideal waveguide^{11} missed the connection between the two invariants (*χ* and *β*), but remains valid since *β* = 1 for small angles (e.g., $\theta <20\xb0$) in ideal waveguides, whose analytical expression is $\beta =cos2\theta $. It turns out that Søstrand^{21} actually applied the waveguide invariant (*β*) and derived Eq. (1) about a year earlier to estimate the range of small explosive charges using a horizontal array, without explicitly defining the array invariant (*χ*). To our delight, the dispersion-based AI is seamlessly applicable to many realistic shallow-water environments with a non-uniform sound speed profile (SSP) and seabed for surface/bottom-reflected arrivals, because their dispersion characteristics are similar to those in ideal waveguides. Moreover, a closed-form expression *β* derived for a Pekeris waveguide with a penetrable bottom^{22} that is more descriptive of shallow-water environments also can be approximated by $\beta =cos2\theta $, further validating the use of an ideal waveguide invariant in shallow-water environments.

In theory, Eq. (1) implies that AI can be extended to a general refracting environment with a different value of *β* (e.g., $\beta =\u22123$ for a uniform sound speed gradient).^{12} However, *β* can be unstable in the presence of small variations in the profile.^{23} In deep water with a refractive SSP and no bottom interaction, *β* varies widely between negative and positive values and also depends on the source range and depth,^{16,24–26} making it impractical for source-range estimation. Thus, the content of this paper is limited to shallow waters endowed with a small variation of *β* around unity for practical applications, where *β* is approximately independent of frequency, range, or source and receiver depths in a range-independent environment.^{16,18}

The objective of this paper is to provide an overview of AI for readers who would be interested in applying this simple and robust approach to source localization in shallow water. Section II presents the definition of AI and derivation of Eq. (1) from two different perspectives, one from the dispersion in the ideal waveguide and the other from the definition of the waveguide invariant. The connection between the two invariants (*χ* and *β*) is then discussed in detail, including the basic physics of the waveguide invariant. Section III investigates the beam-time migration for horizontal and vertical arrays and provides the least squares (LS) estimate of *χ* for source-range estimation using Eq. (1). In Sec. IV, the two invariants (*χ* and *β*) are combined into a single invariant referred to as the adaptive AI ($\chi \beta =\chi /\beta $), which can handle the dependence of *β* on the propagation angle or steep-angle arrivals. Since obtaining the beam-time migration requires knowledge of the time-domain Green's function, Sec. V reviews the ray-based blind deconvolution^{27} for estimation of the Green's function from an unknown broadband source such as a ship. Section VI compares the performance of MFP and AI, using a ship of opportunity and a large-aperture vertical array in a relatively range-independent shallow-water environment. While both methods are sensitive to array tilt, it is shown that AI is capable of self-calibrating the array tilt. Last, Sec. VII extends the AI to a mildly (adiabatic) range-dependent coastal environment with variable bathymetry, using an iterative approach and a short-aperture array, followed by a summary in Sec. VIII.

## II. FORMULATION OF ARRAY INVARIANT (AI)

This section derives Eq. (1) from two different perspectives through the normal modes of a waveguide. One is based on the dispersion property in ideal waveguides, and the other is based on the definition of waveguide invariant valid for general range-independent environments.

### A. Dispersion in ideal waveguides, *χ*

Consider an isovelocity ideal waveguide with a pressure-release surface and bottom (see Fig. 2). The basic ingredient is the modal group speed represented by $vg=c\u2009\u2009cos\u2009\theta $, where *θ* is the modal propagation angle with respect to the horizontal (i.e., grazing angle). The modal index has been suppressed for simplicity. For a source at range *r*_{0}, we substitute $vg=r0/t$ with *t* denoting the travel time and obtain

Note that the cosine of the grazing angle (rather than the angle itself) is inversely proportional to the travel time (*t*).

The array invariant *χ*, which involves beamforming using an array, is defined as the derivative of the cosine of the beam angle with respect to the travel time, equivalent to the slope of the beam-time migration curve in $(t,\u2009cos\u2009\theta )$ coordinates,

where $t=r0/(c\u2009cos\u2009\theta )$ is substituted from Eq. (2). The term $\u2009cos2\theta $ on the right-hand side of Eq. (3) is recognized as the analytic expression for *β* in ideal waveguides (see Sec. II C),^{19,23} i.e., $\beta =cos2\theta $. The substitution completes a concise derivation of Eq. (1).

### B. Waveguide invariant, *β*

In this section, Eq. (1) is derived from the waveguide invariant assuming a range-independent waveguide (not limited to ideal waveguides), which is extended to a mildly (adiabatic) range-dependent environment^{28,29} in Sec. VII. The waveguide invariant *β* is defined as^{18,19}

where *v _{g}* and

*v*are the group and phase speed, and their inverses are the group and phase

_{p}*slowness*,

*s*and

_{g}*s*, respectively. Substituting $sg=t/r0$ and $sp=cos\u2009\theta /c$ with

_{p}*c*being the local sound speed at the receiver depth, Eq. (4) can be written as

The term in brackets on the right-hand side of Eq. (5) is recognized as the inverse of *χ* defined in Eq. (3), making Eq. (5) identical to Eq. (1). While *c* normally depends on the sound speed at the local receiver depth, even a deviation of $\Delta c=15$ m/s from the nominal sound speed of 1500 m/s leads to a range error of just 1%, and so the impact of using *c* = 1500 m/s in Eq. (1) is minuscule for practical purposes.

### C. Connection between the two invariants, *χ* and *β*

The above derivation clearly reveals the link between the two invariants, *χ* and *β*, as follows. The beam-time migration in $(t,\u2009cos\u2009\theta )$ coordinates (see Fig. 3) corresponds to the dispersion curve in (*s _{g}*,

*s*) coordinates, and the slope of the beam-time migration (

_{p}*χ*) is proportional to the slope of the dispersion curve (

*β*) with an opposite sign and the proportional coefficient containing the source-range information ($\u2212c/r0$). In other words,

*χ*and

*β*are two sides of the same coin where the beam-time migration (data) simply reflects the theoretical dispersion in the waveguide for a specific source range (

*r*

_{0}). Thus, it is imperative to understand the physics of the waveguide invariant

^{16–19,23}for the proper use of AI.

In ideal waveguides, all rays or modes interact with both the surface and the bottom, and we have the relationship of $vgvp=c2$ with $vg=c\u2009cos\u2009\theta $ and $vp=c/\u2009cos\u2009\theta $. Then Eq. (4) provides the analytical expression of $\beta =cos2\theta $, which is a function of the grazing angle *θ* (so is *χ*). The Pekeris waveguide with a penetrable bottom^{22} is also shown to have a closed-form expression that is approximated by $\beta =cos2\theta $. The far-reaching benefit is that $\beta =cos2\theta $ is applicable to many realistic shallow-water environments with a non-uniform SSP and penetrable seabed for rays/modes interacting with both the surface and the bottom, because they behave similarly to those in ideal waveguides in terms of the dispersion. Obviously, this happens when the SSP is nearly constant over depth in winter. Since typical shallow-water environments^{19} have bottom critical angles of less than $20\xb0$, we can assume *β* = 1, effectively being treated as invariant (so is *χ*).

When the SSP is downward refracting due to warm surface water in summer,^{30} the early arrivals are refracted below the thermocline and bottom-reflected, followed by surface/bottom-reflected arrivals. The early refracted arrivals usually have a value of *β* greater than unity (e.g., 1.5),^{31–33} so it is not surprising that *β* has been modeled as a distribution rather than a single scalar.^{23,34} In such cases, we can selectively choose for source localization the later arrivals whose dispersion characteristics are similar to those in ideal waveguides.^{13,30} When a source is close to the surface (e.g., ships), it excites only high-order modes that are reflected at both boundaries, and so all arrivals [e.g., Fig. 1(b)] can be exploited assuming $\beta =cos2\theta $ (or *β* = 1 for small angles).^{14,35} It should be pointed out that the waveguide invariant typically has been used to describe the slope of the intensity striations in the frequency and range domain.^{23,28,33,36}

Now that *β* (and *χ*) is invariant to the details of the shallow-water acoustic waveguide, such as the SSP in the water column and geoacoustic bottom properties, a robust source-range estimation can be achieved through Eq. (1), together with the stability and small variation of *β* around unity in many realistic shallow-water environments ($\beta \u22481$). In case steep rays are involved, one may adopt $\beta =cos2\theta \xaf$, where $\theta \xaf$ is the averaged absolute angle of arrivals, and *β* remains invariant (so does *χ*). Alternatively, a more logical approach is to combine the two invariants ($\chi /\beta $), which will be discussed in Sec. IV.

## III. BEAM-TIME MIGRATION

The array invariant defined in Eq. (3) was the slope of the beam-time migration curve in $(t,\u2009cos\u2009\theta )$ coordinates,

To examine beam-time migration, consider an isovelocity ideal waveguide and two common arrays, a horizontal line array (HLA) and vertical line array (VLA), as illustrated in Fig. 2. For simplicity, the HLA is aligned with the source (i.e., endfire),^{21} although it can easily be extended to more general cases with an azimuth angle ($\varphi $) in the horizontal plane.^{11,37,38} The main distinction is that the HLA measures the horizontal component of the spatial wavenumber ($cos\u2009\theta $), whereas the VLA measures the vertical component ($sin\u2009\theta $). Note that beam-time migration does not require the absolute travel times of arrival because Eq. (6) involves the relative travel times, i.e., the differences of travel time (*dt*) (see Figs. 3 and 4).

### A. Horizontal line array (HLA)

Let us assume that *χ* is a constant, which occurs when *β* = 1 for small angles of arrival in shallow water or more generally when *β* is invariant (e.g., $\beta =\u22123$).^{12} The beam-time migration for a horizontal array is a straight line in $(t,\u2009cos\u2009\theta )$ coordinates from Eq. (6),

where *d _{s}* is a constant intercept. For source-range estimation via Eq. (1), the invariant parameter

*χ*needs to be estimated directly from the beam-time migration data. The LS estimate of $\chi \u0302$ can be obtained by fitting the beam angles ($cos\u2009\theta k$) and relative travel times (

*t*) of

_{k}*K*arrivals ($k=1,\u2026,K$),

where $S=[cos\u2009\theta 1,\u2009cos\u2009\theta 2,\u2026,\u2009cos\u2009\theta K]T$ and $T=[(t1,t2,\u2026,tK)T\u20091T]$, with the superscript “*T*” denoting a matrix transpose, and **1** being a *K* row vector of ones.

Figure 3 illustrates simulated beam-time migration in $(t,\u2009cos\u2009\theta )$ coordinates for two different groups of arrivals (squares) using the HLA in the ideal waveguide shown in Fig. 2: (i) $\theta \u226420\xb0$ and (ii) $\theta >20\xb0$. A broadband point source (200–900 Hz) is at a range of 2 km and a depth of 5 m, and the horizontal array at a depth of 50 m consists of 51 elements with 2-m spacing. As expected, the dashed line [Eq. (7)] with LS estimates of $\chi \u0302$ [Eq. (8)] matches well the group (i) arrivals in Fig. 3(a), since *β* is invariant (unity) for small angles.

### B. Vertical line array (VLA)

For the VLA commonly deployed in underwater acoustic experiments [see Fig. 6(a)], the beamformer output is obtained instead as a function of $sin\u2009\theta $, as presented in Fig. 4. Similar to Fig. 3, the multiple arrivals (eigenrays), circles (up-going) and squares (down-going), are selected according to their group, (i) and (ii). Again assuming *χ* is invariant, Eq. (6) leads to an ellipse in $(t,\u2009sin\u2009\theta )$ coordinates (refer to Appendix A),

whose center is located on the *t*-axis at $(tc,0)$ and $1/|\chi |$ is the horizontal semi-major axis. As shown in Fig. 4(a), the elliptic curve (dashed) passes through the group (i) small-angle arrivals with *β* = 1. The LS estimate of $\chi \u0302$ is determined from Eq. (8) with $S=[(1\u2212sin2\theta 1)1/2,(1\u2212sin2\theta 2)1/2,\u2026,(1\u2212sin2\theta K)1/2]T$.

## IV. ADAPTIVE AI, $\chi \beta $

As discussed in Sec. II C, the two invariants (*χ* and *β*) are proportional to each other and both depend on the propagation angle in ideal waveguides (i.e., $\beta =cos2\theta $). Nevertheless, up until now they have been treated as two independent invariant parameters, and the resulting beam-time migration involves *χ* alone [Eqs. (7) and (9)], referred to as the standard AI. This treatment is valid when the propagation angles are small ($<20\xb0$) with *β* = 1, as confirmed in Figs. 3(a) and 4(a) in which the standard migration (dashed) matched well the group (i) arrivals. However, the standard migration (linear or elliptic) clearly failed to match all the arrivals involving the group (ii) steep-angle arrivals in Figs. 3(b) and 4(b).

The idea of adaptive AI is to provide the beam-time migration that matches all the arrivals regardless of the propagation angle.^{38} Since *χ* and *β* are proportional to each other, we can accomplish this by using a ratio of the two parameters ($\chi /\beta $) to cancel out the angle dependence. Denoted by $\chi \beta $, the adaptive AI can be truly invariant such that

The solution to this differential equation can be written as

where $d\beta $ is a constant intercept.

For the HLA, the beam-time migration in $(t,\u2009cos\u2009\theta )$ coordinates is no longer linear, but inversely proportional to each other, consistent with Eq. (2). As evident in Fig. 3(b), the adaptive migration [solid, Eq. (11)] passes through all the arrivals (squares). The LS estimate of $\chi \u0302\beta $ is obtained from Eq. (8) with $S=\u2212[\u2009cos\u22121\theta 1,\u2009cos\u22121\theta 2,\u2026,\u2009cos\u22121\theta K]T$. The source range *r*_{0} is simply determined by $r0=\u2212c/\chi \u0302\beta $, according to Eq. (1).

For the VLA, the beam-time migration in $(t,\u2009sin\u2009\theta )$ coordinates is derived in Appendix B:

Similarly, the above adaptive migration (solid) matches all the arrivals in Fig. 4(b). The LS estimate of $\chi \u0302\beta $ is obtained from Eq. (8) with $S=\u2212[(1\u2212sin2\theta 1)\u22121/2,(1\u2212sin2\theta 2)\u22121/2,\u2026,(1\u2212sin2\theta K)\u22121/2]T$.

Table I summarizes the beam-time migration of standard AI (*χ*) with *β* = 1 and adaptive AI ($\chi \beta $) with $\beta =cos2\theta $, for the HLA (endfire) and VLA in the ideal waveguide shown in Fig. 2. For the HLA with an azimuth angle $\varphi $, the AI parameter becomes $\chi \varphi =\chi \u2009sin\u2009\varphi $.^{37} Also, $\chi \u0302$ and $\chi \u0302\beta $ are LS estimates, $ds=\u2212\chi tc$, and $d\beta =\u2212\chi \beta tc$. The effectiveness of the adaptive AI was demonstrated in a range-independent environment even with small-angle arrivals^{38} and then in a range-dependent coastal environment with a sloping bottom involving steep-angle arrivals,^{29} as described in Sec. VII. Next, we describe how to obtain the beam-time migration when the source waveform is unknown.

AI parameter . | $\chi \u2261d(cos\u2009\theta )dt$ . | $\chi \beta \u2261d(cos\u2009\theta )dt1\u2009cos2\theta $ . |
---|---|---|

HLA | $\chi t+ds=cos\u2009\theta $ | $\chi \beta t+d\beta =\u22121/\u2009cos\u2009\theta $ |

VLA | $(t\u2212tc1/\chi )2+sin2\theta =1$ | $(t\u2212tc1/\chi \beta )2\u221211\u2212sin2\theta =0$ |

Range | $r0\u0302=\u2212(c/\chi \u0302)$ | $r0\u0302=\u2212(c/\chi \u0302\beta )$ |

AI parameter . | $\chi \u2261d(cos\u2009\theta )dt$ . | $\chi \beta \u2261d(cos\u2009\theta )dt1\u2009cos2\theta $ . |
---|---|---|

HLA | $\chi t+ds=cos\u2009\theta $ | $\chi \beta t+d\beta =\u22121/\u2009cos\u2009\theta $ |

VLA | $(t\u2212tc1/\chi )2+sin2\theta =1$ | $(t\u2212tc1/\chi \beta )2\u221211\u2212sin2\theta =0$ |

Range | $r0\u0302=\u2212(c/\chi \u0302)$ | $r0\u0302=\u2212(c/\chi \u0302\beta )$ |

## V. BLIND DECONVOLUTION

The dispersion-based AI in the beam-time domain requires either an impulsive source (e.g., an explosive)^{21} or a time-domain Green's function, which is usually provided by a probe signal (i.e., a cooperating source) such as linear frequency modulation (LFM) chirp after matched-filtering.^{11,30,39} In practice, we do not know the source waveform. For broader applications of AI, we adopt the blind deconvolution^{27} to estimate the Green's function from an unknown broadband source such as a ship.^{14,35,37,40} The theoretical formulation of the ray-based blind deconvolution is reviewed here for a VLA.

Consider a point source located at $r0=(r0,z0)$ in an acoustic waveguide that broadcasts a signal *s*(*t*) whose Fourier transform is $S(\omega )=|S(\omega )|ei\Phi 0(\omega )$. The Fourier transform of the received signal on the *n*th element of a vertical array at the origin $(0,zn)$ is given by

In the ray formulation, the complex Green's function can be represented as the sum of *K* amplitudes *A _{k}* and travel times $T(\theta k)$,

where *θ _{k}* is the

*k*th angle of arrival and $\tau n(\theta k)$ is the local time delay of the

*k*th ray at

*n*th element with respect to the array center $(zc)$. To estimate the source signal and its unknown phase component $\Phi 0(\omega )$, conventional plane wave beamforming can be applied to a specific beam direction

*θ*, whose output can be written in the form

_{k}where *N* is the number of array elements. The phase of the beamformer output, $\psi (\omega ,\theta k)=\Phi 0(\omega )\u2212\omega T(\theta k)$, contains the source phase component $\Phi 0(\omega )$, which can be removed from the data $Pn(\omega )$ by a phase rotation,

The normalization of data $Pn(\omega )$ is adopted here to eliminate the source spectrum $|S(\omega )|$, which generally is not flat. The output of the blind deconvolution, $G\u0302(zn,r0,\omega )$, provides an estimate of the Green's function by an arbitrary amplitude and with a linear frequency shift $\omega T(\theta k)$, which corresponds to a time delay of $T(\theta k)$. Finally, the time-domain Green's function, $g\u0302(zn,r0,t)$, is obtained by Fourier synthesis.

As an illustrative example, a typical spectrogram of ship-radiated noise^{35} on a single element of the VLA during SAVEX15 is shown in Fig. 5(a),^{41} which will be described in Sec. VI. A 2-s window of noise in the frequency band of 200–900 Hz (denoted by a box) is selected across the array to extract the Green's function via blind deconvolution. The output of a conventional beamformer averaged over the frequency band is shown in Fig. 5(b), with a positive angle representing an up-going ray path. Six local maxima including the weak one at around $25\xb0$ are identified with the maximum-intensity angle of arrival at $\theta =\u221214.1\xb0$ (circle), to which blind deconvolution steers the beam in order to estimate the phase of the source waveform. The Green's function is then obtained by removing the source phase component along with normalization, which is shown in Fig. 5(c). It exhibits six wavefronts including the weak later arrival, each corresponding to one of the six local maxima in Fig. 5(b). In fact, the beam-time migration shown in Fig. 1(b) is the output of conventional beamforming of the Green's function in Fig. 5(c).

## VI. EXPERIMENTAL EXAMPLES: SAVEX15

This section provides examples of AI applications to source-range estimation during SAVEX15 conducted in the northeastern East China Sea^{41} using the R/V *Onnuri*, which serves as a source of opportunity. The experimental site had a nearly flat sandy bottom with a water depth of approximately 100 m (see Fig. 6). A 16-element, 56.25-m aperture vertical array (VLA) was deployed on the ocean bottom, covering about half the water column (from 25 to 81.25 m), and substantial buoyancy was provided at the top of the array (20 m) to keep it straight. The hydrophones were spaced 3.75 m apart corresponding to a half-wavelength at 200 Hz (i.e., design frequency). Interestingly, the SSP features an asymmetrical underwater sound channel with the channel axis at around 40-m depth.

During the experiment, the R/V *Onnuri* was towing two broadband sources (3–10 and 12–26 kHz) simultaneously at a depth of 50 m. However, our interest is in the ship itself radiating low-frequency noise (<1 kHz), which did not overlap with the high-frequency waveforms broadcast by the two controlled sources. The GPS ship track during a source-tow run is shown in Fig. 6(b), where the ship encircled the VLA at the origin (blue dot) counterclockwise for the 4-h period (15:45–19:35 UTC) on JD 145 (May 25). Note that the bathymetry in Fig. 6(b) reveals a diagonal strip of sand dune (reddish) and pits (bluish) east of the VLA, with the water depth varying from 95 to 110 m. Since the VLA was located on the sand dune, the sound propagation from the ship to the VLA is considered slightly upslope, except towards the end of the track, which also was on top of the sand dune.

### A. Array aperture and element spacing

AI employs conventional plane wave beamforming for beam-time migration in Sec. III and blind deconvolution in Sec. V, and used a 16-element large-aperture VLA (56.25 m) placed in an inhomogeneous sound field during SAVEX15. This is distinguished from other experiments conducted in shallow water, which used a dense short-aperture VLA for a mid- or high-frequency source: (1) 64-element 12 m for 2–3 kHz (HF97),^{30} (2) 32-element 2.8 m for 2–3.5 kHz (RADAR07),^{42} and (3) 24-element 1.2 m for 7–19 kHz (KAM11).^{39} While the large aperture allows for comparison between MFP and AI in Sec. VI B, it is appropriate to discuss the impact of the array aperture and element spacing on signal processing.^{13}

The classical beamforming usually considers a monochromatic wave, requiring the spatial sampling at a rate greater than half the wavelength of the monochromatic wave to avoid spatial aliasing. The ship noise processed in Fig. 1(b) and Fig. 5(b) was a broadband signal (200–900 Hz) above the design frequency of 200 Hz. This means the aliasing effect is present across the entire frequency band (i.e., sparse). However, the spatial aliasing occurs at different angles for different frequency components, and the integrated broadband beampattern tends to average out the sidelobes (i.e., destructive summation).^{43} Thus, the classical aliasing effect from the sparse array is diminished as the signal bandwidth increases, coupled with typically small grazing angles of arrival (i.e., towards the broadside), as evident in Fig. 5(b). In turn, the large aperture leads to the narrowing of the mainlobe and enables the separation of closely spaced beams. In summary, the spatial aliasing criterion has little importance for broadband arrays [see Fig. 7 (top) for both MFP and AI]. The signal bandwidth also determines the resolution in time in Fig. 1(b).

Moreover, conventional plane wave beamforming is applied assuming a constant sound speed (*c* = 1500 m/s), although the SSP in Fig. 6(a) varies significantly across the array aperture by as much as $\Delta c=18$ m/s. Recall from Sec. II C that we can use *β* = 1 in shallow-water environments with a non-uniform SSP for surface/bottom-reflected arrivals because their dispersion characteristics are similar to those in ideal waveguides. In fact, Fig. 5(c) reveals linear ray paths as in an ideal waveguide with a uniform sound speed averaged between the surface and the bottom. The beam-time migration in Fig. 1(b) is simply mirroring the dispersion in the waveguide, justifying the use of plane wave beamforming in a homogeneous field (free space) for AI applications in an inhomogeneous field.

### B. AI versus MFP

Figure 1 presents the case when the ship was at a range of 1.9 km, denoted by $#1$ in Fig. 6(b). First, AI needs to extract the Green's function via blind deconvolution using a 2-s window of ship noise (200–900 Hz), as illustrated in Fig. 5. It should be mentioned that the 2-s window was chosen because it corresponds to a travel distance of 3 m for the ship moving at 1.5 m/s, equivalent to the wavelength at 500 Hz (center frequency), so that the impact of source motion is negligible over the 2-s period (i.e., quasi-stationary).^{44} Conventional plane wave beamforming of the Green's function produced the beam-time migration in Fig. 1(b), which exhibited three up-going (circles) and three down-going (squares) paths constituting an elliptic curve (solid line). The LS estimate of $\chi \u0302$ then led to the source range $r1\u0302=1.9$ km by applying Eq. (1) and *β* = 1, remarkably close to the ship range.

While the differences between MFP and AI were described in Sec. I, the task of directly comparing the performance of AI and MFP side by side is quite challenging. A major obstacle to achieving a fair comparison is that accurate environmental information is required for the model-based MFP, whereas no such knowledge is required for the dispersion-based AI. Therefore, we employ conventional Bartlett MFP, which is most robust to the environmental mismatch, and incoherent averaging of single-frequency ambiguity surfaces over the frequency band, which generally works well.^{1} Here, the replica vectors are generated using the normal mode program kraken,^{45} assuming a range-independent environment with the SSP and geoacoustic parameters indicated in Fig. 6(a).

The resulting ambiguity surface of MFP in range and depth is shown in Fig. 1(a) with the main peak ($\xb0$) at $r1\u0302=1.8$ km range and 5-m depth, suggesting that our environmental model reasonably captured the actual propagation for comparison purposes. In addition, we are mainly interested in the source range since the ship is near the surface. The same 2-s ship noise in Fig. 5(a) was processed to construct the sample cross-spectral density matrix (CSDM) for individual frequency components and then incoherent averaging of the ambiguity surfaces at 20-Hz intervals across the bandwidth (200–900 Hz).

To evaluate the impact of array aperture and spacing on the performance, Fig. 1 is redrawn in Fig. 7 using only eight elements of the VLA for two different configurations as indicated on the leftmost column: (top) every other element with a similar aperture and (bottom) lower eight elements with a half aperture. Although the main peak still occurs near the source location (left column), it is clear that MFP with a smaller aperture (bottom) suffers a performance degradation with numerous peaks (red) at spurious locations close to the boundaries. This also confirms that the spatial aliasing effect from the array sparsity (top) is less significant for broadband signals. Note that both ambiguity surfaces (left) exhibit an elevated background level (greenish). On the other hand, the beam-time migration of AI (right column) both remains similar to Fig. 1(b), indicating that AI based on the plane wave beamforming is viable with a smaller aperture array for source-range estimation.

### C. Impact of array tilt

One thing AI has in common with MFP is its sensitivity to array tilt induced by ocean currents, which is a surprisingly neglected source of mismatch in the literature.^{46,47} It is not feasible for MFP to distinguish the array tilt from other environmental mismatch, let alone recognize it. Recently a robust multiple-constraint MFP^{9} was developed, which is tolerant to array tilt uncertainty as well as environmental mismatch. For AI, it was found early on that shifted beam angles could significantly impact the source range estimation, even with a small array tilt (e.g., $\Delta \theta <2\xb0$) for a 1.2-m aperture array and a high-frequency source (7–19 kHz).^{39} The sensitivity was then exploited to infer the array tilt inversely for a known source range. As it turns out,^{14} AI is capable of detecting the array tilt from its own beam-time migration, if any exists, and at the same time estimating the source range (i.e., *χ*) as if no array tilt were present.

Displayed in Fig. 8(a) is the beam-time migration when the ship was at a range of 2.2 km, denoted by $#2$ in Fig. 6(b). Similar to Fig. 1(b), it also exhibits six distinct up- and down-going arrivals (circles and squares). Clearly, they are offset downward from the elliptic curve (dashed line) assuming no array tilt, indicating the existence of array tilt. The estimated ship range is $r2\u0302=1.8$ km, which has a large error of −17%. On the other hand, the ambiguity surface in Fig. 8(b) completely fails to localize the source with a chaotic sidelobe structure. This example suggests that the array tilt is far more devastating to MFP.

To estimate the array tilt ($\Delta \theta $), the basic idea is to treat the array tilt as a nuisance parameter that should be adjusted to best fit the beam-time migration, similar to the concept of focalization.^{7} A simple optimal approach^{14} is to introduce a cost function

where *K* is the number of arrivals identified in the beam-time migration [e.g., *K* = 6 in Fig. 8(a)]. The angle $\theta \u0303k=\theta k\u2212\Delta \theta $ is the *k*th beam data (*θ _{k}*) corrected for a presumed array tilt $\Delta \theta $, from which a tilt-incorporated AI parameter $(\chi \u0303)$ can be estimated. On the other hand, $\theta \u0302k$ is the expected

*k*th beam angle along the ellipse for $\chi \u0303$.

The cost function $J(\Delta \theta )$ is the root-mean-squared (RMS) error between the two beam angles, $\theta \u0303k$ and $\theta \u0302k$, and the array tilt is estimated at the value that minimizes *J*, e.g., $\Delta \theta =3.3\xb0$ in Fig. 8(c). Here, a positive angle is defined as the array being tilted in the opposite direction of the source. The tilt-incorporated elliptic curve for $\chi \u0303$ (solid line) is shown in Fig. 8(a), which passes through all six arrivals. In addition, the horizontal line is shifted downward by the amount equivalent to the array tilt, i.e., $sin\u2009\Delta \theta $ with $\Delta \theta =3.3\xb0$. The corresponding range estimate for $\chi \u0303$ is $r2\u0302=2.2$ km, the ship GPS range. In comparison, Fig. 8(d) displays the ambiguity surface with the tilt-incorporated replicas, which also localizes the source unambiguously with the main peak ($\xb0$) at $r2\u0302=2.1$ km with −3% range error, similar to Fig. 1(a) where no array tilt was present.

### D. Performance of ship tracking

The overall performance of tracking the ship is presented in Fig. 9 for the entire 4 h during the source-tow run using AI ($\xb0$) and MFP (+): (a) in the presence of array tilt and (b) after compensation of the array tilt. The solid line represents the GPS ship range. The 2-s window of ship-radiated noise (200–900 Hz) was selected every minute to generate a total of 231 discrete time samples. Each MFP sample (+) represents the range of the main peak. The array tilt $\Delta \theta $ estimated from AI is less than $4\xb0$ as depicted in Fig. 10(a). In the presence of array tilt, both methods in Fig. 9(a) show poor performance in the region outside the gray-shaded area where $|\Delta \theta |>1.3\xb0$. It is also noticeable that the impact of array tilt is far more devastating to MFP (+), which exhibits numerous outliers that are not meaningful given the chaotic sidelobe structure observed in Fig. 8(b).

After compensation of the array tilt, the performance of both methods has drastically improved in Fig. 9(b). A few interesting observations can be made. First, the two methods are mostly comparable, provided that accurate environmental information as well as the array tilt information are fed into MFP. Second, MFP estimates (+) are consistently below the solid line (GPS) as they underestimate the ship range, except towards the end (i.e., last 30 min). Looking at the bathymetry in Fig. 6(b), the underestimation is likely attributable to the slightly upslope propagation from the ship to the VLA residing on the sand dune (reddish).^{47} This finding is supported by a slightly enhanced performance (falling on the GPS curve) over the last 30 min, denoted by a horizontal line $(\u2194)$, when the ship transited above the sand dune with a similar water depth (i.e., range-independent). Third, the MFP outliers at around 0.3 h (+) are apparently due to the pits (bluish) along the propagation path shown in Fig. 6(b), confirming its sensitivity to environmental mismatch. Last, AI equipped with self-calibration of the array tilt ($\xb0$) demonstrates a remarkable and robust performance throughout the event without using any environmental knowledge in the relatively range-independent environment, in conjunction with the blind deconvolution to estimate the Green's function from the ship noise.

### E. Time-evolving array tilt

Figure 10(a) shows the array tilt ($\Delta \theta $) estimated by AI over the entire 4-h period while the ship circled around the VLA at the origin, exhibiting an approximately smooth sinusoidal variation ($|\Delta \theta |<4\xb0$). To be precise, the apparent array tilt (treated as a rigid line) is a projection of the absolute array tilt, whose direction (heading) and inclination are unknown, onto the vertical plane along the line connecting the ship and the origin, either in the direction of the ship (negative) or away from the ship (positive). Nevertheless, at local maxima (labeled by A and E) and minimum (C), the estimated array tilt becomes the absolute tilt. At zero crossings (B and D) between them, we can at least point to the tilt direction (not the magnitude), which must be orthogonal to the vertical plane.

Figures 10(b)–10(f) include the ship's position (black circles) corresponding to each letter, where the red vectors at the origin (A, C, and E) represent both the direction and inclination of the array tilt. While the blue vectors at B and D are correct only in their direction, it is not difficult to choose one from the two orthogonal directions because the array tilt responding to the current forces (mostly tidal) changes slowly.^{48} When all five tilt vectors (three reds and two blues) are put together in Fig. 10(g), it is clear that the orientation of the array tilt is gradually turning clockwise toward the northwest, along with an increase in magnitude (from $1.7\xb0$, to $3.2\xb0$, to $3.5\xb0$). The implication is that the time-evolving array tilt estimated by AI potentially can be utilized for ocean acoustic tomography.^{49}

## VII. RANGE-DEPENDENT ENVIRONMENTS

In this section, the adaptive AI in Sec. IV is extended to mildly (adiabatic) range-dependent environments with variable bathymetry such as coastal waters,^{29} where the waveguide invariant is a complex function of the bathymetry between the source and receiver as well as the propagation angle. Assuming an ideal waveguide with variable bathymetry that is known *a priori*, a high-order approximation for $\beta (r,\theta )$ is provided in Appendix C.

Recalling Eq. (10), the adaptive AI now involves $\beta (r,\theta )$,

The solution to this differential equation can be written in the form

### A. Iterative approach

Since $\beta (r,\theta )$ requires knowledge of the source range (*r*), which is yet to be determined, the adaptive AI must be implemented iteratively with initialization of $\beta =cos2\theta $ and successive updates of *β* and *r* until the convergence of the source range.^{29,42}

- (1)
**Initialization**: We start by initializing $\beta 0=cos2\theta $. Then Eq. (19) becomes(20)$\chi \beta 0t+d\beta 0=\u22121\u2009cos\u2009\theta =F0(\theta ).$The LS estimate of $\chi \u0302\beta 0$ is obtained from Eq. (8) with $S0=[F0(\theta 1),F0(\theta 2),\u2026,F0(\theta K)]T$. The initial estimate of the source range is $r0=\u2212c/\chi \u0302\beta 0$. This completes the initialization: $\beta 0,\u2009\chi \u0302\beta 0$, and

*r*_{0}.**Recursive iteration**: At*m*th iteration, we evaluate $\beta m=\beta (rm\u22121,\theta )$ from Eq. (C1). Plugging this into the integral of Eq. (19) leads to(21)$\chi \beta mt+d\beta m=\u2212\u222b\u2009sin\u2009\theta \beta (rm\u22121,\theta )d\theta =Fm(rm\u22121,\theta ),$where the complex indefinite integral $F(r,\theta )$ can be evaluated using the matlab function,

^{50}“int().” The LS estimate of $\chi \u0302\beta m$ is obtained from Eq. (8) with $Sm=[Fm(rm\u22121,\theta 1),Fm(rm\u22121,\theta 2),\u2026,Fm(rm\u22121,\theta K)]T$. The new range estimate is then $rm=\u2212c/\chi \u0302\beta m$. This completes the*m*th iteration: $\beta m,\u2009\chi \u0302\beta m$, and*r*._{m}**Termination criterion**: The iteration stops when a further iteration provides minimal improvement in source range (typically 3 to 5 iterations).^{29}

### B. Experiment: RADAR07

RADAR07 was conducted in July 2007 on a continental shelf off the west coast of Portugal.^{42,51} To investigate the adaptive AI, we analyzed a segment of data collected during a source-tow run carried out on July 13 (JD 194), 17:37–18:18 UTC (42 min). The schematic of the experiment is shown in Fig. 11 along with the bathymetry. A VLA was moored in 87.5-m water depth, consisting of three nested subsets of hydrophones spaced for various frequency bands of interest. Here, we use a 32-element line array with a 2.8-m aperture centered around 70-m depth. A broadband source (0.5–3.5 kHz) deployed to a nominal depth of 6 m was towed by the R/V *NRP D. Carlos I* at 3 kn (1.4 m/s), which was approaching the VLA (down arrow) in Fig. 11(b). The actual bathymetry was in a bowed shape (solid line), varying from 87.5 m at the receiver to 55 m over a 5-km range. The dashed line is a simplified linear bottom (i.e., chord) whose slope is about $\varphi =0.37\xb0$. The SSP displayed in Fig. 11(a) indicates a downward-refracting environment with the mixed layer down to 10 m and a steep thermocline below. Since the source is near the surface (6 m), the surface-reflected and bottom-reflected arrivals are expected to be dominant.

During the 42-min period, the point source broadcast a 100-ms, 0.5–3.5 kHz, linear frequency modulated (LFM) chirp at 200-ms intervals continuously for 4 s (i.e., 20 chirp transmissions), which was then repeated every 30 s along the track from 4.5 km down to 0.5 km. As representative examples, the channel impulse response (CIR) obtained after matched filtering is shown in Fig. 12 (top) for two ranges: (a) 4 km (17:41 UTC) and (b) 0.5 km (18:18 UTC). The corresponding beam-time migration is displayed in Fig. 12 (bottom). The distinct multiple arrivals, marked by circles (up-going) and squares (down-going), are selected to obtain the LS estimate of $\chi \u0302$ or $\chi \beta \u0302$ iteratively. At the close range in Fig. 12(d), we notice several steep-angle arrivals (e.g., $|sin\u2009\theta |>0.3$).

The performance of source tracking during the source-tow run (42 min) is shown in Fig. 13. Each sample is an average of the 10 best results out of 20 LFM transmissions broadcast over a 4-s period, whereas the remaining 10 results include outliers due to poor matched filtering, low SNR, or possibly horizontal refraction.^{19} Solid circles (●) represent the adaptive AI ($\chi \beta \u0302$) with $\beta (r,\theta )$ computed for the actual bathymetry, Eq. (C1), whereas open circles (○) denote the standard AI ($\chi \u0302$) with $\beta (r)$ for the simplified sloping bottom (chord), Eq. (C4).

Several interesting observations can be made from Fig. 13. The performance of adaptive AI incorporating the actual bathymetry (●) is outstanding with an RMS error of only 3% using such a short-aperture array (2.8 m), without the need for environmental information except the bathymetry. In contrast, the standard AI (○) provides an RMS error of 22% which is reduced by half (11%) when we use the actual bathymetry and Eq. (C2), indicating that AI (also *β*) is sensitive to the bathymetry. In addition, the consistent overestimation of AI (above the solid line) is due to the overestimation of $\beta (r)$ in Eq. (C4) because the linearized water depth at the source (*H _{s}*) is deeper than the actual one in a bowed shape [see Fig. 11(a)]. Finally, the range error for the standard AI (○) tends to increase as the source-receiver separation decreases (or the time increases) in Fig. 13(b). On the other hand, the adaptive AI (●) maintains smaller errors with a slightly better performance toward the end, confirming that it can effectively handle the steep-angle arrivals at close ranges observed in Fig. 12(d).

## VIII. SUMMARY

This paper has reviewed the AI approach for source-range estimation in shallow water, represented by Eq. (1). Although initially thought to be something new, AI (*χ*) defined in the beam-time domain (data) is just a reflection of the modal dispersion (theory) in ocean waveguides, which is elegantly summarized by the waveguide invariant (*β*). More precisely, the two parameters (*χ* and *β*) are proportional to each other with the proportional coefficient containing the source-range (*r*) information. Since *χ* is estimated directly from the beam-time migration data collected by an array, it is crucial to understand the physics of the waveguide invariant for successful applications of AI.

The analytical expression for an ideal waveguide is $\beta =cos2\theta $ and can be approximated by *β* = 1 for small angles ($\theta <20\xb0$). Experimental data and numerical simulations have also shown that many realistic shallow-water waveguides can be well approximated by an ideal waveguide for surface/bottom-reflected arrivals because their dispersion characteristics are similar. With a small variation of *β* around unity, the robust dispersion-based AI has been successfully demonstrated in relatively range-independent shallow-water environments (HF97,^{30} KAM11,^{39} and SAVEX15^{15}) without the need for environmental knowledge such as a SSP in the water column or bottom properties, spanning various frequencies (0.2–19 kHz), ranges (0.5–6 km), and sediment types. For a range-dependent environment with variable bathymetry (RADAR07),^{29} an accurate bathymetry or bottom slope information was necessary to evaluate $\beta (r)$ using Eq. (C1) or Eq. (C4), respectively. However, AI is not practical in deep water with a refractive SSP and no bottom interaction because *β* varies widely between negative and positive values and also depends on the source range and depth.

The standard AI treats the two parameters (*χ* and *β*) as independent, which is valid as long as *β* is assumed constant, for example, in ideal waveguides for small-angle arrivals (*β* = 1). Then the beam-time migration solely based on a constant *χ* becomes linear for a horizontal array and elliptic for a vertical array. For steep-angle arrivals, the angle dependence is no longer negligible with $\beta =cos2\theta $. In such cases, we can apply the adaptive AI ($\chi \beta =\chi /\beta $), which cancels out the angle dependence and leads to a truly invariant parameter for source-range estimation from Eq. (1), $\chi \beta =\u2212c/r$. It was shown that the adaptive AI can improve the performance in a range-independent environment even with small-angle arrivals (SAVEX15) and effectively handle the steep-angle arrivals in a sloping environment at a close range (RADAR07).

AI involves simple plane wave beamforming and thus requires an array just long enough to separate multiple arrivals in the beam domain. For example, a small-aperture VLA (e.g., 2.8 m) was employed for mid-frequency signals (0.5–3.5 kHz) in Sec. VII B. On the other hand, a 16-element large-aperture VLA (56.25-m) with 3.75-m spacing (a half-wavelength at 200 Hz) was applied to ship noise (200–900 Hz) during SAVEX15, which allowed for the performance comparison between AI and MFP in Sec. VI. While the spatial aliasing effect was diminished for the broadband signal, it was also shown that a half-aperture array could provide similar performance as the full-aperture array. What was really interesting is that we were able to use plane wave beamforming even though the SSP varied significantly across the array aperture. Our rationale is that the surface/bottom-reflected arrivals exploited for AI can be viewed as propagating approximately at a depth-averaged (uniform) sound speed between the surface and bottom, since their dispersion characteristics are similar to those in ideal waveguides (*β* = 1).

The additional benefits of AI include the self-calibration of the array tilt from its own beam-time migration data, whereas we observed that the array tilt, if not compensated for, had a devastating impact on MFP. Moreover, AI was applicable to a mildly (adiabatic) range-dependent coastal environment with variable bathymetry through iterations, where the analytic expression $\beta (r,\theta )$ for a known bathymetry exists ( Appendix C).

Last, obtaining the beam-time migration requires knowledge of the time-domain Green's function between the source and the array, which can be obtained from either an explosive source or a known broadband signal broadcast by a cooperating source. However, ray-based blind deconvolution enabled the estimation of the Green's function from an unknown broadband source such as a ship, opening the door to broad applications of AI. When a source is too far away to estimate the Green's function via blind deconvolution, we can use a guide source at a close range as an intermediary by exploiting the virtual receiver concept and waveguide invariant.^{52} In the presence of multiple sources, AI can process each source independently if they occupy different frequency bands or their time-domain Green's functions are separable in time. Even when neither condition is met, i.e., similar frequency-band signals overlap in time, AI still can localize multiple sources simultaneously as long as individual Green's functions can be extracted by blind deconvolution.^{53,54}

## ACKNOWLEDGMENTS

This research was supported by the U.S. Office of Naval Research under Award No. N00014-20-1-2213.

### APPENDIX A: DERIVATION OF EQ. (9)

The integration of Eq. (6) gives

where *t*_{0} and *θ*_{0} are arbitrary constants. Rearranging the equation, we have

where

Note that $ds=\u2212\chi tc$ in Eq. (7). The above equation has the form of an ellipse whose major and minor axes coincide with the $(t,\u2009sin\u2009\theta )$ coordinates

whose center is located on the *t*-axis at $(tc,0)$ and $1/|\chi |$ is the horizontal semi-major axis.

### APPENDIX B: DERIVATION OF EQ. (12)

The integration of Eq. (10) gives

where *t*_{0} and *θ*_{0} are arbitrary constants. Rearranging the equation, we have

where

Note that $d\beta =\u2212\chi \beta tc$ in Eq. (11). Similar to Eq. (A4), Eq. (B2) can be written as

Using an infinite geometric series with $\u2009sin2\theta <1$, we can rewrite the above equation

which reduces to a parabola in $(t,\u2009sin\u2009\theta )$ coordinates when only the first term (*n* = 1) is included in the summation for $\u2009sin2\theta \u226a1$.

### APPENDIX C: BETA FOR A HOMOGENEOUS WAVEGUIDE WITH VARIABLE BATHYMETRY

Assuming that the mode propagation is adiabatic,^{19} an eighth-order approximation of the waveguide invariant $\beta (r,\theta )$ was derived in the Appendix of Ref. 29 as follows:

where *h*(*x*) is the range-dependent water depth and *H* denotes the water depth at receiver (*x* = 0). For small angles (i.e., powers of $\u2009sin2\theta \u22480$), the zeroth-order approximation is

When the ocean bottom has a constant slope such that $h(x)=H\u2212(tan\u2009\varphi )x$ with an upslope angle of $\varphi $ (see Fig. 11), Eq. (C1) can be further simplified to

where $Hs=h(r)$ is the water depth at the source range. The zeroth-order approximation Eq. (C2) also reduces to

This expression was used for a standard AI in a simplified sloping bottom.^{42}