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

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

引 言

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

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

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

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

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

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

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

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

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

2 數(shù)據(jù)來(lái)源

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

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

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

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

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

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

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

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

進(jìn)一步地，對(duì)2014 年和2018 年第1～10 年積日全部的TEC 數(shù)據(jù)進(jìn)行統(tǒng)計(jì)比較，給出兩種產(chǎn)品的偏差dTEC對(duì)比結(jié)果：

式 中：TECmad表 示Madrigal TEC 數(shù) 據(jù)；TECref表 示CODE 數(shù)據(jù)或DCT-PLS 算法重構(gòu)結(jié)果. TEC 偏差的統(tǒng)計(jì)結(jié)果如圖5所示. 可以看出，DCT-PLS 算法重構(gòu)的TEC 數(shù)據(jù)與Madrigal 數(shù)據(jù)的一致性和穩(wěn)定性有明顯增加. 從分析結(jié)果來(lái)看：CODE 電離層TEC 數(shù)據(jù)相對(duì)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)定性評(píng)對(duì)比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 算法的原因有兩個(gè)方面：1) CODE 采用15 階的球諧函數(shù)重構(gòu)全球TEC，數(shù)據(jù)產(chǎn)品的空間分辨率較低；2)對(duì)于同一區(qū)域，CODE 重構(gòu)TEC 地圖使用的GNSS 臺(tái)站數(shù)量也遠(yuǎn)低于本文DCT-PLS 算法.

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

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

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

本文基于DCT-PLS 算法，實(shí)現(xiàn)了區(qū)域高時(shí)空分辨率TEC 地圖的重構(gòu). 通過(guò)與2014 年和2018 年部分Madrigal 高精度TEC 數(shù)據(jù)的對(duì)比結(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. 對(duì)2017-09 一次磁暴期間的電離層重構(gòu)結(jié)果表明，本算法能夠較好地再現(xiàn)磁暴期間電離層精細(xì)化的擾動(dòng)結(jié)構(gòu)特征，相關(guān)研究結(jié)果可以為實(shí)現(xiàn)區(qū)域高分辨率電離層監(jiān)測(cè)和理論研究提供較好的參考.致 謝本文GNSS 觀測(cè)數(shù)據(jù)從國(guó)際GNSS 服務(wù)(international GNSS service, IGS)網(wǎng)站和中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(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/下載，作者在此表示感謝.