中国气象学会主办。
文章信息
- 乔小湜, 闵锦忠, 王世璋 . 2016.
- QIAO Xiaoshi, MIN Jinzhong, WANG Shizhang . 2016.
- 集合卡尔曼滤波同化中雷达位置的敏感性研究
- A sensitivity study of radar location in data assimilation using the ensemble Kalman filter
- 气象学报, 74(5): 796-814.
- Acta Meteorologica Sinica, 74(5): 796-814.
- http://dx.doi.org/10.11676/qxxb2016.050
-
文章历史
- 2016-03-07 收稿
- 2016-05-17 改回
2. 沈阳中心气象台, 沈阳, 110016
2. Shenyang Central Meteorological Observatory, Shenyang 110016, China
强对流天气是造成人民生命财产损失的重要气象灾害之一,据中国气象局官网(http://www.cma.gov.cn/)报道,2001—2004年天气灾害损失中约一半是由强对流造成的。强对流天气因其局地性和突发性强成为当前天气预报的一个难点(郑永光等,2010)。数值预报是提高强对流天气预报准确率的重要技术之一,而数值预报问题是一个初边值问题,其中初值准确性的提高主要通过同化观测资料来实现(Kalnay, 2003)。当前,强对流研究使用的时空分辨率最高的观测资料来自多普勒天气雷达,因此,如何有效地同化雷达观测资料是数值预报研究的前沿热点问题。
目前雷达资料同化的主要方法有三维变分(3DVAR)(Hu et al, 2007a, 2007b)和集合卡尔曼滤波(EnKF)(Evensen, 1994)。3DVAR使用的是静态背景误差协方差,不随时间变化,考虑到强对流天气系统快速发生、发展的特征,静态协方差在循环同化中并不能合理地代表误差的真实结构。同时,变分方法需要使用切线性算子和伴随算子,对其分析的精度也会造成一定影响(Sun et al, 1997; Whitaker et al, 2004; Kawabata et al, 2011)。集合卡尔曼滤波通过流依赖的背景误差协方差来估计误差结构随时间的变化,可以使用非线性观测算子,尤其是非线性的反射率观测算子,因此,集合卡尔曼滤波成为对流尺度资料同化研究的一个热点,在模拟观测和实际观测的研究中都得到了大量的应用(Tong et al, 2005; Coniglio et al, 2006; Xue et al, 2006; Anderson et al, 2009; Dowell et al, 2009; Raeder et al, 2012; Sobash et al, 2013; Chang et al, 2014; Snook et al, 2015)。
尽管一系列的研究已经显著改进了集合卡尔曼滤波同化雷达资料的效果,但由于强对流天气具有较短生命史及较强突发性和局地性的特点,所以,集合卡尔曼滤波同化系统要在发现强对流天气后能够在较短的时间内完成同化并进行预报。针对该问题,闵锦忠等(2012)、王世璋等(2013)、Wang等(2013a)根据Kalnay等(2010)的RIP(running in place)方案,设计了针对对流尺度同化的迭代集合均方根滤波(iEnSRF)方法,有效地降低了集合卡尔曼滤波短时雷达资料同化的分析误差,该方案对同化后的短时确定性预报有较为显著的正贡献。该迭代方案将一个同化循环内的单次同化调整为多次同化,在每一步迭代中,使用分析时刻获得的观测更新前一时刻的集合,然后再进行集合预报到分析时刻,重复这一“预报—分析”过程直到最后一步迭代,更新分析时刻的集合完成整个迭代步骤,并进入下一同化循环。
当前业务化天气雷达网使用的都是固定地基雷达,雷达位置不能移动,因而强对流天气系统有可能出现在雷达的任意一个方位。强对流天气较强的局地性要求集合卡尔曼滤波能够有效地分析出现在雷达各个方位的对流系统,而不是如车载雷达一样先移动到有利位置再进行观测和同化。已有研究大都假设同化系统(包括3DVAR、4DVAR和EnKF)具备这种能力,因此,目前关于雷达位置对同化影响的讨论相对较少。然而,部分文献的研究表明以上假设并不总是成立。Sun等(1994)在研究4DVAR反演阵风锋上的雷达径向风观测资料时发现,当雷达所处位置垂直于阵风锋的锋面时,雷达径向风并不能很好地反映阵风锋系统的流场结构特征,利用该雷达观测的径向风反演出来的风场与阵风锋的结构存在较大差距。Sun等(1997)通过云模式进行4DVAR反演雷达观测时发现,当使用双雷达进行反演时,双雷达之间的相对位置对4DVAR的风场反演结果并没有显著的影响。但是许小永等(2004)关于雷达位置对4DVAR双雷达反演影响的讨论表明,当水平风速的某一个分量在两部雷达的径向上的投影较小时,4DVAR对该分量的反演误差就较大。随后,许小永等(2006)在集合卡尔曼滤波同化雷达资料的研究中发现,雷达位置对集合卡尔曼滤波同化的前几个循环有较为明显的影响,但是在进行了更多的同化循环后,影响基本消失。
上述研究表明,雷达资料同化初期对雷达位置有一定的敏感性,但是讨论多集中于4DVAR反演雷达资料方面,与集合卡尔曼滤波雷达资料直接同化相关的讨论较少。许小永等(2006)仅分析了位于同化区域西侧和南侧的雷达对集合卡尔曼滤波同化效果的影响,并未对影响原因以及其他方位的雷达影响进行讨论。同时,由于当前大多数集合卡尔曼滤波雷达资料同化研究采用了长时间同化方案,如连续同化若干小时(Tanamachi et al, 2013; Putnam et al, 2014; Chang et al, 2014; James et al, 2014),因此,对短时集合卡尔曼滤波同化中可能存在的问题研究较少。但是强对流系统快速发展变化的特点对短时集合卡尔曼滤波同化提出了需求。那么雷达位置是否会给短时同化和随后预报带来影响?对短时集合卡尔曼滤波同化改进的迭代集合均方根滤波方法能否降低集合卡尔曼滤波对雷达位置的敏感性?这些都需要进行研究。已有研究表明雷达位置对水成物同化的影响较小,本研究将主要讨论对风场的影响。为了易于根据需要调整观测雷达的位置,并知道真值和真实误差的分布,本研究将基于模拟观测试验(OSSE,Observation Simulation System Experiment)进行。同时,考虑到在当前业务雷达布网密度条件下对流系统有可能被多个雷达观测到,将分别讨论单雷达同化和双雷达同化。
2 集合卡尔曼滤波同化方法集合卡尔曼滤波最早由Evensen(1994)引入海洋同化,经过20年的发展已经衍生出多个分支,包括集合调整卡尔曼滤波(EAKF)(Anderson, 2001)、集合转换卡尔曼滤波(ETKF)(Bishop et al, 2001)和集合均方根滤波(EnSRF)(Whitaker et al, 2002)。其中集合均方根滤波在对流尺度雷达资料同化中得到了广泛应用(Xue et al, 2006),文中将基于该算法进行研究。
2.1 集合均方根滤波算法传统的集合卡尔曼滤波需要对观测进行扰动才能使集合卡尔曼滤波的分析误差协方差与卡尔曼滤波的理论值相符。但是Whitaker等(2002)指出,对观测进行扰动会引入额外的采样误差,使得集合卡尔曼滤波实际的分析误差协方差被低估,为此提出了一个确定性的算法——集合均方根滤波。集合均方根滤波采用顺序同化,即对当前同化循环内的观测一个接一个进行分析。对于第j个观测,集合均方根滤波的分析方程可以表述为由集合平均(x)的更新方程和集合扰动(x′)的更新方程(以第i个集合成员为例)两部分组成
(1) |
(2) |
其中,
(3) |
为卡尔曼增益。式中,上标a和b分别表示分析场和背景场;Pb是模式空间的背景误差协方差矩阵;H是观测算子;H是H的切线性算子;R在只同化一个观测的条件下表示该观测的误差方差;α是集合均方根滤波引入的小参数,其在单一观测条件下的解(Maybeck,1979)为
(4) |
α的引入使得集合均方根滤波在不扰动观测的情况下不会低估分析误差协方差。
2.2 迭代集合均方根滤波算法迭代集合均方根滤波算法是通过在两个时刻(分析时刻tn和一个较早的没有观测的时刻tn-1*)之间反复进行“预报—分析”(更新较早时刻的背景场;即迭代运算)来改进背景场以及背景误差协方差较差情况下集合卡尔曼滤波同化分析的效果。迭代集合均方根滤波与传统集合均方根滤波方法的主要差别是引入了Sakov等(2010)的AEnKF(Asynchronous EnKF)算法,实现了分析时刻观测对另一个时刻模式背景场的更新。迭代集合均方根滤波在使用tn时刻的观测对tn-1*的背景场进行更新时,使用的分析方程为
(5) |
(6) |
其中, 异步卡尔曼增益矩阵的计算方程为
(7) |
(8) |
(9) |
式中,N为集合成员的个数,其余与前面说明相同。
3 同化模式和试验设计观测模拟试验选用基于WRF(Weather Research and Forecasting Model)研发的集合均方根滤波同化系统,该系统已在Wang等(2013b)、王世璋等(2013)、陈杰等(2012)及闵锦忠等(2013)的研究中得到了检验和应用,具有一定的可靠性。
3.1 模式设置和真实风暴模拟进行观测模拟试验的第一步是模拟一个“真实风暴”作为生成模拟观测和同化试验对比的真值。预报模式采用WRF ARW V2.2.1,模拟区域水平网格为30×30,垂直40层。网格分辨率在水平方向为2 km,在垂直方向为0.5 km。理想风暴的环境场通过一个探空资料生成,该探空数据为WRF软件包自带,垂直结构与Weisman等(1982)所用探空相似。该环境场中各个变量在水平方向上均匀分布,例如在同一高度上各网格点上风等要素是相同的。理想风暴是由在上述环境场中加入一个中心强度为3 K的椭圆形热泡(Yussouf et al,2010)触发而成,热泡在水平方向上的最大半径为6 km,在垂直方向上的最大半径为1.5 km。热泡中心位于x=14 km、y=28 km和z=1.5 km处。
模式的积分方案为Runge-Kutta 3阶时间积分方案,积分步长为12 s,“真实风暴”模拟时间为1.5 h。模式的各参数化方案包括WSM6微物理方案、Dudhia短波辐射方案、RTTM长波辐射方案和1.5阶湍流动能(TKE)次网格扰动方案。模拟区域的侧边界为开放边界,上、下边界为刚性,下垫面平坦无地形。设置的详细介绍见Wang等(2013a)。该超级单体风暴位置随时间的变化如图 1所示。从图中可以看出该风暴处于西风环境气流中,沿着气流自西向东移动,并在移动过程中不断增强并出现分裂,其中向风暴整体前进方向右侧移动的单体在后期发展出类似于钩状回波的结构,表明这是一个较强的超级单体风暴。同时该风暴的移动速度较低,在整个模拟过程中都处于模拟区域的中间位置。
3.2 模拟雷达观测WRF-EnSRF同化系统能够同化雷达径向风(Vr)和反射率(Z)观测,使用与Tong等(2005)相同的观测算子。雷达径向风的观测算子为
(10) |
式中,ug、vg和wg分别为模式预报风速的3个分量在雷达观测点上的插值;wt为降水粒子的下落末速度;α和β分别为雷达观测的仰角和方位角。雷达反射率的观测算子包括雨水(qr)、雪(qs)和雹霰(qg)3部分
(11) |
式中,Zr、Zs和Zh分别是雨水、雪和雹霰的相当反射率因子,具体见Tong等(2005)的式(4)—(7)。
模拟雷达观测的生成方式采用美国业务WSR-88D雷达降水扫描模式(VCP11),共有14层仰角,最低0.5°,最高19.5°(Xu et al, 2008)。Vr和Z在水平方向上位于模式格点,垂直方向上位于雷达体扫锥面,即相当于在垂直方向上保留了雷达观测的仰角信息,在水平方向上进行从雷达观测到模式网格的水平插值。自模拟20 min开始,每5 min生成一个体扫数据。同化时使用的观测数据首先通过式(10)和(11)由观测时刻的真实场计算得到,然后加入随机扰动模拟有误差的观测。随机扰动满足均值为0的高斯分布,径向风和反射率的误差均方差分别为2 m/s和2 dBz。
为了全面考察不同位置雷达探测的信息对集合卡尔曼滤波同化的影响,共设置了8个模拟雷达。雷达距离地面10 m,水平位置相对于模拟区域的位置如图 2所示。可见各个雷达相对于模拟区域中心位置的距离为70—80 km,是比较理想的雷达观测距离。同时8个雷达环绕模拟区域均匀分布,意味着能够考察风暴出现在雷达的任意象限时集合卡尔曼滤波同化雷达资料的能力以及对随后预报的效果。
3.3 同化系统设置和试验设计集合卡尔曼滤波同化循环开始于模式时间第25 min,结束于第90 min,每5 min同化一次。为研究不同位置雷达观测对短时同化后确定性预报的影响,所有同化试验均从第4次同化循环(第40 min)的分析集合平均场开始进行50 min的确定性预报。40个初始集合成员通过在第20 min的水平均一环境场(见3.1节)上加入位温θ和水汽混合比qv的随机扰动得到。这些随机扰动满足均值为0的高斯分布,为了使扰动激发出的风场扰动特征与热泡刚刚触发的对流风场特征类似,文中仅对位温和水汽混合比加上扰动,扰动的均方差分别为3 K和0.5 g/kg,并且扰动仅加在2 km内有大于10 dBz的反射率观测的网格上,这样可以避免在对流区域外引入虚假的对流。同化系统的协方差膨胀方案采用松弛膨胀方案(Zhang et al, 2004),预报扰动和分析扰动的权重相等。此外,为了防止出现离散度不足的情况,还采用了额外的协方差膨胀方案,即当位温的分析离散度小于2.0 K时,将其放大至2.0 K(Wang et al, 2013a)。空间局地化采用Gaspari等(1999)的5阶距离相关函数,水平局地化距离为8 km,垂直局地化距离为4 km。
相对于传统的集合卡尔曼滤波,迭代集合均方根滤波增加了3处需要设置的参数:时间局地化系数、迭代参数、更新tn-1*时刻背景场的同化参数。时间局地化的目的是使tn时刻的观测对tn-1*时刻背景场的影响随时间间隔的增大而减小,文中设为6 min,略长于同化循环的间隔;迭代参数包括两部分,迭代步数和迭代时间间隔,分别设为3步和3 min,即每个同化循环中进行3次迭代,并且,tn-tn-1*=3 min;tn-1*时刻背景场的更新参数与tn时刻的基本相同,只是额外协方差膨胀方案中的2.0 K减为1.8 K,以减小膨胀后给模式积分带来的影响。
同化试验分为3组:集合卡尔曼滤波同化单雷达试验;集合卡尔曼滤波同化双雷达试验;对比试验(即由迭代集合均方根滤波分别同化第1组的单雷达和第2组的双雷达)。已有的反演和同化研究结果表明,双雷达分析的精度更高,尤其是在集合卡尔曼滤波早期的同化循环中(许小永等,2004;王世璋等,2013)。因此,设计此试验可以检验是否会因为双雷达的使用而减弱雷达方位对集合卡尔曼滤波同化的影响。单雷达试验的名称及试验使用的雷达由表 1列出,命名规则是“radar”加上雷达所在网格的方位。双雷达试验的命名规则与单雷达试验相同,只是增加了第2个雷达的网格方位。仅选取两个雷达方位与区域中心连线成90°夹角的双雷达组合,该组合也是双雷达反演中比较理想的雷达配置方式,组合方式如表 2。
试验名称 | x方向网格序号 | y方向网格序号 | 相对区域中心的位置 |
radar_E | 50 | 15 | 正东 |
radar_N | 15 | 50 | 正北 |
radar_W | -20 | 15 | 正西 |
radar_S | 15 | -20 | 正南 |
radar_NW | -15 | 45 | 西北 |
radar_SW | -15 | -15 | 西南 |
radar_NE | 45 | 45 | 东北 |
radar_SE | 45 | -15 | 东南 |
试验名称 | 雷达1 | 雷达2 |
radar_E_S | 正东 | 正南 |
radar_N_E | 正北 | 正东 |
radar_W_N | 正西 | 正北 |
radar_S_W | 正南 | 正西 |
radar_NW_NE | 西北 | 东北 |
radar_SW_NW | 西南 | 西北 |
radar_NE_SE | 东北 | 东南 |
radar_SE_SW | 东南 | 西南 |
根据Tong等(2005)、Xue等(2006)在对流尺度集合卡尔曼滤波同化中的研究,文中只计算真实场反射率大于10 dBz格点上的均方根误差(RMSE)。本节将通过误差相关系数分析雷达位置造成集合卡尔曼滤波分析差异和预报差异的原因。误差相关系数与背景误差协方差有一致的正、负符号,同时集合卡尔曼滤波的分析效果依赖于背景误差协方差,因此,对误差相关系数的分析有助于寻找以上问题的原因和解决方案。
4.1 单雷达观测的位置对同化的影响图 3给出了各个单雷达试验的风场和温度场的分析误差随同化次数的变化。在第1次同化循环中,试验radar_N和radar_S的u分量的分析误差最大,接近7 m/s,而radar_E和radar_W的u分量的分析误差最小,小于5 m/s,其他试验的u分量的分析误差相差不大,约为6 m/s。对于v分量,同化radar_N和radar_S试验的分析误差最小,约为4.5 m/s,而其他试验的v分量的分析误差均大于5 m/s。从首次同化的结果来看,处于模式区域对角线的雷达,即西北、东北、西南和东南的雷达的u和v分量的分析误差水平基本处于所有试验的中间,相互差别也较小,而处于正东、正南、正北和正西的雷达的分析误差则差异较大,其中,正南、正北雷达的v分量的分析误差最小,而u分量的分析误差最大,正东和正西雷达的效果则正好相反。该特征与Sun等(1994)和许小永等(2004)的反演结果相似,即在雷达观测径向投影较大的风速分量的分析效果最好,而与雷达径向垂直的分量的分析效果较差。w分量垂直于水平面而θ是标量,因此,这两个分量在首次分析中的误差受雷达位置的影响较小。在进行4次同化循环后,同化各方位雷达观测都能够显著降低各个变量的分析误差,但是不同雷达分析误差的水平仍有较为明显的差距。雷达观测位于模拟区域对角线方向上时,同化效果相对较好,而处于正南、正北方向时同化效果较差。
当最后一次同化循环(第14次)完成时,所有试验各变量的分析误差均达到同一水平:各试验u分量的分析误差约为2 m/s,最大偏差小于0.4 m/s;v分量的分析误差约为1.7 m/s,最大偏差小于0.3 m/s;w分量的误差约为1.5 m/s,最大偏差小于0.3 m/s;位温θ的误差偏差约为0.2 K。该结果与许小永等(2004)的结果一致,即雷达位置对集合卡尔曼滤波(包括集合均方根滤波)同化的影响主要出现在同化初期,在进行10多次同化循环后该影响即基本消失。
图 4给出了从第4次同化循环的分析场起报的确定性预报的误差。至第90 min预报结束,u分量的预报误差较大的4个试验分别是radar_S、radar_N、radar_W和radar_SW,均超过6.6 m/s,而在u分量的预报误差最小的radar_NW试验中,仅为4.6 m/s;v分量的误差排序和u分量基本一致,效果最好(radar_NW)与最差(radar_N)的预报误差相差2.5 m/s。对于w和θ,radar_S和radar_N的预报误差显著大于其他试验。上述结果表明,同化不同位置的雷达观测对短时同化后的确定性预报具有较明显的影响。
图 4b中radar_S和radar_N试验v分量尽管在同化结束时分析误差最小,但是在随后的确定性预报中,误差快速增大并在预报结束时较其他试验的预报误差大。该现象表明,u分量的误差在这个理想个例中对误差发展的影响超过v分量。考虑到环境场中高层的风向基本以西风为主且风暴基本向东移动,因此,同化对该方向观测的雷达径向风数据有可能不能在较短时间内得到较为理想的分析结果,这与Sun等(1994)的结论基本一致。此外,值得注意的是,尽管radar_SW试验的u、v、w和θ的分析误差均较小,但在最终的预报误差方面,radar_SW试验没有那些在第4次同化循环中分析误差比其大的试验表现好。这可能是由于误差的非线性增长所致。尽管该现象的存在给集合卡尔曼滤波同化和预报的效果带来了一定的不确定性,但分析误差较大时,预报误差通常也相对较大,例如radar_S、radar_N试验。
4.2 双雷达观测的位置对同化的影响图 5给出了多种组合的双雷达观测试验同化窗口内分析误差的变化。图 5具有两个明显的特征:一是早期同化循环中各个试验的分析误差的差异较大(u分量试验结果最明显),可分为两组;二是各个试验的最终分析误差的差异很小,小于第1组单雷达试验的差异。这两个特征表明,即使在双雷达条件下,雷达位置对同化初期的效果仍然有较为明显的影响,而在十几次同化循环之后,雷达位置的影响基本消失。对u分量的结果进行归类可以发现,误差较大的4个试验都是使用模拟区域非对角线上的雷达,即正东、正西、正南、正北雷达之间的组合。考虑到分析误差通常小于预报误差,因此,分析误差的快速增大通常对应预报误差的快速增大,而快速增大的预报误差通常是由于分析场不平衡所致。类似的情况在w和θ的分析误差变化上也可以看到,尤其是w,分析误差在第2次同化循环中急剧增大(超过4 m/s),最大误差超过10 m/s,增加幅度超过了单雷达试验对应的最大幅度(2 m/s)。导致不平衡的可能原因将在下节中讨论。总体而言,并不是所有的与观测对象连线成90°夹角的双雷达组合都能在短时同化中得到理想的分析结果。如果双雷达都处在模拟区域的对角线上,并与模拟区域中心的连线成90°夹角,则集合卡尔曼滤波同化循环早期的分析效果都比较理想。而当双雷达径向方向分别与u和v平行时,集合卡尔曼滤波早期同化的效果则相对不理想。考虑到图 1所示的西风环境气流,上述结论也可以表述为当双雷达与风暴连线成90°并且连线与环境气流夹角大致成45°时,集合卡尔曼滤波的同化分析效果最好。
从多种组合的双雷达短时同化对随后确定性预报的误差影响(图 6)可以看出,对于u分量,分析误差较大的4个试验在第90 min的预报误差均较大(超过5 m/s),而分析误差较小的4个试验,除1个试验的预报误差为5 m/s外,剩余3个试验的预报误差均小于4.5 m/s,其中,radar_NW_NE试验的效果最好,预报误差小于4 m/s。在第90 min各试验v分量误差的排序与u分量的情况类似,只是v预报误差的差异更为明显。几乎相同的预报误差的排序情况也可以在w分量和θ上看到,其中,各个试验的w分量在第90 min的预报误差的最大差值达2 m/s,θ的最大差值达1 K。上述结果表明,即使使用成90°夹角的双雷达进行同化,雷达的位置仍然会对短时同化及其随后的确定性预报产生显著的影响。
4.3 背景误差协方差分析鉴于无论是单雷达还是双雷达同化,雷达位置对集合卡尔曼滤波初期的同化效果都会产生不可忽略的影响,因此有必要探寻其原因,为解决该问题提供依据。为了便于研究,本节将统一使用试验radar_N的模式空间背景误差协方差,对不同位置雷达观测的误差相关空间结构进行分析。从式(7)可以看出,当模式空间的背景误差协方差固定时,同化效果仅与观测及观测算子有关。因此,本节使用同一个模式空间背景误差协方差与不同雷达的观测计算相关系数,这样可以仅考虑雷达位置对同化结果的影响,排除初值的影响。“正北方向雷达”和“正西方向雷达”分别指表 1中雷达位置是“正北”和“正西”的雷达。需要指出的是,不同模式变量是存在相关的,即使是u、v和w这3个互相垂直的变量。
图 7给出了试验radar_N在第4次同化循环中u和w集合平均的预报误差,及其与不同位置的雷达在同一点(红色实心三角处)上的观测的相关系数。由图 7a可见,如果使用来自“正北方向雷达”的观测,预报的径向风与径向风的真值的差值(预报减真值)为正值,意味着如果要通过这一增量修正红色圆圈内正的u的偏差(预报值大于真值),需要正的相关系数。就这一点而言,同化“正北方向雷达”的观测是能够修正u分量的误差的。但从u分量误差的空间分布来看,观测所在点的误差较小,其东北和西北方向的误差较大,而相关系数却从观测点向四周减小,因此,同化这一点的观测并不能合理地修正u的偏差。当这个观测来自“正西方向雷达”时,径向风的差值为负(图 7b),意味着需要负(正)的相关来修正u分量误差为正(负)值的区域。这一相关系数的要求在图 7b中得到了很好的满足。同时,由图 7b还可见,相关系数的中心在空间位置上基本对应于实际误差的中心,因此,同化来自正西方向的雷达观测能够有效修正u分量的误差。
垂直速度w在雷达径向上的投影相对u和v更小,因此,其同化分析效果对背景误差协方差的依赖程度更大。由图 7c可见,径向风的差值为正,因此要修正圆圈内较大的w的负差值需要负的相关系数,但是该区域内的相关系数为正,因此同化来自“正北方向雷达”的径向风并不能修正其周围的w分量的误差。此时如果使用来自“正西方向雷达”的观测,则圆圈内是负径向风、负w差值和正相关系数的组合(图 7d),且相关系数大值中心与误差大值中心相匹配,因此同化该雷达径向风观测能够有效地修正其影响区域内的w的误差。
综合图 7的结果可知,对于同一个模拟空间的背景误差协方差,即在排除初值影响的情况下,不同位置的雷达仍会导致同化时计算的观测空间误差与模式空间误差的协方差的合理性存在较大差别。由来自合适位置的雷达的观测计算出的协方差合理性强,而由来自不理想位置的雷达的观测计算出的协方差合理性较弱。正是此原因导致了同化初期使用不同位置雷达时集合卡尔曼滤波分析效果存在较大的差异。联合图 7c和d能在一定程度上解释图 5中使用正西、正北雷达组合分析得到的w分量在第2次同化循环中快速增长的可能原因:在同化了正西雷达得到合理的w增量(类似图 7c)后,使用正北雷达的分析又引入了错误方向的增量(类似图 7d),这样使用两个雷达对w的分析不仅没有只使用正西雷达的效果好,而且还破坏了变量之间的协调性。同理,图 5c中使用radar_N_E在第2次同化循环时会使w误差快速增长。因此可以认为,当同化两个不同位置的雷达导致出现相互矛盾的分析增量时,有可能在分析场中引入较大的不平衡,使得随后的预报误差快速增大。
由于不同位置的协方差和观测增量关系并不相同,为了确认上述结果不是个别现象,进一步检查了位于图 7观测位置垂直上方2 km的观测(图 8a、b)以及位于东北方向约2 km的观测(图 8c、d)。其中,位于图 7观测点上方2 km的径向风数据对应的雷达仰角更高,径向风包含的垂直运动的分量更大,因此,对w的分析效果应当更合理。图 8a中正北方向的雷达得到的增量大于0(红色三角形处),而该点的w的真值小于预报场的数值,需要一个小于0的协方差以修正图 8a中同化半径内的较大负值误差。但是通过正北雷达先验观测和模拟空间集合算出的协方差在这一区域大于0,同化这一位置上来自正北雷达的观测不能使误差减小。而在图 8b中,这一误差的空间结构得到了误差协方差的较好的代表,负观测的增量在负偏差区的相关系数为正,最大相关系数接近0.6且相关系数中心位于负偏差中心,因此,同化后能够减小w的预报误差。虽然图 8c与图 7观测点位置相差不大,但是相关系数的空间分布出现了显著差异。可见图 8c观测点周围的相关系数较小,并不能有效修正其周围较大的负偏差。而在图 8d中,正的观测增量和较大的负相关系数出现在同一个负偏差区,同化后该区域的误差能够显著减小。使用以上相同的方式检查了其他位置的观测增量与背景误差协方差和误差场的分布特征。结果表明,对于大部分的观测点,使用正西方向的雷达能够得到比正北雷达更合理的背景误差协方差以及分析增量。图 3c中使用正西雷达得到的对流区平均w分析误差小于使用正北雷达的分析误差也验证了这一结果。
根据图 3结果可知,随着同化次数的增多,不同试验间的误差水平逐渐趋于一致。因此,考察中后期同化循环的误差协方差结构有助于探寻降低同化初期雷达位置影响的办法。由图 9a可见,经过7次同化循环后,雷达径向风预报误差与u分量误差的相关系数空间分布已比第4次同化循环中的相应结构更加合理。而图 9c显示,此时误差相关系数的大值区与实际误差的大值区相匹配,使用“正北方向雷达”能够有效修正观测附近w的误差。由图 9b、d可见,通过来自“正西方向雷达”的观测计算的相关系数,在空间分布上与真实误差基本一致,因此,同化来自“正西方向雷达”的观测仍能有效地修正模式空间的误差。上述结果表明,各个试验的中后期同化循环的误差趋于一致的原因是,背景误差协方差得到了显著的改进。
根据集合卡尔曼滤波的理论,其背景误差协方差是流依赖的,因此,误差的空间结构是由其所依赖的“流”来决定的。随着同化次数的增多,分析场中风暴的结构和位置都得到了明显的改进,因此其驱动出的扰动结构也就合理。基于上述分析可以推断,既然迭代集合均方根滤波具有在较短时间内建立更接近真值的风暴的可能,那么在前几次同化循环中迭代集合均方根滤波也有可能减弱雷达位置对集合卡尔曼滤波短时同化及随后预报的影响。
4.4 雷达位置对迭代集合均方根滤波同化的影响为避免重复,本节只考察迭代集合均方根滤波试验与对应集合卡尔曼滤波试验的差值。为了方便比较,用iEnSRF-2和EnSRF-2表示第2个循环(对应模式时间30 min)时迭代集合均方根滤波和集合卡尔曼滤波试验的结果,以此类推。由图 10a可见,使用正北雷达观测时,iEnSRF-4比EnSRF-4试验u分量的分析误差小2 m/s,而使用正南雷达得到的分析误差iEnSRF-4比EnSRF-4约小1 m/s。迭代集合均方根滤波方法会提高大部分单雷达试验u分量的分析。在随后的确定性预报中,第90 min时大部分iEnSRF-14预报误差小于EnSRF-14,尤其是对于图 4a预报误差较大的几个试验,iEnSRF-14比EnSRF-14改进超过1 m/s,使用radar_W和radar_SW时,改进则超过2 m/s。在使用西北和东南雷达时,iEnSRF-14试验u分量的预报误差略大于EnSRF-14,但是幅度也在0.6 m/s以内。考虑到radar_NW和radar_SE试验在图 4a中是误差最小的两个试验,因此,迭代集合均方根滤波方法未能对这两个试验的预报效果产生改进也是可以接受的。对于w分量(图 10c),iEnSRF-14在使用除西北雷达以外的全部单雷达试验,都得到了相对EnSRF-14更小的预报误差,且有3个试验预报误差的改进超过2 m/s。而在使用西北雷达时,iEnSRF-14的预报误差也仅增长约0.6 m/s。
由图 10b、d可见,使用双雷达时,迭代集合均方根滤波同样能够改进分析和预报。所有双雷达观测组合试验中,u分量iEnSRF-4分析误差均小于EnSRF-4的分析误差,其中最大的改进超过2 m/s (radar_S_W),最小的改进也超过0.8 m/s。在随后的预报中,多数iEnSRF-14试验比对应的EnSRF-14试验预报误差小,其中改进幅度最大的接近3 m/s(对应图 6中u预报误差最大的试验radar_E_S)。使用迭代集合均方根滤波方法没有改进的试验radar_NW_NE和radar_SE_SW对应集合卡尔曼滤波试验中预报误差最小的两个,但是iEnSRF-14比EnSRF-14的预报误差也仅大0.6 m/s左右。对于w分量,iEnSRF-2显著降低了EnSRF-2的分析误差。这一改进使得第4次同化循环时,iEnSRF-4的w分量误差比EnSRF-4对应的误差小1—2 m/s。至第90 min预报结束,对于w分量的预报几乎是全部iEnSRF-14试验都取得了相对EnSRF-14试验更小的误差,其中radar_E_S和radar_W_N试验的预报误差减小超过2 m/s。
为了确认迭代集合均方根滤波分析误差的减小是由于背景误差协方差的改进,对其协方差结构也进行了检查。由于迭代集合均方根滤波需要进行迭代,无法使用如前一小节所述方法进行对比,因此,为公平起见,使用第1次同化循环的结果进行对比。从图 11中可以看出,使用正北雷达的集合卡尔曼滤波试验得到的观测增量为负值,而周围误差也为负值,意味着需要正的相关系数才能修正这一偏差,但是其周围相关系数为负,因此使用这一观测的结果是进一步加大了观测附近w的偏差。而迭代集合均方根滤波多次迭代后,背景误差协方差的结构得到了明显的改进。首先,在观测点附近,相关系数的中心与误差中心的匹配程度显著提高,尤其是观测点左侧的正偏差中心处。其次,相关系数的正、负分布也更为合理:负的观测增量在正的相关系数区(观测点右上方)能够修正该区域的负偏差,而在负相关系数区(观测点左侧)能够修正正偏差。因此,背景误差协方差的改进对迭代集合均方根滤波的分析有正贡献。
上述结果表明,使用迭代集合均方根滤波能够在一定程度上减小雷达位置对集合卡尔曼滤波前几次同化分析的影响,也表明第4.3节中关于背景误差协方差的分析是合理的。
5 结论基于已有研究,讨论了雷达位置对对流尺度集合卡尔曼滤波同化的影响,并通过背景误差相关系数分析了造成该影响的可能原因,在此基础上尝试使用迭代集合均方根滤波减小短时集合卡尔曼滤波同化中由于雷达位置造成的同化和预报的不确定性,得到以下结论:
(1) 集合卡尔曼滤波雷达资料直接同化会遇到与雷达观测反演类似的问题,即当雷达径向与所观测的气流方向垂直时难以估计该方向气流的信息。这个问题主要影响集合卡尔曼滤波前期的同化循环,且该影响不可忽略。
(2) 造成集合卡尔曼滤波同化初期对雷达位置较为敏感的一个可能原因是,同化初期背景场偏离真值较多,不能合理地驱动出代表真实误差的误差协方差结构。在此背景误差协方差条件下,使用部分位置的雷达能够得到较理想的分析效果,而使用其他的雷达则难以在短时间内得到合理的结果。
(3) 能够在集合卡尔曼滤波同化初期取得较理想的分析结果的雷达观测,其径向通常与气流的方向成45°夹角(在本例中对应模拟区域的4个对角线方向),而同化初期不理想分析结果的雷达的径向通常垂直于气流方向(在本例中对应正南和正北方向)。
(4) 双雷达集合卡尔曼滤波同化并不一定有利于集合卡尔曼滤波的短时同化,因为当背景场不是很准确时,同化两部雷达的数据可能会向分析场引入互相矛盾的资料,从而导致分析和预报效果变差。通常这种情况会出现在一部雷达与观测对象的连线垂直于环境气流而另一部雷达的相应连线与之垂直的情况下。
(5) 迭代集合均方根滤波能够在一定程度上改进集合卡尔曼滤波同化初期因雷达位置而造成的分析效果较差的情况,并可以改善短时同化后确定性预报的效果。
(6) 当雷达与风暴的连线和环境气流成90°夹角或双雷达中有一个雷达与环境气流成90°时,可以使用迭代集合均方根滤波提高分析效果。而在其他情况,可以使用标准集合均方根滤波方案减小计算代价。
本研究基于理想个例进行讨论,且仅考虑超级单体风暴一种类型,其结论是否适用于其他强对流天气仍需作进一步研究。此外,受篇幅所限,文中主要就雷达方位对同化的影响进行了讨论,而距模拟区域中心不同距离的雷达资料对集合卡尔曼滤波同化有哪些影响尚需进一步研究。
陈杰, 闵锦忠, 王世璋, 等. 2012. WRF-EnSRF系统同化多普勒雷达资料在多类型强对流天气过程的数值试验. 大气科学学报 , 35 (6) : 720–729. Chen J, Min J Z, Wang S Z, et al. 2012. A numerical experiment on WRF-EnSRF for assimilation of Doppler Radar data in multicase strong convective weather processes. Trans Atmos Sci , 35 (6) : 720–729. (in Chinese) |
闵锦忠, 王世璋, 陈杰, 等. 2012. 迭代EnSRF方案设计及在Lorenz96模式下的检验. 大气科学 , 36 (5) : 889–900. Min J Z, Wang S Z, Chen J, et al. 2012. The implementation and test of iterative EnSRF with Lorenz96 model. Chinese J Atmos Sci , 36 (5) : 889–900. (in Chinese) |
闵锦忠, 毕坤, 陈耀登, 等. 2013. 基于物理约束扰动的EnSRF雷达资料同化. 大气科学学报 , 36 (2) : 129–138. Min J Z, Bi K, Chen Y D, et al. 2013. Assimilation of Doppler radar data with EnSRF based on physical constraint perturbation. Trans Atmos Sci , 36 (2) : 129–138. (in Chinese) |
王世璋, 闵锦忠, 陈杰, 等. 2013. 迭代集合平方根滤波在风暴尺度资料同化中的应用. 大气科学 , 37 (3) : 563–578. Wang S Z, Min J Z, Chen J, et al. 2013. Application of iterative ensemble square-root filter in storm scale data assimilation. Chinese J Atmos Sci , 37 (3) : 563–578. (in Chinese) |
许小永, 郑国光, 刘黎平. 2004. 多普勒雷达资料4DVAR同化反演的模拟研究. 气象学报 , 62 (4) : 410–422. Xu X Y, Zheng G G, Liu L P. 2004. Dynamical and microphysical retrieval from simulated Doppler radar observations using the 4DVAR assimilation technique. Acta Meteor Sinica , 62 (4) : 410–422. (in Chinese) |
许小永, 刘黎平, 郑国光. 2006. 集合卡尔曼滤波同化多普勒雷达资料的数值试验. 大气科学 , 30 (4) : 712–728. Xu X Y, Liu L P, Zheng G G. 2006. Numerical experiment of assimilation of Doppler radar data with an ensemble Kalman filter. Chinese J Atmos Sci , 30 (4) : 712–728. (in Chinese) |
郑永光, 张小玲, 周庆亮, 等. 2010. 强对流天气短时临近预报业务技术进展与挑战. 气象 , 36 (7) : 33–42. Zheng Y G, Zhang X L, Zhou Q L, et al. 2010. Review on severe convective weather short-term forecasting and nowcasting. Meteor Mon , 36 (7) : 33–42. (in Chinese) |
Anderson J L, Hoar T, Raeder K, et al. 2009. The data assimilation research testbed: A community facility. Bull Amer Meteor Soc , 90 (9) : 1283–1296. DOI:10.1175/2009BAMS2618.1 |
Anderson J L. 2001. An ensemble adjustment Kalman filter for data assimilation. Mon Wea Rev , 129 (12) : 2884–2903. DOI:10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2 |
Bishop C H, Etherton B J, Majumdar S J. 2001. Adaptive sampling with the ensemble transform Kalman filter.Part I: Theoretical aspects. Mon Wea Rev , 129 (3) : 420–436. DOI:10.1175/1520-0493(2001)129<0420:ASWTET>2.0.CO;2 |
Chang W Q, Chung K S, Fillion L, et al. 2014. Radar data assimilation in the Canadian high-resolution ensemble Kalman filter system: Performance and verification with real summer cases. Mon Wea Rev , 142 (6) : 2118–2138. DOI:10.1175/MWR-D-13-00291.1 |
Coniglio M C, Stensrud D J, Wicker L J. 2006. Effects of upper-level shear on the structure and maintenance of strong quasi-linear mesoscale convective systems. J Atmos Sci , 63 (4) : 1231–1252. DOI:10.1175/JAS3681.1 |
Dowell D C, Wicker L J. 2009. Additive noise for storm-scale ensemble data assimilation. J Atmos Ocean Tech , 26 (5) : 911–927. DOI:10.1175/2008JTECHA1156.1 |
Evensen G. 1994. Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J Geophys Res , 99 (C5) : 10143–10162. DOI:10.1029/94JC00572 |
Gaspari G, Cohn S E. 1999. Construction of correlation functions in two and three dimensions. Quart J Roy Meteor Soc , 125 (554) : 723–757. DOI:10.1002/(ISSN)1477-870X |
Hu M, Xue M. 2007a. Analysis and prediction of 8 May 2003 Oklahoma City tornadic thunderstorm and embedded tornado using ARPS with assimilation of WSR-88D radar data//22nd Conf Wea Anal Forecasting/18th Conf Num Wea Pred. Salt Lake City, Utah: Amer Meteor Soc |
Hu M, Xue M. 2007b. Implementation and evaluation of cloud analysis with WSR-88D reflectivity data for GSI and WRF-ARW. Geophys Res Let , 34 (7) : L07808. DOI:10.1029/2006GL028847 |
James M, Yvette R, Paul M, et al. 2014. An investigation of the Goshen County, Wyoming, tornadic supercell of 5 June 2009 using EnKF assimilation of mobile mesonet and radar observations collected during VORTEX2. Part I: Experiment design and verification of the EnKF analyses. Mon Wea Rev , 142 (2) : 530–554. |
Kalnay E. 2003. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge: Cambridge University Press : 341pp . |
Kalnay E, Yang S C. 2010. Accelerating the spin-up of ensemble Kalman filtering. Quart J Roy Meteor Soc , 136 (651) : 1644–1651. DOI:10.1002/qj.652 |
Kawabata T, Kuroda T, Seko H, et al. 2011. A cloud-resolving 4DVAR assimilation experiment for a local heavy rainfall event in the Tokyo metropolitan area. Mon Wea Rev , 139 (6) : 1911–1931. DOI:10.1175/2011MWR3428.1 |
Maybeck P S. 1979. Stochastic Models, Estimation and Control(Ⅰ). New York: Academic Press : 441pp . |
Putnam B J, Xue M, Jung Y, et al. 2014. The analysis and prediction of microphysical states and polarimetric radar variables in a mesoscale convective system using double-moment microphysics, multinetwork radar data, and the ensemble Kalman filter. Mon Wea Rev , 142 (1) : 141–162. DOI:10.1175/MWR-D-13-00042.1 |
Raeder K, Anderson J L, Collins N, et al. 2012. DART/CAM: An ensemble data assimilation system for CESM atmospheric models. J Climate , 25 (18) : 6304–6317. DOI:10.1175/JCLI-D-11-00395.1 |
Sakov P, Evensen G, Bertino L. 2010. Asynchronous data assimilation with the EnKF. Tellus A , 62 (1) : 24–29. DOI:10.1111/tea.2010.62.issue-1 |
Snook N, Xue M, Jung Y. 2015. Multiscale EnKF assimilation of radar and conventional observations and ensemble forecasting for a tornadic mesoscale convective system. Mon Wea Rev , 143 (4) : 1035–1057. DOI:10.1175/MWR-D-13-00262.1 |
Sobash R A, Stensrud D J. 2013. The impact of covariance localization for radar data on EnKF analyses of a developing MCS: Observing system simulation experiments. Mon Wea Rev , 141 (11) : 3691–3709. DOI:10.1175/MWR-D-12-00203.1 |
Sun J Z, Crook A. 1994. Wind and thermodynamic retrieval from single-Doppler measurements of a gust front observed during Phoenix Ⅱ. Mon Wea Rev , 122 (6) : 1075–1091. DOI:10.1175/1520-0493(1994)122<1075:WATRFS>2.0.CO;2 |
Sun J Z, Crook N A. 1997. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part Ⅰ: Model development and simulated data experiments. J Atmos Sci , 54 (12) : 1642–1661. |
Tanamachi R L, Wicker L J, Dowell D C, et al. 2013. EnKF assimilation of high-resolution, mobile Doppler radar data of the 4 May 2007 Greensburg, Kansas, supercell into a numerical cloud model. Mon Wea Rev , 141 (2) : 625–648. DOI:10.1175/MWR-D-12-00099.1 |
Tong M Q, 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. DOI:10.1175/MWR2898.1 |
Wang S Z, Xue M, Schenkman A D, et al. 2013a. An iterative ensemble square root filter and tests with simulated radar data for storm-scale data assimilation. Quart J Roy Meteor Soc , 139 (676) : 1888–1903. DOI:10.1002/qj.2077 |
Wang S Z, Xue M, Min J Z. 2013b. A four-dimensional asynchronous ensemble square-root filter (4DEnSRF) algorithm and tests with simulated radar data. Quart J Roy Meteor Soc , 139 (672) : 805–819. DOI:10.1002/qj.v139.672 |
Weisman M L, Klemp J B. 1982. The dependence of numerically simulated convective storms on vertical wind shear and buoyancy. Mon Wea Rev , 110 (6) : 504–520. DOI:10.1175/1520-0493(1982)110<0504:TDONSC>2.0.CO;2 |
Whitaker J S, Hamill T M. 2002. Ensemble data assimilation without perturbed observations. Mon Wea Rev , 130 (7) : 1913–1924. DOI:10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2 |
Whitaker J S, Compo G P, Wei X, et al. 2004. Reanalysis without radiosondes using ensemble data assimilation. Mon Wea Rev , 132 (5) : 1190–1200. DOI:10.1175/1520-0493(2004)132<1190:RWRUED>2.0.CO;2 |
Xu Q, Lu H J, Gao S T. 2008. Time-expanded sampling for ensemble Kalman filter: Assimilation experiments with simulated radar observations. Mon Wea Rev , 136 (7) : 2651–2668. DOI:10.1175/2007MWR2185.1 |
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 Tech , 23 (1) : 46–66. DOI:10.1175/JTECH1835.1 |
Yussouf N, Stensrud D J. 2010. Impact of phased-array radar observations over a short assimilation period: Observing system simulation experiments using an ensemble Kalman filter. Mon Wea Rev , 138 (2) : 517–538. DOI:10.1175/2009MWR2925.1 |
Zhang F, Snyder C, Sun J Z. 2004. Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter. Mon Wea Rev , 132 (5) : 1238–1253. DOI:10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2 |