動態地球磁層的時空剖分編碼*

王慈楓 胡曉彥 鄒自明 李云龍 白 曦

1(中國科學院國家空間科學中心 北京 100190)

2(中國科學院大學 北京 100049)

3(國家空間科學數據中心 北京 100190)

0 引言

地球磁層是大多數航天器運行的區域，也是太陽等外源因素影響地球空間的關鍵區域，具有重要的研究與應用價值。隨著數字化、智能化、網絡化傳感器技術等現代探測技術的快速發展與大規模應用，對地球磁層的觀測越發全面、立體、精細。當前，空間物理對地磁資料的組織主要包括兩種方式。第一種是基于數據的語義特征，例如時間、載荷等，構建一系列資源間的層級關系，從而實現資料組織和索引[1]。然而，不同衛星、不同臺站在數據預處理時會建立各自獨立的組織和索引方式，因此這種方式不利于數據的檢索、整合、共享和分發。第二種是基于坐標系的重組和索引，空間物理領域結合學科特征提出了一系列專用坐標系。這些坐標系大致可以分為三類，即地心坐標系、日心坐標系和局地坐標系[2]。這種方式進行數據資料組織的基本思想是對現有資料的時空信息在特定應用場景下進行坐標轉換等處理，然后利用文件系統或傳統的關系型數據模型進行重組和存儲。選擇適當的坐標系可以使資料排列更加合理、計算結果更容易理解。然而，這類模型在不同的應用場景下需要選擇不同的坐標系，且這些坐標系的定義有明顯的學科特征，因此不利于多尺度、跨學科領域的數據組織和計算研究。

在研究空天地一體化地球大數據模型時，提出了地球剖分格網的概念，即對地球球體進行不斷細分，得到多層級的剖分格網，并設計編碼方案實現格網編碼，從而便于格網在計算機中的表達和計算，構建統一、多分辨率、可計算的時空框架。有研究將剖分、編碼的思想應用到了日地空間物理領域的相關研究，提出了一系列涵蓋近地空間部分區域的剖分模型，如GeoSOT-3D 立體剖分模型[3]、球體退化八叉樹剖分格網模型（Sphere Degenerated Octree Grid，SDOG）[4]、層級三角網時空模型（Hierarchical Triangular Mesh –Sphere & Time，HTM-ST）模型[5]等。從這些模型的研究對象出發考慮，對其空間范圍的描述可遵循以地球為中心向外延展至地表、近地空間各圈層乃至外部空間的基本思路。因此，針對規則（橢）球面、（橢）球體提出的剖分方案可以對研究對象的空間位型進行有效擬合。但是，地球磁層是受地球磁場、太陽風等多源因素綜合影響形成的一個非規則化的動態物理空間，因此以（橢）球面、（橢）球體為基礎的空間剖分模型及其編碼方案難以在所有應用場景下對磁層空間進行有效表達。地球磁層時空剖分模型基于粒子運動的漂移殼并結合磁場物理要素，實現了一定時空范圍內的磁層區域剖分。該模型得到的格網形變穩定且可以反映磁場要素的物理特征，適合作為地球磁層區域的通用剖分模型。

格網編碼是剖分得到的各時空格網在計算機中的數字化表達。在剖分模型的基礎上實現剖分格網的編碼是構建時空基礎框架至關重要的內容，支撐著數據的快速索引及高效計算。目前，編碼方案大致可以劃分為層次編碼、填充曲線編碼和整數坐標編碼三類[6]。層次編碼是在給定初始層次單元的碼元后，根據剖分格網之間的層次結構，對其子單元用后綴（前綴）碼元表示，適用于父子單元邊界重合、層次關系明確的剖分模型。例如，HTM-ST 模型[5]基于層次編碼方案設計了模型對應的編碼，完成了衛星軌道面的可計算時空框架構建。填充曲線編碼是一種通過遞歸覆蓋指定區域的一維曲線。SDOG 模型基于退化Morton 曲線編碼實現了數字虛擬球體的構建[4]；圈層格網模型基于Hilbert 填充曲線設計了模型對應的編碼方案，實現了地球圈層空間基礎框架的構建[7]。整數坐標編碼是最簡單直接的格網編碼方案，在格網空間中定義m個坐標軸，則格網的m維整數坐標即為其對應的編碼。填充曲線編碼和整數坐標編碼可以很好地反映格網在指定域中的分布關系，使格網鄰域計算更高效；而層次編碼則可以更便捷、快速地計算格網之間的父子關系。這些傳統的編碼方案適用于規則的、時空關系較為明確的各類格網表達，與地理坐標系等各類傳統坐標系之間的轉換較為簡單，可以支持基于這些編碼的格網間時空關系計算。但是，由于磁層空間位型的特殊性，地球磁層時空剖分模型得到的剖分格網往往是不規則的，且格網之間的時空關系較為復雜，用一維編碼難以實現完整表達。

