電離層VTEC異常分布產(chǎn)出系統(tǒng)設(shè)計與實現(xiàn)

劉 江，王 明，陳 聰，饒?zhí)珖?/p>

（1.四川省地震局，成都 610041；2.四川啟明豪盛軟件科技有限公司，成都 610031）

0 引言

近年來，地基電離層探測技術(shù)已被廣泛應用于地震前兆研究[1]。其中，電離層垂直電子濃度總含量（VTEC）是電離層特性的一個重要參量，對于電離層監(jiān)測具有十分重要的意義。而GNSS觀測由于具有觀測及時、覆蓋范圍廣、觀測精度高、成本低、效率高等優(yōu)勢而被廣泛應用于電離層研究，是監(jiān)測VTEC實時變化的最主要技術(shù)手段之一[2-4]。

隨著地震預報需求的日益迫切及電離層探測水平的不斷提高，地震-電離層預報模式逐漸成為地震預報新的研究方向。自2018年以來，四川省地震局GNSS中心與中國地震局地殼運動監(jiān)測工程研究中心合作，快速獲取陸態(tài)網(wǎng)絡(luò)中國區(qū)域GNSS連續(xù)運行基準站觀測數(shù)據(jù)（包含GPS、GLONASS衛(wèi)星信息），解算中國區(qū)域電離層VTEC[5-7]，開展臨震電離層異常擾動的監(jiān)測和研究。為有效監(jiān)測全球、中國區(qū)域電離層異常擾動變化，研究震前電離層異常短臨特性，本文設(shè)計開發(fā)電離層VTEC異常分布產(chǎn)出系統(tǒng)。該軟件可快速產(chǎn)出全球VTEC、中國區(qū)域VTEC異常分布，自動產(chǎn)出結(jié)果延遲1~2 d，同時對空間環(huán)境參數(shù)F10.7、Kp、Dst及中國區(qū)域格網(wǎng)點VTEC時序變化進行監(jiān)控，極大地提高了電離層VTEC異常監(jiān)測能力，進一步測試中國區(qū)域地震—電離層監(jiān)測預報的應用效能，嘗試為地震監(jiān)測預報提供電離層前兆異常判定依據(jù)。

1 系統(tǒng)設(shè)計

1.1 系統(tǒng)結(jié)構(gòu)

基于Windows 2008/2012操作系統(tǒng)，采用C/S架構(gòu)（Client-Server，即客戶端/服務(wù)器模式），C#作為后臺開發(fā)語言，MySQL作為后臺數(shù)據(jù)庫管理系統(tǒng)，ArcGIS作為圖片成果產(chǎn)出平臺，通過分層技術(shù)體系結(jié)構(gòu)實現(xiàn)系統(tǒng)運行的穩(wěn)定性和交互性。整體系統(tǒng)框架分為基礎(chǔ)層、數(shù)據(jù)層、應用層、業(yè)務(wù)層和用戶層五個層次（見圖1）。其中，基礎(chǔ)層包括操作系統(tǒng)、數(shù)據(jù)庫、網(wǎng)絡(luò)監(jiān)控軟件等基礎(chǔ)設(shè)施；數(shù)據(jù)層包括基礎(chǔ)數(shù)據(jù)和業(yè)務(wù)數(shù)據(jù)，基礎(chǔ)數(shù)據(jù)包含屬性數(shù)據(jù)、地圖數(shù)據(jù)，業(yè)務(wù)數(shù)據(jù)包含星歷及時鐘等控制文件、空間環(huán)境參數(shù)數(shù)據(jù)、VTEC格網(wǎng)數(shù)據(jù)；應用層為該系統(tǒng)提供業(yè)務(wù)支撐，包括ArcGIS平臺、數(shù)據(jù)庫訪問引擎、數(shù)據(jù)下載等；業(yè)務(wù)層按照系統(tǒng)業(yè)務(wù)規(guī)劃和需求，實現(xiàn)所有業(yè)務(wù)功能，包括觀測數(shù)據(jù)解算、VTEC異常提取、全球VTEC異常分布、中國區(qū)域VTEC異常分布、空間環(huán)境參數(shù)及格網(wǎng)點VTEC時序曲線、數(shù)據(jù)查詢等；用戶層面向系統(tǒng)用戶，提供成果查詢、成果展示和成果手動產(chǎn)出等功能。

