鄭敦勇, 姚宜斌, 聶文鋒, 林東方, 梁繼, 陳春花

1 湖南科技大學(xué)地理空間信息技術(shù)國家地方聯(lián)合工程實(shí)驗(yàn)室, 湖南湘潭 411201 2 湖南科技大學(xué)地球科學(xué)與空間信息工程學(xué)院, 湖南湘潭 411201 3 武漢大學(xué)測繪學(xué)院, 武漢 430079 4 廣西空間信息與測繪重點(diǎn)實(shí)驗(yàn)室, 桂林 541006 5 山東大學(xué)空間研究院, 山東威海 264209 6 湖南省測繪科學(xué)技術(shù)研究所, 長沙 410007

0 引言

電離層位于地面上約60～1000 km范圍內(nèi)，是地球高空大氣受太陽高能輻射以及宇宙線的激勵(lì)而被電離的區(qū)域，主要由帶電粒子、中性原子及分子組成，在整個(gè)日地空間環(huán)境中，電離層是直接影響人類空間活動(dòng)的最重要環(huán)節(jié)之一(Liu et al. 2021).GNSS電離層層析技術(shù)(Computerized ionospheric tomography, CIT)作為監(jiān)測和研究電離層的一種重要手段，其重構(gòu)的電子密度分布能同時(shí)反映電離層水平結(jié)構(gòu)和垂直結(jié)構(gòu)的變化，從而克服了二維電離層模型的局限性，因此，在近20余年來得到了突飛猛進(jìn)的發(fā)展.目前，這些模型大致可以分為兩類：一類是函數(shù)基電離層層析模型(Hansen et al., 1997)；另一類是像素基電離層層析模型(Rius et al., 1997; 徐繼生和鄒玉華, 2003).相比之下，函數(shù)基CIT模型可以用少量的模型參數(shù)去反演大范圍的電離層電子密度，且反演結(jié)果較穩(wěn)定，但缺點(diǎn)是反演結(jié)果過于平滑，有時(shí)會(huì)掩蓋一些電離層的小尺度擾動(dòng)結(jié)構(gòu)特征.而在像素基CIT模型的反演過程中，也存在一個(gè)比較棘手的問題：可利用的觀測數(shù)據(jù)不足或分布不均，從而導(dǎo)致待反演的電離層空間離散化后部分像素沒有信號(hào)射線穿過，這也就意味著電離層電子密度的反演問題是一個(gè)不適定問題(聞德保, 2013).為了克服上述問題給像素基CIT模型帶來的不利影響，國內(nèi)外許多學(xué)者提出了解決方法，其主要分為算法改進(jìn)、多源數(shù)據(jù)融合以及神經(jīng)網(wǎng)絡(luò)模型等三類.

針對(duì)層析算法的改進(jìn)，目前主要包括基于奇異值分解(Singular Value Decomposition, SVD)、截?cái)嗥娈愔捣纸?Truncated SVD, TSVD)的非迭代算法，基于代數(shù)重建算法(Algebraic Reconstruction Technique, ART)、乘法代數(shù)重建算法(Multiplicative ART, MART)、同步迭代重構(gòu)算法(Simultaneous Iterative Reconstruction Technique, SIRT)、Kalman濾波以及正則化的改進(jìn)迭代算法，其中，非迭代算法的優(yōu)點(diǎn)是不需要初值，但由于可用的觀測數(shù)據(jù)不足以及離散誤差和觀測噪聲的影響，若單獨(dú)使用會(huì)導(dǎo)致其反演結(jié)果精度不高(聞德保, 2013)，因此，目前主要用于顧及數(shù)據(jù)權(quán)比因子的多源數(shù)據(jù)融合反演(湯俊, 2014)，另外，非迭代算法一般都涉及大規(guī)模矩陣求逆，從而對(duì)硬件要求也較高.

