1 引言

云对地气系统能量收支平衡具有调节作用，云的时空变化和辐射性质是影响全球气候和各个尺度天气系统的重要因子。目前，云顶气压和有效云量等云参数已被广泛应用于数值天气预报云初始化过程中(Wu and Smith, 1992；Hu et al., 2006；Pangaud et al., 2009；丁伟钰等，2010；刘瑞霞等，2013；Zhu et al., 2016)。然而由于大气环流的作用，云在水平和垂直方向瞬息万变，在海洋和陆地的分布范围也有较大差异，因此构建一个快速有效的实时云反演方法是提供实时云产品分析的关键(Auligné et al., 2011)。当前极地轨道卫星和地球静止轨道卫星能联合提供高光谱分辨率和高水平分辨率的数据流(董佩明等，2008)。卫星测量到的云特性可以为有限的地面观测提供有益的补充，进而实现云观测的连续性和全球性(Rossow et al., 1993)。

有关卫星资料云反演研究，国内外学者分别基于云图和卫星辐射率资料已经做了大量卓有成效的工作。关于云图识别，Goodman et al.(1998)总结了三种基本的云检测方法：阈值法、统计学方程以及辐射传输技术。王立志等(1998)和马芳等(2007)分别使用四通道和多通道的云图云检测方法对云量进行了识别。相对于云图资料，卫星辐射率是更为原始的观测资料，能提供细致的云量垂直分布信息。Chahine(1975)提出了一种计算观测和模拟辐射率值之间的均方根误差来确定云顶高度、云量的方法。Menzel et al.(1983)通过CO₂切片法计算云顶气压和有效发射率，进而来计算静止卫星云导风。Smith and Frey(1990)通过CO₂切片法计算云顶气压和有效发射率并对先进TIROS业务垂直探测器大气红外探测器进行云检测。官莉(2007)使用与AIRS(the Advanced InfraRed Sounder)同步的MODIS(MODerate-resolution Imaging Spectro-radiometer)二级云量产品来确定受云影响的视场。陈靖等(2011)基于Goldberg et al.(2003)的AIRS云检测方案，结合我国的GRAPES-三维变分系统和AIRS仪器特征，分别对海洋、陆地的视场进行云检测。Huang et al.(2004)提出了最小局地发射率方差方法，通过计算背景场中特定气压层的云谱发射率的局地方差来得到单层云发射的最优估计。同样基于变分方法，Auligné et al.(2014a，2014b)和Xu et al.(2013, 2015)提出了多元极小残差(Multivariate and Minimum Residual，简称MMR)。MMR云检测方法通过构建代价函数，采用极小化方法来拟合观测，模拟得到各个模式层的云量参数。根据晴空下的辐射率和模拟的全天空条件下的辐射率值差异来确定通道是否被污染。Xu et al.(2015)基于WRFDA[Weather Research and Forecasting (WRF) model data assimilation system；Barker et al., 2012]对MMR反演方法展开了系统性地评估，并验证了该方法的稳定性和合理性。近年来，基于在状态空间传播的随机样本来近似得到概率密度函数的思想[如粒子滤波(Particle Filter，简称PF)和集合卡尔曼滤波(Kalman Filter，简称KF)等]已经被广泛应用到同化(Hamill et al., 2011)和反演大气波导(盛峥等，2012)等多个气象领域。因此本文尝试把概率论中PF算法的思想应用到气象中的卫星资料的云反演领域中。本文基于格点统计插值分析系统(Gridpoint Statistical Interpolation analysis system，简称GSI)，以GOES(Geostationary Operational Environmental Satellites)-Imager资料为例，采用PF算法进行了云反演数值试验，初步评估了PF云反演方法的可行性及其与MMR云反演方法的异同。

2 反演方法、模式与数据 2.1 粒子滤波(Particle Filter)

