何亞軍，陳曉茜

（1.中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司，湖北 武漢 430063；2.湖北省基礎(chǔ)地理信息中心（湖北省北斗衛(wèi)星導(dǎo)航應(yīng)用技術(shù)研究院），湖北 武漢 430074）

地球重力場(chǎng)反映了地球內(nèi)部的物質(zhì)分布和地球旋轉(zhuǎn)運(yùn)動(dòng)信息，作為基本的地學(xué)信息，其在地球物理學(xué)、現(xiàn)代大地測(cè)量學(xué)、地球動(dòng)力學(xué)等學(xué)科研究中有著重要的作用。大地測(cè)量中的地球重力場(chǎng)主要以截?cái)嗟揭欢A次的球諧位系數(shù)來(lái)表示，稱為地球重力場(chǎng)模型。當(dāng)前高精度高分辨率地球重力場(chǎng)模型的確定在國(guó)際上廣受重視，已有多個(gè)國(guó)家發(fā)射專門用于重力探測(cè)的重力衛(wèi)星進(jìn)行地球重力場(chǎng)測(cè)量，如德國(guó)的CHAMP衛(wèi)星、美德的GRACE衛(wèi)星、歐洲的GOCE衛(wèi)星以及美德的GRACE-FO衛(wèi)星等，我國(guó)重力衛(wèi)星也在籌備階段。重力衛(wèi)星主要采用衛(wèi)星地面跟蹤技術(shù)、衛(wèi)星跟蹤衛(wèi)星技術(shù)以及衛(wèi)星重力梯度技術(shù)測(cè)定地球重力場(chǎng)，并結(jié)合地面重力、航空重力觀測(cè)資料和衛(wèi)星測(cè)高資料反演出較高精度的地球重力場(chǎng)模型[1-3]。

目前，國(guó)際地球模型中心ICGEM已收錄靜態(tài)地球重力場(chǎng)模型176個(gè)，這些重力場(chǎng)模型在大地測(cè)量和工程測(cè)量方面有著廣泛應(yīng)用，比如利用地球重力場(chǎng)模型可計(jì)算GNSS點(diǎn)高程異常，從而將GNSS大地高轉(zhuǎn)換為工程測(cè)量所需的正常高，減少水準(zhǔn)測(cè)量的工作；其次利用地球重力場(chǎng)模型可計(jì)算水準(zhǔn)測(cè)量中的正常水準(zhǔn)面不平行改正和重力異常改正，提高水準(zhǔn)測(cè)量精度，滿足海拔變化劇烈地區(qū)的水準(zhǔn)測(cè)量需求[4]。本文利用中國(guó)某區(qū)域的GNSS點(diǎn)數(shù)據(jù)，利用多個(gè)地球重力場(chǎng)模型計(jì)算各點(diǎn)高程異常，并與各點(diǎn)實(shí)測(cè)高程異常相比較，檢驗(yàn)多類地球重力場(chǎng)模型的精度，為工程測(cè)量實(shí)踐工作提供參考和借鑒。

1 原理與方法

1.1 地球重力場(chǎng)模型高程異常計(jì)算

地球重力場(chǎng)的引力位是一個(gè)逼近地球質(zhì)體外部引力位在無(wú)窮遠(yuǎn)處收斂到零值的調(diào)和函數(shù)，通常展開為一個(gè)在理論上收斂的整階次球諧或橢球諧函數(shù)的無(wú)窮級(jí)數(shù)，這個(gè)級(jí)數(shù)展開系數(shù)的集合定義一個(gè)相應(yīng)的地球重力場(chǎng)模型。利用地球重力場(chǎng)模型根據(jù)擾動(dòng)位原理計(jì)算高程異?？杀硎緸閇5-7]

式中，GM為地球引力常數(shù)；r為地心向徑；a為參考橢球長(zhǎng)半軸（a=6 378 137.0 m，采用CGCS2000參考橢球，以下橢球參數(shù)均相同）；和為減去正常橢球重力場(chǎng)后的完全規(guī)格化球諧位系數(shù)；(sinφ)為完全規(guī)格化締合勒讓德函數(shù)（可采用跨階次遞推法）；φ為地心緯度；λ為地心經(jīng)度；γG為GNSS點(diǎn)處的正常重力。

