敖先志，劉四清,2，沈 華，王晶晶，胡駿翔，李 剛

（1.中國科學(xué)院國家空間科學(xué)中心，北京 100190；2.中國科學(xué)院大學(xué)，北京 100190；3.美國阿拉巴馬大學(xué)漢茨維爾分校空間科學(xué)系，漢茨維爾 35899）

引 言

太陽高能粒子事件是影響近地空間（1 AU）以及深空輻射環(huán)境重要的甚至是主要的因素之一，對航天器、人造地球衛(wèi)星，以及宇航員的安全等構(gòu)成巨大的威脅。大規(guī)模的太陽高能粒子事件會造成嚴(yán)重的空間災(zāi)害天氣。因此，太陽高能粒子事件的觀測、產(chǎn)生形成機制以及相應(yīng)的預(yù)報模式是當(dāng)前國際空間物理學(xué)界的最前沿課題之一。

太陽高能粒子事件往往可以分為兩種類型：“脈沖”（impulsive）型和“漸進(jìn)”（gradual）型[1]?！懊}沖”型事件中的高能粒子表現(xiàn)在觀測上的特征是儀器測量到的高能粒子通量有個突然的、急劇的增強，然后又快速下降，持續(xù)時間常常在一個或者數(shù)個小時左右，就好像一個單脈沖一樣；“漸進(jìn)型”事件中的高能粒子其典型觀測特征是升高后的粒子通量緩慢、漸進(jìn)地下降，事件持續(xù)的時間很長，可以達(dá)到2 天以上。這兩種形態(tài)的事件有時候會混合在一起[2-4]。在實際觀測中，“漸進(jìn)”型的太陽高能粒子事件占大多數(shù)[5]。

國外對于太陽高能粒子事件的觀測和理論研究開展得比較早[1,6]。一般認(rèn)為，“漸進(jìn)型”事件中的高能粒子（主要成分為質(zhì)子）是由日冕物質(zhì)拋射驅(qū)動（CME-driven）的激波加速產(chǎn)生[7-11]。這種粒子加速機制通常被稱為激波擴散加速（Diffusive Shock Accelera‐tion，DSA）機制，其本質(zhì)是一階費米加速。這種加速機制可以輕而易舉地把粒子加速到＞100 MeV/nucleon。2000年，Zank 等提出了基于DSA 加速機制的太陽高能粒子加速和傳播模型：PATH（Particle Acceleration and Transport in the Heliosphere）模型。這個模型在過去的10多年里得到了不斷的發(fā)展和完善[12-20]。胡駿翔等進(jìn)一步擴展了前人的研究[21-22]，建立了iPATH（improved PATH）模型，將激波的演化從一維擴展到了二維，從而能夠考慮方位角和垂直擴散系數(shù)對粒子加速和傳播的影響。

本文利用iPATH模型對發(fā)生于2014年04月18日的CME 所引起的太陽高能粒子事件實例進(jìn)行了數(shù)值模擬，并將地球附近（1 AU）的數(shù)值模擬結(jié)果和衛(wèi)星觀測進(jìn)行了對比；在此基礎(chǔ)上加以擴展，進(jìn)一步給出了深空環(huán)境中不同位置處太陽高能粒子通量的模型計算結(jié)果。

1 iPATH模型簡述

iPATH模型是在PATH模型的基礎(chǔ)上進(jìn)行的擴展。iPATH模型的運行主要包括3個步驟：①生成背景太陽風(fēng)環(huán)境；②模擬CME 驅(qū)動的激波并計算粒子在激波附近的加速過程；③計算被激波加速的高能粒子在行星際空間中的傳播。

iPATH對背景太陽風(fēng)的模擬使用了開源數(shù)值模擬程序包ZEUS3D[23]來求解磁流體力學(xué)方程組為

其中：ρ是流體的密度；V是流體的速度；p是流體的熱壓強；B是磁感應(yīng)強度；e是單位體積內(nèi)的能量，其定義包括內(nèi)能和磁能；Φ是引力勢能。

