参考声速剖面误差对主动式声呐定位影响仿真分析

肖圳, 薛树强, 韩保民, 李保金, 孙悦

1 山东理工大学建筑工程与空间信息学院, 山东淄博 2550492 中国测绘科学研究院, 北京 100830

0 引言

海底大地控制网是海底大地测量、水下导航和海底板块运动监测的重要基础设施(杨元喜等,2017).高精度水下定位是实现海底大地基准的关键,主要受GNSS/声呐组合观测系统的观测精度和参考声速误差影响(杨元喜等,2020).GNSS/声呐组合观测技术已在海底精准定位、海底地壳形变监测等方面取得重要研究成果(陈瀚等,2019;Yokota and Ishikawa,2019;乔学军等,2021).水下定位所使用的声速剖面(Sound Velocity Profile, SVP)观测仅代表特定时间和特定空间位置处的声速场环境信息,无法代表声呐观测时刻的声场环境信息,致使海洋环境在时间和空间上的声速扰动对精确定位的影响不容忽视(Yokota et al.,2019).目前,声速误差消除方法主要有两种,一种是建模参数化,另一种是差分模型解算.差分方法可以削弱观测值中的共性误差,但需要考虑差分观测值的相关性(Xue et al.,2022).

参考声速剖面误差是水下定位精度最棘手的误差源.实际信号传播时延不仅包含了几何距离信息,还包含了声速场时空环境对其扰动信息(Yokota et al.,2018).基于此,Fuijita等(2006)将观测时间窗口内平均声速变化以二次多项式或二次曲线进行建模,修正了声速观测误差,显著提高了定位精度; Ikuta等(2008)使用三次B样条时间函数对声速场时变影响导致的参考声速剖面变化进行建模,以削弱声速时间变化误差影响;Yokota等(2019)构造了从GNSS/声呐数据中提取声速梯度的反演方法,其提取的梯度效应不仅代表了广泛的声速结构,而且代表了在非定常扰动中产生的更详细的结构.笔者将上述基于声速变化的参数模型归纳为补偿参考声速剖面误差的声速剖面校正法.Kido等(2008)提出将声速变化影响的传播时延残差归一化到垂直分量,构造了天底延迟估计模型;Watanabe等(2020)构造了信号传播时间校正模型,并通过在观测方程中引入适当的统计特性来提取声速梯度结构的方法.笔者将上述基于声速变化影响的观测时间参数化模型归纳为传播时间校正法.

海底定位精度验证和定位模型优劣对比同样面临很大的挑战,仿真试验验证已成为一种行之有效的方法.国内外学者进行了不少探索和尝试.Xu等(2005)为模拟由声速结构时空变化引起的系统误差,考虑了常数项、短期内波、日潮或半日潮以及区域影响因素4种效应;朱祥娥(2008)根据双曲面定位模型,遵照理想介质中声线传播规律,研究出了深海定位仿真方法.Yusuke等(2016)为检验环境因素对定位精度的影响,专门开发了一个海底大地测量模拟器.为此,笔者结合GNSS/声呐定位原理开展声呐往返双程观测时间和时变声速场仿真研究,提出了使用经验正交函数(Empirical Orthogonal Function, EOF)仿真时变声速场,提出了采用信号传播时间追赶仿真法模拟往返双程观测时间;根据仿真数据进行定量研究参考声速剖面误差对定位的影响,测试对比声速剖面校正法和传播时间校正法这两类海底定位方法的适用性;此外,本文研究也可为海面测线优化设计和海底控制网优化设计提供参考.

1 主动式声呐定位模型及声速误差处理

为降低海底基准站观测设备的复杂性,避免复杂的海底时间同步难题,海底控制网定位往往采用主动式声呐定位模式,即海面测量船安装声呐换能器向水中发射询问声信号,海底基准站应答器接收询问声信号并发送应答声信号,船载换能器接收应答声信号并对信号进行反馈得到声呐信号往返程传播时间,进一步结合声呐测距原理和椭球交会方法以实现海底基准站定位(杨元喜等,2020).因此,水下定位除了受海底基准站声呐信号转发硬件延迟等误差影响外,海洋声速测量误差也是制约水下高精度定位的重要误差源.

1.1 主动式声呐观测模型

主动式声呐定位模型如图1所示,其观测方程可以表示为

图1 测船主动式声呐定位简略示意图

Ti=fi(P(ti+),P(ti-),Xj,Vi(u))+εi,

(1)

其中,Ti为声线传播时间观测值,f(·)为声线传播时间计算值,P(ti-)、P(ti+)为第i次观测换能器的发射位置和接收位置,Xj为海底第j个应答器,Vi(u)表示第i次观测时刻声速剖面,εi为第i次观测偶然误差.

