1 引言

多普勒雷达具有高时空分辨率的特点, 是观测风暴内部三维结构的唯一手段。我国新一代多普勒天气雷达观测网在实时监测及预警对流天气、降水及台风等方面起到了重要作用。多普勒雷达仅能提供反射率、径向速度、谱宽信息, 而为了更好的研究天气系统, 三维风场是非常必要的。多普勒天气雷达资料同化技术能够得到高分辨率、高精度的三维风场, 并且可以为数值预报提供高质量的初始场, 这种方法目前在短临预报及数值预报领域均被采用(张沛源等, 2002; 周海光等, 2002; 盛春岩等, 2006; 陈力强等, 2009; 覃月凤等, 2015; 徐琪等, 2015; 杨银等, 2015; 张桂莲等, 2018; 肖艳姣, 2018)。

雷达径向速度同化可以提供额外的风信息, 尤其在常规数据稀疏的情况下, 雷达径向风数据可以提供重要的信息, 显著改进初始场(王永明等, 2016; Li et al, 2008; 王晓峰等, 2017), 从而加强中尺度对流系统近地层的风场辐合情况, 更好地反映风暴结构(张林等, 2006; Ge et al, 2013; 罗义等, 2014; 马昊等, 2016; 薛谌彬等, 2017; 周荣卫等, 2018)。经过质量控制的雷达径向风、反射率因子资料经同化系统同化后, 可形成合理的分析增量, 使同化风场与实况更接近, 从而提高模拟降雨的准确性(张新忠等, 2015; 文秋实等, 2017)。

雷达资料同化的结果也可以直接应用于短临预警。Gao et al (2013)利用ARPS系统(Xue et al, 2000, 2001)的三维变分同化系统(Gao et al, 2004), 同化多部多普勒雷达形成每5 min、1 km分辨率的分析场, 进行中气旋识别, 该系统在2010年的恶劣天气试验中表现很好, 准确、精细的反映出风暴的强度及涡旋结构。美国全国雷达拼图和多传感器QPE系统(Zhang et al, 2011)及警告决策支持系统(Lakshmanan et al, 2007)都采用雷达资料同化的技术。VDRAS系统快速更新同化雷达资料及地面自动站资料以得到三维分析场(孙娟珍等, 2016), 该系统在2008年北京奥运会期间发挥重要作用(陈明轩等, 2011)。

雷达径向速度资料同化有间接同化和直接同化两种方法(杨毅等, 2007), 现在WRF、ARPS和GSI等主流数值模式大多采用直接同化的方法(Hu et al, 2006b; Li et al, 2012; Abhilash et al, 2012)。在雷达资料预处理方面不可避免的遇到坐标转换问题, 雷达资料呈锥面分布, 因此在对雷达数据进行处理时需要进行从极坐标到笛卡尔坐标的坐标转换。雷达分析的坐标一般可分为圆柱坐标系和笛卡尔坐标系。早期的多普勒雷达分析一般采用圆柱坐标系(Lhermitte et al, 1970; Miller et al, 1974)。进行双多普勒雷达风场分析时需要将雷达数据从球坐标插值到圆柱坐标系, 完成分析后再从圆柱坐标系插值到笛卡尔坐标系。为了提高垂直速度估计, Ziegler(1978)提出了一种变分双多普勒分析方法来处理上边界条件的问题。Sun et al(1997)利用伴随方法进行双雷达分析和物理量反演。与在圆柱坐标系中求解的那些方法相比, 这些采用笛卡尔坐标系的方法仅需要一次插值。

观测资料同化系统中直接同化雷达径向速度资料一般需要计算观测增量, 该过程首先需将雷达径向速度数据处理至模式网格, 其中有两种方法:一是将雷达径向速度资料通过最小二乘法从极坐标映射至模式标量的三维网格点(称为Grid方案); 二是通过双线性插值将雷达径向速度数据水平插值至的标量水平网格, 但垂直高度不进行插值, 保留在雷达径向坐标上(称为Tilt方案)。因为雷达数据在垂直方向分辨率较低, 将数据从仰角面插值到模式坐标会产生较大的误差, 直接同化未做垂直插值的PPI面上的数据可以在三维变化同化框架下进行最优插值(Sun et al, 2001), 在悉尼2000预测示范项目中就采用的Tilt方案(Crook et al, 2002)。江苏省气象局与美国强风暴实验室联合开发的高精度数值分析及预报系统(Precision Weather Analysis and Forecast System, PWAFS)中对雷达径向速度资料的预处理包括Grid方案和Tilt方案, 本文选取了华东地区几次典型天气过程, 用两种方案进行试验, 对比分析它们在不同类型个例中的表现, 对同化效果进行检验, 为数值预报及数值模拟方案的选择提供参考意义。文中地图包括江苏省界, 地图数据来自中国气象局Micaps系统, 底图未改动。

