基于GAMIT系统的GPS/MET数据反演及其在闽南地区获取PWV的应用

杜言霞 ，梁 莺 ，吴勇凯 ，陈州川 ，林添水

（1.泉州市气象局，福建 泉州362000；2.海峡气象开放研究实验室，福建 厦门361012；3.福建省大气探测技术保障中心，福建 福州350008；4.安溪县气象局，福建 安溪362400）

随着全球导航卫星和低轨卫星技术的进步，GPS气象学也得到了迅速发展。GPS技术具有自动、精确、快速、成本低、时空分辨率高和不易受天气变化影响等优势，是一种极具潜力的探测手段[1]，因此得到很多国家和地区的关注。研究表明，利用地基GPS观测数据和气象数据反演出的大气可降水量具有较高的精度，在暴雨预测、数值天气预报和其它极端天气预测中都起到了重要的作用。

GPS气象学（GPS/MET）分为地基GPS气象学和空基GPS气象学。20世纪80年代后期美国的Askne和Nordius最早提出用地基GPS网探测地球大气，他们推导了PWV和大气湿延迟之间的关系[2]。Bevis和Businger进一步分析和讨论了这个转换关系，从而提出地基GPS气象学[3]。

GPS/MET作为一种大气遥感探测新技术，国外发达国家在技术研究及业务应用上，已经取得显著成果，世界上许多国家都建立了GPS观测网。最早是美国在1996年就开始组建覆盖全国的GPS/MET网，德国在地基和空基GPS气象学方面也做了大量的工作，德国气象局和地学研究中心（CFZ&DWD）利用组建的GPS站网数据，每小时提供一次全国GPS/PWV分布图。近年来，我国气象部门也积极开展GPS/MET的研究和业务化试验，取得了一定的成绩，把建设地基GPS探测大气和电离层参数的系统列入了现代化气象探测系统的一项重要内容，并制定了“中国气象局导航卫星接收站网建设规划”。截至目前，我国已组建覆盖全国的实时GPS/MET观测网，提供1 h间隔的全国范围的可降水量图。

1 GAMIT软件简介

GAMIT软件系统是由美国麻省理工学院（MIT）开发的，由于它具有开源性以及版本更新速度快等诸多优点，因此得到了较为广泛的应用。它主要由以下几个程序模块构成：单差自动修复周跳模块、双差自动修复周跳模块、数据文件的格式转换、轨道积分模块、观测值模型、人工交互式修复周跳模块、参数估计、用于创建SOLVE所需的M文件和一些辅助程序。并且在数据编辑完成之后，可利用SOLVE模块求解出相关的未知参数[4]。

GAMIT软件的主要功能包括：估算卫星状态矢量、太阳光压参数、卫星Y方向偏差、地球自转参数、台站坐标和对流层天顶延迟等[5]。

GAMIT采用双差做基本观测量，可以消除站钟和星钟之间的主要误差；轨道误差对测站相对位置不敏感，有利于精密定位；对于台站初始坐标的精度要求也不高[6]。但是减少了有效观测，各站有效观测因其在网内位置的变化而变化，且不能用于单点定位。

2 GAMIT数据处理过程

GAMIT软件数据处理过程分为3个步骤：GPS数据准备、GAMIT参数配置和GAMIT数据解算。

2.1 GPS数据准备

数据处理过程主要用到3种数据文件，分别为观测文件、星历文件和导航文件，这些文件可以从IGS 网站（ftp：//garner/ucsd.edu/pub）下载，但下载后需转换为RINEX格式的文件。若下载的数据文件在WINDOWS系统下进行处理时，可通过dos2unix命令进行格式转换[7]。

在单天目录里用sh_link_rinex命令进行相应文件的链接，也可以直接把原始观测文件、星历文件和导航文件拷贝到相应的单天目录。然后进入tables目录，建立与/gamit/tables目录的连接（%links.table J2000 2001）。

GPS卫星星历的质量对GPS卫星定位的精度起着至关重要的作用，利用GAMIT软件解算GPS数据，对于不同的研究采用不同的星历。如果反演的大气可降水量用于天气预报，那么应采用精密预报星历进行解算，如果用于长期性的气象及气候学研究，则可采用各种精密星历（表1）。

表1 星历一览表

2.2 软件参数配置

为使GAMIT的内置参数能与GPS数据匹配，需重新配置 sittbl.、sestbl.、station.info、session.info 等文件，配置完成后一起放入表目录[8]。