盡管ZEUS3D 具有求解三維磁流體力學(xué)（Mag‐netoHydro Dynamics，MHD）方程組的能力，但是iPATH在當(dāng)前只具有模擬二維條件下的粒子加速和傳播能力，因而在利用ZEUS3D 程序包進(jìn)行計算時，iPATH模型將其中一個維度的大小設(shè)為1。具體來說，ZEUS3D 程序包的設(shè)置采用了球極坐標(biāo)（r，θ，φ）參考系，r，θ，φ分別為徑向距離、仰角和方位角，太陽位于中心，θ=90°的平面是黃道面，iPATH模型的θ角取值被限定于只取90°一個值，從而使得三維模擬程序退化為二維數(shù)值模擬。這種方法的不足在于二維數(shù)值模擬沒有三維數(shù)值模擬更準(zhǔn)確，但是對于太陽高能粒子事件模擬來說可以減少計算資源、大幅提高計算速度，尤其是對于空間天氣預(yù)警和預(yù)報而言，計算時間的減少，帶來的是預(yù)警提前量的提升。

iPATH模型通過在內(nèi)邊界處引入持續(xù)一段時間的擾動來模擬日冕物質(zhì)拋射。內(nèi)邊界處的擾動參數(shù)包括CME 速度、密度、溫度、磁場以及擾動持續(xù)時間，CME的角寬度也是參數(shù)之一，用于描述CME的擾動速度分布的范圍。為簡單起見，iPATH 模型假設(shè)CME的速度分布為高斯分布。

對于發(fā)生于2014年4月18日的CME事件，本文所選擇的初始擾動參數(shù)，盡量使得模擬結(jié)果符合日冕儀和1 AU 處的太陽風(fēng)的實際觀測數(shù)據(jù)。與背景太陽風(fēng)數(shù)值模擬類似，iPATH 模型采用MHD 方程組來描述擾動在太陽風(fēng)中的傳播。隨著CME 擾動的引入，CME驅(qū)動的激波開始在模型中傳播和演化。ZEUS3D程序包采用了交錯網(wǎng)格和迎風(fēng)格式。模型中的初始磁場為簡單的Parker螺旋場。

iPATH模型在計算激波參數(shù)時同步計算粒子的加速過程。對于粒子在激波附近經(jīng)歷加速過程的理論和數(shù)值模型，文獻(xiàn)[12]～[22]有非常詳盡的描述，本文不再重復(fù)敘述。文獻(xiàn)[12]～[20]所發(fā)展的PATH模型在后期已經(jīng)在一定程度上考慮了垂直擴散系數(shù)，但是鑒于PATH模型在處理激波時，采用的是一維模型，因此PATH 無法解決擴散系數(shù)在φ方向上的變化。iPA‐TH 模型從一維擴展為二維空間，擴散系數(shù)進(jìn)而可以表現(xiàn)為φ的函數(shù)，或者更為明確地說，激波波前不同位置具有不同的擴散系數(shù)。

高能粒子的注入率對被加速粒子能譜的形成具有重要作用。注入率受到局地太陽風(fēng)參數(shù)、湍流強度、激波參數(shù)、粒子能量動量等各種因素影響，也是目前高能粒子加速機制研究的一個重要對象?，F(xiàn)階段，在iPATH模型里粒子的注入率是高能粒子能量動量的一個函數(shù)。此外，iPATH模型假設(shè)：激波對粒子有加速作用，而粒子對于激波的反作用可以忽略。iPATH模型利用磁流體力學(xué)方程組來描述和模擬激波傳播，然而磁流體力學(xué)是沒有辦法來刻畫高能粒子加速過程對激波反作用的動力學(xué)特征的。這一點，有待于未來的研究工作來解決。

被加速的高能粒子從激波附近逃逸并沿行星際空間磁力線傳播，傳播過程中被傳播路徑上的湍流所散射。粒子的傳播過程服從聚焦傳輸（focused trans‐port）方程，在iPATH 模型中利用蒙特卡洛（Monte Carlo）方法來求解粒子的傳播過程[21]。

