集成船载重力和GNSS的测线式陆海高程传递方法

摘要

本发明涉及一种集成船载重力和GNSS的测线式陆海高程传递方法，其利用集成船载重力仪和GNSS远距离往返陆地和海岛(礁)一次，在一次船只往返中，沿船只航线进行海洋重力测量和GNSS测量，以完成高精度测线式跨海高程基准传递的数据采集，并将基础数据进行修正和降噪处理，并利用最小二乘配置方法推导出了沿测线重力数据解算高精度航线海洋垂线偏差的公式，并应用天文水准原理，实现了海上测线式高程基准统一和精确高程传递。本发明相对于现有技术，具有所需要观测数据量相对较少、所需费用相对较低；误差小、精度高，适于远距离跨海高程传递等有益效果。

说明书

技术领域

本发明涉及一种远距离跨海高程基准传递方法，尤其涉及一种集成船载重力和GNSS的 测线式陆海高程传递方法。

背景技术

高程作为基础地理空间信息，在国民经济建设和国防建设中具有举足轻重的作用。

海洋经济，港口工程、海洋运输、海洋养殖、海洋捕捞、海洋资源勘探和开发、海底工 程、海岸带工程、海洋环境保护、海洋灾害监测、海岛(礁)开发和保护等都需要丰富的高程 资料，需要统一陆图和海图的垂直基准，以构建统一的陆海空间信息基准。

为了统一陆地和岛屿的高程基准，应进行陆海高程联测。常用的陆海高程联测方法有精 密水准测量、静力水准测量、三角高程测量、海洋动力方法、GPS(Global Positioning System， 全球定位系统)水准法和重力位方法等。其中，精密水准测量虽然精度高，但是测站视距较 短，当进行跨水域测量时，前后视距之差太大，无法实现长距离精密跨海水准测量。

静力水准法是采用连通管进行高程传递，由于跨海距离较长，不但要求连通管具有极高 的质量，而且为了保持流体静力平衡，必须保证连通管中的液体无气泡，还需考虑气压差、 密度差等因素对平衡状态的影响，费用十分昂贵。

尽管世界上多个岛屿的高程基准传递采用了三角高程测量的方法，但众所周知的是，三 角高程测量方法只能测定任意两点之间的几何高差，因此，由三角高程测量直接推求的正高 或正常高将包含大地水准面不平行而引起的误差，同时还受到垂直折光差、垂线偏差等的影 响。在距离超过10千米时照准误差较大，有时甚至无法观测。所以这种方法误差大，又对 天气环境要求高，不适合长距离跨海高程传递。

动力水准法，即验潮法，也称为海洋动力方法，又称潮位观测法，其原理是在两点进行 长期的潮位观测，求出两点的平均海平面，进而由一点的高程推求另一点的正高或正常高。 这种方法由于需要多年的潮汐观测资料，周期长、代价高。一方面，区域平均海平面与大地 水准面有米级的差异；另一方面，海平面存在海面倾斜。

因此，该方法对于远距离海岛(礁)高程传递的精度有待于改善。

GPS水准法的基本原理是首先利用陆地上大量GPS水准点拟合出区域(似)大地水准面， 再把拟合的(似)大地水准面外推到待求点，进而结合待求点的大地高求出正高/正常高来。一 方面拟合的区域(似)大地水准面精度只有分米/厘米量级，另一方面外推方法也带来误差。

这种方法的前提是待求点和陆地上的GPS水准点在同一个拟合(似)大地水准面上，随着 距离增加这个前提逐渐会失去其合理性，这种误差在距离超过20千米时甚至能达到分米级。 另外就是GPS水准点的测量误差以及GPS水准点的分布影响。该方法是利用大地高差与(似) 大地水准面差距传递高程基准，所需经费较少，周期较短，方便实用。但是，这种方法误差 较大，不适合长距离跨海高程传递。

根据重力位理论，在基于重力测量的区域(似)大地水准面精化基础上，利用GNSS技术 可以实现海岛(礁)高程传递。通过确定位于不同陆地和海岛(礁)的水准原点之间的位差，从而 统一各个国家高程基准，建立全球高程框架，提出通过水准测量和重力数据解算位差，这需 要大量的外业工作和复杂的计算，还受到政治因素的严重影响，同时海岛(礁)上也无法直接 联测到陆地。

利用大地测量边值问题解算，可以进行垂直基准连接的方法。

利用线性化固定重力边值问题进行了深圳和香港的高程比较，精度达到了厘米量级，一 方面深圳和香港距离较近，另一方面有较多的高精度重力资料。

可以提出利用重力位差技术实现陆海高程传递，确定位差是一个十分耗时且复杂的工作。

基于重力位技术的高程传递，需要大量统一基准下的重力数据和GNSS测量；

确定两地之间位差的大地测量边值问题(GBVP)的解也应用于高程统一。因为，重力 数据与地形是强相关的，因此，扣掉参考重力场模型之后的剩余重力可用于表述地球重力场 的高频分量。

然而，这些方法需要区域高程(参考于区域高程基准)，计算重力异常，重力归算的基 准就有差异。基于重力位理论，使用固定边值问题的解来进行高程基准统一。利用重力位方 法进行区域(似)大地水准面精化，需要大量的重力数据和复杂的算法，代价高、费用昂贵、 周期较长。

