基于GNSS技术探测大气的方法可以分为地基GNSS网层析技术、单站地基GNSS反演技术和空基GNSS掩星技术。地基单站GNSS接收机在定位中能够获取电磁波传播路径的对流层延迟信息，因此可以用来反演对流层大气折射率剖面。与地基GNSS网层析和空基GNSS掩星技术相比，地基单站GNSS技术具有观测成本低、操作简单、可移动性较强等优点^[1]。

Lowry等^[2]首次提出基于地基单GNSS测站获得的斜路径延迟来反演对流层大气折射率剖面。随后，林乐科等^[3-6]对单站地基GNSS探测大气折射率剖面进行深入研究，实现了基于对流层天顶延迟的折射率剖面反演，并在广东东南沿海地区开展单站地基GNSS探测大气的大型观测实验。目前，单站地基GNSS反演大气折射率的模型仍需进一步改进，以提高大气折射率的反演精度。本文提出一种新的利用天顶对流层延迟反演大气折射率的模型，并选用中国及周边区域8个IGS测站的观测数据进行试算。将测站附近探空产品获取的大气折射率作为真值，统计结果显示，与传统的反演模型相比，新的反演模型提高了大气折射率的反演精度。

1 传统的基于天顶对流层延迟的模型参数反演

利用平差模型，通过GNSS接收机接收到的电磁波信号可以估算出测站的坐标、天顶对流层延迟(ZTD)、接收机钟差等信息。天顶对流层延迟在导航和定位中是一种主要的误差源，而在GNSS气象学中，它却是一种探测对流层大气的信号源^[7-8]，根据是否含有水汽，可划分为天顶对流层干延迟(ZHD)和天顶对流层湿延迟(ZWD)^[9]。实际的反演中，根据不同的模型利用ZHD和ZWD分别反演出大气干折射率N_T和大气湿折射率N_W。

1.1 干项模型参数估计

干项大气折射率比较稳定，近似认为满足大气静力平衡，Hopfield^[10]推导获得满足大气静力平衡时的干项大气折射率为：

$ N_T^h \approx N_T^{{h_0}}{\left( {\frac{{{T_0} + \beta h}}{{{T_0}}}} \right)^{ - \left( {1 + \frac{g}{{\beta {R_d}}}} \right)}} $

(1)

式中, N_T^h₀为地面折射率静力项，T₀为地面温度(K)，h为相对于地面站的高度，β为高度方向的温度梯度，g为重力加速度(一般变化不大)，R_d为干空气的比气体常数。

ZHD可表示为干项折射率的积分^[11]：

$ {\rm{ZHD}} = {10^{ - 6}}\int_{{h_0}}^{{H_h}} {N_T^h} (h){\rm{d}}h $

(2)

由式(2)便可估算出折射率干项在对流层H_Tro内的积分^[11]：

$ \begin{array}{*{20}{c}} {{\rm{ZH}}{{\rm{D}}^{{\rm{Tro}}}} = {{10}^{ - 6}}\int_0^{{H_{{\rm{Tro}}}}} {N_T^h} (h){\rm{d}}h = }\\ {{{10}^{ - 6}}N_T^{{h_0}}\frac{{{{\left( {{T_0} + \beta {H_{{\rm{Tro}}}}} \right)}^{\eta + 1}} - T_\delta ^{\eta + 1}}}{{(\eta + 1)T_0^\eta \beta }}} \end{array} $

(3)

式中，$\eta=-\left(1+\frac{g}{\beta R_{d}}\right)$。可以看出，该积分与β有复杂的函数关系。

由式(3)可以估算出β，将其代入式(1)便可得到不同高度层的大气折射率干项。

1.2 湿项模型参数估计

从ZTD中剔除ZHD便可得到ZWD^[10]：

$ Z{\rm{WD}} = Z{\rm{TD}} - Z{\rm{HD}} $

(4)

ZWD可表示为湿项折射率N_w的积分：

$ Z{\rm{WD}} = {10^{ - 6}}\int_0^\infty {{N_w}} (h){\rm{d}}h $

(5)

N_w随高程呈现近似指数形式的分布^[11]：

$ {N_w}\left( h \right) = N_w^0\exp \left( {\frac{{ - h}}{{{H_{we}}}}} \right) $

(6)

式中，H_we为湿项的特征高度，N_w⁰为起始高度处的湿项折射率。

根据式(5)和式(6)可知^[11]：

$ Z{\rm{WD}} = {10^{ - 6}}{H_{we}}N_w^0 $

(7)

通过式(7)可估算出H_we，将其代入式(6)可计算出不同高程处的湿折射率。

2 新的基于天顶对流层延迟的模型参数反演 2.1 新的干项模型参数估计

干项折射率随高程的变化比较平稳，根据式(1)和式(3)，在获取ZHD和H_Tro后，利用搜索的方法便可计算出β，将其代入式(1)便可获得不同高度层的N_T^h。通常情况下，GNSS测站会配置气象观测文件，能够提供压强、温度和相对湿度。基于压强信息，利用Saastamoinen模型可以精确地估算出ZHD。H_Tro随区域和季节的变化不断改变，而传统上H_Tro通常取一个经验值60 km。本文试图基于Era-Interim产品估算H_Tro。

ERA-Intrim提供了不同高程上的压强、温度、湿度等气象产品，利用这些产品可以计算出不同高程上的大气干项折射率及ZHD：

$ N_T^h = 77.6 \cdot \frac{{{P_h}}}{{{T_h}}} $

