基于IGU星历的单基站CORS大气可降水量估计

2019-02-28吴旭祥郭秋英桑文刚姜英明

山东建筑大学学报 2019年1期

吴旭祥郭秋英桑文刚姜英明

(1.山东建筑大学测绘地理信息学院，山东济南250101；2.千寻位置网络有限公司，上海200438)

0 引言

水汽是大气重要的组成成分，从区域的短期天气变化到全球的气候变化均发挥着至关重要的作用。大气可降水量PWV(Precipitable Water Vapor)反映了大气中水汽含量的丰富程度，是产生降水的基础，也是影响天气变化的关键要素之一，直接关系到短期气象降水及长期气候变化[1]。PWV是形成气候环境的重要特征量，其长期变化特征反映了气候的演变趋势。因此，实时监测PWV在气象学研究和业务领域中具有重要的意义。

目前国际上常规地基气象观测主要依靠探空气球、微波辐射计等手段获取大气水汽含量。由于探空气球发射间隔为12 h且探空站分布稀疏，所以通过探空气球得到的PWV的时空分辨率比较低，难以满足中小尺度数值天气预报的要求。微波辐射计的观测时间分辨率较高，通常5～8 min观测一次，但其在雨雪天气或者浓云等大气环境下会产生较明显的观测误差，不能满足全天候的观测要求[2]。为弥补常规大气水汽探测技术的局限性，将地基GPS观测资料用来反演大气水汽含量，并展现出广阔的应用前景。该技术与常规的大气水汽探测技术相比，具有低成本、全天候稳定实时观测、高时空分辨率等优点，有效弥补了目前常规探测手段探测大气可降水量的不足，能为改善中小尺度数值天气预报和气象学研究提供重要的、较为理想的大气水汽观测数据。

一直以来，欧美等发达国家将GPS气象成果广泛地应用于数值天气预报、气象研究等相关领域。Michaet等[3]推导出天顶湿延迟ZWD(Zenith Wet Delay)和PWV的转换关系，首次将GPS应用于气象研究。Kanda[4]研究发现实际降水和1 h内PWV增量关系密切，观测结果对于短期强降水的预报准确率提升明显。Wayan等[5]研究地基GPS在强对流天气下PWV的变化情况，验证了地基GPS在强对流天气下PWV检测的可靠性。

国内自20世纪90年代中期开始进行地基GPS气象学的研究，利用多种方法对大气可降水量开展了诸多的研究工作。王小亚等[6]首次完成GPS估计对流层延迟实验，得到较高精度的计算结果，初步验证了地面GPS观测为气象服务的可行性和可靠性。陈小雷等[7]分析了地基GPS反演PWV的精度情况，发现GPS反演的PWV具有较高的使用价值。赵静旸等[8]利用ECMWF再分析资料评估了中国地区GPS/PWV能达到的精度，与再分析资料计算的GPS/PWV与气象实测资料相比，年平均均方根偏差RMSD＜1 mm。周茂盛等[9]检验了连续运行(卫星定位服务)参考站 CORS(Continuously Operating Reference Stations)反演PWV的精度，并将其与美国环境预报中心再分析资料对比分析，结果表明通过区域CORS站网反演PWV达到了较高的精度。王洪等[10]比较了济南地区GPS/MET与微波辐射计和L波段探空设备反演PWV的数值偏差特征，发现三者变化特征一致。

我国各地雷暴、冰雹、暴雨等中小尺度的灾害多。中小尺度天气系统时空尺度小、生命周期短，而且水汽等气象要素的水平变化很大，因此对中小尺度天气系统进行监测、数值预报和研究中小尺度天气的发生、发展机理具有重要的理论意义和一定的社会经济效益[11]。随着GPS技术的发展，CORS系统逐渐应用于测量领域，在空间测绘数据信息采集和处理过程中起着重要作用。利用单基站CORS系统反演PWV，可实现对大气水汽的低成本、全天候、稳定和实时监测，即使在各种极端天气情况下，也可实现对大气水汽的高精度连续监测，可为气象和气候研究提供大量的科学数据，用以满足中小尺度天气分析、数值预报等需求。

现阶段，国际GPS服务中心IGS(International GPS Service)发布的最终精密星历产品IGF(IGS Final products)获得的PWV精度可以达到2 mm，能够满足气象研究和数值天气预报的要求[12]。但IGF发布时延达13 d，难以满足中小尺度天气系统近实时、短周期的研究和应用需要。2000年起，IGS中心开始向用户提供可实时获取的超快速星历产品IGU(IGS Ultra-rapid products)，其轨道时长为48 h，前24 h是实测轨道，后24 h是根据前24 h数据外推的预报轨道，IGS中心每天于世界标准时间3∶00、9∶00、15∶00、21∶00 发布超快速星历产品，每 6 h 更新一次[13]。使用当天21∶00发布的超快速星历产品的预报轨道在时间跨度上能够满足之后一天获取PWV的需要，在一定程度上能为小尺度天气系统提供较全面、快速的观测。