中国专利申请CN101957193B公开了一种海岛礁高程传递的优化方法，其通过陆地、海 岸带、海岛礁重力数据的高度集成，通过高程基准系统差改正的空间重力异常数据计算陆地 与所传递海岛礁一致的高精度整体重力似大地水准面数值模型，最后利用计算区域内的已知 高程点，通过测量已知高程点的GPS大地高和待测海岛礁点的GPS大地高，完成海岛礁的 高程传递。

上述技术方案的海岛礁高程传递的优化方法，由于其构建重力似大地水准面数值模型， 需要综合陆地、海岸带和海岛礁的重力数据，因而所需要观测数据量庞大、费用高，各种类 型数据进行统一的过程中会引进新的误差；更为关键的是，这种方法要构建精确的局部重力 似大地水准面范围相对较小，也就是说，其无法实现远距离的跨海高程传递。

发明内容

本发明的目的是，提供一种实现跨海远距离、高精度、低代价远距离跨海高程基准传递 方法。

本发明为实现上述目的需要解决的技术问题是，如何利用集成船载重力仪和GNSS远距 离往返陆地和海岛(礁)一次，在一次船只往返中，沿船只航线进行海洋重力测量和GNSS测 量，以完成高精度测线式跨海高程基准传递的数据采集，并将基础数据进行运算处理，以实 现远距离跨海高程基准传递提高精度、减小系统误差的技术问题。

本发明为解决上述技术问题所采用的技术方案是，一种集成船载重力和GNSS的测线式 陆海高程传递方法，其特征在于，包括以下步骤：

第一步，已知陆地上的高程基准点A的正高，用常用的GNSS方法测量其大地高和用 重力仪测量重力异常；

第二步，用集成船载重力仪和GNSS航船沿陆地已知点A和海岛(礁)待测点B之间的航 线上往返一个航程，测量所经过的测线上的重力数据和GNSS三维坐标；

第三步，在海岛(礁)待测点B用常用的GNSS方法测量其大地高和用重力仪测量重力异 常；

第四步，利用EGM2008地球重力场模型，构建沿船只航线的重力异常自协方差矩阵和 重力—垂线偏差互协方差矩阵；

对所采集到的沿船只航线重力数据和GNSS三维数据进行分析和处理，并利用GNSS三 维数据和处理过的实测沿船只航线的重力数据，优化重力异常方差矩阵和重力—垂线偏差互 协方差矩阵；

第五步，将第四步处理完成的重力数据、重力异常方差矩阵和重力-垂线偏差互协方差矩 阵，运用最小二乘配置原理解算出垂线偏差，将结果代入天文水准公式，计算沿船只航线点 与陆地已知点A点之间的大地水准面差距之差；

在计算的过程中，需要将计算得到的垂线偏差投影到船只测线方向上，并根据大地高、 正高和大地水准面起伏的关系计算船只航线测点的正高；

第六步，根据第五步计算出的船只沿航线测点的正高，对第四步处理过的沿航线重力数 据进行空间重力异常归算；

将归算后的重力数据运用最小二乘配置重新计算沿航线测点垂线偏差，计算结果代入天 文水准公式，重新修正沿船只航线测点大地水准面差距之差，并对A、B两点之间航线的所 有大地水准面差距之差进行积分，得到陆地已知点A和海岛(礁)未知点B之间的大地水准面 差距之差；

第七步，根据第六步算出的陆地已知点A和海岛(礁)未知点B的大地水准面差距之差， 结合第一步测得的A点和B点的大地高，即可求出未知点B点的正高。

作为优选方式，上述的集成船载重力和GNSS的测线式陆海高程传递方法，其利用集成 船载重力仪和GNSS技术实现测线式跨海高程传递的数据采集；

利用EGM2008地球重力场模型，构建重力异常自协方差矩阵和重力-垂线偏差互协方差 矩阵；

结合航线测点重力数据，利用最小二乘配置方法，求取高精度垂线偏差；

最后，利用天文水准原理实现陆海高程基准的远距离传递。

进一步优选，上述重力异常协方差矩阵中还包括有噪声方差矩阵；

所述噪声方差矩阵通过实验探求获取、通过模拟数据或经验数据得到。

上述技术方案，通过已知陆地上的高程基准点A的正高，用常用的GNSS方法测量其大 地高和用重力仪测量重力异常；用集成船载重力仪和GNSS航船沿陆地已知点A和海岛(礁) 待测点B之间的航线上往返一个航程，测量所经过的测线上的重力数据和GNSS三维坐标； 在海岛(礁)待测点B用常用的GNSS方法测量其大地高和重力异常；利用EGM2008地球重力 场模型，构建沿船只航线的重力异常方差矩阵和重力—垂线偏差互协方差矩阵；对所采集到 的沿船只航线重力数据和GNSS三维数据进行分析和处理，并利用GNSS三维数据和处理过 的实测沿船只航线的重力数据，优化重力异常方差矩阵和重力—垂线偏差互协方差矩阵等系 列技术手段的采用，其直接带来的技术效果是：

