Quantifying the influence of source motion on the ray-based blind deconvolution algorithm

(a) Experiment configuration of the SBCEx16 experiment, with schematic of the direct and bottom-reflected ray paths. (b) Depth averaged value of the measured sound speed profile at the experimental site. (c) Simulated ray path trajectory for a near-surface source (depth = 10 m) located 2500 m away from a short bottom-mounted VLA. (d) Simulated CIR using the experimental configuration shown in (c) in the frequency band [100 Hz–1 kHz].

Insights on the specific effects of source motion on RBD can be first gained by considering a free-space environment, where only a single straight ray path of length R exists between the center of the source track and the center of the receiver array. The source is assumed to move horizontally with constant velocity (i.e., an idealized model for a slowly moving surface ship). Here, θ_s denotes the launch angle for the kth ray (k = 1 for free space) for the case of a stationary source, whose location corresponds to the center of the track of the moving source considered hereafter [see Fig. 1(a)]. It is noted that for an arbitrary environment where refraction is present, both the launch angle θ_s and the estimated apparent DOA θ_k [e.g., as estimated by beamforming using Eq. ] depend on the selected ray path. The projection of the source velocity vector in the radial direction along the source-receiver ray path yields the radial velocity component v_r [see Fig. 1(a)]:

v_{r} = | | \vec{v} | | cos (θ_{s}) .

(4)

In free space, an approximation for the received signal on the mth element accounting for the Doppler effect induced by the moving source is¹⁰

p_{m} (t) = s (β [t + \frac{R}{c}] + β [\frac{(m - 1) Δ_{z}}{c}] sin (θ_{s})),

(5)

where the corresponding scaling factor β is determined by the radial velocity v_r and the local sound speed value c,

β = 1 + \frac{v_{r}}{c} .

(6)

Here, Δ_z represents the inter-element spacing of the VLA. Equivalently, Eq. can be expressed in the frequency domain as

P_{m} (ω) = \frac{1}{| β |} S (\frac{ω}{β}) exp (\frac{- i ω}{β c} [R + (m - 1) Δ_{z} sin (θ_{s})]) .

(7)

Under this free-space model for a moving source given by Eq. , the theoretical estimate

θ_{k} (β)

for the apparent DOA of the moving source is given by

sin (θ_{k} (β)) = β sin (θ_{s}) \to θ_{k} (β) = \sin^{- 1} (β sin (θ_{s})) .

(8)

Hence, Eq. indicates that the apparent DOA

θ_{k} (β)

for the moving source (i.e., when

β \neq 1

⁠) differs from the apparent DOA

θ_{k} (1) = θ_{s}

for a stationary source (i.e., β = 1). This difference will be investigated in more detail in Sec. 4.

2.3 Theoretical estimates for ocean waveguides

In this section, a theoretical approximation of the influence of the Doppler effect on the RBD method for the case of a moving surface source in an oceanic waveguide is proposed by analogy with the results obtained in the previous Section for the free-space case. Specifically, for an ocean waveguide, under the ray assumption, the CIR can be expressed as a sum of eigenrays linking the source to a selected receiver of the VLA [see Eq. ]. To the first order in the far field of the source, each kth eigenray, whose integrated path length is denoted

R_{k}

⁠, can be approximated as a plane wave impinging on the VLA with an apparent DOA θ_k. The launch angle for the kth ray is denoted θ_s [see Fig. 1(a)] and is defined based on the reference stationary source position (similarly to the previous subsection). Thus, by analogy with Eq. and under this plane wave approximation, the contribution of the kth ray to the received signal at the mth VLA element is proportional to

P_{m} (ω) = \frac{1}{| β_{k} |} S (\frac{ω}{β_{k}}) exp (\frac{- i ω}{β_{k} c} [R_{k} + (m - 1) Δ_{z} sin (θ_{s})]),

(9)

where the companding factor β_k for the kth ray is defined similarly to Eq. using the associated launch angle θ_s for that ray. By analogy with Eq. , it can be concluded from the construct of Eq. that the same relationship stated in Eq. exists between the launch angle θ_s and the apparent received angle θ_k at the VLA for the moving surface source, that is,

