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

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

引言

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

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

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

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

為解決上述問題，本文提出一種可選的TEC 區(qū)域重構(gòu)方法，從電離層TEC 實(shí)際變化趨勢出發(fā)，采用最小二乘支持向量機(jī)(least squares support vector machine,LS-SVM)擬合實(shí)驗(yàn)變異函數(shù)，實(shí)現(xiàn)電離層TEC 區(qū)域重構(gòu).為驗(yàn)證此方法的準(zhǔn)確性，本文選用中國陸態(tài)網(wǎng)地基GNSS 臺站某時刻三組不同時刻穿刺點(diǎn)垂直總電子含量(vertical TEC,VTEC)值作為實(shí)測數(shù)據(jù)，同時選用普通克里金空間插值方法中的指數(shù)理論變異函數(shù)模型、球狀理論變異函數(shù)模型以及本文模型進(jìn)行實(shí)驗(yàn).結(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ū)域化變量進(jìn)行無偏最優(yōu)估計.區(qū)域化變量以自身具有的隨機(jī)性特征和結(jié)構(gòu)性特征為基礎(chǔ)，對相關(guān)性和連續(xù)性等要素特點(diǎn)進(jìn)行空間描述，模擬地理現(xiàn)象空間分布的相關(guān)性和變異性，因此能夠挖掘區(qū)域化變量的空間結(jié)構(gòu)和空間變化規(guī)律.同時借助變異函數(shù)，既能夠描述其隨機(jī)性變化過程，又能夠反映區(qū)域化變量空間結(jié)構(gòu)性變化過程[15].

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

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

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

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

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

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

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

求解得到如下方程組：

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

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

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

解得

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

式中：x1，...，xk為已知樣本點(diǎn)位；z(x1)，...，z(xk)為對應(yīng)樣本點(diǎn)位的實(shí)際觀測值.

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

1.2 LS-SVM

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

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

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

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

對應(yīng)得到

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

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

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

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

2 計算流程

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

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

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

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

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

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

3 實(shí)驗(yàn)分析

3.1 精度誤差指標(biāo)

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

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

3.2 實(shí)驗(yàn)數(shù)據(jù)

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

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

同時為對比各組數(shù)據(jù)理論函數(shù)變異模型計算得到的MAE 和RMSE 兩類評價指標(biāo)誤差，將不同時刻穿刺點(diǎn)VTEC 值劃分為插值數(shù)據(jù)和測試數(shù)據(jù).隨機(jī)均勻抽取樣本，通常插值數(shù)據(jù)占總樣本數(shù)據(jù)的75%，測試數(shù)據(jù)占總樣本數(shù)據(jù)25%，同時兩者數(shù)據(jù)不重復(fù)，以滿足計算過程的獨(dú)立性[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 實(shí)驗(yàn)結(jié)果

普通克里金空間插值方法中理論變異函數(shù)模型選用指數(shù)理論變異函數(shù)模型和球狀理論變異函數(shù)模型，用于對比本文方法實(shí)現(xiàn)的電離層TEC 區(qū)域重構(gòu)效果.圖2 給出實(shí)驗(yàn)過程中實(shí)驗(yàn)變異函數(shù)的擬合曲線(以第二組數(shù)據(jù)為例).明顯可以看出，LS-SVM 理論變異函數(shù)較好地擬合了所有實(shí)驗(yàn)變異函數(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ù)模型實(shí)現(xiàn)的區(qū)域重構(gòu)效果大致符合通過克里金空間插值方法球狀理論變異函數(shù)模型和指數(shù)理論變異函數(shù)模型重構(gòu)的效果，三者空間變化趨勢相近，空間變化呈現(xiàn)整體變化平緩、低緯地區(qū)偏高、高緯地區(qū)偏低的趨勢，符合電離層TEC 空間分布特征.因此保守地認(rèn)為，本文方法的確可以作為一種可選的電離層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í)驗(yàn)數(shù)據(jù)整體時間變化趨勢基本符合電離層TEC 的空間變化趨勢，與文獻(xiàn)[24]通過改進(jìn)克里金空間插值方法給出的該區(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í)驗(yàn)數(shù)據(jù)精度誤差存在差異性，作者分析認(rèn)為，該差異性主要由數(shù)據(jù)本身的空間屬性變化差異性引起.白天VETC值差異性明顯，同時該研究區(qū)域中低緯地區(qū)處于電離層異常區(qū)域，梯度變化十分劇烈，引起的極值變化較多，因此精度誤差偏大.相反地，夜間VETC 值變化則較為平緩且極值點(diǎn)較少，誤差精度反而較低.

4 結(jié)論

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

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

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

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