从一条轨迹说起。做大气污染溯源的时候光看单条HYSPLIT后向轨迹意义有限——某一口气团只能代表那个时刻那个高度的输送路径没法回答“这个季节污染物主要从哪来”这种大问题。我第一次做污染气团来源分析时以为把几十条轨迹叠在一张图上就能交差结果图一出来自己都懵了几百条彩色线缠成一团根本看不出规律。后来老老实实学了后向轨迹聚类分析才算是把“散点轨迹”变成了“可解释的输送通道”。这里把整套流程从数据准备到出图踩过的坑写出来给要用MeteoInfo和TrajStat做轨迹聚类的朋友省点时间。这套工具链的核心思路其实不复杂用MeteoInfo读入气象再分析数据用TrajStat插件对批量HYSPLIT轨迹做插值、聚类最后把聚类结果叠加在地图上做可视化。关键是每一步都有不少“不说清楚就得自己摸半天”的细节尤其是数据格式校验和聚类数目选择这两块。1. 后向轨迹聚类到底在解决什么问题1.1 单条轨迹与聚类结果的本质区别后向轨迹模拟的是气团在过去一段时间内的移动路径。HYSPLIT模型根据气象场插值计算气团轨迹输出每个时刻气团所在的经纬度和高度。这条轨迹本身是“个体”而聚类分析是把大量轨迹按照空间形态的相似程度进行分组让相似路径归为同一类最终得到几条代表性的“平均轨迹”。打个比方单条轨迹像一个人上班的路线聚类结果像整个城市居民通勤的主流走向。你关心的是城市交通整体规划而不是某个人的具体路线。污染溯源场景里决策者需要知道的是“这个季节气团主导来源是哪几个方向”而不是“1月3日那天气团经过了哪里”。1.2 实际应用场景与适用范围聚类分析结果可以用来支撑以下几类工作污染过程成因分析识别重污染期间气团的主要输送通道判断是本地积累还是区域输送主导。季节传输特征对比按季节分别聚类比较不同季节的气团来源差异对应到沙尘、霾、臭氧等污染类型。潜在源区分析的前置步骤聚类结果可以作为条件概率函数PSCF和浓度权重轨迹CWT分析的分组依据。水汽来源分析在降水同位素研究中聚类后的轨迹用于判断水汽来源方向。适用对象主要是环境监测、气象科研、高校课题组里需要做溯源分析的人。不需要太深的数学基础但要能理解聚类算法的基本逻辑否则容易出现“聚类数设置不合理但结果照用”的问题。2. 数据准备阶段的高频事故区GDAS数据获取与轨迹文件生成2.1 气象数据的选择与下载渠道HYSPLIT模型本身不自带气象数据需要先下载再喂给它。NOAA ARL提供全球数据同化系统GDAS数据是目前最常用的驱动数据。GDAS数据主要分两类数据名称空间分辨率时间分辨率适用场景GDAS11°×1°3小时常规后向轨迹模拟精度足够GDAS0p50.5°×0.5°3小时局地精细化模拟数据量大很多数据从NOAA ARL的FTP站点下载文件名格式类似gdas1.apr18.w1代表2018年4月第1周w1~w6每5天一个文件6个文件覆盖一个月。GDAS0p5文件名是gdas0p5.20180401这种按天命名的格式。需要注意几点FTP下载建议用专门的下载工具断点续传功能很重要GDAS1单月数据量大概2~3GBGDAS0p5单月能到几十GB网络不稳定时中断很常见。数据要覆盖模拟时间窗口。做72小时后向轨迹起始时间是4月18日00时那么气象数据必须覆盖4月15日00时到4月18日00时往前推3天不能只下载起始日当天。跨月、跨年时要把相邻月份的数据都下载完整文件缺了模型运行会直接报错而且报错信息不一定能立即定位到是缺数据。2.2 用HYSPLIT生成批量轨迹文件轨迹来源有两种路径用NOAA READY网站在线生成或者在本地装HYSPLIT PC版批量跑。在线版适合单次少量轨迹本地版适合批量任务。无论哪种方式关键参数设定要提前想清楚起始高度近地面污染分析一般选500m AGL离地高度这个高度能代表边界层内气团的整体输送特征。要做高空输送分析可以选1000m或1500m。不同高度的轨迹不能混在一起聚类必须分开处理。模拟时长常用72小时3天也有做48小时或120小时的。回溯时长越长气团来源范围越大但轨迹误差也越大因为长时间积分累积误差。常规污染过程分析72小时足够。起始时间与频率逐日逐时或逐日固定时刻。做季节聚类通常选择每日00时或08时当地时间起始一天一条一个季节约90条轨迹。做污染过程分析可以加密到逐小时。生成后的轨迹文件通常是tdump格式的文本文件里面按时间顺序记录了每个时刻的轨迹点坐标。文件头部有起始点的经纬度、起始时间、模拟时长等信息。这个文件就是下一步TrajStat要读入的对象。2.3 一个常被忽视的格式问题HYSPLIT在线版生成的tdump文件和第16版PC版生成的tdump文件在表头细节上略有差异。TrajStat在读取时对文件格式有一定要求最容易出问题的是以下两点年份字段位数某些版本的tdump文件里年份只写2位如18有些写4位如2018。TrajStat老版本读2位年份会出现日期解析错误表现为轨迹点的时间轴错乱。缺失值表示轨迹模拟中如果某一时刻气象数据缺失tdump文件里可能出现-9999之类的填充值。插值前如果没有检查这些异常值聚类结果里会出现飞到非洲的诡异轨迹。我的习惯是拿到tdump文件先打开看几行确认表头格式正常、数据行完整再用脚本把异常值轨迹挑出来剔除。这一步虽然不复杂但能省掉后面排查结果的很多麻烦。3. MeteoInfo与TrajStat的环境配置和插值预处理3.1 软件版本匹配问题MeteoInfo是国产开源气象软件由气象科研人员开发在环境领域用得很多。TrajStat是它的一个插件专门做轨迹分析。版本匹配是个比较容易踩坑的地方MeteoInfo 1.x版本对应TrajStat插件需要在软件内通过“插件管理”加载菜单路径是Tools - TrajStat。MeteoInfo 2.x版本Java版中TrajStat的加载方式略有变化需要确认插件包版本与主程序版本一致否则会出现菜单项不显示的问题。新版MeteoInfo3.x以后直接内置了部分轨迹分析功能但很多教程还是基于1.x写的界面路径对不上会让新手很困惑。我自己用下来最稳妥的组合是MeteoInfo 1.4.x加配套的TrajStat插件包。这个版本资料多、稳定遇到问题容易搜到解决方案。装好之后在MeteoInfo主界面能看到TrajStat菜单就说明插件加载成功。3.2 插值处理为什么不能直接用原始轨迹聚类HYSPLIT输出的轨迹点在时间上并不是等间隔的。在路径平直的地方模型输出的轨迹点可能很稀疏在气团转向或地形复杂的区域轨迹点会很密。直接拿原始轨迹做距离计算有两个问题不等间隔的轨迹点会使得密集区域对距离计算的贡献更大导致聚类结果偏向局部特征。不同轨迹之间的对应点高度不一致距离计算时如果包含了高度维度的差异权重分配很难统一。TrajStat的轨迹插值功能就是把每条轨迹重新采样到固定时间间隔。具体操作路径是Tools - TrajStat - Interpolate Traj对话框中设置时间间隔一般选1小时。模拟72小时的轨迹插值后就是73个点含起始点。高度插值方式可以选择固定高度或固定气压。固定高度AGL常用于近地面分析固定气压用于高空分析。插值后生成的文件会保存在你指定的输出目录文件名通常是在原始文件名基础上加上_interp后缀。注意插值过程和聚类过程都要指定同一目录否则后一步找不到前一步的结果文件。3.3 插值前检查轨迹文件数量批量插值前我建议先确认待处理的tdump文件是否都在同一目录下而且目录路径不要包含中文字符。MeteoInfo底层是Java写的某些版本对中文路径的支持不好容易在读取文件时抛出超长路径或编码异常。路径问题看起来低级但实际遇到的人真不少尤其是Windows环境。还有一个问题是轨迹文件特别多的时候插值过程会有点慢。一个季节90条轨迹、每条73个点插值操作一般几秒钟就能跑完如果做逐小时轨迹比如一个月720条内存占用会明显上升。遇到这种情况可以分批插值最后在聚类时把多批结果一起加载不影响后续流程。4. 聚类数目怎么定TSV拐点判读与物理解释校验4.1 TrajStat聚类的算法逻辑TrajStat提供了两种聚类算法K均值K-means先随机指定K个聚类中心然后迭代优化每个轨迹到所属聚类中心的距离总和。计算速度快但需要预先指定聚类数K。Ward最小方差法这是一种层次聚类方法每次合并两个能使类内方差增量最小的类最终形成一个层次树。通过树状图可以直观判断聚类数目不需要预先指定。轨迹聚类的“距离”并不是简单的直线距离而是两条轨迹对应轨迹点之间的空间距离经度差、纬度差、高度差加权的累积。TrajStat默认使用欧氏距离平方。高度差是否纳入距离计算、权重是多少可以在设置里调整对聚类结果有影响。4.2 聚成几类才合适这个问题的标准答案是“看总空间方差的变化”。TrajStat在聚类之前会先计算不同聚类数对应的总空间方差TSV并生成一个TSV随聚类数变化的表格和曲线图。随着聚类数K增大TSV单调下降——因为分得越细组内差异越小。关键是找拐点在拐点之前增加聚类数能显著降低TSV拐点之后增加聚类数收益很小。实际操作中我的做法是先让TrajStat计算K2到K10的TSV。查看TSV曲线找到斜率从陡峭变平缓的位置。把这个位置的K作为候选聚类数。结合物理解释做最终判断。比如分4类时每类对应一个清晰的气团来源方向西北、偏北、偏东、局地分5类时其中两类路径非常接近、物理意义重叠那就取4类。要注意TSV拐点有时候并不明显尤其是气团来源本身就比较分散的情况。这时候优先保证聚类结果的物理可解释性而不是机械地找拐点。4.3 聚类结果文件里有什么聚类完成后TrajStat会输出几类文件聚类信息表每条轨迹属于哪一类。平均轨迹数据每类轨迹的代表性路径平均轨迹。统计汇总各类轨迹数量、占比、平均长度等信息。地图显示文件可以在MeteoInfo中直接查看聚类轨迹的空间分布。这些文件是后续做可视化制图的数据基础建议在聚类前就规划好输出目录避免后期找不到结果散落在哪里。5. 可视化制图把聚类结果画成发表级的轨迹图5.1 在地图上叠加聚类轨迹MeteoInfo提供了灵活的地图绘制功能核心操作逻辑是图层Layers概念。TrajStat聚类结果会作为轨迹图层加载到地图上。基本步骤打开MeteoInfo在工具栏选择添加底图图层海岸线、国界、省界。在Tools - TrajStat菜单下加载聚类结果文件。在图层属性中设置各类轨迹的颜色、线宽、透明度。调整地图投影和显示范围让轨迹覆盖区域合理显示。常用的地图底图数据有两种来源MeteoInfo自带的map资源库或者外部导入Shapefile格式的边界文件。做中国区域的分析我通常会导入省界和市界叠加让轨迹路径和行政区划的关系一目了然。5.2 出图细节规范一幅能放进论文的轨迹聚类图不只是把线画出来就行。几个关键细节颜色区分度每类轨迹用不同颜色注意色盲友好配色红绿搭配尽量避免。图例信息标注各类轨迹占比比如“类别132%”。图例字号要跟图面信息匹配不要过大或过小。起止点标记研究区域的站点位置用醒目的符号标出。比例尺与指北针MeteoInfo的制图布局中可以添加比例尺和指北针具体在Layout下调整。导出格式导出版本要选高分辨率位图TIFF或PNG300dpi以上或矢量图PDF或SVG矢量图文字清晰度最好。5.3 轨迹密度检验老手出图之后还会做一个检查把原始轨迹按类别分别画出来看看聚到同一类的轨迹内部一致性如何。如果某一类内部轨迹非常发散说明这一类本身代表的意义有限可能需要调整聚类数或重新审视高度选择。这个检查在MeteoInfo里操作起来很简单就是把插值后的轨迹文件按类别筛选出来分开加载到地图上查看。看起来很笨但能避免聚出一类平均路径很漂亮、实际却包罗万象的“垃圾类”。6. 实战中的报错盘查与结果解读经验6.1 高频报错及排查链路下面是我这几年用MeteoInfo和TrajStat实际遇到频率最高的几类问题按排查顺序列出来现象可能原因排查步骤轨迹文件加载后地图无显示轨迹起点坐标超出地图范围检查tdump文件头部经纬度核对站点坐标是否写反经度纬度顺序插值过程报错“数据格式错误”文件编码或表头不符合TrajStat要求用文本编辑器打开tdump文件比对表头行数据样例菜单Tools下找不到TrajStat插件未正确加载或版本不匹配检查插件目录是否放置正确重新启动MeteoInfo聚类结果分类数异常所有轨迹分到一类轨迹空间跨度太大或高度权重设置不当检查起始高度是否统一尝试调整距离计算参数出图时中文标注乱码字体配置问题在MeteoInfo中设置中文字体或改用英文标注排查问题有个原则先从数据自身找原因再从软件参数找原因。我见过不少同行遇到聚类结果异常第一反应是软件坏了花时间重装最后发现是轨迹文件本身有缺失值。6.2 结果解读的统计口径聚类结果出来后常规的统计分析是把每个类别的轨迹占比与观测污染物浓度做对比。做法是把污染浓度数据按时间与轨迹起始时间对齐。按轨迹类别分组统计浓度均值、中位数、超标率。比较各类别污染物浓度差异结合输送路径方向给出解释。这里要强调的是相关性不等于因果性。某类轨迹对应高浓度污染只能说这类气团输送条件更容易造成污染积累不能说污染就是来源地直接排出来的。严谨的表述是“该类气团对应的污染物浓度显著高于其他类别可能存在区域输送贡献”。6.3 季节对比的标准化流程做季节聚类分析时建议春3-5月、夏6-8月、秋9-11月、冬12-2月分别聚类而不是把全年数据混在一起聚类。原因很简单不同季节的盛行风向差异很大混合聚类会把季节差异淹没在分类平均中物理意义不清晰。季节对比的图件编排建议采用四联图或2×2网格布局同一研究站点相同聚类设置参数保证结果可比性。四季轨迹方向差异往往一眼就能看出规律这也是这类图在论文里很讨喜的原因。7. 从聚类图到完整溯源结论的最后一公里聚类分析完成后很多分析还会继续往前走一步做潜在源区贡献PSCF或浓度权重轨迹CWT分析。TrajStat本身就支持这两个功能输入聚类结果和污染物浓度序列输出源区贡献网格图。PSCF的核心思想是统计污染轨迹经过每个网格的次数与总轨迹经过次数的比值识别高概率源区CWT则是用浓度值加权能够区分高浓度来源和低浓度来源。这两张图叠加轨迹聚类结果就是一份比较完整的污染溯源图件组合。不过要提醒一句PSCF和CWT的网格分辨率设置需要根据轨迹覆盖范围来确定网格太小会导致很多网格轨迹频数很低、结果噪声很大网格太大又会丢失空间细节。常规做法是先用较粗的网格比如0.5°试跑看覆盖效果再细化。我自己做完一整套流程之后的感受是MeteoInfo和TrajStat这套工具链的学习曲线并不陡峭真正花时间的是理解数据本身的特征和算法参数的含义。很多参数比如起始高度、回溯时间、聚类数目没有绝对正确的答案需要你根据研究问题的物理背景去判断。多试几组参数对比结果比死记硬背教程参数可靠得多。最后分享一个小技巧做完聚类出图后把每次使用的参数记录在分析笔记里包括数据版本、模拟时长、起始高度、插值间隔、聚类方法、聚类数目、TSV拐点截图。这个习惯在写论文方法论部分的时候会救命因为审稿人几乎必问这些细节到时候再回头翻操作记录就很从容。