磁暴事件對高精度TEC二維分布的影響

余侯芳,郭文興,張玉強(qiáng),馮健,甄衛(wèi)民

(1. 中國電波傳播研究所,山東 青島 266107;2中興通訊股份有限公司,廣東 深圳 518057;3. 西安電子科技大學(xué) 物理與光電工程學(xué)院,陜西 西安 710071)

0 引 言

電離層總電子含量(TEC)是表征電離層特性的重要特征參量之一,研究其時(shí)空變化規(guī)律對衛(wèi)星通信、導(dǎo)航定位、空間天氣預(yù)報(bào)等領(lǐng)域有著極其重要的理論意義及應(yīng)用價(jià)值[1-2]．

實(shí)驗(yàn)探測中,電離層TEC值主要由地面接收站接收的衛(wèi)星信號(hào)獲得．在信號(hào)由衛(wèi)星向接收機(jī)傳播時(shí),由于電離層的色散性質(zhì),不同頻率信號(hào)所引起的偽碼時(shí)延和載波相位時(shí)延也不同,因此可以通過不同頻率信號(hào)間的差分求得信號(hào)路徑上的電離層總電子含量[3]．二十世紀(jì)七十年代以后,衛(wèi)星導(dǎo)航定位系統(tǒng)的出現(xiàn),極大地促進(jìn)了電離層探測技術(shù)的發(fā)展,其對電離層TEC的研究有著極其重要的作用及影響．由于全球定位系統(tǒng)(GPS)發(fā)展時(shí)間較早,衛(wèi)星星座多,接收站分布廣泛,因此基于GPS的電離層研究工作較多[4-5]．

目前,國際上已有許多電離層研究機(jī)構(gòu)實(shí)現(xiàn)了基于GPS的TEC地圖重構(gòu)數(shù)據(jù)產(chǎn)品的發(fā)布,包括歐洲定軌中心(CODE)、美國噴氣動(dòng)力實(shí)驗(yàn)室(JPL)以及歐洲空間局(ESA)等機(jī)構(gòu)分別完成了全球電離層地圖重構(gòu)算法研究,并于事后提供了基于IGS的全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)觀測站數(shù)據(jù)的全球電離層TEC地圖產(chǎn)品——GIM[6-7]．此外,日本、歐洲、美國、澳大利亞等國家和地區(qū)都建立精度更高、實(shí)時(shí)性更好的地區(qū)性TEC預(yù)報(bào)系統(tǒng)[8-11]．

近年來,國內(nèi)也有很多單位開展了GPS-TEC的研究與觀測工作．其中,中科院地質(zhì)與地球物理研究所(IGGCAS)開發(fā)了我國首個(gè)覆蓋全境的TEC現(xiàn)報(bào)系統(tǒng)[6]．毛田等[12]用Kriging方法構(gòu)建了中緯度區(qū)域的電離層TEC地圖．雷宵龍等[13-14]利用COSMIC數(shù)據(jù)對全球電離層TEC地圖進(jìn)行了構(gòu)建,證明了利用COMSIC掩星資料構(gòu)建電離層TEC地圖的可行性．

隨著全球電離層TEC地圖的建立,人們對電離層暴的理解與研究越來越深入,例如:Balan等[15-16]對中、低緯電離層暴的產(chǎn)生原因進(jìn)行了深入的探討;鄧忠新等[17]提出了TEC擾動(dòng)指數(shù)DI,借此分析了中國地區(qū)電離層TEC暴擾動(dòng)的產(chǎn)生緣由和分布特性;吳佳姝等[18]通過對2001—2010年間的156次磁暴事件的統(tǒng)計(jì)分析了歐洲扇區(qū)赤道到極光帶不同緯度的電離層暴特征．除此之外,國內(nèi)外學(xué)者對特定磁暴期間不同緯度地區(qū)的電離層的變化情況進(jìn)行了細(xì)致的分析[19-23],例如:Oluwaseyi等[19]重點(diǎn)研究了赤道地區(qū)電離層在磁暴期間的變化情況,劉海濤等[20]通過構(gòu)建北美地區(qū)的二維分布圖像,研究了一次強(qiáng)磁暴期間的電離層擾動(dòng)的變化情形及其原因．