sin (θ_{k} (β_{k})) = β_{k} sin (θ_{s}) \to θ_{k} (β_{k}) = \sin^{- 1} (β_{k} sin (θ_{s})) .

(10)

Furthermore, a first-order approximation of the CIR in an ocean waveguide, under the ray assumption, for the case of a moving surface source can be obtained by analogy with the expression for the stationary CIR given by Eq. ,

G_{m} (\vec{r_{m}}, \vec{r_{s}}, ω) = \sum_{k = 1}^{K} α_{m k} e^{- i ω (T_{η k} + τ_{m k}) / β_{k}} .

(11)

It is noted that the only difference between the CIRs given by Eqs. and is the presence of the companding factor β_k which simply scales the arrival time

(T_{η k} + τ_{m k})

of the kth ray at its reception on the mth VLA element. Hence, by analogy with Eq. , the phase

ϕ_{k} (ω, θ_{k})

of the beamformer output when steered onto the apparent receiving angle θ_k [given by Eq. ] can be approximated as

ϕ_{k} (ω, θ_{k}) \approx \frac{- ω T_{η k}}{β_{k}} + ϕ_{s} (\frac{ω}{β_{k}}) .

(12)

Therefore, by analogy with Eq. , the Fourier transform of the estimated CIR from RBD

{\hat{G}}_{m} (ω, θ_{k})

for the case of moving surface source is simply given by

{\hat{G}}_{m} (ω, θ_{k}) \equiv \frac{G_{m} (ω) e^{i ω T_{η k} / β_{k}}}{\sqrt{\sum_{m} | G_{m} (ω) |^{2}}} .

(13)

It can be noticed that the only difference between the CIRs estimated from RBD given by Eq. (stationary source) and Eq. (moving source) is the presence of the companding factor β_k which simply scales the differential arrival time τ_mk of the kth ray received on the mth VLA element.

In the remainder of this Letter, the influence of the Doppler effect on the RBD algorithm will be quantified by its effect on the estimated apparent receiving angle θ_k and the estimated differential travel times of the estimated CIR from RBD. To this end, a first error metric

Δ Θ_{k}

is defined as the difference between the apparent DOA estimated from the moving source case and the one obtained under the stationary source assumption [i.e.,

{\hat{θ}}_{k} = θ_{k} (β_{k} = 1)

],

Δ Θ_{k} = θ_{k} - {\hat{θ}}_{k} .

(14)

Additionally, for an arbitrary moving source, the relative arrival times for the kth ray of the estimated CIR from RBD, denoted

γ_{m} (| \vec{v} |, R)

for each mth receiver, depend both on the moving source velocity

\vec{v}

and the source-receiver range

R

⁠. In the following section, these arrival times will be computed on the basis of the peak arrival of the envelope of the estimated CIR from RBD. For each receiver, due to the Doppler effect, the relative arrival time

γ_{m} (| \vec{v} |, R)

for the case of a moving source is expected to differ from the corresponding arrival times obtained for a stationary source

γ_{m, s} (R) = γ_{m} (| \vec{v} | = 0, R)

⁠. Therefore, a second error metric,

δ_{t} (| \vec{v} |, R)

⁠, is used to quantify the difference (averaged over all M receivers) for the estimated relative arrival times of RBD between the moving and stationary source cases,

δ_{t} (| \vec{v} |), R) = \frac{1}{M} \sum_{m = 1}^{M} \frac{(γ_{m} (| \vec{v} |, R) - γ_{m, s} (R))}{γ_{m} (| \vec{v} |, R)} .

(15)

Finally, an estimate

γ_{m, β} (| \vec{v} |, R)

of the actual Doppler-influenced relative arrival times

γ_{m} (| \vec{v} |, R)

for the moving source case can be given simply, based on Eq. , as a scaled value of the corresponding stationary source case

γ_{m, s} (R)

⁠,

γ_{m, β} (| \vec{v} |, R) = β_{k} γ_{m, s} (R) .

(16)

Hence, a third error metric,

{\tilde{δ}}_{t} (| \vec{v} |, R)

⁠, will be used in Sec. 4 to quantify the difference (averaged over all M receivers) between the estimated relative arrival times estimated for the moving source case and its theoretical approximation given by Eq. ,