本文對地球磁層時空剖分模型格網的基本特征進行分析。基于剖分格網的時空特性，設計編碼方案對其進行數字化表達，從而為粒子輸運過程、輻射帶動態變化等研究場景提供一種新的時空框架。最后，設計實驗分析編碼效率，驗證編碼方案的高效性。在此基礎上，設計基礎相鄰關系計算效率對比實驗，驗證采用不同編碼方案進行高效時空計算的可行性，為后續復雜計算與科學分析奠定基礎。

1 地球磁層時空剖分模型

1.1 漂移殼的構造

基于漂移殼實現地球磁層動態剖分，首先要根據磁場特征，實現漂移殼的構造。漂移殼參數L為描述漂移殼的重要參數之一。McIlwain[8]將L定義為磁鏡點磁場強度Bm和積分不變量I之間的關系，并提出在偶極磁場條件下，L的幾何含義為磁力線與磁赤道面交點到地心的距離與地球半徑的比值。一般而言，L≤3 的范圍內，地球磁場可以近似為偶極磁場，故而可以將磁層空間分為L≤3 的偶極磁場區域和L ＞3的非偶極磁場區域[9]。

漂移殼由一系列磁力線構成，偶極磁場區域的磁力線方程可以表示為[10]其中，θ和φ分 別為磁緯和磁經，r為徑向距離（其單位為地球半徑Re），C為常數。此時，可以直接利用上述方程構造漂移殼。偶極磁場區域L={1.5,2,2.5}的各漂移殼如圖1 所示。

圖1 偶極磁場區域L ={1.5,2,2.5}的各漂移殼構造Fig. 1 Construction of the drift shells with L={1.5,2,2.5}in the dipole field

非偶極磁場區域以Galperin 提出的L值計算方法為基本思想，漂移殼構造方法如下。

（1）利用式（1）計算L=L0的面，稱為參考L面（ReferenceLSurface），該曲面與磁赤道面（Magnetic Equatorial Plane）的交線用ΘF表示為

（2）從ΘF上任意一起點P出發，在IGRF 磁場模型[11]下利用龍格庫塔–梅森法[12]，向地磁北極方向追蹤磁力線至地球表面，交點記為Foot Point North。

（3）從Foot Point North 出發，在IGRF+Tsyganenko96（T96）磁場模型[11,13]下向地磁南極方向追蹤磁力線至地球表面，若該磁力線不能在預先設置的有限步驟內閉合，則認為該磁力線為開放磁力線，對應暴露于行星際太陽風的區域，不在本文模型中進行剖分。

（4）針對ΘF上所有選取的起點執行步驟（2）（3），形成的所有磁力線構成的曲面即為非偶極磁場區域漂移殼。

上述方法對其中一條磁力線的構造過程如圖2所示。使用外源場模型T96 時所需的太陽風參數、地磁指數等從OMNI 數據庫獲取.

圖2 基于Galperin 計算方法的磁力線構造方案Fig. 2 Construction of a field line of the drift shell based on the calculation of Galperin L

非偶極磁場區域漂移殼的構造結果與時間密切相關，本文統一用yyyyddd表 示時間，其中ddd為年積日。t=2015169非 偶極磁場區域L={6.5,7,7.5}的各漂移殼構造結果如圖3 所示。

