大气科学  2019, Vol. 43 Issue (2): 285-296   PDF    
X波段双线偏振多普勒雷达共极化差分相移资料的滤波方法研究
赵川鸿1, 周筠珺1,2, 肖辉3,4, 赵鹏国1, 张晓玲5, 胡天洁5     
1 成都信息工程大学大气科学学院高原大气与环境四川省重点实验室, 成都 610225
2 南京信息工程大学气象灾害预警预报与评估协同创新中心, 南京 210044
3 中国科学院大气物理研究所中国科学院云降水物理与强风暴重点实验室, 北京 100029
4 中国科学院大学, 北京 100049
5 北京市气象信息中心, 北京 100089
摘要: 为了提高X波段双线偏振多普勒雷达的共极化差分相移ΨC质量,利用714XDP-A型车载X波段双线偏振多普勒雷达在北京顺义的探测资料,构造了一个综合滤波方法。将其应用于滤波处理,得到结果:(1)综合小波去噪滤波方法在滤波处理及差分传播相移率KDP的估算方面相较于中值滤波、综合滑动平均、综合卡尔曼滤波及综合有限脉冲响应(FIR)滤波都是最优秀的,其能较好地保留原始ΨC平均趋势及细节特征(偏差基本在±2°可信范围内,部分区域由于对原始ΨC细节特征的保留会略有超出),避免探测误差零值对滤波产生的影响(误差零值订正值回归平均趋势水平,偏差在±2°可信范围内),对原始ΨC距离廓线的平滑程度达到了94%~97%[平滑了超出±1.5° Gate-1(Gate表示雷达距离库)的波动,对±1.5° Gate-1内波动保留];(2)估算的KDP整体在±5° km-1边界阈值内,避免了探测误差零值、环境噪声和δKDP估算产生的误差影响,KDP对较强降水指示明显(在回波强度30 dBZZH < 37 dBZ的层积云降水中KDP普遍大于0.25° km-1,峰值可达到4.3° km-1;在回波强度30 dBZZH < 50 dBZ的积雨云降水中KDP普遍大于0.1° km-1,峰值可达到1.95° km-1)。本文工作对于降水估测和水成物粒子识别有着改善效果,这对雷暴天气预警预报具备重要意义。
关键词: 共极化差分相移    滤波    X波段双线偏振多普勒雷达    
A Study of Method for Filtering Copolar Differential Phase of X-Band Dual-Polarimetric Doppler Weather Radar
ZHAO Chuanhong1, ZHOU Yunjun1,2, XIAO Hui3,4, ZHAO Pengguo1, ZHANG Xiaoling5, HU Tianjie5     
1 Plateau Atmosphere and Environment Key Laboratory of Sichuan Province, College of Atmospheric Science, Chengdu University of Information & Technology, Chengdu 610225
2 State Key Laboratory of Meteorological Disaster of Ministry of Education, Nanjing University of Information Science & Technology, Nanjing 210044
3 Key Laboratory of Cloud-Precipitation Physics and Severe Storms, Institute of Atmospheric Physics, Chinese Academy of Sciences(IAP/CAS), Beijing 100029
4 University of Chinese Academy of Sciences, Beijing 100049
5 Meteorological Information Center of Beijing, Beijing 100089
Abstract: Measurements of the 714XDP-A X-band dual-polarimetric Doppler radar located at Shunyi of Beijing are used to develop an integrated filter method that can improve the quality of the copolar differential phase ΨC measured by X-band dual-polarimetric Doppler radar. Application of the integrated filter method shows that:(1) The integrated wavelet de-noising method is the best and better than the median, integrated mean, integrated Kalman and integrated FIR (Finite Impulse Response) filters. It retains the mean trend of raw ΨC range profiles and details[deviations are within ±2°, which indicates the results are reliable, although the values may slightly exceed the range in several gates (the numbers of radar range profile)]. Furthermore, the influence of zero value caused by measurement errors on the filtered results is avoided (the corrected value for zero value in the measurements returns to the level of mean trend with deviations within ±2°, which is regarded as credible). In addition, the fluctuation index of ΨC range profiles filtered by integrated wavelet de-noising can be reduced by 94%-97% compared to raw ΨC range profiles, and the smoothing effect is remarkable (smoothing the fluctuations exceeding ±1.5° Gate-1, and retaining the fluctuations inside ±1.5°Gate-1). (2) The estimated specific differential propagation phase (KDP) are within the range of ±5° km-1, and influences caused by zero errors in the measurements and δ and environment noises are avoided. In addition, KDP can well indicate strong precipitation (for precipitation in stratocumulus clouds, radar reflectivity is usually within 30 dBZZH < 37 dBZ, and KDP generally is greater than 0.25° km-1 with the maximum of 4.3° km-1; for precipitation in cumulonimbus clouds, radar reflectivity usually is within 30 dBZZH < 50 dBZ, and KDP is generally greater than 0.1° km-1 with the maximum of 1.95° km-1). This study is beneficial to rainfall estimation and hydrometeors classification, which is important for thunderstorm weather forecast.
Keywords: Copolar differential phase    Filter    X-band dual-polarimetric Doppler radar    
1 引言

