This paper describes a stochastic model of ray trajectory propagation through a medium—such as the ocean—which has an uncertain sound speed profile. We frame ray propagation as a geometric fractal Brownian motion process on the special Euclidean group of dimension two, SE(2). The framing includes diffusion parameters to describe how the stochastic rays deviate from the expected rays, and these diffusion parameters are a function of the uncertainty in the sound speed profile. We demonstrate this framing for the classical Munk profile and a double-ducted profile in the Beaufort.

The accuracy of the ray tracing method is directly constrained by the environmental uncertainty. In the ocean, a range-independent sound speed profile is often a useful simplification for the real ocean due to (unknowable) spatial and temporal variation. This paper contributes to stochastic ray tracing by modeling propagation as geometric Brownian motion on the special Euclidean group of dimension two SE(2). The special Euclidean group is the set of elements where each element encodes a unique transformation of position and orientation, with a group operator that preserves the L 2 norm and the cross-product. An overview of this new SE(2) approach to ray tracing is shown in Fig. 1 along with a comparison to the more familiar Cartesian approach.

Fig. 1.

Coordinate definitions and associated dynamics of ray tracing. Definitions for the Cartesian approach are from Jensen (2011) and used in Porter (2011); the definitions for the SE(2) approach are described later in this paper. Throughout this initial work, we assume the sound speed profile is range-independent for clarity, but it is likely that the model could be extended to range dependent profiles.

Fig. 1.

Coordinate definitions and associated dynamics of ray tracing. Definitions for the Cartesian approach are from Jensen (2011) and used in Porter (2011); the definitions for the SE(2) approach are described later in this paper. Throughout this initial work, we assume the sound speed profile is range-independent for clarity, but it is likely that the model could be extended to range dependent profiles.

Close modal

Using the new coordinate definitions, this paper reports a model for stochastic ray tracing through a medium with an uncertain sound speed profile. The model includes diffusion parameters to describe how the stochastic rays deviate from the expected rays, where the diffusion parameters are a function of the uncertainty in the sound speed profile. Uncertainty in the sound speed profile of a medium can be a result of uncertainty in measurement instruments, the measurement or interpolation process, or underlying ocean variability in space and time.

We assume the unknown true sound speed profile c(y) is a continuous function of depth y and has a continuous gradient. The expectation of the sound speed and the gradient are expressed in Eq. (1), respectively,
c * ( y ) = E [ c ( y ) ] and d d y ( c * ( y ) ) = E [ c ( y ) ] .
(1)
Section 3.4 describes a method to use Gaussian process regression (GPR) to generate a model with these properties from discrete uncertain measurements. If we consider tracing a ray in the plane, we can parameterize a bathymetric profile b(x) as a continuous function of range x.

Current models of ray propagation typically parameterize the path using planar Cartesian or cylindrical coordinates. For instance, the BELLHOP program numerically solves the eikonal equation using ray coordinates that are parameterized by path length s, and the orientation of the ray path is parameterized via a tangent vector as shown in Fig. 1 (Porter, 2011; Jensen , 2011). However, in the case of planar acoustic ray propagation there are three degrees-of-freedom, although four variables are used. The extra degree-of-freedom is eliminated by scaling the tangent vector so the system of four variables solves the eikonal equation.

A comprehensive study of stochastic ray theory (Colosi, 2016) describes models for uncertainty in the path geometry implicitly in either the traditional coordinate systems or in action-angle coordinates. Although the ray angle is sometimes used as a generalized coordinate, in the case of linearized models the variation in ray angle is considered separately from the variation in ray arrival time. Non-linear models of stochastic rays that use action-angle coordinates are derived from the Hamiltonian equation, and results from simulations of these models indicate that trajectories exhibit fractal anomalous diffusion (Brown, 1998). Anomalous diffusion with the same properties is used in this study, although with a different coordinate representation.

In this work, we use a different approach by forward propagating the ray in time as an element of SE(2). We show this new approach generates the same nominal ray trajectories as the traditional approach, but this model admits a geometric fractal Brownian motion process to represent stochastic perturbations such as the type reported in Brown (1998). We conjecture that by expressing the dynamics of ray tracing as a stochastic differential equation, new tools developed for analysis of differential motion on Lie Groups and stochastic calculus can be used for this problem.

