楊 玲，周春元，蘇小寧，李博峰

（1.同濟大學測繪與地理信息學院，上海 200092；2.蘭州交通大學測繪與地理信息學院，甘肅蘭州 730070）

電離層延遲是全球?qū)Ш叫l(wèi)星系統(tǒng)（Global Navigation Satellite Systems，GNSS）進行導航定位的一種重要誤差源，也是致使一般差分GNSS系統(tǒng)的定位精度隨用戶和基準站間的距離增加而迅速降低的主要原因［1-2］。由電離層引起的GNSS導航信號時延可達數(shù)米甚至百米級，會嚴重削弱導航定位的精度和可靠性。目前GNSS可以長期連續(xù)監(jiān)測電離層活動，基于雙頻觀測數(shù)據(jù)利用載波相位平滑偽距的方法進行電離層總電子含量（total electron content，TEC）的提取被證明是一種精度較高的手段［3-5］。基于提取的高精度電離層TEC觀測數(shù)據(jù)，利用球諧函數(shù)可建立全球電離層模型（Global Ionosphere Map，GIM），反映電離層TEC的時空分布特征和電離層活動規(guī)律，并為單頻GNSS接收機用戶提供電離層延遲改正，從而提高導航定位的精度［6-7］。

從1998年起，IGS（International GNSS Service）電離層工作組的各分支中心基于GNSS雙頻觀測數(shù)據(jù)開始獨立解算全球電離層產(chǎn)品，IGS提供的電離層產(chǎn)品為各分支中心產(chǎn)品的加權(quán)平均值［8］。由于在海洋和南極區(qū)域測站數(shù)目較少，在這些區(qū)域利用球諧函數(shù)進行全球電離層建模時存在異常值現(xiàn)象，已有相關(guān)研究利用不等式約束最小二乘來消除建模區(qū)域的TEC負值［9］。為彌補觀測數(shù)據(jù)的不足，利用前一日的GIM產(chǎn)品或經(jīng)驗?zāi)Ｐ驮跀?shù)據(jù)缺失區(qū)域添加虛擬觀測也被證明是有效的方法［10-11］。本文利用最新的國際參考電離層（International Reference Ionosphere，IRI）模型2016版本在觀測數(shù)據(jù)不足的區(qū)域添加虛擬觀測值進行約束，分析模型的改善效果，并利用單點定位結(jié)果評估添加IRI模型約束的有效性。同時，對于只存在C1、P2觀測值的測站，在C1、P1碼間偏差改正效果不佳的情況下，考慮將C1觀測值作為低精度的P1觀測值，增加全球的觀測數(shù)據(jù)，分析建模效果并在定位中評估數(shù)據(jù)選取策略的必要性。

1 全球電離層建?；驹?/h2>

1.1 求解STEC的方法

GNSS原始觀測量包含偽距和載波相位等，偽距觀測方程可表示為

式中：Pi為載波i（i=1，2）上的偽距觀測值；ρ0為衛(wèi)星和接收機之間的真實距離；dion，i為電離層延遲量；dtrop為對流層延遲量；、δr，t為衛(wèi)星和接收機鐘差；、d r，i分別為衛(wèi)星和接收機的碼延遲量，s和r分別為衛(wèi)星和接收機編號；c為光速；ε為隨機噪聲。

在進行電離層延遲和差分碼偏差（differential code bias，DCB）估計時，通常對P1、P2觀測值進行組合，得到組合的P4觀測值為

式中：dion，1、dion，2分別 為L1、L2載波上的電離層延遲；Ds、Dr分別為衛(wèi)星和接收機差分碼偏差。為獲取更精確的DCB估計，需利用載波L4（L4=L1?L2）觀測值對P4觀測值進行逐弧段的平滑處理。在進行載波相位平滑偽距之前，應(yīng)利用MW組合或電離層殘差組合檢測歷元中的周跳和粗差，并進行剔除［1］。

