朱佳楠， 郭建廣 ， 倪彬彬， 曹興， 項正，付松， 顧旭東， 馬新

1 武漢大學電子信息學院空間物理系, 武漢 430072 2 中國科學院比較行星學卓越創(chuàng)新中心, 合肥 230026 3 中國氣象局國家空間天氣監(jiān)測預警中心, 北京 100081

0 引言

地球輻射帶由地球磁場捕獲的高能帶電粒子構成，被槽區(qū)分割為內(nèi)輻射帶和外輻射帶.其中內(nèi)輻射帶(距地心1.1～2.5RE)主要由質(zhì)子組成，其通量水平相對穩(wěn)定，外輻射帶(距地心3～8RE)主要由電子組成，其通量呈現(xiàn)著復雜、劇烈的動態(tài)演化特征.由于外輻射帶中高能電子對衛(wèi)星儀器和宇航員健康存在潛在威脅，研究輻射帶電子的分布與演化過程有著十分重要的意義.

輻射帶電子的分布和演化與磁暴發(fā)生關系密切.Reeves等(2003)通過分析1989—2000年期間276個磁暴事件得出，磁暴可以增加或減少輻射帶相對論電子的通量，也會使電子通量保持不變.電子通量的變化過程是電子輸運、加速和損失過程的平衡.由ULF波驅(qū)動的向內(nèi)的徑向擴散被認為是電子加速的重要機制(Ozeke et al., 2012, 2014)，向外的徑向擴散則造成電子損失(Shprits et al., 2008; Bortnik et al., 2006).波粒相互作用引起的能量和投擲角擴散也能夠有效地加速和損失電子(Meredith et al., 2002; Horne et al., 2005; Thorne et al., 2005; Cao et al., 2019; Ma et al., 2020).一系列研究表明，合聲波既可以加速電子也能導致電子的散射損失(Summers et al., 2007a,b; Thorne et al., 2007, 2013)，嘶聲波通過投擲角散射可損失幾十千到幾兆電子伏特的電子(Meredith et al., 2007, 2009; Baker et al., 2007; Ni et al., 2013; Cao et al., 2017, 2020).

研究輻射帶電子相空間密度(Phase Space Density，簡記為PSD)分布可以分析潛在的電子動態(tài)變化機制(Xiang et al., 2017, 2018; 劉澤源等, 2020).但是衛(wèi)星觀測受限于軌道的時空位置和有限的投擲角、能量范圍.模擬輻射帶的物理模型存在計算方程簡化處理、初值和邊界值難以確定等不足.數(shù)據(jù)同化技術將有限的衛(wèi)星觀測與物理模型結合起來，得到輻射帶電子分布的最佳估計.

Naehr和Toffoletto等(2005)首先將基于卡爾曼濾波的數(shù)據(jù)同化方法使用于輻射帶粒子分布的模擬與預測上，驗證了該方法的可行性.Koller等(2007)、Shprits等(2007)、Kondrashov等(2007)的工作進一步表明，利用數(shù)據(jù)同化方法可以準確地再現(xiàn)輻射帶電子PSD分布.后續(xù)的工作豐富了數(shù)據(jù)同化方法，Ni等(2009a)表明同化過程中使用多個衛(wèi)星的觀測數(shù)據(jù)較使用單個衛(wèi)星誤差更小，并發(fā)現(xiàn)同化結果對于磁場模型的選擇相對不敏感(Ni et al., 2009b).Daae等(2011)研究了邊界條件和損耗模型，表明卡爾曼濾波對各種模型參數(shù)都具有魯棒性.使用上述只考慮徑向擴散的同化模型，Shprits(2012)和Ni等(2013)細致分析了輻射帶電子的動態(tài)變化過程.考慮徑向擴散、投擲角擴散、能量擴散的三維同化模型也被運用于輻射帶研究工作中，并在此基礎上提出了集成卡爾曼濾波(Bourdarie and Maget, 2012)、分裂算子卡爾曼濾波等方法(Shprits et al., 2013; Kellerman et al., 2014).數(shù)據(jù)同化中一個重要的研究對象是新息矢量，它是同化過程中觀測對物理模型修改程度的度量.Shprits等(2012)利用新息矢量進行統(tǒng)計分析，發(fā)現(xiàn)槽區(qū)存在局部加速機制.Cervantes等(2020)通過對比各種損失機制下的新息矢量，得出了不同損失機制的作用范圍及其與地磁活動間的關系.

