不同全球?qū)α鲗犹祉斞舆t產(chǎn)品在中國區(qū)域的比較

李博峰，王苗苗，沈云中，樓立志

（同濟大學(xué) 測繪與地理信息學(xué)院，上海 200092）

電離層與對流層誤差是GNSS（Global Navigation Satellite System）應(yīng)用的兩大主要誤差源，電離層誤差可采用多頻信號有效消除.盡管對流層誤差較電離層小，但在精密定位中無法通過多頻信號消除.對流層延遲包括干延遲和濕延遲兩部分，其中干延遲約占總延遲的90%，可采用Saastamoinen，Hopfield，UNB3等經(jīng)驗?zāi)Ｐ瓦M(jìn)行精確改正［1］.由水汽引起的濕延遲雖然較小，但其變化率是干延遲的3～4倍，因此對流層延遲是制約未來多頻多模GNSS模糊度固定、精密定位的主要因素［2－3］.

隨著對地觀測網(wǎng)絡(luò)的完善和觀測數(shù)據(jù)的增加，對流層延遲從依賴簡單封閉的數(shù)學(xué)模型向依賴大量外部數(shù)據(jù)的改正模型過渡［4］.目前已有諸多機構(gòu)提供全球?qū)α鲗犹祉斞舆t產(chǎn)品，這些產(chǎn)品采用的數(shù)據(jù)源、投影函數(shù)不同，因此它們的精度和適用區(qū)域都不盡相同.IGS（International GNSS Service）融合多個分析中心的結(jié)果發(fā)布全球?qū)α鲗犹祉斞舆t產(chǎn)品，其精度可達(dá)4mm，應(yīng)用范圍較廣［5］.維也納科技大學(xué)UTV（Vienna University of Technology）基于歐洲中尺度氣象預(yù)報中心的數(shù)值氣象資料ECMWF-NWM結(jié)合 Marini［6］投影函數(shù)發(fā)布了全球?qū)α鲗犹祉斞舆t產(chǎn)品，精度可達(dá)到2cm，其在歐洲的適用性較好.加拿大新布倫斯威克大學(xué)UNB（University of New Brunswick）基于美國國家環(huán)境中心和全球多尺度環(huán)境中心NCEP-GEM再分析資料結(jié)合VMF1［7］投影函數(shù)發(fā)布了全球?qū)α鲗犹祉斞舆t產(chǎn)品，其在北美區(qū)域的適用性較好.

不同產(chǎn)品采用不同的數(shù)據(jù)源，即使采用相同的投影函數(shù)計算對流層天頂延遲，全球不同區(qū)域其精度亦有差異.況且現(xiàn)有產(chǎn)品采用的數(shù)據(jù)源和投影函數(shù)都不盡相同，盡管在各自的數(shù)據(jù)源區(qū)域能取得良好的效果，但是鮮有研究分析這些產(chǎn)品在中國區(qū)域內(nèi)的精度差異.本文分析4種對流層天頂延遲產(chǎn)品及不同對流層高程歸算方法和格網(wǎng)插值方法在中國區(qū)域內(nèi)的精度.

1 對流層天頂延遲產(chǎn)品

1.1 VMF1產(chǎn)品

VMF1（Vienna Mapping Functions）格網(wǎng)產(chǎn)品將地球按照地心經(jīng)度、緯度等間隔地分成2.5°×2.0°的格網(wǎng)，每隔6h提供全球?qū)α鲗犹祉斞舆t格網(wǎng)產(chǎn)品，數(shù)據(jù)包括映射函數(shù)系數(shù)、天頂干延遲（zenith hydrostatic delay，ZHD）和天頂濕延遲（zenith wet delay，ZWD）.除了格網(wǎng)產(chǎn)品，同時還提供全球IGS站點的天頂延遲數(shù)據(jù)，包括映射函數(shù)系數(shù)、ZHD和ZWD以及站點氣象數(shù)據(jù)，其時間分辨率與格網(wǎng)產(chǎn)品相同.

