数值天气预报(Numerical Weather Prediction, NWP)是一个初/边值问题。卫星资料同化则为数值预报模式提供尽可能精确初始状态提供了可能(Wang et al, 2014)。高光谱大气红外探测器AIRS(Atmospheric Infra-Red Sounder)搭载在地球观测系统第二颗卫星Aqua上。AIRS共2 378个通道覆盖650~2 700 cm^-1红外谱区域, 其光谱分辨率远超过之前的高分辨率红外辐射探测仪HIRS(High Resolution Infrared Radiation Sounder), 常被用来探测精细的大气温度、湿度廓线等(Aumann et al, 2003)。AIRS亮温资料已被全球几大数值预报中心进行同化使用, 能提高数值天气预报中温度和水汽分析, 对数值预报模式预报精度有重要的作用(Cameron et al, 2005; Liu et al, 2009), 尤其是云区资料对提高气象敏感区的预报精度有明显的正效果(McNally, 2002)。蔡淼等(2015)基于Cloudsat云检测产品与ECMWF再分析资料, 开展三维云场分布诊断方法研究。马占山等(2016)通过使用卫星云观测资料对全球中期数值预报模式GRAPES_GFS湿物理过程参数化方案中降水场和云宏微观场进行诊断和评估, 以揭示其对云和降水的预报性能。李浩然等(2015)研究中国及周边区域云量和云水路径季节变化并分析大尺度环流对中国云水路径季节性变化的影响。由于红外探测器受云影响, 在探测到云时作为“黑体”处理, 同化AIRS通道亮温时, 首先需进行云检测(王根, 2014)。目前, 云检测主要有基于“晴空视场点”、基于“晴空通道”、“云辐射订正”和“不同仪器间的匹配”四种方法。基于晴空视场点的云检测是寻找绝对晴空的视场点, 当判识出所选通道组合中某个通道在某个视场点(field-of-view, FOV)有云, 则剔除所有通道在此视场点的亮温(English et al, 1999), 常用方法有CO₂分层法(Chahine, 1974; Menzel et al, 1983)、最小剩余法(Eyre et al, 1989)和局部最小方差发射率法(Huang et al, 2004)。基于晴空通道的云检测, 其本质是寻找不受云影响的晴空通道, 代表性方法是McNally and Watts法(McNally et al, 2003)。对AIRS资料来说像元内完全无云的可能性很小, 全球仅约10 %的视场点不受云影响(Huang et al, 2005), 因此经过云检测后, 大量AIRS受云影响通道亮温不能完全被使用(Pangaud et al, 2009)。

经典变分同化基于最小二乘理论要求误差服从高斯分布。最小二乘法对离群值较敏感, 由其数学理论保证, 求得的最终解是离群值和正常值之间的一个“妥协”, 与真实解相差较远(Wang et al, 2014)。因此, 在使用经典变分同化时, 首先需识别数据中离群值并加以处理, 然后再进行同化。因离群值并不总是有害, 且部分离群值往往代表了新的亮点, 如极端天气现象, 这就导致丢失较多有用的数据信息。

目前, 国际上基于云检测方法, 开展了一系列红外亮温资料受云影响同化研究(Pangaud et al, 2009; Pavelin et al, 2008; Martinet et al, 2013; Chevallier et al, 2004)。其基本理论是通过CO₂分层法或最小剩余法等得到视场点云参数, 再把云参数提供给辐射传输模式使用; 或者在每个视场点剔除通道权重函数位于云顶下的单个通道亮温(McNally et al, 2003), 但这些方法都存在大部分资料不被同化使用的局限性。Li et al(2005)提出了同化“晴空辐射”资料, 即使用邻域视场点资料消除云信息重构晴空辐射, 但此方法因假设云是各向同性, 不能满足大多数个例, 导致成功同化无云辐射作用不显著(Pangaud et al, 2009)。董超华等(2013)采用不同仪器云产品匹配得到晴空视场点, 主要通过四维变分直接同化受云影响的亮温辐射, 这对能计算云对亮温辐射值的辐射传输模式和伴随算子要求较高, 超出了当前变分同化系统的计算资源能力(Chevallier et al, 2004)。如果卫星视场点FOV内云参数(有效云量和有效云顶气压)已知或可通过算法进行估算, 那该信息可以被变分同化系统中的辐射传输模式使用, 从而可直接同化受云影响的亮温辐射值。但需要精确的云参数才能有效提高分析场。含云资料的引入使亮温偏差呈现较强的非高斯性, 为充分利用受云影响亮温辐射, 需进行同化算法构建。开展AIRS资料研究为中国未来高光谱大气红外探测资料的同化研究奠定基础, 以期提高NWP预报精度。基于此开展AIRS有云资料同化研究。

