简介这是一套面向电离层研究者、GNSS数据处理工程师及空间天气分析人员的Fortran科学计算源码专注于将标准RINEX格式的GNSS观测数据高精度反演为电离层总电子含量TEC并输出为专用GTEX格式有效支撑定位误差校正、电离层建模与太阳活动影响评估等关键任务。资源共42个文件含15个Fortran源文件.f实现核心算法如TEC反演、轨道读取、周跳修正、15个编译目标文件.o、3个Shell脚本.sh用于自动化流程控制以及Makefile、头文件、参数配置列表和详细README文档结构完整、模块清晰便于二次开发与教学实践。压缩包仅290KB轻量但功能完备。目前已有155人学习下载用户可直接获取可编译运行的完整工程涵盖从RINEX数据解析、卫星几何定位、伪距/相位组合TEC计算到GTEX格式写入的全链路实现并附带Julian日历转换、IGS周数计算等实用工具模块。 从源码读懂GNSS电离层处理这事放到现在依然有吸引力。RNX2GTEX这个名字常做GNSS数据处理的人一看就明白它做的是把RINEX观测文件转成TEX格式输出。这里的TEX不是LaTeX那套排版工具而是电离层TEC交换文件Total Electron Content Exchange的简写。整个程序的核心任务就是从双频GNSS观测值中提取总电子含量TEC再把逐卫星、逐历元的TEC结果写成标准格式供后续电离层建模、单频用户电离层改正、空间天气研究使用。我用这个程序处理过不少测站数据也把源码从头到尾翻过几遍坦白说这不是一个代码风格很“现代”的项目但它把GNSS电离层测量的整条链路用最朴素的方式跑通了。这篇内容我就从一个实际使用者的角度把RNX2GTEX涉及的物理原理、源码结构、编译运行细节、输出格式和常见坑都拆开讲一遍适合刚接触GNSS电离层数据处理、或者想通过老Fortran代码理解观测方程的同学参考。1. 先把名词理清楚RINEX、TEX、TEC各自扮演什么角色很多初学者看到RNX2GTEX这个软件名第一反应是去查怎么编译结果被RINEX文件名、TEX格式、TEC单位这些概念绕晕。我建议先花半小时把下面几个名词的关系理清后面看代码会轻松很多。1.1 TEC是什么GNSS为什么能测到电离层TEC的全称是Total Electron Content翻译过来就是总电子含量指沿信号传播路径上单位截面积柱体内的自由电子总数单位是电子数每平方米。平时我们更常用TECU做单位1 TECU等于10的16次方个电子每平方米。电离层里的自由电子会对GNSS信号产生折射延迟这个延迟的大小和信号频率的平方成反比。对频率为f的信号伪距上的电离层延迟近似为40.3乘STEC除以f的平方这里的STEC就是斜路径上的TEC。因为两个载波频率不一样同一颗卫星发出的L1和L2信号穿过电离层时延迟量不同通过对比两个频率的观测值就能把电离层影响单独分离出来。我在实际给学生讲的时候喜欢用一个类比电离层就像一块有色玻璃不同颜色的光穿过时速度不一样。GNSS接收机同时收到两个频率的信号相当于拿两束不同颜色的光同时穿过这块玻璃通过比较它们的到达时间差就能反推玻璃的“厚度”。这里的“厚度”就是TEC。1.2 从RINEX到TEXRNX2GTEX在流水线中的位置RINEX是GNSS观测数据的标准交换格式接收机厂商输出的原始数据经过转换后会以RINEX格式保存伪距、载波相位、多普勒等观测值。RINEX文件里包含的信息非常丰富但直接拿RINEX文件做电离层研究并不方便因为观测值里混着钟差、轨道误差、对流层延迟等一大堆无关量而且不同接收机输出的文件格式细节也不一样。TEX格式则是专门面向电离层TEC数据设计的交换格式它把处理好的TEC结果按站点、按历元、按卫星组织起来省去了用户重复做预处理的工作。RNX2GTEX就是连接RINEX和TEX的桥梁输入一个或多个RINEX观测文件经过质量检查和组合计算输出TEX格式的TEC时间序列。整个GNSS电离层处理链路大致是RINEX观测文件加精密星历和钟差经过预处理、周跳探测、无几何组合、相位平滑等步骤得到STEC再做硬件延迟校正和映射函数转换输出网格化的VTEC或直接输出STEC序列。RNX2GTEX覆盖的是从RINEX到STEC/TEC这一核心段落。1.3 单位转换和数量级怎么判断算出来的TEC是否合理用这个程序之前最好对TEC的数量级心里有数。中纬度地区平静电离层情况下天顶方向VTEC通常在10到50 TECU之间低纬赤道异常区可以到80甚至100 TECU以上高纬和夜间会低很多夜间经常只有几个TECU。如果算出来的结果在一个测站、一整天内变化几百上千个TECU那基本可以判断数据或处理方法出了问题。从几何关系上也要有概念。斜路径STEC一般是VTEC的好几倍仰角越低路径越长STEC越大。仰角30度左右时如果VTEC是30 TECUSTEC大概在50到60 TECU。程序输出的如果明显偏离这个范围就要回头检查组合系数、单位换算是哪里出了问题。还有一个非常实用的换算关系必须掌握伪距组合P1减P2的1米差异大约对应9.52 TECU具体推导依据是40.3乘以(1/f1的平方减1/f2的平方)的倒数f1是1575.42兆赫兹f2是1227.60兆赫兹。这个系数在验证程序输出时特别有用算完TEC以后可以反算一下P4残差看看是否在合理范围内。2. 源码核心逻辑STEC是怎么从观测值里算出来的RNX2GTEX本质上是把教科书里的双频电离层探测公式翻译成了Fortran代码。所以读这份源码数学上并不难难的是理解代码里各种变量、常量和文件操作背后对应的物理过程。我个人推荐在读源码之前先把下面几个关键逻辑在纸上推导一遍。2.1 无几何组合的推导和实现要点电离层探测的第一个关键组合叫无几何组合geometry-free combination也叫电离层残差组合。对于伪距组合形式是P4等于P1减P2对于载波相位组合形式是L4等于L1减L2。这个组合能消掉卫星钟差、接收机钟差、对流层延迟、几何距离等与频率无关的项剩下的主要就是电离层延迟差异和硬件延迟偏差。在源码里这个组合通常不是直接算两个观测值相减就完事还要考虑P1和P2的观测值类型。有的接收机输出C1和P2有的输出P1和P2还有的只有C1和C2不同组合对应的硬件延迟偏差不一样代码里一般会有对应的分支判断。这也是为什么源码里出现一长串if条件判断的原因。从组合值换算到STEC的时候符号问题特别容易搞晕。如果程序里写的是P4等于P1减P2那么P4大约等于负的40.3乘STEC乘以(1除以f1平方减1除以f2平方)换算成TEC需要乘一个负系数如果程序里用的是P2减P1符号就反过来。我看到不少人在看源码时卡在这里最后发现是符号理解反了。2.2 载波相位平滑伪距为什么需要代码里怎么体现伪距观测的噪声比较大尤其C/A码和P码噪声水平通常是分米级甚至米级直接用伪距组合算STEC结果会非常毛糙。载波相位观测的噪声小得多一个量级的差距但相位观测值存在整周模糊度差分后依然有一个未知常数偏差无法直接给出绝对TEC。解决办法是用相位组合的变化量来平滑伪距组合的绝对值这就是经典的载波相位平滑伪距算法。具体实现思路是先对L4做周跳检测把连续的、没有周跳的弧段找出来在每个弧段内计算(L4减去P4)的时间平均这个平均值包含了模糊度和硬件延迟的综合常数然后用P4观测值加上这个平均值得到平滑后的STEC。程序里通常会有一个循环逐历元处理遇到周跳就重新初始化平滑窗口。我在源码里看到平滑窗口长度设置时不同版本差别比较大。窗口太短平滑效果差窗口太长又容易把电离层的真实变化也抹平了。一般长弧段取20到30分钟比较稳妥但如果电离层活跃、TEC变化剧烈窗口要适当缩短。这个参数值得根据你的数据和研究目的多试几组。2.3 DCB偏差处理源码的边界在哪里这里要特别注意伪距组合P4里面除了电离层延迟还包含卫星差分码偏差和接收机差分码偏差统称DCB。即使做了相位平滑DCB依然保留在结果里。如果不修正STEC会出现一个系统性偏置中纬度地区这个偏置折算下来少则几个TECU多则十几二十个TECU对电离层建模影响很大。RNX2GTEX这个层级的程序通常是不做DCB估计的。它的定位是把原始观测换算成不含几何项的电离层组合量并输出成TEXDCB修正一般留给后续处理链路比如用IGS或者CODE发布的DCB产品做后处理剔除。源码里可能预留了DCB文件的读取接口也可能完全没有要看具体版本。你在把TEX数据用于定量研究之前一定要确认DCB处理在哪一步完成否则结果会整体偏移。这条边界一定要搞清楚。如果你拿到一个RNX2GTEX输出的TEX文件第一件事不是画图看趋势而是要问这份数据是原始STEC还是已经做了DCB修正做了卫星端修正还是接收机端也修正了。我见过有人直接拿未修正DCB的数据做VTEC地图出来的图上整个测区都有一个固定偏置事后排查了半天才发现问题出在数据源头。3. 把老Fortran代码跑起来编译与运行的全过程RNX2GTEX是Fortran写的一般是Fortran 77风格的固定格式代码。这种老代码在今天的Linux环境上编译通常会遇到一些小问题但解决起来也不难关键是要知道几个典型的坑。3.1 源码文件组成与gfortran编译命令从源码库拿到的RNX2GTEX可能是一个单独的.f文件也可能是主程序加若干子程序文件的结构。以常见版本为例主程序文件名一般是RNX2GTEX.f里面包含若干个subroutine比如读取RINEX文件头的子程序、读取观测记录的子程序、计算TEC的子程序、写出TEX文件的子程序每个子程序对应一个独立的处理阶段。用gfortran编译时最简单的命令是gfortran -O2 -ffixed-line-length-132 -o rnx2gtex RNX2GTEX.f如果源码拆成了多个文件就全部列在命令后面gfortran -O2 -ffixed-line-length-132 -o rnx2gtex RNX2GTEX.f SUBRTN.f CONST.f-ffixed-line-length-132这个选项容易忽略但经常是编译报错的关键。很多老Fortran代码的注释和续行标志在固定格式下默认只认72列超过的部分会被忽略如果源码里某一行比较长不调整行长限制的话编译时会报一堆莫名其妙的语法错误。如果是64位Linux系统一般不用加额外选项就能编译但个别版本会用到非标准的库函数这时要根据报错信息去源码里查具体是哪个函数再决定是替换实现还是增加兼容代码。我不建议一开始就改动源码逻辑先试着原样编译遇到问题再对症下药。3.2 RINEX文件命名约定这个坑最容易栽RNX2GTEX对输入文件的命名有约定基本上遵循标准RINEX文件命名规则前四个字符是站名缩写第五到第七个字符是年积日第八个字符是日内文件序号第九第十个字符是年份点后面两位是文件类型标识。比如abmf0010.15o表示abmf站、年积日第001天、序号0、2015年、观测文件。这个命名约定是程序正确运行的前提因为很多老程序不会做太智能的文件名解析而是直接按位置截取字符串来提取站名、年份和年积日。如果你把文件重命名成test_obs.15o之类的名字程序要么报错要么给出完全错误的输出。我建议在运行前把所有输入文件统一改成标准命名并放在同一个目录下文件名全部用小写或者全部用大写不要混用。有的程序在文件系统大小写敏感的环境里对文件名的判断很严格稍微不一致就会找不到文件。3.3 运行、交互输入与TEX输出验证编译成功后运行方式通常有两种取决于你拿到的版本。老版本一般是交互式提示输入运行后程序会问你要RINEX文件名有的版本支持命令行参数直接指定文件名比如./rnx2gtex abmf0010.15o运行结束后目录下会多出一个TEX文件。这时不要急着拿去用先打开文件看一下头部信息对不对站名是否与输入一致历元数是否合理有没有出现大量零值或负值。我一般习惯用head命令看前几十行再用awk统计一下STEC列的数值范围如果最小值是负几十、最大值是正几百说明数据预处理环节大概率有问题。如果输出文件是空的或者程序中途崩溃最优先检查的永远是RINEX文件本身可以先确认它能否被其他常用软件正常读取比如用teqc或者gfzrnx做一下质量检查排除RINEX文件损坏的可能再回头查程序参数设置。4. 输出文件长什么样TEX格式解析与Python后处理TEX格式是RNX2GTEX的输出也是后续处理的起点。虽然不同版本输出格式存在差异但结构上大体一致理解之后用脚本解析并不难。4.1 TEX文件头部和正文结构一个典型的TEX输出文件头部通常包含生成程序的标识、站点名称、站点坐标、数据的时间范围、观测文件的相关信息等内容。正文部分按历元组织每个历元下列出可见卫星的TEC值可能带有卫星编号、仰角等信息。对于源码里的写语句建议逐个对照看。有的版本输出STEC有的版本输出VTEC有的版本会同时输出仰角供你后续自己换算。源码中每个格式描述符对应的内容最好对照RINEX文件和程序内部变量Name来确认不要只看文件后缀就默认是VTEC。4.2 用Python快速解析并画一条VTEC时间序列TEX文件虽然可以直接用文本编辑器打开看但要做时间序列分析或画图还是得写脚本。下面给一个很基础的Python解析示例具体列位置需要根据你那个版本的写语句调整import matplotlib.pyplot as plt records [] with open(abmf0010.tex, r) as f: for line in f: if line.startswith(RNX2GTEX OUTPUT): parts line.split() station parts[2] elif line.strip() and line[0].isdigit(): doy int(line[0:3]) sec float(line[4:14]) prn int(line[15:17]) stec float(line[18:27]) records.append((doy, sec, prn, stec)) for prn in sorted(set(r[2] for r in records)): data [r for r in records if r[2] prn] times [r[1] / 3600.0 for r in data] values [r[3] for r in data] plt.plot(times, values, labelfPRN {prn}) plt.xlabel(Hour of Day) plt.ylabel(STEC (TECU)) plt.legend() plt.show()这段代码把每个卫星的STEC按小时画成一条线可以快速看出各卫星之间的系统偏差是否正常。如果某颗卫星整体比别的卫星高出一截大概率是卫星DCB没有修正如果出现锯齿状跳变说明周跳处理有问题。4.3 数据后处理的几点建议我处理TEX数据时习惯做三步检查。第一步看单颗卫星连续弧段是否平滑有没有突跳第二步把所有卫星的STEC映射到天顶方向做VTEC按站点看日变化曲线是否合理正常情况中午高、夜间低第三步用同一天相邻测站的数据做交叉验证如果两个测站距离很近VTEC应该高度一致。映射STEC到VTEC时经典做法是除以仰角的正弦值近似映射函数。但要注意低仰角时映射函数误差很大一般会把仰角低于10度或15度的数据去掉再换算。RNX2GTEX如果本身不输出仰角你可能还需要从RINEX文件或者星历计算里补上这个信息这也是不少人在后续处理时额外写模块的原因。5. 常见问题与排查技巧实录这个项目我用下来真正运行顺利的情况其实不多大多数时间都在跟各种细节较劲。下面这些问题几乎每个使用RNX2GTEX的人都会碰到我按阶段整理成一张排查表方便你遇到问题时直接对照。5.1 编译阶段的问题编译是第一个拦路虎也是最容易劝退新手的环节。我见过最多的报错有两类一类是固定格式行长问题用-ffixed-line-length-132基本能解决另一类是源码里用了非标准的扩展语法比如Tab开头的代码行、超过72列的续行、或者比较老的Fortran特性gfortran默认模式下不接受。遇到这类问题先把编译器的警告信息完整看一遍重点找第一个报错位置因为后续报错往往是连锁反应。如果某个语法确实是老扩展最简单的处理是把那几行改写成标准Fortran 77语法。我不建议为了编译通过而关掉所有警告容易埋下运行时隐患。5.2 运行阶段的数据问题程序编译通过只是开始运行结果异常的排查才更耗时间。以下几种情况我都在实际数据里遇到过第一种输出全部为0或者全部为负值。最常见的原因是观测文件里没有程序期望的观测值类型比如程序默认读P1和P2但现代接收机只输出C1和C2或者双频数据里有大量L2观测值缺失。这时需要检查RINEX文件头里的观测类型列表确认实际包含哪些信号。第二种输出的TEC在长时间段内整体偏置。这是DCB没修正的典型特征尤其是接收机端DCB在同一台接收机的数据里是常数很容易被误认为是真实电离层变化。如果相邻两天的数据在同一时刻都有固定差异大概率就是DCB问题。第三种单颗卫星数据在某个历元突然跳变之后又恢复正常这通常是观测数据本身存在跳变或者周跳漏检。可以检查平滑窗口的周跳检测阈值是否需要调整阈值设得太松会把小周跳放过去设得太紧又会频繁重置平滑窗口导致结果噪声增大。5.3 结果异常排查速查表现象可能原因处理建议输出全部为零或负数观测值类型不匹配、无P1/P2组合查看RINEX头文件观测类型确认输入数据STEC整体偏置不同卫星各自的基线不同未做卫星端/接收机端DCB修正接入DCB产品做后处理修正个别卫星弧段突跳周跳漏检、平滑窗口参数不合理调低周跳检测阈值缩短平滑窗口处理后VTEC夜间出现负值平滑噪声、低仰角数据污染提高截止仰角检查高度角输入程序运行时提示文件无法打开文件名不符合RINEX命名规则按ssssdddf.yyt格式重命名输入文件编译时出现unclassifiable statement固定格式行长或非标准扩展语法加编译选项修改老式语法输出文件只有头部没有正文RINEX观测记录读取失败检查RINEX文件完整性先用teqc等工具质检表中列出的问题覆盖了我看到的大多数求助。实际上每次排查这些问题都能加深对程序和数据格式的理解。你把这个表存下来遇到问题先按图索骥能省不少时间。6. 源码里值得多读几遍的几个地方RNX2GTEX代码量不大但里面有几个片段非常值得反复读。第一个是RINEX头部解析部分这里能看到程序如何处理不同版本的RINEX格式老代码往往用许多分支来兼容不同年代的文件格式读这部分能学到不少处理历史数据的经验。第二个是周跳检测和平滑窗口的实现。这个模块直接决定输出STEC的质量也是后来人改动最多的部分。有的版本用L4变化量超过固定阈值来判断周跳有的版本用相邻历元差分的统计量动态设定阈值两种方式各有优劣。你完全可以在理解原逻辑后把平滑部分替换成更新的算法比如基于卡尔曼滤波的方式。第三个是TEX文件输出的写语句。这部分能直观看出程序作者对输出格式的考虑哪些信息被保留哪些信息被丢弃背后都有取舍。如果你要做更细致的分析可能需要在这里增加输出内容比如加上方位角、高度角、信噪比等。我在读这份源码过程中最大的体会是老程序虽然界面不友好、代码风格不现代但它的逻辑非常直接几乎没有多余的设计。这种直接性反而让学习变得容易你能清楚地看到每一步算子对应教科书里的哪个公式。对于想理解GNSS电离层处理全流程、或者想写自己的TEC处理工具的人来说RNX2GTEX是一份非常合适的入门源码。如果后续想扩展它的能力我建议优先考虑这几个方向一是增加对Galileo和北斗观测值的支持老程序最初主要是针对GPS设计的二是把DCB修正模块直接集成进去省去后处理环节三是增加输出精度因子和高度角信息方便下游做质量加权。改动过程中注意保持原有输出格式的兼容性因为很多下游工具默认了TEX的旧版结构。最后分享一个实际工作中的小习惯每次处理一个新的测站或新一天的数据时我会保留程序输出的原始TEX文件用脚本自动生成一份包含最大值、最小值、均值、有效数据率的统计报告。这样处理大量数据时能快速挑出异常那天。GNSS电离层数据处理麻烦往往不在程序本身而在数据质量的把控上这个习惯帮我省了不少排查时间。本文还有配套的精品资源点击获取