水汽是地球大气中非常重要的一种气体成分, 在天气演变与气候变化中扮演一个极为重要的角色。在通常的温度范围内, 水汽存在相态变化; 小的时空尺度范围内数量变化也很大, 因而难于捕捉水汽的详细变化信息。经过诸多大地测量学家和大气科学家的共同研究而发展起来的地基GPS探测大气水汽主要包括3个方面的内容:一是观测天顶方向的大气水汽总量PW (precipitable water vapor) ^[1-2]; 二是遥测倾斜的GPS信号路径上的水汽总量SW (slant-path water vapor) ^[3-6]; 三是应用组网的GPS倾斜路径观测反演局地上空的水汽四维时空变化, 即水汽层析 (water vapor tomography) ^[7-8]。目前, 地基GPS观测PW技术方法都很成熟, 精度可达毫米量级, 并且已转向业务应用。尽管PW精度很高, 但由于PW仅是测站上空大气水汽的一维分布, 在大气科学等领域中的应用受到了一定的限制。最近, GPS遥测倾斜路径上水汽总量SW的研究也取得了很大的进展, 研究表明, GPS观测SW的精度在毫米量级^{[6, 9-10]}, SW的可能应用包括探测锋面过境时的水汽变化、深对流的发展、低空急流水汽输送等方面。相比PW而言, SW的观测拓展到了测站上空不同的空间方位, 但它仍然不能给出水汽垂直分布信息。应用区域的GPS观测网进行水汽层析解决了以上PW, SW两种观测的不足, 可获得局地上空的水汽垂直分布状况, 时间分辨率可达30 min或更低。GPS卫星星座包括24颗GPS卫星, 一般一个测站上空可同时观测6~8颗卫星, 如果地面上小范围内有多个测站 (大于3个), 那么在一段时间之内, 多个站点对多个卫星的倾斜路径观测, 将产生相当多的SW观测, 形成对观测网上空低层大气的稠密采样, 每一条信号路径上的SW观测包含了这条路径上的水汽信息, 采用反转方法, 有可能把水汽的三维结构信息解算出来。

水汽层析被应用于数值预报模式。大气中水汽的三维空间分布信息对于中小尺度数值预报模式的研究非常重要, 归根到底, 数值预报模式是一个初/始值问题^[11], 即给定当前大气状态的初估场和合适的边界条件, 模式将能模拟出大气的演变结果。显而易见, 初值的估计越精确, 预报质量就越好, 湿度场正是模式需要的初值之一。目前数值预报模式初始场中所用的水汽分布信息主要依靠探空站提供的水汽廓线, 但是探空站站间距过大 (200 km以上), 而且每天只有两次探测, 这样的时空密度难以捕捉中小尺度水汽变化信息。依靠卫星提供的水汽资料, 精度也是一个主要问题。应用组网的地基GPS获取水汽廓线有望在中小尺度水汽探测中发挥出作用, 为数值天气预报提供初始湿度场。

层析技术主要是从沿不同路径的积分观测数据中提供模型参数的估计值, 这种方法已成功应用到医学、地震学和海洋学的研究当中。1998年, Ruffini等^[12]应用GPS进行地球电离层层析研究, 给出了全球电离层四维图像。尽管GPS应用在电离层四维层析中已比较成功, 但应用到对流层层析时在方法及数据处理等方面有诸多不同, 主要原因有3个方面, 一是对流层对信号的影响量级较小, 处理精度要求较高; 二是信号在对流层中传播时没有频散的特点, 所以需要研究计算倾斜路径延迟的方法; 三是对流层时空变化比电离层更快, 需要尺度较小的观测网进行观测。在过去的一些对流层层析研究中, 反演关注的物理量为对流层湿折射率, 基于水汽在大气科学中的重要性, 本文中反演关注的量为每个网格中的水汽密度 (绝对湿度), 以产生水汽在水平和垂直方向上的分布图像。

将观测网上空的对流层大气在水平和垂直方向上划分为许多网格或盒子, 每个网格中的水汽密度是未知的, 但在一定时间内 (如1 h) 假定每个盒子中的水汽分布均匀且密度为常数, 则每一条倾斜观测路径上的水汽总量为这条射线所穿越网格的水汽量之和:

