1 引言

数值天气预报作为初边值问题，对初始场极为敏感，初始场的质量会极大地影响数值预报的准确率，这在中小尺度天气系统的数值预报中体现得尤为明显。若数值预报的初始场主要依靠常规观测资料得到的背景场提供，将会缺少对中小尺度天气系统的细致描述，从而降低中小尺度系统的预报准确率。

多普勒雷达资料具有高时空分辨率的特征，其中的径向风和反射率资料与模式背景场中的三维风场和水物质场有一定的关系，因此将这两种观测量同化进数值模式能够丰富背景场中的中小尺度信息。国内外学者为开发这两类观测量的同化模块（Gao et al., 1999; Xiao et al., 2007），进行了大量的研究。已有的研究表明，同化雷达径向风资料后，初始场会出现明显的中小尺度特征（盛春岩等，2006；Zhang et al., 2009），这有利于增加初始场的动力特性（Li et al., 2012），减少模式的spin-up时间（Xue et al., 2003），改善降水的位置和强度预报的准确率（Xiao et al., 2005;邱学兴和Zhang，2016）。对雷达反射率资料的同化主要有两种方式，一种是直接同化，结果可以通过改进模式初始场中的温、湿参数，增加云水含量，提高降水预报的效果（徐广阔等，2009）；另一种方式是同化雷达反射率反演的降水率，其结果也能够清晰地模拟出对流系统的发展、消亡过程（张诚忠等，2012；高郁东等，2015），总体而言，雷达反射率同化在各类天气系统的数值模拟中均获得了较好的效果（Xiao and Sun, 2007; Jung et al., 2008; Zhang et al., 2014）。相比于单独同化径向风或反射率因子，联合同化这两种观测量则可以结合两者的优点，使初始风场和雷达回波更接近实况（Hu and Xue, 2007;兰伟仁等，2010；刘寅，2011；Dong and Xue, 2013）。

虽然雷达资料同化能够在一定程度上改进初始场，使初始场更接近大气的真实状态，但由于大气是混沌的，具有随机性和不稳定性（Lorenz, 1963），因此初始场的微小偏差在模式积分过程中可以有不同的演变和发展过程。同化增量或同化不同的雷达观测要素造成的初始场差异可以认为是模式初始场的偏差，这些偏差是否对未来模拟结果产生什么影响？依赖于这些偏差在模式积分过程中的演变，而其演变规律既与偏差本身的分布有关，也与背景场结构有关，还与天气系统发展中的物理过程有关。已有研究表明，湿对流引起的潜热释放会给系统带来更多的发展扰动能量（Spyksma and Bartello, 2008），导致小振幅、小尺度初始误差在湿物理过程作用下快速增长（Zhang et al., 2002, 2003）。资料同化的目的是使初始场达到“最优”，但误差最小的初始场，由于小扰动增长，并不能保证在整个积分过程中的误差始终保持在较小的状态。所以为了进一步提高中小尺度天气系统的数值预报准确率，本文从初始可预报性的角度，将同化不同雷达观测要素产生的初始场偏差作为初始扰动，分析了偏差在模式积分过程中的演变特征，目的是为进一步发展资料同化技术提供理论依据。

2 同化方法和同化资料 2.1 四维集合变分同化方法

本文使用的四维集合变分同化方法（POD- 4DEnVar）（Tian et al., 2008, 2011, 2015）是在四维变分（4DVar）的基础上，通过引入具有流依赖特性的集合预报信息，并利用本征正交分解（POD）技术对观测扰动进行分解，将4DVar和集合卡尔曼滤波（EnKF）两种同化方法相结合，显式给出同化增量。

增量形式的4DVar同化目标泛函（Courtier et al., 1994）可表示为

