气象学报  2014, Vol. 72 Issue (1): 133-151 PDF    
http://dx.doi.org/10.11676/qxxb2014.012
中国气象学会主办。
0

文章信息

秦琰琰, 龚建东, 李泽椿, 盛日锋. 2014.
QIN Yanyan, GONG Jiandong, LI Zechun, SHENG Rifeng. 2014.
集合方根滤波同化多普勒雷达资料在一次飑线过程中的应用研究
Assimilation of Doppler radar observations with an ensemble square root filter:A squall line case study
气象学报, 72(1): 133-151
Acta Meteorologica Sinica, 72(1): 133-151.
http://dx.doi.org/10.11676/qxxb2014.012

文章历史

收稿日期:2013-4-23
改回日期:2013-10-24
集合方根滤波同化多普勒雷达资料在一次飑线过程中的应用研究
秦琰琰1,2, 龚建东3, 李泽椿4, 盛日锋5     
1. 国家海洋环境预报中心, 北京, 100081;
2. 国家海洋局海洋灾害预报技术研究重点实验室, 北京, 100081;
3. 中国气象局数值预报中心, 北京, 100081;
4. 国家气象中心, 北京, 100081;
5. 山东省人民政府人工影响天气办公室, 济南, 250031
摘要:针对2005年7月12日发生在山东省中西部地区的一次飑线天气过程,采用集合方根滤波方法开展基于WRF模式的多普勒雷达资料的同化应用试验,考察了此同化系统对实际雷达资料的同化效果。主要结论如下:(1)集合方根滤波同化系统能有效同化实际雷达资料,雷达资料的加入增加了模式的中小尺度信息,使分析场得到了显著改善,有效缩短了模式起转时间,改进了对地面降水的预报。(2)利用三次同化分析后的集合平均分析场进行的确定性预报表明,与控制试验相比,同化后分析场能更准确地预报飑线系统的微物理量场,预报的流场结构符合风暴的动力特征,动力场和热力场的分布与配置也基本合理。(3)集合平均分析场对飑线系统传播方向的预报与实况一致,但预报的系统传播速度较实况快,由于对流系统的非线性发展迅速,对系统的预报时效为5-6 h。
关键词雷达资料同化     集合方根滤波     飑线    
Assimilation of Doppler radar observations with an ensemble square root filter:A squall line case study
QIN Yanyan1,2, GONG Jiandong3, LI Zechun4, SHENG Rifeng5    
1. National Marine Environment Forecast Center, Beijing 100081, China;
2. Key Laboratory of Research on Marine Hazards Forecasting, State Oceanic Administration, Beijing 100081, China;
3. Center for Numerical Prediction, CMA, Beijing 100081, China;
4. National Meteorological Center, Beijing 100081, China;
5. Weather Modification Office of People's Government of Shandong, Jinan 250031, China
Abstract:The effectiveness of using an Ensemble Square Root Filter (EnSRF) to assimilate real Doppler radar observations is investigated by applying the technique to a case of squall line occurred on 12 July 2005 in the midwestern part of Shandong Province as done by the Weather Research and Forecasting (WRF) model. The experimental results show that, (1) the EnSRF system has the ability to effectively assimilate Doppler radar data: since the information at convective scale is added into the numerical model through radar data assimilation the analyzed field is improved noticeably, the model spin-up time is shortened, and the surface precipitation forecast is improved correspondingly; (2) compared with the forecast of control run, the forecast initiated from the EnSRF analysis can give a more accurate prediction of the microphysical field, and the predicted wind field and thermal field are reasonable and they are in accordance with the characteristic of convective storm; (3) the propagation direction of the squall line forecasted by the ensemble mean analysis is consistent with observations, but the propagation speed is faster than observed speed, because of the nonlinear development of the convective storm, the period of validity of forecast is about 5 to 6 hours.
Key words: Radar data assimilation     Ensemble Square Root Filter     Squall line    

1 引言

数值预报是一个初边值问题,一个尽可能准确、协调的初始场对提高数值预报准确性有至关重要的作用。多普勒天气雷达探测资料具有高时空分辨率的特点,用其改进对流尺度数值预报的初始场,对提高对流尺度灾害性天气预警预报的准确率有重要意义。随着新一代多普勒天气雷达网的布设完善,天气雷达在探测能力、资料覆盖率、资料质量上均有极大的提高,为雷达探测资料同化提供了必要条件。近年来资料同化技术的快速发展和计算机水平的不断提高,也使得在高分辨率非静力模式中同化多普勒雷达探测资料成为可能。

目前,同化雷达资料的方法主要有反演类算法(Gal-Chen,1978Hane,et al,1985)、热动力初始化方法(Lin,et al,1993Xue,et al,1998Zhang,et al,1998)、松弛逼近方法(Jones,et al,1997)、变分类方法(杨毅等,2008Sun,et al,1997许小永等,2004杨艳蓉等,2008),以及集合卡尔曼滤波方法(EnKF)。集合卡尔曼滤波方法采用蒙特卡罗思路,通过一组短期预报集合来估计背景场误差协方差,自Evensen(1994)首次引入到地球物理领域以来,集合卡尔曼滤波方法以其所依赖的背景场误差协方差、模式和观测算子可为非线性、不需要切线性伴随、易于并行等诸多优点,已在多种大气和海洋模式中得到了广泛应用,成为目前最具业务发展前景的同化方法之一。近年来,在对流尺度系统中应用集合卡尔曼滤波方法同化雷达资料的研究取得了较大进展,许多学者通过同化模拟的雷达资料的观测系统模拟试验,得到了较准确的风暴流场、热力场和各种微物理量场,初步验证了集合卡尔曼滤波方法在对流尺度的同化效果,显示其在对流尺度应用的可行性和有效性(Snyder,et al,2003Zhang,et al,2004Tong,et al,2005Xue,et al,2006许小永等,2006秦琰琰等,2012),然而对于实际雷达资料在对流尺度的同化应用研究在中国尚少,而在国际上已有相当多的文献发表。