1、用集成船载重力仪和GNSS航船沿航线往返一次共进行了两次测量，可以通过共线平 差有效提高数据精度，其成本低且数据采集覆盖面宽，因而适合远距离海岛(礁)高程的精确 确定基础数据的采集，且基础数据的采集成本低；

2、利用EGM2008地球重力场模型重力位系数标准差构建重力异常方差矩阵和重力-垂线 偏差互协方差矩阵；对船载实测重力数据进行去噪和平差处理，可以提高其精度；

3、利用处理后的重力数据修正重力异常方差矩阵和重力-垂线偏差协方差矩阵；最后， 利用最小二乘配置原理，可以解算出高精度垂线偏差，保证了海上测线式高程基准统一和精 确高程传递，可以满足远海岛(礁)高程的精确确定的需要；

需要说明的是：

上述技术方案中，由于在海洋动态环境下开展海洋测量，受到船只姿态、加速度、厄缶 效应、仪器噪声等的影响，因此必须对采集到的重力数据和大地高数据进行噪声滤波处理。

上述技术方案中，准确地构建重力异常方差矩阵和重力—垂线偏差互协方差矩阵是解算 高精度垂线偏差的基础。只有计算出高精度的垂线偏差，才能将天文水准原理应用于跨海高 程传递之中，因此准确的构建重力异常方差矩阵和重力—垂线偏差互协方差矩阵是很关键的 一步。

上述技术方案中，利用EGM2008地球重力场模型重力位系数标准差构建重力异常方差 矩阵和重力-垂线偏差互相方差矩阵；

对船载实测重力数据进行去噪和平差处理以提高其精度；

利用处理后的重力数据修正重力异常方差矩阵和重力-垂线偏差协方差矩阵；

最后，利用最小二乘配置原理解算高精度垂线偏差。

其中：

1、构建重力-垂线偏差互相方差矩阵过程如下：

(a)重力异常和垂线偏差的球谐展开公式的推导

根据布隆斯公式可得：$\mathrm{\Delta g}=-\frac{\partial T}{\partial r}-\frac{2}{r}T---\left(1\right)$

扰动位的球谐函数展开式：

$T(r,\theta ,\lambda )=\frac{\mathrm{GM}}{r}\underset{n-1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)---\left(2\right)$

是地球重力场减掉正常椭球重力位系 数。

将扰动位对地球半径求偏导得：

$\left(\begin{array}{c}-\frac{\partial T}{\partial r}=-[-\frac{\mathrm{GM}}{r^2}\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\\ -\frac{\mathrm{GM}}{r^2}\underset{n=1}{\overset{\infty }{\Sigma }}n{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)]\\ =\frac{\mathrm{GM}}{r^2}\underset{n=1}{\overset{\infty }{\Sigma }}(n+1){\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\end{array}\right)---\left(3\right)$

将公式(3)代入公式(1)得：

$\mathrm{\Delta g}=\frac{\mathrm{GM}}{r^2}\underset{n-1}{\overset{\infty }{\Sigma }}(n-1){\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)---\left(4\right)$

根据扰动位与垂线偏差之间的关系得：

$\left(\begin{array}{c}\xi =\frac{1}{\mathrm{\gamma r}}\frac{\partial T}{\partial \theta }\\ \eta =-\frac{1}{\gamma r\sin \theta }\frac{\partial T}{\partial \lambda }\end{array}\right)---\left(5\right)$

对扰动位求偏导得：

$\left(\begin{array}{c}\xi =\frac{1}{\mathrm{\gamma r}}\frac{\partial T}{\partial \theta }=\frac{1}{\mathrm{\gamma r}}·\frac{\mathrm{GM}}{r}\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda })\frac{d\left({\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\right)}{\mathrm{d\theta }}\\ =\frac{\mathrm{GM}}{{\mathrm{\gamma r}}^2}\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda })\frac{d\left({\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\right)}{\mathrm{d\theta }}\end{array}\right)---\left(6\right)$

$\left(\begin{array}{c}\eta =-\frac{1}{\gamma r\sin \theta }\frac{\partial T}{\partial \mathrm{\lambda \theta }}=-\frac{1}{\gamma r\sin \theta }·\frac{\mathrm{GM}}{r}\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}m(-{\bar{C}}_{\mathrm{nm}}^{*}\sin \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\cos \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\\ =\frac{\mathrm{GM}}{{\mathrm{\gamma r}}^2\sin \theta }\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}m({\bar{C}}_{\mathrm{nm}}^{*}\sin \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\cos \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\end{array}\right)---\left(7\right)$

(b)将勒让德函数展开代入重力异常和垂线偏差的球谐展开式

勒让德函数展开式为：

