简介splitr是一套面向大气科学研究者与R用户的HYSPLIT轨迹与扩散建模工具包支持从指定接收位置反向追溯气团来源也能对气相或颗粒物的正向/反向输送过程进行模拟广泛应用于极端降水溯源、远距离粉尘与生物气溶胶关联分析、入侵物种扩散及植物病原体长距离传播等研究场景。资源包共108个文件压缩后约10.9MB主体为38个R脚本和26个Rd帮助文档同时附有若干Windows可执行程序、测试数据及示例图便于直接运行与二次开发。已有581人学习下载。通过完整包目录可快速掌握建模流程省去自行配置HYSPLIT编译环境的麻烦结合自带示例与文档可对照理解参数设置、输出解析及可视化的关键步骤适合从事大气环境、生态保护或农业病害预警的科研人员与中级以上R用户使用。1. 把污染溯源放进 R 的工作流里splitr 到底解决了什么拿到一组污染监测数据浓度超标的时段很明确但来源说不清这是环境数据分析里最常遇到的场景。过去做后向轨迹溯源常规流程是打开网页端或桌面端工具填表单、选日期、等结果下载再把轨迹文件导入其他软件绘图。步骤多中间任何一步口径不一致结果就要推翻重来。splitr 这个 R 包做的事情是把 HYSPLIT 模型的轨迹计算拉进 R 环境内部让数据分析师用一套代码完成从气象数据读取、轨迹计算到聚类分析的全流程。我最初接触 splitr 时也带着疑惑HYSPLIT 本身已经足够成熟为什么还要在 R 里再包一层实际用下来才体会到它的价值不在“算轨迹”本身而在“可复现”——每次分析的时间范围、起始高度、气象场、聚类参数全部固化成代码换一台机器重新跑一遍结果一模一样。对于要做污染溯源、气团传输路径分析或区域输送贡献评估的从业者来说这比反复操作图形界面靠谱得多。本文直接从实际使用角度拆解 splitr 的完整流程从安装到聚类再到你需要提前知道的坑。2. 从 HYSPLIT 到 splitr轨迹计算的分工与最小可跑通案例2.1 HYSPLIT 是什么splitr 在其中扮演什么角色HYSPLITHybrid Single-Particle Lagrangian Integrated Trajectory是用于计算气团运动轨迹和扩散过程的模型系统由研究机构开发并持续维护。它的核心是拉格朗日方法把一个虚拟粒子放置在起始位置然后根据气象场的风场分量逐步推算出它在过去或未来一段时间内的位置序列。后向轨迹backward trajectory是最常见的应用即从受体点出发倒推过去 72 小时或更长时间内气团来过哪里从而辅助判断污染物可能的上游来源区域。splitr 则是在 R 层面封装了与 HYSPLIT 的交互逻辑。它负责几件事拼接输入参数、调用轨迹计算核心、读取返回的轨迹文件、并把结果整理成 R 里的数据框结构。可以这样理解HYSPLIT 是计算引擎splitr 是驱动引擎的遥控器而“对其进行更多处理”指的是 splitr 自带的聚类函数、浓度格点分析函数和绘图函数——这些在传统流程里通常需要依赖第三方工具完成如今全部在 R 里收尾。我一般会把整个过程拆成四段准备气象数据、计算轨迹、后处理分析、可视化输出。splitr 覆盖前三段的大部分工作可视化则留给 R 生态里的绘图包继续扩展。2.2 安装 splitr 与跑通第一个后向轨迹安装 splitr 前需要确认 R 版本在 3.5 以上并且系统里有可用的 RtoolsWindows 场景下编译包时需要。如果从 CRAN 安装直接执行标准安装命令如果安装的是开发版本则从代码仓库拉取。注意一点splitr 依赖的某些包可能同时被其他环境占用版本建议在干净环境里安装避免依赖冲突。# 安装 splitr 的常见方式二选一 install.packages(splitr) # 从 CRAN 安装稳定版 # remotes::install_github(作者名/splitr) # 从开发仓库安装最新版需自行确认可用性 # 加载包 library(splitr) # 第一阶段下载或指定气象数据 met_data - C:/metdata # 本地存放气象文件的目录目录内需有对应的数据文件逻辑说明第一步安装包第二步加载库第三步设置气象数据目录。splitr 不会自动连网下载气象文件它要求本地有一个数据目录里面存放 HYSPLIT 格式的气象数据文件。路径中不要出现中文和空格这是后面最容易出错的地方。气象数据文件通常在官方网站注册后获取也可以借助相关工具自动下载归档数据但本文以本地目录假设为例。参数说明met_data 指向的是目录而不是单个文件。splitr 内部通过日期和时长参数去匹配目录下的具体文件因此目录里的数据必须覆盖你要模拟的时间段否则计算会出错。这一点越早确认越好我最初在这里栽过跟头配好了参数却一直报错最后发现目录里少了当天的数据文件。跑通第一个轨迹的基本代码如下# 定义后向轨迹模拟参数 traj - trajec( met met_data, lat 39.9, # 受体点纬度例如某城市监测站点 lon 116.4, # 受体点经度 height 100, # 起始高度单位米 duration 72, # 后向模拟时长单位小时 daily_hours c(0, 6, 12, 18), # 每天生成轨迹的时刻UTC时间 start_date 2023-01-01, end_date 2023-01-05, direction backward ) # 查看轨迹数据结构 head(traj)代码逻辑trajec() 是 splitr 中最核心的轨迹计算函数。传入气象数据目录、经纬度、起始高度、模拟时长、起止日期和方向包会拼接出 HYSPLIT 需要的控制文件调用底层模型最后返回一个数据框。数据框每一行是轨迹上的一个点包含时间、经纬度、高度等字段。参数说明daily_hours 控制每天从哪几个整点时刻起始计算轨迹。c(0,6,12,18) 意味着每天生成 4 条轨迹分别在 UTC 时间的 0 点、6 点、12 点和 18 点。如果你做的是污染过程分析通常只需要选 1 到 2 个时刻如果做季节尺度的来源统计多时刻设置更稳健。duration 是后向追溯的小时数72 小时3 天是污染溯源常用的配置但具体要看污染事件的持续时间和气象场覆盖范围。height 为 100 米是指轨迹起始点在地面以上的高度一般选择 100 米、500 米或 1000 米来代表近地面、边界层中层和高层。运行成功后traj 对象里已经有了全部轨迹的逐小时位置数据。下一步不是急着画图而是先确认数据完整性检查每个起点时刻对应的轨迹长度是否一致、高度是否出现极大异常值。这些细节会在后续章节展开。3. 轨迹计算的四个关键参数时长、高度、日期窗口和气象数据少一个结果就偏3.1 模拟时长与起始高度不同污染物来源的场景差异后向轨迹的时长设置直接影响结果解读的方向。72 小时是气团输送研究的典型值但沙尘过程经常需要追溯到 120 小时以上因为沙尘源区可能远在数千公里之外而分析局地二次污染时24 到 48 小时就足够时间拉太长反而会引入无关区域的干扰。起始高度同样有讲究。近地面监测点测到的是混合充分的近地层空气100 米高度轨迹代表的是近地面气团走势如果想了解边界层中部的输送路径可以用 500 米分析高空传输、例如沙尘的长距离迁移时1000 米以上更合适。我见过不少分析报告把三个高度混在一起画图结果轨迹乱成一团每条线代表的高度层都不一样读图时很难归因。常见做法是固定一个高度层或者分别计算后分层绘图。# 同一日期、同一受体点、不同起始高度的轨迹计算 heights - c(100, 500, 1000) traj_list - lapply(heights, function(h) { trajec( met met_data, lat 39.9, lon 116.4, height h, duration 72, daily_hours c(6), # 本次只取每天 6 点UTC一条 start_date 2023-01-01, end_date 2023-01-03, direction backward ) }) # 合并后添加高度标识 library(dplyr) traj_all - bind_rows(traj_list, .id height_id)逻辑说明lapply 循环依次对三个高度调用轨迹计算最终结果按高度分组合并。这里 daily_hours 只设置一个时刻是为了在控制变量时减少数据量便于逐条检查每个高度对应的路径差异。参数说明daily_hours 写成 c(6) 而不是 6是因为函数期望接收一个向量。如果你在代码里看到 c(0,6,12,18) 这种写法含义是同一批模拟里每天生成多个起始时刻的轨迹。多高度结果合并后画图时应按 height_id 分组着色并用 facet 或分面展示避免线条重叠导致无法判断高度层贡献。3.2 气象数据文件格式、时间覆盖范围和本地目录组织HYSPLIT 计算依赖的气象场数据有两种常见格式ARL 打包格式和 NetCDF 格式。splitr 对两种格式都有支持但实际使用中 ARL 格式的历史数据更全、兼容性更稳定。如果只有 NetCDF 格式的气象场需要确认文件里的变量命名与包内部解析逻辑一致否则会出现读取出 NaN 的情况。气象数据的组织方式是另一个高频问题点。我第一次使用时按年份分目录存放每个月一个子目录结果 trajec() 在跨月模拟时找不到下个月的文件。原因在于包内部按日期顺序拼接文件跨月场景下需要对应的是同目录下的连续文件而不是按月份切换路径。常见做法是建立一个统一的数据仓库目录内部将文件按月或按日命名并保证目录下所有文件的时间连续。例如一个名为“202301”的目录中存放 2023 年 1 月整月的气象文件1 月 30 日到 2 月 1 日的跨月模拟需要同时存在 1 月和 2 月的目录这时可以创建符号链接或软连接来让包在同一个逻辑路径下找到全部文件。我个人的习惯是建一个“all_met”目录用脚本把需要的历史数据全部拷贝或链接进去这样轨迹模拟代码里只需要维护一个路径变量日期窗口随意调整文件查找逻辑稳定可靠。# 创建软连接Linux/macOS 示例 ln -sf /data/met/202301/gdas1.jan23.w3 ./ ln -sf /data/met/202302/gdas1.feb23.w3 ./逻辑说明将 1 月和 2 月的数据文件链接到当前工作目录或统一的 met 目录splitr 在按日期匹配时就能在同目录下找到连续文件。Windows 环境下没有 ln -sf需要改用管理员模式下的 mklink或直接把文件拷贝进同一目录。参数说明文件名中的 jan23、feb23 是 HYSPLIT 气象文件名的月份缩写习惯具体命名取决于你的数据来源渠道。如果文件命名方式跟包默认的解析规则不匹配可以在 trajec() 调用前后检查返回对象的行数为空即说明文件匹配失败。4. 把“更多处理”落到实处轨迹聚类、浓度格点分析与结果解读4.1 轨迹聚类的运行逻辑从几十条路径里提炼主输送通道跑完一个月的每日轨迹得到 30 条甚至 60 条曲线后直接画在一张图上会变成一团乱麻每条轨迹彼此交叉肉眼根本看不出主要输送路径。聚类分析在这里的作用是把空间形态相似的轨迹归为一类每一类代表一条典型的气流路径类别之间的差异则表征不同的来源方向或输送距离。splitr 中实现轨迹聚类的常见流程是先调用聚类函数设定目标类别数内部通过距离算法反复迭代直到簇内距离不再明显下降。聚类结果会输出每个类别包含的轨迹数量、该类别的平均路径以及每条原始轨迹被归到了哪个类。# 假设已经通过 trajec() 得到多日轨迹 traj_all执行聚类分析 clust - trajectory_cluster( traj traj_all, method ward, # 聚类方法 n_clusters 4 # 目标类别数 ) # 查看聚类结果概要 summary(clust)逻辑说明trajectory_cluster() 接受前面计算出的轨迹数据框按空间位置相似度进行层次聚类。method 参数选 ward 时聚类倾向生成规模相近的类别如果气流方向差异特别大方法可以换为 average对异常轨迹的敏感度不同。n_clusters 是典型的需要人工调试的参数不能一次性定死。参数说明n_clusters 的取值没有标准答案。常规做法是设定一个从 2 到 8 的范围逐个计算并画出一条“类别数与簇内距离”的曲线选择曲线拐点处的类别数作为最终值。这个方法在文献里叫肘部法则elbow method虽然听起来有点像玄学但在轨迹聚类场景下比拍脑袋定数字可靠得多。如果某个类的轨迹条数占比不到 5%说明类别数定多了合并到最近的类里重新跑。聚类完成后还需要把类别标签关联回原始轨迹否则后续画图时无法按类着色。# 将聚类标签合并回原始轨迹数据 traj_labeled - left_join(traj_all, clust$cluster_info, by traj_id) table(traj_labeled$cluster)逻辑说明left_join 按轨迹 ID 关联相当于给每一天、每一时刻的轨迹打上类别标签。table() 输出每个类别的轨迹数量用来判断各类占比是否合理。参数说明by traj_id 要求原始数据框和聚类结果里都有这个字段。trajec() 默认返回的数据框包含每条轨迹的唯一 ID但如果你的代码对数据做过筛选或修改这个字段可能丢失。遇到关联失败时检查列名是否完整是最快的排查路径。4.2 浓度格点分析把轨迹频次转化成格点浓度图另一种常见的“更多处理”是浓度格点分析也叫轨迹频次分析。它的思路是把整个研究区域划成规则的经纬度网格统计每一条后向轨迹经过每个网格的次数次数越多说明气团经过该区域的频率越高该区域被推断为潜在源区的概率也越大。splitr 中这部分处理包括两步把轨迹点映射到网格统计每个格点的命中频次。随后将频次矩阵输出再交给绘图函数生成浓度分布图。# 执行网格化频次计算 grid_res - 0.5 # 网格分辨率单位度 conc - concentration_analysis( traj traj_labeled, lat.min 20, lat.max 45, lon.min 100, lon.max 130, grid.res grid_res ) head(conc)逻辑说明concentration_analysis() 统计研究区域内每个网格点的轨迹频次。lat.min、lat.max、lon.min、lon.max 限定研究区范围grid.res 决定网格粗细。0.5 度网格在省市尺度的分析中比较常用如果是全国范围1 度更合适。参数说明网格分辨率不是越小越好。网格太小单条轨迹的经过点很少频次统计会稀疏到看不出规律还大幅增加计算耗时网格太大空间细节被抹平源区的定位价值降低。0.5 到 1 度之间是大多数分析的平衡区间。浓度分析的结果里附带了经纬度网格中心和对应的频次值可以直接用于绘制热力图或等值线。需要注意这里的“浓度”并不是真实污染物浓度而是轨迹密度的代理指标。如果手上同时有站点污染物浓度数据更进阶的做法是对轨迹做加权例如按浓度值加权后再统计频次这样得到的源区分布更贴近实际污染贡献。splitr 的基础函数不直接做加权浓度分析但处理后的数据框允许你手动乘上一个权重列再聚合这一步用 dplyr 几行就能完成。5. 常见避坑清单splitr 跑不通、结果异常先查这五处5.1 气象数据目录路径含中文或空格加载报错现象trajec() 执行后报错提示找不到文件或文件路径无效但手动打开路径确认文件确实存在。原因splitr 在拼接控制文件路径时如果路径包含中文字符或空格底层模型程序解析时会出现意外截断导致路径失效。这个问题在 Windows 系统上尤为频繁。解决把气象数据目录和项目工作目录都改成纯英文、无空格的路径。例如 C:/metdata 而非 C:/我的气象数据。如果是公司内网共享盘可以先映射一个本地盘符再把数据拷贝到本地工作目录不要直接引用 UNC 路径。5.2 轨迹条数比预期少一半原因是日期时间使用了本地时区现象设置了 daily_hours c(0, 6, 12, 18) 和 3 天的模拟日期期望得到 12 条轨迹结果只返回 6 条。原因splitr 内部的时间处理按 UTC 标准如果你传入的日期字符串被系统解释为本地时间例如 UTC8则每天 0 点到 18 点之间的部分时刻会溢出到下一天导致日期范围边缘的轨迹被过滤掉。解决确认 start_date 和 end_date 的时区语义。常见的做法是统一按 UTC 日期传入计算完成后再把输出时间列转换到本地时区用于画图。更稳妥的方式是显式声明时区as.POSIXct(2023-01-01, tz UTC)。5.3 聚类结果中某一类只包含一条轨迹类别数设置不合理现象clust 输出中某个 cluster 的轨迹数量为 1另一类的轨迹数量接近总数一半图形明显失衡。原因n_clusters 设得偏高或者数据里存在极端长距离轨迹与大部分短轨迹之间的距离差过大导致聚类算法把极端轨迹单独切成一个类。解决不盲目信首次运行结果。把 n_clusters 从 2 到 8 各跑一遍画类别数与平均簇内距离的关系曲线选取拐点。若极端轨迹影响明显可以先把超出合理空间范围的轨迹剔除后再聚类。我在实际分析中遇到过沙尘过程的长距离轨迹把聚类结果带偏的情况处理方式是先按轨迹最大长度过滤再看聚类形态。5.4 轨迹高度字段出现 -9999 或 NA绘制的垂直剖面断裂现象轨迹数据框的 height 列出现大量 -9999 或 NA画高度-时间剖面图时线条中断或水平轨迹图上出现飞到上万米的异常线段。原因气象数据文件在某些气压层缺失数据模型在插值时返回填充值。也可能是在数据处理环节中过滤了某一层的数据后没有保留完整字段。解决将填充值统一转为 NA再做逐时段的线性插值。如果异常出现在个别时段直接剔除该时段的重采样点如果大面积出现说明气象文件本身有问题换数据源或检查文件完整性不要在后处理阶段强行补值。5.5 轨迹起点坐标偏移经纬度顺序写反现象打印轨迹前几行时发现起点经纬度和自己传入的参数不一致例如 lat 等于 39.9lon 等于 116.4但数据框里第一行显示 lat 116.4, lon 39.9。原因部分旧版数据处理函数内部按 lon, lat 的顺序组织坐标而 trajec() 的参数顺序是 lat, lon。如果按直觉把经纬度填反起点就会跑到完全不同的区域后续轨迹路径全部偏移数百公里。解决在 trajec() 调用前打印出传入参数确认顺序trajec() 返回后马上检查 head(traj) 第一行的起点坐标是否与设定一致。养成这个习惯后可以避免一大类因坐标顺序导致的返工。这个问题看起来基础我在交付分析报告前自查时抓到过不止一次。6. 进阶技巧把批量回溯与自定义后处理函数结合形成自己的分析闭环当 splitr 的基本流程跑通后值得投入时间做的是批量化。污染溯源分析很少只跑一两天一跑就是一个月甚至跨季节这时适合把整条流程封装成一个可复用脚本输入是日期序列和受体点坐标输出是聚类结果和浓度格点数据。常见的做法是写一个 R 脚本定义三个层次的函数单日轨迹计算函数、整段日期循环函数、结果汇总导出函数。循环部分用 purrr 或基础 lapply 遍历日期序列每个日期生成指定时刻的轨迹最后把全部结果合并。合并后的数据框保存成 CSV 或 RDS 格式方便后续直接读取不用每次都重新调模型。要注意批量处理时的内存增长问题。如果一次运行 90 天的数据traj 对象会包含几十万行轨迹点后续聚类可能因为距离矩阵过大而卡死。我的习惯是先按月份分批计算并保存中间结果全部跑完后汇总到最终分析环节。这样即使某一天的计算中途翻车也不需要整个月重新来一遍。我自己在做某区域的多季节溯源分析时把脚本流程固定成了三步第一步按周批量生成轨迹并保存 RDS第二步对保存的数据统一做聚类和绘图快速迭代类别数第三步把聚类结果和气象观测数据合并输出各输送路径的占比时序。这套流程让我不再每次从零开始写分析代码新数据进来只需要更新日期范围半小时内就能得到一份基本报告。最后分享一个长期养成的习惯每个轨迹文件的生成过程都记录下当时的 splitr 版本、气象数据文件名、起止日期和参数清单作为分析的元数据保存。版本一换结果可能出现细微差异有了元数据才能解释这些差异的来源。HYSPLIT 模型一直在更新气象数据也在持续补全你半年前算的结果和今天算的结果可能因为底层数据不同而不完全一致这是环境数据分析的常态不是代码错误。希望这篇内容能帮你在 splitr 上少走弯路把时间花在分析本身而非排错上。本文还有配套的精品资源点击获取