We model a single ray of a wavefront as an element of SE(2) which departs from the source with a path parameterized by time t. The planar special Euclidean group G = S E ( 2 ) can be represented with a rotation matrix R S O ( 2 ), where SO(2) is a special orthogonal group of dimension two and a position vector p R 2. Elements of SE(2) can be represented as a 3 × 3 homogeneous transformation matrix,
g = [ R 2 × 2 p 0 1 × 2 1 ] S E ( 2 ) ,
(2)
and the group operator denoted by is matrix multiplication. The Lie algebra se(2) associated with SE(2) is defined as
X = x ̂ = [ 0 α v 1 α 0 v 2 0 0 0 ] ,
(3)
where the vector x = | v 1 v 2 α | T can be mapped to X using the operator. The map from X to x is denoted as , i.e., X = x (Murray , 2017). Every element of the group SE(2) can be obtained by using the matrix exponential map on elements of se(2). Written explicitly, this map is
g ( v 1 , v 2 , α ) = exp ( X ) = [ cos ( α ) sin ( α ) [ v 2 ( 1 + cos ( α ) ) + v 1 sin ( α ) ] / α sin ( α ) cos ( α ) [ v 1 ( 1 cos ( α ) ) + v 2 sin ( α ) ] / α 0 0 1 ] .
(4)
The differential motion of the ray is expressed in these “exponential coordinates,” or basis elements of the Lie algebra of SE(2). We also require the use of the adjoint operator for SE(2),
A d ( g ) = [ R 2 × 2 M p 0 1 × 2 1 ] , where M = [ 0 1 1 0 ] .
(5)
We also make use of a reflection operator R F S E 2 which cannot be represented by the group operator of SE(2). The transformation for elements of SE(2) that are reflected about the line y = tan ( θ g ) x + β is defined as
R F S E 2 ( S E ( 2 ) × R × R : S E ( 2 ) ) = ( [ R ( θ ) 2 × 2 p 0 1 × 2 1 ] × θ g × β : [ R ( θ f ) 2 × 2 p f 0 1 × 2 1 ] ) ,
(6)
where
| p f 1 | = [ cos ( 2 θ g ) sin ( 2 θ g ) β · sin ( 2 θ g ) sin ( 2 θ g ) cos ( 2 θ g ) 2 β · cos 2 ( θ g ) 0 0 1 ] | p 1 | and θ f = θ + 2 θ g .
(7)
The inverse of the reflection transformation is simply the original reflection applied again, e.g., R F S E 2 1 · R F S E 2 = I 3 × 3 where R F S E 2 1 = R F S E 2, and I 3 × 3 is the identity matrix of size 3. Finally, multiple reflections can be represented as a recursive sequence of reflection operations. With a slight abuse of notation, we can write the net effect R F S E 2 Total of a compilation of m reflections as R F S E 2 Total = R F S E 2 m ( R F S E 2 m 1 ( R F S E 2 2 ( R F S E 2 1 ) ) ) and the net inverse R F S E 2 Total 1 can be expressed in the same manner.
This section describes the dynamics of a ray trajectory as it travels through a medium using the new SE(2) model. The location and orientation of the ray as a function of time are described using an element of SE(2) as shown in Fig. 1. Since the sound speed profile does not depend on x(t) the expected change in orientation that is described by the eikonal equation can be reduced to Snell's law, and the differential form is
d θ = cot ( θ ) · d c ( y ) c * ( y ) .
(8)
We assumed that c(y) is continuous and the expected gradient is ( d / d y ) ( c * ( y ) ). We can augment Eq. (8) to write the three degree-of-freedom dynamics of the ray in Cartesian coordinates as
| d x d y d θ | = | c * ( y ) · cos ( θ ) · d t c * ( y ) · sin ( θ ) · d t d d y ( c * ( y ) ) · cos ( θ ) · d t | .
(9)
To express the ray dynamics in SE(2), we define the body-relative velocity of the ray per Wang and Chirikjian (2008) as
( g 1 g ̇ ) d t = | cos ( α ) d x + sin ( α ) d y sin ( α ) d x + cos ( α ) d y d α | .
(10)
We substitute Eq. (9) into Eq. (10) with α = θ to obtain an equation of the expected ray dynamics. In alignment with previous theoretical research, as discussed in Sec. 2, we model the stochastic process with the inclusion of geometric fractal diffusion terms σ 1 ( y ( t ) ) · d B H 1 and σ 2 ( y ( t ) ) · d B H 2 in the ray relative coordinates. The result is
( g 1 g ̇ ) d t = | c * ( y ) 0 d d y ( c * ( y ) ) · cos ( θ ) | d t + [ σ 1 ( y ( t ) ) 0 0 0 0 σ 2 ( y ( t ) ) · cos ( θ ) ] | d B H 1 0 d B H 2 | ,
(11)
where dBH is a fractional Brownian motion processes parameterized by the Hurst exponent H and normalized such that the process has unit variance when H = 1 / 2 (the process is normal Brownian motion). In the numerical simulation study reported later in Sec. 4.1, we set H = 2 D = 0.27, where D = 1.73 is the fractal dimension estimated by Brown (1998), and use the magnitude of the diffusion terms to capture the effect of the uncertainty in the sound speed profile. Since the numerical simulation of these equations is discrete with time step Δ t these diffusion terms are
σ 1 2 ( y ( t ) ) = V [ c * ( y ) ] · Δ t and σ 2 2 ( y ( t ) ) = V [ d d y ( c * ( y ) ) ] · Δ t ,
(12)
where V [ · ] is the variance as defined in Sec. 3.4.
Equation (11) describing ray propagation can be written as
( g 1 g ̇ ) d t = h ( g ( t ) ) d t + H ( g ( t ) ) d B H .
(13)
Hereafter, we drop the explicit dependence of the vectors h and H on g(t) for clarity. For small time steps, the transformation will be near identity, therefore the SE(2)-mean can be approximated per Chirikjian and Kyatkin (2016) as
μ ( t ) = g 0 exp ( 0 t h ( τ ) ̂ d τ ) ,
(14)
where g 0 S E ( 2 ) is the original pose of the ray. Details regarding integration on SE(2) can be found in other studies (Wang and Chirikjian, 2008; Chirikjian and Kyatkin, 2016). Although not considered in this study, if the stochastic diffusion is assumed to be simple geometric Brownian motion ( H = 1 / 2) with diffusion parameters σ1 and σ2 that do not depend on the state g, then an analytical solution for the uncertainty about the SE(2)-mean is reported in Chirikjian and Kyatkin (2016).
When a ray makes contact with the ocean bottom we first determine if any acoustic energy is reflected or if all of the energy is transmitted into the bottom medium. To make this determination we calculate the intromission angle using the physical properties of both mediums at the interface (Jensen , 2011). At time of contact tc, y ( t c ) = b ( x ( t c ) ), and the incident angle θg is defined by the relationship
tan ( θ g ) = d d x ( b ( x ( t c ) ) ) .
(15)
If θg is greater than the intromission angle, then some reflection occurs and we apply the following transformation to model reflections via the method of images.
We use the transformation for a reflection, R F S E 2, where the intercept β can be computed from x(c) and y ( t c ) as
β = y ( t c ) d d x ( b ( x ( t c ) ) ) · x ( t c ) .
(16)
The expectation of the reflection of the ray for t > t c is then
μ ( t > t c ) = R F S E 2 ( [ exp ( 0 t c ( h ( τ ) ̂ ) d τ + t c t ( S h ( τ ) ̂ ) d τ ) ] × θ g × β ) ,
(17)
where S = diag ( 1 , 1 , ( 1 ) m ) R 3 × 3 is a symmetry matrix where m is the number of reflection transformations on the ray. Finally, reflections from the surface are modeled using the same approach with θ g = β = 0.