${\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)=\frac{1}{2^n}{(1-{\cos }^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}---\left(8\right)$

对公式(8)求偏导得：

$\left(\begin{array}{c}\frac{d\left({\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\right)}{\mathrm{d\theta }}=\frac{d\left({\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\right)}{d\left(\cos \theta \right)}\frac{d\left(\cos \theta \right)}{\mathrm{d\theta }}\\ =(\frac{1}{2^n}\frac{m}{2}\frac{{(1-{\cos }^2\theta )}^{\frac{m}{2}}}{1-{\cos }^2\theta }(-2\cos \theta )\underset{k=0}{\overset{n}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}\\ +\frac{1}{2^n}{(1-{\cos }^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{n}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}\frac{(n-m-2k)}{\cos \theta }{\left(\cos \theta \right)}^{n-m-2k})(-\sin \theta )\\ =\frac{1}{2^n}{(1-{\cos }^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}(m\cot \theta -(n-m-2k)\tan \theta )\end{array}\right)---\left(9\right)$

将公式(8)代入公式(4)和公式(7)，将公式(9)代入公式(6)得到重力异常、垂 线偏差关于经纬度和地球重力位系数和的展开式形式：

$\left(\begin{array}{c}\mathrm{\Delta g}\frac{\mathrm{GM}}{r^2}\underset{n=1}{\overset{\infty }{\Sigma }}(n-1){\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\\ =\frac{\mathrm{GM}}{r^2}\underset{n=2}{\overset{\infty }{\Sigma }}(n-1){\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_m^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\\ =\frac{\mathrm{GM}}{r^2}\underset{m=2}{\overset{\infty }{\Sigma }}(n-1){\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda })\frac{1}{2^n}{(1-c{\mathrm{os}}^2\theta )}^{\frac{m}{2}}\\ \underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}\end{array}\right)---\left(10\right)$

$\left(\begin{array}{c}\xi =\frac{\mathrm{GM}}{\gamma r^2}\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}({\bar{C}}_{\mathrm{nm}}^{*}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{nm}}\sin \mathrm{m\lambda })\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\\ \underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}(m\cot \theta -(n-m-2k)\tan \theta )\end{array}\right)---\left(11\right)$

$\left(\begin{array}{c}\eta =\frac{\mathrm{GM}}{{\mathrm{\gamma r}}^2\sin \theta }\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}m({\bar{C}}_{\mathrm{nm}}^{*}\sin \mathrm{m\lambda }-{\bar{S}}_{\mathrm{nm}}\cos \mathrm{m\lambda }){\bar{P}}_{\mathrm{nm}}\left(\cos \theta \right)\\ =\frac{\mathrm{GM}}{\gamma r^2\sin \theta }\underset{n=1}{\overset{\infty }{\Sigma }}{\left(\frac{a}{r}\right)}^n\underset{m=0}{\overset{n}{\Sigma }}m({\bar{C}}_{\mathrm{nm}}^{*}\sin \mathrm{m\lambda }-{\bar{S}}_{\mathrm{nm}}\cos \mathrm{m\lambda })\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\\ \underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}\end{array}\right)---\left(12\right)$

(c)利用EGM2008地球重力场模型，构建重力异常-垂线偏差协方差：

将公式(10)～(12)转化成如下格式：

$\left(\begin{array}{c}\mathrm{\Delta g}=\frac{\mathrm{GM}}{r^2}\left\{\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}\right[(n-1){\left(\frac{a}{r}\right)}^n\cos \mathrm{m\lambda }·\frac{1}{2^n}{(1-{\cos }^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}·{\bar{C}}_m^{*}]+\\ \underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}[(n-1){\left(\frac{a}{r}\right)}^n\sin \mathrm{m\lambda }·\frac{1}{2^n}{(1-c{\mathrm{os}}^2\theta )}^{\frac{m}{2}}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}·{\bar{S}}_{\mathrm{nm}}]\}\end{array}\right)---\left(13\right)$

$\left(\begin{array}{c}\xi =\frac{\mathrm{GM}}{\gamma r^2}\left\{\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}\right[{\left(\frac{a}{r}\right)}^n\cos \mathrm{m\lambda }·\frac{1}{2^n}{(1-{\cos }^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}(m\cot \theta -(n-m-2k)\tan \theta ·{\bar{C}}_m^{*}]+\\ \underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}[{\left(\frac{a}{r}\right)}^n\sin \mathrm{m\lambda }·\frac{1}{2^n}{(1-c{\mathrm{os}}^2\theta )}^{\frac{m}{2}}{\underset{k=0}{\overset{q}{\Sigma }}(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}(m\cot \theta -(n-m-2k)\tan {\theta ·\bar{S}}_{\mathrm{nm}}]\}\end{array}\right)---\left(14\right)$

$\left(\begin{array}{c}\eta =\frac{\mathrm{GM}}{\gamma r^2\sin \theta }\left\{\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}\right[{\left(\frac{a}{r}\right)}^{n}m·\sin \mathrm{m\lambda }·\frac{1}{2^n}{(1-{\cos }^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}·{\bar{C}}_m^{*}]-\\ \underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}[{\left(\frac{a}{r}\right)}^{n}m·\cos \mathrm{m\lambda }·\frac{1}{2^n}{(1-c{\mathrm{os}}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}{·\bar{S}}_{\mathrm{nm}}]\}\end{array}\right)---\left(15\right)$

令$(n-1){\left(\frac{a}{r}\right)}^n\cos \mathrm{m\lambda }·\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}=A_{\mathrm{nm}}$

$(n-1){\left(\frac{a}{r}\right)}^n\sin \mathrm{m\lambda }·\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}=B_{\mathrm{nm}}$