本文首先利用范艾倫A星、B星GOES-13、GOES-15四顆衛(wèi)星的數(shù)據(jù)對2013年3月的輻射帶電子通量進行了同化分析，細致對比了一維、二維、三維不同輻射帶物理模型下的同化結果，并利用三維同化模型重現(xiàn)了2013年一整年電子PSD的長期演化.通過分析同化過程中的新息矢量，本文進一步討論了所使用同化模型的性能.

1 同化方法與數(shù)據(jù)

1.1 衛(wèi)星觀測數(shù)據(jù)

本文使用了范艾倫雙星和GOES-13、GOES-15四顆衛(wèi)星2013年一整年的輻射帶電子觀測數(shù)據(jù).需將衛(wèi)星電子通量數(shù)據(jù)在相空間坐標(μ,K,L*)下轉(zhuǎn)換為PSD作為同化過程的輸入，轉(zhuǎn)化過程借鑒了Ni等(2009a)的方法.并參考Kellerman等(2014)的工作，使用PSD匹配算法對來自不同衛(wèi)星的觀測數(shù)據(jù)進行相互校準.

范艾倫衛(wèi)星包括A星、B星兩顆運行軌道和攜帶儀器完全相同的衛(wèi)星，其軌道傾角約10°，近地點高度約1.1RE，遠地點為約5.9RE，運行軌道周期約9 h.本研究使用范艾倫衛(wèi)星搭載的高能粒子、成分分析與熱等離子體(ECT, Energetic Particle, Composition, and Thermal Plasma)科學探測單元，其中磁化電子離子分光儀(MagEIS, The Magnetic Electron Ion Spectrometer)提供了能量范圍約10 keV～4 MeV的電子通量數(shù)據(jù)，相對論電子質(zhì)子探測器(REPT, Relativistic Electron Proton Telescope)提供了能量范圍約1.8～20 MeV的粒子數(shù)據(jù).衛(wèi)星投擲角分布在均勻網(wǎng)格中插值，步長為5°.GOES系列衛(wèi)星為一系列地球同步軌道衛(wèi)星，用于對近地空間氣象學和空間天氣的監(jiān)測.本研究使用了GOES-13和GOES-15兩顆衛(wèi)星的數(shù)據(jù)，其上搭載的粒子輻射探測儀器EPEAD(Energetic Proton, Electron, and Alpha Detector)可用于探測0.6～4 MeV的電子，MAGED(MAGnetospheric Electron Detector)用于探測30～600 keV的電子.

1.2 輻射帶擴散模型

借鑒Shprits等(2008)和Subbotin(2010)等的工作，本文使用VERB-3D程序，通過求解三維Fokker-Planck擴散方程模擬電子PSD的演化，三維Fokker-Planck擴散方程的表達式為

(1)

其中f表示電子PSD；μ為第一絕熱不變量，與電子的回旋運動有關；J為第二絕熱不變量，與電子的彈跳運動有關；L*與第三絕熱不變量Φ成反比關系，沿粒子漂移路徑保持不變；DL*L*、Dp p、Dα α分別是彈跳平均下徑向、能量和投擲角擴散系數(shù)；T(α0)為與粒子彈跳頻率有關的函數(shù).

求解Fokker-Planck擴散方程需要已知擴散系數(shù).參照Brautigam和Albert(2000)的方法，徑向擴散系數(shù)表示為與Kp指數(shù)相關的函數(shù).本研究只考慮嘶聲波和合聲波引起的電子散射，使用FDC程序(Full Diffusion Code)計算得到(Ni et al., 2008, 2015; Shprits et al., 2009; Cao et al., 2016, 2020)電子投擲角擴散和能量擴散系數(shù)，用以量化磁層波動對輻射帶電子的影響.

1.3 卡爾曼濾波和新息矢量

