COSMIC-2掩星反演数据质量分析

薛子悦，鲍艳松，唐歌实，成巍，朱孟斌，袁帅

(1.南京信息工程大学气象灾害预报预警与评估协同创新中心/中国气象局气溶胶与云降水重点开放实验室/气象环境卫星工程与应用联合实验室，江苏 南京210044；2.南京信息工程大学大气物理学院，江苏 南京210044；3.南京信息工程大学，江苏 南京210044；4.北京应用气象研究所，北京100029)

1 引 言

无线电掩星技术起源于天文学中的掩星现象，最早被用于太阳系中其他行星的大气探测。在宇宙探索初期，研究学者发现无线电波信号在经过行星大气层时会产生系统延迟和偏折，而这种系统性的偏折正是由于大气中介质的垂直折射指数变化引起的，包含着大气层的信息。而利用人造地球卫星的无线电掩星探测技术来获取地球大气信息的思想最早是在上世纪六十年代由美国喷气推进实验室JPL提出[1-2]，经过不断的技术演变，人们通过在近地轨道(Low Earth Orbit,LEO)卫星上安装一个双频接收机来接收GNSS(Global Navigation Satellite System)信号，获取由传播介质的垂直折射指数变化而带来的信号延迟，得到大气参数信息。以此为原理的掩星探测技术诞生并迅速发展起来。作为二十一世纪最先进的空间探测技术之一，无线电掩星探测技术不同程度地克服了气象学领域中传统观测的一些弊端，可获得全球大气折射率、温度、气压等中性大气参数廓线，是一项准实时、全球覆盖、精度高、垂直分辨率高、成本低的创新性技术。

气象、电离层和气候卫星联合观测系统(Constellation Observing System for Meteorology Ionosphere and Climate,COSMIC)于2006年发射六颗近地轨道卫星，经过产品验证，其中性大气产品具有较高的精度水平[3]，可成为研究大气变化的高精度资料[4]和提升数值预报水平的重要输入[5]。在经过十四年的运行后，COSMIC-1任务卫星于2020年5月1日退役。美国NOAA与中国台湾NSPO共同将六颗FORMOSAT-7/COSMIC-2(F7C2)掩星任务卫星送入倾角为24°的低地球轨道，取代退役的COSMIC-1任务卫星。其上主要有效载荷名为TGRS(Tri-GNSS Radio-occultation System)，能接收GPS和GLONASS两个导航系统的L1和L2信号，并且达到了迄今为止掩星探测载荷的最高信噪比，比COSMIC-1的平均值高出约两倍。高信噪比不仅提高了大气廓线的反演质量，令以往产生较大误差的低纬度带、低对流层样本的精度有所提高，也使得观测数量有了大幅的增加，更有利于科学研究和业务应用。目前国际上掩星计划均存在由于水汽影响而造成廓线精度不同程度下降的问题，使得掩星数据在低对流层的应用受限。利用高信噪比载荷观测的COSMIC-2理论上可提升这方面的数据精度。但目前国内对于COSMIC-2掩星数据的反演与分析研究还较少，尚未获得很好的应用。因此，对于COSMIC-2数据反演以及质量验证和分析，尤其是对低纬度带、低对流层的廓线样本精度进行评价，显得尤为重要。

本文基于CDAAC(COSMIC Data Analysis and Archive Center)公布的2019年10月—2020年9月COSMIC-2 conPhs大气附加相位L1b数据，利用无线电掩星数据处理软件ROPP(the Radio Occultation Processing Package)反演大气廓线。首先采用几何光学和波动光学方法反演大气折射率廓线，再通过折射率一维变分同化得到温湿廓线。以探空站点资料作为参照，将反演的折射率、温湿廓线结果分不同纬度、不同季节进行了统计验证和精度分析。同时，与COSMIC-2的wetPf2二级数据进行对比，对反演数据的优劣进行讨论。

2 数据与方法

2.1 数据介绍

本文采用每日发布的conPhs大气附加相位L1b数据作为反演输入数据，日均可获取四千余个掩星观测样本。由于COSMIC-2目前在轨的六颗卫星均处于低倾角(24°)轨道，因此所用的探测数据的地理覆盖范围仅在45°S～45°N纬度带之间。探测的基本参数是GNSS无线电载波信号的振幅、相位以及导航卫星与低轨卫星的位置信息。

