1 引言

提高降水的精确预报水平对数值预报是一个很大的挑战 (Sun, 2005)。除了模式本身的误差外，制约降水精确预报水平的一个重要因素是常规观测得到的初始场很难对中小尺度天气系统进行精细的描述 (Dance, 2004)。在众多观测资料中，多普勒雷达资料具有很高的时空分辨率，其观测资料的空间分辨率达到了1公里，时间分辨率为6分钟，能够对大气内部三维结构进行更为精细化的描述，这对于改善中小尺度天气的初始场具有重要意义 (王洪等，2015)。

多普勒雷达资料的有效应用依赖于同化技术，目前的数据同化方法有以三维变分同化 (three-dimensional variational analysis，简称3DVar) 和四维变分同化 (four-dimensional variational analysis，简称4DVar) 为代表的变分同化方法和以集合卡尔曼滤波 (ensemble Kalman filter，简称EnKF) 为代表的顺序同化方法。已有的研究表明，3DVar (Xiao et al., 2005, 2007; Hu et al., 2006a, 2006b; Xiao and Sun, 2007)、4DVar (Kawabata et al., 2011; Choi et al., 2013; Wang et al., 2013; Sun and Wang, 2013) 或者EnKF (Tong and Xue, 2005; Askoy et al., 2009, 2010; Dowell et al., 2011) 能够有效同化雷达资料，对提高降水的精细预报具有重要作用。相对于3DVar，4DVar在资料同化中更具优势，然而4DVar需要的切线性和伴随模式不仅计算量大，且其编码和升级也十分困难，特别是当数值模式存在高度非线性，涉及的模式物理过程具有不连续过程的情况下 (Xu, 1996a, 1996b; Zou, 1997)。这些都使得4DVar的业务应用具有更大的难度。此外，4DVar的背景误差协方差虽然在同化窗口内可以隐式发展，但其在同化时间窗口初端仍采用静态的背景误差协方差，这也制约了其同化性能的发挥 (Beck and Ehrendorfer, 2005; Cheng et al., 2010; Fairbairn et al., 2014)。而EnKF由于计算过程相对简单 (Evensen, 2003)，更易于实施，它在雷达资料同化中得到了越来越多的应用，许多试验 (Tong and Xue, 2005; Askoy et al., 2009, 2010; Dowell et al., 2011; 兰伟仁等, 2010a, 2010b) 也表明EnKF能够有效同化雷达资料。EnKF利用集合样本统计的背景误差协方差具有随流型的特点，能够对观测信息进行更为合理的传播，这对于变化比较快的降水天气的预报具有重要意义。

综合考虑4DVar和EnKF两种同化方法的特点，国内外许多学者提出了集合四维变分同化方法 (4DEnVar；Tian et al., 2008, 2011; Fairbairn et al., 2014)，4DEnVar结合了4DVar和EnKF两种同化方法的优势，规避了各自的不足。4DEnVar可以同化多个时刻的观测资料，但是却不需要4DVar所需要的切线性和伴随模式，因而使得同化过程变得易于实施；其次，4DEnVar采用集合样本统计的背景误差协方差具有随流型的性质，具有EnKF的优势。利用4DEnVar对各类资料的同化试验也表明了这种同化方法的潜力 (Qiu et al., 2007; Liu et al., 2008, 2009; Zhang et al., 2009; Wang et al., 2010; Zhang and Zhang, 2012; Tian et al., 2014)。在众多4DEnVar同化方法中，Tian et al.(2008, 2011, 2014) 基于本征正交分解，提出了集合四维变分同化方法 (POD4DEnVar)，并将其应用于各类资料的同化，试验结果均表明了这种方法的有效性和优越性。基于POD4DEnVar理论，Zhang et al.(2015)首次将其应用于雷达资料同化中，建立了雷达资料同化系统PRAS (PODEn4DVar-based Radar data Assimilation Scheme)，结果表明PRAS能够有效的改善初始场进而提高降水预报。然而POD4DEnVar假设模式扰动和观测扰动之间为线性关系，事实上数值模式和非常规资料的观测算子往往为高度非线性，无疑这个线性假设会直接影响到同化的效果 (Tian and Feng, 2015)。由于PRAS是基于POD4DEnVar建立的，因而不可避免受到线性化假设的影响。为了解决POD4DEnVar在高度非线性情况下的应用问题，Tian and Feng (2015)提出了利用高斯牛顿迭代格式逼近非线性4DVar分析解的NLS-4DVar (Non-Linear Least Squares-4DVar) 同化方法，作为POD4DEnVar方法的一种改进。

