利用大地邊值問題方法統一全球高程基準

趙德軍，徐新強，歐陽明達

(1. 信息工程大學地理空間信息學院，河南 鄭州 450001； 2. 地理信息工程國家重點實驗室，陜西 西安 710054； 3. 西安測繪總站，陜西 西安 710054)

眾所周知，當地平均海水面并不重合于全球大地水準面，它們的差異稱為海面地形。全球海面地形最大差異可達2 m。例如，我國東海地區各長期驗潮站的平均海面比青島驗潮站代表的國家高程基準要高出25 cm左右，而南海驗潮站確定的平均海面要高出1985國家高程基準35 cm以上。因此建立在當地平均海水面上的各國高程基準就不一樣。目前全世界存在超過100個高程基準，如何統一它們成為大地測量、海洋學界努力的目標。

通常有3種方法統一高程基準：①精密水準測量法。最著名的工程就是歐洲垂直參考框架EVRF2007，將27個國家不同高程系統、不同潮汐系統的高程基準統一為正常高系統、零潮汐系統[1]；以及北美大陸垂直基準NAVD88，通過美國、加拿大和墨西哥水準網重新平差而來。②海面地形法。隨著衛星測高技術的發展，已能以3～5 cm的精度確定平均海水面，因此可嘗試通過海面地形的方法來統一高程基準。文獻[2]通過對北美太平洋、大西洋沿岸、歐洲大西洋沿岸的130多個具有GNSS大地高的驗潮站資料分析發現，平均海面高模型結合GOCE+EGM2008重力大地水準面能以分米級的精度傳遞高程基準。③重力位差法。該方法是確定局部高程零點的重力位，又分為兩類：一種是正高或正常高反算法，即根據正高或正常高的定義反算高程基準零點的重力位，文獻[3—6]利用該方法計算了我國1985國家高程基準零點的重力位，文獻[7]用該方法實現了沿海島礁高程傳遞，文獻[8]用該方法確定了希臘16個島嶼水準原點重力位，但這種方法需要用全球重力位系數模型來計算水準點的重力位，位系數模型的精度是影響該方法的一個重要因素；另一種方法就是基于大地邊值問題方法，該方法利用含有局部高程基準信息的“有偏”重力異常，按邊值問題理論確定重力大地水準面，從而間接解算出局部高程基準位差。

基于邊值問題法最早由Colombo提出，后經Rummel等進一步拓展，發展為利用改化的Stokes積分公式計算不同基準之間的位差[9-10]。文獻[11]應用該方法計算出的瑞典高程與芬蘭高程系統的系統差異與水準測量結果一致。文獻[12—13]將該方法引入中國，并拓展了正常高系統中高程零點重力位的計算方法[12]。但是依賴于局部高程基準信息的“有偏”重力異常會干擾重力位差的解算，引出“間接偏差項”的概念，這也是該方法沒有流行起來的原因。有兩種方法處理間接偏差項，一種是新西蘭采用的迭代法，首先用“有偏”重力異常確定出一個粗略的高程基準位差，然后用這個粗略的高程基準位差去改正“有偏”重力異常，再次計算較精確的高程基準位差，如此反復迭代[14]。我國海島(礁)測繪工程中也采用新西蘭的迭代法來確定大陸高程基準重力位。另一種方法是利用GOCE衛星重力場模型計算中長波的全球大地水準面，以此來減弱間接偏差項的影響。2010年歐空局推出了“GOCE+”項目，其第1個主題就是探索GOCE衛星重力場統一全球高程基準的能力，由于純衛星重力大地水準面不受任何局部高程基準影響，可以作為全球高程基準。特別是文獻[15—16]從理論上證明了利用200階次的GOCE衛星重力場，間接偏差項將小于1 cm，從而可忽略其影響。

文獻[9—16]中確定局部高程基準重力位差的方法，無論是Stokes解，還是Molodensky解，均是基于第三邊值問題。后來我國學者提出采用線性固定重力邊值問題法來確定高程基準位差，并成功應用于深圳和香港的高程基準統一[17]。其基本原理為：與重力異常相比，GNSS大地高和重力確定的擾動引力獨立于任何高程基準，是“無偏”擾動引力，其線性化球近似的解其實就是第二邊值問題解，可得到獨立于任何高程基準的重力水準面，用此獨立的重力水準面作為標尺來統一全球高程基準。由于GNSS技術出現之前的重力數據均以重力異常表示，使得目前基于第二邊值問題的方法還難以推廣。因此，本文還是以基于第三邊值問題來討論如何聯接各局部高程基準。

