Centimeter-level-precision seafloor geodetic positioning model with self-structured empirical sound speed profile
具有自结构经验声速剖面的厘米级精度海底大地测量定位模型

Abstract  抽象

In-field Sound Speed Profile (SSP) measurement is still indispensable for achieving centimeter-level-precision Global Navigation Satellite System (GNSS)-Acoustic (GNSS-A) positioning in current state of the art. However, in-field SSP measurement on the one hand causes a huge cost and on the other hand prevents GNSS-A from global seafloor geodesy especially for real-time applications. We propose an Empirical Sound Speed Profile (ESSP) model with three unknown temperature parameters jointly estimated with the seafloor geodetic station coordinates, which is called the 1st-level optimization. Furthermore, regarding the sound speed variations of ESSP we propose a so-called 2nd-level optimization to achieve the centimeter-level-precision positioning for monitoring the seafloor tectonic movement. Long-term seafloor geodetic data analysis shows that, the proposed two-level optimization approach can achieve almost the same positioning result with that based on the in-field SSP. The influence of substituting the in-field SSP with ESSP on the horizontal coordinates is less than 3 mm, while that on the vertical coordinate is only 2–3 cm in the standard deviation sense.
在当前技术水平下，现场声速剖面 （SSP） 测量对于实现厘米级精度的全球导航卫星系统 （GNSS）-声学 （GNSS-A） 定位仍然是必不可少的。然而，现场 SSP 测量一方面会导致巨大的成本，另一方面又使 GNSS-A 无法进行全球海底大地测量，尤其是对于实时应用。我们提出了一个经验声速剖面 （ESSP） 模型，该模型具有三个未知温度参数，与海底大地测量站坐标联合估计，称为一级优化。此外，关于 ESSP 的声速变化，我们提出了一种所谓的二级优化，以实现厘米级精度定位，以监测海底构造运动。长期海底大地测量数据分析表明，所提出的两级优化方法能够获得与基于现场 SSP 的定位结果几乎相同的定位结果。用 ESSP 代替场内 SSP 对水平坐标的影响小于 3 mm，而在标准差意义上，对垂直坐标的影响仅为 2-3 cm。

Introduction  介绍

Global Navigation Satellite System (GNSS)-Acoustic (GNSS-A) positioning technique has become a vital tool for seafloor geodesy and crustal deformation applications of submarine offshore regions (Bürgmann & Chadwell, 2014; Iinuma et al., 2021). GNSS-A positioning technique was pioneered by Fred Spiess (Spiess, 1985; Spiess et al., 1998) and further developed for more than two decades (Chadwell & Sweeney, 2010; Chadwell et al., 1997). Regional seafloor geodetic networks were established in the past (Poutanen & Rózsa, 2020; Yang et al., 2020) and had provided key observations for constraining tectonic motion, crustal deformation and the earthquake cycle (Gagnon & Chadwell, 2007; Gagnon et al., 2005; Iinuma et al., 2021; Kido et al., 2006; Sato et al., 2011). GNSS-A technique achievements made in the past are on the one hand from observation scheme advances (Kido et al., 2008; Sato et al., 2013; Spiess et al., 1998) and on the other hand from acoustic positioning model improvements (Asada & Yabuki, 2006; Fujita et al., 2006; Spiess, 1980; Watanabe et al., 2020; Yang & Qin, 2020; Yokota et al., 2019). Nowadays, terrestrial geodesy has greatly facilitated low-cost and real-time positioning and near-real-time atmosphere state sensing due to positioning model advances (Torge & Müller, 2012), but this has not been achieved in seafloor geodesy. One of important reasons for this is that the high-precision seafloor geodetic positioning seriously relies on the Sound Speed Profile (SSP) measurement in field.
全球导航卫星系统（GNSS）-声学（GNSS-A）定位技术已成为海底大地测量和海底近海区域地壳变形应用的重要工具（Bürgmann & Chadwell，2014;Iinuma et al.， 2021）。GNSS-A 定位技术由 Fred Spiess （Spiess， 1985;Spiess 等人，1998 年），并进一步发展了二十多年（Chadwell & Sweeney，2010 年;Chadwell et al.， 1997）。区域海底大地测量网络在过去就已经建立起来了（Poutanen & Rózsa， 2020;Yang et al.， 2020），并为限制构造运动、地壳变形和地震周期提供了关键观测结果（Gagnon & Chadwell， 2007;Gagnon et al.， 2005;Iinuma et al.， 2021;Kido et al.， 2006;Sato et al.， 2011）。过去取得的 GNSS-A 技术成就一方面来自观测方案的进步（Kido et al.， 2008;Sato et al.， 2013;Spiess 等人，1998），另一方面来自声学定位模型的改进（Asada & Yabuki，2006;Fujita 等人，2006 年;Spiess，1980 年;Watanabe 等人，2020 年;Yang & Qin， 2020;Yokota 等人，2019 年）。如今，由于定位模型的进步，陆地大地测量学极大地促进了低成本和实时定位以及近实时大气状态传感（Torge & Müller，2012），但这在海底大地测量学中尚未实现。其中一个重要原因是，高精度海底大地测量严重依赖于现场声速剖面 （SSP） 测量。

The sound speed structure of the ocean water prominently varies from the depth, spatial location and the time (Yokota et al., 2022). Although several global ocean environment observation plans have been in operation (Hayes et al., 1991; Stammer et al., 2003; Zeng et al., 2016), e.g., Array for Real-time Geostrophic Oceanography (Argo) and Transparent Ocean Plan (TOP) (Wu et al., 2020), the current resolution of ocean environment observations, e.g., the Argo observations having a temporal resolution of several days (Roemmich et al., 2009), cannot satisfy the high-precision seafloor geodetic positioning demands; and therefore an in-field Reference SSP (RSSP) measurement is still required. The sound speed can be directly measured or derived from the Conductivity-Temperature-Depth (CTD) profiler measurement (Wilson & Wayne, 1960). However, it is still unrealistic to conduct a high-resolution Sound Speed Field (SSF) measurement for GNSS-A positioning because of the huge cost; and therefore GNSS-A positioning technique in the current state of the art adopts a RSSP to perform the seafloor geodetic positioning (Watanabe et al., 2020). For this, the positioning model requires a strong resist ability to remedy the sound speed variation effect of the RSSP.
海水的声速结构因深度、空间位置和时间而异（Yokota et al.， 2022）。尽管已经实施了几个全球海洋环境观测计划（Hayes et al.， 1991;Stammer 等人，2003 年;Zeng et al.， 2016），例如，实时地转海洋学阵列 （Argo） 和透明海洋计划 （TOP）（Wu et al.， 2020），目前海洋环境观测的分辨率，例如，Argo 观测具有几天的时间分辨率（Roemmich et al.， 2009），无法满足高精度海底大地测量定位需求;因此，仍然需要进行现场参考 SSP （RSSP） 测量。声速可以直接测量或从电导率-温度-深度（CTD）分析器测量中得出（Wilson & Wayne，1960）。然而，由于成本巨大，对 GNSS-A 定位进行高分辨率声速场 （SSF） 测量仍然不现实;因此，当前技术水平下的 GNSS-A 定位技术采用 RSSP 来执行海底大地测量定位（Watanabe 等人，2020 年）。为此，定位模型需要强大的光刻胶能力来补救 RSSP 的声速变化效应。

The acoustic positioning is seriously affected by the sound speed variations (Kido et al., 2008), which is very familiar to the atmosphere variation affecting GNSS positioning (Chen & Herring, 1997). Inspired by this, the Nadir Total Delay (NTD) model which is analogous to the zenith total delay in GNSS was proposed recently (Honsho & Kido, 2017; Honsho et al., 2019). This idea developed in the GNSS-A positioning can be also found in analyzing the sound speed variation for GNSS-A measurements (Kido et al., 2008). Besides NTD model, a generalized GNSS-A positioning model was also developed to correct the effect of the sound speed variation (Watanabe et al., 2020; Yokota et al., 2022). However, this kind of methods can extract only the sound speed variation relative to the in-field SSP.
声学定位受到声速变化的严重影响（Kido 等人，2008），这与影响 GNSS 定位的大气变化非常熟悉（Chen & Herring，1997）。受此启发，最近提出了类似于 GNSS 的天顶总延迟的最低点总延迟（NTD）模型（Honsho & Kido，2017 年;Honsho 等人，2019 年）。在 GNSS-A 定位中开发的这一想法也可以在分析 GNSS-A 测量的声速变化中找到（Kido et al.， 2008）。除了 NTD 模型外，还开发了广义 GNSS-A 定位模型来校正声速变化的影响（Watanabe et al.， 2020;Yokota 等人，2022 年）。然而，这种方法只能提取相对于场内 SSP 的声速变化。

In fact, the sound speed variation effect was very early studied by the GNSS-A and SSP observations in Hawaii (Osada et al., 2003). Then, the temporal variation inversion based on the RSSP was developed to eliminate this kind of effect (Fujita et al., 2004, 2006; Matsumoto et al., 2008). For more precisely characterizing the sound speed variation, the inversion method based on a 3-order B-spline model was developed in the past (Ikuta et al., 2008), and then Yokota et al. (2019) advanced this method by extracting the first-order and second-order horizontal sound speed gradients to represent more sound speed variation details (Yasuda et al., 2017). Note that, the above-mentioned inversion methods are based on the RSSP. In fact, the sound speed inversion without using the RSSP was also developed in the past (Chen, 2014). Recently, a self-constraint positioning method without the assistance of in-field SSP has also been developed (Zhao et al., 2022). However, without using the in-field SSP it is still hard to achieve a desirable positioning accuracy, e.g., the positioning error of the above self-constraint positioning method is up to 0.54 m. Precise positioning models without the assistance of in-field SSP still need to be developed to reduce the in-field SSP measurement cost, which is the vitally meaningful not only for facilitating the low-cost large-scale and even global seafloor geodesy but also for achieving the real-time seafloor geodetic applications, like nowadays mature GNSS geodesy.
事实上，夏威夷的 GNSS-A 和 SSP 观测很早就研究了声速变化效应（Osada et al.， 2003）。然后，开发了基于 RSSP 的时间变化反转来消除这种效应（Fujita et al.， 2004， 2006;Matsumoto et al.， 2008）。为了更精确地表征声速变化，过去开发了基于 3 阶 B 样条模型的反演方法（Ikuta et al.， 2008），然后 Yokota et al.（2019）通过提取一阶和二阶水平声速梯度来表示更多的声速变化细节（Yasuda et al.， 2017）. 请注意，上述倒置方法都是基于 RSSP 的。事实上，不使用 RSSP 的声速反演在过去也是开发的（Chen， 2014）。最近，还开发了一种无需场内 SSP 辅助的自约束定位方法（Zhao et al.， 2022）。然而，如果不使用现场 SSP，仍然难以达到理想的定位精度，例如，上述自约束定位方法的定位误差高达 0.54 m。在没有现场 SSP 帮助的情况下，仍然需要开发精确定位模型以降低现场 SSP 测量成本，这不仅对于促进低成本的大规模甚至全球海底大地测量具有重要意义，而且对于实现实时海底大地测量应用也具有重要意义，就像现在成熟的 GNSS 大地测量一样。

