基于數據網格化方法的低軌輻射帶建模技術

常崢，王詠梅 ，田天，張賢國

(1.中國科學院 國家空間科學中心，北京100190; 2.中國科學院大學 地球科學學院，北京100049;3.61741 部隊，北京100094)

地球外層空間存在著一個區域，其中充滿地磁場捕獲的高能質子和電子，這個區域被稱為地球輻射帶(以下簡稱為輻射帶)［1］.輻射帶中的質子和電子能量較高［2］，能夠引起航天器材料和器件性能退化甚至失效［3-4］.在地球空間運行的絕大多數航天器都要或多或少地穿越輻射帶，遭遇高能粒子輻射.因此，在航天任務的設計過程中，準確地定量估算航天器所處的空間輻射環境對于航天器設計的優化、運行安全的保障和任務的完成都至關重要.

目前國際和國內航天界普遍采用NASA 編制的質子模型AP8 和電子模型AE8 來計算輻射帶質子和電子的通量，并以此為基礎評估輻射帶對航天器的影響［3，5-6］.目前我國可獲得的AP8/AE8模型的最新版本分別完成于1976 年和1983 年，這兩個版本使用的都是20 世紀70 年代以前的數據［7-8］，觀測資料比較陳舊，模型計算值與實際觀測值經常存在較大偏差［4，6，9-12］，使得這兩個版本的模型應用到當前航天任務中存在時效局限性.NASA 在后期對AP8/AE8 模型進行了多次修訂和補充，形成了多種修正版本，并且于2006 年聯合美國多家單位開始開發用于替代AP8/AE8 的輻射帶質子和電子的新模型:AP9/AE9［6］，但由于技術封鎖，我國目前還無法獲得AP8/AE8 的修正版本和AP9/AE9.因此，我國自主研發輻射帶模型對于我國航天工程的可持續發展具有現實和戰略的意義.地面高度為300 ～1 000 km 的軌道空間即低地球軌道(Low-Earth Orbit，LEO)空間是大量應用衛星運行的區域［13］.因此，這一區域是輻射帶建模中最重要的區域之一.

1 我國輻射帶監測數據情況

我國從20 世紀80 年代開始對空間環境展開了大規模探測，其中LEO 空間輻射帶高能粒子分布是探測的重點［14］.截至到目前為止，已發射衛星中搭載有高能粒子探測器的衛星共有50 余顆，其中LEO 衛星超過30 顆.部分高能粒子探測儀器的衛星軌道、觀測對象、能量范圍、通道數和觀測時段如表1 所示.這些儀器的設計原理相同，都為半導體型探測器.

不同儀器獲得的探測數據在聯合使用前，必須經過交叉定標［15］.因此在建模過程中，研究人員對表1 所列數據進行了同化處理，并以此為基礎建立了LEO、中地球軌道(Medium-Earth Orbit，MEO)和同步軌道(Geosynchronous Orbit，GEO)輻射環境數據庫，數據庫包含了超過一個完整太陽活動周(11 a)的LEO 空間輻射帶數據，可以用于LEO 空間輻射帶建模.

表1 高能粒子探測儀器參數Table 1 High energy charged particle detector parameters

2 輻射帶建模的方法

本文利用我國輻射帶高能粒子的探測數據建立的自主輻射帶模型和NASA 的AP8/AE8 均為經驗模型，即模型的基礎是具有特定組織結構的數據，文中將這種具有特定組織結構的數據定義為表單.在輻射帶模型中，表單指的是能描述輻射環境某一給定時刻、給定空間環境狀態下在所有空間維度上分布特征的粒子通量數據列表.AP8/AE8僅含有一組代表太陽活動高年和低年時對應高能粒子通量平均分布的表單，不能給出輻射帶隨地磁場長期變化而引起的緩變，也不能反映太陽活動周內的各種周期性變化和各種突發性擾動變化［16］，因此稱其為靜態模型;代替它的AP9/AE9 只能通過輸出結果中的概率間接地反映各種變化對輻射帶的影響［17］，說明AP9/AE9也不是真正意義上的動態模型.

