引言

初始场是影响数值天气预报准确率的关键，随着数值模式分辨率的逐步提高以及非常规资料的应用越来越广，在中尺度模式中加入非常规资料改善初始场，从而提高数值天气预报的准确率成为气象科学家研究的热点。在众多的非常规观测中，具有高时空分辨率的多普勒雷达能够更为精确地捕捉到中小尺度天气系统的实时发展和演变，从而得到广泛的应用(许小永等，2005；智协飞等，2010)。

多普勒雷达基本观测资料包括反射率、径向风速和速度谱宽。现有的雷达径向风资料同化方法主要分为两种，间接同化与直接同化(杨毅等，2007)。间接同化方法根据雷达观测信息，在一定的假设基础上，建立风场模型，通过各种方法推演出较为完整的风场信息。应用较广的方法是利用雷达坐标系与笛卡尔坐标系之间的空间几何关系反演常规风场。Liang(2007)在区域方位角内风场均匀的假定下提出IVAP(Integrating Velocity Azimuth Process)方法，将变分-积分思想引入雷达径向风反演方程组，结果表明IVAP方法在较好描述中小尺度风场的同时对多种传统反演方法具有一定的兼容性, 因此本文对于雷达径向风同化算子的改进工作便是基于IVAP反演方程组。根据分析单元以及假设条件不同能够衍生出大量的反演方法：VAD(Lhermitte et al, 1961; 周小刚等，2015), VAP(陶祖钰，1992)，VPP(郎需兴等，2001)等。尽管间接同化能够将完整的风场信息引入同化系统，但是各种反演方法需要引入局地强约束条件来满足反演方程组闭合问题，从而会产生大量的反演误差，影响分析场的准确性。

另一种是直接同化，利用雷达径向风观测算子将背景常规风场处理成径向风速形式后与观测径向风极小化求解最优风场。Sun等(1997)和Sun(1998)建立VDRAS(Variational Doppler Radar Analysis System)系统，利用四维变分方法以及云尺度数值模式，采用雷达径向风传统算子得到中尺度系统细致结构，由于四维变分方法将预报方程组引入代价函数中作为约束，从而产生的分析风场在热力与动力上更为协调。Xiao等(2005)在MM5_3Dvar系统中添加Richardson方程在校正垂直速度的同时对分析场施加了较强的热力、动力约束条件，模拟结果表明加入Richardson后，分析风场的垂直速度的量级得到校正，提升了短时强降水模拟结果。Hu等(2006)在ARPS_3Dvar系统的代价函数中加入三维质量散度的约束项来产生切向风场信息，结果表明加入三维散度约束项的同化结果对于龙卷、风暴等强对流系统有了更好的模拟，然而在散度项中的权重系数是由经验给定而缺少相应的物理基础。刘红亚(2009)利用GRAPES_3Dvar系统，对比分析采用不同方法引入垂直速度对同化结果的影响，一种是利用Richardson方程直接诊断垂直速度增量，另一种则是将Richardson方程作为观测算子同化由反射率反演得到的垂直速度，结果表明前者得到的垂直速度较小，仅代表了天气尺度系统垂直运动的特征，而后者将反演得到的中小尺度的垂直速度作为观测项，对于中小尺度的动力场有了更为细致的反映。国内大量研究者采用直接同化雷达径向风的方法进行模拟试验(张林等，2006；许小永等，2006；陈力强等，2009; 李媛等，2011；陈锋等，2012；秦琰琰等，2012；2014；陈敏等，2014；罗义，2012；肖现等，2013；张旭斌等，2015)，结果表明直接同化雷达径向风速后增加了初始场的中小尺度信息，增强了局地风场的辐合辐散特征，对降水等要素预报具有明显的改进作用。上述研究均采用径向风传统观测算子，在同化过程中附加了相应的约束项，一方面是为了调节由于中小尺度动力场单一改变而引起其他物理量不协调，另一方面由于传统观测算子仅能引入径向风信息，分析场中切向风增量通过背景误差协方差矩阵调节得到，约束项则利用其物理属性对切向风增量进行相应的调整。虽然约束项减弱了分析风场的欠定性，但是切向增量依旧不是由观测信息得到，以及大部分约束项中依旧存在依照经验值给定权重的干扰因素，因此采用传统观测算子产生的分析风场是存在误差的。