(1)

式 (1) 中, W_S^P代表第P个射线路径上的水汽总量, A _ijk^P为射线所穿越的第 (i, j, k) 个网格的距离, X _ijk^P为水汽密度。这样, 在一段时间内, 将形成大量的上述观测方程, 求解这些方程有可能解得网格中所含的水汽量。这个问题本质上是一个反问题, 即已知观测数据, 如何由这些数据来求得模型参数的估计值。

仔细分析发现, 即使方程数目大于未知量数目, 此方程组也难以解算。这是因为GPS卫星在天空中的分布并不均匀, 而且地面GPS观测网的几何构建并不能保证地面上每个网格点都有测站, 致使某些网格没有观测射线穿过。观测方程组表面上超定, 实质上却是欠定的, 而且许多单个观测方程的系数为0, 因此观测方程组属于混定方程组。为解决没有射线穿过盒子导致水汽密度不能确定的问题, 假定同一水平方向上的网格内的水汽分布相互关联, 称之为“水平平滑条件”, 即

(2)

式 (2) 中w _j为第j个网格的权重。假设距离越相近, 则相关性越强。将上式移项整理, 并用矩阵形式表达, 可写为: wX=0。

包括上述平滑条件, 从而形成如下层析观测方程组:

(3)

式 (3) 中, W_S代表倾斜路径水汽总量, A代表网格路径长度, w为权重。从上面的论述可以看到, GPS水汽层析本质上属于反问题的研究, 即通过最终的观测值, 来反求模型的参数值。以Y代表方程左侧的量, 以X代表待求模型的参数值, 则可将上述方程简写为以下线性方程形式:

(4)

很显然, 目的是反转上述线性系统并求得X值。一般来说, 观测值多于未知量个数 (模型参数), 使下式最小化, 可求得最小二乘解:

(5)

式 (5) 中, N_Y是观测量的协方差矩阵。假设协方差是对角矩阵, 观测是等权的, 且观测误差呈高斯分布, 那么, 最小二乘解等同于最大似然解, 表达式可减少为:

(6)

观测方程及未知量个数非常多, 使得采用通常方法求解上述方程中矩阵的逆非常困难, 为此, 采用奇异值分解法 (singular value decomposition, 简称SVD) ^[7]来求解矩阵的广义逆。这种方法是分析线性系统解的有力工具。SVD计算系统的特征值, 并以零特征值表征解中出现零的位置。一个矩阵A可被分解为如下形式:

(7)

式 (7) 中, U表示N×N阶正交矩阵, N为观测量个数; Λ是N×M阶对角矩阵, M是模型参数个数, V是M×M阶矩阵; 利用上述表达式, A的广义逆可表达为:

(8)

在Λ ^-1中, 与零特征值相应的地方也设为零, 而不是无穷。零特征值的出现意味着模型参数无法确定下来。如果想要得到一个良好的解, 应该尽可能减少零特征值的出现。

为了研究如何利用组网的GPS进行水汽层析反演, 获取水汽场的三维结构, 探索GPS层析的可能应用, 北京大学大气科学系、北京合众思壮科技有限公司、海南省气象局共同合作, 于2005年11月在海口周围进行了小范围的地基GPS水汽观测试验。11月北方已是冬季, 较为干旱, 而海南地处中国南端, 气候较为炎热, 水汽较我国北方充沛, 从地理条件上来看, 是有利于试验开展的。试验一共使用了4台GPS设备, 布设在海口站周围, 站点布设如图 1所示, 包括海口站、临高站、澄迈站、定安站。整个观测区域地势较为低平, 海拔约在20 m左右。海口站GPS为合众思壮科技有限公司的设备, 其余3站为Leica公司设备, 天线为扼流圈天线。从图 1可以看到, 临高、澄迈和定安分布在海口站西南部及南部, 4站间距分布较为均匀。同时, 4个站还是气象观测站, 每天进行标准的温度、湿度、气压等气象数据的观测, 海口站也是探空站点, 每天00:00和12:00(世界时, 下同) 释放探空气球获取湿度廓线, 因此海口站是重点研究的观测站点。