若对m个海底控制点进行n次观测,式(1)的误差方程矩阵形式表示为

(2)

VTPV=min,

(3)

可得如下最小二乘解为

(4)

其中,P为权阵,为避免大入射角对定位的影响,可采用与入射角有关的加权函数或分段指数权函数构建权阵(Liu et al.,2020;王薪普等,2021).

理论上,声呐定位模型可使用高时空分辨率的声速场进行解算,可避免声速时空变化产生的不确定性.然而,现实测量中难以实施高频率的声速剖面观测,大多数情况下只能采取某一参考声速剖面V0(u)代替真实声速场实施定位计算.此时,可将(1)式改写为

Ti=fi(P(ti+),P(ti-),Xj,V0(u))+ρvi+εi

(5)

式中,ρvi表示第i次观测声速误差引起的误差项.

研究表明,采用固定的参考声速剖面只能达到分米级甚至米级精度的定位结果(Obana et al.,2000;王振杰等,2016;辛明真等,2020).为削弱时域声速误差对定位的影响,可采用顾及声速时域变化的声速误差修正方法(陈冠旭等,2022).

1.2 时域声速误差处理的声速剖面校正法

观测时刻声速剖面Vi(u)可以描述为随时间变化的系数αi与参考声速剖面V0(u)的乘积(Ikuta et al.,2008):

Vi(u)=αiV0(u),

(6)

其中αi为修正系数,表示声速随时间变化的函数,使用三次B样条函数对其进行建模,其关系表达式为

(7)

其中ak为第k个三次B样条基函数Φk对应的系数,基函数Φk可由de-Boor Cox递推公式获得(De Boor,1978).因此,式(1)的观测方程可改写为

Ti=fi(P(ti+),P(ti-),Xj,αiV0(u))+εi,

(8)

则其误差方程矩阵形式表示为

(9)

将B样条基函数系数a与海底点坐标位置参数X作为未知参数,使用高斯-牛顿迭代求解.

1.3 时域声速误差处理的传播时间校正法

双程观测传播时间Ti可以描述为随时间变化的校正系数与基于参考声速剖面计算的双程传播时间的乘积(Watanabe et al.,2020):

Ti=exp(-γi)·fi(P(ti+),P(ti-),Xj,V0(u)),

(10)

(11)

(12)

式中,exp(-γi)为校正系数表达式,表征观测时刻实际声速场由时空变化引起的实际传播时间与参考传播时间的差异率;α(t)为时间相关函数,可使用三次B样条建模;ΔP表示海面换能器位置P与海面图形中心位置的水平距离差值;ΔX表示海底应答器位置X与海面图形中心位置的水平距离差值,L*表示特征长度(这里取应答器的平均深度).

在这里笔者仅用三次B样条函数对时变声速场造成的传播时间误差进行建模.因此,式(10)的观测方程可改写为

log(Ti)=log(fi(P(ti+),P(ti-),Xj,V0(u)))

-γi+εi,

(13)

(14)

令yi=log(Ti),Fi=log(fi(P(ti+),P(ti-),Xj,V0(u))),则式(13)可以简化为

yi=Fi-γi+εi,

(15)

则其误差方程矩阵形式表示为

(16)

同样的,将B样条基函数系数a与海底点坐标位置参数X作为未知参数使用高斯-牛顿迭代求解.

数值测试表明,定位结果与平均声速有关(Sato and Fujita,2004),则等式(8)可表示为

Ti=αifi(P(ti+),P(ti-),Xj,V0(u))+εi,

(17)

由此可得,等式(8)中校正系数αi与等式(10)中校正系数表达式exp(-γi)具有相同性质.因此,声速剖面校正法与传播时间校正法在时域声速误差处理方面具有一定等效性(Watanabe et al.,2020;陈冠旭等,2022).

2 主动式声呐观测仿真

主动式声呐观测仿真需要解决两个问题,第一是双程观测值仿真,第二是往返程信号传播确定所需的海洋声速场环境仿真.

2.1 主动式声呐信号传播时间仿真

如图2所示,仿真观测中测船将按照观测策略中航线设计的图形进行走航,它随时间的轨迹可表示为H(t),假定发射时刻为起始时刻,则发射时刻位置P1与接收时刻位置P2可以表示为

P1=H(0),

(18)

P2=H(T1+T2),

(19)

其中,T1、T2分别表示往、返程传播时间.在海底基准点X、发送与接收时刻位置P1与P2、声速剖面都已知的情况下,往、返程传播时间又可根据声线跟踪算法G(P,X)计算得到

T1=G(P1,X),

(20)

