气象学报  2013, Vol. Issue (2): 318-331   PDF    
http://dx.doi.org/10.11676/qxxb2013.026
中国气象学会主办。
0

文章信息

万蓉, 郑国光, 于胜杰, 付志康 . 2013.
WAN Rong, ZHENG Guoguang, YU Shengjie, FU Zhikang . 2013.
基于观测约束的地基犌犘犁三维水汽层析技术研究
A study of the ground-based GPS 3D water vapor tomography with radiosonde vertical constraining
气象学报, (2): 318-331
Acta Meteorologica Sinica, (2): 318-331.
http://dx.doi.org/10.11676/qxxb2013.026

文章历史

收稿日期:2012-07-11
改回日期:2012-12-14
基于观测约束的地基犌犘犁三维水汽层析技术研究
万蓉1,2, 郑国光1,3, 于胜杰2, 付志康2    
1. 南京信息工程大学, 南京, 210044;
2. 中国气象局武汉暴雨研究所暴雨监测预警湖北省重点实验室, 武汉, 430074;
3. 中国气象局, 北京, 100081
摘要:全球定位系统(GPS)卫星信号穿过大气层时发生的偏折和延迟, 可以用来反演信号传播路径上的大气水汽总量。为 获取区域高精度的大气水汽三维分布, 借助分布密集的地基GPS观测网及其斜路径水汽观测, 建立新的观测约束层析模型, 提出以高斯函数为水平约束, 区域逐月多年探空观测为垂直约束, 即以平均量为先验值, 以标准偏差为权重矩阵的新方法;并 在层析算法中引入地面观测, 以提高整体尤其是低层反演精度。三维水汽层析网格模型基于长江中游鄂东区域的22站地基 GPS加密网搭建, 实时解算系统可逐时输出三维水汽产品。三维湿折射度和水汽密度可以分别由斜路径的湿延迟总量和水 汽总量观测反演获得。以2010年开展的汛期联合加密探空观测为参照, 对三维层析的总体反演精度、低层反演精度、层析区 域中心与边缘反演精度进行了对比和分析。结果显示:整体样本检验三维水汽密度平均偏差为-0.63g/m3,标准偏差为 1.22g/m3,与探空相关系数为0.98;水汽密度与探空资料的相关比湿折射度与探空资料的相关好;对于不同层析区域, 中心 区域观测元数量较为丰富, 使得位于中心的层析精度好于整体和边缘;加入低层观测的层析结果与探空的相关比未加低层观 测时的好, 低层观测的加入提高了层析与探空的相关, 减小了低层层析标准差、区域中心和2km 以上层析的均差, 有效提高 了反演精度;低层观测的加入对整体标准差的影响, 可能与加剧观测方程中长度矩阵元素间的量级差异有关。
关键词全球定位系统     三维水汽层析     观测约束     水汽密度    
A study of the ground-based GPS 3D water vapor tomography with radiosonde vertical constraining
WAN Rong1,2, ZHENG Guoguang1,3, YU Shengjie2, FU Zhikang2    
1. Nanjing University of Information & Science, Nanjing 210044, China;
2. Hubei Key Laboratory for Heavy Rain Monitoring and Warning Research, Institute of Heavy Rain, China Meteoeological Administration, Wuhan 430074, China;
3. China Meteoeological Administration, Beijing 100081, China
Abstract:The atmospheric water vapor can be retrieved exactly on certain scales by the deflection and delay of the Navigation Satcllite Timing And Ranging Global Position System(GPS) satellite signal along signal path. Based on the water vapor obscr-vations along the signal paths over the density ground-based GPS network,a new GPS tomography model has been established by using a Gauss function as horizontal constraint and the radiosondc(RS) observation vertical constraint, which includes RS mean values of 4 years as priori information and RS standard deviation as the weight matrix to get the regional high-accuracy three-dimensional (3D) distribution of the atmospheric water vapor. The observations on the ground surface have been used in the tomography formulation to improve the tomographic result accuracy in the overall region, especially the low-level. A 3D water vapor tomographic grid has been designed based on the GPS observation network with 22 GPS sites in the eastern part of Hubei Province in the middle reaches of the Yangtze River. The quasi-real-time retrieval system can provide the hourly output of 3D water vapor density. The 3D wet refraction and water vapor density can be retrieved from the observations of Slant Wet Delay(SWD) and Slant Water Vapor(SWV), respectively. The accuracy of the retrieval results in the 3D tomography domain of the overall, the low-level, the center and the edge has been analyzed by comparing to the RS observations during the flood season of 2010. The results show that for overall sample tests the 3D water vapor density average deviation is -0.63 g/m3, the standard deviation is 1.22 g/m3, and the correlation coefficient with RS is 0.98. The correlation between the retrieval results for water vapor density and soundings is better than that between those for the wet refraction and RS ohservations, and the correlation hctween RS ohscrvations and the results of water vapor density after joining the ground surface observations is better than those without joining the ground surface observations in the tomography. For the rich observations path of GPS, the accuracy at the center of the tomography domain is better than that in the overall and at the marginal domain. Combined with the ground surface observations, the 3D water vapor retrieval accuracy is improved with enhancing the correlation between tomographic results and RS observations, reducing the standards deviation of tomographic results in low-level of the region, and the average deviation of tomographic results in the center and the layer above 2 km of the region. The magnitude order difference between the elements in the signal length matrix due to the ground surface observation being used in tomography could disturb the stability of the observation equations and thus increase the overall standard deviation.
Key words: GPS     3D Water vapor tomography     Obscrvations constraining     Water vapor density    
1 引 言