鉴于以上问题，罗义等(2014)提出了一种新的雷达径向风资料同化观测算子，将多普勒雷达径向风速以及径向风速的空间变化项作为同化变量引入同化系统，利用同化系统的最优化方法迭代产生分析风场。这种方法不仅极大程度地引入了真实完整的观测信息，同时利用同化系统本身作为约束求解最优风场，降低计算误差。更重要的是，区域的分析单元也起到自动稀疏化的作用。

本文将传统算子与改进算子及其切线、伴随算子嵌入中国气象局数值预报中心研发的GRAPES_3Dvar同化系统(薛纪善等，2008)，添加了雷达径向风预处理模块，完善和优化了径向风同化模块。选取华北地区2014年6月6日降水个例，分析了采用不同观测算子同化径向风资料后动力场的改变以及对短时预报的影响。

1 技术方法简介 1.1 雷达径向风观测算子

雷达径向风与三维风矢量的几何关系为：

$V_R^* = u\cos \theta \cos \phi + v\sin \theta \cos \phi + w\sin \phi $

(1)

式中，V_R^*为原始雷达径向风，(u, v, w)为常规风场的三维分量，θ为雷达观测点的方位角，φ为雷达观测点所在仰角。

本文仅同化水平风场，因此将式(1) 中径向风垂直分量扣除，重新整理成为式(2) 形式：

${V_r} = u\cos \theta + v\sin \theta $

(2)

式中，V_r为处理后的水平径向风，其与原始雷达径向风关系如式(3) 所示：

${V_r} = \left( {V_r^* - w\sin \phi } \right)/\cos \phi $

(3)

本文将式(2) 的右端作为传统观测算子进行同化试验。通过式(2) 可以发现，传统径向风观测算子仅能引入单一的径向风增量，但是通过同化系统后会产生风场两维增量，即同化系统在吸收径向风信息后根据背景误差协方差矩阵调节得到风场的切向分量，也就是说，分析风场在径向上吸收了观测资料而体现出观测场的特征，而切向分量则无法反映观测场的切向特征。因而采用传统算子同化得到的分析风场是不准确的。

改进算子以IVAP反演方程组为基础。提取单位分析单元内的平均风信息，在同化系统背景场和物理条件的约束下，求取式(4) 左右两端极小化得到分析风场。

$\begin{array}{l} \frac{{\sum\limits_\mathit{\Omega } {{V_r}\cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\cos }^2}\theta } }} = \bar u + \bar v\frac{{\sum\limits_\mathit{\Omega } {sin\theta \cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\cos }^2}\theta } }}\\ \frac{{\sum\limits_\mathit{\Omega } {{V_r}\cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\sin }^2}\theta } }} = \bar u\frac{{\sum\limits_\mathit{\Omega } {sin\theta \cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\sin }^2}\theta } }}+\bar v \end{array}$

(4)

式中，u, v分别为分析单元Ω内的平均纬向风和经向风，本文中分析单元Ω为10 km×10 km的正方形区域，这样避免了分析单元随距离雷达中心的位置改变而发生形变的情况。

将式(4) 的左端作为同化变量，用Y₁, Y₂表示，代表径向风在纬向和经向上的投影信息

$\begin{array}{l} {Y_1} = \frac{{\sum\limits_\mathit{\Omega } {{V_r}\cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\cos }^2}\theta } }}\\ {Y_2} = \frac{{\sum\limits_\mathit{\Omega } {{V_r}\sin \theta } }}{{\sum\limits_\mathit{\Omega } {{{\sin }^2}\theta } }} \end{array}$

(5)

将式(4) 右端作为相应的改进观测算子，用H₁，H₂表示

$\begin{array}{l} {H_1} = \bar u + \bar v\frac{{\sum\limits_\mathit{\Omega } {sin\theta \cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\cos }^2}\theta } }}\\ {H_2} = \bar u\frac{{\sum\limits_\mathit{\Omega } {sin\theta \cos \theta } }}{{\sum\limits_\mathit{\Omega } {{{\sin }^2}\theta } }} + \bar v \end{array}$

(6)