2 數(shù)值模擬結(jié)果

數(shù)個 GOES（Geostationary Operational Environ‐mental Satellites）衛(wèi)星從 2014年04月18日—2014年04月22日均觀測到了強烈的太陽高能質(zhì)子增強事件。這次太陽高能粒子事件來源于一次大的日冕物質(zhì)拋射事件所驅(qū)動的激波。在2014年04月18日，太陽活動區(qū)AR12036 爆發(fā)了一次M7.3 級的耀斑，X 射線的流量在當(dāng)天13：03UT左右達(dá)到了最大值。耀斑爆發(fā)大約20 min 后，位于日地引力平衡點L1 處的太陽和日球?qū)犹綔y器（Solar and Heliospheric Observatory，SO‐HO）衛(wèi)星所搭載的C2 和C3日冕儀觀測到了高速全暈CME。先進(jìn)的成分探測器（Advance Composition Explorer，ACE）衛(wèi)星在 20日大約 10：20UT 左右觀測到了這一次CME爆發(fā)所驅(qū)動的激波。這次全暈CME事件所引發(fā)的太陽高能粒子事件是一起孤立事件，事件前后高能粒子的通量均無由于其它CME 事件所帶來的明顯擾動。GOES-15 衛(wèi)星的觀測表明在事件之前7 天和之后9 天的時間內(nèi)，能量大于2.5 MeV 的質(zhì)子通量保持在背景大小。因此，這次事件非常適合進(jìn)行數(shù)值模擬研究。實際上，在此次CME 爆發(fā)的同一天大約07：24UT 左右爆發(fā)了一個小的CME，其速度和角寬度經(jīng)估算分別為大約387 km/s 和84°。同時期ACE 衛(wèi)星觀測到的太陽風(fēng)速度逐漸從大約400 km/s增加到700 km/s，因而這一次小的CME 事件不大可能會引起SEP事件。另一方面，STEREO的2顆衛(wèi)星也觀測到了這次小的CME 事件。STEREO-A 與地球的夾角為165°，STEREO-B與地球的夾角為156°，這兩顆衛(wèi)星上搭載的日冕儀觀測圖像結(jié)合SOHO衛(wèi)星的日冕儀圖像，我們推斷07：24UT 左右爆發(fā)的小型CME 事件的拋射方向遠(yuǎn)離黃道面，與日地連線接近垂直。故而，這次事件不會對本文的數(shù)值模擬造成非常顯著的干擾。

此次事件數(shù)值模擬計算域的物理尺度為

其中：r的單位為AU；φ的單位為角度。

網(wǎng)格大小設(shè)置為黃道面上的1 500×360 個格點。在r方向和φ方向上均為均勻網(wǎng)格劃分，因此計算域的內(nèi)邊界是一個半徑為0.1 AU 的圓。在計算背景太陽風(fēng)參數(shù)這一步中，為簡化起見，采用了軸對稱設(shè)置，即內(nèi)邊界所有網(wǎng)格點上同一個物理量的數(shù)值是相同。數(shù)值模擬所使用的內(nèi)邊界背景太陽風(fēng)參數(shù)見表1所示。圖1給出了背景太陽風(fēng)的模擬結(jié)果。圖中橫坐標(biāo)是離太陽的距離（以AU 為單位），（a）是歸一化的太陽風(fēng)數(shù)密度，即質(zhì)子數(shù)密度乘以距離的平方；（b）是太陽風(fēng)的速度。在1 AU 處，數(shù)值模擬的太陽風(fēng)數(shù)密度大約為5.9/cm3，速度大約為468 km/s。同時期的ACE衛(wèi)星在1 AU的觀測數(shù)據(jù)表明，太陽風(fēng)數(shù)密度圍繞6/cm3上下波動，速度圍繞500 km/s上下波動。因此，采用表1中的內(nèi)邊界條件數(shù)值模擬得到的結(jié)果是符合實際太陽風(fēng)狀態(tài)的。

