何秀凤，詹 伟，施宏凯

河海大学地球科学与工程学院，江苏 南京 211100

精确的大气水汽三维时空分布信息对于提升中尺度数值气象预报精度具有非常重要的意义[1]。当前对大气水汽进行探测的常规技术(无线电探空、气象卫星等)存在低时空分辨率、高成本、观测精度受天气影响大、观测范围局限等问题[2-3]。全球导航卫星系统(GNSS)最初被应用于区域大气可降水量[4-6]的估计，此后为了进一步得到区域内水汽的三维分布，GNSS三维水汽层析技术[7]兴起，凭借其低成本、高时空分辨率、全天候的特点被广泛应用。但传统算法[8-9]只考虑从层析区域顶端射出的信号，忽略了侧面穿刺的信号，层析方法存在观测数据利用率低、观测分布不均等问题。有效利用高度角较低的边界信号提供的信息可以增加层析区域内网格尤其是底层及边缘网格被穿过的概率、减小层析网格空格率[10]、提高观测数据利用率。

为解决边界信号未能有效利用的问题，文献[11—13]建立了一种基于水汽密度比例因子的三维水汽层析算法；文献[14]通过缩小层析水平范围构建层顶入射信号上的水汽含量与范围缩小后的层析边界高程间的比例关系，用于恢复边界信号的有效水汽信息；文献[15—16]附加一个辅助层析区域，使从研究区域侧面入射的信号也可被有效利用。以上3种模型在一定程度上克服了观测数据利用率较低的局限性，取得了较高的精度，但是均在单一导航系统中得以实现，在多模组合系统水汽层析中尚未有所研究，观测数据利用率有待进一步提高。在此背景下，文献[10]对以上研究进行了详细分析，通过层析时间前3 d的探空信息拟合出不同高程对应的大气可降水量(precipitable water vapor,PWV)与高程的函数关系，从PWV的角度恢复边界信号水汽信息，在多模组合系统中进行水汽层析试验，并证实了引入多系统信号及引入边界信号均能获取更好的水汽分布结果。

文献[10]基于拟合函数建立的边界信号的精度相对于层顶信号存在较大差别，有必要削弱这部分误差对层析结果的影响。与此同时，以上研究中的层析模型均是根据垂直方向上水汽分布呈指数递减规律[17]建立上下层的关系，得到的传统垂直约束方程与层析区域实际水汽分布符合程度较低，采用传统垂直约束解算的结果精度较差，在底层区域表现得较为明显。为了进一步提高三维水汽层析的精度，有必要对上述两个问题加以考虑。因此，本文基于短期探空信息拟合函数建立了层析区域内水汽分布的区域性垂直约束，采用了一种顾及边界信号的GNSS水汽层析方法，设计和实现了融合边界信号的多模水汽层析试验，采用无线电探空数据进行精度验证，详细对比分析了边界信号及本文建立的垂直约束对水汽层析结果的改善程度。

1 基于探空信息拟合函数的三维水汽层析算法

1.1 基本观测方程

倾斜路径水汽含量(slant water vapor,SWV)指的是GNSS信号传播方向上的水汽柱内所含水汽转化为单位面积上的液态水高度，单位为mm，其表达式为

(1)

式中，S为SWV；L为卫星到接收机的信号路径；ρwater表示水汽密度，单位为g/m3；dL表示单位截面积上的体长度，单位为km。将式(1)方程离散化为线性形式的观测方程

S=∑ijk(dijk·ρijk)=A·X

(2)

式中，dijk为信号ρ在网格(i、j、k)内的截距，单位为km；ρijk为(i、j、k)网格内的水汽密度；A表示斜路径信号在层析网格内穿过的截距向量，未穿过的网格截距为零；X为层析网格内部各个网格的水汽密度值。根据式(2)将所有GNSS观测表达式表示为矩阵形式

Am×n·Xn×1=Sm×1

(3)

式中，Sm×1=[S1S2…Sm]T表示所有GNSS斜路径信号上观测值向量；Am×n=[A1A2…Am]T为对应信号在网格内的截距矩阵；Xn×1为各个网格的水汽密度未知量；m表示卫星信号的条数；n为网格个数。

1.2 边界信号观测方程

