桂曉純，洪振杰，趙運(yùn)超

(溫州大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，浙江溫州 325035)

基于Kalman濾波技術(shù)的電離層電子密度反演研究

桂曉純，洪振杰?，趙運(yùn)超

(溫州大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，浙江溫州 325035)

利用Kalman濾波技術(shù)，實(shí)現(xiàn)了基于6個(gè)衛(wèi)星接收站所獲取三頻信標(biāo)數(shù)據(jù)的電離層電子密度反演，并對(duì)影響反演效果的因素進(jìn)行了分析．結(jié)果表明：初始值、水平相關(guān)度以及垂直相關(guān)度等因素都會(huì)對(duì)電離層電子密度反演效果產(chǎn)生重要影響，其中水平相關(guān)設(shè)為14°、垂直相關(guān)設(shè)為200 km時(shí)反演效果較好．

Kalman濾波；電離層；衛(wèi)星三頻信標(biāo)；總電子含量

三頻信標(biāo)是一種重要的空間環(huán)境探測(cè)星載設(shè)備，其工作機(jī)制是發(fā)射機(jī)在VHF、UHF和L頻段上輸出相對(duì)穩(wěn)定相位的載頻信號(hào)，并經(jīng)全向天線向預(yù)定覆蓋區(qū)域輻射，通過(guò)地面接收設(shè)備在單站或者多站進(jìn)行信號(hào)幅度和相位測(cè)量，再進(jìn)行反演可以獲得電離層電子總含量、電離層閃爍和電子密度及電離層不規(guī)則體的信息[1-3]．

美國(guó)海軍實(shí)驗(yàn)室的Bernhardt等人[4]創(chuàng)造性的工作使得三頻信標(biāo)技術(shù)有了重大突破．他們巧妙地利用數(shù)論知識(shí)，經(jīng)過(guò)嚴(yán)格的數(shù)學(xué)推導(dǎo)，得到了TEC（總電子含量，Total Electron Content，簡(jiǎn)稱(chēng)TEC）的計(jì)算方法，使得其相位模糊系數(shù)值大幅度提高，大大降低了相位積分常數(shù)的求解難度．中國(guó)臺(tái)灣和美國(guó)于 1997年聯(lián)合制定了 COSMIC計(jì)劃（氣象、電離層和氣候的衛(wèi)星觀測(cè)系統(tǒng)，Constellation Observing System for Meteorology, Ionosphere and Climate，簡(jiǎn)稱(chēng)COSMIC），并于2006年成功發(fā)射衛(wèi)星．COSMIC攜帶有美國(guó)海軍實(shí)驗(yàn)室研制的TBB（三頻信標(biāo)，Tri-Band Beacon，簡(jiǎn)稱(chēng)TBB）發(fā)射機(jī)，三個(gè)載波分別為150 MHz、400 MHz和1 067 MHz，它們發(fā)射的信號(hào)可以被地面或者其它低軌衛(wèi)星接收機(jī)接收．通過(guò)不同頻率之間測(cè)量的相位差就可以算出發(fā)射機(jī)到接收機(jī)路徑上的總電子含量．將地面接收機(jī)得到的總電子含量數(shù)據(jù)作為輸入?yún)⒘坷肅IT技術(shù)（電離層層析成像，Computerized Ionospheric Tomography，簡(jiǎn)稱(chēng)CIT）就可以反演電離層電子密度[5-9]．三頻信標(biāo)的精確探測(cè)的特點(diǎn)也逐漸引起我國(guó)地震電磁工作者的關(guān)注，經(jīng)過(guò)幾年的論證和實(shí)驗(yàn)，中國(guó)電波傳播研究所吳健研究員、甄衛(wèi)民研究員等人[10-12]在三頻信標(biāo)探測(cè)原理和 CT（層析成像，Computerized (or Computed) Tomography，簡(jiǎn)稱(chēng)CT）算法上有所突破，他們的工作為我國(guó)三頻信標(biāo)的研究奠定了基礎(chǔ)．雖然三頻信標(biāo)TEC探測(cè)和CT成像技術(shù)在各方面都取得了重大的進(jìn)展，但是我們需要看到的是不管是理論還是算法，三頻信標(biāo)CT技術(shù)還遠(yuǎn)沒(méi)有達(dá)到理想的要求，都需要進(jìn)一步的發(fā)展．