圖1 系統(tǒng)框架結(jié)構(gòu)圖Fig.1 System frame structure diagram

1.2 系統(tǒng)處理流程

系統(tǒng)選取陸態(tài)網(wǎng)絡(luò)中國區(qū)域70個均勻分布的GNSS基準站觀測數(shù)據(jù)（站點分布見圖2），通過自動和手動模式采用FTP方式下載各類數(shù)據(jù)，選擇球諧函數(shù)作為電離層VTEC的擬合模型，自動調(diào)用Bernese數(shù)據(jù)處理軟件，建立中國區(qū)域電離層格網(wǎng)模型，解算中國區(qū)域電離層VTEC格網(wǎng)數(shù)據(jù)（經(jīng)緯度分辨率1°×1°）。結(jié)合CODE全球電離層格網(wǎng)數(shù)據(jù)，采用滑動四分位距法分別提取全球、中國區(qū)域VTEC異常信息，將全球VTEC、中國區(qū)域VTEC異常信息及空間環(huán)境參數(shù)導入數(shù)據(jù)庫，利用ArcGIS平臺分別產(chǎn)出全球VTEC異常分布、中國區(qū)域VTEC異常分布、空間環(huán)境參數(shù)及中國區(qū)域格網(wǎng)點VTEC時序曲線，具體處理流程如下（見圖3）。

圖2 陸態(tài)網(wǎng)絡(luò)選取站點分布圖Fig.2 Distribution map of selected sites in crustal movement observation network of China

圖3 系統(tǒng)處理流程圖Fig.3 System processing flow diagram

1.3 系統(tǒng)分析方法

1.3.1 中國區(qū)域電離層格網(wǎng)數(shù)據(jù)解算

對于中國區(qū)域電離層格網(wǎng)數(shù)據(jù)解算，可采用以下方法：利用GPS雙頻接收機2個頻率載波相位觀測值及其加載的偽距碼觀測值，通過載波相位平滑偽距觀測量，形成電離層殘差組合[7]，公式如下：

式（1）中：P4為載波相位平滑偽距，c為光速，f1、f2為衛(wèi)星載波頻率，Δbs是衛(wèi)星系統(tǒng)硬件延遲，Δbr是接收機硬件延遲，TEC為傾斜路徑總電子含量。

采用球諧函數(shù)對電離層VTEC進行建模，數(shù)學表達式[7]：

式（2）中：VTEC（β，s）為垂直總電子含量，nmax為球諧展開的最大階數(shù)，β為穿刺點地理緯度，s為日固系中穿刺點經(jīng)度，為n階m次歸一化締合勒讓德函數(shù)，anm和bnm為待估計電離層模型參數(shù)。TEC與VTEC數(shù)學表達式[7]：

式（4）中：F（Z）是電離層投影函數(shù)，單層模型高度H取350 km，α為0.978，R為6387 km，Z為衛(wèi)星天頂距。

將式（2）、（3）、（4）代入式（1），形成最終解算方程[7]：

式（5）中：anm，bnm，Δbs，Δbr為待估計未知量，采用每2 h多項式模型與單日固定值的接收機硬件延遲和衛(wèi)星硬件延遲共同建立法方程，利用最小二乘求解，解算出每2 h的模型參數(shù)和單日硬件延遲后，便可解算中國區(qū)域電離層格網(wǎng)數(shù)據(jù)[6]。選擇球諧函數(shù)作為電離層VTEC擬合模型，算法成熟，因此VTEC解算結(jié)果具有較高可信度。

1.3.2 異常檢驗方法

選用滑動四分位距法提取VTEC異常信息，選取前15 d相同時刻的觀測值求取相應四分位距（IQR）及中位數(shù)（M），建立上下邊界閾值M±1.5×IQR，當同一時刻觀測值超出上或下邊界時，該觀測值視為異常，且高于上邊界為正異常，低于下邊界為負異常。滑動四分位距法M±1.5×IQR的判斷閾值，約為標準差的2倍，檢測置信度為95%。

從二維空間分布分析全球及中國區(qū)域電離層VTEC異常分布，分析VTEC異常的ΔTEC空間分布，計算方法如下：

其中，ΔTEC為VTEC的擾動異常，L1和L2分別為VTEC背景上下邊界閾值，K為待分析時刻VTEC值。

1.4 系統(tǒng)功能模塊