水汽约占大气总体积的4%,对地-气系统径向辐射能量平衡、大气的垂直稳定度、云的形成和降水的形成及演变有显著影响。精确的水汽观测数据是有效预报中小尺度灾害性天气(暴雨、雷雨、冰雹、暴雪、龙卷风等)的基本保证。利用全球定位系统(GPS)卫星信号穿过大气层时发生偏折和延迟的现象,可以在一定程度上精确反演大气水汽参数,从20世纪90年代,国际上有学者探讨利用GPS数据估计大气可降水的可行性至今,直接从GPS观测数据中反演出来的大气可降水量精度已达毫米级(Duan et al,1996毛节泰,1993),并被用于数值预报和天气分析研究(袁招洪,2005; BeHaan et al,2003)。GPS气象学正从理论研究步入应用阶段,人们已不满足于从GPS中得到整层水汽总量,更需要获得精确的三维水汽分布。随着中外GPS网建设发展,获取三维水汽的条件已逐渐成熟。

层析也称断层扫描技术是目前从GPS网反演三维水汽的主要方法,即利用一定研究区域内不同方向和位置的观测值积分来反演各部分信息,该技术最初于20世纪50年代被美国医学界采用,后来广泛用于医学、地质、物探等领域。GPS斜路径观测方向上的水汽总量是水汽在观测方向上的积分,包含着水汽在垂直和水平剖面分布的信息,基于多个斜路径水汽量,有可能通过层析方法获得水汽的三维分布。层析技术首先通过网格划分使观测路径上各部分信息的积分问题离散化,假定网格内的水汽分布是均匀的,可由卫星信号所经网格长度和网格内的水汽参量,得到每条信号路径上的水汽积分总量而建立观测方程组,求得方程组的解即得到三维水汽参量的值。由于网格模型上空GPS卫星在天空中的分布不均,而当地面的接收站间距大于网格尺度时,就会出现许多网格通过信号较少,甚至没有信号通过;即使站点足够密集,对于层析区域边界的网格也存在上述现象。这导致观测方程个数少于未知量,出现观测方程组的无定解问题。

为克服观测方程的无定解问题,中外采用了多种方法:利用数值预报结果作为背景场,提供先验值和垂直方向的约束得到观测方程的唯一解;或直接将其他辅助观测作为先验信息,以附加观测方程形式与GPS观测方程一起组成层析模型进行解算。国际上,Flores等(2000)率先使用层析技术得到区域四维湿折射度,Macdonald等(2002)利用间距为40 km的GPS观测网,结合地面微波辐射计水汽观测资料和探空站网的湿度资料,以三维变分(3DVAR)技术为基础,建立了应用GPS斜路径累计水汽总量分析三维水汽分布的系统,Gradinarsky等(2002)曹云昌等(2006)、张双成

等(2008)利用卡尔曼滤波技术解决弱几何站点分布的水汽层析问题,得到三维水汽信息。Gradinarsky等(2002)方法对初值具有一定的依赖性,卡尔曼滤波算法中的噪声误差也较难给定。因为是独立于观测之外的数学计算方法,要提高解算精度需要寻找更好的方法。有必要在三维层析中加入先验信息提高反演精度(Bi et al,2006)。附加约束的三维水汽层析法需要提供较为合理的约束条件,一般分为水平和垂直约束,约束先验值主要来自大气特性和观测,在算法上较为遵循大气的特性,有一定的合理性;MacDonald(2002)利用微波辐射计、Braun等(2003)利用激光雷达、Skone等(2004)利用单个探空观测等作为先验值用于层析算法的研究。宋淑丽等(2005)在层析中采用高斯加权函数进行水平约束,垂直方向上采用随高度变化的非等权约束方法,获得上海地区水汽三维分布;Bi等(2006)利用3 d的探空平均作为垂直先验信息获得三维层析结果。

利用2006—2009年探空辅助观测的逐月平均探空垂直廓线作为垂直约束先验值,用其标准偏差作为垂直约束的权重矩阵,利用高斯函数建立大气水平约束的三维水汽层析模型,解决三维水汽层析观测方程组的不适定性问题;由加密GPS网覆盖区域各站斜路径湿延迟量和斜路径水汽总量建立各自的观测方程组,分别得到三维大气水汽湿折射度和水汽密度结构;并在约束模型基础上,首次引入实时的地面观测,增加了低层网格的观测信息。将鄂东区域加密观测网的三维水汽层析结果与探空进行比对显示,结果表明,探空观测的垂直约束和地面观测加入后的层析算法改善了区域整体反演精度,尤其是低层网格水汽的反演精度。2 GPS资料处理与三维层析网格设计2.1 GPS资料斜路径水汽参量解算

