基于万有引力的普适定位

张涛，邹进贵

( 武汉大学测绘学院, 武汉 430079 )

0 引言

当前，地面上开阔空间的定位手段已经十分成熟，尤其是GNSS 定位技术，无论从精度还是便捷性、可靠性方面，均可以满足大多数场合的需求. 但是，在地下、深海定位领域，GNSS 定位等大多数基于电磁波传播的定位技术难以发挥作用，因此需要采用其他技术，诸如声学定位、光学定位、地磁匹配定位、重力场匹配或者重力梯度匹配定位、惯性导航等. 其中的重力场匹配和重力梯度匹配定位技术，与本文提出的普适定位方法都是基于物理场的被动定位，但是重力场匹配和重力梯度匹配定位技术需要首先得到全球的精确重力场数据模型，然后测量所处位置的重力场或者重力梯度，才可以完成“匹配”[1-4]. 同时，这个数据模型受到地表浅层物体的影响而变化，需要及时更新，因此在具体实现上尚有诸多问题待解决. 而本文提出的基于万有引力的普适定位所测量的，是由地月和地日系统以及太阳系其他行星规律运动而形成的随时间变化的引力变化，无需提前采集重力数据，比当前的重力场匹配和重力梯度匹配定位方法更具有普适性. 从另一个角度说，现有的重力场匹配和重力梯度匹配的定位方法是基于万有引力随空间的变化进行定位[5-9]，而本文提出的基于万有引力的定位方法是根据万有引力随时间的变化进行定位. 由于万有引力随着时间的变化量非常小，而当前万有引力观测仪器(较为常见的是重力仪或者倾斜仪，实际观测到的是万有引力以及其他惯性力的合力的大小或者方向)精度有限，因此，采用本方法实现定位所需时间较长，而且定位精度也容易被观测误差影响. 本文对定位时间和精度做了具体的数值分析.

进一步的，只要有附近星体的精密星历，本文提出的定位方法还可以扩展到地外空间，故可称为“普适定位”.

任何物体都会受到万有引力的作用，而且万有引力的大小遵循如下公式：

式中：F为两个物体之间万有引力的大小；G为万有引力常量；m1、m2为两个物体的质量；r是两个物体质心的距离.

显然，万有引力的大小和双方的质量以及距离相关. 以地球上的物体为例，在其受到的诸多万有引力中，地球产生的万有引力占绝大部分，而且，在同一个地点，短期内来自地球的万有引力基本恒定；除此之外，由于太阳、月亮以及其他近地行星也产生可以观测到的引力，表现为潮汐. 距离地球太远的行星以及太阳系外天体虽然也形成万有引力，但是通常的仪器难以观测. 当前，天体的运动规律已经可以精确掌握，因此，排除重力异常等因素，在任意时刻，地球上任意地点的物体所受到的万有引力大小和方向的变化是可以通过计算得到的[8]. 反之，在一个时间段内，连续测量一个未知地点的物体所受到的万有引力的大小和方向，也可以计算其位置.

地球附近的物体，在任意时刻，其所受到的重力加速度主要是由如下几种力矢量合成：地球、月球、太阳以及其他星球对其的万有引力，地球自转形成的离心力，地月系统中的离心力，日地系统中的离心力，以及其他不能忽略的物体(例如附近的重物)对其的万有引力. 按照这些力的特点，将其分为两大类：

第一类，短期变化量较小的：包括地球形成的万有引力，地球自转形成的离心力等；

第二类，始终随时间变化的：月球、太阳以及其他星球形成的万有引力；地月系统以及日地系统中的离心力.

通过测量和计算可知，通常情况下，第一类力产生的加速度，日变化的量级为10-8m/s2；第二类力产生的加速度，日变化的量级为10-6m/s2[9].

虽然第二类力随着时间在变化， 但其运行遵循一定规律，根据时间可以准确计算其大小和方向.

可见，在短期时间内，物体所受合力的大小和方向变化、时间变化以及其所在的位置之间存在一定的函数关系. 当时间已知，物体所受的合力大小或者方向作为观测量，可以计算其所在的位置. 由于它们之间的函数关系非线性，因此所得的解不唯一，需要利用其他方式排除错误解. 进一步地，为了减少误差，需要将物体所受合力大小或者方向的变化(按照时间间隔观测量之间的差)作为已知量. 这是基于万有引力进行普适定位的基础.

1 通过万有引力进行普适定位的方法和步骤

地球周围天体的质量以及距地球的距离不同，他们对重力的影响也不同. 影响最大的是月球，其次是太阳，再次是金星、火星、木星、水星和土星[2].