This section provides a review of GPR for non-parametric inference. In this work we use GPR to generate a continuous model of the unknowable true sound speed profile c(y) from which we can compute an expectation of the value and the gradient as defined in Eq. (1), and determine diffusion parameters σ1 and σ2 in Eq. (12). In the case of planar ray tracing with a range-independent sound speed profile, the inference only depends on a single variable. Such functions can be modeled as a Gaussian process of the form y ̃ = f ( x ̃ ) + η where y ̃ is the observed response from input x ̃, and η is an independently distributed random variable with a normal distribution of variance σ p 2. Here, f ( ) is the (unknown) function to be modeled. We collect n observations of the function as the matrix Y ̃ R n and the corresponding state inputs as X ̃ R n.

Given a test point x ̃ * we can compute the expected value and the variance of the function at this point, f * as
E [ f * | X ̃ , Y ̃ , x * ] = K ( x * , X ̃ ) [ K ( X ̃ , X ̃ ) + σ p 2 I ] 1 Y ̃
(18)
and
V [ f * | X ̃ , x * ] = K ( x * , x * ) K ( x * , X ̃ ) [ K ( X ̃ , X ̃ ) + σ p 2 I ] 1 K ( X ̃ , x * ) ,
(19)
where K ( x * , X ̃ ) is the covariance between the test point x * and the state inputs in the training data X ̃, and K ( X ̃ , X ̃ ) is the covariance matrix for the state inputs in the training data (Williams and Rasmussen, 2006). I is the identity matrix of appropriate size. Here, we use the radial basis function to compute K, e.g., k ( x i , x j ) = α ¯ 2 exp ( | | x i x j | | 2 / ( 2 l 2 ) ) where α ¯ and l are characteristic scaling factors.
The expectation of the gradient of f * can be computed as
E [ f * | X ̃ , Y ̃ , x * ] = [ K ( x * , X ̃ ) ] D k [ K ( X ̃ , X ̃ ) + σ p 2 I ] 1 Y ̃ ,
(20)
where D k = diag ( ( x 1 x * ) / l 2 , ( x 2 x * ) / l 2 , , ( x n x * / l 2 ) ). The variance of the gradient of f * can be computed as (McHutchon, 2014)
V [ x * ( f * | X ̃ , x * ) ] = α ¯ l 2 [ K ( x * , X ̃ ) ] D k [ K ( X ̃ , X ̃ ) + σ p 2 I ] 1 D k [ K ( x * , X ̃ ) ] T .
(21)

