简介这是一套面向电离层物理、GNSS信号处理及空间天气研究领域的Fortran源码工具专用于将标准RINEX格式的GNSS观测数据如GPS/GLONASS高精度反演为电离层总电子含量TEC并输出为GTEX专用格式服务于定位误差校正、电离层建模与太阳活动影响分析等科研与工程场景。资源共42个文件含15个编译目标文件.o、15个核心Fortran源码.f涵盖RINEX读取、轨道插值、TEC计算、坐标转换、周跳修正等模块、3个Shell脚本含自动化处理流程、2个参数配置列表param.list、filename.list、1个可执行程序rnx2gtex及Makefile、README、头文件Define.h等完整构建组件包体仅290KB结构紧凑、依赖精简。已有155人学习下载读者可直接编译运行获得从原始RINEX观测到GTEX格式TEC序列的端到端处理能力并通过源码深入理解TEC反演中的几何映射、卫星钟差修正、双频组合算法及Julian日历转换等关键技术实现细节。1. RNX2GTEX到底是什么一个老工具背后的硬需求做GNSS数据处理的人大概率绕不开一个词TEC。搞电离层研究、做单频定位误差修正、评估空间天气都需要从GNSS观测数据里提取总电子含量。而RNX2GTEX这个Fortran源码项目干的就是这件事——把RINEX格式的GNSS观测文件转换成电离层TEC值。先说清楚这个工具在整条数据处理链路里的位置。GNSS接收机采集到的原始数据通常是厂商私有格式比如Trimble的.dat、Leica的.m00经过转换后变成标准交换格式RINEX。RINEX文件里记录的是伪距、载波相位、多普勒等观测值但这些观测值本身并不是TEC。TEC需要从L1和L2两个频率的观测值差异中推算出来RNX2GTEX就是在这个环节工作的输入RINEX观测文件输出每条信号路径上的TEC时间序列供后续电离层建模、层析成像或精度评估使用。这个工具适合谁三类人最需要。第一类是电离层物理方向的研究生手里攒了一堆RINEX数据但还没有系统的TEC提取工具链第二类是卫星导航定位领域的工程师需要对观测数据做预处理分析比如评估一条基线上的电离层活跃程度第三类是刚接触Fortran科学计算的开发者想找一个结构清晰、能直接编译运行的GNSS处理源码作为学习范本。无论哪类这个项目都能提供一个完整的、可运行的参考实现。我最初拿到这个源码时有点意外因为现在Python做GNSS处理的生态已经很成熟比如georinex解析RINEX、numpy做矩阵运算为什么还要用Fortran带着这个疑问把源码过了一遍之后我意识到这个选择有它的道理。Fortran在科学计算领域积累了几十年的代码资产尤其是电离层、对流层这类地球物理模型中大量成熟的算法库都是Fortran写的RNX2GTEX作为其中一环用Fortran写可以非常方便地嵌入到既有模型流程中。另外处理长时间序列的RINEX数据时Fortran的数组操作和I/O效率确实比脚本语言高出一截尤其是在没有向量化优化习惯的Python代码面前差距很明显。2. TEC计算原理代码背后的物理与数学2.1 双频观测值如何反演TEC要读懂RNX2GTEX的源码先得理解TEC反演的基本原理。GNSS信号从卫星到接收机的传播路径上会穿过电离层而电离层中的自由电子会对无线电信号产生折射效应。这个效应的核心特征是折射量与信号频率的平方成反比。所以用两个不同频率的信号同时传播两者之间的延迟差就直接反映了传播路径上的总电子含量。具体到数学表达伪距观测方程可以写成P1 ρ c(dt_r - dt_s) I1 T ε1 P2 ρ c(dt_r - dt_s) I2 T ε2其中电离层延迟I与频率的关系是I1 40.3 × TEC / f1² I2 40.3 × TEC / f2²两式相减消去几何距离ρ、钟差和对流层延迟T得到P2 - P1 40.3 × TEC × (1/f2² - 1/f1²)整理后就是代码里最核心的公式TEC (P2 - P1) × f1² × f2² / (40.3 × (f1² - f2²))对于GPS系统f1 1575.42 MHzf2 1227.60 MHz代入后系数大约是9.52也就是TEC ≈ 9.52 × (P2 - P1)。这里的TEC单位是TECU1 TECU 10¹⁶ electrons/m²。RNX2GTEX源码里必然有一段代码在做这个换算但实际工程实现远比公式复杂。因为伪距观测值噪声很大单历元的P2 - P1可能包含几TECU的噪声直接使用效果很差。所以常见做法是用载波相位观测值来平滑伪距或者直接用相位观测值计算相对TEC变化再用伪距TEC来消除相位模糊度。2.2 载波相位TEC与模糊度处理载波相位观测值的噪声远小于伪距大约只有毫米级对应的TEC噪声可以低到0.01 TECU以下。但相位观测值有一个整数模糊度问题导致它只能给出TEC的相对变化无法给出绝对值。相位观测方程类似L1 ρ c(dt_r - dt_s) - I1 T λ1 × N1 ε1 L2 ρ c(dt_r - dt_s) - I2 T λ2 × N2 ε2同样做差得到L1λ1 - L2λ2 -(I1 - I2) λ1N1 - λ2N2整理后TEC (L1λ1 - L2λ2 λ2N2 - λ1N1) × f1² × f2² / (40.3 × (f1² - f2²))这里的λ1N1 - λ2N2是一个常数偏差只要没有周跳它在一段连续观测弧段内保持不变。所以RNX2GTEX这类工具的标准做法是用伪距TEC确定这个常数偏差的初值然后用相位TEC去跟踪高精度的相对变化。这种组合方式兼顾了伪距的绝对性和相位的精密性。源码里实现这个策略时最需要小心的就是周跳检测。一旦发生周跳相位观测值会跳变常数偏差也随之改变如果继续沿用之前的偏差值后续所有TEC估计都会整体偏移。我看到的实现里通常会使用电离层残差组合或者M-W组合来检测周跳碰到周跳就重新初始化偏差。2.3 硬件延迟偏差DCB和斜向TEC转垂直TEC还有一个不可忽视的问题卫星和接收机硬件对两个频率信号的处理延迟不同这种差分硬件延迟(DCBDifferential Code Bias)会直接叠加在TEC估计值上。如果不做修正TEC结果可能偏差几TECU到十几TECU在电离层平静时期这个误差已经相当可观了。处理DCB通常有两条路线。一条是用外部产品修正比如CODE分析中心发布的卫星和接收机DCB文件直接扣除另一条是在数据自身内部估计假设一个区域内所有接收机和卫星的DCB之和在一定时间尺度内恒定通过最小二乘平差同时解算TEC和DCB。RNX2GTEX这类单站处理工具如果源码里没有集成外部DCB文件的接口那它的输出就是非修正TEC也就是包含DCB偏差的原始TEC值。这一点在用的时候要特别注意否则后续建模会引入系统性误差。另外从观测值中提取的TEC是沿卫星到接收机视线方向的斜向TECSTEC而大多数应用场景需要的是垂直TECVTEC。转换需要知道电离层穿刺点IPP的位置和卫星高度角常用的是单层映射函数VTEC STEC × cos(arcsin((Re×sin(90°-el)) / (Reh)))其中Re是地球半径约为6371 kmh是电离层单层高度通常取350-450 km。RNX2GTEX如果输出的是STEC那就需要在后处理中自己加这一层映射如果源码里已经包含了映射函数那它输出的就是VTEC。我建议拿到源码后第一件事就是确认这一点看一下主程序里是否调用了高度角计算和映射函数相关的子程序。3. 源码结构与核心模块拆解3.1 RINEX文件解析模块格式细节决定成败RINEX格式是GNSS数据交换的通用语言但解析它并不轻松。RNX2GTEX的源码首先要解决的就是读取RINEX观测文件通常是RINEX 2.11或2.12版本也可能支持3.x。这个模块的代码量往往占整个项目的三分之一以上属于典型的看着简单写起来琐碎的部分。RINEX观测文件的结构分两部分头部区Header和数据区Body。头部区每一行都有固定的列位置比如第1-60列是描述信息第61-80列是标签。程序需要从头部读取的关键信息包括接收机近似位置APPROX POSITION XYZ、观测类型列表# / TYPES OF OBSERV、观测历元间隔INTERVAL、数据类型标志GGPSRGLONASSEGalileoCBDS等。这里有一个特别容易踩的坑不同版本的RINEX对观测类型的命名规则不同。RINEX 2.x用的是两位字符比如L1、L2、P1、P2、C1、D1、S1RINEX 3.x变成了三位字符比如L1C、L2W、C1C、C2W。RNX2GTEX源码里解析模块通常会维护一个映射表将不同版本的观测码统一映射到内部定义的编号然后根据编号索引对应的观测值数组。如果你自己改源码新增了某种观测类型而没有同步更新映射表程序会大概率在读取数据时直接报错。解析数据区时需要注意历元行的格式。RINEX 2.x的历元行格式为YY MM DD HH MM SS 标志 卫星数然后是每颗卫星的观测值。观测值按头部声明的类型顺序排列缺失值用0.0或空白填充。源码中常见的做法是逐行读取后先判断是否为历元行通过正则或格式检查再按卫星数循环读取观测值每一颗卫星的观测值根据被跳过的空白字段对齐。I/O性能方面RNX2GTEX如果处理的是多天、多站的RINEX数据文件I/O的优化就很重要。我看到一些实现里用了Fortran的直接访问direct access或流式I/O来替代默认的顺序读实测在大文件上能快不少。3.2 观测值预处理与质量控制RINEX数据不等于干净数据。由于多路径效应、接收机异常、卫星故障等原因原始观测值里混入了不少粗差RNX2GTEX的预处理模块就是把这些异常值挡在TEC计算之前。预处理通常包含几个环节。首先是卫星高度角计算利用广播星历或精密星历计算每颗卫星在接收机处的仰角和方位角。高度角低于阈值比如10度的卫星观测值会被剔除因为低高度角信号穿过电离层的路径长多路径效应严重TEC提取精度很难保证。其次是粗差剔除。常用的方法是考察相邻历元之间的观测值变化率如果某个历元的P2-P1值与前一个历元的差值超过设定阈值比如50 TECU就标记为异常并剔除。更严格的做法是结合载波相位平滑伪距后的残差来判断残差超过3倍标准差就剔除。再就是周跳检测。我上面提到过相位观测值的连续性直接决定了TEC跟踪的精度所以RNX2GTEX里大概率会有一个专门检测周跳的子程序。经典的检测手段是使用电离层残差组合LI组合即L1 - L2的线性组合这个组合对周跳非常敏感——一个周跳在LI组合中会产生大约5.4 cm的跳变对GPS远超噪声水平。代码中通常会给LI组合设置一个阈值比如10 cm超过就判定为周跳并触发模糊度重新初始化。3.3 TEC计算与输出格式预处理完成后TEC计算就是标准的公式推演了。RNX2GTEX主体循环伪代码思路通常是打开RINEX观测文件 读取头部信息 循环读取每个历元 循环读取每颗卫星 检查观测值有效性 计算卫星高度角/方位角 如果高度角低于阈值跳过 计算伪距TECP2-P1 计算相位TECL1λ1-L2λ2 用伪距TEC初始化相位TEC常数偏差 检测周跳若发生则重新初始化 输出STEC或VTEC 关闭文件输出部分一般有两种设计。一种是直接输出文本文件每行包含时间、卫星PRN、高度角、方位角、STEC、VTEC等字段方便用绘图工具展示另一种是直接输出成IONEX格式或其他电离层地图产品需要的中间格式方便与CODE、IGS等机构的产品做对比。看RNX2GTEX这个名字GTEX大概率是GNSS TEC EXtraction的意思输出的应该是一种自定义的TEC交换格式后面可以做转换器接到其他可视化工具上。这里我建议拿到代码后先看输出样例确认各个字段的单位和坐标系。比如STEC是正的还是负的相位TEC的符号约定是什么时间是GPS时还是UTC这些细节如果不核对后面做时间序列对比时很容易出问题。4. 编译与运行实战4.1 Fortran环境准备gfortran安装与验证RNX2GTEX是Fortran源码第一步是准备一个可用的Fortran编译器。如果你在Linux或macOS上工作gfortran是最稳妥的选择它是GNU编译器套件的一部分免费、跨平台、兼容性好。大部分Linux发行版都可以通过包管理器直接安装Ubuntu/Debian系用aptCentOS/RHEL系用yum或dnf。macOS上如果装了Homebrewbrew install gcc就会一并安装gfortran。Windows上可以安装MinGW-w64或者使用WSLWindows Subsystem for Linux我个人的建议是直接用WSL跑Linux环境省去一堆环境变量和依赖的麻烦。安装完成后用下面的命令验证编译器是否正常gfortran --version如果输出了版本信息就说明环境就绪。接着写一个hello world验证编译链路program hello print *, Fortran is ready end program hellogfortran -o hello hello.f90 ./hello这里有个小经验老版本的RNX2GTEX源码可能用的是Fortran 77语法.f后缀而新版本可能是Fortran 90/95.f90后缀。gfortran对两者都能编译但处理.f文件时会按F77标准编译如果你的代码里用了F90的自由格式写法记得把文件名改成.f90或者在编译时加-ffree-form选项强制指定自由格式。4.2 源码编译流程从Makefile到手动编译拿到RNX2GTEX源码后先看看项目里有没有Makefile。有Makefile的话直接执行make就能出可执行文件。但很多学术代码的Makefile年代久远编译器选项和现代环境不完全兼容需要手动调整。我的建议是先把Makefile里的FC变量指到gfortran把FFLAGS里的优化选项设置好FC gfortran FFLAGS -O2 -ffixed-line-length-none -fno-range-check LDFLAGS 这里-ffixed-line-length-none很重要因为老Fortran代码默认固定格式下每行不能超过72个字符如果源码里有长表达式不加这个选项会直接编译报错。-fno-range-check是应对某些老代码里整数溢出或类型转换检查过严的问题不过用的时候要谨慎它可能掩盖真实的代码bug。如果项目里没有Makefile也可以手动编译gfortran -O2 -ffixed-line-length-none -o rnx2gtex main.f90 tec_calc.f90 rinex_read.f90注意编译顺序依赖模块一定要在主程序之前编译。如果源码用到了module.mod文件还需要保证module文件在编译时能被找到必要时用-I参数指定module搜索路径。4.3 运行参数与输入输出文件运行RNX2GTEX前需要准备好输入文件。最基本的输入是一份RINEX观测文件文件名通常类似abcd0010.21o这种格式其中abcd是测站名001是年积日0表示该天内第0个会话.21o表示2021年的观测文件。如果你的源码还依赖导航文件比如需要计算卫星位置或高度角那还需要准备对应的.21n或.21p文件。程序运行方式通常是命令行传参./rnx2gtex abcd0010.21o或者交互式输入文件名。运行后应该会生成一个TEC输出文件我建议第一次运行先用小数据量尝试比如单站1小时的RINEX文件确认输出符合预期后再跑完整数据。输出文件的每一行大概长这样2021 001 12 00 00 G01 45.3 182.6 12.34 11.87 2021 001 12 00 00 G02 38.7 75.2 18.92 17.44分别是时间、卫星PRN、高度角、方位角、STEC、VTEC。这些字段名可能因源码版本不同而有差异但核心信息是一致的。5. 常见问题与排查技巧实录5.1 编译期报错的原因和解决RNX2GTEX这类老代码编译报错几乎是必然的不用慌大部分问题集中在几类。最常见的是语法兼容性报错。比如老代码用了PAUSE语句这在现代Fortran标准里已经被废弃gfortran会报告Error: PAUSE statement is not standard可以直接注释掉或者改成CONTINUE。还有OPEN语句里的STATUSNEW如果文件已经存在会报错改成STATUSREPLACE即可。另一个高频问题是隐式类型冲突。老Fortran代码默认I-N规则变量名以I、J、K、L、M、N开头的是整型其他是实型但有些变量名在现代代码里容易引起歧义。比如NMAX被当作整型如果代码里给它赋了浮点值编译时会报Type mismatch。解决办法是在主程序和子程序里统一加上IMPLICIT NONE把所有变量显式声明虽然改起来麻烦但能一次性把这类问题清干净。还有一种情况是整数溢出。处理RINEX数据时如果用到整型变量存储时间戳或数据字节数在老代码里可能是16位或32位整型数据量大了之后会溢出。gfortran默认整型是32位存储文件大小没问题但如果代码里有从UNIX时间戳转换的逻辑建议把相关变量改成INTEGER(8)。5.2 运行期数据异常的类型和应对编译过了不等于能拿到正确结果。我实际跑这类工具时经常遇到几类运行期异常。第一类是输出TEC值全部为零。这通常是输入的RINEX文件里没有找到代码期望的观测类型。比如源码默认读取P1和P2观测值但你的接收机输出的是C1和L2源码里与观测类型匹配相关的逻辑就会失效导致所有观测值被标记为无效。解决方法是查看源码中解析观测类型那段把实际存在的观测码加入匹配列表。第二类是TEC值出现明显的系统性偏差或跳变。系统性偏差大概率是DCB未修正导致前面分析过这是预期内的。跳变则要重点检查周跳检测是否生效可以对照该时段有没有发生磁暴等空间天气事件如果没有就很有可能是周跳漏检。第三类是程序运行到某个历元直接崩溃或死循环。常见原因是RINEX文件数据区某行格式损坏导致解析逻辑进入错误分支。排查时可以打开调试选项重新编译或者直接在读取模块里加打印语句输出出问题历元的时间再回到原始RINEX文件里检查那一行数据。5.3 结果验证与质量评估拿到TEC结果后一定要做质量检查而不是直接进入后续建模。我的习惯做法是三个步骤。第一步画出TEC时间序列图看趋势是否合理。正常情况下白天TEC值高中午前后达到峰值夜间TEC值低日出前达到谷值。一天内STEC变化范围通常在几TECU到几十TECU之间。如果曲线完全没有日变化规律或者数值明显超出合理范围说明处理链路有问题。第二步与外部产品对比。IGS和CODE等机构发布全球电离层地图GIM可以从它们的网站下载对应时刻的TEC数据在自己的测站位置插值再与RNX2GTEX的输出对比。两者差异在几个TECU以内是正常的如果差异巨大需要仔细检查是否漏掉了DCB修正。第三步检查多颗卫星之间的一致性。同一测站同一时刻不同卫星的STEC因为传播路径不同会有差异但VTEC经过映射函数修正后应该趋于一致。如果某颗卫星的VTEC明显偏离其他卫星大概率是这颗卫星数据本身有问题或者处理参数设置不对。以下整理了一份快速排查表方便复现时对照。现象可能原因排查方法编译报错PAUSE/STOP老语法不适配注释或替换为CONTINUE编译报类型不匹配隐式类型冲突加IMPLICIT NONE显式声明输出全零观测类型不匹配检查RINEX头部的观测类型列表TEC日变化无规律卫星钟差/轨道未处理确认是否输入了导航文件TEC整体偏大DCB未修正引入外部DCB改正TEC出现阶跃周跳漏检调整周跳检测阈值运行崩溃RINEX数据行损坏定位并修复损坏行6. 实际应用场景与扩展思路6.1 电离层监测与研究中的典型用法RNX2GTEX提取出的TEC时间序列在科研和工程中用途很广。我自己最近的项目里用这类工具处理了一个中等纬度测站连续三天的数据研究磁暴期间TEC的扰动特征。方法是先提取STEC然后用滑动窗口计算TEC变化率指数ROTIRate of TEC Index能很清楚地看到磁暴期间电离层闪烁活动增强的信号。具体做法是对每个卫星的连续弧段计算每30秒的TEC变化率再在5分钟窗口内统计标准差就得到了ROTI。用RNX2GTEX的输出配合几行Python脚本就能完成import numpy as np import pandas as pd data pd.read_csv(tec_output.txt, delim_whitespaceTrue, names[year,doy,hh,mm,ss,prn,el,az,stec,vtec]) data[time] pd.to_datetime(data[[year,doy,hh,mm,ss]].astype(str).agg(-.join, axis1), format%Y-%j-%H-%M-%S) data data.sort_values([prn,time]) for prn, grp in data.groupby(prn): grp grp.reset_index(dropTrue) dtec grp[stec].diff() / grp[time].diff().dt.total_seconds() * 30 roti dtec.rolling(10, min_periods1).std() data.loc[grp.index, roti] roti这类分析在空间天气研究和导航增强系统性能评估中非常常见。6.2 从单站到网络扩展方向RNX2GTEX的单站版本做出来后下一步很自然的扩展就是处理一个测站网络的TEC数据做区域电离层地图。常见做法是先把每个站的TEC输出统一为VTEC然后选择单层高度比如350 km计算每个观测值的电离层穿刺点经纬度再用克里金插值或球谐函数拟合生成区域TEC地图。这个过程在Fortran源码层面可以加但更便捷的做法是保留RNX2GTEX作为前端提取工具把输出用Python的scipy和pykrige做插值。两条路线的区别在于如果数据量巨大、要求实时处理就在Fortran侧把插值也实现了减少I/O开销如果只是离线研究脚本语言就够用了没必要改Fortran源码。另外一个常见扩展是把TEC输出与接收机DCB估计结合起来。如果你有同一个区域多个测站同时段的数据可以在TEC提取结果的基础上用最小二乘平差同时估计每个测站的接收机DCB和每颗卫星的卫星DCB。这个内容本身可以写一篇单独的博文这里提一句是因为RNX2GTEX输出的STEC文件是这类平差的直接输入做好格式对接能省不少事。6.3 源码改造时的一些个人经验如果你打算在RNX2GTEX基础上做二次开发我分享几个实操中总结的经验。第一先跑通再改造。拿到源码第一件事不是急着读代码而是先找一小段RINEX数据把原始代码完整编译运行一遍确认输出正常。只有基线跑通了你之后做的任何改动才有对照。第二新增输出字段时注意保持原有输出格式稳定。很多下游处理流程依赖固定列数的输出文件如果你加了列下游脚本就会错位。稳妥的做法是新增一个输出文件或者在文件末尾附加新列而不是改中间列的位置。第三RINEX 3.x的多系统支持值得优先考虑。现在越来越多的接收机输出RINEX 3.x格式包含GPS、GLONASS、Galileo、BDS等多系统数据。如果RNX2GTEX只支持GPS可以照葫芦画瓢扩展支持其他系统。注意不同系统的频率不同比如BDS的B1I是1561.098 MHzB2I是1207.14 MHzGalileo的E1是1575.42 MHzE5a是1176.45 MHz这些频率常数必须对应修改直接套GPS的公式会得到错误结果。第四单元测试断不可少。TEC计算函数是整个工具的核心建议给它写一组单元测试输入已知的P1、P2差值验证输出TEC是否与理论值一致。比如P2-P1 1米时TEC应该约为9.52 TECU几个简单用例就能把核心算法钉死后续改代码不会改坏。7. 写在最后的一点实操体会用RNX2GTEX处理TEC数据这条技术路线给我最大的感受是科学计算领域工具不论新旧能解决问题就是好工具。Fortran这门语言被很多人认为过时了但在GNSS电离层处理这种对数值精度和计算效率要求高的场景它依然有不可替代的价值。RNX2GTEX的价值不仅在于它本身能提取TEC更在于它是一个完整的、结构清晰的样例新手可以通过读懂它理解TEC提取的全流程老手可以在它的基础上做二次开发。我自己的做法是将它作为整个处理链路的源头后面接Python做可视化和统计分析各取所长。最后再提醒一句跑任何GNSS处理工具之前先确认你对输入数据的格式和数据质量有充分的了解这是所有后续分析可信度的基础。本文还有配套的精品资源点击获取