要开展GPS三维水汽层析,必须先获取各个站点接收机与不同卫星建立的信号路径上的观测总量,即斜路径湿延迟总量(SWD,Slant Wet Delay)和斜路径水汽总量(SWV,Slant Water Vapor)。在卫星信号传播的倾斜路径上,总延迟量是干延迟量和湿延迟量的共同效应,穿过对流层的倾斜路径高度角越低路线会越长,延迟量也会越大,通过与高度角有关的映射函数将倾斜路径的各类大气延迟变换为等效的天顶总延迟,其中,天顶总湿延迟(ZWD,Zenith Wet Delay)代表水汽的均匀态(单位:mm)。而实际大气水汽分布并不均匀,除了映射函数项,还要充分考虑大气水平梯度项、湿残差项等因素。基于美国麻省理工学院研制的GAMIT软件(Release10.4,2010),对流层大气中GPS信号的斜路径湿延迟量SWD可以表示为

其中,ø为方位角,单位:°,θ为仰角,单位:°;右边第1项为等效天顶湿延迟项,mw(θ)为湿映射函数,无量纲;第2项为湿延迟水平梯度项,GNwGEw分别是南北、东西方向湿梯度分量(单位:mm),mΔ(θ)是水平梯度映射函数,无量纲;第3项Re为残差项(单位:mm)。在GAMIT软件解算参数中,选用全球映射函数(Boehm et al,2006),水平梯度映射函数表示为

延迟的水平梯度项包括大气干延迟梯度和大气湿延迟梯度。大气干延迟梯度主要是由气压和温度的水平分布梯度引起的,具有100 km左右的空间尺度和几天的时间尺度,而大气湿延迟主要是由水汽的水平分布梯度引起的,具有较小的空间尺度(<10 km)和几小时的时间尺度,但在总延迟梯度项中占主要的地位。采用平均处理把干大气延迟梯度和湿大气延迟梯度分离出来,即基于大气干延迟梯度具有较大的时间和空间尺度的特点,可假定在12 h内为一个常数,可将软件解算出来的各时次的梯度项在这一段时间内的平均值作为干大气延迟梯度,而各个时次的梯度项减去干大气延迟梯度即为大气湿延迟梯度。

Re延迟残差项为未被模型化的水汽的各向异性部分。在GPS数据处理中,通常基于双差观测值来构建观测方程,这样可以消除卫星和接收机钟差的影响,并能将双差模糊度固定为整数值,若在该过程中,所有的误差均被很好地模型化,数据处理后生成的残差则主要包括未被模型化的大气延迟部分。可通过GAMIT数据处理后得到相关的残差信息。

通过湿延迟与水汽总量的比例关系,可以得到斜路径的水汽总量

其中,PWV为天顶方向水汽总量,也称可降水总量,(单位mm);Π为本地化后大气平均温度的函数,无量纲(徐桂荣等,2009)。

上述方法计算得到的SWV与微波辐射计WVR观测结果比较,精度达到毫米级。仰角30°—90°观测数据对比,拟合度为0.88,相关系数为0.94,平均偏差为7.03 mm,偏差的标准差为5.55 mm。50°—90°仰角范围的对比显示平均偏差为0.8 mm,偏差的标准差为3.6 mm,相关系数为0.97。2.2 鄂东GPS加密网的层析网格设计

层析模型网格水平大小最理想是25—37 km。由于GPS观测网的几何结构会对层析结果产生一定影响,水平网格的划分应使站点尽可能均匀分布在不同网格中;为消除多路径效应等的影响,一般使用仰角15°以上的信号。假设在对流层中层析的高度为10 km,那么可以计算出GPS站所处的水平网格最大边长74.6 km;再考虑到如果要求区域中心大多数网格至少拥有2—3部GPS接收机的信号,因此,每个网格边长不宜小于25 km(图 1)。如果网格边长超过这个限制,解算出的水汽参量空间分辨率较低,如果小于这个限制,每个网格中可以穿过的GPS信号太少或路径比较短,由此可能导致构造的层析观测方程组的观测系数矩阵稀疏,甚至观测方程少于待求参量,造成求解不适定问题。在垂直方向上,每层厚度取为1 km,如果厚度太小信号噪音会影响层析结果(Bi,2006)。

图 1 GPS层析网格边长计算示意图 Fig. 1 Schematic diagram of the suitable size of the horizontal grid in GPS tomography

基于上述原则,在鄂东22站GPS加密网区域构建水汽三维分布模型(图 2)。层析网格模型坐标系统采用局部大地坐标系,区域范围:(29.5°—31.1°N,113.0°—116.0°E),高度0—10 km。将0—10 km 高度范围平均分为10层,每层1 km。水平方向上沿东西向划分为7个网格,间距约30 km;南北向划分为4个网格,间距约30 km(图 3a),模型区域对流层大气将被分为280个网格或六面体,最低一层网格输出值代表网格中心高度500 m,最顶层网格输出值代表网格中心9500 m。