使用对非高斯分布具有稳健性的广义变分同化进行AIRS有云亮温同化研究。首先获取云参数, 再将云参数和其它变量加入变分同化极小化迭代, 并提供给辐射传输模式参与通道亮温模拟。其次用AIRS模拟亮温去评估广义变分同化云辐射系统的性能。在执行过程中动态选择AIRS通道形成通道子集, 并在广义变分同化基础上加上约束项以实现云参数“数学解”到“物理解”的过渡。最后将带约束项的广义变分用于AIRS有云辐射进行同化效果对比分析。

带约束项的广义变分同化核心思想是把M-估计法, 如L2-估计、Huber-估计、Fair-估计和Cauchy-估计(Wang et al, 2014)的权重函数耦合到经典变分同化中(Bhar, 2007; Liu et al, 2005)。带约束项的广义变分同化考虑了偏差的非高斯性, 对“离群值”具有较强的稳健性。偏差定义为观测亮温O与模拟亮温B的差值(O-B)。在变分极小化过程中耦合了质量控制, 降低广义尺度范围外观测资料对目标泛函的贡献率, 并能得到较合理的云参数。

Pavelin et al(2008)采用一维变分反演得到的云参数比CO₂分层法和最小剩余法精度高, 但由于一维变分是基于高斯理论变分处理非高斯分布数据, 结果可能存在一定问题(王根, 2014)。在结合Pavelin et al(2008)前期研究成果的基础上, 采用非高斯思想进行红外资料有云同化的研究。

变分同化是基于贝叶斯最优估计理论组合观测和背景场先验信息, 通过极小化目标泛函以得到统计尺度上“物理解”(模式信息和物理过程)的最优大气估计状态(Rodgers, 2000)。为有效同化受云影响的亮温, 在Wang et al(2014)的研究基础上采用带约束项广义变分同化进行研究, 目标泛函定义为:

式中: x为控制变量; x_b为背景场, 又称第一猜测场(包括温度、湿度、压力和风场等信息); B和R分别是背景和观测误差协方差矩阵; 上标T和-1分别表示矩阵转置和逆; y_o为观测亮温; w(r)为权重因子, 本文称为观测项对目标泛函贡献率的重估计; H为观测算子, 对卫星辐射资料而言, 也称为快速辐射传输模式, 在文中采用Tiros业务垂直探测器辐射传输模式RTTOV(Radiative Transfer for Tirosn Operational Vertical Sounder)进行计算(Saunders et al, 2005); P为约束项。

得到的云参数为数学意义上的解。当N_e∉[N_{e_min}, N_{e_max}]和P_CTP∉[P_{CTP_low}, P_{CTP_top}], 此时的解非物理解, 为实现云参数从“数学解”到“物理解”过渡, 定义约束项P为:

式中: σ_CTP和σ_Ne分别为气压和云量代价函数“调节参数”; CTP为有效云顶气压估算; N_e为有效云量估算; P_{CTP_top}和P_{CTP_low}分别为高云和低云气压值; N_{e_max}和N_{e_min}分别为最大和最小云量(即云覆盖)。

带约束项的广义变分同化目标泛函梯度为:

文中▽P采用泰勒展开一阶和二阶近似表示, 如(CTP-P_{CTP_top})³⇒3·(CTP-P_{CTP_top})²+6·(CTP-P_{CTP_top}), 依此类推。

权重因子定义为:

式中: ρ(r)为M-估计法的代价函数(Bhar, 2007), 使用的M-估计见表 1; d为求导数标记; w(r)为以w(r_i)为元素的对角矩阵, r_i=[y_i-H_i(x)]/σ_i; y_i和H_i(x)分别表示第i通道的观测和模拟亮温; σ_i为观测误差协方差矩阵对角元素的平方根(假定各通道不相关)。

表 1(Table 1)