圖1 UNBVMFG，UNBVMFGcmc和UNBVMFP與VMF1的ZTD差異分布Fig.1 Comparison of ZTD from UNBVMFG，UNBVMFGcmc and UNBVMFP with VMF1

VMF1產(chǎn)品精度較高，可達(dá)到2cm，可用于精密的靜態(tài)定位或者亞米級的定位.通過插值計算，格網(wǎng)對流層天頂延遲能夠獲得所有站點的天頂對流層延遲，因而比站點對流層延遲產(chǎn)品的應(yīng)用范圍更廣泛.

1.2 UNB產(chǎn)品

作為對VMF1產(chǎn)品的一個很好的補充，UNB提供3種全球?qū)α鲗犹祉斞舆t格網(wǎng)產(chǎn)品，即UNBVMFG， UNBVMFGcmc 和 UNBVMFP.UNBVMFG是利用NCEP-GEM再分析資料得出的產(chǎn)品；UNBVMFGcmc是利用全球確定性預(yù)報系統(tǒng)GDPS（Global Deterministic Prediction System）分析得出的產(chǎn)品，其6h和18h的數(shù)據(jù)以經(jīng)過6h的預(yù)測值為基礎(chǔ)；UNBVMFP是利用GDPS分析得出的預(yù)測產(chǎn)品，分別預(yù)測24，30，36，42h的對流層延遲值.3種產(chǎn)品的延遲分別為7d，1d和0d.圖1分別是3種UNB產(chǎn)品與VMF1產(chǎn)品在2012年6月1日12時的全球ZTD（zenith troposphere delay）差異分布圖.

3種UNB產(chǎn)品與VMF1產(chǎn)品在2012年6月1日12時的全球差異最小值分別為－174.8，－121.0，－133.2mm，最大值分別為191.9，140.5，164.9mm，平均值分別為7.21，7.56，9.69mm.在中國范圍內(nèi)大部分區(qū)域的差異值基本一致，約為20 mm，個別地區(qū)出現(xiàn)50mm以上的差異，說明各對流層天頂延遲產(chǎn)品在中國區(qū)域內(nèi)的差異性不盡相同，且在區(qū)域上呈現(xiàn)不連續(xù)性分布.

2 對流層天頂延遲

2.1 對流層延遲高程改正

格網(wǎng)點ZHD，ZWD參考地形格網(wǎng)點高度hg，站點高度為hs，首先需要將對流層延遲從格網(wǎng)高程歸算至站點高程，ZHD的高程歸算方法一般有2種.第1種方法較為簡單，只考慮對流層隨高程變化的物理機制（記為 HGT-CORR改正方法）［8］，即

式中：DH（hs），DH（hg）分別為站點和格網(wǎng)點的干延遲；hg，hs分別為格網(wǎng)點高度和站點高度，m；T（hg）為格網(wǎng)點對應(yīng)的溫度，K；P（hg）為格網(wǎng)點對應(yīng)的氣壓，hPa，溫度和氣壓可以由GPT（Global Pressure and Temperature）［9］模型或者標(biāo)準(zhǔn)大氣解算；R為氣體常數(shù)，R＝0.289 644KJ·kg－1；g為重力常數(shù)，g＝9.784m·s－2.

另一種ZHD高程歸算方法顧及氣壓的影響（記為PRESS-CORR改正方法），數(shù)學(xué)公式如下：

式中：P為氣壓，hPa；φ為緯度；h為高度，m；已知DH（hg），hg，hs和φ，首先根據(jù)式（3）分別計算氣壓值P（hs），P（hg），再根據(jù)已知DH（hg），φ，hg用式（2）反算格網(wǎng)點氣壓P′（hg），然后計算站點的氣壓修正值P′（hs）＝P（hs）＋P′（hg）－P（hg），最后將P′（hs）代入式（2）計算站點干延遲DH（hs）［10］.

