做大气污染溯源这么多年,我越来越觉得后向轨迹分析是一项“看着简单、做好很难”的工作。尤其是当你面对一场PM2.5快速升高的事件,风场却指向本地源,而浓度场又明显带有传输特征时,光靠“看风”是不够的,你需要把气团过去几天的来路完整拉出来。后向轨迹聚类分析就是解决这个问题的标准手段,而MeteoInfo与TrajStat组合,可能是目前免费工具里最适合业务化批量操作的一套方案。
这篇文章我会按自己的实际操作流程来写,从气象数据获取、站点与参数设置,到轨迹批量计算、聚类参数选择和可视化出图,全程给出可以直接复现的步骤。适合刚接触后向轨迹分析的研究生、环境监测站业务人员,以及想做污染过程复盘却不想啃HYSPLIT命令行的人。你不需要会写代码,只要跟着操作,就能把一张带聚类结果的轨迹图做出来,并且知道每个参数为什么会这样选。
1. 后向轨迹分析在污染物溯源中的应用逻辑与工具选型
1.1 为什么需要看气团“从哪里来”
后向轨迹分析本质上是利用气象场和时间推进,逆向追踪一个空气微团在过去若干小时内的运动路径。它回答的问题是:当前这个受体点上的空气,是从哪些区域过来的、经过什么路径、停留了多长时间。
你可以这样理解:假设你站在一条河边,突然看到河水变浑浊了,你想搞清楚是上游哪个支流带来的泥沙,就需要逆流而上找源头。后向轨迹就是给大气中的“河水”画出的逆流路径,只不过这条“河”流动的是空气,而追踪手段不是放漂流瓶,而是用风场数据反推。
实际案例中,一次典型的重污染过程往往包含本地积累、区域传输和二次生成三个因素。后向轨迹能对传输因素做到定量化:比如某站点72小时后向轨迹显示气团来自西北方向,途经蒙古高原和内蒙古沙地,那么这次PM10的高值就很可能包含沙尘输送的贡献。而如果轨迹显示气团长期滞留在本地或弱风速区,则更多指向本地排放积累。
需要说明的是,轨迹计算基于“准拉格朗日”假设,即空气微团随环境风场运动,没有考虑湍流扩散和化学反应对微团性质的改变。因此轨迹结果是“近似路径”,不是“真实粒子路径”,趋势判断没问题,但不要追求厘米级精度。
1.2 为什么选择MeteoInfo + TrajStat:方案优势
后向轨迹计算的底层引擎是美国NOAA开发的HYSPLIT模式,但HYSPLIT原版桌面应用偏研究型,操作路径较长,而且聚类分析、潜在源分析等模块分散在不同界面。MeteoInfo则把这些能力集成到了插件化框架中,TrajStat就是其中一个专门做轨迹分析的功能模块。
我对这个组合的评价是:它能用一套工具覆盖从轨迹计算到聚类分析,再到PSCF/CWT潜在源分析的完整链路,不需要来回切换软件,也不需要手动处理中间文件。这一点在业务化场景里非常重要,因为业务复盘往往时间紧、批次多,像吉日吉时一样地来回倒腾数据会让人失去耐心。
此外,MeteoInfo本身也是一个开源的气象数据分析与制图平台,支持脚本批处理。这意味着你不但可以交互式操作,还可以把整套流程写成脚本批量跑,这对处理一个月90天、每3小时一条轨迹的数据量来说就非常实用了。
1.3 版本选择与运行环境
MeteoInfo目前有三类版本可以选用:Java跨平台版(MeteoInfo Java)、C#版(仅Windows)以及最新的MeteoInfoLab脚本环境。TrajStat插件在Java版本中支持最完整,包括聚类分析界面和PSCF/CWT的计算模块,所以我个人的建议是优先选择Java版。
安装时需要确保系统已配置相应版本的Java运行环境(JRE或JDK),因为MeteoInfo Java版基于Java环境运行。下载完成后解压即可使用,无需复杂的安装步骤。初次打开会请求确认工作空间目录,建议设置到读写权限正常且路径不含中文的目录,这个细节看上去没什么,但能避免后面很多奇怪的坑。
版本更新方面,MeteoInfo和TrajStat都是持续更新的,新版本会修复一些地图显示问题,也会补充数据格式支持。我在实操中遇到过老版本无法读取新版GDAS数据的案例,所以建议定期去官方项目主页查看更新。
2. 数据准备:气象驱动场与站点参数设置
2.1 气象数据的来源与选择
TrajStat计算后向轨迹需要输入的气象数据,主要是全球或区域再分析格点风场数据,其中最常见的是NOAA ARL分发的GDAS(Global Data Assimilation System)数据集。GDAS数据有两种空间分辨率:一种是1度×1度,时间间隔6小时,通常表示为GDAS1;另一种是0.5度×0.5度,时间间隔3小时,简写为GDAS0p5。
选择哪个数据,主要看你研究区域的范围和目标事件的时间段。如果你分析的是城市尺度、几天内的污染过程,GDAS0p5因为时间分辨率和空间分辨率都更细,往往能更好地刻画近地面气流细节。如果你做的是季节尺度的气候统计,比如分析某个站点冬半年气团来源的总体特征,用GDAS1就足够了,而且数据量更小、处理更快。
这里有个要注意的时间覆盖问题:GDAS0p5数据并不是所有年份都有,它从某个时间点开始提供。如果你研究的是较早年份,需要确认对应时段是否已有0.5度产品,如果没有,就回退到GDAS1。我刚开始做的时候没注意到这点,下载了一个较早年月的数据发现文件目录是空的,处理到一半才发现是数据源问题。
下载地址在NOAA ARL的GDAS数据页面,文件按月份和预报时次组织。文件名大概长这样:gdas1.apr2019.w4,它表示2019年4月的第4周数据,或者gdas1.20190423.t00z,表示2019年4月23日00时次的单时次数据。TRAJSTAT读取时是直接读取一个时间范围对应的多个文件的,所以下载时要保证时间连续性。
建议把数据文件按月份放在单独文件夹中,例如:D:\traj\gdas\201911\,然后每个月一个子目录,方便以后重复使用。同一个时次的GDAS数据可能同时被多个站点或多次实验用到,保存好原始文件能避免重复下载,也方便查证。
2.2 站点表与轨迹参数设置
在TrajStat中开始轨迹计算前,需要先设置“站点文件”和“轨迹参数”。站点文件是一个包含站点名称、经度、纬度、海拔高度的表格,可以是文本文件或数据库文件格式。建议直接用界面中的“新建”功能添加站点,并逐项录入坐标。
站点的地理位置精度直接影响轨迹起始位置,所以经纬度至少要精确到小数点后两位。起始高度可以在单层和多层之间选择。研究近地面污染物受传输影响的时候,我通常设置500m这个高度层,它大致代表边界层中下部气流的平均输送高度,比10m近地面高度更能反映区域传输的特征。但如果你关注的是地表扬尘直接贡献,那10m甚至更低层的轨迹也值得看。
多层轨迹则在研究污染物垂直输送时更有优势,比如设置10m、500m、1000m三个层次,可以对比不同高度上的来向差异。但多层轨迹的聚类分析会相对复杂,因为不同高度层气团来向可能差异很大,合并聚类的解释难度会增加。
时间段和频率设置上,TrajStat支持按固定时间间隔输出轨迹。日常分析建议用6小时间隔(每天00、06、12、18时),这样一个月有120条轨迹,聚类分析样本量足够。如果你只关注某个事件日的演变,可以加密到3小时间隔甚至1小时间隔。
轨迹时长(模拟小时数)的选择取决于你的研究尺度。区域污染过程通常选择48至72小时,跨区域长距离沙尘事件则适合96至120小时。模拟时间过长,轨迹后段的不确定性会累积,而且聚类时不同轨迹之间的分离度可能被过度夸大;模拟时间过短,又可能看不到气团真正的来源区域。我的经验是日常空气质量复盘优先用72小时。
2.3 垂直运动方案与输出设置
TrajStat中与轨迹计算相关的垂直运动方案,通常可以在“Vertical Motion Method”中选择。常见选项包括:Data(使用气象数据中的垂直速度字段)、Isentropic(等熵面)、以及3D(使用三维风场分量)等。
垂直运动方式影响的核心问题是:气块在三维空间中的移动路径到底怎么确定。等熵方案假设气块沿等熵面运动,适用于描述平流层-对流层交换等绝热过程,但对一般的污染物输送来说并不完全适合。我在实际业务中多数情况使用Data或3D方案,因为它们在边界层内的表现更稳定。
输出设置中最重要的是选择“计算轨迹点文件”输出格式。TrajStat会生成每条轨迹的逐小时点位信息,包括经纬度、高度、气压等,这些是后续聚类分析和PSCF计算的输入数据。输出文件的命名规律最好自己能看懂,比如包含站点名和时间段。
另外一个容易忽略的选项是“时区”。如果站点位于东八区,而气象数据是UTC时间,那么你设定的起始时间需要先转换成UTC。TrajStat允许在界面中设定时间格式和时区,建议在计算前先确认站点时区设置与你要分析的本地时间一致,否则轨迹起点时间会整体错位,和污染过程的时间对应就会出问题。
3. 轨迹批量计算与聚类分析实操
3.1 轨迹计算实例复盘
我以一个公开常用的分析场景为例,梳理一下轨迹计算的具体操作步骤。假设研究对象是某城市环境监测站点,位置大约在经度116.4°E,纬度39.9°N,目标是分析2019年11月一次PM2.5污染过程中气团的后向轨迹特征。
首先,在TrajStat中新建站点,输入站点名、经纬度、海拔。然后在“Trajectory”模块中选择站点、起始日期(2019年11月1日00时)、结束日期(2019年11月30日18时)、间隔6小时、模拟时长72小时。气象数据格式选择ARL格式的GDAS数据,路径指向已经下载好的201911数据目录。
运行之后,TrajStat会调用HYSPLIT底层模块完成轨迹计算。这个过程的时间长短取决于轨迹条数和数据分辨率,一般一个月的数据在几分钟内可以完成。输出目录中会生成一系列轨迹文件,每个文件对应一条轨迹,文件名通常包含站点ID和时间信息。
检查输出结果时,我习惯先随机打开几条轨迹看看路径是否合理。比如,是否跨过了明显的地形障碍,是否出现起点和终点位置异常跳跃,是否有轨迹因为气象数据缺失而中断。出现轨迹中断时,需要回到数据下载环节确认对应时段的数据是否完整。
3.2 聚类参数选择与聚类实现
轨迹聚类的目标是把大量轨迹按照空间路径的相似性划分为若干类别,每一类代表一个相对固定的传输路径或气团来源方向。TrajStat内置的是层次聚类方法,这一步是整个流程中概念最密集、也最需要经验的环节。
先看距离度量。TrajStat提供了欧氏距离和角度距离两种选择。欧氏距离直接计算轨迹对应时刻点位之间的直线距离总和,对起点附近的位置差异敏感。角度距离则强调的是两条轨迹在方向上的差异,对传输方向的区分更好。我在做城市站点污染物溯源时更倾向于角度距离,因为不同轨迹如果来自同一个方向,即使末端速度有所差异,往往对应同一种天气系统主导下的输送过程。
再看聚类算法。TrajStat默认支持Ward最小方差法,它通过最小化簇内方差来合并轨迹,聚类形状比较紧凑,适合轨迹分析。实际使用中我不会把全部轨迹直接丢进去就跑,而是先设置聚类数范围,比如2到10,然后观察“总方差随聚类数变化”的结果,在拐点附近选择聚类数。
拐点的判断没有绝对标准,但有一个经验参考:当增加聚类数对减少总方差的贡献明显变缓时,说明继续分下去的信息增益有限。比如聚类数从4增加到5,总方差下降10%,而5增加到6只下降2%,那么5个聚类可能就已经足够。如果聚类数过少,不同来源的气团会被强制合并;过多则可能出现只有一两条轨迹的极小簇,结果解释会很困难。
执行聚类后,TrajStat会生成一个包含聚类结果的数据表,每条轨迹被分配到一个簇ID,并给出每个簇的轨迹数量、百分比和簇平均轨迹。此时可以在地图窗口中叠加显示所有轨迹,按簇染色查看分离效果,如果两个簇在空间上明显交叠,就要考虑调整聚类数或距离度量重新聚类。
3.3 聚类结果解读与输出
聚类完成不代表工作结束,解读才是真正体现分析价值的部分。拿到每个簇的轨迹数占比后,要把它们和观测浓度时间序列结合。比如某个站点11月PM2.5浓度高值时段,对应的是哪一类轨迹?这一类轨迹来自哪个方向、路径上有哪些工业城市或沙源地?这才是聚类分析的核心产出。
簇平均轨迹可以导出为带经纬度坐标的曲线文件,既可以在MeteoInfo中与底图叠加出图,也可以导入ArcGIS进一步制图。TrajStat中直接生成的聚类图常常自带统计标签,比如簇编号、轨迹占比,这些信息在论文或报告中都是必不可少的。
要注意的是,聚类结果受输入站点高度、轨迹时长和气象数据分辨率的影响。对同一时段的数据,把起始高度从10m改到500m,聚类结果就可能不同。所以严谨的做法是记录参数设置,并在报告里说明,让别人能看懂你的结果是在什么条件下得到的。
在输出和保存方面,我建议把每个站点的聚类结果单独存为一个工作目录,里面包含原始轨迹、聚类分配表、簇平均轨迹和聚类设置笔记。这样后续做季度汇总或年际比较时,可以直接复用底层的轨迹文件,而不必重新计算。
4. 可视化制图与成果表达
4.1 轨迹簇分布图:从底图到出图
可视化是后向轨迹分析中特别容易出效果也特别容易出错的一环。TrajStat的地图窗口可以叠加站点、轨迹、聚类结果和多种底图要素,但默认显示效果通常比较朴素,需要手工调整才能达到发表级别。
先设置底图。在MeteoInfo地图菜单中,可以加载自带的基础地理数据,也可以导入OpenStreetMap等在线底图。做报告时我一般用矢量行政区划边界叠加河流湖泊,简洁不花哨。如果在线底图加载不了,就使用本地shapefile底图,效果反而更稳定。
然后是投影设置。轨迹分析常用WGS84地理坐标系,如果研究区域纬度较高,也可以使用Lambert等角圆锥投影来减少变形。但要注意,修改投影后要检查轨迹点位是否与底图对齐,避免因为坐标系不一致导致轨迹偏离到海岸线以外。
轨迹显示上,我通常把每条轨迹设为半透明细线,聚类平均轨迹用粗实线,颜色用色盲友好的配色方案。图例中要同时包含簇编号、轨迹占比和平均高度等信息,这样阅读者不需要翻正文就能看懂每一条线条的含义。
出图尺寸和分辨率方面,业务报告建议输出300dpi以上的PNG或PDF,图片宽度控制在15cm以内。如果要在论文中配图,建议导出为矢量格式,这样放大后不会糊。
4.2 潜在源区分析:PSCF与CWT的实用细节
聚类分析告诉我们“气团从哪里来”,而PSCF(潜在源贡献因子)和CWT(浓度加权轨迹)分析可以进一步回答“可能的污染源区在哪里”。这两个模块在TrajStat中可以直接调用,计算量不大,但参数设置会影响结果可读性。
PSCF的原理比较简单:把研究区域划分为网格,统计落在每个格点内的轨迹端点数与总轨迹端点的比值,再用该格点对应的观测浓度超标次数做加权,得到一个“潜在源贡献概率”的空间分布。TrajStat中需要设置网格分辨率,我常用的经验是研究区域跨度约2000km时,网格设为0.5度左右比较合适。
网格太大会丢失空间细节,网格太小则会导致很多格点轨迹点数量稀少,计算出的概率波动很大,结果看起来像噪声。如果发现PSCF图上出现类似“椒盐噪声”的孤点,就说明网格可能太细了,可以适当放松。
CWT(浓度加权轨迹)是在PSCF基础上进一步结合浓度数值的算法。它将轨迹在格点内停留的时间与轨迹起点观测到的污染物浓度对应,计算每个格点的加权浓度贡献。相比PSCF,CWT给出的结果有浓度量级,解释起来更直观。
在TrajStat中执行PSCF/CWT后,结果通常以网格文件的图形式显示。初学者最容易忽略的是网格范围设定,如果范围设太大,边缘区域几乎没有轨迹覆盖,会出现大面积空白;如果范围设太小,又会把重要的远距离源区截掉。建议先查看研究时段内所有轨迹的覆盖范围,再根据轨迹范围设定网格边界。
4.3 批量出图与脚本化:让工作可重复
如果你只需要做一次分析,交互式操作没问题。但如果每个月都要对多个站点做一遍“轨迹聚类+PSCF”,重复劳动就会消耗大量时间。MeteoInfo真正拉开身位的地方,是它提供了脚本化接口。
你可以用MeteoInfo的脚本语言,把数据加载、站点设置、轨迹计算、聚类分析和出图整条链路串联起来。脚本之间的步骤和交互式操作完全一致,只是把界面点击变成了代码调用。我自己的做法是先交互式调好一张图,确认参数无误后再把对应操作改写为脚本,保存为一个函数,后续只需要替换站点和时间段即可复用。
批量出图时要注意输出文件名的规划,建议包含站点名、时间范围和分析类型,例如“siteA_201911_cluster.png”。另外,对脚本中用到相对路径还是绝对路径要有统一约定,否则换一台电脑运行脚本就会报路径找不到。
5. 常见问题与排查技巧实录
5.1 数据读取失败与轨迹数量不足
这是最容易被卡住的一步。TrajStat读取GDAS数据报错,或某一段时间的轨迹没有生成,大概率是数据文件缺失或文件名不匹配。我曾遇到一次性下载了整月数据,却因为某一天的数据文件只有1KB大小,导致轨迹计算到那天时直接中断。
排查思路其实很简单:先看报错的时间点,再检查该时间点对应的GDAS文件是否存在、大小是否正常。有时候NOAA服务器的文件有延迟,某一天的数据当天下载不到,需要改天补下。还有一个常见问题是文件路径含中文或空格,TrajStat在读取时会解析异常,建议数据目录统一用英文字母加数字的命名。
如果轨迹生成成功但数量远少于预期,检查时间设置是不是把结束时间设成了与开始时间重合或早于开始时间。另一个常见原因是,TrajStat中默认的时间步长单位是小时,如果填入了类似“1440”这样的数值,就会把轨迹模拟时长的单位误当成小时而不是分钟,轨迹就会非常长,且数量异常。
5.2 聚类结果“不合常理”的调试路径
聚类结果出来后,有时会发现某一个簇的轨迹几乎平行于海岸线或跨越了明显不连续的风场特征,这时先不要急着调整聚类数,而是检查垂直运动方案和数据源。不同垂直运动方案导致轨迹在垂直方向上的路径差异很大,这种差异在高空风与地面风不一致时会直接反映到聚类上。
还有一种典型情况:某个簇只含少数几条轨迹,且这些轨迹恰好是某一天气象数据缺失后生成的短轨迹。这种“伪簇”必须剔除。建议在聚类前就对轨迹长度做筛选,比如设定轨迹有效点数必须大于总时长的80%,低于该比例的轨迹不参与聚类。
另外,如果你做了季节性汇总分析,需要保证每个季节的轨迹样本量接近。如果春季只有几天有效轨迹,与冬季120条轨迹合并聚类,春季气候特征几乎会被淹没。这种情况下最好按季节分别聚类,而不是混在一起做总体聚类。
5.3 制图过程中坐标与底图对齐问题
MeteoInfo中底图突然不显示、轨迹漂到海上,或者行政边界与轨迹交错但形状明显错位,这些视觉问题多数是投影坐标与地理坐标不统一导致的。默认情况下轨迹输出是经纬度坐标(WGS84),而底图的投影坐标如果设置成了UTM,就必然会出现偏移。
解决办法是先在MeteoInfo的“Map Properties”中把投影参数统一为经纬度,或者统一到研究的中心经纬度对应的投影参数。修改投影后建议先放大到站点附近确认轨迹起点与站点图标重叠,再检查远处轨迹是否与海岸线吻合。
在线底图的缓存也是一个容易被忽略的问题。如果之前加载了一张旧范围的地图,再次分析时动态范围变化后,缓存可能导致底图显示残影或空白。此时清理地图缓存或重启MeteoInfo即可解决。
5.4 常见问题速查表
| 问题现象 | 可能原因 | 快速排查建议 |
|---|---|---|
| 轨迹计算中断 | 气象数据文件缺失或损坏 | 检查对应时段数据文件大小与完整性 |
| 轨迹明显偏长或偏短 | 时间步长或模拟时长设置错误 | 核对基准时间、模拟小时数与时间间隔 |
| 聚类中出现极少轨迹的簇 | 数据缺失产生短轨迹 | 先过滤短轨迹再聚类 |
| PSCF结果噪声大 | 网格分辨率过高 | 调粗网格,或平滑处理 |
| 轨迹与底图错位 | 投影坐标不一致 | 统一为WGS84经纬度坐标 |
| 在线底图不显示 | 网络问题或缓存异常 | 改用本地矢量底图 |
作为一个长期使用这套工具做业务复盘的人,我最后想分享一个体会:这套流程真正的准入门槛不在软件操作,而在于你是否理解每个参数背后的物理意义。MeteoInfo和TrajStat把HYSPLIT的底层能力包成了一个友好的外壳,缩短了从数据到结论的时间,但如果你不清楚“风场分辨率对轨迹的影响”“等熵方案在这种地形下为什么不合适”,一旦结果看起来不合理,你就很难判断是数据问题、参数问题,还是分析逻辑问题。
我自己的习惯是,在正式批量出图之前,先手工做一条轨迹、聚一次小规模的类,用肉眼验证合理性。这个过程只需要十几分钟,却能帮你避开很多后知后觉的大坑。把这套流程固化下来,你会发现后向轨迹分析不再是某个季节才想起来用的“高深工具”,而是像看风场图一样,成为日常分析里随时可以调用的基本能力。