基于萬有引力的普適定位

張濤，鄒進貴

( 武漢大學測繪學院, 武漢 430079 )

0 引言

當前，地面上開闊空間的定位手段已經(jīng)十分成熟，尤其是GNSS 定位技術，無論從精度還是便捷性、可靠性方面，均可以滿足大多數(shù)場合的需求. 但是，在地下、深海定位領域，GNSS 定位等大多數(shù)基于電磁波傳播的定位技術難以發(fā)揮作用，因此需要采用其他技術，諸如聲學定位、光學定位、地磁匹配定位、重力場匹配或者重力梯度匹配定位、慣性導航等. 其中的重力場匹配和重力梯度匹配定位技術，與本文提出的普適定位方法都是基于物理場的被動定位，但是重力場匹配和重力梯度匹配定位技術需要首先得到全球的精確重力場數(shù)據(jù)模型，然后測量所處位置的重力場或者重力梯度，才可以完成“匹配”[1-4]. 同時，這個數(shù)據(jù)模型受到地表淺層物體的影響而變化，需要及時更新，因此在具體實現(xiàn)上尚有諸多問題待解決. 而本文提出的基于萬有引力的普適定位所測量的，是由地月和地日系統(tǒng)以及太陽系其他行星規(guī)律運動而形成的隨時間變化的引力變化，無需提前采集重力數(shù)據(jù)，比當前的重力場匹配和重力梯度匹配定位方法更具有普適性. 從另一個角度說，現(xiàn)有的重力場匹配和重力梯度匹配的定位方法是基于萬有引力隨空間的變化進行定位[5-9]，而本文提出的基于萬有引力的定位方法是根據(jù)萬有引力隨時間的變化進行定位. 由于萬有引力隨著時間的變化量非常小，而當前萬有引力觀測儀器(較為常見的是重力儀或者傾斜儀，實際觀測到的是萬有引力以及其他慣性力的合力的大小或者方向)精度有限，因此，采用本方法實現(xiàn)定位所需時間較長，而且定位精度也容易被觀測誤差影響. 本文對定位時間和精度做了具體的數(shù)值分析.

進一步的，只要有附近星體的精密星歷，本文提出的定位方法還可以擴展到地外空間，故可稱為“普適定位”.

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

式中：F為兩個物體之間萬有引力的大小；G為萬有引力常量；m1、m2為兩個物體的質(zhì)量；r是兩個物體質(zhì)心的距離.

顯然，萬有引力的大小和雙方的質(zhì)量以及距離相關. 以地球上的物體為例，在其受到的諸多萬有引力中，地球產(chǎn)生的萬有引力占絕大部分，而且，在同一個地點，短期內(nèi)來自地球的萬有引力基本恒定；除此之外，由于太陽、月亮以及其他近地行星也產(chǎn)生可以觀測到的引力，表現(xiàn)為潮汐. 距離地球太遠的行星以及太陽系外天體雖然也形成萬有引力，但是通常的儀器難以觀測. 當前，天體的運動規(guī)律已經(jīng)可以精確掌握，因此，排除重力異常等因素，在任意時刻，地球上任意地點的物體所受到的萬有引力大小和方向的變化是可以通過計算得到的[8]. 反之，在一個時間段內(nèi)，連續(xù)測量一個未知地點的物體所受到的萬有引力的大小和方向，也可以計算其位置.

地球附近的物體，在任意時刻，其所受到的重力加速度主要是由如下幾種力矢量合成：地球、月球、太陽以及其他星球?qū)ζ涞娜f有引力，地球自轉(zhuǎn)形成的離心力，地月系統(tǒng)中的離心力，日地系統(tǒng)中的離心力，以及其他不能忽略的物體(例如附近的重物)對其的萬有引力. 按照這些力的特點，將其分為兩大類：

第一類，短期變化量較小的：包括地球形成的萬有引力，地球自轉(zhuǎn)形成的離心力等；

第二類，始終隨時間變化的：月球、太陽以及其他星球形成的萬有引力；地月系統(tǒng)以及日地系統(tǒng)中的離心力.

通過測量和計算可知，通常情況下，第一類力產(chǎn)生的加速度，日變化的量級為10-8m/s2；第二類力產(chǎn)生的加速度，日變化的量級為10-6m/s2[9].

雖然第二類力隨著時間在變化， 但其運行遵循一定規(guī)律，根據(jù)時間可以準確計算其大小和方向.

