简介这份资源面向电离层物理、空间天气与无线电传播方向的研究人员和学生提供名为 codeTEC_L 的 MATLAB 脚本用于处理与可视化 CODEtec 相关的电离层数据。压缩包内仅含 1 个 .m 文件体积约 1KB属于轻量级工具脚本便于直接嵌入既有分析流程。脚本围绕电子密度与总电子含量展开涉及数据读取与清洗、缺失值和噪声处理、电离层模型参数校准、全球分布二维或三维绘图、时间序列与周期性特征分析以及频谱滤波等信号处理环节同时兼顾大规模数据下的算法优化与代码模块化组织。已有 487 人学习适合希望快速上手电离层数据解析、复现可视化结果并理解 CODEtec 数据组织方式的读者参考也可作为空间天气对通信与导航影响研究的入门脚本。1. 电离层 TEC 数据处理从 codeTEC_L.m 看 CODEtec 工具链的落地方式做 GNSS 电离层监测的同行大概率都遇到过这种局面手头攒了一批观测站的 TEC 序列想快速看全球分布、看单站日变化、看磁暴期间的扰动响应结果卡在数据解析和绘图这两步上。codeTEC_L.m就是冲着这个场景来的——它是 CODEtec 工具链里负责电离层 TEC 数据处理与可视化的 MATLAB 脚本配套的codeTEC_L.zip里打包了脚本本体和运行所需的辅助文件。CODEtec 这套东西的核心思路是把电子密度和总电子含量TEC的编码、解码、建模、成图串成一条流水线codeTEC_L.m是其中偏L 波段 / 单站或局域处理的那一环。它适合两类人一是刚接触电离层 TEC 数据、需要一套能跑通的最小可用脚本的研究生和工程师二是已经有自己的数据源、但想借一套现成的预处理加绘图框架来省掉重复造轮子的人。下面按这东西怎么用 → 参数怎么调 → 坑在哪的顺序拆开讲。2. 把 codeTEC_L.m 跑起来环境、数据接口与首次出图2.1 运行环境与依赖确认codeTEC_L.m是纯 MATLAB 脚本不依赖 Simulink也不强制要求额外工具箱但有几处会用到 MATLAB 的基础绘图和矩阵运算能力建议 R2018b 及以上版本。如果你的数据量到了全球多站、长时间序列的级别Parallel Computing Toolbox 能明显缩短循环耗时但这不是跑通脚本的前提。拿到codeTEC_L.zip之后先解压到一个不含中文和空格的路径下这一点在 Windows 上尤其重要——MATLAB 对含中文的路径处理在部分版本上会出问题属于典型的还没开始就算翻车。# 解压后建议的目录结构 codeTEC_L/ ├── codeTEC_L.m # 主脚本 ├── data/ # 放你的 TEC 观测数据 └── output/ # 脚本输出的图和数据解压完先别急着改代码直接在 MATLAB 里cd到脚本所在目录运行一次看报错信息。首次运行大概率会因为找不到数据文件而中断这是正常的——脚本里读取数据的路径需要你按自己的实际情况改。2.2 数据读取与格式适配CODEtec 体系下的 TEC 数据常见来源有两类一类是 GNSS 双频观测反演出来的 TEC 序列通常存成文本或 MAT 文件另一类是电离层模型输出的格网数据。codeTEC_L.m里读取数据的部分一般集中在脚本开头你需要关注的是文件路径、列分隔符和时间格式这三处。% 读取 TEC 数据示例按实际文件格式调整 dataFile fullfile(data, tec_obs.txt); % 假设文件为逗号分隔第一列是时间年 月 日 时 分后续列是各站 TEC rawData readmatrix(dataFile, Delimiter, ,); % 提取时间列并转换为 MATLAB datetime timeVec datetime(rawData(:,1:5)); % 提取 TEC 数值列 tecMatrix rawData(:, 6:end);这里readmatrix是 R2019a 之后推荐的读法老版本用dlmread或textscan也行。关键参数是Delimiter如果你的文件是空格或制表符分隔要相应改成 或\t。datetime那一步把年月日时分拼成时间戳后面画时间序列图时直接拿它当横轴比手动拼字符串省事得多。如果你的数据本身已经是 MATLAB 的.mat文件那就跳过解析直接load进来但要注意变量名和脚本里后续引用的名字对得上。2.3 预处理去野值、补缺失、平滑真实观测数据里出现负 TEC、突跳、整段缺失是家常便饭。codeTEC_L.m通常会带一段预处理逻辑但默认参数未必适合你的数据。常见的做法是先用中位数绝对偏差MAD做粗筛再用滑动窗口做平滑。% 基于 MAD 的野值剔除 medVal median(tecMatrix, omitnan); madVal mad(tecMatrix, 1); % 超过 3 倍 MAD 的判为野值置为 NaN outlierMask abs(tecMatrix - medVal) 3 * madVal; tecMatrix(outlierMask) NaN; % 滑动窗口平滑窗口长度按采样率调整 winLen 5; tecSmooth movmean(tecMatrix, winLen, omitnan);mad的第二个参数设为 1 表示按中位数归一化这是电离层数据处理里比较稳妥的选择。movmean的窗口长度winLen要根据你的采样间隔来定——如果是 30 秒采样5 点窗口对应 2.5 分钟对 TEC 这种缓变量来说够用如果是 1 秒高频采样窗口可以拉到 30 甚至 60。注意omitnan这个选项不加的话一旦窗口里混进 NaN整段平滑结果都会变成 NaN这是很多人第一次跑会踩的坑。2.4 首次出图与结果核对预处理完就可以出图了。codeTEC_L.m的绘图部分一般会生成两类图单站 TEC 时间序列和如果有经纬度信息区域或全球 TEC 分布图。% 单站 TEC 日变化曲线 figure; plot(timeVec, tecSmooth(:, 1), LineWidth, 1.2); xlabel(UT); ylabel(TEC (TECU)); title(单站 TEC 日变化); grid on; % 多站对比取前 5 个站 figure; plot(timeVec, tecSmooth(:, 1:5), LineWidth, 1); legend(STA01,STA02,STA03,STA04,STA05); xlabel(UT); ylabel(TEC (TECU)); title(多站 TEC 对比);出图之后别只看曲线好不好看要核对三件事一是 TEC 量级是否合理中纬度白天一般在 10 到 50 TECU 之间如果动辄几百多半是单位或反演系数出了问题二是日变化形态是否正常白天高、夜间低是基本规律如果反了或者完全没规律检查时间列有没有解析错三是缺失段是否被正确断开了如果 NaN 被当成 0 画进去曲线会出现不自然的掉底。3. 参数调优与建模扩展让 codeTEC_L.m 适配你的数据3.1 时间分辨率与重采样策略不同来源的 TEC 数据时间分辨率差别很大地面 GNSS 站常见 30 秒低轨卫星掩星可能到 1 秒而一些模型输出是 1 小时。codeTEC_L.m默认按输入数据的原始分辨率处理但如果你要拿它跟其他数据源做对比重采样就绕不开。% 将 30 秒采样重采样为 5 分钟均值 tecResampled retime(... timetable(tecSmooth, RowTimes, timeVec), ... regular, mean, TimeStep, minutes(5));retime是处理时间表timetable的利器regular表示按固定步长输出mean是聚合方式。这里要注意重采样之前最好先把野值处理干净否则均值会被个别极端值拉偏。如果你的数据里有明显的周期性比如与太阳活动相关的 27 天周期重采样步长不要设得太大否则会把周期信号抹平。3.2 电离层建模接口IRI 与经验模型codeTEC_L.m本身是一个数据处理脚本不是完整的电离层物理模型但它通常会预留跟经验模型对接的接口。最常见的做法是调用 IRI国际参考电离层模型输出作为背景场然后拿观测 TEC 去减背景看残差扰动。% 调用 IRI 模型获取背景 TEC需 IRI MATLAB 接口 % 假设 iriTEC 已经通过 IRI 接口算好维度与 tecSmooth 一致 tecResidual tecSmooth - iriTEC; % 绘制残差 figure; plot(timeVec, tecResidual(:, 1), LineWidth, 1); xlabel(UT); ylabel(TEC 残差 (TECU)); title(观测 TEC 减 IRI 背景); grid on;这里的关键是 IRI 接口的输入参数要对齐时间、经纬度、高度范围、太阳活动指数F10.7和地磁指数Ap 或 Kp。F10.7 和 Ap 的取值日期必须跟观测时间匹配用错日期的指数是新手常犯的错误出来的背景场会整体偏移。如果你没有 IRI 接口也可以用多项式拟合或滑动中位数当背景虽然粗糙但对识别大尺度扰动够用。3.3 并行加速多站大数据量的处理当站点数量上百、时间跨度几个月时逐站循环会明显拖慢速度。MATLAB 的parfor可以把循环分配到多核上但有几个限制循环体内不能有跨迭代的依赖不能直接往同一个文件写数据。% 多站并行预处理 numStations size(tecMatrix, 2); tecProcessed nan(size(tecMatrix)); parfor s 1:numStations col tecMatrix(:, s); medVal median(col, omitnan); madVal mad(col, 1); col(abs(col - medVal) 3 * madVal) NaN; tecProcessed(:, s) movmean(col, 5, omitnan); endparfor启动本身有开销如果每个站的数据量很小、站数也不多用普通for反而更快。判断标准很简单单次循环体耗时超过几毫秒、循环次数上百才值得上并行。另外parfor里用的变量要提前声明好大小tecProcessed先nan初始化就是干这个的否则 MATLAB 会在每次迭代里动态扩数组速度反而更慢。3.4 输出与批量成图处理完的数据和图片需要批量导出时建议把文件命名规则固定下来方便后续检索和对比。% 批量导出图片 outDir output; if ~exist(outDir, dir); mkdir(outDir); end for s 1:numStations f figure(Visible, off); plot(timeVec, tecProcessed(:, s)); xlabel(UT); ylabel(TEC (TECU)); title(sprintf(Station %02d, s)); saveas(f, fullfile(outDir, sprintf(sta_%02d.png, s))); close(f); endVisible, off让图在后台生成不弹窗批量跑的时候不会满屏窗口。saveas的格式由扩展名决定要矢量图就写.eps或.pdf要位图就.png。如果图很多建议把figure和close配对写不然内存会越吃越多跑到后面可能直接卡死。4. 避坑与排查codeTEC_L.m 实际使用中的五类问题4.1 现象脚本运行报未定义函数或变量原因通常有两个一是codeTEC_L.zip解压后目录结构被改动辅助函数不在 MATLAB 搜索路径里二是脚本里引用了某个工具箱函数但你的 MATLAB 没装那个工具箱。解决方法是先addpath(genpath(codeTEC_L))把整个目录加进路径再which 函数名确认能不能找到。如果提示某个函数找不到查一下它属于哪个工具箱要么装工具箱要么找替代实现。4.2 现象画出来的 TEC 曲线全是 NaN 或断成碎片这多半是预处理阶段野值判据设得太严。3 * madVal这个阈值对干净数据合适但如果你的数据本身噪声大大量正常值会被误判成野值。解决办法是先画一下原始数据的直方图看看分布范围再决定阈值。也可以把判据从固定倍数改成百分位数比如超过 99.5 百分位的才置 NaN这样对不同噪声水平的数据更鲁棒。4.3 现象时间轴对不上日变化出现在错误的地方检查时间解析那一步。常见错误是把年积日DOY当成月日来解析或者把 UT 和 LT 搞混。如果你的数据是年积日格式用datetime(year, 1, 1) days(doy - 1)来转换别直接往datetime里塞。另外确认一下数据文件里的时间是整点还是带分钟readmatrix读进来的数值列如果被截断时间就会整体偏移。4.4 现象并行跑完结果跟串行不一致parfor里如果用了随机数、或者循环体里修改了共享变量结果就会跟串行不一样。电离层数据处理里常见的是在循环内做随机抽样或加噪声这种操作在parfor里必须用RandStream给每个 worker 设独立种子。更稳妥的做法是把随机操作挪到循环外面先生成好随机序列再传进去。4.5 现象批量出图跑到一半内存爆掉MATLAB 的figure对象不手动关的话会一直占内存循环里每画一张就多一份。解决办法就是前面说的close(f)或者用close all在每轮结束时统一清。如果图特别多还可以考虑把数据先存成 MAT另起一个脚本专门读数据出图把处理和绘图拆成两个阶段内存压力会小很多。5. 进阶技巧用 codeTEC_L.m 做扰动指数与事件分析跑通基本流程之后codeTEC_L.m真正有价值的地方在于它能当扰动分析的起点。电离层 TEC 的短期扰动比如磁暴期间的 TEC 增强或耗空是空间天气监测的核心指标之一而扰动指数如 DIX、ROTI的计算并不复杂关键是把背景场和观测场的对齐做扎实。一个我常用的做法是先用滑动中位数窗口取 1 到 2 小时算背景再用观测减背景得到残差最后对残差取滑动标准差就得到一个简化的扰动指数。% 滑动中位数背景窗口 2 小时按 30 秒采样即 240 点 bgWin 240; tecBackground movmedian(tecSmooth, bgWin, omitnan); % 残差 tecRes tecSmooth - tecBackground; % 滑动标准差作为扰动指数 distWin 60; % 30 分钟窗口 tecIndex movstd(tecRes, distWin, omitnan); % 画扰动指数 figure; plot(timeVec, tecIndex(:, 1)); xlabel(UT); ylabel(扰动指数); title(单站 TEC 扰动指数); grid on;这段代码里bgWin和distWin两个窗口长度是核心参数。背景窗口太短会把真实的扰动当成背景减掉太长背景跟不上电离层的快速变化。经验上背景窗口取 1 到 2 小时比较平衡扰动窗口取 30 到 60 分钟。如果你关注的是磁暴主相期间的快速变化两个窗口都可以缩短到 15 到 30 分钟但代价是噪声也会被放大。另一个实用技巧是把扰动指数跟地磁指数Kp、Dst叠在一起看。把 Kp 按时间插值到 TEC 的时间轴上然后双 y 轴画出来扰动响应的时间滞后一目了然。% 假设 kpTime 和 kpVal 是地磁 Kp 数据 kpInterp interp1(kpTime, kpVal, timeVec, linear, extrap); yyaxis left; plot(timeVec, tecIndex(:, 1)); ylabel(TEC 扰动指数); yyaxis right; plot(timeVec, kpInterp); ylabel(Kp); xlabel(UT); title(TEC 扰动与 Kp 对比); grid on;interp1的extrap选项让时间轴两端超出 Kp 数据范围的部分也能有值不然两头会断掉。双 y 轴图在 MATLAB 里用yyaxis实现注意左右轴的plot要分别调用颜色会自动区分。验证这套流程是否靠谱最简单的办法是找一个已知的磁暴事件比如 Kp 大于 6 的日期拿当天的数据跑一遍看扰动指数是不是在磁暴主相期间明显抬升。如果抬升时间跟 Kp 峰值对得上说明窗口参数选得合理如果扰动指数在平静日也一直很高那就是背景窗口太短把噪声当成了扰动。从那以后我每次拿到一批新数据都强制先跑一遍平静日的扰动指数确认基线水平正常再去看事件日的结果。这个习惯帮我挡掉过好几次因为数据段拼接错误导致的假扰动。希望帮到你。本文还有配套的精品资源点击获取
