简介codeTEC_L.m 是一个面向电离层数据分析的 MATLAB 脚本适用于空间物理、大地测量与无线通信领域的师生和研究人员。脚本围绕 CODEtec 数据源实现了从电离层参数读取、噪声清洗、缺失值插补到电子密度与 TEC 变化规律可视化的完整链路可辅助解答电离层对短波通信、卫星导航信号传播的影响等实际问题。压缩包内共 1 个文件类型为 .m 源码整体仅 1KB轻量易读适合直接导入 MATLAB 运行或作为二次开发基底。脚本内容覆盖数据预处理、国际参考电离层模型对照、二维/三维地图投影绘图、时间序列与频谱分析等常用方法并包含模块化的处理思路与简单注释虽然代码量小却清晰展示了处理电离层观测数据的典型流程便于初学者按步骤理解也方便研究者根据需求替换输入数据或扩展绘图功能。目前已有 487 人学习下载适合想通过少量代码快速入门电离层数据可视化与建模分析、并进一步理解空间天气影响的读者。1. 拿到 CODEtec 的 IONEX 文件codeTEC_L 就是那个替你省事的解析器干电离层数据处理的人十有八九第一次拿到 CODEtec 的 IONEX 文件时都会对着那堆 ASCII 格网发怵。codeTEC_L 就是在这个场景里冒出来的它把 CODE 发布的全球电离层 TEC 格网从下载、解析、插值到出图收敛成几条命令和函数调用。我不打算泛泛讲电离层理论直接按自己走过的流程把 codeTEC_L 是什么、怎么跑通、参数怎么调、哪些地方最容易翻车从头到尾说一遍。刚接触 TEC 数据的定位算法工程师、做空间天气数据处理的同学还有要替单频接收机做电离层改正的测绘从业者应该都能在这里找到能直接抄走的部分。2. 拆开 IONEX 黑匣子CODEtec 格网数据里到底有什么先把数据本身看明白。CODEtec 这个名字业内通常指欧盟定轨中心发布的全球电离层 TEC 格网产品。它把全球按经纬度切成格网每小时或每两小时给出一张垂直总电子含量切片所有切片合进一个 IONEX 格式的文本文件里。codeTEC_L 的核心功能就是读这种文件所以如果你连 IONEX 的构成都还不清楚后面调参数就是瞎猜。2.1 CODEtec 不是一块完整地图而是每小时一张的 VTEC 切片一张典型的 CODEtec 地图覆盖经度 -180 到 180、纬度 -87.5 到 87.5纬度步长 2.5 度经度步长 5 度所以每张图是 71×73 个网格点。这里的数值是垂直总电子含量单位是 TECU1 TECU 等于每平方米 1×10^16 个电子。CODEtec 最终产品通常以一天为单位发布一个文件里包含多张按时间排列的切片快速产品延迟小但精度和稳定度不如最终产品。codeTEC_L 在解析时会把每个历元的地图单独切出来放进一个三维数组里第一维是时间第二维是纬度第三维是经度。这个设计让后续按时间抽取和空间插值都变得很直接。如果你以前用普通文本编辑器打开过 IONEX会看到一大串排列整齐的数字但没有任何一眼能看懂的坐标标识——坐标信息全在文件头里不解析头文件根本分不清哪一行对应哪个纬度。这也是为什么很多人第一次手写 IONEX 解析器会翻车只按固定宽度切了数字却忘了文件头里有一堆说明行占了不同的行宽和列数。codeTEC_L 的做法是先定位固定的起始列再按字段宽度读取并且把不同版本 IONEX 的差异收敛到几个配置项里。2.2 IONEX 头文件里的关键字段不读懂就等着读错数据拿到一个 IONEX 文件我先做的第一件事永远是看头文件而不是直接跑脚本。下面的命令能快速打印前 35 行head -35 codeTEC_L_sample.ION输出里需要重点确认的字段有这么几个LAT1 / LAT2 / LON1 / LON2定义了格网的边界DLAT / DLON是纬度和经度步长# OF MAPS IN FILE告诉你有多少张时间切片EXPONENT是数据缩放指数。很多 IONEX 文件里的 TEC 值不是真实 TEC而是乘以 10 的某个次方后存储的整数实际值要乘上 10^EXPONENT 才算对。codeTEC_L 的解析器会主动读这个指数但你自己写脚本时很容易忘。再往后文件的历元信息以EPOCH开头一行一个时间标签标注的是 UTC 或者 GPS 时间取决于发布方。CODEtec 的老文件通常写 UTC但近些年的自动处理文件有混用可能。codeTEC_L 默认按 UTC 解析同时也允许你通过参数强制换到 GPS 时间系统这一步很重要因为时间错一个小时TEC 在高电离层活动期可能差出 5 TECU。对文件完整性做快速检查可以用一个简单的命令grep -c EOF codeTEC_L_sample.IONIONEX 文件里每个历元的数据块以EOF行收尾文件头也有一行EOF。所以这个命令返回的数字应该等于“文件头 EOF 一次 地图数”。如果你算出来的结果比# OF MAPS IN FILE多一个正常少一个或者多很多那就说明文件在传输或解压过程中坏了。这个习惯帮我省过好几次排查时间。2.3 codeTEC_L 在数据链路里的位置解析、插值、输出codeTEC_L 不是一个大而全的电离层建模平台它只负责 IONEX 这条链路上的三件事解析、插值、输出。输入是 CODEtec 的 IONEX 文件输出通常是某个经纬度在某个时刻的 VTEC、某个站星视线方向的 STEC或者一张 TEC 图。它不自己去解算卫星信号也不做层析这些需求得交给其他工具。正因为定位这么窄codeTEC_L 的代码量很小依赖也就 numpy、scipy 和 matplotlib。它把 IONEX 解析器、空间双线性插值、时间线性插值和投影函数分开放在不同模块里你完全可以只借它的解析器插值和投影都换成自己的实现。我见过不少团队就是这么干的把 codeTEC_L 当成一个“IONEX 解码库”后面接自己写的区域电离层模型。在设计层面codeTEC_L 选择用“网格对象 插值方法”的核心抽象解析完文件后得到一个包含时间轴、纬度轴、经度轴和三维 TEC 数组的对象随后所有功能都围绕这个对象展开。这个设计的好处是你不必为了一个单点 TEC 反复重新解析文件一次加载后可以任意取点、取时间、画图、批处理。3. 用 codeTEC_L 跑通第一次 TEC 提取最小命令与三个核心参数现在假设你已经有一份 IONEX 文件codeTEC_L 脚本也放在本地目录里。我用的版本是轻量级 Python 实现Python 3.9 以上就能跑依赖只有 numpy、scipy 和 matplotlib。这个特点很实用因为电离层数据经常在离线环境里处理装一个厚重的依赖树会让现场同事很痛苦。3.1 环境准备和 codeTEC_L 的常见目录结构拿到手先看目录里有什么。常见的 codeTEC_L 目录结构大概是这样的codeTEC_L/ ├── ionex_parser.py # IONEX 文件解析 ├── tec_interp.py # 空间和时间插值 ├── stec.py # 垂直 TEC 到斜 TEC 的投影 ├── cli.py # 命令行入口 └── config.yaml # 默认参数配置这种拆分方式很直观解析器只读文件插值器只算网格投影模块负责几何关系。如果你只想要某个点的 VTEC整个链路只走ionex_parser和tec_interp两步。config.yaml里会存默认的插值方式、默认等效电离层高度、默认时间系统我建议你在跑批处理前先打开看一眼而不是直接信任内置值。配置里的常见参数有mapping_height、interp_order、time_system。不同衍生版本叫法可能略有差异但含义一致。把这几个值统一调好后面所有命令行调用都会稳定很多。3.2 最小命令从一个 IONEX 文件里抓出某个经纬度的 VTEC先跑一个最小命令从 CODEtec 文件里提取某个点的垂直 TECpython cli.py --input CODG2024001.ION --lat 30.5 --lon 116.4 --time 2024-01-01 12:00 --output vtec.txt--input指向 CODEtec 的 IONEX 文件--lat和--lon是你要提取的地理坐标--time是 UTC 时间。命令会在终端打印 VTEC 值同时把结果写到vtec.txt里。如果你更习惯写 Python用接口的方式也差不多from codeTEC_L import CodeTecGrid grid CodeTecGrid(CODG2024001.ION) vtec grid.extract(lat30.5, lon116.4, time2024-01-01 12:00) print(fVTEC at test point: {vtec:.2f} TECU)这里CodeTecGrid是 codeTEC_L 里的核心对象加载文件时就完成了头文件解析和三维 TEC 数组构造。extract方法先找离目标时刻最近的两个时间切片再做时间线性插值然后在每张切片上做空间双线性插值。如果文件内只有一个时刻时间插值会自动跳过只做空间插值这很符合单历元的分析场景。3.3 三个核心参数时间、坐标、插值窗口主要有三个参数决定提取质量时间系统、插值顺序、目标高度。下面这张表是我一般在命令行里最常对照的参数取值示例影响--time-systemutc/gps决定 EPOCH 时间标签的解析方式匹配错误会引入整小时偏差--interp-orderlinear/nearest空间插值方式格网稀疏区域用 nearest 更稳日常分析用 linear--mapping-height350/450等效电离层高度影响斜 TEC 计算和穿刺点位置先说时间系统。CODEtec 文件里的 EPOCH 标签不一定全是 UTC有些快速产品会直接用 GPS 时间。codeTEC_L 的--time-system就是用来切换这个解析基准的。如果你把 UTC 的时间当成 GPS 时间去读出来的结果会整体偏移约 18 秒这个误差对单点 VTEC 影响不大但和载波相位实测对比时肉眼可见。插值顺序也值得注意。CODEtec 的全球格网在低纬度经度方向比较稀疏每 5 度一个点双线性插值在格网内表现没问题但目标点落在格网边缘时容易外插出离谱的值。此时设成nearest更安全。我通常先跑一次默认配置再观察目标点周围四个格网点的 TEC 差异如果差异超过 5 TECU就要警惕附近是不是电离层梯度很大的区域。mapping-height 这个参数在纯 VTEC 提取里不参与计算但一旦你要算斜 TEC 或者穿刺点坐标它就直接决定结果。后一章会展开聊这里先记住默认值别乱改。4. 从垂直 TEC 到斜 TECcodeTEC_L 的投影换算与 TEC 图输出很多用户拿到 VTEC 后直接拿去做单频改正这是不对的。接收机观测到的电离层延迟发生在信号斜穿电离层的路径上而 CODEtec 给的是垂直方向的电子总量。要得到斜路径上的 TEC必须做投影换算。codeTEC_L 的stec.py就是干这件事的。4.1 单频用户为什么不能直接用 VTEC垂直 TEC 与斜 TEC 的关系可以用一个简化的投影函数近似STEC ≈ VTEC / cos(zenith_ipp)。这里的 zenith_ipp 不是测站处的卫星天顶角而是信号穿过等效电离层薄壳时穿刺点位置的天顶角。由于测站和穿刺点并不重合直接拿测站天顶角代入会引入系统偏差尤其是在低高度角时候。等效电离层高度是这套近似里最敏感的参数。欧洲定轨中心的全球格网模型通常默认把电离层压缩在 350km 到 450km 之间的一个薄壳上薄壳高度决定了穿刺点的经纬度也就决定了你去格网上哪几个网格点插值。codeTEC_L 里mapping_height的默认值设在 450km这符合一部分最终产品的参数设定但如果你在低纬度用 350kmSTEC 的差异能到 1 到 3 TECU。所以严格来说VTEC 到 STEC 的换算不是一个固定系数而是一个依赖几何的投影函数。codeTEC_L 的做法是先根据测站坐标、卫星方位角高度角、等效电离层高度算出穿刺点地理坐标然后在穿刺点处插值 VTEC再除以穿刺点处的投影系数。4.2 codeTEC_L 的斜 TEC 接口方位角、高度角与映射高度要算一个测站到一颗卫星的斜 TEC接口调用是from codeTEC_L import CodeTecGrid grid CodeTecGrid(CODG2024001.ION) stec grid.stec(station(30.5, 116.4, 0.0), satellite(156.3, 43.2), # 方位角 156.3°高度角 43.2° time2024-01-01 12:00, mapping_height450.0) print(fSTEC: {stec:.2f} TECU)station元组分别是测站纬度、经度和海拔海拔在大多数情况下可以填 0因为 IONEX 格网是大地坐标系统下的 VTEC测站海拔对最终结果的影响远小于映射高度的影响。satellite里第一个数是方位角第二个是高度角单位是度。codeTEC_L 内部先按 450km 高度计算 IPP 位置然后在 IPP 处用双线性插值得到 VTEC最后除以 IPP 天顶角余弦。这里要注意如果卫星高度角低于 10 度投影函数会被分母的余弦压得非常大任何 VTEC 误差都会被放大所以实际使用时我一般只处理高度角大于 15 度到 20 度的观测。如果同一时刻有多颗卫星正确做法是循环调用这个接口而不是自己改satellite参数重算整张图。因为整个文件解析只做了一次后面的插值和投影都是纯数值操作速度很快。在实测数据里你可以把测站到星历计算的方位角、高度角每 30 秒更新一次然后逐历元调用stec就能生成一段时间内的斜 TEC 序列。4.3 画一张带格网的 TEC 图用 matplotlib 快速出图除了单点提取codeTEC_L 另一个常用功能是出图。把网格对象里的 TEC 数组直接喂给 matplotlib 的contourf就能出图不需要引入地理底图包import matplotlib.pyplot as plt lat grid.lat_grid lon grid.lon_grid vtec_map grid.vtec_maps[0] # 第 1 个历元的 VTEC 格网 plt.figure(figsize(10, 5)) cs plt.contourf(lon, lat, vtec_map, cmapjet, levels15) plt.colorbar(cs, labelVTEC [TECU]) plt.xlabel(Longitude [deg]) plt.ylabel(Latitude [deg]) plt.title(CODEtec VTEC map at epoch #1) plt.savefig(codeTEC_L_map.png, dpi150)先检查grid.vtec_maps[0]的形状是不是 (71, 73)再传给contourf。如果你的文件不是标准 2.5×5 度格网lat和lon数组长度会变化但绘图代码不需要改因为坐标都是从对象里取的。出图时最容易忽略的坑是经纬度轴的顺序。有的 IONEX 文件按纬度从南到北排列有的按从北到南排列如果画出来图是上下颠倒的不是投影问题而是解析器没有正确读取LAT2 LAT1还是反向。codeTEC_L 在解析时统一把纬度转成从小到大但这个逻辑只在解析器里生效你自己手写读取数据时一定要处理。5. 处理 CODEtec 时最容易翻车的 5 个问题现象、原因、解法这一章记录我在实际使用 codeTEC_L 过程中踩过的几个高频坑每一条都按现象的套路来说方便你对照排查。5.1 历元对齐出错UTC 和 GPS 时间混用导致整小时偏差现象是提取出来的 VTEC 和实测观测始终差一两个 TECU画时间序列时还能看到固定的跳变点。原因大概率是 IONEX 文件里的 EPOCH 标签被 codeTEC_L 按 UTC 解析但文件实际用的是 GPS 时间或者反过来。CODEtec 老版本产品多以 UTC 为准但部分快速文件用 GPS不同年份的文件混着看很容易翻车。解决方法是显式指定时间系统。在命令行里加--time-system gps或--time-system utc在读入文件后打印第一张 map 的 EPOCH 时间和文件头里写的时间核对一次。只要有一次对齐错误后续所有统计都会被污染所以我在批处理脚本里永远会显式写这个参数不依赖默认值。5.2 高纬度和跨日期边界出现 NaN插值窗口没兜底现象是站点在北极圈或者分析时段正好跨越午夜 0 点输出结果里有大量nan。原因是 IONEX 文件虽然覆盖全球但在高纬度投影到站星视线方向时穿刺点可能落到格网纬度边界之外跨日期时部分 codeTEC_L 版本不会自动加载前一天或后一天的文件时间插值缺一侧数据。解决方法是两层一层是在空间插值时对越界坐标做 clamp把坐标限制在格网边界内用最近邻值兜底另一层是在时间插值前检查文件的首末历元覆盖范围如果目标时间在范围外主动在批处理脚本里拼接相邻日期的 IONEX 文件。我在 codeTEC_L 的调用层加了一个很小的预处理函数专门负责把输入文件列表按时间排序和拼接这之后跨日期的 NaN 基本绝迹。5.3 不同文件版本里经纬度步长不一致默认参数失效现象是某一天读文件正常换一批文件后程序报维度不匹配或者坐标全乱。原因是 CODEtec 在不同年份调整过格网分辨率有的文件是 2.5°×5°有的可能是更高分辨率或不同步长如果 codeTEC_L 硬编码了 71×73 格网就会出问题。解决方法是让解析器永远从DLAT、DLON和LAT1等头字段动态生成网格轴而不是写死。我在配置里关闭了“固定格网尺寸”的开关后再没遇到过这类问题。这也是第一条我建议你拿到新版 codeTEC_L 后先确认的代码逻辑。5.4 穿刺点高度选了 350 还是 450结果差多少现象是同一颗卫星、同一个时刻用 350km 和 450km 算出的 STEC 能差出 1 到 3 TECU。原因不是 codeTEC_L 算错而是电离层薄壳模型本身对高度假设很敏感高度变化改变了穿刺点坐标从而改变了插值位置和投影系数。解决办法不是找到一个“正确高度”而是统一口径。短基线单频定位用 350km 或 450km 对最终定位结果影响有限但你做数据对比时必须保证所有时期的处理都用同一个mapping_height。我在处理分析里会把这个值写进输出文件的属性里方便以后追溯。5.5 文件命名与压缩格式误判codeTEC_L 读不进压缩盘现象是从服务器下载回来的文件后缀是.Z直接扔给 codeTEC_L 报解压错误。原因很简单IONEX 发布方常用 Unix 压缩而 codeTEC_L 默认只处理明文.ION和部分.gz。解决方法是不要改文件后缀而是在预处理里判断扩展名对.Z先调系统解压命令或者在 Python 里按流式方式解压到内存再交给解析器。我一般在批处理脚本里写一个小函数把.Z、.gz、.ion统一成明文路径再传给CodeTecGrid这样就不会因为压缩格式打断整条处理链。6. 用 codeTEC_L 批量处理多年 IONEX从脚本到精度验证到这一步你已经能用 codeTEC_L 处理单文件和单站数据。再往前的方向是批处理以及用实测数据验证 TEC 精度这里说一个我实际在用的最小方案。6.1 多文件批处理构造时间序列并自动重命名输出批处理的核心是让文件列表和时间轴对齐。我通常这样写import glob from codeTEC_L import CodeTecGrid files sorted(glob.glob(COD*.ION)) for f in files: grid CodeTecGrid(f) vtec grid.extract(lat30.5, lon116.4, time12:00) print(f{f}: {vtec:.2f} TECU)这段代码会把所有匹配COD*.ION的文件都读一遍然后提取正午时刻的 VTEC。这里有个容易忽略的地方time12:00只指定了当天时间具体日期取自文件名所以你需要确保文件名里有日期信息比如CODG2024001.ION里的2024001是年积日。如果文件名不标准就改成从文件头 EPOCH 里读取日期。批量处理时建议每读一个文件就打印一行日志不要等到全部跑完再看结果。这样能在中途快速发现某个文件损坏或者文件版本不一致的问题。6.2 和实测载波相位 TEC 对比算 RMS 时要注意边缘效应用双频载波相位观测可以算出信号路径上的高精度 STEC把它和 codeTEC_L 的 STEC 相减就能评估模型偏差。计算公式很简单residual stec_measured - stec_model但在统计 RMS 之前我强烈建议先设置高度角掩码。低高度角的斜路径穿过电离层很长投影函数和模型误差都会被放大残差会明显偏大把这些数据放进 RMS 里会拉高整体指标掩盖真实精度水平。我一般只统计高度角大于 30 度的观测残差并且把残差超过 10 TECU 的点标记为异常值单独检查。我现在拿到新的 IONEX 文件第一件事不是画图而是先打印文件头里的历元数和经纬度范围。这个习惯帮我躲过了很多次文件版本不一致的坑。希望帮到你。本文还有配套的精品资源点击获取