双线偏振雷达探测降水系统获得的水平反射率因子ZH、差分反射率ZDR、差分传播相移ΦDP、差分传播相移率KDP和共极化相关系数ρHV等参数有利于降水系统的粒子相态、滴谱分布(DSD)等云微物理结构变化规律研究,从而明显提高雷达对于降水估测以及降水粒子相态识别的能力(Seliga and Bringi, 1976Zrnić and Ryzhkov,1999Doviak et al., 2000Straka et al., 2000Bringi and Chandrasekar, 2001Park et al., 2009, 2015Chang et al., 2014)。

为了提高降水估测水平,首要条件是提升差分传播相移率KDP的估算质量,因为通过建模研究证明KDP(或者综合KDPZDR)对于降水估算来说是优越的估算量(Sachidananda and Zrnić,1986Bringi et al., 1990Jameson and Caylor, 1994Ryzhkov et al., 2003)。KDP是差分传播相移ΦDP的径向距离廓线平均趋势的斜率,而实际探测过程中,雷达探测到的共极化差分相移ΨC包含了差分传播相移ΦDP和差分后向散射相移δMueller(1984)提出算法计算ΨC

${\mathit{\Psi }_{\rm C}} = \arg (\left\langle {{S_{{\rm{vv}}}}S_{{\rm{hh}}}^*} \right\rangle) + 2(k_{\rm{h}}^{\rm{r}} - k_{\rm{v}}^{\rm{i}})r = \delta + {\mathit{\Phi }_{{\rm{DP}}}}$ , (1)

其中,“arg”表示取幅角计算,“*”表示取共轭,Svv是垂直散射幅度,Shh是水平散射幅度,kh是水平传播常数,kv是垂直传播常数,上标r、i分别代表实部和虚部,尖角括号代表整体平均,δ是差分后向散射相移(短波长雷达在探测较大水成物粒子时由于非瑞利散射导致δ的产生),ΦDP是差分传播相移。KDP通过有限差分估算(Hubbert and Bringi, 1995):

${K_{{\rm{DP}}}} = \frac{{{\mathit{\Phi }_{{\rm{DP}}}}({r_2}) - {\mathit{\Phi }_{{\rm{DP}}}}({r_1})}}{{2({r_2} - {r_1})}}$ , (2)

其中,r1r2是雷达径向距离。实际雷达探测资料中原始ΨC存在探测误差:(1)探测误差零值;(2)随距离的增加,径向距离廓线波动幅度也增大;(3)在雷达附近区域的ΨC资料质量较差。原始ΨC资料存在的探测误差与环境噪声(包括系统噪声、随机波动以及杂波等)和δ有关(Ryzhkov and Zrnić,1998Gosset,2004Ryzhkov,2007Gorgucci and Baldini, 2015)。Ryzhkov(2007)还指出,中尺度对流系统中,当非均匀波束填充存在,波束的展宽效应是造成ΦDP误差显著的首要原因(比起纯统计误差和δ造成的误差),并提出一个简单模型用以定量非均匀波束填充对ΦDP造成的误差:

${\rm{\Delta }}{\mathit{\Phi }_{{\rm{DP}}}} = 0.02{\mathit{\Omega }^2}(\frac{{{\rm{d}}{{\mathit{\Phi '}}_{{\rm{DP}}}}}}{{{\rm{d}}\theta }}\frac{{{\rm{d}}{Z_{{\rm{HV}}}}}}{{d\theta }} + \frac{{{\rm{d}}{{\mathit{\Phi '}}_{{\rm{DP}}}}}}{{{\rm{d}}\phi }}\frac{{{\rm{d}}{Z_{{\rm{HV}}}}}}{{{\rm{d}}\phi }})$ , (3)

其中,Ω表示单程3 dB天线模式宽度,θ表示仰角,ϕ表示方位角,ΔΦDP表示ΦDP误差。由公式(3)可见,当雷达波束宽度越大,将会导致ΦDP误差过大。想要精确地识别以及订正探测误差是不实际的,但是可以通过判断其大致的偏差范围(-2°~2°)确定资料是否可用于降水估测(Ryzhkov,2007)。因此,由于环境噪声和δ的存在,使得由ΨC近似估算的ΦDP存在误差,进而对KDP的估算造成严重的干扰,从而影响降水估测,所以对ΨC进行滤波处理是很有必要的。国外对于C波段以及X波段的双偏振雷达ΨC资料滤波处理主要应用Hubbert和Bringi(1995)提出的有限冲击响应(FIR)迭代滤波方法(Hubbert and Bringi, 1995Park et al., 2005, 2015Chang et al., 2014)。国内学者对于ΨC资料的滤波处理研究还相对较少,近年来,何宇翔等(2009)引进了卡尔曼滤波,杜牧云等(2012a, 2012b)利用信噪比和ρHV将回波信号分为较好、较差和差,并根据不同需求进行处理,形成了一套分类处理方法,Hu and Liu(2014)利用小波去噪对C波段双偏振雷达探测的台风和飑线两种降水系统的ΨC资料进行滤波处理,并与滑动平均、中值滤波、卡尔曼滤波以及FIR迭代滤波等滤波方法对比,指出利用小波滤波在有效平滑原始ΨC径向距离廓线的同时能较好地保留原始ΨC径向距离廓线平均趋势和细节特征,这些细节特征对于杂波的识别以及定位降水大值区有着很好的作用。

随着对降水估测质量以及粒子相态识别的要求越来越高,由环境噪声和δ导致的ΨC原始探测误差不容忽视,为了更好的估算KDP,研究应用于X波段双线偏振雷达的ΨC资料的滤波方法是很有必要的,这对于提高雷达对降水估测以及降水粒子相态识别的能力具备重要的意义。

