一种基于水平张量重力梯度的垂线偏差实时计算方法与流程

本发明涉及重力测量技术领域，尤其是一种基于水平张量重力梯度的垂线偏差实时计算方法。

背景技术：

在几何大地测量中，以大地椭球面为水平基准将观测数据投影到该面上进行计算。椭球面上的铅垂线方向被定义为重力方向，如果地球的质量分布是均匀的，则椭球面上各点的法线与铅垂线方向一致。然而，地球的质量分布是不均匀的，大地椭球面上各点铅垂线方向同法线方向存在偏差，垂线偏差就是大地水准面垂线和参考椭球体法线方向间的差异，其大小和方向不随点位的不同规则变化。在对大地水平基准有高度依赖性的特殊应用领域中，对垂线偏差的影响十分敏感。例如：在弹道导弹发射时，采用经纬仪等光学设备将已知的基准方位传递给弹上惯性制导设备，此时需将经纬仪等设备相对大地参考椭球面法线精确调平，调平的依据是设备垂线轴与重力方向重合。但由于垂线偏差的存在，该方法并不能保证所需的精确调平。在高精度惯性导航系统中，重力偏差是重要的误差源，使重力加速度按该角度投影进入到水平加速度计中，积分形成水平速度和位置误差。

垂线偏差在子午面上的分量称为垂线偏差子午分量(南北垂线偏差分量)，记作ξ；在卯酉面上的分量称为垂线偏差卯酉分量，记作η。垂线偏差的大小一般在3″～5″之间，最大可达20″～30″。目前，计算垂线偏差的方法主要有四种：天文大地测量方法、重力测量方法、天文重力测量方法和gps测量方法。天文大地测量法是定义意义下的直接测量方法，需对每一个位置点长时间静态天文观测来实现，对于幅员辽阔的国家和地区来说，这种方法求解较为密集的垂线偏差虽然精度高、数据直接，但效率低、实用性差，；重力测量方法是一种间接测量方法，根据全球卫星、航空、海洋和地面重力数据进行地球模型的构建，从而得到地球上不同位置大地水准面与地球椭球模型之间的差异来求解垂线偏差。这种方法需要的数据量巨大，对于测点的位置信息精度要求较高，而计算得到的垂线偏差精度相对较低，计算效率非常高；天文重力测量是将两种方法相结合，折中了两种方法的特点；gps技术测定垂线偏差，须应用水准测量技术获取大地水准面之差或高程异常差，也称gps水准法。这种方法易受地形限制。这四种算法都无法实现垂线偏差实时计算，为惯性导航系统提供实时参数。

为此，我们研究出利用重力梯度水平张量进行垂线偏差计算的理论方法，实现垂线偏差实时和自主获取，以满足潜器长时间自主导航需求。

技术实现要素：

本发明的目的在于避免现有技术的不足，提出一种的基于水平张量重力梯度的垂线偏差实时计算方法，该算法配合重力梯度仪测量可实时输出高精度垂线偏差值，进而实现对高精度惯性导航系统的重力扰动补偿，提高导航精度。

本发明解决其技术问题是采取以下技术方案实现的：

一种基于水平张量重力梯度的垂线偏差实时计算方法，其特征在于：包括以下步骤：

步骤一，在笛卡尔坐标系统下，建立重力位场与位置偏差之间的对应关系模型；

步骤二，在标准椭球体坐标体系下，建立重力位场与位置偏差之间的对应关系模型；

步骤三，根据垂线偏差定义的数学模型，利用坐标系转换的方式推导出垂线偏差与重力梯度水平张量之间的对应关系。

而且，所述步骤一的具体过程为：

建立一个初始坐标点p0，x坐标方向指向北向，y坐标方向指向东向，而z坐标方向与p0点垂线方向一致，并指向地心；研究区域内的另外任意一点pi，pi点的z坐标方向与初始点p0的垂直方向平行，而pi点的重力方向并不与z方向一致，向量piv是i点重力gi在xz平面上的投影，而pih是i点重力gi部分分量gxi在同样平面的投影；

