0.   引　言

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

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

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

式中：F 为两个物体之间万有引力的大小；G为万有引力常量；$ {m}_{1} $、$m_2 $为两个物体的质量；r是两个物体质心的距离.

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

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

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

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

通过测量和计算可知，通常情况下，第一类力产生的加速度，日变化的量级为10⁻⁸ m/s²；第二类力产生的加速度，日变化的量级为10⁻⁶ m/s^2[9].

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

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

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

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

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

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

1) 设定自身坐标为$ \begin{bmatrix}x\\ y\\ z\end{bmatrix}$， 是需要求解的未知数；

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

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

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

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

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

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

8) 根据惯性导航或者精密磁力计得到当前的概略位置$ \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $，如果没有当前概略位置，则可将当前概略位置设为地心$ \begin{bmatrix}0\\ 0\\ 0\end{bmatrix} $，将未知数中的坐标初始化为概略位置，此即为初始解；

9) 根据自身坐标初始解$ \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $，和天文星历，由$ \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $以及观测时刻t_i，可以求得各个时刻的天体坐标，(i=1，···，k，k+1，k≥3取整数)，可以求得在$ \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $处，物体所受在各个时刻的合力方向和大小的近似值A_i，V_i，(i=1，···，k，k+1，k≥3 取整数)；

10) 计算4个时刻(或以上)的万有引力(同时也是加速度)的大小和方向并分别求差值，即重力加速度的大小变化估测数据和重力加速度的方向变化估测数据，然后基于步骤3)~6)得到的重力加速度的大小变化观测数据和重力加速度的方向变化观测数据，根据设定的误差限，进行迭代直至收敛到误差小于误差限，即可求得$ \begin{bmatrix}x\\ y\\ z\end{bmatrix}$. 具体过程实现如下：

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

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

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

d) 步骤a)均是计算以$ \begin{bmatrix}x\\ y\\ z\end{bmatrix} $和时间为自变量的函数，其中时间精确已知，而步骤b)和c)的结果基于步骤a)，因此ΔA_i 和ΔV_i 均为$ \begin{bmatrix}x\\ y\\ z\end{bmatrix} $的函数，故有如下的方程组：

以及

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

只有当前解为真实位置时，残差$ \omega i\mathrm{和}\lambda i $为0，(i=1，···k，k≥3取整数)，否则残差不为0.

f₁是加速度方向的残差$ \omega i $关于$ \begin{bmatrix}x\\ y\\ z\end{bmatrix} $的函数，而f₂是加速度的大小的残差$ \lambda i $ 关于$ \begin{bmatrix}x\\ y\\ z\end{bmatrix} $的函数.

将式(2)与式(3)在$ \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $处进行一阶泰勒级数展开以线性化后，可以得到：

式中，Gv和Ga分别对应于f₁和f₂的雅各比矩阵，而$ \Delta \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $是当前解与上一次解的(第一次迭代则为初始解)坐标变化量，$ \Delta \left(\Delta V\right) $为当前解与上一次解的(第一次迭代则为初始解)加速度方向值的变化量，$ \Delta \left(\Delta A\right) $为当前解与上一次解的(第一次迭代则为初始解)加速度大小的变化量.

将$ \omega i $作为$ \Delta \left(\Delta V\right) $，代入式(4)，或者将$ \lambda i $作为$ \Delta \left(\Delta A\right) $，代入式(5)，单独解这两个方程或者将这两个方程联立，均可以得到$ \Delta \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $.

而当前解$ \begin{bmatrix}{x}_{n}\\ {y}_{n}\\ {z}_{n}\end{bmatrix}=\begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix}+\Delta \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $，此时已经将位置解进行了更新，$ \begin{bmatrix}{x}_{n}\\ {y}_{n}\\ {z}_{n}\end{bmatrix} $比$ \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $ 更接近正确位置.

如此迭代，直到$ \Delta \begin{bmatrix}{x}_{0}\\ {y}_{0}\\ {z}_{0}\end{bmatrix} $小于预定限值，则结束迭代，将当前解作为最终解，完成定位.

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

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

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

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

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

Figure  1.  The magnitude of gravitational acceleration during a day

下载: 全尺寸图片 幻灯片

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

Figure  2.  The direction of gravitational acceleration during a day