根據(jù)電磁波信號穿過電離層介質(zhì)時群速、相速的關(guān)系，GPS測距碼偽距和載波相位觀測值所受到的電離層延遲誤差大小相等、符號相反［12］。僅顧及f2項時，電磁波在電離層中傳播時所受到的電離層延遲可表示為

式中：fi表示載波i的頻率表示沿著信號傳播路徑s對電子密度Ne進行積分，并將其稱為傾斜路徑上的總電子含量I，I以TECU（total electron content unit）為單位，1 TECU表示每平方米含有1016個電子。將式（3）代入式（2）中，經(jīng)過變換可得到

1.2 球諧函數(shù)模型

電離層分布在距離地面60～2 000km高度范圍內(nèi)，信號傳播路徑的高度角和方位角各不相同，為將傾斜路徑的總電子含量轉(zhuǎn)化到天頂方向，利用單層模型將路徑上的所有自由電子都集中在F2層高度處。轉(zhuǎn)換關(guān)系可描述為［1-2］

式中：R為地球半徑，取6 371km；H為近似單層高度；α為常數(shù)，取0.978 2。

為適應(yīng)電離層的區(qū)域特征，通常采用球諧函數(shù)模型進行全球電離層建模。球諧函數(shù)模型可表示為［13］

式中：V(β，s)為地面上空某一高度的球殼上緯度為β、經(jīng)度為s處的垂直總電子含量；nmax為擬合時采用的最大階數(shù)；β為電離層穿刺點地心緯度；s=λ?λ0是日固坐標系中穿刺點經(jīng)度，λ、λ0分別為穿刺點的經(jīng)度和真太陽時；A nm、Bnm為全球電離層模型系數(shù)；是標準化的勒讓德多項式。

將式（4）、（5）和（6）進行結(jié)合并變換，可得到組合觀測方程為

式中：A nm、Bnm為待求電離層模型參數(shù)；Ds、Dr即為待求衛(wèi)星和接收機的差分碼偏差。由于組合觀測方程（7）是奇異的，通常添加所有GPS衛(wèi)星DCB和為零的基準約束，實現(xiàn)觀測方程的滿秩，從而分離衛(wèi)星和接收機的DCB。

1.3 IRI模型約束

IRI模型是根據(jù)大量電離層探空數(shù)據(jù)得到的電離層經(jīng)驗?zāi)Ｐ?，可計算海?0～2 000 km范圍內(nèi)的電子密度等參數(shù)，2 000 km以上至GPS衛(wèi)星軌道高度范圍內(nèi)的電子密度可根據(jù)IRI模型外推獲得［14-15］。2 000 km以下電子密度的方程式可描述為

式中：f(h)為高度h對應(yīng)的電子密度；yj為高度hj對應(yīng)的電子密度；w為待求系數(shù)。根據(jù)最小二乘解算系數(shù)w后可外推2 000 km高度以上的電子密度，累加50 km至GPS衛(wèi)星軌道高度范圍內(nèi)的電子密度即可得到總電子含量。

由于IGS測站在大西洋、南半球部分海洋地區(qū)分布較為稀疏，導致該地區(qū)無穿刺點覆蓋。在穿刺點空白區(qū)域逐歷元利用IRI模型計算虛擬穿刺點處的VTEC作為虛擬觀測值參與球諧函數(shù)建模，實現(xiàn)對全球電離層模型的約束。由于IRI模型計算的VTEC精度較低，通常降低虛擬觀測值的權(quán)重，而平滑后的P4計算所得VTEC權(quán)重默認為單位權(quán)［16］。虛擬觀測方程為

式中：βIPP、sIPP為所添加虛擬觀測點的經(jīng)緯度坐標。在進行模型約束時，通常需要經(jīng)過一輪的解算，求得衛(wèi)星和接收機的DCB參數(shù)，將式（7）轉(zhuǎn)化為式（6）。考慮到衛(wèi)星和接收機DCB短期內(nèi)的穩(wěn)定性，亦可利用前一日的衛(wèi)星和接收機DCB參數(shù)代入式（7），得到式（6）中純凈的VTEC觀測方程。最終聯(lián)立式（6）和式（9）進行解算，便可得到附加IRI模型約束的全球電離層模型。