${\left(\frac{a}{r}\right)}^n\cos \mathrm{m\lambda }·\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}(m\cot \theta -(n-m-2k)\tan \theta =D_{\mathrm{nm}}$

${\left(\frac{a}{r}\right)}^n\sin \mathrm{m\lambda }·\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}(m\cot \theta -(n-m-2k)\tan \theta =E_{\mathrm{nm}}$

${\left(\frac{a}{r}\right)}^{n}m·\sin \mathrm{m\lambda }·\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}=M_{\mathrm{nm}}$

${\left(\frac{a}{r}\right)}^{n}m·\cos \mathrm{m\lambda }·\frac{1}{2^n}{(1-\mathrm{co}{s}^2\theta )}^{\frac{m}{2}}\underset{k=0}{\overset{q}{\Sigma }}{(-1)}^k\frac{1}{k!(n-k)!}\frac{(2n-2k)!}{(n-m-2k)!}{\left(\cos \theta \right)}^{n-m-2k}=N_{\mathrm{nm}}$得到：

$\mathrm{\Delta g}=\frac{\mathrm{GM}}{r^2}(\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{A}_{\mathrm{nm}}·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{B}_{\mathrm{nm}}·{\bar{S}}_{\mathrm{nm}})$

$\xi =\frac{\mathrm{GM}}{\gamma r^2}(\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{D}_{\mathrm{nm}}·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{E}_{\mathrm{nm}}·{\bar{S}}_{\mathrm{nm}})---\left(17\right)$

$\eta =\frac{\mathrm{GM}}{\gamma r^2\sin \theta }(\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{M}_{\mathrm{nm}}·{\bar{C}}_{\mathrm{nm}}^{*}-\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{N}_{\mathrm{nm}}·{\bar{S}}_{\mathrm{nm}})---\left(18\right)$

设P，Q是地面上两个点，其重力异常Δg^P、Δg^Q和垂线偏差Δξ^P、Δξ^Q、Δη^P、Δη^Q分 别是：

$\Delta g^P=\frac{\mathrm{GM}}{r_P^2}(\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{A}_{\mathrm{nm}}^P·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{B}_{\mathrm{nm}}^P·{\bar{S}}_{\mathrm{nm}})---\left(19\right)$

$\Delta g^Q=\frac{\mathrm{GM}}{r_Q^2}(\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{A}_{\mathrm{nm}}^Q·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{B}_{\mathrm{nm}}^Q·{\bar{S}}_{\mathrm{nm}})---\left(20\right)$

${\xi }^P=\frac{\mathrm{GM}}{{{\gamma }_{P}r}_P^2}(\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{D}_{\mathrm{nm}}^P·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{E}_{\mathrm{nm}}^P·{\bar{S}}_{\mathrm{nm}})---\left(21\right)$

${\xi }^Q=\frac{\mathrm{GM}}{{{\gamma }_{P}r}_P^2}(\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{D}_{\mathrm{nm}}^Q·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{E}_{\mathrm{nm}}^Q·{\bar{S}}_{\mathrm{nm}})---\left(22\right)$

${\eta }^P=\frac{\mathrm{GM}}{{{\gamma }_{P}r}_P^2\sin {\theta }_P}(\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{M}_{\mathrm{nm}}^P·{\bar{C}}_{\mathrm{nm}}^{*}+\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{N}_{\mathrm{nm}}^P·{\bar{S}}_{\mathrm{nm}})---\left(23\right)$

${\eta }^Q=\frac{\mathrm{GM}}{{{\gamma }_{Q}r}_Q^2\sin {\theta }_Q}(\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{M}_{\mathrm{nm}}^Q·{\bar{C}}_{\mathrm{nm}}^{*}-\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{N}_{\mathrm{nm}}^Q·{\bar{S}}_{\mathrm{nm}})---\left(24\right)$

根据EGM2008地球重力场模型获取地球重力位系数误差标准偏差为和根据误 差传播定律可以得到P、Q两个点之间的协方差矩阵：

$\mathrm{COV}({\mathrm{\Delta g}}^P,{\mathrm{\Delta g}}^Q)=\frac{\mathrm{GM}}{r_P^2}\frac{\mathrm{GM}}{r_Q^2}(\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{A}_{\mathrm{nm}}^P·{\sigma }_{{\bar{C}}_{\mathrm{nm}}}^2·A_{\mathrm{nm}}^Q+\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{B}_{\mathrm{nm}}^P·{\sigma }_{{\bar{S}}_{\mathrm{nm}}}^2{B}_{\mathrm{nm}}^Q)---\left(25\right)$

$\mathrm{COV}({\xi }^P,{\xi }^Q)=\frac{\mathrm{GM}}{{\gamma }_P{r}_P^2}\frac{\mathrm{GM}}{{{\gamma }_{Q}r}_Q^2}(\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{D}_{\mathrm{nm}}^P·{\sigma }_{{\bar{C}}_{\mathrm{nm}}}^2·D_{\mathrm{nm}}^Q+\underset{n=0}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{E}_{\mathrm{nm}}^P·{\sigma }_{{\bar{S}}_{\mathrm{nm}}}^2{E}_{\mathrm{nm}}^Q)---\left(26\right)$