在確定赤道處的正常重力γEQ后，根據(jù)熟知的Somigliana公式，即可計(jì)算橢球表面大地緯度為B處的正常重力γSF，可表示為[8]

式中，K為理論正常重力公式常數(shù)（K=0.101 931 852 619 31）；e為參考橢球第一偏心率（e=0.081 819 191 042 816）；γEQ代表赤道處正常重力。

式中，e′為第二偏心率（e′=0.082 094 438 151 917）。q0、q′0、m的計(jì)算公式分別為

式中，x、y、z代表計(jì)算點(diǎn)的地固直角坐標(biāo)；u為橢球坐標(biāo)中的共焦橢球短半軸；b為參考橢球短半軸（b=6 353 752.314 1 m）；E=（a2-b2）為線偏心率；ω為參考橢球角速度。

GNSS點(diǎn)處正常重力γG計(jì)算公式為：

式中，H為GNSS點(diǎn)的大地高；f為參考橢球扁率。m′的計(jì)算公式為：

1.2 GNSS點(diǎn)實(shí)測(cè)高程異常計(jì)算

高程異常是似大地水準(zhǔn)面與參考橢球面之間的高差。目前GNSS測(cè)量技術(shù)可以測(cè)定高精度的大地高，水準(zhǔn)或三角高程測(cè)量精確測(cè)定點(diǎn)位正常高，即可計(jì)算該點(diǎn)的實(shí)測(cè)高程異常。GNSS實(shí)測(cè)點(diǎn)高程異?？杀硎緸椋?/p>

式中，H為該點(diǎn)的大地高；h為該點(diǎn)的正常高。

GNSS測(cè)定的大地高及水準(zhǔn)測(cè)量的正常高已能達(dá)到毫米級(jí)精度，計(jì)算的高程異常即可達(dá)到較高精度。而直接利用國(guó)際地球重力場(chǎng)模型解算GNSS點(diǎn)高程異常并未能達(dá)到如此高的精度，首先其所基于的橢球面或似大地水準(zhǔn)面與我國(guó)高程系統(tǒng)存在一定的系統(tǒng)偏差，且重力場(chǎng)模型在計(jì)算中需截?cái)嘀烈欢A次，本身也存在一定誤差。因此，為了提高重力場(chǎng)模型解算高程異常的精度，應(yīng)加入系統(tǒng)偏差加以改善，系統(tǒng)偏差計(jì)算采用分區(qū)常系統(tǒng)偏差矯正法和直接計(jì)算法，計(jì)算公式分別為式（11）和式（12）。

式中，ξi0為i分區(qū)中的計(jì)算的常系統(tǒng)偏差值；ξij為i分區(qū)中的第j個(gè)實(shí)測(cè)高程異常值與模型值的差值；n為i分區(qū)中的GNSS點(diǎn)數(shù)。

式中，ξ0為GNSS點(diǎn)實(shí)測(cè)高程異常；ξM為GNSS點(diǎn)重力場(chǎng)模型高程異常；i為點(diǎn)號(hào)，n為GNSS點(diǎn)的總數(shù)目。

1.3 精度評(píng)定標(biāo)準(zhǔn)

采用相關(guān)系數(shù)（r）、均方根誤差（RMSE）、納什系數(shù)（NSE）作為評(píng)定標(biāo)準(zhǔn)，對(duì)模型高程異常和實(shí)測(cè)高程異常差值的精度和兩者的相關(guān)性以及重力場(chǎng)模型的解算質(zhì)量進(jìn)行檢驗(yàn)。相關(guān)系數(shù)用于檢驗(yàn)兩者的相關(guān)程度；均方根誤差確定未加入系統(tǒng)偏差和加入系統(tǒng)偏差的高程異常差值精度，檢驗(yàn)系統(tǒng)偏差對(duì)結(jié)果精度的改善效果；納什系數(shù)檢驗(yàn)重力場(chǎng)模型解算高程異常值的整體質(zhì)量，具體計(jì)算公式為：

式中，Mi為模型高程異常值；Oi為實(shí)測(cè)高程異常值；和分別代表所有GNSS點(diǎn)模型高程異常和實(shí)測(cè)高程異常的平均值。

2 實(shí)例分析

選擇安徽省至江蘇省某鐵路安徽段的GNSS點(diǎn)數(shù)據(jù)作為研究資料，共有70個(gè)GNSS/水準(zhǔn)點(diǎn)，均采用CGCS2000國(guó)家大地坐標(biāo)系參考橢球體，高程系統(tǒng)采用1985國(guó)家高程系統(tǒng)，如圖1所示。

