短波通信選頻中的電離層同化短期預報方法

湯 軍 姬生云 王 健 王先義

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

引 言

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

1 基于同化技術的短期預報方法

1.1 數據同化方法

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

1.2 Kalman濾波同化原理

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

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

Kalman濾波同化預報方程可表示為

xa=xb+K(x0-Hxb).

(1)

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

K=PbHT(HPbHT+R)-1.

(2)

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

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

(3)

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

1.3 foF2短期預報

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

圖1 foF2短期預報流程圖

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

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

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

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

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

(4)

(5)

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

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

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

(6)

令b=r·Pb，進而可得

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

(7)

(8)

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

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

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

(9)

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

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

2 預報方法評估

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

表1 太陽活動期選擇

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

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

圖2 特征尺度對預報結果的影響

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

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

圖3 太陽活動高年預報結果與實測數據對比圖

圖4 太陽活動中期預報結果與實測數據對比圖

圖5 太陽活動低年預報結果與實測數據對比圖

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

(10)

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

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

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

(11)

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

表3 相對誤差的統計分析

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

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

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

表4 不同預報尺度的誤差分析統計

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

3 結 論

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

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

[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] 劉瑞源, 劉順林, 徐中華, 等. 自相關分析法在中國電離層短期預報中的應用[J]. 科學通報, 2002, 50(24): 2781-2785.

[4] 馮 靜, 柳 文, 焦培南等. 電離層特征參量的自相關原理插值方法[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] 陳 春, 吳振森, 趙振維, 等. 基于神經網絡技術的foF2短期預報方法[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] 陳 春, 吳振森, 孫樹計, 等. 利用神經網絡提前24小時預報電離層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短期預報電離層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] 樂新安, 萬衛星, 劉立波, 等. 基于Gauss-Markov 卡爾曼濾波的電離層數值同化現報預報系統的構件-以中國及周邊地區為例德觀測系統模擬實驗[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濾波電離層短期預報[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] 柳 文, 焦培南, 馮 靜, 等. 電離層參數的相似日短期預測方法[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] 馮 靜, 柳 文, 凡俊梅, 等. 一種電離層短期預報的新方法[J]. 裝備環境工程, 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]. 電波科學學報, 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] 劉成思, 薛紀善. 關于集合Kalman濾波的理論和方法的發展[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] 王 健, 惠守強, 付 煒, 等. 高頻單站定位誤差特性分析及精度優化構想[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)