下载: 全尺寸图片 幻灯片

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

Figure  3.  Relationship between observation times and horizontal positioning error

下载: 全尺寸图片 幻灯片

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

Figure  4.  Relationship between observation times and vertical positioning error

下载: 全尺寸图片 幻灯片

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

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

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

Figure  5.  Horizontal positioning error of 4 observations at different intervals

下载: 全尺寸图片 幻灯片

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

Figure  6.  Vertical positioning error of 4 observations at different intervals

下载: 全尺寸图片 幻灯片

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

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

Figure  7.  Horizontal positioning error of 5 observations at different intervals

下载: 全尺寸图片 幻灯片

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

Figure  8.  Vertical positioning error of 5 observations at different intervals

下载: 全尺寸图片 幻灯片

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

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

Figure  9.  Horizontal positioning error of 6 observations at different intervals

下载: 全尺寸图片 幻灯片

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

Figure  10.  Vertical positioning error of 6 observations at different intervals

下载: 全尺寸图片 幻灯片

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

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

Figure  11.  Horizontal positioning error of 7 observations at different intervals

下载: 全尺寸图片 幻灯片

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

Figure  12.  Vertical positioning error of 7 observations at different intervals

下载: 全尺寸图片 幻灯片

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

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

Table  1.  The number of positioning iterations and positioning errors of seven observations at different observation intervals

间隔时间/ min	迭代 次数	水平定位误差 (NS)/(°)	水平定位误差 (EW)/(°)	竖直定位 误差/m
3	683	0.027 288 97	0.007 576 43	2 093.994
5	289	0.010 477 68	0.002 833 73	784.955
7	242	0.001 398 91	0.000 368 09	102.327
9	39	0.002 025 03	0.000 516 17	144.329
11	39	0.000 615 99	0.000 152 16	42.826
13	216	0.000 674 43	0.000 160 97	45.728
15	141	0.000 298 18	0.000 068 23	19.655
17	75	0.000 000 02	0.000 000 00	0.001
19	44	0.000 000 07	0.000 000 02	0.005
21	38	0.000 000 04	0.000 000 01	0.003

下载: 导出CSV 

| 显示表格

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

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

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

选择加入的噪声级别为0、10⁻¹⁶ 、10⁻¹⁵、10⁻¹⁴、10⁻¹³、10⁻¹²共6个级别，得到的结果如表2所示.

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

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

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

Table  2.  Comparison of influence of different levels of observation noise on positioning results

噪声 方差	迭代 次数	水平定位误差 (NS)/(°)	水平定位误差 (EW)/(°)	竖直定位 误差/m
0	118	0.000 000 17	0.000 000 03	0.010
10−16	93	0.000 007 32	0.000 001 41	0.447
10−15	120	0.000 010 25	0.000 002 64	0.698
10−14	95	0.000 154 65	0.000 030 09	9.444
10−13	98	0.003 068 97	0.000 701 79	199.458
10−12	61	0.005 795 25	0.000 992 08	351.853

下载: 导出CSV 

| 显示表格

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

Table  3.  Comparison of influence of different observation times on positioning results

观测时刻 (UTC)	迭代 次数	水平定位误差 (NS)/(°)	水平定位误差 (EW)/(°)	竖直定位 误差/m
00：00	118	0.000 000 17	0.000 000 03	0.010
02：00	54	0.000 000 02	0.000 000 00	0.001
04：00	42	0.000 000 39	0.000 000 06	0.031
06：00	44	0.000 831 12	0.000 122 98	80.545
08：00	30	0.000 000 01	0.000 000 00	0.001
10：00	28	0.000 326 02	0.000 114 48	36.780
12：00	33	0.000 005 29	0.000 002 20	0.939
14：00	61	0.000 003 48	0.000 000 61	0.780
16：00	26	0.000 123 04	0.000 005 56	24.473
18：00	25	0.000 004 89	0.000 000 28	0.890
20：00	29	0.000 003 74	0.000 000 56	0.719
22：00	29	0.000 015 17	0.000 004 78	2.512

下载: 导出CSV 

| 显示表格

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)贾剑钢老师(武汉大学测绘学院)的帮助，他们提供的数据对于验证算法很有帮助. 在此表示诚挚的感谢！