2017, Vol. 37
2. 武汉大学地球空间环境与大地测量教育部重点实验室，武汉市珞喻路129号，430079
常用的水汽层析算法可分为两大类：一类是迭代重构算法[1-4]，另一类是非迭代重构算法[5-9]。层析观测方程通常是秩亏的，故利用非迭代重构算法解算时需要添加额外的约束方程[4]。约束方程一般由经验大气分布信息给出[7]，因此估算出的水汽值精度有限。而利用迭代重构算法要事先对研究区域的所有网格赋予一个水汽密度初值，其反演精度很大程度上取决于所赋初值的精度。
鉴于迭代重构算法和非迭代重构算法存在的不足，本文提出组合重构算法。先将利用非迭代重构算法得到的解作为迭代初值，然后再通过迭代重构来逐步改善初值精度，从而提高所反演的水汽密度的精度。在实际试算过程中，我们在经典的对流层层析模型基础上提出一个新的层析模型，选用高斯指数模型作为非迭代重构算法的垂直约束，利用组合重构算法估算水汽三维分布信息。通过香港实测数据验证了新算法的可行性。
1 基于组合重构算法的对流层层析技术 1.1 新层析模型的确定Flores[5]提出基于像素基的水汽层析方法进行非迭代重构解算，假设每个层析格网中的水汽密度值均匀不变，单元网格内的斜路径水汽含量SWV表示为单一的水汽密度值与信号斜路径距离的乘积。但是实际情况并非如此，独立层析网格的水汽密度在水平和垂直方向都存在变化。为提高层析精度，作一新的假设，即每个独立网格8个顶点的水汽密度不同，通过这8个顶点的水汽密度值去估算出一个最佳的水汽密度值作为整个独立网格内的水汽密度值(图 1)。
|
图 1 独立层析网格 Fig. 1 Separate tomography voxel |
图 1展示了一个独立的层析网格。首先在水平方向用插值方法分别估算出p1和p2点的水汽密度值m(p1)和m(p2)，然后再内插出p点的水汽密度值m(p)，并把m(p)作为整个网格内的最佳水汽密度值。在水平方向，选用贝叶斯插值[10]分别估算m(p1)和m(p2):
|
(1) |
|
(2) |
式中，C(aik)和C(bik)分别表示第k个层析网格底部和顶部的第i个顶点的水平内插系数，其中i=1, 2, 3, 4；m为层析网格顶点的水汽密度向量；m(h(p1k))和m(h(p2k))中的h表示第k个层析网格的高度。
水汽密度值在垂直方向呈现高斯指数变化：
|
(3) |
式中，M为地表处水汽密度值，h为距地面的高度，hsc为水汽标高，Δh为水汽标高变化量。hsc和Δh可以通过气象产品求得。
利用高斯指数模型，可以最终估算出p点的水汽密度m(h(pk)):
|
(4) |
新的层析模型需要估算层析网格8个顶点的水汽密度值，因此会增加层析方程的参数个数。在实际的非迭代重构算法中，我们将层析网格上边界的水汽密度强制约束为一个定值。如果假设传统方法某一层待估参数的个数为nl×ml，而新的模型中该层的待估参数的个数为(nl+1)×(ml+1)，新方法并不会增加过多的待估参数，不会加重层析系数方程组的秩亏情况[10]。其中，nl和ml分别为层析格网在经度和纬度方向投影平面上待估参数的个数。
1.2 反演水汽密度的非迭代重构算法GPS信号的对流层延迟可以通过双差、单差、非差或精密的对流层模型来估算[11]。利用精密单点定位方法可求取天顶路径总延迟ZPD，根据测站气象文件(M文件)提供的气象产品，选用Saastamoinen对流层干延迟模型估算对应的天顶干延迟ZDD，进一步可得到天顶湿延迟ZWD，获得斜路径水汽含量SWV。用矩阵表示为：
|
(5) |
式中，Sq表示第q条GPS射线穿越网格的距离矩阵。受卫星星座及GPS测站几何分布以及层析时间步长的限制，部分层析网格没有观测值通过而造成层析观测矩阵秩亏。为了克服秩亏，需要在水汽层析反演中附加一定的约束条件。
水平方向上可以采用高斯距离加权函数进行约束[9]，在垂直方向附加先验信息值进行约束[10-11]。本文选用COSMIC掩星产品wetPrf计算的水汽密度值作为反演的先验值。此外，为了提高垂直方向的层析反演精度，选用式(3)的高斯指数模型作为垂直方向的约束条件。
将观测方程与附加约束条件联合解算：
|
(6) |
式中，B为水平约束方程系数矩阵，H为单位矩阵，L为水汽密度的先验值，V为垂直约束系数矩阵，Δi(i=1, 2, 3, 4)为误差阵。在实际估算中，上边界层水汽密度值强制约束为0.01 g/m3，各类观测方程之间采用等权进行解算。
将式(6)通过奇异值分解的方法求解未知的水汽密度m。然后通过式(4)估算每个格网内的水汽密度值m(h(pk))。
1.3 反演水汽密度的迭代重构算法本文对经典的代数重构算法进行改进，可进一步提高迭代收敛速度以及重构精度。改进后的迭代公式为[11]：
|
(7) |
式中，λ为松弛因子，其取值范围为0~1。经过试验，取改进的代数重构算法的最佳松弛因子为0.008。将式(6)和(4)解得的水汽密度值作为式(7)迭代的初始值，然后依|SWVq-S·m |=min作为迭代终止条件，其中min的确定方法是第k次和第(k+1)次迭代结果之差小于0.01 mm。
1.4 反演水汽密度的组合重构算法非迭代重构算法反演精度不高，而迭代重构算法的反演精度又很大程度上取决于初值的精度。基于此，本文提出联合使用迭代重构算法和非迭代重构算法进行水汽反演的组合重构算法。首先利用非迭代重构算法进行解算，由于存在观测噪声以及其他误差，所得结果只是一种近似解。再将反演结果作为迭代重构的初始值进行迭代，不仅可以进一步提高水汽密度值的整体反演精度，而且还克服了迭代重构算法在初始值选取方面对数值预报模型以及气象文件的依赖。
2 试验结果与分析为了检验提出的组合重构算法的有效性，选用香港地区12个地基GPS气象监测站2011-11-13及2012-02-10、02-17、03-20、04-16共5 d的数据进行计算。并选取5个与香港地区共址的COSMIC掩星数据反演结果作为水汽层析约束的先验值(表 1)。共址的标准为：水平经度和纬度的差距都小于1°。
|
表 1 与香港地区共址的COSMIC掩星事件 Tab. 1 COSMIC occultation events in Hong Kong |
COSMIC的wetPrf产品提供了0.1 km垂直分辨率的可计算大气水汽密度值的大气温度、气压、水汽压以及折射指数等大气产品。由下式求得大气水汽密度值：
|
(8) |
式中，mw为水汽密度值；Pw为水汽分压，单位为hPa；uw为水汽的摩尔质量，其值为18.02 kg/kmol；R是普适气体常数, 且R=8 314 Pa·m3·K-1·kmol-1；T为大气温度，单位为K。利用式(3)拟合水汽密度随高度变化的函数。限于篇幅，本文仅给出C005.2012.048.11.10.G20掩星事件的拟合结果。其拟合公式为：
|
(9) |
此外，从图 2可以看出，高斯指数拟合的水汽密度值与原始的水汽密度分布基本趋于一致。同时也可以看出，水汽主要集中在4 km以下，10 km以上高度的水汽密度值接近为0。受GPS掩星探测深度、近地面区域掩星产品精度有限及近地面区域水汽的逆增层分布等影响，使得近地面区域高斯指数拟合值出现较大偏差。通过对其他4 d的掩星产品进行分析，得到同样的结论，故我们确定香港地区水汽层层顶的高度为10 km。因此，我们的研究区域可以定位在香港地区0.5°×0.3°×10 km范围内。空间格网经度方向上的分辨率为0.062 5°，纬度方向上的分辨率为0.06°，垂直分辨率为500 m，东西和南北方向分别划分了8和5个网格，垂直方向分为20层，空间区域内网格总数为8×5×20=800个。所研究区域内的12个GPS测站从水平方向上看分布在不同的网格内(图 3)。
|
图 2 利用高斯指数模型拟合水汽密度 Fig. 2 Fitting water vapor using Gaussian exponential model |
|
图 3 水平方向网格划分及香港地区GPS气象监测站分布 Fig. 3 Horizontal grid division and GPS stations distribution in Hong Kong |
下载与掩星事件发生时间同一天的GPS观测数据及M气象数据，并选取掩星事件发生时刻前0.5 h的数据进行水汽层析反演。试验所用观测数据的时间间隔为5 s，经过精密单点定位处理后，每隔30 s输出一个结果。限于篇幅，仅给出2012-02-17及04-16的计算结果，如图 4和图 5所示。
|
图 4 2012-02-17不同经度区域层析水汽的垂直剖面 Fig. 4 Vertical profiles of tomography water vapor in different longitude regions on February 17, 2012 |
|
图 5 2012-04-16不同经度区域层析水汽的垂直剖面 Fig. 5 Vertical profiles of tomography water vapor in different longitude regions on April 16, 2012 |
可以看出，在高程方向上5 km以下的水汽密度值变化较大，5 km以上的水汽密度值变化趋于平稳。还可以看出，水汽在水平方向上分布比较均匀。为了与单个重构算法进行对比，分别统计了改进的代数重构算法(IART)、非迭代重构算法(NIRA)以及组合重构算法(CRA)反演水汽密度值的精度。其中，IART算法选用COSMIC掩星数据反演的水汽密度值作为初值。精度统计时，利用其中10个测站数据反演得到每个网格的水汽密度值，剩下的2个测站(HKMW，HKWS)由精密单点定位方法估计出SWV作为真值来检验反演获得的水汽密度值，统计结果如表 2(单位mm)所示。
|
|
表 2 SWV反演值与实测值差值的统计结果 Tab. 2 RMS statistical results for SWV differences between inversion and observation |
由表 2可以看出，组合重构算法(CRA)不仅优于非迭代算法而且优于改进的迭代算法，且由组合重构算法估算的SWV与真值之差的RMS值优于5.7 mm。此外，由于迭代算法的精度取决于迭代初始值的精度，因而精度最差。
为检验组合重构算法反演水汽分布的质量及新层析模型的优势，将计算结果(Tomo2)分别与Flores等[5]提出的层析模型计算结果(Tomo1)以及相应时段探空测站(114.16°E, 22.31°N)反演的水汽密度值(Radwv)进行比较。图 6仅给出2012-02-10及02-17两天的试算结果。
|
图 6 层析水汽值与探空水汽值的比较 Fig. 6 Comparisons of the water vapor values obtained from tomography and radiosonde data |
由图 6(a)和图 6(b)可以看出，Tomo2与Tomo1、Radwv的变化趋势基本一致，只是在地面到4 km高度之间出现较大差异。这种情况一方面是由于层析反演与探空数据在时间和空间上的差异造成的；另外一方面，层析反演中选择了COSMIC掩星事件反演得到的水汽密度值作为先验约束条件，但掩星剖面在底层的精度较差，从而造成层析反演与探空数据结果在低空位置出现差异。统计了4 km以下Tomo1与Radwv差值的RMS为1.67 g/m3；Tomo2与Radwv差值的RMS为1.23 g/m3；IART与Radwv差值的RMS为1.41 g/m3及NIRA与Radwv差值的RMS为1.53 g/m3。而在2 km以下，Tomo1与Radwv差值的RMS为3.02 g/m3；Tomo2与Radwv差值的RMS为2.76 g/m3；IART与Radwv差值的RMS为2.85 g/m3及NIRA与Radwv差值的RMS为2.94 g/m3。统计结果再次验证了新的层析模型相对于传统层析模型的优势以及组合重构算法相对于单一重构算法的优势。
3 结语针对单个重构算法在对流层层析水汽密度值反演中的局限性，提出一种组合重构水汽层析方法。该算法首先利用奇异值分解的方法(附加约束条件)求得水汽密度值的近似值，然后将其作为迭代算法的初始值，通过改进的代数重构算法进行迭代。数值分析表明，组合重构方法得到的水汽密度值优于单一重构算法反演的水汽密度值。此外，该方法能够为迭代重构算法提供高质量、高可靠性的初始值。通过与Radwv进行比较表明，Tomo2与其具有较好的一致性；又将Tomo2与Tomo1进行比较，验证了新层析模型的优越性，并验证了掩星产品作为先验信息值反演水汽三维分布的可行性。如何提高水汽在垂直方向上存在“逆增层”现象时的反演精度，是下一步将要研究的内容。
| [1] |
于胜杰, 柳林涛. 利用选权拟合法进行GPS水汽层析解算[J]. 武汉大学学报:信息科学版, 2012, 37(2): 183-186 (Yu Shengjie, Liu Lintao. Application of Fitting Method by Selection of the Parameter Weights on GPS Water Vapor Tomography[J]. Geomatics and Information Science of Wuhan University, 2012, 37(2): 183-186)
(  0) |
| [2] |
Jin S G, Luo O F, Park P. GPS Observations of the Ionospheric F2-Layer Behavior during the 20th November 2003 Geomagnetic Storm over South Korea[J]. Journal of Geodesy, 2008, 82(12): 883-892 DOI:10.1007/s00190-008-0217-x
( 0) |
| [3] |
Wen D B, Liu S Z, Tang P Y. Tomographic Reconstruction of Ionospheric Electron Density Based on Constrained Algebraic Reconstruction Technique[J]. Journal of GPS Solution, 2010, 14(4): 251-258
( 0) |
| [4] |
宋淑丽.地基GPS网对水汽三维分布的监测及其在气象学中的应用[D].上海: 中国科学院上海天文台, 2004 (Song Shuli. Sensing Three Dimensional Water Vapor Structure with Ground-Based GPS Network and the Application in Meteorology[D]. Shanghai: Shanghai Astronomical Observatory, CAS, 2004) http://cdmd.cnki.com.cn/Article/CDMD-80022-2004139591.htm
( 0) |
| [5] |
Flores A, Rius A, Rufffini G. 4D Tropospheric Tomography using GPS Slant Wet Delays[J]. Annales Geophysicae, 2000, 18: 223-224 DOI:10.1007/s00585-000-0223-7
( 0) |
| [6] |
Skone S, Hoyle V. Troposphere Modeling in a Regional GPS Network[J]. Journal of Global Positioning Systems, 2005, 4(1-2): 230-239
( 0) |
| [7] |
于胜杰, 柳林涛, 梁星辉. 约束条件对GPS水汽层析解算的影响分析[J]. 测绘学报, 2010, 39(5): 491-496 (Yu Shengjie, Liu Lintao, Liang Xinghui. Influence Analysis of Constraint Conditions on GPS Water Vapor Tomography[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(5): 491-496)
( 0) |
| [8] |
Riccardo N, Cucca M, Gabella M. Tomographic Reconstruction of Wet and Total Refractivity Fields from GNSS Receiver Networks[J]. Advances in Space Research, 2011, 47(5): 898-912 DOI:10.1016/j.asr.2010.12.025
( 0) |
| [9] |
Notarpietro R, Gabella M, Perona G. Tomographic Reconstruction of Neutral Atmospheres using Slant and Horizontal Wet Delays Achievable through the Processing of Signal Observed from Small GPS Networks[J]. European Journal of Remote Sensing, 2008, 40(2): 63-74
( 0) |
| [10] |
Perler D, Geiger A, Hurter F. 4D GPS Water Vapour Tomography: New Parameterized Approaches[J]. J Geophys, 2011, 85: 539-550
( 0) |
| [11] |
Xia P, Cai C, Liu Z. GNSS Troposphere Tomography Based on Two-Step Reconstructions using GPS Observations and COSMIC Profiles[J]. Ann Geophys, 2013, 31: 1-11 DOI:10.5194/angeo-31-1-2013
( 0) |
2. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, 129 Luoyu Road, Wuhan 430079, China