T2=G(P2,X),

(21)

其中,声线跟踪算法G(P,X)采用p阶割线法以提高计算效率(Yang et al.,2023).整理式(18—21),可得

G(H(G(H(0),X)+T2),X)=T2.

(22)

将T2作为未知参数求解式(22)即可得到返程传播时间,进而确定接收时刻位置P2,但由于不同的观测图形H(t)表示不同,解析表达式具有很强的不确定性.

为解决这个问题,笔者使用一个类似逐次逼近方法,算法模型步骤示意如图2所示,计算流程如下:

(1)如图2a所示,根据已知的海底基准点X、发射时刻位置P1和声速剖面通过式(20)确定T1,同时假定发射时刻位置P1与接收时刻位置P2相同,通过将式(21)确定T2;

(2)将往、返程传播时间T1、T2代入到式(19)使图2a中接收位置更新至如图2b所示的接收时刻位置P2;

(3)如图2c所示,将更新后的接收时刻位置P2代入到式(21)更新T2;

上述算法的几何意义表示为不断追赶距离差的过程,首先假定发射位置与接收位置相同,即往程传播时间与返程传播时间相等,而此时对应的位置与假设位置之间存在位置差,因此,需要根据新位置重新进行传播时间计算,继续追赶,直至双程传播时间小于一定的阈值,则可确定测船的接收位置,完成主动式声呐信号传播时间和位置仿真.为了讨论方便,笔者将上述算法称为信号传播时间追赶仿真法.

2.2 时变声速误差仿真

为了刻画海洋动态变化引起的声速场时变特性,笔者建议结合测区的实测SVP数据,采用经验正交函数的方式重构声呐观测时刻的声速剖面.假设已有M条历史观测声速剖面,声速剖面矩阵V可表示为(禹小康,2021)

(23)

(24)

声速扰动矩阵ΔV可以表示为

Δvi(uj)=vi(uj)-v1(uj).

(26)

由此,可计算扰动矩阵的协方差R:

(27)

对协方差矩阵R进行特征值分解,可得特征值矩阵Λ和特征向量矩阵S:

R·S=Λ·S,

(28)

Λ=diag(λ1,λ2,…,λN),

(29)

S=[s1,s2,…,sN],

(30)

式中Λ的特征值按从大到小的顺序排列,即λ1>λ2>…>λN,si为特征值λi对应的特征向量.将声速扰动投影到特征向量上获得EOF系数矩阵C:

(31)

式中,cij为第i条声速剖面对应的第j个EOF系数.

通过已有的先验信息(例如历史浮标数据等)得知声速场是存在周期变化的,例如:日周期,半日周期等,所以采用与时间变量有关的三角函数构建EOF系数函数模型如下所示:

cij=e1,jsin(2π·ti+φ1)+e2,jsin(4π·ti+φ2),(32)

式中,为符合海洋声速场短时间内的变化规律,在这里笔者将第j个EOF系数对应模型参数e1,j、e2,j、φ1、φ2通过先验信息进行人为指定;t以天为基准单位.

(33)

式中cw(t)表示由指定的EOF系数函数模型参数带入式(32),再由式(32)根据时间变量求解的EOF系数.将仿真中的观测时间代入上式即可求得任意观测时刻声速剖面,实现海洋时变声速剖面的仿真.上述声速场模型仿真中未考虑声速场水平梯度影响,对于较大范围水下定位,声速场的水平梯度影响不容忽视,这是本文的后续研究方向.

3 仿真试验分析

3.1 仿真观测数据生成流程

考虑到水下定位的误差源主要包括(辛明真,2020):(1)随机性误差:GNSS定位误差、载体姿态测量误差、水声时延测量误差;(2)系统性误差:声速剖面代表性误差、水声系统响应误差、设备安装校准偏差.声呐观测仿真流程如图3所示,仿真数据包括实时声速剖面仿真、参考声速剖面仿真以及声速传播时间仿真等.

对于海面和海底控制网型仿真,借鉴Nakamura对测线配置的最佳尺寸相关研究,在仿真3000 m水深定位条件下,海面控制网和海底控制网尺寸均采用3000 m左右(Nakamura et al.,2021),详见图4.研究表明,将多个海底基准站的中心点作为虚拟基准点,可以将声速误差对水平位置的影响减小到几厘米左右(Sato and Fujita,2004).因此,采用多个海底基准站进行观测,同时在仿真数据的解算部分,以虚拟基准点解算坐标偏差均值(MEAN)、标准差(STD)以及均方根误差(RMSE)作为各种图形控制网和解算方法的精度评定指标.

图4 仿真试验海面海底控制网设计图

水下定位观测仿真相关参量见表1、2.