2 實驗及分析

實驗選取2019年10月2日（年積日為275）全球420個測站的觀測數(shù)據(jù)，數(shù)據(jù)采樣率為30 s，其中205個測站包含P1、P2偽距雙頻觀測數(shù)據(jù)，其余215個測站P1觀測值缺失，但存在C1、P2雙頻觀測數(shù)據(jù)。實驗中利用GPS單系統(tǒng)進行建模，采用4種建模方案：①基于P1、P2觀測值的全球電離層建模，相關(guān)結(jié)果記為P1-P2；②在方案1的基礎(chǔ)上，在RMS值超限區(qū)域利用IRI模型進行約束，相關(guān)結(jié)果記為P1-P2-Limited；③將C1觀測值作為低精度數(shù)據(jù)取代缺失的P1觀測值，相關(guān)結(jié)果記為C1-P2；④在方案3的基礎(chǔ)上，在RMS超限區(qū)域利用IRI模型進行約束，相關(guān)結(jié)果記為C1-P2-Limited。

通過數(shù)據(jù)預(yù)處理，可獲得方案1和方案3測站上空的電離層穿刺點分布情況如圖1a和1b所示。可見，由于東太平洋和南半球海洋區(qū)域測站分布較少，其大氣上方無電離層穿刺點分布。經(jīng)過首輪解算，RMS值超限區(qū)域基本位于此范圍內(nèi)，因此方案2重點考慮在A、B、C、D 4個區(qū)域添加虛擬觀測值，如圖1a所示；方案4重點考慮在E、F、G 3個區(qū)域添加虛擬觀測值，如圖1b所示。

圖1 電離層穿刺點分布Fig.1 Distribution of IPP

采用如表1中的設(shè)置及處理策略進行解算，得到衛(wèi)星、接收機DCB和電離層參數(shù)，并利用電離層參數(shù)建立全球的電離層模型。

表1 設(shè)置及處理策略Tab.1 Strategies for setting and processing

2.1 模型解算精度分析

實驗中，方案2和方案4加入了虛擬觀測值，但方案1與方案2、方案3與方案4的DCB解算結(jié)果相同，僅在添加虛擬觀測值后電離層參數(shù)的求解中存在差異。將解算的DCB結(jié)果與IGS提供的DCB產(chǎn)品進行對比，評定方案1與方案3中模型參數(shù)的解算精度。統(tǒng)計2種方案下衛(wèi)星DCB的解算偏差及偏差分布如圖2a和2b所示，接收機DCB的解算偏差及偏差分布如圖3a和3b所示。

圖2 方案1、3衛(wèi)星DCB解算偏差對比Fig.2 Comparison of satellite DCB solution deviations of schemes 1 and 3

圖3 方案1、3接收機DCB解算偏差對比Fig.3 Comparison of receiver DCB solution deviation of schemes 1 and 3

可以看出，利用P1、P2觀測值時，由于精碼偽距具有較高的精度，雖然參與解算的測站數(shù)較少，但衛(wèi)星DCB的解算偏差均在0.1 ns之內(nèi)，90%的接收機DCB解算偏差在0.6 ns之內(nèi)；P1觀測值缺失時，將C1觀測值作為低精度的P1觀測值后，所有測站的觀測數(shù)據(jù)均可參與解算，90%的衛(wèi)星DCB解算偏差在1 ns之內(nèi)，90%的接收機DCB解算偏差在0.8 ns之內(nèi)。對比C1作為低精度P1觀測值前后，DCB解算精度有所下降，分析原因可能為：C1碼本身測距精度較低，且C1碼與P1碼之間存在碼間偏差，將C1作為低精度的P1觀測值后，雖然觀測數(shù)據(jù)增加，但由于添加的觀測值精度較低，導致參數(shù)解算精度整體略有下降，但大部分的DCB解算偏差仍在1 ns之內(nèi)，這也驗證了模型參數(shù)解算的有效性。