PF方法(Arulampalam et al., 2002)的核心思想是基于蒙特卡洛方法，即利用粒子样本来表示概率分布。粒子滤波方法通过从后验概率中抽取的随机状态粒子来表达其分布，是一种顺序重要性采样方法。这里的样本即指粒子，当样本数量n趋近于无穷大时可以逼近任何形式的概率密度分布。尽管算法中的概率分布只是对云真实分布的近似，但由于其非参数化的特点，有效避免了在处理非线性滤波问题时随机量必须满足高斯分布的条件制约，进而能表达比高斯模型更广泛的分布，也对变量参数的非线性特征有更强的建模能力。

在云反演中使用PF算法，首先需要构建初始的云量廓线粒子。如图 1所示，c = c¹, c², ..., cⁿ代表一个视场中各个模式层的有效云覆盖比例(cloud fraction)，c⁰为晴空比例，n为模式层数。针对某一个卫星视场，所有模式层的云覆盖比例以及晴空区域所占比例之和应该为1，即c¹+c²+...+cⁿ+c⁰=1。对于该视场的任意高度面，云覆盖的区域阻挡了来自该平面以下区域比例为c^k的向上辐射，剩余的1－c^k的向上辐射有效。

在PF云反演过程中，定义每个粒子为一组完整的$ {c^1}, {c^2}, \ldots, {c^n} $组合。首个时次反演所选用的粒子总数目为n+1。其中，第一个粒子c¹=1，即所有云均在第一个模式层；第二个粒子c²=1，即所有云均在第二个模式层；以此类推，第n个粒子cⁿ=1，而第n+1个粒子c⁰=1，即完全晴空。某个模式层上的粒子c^k的后验概率为$ \exp \left({ - {{\left({\frac{{R_v^{{\rm{obs}}} - R_v^k}}{\sigma }} \right)}^2}} \right) $，因而粒子c^k的更新方程为

$ {c^k} = {c^k} \cdot \exp \left({ - {{\left({\frac{{R_v^{{\rm{obs}}} - R_v^k}}{\sigma }} \right)}^2}} \right), $

(1)

其中，$ \forall k \in \left[ {1, n} \right] $，σ为观测误差，R_v^obs为在波数v处观测的辐射率，R_v^k为假设在模式层第k层处放入黑体云时，辐射传输模式在波数v处模拟的辐射率值。$ \exp \left({ - {{\left({\frac{{R_v^{{\rm{obs}}} - R_v^0}}{\sigma }} \right)}^2}} \right) $为晴空比例c⁰的后验概率，因此c⁰的更新方程为

$ {c^0} = {c^0} \cdot \exp \left({ - {{\left({\frac{{R_v^{{\rm{obs}}} - R_v^0}}{\sigma }} \right)}^2}} \right), $

(2)

其中，R_v⁰是晴空条件下，辐射传输模式计算得到的辐射率值。为了使各个模式层的云覆盖比例和晴空比例总和为1，在对粒子进行更新后，对所有粒子采用如下公式进行归一化处理：

$ {c^k} = \frac{{{c^k}}}{{\sum\limits_{k = 1}^n {{c^k}R_v^k} }}, $

(3)

其中，$ \forall k \in \left[ {0, n} \right] $。对于首个时次的反演(同化的背景场中没有前一个时次的云反演结果，即冷启动)，设置云等概率均匀地分布在各个模式层，即

$ {c^k} = \frac{1}{{n + 1}}, $

(4)

其中，$ \forall k \in \left[ {0, n} \right] $。试验其他时次的背景场由前一个时次反演结果经过模式动力模块中的平流扩散过程得到。

