一種基于離散余弦變換-懲罰最小二乘回歸的區域高分辨率電離層TEC 地圖重構方法

歐明 陳龍江 吳家燕 甄衛民 呂夢海

（1. 中國電波傳播研究所，青島 266107；2. 中電科(青島)電波技術有限公司，青島 266107）

引 言

電離層是日地空間環境的重要組成，會對穿越其中的無線電波產生折射、反射、散射和吸收等效應，是衛星導航、通信、雷達等諸多無線電信息系統最嚴重的誤差來源之一. TEC 是表征電離層等離子體數密度的一個重要參數[1]，掌握其時空變化特征和規律，對提升無線電信息系統性能具有重要的實用價值，對空間科學研究及空間天氣預報也同樣具有重要的科學意義[2-4].

電離層觀測手段通常是離散不均勻分布的，獲取整個空間高精度電離層TEC 地圖一般需要通過特定的擬合建模或插值方法實現. 隨著衛星技術的發展，地基GNSS 已經成為全球TEC 地圖構建最為重要的數量來源. 1998-06 以來，在國際GNSS 服務(international GNSS service, IGS)電離層關聯分析中心五位成員的努力下，IGS 開始公開發布全球TEC 地圖(global ionospheric map, GIM)數據[5]. 多年以來，隨著GNSS 站網數量的增加及數據處理技術的發展，目前已有多種生成GIM 的建模策略，包括基于三角網格雙三次曲線插值法、球諧函數法、層析成像(computerized tomography, CT)法、球諧函數和廣義三角級數組合法、深度卷積生成對抗網絡法、Kalman濾波同化等多種算法，實現了全球TEC 地圖的重構并在線提供數據產品[6-9]. 現有的GIM 數據為分析電離層年度、季節、晝夜和半晝夜周期及磁緯度、經度等時空的大尺度變化提供了重要的支撐[10].

地震、海嘯、臺風等自然災害可能會導致小而明顯的電離層擾動[11-12]；太陽耀斑或日冕物質拋射(coronal mass ejections, CME)等空間天氣事件中短時間內電離層劇烈的擾動變化也時常發生[13-15]. 更小尺度電離層擾動的分析和建模亟需高時間和空間分辨率的電離層TEC 地圖產品. 但目前大多數GIM 數據產品基于球諧函數及其改進方法實現，其最大擴展階數為15，時間分辨率通常為1～2 h，地理緯度和經度網格大小為2.5°×5°，要實現1°×1°的空間分辨率，則需要高達180 階的球諧函數級數擴展[16]，這對現有全球GNSS 觀測臺站的分布而言是難以實現的.

隨著美國、歐洲、日本、中國等區域范圍內密集GNSS 臺站的建立，采用各類插值方法獲取高分辨率的區域TEC 地圖成為可能. 毛田等基于Kriging 插值算法實現了中國中低緯度TEC 地圖的重構，有效提高了電離層TEC 地圖的準確性[17]，但面對密集的GNSS 觀測數據，Kriging 算法存在計算量大的缺點.Aa 等使用國際參考電離層作為背景場，采用三維變分同化方法構建了中國及周邊區域電離層TEC 地圖[18]，但是數據同化技術實現過程復雜且需要占用較大內存，計算效率不高. 朱永興等提出一種顧及經緯度方向異性的電離層TEC 反距離加權插值算法，插值精度提升明顯，但反距離加權插值算法實現過程中設置的距離衰減指數需要依賴經驗獲得[19]. Goss 等基于多項式和三角B 樣條的張量積等局部化基函數，結合卡爾曼濾波技術，實現了近實時的高分辨率TEC 地圖建模[16]. Foster 等比較了自適應歸一化卷積(adaptive normalized convolution，ANC)、三次樣條插值、克里金(kriging) 法和自然鄰域法等多種TEC 地圖插值技術，分析結果表明：基于三次插值(Cubic)具有相對較低的復雜性和性能，許多情況下性能接近于其他更復雜的方法；Kriging 法在某些稀疏情況下表現良好，但偶爾會產生尖峰和大誤差；ANC 比Kriging 插值算法具有更好的誤差一致性；自然鄰域法具有出色的插值效果，但實現方式較為復雜[20].