相較于非迭代算法，目前在像素基CIT研究中使用改進(jìn)迭代算法相對(duì)較多.如Yavuz等(2005)提出了聯(lián)合整體最小二乘法和ART的混合算法，其中利用整體最小二乘法為ART算法提供初值；肖宏波等(2008)發(fā)展了一種聯(lián)合最速下降法和MART的改進(jìn)算法，該算法有效降低了MART算法對(duì)迭代初值的敏感度，提高了電離層層析的可靠性；Liu等(2010)提出了附加約束的同步迭代代數(shù)重構(gòu)(Constrained CSIRT, CSIRT)算法；Kunitsyn等(2011)提出了基于Sobolev正則化約束的同步迭代重構(gòu)算法；聞德保(2013), 聞德保等(2014)以及Wen等(2015)先后提出了改進(jìn)ART算法、基于Phillips平滑約束和MART的二步算法、基于TSVD和ART的組合算法、約束MART算法、約束ART算法、約束自適應(yīng)聯(lián)合迭代重構(gòu)算法以及選權(quán)擬合重構(gòu)算法；Yao等(2013a, b, 2018)先后提出了基于函數(shù)基與像素基模型的組合算法、混合正則化迭代算法以及附加側(cè)面射線的CIT算法；湯俊(2014)提出了基于多源數(shù)據(jù)融合的CIT算法以及基于總變差最小化的CIT算法；霍星亮等(2016)通過改進(jìn)迭代算法及其松弛因子來提高CIT的反演精度，相對(duì)于傳統(tǒng)ART算法，該方法反演的電子密度剖面更接近于實(shí)測數(shù)據(jù)；Wang等(2016)提出一種基于Tikhonov正則化的CIT方法，利用模型函數(shù)方法去確定最優(yōu)的正則化參數(shù)，以此來平衡觀測數(shù)據(jù)和背景模型值之間的權(quán)重，并提高算法的全局收斂速度；Yin等(2017)發(fā)展了多程層析迭代算法，該算法通過構(gòu)建不同的像素分辨率來實(shí)現(xiàn)多次反演，從而提高反演精度；趙海山等(2018)提出自適應(yīng)層析算法(Adaptive MART, AMART)，其中根據(jù)有信號(hào)射線像素中的截距和電子密度值的綜合影響，合理分配算法的迭代差值，并提出與電子密度值相關(guān)的自適應(yīng)松弛因子，有效克服傳播噪聲對(duì)電子密度反演的影響；Jin和Li(2018)提出了改進(jìn)的二步算法，首先用NeQuick模型對(duì)背景模型進(jìn)行精化，然后結(jié)合MART算法進(jìn)行反演；Chen等(2019)提出了一種參數(shù)化電離層層析方法，根據(jù)像素八個(gè)角處的電子密度值(Ionospheric electron density, IED)，利用垂直指數(shù)插值和水平反向距離加權(quán)插值獲取像素內(nèi)任意位置的IED；除此之外，Zheng等(2015, 2016, 2017, 2018, 2021b)也先后提出了基于多尺度剖分、可變像素高度、可變像素尺度、松弛因子自動(dòng)搜索技術(shù)以及有序子集約束ART的電離層層析算法，均表現(xiàn)出良好的層析效果和效率.總的來說，通過增加約束條件或組合其他算法得到的這類改進(jìn)迭代算法，其性能較之前有大幅度提高，但離完全克服不適定性問題仍有較大差距，主要表現(xiàn)在無信號(hào)射線像素的反演電子密度對(duì)初值比較依賴，而初值通常是由精度不高的經(jīng)驗(yàn)?zāi)Ｐ徒o出，從而拉低了電離層層析的整體精度.