(8)

$ {\rm{ZH}}{{\rm{D}}_h} = \frac{{{{10}^{ - 6}} \cdot \sum\limits_{i = 0}^{i = n - 1} {\left( {N_T^{{h_i}} + N_T^{{h_{i + 1}}}} \right)} }}{2} $

(9)

将起始高度的ZHD_h(即ZHD₁)视为总的天顶干延迟，而不同高程h上估算的ZHD_h占ZHD₁的比例为RZHD：

$ {\rm{RZH}}{{\rm{D}}_h} = \frac{{{\rm{ZH}}{{\rm{D}}_h}}}{{{\rm{ZH}}{{\rm{D}}_1}}} $

(10)

若高程h上RZHD小于等于0.01，则定义h所对应的高度为H_Tro。

2.2 新的湿项模型参数估计

研究无线探空、无线掩星等高垂直分辨率、高精度的气象观测值发现，N_h^w与P_w^h之间的关系可近似表示为：

$ \left( {\frac{{N_h^w}}{{N_{{h_0}}^w}}} \right) = {\left( {\frac{{P_w^h}}{{P_w^{{h_0}}}}} \right)^\gamma } = {\left( {\frac{{{T_0} + \beta \left( {h - {h_0}} \right)}}{{{T_0}}}} \right)^{\frac{{\gamma {g_0}}}{{\beta {R_d}}}}} $

(11)

式中，γ为大气混合比指数。

通过式(5)可知：

$ \begin{array}{*{20}{c}} {Z{\rm{WD}} = {{10}^{ - 6}} \cdot N_{{h_0}}^w\int_0^\infty {{{\left( {\frac{{{T_0} + \beta \left( {h - {h_0}} \right)}}{{{T_0}}}} \right)}^\gamma }} {\rm{d}}h = }\\ {{{10}^{ - 6}} \cdot N_{{h_0}}^w \cdot \frac{{{T_0}\left( {E + 1} \right)}}{\beta } \cdot }\\ {\left[ {1 - {{\left[ {\frac{{{T_0} + \beta \left( {{H_{{\rm{top}}}} - {h_0}} \right)}}{{{T_0}}}} \right]}^{\left( {E + 1} \right)}}} \right]} \end{array} $

(12)

式中，$E=\frac{\gamma g_{0}}{\beta R_{d}}$，H_top为对流层顶。在已知γ和H_top的情况下，利用式(12)便可估算出β，将其代入式(11)，便可得到不同高程处的N_h^w。通常情况下，γ取值为4^[12-13]，而H_top可以基于Era-Interim产品，利用CPT (cold point temperature)法或WMO(world meteorological organization)法来确定。

3 结果分析 3.1 数据来源

选择中国及周边区域8个IGS站(图 1)2016-07的观测数据及气象观测值，使用Bernese5.2软件及IGS事后轨道产品和钟差产品，卫星截止高度角设置为7°，然后选用精密单点定位的计算模型估算ZTD。基于Saastamoinen模型可以获得ZHD^[12]，然后从ZTD中分离出ZWD。

3.2 确定H_Tro和H_top

Era-Interim提供高垂直分辨率的大气压强、大气温度和湿度等大气产品，可从网站http://apps.ecmwf.int/datasets/data/interim-full-daily/下载。选择中国地区2006~2016-07的Era-Interim产品，分别利用式(10)和MWO方法估算出H_Tro和H_top，最后分别利用传统和新的模型参数反演出大气折射率。表 1给出中国及周边区域8个IGS测站2016-07的H_Tro和H_top。

从表 1看出，不同地区的H_Tro最大差距为16 km，而H_top的最大差距为5.09 km；LHAZ和URUM站等偏干旱区域的H_Tro和H_top要小于偏湿润的TNML站。

3.3 反演结果与探空数据的比较

从图 1中可以看出，每一个IGS测站附近都有一个探空测站，每天提供2次探测结果。探空数据可以提供高精度的、随高程变化的大气参数，如温度、压强、相对湿度和混合比等。利用探空产品获得的对流层折射率作为真值来检验本文方法获得的结果。图 2和图 3分别给出URUM和WUHN站单个反演剖面的结果。

从图 2和图 3可以看出，与传统的模型参数法反演的大气折射率相比，新的模型参数法反演的大气折射率更接近真值。当大气折射率在空间分布中存在“逆增层”现象时，新的模型参数法及传统的模型参数法反演的结果与真值存在较大的偏差。此外，给出2个测站模型参数法获得的结果与真值之间1个月的偏差分布图(图 4、5)。

图 4、5显示，大部分情况下新的模型参数法反演得到的大气折射率更接近真值。此外，统计8个测站模型参数法获得的结果与真值之间的偏差(表 2)。

从表 2可以看出，新的模型参数法获得的结果更接近真值，从而证实了其有效性。

4 结语

GNSS气象学使得对流层延迟成为反演大气参数的有效信号源。本文利用天顶对流层延迟值，基于模型参数法反演不同高程处的大气折射率，并利用探测产品获得的大气折射率作为真值来评估新的模型参数法的可行性。选择中国及周边地区8个IGS测站2016-07的观测数据进行试算，结果显示，与传统的模型参数法反演的结果相比，新的模型参数法反演的大气折射率更接近真值。然而，当大气折射率的空间分布中存在“逆增层”现象时，2种模型参数法获得的结果都与真值存在较大的偏差。如何反演大气折射率“逆增层”现象将是下一步重点研究的内容。