湯 軍 姬生云 王 健 王先義

(1．海司信息化部，北京 100841；2．中國電波傳播研究所，山東 青島 266107)

引 言

電離層短期預(yù)報技術(shù)被廣泛應(yīng)用于高頻頻率管理，與可用頻率預(yù)測、頻率分配等方法成為短波鏈路通信選頻不可或缺的重要組成部分，為短波通信選頻提供重要支撐.當前，電離層垂直探測獲得的foE、foF2、M(3000)F2等電離層特性參數(shù)是短波鏈路通信選頻重要的基礎(chǔ)數(shù)據(jù)[1-2]；其中，foF2是短波可用頻率預(yù)測最主要的參數(shù)，對foF2預(yù)測的準確性直接影響了短波頻率預(yù)報的準確性，進而影響了短波頻率資源劃分、資源管理的可靠性.目前采用的方法主要有自相關(guān)函數(shù)法[3-4]，多元線性回歸法[5]，人工神經(jīng)網(wǎng)絡(luò)法[6-7]，等效太陽黑子數(shù)法[8]， 卡爾曼濾波法[9-10]、 相似日法[11-12]、暴時電離層預(yù)報方法[13]等等.上述方法各有優(yōu)勢，如自相關(guān)法在低于3 h的預(yù)報時間內(nèi)能夠取得高的預(yù)測精度，相似日在1～24 h情況下均取得較高的預(yù)測精度，且穩(wěn)定性高.然而這些方法都需要較長時間的數(shù)據(jù)積累以及一些其他的參數(shù)支撐，例如自相關(guān)方法需要20～30 d的數(shù)據(jù)，神經(jīng)網(wǎng)絡(luò)法需要5～7 d的電離層觀測數(shù)據(jù)以及地磁K指數(shù)和Ap指數(shù)等信息，而相似日法則需要預(yù)先利用自相關(guān)函數(shù)判斷相似性，再進行預(yù)測，所需數(shù)據(jù)量與自相關(guān)法接近.本文將卡爾曼同化技術(shù)引入電離層參數(shù)短期預(yù)報中，基于電離層的日變化規(guī)律利用同化技術(shù)，利用1～2 d的電離層觀測數(shù)據(jù)，實現(xiàn)電離層特征參數(shù)foF2的1～24 h的短期預(yù)報.

1 基于同化技術(shù)的短期預(yù)報方法

1.1 數(shù)據(jù)同化方法

數(shù)據(jù)同化方法的出現(xiàn)解決了兩類問題：一是如何把觀測和模式所帶來的兩種不同但又“互補”的信息融合起來[14]，二是如何利用同化思想進一步提高資料分析質(zhì)量和數(shù)據(jù)預(yù)報的準確性[15].方法有早期的統(tǒng)計訂正法[16](Successive Correction Method，SCM)、最優(yōu)插值法[17](Optimal Interpolation，OI)，以及現(xiàn)在較為成熟的變分同化法[18](Variational Analysis)，卡爾曼濾波(Kalman Filter，KF)和拓展卡爾曼濾波(Extended Kalman Filter，EnKF)[19]等等.

1.2 Kalman濾波同化原理

鑒于KF方法已經(jīng)成熟，且易于工程實現(xiàn)[20]的特點，本文選用了該方法進行數(shù)據(jù)同化分析.KF最早由Kalman于1960年首先提出，其基本思想是：首先進行模式狀態(tài)的預(yù)報，在此基礎(chǔ)上，根據(jù)觀測數(shù)據(jù)對模式狀態(tài)進行修正.Kalman濾波是一種有嚴格數(shù)學理論基礎(chǔ)的預(yù)報—校正統(tǒng)計優(yōu)化方法，這種優(yōu)化的特點是能夠不斷地將觀測資料同化到動態(tài)系統(tǒng)中，以數(shù)值模式為動力約束條件，在使觀測和模式結(jié)果均方誤差達到最小的條件下，得到基于以前所有觀測和當前觀測的系統(tǒng)變量的最優(yōu)線性估計.而且，該方法具有遞推特性，可以借助前一時刻的濾波結(jié)果推出當前時刻的狀態(tài)估計值，從而大大減少計算量和存儲量.