通过式(4) 可以发现，同化变量由传统单一的径向风变量转变为径向风在经向和纬向上的投影项，在引入真实径向风的同时引入了部分切向信息，使得进入同化系统的信息是完整的风场观测信息，也就是说，引入径向风在经向和纬向投影上的两个信息后通过同化系统得到常规风场的两维增量，即两个方程(式4) 求解得到两个未知变量(常规风场的两维分量增量)满足了计算方程租的闭合问题。改进算子另一个优势则体现在分析单元的改变。由于雷达资料空间分辨率远高于模式分辨率，为了降低计算代价和避免观测相关性问题引入同化系统(赖安伟，2008)，稀疏化则成为同化前必须进行的步骤。传统算子的分析单元为格点，传统稀疏化(Super-Obbing，Liu et al，2005)方法会人为地去除大量的有效雷达观测信息以及引入大量的插值误差而污染原始观测。而改进算子采用人为可调节的区域作为分析单元，稀疏化的过程能够通过观测算子自身来实现。区域的分析单元一方面保证了大量原始观测数据能够被同化系统吸收，另一方面，同化变量的大小取决于分析单元内所有观测信息的累加，这样降低的单个观测信息或特殊位置零散观测对于同化系统的影响权重，保证了数据的连续性、准确性和稳定性。

改进算子不再单纯地引入传统意义的风场分量，而是将径向风在经向和纬向上的投影信息作为同化变量引入同化系统，利用同化系统的背景场和物理平衡关系作为弱约束条件，通过迭代方法得到全局最优风场。从本质上来看，采用改进算子的同化方法不再是一种强约束求解风场的过程，它利用同化系统的自身属性来得到最优风场，也就是说均匀风假设的强约束不一定完全满足，分析场仅需满足代价函数值极小即可，这样极大程度上削弱了由于引入局地假设条件而带来的反演误差。

1.2 雷达径向风同化系统检验

图 1为GRAPES_3Dvar雷达径向风同化流程图以及相应程序模块。在进行同化试验前，对GRAPES_ 3Dvar系统进行切线性检验与伴随检验(Kalney，2003)。

切线性模式检验判据：

