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.

## 1. Introduction

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.

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.

*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,

*b*(

*x*) as a continuous function of range

*x*.

## 2. Related work

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.

## 3. Methods

### 3.1 Mathematical definitions

*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 \u2208 S O ( 2 )$, where

*SO*(2) is a special orthogonal group of dimension two and a position vector $ p \u2192 \u2208 R 2$. Elements of

*SE*(2) can be represented as a 3 × 3 homogeneous transformation matrix,

*se*(2) associated with

*SE*(2) is defined as

*X*using the $\u2227$ operator. The map from

*X*to $ x \u2192$ is denoted as $\u2228$, i.e., $ X \u2228 = x \u2192$ (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

*SE*(2). We also require the use of the adjoint operator for

*SE*(2),

*SE*(2). The transformation for elements of

*SE*(2) that are reflected about the line $ y = tan ( \theta g ) x + \beta $ is defined as

*m*reflections as $ R F S E 2 Total = R F S E 2 m ( R F S E 2 m \u2212 1 ( \cdots \u2009 R F S E 2 2 ( R F S E 2 1 ) ) )$ and the net inverse $ R F S E 2 Total \u2212 1$ can be expressed in the same manner.

### 3.2 Ray propagation in SE*(2)*

*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

*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

*SE*(2), we define the body-relative velocity of the ray per Wang and Chirikjian (2008) as

*dB*is a fractional Brownian motion processes parameterized by the Hurst exponent

_{H}*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 \u2212 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 $ \Delta t$ these diffusion terms are

*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

*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).

### 3.3 Reflections

*t*, $ y ( t c ) = b ( x ( t c ) )$, and the incident angle

_{c}*θ*is defined by the relationship

_{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.

_{g}*β*can be computed from

*x*(

*c*) and $ y ( t c )$ as

*m*is the number of reflection transformations on the ray. Finally, reflections from the surface are modeled using the same approach with $ \theta g = \beta = 0$.

### 3.4 Review of GPR

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 \u0303 = f ( x \u0303 ) + \eta $ where $ y \u0303$ is the observed response from input $ x \u0303$, and *η* is an independently distributed random variable with a normal distribution of variance $ \sigma p 2$. Here, $ f ( )$ is the (unknown) function to be modeled. We collect *n* observations of the function as the matrix $ Y \u0303 \u2208 R n$ and the corresponding state inputs as $ X \u0303 \u2208 R n$.

*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 ) = \alpha \xaf 2 \u2009 exp ( \u2212 | | x i \u2212 x j | | 2 / ( 2 l 2 ) )$ where $ \alpha \xaf$ and

*l*are characteristic scaling factors.

## 4. Numerical simulation results

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.

### 4.1 Geometric stochastic diffusion of ray trajectories

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 \u2212 1$ with a variance of 1 × 10^{−4} $ ( s \u2212 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.

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.

## 5. Conclusion

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).

## Supplementary Material

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

## Acknowledgments

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.

## Author Declarations

### Conflict of Interest

The authors have no conflicts of interest to declare.

## Data Availability

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.

## REFERENCES

*Harmonic Analysis for Engineers and Applied Scientists: Updated and Expanded Edition*

*Sound Propagation through the Stochastic Ocean*

*Computational Ocean Acoustics*

*A Mathematical Introduction to Robotic Manipulation*

*Gaussian Processes for Machine Learning*