{\tilde{δ}}_{t} (| \vec{v} |), R) = \frac{1}{M} \sum_{m = 1}^{M} \frac{[γ_{m} (| \vec{v} |, R) - γ_{m, β} (| \vec{v} |, R)]}{γ_{m} (| \vec{v} |, R)} .

(17)

2.4 Numerical implementation of the Doppler effect

An exact formulation for the Fourier transform of the pressure field recorded at a fixed location

\vec{r_{m}}

(representing here the location of the mth element of the VLA) by a moving source along an arbitrary track

\vec{r_{s}} (t)

while broadcasting a signal s(t) is given by the piece-wise convolution operation,⁸^,

P_{m} (ω) = \int_{- \infty}^{\infty} G_{m} (\vec{r_{m}}, \vec{r_{s}} (t), ω) \cdot s (t) e^{- i ω t} d t .

(18)

Numerically, Eq. can be implemented using a discretized formulation for the time variable,⁸

P_{m} (ω) = \sum_{n = 1}^{N} G_{m} (\vec{r_{m}}, \vec{r_{s}} (t_{n}), ω) \cdot s (t_{n}) e^{- i ω t_{n}},

(19)

where N represents the amount of sampled points of the recording. The operation in Eq. can be used to obtain an exact expression for the estimated CIR from RBD when the source of opportunity is in motion (e.g., a shipping vessel).⁹

3. Simulations of moving sources in the SBCEx16 Environment

Ray-based numerical simulations using the BellHop¹¹ program were performed to generate synthetic random waveform recordings P_m from moving sources of varying velocity. Using the environmental parameters of the SBCEx16 experiment [Fig. 1(b)], synthetic CIRs between source and receiver were generated [Fig. 1(d)] and utilized in methods as described in Sec. 2.4. In contrast to the experimental shipping noise data gathered during the experiment, ray-based simulations are free from some specific confounding factors which were preferable to avoid. Such factors include changes in source level and radiated noise characteristics due to increasing source velocity, as well as inaccurate source positioning data. During the experiment, large vessels had observed velocities (using AIS) rarely exceeding 10 m/s,⁹ an effective noise radiation depth of 10 m, and were processed in the 100–1000 Hz band with a recording duration of T = 4 s. To capture the effects of prolonged source motion, the velocity of the source was modeled up to $| \vec{v} | = 20$ m/s moving both towards and away from the receiver, with recording time up to T = 10 s. A bottom-mounted VLA containing 32 variably spaced elements was positioned to match the third VLA used during SBCEx16,⁹ spanning from 500 to 577 m depth [Fig. 1(c)]. The simulations generated CIRs between the moving source and VLA for four different source-receiver ranges: 1500, 2000, 2500, and 3000 m. These four reference source-receiver ranges served as center position for the tracks of the horizontally moving sources [Fig. 1(a)]. The results stemming from 2000 and 2500 m range simulations are emphasized in this study, as these ranges best fit the plane wave assumption for the direct ray path and clearly showed the direct path in previous experimental RBD results.⁹

4. Results

4.1 Source-motion induced error

When comparing the difference in apparent DOA for direct arrival, calculated between the moving source case and the stationary source case [denoted $Δ Θ_{k}$ according to Eq. (14)], an increase in source velocity is associated, as expected, with a nearly proportional change in the apparent DOA for all simulated ranges [Figs. 2(a), 2(b), and 3]. In particular, the greatest difference of $Δ Θ_{k}$ (obtained for a source velocity of 20 m/s), is $0.19 °$ at 2000 m and $0.12 °$ at 2500 m [Figs. 3(c) and 3(d)]. Furthermore, this change in θ_k is consistent with the change in the arrival times of the CIR estimated from RBD: higher source velocities are associated with larger differences in arrival times [as quantified by the metric δ_t, see Eq. (15)] between the stationary and moving source cases [Figs. 1(c) and 1(d)]. For example, the respective relative variation in arrival times is around 0.018 and 0.014 for a source velocity of 20 m/s in the source-receiver range of 2000 and 2500 m. Extending this analysis to greater ranges for all modeled positive velocities shows the increasing impact of source velocity and closer range on the estimated arrival times [Figs. 2(c), 2(d), and 4(e)]. At moving vessel velocities⁹ observed in the SBCEx16 experiment, the Doppler induced offset is shown to be on the order of $\approx 1 %$ [Fig. 4(e)].