2 资料和方法 2.1 数据来源

本文数据来自北京顺义(40.18°N,116.68°E)的714XDP-A型车载X波段双线偏振多普勒天气雷达。该雷达海拔高度为26.1 m,最大探测范围150 km,采用双发双收模式,可探测8个偏振参量,包含水平反射率因子(ZH)、垂直反射率因子(ZV)、多普勒径向速度(V)、速度谱宽(W)、差分反射率(ZDR)、共极化差分相移(ΨC)、差分传播相移率(KDP)以及共极化相关系数(ρHV),主要性能指标如表 1所示。

表 1 714XDP-A型X波段双线偏振多普勒天气雷达主要性能指标 Table 1 The main parameters of 714XDP-A X-band dual- polarimetric Doppler weather radar
2.2 数据预处理

对714XDP-A型车载X波段双线偏振多普勒天气雷达探测的ΨC资料先进行地物杂波抑制及退折叠处理(如有必要的情况下进行,高仰角3°~5°的探测资料被认为基本不受地物杂波影响)。地物杂波识别采用模糊逻辑识别(Cho et al., 2006),对于ΨC资料折叠区域的退折叠,沿用径向连续性检查的方法(肖艳姣等,2012)。因为X波段雷达存在的衰减现象,还需对ZH以及ZDR进行衰减订正,采用自适应约束算法(Park et al., 2005毕永恒等,2012)。以下资料都是经过预处理之后的资料。

2.3 原始ΨC探测误差

选取2015年8月7日19:50(北京时,下同)探测资料对原始ΨC探测误差进行说明,图 1ab中,ΨC红色矩形框内出现的零值区域(距离库数264~266),但是对应ZH距离廓线(图 1a)在此区域仍有有效回波,这种ZH不为零而ΨC出现的零星零值是探测误差零值。同时,该资料随着距离库数的增加,ΨC距离廓线的抖动幅度也逐渐增大,并且在距离库数6以内ΨC探测值极不稳定。其中在距离库数217~218的ΨC径向增量达到9°,考虑为明显的δ效应。利用中值滤波、滑动平均、小波去噪、卡尔曼滤波及FIR对该资料进行滤波处理,并选取方位角37°径向距离廓线(图 2)。

图 1 2015年8月7日19:50(北京时,下同)探测资料:(a)反射率因子ZH在方位角37°径向距离廓线;(b)共极化差分相移ΨC在方位角37°径向距离廓线(红色矩形框内是ΨC零值) Fig. 1 Measurements at 1950 BJT (Beijing time) 7 August 2015: (a) The radial range profile at azimuth 37° of reflectivity ZH; (b) the radial range profile at azimuth 37° of copolar differential phase ΨC (the values of ΨC inside the red rectangle are 0°)

图 2 2015年8月7日19:50 ΨC在方位角37°径向距离廓线:(a)原始资料;(b–f)对资料分别进行小波去噪、滑动平均、中值滤波、卡尔曼滤波和FIR滤波处理。红色虚线是原始ΨC径向距离廓线平均趋势 Fig. 2 The radial range profiles of ΨC at azimuth 37° at 1950 BJT 7 August 2015: (a) Raw data, (b–f) filtered by wavelet de-noising, running mean, median, Kalman, and FIR (Finite Impulse Response) filters. The red dashed line in each panel denotes the mean trend of raw ΨC range profile

图 2可见:(1)滑动平均、小波去噪、卡尔曼滤波以及FIR迭代滤波都对探测误差零值敏感,不能对其进行有效订正,尤其是FIR滤波后由于探测误差零值的影响造成距离库数266~282上偏离平均趋势严重,数据严重失真;(2)滑动平均、小波去噪以及卡尔曼滤波对δ(如距离库数217~218)和环境噪声(如距离库数220~225)特征进行了保留;(3)FIR迭代滤波对于原始ΨC径向距离廓线的平均趋势保留不如中值滤波、滑动平均、小波去噪和卡尔曼滤波,其在距离库数213~238段失真较为严重。综上,在利用滑动平均、小波去噪、卡尔曼滤波以及FIR迭代滤波时,需要对原始ΨC探测误差零值进行订正避免对滤波结果进行干扰造成KDP估算误差;对于滑动平均、小波去噪和卡尔曼滤波还需要对原始ΨCδ及环境噪声进行去噪处理。

2.4 方法

对于原始ΨCδ及环境噪声的去噪处理本文提出资料识别方法,对于原始ΨC探测误差零值的订正本文提出逐步订正方法。

(1)资料识别

Hubbert and Bringi(1995)指出差分传播相移ΦDP是距离累积量,即使在大雨情况下其增量也应小于8° km−1~10° km−1(根据雷达频率不同而有所不同),ΦDP是距离的平滑函数并且其径向距离廓线变化缓慢。因此从理论上讲,ΦDP的径向距离廓线整体上应该呈随距离的增大而增加的趋势(Hubbert and Bringi, 1995Ryzhkov and Zrnić,1998Wang and Chandrasekar, 2009),并且其径向距离增量在一定范围内。Ryzhkov(2007)指出,在X和C波段双偏振雷达实际应用中,1°~2°的统计波动范围被认为是ΦDP估算存在的典型边界。这说明由δ和环境噪声引起的不可取的波动超出±1°~2°范围,并且结合Hubbert and Bringi(1995)指出的ΦDP径向增量理论阈值应在8° km−1~10° km−1(对应本文1.2° Gate−1~1.5° Gate−1)。