现有技术在只需准确时刻的条件下，已经能够精确计算这些天体的运行轨迹. 而且目前精密时钟技术也非常先进，铯钟精度很容易达到10-14/d 的水平，可以保证1 年的误差不超过100 ps. 现代的倾斜仪对于倾斜的观测可以轻易达到0.0002 arcsec 的精度[10]，精密相对重力仪的精度更是达到了10-11m/s2[11]. 随着技术的进步，这些仪器的精度仍在不断提高[6]. 由于天体造成的万有引力变化随时间的变化量极其微小，因此，观测仪器的精度越高，越有可能在短期内获取定位结果. 虽然使用精度低的仪器进行长时间观测能够测得万有引力的变化，但是在观测量中也引入了其他干扰项，会对定位结果造成不利影响.

以下是根据倾斜仪或者重力仪的间隔观测量来得到当前位置的具体步骤：

2) 使用倾斜仪测量当前重力的方向，或者使用重力仪测量重力的大小，两者只需其一即可，如果两者都有，可以增加约束条件，减少误差；

3) 间隔一段时间，具体时长取决于仪器的精度，如果精度高则可以间隔短些，否则测量不出引力的变化量. 但是间隔时间长的话，诸如重力异常之类的非规律性引力变化幅度变大，会严重影响定位结果，文章后面会有定量分析；

4) 从精密时钟获取观测时刻，计算观测时刻的月球、太阳等天体位置；

5) 由于定位需要求的未知数是三维坐标，而且需要根据观测量的差进行计算定位，因此重复步骤1)～3)，执行至少4 个循环，为了能组成超定方程增加约束条件，应适当增加观测次数，从而得到至少4 个不同时刻的数据：观测时刻(t1，t2，t3，t4，···)，以及该时刻对应的倾斜角度(θ1，θ2，θ3，θ4，···)，或者重力值(g1，g2，g3，g4，···)；

6) 由于经过了航行、晃动等过程，不能保证倾斜仪观测到的角度是绝对倾斜角度，因此需要根据倾斜角度(θ1，θ2，θ3，θ4，···)计算倾斜角度的变化量Δθi=θi+1-θi作为观测数据，即重力加速度的方向变化观测数据，(i=1，···，k，k≥3 取整数)，从而得到至少3 个倾斜角度变化量Δθ1，Δθ2，Δθ3；

7) 如果有重力数据，计算重力差值：Δgi=gi+1-gi作为观测数据，(i=1,…,k，k≥3 取整数)，即重力加速度的大小变化观测数据，如果是相对重力仪则可以省略执行步骤4)，Δgi由相对重力仪装置直接测量得到；

a) 根据自身坐标初始解和时间计算加速度的大小和方向Ai，Vi，(i=1，···，k，k+1，k≥3 取整数)；

b) 计算加速度大小的差值ΔAi=Ai+1-Ai，(i=1，···，k，k≥3 取整数)；

c) 计算加速度的方向的差值ΔVi=Vi+1-Vi，(i=1，···，k，k≥3 取整数)；

以及

其中，根据步骤9)所得结果，以及步骤6)的观测结果，计算加速度方向残差或者加速度大小的残差：

只有当前解为真实位置时，残差 ωi和λi为0，(i=1，···k，k≥3 取整数)，否则残差不为0.

f1是加速度方向的残差 ωi关于的函数，而f2是加速度的大小的残差 λi关于的函数.

式中，Gv和Ga分别对应于f1和f2的雅各比矩阵，而是当前解与上一次解的(第一次迭代则为初始解)坐标变化量， Δ(ΔV) 为当前解与上一次解的(第一次迭代则为初始解)加速度方向值的变化量，Δ(ΔA)为当前解与上一次解的(第一次迭代则为初始解)加速度大小的变化量.

将 ωi作为 Δ(ΔV) ，代入式(4)，或者将 λi作为Δ(ΔA)，代入式(5)，单独解这两个方程或者将这两个方程联立，均可以得到.

2 通过万有引力进行普适定位的正演和反演模拟

根据上述步骤，采用生成模拟重力观测值的方法，验证根据万有引力进行定位的可行性. 假设只使用重力仪进行观测. 定位结果以经纬度和高程体现.

以武汉市的概略坐标(30.5°N ， 114.3°E ，海拔高度为30 m)作为模拟实验地点，模拟时间为2022 年12 月1 日(UTC)，观测间隔为1 min[12]，可得该处的全天重力加速度的大小和方向变化曲线分别如图1和图2 所示.

