基于穩健卡爾曼濾波的傾斜探測電離層MUF短期預報方法

李國軍，鄭廣發，葉昌榮，周銀萍

（重慶郵電大學超視距可信信息傳輸研究所，重慶 400065）

1 引言

電離層是隨時間、空間高度變化的復雜系統，它對短波天波的傳播、地空和星間的信息鏈路有著至關重要的影響，是短波通信鏈路的中間介質。由于受太陽輻射、大氣環流和地磁環境等多種因素的影響，電離層具有時變色散特性，容易發生擾動，影響短波通信質量，甚至導致通信中斷。因此準確預報電離層能反射的最大可用頻率（MUF,maximum usable frequency）[1]，在短波通信[2-3]以及電離層探測[4]中具有重要的意義。

電離層的頻率預報主要分為2 種，即長期預報[5]和短期預報[6]。長期預報是通過建立電離層參考模型（IRI,international reference ionosphere）[7]，根據太陽活動指數以及電離層月中值參數，預測未來時間電離層的變化情況。正是由于這種月中值特性，長期預報方法在具體時刻的預報存在較大的偏差，往往僅適用于頻率的初選[8]。相比于長期預報，短期預報根據電離層短期內會保持相對穩定的電離特性，結合實測數據，通過算法得出短期電離層的變化情況，預報準確性更高，適用性更強。因此，在過去的幾十年里，電離層預報的發展趨勢也逐步由長期預報向短期預報發展[9]。目前，電離層短期預報方法主要有等效太陽黑子數法[10]、自相關分析法[11]、多元線性回歸法[12]、神經網絡法[13-14]、卡爾曼濾波（KF,Kalman filter）法[15-16]、相似日法[17]等。其中，等效太陽黑子數法和神經網絡法需要大量的歷史數據作為輸入，而自相關分析法和卡爾曼濾波法需要的歷史數據相對較少，因此更常用。文獻[18]對比了卡爾曼濾波法和自相關分析法在電離層短期預報的精度，得出卡爾曼濾波法短期預報的性能總體上優于自相關分析法。但對于卡爾曼濾波法，當觀測噪聲為非高斯干擾噪聲或者存在大量觀測野值的情況下，濾波精度和穩定性會下降，甚至可能會導致不收斂。本文旨在提出一種強跟蹤的穩健卡爾曼濾波（RKF,robust Kalman filter）頻率預報方法，改善傳統電離層頻率預報的準確性。

2 卡爾曼濾波

卡爾曼濾波是1960 年由Kalman 提出的，其基本思想是以分析誤差的最小方差為標準，假設系統為線性系統，噪聲為高斯白噪聲，采用信號與噪聲的狀態空間模型，根據預測方程和量測方程對系統狀態進行評估與預測，是一種最優化自回歸數據處理方法。具體地，利用前一時刻的估計值和當前時刻的觀測值來更新對狀態變量的估計，得到最佳的狀態估計值。由于卡爾曼濾波的基本方程在時間域中是遞推形式的，遞推過程是一個“預測?修正?再預測”的過程，在求解時不需要大量的已知數據，一旦觀測到新的數據，就可以根據算法得到新的濾波值，保持結果最優，因此非常便于數據的實時處理和計算機實現，得到廣泛的應用[19]。

2.1 傳統卡爾曼濾波預測模型

考慮到電離層MUF 的時間離散性，系統的預測過程可看作簡單的時間離散、線性隨機的動態系統，模型框架如圖1 所示，其中，IZ 表示正向推導（濾波），IZ?1表示反向推導的上一時刻的狀態。

圖1 模型框架

模型表達式為[20]

式(1)是描述隨機動態系統狀態隨時間變化的狀態方程，式(2)是描述狀態矢量與觀測矢量之間關系的量測方程。通過狀態方程，可以利用上一時刻（t時刻）的狀態隨機變量Xt估計出下一時刻（t+1時刻）的狀態隨機變量Xt+1，估計過程受到隨機噪動wt的影響，這種隨機噪動也稱為動態噪聲，其均值為0、方差為Q，服從wt≈N(0,Q)。對于量測方程，狀態變量與量測變量之間也受到量測噪聲的影響，其均值為0、方差為R，服從vt≈N(0,R)。