p0点处的天文纬度为φ0，δφi表示piv方向与pi点z轴方向的夹角，因此，pi点的纬度为：

φi＝φ0+δφi.(1)

-gxi＝gisinδφi，认为δφi＝sinδφi，则有：

在yz平面的天文经度的变化的具体表达式为：

方程(2)和(3)表明大地水准面上p0和pi点处产生的北向和东向分量，同样在任意pk点处得到类似的δφk和δλk表达式，通过两个点之间的经纬度关系，得出：

同时，

其中，φi和φk为pi和pk点的纬度值，类似于方程(2)和方程(3)，方程(4)和方程(5)中的δφi、δλi和δφk、δλk分别在pi和pk点处的北向和东向分量。

而且，用w表示地球真实重力位场数据：和方程(4)和(5)简化表述为：

而且，所述步骤二的具体过程为：

φ和λ表示标准椭球体下的大地纬度和和大地经度，u为标准椭球体下重力位，与方程(6)和(7)类似的表述方式在标准椭球体下将变化为：

其中，γi和γk和是在pi和pk点处的两个标准椭球绝对重力值，利用方程(6)减去方程(8)和利用方程(7)减去(9)得到，

而且，所述步骤三的具体过程为：

大地纬度和大地经度等于标准椭球体下的经纬度，其线性对应关系为：φ＝nφ-κ和λ＝nλ，κ值代表测量点地球表面和标准椭球表面沿着铅垂线方向的差异；

垂线偏差定义的数学表达为：

其中，φ和λ是天文坐标系下的经纬度定义，φ和λ是大地坐标系下的经纬度定义，将方程(11)和(12)带入方程(10)，则公式(10)变化为：

其中，方程(13)和(14)中ξi和ξk、ηi和ηk表示pi和pk点处的垂线偏差；

在0.5°ⅹ0.5°范围内，两个坐标系下的正常重力值相等，因此gk＝γk，gi＝γi，方程(13)和(14)则变化为：

在从公式(15)和(16)中，以重力在x和y方向的分量作为中间变量，推导出垂线偏差与重力梯度之间的对应关系，引入新的表达式：

δw＝w-u(17)

得到如下简化公式：

重力梯度的水平分量wδ和wxy是地球的曲率效应引起的，它们与垂线偏差之间的关系直接，其具体的表达式为：

对公式(18)和(19)进行转化，建立并引入一个新的坐标系反映重力梯度水平分量与垂线偏差变化之间的对应关系，从而实现垂线偏差的实时计算，

根据新的坐标参考系统进行坐标转换，以任意点i和k连线方向作为n轴，与其垂线的方向轴作为s轴，将这个两个观测点的数据转换到新的坐标体系中，具体的变换方式如下：

则相应的坐标转换结果为

则新坐标系下重力梯度值为：

在新坐标系下建立的新的重力梯度水平张量，其与原来坐标系下的重力梯度水平张量之间的关系时建立在两个测量点i和k连续的方位角条件下进行的，因此利用这种坐标转换方式作为中间变量环节，来获取新的方程表达式；

目前，需要对方程(22)左侧进行积分，则有

如果点pi和pk之间的距离很近的话，将二阶导数wns变化看做是一个线性的，因此方程(23)的积分表达则变为：

其中，nik＝nk-ni来代表点pi和pk之间的距离；

根据(23)的积分结果以及方程(21)进行有效地合并转换则有，

wsk-wsi＝-(wxk-wxi)sinαik+(wyk-wyi)cosαik(25)

类似的情况，在大地坐标系下的展开公式为：

usk-usi＝-(uxk-uxi)sinαik+(uyk-uyi)cosαik(26)

利用新的变量公式来代表方程(25)减去(26)的差，gik/g＝δυik，并且表示为两个点pi和pk在n方向上的水平变量，得到如下表达式：

根据垂线偏差的分析结果则有：

gik/g＝(ξk-ξi)sinαik-(ηk-ηi)cosαik(28)