本文試圖建立的自主輻射帶模型的建模主要步驟如圖1 所示.從圖1 中可以看出，模型輸出結果是通過對模型數據中的4 組表單進行查詢、計算并將結果進行疊加后得到的，因此這4 組表單是模型的核心內容.這4 組表單分別給出輻射帶的寧靜狀態、長期緩變狀態、周期變化狀態和擾動變化狀態，因此自主輻射帶模型屬于動態模型，這也是該模型相比于AP8/AE8 的最大改進之處.表單的建立流程如圖2 所示.同化后的數據在空間中離散分布，必須先經過空間網格化處理，使數據在空間中均勻分布后，再進行分類處理，生成各類數據庫.因此，空間網格化是建模中的重要步驟.3 種軌道空間中，LEO 空間衛星最多，且LEO 空間包含了多個不同高度的軌道，而MEO(特指高度為22 000 km、傾角為55°的衛星軌道)空間和GEO 空間均僅包含單一軌道的若干顆衛星.因此，LEO 空間的數據網格化成為建模中的重點.

圖1 自主輻射帶模型框架Fig.1 Framework of own radiation belts model

圖2 自主輻射帶模型表單建立流程Fig.2 Procedure of creating data table in own radiation belts model

3 LEO 空間數據網格化

3.1 空間數據網格

鑒于LEO 空間觀測數據的分布特點，為了得到如圖3 所示的空間數據網格，本文采取了以下處理方法進行數據網格化:

1)利用11 個軌道高度上的衛星觀測資料，完成對應11 個空間球面上的數據網格化.對于同系列衛星觀測數據的使用原則是，利用其中部分衛星的觀測資料完成插值，剩余衛星的觀測資料用于對插值結果進行驗證.

2)11 個空間球面以外的空間，利用11 個空間球面的數據網格，結合大橢圓軌道衛星的觀測資料通過插值獲得;11 個空間球面之間的區域采用內插，高度范圍為300 ～340 km、860 ～1 000 km的區域采用外插.

圖3 LEO 空間數據網格Fig.3 Data grid of LEO space

3.2 可選用的插值方法

空間網格點上的粒子通量大小是通過對觀測數據插值獲得的，而插值中最關鍵的步驟是利用同一空間球面上觀測數據進行插值得到該球面上的網格點通量.常用的插值方法包括反距離加權法(Inverse Distance Weighting，IDW)［18］、克里格法(Kriging)［19］、自然鄰點法(natural neighbor)［20］、最近鄰點法(nearest neighbor)［21］、多元回歸法(polynomial regression)［22］、最小曲率法(minimum curvature)［23-24］、C1 連續5 次雙變量插值法(quintic)［25-26］、改進謝別德法(modified Shepard’s)［27-29］徑向基函數法(Radial Basis Function，RBF)［30-31］、三角形剖分法(triangulation)［32-34］等.這些方法在基本原理、數學模型、計算效率等方面各不相同，但在使用中一般都需要輸入輔助參數，例如:克里格法需要輸入所采用的變異函數的類型、變程、塊金值、基臺值［35］;改進謝別德法需要輸入節函數的兩個影響半徑［36］;徑向基函數法需要確定插值使用的徑向基函數［37］.

在應用插值方法時，由于網格點和觀測數據分布在球面上，因此插值形式應為球面插值;同時由于輻射帶內的粒子通量分布服從冪律譜［38］，因此計算中應采用對數插值.

使用空間球面上的網格點通量可以反演出該球面上任意位置的通量.通過對反演結果的誤差分析可以確定網格點通量的精度，進而確定選取的插值方法是否合適.對反演結果的誤差分析主要是對反演值(P)和觀測值(O)之差(D=P －Q)進行分析.誤差分析包括兩部分:計算各類誤差成分在總誤差中的比例和計算誤差大小［39］.本文采用均方誤差(Mean Square Error，MSE)來描述誤差比例.MSE 包括系統均方誤差(systematic MSE)MSEs和非系統均方誤差(unsystematic MSE)MSEu，其中MSEs由3 部分組成，包括MSE 中的偏移成分(additive MSE)MSEa、縮放成分(proportional MSE)MSEp和MSEa、MSEp互相依賴成分(interdependence MSE)MSEi.以上誤差的計算方法由式(1)～式(6)給出:

式中:N 為觀測點數;Pi和Oi分別為第i 點的計算值和觀測值;為觀測值的平均值;a 和b 分別為以O 為自變量、以P 為因變量采用最小二乘法擬合出的直線的截距和斜率.

本 文 采 用RMSE、RMSEs、RMSEu、RMSEa、RMSEp、RMSEi來描述誤差大小，這6 個參數值的量綱與粒子通量相同，大小分別等于MSE、MSEs、MSEu、MSEa、MSEp、MSEi的平方根，符號與MSE、MSEs、MSEu、MSEa、MSEp、MSEi相同.同時，還采用了一致性系數d 來描述P 和O 的一致程度，d的計算方法為

d 的取值范圍在0.0 ～1.0 之間，d 越接近1.0表示計算值和觀測值的一致性越高.

3.3 插值方法的應用

本文通過對軌道高度為630 km 的兩顆衛星的觀測數據的處理來說明空間球面上的數據網格化方法及精度評估標準.這兩顆衛星為同系列衛星，各自搭載了一臺由中國科學院空間中心空間環境探測研究室研制的高能粒子探測器，兩臺探測器的工作原理、制造工藝和技術指標完全相同.該探測器可以對6 個能道的高能質子通量和1 個能道的高能電子通量同時進行觀測.選定的觀測數據為2010 年8 月29 日至10 月10 日期間能量范圍為3 ～5 MeV 的質子通量，該時間段內空間環境始終保持寧靜，得到的觀測數據中不包含突發性空間環境擾動引起的質子通量變化.

數據處理的具體步驟如下:

1)利用其中一顆衛星的數據，選取一種插值方法，給定該方法所需的各項參數值，得到衛星所處空間球面上的網格點通量，網格劃分以經度、緯度為單位，每1°作為一個網格.

2)使用網格點通量反演另一顆衛星各觀測位置的通量，反演采用趨勢面插值法，即在網格點中選取3 個點，此3 個點形成的三角形是可以包圍采樣位置的最小三角形，將網格點的經度、緯度作為x、y 坐標，通量作為z 坐標，擬合出一個趨勢面方程，利用該方程計算出觀測位置的通量.

3)將反演值與觀測值進行比對，得到所選定插值方法在采用選定參數時的誤差分析結果.

在第3.2 節給出的各種插值方法中，反距離加權法、自然鄰點法、最近鄰點法、克里格法和徑向基函數法支持球面插值，克里格法和徑向基函數法在使用中需要較多的人工干預才能得到精度較高的結果，不適于工程應用.因此，在實際建模中本文選擇了反距離加權法、自然鄰點法和最近鄰點法進行對比.表2 給出了這3 種插值方法處理結果的誤差分析對比，改變反距離加權法的輸入參數會得到不同的處理結果，表2 給出的是其誤差分析最優的結果.從表2 給出的結果可以看出，反距離加權法在最優狀態下得到的處理結果的精度高于自然鄰點法和最近鄰點法，原因在于其輸入參數優化處理能夠大幅減小系統偏差(RMSEs)中的縮放成分(RMSEp)，自然鄰點法和最近鄰點法沒有可調整的輸入參數，反距離加權法有兩個輸入參數可以調整，其計算原理為

表2 采用反距離加權法、自然鄰點法和最近鄰點法處理結果的誤差分析對比Table 2 Error analysis comparison of processing results of IDW，natural neighbor and nearest neighbor

在計算中，將r 的范圍確定為1 ～22013.14 km，其中1 km 為觀測點至網格點距離的最小值，22 013.14 km為衛星軌道球面上大圓的半周長.在此范圍下，本文獲得了p 從1 ～10 的處理結果.圖4和圖5 所示分別為不同p 取值下的d 和RMSE 的對比.從圖4 可以看出，當r ＜1 000 km時，不同p 取值下的結果基本相同;當r ＞1000 km時，在p=1 和p=2 的情況下d 會隨著r 的增加迅速減小，在p≥3 的情況下d 基本保持不變，仍然維持較高水平.圖5 所示的RMSE 的變化規律與d 類似，在p=1、2 和3 的情況下最小RMSE 值均小于p≥4，但獲得最小RMSE 的r 各不相同.圖4和圖5 的結果說明，在p ≥4 的情況下，當r ＞65 km時，不同p 值下d 和RMSE 的變化規律幾乎完全相同，說明此種情況下單純增大r 不能影響插值結果的精度，這是由反距離加權法的原理決定的.從式(10)中可知，p 增大會增加距離網格點較近的實測點在插值結果中的權重.因此，當r 取值使得搜索范圍包含全部對插值結果有主要作用的實測點后，增大r 基本不會改變插值結果的精度.

