綜合半參數核估計和自回歸補償的全球電離層總電子含量預報模型

邱封欽， 潘雄， 羅小敏， 賴祖龍， 蔣可， 王廣興

中國地質大學(武漢)地理與信息工程學院， 武漢 430078

0 引言

電離層延遲是衛星導航定位系統的一項重要誤差來源，電離層對衛星偽距和載波信號產生時延效應，最大可達數十米，在實時單頻導航定位中必須對電離層延遲加以改正.電離層總電子含量(Total Electron Content ,TEC)直接影響電離層延遲的大小，為了降低電離層對于衛星信號的影響，需要對電離層TEC進行精確預報(張強等，2014；章紅平等，2006).

電離層預報模型分為經驗模型和非經驗模型.經驗模型有Klobuchar、IRI、Bent等模型(袁運斌等，2017)；非經驗模型有球諧函數模型、多項式模型等(李志剛等，2007；安家春等，2015；蔡成輝等，2015).這些模型計算簡單，使用方便，但電離層受到太陽活動、災害性天氣等很多因素影響，模型無法精確反映電離層不規則擾動，且模型本身具有不合理性，例如Klobuchar模型將夜間TEC視為常數，導致經驗模型預報精度不高(Pongracic et al.,2019；周仁宇等，2019).球諧函數模型通過預報電離層球諧函數系數，將預報得到的球諧系數代入球諧函數模型中，計算出格網點的TEC預報值(Li et al.,2019；Schaer, 1999).該模型可計算全球任意經緯度的電離層TEC值，但是模型觀測數據存在系統誤差以及模型誤差的影響.對此，王喜江等(2020)提出了一種基于半參數核估計的電離層球諧函數模型(Semiparametric-Spherical Harmonic, Semi-SH)，將系統誤差及小周期項歸入到非參數分量部分，有效地削弱了系統誤差對預報精度的影響.

在半參數核估計理論中，窗寬參數是一個非常關鍵的平滑因子，窗寬參數的選取直接影響參數估計和非參數估計的精度，作為權衡核估計偏差和方差之間的參數，它不能使參數估計偏差和非參數估計偏差同時變小(趙春茹，2011；李文濤等，2020；Yan et al.,2020).顏雄(2019)在利用半參數核估計模型預報北斗鐘差時，分析了不同核函數和窗寬參數選取對模型誤差的影響，發現選取不同核函數及窗寬參數建模對鐘差數據的有用信息均存在提取不充分的情況.所以半參數球諧函數模型對電離層TEC預報依然存在偏差.鑒于此，本文利用半參數球諧函數模型對15階球諧系數預報，擬合殘差建立自回歸(Autoregressive model, AR)模型進行補償預報，提出一種組合預報模型，減小因窗寬參數選取導致的估計偏差，提高電離層預報精度.

1 模型構建

1.1 半參數球諧函數模型

電離層球諧函數各階次球諧系數隨時間具有周期性變化特性，利用頻譜分析法對各階次球諧系數時間序列三角級數分解，建立球諧函數頻譜分析模型，即可對球諧系數擬合和預報，進而預報全球TEC值.

電離層球諧系數時間序列可表示為

L(t)=F(t)+U(t)+S(t)+Δ，

(1)

其中,L(t)為球諧系數時間序列；F(t)為趨勢項，可表示為式(2)的多項式回歸方程；U(t)為周期項，可表示為式(3)的三角級數形式；S(t)為模型的隨機信號，與球諧系數無確定關系，是非參數分量；Δ為觀測噪聲.

F(t)=a0+a1t+a2t2,

(2)

(3)

(4)

其中

V=CX+S-(L-F).

(5)

(6)

(7)

MSE(h)=tr((H(h)C)-1)H(h)(VVT+SST)

2.7 保持樹勢中庸 果園平衡施肥，減少氮肥施入，減少激素噴施，調節葉果比，前期適度低溫管理等措施均可保持樹體穩健中庸，協調生長，從而有效防止棗果萎蔫。

·H(h)T(H(h)C)-1),

(8)

式(8)中，tr(·)為矩陣的跡；H(h)為影響矩陣，H(h)=CT(I-W)T(I-W)；當MSE(h)取值最小時，確定對應的窗寬參數.

結合式(5)和式(6)，利用最小二乘原則可求得待估參數X和非參數分量S的估值：

(9)

(10)

最后將擬合的趨勢項、解算的待估參數和隨機信號代入式(1)，建立半參數球諧函數模型，實現對球諧系數的擬合和預報.

1.2 綜合半參數核估計與自回歸補償預報模型