目前，美國麻省理工學院(Massachusetts Institute of Technology, MIT)利用全球超過2 000 個GNSS 監測站，發布了1°×1°×5 min 分辨率的散點TEC 數據產品[21]；比利時皇家天文臺(Royal Observatory of Belgium, ROB)利用歐洲永久GNSS網絡(EUREF)非常 密 集 的 監 測 數 據，實 現 了 該 區 域0.5°×0.5°×15 min 分辨率的近實時TEC 地圖產品的在線發布(ftp://gnss.oma.be/ gnss/products)；國外美國JPL、日本NICT、德國DLR 及國內武漢大學、中科院等單位也正在開展高分辨率電離層TEC 地圖重構技術研究[18].

本文利用中國大陸構造環境監測網絡和IGS 的地基GNSS 監測數據，基于離散余弦變換-懲罰最小二乘回歸(discrete cosine transform and penalized least square regression, DCT-PLS)算法，實現了區域高時空分辨率TEC 地圖的重構. 與MIT Madrigal 數據庫高精度全球垂直TEC 數據的對比結果表明，DCTPLS 算法給出的TEC 誤差相比歐洲定軌中心(Center for Orbit Determination in Europe, CODE)的GIM 數據，一致性和穩定性有明顯提升，且能夠較好地給出磁暴期間電離層精細化的擾動結構特征.

1 TEC 地圖重構方法

電離層TEC 地圖重構引入了DCT-PLS 算法.DCT-PLS 算法最初被提出用于自動平滑多維不完整數據[22]，在此對其進行調整以填補地基GNSS 電離層穿刺點(ionospheric pierce point, IPP) TEC 數據集無法覆蓋區域的數據空白. PLS 是一種常用于一維數組的平滑器，該方法由Whittaker 于1923 年提出，它在數據的保真度與均值函數的粗糙度之間進行權衡. Garcia 證明了PLS 可以通過DCT 來表達，DCT 用在不同頻率振蕩的余弦函數的總和來表示數據[22]. 由于DCT 可以是多維的，因此基于DCT 的PLS 可以立即擴展到多維數據集.

DCT-PLS 算法通過使誤差E最小，計算地基GNSS觀測無法覆蓋網格點(數據空白區)的TEC 值，從而實現區域電離層TEC 地圖的重構.

2 數據來源

采用中國大陸構造環境監測網絡 (Crustal Movement Observation Network of China, CMONOC)和IGS 的地基GNSS 監測網2014 年和2018 年每年第1～10 年積日約290 個GNSS 臺站的觀測數據進行區域TEC.

首先利用Bernese 軟件對RINEX 格式文件中GNSS 的觀測數據進行周跳和異常值的探測，并采用雙頻的載波相位觀測數據平滑碼偽距，平滑處理后的觀測數據仍以RINEX 格式的文件輸出，這樣處理可較好剔除GNSS 數據中的部分異常觀測值[7]；然后利用意大利Gran Sasso 太空研究所研制的GNSS TEC 數據處理軟件對各臺站的RINEX 觀測數據進行處理，從而計算輸出各臺站觀測的GNSS 衛星編號、觀測時刻、衛星仰角、方位角、IPP 經緯度及傾斜路徑上的TEC 信息[24]. TEC 求解過程中，電離層球殼高度是400 km，截止仰角是10°，硬件延遲估計采用了內置的Arcs 方法.

通過計算衛星-接收機間的傾斜因子，傾斜TEC 可轉換為對應IPP 處的垂直TEC 值. 劃定緯度5～60°N、經度70～140°E 為電離層TEC 地圖重構區域. 圖1 所示為2014-01-01T04:00UT 地基GNSS IPP 及對應垂直電離層TEC 分布. 從圖1 中可以看出，IPP 基本覆蓋了設定的重構區域，但仍然存在不少的數據空白區域，這些區域的電離層TEC 值需要通過DCT-PLS 重構算法計算得到.

圖1 2014-01-01T04:00UT 地基GNSS IPP 及對應垂直電離層TEC 分布Fig. 1 The distribution map of GNSS ionospheric pierce points and corresponding vertical TEC at 04:00UT on January 1, 2014

選擇美國麻省理工大學Madrigal 數據庫中高精度垂直TEC 數據作為“真值”對TEC 地圖的重構結果進行評估；為了對比TEC 地圖的重構效果，DCTPLS 算法重構結果、CODE 的全球電離層TEC 也同時與Madrigal 數據進行比較分析.

3 TEC 地圖重構結果分析

3.1 TEC 地圖重構誤差分析

基于地基GNSS IPP 垂直TEC 數據，通過DCTPLS 算法對設定區域TEC 地圖進行重構，重構數據產品的空間地理經緯度分辨率為1°×1°，時間分辨率為15 min. 為顯示方便，圖2～4 分別按照1 h 時間間隔給出了2014-01-01T00:00―23:00UT 緯度5～60°N、經 度70～140°E 對 應 區 域Madrigal 高 精 度TEC、CODE 的全球電離層TEC 以及DCT-PLS 算法重構TEC 的時空分布特征. 從圖2～4 的對比可以看出：相比于Madrigal 的TEC 產品，CODE 給出的TEC 分布與觀測值基本形態較為符合，但整體分布過于平滑，忽略了絕大部分局部的TEC 變化特征；而DCT-PLS 算法重構的TEC 則與Madrigal 數據更為一致.

圖2 2014-01-01T00:00―23:00UT 指定區域Madrigal 垂直電離層TEC 分布Fig. 2 Distribution of Madrigal vertical TEC in corresponding region during 00:00―23:00UT on January 1, 2014

圖3 2014-01-01T00:00―23:00UT 指定區域CODE 全球電離層TEC 分布Fig. 3 Ionospheric TEC distribution of CODE-GIM in corresponding region during 00:00―23:00UT on January 1, 2014

圖4 2014-01-01T00:00―23:00UT 指定區域DCT-PLS 算法重構電離層TECFig. 4 Ionospheric TEC map reconstruction results in corresponding region using DCT-PLS algorithm during 00:00―23:00UT on January 1, 2014

進一步地，對2014 年和2018 年第1～10 年積日全部的TEC 數據進行統計比較，給出兩種產品的偏差dTEC對比結果：

式 中：TECmad表 示Madrigal TEC 數 據；TECref表 示CODE 數據或DCT-PLS 算法重構結果. TEC 偏差的統計結果如圖5所示. 可以看出，DCT-PLS 算法重構的TEC 數據與Madrigal 數據的一致性和穩定性有明顯增加. 從分析結果來看：CODE 電離層TEC 數據相對Madrigal 數據的平均誤差為3.9 TECU，標準差為3.7 TECU；而DCT-PLS 算法的平均誤差和標準差則分別下降為2.0 TECU 和2.7 TECU，平均誤差下降了約46%，標準差下降了約31%.

圖5 CODE 數據和DCT-PLS 算法重構的電離層TEC 地圖與Madrigal 數據一致性和穩定性評對比Fig. 5 Comparison of consistency and stability between TEC map by GIM and DCT-PLS algorithm and Madrigal data

應該指出的是，CODE 數據給出的TEC 平均誤差和標準差大于DCT-PLS 算法的原因有兩個方面：1) CODE 采用15 階的球諧函數重構全球TEC，數據產品的空間分辨率較低；2)對于同一區域，CODE 重構TEC 地圖使用的GNSS 臺站數量也遠低于本文DCT-PLS 算法.