圖4 反距離加權法取不同p 時處理結果的d 對比Fig.4 Comparison of processing results of d with different p in IDW

圖5 反距離加權法取不同p 時處理結果的RMSE 對比Fig.5 Comparison of processing results of RMSE with different p in IDW

圖6 給出了反演結果和觀測結果的對比，反演結果采用反距離加權法處理獲得，p 從1 取到10，選定的r 使當前p 下反演結果的RMSE 最小.

圖6 質子通量觀測結果和網格數據反演結果對比Fig.6 Proton flux comparison of observation and inversion result by data grid

從圖6 可以發現觀測結果具有以下缺陷:

1)本底較高，基本維持在大于10/(cm2·sr·s)的水平，這是由探測器幾何因子較小決定的.

2)通量范圍在10 ～100/(cm2·sr·s)的區間內，觀測結果抖動比較嚴重，這是受儀器本底噪聲和觀測數據采用遙測分層值方式采集的共同影響造成的.

理想的網格數據給出的反演結果，應該可以對觀測數據中出現的以上缺陷予以修正.對比圖6給出的各階反演結果，可以發現它們具有以下特點:

1)通量大于100/(cm2·sr·s)時，各階反演結果與觀測結果一致性均較高.

2)通量小于100/(cm2·sr·s)時，低階(p=1，p =2)反演結果具有較好的消除抖動的效果，對于本底以下的通量具有較好的外推效果.

因此，使用反距離加權法進行處理，當距離階數p 選定為1 或2 時，得到的反演結果合理程度更高;同時，根據圖4 和圖5 給出的誤差分析結果，選擇合適的搜索半徑r，可以使反演結果精度最高.

4 結 論

本文說明了自主輻射帶模型屬于動態模型的原因，指出了自主模型相比于AP8/AE8 的改進之處，即可以給出各種輻射帶周期性變化和擾動變化;討論了輻射帶自主建模的目標和思路，針對LEO 空間數據網格化問題，根據目前我國自主探測數據的實際情況，使用了空間分層的方式對每一個衛星所處空間球面上的數據分別進行處理;在生成數據網格的插值計算中，根據各常用插值方法的特點，選擇反距離加權法、自然鄰點法和最近鄰點法作為考察的插值方法，通過對網格數據反演結果的誤差分析和與觀測結果的比對，確定反距離加權法具有更高的精度，且反距離加權法在采用低階距離的情況下生成的數據網格得到的反演結果具有更高的合理性，并可以對實測數據中的一些缺陷進行彌補.

本文討論的LEO 空間數據網格化方法的評價標準建立在誤差分析的基礎上，是一種數值分析的方法.實際上，輻射帶的形成有著復雜的物理機制，其內部高能帶電粒子的分布也服從相應的物理規律，因此輻射帶建模不能脫離物理理論的指導.在后續工作中需要重點改進的地方是將物理理論與數值分析方法相結合，利用物理理論對數值分析方法進行完善.

References)

［1］ Bothmer V，Daglis I A.Space weather:Physics and effects［M］.New York:Springer，2007:154-159.

［2］ 徐文耀.地球電磁現象物理學［M］.合肥:中國科學技術大學出版社，2009:388-389.

Xu W Y.Physics of electromagnetic phenomena of the earth［M］.Heifei:University of Science and Technology of China Press，2009:388-389(in Chinese).

［3］ 沈自才.空間輻射環境工程［M］.北京:中國宇航出版社，2013:4-29.

Shen Z C.Space radiation environmental engineering［M］.Beijing:China Astronautic Publishing House，2013:4-29(in Chinese).

［4］ 姜景山，王文魁，都亨，等.空間科學與應用［M］.北京:科學出版社，2000:649-652.