由Kalman濾波同化基本思路可知：隨著模式狀態(tài)預(yù)報的持續(xù)進行和新的觀測數(shù)據(jù)的陸續(xù)輸入，同化過程可不斷向前推進.

Kalman濾波同化預(yù)報方程可表示為

xa=xb+K(x0-Hxb).

(1)

式中：xb是模式預(yù)報向量(即背景場)；xo是觀測向量(即觀測場)；H是觀測算子，即模式預(yù)報向量與觀測向量轉(zhuǎn)換因子；K是增益矩陣，可表示為

K=PbHT(HPbHT+R)-1.

(2)

式中：R是觀測誤差的協(xié)方差矩陣；Pb是背景場誤差協(xié)方差. 通過這兩個參數(shù)就可以求出分析誤差協(xié)方差矩陣Pa.

Pa=Pb-PbHT(HPbHT+R)-1HPb.

(3)

簡單地講，同化預(yù)報的狀態(tài)xa是通過背景場xb和觀測增長(觀測場與背景場的差異)xo-Hxb之間的加權(quán)平均來估計的，權(quán)重是K，它起到用觀測增長在相對附近的模式預(yù)報點上來調(diào)整背景場的作用.一般而言，離觀測點越近的模式預(yù)報點要比遠離觀測點的模式預(yù)報點受到的影響要大.

1.3 foF2短期預(yù)報

利用上述思路，基于Kalman濾波同化理論的foF2預(yù)報過程如圖1所示.

圖1 foF2短期預(yù)報流程圖

1) 讀取觀測數(shù)據(jù)：讀取預(yù)報前8 h的數(shù)據(jù)作為觀測場xo；預(yù)報前9～32 h(觀測場前24 h)的數(shù)據(jù)為初始猜測場；讀取中國參考電離層預(yù)報所得的24 h月中值數(shù)據(jù)作為背景場xb.

2) 選擇特征尺度：特征尺度l是計算背景場相關(guān)函數(shù)的一個重要參數(shù)，其大小直接影響到了預(yù)報準確性.l的選擇需要根據(jù)實際的情況，通過對前期預(yù)報結(jié)果進行比較，選擇一個誤差相對較小的數(shù)值，作為最終確定的特征尺度.

3) 計算誤差協(xié)方差：分別計算背景場和觀測場的誤差協(xié)方差Pb和觀測誤差的協(xié)方差矩陣R.

4) 計算背景場的相關(guān)矩陣：利用歐氏距離公式計算背景場各個時刻數(shù)據(jù)之間的距離d，利用高斯相關(guān)函數(shù)計算相關(guān)矩陣.

(i,j= 1,…,4)，

(4)

(5)

式中，fi和fj分別為i和j時刻的foF2，單位為MHz.

5) 計算增益矩陣：其中關(guān)鍵是計算觀測算子H，由于電離層的日變化規(guī)律特性，同時背景場數(shù)據(jù)與觀測場數(shù)據(jù)都是foF2，變量參數(shù)也一致，因此，觀測算子是“完美的”，即

H=E?H(i,j)=1.

(6)

令b=r·Pb，進而可得

(i=1,…,8;j=1,…,24);

(7)

(8)

至此，可以根據(jù)式(2)計算卡爾曼系數(shù)矩陣，即增益矩陣.

6) 同化分析：因為觀測算子是完美的，所以有Hxb=xb，根據(jù)式(1)，獲得同化后的預(yù)報結(jié)果.在計算增益矩陣時，部分參數(shù)都已經(jīng)獲得，但在計算分析誤差協(xié)方差時還有一個未知量HPb，該變量的計算方法如下

(i=1,…,8;j=1,…,24),

(9)

就可以利用式(3)計算分析誤差協(xié)方差矩陣.

7) 判斷同化過程是否結(jié)束：設(shè)定同化過程結(jié)束條件，約束條件有兩個：一是設(shè)定一個合理的分析誤差協(xié)方差，每次同化過程結(jié)束時都判斷這個條件是否滿足，滿足則終止同化，不滿足則繼續(xù)同化；二是設(shè)定一個最多的同化次數(shù)，當超出這個同化次數(shù)時，自動終止同化預(yù)報過程,將同化的結(jié)果作為預(yù)報結(jié)果進行輸出，如果兩個條件都不滿足則將同化結(jié)果作為初始猜測場，跳入第一步，重新開始同化，直到兩個條件滿足一個為止.