2020年5月之后的conPhs数据文件中不再存放以100 Hz频率记录的卫星位置信息，而是每秒仅记录一次。因此，为了对COSMIC-2进行数据反演和分析，本文将文件中的卫星位置坐标分别插值到了100 Hz的频率上，并使用插值后的位置坐标的斜率计算出卫星运行速度，重新创建完整的conPhs文件。

wetPf2数据为COSMIC-2官方公布的二级产品，CDAAC以COSMIC-2掩星探测结果(弯曲角和折射率)作为观测量，以欧洲中期天气预报中心ECMWF发布的再分析数据ERA-Interim作为背景先验信息[6]，使用一维变分方法来反演大气温度、湿度和压强[7]。Chen等[8]以无线电探空仪数据为标准进行精度验证，得出wetPf2产品在0～15 km范围内温度和水汽压标准差分别为1.3～1.5 K和0～2.2 hPa的评价。

2.2 反演流程

当GNSS信号发射机从高空发射的无线电信号穿过大气层时，由于大气中介质垂直密度分布的变化，信号传播路径会发生一定的偏折。偏折前后产生一个弯曲角α，根据这一弯曲角可利用积分变换计算大气折射率。本文所用的conPhs大气附加相位L1b数据可提供单次掩星事件的大气附加相位、信噪比、卫星位置和速度等辅助数据，在经过对附加相位数据的滤波、质控后，采用几何光学和波动光学反演方法进行弯曲角的计算。在局部球对称的假设下，该弯曲角可与碰撞参数a(定义为折射中心点到信号射线的垂直距离)建立起唯一的对应关系。在去除电离层效应并统计优化后，使用Abel积分变化将弯曲角转换为大气折射率。再利用ECMWF再分析数据作为背景场数据，进行掩星折射率一维变分同化反演，得到大气温湿廓线。本文使用ROPP掩星数据处理软件实现上述反演流程。

ROPP掩星数据处理软件由欧洲气象局无线电掩星气象卫星应用机构(ROM SAF)开发，用于掩星观测数据处理与反演，可适用于世界上几乎所有掩星观测计划发布的数据。王大钊等[9]利用ROPP成功地进行了无线电掩星数据反演，与二级产品进行对比，得到的折射率和湿度相对误差在2%以内，温度误差不超过2 K。参照其反演流程进行实验设计，给出所用软件模块进行反演的主要流程图(图1)。

图1 反演流程图

在大气折射率反演部分，主要使用处理模块中的ropp_pp_occ_tool工具[10]得到折射率结果。该部分首先对L1和L2双频信号观测的附加相位和振幅进行滤波降噪等质量控制[11]，在信号跟踪误差非常大的情况下，会对数据进行截断并舍去。相比于L1信号，L2信号的信号跟踪误差更大，利用MSIS全球气候学模型计算得到的模型信号值来代替L2信号[12]。对预处理后的双频信号采用几何光学[13]和波动光学(正则变换CT2算法[14])方法反演弯曲角廓线，两种方法反演的结果在25 km左右高度处拼接。25 km以下采用波动光学反演结果，可更好地减小在低对流层由于水汽增多引起的多路径效应影响。为了去除弯曲角廓线结果中存在的电离层噪声贡献，应用电离层校正算法，结合实测和气候弯曲角数据，通过统计优化来降低噪声[11]。最后，利用Abel积分变换[15]的逆变换公式将弯曲角α转换为折射率N。

在大气温湿反演部分，主要使用折射率一维变分模块中的ropp_1dvar_refrac工具[16]得到大气温湿结果。对以下代价函数进行求解：

其中，x为大气状态向量的最大似然解；观测场yo为反演得到的大气折射率数据；背景场x b选择137层的ERA5模式层数据；背景误差协方差B和观测误差协方差O选用ROPP提供的对应层数的误差信息文件；正演算子h(x)的基本表达式采用Smith公式[17]。尽管许多学者基于Smith公式进行了更多研究，以得到精度更高的表达式，但Smith公式仍然是GPS气象领域广泛采用的折射率“标准”表达式[18]。利用基于拟牛顿法[19]的ROPP极小化工具对代价函数进行极小化求解，得出最优大气状态参量。

