范昊鵬，孫中苗，張麗萍，劉曉剛

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

對(duì)于甚長(zhǎng)基線干涉測(cè)量(very long baseline interferometry，VLBI)、全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system，GNSS)等空間大地測(cè)量技術(shù)，對(duì)流層延遲是其主要誤差源之一。在利用這些測(cè)量技術(shù)求解大地測(cè)量參數(shù)(包括測(cè)站坐標(biāo)、速度、地球定向參數(shù)等)的過程中，一般做法是：將對(duì)流層延遲中的主項(xiàng)，即天頂干延遲部分作為延遲的先驗(yàn)值；因干延遲計(jì)算精度較高，可達(dá)毫米級(jí)，遠(yuǎn)小于天頂濕延遲值，所以可將天頂濕延遲作為干延遲的改正值估計(jì)得出；最終，干、濕天頂延遲分別與對(duì)應(yīng)的映射函數(shù)相乘并求和得到斜路徑延遲[1-2]。以下如不作特殊說明，統(tǒng)一將天頂干、濕延遲簡(jiǎn)稱為干、濕延遲。

這種估計(jì)策略并未考慮映射函數(shù)的誤差，如果選用GMF(global mapping function)、NMF(Niell mapping function)等經(jīng)驗(yàn)映射函數(shù)模型，則干、濕分量有可能存在10-2、10-1量級(jí)的誤差[3-4]，此時(shí)與干、濕延遲相乘會(huì)產(chǎn)生厘米級(jí)甚至分米級(jí)的斜路徑延遲誤差，進(jìn)而對(duì)測(cè)站位置等大地測(cè)量參數(shù)的解算造成不可忽略的影響[5]。雖然VMF(Vienna mapping function)系列模型[6-7]相比GMF和NMF具有更高的精度和時(shí)間分辨率，但VMF模型需要利用氣象數(shù)據(jù)事后反算得到，因而在實(shí)際觀測(cè)中無法實(shí)時(shí)獲取。這對(duì)未來大地測(cè)量產(chǎn)品的快速、高精度發(fā)布較為不利[8]。

本文以GMF經(jīng)驗(yàn)?zāi)Ｐ蜑槔?，分析了該映射函?shù)的誤差特點(diǎn)，利用統(tǒng)計(jì)方法構(gòu)建了該誤差隨高度角變化的函數(shù)模型，進(jìn)而提出了顧及映射函數(shù)誤差的對(duì)流層延遲兩步估計(jì)法，并用試驗(yàn)證明了該方法的有效性。

1 傳統(tǒng)估計(jì)方法

傳統(tǒng)估計(jì)方法中，對(duì)流層延遲的數(shù)學(xué)模型如下

L=ZHD·mfh(el)+ZWD·mfw(el)+

maz(el)·[Gnscos(α)+Gewsin(α)]

(1)

式中，L表示斜路徑對(duì)流層延遲；ZHD、ZWD分別表示天頂對(duì)流層干、濕延遲；mfh(el)、mfw(el)表示高度角為el時(shí)干、濕分量的映射函數(shù)；maz(el)表示高度角為el時(shí)對(duì)流層水平梯度映射函數(shù)，具體形式為[9]

(2)

式中，α為方位角；Gns、Gew分別被定義為水平梯度的南北、東西分量。

此處需說明的是：①一般情況下，干延遲可按照Saastamonien模型計(jì)算[10]，當(dāng)已知地表氣壓時(shí)，認(rèn)為干延遲誤差可忽略；其中，地表氣壓可通過氣象儀器觀測(cè)獲取，或利用數(shù)值氣象模型等方法得到。本文采用的是后者，具體是歐洲中期天氣預(yù)報(bào)中心(European Centre for Medium-Range Weather Forecasts，ECMWF)提供的分層格網(wǎng)氣象數(shù)據(jù)。除此之外，后文用于評(píng)定精度的映射函數(shù)、濕延遲等參數(shù)的參考值，也根據(jù)該數(shù)據(jù)計(jì)算得到。ECMWF數(shù)據(jù)根據(jù)多源氣象信息綜合得到，在國(guó)際氣象學(xué)研究中具有十分權(quán)威的地位，常被學(xué)者用作評(píng)判對(duì)流層延遲精度的參考數(shù)據(jù)。②映射函數(shù)通常是將待求地區(qū)的位置和年積日代入GMF等全球經(jīng)驗(yàn)?zāi)Ｐ陀?jì)算得到(具體公式可參照文獻(xiàn)[11—12])，并忽略其誤差。此時(shí)觀測(cè)方程為

Li=[mfw(eli)maz(eli)cos(αi)maz(eli)sin(αi)]

