GARPOS: Analysis Software for the GNSS‐A Seafloor Positioning With Simultaneous Estimation of Sound Speed Structure

Global Navigation Satellite System–Acoustic ranging combined seafloor geodetic technique (GNSS-A) has extended the geodetic observation network into the ocean. The key issue for analyzing the GNSS-A data is how to correct the effect of sound speed variation in the seawater. We constructed a generalized observation equation and developed a method to directly extract the gradient sound speed structure by introducing appropriate statistical properties in the observation equation, especially the data correlation term. In the proposed scheme, we calculate the posterior probability based on the empirical Bayes approach using the Akaike’s Bayesian Information Criterion for model selection. This approach enabled us to suppress the overfitting of sound speed variables and thus to extract simpler sound speed field and stable seafloor positions from the GNSS-A dataset. The proposed procedure is implemented in the Python-based software “GARPOS” (GNSS-Acoustic Ranging combined POsitioning Solver).

Introduction

Basic Configurations of the GNSS‐A Observation

Precise measurements of seafloor position in the global reference frame opens the door to the “global” geodesy in the true sense of the word. It extended the observation network for crustal deformation into the ocean and has revealed the tectonic processes in the subduction zone including megathrust earthquakes (e.g., Bürgmann and Chadwell, 2014; Fujimoto, 2014, for review). Many findings have been reported especially in the northwestern Pacific along the Nankai Trough (e.g., Yokota et al., 2016; Yasuda et al., 2017; Yokota and Ishikawa, 2020), and the Japan Trench (e.g., Kido et al., 2011; Sato et al., 2011; Watanabe et al., 2014; Tomita et al., 2015; Tomita et al., 2017). These achievements owe to the development of GNSS-A (Global Navigation Satellite System–Acoustic ranging combined) seafloor positioning technique, proposed by Spiess (1980).

Observers can take various ways to design the GNSS-A observation for the positioning of the seafloor benchmark. They have to solve the difficulties not only in the technical realizations of GNSS-A subcomponents such as the acoustic ranging and the kinematic GNSS positioning, but also in designing the observation configurations and analytical models to resolve the strongly correlated parameters. For example, because the acoustic ranging observations are performed only on the sea surface, the sound speed perturbations and the depth of the benchmark are strongly correlated.

In the very first attempt for the realization, Spiess et al. (1998) derived horizontal displacement using a stationary sea-surface unit which was approximately placed on the horizontal center of the array of multiple seafloor mirror transponders. They determined the relative positions and depths of the transponders in advance. The relative horizontal positions of the sea-surface unit to the transponder array can be determined by acoustic ranging data, to be compared with its global positions determined by space geodetic technique. In this “stationary” GNSS-A configuration, the temporal variation of sound speed is less likely to affect the apparent horizontal position under the assumption that the sound speed structure is horizontally stratified. Inversely, comparing the residuals of acoustic travel time from multiple transponders, Osada et al. (2003) succeeded in estimating the temporal variation of sound speed from the acoustic data. Kido et al. (2008) modified the expression to validate the stationary configuration for a loosely tied buoy even in the case where the sound speed has spatial variations. The stationary GNSS-A configuration is applied mainly by the groups in the Scripps Institution of Oceanography (e.g., Gagnon et al., 2005; Chadwell and Spiess, 2008) and in the Tohoku University (e.g., Fujimoto, 2014; Tomita et al., 2015; Tomita et al., 2017).

On the other hand, Obana et al. (2000) and Asada and Yabuki (2001) took a “move-around” approach where the 3-dimensional position of single transponder can be estimated by collecting the acoustic data from various relay points on the sea surface. Figure 1 shows the schematic image of move-around configuration. The move-around GNSS-A configuration is developed and practicalized mainly by the collaborative group of the Japan Coast Guard and the University of Tokyo, and the Nagoya University. Unlike the stationary configuration, the horizontal positions of transponders are vulnerable to bias errors of sound speed field. Fujita et al. (2006) and Ikuta et al. (2008) then developed the methods estimating both the positions and the temporal variations of sound speed.

FIGURE 1

FIGURE 1. Schematic image of the GNSS-A system in the move-around configuration.

Similar to the effects of distribution of the GNSS satellites on the positioning, well-distributed acoustic data is expected to decrease the bias errors of the estimated transponders’ positions in the move-around configuration. By implementing the sailing observations where the sea-surface unit sails over the transponder array to collect geometrically symmetric data, positioning accuracy and observation efficiency have improved (Sato et al., 2013; Ishikawa et al., 2020).

In order to enhance the stability of positioning, an assumption that the geometry of transponder array is constant over whole observation period is usually adopted (e.g., Matsumoto et al., 2008; Watanabe et al., 2014; Chen et al., 2018; Yokota et al., 2018). Misestimates of sound speed cause the positional biases parallel to the averaged acoustic-ray direction, which results in the distortion of the estimated array geometry. Constraining the array geometry contributes to reducing the bias error in the sound speed estimates and the transponders’ centroid position.

It should be noted that these two configurations are compatible under the adequate assumptions and constraints. Recently, the group in the Tohoku University uses not only the stationary but also the move-around observation data collected for determining the array geometry (Honsho and Kido, 2017).

Recent Improvements on GNSS-A Analytical Procedures

In the late 2010s, analytical procedures with the estimation of the spatial sound speed gradient for the move-around configuration have been developed. In the earlier stage of the move-around GNSS-A development, the spatial variations of sound speed were approximated as the temporal variations, because most of the sound speed change are confined in the shallowest portion along the acoustic ray paths (e.g., Watanabe and Uchida, 2016). Actually, Yokota et al. (2019) extracted the spatial gradient of the sound speed in the shallow layer from the estimated temporal sound speed variation. However, the smoothly modeled temporal variations cannot represent the transponder-dependent variation which is caused by the sound speed gradient in the relatively deeper portion. Therefore, Yokota et al. (2019) extracted the transponder-dependent correction term from the residuals of the results derived by the conventional method of Fujita et al. (2006).

Yasuda et al. (2017) took a different approach where the sound speed structure shallower than 1,000 m is assumed to be inclined in one direction due to the Kuroshio current flowing near their sites in the offshore region south of Kii Peninsula, Japan. Because their model reflects the specific oceanographic feature, the estimated parameters are easier to be interpreted than that of Yokota et al. (2019) which has higher degree of freedom to extract the oceanographic features.

Meanwhile, Honsho et al. (2019) showed a more general expression for one-directional sound speed gradient. As they mentioned, the gradient terms in their formulation correspond to the extracted features in Yokota et al. (2019). The work by Honsho et al. (2019) showed the possibility to connect all the GNSS-A configurations into a unified GNSS-A solver. However, due to the limitation in resolving the general gradient structure, an additional constraint was taken for the practical application, which is essentially the same formulation as Yasuda et al. (2017).

In this study, to overcome the limitation above, we propose a method to directly extract the gradient sound speed structure by introducing appropriate statistical properties in the observation equation. This paper first shows the reconstructed general observation equation for GNSS-A, in which the continuity of the sound speed field in time and space is assumed. The generalized formulation approximately includes the practical solutions in the previous studies by Yasuda et al. (2017), Yokota et al. (2019), and Honsho et al. (2019) as special cases. We then describe the analytical procedure to derive the posterior probability based on the empirical Bayes approach using the Akaike’s Bayesian Information Criterion (ABIC; Akaike, 1980) for model selection. We obtain the solution which maximizes the posterior probability under the empirically selected prior distribution. This is implemented in the Python-based software “GARPOS” (GNSS-Acoustic Ranging combined POsitioning Solver; Watanabe et al., 2020a, available at https://doi.org/10.5281/zenodo.3992688).

Methodology

Positioning of Sea-Surface Transducer

The key subcomponent of the GNSS-A is the global positioning of the transducer, generally realized by GNSS observation. Whereas acoustic measurement determines the relative position of the seafloor transponder and the sea-surface transducer, GNSS plays a role to align them to the earth-centered, earth-fixed (ECEF) coordinates such as the International Terrestrial Reference Frame (ITRF). When the transducer’s position, $\mathit{P}\left(t\right)$, is determined in the GNSS’s reference frame, a realization of the ITRF, the global positions of transponders can be estimated.

It should be noted that the transponders’ positions are generally a function of time, including the solid earth tide as well as global and local crustal deformation (e.g., IERS Conventions, 2010). For the purpose of detecting crustal deformation, it is better to determine the seafloor positions in the solid-earth-tide-free coordinates. Because the observation area is limited to several-kilometers-width, solid-earth-tide-free solutions can be obtained when the trajectory of the transducer is determined in the solid-earth-tide-free coordinates. Hereafter, the positions are expressed in solid-earth-tide-free coordinates in this paper.

In order to determine $\mathit{P}\left(t\right)$ in the ECEF coordinates, a set of GNSS antenna/receiver and a gyro sensor should be mounted on the sea-surface unit. The positions of GNSS antenna, $\mathit{Q}\left(t\right)$, can be determined using any of appropriate kinematic GNSS solvers. The gyro sensor provides the attitude of the sea-surface platform, $\mathit{\Theta }\left(t\right)={\left[\begin{array}{ccc}{\theta }_r& {\theta }_p& {\theta }_h\end{array}\right]}^T$, i.e., roll, pitch, and heading (Figure 2). Because the attitude values are aligned to the local ENU coordinates, it is convenient to transform $\mathit{Q}\left(t\right)$ from ECEF to local ENU coordinates, i.e., $\mathit{Q}\left(t\right)={\left[\begin{array}{ccc}{Q}_e& Q_n& Q_u\end{array}\right]}^T$. Using the relative position of the transducer to the GNSS antenna in the gyro’s rectangular coordinate (called “ATD offset” hereafter; Figure 2), $\mathit{M}={\left[\begin{array}{ccc}{M}_r& M_p& M_h\end{array}\right]}^T$, we obtain the transducer’s position in the local ENU coordinates as,

$\begin{array}{c}\mathit{P}\left(t\right)=\mathit{Q}\left(t\right)+\mathit{R}\left(\mathit{\Theta }\left(t\right)\right)\mathit{M},\end{array}(1.1)$

with,

$\begin{array}{c}R\left(\mathit{\Theta }\right)=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& -1\end{array}\right]\left[\begin{array}{ccc}\cos {\theta }_h& -\sin {\theta }_h& 0\\ \sin {\theta }_h& \cos {\theta }_h& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos {\theta }_p& 0& \sin {\theta }_p\\ 0& 1& 0\\ -\sin {\theta }_p& 0& \cos {\theta }_p\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\theta }_r& -\sin {\theta }_r\\ 0& \sin {\theta }_r& \cos {\theta }_r\end{array}\right]\end{array}(1.2)$

