湿度分析效果对全球数值天气预报系统而言十分重要，它不仅直接影响湿度初值场的三维模拟，还会对中高纬度1—2 d的形势场预报效果产生影响，并对赤道地区2—3 d的风场预报影响显著(Andersson et al，2005，2007; Bengtssonet et al，2005)。然而，目前的湿度分析仍存在较多问题，如英国气象局全球数值天气预报系统的湿度分析在对流层中高层存在偏干现象(Inglebyet et al，2013)，欧洲中期数值预报中心(ECMWF)全球数值天气预报系统也出现过湿度分析偏干的现象(Tompkins et al，2007)。中国自主发展的GRAPES全球数值天气预报系统(水平分辨率0.25°，垂直60层)的湿度分析则存在明显偏湿现象，这在与中国、欧洲、北美等观测密集地区的探空观测资料以及ECMWF全球再分析资料(ERA-Interim)的对比中均表现明显。进一步分析表明，偏湿主要发生在模式背景相对湿度大于60%的地区，并且当相对湿度超过80%时最为严重。这也导致了GRAPES全球模式对中雨以下量级降水的空报增多。

造成GRAPES湿度分析效果不好的原因多样。一方面是模式预报本身具有湿偏差并反馈到湿度分析中，模式预报湿偏差的具体原因还在诊断中。另一方面与湿度观测存在系统性偏差以及同化所使用的湿度观测数量偏少有关。探空湿度观测普遍存在偏差，需要进行偏差订正(Christian，2006; Sun et al，2010; Bian et al，2011; 郝民等，2015)。ECMWF等业务数值天气预报中心对探空湿度资料进行偏差订正后，明显改善了湿度分析质量(Lorenc et al，1996; Agusti-Panareda et al，2009)。在ERA-Interim再分析中，还采用黑名单形式剔除质量不高探空的湿度观测(Dee et al，2011)。目前GRAPES全球三维变分同化(GRAPES-3DVar)所使用的探空湿度观测还没有进行偏差订正，也没有黑名单，影响了分析效果。此外，相较于ECMWF等先进业务中心综合使用常规、非常规多种湿度观测(Andersson et al，2007)，GRAPES-3DVar采用的湿度观测类型尚有欠缺，仅限于探空、地面、船舶等常规湿度观测。

而从变分同化系统的设计和发展角度来看，湿度分析自然还受到同化框架中湿度分析方案的影响，这也是本研究的着眼点。从该角度考察，湿度分析涉及的关键科学问题是如何构造合理的表征湿度信息的控制变量。该变量的构造又涉及两方面的问题:如何选择湿度分析变量以及如何构建湿度分析变量和其他动力学分析变量间的动力平衡约束，以便抽取出平衡部分，得到误差独立的控制变量。

为了计算的可行性，业务运行的变分同化框架均假设分析变量的背景误差满足高斯无偏特征。对于动力学变量而言，其误差特征能基本满足要求，而对于常用的湿度变量(例如比湿(q)、相对湿度等)，它们的误差不仅存在非高斯分布特征，而且系统性偏差明显。Dee等(2003)分别比较了相对湿度、比湿q、lgq、拟相对湿度作为分析变量的优劣，指出相对湿度和拟相对湿度的统计特征更为合适。基于此，美国GSI(Gridpoint Statistical Interpolation)系统以及研究模式WRF(Weather Research Forecasting)变分同化系统的湿度分析变量均采用了拟相对湿度(Kleist et al，2009; Huang et al，2009)。

对于湿度和动力学变量间的动力平衡约束，由于之间物理关系十分复杂，难以通过方程显式给出，早期发展的同化系统均未予以考虑，即直接将所选择的湿度分析变量作为控制变量。之后，在区域变分同化系统的发展中有采用统计方法直接构建湿度和动力学变量间的平衡约束的研究，但统计过程中未考虑实际的物理意义以及具体天气形势(Berre，2000; Chen et al，2013)。Hólm等(2002)分析指出，在考虑背景湿度信息的前提下，引入温度对湿度的统计解释来定义湿度控制变量，并采用对称变换处理，能明显改善湿度控制增量的误差统计特征，改进湿度分析质量，进而改善卫星湿度观测资料同化应用效果。这一成果也在ERA-Interim中使用(Dee et al，2011)。温度与湿度的物理约束关系的实质是：当相对湿度接近饱和时，云中水汽凝结加热带来的温度升高会减缓相对湿度进一步升高，而水汽蒸发吸收热量带来的温度降低会减缓相对湿度降低。若不考虑凝结加热影响，在接近饱和的相对湿度上继续叠加正的相对湿度分析增量可能出现超饱和现象。考虑温度与湿度统计约束关系后，对全球湿度分析质量的改善有正效果(Andersson et al，2005; Ingleby et al，2013)，对区域湿度分析也有一定作用(Gustafsson et al，2011)。英国气象局的研究表明，考虑温度与湿度约束关系后，对流层中高层湿度分析偏干的情况有所减缓(Ingleby et al，2013)。