Jiang J S，Wang W K，Du H，et al.Space science and applications［M］.Beijing:Science Press，2000:649-652(in Chinese).

［5］ 褚桂柏.空間飛行器設計［M］.北京:航空工業出版社，1996:49-54.

Chu G B.Space vehicle design［M］.Beijing:Aviation Industry Press，1996:49-54(in Chinese).

［6］ Ginet G P，O’Brien T P.Trapped radiation and plasma models requirements specification，AE-9/AP-9［R］.Washington，D.C.:NASA，2009.

［7］ Vette J I.The AE-8 trapped electron model environment，NASA STI/Recon Technical Report N，1991，92:24228［R］.Washington，D.C.:NSSDC/WDC-A-R＆S，1991.

［8］ Sawyer D M，Vette J I.AP-8 trapped proton environment for solar maximum and solar minimum，NASA STI/Recon Technical Report N，1976，77:18983［R］.Washington，D.C.:NSSDC/WDC-A-R＆S，1976.

［9］ 王世金，葉宗海.實踐四號衛星高能質子重離子探測器探測數據初步分析［J］.航天器工程，1995，4(3):1-7.

Wang S J，Ye Z H.Initial data analysis from the energetic proton and heavy ion detector on-board SJ-4 satellite［J］.Spacecraft Engineering，1995，4(3):1-7(in Chinese).

［10］ 師立勤，王世金，葉宗海，等.實踐五號高能粒子環境探測結果［C］∥實踐五號衛星空間探測與試驗成果學術會議論文集.北京:中國空間科學學會空間探測專業委員會，2002:51-61.

Shi L Q，Wang S J，Ye Z H.et al.Energetic particle environment detection results by SJ-5 satellite［C］∥Symposium of the Scholar Meeting on Space Exploration and Test Results by SJ-5 Satellite.Beijing:Space Exploration Committee of Chinese Society of Space Research，2002:51-61(in Chinese).

［11］ 王世金，朱光武，梁金寶，等.FY-1C 衛星空間粒子成分監測器及其探測結果［J］.上海航天，2001，18(2):24-28.

Wang S J，Zhu G W，Liang J B，et al.FY-1C space particle composition monitor and the results detected［J］.Aerospace Shanghai，2001，18(2):24-28(in Chinese).

［12］ Leonov A，Cyamukungu M，Cabrera J，et al.Pitch angle distribution of trapped energetic protons and helium isotope nuclei measured along the Resurs-01 No.4 LEO satellite［J］.Annales Geophysicae，2005，23(9):2983-2987.

［13］ 王鵬，徐青，李建勝，等.空間環境建模與可視化仿真技術［M］.北京:國防工業出版社，2012:92-95.

Wang P，Xu Q，Li J S，et al.Space environment modeling and visual simulation techniques［M］.Beijing:National Defence Industry Press，2012:92-95(in Chinese).

［14］ 王世金，朱光武，梁金寶.我國空間環境天基監測網建設［C］∥中國空間科學學會空間探測專業委員會第十六次學術會議論文集(上).北京:中國空間科學學會空間探測專業委員會，2003:65-68.

Wang S J，Zhu G W，Liang J B.Space-based space weather monitoring network building in China［C］∥Symposium of the 16th Scholar Meeting Held by Space Exploration Committee of Space Science Society in China(I).Beijing:Space Exploration Committee Chinese Society of Space Research，2003:65-68(in Chinese).

［15］ Friedel R H W，Bourdarie S，Cayton T E.Intercalibration of magnetospheric energetic electron data［J］.Space Weather，2005，3(9):2561-2570.

［16］ 莊洪春.宇航空間環境手冊［M］.北京:中國科學技術出版社，2000:140-148.

Zhuang H C.Aerospace environmental handbook［M］.Beijing:China Science and Technology Press，2000:140-148(in Chinese).

［17］ O’Brien T P.A framework for next-generation radiation belt models［J］.Space Weather，2005，3(7):1891-1904.

［18］ Serón F J，Badal J I，Sabadell F J.Spatial prediction procedures for regionalization and 3-D imaging of Earth structures［J］.Physics of the Earth and Planetary Interiors，2001，123(2):149-168.

