CN113310486A - 一种万有引力定位方法及系统 - Google Patents

Abstract

本发明提供一种万有引力定位方法及系统，基于天体产生的万有引力、地球产生的万有引力以及其他惯性力影响的合成，导致待定位位置的重力加速度发生变化，通过定期监测待定位位置的重力加速度变化，反推得到该待定位位置的定位结果；所述重力加速度变化，为重力加速度的方向变化和/或大小变化；根据自身坐标初始解计算重力加速度的大小变化和/或方向变化估测数据，相应结合重力加速度的大小变化和/或方向变化观测数据，进行迭代线性化求解直至收敛到误差小于误差限，得到定位结果。本发明通过定期监测待定位位置的重力加速变化实现定位，这一过程完全被动，无需发射任何信号，隐蔽性好，也不易受到干扰，便于实施。

Description

一种万有引力定位方法及系统

技术领域

本发明涉及测量定位技术领域，具体涉及一种万有引力定位技术方案。

背景技术

在某些特定场合下，常用的定位方法(例如GNSS,天文观测)会失效，例如在GNSS信号受到干扰，或者在深海中，无法进行GNSS和天文观测。而惯性导航会有累积误差，在长时间运行后误差无法接受，海底地形图、地磁导航和现有正在研究的重力导航都需要事先采集大量的数据，这一工作极其耗费时间，而且无法及时更新。因此探究特殊环境下的新型定位方法非常重要。这种需求在潜艇、海底探测器中普遍存在。

发明内容

为了解决上述现有技术中存在的不足，本发明提出了一种万有引力定位技术方案。

为了实现上述目的，本发明提出一种万有引力定位方法，基于天体产生的万有引力、地球产生的万有引力以及其他惯性力影响的合成，导致待定位位置的重力加速度发生变化，通过定期监测待定位位置的重力加速度变化，反推得到该待定位位置的定位结果；所述重力加速度变化，为重力加速度的方向变化和/或大小变化；

监测重力加速度的方向变化时，测量各个时刻t₁,t₂,…,t_k+1重力的方向θ₁,θ₂,…,θ_k+1，得到重力加速度的方向变化观测数据Δθ_i＝θ_i+1–θ_i，i＝1…k；设根据坐标

时刻、相关天体位置，能够求出各个时刻t₁,t₂,…,t_k+1各个力的合成形成的加速度的方向为V₁,V₂,…,V_k+1，则重力加速度的方向变化估测数据ΔV_i＝V_i+1-V_i，i＝1…k，k为不小于3的整数，建立如下的重力加速度的方向变化相应方程组，

其中，f1()是加速度的方向的差值ΔV_i关于

的函数；

监测重力加速度的大小变化时，测量各个时刻t₁,t₂,…,t_k+1重力的大小g ₁,g₂,…,g _k+1，得到重力加速度的大小变化观测数据Δg_i＝g_i+1–g_i，i＝1…k；设根据坐标

时刻、相关天体位置，能够求出各个时刻t₁,t₂,…,t_k+1各个力的合成形成的加速度的大小为A₁,A₂,…,A_k+1，则重力加速度的大小变化估测数据ΔA_i＝A_i+1-A_i，i＝1…k，k为不小于3的整数，建立如下的重力加速度的大小变化方程组，

其中，f2()是加速度的大小的差值ΔA_i关于

的函数；

反推得到待定位位置的定位结果实现方式为，根据自身坐标初始解计算重力加速度的大小变化和/或方向变化估测数据，相应结合重力加速度的大小变化和/或方向变化观测数据，相应地通过重力加速度的方向变化相应方程组和/或大小变化方程组，进行迭代线性化求解直至收敛到误差小于误差限，得到定位结果。

而且，监测重力加速度的方向变化采用倾斜仪测量实现。

而且，监测重力加速度的大小变化采用重力仪测量实现。

而且，定期监测待定位位置的重力加速度变化时，采用原子钟获取当前时刻。