为改进GRAPES-3DVar的湿度分析，本研究借鉴Hólm等(2002)的思想，采用偏差更小、更为接近高斯分布的湿度分析变量，同时利用温度与湿度的统计约束关系定义新的湿度控制变量，减缓GRAPES-3DVar中湿度分析偏湿现象具有重要的现实意义。另外，之前的研究主要集中于新的湿度控制变量构造的方案上，而对构造过程中温湿统计平衡约束关系的具体特征，特别是随纬度、高度变化的空间特征研究不足，在此将结合GRAPES-3DVar自身的框架特点对这些问题做进一步分析，完善该领域的科学认识。

根据贝叶斯后验概率理论，增量形式的变分同化目标函数可以表示为(Lorenc，1986; Courtier et al，1994)

点击浏览公式

式中，δx是增量形式的分析场(分析场减去背景场)，d表示观测增量(观测场减去背景场在观测空间的相当量)，H表示线性化的观测算子，而B和R分别表示背景场和观测场的误差协方差矩阵。在式(1)的推导过程中，一个基本假设就是背景和观测的误差分布满足高斯无偏特征。基于这样的假设，极小化式(1)得到的δx也满足高斯无偏分布。

在极小化求解式(1)给出的目标函数时，一个关键难点在于背景误差协方差矩阵B的表达和求逆运算。这是因为B矩阵是超大矩阵，并且接近病态。为解决该问题，业务运行的变分同化框架均采用变分预处理技术对B矩阵进行简化。变分预处理的第一步就是将分析变量δx中不同变量间的平衡约束抽取出来，构造相互独立的控制变量(记作δx_u)

点击浏览公式

式中，U_p常被称作物理变换算子。经式(2)变换得到的控制变量δx_u之间相互独立，它们对应的背景误差协方差矩阵退化为块对角阵，变分问题得以简化。

在GRAPES-3DVar中，采用标量场流函数ψ和势函数χ表征水平风场，采用无量纲气压π表征大气质量场，比湿q表征湿度信息，因而其增量形式的分析变量δx=[δψ，δχ，δπ，δq]^T(薛纪善等，2008)。在控制变量δx_u的构造中，主要考虑了风压场之间的平衡约束关系，将δχ和δπ中与δψ相平衡的部分抽取出来，剩余的δχ_u和δπ_u作为描述非平衡散度风和质量场信息的控制变量(王瑞春等，2014，2015)。在这一过程中忽略了湿度和动力学变量的动力平衡约束关系，也即直接将δq全量作为湿度控制变量。

在上述框架中，湿度变量的选取和构造存在两个主要问题：(1)引言中提到的分析变量q的误差特征并不满足高斯无偏分布特征；(2)忽略湿度变量和动力学变量间的平衡约束关系，与实际情形差异较大。为此，文中同样分两步出发构造新的湿度的控制变量：(1)构造满足高斯无偏分布的湿度分析变量；(2)提取湿度分析变量和动力学变量间的平衡约束。下面先讨论这两步的实施方案，获取新的湿度控制变量，然后再给出对应的湿度变换方案。

高斯无偏分析变量的具体含义指：分析变量在所有格点上的背景误差的总样本满足高斯无偏特征(Hólm et al，2002)。由于背景误差样本无法获取，一般采用同一时刻不同预报场的差值近似背景误差样本。这里不同的预报场既可以是从不同时刻起报，预报到同一时刻的预报场(NMC方法)；也可以是同一时刻起报的不同集合预报成员(集合方法)(王金成等，2014)。这里选用NMC方法获取变量背景误差的替代样本，有关样本的具体获取方案参见第3节参数统计部分。

对于常用的湿度变量，例如比湿和相对湿度，Dee等(2003)研究表明相对湿度增量(δf)的统计特征更满足需求。即便如此，相对湿度增量的概率密度函数(PDF)接近指数函数，并存在明显的偏差，仍需进一步处理(Hólm et al，2002)。进一步研究表明，如果依据相对湿度的背景场(f_b)的大小将δf进行分组，每一组的δf的统计特征更加接近高斯分布。也即相对于总的p(δf)，条件概率函数p(δf|f_b)更接近高斯分布。而上述两个概率密度函数的关系可以表述为

点击浏览公式

式中，φ是所选择的参考变量，这里为f_b。之所以出现p(δf|f_b)满足高斯分布，而积分得到的p(δf)接近指数函数，主要原因在于δf的标准差和偏差是随着参考变量f_b变化的。如果能找到一个参考量φ，使得δf的方差和偏差不随φ变化，那么δf整体的概率密度函数也是高斯型的。

为实现上述目标，可以对δf进行如下的标准化操作：δ=[δf-b(φ)]/σ(φ)，这里b(φ)和σ(φ)指针对参考变量φ的不同数值δf的偏差和标准差。但这会带来一个新的问题，由于b(φ)≠0，在没有观测资料进入同化系统时δf也不为0。为此，更为合适的方案是寻找一个参考量φ，使得b(φ)始终等于或接近于0，并将标准化操作简化为δ=δf/σ(φ)。

δf出现明显偏差的主要原因在于自然界中湿度变量是有界的，当f_b接近0时，δf有明显正偏差；而当f_b接近饱和时，δf有明显负偏差。上述现象可以从图 1中看出，对于f_b最小的一组中，δf样本出现正值“长尾”型分布(图 1a)；对于f_b最大的一组中，δf样本出现负值“长尾”型分布。为了克服这样的问题，将式(3)中的参考变量选为背景场和分析场的均值，也即φ==(f_b+f_a)/2，此时p(δf|f)=p(f_b-f_a|(f_b+f_a)/2)。考虑到f_a和f_b的概率密度函数接近一致，两者位置可以互换，置换后可得如下等式