图 2 长江中游暴雨外场观测试验基地的高密度水汽监测网(方框为层析区域) Fig. 2 Distribution of the GPS/MET network in the middle reach of the Yangtze River(GPS tomographic region marked with box is in eastern Hubei Province)
图 3(a)层析区域水平网格与站点分布,(b)斜路径信号穿越层析网格示意图 Fig. 3 Distribution of the GPS sites and the horizontal grids in the tomographic region(a),and schematic diagram of the signals in slant way passing through the tomographic grids(b)
3 SWDSWV的观测方程组建立

以GPS斜路径湿延迟SWD为观测值,给出层析湿折射度网格模型的构建方法。设GPS斜路径湿延迟观测值是沿GPS信号方向上湿折射率的积分(Gradinarsky,2002)

其中,nw为湿折射率(无量纲),dl为GPS信号斜路径(单位:km)。对流层的大气湿折射度表示成Nw=106(nw-1)(单位:mm/km),式(4)可改写如下:

基于斜路径的湿延迟总和与折射度是积分关系,通过建立和求解多个斜路径上延迟量积分观测方程组,可以获取湿折射度Nw的三维分布。大气湿折射度是水汽压、温度的函数,一定程度反映大气水汽特征(Rüeger,2002)

其中,e为水汽压(单位:hPa);T为温度(单位:K);大气物理常数k2=71.2952 K/hPa,k3=375463 K2/hPa。

若将获取的三维层析湿折射度转换为水汽密度,尚需要同一网格中的温度等气象参数,似乎不太可能。可建立基于SWV的观测方程,直接获取水汽密度。将斜路径方向上的水汽总量定义为该方向上积分水汽量的等效液水高度,观测方程为

其中,SWV代表倾斜路径方向上测站与卫星之间的水汽总量,ρ为液水密度(单位:g/m3),ρw为水汽密度(单位:g/m3)。将每一条观测路径上的水汽总量分配到这条射线所穿越的每个网格中,斜路径水汽总量与网格内的分段水汽关系为

由此建立基于SWV观测的层析模型,通过层析算法求解,获取斜路径上的逐个网格的水汽密度ρw,更便于天气研究。

对划分的网格逐一编号,设在一定时段内(如1 h)每个网格的湿折射度(或水汽密度)是一未知常数且分布均匀,则式(4)(或式(7))可以简化为

式中,SWVjWD为斜路径j上的湿延迟(或水汽总量),j=1,……,m,m为斜路径湿延迟(或水汽总量)参数的个数,sji为斜路径j在网格i中的长度,xρsNwi是处于网格i中的湿折射度(或水汽密度),i=1,……,n,n为划分的总网格数。SWD(或SWV)等于GPS信号射线穿过各网格的未知的湿折射度(或水汽密度)的总和(图 3b)。所以通常只选用从网格模型顶层穿出的观测值。在1 h内,根据穿过网格模型的多个斜路径观测值由式(9)可得到观测方程组
式中,m为斜路径湿延迟(或水汽总量)参数的个数,n为划分的总网格数。某一时刻整个GPS网的斜路径湿延迟(或斜路径水汽总量)可以合并为线性方程
式中,SWVWD为观测值斜路径湿延迟(或斜路径水汽总量)矩阵,包含一定时间段内所有的斜路径参数组成的观测值向量,S为穿过网格长度的设计矩阵,其第j行第i列的元素表示第j条射线穿过第i个网格的长度;XρsNw为未知参数湿折射度(或平均水汽密度数)的列矩阵。

由于接收机的数目有限,且它们在层析区域内分布不均匀,不能保证每一个网格中都有斜路径射线通过,致使线性方程中有一些观测元不会出现,直接导致了水汽层析观测方程的不适定性,无法获得三维湿折射度或水汽密度廓线的唯一解。有必要建立约束和辅助其他观测信息进行求解。4 附加约束条件的层析算法构建与改进4.1 基于高斯平滑的水平方向上的约束方程

大气水汽在水平空间分布上存在连续性,对于静稳天气可以根据距离越近大气特性相关性越强的原则,建立起某网格待求参量xρsNwi(后面简称xi)与其周围等高网格中参量的数学关系。该约束方程的基本形式可以表示为

式中,xi为计算网格i的未知待求参量;wk,i为约束权重系数,k≠i,n为总的网格数,k=1,……,n-1。设定参与加权的网格位于不同高度层时,wk值为0;同一高度的网格利用高斯加权函数的方法来确定权重系数的值,权重与网格间的距离大小有关。可表示为

式中,Hi为计算网格的海拔高度(单位:m),Hk为参与加权的网格高度(单位:m)。dik为所有网格到计算网格的距离(单位:km),σ 为平滑因子(单位:km)。增加(0,1)函数M(Hi,Hj),以避免其他高度网格参量参与加权,扰乱水平约束。

