中国分区大气加权平均温度模型的建立与应用研究

高德恒，张 伟，沈清华

(中水珠江规划勘测设计有限公司，广东 广州 510610)

GNSS(全球卫星导航定位系统，Global Navigation Satellite System)气象学概念最早由Bevis[1]提出，是GNSS技术的一个重要应用领域，通过估算卫星信号在对流层中传播的总延迟量，剔除可精确计算的静力延迟部分，再乘以转换系数(Π)即可得到高精度的PWV(可降雨量,Precipitable Water Vapor),可用于气候研究及短临天气预报等。Rocken等[2]基于Bernese软件、精密轨道产品和实测地面气象参数，得到GNSS水汽值序列，其与同期水汽辐射计观测结果相比，RMS优于1 mm；Tregoning等[3]分析了基于GAMIT解算的GPS水汽、微波辐射计水汽及雷达探空水汽结果，两两间RMS优于1.5 mm；叶世榕等[4]应用精密单点定位方法估计对流层延迟，与IGS(国际GNSS服务，International GNSS Service)相应产品相比，精度优于6 mm，意味着基于该方法进行水汽反演精度能够达到1 mm，江鹏[5]基于GAMIT网解和数值气象资料插值的地面气象参数反演得到GPS水汽，其与探空水汽相比RMS小于3 mm。以上研究表明，采用精密星历和钟差产品，可以获取亚厘米乃至毫米级精度的对流层延迟，进而反演得到毫米级精度水汽结果。

在GNSS水汽反演时，除地表气象参数精度影响，求解Π的主要变量大气加权平均温度(Tm)也是影响水汽含量转换精度的关键参数。Tm是测站上空水汽压和绝对温度沿天顶方向的积分函数，利用无线电探空资料进行数值积分是公认的最为精确的计算方法，但是单次探空成本高，无法满足高时空分辨率的研究应用。因此Bevis[1]首次使用美国地区8 718次探空资料，建立了Tm与地面温度(Ts)的线性回归模型，并被广泛应用于全球其他区域的GNSS水汽反演研究；Wang等[6]基于ECMWF和NCEP中心发布的重分析数值气象产品及IGRA探空数据等研究全球Tm变化情况，发现Bevis公式在不同地区的适用性不一致；刘焱雄等[7]探讨了在香港地区使用Bevis经验公式引起的系统误差，进而建立了香港与地表温度回归模型；谷晓平等[8]利用广东清远探空站资料建立了适合该区域的经验公式；Raju等[9]使用8个探空站的观测资料，建立了印度经验模型；Sapucci等[10]基于ERA-40数值气象产品及探空资料建立了巴西Tm模型；谢劭峰等[11]选择广西地区4个探空站，开展了Tm建模与GPS水汽反演研究；Li等[12]基于3个探空站观测记录，建立了湖南省Tm模型。这些研究成果均表明Tm模型参数具有地域差异性，建立局部区域模型有助于提升水汽反演精度，但是在模型的适用范围确定方面，更多依据于行政区域的划分，具有较强的主观经验性。

在实际应用中，绝大多数GNSS站点没有并址的探空观测数据，只能利用临近的探空数据建立大气加权平均温度与地面温度的回归模型，再根据该模型和地面温度观测值计算大气加权平均温度。中国地域辽阔，地形复杂，气候类型多变，在开展大范围GNSS站点的水汽反演研究时，若建立统一的Tm模型，必定会出现在不同区域的适用性不一致问题；若根据各探空站资料分别建立局地Tm模型，不仅工作量增加，而且模型的适用范围界定标准亦不明确。因此，本文探讨了根据历史地面气温资料实现对中国的合理分区，再基于各分区内的探空资料建立分区Tm模型，最后应用这些模型开展各分区内GNSS站点的水汽反演工作。

1 数据集

1.1 中国地面气象年值日数据集

地表气温数据可以在中国气象数据网申请，获取中国地面累年值日值数据集(1981—2010年)，该数据集为中国基本、基准和一般地面气象观测站观测结果，包括气压、气温、降水、风要素等日气候标准值数据。地面基础资料专项工作对历史信息化文件重新进行了质量检测，并对存在问题和分歧的站点文件进行了修订，其站点地理分布见图1。

图1 站点分布

1.2 无线电探空数据