可見，在短期時間內(nèi)，物體所受合力的大小和方向變化、時間變化以及其所在的位置之間存在一定的函數(shù)關系. 當時間已知，物體所受的合力大小或者方向作為觀測量，可以計算其所在的位置. 由于它們之間的函數(shù)關系非線性，因此所得的解不唯一，需要利用其他方式排除錯誤解. 進一步地，為了減少誤差，需要將物體所受合力大小或者方向的變化(按照時間間隔觀測量之間的差)作為已知量. 這是基于萬有引力進行普適定位的基礎.

1 通過萬有引力進行普適定位的方法和步驟

地球周圍天體的質(zhì)量以及距地球的距離不同，他們對重力的影響也不同. 影響最大的是月球，其次是太陽，再次是金星、火星、木星、水星和土星[2].

現(xiàn)有技術在只需準確時刻的條件下，已經(jīng)能夠精確計算這些天體的運行軌跡. 而且目前精密時鐘技術也非常先進，銫鐘精度很容易達到10-14/d 的水平，可以保證1 年的誤差不超過100 ps. 現(xiàn)代的傾斜儀對于傾斜的觀測可以輕易達到0.0002 arcsec 的精度[10]，精密相對重力儀的精度更是達到了10-11m/s2[11]. 隨著技術的進步，這些儀器的精度仍在不斷提高[6]. 由于天體造成的萬有引力變化隨時間的變化量極其微小，因此，觀測儀器的精度越高，越有可能在短期內(nèi)獲取定位結果. 雖然使用精度低的儀器進行長時間觀測能夠測得萬有引力的變化，但是在觀測量中也引入了其他干擾項，會對定位結果造成不利影響.

以下是根據(jù)傾斜儀或者重力儀的間隔觀測量來得到當前位置的具體步驟：

2) 使用傾斜儀測量當前重力的方向，或者使用重力儀測量重力的大小，兩者只需其一即可，如果兩者都有，可以增加約束條件，減少誤差；

3) 間隔一段時間，具體時長取決于儀器的精度，如果精度高則可以間隔短些，否則測量不出引力的變化量. 但是間隔時間長的話，諸如重力異常之類的非規(guī)律性引力變化幅度變大，會嚴重影響定位結果，文章后面會有定量分析；

4) 從精密時鐘獲取觀測時刻，計算觀測時刻的月球、太陽等天體位置；

5) 由于定位需要求的未知數(shù)是三維坐標，而且需要根據(jù)觀測量的差進行計算定位，因此重復步驟1)～3)，執(zhí)行至少4 個循環(huán)，為了能組成超定方程增加約束條件，應適當增加觀測次數(shù)，從而得到至少4 個不同時刻的數(shù)據(jù)：觀測時刻(t1，t2，t3，t4，···)，以及該時刻對應的傾斜角度(θ1，θ2，θ3，θ4，···)，或者重力值(g1，g2，g3，g4，···)；

6) 由于經(jīng)過了航行、晃動等過程，不能保證傾斜儀觀測到的角度是絕對傾斜角度，因此需要根據(jù)傾斜角度(θ1，θ2，θ3，θ4，···)計算傾斜角度的變化量Δθi=θi+1-θi作為觀測數(shù)據(jù)，即重力加速度的方向變化觀測數(shù)據(jù)，(i=1，···，k，k≥3 取整數(shù))，從而得到至少3 個傾斜角度變化量Δθ1，Δθ2，Δθ3；

7) 如果有重力數(shù)據(jù)，計算重力差值：Δgi=gi+1-gi作為觀測數(shù)據(jù)，(i=1,…,k，k≥3 取整數(shù))，即重力加速度的大小變化觀測數(shù)據(jù)，如果是相對重力儀則可以省略執(zhí)行步驟4)，Δgi由相對重力儀裝置直接測量得到；

a) 根據(jù)自身坐標初始解和時間計算加速度的大小和方向Ai，Vi，(i=1，···，k，k+1，k≥3 取整數(shù))；

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

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

以及

其中，根據(jù)步驟9)所得結果，以及步驟6)的觀測結果，計算加速度方向殘差或者加速度大小的殘差：

只有當前解為真實位置時，殘差 ωi和λi為0，(i=1，···k，k≥3 取整數(shù))，否則殘差不為0.

f1是加速度方向的殘差 ωi關于的函數(shù)，而f2是加速度的大小的殘差 λi關于的函數(shù).

式中，Gv和Ga分別對應于f1和f2的雅各比矩陣，而是當前解與上一次解的(第一次迭代則為初始解)坐標變化量， Δ(ΔV) 為當前解與上一次解的(第一次迭代則為初始解)加速度方向值的變化量，Δ(ΔA)為當前解與上一次解的(第一次迭代則為初始解)加速度大小的變化量.