根據(jù)業(yè)務(wù)需求，該系統(tǒng)分為系統(tǒng)配置、自動產(chǎn)出、手動產(chǎn)出、數(shù)據(jù)服務(wù)、任務(wù)查詢等5個主要部分，具體功能模塊如下（見圖4）。

圖4 系統(tǒng)功能圖Fig.4 System function diagram

2 系統(tǒng)實現(xiàn)

2.1 系統(tǒng)界面

系統(tǒng)基于Microsoft Visual Studio開發(fā)環(huán)境，C#作為后臺開發(fā)語言。系統(tǒng)登錄界面（見圖5），輸入用戶名、密碼登錄進入系統(tǒng)運行主界面。主界面分為功能菜單、用戶操作及圖形展示三部分，根據(jù)用戶需求分別展示全球VTEC異常分布（見圖6）、中國區(qū)域VTEC異常分布（見圖7）。

圖5 系統(tǒng)登錄界面Fig.5 System login interface

功能菜單位于主界面左上側(cè)，包含系統(tǒng)配置、手動產(chǎn)出、時序曲線、數(shù)據(jù)轉(zhuǎn)換、數(shù)據(jù)導出、數(shù)據(jù)下載、任務(wù)列表等功能。

用戶操作位于功能菜單下方，包含查詢條件和查詢結(jié)果兩部分。其中，查詢條件需要設(shè)置日期、檢驗倍數(shù)、滑動天數(shù)、產(chǎn)出方式和圖片分類。檢驗倍數(shù)為四分位距異常檢驗倍數(shù)，滑動天數(shù)為背景觀測數(shù)據(jù)滑動窗口長度，默認值分別為1.5倍和15 d；產(chǎn)出方式分為自動和手動，圖片分類分為中國和全球，用戶根據(jù)需要選擇設(shè)置。查詢結(jié)果分為圖片文件和地圖文件，圖片文件為電離層VTEC異常分布圖，全球VTEC異常分布圖時間間隔為1 h，中國區(qū)域VTEC異常分布圖時間間隔為2 h（見圖6、圖7）；地圖文件為圖片文件對應的矢量圖，用戶根據(jù)需要選擇文件并利用ArcGIS軟件修改圖片底圖、圖片標注、VTEC異常分布等設(shè)置（見圖8）。

圖6 系統(tǒng)運行主界面（全球）Fig.6 Main interface of system operation（Global）

圖7 系統(tǒng)運行主界面（中國）Fig.7 Main interface of system operation（China）

圖8 ArcGIS地圖文件處理界面（中國）Fig.8 Processing interface of ArcGIS map file（China）

圖形展示位于用戶操作界面的右側(cè)，根據(jù)查詢條件和查詢結(jié)果可分別查看全球VTEC異常分布、中國區(qū)域VTEC異常分布。

2.2 系統(tǒng)設(shè)置

系統(tǒng)設(shè)置包括TEC分布圖參數(shù)，GNSS數(shù)據(jù)下載，CODE數(shù)據(jù)下載三部分（見圖9）。其中，TEC分布圖參數(shù)需要設(shè)置檢驗倍數(shù)和滑動天數(shù)，默認值分別為1.5倍和15 d；GNSS數(shù)據(jù)下載通過FTP方式每天自動下載陸態(tài)網(wǎng)絡(luò)中國區(qū)域GNSS基準站觀測數(shù)據(jù)；CODE數(shù)據(jù)下載通過FTP方式每天自動下載星歷、時鐘等數(shù)據(jù)解算控制文件及全球電離層格網(wǎng)數(shù)據(jù)文件。

圖9 系統(tǒng)參數(shù)配置界面Fig.9 System parameters configuration interface

2.3 自動模式

自動模式為系統(tǒng)默認工作模式，按照系統(tǒng)設(shè)置的默認參數(shù)，每天自動下載中國區(qū)域GNSS觀測數(shù)據(jù)、星歷及時鐘等控制文件、CODE全球電離層格網(wǎng)數(shù)據(jù)和空間環(huán)境參數(shù)等數(shù)據(jù)，調(diào)用Bernese軟件完成數(shù)據(jù)自動解算，獲取中國區(qū)域電離層格網(wǎng)數(shù)據(jù)。結(jié)合CODE全球電離層格網(wǎng)數(shù)據(jù)，采用滑動四分位距法提取異常信息，自動產(chǎn)出全球VTEC、中國區(qū)域VTEC異常分布。用戶設(shè)置查詢條件，產(chǎn)出方式選擇自動，點擊查詢功能，選擇查詢結(jié)果中的圖形文件，即可查看電離層VTEC異常分布結(jié)果；點擊播放功能，可動態(tài)展示電離層VTEC異常分布變化；修改查詢條件，點擊生成功能，可重新生成電離層VTEC異常分布結(jié)果（見圖6，圖7）。