$\mathrm{COV}({\eta }^P,{\eta }^Q)=\frac{\mathrm{GM}}{{\gamma }_P{r}_P^2\sin {\theta }_P}\frac{\mathrm{GM}}{{\gamma }_Q{r}_Q^2\sin {\theta }_Q}(\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{M}_{\mathrm{nm}}^P·{\sigma }_{{\bar{C}}_{\mathrm{nm}}}^2·M_{\mathrm{nm}}^Q+\underset{n=1}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{N}_{\mathrm{nm}}^P·{\sigma }_{{\bar{S}}_{\mathrm{nm}}}·N_{\mathrm{nm}}^Q)---\left(27\right)$

$\mathrm{COV}({\mathrm{\Delta g}}^P,{\xi }^Q)=\frac{\mathrm{GM}}{r_P^2}\frac{\mathrm{GM}}{{\gamma }_Q{r}_Q^2}(\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{A}_{\mathrm{nm}}^P·{\sigma }_{{\bar{C}}_{\mathrm{nm}}}^2·D_{\mathrm{nm}}^Q+\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{B}_{\mathrm{nm}}^P·{\sigma }_{{\bar{S}}_{\mathrm{nm}}}^2{E}_{\mathrm{nm}}^Q)---\left(28\right)$

$\mathrm{COV}({\mathrm{\Delta g}}^P,{\eta }^Q)=\frac{\mathrm{GM}}{r_P^2}\frac{\mathrm{GM}}{{\gamma }_Q{r}_Q^2\sin {\theta }_Q}(\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{A}_{\mathrm{nm}}^P·{\sigma }_{{\bar{C}}_{\mathrm{nm}}}^2·M_{\mathrm{nm}}^Q+\underset{n=2}{\overset{\infty }{\Sigma }}\underset{m=0}{\overset{n}{\Sigma }}{B}_{\mathrm{nm}}^P·{\sigma }_{{\bar{S}}_{\mathrm{nm}}}^2{N}_{\mathrm{nm}}^Q)---\left(29\right)$

2、根据最小二乘配置方法求沿测线测点垂线偏差

采集并处理后的沿测线船载重力数据为Δg₁,Δg₂,…Δg_i…，所求得的垂线偏差子午分量和 卯酉分量分别为ξ₁,ξ₂,…ξ_i…和η₁,η,…η_i…，最小二乘配置原理为：

$S=C_{\mathrm{sl}}(C_{\mathrm{ll}}^{-1}+D)L---\left(30\right)$

式中：

$S=\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\\ ...\\ {\xi }_i\\ ...\end{array}\right)$或$\left(\begin{array}{c}{\eta }_1\\ {\eta }_2\\ ...\\ {\eta }_i\\ ...\end{array}\right)$$C_{\mathrm{sl}}=\left(\begin{array}{ccccc}\mathrm{COV}({\mathrm{\Delta g}}^1,{\xi }^1)& \mathrm{COV}({\mathrm{\Delta g}}^1,{\xi }^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^1,{\xi }^i)& ...\\ \mathrm{COV}({\mathrm{\Delta g}}^2,{\xi }^1)& \mathrm{COV}({\mathrm{\Delta g}}^2,{\xi }^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^2,{\xi }^i)& ...\\ ...& ...& ...& ......& \\ \mathrm{COV}({\mathrm{\Delta g}}^i,{\xi }^1)& \mathrm{COV}({\mathrm{\Delta g}}^i,{\xi }^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^i,{\xi }^i)& ...\\ ...& ...& ...& ......& \end{array}\right)$或

$C_{\mathrm{sl}}=\left(\begin{array}{ccccc}\mathrm{COV}({\mathrm{\Delta g}}^1,{\eta }^1)& \mathrm{COV}({\mathrm{\Delta g}}^1,{\eta }^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^1,{\eta }^i)& ...\\ \mathrm{COV}({\mathrm{\Delta g}}^2,{\eta }^1)& \mathrm{COV}({\mathrm{\Delta g}}^2,{\eta }^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^2,{\eta }^i)& ...\\ ...& ...& ...& ......& \\ \mathrm{COV}({\mathrm{\Delta g}}^i,{\eta }^1)& \mathrm{COV}({\mathrm{\Delta g}}^i,{\eta }^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^i,{\eta }^i)& ...\\ ...& ...& ...& ......& \end{array}\right)$$L=\left(\begin{array}{c}\Delta g_1\\ {\mathrm{\Delta g}}_2\\ ...\\ {\mathrm{\Delta g}}_i\\ ...\end{array}\right)$

$C_{\mathrm{ll}}=\left(\begin{array}{ccccc}\mathrm{COV}({\mathrm{\Delta g}}^1,{\mathrm{\Delta g}}^1)& \mathrm{COV}({\mathrm{\Delta g}}^1,{\mathrm{\Delta g}}^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^1,{\mathrm{\Delta g}}^i)& ...\\ \mathrm{COV}({\mathrm{\Delta g}}^2,{\mathrm{\Delta g}}^1)& \mathrm{COV}({\mathrm{\Delta g}}^2,{\mathrm{\Delta g}}^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^2,{\mathrm{\Delta g}}^i)& ...\\ ...& ...& ...& ......& \\ \mathrm{COV}({\mathrm{\Delta g}}^i,{\mathrm{\Delta g}}^1)& \mathrm{COV}({\mathrm{\Delta g}}^i,{\mathrm{\Delta g}}^2)& ...& \mathrm{COV}({\mathrm{\Delta g}}^i,{\mathrm{\Delta g}}^i)& ...\\ ...& ...& ...& ......& \end{array}\right),$D为噪声方差矩阵。