The ATD offset values should be measured before the GNSS-A observation.

FIGURE 2

FIGURE 2. Definitions of the attitude parameters and the ATD offset vector for the sea-surface platform. Heading is zero when the roll axis directs to the north. The roll and pitch axes direct forward and rightward (portside) of the vessel, respectively.

Underwater Acoustic Ranging

Another key subcomponent is the technique to measure the acoustic travel time between the sea-surface transducer and the seafloor transponders. The techniques for the precise ranging using acoustic mirror-type transponders had been developed and practicalized in early studies (e.g., Spiess, 1980; Nagaya, 1995). Measuring round-trip travel time reduces the effect of advection of the media between the instruments.

The round-trip travel time for the $i$th acoustic signal to the $j$th transponder, $T_i$, is calculated as a function of the relative position of the transponder to the transducer and the 4-dimensional sound speed field, $V\left(e,n,u,t\right)$, i.e.,

$\begin{array}{c}{T}_i=T_i^c\left(\mathit{P}\left(t_{i_{+}}\right),\mathit{P}\left(t_{i_{-}}\right),{\mathit{X}}_j,V\left(e,n,u,t\right)\right),\end{array}(2)$

where $t_{i_{+}}$, $t_{i_{-}}$, and ${\mathit{X}}_j$ are the transmitted and received time for the $i$th acoustic signal, and the position of seafloor transponder numbered $j$, respectively. Note that $j$ is a function of $i$.

Although the concrete expression is provided as the eikonal equation (e.g., Jensen et al., 2011; Sakic et al., 2018), it requires much computational resources to numerically solve. When the sound speed structure is assumed to be horizontally stratified, we can apply a heuristic approach based on the Snell’s law (e.g., Hovem, 2013), which has an advantage in computation time (e.g., Chadwell and Sweeney, 2010; Sakic et al., 2018).

Therefore, we decomposed the 4-dimensional sound speed field into a horizontally stratified stational sound speed profile and a perturbation to obtain the following travel time expression:

$\begin{array}{c}{T}_i^c\left(\mathit{P}\left(t_{i_{+}}\right), \mathit{P}\left(t_{i_{-}}\right),{\mathit{X}}_j,V\left(e,n,u,t\right)\right)=\exp \left(-{\gamma }_i\right)·{\tau }_i\left(\mathit{P}\left(t_{i_{+}}\right), \mathit{P}\left(t_{i_{-}}\right),{\mathit{X}}_j,V_0\left(u\right)\right), \end{array}(3)$

where ${\tau }_i$ and $V_0\left(u\right)$ denote the reference travel time and the reference sound speed profile, respectively. $V_0\left(u\right)$ is given as a piecewise linear function of height, so that the propagation length along the radial component and the propagation time can be calculated for the given incidence angle according to the Snell's law (e.g., Hovem, 2013; Sakic et al., 2018). The expression of the correction coefficient, $\exp \left(-{\gamma }_i\right)$, is selected for the simplification in the following expansion. It represents the discrepancy ratio of the actual travel time to the reference, which caused by the spatial and temporal perturbations of the sound speed field.

In the right-hand side of Eq. 3, ${\gamma }_i$ and ${\mathit{X}}_j$ are assigned as the estimator. Eq. 1 gives the transducer’s position $\mathit{P}\left(t\right)$ as a function of the GNSS antenna’s position $\mathit{Q}\left(t\right)$, the attitude vector $\mathit{\Theta }\left(t\right)$, and the ATD offset $\mathit{M}$. The time-independent parameter $\mathit{M}$ can be also assigned as the estimator when the variation of the attitude value is large enough to resolve the parameter. Hence, the reference travel time can be rewritten as ${\tau }_i={\tau }_i\left({\mathit{X}}_j,\mathit{M}\text{|}\mathit{Q}\left(t\right),\mathit{\Theta }\left(t\right),V_0\left(u\right)\right)$, where the variables on the left and right sides of the vertical bar indicate the estimators and the observables, respectively.

Sound Speed Perturbation Model

In seawater, sound speed is empirically determined as a function of temperature, salinity, and pressure (e.g., Del Grosso, 1974). Because these variables strongly depend on the water depth, the vertical variation of the sound speed is much larger than the horizontal variation in the observation scale. Thus, $\left|{\gamma }_i\right|\ll 1$ will be satisfied in most cases where the reference sound speed appropriately represents the sound speed field. In such cases, the average sound speed along the actual ray path is expressed as $\overline{V_0}+\delta V_i \sim \overline{V_0}+{\gamma }_i\overline{V_0}$, where $\overline{V_0}$ denotes the average sound speed of the reference profile.

Recalling that the sound speed field is continuous and usually smooth in time and space within the observation scale, we can introduce a scalar field which is continuous with time and acoustic instruments’ positions, i.e., $\mathit{\Gamma }\left(t, \mathit{P},\mathit{ X}\right)$, from which the correction coefficient is extracted. Because the temporal variation of the sound speed structure is small during the travel of the acoustic signal and is usually concentrated in the shallower portion of the sea, ${\gamma }_i$ is approximated by the average of the transmission and the reception times, i.e., ${\gamma }_i\equiv \frac{1}{2}{\displaystyle \sum }_{l=i_{+},i_{-}}\Gamma \left(t_l, \mathit{P}\left(t_l\right),\mathit{ }{\mathit{X}}_j\right)$. The function $\Gamma \left(t, \mathit{P},\mathit{ X}\right)$ can be called the sound speed perturbation model.

For simplification, we put the sound speed perturbation model as time-varying linear spatial function in space as follows:

$\begin{array}{c}\Gamma \left(t, \mathit{P},\mathit{ X}\right)\equiv {\alpha }_0\left(t\right)+{\mathit{\alpha }}_1\left(t\right)·\frac{\mathit{P}}{L^{\ast }}+{\mathit{\alpha }}_2\left(t\right)·\frac{\mathit{X}}{L^{\ast }}, \end{array}(4)$

where $L^{\ast }\mathit{ }$ indicates the characteristic length of the observation site (typically in several kilometers). ${\alpha }_0\left(t\right)$, ${\mathit{\alpha }}_1\left(t\right)$ and ${\mathit{\alpha }}_2\left(t\right)$ are the time-dependent coefficients for each term. Because the vertical variation of $\mathit{P}$ and $\mathit{X}$ are much smaller than the horizontal variation, we can practically ignore the vertical component of ${\mathit{\alpha }}_1\left(t\right)$ and ${\mathit{\alpha }}_2\left(t\right)$. Thus, ${\mathit{\alpha }}_1\left(t\right)$ and ${\mathit{\alpha }}_2\left(t\right)$ are reduced to a 2-dimensional vector to denote the horizontal gradient.

Each coefficient can be represented by a linear combination of basis functions ${\mathit{\Phi }}_k\left(t\right)$:

$\begin{array}{c}\{\begin{array}{c}{\alpha }_0\left(t\right)={\displaystyle \sum }_{k=0}^{K_a}{a}_k^0{\mathit{\Phi }}_k^0\left(t\right)\\ {\mathit{\alpha }}_1\left(t\right)={\displaystyle \sum }_{k=0}^{K_b} \left(a_k^{1E}{\Phi }_k^{1E}\left(t\right),\mathit{ }{a}_k^{1N}{\mathit{\Phi }}_k^{1N}\left(t\right), 0\right)\\ {\mathit{\alpha }}_2\left(t\right)={\displaystyle \sum }_{k=0}^{K_c} \left(a_k^{2E}{\Phi }_k^{2E}\left(t\right),\mathit{ }{a}_k^{2N}{\mathit{\Phi }}_k^{2N}\left(t\right), 0\right)\end{array},\end{array}(5)$

where $a_k^{\langle ·\rangle }$ are the coefficients of the $k$th basis function, ${\mathit{\Phi }}_k^{\langle ·\rangle }\left(t\right)$, for each term named $\langle ·\rangle$. $E$ and $N$ in $\langle ·\rangle$ denote the eastward and northward components of the vector, respectively. For simplification, we compile these coefficients into vector $\mathit{a}$, hereafter.

Because the initial values for $\mathit{M}$ and ${\mathit{X}}_j$ are usually obtained in the precision of less than meters prior to the GNSS-A analysis, we approximate $\mathit{P}$ and ${\mathit{X}}_j\mathit{ }$in $\Gamma$ substituting the initial values, i.e., ${\mathit{M}}^0$ and ${\mathit{X}}_j^0$, and not updating them with the iteration. This reduces the number of estimation parameters in the correction term, i.e., ${\gamma }_i={\gamma }_i\left(\mathit{a}\text{|}{\mathit{X}}_j^0,{\mathit{M}}^0,\mathit{Q}\left(t\right),\mathit{\Theta }\left(t\right)\right)$.

Rigid Array Constraints

Usually, the local deformation within the transponders’ array is assumed to be sufficiently small, so that the same array geometry parameters can be used throughout all visits. Because the relative positions of the transponders in the array are strongly coupled with the sound speed estimates and the position of array centroid, constraining the array geometry is expected to stabilize the GNSS-A solutions. Matsumoto et al. (2008) developed the rigid-array constraint method, which has been adopted in the subsequent studies (e.g., Watanabe et al., 2014; Yokota et al., 2016) except in the cases where the rigid-array assumption is inadequate (e.g., Sato et al., 2011).

To implement the rigid-array constraint, slight change in the observation equation is needed. We divide the transponders’ positions as ${\mathit{X}}_j=\overline{{\mathit{X}}_j}+\mathit{\Delta X}$, where $\overline{{\mathit{X}}_j}$ and $\mathit{\Delta X}$ denote the relative positions of each transponder and the parallel translation of the transponder array, respectively. The array geometry, $\overline{{\mathit{X}}_j}$, should be determined prior to the analytical procedure, using the data of multiple observation visits.

Meanwhile, $\overline{{\mathit{X}}_j}$ can also be determined simultaneously with the positioning procedure by combining the data vectors, model parameter vectors, and observation equation for all series of the observation visits, as the original formulation of Matsumoto et al. (2008). However, it requires huge computational resources to solve all the parameters, as the number of observations increases. Therefore, we are not concerned with the simultaneous determination of the array geometry in the present paper.

Analytical Procedures

Observation Equation

In the GNSS-A analysis, observed travel time, $T_i^o$, are compared with the model, $T_i^c$. For the interpretability of variables and the simplification in the expansion, we took the logarithms of travel time. Summarizing the above expansion, we put the following observation equation for$i$th acoustic round-trip travel time:

$\begin{array}{c}\log \left(T_i^o/T^{\ast }\right)=\log \left({\tau }_i\left({\mathit{X}}_j,\mathit{M}\text{|}\mathit{Q},\mathit{\Theta },V_0\right)/T^{\ast }\right)-{\gamma }_i\left(\mathit{a}\text{|}{\mathit{X}}_j^0,{\mathit{M}}^0,\mathit{Q},\mathit{\Theta }\right)+e_i, \end{array}(6.1)$

or in the form with the rigid-array constraint,

$\begin{array}{c}\log \left(T_i^o/T^{\ast }\right)=\log \left({\tau }_i\left(\mathit{\Delta X},\mathit{M}\text{|}\overline{{\mathit{X}}_j},\mathit{Q},\mathit{\Theta },V_0\right)/T^{\ast }\right)-{\gamma }_i\left(\mathit{a}\text{|}{\mathit{X}}_j^0,{\mathit{M}}^0,\mathit{Q},\mathit{\Theta }\right)+e_i,\end{array}(6.2)$

where $T^{\ast }$ is the characteristic travel time and $e_i$ is the observation error vector. Figure 3 indicates the summary for constructing the observation equation. It should be noted that, in this formulation, the continuity of sound speed field is assumed.

FIGURE 3

FIGURE 3. Flow chart to construct the observation equation.

This section shows the algorithm to estimate the model parameters from the nonlinear observation equation (Eq. 6). We took a Bayesian approach because of its simple expression when incorporating prior information. Furthermore, it provides a well-defined index for the model selection, i.e., the Akaike’s Bayesian Information Criterion (ABIC; Akaike, 1980). The expansion shown in this section is based on Tarantola and Valette (1982) and Matsu’ura et al. (2007).

Prior Information

The observation equation can be rewritten as,

$\begin{array}{c}\mathit{y}=\mathit{f}\left(\mathit{x}\right)+\mathit{e},\end{array}(7)$

where $\mathit{x}={\left[\begin{array}{ccc}{\mathit{X}}_j^T& {\mathit{M}}^T& {\mathit{a}}^T\end{array}\right]}^T$, $y_i=\log \left(T_i^o/T^{\ast }\right)$, and $f_i=\log \left({\tau }_i/T^{\ast }\right)-{\gamma }_i$. Let us consider the direct prior information for the model parameters ${\mathit{X}}_j$ and $\mathit{M}$ written as,

$\begin{array}{c}\left[\begin{array}{c}{\mathit{X}}_j^0\\ {\mathit{M}}^0\end{array}\right]=\left[\begin{array}{c}{\mathit{X}}_j\\ \mathit{M}\end{array}\right]+\left[\begin{array}{c}{\mathit{d}}_X\\ {\mathit{d}}_M\end{array}\right],\end{array}(8)$

where ${\mathit{X}}_j^0$, ${\mathit{M}}^0$, and $\mathit{d}={\left[\begin{array}{cc}{\mathit{d}}_X^T& {\mathit{d}}_M^T\end{array}\right]}^T$ denote the predicted model parameter vectors and the error vector, respectively. Let us assume that ${\mathit{d}}_X$ and ${\mathit{d}}_M$ follow a normal distribution with a variance-covariance of $D_X\left({\rho }^2\right)$ and $D_M\left({\rho }^2\right)$, whose scale can be adjusted by a hyperparameter ${\rho }^2$, i.e., $D_X={\rho }^2\widetilde{D_X}$ and $D_M={\rho }^2\widetilde{D_M}$, respectively. The prior probability density function (pdf) for the constraints can be written as,

$\begin{array}{c}p\left({\mathit{X}}_j,\mathit{M};{\rho }^2\right)=c·\exp \left[-\frac{1}{2}{\left(\left[\begin{array}{c}{\mathit{X}}_j^0\\ {\mathit{M}}^0\end{array}\right]-\left[\begin{array}{c}{\mathit{X}}_j\\ \mathit{M}\end{array}\right]\right)}^T{\left[\begin{array}{cc}{D}_X\left({\rho }^2\right)& 0\\ 0& D_M\left({\rho }^2\right)\end{array}\right]}^{-1}\left(\left[\begin{array}{c}{\mathit{X}}_j^0\\ {\mathit{M}}^0\end{array}\right]-\left[\begin{array}{c}{\mathit{X}}_j\\ \mathit{M}\end{array}\right]\right)\right],\end{array}(9)$

where $c$ denotes the normalization constant.

For the model parameter $\mathit{a}$, an indirect prior information can be applied that the temporal change of sound speed perturbation model $\Gamma$ is small. Specifically, the roughness which can be defined by the derivatives of each term in Eq. 4 should be small. In this study, we use the square of second derivative as roughness $\varphi$, whereas Ikuta et al. (2008) used the first derivative. When using the B-spline functions ${\mathit{\Phi }}_k^{\langle ·\rangle }\left(t\right)$ (e.g., de Boor, 1978) as the basis of temporal sound speed variation, the roughness can be written in a vector form, i.e.,

$\begin{array}{c}\varphi ={\displaystyle \underset{t}{\int }{\left(\frac{{\partial }^2{\alpha }_{\langle ·\rangle }\left(t\right)}{\partial t^2}\right)}^{2}dt={\mathit{a}}^{{\langle ·\rangle }^T}{H}_{\langle ·\rangle }{\mathit{a}}^{\langle ·\rangle }},\end{array}(10.1)$

where,

$\begin{array}{c}{H}_{{\langle ·\rangle }_{ij}}={\displaystyle \int \frac{{\partial }^2{\mathit{\Phi }}_i^{\langle ·\rangle }\left(t\right)}{\partial t^2}\frac{{\partial }^2{\mathit{\Phi }}_j^{\langle ·\rangle }\left(t\right)}{\partial t^2}dt.}\end{array}(10.2)$

Then, the prior pdf can be written using the hyperparameter ${\lambda }_{\langle ·\rangle }$ as,

$\begin{array}{c}p\left({\mathit{a}}^{\langle ·\rangle };{\lambda }_{\langle ·\rangle }^2\right)=c·\exp \left[-\frac{1}{2{\lambda }_{\langle ·\rangle }^2}{\mathit{a}}^{{\langle ·\rangle }^T}{H}_{\langle ·\rangle }{\mathit{a}}^{\langle ·\rangle }\right],\end{array}(11)$

where $c$ denotes the normalization constant.

Combining these prior informations, we obtain the following prior pdf:

$\begin{array}{c}p\left(\mathit{x};{\rho }^2,{\mathit{\lambda }}^2\right)={\left(2\pi \right)}^{-\frac{g}{2}}{\Vert {\Lambda }_G\Vert }^{\frac{1}{2}}\exp \left[-\frac{1}{2}{\left({\mathit{x}}^0-\mathit{x}\right)}^{T}G\left({\rho }^2, {\mathit{\lambda }}^2\right)\left({\mathit{x}}^0-\mathit{x}\right)\right], \end{array}(12.1)$

with ${\mathit{\lambda }}^2=\left[\begin{array}{ccc}{\lambda }_0^2& {\lambda }_{1E}^2& \begin{array}{ccc}{\lambda }_{1N}^2& {\lambda }_{2E}^2& {\lambda }_{2N}^2\end{array}\end{array}\right]$, ${\mathit{x}}^0={\left[\begin{array}{ccc}{\mathit{X}}_j^{0^T}& {\mathit{M}}^{0^T}& {\mathbf{0}}^T\end{array}\right]}^T$, and,

$\begin{array}{c}G\left({\rho }^2, {\mathit{\lambda }}^2\right)=\left[\begin{array}{ccc}\begin{array}{c}{D}_X{\left({\rho }^2\right)}^{-1}\\ \\ \begin{array}{c}\\ \\ \begin{array}{c}\\ \\ \end{array}\end{array}\end{array}& \begin{array}{c}\\ D_M{\left({\rho }^2\right)}^{-1}\\ \begin{array}{c}\\ \\ \begin{array}{c}\\ \\ \end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\\ \\ \begin{array}{c}{H}_0/{\lambda }_0^2\\ \\ \begin{array}{c}\\ \\ \end{array}\end{array}\end{array}& \begin{array}{c}\\ \\ \begin{array}{c}\\ H_{1E}/{\lambda }_{1E}^2\\ \begin{array}{c}\\ \\ \end{array}\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}\\ \\ \begin{array}{c}\\ \\ \begin{array}{c}{H}_{1N}/{\lambda }_{1N}^2\\ \\ \end{array}\end{array}\end{array}& \begin{array}{c}\\ \\ \begin{array}{c}\\ \\ \begin{array}{c}\\ H_{2E}/{\lambda }_{2E}^2\\ \end{array}\end{array}\end{array}& \begin{array}{c}\\ \\ \begin{array}{c}\\ \\ \begin{array}{c}\\ \\ H_{2N}/{\lambda }_{2N}^2\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right],\end{array}(12.2)$

where $g$ and $\left|\right|{\Lambda }_G\left|\right|$ represent the rank of $G$ and the absolute value of the product of non-zero eigenvalues of $G$, respectively.

Variance-Covariance of Data Error

Now for the observed data, we take the assumption that $\mathit{e}$ also follows a normal distribution with a variance-covariance of ${\sigma }^{2}E$, i.e.,

$\begin{array}{c}p\left(\mathit{y}\text{|}\mathit{x};{\sigma }^2\right)={\left(2\pi {\sigma }^2\right)}^{-\frac{n}{2}}{\left|E\right|}^{-\frac{1}{2}}\exp \left[-\frac{1}{2{\sigma }^2}{\left(\mathit{y}-\mathit{f}\left(\mathit{x}\right)\right)}^T{E}^{-1}\left(\mathit{y}-\mathit{f}\left(\mathit{x}\right)\right)\right],\end{array}(13)$

where $n$ is the number of data and $|·|$ denotes the determinant of the matrix.

The major error sources for the measurement and calculation of travel time are 1) measurement error when reading the return signal, 2) transducer’s positioning error, and 3) modeling error of the sound speed field. Non-diagonal components of $E$ are caused not by measurement error, but by transducer’s positioning error and sound speed modeling error. The transducer’s positioning error may have temporal correlation which comes from the kinematic GNSS noise. The modeling error has spatio-temporal correlation because the sound speed variation is modeled by a smooth function of space and time. Thus, we assumed the following covariance terms using two hyperparameters, i.e., ${\mu }_t$ and ${\mu }_{MT}$, to adjust the non-diagonal component of $E$:

$\begin{array}{c}{E}_{ij}\left({\mu }_t,{\mu }_{MT}\right)=\{\begin{array}{c}\sqrt{E_{ii}{E}_{jj}}\exp \left(-\frac{\left|t_i-t_j\right|}{{\mu }_t}\right) if the transponders for i and j are the same \\ {\mu }_{MT}\sqrt{E_{ii}{E}_{jj}}\exp \left(-\frac{\left|t_i-t_j\right|}{{\mu }_t}\right) for others\end{array},\end{array}(14)$

whose formulation refers to Fukahata and Wright (2008). Eq. 14 means that the densely sampled data would have smaller weights in the model than the isolated data. A factor ${\mu }_{MT}\in \left[0, 1\right]$ was introduced to surpress the error correlation between the different transponders because the acoustic rays for different transponders take separate paths as the depths increases. Consideration of the non-diagonal components of the data variance-covariance contributes to reduce the complexity of the model against the excessively high-rate data sampling.

Subsequently, we consider the diagonal component of $E$ which controls the weight of individual data. Because the measurement errors of acoustic travel time are caused by mis-reading of the return signal, it is independent on the travel time value. Therefore, it is reasonable to put the assumption that the error of $T_i^o-T_i^c$ follows a normal distribution. Nonetheless, because the GNSS-A typically gives the precision of $T_i^o/T_i^c \sim 1\pm {10}^{-4}$, we can suppose that $T_i^o/T_i^c$ approximately follows a log-normal distribution as assumed in Eq. 13. In order to put the same weight for all measured travel time in the real scale, we applied $E_{ii}={\left(T^{\ast }/T_i^o\right)}^2$ for scaling the diagonal component.

Posterior Probability

The posterior pdf after the data acquisition, which can be defined to be equal to the likelihood of the model parameter with the given data, can be written as,

$\begin{array}{c}p\left(\mathit{x};{\sigma }^2,{\mu }_t,{\mu }_{MT},{\rho }^2,{\mathit{\lambda }}^2\text{|}\mathit{y}\right)=c·{\left(2\pi {\sigma }^2\right)}^{-\frac{\left(n+g\right)}{2}}{\left|E\right|}^{-\frac{1}{2}}{\Vert \widetilde{{\Lambda }_G}\Vert }^{\frac{1}{2}}\exp \left[-\frac{1}{2{\sigma }^2}s\left(\mathit{x}\right)\right],\end{array}(15.1)$

with,

$\begin{array}{c}s\left(\mathit{x}\right)={\left(\mathit{y}-\mathit{f}\left(\mathit{x}\right)\right)}^T{E}^{-1}\left(\mathit{y}-\mathit{f}\left(\mathit{x}\right)\right)+{\left({\mathit{x}}^0-\mathit{x}\right)}^T\tilde{G}\left({\mathit{x}}^0-\mathit{x}\right),\end{array}(15.2)$

where $\tilde{G}={\sigma }^{2}G\left({\rho }^2, {\mathit{\lambda }}^2\right)$ and $\Vert \widetilde{{\Lambda }_G}\Vert$ represents the absolute value of the product of non-zero eigenvalues of $\tilde{G}$.

Defining $\widehat{\mathit{x}}\left({\sigma }^2,{\mu }_t,{\mu }_{MT},{\rho }^2,{\mathit{\lambda }}^2\right)$ as $\mathit{x}$ that maximizes the posterior probability (Eq. 15) under the given hyperparameters, the partial derivative of $p\left(\mathit{x}\text{|}\mathit{y}\right)$ with respect to $\mathit{x}$ should be zero for $\mathit{x}=\widehat{\mathit{x}}$. Hence, $\widehat{\mathit{x}}$ should satisfy the following equation:

$\begin{array}{c}A{\left(\widehat{\mathit{x}}\right)}^T{E}^{-1}\left(\mathit{y}-\mathit{f}\left(\widehat{\mathit{x}}\right)\right)+\tilde{G}\left({\mathit{x}}^0-\widehat{\mathit{x}}\right)=0,\end{array}(16.1)$

where $A\left(\mathit{x}\right)$ is the Jacobian matrix at point $\mathit{x}$ defined as,

$\begin{array}{c}A\left(\mathit{x}\right)=\left[\begin{array}{ccc}\frac{\partial f_1}{\partial x_{k1}}\left(\mathit{x}\right)& \cdots & \frac{\partial f_1}{\partial x_{km}}\left(\mathit{x}\right)\\ ⋮& \ddots & ⋮\\ \frac{\partial f_n}{\partial x_{k1}}\left(\mathit{x}\right)& \cdots & \frac{\partial f_n}{\partial x_{km}}\left(\mathit{x}\right)\end{array}\right].\end{array}(16.2)$

We can solve the nonlinear equation (Eq. 16) numerically by performing an iterative method, where ${\mathit{x}}_k$ is corrected in each step with the following algorithm:

$\begin{array}{c}{\mathit{x}}_{k+1}={\mathit{x}}_k+{\left(A{\left({\mathit{x}}_k\right)}^T{E}^{-1}A\left({\mathit{x}}_k\right)+\tilde{G}\right)}^{-1}\left(A{\left({\mathit{x}}_k\right)}^T{E}^{-1}\left(\mathit{y}-\mathit{f}\left({\mathit{x}}_k\right)\right)+\tilde{G}\left({\mathit{x}}^0-{\mathit{x}}_k\right)\right),\end{array}(17)$

to satisfy the following convergence criteria:

$\begin{array}{c}A{\left({\mathit{x}}_k\right)}^T{E}^{-1}\left(\mathit{y}-\mathit{f}\left({\mathit{x}}_k\right)\right)+\tilde{G}\left({\mathit{x}}^0-{\mathit{x}}_k\right)\ll 1.\end{array}(18)$

Ignoring the term $O\left({\left(\mathit{x}-\widehat{\mathit{x}}\right)}^2\right)$ in $\mathit{f}\left(\mathit{x}\right)$ around $\widehat{\mathit{x}}$, $s\left(\mathit{x}\right)$ can be rewritten as,

$\begin{array}{c}s\left(\mathit{x}\right) \sim s\left(\widehat{\mathit{x}}\right)+{\left(\mathit{x}-\widehat{\mathit{x}}\right)}^T\left(A{\left(\widehat{\mathit{x}}\right)}^T{E}^{-1}A\left(\widehat{\mathit{x}}\right)+\tilde{G}\right)\left(\mathit{x}-\widehat{\mathit{x}}\right).\end{array}(19)$

Therefore, the linearized variance-covariance matrix around $\widehat{\mathit{x}}$ can be obtained as,

$\begin{array}{c}\widehat{C}={\sigma }^2{\left(A{\left(\widehat{\mathit{x}}\right)}^T{E}^{-1}A\left(\widehat{\mathit{x}}\right)+\tilde{G}\right)}^{-1}.\end{array}(20)$

Hyperparameter Tuning

The appropriate values of the hyperparameters can be determined by minimizing Akaike’s Bayesian Information Criteria (ABIC; Akaike, 1980),

$\begin{array}{c}ABIC=-2\log {\displaystyle \int p\left(\mathit{y}\text{|}\mathit{x};{\sigma }^2,{\mu }_t,{\mu }_{MT}\right)p\left(\mathit{x};{\rho }^2,{\mathit{\lambda }}^2\right)d\mathit{x}+2N_{\text{HP}}}\end{array},(21)$

where $N_{\text{HP}}$ denotes the number of hyperparameters. Although it is difficult to analytically calculate the integral for the marginal likelihood because of the nonlinearity in $\mathit{f}\left(\mathit{x}\right)$, the Laplace’s method can be applied in this case where the degree of freedom is sufficiently large and $s\left(\mathit{x}\right)$ can be almost unimodal. Thus, an approximated form for ABIC is obtained as follows:

$\begin{array}{c}ABIC\cong \left(n+g-m\right)\log s\left(\widehat{\mathit{x}}\right)-\log \left|E^{-1}\right|-\log \Vert {\Lambda }_G\Vert +\log \left|A{\left(\widehat{\mathit{x}}\right)}^T{E}^{-1}A\left(\widehat{\mathit{x}}\right)+\tilde{G}\right|+const.\end{array}(22)$

where $m$ is the number of model parameters. For the derivation, we used the following relationship:

$\begin{array}{c}{\sigma }^2=\frac{s\left(\widehat{\mathit{x}}\right)}{n+g-m},\end{array}(23)$

which is derived from the condition that the partial derivative of ABIC with respect to ${\sigma }^2$ should be zero. We can tune the hyperparameters to minimize the approximated ABIC value defined in Eq. 22, to obtain the solution ${\mathit{x}}^{\mathit{*}}=\widehat{\mathit{x}}\left({\sigma }^{2^{\ast }},{\mu }_t^{\ast },{\mu }_{MT}^{\ast },{\rho }^{2^{\ast }},{\mathit{\lambda }}^{2^{\ast }}\right)$, where $\ast$ denotes the selected hyperparameters.

Features of “GARPOS”

GARPOS (Watanabe et al., 2020a; available at https://doi.org/10.5281/zenodo.3992688) has been developed to implement the GNSS-A analysis procedure. GARPOS is compatible with Python 3, with other packages NumPy, SciPy, pandas, and matplotlib. These packages are pre-installed in most of the Python distributions such as Anaconda. Sample scripts and data for testing GARPOS are also stored in the repository.

GARPOS is distributed as a series of files, which requires a driver script to run. The toolset consists of multiple Python files and a Fortran 90 library for ray tracing. GARPOS requires the following input files:

(I-1) Initial site parameter file (in Python’s configuration format),

(I-2) Acoustic observation data file (in csv format),

(I-3) Reference sound speed data file (in csv format),

(I-4) Setting file (in Python’s configuration format).

Initial site parameter file (I-1) contains the initial values of the transponders’ positions, the ATD offset and the relevant prior covariance information, as well as the metadata for the observation site and conditions. Acoustic observation data file (I-2) contains the list of the observation data associated with each acoustic ranging, such as travel time, positions, attitude and other metadata. Reference sound speed data file (I-3) contains the reference sound speed profile approximated into a polygonal curve. Setting file (I-4) contains the parameters to control the analysis procedures including the hyper parameters. Users can put the lists of candidates of hyperparameters in which the best combination may be within. The parameters $nmp0$, $nmp1$, and $nmp2$ in the setting file control the number of basis functions, $K_a$, $K_b$, and $K_c$ in Eq. 5.

The results are written in the following output files:

(O-1) Estimated site parameter files (in Python’s configuration format),

(O-2) Modified acoustic observation data file (in csv format),

(O-3) Model parameter list file (in csv format),

(O-4) Posterior variance-covariance matrix file (in csv format).

Estimated site parameter files (O-1) is written in the same format as the file (I-1). Modified acoustic observation data file (O-2) contains the calculated travel time data and the coefficients of sound speed perturbation model, as well as the original data/metadata set in (I-2). Model parameter list file (O-3) and posterior variance-covariance matrix file (O-4) contain the whole estimated model parameter vector and its variance-covariance, respectively.

Major input/output parmeters and hyperparameters for GARPOS are listed in Tables 1 and 2, respectively.

TABLE 1

TABLE 1. List of observable and estimation parameters used in GARPOS.

TABLE 2

TABLE 2. List of hyperparameter in GARPOS.

We developed GARPOS to be compatible with both observation configurations. When handling the GNSS-A data collected in the stationary configurations, we should process data with some constraints on model parameters. Specifically, 1) upward components of transponders’ positions should be fixed to zero, and 2) spatial gradient components of the sound speed perturbation model should not be solved, i.e., $nmp1=nmp2=0$, because these parameters cannot be well resolved in the stationary configuration. Although further parameter tuning may be required for optimization, users can solve the seafloor position by GARPOS with the stationary data in addition to the move-around data.

