孫 帆，寧一鵬，邢建平*，王 森，代培培

(1.山東大學(xué)微電子學(xué)院，山東 濟(jì)南 250100；2.山東建筑大學(xué)測(cè)繪地理信息學(xué)院，山東 濟(jì)南 250101)

電離層位于地面上空，影響頻率信號(hào)傳播，是導(dǎo)航定位結(jié)果出現(xiàn)偏差的重要原因，制約著導(dǎo)航定位精度，對(duì)天基、地基系統(tǒng)的穩(wěn)定運(yùn)行及可靠性，航天器在電離層的運(yùn)行安全，無(wú)線電波的傳輸，以及全球?qū)Ш叫l(wèi)星系統(tǒng)具有很大影響，例如我國(guó)在電離層高發(fā)期多次觀測(cè)到電離層閃爍導(dǎo)致全球衛(wèi)星導(dǎo)航系統(tǒng)(Global NavigaTIon Satellite System，GNSS)信號(hào)失鎖[1]，鑒于此，研究精確度高、實(shí)時(shí)性優(yōu)良的電離層模型在民用、軍事等領(lǐng)域的應(yīng)用非常必要且具有重要的作用。

目前電離層模型研究進(jìn)展可分為三類。第一類為建立電離層延遲經(jīng)驗(yàn)?zāi)Ｐ?，如?guó)際GPS 服務(wù)(International GNSS Service，IGS)提供的全球電離層模型Klobuchar[2]、Abdus Salam 發(fā)布的NeQuick 模型[3]等經(jīng)驗(yàn)?zāi)Ｐ停诙悶橛?jì)算得到的數(shù)學(xué)擬合模型，如文獻(xiàn)[4-5]分別基于球諧函數(shù)和球冠諧函數(shù)擬合方法建立區(qū)域電離層(Total Electron Content，TEC)模型；周晨等[6]采用雙頻載波偽距平滑方法建立動(dòng)態(tài)電離層模型；第三類為利用神經(jīng)網(wǎng)絡(luò)、時(shí)間序列等方法得到預(yù)測(cè)模型，如陳春，吳振森等[7]提出通過(guò)神經(jīng)網(wǎng)絡(luò)ANN 方法提前1 h 預(yù)測(cè)電離層TEC 模型，謝劭峰等[8]利用IGS 中心提供的數(shù)據(jù)采用Holt-Winters’加法模型和乘法模型建立TEC 預(yù)報(bào)模型。

以上為傳統(tǒng)電離層建模方法，如文獻(xiàn)[9]依據(jù)最新發(fā)布數(shù)據(jù)建立格網(wǎng)點(diǎn)利用時(shí)間序列法得到所求位置點(diǎn)具有的電離層數(shù)據(jù)發(fā)布延遲誤差，誤差較大且不能反映電離層異常時(shí)的情況；文獻(xiàn)[10]利用偽距解算得到實(shí)時(shí)電離層延遲，受限于所用衛(wèi)星偽距數(shù)據(jù)質(zhì)量及數(shù)目，且高精度計(jì)算方法需考慮整周模糊度的消除，計(jì)算較為繁瑣。

本文將發(fā)布的格網(wǎng)數(shù)據(jù)進(jìn)行預(yù)測(cè)建模，結(jié)合多頻載波偽距法經(jīng)電離層異常條件判定改正，建立改正擬合模型。通過(guò)預(yù)測(cè)模型降低因數(shù)據(jù)質(zhì)量不良造成的偽距模型誤差，通過(guò)實(shí)時(shí)偽距模型提高預(yù)測(cè)模型與時(shí)間的關(guān)聯(lián)度，經(jīng)異常判定建立模型，有效提高模型實(shí)時(shí)性及模型擬合精度。本文對(duì)比分析建立模型與多頻載波偽距法、Holt-Winters’對(duì)比，驗(yàn)證模型的正確性及可行性。

1 電離層模型

1.1 單站區(qū)域電離層球諧級(jí)數(shù)模型

1.1.1 載波相位平滑偽距原理

基于GNSS 觀測(cè)數(shù)據(jù)的偽距觀測(cè)值建立電離層模型，針對(duì)偽距法建立TEC 模型精確度低的問(wèn)題，采用載波相位平滑偽距法可有效平滑偽距誤差，削弱觀測(cè)噪聲，提高偽距觀測(cè)值精度。GNSS 系統(tǒng)[11]雙頻偽距觀測(cè)方程和載波相位觀測(cè)方程如下:

式中:j表示衛(wèi)星編號(hào)，k表示頻率編號(hào)，ρ表示衛(wèi)星到觀測(cè)站接收機(jī)的偽距觀測(cè)值，L表示載波相位觀測(cè)值，c表示光速，dt1表示接收機(jī)鐘差，dt2表示衛(wèi)星鐘差，T與I分別表示信號(hào)傳播路徑對(duì)流層與電離層延遲誤差，br 和bs 分別表示接收機(jī)與衛(wèi)星硬件延遲誤差；ε表示偽距碼觀測(cè)噪聲，λ表示波長(zhǎng)，N表示載波相位整周模糊度[12]。

式中:f1和f2分別表示和Lj2 載波頻率，由于頻間差分碼偏差(Differential Code Bias，DCB)參數(shù)通常一天之內(nèi)幾乎無(wú)變化，代入IGS 發(fā)布的DCB 文件中觀測(cè)站接收機(jī)和衛(wèi)星的頻間DCB1、DCB2 參數(shù)值，得到信號(hào)路徑上的TEC 值(Slant Total Electron Content，STEC)[13]。

1.1.2 穿刺點(diǎn)經(jīng)緯度計(jì)算

式中:φPP、λPP表示衛(wèi)星到觀測(cè)站的穿刺點(diǎn)坐標(biāo)，φ、λ表示觀測(cè)點(diǎn)坐標(biāo)，H表示電離層高度值(一般情況為300 km～450 km，本文選用450 km)，RE表示地球平均半徑，為6 371 km，El、Az 分別表示衛(wèi)星相對(duì)于觀測(cè)站的高度角和方位角，z和z＇分別表示觀測(cè)站和穿刺點(diǎn)的天頂角[14]。

圖1 所示為穿刺點(diǎn)原理圖，可見(jiàn)衛(wèi)星穿刺點(diǎn)坐標(biāo)與地面坐標(biāo)正上空不是同一個(gè)點(diǎn)。

圖1 穿刺點(diǎn)原理圖

1.1.3 VTEC 計(jì)算

擬合兩個(gè)組合得到的測(cè)點(diǎn)上空電離層(Vertical Total Electron Content，VTEC)值，若雙頻時(shí)P1＝1，P2＝0，或者P1＝0，P2＝1，通過(guò)計(jì)算得到的穿刺點(diǎn)經(jīng)緯度坐標(biāo)，得到投影函數(shù)，從而計(jì)算得到VTEC[12]。

將STEC 通過(guò)單層模型投影函數(shù)轉(zhuǎn)化為觀測(cè)站上空垂直電離層VTEC，轉(zhuǎn)化公式為:

式中:F(z)表示穿刺點(diǎn)處的投影函數(shù)，z表示接收器天頂角，RE表示地球半徑，Hion表示薄層高度。

1.1.4 球諧函數(shù)模型

經(jīng)過(guò)式(5)計(jì)算后，將不同VTEC 值及其穿刺點(diǎn)代入到電離層球諧函數(shù)模型，表示如下:

式中:Cnm，Snm為球諧函數(shù)的參數(shù)，Pnm(sinβ)為締合Legendre 函數(shù)，nmax為模型階數(shù)，ms為過(guò)穿刺點(diǎn)的經(jīng)線與過(guò)地心和太陽(yáng)連線的經(jīng)線構(gòu)成的夾角，β為地磁緯度。

1.2 Holt-Winters’加法預(yù)測(cè)模型

選取時(shí)間序列法中的周期季節(jié)性時(shí)間序列法更符合電離層變化趨勢(shì)，提高電離層預(yù)測(cè)模型準(zhǔn)確度。使用歐洲定軌中心(CODE)發(fā)布的前10 天數(shù)據(jù)，預(yù)測(cè)得到當(dāng)天數(shù)據(jù)，構(gòu)建電離層預(yù)測(cè)模型。

加法模型初始值計(jì)算公式如下:

式中:Xt表示t時(shí)刻觀測(cè)值，It表示t時(shí)刻季節(jié)成分，bi表示t時(shí)刻趨勢(shì)成分，F(xiàn)t+m表示m時(shí)刻趨勢(shì)成分，m表示預(yù)測(cè)的時(shí)刻數(shù)，L表示季節(jié)長(zhǎng)度，α、β、γ表示平滑參數(shù)[9]。

1.3 擬合改正模型

1.3.1 初步擬合模型