由式(8)、式(9)和式(10)可知，窗寬參數的合理選擇對半參數球諧函數模型的預報結果影響較大.半參數球諧函數模型中，窗寬參數作為一個懲罰因子，在分離參數分量和非參數時，會出現估計偏差，需對擬合殘差中的有用信號進一步分離.將球諧系數觀測值減去半參數球諧函數模型擬合值，得到的球諧系數殘差序列Y采用AR模型進行擬合和預報.

(11)

Yi=b0+b1Yi-1+b2Yi-2+…bnYi-n+εt,

(12)

其中，b0為常數項，b1,b2,…,bn為自回歸模型系數；i=1,2,…,N；N為殘差序列長度；n為AR模型階數；εt是均值為0，方差為σ2的白噪聲.令

Dn=[1Yi-1Yi-2…Yi-n],

Xbn=[b0b1b2…bn]T,

將式(12)轉化為矩陣形式：

Yi=DnXbn.

(13)

i時刻的殘差值是由該時刻之前的n個時刻的殘差值進行估計，將殘差序列代入式(13)組成觀測方程,根據式(14)最小二乘原則可求得AR模型參數，建立殘差AR模型，實現對半參數球諧函數模型殘差的補償預報.

Xb=(DTD)-1DTY.

(14)

改進的組合模型預測值采用串聯式組合，如式(15)所示：

R(ti)=L(ti)+Y(ti),

(15)

式中，R(ti)為組合預報模型預測值；L(ti)和Y(ti)分別為半參數球諧函數預測值和AR模型殘差預測值.

將半參數球諧函數模型預報值與AR模型殘差預報值求和，即得到組合預報模型的球諧系數預報值.將各階次球諧系數預報值代入全球電離層電子含量球諧函數模型，即可求得全球電離層TEC格網值.模型方程如式(16)所示：

(16)

2 算例分析

2.1 數據來源

利用歐洲定軌中心(Center for Orbit Determination in Europe, CODE)發布的高精度全球電離層球諧系數數據作為樣本序列，選取2018年1月1日至2018年2月4日共計35天的15×15階共計256個球諧函數系數作為觀測值，采樣頻率為1 h.為了驗證組合模型的有效性，分別利用Semi-SH模型和組合模型進行預報，將預報結果與CODE所發布的高精度全球電離層TEC格網值(CODG值)進行對比.CODG值是經差為5°、緯差為2.5°(87.5°N—87.5°S)共計5183個高精度的全球格網值.本文采用的精度評定指標為均方根誤差(RMS)和平均相對精度(P)

(17)

(18)

2.2 單天預報

利用2018年1月1日—1月30日共計30天的球諧系數，預報2018年1月31日1天內逐小時的球諧函數系數，根據球諧函數模型計算預報的總電子含量，實現全球TEC單天預報.

為了檢驗預報精度，將CODE發布的一天預報產品(C1PG 值)、Semi-SH 模型預報值和組合模型預報值與CODG值24 h的全球格網值進行比較，其差值圖分布如圖 1(a、b、c)所示.

圖1可以看出，Semi-SH模型和組合模型相比C1PG模型精度有明顯提升，在0～16 h精度提升較大，兩種模型的預報精度均隨預報時間長度降低.在南極地區，C1PG模型預報精度較低.組合模型預報殘差大部分在±3 TECU以內，說明了組合模型的有效性.從圖1b和1c可以看出，組合模型在前半天對單一的Semi-SH模型的改正效果較顯著，主要表現在赤道附近殘差較大區域；在后半天改正較小，這與自回歸模型相關性逐漸減弱，模型趨穩有關.預報精度較低地區主要分布在南北緯20°左右，在赤道兩側呈駝峰狀，與電離層赤道異常現象相關(田耀宇等，2019；黃江等，2013；Balan et al.,2018).

圖1 C1PG值(a)、Semi-SH模型(b)和組合模型(c)單天預報殘差全球分布圖Fig.1 The global map of residuals of C1PG (a), Semi-SH (b) and combined model (c) value with respect to CODG value

表1是C1PG模型、Semi-SH模型、組合模型三種模型的預報殘差分析.由表1可知，組合模型預報精度高于單一的Semi-SH模型，預報均方根誤差為1.59 TECU，C1PG模型精度最低；在北半球地區，三種模型的預報精度都比南半球地區的預報精度高，這可能是南半球處于夏季，電離層電子活躍度高于北半球，以及南半球海洋面積廣闊，地基觀測站觀測數據密度小，導致預報精度降低.

表1 C1PG模型、Semi-SH模型和組合模型預報值殘差分析Table 1 C1PG model, Semi-SH model and combined model prediction residual analysis