（1）台站坐标文件。获得测站概略坐标的方法主要有：利用原始观测RINEX文件获得、利用命令gapr_to_l生成L文件、利用命令sh_rx2apr生成单测站的L文件，得到的L文件需进行手工编辑。

在GAMIT运行过程中，初始坐标的误差必须在10 m以内，否则可能会导致最后迭代不收敛。初始坐标的误差在1 m以内，才能保证最终的台站坐标精度好于1 mm。但有时因为接收机性能不稳定，或者是周围观测条件不佳，会影响某些时段的观测，导致单点定位精度不好或者定位时不收敛，那么需要重新选择其他时段的RINEX观测文件。为了获得较精准的初始坐标，通常可以多选取几天的RINEX观测文件，进行单点定位，对多次定位结果取均值。

（2）台站信息文件。station.info文件包含了所有需要解算的台站信息，包括台站的天线相位中心对于地球标准椭球面在X、Y、Z方向的偏差以及接收机型号、天线型号等信息。这些信息都可从每个台站的观测文件的头文件中获取，其中所谓天线相位中心在X、Y、Z方向的偏差，实际只有天线相位中心的高度，因为天线相位中心在Y和X方向的偏差一般为零。可以参考/templates/station.info对文件按照规定的格式进行手工编辑，并根据研究需要引入相应测站。另外，工程文件的名字也在这里定义。

值得注意的是，有些台站的信息因为更换接收机或迁址等原因，会进行信息更新，在station.info文件中同一个台站会有几行信息，必须按时间先后顺序排列，将最新的信息放在最后一行，GAMIT默认读取最后一行信息。

（3）运行控制文件。sestbl.文件是数据处理方案的核心控制文件，很多参数都需要在该文件里配置。具体要根据GPS网类型、所要解决的问题等采取相应的配置方案。主要对sestbl.文件中的几项进行设置：①choice of experiment。RELAX：定轨、定位、解ERP；BASELINE：仅仅定位；ORBIT：仅仅定轨。一般为解算水汽获取大气总延迟的时候，仅仅需要选择BASELINE，不需要对卫星进行定轨，广播星历的精度已经足够。

②choice of observable。LC_HELP：用 LC 观测解模糊度；LC_RANGE：用于LC_HELP，但更强调伪距的作用；LC_ONLY：用 LC观测不解模糊度；L1_ONLY：仅使用L1解模糊度；L2_ONLY：仅使用L2 解模糊度；L1，L2_INDEPEND：分别使用 L1，L2解模糊度。

③type of analysis。QUICK：每一个缺口加周跳标志，对每一个周跳标志（在连续资料中或在缺口处）求解bias，这样解得参数多了，解的精度也就低了，但对周跳不敏感；0-ITER：ARC（optional），MODEL，AUTCLN （optional），SOLVE，SCANDD；1-ITER：两个0-ITER序列，但第一个是QUICK解；2-ITER：在1-ITER基础上再加一序列，用于确定轨道；对于仅仅解算基线，获取大气延迟量，只需要选取0-ITER即可。

（4）台站约束文件。通过设置sittbl.文件可以完成测站约束，另外大气模型等参数也在这个文件里设置。首先要检查是否包含所有需要解算的台站的信息，设置的约束值的大小要依据L文件的精度，如果L文件中初始台站坐标比较准确，约束值可以设置比较小；如果L文件误差大，约束值就要比较大，一般为5~10 m，否则在最后解算时会导致迭代不收敛。

（5）时段信息文件。可以手工编辑session.info文件，也可以用sh_makexp命令创建，采样间隔和卫星号等参数都在这个文件里设置。

2.3 数据处理

GAMIT数据处理过程主要分为3个阶段：首先利用GAMIT软件处理下载的观测数据，然后利用卡尔曼滤波进行综合解算，最后求出对流层天顶延迟[9]。

用户可以根据研究情况需要，对约束条件和参数进行设置。卫星截止高度角、数据采样率、星历等参数也可以根据研究需要进行选择。GPS数据处理流程如图1所示。

首先，利用GAMIT软件解算出载波相位的理论值，然后根据理论值和观测值求出相位残差，最终求出坐标以及各个大气参数等。

以外文期刊为例，尝试构建一个基于学科与用户行为的外文期刊的资源导航，该模型如何搭建？需要从资源保障和绩效评估两方面构建各种指标体系。