而且，所述自身坐标初始解采用根据惯性定位设备或者地磁场定位设备得到的当前的概略位置。

而且，所述天体包括太阳和月亮。

而且，根据自身坐标初始解计算重力加速度的大小变化和/或方向变化估测数据，实现方式如下，

根据自身坐标初始解、地球质量、万有引力常数计算地球造成的加速度g_e以及其方向v_e；

根据自身坐标初始解、地球自转速度计算地球自转造成的加速度a_c及其方向v_c；

根据时刻，月球质量、月球坐标、万有引力常数、地心坐标、自身坐标初始解计算在地月系统中的加速度a_em及其方向v_em；

根据时刻、太阳质量、太阳坐标、万有引力常数、地心坐标、自身坐标初始解计算在日地系统中的加速度a_se及其方向v_se；

将以上骤得到的四个加速度矢量进行合成，形成加速度的大小和方向A_i,V_i,i＝1…k+1，k为不小于3的整数，然后计算相应差值。

另一方面，本发明还提供一种万有引力定位系统，用于实现如上任一项所述的一种万有引力定位方法。

而且，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种万有引力定位方法。

或者，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种万有引力定位方法。

本发明利用了天体产生的万有引力与地球产生的万有引力以及其他惯性力的关系不断在变化，导致待定位位置的重力加速度发生极其微弱的变化，通过定期监测该位置的重力加速度方向变化以及重力加速度大小的变化，即可反推出自身位置。这一过程完全被动，无需发射任何信号，隐蔽性好，也不易受到干扰。通过倾斜仪即可测量出重力方向的变化。

本发明方案实施简单方便，实用性强，解决了相关技术存在的实用性低及实际应用不便的问题，能够提高用户体验，具有重要的市场价值。

具体实施方式

以下结合实施例具体说明本发明的技术方案。

万有引力是一种难以被干扰的力，其特点是穿透性、稳定性都很好。本发明即是利用了天体产生的万有引力与地球产生的万有引力以及其他惯性力的关系不断在变化，导致待定位位置的重力加速度发生极其微弱的变化，力的方向就是加速度的方向，加速度大小和力的大小成正比，因此测量加速度就是测量力。而每一个位置上随着时间的变化都是不同的，通过定期监测该位置的重力加速度方向变化以及重力加速度大小的变化，即可反推出自身位置。这一过程完全被动，无需发射任何信号，隐蔽性好，也不易受到干扰。通过倾斜仪即可测量出重力方向的变化。

设在地球某处，其在地心地固坐标系下的坐标为

那么在任意时刻，其所受到的重力加速度是由这么几种力的矢量合成而形成的：地球、月球、太阳以及其他星球对其的万有引力，地球自转形成的离心力，地月系统中的离心力，日地系统中的离心力，以及其他不能忽略的物体(例如附近的重物)对其的万有引力。按照短期恒定不变与始终在变化对这些力分为两大类：

第一类，短期恒定不变的：地球形成的万有引力，地球自转形成的离心力，通常来说，附近的重物短期内也不会有很大的变化；

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

对于第二类力，虽然随着时间在变化，但是其运行规律早已经被人们熟知，只需要知道时间，就可以准确计算其大小和方向。

因此，理论上来说，只需要精确知道时间，又可以精确测量出重力的方向和大小，就可以知道自身的位置。(表面上看，这个位置解有两个，即对称分布在地球两端，但是实际上通过时间可以进一步排除掉一个不合理解)。然而由于目前还不可能知道准确的水平方向，地球的引力并不能直接由公式准确算得，也就是说无法测量出绝对的倾斜角度和知道引起重力加速度的最大因素---地球引力的大小，因此通过重力方向和大小的变化来定位更为现实。这样既可以不被水平方向所干扰，又能通过差的方法减弱甚至消除掉大多数其他因素的误差，例如长期温度漂移、变形、重力异常等。