因此本文对ΨC求取径向增量(∆ΨC),其单个径向距离库单元分辨率为0.15 km,利用模糊逻辑方法对径向增量(∆ΨC)进行识别,从而识别出δ和环境噪声引起的增量。∆ΨC计算公式如下:

$\Delta {\mathit{\Psi }_{\rm{C}}} = {\mathit{\Psi }_{{\rm{C}}(i + 1)}} - {\mathit{\Psi }_{{\rm{C}}(i)}}$ , (4)

其中,i代表径向距离库数。在计算∆ΨC前,由于本部雷达在距离库数6以内,雷达附近ΨC质量很差,其径向增量超过±50°,ΨC变化趋势不可信,就算其中参杂有降水回波信号,也无法参考,所以不考虑这部分不可信资料(杜牧云等,2012a)。

综合前人对于ΨC径向增量的统计研究建立了基于不对称梯形函数的特征参量∆ΨC的隶属函数模糊基。其中,不对称梯形函数:

${\rm{MF}}(x, {\rm{ }}x1, {\rm{ }}x2, {\rm{ }}x3, {\rm{ }}x4) = \left\{ {\begin{array}{*{20}{c}} {0, \;\;\;\;\;\;\;\;\;\;\;x < x1}\\ {\frac{{x - x1}}{{x2 - x1}}, {\rm{ }}x1 \le x \le x2}\\ {\frac{{x4 - x}}{{x4 - x3}}, {\rm{ }}x3 \le x \le x4}\\ {0.\;\;\;\;\;\;\;\;\;\;x4 < x} \end{array}} \right.$ (5)

隶属函数值MFtotal采用等权重求和集成,阈值定义为0.5,MFtotal≥0.5的ΨC资料为差资料(表示ΨC径向距离库增量超出±1.5°,受δ和噪声污染),0<MFtotal<0.5的ΨC资料为一般资料(表示ΨC径向距离库增量在±1.5°内,有可能存在轻微的δ和噪声污染),MFtotal=0的ΨC资料为较好资料(表示ΨC径向距离库增量在±1°内,可认为不存在δ和噪声污染)。

(2)逐步订正

将识别后的资料进行处理,重构ΨC数据。首先因为差资料(MFtotal≥0.5)存在δ以及环境噪声明显,因此滤除,对于一般资料(0<MFtotal<0.5),即使其中存在噪声也是轻微的,对于滤波结果不会有太大影响,为了尽可能保留原始ΨC数据特征,因此选取一般资料以及较好资料一起重构ΨC。利用识别后得到的较好资料和一般资料作为可信值,而后利用ZH识别资料序列中的零值是否为误差值(对比ZHΨC,如果在ZH不为零时,ΨC出现零值,则判定为误差值;如果两者都为零,则是合理零值),如果是误差零值则用向前相邻的可信值代替,合理零值则继续保留(将其称为逐步订正法),因为ΦDP的径向距离廓线变化是缓慢的,而且整体上应该呈随距离的增大而增加的趋势(Hubbert and Bringi, 1995Ryzhkov and Zrnić,1998Wang and Chandrasekar, 2009),因此逐步订正的结果是较为合理的并且近似满足ΦDP的理想变化趋势,同时,从理论上讲,这样的订正可以修正原始ΨC的误差零值,并且不改变合理零值。经过逐步订正法处理能简单且较为合理地订正原始ΨC资料探测误差零值,并以原始ΨC径向距离廓线平均趋势作为基准判定滤波后ΨC径向距离廓线与原始ΨC径向距离廓线平均趋势偏差是否在±2°可信范围内。

因此,利用滑动平均、小波去噪和卡尔曼滤波时先对原始ΨC进行资料识别和逐步订正(称为A类处理),而利用FIR迭代滤波时只需要对原始ΨC探测误差零值进行订正(称为B类处理),对于中值滤波不做处理,如表 2所示。

表 2 ΨC的综合处理方法 Table 2 The composite processing of ΨC
3 结果分析 3.1 滤波对比

由于层云降水和对流云降水的不同,本文分别选取了一次层积云降水过程和一次积雨云降水过程进行讨论。2015年8月7日19:50(仰角5°)雷达平面位置显示器(PPI)显示(图 3a)这是一次层积云降水过程,在方位角37°径向上主要为层云降水,混杂小范围积云对流降水,降水范围大致在距离库数93~311,其中在距离库数264~266存在探测误差零值,该径向上存在30 dBZZH<37 dBZ的较强降水区域[张培昌等(2000)指出,如果ZH≥30 dBZ,就可能产生较强的降水],其范围大致在距离库数198~259(29.7~38.85 km),云高大约在3 km左右。

图 3 2015年(a)8月7日19:50(仰角5°)反射率因子ZH(红线代表方位角37°径向)和(b)6月26日22:26(仰角3°)反射率因子ZH(红线代表方位角210°径向) Fig. 3 (a) The reflectivity ZH at the elevation angle 5° at 1950 BJT 7 August 2015 (the red solid line denotes the radial at azimuth 37°); (b) the reflectivity ZH at the elevation angle 3° at 2226 BJT 26 June 2015 (the red solid line denotes the radial at azimuth 210°)

同时为了探究对流降水系统,选取了2015年6月26日22:26(仰角3°)探测资料,如图 3b所示,在雷达站西南方向,方位角210°径向上存在雷暴单体,为积雨云降水,降水范围大致在距离库数533~690,其中距离库数586、648~649存在探测误差零值,该径向上存在回波强度为30 dBZZH<50 dBZ的较强降水区域,大致范围在距离库数623~682(93.45~102.3 km),其中心回波强度45 dBZZH<50 dBZ,所在高度在5.24 km左右,所在范围667~674。

针对降水区域应用中值滤波、综合滑动平均、综合小波去噪、综合卡尔曼滤波和综合FIR滤波五种方法对ΨC分别进行滤波处理(已经过A、B类综合处理;图 4),并统计五种方法滤波后ΨC径向距离廓线与原始ΨC径向距离廓线平均趋势的偏差(图 5)以及滤波后的ΨC径向波动指数FIX(Fluctuation Index;表 3),FIX算法如下:

图 4 2015年8月7日19:50时层积云(Sc)和2015年6月26日22:26时积雨云(Cb)中(a1、a2)原始ΨC、(b1、b2)综合小波去噪、(c1、c2)综合滑动平均、(d1、d2)中值滤波、(e1、e2)综合卡尔曼滤波和(f1、f2)综合FIR滤波(其中红色虚线是原始ΨC径向距离廓线平均趋势) Fig. 4 (a1, a2) The raw ΨC range profiles for Stratocumulus (Sc) (1950 BJT on 7 August 2015) and Cumulonimbus (Cb) (2226 BJT on 26 June 2015) respectively, and ΨC filtered by (b1, b2) integrated wavelet de-noising, (c1, c2) integrated running mean, (d1, d2) median, (e1, e2) integrated Kalman, and (f1, f2) integrated FIR filters (the red dashed line in each panel denotes the mean trend of raw ΨC range profile)

图 5 层积云(Sc)和积雨云(Cb)中(a1、a2)原始ΨC偏差统计、(b1、b2)综合小波去噪偏差统计、(c1、c2)综合滑动平均偏差统计、(d1、d2)中值滤波偏差统计、(e1、e2)综合卡尔曼滤波偏差统计、(f1、f2)综合FIR滤波偏差统计(其中红色实线表示±2°),时间与图 4一致 Fig. 5 Statistics of deviations for (a1, a2) raw ΨC range profiles at Sc and Cb respectively, ΨC filtered by (b1, b2) integrated wavelet de-noising, (c1, c2) integrated running mean, (d1, d2) median, (e1, e2) integrated Kalman, and (f1, f2) integrated FIR filters (the red solid lines in each panel denote ±2°), the time is the same as Fig. 4

表 3 波动指数FIX统计 Table 3 The statistics of FIX (fluctuation index)
${\rm{FIX}} = \frac{1}{{n - 1}}\sum {_{i = 2}^n\left| {\mathit{\Phi }_{{\rm{DP}}}^i - \mathit{\Phi }_{{\rm{DP}}}^{i - 1}} \right|} $ , (6)

其中,i表示距离库数,n表示总的距离库数,ΦDPΨC滤波后得到的近似值。

图 4图 5可见:(1)经过A、B类处理后的综合滤波方法较好地抑制了探测误差零值(图 4a1距离库数264~266,图 4a2距离库数586、648~649)及δ(如图 4a1中距离库数217~218)和环境噪声(如图 4a1中距离库数220~225)对滤波结果的影响。(2)综合卡尔曼滤波、综合小波去噪和综合滑动平均都对原始ΨC细节特征有所保留,其中综合卡尔曼滤波保留效果最好,综合小波去噪其次,综合滑动平均不及前两种方法,如图 5c1中综合滑动平均后ΨC在距离库数216~229、245~248略有超出±2°可信偏差范围是因为对原始ΨC细节特征的保留(对应原始ΨC距离库数220~224、244~248),值得指出的是距离库数262~268上超出±2°可信偏差范围是由于未能彻底避免探测误差零值(264~266)的干扰;图 5b1中综合小波去噪后ΨC在距离库数219~229、244~246略有超出±2°可信偏差范围也是由于对原始ΨC细节特征的保留,就其保留特征的范围比滑动平均要更为接近原始ΨC,而且避免了探测误差零值(264~266)产生的干扰,其误差零值的订正值与原始ΨC平均趋势保持在±2°可信偏差范围内;综合卡尔曼滤波虽然对原始ΨC细节特征保留最好,但是无法完全避免探测误差零值干扰,而且平滑程度较其他滤波方法来说最低(波动指数FIX可表征径向距离廓线波动程度,表 3)。(3)中值滤波虽然对于原始ΨC平均趋势保留最好,但是其对于原始ΨC细节特征保留较差(如图 4a1中距离库数220~224、244~248,图 4a2中距离库数540~560、580~590、660~680)。(4)综合FIR对于原始ΨC径向距离廓线的平滑效果最好(表 3),但是此方法滤波后对于原始ΨC细节特征保留较差,并且造成ΨC形变,失真较为严重,如图 4f1中距离库数216~239、253~258,图 5f1中距离库数270~274上探测误差零值(264~266)导致其超出±2°可信偏差范围。

因此,从保留原始ΨC平均趋势和细节特征以及平滑原始ΨC径向距离廓线三方面综合考虑,综合小波去噪方法相比较其他四种滤波方法是最优秀的,其能较好地保留原始ΨC平均趋势及细节特征(偏差基本在±2°可信范围内,部分区域由于对原始ΨC细节特征的保留会略有超出),避免探测误差零值对滤波结果产生的影响(误差零值订正值回归平均趋势水平,偏差在±2°可信范围内),对原始ΨC距离廓线的平滑程度达到了94%~97%(平滑了超出±1.5° Gate−1的波动,对±1.5° Gate−1内的波动保留)。

利用滤波后的近似ΦDP估算KDP图 6),由于1°~2°内的统计波动被认为是ΦDP估算典型存在的边界,并且结合Hubbert and Bringi(1995)指出的ΦDP理论阈值应在8° km−1~10° km−1(对应本文1.2° Gate−1~1.5° Gate−1),本文将±1.5°看作ΦDP波动边界,因此从理论上讲[公式(2)],估算出的KDP值在±5° km−1内是可信的,所以将±5° km−1作为KDP的边界阈值。