雖然在電離層TEC領(lǐng)域已經(jīng)開展了大量的研究工作,也取得了不錯(cuò)的研究成果,但隨著科技的進(jìn)步,研究的發(fā)展,必然對電離層TEC的精度和分辨率提出更高的要求．本文提出了一種高精度TEC地圖重構(gòu)方法,利用亞大地區(qū)60個(gè)GPS臺(tái)站的實(shí)測數(shù)據(jù),基于Kriging內(nèi)插方法實(shí)現(xiàn)了亞大區(qū)域高分辨率TEC地圖重構(gòu),通過與實(shí)測數(shù)據(jù)進(jìn)行對比并對成像精度進(jìn)行了驗(yàn)證;在此基礎(chǔ)上,詳細(xì)分析了南北半球不同緯度的電離層TEC值在磁暴期間的不同表現(xiàn)．

1 Kriging插值方法

1.1 半變異函數(shù)(semi-variogram)

要使用Kriging法進(jìn)行插值,必須先求出半變異函數(shù),這是用于描述插值中所用數(shù)據(jù)之間空間相關(guān)性的函數(shù),它之所以重要是由于它是Kriging算法的主要輸入量．半變異函數(shù)可通過在固定距離dl±Δdl/2處對不同觀測數(shù)據(jù)做方差求得:

(1)

式中:γ*為經(jīng)驗(yàn)半變異函數(shù);m(dl)為在距離dl處的不同觀測對的數(shù)量,Zi和Zj是相應(yīng)點(diǎn)xi、xj在距離|xi-xj|=dl處的觀測值．

一旦求得經(jīng)驗(yàn)半變異函數(shù),下一步便是通過滿足各種數(shù)學(xué)條件轉(zhuǎn)化為可以應(yīng)用在Kriging方程中．

1.2 Kriging方法

Kriging是一種屬于最優(yōu)線性無偏估計(jì)的線性內(nèi)插方法,因此,Kriging方法的主要目的就是利用已知量的線性組合來估計(jì)未知變量Z*:

(2)

式中:ωi是權(quán)重因子,對應(yīng)著每個(gè)已知量Zi．

(3)

式中:E[(Z*-Z)2]是Z*的方差,能夠用半變異函數(shù)來表示:

(4)

減去方程(3),可得簡化后的方程:

(5)

因此,為了求得權(quán)重因子ωi,必須先計(jì)算出下列參數(shù):

Ω=Γ-1Γ0,

(6)

式中：Ω為包含權(quán)重因子ωi和拉格朗日乘子λ的向量；Γ為包含已知量和方位半變異的矩陣；Γ0為包含未知量和方位半變異的矩陣．

2 TEC的二維分布結(jié)果與分析

2.1 數(shù)據(jù)來源

本文使用的衛(wèi)星數(shù)據(jù)全部來自于IGS的全球GPS臺(tái)網(wǎng),從ftp://cddis.gsfc.nasa.gov/pub/gps數(shù)據(jù)中心下載而得,選取觀測數(shù)據(jù)的經(jīng)度和緯度范圍為:90°E～150°E、50°S～50°N．此區(qū)域涵蓋了約60個(gè)觀測臺(tái)站,觀測臺(tái)站接收數(shù)據(jù)的時(shí)間分辨率均為30 s．

2.2 TEC的二維分布與精度驗(yàn)證

由于測量精度、儀器偏差等問題的存在,在解算TEC的過程中,不同臺(tái)站獲得的TEC值可能會(huì)存在一定的誤差,這會(huì)對TEC二維分布重構(gòu)產(chǎn)生較大的影響．為此需要對各臺(tái)站的TEC數(shù)據(jù)進(jìn)行預(yù)處理,具體方法為:首先把同一段時(shí)間內(nèi)(本文取15 min)的垂直TEC值按經(jīng)度進(jìn)行劃分,經(jīng)度間隔取2.5°;然后畫出相同經(jīng)度區(qū)間內(nèi)垂直TEC隨緯度的變化曲線,采用數(shù)據(jù)線性擬合的方法對曲線中相比其它數(shù)據(jù)點(diǎn)偏差過大的點(diǎn)進(jìn)行濾除,以保證所有的垂直TEC曲線都較為平滑．