基于雙頻載波相位平滑偽距法得到的球諧函數(shù)模型受數(shù)據(jù)質(zhì)量影響，精確度不穩(wěn)定，圍繞真實(shí)值跳動(dòng)較大，造成誤差較大問(wèn)題；依據(jù)CODE 發(fā)布數(shù)據(jù)進(jìn)行的時(shí)間序列法得到的預(yù)測(cè)值模型，因其實(shí)時(shí)性關(guān)聯(lián)不夠緊密，存在延遲問(wèn)題，造成預(yù)測(cè)模型雖然整體符合真實(shí)值變化，但當(dāng)電離層出現(xiàn)異常擾動(dòng)變化時(shí)，該模型無(wú)法及時(shí)預(yù)測(cè)，從而產(chǎn)生誤差。基于此，本文將二者進(jìn)行改正擬合，使其兩種方法優(yōu)勢(shì)互補(bǔ)，彌補(bǔ)傳統(tǒng)模型的不足。

通過(guò)加權(quán)最小二乘法對(duì)10 天內(nèi)3 萬(wàn)歷元數(shù)據(jù)處理，初步得到確定系數(shù)a，b。

1.3.2 電離層異常判定標(biāo)準(zhǔn)

初步擬合模型在精準(zhǔn)度上較兩種傳統(tǒng)模型略有提升，但在電離層出現(xiàn)異常擾動(dòng)情況時(shí)，受時(shí)間序列法原理制約，誤差仍然較大，拉低初步擬合模型整體精準(zhǔn)度，針對(duì)這一問(wèn)題，本文使用電離層異常判斷標(biāo)準(zhǔn)判定是否發(fā)生異常，為后面異常情況下改正擬合提供依據(jù)。

電離層發(fā)生異常時(shí)，增加使用載波相位平滑偽距法計(jì)算得到的球諧函數(shù)模型在擬合模型中占比，降低預(yù)測(cè)得到的時(shí)間序列模型占比。異常符合判斷時(shí)刻前后5 天的VTEC 值服從均值為μ標(biāo)準(zhǔn)差為σ的正態(tài)分布[15]，異常判定上下界表示如下:

C1iup表示判定上界，C2idown表示判定下界，下標(biāo)i表示時(shí)刻，VTECmedian表示預(yù)測(cè)值中值，VTEC(Projection)表示預(yù)測(cè)模型的VTEC。

2 算例分析

2.1 建模方法

考慮太陽(yáng)磁暴和地理位置，本文選取MIZU 參考站，時(shí)間為2013 年1 月10 日(BRDC0100)到20日(BRDC0200)，中心經(jīng)緯度范圍為10°×10°以內(nèi)，半徑＜600 km[9]的GPS 數(shù)據(jù)作為分析數(shù)據(jù)構(gòu)建模型。使用載波相位平滑偽距法提高偽距觀測(cè)精度，衛(wèi)星截止仰角設(shè)為15°，投影函數(shù)采用標(biāo)準(zhǔn)投影函數(shù)，時(shí)間采樣率為30 s，電離層薄層高度設(shè)置為450 km[16]。對(duì)比分析載波相位平滑偽距法得到的球諧函數(shù)模型[11]，利用CODE 發(fā)布數(shù)據(jù)構(gòu)建的周期季節(jié)性時(shí)間序列法Holt-Winters 加法預(yù)測(cè)模型，以及電離層異常改正擬合后的建立模型進(jìn)行區(qū)域VTEC 建模，以0.15°間隔劃分網(wǎng)格點(diǎn)，求解球諧級(jí)函數(shù)得到模型系數(shù)，反代入網(wǎng)格點(diǎn)，計(jì)算得出各網(wǎng)格點(diǎn)處的VTEC 值，預(yù)測(cè)值法使用CODE 發(fā)布的數(shù)據(jù)利用Holt-Winters 預(yù)測(cè)得到各網(wǎng)格點(diǎn)處的VTEC 值，基于以上兩種模型進(jìn)行最小二乘法擬合得到初步擬合模型，在初步擬合模型基礎(chǔ)上通過(guò)異常判斷改正異常點(diǎn)，構(gòu)建改正擬合模型。

2.2 結(jié)果與分析

2.2.1 初步加權(quán)擬合

通過(guò)最小二乘法如表1 所示得到未發(fā)生異常時(shí)刻時(shí)載波相位平滑偽距法與預(yù)測(cè)模型的加權(quán)系數(shù)，當(dāng)a＝0.642，b＝0.358 時(shí)，擬合模型最貼近真實(shí)值。

表1 時(shí)間序列法模型系數(shù)及其對(duì)應(yīng)誤差

2.2.2 異常判定

