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.
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,  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,  extracted the transponder-dependent correction term from the residuals of the results derived by the conventional method of Fujita et al. (2006). Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532 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  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 . 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), 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).

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 earthcentered, earth-fixed (ECEF) coordinates such as the International Terrestrial Reference Frame (ITRF). When the transducer's position, P(t), 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 severalkilometers-width, solid-earth-tide-free solutions can be obtained when the trajectory of the transducer is determined in the solidearth-tide-free coordinates. Hereafter, the positions are expressed in solid-earth-tide-free coordinates in this paper.
In order to determine P(t) 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, Q(t), can be determined using any of appropriate kinematic GNSS solvers. The gyro sensor provides the attitude of the sea-surface platform, Θ(t) θ r θ p θ h 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 Q(t) from ECEF to local ENU coordinates, i.e., Q(t) Q e Q n Q u T . Using the relative position of the transducer to the GNSS antenna in the gyro's rectangular coordinate (called "ATD offset" hereafter; Figure 2), M M r M p M h T , we obtain the transducer's position in the local ENU coordinates as, The ATD offset values should be measured before the GNSS-A observation. Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532

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 ith acoustic signal to the jth 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(e, n, u, t), i.e., where t i+ , t i− , and X j are the transmitted and received time for the ith 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: where τ i and V 0 (u) denote the reference travel time and the reference sound speed profile, respectively. V 0 (u) 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(−c i ), 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, c i and X j are assigned as the estimator. Eq. 1 gives the transducer's position P(t) as a function of the GNSS antenna's position Q(t), the attitude vector Θ(t), and the ATD offset M. The time-independent parameter 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 τ i τ i (X j , M|Q(t), Θ(t), V 0 (u)), 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, c i ≪ 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 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., Γ(t, P, X), 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, c i is approximated by the average of the transmission and the reception times, i.e., c i ≡ 1 The function Γ(t, P, X) 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: where L p indicates the characteristic length of the observation site (typically in several kilometers). α 0 (t), α 1 (t) and α 2 (t) are the time-dependent coefficients for each term. Because the vertical variation of P and X are much smaller than the horizontal variation, we can practically ignore the vertical component of α 1 (t) and α 2 (t). Thus, α 1 (t) and α 2 (t) are reduced to a 2dimensional vector to denote the horizontal gradient. Each coefficient can be represented by a linear combination of basis functions Φ k (t): where a 〈·〉 k are the coefficients of the kth basis function, Φ 〈·〉 k (t), for each term named 〈 · 〉. E and N in 〈 · 〉 denote the eastward and northward components of the vector, respectively. For simplification, we compile these coefficients into vector a, hereafter.
Because the initial values for M and X j are usually obtained in the precision of less than meters prior to the GNSS-A analysis, we approximate P and X j in Γ substituting the initial values, i.e., M 0 and X 0 j , and not updating them with the iteration. This reduces the number of estimation parameters in the correction term, i.e., c i c i (a|X 0 j , M 0 , Q(t), Θ(t)).

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 X j X j + ΔX, where X j and ΔX denote the relative positions of each transponder and the parallel translation of the transponder array, respectively. The array geometry, X j , should be determined prior to the analytical procedure, using the data of multiple observation visits.
Meanwhile, 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 o i , are compared with the model, T c i . 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 ith acoustic round-trip travel time: (6.1) or in the form with the rigid-array constraint, where T p 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. 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, FIGURE 3 | Flow chart to construct the observation equation.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532 Let us consider the direct prior information for the model parameters X j and M written as, T denote the predicted model parameter vectors and the error vector, respectively. Let us assume that d X and d M follow a normal distribution with a variance-covariance of D X (ρ 2 ) and D M (ρ 2 ), whose scale can be adjusted by a hyperparameter ρ 2 , i.e., D X ρ 2 D X and D M ρ 2 D M , respectively. The prior probability density function (pdf) for the constraints can be written as, where c denotes the normalization constant.
For the model parameter a, an indirect prior information can be applied that the temporal change of sound speed perturbation model Γ 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 ϕ, whereas Ikuta et al. (2008) used the first derivative. When using the B-spline functions Φ 〈·〉 k (t) (e.g., de Boor, 1978) as the basis of temporal sound speed variation, the roughness can be written in a vector form, i.e., where, Then, the prior pdf can be written using the hyperparameter λ 〈·〉 as, where c denotes the normalization constant. Combining these prior informations, we obtain the following prior pdf: (12.2) where g and ||Λ G || 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 e also follows a normal distribution with a variance-covariance of σ 2 E, i.e., 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 spatiotemporal 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., μ t and μ MT , to adjust the non-diagonal component of E: if the transponders for i and j are the same 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 μ MT ∈ [0, 1] 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 o i − T c i follows a normal distribution. Nonetheless, because the GNSS-A typically gives the precision of T o i /T c i ∼ 1 ± 10 − 4 , we can suppose that T o i /T c i 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 (T p /T o i ) 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, whereG σ 2 G(ρ 2 , λ 2 ) and Λ G represents the absolute value of the product of non-zero eigenvalues ofG. Defining x(σ 2 , μ t , μ MT , ρ 2 , λ 2 ) as x that maximizes the posterior probability (Eq. 15) under the given hyperparameters, the partial derivative of p(x|y) with respect to x should be zero for x x. Hence, x should satisfy the following equation: where A(x) is the Jacobian matrix at point x defined as, We can solve the nonlinear equation (Eq. 16) numerically by performing an iterative method, where x k is corrected in each step with the following algorithm: to satisfy the following convergence criteria: can be rewritten as, Therefore, the linearized variance-covariance matrix around x can be obtained as,