2 預(yù)報方法評估

因電離層受太陽周期活動、季節(jié)等影響明顯，故選取了太陽活動高年、中期、低年三個典型時期進行分析，如表1所列.

表1 太陽活動期選擇

選取中國境內(nèi)北京、長春、廣州、海南、拉薩、蘭州、滿洲里、青島、烏魯木齊和重慶電離層觀測站點數(shù)據(jù).基于這三組實測數(shù)據(jù)利用Kalman同化技術(shù)對未來1～24 h的foF2進行預(yù)報.

特征尺度的差異會改變預(yù)報精度，圖2顯示了特征尺度從0.001到100變化，不同年份預(yù)報精度的相對誤差.

圖2 特征尺度對預(yù)報結(jié)果的影響

從圖2中可以看出，特征尺度對預(yù)報結(jié)果的影響較大，以2000年為例，預(yù)報誤差接近5%.因此，為了獲得較好的預(yù)報精度，2000年的特征尺度l=10，2003年的特征尺度l=10，2008年的特征尺度l=0.001.

根據(jù)上面確定的特征尺度，分別對太陽活動高年、中期以及低年三個不同月份進行1～24 h的預(yù)報，并與實測結(jié)果進行比較，如圖3、圖4和圖5所示：為了便于觀察分析，所有的預(yù)報都是從北京時間0點開始.

圖3 太陽活動高年預(yù)報結(jié)果與實測數(shù)據(jù)對比圖

圖4 太陽活動中期預(yù)報結(jié)果與實測數(shù)據(jù)對比圖

圖5 太陽活動低年預(yù)報結(jié)果與實測數(shù)據(jù)對比圖

對上述不同太陽活動時期的1～24 h的預(yù)報結(jié)果進行均方根誤差統(tǒng)計，如表2所列，均方根誤差的統(tǒng)計公式為

(10)

式中，Δfi為每個時刻預(yù)測值與實測值之間的絕對誤差，單位為MHz.

表2 均方根誤差分析表(單位MHz)

對上述不同太陽活動時期的1～24 h的預(yù)報結(jié)果進行相對誤差統(tǒng)計，如表3所列，相對誤差的統(tǒng)計公式為

(11)

式中，fi為每個時刻的實測值.

表3 相對誤差的統(tǒng)計分析

通過統(tǒng)計，太陽活動高年預(yù)報的均方根誤差最大，均值為0.82 MHz，中期次之，均值為0.70 MHz，低年最小，均值為0.50 MHz.結(jié)合表2可以看出，在2000年，太陽活動高年的平均均方根誤差最大，拉薩、海南和長春三站預(yù)報的均方根誤差尤其高，最高達到了1.74 MHz，其余各站則較小，滿洲里站最小為0.29 MHz；在太陽活動低年，廣州和海南站預(yù)報的均方根誤差較大，分別為0.99 MHz和1.27 MHz，其余各站誤差基本上都小于0.4 MHz；太陽活動中期預(yù)報的均方根誤差較為平穩(wěn)，除了海南、重慶和拉薩站外，其他各站的誤差都在0.5 MHz左右.

由表3可以看出，海南、拉薩和廣州站的預(yù)報精度稍差，相對誤差在15%左右，北京、青島和烏魯木齊的預(yù)報精度較高，相對誤差在8%左右，其他站相對誤差大體上10%左右，總體上的預(yù)報精度達到了10%.綜合比較均方根誤差和相對誤差可以發(fā)現(xiàn)，中高緯度地區(qū)的預(yù)報準確性要高于中低緯度地區(qū).

下面，對不同預(yù)報尺度下的預(yù)報誤差分析，預(yù)報尺度間隔分別為1～24 h，25～48 h以及49～72 h，如表4所示.

表4 不同預(yù)報尺度的誤差分析統(tǒng)計