Applications to the Actual Data

Data and Settings

In order to verify the proposed analytical procedure, we reanalyzed the GNSS-A data at the sites named “TOS2” and “MYGI” (Table 3; Figure 4) in 2011–2019. The test sites were selected for several reasons: 1) whereas TOS2 is expected to move at almost constant rate, MYGI will show the transient displacement due to the postseismic crustal deformation of the 2011 Tohoku-oki earthquake; 2) the oceanographic environments are different, i.e., the effect of the Kuroshio current is dominant at TOS2; but 3) the depths of both sites are almost the same. The observation epochs used in this study is listed in Supplementary Tables S1 and S2. The datasets used in this study are available at https://doi.org/10.5281/zenodo.3993912 (Watanabe et al., 2020b).

TABLE 3

TABLE 3. Locations and observation periods of the GNSS-A observation sites used in this study.

FIGURE 4

FIGURE 4. Locations of the tested GNSS-A sites TOS2 and MYGI.

Acoustic round-trip travel times were measured on the survey vessel using the hull-mounted acoustic transducer (e.g., Ishikawa et al., 2020). Processing delays in the acoustic devices were subtracted from the acoustic data beforehand.

Solid-earth-tide-free positions of GNSS antenna $\mathit{Q}\left(t\right)$ were determined at 2 Hz by the open source software RTKLIB version 2.4.2 (Takasu, 2013) in post-processing kinematic Precise Point Positioning (PPP) mode, using the precise satellite orbit and the 30-s satellite clock solutions (final products) provided by the International GNSS Service (International GNSS Service, a; International GNSS Service, b), in the same procedures as Watanabe et al. (2020c). The ATD offset values for each vessel, $\mathit{M}$, were measured by leveling, distance, and angle surveys before the first GNSS-A observation cruise, to be used as ${\mathit{M}}^0$.

Along with the acoustic observations, the profiles of temperature and/or conductivity were measured by CTD, XCTD or XBT probes several times. The reference sound speed profile, $V_0\left(u\right)$, was calculated from the observed temperature and salinity profiles using the empirical relationship proposed by Del Grosso (1974). To save the computational cost for ray tracing, the profile was approximated into a polygonal curve with several tens of nodes (Figure 5).

FIGURE 5

FIGURE 5. Reference sound speed profiles (blue lines) for epochs (A) TOS2.1301 (January 2013), (B) TOS2.1508 (August 2015), (C) MYGI.1302 (February 2013), and (D) MYGI.1508 (August 2015). Red lines indicate 1-m sound speed profiles obtained from the 1-m layered XBT/XCTD data.

During a GNSS-A survey, the vessel sails on a pre-determined track over the seafloor transponder array to collect geometrically balanced acoustic data (e.g., Figure 1). The along-track observation (called “subset”, hereafter) is repeated several times by reversing the sailing direction in order to reduce the bias due to the errors in the ATD offset.

During an observation cruise, it occasionally took more than a few weeks to collect sufficient acoustic data at a single site due to weather conditions or other operational restrictions. Even so, we compiled a single dataset per site per cruise for the static seafloor positioning in practice, because the positional changes should be too small to detect. We call the collection of a single GNSS-A dataset “observation epoch” or “epoch”, hereafter.

We set the parameters for the numbers of basis functions, $K_a$, $K_b$, and $K_c$, in Eq. 5, as $nmp0=nmp1=nmp2=15$ for both preprocess and main process. Knot intervals of B-spline basis functions were approximately 10–20 min. for most epochs.

Array Geometry Determination

In order to calculate the proper array geometry $\overline{{\mathit{X}}_j}$ for the rigid-array constraint, we first determined the positions of each transponder for all observations. Note that not all transponders are used in each observation, for example, because of additional installation of transponders for replacing transponders which were decommissioned due to battery outage. $\overline{{\mathit{X}}_j}$ and the positional difference of the array center for $n$th observation, ${\mathit{c}}^{\left(n\right)}$ were calculated by solving the following simultaneous equations:

$\begin{array}{c}\{\begin{array}{c}{\mathit{X}}_j^{\left(n\right)}={\delta }_j^{\left(n\right)}\overline{{\mathit{X}}_j}+{\delta }_j^{\left(n\right)}{\mathit{c}}^{\left(n\right)} \left(\text{for} j=1\dots J \text{and} n=1\dots N\right) \\ 0={\displaystyle \sum }_{n=1}^N{\mathit{c}}^{\left(n\right)}\end{array},\end{array}(24.1)$

with,

$\begin{array}{c}{\delta }_j^{\left(n\right)}=\{\begin{array}{c}1 \mathrm{i}\mathrm{f} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r} j \mathrm{i}\mathrm{s} \mathrm{u}\mathrm{s}\mathrm{e}\mathrm{d} \mathrm{i}\mathrm{n} n \mathrm{t}\mathrm{h} \mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \\ 0 \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s}\end{array},\end{array}(24.2)$

where $J$ and $N$ are the number of transponders and observations, respectively, and ${\mathit{X}}_j^{\left(n\right)}$ denotes the predetermined transponders’ positions for the $n$th observation.

The preliminary array-free positioning was also used for the verification of the collected data. We eliminated the outliers whose discrepancies from the preliminary solution were larger than the arbitrary threshold. We set the threshold to be 5 times as large as the root mean square value (RMS) of the travel time residuals.

Hyperparameter Search

In order to get the solution ${\mathit{x}}^{\mathit{*}}$, we should determine the appropriate values for the various hyperparameters, i.e., ${\sigma }^2$, ${\mu }_t$, ${\mu }_{MT}$, ${\rho }^2$, ${\lambda }_0^2$, ${\lambda }_{1E}^2$, ${\lambda }_{1N}^2$, ${\lambda }_{2E}^2$, and ${\lambda }_{2N}^2$. In the scheme of the ABIC minimization, ${\sigma }^2$ can be determined analytically by Eq. 23. It is reasonable to assume ${\lambda }_{1E}^2={\lambda }_{1N}^2={\lambda }_{2E}^2={\lambda }_{2N}^2$ because these hyperparameters control the smoothness of the spatial sound speed structure. We hereafter use a variable ${\lambda }_g^2$ to represent the value of these hyperparameters. For the purpose of single positioning, $\rho$ should be a large number, for example in meter-order. The large $\rho$ hardly changes the ABIC value and thus the solution.

In order to save the computational resources, we should further reduce the number of hyperparameters. We tentatively put ${\mu }_{MT}=0.5$. For the sound speed variations, we had to assume the strong constancy of spatial sound speed structure to resolve them with the single transducer GNSS-A. For this reason, we selected the ratio of ${\lambda }_0^2$ and ${\lambda }_g^2$, as ${\lambda }_g^2=0.1 {\lambda }_0^2$. The last two hyperparameters, ${\mu }_t$ and ${\lambda }_0^2$, were determined with the grid search method. The tested values for ${\mu }_t$ and ${\lambda }_0^2$ are ${\mu }_t=\left(0 \text{min}., 0.5 \text{min}., 1 \text{min}., 2 \text{min}.,3 \text{min}.\right)$ and ${\lambda }_0^2=\left({10}^{-3}, {10}^{-2}, {10}^{-1}, {10}^0, {10}^1, {10}^2\right)$, respectively.

Results

Figure 6 shows the time series of the estimated positions at sites TOS2 and MYGI. The positions are aligned to the ITRF2014 (Altamimi et al., 2016) and transformed into local ENU coordinates. Comparing the time series derived by the existing scheme (SGOBS version 4.0.2; used in Yokota et al., 2019), GARPOS reproduced almost the same trends for both sites. GARPOS might have succeeded in slightly suppressing the dispersion around the averaged velocity extracted from the neighboring epochs. Whereas the previous method corrected the sound speed gradient structure with step-by-step procedure, the proposed method has an advantage in directly extracting the structure by simultaneous estimation of all parameters.

FIGURE 6

FIGURE 6. Time series of displacement at (A) TOS2 and (B) MYGI solved by GARPOS (orange circles) and SGOBS version 4.0.2 (blue squares). The positions are aligned to the ITRF 2014.

TOS2 is located offshore in the south of Shikoku Island, southwestern Japan, above the source region of the 1946 Nankaido earthquake (e.g., Sagiya and Thatcher, 1999) along the Nankai Trough. According to Yokota and Ishikawa (2020), who investigated the transient deformations at the GNSS-A sites along the Nankai Trough, no significant signal was detected at TOS2. The results by the proposed method show the same trends as the conventional results. Although the trend of horizontal displacement seems to be changed in 2018 or 2019, careful inspection is needed because the transponders had been replaced during this period.

MYGI is located in the offshore east of Miyagi Prefecture, northeastern Japan, which experienced the 2011 Tohoku-oki earthquake (Sato et al., 2011). After the earthquake, significant westward postseismic movement and subsidence due to the viscoelastic relaxation has been observed at MYGI (Watanabe et al., 2014). The postseismic movements continue but appear to decay. It is true that the changes in the displacement rate at these sites are crucial in seismic and geodetic researches, but discussing these matters is beyond the scope of the present paper. The point is that the seafloor positioning results were well reproduced by the proposed method.

Discussions

Interpretations for the Correction Coefficient

As mentioned in Sound Speed Perturbation Model, it is convenient to relate the correction coefficient to the sound speed perturbation by assuming the case for $\left|{\gamma }_i\right|\ll 1$ for better understanding, though observation equation (Eq. 6) is valid for arbitrary value of ${\gamma }_i$. For the relationship $\delta V_i \sim {\gamma }_i\overline{V_0}$, we can convert each term of $\Gamma$ into the dimension of speed and speed gradient as, $\delta V_0\left(t\right)\equiv \overline{V_0}{\alpha }_0\left(t\right)$, ${\mathit{g}}_1\left(t\right)\equiv \overline{V_0}{\mathit{\alpha }}_1\left(t\right)$, and ${\mathit{g}}_2\left(t\right)\equiv \overline{V_0}{\mathit{\alpha }}_2\left(t\right)$.