从目前掌握的知识，地球周围有很多天体，但是因为这些天体的质量以及和地球之间的距离不同，他们对重力的影响是不同的。影响最大的是月球，其次是太阳，再其次是金星、火星、木星、水星和土星。由于目前的仪器水平尚不能测量出金星等星体对地球的万有引力的影响，因此在此暂不予考虑。待以后仪器的精度水平提高后，再进行考虑，从而进一步提高观测精度。

通过计算，有如下结果：

每间隔1秒钟，月球对地面物体造成的重力方向的影响变化量为10^-6角秒量级，而如果要准确提取其变化量，需要的倾斜仪的精度至少需要达到10^-10角秒。

每间隔1分钟，月球对地面物体造成的重力方向的影响变化量为10^-5角秒量级，而如果要准确提取其变化量，需要的倾斜仪的精度至少需要达到10^-7角秒。

每间隔10分钟，月球对地面物体造成的重力方向的影响变化量为10^-4角秒量级，而如果要准确提取其变化量，需要的倾斜仪的精度至少需要达到10^-5角秒。

每间隔60分钟，月球对地面物体造成的重力方向的影响变化量为10^-3角秒量级，而如果要准确提取其变化量，需要的倾斜仪的精度至少需要达到10^-4角秒。

但是如果要区分所在位置100米差异对重力方向造成的影响变化量，倾斜仪的精度需要在以上的基础上再提高两到三个数量级。

太阳造成的重力方向变化大约为月球的1/3。其他星球造成的影响更小，如果需要高精度测量，则建议考虑金星、火星、木星、水星和土星等行星的影响。

以上均考虑了地月系统、日地系统是非惯性系的问题。

目前，现有技术已经能够对这些天体的运行轨迹进行精确计算，只需要知道准确的时刻即可。而目前精密时钟技术也非常先进，现在的铯钟，精度可以轻易达到10^-14/5Day的水平，可以保证1年的误差不超过100ps。对于倾斜的观测可以轻易达到10^-5角秒的精度。基于以上设备，已经具备了根据万有引力变化计算自身位置的条件，但是测量精度和测量频率距离并不十分理想。显然，倾斜仪的精度越高，则间隔观测时间可以越短，定位的时间也就越短，定位的精度也越高，重力仪的精度提高也可以达到相同的效果。从目前的技术来说，倾斜仪精度的提高要远比提高重力仪精度容易，且成本也低得多，例如电容式倾斜仪，只需要极板面积100平方毫米左右，摆长50毫米，间距0.25毫米，其大小如同保温杯，即可达到10^-5角秒的精度。光学法倾斜仪的也在相近的水平。具体实施时，可以采用相应的定制倾斜仪实现本发明技术方案。并且可以预见，本发明所述的定位方法可能会促进倾斜仪市场产品的精度进一步提高。

本发明实施例一通过定期监测待定位位置的重力加速度的方向变化，得到该位置的定位结果；提出测量的方法实现方式为：

1)根据使用的倾斜仪的精度水平，选择不同的时间间隔，在四个时刻t₁,t₂,t₃,t₄观测倾斜仪的角度θ₁,θ₂,θ₃,θ₄，然后将相邻的两个观测值作差，可以得到三个观测值差Δθ₁,Δθ₂,Δθ₃。

Δθ₁＝θ₂-θ₁,

Δθ₂＝θ₃-θ₂,

Δθ₃＝θ₄-θ₃。

2)再假设根据坐标

时刻、星体位置，可以准确求出各个时刻前面提到的各个力的合力方向V₁,V₂,V₃,V₄。该方向为一个角度，实际是X,Y,Z三个方向的合成，主要贡献来自X方向和Y方向。当连续观测4个时刻后，将相邻的两个方向作差，即可以求得三个计算方向差ΔV₁,ΔV₂,ΔV₃，

3)如果

就是倾斜仪所在的真实的位置，那么Δθ₁,Δθ₂,Δθ₃应该分别等于ΔV₁,ΔV₂,ΔV₃。