圖 3 非偶極磁場區域 t =2015169， L={6.5,7,7.5}的各漂移殼構造Fig. 3 Construction of the drift shells witht=2015169 and L ={6.5,7,7.5} in the non-dipole field

1.2 剖分模型

以上述漂移殼構造方法為基礎，提出一種地球磁層時空剖分模型。剖分后得到一系列具有時空內在聯系的時空格網。為方便論述，對地球磁層時空格網相關術語定義如下。

時間間隔 Δt(N)： 剖分N次時兩個時間格網之間的時間差。

徑距ΔL： 兩個漂移殼的參數值L之差。

緯距 Δθ(N)： 同一時間下，空間剖分N次時兩個空間格網之間最大緯線與最小緯線的緯度差。

磁力線距ΔφF(N)： 同一時間下，空間剖分N次時空間格網左下頂點與右下頂點的經度差。

經距 Δφ(Nφ)： 在漂移殼上進行Nφ次等經度剖分，兩條經線之間的經度差。

在空間剖分方面，該模型首先在徑向上選擇一系列不同的L值，將磁層空間離散為一組漂移殼。針對不同場景，可選擇等徑距離散或根據粒子的分布特征在徑向上進行適應性離散，隨后對每個漂移殼曲面進行剖分。漂移殼由磁力線段構成，沿磁力線對漂移殼進行剖分可以使剖分格網與磁場條件相關聯，以支持基于磁場的時空關系計算。沿磁緯方向采用等緯度間隔剖分，使計算量最小。

在非偶極磁場條件下，漂移殼的構造與時間密切

由剖分模型可知，偶極磁場區域的磁力線與磁經線保持一致，即在同一條磁力線上，所有點的磁經度相同。因此，有N=Nφ且

（1） 磁力線距不確定性。格網間的磁力線距ΔφF(N)是 磁場決定的，不同格網之間的ΔφF(N)不同。

（2） 層次關系不確定性。在進行下一層級的剖分時，是對Θeq進行遞歸二等分，得到起點后構造磁力線，因此剖分層級的格網之間的層次關系難以確定。

（3） 時空關系復雜性。在不同時間格網下，漂移殼的空間位型不同。因此不同時間下，同一L值對應的漂移殼之間的時空關系極為復雜。本文暫不討論時空關系計算，只在同一時間的條件下討論基礎空間關系計算。

2 時空剖分模型編碼方案

為了構建支持高效檢索與計算的時空框架，需要對剖分得到的各時空格網進行編碼表達，以便于計算機對其存儲和處理。編碼過程中用到的符號標記定義如下：N表示該編碼對應的格網剖分層級，上標表示在某一維度的編碼，下標表示編碼進制，如(idθ(N))2表 示剖分N次時，格網緯度對應的二進制編碼。在下文出現的二進制計算中，乘法和加法都是按位計算。

圖4 L =1.5的漂移殼剖分Fig. 4 Subdivision of the drift shell with L=1.5

圖5 非偶極磁場區域L =7.5， t =2015169漂移殼剖分Fig. 5 Subdivision of the drift shell with L =7.5,t=2015169

利用一個5 位二進制編碼(idL)2表達格網所在的漂移殼，從L=1.1向外進行順序編碼，即

每個時空格網的編碼由其所在的漂移殼順序編碼和該格網在漂移殼上的時空位置編碼串聯構成，即

其中，⊕表示編碼串聯。

2.1 Morton 編碼

基于Morton 曲線對剖分模型得到的時空格網進行編碼，基本思想為分別在磁緯線（θ）、磁力線（F）、和時間段（t）三個維度上進行順序編碼，記為(idθ(N))10，(idF(N))10，(idt(N))10，再耦合得到對應的時空Morton 編碼。

在偶極磁場區域，具體編碼步驟如下。

步驟1 分別計算三個維度的十進制順序編碼，有

其中， (idST(N))2的 編碼位數為 ( 3N)，其他二進制編碼的位數為N。