图 1 在fb最小(a)和最大组(b)中，δf样本的概率分布特征(图中虚线为标准高斯分布) Fig. 1 Probability distributions of δf for groups with the lowest (a) and the highest (b) values of fb respectively (dashed lines denote the corresponding standard normal distribution)

图选项

点击浏览公式

根据式(4)，依据的值将δf进行分组，理论上可以保证样本的无偏性。图 2给出了上述操作得到的结果，可以看到无论是在f的最大组还是最小组，δf的偏差均远小于图 1，并且分布特征也更加接近高斯形态。

图 2 在最小(a)和最大组(b)中，δf样本的概率分布特征(图中虚线为标准高斯分布) Fig. 2 Probability distributions of δf for groups with the lowest (a) and the highest (b) values of respectively (dashed lines denote the corresponding standard normal distribution)

图选项

综上，对于δf而言，可以选择作为参考变量对其进行分组，并做标准化处理，也即δ=δf/σ()。这样的处理方案不仅适用于相对湿度全量，也可以应用于扣除平衡部分后的湿度控制变量。下面给出如何将湿度和动力学变量间的平衡约束考虑进来。

根据定义，相对湿度f等于一定温度T和压强p下水汽压e和饱和水汽压e_s(T)的比值，近似等于比湿q和饱和比湿q_s(T)的比值，也即f=q/q_s(T)。进一步的，饱和比湿q_s(T)≈εe_s(T)/p，其中ε(=0.622)是一个常数。那么，基于背景场，相对湿度的增量可以表示为

点击浏览公式

式(5)右端第1项就是Dee等(2003)给出的拟相对湿度的增量；第2项是气压变化造成的相对湿度变化，该项数值较小，可忽略；第3项是饱和水汽压变化造成的相对湿度变化。引入克拉伯龙-克劳修斯方程，将饱和水汽增量转换为温度增量，式(5)可以变换为

点击浏览公式

式中，L是混合相态下的凝结加热率，R_v是湿空气常数。现在，式(6)右端第2项是由于温度变化产生的相对湿度变化，在背景湿度较大、温度较低时该项与第1项可比。Hólm等(2002)的研究表明这两项在f_b接近饱和时存在密切联系，可以用来引入湿度与温度的平衡约束关系。

为简洁起见，定义变量q′、T′

点击浏览公式

式中，系数是背景温度和相对湿度的函数。进一步地，将q′分解为与T′相平衡的部分和非平衡部分(或称独立部分)

点击浏览公式

式中，Q_T是描述湿度和温度间平衡约束关系的算子。通过增量样本，采用线性回归的方法来求解算子Q_T，其计算公式为

点击浏览公式

式中，〈·，·〉是协方差求解算子。Q_T的计算和具体特征将在第3节给出。

根据式(8)，可以定义一个相对于温度独立的湿度变量q′u，并在此基础上通过2.2节给出的对称变换对q′_u做进一步处理，以保证其满足高斯无偏特征。最终，新的湿度控制变量可以定义为

点击浏览公式

式中，σ_u()是基于的大小分组求得的q′_u的标准差。

由式(10)获取了一个与温度分析相独立，并且满足高斯无偏分布的湿度控制变量。也即从与温度明显相关的δf转换到与温度变化相独立的q′_u，再通过高斯无偏处理获取。这样的变换过程在同化系统中并不显式出现，框架中实际出现的是上述变换的逆变换。同化框架以作为湿度控制变量进行极小化迭代，而在迭代过程中为了与观测信息相比较，需要将转换到δf、δq等常用湿度变量的增量。下面就给出这部分变换的实施方案。

为了导出常用湿度变量，首先将<转换为q′，根据式(8)和(10)，该变换可以表达为

点击浏览公式

由于σ_u和Q_T均是基于f的大小分组统计得到，因而都是f的函数，这里f=f_b+12δf 。根据式(11)计算得到q′后，再根据式(6)和(7)导出δf。此时，可以发现由向δf的转换是一个非线性过程，因为计算过程中就需要用到δf的值。在同化框架内循环的线性化迭代过程中，此处与ECMWF以及英国气象局一样，采用f_b替代式(11)中的f，也即采用线性计算近似非线性转换。

在内循环极小化迭代求得的收敛值之后，在外循环更新时采用非线性变换方案重新求解q′。重组式(11)，将其表达为如下非线性代数方程

点击浏览公式

在上式计算中，如果Q_T也随δf变化，表达式将会非常复杂，也会影响到非线性代数方程的单调性，因在实际计算中，Q_T仍取为随f_b变化的函数。进一步将式(12)表达为有关δf的非线性方程

点击浏览公式

由于σ_u 是由样本拟合出的函数，由式(13)直接计算函数F(δf)的导数不易，可采用割线法求解该非线性方程，一般迭代几步就可以获得收敛值δf^*。在计算得到δf^*之后，可以根据式(7)中q′的定义导出相应的δq^*。最后，δq^*的值被用于更新q的背景场值，并提供给模式作为湿度初值。

