GNSS-A Positioning Derivation

Consider a transponder sitting on the seafloor at some known position \(X_1\) and a transducer at the sea surface at some known position \(X_G\). Let the line-of-sight distance between the transponder and transducer be denoted by \(\vec{D}_1=X_1-X_G\), and assume the water column consists of a single layer with sound velocity \(c\). Under such conditions, the one-way travel time of an acoustic pulse \(a_1\) between the transponder and the transducer may be written as

\[\frac{||\vec{D}_1||}{c}=a_1,\]

\[\frac{\hat{D}_1}{c} \cdot \vec{D}_1=a_1\]

As a quick unit check, note that \(c\) has units of velocity (nominally m/s), \(\vec{D}_1\) has units of length (nominally meters), and \(\hat{D}_1\) is unitless, so \(a_1\) has units of time (nominally s) as expected.

Now let us assume that the seafloor transponder remains stationary but is offset from \(X_1\), while holding the transducer position \(X_G\) constant. In this case, let \(a_{mod}\) be the expected travel time had the acoustic pulse traveled along \(\vec{D}_1\) and \(a_{meas}\) be the measured travel time. If \(a_{meas} \neq a_{mod}\), then the “true” raypath \(\vec{D}_{1t}\) must differ from \(\vec{D}_1\) by offset \(\Delta X_1\). Assuming that \(||\Delta X_1|| << ||\vec{D}_1||,\) we may approximate \(\hat{D}_{1t} \approx \hat{D}_1\). This lets us write a simplified equation solving for the travel time residual \(\Delta a_1 = a_{meas}-a_{mod}\),

\[\frac{\hat{D}_1}{c} \cdot \Delta X_1 = \Delta a_1\]

If we instead consider the two-way travel time, the only difference is that the raypath is multiplied by \(2/c\) instead of \(1/c\).

Now, let us relax the constraint that \(X_G\) is constant and allow the transducer at the sea surface to move, defining the transducer position when it sends an acoustic pulse to the transponder as \(X_S\) and its position when it receives a reply as \(X_R\). During this time, the transponder on the seafloor is stationary but may still be offset from its assumed position by \(\Delta X_1\). In this case the travel time residual is the sum of the residual from the raypath \(\vec{D}_{1S}\) from the initial transducer position to the transponder and the residual from the return raypath \(\vec{D}_{1R}\) from the transponder to the final transducer position:

\[\left( \frac{\vec{D}_{1S} + \vec{D}_{1R}}{c} \right) \cdot \Delta X_1 = \Delta a_{1S} + \Delta a_{1R} = \Delta a_1\]

Until now we have only considered a single transponder, but we may generalize this to an array with \(i\) transponders. Assuming that the array moves as a block, each transponder has the same offset, so \(\Delta X_1 = \Delta X_2 = \cdots = \Delta X_i = \Delta X\). Furthermore, since there may be subtle oceanographic variations across the array, let us assume that the raypath from the transducer to each transponder travels through a different mean sound velocity \(c_i\). With this assumption, we can define the variable \(\vec{P}_i\) as

\[\vec{P}_i = \left( \frac{\vec{D}_{iS} + \vec{D}_{iR}}{c_i} \right)\]

Now, we can write a generalized equation relating the change in array position to the travel time residuals:

\[ A \cdot \Delta X = \Delta \vec{a}, \]

\[\begin{split} A = \begin{bmatrix} \vec{P}_1 \\ \vec{P}_2 \\ \vdots \\ \vec{P}_i \end{bmatrix}, \end{split}\]

\[\begin{split} \Delta \vec{a} = \begin{bmatrix} \Delta a_1 \\ \Delta a_2 \\ \vdots \\ \Delta a_i \end{bmatrix} \end{split}\]

Note that in the equations above, \(A\) is an \(i \times 3\) matrix and a has a length of \(i\).

You will no doubt have noticed at this point that I have very deliberately structured the relationship between the transponder and array offsets and the travel time residuals as a system of linear equations. This is to imply that if we measure the two-way travel time between a sea surface transducer and an array of seafloor transponders, we may infer the offset of the array from an assumed a priori position as long as we have a good estimate of the sound velocity profile across the array and the transducer position when sending and receiving acoustic pulses. This exercise forms the backbone of the GNSS-Acoustic technique, although in practice there are some more steps we must take when interpreting some data.