Stolle等(2003)在基于GPS數(shù)據(jù)的CIT反演中，提出融合掩星數(shù)據(jù)來提高電子密度的垂直分辨率；Bust等(2004)基于三維變分同化技術(shù)3DVAR，發(fā)展了三維電離層數(shù)據(jù)同化算法(IDA3D)，IDA3D可以同化GNSS、低軌信標(biāo)以及測高儀等多類觀測數(shù)據(jù)；Garca-Fernández等(2005)利用日本密集的GPS網(wǎng)數(shù)據(jù)，融合SAC-C掩星和測高儀數(shù)據(jù)，反演了日本上空的三維電離層電子密度，其中掩星和測高儀數(shù)據(jù)以及改進(jìn)的Abel算法被用來構(gòu)造電子密度剖面的背景分布，以此來提高電子密度的垂直分辨率；徐繼生等(2005)建立了融合GPS和掩星數(shù)據(jù)的時(shí)變CIT模型，并利用實(shí)測數(shù)據(jù)對(duì)模型進(jìn)行了仔細(xì)驗(yàn)證；Lee等(2007)融合GPS和測高儀數(shù)據(jù)，基于Tikhonov和總變分TV算法，得到兩種不同的CIT模型，并利用測高儀數(shù)據(jù)和Chapman函數(shù)構(gòu)造電子密度垂直剖面模型來驗(yàn)證反演精度；Angling(2008)基于最優(yōu)化線性無偏估計(jì)理論，融合GPS和掩星數(shù)據(jù)，發(fā)展了電離層電子密度同化模型；Chartier等(2012)發(fā)展了融合GPS和電離層測高儀數(shù)據(jù)的CIT方法，其中利用測高儀數(shù)據(jù)構(gòu)建垂直方向的基函數(shù)，自適應(yīng)地在反演中提供電子密度信息；Li等(2012)發(fā)展了一種融合地基GPS和COSMIC掩星數(shù)據(jù)的CIT方法；Zhao等(2013)提出了一種融合三頻GPS、地面垂測和斜測數(shù)據(jù)以及頂部探測儀數(shù)據(jù)的CIT方法；Yu等(2014)發(fā)展了電離層數(shù)據(jù)同化系統(tǒng)IDAAS，并對(duì)東南亞上空的電離層三維電子密度結(jié)構(gòu)進(jìn)行了重構(gòu)，被同化的數(shù)據(jù)包括地基GPS和測高儀數(shù)據(jù)；Semeter等(2016)提出了融合GPS和非相干散射雷達(dá)ISR數(shù)據(jù)的CIT方法，其中把插值的ISR密度場作為反演初值，然后使用最大熵正則化方法求解；Aa等(2016)利用三維變分?jǐn)?shù)據(jù)同化算法以及地基GPS和掩星數(shù)據(jù)，得到中國及相鄰區(qū)域上空電離層的三維電子密度結(jié)構(gòu)，其結(jié)果證明了該同化技術(shù)在提高區(qū)域電離層形態(tài)規(guī)范方面的有效性；歐明等(2016)提出了一種融合天地基多源數(shù)據(jù)的CIT算法，其多源數(shù)據(jù)包括地基GNSS、LEO掩星、衛(wèi)星信標(biāo)及測高儀等數(shù)據(jù)，并利用改進(jìn)的Kriging插值與MART算法實(shí)現(xiàn)多源數(shù)據(jù)的有效融合；She等(2017)提出了一種基于多源數(shù)據(jù)同化技術(shù)的CIT算法，該方法在水平和垂直方向上分別采用球諧函數(shù)和經(jīng)驗(yàn)正交基函數(shù)來進(jìn)行電子密度估計(jì)，其同化數(shù)據(jù)包括地基GPS、LEO掩星、海洋測高衛(wèi)星數(shù)據(jù)等；dos Santos Prol等(2019)提出融合測高儀、掩星、以及垂直方向總電子含量(Vertical Total Electron Content, VTEC)數(shù)據(jù)的層析方法，其中利用測高儀和掩星數(shù)據(jù)對(duì)背景模型進(jìn)行精化，然后把通過投影函數(shù)得到的VTEC數(shù)據(jù)加入到層析方程中，從而提高層析精度.盡管上述研究表明多源數(shù)據(jù)融合技術(shù)的發(fā)展使陸上大范圍的層析精度有了一定提高，但仍存在兩個(gè)主要問題：①部分研究中，直接將各類觀測數(shù)據(jù)按等權(quán)處理，這顯然是不合理的，而另一部分研究中，利用方差分量估計(jì)對(duì)觀測數(shù)據(jù)的權(quán)重進(jìn)行了確定，但其過程略顯復(fù)雜，且實(shí)施起來有一定難度，如很多學(xué)者利用Helmert方差分量估計(jì)法來確定權(quán)重，但在實(shí)際應(yīng)用時(shí)，各類方差因子相等的收斂條件是難以滿足的；②由于受時(shí)空分辨率的限制，多源數(shù)據(jù)融合技術(shù)并不能完全克服不適定性問題，特別是針對(duì)陸上的小范圍層析以及海域?qū)游鰰r(shí)，可利用的多源數(shù)據(jù)不多甚至缺失，因此，不適定性問題依然會(huì)較嚴(yán)重.