Fig. 2.

(a) Difference ΔΘk [see Eq. (14)] for a range R = 2000 m between the simulated apparent DOA of for a stationary source ( θ̂k) and moving source (θk), as a function of the magnitude of the source velocity and the integration time used for RBD. (b) Same as (a) but for a of range R = 2500 m. (c) Respective offset between the averaged relative arrival time structures of the RBD CIR or a range R = 2000 m between the stationary source case and the moving source recording (δt). (d) Same as (c) but for a of range R = 2500 m.

(a) Difference $Δ Θ_{k}$ [see Eq. (14)] for a range R = 2000 m between the simulated apparent DOA of for a stationary source (⁠ ${\hat{θ}}_{k}$ ⁠) and moving source (θ_k), as a function of the magnitude of the source velocity and the integration time used for RBD. (b) Same as (a) but for a of range R = 2500 m. (c) Respective offset between the averaged relative arrival time structures of the RBD CIR or a range R = 2000 m between the stationary source case and the moving source recording (δ_t). (d) Same as (c) but for a of range R = 2500 m.

Fig. 3.

(a) Angular variations of the plane wave beamformer output for increasing source velocity for a range R = 2 km. The apparent DOA θ_k is marked by cyan dots. For comparison, the theoretically estimated DOA $θ_{k} (β_{k})$ is shown by magenta dots. (b) Same as (a) but for R = 2.5 km. (c) Variations of the error in DOA $Δ Θ_{k}$ [see Eq. (14)] between the moving source case and the stationary source case for increasing source velocity. (d) Variations of the difference between the estimated DOA θ_k and the predicted DOA $θ_{k} (β_{k})$ ⁠; i.e., the difference between the values shown by the cyan dots and red dots on (a-b)- for increasing source velocity.

Fig. 4.

$(a) Variations of the relative arrival time for the direct arrival along the VLA aperture as a function of the depth of the element, for a range R = 2500 m and a source velocity of 2 m/s. (b) Same as in (a) but for a source velocity of 10 m/s. (c) Difference of the relative arrival times of the direct path between a stationary source ( γm,s) and a moving source (γm) with a velocity of 2 m/s shown across all M = 16 VLA elements (the black line); and similarly between the same non-moving source and the theoretically predicted arrival times [ γm,β, see Eq. (16)] red line. (d) Same as in (c) but for a source velocity of 10 m/s. (e) Variations of the absolute fractional change in arrival times averaged over all M = 16 receivers δt [see Eq. (15)] between the stationary and moving source cases for the four modeled source-receiver ranges versus the increase in source velocity magnitude. (f) Same as in (e) but showing instead the variations of the fractional change δ̂t [see Eq. (17)] of arrival times between the values obtained for the moving source case and the corresponding predicted values.$

https://doi.org/10.1121/1.5003327

(a) Variations of the relative arrival time for the direct arrival along the VLA aperture as a function of the depth of the element, for a range R = 2500 m and a source velocity of 2 m/s. (b) Same as in (a) but for a source velocity of 10 m/s. (c) Difference of the relative arrival times of the direct path between a stationary source (⁠ $γ_{m, s}$ ⁠) and a moving source (γ_m) with a velocity of 2 m/s shown across all M = 16 VLA elements (the black line); and similarly between the same non-moving source and the theoretically predicted arrival times [ $γ_{m, β}$ ⁠, see Eq. (16)] red line. (d) Same as in (c) but for a source velocity of 10 m/s. (e) Variations of the absolute fractional change in arrival times averaged over all M = 16 receivers δ_t [see Eq. (15)] between the stationary and moving source cases for the four modeled source-receiver ranges versus the increase in source velocity magnitude. (f) Same as in (e) but showing instead the variations of the fractional change ${\hat{δ}}_{t}$ [see Eq. (17)] of arrival times between the values obtained for the moving source case and the corresponding predicted values.

4.2 Evaluation of the theoretical predictions from source motion