由于觀測(cè)站所處區(qū)域在1 月10 日至20 日均無(wú)太陽(yáng)磁暴現(xiàn)象，但17/18/20 均有小規(guī)模地震發(fā)生，考慮其對(duì)電離層擾動(dòng)影響，分別分析地震前后VTEC 值，取當(dāng)前時(shí)刻前后5 天的預(yù)測(cè)模型VTEC(或foF2)中值和標(biāo)準(zhǔn)差作為背景值，以2 倍標(biāo)準(zhǔn)差作為判斷標(biāo)準(zhǔn)，如果載波相位平滑偽距法的VTEC計(jì)算值超出背景值上下邊界，則認(rèn)為該時(shí)刻電離層出現(xiàn)異常，經(jīng)電離層異常閾值判定，增加載波相位平滑偽距法所占權(quán)重，降低預(yù)測(cè)模型所占權(quán)重，通過(guò)最小二乘法修正異常時(shí)刻系數(shù)a＝0.67，b＝0.33，異常修正后均方差(Mean Square Error，RMSE)提高0.5。

如圖2 所示，c1、c2為異常判定的上界和下界，y3計(jì)算值為載波相位平滑偽距法建立的球諧級(jí)函數(shù)模型得到的測(cè)站點(diǎn)上空VTEC 值，當(dāng)計(jì)算值超出上下閾值邊界時(shí)，判定電離層發(fā)生異常擾動(dòng)。對(duì)此時(shí)的加權(quán)擬合值進(jìn)行重新加權(quán)，增大計(jì)算值的所占比重。

圖2 計(jì)算值上下邊界異常判定

2.2.3 異常判定改正前后的模型與真實(shí)值對(duì)比

相比初步擬合模型，改正加權(quán)系數(shù)比重后，如圖3所示，在28-34、53-56、148-153、155-167 時(shí)刻精度均有提高，尤其在53-56 時(shí)刻提升最為明顯。

圖3 異常點(diǎn)處改正前后與真實(shí)值的對(duì)比

2.2.4 模型比較

如圖4 所示，y1代表VTEC 真實(shí)值，y2代表VTEC初步擬合值。預(yù)測(cè)模型在192 個(gè)時(shí)刻內(nèi)59 個(gè)時(shí)刻精度較高，相比于真實(shí)值，誤差低于0.2VTEC；32 個(gè)時(shí)刻預(yù)測(cè)值與真實(shí)值相差較大，誤差高于2VTEC，其他時(shí)刻計(jì)算值與真實(shí)值誤差均保持在(0.2，2)以內(nèi)，整體日RMSE 為22.441。由此可知，預(yù)測(cè)模型總體滿足2VTEC 精度，且精度較高的時(shí)刻較多，但某些時(shí)刻存在較大波動(dòng)，且當(dāng)遇到由于電離層擾動(dòng)造成的真實(shí)值突變時(shí)，預(yù)測(cè)模型往往無(wú)法有效預(yù)測(cè)，使得預(yù)測(cè)模型產(chǎn)生較大誤差，若能有效探測(cè)電離層是否發(fā)生異常，并進(jìn)行改正，可大大提高預(yù)測(cè)模型精度。

圖4 預(yù)測(cè)值與真實(shí)值對(duì)比

如圖5 所示，載波相位平滑偽距法得到的球諧函數(shù)模型在192 個(gè)時(shí)刻內(nèi)共39 個(gè)時(shí)刻精度較高，相比于真實(shí)值，誤差低于0.2VTEC；10 個(gè)時(shí)刻預(yù)測(cè)值與真實(shí)值相差較大，誤差高于2VTEC，其他時(shí)刻預(yù)測(cè)值與真實(shí)值誤差均保持在(0.2，2)以內(nèi)，日RMSE在14.618。由此可知，球諧函數(shù)模型，誤差較大時(shí)刻遠(yuǎn)遠(yuǎn)少于預(yù)測(cè)值模型，整體可以描述真實(shí)值變化情況。但受數(shù)據(jù)影響造，與真實(shí)值間的差值波動(dòng)較大，高精度時(shí)刻較少。

圖5 球諧模型值與真實(shí)值

