搞大气水汽反演的人应该都有过这种纠结探空站一天就放两次气球时间分辨率低到让人抓狂再分析资料空间分辨率虽然过得去但时效性和精度没法保证。GNSS连续运行参考站网把这个缺口补上了大半——卫星信号穿过大气层时受到的对流层延迟这个定位中的“噪声”被当成信号来用反演出的GNSS PWV时间分辨率能做到5分钟甚至更高。而要把这件事做准天顶总延迟ZTD的解算精度是第一道门槛。我这两年一直在用武汉大学PNT中心开源的PRIDE PPP-AR来处理ZTD解算和PWV反演踩过不少坑也积累了一套比较顺手的全流程整理出来给同行参考。这篇文章不会只贴命令会把参数背后的逻辑、输出的含义、精度评估里容易出问题的地方都说清楚适合刚接触GNSS气象学、准备用PPP-AR做水汽时间序列的研究生和工程师。1. 为什么是PRIDE PPP-AR这条技术链路的价值与选型逻辑1.1 GNSS水汽反演到底解决了什么痛点水汽是大气中最活跃的成分虽然只占大气总量的0.1%到3%但它的时空变化直接决定了强对流天气的发生发展。传统观测手段里探空站一天只放两次气球时间分辨率12小时遇到天气系统快速演化的过程基本抓不住微波辐射计精度尚可但设备贵、站点稀少再分析资料像ERA5虽然时空分辨率不错但它是模式输出结果时空分辨率受限于同化系统和观测密度在局地对流过程面前仍然会“钝化”。地基GNSS反演PWV的出现相当于给水汽观测网格里补上了高时间分辨率的一环。只要有一个连续运行的GNSS测站卫星信号穿过大气层时产生的天顶总延迟ZTD就可以按分钟级解算出来从中剥离出湿延迟再乘上一个转换系数就能得到大气可降水量PWV。这个值与探空、微波辐射计的结果有很好的一致性而且全天候、不受云雨影响、设备维护成本低。现在很多气象业务单位把GNSS PWV纳入短临预报的参考因子尤其是强降水过程前水汽快速增加的特征GNSS时间序列能给出非常有价值的信息。1.2 为什么选PRIDE PPP-AR而不是GAMIT和Bernese做ZTD解算的软件不少常见的路线大概有三类。我只说我实际用下来的感受不做盲目吹捧。GAMIT/GLOBK是相对定位模式需要多站联测做网解解算前要挑选基准站、做基线处理、处理周期模糊度整个流程比较重。如果只是想要单站或者十几个站的水汽时间序列用GAMIT会有“杀鸡用牛刀”的感觉而且计算量大、学习曲线陡。Bernese是功能最全面的学术软件之一从PPP到网解、从对流层到电离层都能做但操作门槛更高命令行体系、脚本配置、文件组织方式都很“学院派”新手从零到跑通第一个解算流程通常要花一到两周。PRIDE PPP-AR的特点在于它走的是精密单点定位加模糊度固定路线单站就能解算不需要联测基准站支持GPS、GLONASS、Galileo、BDS多系统联合解算数据利用率高编译和使用相对轻量批量处理测站很方便最关键的是它对ZTD这类对流层参数的估计精度与GAMIT和Bernese的结果在毫米级能够保持一致。对于做水汽研究的人来说这意味着不用维护一套复杂的网解环境也能拿到高质量的对流层延迟产品。另外这几年PRIDE PPP-AR一直在更新支持RINEX 3.x格式、支持多系统多频率模糊度固定对北斗数据也有比较好的处理能力。国内做GNSS气象应用的研究组用得越来越多社区里的经验积累也相对多一些遇到问题容易找到人讨论。1.3 整条技术链路长什么样从原始数据到最终的PWV精度评估完整链路可以拆成四个环节数据准备GNSS观测文件RINEX格式、精密星历、精密钟差、DCB文件、天线相位中心文件。ZTD解算用PRIDE PPP-AR做精密单点定位解算逐历元输出天顶总延迟ZTD。PWV反演从ZTD中减去干延迟ZHD得到湿延迟ZWD再通过加权平均温度Tm和转换系数Π把ZWD换算成PWV。精度评估把GNSS PWV与探空资料或ERA5再分析资料的PWV做对比计算偏差、标准差、均方根误差和相关系数。这个链路每一步都有坑数据下载缺文件、配置参数选错、干延迟扣除精度不足、Tm模型偏差过大、时间系统没统一等等。后面我会一个一个展开说把我实际遇到过的情况和解决办法都写出来。2. 环境与数据从源码编译到观测文件准备2.1 编译前的系统环境配置PRIDE PPP-AR目前主要面向Linux环境Windows下虽然可以用WSL或者虚拟机跑但原生Linux还是最省心的。我用的是Ubuntu 20.04长期支持版内核和编译工具链都比较稳定。新装完系统后先装好基础依赖再编译不然会在configure阶段反复报错。sudo apt update sudo apt install build-essential gfortran make wget curl sudo apt install libcurl4-openssl-dev libxml2-dev这里有几个容易忽略的点。第一gfortran版本不要太老也不要太激进Ubuntu 20.04自带的gfortran 9处理PRIDE PPP-AR的代码没问题但如果你用的是CentOS 7自带的老版本gfortran某些新语法会编译失败。第二如果源码里使用了网络下载模块libcurl开发库是必须的缺了这个头文件会在编译中报curl/curl.h: No such file or directory这种错误。第三建议装一个nco和cdo工具包后面处理ERA5 NetCDF数据会用到虽然不属于PRIDE PPP-AR的依赖但在精度评估阶段几乎是刚需。2.2 获取源码与编译PRIDE PPP-AR的源码托管在GitHub上搜索PRIDE-PPP-AR就能找到也可以从武汉大学PRIDE实验室官网的链接跳转。下载后进入源码目录按标准流程编译git clone https://github.com/PRIDE-Lab/PRIDE-PPP-AR.git cd PRIDE-PPP-AR ./configure --prefix$(pwd)/install make -j4 make install--prefix参数指定安装路径这是我一直建议的做法把编译产物和源码分开后续更新版本时不用重新编译旧目录只需要切换路径就行。-j4表示四核并行编译如果是八核机器可以改成-j8能省不少时间。编译成功的标志是src目录下生成了主程序pridepppar或者你设置了--prefix后在install/bin下能看到它。如果遇到gfortran: error: unrecognized command-line option之类报错基本可以判定是编译器版本问题换gfortran 9或10重新编译即可。一个值得注意的经验PRIDE PPP-AR的代码对编译选项比较敏感如果configure阶段出现奇怪的错误建议先make distclean清理干净再来一遍不要硬着头皮在报错状态下反复修改makefile。2.3 数据准备观测文件、精密星历与钟差下载ZTD解算需要三类核心输入测站观测文件、精密星历、精密钟差。观测文件从IGS数据中心下载CDDIS是常用的来源但国内访问有时不太顺畅中国这边的替代方案是武汉大学IGS数据中心igs.gnsswhu.cn下载速度明显更好。RINEX文件的命名遵循统一规则例如BJFS00CHN_R_20240010000_01D_30S_MO.rnx表示2024年第1天北京房山站BJFS的30秒采样、24小时观测文件。精密星历和钟差产品按时效分为最终、快速、超快速三类。做研究用最终产品精度最高但发布时间有一到两周延迟近实时应用用快速产品钟差精度在亚纳秒量级实时应用用超快速预报产品ZTD精度会差一些。IGS的最终产品命名类似igs21927.sp3和igs21927.clk前者是SP3格式的精密轨道后者是精密钟差文件这个21927对应GPS周和星期几解算前要确保把当天对应的所有历元覆盖。DCB文件差分码偏差也是容易漏掉的一项。多系统PPP中不同频率、不同系统间的码偏差如果不改正会影响模糊度固定和时间参数估计。CODE分析中心发布DCB产品格式是CODE20240010000.DCB这种命名。同时还需要天线相位中心文件igs20.atx这些文件可以和精密产品一起放在一个公共目录里各测站解算时共用。我建议在本地建一个清晰的数据目录结构至少包括obs/、orbit/、clock/、dcb/、atx/这几个子目录名称固定下来不要变。PRIDE PPP-AR运行时会按配置文件里的路径去读文件目录结构越规律写批处理脚本越省事。3. ZTD解算实操配置参数背后的原理与输出文件解读3.1 解算前先理解ZTD是怎么被估计出来的在PPP解算里ZTD不是直接观测量而是作为未知参数参与平差。卫星信号穿过对流层时产生的延迟大致分成两部分静力学延迟干延迟ZHD和湿延迟ZWD。干延迟变化比较有规律可以先通过模型扣除一部分湿延迟变化快、与空间位置关系密切是残余估计的主角。PRIDE PPP-AR的处理逻辑是先用Saastamoinen模型给出一个先验ZTD再用VMF1或GMF映射函数把天顶方向的延迟映射到卫星信号的实际传播路径上路径上的误差被分解成水平梯度和残余天顶延迟两部分。残余天顶延迟按随机游走过程进行逐历元估计后面再加上模糊度固定把载波相位模糊度固定到整数从而显著提高ZTD估计的稳定性和精度。这就解释了为什么模糊度固定对ZTD这么重要不固定模糊度时模糊度参数和天顶延迟参数之间的相关性会通过观测方程传染ZTD序列容易出现小幅漂移模糊度固定后参数之间的耦合被切断ZTD时序就变得干净很多。所以用PRIDE PPP-AR时我的建议是务必开启模糊度固定功能否则软件最大的优势没有发挥出来。映射函数的选择也要留意。VMF1和GMF各有优劣VMF1使用的是ECMWF数值天气预报产品计算的对流层映射函数精度略好但需要外部下载VMF1格网产品GMF是经验模型不依赖外部数据穿透性和稳定性好。如果你做的是高精度水汽研究我推荐VMF1如果只是批量处理常规站GMF完全够用。3.2 配置文件逐项解读PRIDE PPP-AR通过一个文本配置文件控制解算行为不同版本的配置格式略有差异但核心参数是相通的。以下是我自己常用的一套配置内容标注了每个参数的实际意义# 测站与数据路径 site_name BJFS obs_file ./obs/BJFS00CHN_R_20240010000_01D_30S_MO.rnx nav_file ./obs/BJFS00CHN_R_20240010000_01D_MN.rnx sp3_file ./orbit/igs21927.sp3 clk_file ./clock/igs21927.clk dcb_file ./dcb/CODE20240010000.DCB atx_file ./atx/igs20.atx # 解算策略 mode ppp-ar interval 30 # 观测值采样间隔单位秒 elevation_mask 7 # 截止高度角单位度 system grce # 参与解算的系统gps/glonass/galileo/bds # 对流层估计 trop_model saastamoinen mapping_function vmf1 trop_interval 30 # 天顶延迟估计间隔单位分钟 gradient_est yes # 估计水平梯度interval 30对应观测文件里数据本身的采样率如果你的RINEX文件是30秒采样这里就必须写30不能随意改。elevation_mask 7是我批量处理时比较折中的方案低于7度的观测值受多路径和大气折射影响较大太高了又会损失低仰角卫星带来的水汽信息。低仰角卫星信号穿过大气的路径更长携带的水汽信息也更多对ZTD估计是有利的但要付出噪声增大的代价。gradient_est yes建议一直打开。天顶延迟是一维的但大气水汽在水平方向往往不均匀尤其是遇到锋面过境时东西方向和南北方向的湿延迟差异很大。如果忽略水平梯度天顶延迟估计会被污染表现就是解出的ZTD序列在天气系统过境时出现系统性偏差。trop_interval 30表示每30分钟估计一个天顶延迟参数。这个参数决定ZTD时间分辨率想要更高分辨率的PWV可以改成15甚至5但要注意时间分辨率越高单个历元的观测冗余越少ZTD噪声会增大。做常规水汽研究15到30分钟是平衡点。3.3 运行解算与输出文件解析配置写好后运行主程序解算单个测站./pridepppar config_BJFS.in程序运行时会打印每个历元的处理进度同时输出一个后缀为.ztd或类似命名的文件里面按时间顺序列出每个历元的ZTD估计值。典型输出列的格式大致是年、年积日、秒、天顶总延迟ZTD单位米、估计标准差、水平和垂直方向的梯度参数等。拿到输出后第一件事不是直接算PWV而是看ZTD时间序列是否平滑连续、有没有跳变。正常情况下ZTD是缓慢变化的大气信号数值一般在1.5米到2.7米之间取决于海拔和季节如果出现单历元突然跳变几厘米、甚至到负值或者0的情况多半是观测文件中存在粗差或者模糊度固定失败需要把这一段标记为无效不要进入后续反演。我习惯在批量处理完后写一个质量控制脚本把ZTD序列中偏离前后中位数超过3倍标准差的历元剔除掉同时删除解算状态标志不是正常收敛的历元。这一步看似简单但它决定了PWV时间序列的“干净程度”不能省。4. PWV反演从湿延迟到可降水量的完整脚本4.1 干延迟与湿延迟分离ZTD解算出来后得到的是天顶方向的总延迟要得到湿延迟ZWD需要先扣除干延迟ZHD。干延迟可以用地面气压和测站位置计算最常用的是Saastamoinen模型ZHD 0.0022768 × Ps / (1 - 0.00266 × cos(2φ) - 0.00028 × H)其中Ps是测站处地面气压单位hPaφ是测站纬度H是测站椭球高单位km。从公式可以看出干延迟与气压几乎成正比气压1百帕的误差会引起大约2.3毫米的干延迟误差进而影响同量级的PWV误差所以地面气压的精度很关键。最佳做法是在测站附近布设气象传感器实时获取气压、温度、湿度。但现实往往是测站没有气象观测设备这种情况下可以用数值预报模型的气压产品。欧洲中心ERA5的逐小时2米气压或者GPT2w、GPT3经验模型内插的气压值都能用。我自己实测下来GPT3模型给出的气压与真实气压差通常在1到2百帕以内换算成PWV误差大约在0.5到1毫米对大多数研究场景可以接受。4.2 加权平均温度Tm的获取Bevis公式与GPT2w/GPT3从湿延迟换算可降水量核心转换系数里有一个变量叫加权平均温度Tm它的定义是大气柱内水汽压和温度的加权平均。Tm不能直接观测但可以经验估算。最经典的方法是美国学者Bevis提出的线性回归公式Tm 70.2 0.72 × TsTs是测站处地表温度单位K这个公式是从美国地区探空资料拟合出来的在北半球中纬度地区表现很好。但如果你做的区域是热带、高原或者高纬度地区直接用Bevis公式可能带来1到3开的偏差而Tm每偏差1开PWV大约有1%的相对误差对于20毫米的PWV就是0.2毫米做高精度研究时不可忽视。更好的方案是用GPT2w或GPT3经验模型它们按经纬度和高程给出Tm的长期平均值和周期变化项输入测站坐标和年积日就能插值出来空间代表性比单点的Bevis公式更好。用GPT系列模型不需要外部实时数据而且它在全球范围都做过验证。我个人的处理流程是优先用测站气象观测数据按Bevis公式计算Tm没有气象数据就用GPT3模型获取表面温度和Tm如果做批量历史数据处理也可以直接插值ERA5气压层温度和水汽剖面来计算Tm精度更高但计算量也更大。4.3 PWV计算与Python脚本实现转换系数Π的公式如下Π 10^6 / (ρ_w × R_v × (k3 / Tm k2))其中ρ_w是液态水密度取1000 kg/m³R_v是水汽气体常数取461.495 J/(kg·K)k2和k3是大气折射率常数常用的一组取值是k2 22.1 K/hPak3 373900 K²/hPa。注意不同文献里常数取值略有差异选一组并保持全程一致就好换成不同常数组会让PWV产生约0.5%到1%的差异。PWV Π × ZWD下面是完整的Python计算脚本注释已经写得比较详细直接抄作业即可import numpy as np import pandas as pd # 读取ZTD解算结果 # 假设文件列为year doy sod ztd ztd_std ... df pd.read_csv(BJFS.ztd, delim_whitespaceTrue, names[year, doy, sod, ztd, std, grad_n, grad_e], skiprows1) # 测站参数 lat 39.6 # 纬度单位度 h_ell 80.0 # 椭球高单位 km北斗站约80m 0.08km P_s 1013.25 # 地面气压单位 hPa来自气象数据或GPT2w T_s 288.15 # 地表温度单位 K # 1. 干延迟 ZHDSaastamoinen 模型 phi np.deg2rad(lat) ZHD 0.0022768 * P_s / (1 - 0.00266 * np.cos(2 * phi) - 0.00028 * h_ell) # 2. 湿延迟 ZWD df[ZWD] df[ztd] - ZHD # 3. 加权平均温度 Tm这里用Bevis公式可替换为GPT2w插值结果 Tm 70.2 0.72 * T_s # 单位 K # 4. 转换系数 Π rho_w 1000.0 Rv 461.495 k2_prime 22.1 k3 373900.0 Pi 1e6 / (rho_w * Rv * (k3 / Tm k2_prime)) # 5. 可降水量 PWV单位 mm df[PWV] df[ZWD] * Pi * 1000 # ZWD单位m乘1000转为mm # 6. 质量控制剔除负值和不合理的大值 df.loc[df[PWV] 0, PWV] np.nan df.loc[df[PWV] 100, PWV] np.nan df.to_csv(BJFS_PWV_5min.csv, indexFalse)这段脚本对每个历元的ZTD做一个固定的干延迟扣除。严格来说气压和温度随时间变化ZHD也应该是时变的但逐时气压变化引起的ZHD变化一般只有几毫米折合PWV误差约1毫米。如果你的精度目标高于1毫米建议用测站逐小时气象数据或者ERA5气压产品按时间逐历元做干延迟扣除脚本逻辑是一样的只是把ZHD从标量改成向量。5. 精度评估与探空、ERA5对比的核心指标与结果解读5.1 参考数据的获取与时间匹配PWV反演出来后引用审稿人或者导师最喜欢问的一句话就是“精度怎么样”。要回答这个问题必须有第三方参考数据。常见的参考数据源是探空资料和ERA5再分析产品。探空资料推荐用IGRA2数据集它整合了全球探空站自1940年代以来的数据温度和湿度剖面质量经过严格检验。从IGRA2计算PWV的做法是从地面开始沿气压层积分湿度廓线公式为PWV -(1/g) ∫ q dp实际用数值积分实现。这里容易踩的一个坑是探空仪的湿度传感器存在“迟滞效应”尤其在湿度剧烈变化时上升过程的高层湿度往往偏干导致探空PWV在某些天气系统过境时偏低这是系统性误差对比时要有心理预期。ERA5数据可以从CDSCopernicus Climate Data Store下载单层变量里面有一个直接的变量叫做total_column_water_vapour单位是kg/m²数值上就等于PWV的毫米数拿过来直接用就行不需要自己积分。如果你下载的是气压层变量也可以按比湿和气压计算PWV但没必要直接拉这个变量更省事。时间匹配上要注意探空资料一天只有0点和12点两个时次GNSS PWV要取相同时刻前后30分钟的平均值再和探空比较不能直接把历元值和探空值硬套。ERA5是逐小时输出可以比较GNSS PWV与整点时刻的值匹配容易得多。5.2 评估指标与Python计算精度评估通常报告四个指标偏差Bias、标准差STD、均方根误差RMSE和相关系数R。计算公式如下Bias Σ(GNSS - REF) / nSTD sqrt( Σ(Δ - Bias)² / (n - 1) )RMSE sqrt( Σ(GNSS - REF)² / n )R cov(GNSS, REF) / (std(GNSS) × std(REF))计算脚本也不复杂import numpy as np import pandas as pd from scipy import stats # 读入GNSS PWV和参考PWV且已经完成时间匹配 gnss pd.read_csv(BJFS_PWV.csv, parse_dates[time]) ref pd.read_csv(BJFS_ref.csv, parse_dates[time]) merged pd.merge(gnss, ref, ontime, howinner) valid merged.dropna(subset[PWV_gnss, PWV_ref]) d valid[PWV_gnss] - valid[PWV_ref] bias np.mean(d) std np.std(d, ddof1) rmse np.sqrt(np.mean(d ** 2)) r, p stats.pearsonr(valid[PWV_gnss], valid[PWV_ref]) print(fBias {bias:.2f} mm, STD {std:.2f} mm, RMSE {rmse:.2f} mm, R {r:.3f})实际处理中我还会额外画两幅图一个是GNSS PWV与参考PWV的时间序列对比曲线另一个是散点密度图加1:1线。散点图能直观显示偏差的方向性比如如果散点整体偏到1:1线下方说明GNSS PWV系统性偏低就要回头检查气压输入或者Tm是否偏大。5.3 结果解读与常见偏差来源从大量已发表研究和我的实测经验来看GNSS PWV与探空的RMSE通常在2到4毫米与ERA5的RMSE在2到5毫米这个区间Bias一般能控制在±1毫米以内相关系数在0.95以上。如果结果明显差于这个范围优先排查以下几个问题。气压误差是最大嫌疑。前面说过1百帕的气压误差直接带来约2.3毫米的ZHD误差等价于大约0.3毫米的PWV误差。检查一下你用的气压是测站实测还是模型内插内插值的准确性会直接影响干延迟扣除。Tm误差排在第二。很多人在热带或者高海拔地区直接用Bevis公式这种地区Tm与探空实际值的偏差可能达到2到4开对应PWV相对误差2%到4%。如果发现偏差随季节呈正弦变化基本就是Tm模型没选对。还有一个容易被忽略的是测站高程系统的不一致。PRIDE PPP-AR解算时使用的测站坐标如果是CGCS2000或ITRF框架下的椭球高那Saastamoinen公式里H用的就是椭球高但如果测站资料里混入了正常高ZHD计算就会差出十几米高程对应的气压差这在山区站点尤其明显。每次处理新测站时先确认高程是椭球高还是正高两者差一个大地水准面差距千万别搞混。6. 实测下来最容易被忽视的六个细节第一时间系统必须统一。PRIDE PPP-AR输出文件里通常用的是GPS时或者UTC但探空资料、ERA5用的是世界时UTC如果两边时间标签差出十几秒甚至几分钟短时强降水过程中的水汽快速变化特征的对比会被明显削弱。数据预处理阶段要先把时间系统转成同一个我一般在脚本里统一转成UTC以GNSS输出的时标为基准核对。第二ZTD解算的采样间隔不要盲目追求高分辨率。5分钟ZTD看起来很诱人但每次估计的观测冗余少噪声大最后算出的PWV曲线毛刺会很明显。如果是业务预报需求10到15分钟已经很够用研究用的历史序列我通常生成30分钟产品稳。第三DCB文件不要漏但也别用错版本。用CODE发布的DCB文件时要确保DCB对应的卫星、频率和你的解算策略一致。我遇到过一次把GPS和BDS的DCB搞混结果BDS模糊度固定成功率暴跌ZTD序列出现系统性偏差排查了很久才定位到问题。第四天线相位中心文件igs20.atx要用新版本。随着IGS框架更新天线PCO/PCV参数也在更新老版本的atx会让坐标解算有毫米到厘米级误差但ZTD受影响相对小。不过为了严谨建议每半年检查一次IGS发布的新版本。第五ZTD质量控制要放在PWV反演之前不要放到PWV之后。先根据ZTD连续性和解算标志剔除坏历元再做干延迟扣除和PWV计算。如果先算PWV再剔除坏值会通过插值污染相邻历元尤其是靠近坏值的PWV值会在时间序列图上出现明显缺口或尖刺。第六与探空做精度评估时要把探空气球漂移纳入考量。探空气球不是垂直上升的释放后会在高空风作用下水平漂移几十公里PWV对比隐含了空间不一致性。对于天气系统尺度小、地形复杂的区域这种漂移导致的代表性误差可能到1到2毫米解释结果时要留有余地。做GNSS水汽反演这件事核心不是把软件跑通而是每一步都知道自己在算什么、误差从哪里来。PRIDE PPP-AR把ZTD解算的门槛降下来了但精度上限还是由你对数据处理细节的理解程度决定的。如果这篇文章能让你少走一些弯路少浪费几周调试时间目的就达到了。
PRIDE PPP-AR实现GNSS水汽反演全流程：从ZTD解算到PWV精度评估
搞大气水汽反演的人应该都有过这种纠结探空站一天就放两次气球时间分辨率低到让人抓狂再分析资料空间分辨率虽然过得去但时效性和精度没法保证。GNSS连续运行参考站网把这个缺口补上了大半——卫星信号穿过大气层时受到的对流层延迟这个定位中的“噪声”被当成信号来用反演出的GNSS PWV时间分辨率能做到5分钟甚至更高。而要把这件事做准天顶总延迟ZTD的解算精度是第一道门槛。我这两年一直在用武汉大学PNT中心开源的PRIDE PPP-AR来处理ZTD解算和PWV反演踩过不少坑也积累了一套比较顺手的全流程整理出来给同行参考。这篇文章不会只贴命令会把参数背后的逻辑、输出的含义、精度评估里容易出问题的地方都说清楚适合刚接触GNSS气象学、准备用PPP-AR做水汽时间序列的研究生和工程师。1. 为什么是PRIDE PPP-AR这条技术链路的价值与选型逻辑1.1 GNSS水汽反演到底解决了什么痛点水汽是大气中最活跃的成分虽然只占大气总量的0.1%到3%但它的时空变化直接决定了强对流天气的发生发展。传统观测手段里探空站一天只放两次气球时间分辨率12小时遇到天气系统快速演化的过程基本抓不住微波辐射计精度尚可但设备贵、站点稀少再分析资料像ERA5虽然时空分辨率不错但它是模式输出结果时空分辨率受限于同化系统和观测密度在局地对流过程面前仍然会“钝化”。地基GNSS反演PWV的出现相当于给水汽观测网格里补上了高时间分辨率的一环。只要有一个连续运行的GNSS测站卫星信号穿过大气层时产生的天顶总延迟ZTD就可以按分钟级解算出来从中剥离出湿延迟再乘上一个转换系数就能得到大气可降水量PWV。这个值与探空、微波辐射计的结果有很好的一致性而且全天候、不受云雨影响、设备维护成本低。现在很多气象业务单位把GNSS PWV纳入短临预报的参考因子尤其是强降水过程前水汽快速增加的特征GNSS时间序列能给出非常有价值的信息。1.2 为什么选PRIDE PPP-AR而不是GAMIT和Bernese做ZTD解算的软件不少常见的路线大概有三类。我只说我实际用下来的感受不做盲目吹捧。GAMIT/GLOBK是相对定位模式需要多站联测做网解解算前要挑选基准站、做基线处理、处理周期模糊度整个流程比较重。如果只是想要单站或者十几个站的水汽时间序列用GAMIT会有“杀鸡用牛刀”的感觉而且计算量大、学习曲线陡。Bernese是功能最全面的学术软件之一从PPP到网解、从对流层到电离层都能做但操作门槛更高命令行体系、脚本配置、文件组织方式都很“学院派”新手从零到跑通第一个解算流程通常要花一到两周。PRIDE PPP-AR的特点在于它走的是精密单点定位加模糊度固定路线单站就能解算不需要联测基准站支持GPS、GLONASS、Galileo、BDS多系统联合解算数据利用率高编译和使用相对轻量批量处理测站很方便最关键的是它对ZTD这类对流层参数的估计精度与GAMIT和Bernese的结果在毫米级能够保持一致。对于做水汽研究的人来说这意味着不用维护一套复杂的网解环境也能拿到高质量的对流层延迟产品。另外这几年PRIDE PPP-AR一直在更新支持RINEX 3.x格式、支持多系统多频率模糊度固定对北斗数据也有比较好的处理能力。国内做GNSS气象应用的研究组用得越来越多社区里的经验积累也相对多一些遇到问题容易找到人讨论。1.3 整条技术链路长什么样从原始数据到最终的PWV精度评估完整链路可以拆成四个环节数据准备GNSS观测文件RINEX格式、精密星历、精密钟差、DCB文件、天线相位中心文件。ZTD解算用PRIDE PPP-AR做精密单点定位解算逐历元输出天顶总延迟ZTD。PWV反演从ZTD中减去干延迟ZHD得到湿延迟ZWD再通过加权平均温度Tm和转换系数Π把ZWD换算成PWV。精度评估把GNSS PWV与探空资料或ERA5再分析资料的PWV做对比计算偏差、标准差、均方根误差和相关系数。这个链路每一步都有坑数据下载缺文件、配置参数选错、干延迟扣除精度不足、Tm模型偏差过大、时间系统没统一等等。后面我会一个一个展开说把我实际遇到过的情况和解决办法都写出来。2. 环境与数据从源码编译到观测文件准备2.1 编译前的系统环境配置PRIDE PPP-AR目前主要面向Linux环境Windows下虽然可以用WSL或者虚拟机跑但原生Linux还是最省心的。我用的是Ubuntu 20.04长期支持版内核和编译工具链都比较稳定。新装完系统后先装好基础依赖再编译不然会在configure阶段反复报错。sudo apt update sudo apt install build-essential gfortran make wget curl sudo apt install libcurl4-openssl-dev libxml2-dev这里有几个容易忽略的点。第一gfortran版本不要太老也不要太激进Ubuntu 20.04自带的gfortran 9处理PRIDE PPP-AR的代码没问题但如果你用的是CentOS 7自带的老版本gfortran某些新语法会编译失败。第二如果源码里使用了网络下载模块libcurl开发库是必须的缺了这个头文件会在编译中报curl/curl.h: No such file or directory这种错误。第三建议装一个nco和cdo工具包后面处理ERA5 NetCDF数据会用到虽然不属于PRIDE PPP-AR的依赖但在精度评估阶段几乎是刚需。2.2 获取源码与编译PRIDE PPP-AR的源码托管在GitHub上搜索PRIDE-PPP-AR就能找到也可以从武汉大学PRIDE实验室官网的链接跳转。下载后进入源码目录按标准流程编译git clone https://github.com/PRIDE-Lab/PRIDE-PPP-AR.git cd PRIDE-PPP-AR ./configure --prefix$(pwd)/install make -j4 make install--prefix参数指定安装路径这是我一直建议的做法把编译产物和源码分开后续更新版本时不用重新编译旧目录只需要切换路径就行。-j4表示四核并行编译如果是八核机器可以改成-j8能省不少时间。编译成功的标志是src目录下生成了主程序pridepppar或者你设置了--prefix后在install/bin下能看到它。如果遇到gfortran: error: unrecognized command-line option之类报错基本可以判定是编译器版本问题换gfortran 9或10重新编译即可。一个值得注意的经验PRIDE PPP-AR的代码对编译选项比较敏感如果configure阶段出现奇怪的错误建议先make distclean清理干净再来一遍不要硬着头皮在报错状态下反复修改makefile。2.3 数据准备观测文件、精密星历与钟差下载ZTD解算需要三类核心输入测站观测文件、精密星历、精密钟差。观测文件从IGS数据中心下载CDDIS是常用的来源但国内访问有时不太顺畅中国这边的替代方案是武汉大学IGS数据中心igs.gnsswhu.cn下载速度明显更好。RINEX文件的命名遵循统一规则例如BJFS00CHN_R_20240010000_01D_30S_MO.rnx表示2024年第1天北京房山站BJFS的30秒采样、24小时观测文件。精密星历和钟差产品按时效分为最终、快速、超快速三类。做研究用最终产品精度最高但发布时间有一到两周延迟近实时应用用快速产品钟差精度在亚纳秒量级实时应用用超快速预报产品ZTD精度会差一些。IGS的最终产品命名类似igs21927.sp3和igs21927.clk前者是SP3格式的精密轨道后者是精密钟差文件这个21927对应GPS周和星期几解算前要确保把当天对应的所有历元覆盖。DCB文件差分码偏差也是容易漏掉的一项。多系统PPP中不同频率、不同系统间的码偏差如果不改正会影响模糊度固定和时间参数估计。CODE分析中心发布DCB产品格式是CODE20240010000.DCB这种命名。同时还需要天线相位中心文件igs20.atx这些文件可以和精密产品一起放在一个公共目录里各测站解算时共用。我建议在本地建一个清晰的数据目录结构至少包括obs/、orbit/、clock/、dcb/、atx/这几个子目录名称固定下来不要变。PRIDE PPP-AR运行时会按配置文件里的路径去读文件目录结构越规律写批处理脚本越省事。3. ZTD解算实操配置参数背后的原理与输出文件解读3.1 解算前先理解ZTD是怎么被估计出来的在PPP解算里ZTD不是直接观测量而是作为未知参数参与平差。卫星信号穿过对流层时产生的延迟大致分成两部分静力学延迟干延迟ZHD和湿延迟ZWD。干延迟变化比较有规律可以先通过模型扣除一部分湿延迟变化快、与空间位置关系密切是残余估计的主角。PRIDE PPP-AR的处理逻辑是先用Saastamoinen模型给出一个先验ZTD再用VMF1或GMF映射函数把天顶方向的延迟映射到卫星信号的实际传播路径上路径上的误差被分解成水平梯度和残余天顶延迟两部分。残余天顶延迟按随机游走过程进行逐历元估计后面再加上模糊度固定把载波相位模糊度固定到整数从而显著提高ZTD估计的稳定性和精度。这就解释了为什么模糊度固定对ZTD这么重要不固定模糊度时模糊度参数和天顶延迟参数之间的相关性会通过观测方程传染ZTD序列容易出现小幅漂移模糊度固定后参数之间的耦合被切断ZTD时序就变得干净很多。所以用PRIDE PPP-AR时我的建议是务必开启模糊度固定功能否则软件最大的优势没有发挥出来。映射函数的选择也要留意。VMF1和GMF各有优劣VMF1使用的是ECMWF数值天气预报产品计算的对流层映射函数精度略好但需要外部下载VMF1格网产品GMF是经验模型不依赖外部数据穿透性和稳定性好。如果你做的是高精度水汽研究我推荐VMF1如果只是批量处理常规站GMF完全够用。3.2 配置文件逐项解读PRIDE PPP-AR通过一个文本配置文件控制解算行为不同版本的配置格式略有差异但核心参数是相通的。以下是我自己常用的一套配置内容标注了每个参数的实际意义# 测站与数据路径 site_name BJFS obs_file ./obs/BJFS00CHN_R_20240010000_01D_30S_MO.rnx nav_file ./obs/BJFS00CHN_R_20240010000_01D_MN.rnx sp3_file ./orbit/igs21927.sp3 clk_file ./clock/igs21927.clk dcb_file ./dcb/CODE20240010000.DCB atx_file ./atx/igs20.atx # 解算策略 mode ppp-ar interval 30 # 观测值采样间隔单位秒 elevation_mask 7 # 截止高度角单位度 system grce # 参与解算的系统gps/glonass/galileo/bds # 对流层估计 trop_model saastamoinen mapping_function vmf1 trop_interval 30 # 天顶延迟估计间隔单位分钟 gradient_est yes # 估计水平梯度interval 30对应观测文件里数据本身的采样率如果你的RINEX文件是30秒采样这里就必须写30不能随意改。elevation_mask 7是我批量处理时比较折中的方案低于7度的观测值受多路径和大气折射影响较大太高了又会损失低仰角卫星带来的水汽信息。低仰角卫星信号穿过大气的路径更长携带的水汽信息也更多对ZTD估计是有利的但要付出噪声增大的代价。gradient_est yes建议一直打开。天顶延迟是一维的但大气水汽在水平方向往往不均匀尤其是遇到锋面过境时东西方向和南北方向的湿延迟差异很大。如果忽略水平梯度天顶延迟估计会被污染表现就是解出的ZTD序列在天气系统过境时出现系统性偏差。trop_interval 30表示每30分钟估计一个天顶延迟参数。这个参数决定ZTD时间分辨率想要更高分辨率的PWV可以改成15甚至5但要注意时间分辨率越高单个历元的观测冗余越少ZTD噪声会增大。做常规水汽研究15到30分钟是平衡点。3.3 运行解算与输出文件解析配置写好后运行主程序解算单个测站./pridepppar config_BJFS.in程序运行时会打印每个历元的处理进度同时输出一个后缀为.ztd或类似命名的文件里面按时间顺序列出每个历元的ZTD估计值。典型输出列的格式大致是年、年积日、秒、天顶总延迟ZTD单位米、估计标准差、水平和垂直方向的梯度参数等。拿到输出后第一件事不是直接算PWV而是看ZTD时间序列是否平滑连续、有没有跳变。正常情况下ZTD是缓慢变化的大气信号数值一般在1.5米到2.7米之间取决于海拔和季节如果出现单历元突然跳变几厘米、甚至到负值或者0的情况多半是观测文件中存在粗差或者模糊度固定失败需要把这一段标记为无效不要进入后续反演。我习惯在批量处理完后写一个质量控制脚本把ZTD序列中偏离前后中位数超过3倍标准差的历元剔除掉同时删除解算状态标志不是正常收敛的历元。这一步看似简单但它决定了PWV时间序列的“干净程度”不能省。4. PWV反演从湿延迟到可降水量的完整脚本4.1 干延迟与湿延迟分离ZTD解算出来后得到的是天顶方向的总延迟要得到湿延迟ZWD需要先扣除干延迟ZHD。干延迟可以用地面气压和测站位置计算最常用的是Saastamoinen模型ZHD 0.0022768 × Ps / (1 - 0.00266 × cos(2φ) - 0.00028 × H)其中Ps是测站处地面气压单位hPaφ是测站纬度H是测站椭球高单位km。从公式可以看出干延迟与气压几乎成正比气压1百帕的误差会引起大约2.3毫米的干延迟误差进而影响同量级的PWV误差所以地面气压的精度很关键。最佳做法是在测站附近布设气象传感器实时获取气压、温度、湿度。但现实往往是测站没有气象观测设备这种情况下可以用数值预报模型的气压产品。欧洲中心ERA5的逐小时2米气压或者GPT2w、GPT3经验模型内插的气压值都能用。我自己实测下来GPT3模型给出的气压与真实气压差通常在1到2百帕以内换算成PWV误差大约在0.5到1毫米对大多数研究场景可以接受。4.2 加权平均温度Tm的获取Bevis公式与GPT2w/GPT3从湿延迟换算可降水量核心转换系数里有一个变量叫加权平均温度Tm它的定义是大气柱内水汽压和温度的加权平均。Tm不能直接观测但可以经验估算。最经典的方法是美国学者Bevis提出的线性回归公式Tm 70.2 0.72 × TsTs是测站处地表温度单位K这个公式是从美国地区探空资料拟合出来的在北半球中纬度地区表现很好。但如果你做的区域是热带、高原或者高纬度地区直接用Bevis公式可能带来1到3开的偏差而Tm每偏差1开PWV大约有1%的相对误差对于20毫米的PWV就是0.2毫米做高精度研究时不可忽视。更好的方案是用GPT2w或GPT3经验模型它们按经纬度和高程给出Tm的长期平均值和周期变化项输入测站坐标和年积日就能插值出来空间代表性比单点的Bevis公式更好。用GPT系列模型不需要外部实时数据而且它在全球范围都做过验证。我个人的处理流程是优先用测站气象观测数据按Bevis公式计算Tm没有气象数据就用GPT3模型获取表面温度和Tm如果做批量历史数据处理也可以直接插值ERA5气压层温度和水汽剖面来计算Tm精度更高但计算量也更大。4.3 PWV计算与Python脚本实现转换系数Π的公式如下Π 10^6 / (ρ_w × R_v × (k3 / Tm k2))其中ρ_w是液态水密度取1000 kg/m³R_v是水汽气体常数取461.495 J/(kg·K)k2和k3是大气折射率常数常用的一组取值是k2 22.1 K/hPak3 373900 K²/hPa。注意不同文献里常数取值略有差异选一组并保持全程一致就好换成不同常数组会让PWV产生约0.5%到1%的差异。PWV Π × ZWD下面是完整的Python计算脚本注释已经写得比较详细直接抄作业即可import numpy as np import pandas as pd # 读取ZTD解算结果 # 假设文件列为year doy sod ztd ztd_std ... df pd.read_csv(BJFS.ztd, delim_whitespaceTrue, names[year, doy, sod, ztd, std, grad_n, grad_e], skiprows1) # 测站参数 lat 39.6 # 纬度单位度 h_ell 80.0 # 椭球高单位 km北斗站约80m 0.08km P_s 1013.25 # 地面气压单位 hPa来自气象数据或GPT2w T_s 288.15 # 地表温度单位 K # 1. 干延迟 ZHDSaastamoinen 模型 phi np.deg2rad(lat) ZHD 0.0022768 * P_s / (1 - 0.00266 * np.cos(2 * phi) - 0.00028 * h_ell) # 2. 湿延迟 ZWD df[ZWD] df[ztd] - ZHD # 3. 加权平均温度 Tm这里用Bevis公式可替换为GPT2w插值结果 Tm 70.2 0.72 * T_s # 单位 K # 4. 转换系数 Π rho_w 1000.0 Rv 461.495 k2_prime 22.1 k3 373900.0 Pi 1e6 / (rho_w * Rv * (k3 / Tm k2_prime)) # 5. 可降水量 PWV单位 mm df[PWV] df[ZWD] * Pi * 1000 # ZWD单位m乘1000转为mm # 6. 质量控制剔除负值和不合理的大值 df.loc[df[PWV] 0, PWV] np.nan df.loc[df[PWV] 100, PWV] np.nan df.to_csv(BJFS_PWV_5min.csv, indexFalse)这段脚本对每个历元的ZTD做一个固定的干延迟扣除。严格来说气压和温度随时间变化ZHD也应该是时变的但逐时气压变化引起的ZHD变化一般只有几毫米折合PWV误差约1毫米。如果你的精度目标高于1毫米建议用测站逐小时气象数据或者ERA5气压产品按时间逐历元做干延迟扣除脚本逻辑是一样的只是把ZHD从标量改成向量。5. 精度评估与探空、ERA5对比的核心指标与结果解读5.1 参考数据的获取与时间匹配PWV反演出来后引用审稿人或者导师最喜欢问的一句话就是“精度怎么样”。要回答这个问题必须有第三方参考数据。常见的参考数据源是探空资料和ERA5再分析产品。探空资料推荐用IGRA2数据集它整合了全球探空站自1940年代以来的数据温度和湿度剖面质量经过严格检验。从IGRA2计算PWV的做法是从地面开始沿气压层积分湿度廓线公式为PWV -(1/g) ∫ q dp实际用数值积分实现。这里容易踩的一个坑是探空仪的湿度传感器存在“迟滞效应”尤其在湿度剧烈变化时上升过程的高层湿度往往偏干导致探空PWV在某些天气系统过境时偏低这是系统性误差对比时要有心理预期。ERA5数据可以从CDSCopernicus Climate Data Store下载单层变量里面有一个直接的变量叫做total_column_water_vapour单位是kg/m²数值上就等于PWV的毫米数拿过来直接用就行不需要自己积分。如果你下载的是气压层变量也可以按比湿和气压计算PWV但没必要直接拉这个变量更省事。时间匹配上要注意探空资料一天只有0点和12点两个时次GNSS PWV要取相同时刻前后30分钟的平均值再和探空比较不能直接把历元值和探空值硬套。ERA5是逐小时输出可以比较GNSS PWV与整点时刻的值匹配容易得多。5.2 评估指标与Python计算精度评估通常报告四个指标偏差Bias、标准差STD、均方根误差RMSE和相关系数R。计算公式如下Bias Σ(GNSS - REF) / nSTD sqrt( Σ(Δ - Bias)² / (n - 1) )RMSE sqrt( Σ(GNSS - REF)² / n )R cov(GNSS, REF) / (std(GNSS) × std(REF))计算脚本也不复杂import numpy as np import pandas as pd from scipy import stats # 读入GNSS PWV和参考PWV且已经完成时间匹配 gnss pd.read_csv(BJFS_PWV.csv, parse_dates[time]) ref pd.read_csv(BJFS_ref.csv, parse_dates[time]) merged pd.merge(gnss, ref, ontime, howinner) valid merged.dropna(subset[PWV_gnss, PWV_ref]) d valid[PWV_gnss] - valid[PWV_ref] bias np.mean(d) std np.std(d, ddof1) rmse np.sqrt(np.mean(d ** 2)) r, p stats.pearsonr(valid[PWV_gnss], valid[PWV_ref]) print(fBias {bias:.2f} mm, STD {std:.2f} mm, RMSE {rmse:.2f} mm, R {r:.3f})实际处理中我还会额外画两幅图一个是GNSS PWV与参考PWV的时间序列对比曲线另一个是散点密度图加1:1线。散点图能直观显示偏差的方向性比如如果散点整体偏到1:1线下方说明GNSS PWV系统性偏低就要回头检查气压输入或者Tm是否偏大。5.3 结果解读与常见偏差来源从大量已发表研究和我的实测经验来看GNSS PWV与探空的RMSE通常在2到4毫米与ERA5的RMSE在2到5毫米这个区间Bias一般能控制在±1毫米以内相关系数在0.95以上。如果结果明显差于这个范围优先排查以下几个问题。气压误差是最大嫌疑。前面说过1百帕的气压误差直接带来约2.3毫米的ZHD误差等价于大约0.3毫米的PWV误差。检查一下你用的气压是测站实测还是模型内插内插值的准确性会直接影响干延迟扣除。Tm误差排在第二。很多人在热带或者高海拔地区直接用Bevis公式这种地区Tm与探空实际值的偏差可能达到2到4开对应PWV相对误差2%到4%。如果发现偏差随季节呈正弦变化基本就是Tm模型没选对。还有一个容易被忽略的是测站高程系统的不一致。PRIDE PPP-AR解算时使用的测站坐标如果是CGCS2000或ITRF框架下的椭球高那Saastamoinen公式里H用的就是椭球高但如果测站资料里混入了正常高ZHD计算就会差出十几米高程对应的气压差这在山区站点尤其明显。每次处理新测站时先确认高程是椭球高还是正高两者差一个大地水准面差距千万别搞混。6. 实测下来最容易被忽视的六个细节第一时间系统必须统一。PRIDE PPP-AR输出文件里通常用的是GPS时或者UTC但探空资料、ERA5用的是世界时UTC如果两边时间标签差出十几秒甚至几分钟短时强降水过程中的水汽快速变化特征的对比会被明显削弱。数据预处理阶段要先把时间系统转成同一个我一般在脚本里统一转成UTC以GNSS输出的时标为基准核对。第二ZTD解算的采样间隔不要盲目追求高分辨率。5分钟ZTD看起来很诱人但每次估计的观测冗余少噪声大最后算出的PWV曲线毛刺会很明显。如果是业务预报需求10到15分钟已经很够用研究用的历史序列我通常生成30分钟产品稳。第三DCB文件不要漏但也别用错版本。用CODE发布的DCB文件时要确保DCB对应的卫星、频率和你的解算策略一致。我遇到过一次把GPS和BDS的DCB搞混结果BDS模糊度固定成功率暴跌ZTD序列出现系统性偏差排查了很久才定位到问题。第四天线相位中心文件igs20.atx要用新版本。随着IGS框架更新天线PCO/PCV参数也在更新老版本的atx会让坐标解算有毫米到厘米级误差但ZTD受影响相对小。不过为了严谨建议每半年检查一次IGS发布的新版本。第五ZTD质量控制要放在PWV反演之前不要放到PWV之后。先根据ZTD连续性和解算标志剔除坏历元再做干延迟扣除和PWV计算。如果先算PWV再剔除坏值会通过插值污染相邻历元尤其是靠近坏值的PWV值会在时间序列图上出现明显缺口或尖刺。第六与探空做精度评估时要把探空气球漂移纳入考量。探空气球不是垂直上升的释放后会在高空风作用下水平漂移几十公里PWV对比隐含了空间不一致性。对于天气系统尺度小、地形复杂的区域这种漂移导致的代表性误差可能到1到2毫米解释结果时要留有余地。做GNSS水汽反演这件事核心不是把软件跑通而是每一步都知道自己在算什么、误差从哪里来。PRIDE PPP-AR把ZTD解算的门槛降下来了但精度上限还是由你对数据处理细节的理解程度决定的。如果这篇文章能让你少走一些弯路少浪费几周调试时间目的就达到了。