2.3 反演结果检验方法

利用美国怀俄明大学网站提供的探空站点资料作为参考真值，对反演结果进行对比验证[20]。由于探空站点数据时间相对固定(每日00 UTC和12 UTC)，而掩星观测数据的时空分布较为随机，所以对两者进行统计比较时，观测时间和位置必须匹配至合理范围内。匹配原则参考Kuo等[20]的做法，首先将两者观测时间差控制在±1 h内，再以探空站点的位置为中心，周围100 km范围内匹配掩星样本点。图2显示了2019年12月21日单日样本的时空匹配情况，将匹配好的数据对进行位势高度上的插值。设定高度层范围为0～35 km，分辨率为0.1 km，共350层。对于在整个高度层上存在部分缺测的样本数据，仍保留廓线，利用有数据的部分高度层加入统计。2019年10月—2020年9月共得到3 689对匹配样本。

图2 2019年12月21日COSMIC-2与探空站点数据全球匹配情况

进行对比验证的大气参量为大气折射率、大气温度和相对湿度廓线。由于探空站点没有直接的折射率数据提供，故在对COSMIC-2反演折射率数据进行质量分析之前，先将探空资料由提供的温度、气压、水汽混合比转换为折射率参数：

其中，r为水汽混合比(g/kg)，P为大气压强(hPa)，e为水汽压强(hPa)，T为温度(K)，N为大气折射率。

同时，将反演的比湿廓线结合温度和气压，转换为相对湿度廓线：

其中，T为温度(K)，E为温度T时的饱和水汽压(hPa)，P为大气压强(hPa)，q为比湿(g/kg)，e为水汽压强(hPa)，RH为相对湿度。

为保证数据质量，在对比分析前对探空站点和反演数据分别进行质控。对折射率数据的质控手段主要有两步。第一步进行极值检查：根据大气折射率极值范围，将存在折射率值小于0或大于400的廓线样本剔除；第二步进行单调性检查：根据大气折射率与高度的指数关系[21]，若单个廓线样本中，折射率随高度增加而不单调递减的高度层超过总高度层数的20%，则剔除该廓线。经过上述质控，剔除约14%的样本量。大气温度和相对湿度采用折射率质控后的对应样本进行统计。

在对比分析中，折射率采用平均相对偏差MRE和相对偏差标准差STD来验证精度，分别由公式(7)、公式(8)和公式(9)得到；大气温度和相对湿度采用平均偏差ME和均方根误差RMSE来验证精度，分别由公式(10)和公式(11)得到。i表示第i个廓线样本，n为统计样本数，x代表统计量，即大气温度和相对湿度，RO和Sound分别代表掩星反演得到的和探空站点的数据。

在分纬度带对比分析中，将反演结果样本分为南北纬0°～15°、15～30°、30～45°三个纬度带；在分季节对比分析中，将反演结果样本根据掩星点的纬度区分为南北半球，并以不同月份划分四季：北半球3—5月和南半球9—11月为春季，北半球6—8月和南半球12—2月为夏季，北半球9—11月和南半球3—5月为秋季，北半球12—2月和南半球6—8月为冬季。

3 反演结果与分析

3.1 大气折射率廓线

3.1.1 整体验证

图3和图4从整体上验证了COSMIC-2数据反演的大气折射率的精度。图3为探空数据与反演折射率数据在0～35 km范围内选取的九个高度层上的散点对比图，各层的数据点分布在y=x参考直线两侧，大部分高度层上反演折射率比探空折射率略小，相关系数在0.663～0.967之间，其中15 km高度层最优，在2 km和35 km高度附近负偏差较明显。2 km及以下高度大气中超折射现象增多，应用算法中的Abel积分变换导致了负偏差[22]；而35 km高空则可能在电离层校正和降噪过程中对弯曲角的不当加权，导致优化的弯曲角存在一定偏差，这一偏差最终进入到了折射率结果中[23]，这与折射率相对偏差在平流层中通常随高度增加而增大的事实[22]相符。总体上各高度层偏差不大，说明反演折射率与探空折射率之间有很好的一致关系。