This contribution is to develop a Self-structured Empirical SSP (SESSP) approach to achieve centimeter-precision-level seafloor geodetic three-dimensional positioning. In "Three-parameter empirical SSP" section, a three-parameter Empirical Temperature Profile (ETP) model is proposed to structure an Empirical Sound Speed Profile (ESSP) by using the Del Grosso’s sound speed formula. In "Two-level optimizations on ESSP" section, a novel GNSS-A positioning model based two-level optimizations is proposed, of which the 1st-level optimization is to fix the ETP model parameters, and the 2nd-level optimization is to finally achieve high-precision location regarding the sound speed variations relative to ESSP. In "Experiment results and tests" section, the proposed models are verified by the long-term seafloor geodetic array observations.
这一贡献是开发了一种自结构经验 SSP （SESSP） 方法，以实现厘米级精度的海底大地三维定位。在“三参数经验 SSP”部分中，提出了一个三参数经验温度曲线 （ETP） 模型，以使用德尔格罗索声速公式构建经验声速曲线 （ESSP）。在“ESSP 的两级优化”部分中，提出了一种基于 GNSS-A 定位模型的新型两级优化，其中第一级优化是固定 ETP 模型参数，第二级优化是最终实现相对于 ESSP 的声速变化的高精度定位。在“实验结果和测试”部分，通过长期海底大地测量阵列观测验证了所提出的模型。

Three-parameter empirical SSP
三参数经验 SSP

The horizontal sound speed stratification can be expressed as an exponential SSP with unknowns for applying much of the World’s oceans (Munk, 1974). Besides, an empirical bilinear SSP $c_{\text{BP}}(u,{pp}_{\text{BP}})$ with four unknown parameters was also proposed and applied in the seafloor geodetic positioning, that is Chen (2014)
水平声速分层可以表示为指数 SSP，对于应用世界上大部分海洋来说，这是未知数（Munk，1974 年）。此外，还提出了一个具有四个未知参数的经验双线性 SSP $c_{\text{BP}}(u,{pp}_{\text{BP}})$ ，并将其应用于海底大地测量定位，即 Chen （ 2014）

$$c_{\text{BP}}(u,{pp}_{\text{BP}})=\{\begin{array}{ll}{v}_{\text{s}}+g_{\text{u}}u& 0\le u<u_{\text{b}}\\ v_{\text{b}}+g_{\text{d}}(u-u_{\text{b}})& u_{\text{b}}\le u\end{array}$$

(1)

where $u$ is the depth, ${pp}_{\text{BP}}=\left(\begin{array}{cccc}{v}_{\text{s}}& g_{\text{u}}& g_{\text{d}}& u_{\text{b}}\end{array}\right)$ is the unknown parameter vector, $v_{\text{s}}$ is the surface sound speed, $u_{\text{b}}$ is the bilinear break depth, $v_{\text{b}}$ is the sound speed corresponding to the depth $u_{\text{b}}$, and $g_{\text{u}}$ and $g_{\text{d}}$ are the piecewise gradients of the bilinear function. The unknown parameters can be estimated jointly with the seafloor geodetic station coordinates.
其中 $u$ 是深度， ${pp}_{\text{BP}}=\left(\begin{array}{cccc}{v}_{\text{s}}& g_{\text{u}}& g_{\text{d}}& u_{\text{b}}\end{array}\right)$ 是未知参数向量， $v_{\text{s}}$ 是表面声速， $u_{\text{b}}$ 是双线性断裂深度， $v_{\text{b}}$ 是对应于深度 $u_{\text{b}}$ 的声速， $g_{\text{u}}$ 和 $g_{\text{d}}$ 是双线性函数的分段梯度。未知参数可以与海底大地测量站坐标联合估计。

Chen (2014) pointed out that $v_{\text{s}}$ and $u_{\text{b}}$ were treated knowns for the facilitation of surface sound speed measurement, but this is not so practical in some cases. Note that there are a series of equivalent profiles that can be used to obtain almost the same positioning result (Sun et al., 2019; Zielinski & Geng, 1999), i.e., this is an ill-posed problem widely existed in nonlinear inversion. To obtain a meaningful solution, we should impose a prior information on ${pp}_{\text{BP}}$ with a certain uncertainty, see the Table 3. For solving this problem, we however present a novel ESSP by Del Grosso sound speed formula with an ETP as follow
Chen （ 2014） 指出， $v_{\text{s}}$ 和 $u_{\text{b}}$ Treated Known Known for facilitate of Surface Sound Speed Measurement，但这在某些情况下并不那么实用。请注意，有一系列等效的轮廓可用于获得几乎相同的定位结果（Sun 等人，2019 年;Zielinski & Geng， 1999），即这是一个广泛存在于非线性反演中的病态问题。为了获得有意义的解决方案，我们应该对 ${pp}_{\text{BP}}$ 具有一定不确定性的先验信息施加，见表 3。然而，为了解决这个问题，我们提出了一种新颖的 ESSP by Del Grosso 声速公式，其 ETP 如下

$$T_{\text{ETP}}(u,{pp}_{\text{ETP}})=\tau p{\tau p}_{\text{ETP}}=T_{\text{m}}+\mathrm{\Delta }T{\text{e}}^{-\frac{u}{u_0}}$$

(2)

where ${pp}_{\text{ETP}}={\left(T_{\text{m}}\mathrm{\Delta }T\right)}^{\text{T}}$ is an unknown parameter vector of ETP to be jointly estimated with the seafloor geodetic station coordinates, $\tau \tau =\left(T_{\text{m}}{\text{e}}^{-u/u_0}\right)$ is the corresponding coefficient matrix, $T_{\text{m}}$ represents the intermediate overmeasurement, $\mathrm{\Delta }T$ represents the temperature difference between the sea-surface and sea-bottom, $u_0$ represents the depth of the thermocline, see Fig. 1. The $u_0$ and $T_{\text{m}}$ can be statistically induced from the long-term ocean environment observations and thereby we can impose prior constraints on the parameter estimation.
其中 ${pp}_{\text{ETP}}={\left(T_{\text{m}}\mathrm{\Delta }T\right)}^{\text{T}}$ 是要与海底大地测量站坐标联合估计的 ETP 的未知参数向量， $\tau \tau =\left(T_{\text{m}}{\text{e}}^{-u/u_0}\right)$ 是相应的系数矩阵， $T_{\text{m}}$ 表示中间超测， $\mathrm{\Delta }T$ 表示海面和海底的温度差， $u_0$ 表示温跃层的深度，见图 1。和 $u_0$ $T_{\text{m}}$ 可以从长期海洋环境观测中统计推断出来，因此我们可以对参数估计施加先验约束。

Fig. 1  图 1

Illustration on the empirical temperature profile
经验温度曲线图

As shown in Fig. 1, the three parameters $T_{\text{m}}$, $\mathrm{\Delta }T$ and $u_0$ to be estimated are the control parameters of the ETP shape. Note that the deep-sea bottom water temperature is generally very stable, for which it is convenient to impose a priori knowledge or constraint on the model parameter $T_{\text{m}}$. Besides, the deep-sea water temperature will decrease 1–2 °C per 1000 m, which might also be useful to characterize the ETP. Then, substituting ETP (2) and the average salinity S = 35‰ into the Del Grosso formula (Grosso, 1974; Wong & Zhu, 1995), we can establish an ESSP as
如图 1 所示，三个参数 $T_{\text{m}}$ 和 $\mathrm{\Delta }T$ $u_0$ 待估计 是 ETP 形状的控制参数。注意，深海海底水温一般都很稳定，方便对模型参数 $T_{\text{m}}$ 施加先验知识或约束。此外，海洋深水温度每 1000 m 将降低 1–2 °C，这也可能有助于表征 ETP。然后，将 ETP （ 2） 和平均盐度 S = 35‰ 代入 Del Grosso 公式（Grosso，1974 年;Wong & Zhu， 1995），我们可以将 ESSP 建立为

$$c_{\text{ESSP}}(u,{pp}_{\text{ETP}})={\beta }_{\text{C}}+{\beta }_{\text{P}}+{\beta }_{\text{T}}({pp}_{\text{ETP}})$$

(3)

where $c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ is ESSP,
其中 $c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ 是 ESSP，

$${\beta }_{\text{C}}={\text{C}}_{\text{c}}+{\text{C}}_{\text{s}(1)}S+{\text{C}}_{\text{s}(2)}{S}^2$$

(4)

is a constant associated with the salinity S, ${\text{C}}_{\text{c}}$ is the constant term in the Del Grosso formula,
是与盐度 S 相关的常数， ${\text{C}}_{\text{c}}$ 是 Del Grosso 公式中的常数项，

$$\begin{array}{rl}{\beta }_{\text{P}}& ={\text{C}}_{\text{s(2)p(2)}}{S}^2{P}^2+{\text{C}}_{\text{p(1)}}P+{\text{C}}_{\text{p(2)}}{P}^2+{\text{C}}_{\text{p(3)}}{P}^3\\ & {\text{ = C}}_{\text{p(1)}}P+{\text{(C}}_{\text{s(2)p(2)}}{S}^2+{\text{C}}_{\text{p(2)}})P^2+{\text{C}}_{\text{p(3)}}{P}^3\\ & ={\text{C}}_{\text{p(1)}}P+{\text{C}}_{\text{d(1)}}{P}^2+{\text{C}}_{\text{p(3)}}{P}^3\end{array}$$

(5)

is the sound speed variation associated with the pressure P which is recommended to use the Leroy’s formula as Leroy and Parthiot (1998)
是与压力 P 相关的声速变化，建议使用 Leroy 公式，如 Leroy 和 Parthiot （ 1998）

$$P=1.0052405(1+5.28\times {10}^{-3}{\sin }^2\text{φ})u+2.36\times {10}^{-6}{u}^2$$

(6)

where $\text{φ}$ is the latitude,
其中 $\text{φ}$ 是纬度，

$$\begin{array}{rl}{\beta }_{\text{T}}({pp}_{\text{ETP}})& ={\text{C}}_{\text{t(1)}}{T}_{\text{ETP}}({pp}_{\text{ETP}})+{\text{C}}_{\text{t(2)}}{T}_{\text{ETP}}^2({pp}_{\text{ETP}})+{\text{C}}_{\text{t(3)}}{T}_{\text{ETP}}^3({pp}_{\text{ETP}})\\ & +{\text{C}}_{\text{tp}}P{T}_{\text{ETP}}({pp}_{\text{ETP}})+{\text{C}}_{\text{t(3)p}}P{T}_{\text{ETP}}^3({pp}_{\text{ETP}})+{\text{C}}_{\text{tp(2)}}{P}^2{T}_{\text{ETP}}({pp}_{\text{ETP}})+{\text{C}}_{\text{t(2)p(2)}}{P}^2{T}_{\text{ETP}}^2({pp}_{\text{ETP}})+{\text{C}}_{\text{tp(3)}}{P}^3{T}_{\text{ETP}}({pp}_{\text{ETP}})\\ & +{\text{C}}_{\text{st}}S{T}_{\text{ETP}}({pp}_{\text{ETP}})+{\text{C}}_{\text{st(2)}}S{T}_{\text{ETP}}^2({pp}_{\text{ETP}})+{\text{C}}_{\text{stp}}SP{T}_{\text{ETP}}({pp}_{\text{ETP}})+{\text{C}}_{\text{s(2)tp}}{S}^{2}P{T}_{\text{ETP}}({pp}_{\text{ETP}})\\ & =({\text{C}}_{\text{t(1)}}+{\text{C}}_{\text{st}}S)T_{\text{ETP}}({pp}_{\text{ETP}})+(({\text{C}}_{\text{tp}}+{\text{C}}_{\text{stp}}S+{\text{C}}_{\text{s(2)tp}}{S}^2)P+{\text{C}}_{\text{t(1)p(2)}}{P}^2+{\text{C}}_{\text{t(1)p(3)}}{P}^3)T_{\text{ETP}}({pp}_{\text{ETP}})\\ & +{\text{ (C}}_{\text{st(2)}}S+{\text{C}}_{\text{t(2)}})T_{\text{ETP}}^2({pp}_{\text{ETP}})+{\text{C}}_{\text{t(2)p(2)}}{P}^2{T}_{\text{ETP}}^2({pp}_{\text{ETP}})\\ & +{\text{C}}_{\text{t(3)}}{T}_{\text{ETP}}^3({pp}_{\text{ETP}})+{\text{C}}_{\text{t(3)p}}P{T}_{\text{ETP}}^3({pp}_{\text{ETP}})\\ & ={\text{C}}_{\text{d(2)}}{T}_{\text{ETP}}({pp}_{\text{ETP}})+({\text{C}}_{\text{d(3)}}P+{\text{C}}_{\text{t(1)p(2)}}{P}^2+{\text{C}}_{\text{t(1)p(3)}}{P}^3)T_{\text{ETP}}({pp}_{\text{ETP}})\\ & +{\text{C}}_{\text{d(4)}}{T}_{\text{ETP}}^2({pp}_{\text{ETP}})+{\text{C}}_{\text{t(2)p(2)}}{P}^2{T}_{\text{ETP}}^2({pp}_{\text{ETP}})\\ & +{\text{C}}_{\text{t(3)}}{T}_{\text{ETP}}^3({pp}_{\text{ETP}})+{\text{C}}_{\text{t(3)p}}P{T}_{\text{ETP}}^3({pp}_{\text{ETP}})\end{array}$$

(7)

is the sound speed variation determined by the water temperature T and the temperature-depth mixture terms.
是由水温 T 和温度-深度混合物项确定的声速变化。

With S = 35‰, we can then figure out the coefficients of Eqs. (3), (4), (5), (7) and they are given in Table 1 for facilitating the calculation. Finally, for an arbitrary latitude specified we can figure out the ESSP by employing the data in Table 1.
当 S = 35‰ 时，我们可以计算出 Eqs 的系数。（ 3）、（ 4）、（ 5）、 （ 7） 中给出，它们列于表 1 中，以便于计算。最后，对于指定的任意纬度，我们可以通过使用表 1 中的数据来计算 ESSP。

Table 1 ESSP coefficients
表 1 ESSP 系数

As the proposed ESSP regards the sound speed variation with hydrostatic pressure of the water that implied in the Leroy’s formula and the water temperature’s exponential decaying characteristic with the depth that implied in the proposed empirical temperature profile, it is more meaningful and accurate to perform the SSP inversion result. We can also use other sound speed formulae to structure the ESSP complying with their applicable conditions, e.g., Wilson formula (Wilson, 1960) and Chen-Millero formula (Chen & Millero, 1977). In the following section, we will propose an optimization approach to determine the model parameters of ETP (2).
由于所提出的 ESSP 考虑了 Leroy 公式中隐含的水的声速随静水压力的变化以及所提出的经验温度剖面中隐含的水温与深度的指数衰减特性，因此执行 SSP 反演结果更有意义和准确。我们还可以使用其他声速公式来构建符合其适用条件的 ESSP，例如，Wilson 公式（Wilson，1960）和 Chen-Millero 公式（Chen & Millero，1977）。在下一节中，我们将提出一种优化方法来确定 ETP 的模型参数 （2）。

Two-level optimizations on ESSP
ESSP 上的两级优化

Overall Scheme  整体方案

Next, we will conduct a joint estimation of ESSP parameters and the seafloor geodetic station coordinates by employing the 1st-level optimization as shown in Fig. 2, which is called as sound speed self-structured empirical SSP approach for its independence on the in-field SSP measurement.
接下来，我们将采用如图 2 所示的一级优化，对 ESSP 参数和海底大地测量站坐标进行联合估计，这被称为声速自结构化经验 SSP 方法，因为它独立于现场 SSP 测量。

Fig. 2  图 2

Overall scheme of the two-level optimization
两级优化总体方案

As shown in Fig. 2, the 1st-level optimization is performed in the context of the ray-tracing positioning procedure while the 2nd-level one is performed by B-splines for characterizing acoustic delay caused by the sound speed variations relative to ESSP.
如图 2 所示，第一级优化是在射线追踪定位程序的背景下进行的，而第二级优化是由 B 样条执行的，用于表征由相对于 ESSP 的声速变化引起的声学延迟。

The 1st-level optimization
第一级优化

GNSS-A positioning is achieved by a combination of the sea-surface GNSS positioning with the precise acoustic round-trip time measurement from the sea-surface acoustic transducer to the seafloor acoustic transponder (Asada & Yabuki, 2006). Due to the horizontal density stratification of the ocean, the seafloor geodetic positioning generally uses the ray-tracing positioning model. For taking place of the high-cost in-field SSP measurement, we use ESSP $c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ to perform a joint estimation of ${pp}_{\text{ETP}}$ and the seafloor geodetic coordinates. This time that the GNSS-A positioning model reads
GNSS-A 定位是通过结合海面 GNSS 定位和从海面声学换能器到海底声学转发器的精确声学往返时间测量来实现的（Asada & Yabuki， 2006）。由于海洋的水平密度分层，海底大地定位一般采用射线追踪定位模型。为了进行高成本的现场 SSP 测量，我们使用 ESSP $c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ 对海底大地坐标进行联合估计 ${pp}_{\text{ETP}}$ 。这一次，GNSS-A 定位模型读取

$$T_{\text{obs},i}=T_i\left(XX,c_{\text{ESSP}}(u,{pp}_{\text{ETP}})\right)+{\text{ε}}_{\text{T},i}$$

(8)

where $T_{\text{obs},i}$ is the ith round-trip time observation, $XX$ is the seafloor transponder coordinate vector to be estimated, ${\text{ε}}_{\text{T},i}$ is the random error of observation. $T_i=T_{\text{s},i}+T_{\text{r},i}$ represents calculated the ith round-trip time,
其中 $T_{\text{obs},i}$ 是第 i 个往返时间观测， $XX$ 是要估计的海底应答器坐标向量， ${\text{ε}}_{\text{T},i}$ 是观测的随机误差。 $T_i=T_{\text{s},i}+T_{\text{r},i}$ 表示计算的第 i 次往返时间，

$$T_{J,i}={\int }_{u_{\text{X}}}^{u_{\text{x}}^J}\frac{1}{\cos {\beta }_i(u)}\frac{1}{c_{\text{ESSP}}(u,{pp}_{\text{ETP}})}\text{d}uJ\in \{\text{s},\text{r}\}$$

(9)

is one half of the round-trip travelling time calculated by the Two Dimensional (2D) ray-tracing model (Spiess, 1980), $u_{\text{x}}^J,u_{\text{X}}$ are depths of the transducer and transponder, respectively. Note that, horizontal coordinates of transponder are implied in the incident angle ${\beta }_i(u)$ of the ith ray calculated by the inversion of the eigenray (Yang et al., 2021).
是由二维 （2D） 光线追踪模型 （Spiess， 1980） 计算的往返旅行时间的一半， $u_{\text{x}}^J,u_{\text{X}}$ 分别是换能器和应答器的深度。请注意，应答器的水平坐标隐含在由特征射线反转计算的第 i 条射线的入射角 ${\beta }_i(u)$ 中（Yang 等人，2021 年）。

For n observations we have an over-determined vector-form observation equation ${TT}_{\text{obs}}=TT\left(XX,c_{\text{ESSP}}(u,{pp}_{\text{ETP}})\right)+{\text{ε}\text{ε}}_{\text{T}}$ solved by the nonlinear least squares (LS) criterion reads:
对于 n 个观测值，我们有一个由非线性最小二乘法 （LS） 准则 ${TT}_{\text{obs}}=TT\left(XX,c_{\text{ESSP}}(u,{pp}_{\text{ETP}})\right)+{\text{ε}\text{ε}}_{\text{T}}$ 求解的超定向量形式观测方程 ，如下所示：

$$\underset{{pp}_{\text{ETP}}}{min}{g}_1(XX,{pp}_{\text{ETP}}):={VV}^{\text{T}}(XX,{pp}_{\text{ETP}}){\text{Σ}\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})$$