神經(jīng)網(wǎng)絡(luò)以其自身的自組織、自適應(yīng)和自學(xué)習(xí)等特點(diǎn)被廣泛應(yīng)用于各個(gè)領(lǐng)域，Ma等(2005)首次提出利用BP神經(jīng)網(wǎng)絡(luò)重構(gòu)區(qū)域電離層電子密度的思想，并融合測高儀數(shù)據(jù)來提高電子密度垂直分辨率；陳必焰等(2012)基于傳統(tǒng)層析模型的反演結(jié)果，利用神經(jīng)網(wǎng)絡(luò)對(duì)電子密度進(jìn)行預(yù)測，在建模過程中，用當(dāng)前反演的全部電子密度及其像素特征構(gòu)建學(xué)習(xí)樣本，然后用訓(xùn)練好的模型預(yù)測下一時(shí)段的電子密度，其預(yù)測精度明顯高于IRI(International Reference Ionosphere)模型；Razin和Voosoghi(2016)設(shè)計(jì)了一種多隱含層神經(jīng)網(wǎng)絡(luò)CIT算法，并基于IRI2012模型構(gòu)造了一個(gè)經(jīng)驗(yàn)正交基函數(shù)作為神經(jīng)網(wǎng)絡(luò)訓(xùn)練的目標(biāo)函數(shù)，以此來提高電離層垂直方向的分辨率，在此基礎(chǔ)上，Razin和Voosoghi(2019)發(fā)展了基于粒子群優(yōu)化訓(xùn)練的小波神經(jīng)網(wǎng)絡(luò)CIT算法，該算法不僅能改善層析模型的不適定問題，且能保證效率和全局收斂.盡管這些基于神經(jīng)網(wǎng)絡(luò)的層析算法對(duì)促進(jìn)CIT技術(shù)的發(fā)展起到了一定的推動(dòng)作用，但仍存在三個(gè)主要問題：①把精度不高的經(jīng)驗(yàn)?zāi)Ｐ?IRI)值當(dāng)做真值去構(gòu)建學(xué)習(xí)樣本，這勢必會(huì)導(dǎo)致最終神經(jīng)網(wǎng)絡(luò)模型的預(yù)測精度也不高；②用傳統(tǒng)層析算法的全部反演結(jié)果構(gòu)建學(xué)習(xí)樣本，沒有顧及那些無射線穿過像素的反演電子密度對(duì)初值的依賴問題；③學(xué)習(xí)樣本的輸入?yún)?shù)和時(shí)間跨度沒有考慮全面.鑒于此，Zheng等(2021a)提出了基于AdaBoost-BPNN算法的電離層層析模型，該模型利用有信號(hào)射線像素的特征及其電子密度反演值對(duì)機(jī)器學(xué)習(xí)模型進(jìn)行訓(xùn)練，然后預(yù)測無信號(hào)射線像素的電子密度，有效避免傳統(tǒng)層析算法中仍然存在的不適定性對(duì)層析結(jié)果的影響，進(jìn)而提高電離層層析的整體精度.然而在該模型中，仍然存在未完全顧及層析初值的影響，且其機(jī)器學(xué)習(xí)算法有待進(jìn)一步優(yōu)化等問題，為此，本文提出使用機(jī)器學(xué)習(xí)集成算法XGBoost對(duì)層析初值進(jìn)行預(yù)測，最終建立一種高精度的機(jī)器學(xué)習(xí)電離層層析模型(簡稱XGB-CIT模型)，并以獨(dú)立的電離層測高儀數(shù)據(jù)為參考值，對(duì)XGB-CIT模型進(jìn)行有效驗(yàn)證.

1 XGB-CIT模型構(gòu)建

GNSS信號(hào)路徑上的總電子含量(Total Electron Content, TEC)可表示為電子密度沿信號(hào)傳播路徑的線性積分，在電離層層析過程中，首先利用相位平滑偽距方法(Jin et al., 2012)獲取信號(hào)路徑上的TEC，然后對(duì)上述線性積分進(jìn)行離散化后并利用相應(yīng)的反演算法獲取其像素的IED(聞德保, 2013).

如Zheng等(2018)所示的改進(jìn)IART算法在本文中被使用，一方面可以為機(jī)器學(xué)習(xí)模型提供輸入數(shù)據(jù)，另一方面把機(jī)器學(xué)習(xí)模型預(yù)測的電子密度作為初值，然后反演得到最終的電離層電子密度，其具體的流程圖如圖2所示，同時(shí)IART算法還可以作為參考方法進(jìn)行對(duì)比分析.本文中對(duì)IART算法施加水平和垂直平滑約束，以避免相鄰像素之間電子密度值的突然變化，其中水平約束采用高斯加權(quán)函數(shù)構(gòu)造(Wen et al., 2015)，垂直約束利用Chen等(2019)中的方法，并結(jié)合機(jī)器學(xué)習(xí)模型預(yù)測的電子密度得到.

1.1 XGBoost算法

集成學(xué)習(xí)(Ensemble learning)是通過構(gòu)建并結(jié)合多個(gè)個(gè)體學(xué)習(xí)器來完成學(xué)習(xí)任務(wù)(圖1)，從而獲得比單一學(xué)習(xí)器顯著優(yōu)越的泛化性能(周志華, 2016).根據(jù)個(gè)體學(xué)習(xí)器的生成方式，目前的集成學(xué)習(xí)方法大致可分為兩大類：基于Boosting的和基于Bagging的，前者的代表算法有AdaBoost、GBDT、XGBoost、LightGBM以及CatBoost，后者的代表算法主要是隨機(jī)森林(RF).通過對(duì)比分析，本文選擇XGBoost算法并使用Python軟件來實(shí)現(xiàn)電離層層析重構(gòu).

圖1 集成學(xué)習(xí)示意圖Fig.1 The schematic diagram of ensemble learning

XGBoost(Extreme Gradient Boosting)是由Chen和Guestrin于2011年首次提出的機(jī)器學(xué)習(xí)模型，并在眾多學(xué)者的后續(xù)研究中不斷優(yōu)化和改進(jìn)(Chen and Guestrin, 2016; Li et al., 2019).XGBoost對(duì)損失函數(shù)進(jìn)行二階泰勒展開，可以自動(dòng)使用CPU的多線程進(jìn)行并行計(jì)算，且使用多種方法避免過度擬合.下面簡要介紹XGBoost算法，更多詳細(xì)信息可參考文獻(xiàn)Chen和Guestrin (2016).

XGBoost中用加法將所有樹模型集成，假設(shè)總共有K個(gè)樹，則有

(1)