2.4 手動模式

手動模式為用戶選擇模式，手動導入中國區(qū)域GNSS數(shù)據(jù)，設(shè)置圖形結(jié)果保存目錄、檢驗倍數(shù)、滑動天數(shù)、開始和結(jié)束時間等參數(shù)，勾選數(shù)據(jù)解算和異常分布圖選項，點擊確定。系統(tǒng)根據(jù)用戶需要自動調(diào)用Bernese軟件完成數(shù)據(jù)自動解算、異常信息提取、圖形結(jié)果產(chǎn)出等功能（見圖10）。用戶設(shè)置查詢條件，產(chǎn)出方式選擇手動，點擊查詢功能，選擇查詢結(jié)果中的圖形文件，即可查看電離層VTEC異常分布結(jié)果，播放和生成功能與自動模式相同（見圖6，圖7）。

圖10 手動產(chǎn)出參數(shù)配置界面Fig.10 Configuration interface of manual output parameters

2.5 時序曲線

時序曲線功能基于數(shù)據(jù)庫中保存的空間環(huán)境參數(shù)（F10.7、Kp、Dst）及中國區(qū)域電離層格網(wǎng)數(shù)據(jù)，用戶通過設(shè)置經(jīng)緯度、異常檢驗方法、檢驗倍數(shù)、滑動天數(shù)、開始及結(jié)束時間等基本參數(shù)，設(shè)置F10.7、Kp、Dst、VTEC、ΔTEC、時間刻度等參數(shù)，點擊生成功能，系統(tǒng)根據(jù)用戶需要繪制空間環(huán)境參數(shù)及中國區(qū)域格網(wǎng)點VTEC時序曲線。其中，空間環(huán)境參數(shù)時序圖中的紅色直線為參數(shù)閾值；VTEC時序圖中的紅色曲線為觀測值，灰色曲線為VTEC觀測背景上下邊界閾值；ΔTEC為VTEC的擾動異常。同時，系統(tǒng)具備圖形導出、異常數(shù)據(jù)導出功能（見圖11）。

圖11 時序曲線參數(shù)配置界面Fig.11 Configuration interface of time sequence curve parameters

2.6 數(shù)據(jù)服務(wù)

系統(tǒng)根據(jù)用戶需求，提供數(shù)據(jù)轉(zhuǎn)換、數(shù)據(jù)導出和數(shù)據(jù)下載等功能。

數(shù)據(jù)轉(zhuǎn)換：用于將原始觀測數(shù)據(jù)壓縮文件（D.Z文件）轉(zhuǎn)換為原始觀測數(shù)據(jù)文件（O文件）。選擇觀測數(shù)據(jù)文件目錄，點擊確定，數(shù)據(jù)自動轉(zhuǎn)換（見圖12 a）。

數(shù)據(jù)導出：用于將起止時間內(nèi)的空間環(huán)境參數(shù)F10.7、Dst、Kp以文本文件格式導出，便于用戶查找和讀取。設(shè)置開始及結(jié)束時間，選擇數(shù)據(jù)保存目錄，勾選空間環(huán)境參數(shù)，點擊確定，完成數(shù)據(jù)導出（見圖12 b）。

數(shù)據(jù)下載：用于下載起止時間內(nèi)的星歷、時鐘等數(shù)據(jù)解算控制文件及CODE全球電離層格網(wǎng)數(shù)據(jù)文件。設(shè)置開始及結(jié)束時間，選擇數(shù)據(jù)保存目錄，點擊確定，完成數(shù)據(jù)下載（見圖12 c）。

圖12 數(shù)據(jù)服務(wù)界面Fig.12 Data service interface

2.7 任務(wù)查詢

用戶選擇查詢條件，在任務(wù)列表中可查詢自動或手動方式下任務(wù)執(zhí)行進程，若任務(wù)運行失敗，可查詢錯誤信息，方便用戶及時查看并處理（見圖13）。