1 大地邊值問題法

如上文所述，無論基于第二、還是第三大地邊值問題，統一高程基準的關鍵是要確定獨立于任何高程基準的全球統一的大地水準面或似大地水準面。

1.1 (似)大地水準面之間的平行性

大地水準面是重力等位面，其上的重力位處處相等。任何一個局部高程基準l與全球高程基準W0之間的位差Cl0也處處相等，它們之間的差距為

(1)

大地水準面具有平行性，但似大地水準面不是重力等位面，似大地水準面之間是否平行文獻[16,18]給出了證明：似大地水準面雖不是重力等位面，但重力位差Cl0卻處處相等，因此似大地水準面之間的差距為

(2)

1.2 純地面重力的邊值問題解

本文以Stokes理論為例來解第三邊值問題，Molodensky理論參見文獻[12,18]。根據廣義布隆斯公式，局部高程基準的大地水準面Nl與全球大地水準面上擾動位T0的關系為[9-11]

(3)

式中，ΔW0=W0-U0，為大地水準面重力位W0與參考橢球正常重力位U0之差；γ為參考橢球面上的正常重力。在Stokes理論中，要確定全球大地水準面需要全球大地水準面上的重力異常，但是只能得到局部高程基準l中的重力異常，全球大地水準面上和局部大地水準面重力異常關系為[9-11]

(4)

式中，Δg為從地面歸算到全球大地水準面上的重力異常；Δgl為從地面歸算到局部大地水準面l上的重力異常，因其嚴重依賴于局部高程基準，也稱“有偏”重力異常；R為地球平均半徑。式(4)表明，全球大地水準面上的重力異常Δg，可用局部高程基準上的重力異常Δgl來表示。第三邊值問題的廣義Stokes解為

(5)

式中，δGM=GM-GM0，為實際地球質量引力常數GM和參考橢球質量引力常數GM0之差；S(ψ)為Stokes函數。進一步整理得到

(6)

若存在K個高程基準，每個高程基準中重力位差Ck0只有一個，則式(6)右端第3項可以用和的形式表示，即

(7)

將式(6)和式(7)代入式(3)，得

(8)

令

(9)

式中，Nzero為大地水準面0階項，目前能精確確定，可作為已知量。令

(10)

式中，NStokes為地面重力異常計算的大地水準面，可作為已知量。令

(11)

式中，Nind稱為間接偏差項，其中

(12)

則式(8)形式上可簡化為

(13)

式中，Nl=h-Hl，h為GNSS獲取的大地高，Hl為局部高程基準中的正高，Nl為已知量。

式(14)中未知量存在直接位差Cl0，以及其他高程基準中的位差Ck0。可按式(13)建立觀測方程，平差解算出未知量，文獻[9]建立了平差模型，文獻[19]根據已有觀測資料將全世界分成亞洲、歐洲、非洲等7個高程基準進行模擬計算。由于Stokes函數收斂速度慢，因此其積分需要跨越多個高程基準，使得觀測方程變得復雜，而不實用。

1.3 聯合衛星和地面重力的邊值問題解

(14)

(15)

式中，Pn(cosψ)為勒讓德級數。改化Stokes核函數的作用為加速積分收斂，減小球冠Ω0積分面積，同時間接偏差項變為

(16)

對比式(16)和式(12)，區別僅僅在于改化Stokes核函數，但是卻起到了驚人的效果。文獻[15]證明使用50階次的衛星重力，即改化Stokes核函數nmax=50時，間接偏差項ISk影響小于3 cm；而使用200階次的衛星重力(nmax=200)，間接偏差項ISk影響小于1 cm，因而可直接忽略間接偏差項的影響。聯合衛星重力場和地面重力的全球大地水準面為

(17)

最后給出高程基準重力位差的計算式為

Cl0=γ(h-Hl-N0)

(18)

正常高系統中，重力位差的計算式為

(19)

2 數值試驗

2.1 GNSS/水準數據

