楚偉, 秦剛, 黃建平, 許嵩, 澤仁志瑪, 申旭輝

1 應(yīng)急管理部國家自然災(zāi)害防治研究院， 北京 100085 2 哈爾濱工業(yè)大學(xué)(深圳)， 深圳 518055

0 引言

太陽高能粒子是由太陽耀斑和日冕物質(zhì)拋射爆發(fā)時(shí)導(dǎo)致行星際和磁層空間粒子環(huán)境增強(qiáng)的事件(Le and Zhang, 2017; Le et al., 2017; Zhao et al., 2018).如果太陽質(zhì)子事件中的質(zhì)子能量達(dá)到相對(duì)論的能量，這些相對(duì)論能量的質(zhì)子可以穿越磁層，并與地球大氣發(fā)生核反應(yīng)，從而使得地面儀器探測到相對(duì)論質(zhì)子與地球大氣發(fā)生核反應(yīng)的次級(jí)成分，這種現(xiàn)象稱為地面水平增強(qiáng)事件，簡稱Ground Level Enhancements(GLE)事件.這種事件主要出現(xiàn)在太陽活動(dòng)高年(Le and Liu, 2020; Zhao and Le,2020).太陽高能粒子事件是偶發(fā)的，而銀河宇宙線是恒定存在的，只是其強(qiáng)度受太陽活動(dòng)的調(diào)制.銀河宇宙線和太陽高能粒子進(jìn)入磁層時(shí)會(huì)導(dǎo)致磁層高能粒子環(huán)境更惡劣，將對(duì)磁層空間的高能粒子環(huán)境產(chǎn)生重要的影響.高能質(zhì)子能夠穿透電子元件防護(hù)層導(dǎo)致單粒子反轉(zhuǎn)，對(duì)宇航員的健康構(gòu)成威脅，同時(shí)也會(huì)降低太陽能電池的壽命，影響極區(qū)短波通訊等，嚴(yán)重情況下將造成電子元器件損傷引發(fā)衛(wèi)星故障，甚至將衛(wèi)星徹底摧毀(Baker, 1998, 2001, 2002).因此研究行星際和磁層高能粒子耦合機(jī)制是日地空間物理的重點(diǎn)問題.

從上述公式中可以發(fā)現(xiàn)，當(dāng)磁緯度為90°時(shí)，截止剛度達(dá)到最小.但上述公式中沒有給出截止剛度與粒子投擲角的相關(guān)關(guān)系，同時(shí)后續(xù)研究人員更多地將研究內(nèi)容集中在垂直方向上，即天頂方向.Smart和Shea(2005)指出，在計(jì)算真實(shí)的截止剛度時(shí)，由于行星際高能粒子的各項(xiàng)同性，需要考慮來自不同方向粒子的截止剛度(Masarik and Reedy, 1995).同時(shí)行星際高能粒子的投擲角作為一個(gè)重要參量(可以研究粒子的加速以及絕熱聚焦效應(yīng))，對(duì)研究行星際和磁層高能粒子耦合具有至關(guān)重要的作用，因此針對(duì)行星際高能粒子的傳輸研究都需要考慮投擲角的擴(kuò)散狀況.