將每種方案解算的12組電離層參數(shù)利用球諧函數(shù)模型展開，加入歷元、經(jīng)緯度信息可得到一天內(nèi)任意時刻、任意經(jīng)緯度的電離層VTEC值。圖4分別統(tǒng)計了2019年10月2日0：00、6：00、12：00和18：00時刻不同方案的電離層建模結(jié)果與IGS產(chǎn)品的差值，在默認IGS提供的電離層產(chǎn)品精度較高的情況下，分析不同方案建模結(jié)果的優(yōu)劣。

圖4 附加IRI模型約束前后全球電離層VTEC對比（單位：TECU）Fig.4 Comparison of global ionospheric VTEC after addition of IRI model constraints(unit:TECU)

該圖顯示，方案1中，由于含有P1、P2觀測值的測站在海洋區(qū)域數(shù)量較少且分布不均勻，導致在太平洋東部、大西洋中部和印度洋南部海域建模結(jié)果存在異常值。方案2中添加虛擬觀測值后，6：00和12：00時刻在太平洋南部和印度洋南部海域的異常值范圍顯著減小，建模效果顯著提升；0：00時刻在印度洋西南部的異常值改善效果顯著，而太平洋大部分異常值改善效果甚微；18：00時刻2種方式的建模效果無顯著差異。結(jié)合對比其他時刻的建模結(jié)果表明，添加虛擬觀測值后，觀測數(shù)據(jù)在全球得到充分補充，建模效果在部分區(qū)域提升明顯，但較為稀疏的虛擬觀測數(shù)據(jù)對模型的約束力并不強，使得改善效果存在不確定性，并不能從根本上解決觀測數(shù)據(jù)缺失的問題。

方案3中將C1觀測值作為低精度的P1觀測值后，觀測數(shù)據(jù)明顯增多，與方案1相比建模異常值的范圍顯著縮小，全球整體建模效果提升顯著，但在東太平洋無測站的海域依然存在異常值現(xiàn)象。方案4中添加虛擬觀測值后，由于添加虛擬觀測的范圍有限且精度較低，與方案3相比各歷元的建模效果無顯著差異。

為客觀分析建模效果，統(tǒng)計了4種方案各歷元時刻電離層VTEC值與IGS參考值差值的平均值（Mean）、均方根（RMS）和標準差（STD）如圖5所示，4種方案所有歷元各統(tǒng)計量的均值如表2所示。

表2 不同方案全部歷元統(tǒng)計量均值統(tǒng)計Tab.2 Mean value of all epoch statistics in differ?ent schemes（單位：TECU）

圖5 附加IRI模型約束前后各統(tǒng)計量對比Fig.5 Comparison of statistics after addition of IRI model constraints

可以看到，與方案1相比，方案2中添加IRI模型約束后，8：00和20：00時刻的各統(tǒng)計量絕對值都有所減小，說明模型的內(nèi)外符合精度都有提升，模型改善效果顯著；其他歷元時刻的Mean的絕對值均有所增大，RMS值均有不同程度的減小。統(tǒng)計所有歷元時刻的Mean、RMS和STD平均值后發(fā)現(xiàn)，Mean平均值的絕對值由0.41增大到0.69，RMS和STD的平均值分別由4.22、4.14減小為2.99、2.85，模型內(nèi)符合精度降低而外符合精度有所提高，分析原因可能為添加虛擬觀測值后在觀測數(shù)據(jù)較多的區(qū)域模型精度略有降低，但在觀測數(shù)據(jù)缺失區(qū)域建模效果提升顯著。所有歷元附加約束前后STD值均有不同程度的減小，說明附加IRI模型約束可顯著降低模型中異常值的概率，改善了建模效果。