PF反演方法是解决概率问题的有效手段，然而该算法仍然存在着一些问题。其中最主要的问题是需要采用大量的样本数量才能很好地近似系统的后验概率密度。概率问题越复杂，描述后验概率分布所需要的样本数量就越多，算法的复杂程度就越高。另外，重采样阶段会造成样本有效性和多样性的损失，导致样本贫化现象的发生。如何保持粒子的有效性和多样性，克服样本贫化，也是该算法研究重点。前面介绍的粒子滤波的云反演方法采用的粒子数目为模式的垂直层数加1。试验中粒子样本数目为40(模式的层次)，故粒子数为41。粒子滤波容易出现滤波发散问题。针对PF反演方法有限粒子样本数的缺陷，进一步提出并设计了一种新的基于扰动的PF方法。即考虑在云反演过程中在背景场的云量廓线的基础上，通过上下移动和增大缩小云覆盖比例等扰动方法得到更多的粒子样本。其中，q和p分别表示上下移动幅度和缩大或放小的系数，如：

$ {c^k}\left({p, q} \right) = p \cdot {c^{k \pm q}}, $

(5)

其中，$ \forall k \in \left[ {1, n} \right] $。因此扰动后的新样本数目为n+1+2p·q。本文试验中移动幅度q为上下5个模式层，p为缩小到50%、55%、…、95%，一直到扩大到150%(以5%为间隔)。增加样本后的粒子数目为251。这种通过上下移动和缩放云量的粒子滤波方法，以下简称为APF(Advanced Particle Filter)。

2.2 MMR多元极小残差法

多元极小残差云检测方法MMR(Auligné, 2014a，2014b；Xu et al., 2015)，通过逐步拟合辐射率观测值来得到每个模式层(假定n个模式层)的有效云覆盖比例，其拟合的全天空辐射率值R_v^cloud为

$ R_v^{{\rm{cloud}}}\left({{c^0}, {c^1}, {c^2}, \ldots, {c^n}} \right) = {c^0}R_v^0 + \sum\limits_{k = 1}^n {{c^k}R_v^k}, $

(6)

其中，