针对传统模型中未能有效利用许多侧边穿刺信号的问题，文献[10]提出直接从PWV与高程的函数关系恢复边界信号水汽值的方法。其思路是：PWV作为求解SWV的关键元素，与高程存在函数关系，而与梯度改正项无关，因而利用五阶多项式拟合PWV与高程的相互关系，直接恢复边界信号的SWV。本文考虑到水汽在垂直方向上呈指数递减的规律，利用指数函数拟合两者的关系以恢复边界信号信息，并为层析高度的选取提供参考。其中，PWV与高程相互关系拟合过程如下。

(1)利用层析区域内探空站前3 d的探空数据，在垂直方向上计算出采集点累计PWV占整个垂直方向上PWV的比例因子γ。

(2)由步骤(1)计算的γ及其对应的高程H进行指数函数拟合，建立函数关系式γ=a·eb·H+c·ed·H。

(3)通过步骤(2)得到的拟合函数，只需要代入每条边界信号穿刺点处的高程H，便可由测站处PWV计算得到对应的位于层析区域内的有效PWVH=PWV·γ。然后采用SWV计算公式[18]计算出边界信号水汽值SWVH

SWVH=mw(θ)·PWVH+Π·{mΔ(θ)·

(GN·cosφ+GE·sinφ)+ε}

(4)

式中，mw(θ)、mΔ(θ)分别为湿映射函数和梯度映射函数，与卫星高度角θ有关；Π为水汽转换系数；GN、GE分别为南北方向和东西方向梯度参数；φ为方位角。综上建立起边界入射信号的观测方程，矩阵形式如下

AHm×n·Xn×1=SHm×1

(5)

式中，SHm×1为边界信号位于层析区域内的有效水汽值；AHm×n为边界信号的系数矩阵；m为信号条数。

1.3 水平约束与垂直约束

GNSS观测方程系数矩阵是一个绝大多数为零的稀疏矩阵，通过引入额外的水平和垂直约束条件[19]对观测方程进行改善。

水平约束条件根据同层网格之间的相关性，基于高斯加权法[11]对每一个网格建立约束条件，构成如下水平约束方程

H·X=0

(6)

式中，H为水平约束系数矩阵。

文献[20]通过大量试验指出，对流层低层的垂直分辨率相对较高，采用垂直不均匀分层方法得到的大气水汽更加精确，垂直不均匀分层的最佳高度为500～1200 m。传统方法中垂直约束根据垂直方向上水汽分布呈指数递减规律建立上下层的关系

xk+1=xk·e(-Z/H)

(7)

式中，xk+1、xk分别为上下两层的水汽密度参数；Z为上下层的高程差；H为水汽标高，通常取值1～2 km。

本文垂直约束由1.2小节中得到的拟合函数构建上下两层之间的函数关系

xk+1=α·xk

(8)

V·X=0

(9)

式中，V为垂直约束方程系数矩阵。

1.4 层析方程组解算

综合上述内容，引入GNSS多模观测数据得到顾及边界信号的层析方程组

(10)

式中，m=g+r+e+c；Hm=Hg+Hr+He+Hc；m、g、r、e、c分别代表GNSS组合和G、R、E、C各单系统的层顶入射信号数目；Hm、Hg、Hr、He、Hc分别代表GNSS组合和G、R、E、C各单系统的边界入射信号数目，即

添加约束条件后的层析方程组仍然是一个绝大部分元素为零的稀疏矩阵，如何对方程组进行解算将影响层析结果的准确性及稳定性，一般采用代数重构法[21-25](algebraic reconstruction techniques,ARTs)求解水汽参数，可以避免系数矩阵求逆的运算。本文采用加型ART进行方程组解算，其公式如下

(11)

式中，k表示第k次迭代；xk为第k次迭代的结果；Ai表示观测矩阵的第i行；Si为对应的SWVi；〈·,·〉表示向量内积；λ为松弛因子。

2 GNSS三维水汽层析

2.1 试验区域及数据

选择香港地区作为试验区域，获取了香港CORS的10个站点2016年DOY 129—DOY 135一周的GNSS观测数据和气象数据，采用CODE分析中心提供的精密星历、钟差、轨道等产品，基于数据处理软件RTKLIB进行GPS/GLONASS/Galileo/BDS四系统组合精密单点定位(precise point positioning,PPP)解算，估算对流层延迟参数，数据处理策略见表1。