步驟3 與漂移殼順序編碼串聯得到完整編碼(id(N))2。例如，T0=2014005，對于一個坐標為(t,L,θ,φ)=(2015169,1.5,0,0)的時空格網，令剖分次 數N=2 時， ΔT=273 d，則 有 (idL)2=00000，(idt(2))10=2， (idθ(2))10=2， (idF(2))10=0，耦 合得(idST(2))2=110000，最終計算得到該格網的(id(2))2=00000110000。

非偶極磁場區域的編碼方式與偶極磁場區域類似。但是，對于同一條磁力線上的點可能處于不同磁經度，因此(idF(N))10的計算方法需進行修正。通過計算每條磁力線與磁赤道面（θ=0）的交點經度，再按照交點經度大小對磁力線進行排序得到順序編碼(idF(N))10。 然后將 (idF(N))10代入步驟2 和3，得到非偶極磁場區域Morton 編碼(id(2))2。

2.2 整數坐標編碼

整數坐標編碼是最簡單直接的編碼方案，其基本思想是在空間中定義m個坐標軸，將剖分格網按照其沿著坐標軸前進的步長進行編碼。偶極磁場區域的整數坐標編碼方案步驟如下。

步驟1 根據輸入的時空坐標 (t,L0,θ,φ)分別計算其在磁緯（θ）、磁經（φ）、時間（t）三個維度上的整數坐標，有

2.3 漂移殼剖分格網編碼方案

偶極磁場區域磁力線均勻對稱， 有(idF(N))10=(idφ(N))10。因此，任意一種編碼(id(N))2都可以較為完整地反映格網的時空信息及其之間的時空關系。但是在非偶極磁場區域，Morton 編碼缺乏(idφ(N))10的 信息，即在θ／=0的緯線上，(id(N))2無法直接反映格網的空間左右相鄰關系。而整數坐編碼則缺少(idF(N))10的信息，因此無法反映兩個空間格網是否在一條磁力線上，較難計算空間上下相鄰關系。針對非偶極磁場區域，本文結合Morton 編碼和整數坐標編碼，設計三種能完整表達格網時空信息的漂移殼剖分格網時空編碼（Drift Shell Grids Coding, DSGC），方案如下。

方案1

其 中： (idF(N))2的 編 碼 位 數 為N； (idt(N))2，(idθ(N))2及 (idφ(N))2的 編 碼 位 數 為 max(N,Nφ)。則(idS1T(N))2的 編碼位數為3 ×max(N,Nφ)。最后串聯漂移殼編碼，即

3 實驗及結果分析

3.1 編碼效率對比

實驗使用的軟硬件配置列于表1。這里使用了中國科學院國家空間科學中心公共技術服務中心空間科學數據融合計算平臺的計算服務。

表1 集群環境Table 1 Cluster environment

圖6 非偶極磁場區域漂移殼格網編碼Fig. 6 Diagram of the DSGC for the grids in the non-dipole field

在偶極磁場和非偶極磁場兩個特征區域內，選擇典型的漂移殼對其進行N次剖分，然后對剖分得到的所有時空格網進行編碼，計算并分析各編碼方案的編碼效率（ηc），則有

式中，nc表示總時空格網數，Tc表示編碼時間。

偶極磁場區域選擇L=1.5，各剖分層級上的編碼效率列于表2。非偶極磁場區域選擇L={7.5,8}，時間t={2014005,2015169}這4 個漂移殼，計算各剖分層級上的平均編碼效率。

從表2 可以看出，隨著剖分次數的增多，剖分得到的格網數成倍增加，但編碼效率較為穩定，且Morton 曲線編碼和整數坐標編碼方案的效率相差不大。由表3 可知，漂移殼剖分格網編碼方案與傳統的Morton 編碼方案和整數坐標編碼方案的效率相差不大，且當剖分到12 層級時仍能維持在4500 s–1以上。可以認為，本文提出的編碼方案編碼效率較高。并且對于非偶極磁場區域的格網而言，Morton 編碼方案與整數坐標編碼方案難以完整表達格網的時空信息，給后續時空關系的計算帶來一定的困難。