至此，完整给出了考虑温湿平衡约束关系的新湿度控制变量设计的科学思路以及基于该变量的湿度变换方案。在上述变换过程中，参数σ_u和Q_T是事先统计好的，下面给出有关这两个参数的统计方案和空间分布特征。

使用NMC方法获取误差样本来统计σ_u 、Q_T，两个预报时效分别为24与48 h。选取GRAPES全球模式1°分辨率24与48 h的气压、温度、比湿预报样本p²⁴、T²⁴、q²⁴与p⁴⁸、T⁴⁸、q⁴⁸，并把24与48 h的样本分别视为背景场与分析场，计算温度与比湿的样本平均(T、q)与样本差异(δT、δq)。样本时段为2015年5月31日—6月30日，共计获取61个误差样本。考虑到实际大气存在冰水混合相态，参照ECMWF的方法来计算冰水混合相态下的饱和水汽压，并由气压、比湿和饱和水汽压计算相对湿度f²⁴与f⁴⁸，求得样本平均f与样本差异δf。饱和比湿q_s与温度T是强非线性关系，不直接使用平均温度计算平均饱和比湿，而是先分别计算24与48 h的饱和比湿，再做平均得到q_s(T)。

对所有样本按f大小进行排序，并将样本按[f，f+Δf]分成多个组，统计σ_u(f)、Q_T(f)等参数随f的变化。考虑到不同组对应的样本数差异较大，为避免f很小和接近饱和情况下样本较少对统计结果造成的干扰，结合等Δf间隔和等样本数分组两个方案来确定Δf间距。具体实施中，将f(范围为[0.05，1.05])划分为51个样本间隔来归类样本，使得每一组中有足够多的样本数用于统计，以保证统计结果的稳健性。

分组统计得到离散形式的σ_u、Q_T不能直接用于计算，需要给出它们的连续形式。对离散σ_u、Q_T值沿着f方向进行多项式拟合，并考虑每组σ_u、Q_T数值的误差方差作为多项式拟合的权重，可以获得随f变化的连续σ_u、Q_T拟合函数。试验表明采用5阶多项式函数就能取得较高的精度。温度对湿度的统计平衡关系Q_T在f比较大时有物理意义，在f小于80%时出现无物理意义的负相关，需要对Q_T进行非负值限制，仅给出Q_T大于0的值。

下面分别给出统计参数的全球平均特征以及随纬度和高度变化的情况。并且为了进一步理解2.2节选用f替代f_b作为参考变量构造无偏的湿度控制变量的影响，还沿f_b做了相应统计计算，以便与沿f统计得到的结果作比较。

分别考察沿f_b与f分布的拟相对湿度q′的均方差σ_q 、扣除与温度平衡部分后的非平衡拟相对湿度q′_u的均方差σ_u以及温度对湿度解释Q_T参数的全球平均特征。如图 3所示，沿f_b分布的q′样本在不同f_b值下的偏差明显，σ_q在f_b值较小时相应较小，当f_b增大时稳定在18%左右；而沿f分布的q′样本偏差基本为0，σ_q 在f值居中时最大，在f值较大或较小时明显减小。根据式(6)和(11)，对于同一湿度控制变量，σ_q愈大，对应的δf也愈大。根据σ_q 沿f统计得到的特征(图 3b)，当f减小时，σ_q也减小，负δf增量幅度受到抑制，分析出现负水汽的现象得到抑制；反之在f接近饱和时，σ_q同样减小，正δf增量幅度受到抑制，分析出现过饱和水汽的现象得到抑制。传统的水汽分析变量容易出现负水汽或超饱和分析值，需要增加额外的负水汽和超饱和分析值的限制。而通过对不同f给定不同误差σ_q，隐含的对δf的分析增加了上下边界条件的约束，从而不再需要额外的约束，这是该方法的一个优点。

图 3 分别采用fb(a)和f(b)作为参考场时，参数σq、σu、QT的数值以及q′偏差的分布(σq,σu以及q′偏差的数值由左侧坐标轴表示，QT的数值由右侧坐标轴表示;单位：1) Fig. 3 σq ,σu ,the bias of q′(specified by left y-axes，units:1) and QT (specified by right y-axes,units:1) as a function of fb (a) or f (b)

图选项

Q_T表示温度与湿度间的统计平衡约束。具体地，根据式(7)和(8)，Q_T反映了δq和δT间的约束关系。根据Q_T 统计结果，在湿度接近饱和时，δq和δT间存在明显正相关(Q_T为正值)。其内在物理机制为：当f接近饱和时，继续增加水汽，水汽凝结加热将会带来温度升高；而减少水汽时，云水蒸发吸热会使温度降低。这样的相关关系主要体现在湿度接近饱和的情形下，而当f降至80%时，湿度变化和温度变化间的相关就已经很小，Q_T的数值接近为0。上述情形可以通过f的背景误差中能被温度解释的比例看出(图 4)：当f接近饱和时，温度对总相对湿度误差的解释可以达到80%，而当f减小到80%以下时，其对相对湿度误差的解释接近0，可以忽略不计。

图 4 600 hPa上，f的背景误差能被温度变化解释的比例 Fig. 4 The ratio of f variance explained by temperature at 600 hPa

图选项