方案3中利用部分C1觀測值作為低精度P1觀測值后，與方案1相比各歷元Mean值的絕對值均有所增大，而RMS值除14：00時刻外均減小，說明加入額外的低精度觀測數(shù)據(jù)后，模型內(nèi)符合精度降低而外符合精度提高。方案4中添加IRI模型約束后，整體建模效果略好于方案3，但并無顯著差異?？傮w而言，利用IRI模型添加虛擬觀測值和利用C1觀測值作為低精度P1觀測值的目的都是為了增加觀測量，擴大電離層穿刺點的覆蓋范圍，提升全球整體建模的效果，但是由于添加的多余觀測數(shù)據(jù)精度往往更低，因此會降低電離層模型的內(nèi)符合精度。

2.2 定位精度分析

為分析不同方案電離層建模效果的優(yōu)劣，采用偽距單點定位的方式評定不同方案所得電離層產(chǎn)品的定位效果?？紤]到全球IGS測站的分布及赤道附近電離層雙駝峰結(jié)構(gòu)的影響，以經(jīng)度0°、緯度30°N和30°S為分界線將全球劃分為6大區(qū)域并標記為A、B、C、D、E和F，每一區(qū)域選取6個測站的觀測數(shù)據(jù)用于定位精度分析，區(qū)域劃分及參與測試的測站分布如圖6所示。

圖6 區(qū)域劃分及測試測站分布Fig.6 Area division and distribution of test stations

在定位中，分別加入4種方案所得到的電離層產(chǎn)品用于修正電離層延遲量，同時添加Klobuchar模型和IGS提供的GIM產(chǎn)品作為對比。為客觀地分析不同方案電離層產(chǎn)品的定位效果，統(tǒng)計了各區(qū)域6個測站定位偏差的平均RMS分布如圖7所示，各方向的平均定位誤差（95%顯著水平）如圖8所示?？梢钥闯?，4種方案電離層產(chǎn)品的定位效果普遍優(yōu)于Klobuchar模型，在觀測數(shù)據(jù)較多的A、B、D和F區(qū)域，利用不同方案電離層產(chǎn)品的定位效果與IGS提供的GIM產(chǎn)品效果相當，而在缺少觀測數(shù)據(jù)的C和E區(qū)域，不同方案電離層產(chǎn)品的定位效果表現(xiàn)出較大的差異。

圖7 不同方案各區(qū)域平均RMS分布Fig.7 Average RMS distribution of different schemes in each area

可以發(fā)現(xiàn)，方案2中添加虛擬觀測值之后，與方案1相比在測站數(shù)據(jù)較多的A、B、D和F區(qū)域各方向的RMS整體有略微增大的趨勢，N、E、U方向RMS值增大的平均值為0.015 m、0.001 m和0.002 m。而在測站數(shù)據(jù)較少的C和E區(qū)域，各方向上的RMS顯著減小，其中以E區(qū)域N、E方向最為明顯，RMS分別減小0.109 m、0.073 m，U方向減小0.045 m。圖8顯示，添加虛擬觀測值后，在A、B、D和F區(qū)域各方向定位誤差整體略微增大，N、E、U方向定位誤差平均增加0.016 m、0.003 m和0.012 m，在觀測數(shù)據(jù)較少的C區(qū)域定位誤差略微減小，而觀測數(shù)據(jù)最少的E區(qū)域定位誤差減小更為顯著，N、E、U方向上分別減小0.260 m、0.146 m和0.103 m。說明在加入低精度虛擬觀測值之后，對于觀測數(shù)據(jù)較多的區(qū)域，重新解算電離層模型時參數(shù)精度有所下降，影響了該區(qū)域電離層模型的精確建立，最終導致定位效果略有下降。而在觀測數(shù)據(jù)不足的南半球海洋區(qū)域，添加低精度虛擬觀測值一定程度上彌補了觀測數(shù)據(jù)不足的問題，改善了電離層參數(shù)中異常值的出現(xiàn)，增加了該區(qū)域電離層模型的整體符合度，并一定程度上提升了該地區(qū)的定位效果。

