重力資料處理中的坐標(biāo)轉(zhuǎn)換問(wèn)題研究

鄭慶敏

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

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

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

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

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

如圖1所示，在大地坐標(biāo)系O-XYZ中，已知P點(diǎn)大地坐標(biāo)(L,B,H)，要獲得其精密地心坐標(biāo)(相對(duì)于大地坐標(biāo)系參考橢球的質(zhì)心)，一般需要按照如下兩個(gè)步驟來(lái)實(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ū)間的變化趨勢(shì)(取大地高為0)，如圖2和表1所示，此差值在大地緯度為45°時(shí)達(dá)到最大值11′32″.723 6，而在赤道和兩極時(shí)均為0。這也意味著在我國(guó)范圍內(nèi)進(jìn)行的重力場(chǎng)計(jì)算中，是不能使用大地緯度來(lái)代替地心緯度的。

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

表1B-φ在B[0, 90]區(qū)間的變化趨勢(shì)(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ū)間的變化趨勢(shì)，隨著大地高的增加，大地緯度與地心緯度之差呈現(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ū)間的變化趨勢(shì)(B=45°)

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

(5)

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

(6)

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

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

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

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

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

根據(jù)前述，已知重力測(cè)點(diǎn)的高程和大地經(jīng)緯度信息，并同時(shí)能獲得測(cè)點(diǎn)所在區(qū)域的(似)大地水準(zhǔn)面數(shù)值模型，為計(jì)算點(diǎn)位的地心經(jīng)緯度，可以利用數(shù)值內(nèi)插方法計(jì)算出測(cè)點(diǎn)的正常高或大地水準(zhǔn)面高，再進(jìn)一步計(jì)算出測(cè)點(diǎn)的大地高信息，然后利用第1節(jié)的嚴(yán)密轉(zhuǎn)換過(guò)程計(jì)算出測(cè)點(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ì)于位于參考橢球面上的測(cè)點(diǎn)來(lái)說(shuō)非常簡(jiǎn)便，然而根據(jù)前述的條件，要確定測(cè)點(diǎn)是否位于參考橢球面上則是非常困難的，如果測(cè)點(diǎn)并不位于參考橢球面上，那么這種做法會(huì)帶來(lái)較大的誤差。

2.3 兩步法

根據(jù)圖1，P0點(diǎn)為P點(diǎn)沿參考橢球面法線在橢球面上的投影點(diǎn)，由于大地水準(zhǔn)面高相對(duì)于地球半徑來(lái)說(shuō)是一小量，認(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)就是顧及了測(cè)點(diǎn)的海拔高。

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

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

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

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

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

在分析過(guò)程中，首先利用上述測(cè)點(diǎn)的CGCS2000大地坐標(biāo)計(jì)算出其地心緯度的真實(shí)值，然后在假設(shè)缺乏點(diǎn)位大地高信息的情況下，分別利用上述的兩種方法計(jì)算出測(cè)點(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可知：在不具備測(cè)點(diǎn)大地高信息的情況下，利用大地水準(zhǔn)面模型法轉(zhuǎn)換得到的地心緯度與真實(shí)值的差別最小、其次是兩步法的結(jié)果、一步法的轉(zhuǎn)換結(jié)果與真實(shí)值的差別最大；兩步法與模型法在轉(zhuǎn)換得到地心緯度的精度上非常接近，因此建議在無(wú)法獲得可用的大地水準(zhǔn)面數(shù)值模型的情況下，可用兩步法或模型法求取測(cè)點(diǎn)的地心緯度。

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

測(cè)定重力點(diǎn)的目的在于得到測(cè)點(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)空間異常中誤差與觀測(cè)重力值、點(diǎn)位和高程這三個(gè)量的測(cè)定誤差之間的關(guān)系。根據(jù)對(duì)空間異常的精度要求，按式(15)可對(duì)上面三個(gè)量的中誤差做出一定的安排?？紤]到點(diǎn)位誤差不能過(guò)大，因?yàn)辄c(diǎn)位誤差大了必然影響異常值的質(zhì)量，并顧及到測(cè)定重力值、點(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)換誤差較大的一步法獲得測(cè)點(diǎn)的地心緯度，仍然不會(huì)對(duì)最終的空間異常造成較大的影響。然而，隨著未來(lái)測(cè)量手段的進(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é) 語(yǔ)

本文提出利用已有信息獲得地心緯度的幾種方法，并通過(guò)數(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)換問(wèn)題，我們推薦使用兩步法。另外，在計(jì)算出地心緯度后，可進(jìn)一步利用式(6)算出大地水準(zhǔn)面高，可得到更加精確的大地高，以及更加精確的地心距，從而使基于地心坐標(biāo)的重力場(chǎng)計(jì)算更加便利。