Shea等(1965)使用數(shù)值模擬方法研究了地球表面超過300個(gè)點(diǎn)的垂直截止剛度，發(fā)現(xiàn)在南非、南大西洋和加那利群島附近，其截止剛度和地面觀測值具有高達(dá)15%的偏差.Daniel和Stephens(1966)研究了印度的HYDERABAD臺(tái)站附件的高能粒子截止剛度的方向依賴性.Bland和Cioni(1968)研究了有限方向(方位角方向15°間隔和極角方向10°間隔)的非垂直方向截止剛度，并利用兩個(gè)不同觀測位置(Trapani和Aire-sur-Adour)氣球測量試驗(yàn)發(fā)現(xiàn)，測量和計(jì)算得到的截止剛度在上述兩個(gè)不同的觀測位置的均方差分別為4.8%和4.6%.Fanselow和Stone(1972)研究表明衛(wèi)星實(shí)際觀測結(jié)果比理論計(jì)算結(jié)果得到的垂直方向截止緯度范圍通常要低3°～5°，也就意味著垂直方向的截止剛度應(yīng)該不是該點(diǎn)的最小地磁截止剛度.Cooke等(1991)針對(duì)變化磁場引起的截止剛度變化進(jìn)行了研究，并針對(duì)計(jì)算過程中只考慮粒子剛度而忽略粒子方向的做法進(jìn)行了指正.Bhattacharyya和Mitra(1997)研究了截止剛度隨著變化磁場的變化關(guān)系.Smart和Shea兩人針對(duì)截止剛度問題做了一系列的工作(Shea et al. 1965，1968； Smart，1999; Smart et al.,1999, 2000; Shea and Smart， 2001; Smart and Shea，2003a，b，2005)，針對(duì)地面、衛(wèi)星觀測和數(shù)值模擬存在的偏差展開了論述，并得到在進(jìn)行衛(wèi)星觀測和數(shù)值模擬結(jié)果精確匹配過程中需要考慮不同方向粒子入射狀況的結(jié)論.Derome和Buénerd(2001)發(fā)現(xiàn)衛(wèi)星觀測到的剛度小于模擬獲得的截止剛度.Shimazu等(2006)研究了行星際磁場方向?qū)刂咕暥鹊挠绊?Dorman等(2008)研究了非垂直方向的截止剛度，發(fā)現(xiàn)當(dāng)剛度為7～8 GV時(shí)，垂直方向和非垂直方向的截止剛度相差可能達(dá)到1 GV；Dmitriev等(2010)通過POES衛(wèi)星的觀測數(shù)據(jù)，確定了高能粒子截止緯度的橢圓形表達(dá)式.國內(nèi)研究者通過數(shù)值模擬(甄杰和楚偉，2013； Chu and Qin, 2016)和觀測研究(朱邦耀，1979；樂貴明等，2003)發(fā)現(xiàn)地磁活動(dòng)對(duì)截止剛度的影響.

由于使用解析表達(dá)式的垂直方向高能粒子截止剛度不能真實(shí)反映實(shí)際狀態(tài)的磁層高能粒子截止剛度，為了盡可能準(zhǔn)確反映真實(shí)磁層的地磁截止剛度，需要研究真實(shí)磁場狀態(tài)下不同方向(投擲角)的截止剛度.

通過測試粒子的方式對(duì)不同能量(剛度)的粒子進(jìn)行倒向追蹤計(jì)算，得到給定點(diǎn)的截止剛度.由于地球磁場的結(jié)構(gòu)，使用此種方法需要計(jì)算包括回旋運(yùn)動(dòng)在內(nèi)的高能粒子運(yùn)動(dòng)，同時(shí)需要追蹤大量的粒子運(yùn)動(dòng)，因此往往需要耗費(fèi)大量的計(jì)算資源和計(jì)算時(shí)間，因此本文只給出地磁平靜期間的典型高能粒子截止剛度狀態(tài).

1 數(shù)值模擬方法

1.1 高能粒子剛度

高能粒子剛度作為描述高能粒子能量大小的物理量，其定義如下：R=mvc/Ze，其中R表示粒子剛度，m為高能粒子質(zhì)量，v為粒子速度，c為光速，Z為電荷數(shù)，e為元電荷電量.剛度在數(shù)值上等于單位單核的能量，單位為V.截止剛度表征到達(dá)給定點(diǎn)的最小剛度.地磁粒子截止剛度可以描述地磁場對(duì)高能粒子的屏蔽作用.Le等(2006)推導(dǎo)了粒子剛度和能量的關(guān)系式，依據(jù)該關(guān)系式，我們可計(jì)算不同剛度粒子的能量.圖1為剛度與質(zhì)子能量的關(guān)系圖，從圖上可以發(fā)現(xiàn)，兩者不是簡單的線性關(guān)系.

圖1 質(zhì)子剛度Rigidity與能量Energy的關(guān)系圖Fig.1 The relationship between rigidity and energy for proton

Weygand和Raeder(2005)研究表明，使用目前通用的背景場模型都可以產(chǎn)生比較可靠的截止剛度模式，其中使用T96模型(Tsyganenko，1995； Tsyganenko and Stern，1996)作為外源場計(jì)算得到的截止剛度對(duì)應(yīng)的緯度與衛(wèi)星觀測結(jié)果平均偏差為4.0°±1.4°，而使用T89模型(Tsyganenko，1989)作為外源場計(jì)算得到的數(shù)值為3.9°±2.4°，兩者的偏差并不明顯.由于本文研究平靜期間的截止剛度狀況，不采用T05(Tsyganenko and Sitnov，2005)模式，同時(shí)為了計(jì)算快捷方便，采用T89模式(Tsyganenko，1989)作為外源場.使用單粒子運(yùn)動(dòng)理論，結(jié)合地球磁層的內(nèi)源背景場IGRF12模型以及T89外源場模型作為背景磁場，倒向追蹤粒子運(yùn)動(dòng)狀態(tài)，得到不同投擲角的粒子截止剛度，從中選擇最小剛度作為給定點(diǎn)的截止剛度.在計(jì)算過程中，只考慮背景磁場作用，忽略電場影響，因此運(yùn)動(dòng)過程中粒子的能量/剛度不會(huì)發(fā)生改變.同時(shí)由于是計(jì)算行星際粒子進(jìn)入磁層，不考慮回旋運(yùn)動(dòng)和彈跳運(yùn)動(dòng)的投擲角損失錐造成的影響.