The early models by Fujita et al. (2006) and Ikuta et al. (2008) took only the term $\delta V_0\left(t\right)$ into account. Whereas Ikuta et al. (2008) used the cubic B-spline functions as basis functions, Fujita et al. (2006) applied the multiple second degree polynomial functions with 10–20-min. time windows. Although these models do not include any transponder dependent term ${\mathit{g}}_2\left(t\right)$, the transponder independent spatial gradient ${\mathit{g}}_1\left(t\right)$ can be indirectly extracted as shown by Yokota et al. (2019).

In addition to estimating the term identical to $\delta V_0$, Yokota et al. (2019) implemented the additional process to estimate ${\mathit{g}}_2$ from the residuals of the solution by the method of Fujita et al. (2006). Strictly, the derived parameters in their scheme, i.e., $\Delta V_1$ and $\Delta V_2$ in Yokota et al. (2019), are the same as ${\mathit{g}}_1+{\mathit{g}}_2$ and ${\mathit{g}}_2$ in this study, respectively. For these parameters, our team have already made a qualitative interpretation in Yokota and Ishikawa (2019).

In order to show the relationship with other conventional models, we expand the proposed formulation to those by Yasuda et al. (2017), Honsho et al. (2019), and Kinugasa et al. (2020). Because Honsho et al. (2019) practically assumed 1-directional sound speed gradient, they constructed the model basically in the 2-dimensional plane spanned by the gradient direction and vertical direction.

For simplification, we assume that the ray path is a straight line connecting both ends. Putting $L^{\ast }$ equal to the depth of the observation site, the emission angle $\theta$ defined in Figure 3 of Honsho et al. (2019) can be expressed as,

$\begin{array}{c}\frac{X_j}{L^{\ast }}-\frac{P\left(t\right)}{L^{\ast }}=\tan \theta .\end{array}(25)$

Furthermore, assuming that the transmit/reception positions are the same and that the difference between transmit/reception time is so small that ${\alpha }_0\left(t\right)$, ${\alpha }_1\left(t\right)$, and ${\alpha }_2\left(t\right)$ hardly change, ${\gamma }_i$ can be written as,

$\begin{array}{c}{\gamma }_i={\alpha }_0\left(t\right)+\left({\alpha }_1\left(t\right)+{\alpha }_2\left(t\right)\right)\frac{P\left(t\right)}{L^{\ast }}+{\alpha }_2\left(t\right)\tan \theta .\end{array}(26)$

Because $\delta T$ defined in Eqs 2 and 5 of Honsho et al. (2019) is equivalent to $T_i^c-{\tau }_i$ in our formulation, we have,

$\begin{array}{c}\left(\exp \left(-{\gamma }_i\right)-1\right){\tau }_i=\frac{1}{\cos \theta }\left(c_0\left(t\right)+g\left(t\right)x_0+w\left(t\right)\tan \theta \right),\end{array}(27)$

where $c_0\left(t\right)$, $g\left(t\right)$, $w\left(t\right)$, and $x_0=P$are defined in Eqs. 6–8 of Honsho et al. (2019) and the transducer’s position in their formulation, respectively. Recalling that the slant range of acoustic ray path is $2L^{\ast }/\cos \theta$, the reference round-trip travel time can be written as,

$\begin{array}{c}{\tau }_i=\frac{2L^{\ast }}{\overline{V_0\left(u\right)}\cos \theta }.\end{array}(28)$

Considering the case where $\left|{\gamma }_i\right|\ll 1$, Eq. 27 is approximated to,

$\begin{array}{c}-\frac{2L^{\ast }}{\overline{V_0}}{\gamma }_i=c_0\left(t\right)+g\left(t\right)x_0+w\left(t\right)\tan \theta .\end{array}(29)$

From Eqs 26 and 29, the following relationships are derived:

$\begin{array}{c}{c}_0\left(t\right)=-\frac{2L^{\ast }}{\overline{V_0}}{\alpha }_0\left(t\right),\end{array}(30.1)$

$\begin{array}{c}g\left(t\right)=-\frac{2}{\overline{V_0}}\left({\alpha }_1\left(t\right)+{\alpha }_2\left(t\right)\right),\end{array}(30.2)$

$\begin{array}{c}w\left(t\right)=-\frac{2L^{\ast }}{\overline{V_0}}{\alpha }_2\left(t\right).\end{array}(30.3)$

In Honsho et al. (2019), $w\left(t\right)$ is extended to a 2-dimensional vector, i.e.,

$\begin{array}{c}\mathit{w}\left(t\right)=-\frac{2L^{\ast }}{\overline{V_0}}{\mathit{\alpha }}_2\left(t\right).\end{array}(31.1)$

Similarly, when extending $g\left(t\right)$ to the 2-dimensional vector, we can use the following vector form:

$\begin{array}{c}\mathit{g}\left(t\right)=-\frac{2}{\overline{V_0}}\left({\mathit{\alpha }}_1\left(t\right)+{\mathit{\alpha }}_2\left(t\right)\right),\end{array}(31.2)$

though they consequently use the assumption that $\mathit{g}\left(t\right)$ is parallel to $\mathit{w}\left(t\right)$. It is equivalent to the case that ${\mathit{\alpha }}_1$ is parallel to ${\mathit{\alpha }}_2$ in the proposed formulation.

Honsho et al. (2019) supposed the physical model where a spatially homogeneous 1-directional gradient of slowness lies in the certain layer, from sea-surface to the depth $D$, in the water. In such cases, $\mathit{w}\left(t\right)$ is proportional to $\mathit{g}\left(t\right)$, as $\mathit{w}=\left(D/2\right)\mathit{g}$. This is exactly the same assumption as the model by Yasuda et al. (2017). The model of Kinugasa et al. (2020) is the special case of those models where $D$ equals to the water depth.

In the proposed method, the sound speed field is approximately interpreted by their models when the unit vector of ${\mathit{\alpha }}_1$ is supposed to be same as that of ${\mathit{\alpha }}_2$ and $\left|{\mathit{\alpha }}_1\right|\ge \left|{\mathit{\alpha }}_2\right|$. The depth of the gradient layer is calculated as,

$\begin{array}{c}D=\frac{2L^{\ast }}{1+{\mathit{\alpha }}_1/{\mathit{\alpha }}_2}.\end{array}(32)$

When ${\mathit{\alpha }}_1={\mathit{\alpha }}_2$, it concludes to the model of Kinugasa et al. (2020). Conversely, when $\left|{\mathit{\alpha }}_2\right|\ll \left|{\mathit{\alpha }}_1\right|$, sound speed gradient lies in the thin layer near the surface.

In addition to the simple model above, the proposed method can extract more complicated sound speed field, which partly described by Yokota and Ishikawa (2019). Extracted parameters for the sound speed perturbation indicate the complicity of oceanographic structure, as shown in the next section.

Validity of Extracted Sound Speed Perturbation Model

Typical examples for the estimation results for each observation, i.e., the time series of travel time residuals, and sound speed perturbation interpreted from the correction coefficient, are shown in Figure 7. Results for all the datasets are available in Supplementary Figure S1.

FIGURE 7

FIGURE 7. Estimated results of the most preferred model for epochs (A) TOS2.1301.kaiyo_k4, (B) TOS2.1508.meiyo_m5, and (C) TOS2.1711.kaiyo_k4, (D) MYGI.1211.kaiyo_k4, (E) MYGI.1508. kaiyo_k4, and (F) MYGI.1802. kaiyo_k4. The top panels show the model residuals of the round-trip travel time. The second panels show the rejected acoustic data in the preprocessing step for determining the array geometry. The third panels indicate the sound speed perturbations, i.e.,${\gamma }_i\overline{V_0}$ (the crosses), and $\delta V_0\left(t\right)\equiv \overline{V_0}{\alpha }_0\left(t\right)$ (black line). The colors of the symbols in these panels identify the target transponders. The blue and purple arrows on the bottom panels indicate the spatial gradient of the sound speed perturbations in north-up expression, i.e., ${\mathit{g}}_1\left(t\right)\equiv \overline{V_0}{\mathit{\alpha }}_1\left(t\right)$, and ${\mathit{g}}_2\left(t\right)\equiv \overline{V_0}{\mathit{\alpha }}_2\left(t\right)$, respectively. Dotted lines and solid lines show the temporal variations of eastward and northward components, respectively. The colored horizontal lines denote the ranges of the observation subsets.

In the most cases for site TOS2, both terms of the estimated sound speed gradient vector stably direct south to southeast. Because the sound speed increases with the water temperature, it means that the water temperature is higher in the southern region. The results that ${\mathit{g}}_2$ is comparable with ${\mathit{g}}_1$ in many cases indicate that the gradient of water temperature continues to the deeper portion, as discussed in the previous section. This is consistent with the fact that the Kuroshio current continuously flows on the south of TOS2.

In contrast, the directions of gradient terms at MYGI have less constancy than TOS2. Unlike the area around TOS2 where the Kuroshio current dominantly affects the seawater structure, MYGI is located in an area with a complicated ocean current system (e.g., Yasuda, 2003; Miyazawa et al., 2009). Watanabe and Uchida (2016) have also shown that the temperature profiles at MYGI vary widely with observation epochs. These features cannot be resolved by the simpler model with single sound speed gradient parameter.

The complexity in the sound speed variation at MYGI tends to lead to large variations in the residual travel time. Nevertheless, the proposed method successfully extracted the smooth sound speed structure for many observation epochs, except a few epochs such as June 2013 (MYGI.1306.kaiyo_k4) and June 2019 (MYGI.1906.meiyo_m5) shown in Supplementary Figure S1. In these epochs, relatively larger values for the hyperparameter ${\lambda }_0^2$ were adopted and caused larger variations in each term of $\Gamma$. Possible causes of this include the systematic errors in other observation subcomponents such as the random walk noise in GNSS positioning, the drifts of gyro sensor, or the time synchronization error between the devices.

Preferred models for all the tested epochs had positive values for data correlation length, ${\mu }_t$. It is considered that the plausible estimation of sound speed is realized by introducing the statistic information criteria and the information of data covariances.