2 同化系统

本研究采用江苏省气象局与美国强风暴实验室联合开发的PWAFS系统, 该系统的同化系统采用ARPS系统的三维变分系统(Gao et al, 2004; Xue et al, 2003; Hu et al, 2006a, 2006b)。其代价函数写为(Gao et al, 2004):

$ \begin{aligned} J(x)=& \frac{1}{2}\left(x-x^{b}\right)^{T} B^{-1}\left(x-x^{b}\right) \\ &+\frac{1}{2}\left[H(x)-y^{o}\right]^{T} R^{-1}\left[H(x)-y^{o}\right]+J_{c}(x), \end{aligned} $

(1)

式中: x为分析因子, 包含三维风分量(u, v, w), 位温θ, 气压p, 水汽混合比q_v; x^b代表背景场; B是背景误差协方差矩阵; H(x)是前向算子, 它将模式变量转换到观测点坐标; R是观测误差协方差矩阵, 包括仪器误差和代表性误差; y^o为观测向量。代价函数包括三项, 右侧第一项是背景项, 第二项是观测项, 第三项J_c(x)为弱连续方程约束。

本文中y^o只包括雷达径向速度, 雷达前向观测算子可写为(Doviak et al, 1993):

$ v_{r}=\frac{\mathrm{d} h}{\mathrm{d} r} w+\frac{\mathrm{d} s}{\mathrm{d} r}(u \sin \phi+v \cos \phi), $

(2)

式中: v_r是雷达径向速度投影; r是斜距; h是高度(需考虑地球曲率半径); s是沿地球表面的距离; φ是雷达方位角。根据典型雷达设备误差(Miller et al, 2003), 观测误差设为1 m·s^-1, 则式(1)中第二项简化为:

$ J_{o}(x)=\frac{1}{2}\left[H\left(v_{r}\right)-v_{r}^{o b}\right]^{2}, $

(3)

利用式(2)将分析场数据投影到观测空间, 即通过一个三线性插值算子将背景场数据从模式格点映射到雷达观测点。

其中, 由于雷达风场只有径向观测, 为约束风场, 并产生合理的切向和垂直风场, 使用了滞弹性近似的连续方程约束, 从而建立了风场三个分量之间的关系, 并在一系列研究中取得较好的分析和预报结果(Gao et al, 2004; Hu et al, 2006b; Schenkman et al, 2011)。

$ J_{c}=\frac{1}{2} \lambda_{c}\left(\frac{\partial \overline{\rho} u}{\partial x}+\frac{\partial \overline{\rho} v}{\partial y}+\frac{\partial \overline{\rho} w}{\partial z}\right)^{2}, $

(4)

式中: $\overline{\rho} $是给定水平高度的平均空气密度; 权重系数λ_c控制弱约束项在代价函数中的重要程度, 本研究中λ_c=5. 0×10^-4。

3 两种方案描述

PWAFS系统的同化系统采用ARPS3DVAR, 预报系统采用ARW-WRF。WRF系统以GFS预报场为背景场冷启动4次, 输出9 km格距预报结果, 水平格点数240×240, 垂直方向27层, 覆盖我国东部及近海地区。微物理参数化方案采用Single-Moment 6-class, 长波辐射及短波辐射均采用RRTMG方案, 地表层采用改进的MM5 Monin-Obukhov方案, 陆面过程采用统一的Noah方案, 积云参数化方案采用Kain-Fritsch方案, 边界层采用YSU方案。该预报系统逐12 min输出预报结果, 将此预报结果从9 km格距的WRF网格插值到1 km格距的ARPS网格作为同化分析系统的背景场。同化分析系统水平格距1 km, 格点数720×720, 垂直方向50层, 覆盖江苏省及周边地区, 每12 min运行一次, 输出同化分析场数据。该系统可同化江苏及周边地区多普勒天气雷达资料、地面自动站资料及风廓线资料, 本研究中仅同化多普勒天气雷达资料。本文中对比的两种方案主要差别在雷达资料预处理中对雷达径向速度资料插值方案的选择, 两种方案对反射率资料的处理是一样的, 均是插值到模式三维网格后进行同化。在PWAFS系统中利用88d2arps(Brewster et al, 2005)模块对雷达资料进行数据质量控制、坐标转换等预处理。