1.2 模擬條件的設(shè)定

模擬過程中采用四階龍格庫塔積分，初始高度設(shè)定為450 km高度.內(nèi)邊界設(shè)定為1個(gè)地球半徑，外邊界由TS模式模型給出的磁層頂為邊界.同時(shí)為了保證計(jì)算效率最大化，排除捕獲粒子長時(shí)間做彈跳運(yùn)動(dòng)的影響，模擬過程中最大運(yùn)動(dòng)時(shí)長設(shè)定為1小時(shí)，最大運(yùn)動(dòng)距離設(shè)定為1000個(gè)地球半徑，最大運(yùn)動(dòng)步數(shù)為106，這樣既能確保計(jì)算結(jié)果的準(zhǔn)確性，又能有效提高計(jì)算效率.

積分步長的設(shè)定將在一定程度上決定計(jì)算速度以及計(jì)算準(zhǔn)確性.對(duì)于積分步長的設(shè)定目前主要存在兩種方法.第一種方法是使用龍格庫塔積分程序rkqs(Levin，1998)，設(shè)定誤差范圍內(nèi)自適應(yīng)的調(diào)節(jié)積分時(shí)間步長；第二種方法為使用粒子回旋周期的一定比例(通常情況下設(shè)定為粒子局地回旋周期的1/100，詳見Smart et al., 2000；Kress et al., 2010).在本文模擬過程中，在準(zhǔn)確描述粒子運(yùn)動(dòng)前提下，盡量提高粒子積分效率.我們采用第二種方法，模擬過程中時(shí)間步長設(shè)定為粒子回旋周期的1/100.由于緯度對(duì)磁場的影響遠(yuǎn)遠(yuǎn)大于經(jīng)度，因此模擬過程中緯度采用1°的間隔，經(jīng)度采用10°的間隔.同時(shí)鑒于目前近地空間探測器的投擲角分辨率普遍大于5°，因此投擲角采用5°的間隔，粒子初始垂直平面方向運(yùn)動(dòng)速度進(jìn)行均勻分配(Weygand and Raeder，2005).

我們?cè)诖耸褂玫膬?nèi)部場模型是國際地磁參考場(IGRF)模型，該模型是地球主要磁場及其長期變化的標(biāo)準(zhǔn)數(shù)學(xué)描述.IGRF模型針對(duì)特定年份進(jìn)行了標(biāo)準(zhǔn)化處理，反映了當(dāng)時(shí)可用的最精確值.為了計(jì)算內(nèi)源場和外源場貢獻(xiàn)，并進(jìn)行計(jì)算過程中的坐標(biāo)轉(zhuǎn)換以及磁層頂外邊界確認(rèn)，我們使用Tsyganenko提供的軟件包GEOPACK-2008.以上所有有關(guān)磁場和輸入的參數(shù)均可從http:∥geo.phys.spbu.ru/～tsyganenko/modeling.html網(wǎng)站下載.

為了進(jìn)一步提高計(jì)算效率，背景磁場使用插值進(jìn)行計(jì)算，首先計(jì)算每一個(gè)格點(diǎn)的背景磁場，然后通過線性插值方法得到粒子運(yùn)動(dòng)位置的磁場強(qiáng)度.

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

圖2 計(jì)算過程中的截止剛度的半影區(qū)示意圖允許(白色)和禁止(黑色)在西經(jīng)40°、北緯30°垂直方向，距離地面450 km高度處，剛度計(jì)算間隔為1 MeV的最大截止剛度、最小截止剛度和有效截止剛度示意圖.Fig.2 Schematic diagram of penumbra area of cut-off rigidityThe maximum, minimum and effective cut-off rigidity at intervals of 1 MeV for the calculation of the rigidity at a height of 450 km from the ground in the vertical direction of 40 degrees west longitude, 30 degrees north latitude for allowed (white) and prohibited (black).

2.1 不同投擲角截止剛度的對(duì)比

為了使研究結(jié)果具有普遍性，此研究中我們選取地磁平靜期間進(jìn)行模擬.