將 ωi作為 Δ(ΔV) ，代入式(4)，或者將 λi作為Δ(ΔA)，代入式(5)，單獨解這兩個方程或者將這兩個方程聯(lián)立，均可以得到.

2 通過萬有引力進行普適定位的正演和反演模擬

根據(jù)上述步驟，采用生成模擬重力觀測值的方法，驗證根據(jù)萬有引力進行定位的可行性. 假設只使用重力儀進行觀測. 定位結果以經(jīng)緯度和高程體現(xiàn).

以武漢市的概略坐標(30.5°N ， 114.3°E ，海拔高度為30 m)作為模擬實驗地點，模擬時間為2022 年12 月1 日(UTC)，觀測間隔為1 min[12]，可得該處的全天重力加速度的大小和方向變化曲線分別如圖1和圖2 所示.

圖1 一天內(nèi)重力加速度大小的變化

圖2 一天內(nèi)重力加速度方向的變化

根據(jù)第2 節(jié)描述的方法，按照間隔1 min 觀測，分別取開始的15～60 個觀測量進行反演，定位結果歸化為經(jīng)緯度和高程. 定位結果精度如圖3 和圖4 所示.

圖3 觀測次數(shù)和水平定位誤差的關系

圖4 觀測次數(shù)和豎直定位誤差的關系

可以看出，觀測間隔為1 min，即使觀測量達到60 個時，定位精度也僅達到數(shù)百米，此時的方程數(shù)量龐大，迭代時間較長. 因此，如果觀測間隔時間過短，只通過增加觀測量來提高定位精度的話，效果并不明顯. 下面進行固定觀測次數(shù)、改變觀測間隔的方法進行測試.

首先使用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 個觀測量，不同觀測間隔時間的水平定位誤差

為節(jié)省篇幅，如表1 所示，除了提供7 次觀測量的數(shù)據(jù)表格外，其余只留曲線圖.

表1 7 個觀測量，不同觀測間隔時間的定位迭代次數(shù)和定位誤差

由表1 可知，當前使用基于萬有引力的定位算法，在沒有任何誤差參與的情況下，當觀測總時長達到20 min 時，定位精度約為1 km. 要想得到亞米級的定位結果，觀測時間至少需要1 h；當觀測時間達到2 h 時，結果會比較穩(wěn)定，維持在厘米級.

實際情況下，觀測量會有各種誤差，包括儀器自身的誤差，外界的干擾. 下面驗證誤差和定位精度的關系.

采用7 個觀測量，觀測間隔為20 min. 對模擬觀測量加入不同等級的高斯噪聲，然后驗證定位結果.

選擇加入的噪聲級別為0、10-16、10-15、10-14、10-13、10-12共6 個級別，得到的結果如表2 所示.

表2 不同級別的觀測噪聲對定位結果的影響對比

可見本文于誤差的容忍度仍然太低，定位算法比較脆弱.

此外，由于基于萬有引力的普適定位方法是根據(jù)不同事件的重力觀測值的差來計算位置，由圖1 可知，不同時間的重力觀測值的變化率是不一樣的，因此在不同時間定位的精度也不盡相同. 例如，同樣是7 個觀測量，觀測間隔為20 min，在一天中不同時刻的定位結果誤差如表3 所示.

表3 2022 年12 月1 日不同觀測時刻對定位結果的影響對比

3 總結與展望

通過上述分析和計算證明：通過萬有引力進行普適定位是存在可行性的. 在理想狀況下，當總觀測時長達到20 min 的時候，定位精度約為1 km，在觀測總時長達到2 h 時，定位精度可以達到厘米級. 但是目前本文的算法對測量誤差極為敏感. 下一步的工作是對定位方法進行優(yōu)化和完善，提高其穩(wěn)健度，并結合高精度儀器進行實際測量，進行實際定位實驗，進一步對定位方法進行驗證. 如果證實可行，這種基于萬有引力的普適定位方法可以應用于地下、水下定位，由于完全被動無源，而且萬有引力難以被干擾，因此安全性高. 進一步的，這種定位方法還可以拓展應用到其他星球上的定位. 另外，由于誤差傳播放大的作用，有可能對引力傳播速度進行驗證.

致謝：在研究過程中，得到了吳立恒老師(應急管理部國家自然災害防治研究院)、邢樂林老師(中國科學院精密測量科學與技術創(chuàng)新研究院)、張燕老師(湖北省地震局)、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)賈劍鋼老師(武漢大學測繪學院)的幫助，他們提供的數(shù)據(jù)對于驗證算法很有幫助. 在此表示誠摯的感謝！