2010, Vol. 53
2. 中国科学院上海天文台, 上海 200030
2. Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China
从20世纪80年代后期开始，随着GPS观测精度的提高，和通讯、计算机和卫星技术的发展，兴起了利用GPS观测地球大气的研究，即通过测量穿过大气层的GPS信号延迟，反演大气中的水汽、温度、压强等信息，这种方法称为GPS气象学[1].根据观测模式的差别，它可分为地基和空基GPS气象学.地基GPS气象学利用固定在地球表面上的GPS观测网络接收来自GPS卫星星座的信号，连续地对地球可降水汽进行监测，获得测站上空高时间分辨率的可降水汽垂直积分序列.与常规探空气球、水汽辐射仪和遥感卫星观测相比较，地基GPS观测可以获得高时间分辨率、高精度的可降水汽序列.
无线电信号在地球大气的传播过程中，产生了信号的延迟和弯曲.当地基GPS在低高度角观测时，大气延迟和弯曲会变得显著，从GPS观测资料可以获得高精度的大气延迟和弯曲，如3°高度角时，大气延迟和弯曲的精度分别在厘米和10-5rad量级，其中主要的贡献来自于低层大气的影响[2].当存在较大大气折射率梯度时（通常在大气边界层顶附近），低高度角的大气延迟和弯曲可能会出现急剧的增大.大气层内的大气折射率垂直梯度小于-157（N-unit）km-1时，该大气层称为陷获层（也称波导层）.在大气陷获层内，无线电信号小于临界发射角的初始发射角（绝对值）传播时，在陷获层的顶层和底层之间发生全反射，这称为大气波导.这一物理现象提供了利用高精度的地基GPS低高度角观测资料，探测低层大气特性（如大气波导、水汽分布等）的可能性.Sokolovskiy等提出利用低高度角的地基GPS资料计算大气弯曲角的方法，探索了GPS观测与气象分析场得到的弯曲角之间的差异[2].Lowey等提出GPS斜路径大气延迟反演大气折射率的垂直廓线的方法.他们认为，仅用地基GPS斜路径大气延迟反演大气折射率廓线的效果不如弯曲场的反演[3]效果好.
GPS观测得到的大气延迟和弯曲是大气折射同一物理过程的两个不同的观测效应[4].在低高度角观测中包含了更多的低层大气信息.本文发展了一种利用地基GPS大气延迟和弯曲角观测联合反演大气折射率廓线的算法；考虑了弯曲角和延迟空间之间观测协方差矩阵的转换关系；通过较为合理的三折线大气廓线模型，构造出球对称大气假设下的射线追踪算子，同时得到GPS观测大气延迟和弯曲；利用St.Helena台站的探空资料作为模拟观测，采用变分同化技术，重构大气折射率廓线.
本文的结构如下：第2节介绍球对称大气假设下的射线追踪正演算子，比较了用高分辨率的探空资料和大气NCEP（National Center for Environmental Prediction）分析场得到大气延迟和弯曲之间的差异；在第3节引入了大气三折线模型；第4节描述变分同化反演算法；第5节利用探空资料进行模拟仿真实验，并比对、验证了反演结果；第6节包含了本文的一个简单的小结和讨论.
2 射线追踪技术在地球大气中，由于气体分子、水汽、自由电子等介质的存在，电磁波在大气中的传播速度会发生变化，从而产生了大气延迟.介质的梯度使电波发生了弯曲.
在无线电波段上，大气折射率（refractivity）N与大气的密度、温度、压力、水汽等物理参量满足[5]：
|
(1) |
其中，n是大气折射指数（refractiveindex），Pd和Pw分别是大气的干分压（hPa）和湿分压（hPa），T是绝对温度（K）.
在几何光学近似和大气球对称假设下，Snell公式成立.如图 1所示，用路径积分方法可以得到GPS信号路径上的大气延迟和弯曲[3, 6].
|
图 1 GPS信号路径示意图 Fig. 1 Geometry of GPS signal |
图 1中，R和T分别对应于接收机和发射机的位置.在地球大气局部球对称假设下，同时假设原点0为地球参考椭球在R点的局部曲率中心.信号路径的相位延迟和弯曲可以分别表示为
|
(2) |
|
(3) |
式中r是信号到原点的距离，
|
(4) |
式中r′=dr/dθ，利用Snell定律
|
(5) |
式中ϕ是信号路径的天顶角，a是信号的碰撞参数.相位和弯曲可以重新改写为
|
(6) |
|
(7) |
式中x=rn（r）是折射率半径，m（x）=ln[n（x）].对于GPS卫星，点2在大气层以外，n（r2）=1.给定α，φ1，r1和r2，点1和点2之间的中心夹角θ等于
|
(8) |
点1和点2之间的距离S0为：
|
(9) |
大气延迟为ΔS=S-S0.
图 2显示用射线追踪技术得到大气延迟和弯曲序列的两个例子，其中实线和虚线分别表示利用高分辨率的探空资料（St.Helena台站坐标[15.9°N，56.6°W]）与相应的NCEP气象分析场资料得到的结果.图 2的a1和a2分别表示个例1和个例2的高度角与弯曲之间的关系，图 2的b1和b2表示个例1和个例2的高度角与大气延迟之间的关系.当高度角大于2°时，由探空（实线）和NCEP分析场（虚线）得到的大气延迟和弯曲，两者之间的差异很小，但当高度角小于2°时，探空得到的大气延迟和弯曲随高度角的减小急剧增大（弯曲角更明显），而NCEP分析场得到的弯曲和大气延迟增加的速度明显小，并且高度角越低差异越明显.造成这种现象是由于St.Helena台站位于热带海洋上，存在大气波导现象，易被高分辨率的探空资料较准确反映出来，而NCEP分析场垂直分辨率低，一般不能准确地反映出低层大气的变化情况.在低高度角时，由探空资料与气象分析场资料根据射线追踪法分别得到的弯曲和大气延迟的比对，可以看出较明显的差异.因此，低高度角的地基GPS观测可以获得低层大气中的信息.
|
图 2 探空资料与NECP资料得到的弯曲和延迟与高度角关系图 Fig. 2 Relation between elevation and bending angle (up) or phase (down) from radiosonde and NECP respectively |
在大气对流层中，可以用简单的四参数三折线模型对低层大气折射率进行模型参数化，如图 3所示，其中横坐标为折射率，纵坐标为高度.折射率可以表示成分段函数：
|
(10) |
|
图 3 四参数形式的大气折射率模型 Fig. 3 The refractivity model with four parameters |
其中，Z1和N1为GPS接收机位置的高度和大气折射率，它们从台站的气象资料得到；Z2、N2和Z3、N3为待估参数，分别为点2和点3处的高度和大气折射率；6km以上的大气折射率可以从气候模型CIRA+Q[7]中获取，N4是Z4≈6km处的折射率.
4 同化方法优化估计三折线模型参数设解向量x=（N2，N3，Z2，Z3）T.价值函数f（x）可以写成
|
(11) |
其中，ys是GPS大气延迟观测矢量，yα是GPS弯曲观测矢量，Hs{x}为大气延迟射线追踪算子，它是从解矢量x到大气延迟的映射，Hα{x}为弯曲射线追踪算子，它是从解矢量x到弯曲的映射，Cs=IE为大气延迟观测误差的协方差矩阵，I为单位矩阵，E为常数，Cα为弯曲观测误差的协方差矩阵.它们可以从弯曲角与大气延迟的关系中得到[8, 9].
采用统计搜索法，求价值函数方程（11）的极小值.高度的搜索步长为25m，大气折射率的步长为5N-unit.为了使大气折射率三折线模型更符合低层大气的特性，采用了如下约束：线段a和线段c的斜率皆小于中间线段b的斜率；Z3 -Z2小于800m.
5 模拟仿真结果采用St.Helena台站4次探空资料：07010211、07012511、07012611、07012911作为真值，利用第2节中的射线追踪技术得到GPS大气延迟和弯曲的模拟观测序列.将不同高度角序列的大气延迟和弯曲，采用上节的优化估计的算法，反演出大气折射率廓线.
图 4是观测资料在高度角为1°时的反演结果.图中实线表示大气折射率的反演结果，虚线表示探空资料的大气折射率.低高度角的地基GPS观测模拟仿真的反演结果能较好地反映低层真实大气折射率的廓线，同时也可以探测到低空的大气波导现象.根据大气波导发生条件，即大气折射率垂直梯度小于-157（N-unit）km-1，参照图 3，对图 4中CaseA、B、C、D进行处理，分别计算在点3和点2之间的大气折射率垂直梯度，结果依次为-91.4285、-980.0、-1080.0、-70.0（N-unit）km-1.可以判断只有在CaseB和CaseC中出现明显的大气波导现象.
|
图 4 反演结果示意图 Fig. 4 Results of inversion |
为了检验反演结果，分别计算观测资料在高度角为1°时的四个个例的反演结果与真值的残差
|
(12) |
式中的Ni rad为探空资料中点i对应的折射率，Ni mod为大气三折线模型反演结果对应的折射率，n为比较的个数.结果分别为8.0622、8.5124、4.3011、7.1063（相当于2%~3%折射率相对误差）.同时我们也给出了四个个例在不同高度角下面折射率残差的变化情况，如表 1所示.
|
表 1 不同高度角下的RMS Table 1 The RMS under different elevations |
从表 1中的数据可以看出，当高度角为1°时，采用的4个个例的反演结果的比较残差RMS值最小.随着高度角的变大，反演的精度也随之下降.当高度角大于3°的时候，反演结果与真值差异很大.因此对于地基GPS反演大气折射率廓线，应选取高度角低于3°的观测资料.
6 结语本文将地基观测得到的弯曲和延迟作为观测量，发展了地基GPS观测联合反演大气折射率的算法，首次比较合理地解决了价值函数中弯曲角和延迟的转换关系，将不同观测空间更加合理地结合起来，并进行模拟仿真反演的试验：采用高分辨率的探空资料为真值，利用射线追踪技术，获得仿真的GPS大气延迟和弯曲的观测序列；将低层大气的折射率廓线用三折线模型参数化，可以获得大气延迟和弯曲的序列；利用变分同化技术，优化估计三折线模型参数，反演出低层大气折射率廓线.
GPS观测技术具有高精度、高分辨率、低成本、全天候、灵活运作等特点，在处理地基GPS实测资料中，可以用本文的方法得到较低高度角的大气斜路径延迟和弯曲序列，并用优化估算方法反演出大气折射率垂直廓线.三折线模型具有它的简单性和一定程度上的合理性，但是仍然存在一定的局限性；采用单高度角获得模型参数的解稳定性是否很好有待进一步研究.改进低层大气廓线模型和本文方法结果的稳定性是将来的主要研究方向之一.在研究中所采用的搜索方法为传统的穷举搜索法，准确率虽高但搜索量大，且花费的时间较长，不宜采用进行大规模的计算.需进一步考虑采用更好的算法，如蚁群算法进行搜索.本文提出的方法可应用于大气边界层的探测，也为低层大气的结构特性（例如波导现象）的研究提供新的思路[10].
致谢衷心感谢匿名评审专家提出的宝贵修改意见.
| [1] | Yuan L L, Anthes R A, Ware R H, et al. Sensing climate change using the Global Positioning System. J.Geophys.Res. , 1993, 98(D8): 14925-14937. DOI:10.1029/93JD00948 |
| [2] | Sokolovskiy S V, Rocken C, Lowry A R. Use of GPS for estimation of bending angles of radio waves at low elevations. Radio Sci. , 2001, 36(3): 473-482. DOI:10.1029/2000RS002541 |
| [3] | Lowry A R, Rocken C, Sokolovskiy S V, et al. Vertical profiling of atmospheric refractivity from ground-based GPS. Radio Sci. , 2002, 37(3). |
| [4] | 严豪健, 符养, 洪振杰. 现代大气折射引论. 上海: 上海科技教育出版社, 2006 . Yan H J, Fu Y, Hong Z J. Introduction of the Modern Atmospheric Refraction (in Chinese). Shanghai: Shanghai Science and Technology Education Press, 2006 . |
| [5] | Bean B R, Dutton E J. Radio Meteorology. New York: Dover Publications, 1968 . |
| [6] | 洪振杰, 叶帆, 郭鹏. 非对称大气模式下二维射线追踪算子. 自然科学进展 , 2007, 17(3): 405–409. Hong Z J, Ye F, Guo P. The Ray-tracing mapping operator in an asymmetric atmosphere. Progress in Natural Science (in Chinese) , 2007, 17(3): 405-409. |
| [7] | Kirchengast G, Hafner J, Poetzi W.The CIRA86aQ_UoG model: An extension of the CIRA-86 monthly tables including humidity tables and a Fortran 95 global moist air climatology model.IMG/UoG technical report for ESA/ESTEC, number 8, Inst.für Meteorol.und Geophys., Univ.of Graz, Graz, Austria, 1999 |
| [8] | Yan H J, Ping J S. The generator function method of the tropospheric refraction corrections. A J , 1995, 110(2): 934-939. |
| [9] | Yan H J. A new expression for astronomical refraction. A J , 1996, 112(3): 1312-1316. |
| [10] | Sokolovskiy S V, Kuo Y H, Rocken C, et al. Monitoring the atmospheric boundary layer by GPS radio occultation signals recorded in open-loop mode. Geophys.Res.Lett. , 2006, 33: L12813. DOI:10.1029/2006GL025955 |