文章利用GAMIT软件，分别使用IGU和IGF星历产品参与解算单基站CORS的观测数据，验证解算精度后将估计的PWV与欧洲中期数值预报中心ECMWF(European Centre for Medium-Range Weather Forecasts)再分析资料产品ERA-Interim插值计算得到的PWV进行比较分析，验证IGU星历在单基站CORS估计PWV的可行性。

1 单基站CORS反演大气可降水量原理

在中纬度地区，高度大约＜12 km的大气层称为对流层，其集中了大气总质量的75%和几乎全部的水汽。GPS信号在大气中传播时，受到对流层中大气成分的影响，其传播速度减慢，传播路径发生弯曲，从而产生信号延迟。GPS气象学利用GPS信号在对流层延迟的某些特性反演PWV。

对流层中的GPS信号延迟，是由90%的干燥空气和10%的水汽造成的[14]。其中，由干燥大气引起的信号延迟称为干延迟或天顶静力学延迟ZHD(Zenith Hydrostatic Delay)，由水汽引起的信号延迟ZWD，ZHD和ZWD共同组成对流层总延迟ZTD(Zenith Troposheric Delay)。PWV可由ZWD计算得到，其关系由式(1)表示为

式中:Π为与对流层温度有关的转换系数。

ZWD受水汽影响变化较为复杂，难以用数学模型精确计算。天顶方向上，可以用ZTD减去ZHD得到天顶湿延迟ZWD，由式(2)表示为

ZHD可以通过 Saastamoinen模型[15]、Hopfield模型[16]推算得到。选择Saastamoinen模型计算静力学延迟，计算方法由式(3)表示为

式中:ZHD为天顶静力学延迟，mm；Ps为测站地表气压，hPa；λ为测站地理纬度；H为测站大地高，km。

ZTD可以使用GAMIT/GLOBK等软件计算出来，精度可以达到mm级，计算方法由式(4)表示为

式中:EL为卫星高度截止角；ATDEL为与卫星高度截止角EL有关的大气延迟量；DRYMAP为与EL有关的干延迟模型；WETMAP为与EL有关的湿延迟模型。

式中:k′2、k3分别为大气折射率试验常数，一般取固定值；Rv为水汽的气体常数，J/(mol·K)；Tm为加权平均温度，K。Tm由式(6)表示为

式中:Pv为水汽分压，hPa；T为同一大气断面上的绝对温度，K。

Tm是影响GPS反演PWV精度的一个重要参数，确定加权平均温度常用的方法有2种:(1)使用探空资料或数值分析资料等获得，该方法精度高但不能满足GPS反演PWV实时性的要求，因此通常用于事后分析；(2)通过地面气温Ts的经验公式计算得到，关系式由式(7)表示为

Bevis等[17]结合无线电探空资料给出了经验系数a为70.2、b为0.72。由于地理环境上的差异，Bevis模型应用于中国地区会产生较大的模型误差，为了获得最优化的PWV，不同地区可以通过适用于当地的回归模型进行修正。李建国等[18]对中国东部地区Tm进行回归分析，得出适合我国东部使用的经验公式，Tm与Ts相关系数＞0.9，由式(8)表示为

2 实验数据与方法

2.1 实验数据

2.1.1 单基站CORS数据

山东建筑大学CORS站为单基站CORS，站点名标记为 JZDX，使用的仪器为广州中海达提供的VNet 9北斗参考站接收机和配套扼流圈天线。观测天线放置于校内教学楼顶，四周无高大物体遮挡，对空通视良好。采集的GPS观测数据，单天观测时长为24 h，采样间隔为30 s，选用观测时间为2018年4月30日至5月30日，即年积日120～150 d，共31 d的观测数据作为实验的数据来源。

按照地基导航卫星遥感水汽观测规范标准，需对观测文件的数据质量进行检核，规范要求高度角＞10°的观测量中应有＞95%的有效观测量，平均多路径误差MP1、MP2应≤0.5 m，信号比率CSR(Clutter-to-Signal Ratios)值应≤0.5[19]。其中，数据利用率是实际采集数据量与应采集数据量的比值，反映了GPS接收机采集数据的能力，利用率越高表明GPS接收机接收数据的能力越强；MP1、MP2反映扼流圈天线对多路径效应的抑制效果，其值越小表明抗多路径效应的能力越强；CSR能够反映出接收机锁星的能力，其值越小说明GPS接收机的周跳修复能力越强。JZDX测站点的单天数据质量均符合规范要求，JZDX站点观测数据120～150 d统计结果的平均值见表1。

表1 JZDX测站观测质量统计表