图 6 层积云(Sc)和积雨云(Cb)中KDP:(a1、a2)原始、(b1、b2)综合小波去噪、(c1、c2)综合滑动平均、(d1、d2)中值滤波、(e1、e2)综合卡尔曼滤波和(f1、f2)综合FIR滤波(其中红色虚线表示±5° km-1, 蓝色虚线表示未经过综合处理),时间与图 4一致 Fig. 6 The KDP range profiles for Sc and Cb respectively, (a1, a2) raw, (b1, b2) integrated wavelet de-noising, (c1, c2) integrated running mean, (d1, d2) median, (e1, e2) integrated Kalman, and (f1, f2) integrated FIR filters (the red solid lines denote ±5° km-1, blue dashed lines denote profiles untreated by integrated process), the time is the same as Fig. 4

图 6可知:(1)经过A、B类处理后的综合滤波方法很好地避免了探测误差零值以及δ和环境噪声引起的KDP估算误差;(2)中值滤波、综合卡尔曼滤波由于平滑效果的影响(中值滤波后虽然ΨC波动指数FIX较小,但是其距离廓线并不光滑,毛刺较多)导致估算的KDP抖动明显,其结果不利于降水估测,并且两种方法滤波后估算的KDP存在超出边界阈值5° km−1的现象,不可取;(3)综合滑动平均、综合小波去噪以及综合FIR滤波后估算的KDP都在边界阈值5° km−1内,但是综合FIR滤波方法由于失真较为严重,导致估算的KDP对于原始ΨC细节波动特征保留较差(尤其是对降水大值区的指示效果较差)。综合滑动平均滤波方法在保留原始ΨC细节特征和平滑效果上都不如综合小波去噪。并且综合小波去噪方法在层积云降水以及积雨云降水系统中,对于较强降水区域的指示作用很明显[层积云中较强降水区域(29.7~38.85 km)KDP普遍大于0.25° km−1,峰值可达到4.3° km−1;积雨云中较强降水区域(93.45~102.3 km)KDP普遍大于0.1° km−1,峰值可达到1.95° km−1]。