Hyperparameter Tuning
The appropriate values of the hyperparameters can be determined by minimizing Akaike's Bayesian Information Criteria (ABIC; Akaike, 1980), where N HP denotes the number of hyperparameters. Although it is difficult to analytically calculate the integral for the marginal likelihood because of the nonlinearity in f (x), the Laplace's method can be applied in this case where the degree of freedom is sufficiently large and s(x) can be almost unimodal. Thus, an approximated form for ABIC is obtained as follows: where m is the number of model parameters. For the derivation, we used the following relationship: which is derived from the condition that the partial derivative of ABIC with respect to σ 2 should be zero. We can tune the hyperparameters to minimize the approximated ABIC value defined in Eq. 22, to obtain the solution x * x(σ 2 p , μ p t , μ p MT , ρ 2 p , λ 2 p ), where p 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 Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532 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: Major input/output parmeters and hyperparameters for GARPOS are listed in Tables 1 and 2, respectively.
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 Attitude of platform at t i+ roll0 I-4 Setting - Scale of a priori positioning error N/A N/A m 2 *Note that σ 2 is calculated analytically, and that ρ 2 is set in (I-2).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532 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.

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). 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 Q(t) 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, M, were measured by leveling, distance, and angle surveys before the first GNSS-A observation cruise, to be used as 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 (u), 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).
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 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. X j and the positional difference of the array center for nth observation, c (n) were calculated by solving the following simultaneous equations: where J and N are the number of transponders and observations, respectively, and X (n) j denotes the predetermined transponders' positions for the nth 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 x * , we should determine the appropriate values for the various hyperparameters, i.e., σ 2 , μ t , μ MT , ρ 2 , λ 2 0 , λ 2 1E , λ 2 1N , λ 2 2E , and λ 2 2N . In the scheme of the ABIC minimization, σ 2 can be determined analytically by Eq. 23. It is reasonable to assume λ 2 1E λ 2 1N λ 2 2E λ 2 2N because these hyperparameters control the smoothness of the spatial sound speed structure. We hereafter use a variable λ 2 g to represent the value of these hyperparameters. For the purpose of single positioning, ρ should be a large number, for example in meter-order. The large ρ 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 μ 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 λ 2 0 and λ 2 g , as λ 2 g 0.1 λ 2 0 . The last two hyperparameters, μ t and λ 2 0 , were determined with the grid search method. The tested values for μ t and λ 2 0 are μ t (0 min., 0.5 min., 1 min., 2 min., 3 min.) and λ 2 0 (10 − 3 , 10 − 2 , 10 − 1 , 10 0 , 10 1 , 10 2 ), respectively. 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 , 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.

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

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 c i ≪ 1 for better understanding, though observation equation (Eq. 6) is valid for arbitrary value of c i . For the relationship δV i ∼ c i V 0 , we can convert each term of Γ into the dimension of speed and speed gradient as, δV 0 (t) ≡ V 0 α 0 (t), g 1 (t) ≡ V 0 α 1 (t), and g 2 (t) ≡ V 0 α 2 (t).
The early models by Fujita et al. (2006) and Ikuta et al. (2008) took only the term δV 0 (t) 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 g 2 (t), the transponder independent spatial gradient g 1 (t) can be indirectly extracted as shown by .
In addition to estimating the term identical to δV 0 ,  implemented the additional process to estimate 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., ΔV 1 and ΔV 2 in , are the same as g 1 + g 2 and g 2 in this study, respectively. For these parameters, our team have already made a qualitative interpretation in .
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 p equal to the depth of the observation site, the emission angle θ defined in Figure 3 of Honsho et al. (2019) can be expressed as, Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532 Furthermore, assuming that the transmit/reception positions are the same and that the difference between transmit/reception time is so small that α 0 (t), α 1 (t), and α 2 (t) hardly change, c i can be written as, Because δT defined in Eqs 2 and 5 of Honsho et al. (2019) is equivalent to T c i − τ i in our formulation, we have, where c 0 (t), g(t), w(t), 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 p /cos θ, the reference round-trip travel time can be written as, Considering the case where c i ≪ 1, Eq. 27 is approximated to, From Eqs 26 and 29, the following relationships are derived: In Honsho et al. (2019), w(t) is extended to a 2-dimensional vector, i.e., Similarly, when extending g(t) to the 2-dimensional vector, we can use the following vector form: though they consequently use the assumption that g(t) is parallel to w(t). It is equivalent to the case that α 1 is parallel to α 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, w(t) is proportional to g(t), as w (D/2)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 α 1 is supposed to be same as that of α 2 and |α 1 | ≥ |α 2 |. The depth of the gradient layer is calculated as, When α 1 α 2 , it concludes to the model of Kinugasa et al. (2020). Conversely, when |α 2 | ≪ |α 1 |, 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 . 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.
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 g 2 is comparable with 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 λ 2 0 were adopted and caused larger variations in each term of Γ. 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, μ 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., μ t 0. Figure 8 shows the preferred models selected from λ 2 0 (10 − 3 , 10 − 2 , 10 − 1 , 10 0 , 10 1 , 10 2 , 10 3 , 10 4 ) and μ t 0. It is clear that the preferred models without assuming the data correlation have larger λ 2 0 . Although the residuals of travel time were reduced in these models, overfittings occurred for each term of Γ. 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.
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 λ 2 0 . It indicates that finer hyperparameter tuning is not required when considering the application to seafloor positioning.
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 g 1 in the epoch of MYGI.1802.kaiyo_k4 ( Figure 7F) suggests that the dominant show the rejected acoustic data in the preprocessing step for determining the array geometry. The third panels indicate the sound speed perturbations, i.e., c i V 0 (the crosses), and δV 0 (t) ≡ V 0 α 0 (t) (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., g 1 (t) ≡ V 0 α 1 (t), and g 2 (t) ≡ V 0 α 2 (t), 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.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 597532 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 FIGURE 8 | Same as Figure 7, but for the most preferred model in the models with μ t 0.