表1 PPP数据处理策略Tab.1 PPP calculation strategy

研究区域网格划分：纬度范围为22.16°N—22.52°N，经度范围为113.90°E—114.35°E，经纬度网格密度为8×8均匀划分；垂直方向上由探空信息拟合函数计算得到8000 m高度(DOY 129—DOY 135)PWV比例因子均值为0.982，最小值为0.973，8000 m高度以上分布的水汽均值低于0.9 mm，最大低于1.6 mm，因而将层析高度定为8000 m。采用不均匀分层(0、500、1000、1700、2400、3100、4000、4900、5800、6900、8000，单位为m)，总计640个网格。利用KingsPark探空站无线电探空数据作为水汽层析结果的验证数据。研究区域及测站分布如图1所示。

图1 层析区域及GNSS测站分布Fig.1 Tomography area and distribution of GNSS stations

2.2 探空信息拟合函数精度分析

图2为2016年DOY 129—DOY 134根据探空站数据计算的PWV比例因子和高程相互关系的拟合结果。由图2可知：DOY 129—DOY 131 3 d采样点数据分布相对密集，对应的拟合结果均方根误差(RMSE)分别为0.019 2、0.017 3、0.021 1；DOY 132—DOY 134 3 d则相对分散一些，RMSE分别为0.039 1、0.049 9、0.052 8。图2结果表明，基于探空信息得到的拟合函数存在一定的模型误差，表现在PWVH上1～3 mm，由此模型计算的边界信号的SWV观测值相对于层顶信号的SWV精度要低一些，因此，有必要在层析方程组解算时采取合理的方法对这部分误差进行削弱。

2.3 层析区域内信号数量与空格率分析

为研究边界信号对层析模型的改善程度，本文首先统计了仅考虑层析区域天顶信号(图3(a)，以GNSS表示)和顾及边界信号(图3(b)，以GNSS+表示)的两种层析模型在30 min内的有效GNSS信号数量和空格率，其中数据采样间隔为1 min。

图2 DOY 129—DOY 134指数函数拟合曲线Fig.2 Curve fitted by exponential function of DOY 129—DOY 134

图3 两种层析模型Fig.3 Diagrammatic drawing of two tomographic models

融合边界信号后层析区域内的观测量增多，同时低高度角的边界信号使得层析区域底层及边缘网格被穿过的机率大大提高。为此，本文对DOY 129—DOY 135每天UTC 12 h层析模型的观测量、网格空格率及边缘网格(层析区域侧面最外围网格，不包括顶层)空格率进行了对比分析。其中层析区域内网格总量为640、边缘网格为198，统计结果见图4和表2。

图4和表2的结果表明：顾及边界信号的GNSS+模型在观测量、网格空格率，以及边缘网格空格率3个方面上相对于舍弃边界信号的GNSS模型都存在优势。在层析解算时刻：GNSS方案能够利用的平均有效观测量只有4637，而GNSS+为7045，观测量提高了51.9%，并且GNSS+模型的最低观测量都要高于GNSS模型的最高观测量；GNSS+方案的平均空格率为14.3%，相对于GNSS的26.1%降低了11.8%；边缘空格率则表现得更为明显，相对于GNSS方案47.1%的平均空格率，GNSS+降低了24.1%，只有23.0%。

2.4 三维层析结果精度分析

本文利用2016年DOY 129—DOY 135一周的数据进行三维水汽层析精度分析，Kingspark探空站位于层析区域内，并且其配备的无线电探空仪能够提供垂直方向上精确的水汽密度廓线图，因而将该探空站所在网格对应的层析垂直廓线与探空数据进行对比。本文采用平均绝对偏差(MAE)和均方根误差(RMS)作为评定层析方法的指标，共设计以下4种方案进行试验。

图4 有效信号数、空格率、边缘空格率统计Fig.4 Statistical results of effective signal number,space rate and edge space rate

表2 有效信号数、空格率、边缘空格率统计Tab.2 Statistical results of effective signal number,space rate and edge space rate

方案1：仅考虑层析区域天顶信号的层析模型(GNSS)。

方案2：顾及边界信号的层析模型(GNSS+)。

方案3：顾及边界信号但不考虑PWV拟合函数模型误差的层析模型(GNSS-)。