(3)

式中，下標(biāo)i表示第i次觀測(cè)的量。

但事實(shí)上，映射函數(shù)并非準(zhǔn)確。以GMF映射函數(shù)模型為例，繪出西安地區(qū)(大地坐標(biāo)為34°10′41.74″N，108°59′8.16″E，504.427 m)2016年5°高度角時(shí)，GMF與映射函數(shù)參考值如圖1(此時(shí)暫不考慮方位角對(duì)映射函數(shù)的影響)。其中映射函數(shù)參考值采用ECMWF分層格網(wǎng)氣象數(shù)據(jù)并根據(jù)射線追蹤方法反算得到[13-14](以下試驗(yàn)的參考值計(jì)算方法與之一致)。

圖1 西安地區(qū)2016年5°高度角時(shí)干、濕分量映射函數(shù)Fig.1 Hydrostatic and wet mapping function of Xi’an area during 2016 at elevation of 5°

統(tǒng)計(jì)映射函數(shù)各分量誤差見表1。

表1西安地區(qū)5°高度角的映射函數(shù)各分量誤差

Tab.1ErrorsofeachmappingfunctioncomponentsofXi’anareaatelevationof5°

映射函數(shù)分量最小值最大值平均值標(biāo)準(zhǔn)差干分量-0.00730.03560.01980.0061濕分量-0.72211.87320.2800.3541

統(tǒng)計(jì)一年當(dāng)中誤差分布頻數(shù)直方圖見圖2。

圖2 映射函數(shù)誤差頻數(shù)直方圖Fig.2 Histogram of mapping function errors

從圖1、表1可知，高度角為5°時(shí)，西安地區(qū)干、濕分量映射函數(shù)均存在明顯的系統(tǒng)偏差，濕分量的RMS超過0.3，且最大誤差超過1.8；另外由圖2可知，當(dāng)高度角一定時(shí)，一年當(dāng)中的干、濕分量映射函數(shù)誤差大體呈正態(tài)分布。

2 顧及映射函數(shù)誤差的估計(jì)方法

顧及映射函數(shù)誤差的對(duì)流層延遲數(shù)學(xué)模型如下

maz(el)·[Gnscos(α)+Gewsin(α)]

(4)