對所有經(jīng)度區(qū)間的數(shù)據(jù)點(diǎn)處理過后,即可得到所測電離層垂直TEC數(shù)據(jù)的分布圖,如圖1所示．得到垂直TEC值的分布以后便可進(jìn)行插值處理,本文使用Kriging法進(jìn)行插值處理．由于數(shù)據(jù)量過大,計(jì)算過程中會(huì)出現(xiàn)超出計(jì)算機(jī)內(nèi)存的情況,需要對原始數(shù)據(jù)進(jìn)行取中值處理:對于每一對衛(wèi)星接收站而言,在15 min內(nèi)大約可以獲得30個(gè)位于電離層穿刺點(diǎn)處的垂直TEC值,按每十組數(shù)據(jù)對垂直TEC取中值的方法對衛(wèi)星信號(hào)進(jìn)行處理,從而每對衛(wèi)星接收站之間可以獲得的數(shù)據(jù)量可以減少到約3組,這樣可以大大降低運(yùn)算內(nèi)存且能優(yōu)化計(jì)算結(jié)果．經(jīng)過插值后便可以得到TEC的二維分布圖,如圖2所示．可以看出:去掉測量過程中一些偏差較大的數(shù)據(jù)可明顯提高TEC值的重構(gòu)質(zhì)量．

圖1 數(shù)據(jù)處理前后觀測值分布情況

圖2 由Kriging內(nèi)插得到的TEC二維分布圖

在獲得高分辨率的TEC二維分布圖后,采用以下方法對精度進(jìn)行分析:在60個(gè)接收站隨機(jī)選取59個(gè)接收站的數(shù)據(jù)進(jìn)行TEC重構(gòu),剩下1個(gè)站的數(shù)據(jù)作為真實(shí)值對我們的成像結(jié)果進(jìn)行驗(yàn)證,以評估本文方法的TEC重構(gòu)誤差．

設(shè)重構(gòu)得到的TEC值為TECsimulation,TEC真實(shí)值為TECstation,則誤差為

圖3是TEC重構(gòu)誤差分布圖,選取不同日期不同時(shí)段的共611個(gè)TEC值進(jìn)行對比分析,從圖中我們可以看出,經(jīng)過插值得到的TEC地圖的誤差主要集中在10%以下,精度還是比較理想的,驗(yàn)證了本文方法的有效性和可靠性．

圖3 TEC成像結(jié)果的誤差分布圖

2.3 實(shí)測數(shù)據(jù)成像結(jié)果分析

圖4示出了東經(jīng)120°,不同緯度下TEC值在一天內(nèi)隨時(shí)間的變化情況,很明顯地出現(xiàn)了一天中TEC值先升高后降低的變化趨勢,而且緯度越低變化趨勢越明顯;在白天可以很清晰地看出TEC值從赤道向緯度升高的方向變小的趨勢,而晚上各緯度TEC值相差不大;且從圖中可以看出白天TEC值在20°緯度附近有一個(gè)反常變大的過程,這是由于赤道異常區(qū)的影響．一般來講,赤道異常在日出后形成于地方時(shí)10:00并逐漸增強(qiáng),在地方時(shí)14:00左右達(dá)到最大值,然后又逐漸衰退,在地方時(shí)04:00左右有極小值．一天之中,異常峰先在低緯度發(fā)展出現(xiàn),并于增強(qiáng)時(shí)逐漸移向較高緯度,在其之后衰退時(shí),又逐漸返回低緯度,最后消失[24]．圖5所示為2014年02月21日00:00-24:00UT時(shí)刻的TEC二維分布重構(gòu)結(jié)果．在圖中可以很明顯地看出TEC極大值區(qū)域在隨著時(shí)間自東向西偏移,這是因?yàn)榈厍蜃晕飨驏|自轉(zhuǎn)因素的存在;由于白天電子密度會(huì)升高,白天各地區(qū)的TEC值比晚上大很多．本文重構(gòu)的TEC結(jié)果較好反映了電離層的日變化和區(qū)域變化特征．