q(x)表示將樣本x分到了某個(gè)葉子節(jié)點(diǎn)上，在本文中，針對(duì)如表1所示的樣本數(shù)據(jù)，若我們以輸入變量Lon<112.75°E為分裂點(diǎn)，則可得到左右兩個(gè)葉子節(jié)點(diǎn)，其中左葉子節(jié)點(diǎn)中的樣本有3，右葉子節(jié)點(diǎn)中的樣本有1、2、4、5；w是葉子節(jié)點(diǎn)的權(quán)重，在本文中即為電子密度預(yù)測值，最終根據(jù)(5)式計(jì)算得到.

表1 樣本數(shù)據(jù)示例Table 1 Sample data example

XGBoost的目標(biāo)函數(shù)如下：

(2)

(3)

Ω是防止過擬合以及控制模型復(fù)雜度的正則化函數(shù)：

(4)

其中T是樹的葉子總數(shù)，γ和λ分別為L1和L2正則的懲罰項(xiàng).為了使得目標(biāo)函數(shù)最小，正則項(xiàng)越小越好，因此，一方面要使葉子節(jié)點(diǎn)個(gè)數(shù)T最少，這樣樹就簡單，從而提高計(jì)算效率，實(shí)際應(yīng)用時(shí)，T由樹的深度max_depth決定；另一方面要對(duì)葉子節(jié)點(diǎn)的值w進(jìn)行L2正則，這樣才會(huì)把w的解約束在一個(gè)范圍內(nèi).

通過對(duì)目標(biāo)函數(shù)(4)的改寫、泰勒展開以及求導(dǎo)，可分別得到第j個(gè)葉子的最優(yōu)權(quán)重(式(5))和相應(yīng)的最優(yōu)目標(biāo)函數(shù)值(式(6))：

(5)

(6)

這里Ij={i|q(xi)=j}表示i樣本落在第j個(gè)葉子節(jié)點(diǎn)上，

通常不可能枚舉所有可能的樹結(jié)構(gòu)q，取而代之的是一種貪婪算法，該算法迭代地添加樹，一直到K顆樹停止，實(shí)際應(yīng)用中，弱分類器數(shù)目n_estimators即為樹的顆數(shù)，具體迭代過程如下：

(7)

另外，假設(shè)IL和IR是拆分后左右節(jié)點(diǎn)的實(shí)例集，設(shè)I=IL∪IR，則分割后的損失函數(shù)由下式給出：

(8)

1.2 XGB-CIT模型流程

如圖2所示，XGB-CIT的流程大致可分為3部分.

圖2 XGB-CIT流程圖Fig.2 The flow chart of XGB-CIT

(1)訓(xùn)練樣本構(gòu)建

基于傳統(tǒng)層析IART算法，利用多個(gè)連續(xù)時(shí)段(滑動(dòng)窗口內(nèi))中有信號(hào)射線像素的特征及其電子密度反演值(IED1)構(gòu)建XGBoost的訓(xùn)練樣本，其輸出參數(shù)即為IED1，輸入?yún)?shù)來源于與電子密度變化相關(guān)的幾類變量，結(jié)合Zheng等(2021a)的介紹以及大量試驗(yàn)，本文最終確定12個(gè)輸入?yún)?shù)，包括像素中心的緯度(Lat)、經(jīng)度(Lon)以及像素中心距地面的高度H，世界時(shí)的正弦(UTS)和余弦(UTC)分量(sin(2π·UT/24), cos(2π·UT/24))，太陽仰角的正弦(ELES)和余弦(ELEC)分量(sin(2π·ELE/360), cos(2π·ELE/360))，太陽方位角的正弦(AZIS)和余弦(AZIC)分量(sin(2π·AZI/360), cos(2π·AZI/360))，地磁指數(shù)Dst、ap和AE值.

(2)機(jī)器學(xué)習(xí)模型訓(xùn)練

目標(biāo)函數(shù)是判斷機(jī)器學(xué)習(xí)模型質(zhì)量的基礎(chǔ)，因此構(gòu)造合理的目標(biāo)函數(shù)是XGB-CIT成功的關(guān)鍵，本文的總目標(biāo)函數(shù)(total objective function, OF)定義如下：

OF=model_gs.best_score+rmse_score,

(9)

式中model_gs.best_score與rmse_score分別是交叉檢驗(yàn)的最優(yōu)得分和檢驗(yàn)樣本(未參與訓(xùn)練)的TEC正演誤差.

由于XGBoost在使用時(shí)需要調(diào)整的參數(shù)較多，根據(jù)相關(guān)研究(Li et al., 2019)和實(shí)踐經(jīng)驗(yàn)，本文調(diào)整了如表2中所示的參數(shù)以使模型發(fā)揮最佳性能，其具體的調(diào)整流程如圖3所示，即針對(duì)不同參數(shù)，依次找到最優(yōu)值；同時(shí)，在參數(shù)調(diào)整過程中使用了交叉驗(yàn)證(GridSearchCV)方法，其主要參數(shù)設(shè)置為：模型評(píng)估方法scoring=‘neg_mean_squared_error’，分割次數(shù)cv=5.