We confirm the theoretical analysis using numerical simulation of an acoustic ray in an ocean environment. A comparison of the SE(2) method and BELLHOP program output is shown in Fig. 2. The Munk sound speed profile and a synthetic bathymetry were used as the nominal medium parameters in both cases. The BELLHOP program was set to use curvilinear interpolation and the SE(2) method used the GPR interpolation. The differences in the ray traces can be attributed to the difference in the interpolation method, and this was confirmed by numerically integrating the “traditional” ray parameters as described in Fig. 1 (Jensen , 2011) alongside the SE(2) method; both methods computed the same trajectories for each ray.

Fig. 2.

(a) Ray trace output from SE(2) model using the Munk profile. (b) Bellhop output with same inputs and Munk profile.

Fig. 2.

(a) Ray trace output from SE(2) model using the Munk profile. (b) Bellhop output with same inputs and Munk profile.

Close modal

We demonstrated the ability of this SE(2) model to capture geometric stochastic diffusion in ray trajectories. For these simulations, we used a dataset of 920 measured profiles collected and made available by the Ice-Tethered Profiler Program (Toole , 2011). The data were collected from June through September 2014 within a grid bounded by 73.5°N, 145°W to 75°N, 135°W. The uncertainty in the composite profile was assumed to be Gaussian and a GPR model was fit to reflect the mean and standard deviation of the sound speed profiles at each depth. It should be noted that the expectation of the SSP gradient below 1000 m was assumed to be 0.0167  s 1 with a variance of 1 × 10−4 ( s 1 ) 2. Figure 3(a) shows the expected value for sound speed as well as the effective sound speeds—the nominal speed plus the perturbations that would have caused the stochastic ray motions—for each Monte-Carlo simulation. The error between the effective sound speeds and the nominal sound speed across all depths is shown in Fig. 3(c). Figure 3(b) shows trajectories of stochastic rays in the water column originating at two points: one trapped inside the duct and the second below. All rays have the same number of time steps (same arrival time) and have an approximate travel distance of 40 km.

Fig. 3.

(a) Arctic sound speed profiles. (b) Computed trajectories of rays from source at x = 0 with two starting conditions: {y = −20 m, θ = 20 deg}, {y = −100 m, θ = 1 deg}. Expected ray positions are in blue, and Monte-Carlo stochastic rays are in gray. (c) Histogram of error between nominal sound speed and effective sound speed for stochastic rays.

Fig. 3.

(a) Arctic sound speed profiles. (b) Computed trajectories of rays from source at x = 0 with two starting conditions: {y = −20 m, θ = 20 deg}, {y = −100 m, θ = 1 deg}. Expected ray positions are in blue, and Monte-Carlo stochastic rays are in gray. (c) Histogram of error between nominal sound speed and effective sound speed for stochastic rays.