卡爾曼濾波是一種對線性系統(tǒng)狀態(tài)量進行最小方差無偏估計的算法，被用于數(shù)據(jù)同化過程.假設模型和觀測誤差是無偏的，且服從高斯分布.在更新時刻，卡爾曼濾波對當前時刻的觀測值和擴散模型的預測值進行比較，在知道前一時刻分析誤差和前一時刻模型誤差的情況下最優(yōu)地估計出當前時刻的預測誤差，從而對預測值進行調(diào)整.卡爾曼濾波的計算方法如下：

(2)

(3)

(4)

國際文憑(International Baccalaureate， IB)是由國際文憑組織(International Baccalaureate Organization， IBO)推出的課程和考試項目，大學預科國際文憑課程(Diploma Program， DP)則是面向全球16～19歲優(yōu)秀高中生開展的課程[1]。IBDP生物學課程的總體特色是“輕知識、重技能”，強調(diào)實驗教學的過程性和方法，旨在關注學生知識的生成過程，培養(yǎng)學生解決問題的能力，全方面提高生物學素養(yǎng)，在情感教育上則充分體現(xiàn)了以人為本、尊重生命的課程目標。

(5)

與三維Fokker-Planck擴散方程相對應，本文將考慮一種、兩種及三種擴散機制的輻射帶物理模型描述為一維、二維、三維擴散模型.三組輻射帶擴散模型的計算結果被用做本文同化過程的輸入，分別為一維純徑向擴散模型，二維徑向擴散加投擲角擴散模型，三維徑向擴散、投擲角擴散加能量擴散模型.經(jīng)過卡爾曼濾波的同化處理，分別對應了下文所述的一維、二維、三維同化模型.

2 同化結果與分析

2.1 2013年3月不同擴散模型下輻射帶電子通量的同化結果

采用上述數(shù)據(jù)同化的方法，本文首先計算了2013年3月輻射帶電子通量的觀測及多維同化結果，如圖1所示.同化過程中設定的觀測誤差和模型誤差均為0.5，即觀測誤差與模型誤差的比值為100%.圖1a展示了在L*=3～7范圍內(nèi)Ek=0.6 MeV，αeq=30°電子通量的衛(wèi)星觀測結果.圖1b—d展示了該能量和投擲角下一維、二維、三維不同擴散模型的同化結果.與上述類似，圖1e—h展示了Ek=1 MeV，αeq=55°電子通量的衛(wèi)星觀測及同化結果.圖1i給出了該時間段內(nèi)Kp和Dst指數(shù)的變化情況.圖中的黑線表示由Carpenter和Anderson(1992)模型計算得到的等離子體層頂位置.如Kp指數(shù)所示，3月1日和3月17日分別發(fā)生了兩次磁暴，其中3月17日的磁暴較強.圖1a的觀測結果顯示，磁暴主相期間輻射帶電子通量快速損失，在恢復相通量有所增強并形成峰值.圖1b—d所示三組不同擴散模型下的同化結果均模擬出了磁暴前后電子通量的變化過程，且填補了L*=5.5附近和L*>6.5時觀測值的空缺，較為精確地重現(xiàn)了輻射帶電子通量的徑向分布.圖1e—h所示Ek=1 MeV，αeq=60°電子的通量值略小于Ek=0.6 MeV，αeq=30°的電子通量值，同化模型在該能量和投擲角下也可以模擬出電子通量的分布和變化.

圖1 設定觀測誤差和模型誤差均為0.5時，2013年3月輻射帶電子通量的觀測結果和一維、二維、三維不同輻射帶擴散模型下隨時間和L*分布變化的同化結果.(a—d) 對應Ek=0.6 MeV，αeq=30°的輻射帶電子； (e—h) 對應Ek=1 MeV，αeq=55°的輻射帶電子； (i) Kp和Dst地磁指數(shù)隨時間的變化.(a—h)中黑線表示等離子體層頂位置Fig.1 For the case that both the observation error and model error are 100%, observations of radiation belt electron fluxes and the corresponding data assimilation results using the 1-D, 2-D and 3-D radiation belt diffusion models over the month of March 2013 for radiation belt electrons with (a—d) Ek= 0.6 MeV, αeq=30° and (e—h) Ek=1 MeV, αeq=55°. (i) The time series of Kp and Dst indices. Black curves in (a—h) show the plasmapause location