$Ratio = \frac{{F\left( {{x_0} - \alpha \times {\rm{d}}x} \right) - F\left( {{x_0}} \right)}}{{\alpha \times F'\left( {{\rm{d}}x} \right)}}$

(7)

伴随模式检验判据：

$F'\left( {{\rm{d}}x} \right) \times F'\left( {{\rm{d}}x} \right) \approx {\rm{d}}x \times {F^{\rm{T}}}\left( {{\rm{d}}x} \right)$

(8)

式(7)、(8) 中F(x)为正演算子，F′(x)为切线性算子，F^T(x)为伴随算子。理论认为，若式(7) 中Ratio值趋近于1，同时式(8) 中左右两端计算结果趋于相等时，同化系统通过正确性检验。表 1为切线性模式的结果，通过表 1可以发现，随着α由0.1到1，检验结果得到的Ratio值稳定，均接近于1。表 2为式(8) 左右两端计算结果以及两者的差值。在双精度的情况下，式(8) 左右两端差值较小，表明同化系统的伴随模式通过了检验。以上数据皆表明了加入雷达径向风同化模块后的GRAPES_3Dvar同化系统通过了正确性检验，能够进行后续试验研究。

2 雷达径向风质量控制

由于雷达探测到的信息通常包含大量的气象噪音，使得雷达回波中存在大量的错误数据，因此雷达资料预处理决定了进入同化系统观测信息的质量。

地物杂波在径向风速图上的特征是低仰角扫描径向速度值接近零的区域内出现的非零观测，本文统一将径向风速度小于0.5 m·s^-1的观测点去除(俞小鼎，2006)。

为了保证进入同化系统的观测资料具有较好的连续性，因此剔除孤立观测值，若有效观测点周围8个点中有多于4个点为缺测值，将该点判定为缺测值。

中小尺度系统中水平风速特征值是垂直速度特征值的100倍，而垂直速度特征值与降水粒子末速度特征值量级相近(万齐林等，2005)。因此本文利用雨滴下落末速度替代径向风垂直分量后直接按照式(3) 扣除，仅进行水平径向风的同化试验。

传统稀疏化方法将原始径向风在水平方向粗化为与模式分辨率相同的观测点。本文中传统算子同化试验时采用的稀疏化方法是选取距离模式背景场格点最近的有效观测值作为该格点的观测信息。同时稀疏化过程中保留观测点的风速、方位角、距离库和经纬度等原始属性，在进行观测变量计算时采用观测点原始信息，不引入插值误差。

改进算子自身具有稀疏化功能，因此仅需将观测资料按分析区域大小进行分组。按照同化试验中确定的区域大小(本研究中选取10 km×10 km的正方形区域)将完整的观测场划分为若干个分析区域，将各个分析区域内的观测资料按照式(5) 计算后作为一组同化变量进入同化系统后将该分析区域内观测赋为缺省。一方面避免了观测资料多次利用的同时也达到了“稀疏化”的目的，另一方面，克服了雷达不同扫描半径上空间分辨率不同的问题。同时，为了避免分析单元内观测数量较少而引入较大误差，对分析单元设定阈值，这样保证了分析区域内有足够的观测资料用来计算同化变量，也保证了分析单元中观测信息的连续性。

图 2为2014年6月6日03时北京南郊S波段雷达2.5°(第三个仰角)的原始径向风观测与经过稀疏化与分组过程后的观测分布示意图。通过图 2可以发现，传统稀疏化处理得到的观测场(图 2b)的信息以单点形式存在，而经过改进分组得到的观测场(图 2c)信息则是以正方形的区域为代表，每个方框为一组同化变量，而方框中的所有有效观测信息均参与了相应分析区域的同化变量的计算，相比于传统稀疏化方法仅选用单个观测点作为同化变量，改进方案一方面能够吸收更多的观测信息，另一方面也降低了单个观测信息对于分析场的影响比重，确保了雷达完整观测信息的可靠性。

3 个例分析 3.1 资料与模式配置

本文选取2014年6月6日华北地区的一次降水过程进行对比分析。03时雷达强回波带(图 3a)恰好位于北京西北部边界区域，之后雷达回波不断向东移动，逐步影响北京城区以及河北东北部，因而选取03时作为积分初始时刻。雷达径向风观测资料来自北京南郊观象台S波段雷达(39.88°N、116.47°E)，采用雷达9个仰角的观测，径向速度的最大观测半径为230 km，为了有效抑制二程回波，本文仅选用150 km(600距离库)以内的观测(陈敏等，2014)。同化试验背景场采用美国NCEP中心GFS(Global Forecast System)的再分析资料，水平分辨率为1°×1°，垂直层次为26层。模式分析区域为(36.88°~42.88°N、113.47°~119.47°E)，模式中心位于39.88°N、116.47°E，模式水平方向格点数为201×201，水平格距为0.03°×0.03°，垂直方向为31层。设计以下三组试验，验证不同观测算子对于数值模拟结果的影响。

Control试验：采用GFS资料生成初始场与边界场，进行6 h积分预报。

Tradition试验：采用GFS资料生成背景场，利用传统算子同化雷达径向风后，积分6 h。

Advance试验：采用GFS资料生成背景场，利用改进算子同化雷达径向风后，积分6 h。

3.2 分析风场与散度场

图 3为模式第九层(高度约为2.5 km)背景风场与两种不同观测算子得到初始场的分析风场与散度场的合成图。图 3a中青色线条为组合反射率大于35 dBz的区域，即该区域是对流系统中心所在区域。从图 3a可以发现，分析区域受10 m·s^-1左右的均匀西南风气流控制，仅在41°N、115°E处存在少量的局地辐散信息(正值阴影)，而在对流较为旺盛的区域缺少与对流系统对应的中小尺度信息。加入雷达径向风资料后，两个试验的分析风场在北京西北部边界由之前均匀的西南风转变为一条西风与西南风的辐合带，风速也有所增加。同时，同化雷达资料之后的散度场相较于背景场的散度信息也有了大幅度的增加。Tradition试验结果(图 3b)表现为一条狭窄的辐合带(负值阴影)，辐合中心位于雷达中心东北部，而北京西北边界以外的辐散区域(正值阴影)较背景场增加幅度相对较小，范围也相对狭窄。Advance试验的散度特征(图 3c)在北京西北边界内部表现为一条量级相近的辐合带(负值阴影)，范围较Tradition试验结果有所增加，同时，北京西北边界以外的辐散场较Tradition试验结果增幅较大，接近带状分布。

由于对流系统前部较强的斜升气流与对流系统低层的冷池出流相互作用，因而会加强对流系统前部低层大气的辐合(陈明轩等，2012)。加入雷达径向风后，两种算子得到的散度场在北京西北边界内均表现为局地辐合而在对流系统后部表现为局地辐散。这些均表明加入的雷达径向风产生的局地动力信息是准确合理的。然而Tradition试验在雷达中心东北部产生了一个辐合中心，根据图 3a可以发现，对流系统在相应位置并无明显的回波中心，仅在附近存在一些零散的回波信息。而Advance试验散度场得到的散度场分布相对均匀，仅在40°N、116°E存在一个范围较小的辐合中心，与强回波带位置对应较好。

3.3 径向风场

为了更加清晰地说明不同算子对于风场的影响，以下将针对风场径向分量与切向分量进行对比分析。图 4为模式第九层观测、背景以及不同试验得到的径向风场的对比图。图 4a为观测径向风场，正径向风速中心位于北京中部而负径向风速中心位于北京以西区域，同时北京西北部边界存在一条约为5 m·s^-1的负径向风速带，向东北延伸至41°N，在北京西北部边界处造成正负径向风速的辐合区。图 4b为背景径向风场，表现为正负径向风速的均匀过渡，缺少大量的中小尺度动力场的辐合辐散的信息，同时正径向风中心位置偏北，负径向风中心量级偏小，位置偏南。加入雷达径向风资料后，不同试验结果得到的分析径向风场相比于背景径向风场均接近于观测场，在北京西北边界区域产生一条向东北方向延伸的负径向风速带，同时，观测场中的正径向风速中心的位置、量级也得到了校正。然而不同试验得到的分析径向风场依然存在差异，Tradition试验(图 4c)产生的负值中心的范围与量级相比与Advance试验(图 4d)的结果偏小，而后者更接近实况观测。由于两个观测算子中均引入了V_r项，因而得到的分析场均接近雷达观测。然而采用改进算子的同化系统能够吸收的径向风观测的数量远多于采用传统算子的同化系统，大量的观测信息使得改进算子得到的径向风场在极值中心量级以及范围上更接近于观测。

3.4 切向风场

由于无法获取对流系统的观测切向风信息，因而本文将IVAP反演(Liang，2007)得到的切向风场作为参考进行对比分析。图 5为模式第九层反演、背景以及不同试验得到的切向风场的对比图。图 5a为由IVAP方法反演得到的切向风场，整个观测区域内存在明显的正负切向风辐合，北京西北部边界以内为正的切向风区域且量级较大，而在参考场的西南部以及北部则存在切向风负值区域。图 5b为背景切向风场，整个区域内风速基本为10 m·s^-1，在41°N、115.5°E区域存在一正值中心，与IVAP反演得到的切向风场差异较大。Tradition试验的切向风场(图 5c)与背景切向风场的分布相近，仅削弱了位于41°N、115.5°E区域的切向风中心范围。采用传统算子得到的切向风增量是由背景误差协方差矩阵调节得到，因而切向风场相对于背景切向风场变化相对较小。而改进算子将径向风速随空间的变化信息(切向风信息)引入同化系统，因此Advance试验(图 5d)的切向风相对于背景切向风场变化较大，在北京西部形成了一条西南—东北走向的切向风正值带，同时分析风场的西南部以及北部的正值区域的量级有所削减，这都表明Advance试验的分析切向风场接近IVAP方法得到的反演切向风场。然而改进算子的切向风场并未出现反演风场(图 5a)中量级较大的极值中心，主要是由于IVAP方法在局地强约束下得到的风场存在较大的反演误差，而改进算子以背景误差协方差为弱约束，在求解切向风的同时一定程度上保留了背景切向风场的特征。

以上分析表明，在同化雷达径向风资料后，径向风场与切向风场相对于背景风场均产生了大量中小尺度信息。由于两个观测算子均引入了径向风观测项，因而均能有效准确地吸收雷达观测到的径向风信息，其中由于采用改进算子的同化系统利用了较多的雷达观测信息，因而得到的径向风场更为准确。而针对于切向风场，传统算子得到的切向分量与背景切向分量差别相对较小，而改进算子引入了径向风随空间的变化特征，相比于背景风场的切向分量表现出较大的改变，更接近IVAP反演得到的切向风场，同时在背景误差协方差约束下求解得到的切向风场也削减了由于局地强约束求解而带来的反演误差。

3.5 降水模拟

图 6为2014年6月6日03—09时华北地区自动站累计降水以及不同试验得到的累计降水模拟。从图 6a可以看出，降水主要集中在北京以及河北东北部地区，较大量级的降水雨带主要集中在北京西北部边界以及天津北部。Control试验结果(图 6b)较好地模拟出降水雨带的总体范围，对于北京西北部的降水中心模拟较为准确，然而在河北东北部则产生了一条西南—东北走势的较大量级的虚假降水雨带。加入雷达径向风资料后，Tradition试验(图 6c)对于北京东北部虚假雨带存在一定的削弱，但是对于北京西北部边界的降水雨带也产生了一定的削减，而Advance试验结果(图 6d)相比与Tradition试验结果，前者在较好模拟出北京西北部降水雨带范围、量级的同时，也极大地削弱了河北东北部虚假雨区的范围，同时对于北京南部大雨量级以上的降水有所改善，与实况观测较为接近。

对三组试验结果进行TS评分与Bias评分检验。图 7为三组不同试验6 h累计降水的TS评分与Bias评分。通过图 7a，可以看出加入雷达资料后，TS评分值相对于控制试验有了一定的增加，其中Advance试验的TS评分值相比与Tradition试验的结果更大，改进效果尤其体现在较大量级的降水上。结合Bias评分(图 7b)来看，所有量级降水的Bias评分值均接近于1，表明所有模拟试验的干湿偏差均在可接受范围内。加入雷达资料后，Bias评分均有所增加，其中以Advance试验结果最为明显。综合降水评分可以看出，加入雷达资料后的初始场在积分过程中有效地抑制了较高量级的虚假降水，其中以Advance试验的效果最好。

4 结论

本文对加入雷达资料预处理模块以及雷达径向风同化模块后的GRAPES_3Dvar系统进行了切线性模式检验与伴随模式检验，然后选取2014年6月6日华北地区的降水过程进行同化模拟试验，对比分析采用不同观测算子对初始场的动力调整以及对中小尺度降水的影响。

切线性模式与伴随模式检验结果均表明，加入雷达径向风同化模块后的GRAPES_3Dvar同化系统通过了正确性检验，能够进行雷达径向风变分同化试验。

从个例的试验结果来看，雷达径向风资料同化加强了中小尺度对流系统低层的风场辐合，同化后的分析场具有更为明显的中尺度信息，其中Advance试验所产生的散度场与雷达回波信息对应较好。从初始风场径向分量来看，两个试验均能有效地吸收观测的径向风信息，由于改进算子能够利用更多的观测信息，因而Advance试验结果在极值中心的范围与强度上与观测更为接近；从切向分量的对比来看，Tradition试验得到的分析风场的切向分量相对于背景切向风改变较小，而改进算子则充分吸收了切向信息，与背景切向风场差异较大，在更为接近IVAP反演切向风场特征的同时也保留了部分背景场的切向特征，即在保证吸收切向观测信息的同时也抑制了局地强约束导致的反演误差。因而，就动力场的调整而言，改进算子的同化结果优于传统算子的同化结果。

对不同观测算子得到降水模拟预报结果来看，加入雷达径向风的试验结果均提高了预报的准确率，在有效模拟出降水雨带范围与量级的同时抑制了控制试验中的大面积虚假降水，其中以Advance试验效果较好。而降水评分也验证了相应的结论。

雷达径向风同化对于短时强降水预报有着较大的贡献，本文基于均匀风假设引入IVAP反演方程组作为改进观测算子的基础，验证了改进算子相对于传统算子的优越性。然而不同的反演方法能够提供不同的风场信息，改进算子能够进一步扩展至线性风场假设甚至非线性风场假设，下一步研究会将风场假设进一步拓展，引入线性风假设，扩展观测算子进行相应的同化试验。