GNSS/水準數據采用美國大地測量局NGS的GNSS/水準數據集——GPSBM2009，該數據集參與了美國混合大地水準面GEOID12A的構建。GPSBM2009數據的大地基準采用北美大地基準NAD83，高程基準采用北美高程基準NAVD88。NAD83坐標系原點與現代ITRF地心存在約2 m的偏差，定向也存在一定的系統差，因此數據預處理時采用NGS提供的Htdp坐標轉換軟件，將大地坐標轉換到ITRF2008框架[22]。NAVD88屬于單驗潮站確定的高程基準，通過重新平差美國、加拿大、墨西哥歷史上幾十年的水準測量數據獲得，并非開展北美全洲統一的水準測量，存在一定的系統差，尤其在美加邊境。因此加拿大最終放棄了NAVD88高程基準，仍采用CGVD28作為法定高程基準，只有美國、墨西哥采用NAVD88作為法定高程基準。

2.2 重力大地水準面

從加拿大自然資源部收集了加拿大重力水準面CGG2010，其采用200階次的GOCE重力場模型goco01s，覆蓋范圍包括加拿大和美國[23]。美國大地測量局收集了美國重力水準面USGG2012，其采用200階的GOCE模型goco02s，覆蓋范圍包括美國和加拿大大部分國土。250階次的GOCE和251～2160階次的EGM2008組合的重力大地水準面GOCE+EGM08，GOCE有多種版本，經比較第5代直接法解算的模型精度較高[22]。

2.3 重力位差

NAVD88高程基準在美國境內與GOCE重力水準面存在明顯的東西方向和南北方向的傾斜(如圖1 所示)，因此采用2參數曲面模型來吸收系統誤差，即

δN=x1(φp-φ0)+x2(λp-λ0)cosφp

(20)

式中，δN=hNAD83-HNAVD88-NGG為大地水準面差異，hNAD83為GNSS大地高，HNAVD88為赫爾默特正高，NGG為重力大地水準面；φp、λp分別為GNSS/水準點的大地緯度、經度；φ0、λ0分別為水準網幾何中心的大地緯度、經度；x1、x2分別為南北、東西方向水準面的傾斜度，是未知量。結合式(18)和式(20)計算了美國高程基準位差。表1中第2列為NAVD88高程基準與W0的重力位差，第3列為高程基準垂直差距，第4、5列分別為南北、東西方向的傾斜度。可以看出，CGG2010、GOCE+EGM08重力水準面之間較吻合，但USGG2012與其他二者相差了約3 cm。因此取三者均值，最終確定美國高程基準與全球高程基準重力位差為-4.82±0.05 m2s-2，垂直偏差為-49.2±0.5 cm。

圖1 NAVD88高程基準面系統偏差

重力水準面Ci0/m2s-2δN/cmx1/(m/°)x2/(m/°)CGG2010-4.71±0.05-48.08±0.5-0.0360.013USGG2012-5.01±0.05-51.13±0.5-0.0360.013GOCE+EGM08-4.73±0.05-48.29±0.5-0.0350.013

3 結論與建議

本文利用含有局部高程基準的“有偏”重力異常、局部高程基準的正高，按第三邊值問題推算了局部高程基準重力位差的計算方法。針對間接偏差項，利用衛星重力場模型和改化積分核減弱其影響。得出如下結論與建議：

(1) 大地邊值問題法確定高程異常差會引入間接偏差項。充分利用200階的GOCE衛星重力場，間接偏差項小于1 cm，可以忽略不計。

(2) 利用實測GNSS/水準數據確定美國高程基準NAVD88與全球高程基準差為-49.2±0.5 cm。

(3) EGM08重力場模型構建時，采用60階GRACE衛星重力場模型GGM02S作為其低階項，對無重力數據區域或數據專利區域，采用剩余地形模型來填充。因此EGM08在許多地方精度較低，而且60階的GRACE衛星重力場模型并不能完全消除間接偏差項的影響。因此建議，對于有地面重力數據的地區，應聯合GOCE和地面重力按邊值問題解算高精度的重力大地水準面，從而確定高程基準差；對于無地面重力數據的地區，可用GOCE+EGM08組合的重力大地水準面來確定高程基準差。

參考文獻：

[1] SACHER M，IHDE J，LIEBSCH G．EVRF2007 as Realization of the European Vertical Reference System[J]．Bollettino Di Geodesia E Scienze Affin，2009，68：35-50．

[2] WOODWORTH P，HUGHES C W，BINGHAM R J．Towards Worldwide Height System Unification Using Ocean Information[J]．Journal of Geodetic Science，2012，2(4)：302-318.DOI:10.2478/v10156-012-0004-8.