采用高斯加权的水平约束矩阵方程可表达为

HX=0

包含水平约束项的线性系统方程为

4.2 基于探空观测的垂直方向的约束方程

根据武汉、宜昌等探空2006—2009年的资料,计算该地的各高度层的湿折射度式(6)、或水汽密度指数式(18),并进行统计分析,得到其均值和标准差。各高度层的湿折射度(或水汽密度指数)的多年平均值作为水汽垂直分布结构的先验信息,各个高度上的标准差可建立起该约束方程的权重矩阵,并建立约束方程。基于探空信息建立的垂直约束方程的数学模型可表示为

式中,k为垂直高度划分的层数,L1L2,…,Lk均为28维行向量(28为同一层内水平网格的数目),Lk=[xk0,xk0,…,xk0]。将探空资料转换的湿折射度(或水汽密度指数)统计结果为每个网格赋值,B是单位矩阵,xk0为对探空资料第k层高度处的水汽参量进行统计得到的均值。采用2006—2009年武汉探空站每年12个月的所有的探空为总样本,得到垂直约束需要的均值和标准差(图 4a),均值为先验值,标准差为权重矩阵。也分别采用每月4 a的探空为样本,得到该月份期间数据层析中垂直约束需要的均值和标准差,图 4b给出了7月的两个垂直约束条件。图 4a和b比较,对于全年和单月统计出来的先验值差异不大,但权重矩阵有较大差异,采用单月的统计结果建立垂直约束。
图 4 由2006—2009年武汉探空统计的湿折射度廓线的均值(mean)和标准差(std)(a.年统计,b. 7月统计) Fig. 4 Statistical average value(mean) and st and ard deviation(std)of the wet refraction by the radiosonde data at Wuhan site during 2006-2009(a. annual statistics results,b. monthly statistics results in July)

水汽密度(ρw)公式用于计算探空的水汽密度先验值。由水汽状态方程:e=ρwRvT得出

式中,e为水汽压(单位:hPa),可由露点温度定义建立与水汽和饱和水汽压的关系,结合饱和水汽压Tetens经验公式(盛裴轩等,2003)得到水汽压
式中,Td为露点温度(单位:K),es为饱和水汽压,a0b0的取值满足如下条件:

a0=7.5,b0=237.3(T≥273.15 K)

a0=9.5,b0=265.7(T<273.15 K)

结合GPS水汽层析的观测方程及水平、垂直约束方程,可得到如下方程组:
其中,HX=0表示水平约束方程,BX=V表示利用先验信息构建的垂直约束方程。根据最小二乘法,从协方差矩阵PSWVSWDPhPv组成的约束式中就可以求出层析解
式中,PSWVSWD为观测方程的权重矩阵,PhPv分别为水平约束方程和垂直约束方程的权重矩阵。4.3 地面观测的引入

基于上述水平约束和垂直约束,基本可以对三维水汽层析方程组进行求解,为充分利用观测资料,提高三维水汽层析精度,在观测方程(10)中增加一层近地面网格(图 5),即将低层1000 m高的网格分为两层,0—500 m的一层网格中心的水汽密度由地面观测获取。利用水汽密度随高度负指数变化函数(Tomasi,1981),由GPS接收站同址的地面观测推导出,表示为

图 5 地面观测引入层析算法的示意图 Fig. 5 LocationofthenetwSchematic diagram for adding the ground surface observations in the tomography algorithmorking
式中,ρwzρw分别表示高度Z和地面的水汽密度,水汽标高为Hw(单位:mm),可表示为
引入地面观测后建立新的方程
式中,sm0为第m号斜路径射线近地面的250 m高度以下的射线长度,xρsNwm0m号射线在250 m高度的湿折射度或水汽密度。设

加入地面观测后,约束层析观测方程为:

根据最小二乘法,从协方差矩阵PSWVSWD′PhPv组成的约束式中就可以求出层析解

5 三维水汽层析结果检验5.1 三维层析产品

采用2010年6、7月鄂东区域加密GPS网的观测开展了三维水汽层析研究,试验利用美国麻省理工学院(MIT)提供的GAMIT软件,通过联测国际GPS服务站解算出各站天顶湿延迟,天顶延迟分辨率30 min,水平梯度分辨率1 h,映射函数采用动态映射函数;然后采用顾及双差残差算法来进一步获取穿过层析网格模型顶层的GPS斜路径湿延迟观测量。基于获得的湿延迟观测量,层析窗口取1 h,采用约束的层析算法得到逐时的区域湿折射度和水汽密度的三维分布。图 6、7分别为基于2010年7月10日02时(北京时)22个GPS站观测资料层析的三维湿折射度和三维水汽密度。比较可见同一高度上反演的湿折射度(图 6)和水汽密度(图 7)分布较为一致。