而相应的δwns可以计算为

其中，δwδ＝wδ-uδ和δwxy＝wxy-uxy，其中，wδ和wxy可以通过重力梯度仪进行测量，而uδ和uxy是大地坐标系下的重力梯度值可以利用公式进行求解，将方程(30)带入方程(29)后则有

将方程(31)带入方程(28)，得到最终垂线偏差与重力梯度水平张量之间的对应关系：

从公式(32)中，能够准确地分析出各个测量点的重力梯度水平分量值与观测量点垂线偏差之间的对应关系，利用该公式实现垂线偏差实现计算。

本发明的优点和积极效果是：

本方法建立重力梯度水平张量与垂线偏差之间的对应关系，利用重力梯度仪实时测量的重力梯度水平张量数据实时获取垂线偏差，较常规垂线偏差测量方法，本方法可以实现垂线偏差实时准确计算，并无需依赖天文和卫星导航数据，具备自主测量能力。

附图说明

图1为任意观测点在天文坐标系统之间的对应关系；

图2为x和y向n和s坐标转换形成新的坐标系统。

具体实施方式

下面结合附图并通过具体实施例对本发明作进一步详述，以下实施例只是描述性的，不是限定性的，不能以此限定本发明的保护范围。

一种基于水平张量重力梯度的垂线偏差实时计算方法，建立不同坐标体系下重力梯度与经纬度之间的线性关系，利用2～3个已知点位的垂线偏差作为基准点进行垂线偏差迭代求解。

具体过程如下：

根据垂线偏差定义可知，垂线偏差的实质就是不同坐标系(笛卡尔坐标和标准椭球坐标系)下的法线方向的差异。

在笛卡尔坐标系统下，假定任意一个初始坐标点p0，其中，x坐标方向指向北向，y坐标方向指向东向，而z坐标方向与p0点垂线方向一致，并指向地心。因此，在研究区域有另外任意一点pi，它的z方向平行于p0点的z坐标轴方向，xi方向平行于p0点的天文子午面的切线，这两个点之间的具体关系如图1所示。假设，pi点的z坐标方向与初始点p0的垂直方向平行，而pi点的重力方向并不与z方向一致。在图1中，向量piv是i点重力gi在xz平面上的投影，而pih是i点重力gi部分分量gxi在同样平面的投影。

p0点处的天文纬度为φ0，δφi表示piv方向与pi点z轴方向的夹角，因此，pi点的纬度为：

φi＝φ0+δφi.(1)

根据图1所示，-gxi＝gisinδφi；实际情况下的δφi非常小，近似认为δφi＝sinδφi，则有：

在y方向也具有同样的特性，在yz平面的天文经度的变化的具体表达式为：

方程(2)和(3)表明大地水准面上p0和pi点处产生的北向和东向分量，同样在任意pk点处得到类似的δφk和δλk表达式，通过两个点之间的经纬度关系，得出：

同时，

其中，φi和φk为pi和pk点的纬度值，类似于方程(2)和方程(3)，方程(4)和方程(5)中的δφi、δλi和δφk、δλk分别在pi和pk点处的北向和东向分量。

为了更简单明确地表述这些公式，用w表示地球真实重力位场数据，然后给出如下定义：和那么方程(4)和(5)可以用以下公式进行表述：

在标准椭球体坐标体系下，我们建立同样类似的重力位场与位置偏差之间的对应关系，其中，φ和λ表示标准椭球体下的大地纬度和和大地经度，u为标准椭球体下重力位，与方程(6)和(7)类似的表述方式在标准椭球体下将变化为：

其中，γi和γk和是在pi和pk点处的两个标准椭球绝对重力值，利用方程(6)减去方程(8)和利用方程(7)减去(9)得到，