Dowell等(2004)最早将集合卡尔曼滤波方法用于对流尺度同化实际雷达资料,使用Sun等(1997)建立的三维暖云模式和基于集合方根滤波(EnSRF)的集合卡尔曼滤波同化系统,对一次龙卷超级单体个例,尝试了真实的单多普勒雷达反射率和径向风资料的同化,试验中同时同化了单多普勒雷达10个连续体扫47 min的反射率和径向风观测,较好地反演了龙卷超级单体风场结构,但由于模式误差和观测的限制,没能成功反演风暴低层冷池。许小永等(2006)使用相同的暖云模式,同时同化了宜昌和荆州双多普勒雷达6个连续体扫的径向风和由反射率导出的雨水观测,结果显示,集合卡尔曼滤波能准确分析暴雨中尺度对流系统的风场结构,反演的热力场和微物理场也基本合理。上述研究使用的均为暖云模式,不包含冰相微物理过程,同化分析中只对雨水混合比进行了更新,不包含对冰相水凝物的反演和分析,而冰相水凝物在发展旺盛的对流单体中是实际存在的,并且,冰相过程在其中所起的作用具有一定的重要性(胡志晋等,1998)。此外,其研究都只关注了集合卡尔曼滤波同化分析对风场、热力场和微物理量场的反演能力,并未关注利用同化后分析场进行确定性预报的效果。

本研究使用包含复杂冰相微物理过程的WRF模式,针对2005年7月12日发生在山东省中西部的一次飑线天气过程,采用集合方根滤波同化方法开展了实际雷达资料的同化试验,并利用三次同化循环后的集合平均分析场进行确定性预报,考察了集合方根滤波同化系统对实际雷达资料的同化效果。同化分析中对包含冰相水凝物在内的所有变量都进行了更新,不但关注同化时段内集合方根滤波分析对大气状态的反演能力,更考察了雷达资料同化对数值模式的初始化能力及对热力场、流场、水凝物场以及降水的预报能力。 2 集合方根滤波方法及同化系统介绍

传统的集合卡尔曼滤波算法在实施时需要扰动观测(Houtekamer,et al,1998),这会产生由观测扰动带来的样本误差,为此Whitaker等(2002)提出了不扰动观测的集合方根滤波方法,并证明对于洛伦茨模式和两层大气环流模式,即使在相同集合数的情况下,集合方根滤波也能得到比集合卡尔曼滤波更低的分析场集合平均误差。

与其他集合卡尔曼滤波算法一样,集合方根滤波算法也包含预报积分和分析更新两个步骤。在预报步骤实施一组短期集合预报,在分析步骤利用这组短期集合预报来估计背景场误差协方差,进而对预报场进行分析更新;更新后的分析场接着积分到下一个有观测的时刻,再进行分析更新,这样不断循环。与其他集合卡尔曼滤波算法不同的是,集合方根滤波将对预报场的更新分为对预报集合平均场 a进行更新和对相对于集合平均的扰动场x′a进行更新两部分,在观测误差独立且与背景场误差不相关的假设下,同化观测的先后顺序不影响同化结果,观测被逐个顺序同化,用同化了这个观测后的平均分析场和扰动分析场作为下个观测同化前的背景场,直到同化完所有的观测。

集合平均场的更新方程为

式中,上标a和f分别表示分析场和预报场,yo为同化的观测,H为观测算子,包含空间变换和物理转换,将变量从模式空间转换到观测空间,成为空间位置和要素类型都与观测对应的观测相当量。 K 为卡尔曼增益矩阵式中,R 为观测误差协方差,H 为H的雅可比矩阵,P f为从预报集合估计的背景场误差协方差,在集合方根滤波算法的实施中不需要直接将其计算出来,而是计算 P f H T和HP f H T 式中,下标k为集合成员序号,N为集合成员数。

由有限集合样本数引入的样本误差会导致背景场误差协方差的估计误差,尤其会在观测较远处产生虚假的分析增量。本研究采用Houtekamer等(2001)提出的Schur乘积方法解决这个问题,即观测空间的协方差局地化。如式(3)所示,ρ表示对背景场误差协方差矩阵做Schur乘积,ρ为准高斯型的局地相关系数,在观测点上ρ=1,而在观测的影响范围之外ρ=0。采用Gaspari等(1999)的五阶分段有理函数来定义局地相关函数ρ,五阶分段有理函数的形式为

式中,r为观测与背景场上分析格点之间的距离,c为给定的长度尺度参数。

集合方根滤波方法中,观测被依次逐个同化,这时 HP f H T和R 都将成为标量,这样对(HP f H T+ R)-1的计算就成为求标量的倒数。顺序同化观测技术和协方差局地化技术的使用,使得求算 K 过程中矩阵求逆的巨大计算量问题得以较好的解决,使集合方根滤波方法在实际大气模式的同化应用中可行。

相对于集合平均场的第k个成员的扰动场的更新公式为

式中,γ为协方差扩大因子,用以调节集合成员的发散度,从而缓解由于有限集合数带来的样本误差而造成的滤波发散问题。集合的离散度应该代表集合平均场的误差统计特征,然而有限集合数带来的样本误差会造成集合离散度与集合平均场误差的不一致,离散度过小会导致集合卡尔曼滤波分析不吸收观测,而离散度过大则会使分析过度偏离背景场。实施中,在同化第一个观测前在每个成员相对于集合平均的偏差场上乘以此系数。

同化的观测包括雷达径向风和反射率因子,雷达径向风和反射率因子的观测算子与Tong等(2005)中相同。雷达径向风的观测算子为

式中,Vr为雷达径向风观测,α为雷达波束的仰角,β为雷达波束的方位角,u、v、w为模式输出的风分量。

雷达反射率因子的观测算子为

式中,Z是反射率因子Ze的对数表示,其单位分别为dBz和mm6/m3

总的雷达反射率因子由雨、雪、霰的反射率因子Zer、Zes、Zeh相加得到

式中,各水凝物反射率因子的计算基于Smith等(1975)Lin等(1983)的算法求得,雨滴的反射率因子为

温度低于0°时,干雪的反射率因子为

温度高于0°时,湿雪的反射率因子为

霰的反射率因子