In order to discuss the effects of the data covariance, we tested the cases for the models without assuming the data correlation, i.e., ${\mu }_t=0$. Figure 8 shows the preferred models selected from ${\lambda }_0^2=\left({10}^{-3}, {10}^{-2}, {10}^{-1}, {10}^0, {10}^1, {10}^2, {10}^3, {10}^4\right)$ and ${\mu }_t=0$. It is clear that the preferred models without assuming the data correlation have larger ${\lambda }_0^2$. Although the residuals of travel time were reduced in these models, overfittings occurred for each term of $\Gamma$. Comparing the preferred and less-preferred results, the existence of data covariance components contributes to suppressing the overfitting and to selecting a model with less perturbation by decreasing the impact of individual data on model parameters.

FIGURE 8

FIGURE 8. Same as Figure 7, but for the most preferred model in the models with ${\mu }_t=0$.

To confirm the stability of the seafloor positioning results, the differences of seafloor position for the tested models from the most preferred models are summarized in Figure 9. The differences in estimated positions for most of the tested models converged in several centimeters. For both sites, variations in the vertical component tend to be larger for larger values of ${\lambda }_0^2$. It indicates that finer hyperparameter tuning is not required when considering the application to seafloor positioning.

FIGURE 9

FIGURE 9. Distributions of differences of positions of the tested models from the preferred ones at (A) TOS2 and (B) MYGI for northward-eastward (left), northward-upward (center), and upward-eastward (right) components. The colors of circles indicate the value of ${\lambda }_0^2/{\lambda }_0^{2^{\ast }}$.

As another application of GNSS-A to oceanography, temporal changes of the oceanographic structure within the observation epoch can be extracted using the proposed method. For example, the estimated sound speed gradient ${\mathit{g}}_1$ in the epoch of MYGI.1802.kaiyo_k4 (Figure 7F) suggests that the dominant oceanographic structure had changed at 01:00–03:00 UTC. On the other hand, a temporal variation with a relatively short period of several tens of minutes remains in the travel time residuals, which might be caused by the internal gravity wave. To improve the detectability of relatively short-period perturbations, further adjustments and verifications of the proposed model will be required.

Conclusion

We reconstructed the GNSS-A observation equation and developed the Python-based software GARPOS to solve the seafloor position as well as the sound speed perturbations using the empirical Bayes approach. It provides a stable solution for a generally ill-posed problem caused by the correlation among the model parameters, by introducing the hyperparameter tuning based on the ABIC minimization and data covariance to rationalize the normalization constant of the posterior pdf.

The most important point is that the proposed method succeeded in directly extracting the time-dependent sound speed field with two end members of spatial gradient terms, which are roughly characterized by depths, even when the observers used only one sea-surface unit. Statistical approach allowed us to suppress the overfitting and thus to obtain simpler sound speed field from densely collected dataset. It successfully reproduced the stationary southward sound speed gradient at TOS2, which is consistent with the Kuroshio current.

On the other hand, model overfits were shown in several epochs. These overfits can be caused not only by the actually complicated sound speed field but also by other error sources which were not well included in the model. It means that the hyperparameter tuning also plays a role in the verification of dataset and model. Error analyses in such cases might rather help improving the GNSS-A accuracy and methodology.

We suggested a simplified formatting for the GARPOS input files. Researchers can enter into the field of seafloor geodesy by collecting the listed data with adequate precision. Since each subcomponent of GNSS-A technique, i.e., GNSS positioning, acoustic ranging, and so on, has been well established, observers can combine them on their platform. Especially, GNSS-A is expected to be practicalized in the near future with an unmanned surface vehicle (Chadwell, 2016) or a buoy (e.g., Kinugasa et al., 2020; Tadokoro et al., 2020). Even in the case of the stationary observation due to small cruising speed, GARPOS may provide the solutions by making a slight modification in the prior variance-covariance matrix.