本文基于NLS-4DVar理论，对PRAS进行了改进，建立新的雷达资料同化系统NRAS (NLS-4DVar-based Radar data Assimilation Scheme)。为了考察NRAS的性能，以WRF (Weather Research and Forecasting) 模式作为预报模式，设计了观测系统模拟试验 (Observing System Simulation Experiments，简称OSSEs) 和实际雷达资料同化试验。

2 NRAS同化系统介绍 2.1 NLS-4DVar方法简介

由于NRAS是基于NLS-4DVar方法构建的，因而首先对NLS-4DVar方法作一简单介绍，详细的推导可参阅Tian and Feng (2015)。NLS-4DVar的出发点仍基于如下增量形式的4DVar代价函数：

$ \begin{array}{l} J(\mathit{\boldsymbol{x}}') = \frac{1}{2}{\left({\mathit{\boldsymbol{x}}'} \right)^{\rm{T}}}{\mathit{\boldsymbol{B}}^{ - 1}}\left({\mathit{\boldsymbol{x}}'} \right) + \frac{1}{2}\sum\limits_{k = 1}^s {[{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{M}}_{0 \to k}}({\mathit{\boldsymbol{x}}_{\rm{b}}} + \mathit{\boldsymbol{x'}}) - } \\ {\mathit{\boldsymbol{y}}_{{\rm{obs}}, k}}{]^{\rm{T}}}\mathit{\boldsymbol{R}}_k^{ - 1} \times \left[ {{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{M}}_{0 \to k}}({\mathit{\boldsymbol{x}}_b} + \mathit{\boldsymbol{x'}}) - {\mathit{\boldsymbol{y}}_{{\rm{obs}}, k}}} \right], \end{array} $

(1)

其中，x'表示待求的分析增量，T表示矩阵转置，x_b表示背景场，k表示同化窗口内第k个观测时刻，s表示同化时间窗口内总的观测时刻数，H_k是第k个观测时刻的观测算子，M_{0→ k}是预报模式，y_{obs, k}表示第k个时刻的观测值，B表示背景误差协方差，R_k表示和第k个时刻的观测误差协方差。

NLS-4DVar在极小化求解x'时，首先给出一组集合样本，假定集合样本数为N个，x_j(j=1, …N) 表示各集合成员，其相对于背景场x_b可以产生N个模式扰动$\mathit{\boldsymbol{x'}}_j^{} = {\mathit{\boldsymbol{x}}_j} - {\mathit{\boldsymbol{x}}_b} $(j=1, …N)，从而构成模式扰动空间 ($\mathit{\boldsymbol{x'}}_j^{}, j = 1, \cdots, N$)。假定分析增量x'线性蕴含在模式扰动空间中，即x'可以表示为如下模式扰动的线性组合形式：

$ \mathit{\boldsymbol{x'}} = {\mathit{\boldsymbol{P}}_x}\mathit{\boldsymbol{\beta }} $

(2)

其中，${\mathit{\boldsymbol{P}}_x} = (\mathit{\boldsymbol{x'}}_1^{}, \mathit{\boldsymbol{x'}}_2^{}, \cdots, \mathit{\boldsymbol{x'}}_N^{}) $表示模式扰动空间；$\boldsymbol{\beta} = {({\beta _1}, {\beta _2}, \cdots, {\beta _N})^T} $，表示模式扰动的线性组合系数。

利用模式扰动可以得到相应的集合背景误差协方差$\boldsymbol{B}={{\boldsymbol{P}_{x}}\boldsymbol{P}_{x}^{\text{T}}}/{(N-1)}\;$，将其和公式 (2) 代入到公式 (1) 得到如下形式的代价函数：

$ \begin{array}{l} J(\mathit{\boldsymbol{\beta }}) = \frac{1}{2}(N - 1){\mathit{\boldsymbol{\beta }}^{\rm T}}\mathit{\boldsymbol{\beta }} + \frac{1}{2}\sum\limits_{k = 1}^s {[{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{M}}_{0 \to k}}({\mathit{\boldsymbol{x}}_{\rm{b}}} + } {\mathit{\boldsymbol{P}}_x}\mathit{\boldsymbol{\beta }}) - {\mathit{\boldsymbol{y}}_{{\rm{obs}}, k}}{]^{\rm{T}}} \cdot \\ \mathit{\boldsymbol{R}}_k^{ - {\rm{1}}} \times \mathit{\boldsymbol{H}}\left[ {{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{M}}_{0 \to k}}({\mathit{\boldsymbol{x}}_{\rm{b}}} + {\mathit{\boldsymbol{P}}_x}\mathit{\boldsymbol{\beta }}) - {\mathit{\boldsymbol{y}}_{{\rm{obs}}, k}}} \right] \end{array} $

(3)

从公式 (3) 可见，新形式的代价函数将对x'的求解转化为对线性组合系数$ \boldsymbol{\beta} = {({\beta _1}, {\beta _2}, \cdots, {\beta _N})^T} $的求解，由于集合成员数目N一般远小于模式变量的维数，因而大大降低了计算量。POD4DEnVar对公式 (3) 直接最小化求解，得到线性组合系数β(Tian et al., 2008, 2011)。相对于POD4DEnVar，NLS-4DVar则进一步将公式 (3) 转化为如下的非线性最小二乘法形式：

从公式 (3) 可见，新形式的代价函数将对x'的求解转化为对线性组合系数$ \boldsymbol{\beta} = {({\beta _1}, {\beta _2}, \cdots, {\beta _N})^T} $的求解，由于集合成员数目N一般远小于模式变量的维数，因而大大降低了计算量。POD4DEnVar对公式 (3) 直接最小化求解，得到线性组合系数β(Tian et al., 2008, 2011)。相对于POD4DEnVar，NLS-4DVar则进一步将公式 (3) 转化为如下的非线性最小二乘法形式：

$ J(\mathit{\boldsymbol{\beta }}) = \frac{1}{2}Q{(\mathit{\boldsymbol{\beta }})^{\rm{T}}}Q(\mathit{\boldsymbol{\beta }}) $

(4)

其中，

$ {\rm{Q}}\left(\mathit{\boldsymbol{\beta }} \right){\rm{ = }}\left({\begin{array}{*{20}{c}} {\mathit{\boldsymbol{R}}_1^{ - 1/2}\left[ {{\rm{L'}}_1^{}({\mathit{\boldsymbol{P}}_{\rm{x}}}\mathit{\boldsymbol{\beta }}) - \mathit{\boldsymbol{y'}}_{{\mathop{\rm obs}\nolimits} {\rm{, }}1}^{}} \right]}\\ \vdots \\ {\mathit{\boldsymbol{R}}_s^{ - 1/2}\left[ {{\rm{L'}}_{\rm{s}}^{}({\mathit{\boldsymbol{P}}_{\rm{x}}}\mathit{\boldsymbol{\beta }}) - \mathit{\boldsymbol{y'}}_{{\mathop{\rm obs}\nolimits}, {\rm{s}}}^{}} \right]}\\ {\sqrt {{\rm{N}} - 1} \mathit{\boldsymbol{\beta }}} \end{array}} \right) $

(5)

$ {\mathit{\boldsymbol{R}}_{\rm{k}}}{\rm{ = }}\mathit{\boldsymbol{R}}_{\rm{k}}^{1/2}{(\mathit{\boldsymbol{R}}_{\rm{k}}^{1/2})^T} $

(6)

$ \begin{align} & L_{k}^{'}(\mathit{\boldsymbol{{x}'}})={{\mathit{\boldsymbol{H}}}_{k}}{{\mathit{\boldsymbol{M}}}_{0\to k}}({{\mathit{\boldsymbol{x}}}_{\rm{b}}}+\mathit{\boldsymbol{{x}'}})-{{\mathit{\boldsymbol{H}}}_{k}}{{\mathit{\boldsymbol{M}}}_{0\to k}}({{\mathit{\boldsymbol{x}}}_{\rm{b}}})\ \\ & (k=1, ..., s) \\ \end{align} $

(7)

其中，$ L_{k}^{'}({\boldsymbol{x}}') $表示各时刻的观测扰动。

$ \mathit{\boldsymbol{{y}'}}_{\operatorname{obs}, k}^{{}}={{\mathit{\boldsymbol{y}}}_{\operatorname{obs}, k}}-{{\mathit{\boldsymbol{H}}}_{k}}{{\mathit{\boldsymbol{M}}}_{0\to k}}({{\mathit{\boldsymbol{x}}}_{\rm{b}}})\ (k=1, ..., s) $

(8)

由于采用有限数量的集合样本构造背景误差协方差，不可避免产生远距离虚假相关，为使观测的影响保持在一定范围内，实施了局地化方案。最终得到迭代形式的分析增量解$ {\boldsymbol{x}'^{, i + 1}} $($i = 0, 1, 2, \cdots, {I_{\max }}$)，I_max为最大迭代次数，如下式：

$ \begin{align} & {{{\mathit{\boldsymbol{{x}'}}}}^{, i+1}}={{{\mathit{\boldsymbol{{x}'}}}}^{, i}}+\sum\limits_{k=1}^{s}{\left(\mathit{\boldsymbol{ }}\!\!\rho\!\!\rm{ }\circ {{\mathit{\boldsymbol{P}}}_{x}}\overset{\mathit{\boldsymbol{ }}\!\!\hat{\ }\!\!\rm{ }}{\mathop{{{\mathit{\boldsymbol{P}}}_{\mathit{\boldsymbol{y, k}}}}}}\, \right)L_{k}^{'}\left({{{\mathit{\boldsymbol{{x}'}}}}^{, i}} \right)}+ \\ & \sum\limits_{k=1}^{s}{\left(\mathit{\boldsymbol{ }}\!\!\rho\!\!\rm{ }\circ {{\mathit{\boldsymbol{P}}}_{x}}\overset{\mathit{\boldsymbol{ }}\!\!\tilde{\ }\!\!\rm{ }}{\mathop{{{\mathit{\boldsymbol{P}}}_{\mathit{\boldsymbol{y, k}}}}}}\, \right)\left[ \mathit{\boldsymbol{{y}'}}_{\operatorname{obs}, k}^{{}}-L_{k}^{'}\left({{{\mathit{\boldsymbol{{x}'}}}}^{, i}} \right) \right]} \\ \end{align} $

(9)

从公式 (9) 可以看出，NLS-4DVar方法通过迭代的形式求得最终的分析增量，其中，ρ表示局地化系数，$ \circ $表示两个矩阵的舒尔积，${\boldsymbol{\widehat P}_{y, k}} $和$ {\boldsymbol{\widetilde P}_{y, k}} $的定义可以参阅Tian and Feng (2015)。公式 (9) 中观测扰动L_k(x_a^{′, i}) 是利用非线性模式M_{0→ k}和非线性观测算子H_k直接计算的，因而很大程度上保持了原有代价函数的非线性性质。当预报模式M_{0→ k}和观测算子H_k为线性情况下，NLS-4DVar可退化为POD4DEnVar方法。NLS-4DVar的第一步迭代形式也是POD4DEnVar分析增量解，因而在上述提到的高度非线性情况下，NLS-4DVar通过增加迭代的次数，可以进一步优化POD4DEnVar的同化结果。

2.2 NRAS同化系统构建

基于NLS-4DVar理论，本文引入雷达资料的观测算子 (Sun and Crook, 1997, 1998)，对PRAS进行了改进，以WRF模式为预报模式，建立了雷达资料同化系统NRAS，实现了对雷达资料的反射率和径向风同化，参照Wang et al.(2012)和Zhang et al.(2015)的文章，NRAS中选取1小时作为雷达资料的同化时间窗口长度。为了验证该系统的同化效果，设计了OSSEs试验和针对两次暴雨个例的实际雷达资料同化试验。

3 雷达资料同化系统NRAS的验证 3.1 OSSEs验证

OSSEs试验是检验同化系统的有效手段，本节首先利用OSSEs试验对基于NLS-4DVar方法建立的雷达资料同化系统NRAS进行验证，并分析其相对于PRAS的优势。

(1) 资料和模式设置

WRF模式中使用的资料有NCEP的1°×1°的再分析资料，同化的雷达资料为武汉站的径向风和反射率的“模拟观测”。模式水平方向的格点数是400×400，格距为6 km，垂直方向为不等距27层，模式层顶为50 hPa。模式物理过程选用RRTM长波辐射方案 (Mlawer et al., 1997)、Dudhia短波辐射方案 (Dudhia, 1989)、YSU边界层方案 (Hong et al., 2006)、WSM6微物理方案 (Hong and Lim, 2006)，Noah LSM的陆面方案 (Ek et al., 2003) 以及无积云参数化方案。

(2) 试验设计

试验中的分析时间为2010年7月8日00时 (协调世界时，下同)，同化时间窗为7月8日00时到7月8日01时。为了详细考察NRAS的同化效果，OSSEs中设计了4个试验，分别为CTRL、PRAS、NRAS2和NRAS3，其中，CTRL为不同化雷达资料的控制试验，试验PRAS、试验NRAS2和试验NRAS3分别为NRAS迭代一次、两次和三次的同化试验。通过对比NRAS2、NRAS3和PRAS的同化效果，可以考察NRAS相对于PRAS的优势。

为了产生“模拟观测”。在OSSEs试验中，通过WRF-3DVar的三维变分控制变量的扰动方法 (Barker et al., 2004; Barker, 2005) 模块得到扰动的背景场作为“真实初始场”，并积分1小时得到同化窗口内的“真实场”。由于我国多普勒雷达的扫描频率为6分钟一次，因而模式产生的“真实场”也是6分钟输出一次。在此基础上，选取武汉站的实际雷达位置作为观测点，并通过观测算子和同化窗口内的“真实场”构造相应的雷达“模拟观测”，最后在径向风和反射率“模拟观测”上分别加入均值为0，标准差分别为1 m s⁻¹和1 dBZ的随机误差，形成同化的“模拟观测”。雷达体扫的9个仰角分别为0.5°、1.5°、2.4°、3.4°、4.3°、6.0°、9.9°、14.6°和19.5°。

无论是NRAS或者PRAS，均需要产生模式扰动和相应的观测扰动，按照Zhang et al.(2015)的试验，围绕同化窗口进行两个较长时间积分过程，利用滑动采样的方式 (Tian et al., 2008, 2011, 2014) 得到相应的集合样本。最后基于背景场和集合样本，产生相应的模式扰动和观测扰动。

(3) 试验结果

为了衡量各试验的同化效果，尤其是对NRAS和PRAS的同化性能进行对比，计算了各试验初始时刻的分析场与“真实场”的均方根误差，均方根误差越小，表明同化结果越接近于“真实场”，即同化效果越好。图 1给出了各变量均方根误差的垂直分布图。从图中可以就看出，相对于试验CTRL (黑色方块线)，在大部分垂直层中，试验PRAS (红色三角线)、NRAS2(绿色圆圈线) 和NRAS3(蓝色加号线) 的均方根误差均明显减小，尤其是中层。这表明NRAS通过同化雷达资料有效提高了初始场的精度。进一步分析同化雷达资料的3个试验发现，相对于试验PRAS，NRAS2和NRAS3均方根误差更小，这表明相对于PRAS，NRAS能够进一步提高同化效果，有效解决高度非线性情况下PRAS的应用问题。同时还可见试验NRAS3的均方根误差小于试验NRAS2，但二者的差异明显小于试验PRAS和NRAS2，这表明随着迭代次数的增加，NRAS对同化效果的改善越来越小，一定次数的迭代便可取得近似收敛的结果。通过两次迭代，NRAS即可得到近似收敛的结果。

以上的OSSEs试验表明相对于PRAS，NRAS能够进一步提高雷达资料的同化效果，改善高度非线性情况下PRAS的同化效果，且有限的迭代次数即可得到近似收敛的结果。

3.2 实际雷达资料同化试验

为了考察NRAS对实际雷达资料的同化性能，本文选择了两个强降水事件作为研究个例，分别为2010年7月8日发生在我国中部地区的24小时降水过程和2014年3月30日发生在我国华南地区的6小时降水过程，针对这两个降水事件进行了同化及预报，验证相对于PRAS，NRAS是否能够进一步改善初始场和提高降水预报效果。

3.2.1 2010年7月8日的降水事件

(1) 降水事件简介

2010年7月8日，我国中部地区发生了强降水过程，图 2a为24小时累积降水实况。其中湖北东部的降水最强，最强降水24小时达到了287 mm，其次在湖北和湖南交界、安徽和江西交界也存在强降水。大部分地区，24小时降水量均超过了90 mm。

(2) 试验设计

试验中的使用资料和模式设置同OSSEs。为了考察NRAS对实际雷达资料的同化效果，同OSSEs试验一样，设计了四个实际雷达资料同化试验，分别为试验HCTRL、试验HNLS1、试验HNLS2和试验HNLS3，其中试验HCTRL为不同化雷达资料的控制试验，以2010年7月8日00时的背景场积分24小时得到预报场，试验HNLS1、试验HNLS2和试验HNLS3分别为NRAS迭代一次、两次和三次的同化试验，考虑到模式的spin-up问题 (张斌等，2014)，模式首先积分6小时，然后将积分6小时后的预报场作为同化试验的背景场进行雷达资料的同化，然后将同化后的分析场再进行24小时的预报。试验中同化的雷达资料为武汉站的实际雷达观测的径向风和反射率资料。

(3) 降水预报结果

对比24小时观测降水 (图 2a)，可以看出试验HCTRL (图 2b) 没有预报出湖北东部的强降水，而在湖北东南部却预报出了最大降水中心超过240 mm大量虚假降水。试验HNLS1、HNLS2和HNLS3的预报结果如图 3所示。图 3a为HNLS1的预报降水，即PRAS同化雷达资料后的预报降水，相对于试验HCTRL，其增强了湖北东部降水和削弱了东南部虚假降水，预报结果更接近于实况，说明HNLS1明显的提高了降水预报效果，详细的分析可见Zhang et al.(2015)。对比试验HNLS1，试验HNLS2和HNLS3加大了降水强度，使预报结果更加接近实况。这与前面OSSEs的试验结果一致，相对于PRAS，NRAS能够更好的改善降水预报效果。同时，相对于试验HNLS2，试验HNLS3使得安徽西部的降水略微西移，更加接近观测，但是改进很小，试验HNLS2和试验HNLS3的预报结果基本相同。这表明，随着迭代次数的增加，对预报的改善逐步减小，NRAS两步迭代已可取得近似收敛的结果。

为了定量客观的评价降水预报效果，采用SAL评分 (Wernli et al., 2008) 对降水预报结果进行了评估。在SAL评分中，S为衡量降水结构的指标，A为衡量降水强度的指标，L为衡量降水位置的指标。总的来说，S、A和L的绝对值越小，表明降水预报效果越好，其中L指标最重要 (公颖，2010)。各试验的SAL降水评分结果见图 4，可以看出相对于HCTRL，试验HNLS1、HNLS2和HNLS3的L值逐步减小，这表明同化试验使得降水落区预报更加接近于实况。所有试验A值均为负值，表明各试验预报降水均偏弱，然而相对于HCTRL，试验HNLS1、HNLS2和HNLS3的A值逐步增大，更接近于0，这也表明各同化试验对于降水强度有了改善。大部分同化试验的S值相对HCTRL也更接近于0。上述结果表明，NRAS的同化效果要优于PRAS，随着迭代次数的增加，对于降水预报的改善越来越小，趋于收敛，这些结论与上述OSSEs是一致的。

3.2.2 2014年3月30日的降水事件

(1) 降水事件简介

图 5a为2014年3月30日00时到06时的降水实况，可以看到强降水主要集中在广东地区，在23°N附近出现了一条较强的准纬向雨带，其中最强降水发生在广州市附近，中心最大降水量达到140mm以上，因而此次降水属于短时间的强降水事件。

(2) 资料和模式设置

模式中的模拟区域为两层双向嵌套，其中粗网格的水平格点数为300×200，水平分辨率为18 km，细网格的水平格点数为424×241，水平分辨率为3 km，粗细网格的垂直方向采用30层地形跟随坐标，模式层顶为50 hPa。粗网格的物理方案选择为：RRTM长波辐射方案 (Mlawer et al., 1997)、Dudhia短波辐射方案 (Dudhia, 1989)、YSU边界层方案 (Hong et al., 2006)、WSM6微物理方案 (Hong and Lim, 2006)、Noah LSM的陆面方案 (Ek et al., 2003) 和Kain-Fritsch积云参数化方案 (Kain，2004)。细网格没有使用积云参数化方案，其他方案同粗网格。雷达资料为广州站和梧州站的实际雷达观测的径向风和反射率资料。

(3) 试验设计

试验中的分析时间为2014年3月30日00时，同化时间窗为3月30日00时到3月30日01时。同样设计了四个同化试验，分别为不同化雷达资料的控制试验GCTRL和分别为NRAS迭代一次、两次和三次的同化试验GNLS1、试验GNLS2和试验GNLS3。考虑到模式的spin-up问题 (张斌等，2014)，模式首先积分6小时，然后将积分6小时后的预报场作为背景场进行控制试验和同化试验，所有同化试验均在细网格中进行，并对细网格结果进行分析。

(4) 降水预报结果

对比6小时降水实况 (图 5a)，可以看出试验GCTRL (图 5b) 预报的降水位置过于偏北，降水强度偏弱，雨带范围偏小，另外在广西和湖南交界处还出现了局地性的虚假强降水。图 6a为试验GNLS1的预报降水，即PRAS同化雷达资料后的预报降水，相对于试验GCTRL，广东境内的预报降水整体南移，且降水在纬向范围内变大，同时广西和湖南交界处出现的虚假降水得到了明显的削弱，试验GNLS1的预报降水与实况更为接近。对比试验GNLS 1，试验GNLS2(图 6b) 和试验GNLS3(图 6c) 使得广东境内的雨带进一步向南移动，同时使得GNLS1在 (24°N，113°E) 附近区域预报的强降水范围减小，更加接近实况。这再次说明相对于PRAS，NRAS能够更好的改善降水预报效果。相对于试验GNLS2，试验GNLS3使雨带进一步南移，从而更加接近观测，但是改善的幅度逐步减小。

图 7为各试验的SAL降水评分结果。相对于试验GCTRL，试验GNLS1、GNLS2和GNLS3的S、A和L值更为接近于0，这表明同化试验有效的提高了降水的预报结果。试验GNLS1、GNLS2和GNLS3的S、A和L值逐步减小，这再次表明NRAS的同化效果要优于PRAS，随着迭代次数的增加，能够逐步提高降水的预报。

3.2.3 初始场增量

上述同化试验对于降水预报的改善直接源于对初始场的调整，下面以2010年7月8日的降水事件为例 (2014年3月30日的降水事件结论类似)，分析各同化试验对于初始场的调整。

图 8为模式第8层 (约850 hPa)，试验HNLS1相对于试验HCTRL (图 8a和b)、试验HNLS2相对于试验HNLS1(图 8c和d)、试验HNLS3相对于试验HNLS2(图 8e和f) 的风场 (图 8a、c、e) 和雨水混合比 (图 8b、d、f) 的增量场。由于试验HNLS1的同化等价于PRAS的同化，因而图 8a和b即为PRAS同化的增量场，具体分析可参阅Zhang et al.(2015)。下面着重分析图 8c-f。相对于图 8a和b，可以看到图 8c、e的风场增量和图 8d、f的雨水混合比增量均明显减小，这表明NRAS的第一步迭代，即PRAS对同化产生的影响最大；其次从图 8c和d中可以看出HNLS2相对于HNLS1的增量仍较大，例如湖北东南部的风场辐散进一步增强 (图 8c)，垂直速度 (图 8c) 和雨水混合比 (图 8d) 进一步减弱，这与上述同化预报试验中湖北东南部虚假降水的进一步减弱有着很好的对应；安徽西部的风场辐合增强 (图 8c)、垂直速度 (图 8c) 和雨水混合比 (图 8d) 增大也与安徽西部降水的增强对应。这表明HNLS2能够进一步调整初始场，NRAS第二步迭代对初始场进一步的调整是具有物理意义的，能够进一步提升降水预报精度。从图 8e和f可以看出试验HNLS3相对于试验HNLS2的增量场变得非常小，尤其是垂直速度 (图 8e) 和雨水混合比 (图 8f)。这表明随着迭代次数的增加，对初始场的改善越来越小，有限次数的迭代即可得到近似收敛的结果。

4 结论

本文基于NLS-4DVar方法，对雷达资料同化系统PRAS进行改进，建立了新的雷达资料同化系统NRAS。为了考察NRAS相对于PRAS是否能够进一步提高同化效果，首先设计了OSSEs试验。OSSEs试验结果表明：相对于PRAS，NRAS能够进一步提高同化效果，随着迭代次数增加，改善越来越小，有限次数的迭代便可以取得近似收敛的结果。这表明了NRAS能够有效解决PRAS在高度非线性情况下的应用问题。在此基础上，针对两次次实际暴雨过程，进一步开展了实际雷达资料的同化试验，考察NRAS在实际雷达资料同化应用中的性能。结果同样表明，相对于PRAS，NRAS能够进一步改善初始场，通过对初始场中的风场和水汽场的有效调整，提高了降水强度和位置的预报精度。但随着迭代次数的增加，对于初始场的调整变小，进而对降水预报效果的改进也减小，NRAS通过有限次数的迭代，即可得到近似收敛的结果。这与OSSEs试验的结论是一致的。

中小尺度的天气发展过程往往具有较强的非线性，NRAS在同化雷达资料时考虑到了非线性问题，其通过迭代同化的过程逼近非线性的分析解，因而NRAS在数值预报中能够更为有效的同化雷达资料，从而更好地改善初始场的中小尺度结构，这对于提高中小尺度天气的精确预报水平具有重要意义。