表1 數(shù)值模擬內(nèi)邊界參數(shù)Table1 Inner boundary conditions

圖1 背景太陽風(fēng)參數(shù)數(shù)值模擬結(jié)果。上：歸一化的太陽風(fēng)質(zhì)子數(shù)密度；下：太陽風(fēng)速度。橫坐標(biāo)是r方向上的距離Fig.1 Numerical simulation of the background solar wind conditions

中國科學(xué)院國家空間科學(xué)中心所屬空間環(huán)境預(yù)報與研究中心（Space Environment Prediction Center，SEPC）部署了一套日冕物質(zhì)拋射自動識別軟件程序。這套軟件利用Hough 變換和J-maps[24-25]方法從SOHO衛(wèi)星日冕儀的白光日冕圖像中自動檢測和識別CME，并進(jìn)一步采用冰激凌錐模型[26]來反演CME 的速度和角寬度等參數(shù)。在某些情況下，例如全暈CME，自動識別和反演程序可能會把1個CME識別成2個或者多個CME，SEPC 的CME 識別軟件同時提供人工干預(yù)接口以除去明顯不合理的結(jié)果。

圖2是計算機程序自動識別和人工干預(yù)識別對2014年04月18日的CME爆發(fā)進(jìn)行識別的對比圖。圖2（a）是自動識別的結(jié)果，圖2（b）是人工干預(yù)識別的結(jié)果。自動識別將此次全暈CME 識別為2 個CME，本文在識別CME時進(jìn)行了人工修正。將CME識別結(jié)果利用冰激凌錐模型進(jìn)行反演，結(jié)果見圖3。圖3（a）是將圖2（b）的結(jié)果輸入冰激凌錐模型反演的結(jié)果，圖3（b）是全自動反演的結(jié)果。顯然，圖3（b）結(jié)果不是真實的。因此，本文在后續(xù)的數(shù)值模擬中采用圖3（a）的結(jié)果：速度1 290 km/s，角寬度120°。此外，數(shù)值模擬所采用的CME 持續(xù)時間為1.5 h，其它的CME參數(shù)見表1。

圖2 SOHO衛(wèi)星白光日冕圖像CME識別與提取Fig.2 CME detection and recognition based on SOHO coronagraph images

圖3 冰激凌錐模型反演CME參數(shù)Fig.3 CME parameters derived by Cone model fitting method

圖4所示為CME/激波在黃道面內(nèi)傳播的數(shù)值模擬結(jié)果。每個小圖中，太陽位于中心，以白色表示；綠色代表地球所在的位置，藍(lán)色、紅色和黃色代表另外3個假想的觀測點，其中紅色和黃色觀測點的方位角分別為30°和5°，地球和藍(lán)色觀測點的方位角為90°；黑色曲線代表磁力線；顏色棒表示的是速度的大小，單位為km/s。觀測點和激波波前的磁力線連接位置（Connecting with the OBserver Point，COB‐Point）在太陽高能粒子事件中起著非常重要的作用。這幾個觀測點在研究磁力線聯(lián)通與SEP事件之間的關(guān)系具有相當(dāng)?shù)拇硇浴＠?，紅色觀測點在爆發(fā)開始時，磁力線連接在激波最強的中心區(qū)域，隨著時間的推移，連接點（COBPoint）逐漸向激波右側(cè)邊緣移動；黃色觀測點的磁力線連接點（COBPoint）起初位于激波左側(cè)邊緣，隨激波傳播逐漸向激波中心區(qū)域移動；藍(lán)色觀測點在開始時刻幾乎不與激波波前具有磁力線連接。圖3中徑向距離的單位為AU，黑色實線圓為1 AU 的地球軌道。圖4中可以直觀地看到，激波法向和磁力線的夾角隨激波波前位置不同而變化。圖4中每一個小圖，激波中間部分壓縮比大，向兩側(cè)，壓縮比逐漸減小。因此，不同位置處被加速的粒子也具有不一樣的能譜。