图 6 2010年7月10日02时鄂东区域上空不同高度层的湿折射度水平分布(a. 500 m,b. 2500 m,c. 5500 m,d. 9500 m) Fig. 6 Horizontal distributions of the wet refraction over the Eastern Hubei at 02:00 BT 10 July 2010 at the various height layers of 500 m(a),2500 m(b),5500 m(c) and 9500 m(d)
图 7 2010年7月10日02时鄂东区域上空不同高度层的水汽密度水平分布(a.500 m,b.2500 m,c.5500 m,d.9500 m) Fig. 7 Horizontal distributions of the water vapor density over the Eastern Hubei at 02:00 BT 10 July 2010 at the various height layers of 500 m(a),2500 m(b),5500 m(c) and 9500 m(d)
5.2 总体反演精度

约束层析算法中,可以由SWDSWV分别推导出湿折射度和水汽密度的三维分布。区域内以武汉(btcd)、咸宁站(bfxp)的探空观测作为参照,与层析结果进行比对检验。根据式(6)和式(19)可以由探空的露点温度和温度推导湿折射度,根据式(18)—(19)可由探空的温度、露点得到水汽密度。

选取咸宁2010年7月9日20时、10日02时两个时次的探空为参照,与层析结果比较。咸宁GPS移动探空每2 s返回一组数据,垂直分辨率最高可达10 m量级。从图 8可见,两种探测获得的湿折射度廓线变化趋势一致,三维水汽层析得到的廓线值偏小。02时两种探测的廓线更接近。从图 9可见,两种探测获得的水汽密度线变化趋势一致,三维水汽层析得到的廓线值偏小。比较图 8和9,水汽密度廓线与探空的相关性要强。

图 8 2010年7月9日20时(a)、10日02时(b)咸宁站层析(GPS)与探空观测(RS)的湿折射度对比 Fig. 8 Wet refraction profile comparison between tomography(GPS)results and radiosonde observations(RS)over Xianning site at 20:00 BT 9(a) and 02:00 BT 10(b)July 2010
图 9 2010年7月9日20时(a)、10日02时(b)咸宁站层析(GPS)与探空观测(RS)的水汽密度对比 Fig. 9 Water vapor density profile comparison between tomography(GPS)results and radiosonde observations(RS)over Xianning site at 20:00 BT 9 and 02:00 BT 10 July 2010

整理咸宁和武汉两站的探空资料得到1161个样本用于检验层析结果,代表层析的总体精度。三维层析反演湿折射度廓线、水汽密度与探空比较的平均偏差、标准差和相关性,由图 10可见:总体上层析结果小于探空观测,平均偏差都为负值(表 1);水汽密度层析结果与观测的相关优于湿折射度,湿折射度廓线与探空相关系数为0.9783,水汽密度与探空相关系数为0.9830(表 1);加低层观测的层析结果显示,层析与探空的相关提高了1%,水汽密度整体均差不变,但标准差增大。层析结果整体较探空观测偏小,可能与模型网格较大以及假设网格内水汽均匀稳定有关,此外,在由地面观测推导低层水汽算法中,假设近地面高度以下水汽密度随高度负指数变化,虽然是基于统计规律,而实际情况要复杂得多,低层的水汽偏差也会影响反演结果。低层观测的加入增强区域水汽的非线性特征,若方程组使用线性求解算法,引起解的离散度增大。

图 10 1161个样本三维层析结果与探空观测的平均偏差和相关性(a、b、c.层析湿折射度和加入地面观测前、后层析的水汽密度与探空的平均偏差;d、e、f.层析湿折射度和加入地面观测前、后层析的水汽密度与探空的相关散点图) Fig. 10 10 Average deviation and correlation of the 1161 samples between 3D tomographic results and radiosonde observations,a.-c.(d.-f.)are the average deviation(correlation)between the radiosonde observations and respectively the wet fraction,the water vapor densities before and after adding the ground surface observations to the tomographic algorithm
表 1 层析结果与探空观测对比 Table 1 Comparison between tomographic results and radiosonde observations
样本 湿折射度 水汽密度 加低层观测后的水汽密度
均差
(mm/km)
标准差
(mm/km)
相关系数均差
(g/m3)
标准差
(g/m3)
相关系数 均差
(g/m3)
标准差
(g/m3)
相关系数
所有探空 -1.74 8.48 0.978-0.63 1.22 0.983 -0.63 1.32 0.985
2 km以下 -1.81 1.47 0.865-0.20 1.89 0.742 -1.03 1.86 0.777
武汉站 1.18 8.43 0.978 -0.16 1.18 0.983 -0.14 1.25 0.985
咸宁站 -4.50 7.57 0.983 -1.06 1.10 0.988 -1.09 1.21 0.990
5.3 低层反演精度