(10)

where $g_1$ is the weighted sum of the squared residuals, $VV(XX,{pp}_{\text{ETP}}):={TT}_{\text{obs}}-TT\left(XX,c_{\text{ESSP}}(u,{pp}_{\text{ETP}})\right)$ is the residual vector, and ${\text{Σ}\text{Σ}}_{\text{L}}^{}={PP}^{-1}{\text{σ}}_0^2$ and ${\text{σ}}_0^2$ are the variance of observations and prior unit weight variance, respectively. $PP=\text{diag}(p_1,p_2,\dots ,p_n)$ is the weight matrix of observations where $p_i=(\sin (e_i){)}^2$ is the weight of the ith observation and $e_i$ is the corresponding elevation angle, respectively. Hereafter, we however adopt an equal weight matrix to simplify the discussion, i.e., $p_i=1$. To obtain LS solution we get the following first-order partial derivatives of $g_1\left(XX,{pp}_{\text{ETP}}\right)$, that is
其中 $g_1$ 是残差平方和的加权和， $VV(XX,{pp}_{\text{ETP}}):={TT}_{\text{obs}}-TT\left(XX,c_{\text{ESSP}}(u,{pp}_{\text{ETP}})\right)$ 是残差向量， ${\text{Σ}\text{Σ}}_{\text{L}}^{}={PP}^{-1}{\text{σ}}_0^2$ 和 ${\text{σ}}_0^2$ 分别是观测值的方差和先验单位权重方差。 $PP=\text{diag}(p_1,p_2,\dots ,p_n)$ 是观测值的权重矩阵，其中 $p_i=(\sin (e_i){)}^2$ 是第 i 个观测值的权重，分别 $e_i$ 是相应的仰角。然而，在下文中，我们采用一个等权重矩阵来简化讨论，即 $p_i=1$ 。为了获得 LS 解，我们得到以下 $g_1\left(XX,{pp}_{\text{ETP}}\right)$ 的一阶偏导数 ，即

$$hh(XX,{pp}_{\text{ETP}})=\left(\begin{array}{c}\mathrm{\partial }{g}_1/\mathrm{\partial }XX\\ \mathrm{\partial }{g}_1/\mathrm{\partial }{pp}_{\text{ETP}}\end{array}\right)=\left(\begin{array}{l}{AA}^{\text{T}}(XX){\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})\\ {BB}_{}^{\text{T}}({pp}_{\text{ETP}}){\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})\end{array}\right)$$

(11)

where  哪里

$$AA(XX)=\frac{\mathrm{\partial }TT}{\mathrm{\partial }XX}=\left(\begin{array}{c}\sum _{J\in \{\text{s},\text{r}\}}\frac{\mathrm{\partial }{T}_{J,1}}{\mathrm{\partial }XX}\\ \sum _{J\in \{\text{s},\text{r}\}}\frac{\mathrm{\partial }{T}_{J,2}}{\mathrm{\partial }XX}\\ ⋮\\ \sum _{J\in \{\text{s},\text{r}\}}\frac{\mathrm{\partial }{T}_{J,n}}{\mathrm{\partial }XX}\end{array}\right)=c_{\text{X}}^{-1}\left(\begin{array}{c}\sum _{J\in \{\text{s},\text{r}\}}\left(\begin{array}{ccc}\sin {\alpha }_{J,1}\sin {\beta }_{J,1}& \cos {\alpha }_{J,1}\sin {\beta }_{J,1}& \cos {\beta }_{J,1}\end{array}\right)\\ \sum _{J\in \{\text{s},\text{r}\}}\left(\begin{array}{ccc}\sin {\alpha }_{J,2}\sin {\beta }_{J,2}& \cos {\alpha }_{J,1}\sin {\beta }_{J,2}& \cos {\beta }_{J,2}\end{array}\right)\\ ⋮\\ \sum _{J\in \{\text{s},\text{r}\}}\left(\begin{array}{ccc}\sin {\alpha }_{J,n}\sin {\beta }_{J,n}& \cos {\alpha }_{J,n}\sin {\beta }_{J,n}& \cos {\beta }_{J,n}\end{array}\right)\end{array}\right)$$

(12)

is the Jacobian matrix of the nonlinear observation about the coordinate vector, $c_{\text{X}}^{}$ is the sound speed at the depth of the transponder, ${\alpha }_{J,i},{\beta }_{J,i}$ are the azimuth angle and incident angle of the ith ray, respectively.
是非线性观测关于坐标矢量的雅可比矩阵， $c_{\text{X}}^{}$ 是应答器深度处的声速， ${\alpha }_{J,i},{\beta }_{J,i}$ 分别是第 i 条射线的方位角和入射角。