表1 观测误差仿真参数设置表

为分析走航轨迹对水下定位精度的影响,笔者设计了不同的航迹图形,如表3所示.

表3 控制网型海面航迹图形设计方案

笔者选取声速剖面观测过程中的中间时刻作为声速剖面的代表时间,则南海8条声速剖面各代表时间如表4所示.如图5a为南海8条实测声速剖面,图5b为各声速剖面相对于参考观测时间所用声速剖面之间的声速变化.如下图6所示,单个海底站点仿真所用到的声速剖面,其中图6a为时变声速剖面的仿真,图6b为各声速剖面相对于参考观测时间所用声速剖面之间的声速的变化.对比图5a和图6a,可以了解到仿真生成的声速剖面与实际观测的声速剖面具有近似的趋势形状,由于仿真的基准参考时间与实测SVP1的代表时间相同,且仿真过程中涵盖了SVP2,从图5b可以看出SVP2相对于SVP1的声速变化范围为-3.1799～1.5313 m·s-1,仿真声速剖面随时间变化在此范围内.

表4 实测SVP代表时刻

图5 实测声速剖面(a)及其与基准剖面(b)的声速变化

图6 仿真声速剖面(a)及其与基准剖面(b)的声速变化

3.2 仿真试验解算及其结果分析

为验证参考声速剖面误差影响,增添了使用参考声速剖面进行定位(方法1).关于方法2、3的参数个数的选取的问题,在大多数情况下,B样条基函数的节点间隔约为10～20 min(Watanabe et al.,2020).因此,本文B样条基函数系数参数的个数定为12个.

根据表5的解算方法设计对表3的不同控制网的海底虚拟基准点进行解算,控制网高度角变化范围如图7所示,海底定位的精度衰减因子(Position Dilution of Precision,PDOP)、Fisher信息量(法方程系数矩阵行列式)如图8所示,方法2、3反演的加权声速变化如图9所示,虚拟基准点解算精度结果如表6所示.

表5 解算方法设计

表6 五种网型定位虚拟中心点解算精度(单位:cm)

图8 五种控制网PDOP和Fisher信息量

图9 五种网型声速变化解算结果

五种网型反演声速变化解算结果和虚拟中心点解算结果精度分析如下:

(1)图9表明,两种解算方法都可较好拟合实际加权声速变化趋势,且网型3的拟合最优,相对于真值RMSE分别为0.438 cm·s-1和0.452 cm·s-1,网型4的曲线拟合最差,RMSE分别为1.063 cm·s-1和1.033 cm·s-1.

(2)表6表明,参考声速剖面误差可产生分米级的误差影响.两种解算方法在解算偏差、标准差和均方根误差仅存在毫米级的微小差异,认为两种解算方法在厘米级精度量级是等价的.

(3)对比表6中不同网型的定位精度可以看出,网型1与网型2的区别主要在于海面航迹尺寸不同,如表3所示,相比于网型2,网型1的PDOP更优,Fisher量更大,定位精度更好,如图8所示;网型3具有很强的垂向几何结构如图7所示,这也导致网型3的解算结果中高程方向具有较高精度的原因;网型4的图形设计是对网型1和网型2图形上的叠加,从解算结果上来看,网型4中方法1的解算结果优于网型1和网型2,解算结果偏差次于网型1、2,标准差优于网型1、2,其原因可能是嵌套形图形设计可以消除一部分系统误差的影响;网型5的图形设计是网型1、2、3总体上的叠加,其中既有丰富的高度角变化信息如图7所示,又有很好的几何构型如图8所示,具有最佳定位结果.从网型1、网型4、网型5这三种网型整体来看,随着网型设计中海面观测图形逐渐完备,三个分量的定位精度逐渐提高,即良好的海面测线设计有利于削弱系统误差的影响.

4 结论

本文提出了“信号传播时间追赶仿真法”和“带周期函数约束的声速场建模”用于更加真实模拟主动式声呐观测模式中声学信号双程传播过程和声速场周期时间变化.信号传播时间追赶仿真法解决了主动式声呐双程观测时间和接收位置难以确定的问题,而带周期函数约束的声速场建模利用了测区海洋环境观测先验信息,可获得更为逼真的海洋声速场变化.仿真试验结果表明,声速剖面校正法与传播时间校正法都可以较好的反演出实际加权平均声速的变化趋势,且两种方法解算结果相差毫米级,即在目前水下定位厘米级精度水平下是等价的;此外,对比五种控制网型解算结果发现,良好的海面测线图形设计有利于削弱参考声速剖面误差影响.

致谢感谢两名匿名评审专家提出的宝贵建议,感谢编辑提供的帮助.