由水汽引起的濕延遲ZWD變化率大、隨機性較強，為了完整性和便于直接比較，類似于干延遲在高程上的改正，選擇指數(shù)衰減函數(shù)對ZWD在高程上進(jìn)行改正［10］.

式中，DW（hg），DW（hs）分別是格網(wǎng)點和站點的天頂濕延遲.

2.2 對流層延遲格網(wǎng)插值

格網(wǎng)插值常采用插值方法簡單、以格網(wǎng)值近似地遵守線性變化為基礎(chǔ)的雙線性內(nèi)插法（BIL）.如圖2a所示，設(shè)

式中：0≤p＜1，0≤q＜1；λ和φ分別為插值點所在的經(jīng)度和緯度；λ00和φ00分別是插值點所在格網(wǎng)左下角點的經(jīng)度和緯度.可以先線性內(nèi)插出點A，B或者點C，D，然后再線性內(nèi)插出點P.點P的Z（λ，φ）值計算公式為

式中：Z0，0，Z1，0，Z0，1，Z1，1分別是點P所在格網(wǎng)各角點的值；Q00，Q10，Q01，Q11分別是對應(yīng)的因子.

圖2 周圍4點與8點插值Fig.2 Interpolation with nearby 4points and 8points

另一種常用的計算點P的Z（λ，φ）值的方法為 距離加權(quán)格網(wǎng)插值（WGI）［11］.不同于BIL方法，定義函數(shù)Q（p，q）為

則Q00，Q10，Q01，Q11分別為

還可采用雙二次多項式內(nèi)插方法（BIQ），其對由單個格網(wǎng)拼接成整體曲面造成的不光滑以及格網(wǎng)點處形成“尖點”和格網(wǎng)邊上形成“折痕”有所改善［12］.設(shè)平面上3點滿足二次多項式y(tǒng)＝ax2＋bx＋c，由3個已知點計算出多項式系數(shù)，從而可計算出任意點的對流層延遲.如圖2b所示，找到與插值點P距離最近的中心格網(wǎng)點（圖中點5），用中心點周圍的點1～9對點P進(jìn)行插值.由點1～3插值點A，點4～6插值點B，點7～9插值點C，然后由A，B，C三點插值計算點P.

3 試驗分析

綜合格網(wǎng)插值方法和干延遲在高程上的歸算方法，組成計算方案1～6如下：BIL與HGT-CORR組合；BIL 與 PRESS-CORR 組 合；WGI與 HGTCORR組合；WGI與PRESS-CORR組合；BIQ與HGT-CORR組合；BIQ與PRESS-CORR組合.

各方案采用的濕延遲改正方式相同.采用中國區(qū)域內(nèi)9個IGS站從2012年6月1日至9月1日共92d的數(shù)據(jù)，用6種計算方案計算比較4種全球?qū)α鲗犹祉斞舆t格網(wǎng)產(chǎn)品和干延遲高程歸算方法及格網(wǎng)插值方法在中國不同區(qū)域的差異性.試驗中，內(nèi)插站點的ZTD以IGS提供的產(chǎn)品值為參考，ZHD，ZWD分別以VMF1站點產(chǎn)品的ZHD，ZWD為參考.表1是9個IGS站點的基本信息.

表中各站點的高程離散度Δh按式（10）計算：

式中：hi為站點所在格網(wǎng)第i個格網(wǎng)點的高程；h為站點的高程.站點高程離散度用于反映內(nèi)插站點高程相對于格網(wǎng)點高程的相對變化幅度.

首先分析VMF1對流層格網(wǎng)產(chǎn)品.對VMF1對流層格網(wǎng)產(chǎn)品分別采用6種方案計算9個IGS站92 d的ZTD，然后分別統(tǒng)計它們的內(nèi)符合精度（STD）和外符合精度（RMS），如圖3所示.除此之外，還分析了不同方案ZHD和ZWD的內(nèi)插效果，其統(tǒng)計結(jié)果如表2，其中方案1，3，5的RMS統(tǒng)計不包括LHAZ站.

