趙德軍，徐新強(qiáng)，歐陽明達(dá)

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

眾所周知，當(dāng)?shù)仄骄Ｋ娌⒉恢睾嫌谌虼蟮厮疁?zhǔn)面，它們的差異稱為海面地形。全球海面地形最大差異可達(dá)2 m。例如，我國東海地區(qū)各長期驗(yàn)潮站的平均海面比青島驗(yàn)潮站代表的國家高程基準(zhǔn)要高出25 cm左右，而南海驗(yàn)潮站確定的平均海面要高出1985國家高程基準(zhǔn)35 cm以上。因此建立在當(dāng)?shù)仄骄Ｋ嫔系母鲊叱袒鶞?zhǔn)就不一樣。目前全世界存在超過100個高程基準(zhǔn)，如何統(tǒng)一它們成為大地測量、海洋學(xué)界努力的目標(biāo)。

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

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

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

1 大地邊值問題法

如上文所述，無論基于第二、還是第三大地邊值問題，統(tǒng)一高程基準(zhǔn)的關(guān)鍵是要確定獨(dú)立于任何高程基準(zhǔn)的全球統(tǒng)一的大地水準(zhǔn)面或似大地水準(zhǔn)面。

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

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

(1)

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

(2)

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

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

令

(9)

式中，Nzero為大地水準(zhǔn)面0階項(xiàng)，目前能精確確定，可作為已知量。令

(10)

式中，NStokes為地面重力異常計(jì)算的大地水準(zhǔn)面，可作為已知量。令

(11)

式中，Nind稱為間接偏差項(xiàng)，其中

(12)

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

(13)

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

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

1.3 聯(lián)合衛(wèi)星和地面重力的邊值問題解

(14)

(15)

式中，Pn(cosψ)為勒讓德級數(shù)。改化Stokes核函數(shù)的作用為加速積分收斂，減小球冠Ω0積分面積，同時間接偏差項(xiàng)變?yōu)?/p>

(16)

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

(17)

最后給出高程基準(zhǔn)重力位差的計(jì)算式為

Cl0=γ(h-Hl-N0)

(18)

正常高系統(tǒng)中，重力位差的計(jì)算式為

(19)

2 數(shù)值試驗(yàn)

2.1 GNSS/水準(zhǔn)數(shù)據(jù)

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

2.2 重力大地水準(zhǔn)面

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

2.3 重力位差

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

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

(20)

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

圖1 NAVD88高程基準(zhǔn)面系統(tǒng)偏差

重力水準(zhǔn)面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 結(jié)論與建議

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

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

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

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

參考文獻(xiàn)：

[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國家高程基準(zhǔn)相對于大地水準(zhǔn)面的垂直偏差[J]．測繪學(xué)報，2002，31(3):196-200．

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

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

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

[7] 吳富梅，魏子卿．利用GNSS和EGM2008模型進(jìn)行跨海高程傳遞[J]．武漢大學(xué)學(xué)報(信息科學(xué)版)，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] 莫登沅．高程基準(zhǔn)統(tǒng)一問題的研究[D]．鄭州：信息工程大學(xué)，2000.

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

[14] AMOS M J，F(xiàn)EATHERSTONE 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] 傅露，褚永海．區(qū)域大地水準(zhǔn)面確定中Stokes核函數(shù)的應(yīng)用[J]．大地測量與地球動力學(xué)，2013，33(2):110-113．

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

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

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