垂线偏差和子午分量卯酉分量之间的关系式为：

μ_i＝ξcosA_i+ηsinA_i     (31)

A_i为船只沿测线i,i+1两测点的大地方位角。

综上所述，本发明相对于现有技术，具有所需要观测数据量相对较少、所需费用相对较 低；误差小、精度高，适于远距离跨海高程传递等有益效果。

附图说明

图1测线式跨海高程传递数据总体技术路线示意图；

图2天文水准高程传递原理图。

具体实施方式

下面结合附图，对本发明进行详细说明。

本发明的集成船载重力和GNSS的测线式陆海高程传递方法，利用集成船载重力仪和 GNSS技术实现测线式跨海高程传递的数据采集，其包括以下步骤：

第一步，已知陆地上的高程基准点A的正高h_A，用常用的GNSS方法测量其大地高H_A和 用重力仪测量重力异常Δg_A；

第二步，用集成船载重力仪和GNSS航船沿陆地已知点A和海岛(礁)待测点B之间的航 线上往返一个航程，测量所经过的测线上的重力数据Δg_i(i＝1,2,…)和GNSS三维坐标

第三步，在海岛(礁)待测点B用常用的GNSS方法测量其大地高H_B和用重力仪测量重力 异常Δg_B；

如图1所示，本发明的集成船载重力和GNSS的测线式陆海高程传递方法，其中的数据 处理共分为三大块：

1、船测重力数据处理流程：在海洋动态环境下，将重力仪安置于船只上，开展海洋重力 测量，受到船只姿态、加速度、厄缶效应、仪器噪声等的影响，对船载重力数据进行处理以 改善船载重力精度。船载重力测量数据处理流程其误差改正包括：迟后效应改正、零点漂移 改正、垂直扰动加速度剩余影响改正、水平扰动加速度剩余影响改正、厄缶效应改正等。计 算出海面高后进行航线测点空间重力异常改正。改正涉及的方法有沿航线的高斯滤波算法和 共线平差方法等。

2、航线海面高计算流程：改正包括姿态改正、迟后改正以及由吃水深度模型和浮子GNSS 确定的高度改正。历元海面高减去DTU10平均海面高得到海面高残差，采用高斯方法对残 差海面高进行滤波。采用卫星测高+验潮+潮汐模型的综合方法进行航线海面高精度验证。

船载GNSS处理采用单历元精密单点定位方法。

首先，构建船载GNSS数据质量评价和控制体系。

其次，对船载GNSS数据进行预处理，进行周跳探测与修复以及时间同步。进行初始化， 确定模糊度。以NCEP、ECMWF或UCAR模型计算对流层延迟作为初值，利用GNSS精密 星历，未知参数包括位置、接收机钟差和对流层延迟残差，进行最小二乘批处理。在优化配 置船载GNSS接收机和天线基础上，解算精确历元位置、速度、加速度和船姿。

当无法获得GNSS精密星历时，采用历史精密星历和广播星历，构建星历残差序列， 采用谱分析方法，建立残差星历预报模型，实时广播星历加上残差星历预报模型代替精密星 历。

3、沿航线的垂线偏差：以EGM2008为参考重力场，采用移去-恢复技术，利用沿航线 的重力残差数据由最小二乘配置方法解算垂线偏差。

通过船测重力数据处理，就可以确定沿航线重力数据方差矩阵和噪声矩阵。

根据垂线偏差和重力异常与位系数之间的函数关系，由EGM2008模型建立垂线偏差与重力 之间的互协方差矩阵，通过实测数据处理和分析，不断改善该互协方差矩阵。

上述数据处理方法的具体流程，包括以下步骤：

第一步，利用EGM2008地球重力场模型，构建沿船只航线的重力异常协方差矩阵 cov(Δg_i,Δg_j)和重力-垂线偏差互协方差矩阵cov(Δg_i,ξ_j)和cov(Δg_i,η_j)；

第二步，对所采集到的沿船只航线重力数据和GNSS三维大地坐标数据进行分析和处理； 对实测船载重力数据进行去噪和平差处理以提高其精度；船载GNSS测量受到船姿、船只吃 水深度、风浪场等影响，构建相应改正模型，采用浮子GNSS技术，进行单历元精密单点定 位/相对定位，通过共线平差和滤波方法改善海面高精度；优化配置船载GNSS，构建船载GNSS 数据质量评价和控制体系，实现了海洋环境下的单历元相对定位/精密单点定位，达到位置、 速度、加速度和姿态的高精度测量；利用精确测量获得的测线测点海面大地高及吃水深度对 船测重力数据进行空间改正，检验和评估海面重力数据及精度；

第三步，利用GNSS三维数据和处理过的实测沿船只航线的重力数据，优化重力异常协 方差矩阵cov(Δg_i,Δg_j)和重力-垂线偏差互协方差矩阵cov(Δg_i,ξ_j)和cov(Δg_i,η_j)；