ARC模块的作用是得到星历表T文件；MODEL模块的作用是将理论值与观测值的偏差存入C文件，并对其进行编辑和估算；AUTCLN模块由两个模块组成，分别为单差修复周跳模块（SINCLN）和双差修复周跳模块（DBCLN），在GPS数据处理中，正确探测与修复周跳是一个非常重要的问题[10]；CFMRG模块的作用是生成观测方程M文件；SOLVE模块的作用是利用最小二乘方法解算未知参数。H文件是GLOBK和GAMIT的数据交换中介。

图1 GPS数据处理流程图

GAMIT求解天顶延迟量的数据处理过程主要分为7个步骤[11]，具体流程见图2。

（1）建立测站先验坐标文件L-file。值得注意的是，当用于水汽反演研究时，要求初始坐标必须为精确的坐标。

（2）执行makexp命令。这一步主要是给定数据处理的时间，包括开始时刻、采样频率和计算历元，查看所需的文件是否齐全，建立所有准备文件的输出及一些模块的输入文件。这一步结束以后，会生成一个session.info文件，其中包含了数据处理相关的时间信息和卫星信息。

（3）执行sh_sp3fit命令。这一步是使用伪观测（从观测文件所得的位置和速度）获得精确的卫星初始状态，得到T文件和Y文件。

（4）通过执行makej命令，得到J文件。

（5）执行makex命令。通过执行这条命令可以生成两个文件，分别为X文件和K文件。X文件主要是对观测文件分析后按GAMIT需要的格式重新写过，K文件是利用台站的初始坐标、卫星广播星历和伪距信息计算出的接收机钟差。

（7）执行csh b*.bat。批处理工作由ARC（轨道积分）、MODEL（组成观测方程）、AUTCLN（自动修复周跳）和SOLVE（最小二乘法求解参数）等模块组成。

图2 GAMIT求解天顶延迟量的数据处理过程

3 实例分析

3.1 GPS站的选取

在局域GPS观测网（小网）反演大气可降水量时，如果想提高反演结果的精度，则应该适当的添加一些辅助站，组成一个监测网进行解算，这样就能减小监测网内水汽的相关性[12]。本文以闽南地区的厦门站（站号：59134，下同）为基准站，另外选取了该地区的其他5个站点联合进行解算，分别为：莆田（58946）、平和（59125）、平潭（58944）、南靖（59124）、福清（58942）。

3.2 GAMIT软件数据处理参数配置

本文采用GAMIT软件的版本号为10.34版。选取6个GPS永久跟踪站2017年5月BT的观测数据、以及精密星历和导航文件等数据。数据处理是在惯性坐标系J2000中进行的，数据处理观测量为无电离层延迟LC观测值的双差组合观测量[13]，处理时对厦门跟踪站的坐标进行强约束。资料处理的参数约束和方案见表2。

表2 资料处理中的参数约束和方案

相关的文件配置好以后，就可以进行后面的数据解算过程。主要的数据处理过程如下：

（1）利用厦门及其它GPS基准站组成的观测网来综合解算厦门地区的对流层天顶总延迟（Zenith Total Delay，ZTD）；

（2）利用相对精度较高Sastamonien模型求出厦门地区的天顶干延迟量（Zenith Hydrostatic Delay，ZHD），由ZTD减去ZHD，得到对流层湿延迟（Zenith Wet Delay，ZWD）[14]；

（3）利用厦门地区2016年的实际气象资料统计拟合出当地的加权平均温度模型，从而得到加权平均温度Tm和反演算法中的转换因子Π[15]；

Π称为无量纲比例因子，即水汽转换因子。k1、k2为常数，ρw表示湿大气的摩尔质量，Tm为大气加权平均温度（K），Rv为水汽常数。

（4）湿延迟量乘上由加权平均温度求出的转换因子，即PWV=Π·ZWD，便得出厦门地区的大气水汽含量。

3.3 GAMIT解算结果分析

判断结果是否正确主要有2个判断标准：（1）是否具有可进行合理估计的足够数据；（2）数据的噪声水平是否与选用的模型匹配[16]。第一条标准的首要指标是基线分量误差的大小，如果比预先设定的站点坐标和轨道参数的先验约束还要大，则会发现Q-文件或autcln.sum文件中大量数据被autcln丢弃。对于第二条标准，主要的指标是解决nrms（标准化的均方根误差），也就是chi-square x2的自由度。在Q文件中可以查看解算出的结果是否符合标准，本例的nrms为0.3，符合要求结果低于0.5的标准。6个测站的基线边长度及精度见表3。

