鄭慶敏

(海南省海洋地質(zhì)調(diào)查研究院,海南 ?？?570100)

在地球外部重力場的相關(guān)研究領(lǐng)域，如大地水準(zhǔn)面的精化、外部重力場的逼近等都需要處理大量的重力觀測資料以得到具有重要使用價(jià)值的重力場元素[1-3]。地球外部重力場的大量計(jì)算均需要重力點(diǎn)精確地心坐標(biāo)，即地心經(jīng)緯度和地心距，由于各種重力資料的來源、時(shí)期等的差異，致使已有重力測點(diǎn)的坐標(biāo)信息良莠不齊。根據(jù)相關(guān)知識(shí)，只有測點(diǎn)具備了的全部精確大地坐標(biāo)，即大地經(jīng)緯度和大地高的情況下才能轉(zhuǎn)換得到其精確地心坐標(biāo)[3-5]，而目前我國大部分物理大地測量試驗(yàn)區(qū)的大量重力點(diǎn)的坐標(biāo)是大地經(jīng)緯度和水準(zhǔn)高程，并不具備完整的三維空間坐標(biāo)，顯然，重力資料的現(xiàn)實(shí)情況和實(shí)際要求是矛盾的。

上述問題產(chǎn)生的原因是我國的重力網(wǎng)建設(shè)和三維坐標(biāo)控制網(wǎng)的建設(shè)不同步造成的。我國重力網(wǎng)的建設(shè)，從1957國家重力基本網(wǎng)(簡稱“57網(wǎng)”)到1985國家重力基本網(wǎng)(簡稱“85網(wǎng)”)，再到最新的“2000國家重力基本網(wǎng)”(簡稱“2000網(wǎng)”)，基準(zhǔn)點(diǎn)、基本點(diǎn)和引點(diǎn)的數(shù)量、點(diǎn)位分布合理度、覆蓋面和精度都有了很大的進(jìn)步[6-7]。相較之下，我國的三維坐標(biāo)基準(zhǔn)的建設(shè)則是始于上世紀(jì)90年代初，如國家GPS A級(jí)網(wǎng)于1992年布測，GPS B級(jí)網(wǎng)于1991～1995年布測，全國GPS一、二級(jí)網(wǎng)則是于1991～1997年布測，中國地殼運(yùn)動(dòng)觀測網(wǎng)絡(luò)則是從1998年開始布測。直到2008年7月1日我國正式啟用2000國家大地坐標(biāo)系，才標(biāo)志著實(shí)際測量工作進(jìn)入了具備三維空間坐標(biāo)的時(shí)期。然而，在此之前所完成的多數(shù)重力測量工作因?yàn)榉N種原因而與三維空間坐標(biāo)的關(guān)系沒有那么密切。在我國的重力控制網(wǎng)中，除了高等級(jí)的基準(zhǔn)點(diǎn)、基本點(diǎn)、引點(diǎn)及其他重要重力控制點(diǎn)具備完整的三維空間坐標(biāo)信息，甚至具備高等級(jí)的天文坐標(biāo)信息，其他大多數(shù)的低等級(jí)重力點(diǎn)的高程坐標(biāo)多是通過導(dǎo)線測量、水準(zhǔn)測量、三角高程導(dǎo)線等確定的，只有少數(shù)的重力點(diǎn)具備較可靠的GPS測量成果。再考慮到目前相對(duì)重力測量中使用的精度最高的CG-5相對(duì)重力儀，其測點(diǎn)的GPS坐標(biāo)也只是依靠一個(gè)粗定位的單天線而獲得的結(jié)果[8-9]，坐標(biāo)精度可想而知。

顯然，對(duì)我國重力網(wǎng)中坐標(biāo)信息不完備的大量低等級(jí)重力點(diǎn)進(jìn)行三維空間坐標(biāo)的補(bǔ)測是不現(xiàn)實(shí)的，因此需要尋找其它途徑來改善這一問題。本文針對(duì)具有大地經(jīng)緯度和水準(zhǔn)高信息的重力點(diǎn)，研究采用何種方法才能獲得比較精確的地心坐標(biāo)，其核心則是地心緯度的計(jì)算[10]。目的是在重力點(diǎn)位坐標(biāo)的確定過程中，將坐標(biāo)轉(zhuǎn)換導(dǎo)致的誤差控制在較低的程度，從而可靠地用于地球外部重力場的各種計(jì)算過程中。