圖3 VMF1產(chǎn)品的內(nèi)符合精度與外符合精度Fig.3 STD and RMS of ZTD from VMF1

除了拉薩站LHAZ外，各站點6種方案獲得的ZTD，ZHD和ZWD的精度基本一致，外符合精度分別小于50，20，40mm；方案間相差最大分別不超過10，5，8mm；所有站點平均偏差值均分別在－40～30mm，－15～20mm，－20～25mm之間；所有站點的內(nèi)符合精度均分別小于35.0，3.2，32.0mm，即對流層延遲在高程上的歸算方法和插值方法在站點上差異不明顯.相同的方案在不同的區(qū)域計算的精度有差異，外符合精度差異最大分別約為30，16，30 mm，說明不同的區(qū)域影響對流層延遲在高程上的歸算方法和插值方法的精度.站點內(nèi)符合精度之間的差異說明產(chǎn)品在不同區(qū)域的穩(wěn)定性不同.

表2 6種方案計算的內(nèi)外符合精度統(tǒng)計Tab.2 STD and RMS from six interpolation schemes mm

對流層天頂干延遲的歸算精度與內(nèi)插點的高程離散度相關(guān)，對LHAZ站而言，其高程離散度為1 614.57m，遠(yuǎn)大于其他IGS站，所以干延遲高程歸算方法在LHAZ產(chǎn)生約145mm的精度差異，而在其他站點差異不明顯.因此，LHAZ站較大的高程離散度是導(dǎo)致2種干延遲高程歸算方法（HGT-CORR和PRESS-CORR）產(chǎn)生較大差異的原因，并且HGT-CORR方法對高程離散度更為敏感，在高程離散度較大時，改正精度較低；而PRESS-CORR方法能有效克服HGT-CORR的缺陷，即使對高程離散度較大的LHAZ站也能取得較好效果.當(dāng)高程離散度較小時，2種干延遲高程歸算方法都能取得良好的效果.

類似于對VMF1產(chǎn)品的分析，對3種UNB對流層產(chǎn)品也做了分析，得到類似的結(jié)論.

由上述分析知，6種方案在各站點獲得的對流層延遲精度基本一致，現(xiàn)選擇方案4分析4種對流層產(chǎn)品.采用方案4分別計算4種格網(wǎng)產(chǎn)品獲得各站點的對流層延遲，統(tǒng)計其內(nèi)外符合精度，其中ZTD的內(nèi)外符合精度如圖4，對應(yīng)的統(tǒng)計結(jié)果如表3.

圖4 4種產(chǎn)品計算ZTD的內(nèi)符合精度和外符合精度Fig.4 STD and RMS of ZTD from four different products

表3 4種產(chǎn)品計算的ZTD的內(nèi)符合精度與外符合精度Tab.3 STD and RMS of ZTD from four different products mm

采用方案4計算UNB提供的3種對流層天頂延遲產(chǎn)品獲得的ZTD，ZHD和ZWD的精度在區(qū)域上的變化基本一致，差值最大分別約為25，3，25mm，UNBVMFG產(chǎn)品在中國區(qū)域精度相對較低（除BJFS站外），而UNBVMFGcmc產(chǎn)品在中國區(qū)域的精度相對較高.在GUAO，BJFS，WUHN站處，VMF1產(chǎn)品與3種UNB產(chǎn)品在區(qū)域上的變化有差異，并且各產(chǎn)品在區(qū)域上的精度變化也有不同.分析原因可能是由于2個機構(gòu)提供的對流層天頂延遲產(chǎn)品采用的數(shù)據(jù)源、映射函數(shù)和計算方法不同引起產(chǎn)品的精度變化在區(qū)域上呈現(xiàn)差異性.

