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

歐明 陳龍江 吳家燕 甄衛(wèi)民 呂夢海

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

引 言

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

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

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

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

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

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

1 TEC 地圖重構(gòu)方法

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

DCT-PLS 算法通過使誤差E最小，計算地基GNSS觀測無法覆蓋網(wǎng)格點(數(shù)據(jù)空白區(qū))的TEC 值，從而實現(xiàn)區(qū)域電離層TEC 地圖的重構(gòu).

2 數(shù)據(jù)來源

采用中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò) (Crustal Movement Observation Network of China, CMONOC)和IGS 的地基GNSS 監(jiān)測網(wǎng)2014 年和2018 年每年第1～10 年積日約290 個GNSS 臺站的觀測數(shù)據(jù)進行區(qū)域TEC.

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

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

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

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

3 TEC 地圖重構(gòu)結(jié)果分析

3.1 TEC 地圖重構(gòu)誤差分析

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

圖2 2014-01-01T00:00―23:00UT 指定區(qū)域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 指定區(qū)域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 指定區(qū)域DCT-PLS 算法重構(gòu)電離層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 數(shù)據(jù)進行統(tǒng)計比較，給出兩種產(chǎn)品的偏差dTEC對比結(jié)果：

式 中：TECmad表 示Madrigal TEC 數(shù) 據(jù)；TECref表 示CODE 數(shù)據(jù)或DCT-PLS 算法重構(gòu)結(jié)果. TEC 偏差的統(tǒng)計結(jié)果如圖5所示. 可以看出，DCT-PLS 算法重構(gòu)的TEC 數(shù)據(jù)與Madrigal 數(shù)據(jù)的一致性和穩(wěn)定性有明顯增加. 從分析結(jié)果來看：CODE 電離層TEC 數(shù)據(jù)相對Madrigal 數(shù)據(jù)的平均誤差為3.9 TECU，標(biāo)準(zhǔn)差為3.7 TECU；而DCT-PLS 算法的平均誤差和標(biāo)準(zhǔn)差則分別下降為2.0 TECU 和2.7 TECU，平均誤差下降了約46%，標(biāo)準(zhǔn)差下降了約31%.

圖5 CODE 數(shù)據(jù)和DCT-PLS 算法重構(gòu)的電離層TEC 地圖與Madrigal 數(shù)據(jù)一致性和穩(wěn)定性評對比Fig. 5 Comparison of consistency and stability between TEC map by GIM and DCT-PLS algorithm and Madrigal data

應(yīng)該指出的是，CODE 數(shù)據(jù)給出的TEC 平均誤差和標(biāo)準(zhǔn)差大于DCT-PLS 算法的原因有兩個方面：1) CODE 采用15 階的球諧函數(shù)重構(gòu)全球TEC，數(shù)據(jù)產(chǎn)品的空間分辨率較低；2)對于同一區(qū)域，CODE 重構(gòu)TEC 地圖使用的GNSS 臺站數(shù)量也遠低于本文DCT-PLS 算法.

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

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

圖6 2017-09-08T12:00―17:45UT 磁暴期間指定區(qū)域高分辨率垂直電離層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 磁暴期間指定區(qū)域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 小 結(jié)

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

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