Grid方案将雷达径向速度资料通过最小二乘法从极坐标映射到模式三维网格。这种方案选用目标格点相邻两层仰角、每层仰角上目标格点周围小范围内的数据进行计算。首先, 径向速度从平均风转换为增量, 其中平均风廓线由雷达周围的背景风场中的数据点的平均值提供。其次, 通过检测距离库之间切变的方法进行水平风一致性检查(Ellis et al, 1990)。最后, 利用拟合插值法将数据从极坐标重新映射到模式的笛卡尔地形伴随坐标, 使雷达数据有相同的空间分辨率。这是通过最小二乘拟合到局部多项式函数来实现的:

$ A=a_{0}+a_{1} x+a_{2} x^{2}+a_{3} y+a_{4} y^{2}+a_{5} x y+a_{6} z, $

(5)

可见其在水平方向是二次函数, 垂直方向是线性的。其中: A是分析变量; a_i是多项式系数。(x, y, z)是模式坐标值偏差。为了避免不合理的外推, 结果被限制在符合最小二乘计算的数据范围内。当模式网格分辨率为几公里时, 这个算法在近雷达处进行平滑, 在远雷达处进行插值, 使数据分布在水平方向上更均匀, 以避免在近距离范围内同化过多数据带来的计算成本。

Tilt方案中通过双线性插值将雷达径向速度数据水平映射至的标量水平网格, 但垂直高度不进行插值, 保留在雷达径向坐标上。该方案仅用目标格点对应水平位置两侧相邻径向上、前后两个库的四个数据进行双线性插值, 得到的数据保留在雷达仰角面上。由于模式分辨率(一般为几公里)往往较雷达径向速度分辨率(250 m)低, 该方案减少了雷达径向风观测垂直插值引发的误差, 并在一定程度上起到了数据稀薄的作用(Tong et al, 2008; Dong et al, 2013)。

为验证方案的不同效果, 本文选取华东地区几次典型个例进行敏感性试验, 并对低层(< 1000 m)流场进行详细分析。选取的个例包括2016年6月23日江苏阜宁龙卷个例, 2017年3月1日江苏阜宁大风个例及2017年6月10日长江下游地区梅雨个例。

4 试验结果分析和对比 4.1 2016年6月23日阜宁龙卷个例

2016年6月23日, 江苏省盐城市阜宁县和射阳县发生由龙卷导致的特大灾害, 导致99人死亡、800多人受伤, 大量基础设施损毁, 基于详细的现场调查资料, 评估江苏阜宁龙卷为EF4级(郑永光等, 2016)。利用当天盐城雷达资料分析发现23日14 : 00(北京时, 下同)—15 : 00中气旋影响区域及其西侧共计有5个自动站的瞬时风速超过8级, 大风范围非常小, 仅出现在阜宁县西南部长25 km、宽10 km的范围内, 最大风速在阜宁县新沟镇为34. 6 m·s^-1 (12级以上, 时间为14 : 29)(张小玲等, 2016)。新沟镇距离盐城雷达52 km, 雷达0. 5°仰角在此处高度700 m。