[3] 焦文海，魏子卿，馬欣，等．1985國家高程基準相對于大地水準面的垂直偏差[J]．測繪學報，2002，31(3):196-200．

[4] 郭海榮，焦文海，楊元喜．1985國家高程基準的系統差[J]．武漢大學學報(信息科學版)，2004，29(8)：715-717.

[5] 翟振和，魏子卿，吳富梅．利用EGM2008位模型計算中國高程基準與大地水準面間的垂直偏差[J]．大地測量與地球動力學，2011，31(4):116-118.

[6] 赫林，李建成．褚永海．1985國家高程基準與全球高程基準之間的垂直偏差[J]．測繪學報，2016，45(7):768-774．DOI:10.11947/j.AGCS.2016.20160029．

[7] 吳富梅，魏子卿．利用GNSS和EGM2008模型進行跨海高程傳遞[J]．武漢大學學報(信息科學版)，2016，41(5)：698-703.DOI:10.13203/j.whugis20140392.

[8] KOTSAKIS C．Estimation of the Zero-height Geopotential Level WLVD0 in a Local Vertical Datum From Inversion of Co-located GPS，Leveling and Geoid Heights:A Case Study in the Hellenic Islands[J]．Journal of Geodesy，2012，86(6):423-439．DOI：10.1007/s00190-011-0530-7.

[9] RUMMEL R，TEUNISSEN P．Height Datum Definition，Height Datum Connection and the Role of the Geodetic Boundary Value Problem[J]．Bulletin Géodésique，1988，62(4):477-498.

[10] SANSO F，USAI S．Height Datum and Local Geodetic Datums in the Theory of Geodetic Boundary Value Problems[J]．Allgemeine Vermessung Nachrichten，1995，102:343-355.

[11] PAN M，SJOBERG L E．Unification of Vertical Datums by GPS and Gravimetric Geoid Models With Application to Fennoscandia[J]．Journal of Geodesy，1998，72(2):64-70.

[12] 莫登沅．高程基準統一問題的研究[D]．鄭州：信息工程大學，2000.

[13] 徐菊生，王建華，賴錫安，等．垂直基準定義及聯接[J].地殼形變與地震，2001，21(1)：24-34.

[14] AMOS M J，FEATHERSTONE W E．Unification of New Zealand’s Local Vertical Datums: Iterative Gravimetric Quasigeoid Computations[J]．Journal of Geodesy，2009，83(1)：57-68.

[15] GERLACH C，RUMMEL R．Global Height System Unification with GOCE:A Simulation Study on the Indirect Bias Term in the GBVP Approach[J]．Journal of Geodesy，2013，87(1):57-67．DOI:10.1007/s00190-012-0579-y．

[16] GATTI A，REGUZZONI M，VENUTI G．The Height Datum Problem and the Role of Satellite Gravity Models[J]．Journal of Geodesy，2013，87(1):15-22．DOI：10.1007/s00190-012-0574-3.

[17] ZHANG L，LI F，CHEN W，et al．Height Datum Unification Between Shenzhen and Hong Kong Using the Solution of the Linearized Fixed-gravimetric Boundary Value Problem[J]．Journal of Geodesy，2008，83(5)：411-417．DOI:10.1007/s00190-008-0234-9．

[18] JEKELI C．Heights，the Geopotential，and Vertical Datums[C]∥Report 459.[S.l.]：The Ohio State University，2000．

[19] XU P．A Quality Investigation of Global Vertical Datum Connection[J]．Geophysical Journal International，1992，110(2)：361-370．DOI:10.1111/j.1365-246X.1992.tb00880.x．

[20] 傅露，褚永海．區域大地水準面確定中Stokes核函數的應用[J]．大地測量與地球動力學，2013，33(2):110-113．

[21] 赫林，李建成．褚永海．1985國家高程基準與全球高程基準之間的垂直偏差[J]．測繪學報，2016，45(7)：768-774．DOI:10.11947/j.AGCS.2016.20160029．

[22] 趙德軍．GOCE重力場模型的精度評估[J]．大地測量與地球動力學，2015，35(2)：21-25．

[23] HUANG J，VERONNEAU M．Canadian Gravimetric Geoid Model 2010[J]．Journal of Geodesy，2013，87(8)：771-790．