徐 彬 吳永宏 劉毅敏

（1.海軍電磁頻譜管理辦公室，北京100841；2.中國電波傳播研究所，山東 青島266107）

1.引 言

最高可用頻率指實際通信中能被電離層反射回地面的電波的最大頻率[1]，若選用的工作頻率超過它，則電波穿出電離層，不再返回地面；若選用的工作頻率太低，電波傳輸損耗較大，通信效果較差[2]。將工作頻率選定在最高可用頻率附近通常會帶來較好的通信效果[3]，因此預(yù)報最高可用頻率對短波通信系統(tǒng)合理選取工作頻率具有重要的意義。

通常，可以基于電離層模式進(jìn)行最高可用頻率的長期預(yù)報[4－6]，也可以基于電離層實際測量數(shù)據(jù)進(jìn)行最高可用頻率的短期預(yù)報。Zolesi等人利用有效太陽黑子數(shù)代替R12，進(jìn)而確定F2層臨界頻率f0F2，然后結(jié)合電離層傳播因子 M（3000）F2對鏈路的最高可用頻率進(jìn)行短期預(yù)報[7－10]。Pietrella等人 提 出 了ISWIRM 模 型 （Instantaneous Space Weighted Ionospheric Regional Model），該模型利用歐洲地區(qū)的4個探測站的實時數(shù)據(jù)預(yù)報F2層臨界頻率f0F2，僅適用于北緯35°到70°西經(jīng)5°到東經(jīng)40°所覆蓋的區(qū)域[11]。黃德耀[12]和張訓(xùn)械[13]分別采用線性相關(guān)和神經(jīng)網(wǎng)絡(luò)方法預(yù)報最高可用頻率。

Kalman濾波技術(shù)作為一種數(shù)學(xué)工具，從1960年提出以來得到了快速發(fā)展，在衛(wèi)星導(dǎo)航、飛行器控制、數(shù)據(jù)通信等領(lǐng)域得到了廣泛應(yīng)用[14－17]。為了根據(jù)鏈路當(dāng)前時刻的探測數(shù)據(jù)預(yù)報未來時刻的最高可用頻率，本文在前人工作的基礎(chǔ)上引入了Kalman濾波技術(shù)，設(shè)計了最高可用頻率自適應(yīng)預(yù)報算法，并用中國電波觀測網(wǎng)的實測數(shù)據(jù)對該算法的性能進(jìn)行了檢驗，取得了較好的效果。

2.Kalman濾波技術(shù)[9]

Kalman濾波模型用狀態(tài)空間來描述數(shù)學(xué)公式，而且它的解也是遞歸計算的，可以不加修改地應(yīng)用于平穩(wěn)和非平穩(wěn)環(huán)境，其狀態(tài)方程為

式中：x（n）為M×1維狀態(tài)向量；A（n＋1，n）為M×M維轉(zhuǎn)移矩陣；ν1（n）為M×1維高斯白噪聲；n表示時刻。

測量方程為

式中：y（n）為N×1維系統(tǒng)輸出向量；C（n＋1，n）為N×M維測量矩陣；ν2（n）為N×1維高斯白噪聲；n表示時刻。

令Yn為y（1），y（2），…，y（n）張成的空間，新息α（n）＝y(tǒng)（n）－（n｜Yn－1），狀態(tài)預(yù)測誤差向量為ε（n，n－1）＝x（n）－（n｜Yn－1），Q2（n）為測量噪聲向量ν2（n）的相關(guān)矩陣，Q1（n）為狀態(tài)噪聲向量ν1（n）的相關(guān)矩陣，狀態(tài)預(yù)測誤差相關(guān)矩陣為K（n，n－1）＝E[ε（n，n－1）εH（n，n－1）] ，對式（1）和（2）進(jìn)行推導(dǎo)，可得狀態(tài)向量的估計值為

式中Kalman增益矩陣為

其中

利用Riccati差分方程可得狀態(tài)預(yù)測誤差相關(guān)矩陣的遞歸計算表達(dá)式為

式中

