一種可選的電離層TEC 區(qū)域重構(gòu)方法

車磊 歐明 陳奇東 蔡紅濤 甄衛(wèi)民 陳龍江 靳睿敏

（1.中國電波傳播研究所,青島 266107；2.西安電子科技大學，西安 710071；3.武漢大學電子信息學院，武漢 430072）

引言

作為表征電離層變化的一個重要特征參量，電離層總電子含量(total electron content,TEC)反映了眾多電離層不同空間的變化特性[1-2].通過空間插值(spatial interpolation)的方法得到整個區(qū)域TEC 的預測值是TEC 區(qū)域重構(gòu)常用的技術(shù)方法[3].由于電離層暴、太陽活動等自然現(xiàn)象的頻發(fā)，電離層TEC 的非線性、非平穩(wěn)變化可能存在波動和異常，影響無線電波的傳播特性，不可避免地引起包括衛(wèi)星導航、通訊、雷達和定位異常等問題[4]，因此，如何清晰地描述電離層TEC 的空間分布特征，準確地實現(xiàn)特定區(qū)域的電離層TEC 區(qū)域重構(gòu)，完成電離層監(jiān)測與預警，一直是相關(guān)領(lǐng)域研究的熱點和難點問題[4-5].

當前，對于區(qū)域電離層TEC 區(qū)域重構(gòu)主要通過空間插值的方法來估算一定區(qū)域范圍內(nèi)電離層TEC.這里提及的空間插值方法主要是基于已知觀測站點獲取的電離層TEC 真實值，通過插值的原理來估計其他位置站點的電離層TEC 預測值，其根本原理是通過基于構(gòu)建的函數(shù)關(guān)系理論模型，綜合已知監(jiān)測站點的電離層TEC 空間位置關(guān)系以及空間相關(guān)性，從而估算其他任意點的電離層TEC[6].空間插值方法本質(zhì)上追求構(gòu)建盡可能符合原始觀測數(shù)據(jù)的函數(shù)關(guān)系理論模型[7].

空間插值方法種類眾多，應用也十分廣泛.常用的插值方法包括反距離加權(quán)法、線性內(nèi)插法、泰森多邊形法、樣條函數(shù)法、移動擬合法、趨勢分析法、克里金插值法等.克里金空間插值方法(Kriging interpolation)適用于樣本數(shù)據(jù)存在隨機性和結(jié)構(gòu)性特征的場景，應用于空氣污染、降雨、環(huán)境監(jiān)測等領(lǐng)域，并結(jié)合電離層TEC 空間分布的特殊性及相關(guān)性，得到了廣泛應用[8].克里金空間插值方法又稱空間自協(xié)方差最佳插值法[9]，基于區(qū)域化變量(regionalized variable)自身具有的隨機性特征和結(jié)構(gòu)性特征為基礎(chǔ)，同時通過變異函數(shù)(variogram)對區(qū)域化變量進行空間描述，模擬地理現(xiàn)象空間分布的相關(guān)性和變異性，因此能夠挖掘區(qū)域化變量的空間結(jié)構(gòu)和空間變化規(guī)律[10].電離層TEC 正是具有這種隨機性(不確定性)與結(jié)構(gòu)性(相關(guān)性)雙重特征的區(qū)域化變量，應用克里金空間插值實現(xiàn)對電離層TEC 區(qū)域重構(gòu)，其實質(zhì)在于通過已經(jīng)位置點的電離層TEC 內(nèi)插或外推的方式，對待估位置點電離層TEC 的取值進行無偏、最優(yōu)估計[11].

諸多專家學者通過克里金空間插值方法實現(xiàn)了電離層TEC 區(qū)域重構(gòu).Stanislawska 等人改進了克里金空間插值方法，通過加入電離層空間距離的影響因素，實現(xiàn)對歐洲區(qū)域的電離層TEC 區(qū)域重構(gòu)[12-13]；陳春等人根據(jù)foF2時間和空間相關(guān)性，通過克里金空間插值方法引入電離層空間距離、經(jīng)度因子和緯度因子等參數(shù)實現(xiàn)了電離層TEC 區(qū)域重構(gòu)[3]；劉瑞源等人提出了一種適用于中國地區(qū)電離層TEC 的短期預報方法，并定量分析了低緯站和邊緣站對區(qū)域重構(gòu)的精度誤差影響[14]；然而通過克里金空間插值方法實現(xiàn)電離層TEC 區(qū)域重構(gòu)過程中，區(qū)域重構(gòu)的精度取決于模型對待估算點位和已知樣本點位空間位置及其空間相關(guān)性這兩者的反映程度[7].但隨之帶來的問題是，當克里金空間插值方法擬合變異函數(shù)時，傳統(tǒng)理論變異函數(shù)模型面臨函數(shù)曲線固定、空間細節(jié)變化無法反映以及模型選取人為主觀等問題[15].