式(1)和式(2)所組成的系統模型，可用于預測電離層的最大可用頻率，進而得出傳統的卡爾曼濾波方法，遞推預測過程如下。

協方差矩陣量測更新為

KF 狀態和遞推轉移流程如圖2 所示。

圖2 KF 狀態和遞推轉移流程

2.2 Huber-M 穩健卡爾曼濾波方法

在傳統卡爾曼濾波中，w t和vt為互不相關的高斯白噪聲，實現的功能也是最小均方根誤差意義下的最優估計。在標準模型里，量測方程中噪聲是一種理想的情況，而在實際的量測中，測量值可能會存在野值，一旦量測中含有大量的野值，標準卡爾曼濾波的性能就會大幅退化，甚至是不收斂的。因此本文提出一種基于Huber-M 估計的穩健式卡爾曼濾波方法，來消除量測野值對預報精度的影響。

考慮到實際測量中包含野值的情況，本文改進量測噪聲的分布函數，具體描述為

其中，FN為量測噪聲的主體分布，保持和原模型一致的高斯分布；FP為污染部分，它可以是方差遠大于主體分布的其他厚尾分布；α為污染率，其值越大，代表在實際的測量過程中存在野值的數量就越多。

Huber 最先提出M 估計[21]，其核心是使定義的代價函數取得最小值。本文將該方法應用于傳統卡爾曼濾波的觀測更新方程中，并對其進行改進。推導過程如下。

首先，定義了t+1 時刻的狀態誤差為

其中，xt+1為t+1 時刻的狀態真實值，為t時刻對t+1 時刻狀態的先驗預測值。

聯合式(2)構造線性回歸模型為

定義矩陣

則式(10)可以改寫為

定義Huber-M 估計的代價函數為

其中，ui為殘差向量u=Lt+1?M t+1xt+1的第i個分量；ρ為Huber 函數，表達式為

其中，γ為可調節因子，取值為1～2，本文中，γ=1.345。因為在純高斯分布的條件下，該值的效率是基于L2范數估計的95%。

當代價函數存在極小值時，本文可以通過對代價函數求偏導的方式對其進行求解。因此對于式(16)，在等式左右兩邊分別同時對xt+1求偏導，并令其結果等于0。

將上述量測數據值和預測狀態值之間的線性回歸問題添加到卡爾曼濾波遞推公式的量測更新方程中，就得到了Huber-M 估計的穩健卡爾曼濾波，簡稱RKF。圖3 為RKF 狀態模型與遞推轉換的流程框架。

圖3 RKF 狀態模型與遞推轉換的流程框架

2.3 遞推參數初值確定

由于短波通信頻率的日變化性，因此需要對預報前一個月的數據進行日變換統計分析。參考文獻[22]中遞推參數初值的確定方法，初始化P0、Q、R這3 個參數的初始值。

P0為初始系統狀態的噪聲協方差矩陣，由于X0是從實測的樣本數據獲取的，嚴格為系統真值，因此認為P0=0，即為零矩陣。

對于動態隨機噪聲和量測隨機噪聲，如果采用經驗進行分析估計，則具有隨機性，不便于實際模型的輸入。于是本文采用回歸系數估計的方法，建立傾斜鏈路探測站8 月和9 月的回歸方程，根據2 個月所得回歸方程的系數變化進行估算，確定Q矩陣為

3 數據來源和誤差分析方法

分析用到的電離層MUF 樣本數據采用重慶、海南、昆明3 個電離層觀測站的傾斜探測數據，該數據由中國電波傳播研究所計算中心提供。由于短波在單個頻率點進行通信時，通信時間具有連續性。結合該特性，本文在選擇探測樣本數據時，每小時內選取2 個探測頻率點作為驗證數據，使檢驗結果更具代表性。

采用均方根誤差（ERMS）和相對誤差（ΔPD）這2 個指標[23]來評估與分析預報結果的好壞，指標定義分別為

其中，MUFob為某個時刻的探測值，MUFpre為相應時刻的預報值，N為所選樣本數據的總數。

4 預報結果比對與分析

圖4～圖7 分別為2019 年10 月4 日—6 日，重慶到海南、昆明探測站（正向、反向）通過傳統卡爾曼濾波（KF）、穩健卡爾曼濾波（RKF）和國際電離層參考模型（IRI）3 種預報方法得到的MUF 預報結果，同時給出不同方法在相同時刻頻率預報的誤差。