式中,ρ为空气密度;ρr=1000 kg/m3,为雨水密度;ρs=100 kg/m3,为雪的密度;ρi=917 kg/m3,为冰晶的密度;ρh=913 kg/m3,为霰的密度;qrqsqh分别为雨、雪、霰的比含水量;Nr=8.0×10-6 m-4Ns=3.0×10-6 m-4Nh=4.0×10-6 m-4分别为雨、雪、霰的Marshall-Palmer指数型滴谱分布的截断参数,K2i=0.176、K2r=0.93分别为冰晶和雨滴的介电系数。

秦琰琰等(2012)通过观测系统模拟试验,对此同化系统的可行性进行了验证。本研究采用此系统进行实际雷达资料的同化试验研究,验证其在实际天气过程中的同化效果。 3 个例介绍 3.1 天气概况

2005年7月12日,受东北低涡横槽转竖的影响,山东省中西部出现了一次以雷雨大风和冰雹为主的飑线过程。11时(北京时,下同)至15时,高唐、禹城、齐河、济南、莱芜、博山、沂源等地先后出现了雷雨大风或冰雹。降雹地点主要出现在高唐、禹城、莱芜、沂源等地,其中,莱芜市辛庄镇的冰雹直径达到25 mm,地面积雹最多处达30-40 mm;大风区位于济南至沂源之间,阵风风力都在8级以上,13时30分莱芜最大风速为22.4 m/s,14时40分沂源瞬时最大风速为29.2 m/s(11级)。这次过程给受灾地区带来了巨大的经济损失,仅沂源县的直接经济损失即达2.55亿元

淄博市气象局.2005.淄博市2005年7月12日雷雨大风冰雹灾害评估报告 3.2 形势场分析

2005年7月8-11日,东北地区一直维持一个低涡,在该低涡后部,西风带高压脊逐渐加强并向东北方向发展,低涡中心南移且减弱,低涡中心西部有横槽生成。12日08时,低涡中心减弱消失,横槽将转竖,山东500(图 1a)和700 hPa(图略)处于西北偏西气流控制之下,850 hPa(图 1b)及以下对流层低层,在河北、河南与山东交界处有一气旋式的风向切变线,山东的中西部地区为弱的西南风。分析济南章丘探空站变温发现,12日08时850 hPa 24 h变温为+2℃,500和700 hPa为负变温(-1℃),表明对流层低层增温明显,对流层中高层有冷空气扩散侵入。在这种高低空配置形势下,当横槽转竖时高层冷空气向下扩散,可形成上干冷下暖湿的不稳定层结,从而导致强对流天气的产生。

图 1 2005年7月12日08时500 hPa(a)和850 hPa(b)形势场 Fig. 1 Synoptic chart of 500 hPa(a) and 850 hPa(b)at 08:00 BT 12 July 2005
3.3 降水过程分析

本次飑线过程自西向东影响山东中西部地区,降水具有自西向东不断移动发展的特征(图略)。7 月12日08时在鲁西北地区开始出现降水,以雷雨或阵雨为主,量级均在1 mm以下;12-17时降水强度和范围逐渐增强增大,局地降水增大明显。12时雨强达18 mm/h,14时降水范围继续扩大且强度达到最强的43 mm/h,中心位于泰安附近,16时,降水范围继续东移,中心位于沂源附近,雨强为12 mm/h。18时,降水强度逐渐减弱,为一般性的雷阵雨,降水区分布在地面气压的高值区,18-20时,降水强度减弱、范围变小,系统逐渐移出山东。

3.4 雷达观测分析

位于济南齐河的CINRAD/SA多普勒天气雷达(36.81°N,116.78°E,雷达天线海拔高度72.9 m)捕捉到了这次中尺度系统的发展演变过程。该雷达采用VCP21观测模式进行连续体扫观测,每5-6 min一次完整体扫,共9个仰角,仰角为0.5°-19.5°,每个仰角体扫观测的反射率和径向速度的最大探测半径分别为460和230 km,分辨率分别为1和0.25 km。

从组合反射率因子演变可知,12日凌晨以后,在山西阳泉附近有对流回波生成,随着对流回波的东移生消、发展,于10时56分在德州夏津(回波中心强度55 dBz,回波顶高12 km)和河北清河附近各发展成两短对流带,这两处对流回波不断加强,逐渐首尾相连为带状(图 2a)。11时26分高唐附近出现强回波(图 2a)并逐渐南移加强,于12时02分在其西端出现弯曲的强中心(图 2b),回波中心强度达68 dBz,弓状回波开始发展。同时,在回波南面新生回波的后部形成明显的出流边界,出流边界的后侧不断有对流回波产生,于12时32分新生的对流回波与之前的回波带基本连在一起,形成 不连续的弓状回波带(图 2c),之后在弓状回波头颈部的后方出现入流槽口,意味着弓状回波的后部有强大的后侧入流。11-13时,弓状回波经过的高唐、禹城出现了冰雹;聊城、济南出现雷雨大风,济南的阵风达18 m/s。弓状回波不断向东南方向移动,并与其前方新的对流单体不断合并,13时57分之后,其北部逐渐减弱,回波强度降低。12时38分莱芜附近开始有孤立的块状回波单体生成,该回波在向东北移动的过程中发展加强,13时21分达到63 dBz(图 2d),此后逐渐减弱并于16时09分消亡。同时在该对流风暴的右侧又有新对流单体发展,于13时39分发展成超级单体,最大回波强度为67 dBz,回波顶高为13 km。受超级单体影响,莱芜出现冰雹和雷雨大风。14时03-21分,弓状回波演变成逗点云系,此后回波向山东东南方向移动,强度逐渐减弱。

图 2 2005年7月12日济南雷达的组合反射率因子(a.11时26分,b.12时02分,c.12时32分,d.13时21分) Fig. 2 Combination reflectivity from the Jinan radar at 11:26(a),12: 02(b),12:32(c) and 13:21 BT(d)12 July 2005
4 同化试验设计 4.1 试验方案

图 3为试验方案的示意图,包括实际雷达资料同化试验和控制试验两部分。控制试验只模拟不同化,采用NCEP再分析资料提供初始场和侧边界,从2005年7月12日08时开始模拟,积分12 h,模拟时段包含了这次飑线天气的整个发展过程;同化试验也同样自7月12日08时开始,在第1次同化分析前先做1 h集合预报,以初步发展一个对应于此次天气过程的、随背景场流型演变的背景场误差协方差结构,09时进行第1次同化分析,之后每30 min进行一次同化分析,如此循环,共进行3次同化分析,然后用3次同化后的集合平均分析场做10 h确定性预报。