圖4 垂直TEC隨緯度和時(shí)間的變化結(jié)果

圖5 同一天內(nèi)不同時(shí)刻的TEC二維分布圖

在得到亞大地區(qū)的高精度TEC二維分布圖后,便可以借此分析TEC隨緯度、時(shí)間等的變化情況以及在各種突發(fā)事件中的變化趨勢．相比于磁靜日,磁暴期間電離層會(huì)發(fā)生顯著的變化．如圖6所示,給出了2014年2月25日-3月5日一次磁暴事件期間電離層TEC的變化結(jié)果．從圖中可以看出Dst變化情況和東徑120°不同緯度處TEC變化情況，本次磁暴發(fā)生在2月27日夜晚,最小Dst值為-94 nT;其中黑色實(shí)線表示實(shí)時(shí)的TEC值,紅色虛線表示靜日TEC平均日變化．本文以Dst指數(shù)作為靜日的選取標(biāo)準(zhǔn),期間所有Dst指數(shù)均不小于-20 nT,以此來選擇磁暴前后共7天的TEC值進(jìn)行平均以作為本文的靜日參考．可以看出隨著磁暴的發(fā)生不同緯度的TEC值都出現(xiàn)了不同程度的變化,在磁暴恢復(fù)相期間,各緯度TEC值出現(xiàn)不同程度的減小,呈現(xiàn)負(fù)暴效應(yīng),尤其在低緯赤道異常區(qū)的負(fù)暴效應(yīng)十分明顯,這應(yīng)該是由于磁暴期間出現(xiàn)的擾動(dòng)發(fā)電機(jī)電場阻礙了赤道異常駝峰的形成所造成的[16]．由于磁暴主相持續(xù)時(shí)間較短,因此主相期間TEC的變化情況我們以圖7進(jìn)行詳細(xì)分析．

圖6 一次磁暴期間的Dst指數(shù)以及不同緯度TEC的變化情況

圖7 磁暴主相期間TEC與靜日TEC對比圖

圖7示出了磁暴主相期間TEC值與寧靜日里相對應(yīng)地方時(shí)的TEC值之間的差異,本次磁暴急始發(fā)生于27日17:00UT時(shí),并于27日23:00UT時(shí)達(dá)到主相極大．26日作為我們的靜日參考,從中可以看出在主相期間,大部分地區(qū)TEC值都出現(xiàn)了明顯的增大,整體呈現(xiàn)正暴相的特征．不過由于磁暴期間的物理機(jī)制極其復(fù)雜,不同區(qū)域內(nèi)TEC值的變化情況也有所不同．例如南北半球TEC值的變化差異極大,TEC值在北半球東亞地區(qū)升高幅度較為明顯,最高可達(dá)100%,而在澳大利亞南部區(qū)域則出現(xiàn)了下降的趨勢,下降最大百分比接近60%,這種南北半球的差異化表現(xiàn)可能是由于磁暴帶來的帶電粒子由南向北的遷移現(xiàn)象造成的．

圖8示出了磁暴恢復(fù)相期間TEC值與寧靜日里相對應(yīng)地方時(shí)的TEC值之間的差異,可以看出在恢復(fù)相期間,大部分地區(qū)TEC值都出現(xiàn)了不同程度的下降,呈現(xiàn)負(fù)暴的特征,而且南半球TEC的下降幅度要比北半球大很多,這也符合圖6中顯示的TEC變化趨勢．

圖8 磁暴恢復(fù)相期間與寧靜日TEC對比

由于磁暴期間復(fù)雜的物理機(jī)制,不同經(jīng)緯度地區(qū)的TEC值變化情況也有所不同．例如在北半球中低緯區(qū)域有著明顯的先正暴后負(fù)暴的現(xiàn)象,而在澳大利亞南部區(qū)域卻只有負(fù)暴產(chǎn)生,全程未見正暴的存在．這一南北半球的差異化表現(xiàn)對于深入了解電離層暴的特征及其動(dòng)力學(xué)過程無疑是非常有益的．