图3 探空数据(横坐标)与COSMIC-2反演折射率(纵坐标)不同高度层处散点图 单位：N-units。

图4为反演大气折射率与wetPf2二级产品折射率数据的平均相对偏差(MRE)和相对偏差标准差(STD)图，retri代表反演折射率，wetL2代表二级产品折射率，样本数量为掩星样本和探空样本在不同高度上匹配后得到的样本分布情况。从反演折射率的整体统计结果来看，MRE在-2.5%～0.5%以内，整体上反演折射率略低于探空站点折射率，5 km以下和25 km以上均以负偏差为主，低层最大负偏差为1 km处的-2.2%，随高度升高至5 km处的0.3%，在5～25 km高度层上始终稳定保持在±0.5%以内，25 km以上又显示出负偏差的趋势；STD整体上在1.5%～4.3%，最大出现在1 km处为4.3%，随高度升高逐渐减小，至35 km都维持在3%以内，最小值出现在25 km处为1.5%；不同高度层的统计样本数量在5～10 km内最多且稳定，5 km以下样本数量减小是因为掩星探测高度有限，在一些多水汽的情况下无法探测到低层；而由于探空折射率数据在质控后越往高层越容易出现缺测误测的情况，导致了10 km以上样本匹配数量的减少。

图4 折射率廓线平均相对偏差和相对偏差标准差图

从反演折射率与产品折射率的对比上来看，两者都保持在较高精度水平且都存在一定的负偏差，但反演折射率5 km以下的负偏差比产品折射率更大，在5～30 km内更小；在24 km以上反演折射率STD明显小于产品折射率STD，而在24 km以下差距不大。造成这种结果的原因可能是对于L2信号的处理不同：在ROPP中采用MSIS模型作为L2信号的模型代替值，CDAAC则采用CIRA+Q模型[24]，两种模型之间的差异可能导致了折射率STD的差异。同时，25 km附近是反演折射率两种方法结果的拼接高度，这一高度以下当L2信号噪声过大时将不再使用L2信号加入反演过程中，因此仅在24～35 km中差异较大。两者统计样本数量在5 km以上趋近一致，但在5 km以下产品折射率样本的数量更多。

3.1.2 分纬度验证

图5为反演的大气折射率按纬度划分的平均相对偏差和相对偏差标准差图，从分纬度带对比平均相对偏差中可以看出，在5 km以下，各纬度带样本均有不同程度的负偏差，且纬度越低，负偏差越大，最大负偏差出现在15°纬度带样本中，为-4%；在5～25 km内各纬度带样本差距不大，三类样本的MRE均在±0.5%以内；25 km以上的负偏差也以15°纬度带样本最为明显。从分纬度带对比相对偏差标准差中可以看出，5 km以下15°纬度带样本的STD最小，30°纬度带样本的STD最大；在5 km以上均显示出纬度越低，STD越小的趋势，特别是在8～16 km高度层上，三类样本STD的差距相对较大。因此可以看出，COSMIC-2低纬度样本的折射率精度都一定程度地优于中纬度样本，且在不同高度层上呈现不同差距。

图5 折射率廓线按纬度划分的平均相对偏差和相对偏差标准差图

3.1.3 分季节验证

图6为反演的大气折射率的按季节划分的平均相对误差和相对偏差标准差图，从分季节对比平均相对偏差中可见，夏季在5 km以下的低对流层负偏差最明显；5～15 km四类样本差距不大，均在±0.3%内；25 km以上的负偏差又以夏季最明显。从分季节对比相对偏差标准差中可见，18 km以下夏季折射率STD始终为四类样本中最大，秋季样本次之，冬季和春季最小；尤其在3 km以下，夏季折射率STD最大可达5.5%，冬季样本最优；而在20 km以上，四类样本STD变化则显示出较高的一致性。出现以上现象的主要原因是，在夏季低层尤其是低对流层水汽较多，造成了较大的折射率梯度，对掩星无线电信号产生较大影响，导致了相对较大的误差。

图6 折射率廓线按季节划分的平均相对偏差和相对偏差标准差图

3.2 大气温度廓线

3.2.1 整体验证

