On the principle and applications of the joint adjustment method for GNSS/A positioning of the underwater geodetic control point-(Invited)
-
摘要: 为分析海面换能器位置误差对海底大地控制点定位精度的影响规律，我们提出了GNSS声学(GNSS/Acoustic，GNSS/A )定位联合平差(joint adjustment, JA)方法. 针对联合平差方法理论证明不足及无深海数据验证的问题，本文进一步阐述了联合平差方法的原理，给出了联合平差方法严密的优越性证明和精度评定公式. 最后采用松花湖数据和南海数据进行验证，联合平差方法可以提高传统GNSS/A方法2%~26%的定位精度，展示联合平差方法在湖区和深海区的应用效果.Abstract: On the impact of the coordinate errors of the acoustic transducer on the positioning accuracy of the undersea geodetic control point, we proposed a joint adjustment method of the seafloor geodetic control point for GNSS/A underwater precise positioning, whose crux was treating the positions of both transducer and transponder as unknown parameters in the acoustic ranging equation. The paper further expounds the principle of the joint adjustment method, for the shortcomings in the proof of the superiority of the joint adjustment method and the lack of verification of deep-sea data, the strict superiority proof and the precision evaluation formula of the joint adjustment method are given in this paper. Using Songhua Lake data and South China Sea data for verification, the precision of seafloor geodetic control point positioning with the proposed method is improved by 2% to 26%. The application effect of the joint adjustment method in the lake area and the deep-sea area is shown.
-
0. 引 言
海底大地控制网是海洋导航定位技术的基础设施[1]. 目前，主要以GNSS声学(GNSS/Acoustic，GNSS/A)定位技术为主，根据非差/差分定位模型实现海底大地控制点的精确标定. 由于复杂海洋环境的影响和水下测距技术的限制，海底大地控制点的精确标定需要发展更加完备且高精度的GNSS/A定位数据处理方法[1]. 国际大地测量协会(International Association of Geodesy，IAG)于2019年设立了海洋大地测量工作组，目的是推动海洋空间基准与水下定位技术的发展[2-3].
为了提高海底大地控制点应答器的定位精度，研究学者先后提出了航迹规划、差分定位和附加深度约束等技术手段. 航迹规划通过增加观测值的冗余性和几何对称性，提高海底应答器平面坐标的定位精度，但高程方向的定位精度还有待提高[4-5]. 历元间差分能够削弱观测值的共模误差，比如硬件延迟误差和由海洋长周期误差引起的声速代表性误差等，提高了海底应答器的定位精度，但无法消除由海洋短周期误差引起的声速代表性误差[6-7]. 基于压力传感器提供的深度值，提出了深度差约束的水下定位方法，提高了海底应答器垂直方向的定位精度[8-11]. 但是，这些方法均是将海面换能器的位置作为没有误差的已知值，结合声学测距来计算海底应答器的位置[12]. 海面换能器位置不可避免的存在定位误差，由GNSS定位提供的海面换能器位置在各个历元具有不同的定位精度[13]. 因此，应精确评定由GNSS定位提供的海面换能器位置精度信息，采用更为合理的数据处理模型来减弱声学换能器位置误差对海底大地控制点定位的影响.
1. 传统GNSS/A定位平差方法
如图1所示，GNSS/A定位技术通常借助海面平台(如测量船)实施水下声呐测距来反演海底大地控制点的位置. 测量船上装载了GNSS接收机、姿态测量仪和声学换能器等设备. 通常采用GNSS定位技术，如实时动态(real time kinematic，RTK)技术或者精密单点定位(precise point positioning，PPP)技术来获取GNSS天线位置，然后根据精密工程测量方法和船载姿态仪数据，将GNSS天线位置归算到换能器信标中心位置[14]. 声学换能器向海底应答器发送声学信号并记录海底应答器的回波信息，可得到声学换能器至海底应答器的传播时间，根据声速和传播时间计算可得到声学换能器至海底应答器的声学距离.
$$ {\rho }_{i}={f}_{i}({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r})+{\varepsilon }_{i} $$ (1) 式中：
$ i=\mathrm{1,2},\cdots ,n $ 为观测历元；$ {\;\rho }_{i}=c{t}_{i} $ 为海底应答器到换能器的观测距离，$ c $ 为声速，$ {t}_{i} $ 为声信号单向传播时间；$ f({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}) $ 为声学换能器至海底应答器的几何距离；$ {\boldsymbol{X}}_{i}^{s}=[{x}_{i}^{s},{y}_{i}^{s},{z}_{i}^{s}{]}^{{\rm{T}}} $ 为声学换能器的坐标；$ {\boldsymbol{X}}^{r}=[{x}^{r}, {y}^{r},{z}^{r}{]}^{\rm{T}} $ 是海底应答器的坐标；$ {\varepsilon }_{i} $ 为观测误差.忽略海面换能器的位置误差，将式(1)进行线性化，线性化观测方程的向量表达式为
$$ \boldsymbol{L}=\boldsymbol{A}{\rm{d}}{\boldsymbol{X}}^{r}+\boldsymbol{\varepsilon } $$ (2) 式中：
$ \boldsymbol{A} $ 为设计矩阵；${\rm{d}}{\boldsymbol{X}}^{r}=[{{\rm{d}}x}^{r},{{\rm{d}}y}^{r},{{\rm{d}}z}^{r}{]}^{\mathrm{T}}$ 为海底应答器坐标的改正数；$ {{L}_{i}=\rho }_{i}-{f}_{i}({\boldsymbol{X}}_{i}^{s},{X}_{0}^{r}) $ ,$ {\boldsymbol{L}=({L}_{1},{L}_{2},\cdots ,{L}_{n})}^{\rm{T}} $ 为观测值；$ \boldsymbol{\varepsilon }={({\varepsilon }_{1},{\varepsilon }_{2},\cdots ,{\varepsilon }_{n})}^{{\rm{T}}} $ 为观测误差；$ {\boldsymbol{X}}_{0}^{r}=[{x}_{0}^{r},{y}_{0}^{r}, {z}_{0}^{r}{]}^{\mathrm{T}} $ 为海底应答器的线性化初值.$$ {\boldsymbol{A}}=\begin{bmatrix}{\begin{array}{*{20}{c}} {\left(\dfrac{\partial {f}_{1}}{\partial {x}_{r}}\right)}_{0}&{\left(\dfrac{\partial {f}_{1}}{\partial {y}_{r}}\right)}_{0}& {\left(\dfrac{\partial {f}_{1}}{\partial {z}_{r}}\right)}_{0}\\ {\left(\dfrac{\partial {f}_{2}}{\partial {x}_{r}}\right)}_{0}&{\left(\dfrac{\partial {f}_{2}}{\partial {y}_{r}}\right)}_{0}&{\left(\dfrac{\partial {f}_{2}}{\partial {z}_{r}}\right)}_{0}\\ \vdots &\vdots &\vdots \\ {\left(\dfrac{\partial {f}_{n}}{\partial {x}_{r}}\right)}_{0}&{\left(\dfrac{\partial {f}_{n}}{\partial {y}_{r}}\right)}_{0}&{\left(\dfrac{\partial {f}_{n}}{\partial {z}_{r}}\right)}_{0}\end{array}}\end{bmatrix} $$ (3) 根据最小二乘原理，由式(2)可得海底应答器的坐标为
$$ {\boldsymbol{X}}^{{r}}={\boldsymbol{X}}_{0}^{{r}}+{\rm{d}}{\boldsymbol{X}}^{{r}}={\boldsymbol{X}}_{0}^{{r}}+({\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{1}\boldsymbol{A}{)}^{-1}{\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{1} \boldsymbol{L}$$ (4) 式中，
$ {\boldsymbol{P}}_{1} $ 声学观测值的权阵.传统GNSS/A定位方法是将海面换能器的位置当作没有误差的已知值，结合声学测距来计算海底应答器的位置. 然而，海面换能器GNSS定位在各个历元具有不同的精度，在平差过程中，海面换能器的位置误差会传递给海底待估参数，降低海底大地控制点的定位精度. 为了分析海面位置误差对海底应答器定位精度的影响，将海面换能器和海底应答器的位置参数均作为待估参数，我们提出了GNSS/A定位联合平差(joint adjustment, JA)方法[13]，并采用模拟数据和浅水区实测数据进行了验证，但是该方法存在优越性证明方面的不足和无深海数据验证的问题. 本文进一步阐述联合平差方法的原理，给出了联合平差方法严密的优越性证明和精度评定公式，并采用松花湖浅水数据和南海深水数据进行验证.
2. 顾及海面换能器位置误差的GNSS/A联合平差方法
2.1 GNSS/A联合平差的函数模型
将海面换能器位置和海底应答器位置均作为未知参数，将式(1)线性化可得
$$ \begin{aligned} {\rho }_{i} - {f}_{i}\left({\boldsymbol{X}}_{i,0}^{s},{\boldsymbol{X}}_{0}^{r}\right) =&{\left(\frac{\partial {f}_{i}\left({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}\right)}{\partial {x}_{i}^{s}}\right)}_{0}{\rm{d}}{x}_{i}^{s}\\& +{\left(\frac{\partial {f}_{i}\left({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}\right)}{\partial {y}_{i}^{s}}\right)}_{0}{\rm{d}}{y}_{i}^{s}\\& +{\left(\frac{\partial {f}_{i}\left({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}\right)}{\partial {z}_{i}^{s}}\right)}_{0}{\rm{d}}{z}_{i}^{s}-{\left(\frac{\partial {f}_{i}\left({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}\right)}{\partial {x}^{r}}\right)}_{0}{\rm{d}}{x}^{r} \\& - {\left(\frac{\partial {f}_{i}\left({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}\right)}{\partial {y}^{r}}\right)}_{0}{\rm{d}}{y}^{r} - {\left(\frac{\partial {f}_{i}\left({\boldsymbol{X}}_{i}^{s},{\boldsymbol{X}}^{r}\right)}{\partial {z}^{r}}\right)}_{0}{\rm{d}}{z}^{r}\\&+{\varepsilon }_{i}\end{aligned} $$ (5) 式中：
${\boldsymbol{X}}_{i,0}^{s}=[{x}_{i,0}^{s},{y}_{i,0}^{s},{z}_{i,0}^{s}{]}^{{\rm{T}}}$ 为声学换能器的线性化初值；${\boldsymbol{X}}_{0}^{r}$ 含义同上.假设测量船围绕海底应答器做走航运动，共采集了
$ n $ 个历元的观测数据，则线性化观测方程的向量表达式为$$ \boldsymbol{L}=\boldsymbol{B}{\rm{d}}{\boldsymbol{X}}^{s}+{\boldsymbol{A}}{\rm{d}}{\boldsymbol{X}}^{r}+\boldsymbol{\varepsilon } $$ (6) 式中：
$ \boldsymbol{L} $ 、$ \boldsymbol{A} $ 、$ {\rm{d}}{\boldsymbol{X}}^{r} $ 、$ \boldsymbol{\varepsilon } $ 含义同上；$ {\rm{d}}{\boldsymbol{X}}^{s}=[{{\rm{d}}x}^{s},{{\rm{d}}y}^{s},{{\rm{d}}z}^{s}{]}^{\mathrm{T}} $ 为海面换能器的坐标改正数.$$ {\boldsymbol{B}}=\begin{bmatrix}{\begin{array}{*{20}{c}}{\left(\dfrac{\partial {f}_{1}}{\partial {x}_{1}^{s}}\right)}_{0}& {\left(\dfrac{\partial {f}_{1}}{\partial {y}_{1}^{s}}\right)}_{0}& {\left(\dfrac{\partial {f}_{1}}{\partial {z}_{1}^{s}}\right)}_{0}\;&\;&\;&\;&\;&\;&\;&\\ \;&\;&\;&{\left(\dfrac{\partial {f}_{2}}{\partial {x}_{2}^{s}}\right)}_{0}& {\left(\dfrac{\partial {f}_{2}}{\partial {y}_{2}^{s}}\right)}_{0}& {\left(\dfrac{\partial {f}_{2}}{\partial {z}_{2}^{s}}\right)}_{0}\\ \;&\;&\;&\;&\;&\;&\ddots \\ \;&\;&\;&\;&\;&\;&\;&{\left(\dfrac{\partial {f}_{n}}{\partial {x}_{n}^{s}}\right)}_{0}& {\left(\dfrac{\partial {f}_{n}}{\partial {y}_{n}^{s}}\right)}_{0}& {\left(\dfrac{\partial {f}_{n}}{\partial {z}_{n}^{s}}\right)}_{0}\end{array}}\end{bmatrix} $$ (7) 将海面换能器位置作为虚拟观测，则虚拟观测方程为
$$ {\boldsymbol{L}}_{{{x}}^{s}}={\boldsymbol{E}}_{3n}{\rm{d}}{\boldsymbol{X}}^{s}+{{\boldsymbol{\varepsilon}} }_{{{x}}^{s}} $$ (8) 式中：
$ {\boldsymbol{L}}_{{{x}}^{s}} $ 为虚拟观测值；$ {\boldsymbol{E}}_{3n} $ 为$ 3n\times 3n $ 的块对角矩阵，$ {\boldsymbol{E}}_{3n}={\rm{blkdiag}}(\boldsymbol{E},\boldsymbol{E},\cdots ,\boldsymbol{E}) $ ；$ \boldsymbol{E} $ 为$ 3\times 3 $ 的单位阵；$ {\boldsymbol{\varepsilon }}_{{{x}}^{{s}}} $ 为$ 3n\times 1 $ 的虚拟观测误差. 根据式(6)和式(8)，可建立GNSS/A联合平差的函数模型$$ \left\{\begin{array}{l}\boldsymbol{L}=\boldsymbol{B}{\rm{d}}{\boldsymbol{X}}^{s}+\boldsymbol{A}{\rm{d}}{\boldsymbol{X}}^{r}+{\boldsymbol{\varepsilon}} \\ {\boldsymbol{L}}_{{{x}}^{s}}={\boldsymbol{E}}_{3n}{\rm{d}}{\boldsymbol{X}}^{s}+{{\boldsymbol{\varepsilon}} }_{{{x}}^{s}}\end{array}\right. $$ (9) 2.2 GNSS/A联合平差的参数估计方法
根据最小二乘原理，由式(9)可得GNSS/A联合平差的法方程为
$$ \left\{\begin{array}{l}\left({\boldsymbol{B}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{B}+{\boldsymbol{P}}_{{{x}}_{{s}}}\right){\rm{d}}{\boldsymbol{X}}^{{s}}+{\boldsymbol{B}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{A}{\rm{d}}{\boldsymbol{X}}^{{r}}={\boldsymbol{B}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{L}+{\boldsymbol{P}}_{{{x}}_{{s}}}{\boldsymbol{L}}_{{{x}}_{{s}}}\\ {\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{B}{\rm{d}}{\boldsymbol{X}}^{{s}}+{\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{A}{\rm{d}}{\boldsymbol{X}}^{{r}}={\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{L}\end{array}\right. $$ (10) 式中：
$ {\boldsymbol{P}}_{{{1}}} $ 含义同上；$ {\boldsymbol{P}}_{{{x}}^{\boldsymbol{s}}} $ 为虚拟观测值的权阵，由换能器GNSS定位的方差-协方差信息确定.若根据式(10)直接进行求解，会面临法矩阵计算量大的问题. 假设海面调查船共采集了800个历元的观测值，法矩阵维数为2403×2403，且法矩阵的维数随着观测历元数而增加，将会占用很大的储存空间和计算资源. 相对于海面换能器的位置，我们更关注海底应答器的位置. 为此，可将法方程进行等价转换，消除海面换能器的位置参数[15]. 将式(10)第一式左乘以
${\boldsymbol{M}=-\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{B}({\boldsymbol{B}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{B}+{\boldsymbol{P}}_{{{x}}_{{s}}}{)}^{-1}$ 再与式(10)第二式相加，化简可得$$\begin{aligned} &\left({\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{A}+\boldsymbol{M}{\boldsymbol{B}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{A}\right){\rm{d}}{\boldsymbol{X}}^{{r}}=\\&{\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{L}+ \boldsymbol{M}{\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{L}+\boldsymbol{M}{\boldsymbol{P}}_{{{x}}_{{s}}}{\boldsymbol{L}}_{{{x}}_{{s}}} \end{aligned} $$ (11) 根据式(11)，可得到海底应答器的坐标改正数为
$$\begin{aligned} {\rm{d}}{\boldsymbol{X}}^{{r}}=&({\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{P}}_{1}\boldsymbol{A}+\boldsymbol{M}{\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{A}{)}^{-1}\\&({\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{L}+{\boldsymbol{M}\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{P}}_{1}\boldsymbol{L}+\boldsymbol{M}{\boldsymbol{P}}_{{{x}}_{{s}}}{\boldsymbol{L}}_{{{x}}_{{s}}}) \end{aligned} $$ (12) 由式(12)可以发现，通过对法方程进行等价转换，消除了海面换能器的位置参数，直接得到海底应答器的位置参数，解决了由于待估参数过多引起的计算量大的问题.
2.3 联合平差优越性的证明
根据协方差传播律[16]，结合式(4)，传统GNSS/A定位方法计算得到的海底应答器的协方差阵
${\boldsymbol{D}}_{\rm{L}\rm{S}}$ 为$$ {\boldsymbol{D}}_{\rm{LS}}={\sigma }_{0}^{2}({\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{{{1}}}\boldsymbol{A}{)}^{-1} $$ (13) 式中，
$ {\sigma }_{0}^{2} $ 为单位权方差. 根据式(9)，联合平差的函数模型可表达为$$ \boldsymbol{S}=\boldsymbol{G}{\rm{d}}{\boldsymbol{X}}^{r}+\boldsymbol{H}{\rm{d}}{\boldsymbol{X}}^{s}+{\boldsymbol{\varepsilon }}_{{l}}=\left[\begin{array}{cc}\boldsymbol{G}& \boldsymbol{H}\end{array}\right]\left[\begin{array}{c}{\rm{d}}{\boldsymbol{X}}^{r}\\ {\rm{d}}{\boldsymbol{X}}^{s}\end{array}\right]+{\boldsymbol{\varepsilon }}_{l} $$ (14) 式中：
$ \boldsymbol{S} $ 为观测值向量；$\boldsymbol{S}={\left[\begin{array}{cc}{\boldsymbol{L}}^{\rm{T}}& {{\boldsymbol{L}}_{{{x}}^{s}}^{\rm{T}}}\end{array}\right]^{\mathrm{T}}}$ ；$\boldsymbol{G}={\left[\begin{array}{cc}{\boldsymbol{A}}^{\rm{T}}& {\bf{0}}\end{array}\right]}^{\mathrm{T}}$ ；$\boldsymbol{H}={\left[\begin{array}{cc}{\boldsymbol{B}}^{\rm{T}}& {\boldsymbol{E}}_{3n}\end{array}\right]}^{\rm{T}}$ ；$ {\boldsymbol{\varepsilon }}_{l} $ 为误差向量；$ {\boldsymbol{\varepsilon }}_{l}={\left[\begin{array}{cc}{\boldsymbol{\varepsilon }}^{\rm{T}}& {\boldsymbol{\varepsilon }}_{{x}^{s}}^{\rm{T}}\end{array}\right]}^{\rm{T}} $ ；$ \boldsymbol{P} $ 为$ 4n\times 4n $ 的块对角矩阵；$ \boldsymbol{P}= {\rm{blkdiag}}\left({\boldsymbol{P}}_{1},{\boldsymbol{P}}_{{{x}}_{\boldsymbol{s}}}\right) $ . 根据最小二乘原理，由式(14)可得GNSS/A联合平差的法方程为$$ \left[\begin{array}{cc}{\boldsymbol{M}}_{11}& {\boldsymbol{M}}_{12}\\ {\boldsymbol{M}}_{21}& {\boldsymbol{M}}_{22}\end{array}\right]\left[\begin{array}{c}{\rm{d}}{\boldsymbol{X}}^{r}\\ {\rm{d}}{\boldsymbol{X}}^{s}\end{array}\right]=\left[\begin{array}{c}{\boldsymbol{W}}_{1}\\ {\boldsymbol{W}}_{2}\end{array}\right] $$ (15) 式中：
${\boldsymbol{M}}_{11}={\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\boldsymbol{G}$ ；${\boldsymbol{M}}_{12}={\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\boldsymbol{H} $ ；${\boldsymbol{M}}_{21}={\boldsymbol{H}}^{\rm{T}}\boldsymbol{P}\boldsymbol{G} $ ；${\boldsymbol{M}}_{22}={\boldsymbol{H}}^{\rm{T}}\boldsymbol{P}\boldsymbol{H} $ ；${\boldsymbol{W}}_{1}={\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\boldsymbol{S} $ ；${\boldsymbol{W}}_{2}={\boldsymbol{H}}^{\rm{T}}\boldsymbol{P}\boldsymbol{S} $ . 对式(15)做等价变换[15]，得：$$ \boldsymbol{C}\left[\begin{array}{cc}{\boldsymbol{M}}_{11}& {\boldsymbol{M}}_{12}\\ {\boldsymbol{M}}_{21}& {\boldsymbol{M}}_{22}\end{array}\right]\left[\begin{array}{c}{\rm{d}}{\boldsymbol{X}}^{r}\\ {\rm{d}}{\boldsymbol{X}}^{s}\end{array}\right]=\boldsymbol{C}\left[\begin{array}{c}{\boldsymbol{W}}_{1}\\ {\boldsymbol{W}}_{2}\end{array}\right] $$ (16) $$ \left[\begin{array}{cc}{\boldsymbol{M}}_{11}& {\boldsymbol{M}}_{12}\\ {\boldsymbol{M}}_{1}& {\boldsymbol{0}}\end{array}\right]\left[\begin{array}{c}{\rm{d}}{\boldsymbol{X}}^{r}\\ {\rm{d}}{\boldsymbol{X}}^{s}\end{array}\right]=\left[\begin{array}{c}{\boldsymbol{W}}_{1}\\ {\boldsymbol{R}}_{1}\end{array}\right] $$ (17) $$ \boldsymbol{C}=\left[\begin{array}{cc}{\boldsymbol{E}}_{3}& {\boldsymbol{0}}\\ {\boldsymbol{E}}_{3}& -\boldsymbol{Z}\end{array}\right] $$ (18) $$ \boldsymbol{Z}={\boldsymbol{M}}_{12}{\boldsymbol{M}}_{22}^{-1} $$ (19) $$ {\boldsymbol{M}}_{1}={\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}({\boldsymbol{E}}_{4n}-\boldsymbol{H}{\boldsymbol{M}}_{22}^{-1}{\boldsymbol{H}}^{\rm{T}}\boldsymbol{P})\boldsymbol{G} $$ (20) $$ {\boldsymbol{R}}_{1}={\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}({\boldsymbol{E}}_{4n}-\boldsymbol{H}{\boldsymbol{M}}_{22}^{-1}{\boldsymbol{H}}^{\rm{T}}\boldsymbol{P})\boldsymbol{S} $$ (21) 式中：
$ {\boldsymbol{0}} $ 是零矩阵；$ {\boldsymbol{E}}_{4n} $ 是$ 4n\times 4n $ 的单位阵；$ {\boldsymbol{E}}_{3} $ 是$ 3\times 3 $ 的单位阵. 式(17)的第二式可表达为$$ {\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\left({\boldsymbol{E}}_{4n}-\boldsymbol{H}{\boldsymbol{M}}_{22}^{-1}{\boldsymbol{H}}^{\rm{T}}\boldsymbol{P}\right)\boldsymbol{G}{\rm{d}}{\boldsymbol{X}}^{r} ={\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}({\boldsymbol{E}}_{4n}-\boldsymbol{H}{\boldsymbol{M}}_{22}^{-1}{\boldsymbol{H}}^{\rm{T}}\boldsymbol{P})\boldsymbol{S} $$ (22) 令
$$ \boldsymbol{J}=\boldsymbol{H}{\boldsymbol{M}}_{22}^{-1}{\boldsymbol{H}}^{\rm{T}}\boldsymbol{P} $$ (23) 则：
$$ {({\boldsymbol{E}}_{4n}-\boldsymbol{J})}^{2}={\boldsymbol{E}}_{4n}-\boldsymbol{J} $$ (24) $$ {({\boldsymbol{E}}_{4n}-\boldsymbol{J})}^{{\rm{T}}}\boldsymbol{P}=\boldsymbol{P}({\boldsymbol{E}}_{4n}-\boldsymbol{J}) $$ (25) $ {\boldsymbol{E}}_{4n}-\boldsymbol{J} $ 是幂等矩阵，$ {({\boldsymbol{E}}_{4n}-\boldsymbol{J})}^{\rm{T}}\boldsymbol{P} $ 是对称矩阵.根据式(23)~(25)，式(22)可改写为
$$ {\boldsymbol{G}}^{\rm{T}}{\left({\boldsymbol{E}}_{4n}-\boldsymbol{J}\right)}^{\rm{T}}\boldsymbol{P}\left({\boldsymbol{E}}_{4n}-\boldsymbol{J}\right)\boldsymbol{G}{\rm{d}}{\boldsymbol{X}}^{r} ={\boldsymbol{G}}^{\rm{T}}{\left({\boldsymbol{E}}_{4n}-\boldsymbol{J}\right)}^{\rm{T}}\boldsymbol{P}\boldsymbol{S} $$ (26) 令
$$ \boldsymbol{F}=({\boldsymbol{E}}_{4n}-\boldsymbol{J})\boldsymbol{G} $$ (27) 根据式(26)可得海底应答器的坐标改正数为
$$ {\rm{d}}{\boldsymbol{X}}^{r}={{(\boldsymbol{F}}^{\rm{T}}\boldsymbol{P}\boldsymbol{F})}^{-1}{\boldsymbol{F}}^{\rm{T}}\boldsymbol{P}\boldsymbol{S} $$ (28) 则海底应答器的坐标为
$$ {\boldsymbol{X}}^{{r}}={\boldsymbol{X}}_{{{0}}}^{{r}}+{\rm{d}}{\boldsymbol{X}}^{{r}}={\boldsymbol{X}}_{{{0}}}^{{r}}+{{(\boldsymbol{F}}^{\rm{T}}\boldsymbol{P}\boldsymbol{F})}^{-1}{\boldsymbol{F}}^{\rm{T}}\boldsymbol{P}\boldsymbol{S} $$ (29) 根据式(28)，联合式(4)可得
$$\begin{aligned} {\rm{tr}}\left({\boldsymbol{D}}_{\rm{JA}}\right)=&{\sigma }_{0}^{2}{\rm{tr}}{{(\boldsymbol{F}}^{\rm{T}}\boldsymbol{P}\boldsymbol{F})}^{-1} \\=&{\sigma }_{0}^{2}{\rm{tr}}\left({\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\boldsymbol{G}-{\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\boldsymbol{H}{\boldsymbol{M}}_{22}^{-1}{\boldsymbol{H}}^{\rm{T}}\boldsymbol{P}\boldsymbol{G}\right) \\&< {\sigma }_{0}^{2}{\rm{tr}}\left({\boldsymbol{G}}^{\rm{T}}\boldsymbol{P}\boldsymbol{G}\right)={\sigma }_{0}^{2}{\rm{tr}}({\boldsymbol{A}}^{\rm{T}}{\boldsymbol{P}}_{1}\boldsymbol{A}{)}^{-1} \\=&{\rm{tr}}\left({\boldsymbol{D}}_{\rm{LS}}\right)\end{aligned} $$ (30) 式中，
$ {\sigma }_{0}^{2} $ 为单位权方差.由式(30)可以看出，由联合平差计算得到的海底大地控制点坐标的方差-协方差矩阵的迹小于传统的GNSS/A方法. 这表明在GNSS/A定位中，联合平差方法能够提高传统GNSS/A定位方法的精度. 这种证明方法将联合平差问题转化成了常规的间接平差问题，证明简单，易于理解，不同于文献[13]采用的拉格朗日极值证明方法.
另外，由式(28)可以看出，海底应答器的坐标改正数的计算公式和间接平差的计算公式形式上完全类似，这样就把复杂的联合平差问题转化成了常规的间接平差问题，根据式(28)很容易得到严密的精度评定公式为[16]
$$ {\boldsymbol{D}}_{{\boldsymbol{X}}^{r}}={{{\sigma }_{0}^{2}(\boldsymbol{F}}^{\rm{T}}\boldsymbol{P}\boldsymbol{F})}^{-1} $$ (31) 3. 实 验
3.1 松花湖数据测试
采用松花湖实测数据验证新方法的有效性. 测量船配备了康斯伯格Simrad EM3000多波束系统、徕卡GNSS-1200接收机、Seatex MRU-05姿态测量装置、AML声速计和压力计等. 图2(a)给出了多波束测深数据，湖底平均深度约为60 m. 声速剖面仪获得的声速结构如图2(b)所示，表面声速为1465.34 m/s. 采用分层等梯度射线追踪法进行声线改正. 三个海底应答器放置在湖床中，走航船分别对三个海底应答器进行圆走航观测，三个海底应答器编号为C02、C04和C08.
分别采用传统GNSS/A定位方法(简写为LS方法)和联合平差法(简写为JA 方法)进行计算，结果如表1所示. 对于编号为C02、C04和C08的湖底应答器，传统GNSS/A定位方法东(east，E)、北(north， N)、天顶(up，U)三个方向的均方根(root mean square，RMS)分别为0.024 1 m、0.022 3 m和0.022 2 m，联合平差E、N、U三个方向的RMS分别为0.0232 m、0.014 6 m和0.014 1 m，可见传统GNSS/A方法的内符合精度低于联合平差方法. 这主要是因为该方法将海面换能器的位置作为没有误差的已知值，在平差过程中，海面换能器的位置误差会传递给湖底待估参数，降低了湖底大地控制点的定位精度. 将压力传感器提供的深度值作为真值来验证两种方法的外符合精度. 传统GNSS/A定位方法U方向的定位偏差分别为0.029 1 m，0.264 2 m和0.105 1 m，联合平差U方向的定位偏差分别为0.028 4 m、0.210 5 m和0.082 3 m，相比于传统GNSS/A定位方法，湖底大地控制点的高程定位精度分别提高了2%、20%和22%. 由此可知，新方法的外符合精度优于传统GNSS/A定位方法.
表 1 海底应答器定位结果与精度统计m PRN 方法 海底应答器定位结果 偏差_U RMS E N U E N U C02 LS 315690.3981 4841955.1130 −60.7708 0.0291 0.0164 0.0166 0.0061 JA 315690.3985 4841955.1128 −60.7718 0.0284 0.0157 0.0160 0.0058 C04 LS 315697.7689 4841835.0949 −59.8797 0.2642 0.0148 0.0150 0.0073 JA 315697.7938 4841835.0788 −59.9355 0.2105 0.0096 0.0098 0.0051 C08 LS 315736.3626 4841896.6657 −60.2694 0.1051 0.0158 0.0145 0.0057 JA 315736.4328 4841896.6636 −60.2922 0.0823 0.0096 0.0096 0.0038 图3给出了三个湖底应答器平差后的声学测距残差值. 从图3可以看出，联合平差方法的声学测距残差比传统GNSS/A定位方法的残差更集中. 计算残差的均方根误差以评估湖底大地控制点的精度. 对于编号为C02、C04和C08的湖底大地控制点，联合平差方法残差的RMS分别为0.084 6 m、0.071 1 m和0.050 6 m，传统GNSS/A定位方法的残差RMS分别为0.098 2 m、0.083 4 m和0.059 7 m，新方法可以获得更加稳定的解算结果.
3.2 南海数据测试
实验区域海底地形平坦，水深约为2 000 m. 在海底布设了1个海底大地控制点，在实验区域内以圆形航迹对海底大地控制点进行持续性观测，如图4(a)所示. 测量期间进行了声速剖面测量，如图4(b)所示. 严格测定船载声学换能器、GNSS天线在船体坐标系下的坐标，利用姿态测量设备获得测量船姿态，进行坐标转换获取声学换能器的位置. 采用分层等梯度射线追踪法进行声线改正.
分别采用传统GNSS/A定位方法(简写为LS方法)和联合平差法(简写为JA 方法)进行解算. 表2给出了传统GNSS/A定位方法和联合平差定位的计算结果及其精度统计表. 由表2可知，两种方法在E方向的定位差值为0.039 2 m，在N方向的定位差值为0.028 3 m，在U方向的定位差值为0.083 0 m. 传统方法三个方向RMS的平方根为0.181 6 m，联合平差方法三个方向RMS的平方根为0.134 1 m，相对于传统GNSS/A定位方法，新方法的RMS降低了26%.
表 2 海底应答器定位结果及其精度统计表m 方法 E N U RMS E N U LS 2438927.5608 491999.1574 −1960.1448 0.1218 0.1235 0.0538 JA 2438927.5216 491999.1857 −1960.2278 0.0899 0.0912 0.0397 图5给出了两种定位方法的残差比较图和残差直方图. 残差直方图表明，两种定位方法的声测距残差概率密度都服从正态分布. 根据残差直方图和残差比较图可知，联合平差方法中的声学测距残差值比传统GNSS/A定位方法的声学测距残差值更集中. 计算残差的RMS来进一步评估海底大地控制点的定位精度. 联合平差方法声学测距残差的RMS值为0.495 9 m，传统的GNSS/A方法声学测距残差的RMS值为0.901 6 m，与传统的GNSS/A定位方法相比，联合平差方法的声学测距残差RMS值降低了45%.
4. 结 论
综合本文的研究结果，可得到如下结论：
1)针对联合平差方法优越性证明方面的不足和无深海数据验证的问题，本文对顾及海面换能器位置误差的GNSS/A联合平差方法的原理进行了深化和分析，给出了间接平差框架下的计算和精度评定公式，在间接平差框架下证明了新方法相比于传统GNSS/A定位方法能够获取更加稳定可靠的参数估计解.
2)通过松花湖实测数据和南海实测数据进行了验证，结果表明，联合平差方法考虑了海面换能器位置误差对海底大地控制点定位精度的影响，相比于传统GNSS/A定位方法，定位精度提高了2%~26%，声学测距残差更聚集，验证了联合平差方法在湖区和深海区的应用效果.
致谢：特别感谢崂山实验室“问海计划”项目组及深海试验全体科研人员的合作与支持.
-
表 1 海底应答器定位结果与精度统计
m PRN 方法 海底应答器定位结果 偏差_U RMS E N U E N U C02 LS 315690.3981 4841955.1130 −60.7708 0.0291 0.0164 0.0166 0.0061 JA 315690.3985 4841955.1128 −60.7718 0.0284 0.0157 0.0160 0.0058 C04 LS 315697.7689 4841835.0949 −59.8797 0.2642 0.0148 0.0150 0.0073 JA 315697.7938 4841835.0788 −59.9355 0.2105 0.0096 0.0098 0.0051 C08 LS 315736.3626 4841896.6657 −60.2694 0.1051 0.0158 0.0145 0.0057 JA 315736.4328 4841896.6636 −60.2922 0.0823 0.0096 0.0096 0.0038 
下载: 导出CSV
表 2 海底应答器定位结果及其精度统计表
m 方法 E N U RMS E N U LS 2438927.5608 491999.1574 −1960.1448 0.1218 0.1235 0.0538 JA 2438927.5216 491999.1857 −1960.2278 0.0899 0.0912 0.0397
下载: 导出CSV
-
[1] 杨元喜, 徐天河, 薛树强. 我国海洋大地测量基准与海洋导航技术研究进展与展望[J]. 测绘学报, 2017, 46(1): 1-8. DOI: 10.11947/j.AGCS.2017.20160519 [2] 杨元喜, 刘焱雄, 孙大军, 等. 海底大地基准网建设及其关键技术[J]. 中国科学:地球科学, 2020, 50(7): 936-945. [3] 刘经南, 陈冠旭, 赵建虎, 等. 海洋时空基准网的进展与趋势[J]. 武汉大学学报(信息科学版), 2019, 44(1): 17-37. [4] CHEN G X, LIU Y, LIU Y X, et al. Improving GNSS-acoustic positioning by optimizing the ship’s track lines and observation combinations[J]. Journal of geodesy, 2020, 94(6): 83-96. DOI: 10.1007/s00190-020-01389-1
[5] SATO M, FUJITA M, MATSUMOTO Y, et al. Improvement of GPS/acoustic seafloor positioning precision through controlling the ship’s track line[J]. Journal of geodesy, 2013, 87(9): 105-114. DOI: 10.1007/s00190-013-0649-9
[6] 赵爽, 王振杰, 吴绍玉, 等. 基于选权迭代的走航式水声差分定位方法[J]. 石油地球物理勘探, 2017, 52(6): 1137-1145. [7] XU P L, ANDO M, TADOKORO K. Precise, three-dimensional seafloor geodetic deformation measurements using difference techniques[J]. Earth planets and space, 2005, 57(9): 795-808. DOI: 10.1186/bf03351859
[8] 王振杰, 刘慧敏, 杨慧良, 等. 基于垂直约束的深海拖曳系统USBL/DVL组合导航算法[J]. 中国惯性技术学报, 2019, 27(5): 670-676. [9] 赵建虎, 邹亚靖, 吴永亭, 等. 深度约束的海底控制网点坐标确定方法[J]. 哈尔滨工业大学学报, 2016, 48(10): 137-141. DOI: 10.11918/j.issn.0367-6234.2016.10.020 [10] 赵建虎, 陈鑫华, 吴永亭, 等. 顾及波浪影响和深度约束的水下控制网点绝对坐标的精确确定[J]. 测绘学报, 2018, 47(3): 413-421. DOI: 10.11947/j.AGCS.2018.20170246 [11] 刘慧敏, 王振杰, 赵爽. 深度约束的浅海多目标声学定位方法[J]. 石油地球物理勘探, 2019, 54(6): 1181-1187. [12] 王毅. 石油勘探中水下高精度定位算法研究[D]. 青岛: 中国石油大学(华东), 2014. [13] ZHAO S, WANG Z J, NIE Z X, et al. Investigation on total adjustment of the transducer and seafloor transponder for GNSS/Acoustic precise underwater point positioning[J]. Ocean engineering, 2021, 221(1): 108533. DOI: 10.1016/j.oceaneng.2020.108533
[14] NIE Z X, WANG B Y, WANG Z J, et al. An offshore real-time precise point positioning technique based on a single set of BeiDou short-message communication devices[J]. Journal of geodesy, 2020, 94(9): 78. DOI: 10.1007/s00190-020-01411-6
[15] 周江文. 误差理论[M]. 北京: 测绘出版社, 1979. [16] 黄维彬. 近代平差理论及其应用[M]. 北京: 解放军出版社, 1992.
计量
- 文章访问数: 330
- HTML全文浏览量: 101
- PDF下载量: 42
