于錦海，徐 煥，萬曉云

1. 中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100049; 2. 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京 100049; 3. 中國地質(zhì)大學(xué)(北京)土地科學(xué)與技術(shù)學(xué)院，北京 100083

地球重力場(或引力場)的研究工作目前已經(jīng)取得了較大進(jìn)展。例如：綜合利用各類重力觀測數(shù)據(jù)建立了EGM2008模型，重力場的分辨率有了極大的提高。又如：衛(wèi)星重力技術(shù)的應(yīng)用，特別是GRACE衛(wèi)星計(jì)劃的實(shí)施，應(yīng)用時(shí)變重力場信息可以解釋地球質(zhì)量的遷移特征[1-2]。此外，重力測量技術(shù)也獲得了較大的發(fā)展，特別是衛(wèi)星測高技術(shù)與重力梯度測量技術(shù)的應(yīng)用[3-6]，為解算重力場提供了豐富的觀測數(shù)據(jù)。所有這些關(guān)于重力場方面的進(jìn)展，使得進(jìn)一步地拓展重力場應(yīng)用成為可能。

重力輔助導(dǎo)航是重力場理論的潛在應(yīng)用方向之一[7-9]。自20世紀(jì)90年代開始，便有國外的學(xué)者進(jìn)行了研究與模擬測試[10-14]，而我國也有學(xué)者開展了相應(yīng)的研究工作[15-16]。若將上述研究工作進(jìn)一步延伸，能否利用重力理論來實(shí)現(xiàn)自主定位與導(dǎo)航(即單純使用重力測量數(shù)據(jù)實(shí)現(xiàn)目標(biāo)的定位與導(dǎo)航)。顯然，這是一項(xiàng)極有挑戰(zhàn)的工作，其中包含了利用重力進(jìn)行定位與導(dǎo)航的理論問題、重力測量儀器的研制、高精度與高分辨率的重力場模型研發(fā)等諸多因素[17]。

如何建立重力場(引力場)數(shù)據(jù)與相應(yīng)空間點(diǎn)位之間的直接關(guān)系，以及建立的關(guān)系有何應(yīng)用價(jià)值，是能否利用重力數(shù)據(jù)進(jìn)行自主定位與導(dǎo)航的基礎(chǔ)理論問題。本文的目標(biāo)從引力場對應(yīng)的不變量入手，引入一組曲線坐標(biāo)系，并研究其性質(zhì)和探討潛在的應(yīng)用。因?yàn)橐氲那€坐標(biāo)系是由不變量構(gòu)成的，因此該坐標(biāo)系直接描述引力場與空間點(diǎn)位之間直接的關(guān)系。

1 引力場的不變量

(1)

(2)

經(jīng)計(jì)算，這些不變量可寫為[19-20]

(3)

由于引力位v滿足Laplace方程，即I1=0，所以從引力梯度數(shù)據(jù)可構(gòu)造出函數(shù)I2和I3。因?yàn)镮2和I3是不變量，這意味著I2和I3也僅是點(diǎn)位的函數(shù)，與引力梯度分量的方向無關(guān)。

至此，就從引力場自身引入了3個(gè)僅與點(diǎn)位相關(guān)的函數(shù)g、I2和I3。由于這3個(gè)函數(shù)都有各自的物理量綱，所以需進(jìn)行去量綱化處理。

為了確保不變量坐標(biāo)系(f1,f2,f3)構(gòu)成了曲線坐標(biāo)系，需要在理論上證明相應(yīng)的Jacobi行列式是非零的，即

(4)

事實(shí)上，因?yàn)镴acobi行列式J是不恒為零的，所以從三維空間分布看面，J=0僅可能構(gòu)成空間的曲面，圖1的結(jié)果表明，該曲面大致是赤道沿徑向向外延伸的曲面。

圖1 利用EGM2008前180階次模型計(jì)算的Jacobi行列式在地球平均球面上的分布Fig.1 The distribution of Jacobi’s determinant for EGM2008 model with the former 180 degree and order on the average sphere of the earth

2 不變量坐標(biāo)系與空間點(diǎn)位之間的關(guān)系