組合模型預報差值在1 TECU以內占比54.30%、在3 TECU以內占比93.40%，普遍高于C1PG模型和單一的Semi-SH模型；三種模型預報差值在5 TECU以內占比均達到98%以上，其中組合模型預報差值在5 TECU以內占比98.61%，優于CIPG模型和單一的Semi-SH模型.

2.3 多天預報

利用2018年1月1日—1月30日共計30天的球諧系數，單次預報一天，滑動預報2018年1月31日—2月4日共計5天的球諧系數，根據球諧函數模型計算預報的總電子含量，實現全球TEC預報，預報結果與COGD值進行檢驗.

圖2為Semi-SH模型和組合模型滑動5天預報值與CODG值對比的均方根誤差統計.因C1PG為CODE發布的一天預報產品模型，與滑動5天預報結果可比性不強，故未引入比較.從圖中可以看出，組合模型預報精度均高于Semi-SH模型.隨著預報時長增加，兩種模型的預報精度都逐步降低.

圖2 Semi-SH模型和組合模型預報RMS統計Fig.2 RMS statistics of Semi-SH value and combined model value relative to CODG value

取滑動預報5天每天00∶00 h共計25915(71×73×5)個全球格網點預報殘差數據進行統計，圖3a和3b分別是Semi-SH模型和組合模型滑動5天預報殘差概率密度統計圖.從圖3中可以看出，兩種模型預報殘差主要分布在±3 TECU以內，大于5 TECU的殘差值所占比例較小，說明兩種模型用于全球電離層TEC預報的可靠性.組合模型預報殘差值小于1 TECU占比53.43%， Semi-SH模型小于1 TECU占比49.83%，組合模型有明顯提高.

圖3 Semi-SH模型(a)和組合模型(b)預報殘差密度統計Fig.3 Semi-SH model (a) and combined model (b) forecast residual probability density statistics

為了分析不同緯度下組合模型的預報精度，根據滑動預報5天的總電子含量數據，分別選取不同緯度TEC值進行分析.圖4a—4e分別為Semi-SH模型和組合模型滑動預報5天00∶00 h的總電子含量與CODG格網值比較，從上至下依次為70°N、50°N、15°N、15°S、50°S和70°S六個緯度數據.

從圖4可以看出Semi-SH模型和組合模型預報值都能較好反映出真實總電子含量及其變化趨勢，組合模型的改進效果隨著預報時間的增長而下降，這與預報誤差滑動積累有關.白天中低緯度預報殘差較大，這與白天中低緯度地區電離層電子活躍度高有關.在中低緯度地區，第2天的CODE-TEC值與其他天有一定區別，導致該天Semi-SH模型和組合模型的預測誤差較大，這可能是當天該時刻地磁環境出現擾動造成的.組合模型在中高緯度地區以及預報前期比Semi-SH模型更加接近CODG值，后期與Semi-SH模型預報結果相近，改正效果減弱.

圖4 Semi-SH模型、組合模型預報值與CODG值對比Fig.4 Comparison of Semi-SH model, combined model forecast value with CODG value

受季節因素影響，北半球高緯度地區夜間電離層電子濃度較低，該區域CODG值含有較多零值，相對精度失去評定意義，故只統計南半球15°S、50°S和70°S三個緯圈兩種模型滑動預報5天的預報精度，結果見表2.從表中可以看出，組合模型滑動預報5天精度在低、中、高緯度下均高于Semi-SH模型.組合模型在高緯度地區的改正效果最好，相對精度提高了2.32%；其次是低緯度地區，相對精度提高了2.14%；在中緯度地區改正效果較小，相對精度提高了1.32%.兩種模型預報結果的均方根誤差大致隨著緯度的提高而減小.

表2 Semi-SH模型和組合模型預報精度對比Table 2 Comparison between the prediction accuracies of the Semi-SH model and the combined model

3 結論

本文對發布的15階電離層球諧系數進行建模與預報，利用自回歸模型對殘差進行建模，實現對窗寬參數選取產生的估計偏差的補償修正，提出一種全球電離層TEC組合預報模型，并通過實驗，得出以下結論：

1)組合預報模型在2018年1月31日當天的預報均方根誤差為1.59 TECU，預報差值小于3 TECU的格網點占90%以上，表明組合預報模型的可行性；

2)與Semi-SH模型相比，組合模型單天預報差值小于1 TECU、3 TECU和5 TECU的格網點分別占54.30%、93.40%和98.61%，均高于Semi-SH模型，且組合預報模型的預報精度優于單一Semi-SH模型；

3)在滑動預報5天中，組合模型預報精度優于Semi-SH模型，組合模型對高緯度和低緯度地區改正效果明顯，在中緯度地區改正較小.

致謝感謝CODE提供的電離層數據產品.