表 1的数据可见,2 km以下,层析结果还是小于探空,平均偏差为负值,2 km以下水汽密度的均差占所有样本的31%;加入低层观测后总体均差不变,2 km以下水汽密度均差增大,为总体均差的1.68倍,说明低层观测的加入,尽管使低层层析结果偏小,但高层的层析结果偏大,更接近于真实值;2 km以下湿折射度与探空相关系数为0.87,水汽密度与探空相关系数为0.74,与整体样本显示的情况相反,说明2 km以下的湿折射度和2 km以上的水汽密度更能接近大气水汽的特征;加入地面观测约束后,2 km以下水汽密度与探空的相关提高了4.3%,并且,标准偏差减小2%,对2 km以下水汽密度的精度有明显提高;2 km以下的水汽密度标准差比总体样本的要大,说明观测方程中对于低层的未知量的求解还存在不稳定,加地面观测后2 km以下标准差减小,总体标准差增大,说明加地面观测对高层的求解有影响,近地层水汽的垂直变化较为复杂,在引入低层观测时假设近地面高度以下水汽密度随高度负指数变化,基于统计规律的假设与实际状况的差异对反演结果多少会有影响。5.4 层析区域边缘站与中心站水汽反演精度

根据层析区域GPS站点分布(图 3a),武汉站位于区域中心网格,咸宁站在网格靠近区域边缘位置,选取武汉和咸宁分别作为GPS层析区域的中心和边缘代表站,分析区域中心和边缘的层析精度。武汉站探空有L波段的无线电探空仪,咸宁站有GPS移动探空,对比检验资料来源于2010年汛期探空加密观测试验。

(1)层析中心武汉站上空的层析精度

武汉站累计2010年汛期所有可以对比的数据作为统计样本共571个。如图 11显示层析中心湿折射度与探空的平均偏差为正值,整体层析结果大于探空观测,可能与处于层析中心观测元较为丰富有关;湿折射度、水汽密度与加地面观测后层析的水汽密度和探空的相关逐渐增强,并都接近各自同类整体样本的相关系数,说明区域中心站在相关性方面有代表性,加低层观测后提高了中心层析结果的相关性;加低层观测后的层析中心水汽密度与探空的平均偏差减小了1.5%(表 1),可见精度有所提高;加了低层观测后,武汉站层析与探空的标准差增加了6%,整体的标准差增加了8%,中心层析结果受低层观测加入后的扰动小于整体水平。

图 11 武汉站571样本三维层析结果与探空观测的平均偏差和相关性(a、b、c.层析湿折射度和加入地面观测前、后层析的水汽密度与探空的平均偏差;d、e、f.层析湿折射度和加入地面观测前、后层析的水汽密度与探空的相关散点图) Fig. 11 Average deviation and correlation of the 571 samples over Wuhan site between 3D tomographic results and radiosonde observations,a-c(d-f)are the average deviation(correlation)between the radiosonde observations and respectively the wet fraction,the water vapor densities before and after adding the ground surface observations to the tomographic algorithm

(2)层析边缘咸宁站上空水汽层析精度

位于层析区域边缘的咸宁站累计样本共590个。如图 12显示,湿折射度、水汽密度与加地面观测后的水汽密度和探空的相关系数依次提高,显示低层观测的加入提高了层析与探空的相关。而边缘各类层析结果与探空平均偏差为负值,整体显示比探空观测要小,比层析中心各类层析结果与探空的均差要高出一个量级,是整体层析均差的两倍,因为层析边缘观测元数量不够多,层析精度有待提高;由图 12表 1也可知,加低层观测后,扰动对咸宁站层析结果影响较武汉站明显,标准差增大了10%。

图 12 咸宁站590样本三维层析结果与探空观测的平均偏差和相关性(a、b、c.层析湿折射度、加入地面观测前、后层析的水汽密度与探空的平均偏差;d、e、f.层析湿折射度、加入地面观测前、后层析的水汽密度与探空的相关散点图) Fig. 12 As in Fig. 11 but for the 590 samples over Xianning site

无论对于层析中心还是边缘,层析水汽密度与探空的相关大于湿折射度,加地面观测后的层析能进一步提高与探空的相关。低层观测的加入,明显减小中心层析结果的均差。由于所处层析区域不同,中心观测元数量较多,使得中心层析的精度好于整体和边缘层析的精度。低层观测的加入对整体标准差有影响,尤其对边缘层析结果影响明显,使得层析标准差增大,层析解出现不稳定。

在层析技术数值计算中由于穿过网格长度长短不一致,长的达3.5×104 m,短的可以到10 m,使得观测方程(24)中穿过网格长度的设计矩阵S各元间存在较大的量级差异,将会在矩阵求解过程中造成大的扰动,引起大的舍入误差,导致方程解的不稳定。而低层观测的加入,增加了0—500 m一层的观测值,将原有的第1层网格平分为两层,虽然穿越网格的射线增多,但网格内出现较多的短的线段,由此更加剧了观测方程的病态特征和解的不稳定。如何求解该类病态方程是层析算法中面临的另一个问题,正在进一步解决中。6 结论和讨论

为提高地基GPS三维大气水汽反演精度,借助密集的地基GPS斜路径水汽观测、多年探空观测和地面观测,开展了基于约束的层析算法研究,得到以下结论:

(1) 以4 a的探空观测为垂直约束,即逐月的探空平均为先验值,标准偏差为权重矩阵,以高斯函数为水平约束,构建了基于探空观测的约束层析算法,基于该算法不仅提供大气三维湿折射度,并能提供三维的水汽密度。为改善层析精度,算法中引入了地面观测信息。

(2) 建立了湖北省鄂东区域GPS加密大气水汽层析网格模型,基于上述观测约束层析解算技术实现了准实时三维水汽层析自动解算系统,实现逐时的层析结果输出。

(3) 基于探空观测约束条件的层析得到水汽密度与探空相关性好于湿折射度,整体样本检验三维水汽密度平均偏差为-0.63 g/m3,标准偏差为1.22 g/m3,与探空相关系数为0.98。

(4) 由于所处层析区域不同,中心观测元数量较多,使得中心层析的精度好于整体和边缘层析的精度。

(5) 加入低层观测的层析结果与探空的相关性好于未加低层观测时的状况。低层观测的加入提高了层析与探空的相关,减小了中心层析均差、低层层析标准差和2 km以上的层析均差,对提高层析精度有改进作用

低层观测的加入在对层析相关性和均差方面有改进的同时,对整体标准差有负影响,尤其对边缘层析标准差负影响明显。初步分析与加剧观测方程中长度矩阵元素间的量级差异有关,由于元素量级的较大差异将会在矩阵求解过程中造成大的扰动,引起大的舍入误差,导致方程解的不稳定性问题。

三维水汽层析解算是寻求唯一解和稳定解的过程,通过观测约束已初步获得了三维水汽的唯一解;为进一步提高层析精度,需要探索一种合适的算法解决解的稳定性问题。

致谢: 感谢毛节泰教授、曹云昌研究员在层析算法上给予的有益讨论和宝贵的建议;感谢麻省理工学院地球大气与行星科学系无偿提供的GAMIT软件。

参考文献
曹云昌, 陈永奇, 李炳华等. 2006. 利用地基GPS测量大气水汽廓线的方法. 气象科技, 34(3): 241-245
毛节泰. 1993. GPS的气象应用. 气象科技, (4): 45-49
盛裴轩, 毛节泰, 李建国等. 2003. 大气物理学. 北京: 北京大学出版社,522pp
宋淑丽, 朱文耀, 丁金才等. 2005. 上海GPS网层析水汽三维分布改善数值预报湿度场. 科学通报, 50(20): 2271-2277
徐桂荣, 万蓉, 李武阶等. 2009. 地基GPS反演大气可降水量方法的改进. 暴雨灾害, 28(3): 203-209
袁招洪. 2005. GPS可降水量资料应用于MM5模式的变分同化试验. 气象学报, 63(4): 391-404
张双成, 叶世榕, 万蓉等. 2008. 基于Kalman滤波的断层扫描初步层析水汽湿折射率分布. 武汉大学学报(信息科学版), 33(8): 797-799, 809
Bi Y M, Mao J T, Li C C. 2006. Preliminary results of 4-D water vapor tomography in the troposphere using GPS. Adv Atmos Sci, 23(4): 551-560
Boehm J, Niell A, Tregoning P, et al. 2006. Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data. Geophys Res Lett, 33(L07304), doi: 10.1029/2005GL025546
Braun J, Rocken C, Liljegren J. 2003. Comparisons of line-of-sight water vapor observations using the Global Positioning System and a pointing microwave radiometer. J Atmos Oceanic Technol, 20: 606-612
de Haan S, Barlag S, Klein-Baltink H, et al. 2003. Synergetic use of GPS water vapor and meteosat images for synoptic weather forecasting. J Appl Meteor, 43(3): 514-518
Duan J P, Bevis M, Fang P, et al. 1996. GPS meteorology: Direct estimation of the absolute value of precipitable water. J Appl Meteor, 35(6): 830-838
Flores A, Ruffini G, Rius A. 2000. 4D Tropospheric tomography using GPS slant wet delays. Ann Geophys, 18(2): 223-234
GAMIT Release10.4. 2010. http://www-gpsg.mit.edu/-simon/gtgk/index.htm
Gradinarsky L P, Jarlemark P. 2002. GPS tomography using the permanent network in Goteborg: simulations//Proceedings of IEEE Positioning Location and Navigation Symposium. Palms Springs, CA: IEEE, 128-133
Gradinarsky L P. 2002. Sensing Atmospheric Water Vapor Using Radio Waves. Department of Radio and Space Science, School of Electrical Engineering, Chalmers University of Technology
MacDonald A E, Xie Y F, Ware R H. 2002. Diagnosis of three-dimensional water vapor using a GPS network. Mon Wea Rev, 130: 386-397
Rüeger J M. 2002. Refractive index formulae for radio waves. Washington, DC: FIG XXII International Congress, 13pp
Skone S, Hoyle V. 2004. Troposphere Modeling in a Regional GPS Network. Sydney, Australia: International Symposium on GNSS/GPS, 15pp
Tomasi C. 1981. Determination of the total precipitable water by varying the intercept in Reitian’s relationship. J Appl Meteor, 20: 1058-1069