表 1 M-估计法 Table 1 M-estimator M-估计法 代价函数ρ(r) 影响函数φ(r) 权重函数w(r) 广义尺度 L2 r2/2 r 1 ${\rm{Huber}}\left\{ \begin{array}{l}{\rm{if}}\left| r \right| \le c\\{\rm{if}}\left| r \right| > c\end{array} \right.$ $\;\left\{ \begin{array}{l}{r^2}/2\\c\left( {\left| r \right| - r/2} \right)\end{array} \right.$ $\left\{ \begin{array}{l}r\\c{\mathop{\rm sgn}} \left( r \right)\end{array} \right.$ $\;\left\{ \begin{array}{l}1\\c/\left| r \right|\end{array} \right.$ c=1. 345 Cauchy (c2/2){ln[1+(r/c)2]} r/[1+(r/c)2] 1/[1+(r/c)2] c=2.384 9 Fair c2[|r|/c-ln(1+|r|/c)] r/(1+|r|/c) 1/(1+|r|/c) c=1.399 8 广义尺度的取值可结合经验给定或通过稳健统计自动估计

表 1 M-估计法 Table 1 M-estimator

经典变分同化在同化通道亮温前需获得误差分布信息, 若存在离群值, 则最优参数化得到的结果没有意义。目前, 质量控制大多基于阈值给定, 主观性强, 导致离群值存在的可能。带约束项广义变分同化把M-估计法的权重函数耦合到经典变分同化中, 即使误差分布假设不满足高斯分布或存在离群值, 其变分同化极小化迭代仍可得到有意义的结果。带约束项广义变分依赖M-估计法权重函数的选取。因权重函数在目标泛函极小化迭代过程中会更新, 能实现观测项对目标泛函贡献率的动态重估计。

在选择M-估计权重因子时应满足一些约束, 其中最重要的约束是在变量空间要求单个函数为凸函数, 以保证最终结果收敛到唯一解。

带约束项广义变分同化不仅依赖权重因子的选取, 还依赖M-估计法中广义尺度的选取(Black et al, 1998)。此处广义尺度值也称为“离群值”阈值。广义尺度定义为c:

式中: y^obs和y^sim分别表示通道的观测和模拟亮温; MAD表示绝对偏差中值; median为中位数。

采用辐射传输模式RTTOV模拟AIRS通道亮温。RTTOV是快速回归模型在数据同化和反演系统中用于处理星下点红外和微波被动遥感。RTTOV垂直层从0.1~1 013.5 hPa共分为43层, 除地表温度、地表湿度、地表类型和温度、湿度、云中液态水等廓线外, 还包括云参数。具有正演模式、切线性模式、伴随模式以及雅克比(Jacobi)模式。亮温模拟使用RTTOV8.5扩展版本, 该版本包含云辐射传输模块(RTTOVCLD), 此时云被假定为单层灰体且不考虑光学厚度(Saunders et al, 2005)。

由辐射传输理论可知, 卫星探测器通道k观测到的辐射值可写成:

式中: N_e为有效云量; R_k^clear为晴空辐射; R_k^cld(p_{c, k})为云顶气压在p_{c, k}的有云辐射。

对于不透明云, N_e是有效云量; 对于半透明云, N_e则代表云的覆盖和云的发射率, 文中只考虑发射和吸收, 且假定不同通道N_e值不同, 半透明云的发射率独立于波长。给定大气温度、湿度和地表信息, 使用辐射传输模式进行模拟晴空和给定气压层范围的有云辐射。

在采用带约束项的广义变分同化前, 所有视场点FOV被视为潜在的有云视场点, 在执行变分同化极小化迭代时, 需输入初始背景场以得到云参数。

ECMWF廓线集中包含了云量的信息, 原始廓线共60层。借鉴Pavelin et al(2008)方法, 对廓线进行插值和加高斯随机误差前, 定义“有效云量”和“有效云顶气压”。廓线分层中云量最大值为“有效云量”。从顶层(0.1 hPa)向近地面层(998.82 hPa)算起, 当某层云量超过廓线有效云量10 %时, 对应气压层为该条廓线的“有效云顶气压”。该方法能有效处理云的分层结构。将得到的有效云量和有效云顶气压提供给辐射传输模式使用, 并将云参数和其他变量引入带约束项的同化极小化迭代过程中。带约束项广义变分同化有足够的自由度将云参数调整到实际值且拒绝不收敛视场点FOV。

从高云、中云和低云的有效云量和有效云顶气压示意图(资料选自ECMWF廓线集)(图 1)中可以看出, 高云视场点云量廓线较复杂, 最大云量为0.898, 对应的气压为787.23 hPa。由前面对于云量的定义可知, 有效云量定义为0.898;由顶层(0.1 hPa)向近地面层(998.82 hPa)算起, 第一次超过有效云量10 %的云量为0.35, 此时对应气压值为228.27 hPa, 即有效云顶气压。以此类推, 可以得到中云和低云的有效云量和有效云顶气压值。

图 1(Figure 1)

图 1 高云、中云和低云的有效云量和云顶气压定义解释 Figure 1 Illustrations of the definition of cloud-top pressure and effective cloud fraction for high, mid and low clouds

使用ECMWF数值预报模式的40年再分析资料共13 495条廓线集进行AIRS有云亮温带约束项广义变分同化研究(Chevallier, 2001)。该廓线集被设计为不同样本, 分为不同大气状态、季节和地形。廓线包括温度、湿度、臭氧、云中液态水含量、云冰含量、云量、地表温度、地表湿度和地表类型。基于气压对数坐标进行线性插值, 将原始13 495条廓线的60层信息插值到RTTOV的43层。插值后的廓线集为变分同化的背景廓线, 并将由上述方法得到的云参数提供给辐射传输模式RTTOV进行AIRS通道有云或晴空亮温模拟。

广义变分同化需进行大气状态的背景估计, 在实际变分同化系统中, 背景场是通过早期的模式预报得到。温度和比湿误差值的给定分别参考Eyre et al(1993)和Rabier et al(1998)的方法。因陆地地表类型和地表发射率的不确定, 选用ECMWF廓线集中海洋视场点作为“真实廓线”开展AIRS理想数据研究。真实温度和湿度廓线通过添加高斯随机误差作为背景廓线(Pavelin et al, 2008)。

式中: x_b为扰动后的背景廓线; x_t为真实廓线; L_i和λ_i分别为背景误差协方差矩阵B的特征向量和特征值; ε_i为均值为0方差为1的高斯随机数。

使用辐射传输模式RTTOV进行AIRS受云影响亮温模拟, 包含混合相位和多层云的影响。对选取的每个ECMWF数据集中的大气廓线, 模拟324个通道亮温。324个通道说明和分布见网址http://cimss.ssec.wisc.edu/~paulv/AIRS/AIRS_324_Channel_Subset.html。“真实廓线”相应的模拟亮温作为“真实亮温”。在AIRS的模拟亮温基础上加入随机测量误差, 记为“观测值”y:

式中: y为扰动后的观测亮温; H(x_t)为真实廓线模拟得到的亮温; ε₀为高斯随机误差; R为观测误差协方差矩阵。

涉及快速辐射传输模式RTTOV(Saunders et al, 2005)和变分同化系统2个模型。其中变分同化系统是基于欧洲卫星数值天气预报应用研究小组开发的一维变分同化模型基础上进行修改(杜华栋等, 2008)。具体修改如下: (1)关闭其他变量, 只考虑云参数、温度、湿度和O₃(气候态)等廓线及2 m温度、2 m湿度、u、v风场、地表气压和地表温度。(2)添加带约束项的广义变分同化和相应M-估计法代码。(3)改写目标泛函和相应泛函梯度的程序接口。

设定带约束项广义变分同化最大迭代次数为10步。在实际分析过程中大多数视场点FOV在第3或第4步迭代时收敛。具体算法步骤如下: (1)输入背景场, 采用RTTOV模拟卫星通道亮温。(2)统计各通道亮温偏差(O-B), 由公式(5)计算每个通道离群值阈值。(3)计算带约束项广义变分同化目标泛函和梯度, 包括背景、观测和云代价函数及梯度3个部分。该步需参考表 1选取相应M-估计法权重函数。(4)调用极小化算法进行求解(Limited-Memory BFGS, L-BFGS), 得到控制变量x_k=x_b+δx_k, k表示迭代次数。(5)若‖x_k-x_k-1‖₂ < ε或k > n_max, 则停止迭代, 否则转回(3)进行下一次迭代。ε是迭代停止阈值; n_max是极小化迭代中最大迭代次数。

在同化时对低云(800 hPa以下)、中云(800~500 hPa)和高云(500 hPa以上)分别考虑。

AIRS光谱区域覆盖了长波CO₂、短波CO₂、H₂O、O₃、CH₄和N₂O等吸收带。不同吸收带通道对云的敏感性各不同, 且AIRS各个通道的权重函数峰值层也不同。若只同化晴空视场点的通道亮温, 且不考虑地表发射率等对窗区通道的影响, 则不需进行通道选择。为更有效地同化受云影响的通道亮温, 在同化过程中有必要进行通道的动态选择。具体步骤如下: (1)计算AIRS通道子集中各通道在每个视场点的温度雅可比矩阵。温度雅可比矩阵表示模拟或观测亮温对模式空间温度变量的敏感程度(王根, 2014)。(2)在同化过程中, 当通道的温度雅可比垂直积分超过10 %位于云顶之下, 则拒绝视场点FOV中的通道亮温(Pavelin et al, 2008)。

图 2分别给出了高、中、低云的温度雅可比矩阵和温度廓线分布。该个例中, 高云、中云和低云分别选用324个通道中的100个、130个和175个通道。

图 2(Figure 2)

图 2 1Dvar中324通道的温度雅可比(左), 随机选取高、中和低云个例峰值在云顶之上的通道分布(中), 温度廓线(右) Figure 2 Temperature Jacobian for all channels used in the 1DVar (left) and selected channels peaking above cloud (middle) for randomly selected cases with high, mid and low clouds, and temperature profile (right)

为验证带约束项广义变分同化红外探测器受云影响亮温的可行性, 进一步进行温度和湿度同化反演。采用M-估计法(L2-估计、Huber-估计、Fair-估计和Cauchy-估计)的权重函数进行观测项对目标泛函的“动态”重估计。因湿度时空变化大且与亮温具有较强的非线性, 微小的湿度扰动, 都会引起亮温大的变化。故同化前对AIRS资料进行了初步质量控制, 当通道亮温偏差超过20 K时(Lsaksen, 2010), 剔除此通道亮温。形成“通道子集”后, 执行广义变分同化时极小化迭代时, 当通道亮温贡献率小于25 %时, 该通道亮温在该次迭代过程中不使用。

同化AIRS 324个通道亮温组合, 随机选取ECMWF廓线集中100个海洋视场点进行统计研究。通道的“动态选择”基于温度雅克比进行。图 3为采用不同方法估计的温度和湿度误差均方根统计。

图 3(Figure 3)

图 3 不同方法下高云、中云和低云的温度、湿度误差均方根统计 Figure 3 Root-mean-square errors (RMSE) statistics of temperature and humidity of high cloud, mid-level and low-level in different methods

图 4进一步给出了L2-估计与Huber-估计、Cauchy-估计和Fair-估计变分同化得到的温度和湿度偏差均方根误差统计。

图 4(Figure 4)

图 4 L2估计与不同方法估计的温度、湿度偏差均方根误差统计 Figure 4 Temperature and humidity RMSE statistics difference between L2 estimator and other different methods

结合图 3和图 4可以看出, 带约束项广义变分同化高、中、低云时温度、湿度同化反演的效果整体优于经典变分法。对于高云, 基于Cauchy-估计的变分同化反演结果最好, 而对于中云和低云, 基于Huber-估计得到了较好的同化反演结果。然而, 在反演模式高层温度时, 基于Fair-估计反而得到了较差的结果, 但其对于湿度反演效果较为理想, 其结果可能与Fair-估计分布的固有特点有关。说明采用广义变分同化部分有云亮温时需根据不同问题, 选取相应权重因子。

广义变分同化考虑了误差非高斯性, 在变分同化的过程中耦合了质量控制, 每次被“拒绝”的通道亮温能重新使用。“约束项”能实现云参数从“数学解”到“物理解”的过渡, 使结果具有“合理性”。本文分析和试验表明Huber-估计效果较好, 初步验证了Huber-估计是众所周知的“健壮统计”(Huber, 1981)。考虑到误差服从高斯分布只是一种理论上的假设, 受云影响的亮温偏差往往具有较强的非高斯, 因此研究误差服从非高斯变分同化具有前瞻性和必要性。

针对目前国际上高光谱大气红外探测器AIRS受云影响同化研究热点, 基于非高斯模型开展AIRS受云影响亮温同化初步研究, 将云参数(有效云量和有效云顶气压)用于变分极小化迭代和卫星通道亮温模拟, 得到以下主要结论:

(1) 区别于经典变分只能同化晴空视场点或晴空通道亮温资料, 带约束项广义变分同化能直接同化受云影响亮温, 使用资料信息量较多。

(2) 动态形成“通道子集”, 充分利用不同光谱带通道的亮温信息, 使温度雅克比峰值在云顶之上的通道亮温能较好地被同化使用。

(3) 区别于经典变分同化误差服从高斯分布, 广义变分同化耦合了M-估计法, 具有非高斯性, 带约束项能合理得到云参数信息, 从而更好地用于变分极小化迭代。

(4) 不同M-估计权重因子对变分同化结果影响较大, 耦合Huber-估计权重因子的广义变分同化比经典变分同化更能得到较好的结果。

本文只进行了AIRS模拟亮温的扰动理想试验。在后续科研工作中, 将根据不同数据特性给定不同权重因子, 并将带约束项广义变分同化思想用于WRFDA(Weather Research and Forecasting Model Data Assimilation)和GRAPES(Globe and Regional Assimilation and Prediction System)以开展大量实际个例的应用研究。