图 1

图 1. 海南地基GPS观测网及水平网格划分示意图 Fig 1. Location of the ground-based GPS network in Hainan with the horizontal grids

采用GAMIT软件处理GPS原始数据, 首先计算出精确的站坐标, 然后依次估计天顶延迟, 估算倾斜路径延迟及所有GPS站到所有可见GPS卫星的SW。由于试验站点间距较小, 直接处理会得到相对天顶延迟, 为了得到绝对的对流层延迟量, 需要引入较远距离的观测站点 (本次试验中引入的是北京房山、昆明等IGS的站点)。数据处理时, 为消除多路径效应等的影响, 卫星截止角为15°, 天顶延迟每30 min计算一次, 大气梯度每2 h计算一次, 映射函数采用NMF函数 (Niell映射函数)。获得层析所需的观测量后, 进行层析水汽反演。

综合考虑4站GPS数据采集情况以及匹配气象数据的观测情况, 选取2005年11月17-21日 (年积日321~325) 的观测数据进行了层析反演试验。层析水平区域网格划分如图 1所示, 经度为109.17°~110.83°E, 格距大小为0.33°, 纬度为19.33° 20.17°N, 格距大小为0.13°, 垂直方向从地面开始, 每1 km一层, 至10 km高度, 网格共5×5×10。在这次试验中, 首先用5 d时段的探空平均廓线作层析先验值, 按照第1章所述的数据处理方法, 进行了层析反演, 获取了每小时1次的水汽廓线。由于海口站同时也是探空站, 所以海口站上空GPS层析结果可与单次探空探测结果作一比对。从整个数据集来看, 平均差0.11g·m^-3, 均方根误差 (RMSE) 为0.53g·m^-3, 部分对比结果示于图 2至图 3中。在5 km以下水汽较多的区域, GPS与探空符合较好, 在观测后期大气变干过程中, GPS水汽明显偏离先验廓线, 几乎与反演时段内释放的单次探空廓线完全重合, 显示出GPS观测的作用。在5 km以上, 由于水汽量较少, 扰动因素明显, 部分反演时刻层析结果与探空相比, 出现差异较大的点。在此应该说明的是, GPS与探空二者之间的偏差并不全是GPS误差造成的, 探空水汽测量本身也有误差, 而且在匹配空间和时间上不同也会带来误差, GPS是对一段时间内大气的平均状况的反映, 探空反映的是大气瞬时状况。综合试验期间探空探测的水汽变化情况来看, 试验期间大气水汽逐渐变少, GPS与单次探空符合较好, 反映出了此次大气变干过程。

图 2

图 2. 2005年11月20日 (年积日324) 海口站层析廓线与探空对比图 Fig 2. Comparision of water vapor vertical profiles between GPS and radiosonde at Haikou station at 00:00-01:00 and 12:00-13:00 on Nov 20, 2005

图 3

图 3. 2005年11月21日 (年积日325) 海口站层析廓线与探空对比图 Fig 3. Comparision of water vapor vertical profiles between GPS and radiosonde at Haikou station at 00:00-01:00 and 12:00-13:00 on Nov 21, 2005

通过上述与探空的比较, 再一次可以看出GPS观测的明显优势, 在此GPS网观测区域内, 探空只能在00:00和12:00给出一个观测站点的探测结果, 时间和空间分辨率不足以捕捉中小尺度水汽变化信息, 而GPS时间分辨率远高于探空, 可达每小时1次, 给出了更详细的水汽变化信息。

另外, 分别测试了600 m和800 m垂直分辨率下层析情况, 结果与1 km分辨率相近, 数据对比情况示于表 1。3种分辨率下, 与探空相比平均差约为0.1g·m^-3, 均方根误差约为0.5g·m^-3。因此, 可以认为GPS水汽层析垂直分辨率可选择在500~1000 m之间。

表 1

表 1 不同分辨率下的层析结果统计分析 Table 1 Tomographic statistics results of different resolutions

表 1 不同分辨率下的层析结果统计分析 Table 1 Tomographic statistics results of different resolutions

