1 引 言
水汽是大气的重要组成部分，在众多天气现象中扮演着重要的角色。随着数值天气预报技术的发展，水汽三维分布信息作为数值天气预报初始场输入量以及各种天气过程的诊断分析的重要参考资料的重要性也日渐凸显，特别的，高空间分辨率的定量降水预报的预报结果对于温度场及湿度场尤其敏感[1, 2, 3]。目前，在实际应用中，数值天气预报湿度场信息主要取自无线电探空资料，该方法空间及时间分辨率都很低，难以满足实时获取大气水汽时空分布变化的要求。利用GPS探测大气水汽技术则可以快速、准确、全天候的连续获取大气水汽分布信息[4]，其中，地基GPS网层析水汽三维分布技术可以进一步获取水汽的垂直结构，具备其他探测技术所不具备的优势，因此国内外众多研究人员对此进行了大量的研究与验证，取得了一系列研究成果[5, 6, 7, 8, 9, 10, 11, 12, 13]。然而，众多研究在设计观测方程时一般假设单位网格内水汽分布均匀不变，这与实际情况显然差别较大，尤其是当网格边长较长时，对层析结果的精度会造成不利影响。
针对如上所述地基GPS网层析水汽技术中假设的局限性，文献[14]提出了基于数值积分的新参数化方法，摆脱了单元网格水汽均匀分布这一不合理假设。本文研究并发展了这一方法的基本思想，应用了更符合大气实际分布，运算量更小的内插计算方法。另外，传统方法中采用区域平面网格模型，而本文方法顾及了地球曲率的影响。
2 地基GPS网层析水汽三维分布方法 2.1 GPS倾斜路径湿延迟的获取GPS无线电信号在从卫星到达地面接收机的传播过程中会受到地球大气的延迟效应影响，大气延迟又可以分为电离层延迟与对流层延迟。电离层延迟可以通过GPS双频组合观测值加以消除，而对流层大气对于直至大约15 GHz的射电频率呈中性，信号传播产生非色散延迟[15]，即其折射量与信号频率无关。因此在高精度GPS数据处理中，对流层延迟必须作为未知参数加以估计。另外，考虑到实际大气的非球对称分布，可将表示大气在水平方向不均匀分布的水平梯度GNS、GEW作为未知参数加以估计,并假设GPS数据处理后的相位残差ξ主要由对流层延迟湿分量造成。因此，信号传播路径对流层湿折射量SWD可以表示为
式中，STD、SHD分别为传播路径对流层总延迟及其干分量，可以分别通过天顶方向对流层总延迟及其干分量与相应的映射函数计算得到，天顶方向对流层延迟干分量可以利用测站地面气压通过Saastamoinen模型计算得到[16]；ε、φ为卫星高度角与方位角；MH(ε)为水平梯度映射函数。
当利用GPS双差观测值进行数据处理时，由于观测值之间存在相关性，估计得到的对流层延迟为测站间相对量。为了得到对流层延迟绝对量，通常有两种方法：一种是加入精确的气象观测数据(如探空、水汽辐射计等)，以确定相对量与实际绝对量之间的偏差；另一种是解算时引入多个500 km以上的长距离测站以减少观测数据间的相关性[17]。本文采用的是第2种方法。
2.2 基于数值积分的层析原理与方法 2.2.1 观测方程的建立由对流层延迟的定义可以将GPS信号传播路径上的湿折射量表示成大气湿折射率Nw沿传播路径S的积分
式中
式中，nw为大气折射指数湿分量，其大小与水汽有关，可以通过水汽压e(单位为hPa)以及温度(单位为K)计算得到。
传统的层析算法都是将观测区域上空划分成若干规则形状的网格单元，将信号传播路径分割成若干通过不同网格的区段，并假设各网格内水汽分布均匀不变，进而将传播路径湿折射量离散表示为各网格湿折射量之和，网格湿折射量等于该网格湿折射率与相应传播路径长度的乘积。然而该假设显然与实际情况不符，特别是在垂直方向上，不同高度处湿折射率变化非常明显。针对这一局限，如图 1所示，本文将观测区域上空划分成若干层，将传播路径分割成若干段，对每一段传播路径的湿折射量利用Newton-Cotes 积分公式计算[14]，即式(3)可以表达为
式中，pi=a+i·b－a/4(i=1,2,3)。
|
| 图 1 层析网格模型示意图 Fig. 1 Voxel model of new tomographic approach |
通过Newton-Cotes公式，可以将一段传播路径湿折射量表示成5个离散点湿折射率与该路径长度的线性表达式，而各离散点湿折射率则可以由网格顶点的湿折射率内插计算得到。大气中水汽分布在水平和垂直方向上有着不同的特性，因此内插时应分别在水平方向与垂线方向上进行。首先，将离散点沿垂线方向投影到各层层面上，依据Barnes插值公式[18]，将各层投影点湿折射率表示成所在网格的顶点湿折射率的线性表达式
式中，pij为pi点在第j层的投影点，j=1,2,…,n； Nw(xmj)(m=1，2，3，4)为第j层pij所在网格4个顶点的湿折射率，ωmj为相应的权重系数；X为所有网格顶点湿折射率向量，Aj为相应的系数矩阵
式中，rmj为网格顶点至pij点的距离；E为经验常数，一般取为0.25。
然后，在垂线方向上，离散点湿折射率Nw(pi)可以由上述各投影点湿折射率Nw(pij)内插得到。边界条件为两端点处二阶导数值为0的三次样条插值方法在文献[14]中已有论述，该方法较为复杂，当观测值较多时会显著增加矩阵运算量，因此本文提出另外一种简化的指数模型内插的方法。依据对流层水汽分布在垂线方向上基本呈现指数分布的特性，位于第j层与j+1层之间的pi点湿折射率可通过其在j、j+1层的投影点的湿折射率Nw(pij)、Nw(pij+1)内插得到
式中，hj表示第j层层面高程； hsc为湿折射率标称高度，可以通过统计当地探空数据获得。
最后，得到各离散点湿折射率与X的关系式，通过式(4)可以得到各层传播路径湿折射量与X间转化矩阵，将所有转化矩阵叠加即为GPS传播路径总湿折射量转换矩阵，若t时刻有k个传播路径湿折射量观测值由网格模型顶部刺入网格，则
式中，Δ(t)为观测值噪声向量。
过去的研究在划分网格过程中近似地认为地面为水平面，建立的网格模型为平面围成的立方体。而事实上地球是一个近似椭球体，因此，同一点在区域平面网格模型中的高程与其实际高程并不一致，这就会导致在相同的对流层顶层高度设定下，相同的观测路径在区域平面网格模型中和实际穿过对流层的路径长度并不一致，两者差值会随着路径高度角的降低而增大，因此在层析中有必要考虑地球曲率以避免上述因素造成的误差。另外，使用顾及地球曲率的网格模型也可以降低实际观测方程中零系数参数的比例。在地面测站难以实现高密度布设的条件下，受卫星星座的空间分布以及观测值截止高度角的影响，单纯的增加观测卫星数量(例如增加GLONASS、北斗卫星等)实际上难以实现上述目的，但却可以减弱直至消除实际观测量组成的观测方程的不适定性。在多种导航卫星系统同时发展，可观测卫星数量将大大增加的情况下，减少零系数参数比例对于层析来说是有利的，能够使得层析的结果反映更多区域的实际水汽分布情况。
为计算方便，假设地球是半径为R的球体，由信号传播路径和地心O可以确定出一个地球大圆横切面，如图 2所示，a点为测站相位中心点，b点为传播路径与模型一层面的交点，ha、hb分别为a、b两点的高程，ε为传播路径高度角，易知
|
| 图 2 过地心及GPS信号传播路径的地球横切面图 Fig. 2 Cross-section of earth with earth center and a GPS signal’s propagation path on it |
求得Sab后，可以进一步得到各离散点Pi到站心a的距离Rpi，结合传播路径高度角ε、方位角φ组成以a点为原点的站心极坐标，通过坐标转换求得离散点在地心空间直角坐标系中的坐标
式中，Ta－1为相应的坐标旋转矩阵。再进一步将得到的pi点地心空间直角坐标转换成为大地坐标，本文不再赘述。
2.2.2 约束条件由于GPS星座及地基GPS观测站空间分布和数量的限制，观测方程在多数情况下是不适定的 ，因此需要加入额外的约束方程以解决这一问题。本文采用3种约束条件：水平约束、垂直约束和边界值约束。
(1) 水平约束。即某一网格顶点的湿折射率值可以由同一层面与其相邻的其他4个顶点湿折射率按照式(5)计算得到。
(2) 垂直约束。可以由统计信息得到下列公式对垂直方向相邻两层湿折射率进行垂直约束[19]
各参数意义如前所述。
(3) 边界值约束。即认为网格模型最顶层(高度=15 km)湿折射率为0。
以上约束将和实际观测方程共同组成新的最终观测方程。
2.2.3 利用Kalman滤波估计参数假设在短时间内各网格的湿折射率符合高斯-马尔科夫的随机游走平稳过程，即
w(t)为状态噪声。由式(8)及约束条件方程组成的观测方程以及如式(12)的状态方程，通过Kalman滤波有
在t时刻，H(t)为观测方程系数矩阵；DΔ(t)为观测值方差协方差阵；Q(t)为过程噪声方差协方差阵。由式(13)即可求出所有网格顶点湿折射率及其方差。
3 实测数据试验 3.1 数据来源及数据处理策略为了验证新方法的有效性与精度，本文以美国Dallas市及其周边地区为试验区域，通过处理该区域内CORS实测数据，在此基础上层析得到了目标地区上空水汽垂直结构。试验区域共计有17个GPS站以及相应的地面气象观测站，测站间距40～80 km，GPS观测数据采样率取30 s，地面气象观测站每隔一小时左右观测记录一次。垂直方向网格划分采用(300 m，600 m，900 m，1300 m，1900 m，…)的非均匀分层。水平网格划分及各测站空间分布如图 3所示。
|
| 图 3 GPS网测站分布及网格的划分 Fig. 3 Distribution of GPS stations and tomographic voxels |
本文利用GAMIT高精度GPS数据处理软件处理实测数据，由于GAMIT处理的是GPS双差观测值，因此本文引入了AMC2、BLYT、INEG、PIE1 4个远距离IGS跟踪站以削弱观测资料之间的相关性。天顶延迟和水平梯度参数时间分辨率分别为1 h和2 h；投影函数为GMF；截止高度角取10°；再结合测站地面气象观测数据及GPS数据后验残差，本文每隔5 min由2.1节所述方法得到所有有效传播路径的湿折射量，作为Kalman滤波的观测值；经过多次试验，本文取实际观测方程中非零系数参数对应的水平约束方差为602，而零系数参数水平约束方差为102；垂直约束方差均取为602。
为了获得较好的Kalman滤波初值，本文获取了试验时间段内FWD探空站所有探空数据，湿折射率并不是探空的直接观测值，需要通过式(3)计算得到，结果如图 4所示。对得到的所有湿折射率观测值依据高度进行指数拟合，状态方程的Q(t)可以通过分析各高度层探空结果时间序列计算得到[20]，滤波状态初值如式(14)所示
式中，X0(h)、DX0(h)为待求点初始状态及其方差，h为待求点高度，单位为 km。
|
| 图 4 FWD探空站2010年6月份探空湿折射率垂直轮廓 Fig. 4 Wet refractivity vertical profile from radio soundings of FWD RS station in June，2010 |
层析过程起始于UTC 2012年第184天00时，结束于UTC 2012年第193天23时。从起始时刻开始，FWD探空站(32.83°N，262.70°E)间隔12 h释放一次探空气球，共进行了20次观测。
前文提到，当网格模型顾及地球曲率影响时，将会减少实际观测方程中零系数参数的比例，图 5统计了在相同的观测值与相同的网格划分情况下，新层析方法和传统方法零系数参数比例的时间序列，可以看到，两者间有着较强的相关性，但是传统层析方法的比例明显高于新方法。零系数参数越少，意味着更多的参数是由实际观测方程而非约束方程计算得到，在有限的地面测站分布情形下，增加观测卫星数量的方法对于新方法的网格模型来说更加有效，这也从一个方面反映了新方法的优越性。
|
| 图 5 各历元空白网格比例 Fig. 5 Percentage of zero coefficient parameters at each epoch |
以层析时间段内的20次无线电探空值为真值，3种层析方法结果的误差分布区间及其RMS统计如图 6所示，3种层析方法的误差在-1×10-6～3×10-6之间分布最集中；与此同时，按图 6顺序从左至右的3种方法分别有86.83%、87.29%、87.37%的绝对误差小于10×10-6，比例依次递增；3种层析方法结果RMS分别为9.622 2×10-6、9.524 9×10-6和 9.518 7×10-6，指数垂直内插法略小于三次样条内插法，两者都小于传统方法；笔者发现最底层(即高度小于300 m)的层析结果误差相比其他层偏大，如果除去最底层的数据进行统计，则3种方法的RMS分别变为8.171 5×10-6、8.008 3×10-6、7.993 4×10-6，都明显小于前述结果，但是三者之间的大小关系并未改变，最底层误差相对较大可能与测站站间距较大有关，另外近地空间干扰因素较多，探空数据也可能存在较大误差。以上统计信息都说明新的层析方法的统计精度要高于传统方法。图 7给出了其中4个时刻无线电探空得到的湿折射率垂直轮廓线与3种层析方法得到的距探空站最近的垂直轮廓线(32.877 6°N，262.736 6°E)的比较图。可以看到，3种方法的层析结果大体上都与探空结果较为吻合，而两种新的层析方法结果则较为接近。
|
| 图 6 3种层析方法误差分布区间统计 Fig. 6 Distribution of wet refractivity results’ errors by three tomographic methods |
|
| 图 7 探空湿折射率垂直轮廓与层析湿折射率垂直轮廓线的比较 Fig. 7 Comparison of radio soundings (circular dots) and three tomographic methods’ results |
利用美国Dallas市地区CORS网GPS站及其相应的气象站观测数据，本文验证分析了顾及地球曲率的GPS层析水汽垂直轮廓线数值积分方法的有效性及精度。结果表明：顾及地球曲率的数值积分方法避免了传统层析方法的不合理假设，有效降低了实际观测方程中零系数参数的比例；通过与探空水汽垂直轮廓线的比较，表明数值积分方法结果的RMS值以及误差分布均优于传统层析方法；不同的垂直内插公式对于数值积分方法的层析结果也有影响，简单的垂直指数内插公式结果要略优于复杂的三次样条内插公式。
本文在研究中也发现若干问题有待解决，试验虽然采用了较高精度的垂直约束，但是仍然难免部分层析值与探空值相差较大，下一步应该研究如何突破GPS观测站分布的限制，引入其他方法(如雷达、数值预报产品等)垂直廓线信息作为约束以提高GPS层析垂直轮廓线的精度。另外，本文在评价层析精度时以探空结果为真值，但事实上探空数据本身也存在误差，如何更好地评价层析精度有待进一步研究与讨论。
致谢：感谢美国MIT授权使用GAMIT/GLOBK软件，感谢NOAA提供的试验数据
| [1] | DUCROCQ V,RICARD D,LAFORE J P,et al.Storm-scale Numerical Rainfall Prediction for Five Precipitating Events over France:on the Importance of the Initial Humidity Field [J].Weather and Forecasting,2002,17(6):1236-1256. |
| [2] | PARK S K.Nonlinearity and Predictability of Convective Rainfall Associated with Water Vapour Perturbations in a Numerically Simulated Storm[J].Journal of Geophysical Research. (Atmospheres),1999, 104:31575-31587. |
| [3] | BASTIN Sophie,CHAMPOLLION Cedric,BOCK Olivier,et al.On the Use of GPS Tomography to Investigate Water Vapor Variability during a Mistral/sea Breeze Event in Southeastern France [J].Geophysical.Research Letters,2005,32(5):15787-15801. |
| [4] | BEVIS M,BUSINGER S,HERRING T,et al.GPS Meteorology:Remote Sensing of Atmospheric Water Vapor Using the Global Positioning System [J]. Journal of Geophysical Research,1992 ,97 (14):15787-15801. |
| [5] | FLORES A,RUFFINI G,RIUS A.4D Tropospheric Tomography Using GPS Slant Wet Delays [J].Annales Geophysicae,2000,18(2):223-234. |
| [6] | FLORES A,J VILà-GUERAU DE ARELLANO,GRADINARSKY L P, et al.Tomography of the Lower Troposphere Using a Small Dense Network of GPS Receivers [J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(2):439-447. |
| [7] | NILSSON T,GRADINARSKY L.Water Vapor Tomography Using GPS Phase Observations:Simulation Results[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(10):2927-2941. |
| [8] | TROLLER M,GEIGER A,BROCKMANN E,et al.Tomographic Determination of the Spatial Distribution of Water Vapor Using GPS Observations [J].Advances in Space Research,2006,37:2211-2217. |
| [9] | LIU Haixia,XUE Ming,PURSER R J,et al.Retrieval of Moisture from Simulated GPS Slant-path Water Vapor Observations Using 3DVAR with Anisotropic Recursive Filters [J].Monthly Weather Review,2007,135(4):1506-1521. |
| [10] | BENDER M,GALINA D,GE Maorong,et al.Development of a GNSS Water Vapour Tomography System Using Algebraic Reconstruction Techniques [J].Advances in Space Research,2011,47(10):1704-1720. |
| [11] | NOTARPIETRO R,CUCCA M,MARCO G,et al.Tomographic Reconstruction of Wet and Total Refractivity Fields from GNSS Receiver Networks[J].Advances in Space Research,2011,47(5):898-912. |
| [12] | 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.(宋淑丽.地基GPS网对水汽三维分布的监测及其在气象学中的应用[D].上海:中国科学院上海天文台,2004.) |
| [13] | ZHANG Shuangcheng.Research on Technology and Application of Remote Sensing Water Vapor Using Ground-Based GPS/Met [D].Wuhan:Wuhan University,2009.(张双成.地基GPS遥感水汽空间分布技术及其应用的研究[D].武汉:武汉大学,2009.) |
| [14] | DONAT P,ALAIN G,FABIAN H.4D GPS Water Vapor Tomography: New Parameterized Approaches [J].Journal of Geodesy,2011,85:539-550. |
| [15] | WEI Ziqing,GE Maorong. GPS Relative Positioning Model [M].Beijing:Surveying and Mapping Press,1998.(魏子卿,葛茂荣.GPS相对定位的数学模型 [M].北京:测绘出版社.1998.) |
| [16] | CHEN Yongqi,LIU Yanxiong,WANG Xiaoya, et al.GPS Real-time Estimation of Precipitable Water Vapor:HongKong Experience[J].Acta Geodaetica et Cartographica Sinica,2007,36(1): 9-12.(陈永奇,刘焱雄,王晓亚,等.香港实时GPS水汽监测系统的若干关键技术[J].测绘学报,2007,36(1):9-12.) |
| [17] | DUAN Jingping, BEVIS M, FANG Peng, et al.GPS Meteorology:Direct Estimation of the Absolute Value of Precipitable Water Vapor [J].Journal of Applied Meteorology,1996,35:830-838. |
| [18] | SHOU Shaowen.Synoptic Meteorology [M].Beijing:China Meteorological Press,2008.(寿绍文.天气学 [M].北京:气象出版社,2008.) |
| [19] | 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.(于胜杰,柳林涛,梁星辉.约束条件对GPS水汽层析解算的影响分析[J].测绘学报, 2010,39(5):491-496.) |
| [20] | GRADINARSKY L P,JARLEMARK P.Ground-Based GPS Tomography of Water Vapor:Analysis of Simulated and Real Data[J].Journal of the Meteorological Society of Japan,2004,82(1B):551-560. |