统计结果表明，Q_T表征的温湿统计约束关系会随纬度产生明显变化，赤道地区与中高纬度地区的差异尤其明显。这与水汽分布特征和大尺度凝结发生区域有关。数值试验结果表明，不考虑这种差异会对赤道地区温度与湿度分析结果产生严重的影响，异常强的温湿平衡约束关系会导致温度与湿度分析产生系统性的负偏差。

为引入Q_T随纬度的变化，同时兼顾统计样本的数量，将全球分为8个纬度带，北半球被划分为0—15°，15°—30°，30°—60°，60°—90°N 4个纬度带，南半球按同样方式剖分。每个纬度带的统计结果赋值给该纬度带的中心纬度，其余纬度上的值由中心纬度上的数值线性插值得到。每个纬度带统计时，采用了相同数值的[f，f+Δf]分组方式。

图 5给出北半球4个纬度带统计得到的σ_u，南半球与之类似。从图中可以看到，赤道地区300 hPa附近出现的σ_u极大值中心随着纬度升高而逐渐降低，这与湿度相关物理过程发生的垂直位置一致。此外在对流层顶，也出现了另一个σ_u 大值区域，与该区域的变率较大有关。此外，σ_u 极大值位于f约为50%的区段，并随f的增大或减小而减小，这与图 3b的结论一致。

图 5 σu(等值线)随纬度带和高度的变化特征(填色是参考场f的值，横坐标是样本分组的组数) (a.0—15°N,b.15°—30°N,c.30°—60°N,d.60°—90°N) Fig. 5 σu (counter lines) as a function of latitude band and vertical pressure (color shaded: reference f, x-axes: the group number of the samples) (a.0—15°N,b.15°—30°N,c.30°—60°N,d.60°—90°N)

图选项

图 6给出不同纬度带上Q_T的统计结果。中高纬度400 hPa附近，Q_T较大，在f＜70%的情况下，温度变化与湿度变化仍然有统计相关，这种相关随着纬度减小而逐步减弱。赤道地区，中低层温度与湿度变化基本不相关，而在200 hPa以上高度，温度与湿度变化也有一定相关。图 7a进一步给出了经插值和纬圈平均得到的Q_T随纬度的连续变化特征，结论与图 6一致，Q_T大值区主要出现在中高纬度地区。

图 6 QT(等值线，单位：1)随纬度带和高度的变化特征(说明同图 5) Fig. 6 Same as Fig.5 but for QT (units:1)

图选项

图 7 QT(单位：1)(a)、温度对相对湿度背景误差均方差的贡献()(b)随纬度和高度变化的特征 Fig. 7 QT (units:1) (a) and the contribution of temperature to the standard deviation of relative humidity(, make σT =1 K) (b) as a function of latitude and vertical pressure

图选项

根据式(6)—(8)，新方案中f的背景误差均方差σ_f可由σ_u、Q_T的值计算得到

点击浏览公式

式中，σ_f 随f演变，而在线性变换中仅随f_b 演变。在式(14)中，包含Q_T 的项反映的是温湿平衡约束关系对σ_f 的贡献，Q_T 越大，σ_T 对σ_f的贡献也越大。为方便分析，在所有格点取σ_T =1 K，计算其对σ_f的贡献，结果见图 7b。从图中可以看出温湿平衡约束的贡献在对流层低层较为微弱，而在对流层中高层更为明显。这是因为α不但与背景相对湿度有关，也与背景温度有关，饱和情形下α的取值随温度降低而增大。同样的，由于背景温度随纬度升高而减小，因而温湿平衡约束对σ_f的贡献在高纬度要更大一点。ECMWF在评估线性化物理过程对12 h预报误差减小的影响时发现：大尺度水汽凝结物理过程对赤道外对流层中层、赤道地区对流层顶附近的贡献大(Mahfouf，1998)。其大尺度水汽凝结的影响区域与图 7所示的温湿平衡约束贡献较大的区域对应较好，说明Q_T主要代表了大尺度水汽凝结的作用。而赤道地区主要以对流性降水活动为主，需采用次网格对流参数化过程描述，从统计上不能给出温度对湿度的解释部分。这也是温湿统计平衡约束在赤道地区贡献较小的原因之一。

为了考察引入的新的湿度控制变量的具体效果，分析温湿统计平衡约束关系在其中起到的作用，进行了一系列理想试验。首先给出的是线性湿度变换的结果，也即在式(11)计算过程中采用f_b替代f，该变换在内循环中使用，然后给出非线性湿度变换情形下的结果，并对比两者的差异。

σ_u、Q_T随f_b演变的特点，决定了分析增量δf 也具有随流型演变的特征。为考察湿度分析的流依赖特征，采用单点理想观测试验予以检验。任意选取(20°—45°N，110°—80°W)作为研究区域，该区域是2014年5月一次冷锋锋面过境情况。如图 8a所示在95°W经线上，相对湿度从低于10%快速增大到100%，温度和气压也有较大的变化。在该区域，σ_u随f_b的不同而变化，呈现出随流型演变的特征(图 8b)。而在相对湿度较大的区域，由估计出的温度对σ_f贡献部分也较明显，量值与σ_u可比(图 8c)。