3.2 TEC 地圖在磁暴期間電離層擾動分析中的應用

2017-09-07-08 期間，兩次CMEs 先后經過地球，并產生了具有雙主相的強地磁暴，第二主相期間，中國及鄰近區域出現了電離層擾動觀測結果[25]. 基于高分辨率電離層TEC 地圖的重構結果，圖6 給出了2017-09-08T12:00—17:45UT 期間對應區域的垂直電離層TEC 變化. 可以看出，從當地夜間大約13:00UT左右開始，TEC 耗空開始出現并覆蓋了廣泛的區域；然后向西北方向擴展，耗空區域主要在15～45°N 和80～110°E 區域以內，持續近5 h(約13:00―18:00UT)，TEC 耗空的深度在5～15 TECU，耗空達到的最高緯度約為50°N(磁緯約45.5°N). 進一步地，選擇緯度25°N 進行分析，基于高分辨率TEC 構建了一個TEC 概要圖(Keogram)[26]，以便于更詳細地顯示磁暴期間(12:00―19:00UT)垂直電離層TEC 的變化，具體結果如圖7 所示. 圖7 中TEC 耗空區域由黑色虛線突出顯示，結合圖6 分析發現TEC 耗空區能以較完好的兩個單獨的平行結構出現，TEC 耗空呈現從東南到西北的變化趨勢，估算平均緯向漂移速度約105 m/s，相關結果與文獻[26]基本一致. 說明本文的高分辨率TEC 地圖重構算法可以有效應用于中小尺度的電離層擾動分析中.

圖6 2017-09-08T12:00―17:45UT 磁暴期間指定區域高分辨率垂直電離層TEC 變化Fig. 6 Variation of vertical TEC in corresponding region during 12:00-17:45UT on September 8, 2017

圖7 2017-09-08T12:00-19:00UT 磁暴期間指定區域25°N 垂直電離層TEC Keogram 變化Fig. 7 TEC Keogram variation analysis of 25°N in corresponding region during 12:00—19:00UT on September 8, 2017

4 小 結

為滿足更小尺度電離層擾動的分析和建模的需求，亟需獲得高時間和高空間分辨率的電離層TEC 地圖產品. 隨著區域范圍內密集的GNSS 臺站的建立，采用各類插值方法獲取高分辨率的區域TEC 地圖成為可能.

本文基于DCT-PLS 算法，實現了區域高時空分辨率TEC 地圖的重構. 通過與2014 年和2018 年部分Madrigal 高精度TEC 數據的對比結果表明，DCTPLS 算法給出的垂直TEC 相比CODE 的GIM 數據的一致性和穩定性有明顯提升，其中TEC 平均誤差由3.9 TECU 下降為2.0 TECU，標準差由3.7 TECU下降為2.7 TECU. 對2017-09 一次磁暴期間的電離層重構結果表明，本算法能夠較好地再現磁暴期間電離層精細化的擾動結構特征，相關研究結果可以為實現區域高分辨率電離層監測和理論研究提供較好的參考.致 謝本文GNSS 觀測數據從國際GNSS 服務(international GNSS service, IGS)網站和中國大陸構造環境監測網絡(Crustal Movement Observation Network of China，CMONOC)獲取. 全球電離層地圖CODE數據從http://ftp.aiub.unibe.ch/CODE/網站獲取. 麻省理工學院高精度Madrigal 數據由http:// cedar.openmadrigal.org/下載，作者在此表示感謝.