ΔV₁＝V₂-V₁,

ΔV₂＝V₃-V₂,

ΔV₃＝V₄-V₃。

而事实上，

是需要求的未知数。由于ΔV₁,ΔV₂,ΔV₃是

以及时间的函数，因此以上的三个计算方向差的算式可以表达为：

其中，t₁,t₂,t₃,t₄均为已知的精确观测值，而ΔV₁,ΔV₂,ΔV₃是已知的观测值Δθ₁,Δθ₂,Δθ₃，f()是方向的差值关于

的函数。因此上式成为了一个有三个未知数的方程组，包含了三个公式，因此不出意外的话，其有唯一解。该方程是非线性方程，可以通过泰勒展开、迭代求解。

具体实施时，为了能够缩短收敛时间，并提高定位的成功率，最好有辅助定位设备，例如惯性导航、精密磁力计等。

为了达到更精确的解，可以增加观测次数，以得到超定方程组，然后使用最小二乘等方法求解。

如果场景特殊，则需要考虑其他物体产生的引力。例如潜艇中，潜艇自身的质量巨大，其影响不可忽略，应当计入，通常这个数值是个常量。

由于根据以上条件，不仅可以计算出重力加速度的方向，也可计算出重力加速度的大小，因此，如果另有一台精密重力仪，在测量倾斜仪的同时测量重力值，同样进行两次观测值之间求差，即可得到另一组方程组，将两个方程组联合得到新的方程组，可以在一定程度上加强约束条件，从而减少误差。

以上涉及的计算方法均为成熟的方法，具体过程不再赘述。

实施例二提出的一种万有引力定位方法，通过定期监测待定位位置的重力加速度的方向变化和大小变化，或者只监测待定位位置的重力加速度的方向，得到该位置的定位结果；实现具体如下：

(为简化表述，该实施例仅考虑月球、太阳造成的引力影响。如果要考虑其它天体等因素，只需要仿照该例加入相应的条件即可。)

1，所需的设备：高精密时钟(一般用原子钟)，计算设备(可采用计算机或者其他带处理器和存储器的设备)，精密倾斜仪，精密重力仪(可选)，惯性定位设备(可选)，地磁场定位设备(可选，例如精密磁力计)；

2，从高精密时钟获取当前时刻(年月日时分秒)，根据当前时刻计算月球、太阳等天体的位置。由于最终要计算自身在地球上的位置，因此一般使用地心地固坐标系。该计算方法是既有方法，不再详述；

3，测量当前重力的方向(倾斜仪的倾斜角度)；

4，测量当前重力的大小(如果有重力仪则执行本步骤，否则直接进入步骤5并在后续流程中忽略步骤8以及步骤11中的子步骤f，而且在步骤11中的子步骤h中，得不到方程2；也即如果能测量当前重力，则可以通过解方程1与方程2的联立方程得到位置，如果不能则只通过方程1得到位置，实际上也可以只通过解方程2得到位置，但是成本会高出很多)；

5，间隔一段时间，根据现有倾斜仪的精度，建议时间间隔为1分钟；

6，重复步骤2,3,4,5一共执行至少4个循环，从而得到至少4组数据：

观测时刻：t₁,t₂,t₃,t₄…

月球坐标：

太阳坐标：

倾斜角度：θ₁,θ₂,θ₃,θ₄…

重力值：g₁,g₂,g₃,g₄…；

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

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

9，设定自身坐标为

则存在：

由

以及观测时刻t_i，月球坐标

太阳坐标

(i＝1…k,k+1，k为不小于3的整数)，可以求得在

处，每个观测时刻的地球引力大小Ae和方向Ve，月球引力大小Am和方向Vm，太阳引力大小As和方向Vs，以及地球自转、地月系统中、日地系统中的离心力的大小和方向。由这些可以进一步计算合力的大小和方向。这个合力的大小和方向应该与观测值一致，但因为

是未知数，而目前只有观测值，所以可以通过观测值求得

