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)共享的数据支持．