圖8 不同方案各區(qū)域平均定位誤差（95%顯著水平）Fig.8 Average positioning error of different schemes in each area(95%significant level)

方案3中將C1觀測值作為低精度P1觀測值之后觀測數(shù)據(jù)顯著增多，與方案1相比，在觀測數(shù)據(jù)較多的A、B、D和F區(qū)域RMS和定位誤差有略微增大趨勢，而在觀測數(shù)據(jù)不足的C和E區(qū)域RMS和定位誤差的增大趨勢則更為顯著，其中E區(qū)域的RMS在N、E、U方向上分別增加0.200 m、0.024 m和0.384 m，定位誤差在N、E、U方向上分別增加0.524 m、0.111 m和1.549 m，方案3的定位效果整體不如方案1，說明利用C1觀測值作為低精度的P1觀測值后，雖然圖4a和4c直觀顯示全球整體建模效果有所提升，但是一定程度上損失了模型的精度，導致在測試區(qū)域采用方案3的電離層模型的定位效果并不如方案1。

以上分析表明，4種方案中方案2整體上最優(yōu)，與方案1相比，雖然在測站數(shù)目較多區(qū)域的定位效果有所下降，但下降幅度僅在厘米級，而在缺少觀測數(shù)據(jù)的區(qū)域，定位效果的提升能達到分米級，且在嚴重缺少觀測數(shù)據(jù)的E區(qū)域最為明顯。因此，在利用載波相位平滑偽距的方法進行全球電離層建模時，選取利用包含P1、P2雙頻觀測值的觀測數(shù)據(jù)是非常有必要的，在觀測數(shù)據(jù)嚴重不足的區(qū)域，利用IRI模型添加虛擬觀測值能在保證定位精度的基礎(chǔ)上提升模型整體的建模效果，而存在C1、P2觀測值時，則不建議將C1觀測值作為低精度的P1觀測值參與電離層模型的建立。

3 結(jié)論

采用載波相位平滑偽距的全球電離層建模方法，分析對比4種電離層建模方案的效果差異，評定了不同方案所得到的電離層產(chǎn)品在單頻偽距單點定位中的改善效果，最終確定了不同數(shù)據(jù)選取策略和添加IRI模型約束的必要性和有效性。

基于P1、P2原始雙頻觀測數(shù)據(jù)添加IRI模型約束后，部分歷元電離層建模整體效果改善顯著，異常值出現(xiàn)范圍顯著縮小，但是該方法降低了模型的內(nèi)符合精度，且無法從根本上解決觀測數(shù)據(jù)缺失的問題。在定位效果上，添加IRI模型約束后，在觀測數(shù)據(jù)不足的區(qū)域定位誤差顯著減小，N、E、U方向上分別減小0.260 m、0.146 m和0.103 m。因此在進行全球電離層建模時，利用IRI模型在觀測數(shù)據(jù)不足區(qū)域進行約束是非常有效的。

將C1觀測值作為低精度P1觀測值后，電離層整體建模效果明顯提升，但模型精度損失嚴重，數(shù)據(jù)不足區(qū)域的定位誤差在N、E、U方向上分別增加0.524 m、0.111 m和1.549 m。因此進行全球電離層建模時，利用包含P1、P2雙頻觀測值的觀測數(shù)據(jù)是非常有必要的，而存在C1、P2觀測值時，則不建議將C1觀測值作為低精度的P1觀測值參與電離層模型的建立。

作者貢獻聲明：

楊 玲：統(tǒng)籌論文的研究工作，指導論文的研究方向并修改論文。

周春元：查找文獻、分析數(shù)據(jù)，負責論文的撰寫。

蘇小寧：提供部分測站數(shù)據(jù)及相關(guān)問題的咨詢。

李博峰：參與部分數(shù)據(jù)分析，提供修改意見。