2016年6月23日14 : 24盐城雷达观测到的反射率因子及径向速度图(图 1)显示, 新沟镇位于龙卷勾状回波位置, 从径向速度图上显示明显的正负速度对, 表示此处有很强的气旋, 这与典型的龙卷回波结构特征一致。400 m及1000 m高度上背景场及两种方案的同化结果(图 2)表明, 400 m高度低于雷达最低探测仰角在新沟镇的探测高度, 因此在400 m高度、龙卷发生位置附近是没有雷达探测数据的。从图 2(a)可以看出, 背景场中无龙卷的勾状回波特征和涡旋结构, 风场较均匀。而雷达资料同化之后可以清楚的看到龙卷的勾状回波[图 2(b), (c)]。径向速度预处理的两种方案在风场同化结果上表现出明显的不同。图 2(b)为Grid方案得到的同化结果, 这种方案在龙卷勾状回波西侧区域得到的水平风速超过12 m·s^-1, 在勾状回波南侧同化得到的水平风速超过18 m·s^-1, 在靠近雷达中心西边方向得到了超过16 m·s^-1的水平风速区域, 在龙卷北部的大部分区域也得到超过14 m·s^-1的水平风速。结果显示由于垂直插值的影响, 在部分没有雷达观测资料或雷达观测到径向速度较小的地方也得到了较大的同化风场, 如在龙卷西侧, 雷达观测到径向速度仅5 m·s^-1, 而同化风场超过12 m·s^-1。图 2(c)为Tilt方案得到的同化结果。对比图 2(b)可以看出, 该方案在龙卷勾状回波东南侧得到的风场结构及水平风速与Grid方案几乎一致, 在整个区域内水平风速超过10 m·s^-1的范围及最大水平风速均小于Grid方案, 该方案的结果与雷达观测径向速度结果较为一致。1000 m高度同化结果[图 2(d)~(f)]也可以看到相似的情况, 经过雷达资料同化后, 反射率及风场都比背景场有明显的调整, Grid方案对水平风速的调整范围比Tilt方案大, 在新沟镇西侧得到一个在雷达径向速度观测中未出现的大风速区。通过对比分析可知, Tilt方案得到更为贴近观测、合理的分析结果。图 3为400 m高度上两种方案同化结果与背景场的差, 亦即分析增量。从图 3可以看出, Grid方案分析增量区明显大于Tilt方案, Grid方案在超出雷达回波区域也得到很多分析增量, 而Tilt方案的分析增量大部分局限于回波覆盖范围内。1000 m高度上也有相同的现象(图略)。

图 4为Tilt方案的回波强度及风场结果。图 4(a)显示低层(400 m)有明显的强回波中心、钩状回波、涡旋中心及入流。龙卷涡旋中心位于钩状回波尾部, 强回波中心位于涡旋中心东北侧, 最大回波强度超过60 dBz, 钩状回波位于涡旋中心东南侧, 钩状回波东侧是很强的水平入流, 最大风速超过18 m·s^-1。在龙卷涡旋中心沿西南-东北方向线段及西北-东南方向线段[图 4(a)]进行垂直剖面分别得到图 4(b)及图 4(c)。图 4(b)中低层有明显的辐合, 后向气流在强回波中心后侧下沉, 前侧入流一直延伸到强回波后侧。强回波后侧低层有较弱的上升运动, 强回波整体呈强的下沉, 最大下沉速度超过10 m·s^-1, 7 km以上区域有上升运动。图 4(c)显示从龙卷回波前侧(钩状回波一侧)入流区到后部, 低层气流呈现下沉-上升-下沉的分布特征, 即在对流系统前侧入流区下沉, 强回波前部上升、后部下沉, 强回波后部的下沉区在地面产生了明显的辐散性流场。该结果清楚展示了龙卷发生位置的回波及流场的精细结构, 这些分析结果与经典龙卷分析结果一致(Lemon et al, 1979), 证明了方案的合理性。

4.2 2017年3月1日盐城大风个例

2017年3月1日16 : 00左右, 江苏省盐城市阜宁县境内因受较强冷空气影响, 突发雷雨大风天气, 造成人员伤亡和财产损失。根据阜宁气象观测站控测到的数据显示, 瞬间极大风速达到26. 9 m·s^-1(10级), 另根据区域自动观测站(三灶镇)探测数据显示为30. 7 m·s^-1(11级)。根据当日气象观测记录, 1日下午江苏省大部分地区出现7级以上大风, 截止17 : 00, 全省有26个基本站超8级、6个站超过9级、最大盐城阜宁达10级。

从2017年3月1日16 : 46分盐城雷达观测到的0. 5°仰角雷达反射率及径向速度可以看到, 雷达站西北侧有明显的雷暴区域[图 5(a)], 雷暴主体及其前侧有很强的朝向雷达中心的负速度区, 最大径向速度超过20 m·s^-1[图 5(b)], 表明从雷暴主体到前侧有很强的低层出流。分析认为这次风灾是由雷暴产生的强出流所致。

从两种同化方案得到的400 m及1000 m高度雷达反射率、水平风场及水平风速(图 6)可以看出, 两种方案得到的风场结构较一致, 水平风速量级及分布也较一致。Grid方案在雷暴区得到的最大风速为26. 96 m·s^-1, Tilt方案在雷暴区得到的最大水平风速为29. 64 m·s^-1。两种方案得到的结果与观测相比都很吻合, 对比来看Tilt方案的结果更接近观测最大风速值。Grid方案中雷达中心东侧分析得到的风速超过24 m·s^-1, 对比雷达观测径向速度图, 此处得到的同化速度偏强。

