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

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

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

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

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

1 對(duì)流層天頂延遲產(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)品，同時(shí)還提供全球IGS站點(diǎn)的天頂延遲數(shù)據(jù)，包括映射函數(shù)系數(shù)、ZHD和ZWD以及站點(diǎn)氣象數(shù)據(jù)，其時(shí)間分辨率與格網(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)定位或者亞米級(jí)的定位.通過(guò)插值計(jì)算，格網(wǎng)對(duì)流層天頂延遲能夠獲得所有站點(diǎn)的天頂對(duì)流層延遲，因而比站點(diǎn)對(duì)流層延遲產(chǎn)品的應(yīng)用范圍更廣泛.

1.2 UNB產(chǎn)品

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

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

2 對(duì)流層天頂延遲

2.1 對(duì)流層延遲高程改正

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

式中：DH（hs），DH（hg）分別為站點(diǎn)和格網(wǎng)點(diǎn)的干延遲；hg，hs分別為格網(wǎng)點(diǎn)高度和站點(diǎn)高度，m；T（hg）為格網(wǎng)點(diǎn)對(duì)應(yīng)的溫度，K；P（hg）為格網(wǎng)點(diǎn)對(duì)應(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）分別計(jì)算氣壓值P（hs），P（hg），再根據(jù)已知DH（hg），φ，hg用式（2）反算格網(wǎng)點(diǎn)氣壓P′（hg），然后計(jì)算站點(diǎn)的氣壓修正值P′（hs）＝P（hs）＋P′（hg）－P（hg），最后將P′（hs）代入式（2）計(jì)算站點(diǎn)干延遲DH（hs）［10］.

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

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

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

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

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

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

圖2 周?chē)?點(diǎn)與8點(diǎn)插值Fig.2 Interpolation with nearby 4points and 8points

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

則Q00，Q10，Q01，Q11分別為

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

3 試驗(yàn)分析

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

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

表中各站點(diǎn)的高程離散度Δh按式（10）計(jì)算：

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

［1］ 楊玲，李博峰，樓立志.不同對(duì)流層模型對(duì)GPS定位結(jié)果的影響［J］.測(cè)繪通報(bào)，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線(xiàn)性模型的參數(shù)與方差－協(xié)方差分量估計(jì)理論與方法［D］.上海：同濟(jì)大學(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］ 鄧興升，郭云開(kāi)，花向紅.似大地水準(zhǔn)面格網(wǎng)雙二次多項(xiàng)式插值方法［J］.測(cè)繪學(xué)報(bào)，2009，38（1）：35.

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