從預報結果來看，IRI 的預報值與實際觀測值相比偏差較大，但是在部分時間節點上也能取得不錯的預測效果；KF 和RKF 在整個預報過程中都保持著較高的精度和穩定性。

圖4 重慶站—海南站的MUF 預報結果

圖5 海南站—重慶站的MUF 預報結果

圖6 重慶站—昆明站的MUF 預報結果

圖7 昆明站—重慶站的MUF 預報結果

表1 3 種預報方法的結果分析

表1 分別給出了RKF、KF、IRI 這3 種預報方法在不同探測鏈路上預報結果的誤差分析。在重慶到海南（正、反向）沒有觀測野值的探測鏈路上，相比IRI，KF 的平均相對誤差降低了6.665%，平均均方根誤差減少了1.100 85 MHz；RKF 的平均相對誤差降低了6.525%，平均均方根誤差減少了1.117 35 MHz。對比均方根誤差和相對誤差可以得出，RKF 和KF 預報精度均高于IRI。將KF 和RKF兩者相比，平均的 ΔRMS=0.0165、ΔPD=0.14%，說明在無探測野值的鏈路上KF 和RKF 的預報精度相差不大。

在重慶到昆明探測鏈路上，相比IRI，KF 的相對誤差降低了2.48%，均方誤差減少了0.383 1 MHz；RKF 的相對誤差降低了2.08%，均方根誤差減少了0.331 5 MHz，ΔRMS=0.0516、ΔPD=0.4%。但從圖7(a)的探測數據看出，實際觀測的第一天和第三天在21:00—24:00（夜間）時間段，探測到的最大可用頻率異常偏高，根據實測經驗夜間短波通信頻率較低，該時間段存在探測野值。對于該條鏈路，若采用KF，均方根誤差為1.373 4 MHz，比IRI 高出0.298 7 MHz，濾波性能下降；若采用RKF，則比IRI 的相對誤差降低了5.96%，均方根誤差也減少了0.468 3 MHz，ΔRMS=0.767、ΔPD=3.59%。結合重慶到昆明的數據分析結果可知，RKF 在有觀測野值的情況下性能明顯優于KF，穩定性和抗干擾能力更強，具有更好的穩健性。通過圖7(b)誤差結果分析也可以看出，在存在觀測野值的情況下，RKF 的濾波性能最好，KF 次之，IRI 稍差。

最后，根據2019 年10 月多個電離層探測站點探測到的數據，采用本文提出的RKF 短期預報方法，獲得MUF 三維曲面，如圖8 所示。

選取MUF 頻點的85%作為測試頻點，在重慶到廣州（超過1 000 km）短波通信鏈路上進行測試分析。通過短波電臺，每30 min 呼叫一次，重慶到廣州24 h 可通率高于90%；每10 min 呼叫一次，重慶到廣州24 h 可通率高于87%。

5 結束語

本文提出的穩健卡爾曼濾波短期預報方法，進一步提高了電離層頻率預報方法的準確性。通過對實測MUF 數據變化規律及觀測誤差產生的成因與特征進行回歸分析，利用電離層參考模型作為先驗信息，改進卡爾曼濾波模型的狀態估計方程；引入代價函數機制，通過Huber-M 估計實現對預報狀態的量測更新，消除觀測野值帶來的污染。驗證結果表明，本文提出的穩健卡爾曼濾波方法在實際鏈路的頻率預報上均方根誤差和相對誤差均優于傳統的卡爾曼濾波方法，在昆明到重慶（存在觀測野值）的探測鏈路上，所提方法比KF 的相對誤差降低了3.59%，均方根誤差也減少了0.767 MHz。能夠對電離層的最大可用頻率實現更好的追蹤，具有更好的穩健性和穩定性。此外，根據多個電離層探測站點探測到的數據，采用本文提出的穩健卡爾曼濾波短期預報方法，獲得最高可用頻率曲面圖，最后通過短波電臺，進行了重慶到廣州數千米鏈路的實際通信測試，可通率達87%以上，具有良好的連通率。考慮到短波信道的特點，后續工作將結合不同季節電離層的信道特性、磁暴現象等因素展開進一步的分析與研究。

圖8 不同時間、不同距離最高可用頻率三維曲面