图 8 单点试验附近区域背景场与背景误差分布特征 (a. 背景场的气压(白色线，单位：hPa)、温度(黑色线，单位：K)、 相对湿度(填色，单位：%)在模式13层(近500 hPa)的分布情况， b.σu的分布， c. 温度对σf的贡献部分，等值线，单位：%) Fig. 8 Background states and error distribution around the single pseudo-observation (a. background pressure (white contour lines, units: hPa), temperature (black contour lines, units: K), and relative humidity (color shaded, units: %, also shown in plot (b) and (c)) at model level 13 (≈500 hPa); b. the distribution of σu; c. the contribution of temperature to σf ，contour lines, units: %)

图选项

在该区域500 hPa高度的锋面附近，选取3个单点(经度为95°W，纬度分别为25.5、31.5、35.5°N，记为单点1、2、3)，其背景相对湿度f_b分别为9%、53%、100%，给定新息向量f_o-f_b为15%、15%、-15%，逐点进行单点分析试验。在选择的3个单点位置，其对应背景场和背景误差特征如表 1所示。

根据表 1，对于单点1，f_b=9%，Q_T=0，且α的数值较小，因而温度对相对湿度背景误差σ_f的贡献很小，σ_f的值接近σ_u。由于该点所处区域的f_b相对均匀，σ_u分布也较为均匀，分析增量基本呈现出各向同性的圆形分布(图 9a)。对单点2，f_b=53%，对应的Q_T也为0。该点所处的区域f_b变化幅度大，从10%变化到100%，σ_f变化很大，分析增量不再呈现各向同性的圆形分布，而是呈现出一定的随流型演变的特点(图 9b)。对单点3，其f_b=100%，背景相对湿度饱和，温度分析增量对σ_f的贡献变大。f_b的快速变化使得f分析增量的极值中心偏离观测所在位置，位于f_b≈80%处(图 9c)。分析表 1的数据表明，500 hPa上σ_f的贡献主要来自湿度本身，而温度分析增量的贡献在高湿区也占一定比例。

图 9 采用线性湿度变换时，单点f观测产生的相对湿度分析增量(等值线，单位：%) (a. 单点1 (25.5°N，95°W)，b. 单点2(31.5°N，95°W)， c. 单点3 (35.5°N，95°W)；图中填色为背景相对湿度的值) Fig. 9 Analysis increments (contour lines, units: %) of relative humidity generated by a pseudo f observation at 95°W and at 25.5°N (a), 31.5°N (b), 35.5°N (c) when the linear humidity transform is used (the color shaded show the background f displayed in Fig. 8)

图选项

首先在单点上模拟非线性湿度变换的作用。取背景温度T_b为273 K，温度增量δT为0.5 K，将f_b从1%到100%等间距分为101个组。针对不同f_b值，令控制变量在-6到6(间距为1)范围内变化，分别采用线性和非线性变换两种方式计算δf，结果如图 10所示。对于线性变换，在f_b较大或较小时，相对湿度的正负分析增量呈现非对称分布，并且由于其对应的σ_u较小，分析增量也比较小，而当f_b居中时对应的分析增量最大。而非线性变换计算得到的δf沿f(垂直于f_a=f_b所在的45°斜线方向)成对称分布。

图 10 针对一理想温度观测，分别采用线性(a) 和非线性(b)湿度变换计算得到的相对湿度增量 以及两个方案计算值的差异 (c) Fig. 10 Different f analysis responses to a perfect temperature observation by using linear (a) or nonlinear (b) humidity transform, and their deviations are (c)

图选项

从非线性与线性变换计算得到的δf 差异(图 10c)中可以看到，当f较小时，对于正分析增量，因σ_u随f增大而变大(图 3b)，正分析增量在非线性调整后变大，越大非线性调整越大；对于负分析增量，因σ_u随f减小而变小，非线性调整作用很弱。当f较大接近饱和时，因σ_u随f增大而减小，负分析增量在非线性调整后产生更大的分析增量，反之则效果不明显。非线性调整带来的这些变化可以帮助限制负水汽和超饱和的出现。

非线性湿度变换的作用也可以从单点试验(单点设置同图 9)中看到，对单点1(图 11a)，正值δf 在非线性调整后进一步增大，与图 10c低相对湿度的结果一致。对单点2(图 11b)，在低湿度区，正值δf 在调整后数值变大，而在高湿度区，正值δf 的数值变小。在线性变换结果中(图 9b)，受水平误差相关模型的影响，单点2附近饱和湿度区出现了数值较小的正分析增量，使得湿度分析超饱和。而引入非线性调整后，负水汽和超饱和现象的出现受到限制，因而调整后的差异出现负值。对单点3(图 11c)，负值δf 在调整后数值进一步变小(绝对值变大)，与图 10c高相对湿度的结果一致。值得注意的是，单点试验结果也表明，引入非线性湿度变换对分析增量的影响并不是很大，后面的批量试验结果也显示这一现象。其原因可能与目前GRAPES变分框架只采用了一次外循环更新有关，Hólm等(2002)研究表明非线性变换的作用在采用多重外循环更新时较为明显。

图 11 引入非线性湿度变换对图 9单点试验中相对湿度分析增量δf的影响 (图中等值线为非线性变换得到δf与线性变换条件下δf的差值；其余设置同图9) Fig. 11 Same as Fig. 9, but the contour lines show the difference in f increments between using nonlinear humidity transform and using linear humidity transform