$$BB({pp}_{\text{ETP}})=\frac{\mathrm{\partial }TT}{\mathrm{\partial }{pp}_{\text{ETP}}}=\frac{\mathrm{\partial }TT}{\mathrm{\partial }{c}_{\text{ESSP}}}\frac{\mathrm{\partial }{c}_{\text{ESSP}}}{\mathrm{\partial }{pp}_{\text{ETP}}}$$

(13)

is the first-order derivative of the round-trip time $TT$ about ${pp}_{\text{ETP}}$. $BB(XX)$ is hardly analytically given but it can be numerically solved. Then, vanishing $hh(XX,{pp}_{\text{ETP}})$, i.e., we have
是往返时间 $TT$ 关于 ${pp}_{\text{ETP}}$ 的一阶导数。 $BB(XX)$ 很难分析给出，但可以用数值求解。然后，消失 $hh(XX,{pp}_{\text{ETP}})$ ，即我们有

$$\{\begin{array}{c}{AA}^{\text{T}}(XX){\text{Σ}\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})=0\\ {BB}_{}^{\text{T}}({pp}_{\text{ETP}}){\text{Σ}\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})=0\end{array}$$

(14)

that can be used to obtain the LS solution. With initial values ${XX}_0,{pp}_{\text{ETP}(0)}^{}$, if $g_1^{\mathrm{\prime }\mathrm{\prime }}({XX}_0,{pp}_{\text{ETP}(0)}^{})$ is positively defined, then LS solution can be locally solved by Newton’s method with the second-order partial derivative of $g_1(XX,{pp}_{\text{ETP}})$, that reads (Xue et al., 2014):
可用于获取 LS 解。对于初始值 ${XX}_0,{pp}_{\text{ETP}(0)}^{}$ ，如果 $g_1^{\mathrm{\prime }\mathrm{\prime }}({XX}_0,{pp}_{\text{ETP}(0)}^{})$ 是正定义的，则 LS 解可以通过牛顿方法局部求解，该方法具有 $g_1(XX,{pp}_{\text{ETP}})$ 的二阶偏导数 ，如下所示（Xue et al.， 2014）：

$$\left(\begin{array}{c}{XX}_{k+1}\\ {pp}_{\text{ETP},k+1}\end{array}\right)=\left(\begin{array}{c}{XX}_k\\ {pp}_{\text{ETP},k}\end{array}\right)+{\left(\begin{array}{ll}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}+{\Gamma \Gamma }_{\text{X}}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}+{\Gamma \Gamma }_{\text{p}}\end{array}\right)}^{-1}\left(\begin{array}{c}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\end{array}\right)$$

(15)

where ${AA}_k^{}:=AA({XX}_k)$, ${BB}_k^{}:=BB({pp}_{\text{ETP},k})$, ${VV}_k:=VV({XX}_k,{pp}_{\text{ETP},k})$ and k is the iteration index,
其中 ${AA}_k^{}:=AA({XX}_k)$ ， ${BB}_k^{}:=BB({pp}_{\text{ETP},k})$ ， 和 ${VV}_k:=VV({XX}_k,{pp}_{\text{ETP},k})$ k 是迭代索引，

$${\Gamma \Gamma }_{\text{X}}=\sum _{i=1}^n{V}_i(XX,{pp}_{\text{ETP}})p_i{\text{σ}}_0^{-2}{SS}_i(XX)$$

(16)

in which ${SS}_i(XX)=\mathrm{\partial }{aa}_i^{\text{T}}(XX)/\mathrm{\partial }XX$ is the first-order derivative of the ith row ${aa}_i(XX)$ of $AA(XX)$, i.e., the Hessian matrix of the round-trip time $T_i(XX)$ about $XX$,
其中 ${SS}_i(XX)=\mathrm{\partial }{aa}_i^{\text{T}}(XX)/\mathrm{\partial }XX$ 是 的第 i 行 ${aa}_i(XX)$ 的一阶导数 $AA(XX)$ ，即往返时间 $T_i(XX)$ 约为 $XX$ 的 Hessian 矩阵

$${\Gamma \Gamma }_{\text{p}}={\sum }_{i=1}^n{V_i(XX,{pp}_{\text{ETP}})p_i{\text{σ}}_0^{-2}QQ}_i({pp}_{\text{ETP}})$$

(17)

in which ${QQ}_i({pp}_{\text{ETP}})=\mathrm{\partial }{bb}_i^{\text{T}}({pp}_{\text{ETP}})/\mathrm{\partial }{pp}_{\text{ETP}}$ is the first-order derivative of the ith row ${bb}_i(XX)$ of $BB({pp}_{\text{ETP}})$, i.e., the Hessian matrix of the round-trip time ${TT}_i({pp}_{\text{ETP}})$ about ${pp}_{\text{ETP}}$. Note that, ${\Gamma \Gamma }_{\text{X}}$ and ${\Gamma \Gamma }_{\text{p}}$ can be ignored for long-distance or small-residual cases, and at this time we recommend using Gauss–Newton iterative formula as
其中 ${QQ}_i({pp}_{\text{ETP}})=\mathrm{\partial }{bb}_i^{\text{T}}({pp}_{\text{ETP}})/\mathrm{\partial }{pp}_{\text{ETP}}$ 是 的第 i 行 ${bb}_i(XX)$ 的一阶导数 $BB({pp}_{\text{ETP}})$ ，即往返时间 ${TT}_i({pp}_{\text{ETP}})$ 约为 ${pp}_{\text{ETP}}$ 的 Hessian 矩阵。请注意， ${\Gamma \Gamma }_{\text{X}}$ 对于长距离或小残差情况，可以忽略 and ${\Gamma \Gamma }_{\text{p}}$ ，此时我们建议使用高斯-牛顿迭代公式 as

$$\left(\begin{array}{c}{XX}_{k+1}\\ {pp}_{\text{ETP},k+1}\end{array}\right)=\left(\begin{array}{c}{XX}_k\\ {pp}_{\text{ETP},k}\end{array}\right)+{\left(\begin{array}{cc}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}\end{array}\right)}^{-1}\left(\begin{array}{c}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\end{array}\right)$$

(18)

which is recommended to be terminated under the condition $\underset{i=1,2,\dots ,S}{max}\left(\Vert {XX}_{i,k+1}-{XX}_{i,k}\Vert \right)<\delta$ where $\delta$ is a sufficient small positive value, i is the index of the S seafloor geodetic stations. It is generally necessary to introduce a certain number of constraints on partial parameters to stabilize the iteration or to obtain a meaningful solution.
建议在足够小的正值的情况下终止 $\underset{i=1,2,\dots ,S}{max}\left(\Vert {XX}_{i,k+1}-{XX}_{i,k}\Vert \right)<\delta$ $\delta$ ，i 是 S 海底大地测量台站的索引。通常需要对 partial 参数引入一定数量的约束，以稳定迭代或获得有意义的解决方案。

Next, we impose a priori constraints on ${pp}_{\text{ETP}}$. The first a priori constraint can be structured by the temperature stability of the seawater-bottom, that is
接下来，我们对 施加先验约束 ${pp}_{\text{ETP}}$ 。第一个先验约束可以由海水底部的温度稳定性构成，即

$$T_{\text{b}}=T_{\text{ETP}}(u_{\text{X}},{pp}_{\text{ETP}})+{\text{ε}}_{\text{b}}^{}$$

(19)

where $u_{\text{X}}$ is the depth of the seafloor geodetic station, the error ${\text{ε}}_{\text{b}}^{}$ represents the uncertainty of the a priori seawater-bottom temperature. Note that because deep-sea temperature is generally quite stable, the average temperature of the seawater bottom can be easily obtained by global ocean environment observations.
其中 $u_{\text{X}}$ 是海底测地测站的深度，误差 ${\text{ε}}_{\text{b}}^{}$ 表示先验海水底部温度的不确定性。请注意，由于深海温度通常相当稳定，因此可以通过全球海洋环境观测轻松获得海水底部的平均温度。

Let ${\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}$ be the variance of $T_{\text{b}}$, we can structure the following LS criterion
设 ${\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}$ 是 的方差 $T_{\text{b}}$ ，我们可以构建以下 LS 准则

$$min{g}_1(ZZ)\mathbf{=}{VV}^{\text{T}}(XX,{pp}_{\text{ETP}}){\text{Σ}\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})+(T_{\text{b}}-\tau p{\tau p}_{\text{ETP}}{)}^{\text{T}}({\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}{)}^{-1}(T_{\text{b}}-\tau p{\tau p}_{\text{ETP}})$$

(20)

where $ZZ=({XX}^{\text{T}}{pp}_{\text{ETP}}^{\text{T}}{)}^{\text{T}}$ is the parameter vector to be estimated. Omitting the deduction, we can obtain the Gauss–Newton iterative formula as
其中 $ZZ=({XX}^{\text{T}}{pp}_{\text{ETP}}^{\text{T}}{)}^{\text{T}}$ 是要估计的参数向量。省略推导，我们可以得到高斯-牛顿迭代公式为

$$\left(\begin{array}{c}{XX}_{k+1}\\ {pp}_{\text{ETP},k+1}\end{array}\right)=\left(\begin{array}{c}{XX}_k\\ {pp}_{\text{ETP},k}\end{array}\right)+{\left(\begin{array}{ll}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}+{\tau \tau }^{\text{T}}({\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}{)}_{}^{-1}\tau \tau \end{array}\right)}^{-1}\left(\begin{array}{l}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k+({\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}{)}_{}^{-1}(T_{\text{b}}-\tau p{\tau p}_{\text{ETP},k})\end{array}\right)$$

(21)

which is recommended to use the same termination condition. It is generally hard to obtain a global or meaningful solution of the problem without an effective initial value.
建议使用相同的终止条件。如果没有有效的初始值，通常很难获得问题的全局或有意义的解决方案。

For most cases an arbitrary guess value of the unknown ${pp}_{\text{ETP}}$ is sufficient to start the iteration, but we still recommend using a coarse grid search approach to obtain a robust initial value. Let ${pp}_{\text{ETP}(0)}$ and ${\text{Σ}\text{Σ}}_{\text{p}}^{\text{ETP}}$ be the initial value and the variance of ${pp}_{\text{ETP}}$, respectively, at this time we can structure the following LS criterion
在大多数情况下，未知 ${pp}_{\text{ETP}}$ 数的任意猜测值足以开始迭代，但我们仍然建议使用粗略网格搜索方法来获得稳健的初始值。设 ${pp}_{\text{ETP}(0)}$ 和 ${\text{Σ}\text{Σ}}_{\text{p}}^{\text{ETP}}$ 分别为 的初值和方差 ${pp}_{\text{ETP}}$ ，此时我们可以构建以下 LS 准则

$$\begin{array}{rl}min{g}_1(ZZ)& \mathbf{=}{VV}^{\text{T}}(XX,{pp}_{\text{ETP}}){\text{Σ}\text{Σ}}_{\text{L}}^{-1}VV(XX,{pp}_{\text{ETP}})+(T_{\text{b}}-\tau p{\tau p}_{\text{ETP}}{)}^{\text{T}}({\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}{)}_{}^{-1}(T_{\text{b}}-\tau p{\tau p}_{\text{ETP}})\\ & +({pp}_{\text{ETP}(0)}-{pp}_{\text{ETP}}{)}^{\text{T}}({\text{Σ}\text{Σ}}_{\text{p}}^{\text{ETP}}{)}_{}^{-1}({pp}_{\text{ETP}(0)}-{pp}_{\text{ETP}})\end{array}$$

(22)

Omitting the deduction, we can write out the Gauss–Newton iterative formula as
省略推导，我们可以将高斯-牛顿迭代公式写为