為解決上述問題，本文提出一種可選的TEC 區(qū)域重構(gòu)方法，從電離層TEC 實際變化趨勢出發(fā)，采用最小二乘支持向量機(least squares support vector machine,LS-SVM)擬合實驗變異函數(shù)，實現(xiàn)電離層TEC 區(qū)域重構(gòu).為驗證此方法的準確性，本文選用中國陸態(tài)網(wǎng)地基GNSS 臺站某時刻三組不同時刻穿刺點垂直總電子含量(vertical TEC,VTEC)值作為實測數(shù)據(jù)，同時選用普通克里金空間插值方法中的指數(shù)理論變異函數(shù)模型、球狀理論變異函數(shù)模型以及本文模型進行實驗.結(jié)果表明，本文提供的電離層TEC區(qū)域重構(gòu)方法計算的均方根誤差(root mean square error,RMSE)和平均絕對誤差(mean absolute error,MAE)均小于其他兩種理論變異函數(shù)模型，插值精度最好，為電離層TEC 區(qū)域重構(gòu)提供了一種可選的思路.

1 相關(guān)理論原理

1.1 克里金空間插值方法

克里金空間插值方法在有限的區(qū)域范圍內(nèi)對區(qū)域化變量進行無偏最優(yōu)估計.區(qū)域化變量以自身具有的隨機性特征和結(jié)構(gòu)性特征為基礎(chǔ)，對相關(guān)性和連續(xù)性等要素特點進行空間描述，模擬地理現(xiàn)象空間分布的相關(guān)性和變異性，因此能夠挖掘區(qū)域化變量的空間結(jié)構(gòu)和空間變化規(guī)律.同時借助變異函數(shù)，既能夠描述其隨機性變化過程，又能夠反映區(qū)域化變量空間結(jié)構(gòu)性變化過程[15].

變異函數(shù) γ(x,h)定義為區(qū)域化變量z(x) 在x軸 方向上，z(x)在點位x和x+h處變量值之差的方差一半，如式(1)計算：

式中：Var[·]表 示方差；E[·]表示期望.

在二階平穩(wěn)假設(shè)的情況下，對任意樣本點對的距離h有

式(3)可理解為變異函數(shù) γ(x,h)依賴于方向和距離兩個變量變化，倘若變異函數(shù)僅僅依賴于距離變化時，則 γ(x,h)可以寫為 γ(h)，同時稱 γ(h)為 各向同性.此時，離散樣本數(shù)據(jù)的實驗變異函數(shù)可以通過式(4)所示：

式中：h代表樣本點對的空間距離；N(h)代表當樣本點對距離為h時，所有樣本點對的總數(shù)量；z(xi)和z(xi+h)分別表示z(x)在點位xi和點位xi+h處的實際觀測值，即真實值.

克里金空間插值方法提供了包括線性模型、高斯模型、球狀模型、指數(shù)模型等在內(nèi)的幾種常用理論變異函數(shù)模型[16].接下來以普通克里金空間插值方法為例，描述其實現(xiàn)的基本原理.

該方法實現(xiàn)空間估計主要滿足兩個基本條件：無偏性和估計方差最小，即：

求解得到如下方程組：

式中：λi為 空間權(quán)重系數(shù)，表示點位xi處的區(qū)域化變量z(xi)值對待估點位x0的貢獻程度；j=1,...,k，k為監(jiān)測站點的總個數(shù)；為點位xi與點位xj距離下的實驗變異函數(shù)值；μ為拉格朗日乘子；為待估點位x0與點位xj距離下的實驗變異函數(shù)值.

將式(7)用矩陣展開，得到：

普通克里金空間插值方程為

解得

最終得到待估位置點x0處的區(qū)域化變量估計值z*(x0)：

式中：x1，...，xk為已知樣本點位；z(x1)，...，z(xk)為對應樣本點位的實際觀測值.

綜上所述，普通克里金空間插值方法基于無偏性和估計方差最小兩個原則，在此基礎(chǔ)上建立了含有約束條件的拉格朗日函數(shù)，插值結(jié)果的好壞完全取決于權(quán)重系數(shù)[8].通過約束條件和求極值問題解決待估位置點位區(qū)域化變量的無偏、最優(yōu)估計[17].