2.1.2 ECMWF再分析资料

以“真人”或“神人”为理想人格的道家人格。道家主张清静无为、无为不争、少思寡欲、淡泊名利，来恢复人性的原始和质朴状态，把自己看作客观世界的一部分。正如老子认为的“圣人处无为之事，行不言之教”，体现了道家自然朴素、身心合一的理想人格境界。道家在提倡自然心性的同时，也强调了人们的精神自由和意志独立。道家理想人格的修养方法有“坐忘”“心斋”等，在道家看来，泯灭物我和人我是人格修养的重要阶段，最终消除自我与外物、外界的矛盾。道家的返朴归真和自然无为，对当代一部分中国人，特别是一些文人墨客的人生态度和生活情趣产生了深刻的影响，使很多人以开朗豁达的心态去面对人际交往和处理社会矛盾，不断提高自身人格修养。

在利用GPS反演PWV时，GPS测站相应的气象数据是反演PWV必不可少的参数。山东建筑大学单基站CORS用于大地测量学的应用和研究，并未配备相应的气象传感器，且该站15 km内无可用气象站和探空站，很难获取GPS反演PWV所需要的气象观测数据。针对这一问题，采用ECMWF再分析资料作为GPS反演PWV的气象资料。

ECMWF是国际上权威的天气预报研究和业务机构，提供高精度的再分析模式和预报产品。自1979年开始运行制作中期天气预报，提供的再分析产品具有空间分布广、资料完善、精度高等优势。

采用ECMWF发布的ERA-Interim再分析资料，从地面到最上层0.1 hPa(高度约64 km)，空间分辨率为0.125°×0.125°，分析产品为逐日4次(0、6、12、18 UTC)，每12 h 进行一次误差校正，保证产品精度与可靠性。ERA-Interim再分析资料中的地面层分析场提供18种观测要素，其中要素总水汽含量对应GPS气象学中的PWV，要素2 m气温对应近地面气温Ts。ERA-Interim再分析资料为格网数据，为了获得精度较高的数据插值结果，通过与JZDX站点相邻的4个格网点到测站的距离，采用双线性插值的方法得到JZDX站点的总水汽含量值及地面气温值。

2.1.3 IGS站数据

PWV通常被视为一个绝对量，因此必须求得绝对对流层延迟。为获得绝对的ZTD值，使用GAMIT/GLOBK软件解算时基线长度应＞500 km，确保得到的ZTD是独立的估计值。因此在数据处理中引入IGS站进行联合解算，引用JZDX测站周边的6个 IGS站参与数据解算，分别是房山站(BJFS)、长春站(CHAN)、黄石站(HKWS)、上海站(SHAO)、桃园站(TWTF)和乌兰巴托站(ULAB)。

2.2 实验方法

GAMIT/GLOBK的解算方法比较复杂，由多个子模块和相应的控制程序配合以实现数据的解算。GAMIT/GLOBK软件中测段信息控制文件Sestbl定义了常用的解算控制设置，是数据处理的核心控制文件，具体参数设置见表2。

表2 GAMIT/GLOBK解算主要参数设置表

在将解算数据进行预处理之后，(1)执行makex命令，生成后续命令的控制文件。其次执行sh＿sp3fit命令，分别根据IGF和IGU星历生成轨道初值文件。再次执行makej命令，由广播星历生成卫星钟文件。(2)执行makex命令，将JZDX站及其周边的6个IGS站的观测数据转换成双差观测文件。另外执行fixdrv命令，根据解算控制文件生成一系列的批处理脚本。(3)执行csh b*.dat启动解算脚本及生成成果文件。(4)执行sh＿metutil命令，处理包含对流层参数解算结果的文件并生成含有逐小时PWV的文件。

3 基于IGU星历大气可降水量解算

3.1 解算质量评价

3.1.1 基线解算质量评价

(1)单天解标准化均方差NRMS

基线解算的精度与可靠性是判断解算可降水量精度的标志之一[20]。其中，单天解标准化均方差NRMS表示基线值与加权平均值的偏离程度，是评价基线解算质量的重要指标。NRMS值越小，表示基线解算精度越高，通常情况下NRMS值＜0.3 m可视为质量合格，而＞0.3 m时则表示解算过程中存在参数设置问题或周跳未完全修复[21]。

使用IGS最终精密星历IGF和超快速星历IGU分别进行基线解算后得到62个单天基线正规解，其NRMS统计结果见表3。使用IGU和IGF解算的基线解，其偏倚值Bias分别为3.5和3.7 mm，均方根误差RMSE分别约为4.9和5.2 mm，表明2种星历解算的结果合格且稳定性较好，可以用于实验分析。

表3 NRMS统计结果表

(2)基线重复率