图 3 同化试验方案流程 Fig. 3 Flow chart of the assimilation experiments scheme
4.2 预报模式及其设置

本试验采用完全可压缩欧拉非静力中尺度WRF v3.2模式,包含6个预报方程和一个静力诊断方程,预报变量为三维风场、扰动位温、扰动位势、干空气地面扰动气压,以及水汽、云水、雨水、冰晶、雪和霰的混合比等可选择的预报变量。采用两重双向嵌套,外层区域格点数为201×201,格距9 km,内层区域格点数为280×247,格距3 km,模式垂直28层,模式顶层为50 hPa,外层区域积分步长60 s。外层区域采用Kain-Fritsch(Kain,et al,19901993)积云对流方案,内层嵌套区域关闭积云对流方案,显式云方案采用WSM6方案(Hong,et al,2004),输出水汽含量和云水、雨水、冰晶、雪、霰共5种水凝物的混合比,初始场和侧边界均采用每6 h一次1°×1°的NCEP再分析资料。

4.3 同化系统设置

本试验同时同化了雷达反射率因子和径向风观测,并对所有变量进行了分析更新,但只更新内层区域。协方差局地化半径水平方向取为10个格距,垂直方向半径为5个格距;协方差扩大因子取2.0;观测算子为式(9)-(15);受计算机资源限制,集合成员取20个,初始集合由初始扰动场叠加NCEP再分析场而形成。基于理想试验中同化3个循环后就能较好重建模式风暴的结果(秦琰琰等,2012),本同化试验进行了3次同化分析。

4.4 初始扰动场生成方法

初始扰动场由WRF模式三维变分同化系统的cv3背景场误差协方差生成,与Zhang等(2006)用MM5模式三维变分同化系统来产生集合卡尔曼滤波初始集合扰动场的方法相同。首先挑选20个与WRF-3DVar系统中的背景场误差协方差一致的流函数扰动场,即控制变量场;然后采用递归滤波、经验正交函数进行水平变换和垂直变换,再通过物理转换实现从控制变量所在空间向模式空间的转换。由此导出地转平衡的扰动变量场,包括水平风场、温度场、水汽混合比、气压扰动场等,而垂直速度及云水、雨水、冰晶、雪、霰比含水量等其他预报变量不做扰动。 5 雷达资料质量控制及预处理 5.1 雷达资料质量控制

对雷达资料做了如下预处理和质量控制:

(1)使用三维反射率结构特征对地物回波进行抑制:噪声点回波过滤,通过计算反射率水平纹理、反射率垂直梯度及限定反射率垂直梯度的高度滤除地物回波,最后对降水回波洞进行填补,具体做法见肖艳姣(2007)文章。

(2)雷达资料稀疏化:径向风和反射率因子观测都进行了稀疏化处理,均为径向每4 km取一个观测,稀疏化后的观测空间分辨率与模式格距相当。

(3)背景场检查:雷达反射率因子观测误差假设为5 dBz,雷达径向风观测误差假设为3 m/s,当观测与背景场的偏差大于3倍观测误差时,剔除该观测。

5.2 进入同化分析的雷达资料

经过资料质量控制和预处理,进入同化分析的雷达资料被剔除了很大一部分(表 1),反射率只占一次完整体扫的2%-3%,径向风约只占一次完整体扫的1%。于09时、09时30分、10时共进行3次同化分析,分别对应采用09时01分、31分、10时01分的雷达体扫资料。图 4为经过质量控制后进入第一次同化分析的观测数值及分布,图中将不同高度的所有观测均投影到同一水平面上,用不同颜色代表观测的数值大小,第2、第3次同化的观测分布及数值图略。

表 1给出了3次同化分析的总观测数目、经过质量控制后进入分析的观测数目及其占总观测的百分比。

表 1 同化时刻的雷达观测数目 Table 1 The number of radar observations assimilated at the scanning times
同化循环雷达体扫时刻总观测数目进入分析的观测数目占总观测的百分比
反射率径向风反射率径向风反射率径向风
第1次09时01分764898232717322892199652.99﹪0.86﹪
第2次09时31分764606232641321194199102.77﹪0.86﹪
第3次10时01分765659232837416971190162.22﹪0.82﹪

图 4 09时进入第1次同化分析的观测数值及分布(a.反射率因子,b.径向风) Fig. 4 Observations of reflectivity(a) and radial velocity(b)be assimilated in the first assimilation cycle at 09:00 BT
6 同化试验结果 6.1 集合方根滤波对飑线系统的分析能力

将控制试验的预报场、同化试验的集合平均预报场、集合平均分析场与观测进行对比,通过对3次同化分析前后的风场、水凝物场、降水量等进行比较,评估同化分析对系统流场结构和微物理量场的分析反演能力。

6.1.1 对风场结构的分析

先通过对比0.5°仰角的径向风PPI观测与控制试验、同化试验的风场,检验集合方根滤波同化系统对低层水平风场的分析能力。2005年7月12日09-10时,济南雷达的西北方有离开雷达的径向风大值区,径向风速达19 m/s,东南方向及南部有朝向雷达5-10 m/s的径向风,可见地面至850 hPa左右的低层,雷达西北方向的清河、夏津、临清、武城等地主要为较强的东南风及偏南风,雷达西南部的阳谷县附近为弱的西南风,南部的东平县附近为弱的偏南风,在山东西部、河北东部两省交界处存在南西南风转东南风的气旋式变化。

3次同化分析后集合平均分析场的雷达径向风可以通过径向风观测算子(式(9))计算得到,图 6为相应时间模式第5层(约900 hPa)雷达径向风分布。对比图 5可见,同化分析后径向风的大小和方向与实况基本相符。

图 5 2005年7月12日0.5°仰角雷达径向风观测PPI(a.09时01分,b.09时31分,c.10时01分) Fig. 5 Radial velocity PPI observations at 0.5 elevation at 09:01(a),09:31(b),and 10:01 BT(c)12 July 2005