本節(jié)將通過某些實(shí)例計(jì)算來驗(yàn)證：假若事先給出了不變量坐標(biāo)系的具體表達(dá)式，則利用該表達(dá)式就可以進(jìn)行空間點(diǎn)的定位。正如地圖一樣，通過對若干參照物的觀測可以利用地圖作為工具來確定點(diǎn)位的坐標(biāo)，因此在某種意義上看，不變量坐標(biāo)系起著地圖的作用。

本節(jié)將以EGM2008引力場模型的前360階次作為實(shí)際地球引力場，由此便可以得到相應(yīng)的引力不變量曲線坐標(biāo)系(f1,f2,f3)。假設(shè)在地面上某點(diǎn)P進(jìn)行引力和引力梯度觀測，并進(jìn)行數(shù)據(jù)處理后得到了P點(diǎn)在不變量坐標(biāo)系下的值(F1,F2,F3)，其目標(biāo)是解算出P點(diǎn)的空間坐標(biāo)(r,θ,λ)。事實(shí)上，從P點(diǎn)的觀測值(F1,F2,F3)，可以建立下列方程組

(5)

由于不變量坐標(biāo)系(f1,f2,f3)是已知的，故通過解算上述方程組便能得到P點(diǎn)對應(yīng)的球坐標(biāo)。顯然方程組(5)關(guān)于待解算的變量(r,θ,λ)是非線性的，因此直接進(jìn)行求解是難以實(shí)現(xiàn)的。理論上講，對于形式如方程組(5)的非線性方程而言，只要在某點(diǎn)處的Jacobi行列式非零，那么在該點(diǎn)附近是局部可解的，用數(shù)學(xué)語言描述，就是存在該點(diǎn)的某個(gè)鄰域，在該鄰域內(nèi)方程組是唯一可解，而且解還具有穩(wěn)定性[21]。至于求解方法通常是采用Newton迭代法，即線性化迭代解法。

下面就方程組(5)，給出具體的線性化迭代過程如下

(6)

這里j=0,1,2,…，而初始值(r0,θ0,λ0)的選擇與待計(jì)算的P點(diǎn)要盡可能地接近，其理由就是因?yàn)榉蔷€性方程組的解僅在局部具有存在性、唯一性、穩(wěn)定性。在方程組(6)中，待解的量是(rj+1,θj+1,λj+1)，因此方程組(6)關(guān)于解算量是線性的，所以方程組(6)稱為方程組(5)的線性化迭代形式。由于方程組(6)是線性的，所以其可解性完全取決于在(rj,θj,λj)處的Jocobi行列式J(rj,θj,λj)。如果能通過實(shí)際解算驗(yàn)證方程組(6)關(guān)于(rj+1,θj+1,λj+1)是唯一可解的，那么便間接地證明了Jocobi行列式是非零的。此外，至于需要進(jìn)行多少次迭代，這取決于事先給出的誤差。事實(shí)上，只要迭代前后解的差小于誤差要求，便可以停止迭代，從而得到所需的解。

為了驗(yàn)證方程組(5)以及線性化迭代形式(6)的可解性，本文通過3個(gè)低、中、高緯度實(shí)際點(diǎn)位上的不變量坐標(biāo)值對其進(jìn)行了驗(yàn)證，其解算結(jié)果匯總見表1—表3。

表1 低緯度點(diǎn)驗(yàn)證結(jié)果

表2 中緯度點(diǎn)驗(yàn)證結(jié)果

表3 高緯度點(diǎn)驗(yàn)證結(jié)果