GPS层析中水汽分布的先验值有助于解决水汽的垂直结构信息。这就存在一个问题, 不同的先验信息会对层析结果造成什么样的影响?为了回答这一问题, 应用前述的层析方法, 进行了不同的先验信息下的层析测试。从以往的层析研究来看, 大多采用1000 m的分辨率^[13-14]; 就海南地区来看, 从层析结果可以看出, 不同垂直分辨率下的层析结果相差不大, 因此选定1000 m的垂直分辨率进行测试。在这一垂直分辨率下, 再采用3种不同的先验信息方案来测试结果。第1种是所有的层析先验信息采用最早的一次探空, 即2005年11月17日00:00的探空结果; 第2种是在进行此次层析反演之前, 采用距离此次反演时间最近的前一次探空结果; 第3种是采用5 d平均的探空廓线方案。层析的时间段为11月18-21日。在第1种和第3种方案中, 先验信息一直保持不变, 而在第2种方案中, 先验信息是不断变化的。

图 4给出了11月21日采用不同先验信息的层析结果与实测探空的对比, 从图中可以看出, 在这一天的两次反演对比中, 采用11月17日00:00的探空为先验值的第1种方案得到的层析结果在低层与实测探空差异较大, 采用前一次探空为先验值的第2种层析结果稍好于第1种方案, 以平均廓线为先验值的第3种层析方案的结果与实测探空符合最好。图 5给出了5 d试验期间不同方案整个数据集上所有结果的对比散点图, 在3种层析与探空的对比中, 第1种方案的RMSE为0.71, 第2种方案为0.53, 第3种方案为0.47。从上述对比可以看出, 先验值会对层析结果造成一定的影响, 先验值较为准确的层析结果优于先验值较差的结果。第2种方案和第3种方案的结果相近, 差异不大。因此, 在实际层析应用中可采用邻近探空为先验值或最近几天的平均探空为先验值。

图 4

图 4. 2005年11月21日 (年积日325) 采用不同先验信息的层析结果与实测探空对比 (十字形:方案一, 即采用第321天00:00探空为先验值的层析结果; 圆形:方案二, 即采用时间最近的上次探空为先验值; 正方形:方案三, 即用平均探空为先验值; 菱形:第325天的反演时刻的单次探空) Fig 4. Tomographic results from different a priori information versus observations by radiosonde (line with cross-shaped symbol is the first plan using radiosonde in 321 day of year as a priori values; line with circle symbol represents the second plan using the latest radiosonde before this time as a priori information; line with square symbol is the third plan using a priori information from average radionsonde; line with diamond symbol shows the single radiosonde in day 325)

图 5

图 5. 试验期间3种先验值层析方案与探空对比的散点图 (a) 采用2005年11月17日00:00探空为先验值的层析结果, (b) 采用反演时刻前一次探空为先验值的层析结果, (c) 以平均探空廓线为先验值的结果 Fig 5. Scatterplots of results from radiosonde and tomography during the experiment (a) tomographic results using radiosonde at 00:00 as a priori information on Nov 17, 2005, (b) using latest radiosonde before the time of retrieval, (c) using average radiosonde

由上述分析可见, 给出水汽的垂直分布先验信息, 应用层析方法, 利用GPS观测量可反演出更为精确的水汽分布信息。图 6给出了2005年11月21日单个时次GPS层析、先验信息分别与探空的对比散点图。由图 6可以看出, 先验信息与所给层析时刻的探空相比, 差异较大; 经过GPS层析后, 水汽分布与探空更为接近, 表现在通过原点的拟合直线更趋近于y=x或者均方根误差得到了明显降低。00:00和12:00两次与探空的对比表明, 先验信息与探空相比均方根误差分别为1.02和0.84, 但经过层析反演后, 均方根误差分别降为0.24和0.27。

图 6

图 6. 2005年11月21日GPS层析、先验信息分别与探空的对比散点图 (菱形:先验信息与探空的对比; 实线:通过原点的拟合直线; 三角形: GPS层析与探空的对比; 虚线:通过原点的拟合直线) Fig 6. Scatterplots of GPS tomography, a priori information and radiosonde on Nov 21, 2005 (diamond symbol denotes comparison between a priori information and radiosonde, solid line is fitting line passing origin; triangle symbol denotes the results from GPS and radiosonde; dashed line is fitting line passing origin)