如圖6 所示，初步擬合值模型在192 個(gè)時(shí)刻內(nèi)41 個(gè)時(shí)刻精度較高，相比于真實(shí)值，誤差低于0.2VTEC，3 個(gè)時(shí)刻精度較低，與真實(shí)值誤差高于2VTEC，其他時(shí)刻預(yù)測(cè)值與真實(shí)值誤差均保持在(0.2，2)以內(nèi)，整體日RMSE 為9.672。由此可知，初步擬合模型，中和了以上兩種模型的優(yōu)點(diǎn)，彌補(bǔ)了兩種模型的缺點(diǎn)，與預(yù)測(cè)值模型相比，低精度時(shí)刻變少，與球諧函數(shù)模型相比高精度時(shí)刻增多，并且總體擬合精度提高33.83%。

圖6 初步擬合模型與真實(shí)值對(duì)比

如圖7 所示，異常判斷改正后的加權(quán)擬合模型在192 個(gè)時(shí)刻內(nèi)61 個(gè)時(shí)刻精度較高，相比于真實(shí)值，誤差低于0.2VTEC，5 個(gè)時(shí)刻精度較低，相比于真實(shí)值，誤差高于2VTEC，其他時(shí)刻預(yù)測(cè)值與真實(shí)值誤差均保持在(0.2，2)以內(nèi)，日RMSE 在9.111。由此可見(jiàn)，改正后的擬合模型更好地將以上兩種模型優(yōu)點(diǎn)繼承下來(lái)，在保留高精度時(shí)刻不變甚至略有增加的情況下，低精度時(shí)刻略有降低，且總體擬合模型精度提高了37.67%。

圖7 異常判定后擬合模型與真實(shí)值對(duì)比

表2 顯示加權(quán)擬合后的模型在精度上優(yōu)于僅使用載波相位平滑偽距以及季節(jié)性時(shí)間序列法建立的球諧模型以及預(yù)測(cè)模型；隨著利用數(shù)據(jù)的增多，載波相位平滑偽距法和季節(jié)性時(shí)間序列法的日RMSE 值均在降低，而異常改正后的擬合模型在精度上更高，且在發(fā)生地震等影響電離層，使其發(fā)生異常擾動(dòng)時(shí)，該種方法精確度相比較為判斷改正的精度更高。

表2 10 日內(nèi)4 種方法與真實(shí)值的均方根誤差

測(cè)站及其周邊接收機(jī)對(duì)20 日24 h 內(nèi)上空經(jīng)度135°到145°，北緯35°到45°的電離層誤差累積渲染圖像如圖8～圖11 所示，初步加權(quán)擬合后的模型在精度上已經(jīng)在大部分區(qū)域優(yōu)于僅使用載波相位平滑偽距以及季節(jié)性時(shí)間序列法；而經(jīng)過(guò)異常改正后的擬合模型更是對(duì)初步擬合模型誤差較大區(qū)域進(jìn)行了明顯的改善。

圖8 預(yù)測(cè)值模型

圖9 球諧函數(shù)模型

圖10 初步擬合模型

圖11 異常判定后擬合模型

3 分析與總結(jié)

本文提出了一種電離層異常改正擬合模型，對(duì)比分析區(qū)域電離層建模的2 種傳統(tǒng)模型，選取10 天內(nèi)3 萬(wàn)個(gè)歷元數(shù)據(jù)作為分析數(shù)據(jù)，以0.15°為單位的網(wǎng)格點(diǎn)地圖進(jìn)行日誤差渲染，研究結(jié)果如下:

①載波相位平滑偽距法構(gòu)建球諧函數(shù)模型與預(yù)測(cè)模型經(jīng)異常判定條件后，建立的改正擬合模型，可有效提高預(yù)測(cè)模型實(shí)時(shí)性，減輕載波相位平滑偽距法構(gòu)建球諧函數(shù)模型的不穩(wěn)定性，提高電離層模型精度，支撐電離層建模研究。

②文中的試驗(yàn)是基于日本測(cè)站MIZU 的周邊區(qū)域10 天內(nèi)的數(shù)據(jù)來(lái)分析驗(yàn)證的，沒(méi)有對(duì)其他不同地區(qū)測(cè)站和更長(zhǎng)時(shí)間段的區(qū)域進(jìn)行分析。并且該種模型加權(quán)改正系數(shù)受地理分布位置、磁暴、地震和海嘯等電離層擾動(dòng)發(fā)生頻率規(guī)模差異和觀測(cè)數(shù)據(jù)質(zhì)量和密集程度等多種因素的影響，未來(lái)為了達(dá)到一個(gè)更為詳細(xì)的試驗(yàn)結(jié)論和更精細(xì)的模型，需要在時(shí)間長(zhǎng)度上、電離層異常擾動(dòng)的標(biāo)準(zhǔn)、不同測(cè)站地區(qū)上進(jìn)行更加完善的試驗(yàn)。