$$\left(\begin{array}{c}{XX}_{k+1}\\ {pp}_{\text{ETP},k+1}\end{array}\right)=\left(\begin{array}{c}{XX}_k\\ {pp}_{\text{ETP},k}\end{array}\right)+{\left(\begin{array}{ll}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{BB}_k^{}+{\tau \tau }^{\text{T}}({\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}{)}_{}^{-1}\tau \tau +({\text{Σ}\text{Σ}}_{\text{p}}^{\text{ETP}}{)}_{}^{-1}\end{array}\right)}^{-1}\left(\begin{array}{l}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\\ {BB}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k+({\text{Σ}\text{Σ}}_{\text{T}}^{\text{b}}{)}_{}^{-1}(T_{\text{b}}-\tau p{\tau p}_{\text{ETP},k})+({\text{Σ}\text{Σ}}_{\text{p}}^{\text{ETP}}{)}_{}^{-1}({pp}_{\text{ETP}(0)}-{pp}_{\text{ETP},k})\end{array}\right)$$

(23)

The sea-surface temperature has a great fluctuation throughout the year, and thereby we recommend a very loose constraint on the parameter $\mathrm{\Delta }T$.
海面温度全年波动很大，因此我们建议对 parameter $\mathrm{\Delta }T$ 进行非常宽松的约束。

It is notable that Newton-type methods are of local convergence and may seriously suffer from ill-posedness of the problem, e.g., the non-uniqueness of the nonlinear parameter $u_0$, and therefore it is recommended to be found out by a grid search method. In fact, the LS solution $\widehat{ZZ}$ of (23) can be treated as a function about the variable $u_0$ because of $\tau (u_0)=\left(T_{\text{m}}{\text{e}}^{-u/u_0}\right)$ and it is denote as $\widehat{ZZ}(u_0)$.
值得注意的是，牛顿型方法具有局部收敛性，可能会严重受到问题的病态性影响，例如非线性参数 $u_0$ 的非唯一性，因此建议通过网格搜索方法找出。事实上，（ 23） 的 LS 解 $\widehat{ZZ}$ 可以被视为关于变量 $u_0$ 的函数，因为 $\tau (u_0)=\left(T_{\text{m}}{\text{e}}^{-u/u_0}\right)$ ，它表示为 $\widehat{ZZ}(u_0)$ 。

The algorithm is given in Fig. 3.
该算法如图 3 所示。

Fig. 3  图 3

1st-level optmization algorithm
1 级优化算法

Note that $T_{\text{m}},\mathrm{\Delta }T$ are linear parameters for (2) but they are nonlinear parameters both for (3) and (8), and therefore must be iteratively solved.
请注意， $T_{\text{m}},\mathrm{\Delta }T$ （ 2） 是线性参数，但它们是 （ 3） 和 （ 8） 的非线性参数，因此必须迭代求解。

The 2nd-level optimization
第 2 级优化

Next, we will take the ESSP $c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ derived from the 1st-level optimization as a RSSP $c_0(u)$ to conduct the precise positioning. Regarding the N and E direction variations and temporal variation of the sound speed, we can express the Four Dimensional (4D) SSF $c(n,e,u,t)$ to be the form as
接下来，我们将 1st level 优化 $c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ 得出的 ESSP 作为 RSSP $c_0(u)$ 进行精确定位。关于声速的 N 和 E 方向变化和时间变化，我们可以将四维 （4D） SSF $c(n,e,u,t)$ 表示为以下形式

$$c(n,e,u,t)=c_0(u)+k_{\text{t}}(u,t)t+k_{\text{n}}(u,t)n+k_{\text{e}}(u,t)e$$

(24)

where $c_0(u):=c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ is the ESSP obtained from the 1st-level optimization, $k_{\text{n}}(u,t),k_{\text{e}}(u,t),k_{\text{t}}(u,t)$ are sound speed gradients of the SSF relative to RSSP.
其中 $c_0(u):=c_{\text{ESSP}}(u,{pp}_{\text{ETP}})$ 是从第一级优化 获得的 ESSP， $k_{\text{n}}(u,t),k_{\text{e}}(u,t),k_{\text{t}}(u,t)$ 是 SSF 相对于 RSSP 的声速梯度。

It is very familiar to the sound speed inversion in the 1st-level optimization that the spatiotemporal sound speed gradients can be jointly estimated with the seafloor geodetic coordinates. With “Appendix 1” we can directly write out the positioning model regarding the three sound speed gradients, that is
与一级优化中的声速反演非常熟悉的是，时空声速梯度可以与海底大地坐标联合估计。使用“附录 1”，我们可以直接写出关于三个声速梯度的定位模型，即

$$\begin{array}{rl}{T}_{\text{obs}}& ={\int }_{u_{\text{X}}}^{u_{\text{x}}}{(\cos \beta (u))}^{-1}{c}_0^{-1}(u)\text{d}u\\ & -m_{\text{t}}{Z}_{\text{t}}(t,p_{\text{t}})-m_{\text{n}}{Z}_{\text{n}}(t,{pp}_{\text{n}})-m_{\text{e}}{Z}_{\text{e}}(t,{pp}_{\text{e}})-m_{{\text{n}}^{\mathrm{\prime }}}{Z}_{{\text{n}}^{\mathrm{\prime }}}(t,p_{{\text{n}}^{\mathrm{\prime }}})-m_{{\text{e}}^{\mathrm{\prime }}}{Z}_{{\text{e}}^{\mathrm{\prime }}}(t,{pp}_{{\text{e}}^{\mathrm{\prime }}})+{\text{ε}}_{\text{T}}\end{array}$$

(25)

where t is the travel time observation,
其中 t 是旅行时间观测值，

$$\{\begin{array}{l}{m}_{\text{t}}=(\cos z{)}^{-1}\\ m_{\text{n}}=(\cos z{)}^{-1}\tan z\cos \alpha \\ m_{\text{e}}=(\cos z{)}^{-1}\tan z\sin \alpha \\ m_{{\text{n}}^{\mathrm{\prime }}}=(\cos z{)}^{-1}\tan z^{\mathrm{\prime }}\cos {\alpha }^{\mathrm{\prime }}\\ m_{{\text{e}}^{\mathrm{\prime }}}=(\cos z{)}^{-1}\tan z^{\mathrm{\prime }}\sin {\alpha }^{\mathrm{\prime }}\end{array}$$

(26)

is the mapping function, $z$ and $\alpha$ are the zenith angle and azimuth angle of the sea-surface platform observed at the seafloor geodetic station, respectively; $z^{\mathrm{\prime }}$ and ${\alpha }^{\mathrm{\prime }}$ are the zenith angle and azimuth angle observed at the seafloor geodetic network array center, respectively.
是制图函数， $z$ $\alpha$ 分别是在海底大地测量站观测到的海面平台的天顶角和方位角; $z^{\mathrm{\prime }}$ 和 ${\alpha }^{\mathrm{\prime }}$ 分别是在海底大地测量网络阵列中心观测到的天顶角和方位角。

$$\{\begin{array}{l}{Z}_{\text{t}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{t{k}_{\text{t}}(u,t)}{c_0^2(u)}\text{d}u\\ Z_{\text{n}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{n}}(u,t)u}{c_0^2(u)}\text{d}u\\ Z_{\text{e}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{e}}(u,t)u}{c_0^2(u)}\text{d}u\\ Z_{{\text{n}}^{\mathrm{\prime }}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{n}}(u,t)u_{\text{X}}}{c_0^2(u)}\text{d}u\\ Z_{{\text{e}}^{\mathrm{\prime }}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{e}}(u,t)u_{\text{X}}}{c_0^2(u)}\text{d}u\end{array}$$

(27)

are the five zenith delays, $\lambda (u)=(\cos z{)}^{-1}\cos (\beta (u))$ can be defined as the measure of the ray bending degree. To characterize the five time-varying zenith delay parameters $Z_K(t),K\in \left\{t,n,e,n^{\mathrm{\prime }},e^{\mathrm{\prime }}\right\}$, the following B-splines
是 5 个 Zenith 延迟， $\lambda (u)=(\cos z{)}^{-1}\cos (\beta (u))$ 可以定义为光线弯曲度的度量。为了表征 5 个时变天顶延迟参数 $Z_K(t),K\in \left\{t,n,e,n^{\mathrm{\prime }},e^{\mathrm{\prime }}\right\}$ ，以下 B 样条曲线

$$Z_K(t,{pp}_K)=\sum _{q=0}^{Q_K}{p}_{K,q}{\text{Φ}}_{K,q,k}(t)K\in \{t,n,e,n^{\mathrm{\prime }},e^{\mathrm{\prime }}\}$$

(28)

are recommended, ${pp}_K=\left(p_{K,1},p_{K,2},\dots ,p_{K,Q(K)}\right)$ is the unknown model parameter vector to be estimated, ${pp}_{K,q}$ is the coefficient of the qth k-degree B-spine basis function ${\mathrm{\Phi }}_{K,q,k}(t)$. This indicates the observation time span is split into $L_J=(Q_J-k+1)$ intervals of which the knots span is $S_K=S_{\text{all}}/L_J$ where $S_{\text{all}}$ is the observation time span. If it is not especially specified, the subscript $k$ will be omitted in the following discussion.
是要 ${pp}_K=\left(p_{K,1},p_{K,2},\dots ,p_{K,Q(K)}\right)$ 估计的未知模型参数向量， ${pp}_{K,q}$ 是第 q k 度 B-spine 基函数 ${\mathrm{\Phi }}_{K,q,k}(t)$ 的系数。这表示观测时间跨度被划分为 $L_J=(Q_J-k+1)$ 多个间隔，其中 knots 跨度是 $S_K=S_{\text{all}}/L_J$ $S_{\text{all}}$ 其中 是观测时间跨度。如果没有特别指定，在下面的讨论中将省略下标 $k$ 。

At this stage, considering the smooth nature of the physical ocean process signal, we can conduct a series of 2-order smooths on the sound speed variations, and therefore we can alternatively use the optimization criterion as follow
在这个阶段，考虑到物理海洋过程信号的平滑性质，我们可以对声速变化进行一系列的 2 阶平滑，因此我们可以交替使用优化准则，如下所示

$$\underset{XX,p{p}_K}{min}g(XX,{pp}_K):={VV}^{\text{T}}(XX,{pp}_K)PVPV(XX,{pp}_K)+\sum _{K\in \{t,n,e,n^{\mathrm{\prime }},e^{\mathrm{\prime }}\}}{\Vert Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)\Vert }_{}^2$$

(29)

where $Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)$ is the second-order derivative of $Z_K(t,{pp}_K)$ about the time,
其中 $Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)$ 是大约时间的 $Z_K(t,{pp}_K)$ 二阶导数，

$${\Vert Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)\Vert }_{}^2={\lambda }_K^2{\int }_{t_1}^{t_2}(Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K){)}^2\text{d}t$$

(30)

is the squared norm of the estimated signal $Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)$, where $t_1$ and $t_2$ are the starting time and end time of the observation, respectively, which is defined by the inner product $⟨f_i(t),f_j(t)⟩=\int f_i(t){\lambda }_K{f}_j(t)\text{d}t$ of $f_i(t)$ and $f_j(t)$, ${\lambda }_K^2$ is the scale factor of the space specified by the above inner produce, which can be defined as
是估计信号 $Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)$ 的平方范数 ，其中 $t_1$ 和 $t_2$ 分别是观测的开始时间和结束时间，由 $f_i(t)$ 和 $f_j(t)$ 的内积定义， ${\lambda }_K^2$ 是上述内积 $⟨f_i(t),f_j(t)⟩=\int f_i(t){\lambda }_K{f}_j(t)\text{d}t$ 指定的空间的比例因子，可以定义为

$${\lambda }_K^2={\text{σ}}_{0,K}^{-2}$$

(31)