1 具備大地高條件下的嚴(yán)密轉(zhuǎn)換過程

如圖1所示，在大地坐標(biāo)系O-XYZ中，已知P點(diǎn)大地坐標(biāo)(L,B,H)，要獲得其精密地心坐標(biāo)(相對(duì)于大地坐標(biāo)系參考橢球的質(zhì)心)，一般需要按照如下兩個(gè)步驟來實(shí)現(xiàn)[6,10]

圖1 大地坐標(biāo)與地心坐標(biāo)

1.1 根據(jù)大地坐標(biāo)(L,B,H)計(jì)算地心直角坐標(biāo)(X,Y,Z)如式(1)所示：

X=(N+H)cosBcosL

Y=(N+H)cosBsinL

Z=(N(1-e2)+H)sinB

(1)

1.2 由地心直角坐標(biāo)(X,Y,Z)計(jì)算地心球坐標(biāo)(λ,φ,r)如式(2)所示。

(2)

式(2)中的地心經(jīng)緯度的公式還可以表示為：

λ=L

(3)

(4)

顯然，經(jīng)度在轉(zhuǎn)換前后是不變的，而地心緯度φ則是大地緯度和大地高的函數(shù)，φ和B之間存在著明顯的差異，這里給出B-φ在B∈[0，90]區(qū)間的變化趨勢(取大地高為0)，如圖2和表1所示，此差值在大地緯度為45°時(shí)達(dá)到最大值11′32″.723 6，而在赤道和兩極時(shí)均為0。這也意味著在我國范圍內(nèi)進(jìn)行的重力場計(jì)算中，是不能使用大地緯度來代替地心緯度的。

圖2 B-φ在B[0, 90]區(qū)間的變化趨勢(H=0)

表1B-φ在B[0, 90]區(qū)間的變化趨勢(H=0)

B(°)B-φ ('.″)B(°)B-φ ('.″)00.000 0005011.225 97751.598 9415510.516 961103.561 8086010.009 259155.453 583658.518 059207.241 325707.264 236258.495 148755.473 731309.589 111803.576 7633510.502 007852.006 8984011.218 020900.000 0004511.327 236

圖3表示的是在B=45°處，B-φ在H[0, 3 000]m區(qū)間的變化趨勢，隨著大地高的增加，大地緯度與地心緯度之差呈現(xiàn)出線性衰減的特點(diǎn)，對(duì)應(yīng)的地心緯度從44°48′27″.2 764(H=0)增加至44°48′27″.6 026(H=3 000)。

圖3 B-φ在H[0, 3 000]m區(qū)間的變化趨勢(B=45°)

在涉及重力資料處理的背景下，大地緯度是已知量，而大地高則是含有一定的誤差的，進(jìn)一步得到地心緯度對(duì)大地高的偏導(dǎo)數(shù)如下式：

(5)

從式(5)可知：地心緯度、大地高、大地緯度之間的關(guān)系比較復(fù)雜，因此對(duì)缺乏精確大地高信息的重力資料而言，若要獲得重力資料的精密地心緯度是比較困難的。另外，由于無法獲得精確的地心坐標(biāo)，亦無法直接利用下式計(jì)算大地水準(zhǔn)面高Ng，

(6)

2 不具備大地高情況下的近似轉(zhuǎn)換過程

大地高的定義如圖所示，H=Ng+h。式中，Ng為大地水準(zhǔn)面高，h為水準(zhǔn)高程。亦可根據(jù)似大地水準(zhǔn)面來定義大地高，H=Hγ+ξ，式中，Hγ為正常高，ξ為高程異常。

為解決不具備大地高情況下獲取地心緯度的突出問題，本文提出如下的幾種近似方法來獲得重力資料的地心緯度。

圖4 大地高構(gòu)成圖

2.1 借助于已知大地水準(zhǔn)面數(shù)值模型的轉(zhuǎn)換

根據(jù)前述，已知重力測點(diǎn)的高程和大地經(jīng)緯度信息，并同時(shí)能獲得測點(diǎn)所在區(qū)域的(似)大地水準(zhǔn)面數(shù)值模型，為計(jì)算點(diǎn)位的地心經(jīng)緯度，可以利用數(shù)值內(nèi)插方法計(jì)算出測點(diǎn)的正常高或大地水準(zhǔn)面高，再進(jìn)一步計(jì)算出測點(diǎn)的大地高信息，然后利用第1節(jié)的嚴(yán)密轉(zhuǎn)換過程計(jì)算出測點(diǎn)的地心緯度。