图 6 2005年7月12日3次同化分析后的低层雷达径向风分布(a.09时,b.09时30分,c.10时;黑点为济南雷达天线位置) Fig. 6 Radial velocity around 900 hPa calculated from the ensemble mean analysis field at 09:00(a),09:30(b),and 10:00 BT(c)12 July 2005(the thick black dot denotes the location of the Jinan radar station)

图 7为09时30分控制试验预报场、同化试验集合平均预报场和分析场的900 hPa水平风矢量分布,09时、10时的图略。可见在控制试验模拟的900 hPa水平风场,山东西部及中部都为均匀南向风场,风速较小,与实况观测到的鲁西北低层较强东南气流及气旋式切变等风场分布特征不一致;而同化分析后,集合平均分析场低层的径向风从大小到方向都与实况对应较好,清河、夏津、临清等地出现了大于19 m/s的背离雷达的径向风(图 6);同化分析后的900 hPa预报场和分析场,在山东西部与河北交界处,有一明显的南西南风转东南风的风向转变,这与雷达低层径向风观测一致;同时,同化分析后,雷达北部陵县附近南风比控制试验也有所增强,与实况更为接近。低层水平风切变,清河、夏津、临清等地出现东南风增强及陵县附近南风增强等现象,引起了鲁西北地区低层气流辐合,有助于模式大气产生或增强上升气流,为后期强对流天气过程的发生提供了准备。

图 7 2005年7月12日09时30分900 hPa水平风矢量分布(a.控制试验预报场,b.第2次同化分析前的集合平均预报场,c.同化后的集合平均分析场) Fig. 7 Horizontal wind vector at 900 hPa from the forecast of control run(a),ensemble mean forecast(b) and ensemble mean analysis(c)of the second assimilation cycle at 9:30 BT 12 July 2005

将控制试验、同化试验高层的径向风与1.5°、2.4°仰角的径向风观测进行对比发现,同化后分析场的高层风比控制试验更接近观测实况。从09时01、31分、10时01分的1.5°、2.4°仰角的雷达径向风PPI观测可见(图略),山东西部及西北部、河北东南部的500 hPa高空附近有较强的西北气流;山东中部则为弱的西南偏西气流,雷达附近的高空对应西北风转西南风的槽线,雷达南部的鲁西南地区的500 hPa高空附近为弱的西南偏西气流。对比从集合平均分析场计算出的模式第10层(约550 hPa)雷达径向风可见,同化分析后,高层径向风的方向和大小都与雷达观测吻合较好,西北象限为朝向雷达的径向风,入流速度达19 m/s,西南象限径向风朝向雷达,东面径向风离开雷达。

对比09时、09时30分、10时控制试验预报场、同化试验集合平均预报和分析场的500 hPa水平风矢量分布可见,控制试验虽然模拟出了500 hPa槽线的主要分布形势,但流场相对均匀平滑,没有反应出雷达西北方向德州-衡水-辛集-无极一带出现的强西北气流,而经过同化后,这些区域的西北风得到加强,且在冀东南河间、献县的风场也得到改进。

下面考察集合方根滤波系统对垂直风场的分析能力。从上述水平风分析可知,09-10时在37.2°N附近的底层有风场的辐合,此处有对流云团发展旺盛,对应有较强的垂直运动。对比控制试验预报场、同化试验预报场和分析场在09时、09时30分、10时沿37.2°N剖面的垂直速度、总水凝物(即云水、雨水、冰晶、雪、霰比含水量的总和)及温度可见,控制试验中垂直速度增加缓慢,直到10时才达到1 m/s,控制试验的总水凝物含量也较小(图 8b),如此弱的上升运动和低含量的水凝物与实况的强对流不符;同化试验第1次同化前的集合平均预报场垂直上升运动很弱,不到0.5 m/s(图 8c),但经过第1次同化后,上升运动明显增强,中心垂直速度达4 m/s(图 8d)。模式预报水凝物也因为上升运动迅速发展起来,可以看到,总水凝物含量与上升运动有较好的对应,10时的垂直运动较09时30分有所减弱,这与实况观测到的这一时期对流云团处于逐渐衰弱阶段相吻合。

图 8 2005年7月12日沿37.2°N的垂直速度(色阶)、温度(绿线,℃)及总水凝物(黑实线,g/kg)剖面(a.09时控制试验预报场,b.10时的控制试验预报场,c.09时同化试验集合平均预报场,d.09时同化试验集合平均分析场) Fig. 8 Cross-section diagram of vertical velocity(shaded),temperature(green line) and total hydrometeor(black solid line,g/kg)along 37.2°N from the forecast of control run(a)at 09:00,forecast of control run(b)at 10:00,ensemble mean forecast(c)at 09:00 and ensemble mean analysis(d)of the second assimilation cycle at 09:00 BT 12 July 2005

可见,通过集合方根滤波雷达资料同化,模式能更准确地分析出低层流场的气旋式切变、强风速区及夏津等地的低层风场辐合等细节特征,使垂直运动得到增强,促进了对流的发生、发展,雷达观测资料同化分析使流场与实况更为吻合。

6.1.2 微物理量场的反演

从0.5°仰角的雷达回波PPI(图 9)和1.5°仰角的雷达回波PPI(图略)可见,雷达西北部出现发展旺盛的对流云团,中心最大回波均达到60 dBz,09-10时此对流云团逐渐衰弱,由回波位置及仰角度数推算,回波高值中心所在高度即2-5.5 km处(700-500 hPa)应有水凝物的大值区。但控制试验结果表明,模式由于起转的原因,对流尚未发 展起来,500和700 hPa上水凝物仅在局地出现,且数值较小,只有局部地区达到1.0 g/kg(图 10a1、b1、c1;控制试验700 hPa总水凝物分布图略)。经过同化分析后,水凝物信息被加入到了分析场,大值中心得到明显加强,大部分区域达到了1.0 g/kg,水凝物分布得到显著改进,在500和700 hPa的山东西北地区都出现了与实况更吻合的水凝物分布,水凝物大值区也与雷达回波高值区对应较好(图 10a2、b2、c2;同化试验700 hPa总水凝物分布图略)。

