李國軍，鄭廣發(fā)，葉昌榮，周銀萍

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

1 引言

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

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

2 卡爾曼濾波

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

2.1 傳統(tǒng)卡爾曼濾波預測模型

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

圖1 模型框架

模型表達式為[20]

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

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

協(xié)方差矩陣量測更新為

KF 狀態(tài)和遞推轉(zhuǎn)移流程如圖2 所示。

圖2 KF 狀態(tài)和遞推轉(zhuǎn)移流程

2.2 Huber-M 穩(wěn)健卡爾曼濾波方法

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

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

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

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

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

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

聯(lián)合式(2)構(gòu)造線性回歸模型為

定義矩陣

則式(10)可以改寫為

定義Huber-M 估計的代價函數(shù)為

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

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

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

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

圖3 RKF 狀態(tài)模型與遞推轉(zhuǎn)換的流程框架

2.3 遞推參數(shù)初值確定

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

P0為初始系統(tǒng)狀態(tài)的噪聲協(xié)方差矩陣，由于X0是從實測的樣本數(shù)據(jù)獲取的，嚴格為系統(tǒng)真值，因此認為P0=0，即為零矩陣。

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

3 數(shù)據(jù)來源和誤差分析方法

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

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

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

4 預報結(jié)果比對與分析

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

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

圖4 重慶站—海南站的MUF 預報結(jié)果

圖5 海南站—重慶站的MUF 預報結(jié)果

圖6 重慶站—昆明站的MUF 預報結(jié)果

圖7 昆明站—重慶站的MUF 預報結(jié)果

表1 3 種預報方法的結(jié)果分析

表1 分別給出了RKF、KF、IRI 這3 種預報方法在不同探測鏈路上預報結(jié)果的誤差分析。在重慶到海南（正、反向）沒有觀測野值的探測鏈路上，相比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)的探測數(shù)據(jù)看出，實際觀測的第一天和第三天在21:00—24:00（夜間）時間段，探測到的最大可用頻率異常偏高，根據(jù)實測經(jīng)驗夜間短波通信頻率較低，該時間段存在探測野值。對于該條鏈路，若采用KF，均方根誤差為1.373 4 MHz，比IRI 高出0.298 7 MHz，濾波性能下降；若采用RKF，則比IRI 的相對誤差降低了5.96%，均方根誤差也減少了0.468 3 MHz，ΔRMS=0.767、ΔPD=3.59%。結(jié)合重慶到昆明的數(shù)據(jù)分析結(jié)果可知，RKF 在有觀測野值的情況下性能明顯優(yōu)于KF，穩(wěn)定性和抗干擾能力更強，具有更好的穩(wěn)健性。通過圖7(b)誤差結(jié)果分析也可以看出，在存在觀測野值的情況下，RKF 的濾波性能最好，KF 次之，IRI 稍差。

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

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

5 結(jié)束語

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

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