2.2 一步法

參考圖1，由于P點(diǎn)的高程與地球半徑相比較為一小量，忽略P點(diǎn)的高程，即認(rèn)為P點(diǎn)位于參考橢球面上，直接利用公式

φ=arctan((1-f)2tanB)

(7)

將大地緯度轉(zhuǎn)換為地心緯度。其中,f為參考橢球的扁率。這種做法對(duì)于位于參考橢球面上的測點(diǎn)來說非常簡便，然而根據(jù)前述的條件，要確定測點(diǎn)是否位于參考橢球面上則是非常困難的，如果測點(diǎn)并不位于參考橢球面上，那么這種做法會(huì)帶來較大的誤差。

2.3 兩步法

根據(jù)圖1，P0點(diǎn)為P點(diǎn)沿參考橢球面法線在橢球面上的投影點(diǎn)，由于大地水準(zhǔn)面高相對(duì)于地球半徑來說是一小量，認(rèn)為P點(diǎn)的水準(zhǔn)高程即P點(diǎn)的大地高，

h≈H

首先借助于上述一步法的公式計(jì)算出P0點(diǎn)的地心緯度

φS=arctan((1-f)2tanB)

(8)

即將地球視為半徑為R的球體時(shí)對(duì)應(yīng)地面點(diǎn)的地心緯度，接著用下式計(jì)算出地心緯度

(9)

式中,h為海拔高(到平均海平面的距離)，f為參考橢球的扁率，R為地球的平均半徑，可取6 371 km。

顯然，與一步法相比較，兩步法的優(yōu)點(diǎn)就是顧及了測點(diǎn)的海拔高。

3 數(shù)值分析與討論

3.1 特征區(qū)域檢查點(diǎn)的計(jì)算

為了檢驗(yàn)上述三種坐標(biāo)轉(zhuǎn)換方法的精度，選取了中國范圍內(nèi)五個(gè)特征地區(qū)的已知重力點(diǎn)共計(jì)71點(diǎn)，分別位于西南、東南、中部、西北和東北等地。這些重力點(diǎn)的坐標(biāo)和其它成果齊全，不僅具有CGCS2000系統(tǒng)的大地坐標(biāo)，而且具有1985國家高程基準(zhǔn)的水準(zhǔn)高程信息，其重力系統(tǒng)為2000國家重力基本網(wǎng)系統(tǒng)。其統(tǒng)計(jì)信息如表2、表3所示。

表2 已知重力點(diǎn)數(shù)量

表3 區(qū)域內(nèi)重力點(diǎn)的高程特點(diǎn)

在分析過程中，首先利用上述測點(diǎn)的CGCS2000大地坐標(biāo)計(jì)算出其地心緯度的真實(shí)值，然后在假設(shè)缺乏點(diǎn)位大地高信息的情況下，分別利用上述的兩種方法計(jì)算出測點(diǎn)的地心緯度，分別與真實(shí)值作差比較，并進(jìn)一步將比較結(jié)果轉(zhuǎn)化為對(duì)應(yīng)的地面距離，結(jié)果如圖5、圖6和表4所示。

圖5 模型法和兩步法結(jié)果與真實(shí)值之差

圖6 三種方法結(jié)果與真實(shí)值之差

表4三種方法引起的坐標(biāo)誤差統(tǒng)計(jì)/m

方法最大值最小值均值中誤差模型法0.029 1-0.133 2-0.053 30.062 0兩步法0.029 7-0.132 6-0.053 40.062 2一步法9.351 90.059 73.138 83.213 8

綜合分析圖5、圖6及表4可知：在不具備測點(diǎn)大地高信息的情況下，利用大地水準(zhǔn)面模型法轉(zhuǎn)換得到的地心緯度與真實(shí)值的差別最小、其次是兩步法的結(jié)果、一步法的轉(zhuǎn)換結(jié)果與真實(shí)值的差別最大；兩步法與模型法在轉(zhuǎn)換得到地心緯度的精度上非常接近，因此建議在無法獲得可用的大地水準(zhǔn)面數(shù)值模型的情況下，可用兩步法或模型法求取測點(diǎn)的地心緯度。