There is a room for improvement in setting the prior information for transponders’ positions, ${\mathit{X}}_j^0$. For instance, the displacement of transponder array from the previous epoch is predicted as small as several centimeters when the interval of observation visits is short. Such assumption leads to the application of the inter-epoch filtering. Furthermore, it has a possibility to progress to the kinematic seafloor positioning, as shown by Tomita et al. (2019). We expect that the publication of GARPOS on the open-access repository will enhance the researchers’ engagement and the future development on the GNSS-A technique.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GARPOS v0.1.0 (https://zenodo.org/record/3992688); GNSS‐A data obtained at the sites “TOS2” and “MYGI” in 2011-2019 (https://doi.org/10.5281/zenodo.3993912).

Author Contributions

SW designed the study and wrote the manuscript. SW developed “GARPOS” and processed the data. SW, TI, YY, and YN discussed about the methodology and commented to improving the manuscript.

Funding

The submission of the manuscript was funded by the Japan Coast Guard.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank many staff members from the Hydrographic and Oceanographic Department, Japan Coast Guard, including the crew of the survey vessels Takuyo, Shoyo, Meiyo, and Kaiyo for their support in our observations and technological developments. We especially thank the active senior staff members from the Geodesy and Geophysics Office, Hydrographic and Oceanographic Department, Japan Coast Guard, for their devoted maintenance and management of the equipment. We also thank the reviewers for their comments and suggestions for improving the manuscript.

Supplementary Materials

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.597532/full#supplementary-material

References

Akaike, H. (1980). “Likelihood and the Bayes procedure”, In Bayesian statistics. Editors J. M. Bernardo, et al. , (Valencia, Spain: University Press), 143–166.

Google Scholar

Altamimi, Z., Rebischung, P., Métivier, L., and Collilieux, X. (2016). ITRF2014: a new release of the International Terrestrial Reference Frame modeling nonlinear station motions. J. Geophys. Res. Solid Earth 121, 6109. doi:10.1002/2016JB013098.

CrossRef Full Text | Google Scholar

Asada, A., and Yabuki, T. (2001). Centimeter-level positioning on the seafloor. Proc. Jpn. Acad. Ser. B Phys. Biol. Sci. 77, 7–12. doi:10.2183/pjab.77.7.

CrossRef Full Text | Google Scholar

Bürgmann, R., and Chadwell, D. (2014). Seafloor geodesy. Annu. Rev. Earth Planet Sci. 42(1), 509–534. doi:10.1146/annurev-earth-060313-054953.

CrossRef Full Text | Google Scholar

Chadwell, C. D. (2016). “Campaign-style GPS-Acoustic with wave gliders and permanent seafloor benchmarks”, in Proceedings of the subduction zone observatory workshop. (Boise, ID:Boise Center), Sep. 29 – Oct. 1 2016.

Google Scholar

Chadwell, C. D., and Spiess, F. N. (2008). Plate motion at the ridge-transform boundary of the south Cleft segment of the Juan de Fuca Ridge from GPS-Acoustic data. J. Geophys. Res. 113, B04415. doi:10.1029/2007JB004936.

CrossRef Full Text | Google Scholar

Chadwell, C. D., and Sweeney, A. D. (2010). Acoustic ray-trace equations for seafloor geodesy. Mar. Geodes. 33(2–3), 164–186. doi:10.1080/01490419.2010.492283.

CrossRef Full Text | Google Scholar

Chen, H.-Y., Ikuta, R., Lin, C.-H., Hsu, Y.-J., Kohmi, T., Wang, C.-C., et al. (2018). Back-arc opening in the western end of the Okinawa Trough revealed from GNSS/Acoustic Measurements. Geophys. Res. Lett. 45, 137–145, doi:10.1002/2017GL075724.

CrossRef Full Text | Google Scholar

de Boor, C. (1978). A practical guide to splines. New York:Springer-Verlag, Vol. 27.

Google Scholar

Del Grosso, V. A. (1974). New equation for the speed of sound in natural waters (with comparisons to other equations). J. Acoust. Soc. Am. 56, 1084–1091. doi:10.1121/1.1903388.

CrossRef Full Text | Google Scholar

Fujimoto, H. (2014). Seafloor geodetic approaches to subduction thrust earthquakes, Monogr. Environ. Earth Planets 2, 23–63. doi:10.5047/meep.2014.00202.0023.

CrossRef Full Text | Google Scholar

Fujita, M., Ishikawa, T., Mochizuki, M., Sato, M., Toyama, S., Katayama, M., et al. (2006). GPS/acoustic seafloor geodetic observation: method of data analysis and its application. Earth Planets Space 58, 265–275. doi:10.1186/BF03351923.

CrossRef Full Text | Google Scholar

Fukahata, Y., and Wright, T. J. (2008). A non-linear geodetic data inversion using ABIC for slip distribution on a fault with an unknown dip angle. Geophys. J. Int. 173, 353–364. doi:10.1111/j.1365-246X.2007.03713.x.

CrossRef Full Text | Google Scholar

Gagnon, K., Chadwell, C. D., and Norabuena, E. (2005). Measuring the onset of locking in the Peru-Chile trench with GPS and acoustic measurements. Nature 434(7030), 205–208. doi:10.1038/nature03412.

CrossRef Full Text | Google Scholar

Honsho, C., and Kido, M. (2017). Comprehensive analysis of traveltime data collected through GPS-acoustic observation of seafloor crustal movements. J. Geophys. Res. Solid Earth 122, 8583–8599. doi:10.1002/2017JB014733.

CrossRef Full Text | Google Scholar

Honsho, C., Kido, M., Tomita, F., and Uchida, N. (2019). Offshore postseismic deformation of the 2011 Tohoku earthquake revisited: application of an improved GPS‐acoustic positioning method considering horizontal gradient of sound speed structure. J. Geophys. Res. Solid Earth 124, 5990. doi:10.1029/2018JB017135.

CrossRef Full Text | Google Scholar

Hovem, J. M. (2013). “Ray trace modeling of underwater sound propagation”, in Modeling and measurement methods for acoustic waves and for acoustic microdevices. M. G. Editor Beghi. (Rijeka, Croatia: IntechOpen). doi:10.5772/55935.

CrossRef Full Text | Google Scholar

IERS Conventions (2010). Editors. G. Petit and B. Luzum. (IERS technical note; 36) frankfurt am main: Verlag des Bundesamts fr Kartographie und Geodsie, 179, ISBN 3-89888-989-6

Google Scholar

Ikuta, R., Tadokoro, K., Ando, M., Okuda, T., Sugimoto, S., Takatani, K., et al. (2008). A new GPS-acoustic method for measuring ocean floor crustal deformation: application to the Nankai Trough. J. Geophys. Res. 113, B02401. doi:10.1029/2006JB004875.

CrossRef Full Text | Google Scholar

International GNSS Service (a). GNSS final combined orbit solution product, Greenbelt, MD, USA: NASA Crustal Dynamics Data Information System (CDDIS). doi:10.5067/gnss/gnss_igsorb_001.

CrossRef Full Text | Google Scholar

International GNSS Service (b). GNSS final combined satellite and receiver clock solution (30 second) product, Greenbelt, MD, USA: NASA Crustal Dynamics Data Information System (CDDIS). doi:10.5067/GNSS/gnss_igsclk30_001.

CrossRef Full Text | Google Scholar

Ishikawa, T., Yokota, Y., Watanabe, S., and Nakamura, Y. (2020). History of on-board equipment improvement for GNSS-A observation with focus on observation frequency, Front. Earth Sci. 8, 150. doi:10.3389/feart.2020.00150.

CrossRef Full Text | Google Scholar

Jensen, F. B., Kuperman, W. A., Porter, M. B., and Schmidt, H. (2011). Computational ocean acoustics, Vol. 97. New York, NY: Springer New York, ISBN:978-1-4419-8677-1.

Google Scholar

Kido, M., Osada, Y., Fujimoto, H., Hino, R., and Ito, Y. (2011). Trench-normal variation in observed seafloor displacements associated with the 2011 Tohoku-Oki earthquake. Geophys. Res. Lett. 38, L24303. doi:10.1029/2011GL050057.

CrossRef Full Text | Google Scholar

Kido, M., Osada, Y., and Fujimoto, H. (2008). Temporal variation of sound speed in ocean: a comparison between GPS/acoustic and in situ measurements. Earth Planets Space 60(3), 229–234, doi:10.1186/BF03352785.

CrossRef Full Text | Google Scholar

Kinugasa, N., Tadokoro, K., Kato, T., and Terada, Y. (2020). Estimation of temporal and spatial variation of sound speed in ocean from GNSS-A measurements for observation using moored buoy. Prog. Earth Planet. Sci. 7, 21. doi:10.1186/s40645-020-00331-5.

CrossRef Full Text | Google Scholar

Matsu'ura, M., Noda, A., and Fukahata, Y. (2007). Geodetic data inversion based on Bayesian formulation with direct and indirect prior information Geophys. J. Int. 171(3), 1342–1351. doi:10.1111/j.1365-246X.2007.03578.x.

CrossRef Full Text | Google Scholar

Matsumoto, Y., Fujita, M., and Ishikawa, T. (2008). Development of multi-epoch method for determining seafloor station position [in Japanese], Rep. Hydrogr. Oceanogr. Res. 26, 16–22| .

Google Scholar

Miyazawa, Y., Zhang, R., Guo, X., Tamura, H., Ambe, D., Lee, J.-S., et al. (2009). Water mass variability in the western North Pacific detected in a 15-year eddy resolving ocean reanalysis. J. Oceanogr. 65, 737–756. doi:10.1007/s10872-009-0063-3.

CrossRef Full Text | Google Scholar

Nagaya, Y. (1995). Basic study on a sea floor strain measurement using acoustic techniques [in Japanese with English abstracts]. Rep. Hydrogr. Res. 31, 67–76.

Google Scholar

Obana, K., Katao, H., and Ando, M. (2000). Seafloor positioning system with GPS-acoustic link for crustal dynamics observation-a preliminary result from experiments in the sea-. Earth Planets Space. 52, 415–423. doi:10.1186/BF03352253.

CrossRef Full Text | Google Scholar

Osada, Y., Fujimoto, H., Miura, S., Sweeney, A., Kanazawa, T., Nakao, S., et al. (2003). Estimation and correction for the effect of sound velocity variation on GPS/Acoustic seafloor positioning: an experiment off Hawaii Island. Earth Planets Space 55, e17–e20. doi:10.1186/BF03352464.

CrossRef Full Text | Google Scholar

Sagiya, T., and Thatcher, W. (1999). Coseismic slip resolution along a plate boundary megathrust: the Nankai Trough, southwest Japan. J. Geophys. Res. 104(B1), 1111–1129. doi:10.1029/98JB02644.

CrossRef Full Text | Google Scholar

Sakic, P., Ballu, V., Crawford, W. C., and Wöppelmann, G. (2018). Acoustic ray tracing comparisons in the context of geodetic precise off-shore positioning experiments. Mar. Geodes. 41(4), 315–330. doi:10.1080/01490419.2018.1438322.

CrossRef Full Text | Google Scholar

Sato, M., Fujita, M., Matsumoto, Y., Saito, H., Ishikawa, T., and Asakura, T. (2013). Improvement of GPS/acoustic seafloor positioning precision through controlling the ship's track line. J. Geodyn. 87, 825–842. doi:10.1007/s00190-013-0649-9.

CrossRef Full Text | Google Scholar

Sato, M., Ishikawa, T., Ujihara, N., Yoshida, S., Fujita, M., Mochizuki, M., et al. (2011). Displacement above the hypocenter of the 2011 Tohoku-oki earthquake. Science 332, 1395. doi:10.1126/science.1207401.

CrossRef Full Text | Google Scholar

Spiess, F. N. (1980). Acoustic techniques for marine geodesy, Mar. Geodes. 4(1), 13–27. doi:10.1080/15210608009379369.

CrossRef Full Text | Google Scholar

Spiess, F. N., Chadwell, C. D., Hildebrand, J. A., Young, L. E., Purcell, G. H., and Dragert, H. (1998). Precise GPS/Acoustic positioning of seafloor reference points for tectonic studies. Phys. Earth Planet. In. 108(2), 101–112. doi:10.1016/s0031-9201(98)00089-2.

CrossRef Full Text | Google Scholar

Tadokoro, K., Kinugasa, N., Kato, T., Terada, Y., and Matsuhiro, K. (2020). A marine-buoy-mounted system for continuous and real-time measurment of seafloor crustal deformation. Front. Earth Sci. 8, 123. doi:10.3389/feart.2020.00123.

CrossRef Full Text | Google Scholar

Takasu, T. (2013). RTKLIB ver. 2.4.2: an open source program package for GNSS positioning, http://www.rtklib.com/

Google Scholar

Tarantola, A., and Valette, B. (1982). Generalized nonlinear inverse problems solved using the least squares criterion. Rev. Geophys. 20(2), 219–232. doi:10.1029/RG020i002p00219.

CrossRef Full Text | Google Scholar

Tomita, F., Kido, M., Honsho, C., and Matsui, R. (2019). Development of a kinematic GNSS-Acoustic positioning method based on a state-space model. Earth Planets Space 71, 102. doi:10.1186/s40623-019-1082-y.

CrossRef Full Text | Google Scholar

Tomita, F., Kido, M., Ohta, Y., Iinuma, T., and Hino, R. (2017). Along-trench variation in seafloor displacements after the 2011 Tohoku earthquake. Sci. Adv. 3(7), e1700113. doi:10.1126/sciadv.1700113.

CrossRef Full Text | Google Scholar

Tomita, F., Kido, M., Osada, Y., Hino, R., Ohta, Y., and Iinuma, T. (2015). First measurement of the displacement rate of the Pacific Plate near the Japan Trench after the 2011 Tohoku-Oki earthquake using GPS/acoustic technique. Geophys. Res. Lett. 42, 8391–8397, doi:10.1002/2015GL065746.

CrossRef Full Text | Google Scholar

Watanabe, S., Sato, M., Fujita, M., Ishikawa, T., Yokota, Y., Ujihara, N., et al. (2014). Evidence of viscoelastic deformation following the 2011 Tohoku-oki earthquake revealed from seafloor geodetic observation. Geophys. Res. Lett. 41, 5789–5796. doi:10.1002/2014GL061134.

CrossRef Full Text | Google Scholar

Watanabe, S., Ishikawa, T., Yokota, Y., and Nakamura, Y. (2020a). GARPOS v0.1.0: analysis tool for GNSS-Acoustic seafloor positioning. Zenodo. Version 0.1.0. doi:10.5281/zenodo.3992688.

CrossRef Full Text | Google Scholar

Watanabe, S., Ishikawa, T., Yokota, Y., and Nakamura, Y. (2020b). GNSS-A data obtained at the sites “TOS2” and “MYGI” in 2011-2019. Zenodo. doi:10.5281/zenodo.3993912.

CrossRef Full Text | Google Scholar

Watanabe, S., and Uchida, T. (2016). Stable structures of temperature and salinity validated by the repeated measurements in the few-miles square regions off Japan coast in the western Pacific [in Japanese with English abstract]. Rep. Hydro. Ocean. Res. 53, 57–81.

Google Scholar

Watanabe, S., Yokota, Y., and Ishikawa, T. (2020c). Stability test to validate the GNSS-A seafloor positioning with kinematic precise point positioning [in Japanese with English abstract and captions]. J. Geod. Soc. Japan 66, 1–7. doi:10.11366/sokuchi.66.1.

CrossRef Full Text | Google Scholar

Yasuda, I. (2003). Hydrographic structure and variability in the kuroshio-oyashio transition area. J. Oceanogr. 59, 389–402. doi:10.1023/A:1025580313836.

CrossRef Full Text | Google Scholar

Yasuda, K., Tadokoro, K., Taniguchi, S., Kimura, H., and Matsuhiro, K. (2017). Interplate locking condition derived from seafloor geodetic observation in the shallowest subduction segment at the Central Nankai Trough, Japan. Geophys. Res. Lett. 44, 3572–35790. doi:10.1002/2017GL072918.

CrossRef Full Text | Google Scholar

Yokota, Y., and Ishikawa, T. (2019). Gradient field of undersea sound speed structure extracted from the GNSS-A oceanography: GNSS-A as a sensor for detecting sound speed gradient. SN Appl. Sci. 1, 693. doi:10.1007/s42452-019-0699-6.

CrossRef Full Text | Google Scholar

Yokota, Y., and Ishikawa, T. (2020). Shallow slow slip events along the Nankai Trough detected by GNSS-A, Sci. Adv. 6(3), eaay5786. doi:10.1126/sciadv.aay5786.

CrossRef Full Text | Google Scholar

Yokota, Y., Ishikawa, T., and Watanabe, S. (2019). Gradient field of undersea sound speed structure extracted from the GNSS-A oceanography. Mar. Geophys. Res. 40(4), 493–504. doi:10.1007/s11001-018-9362-7.

CrossRef Full Text | Google Scholar

Yokota, Y., Ishikawa, T., and Watanabe, S. (2018). Seafloor crustal deformation data along the subduction zones around Japan obtained by GNSS-A observations. Sci. Data 5, 180182. doi:10.1038/sdata.2018.182.

CrossRef Full Text | Google Scholar

Yokota, Y., Ishikawa, T., Watanabe, S., Tashiro, T., and Asada, A. (2016). Seafloor geodetic constraints on interplate coupling of the Nankai Trough megathrust zone. Nature 534, 374–377. doi:10.1038/nature17632.

CrossRef Full Text | Google Scholar