［19］ Isaaks E H，Srivastava R M.Applied geostatistics［M］.New York:Oxford University Press，1989:279-322.

［20］ Watson D F.Contouring:A guide to the analysis and display of spatial data［J］.Computer Methods in the Geosciences，1993，19(10):1571-1572.

［21］ Aurenhammer F.Voronoi diagrams:A survey of a fundamental geometric data structure ［J］.ACM Computing Surveys(CSUR)，1991，23(3):345-405.

［22］ PIZOR P J.Principles of geographical information systems for land resources assessment［J］.Soil Science，1987，144(4):306.

［23］ Barrodale I，Skea D，Berkley M，et al.Warping digital images using thin plate splines［J］.Pattern Recognition，1993，26(2):375-376.

［24］ Powell M J D.Tabulation of thin plate splines on a very fine two-dimensional grid，DAMTP 1992/NA2［R］.Cambridge:University of Cambridge，1992.

［25］ Akima H.Algorithm 761:Scattered-data surface fitting that has the accuracy of a cubic polynomial［J］.ACM Transactions on Mathematical Software(TOMS)，1996，22(3):362-371.

［26］ Renka R J，Cline A K.A triangle-based C1interpolation method［J］.Rocky Mountain Journal，1984，14(1):223-238.

［27］ Franke R，Nielson G.Smooth interpolation of large sets of scattered data［J］.International Journal for Numerical Methods in Engineering，1980，15(11):1691-1704.

［28］ Renka R J.Algorithm 790:CSHEP2D:Cubic shepard method for bivariate interpolation of scattered data［J］.ACM Transactions on Mathematical Software (TOMS)，1999，25(1):70-73.

［29］ Shepard D.A two-dimensional interpolation function for irregularly-spaced data［C］∥Proceedings of the 1968 23rd ACM National Conference.New York:ACM，1968:517-524.

［30］ Franke R.A critical comparison of some methods for interpolation of scattered data，NPS-53-79-003［R］.Monterey，CA:Naval Postgraduate School，1979.

［31］ Hardy R L.Theory and applications of the multiquadric-biharmonic method 20 years of discovery 1968-1988［J］.Computers＆ Mathematics with Applications，1990，19(8-9):163-208.

［32］ Peucker T K，Fowler R J，Little J J，et al.The triangulated irregular network［C］∥International Symposium on Cartography and Computing:Applications in Health and Environment.Gaithersburg，Maryland:AmericanCongress on Surveying and Mapping，1978:516-532.

［33］ Krajewski S A，Gibbs B L.Understanding contouring:A practical guide to spatial estimation and contouring using a computer and basics of using variograms［M］.Boulder，CO:Gibbs Associates，2001:20-25.

［34］ Lee D T，Schachter B J.Two algorithms for constructing a Delaunay triangulation［J］.International Journal of Computer ＆Information Sciences，1980，9(3):219-242.

［35］ 侯景儒，尹鎮南，李維明，等.實用地質統計學［M］.北京:地質出版社，1998:45-50.

Hou J R，Yin Z N，Li W M，et al.Practical geology statistics［M］.Beijing:Geological Publishing House，1998:45-50(in Chinese).

［36］ 張維娜，呂云霄，吳美平.改進謝別德插值方法在地磁圖構建中的應用［J］.導航與控制，2011，10(1):28-32.

Zhang W N，Lv Y X，Wu M P.Application of modified shepard interpolation in geomagnetic map［J］.Navigation and Control，2011，10(1):28-32(in Chinese).

［37］ 魏義坤，楊威，劉靜.關于徑向基函數插值方法及其應用［J］.沈陽大學學報，2008，20(1):7-9.

Wei Y K，Yang W，Liu J.Interpolation method of radial basis function and its application［J］.Journal of Shenyang University，2008，20(1):7-9(in Chinese).

［38］ 涂傳詒.日地空間物理學:行星際與磁層.下冊［M］.北京:科學出版社，1988:118-127.

Tu C Y.Sun-earth space physics:Interplanetary and magnetosphere (II)［M］.Beijing:Science Press，1988:118-127(in Chinese).

［39］ Willmott C J.On the validation of models［J］.Physical Geography，1981，2(2):184-194.