表1—表3中，第1列是點(diǎn)位坐標(biāo)與實(shí)測不變量坐標(biāo)值的記號，第2列是實(shí)際的點(diǎn)位與不變量坐標(biāo)值，第3列是迭代形式的初始值(r0,θ0,λ0)以及對應(yīng)的不變量坐標(biāo)值，第4列是通過迭代解算式(6)后得到結(jié)果，第5與第6列則是初始結(jié)果和最終計(jì)算結(jié)果分別與實(shí)際值的差，而最下面的行(距離)則是誤差。表1表明，對于低緯度點(diǎn)，迭代初始值的點(diǎn)位距實(shí)際解算點(diǎn)位17.676 km，通過迭代計(jì)算后得到的最終點(diǎn)位坐標(biāo)與實(shí)際點(diǎn)位的差僅有3.87 mm。表2表明，對于中緯度點(diǎn)，迭代初始值的點(diǎn)位距實(shí)際解算點(diǎn)位32.634 km，通過迭代計(jì)算后得到的最終點(diǎn)位坐標(biāo)與實(shí)際點(diǎn)位的差僅有1.34 mm。表3表明，對于高緯度點(diǎn)，迭代初始值的點(diǎn)位距實(shí)際解算點(diǎn)位42.519 km，通過迭代計(jì)算后得到的最終點(diǎn)位坐標(biāo)與實(shí)際點(diǎn)位的差僅有7.69 mm。為了清晰地反映迭代求解過程，本文將表1—表3中的計(jì)算迭代過程與次數(shù)匯制成圖2。從圖2能清晰可見關(guān)于低、中、高緯度點(diǎn)位的迭代解算過程，此外迭代次數(shù)分別是4、3、9次。

圖2 迭代求解過程Fig.2 Graphical processes of iteration

通過上述關(guān)于低、中、高緯度點(diǎn)的驗(yàn)證，可得到如下結(jié)論，若事先能給出不變量坐標(biāo)系，那么通過重力與重力梯度觀測是可以反解出空間點(diǎn)的位置坐標(biāo)的，這為單純利用重力數(shù)據(jù)進(jìn)行自主定位導(dǎo)航提供了新的途徑。

需要說明的是：①若出現(xiàn)Jacobi行列式為零的情況，那么本文給出的方法將無法反演出點(diǎn)位坐標(biāo)，如圖1所示，在赤道附近有可能會出現(xiàn)這樣的情況；②諸如引力(重力)和梯度測量是可以連續(xù)實(shí)施的，因此在求解線性化方程組(6)時(shí)其初始值通常應(yīng)該取為上一個(gè)點(diǎn)位的坐標(biāo)，這樣就能保證初始值與定位點(diǎn)的距離差距不大。

3 結(jié)論與潛在的應(yīng)用

本文利用地球引力場(或重力場)的性質(zhì)，描述了相應(yīng)的引力不變量曲線坐標(biāo)系。利用該坐標(biāo)系，可以進(jìn)行空間點(diǎn)位的反演計(jì)算。事實(shí)上，不變量坐標(biāo)系在定位時(shí)就像地圖一樣，是需要事先確定的。不同的是，地圖是紙質(zhì)或數(shù)字形式表現(xiàn)出來的，而不變量坐標(biāo)系是以函數(shù)形式表現(xiàn)出來的。若要將本文論述不變量曲線坐標(biāo)系用于實(shí)際定位與導(dǎo)航應(yīng)用，大致需要進(jìn)行下列幾項(xiàng)工作：①通過各種引力(重力)測量方法來獲取相關(guān)數(shù)據(jù)，對其進(jìn)行處理后給出不變量坐標(biāo)系的表達(dá)形式；②在進(jìn)行實(shí)際定位或?qū)Ш綍r(shí)，需要進(jìn)行實(shí)時(shí)的引力與引力梯度測量，即需要精度滿足要求的加速度計(jì)和梯度測量儀器；③采用理論上可行的算法，利用事先建立的不變量坐標(biāo)系與實(shí)時(shí)測量值來反演點(diǎn)位的坐標(biāo)。本文工作主要是論證上述方法在理論上的可行性，并針對第3部分內(nèi)容進(jìn)行分析與討論。