圖1 GNSS/水準(zhǔn)點(diǎn)線路分布

選取了10個(gè)精度較好的地球重力場(chǎng)模型進(jìn)行模型高程異常的計(jì)算，各個(gè)模型的公布年份、數(shù)據(jù)來(lái)源、研制者以及最高階數(shù)見表1。

表1 地球重力場(chǎng)模型對(duì)比

根據(jù)圖1顯示，GNSS/水準(zhǔn)點(diǎn)經(jīng)緯度變化范圍較小，且呈現(xiàn)規(guī)律性，故采用分區(qū)常系統(tǒng)偏差矯正法時(shí)，將其按照一個(gè)區(qū)域進(jìn)行分區(qū)進(jìn)行系統(tǒng)偏差，具體精度結(jié)果見表2。采用分區(qū)常系統(tǒng)偏差矯正法消除系統(tǒng)誤差后，GECO模型解算精度最高，高達(dá)3.94 cm，相比較未降低系統(tǒng)偏差的結(jié)果精度提高了81.32%。整體上看，各個(gè)模型的高程異常精度提高顯著，系統(tǒng)偏差改善效果明顯。

表2 分區(qū)常系統(tǒng)偏差矯正法降低系統(tǒng)偏差精度統(tǒng)計(jì)結(jié)果/cm

表3給出了各模型采用直接計(jì)算法解算的高程異常值與GNSS/水準(zhǔn)點(diǎn)實(shí)測(cè)高程異常的比較結(jié)果，在降低系統(tǒng)偏差后，XGM2019e_2159、EGM2008、EIGEN-6C4、GECO、SGG-UGM-1、EIGEN-6C3stat 6個(gè)高階重力場(chǎng)模型中，EIGEN-6C4模型解算精度最高，為15.52 cm；XGM2016、EGM96、GGM05C、PGM2000a4個(gè)低階重力場(chǎng)模型中，PGM2000a模型解算精度最高，高達(dá)7.22 cm。同時(shí)整體上，在降低系統(tǒng)偏差后，各個(gè)模型的高程異常精度均有所提高，其中高階重力場(chǎng)模型提高幅度相較低階重力場(chǎng)模型更為顯著，提升范圍為16%～21%。但直接計(jì)算法降低系統(tǒng)偏差后的精度結(jié)果相比較分區(qū)常系統(tǒng)偏差矯正法的精度結(jié)果，解算精度數(shù)值低，系統(tǒng)誤差并沒有完全消除，仍然需要改善。

表3 直接計(jì)算法降低系統(tǒng)偏差精度統(tǒng)計(jì)結(jié)果/cm

對(duì)比兩種系統(tǒng)偏差計(jì)算法的精度結(jié)果，分區(qū)常系統(tǒng)偏差矯正法的精度提升效果更為明顯，模型解算值更為接近真實(shí)高程異常值。

利用相關(guān)系數(shù)（r）、納什系數(shù)（NSE）可計(jì)算出各個(gè)模型的精度比較結(jié)果，檢驗(yàn)?zāi)Ｐ透叱坍惓Ｅc實(shí)測(cè)高程異常的符合程度，具體計(jì)算結(jié)果見圖2、3。圖2、3結(jié)果表明采用兩種降低系統(tǒng)偏差方法下，各個(gè)模型的相關(guān)系數(shù)均在0.95以上，說明模型高程異常與實(shí)測(cè)高程異常具有高度相關(guān)關(guān)系。此外，各個(gè)模型的納什系數(shù)均在0.8以上，說明重力場(chǎng)模型解算結(jié)果可信度高，模型質(zhì)量好。對(duì)比圖2、3，可比較出兩種降低系統(tǒng)偏差方法計(jì)算出的各模型結(jié)果的相關(guān)系數(shù)極為接近，而納什系數(shù)在數(shù)值上，分區(qū)常系統(tǒng)偏差矯正法明顯優(yōu)于直接計(jì)算法，更適合于重力場(chǎng)模型解算高精度高程異常。

圖2 分區(qū)常系統(tǒng)偏差矯正法精度評(píng)定結(jié)果圖