第四步，采用移去-恢复技术，以EGM2008重力场模型为参考，根据处理后的实测GNSS 三维大地坐标数据解算模型重力异常Δg_i^model和模型垂线偏差(ξ^model_i,η^model_i)；根据处理后的实 测重力异常数据Δg_i和模型重力异常Δg_i^model计算残余重力异常Δg_i^res，

Δg_i^res＝Δg_i-Δg_i^model     (32)

第五步，将第四步处理完成的残余重力数据Δg_i^res、重力异常协方差矩阵cov(Δg_i,Δg_j)和重 力-垂线偏差互协方差矩阵cov(Δg_i,ξ_j)和cov(Δg_i,η_j)，代入最小二乘配置公式(30)解算出残 余垂线偏差(ξ^res_i,η^res_i)；模型垂线偏差(ξ^model_i,η^model_i)和残余垂线偏差(ξ^res_i,η^res_i)恢复沿测线测 点的垂线偏差(ξ_i,η_i)；

${\xi }_i={\xi }_i^{\mathrm{res}}+{{\xi }_i}^{\mathrm{mode}l}$

${\eta }_i={\eta }_i^{\mathrm{res}}+{{\eta }_i}^{\mathrm{mode}l}---\left(33\right)$

将垂线偏差(ξ_i,η_i)计算结果代入公式(31)计算垂线偏差μ_i；

第六步，根据天文水准的基本原理，有无限接近的两个点a，b，如图2所示，在a点上 n₁n₁'表示参考椭球体的法线，a₁a₁'表示大地水准面的法线，即重力方向，在ab方向上的天文 大地垂线偏差分量为μ_ab。

dN_ab表示ab两点相对于参考椭球体的大地水准面差距之差，用ds表示ab两点之间的 距离。

dN_ab＝-μ_abds     (34)

则陆地已知点A和海岛(礁)未知点B之间的大地水准面差距之差ΔN_AB可表示为：

$\Delta N_{\mathrm{AB}}=-{\int }_A^B{\mu }_i\mathrm{ds}---\left(35\right)$

上式离散表示为，

$\Delta N_{\mathrm{AB}}=-\underset{i=1}{\overset{n}{\Sigma }}{\mu }_i\Delta s_i---\left(36\right)$

将第五步垂线偏差μ_i计算结果代入天文水准公式公式(35)，计算沿船只航线i(i＝1,2,…) 点与陆地已知点A点之间的大地水准面差距之差

大地高、正高和大地水准面起伏的关系为：

H_i＝h_i+N_i     (37)

在计算的过程中，需要将计算得到的垂线偏差投影到船只测线方向上，并根据大地高、 正高和大地水准面起伏的关系计算船只航线i点的正高；

第七步，根据第六步计算出的船只沿航线测点i的正高，对第四步处理过的沿航线重力数 据Δg_i^res进行空间重力异常归算；

将归算后的残余重力数据运用最小二乘配置重新计算沿航线测点残余垂线偏差，并重新 恢复垂线偏差，计算结果代入天文水准公式，重新修正沿船只航线测点大地水准面差距之差， 并对A、B两点之间航线的所有大地水准面差距之差进行积分得到陆地已知点 A和海岛(礁)未知点B之间的大地水准面差距之差ΔN_BA；

第八步，根据第七步算出的陆地已知点A和海岛(礁)未知点B的大地水准面差距之差， 结合第一步测得的A点和B点的大地高，即可求出未知点B点的正高，

A、B两点之间的正高差为，

$\left(\begin{array}{c}{h}_{\mathrm{AB}}=-H_A+N_A+(H_B-N_B)=\Delta H_{\mathrm{AB}}-\Delta N_{\mathrm{AB}}\\ =\Delta H_{\mathrm{AB}}+\underset{i=1}{\overset{n}{\Sigma }}{\mu }_i\Delta s_i\end{array}\right)---\left(38\right)$

未知点B点的正高为：

$h_B=h_A+h_{\mathrm{AB}}=h_A+H_B-H_A+\underset{i=1}{\overset{n}{\Sigma }}{\mu }_i\Delta s_i---\left(39\right)$

上述的跨海高程基准传递方法，优选为，利用集成船载重力仪和GNSS技术实现测线式 跨海高程传递的数据采集；

高精度的垂线偏差数据是高程传递方法的基础资料，在船载重力数据处理基础上，通过 EGM2008地球重力场位系数标准差建立重力数据方差矩阵和重力-垂线偏差互协方差矩阵， 结合实测重力数据优化协方差矩阵，通过最小二乘配置方法解算高精度沿航线的垂线偏差；

如图2所示，最后，利用天文水准原理实现陆海高程基准的远距离传递；

上述的集成船载重力和GNSS的测线式陆海高程传递方法，还可以优选为，在所述重力 异常协方差矩阵中还包括有噪声方差矩阵；为了更进一步提高数据处理的精度，可以在重力 异常协方差矩阵中加入噪声方差矩阵，进行修正；

所述噪声方差矩阵的获取可以通过实验探求，也可通过模拟数据或经验数据得到； 需要说明的是，由于噪声的影响具有随机性和偶然性，其对数据影响的大小，需要通过实验 进行探索和验证。