通過分析式（1）～（7），可以看出轉(zhuǎn)移矩陣A（n＋1，n）、測量矩陣C（n＋1，n）、測量噪聲相關(guān)矩陣Q2（n）、狀態(tài)向量噪聲相關(guān)矩陣Q1（n）是已知的，配置好初始條件＾x（1｜Y0）、K（1，0）即可利用式（3）～（7）進(jìn)行遞歸計算，估算狀態(tài)向量x（n）。

3.最高可用頻率自適應(yīng)預(yù)報

電離層斜向探測設(shè)備每半個小時測量到一幅電離圖，從該電離圖可以讀出最高可用頻率，把讀出的最大可用頻率按時間進(jìn)行排序，就形成了一個時間序列，實踐證明該序列是一個平穩(wěn)隨機過程。因此，可以根據(jù)電波觀測網(wǎng)測量的已知數(shù)據(jù)預(yù)報特定鏈路在未來時刻的最高可用頻率。

若最高可用頻率觀測數(shù)據(jù)用y（n）來表示，利用平穩(wěn)時間序列自回歸模型可得

式中：φ1，φ2，…，φp為待定系數(shù)；an為高斯白噪聲。

令Y（n－1）＝[y（n－1），y（n－2），…，y（n－p）] ，φ（n－1）＝[φ1，φ2，…，φp]T，則式（8）可變?yōu)?/p>

而且

將式（9）、（10）和式（1）、（2）進(jìn)行對比，可得

把式（11）～（13）代入式（3）～（7），可得這樣式（14）～（17）就構(gòu)成了遞歸求解算法，其計算流程如圖1所示，圖中測量噪聲相關(guān)矩陣Q2和測量數(shù)據(jù)y（n）為已知參數(shù)（0）、K（0）為遞歸初始條件。也就是說，通過遞歸計算可以獲得較為準(zhǔn)確的向量＾φ（n），那么未來時刻最高可用頻率的預(yù)測值為

圖1 自適應(yīng)預(yù)報算法流程

把上述算法稱為自適應(yīng)預(yù)報算法，該算法可以根據(jù)實際測量數(shù)據(jù)的變化自動修正自回歸模型的待定系數(shù)，然后對未來時刻的最高可用頻率進(jìn)行預(yù)報。

4.實測數(shù)據(jù)檢驗

我們收集了2010年5月到2010年10月中國電波觀測網(wǎng)電離層斜向探測設(shè)備測量的北京－青島鏈路電離圖，從這些圖中讀取最高可用頻率，然后用這些實測數(shù)據(jù)對自適應(yīng)預(yù)報算法的性能進(jìn)行檢驗。

圖2 自適應(yīng)預(yù)報算法跟蹤過程