基线重复率是衡量GAMIT/GLOBK基线解算质量的重要指标之一，它反映了基线内符合精度，其值越小，基线内符合精度越高。NRMS满足要求后，应提取结果文件计算基线重复率的情况。在基线解算中，与JZDX测站相关基线共3条，单天解基线重复性的均值统计情况见表 4。可以看出，GAMIT/GLOBK解算的基线质量较好，绝对基线重复性均＜3.5 mm，相对基线重复性在10-9量级，满足高精度GPS基线处理要求。

表4 单天解基线重复性统计表

3.1.2 GAMIT解算的ZTD精度评价

由于地基GPS反演PWV的精度很大程度上取决于ZTD，因此将GAMIT/GLOBK使用IGU解算的ZTD结果与IGS发布的ZTD产品进行对比分析，由此验证GAMIT/GLOBK解算ZTD的精度。6个IGS站的GAMIT解算值与IGS发布的ZTD产品的数值偏差情况如图1所示。可以看出6个IGS站的ZTD偏差基本上＜10 mm，异常值较少，整体较为稳定。

图1 IGS站ZTD的解算值和IGS/ZTD产品的数值差异情况图(ULAB站年积日第136、137、139 d数据缺失)

6个IGS测站ZTD解算值和IGS/ZTD产品数值差异统计情况见表5。可以看出，Bias＜3 mm、RMS＜6.28 mm、RMSE＜0.22 mm，GAMIT/GLOBK 解算的ZTD与IGS/ZTD相关系数＞0.9。考虑到个别数值存在粗差，将＞2倍中误差(2σ)数值剔除掉之后，RMS明显减小，整体平均值＜5.11 mm。综上所述GAMIT/GLOBK使用超快速星历产品对ZTD影响在精度范围内且两者吻合程度较高，一致性较好。

表5 GAMIT/GLOBK解算得到的ZTD与IGS/ZTD产品数值分析统计表

3.2 解算结果对比分析

3.2.1 基于IGU的PWV与基于ECMWF的PWV对比分析

根据JZDX测站点的经纬度坐标，获取到ERAInterim再分析资料中与JZDX站点相邻的4个格网点到测站点的距离，采用双线性插值方法得到JZDX站点的PWV序列(ECMWF/PWV)。由于 ERAInterim再分析资料中PWV数据为每6 h一个值，所以将单基站CORS观测得到的IGU/PWV序列再按6 h抽取，以便与ECMWF/PWV作对比分析。

JZDX测站的ECMWF/PWV与IGU/PWV的Bias为-0.1 mm，RMS为3.31 mm，在剔除掉＞2倍中误差的数值(占比 4.07%)之后，其RMS减小为2.89 mm，接近于数值天气预报对PWV精度在2～2.5 mm范围内的要求。引起该误差的原因可能是样本容量有限，少量误差会对整体精度产生较大影响，还可能是进行空间双线性插值及时间抽取过程中产生的精度损失。

JZDX单基站CORS观测得到的ECMWF/PWV与IGU/PWV的时间序列对比情况如图2所示。单基站CORS观测得到的IGU/PWV和ECMWF/PWV序列在趋势上有较好的一致性，两者之间没有表现出明显的系统性偏差，说明单基站CORS使用超快速星历获得的PWV与ECMWF再分析资料成果具有很好的一致性。

图2 JZDX测站IGU/PWV与ECMWF/PWV时间序列图

JZDX测站ECMWF/PWV与IGU/PWV序列的线性拟合情况如图3所示。可以看出，IGU/PWV序列和ECMWF/PWV序列间的相关系数为0.96，两者表现为正相关且强相关；P值为0，表示差异具有高度统计意义。因此，通过 IGU/PWV与 ECMWF/PWV进行外部检核，可以认为基于超快速星历的单基站CORS估计的大气可降水汽含量稳定可靠，无明显的系统性偏差。

图3 JZDX测站ECMWF/PWV与IGU/PWV序列的线性拟合图

3.2.2 基于IGU的PWV与基于IGF的PWV对比分析

为了验证IGU超快速星历产品解算PWV的可靠性和精度，将其与精度最高的IGF最终精密星历产品进行对比分析。IGU与IGF 2种星历产品分别解算得到的PWV及其差异值如图4所示。可以看出，2种星历产品的Bias为0.005 mm、平均偏差为0.137 mm、RMS为0.178 mm，两者在数值上非常的接近，总体表现为很高的可靠性。

JZDX测站点IGF/PWV和IGU/PWV的线性拟合如图5所示。IGF/PWV和IGU/PWV序列在数值上相当接近，其相关系数达到0.99，表现为强相关性；P值为0，表示差异具有高度统计意义。综合以上可以认为单基站CORS使用IGU星历产品与事后发布的IGF星历产品解算PWV的精度相当，单基站CORS使用前一天发布的IGU星历产品估计当天PWV是切实可行的，精度是稳定可靠的。