$ \begin{array}{l} J(\mathit{\boldsymbol{x'}}) = \frac{1}{2}{(\mathit{\boldsymbol{x'}})^{\rm{T}}}{\mathit{\boldsymbol{B}}^{ - 1}}(\mathit{\boldsymbol{x'}}) + \frac{1}{2}{\left[ {\mathit{\boldsymbol{y'}}(\mathit{\boldsymbol{x'}}) - ({\mathit{\boldsymbol{y}}^{{\rm{obs}}}})'} \right]^{\rm{T}}} \cdot \\ {\mathit{\boldsymbol{R}}^{ - 1}}\left[ {\mathit{\boldsymbol{y'}}(\mathit{\boldsymbol{x'}}) - ({\mathit{\boldsymbol{y}}^{{\rm{obs}}}})'} \right], \end{array} $

(1)

式中，$ J(\mathit{\boldsymbol{x'}})$表示增量$\mathit{\boldsymbol{x'}}{\mathbf{ = }}\mathit{\boldsymbol{x}}-{\mathit{\boldsymbol{x}}_{\mathbf{b}}}$的代价函数，$\mathit{\boldsymbol{x'}} \in {\mathfrak{R} ^{n \times n}} $表示要求取的分析增量，n为模式变量的维数，x和x_b分别表示分析变量和背景场向量，$ {\mathit{\boldsymbol{B}}^{ - 1}} \in {\mathfrak{R} ^{n \times n}}$表示背景误差协方差矩阵的逆，$ {\mathit{\boldsymbol{R}}^{ - 1}} \in {\mathfrak{R} ^{n \times n}}$表示观测误差协方差矩阵的逆。$\mathit{\boldsymbol{y'}}(\mathit{\boldsymbol{x'}}) $和$({\mathit{\boldsymbol{y}}^{{\rm{obs}}}})' $表示模拟观测扰动和观测新息，其表达式为

$ \mathit{\boldsymbol{y'}}{\mathbf{(}}\mathit{\boldsymbol{x'}}{\mathbf{)}} = \left[ \begin{array}{l} {{\mathit{\boldsymbol{y'}}}_1}{\mathbf{(}}\mathit{\boldsymbol{x'}}{\mathbf{)}}\\ {{\mathit{\boldsymbol{y'}}}_2}{\mathbf{(}}\mathit{\boldsymbol{x'}}{\mathbf{)}}\\ {\rm{ }} \vdots \\ {{\mathit{\boldsymbol{y'}}}_s}{\mathbf{(}}\mathit{\boldsymbol{x'}}{\mathbf{)}} \end{array} \right], $

(2)

$ {\mathit{\boldsymbol{y'}}_t}(\mathit{\boldsymbol{x'}}) = {\mathit{\boldsymbol{y}}_t}({\mathit{\boldsymbol{x}}_b} + \mathit{\boldsymbol{x'}}) - {\mathit{\boldsymbol{y}}_t}({\mathit{\boldsymbol{x}}_b}), {\rm{ }}t = 1, 2, \cdots, s, $

(3)

$ ({\mathit{\boldsymbol{y}}^{{\rm{obs}}}})' = \left[ \begin{array}{l} (\mathit{\boldsymbol{y}}_{\rm{1}}^{{\rm{obs}}})'\\ (\mathit{\boldsymbol{y}}_{\rm{2}}^{{\rm{obs}}})'\\ {\rm{ }} \vdots \\ (\mathit{\boldsymbol{y}}_s^{{\rm{obs}}})' \end{array} \right], $

(4)

$(\mathit{\boldsymbol{y}}_t^{{\rm{obs}}})' = \mathit{\boldsymbol{y}}_t^{{\rm{obs}}} - {\mathit{\boldsymbol{y}}_t}({\mathit{\boldsymbol{x}}_b}), {\rm{ }}t = 1, 2, \cdots, s, $

(5)

${\mathit{\boldsymbol{y}}_t}(\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{H}}_t}\left[ {{\mathit{\boldsymbol{M}}_{1 \to t}}(\mathit{\boldsymbol{x}})} \right], $

(6)

式中，H_t表示$t $时刻的观测算子，${\mathit{\boldsymbol{M}}_{1 \to t}}$表示预报模式，$\mathit{\boldsymbol{y}}_t^{{\rm{obs}}} $表示t时刻的观测值，s表示同化时间窗内总的观测时次，t表示同化时间窗内的观测时刻。该方法的具体流程如下：

首先，制作集合预报。在分析时刻利用N个集合成员$({\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}, \cdots, {\mathit{\boldsymbol{x}}_N})$构造相对于背景场${\mathit{\boldsymbol{x'}}_{\rm{b}}}$的N个集合扰动$({\mathit{\boldsymbol{x'}}_{1}}, {\mathit{\boldsymbol{x'}}_{2}}, \cdots, {\mathit{\boldsymbol{x'}}_{N}})$，记为模式扰动${\mathit{\boldsymbol{x}}^p} $，其中${\mathit{\boldsymbol{x'}}_i} = {\mathit{\boldsymbol{x}}_\mathit{\boldsymbol{i}}} - {\mathit{\boldsymbol{x}}_b}(i = 1, 2, \cdots, N)$。

其次，对观测扰动和模式扰动进行POD分解。在同化时间窗口内的s个观测时刻，利用每个时刻的$N $个集合成员，按照公式（2）产生$N \times s$个观测新息$({\mathit{\boldsymbol{y'}}_{1, t}}, {\mathit{\boldsymbol{y'}}_{2, t}}, \cdots, {\mathit{\boldsymbol{y'}}_{N, t}}), t = 1, 2, \cdots, s$，记为观测扰动${\mathit{\boldsymbol{y}}^p}$，对${\mathit{\boldsymbol{y}}^p}$作POD分解：

${({\mathit{\boldsymbol{y}}^p})^{\rm{T}}}{\mathit{\boldsymbol{y}}^p} = \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^2}{\mathit{\boldsymbol{V}}^{\rm{T}}}, $

(7)

其中，V表示分解后的特征向量，${\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^2}$表示相应的特征值。在假设模式扰动和观测扰动为线性关系的基础上，对观测扰动${\mathit{\boldsymbol{y}}^p}$和模式扰动${\mathit{\boldsymbol{x}}^p}$作如下转换：${\mathit{\boldsymbol{P}}_y} = $ ${\mathit{\boldsymbol{y}}^p}\mathit{\boldsymbol{V}}$，${\mathit{\boldsymbol{P}}_x} = {\mathit{\boldsymbol{x}}^p}\mathit{\boldsymbol{V}}$，则分析增量$\mathit{\boldsymbol{x'}}$和对应的观测扰动$\mathit{\boldsymbol{y'}}$可以表示成${\mathit{\boldsymbol{P}}_x}$和${\mathit{\boldsymbol{P}}_y}$的线性权重相加，即如下形式：

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

(8)

$\mathit{\boldsymbol{y'}} = {\mathit{\boldsymbol{P}}_y}\beta, $

(9)

其中，$\beta = {({\beta _1}, {\beta _2}, \cdots, {\beta _N}{\rm{)}}^{\rm{T}}}$为线性组合系数。

最后，计算分析增量。将公式（8）、（9）带入公式（1），可得分析增量的显式表达式为

$\mathit{\boldsymbol{x'}} = {\mathit{\boldsymbol{P}}_x}{\left[ {(N - 1)\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{P}}_y^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{P}}_y}} \right]^{ - 1}}\mathit{\boldsymbol{P}}_y^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}({\mathit{\boldsymbol{y}}^{{\rm{obs}}}})', $

(10)

其中，背景误差协方差B的表达式如下：

$\mathit{\boldsymbol{B}} = \frac{{{\mathit{\boldsymbol{P}}_x}\mathit{\boldsymbol{P}}_x^{\rm{T}}}}{{N - 1}}.$

(11)

集合样本的构造是POD-4DEnVar方法实施的基础，为了产生更多的集合成员，本文使用历史采样法（Wang et al., 2010; Shen et al., 2015）构造，得到142组集合成员。

2.2 雷达资料观测算子

在进行雷达资料同化时，由于同化的雷达径向风和反射率不是模式变量，因此直接同化这两种观测变量需构建观测算子，将模式变量转化到观测空间。雷达径向风V_r和反射率Z的观测算子（Sun and Crook, 1997, 1998）如下：

${V_r} = u\frac{{x - {x_{{\rm{radar}}}}}}{d} + v\frac{{y - {y_{{\rm{radar}}}}}}{d} + (w - {V_{{\rm{Tm}}}})\frac{{z - {z_{{\rm{radar}}}}}}{d}, $

(12)

$z = 43.1 + 17.5{\rm{log}}(\rho q), $

(13)

式中，$\left({{\rm{ }}u{\rm{, }}v{\rm{, }}w{\rm{ }}} \right) $表示大气三维风场，d表示坐标位置为$\left({{\rm{ }}x, y{\rm{ }}, z{\rm{ }}} \right) $的模式格点与坐标位置为$({x_{{\rm{radar}}}}, {y_{{\rm{radar}}}}, $ ${z_{{\rm{radar}}}})$雷达站的距离，$ \rho $表示空气密度，$ q$为雨水混合比，${V_{{\rm{Tm}}}}$为雨滴末端下降速度，计算如下：

${V_{{\rm{Tm}}}} = 5.40a{q^{0.125}}, $

(14)

$a = {(\frac{{{p_0}}}{{\overline p }})^{0.4}}, $

(15)

式中，$ a$表示修正系数，${\bar p} $表示基态气压，${p_0} $表示地表气压。

2.3 同化资料和模式设置

试验选择一次冬季暴雨个例，该次暴雨发生在2015年12月8日20时（协调世界时，下同）至9日20时，降水区域覆盖华东、华南多个省份。在广东省东部、福建省南部以及浙江省东南部出现了暴雨，部分站点如佛冈、新丰的24 h降水量超过100 mm，达到大暴雨量级，降水强度远超历史平均水平（杨雨轩等，2016）。

试验基于WRFV3.6模式，模拟区域如图 1所示。采用两重双向嵌套，内、外区域的水平分辨率分别为5 km和15 km，水平格点数为433×394和200×200，垂直方向分层为30层，模式层顶气压为50 hPa。模式微物理过程采用WSM6方案，积云对流参数化采用Betts-Miller-Janjic方案，边界层采用YSU方案，陆面过程采用Noah方案，近地面采用Monin-Obukhov方案，长波辐射采用RRTM方案，短波辐射采用Dudhia方案。两层区域的物理参数化方案相同。

多普勒雷达资料同化在内层区域进行，同化起始时刻为2015年12月8日19时，同化窗口为1 h，8日20时结束，同化时间窗内每6分钟同化一次雷达资料。同化结束后，将同化得到的8日19时分析场作为模式初始场继续向前积分至10日00时。

模式初始场和侧边界条件采用美国国家环境预报中心（NCEP）提供的每日四次的再分析资料，水平分辨率为1°×1°。降水资料采用中国地面降水资料与CMORPH降水产品融合的逐小时降水量网格数据集（1.0版本），水平分辨率为0.1°×0.1°。

雷达资料使用S波段多普勒天气雷达观测的径向风和反射率因子，雷达位于梅州站（位置：24.26°N，115.99°E；天线海拔高度：416 m）。

2.4 SAL降水评估

为了定量评估数值试验对降水的模拟效果，采用SAL（Structure, Amplitude, Location）评估方法（Wernli et al., 2008）对降水的结构（S）、强度（A）以及位置（L）进行检验。

2.4.1 结构S的检验

$S = \frac{{V({R_{{\rm{mod}}}}) - V({R_{{\rm{obs}}}})}}{{0.5\left[ {V({R_{{\rm{mod}}}}) - V({R_{{\rm{obs}}}})} \right]}}, $

(16)

$ V(R) = \frac{{\sum\limits_{n = 1}^m {{R_n}{V_n}} }}{{\sum\limits_{n = 1}^m {{R_n}} }}, $

(17)

$ {V_n} = \frac{{\sum\limits_{(i, j) \in {R_n}} {{R_{i, j}}} }}{{R_n^{{\rm{max}}}}}{\rm{ = }}\frac{{{R_n}}}{{R_n^{{\rm{max}}}}}, $

(18)

其中，${R_{i, j}}$表示格点$(i, j)$处的降水，$R $表示降水个体内的最大降水值，$R_n^{{\rm{max}}}$表示降水个体$ n$内的总降水量，${R_n}$表示降水个体$ n$内的总降水量，V_n为第n个降水个体内总降水量与最大降水量之比，V(R)为m个降水个体以该个体内总降水量为权重的V_n的加权平均，下标mod表示模式预报场，obs表示实况场。

S＞0表示所取区域内的预报范围相对实况较大或预报降水中心降水量小于实况，或两种情况同时存在，S＜0表示相反的情况。

2.4.2 强度A的检验

$A = \frac{{D({R_{{\rm{mod}}}}) - D({R_{{\rm{obs}}}})}}{{0.5[D({R_{{\rm{mod}}}}) + D({R_{{\rm{obs}}}})]}}, $

(19)

$D(R) = \frac{1}{N}\sum\limits_{(i, j) \in D} {({R_{i, j}})}, $

(20)

其中，D(R)为所取区域内非缺省格点上降水的平均值。A＞0表示所取区域内的平均降水预报强度偏强，A＜0表示平均降水预报强度偏弱。

2.4.3 位置L的检验

$L = {L_1} + L{}_2, $

(21)

${L_1} = \frac{{\left| {x({R_{{\rm{mod}}}}) - x({R_{{\rm{obs}}}})} \right|}}{d}, $

(22)

${L_2} = 2\frac{{\left| {r({R_{{\rm{mod}}}}) - r({R_{{\rm{obs}}}})} \right|}}{d}, $

(23)

$r = \frac{{\sum\limits_{n = 1}^m {{R_n} \cdot \left| {x - {x_n}} \right|} }}{{\sum\limits_{n = 1}^m {{R_n}} }}, $

(24)

其中，$x\left({{\rm{ }}R} \right) $为降水主体的重心位置，d为所取区域内非缺省值上格点间的最大距离。L的范围在0~2之间，L值越小表示预报结果越好。

2.5 试验设计

为了考察POD-4DEnVar方法同化不同的雷达观测要素对暴雨模拟的影响，设计了如表 1所示的4组试验。试验1不同化雷达资料，记为控制试验（CNTL试验）；试验2只同化雷达径向风，记为径向风同化试验（Vr试验）；试验3只同化雷达反射率，记为反射率同化试验（R试验）；试验4同时同化雷达径向风和反射率，记为径向风和反射率同化试验（Vr+R试验）。

3 不同试验模拟结果分析 3.1 24 h累积降水

图 2给出了12月8日20时至9日20时降水实况、控制试验以及3组同化试验模拟的24 h累积降水分布。与实况（图 2a）对比可见，4组试验大致都模拟出了降水落区的分布，但对降水强度的模拟却存在较大差异。其中，CNTL试验模拟的降水强度较实况明显偏强，在广东和福建两省沿海交界处，降水强度达到了大暴雨量级，在福建和浙江两省交界处也模拟出了暴雨量级的降水，明显强于该地的实况降水（图 2b）。Vr试验对广东和福建两省交界处暴雨区的模拟效果有一定改进，大暴雨区范围比CNTL试验明显缩小，福建和浙江两省交界处的降水模拟强度也有所减弱，但在广东、广西和湖南三省交界处模拟出一个大于50 mm的虚假暴雨中心（图 2c）。R试验基本上没有模拟出广东和福建两省交界处的暴雨，该地区的模拟降水只达到了大雨量级，在广东、广西两省交界处模拟出了一条东北—西南向的虚假雨带，虚假降水中心达到了暴雨量级（图 2d）。Vr+R试验在实况暴雨区所模拟的降水强度介于Vr试验和R试验之间，但该试验同样在主体暴雨区的西北部模拟出一条东北—西南向的虚假雨带（图 2e）。

为了定量评价降水模拟结果，对12月8日20时至9日20时的24 h累积降水进行了SAL检验，检验结果体现为S、A、L三个指标，分别表示对降水结构、强度和位置的模拟效果。在L值相近的情况下，S和A的绝对值越小，则表明降水模拟的效果越好，图 3给出了各试验的SAL评分。从反映降水位置模拟效果的L值看，CNTL试验的L值最小，R试验的L值最大，虽然同化试验没有改善降水落区的模拟效果，但4组试验对降水位置的模拟差异并不是很大；从反映降水结构模拟效果的S值看，CNTL试验的|S|值最大，之后依次为R试验和Vr试验，Vr+R试验的|S|值最接近0，这说明经过同化，降水结构的模拟效果得到了改善，其中Vr+R试验对降水结构的模拟效果最佳；从反映降水强度模拟效果的A值看，同化试验的A值均小于CNTL试验，其中R试验对CNTL试验偏强降水的削弱作用最明显。

从SAL降水定量评估结果可以看出，同化多普勒天气雷达资料能够改善降水模拟效果，但不同雷达观测要素的同化效果和改善内容也不同，对S、A、L这三个指标进行综合评估，可以认为Vr+R试验的降水模拟结果最好。

3.2 同化分析场与背景场的偏差

同化试验与控制试验的不同在于同化试验的初始场是在控制试验的初始场上叠加了同化增量，即控制试验的初始场是同化试验的背景场，所以可以认为初始场的差异是导致降水模拟结果不同的原因，为了揭示不同雷达观测要素的同化试验与控制试验的初始场差异，图 4给出了8日19时850 hPa控制试验的初始风场、水汽混合比，以及3组同化试验的相应同化增量场。由图可见，在背景场上，广东省、福建省和浙江省沿海一带以东南气流为主，在广东省和福建省交界处有偏南风和东南风的辐合，偏南气流的辐合不仅有利于激发上升运动，同时还为暴雨区提供了充沛的水汽，使降水区域内水汽混合比达到9 g kg^-1以上（图 4a）。同化雷达资料后（图 4b、c、d），广东省东部暴雨区呈现较强的偏西风增量，福建省南部暴雨区则表现为较强的西南风增量，这些偏西风增量一方面减弱了背景场的偏南风和偏东风在广东和福建两省交界处的辐合，另一方面有利于削弱来自海上的暖湿气流对水汽的输送。水汽混合比增量则表现为以雷达观测为中心的负增量，这表明同化后暴雨区的水汽条件明显减弱。

比较3组同化试验的增量场（图 4b、c、d）可以看出，同化不同的雷达观测要素均能修正背景风场和水汽场，进而改进初始分析场和降水模拟结果。

3.3 同化分析场之间的偏差

同化试验相对于控制试验的同化增量反映了同化分析场与背景场的差异，而不同雷达观测要素的同化对初始场也有不同的修正，从而导致了不同同化试验的降水模拟结果不同。为了揭示不同雷达观测要素的同化对降水模拟的影响，以下比较了3组同化试验的同化分析场之间的偏差，图 5给出了Vr+R试验分别与Vr试验和R试验之间的850 hPa水汽混合比、风场的偏差。分析混合比的偏差发现，相对于Vr+R试验，Vr试验的混合比增量表现为以雷达观测为中心的负值，而R试验为正值，这说明在暴雨区，Vr试验的混合比最大，R试验最小，Vr+R试验则介于Vr试验和R试验之间。对比同化背景场（图 4a），同化雷达资料后对背景场的水汽条件都有削弱作用，在3组同化试验中，Vr试验对背景场的水汽混合比削弱最小，R试验的削弱最大，Vr+R试验对背景水汽条件的削弱程度介于Vr试验和R试验之间。分析风场偏差可见，在暴雨区Vr+R试验与Vr试验的差异表现为偏西风增量，即Vr+R试验的海上暖湿气流弱于Vr试验；Vr+R试验与R试验的差异表现为偏东风增量，即Vr+R试验的海上暖湿气流强于R试验。这表明在3组同化试验中，Vr试验对背景风场的调整最小，R试验对背景风场的调整最大，Vr+R试验对背景风场的调整介于Vr试验和R试验之间。

通过比较3组同化试验相对于背景场的同化增量场以及3组同化试验之间的同化分析场偏差可以看出，3组同化试验的初始风场和水汽条件存在较小的差异，3组同化试验相对于控制试验的同化增量要远大于3组同化试验之间的差异，这表明同化不同种类的雷达观测要素对初始风场和混合比的调整呈现类似的特征。

4 不同试验初始偏差的演变

由于不同雷达观测要素的同化对模式初始场的修正不同，因此3组同化试验的初始场差异是导致降水模拟结果不同的直接原因。3组同化试验对背景场的修正，或3组同化试验之间的同化分析偏差对于同化背景场来说都可以认为是叠加在背景场上的扰动，其在预报时效内的影响决定于扰动（即增量或偏差）的演变，而其演变规律既与偏差的结构（分布）有关，又与背景环流有关，也与降水的演变有关。为了探究不同雷达观测要素的同化通过改变初始场，进而影响降水模拟结果的原因，以Vr+R试验作为“真实”试验，分析其与Vr试验和R试验之间的初始偏差场的演变。

4.1 逐小时降水偏差的演变

图 6给出了Vr+R试验分别与Vr试验和R试验模拟的逐小时降水偏差及偏差增长率随时间的演变，从图 6可以看出，在模式积分过程中，Vr试验和R试验相对于Vr+R试验的降水偏差随时间逐渐增大，且Vr试验的降水偏差及其增长率要大于R试验。值得注意的是，在模式积分至16 h之前，Vr试验和R试验的降水偏差增长率较小，降水偏差增长比较平稳；模式积分至16~20 h之间时，降水偏差增长率增大，降水偏差增长加快；模式积分至20~24 h，降水偏差增长再次变缓，降水偏差增长率减小。因此可以认为，3组同化试验的24 h累积降水模拟的差异源于初始差异在不同阶段的不同演变。

4.2 垂直速度和雨水混合比

为了进一步分析不同同化试验降水偏差演变的原因，图 7给出了Vr试验和R试验相对于Vr+R试验的垂直速度和雨水混合比偏差在垂直方向上的分布。分析图 7a和图 7b可以看出，经过模式的4 h积分，垂直速度场完成了模式的spin-up调整；在4~16 h的阶段，在500 hPa以下，偏差增幅较小；但当模式积分至16 h后，垂直速度偏差开始迅速增大，且在850~700 hPa附近增长最为明显。雨水混合比偏差出现在600 hPa以下，其垂直廓线演变也具有类似的特征（图 7c、d）。

由于垂直速度偏差和雨水混合比偏差的大值都出现在低层，为了更加直观定性的分析变量偏差的增长特征，图 8给出了850~700 hPa平均垂直速度偏差和雨水混合比偏差的逐时演变，分析图 8发现，在整个积分过程中，Vr试验的变量偏差始终大于R试验。此外，无论是Vr试验还是R试验，垂直速度偏差和雨水混合比偏差在模式积分至16 h后均大幅增长。这两点与降水偏差的逐时演变特征类似。这说明降水偏差的演变与850~700 hPa垂直速度偏差以及雨水混合比偏差的逐时演变有关，也就是说，在模式积分过程中，垂直速度偏差和雨水混合比偏差存在急剧增长，这种急剧增长正是逐时降水偏差迅速增加的主要原因，进而最终导致3组同化试验的24 h降水模拟结果存在较大的差异。

4.3 偏差能量演变

为了综合分析Vr试验和R试验相对于Vr+R试验变量偏差增长特征，参照文献（Zhang et al., 2003），定义了偏差能量（Difference Total Energy，简称DTE）：

${{E}_{\text{DT}}}=\frac{1}{2}\sum\limits_{i}{\sum\limits_{j}{\sum\limits_{k}{\left( u_{ijk}^{'2}+v_{ijk}^{'2}+\frac{{{c}_{p}}}{R}T_{ijk}^{'2} \right),}}}$

(25)

其中，i、j、k为x、y、z方向的格点数，${u'_{ijk}}$、${v'_{ijk}}$、${T'_{ijk}}$分别为Vr试验和R试验相对于Vr+R试验的风场和温度场的偏差，c_p=1004 J (kg K)^-1，R=287 J (kg K)^-1。

图 9给出了850~700 hPa的平均偏差能量及其增长率的时间演变，从图中可以看出，在初始时刻，Vr试验和R试验的偏差能量较小，但随着时间的演变，偏差能量持续增加，根据其演变特征，可分为4个阶段：第一阶段，0~4 h，由于模式积分处于spin-up阶段，因此偏差能量主要以调整为主；第二阶段，4~12 h，偏差能量随时间缓慢增长，相应的增长率波动上升；第三阶段，12~19 h，偏差增长率迅速增大，偏差能量以指数形式迅速增长；第四阶段，19~24 h，偏差增长率急剧减小，偏差能量增速减缓，表明偏差能量逐渐达到饱和状态。结合逐小时降水偏差演变（图 6）、垂直速度和雨水混合比偏差演变（图 8）发现，变量偏差迅速增长的阶段（16~20 h）恰好为偏差能量增速最大的时段，而且偏差能量迅速增长较降水偏差的迅速增长早4 h。

图 10为不同时刻Vr试验的850~700 hPa平均偏差能量、垂直速度和雨水混合比偏差分布。可以看出，广东和福建两省交界处的暴雨区为偏差能量大值区，且偏差能量在15~16 h增幅明显，大值中心从800×10⁵ m² s^-2迅速增长到1600×10⁵ m² s^-2，与此同时，暴雨区的垂直速度偏差和雨水混合比偏差也明显增大并沿雨带传播。当模式积分至19 h时，偏差能量增速达到峰值（图 9b），暴雨区内形成一条东北—西南走向的偏差能量大值带，垂直速度偏差和雨水混合比偏差进一步增长，其大值中心位于偏差能量梯度大值区。对R试验的变量偏差演变进行分析也可以得到类似的演变特征（图略）。

综合物理量偏差的时空演变可以看出，在偏差迅速增长阶段（16~20 h），偏差能量和变量偏差在时空分布上有较好的对应，垂直速度偏差和雨水混合比偏差随着偏差能量的增大而增大，变量偏差增长最明显的时段为偏差能量增幅最大的时段，变量偏差增长最明显的区域为偏差能量梯度较大的区域。这说明数值模式对初始场十分敏感，初始时刻的微小偏差会随着积分时间的延长不断放大，而且偏差在模式积分过程中可能存在急剧增长的过程，这可能就是初始时刻的微小差异带来巨大预报差异的直接原因。

5 结论与讨论

本文基于WRF模式，采用POD-4DEnVar同化方法，对一次冬季暴雨过程进行数值模拟，探讨了同化不同的雷达观测要素对暴雨模拟的影响，得到以下主要结论：

（1）同化多普勒天气雷达资料有利于削弱控制试验对降水偏强的模拟结果，改善降水分布结构；同化不同的雷达观测要素所得到的降水模拟结果不同，SAL结果表明同时同化径向风和反射率的降水模拟结果最好。

（2）同化多普勒天气雷达资料有利于调整初始风场和水汽条件，一方面减弱背景场的偏南风和偏东风在暴雨区的辐合，阻碍海上暖湿气流对暴雨区的水汽输送，另一方面直接削弱暴雨区的水汽条件，大幅降低水汽条件。通过修正初始动力场和水汽场达到改进降水模拟结果的目的。

（3）比较3组同化试验相对于控制试验的同化增量场以及3组同化试验之间的同化分析场差异发现，同化增量场要远大于同化试验之间的差异，这表明同化不同的雷达观测要素对初始风场和湿度场的调整呈现类似的特征，但3组同化试验的初始风场和水汽混合比存在较小的差异。

（4）分析初始场差异较小的3组同化试验降水模拟结果差异明显的原因发现，垂直速度偏差和雨水混合比偏差的最大值出现在低层，在模式积分至16 h，850~700 hPa的平均垂直速度偏差和雨水混合比偏差开始急剧增长，这种变量偏差的急剧增长与逐时降水偏差的迅速增加一致，这是降水偏差增长的直接原因。

（5）在偏差迅速增长阶段（16~20 h），垂直速度偏差和雨水混合比偏差随着偏差能量的增大而增大，变量偏差增长最明显的时段为偏差能量增幅最大的时段，变量偏差增长最明显的区域为偏差能量梯度较大的区域，偏差能量迅速增长早于变量偏差和降水偏差的迅速增长。

基于POD-4DEnVar的单部多普勒雷达资料同化已应用于夏季的梅雨锋暴雨（Zhang et al., 2015）并取得了较好的效果，为了检验该方法对冬季暴雨的适用性，本文也使用了单部雷达的资料。但是从理论上讲，对雨区及附近的多部雷达资料进行同化可以应用更多的观测信息，效果应该会更好。目前也有文献对多部雷达资料同化进行了研究，且取得了优于单部雷达资料同化的效果。未来可以考虑基于POD-4DEnVar方法的多部雷达资料同化，以期进一步提高同化效果。