对于整个层积云降水系统以及积雨云降水系统,根据图 3图 7可知:(1)原始ΨC由于受环境噪声和δ影响严重,估算的KDP不可取,与降水回波不能很好的对应;(2)中值滤波后由于近似估算的ΦDP不够平滑,所以导致估算的KDP质量较差,不可用于降水估测;(3)综合FIR滤波后估算的KDP最均匀,但是由于其大尺度的平滑效果导致对于原始细节特征的保留效果较差,并且在层积云降水系统中东北东方向约24~36 km处存在超出5° km−1边界阈值的现象较其他滤波方法严重,同时边界效应十分明显(Hu and Liu, 2014),在降水系统边缘原始信息失真严重;(4)综合卡尔曼滤波、综合滑动平均和综合小波去噪方法对于原始细节特征的保留较为完整,其中综合卡尔曼滤波和综合滑动平均估算的KDP不如综合小波去噪估算的均匀,并且综合滑动平均估算的KDP也存在超出±5° km−1边界阈值的现象较综合小波去噪严重;(5)整体而言,综合小波去噪方法估算的KDP相对均匀并且较好地保留了降水系统原始细节特征,取值也在±5° km−1边界阈值内,能较好地对应降水系统。

图 7 层积云(Sc)和积雨云(Cb)KDP PPI(雷达平面位置显示器), 原始(a1、a2)、综合小波去噪(b1、b2)、综合滑动平均(c1、c2)、中值滤波(d1、d2)、综合卡尔曼滤波(e1、e2)和综合FIR(f1、f2) Fig. 7 PPI (Plan Position Indicator) of KDP for Sc and Cb respectively, (a1, a2) raw, (b1, b2) integrated wavelet de-noising, (c1, c2) integrated running mean, (d1, d2) median, (e1, e2) integrated Kalman, and (f1, f2) integrated FIR filters after integrated process

综上,综合小波去噪方法对于2015年8月7日的层积云降水资料以及2015年6月22日的积雨云降水资料来说,不仅在滤波效果上优于中值滤波、综合滑动平均、综合卡尔曼滤波及综合FIR滤波,而且滤波后估算的KDP也优于以上四种方法。同时,对于积雨云较强降水区域KDP取值较层积云降水小,并且负值较多,主要考虑是由于该降水系统为强盛对流单体,对流云发展较高,冷云过程较强,云中存在较多固态粒子导致。

4 总结与讨论