Li=[maz(eli)cosαimaz(eli)sinαiZHD

(5)

此處需說明兩點(diǎn)：①為避免觀測(cè)方程中出現(xiàn)兩種不同待估參數(shù)相乘的現(xiàn)象，式(4)中引入濕延遲先驗(yàn)值，并在式(5)中忽略濕延遲改正與映射函數(shù)改正相乘的項(xiàng)；②觀測(cè)方程待估參數(shù)中的映射函數(shù)改正項(xiàng)，即δmfh(eli)、δmfw(eli)包含了隨高度角而變的量，需將高度角相關(guān)項(xiàng)剝離，方可按照最小二乘原理求解。故下文首先對(duì)δmfh(eli)和δmfw(eli)(也可將其視為映射函數(shù)誤差)的變化特性進(jìn)行分析。

2.1 映射函數(shù)誤差特性分析

在我國(guó)境內(nèi)選取西安、武夷、烏魯木齊、日喀則、?？?個(gè)地區(qū)，分別代表溫帶季風(fēng)、亞熱帶季風(fēng)、溫帶大陸、高原和熱帶季風(fēng)5種主要的氣候類型區(qū)域進(jìn)行研究。各地具體位置見表2。

表2 各地區(qū)大地坐標(biāo)

統(tǒng)計(jì)2016年各地區(qū)映射函數(shù)誤差的年均值和標(biāo)準(zhǔn)差隨高度角變化的情況見圖3。

從圖3可知，該誤差的平均值和標(biāo)準(zhǔn)差隨高度角增大而迅速減小，當(dāng)高度角大于40°時(shí)，映射函數(shù)誤差對(duì)斜路徑延遲的影響將降至亞毫米量級(jí)，基本可忽略不計(jì)；同時(shí)，該誤差呈現(xiàn)出較強(qiáng)的地域性，即隨地理位置不同而變，尤其是濕分量誤差平均值的變化最明顯。以下以西安地區(qū)為例進(jìn)一步研究。

統(tǒng)計(jì)西安2013—2016年各年映射函數(shù)誤差的平均值和標(biāo)準(zhǔn)差隨高度角的變化情況見圖4。

圖3 各地區(qū)干、濕分量映射函數(shù)誤差的平均值和標(biāo)準(zhǔn)差隨高度角的變化Fig.3 Variation of average and standard deviation of hydrostatic and wet mapping function errors with elevation in different areas

從圖4可知，相比地域性特點(diǎn)，映射函數(shù)誤差的時(shí)域性表現(xiàn)出較強(qiáng)的一致性，即該誤差的平均值和標(biāo)準(zhǔn)差各年縱向相比變化趨勢(shì)十分接近；同時(shí)觀察各曲線的形狀特征，可考慮利用指數(shù)函數(shù)對(duì)其擬合，具體形式如下，其中a1、a2、b1、b2為待擬合系數(shù)

y=a1exp(a2x)+b1exp(b2x)

(6)

統(tǒng)計(jì)各年份、各分量的平均值、標(biāo)準(zhǔn)差曲線的擬合系數(shù)見表3?？煽闯觯髂陻M合系數(shù)十分接近，與上文結(jié)論吻合；故考慮利用前3年擬合系數(shù)的平均值代表后續(xù)年份的情況，即對(duì)2013—2015年的a1、a2、b1、b2分別取平均，以其平均值代表2016年的系數(shù)，并統(tǒng)計(jì)這種方法的擬合精度，結(jié)果見表4。

表3 映射函數(shù)誤差擬合系數(shù)及殘差標(biāo)準(zhǔn)差統(tǒng)計(jì)

續(xù)表3

圖4 2013—2016年西安地區(qū)干、濕分量映射函數(shù)誤差的平均值和標(biāo)準(zhǔn)差隨高度角的變化Fig.4 Variation of average and standard deviation of hydrostatic and wet mapping function errors with elevation in Xi’an during 2013—2016

表4利用2013—2015年擬合系數(shù)平均值代替2016年系數(shù)后的擬合精度

Tab.4Statisticsoffittingaccuracywhenthecoefficientsof2016arechangedwith3-yearaverageonesfrom2013to2015

映射函數(shù)分量精度最小值最大值平均值標(biāo)準(zhǔn)差干分量平均值-7.1397e-48.6899e-5-9.8463e-51.9445e-4標(biāo)準(zhǔn)差-1.1015e-43.5432e-44.0411e-51.3045e-4濕分量平均值-0.0044 0.0017 -5.2254e-40.0014 標(biāo)準(zhǔn)差-0.00929.0738e-4-0.00240.0025

(7)

(8)

(9)

(10)

(11)

表5 統(tǒng)計(jì)參數(shù)各項(xiàng)系數(shù)

2.2 觀測(cè)方程改化

根據(jù)上文分析，可將觀測(cè)方程式(5)寫作如下形式

Li=[maz(eli)cos(αi)maz(eli)sin(αi)

(12)

此處需說明的是，在利用式(12)求解之前，應(yīng)首先根據(jù)式(3)，即傳統(tǒng)估計(jì)方法解算得到濕延遲先驗(yàn)值，也即ZWD0，同時(shí)得ZWD0的權(quán)逆陣；進(jìn)而采用部分參數(shù)加權(quán)平差[15-16]進(jìn)行解算。最終參數(shù)平差值的計(jì)算公式為

(13)

(14)

式中，V為殘差向量；n為觀測(cè)向量維數(shù)；t1為非隨機(jī)參數(shù)維數(shù)，此處為4，即除去具有先驗(yàn)信息的濕延遲后的待估參數(shù)個(gè)數(shù)。

此即為顧及映射函數(shù)誤差的兩步估計(jì)法，即第一步利用傳統(tǒng)估計(jì)方法(式(3))對(duì)濕延遲進(jìn)行估計(jì)；隨后根據(jù)部分參數(shù)加權(quán)平差，利用式(12)得到最終結(jié)果。為敘述方便，后文統(tǒng)一將其簡(jiǎn)稱為兩步法。

3 在新一代VLBI網(wǎng)觀測(cè)中的應(yīng)用效果仿真

新一代VLBI全球觀測(cè)系統(tǒng)(VLBI global observing system，VGOS)建成后，未來大地測(cè)量參數(shù)的解算精度有望達(dá)到1 mm量級(jí)[8,17]，將成為毫米級(jí)地球科學(xué)領(lǐng)域的重要技術(shù)基礎(chǔ)，然而大氣延遲誤差依舊會(huì)成為制約VGOS測(cè)量精度的重要因素[18-19]。故下文通過仿真，研究傳統(tǒng)方法與兩步法的估計(jì)效果。

根據(jù)未來VGOS網(wǎng)的指標(biāo)[20-21]，每臺(tái)站單天的觀測(cè)將達(dá)到上千次。故本文首先模擬生成1440個(gè)觀測(cè)方向(平均每分鐘1次)，即1440組方位角(0°～360°)、高度角(5°～90°)序列，各觀測(cè)方向具體在測(cè)站上空的分布及其在不同高度角、方位角區(qū)間的占比可參見圖5、圖6；隨后以西安為例，根據(jù)該地區(qū)2016年年積日185 d的分層格網(wǎng)氣象數(shù)據(jù)(數(shù)據(jù)來自ECMWF，原始分辨率為6 h，此處通過線性插值生成時(shí)間分辨率為1 min的數(shù)據(jù))，利用射線追蹤方法計(jì)算各方向斜路徑延遲值，并將其作為真實(shí)值；接著在該值中加入標(biāo)準(zhǔn)差為4 ps的白噪聲和方差為0.3 ps2/sec的隨機(jī)漫步噪聲(分別代表時(shí)延測(cè)量誤差和時(shí)鐘誤差[22]，該誤差量級(jí)與VGOS設(shè)備精度一致)，以此作為虛擬觀測(cè)量；最后分別采用傳統(tǒng)方法和兩步法每60歷元估計(jì)一次結(jié)果，統(tǒng)計(jì)兩種方法得到的斜路徑延遲殘差(見圖7，精度統(tǒng)計(jì)見表6)和濕延遲估計(jì)值(見圖8，精度統(tǒng)計(jì)見表7)。

表6 西安地區(qū)斜路徑延遲殘差

表7 西安地區(qū)濕延遲估計(jì)誤差

圖5 模擬觀測(cè)方向在測(cè)站上空的分布外圈為方位角軸，向徑為高度角軸Fig.5 Distribution of simulated observing orientations over the station outer circle represents the azimuthal axis, while inside diameter represents the elevation angle axis

圖6 模擬觀測(cè)方向在各高度角、方位角區(qū)間的占比Fig.6 Percentages of simulated observing orientations over various elevation angle and azimuth ranges

從圖7、表6可知，相比傳統(tǒng)法，兩步法的斜路徑延遲殘差被大大削減，RMS精度提升至6 mm以內(nèi)，改善幅度達(dá)46.22%，同時(shí)最大誤差降低54.25%；從圖8、表7知，濕延遲估計(jì)精度有小幅提升，但提升效果十分有限，改善程度在亞毫米級(jí)，這主要與第二次平差中濕延遲初值較多依賴第一次平差結(jié)果有關(guān)。

圖7 兩種方法計(jì)算的斜路徑延遲誤差Fig.7 Slant path delay errors of two different methods

圖8 濕延遲參考值與兩種方法的估計(jì)值Fig.8 Reference wet delay values and estimated ones by two methods

同理，將上文方法用于其他地區(qū)，得到的結(jié)果見表8、表9。可見該方法在不同地區(qū)均有相近的改善效果。

4 結(jié) 語

本文首先對(duì)對(duì)流層映射函數(shù)的誤差進(jìn)行研究，利用以往若干年的分層氣象格網(wǎng)數(shù)據(jù)構(gòu)建了映射函數(shù)誤差隨高度角變化的函數(shù)模型，并提出了顧及映射函數(shù)誤差的對(duì)流層延遲兩步估計(jì)法。試驗(yàn)表明：①本文構(gòu)建的映射函數(shù)誤差模型，能精確逼近研究地區(qū)映射函數(shù)誤差在各高度角的情況，可為分析各地區(qū)映射函數(shù)誤差提供依據(jù)；②兩步法可有效削減斜路徑延遲殘差，改善幅度達(dá)50%左右，這意味著采用該方法可將對(duì)流層延遲對(duì)VLBI觀測(cè)結(jié)果的影響削弱50%，同時(shí)對(duì)對(duì)流層濕延遲的估計(jì)結(jié)果也有一定程度的改善。另外，雖然目前VMF系列映射函數(shù)精度較高，但模型系數(shù)的發(fā)布往往滯后期較長(zhǎng)，而本文提出的方法可將映射函數(shù)改正值和濕延遲一并估計(jì)，有利于近實(shí)時(shí)性地獲取映射函數(shù)結(jié)果。然而本文方法尚存在不足之處，如構(gòu)建兩步法的觀測(cè)方程時(shí)，認(rèn)定濕延遲改正的量級(jí)足夠小，故忽略了該改正與映射函數(shù)改正相乘的項(xiàng)；但如果先驗(yàn)濕延遲精度不夠高(如大于1 cm)，則可能會(huì)引入毫米級(jí)甚至1 cm以上的誤差。因此如何選用更合理的濕延遲先驗(yàn)值是進(jìn)一步提高估計(jì)精度的關(guān)鍵。

表8 各地區(qū)斜路徑延遲殘差

表9 各地區(qū)濕延遲估計(jì)誤差