從表4可以看出，均方根誤差和相對誤差的變化具有一致性，前24 h的預(yù)報誤差最小，25～48 h的預(yù)報誤差次之，49～72 h的預(yù)報誤差最大，即隨著預(yù)報時間的增長，預(yù)報準確性會下降.1～24 h預(yù)報的均方根誤差在0.57～0.72 MHz之間，相對誤差在9%～12%之間.根據(jù)相對誤差和預(yù)報結(jié)果的穩(wěn)定性，本文推薦的預(yù)報時間為1～24 h.與文獻[11]介紹的相似日分析方法相比，本方法對1～24 h的預(yù)報相對誤差為10%，與相似日方法的相對誤差基本相同.

3 結(jié) 論

本文以探測數(shù)據(jù)為基礎(chǔ)，利用Kalman濾波同化方法，提出了一種1～24 h的電離層參量foF2短期預(yù)報方法.并利用該方法對太陽活動高、低、中期年分不同季節(jié)不同地區(qū)預(yù)報結(jié)果進行分析.分析結(jié)果顯示了該方法在不同太陽活動周期時的預(yù)報準確度相近，對中高緯度地區(qū)的預(yù)報準確性較高，如青島、北京和滿洲里等地，同時，預(yù)報準確度隨著預(yù)報尺度的增加而下降，例如1～24 h的預(yù)報準確度優(yōu)于25～48 h，該方法對數(shù)據(jù)積累的時間要求低，更利于實現(xiàn)短波頻率預(yù)報工程化，為短波頻率管理提供可靠的支撐[1-2]，同時可為相關(guān)領(lǐng)域的擴展應(yīng)用提供支撐[21].

通過對數(shù)據(jù)誤差進行分析，該方法還有需要改進的地方，如中低緯地區(qū)的預(yù)報準確性較低.因此，借助于逐步成熟的集合Kalman同化以及4維變分同化等技術(shù)，提升中低緯度電離層短期預(yù)報準確性是下一階段研究工作的重點.

[1] ITU. P. 533-9 Method for the Prediction of the Performance of HF Circuits[S]. Geneva: ITU-R, 2008.

[2] ITU. P.1240-1 Methods of Basic MUF, Operational MUF and Ray-path Prediction[S]. Geneva: ITU-R, 2008.

[3] 劉瑞源, 劉順林, 徐中華, 等. 自相關(guān)分析法在中國電離層短期預(yù)報中的應(yīng)用[J]. 科學通報, 2002, 50(24): 2781-2785.

[4] 馮 靜, 柳 文, 焦培南等. 電離層特征參量的自相關(guān)原理插值方法[J]. 空間科學學報, 2009, 29(2): 195-201.

FENG Jing, LIU Wen, JIAO Peinan, et al. Autocorrelation method for interpolation of ionospheric characteristic parameters[J]. Chinese Journal of Space Science, 2009, 29(2): 195-201. (in Chinese)

[5] MARIN D, G MIRO A V MIKHAILOV. A method forfoF2short-term prediction, 4th COST 251 Workshop proceedings [C]// 4th COST 251 Workshop proceedings, Madeira, 1999: 214-222.

[6] 陳 春, 吳振森, 趙振維, 等. 基于神經(jīng)網(wǎng)絡(luò)技術(shù)的foF2短期預(yù)報方法[J]. 電波科學學報, 2008, 23(4): 708-712.

CHEN Chun, WU Zhensen, ZHAO zhenwei, et al. A short-termfoF2forecasting method based on neural network techniques[J]. Chinese Journal of Radio Science, 2008, 23(4): 708-712. (in Chinese)

[7] 陳 春, 吳振森, 孫樹計, 等. 利用神經(jīng)網(wǎng)絡(luò)提前24小時預(yù)報電離層foF2[J]. 電波科學學報, 2009, 24(1): 152-156.

CHEN Chun, WU Zhensen, SUN Shuji, et al. Forecasting the ionosphericfoF224 hour in advance by neural network techniques[J]. Chinese Journal of Radio Science, 2009, 24(1): 152-156. (in Chinese)

[8] 陳 春, 吳振森, 孫樹計, 等. 利用F10.7短期預(yù)報電離層foF2[J]. 空間科學學報, 2009, 29(4): 383-388.

CHEN Chun , WU Zhensen, SUN Shunji, et al. Short-term Forecasting of the IonosphericfoF2With F10.7[J]. Chin J Space Sci, 2009, 29(4): 383-388. (in Chinese)

