全球卫星导航系统(GNSS)可为地表和近空任意地方提供高精度定位服务，但其信号载体(电磁波)在海水中衰减严重，无法直接服务于海洋内部。声波在水体中衰减小、传播距离远，使得水声定位成为最主要的水下定位方法。为了支撑海洋内部的复杂人类活动，如工程建设、资源开发、环境监测，需利用水声定位技术，通过信号接力的方式，将GNSS维持的全球位置框架引入到海洋内部^[1-3]。虽然国外学者通过联合GNSS与水声定位在海底大地测量等方面取得了丰硕成果，但采用的是高度集成的海面平台和硬件设备，成本和技术要求较高^[4-5]。国内GNSS-声学海底定位实验则主要采用多种硬件设备现场安装，以相互配合的方式完成实验任务。

现场安装的水声定位实验中，不同硬件设备会生成不同类型的观测数据，主要包括GNSS观测数据、声学时延观测数据和惯性导航数据，分别反映测船位置、测船与海底目标点相对位置和测船姿态的时间序列信息，需依靠各自的时间标签进行匹配，进而反演出与真实情况相符的实验过程，以支撑水声定位数据处理。虽然利用GNSS技术驯服本地时钟可以完成ns级的时间同步，但实验过程中并非所有设备都能实时与GNSS进行时间同步，这会导致观测数据时间标签出现跳变和钟漂问题^[6]。Ikuta等^[7]通过对比发现，声学设备采用的姿态数据的时间标签存在约1 s/半天的漂移。Xu等^[8]利用高频PPP测量地震波，通过互相关搜索发现，惯性测量单元的观测数据时标中存在约0.01 s的偏差。

由于时标跳变会导致观测数据匹配失败，进而大幅降低水声定位的精度，因此针对观测数据无冗余的情况，需要考虑数据时标跳变的问题。Kuang等^[9]针对时标跳变问题，基于测船以固定速率沿直线航行的假设，推导出时标不匹配误差的表达方式，并利用最小二乘准则进行参数化解算。然而，当测船以固定速率进行圆形航迹或者以变速率进行线性航迹时，时标不匹配问题造成的影响便不再是固定值，而是与速度相关的变量，导致无法通过增设参数对时标偏差进行代数求解。为此，文本基于水声定位原理，提出基于参数搜索的时标跳变改正方法，并利用实测数据对提出的方法进行验证。

1 水声定位原理

GNSS定位技术已经相当成熟和完善，可通过精密单点定位或者事后差分定位等手段获取海面平台cm级至dm级的定位结果，因此本文研究聚焦于水声定位。对于水声定位部分，可分为同步模式和应答模式。同步模式直接由海底设备定时播发声信号，海面平台只负责接收，对应观测值为单程传播时间，可构建单程定位模型。应答模式由海面平台播发声信号，海底控制点收到信号后回传声信号给海面平台，对应观测值为双程传播时间，可构建双程定位模型。

单程模型与双程模型的定位解算过程一致，本文以单程模型为例进行说明。单程模型观测方程类似于GNSS非差定位中采用的观测方程，可表示为：