图 9 2005年7月12日0.5°仰角雷达反射率PPI(a.09时01分,b.09时31分,c.10时01分) Fig. 9 Reflectivity PPI observations at 0.5 elevation at 09:01(a),09:31(b),and 10:01 BT(c)12 July 2005

图 10 2005年7月12日09时(a)、 09时30分(b)、10时(c)500 hPa总水凝物的分布(a1、b1、c1.控制试验预报场,a2、b2、c2.同化后集合平均分析场) Fig. 10 Forecasted total hydrometeor(g/kg)by the control run at 09:00(a1),09:30(b1),10:00 BT(c1),and the ensemble mean analysis at 09:00(a2),09:30(b2),10:00 BT(c2)12 July 2005 at 500 hPa
6.2 同化分析场对飑线系统的预报能力

基于同化分析对飑线系统动力场、热力场等的准确反演,将经过3个同化循环后的分析场预报至20时,与控制试验对应时刻的预报场及雷达观测进行比较,检验集合方根滤波同化分析场对飑线系统的预报能力。

6.2.1 对微物理量场的预报能力

首先通过比较模式预报的雷达反射率因子与实况雷达反射率因子,来检验集合平均分析场对微物理量场的预报能力。

雷达1.5°仰角时其径向150 km处高度约为4 km,即约600 hPa。由图 11a可见,济南南部上空有一块弓状回波,强度约为50 dBz。控制试验预报对应时刻的反射率因子分布与观测差别较大,回波中心强度基本为55-65 dBz,在数值上也整体偏大(图 11b);而由集合平均分析场3 h预报的反射率因子分布与雷达反射率观测更为吻合,强回波中心和弱回波区域的分布也基本一致,在数值上与控制试验相比也更接近观测(图 11c)。将同化试验和控制试验预报的14时雷达反射率因子与观测进行对比可见,集合平均分析场在经过4 h预报后仍然有较好的预报效果(图略)。

图 11 2005年7月12日1.5°仰角13时03分雷达反射率因子PPI(a)、控制试验(b)和同化试验分析场(c)预报13时600 hPa雷达反射率因子 Fig. 11 Reflectivity PPI observations at 1.5 elevation at 13:03(a),forecasted reflectivity at 600 hPa by control run(b) and ensemble mean analysis of assimilation experiment(c)at 13:00 BT 12 July 2005

5 h后,集合平均分析场对水凝物的预报准确率逐渐下降,主要表现为预报的飑线系统发展得比实况更快,系统移动速度较实况更快(图略)。这可能与同化进来的西风偏大有关。但与控制试验相比,同化后的分析场的预报效果还是与实况更为接近。经过8 h的预报后,不论是控制试验还是同化试验,其预报结果都与实况相差较大。

下面从反射率垂直剖面来考察集合平均分析场对系统垂直结构的预报能力。图 12a为14时09分反射率因子1.5°仰角PPI,图 12b为沿图a所示线段AB所做的剖面,剖面位置约在35.7°N,图 13为模拟的反射率沿35.7°N所做的剖面。由反射率因子剖面图可见,飑线系统呈明显带状排列,与控制试验相比同化后的分析场更准确地预报出了飑线系统多个对流单体的带状分布,各单体发展的高度及中心强度都与实况更加吻合。

图 12 2005年7月12日14时09分的反射率因子观测(a.1.5°仰角PPI,b.沿图a所示线段AB的剖面) Fig. 12 Reflectivity observation at 1.5° elevation as shown by the PPI(a)with its cross-section(b)along line AB in(a)at 14:09 BT 12 July 2005

图 13 控制试验(a)和同化试验分析场(b)预报的2005年7月12日14时反射率因子沿35.7°N的剖面 Fig. 13 Forecasted cross-section reflectivity along 35.7°N by control run(a) and ensemble mean analysis of assimilation experiment(b)at 14:00 BT 12 July 2005

由MODIS卫星7月12日10时50分观测的云可见光信息可见,在山东西部地区有发展深厚的云系,系统范围小,外围云系较中心浅薄(图 14)。控制试验直接预报的500 hPa总水凝物比MODIS观测范围小,同化雷达资料后的集合平均预报场所做的预报改进了水凝物的分布,与观测更为一致,云系的带状结构分布也更加合理(图 15)。

图 14 MODIS卫星7月12日10时50分观测的可见光云图 Fig. 14 Visible cloud picture by the MODIS at 10:50 BT 12 July

图 15 控制试验(a)和同化试验分析场(b)预报的2005年7月12日11时500 hPa总水凝物分布 Fig. 15 Total hydrometeor at 500 hPa forecasted by control run(a) and ensemble mean analysis of assimilation experiment(b)at 11:00 BT 12 July 2005
6.2.2 对流场及风暴结构的预报能力

图 16为用3次同化循环后的集合平均分析场预报4 h的500和950 hPa的水平风场和总水凝物,可见从低层到高层环境风有较强的垂直切变,800 hPa以上环境风主要为西北偏北风转偏西风,800 hPa以下环境风切变为西南偏南风转偏东风,高、低层风切变达90°-180°,如图 16中红圈标注所示位置,这种风的垂直切变符合风暴的流场结构特点,有利于风暴的维持和发展(俞小鼎等,2006)。

图 16 同化试验分析场预报的2005年7月12日14时500(a)、950 hPa(b)水平风场(矢量)和总水凝物(色阶) Fig. 16 Forecasted horizontal wind field(vector) and total hydrometeor(shaded)by ensemble mean analysis at 500 hPa(a) and 950 hPa(b)at 14:00 BT 12 July 2005

图 16b可见,在低层风暴的前沿为下沉气流的出流区,水平风速很大,达20 m/s,这是降水粒子在湿绝热下沉过程中蒸发冷却和融化,其拖曳作用形成了冷性下沉气流,在近地面层向外扩散,冷性出流与南面的暖湿西南气流在风暴的前沿强烈辐合交汇,形成阵风锋,这与此处14时地面的大风降温天气实况一致,阵风锋附近强烈的辐合上升运动和不稳定能量驱动了风暴的继续发展。由14时预报的总水凝物场和温度场沿35.4°N的剖面可见,对流单体在边界层附近存在冷池(图略),这是由降水粒子的冷性下沉气流所致,同时单体的中部对应暖的温度扰动,可能与单体中部水凝物的凝结潜热释放有关,符合对流单体的动力和热力结构特征。