為進一步展示采用不同維度擴散模型對同化結果的影響，我們在圖1結果的基礎上重新設置觀測誤差為1，模型誤差為0.5，即觀測誤差與模型誤差的比值為200%，其余參數(shù)與圖1相同，結果如圖2所示.需要說明的是，觀測誤差與模型誤差的比值越大，同化結果越忽視觀測對擴散模型的修改，而更關注擴散模型本身的模擬效果.由圖2所示同化結果可以看出，磁暴發(fā)生前后電子通量的變化趨勢與圖1基本一致，但三組擴散模型的模擬效果表現(xiàn)出明顯的差異.與觀測相比，圖2b和圖2f所示的一維同化結果在磁暴主相期間高L*區(qū)域的電子通量值遠小于觀測值，在磁平靜期低L*的電子通量遠高于觀測值.二維擴散模型在純徑向擴散的基礎上加入了投擲角擴散，可以從圖2c和和2g所示的同化結果觀察到低L*電子的顯著損失，即嘶聲波可以對等離子體層頂以內(nèi)的輻射帶電子有效散射.圖2d和圖2h在二維擴散模型的基礎上加入了能量擴散，可以在L*>4的區(qū)域觀察到磁暴后電子通量的迅速增加，即在等離子體層頂以外合聲波能對輻射帶電子有效加速.對比三組維度擴散模型下的同化結果，三維同化模型提供了有效的加速和損失機制，其同化結果能最優(yōu)地模擬出輻射帶電子通量的徑向分布.

圖2 同圖1類似，但是設定觀測誤差為1，模型誤差為0.5Fig.2 Same as in Figure 1, except for the setting that the observation error is 1 and the model error is 0.5

對應于圖1和圖2，圖3展示了不同維度擴散模型下的同化結果在不同L*的模擬效果.從左至右對應兩組能量和投擲角，從上至下對應兩組觀測和模型誤差，子圖自上而下分別為L*= 4、5、6時電子通量的觀測和三組同化結果隨時間的變化.與圖1、圖2的結論一致，L*= 4時，一維同化結果在3月1日磁暴后的磁平靜期高于觀測結果，L*= 5、6時，一維、二維同化結果在磁暴后均低于觀測結果.三維同化結果在兩組能量、投擲角和三組L*下均與觀測結果很好地吻合，有效修正了一維和二維同化結果的誤差，說明除了徑向擴散和投擲角擴散，能量擴散對輻射帶電子演化也發(fā)揮著重要作用.還可以看出，三組同化結果與觀測的吻合度隨時間的推移變高，由于同化過程會在每個時間步長更新誤差矩陣，同化的時間越長模型的同化效果越好.

圖3 對應于圖1和圖2，當L*=4、5、6時，輻射帶電子通量的觀測和不同維度輻射帶擴散模型下的同化結果的比較Fig.3 Corresponding to Figures 1 and 2, comparisons between the observations of radiation belt electron fluxes and the assimilation results in March 2013 for the indicated two at L*=4, 5 and 6

圖4 偶極地磁場下絕熱不變量K、μ與赤道投擲角αeq、電子能量Ek隨L*的對應關系Fig.4 Dependence of equatorial pitch angle and electron kinetic energy on L* in a dipole geomagnetic field for the indicated four pairs of (K, μ)

增大觀測誤差與模型誤差的比值可以展示不同維度擴散模型下同化結果的顯著差異，減小觀測誤差與模型誤差的比值可以得到與觀測更為吻合的同化模型.通過一系列測試，當觀測誤差與模型誤差均為0.5時，三維同化結果已經(jīng)可以精確地重現(xiàn)輻射帶電子的分布和變化.

2.2 2013年全年輻射帶電子PSD的同化結果分析