图7和图8从整体上验证了COSMIC-2数据反演的大气温度的精度。图7为探空数据与反演温度数据在0～35 km范围内选取的九个高度层上的散点对比图。反演温度与探空温度在大部分高度层上偏差不大，相关系数总体在0.677～0.939之间，2 km处最优，25 km处最差，相关系数的下降可能与掩星观测和探空观测越往高层水平漂移越大有关。对流层温度误差相对更小，RMSE可达到1.5 K以内；平流层温度误差相对较高，但整体上RMSE均不超过2.5 K。

图7 探空数据(横坐标)与COSMIC-2反演温度(纵坐标)不同高度层处散点图 单位:K。

图8为反演大气温度与wetPf2二级产品温度数据的平均偏差(ME)和均方根误差(RMSE)图。从反演温度的整体统计结果来看，ME整体上在±0.5 K以内，低层0～2 km以负偏差为主，最大为0.3 km处的-0.4 K，随高度升高至4 km处的0.2 K，在4～12 km高度层上稳定在0.2 K左右，之后开始减小，从15 km开始往上到30 km再次变为负偏差，最大达-0.3 K；RMSE整体上在1.3～2.2 K，0～14 km高度处稳定保持在1.3～1.5 K，随后增大至25 km处的2.1 K；不同高度层的统计样本在2～16 km内数量最多且稳定，这个范围也是温度RMSE最小的高度范围。从反演温度与产品温度的对比上来看，两种数据在各高度层上误差趋势很接近，反演温度的ME绝对值比产品温度的更小，在0～13 km和27～35 km高度上反演温度的RMSE比产品温度更小，且反演温度的样本数量在5 km以下的底层更多。因此，反演温度的精度在大多数高度层上略优于二级产品温度，这可能是两种数据进行一维变分同化时所选取的背景场信息不同所导致的。

图8 大气温度廓线平均偏差和均方根误差图

3.2.2 分纬度验证

图9为反演的大气温度按纬度划分的平均偏差和均方根误差图。从分纬度带对比的平均偏差中可以看出，样本点的纬度越高，ME的绝对值越大；而在17 km以下，15°纬度带的样本以负偏差为主，30°和45°纬度带的样本以正偏差为主，17 km以上则大致相反。从分纬度带对比均方根误差图中也可以更明显地看出，样本点的纬度越高，RMSE越大，廓线精度也越差；三类样本的RMSE在17 km以上走势逐渐趋同，而17 km以下则显示出较大差距；三类样本在同一高度层的RMSE最大差距在10 km左右处，15°纬度带样本RMSE比45°纬度带样本小0.8 K左右。因此可以看出，COSMIC-2低纬度样本的大气温度精度优于中纬度样本，且在对流层差距明显。

图9 大气温度廓线按纬度划分的平均偏差和均方根误差图

3.2.3 分季节验证

图10为反演的大气温度按季节划分的平均相对误差图。从分季节对比中可以看出，24 km以上夏季与秋季的温度ME绝对值略大一些，24 km以下四季样本ME差异较小；15 km以下秋季样本RMSE略小于其他三季，15 km以上四季样本RMSE差异较小。但可以看出，各类样本在各高度层上的差异较小，温度的季节性差异不显著。

图10 大气温度廓线按季节划分的平均偏差和均方根误差图

3.3 相对湿度廓线

3.3.1 整体验证

图11和图12从整体上验证了COSMIC-2数据反演湿度的精度。由于大气中的水汽主要存在于低层，因此本文只统计0～12 km高度范围内的湿度情况。图11为探空数据与反演湿度数据在0～12 km范围内选取的六个高度层上的散点对比图。反演相对湿度与探空相对湿度各高度层上的相关系数在0.784～0.903，精度随高度升高逐渐变差，但最大均方根误差不超过15%。

图11 探空数据(横坐标)与COSMIC-2反演湿度(纵坐标)不同高度层处散点图 单位:%。