通过上述分析可见,由集合平均分析场预报的系统流场及结构符合风暴的动力特征,动力场和热力场的分布、配置基本合理。

6.2.3 对地面降水的预报能力

由降水实况(图 17)可见,08-10时河北省东南部、河北省与山东省交界处、山东省北部出现了降水,局部地区小时降水量超过30 mm。控制试验从08时开始起报,由于起转原因,模式积分1.5 h仍无降水产生(图 18a);而同化试验中,在经过09时第1次同化分析后,模式经过半小时的继续积分,很快响应出了降水(图 18b),这也与同化对水凝物场的准确分析有关,对比09-10时实况可见,同化试验09时00-30分对山东省西部与河北省交界处的降水预报比较成功,虽然量值偏大,但降水的位置与实况对应很好。可见,集合方根滤波系统对雷达资料的同化有效缩短了模式起转时间,明显改进了降水预报。

图 17 2005年7月12日08-09时(a)、09-10时(b)地面逐时降水(单位:mm)实况 Fig. 17 Hourly surface precipitation observations(mm)during 08:00-09:00 BT(a),09:00-10:00 BT(b)12 July 2005

图 18 控制试验(a)、同化试验(b)预报的2005年7月12日08时至09时30分累积降水 Fig. 18 Forecasted accumulated precipitation by control run(a) and assimilation experiment(b)during 08:00 to 09:30 BT 12 July 2005

由08-20时的累积降水实况观测可见,降水主要位于山东省中西部,呈西北-东南向分布,强降水中心位于(36.3°N,116.8°E),强度超过60 mm(图 19a)。控制试验预报的降水都普遍高于实况,特别是在山东省北部出现了虚假降水(图 19c)。采用集合方根滤波同化后的雷达资料分析场做确定性预报,则明显改进了降水的预报能力:减弱了山东省北部和半岛地区的虚假降水,较好地预报出飑线附近降水的位置和强度;同时,同化后的集合平均分析场还较准确预报出半岛、日照等地的局地降水(图 19b)。

图 19 08-20时累积降水(a.实况;b.同化试验分析场12 h降水预报;c.控制试验12 h降水预报) Fig. 19 Accumulated precipitation from the observation(a),the forecast by control run(b) and the forecast by assimilation experiment(c)during 08:00 to 20:00 BT 12 July 2005
7 结论与讨论

本研究基于WRF模式采用集合方根滤波方法同化多普勒雷达资料,对山东省中西部地区一次飑线天气过程进行了同化应用试验,并结合雷达、卫星、实况降水等观测资料对同化后的分析和预报效果进行了对比检验,主要结论如下:

(1)基于WRF模式的集合方根滤波同化系统能有效同化实际雷达资料,30 min间隔的雷达资料同化有效增加了模式的中小尺度信息,较准确地分析出了飑线系统流场结构的细致特征,如低层流场的气旋式切变、强风速区及低层风场辐合等特征,对微物理量场也能进行较准确的反演,并由此有效缩短了模式起转时间,明显改进了地面降水的预报。

(2)基于同化分析对飑线系统动力场、热力场等的准确反演,利用3次同化循环后的集合平均分析场开展了预报试验。结果表明,与控制试验相比,同化后分析场能更准确地预报出对流系统的微物理量场,尤其是在前3 h效果非常明显,预报的流场结构符合风暴的动力特征,动力场和热力场的分布、配置也基本合理。

(3)集合平均分析场对飑线系统传播方向的预报与实况一致,但预报的系统传播速度较实况快,由于对流系统的非线性发展迅速,对系统的预报时效只有5-6 h。

发展旺盛的中小尺度对流系统(如飑线、龙卷等),云体发展深厚,云顶高度多在0℃等温线以上,云内不仅存在着液态水凝物(云水、雨滴),也存在冰晶、雪、霰或雹等冰相粒子水凝物。这些冰相水凝物及其微物理过程在对流系统的发生和发展中具有重要的作用。雷达反射率因子反映了大气中云雨粒子的尺度和数密度,既包括了液态水凝物也包括了冰相水凝物,若数值模式中的云方案仅包含液态水凝物,那么在同化更新中将不能分析反演出天气系统中水凝物的真实状况,并给同化后的云场带来虚假信息,也将直接影响其他要素的预报准确性。由6.1.2和6.2.1节结果可见(图 10图 15),同化分析和预报的500 hPa水凝物都与实况有较一致的分布,而500 hPa总水凝物贡献最大的主要来自于雪、霰等冰相水凝物,能获得大小与分布都与实况较一致的水凝物的分析和预报,与同化雷达反射率观测时同化方案对冰相水凝物的准确分析,及预报模式中使用了包含复杂冰相微物理过程的云方案有密切关系。

与假设模式完美的观测系统模拟试验相比,实际雷达资料同化试验存在很多困难,最大的困难在于模式误差的引入,这种模式误差既包括了数值模式对实际大气描述的误差(格点离散化、物理过程参数化、模式边界条件等),也包括了观测算子将观测解释为模式变量中的误差。对于对流尺度系统,微物理过程参数化可能是模式误差和不确定性的主要来源,它能较显著地影响风暴的结构和演变,进一步的工作可以尝试在集合成员中使用多种物理过程方案,以评估由微物理过程带来的模式误差。格点离散化是模式误差的另一来源,可以尝试通过提高分辨率部分减轻由分辨率带来的模式误差。针对模式误差引起的滤波发散问题,在实际应用中,对协方差扩大因子进行调整是解决此问题简单而有效的手段,在本研究的同化试验中就是采用这个方法来解决滤波发散问题,进一步的工作可以尝试对同化系统进行改进,在每次集合方根滤波分析中动态调整协方差扩大因子,从而使集合发散度在每次分析时都能自动调整到对集合平均误差的最好估计。

集合方根滤波初始扰动场的生成有多种方法,在实际雷达资料同化时,分析和预报的质量受到各方面不确定性因素的影响,恰当的初始扰动场能减少来自环境场的不确定性。本研究采用WRF模式的三维变分系统来产生初始扰动场,由此得到的初始扰动场是地转平衡的,而积云尺度对流系统变量间的关系是高度非地转的,在初始时刻生成适合于对流系统的扰动场也许能一定程度上改进同化效果。