10，根据惯性导航或者精密磁力计得到当前的概略位置

如果没有当前概略位置，则可将当前概略位置设为地心

将未知数中的坐标初始化为概略位置，此即为初始解；具体实施时，如果设置了惯性定位设备，可以根据惯性导航得到当前的概略位置，如果设置了地磁场定位设备，可以根据精密磁力计得到当前的概略位置。

11，根据自身坐标初始解

以及步骤6中得到的天体坐标数据(月球坐标和太阳坐标)，计算4个时刻(或以上)的万有引力(同时也是加速度)的大小和方向并分别求差值(即重力加速度的大小变化估测数据和重力加速度的方向变化估测数据)，然后基于步骤7、8得到的重力加速度的大小变化观测数据和重力加速度的方向变化观测数据，根据设定的误差限，进行迭代直至收敛到误差小于误差限，即可求得

具体技术过程实现如下：

a)根据自身坐标初始解、地球质量、万有引力常数计算地球造成的加速度g_e以及其方向v_e；

b)根据自身坐标初始解、地球自转速度计算地球自转造成的加速度a_c及其方向v_c；

c)根据时刻、月球质量、月球坐标、万有引力常数、地心坐标、自身坐标初始解计算在地月系统中的加速度a_em及其方向v_em；

d)根据时刻、太阳质量、太阳坐标、万有引力常数、地心坐标、自身坐标初始解计算在日地系统中的加速度a_se及其方向v_se；

e)将以上a,b,c,d四个步骤得到的四个(或以上)加速度矢量进行合成，形成加速度的大小和方向A_i,V_i,(i＝1…k,k+1，k为不小于3的整数)；

f)计算加速度大小的差值ΔA_i＝A_i+1-A_i，(i＝1…k，k为不小于3的整数)；

g)计算加速度的方向的差值ΔV_i＝V_i+1-V_i，(i＝1…k，k为不小于3的整数)；

h)其中a,b两个步骤均是计算以

为自变量的函数，而c,d两个步骤是计算以

和时间为自变量的函数，而时间是已知量。而步骤e,f的结果是基于a,b,c,d四个步骤的，因此ΔA_i和ΔV_i均为

的函数，因此有如下的方程组：

以及：

其中，根据步骤7和步骤8所得结果，以及步骤6的观测结果，计算加速度方向

残差和加速度大小的残差：

ωi＝ΔVi-Δθi

λi＝ΔAi-Δgi

(i＝1…k，k为不小于3的整数)

只有当前解为真实位置的时候，残差ωi和λi为0，(i＝1…k，k为不小于3的整数)，否则残差不为0。

f1()是加速度的方向的残差ωi关于

的函数，而f2()是加速度的大小的残差λi关于

的函数。

将方程1与方程2在

处进行一阶泰勒级数展开以线性化后，可以得到：

和

其中，Gv和Ga分别对应于f1()和f2()的雅各比矩阵，而

是当前解与上一次的解(第一次迭代则为初始解)坐标的变化量，Δ(ΔV)为当前解与上一次的解(第一次迭代则为初始解)加速度方向值的变化量，Δ(ΔA)为当前解与上一次的解(第一次迭代则为初始解)加速度大小的变化量。

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

而当前解

此时已经将位置解进行了更新，

比

更接近正确位置。

如此迭代，直到