令p＝2，測量誤差相關(guān)矩陣Q2＝1，＾φ（0）＝[1，，利用自適應(yīng)預(yù)報算法給出的結(jié)果如圖2所示，圖中直線為2010.5～2010.10期間每日13點對應(yīng)的最高可用頻率的觀測值，點線為預(yù)報值，可以看出，自適應(yīng)預(yù)報算法在第12個數(shù)值點附近就跟蹤上了觀測值，而且隨后較準(zhǔn)確地重現(xiàn)了觀測值。

定義預(yù)報誤差的均方差為

式中：i為預(yù)報值下標(biāo)；N為預(yù)報值總數(shù)量；e（i）為預(yù)報相對誤差，其表達(dá)式為

且y（i）為觀測值（i）為預(yù)報值。

表1為根據(jù)自適應(yīng)預(yù)報算法給出的一天中每個小時預(yù)報誤差的均方差，表2為根據(jù)中國參考電離層模型給出的一天中每個小時預(yù)報誤差的均方差，圖3為這兩個均方差的圖形化顯示。可以看出，與長期預(yù)報算法相比，基于實測數(shù)據(jù)的自適應(yīng)預(yù)報算法可以給出更精確的結(jié)果。

表1 自適應(yīng)預(yù)報算法給出的均方差

表2 長期預(yù)報算法給出的均方差

圖3 一天中均方差變化曲線

5.結(jié) 論

本文從短波通信系統(tǒng)和短波鏈路需要預(yù)報或指配優(yōu)質(zhì)頻率的工程實際出發(fā)，引入了Kalman濾波技術(shù)和最高可用頻率觀測數(shù)據(jù)自回歸模型，利用Kalman濾波技術(shù)推導(dǎo)自回歸模型中待定系數(shù)的求解方法，進(jìn)而給出了最高可用頻率自適應(yīng)預(yù)報算法。利用中國電波觀測網(wǎng)測量的北京－青島鏈路數(shù)據(jù)對該算法的性能進(jìn)行檢驗，結(jié)果表明最高可用頻率自適應(yīng)預(yù)報算法能夠從人為設(shè)定的初值較快跟蹤上自回歸模型的待定系數(shù)，此后基于實際觀測數(shù)據(jù)較準(zhǔn)確地預(yù)報了未來時刻的最高可用頻率，其預(yù)報誤差的均方差遠(yuǎn)低于中國參考電離層模型給出的結(jié)果。

[1] 沈琪琪，朱德生.短波通信[M] .西安：西安電子科技大學(xué)出版社，1989.

[2] 熊年祿，唐存琛，李行健.電離層物理概論[M] .武漢：武漢大學(xué)出版社，1999.

[3] GOODMAN J M.HF communication：Science and Technology[M] .Van Nostrand，New York，1992.

[4] WANG Jian，F(xiàn)ENG Xiaozhe，CHEN Li.Basic MUF observation and comparison of HF radio frequency prediction based on different ionosphere models[C] //9th International Symposium on Antenna，Propagation and EM Theory（ISAPE），2010，403－406.

[5] HF propagation prediction method[R] .ITU－R 533－7，2001.

[6] 中國參考電離層[M] .GJB1925－1994，1994.

[7] ZOLESI B，BELEHAKI A，TSAGOURI I，et al.Real－time updating of the simplified ionosphere regional model for operational applications[J] .Radio Science，2004，139（2），doi：10.1029/2003RS002936.

[8] TSAGOURI I，ZOLESI B，BELEHAKI A，et al.E－valuation of the performance of the real－time updated simplified ionospheric regional model for the European area[J] .Journal of Atmospheric and Solar－Terrestrial Physics，2005，67（12）：1137－1146.

[9] ZOLESI B，CANDER L R，De and FRANCESCHI G.Simplified ionospheric regional model for telecommunication applications[J] .Radio Science，1993，28（4）：603－612.

[10] LOCKWOOD M.A simple M－factor algorithm for improved estimation of the basic maximum usable frequency of radio waves reflected from the ionospheric F region[J] .IEE Proceedings－F：Communication，Radar and Signal Processing，1983，130（4）：296－302.

[11] PIETRELLA M，PERRONE L.Instantaneous space weighted ionospheric regional model for instantaneous mapping of the critical frequency of the F2 layer in the European region[J] .Radio Science，2005，40（1），doi：10.1029/2003RS003008.

[12] 黃德耀.短波電路最高可用頻率準(zhǔn)實時預(yù)報方法[J] .通信學(xué)報，1991，12（6）：105－107.HUANG Deyao.Quasi－real－time forecast method of the shortwave circuit MUF[J] .Journal on Communications，1991，12（6）：105－107.（in Chinese）

[13] ZENG Wen，ZHANG XunJie.Prediction of HF communication MUF in the region of south china sea[J] .IEEE Antennas and Propagation Magazine，1999，41（4）：35－38.

[14] CARMI A，GURFIL P，KANEVSKY D.Methods for sparse signal recovery using Kalman filtering with embedded pseudo－measurement norms and quasinorms[J] .IEEE transactions on Signal Processing，58（4）：2405－2409，2010.

[15] MURALI R.Data－based Techniques to Improve State Estimation in Model Predictive Control[D] .University of Wisconsin－Madison，October 2007.

[16] 鄭寶玉，等譯，Simon Haykin著.自適應(yīng)濾波器原理（第四版）[M] .北京：電子工業(yè)出版社，2006.

[17] KALMAN R E.A new approach to linear filtering and prediction problems[J] .Journal of Basic Engineering，1960，82（1）：35－45.