在对同化效果进行评估时,所能获取到用于检验的观测主要包括探空站资料、常规地面观测、地面自动站观测和卫星观测,本研究还使用了雷达资料进行检验,由于这些资料的时空分辨率不足以全面描述系统的对流尺度特征,雷达观测虽然有足够的分辨率,但所探测的物理量都不是模式变量,所以,在对同化效果的评估上存在一定困难。此外,由于观测资料本身也包含观测误差,所以对集合方根滤波同化实际雷达资料效果的评估具有一定的不确定性。

参考文献
胡志晋, 楼小凤, 包绍武等. 1998. 一个简化的混合相云降水显式方案. 应用气象学报, 9(3): 257-264
秦琰琰, 龚建东, 李泽椿. 2012. 集合卡尔曼滤波同化多普勒雷达资料的观测系统模拟试验. 气象, 38(5): 513-525
肖艳姣. 2007. 新一代天气雷达三维组网技术及其应用研究. 南京: 南京信息工程大学, 128pp
许小永, 郑国光, 刘黎平. 2004. 多普勒雷达资料4DVAR同化反演的模拟研究. 气象学报, 62(4): 410-422
许小永, 刘黎平, 郑国光. 2006. 集合卡尔曼滤波同化多普勒雷达资料的数值试验. 大气科学, 30(4): 712-728
杨艳蓉, 王振会, 杨洪平等. 2008. 多普勒雷达反射率与径向风资料在数值模式中的应用试验. 气象, 34(6): 26-34
杨毅, 邱崇践, 龚建东等. 2008. 利用3维变分方法同化多普勒天气雷达资料的试验研究. 气象科学, 28(2): 124-132
俞小鼎, 姚秀萍, 熊廷南等. 2006. 多普勒天气雷达原理与业务应用. 北京: 气象出版社, 314pp
Dowell D C, Zhang F Q, Wicker L J, et al. 2004. Wind and temperature retrievals in the 17 May 1981 Arcadia, Oklahoma supercell: Ensemble Kalman filter experiments. Mon Wea Rev, 132(8): 1982-2005
Evensen G. 1994. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J Geophys Res, 99(C5): 10143-10162
Gal-Chen T. 1978. A method for the initialization of the anelastic equations: Implications for matching models with observations. Mon Wea Rev, 106(5): 587-602
Gaspari G, Cohn S E. 1999. Construction of correlation functions in two and three dimensions. Quart J Roy Meteor Soc, 125(554): 723-757
Hane C E, Ray P S. 1985. Pressure and buoyancy fields derived from Doppler radar data in a tornadic thunderstorm. J Atmos Sci, 42(1): 18-35
Hong S Y, Dudhia J, Chen S H. 2004. A revised approach to ice microphysical processes for the bulk parameterization of clouds and precipitation. Mon Wea Rev, 132(1): 103-120
Houtekamer P L, Mitchell H L. 1998. Data assimilation using an ensemble Kalman filter technique. Mon Wea Rev, 126(3): 796-811
Houtekamer P L, Mitchell H L. 2001. A sequential ensemble Kalman filter for atmospheric data assimilation. Mon Wea Rev, 129(1): 123-137
Jones C D, Macpherson B. 1997. A latent heat nudging scheme for the assimilation of precipitation data into an operational mesoscale model. Meteor Appl, 4(3): 269-277
Kain J S, Fritsch J M. 1990. A one-dimensional entraining/detraining plume model and its application in convective parameterization. J Atmos Sci, 47(23): 2784-2802
Kain J S, Fritsch J M. 1993. Convective parameterization for mesoscale models: The Kain-Fritcsh scheme. The representation of cumulus convection in numerical models. Amer Meteor Soc, 84pp
Lin Y, Ray P S, Johnson K W. 1993. Initialization of a modeled convective storm using Doppler radar-derived fields. Mon Wea Rev, 121(18): 2757-2775
Lin Y L, Farley R D, Orville H D. 1983. Bulk parameterization of the snow field in a cloud model. J Climate Appl Meteor, 22(6): 1065-1092
Smith P L Jr, Myers C G, Orville H D. 1975. Radar reflectivity factor calculations in numerical cloud models using bulk parameterization of precipitation. J Appl Meteor, 14(6): 1156-1165
Snyder C, Zhang F. 2003. Assimilation of simulated Doppler radar observations with an ensemble Kalman filter. Mon Wea Rev, 131(8): 1663-1677
Sun J Z, Crook N A. 1997. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint: Part I. Model development and simulated data experiments. J Atmos Sci, 54(12): 1642-1661
Tong M J, Xue M. 2005. Ensemble Kalman filter assimilation of Doppler radar data with a compressible nonhydrostatic model: OSS experiments. Mon Wea Rev, 133(7): 1789-1807
Whitaker J S, Hamill T M. 2002. Ensemble data assimilation without perturbed observations. Mon Wea Rev, 130(7): 1913-1924
Xue M, Wang D, Hou D, et al. 1998. Prediction of the 7 May 1995 squall line over the central U. S. with intermittent data assimilation//Preprints 12th Conference on Numerical Weather Prediction, Phoenix, AZ: Amer. Meteor. Soc., 191-194
Xue M, Tong M J, Droegemeier K K. 2006. An OSSE framework based on the ensemble square root Kalman filter for evaluating the impact of data from radar networks on thunderstorm analysis and forecasting. J Atmos Ocean Technol, 23(1): 46-66
Zhang J, Carr F, Brewster K. 1998. ADAS cloud analysis//Preprints 12th Conference on Numerical Weather Prediction. Phoenix, AZ: Amer Meteor Soc, 185-188
Zhang F, Snyder C, Sun J. 2004. Tests of an ensemble Kalman filter for convective-scale data assimilation: Impact of initial estimate and observations. Mon Wea Rev, 132(5): 1238-1253
Zhang F Q, Meng Z Y, Aksoy A. 2006. Tests of an ensemble Kalman filter for mesoscale and regional-scale data assimilation. Part Ⅰ: Perfect-model experiments. Mon Wea Rev, 134(2): 722-736