4.3 2017年6月10日梅雨降水个例

2017年6月10日, 受江淮气旋影响, 江苏的沿江苏南地区出现大暴雨, 南京站累计日降水量超过230 mm, 破历史纪录。南京市气象台2017年6月10日10 : 11升级暴雨黄色预警信号为暴雨橙色预警信号。全省多地区伴有强雷电、短时强降水和雷暴大风等强对流天气, 降雨最大的是常州金坛和镇江句容, 24 h累计雨量超250 mm, 达特大暴雨级别。

2017年6月10日10 : 00南京及常州雷达0. 5°仰角的反射率及径向速度图(图 7)显示, 南京及周边地区被降水回波覆盖, 在雷达站的西南方向吹西南风, 西南侧回波边缘处最大径向速度超过20 m·s^-1, 雷达站连线方向吹偏东风, 最大径向速度超过10 m·s^-1。图 8为两种方案在400 m及1000 m高度上得到的同化结果。两种方案得到的风场结构相似, 呈现气旋性结构的流场, 400 m高度西南侧吹偏南风, 雷达站连线吹偏东风, 1000 m高度南侧吹西南风, 雷达站连线吹偏东风。在400 m高度上, Grid方案[图 8(b)]未能较好的分析出回波西南边界附近的大风速区, 对回波东北部的径向速度超过10 m·s^-1的大片区域也未能分析出来, 在南京雷达站西北侧的风速区有较好的表现。Tilt方案[图 8(c)]在回波西南边界及东北侧都得到明显的水平风速大值, 最大风速23. 68 m·s^-1, 这种结果比Grid方案[图 8(b)]更符合观测结果。1000 m高度上, Tilt方案得到的水平风速大值区范围明显优于Grid方案, 两者最大水平风速差别不大。通过对比雷达观测, 该个例中Tilt方案得到的结果优于Grid方案。

5 结论与讨论

利用龙卷、强对流及梅雨天气个例对比分析两种雷达径向速度坐标转换方案对同化结果的影响, 得到如下主要结论:

(1) Grid方案将雷达反射率及径向速度资料从极坐标映射到模式三维网格, 该方案计算简便, 结果可靠, 并已被长期应用于三维变分同化系统中; Tilt方案中对反射率资料的处理与Grid方案一致, 仅将径向速度资料水平插值到模式水平网格, 但垂直方向不插值, 尽可能多的将径向速度资料保留在雷达径向上, 该方案更多的保留了观测信息, 避免插值引起的平滑。这种方案近年来被多个业务预报系统采用。

(2) 针对不同个例同化分析结果的检验发现, 2016年6月23日阜宁龙卷个例中, Grid方案在部分没有雷达观测资料或雷达观测到径向速度较小的地方也得到了较大的同化风场。Tilt方案结果清楚展示了龙卷发生位置的回波及流场的精细结构, 这些分析结果与经典龙卷分析结果一致。2017年3月1日阜宁大风个例中, 两种方案得到的最大风速值差3 m·s^-1, Tilt方案的结果更接近观测最大风速值。Grid方案在雷达中心东侧分析得到的风速比雷达观测偏强。2017年6月10日梅雨锋暴雨个例中, Grid方案未能较好的分析出回波西南边界附近的大风速区, 对回波东北部的径向速度超过10 m·s^-1的大片区域也未能分析出来, Tilt方案得到的水平风速大值区范围明显优于Grid方案。

低层大气结构对强天气触发及发展的作用非常重要, 因此获取低层信息对强天气的分析及研究很有帮助。Crook et al(2002)论述, 当近地面层的水平方向上相邻模式格点分别位于雷达最低层仰角上下方时, 如果两个格点处有显著的垂直风切变, 将导致同化产生异常的水平风切变。虽然可以通过改进观测及背景场的误差减轻这种现象, 但这是坐标转换的自然结果, 无法消除。Grid方案中得到的“最低点”高于雷达探测最低仰角, 而Tilt方案的“最低点”直接在最低仰角面上。在靠近雷达中心的低层, 观测资料密集, Tilt方案能够更好的反应实际大气状态。从这个角度来说, Tilt方案优于Grid方案。本文研究区域内雷达密集, 低层观测资料丰富, 采用Tilt方案有更好的优势。对比三个个例的同化分析结果, Tilt方案得到的水平风速分布比Grid方案更接近观测, 但是因为缺乏其他观测资料进行验证, 两种方案的效果还需要利用数值预报或其他方法进行对比。