小于预定限值(例如三个量均小于10米，或者当前解位置与上一次的解之间的距离小于17米

则结束迭代，将当前解作为最终解，完成定位。

具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。

在一些可能的实施例中，提供一种万有引力定位系统，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种万有引力定位方法。

在一些可能的实施例中，提供一种万有引力定位系统，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种万有引力定位方法。

在一些可能的实施例中,可以使用多台倾斜仪做差分以消除某些误差，得到更精确和可靠的结果。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (10)

1.一种万有引力定位方法，其特征在于：基于天体产生的万有引力、地球产生的万有引力以及其他惯性力影响的合成，导致待定位位置的重力加速度发生变化，通过定期监测待定位位置的重力加速度变化，反推得到该待定位位置的定位结果；所述重力加速度变化，为重力加速度的方向变化和/或大小变化；

监测重力加速度的方向变化时，测量各个时刻t₁,t₂,…,t_k+1重力的方向θ₁,θ₂,…,θ_k+1，得到重力加速度的方向变化观测数据Δθ_i＝θ_i+1–θ_i，i＝1…k；设根据坐标

时刻、相关天体位置，能够求出各个时刻t₁,t₂,…,t_k+1各个力的合成形成的加速度的方向为V₁,V₂,…,V_k+1，则重力加速度的方向变化估测数据ΔV_i＝V_i+1-V_i，i＝1…k，k为不小于3的整数，建立如下的重力加速度的方向变化相应方程组，

其中，f1()是加速度的方向的差值ΔV_i关于

的函数；

监测重力加速度的大小变化时，测量各个时刻t₁,t₂,…,t_k+1重力的大小g₁,g₂,…,g_k+1，得到重力加速度的大小变化观测数据Δg_i＝g_i+1–g_i，i＝1…k；设根据坐标

时刻、相关天体位置，能够求出各个时刻t₁,t₂,…,t_k+1各个力的合成形成的加速度的大小为A₁,A₂,…,A_k+1，则重力加速度的大小变化估测数据ΔA_i＝A_i+1-A_i，i＝1…k，k为不小于3的整数，建立如下的重力加速度的大小变化方程组，

其中，f2()是加速度的大小的差值ΔA_i关于

的函数；

反推得到待定位位置的定位结果实现方式为，根据自身坐标初始解计算重力加速度的大小变化和/或方向变化估测数据，相应结合重力加速度的大小变化和/或方向变化观测数据，相应地通过重力加速度的方向变化相应方程组和/或大小变化方程组，进行迭代线性化求解直至收敛到误差小于误差限，得到定位结果。

2.根据权利要求1所述万有引力定位方法，其特征在于：监测重力加速度的方向变化采用倾斜仪测量实现。

3.根据权利要求1所述万有引力定位方法，其特征在于：监测重力加速度的大小变化采用重力仪测量实现。

4.根据权利要求1所述万有引力定位方法，其特征在于：定期监测待定位位置的重力加速度变化时，采用原子钟获取当前时刻。

5.根据权利要求1所述万有引力定位方法，其特征在于：所述自身坐标初始解采用根据惯性定位设备或者地磁场定位设备得到的当前的概略位置。

6.根据权利要求1或2或3或4或5所述万有引力定位方法，其特征在于：所述天体包括太阳和月亮。

7.根据权利要求6所述万有引力定位方法，其特征在于：根据自身坐标初始解计算重力加速度的大小变化和/或方向变化估测数据，实现方式如下，

根据自身坐标初始解、地球质量、万有引力常数计算地球造成的加速度g_e以及其方向v_e；

根据自身坐标初始解、地球自转速度计算地球自转造成的加速度a_c及其方向v_c；

根据时刻，月球质量、月球坐标、万有引力常数、地心坐标、自身坐标初始解计算在地月系统中的加速度a_em及其方向v_em；

根据时刻、太阳质量、太阳坐标、万有引力常数、地心坐标、自身坐标初始解计算在日地系统中的加速度a_se及其方向v_se；

将以上骤得到的四个加速度矢量进行合成，形成加速度的大小和方向A_i,V_i,i＝1…k+1，k为不小于3的整数，然后计算相应差值。

8.一种万有引力定位系统，其特征在于：用于实现如权利要求1-7任一项所述的一种万有引力定位方法。

9.根据权利要求8所述万有引力定位系统，其特征在于：包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如权利要求1-7任一项所述的一种万有引力定位方法。

10.根据权利要求8所述万有引力定位系统，其特征在于：包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如权利要求1-7任一项所述的一种万有引力定位方法。