where ${\text{σ}}_{0,K}^2$ represents the variance of the random signal $Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)$.
其中 ${\text{σ}}_{0,K}^2$ 表示随机信号 $Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)$ 的方差 。

To connect the minimization (29) with the normal form (20) of LS with constraints, we can rewrite (30) into the form as
要将 LS 的最小化 （ 29） 与带有约束的 LS 的范式 （ 20） 连接起来，我们可以将 （ 30） 改写为

$${\Vert Z_K^{\mathrm{\prime }\mathrm{\prime }}(t,{pp}_K)\Vert }_{}^2=({\mathbf{0}}^{\text{T}}-{pp}_K^{\text{T}})({\text{Σ}}_{\text{p}}^K{)}_{}^{-1}(\mathbf{0}-{pp}_K^{})$$

(32)

where  哪里

$$({\text{Σ}\text{Σ}}_{\text{p}}^K{)}_{}^{-1}={PP}_{\text{p}}^K{\text{σ}}_{0,K}^{-2}$$

(33)

of which  其中

$${PP}_{\text{p}}^K=\left(\begin{array}{cccc}\left({\text{Φ}}_{K,1}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,1}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)& \left({\text{Φ}}_{K,1}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,2}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)& \cdots & \left({\text{Φ}}_{K,1}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,Q(K)}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)\\ \left({\text{Φ}}_{K,2}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,1}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)& \left({\text{Φ}}_{K,2}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,2}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)& \cdots & \left({\text{Φ}}_{K,2}^{\mathrm{\prime }\mathrm{\prime }}(u),{\text{Φ}}_{K,Q(K)}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)\\ ⋮& ⋮& & ⋮\\ \left({\text{Φ}}_{K,Q(K)}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,1}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)& \left({\text{Φ}}_{K,Q(K)}^{\mathrm{\prime }\mathrm{\prime }}(t),{\text{Φ}}_{K,2}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)& \cdots & \left({\text{Φ}}_{K,Q(K)}^{\mathrm{\prime }\mathrm{\prime }}(u),{\text{Φ}}_{K,Q(K)}^{\mathrm{\prime }\mathrm{\prime }}(t)\right)\end{array}\right)$$

(34)

is the weight matrix. Then, Omitting the deduction, we can immediately write out the Gauss–Newton solution of (29). To save the space, we write the Gauss–Newton solution only with the first three zenith delays, that is
是权重矩阵。然后，省略推导，我们可以立即写出 （ 29） 的高斯-牛顿解。为了节省空间，我们只编写前三个天顶延迟的高斯-牛顿解，即