4種對流層天頂延遲產(chǎn)品獲得的站點ZTD，ZHD和ZWD的內(nèi)符合精度STD在區(qū)域上的變化基本一致（UNBVMFG在區(qū)域上的變化幅度稍大），差值最大分別約為25，3，30mm.UNBVMFG產(chǎn)品在中國區(qū)域的穩(wěn)定性相對較差，VMF1產(chǎn)品的穩(wěn)定性相對較高.

4 結(jié)論

選取中國區(qū)域內(nèi)9個IGS站點，采用6種方案分析了4種全球?qū)α鲗犹祉斞舆t格網(wǎng)產(chǎn)品VMF1，UNBVMFG，UNBVMFGcmc，UNBVMFP，分析結(jié)果表明，當(dāng)內(nèi)插點的高程離散度較大時，單純以修正高程歸算對流層干延遲的方式對高程離散度更加敏感，高程離散度大時獲得的精度較低，而氣壓修正的方式可以克服該缺陷.區(qū)域差異影響對流層延遲在高程上的歸算方法和插值方法的精度.各產(chǎn)品在中國區(qū)域內(nèi)的精度和穩(wěn)定性及其在區(qū)域上的變化不盡相同，但是各產(chǎn)品都可以靈活地為導(dǎo)航、定位、水汽反演等應(yīng)用提供一定精度的對流層天頂延遲.

［1］ 楊玲，李博峰，樓立志.不同對流層模型對GPS定位結(jié)果的影響［J］.測繪通報，2009，4（3）：9.

YANG Ling，LI Bofeng，LOU Lizhi.Effects of different troposphere correction models on GPS positioning［J］.Bulletin of Surveying and Mapping，2009，4（3）：9.

［2］ 李博峰.混合整數(shù)GPS線性模型的參數(shù)與方差－協(xié)方差分量估計理論與方法［D］.上海：同濟大學(xué)，2010.

LI Bofeng.Theory and method for parameter and variancecovariance component estimation in mixed integer linear GPS model［D］.Shanghai：Tongji University，2010.

［3］ LI Bofeng，F(xiàn)ENG Yanming，SHEN Yunzhong，etal.Geometry-specified troposphere decorrelation for subcentimeter real-time kinematic solutions over long baselines［J］.Journal of Geophysical Research，2010：DOI：10.1029／2010JB007549.

［4］ Urquhart L，Santos M.Development of a VMF1-like service at UNB［D］：Fredericton：Department of Geodesy and Geomatics Engineering of University of New Brunswick，2011.

［5］ Byun S H，Bar-Sever Y E.A new type of troposphere zenith path delay product of the international GNSS service［J］.Journal of Geodesy，2009，83（314）：1.

［6］ Marini J.Correction of satellite tracking data for an arbitrary tropospheric profile［J］.Radio Science，1972，7（2）：223.

［7］ Boehm J，Werl B，Schuh H.Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for medium-range weather forecasts operational analysis data［J］.Journal of Geophysical Research，2006：DOI：10.1029／2005JB003629.

［8］ Steigenberger P，Boehm J，Tesmer V.Comparison of GMF／GPT with VMF1／ECMWF and implications for atmospheric loading［J］.Journal of Geodesy，2009，83（10）：943.

［9］ Boehm J，Heinkelmann R，Schuh H.Short note：aglobal model of pressure and temperature for geodetic applications［J］.Journal of Geodesy，2007，81（10）：679.

［10］ Kouba J.Implementation and testing of the gridded Vienna Mapping Function 1（VMF 1）［J］.Journal of Geodesy，2008，82（4）：193.

［11］ Chao YC. Real Time Implemental of the wide area augumention system for the Global Position System with an emphasis on ionospheric modeling［D］.Stanford：Stanford University，1997.

［12］ 鄧興升，郭云開，花向紅.似大地水準(zhǔn)面格網(wǎng)雙二次多項式插值方法［J］.測繪學(xué)報，2009，38（1）：35.

DENG Xingsheng， GUO Yunkai， HUA Xianghong.Quasigeoid grid network bi-quadratic interpolation method［J］.ActaGeodaeticaetCartographicaSinica，2009，38（1）：35.