The theoretical predictions for the shift in estimated DOA and relative arrival times from RBD caused by source motion are shown to have good corrective capabilities in the same set of simulations. Under the recording time of 4 s used and the moving source velocities observed experimentally (rarely exceeding 10 m/s), the predicted DOA $θ_{k} (β_{k})$ is close in value to θ_k [Figs. 3(a) and 3(b)], and trends with the peak angular energy from the beamformer output of Eq. (2). For the highest simulated source velocity at the source-receiver range of 2000 m, their peak offset is $\pm 0.04 °$ [Fig. 3(d)], while the raw difference with an unshifted DOA stands at $\pm 0.2 °$ [Fig. 3(c)]. Similarly, peak offset at the 2500 m range stands at $\pm 0.04 °$ for the predicted DOA and $\pm 0.1 °$ for the unshifted DOA offset.

The predicted relative arrival times from the RBD CIR also show good fidelity with those obtained by performing RBD on a Doppler signal recording. The predicted peak arrival times, $γ_{m, β_{k}}$ ⁠, are shown to have a smaller discrepancy with those of the moving source CIR (γ_m) compared to arrival times without any source-motion related adjustment. Specifically, the predicted arrival times have a relative change to the moving source CIR not exceeding 0.005, or about 0.5% for the shortest range. Without a shift induced by β_k [see Eq. (13)], the relative difference of the arrival times at the highest simulated source velocities is near or even exceeds 0.01, or a 1% difference [Fig. 4(e)]. It is noted that the source-receiver range of 1500 m is an outlier, as it suffers from decreased accuracy due to residual curvature of the actual CIR wavefronts. This is caused by the short bottom-mounted VLA geometry of the SBCEx16 experiment. The three other ranges have significantly more accurate arrival time predictions of less than 0.001.

4.3 Implications for source localization and tomography

Source localization schemes determine the position of an acoustic source based on the received signals at various sensors. Such methods utilizing RBD and simulation rely on a constructed library of synthetic RBD estimated CIR samples. Each sample is assigned to the respective location points of the simulated point sources that provided the recordings. The errors for both experimental and synthetic libraries in the SBCEx16 environment are established and understood³ for source-receiver ranges with the direct path present. In particular, a mean-absolute percent error bound of $\approx 3 %$ ³ of localization is observed when synthetic CIR libraries without Doppler corrections are used to match CIRs from experimental recordings. Such error may be reduced by leveraging AIS information to shift ${\hat{G}}_{m}$ per Eq. (13), although the improvement would be limited and unable to overcome other causes of error.

The change in the CIRs estimated from RBD due to Doppler is also relevant to tomography, as small changes in arrival times are reflected by significant changes in ocean parameters such as temperature. The fractional change in the relative arrival times of the VLA due to Doppler (δ_t) as highlighted by results in Figs. 4(e) and 4(f) infer an average sound speed change of up to 1%. With a reference sound speed of 1500 m/s, the change in interpreted sound speed then leads to a change of at most $\approx 3.75 ° C$ ¹² via the TEOS-10 standard sound speed formulation, which is a significant difference.

5. Conclusions

The effects of source motion on the RBD algorithm in the direct path regime of the SBCEx16 experiment environment were investigated theoretically, and numerically through a series of moving-source simulations. An analysis of RBD performed on various modeled moving source recordings established a quasilinear relationship between the Doppler induced change in the apparent DOA, and subsequently in relative arrival times of the CIRs. Furthermore, the Doppler induced changes in the estimated CIR arrival times were found to approach a 1% bound of error at closer ranges for typical vessel speeds.

Predictions for the estimated change in apparent DOA and arrival times of the estimated CIRs due to source motion were presented in free space and generalized to ocean waveguides. They encompass a shifting term driven by the radial velocity between the source and receiver and provide a means for adjusting the DOA and CIR by utilizing known vessel velocities (via AIS or other information). In the simulations performed, the estimated DOA and CIR through RBD for large shipping source velocities observed in the SBCEx16 experiment were predicted with high accuracy. For ranges of 2 km and beyond, the motion-induced error of the CIR can subsequently be reduced to $\approx 0.1 %$ ⁠.