$ \left\{\begin{array}{l} c t=\Delta \rho_{\text {vos }}+\Delta \rho_{\text {hardware }}+\varepsilon+\rho \\ \rho=\sqrt{\left(x_{\text {trnspndr }}-x_{\text {trnscvr }}\right)^2+\left(y_{\text {trnspndr }}-y_{\text {trnscvr }}\right)^2+\left(z_{\text {trnspndr }}-z_{\text {trnscvr }}\right)^2} \end{array}\right. $

(1)

式中，t为声信号的单程传播时间；c为目标海域实验期间的声速；ρ为船载换能器与海底目标点的几何距离；Δρ_vos和Δρ_hardware分别为声速误差和硬件延迟引起的距离误差；ε为距离观测值中的随机误差；x_trnscvr、y_trnscvr和z_trnscvr为船载换能器在大地坐标系下的3个坐标分量，可基于GNSS定位结果、测船姿态数据和船载设备间的杆臂矢量进行计算；x_trnspndr、y_trnspndr和z_trnspndr为待求海底目标点的三维坐标。

对式(1)进行线性化处理，可得：

$ v=a_1 \mathrm{~d} x_{\text {trnspndr }}+a_2 \mathrm{~d} y_{\text {trnspndr }}+a_3 \mathrm{~d} z_{\text {trnspndr }}-\left(c t-\Delta \rho_{\text {vos }}-\Delta \rho_{\text {hardware }}-\rho_0\right) $

(2)

$ \rho_0=\sqrt{\left(x_{\text {trnspndr }_0}-x_{\text {trnscvr }}\right)^2+\left(y_{\text {trnspndr }_0}-y_{\text {trnscvr }}\right)^2+\left(z_{\text {trnspndr }_0}-z_{\text {trnscvr }}\right)^2} $

(3)

$ \left\{\begin{array}{l} a_1=\left.\frac{\partial \rho}{\partial x_{\text {trnspndr }}}\right|_{\boldsymbol{X}_{\text {trnspndr }_0}}=\frac{x_{\text {trnspndr }_0}-x_{\text {trnscvr }}}{\rho_0} \\ a_2=\left.\frac{\partial \rho}{\partial y_{\text {trnspndr }}}\right|_{\boldsymbol{X}_{\text {trnspndr }_0}}=\frac{y_{\text {trnspndr }_0}-y_{\text {trnscvr }}}{\rho_0} \\ a_3=\left.\frac{\partial \rho}{\partial z_{\text {trnspndr }}}\right|_{\boldsymbol{X}_{\text {trnspndr }_0}}=\frac{z_{\text {trnspndr }_0}-z_{\text {trnscvr }}}{\rho_0} \end{array}\right. $

(4)

式中，ρ₀为用待求参数初始值计算的几何距离；a₁、a₂和a₃为几何距离相对于待求参数的偏导数；X_trnspndr0为船载换能器在大地坐标系下的三维初始坐标；v为观测值残差。

当有n个历元观测值时，可组成误差方程组式(5)，其相应的矩阵形式见式(6)：

$ \left[\begin{array}{c} v_1 \\ v_2 \\ \vdots \\ v_n \end{array}\right]=\left[\begin{array}{ccc} \frac{x_{\text {trnspndr }_0}-x_{\text {trnscvr }_1}}{\rho_{01}} & \frac{y_{\text {trnspndr }_0}-y_{\text {trnscvr }_1}}{\rho_{01}} & \frac{z_{\text {trnspndr }_0}-z_{\text {trnscvr }_1}}{\rho_{01}} \\ \frac{x_{\text {trnspndr }_0}-x_{\text {trnscvr }_2}}{\rho_{02}} & \frac{y_{\text {trnspndr }_0}-y_{\text {trnscvr }_2}}{\rho_{02}} & \frac{z_{\text {trnspndr }_0}-z_{\text {trnscvr }_2}}{\rho_{02}} \\ \vdots & \vdots & \vdots \\ \frac{x_{\text {trnspndr }_0}-x_{\text {trnscvr }_n}}{\rho_{0 n}} & \frac{y_{\text {trnspndr }_0}-y_{\text {trnscvr }_n}}{\rho_{0 n}} & \frac{z_{\text {trnspndr }_0}-z_{\text {trnscvr }_n}}{\rho_{0 n}} \end{array}\right]\left[\begin{array}{c} \mathrm{d} x_{\text {trnspndr }} \\ \mathrm{d} y_{\text {trnspndr }} \\ \mathrm{d} z_{\text {trnspndr }} \end{array}\right]-\left[\begin{array}{c} l_1 \\ l_2 \\ \vdots \\ l_n \end{array}\right] $

(5)

$ \boldsymbol{V}=\boldsymbol{A} \mathrm{d} \boldsymbol{X}-\boldsymbol{L} $

(6)

基于最小二乘准则，利用式(6)解算得到待求参数的改正数：

$ \mathrm{d} \boldsymbol{X}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P A}\right)^{-1}\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{L}\right) $