根据垂线偏差的定义：垂线偏差就是大地水准面坐标系统下的法线方向与标准地球椭球体坐标体系下的法线方向的偏差，此处近似地认为大地纬度和大地经度等于标准椭球体下的经纬度，其实它们之间存在一定的线性对应关系：φ＝nφ-κ和λ＝nλ，其中的主要区别就是κ值，它代表测量点地球表面和标准椭球表面沿着铅垂线方向的差异。在实际应用中，并不区别这两个概念，这个值一般情况下可被忽略。

垂线偏差定义的数学表达为：

其中，φ和λ是天文坐标系下的经纬度定义，φ和λ是大地坐标系下的经纬度定义，将方程(11)和(12)带入方程(10)，则公式(10)变化为：

其中，方程(13)和(14)中ξi和ξk、ηi和ηk表示pi和pk点处的垂线偏差。

查阅各类文献资料显示，0.5°ⅹ0.5°范围内，可认为两个坐标系下的正常重力值近似相等，因此gk＝γk，gi＝γi。方程(13)和(14)则变化为：

在从公式(15)和(16)中，可以看出垂线偏差与重力位的x和y方向导数相关。目前，并没有仪器可以将这两个分量准确测量出来，因此，只能以重力在x和y方向的分量作为中间变量，推导出垂线偏差与重力梯度之间的对应关系，引入新的表达式：

δw＝w-u(17)

得到如下简化公式：

重力梯度的水平分量wδ和wxy是地球的曲率效应引起的，它们与垂线偏差之间的关系直接，其具体的表达式为：

为了找出重力梯度水平张量与垂线偏差之间的内在联系，对公式(18)和(19)进行转化，建立并引入一个新的坐标系反映重力梯度水平分量与垂线偏差变化之间的对应关系，从而实现垂线偏差的实时计算，引入的的新坐标系如图2所示。

根据新的坐标参考系统进行坐标转换，以任意点i和k连线方向作为n轴，与其垂线的方向轴作为s轴，将这个两个观测点的数据转换到新的坐标体系中，具体的变换方式如下：

则相应的坐标转换结果为

则新坐标系下重力梯度值为：

在新坐标系下建立的新的重力梯度水平张量，其与原来坐标系下的重力梯度水平张量之间的关系时建立在两个测量点i和k连续的方位角条件下进行的，因此利用这种坐标转换方式作为中间变量环节，来获取新的方程表达式。目前，需要对方程(22)左侧进行积分，则有

如果点pi和pk之间的距离很近的话，我们可以近似将二阶导数wns变化看做是一个线性的，因此方程(23)的积分表达则变为：

其中，nik＝nk-ni来代表点pi和pk之间的距离。根据(23)的积分结果以及方程(21)进行有效地合并转换则有，

wsk-wsi＝-(wxk-wxi)sinαik+(wyk-wyi)cosαik(25)

类似的情况，在大地坐标系下的展开公式为：

usk-usi＝-(uxk-uxi)sinαik+(uyk-uyi)cosαik(26)

利用新的变量公式来代表方程(25)减去(26)的差，gik/g＝δυik，并且表示为两个点pi和pk在n方向上的水平变量。则可以得到如下表达式：

根据垂线偏差的分析结果则有：

gik/g＝(ξk-ξi)sinαik-(ηk-ηi)cosαik(28)

而相应的δwns可以计算为

其中，δwδ＝wδ-uδ和δwxy＝wxy-uxy。其中，wδ和wxy可以通过重力梯度仪进行测量，而uδ和uxy是大地坐标系下的重力梯度值可以利用公式进行求解。则将方程(30)带入方程(29)后则有

将方程(31)带入方程(28)则可以得到最终垂线偏差与重力梯度水平张量之间的对应关系：

利用坐标系转换的方式推导出垂线偏差与重力梯度水平张量之间的对应关系，从公式(32)中，我们可以清晰准确地分析出各个测量点的重力梯度水平分量值与观测量点垂线偏差之间的对应关系，利用该公式可以实现垂线偏差实现计算。

尽管为说明目的公开了本发明的实施例和附图，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换、变化和修改都是可能的，因此，本发明的范围不局限于实施例和附图所公开的内容。