圖3 XGBoost的參數(shù)調(diào)整流程圖Fig.3 Parameter adjustment flow chart of XGBoost

表2 XGBoost的參數(shù)及其調(diào)整范圍Table 2 Parameters of XGBoost and its adjustment range

(3)電離層電子密度最終輸出

利用訓(xùn)練好的XGBoost預(yù)測下一時(shí)段中所有像素的電子密度IED0，并以此作為迭代初值，再次結(jié)合IART算法進(jìn)行層析反演，從而得到最終的電子密度IED2，同時(shí)，IED2又參與下一個(gè)滑動(dòng)窗口內(nèi)的訓(xùn)練樣本構(gòu)建.

總的來說，如表3所示，相較于傳統(tǒng)層析方法，盡管XGB-CIT的建模過程相對(duì)復(fù)雜，但其優(yōu)點(diǎn)也比較明顯.

表3 XGB-CIT與IART的構(gòu)建方法比較Table 3 Comparison of construction methods between XGB-CIT and IART

2 數(shù)值實(shí)驗(yàn)

利用湖南省CORS網(wǎng)的GNSS觀測數(shù)據(jù)對(duì)XGB-CIT的性能進(jìn)行驗(yàn)證，如圖4所示，重構(gòu)區(qū)域的經(jīng)緯度范圍分別為108.5°E—114.5°E和24.5°N—30.5°N，共包含73個(gè)GNSS測站和2個(gè)電離層測高儀站.重構(gòu)區(qū)域的高度范圍設(shè)為70～1000 km，并采用0.5°×0.5°×15 km的像素尺度對(duì)電離層區(qū)域進(jìn)行剖分，最終共有8928個(gè)像素.實(shí)驗(yàn)中使用了2015-12-30—2015-12-31連續(xù)兩天(每小時(shí)中只采用半小時(shí)數(shù)據(jù)，且每天23∶00—23∶59都沒有數(shù)據(jù))共46個(gè)時(shí)段(UT1, …, UT46)的數(shù)據(jù).同時(shí)，圖5給出了2015年12月30日和12月31日的ap、Dst和AE等地磁指數(shù)的時(shí)間序列，由圖可知，12月30日電離層處于平靜狀態(tài)，而12月31日則相反.

圖4 重構(gòu)區(qū)域及測站分布圖Fig.4 The reconstructed area and observation stations distribution

圖5 連續(xù)兩天的ap (a)、Dst (b)和AE (c)指數(shù)序列圖Fig.5 ap (a), Dst (b) and AE (c) index sequence diagram of two consecutive days

在討論XGB-CIT的性能之前，利用Jin等(2012)中的載波相位平滑偽距方法得到GNSS信號(hào)路徑上的TEC，反演區(qū)域的每個(gè)像素由IRI2016模型初始化.基于XGB-CIT的流程，利用每個(gè)反演周期的TEC數(shù)據(jù)，首先采用IART算法對(duì)重建區(qū)域的IED進(jìn)行反演，然后利用有觀測信息像素的反演IED以及像素特征，得到XGBoost的學(xué)習(xí)樣本，隨后對(duì)2015年12月31日23個(gè)時(shí)段內(nèi)每個(gè)像素的IED進(jìn)行滾動(dòng)預(yù)測，兩種樣本的結(jié)構(gòu)如表4所示.

表4 學(xué)習(xí)和預(yù)測樣本的時(shí)間跨度Table 4 Time span of learning and prediction samples

針對(duì)12月31日的每個(gè)反演時(shí)段，利用XGBoost與IRI2016預(yù)測的IED重新計(jì)算TEC并進(jìn)行比較，表5分別列出了23個(gè)反演時(shí)段內(nèi)TEC的均方根誤差(RMSE)和平均絕對(duì)誤差(ΔE)，其計(jì)算公式如下：

(10)

由表5可知，XGBoost的兩項(xiàng)誤差在所有的時(shí)段中都遠(yuǎn)低于IRI2016模型的，其中IRI2016與XGBoost的平均RMSE分別為6.97 TECU和2.23 TECU、平均ΔE分別為6.58 TECU和1.82 TECU，由此可知，基于RMSE統(tǒng)計(jì)，XGBoost的IED預(yù)測精度相較于IRI2016平均提高了68%，因此，初步驗(yàn)證了XGBoost算法結(jié)果用于XGB-CIT模型迭代初值的可行性.

由圖2可知，在XGB-CIT模型中，把XGBoost算法得到的IED作為迭代初值，并再次與IART算法相結(jié)合，得到最終的IED.結(jié)合(10)式，圖6分別