本文为了提高X波段双线偏振多普勒雷达的共极化差分相移ΨC质量,构造了对于ΨC的综合滤波方法,在滑动平均、小波去噪和卡尔曼滤波方法滤波前先对原始ΨC进行资料识别和逐步订正的综合处理,在FIR滤波前进行了逐步订正的综合处理,对中值滤波不做处理。利用位于北京顺义(40.18°N,116.68°E)的714XDP-A型车载X波段双线偏振多普勒天气雷达探测资料对综合滤波方法进行了应用与评估,得到以下结论:

(1)综合小波去噪滤波方法不论是滤波处理还是对于KDP的估算相较于其他方法都是最优秀的,其能较好地保留原始ΨC平均趋势及细节特征(偏差基本在±2°可信范围内,部分区域由于对原始ΨC细节特征的保留会略有超出),避免探测误差零值对滤波结果产生的影响(误差零值订正值回归平均趋势水平,偏差在±2°可信范围内),对原始ΨC距离廓线的平滑程度达到了94%~97%(平滑了超出±1.5° Gate−1的波动,对±1.5° Gate−1内的波动保留)。

(2)综合小波去噪方法滤波后估算的KDP整体在±5° km−1边界阈值内,避免了原始ΨC不合理零值、环境噪声和δKDP估算产生的误差影响,KDP对较强降水指示明显(在回波强度30 dBZZH<37 dBZ的层积云降水中KDP普遍大于0.25° km−1,峰值可达到4.3° km−1;在回波强度30 dBZZH<50 dBZ的积雨云降水中KDP普遍大于0.1° km−1,峰值可达到1.95° km−1)。

综上所述,综合小波去噪滤波方法对于X波段双偏振雷达的ΨC资料滤波效果优良,滤除了δ和环境噪声,有效抑制了原始ΨC零值误差对滤波结果以及KDP估算的影响,其滤波后近似ΦDP估算的KDP在±5° km−1边界阈值内,KDP对于较强降水指示明显,但是在回波强度30 dBZZH<50 dBZ的积雨云降水中KDP取值相对于回波强度30 dBZZH<37 dBZ的层积云降水中较小,本文考虑是因为积雨云系统中固态粒子较多导致。综合小波去噪滤波方法对于KDP估算结果的改善将有利于提高降水估测和粒子识别的精准度,对于雷暴天气预警预报具备重要意义。但是由于观测资料有限,综合小波去噪滤波方法对于降水估测及粒子识别的改善作用还需要大量的观测资料以及借助地面观测站数据来验证,这将会是本文后续进行的工作。