无线电探空是将各种气象要素测量元件和通讯设备组装在一起，称作无线电探空仪，将其固定在气象气球上，随着气球上升，测量沿途各高度上的气象要素，并将观测结果发送至地面接收站。由于采用一次性观测设备，成本较高，难以实现高时空分辨率的观测。多数站点每天观测2次，于UTC(协调世界时,Coordinated Universal Time)0时和12时开始；部分站点每天观测4次，分别于UTC0时、6时、12时和18时进行。美国俄怀明州立大学网站(http://weather.uwyo.edu/upperair/sounding.html)提供覆盖全球的无线电探空站数据下载，包含无线电探空仪沿途不同高度上的多个气象要素值。以2015年1月1日0时沈阳探空站观测结果为例，其各层的气压、气温、露点温度、相对湿度观测结果见图2。本文选取了国内及周边区域共计90个探空站点，站点分布见图1。根据观测结果网页规则编程并下载了这些站点2015年观测记录用于模型参数估计。

a)气压廓线

2 理论与方法

2.1 地表气温数据处理和插值

Tm和地面气象参数如气温、气压、露点温度都具有较高地相关性，王晓英等[13]、李国翠等[14]尝试建立多因子的Tm回归模型，并与单气温因子Tm模型进行了比较，认为使用多因子建模未能显著提升精度。因此，本文只考虑Tm与Ts的单因子回归模型，而Tm与Ts间为线性关系，Tm分区与Ts分区两者结果近似等价。同时，宋艳华等[15]、金巍等[16]、韩微等[17]在对中国地面气温空间分布和温度区划分的研究结果均表明，中国的气温分布存在明显的空间聚集现象。

综上，基于中国地面累年值日值数据集(1981—2010年)提供的各气象站气象要素累年逐日均值，剔除数据异常站点，提取各站点相应气温累年日均值，计算均值，得到共计2 128个气象站30年日均气温。应用ArcMAP地理统计分析工具集中含障碍的扩散插值方法，插值得到全中国区域的30年日均气温。根据刘焱雄等[5]的分析，假定湿延迟量取极大值500 mm，若要保证水汽转换精度小于1 mm，则精度必须优于3.4 K，选定分区间隔为3.4 K，进而对插值后的结果按照等间隔进行分区。

2.2 计算与建模

Tm定义为：

(1)

式中Z——沿天顶方向的高度，m；Pv——水汽压，hPa；T——绝对温度，K。

由于探空资料观测记录非连续值，故需要将上式离散化后，进行数值积分：

(2)

Pvi=PW·rh=6.112e(17.62t/(243.12+t))·rh

(3)

如式(4)，使用探空数据首层地面气温值以及求取的Tm值，采用最小二乘方法求解模型参数。

Tm=Ts·X+ε

(4)

由最小二乘原理可以求取模型参数：

X=(TsT·Ts)-1·TsT·Tm

(5)

则残差向量可表示为：

V=Tm-Ts·X

(6)

进而可以计算残差均值：

(7)

残差RMS：

VRMS=(VT·V)/(n-1)

(8)

将残差大于3倍残差RMS的观测记录剔除，重复上述过程，直至全部样本残差值小于3倍RMS，得到最终模型参数。根据残差均值和RMS大小，判定模型值与真值的偏移和离散程度，其值越小，则参数越优。

3 研究结果

3.1 分区结果

建立的分区结果见图4，分区间隔为3.4K，可以看出，地面气温均值具有很好地集群性，全国共分为8个级别，13个区域。在中国东部地区主要为平原和丘陵地形，地面平均气温分布受纬度影响较大，受地形影响相对较小，从低纬度地区向高纬度地区，地表平均气温逐渐降低，与纬度高低一致；相反，在中国西部地区，多以山地、高原和盆地地形为主，其地面平均气温分布主要受地形影响，因而，相较于东部地区分区结果，没有呈现出自下而上的规律。可认为该分区结果是合理的。

图4 基于30年平均气温分区结果

3.2 基于探空数据的Tm建模分析

3.2.1层顶高度与观测层数对Tm计算影响

统计2015年的探空数据发现实有观测记录64 131次，若按每天2次观测记录计算，应有探空记录65 700次，探空成功率大于97.6%。图5为探空记录的观测层数和层顶高度的统计结果，可以看出观测层数在集中于5～50层，层顶高度9～15 km的观测记录的比重最高，约占总观测记录的80%。同时发现有层顶高度仅数百米及负值、气球升至空中后才开始采集观测数据、气压与高程观测值不符等情况存在，为保证计算结果的精度，在开展建模工作前，已剔除这些存在质量问题的观测记录。

a)观测层数统计

一般而言，在各探空气象观测仪器能够正常工作的情况下，探空层数越多，层顶高度越高，获取的站点上空的垂直剖面信息越详细，Tm计算结果也就越精确。为了验证探空层数和层顶高度对Tm计算结果的影响，分别进行了实验分析。