$$\left(\begin{array}{c}{XX}_{k+1}\\ {pp}_{\text{N},k+1}\\ {pp}_{\text{E},k+1}\\ {pp}_{\text{t},k+1}\end{array}\right)=\left(\begin{array}{c}{XX}_k\\ {pp}_{\text{N},k}\\ {pp}_{\text{E},k}\\ {pp}_{\text{t},k}\end{array}\right)+{\left(\begin{array}{cccc}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{N}}^{}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{E}}^{}& {AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{t}}^{}\\ {MM}_{\text{N}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {MM}_{\text{N}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_N^{}+({\text{Σ}\text{Σ}}_{\text{p}}^{\text{N}}{)}_{}^{-1}& {MM}_{\text{N}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{E}}^{}& {MM}_{\text{N}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{t}}^{}\\ {MM}_{\text{E}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {MM}_{\text{E}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{N}}^{}& {MM}_{\text{E}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{E}}^{}+({\text{Σ}\text{Σ}}_{\text{p}}^{\text{E}}{)}_{}^{-1}& M_{\text{E}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{t}}^{}\\ {MM}_{\text{t}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{AA}_k^{}& {MM}_{\text{t}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{N}}^{}& {MM}_{\text{t}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{E}}^{}& {MM}_{\text{t}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{MM}_{\text{t}}^{}+({\text{Σ}\text{Σ}}_{\text{p}}^{\text{t}}{)}_{}^{-1}\end{array}\right)}^{-1}\left(\begin{array}{c}{AA}_k^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k\\ \begin{array}{c}{MM}_{\text{N}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k-({\text{Σ}\text{Σ}}_{\text{p}}^{\text{N}}{)}_{}^{-1}{pp}_{\text{N},k}\\ {MM}_{\text{E}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k-({\text{Σ}\text{Σ}}_{\text{p}}^{\text{E}}{)}_{}^{-1}{pp}_{\text{E},k}\\ {MM}_{\text{t}}^{\text{T}}{\text{Σ}\text{Σ}}_{\text{L}}^{-1}{VV}_k-({\text{Σ}\text{Σ}}_{\text{p}}^{\text{t}}{)}_{}^{-1}{pp}_{\text{t},k}\end{array}\end{array}\right)$$

(35)

where ${MM}_K^{}$ is the design matrix of the parameter ${pp}_K$. The above formula can be easily extended for estimating the five zenith delays. Note that hyperparameter ${\text{σ}}_{0,K}^2$ might be optimally determined by the cross-test method or AIC criterion, but hereafter we set them to be zeros.
其中 ${MM}_K^{}$ 是参数 ${pp}_K$ 的设计矩阵 。上述公式可以很容易地扩展以估计五个天顶延迟。请注意，超参数 ${\text{σ}}_{0,K}^2$ 可能由交叉测试方法或 AIC 标准最佳确定，但此后我们将它们设置为零。

Experiment results and tests
实验结果和测试

Experimental data  实验数据

The test adopts Japanese seafloor geodetic network array MYGI to verify the proposed GNSS-A positioning models. The observation span of the opened MYGI station data is from 2011 to 2020, having 35 repeated observations. The long-term displacement time series solution facilitates the positioning accuracy verification based on the fact the station has a linear motion trend with the tectonic movement of the plate (Fig. 4).
该测试采用日本海底大地测量网络阵列 MYGI 对所提出的 GNSS-A 定位模型进行了验证。已打开的 MYGI 站数据的观测跨度为 2011 年至 2020 年，共 35 次重复观测。长期位移时间序列解决方案有助于根据站与板块构造运动呈线性运动趋势的事实进行定位精度验证（图 4）。

Fig. 4  图 4

MYGI station location  MYGI 车站位置

We take the GNSS-Acoustic Ranging combined POsitioning Solver (GARPOS V1.0.0) software output as an external reference solution to evaluate the propose models. The GARPOS software adopts the recommended default parameter settings.
我们采用 GNSS-Acoustic Ranging 组合定位求解器 （GARPOS V1.0.0） 软件输出作为外部参考解决方案来评估所提出的模型。GARPOS 软件采用推荐的默认参数设置。

1st-level optimization results
第一级优化结果

Prior temperatures from Argo observation
Argo 观测的先前温度

We collected Argo observations around MYGI within an area of 110,889 km², lying within the interval [(142.9167° ± 1.5°) E, (38.0833° ± 1.5°) N], from March 28, 2011 to June 15, 2020, see Fig. 5. It shows that the sea-surface temperature has a very large fluctuation, but the seawater bottom temperature is very stable. It is also hard to precisely fix the thermocline depth which lies within the interval (300 m 500 m).
我们收集了 2011 ² 年 3 月 28 日至 2020 年 6 月 15 日期间 [（142.9167° ± 1.5°） E， （38.0833° ± 1.5°） N] 区间内的 MYGI 周围 Argo 观测数据，见图 5。它表明海面温度有很大的波动，但海水底部温度非常稳定。也很难精确确定位于间隔（300 m 500 m）内的温跃层深度。

Fig. 5  图 5

MYGI surrounding Argo observations from 2011 to 2020
2011 年至 2020 年 Argo 观测的 MYGI

The depth of MYGI station plots about 1727.8 m. the probability distribution of seawater bottom temperature at this depth in Fig. 6. It shows that, the mean of the seawater bottom temperature is 2.15 ${}^{\circ }\mathrm{C}$ and STD is 0.09 ${}^{\circ }\mathrm{C}$. Therefore, we adopt the 1st-level optimization parameter settings in Table 2. The same method is used to analyze the seawater bottom sound speed.
MYGI 站的深度约为 1727.8 m。图 6 中该深度的海水底部温度概率分布。它表明，海水底部温度的平均值为 2.15，STD ${}^{\circ }\mathrm{C}$ 为 0.09 ${}^{\circ }\mathrm{C}$ 。因此，我们采用表 2 中的第一级优化参数设置。同样的方法用于分析海水海底声速。

Fig. 6  图 6

Seawater temperature distribution at the depth 1727.8 m
深度 1727.8 m 处的海水温度分布

Table 2 Parameter settings of 1st-level optimization on ESSP
表 2 ESSP 一级优化参数设置

The parameter settings of 1st-level optimization on bilinear SSP are given in Table 3. We impose lose constraints on the surface sound speed and the two sound speed gradients. Like constraining the bottom water temperature of ESSP, the bottom sound speed is imposed with a tight constraint 1486.62 m/s with the STD 0.01 m/s.
表 3 给出了双线性 SSP 上一级优化的参数设置。我们对表面声速和两个声速梯度施加了损失约束。与限制 ESSP 的底部水温一样，底部声速受到 1486.62 m/s 的严格约束，STD 为 0.01 m/s。

Table 3 Parameter settings of 1st-level optimization on bilinear SSP
表 3 双线性 SSP 一级优化参数设置

SSP inversion precision analysis
SSP 反演精度分析

The SSP inversion results are compared with the in-field SSP and they are given in Fig. 7. It shows that the SSP inversions can overall characterize the shape of the in-field SSP, but it is hard to precisely fit the sea-surface sound speed. To precisely reflect the SSP inversion precision we adopt three piece-wise statistics, see Fig. 8.
将 SSP 反演结果与场内 SSP 进行了比较，如图 7 所示。结果表明，SSP 反演可以总体表征场内 SSP 的形状，但很难精确拟合海面声速。为了精确反映 SSP 反演精度，我们采用三分段统计，见图 8。

Fig. 7  图 7

SSP inversion results from 1st-level optimization
第一级优化的 SSP 反转结果

Fig. 8  图 8

Piece-wise statistics for SSP inversion precision
SSP 反转精度的分段统计

Figure 8 shows that both for the bilinear SSP and ESSP the sound speed inversion precision increases with the depth, e.g., for the proposed ESSP, the mean bias and STD of the inversion sound speed in deepwater are 0.085 m/s and 2.448 m/s respectively, but those in the shallow water are up to 2.455 m/s and 10.862 m/s respectively. Table 4 shows that the inversion precision of proposed ESSP is overall better than that of the bilinear SSP even though the bilinear SSP has a relatively small bias for whole water column. Note that it is hard to avoid a bias in the inversion SSP relative to the in-field SSP.
图 8 表明，对于双线性 SSP 和 ESSP 的声速反演精度都随着深度的增加而增加，例如，对于所提出的 ESSP，深水反演声速的平均偏差和 STD 分别为 0.085 m/s 和 2.448 m/s，而浅水的反演声速分别高达 2.455 m/s 和 10.862 m/s。表 4 显示，所提出的 ESSP 的反演精度总体上优于双线性 SSP，尽管双线性 SSP 对整个水柱的偏差相对较小。请注意，很难避免反转 SSP 相对于场内 SSP 的偏差。

Table 4 Mean bias and STD of the inversion SSP
表 4 反转 SSP 的平均偏差和 STD

Positioning precision analysis
定位精度分析

We adopt the Argo SSP nearby MYGI station and the inversion SSPs to perform the positioning. Taking the positioning result based on the in-field SSP as a reference, we can then figure out the positioning errors caused by the substitutions of the in-field SSP and they given in Fig. 9. It shows that the Argo SSP solution has a relatively large positioning error especially in the vertical direction, but the solutions based on the inversion SSPs are almost same with each other.
我们采用 MYGI 站附近的 Argo SSP 和反转 SSP 进行定位。以基于场内 SSP 的定位结果为参考，我们可以计算出场内 SSP 替换引起的定位误差，如图 9 所示。结果表明，Argo SSP 解具有比较大的定位误差，尤其是在垂直方向上，但基于反演 SSP 的解彼此几乎相同。

Fig. 9  图 9

Solutions of Argo SSP, bilinear SSP and ESSP
Argo SSP、双线性 SSP 和 ESSP 的解决方案

The bias and STD of the positioning error series in Fig. 9 are figured out and they are given in Table 5. It shows that both the proposed ESSP solution and bilinear SSP can achieve decimeter-level-precision positioning, e.g., the positional error of the proposed ESSP solution for each coordinate component doesn’t exceed 0.4 m, but there are biases 0.003 m and − 0.016 m in E direction and N direction, respectively. The vertical bias of proposed ESSP solution is smaller than that of the bilinear SSP solution and they are − 0.219 m and − 0.260 m, respectively, but their horizontal biases can be considered at the same order for centimeter-precision-level positioning. Although the bilinear SSP inversion precision is significantly lower than that the proposed ESSP, both the bilinear SSP solution and the proposed ESSP solution possess almost the same residual with the in-field SSP solution, see Fig. 10. This indicates that the two inversion SSPs are approximately equivalent to each other for the geodetic positioning application. This also shows that imposing a prior knowledge about the sound speed on the inversion is vitally important to obtain a meaningful SSP inversion solution.
图 9 中定位误差序列的偏置和 STD 已计算出来，它们在表 5 中给出。结果表明，所提出的 ESSP 解和双线性 SSP 都可以实现分米级精度定位，例如，所提出的 ESSP 解对每个坐标分量的位置误差不超过 0.4 m，但在 E 方向和 N 方向分别存在 0.003 m 和 − 0.016 m 的偏差。所提出的 ESSP 解决方案的垂直偏差小于双线性 SSP 解决方案的垂直偏差，它们分别为 -0.219 m 和 -0.260 m，但它们的水平偏差可以以相同的顺序进行厘米精度级定位。尽管双线性 SSP 反演精度明显低于所提出的 ESSP，但双线性 SSP 解和所提出的 ESSP 解都具有与场内 SSP 解几乎相同的残差，见图 10。这表明对于大地测量定位应用程序，两个反演 SSP 彼此大致等效。这也表明，在逆温上施加有关声速的先验知识对于获得有意义的 SSP 逆温解至关重要。

Table 5 Positioning error comparison of the solutions based on different SSPs
表 5 基于不同 SSP 的解的定位误差比较

Fig. 10  图 10

1st-level optimization residuals of different SSP solutions
不同 SSP 解的 1 级优化残差

2nd-level optimization results
第 2 级优化结果

We use the parameter settings in Table 6 for performing the 2nd-level optimization. Note that the positioning accuracy may be further improved by optimally selecting the knots span $S_{\text{t}},S_{\text{NEU}}$ and the smooth factor ${\text{σ}}_K^{-2}$.
我们使用表 6 中的参数设置来执行第二级优化。请注意，通过最佳选择节长 $S_{\text{t}},S_{\text{NEU}}$ 和平滑因子 ${\text{σ}}_K^{-2}$ ，可以进一步提高定位精度。

Table 6 Parameter settings for 2nd-level optimization
表 6 二级优化参数设置

The above empirical parameter settings listed in Table 6 will be further validated in following page.
表 6 中列出的上述经验参数设置将在下一页中进一步验证。

The positioning results based on the 2nd-level optimization are given in Fig. 11. It shows that the E and N coordinate biases existing in the 1st-level optimization have been removed at a great extent.
基于 2 级优化的定位结果如图 11 所示。这表明 1 级优化中存在的 E 和 N 坐标偏差在很大程度上被去除了。

Fig. 11  图 11

Comparison of 1st-level and 2nd-level optimizations
第 1 级和第 2 级优化的比较

Table 7 shows that the horizontal precision of 2nd-level optimization for ESSP is better than 3 mm, while the vertical precision is better than 3 cm. The horizontal positioning accuracy of ESSP is very close to that of bilinear SSP, but vertical STD of ESSP solution is 0.0238 m which is significantly smaller than that 0.0433 m of bilinear SSP solution. Further considering the positioning error of the reference solution based on the in-field SSP, we can draw the conclusion that the proposed two-level optimization approach can achieve almost the same horizontal positioning precision with that based on in-field SSP. This conclusion becomes more solid when making a comparison among the 2nd-level optimization residuals of different SSP solutions as shown in Fig. 12. It also shows that the residuals are sharply shrunk after applying the 2nd optimization, compare to Fig. 9. A special attention on an existing bias about 1 cm in the vertical direction needs to be paid in future studies.
表 7 显示，ESSP 的二级优化的水平精度优于 3 mm，而垂直精度优于 3 cm。ESSP 的水平定位精度与双线性 SSP 非常接近，但 ESSP 解的垂直 STD 为 0.0238 m，明显小于双线性 SSP 解的 0.0433 m。进一步考虑基于场内 SSP 的参考解的定位误差，我们可以得出结论，所提出的两级优化方法可以实现与基于场内 SSP 的水平定位精度几乎相同的水平定位精度。当对不同 SSP 解决方案的 2 级优化残差进行比较时，这个结论变得更加可靠，如图 12 所示。它还显示，与图 9 相比，在应用第二次优化后，残差急剧缩小。在未来的研究中，需要特别注意垂直方向上约 1 cm 的现有偏差。

Table 7 Accuracy comparison of 1st-level optimization and 2nd-level optimization
表 7 一级优化和二级优化精度对比

Fig. 12  图 12

2nd-level optimization residuals of different SSP solutions
不同 SSP 解的 2 级优化残差

Long-term displacement time series analysis
长期位移时间序列分析

Let ${\overline{XX}}_j$ be the undetermined array geometry to perform the rigid-array fixed solution, $\mathrm{\Delta }{XX}^{\left(\kappa \right)}$ be the positional difference of the array center for kth-epoch observation, and they can be determined by solving the equation as (Watanabe et al., 2020):
设 ${\overline{XX}}_j$ 为执行刚阵列固定解的未确定数组几何结构， $\mathrm{\Delta }{XX}^{\left(\kappa \right)}$ 为 k 个纪元观测的数组中心的位置差，它们可以通过求解方程来确定，如下所示（Watanabe et al.， 2020）：

$$\{\begin{array}{ll}{XX}_j^{\left(\kappa \right)}={\delta }_j^{\left(\kappa \right)}{\overline{XX}}_j+{\delta }_j^{\left(\kappa \right)}\mathrm{\Delta }{XX}^{\left(\kappa \right)}& \left(j=1,\dots ,w\text{ ; }\kappa =1,\dots ,\eta \right)\\ 0=\sum _{\kappa =1}^{\eta }\mathrm{\Delta }{XX}^{\left(\kappa \right)}& \end{array}$$

(36)

where, if the transponder j is used in $\kappa$th observation, ${\delta }_j^{\left(\kappa \right)}=1$. In addition, ${\delta }_j^{\left(\kappa \right)}=0$. where w and $\eta$ are the number of transponders and epochs, respectively, and ${XX}_j^{\left(\kappa \right)}$ denotes the transponders’ position at the kth-epoch.
其中，如果应答器 j 用于 $\kappa$ th 观测，则 ${\delta }_j^{\left(\kappa \right)}=1$ 。此外， ${\delta }_j^{\left(\kappa \right)}=0$ .其中 w 和 $\eta$ 分别是应答器和纪元的数量，表示 ${XX}_j^{\left(\kappa \right)}$ 应答器在第 k 纪元的位置。

The long-term displacement time series obtained by (36) is plotted in Fig. 13. The GARPOS solution time series based on in-field SSP as an external reference is also plotted in Fig. 13 for conducting the comparison. Then, we can use the linear model x(t) = x₀ + vt + e where x is the displacement time series and v is the velocity to fit the displacement time series to verify the proposed approach. The LS fitting residual as an estimation of the observation error e can be then used to evaluate the positioning precision.
由 （ 36） 获得的长期位移时间序列绘制在图 13 中。图 13 中还绘制了基于现场 SSP 作为外部参考的 GARPOS 解决方案时间序列，以便进行比较。然后，我们可以使用线性模型 x（t） = x ₀ + vt + e，其中 x 是位移时间序列，v 是拟合位移时间序列的速度来验证所提出的方法。然后，可以使用 LS 拟合残差作为观测误差 e 的估计值来评估定位精度。

Fig. 13  图 13

GARPOS solution and the proposed 2nd-level optimization solution
GARPOS 解决方案和建议的第 2 级优化解决方案

Figure 13 shows that the proposed 2nd-level optimization approach can produce almost the same station movement trend as the GARPOS solutions, and more detailed information about the time series residuals is given in Table 8.
图 13 显示，所提出的 2 级优化方法可以产生与 GARPOS 解决方案几乎相同的台站移动趋势，表 8 中给出了有关时间序列残差的更多详细信息。

Table 8 Displacement time series analysis results
表 8 位移时间序列分析结果

Table 8 shows that the largest difference among the station velocities along E, N, U directions is about 5.5 mm/a and it happens in the U direction. The horizontal displacement time series residual STDs of proposed models are close to that of GARPOS model. The influence of the substitution of the in-field SSP with the self-structured ESSP on the station velocity estimation are 0.0, 0.4 and − 1.4 mm/a in E, N, U directions, respectively, which can be ignored in applying GNSS-A to the seafloor geodesy in the current state of the art. We can draw almost the same conclusion for apply the bilinear SSP to the seafloor geodetic positioning. It is a very interesting but important that for seafloor geodetic positioning at centimeter precision and long-term tectonic displacement monitoring in-field SSP measurements might be unnecessary. However, as an extra product the SSP inversion should keep its physical meaning by imposing a prior knowledge, such as the hydrostatic pressure of the water and the water temperature’s exponential decaying characteristic with the depth.
表 8 显示，沿 E、N、U 方向的站点速度之间的最大差异约为 5.5 mm/a，并且发生在 U 方向。所提模型的水平位移时间序列残差 STD 与 GARPOS 模型的水平位移时间序列残差 STD 接近。用自结构 ESSP 替换场内 SSP 对台站速度估计的影响在 E、N、U 方向分别为 0.0、0.4 和 − 1.4 mm/a，在当前技术水平下将 GNSS-A 应用于海底大地测量时可以忽略不计。对于将双线性 SSP 应用于海底大地测量定位，我们可以得出几乎相同的结论。对于厘米级精度的海底大地定位和长期的构造位移监测，现场 SSP 测量可能没有必要，这是一个非常有趣但很重要的事情。然而，作为一个额外的产品，SSP 反演应该通过施加先验知识来保持其物理意义，例如水的静水压力和水温随深度的指数衰减特性。

Remarks and conclusions  备注和结论

High-cost in-field SSP measurements have prevented GNSS-A from global seafloor geodesy, especially for real-time applications. The proposed self-structured SSP (SSSP) approach is a useful way to facilitate GNSS-A for conducting the large-scale and even global seafloor geodetic positioning.
高成本的现场 SSP 测量使 GNSS-A 无法进行全球海底大地测量，尤其是对于实时应用。所提出的自结构 SSP （SSSP） 方法是促进 GNSS-A 进行大规模甚至全球海底大地定位的有用方法。

The overall shape of the in-field SSP can be characterized by the proposed three-parameter empirical temperature profile by performing a joint estimation of the three parameters with the seafloor geodetic coordinates. However, the global optimal solution of the thermocline depth parameter is hard to be obtained in the context of the Gauss–Newton method because of its non-uniqueness and local convergence, and therefore the grid search method is recommended to be used. The seawater bottom temperature might also face with ill-posed problems, but fortunately the seawater bottom temperature prior constraint can be easily structured by the history ocean environment observations or by the current ocean temperature knowledge because of its stability.
通过与海底大地坐标对三个参数进行联合估计，可以通过提出的三参数经验温度剖面来表征场内 SSP 的整体形状。然而，由于温跃层深度参数的非唯一性和局部收敛性，在高斯-牛顿方法的背景下很难获得温跃层深度参数的全局最优解，因此建议使用网格搜索方法。海水海底温度也可能面临病态问题，但幸运的是，由于其稳定性，可以很容易地通过历史海洋环境观测或当前海洋温度知识来构建海水海底温度先验约束。

The inaccuracy of the self-structured SSP can be almost completely absorbed by the proposed 2nd-level optimization such that it can achieve almost the same positioning result as that based on the in-field SSP. The influence of substituting the in-field SSP with the proposed SSSP on the horizontal positioning is less than 3 mm while that on the vertical positioning is better than 3cm in the STD sense. The influence of the substitution of the in-field SSP with the self-structured SSP on the station velocity estimation are further reduced to be omitted for applying GNSS-A to the seafloor geodesy in current state of the art. Sound speed inversion accuracy of the proposed SSSP is more accurate than the bilinear SSP and this leads to a more accurate vertical positioning precision.
自结构 SSP 的不精确性几乎可以完全被所提出的第二级优化所吸收，从而可以获得与基于场内 SSP 几乎相同的定位结果。用建议的 SSSP 代替场内 SSP 对水平定位的影响小于 3 mm，而对垂直定位的影响在 STD 意义上优于 3 cm。在当前技术水平下，将 GNSS-A 应用于海底大地测量时，用自结构 SSP 替换场内 SSP 对台站速度估计的影响进一步减少，从而省略。所提出的 SSSP 的声速反转精度比双线性 SSP 更准确，这导致了更准确的垂直定位精度。

Appendix 1: Seafloor geodetic positioning models with five acoustic zenith delays

As shown in Fig. 14, without loss of generality, let the origin of NEU coordinate system be the centroid of the seafloor geodetic network array, let the RSSP $c_0(u)$ be at $(n_0=0,e_0=0,t_0=0)$, the first-order Taylor series expansion of SSF $c(n,e,u,t)$ at $(n_0=0,e_0=0,t_0=0)$ reads

$$c(n,e,u,t)=c_0(u)+k_{\text{n}}(u,t)n+k_{\text{e}}(u,t)e+k_{\text{t}}(u,t)t$$

(37)

where $k_{\text{n}}(u,t):=k_{\text{n}}(0,0,u,t)$, $k_{\text{e}}(u,t):=k_{\text{e}}(0,0,u,t)$, $k_{\text{t}}(u,t):=k_{\text{t}}(0,0,u,t)$ are the partial derivatives about the variable n, e and t, respectively, which represent the gradients of SSF.

Fig. 14

Seafloor geodetic network geometry

As shown in Fig. 14, let ${\alpha }^{\mathrm{\prime }}$ and $z^{\mathrm{\prime }}$ are the azimuth angle and zenith angle of the surface transducer position x, $XX$ is the seafloor transponder coordinate vector to be estimated, $u_{\text{X}}$ be the depth of the seafloor geodetic station, then SSP at the position x can be expressed as

$$\begin{array}{rl}{c}_0^{\mathrm{\prime }}(u,t)& =c_0(u)+k_{\text{t}}(u,t)t+n_{\text{x}}{k}_{\text{n}}(u,t)+e_{\text{x}}{k}_{\text{e}}(u,t)\\ & =c_0(u)+k_{\text{t}}(u,t)t+u_{\text{X}}\tan (z^{\mathrm{\prime }})(\cos ({\alpha }^{\mathrm{\prime }})k_{\text{n}}(u,t)+\sin ({\alpha }^{\mathrm{\prime }})k_{\text{e}}(u,t))\end{array}$$

(38)

The acoustic range signal travelling time is determined by the SSP along the ray. With SSP $c_0^{\mathrm{\prime }}(u,t)$ at the surface transducer position, we can figure out SSP along the ray. The horizontal gradient of SSF is very small and therefore the sound speed along the ray can be approximated by that along the Line of Sight (LOS). Then, let $\alpha$ and $z$ are the azimuth angle and zenith angle of the ray observed at the seafloor transponder position x, at arbitrary depth u, we have

$$c(LOS)=c_0^{\mathrm{\prime }}(u,t)+u\tan (z)(\cos (\alpha )k_{\text{n}}(u,t)+\sin (\alpha )k_{\text{e}}(u,t))$$

(39)

For approximating the ray travelling time function $T={\int }_L{c}^{-1}(L)\text{d}s$ where L is the ray trajectory, with the Taylor series approximation $(y+\text{d}y{)}^{-1}=y^{-1}-y^{-2}\text{d}y$ we have

$$\begin{array}{rl}T& ={\int }_L{c}^{-1}(L)\text{d}s\approx {\int }_L{c}^{-1}(LOS)\text{d}s\\ & ={\int }_{u_{\text{X}}}^{u_{\text{x}}}(\cos \beta (u){)}^{-1}{c}^{-1}(LOS)\text{d}u\\ & \approx {\int }_{u_{\text{X}}}^{u_{\text{x}}}(\cos \beta (u){)}^{-1}{c}_0^{-1}(u)\text{d}u-D_{\text{t}}(t)-\tan (z)(\cos (\alpha )D_{\text{n}}(t)+\sin (\alpha )D_{\text{e}}(t))\\ & -\tan (z^{\mathrm{\prime }})(\cos ({\alpha }^{\mathrm{\prime }})D_{{\text{n}}^{\mathrm{\prime }}}(t)+\sin ({\alpha }^{\mathrm{\prime }})D_{{\text{e}}^{\mathrm{\prime }}}(t))\end{array}$$

(40)

where

$$\{\begin{array}{c}{D}_{\text{t}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}\frac{1}{\cos \beta (u)}\frac{t{k}_{\text{t}}(u,t)}{c_0^2(u)}\text{d}u\\ D_{\text{n}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}\frac{1}{\cos \beta (u)}\frac{u{k}_{\text{n}}(u,t)}{c_0^2(u)}\text{d}u\\ D_{\text{e}}\text{(}t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}\frac{1}{\cos \beta (u)}\frac{u{k}_{\text{e}}(u,t)}{c_0^2(u)}\text{d}u\\ D_{{\text{n}}^{\mathrm{\prime }}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}\frac{1}{\cos \beta (u)}\frac{u_{\text{X}}{k}_{\text{n}}(u,t)}{c_0^2(u)}\text{d}u\\ D_{{\text{e}}^{\mathrm{\prime }}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}\frac{1}{\cos \beta (u)}\frac{u_{\text{X}}{k}_{\text{e}}(u,t)}{c_0^2(u)}\text{d}u\end{array}$$

(41)

are the acoustic delays along LOS. To counterbalance the ray incident angle $\beta (u)$ involved in (41), introducing the azimuthal auqengle z and then we can rewrite the model (40) into the form as follow

$$T_{\text{obs}}={\int }_{u_{\text{X}}}^{u_{\text{x}}}(\cos \beta (u){)}^{-1}{c}_0^{-1}(u)\text{d}u-m_{\text{t}}{Z}_{\text{t}}(t)-m_{\text{n}}{Z}_{\text{n}}(t)-m_{\text{e}}{Z}_{\text{e}}(t)-m_{{\text{n}}^{\mathrm{\prime }}}{Z}_{{\text{n}}^{\mathrm{\prime }}}(t)-m_{{\text{e}}^{\mathrm{\prime }}}{Z}_{{\text{e}}^{\mathrm{\prime }}}(t)+{\text{ε}}_{\text{T}}$$

(42)

where $T_{\text{obs}}$ is the observed value of round-trip time, ${\text{ε}}_{\text{T}}$ is the random error of observation.

$$\{\begin{array}{l}{m}_{\text{t}}=(\cos z{)}^{-1}\\ m_{\text{n}}=(\cos z{)}^{-1}\tan z\cos \alpha \\ m_{\text{e}}=(\cos z{)}^{-1}\tan z\sin \alpha \\ m_{{\text{n}}^{\mathrm{\prime }}}=(\cos z{)}^{-1}\tan z^{\mathrm{\prime }}\cos {\alpha }^{\mathrm{\prime }}\\ m_{{\text{e}}^{\mathrm{\prime }}}=(\cos z{)}^{-1}\tan z^{\mathrm{\prime }}\sin {\alpha }^{\mathrm{\prime }}\end{array}$$

(43)

is the zenithal mapping function,

$$\{\begin{array}{c}{Z}_{\text{t}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{t{k}_{\text{t}}(u,t)}{c_0^2(u)}\text{d}u\\ Z_{\text{n}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{n}}(u,t)u}{c_0^2(u)}\text{d}u\\ Z_{\text{e}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{e}}(u,t)u}{c_0^2(u)}\text{d}u\\ Z_{{\text{n}}^{\mathrm{\prime }}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{n}}(u,t)u_{\text{X}}}{c_0^2(u)}\text{d}u\\ Z_{{\text{e}}^{\mathrm{\prime }}}(t)={\int }_{u_{\text{X}}}^{u_{\text{x}}}{\lambda }^{-1}(u)\frac{k_{\text{e}}(u,t)u_{\text{X}}}{c_0^2(u)}\text{d}u\end{array}$$

(44)

where $\lambda (u)=(\cos z{)}^{-1}\cos \beta (u)\approx 1$ such that $Z_{\ast }(t)$ can represent the zenith delay unrelated to the time-varying of LOS. The above proposed model is different from the GARPOS model which adopts a default multiplicative compensation on the RSSP as

$$\begin{array}{rl}T(xx)& ={\int }_{u_{\text{X}}}^{u_{xx}}\frac{1}{c(u,t)}\text{d}s\\ & ={\int }_{u_{\text{X}}}^{u_{xx}}\frac{1}{\Gamma (t)c_0(u)}\text{d}s={\Gamma }^{-1}(t){\int }_{u_{\text{X}}}^{u_{xx}}\frac{1}{c_0(u)}\text{d}s={\Gamma }^{-1}(t)T_0(xx)\end{array}$$

(45)

where $\Gamma (t)$ is a time-varying factor to characterize SSP variation, $T_0(xx)$ is the travel time calculated by the RSSP, $T(xx)$ is the travel time calculated by the time-varying SSP. To deal with the multiplicative compensation, we can follow the GARPOS approach taking logarithms on both sides of the observation $T_{\text{obs}}={\Gamma }^{-1}\left(t\right)T_0\left(xx\right)+\text{ε}$, that is

$$\ln T_{\text{obs}}=\ln T_0(xx)+\ln {\Gamma }^{-1}(t)+{\text{ε}}^{\mathrm{\prime }}$$

(46)

Assuming ${\Gamma }^{-1}(t)=e^{\gamma (t)}$, we can have

$$\ln T_{\text{obs}}=\ln T_0(xx)+\gamma (t)+{\text{ε}}^{\mathrm{\prime }}$$

(47)

where ${\text{ε}}^{\mathrm{\prime }}$ is the observation error vector, $\gamma (t)$ is the correction coefficient, which can be parameterized by a B-spline about the time. Further considering the horizontal gradient of sound speed, GARPOS finally adopts

$$\gamma (t)={\gamma }_0\left(t\right)+{\gamma }_1\left(t\right)\cdot P^{\mathrm{\prime }}{P}^{\mathrm{\prime }}+{\gamma }_2\left(t\right)\cdot X^{\mathrm{\prime }}{X}^{\mathrm{\prime }}$$

(48)

where $P^{\mathrm{\prime }}{P}^{\mathrm{\prime }}$ is horizontal coordinate under the NEU coordinate system centered on specific original point which is the mean coordinates of transducer positions during the whole surveying session. $X^{\mathrm{\prime }}{X}^{\mathrm{\prime }}$ is horizontal coordinate under the NEU coordinate system centered on specific original point which is the mean coordinates of seafloor transponder array positions.