[9] 樂新安, 萬衛(wèi)星, 劉立波, 等. 基于Gauss-Markov 卡爾曼濾波的電離層數(shù)值同化現(xiàn)報預(yù)報系統(tǒng)的構(gòu)件-以中國及周邊地區(qū)為例德觀測系統(tǒng)模擬實驗[J]. 地球物理學報, 2010, 53(4): 787-795.

YUE Xinan, WAN Weixing, LIU Libo, et al. Development of an ionospheric numerical assimilation nowcast and forecast system based on Gauss-Markov Kalman filter-an observation system simulation experiment taking example for China and its surrounding area[J]. Chinese Journal of Geophysics, 2010, 53(4): 787-795. (in Chinese)

[10] 徐 彤, 吳振森, 吳 健, 等. Kalman濾波電離層短期預(yù)報[J]. 電波科學學報, 2007, 22(增刊): 146-149.

XU Tong, WU Zhensen, WU Jian, et al. Ionospheric short-time forecasting using the Kalman filter[J]. Chinese Journal of Radio Science, 2007, 22(Sup): 146-149. (in Chinese)

[11] 柳 文, 焦培南, 馮 靜, 等. 電離層參數(shù)的相似日短期預(yù)測方法[J]. 電波科學學報, 2010, 25(2): 240-247.

LIU Wen, JIAO Peinan, FENG Jing, et al. Method of similar days for ionospheric parameter short-term forecasting[J]. Chinese Journal of Radio Science, 2010, 25(2): 240-247. (in Chinese)

[12] 馮 靜, 柳 文, 凡俊梅, 等. 一種電離層短期預(yù)報的新方法[J]. 裝備環(huán)境工程, 2009, 6(3): 15-20.

FENG Jing, LIU Wen, FAN Junmei, et al. A new method for short-term ionospheric forecast[J]. Equipment Environmental Engneering, 2009, 6(3): 15-20. (in Chinese)

[13] 孫樹計, 丁宗華, 陳 春, 等. 一種爆時電離層foF2的經(jīng)驗預(yù)報方法[J]. 電波科學學報, 2009, 24(4): 637-644.

SUN Shuji, CHEN Chun, DING Zonghua, et al. A storm-time empirical ionosphericfoF2prediction method[J]. Chinese Journal of Radio Science, 2009, 24(4): 637-644. (in Chinese)

[14] 劉成思, 薛紀善. 關(guān)于集合Kalman濾波的理論和方法的發(fā)展[J]. 熱帶氣象學報, 2005, 21(6): 628-633.

LIU Chengsi, XUE Jishan. The ensemble kalman filter theory and method development[J]. Journal of Tropical Meteorology, 2005, 21(6): 628-633.(in Chinese)

[15]BOUTTIER F, COURTIER P. Data Assimilation Concepts and Methods, Meteorological[R/OL]. 1999[2012-08-31]. http://www.ecmwf.int/newsevents/training/lecture_notes/pdf_files/ASSIM/Ass_cons.pdf

[16] BERGTHORSSON P, DONS B. Numerical weather map analysis. [J]. Tellus, 1955, 7(3): 329-240．

[17]GANDIN L S. Objective Analysis of Meteorology Fields[M]. Leningrad: Leningrad Hydromet Press, 1963.

[18] SASAKI Y. Some basic formalisms in numerical variational analysis[J]. Mon Wea Rev, 1970, 98: 875-883.

[19] KALMAN R E. A new approach to linear filter and prediction problems[J]. J Basic Eng, 1960, 82(3): 35-45.

[20] EVENSEN G. The ensemble Kalman filter: theoretical formulation and practical implementation[J]. Ocean Dymamics, 2003, 53: 343-367.

[21] 王 健, 惠守強, 付 煒, 等. 高頻單站定位誤差特性分析及精度優(yōu)化構(gòu)想[J]. 電波科學學報, 2010, 25(5): 925-933.

WANG Jian, HUI Shouqiang, FU Wei, et al. Error characteristic analysis and accuracy optimizing idea of HF single site location[J]. Chinese Journal of Radio Science, 2010, 25(5): 925-933. (in Chinese)