為了展示三維輻射帶同化模型對較長時間輻射帶電子動態(tài)演變的模擬效果，本文對2013年一整年輻射帶電子PSD進行了同化分析.通過研究以三個絕熱不變量(μ,K,L*)構成相空間坐標系下電子PSD的長期演化，可以區(qū)分非絕熱效應和絕熱效應、局部加速和徑向傳輸過程，有助于理解輻射帶電子加速、輸運和損耗過程之間的平衡.如章節(jié)1.2所述，μ是第一絕熱不變量；第二絕熱不變量J通常用K來代替，K只與背景磁場的強度和結構有關，與帶電粒子質(zhì)量和電荷數(shù)無關.我們選取兩組固定K值下不同μ值的四個有代表性的絕熱不變量結果進行分析，其中第一組為K=0.5 G0.5RE，μ=100 MeV/G和μ=400 MeV/G，第二組為K=0.11 G0.5RE，μ=300 MeV/G和μ=1300 MeV/G，偶極磁場下它們與αeq和Ek隨L*的對應關系如圖4所示.通過選取這四組絕熱不變量，可以獲得在較高和較低投擲角、較大能量范圍下的電子PSD.

圖5 設定K=0.5 G0.5RE，2013年全年輻射帶電子相空間密度隨L*和時間分布的觀測和三維同化結果(a—b) 對應μ=100 MeV/G的輻射帶電子； (c—d) 對應μ=400 MeV/G的輻射帶電子； (e) 地磁指數(shù)Kp和Dst.(a—d)中黑線表示等離子體層頂位置.Fig.5 For the case that K=0.5 G0.5RE, observations of radiation belt electron PSD and the corresponding data assimilation results using 3-D radiation belt diffusion models during 2013 with (a—b) μ=100 MeV/G and (c—d) μ=400 MeV/G. (e) Kp and Dst indices. Black curves in (a—d) show the plasmapause location

同圖5類似，圖6顯示了設定K=0.11 G0.5RE，μ=300 MeV/G和μ=1300 MeV/G時輻射帶電子PSD的觀測和同化結果.電子的投擲角更高，其PSD值也更高，但PSD的變化趨勢基本不變.在該組絕熱不變量下，三維同化結果仍可以很好地模擬出電子PSD的徑向分布和長期演化過程.四組絕熱不變量下的同化結果表明，磁暴造成的電子PSD損失與增加可以發(fā)生在廣泛的能量和投擲角范圍中.

圖6 同圖5類似，但設置K=0.11 G0.5RE，μ=300 MeV/G和μ=1300 MeV/GFig.6 Same as in Fig.5, except for the setting that K=0.11 G0.5RE, μ=300 MeV/G and μ=1300 MeV/G

本文計算了2013年一整年的平均新息矢量，結果如圖7所示：從左至右分別代表三個不同維度的同化模型，從上至下對應不同的μ、K.如章節(jié)1.3所述，新息矢量反映了同化過程在多大程度上修改了模擬值，新息矢量為正(負)值表示衛(wèi)星觀測值大于(小于)擴散模型模擬值，即模擬值與觀測值相比缺少了一些加速(損失)機制，同化過程通過從模擬值中加入(減去)附加PSD值修正模擬值與觀測值之間的誤差.新息矢量還可以作為模型中電子損耗與源的指標，當模型中加入不同擴散機制后，同一絕熱不變量條件下的新息矢量由正值變?yōu)樨撝担砻髂Ｐ图尤肓擞行У募铀贆C制；新息矢量由負值變?yōu)檎担砻骷尤肓擞行У膿p失機制.

從圖7可以看出，一維同化結果的新息矢量在L*>5.8的區(qū)域為正，在低L*趨于零，同化過程未對模型值有效修正.比較一維和二維的結果，L*=4～5的新息矢量由趨于零變?yōu)檎Ｐ椭屑尤肓藫p失機制，表明投擲角擴散對電子有效損耗；比較二維和三維的結果，L*>5.2的新息矢量由正變?yōu)樨?，模型中加入了加速機制，表明能量擴散對電子有效加速.三維同化結果的新息矢量在L*=4～5為正，即觀測值高于模擬值，因此同化過程在模擬值上加入了一定PSD；新息矢量在L*>5.2為負，觀測值低于模擬值，因此同化過程從模擬值中減去了一定PSD.觀測值與模型之間的誤差表明本文數(shù)值模擬階段使用的波動經(jīng)驗模型可能不夠精確，加入了過多的損失和加速，也可能是模擬過程中忽略了一些物理機制.高L*處三維同化結果的新息矢量值隨著Kp指數(shù)增大而增大，說明同化過程在地磁活動劇烈時對模型值的修正更多，此時波動經(jīng)驗模型的模擬效果較差.