圖6 GNSS觀測值和層析模型獲取的TEC在23個(gè)層析時(shí)段的誤差比較Fig.6 Errors comparison between GNSS observations and TEC obtained by tomography model in 23 tomography periods

給出了傳統(tǒng)層析算法IART與XGB-CIT模型在23個(gè)反演時(shí)段內(nèi)TEC的反演誤差柱狀圖.由圖可知，除了UT26、UT37、UT38時(shí)段，XGB-CIT的RMSE和ΔE均小于或等于IART，其中IART與XGB-CIT的平均RMSE分別為0.55 TECU和0.53 TECU、平均ΔE分別為0.33 TECU和0.32 TECU，盡管基于平均RMSE和ΔE，XGB-CIT相較于IART的提高幅度不大，但由圖7可知，由于迭代初值的改

圖7 兩種層析方法的迭代次數(shù)比較Fig.7 Iteration times comparison of two tomography methods

進(jìn)，XGB-CIT的收斂速度比IART提高了20%，因此，這初步證明了XGB-CIT模型的優(yōu)良性能，同時(shí)進(jìn)一步驗(yàn)證了XGBoost算法的可行性.

圖8中給出了6個(gè)歷元中不同反演方法的電子密度廓線比較情況，其子圖從上往下分別代表UT26、UT27、UT33、UT34、UT41、UT42，每個(gè)子圖的左側(cè)是邵陽測高儀站(111.3°E, 27.1°N)，右側(cè)是武漢測高儀站(114.4°E, 30.5°N)，總的來說，相對(duì)于IRI2016和IART，XGB-CIT的電子密度廓線與測高儀的真實(shí)數(shù)據(jù)更加符合.由圖4可知，由于邵陽站位于研究區(qū)域的中心位置，其垂直方向的像素都有信號(hào)射線穿過，因此，除了子圖(f)，IART和XGB-CIT的電子密度廓線在邵陽站幾乎是重合的，這也說明，當(dāng)像素有信號(hào)射線穿過時(shí)，即使提高迭代初值精度，但其最終的反演精度也不會(huì)大幅提高，這也是對(duì)圖6中反演精度的一個(gè)直觀反映；然而武漢站位于研究區(qū)域的邊緣位置，其垂直方向的像素幾乎都沒有信號(hào)射線穿過，因此，當(dāng)?shù)踔稻忍岣吆?，似乎?huì)直接導(dǎo)致反演電子密度的精度提高幅度較大，但實(shí)際上并非如此.圖9中顯示了不同歷元中無信號(hào)射線像素的反演電子密度比較，每個(gè)子圖中從左往右，分別是IART與IRI2016的比較、XGBoost獲得的層析初值與IRI2016的比較、以及XGB-CIT與IRI2016的比較，其中子圖(a)、(b)、(c)和(d)分別代表歷元UT24、UT29、UT36和UT39，且他們的無信號(hào)射線穿過的像素比例分別為78.47%、63.89%、37.67%和32.46%.由圖可知，無信號(hào)射線像素的比例越大，其反演結(jié)果與初值越接近，而UT36和UT39中，盡管反演結(jié)果在初值基礎(chǔ)上有所改變，但其幅度也不大，因此，我們認(rèn)為對(duì)于無觀測區(qū)域的電子密度反演，其本質(zhì)是電子密度在時(shí)間和空間上的預(yù)測問題.總的來說，在有觀測信息像素的反演精度提高空間有限的情況下，能大幅提高無觀測信息像素的初值精度，這也是我們利用機(jī)器學(xué)習(xí)對(duì)迭代初值進(jìn)行精化的一個(gè)主要目的.

圖8 不同層析歷元的電子密度廓線對(duì)比Fig.8 Comparison of electron density profiles in different tomography epochs

圖9 無信號(hào)射線像素的反演電子密度比較Fig.9 Comparison of inversion ionospheric electron density of pixels without signal ray