圖13 任務(wù)查詢界面Fig.13 Task query interface

3 系統(tǒng)應用實例

該系統(tǒng)于2020年6月設(shè)計完成，并安裝部署試運行，系統(tǒng)運行情況良好。2019年6月17日四川長寧發(fā)生MS6.0地震，利用該系統(tǒng)手動模式產(chǎn)出全球VTEC、中國區(qū)域VTEC異常分布。結(jié)果顯示，2019年6月14日震中區(qū)域附近均出現(xiàn)明顯的電離層異常擾動，其中6月14日12：00（UT）異常最為明顯，異常持續(xù)時間長達8 h，電離層異常耦合的局部性、集中性和顯著性特征明顯（見圖14、圖15）。

圖14 2019年6月14日12：00 UT全球VTEC異常分布圖Fig.14 Distribution map of global VTEC anomalies at 12：00 UT on June 14th，2019

圖15 2019年6月14日12：00 UT中國區(qū)域VTEC異常分布圖Fig.15 Distribution map of VTEC anomalies in China at 12：00 UT on June 14th，2019

此外，利用時序曲線功能產(chǎn)出2019年6月1日至6月20日空間環(huán)境參數(shù)及震中附近格網(wǎng)點（28°N，104°E）電離層VTEC時序曲線（見圖16）。其中，空間環(huán)境參數(shù)時序曲線顯示6月14日空間環(huán)境參數(shù)（F10.7、Kp、Dst）均未超過閾值，排除長寧地震前太陽活動、磁暴等空間環(huán)境變化對電離層異常擾動的影響；VTEC時序曲線顯示6月14日震中附近格網(wǎng)點電離層VTEC異常擾動明顯，并從時間尺度驗證了該異常擾動與長寧地震發(fā)生的相關(guān)性。

圖16 2019年6月1日至6月20日空間環(huán)境參數(shù)及格網(wǎng)點（28°N，104°E）電離層VTEC時序曲線變化Fig.16 Variation of space environmental parameters and ionospheric VTEC time series curve at grid point（28 ° N，104 ° E）from June 1st to June 20th，2019

4 討論及結(jié)論

本文設(shè)計開發(fā)了電離層VTEC異常分布產(chǎn)出系統(tǒng)，實現(xiàn)了全球VTEC、中國區(qū)域VTEC異常變化實時監(jiān)測。該系統(tǒng)具有以下特點：

（1）系統(tǒng)自動產(chǎn)出全球VTEC、中國區(qū)域VTEC異常分布圖，產(chǎn)出結(jié)果延遲1~2 d，有助于全球、中國區(qū)域電離層VTEC異常實時監(jiān)測；

（2）系統(tǒng)可根據(jù)用戶需求手動產(chǎn)出指定時間段全球VTEC、中國區(qū)域VTEC異常分布圖，有助于分析歷史地震發(fā)生前后電離層VTEC異常擾動變化。

（3）系統(tǒng)提供空間環(huán)境參數(shù)、數(shù)據(jù)解算控制及電離層VTEC格網(wǎng)數(shù)據(jù)等文件下載功能，同時提供圖形產(chǎn)出結(jié)果修改編輯功能，便于用戶使用。

目前，電離層VTEC異常分布產(chǎn)出系統(tǒng)已成功應用到四川省地震局監(jiān)測預報日常工作，其圖形產(chǎn)出結(jié)果已參與地震監(jiān)測預報會商，實現(xiàn)了全球VTEC、中國區(qū)域VTEC異常變化實時監(jiān)測。為了測試地震-電離層監(jiān)測在川滇地區(qū)的應用效能，該系統(tǒng)仍需完善。未來，基于地基和空基電離層探測技術(shù)，開展天地一體化聯(lián)合觀測，進一步加強川滇地區(qū)電離層異常擾動監(jiān)測能力，分析地震前后電離層多參量異常變化，深化不同參量之間的關(guān)聯(lián)，完善地震-電離層耦合理論模型，嘗試將電離層地球物理異常信息用于川滇區(qū)域臨震監(jiān)測預報，這對于川滇地區(qū)震情判定具有重要意義。

致謝：CODE為該系統(tǒng)提供全球電離層格網(wǎng)數(shù)據(jù)；中國地震局地殼運動監(jiān)測工程研究中心、中國地震局地震預測研究所為該系統(tǒng)提供觀測數(shù)據(jù)。