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

余侯芳,郭文興,張玉強,馮健,甄衛民

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

0 引 言

電離層總電子含量(TEC)是表征電離層特性的重要特征參量之一,研究其時空變化規律對衛星通信、導航定位、空間天氣預報等領域有著極其重要的理論意義及應用價值[1-2]．

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

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

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

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

雖然在電離層TEC領域已經開展了大量的研究工作,也取得了不錯的研究成果,但隨著科技的進步,研究的發展,必然對電離層TEC的精度和分辨率提出更高的要求．本文提出了一種高精度TEC地圖重構方法,利用亞大地區60個GPS臺站的實測數據,基于Kriging內插方法實現了亞大區域高分辨率TEC地圖重構,通過與實測數據進行對比并對成像精度進行了驗證;在此基礎上,詳細分析了南北半球不同緯度的電離層TEC值在磁暴期間的不同表現．

1 Kriging插值方法

1.1 半變異函數(semi-variogram)

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

(1)

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

一旦求得經驗半變異函數,下一步便是通過滿足各種數學條件轉化為可以應用在Kriging方程中．

1.2 Kriging方法

Kriging是一種屬于最優線性無偏估計的線性內插方法,因此,Kriging方法的主要目的就是利用已知量的線性組合來估計未知變量Z*:

(2)

式中:ωi是權重因子,對應著每個已知量Zi．

(3)

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

(4)

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

(5)

因此,為了求得權重因子ωi,必須先計算出下列參數:

Ω=Γ-1Γ0,

(6)

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

2 TEC的二維分布結果與分析

2.1 數據來源

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

2.2 TEC的二維分布與精度驗證

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

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

圖1 數據處理前后觀測值分布情況

圖2 由Kriging內插得到的TEC二維分布圖

在獲得高分辨率的TEC二維分布圖后,采用以下方法對精度進行分析:在60個接收站隨機選取59個接收站的數據進行TEC重構,剩下1個站的數據作為真實值對我們的成像結果進行驗證,以評估本文方法的TEC重構誤差．

設重構得到的TEC值為TECsimulation,TEC真實值為TECstation,則誤差為

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

圖3 TEC成像結果的誤差分布圖

2.3 實測數據成像結果分析

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

圖4 垂直TEC隨緯度和時間的變化結果

圖5 同一天內不同時刻的TEC二維分布圖

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

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

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

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

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

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

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

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

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

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

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

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

3 結 論

本文提出了一種高精度TEC二維分布重構方法,基于亞大地區60個GPS臺站的實測數據,利用Kriging內插方法實現了亞大區域高分辨率TEC地圖重構．基于高分辨率TEC地圖重構結果,本文對亞大區域幾次磁暴事件期間的電離層TEC變化進行了分析,研究結果如下:

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

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

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