圖10給出了IRI2016(左)、IART(中)、XGB-CIT(右)三種反演方法在UT34(上)和UT42(下)的電子密度三維示意圖，其中UT34和UT42的地方時(shí)分別是18∶00—18∶30 LT和02∶00—02∶30 LT，由圖5可知，UT34處于磁暴發(fā)生的上升階段，而UT42處于磁暴的短暫下降階段，因此導(dǎo)致兩個(gè)歷元中電子密度的差別較大.在兩個(gè)子圖中，由IART獲取的電子密度剖面圖在邊緣位置幾乎都與它的迭代初值IRI2016相同，而XGB-CIT的電子密度剖面圖中則不存在這個(gè)問題，結(jié)合圖8的分析可知，邊緣地帶的像素中由于缺少信號(hào)射線，是導(dǎo)致這種現(xiàn)象的主要原因.因此，這充分證明了利用機(jī)器學(xué)習(xí)進(jìn)行迭代初值精化的必要性和可行性.另外，圖11中顯示了四個(gè)反演歷元中由XGB-CIT獲取的VTEC與CODE的電離層格網(wǎng)產(chǎn)品(http:∥ftp.aiub.unibe.ch/CODE/IONO/2015/)對(duì)比情況，其中子圖(a)、(b)、(c)和(d)分別代表歷元UT24、UT31、UT34和UT37，但由于CODE提供的格網(wǎng)產(chǎn)品分辨率是2.5°×5°，而我們的層析像素分辨率是0.5°×0.5°，因此能覆蓋我們研究區(qū)域(24.5°N—30.5°N, 108.5°E—114.5°E)的CODE電離層格網(wǎng)點(diǎn)只有15個(gè)，特別是在有些歷元中，CODE提供的這15個(gè)格網(wǎng)值都是相同的，這導(dǎo)致CODE-VTEC圖并不能完全反映出我們研究區(qū)域內(nèi)電離層的真實(shí)變化，毫無疑問，我們利用XGB-CIT反演方法得到的電離層精度肯定高于CODE產(chǎn)品，這也充分說明在局部區(qū)域建立電離層模型的必要性和利用層析模型獲取高精度電離層信息的可行性.

圖10 IRI2016(a)、IART(b)和XGB-CIT(c)三種反演方法在兩個(gè)歷元UT34(A)和UT42(B)中的電子密度三維示意圖Fig.10 Three-dimensional schematic diagram of the electron density in two epochs UT34 (A) and UT42 (B) by three inversion methods, IRI2016 (a), IART (b) and XGB-CIT (c)

圖11 XGB-CIT獲取的VTEC (右)與CODE (左)的電離層格網(wǎng)產(chǎn)品比較Fig.11 Comparisons of the VTEC from XGB-CIT (right) with the CODE ionospheric grid product (left)

圖12、13分別給出了2015年12月31日XGB-CIT方法在120°E經(jīng)度面和25.5°N緯度面上的電子密度剖面圖，總的來說，兩個(gè)圖都很好地展現(xiàn)了電子密度的日變化，且電子密度在經(jīng)度面和緯度面上的變化比較一致，都在UT33(17∶00—17∶30 LT)達(dá)到峰值，一方面UT33處于磁暴的上升階段，另一方面本地時(shí)17∶00—17∶30 LT也是電離層比較活躍的階段，因此圖12—13中電子密度峰值出現(xiàn)的時(shí)間與實(shí)際情況也是比較吻合的.由圖12可知，在該磁暴期間，其南部的電子密度明顯大于北部的，但在緯度面上，電子密度有比較明顯的從東向西的轉(zhuǎn)移現(xiàn)象，另外，兩個(gè)圖中的電子密度從UT41(01∶00—01∶30 LT)開始又出現(xiàn)了短暫的回升，我們認(rèn)為這是受磁暴影響導(dǎo)致的異常情況.綜上所述，本文提出的新方法在磁暴期間能很好地重構(gòu)電離層電子密度三維結(jié)構(gòu)，且其三維結(jié)構(gòu)能精確地反映出電子密度的日變化及其異常情況，這有助于后續(xù)的電離層探測以及電離層異常擾動(dòng)研究.

圖12 XGB-CIT在120°E經(jīng)度面上的電子密度剖面圖Fig.12 Electron density profile of XGB-CIT on 120°E longitude plane

圖13 XGB-CIT在25.5°N緯度面上的電子密度剖面圖Fig.13 Electron density profile of XGB-CIT on 25.5°N latitude plane

3 結(jié)論

針對(duì)因觀測數(shù)據(jù)不足或分布不均導(dǎo)致的反演不適定性問題，本文提出了一種基于機(jī)器學(xué)習(xí)集成算法(XGBoost算法)的電離層層析新方法(XGB-CIT)，其中利用XGBoost算法獲取的迭代初值，相較于常用的IRI2016模型值，其精度提高了68%，而基于該初值構(gòu)建的XGB-CIT模型，在最終的反演精度上也展現(xiàn)了優(yōu)良的性能，尤其對(duì)無觀測信息像素的電子密度反演，其精度有了大幅提高，同時(shí)，根據(jù)其電子密度的二維和三維剖面圖分析可知，該方法也能精準(zhǔn)地反映出電子密度的日變化，另外，由于該方法可以對(duì)電離層未來的變化提供相對(duì)精準(zhǔn)的信息，因此其預(yù)測的電子密度還可以為電離層異常擾動(dòng)的探測提供準(zhǔn)確的背景值，因此，我們認(rèn)為XGB-CIT是一種新的、可靠的電離層探測手段.

致謝感謝湖南省測繪科技研究所提供的GNSS觀測數(shù)據(jù)，感謝地球物理數(shù)據(jù)中心——地球系統(tǒng)科學(xué)數(shù)據(jù)共享平臺(tái)提供電離層測高儀數(shù)據(jù).