图选项

试验使用2015年集成的GRAPES全球同化预报系统，模式垂直分辨率60层，模式层顶高度为36.4 km，模式水平分辨率0.5°，三维变分内循环分辨率1.0°，外循环分辨率0.5°，同化系统每6 h做一次分析。使用的观测资料包括探空、小球测风、地面、船舶浮标、飞机报、云导风等常规资料，卫星GPS掩星资料，NOAA16/18/19及METOP-A 温度垂直探测仪资料AMSUA。这些观测中的湿度信息仅包括：探空湿度廓线、地面和船舶湿度观测，湿度观测数量明显偏少，且全球分布很不均匀。为方便对比，选取GRAPES连续试验常用的2013年5月进行不同湿度控制变量的分析与预报对比数值试验。将采用本研究发展的非平衡拟相对湿度作为湿度控制变量的连续试验(下文简称为“方案”)与原有采用比湿作为湿度控制变量的试验(下文称为“q方案”)进行比对。Bengtsson等(2005)研究表明，初值中水汽分析的作用主要影响1—3 d的预报，因此对比试验仅比较84 h内的预报差异。

如图 12a，在东亚区域约400 hPa高度，q方案得到的δf 普遍比较大，且出现超过-100%的异常值。这样的异常在其他区域也有所表现，主要原因是q方案中温度与湿度各自独立分析，且湿度分析没有考虑不能出现负水汽与超饱和水汽的物理限制。而在考虑了温湿统计平衡约束关系的′_u方案中，相对湿度分析增量比较均匀，分析增量幅度在-30%—30%，没有出现负水汽现象，且温度分析与湿度分析耦合较好(图 12b)。对比两种方案在图 12中的差异，在q方案δf 超过-100%的位置，方案出现明显的正温度分析增量，同时比湿q的负分析增量非常小(图略)。正温度分析增量结合非常小的负比湿分析增量，会限制δf出现异常负值，从而使得湿度分析更为合理。

图 12 q方案(a)和方案(b)给出的相对湿度(填色)和温度(等值线，单位：K)分析增量在近400 hPa模式层上的分布 Fig. 12 Analysis increments of relative humidity (color shaded) and temperature (contour lines, units: K) at model level (≈400 hPa) for the experiments with q scheme (a) and scheme (b)

图选项

首先采用探空湿度观测检验两个方案得到的湿度分析的差异。使用中国、欧洲、北美3个探空测站密集、探空仪型号相对统一的区域进行检验(图 13)。与探空观测相比，q方案得到的湿度分析存在明显的湿偏差，600 hPa以下偏湿1%—4%，500—200 hPa偏差达到10%—15%。偏湿现象主要发生在背景相对湿度f_b大于60%的情形下，并在f_b超过80%时表现最为显著。这一特征在3个检验区域较为一致。与之相对比，方案的湿度分析偏差在600 hPa以上有所减小，均方根误差也有所减小。而如果将探空观测按北半球、赤道、南半球予以划分，分别统计两个方案的分析结果与探空的差异，得到的结论也与图 13一致(图略)，表明方案略有优势。

图 13 q方案(exp1)和方案(exp2)的相对湿度分析 相对于中国(a)、欧洲(b)、北美(c)探空观测的偏差与均方根误差 Fig. 13 Biases and RMSEs of humidity analysis in the experiments with q scheme (exp1) and scheme (exp2) verified against radiosonde f observationsin China (a), Europe (b), and North American (c)

图选项

考虑到探空湿度观测本身可能存在明显偏差(Christian，2006; Sun et al，2010; Bian et al，2011; 郝民等，2015)，ECMWF等业务数值预报中心在使用探空湿度资料之前进行了偏差订正(Lorenc et al，1996; Agusti-Panareda et al，2009)，并且在ERA-Interim分析中，采用黑名单的形式剔除少数质量很差的观测(Dee et al，2011)，以提高湿度分析质量。而目前GRAPES-3DVar在湿度分析时，尚未对探空湿度观测进行偏差订正，因而需要进一步补充检验。郝民等(2015)的研究表明，在对流层中低层，ERA-Interim的分析是比较可信的。因而，这里选用ERA-Interim再分析资料对两个方案得到的湿度分析做进一步检验。检验发现，与q方案相比，′_u方案的相对湿度偏差在北半球中下层(700—500 hPa)以及赤道地区略有减小，而在南半球显著减小。同样的，相对湿度的均方根误差在北半球有减少，南半球显著减少，在赤道略有减少(图 14)。

图 14 q方案(exp1)和方案(exp2)的相对湿度分析相对于ERA-Interim再分析资料 在北半球(a)、赤道(b)和南半球(c)的均方根误差 Fig. 14 RMSEs of relative humidity analysis in the northern hemisphere (a), the tropics (b), and the southern hemisphere (c) for the experiments with q scheme (exp1) and scheme (exp2) (The verification is done against ERA-Interim reanalysis data)

图选项

两个方案湿度分析的差异也可以从两者背景误差幅度的不同给出一些解释。方案实际使用的相对湿度背景误差(由式(14)隐式给出)在南北半球中高纬度的大部分区域要大于q方案的值(图略)。当观测误差一致而背景误差更大时，观测的作用将更为明显，分析增量幅度将加大。