圖9示出了2015年與2016年的兩次磁暴期間的電離層TEC變化情況．在第一次磁暴中,磁暴急始發(fā)生于12月20日04:00UT時(shí)左右,并于22:00UT時(shí)達(dá)到主相極大,Dst指數(shù)達(dá)-170 nT．從圖中可以看出,本次磁暴期間大體呈現(xiàn)了先正暴后負(fù)暴的現(xiàn)象,而不同緯度的具體表現(xiàn)也有著細(xì)微的不同．圖10(a)示出了本次磁暴期間TEC與寧靜時(shí)的差值,在主相期間,各緯度的TEC值都有明顯的增大,而且有著增大幅度從赤道向兩極逐漸變強(qiáng)的趨勢;而恢復(fù)相期間不同緯度的TEC變化情況有著較大的差異,赤道與中低緯地區(qū)的TEC值有著明顯的降低,而高緯區(qū)的TEC值繼續(xù)保持著增大的趨勢．

圖9 磁暴期間TEC的變化情況

第二次磁暴發(fā)生于2016年5月7日-2016年5月9日．在本次磁暴中,磁暴急始發(fā)生于8日01：00UT時(shí)左右,并于08：00UT時(shí)達(dá)到主相極大,Dst指數(shù)達(dá)到-92 nT．圖10(b)示出了本次磁暴期間TEC與寧靜時(shí)的差值,在本次磁暴期間不同緯度的TEC變化趨勢有較大的差異,在主相期間,赤道地區(qū)TEC變化程度不太明顯,而高緯地區(qū)出現(xiàn)了明顯的正暴現(xiàn)象;而在恢復(fù)相期間,南北半球出現(xiàn)了較大的差異,北半球中低緯區(qū)域主要出現(xiàn)負(fù)暴,而南半球相應(yīng)緯度出現(xiàn)了明顯的正暴,這與2014年2月27日磁暴期間北半球TEC值升高,南半球TEC值減小的現(xiàn)象正好相反．

由這幾次磁暴事件的對比可以看出,亞大區(qū)域在磁暴期間出現(xiàn)先正暴后負(fù)暴的現(xiàn)象比較普遍,不過不同緯度區(qū)域的具體表現(xiàn)情況有所不同:例如主相期間高緯TEC值的升高程度要普遍大于低緯地區(qū),而且南北半球相同緯度的TEC變化趨勢也有所不同,經(jīng)常會(huì)出現(xiàn)相反的變化趨勢．

(a) 第一次磁暴 (b) 第二次磁暴圖10 磁暴期與寧靜時(shí)的TEC差值

3 結(jié) 論

本文提出了一種高精度TEC二維分布重構(gòu)方法,基于亞大地區(qū)60個(gè)GPS臺(tái)站的實(shí)測數(shù)據(jù),利用Kriging內(nèi)插方法實(shí)現(xiàn)了亞大區(qū)域高分辨率TEC地圖重構(gòu)．基于高分辨率TEC地圖重構(gòu)結(jié)果,本文對亞大區(qū)域幾次磁暴事件期間的電離層TEC變化進(jìn)行了分析,研究結(jié)果如下:

1)通過與實(shí)測結(jié)果的對比可以看出,本文的TEC地圖繪制結(jié)果精度較為理想,經(jīng)計(jì)算誤差普遍在10%以下,與實(shí)測值平均相差2.6 TECU左右,驗(yàn)證了本文方法的有效性．

2)基于TEC的二維分布分析了亞大地區(qū)在不同年份的幾次磁暴事件的影響,分析了不同緯度地區(qū)在磁暴主相與恢復(fù)相期間的差異化表現(xiàn),大部分地區(qū)會(huì)出現(xiàn)先正暴后負(fù)暴的現(xiàn)象,高緯地區(qū)的TEC變化幅度較之低緯地區(qū)要更為明顯,而且可以看出南北半球相同緯度地區(qū)在磁暴期間的表現(xiàn)也存在極大不同,這對進(jìn)一步了解磁暴的物理機(jī)制是極其有意義的．

致謝：本文所用實(shí)測數(shù)據(jù)均下載自IGS組織提供的全球觀測資料,在此感謝他們對電離層觀測研究工作做出的無私貢獻(xiàn)．