GARPOS: analysis software for the GNSS-A seafloor positioning with simultaneous estimation of sound speed structure

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


Basic configurations of the GNSS-A observation 24
Precise measurements of seafloor position in the global reference frame opens the door to the 25 "global" geodesy in the true sense of the word. It extended the observation network for crustal 26 deformation into the ocean and has revealed the tectonic processes in the subduction zone including 27 megathrust earthquakes (e.g., Bürgmann and Chadwell, 2014;Fujimoto, 2014, for review). Many 28 findings have been reported especially in the northwestern Pacific along the Nankai Trough (e.g., 29 Yokota combined) seafloor positioning technique, proposed by Spiess (1980). 33 Observers can take various ways to design the GNSS-A observation for the positioning of the 34 seafloor benchmark. They have to solve the difficulties not only in the technical realizations of 35 GNSS-A subcomponents such as the acoustic ranging and the kinematic GNSS positioning, but also 36 in designing the observation configurations and analytical models to resolve the strongly correlated 37 parameters. For example, because the acoustic ranging observations are performed only on the sea 38 surface, the errors in sound speed perturbations are strongly correlated with the relative distance, 39 typically the depths of the benchmark. 40 In the very first attempt for the realization, Spiess et al. (1998) derived horizontal displacement using 41 a stationary sea-surface unit which was approximately placed on the horizontal center of the array of 42 multiple seafloor mirror transponders. They determined the relative positions and depths of the 43 transponders in advance. The relative horizontal positions of the sea-surface unit to the transponder 44 array can be determined by acoustic ranging data, to be compared with the global positions 45 determined by space geodetic technique. In this "stationary" GNSS-A configuration, the temporal 46 variation of sound speed is less likely to affect the apparent horizontal position under the assumption 47 that the sound speed structure is horizontally stratified. Inversely, comparing the residuals of acoustic 48 travel time from multiple transponders, Osada et al. (2003) succeeded in estimating the temporal 49 variation of sound speed from the acoustic data. Kido et al. (2008) modified the expression to 50 validate the stationary configuration for a loosely tied buoy even in the case where the sound speed 51 has spatial variations. The stationary GNSS-A configuration is applied mainly by the groups in the 52 Scripps Institution of Oceanography (e.g., Gagnon et al., 2005;Chadwell and Spiess, 2008) and in 53 the Tohoku University (e.g., Fujimoto et al., 2014;Tomita et al., 2015). 54 On the other hand, Obana et al. (2000) and Asada and Yabuki (2001) took a "move-around" 55 approach where the 3-dimensional position of single transponder can be estimated by collecting the 56 acoustic data from various relay points on the sea surface. Figure  In order to enhance the stability of positioning, an assumption that the geometry of transponder array 68 is constant over whole observation period is usually adopted (e.g., Matsumoto  the averaged acoustic-ray direction, which results in the distortion of the estimated array geometry. 71 Constraining the array geometry contributes to reducing the bias error in the sound speed estimates 72 and the transponders' centroid position. 73 It should be noted that these two configurations are compatible under the adequate assumptions and 74 constraints. Recently, the group in the Tohoku University uses not only the stationary but also the 75 move-around observation data collected for determining the array geometry (Honsho and Kido,76 2017). 77

Recent improvements on GNSS-A analytical procedures 78
In the late 2010s, analytical procedures with the estimation of the spatial sound speed gradient for the 79 move-around configuration have been developed . In the earlier stage of the move-around GNSS-A  80  development, the spatial variations of sound speed were approximated as the temporal variations,  81  because most of the sound speed change are confined in the shallowest portion along the acoustic ray  82 paths (e.g., Watanabe and Uchida, 2016  The ATD offset values should be measured before the GNSS-A observation. 140

Underwater acoustic ranging 141
Another key subcomponent is the technique to measure the acoustic travel time between the sea-142 surface transducer and the seafloor transponders. The techniques for the precise ranging using 143 acoustic mirror-type transponders had been developed and practicalized in early studies (e.g., Spiess, 144 1980;Nagaya, 1995). Measuring round-trip travel time reduces the effect of advection of the media 145 between the instruments. 146 structure is assumed to be horizontally stratified, we can apply a heuristic approach based on the 155 Snell's law (e.g., Hovem, 2013), which has an advantage in computation time (e.g., Chadwell and 156 Sweeney, 2010; Sakic et al., 2018). 157 Therefore, we decomposed the 4-dimensional sound speed field into a horizontally stratified stational 158 sound speed profile and a perturbation to obtain the following travel time expression: 159

Sound speed perturbation model 174
In seawater, sound speed is empirically determined as a function of temperature, salinity, and 175 pressure (e.g., Del Grosso, 1974). Because these variables strongly depend on the water depth, the 176 vertical variation of the sound speed is much larger than the horizontal variation in the observation 177 scale. Thus, | | ≪ 1 will be satisfied in most cases where the reference sound speed appropriately 178 represents the sound speed field. In such cases, the average sound speed along the actual ray path is 179 expressed as 0 ̅ + ~ 0 ̅ + 0 ̅ , where 0 ̅ denotes the average sound speed of the reference 180 profile. 181 Recalling that the sound speed field is continuous and usually smooth in time and space compared to 182 the sampling rates of acoustic data, the acoustic ray path also has continuity in time and positions of 183 both ends, within the observation scale. It means that the acoustic rays from/to the neighboring ends 184 transmitted at almost the same time will take almost the same paths. Thus, can be modeled with a 185 smooth function of time and acoustic instruments' positions for the transmission and return paths, 186 i.e.,  vector, respectively. Let us assume that and follow a normal distribution with a variance-248 covariance of ( 2 ) and ( 2 ), whose scale can be adjusted by a hyperparameter 2 , i.e., = 249 2̃ and = 2̃, respectively. The prior probability density function (pdf) for the constraints 250 can be written as, 251 Then, the prior pdf can be written using the hyperparameter 〈•〉 , 261 where denotes the normalization constant. 263 Combining these prior informations, we obtain the following prior pdf: 264 and, 266 where and ‖ ‖ represent the rank of and the absolute value of the product of non-zero 268 eigenvalues of , respectively. 269

Variance-covariance of data 270
Now for the observed data, we take the assumption that also follows a normal distribution with a 271 variance-covariance of 2 ( , ), where and are the hyperparameters which control the 272 non-diagonal component of , i.e., 273 where is the number of data and |•| denotes the determinant of the matrix. temporal correlation which comes from the kinematic GNSS noise. The modeling error has spatio-280 temporal correlation because the sound speed variation is modeled by a smooth function of space and 281 time. Thus, we assumed the following covariance terms: 282 independent on the travel time value. Therefore, we apply = ( * ⁄ ) 2 , so that all measured data, 292 , has the same weight in the real scale. 293

Posterior probability 294
The posterior pdf after the data acquisition, which can be defined to be equal to the likelihood of the 295 model parameter given the data, can be written as, 296  Tables 1 and 2,  364 respectively. 365 We developed GARPOS to be compatible with both observation configurations. When handling the 366 GNSS-A data collected in the stationary configurations, we should process data with some 367 constraints on model parameters. Specifically, (1) upward components of transponders' positions 368 should be fixed to zero, and (2) spatial gradient components of the sound speed perturbation model 369 should not be solved, i.e., 1 = 2 = 0, because these parameters cannot be well resolved in 370 the stationary configuration. Although further parameter tuning may be required for optimization, 371 users can solve the seafloor position by GARPOS with the stationary data in addition to the move-372 around data. 373

5
Applications to the actual data 374

Data and settings 375
In order to verify the proposed analytical procedure, we reanalyzed the GNSS-A data at the sites 376 named "TOS2" and "MYGI" (Table 3 Along with the acoustic observations, the profiles of temperature and/or conductivity were measured 394 by CTD, XCTD or XBT probes several times. The reference sound speed profile, 0 ( ), was 395 calculated from the observed temperature and salinity profiles using the empirical relationship 396 proposed by Del Grosso (1974). To save the computational cost for ray tracing, the profile was 397 approximated into a polygonal curve with several tens of nodes ( Figure 5). 398 During a GNSS-A survey, the vessel sails on a pre-determined track over the seafloor transponder 399 array to collect geometrically balanced acoustic data (e.g., Figure 1). The along-track observation is 400 repeated several times by reversing the sailing direction in order to reduce the bias due to the errors 401 in the ATD offset. The along-track observation (called "subset", hereafter) is repeated several times, 402 with reversed sailing direction in order to reduce the bias due to the errors in the ATD offset. 403 During an observation cruise, it occasionally took more than a few weeks to collect sufficient 404 acoustic data at a single site due to weather conditions or other operational restrictions. Even so, we 405 compiled a single dataset per site per cruise for the static seafloor positioning in practice, because the 406 positional changes should be too small to detect. We call the collection of a single GNSS-A dataset 407 "observation epoch" or "epoch", hereafter. 408 We set the parameters for the numbers of basis functions, , , and , in equation 5, as 0 = 409 1 = 2 = 15 for both preprocess and main process. 410

Array geometry determination 411
In order to calculate the proper array geometry ̅̅̅ for the rigid-array constraint, we first determined The preliminary array-free positioning was also used for the verification of the collected data. We 423 eliminated the outliers whose discrepancies from the preliminary solution were larger than the 424 arbitrary threshold. We set the threshold to be 5 times as large as the root mean square value (RMS) 425 of the travel time residuals. 426

Hyperparameter search 427
In order to get the solution * , we should determine the appropriate values for the various 428 hyperparameters, i.e., 2 , , , 2 , 0 2 , 1 2 , 1 2 , 2 2 , and 2 2 . In the scheme of the ABIC 429 minimization, 2 can be determined analytically by equation 23. It is reasonable to assume 2 ≡ 430 1 2 = 1 2 = 2 2 = 2 2 because these hyperparameters control the smoothness of the spatial 431 sound speed structure. For the purpose of single positioning, should be a large number, for example 432 in meter-order. The large hardly changes the ABIC value and thus the solution. 433 In order to save the computational resources, we should further reduce the number of 434 hyperparameters. We tentatively put = 0.5. For the sound speed variations, we had to assume the 435 strong constancy of spatial sound speed structure to resolve them with the single transducer GNSS-A. 436 For this reason, we selected the ratio of 0 2 and 2 , as 2 = 0.1 0 2 . The last two hyperparameters, 437 and 0 2 , were determined with the grid search method. The tested values for and 0 2 are = 438 (0 min. , 0.5 min. , 1 min. , 2 min., 3 min. ) and 0 2 = (10 −3 , 10 −2 , 10 −1 , 10 0 , 10 1 , 10 2 ), 439 respectively. 440 TOS2 is located offshore in the south of Shikoku Island, southwestern Japan, above the source region 446 of the 1946 Nankaido earthquake (e.g., Sagiya and Thatcher, 1999) along the Nankai Trough. 447 According to Yokota and Ishikawa (2020), who investigated the transient deformations at the GNSS-448

Results
A sites along the Nankai Trough, no significant signal was detected at TOS2. The results by the 449 proposed method show the same trends as the conventional results. Although the trend of horizontal 450 displacement seems to be changed in 2018 or 2019, careful inspection is needed because the 451 transponders had been replaced during this period. 452 MYGI is located in the offshore east of Miyagi Prefecture, northeastern Japan, which experienced the 453 2011 Tohoku-oki earthquake (Sato et al., 2011). After the earthquake, significant westward 454 postseismic movement and subsidence due to the viscoelastic relaxation has been observed at MYGI 455 (Watanabe et al., 2014). The postseismic movements continue but appear to decay. It is true that the 456 changes in the displacement rate at these sites are crucial in seismic and geodetic researches, but 457 discussing these matters is beyond the scope of the present paper. The point is that the seafloor 458 positioning results were well reproduced by the proposed method.  gradient of slowness lies in the certain layer, from sea-surface to the depth , in the water. In such 510 cases, ( ) is propotional to ( ), as = ( 2 ⁄ ) . This is exactly the same assumption as the 511 model by Yasuda et al. (2017). The model of Kinugasa et al. (2020) is the special case of those 512 models where equals to the water depth. 513 In the proposed method, the sound speed field is approximately interpreted by their models when the 514 unit vector of is supposed to be same as that of and | | ≥ | |. The depth of the gradient 515 layer is calculated as, 516 = 2 * 1 + ⁄ (32) 517 When = , it concludes to the model of Kinugasa et al. (2020). Conversely, when | | ≪ | |, 518 sound speed gradient lies in the thin layer near the surface. 519 In addition to the simple model above, the proposed method can extract more complicated sound 520 speed field, which partly described by . Extracted parameters for the 521 sound speed perturbation indicate the complicity of oceanographic structure, as shown in the next 522 section. 523

Validity of extracted sound speed perturbation model 524
Typical examples for the estimation results for each observation, i.e., the time series of travel time 525 residuals, and sound speed perturbation interpreted from the correction coefficient, are shown in 526 Figure 7. Results for all the datasets are available in Supplementary Figure 1. 527 In the most cases for site TOS2, both terms of the estimated sound speed gradient vector stably direct 528 south to southeast. Because the sound speed increase with the water temperature, it means that the 529 water temperature is higher in the southern region. The results that 2 is comparable with 1 in many 530 cases indicate that the gradient of water temperature continues to the deeper portion, as discussed in 531 the previous section. This is consistent with the fact that the Kuroshio current continuously flows on 532 the south of TOS2. 533 In contrast, the directions of gradient terms at MYGI have less constancy than TOS2. Unlike the area 534 around TOS2 where the Kuroshio current dominantly affects the seawater structure, MYGI is located 535 in an area with a complicated ocean current system (e.g., Yasuda, 2003;Miyazawa et al., 2009 in other observation subcomponents such as the random walk noise in GNSS positioning, the drifts of 545 gyro sensor, or the time synchronization error between the devices. 546 Preferred models for all the tested epochs had positive values for data correlation length, . It 547 contributed to avoiding overfitting of the correction coefficient . It is considered that the plausible 548 estimation of sound speed is realized by introducing the statistic information criteria and the 549 information of data covariances. 550 Figure 8 shows the examples of the cases for the models without assuming the data correlation, i.e., 551 = 0. The preferred models were selected from 0 2 = (10 −3 , 10 −2 , 10 −1 , 10 0 , 10 1 , 10 2 , 10 3 , 10 4 ). 552 It is clear that the preferred models without assuming the data correlation have larger 0 2 . Although 553 the residuals of travel time were reduced in these models, overfittings occurred for each term of . 554 Comparing the preferred and less-preferred results, the existence of data covariance components 555 contributes to the selection of a model with less perturbation by decreasing the impact of individual 556 data on model parameters. 557 Finally, we confirm the stability of the seafloor positioning results. The differences of seafloor 558 position for the tested models from the most preferred models are summarized in Figure 9. The 559 differences in estimated positions for most of the tested models converged into several centimeters. 560 For both sites, variations in vertical component tend to be larger for larger values of 0 2 . It indicates 561 that finer hyperparameter tuning is not required when considering the application to seafloor 562 positioning. 563

Conclusions 564
We reconstructed the GNSS-A observation equation and developed the Python-based software 565 GARPOS to solve the seafloor position as well as the sound speed perturbations using the empirical 566 Bayes approach. It provides a stable solution for a generally ill-posed problem caused by the 567 correlation among the model parameters, by introducing the hyperparameter tuning based on the 568 ABIC minimization and data covariance to rationalize the normalization constant of the posterior pdf. 569 The most important point is that the proposed method succeeded in directly extracting the time-570 dependent sound speed field with two end members of spatial gradient terms, which are roughly 571 characterized by depths, even when the observers used only one sea-surface unit. Statistical approach 572 allowed us to suppress the overfitting and thus to obtain simpler sound speed field from densely 573 collected dataset. It successfully reproduced the stationary southward sound speed gradient at TOS2, 574 which is consistent with the Kuroshio current. 575 On the other hand, model overfits were shown in several epochs. These overfits can be caused not 576 only by the actually complicated sound speed field but also by other error sources which were not 577 well included in the model. It means that the hyperparameter tuning also plays a role in the 578 verification of dataset and model. Error analyses in such cases might rather help improving the 579 GNSS-A accuracy and methodology. 580 We suggested a simplified formatting for the GARPOS input files. Researchers can enter into the 581 field of seafloor geodesy by collecting the listed data with adequate precision. Since each 582 subcomponent of GNSS-A technique, i.e., GNSS positioning, acoustic ranging, and so on, has been 583 well established, observers can combine them on their platform. Especially, GNSS-A is expected to 584 be practicalized in the near future with an unmanned surface vehicle (Chadwell, 2016)