顧及映射函數(shù)誤差的對流層延遲兩步估計法

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

1. 信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450052； 2. 西安測繪研究所，陜西 西安 710054； 3. 地理信息工程國家重點實驗室，陜西 西安 710054

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

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

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

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

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

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

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

(1)

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

(2)

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

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

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

(3)

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

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

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

統(tǒng)計映射函數(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)計一年當(dāng)中誤差分布頻數(shù)直方圖見圖2。

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

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

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

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

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

(4)

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

(5)

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

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

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

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

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

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

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

y=a1exp(a2x)+b1exp(b2x)

(6)

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

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

續(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)計參數(shù)各項系數(shù)

2.2 觀測方程改化

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

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

(12)

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

(13)

(14)

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

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

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

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

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

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

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

圖5 模擬觀測方向在測站上空的分布外圈為方位角軸，向徑為高度角軸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 模擬觀測方向在各高度角、方位角區(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%，同時最大誤差降低54.25%；從圖8、表7知，濕延遲估計精度有小幅提升，但提升效果十分有限，改善程度在亞毫米級，這主要與第二次平差中濕延遲初值較多依賴第一次平差結(jié)果有關(guān)。

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

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

同理，將上文方法用于其他地區(qū)，得到的結(jié)果見表8、表9?？梢娫摲椒ㄔ诓煌貐^(qū)均有相近的改善效果。

4 結(jié) 語

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

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

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