表3 基线边解算精度统计

本文以解算厦门2017年4月BT水汽结果为例。图3为利用GAMIT软件解算出的对流层天顶总延迟（Zenith Total Delay，ZTD）时间序列图。图4为利用萨氏模型计算出的对流层干延迟分量（Zenith Hydrostatic Delay，ZTD）时间序列图。图5为利用天顶总延迟与干延迟分量相减间接求出的对流层湿延迟（Zenith Wet Delay，ZTD）时间序列图。

利用经验模型计算湿延迟需要很多的地面气象参数，给很多实际应用造成不便，而且会增加计算湿延迟的误差。因此，地基GPS气象学中通常是从天顶总延迟中减去干延迟分量从而得到天顶湿延迟，而不选用经验模型计算天顶湿延迟[17]。

3.4 厦门地区加权平均温度的局地算式

图3 厦门站ZTD的时间序列图

图4 厦门站ZHD的时间序列图

图5 厦门站ZWD的时间序列图

大气加权平均温度Tm的精度制约着水汽转换因子的精度，因此要想反演出高精度的水汽结果，就必须保证加权平均温度的精度[18]。Tm的获取途径有4 个[19]：（1）常数法，取 Tm=281 K，但是这种方法的精度较低；（2）利用当地探空数据采样积分公式计算；（3）利用数值预报得到加权平均温度的值；（4）利用当地气象数据算出适合当地当时段的回归经验公式。目前，由于气象探空数据的时空分辨率较低，而且高空的数据不易获得，常数法的精度又太低，数值预报的应用还不够成熟，因此一般都采用回归经验公式的方法。本文选取厦门2016年的气象数据资料，经过统计分析，求出当地的加权平均温度模型。图6为厦门地区地面气温与加权平均温度的散点图及趋势线。

图6 厦门地区地面气温与加权平均温度的散点图及拟合线

从两者地面气温和加权平均温度的趋势图可看出，所有的散点基本上分布在一条直线附近，说明加权平均温度与地面气温之间是线性相关的，因此可拟合出Tm和Ts的线性关系：

通过对加权平均温度和地面气温值进行统计分析并拟合，求出厦地区加权平均温度的线性回归关系为Tm=106.7+0.583 1Ts，均方根误差RMS=3.045 K。

3.5 水汽结果求解与分析

局地加权平均温度求出之后就可以确定转换系数Π，将天顶湿延迟ZWD乘上一个转换系数Π可以转换得到大气可降水量PWV[20]。图7为厦门站2017年4月BTGPS反演的PWV与由探空资料反演求得的PWV的序列比较图，图8为2种方法求得的PWV之差的时间序列图。

从图7可以看出，GPS反演的水汽含量较好的反应出水汽的积累和释放过程，水汽总量大致是以迅速增加、相对稳定、迅速降低的规律呈周期变化；并通过与探空资料反演的PWV比较可以看出，两种方法得到的大气可降水量具有良好的相关性。

图7 GPS与探空反演的PWV比较

图8 GPS与探空反演的PWV之差

从图8可以看出，GPS反演的PWV与由探空资料反演的PWV相比，两者偏差最小为0.3 mm，最大达到8 mm，其偏差的均方差为4.3 mm，均在允许的误差范围之内。从而说明，利用GAMIT软件可以很好的解算PWV。

值得一提的是，GPS反演的PWV与探空PWV的相关性受大气水汽含量的影响，当大气水汽含量较大时，GPS反演的PWV与探空PWV的平均偏差会增大，两者的相关系数也会变小，这可能与GPS水汽反演方法有关，而与GAMIT软件求解的天顶延迟无关[21]。

4 结论

基于闽南地区厦门及其他5个GPS站点的观测资料和从IGS下载的精密星历数据，通过GAMIT软件计算出对流层湿延迟，并利用厦门地区的实际气象资料，通过回归的方法得到厦门地区加权平均温度局地算式，进而求出了该地区的大气可降水量，并将最终反演结果与通过探空方法计算出的大气可降水量进行对比分析。

（1）利用GAMIT软件可以很好的解算PWV。

（2）大气加权平均温度进行本地化订正有利于提高GPS反演PWV的准确度。

大气水汽含量的多少对GPS反演PWV准确度的影响是一个值得重视的研究问题。