图1 一天内重力加速度大小的变化

图2 一天内重力加速度方向的变化

根据第2 节描述的方法，按照间隔1 min 观测，分别取开始的15～60 个观测量进行反演，定位结果归化为经纬度和高程. 定位结果精度如图3 和图4 所示.

图3 观测次数和水平定位误差的关系

图4 观测次数和竖直定位误差的关系

可以看出，观测间隔为1 min，即使观测量达到60 个时，定位精度也仅达到数百米，此时的方程数量庞大，迭代时间较长. 因此，如果观测间隔时间过短，只通过增加观测量来提高定位精度的话，效果并不明显. 下面进行固定观测次数、改变观测间隔的方法进行测试.

首先使用4 个观测量，改变间隔时间从11～59 min，每2 min 递增一次，结果如图5～6 所示.

图5 4 个观测量，不同间隔时间的水平定位误差

图6 4 个观测量，不同间隔时间的竖直定位误差

观测5 次，间隔从5 min 开始，每次增加2 min，结果如图7～8 所示.

图7 5 个观测量，不同观测间隔时间的水平定位误差

图8 5 个观测量，不同观测间隔时间的竖直定位误差

观测6 次，间隔从5 min 开始，每次增加2 min，结果如图9～10 所示.

图9 6 个观测量，不同观测间隔时间的水平定位误差

图10 6 个观测量，不同观测间隔时间的水平定位误差

观测7 次，间隔从3 min 开始，每次增加2 min，结果如图11～12 所示，为了能放大纵轴，看清精度，图11 从间隔7 min 开始绘制.

图11 7 个观测量，不同观测间隔时间的水平定位误差

图12 7 个观测量，不同观测间隔时间的水平定位误差

为节省篇幅，如表1 所示，除了提供7 次观测量的数据表格外，其余只留曲线图.

表1 7 个观测量，不同观测间隔时间的定位迭代次数和定位误差

由表1 可知，当前使用基于万有引力的定位算法，在没有任何误差参与的情况下，当观测总时长达到20 min 时，定位精度约为1 km. 要想得到亚米级的定位结果，观测时间至少需要1 h；当观测时间达到2 h 时，结果会比较稳定，维持在厘米级.

实际情况下，观测量会有各种误差，包括仪器自身的误差，外界的干扰. 下面验证误差和定位精度的关系.

采用7 个观测量，观测间隔为20 min. 对模拟观测量加入不同等级的高斯噪声，然后验证定位结果.

选择加入的噪声级别为0、10-16、10-15、10-14、10-13、10-12共6 个级别，得到的结果如表2 所示.

表2 不同级别的观测噪声对定位结果的影响对比

可见本文于误差的容忍度仍然太低，定位算法比较脆弱.

此外，由于基于万有引力的普适定位方法是根据不同事件的重力观测值的差来计算位置，由图1 可知，不同时间的重力观测值的变化率是不一样的，因此在不同时间定位的精度也不尽相同. 例如，同样是7 个观测量，观测间隔为20 min，在一天中不同时刻的定位结果误差如表3 所示.

表3 2022 年12 月1 日不同观测时刻对定位结果的影响对比

3 总结与展望

通过上述分析和计算证明：通过万有引力进行普适定位是存在可行性的. 在理想状况下，当总观测时长达到20 min 的时候，定位精度约为1 km，在观测总时长达到2 h 时，定位精度可以达到厘米级. 但是目前本文的算法对测量误差极为敏感. 下一步的工作是对定位方法进行优化和完善，提高其稳健度，并结合高精度仪器进行实际测量，进行实际定位实验，进一步对定位方法进行验证. 如果证实可行，这种基于万有引力的普适定位方法可以应用于地下、水下定位，由于完全被动无源，而且万有引力难以被干扰，因此安全性高. 进一步的，这种定位方法还可以拓展应用到其他星球上的定位. 另外，由于误差传播放大的作用，有可能对引力传播速度进行验证.

致谢：在研究过程中，得到了吴立恒老师(应急管理部国家自然灾害防治研究院)、邢乐林老师(中国科学院精密测量科学与技术创新研究院)、张燕老师(湖北省地震局)、Clark R Wilson 教授(Professor Emeritus, Department of Geological Sciences, Jackson School of Geosciences，The University of Texas at Austin， Research Professor, Center for Space Research, Cockrell School of Engineering， Dave P.Carlton Centennial Professorship in Geophysics)贾剑钢老师(武汉大学测绘学院)的帮助，他们提供的数据对于验证算法很有帮助. 在此表示诚挚的感谢！