(7)

最终定位结果$\hat{\boldsymbol{X}}$及其中误差$\hat{\sigma}_{\hat{X}}$可表示为：

$ \hat{\boldsymbol{X}}=\boldsymbol{X}_0+\mathrm{d} \boldsymbol{X} $

(8)

$ \hat{\sigma}_{\hat{X}}=\sqrt{\hat{\sigma}_0^2\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P A}\right)^{-1}} $

(9)

式中，$\hat{\sigma}_0^2$为单位权方差，可表示为：

$ \hat{\sigma}_0^2=\frac{\boldsymbol{V}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{V}}{n-3} $

(10)

2 参数搜索法

不同类型的数据产生于不同的硬件设备。其中，GNSS观测数据的时标由GNSS接收机直接标定，包含接收机钟差。考虑到GNSS接收机钟差绝对值被严格控制在1 ms以内，GNSS-声学定位过程中船速不大于4节(约2.06 m/s)，因此GNSS观测数据时标偏差的影响不超过2.06 mm，可以忽略。声学时延数据的时标由声学设备内置时钟标定，若声学设备的时钟未采用GNSS技术进行严格驯服，则时钟累积误差将影响声学观测值时间标签，甚至时延观测值本身的准确度。测船姿态数据由惯导设备直接导出，其自身并不带时间标签，可由计算机软件(如HYPACK)利用GNSS秒脉冲标定时间标签，或使用声学设备读取后与时延数据一起由声学设备内部时钟标定。

GNSS观测数据、声学时延数据和测船姿态数据的采样率不同。GNSS观测数据的采样率通常为1 Hz；声学时延数据的采样率根据实验区水深确定(浅海0.15~0.50 Hz、深海0.08~0.17 Hz)；姿态数据偏航角的采样率通常为1~5 Hz，姿态数据俯仰角和翻滚角的采样率可达50 Hz。利用时间标签对数据进行内插和整理，可将各类数据匹配在一套时标序列下。由于实际船速为变化的矢量，观测数据时标出现跳变时，跳变导致的匹配误差难以在定位观测方程中进行严格的数学表达，故需直接在观测值域利用冗余观测搜索时标跳变，进而对各类数据进行统一。

冯义楷等^[10]基于GNSS-声学海底定位实验中实测姿态数据的2套时标(HYPACK采用GNSS秒脉冲标定的时标、声学设备采用内部时钟标定的时标)，利用互相关搜索，成功发现声学设备时标存在具有稳定物理机理的偏差，并通过拟合方法反演出钟漂约为7 μs/s。但当各类观测数据无冗余观测时，上述观测值域搜索方法便不再适用。因此，本文假设式(2)中声速误差和硬件延迟引起的距离误差对定位结果中误差的影响为固定值，进而提出一种参数搜索方法。主要步骤为：1)选定GNSS观测数据的时标为基准时标，为其他观测数据时标建立以标准值为中心的搜索空间(如[-2 s, 2 s])；2)以时标的最小分辨率(如0.001 s)为间隔，等步长提取搜索空间中的值作为时钟偏差的测试值，对待搜索数据的时标进行整体调整；3)对时标调整后的数据进行匹配，并完成估计过程，得到定位结果及其中误差。由于单位权方差$\hat{\sigma}_0^2$为观测值残差v的函数，可用于表达观测值中未建模误差的大小，因此选择使单位权中误差最小的时标偏差作为搜索值，并利用其调整数据时标，进行后续的平差计算。虽然在实际海试情况下无法通过代数方式将时标偏差作为参数进行解算，但利用上述参数搜索方法可从数值角度确定观测数据的时标偏差。