图12为反演相对湿度与wetPf2二级产品相对湿度数据的平均偏差(ME)和均方根误差(RMSE)图。从反演相对湿度的整体统计情况来看，ME在±5%范围内变化，随着高度的升高先由正偏差减小为负偏差，再变为正偏差，且4 km以下的低层偏差绝对值要大于4 km以上的部分；RMSE随着高度升高逐渐变大，但始终保持在15%以内，具有较高的精度。从反演湿度与产品湿度的对比情况可见，3 km以下两者ME呈相反走势，但都在±5%以内，3 km以上反演相对湿度ME明显小于产品相对湿度；从RMSE上也可看出，反演相对湿度结果优于二级产品，在6 km以上尤为明显，两者差距最大可达7%左右；统计样本数量呈相同趋势分布，且反演样本比产品样本多约50%，但湿度样本数量相对于温度来说更少，这可能与湿度反演过程中的质控手段[7]有关。因此，反演相对湿度的精度在整体上都优于二级产品相对湿度，与温度结果在0～13 km之间反演精度优于二级产品的结论类似。

图12 相对湿度廓线平均偏差和均方根误差图

3.3.2 分纬度验证

图13为反演相对湿度按纬度划分的平均偏差和均方根误差图，2 km以下相对湿度ME在不同纬度样本中趋势类似，在2～5 km高度内15°和30°纬度带样本负偏差更明显，在5 km以上两者ME绝对值比45°纬度带样本都要小；在5 km以下，15°纬度带样本在三类样本中RMSE最小，且在近地面可达10%以内，5 km以上三类样本则显示类似的误差趋势，差异较不明显。因此可见，以往掩星探测结果中因水汽造成的较大误差被缩小，尽管差异不大，但低纬度样本在低对流层显示出了更好的精度。

图13 相对湿度廓线按纬度划分的平均偏差和均方根误差图

3.3.3 分季节验证

图14为反演相对湿度按季节划分的平均偏差和均方根误差图。其中冬季样本的ME绝对值总体上保持最小，夏季最大，其余两季样本ME变化差异不大；在4 km以上冬季样本的RMSE略小于其他三季样本。因此可见，反演相对湿度样本在冬季的误差略小，其余样本误差变化的季节性差异并不显著。

图14 相对湿度廓线按季节划分的平均偏差和均方根误差图

3.4 误差原因分析

本文对COSMIC-2数据的反演结果相比于二级产品，在一定程度上取得了更高的精度，但相比于探空观测资料还有一定的误差存在。造成误差的原因可能有：(1)反演数据与探空站点资料时空匹配时仍存在时间和空间位置上的差异。探空气球上升探测所形成的廓线并不是严格垂直的，存在一定的漂移，越往高空漂移量越大；同时，掩星观测的廓线也存在位置上的水平漂移，实际处理数据中35 km高度范围内最大水平漂移有10～50 km，两者在水平漂移方向上的不确定性也会造成一定的误差；(2)在水汽含量很大的低对流层通常会产生超折射，在此处采用Abel积分变换导致了低层折射率的负偏差。

4 结 论

本文使用COSMIC-2 conPhs大气附加相位L1b数据，基于ROPP软件计算得到大气折射率和大气温湿廓线，对上述大气廓线进行对比验证和精度分析，得出以下主要结论。

(1)COSMIC-2反演大气折射率、大气温度、相对湿度与探空站点数据对比，在不同高度层上均具有较好的一致性，相关系数最高分别可达0.967、0.939、0.903；反演折射率相对偏差标准差为1.5%～4.3%，反演温度和相对湿度均方根误差分别为1.3～2.2 K、10%～15%；反演折射率的精度在24 km以上优于二级产品折射率，反演温度的精度在0～13 km和27～35 km上优于二级产品温度，反演相对湿度在整体上都优于二级产品相对湿度。

(2)以探空站点资料为参考，COSMIC-2反演折射率和温度的精度随纬度变化显著，低纬度样本精度总体上均优于中纬度样本；反演相对湿度在3 km以下低纬度样本精度优于中纬度样本，3 km以上差距不大。表明COSMIC-2数据对低纬度带观测的改善显著，使其精度大大提升，成为各纬度带中精度最高的样本。

(3)以探空站点资料为参考，COSMIC-2反演数据的精度随季节变化不明显，总体上冬、春季样本误差略小于夏、秋季。表明COSMIC-2数据减小了由水汽含量带来的误差影响，缩小了夏季样本与其余季节样本的精度差距。

致 谢：感谢COSMIC数据分析中心(CDAAC)提供的COSMIC-2掩星观测数据，感谢南京信息工程大学高性能计算中心的计算支持和帮助。