高能粒子的注入率在SEP事件數(shù)值模擬中扮演著非常重要的作用，SEP 事件中的粒子通量均與之相關(guān)。假設(shè)粒子的注入率為[20]

圖4 CME/激波在黃道面內(nèi)傳播Fig.4 The simulated propagation of the CME/shock in the ecliptic plane

其中：θ為激波法向和磁力線的夾角；E0是質(zhì)子在平行激波條件下的注入能；是質(zhì)子在斜激波條件下的注入能，它是θ的函數(shù)，( -δ)是種子粒子的能譜分布函數(shù)的冪指數(shù)。

在深空環(huán)境下，很少有機會能夠直接將觀測結(jié)果和模擬結(jié)果進(jìn)行對比。本文利用地球軌道的觀測結(jié)果來確定數(shù)值模型的參數(shù)，再以同樣的參數(shù)來模擬深空環(huán)境下的SEP事件。如前所述，假設(shè)粒子注入率中的χ= 2%，數(shù)值模擬所得的2014年04月18日SEP事件在1 AU地球附近的高能質(zhì)子的單向積分通量見圖5。

圖6是數(shù)值模擬所得到的地球附近的高能質(zhì)子的單向微分通量，圖6給出了6個不同能段粒子的單向微分通量?？偠灾?，對于2014年04月18日的SEP事件在地球附近的表現(xiàn)所進(jìn)行的數(shù)值模擬基本上能反映實際的觀測結(jié)果，本文利用地球觀測的結(jié)果所確定的數(shù)值模擬參數(shù)，將模擬范圍擴大到1 AU 以外，進(jìn)一步對深空環(huán)境下的SEP事件進(jìn)行模擬。

圖5中虛線是地球同步軌道GOES衛(wèi)星的觀測結(jié)果，實線是iPATH模型的模擬結(jié)果。紅色、藍(lán)色和綠色曲線分別代表能量大于10、50、100 MeV 的高能質(zhì)子的單向積分通量。以大于10 MeV的高能質(zhì)子為例，數(shù)值模擬的結(jié)果在事件發(fā)生后的前12 h內(nèi)和觀測較為吻合，之后我們模擬的積分通量開始下降，與之相比觀測值仍然維持了很長時間才開始下降。通過與GOES衛(wèi)星觀測的高能質(zhì)子的單向微分通量以及軟X射線的觀測進(jìn)行對比發(fā)現(xiàn)，在事件發(fā)生的大約前12～24 h，數(shù)值模擬的結(jié)果相對于觀測結(jié)果偏小，其原因可能與行星際太陽風(fēng)結(jié)構(gòu)有關(guān)。在大約27 h后，太陽有一次小的耀斑爆發(fā)，耀斑爆發(fā)產(chǎn)生的入射粒子可能作為種子粒子注入激波附近進(jìn)一步增強了這次SEP事件中＞10 MeV 能量粒子的通量，但是其對更高能的粒子＞50 MeV 的通量影響不大。事件發(fā)生后的前12 h內(nèi)，對于＞50 MeV和＞100 MeV的高能粒子，本文數(shù)值模擬的結(jié)果偏大。

圖5 地球附近SEP事件單向積分通量觀測和模擬對比虛線為同步軌道GOES衛(wèi)星觀測值，實線為數(shù)值模擬結(jié)果Fig.5 Comparison between the observed unidirectional integral flux（dashed lines）by GOES satellite and the simulated results（solid lines）

圖6 地球附近SEP事件單向微分通量的數(shù)值模擬結(jié)果Fig.6 The simulated unidirectional differential flux at near-Earth space

圖7所示為1 AU距離30°處的觀測點高能質(zhì)子單向積分通量的數(shù)值模擬結(jié)果。由于該位置處的磁力線在CME 爆發(fā)開始時刻就連接于激波波前中心區(qū)域，因而高能粒子的通量迅速達(dá)到最高值，隨著激波快速向外傳播，磁力線在激波波前的連接點（COBPoint）也迅速地移動到激波波前邊緣，在事件發(fā)生后大約15 h后，各個能段高能粒子的單向微分通量漸漸回歸到背景大?。▓D8）。