圖3為截止剛度的全球分布圖像，左上圖為垂直方向的截止剛度分布圖像，右上圖為非垂直方向?qū)?yīng)的截止剛度圖像，左下圖為全向的最小截止剛度圖像(包括不同方向的粒子截止剛度)，右下圖為垂直方向截止剛度和最小截止剛度的差值全球分布圖像.對(duì)比三個(gè)全球截止剛度的分布圖像，可以發(fā)現(xiàn)，通常所說的垂直方向(也就是天頂方向)一般意義上不是最小有效截止剛度對(duì)應(yīng)的方向，尤其是在中緯度地區(qū)以及西經(jīng)0～90°范圍內(nèi)的差值最大，差值比例可能高達(dá)60%至70%.因此如果使用垂直方向的截止剛度作為地磁有效截止剛度，將會(huì)出現(xiàn)普遍高估的情形.從全球平均意義來看，垂直方向截止剛度將平均高估13%左右，最大偏差將達(dá)到70%，由此可見研究不同投擲角的高能粒子進(jìn)入磁層的狀態(tài)對(duì)研究地磁截止剛度具有重要意義.從圖4可以發(fā)現(xiàn)，全方向計(jì)算得到的最小截止剛度和垂直方向計(jì)算得到的最小截止剛度偏差隨著緯度呈現(xiàn)蝴蝶狀分布形態(tài)，這也意味著中緯度地區(qū)的偏差比值達(dá)到最大.

2.2 最小截止剛度對(duì)應(yīng)的投擲角的全球統(tǒng)計(jì)

圖5為最小有效截止剛度對(duì)應(yīng)的投擲角的全球分布圖像(為了研究方便，我們使用185°表示垂直方向角度).從圖中可以發(fā)現(xiàn)，行星際高能粒子到達(dá)磁層的投擲角全球分布沒有明顯經(jīng)緯度依賴性.

高能粒子的投擲角將影響高能粒子由行星際進(jìn)入磁層，導(dǎo)致行星際高能粒子在磁層空間的分布投擲角依賴性明顯.從圖6我們可以發(fā)現(xiàn)，高能粒子由行星際空間進(jìn)入磁層的最強(qiáng)有力的方向?yàn)檠刂帕€方向，也就是投擲角為0°或者180°時(shí)，此時(shí)高能粒子平行方向的能量巨大，可以直接沿著磁力線到達(dá)地球.同時(shí)由于研究的高能粒子能量非常高，天頂方向也是高能粒子進(jìn)入磁層的一個(gè)重要通道.其他的投擲角方向分布符合中心點(diǎn)為90°的正態(tài)分布.

圖3 高能粒子截止剛度全球分布圖像(a) 垂直方向的截止剛度分布圖像； (b) 非垂直方向?qū)?yīng)的截止剛度圖像； (c) 全向的最小截止剛度圖像(包括不同方向的粒子截止剛度)； (d) 垂直方向截止剛度和最小截止剛度的差值全球分布圖像.Fig.3 Global distribution image of cut-off rigidity(a) is vertical cut-off rigidity distribution image; (b) is non-vertical cut-off rigidity; (c) is omnidirectional minimum cut-off rigidity (including particle cut-off rigidity in different directions), and (d) is vertical direction Global distribution of the difference between the cut-off rigidity and the minimum cut-off rigidity.

圖4 截止剛度差值隨著經(jīng)緯度變化圖像Fig.4 Cutoff rigidity difference changes with latitude and longitude

圖5 最小有效截止剛度對(duì)應(yīng)的投擲角全球分布圖像(為了研究方便，我們使用185°表示垂直方向角度)Fig.5 The global distribution of the pitch angle corresponding to the minimum effective cut-off rigidity. (For convenience in the study, we use 185 degrees to represent the vertical direction)

圖6 最小截止剛度對(duì)應(yīng)的投擲角占比柱狀圖(為了方便問題研究，我們采用185°表示垂直方向的角度，因此圖中出現(xiàn)185°的情形)Fig.6 The histogram of the pitch angle corresponding to the minimum cut-off rigidity (To facilitate the study of the problem, we use 185 degrees to represent the angle in the vertical direction, so the situation of 185 degrees appears in the figure)