In the scope of tomography and source localization methods relying on RBD, these results have interesting implications. An increased understanding of the error bounds caused by motion (and their prediction for adjusting purposes) allows for more robust matched field localization schemes³ and measurement of Ocean parameters and their variability. Additional research is required to evaluate these effects for longer ranges and to evaluate the enhancements that the CIR and DOA adjustments proposed in this study have on CIR-based methods mentioned in this study.

Acknowledgments

This study was supported by the Office of Naval Research under Grant No. N00014-20-1–2416.

Author Declarations

Conflict of Interest

The authors do not have any conflict of interest to report. The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1.

G.

Byun

,

J. S.

Kim

,

C.

Cho

,

H. C.

Song

, and

S.-H.

Byun

, “

Array invariant-based ranging of a source of opportunity

,”

J. Acoust. Soc. Am.

142

(

3

),

EL286

–

EL291

(

2017

).

https://doi.org/10.1121/1.4976046

2.

S.-H.

Byun

,

C. M. A.

Verlinden

, and

K. G.

Sabra

, “

Blind deconvolution of shipping sources in an ocean waveguide

,”

J. Acoust. Soc. Am.

141

(

2

),

797

–

807

(

2017

).

https://doi.org/10.1121/10.0009083

3.

N. C.

Durofchalk

,

J.

Jin

,

H. J.

Vazquez

,

K. L.

Gemba

,

J.

Romberg

, and

K. G.

Sabra

, “

Data driven source localization using a library of nearby shipping sources of opportunity

,”

JASA Express Lett.

1

(

12

),

124802

(

2021

).

https://doi.org/10.1121/1.5052259

4.

K. L.

Gemba

,

J.

Sarkar

,

B.

Cornuelle

,

W. S.

Hodgkiss

, and

W. A.

Kuperman

, “

Estimating relative channel impulse responses from ships of opportunity in a shallow water environment

,”

J. Acoust. Soc. Am.

144

(

3

),

1231

–

1244

(

2018

).

https://doi.org/10.1121/1.5025844

5.

G.

Byun

,

C.

Cho

,

H. C.

Song

,

J. S.

Kim

, and

S.-H.

Byun

, “

Array invariant-based calibration of array tilt using a source of opportunity

,”

J. Acoust. Soc. Am.

143

(

3

),

1318

–

1325

(

2018

).

https://doi.org/10.1121/10.0000605

6.

X.

Zhang

,

N. C.

Durofchalk

,

H.

Niu

,

L.

Wu

,

R.

Zhang

, and

K. G.

Sabra

, “

Geoacoustic inversion using ray-based blind deconvolution of shipping sources

,”

J. Acoust. Soc. Am.

147

(

1

),

285

–

299

(

2020

).

https://doi.org/10.1121/1.3284548

7.

K. G.

Sabra

,

H.-C.

Song

, and

D. R.

Dowling

, “

Ray-based blind deconvolution in ocean sound channels

,”

J. Acoust. Soc. Am.

127

(

2

),

EL42

–

EL47

(

2010

).

https://doi.org/10.1121/1.1703539

8.

K. G.

Sabra

and

D. R.

Dowling

, “

Broadband performance of a time reversing array with a moving source

,”

J. Acoust. Soc. Am.

115

(

6

),

2807

–

2817

(

2004

).

https://doi.org/10.1121/10.0000919

9.

N. C.

Durofchalk

and

K. G.

Sabra

, “

Analysis of the ray-based blind deconvolution algorithm for shipping sources

,”

J. Acoust. Soc. Am.

147

(

3

),

1927

–

1938

(

2020

).

https://doi.org/10.1121/10.0009268

10.

K. L.

Gemba

,

H. J.

Vazquez

,

J.

Sarkar

,

J. D.

Tippman

,

B.

Cornuelle

,

W. S.

Hodgkiss

, and

W. A.

Kuperman

, “

Moving source ocean acoustic tomography with uncertainty quantification using controlled source-tow observations

,”

J. Acoust. Soc. Am.

151

(

2

),

861

–

880

(

2022

).

https://doi.org/10.1121/1.395269

11.

M. B.

Porter

and

H. P.

Bucker

, “

Gaussian beam tracing for computing ocean acoustic fields

,”

J. Acoust. Soc. Am.

82

(

4

),

1349

–

1359

(

1987

).