3 时标跳变改正算例

利用2018-09-22青岛灵山岛海域(水深约为30 m)的实验数据进行时标跳变算法分析。实验期间，测船按预定的圆形航迹航行，观测数据包括GNSS观测数据、姿态数据、声信号传播时间、声速剖面数据以及杆臂矢量量测数据。其中，GNSS动态定位结果的北、东和高程方向中误差均值分别为1.1 cm、1.0 cm和2.7 cm；惯导系统POSMV测得的姿态数据偏航角、俯仰角和翻滚角中误差均值分别为0.02°、0.01°和0.01°；声学数据采集时间为UTC 02:30:00~06:40:00，时间跨度为250 min。实验期间，利用miniSVP声速剖面仪每隔约1.5 h测量一次声速剖面。

为聚焦于数据时标跳变研究，忽略声速误差和杆臂矢量误差的影响，假设观测方程的不准确性仅由数据时标偏差所致。按照圆形航迹对观测数据进行时间域划分以进行时标改正研究，共提取8个圆形轨迹，时间跨度介于3~4 min，半径介于40~70 m，航迹见图 1。其中，前4个圆形轨迹于02:30:00~02:50:00测得，后4个圆形轨迹于05:00:00~05:50:00测得，两者相隔约2 h，数据提取指标统计见表 1。按照数据的原始时标进行GNSS-声学海底定位，解算结果见表 2。从表中可以看出，前4个圆形轨迹的定位结果在各方向的中误差均约为mm级，单位权方差接近0，残差中误差为cm级；然而，对于后4个圆形轨迹，定位结果在各方向的中误差均约为dm级，单位权方差和残差中误差中均出现异常偏差，最大值分别达到0.806 m²和0.890 m，说明观测数据中存在时标跳变问题。

对各圆形轨迹进行时标偏差的参数搜索，在-2~2 s区间内以1 ms为间隔进行声学数据时标调整，利用时标调整后的观测数据进行GNSS-声学海底定位。假设观测方程中不存在时标偏差以外的误差项，则定位结果单位权中误差最小的时标偏差即为目标偏差。分别对每个圆形轨迹进行参数搜索，结果见图 2。从图中可以看出，前4个圆形轨迹的时标跳变集中在-17~-27 ms，后4个圆形轨迹的时标跳变则集中在967~990 ms，即时标在2组圆之间发生约1 s的跳变。分别对各航迹数据进行时标修正，再次进行GNSS-声学海底定位，结果见表 3。从表中可以看出，各轨迹定位结果中误差和残差中误差均有一定程度的提升，尤其是后4个圆形轨迹，中误差从dm级提升至与前4个圆形一致的mm级，单位权方差和残差中误差也提升至与前4个圆形一致的量级，验证了时标修正的正确性。因此，在无冗余观测信息情况下进行观测值域的时标偏差解算时，可用上述参数搜索法进行观测数据时标跳变的探测和修正。

4 结语

以GNSS定位结果为基准，利用水声定位技术进行信号接力，可将全球性位置基准引入到海洋内部。然而，国内GNSS-声学海底定位实验的设备集成度不够，需多种硬件设备各自运行，并在观测值域进行匹配融合，进而完成水下定位实验任务。但各设备所输出的观测数据时标难免存在时标跳变问题。通过对比设备冗余观测信息，可搜索观测数据中的时标跳变。针对无冗余观测信息的情况，本文提出基于参数搜索的时标改正方法，构建时标偏差的搜索空间，在其中等步长提取测试值，以水声定位估算的单位权方差为指标，确定观测值时标中的跳变。利用实测海试数据对提出的方法进行验证，结果表明，参数搜索方法成功探测到声学设备观测数据中约1 s的时标跳变，并通过改正数据时标将水声定位结果精度从dm级提升至mm级，残差中误差也由0.50~0.90 m降低到0.01~0.04 m。