$ \left\{ \begin{array}{l} 0 \le {c^k} \le 1, \;\;\forall k \in \left[ {0, k} \right]\\ \;\;\;\;\;\;{c^0} + \sum\limits_{k = 1}^n {{c^k} = 1} \end{array} \right.. $

可见，公式(6)是基于云参数$ {c^0}, {c^1}, \ldots, {c^n} $线性化的全天空辐射率观测算子。通过构建以下代价函数, 采用极小化变分方法，逐步拟合辐射率观测值来得到n个模式层的云参数信息$ \mathit{\boldsymbol{c}} = {c^1}, {c^2}, \ldots, {c^n} $：

$ J\left({{c^0}, \mathit{\boldsymbol{c}}} \right) = \frac{1}{2}{\sum\limits_v {\left[ {\frac{{R_v^{{\rm{cloud}}} - R_v^{{\rm{obs}}}}}{{R_v^{\rm{0}}}}} \right]} ^2}. $

(7)

2.3 模式介绍

本文采用一个完全可压缩、非静力中尺度模式WRF(Skamarock et al., 2008)作为系统预报模式。该模式水平方向采用荒川C网格，垂直方向采用随地形的质量坐标。同化系统为格点统计插值(Gridpoint Statistical Interpolation，简称GSI)分析系统。该分析系统最初由美国国家海洋和大气局NOAA(National Oceanic and Atmospheric Administration)、国家环境预报中心NCEP(National Centers for Environmental Prediction)在其早期业务运行的波谱统计插值分析系统的基础上研发的下一代区域及全球数值预报分析系统。GSI分析系统于2006年6月和2007年5月分别被应用在NOAA的北美区域资料同化系统(NDAS)和全球资料同化系统(GDAS)中，并实现了业务化运行。GSI数值模式自带辐射传输模式CRTM(the Community Radiative Transfer Model；Han et al., 2006；Liu and Weng, 2006)，使用的版本为CRTM2.2.3，具有很好的开发性、兼容性和稳定性。由于CRTM模式是从0.1~1013 hPa，而预报模式是从10~1000 hPa，对于预报模式顶层以上(＜10 hPa)及底层以下(＞1000 hPa)，用线性外推法，试验中采用辐射传输模式自带参数文件提供的参考值。对于臭氧廓线采用气候值，地面参数也取自背景场，地面温度取与2 m温度值相同。

2.4 数据介绍

美国NOAA静止轨道业务卫星系列GOES，采用双星运行体制。GOES-East卫星和GOES-West卫星分别定点在75°W和135°W赤道上空，覆盖范围为20°W~165°E，占近地球1/3的面积。GOES卫星有5个通道(一个可见，四个红外)，5个通道的中心波长分别为0.62、3.9、6.5、10.6和13.3 μm。通道2、4、5的分辨率为4 km，而通道1和3的分辨率为8 km。本文使用美国国家航空航天局NASA兰利研究中心发布的GOES产品(http://www-angler.larc.nasa.gov [2016-01-10])中的辐射率资料作为反演的输入资料。选取2013年12月15日00:00(协调世界时，下同)至2013年12月15日06:00进行逐小时的云反演试验。云反演试验中2013年12月15日00:00的气象场由1°×1°分辨率的NCEP再分析资料启动，背景场中云均匀分布在各个模式层上，初始的粒子数为n+1(包含由2.1节中介绍的各个模式层的单层云以及全晴空状态)。其他时次的气象场由前一个时刻的气象场经过1 h预报得到，云图背景信息由前一个时次反演结果经过模式的平流扩散过程得到，初始粒子构造方法同前。具体流程如图 2所示。

3 反演试验与结果分析

红外辐射率观测资料的预处理、质量控制以及云反演试验均基于GSI系统展开。本文选取2013年12月15日00:00至2013年12月15日06:00 7个时刻的不同类型的云天个例进行反演试验。试验的背景场由WRF模式预报得到，网格分辨率为15 km，水平格点数为415×325，模式垂直层数为40层，模式顶层气压为50 hPa。Xu et al.(2013)的研究结果发现大气辐射传输模式基于WRF模式的背景场对长波CO₂波段的辐射值的模拟效果较好，同时15 μm长波波段通道的卫星观测也是云顶高度反演的主要通道。因而本文选取各种辐射率资料的长波波段通道(波长约15 μm)来反演云参数。

使用NASA兰利研究中心发布的1°×1°的GOES云顶气压和云底气压作为参考场来检验反演的云图产品质量。尽管NASA-GOES云产品并不是云分布的实况，但是该产品是基于发展较为成熟的CO₂切片法(Menzel et al., 1983)得到了的，且已被各大业务部门广泛使用。因而NASA-GOES云产品在一定程度上反应了检验时刻的云分布形势，为验证反演结果提供了参考依据。标记NASA-GOES云产品中存在云顶气压和云底气压记录的区域为有云区域为1，标记本文试验反演的整层存在大于1%云覆盖比例的区域为1。NASA-GOES云产品和本文试验反演的云产品同时插值到0.1°×0.1°的网格上(有云标记为1，晴空处为0)，再进行评分。

3.1 云掩膜结果

图 3给出了NASA-GOES云产品和本文中MMR方法(或PF方法)反演云产品插值到统一网格点后的云掩膜(Cloud mask)，其中红色为有云(标记为1的区域)，晴空为白色(标记为0的区域)。通过逐点对比NASA-GOES云产品和本文试验反演的云掩膜进而得到最终统计结果。图 4显示了基于GOES-Imager辐射率资料反演得到的云覆盖量相对于GOES云产品的反演正确(Hits)，误报(False alarms)和漏报(Misses)的分布形势。可以发现，两种方法基于GOES-Imager辐射率资料反演的云产品和GOES反演的云产品在大部分地区保持较好的一致性(Hits在整个区域中占的分布最大)。MMR方法辐射率资料得到的云产品高估了云的反演，而PF反演方法由于较少的粒子样本数，漏报区域较多。

图 4为MMR方法、PF方法以及APF方法基于GOES-Imager辐射率资料反演得到的云覆盖量相对于GOES云产品的反演正确(Hits)，误报(False alarms)和漏报(Misses)的分布形势。由图 4a和b可见，MMR方法相对于PF方法，辐射率资料得到的云产品存在较多的误报区域。从图 4b、c可以发现，而PF反演方法由于较少的粒子样本数，漏报区域较多。APF方法基于GOES-Imager辐射率资料反演的云产品和GOES反演的云产品在大部分地区有很好的一致性。APF方法基于充足的粒子样本，较大程度改善了PF反演方法漏报现象，反演正确的区域增多(如椭圆区域)。

3.2 评分结果

图 5为2013年12月15日00:00至2013年12月15日06:00，3种云反演方法基于GOES-Imager资料反演的云覆盖量相对于GOES云产品的ETS评分(Equitable Threat Scores；Schaefer，1990)。可见，MMR方法反演产品得到的评分随着反演次数的增多，评分逐步提高；PF方法在第一个时次的评分高于MMR方法，但是由于有限的样本数，经过多个时次后出现滤波发散问题，评分逐渐降低。APF方法则通过增加粒子样本数目，在整个反演区间的ETS评分要普遍高于MMR反演方法和PF反演方法得到的ETS评分。

3.3 云量廓线结果

图 6a显示的是生成图 6b、c、d中剖面的剖面线和模式区域地形高度的示意图，图 6b、c、d给出了2013年12月15日06:00 MMR、PF和APF这3种方法云反演产品的剖面图。MMR方法的云反演产品普遍存在云顶过高的现象。PF反演方法相对于MMR反演方法则能更加清晰地反演出云量的廓线。随着有效样本数目的增加，APF反演方法相对PF反演方法则能更好地捕捉微小的云信息，进而有效避免了滤波发散问题。

3.4 云顶气压结果

图 7给出了2013年12月15日06:00 GOES云产品参考场以及分别利用MMR、PF以及APF方法反演得到的云顶气压。云顶气压由存在大于1%的云覆盖比例的最高模式层的气压值得到。由图 7a可见，GOES云产品的在美国西部海域和东部海岸(如图中方框中所示)有较多低云(如950 hPa)。而MMR方法、PF方法以及APF方法反演的结果均存在漏报低云的现象，这与Auligné(2014a)检测出低云较差的观点相一致。MMR反演方法利用大多数传感器辐射率资料得到的云反演产品与GOES云产品，整体对应较差，整场存在云顶过高的现象，这与Xu et al.(2015)基于GOES-Imager资料反演云高结果较差的结论相一致。PF反演方法相对于MMR反演方法反演的云顶气压更加接近参考场。由于粒子数目的增加，APF方法反演的云顶气压整体反演结果最好。另一方面，对于高云的反演(图中椭圆中所示)，PF方法整体较差。APF反演方法对于高云的反演的云高要优于PF反演方法。

4 结论

基于GSI设计了PF云反演方法，并与MMR云反演方法进行了对比研究试验，结果表明：

(1) 两种方法利用静止卫星GOES-Imager辐射率资料反演得到云的分布和云顶高度和现有的NASA利用CO₂切片法得到的GOES云产品较为一致。两种方法的反演产品的优点是云图信息是三维分布的，相对于现有的大多数云产品能提供更全方位立体的云信息。其中粒子滤波PF方法通过近似云分布的后验概率来更新云背景，相对于多元残差MMR反演方法，其计算效率提高了80%~90%，为未来的准业务化运行奠定了良好基础。

(2) 针对PF反演方法样本匮乏的问题，进一步提出并设计了APF反演方法，通过扰动云覆盖比例值的方式来增加现有的粒子样本数目，从而更好地接近后验概率密度，有效避免粒子滤波发散问题。

(3) PF方法相对于多元残差方法对于云顶气压的反演结果更加接近参考场。增加粒子样本后(样本数量约为250)，APF方法反演的云顶气压整体反演结果最好。3种方法均有漏报低云现象，对于高云的反演，MMR方法和PF方法整体较差，APF方法效果最好。