本文基于Kalman濾波技術(shù)，在介紹三頻衛(wèi)星信標(biāo)測(cè)量TEC原理的基礎(chǔ)上，實(shí)現(xiàn)了三頻信標(biāo)數(shù)據(jù)資料的反演，得到了2維的電離層電子密度反演圖像，并對(duì)影響反演效果的一些因素進(jìn)行了討論，本研究對(duì)于三頻信標(biāo)數(shù)據(jù)獲取站位布設(shè)及反演相關(guān)條件設(shè)定具有重要現(xiàn)實(shí)指導(dǎo)意義．

1 Kalman濾波反演算法及原理

Kalman濾波是Kalman在1960年提出的從與被提取信號(hào)有關(guān)的觀測(cè)量中通過(guò)算法估計(jì)出所需要的信號(hào)的一種濾波算法[13]．Kalman濾波方法是一種時(shí)域方法，對(duì)于具有高斯分布噪聲的線性系統(tǒng)可以得到系統(tǒng)狀態(tài)的遞推最小方差估計(jì)（Recursive Minimum Mean-Square Estimation, RMMSE）．

Kalman濾波是一種最優(yōu)估計(jì)理論，與傳統(tǒng)的估計(jì)方法相比，它具有以下幾個(gè)特點(diǎn)：

第一，在時(shí)域內(nèi)設(shè)計(jì)濾波器，并且算法是遞推的，每一步只需要處理一個(gè)時(shí)刻的量測(cè)信息，既綜合利用了t時(shí)刻以前的全部量測(cè)信息又不會(huì)使計(jì)算量隨時(shí)間增大．

第二，設(shè)計(jì)Kalman濾波器采用了狀態(tài)空間法，對(duì)于多維隨機(jī)過(guò)程的估計(jì)也適用．采用動(dòng)態(tài)方程來(lái)描述被估計(jì)量的動(dòng)態(tài)變化規(guī)律，被估計(jì)量的動(dòng)態(tài)統(tǒng)計(jì)信息由激勵(lì)白噪聲的統(tǒng)計(jì)信息和動(dòng)態(tài)方程確定．由于激勵(lì)白噪聲是平穩(wěn)過(guò)程，而且動(dòng)態(tài)方程是已知的，所以被估計(jì)量可以是平穩(wěn)的也可以是非平穩(wěn)的，所以Kalman濾波也適用于非平衡過(guò)程．

假設(shè)Y為電子密度的觀測(cè)值，對(duì)于電離層電子總含量，它和電子密度的關(guān)系為：

其中T(r1,r2)表示從衛(wèi)星位置r1到接收器位置r2的電子總含量，Ne為電子密度．

把電子密度作為一個(gè)三維網(wǎng)格區(qū)域值可以得到一個(gè)離散近似為：

H中的元素Hi,k表示射線i經(jīng)過(guò)網(wǎng)格點(diǎn)k處射線長(zhǎng)度．

由于存在不確定的模糊度，TEC的測(cè)量通?？醋魇窍鄬?duì)測(cè)量，而不是絕對(duì)測(cè)量．通過(guò)將每個(gè)臺(tái)站的測(cè)量值減去一個(gè)特定的測(cè)量值，轉(zhuǎn)換成絕對(duì)測(cè)量量T1：

T0為給定接收器對(duì)指定衛(wèi)星的參考測(cè)量，H0為參考測(cè)量的射線路徑．

為簡(jiǎn)化推導(dǎo)過(guò)程，在接下來(lái)的討論中默認(rèn)H即為H-H0．在本文討論中都是經(jīng)度固定在120°時(shí)討論的，所以H可以看成一個(gè)二維的矩陣．