水汽会对GPS信号产生显著而独特的延迟, 利用这一信息可应用GPS观测网的SW观测量层析反演出对流层大气水汽的垂直廓线结构^[15]。对水汽场的认知不足, 限制了数值天气预报结果的准确性, 层析可能的应用包括为中小尺度的数值天气预报模式提供好的水汽初始场以及为中小尺度的复杂天气过程提供水汽图像, 从而提高数值天气预报尤其是临近预报的准确性。目前, 地基GPS观测在其他领域的应用也正迅速发展, GPS观测网的数目在迅速增加, 这客观上为GPS在气象学上的应用提供了便利的条件^[16]。

GPS层析结果受许多因素的限制, 还有许多问题有待研究。问题之一来源于风的影响, 本研究没有考虑这一点, 而是假设1 h内水汽的变化可以忽略, 也就是说, 背景风的影响忽略不计。而实际上大气总处在不停的水平和垂直运动当中, Noguchi等^[17]曾提出一种考虑水平风场影响的方法, 他们在层析方法中假设整个空间网格以背景风速移动, 但这也只是为解决水平风影响提供了一种可能, 实际大气状况下, 中小尺度的风场变化是很难以预知的。问题之二在于GPS观测网的几何结构对水汽分辨率的影响。若所有GPS观测站点位于同一水平面上, 这不利于获取水汽的垂直结构信息; 另一方面, GPS观测站点之间的距离也会对水汽分辨率产生重要的影响, 合理的观测网应该保证站间距均匀, 尽量使大多数网格有GPS信号射线穿过。为证明这一点, 可以假设一种极端的情况, 即所有观测站距离很近, 以至同处在地面同一个网格或盒子当中, 导致相当于用一个测站来层析水汽, 这显然不能得到水汽的垂直变化信息。稠密的GPS观测站点允许将观测网上空大气分为较小的体积单元, 本文研究只用了4个GPS观测站点, 体积单元的划分不宜太小。在垂直分辨率方面, 前人研究并未给出适宜的分层方法, 在本文的结果当中3种分辨率下的结果相差不大, 因此认为, 垂直分辨率可选择在500~1000 m之间。

因此, 在以后的研究或应用中, 应该综合考虑以下问题: ①地基GPS网的布站, 在实际应用中, 这一问题非常重要, 原因在于GPS站网的合理布局有助于获得水汽廓线信息; ②风的影响, 风的水平分量和垂直分量使得单个网格内水汽在一段时间内保持不变的假定受到影响, 如何考虑风的作用, 值得研究; ③水平约束影响, 水平约束条件的使用是为了解决某些网格无法解算的问题, 但这一假定会对水汽的水平分布产生平滑作用; ④水汽结构的检验问题, 虽然GPS水汽层析得到的小尺度水汽结构的全面检验较为困难, 但与小型飞行探测器、卫星资料或数值模式结果进行比对是可行的, 这一问题有待于继续研究。

应用海南地基GPS试验网, 成功得到了2005年11月17-21日时段内时间分辨率为1 h的水汽廓线。初步得到以下3点结论:

1) 层析以多天无线电探空探测平均结果为先验信息, 得到的结果与单次探空相比, 廓线符合一致性较好, 反映了试验期间水汽减少大气变干的过程。

2) 测试了不同先验信息方案下的层析结果, 结果表明, GPS层析对先验信息的调整明显, 较好的先验信息有助于层析反演。

3) 统计结果显示, 不同垂直分辨率下 (600 m, 800 m, 1000 m) 层析结果均方根误差在0.5g·m^-3左右。层析垂直分辨率可选择在500~1000 m之间。

层析使得地基GPS观测在大气科学等领域拥有广阔的应用前景。GPS层析反演出的水汽廓线信息在数值天气预报领域有重要的应用价值, 也可应用于确认模式输出结果, 描绘复杂天气现象中的水汽结构, 水循环研究, 合成孔径雷达成像、精密定位等方面。