圖7 紅色觀測點位置處的SEP事件單向積分通量數(shù)值模擬結(jié)果Fig.7 The simulated unidirectional integral flux at location of the red spot

圖8 紅色觀測點位置處的SEP事件單向微分通量數(shù)值模擬結(jié)果Fig.8 The simulated unidirectional differential flux at location of the red spot

圖9～10是位于黃道面上90°，1.5 AU 處高能粒子的數(shù)值模擬結(jié)果。由于事件開始時刻觀測點的磁力線并未很好地連接在激波波前，因此起初在該位置較低能段高能粒子的通量沒有抬升，而較高能段的高能粒子由于橫向的擴散，其通量僅有較小的抬升，在最初的5 h 左右，＞50 MeV 的高能質(zhì)子通量的數(shù)量級從10-2上升為10-1。隨著激波向外傳播，該處的磁力線在激波波前的連接點（COBPoint）向激波中間區(qū)域移動，低能段的高能粒子通量也開始抬升，但是由于激波在傳播過程中逐漸減弱，對高能粒子的加速作用也降低，因而高能段的高能粒子通量的抬升也不大，從圖10中可以看到，18 MeV的高能粒子單向微分通量最高到了大約0.4 個/（cm2s sr MeV），而1.1 MeV的高能粒子的單向微分通量最高則上升了到了接近100個/（cm2s sr MeV）。

圖9 藍(lán)色觀測點位置處的SEP事件單向積分通量數(shù)值模擬結(jié)果Fig.9 The simulated unidirectional integral flux at location of the blue spot

圖10 藍(lán)色觀測點位置處的SEP事件單向微分通量數(shù)值模擬結(jié)果Fig.10 The simulated unidirectional differential flux at location of the blue spot

圖11～12是位于黃道面上 5°，1.5 AU 處的 SEP數(shù)值模擬結(jié)果。在事件初始時刻，此處磁力線在激波面上的連接點（COBPoint）靠近左側(cè)，伴隨連接點逐漸向激波中心區(qū)域移動，高能粒子的通量開始抬升，到達(dá)最高點后開始下降。隨著連接點逐漸向激波右側(cè)邊緣移動，激波的加速作用很快降低，在大約25 h 后觀測點處的高能粒子通量基本恢復(fù)到背景水平。

圖11 黃色觀測點位置處的SEP事件單向積分通量數(shù)值模擬結(jié)果Fig.11 The simulated unidirectional integral flux at location of the yellow spot

圖12 黃色觀測點位置處的SEP事件單向微分通量數(shù)值模擬結(jié)果Fig.12 The simulated unidirectional differential flux at location of the yellow spot

3 結(jié) 論

本文利用iPATH模型對發(fā)生于2014年04月18日爆發(fā)的CME 所驅(qū)動的激波引起的行星際空間太陽高能粒子事件進(jìn)行了數(shù)值模擬。數(shù)值模擬局限于2個天文單位距離以內(nèi)的二維黃道面內(nèi)。本文在模擬中選取了4個具有代表性的觀察點，其中一個觀察點是地球所在位置。通過對比模擬結(jié)果和GOES衛(wèi)星的觀測結(jié)果來確定合適的數(shù)值模擬參數(shù)，并以此為基礎(chǔ)擴展到了更遙遠(yuǎn)的深空。結(jié)果表明：黃道面內(nèi)不同距離不同角度的位置，在一次SEP事件中可以經(jīng)歷完全不同的過程。

現(xiàn)階段太陽高能粒子事件的數(shù)值模擬還較為粗略，由于整個SEP事件的物理過程非常復(fù)雜，iPATH模型做了必要的簡化，所采用的方法和結(jié)果不但揭示了深空環(huán)境中SEP的物理過程，而且也為深空環(huán)境的高能粒子輻射分析提供了參考。深空探測計劃必然需要考慮項目實施過程中來自太陽高能粒子事件的影響，本研究工作所獲得的成果可以在這方面得到應(yīng)用。