在數(shù)據(jù)空間，分析電子密度X和初始電子密度X0具有如下的關(guān)系：

其中P為背景協(xié)方差矩陣，X0為初始電子密度，R為觀測(cè)協(xié)方差．（5）式中PHT[R+HPHT]-1(Y-HX)在Kalman濾波中通常稱(chēng)為Kalman增益．

分析協(xié)方差與背景協(xié)方差具有如下關(guān)系：

P為背景協(xié)方差，aP為分析協(xié)方差，在迭代的過(guò)程中為下一時(shí)刻的背景協(xié)方差或者將其作為背景協(xié)方差的一個(gè)重要參考．

網(wǎng)格點(diǎn)k和l之間的背景協(xié)方差可以寫(xiě)為：

其中r表示經(jīng)驗(yàn)誤差估計(jì)，x0k和x0l分別表示k和l處的初始電子密度，Z表示高度，rkl表示點(diǎn)k和l之間緯度之差，Lz表示垂直相關(guān)度，Lα表示水平相關(guān)度．

2 實(shí)驗(yàn)過(guò)程及其結(jié)果分析

2.1 模擬反演實(shí)驗(yàn)準(zhǔn)備

假設(shè)低軌衛(wèi)星在800 km高空沿120°經(jīng)度圈自北向南飛行，在中國(guó)境內(nèi)沿衛(wèi)星飛行軌道（即120° E）設(shè)置6個(gè)衛(wèi)星接收站，站臺(tái)緯度依次設(shè)為2°、5°、8°、11°、14°和17°，因此可以以高度和緯度進(jìn)行網(wǎng)格劃分，得到一個(gè)2維的區(qū)域．高度選取100 km至800 km，間隔20 km．緯度選為0°到20°，間隔1°．R矩陣我們通過(guò)經(jīng)驗(yàn)給出，本實(shí)驗(yàn)取R=k·E，其中，E表示單位陣，其維數(shù)為所有衛(wèi)星信號(hào)射線數(shù)量之和，k= 0.1× dim（E）．

對(duì)于系數(shù)矩陣的算法，此處假定的算法是：從第一條射線的 TEC開(kāi)始，將每一條射線的TEC減去所有射線中最小的TEC，從而得到系數(shù)矩陣A．初始水平相關(guān)設(shè)為10°，垂直相關(guān)設(shè)為30 km，經(jīng)驗(yàn)誤差估計(jì)設(shè)為0.2．衛(wèi)星的高度設(shè)為800 km，衛(wèi)星緯度跨越為南緯3°至北緯23°，TEC采樣間隔約為0.064°．在本文研究中，真實(shí)的電子密度以及初始的電子密度均由NeQuick模型[14]給出．

2.2 初始場(chǎng)對(duì)反演結(jié)果的影響分析

首先我們分析初始值對(duì)反演結(jié)果的影響．真實(shí)值為1月份UT6時(shí)，f 10.7（太陽(yáng)10.7 cm輻射通量）為63；初始值為1月份UT6時(shí)，f 10.7分別設(shè)為68、73和78．得到電離層電子密度真實(shí)值與初始值以及它們的誤差圖像，如圖1所示．從圖1可以看出：初始值的設(shè)置對(duì)反演結(jié)果影響很大，初始值與真實(shí)值偏差越大，其反演結(jié)果的誤差越大．

2.3 誤差協(xié)方差矩陣的水平相關(guān)因子與垂直相關(guān)因子分析

誤差協(xié)方差矩陣通常是根據(jù)經(jīng)驗(yàn)給出，對(duì)反演結(jié)果的影響也比較大，其中水平和垂直相關(guān)因子與地理位置和晝夜變化相關(guān)．地理位置分低、中及高緯度討論，假設(shè)低緯度的范圍是0°到20°，中緯度是20°到40°，高緯度是40°到60°．對(duì)于不同的緯度，分為UT6時(shí)（下午）和UT18時(shí)（晚上）分別討論．