考慮到實(shí)際工作中GNSS/水準(zhǔn)點(diǎn)的數(shù)據(jù)存在不完整情況，故選擇單個(gè)點(diǎn)位模型高程異常值與實(shí)測(cè)高程異常值的差值直接作為系統(tǒng)偏差改善模型解算精度，考察該種方法的解算精度，檢驗(yàn)其是否用于今后的工程實(shí)踐中。由于分區(qū)常系統(tǒng)偏差矯正法降低系統(tǒng)偏差效果更好，故選擇該方法的RMSE精度值作為參考值，限定精度差值為3 cm進(jìn)行計(jì)算，確定滿足條件的GNSS/水準(zhǔn)點(diǎn)數(shù)目，檢驗(yàn)總點(diǎn)位的符合情況，具體計(jì)算結(jié)果如表4。據(jù)表4結(jié)果可知，整體上各個(gè)模型計(jì)算的符合條件的點(diǎn)位數(shù)比較理想，最低個(gè)數(shù)為32，最高個(gè)數(shù)為57，解算精度效果較為理想。在點(diǎn)位實(shí)測(cè)數(shù)據(jù)缺失的情況下，可選擇合適點(diǎn)位模型值與實(shí)測(cè)高程異常值的差值作為系統(tǒng)偏差，提高模型解算精度，減少測(cè)量工作任務(wù)量。特別在項(xiàng)目初步設(shè)計(jì)即草測(cè)或初測(cè)階段，甚至在定測(cè)階段，線路方案不穩(wěn)定，比較方案較多，在只進(jìn)行優(yōu)選方案的水準(zhǔn)測(cè)量及所有方案的GNSS測(cè)量情況下，選取分區(qū)常系統(tǒng)偏差矯正法精度最優(yōu)且該模型下計(jì)算的符合條件的點(diǎn)位數(shù)較多的地球重力場(chǎng)模型，即可計(jì)算出滿足精度要求的比較方案的GNSS點(diǎn)的正常高，大大減少了水準(zhǔn)測(cè)量工作量?；蛟陧?xiàng)目初步設(shè)計(jì)階段，收集線路沿線一定數(shù)量有正常高及CGCS2000國(guó)家大地測(cè)量成果，在只測(cè)量框架網(wǎng)（CP0）及基礎(chǔ)平面控制網(wǎng)（CPI）的情況下，選取分區(qū)常系統(tǒng)偏差矯正法精度最優(yōu)的地球重力場(chǎng)模型，計(jì)算所有GNSS點(diǎn)的正常高，大大縮減了勘測(cè)前期精測(cè)網(wǎng)準(zhǔn)備階段工期及省去了繁重的水準(zhǔn)測(cè)量工作。

表4 各種模型計(jì)算符合條件點(diǎn)位數(shù)統(tǒng)計(jì)結(jié)果/個(gè)

圖3 直接計(jì)算法精度評(píng)定結(jié)果圖

3 結(jié) 語(yǔ)

根據(jù)各個(gè)模型解算的該鐵路GNSS點(diǎn)的高程異常結(jié)果可得，各個(gè)重力場(chǎng)模型計(jì)算的相關(guān)系數(shù)、均方根誤差以及納什系數(shù)三項(xiàng)精度評(píng)定值均符合限值范圍，其中GECO模型高程異常解算值降低系統(tǒng)偏差后精度值高達(dá)3.94 cm，相關(guān)系數(shù)和納什系數(shù)也極為理想，最為符合實(shí)測(cè)高程異常。利用地球重力場(chǎng)模型解算GNSS點(diǎn)高程異常相比較傳統(tǒng)的人力測(cè)量有著快速、簡(jiǎn)便、工作量小等諸多優(yōu)勢(shì)，可大大縮減測(cè)量工期及節(jié)省測(cè)量所需的人力、物力資源。伴隨著我國(guó)未來(lái)也將發(fā)射重力衛(wèi)星，將可得到精度和分辨率更高的地球重力場(chǎng)模型，同時(shí)利用多個(gè)重力場(chǎng)模型進(jìn)行組合計(jì)算，可進(jìn)一步提高高程異常的精度，未來(lái)單獨(dú)使用地球重力場(chǎng)模型解算高程異常的精度將有望達(dá)到四等或三等水準(zhǔn)測(cè)量，在工程測(cè)量方面發(fā)揮出更加重要的作用。