Linux Shell在電離層TEC格網數據提取和分析中的應用

李涌濤,李建文,顧晨鐘,陳晨,張碩,車通宇

(1.西安測繪總站,陜西 西安 710054;2.信息工程大學 地理空間信息學院,河南 鄭州450001;3.西昌衛星發射中心,四川 西昌 615000)

0 引 言

電離層總電子含量(TEC)作為描述電離層形態和結構的重要參數,在電離層延遲改正及電離層變化規律等相關研究中起著非常重要的作用[1-2]．全球電離層格網(GIM)數據,一方面因其將全球按經緯度格網化,并給出了相應格網點的電離層TEC值,易于內插計算其他經緯度位置的TEC值,從而成為目前為用戶提供電離層延遲改正通用的數據格式[3-5]．另一方面,電離層作為空間大氣的組成部分,研究其變化規律對掌握空間大氣變化規律具有重要意義[6]．由于電離層受到太陽及其他日地活動的影響,從而發生復雜的物理變化,致使電離層成為了一個不斷變化且復雜的開放性系統[7-8],而電離層總電子含量TEC則成為了研究全球或區域電離層變化規律的重要依據．GIM數據,在電離層TEC分析研究中的應用主要分為:全球不同電離層格網模型精度分析、不同系統或者多系統全球電離層格網模型分析、不同緯度和不同區域電離層TEC變化特性分析以及不同GIM產品精度分析等[9-12]．

Linux系統因為其具有高穩定性、健壯、安全、高性能和自由開源等特性,再輔以Shell、Python和Perl等腳本語言進行操作運行,從而被世界范圍內大多數科研機構所采用,其中包括國際全球導航衛星系統(GNSS)服務(IGS)分析中心和數據中心以及國際GNSS監測評估系統(iGMAS)數據中心和分析中心[13-16]．Linux Shell是一種具有命令解釋功能的程序,也是連接用戶和Linux內核并提供用戶使用操作系統的接口,Linux Shell自身就是一個解釋型的程序設計語言[17-19]．Shell腳本是集成了多條命令或程序語言的可執行程序文件,且Shell在grep、awk和sed等眾多Linux系統命令的支撐下,使其在進行文本處理時相較于其他中、高級編程語言,具有編寫快速簡單、執行高效、容易修改維護和可移植性好等優勢[20-22]．

本文基于歐洲定軌中心(CODE)發布的2017年全球電離層GIM格網數據,利用Shell腳本展示并實現了對GIM相應數據的提取和分析等,其中包括GIM中電離層TEC小時平均值、日平均值、不同緯度平均值、特定點或區域TEC的提取計算,TEC時間和空間上差值及相應數值的提取計算,以及利用通用制圖工具(GMT)軟件批量繪制全球及區域電離層TEC圖等．

1 數據準備

選取CODE提供的2017年365天GIM數據,數據下載網址為(ftp://ftp.aiub.unibe.ch/CODE/),GIM數據采用國際通用格式[23-24]．GIM格網的緯度范圍為87.5°～-87.5°,緯度差為2.5°,經度范圍為-180°～180°,經度差為5°,即分辨率為2.5°×5°,將全球劃分為71×73個格網點,如圖1所示;采樣時間間隔為1 h,每天從0:00-24:00共25幅GIM格網數據[25]．在Linux系統中,利用Shell腳本對所需的GIM數據進行提取、整理、計算并按用戶指定的格式輸出．

圖1 71×73 GIM格網圖

2 GIM基礎數據整理

2.1 GIM數據格式

以CODE 2017年12月26日,年積日(DOY)360天的GIM為例,GIM數據文件格式如圖2所示．文件每行限定為80個字符,以“END OF HEADER”為分界,“END OF HEADER”之上為文件頭,是GIM的基本信息說明,之下為格網點TEC值信息．“START OF TEC MAP”為GIM圖幅起始的編號,CODE的GIM采樣間隔為1 h,一天共有1～25幅;“EPOCH OF CURRENT MAP”為當前圖幅的歷元“年 月 日 時 分 秒”;“LAT/LON1/LON2/DLON/H”為當前圖幅的格網位置信息,即“緯度/第一個格網的經度/最后一個格網經度/經度差/電離層TEC值高度”;“LAT/LON1/LON2/DLON/H”以下為與經緯度對應格網點的TEC值(綠色矩形框內),由于受每行字符數所限,TEC值共分為5行,前4行每行16個,第5行9個,即同一緯度線上,經度差為5°的73個格網點TEC值．后面的數據為緯度差為2.5°的相應緯度線格網點TEC值．“END OF TEC MAP”為當前TEC圖幅結束,與起始編號一致,按歷元順序第1～25幅以此類推．TEC數據之后是相對應格網點的均方根值(RMS),數據格式與TEC一致,起止標志字符串分別為“START OF RMS MAP”和 “END OF RMS MAP”,以此和TEC值區別,TEC和RMS單位都是0.1 TECU(1TECU=1016個電子/m2)[7]．

圖2 GIM格網數據格式

2.2 GIM格網TEC提取和整理

利用Shell腳本提取GIM格網數據中的TEC值,整理成71×73的數列,并以“格網范圍(全球或區域)-TEC(或RMS)-四字符機構名稱-三字符DOY-小時歷元h-兩位年份i-guocheng.txt”的自定義格式存儲單幅GIM格網TEC值(如:glob-TEC-CODE-360-14 h-17i-guocheng.txt),以備之后進一步計算相關數據,以CODE 2017年DOY 360的GIM格網數據為例,輸入文件名稱為“CODG3600.17I”．Shell腳本如代碼1)～8)所示:

1) starTEC=＄(sed-n'1 (55個空格) START OF TEC MAP/=' CODG3600.17I) #注釋:提取第1幅TEC開始行號賦值給變量＄starTEC．

2) endTEC=＄(sed-n'1 (55個空格) END OF TEC MAP/=' CODG3600.17I) #注釋:提取第1幅TEC結束行號賦值給變量＄ endTEC．

3)sed-n"＄{starTEC},＄((endTEC))p" CODG3600.17I> 001 #注釋:提取＄starTEC和＄endTEC行號之間的數據存儲為臨時文件001．

4) sed-i '/START OF TEC MAP/d' 001 #注釋:刪除文件001中包含字符串“START OF TEC MAP”的行．

5) sed-i '/EPOCH OF CURRENT MAP/d' 001 #注釋:刪除文件001中包含字符串“EPOCH OF CURRENT MAP”的行．

6) sed-i '/LAT/LON1/LON2/DLON/H/d' 001 #注釋:刪除文件001中包含字符串“LAT/LON1/LON2/DLON/H”的行．

7) sed-i '/END OF TEC MAP/d' 001 #注釋:刪除文件001中包含字符串“END OF TEC MAP”的行．

8) awk 'ORS=NR%5?" ":" "{print}' 001> glob-TEC-CODG-360-1h-17i-guocheng.txt #注釋:將文件001中的數據每5行合并為1行,即同一緯度的TEC值為一行．

其中1)～3)從GIM中提取了只含有第1幅TEC的數據,4)～7)刪除了多余的解釋性字符串,只保留了每個緯度5行的TEC數據,8)將以上提取的數據每5行合并為一行,生成對應GIM經緯格網的TEC值71×73的數列,并存儲第1 h的TEC數據文件名為“glob-TEC-CODG-360-1h-17i-guocheng.txt”(以下簡寫為“CODG3601h.txt”)的過程文件,以備后用．不同小時和不同天的過程文件可用for循環進行生成．

3 GIM格網TEC相關數值計算

基于以上提取整理的每幅(即每小時)71×73的TEC值數列,可以分別提取計算全球或者區域內每小時、每天的TEC均值,以分析研究TEC的小時和日變化;計算經度互差、緯度互差,以分析TEC在空間經緯度方向上的變化;計算不同GIM對應格網點TEC的差值,以分析TEC在時間的變化率及變化范圍等．

3.1 TEC小時均值和日均值的計算

3.1.1 全球范圍

以上提取生成的CODG3601h.txt為全球范圍GIM格網按小時劃分的TEC值,因此可以借助直接用于計算全球范圍內的TEC小時均值,其中用到awk命令,如代碼9):

9) awk '{for(i=1;i<=NF;i++) x+=＄i}END{print x/NR/NF/10}' CODG3601h.txt? AVG.txt #注釋:可用“?”符號重定向,將一天25幅按每小時計算的平均值存儲在同一文件名為AVG.txt的文件中．

計算TEC日均值有兩種方法．方法一:用25幅GIM小時平均值求均值;方法二:將一天25幅GIM文件用cat命令合成一個文件,按代碼9)計算全天TEC日平均值．全球電離層TEC小時均值(以DOY 360～365為例)如圖3所示.

圖3 DOY 360～365全球電離層TEC小時均值