表2 偶極磁場區域編碼效率Table 2 Encoding efficiency in the dipole field

表3 非偶極磁場區域平均編碼效率Table 3 Mean encoding efficiency in the non-dipole field

3.2 相鄰關系計算效率對比

基于非偶極磁場區域的三類編碼，設計基礎的空間相鄰關系計算，以分析對應編碼在進行基礎計算時的效率。在同一時間下，一個空間格網的相鄰關系包括上鄰居（T）、下鄰居（B）、左鄰居（L）、右鄰居（R）、左上鄰（TL）、左下鄰（BL）、右上鄰（TR）、右下鄰（BR）。其中，上下左右4 個基礎相鄰關系的計算列于表4。首先分別計算左鄰居的上鄰居以及上鄰居的左鄰居，作為左上鄰兩個候選者，然后選擇距離所求格網較近的候選者作為其左上鄰，左下鄰、右上鄰、右下鄰計算方式與此類似。

表4 剖分格網的基礎相鄰關系Table 4 Basic adjacency of the grids

選擇L={1.5,7.5}， 時間t=2015169漂移殼，分別對其進行N=10次剖分，然后在其上隨機選取10000 個時空點，分別對其編碼并基于三類編碼求其對應的8 個鄰居，得到一次計算的效率

式中，ns表 示所取樣本數，ts表示計算所有樣本的8 個鄰居的總時間。進行10 次計算后求平均值η-s，得到基礎相鄰關系計算效率列于表5 和表6。

由表5 可知，偶極磁場區域兩種編碼方案在計算相鄰關系時的效率較高且兩者相差不大。因此，對偶極磁場區域可以根據物理存儲結構選擇合適的編碼方案。由表6 可知，一維編碼比二維編碼在進行基礎相鄰關系計算時效率低，且計算左右相鄰關系時差距尤為明顯。可以說明，非偶極磁場區域進行簡單的時空計算時，二維編碼有一定的優勢。

表5 偶極磁場區域格網基礎相鄰關系計算效率對比Table 5 Efficiency comparison of the basic adjacency of the grids in the dipole field

表6 非偶極磁場區域格網基礎相鄰關系計算效率對比Table 6 Efficiency comparison of the basic adjacency of the grids in the non-dipole field

4 討論與結論

分析了地球磁層時空剖分模型得到的非偶極磁場區域剖分格網的基本特征，將其概括為磁力線距不確定性、層次關系不確定性以及時空關系復雜性三點。根據這些特征，本文融合整數坐標編碼及Morton 曲線編碼的基本思想，提出了三種漂移殼剖分格網編碼方案，并設計實驗分析了各方案編碼效率及基于對應方案的相鄰關系計算效率。實驗證明，本文提出的編碼方案效率較高且可以完整地表達格網之間的時空關系，因此可以較好地支持相鄰關系計算。此外，探討了結合物理場來構建時空框架的方案，試圖為后續的科學計算與分析決策提供便利。

但模型仍存在若干局限。第一，本文設計的實驗是基于邏輯模型進行的，沒有考慮各編碼在計算機中的物理層存儲設計。幾種編碼方案雖然編碼效率接近，但是在不同的存儲模型下，存取速度存在差異，需要根據編碼特征設計存儲方案。尤其是二維編碼(id3(N))2存儲到計算機中要保證兩個維度的編碼都可以唯一地、確定地標記時空格網，因此其對應的存儲模型需針對這種特性進行優化設計。第二，在空間適用性方面，完成時空框架構建后要進行數據映射才能夠實現該框架在數據組織及高效計算上的應用。而進行數據映射時，對于L值較大的漂移殼（例如，L≥10），可能會引入較大的位置配準誤差。因此，基于特定的應用場景提升模型適用性也是未來研究方向之一。

致謝 計算服務支持由國家科技資源共享服務平臺–國家空間科學數據中心（https://www.nssdc.ac.cn）提供。