方案4：顾及边界信号但使用传统垂直约束的层析模型(TRA)。

为了对比分析4种方案的层析效果，本文首先计算了DOY 129—DOY 134每天UTC 12 h各种方案解算的无线电探空站所在位置不同高度上的水汽密度，并分别与无线电探空仪计算的水汽密度对比。各方案与无线电探空仪对比的垂直廓线结果如图5所示。

图5表明，各方案在层析时刻的水汽结果具有相同的变化趋势，与探空水汽结果的吻合程度均较高，但是各方案的反演结果在底层区域有所差异，使用传统垂直约束的TRA模型反演的底层水汽密度值易出现偏大的情况，这是由于层析区域底层水汽密度并非指数分布，采用按照指数递减规律的传统垂直约束方程在底层会导致较大偏差。为了对层析结果和探空数据计算结果进行直接对比，本文统计了DOY 129—DOY 135每天UTC 12 h各种方案的层析结果，利用层析结果计算出无线电探空站所在位置的水汽密度，分别与无线电探空仪数据计算的结果进行对比。图6为各方案统计的每天的MAE和RMS，表3为对应的统计结果。

图5 DOY 129—DOY 134(UTC 12)时的垂直水汽廓线分布Fig.5 Vertical water vapor profiles distribution at UTC 12 h in DOY 129—DOY 134

图6 4种方案(DOY 129—DOY 135)层析结果与无线电探空仪对比误差统计Fig.6 The statistical result of difference between different tomographic results and radiosonde for DOY 129—DOY 135

表3 4种方案(DOY 129—DOY 135)层析结果与无线电探空仪对比误差统计Tab.3 The statistical result of difference between different tomographic results and radiosonde for DOY 129—DOY 135 g/m3

图6表明，GNSS方案和TRA方案整体上层析结果较差，MAE和RMS偏大，GNSS+方案结果最优，GNSS-方案略低于GNSS+方案。由表4数据可知，GNSS+方案的层析结果平均MAE和RMS最优，分别为1.45 g/m3和1.82 g/m3，相对于GNSS-方案略有降低，分别为1.4%和1.6%；相对于GNSS方案和TRA方案有较大降低，前者分别降低9.4%和12.1%，后者分别降低5.2%和5.7%。这说明本文采用的顾及边界信号并考虑PWV拟合函数模型误差的GNSS+层析模型精度最高，相对于不考虑PWV拟合函数模型误差的GNSS-模型精度略有提升，相对于只考虑层顶信号的GNSS模型和使用传统垂直约束的TRA模型这两个传统的层析模型，精度存在较大提高。

为了比较不同方案下层析结果的精确性和可靠性，本文将DOY 129—DOY 135一周数据下4种方案解算的探空站上廓线水汽密度与无线电探空数据进行相关性分析，结果如图7所示。

图7 4种方案下层析结果与探空结果的线性相关程度Fig.7 Linear correlation between the results of tomography and the results of sounding information of four schemes

图7表明，GNSS+层析模型的相关系数最高，达到0.946。GNSS-方案相关系数略低于GNSS+。GNSS和TRA两种传统模型相关系数较低，GNSS最低。相关系数分析结果与上述MAE和RMS分析结果相一致，GNSS+方案结果最优，精度最高。此外，TRA模型在底层区域解算的水汽密度值(图中方框内区域)误差较大，这与TRA模型在底层区域水汽密度解算结果易偏大的特点相一致。

3 结 论

本文以香港CORS网2016年5月8日至2016年5月14日的观测数据和无线电探空数据，对融合边界信号、基于探空信息拟合函数建立垂直约束的多模层析方法进行了试验论证，详细分析了边界信号及垂直约束方法对水汽反演结果的优化效果。研究结果表明：边界信号的加入使得观测量提高了51.9%，空格率尤其是边缘和底层空格率大为下降，层析结果的平均MAE和RMS分别降低了9.4%、12.1%；通过搜索松弛因子削弱了PWVH拟合函数误差对层析结果的影响；垂直约束方法的改进改善了水汽反演结果，能够改善基于指数递减特性的传统垂直约束在底层区域反演结果易偏大的不足，平均MAE和RMS分别降低了5.2%、5.7%。本文采用的GNSS+方案的层析结果相较其他方案能够获得更高精度的三维水汽分布信息。