將全球南北半球按緯度30°劃分為高緯度帶(60°～90°)、中緯度帶(30°～60°)和低緯度帶(0°～30°),分別計算不同緯度的TEC小時或日均值．

CODG3601h.txt文件中每一行的數據即為緯度線上的TEC值,因此可以根據不同緯度對應不同的行數,提取不同緯度帶的TEC值,并用代碼9)進行均值計算,如代碼10)～14)所示．

10) sed-n "1,12p" CODG3601h.txt> hig-lat.txt

11) sed-n "60,71p" CODG3601h.txt? hig-lat.txt #注釋:CODG3601h.txt文件中1～12、60～71行對應南北緯高緯度帶(60°～90°)．

12) sed-n "13,24p" CODG3601h.txt> mid-lat.txt

13) sed-n "48,59p" CODG3601h.txt? mid-lat.txt #注釋:CODG3601h.txt文件中13～24、48～59行對應南北緯中緯度帶(30°～60°)．

14) sed-n "25,47p" CODG3601h.txt> low-lat.txt #注釋:CODG3601h.txt文件中25～47行對應南北緯低緯度帶(0°～30°)．

2017年365天全球不同緯度帶電離層TEC日均值如圖4所示，全球不同緯度帶電離層TEC小時均值(以DOY 360為例)如圖5所示．

圖4 2017年全球電離層TEC日均值

圖5 DOY 360 全球不同緯度帶電離層TEC小時均值

3.1.2 自定義區域范圍

CODG3601h.txt文件已經將TEC值提取并按經緯度格式化,經緯度與行列一一對應,可以根據該對應關系,提取自定義區域范圍(包括單一格網點)的TEC值,即按以上方法先提取出所需緯度范圍對應行的數據,再提取所需的經度范圍對應列的數據,最后進行相關計算．2017年DOY 360～365天自定義區域(15°N～55°N,75°E～135°E)的電離層TEC日均值如圖6所示,其中經緯度分別對應CODG3601h.txt文件中14 ～30行和52 ～64列．

圖6 DOY 360～365自定義區域電離層TEC日均值

3.2 TEC時空互差計算

3.2.1 TEC經緯度空間互差

在同一小時的GIM格網中,通過不同經度的TEC值互差和不同的緯度的TEC值互差,可研究電離層TEC在經緯度方向上的變化范圍及規律．以DOY 360第1 h數據,經度方向上相鄰經度差為5°的格網點TEC值互差為例,因互差涉及不同行列數據,故用到兩個for嵌套循環,如代碼15)～26)所示．

15) for ((row=1;row<=71;row++)) #注釋:第1個for作行數據循環

16) do

17) for ((col=1;col<=72;col++)) #注釋:第2個for作列數據循環

18) do

19) let nextcol=＄{col}+1

20) a=`awk 'NR==" '＄row' " {print ＄' "＄col" '}' CODG3601h.txt`

21) b=`awk 'NR==" '＄row' " {print ＄' "＄ nextcol " '}' CODG3601h.txt`

22) let c=＄b-＄a #注釋:同一緯度上相鄰經度TEC值作差

23) echo ＄c?lngdiff.txt #注釋:將所有差值作為一列輸出到文件lngdiff.txt．

24) done

25) done

26) awk 'ORS=NR%72?" ":" "{print}' lngdiff.txt> lngdiff-3601h.txt #注釋:將lngdiff.txt文件每72行合并為一行,即整理為71×72的經度差互差數列,命名為lngdiff-3601h.txt．

lngdiff-3601h.txt以備之后從中提取最值和計算均值所用．2017年經度差5°差值最大值、最小值和絕對值的平均值如圖7所示．

圖7 2017年經度差5°差值的最大值、最小值和絕對值的平均值

3.2.2 TEC經緯度時間互差

在不同小時的的GIM格網中,通過對應經緯度的TEC值互差,可研究電離層TEC在不同時間間隔內的變化范圍及規律．

以DOY 360第1 h和第2 h數據,對應格網點TEC值互差為例,因互差涉及不同行列數據,故用到兩個for嵌套循環,如代碼27)～ 37)所示．

27) for ((row=1;row<=71;row++)) #注釋:第1個for作行數據循環

28) do

29) for ((col=1;col<=73;col++)) #注釋:第2個for作列數據循環

30) do

31) a=`awk 'NR==" '＄row' " {print ＄' "＄col" '}' CODG3601h.txt` #注釋:此處CODG3601h.txt做為awk的輸入文件．