The observations used in this inversion are not perfect, so it is a good idea to construct a weighting matrix \(W\) which depends on the uncertainties of the travel time residuals. The travel time residuals have uncertainties dependent on three sources: the acoustic measurement uncertainty of the transducer, the position uncertainty of the transducer when the acoustic pulse is sent, and the position uncertainty of the transducer when the return acoustic pulse is received. Let \(\sigma_a^2\) be the variance of an acoustic measurement, \(C_S\) be the \(3 \times 3\) covariance matrix of the transducer position \(X_S\), and \(C_R\) be the \(3 \times 3\) covariance matrix of the transducer position \(X_R\). Assuming that \(\Delta a\) is a stochastic variable, the error propagation to derive the \(i \times i\) matrix \(W\) may be written out as:

\[ W = \sigma_a^2 I + \left( \left( \frac{\partial \Delta \vec{a}}{\partial X_S} \right)^T \cdot C_S \cdot \left( \frac{\partial \Delta \vec{a}}{\partial X_S} \right) \right) + \left( \left( \frac{\partial \Delta \vec{a}}{\partial X_R} \right)^T \cdot C_R \cdot \left( \frac{\partial \Delta \vec{a}}{\partial X_R} \right) \right) \]

Note that in the above formulation, the partial derivatives simplify to:

\[ \frac{\partial \Delta \vec{a}}{\partial X_S} = \left( \frac{\hat{D}_{1S}}{c_1}, \frac{\hat{D}_{2S}}{c_2}, \cdots , \frac{\hat{D}_{iS}}{c_i} \right), \]

\[ \frac{\partial \Delta \vec{a}}{\partial X_R} = \left( \frac{\hat{D}_{1R}}{c_1}, \frac{\hat{D}_{2R}}{c_2}, \cdots , \frac{\hat{D}_{iR}}{c_i} \right) \]

As another unit check, the elements in the partial derivatives above all have units of slowness (nominally s/m) and the elements of \(C_S\) and \(C_R\) all have units of length squared (nominally m²), so the entire matrix multiplication will result in a covariance matrix whose elements have units of time squared (nominally s²). Likewise, \(\sigma_a^2\) also has units of time squared so the units of \(W\) are consistent.

As one final note for constructing \(W\), it is possible to save some computing time by assuming that \(X_S\) and \(X_R\) are close to each other and that \(C_S \approx C_R\). In this case you only have to compute the above matrix multiplication once and can write \(W\) as

\[ W = \sigma_a^2 I + 2 \left( \left( \frac{\partial \Delta \vec{a}}{\partial X_S} \right)^T \cdot C_S \cdot \left( \frac{\partial \Delta \vec{a}}{\partial X_S} \right) \right) \]

Following this, the weighted inversion may be written out as

\[ A^T W A \cdot \Delta X = A^T W \cdot \Delta \vec{a} \]

This inversion may be solved according to the user’s preference, such as with the method of least squares.

One more consideration is that in the real world we do not collect a single travel time measurement for each transponder but rather a time series of measurements during which the seafloor transponders are stationary but the sea surface transducer moves. In this case, we can construct the \(A\), \(W\), and \(\Delta \vec{a}\) variables by calculating them at each epoch \(j\) as detailed above and then summing them prior to the inversion. Thus,

\[ \left( \sum_j A_j^T W_j A_j \right) \cdot \Delta X = \left( \sum_j A_j^T W_j \Delta \vec{a}_j \right) \]

You may also want to construct a pseudo-constraint matrix \(Q\) in addition to the above variables to perform a constrained inversion. For instance, you could construct a \(Q\) matrix to keep the baselines between transponders constant and force them to resolve the same \(\Delta X\).

In addition, the inversion may not always converge immediately. In this case you may iterate the inversion until its solution converges. Let the solution of the \(k\)th inversion be \(\Delta X_k\). Simply repeat the inversion on subsequent iterations while updating the transponder positions such that

\[ \Delta X_{ki} = \Delta X_i + \sum_0^{k-1} \Delta X_k \]

The final transponder offsets will be

\[ \Delta X = \sum_k \Delta X_k \]