Linux Shell在電離層TEC格網(wǎng)數(shù)據(jù)提取和分析中的應(yīng)用

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

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

0 引 言

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

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

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

1 數(shù)據(jù)準(zhǔn)備

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

圖1 71×73 GIM格網(wǎng)圖

2 GIM基礎(chǔ)數(shù)據(jù)整理

2.1 GIM數(shù)據(jù)格式

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

圖2 GIM格網(wǎng)數(shù)據(jù)格式

2.2 GIM格網(wǎng)TEC提取和整理

利用Shell腳本提取GIM格網(wǎng)數(shù)據(jù)中的TEC值,整理成71×73的數(shù)列,并以“格網(wǎng)范圍(全球或區(qū)域)-TEC(或RMS)-四字符機構(gòu)名稱-三字符DOY-小時歷元h-兩位年份i-guocheng.txt”的自定義格式存儲單幅GIM格網(wǎng)TEC值(如:glob-TEC-CODE-360-14 h-17i-guocheng.txt),以備之后進(jìn)一步計算相關(guān)數(shù)據(jù),以CODE 2017年DOY 360的GIM格網(wǎng)數(shù)據(jù)為例,輸入文件名稱為“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結(jié)束行號賦值給變量＄ endTEC．

3)sed-n"＄{starTEC},＄((endTEC))p" CODG3600.17I> 001 #注釋:提取＄starTEC和＄endTEC行號之間的數(shù)據(jù)存儲為臨時文件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中的數(shù)據(jù)每5行合并為1行,即同一緯度的TEC值為一行．

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

3 GIM格網(wǎng)TEC相關(guān)數(shù)值計算

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

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

3.1.1 全球范圍

以上提取生成的CODG3601h.txt為全球范圍GIM格網(wǎng)按小時劃分的TEC值,因此可以借助直接用于計算全球范圍內(nèi)的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文件中每一行的數(shù)據(jù)即為緯度線上的TEC值,因此可以根據(jù)不同緯度對應(yīng)不同的行數(shù),提取不同緯度帶的TEC值,并用代碼9)進(jìn)行均值計算,如代碼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行對應(yīng)南北緯高緯度帶(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行對應(yīng)南北緯中緯度帶(30°～60°)．

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

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

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

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

3.1.2 自定義區(qū)域范圍

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

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

3.2 TEC時空互差計算

3.2.1 TEC經(jīng)緯度空間互差

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

15) for ((row=1;row<=71;row++)) #注釋:第1個for作行數(shù)據(jù)循環(huán)

16) do

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

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 #注釋:同一緯度上相鄰經(jīng)度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的經(jīng)度差互差數(shù)列,命名為lngdiff-3601h.txt．

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

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

3.2.2 TEC經(jīng)緯度時間互差

在不同小時的的GIM格網(wǎng)中,通過對應(yīng)經(jīng)緯度的TEC值互差,可研究電離層TEC在不同時間間隔內(nèi)的變化范圍及規(guī)律．

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

27) for ((row=1;row<=71;row++)) #注釋:第1個for作行數(shù)據(jù)循環(huán)

28) do

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

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 #注釋:不同文件對應(yīng)GIM格網(wǎng)點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的經(jīng)緯度1 h時間差互差數(shù)列,命名為timediff.txt-3601-2h.txt．

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

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

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

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

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

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

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

39) do

40) let nn=＄ row-1

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

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

43) a=(＄r)

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

45) do

46) let mm=＄col-1

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

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

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

50) done

51) done

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

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

圖10 電離層TEC云圖GMT格式文本數(shù)據(jù)

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

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

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

圖12 2017DOY 360自定義區(qū)域電離層TEC圖

5 結(jié)束語

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

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