参考文献
毕永恒, 刘锦丽, 段树, 等. 2012. X波段双线偏振气象雷达反射率的衰减订正[J]. 大气科学, 36(3): 495-506. Bi Yongheng, Liu Jinli, Duan Shu, et al. 2012. Attenuation correction of reflectivity for X-band dual-polarization radar[J]. Chinese Journal of Atmospheric Sciences (in Chinese), 36(3): 495-506. DOI:10.3878/j.issn.1006-9895.2011.11075
杜牧云, 刘黎平, 胡志群, 等. 2012a. 双线偏振雷达差分传播相移的质量控制[J]. 应用气象学报, 23(6): 710-720. Du Muyun, Liu Liping, Hu Zhiqun, et al. 2012a. Quality control of differential propagation phase shift for dual linear polarization radar[J]. Journal of Applied Meteorological Science (in Chinese), 23(6): 710-720. DOI:10.3969/j.issn.1001-7313.2012.06.008
杜牧云, 刘黎平, 胡志群, 等. 2012b. 双线偏振雷达差分传播相移的小波滤波初探[J]. 暴雨灾害, 31(3): 248-254. Du Muyun, Liu Liping, Hu Zhiqun, et al. 2012b. Preliminary study on wavelet filtering of differential propagation phase shift for dual linear polarization radar[J]. Torrential Rain and Disasters (in Chinese), 31(3): 248-254.
何宇翔, 吕达仁, 肖辉, 等. 2009. X波段双线极化雷达反射率的衰减订正[J]. 大气科学, 33(5): 1027-1037. He Yuxiang, Lü Daren, Xiao Hui, et al. 2009. Attenuation correction of reflectivity for X-band dual polarization radar[J]. Chinese Journal of Atmospheric Sciences (in Chinese), 33(5): 1027-1037. DOI:10.3878/j.issn.1006-9895.2009.05.13
肖艳姣, 王斌, 陈晓辉, 等. 2012. 移动X波段双线偏振多普勒天气雷达差分相位数据质量控制[J]. 高原气象, 31(1): 223-230. Xiao Yanjiao, Wang Bin, Chen Xiaohui, et al. 2012. Differential phase data quality control of mobile X-band dual-polarimetric Doppler weather radar[J]. Plateau Meteorology (in Chinese), 31(1): 223-230.
张培昌, 杜秉玉, 戴铁丕, 等. 2000. 雷达气象学[M]. 北京: 气象出版社: 417-418. Zhang Peichang, Du Bingyu, Dai Tiepi, et al. 2000. Radar Meteorology[M] (in Chinese). Beijing: China Meteorological Press: 417-418.
Bringi V N, Chandrasekar V. 2001. Polarimetric Doppler Weather Radar:Principles and Applications[M]. Cambridge, New York: Cambridge University Press: 636pp.
Bringi V N, Chandrasekar V, Balakrishnan N, et al. 1990. An examination of propagation effects in rainfall on radar measurements at microwave frequencies[J]. J. Atmos. Oceanic Technol., 7(6): 829-840. DOI:10.1175/1520-0426(1990)007<0829:AEOPEI>2.0.CO;2
Chang W Y, Vivekanandan J, Wang T C C. 2014. Estimation of X-band polarimetric radar attenuation and measurement uncertainty using a variational method[J]. Journal of Applied Meteorology and Climatology, 53(4): 1099-1119. DOI:10.1175/JAMC-D-13-0191.1
Cho Y H, Lee G W, Kim K E, et al. 2006. Identification and removal of ground echoes and anomalous propagation using the characteristics of radar echoes[J]. J. Atmos. Oceanic Technol., 23(9): 1206-1222. DOI:10.1175/JTECH1913.1
Doviak R J, Bringi V, Ryzhkov A, et al. 2000. Considerations for polarimetric upgrades to operational WSR-88D radars[J]. J. Atmos. Oceanic Technol., 17(3): 257-278. DOI:10.1175/1520-0426(2000)017<0257:CFPUTO>2.0.CO;2
Gorgucci E, Baldini L. 2015. Influence of beam broadening on the accuracy of radar polarimetric rainfall estimation[J]. Journal of Hydrometeorology, 16(3): 1356-1371. DOI:10.1175/JHM-D-14-0084.1
Gosset M. 2004. Effect of nonuniform beam filling on the propagation of radar signals at X-band frequencies. Part Ⅱ:Examination of differential phase shift[J]. J. Atmos. Oceanic Technol., 21(2): 358-367. DOI:10.1175/1520-0426(2004)021<0358:EONBFO>2.0.CO;2
Hu Z Q, Liu L P. 2014. Applications of wavelet analysis in differential propagation phase shift data de-noising[J]. Adv. Atmos. Sci., 31(4): 825-835. DOI:10.1007/s00376-013-3095-y
Hubbert J, Bringi V N. 1995. An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements[J]. J. Atmos. Oceanic Technol., 12(3): 643-648. DOI:10.1175/1520-0426(1995)012<0643:AIFTFT>2.0.CO;2
Jameson A R, Caylor I J. 1994. A new approach to estimating rainwater content by radar using propagation differential phase shift[J]. J. Atmos. Oceanic Technol., 11(2): 311-322. DOI:10.1175/1520-0426(1994)011<0311:ANATER>2.0.CO;2
Mueller E A. 1984. Calculation procedure for differential propagation phase shift[C]//Preprints 22nd Radar Meteorology Conference. Zurich: AMS, 397-399.
Park H S, Ryzhkov A V, Zrnić D S, et al. 2009. The hydrometeor classification algorithm for the polarimetric WSR-88D:Description and application to an MCS[J]. Wea. Forecasting, 24(3): 730-748. DOI:10.1175/2008WAF2222205.1
Park S G, Bringi V N, Chandrasekar V, et al. 2005. Correction of radar reflectivity and differential reflectivity for rain attenuation at X band. Part Ⅰ:Theoretical and empirical basis[J]. J. Atmos. Oceanic Technol., 22(11): 1621-1632. DOI:10.1175/JTECH1803.1
Park S G, Kim J H, Ko J S, et al. 2015. Identification of range overlaid echoes using polarimetric radar measurements based on a fuzzy logic approach[J]. J. Atmos. Oceanic Technol., 33(1): 61-80. DOI:10.1175/JTECH-D-15-0042.1
Ryzhkov A V. 2007. The impact of beam broadening on the quality of radar polarimetric data[J]. J. Atmos. Oceanic Technol., 24(5): 729-744. DOI:10.1175/JTECH2003.1
Ryzhkov A V, Zrnić DS. 1998. Beamwidth effects on the differential phase measurements of rain[J]. J. Atmos. Oceanic Technol., 15(3): 624-634. DOI:10.1175/1520-0426(1998)015<0624:BEOTDP>2.0.CO;2
Ryzhkov A V, Giangrande S E, Schuur T J. 2003. Rainfall measurements with the polarimetric WSR-88D radar[R]. Norman: National Severe Storms Laboratory, 98pp.
Sachidananda M, Zrnić D S. 1986. Differential propagation phase shift and rainfall rate estimation[J]. Radio Sci., 21(2): 235-247. DOI:10.1029/RS021i002p00235
Seliga T A, Bringi V N. 1976. Potential use of radar differential reflectivity measurements at orthogonal polarizations for measuring precipitation[J]. J. Appl. Meteor., 15(1): 69-76. DOI:10.1175/1520-0450(1976)015<0069:PUORDR>2.0.CO;2
Straka J M, Zrnić D S, Ryzhkov A V. 2000. Bulk hydrometeor classification and quantification using polarimetric radar data:Synthesis of relations[J]. J. Appl. Meteor., 39(8): 1341-1372. DOI:10.1175/1520-0450(2000)039<1341:BHCAQU>2.0.CO;2
Wang Y T, Chandrasekar V. 2009. Algorithm for estimation of the specific differential phase[J]. J. Atmos. Oceanic Technol., 26(12): 2565-2578. DOI:10.1175/2009JTECHA1358.1
Zrnić D S, Ryzhkov A V. 1999. Polarimetry for weather surveillance radars[J]. Bull. Amer. Meteor. Soc., 80(3): 389-406. DOI:10.1175/1520-0477(1999)080<0389:PFWSR>2.0.CO;2