3.2 重力點(diǎn)的重力值、高程和坐標(biāo)的精度要求分析

測定重力點(diǎn)的目的在于得到測點(diǎn)的空間異常值?？臻g異常Δg按下式計(jì)算

Δg=g-γ0+0.308 6h

(10)

式(10)中，g為地面上的重力值，h為該點(diǎn)的高程，γ0為正常橢球面上相應(yīng)點(diǎn)的正常重力值，它可由該點(diǎn)的緯度由正常重力公式算得：

γ0=γα(1+βsin2φ-β1sin22φ)

(11)

由式(10)可以導(dǎo)出空間異常的中誤差為：

(12)

由式(11)得：

mγ0=γαβsin2φmφ

(13)

以γα=978 030 mGal和β=0.005 302代入，取sin2φ=1，得：

mγ0=0.025mφ

(14)

式中，mφ為重力點(diǎn)緯度的誤差，單位S，算得的mγ0的單位為mGal。將式(14)代入式(12)得到

(15)

式(15)給出了重力點(diǎn)空間異常中誤差與觀測重力值、點(diǎn)位和高程這三個(gè)量的測定誤差之間的關(guān)系。根據(jù)對(duì)空間異常的精度要求，按式(15)可對(duì)上面三個(gè)量的中誤差做出一定的安排?？紤]到點(diǎn)位誤差不能過大，因?yàn)辄c(diǎn)位誤差大了必然影響異常值的質(zhì)量，并顧及到測定重力值、點(diǎn)位和高程各自的困難程度，如在要求空間異常中誤差為0.1 mGal時(shí)，可按表5進(jìn)行合理的安排限差。

表5 空間異常誤差對(duì)坐標(biāo)限差的要求

結(jié)合前一小節(jié)及這里得到的限差要求，可知即使在目前對(duì)空間異常誤差要求嚴(yán)格的情況下，前面三種坐標(biāo)轉(zhuǎn)換所致誤差仍在表5中坐標(biāo)限差的范圍內(nèi)。這就意味著，即使我們使用轉(zhuǎn)換誤差較大的一步法獲得測點(diǎn)的地心緯度，仍然不會(huì)對(duì)最終的空間異常造成較大的影響。然而，隨著未來測量手段的進(jìn)步，空間異常的精度將進(jìn)一步提高，對(duì)坐標(biāo)誤差的限差要求也會(huì)更加嚴(yán)格，在此情況下，精度較高而又便利的兩步法將是更好的選擇。另外，綜合表4與表5的結(jié)果可知，大地水準(zhǔn)面模型法和兩步法的轉(zhuǎn)換誤差在坐標(biāo)誤差限差中所占比例非常小，而一步法的轉(zhuǎn)換誤差在坐標(biāo)誤差限差中所占比例稍大。

4 結(jié) 語

本文提出利用已有信息獲得地心緯度的幾種方法，并通過數(shù)值試驗(yàn)對(duì)轉(zhuǎn)換方法進(jìn)行了分析比較，綜合分析結(jié)果可以得到如下結(jié)論：

(1)在大地水準(zhǔn)面模型法、兩步法和一步法這三種方法中，大地水準(zhǔn)面模型法的轉(zhuǎn)換精度最高，其前提是需要轉(zhuǎn)換區(qū)域的大地水準(zhǔn)面數(shù)值模型。兩步法的轉(zhuǎn)換精度稍低，但與大地水準(zhǔn)面模型法的轉(zhuǎn)換精度非常接近，并且此方法非常直接便利。一步法的精度最低。

(2)在某些對(duì)空間異常精度要求嚴(yán)格的情況下，三種方法的轉(zhuǎn)換誤差都在坐標(biāo)誤差限差的范圍內(nèi)。其中大地水準(zhǔn)面模型法和兩步法的轉(zhuǎn)換誤差在坐標(biāo)誤差限差中所占比例非常小，而一步法的轉(zhuǎn)換誤差在坐標(biāo)誤差限差中所占比例稍大。

(3)在重力資料的處理中若遇到坐標(biāo)轉(zhuǎn)換問題，我們推薦使用兩步法。另外，在計(jì)算出地心緯度后，可進(jìn)一步利用式(6)算出大地水準(zhǔn)面高，可得到更加精確的大地高，以及更加精確的地心距，從而使基于地心坐標(biāo)的重力場計(jì)算更加便利。