討論水平相關(guān)性時(shí)，取水平相關(guān)分別為4°、6°、8°、10°、12°、14°、16°、18°和20°；討論垂直相關(guān)性時(shí)，取垂直相關(guān)分別為30 km、50 km、70 km、100 km、130 km、150 km、170 km、200 km、230 km、250 km、270 km和300 km．

圖1 真實(shí)值圖像、f 10.7分別為68、73及78時(shí)初始值的圖像以及二者的誤差圖像

為了直觀的看出水平相關(guān)和垂直相關(guān)對(duì)電子密度反演結(jié)果的影響，截取200 km到500 km這一個(gè)區(qū)域，計(jì)算每個(gè)網(wǎng)格點(diǎn)的反演誤差，并統(tǒng)計(jì)總的平均誤差，用來(lái)比對(duì)反演結(jié)果的好壞，分別給出低、中以及高緯度相關(guān)性誤差分析的圖像，如圖2 – 圖4所示：

圖2 低緯度電子密度相關(guān)性的誤差分析圖像

圖3 中緯度電子密度相關(guān)性的誤差分析圖像

圖4 高緯度電子密度相關(guān)性的誤差分析圖像

對(duì)于低緯度，UT6時(shí)，隨著水平相關(guān)和垂直相關(guān)的增大，誤差整體上都呈現(xiàn)先增后減的趨勢(shì)；UT18時(shí)，誤差整體呈現(xiàn)減小的趨勢(shì)．

對(duì)于中緯度，UT6時(shí)，隨著水平相關(guān)和垂直相關(guān)的增大，誤差整體呈現(xiàn)減小的趨勢(shì)．UT18時(shí)，水平相關(guān)越大，誤差會(huì)隨之先減小后增大；垂直相關(guān)越大，誤差會(huì)隨之減?。?/p>

對(duì)于高緯度，UT6時(shí)，水平相關(guān)越大，誤差整體呈現(xiàn)減小的趨勢(shì)；垂直相關(guān)越大，誤差先減小后增大．UT18時(shí)，隨著水平相關(guān)和垂直相關(guān)的增大，誤差整體呈現(xiàn)先減小后增大的趨勢(shì)．

選取高緯度UT6時(shí)，將水平相關(guān)設(shè)為14°，垂直相關(guān)設(shè)為200 km，由公式（5）可以計(jì)算出每個(gè)網(wǎng)格點(diǎn)的電子密度，將其與真實(shí)電子密度先做差，再取絕對(duì)值，最后算出平均值為2.38%,反演效果圖如圖5所示．由圖5不難看出，初始值與真實(shí)值的誤差較明顯，但反演值與真實(shí)值的差別較小，說(shuō)明將水平相關(guān)設(shè)為 14°，垂直相關(guān)設(shè)為200 km時(shí)反演效果比較好．

圖5 高緯度UT6時(shí), 水平相關(guān)為14°、垂直相關(guān)為200 km的電子密度反演圖像

3 結(jié) 論

本文介紹了三頻信標(biāo)技術(shù)的國(guó)內(nèi)外研究現(xiàn)狀、Kalman濾波原理與算法，利用模擬實(shí)驗(yàn)對(duì)電離層電子密度進(jìn)行了反演，分析了初始值等因子對(duì)反演效應(yīng)的影響．研究發(fā)現(xiàn)初始值、水平相關(guān)度以及垂直相關(guān)度等因素都會(huì)對(duì)電離層電子密度反演效果產(chǎn)生重要影響，其中水平相關(guān)設(shè)為14°、垂直相關(guān)設(shè)為200 km時(shí)反演效果較好．

[1] Na H, Shen J, Lee H. A fourier domain technique for ionospheric tomography [J]. Radio science, 1995, 30(3):747-754.