前文分析表明，与q方案相比，方案可以帮助减小湿度分析的偏差，这一特征在预报中也有所体现，并对降水产生影响。利用中国2400站地面自动站降水观测资料来检验GRAPES全球模式00时(世界时)预报的12—36 h、36—60 h、60—84 h三个时段内的降水预报。分别使用ETS评分、Bias评分给出降水预报技巧的变化，同时根据不同降水阈值的样本数给出检验结果的95%的统计置信(图 15)。

图 15 q方案(exp1)和′u方案(exp2)的36—60 h(a)和60—84 h(b)降水的ETS(a1、b1)和 Bias(a2、b2)评分(两个方案评分差异落在红色方框外代表通过95%的信度检验) Fig. 15 ETS (a1，b1) and Bias (a2，b2) scores of 36-60 h (a) and 60-84 h (b) precipitation forecasts in the experiments with q scheme (exp1) and scheme (exp2) (Differences outside the hallow bars are significant at the 95% confidence level)

图选项

对于12—36 h降水，两个方案的ETS评分与Bias评分均差异很小，均不显著(图略)。对于36—60 h降水，25 mm以下量级的ETS变化并不显著，不过方案的偏差明显减小，0.1—10 mm量级的偏差更接近于1；而对25 mm以上量级，Bias小于1，漏报明显，但两个方案的评分差异不显著(图 15a)。对于60—84 h降水，ETS评分的差异也不显著，方案在10 mm量级的Bias评分更接近1，表现显著，但25 mm以上量级的Bias评分明显小于1，漏报现象更为严重。降水的这种变化与方案带来的湿度分析的湿偏差明显减小相一致。

由于水平分辨率为0.5°的全球模式的小到中雨(0.1—10 mm)降水样本非常多，而大雨以上量级的样本较少，小到中雨的差异意义更为重要一些。在小到中雨的预报中，方案在不显著减小ETS评分的同时，显著改善了空报过大现象，表现出了一定优势。引入新的湿度控制变量对其他预报要素，如风、温度、位势高度等的影响较小。

改进资料同化系统中湿度分析质量对改善全球模式的湿度和降水预报非常重要。GRAPES-3DVar系统的湿度分析与探空观测以及ECMWF再分析相比都存在明显的湿偏差，改善其分析质量是本研究的主要目的。湿度变量的空间变化往往非常剧烈，且存在上下边界，选取合适的湿度控制变量对改善湿度分析质量十分重要。目前GRAPES-3DVar采用比湿作为湿度控制变量，但比湿的误差特征不满足无偏高斯分布假设，这是导致湿度分析存在湿偏差的一个重要原因。为此，借鉴Hólm等(2002)的思想，结合GRAPES-3DVar自身框架特点，采用考虑了温湿统计平衡约束关系的非平衡拟相对湿度作为新的湿度控制变量，并利用2个月的样本统计了相关参数。与Hólm等(2002)工作不同，本研究将全球分为8个纬度带分别统计，使得统计得到的平衡约束关系随纬度变化。在每个纬度带，采用多项式拟合办法来获得温湿统计平衡约束随相对湿度变化的连续函数，同时考虑了不同纬度带间参数的平滑过渡。

参数的统计结果表明，温湿统计平衡约束作用显著的区域主要出现在南北半球中高层相对湿度大于80%的区域，赤道地区作用较弱，具有较强的随纬度变化的特征。理想试验和一个月的批量试验结果表明：对称处理后的拟相对湿度的偏差基本为0，满足高斯分布特征；新方案中相对湿度背景误差会随背景场变化，这使得湿度分析具备流依赖特征，并有效地抑制了负水汽与超饱和水汽的出现；考虑温湿统计平衡约束关系的新方案使得湿度分析的偏差和均方根误差均有所减小。最后，基于中国2400站的降水检验表明，对0.1—10 mm降水预报，新方案在ETS评分没有显著降低的情况下，使得Bias评分更靠近1，空报程度有所减缓，体现了该方案的优势。不过，在60—84 h时段25 mm以上降水的检验中，新方案中漏报现象有所增大，表明未来还需从湿度的同化分析和物理过程两方面协同进行进一步的改进工作。

本研究是系统改进GRAPES全球同化预报系统湿度分析、预报能力工作的第一步，模式中与水物质相关的物理参数化过程的改进也在同步进行。在同化分析方面，除了从框架角度进一步完善文中引入的新的湿度变量变换方案以外，还需不断夯实和完善大量湿度观测资料的有效应用。例如，探空湿度观测普遍存在的干偏差会对湿度分析造成明显影响，对它进行偏差订正非常关键。此外，目前GRAPES-3DVar中湿度观测资料明显偏少，提高卫星湿度垂直探测仪资料的同化能力，是弥补探空、地面等常规观测不足的重要途径。现有资料的完善和新资料的不断引入过程，将为本方案的不断完善以及效果的体现提供重要机遇。

致 谢：感谢与ECMWF Hólm博士就本研究所做的深入探讨，感谢与中国气象局数值预报中心刘坤博士就GRAPES模式预报偏差产生原因的讨论，感谢数值预报中心赵滨博士对降水预报检验提供的帮助。