一般假設(shè)銀河宇宙線為各項(xiàng)同性分布，因此由銀河宇宙線造成的進(jìn)入磁層的高能粒子，其投擲角分布形態(tài)應(yīng)呈現(xiàn)正態(tài)分布，此研究對(duì)磁層中高能粒子的投擲角分布具有重要作用.磁層中高能粒子的投擲角分布主要有三種可能的形態(tài)分布，各向同性分布、蝴蝶型分布以及薄餅型分布.此處獲得截止剛度對(duì)應(yīng)的投擲角分布圖像與薄餅形態(tài)最接近，但是在兩端呈現(xiàn)上翹的形態(tài).

3 小結(jié)

高能粒子，尤其是行星際高能粒子進(jìn)入磁層空間，將對(duì)磁層空間的粒子輻射環(huán)境產(chǎn)生影響.研究高能質(zhì)子的截止剛度，尤其是中高緯度地區(qū)的高能粒子截止剛度，對(duì)研究行星際高能粒子在磁層空間形成的分布具有重要意義.本文使用單粒子數(shù)值模擬方法研究了近地空間高能粒子截止剛度與粒子投擲角的相關(guān)關(guān)系，研究發(fā)現(xiàn)：

(1)天頂方向或者垂直方向的截止剛度通常不是最小地磁截止剛度.

(2)最小地磁截止剛度對(duì)應(yīng)的投擲角方向最大為沿著磁場方向，即0°方向；其次為天頂方向，也就是通常所說的垂直方向；然后為180°方向，即與磁場反方向.這種高能粒子進(jìn)入磁層的投擲角分布狀態(tài)符合理論預(yù)期，因?yàn)檠刂帕€方向是粒子最容易進(jìn)入磁層的方式，因此對(duì)應(yīng)的截止剛度比較小.

(3)全球范圍截止剛度對(duì)應(yīng)的投擲角分布符合兩端上翹的正態(tài)分布形態(tài)，不考慮兩端最大占比，中心位置出現(xiàn)在90°附近.

(4)通過地磁平靜期間的數(shù)值模擬發(fā)現(xiàn)，使用垂直方向的截止剛度對(duì)比最小截止剛度平均將高估13.17%，最大可能高估70%.

(5)不同經(jīng)緯度高能粒子的截止剛度與投擲角不存在明顯關(guān)系.

本文使用了數(shù)值模擬方法進(jìn)行高能粒子的運(yùn)動(dòng)模擬，由于使用的背景磁場和磁層頂邊界模型的不盡精確，得到的結(jié)果可能不能完全反映到達(dá)近地空間最小的高能粒子剛度，但能反映截止剛度的總體趨勢(shì).為了得到更普適的截止剛度模式，在將來工作中，我們將利用數(shù)值模擬結(jié)合衛(wèi)星觀測數(shù)據(jù)，進(jìn)一步修正不同經(jīng)緯度的高能粒子截止剛度，以期做到盡可能精確可靠，為研究近地空間高能粒子環(huán)境變化提供精確的背景模型(Wang et al. 2020).隨著我國第一個(gè)電磁監(jiān)測試驗(yàn)衛(wèi)星的成功發(fā)射，星上搭載的高能粒子探測器將對(duì)近地空間高能粒子展開軌道高度全球范圍、高投擲角分辨率的高能粒子探測數(shù)據(jù)，將有利于對(duì)此問題進(jìn)行更加精確的研究.Rodger(2006)使用地基觀測數(shù)據(jù)和理論預(yù)期進(jìn)行了對(duì)比，發(fā)現(xiàn)在磁擾動(dòng)期間(Kp指數(shù)小于5的情況下)觀測和理論預(yù)期符合得很好，但是當(dāng)Kp指數(shù)大于7，理論預(yù)期要比實(shí)際觀測結(jié)果大很多.我們后續(xù)將針對(duì)擾動(dòng)條件下不同方向的截止剛度展開研究，屆時(shí)將對(duì)比平靜期間與擾動(dòng)期間的變化磁場導(dǎo)致的截止剛度變化狀況.

致謝本研究得到了國家重點(diǎn)研發(fā)專項(xiàng)-地球物理探測衛(wèi)星數(shù)據(jù)分析處理技術(shù)與地震預(yù)測應(yīng)用研究(2018YFC1503502-05)課題、地震活動(dòng)期間ZH-1電離層掩星數(shù)據(jù)異常研究(ZDJ2019-03)課題以及ISSI-BJ(2019IT-33)課題的支持，同時(shí)得到了組內(nèi)各位同事的鼎力幫助，在此感謝同事們的辛勤付出，使得本文得以完成.感謝IGRF提供的內(nèi)源場源代碼以及Tsyganenko提供的外源場模型源代碼.同時(shí)在此對(duì)審稿老師提出的寶貴建議表示感謝.