[2] Garriott O K, Nichol A W. Ionospheric information deduced from the doppler shifts of harmonic frequencies from earth satellites [J]. Journal of Atmospheric and Terrestrial Physics, 1961, 22(1): 50-63.

[3] Titheridge J E. Determination of ionospheric electron content from the faraday rotation of geostationary satellite signals [J]. Planetary and Space Science, 1972, 20(3): 353-369.

[4] Bernhardt P A, Siefring C L. New satellite-based systems for ionospheric tomography and scintillation region imaging [J]. Radio science, 2006, 41(5): 3-14.

[5] 徐繼生, 葉公節(jié). 東亞赤道異常區(qū)電離層CT診斷: 實(shí)驗(yàn)及初步結(jié)果[J]. 地球物理學(xué)報(bào), 1995, 38(5): 553-563.

[6] 許正文, 吳健. Mexico草帽小波與電離層CT反演[J]. 電波科學(xué)學(xué)報(bào), 1999, 14(3): 308-312.

[7] 吳雄斌, 徐繼生. 一種改進(jìn)的電離層層析成像算法[J]. 地球物理學(xué)報(bào), 2000, 43(1): 19-28.

[8] 萬(wàn)衛(wèi)星, 寧百齊, 袁洪, 等. 電離層擾動(dòng)的GPS探測(cè)[J]. 空間科學(xué)學(xué)報(bào), 1998, 18(3): 247-252.

[9] 徐繼生, 鄒玉華, 馬淑英. GPS地面臺(tái)網(wǎng)和掩星觀測(cè)結(jié)合的時(shí)變?nèi)S電離層層析[J]. 地球物理學(xué)報(bào), 2005, 48(4): 759-767.

[10] 吳健. 用三頻衛(wèi)星信標(biāo)測(cè)量電離層天氣新方法[J]. 中國(guó)科學(xué): A輯, 2000, 30(S1): 111-114.

[11] 趙海生, 許正文, 吳健, 等. 三頻衛(wèi)星信標(biāo)測(cè)量TEC方法探討及數(shù)值模擬[J]. 裝備環(huán)境工程, 2010, 7(4): 9-13.

[12] 許正文, 吳健, 曹沖. 電離層CT變換重建方法中的小波外推技術(shù)[J]. 空間科學(xué)學(xué)報(bào), 1999, 19(2): 114-121.

[13] 秦永元, 張洪鉞, 汪叔華. 卡爾曼濾波與組合導(dǎo)航原理[M]. 西安: 西北工業(yè)大學(xué)出版社, 1998: 33-48.

[14] Di Giovanni G, Radicella S M. An analytical model of the electron density profile in the ionosphere [J]. Advances in Space Research, 1990, 10(11): 27-30.

Research of Ionospheric Electron Density Inversion Based on Kalman Filter Technology

GUI Xiaochun, HONG Zhenjie, ZHAO Yunchao
(College of Mathematics and Information Science, Wenzhou University, Wenzhou, China 325035)

This paper probes into the experiment that a tri-band beacon data obtained from six satellites and earth stations has been developed to realize ionospheric electron density inversion based on Kalman filter technology and the factors that affect the inversion effect is also analyzed. It turns out that the factors like initial value, horizontal correlation and vertical correlation have a significant effect on ionospheric electron density inversion. When horizontal correlation is set as 14 degree and vertical correlation is set as 200 kilometer, the inversion effect is the best.

Kalman Filter; Ionosphere; Satellite Tri-band Beacon; TEC (Total Electron Content)

P353

A

1674-3563(2015)04-0056-06

10.3875/j.issn.1674-3563.2015.04.010 本文的PDF文件可以從xuebao.wzu.edu.cn獲得

(編輯：封毅)

2014-12-22

公益性行業(yè)(地震)科研專(zhuān)項(xiàng)基金(201108004)

桂曉純(1990- )，男，安徽宣城人，碩士研究生，研究方向：應(yīng)用分析與最優(yōu)化理論．? 通訊作者，hong@wzu.edu.cn