32) b=`awk 'NR==" '＄row' " {print ＄' "＄col" '}' CODG3602h.txt` #注釋:此處CODG3602h.txt做為awk的輸入文件．

33) let c=＄b-＄a #注釋:不同文件對應GIM格網點TEC值作差

34) echo ＄c?timediff.txt #注釋:將所有差值作為一列輸出到文件timediff.txt．

35) done

36) done

37) awk 'ORS=NR%73?" ":" "{print}' timediff.txt> timediff.txt-3601-2h.txt #注釋:將timediff.txt文件每73行合并為一行,即依舊整理生成71×73的經緯度1 h時間差互差數列,命名為timediff.txt-3601-2h.txt．

timediff.txt可為之后從中提取最值和計算均值所用．2017年時間差為1 h對應GIM格網點互差值的最值如圖8所示，2017年時間差為1 h、2 h、4 h和8 h的對應GIM格網點互差值絕對值的平均值如圖9所示．

圖8 2017年時間差1 h差值的最大值和最小值

圖9 不同時間間隔差值絕對值的平均值

4 基于GMT繪制電離層TEC云圖

通用制圖工具(GMT)是一種開源的地圖繪制工具,其遵循LGPL(GUN Lesser General Public License)開源協議,自帶80多種制圖工具,并支持多種輸出方式,可在Linux系統上輔以Shell高效運行,被廣泛應用于全球地學界各種成果的表達[26-28]．基于GMT繪制電離層TEC圖,最核心的步驟是根據提取整理好的格網數列文件CODG3601h.txt,生成GMT要求的文本數據,數據格式為每行三列,每一列分別為:緯度、經度、TEC值,以空格分隔．

根據CODG3601h.txt文件生成電離層TEC云圖GMT格式文本,如代碼38)～51)所示．

38) for ((row=1; row <=71; row++)) #注釋:第1個for作行數據循環．

39) do

40) let nn=＄ row-1

41) alt=`echo 87.5-＄nn*2.5|bc` #注釋:根據行數生成對應緯度,間隔2.5°．

42) r=`sed-n "＄{row}p" CODG3601h.txt`

43) a=(＄r)

44) for ((col=1; col<=73; col++)) #注釋:第2個for作列數據循環．

45) do

46) let mm=＄col-1

47) lng=`echo-180+＄mm*5.0|bc` #注釋:根據行數生成對應經度,間隔5°．

48)data=`echo"scale=2;＄{a[＄col]}/10"|bc` #注釋:根據CODG3601h.txt文件生成對應經緯度的電離層TEC值．

49) echo ＄alt" "＄lng" "＄data? CODG3601h-gmt.txt #注釋:按“緯度 經度 TEC值”格式輸出．

50) done

51) done

最終生成的CODG3601h-gmt.txt文件,如圖10所示．

按照GMT作圖要求,加入GMT各種作圖參數設置,并結合CODG3601h-gmt.txt文件,最終GMT繪制得到的是PostScript(PS)格式文件,即“.ps”格式的文本文件,PS是一種頁面描述性編程語言,兼具代碼的可讀性和矢量圖的縮放性．但是“.ps”文本文件并不方便閱讀,用GMT的“ps2raster”命令將其轉換為“.png”圖片格式,如代碼52)所示．

圖10 電離層TEC云圖GMT格式文本數據

52) ps2raster CODG3601h-ion.ps-Tg #注釋:g表示轉換為“.png”格式．

以2017年DOY 360 UTC 2:00 、4:00、6:00的GIM格網數據為例,生成全球電離層TEC圖如圖11所示,自定義區域范圍的電離層TEC圖如圖12所示．

圖11 2017DOY 360全球電離層TEC圖

圖12 2017DOY 360自定義區域電離層TEC圖

5 結束語

GIM電離層TEC格網數據對電離層TEC的研究起著很重要的作用．本文基于Linux Shell腳本編寫簡單快速、執行高效和容易修改維護等優勢,結合電離層TEC數據分析需求,利用Linux Shell腳本對GIM電離層TEC格網數據進行了基礎TEC數據提取和整理．在此基礎上實現了全球和自定義區域電離層TEC小時均值和日均值的計算,經緯度TEC值時空互差及最值、絕對值的平均值計算,以及利用GMT軟件實現了不同區域范圍內電離層TEC云圖的繪制．以上工作對全球或區域的電離層TEC的周年變化、季節變化、周日變化規律,電離層TEC的時空變化特性等相關分析研究具有重要的參考意義．

致謝:衷心感謝歐洲定軌中心(CODE)共享的數據支持．