如何事先給出引力(重力)不變量曲線坐標(biāo)系，這取決于地球引力場(重力場)模型的建立。如果建立了精度足夠高的引力場(重力場)模型，則利用該模型自然就可以生成不變量坐標(biāo)系。目前EGM2008重力場模型已經(jīng)達(dá)到了2160階次，分辨率可達(dá)8 km左右，而國內(nèi)和國際上即將推出更高階次的重力場模型，這對于建立高精度的全球不變量曲線坐標(biāo)系無疑有著直接的作用。對于局部問題(例如南海區(qū)域)可以通過在該區(qū)域進(jìn)行高分辨率的重力測量，結(jié)合衛(wèi)星海洋測高數(shù)據(jù)解算出區(qū)域重力場模型，例如，可采用球冠諧級數(shù)模型或移去恢復(fù)法建立局部模型[22-26]，由此便可得到局部區(qū)域相應(yīng)的不變量坐標(biāo)系?？傊?，要實(shí)現(xiàn)事先能給出不變量曲線坐標(biāo)系，需要進(jìn)行大量的實(shí)際測量、數(shù)據(jù)處理等系列工作，正如要繪制出精確的地圖一樣，需要事先完成系列工程性的工作。顯然，這不是本文能夠解決的問題。

為何引入引力梯度不變量作為曲線坐標(biāo)系，而不是采用梯度的分量。這是因?yàn)樘荻确至坎粌H與點(diǎn)位有關(guān)，而且依賴于方向，這會導(dǎo)致更多參數(shù)的引入。單純從數(shù)學(xué)角度講，若使用梯度分量，則對應(yīng)數(shù)組已不再是曲線坐標(biāo)系的概念了。從實(shí)際應(yīng)用看，更多參數(shù)的引入會增加理論與計(jì)算的復(fù)雜性，甚至?xí)?dǎo)致欠定問題，即給出的條件數(shù)少于待解未知數(shù)的個(gè)數(shù)。

在引力不變量曲線坐標(biāo)系中，引力(重力)是可以直接測量的，而引力(重力)梯度不變量I2和I3則是梯度分量的組合。事實(shí)上，文獻(xiàn)[20]曾證明：計(jì)算I2和I3僅需梯度的3個(gè)對角分量vrr、vθθ和vλλ即可，理由是這3個(gè)對角分量是梯度分量的主項(xiàng)，其余分量的計(jì)算可利用重力場模型迭代計(jì)算來解決。對于該3個(gè)對角分量，由于徑向r與引力(重力)的反方向相差不大，所以可以利用垂線方向作為標(biāo)定軸向，其余兩個(gè)軸向只要與標(biāo)定軸垂直即可?？傊?，在計(jì)算不變量I2和I3時(shí)，梯度的3個(gè)對角分量是必不可少的。概括來說，若在某點(diǎn)給出了引力以及3個(gè)梯度對角分量(共4個(gè)觀測量)，則可以算出該點(diǎn)處的不變量坐標(biāo)值，從而解出該點(diǎn)的空間坐標(biāo)位置。如果不采用不變量坐標(biāo)系，需要解算的量是坐標(biāo)位置(3個(gè)分量)和梯度姿態(tài)(3個(gè)歐拉角)共6個(gè)未知數(shù)，因此將導(dǎo)致欠定問題，即4個(gè)觀測量，6個(gè)未知數(shù)。

至于第2項(xiàng)主要工作，即加速度計(jì)與梯度測量儀器的研制與應(yīng)用，顯然是能否將本文提出的方法予以實(shí)施的關(guān)鍵。對于引力(重力)測量，微伽級精度的加速度計(jì)已經(jīng)用于實(shí)際測量，因此引力測量的條件是具備的。關(guān)于引力(重力)梯度測量，目前得到應(yīng)用的最好成果當(dāng)屬GOCE衛(wèi)星上搭載的梯度測量儀[27]，其精度能達(dá)到幾個(gè)mE(10-12s-2)量級[3]。此外，美國Maryland大學(xué)[28]在實(shí)驗(yàn)室里研制出了精度達(dá)到0.14 mE的重力梯度儀，而我國也有研究機(jī)構(gòu)正致力于重力梯度儀的研制。期待不久的將來，能生產(chǎn)出可用于實(shí)際測量的重力梯度儀。

本文的目的是闡述單純利用引力場(或重力場)理論進(jìn)行自主定位與導(dǎo)航的可行性。盡管學(xué)術(shù)界已經(jīng)提出過重力輔助匹配導(dǎo)航的思路，但是以曲線坐標(biāo)系的形式進(jìn)行引力場(重力場)獨(dú)立定位導(dǎo)航的思路尚屬首次提出，這對于拓展重力場的應(yīng)用無疑是有益的。