a)观测层数影响。从图5可以看出，5层以上20层以下的观测记录有39 554次，占比约61.7%，为保证这些观测记录中的大部分能够应用于建模，重点分析了在该层数区间，每增加一层观测值，Tm计算的精度改善情况。过程如下：筛选出观测层数在20以上的全部观测记录，在保留首层和顶层数据的条件下，按层号对中间层均匀重采样至目标层数，并通过重采样后的观测记录计算Tm值，与重采样前观测记录计算出的Tm值比较，并计算RMS，分析观测层数对Tm计算的影响。目标层数从5依次增加至20，得到的RMS结果见图6，随着观测层数增加，RMS值逐渐变小，趋近于0，两者呈指数函数关系。对其进行指数函数拟合，两者关系可用图6中表达式近似描述，当观测层数为14时，Tm计算精度约为1 K。

图6 观测层数对Tm计算精度影响

b)层顶高度影响。与分析观测层数对Tm计算精度方法类似，探空观测记录的层顶高度集中在9～15 km，因此筛选出层顶高度大于15 km的观测记录，并按照不同层高假设对这些观测记录重新采样，然后计算相应的Tm值，与重采样前观测记录计算出的Tm值比较，分析层顶高度对Tm计算精度的影响。层顶高度值从5 km逐步升高至15 km，单次增加1 km，得到的RMS结果见图6，可以看出层顶高度越大，计算出的Tm与真值越接近，随着层高的增加，计算出的Tm与真值RMS逐渐趋于收敛。两者也存在指数函数关系，并可以使用图7中关系式近似表达，当层顶高度为9 km时，Tm计算精度约为1 K。

图7 层顶高程对Tm计算精度影响

根据上述分析，综合考虑观测记录的层顶高程和观测层数分布情况，最终取观测层数大于8且层顶高程大于9 km的55 729个观测记录开展进一步研究。

3.2.2建模结果

基于全部站点探空观测记录得到的Tm-Ts对，建立全国统一的回归模型，记作M1，结果见式(9)，拟合残差RMS为5.04 K，与Bevis模型拟合残差4.7 K相近。

Tm=54.86+0.773Ts

(9)

图8 全国Tm模型拟合结果

根据3.1节全国分区的结果，采用各分区探空站点数据，分别进行参数拟合，建立了12个(Ⅷ-2区域内没有探空资料)分区模型，结果见表2。从残差结果可以看出，在Ⅳ-1、Ⅶ-1和Ⅷ分区，全国模型与分区模型精度相当，但在其他分区，全国模型计算得到值与真值存在系统性的偏差，偏差均值最大达-6.98 K，各分区模型残差RMS均优于全国模型。因此，全国统一模型在中国不同区域的适用性不一致，建立分区大气加权平均温度模型优于使用统一的模型。

表1 分区Tm回归模型

3.3 基于分区Tm模型的GNSSPWV反演

根据GNSS站点和探空站概略坐标，得到两者距离在10 km以内的全部站点对共14对，详细信息见表2，站点分布见图9。

表2 站点对详细信息

图9 研究站点对分布

基于武汉大学精密单点定位软件Panda，计算得到GNSSPWV，探空观测记录中已给出本次PWV值，可直接使用。统计结果见图10(横轴为GNSS站点编码)，除百色(GXBS)和腾冲(YNTC)站外，其他站点的PWV偏差均值均较小，表明二者获取的PWV值较一致；若以RS PWV作为真值，则RMS反映了本文GNSS PWV反演的精度，只有6个站点精度在3 mm以下，而其余8个站点反演精度均大于3 mm，最大的为百色站，达到了5.5 mm。推测可能由以下几种原因造成。

a)偏差的值

a)站点使用的地表气压参数是使用数值气象产品插值得到，精度不足，将干延迟计算的误差引入到水汽反演中。

b)本文的ZTD结果是使用前12 h的GNSS观测数据解算得到，探空获取的则是观测开始后一段时间内的大气垂直剖面信息，在雨季、极端天气等情况下，空气中的水汽变化是十分活跃的，观测时段的差异很可能是造成RMS增加的一个重要因素，图10b可以看出RMS较大的站点多处于中国南方，是空气中水汽含量和变化均较大的区域。

c)探空站点与GNSS站点并非并址，因地形及高程等差异也很可能造成两者水汽不一致，譬如百色、蒙自和腾冲等站点处，两者偏差均值均大于1.9 mm。

4 讨论与结论

本文根据历史地面气温资料实现对中国的合理分区，建立分区Tm模型，并应用这些模型开展各分区内GNSS站点的水汽反演工作。通过回归建模得到了全国统一的以及分区域的大气加权平均温度回归模型，两者对比可以发现，分区Tm模型能有效地削弱全国模型在不同地区适用性不一致的影响，拟合残差RMS平均降低了17.8%；并且，基于该分区模型反演的可降雨量序列，与同期的相邻探空站点的观测数据计算得到的可降雨量序列变化趋势一致、大小相近，进一步证明应用分区Tm模型开展GNSS水汽反演是行之有效的。