Close modal

The distribution of ray locations at approximately 40 km is given in Fig. 4. The distributions of Monte-Carlo rays are governed by differences in the sound speed profile and centered about the expected ray locations. An animation of the ray trajectory in Fig. 4(a) is included as a supplementary material multi-media file, and the rays in Fig. 4(b) include one reflection from the surface.

Fig. 4.

Detail views for end of rays in Fig. 3. (a) Starting location {y = −100 m, θ = 1 deg}. (b) Starting location {y = −20 m, θ = 20 deg}.

Fig. 4.

Detail views for end of rays in Fig. 3. (a) Starting location {y = −100 m, θ = 1 deg}. (b) Starting location {y = −20 m, θ = 20 deg}.

Close modal

In this paper, we presented a new model of stochastic ray tracing using geometric fractal Brownian motion on the special Euclidean group SE(2). The accuracy of this model was demonstrated with numerical simulations and compared against the BELLHOP program and the traditional ray tracing equations for the Munk profile and real Arctic profiles featuring the Beaufort Lens. Monte Carlo simulations show a ray ensemble reflective of diffusion parameters that match observed variance in the sound speed profile. Future work includes searching for a closed-form solution that describes the distribution of uncertainty about the expected ray location as a function of the diffusion. We conjecture this approach could be readily applied for tomographic inversions, extended to mediums with range-dependent sound speed profiles, and extended to three dimensions using SE(3).

See the supplementary material for an animation of the ray trajectory in Fig. 4(a).

This work was unfunded. We gratefully acknowledge publication support from the MIT Libraries Open Access Article Publication Subvention Fund. We thank our reviewers for their helpful feedback.

The authors have no conflicts of interest to declare.

The ITP data were collected and made available by the ITP Program (Toole , 2011) based at the WHOI. The website is http://www2.whoi.edu/site/itp. Relevant software for this project is available from the corresponding author upon reasonable request.

1.
Brown
,
M.
(
1998
). “
Phase space structure and fractal trajectories in 1 1 2 degree of freedom Hamiltonian systems whose time dependence is quasiperiodic
,”
Nonlin. Processes Geophys.
5
(
2
),
69
74
.
2.
Chirikjian
,
G. S.
, and
Kyatkin
,
A. B.
(
2016
).
Harmonic Analysis for Engineers and Applied Scientists: Updated and Expanded Edition
(
Courier Dover Publications
,
Mineola, NY
).
3.
Colosi
,
J. A.
(
2016
).
Sound Propagation through the Stochastic Ocean
(
Cambridge University Press
.
Cambridge, UK
).
4.
Jensen
,
F. B.
,
Kuperman
,
W. A.
,
Porter
,
M. B.
,
Schmidt
,
H.
, and
Tolstoy
,
A.
(
2011
).
Computational Ocean Acoustics
, 2nd ed. (
Springer
,
New York
).
5.
McHutchon
,
A.
(
2014
). “
Nonlinear modelling and control using Gaussian processes
,” Ph.D. thesis (
University of Cambridge
,
Cambridge, UK
).
6.
Murray
,
R. M.
,
Li
,
Z.
, and
Sastry
,
S. S.
(
2017
).
A Mathematical Introduction to Robotic Manipulation
(
CRC Press
,
Boca Raton, FL
).
7.
Porter
,
M. B.
(
2011
). “
The BELLHOP manual and user's guide: Preliminary draft
,” Report No. 260
(
Heat, Light, and Sound Research, Inc
.,
La Jolla, CA
).
8.
Toole
,
J. M.
,
Krishfield
,
R. A.
,
Timmermans
,
M.-L.
, and
Proshutinsky
,
A.
(
2011
). “
The ice-tethered profiler: Argo of the Arctic
,”
Oceanography
24
(
3
),
126
135
.
9.
Wang
,
Y.
, and
Chirikjian
,
G. S.
(
2008
). “
Nonparametric second-order theory of error propagation on motion groups
,”
Int. J. Rob. Res.
27
(
11–12
),
1258
1273
.
10.
Williams
,
C. K.
, and
Rasmussen
,
C. E.
(
2006
).
Gaussian Processes for Machine Learning
(
MIT Press
,
Cambridge, MA
), Vol. 2.

Supplementary Material