1.2 LS-SVM

SVM 是由Vapnik 等人在20 世紀60 年代提出的一種有限樣本機器學習理論，以統(tǒng)計學習理論為基礎(chǔ)，SVM 模型不僅追求模型本身的泛化性能，且追求有限個樣本條件下的最優(yōu)解[18].SVM 模型構(gòu)建過程中，以結(jié)構(gòu)風險最小化為基礎(chǔ)，主要用于處理小樣本、高維數(shù)、非線性、局部最優(yōu)解等問題.LS-SVM作為SVM 演變的一種類型，其原理是計算損失函數(shù)時，以平方和誤差損失函數(shù)取代Vapnik 的ε 不敏感損失函數(shù)，同時構(gòu)建等式約束條件[19-20]，計算原理如下：

給定N個樣本數(shù)據(jù)集 {其中第m個 輸入xm對應的輸出值為ym，回 歸函數(shù)f(x)的基本形式如式(14)所示：

式中：ω為權(quán)系數(shù)向量(列向量)；φ(x)為輸入空間到特征空間的映射函數(shù)，即低維空間向高維空間轉(zhuǎn)化的映射函數(shù)；b為常數(shù)項.

模型優(yōu)化函數(shù)如式(15)所示：

對應得到

同時，上述公式也需滿足以下等式約束條件：

綜合上述條件，構(gòu)建含有約束條件的拉格朗日函數(shù)，基于KKT(Karush-Kuhn-Tucker)條件，求解方程組，過程不再贅述.

最終回歸函數(shù)模型f(x)如下：

式中：αm表示拉格朗日乘子；K(xm,x)表示核函數(shù).

2 計算流程

通過本文方法實現(xiàn)電離層TEC 區(qū)域重構(gòu)方法，步驟如下：

步驟一：通過已知監(jiān)測站點的地理坐標構(gòu)建該區(qū)域網(wǎng)格化坐標點及范圍.根據(jù)已知站點的地理坐標經(jīng)度和緯度，以經(jīng)度的最小值和最大值構(gòu)建網(wǎng)格區(qū)域的長邊，以緯度的最小值和最大值構(gòu)建網(wǎng)格區(qū)域的寬邊，完成區(qū)域網(wǎng)格化坐標點及范圍的創(chuàng)建.

步驟二：借助離散變異函數(shù)式(4)，得到所有站點電離層TEC 樣本點對的實驗變異函數(shù)值，若實驗變異函數(shù)值較多則進行分組操作，便于后續(xù)擬合實驗變異函數(shù).

步驟三：采用LS-SVM 擬合實驗變異函數(shù)值，得到理論變異函數(shù)模型.通過LS-SVM 與克里金空間插值相結(jié)合，采用LS-SVM 擬合實驗變異函數(shù)，能夠更精確地刻畫變異函數(shù)，反映電離層TEC 空間變化趨勢.

步驟四：根據(jù)式(7)建立方程組求解權(quán)重系數(shù) λi.構(gòu)建含有約束條件的拉格朗日函數(shù)，通過步驟三得到理論變異函數(shù)模型及其他參數(shù)，求解權(quán)重系數(shù) λi.

步驟五：根據(jù)式(13)計算待估位置點的電離層總電子含量z*(x0)，即實現(xiàn)了網(wǎng)格內(nèi)所有待估位置點電離層總電子含量無偏、最優(yōu)的估計.

3 實驗分析

3.1 精度誤差指標

為檢驗本文提出的電離層TEC 區(qū)域重構(gòu)方法，借助RMSE 和MAE 兩個精度誤差指標[21].RMSE 反映的是靈敏度變化和可能存在的極值誤差效應，MAE 反映的是總體精度誤差.兩個精度誤差指標數(shù)值越小，表明電離層TEC 區(qū)域重構(gòu)效果越好.RMSE和MAE 的定義如下：

式中：S代表待估位置點的總個數(shù)；待估位置點的電離層TE C 真實值為yi；待估位置點的電離層TEC 估計值為

3.2 實驗數(shù)據(jù)

以中國陸態(tài)網(wǎng)地基GNSS 臺站三組不同時刻穿刺點觀測值作為實驗數(shù)據(jù).采用文獻[22]中方法對數(shù)據(jù)進行預處理得到斜向總電子含量(slant TEC,STEC)，再通過薄層模型(single layer model,SLM)映射函數(shù)[23]轉(zhuǎn)化因子計算得到VTEC值，三組數(shù)據(jù)VTEC 分布如圖1 所示.其中，(a)表示2017-09-05T14:00:00UT(第一組)穿刺點VTEC 空間分布，(b)表示2017-09-07T03:31:00UT(第二組)穿刺點VTEC 空間分布，(c)表示2017-09-10T06:18:30UT(第三組)穿刺點VTEC 空間分布.

圖1 不同時刻穿刺點VTEC 空間分布Fig.1 The schematic of puncture points VTEC spatial distribution at different times

同時為對比各組數(shù)據(jù)理論函數(shù)變異模型計算得到的MAE 和RMSE 兩類評價指標誤差，將不同時刻穿刺點VTEC 值劃分為插值數(shù)據(jù)和測試數(shù)據(jù).隨機均勻抽取樣本，通常插值數(shù)據(jù)占總樣本數(shù)據(jù)的75%，測試數(shù)據(jù)占總樣本數(shù)據(jù)25%，同時兩者數(shù)據(jù)不重復，以滿足計算過程的獨立性[18].各組數(shù)據(jù)分類個數(shù)如表1 所示.

表1 三組數(shù)據(jù)樣本類別及數(shù)量信息Tab.1 Three groups of data sample category and quantity information

區(qū)域重構(gòu)范圍為70～140°E、15～55°N，基本覆蓋中國大陸區(qū)域，重構(gòu)數(shù)據(jù)的網(wǎng)格空間分辨率為0.5°×0.5°.

3.3 實驗結(jié)果

普通克里金空間插值方法中理論變異函數(shù)模型選用指數(shù)理論變異函數(shù)模型和球狀理論變異函數(shù)模型，用于對比本文方法實現(xiàn)的電離層TEC 區(qū)域重構(gòu)效果.圖2 給出實驗過程中實驗變異函數(shù)的擬合曲線(以第二組數(shù)據(jù)為例).明顯可以看出，LS-SVM 理論變異函數(shù)較好地擬合了所有實驗變異函數(shù)值，整體擬合曲線符合電離層TEC 數(shù)據(jù)本身的空間變化特征，展現(xiàn)了數(shù)據(jù)本身的空間變化趨勢.

圖2 變異函數(shù)擬合曲線Fig.2 Variogram fitting curve

圖3、4 和5 分別為三組數(shù)據(jù)電離層TEC 區(qū)域重構(gòu)效果.

圖3 2017-09-05T14:00:00UT 電離層TEC 區(qū)域重構(gòu)效果Fig.3 The effect diagram of ionospheric TEC region reconstruction at 14:00:00UT on September 5,2017

圖4 2017-09-07T03:31:00UT 電離層TEC 區(qū)域重構(gòu)效果Fig.4 The effect diagram of ionospheric TEC region reconstruction at 03:31:00UT on September 7,2017

圖5 2017-09-10T06:18:30UT 電離層TEC 區(qū)域重構(gòu)效果Fig.5 The effect diagram of ionospheric TEC region reconstruction at 06:18:30UT on September 10,2017

從圖3～5 可以看出，對比三種方法重構(gòu)出來的電離層TEC 分布結(jié)果，LS-SVM 理論變異函數(shù)模型實現(xiàn)的區(qū)域重構(gòu)效果大致符合通過克里金空間插值方法球狀理論變異函數(shù)模型和指數(shù)理論變異函數(shù)模型重構(gòu)的效果，三者空間變化趨勢相近，空間變化呈現(xiàn)整體變化平緩、低緯地區(qū)偏高、高緯地區(qū)偏低的趨勢，符合電離層TEC 空間分布特征.因此保守地認為，本文方法的確可以作為一種可選的電離層TEC區(qū)域重構(gòu)方法.

從TEC 區(qū)域重構(gòu)結(jié)果分析，TEC 值變化呈現(xiàn)從低緯地區(qū)向中高緯地區(qū)逐步銳減的趨勢，本文選用的三組數(shù)據(jù)，包括夜間VTEC 和日間VTEC.第一組數(shù)據(jù)選用夜間VTEC，取值范圍[0,35] TECU，變化平緩；第二組和第三組數(shù)據(jù)選用日間VTEC，取值范圍[0,65] TECU，變化劇烈.三組實驗數(shù)據(jù)整體時間變化趨勢基本符合電離層TEC 的空間變化趨勢，與文獻[24]通過改進克里金空間插值方法給出的該區(qū)域電離層TEC 重構(gòu)結(jié)果基本一致.

同時，為定量對比三種方法得到的區(qū)域重構(gòu)結(jié)果精度誤差，分別統(tǒng)計了不同理論變異函數(shù)模型下的RMSE 和MAE，如表2、表3 和表4 所示.

表2 第一組數(shù)據(jù)精度誤差對比Tab.2 Accuracy error comparison (the first group)

表3 第二組數(shù)據(jù)精度誤差對比Tab.3 Accuracy error comparison (the second group)

表4 第三組數(shù)據(jù)精度誤差對比Tab.4 Accuracy error comparison (the third group)

表2～4 直觀地顯示出不同理論變異函數(shù)模型計算得到的電離層TEC 區(qū)域重構(gòu)結(jié)果，傳統(tǒng)克里金插值擬合模型中，指數(shù)理論變異函數(shù)模型誤差小于球狀理論變異函數(shù)模型；LS-SVM 理論變異函數(shù)模型RMSE 誤差分別為1.54 TECU、1.76 TECU 和2.45 TECU，MAE 誤差分別為1.21 TECU、1.23 TECU 和1.62 TECU.因此，這三種區(qū)域重構(gòu)方法中，本文方法效果最好.同時，對比分析表2～4，三組實驗數(shù)據(jù)精度誤差存在差異性，作者分析認為，該差異性主要由數(shù)據(jù)本身的空間屬性變化差異性引起.白天VETC值差異性明顯，同時該研究區(qū)域中低緯地區(qū)處于電離層異常區(qū)域，梯度變化十分劇烈，引起的極值變化較多，因此精度誤差偏大.相反地，夜間VETC 值變化則較為平緩且極值點較少，誤差精度反而較低.

4 結(jié)論

電離層TEC 區(qū)域重構(gòu)常采用克里金空間插值方法，但克里金空間插值方法擬合變異函數(shù)時，傳統(tǒng)理論變異函數(shù)模型面臨函數(shù)曲線固定、空間細節(jié)變化無法反映以及模型選取人為主觀等問題，為解決諸如此類的問題，本文提出了一種可選的電離層TEC區(qū)域重構(gòu)方法.這種方法通過結(jié)合電離層TEC 數(shù)據(jù)的空間變化特征，能夠更精確地刻畫變異函數(shù)，反映電離層TEC 數(shù)據(jù)在空間中的實際變化趨勢，從而在一定程度上提高了電離層TEC 區(qū)域重構(gòu)的精度.

需要與讀者說明的是，本文提出的這種電離層TEC 區(qū)域重構(gòu)方法主要基于兩方面考慮：1)電離層TEC 數(shù)據(jù)本身是一種具有地理空間特性的區(qū)域化變量，在變化過程中與周圍鄰域的位置點及位置數(shù)據(jù)產(chǎn)生作用，體現(xiàn)為空間屬性的相關(guān)性或相似性特征，因此適合采用克里金空間插值方法實現(xiàn)區(qū)域重構(gòu)；2)電離層TEC 區(qū)域重構(gòu)過程中，理論變異函數(shù)模型擬合實驗變異函數(shù)離散點，如何實現(xiàn)更精確的擬合和更優(yōu)良的泛化性能，是追求無偏性和最優(yōu)解的目標.通過真實數(shù)據(jù)進行仿真驗證，結(jié)果表明該方法在一定程度上可以提高插值精度，即提高電離層TEC 區(qū)域重構(gòu)的準確性，為研究此領(lǐng)域的作者提供一種思路供大家參考借鑒.

但同時值得思考的是，實驗過程中也存在一些需要后續(xù)繼續(xù)研究的內(nèi)容：不同樣本數(shù)據(jù)即電離層TEC 數(shù)據(jù)計算得到的區(qū)域重構(gòu)精度不盡相同，作者分析認為一方面是數(shù)據(jù)本身可能存在特殊性，另一方面LS-SVM 模型本身的泛化性能會根據(jù)數(shù)據(jù)的不同而存在差異性.后續(xù)研究不僅需要更多的數(shù)據(jù)參與實驗驗證，同時引入其他人工智能方法實現(xiàn)不同方法之間的對比提高.

致謝：本文GNSS 觀測數(shù)據(jù)從中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(Crustal Movement Observation Network of China，CMNOC)獲取，作者在此表示感謝.本工作得到國家重點研發(fā)計劃(2018YFF01013702 和2018 YFB0505100)的資助.