圖7 2013年一整年的平均新息矢量，橫坐標為L*，縱坐標為Kp指數(shù).從左至右：一維、二維、三維擴散模型.從上至下：K=0.5 G0.5RE，μ=200 MeV/G和400 MeV/G；K=0.11 G0.5RE，μ=300 MeV/G和1300 MeV/G.Fig.7 For the indicated adiabatic invariants (K=0.5 G0.5RE, μ=200 MeV/G and 400 MeV/G; K=0.11 G0.5RE, μ=300 MeV/G and 1300 MeV/G), average innovation vector as a function of L* and Kp using the 1-D, 2-D and 3-D radiation belt diffusion models.

3 總結與討論

本文使用了基于卡爾曼濾波的數(shù)據(jù)同化方法，將范艾倫A、B星和GOES-13、GOES-15四顆衛(wèi)星的觀測數(shù)據(jù)與不同維度的輻射帶物理模型相結合，重構了2013年外輻射帶電子的徑向分布和時空演化過程，得到以下主要結論：

(1) 與一維純徑向擴散數(shù)值模型和二維徑向擴散加投擲角擴散數(shù)值模型相比，加入了能量擴散的三維數(shù)值模型下的數(shù)據(jù)同化結果與衛(wèi)星觀測結果的吻合度最好，這說明在外輻射帶電子動力學變化過程中，除了徑向擴散和投擲角擴散，能量擴散對電子的加速過程也發(fā)揮著重要作用.

(2) 通過對不同絕熱不變量組合的同化模擬，發(fā)現(xiàn)三維輻射帶同化模型可以合理地獲取L*=3～7范圍內(nèi)的地球外輻射帶不同投擲角、不同能量上的電子相空間密度分布，能填補現(xiàn)有觀測數(shù)據(jù)在空間分布上的不足，提供電子相空間的完整徑向分布.

(3)隨著地磁活動的變化，三維輻射帶同化模型較好地重構了外輻射帶電子長期的動態(tài)演化過程.同化結果表明，電子相空間密度的時空變化與磁暴事件具有很強的依賴性：磁暴主相期間會經(jīng)常發(fā)生電子相空間密度在大多數(shù)L*上的快速降低，而在磁暴恢復相期間電子相空間密度一般會增強并達到峰值.

(4)通過分析同化過程中的新息矢量，可以看出輻射帶物理模型常在低L*處低估觀測值而在高L*處高估觀測值.這些差異可能由物理模型中考慮的波動模型不完善造成，也可能由模擬過程中一些物理機制的缺失引起.新息矢量的數(shù)值隨地磁活動增強而增大，說明本文采用的輻射帶物理模型在磁擾期的模擬效果與地磁平靜期相比不確定性增加，有待進一步深入研究.

綜上，本文開展的地球電子輻射帶三維數(shù)據(jù)同化建模具有有效、合理重構外輻射帶電子的徑向分布以及長期時空演化過程的能力.我們會在后續(xù)的研究中依據(jù)新息矢量提供的重要信息對模型加以改進，包括建立更加真實合理的磁層波動全球分布模型、考慮波動對輻射帶電子的混合擴散機制以及加入更多衛(wèi)星數(shù)據(jù)到同化過程.我們希望通過多維度數(shù)據(jù)同化的方法獲取關于地球輻射帶動態(tài)變化的更多詳細信息，進一步加深對近地空間粒子輻射環(huán)境動力學過程與機制的理解認知，服務于空間天氣過程與效應的預測預報.

致謝感謝范艾倫衛(wèi)星團隊提供數(shù)據(jù).范艾倫衛(wèi)星電子通量數(shù)據(jù)來源于(https:∥www. rbsp-ect. lanl.gov/data_pub).地磁活動指數(shù)來自NASA OMNIWEB (https:∥omniweb. gsfc. nasa. gov).感謝德國地學研究中心(GFZ)Yuri Shprits教授給予的幫助.