移动X波段双线偏振雷达数据质量分析及偏差订正*

2021-04-20冯婉悦刘新安

气象 2021年3期

张林李峰冯婉悦刘新安

1 中国气象局气象探测中心，北京 100081 2 新疆维吾尔自治区气象技术装备保障中心，乌鲁木齐 830002 3 南京恩瑞特实业有限公司，南京 211106

提要： 2019年8月，中国电子科技集团第十四研究所南京恩瑞特实业有限公司研制的GLC-12A型移动X波段双线偏振多普勒天气雷达在江苏大丰进行外场观测。根据该雷达外场观测资料，联合盐城新一代SA天气雷达和江苏省6个雨量站资料，对移动X波段双线偏振雷达做数据质量分析。结果表明，针对同一次降水过程而言，距离越远，衰减越大。为了降低衰减造成的观测误差，首先利用信噪比对差分反射率Zdr、零滞后相关系数ρhv(0)进行偏差订正，提高Zdr和ρhv(0)的数据质量，以便进行地杂波、生物回波等非降水回波过滤。然后应用去抖动滤波后的差分传播相移率Kdp或总差分传播相位φdp表达在雨中距离雷达R处的总双向反射率因子衰减订正值。数据分析表明，订正后移动X波段双线偏振雷达数据质量得到了提高。

引言

双线偏振雷达是指可以发射水平和垂直偏振状态电磁波的雷达系统。它可以通过探测到的原始数据反演出一些新的物理量，如差分反射率因子Zdr、差分传播相位φdp、差分传播相移率Kdp、零滞后相关系数ρhv(0)等偏振参数。使用这些新物理量，可以提高测雨精度和识别冰雹的准确率，了解一定区域内降水粒子相态与形状，确定雨滴谱参数等，以及对非球形的云雨等降水粒子群的探测也都具有实际意义(张培昌等，2018)。

目前，美国已全面完成158部布网S波段雷达(WSR-88D)双线偏振雷达升级和业务化应用。我国也正在加快推进双线偏振雷达业务化发展进程，广东全省雷达和江苏省部分雷达已率先完成了S波段WSR-98D双线偏振雷达升级和业务化(张林等，2018)。

关于雷达数据质量分析，王超等(2019)分析了北京市气象局和佛山市气象局的X波段双线偏振雷达资料，统计偏振测量参数的变化规律，从而分析了地物、避雷针和旋转及俯仰关节等对双偏振参量测量精度的问题，并提出对应的质量控制方法。杜牧云等(2013)分析了升级改造前后的C波段双线偏振雷达资料，指出信噪比小于25 dB时，Zdr和ρhv(0) 受噪声影响存在明显测量误差。张志强和刘黎平(2011)、吴翀等(2014)分析提出了基于不同地理位置、不同分辨率的雷达反射率因子匹配方法和观测资料的定量对比方法。张林等(2018)通过时间、空间一致性匹配，将相邻两部雷达重叠区资料进行定量对比分析。本文用新一代SA天气雷达数据和自动雨量站资料与X波段双线偏振雷达定量对比分析，结果表明X波段双线偏振雷达信号衰减为影响其观测误差的重要因素。

针对X波段双线偏振雷达信号衰减问题，国内学者也开展了大量的研究工作，胡志群等(2012)采用Zh-Kdp综合订正法进行衰减订正并评估其对降水估测的影响。王晗等(2016)对比分析了Kdp订正法、Zh订正法和ZPHI降水廓线订正法在X波段双线偏振雷达衰减订正中的使用效果，结果表明，Kdp订正法和ZPHI订正法的订正效果比较好，优于Zh订正法。毕永恒等(2012)采用自适应约束算法对X波段雷达进行衰减订正，分析表明当存在较大范围降水时，订正后的效果显著优于订正前，能提高降水估测精度。

关于双线偏振雷达数据质量控制，美国强风暴实验室(NSSL)研发了美国业务双线偏振雷达(WSR-88D)非降水回波识别算法，并应用于定量估测降水业务中(Tang et al,2014)。张林等(2018)针对中国上海南汇S波段WSR-88D双线偏振雷达，设计了数据质量控制流程。朱轶明等(2019)研究了上海WSR-88D双线偏振生物回波识别算法。肖艳姣等(2012)针对移动X波段双线偏振雷达数据质量问题，提出了一套差分传播相位数据质量控制算法，用来识别和消除差分传播相位数据中的地物杂波、退折叠的相位和滤除正常的随机差分相位波动。杜牧云等(2012)处理C波段双线偏振雷达资料，根据信噪比和相关系数将差分传播相位资料分为较好、较差和差这三类分别处理，提高差分传播相位数据质量。差分传播相移率Kdp是雷达波束穿过雨区时，单位距离内水平与垂直偏振方向相位差引起的，与雨区衰减无关，普遍用于短波长双线偏振雷达型号衰减订正上。而Kdp订正法的前提是数据质量稳定可靠。原始资料中可能存在地杂波等非降水回波，虽然ρhv(0)和Zdr是消除非降水回波的有效偏振参量(Tang et al,2014；张林和杨洪平，2018)，但两者在信噪比较小时容易产生测量偏差，因此需要在衰减订正前对偏振测量参数进行偏差订正。

2008年12月国家重点研发计划“龙卷风探测雷达研制及业务化应用研究”项目落地，项目计划在龙卷风灾害风易发区江苏省大丰、高邮等多个地区布设多波段多体制雷达(包括移动X波段双线偏振雷达、C波段相控阵雷达等)，开展外场观测试验，捕捉龙卷、强对流天气，利用高分辨率的观测资料分析龙卷精细化结构及演变特征。本文分析和处理2019年8月在龙卷风观测区域开展试验的GLC-12A型移动X波段双线偏振雷达观测资料，利用Kdp表达在雨中任意距离的水平反射率因子衰减订正值。为龙卷风观测区域多部雷达组网拼图、研究分析龙卷精细结构和演变特征奠定基础。

1 GLC-12A型移动X波段双线偏振雷达系统和外场试验

1.1 性能参数

表1列出了GLC-12A型移动X波段双线偏振雷达主要性能指标，在本次观测过程中，采用的是单发双收的极化方式，距离库长为60 m，距离范围为75 km。

表1 GLC-12A型移动X波段双线偏振雷达主要性能参数Table 1 The main parameters of GLC-12A mobile X band dual polarization weather radar

1.2 试验场地和位置

国家重点研发计划“龙卷风探测雷达研制及业务化应用研究”第二课题“龙卷探测雷达组网技术及自适应协同观测试验”负责外场观测试验。为尽快取得龙卷、强对流观测资料，由江苏省气象局联合中国电子科技集团第十四研究所南京恩瑞特实业有限公司，于2019年8月在江苏省大丰市进行外场观测试验。图1a为GLC-12A型移动X波段双线偏振雷达样机。雷达架设在一栋楼顶上(图1b)。图1c为高德地图上显示该X波段双偏振雷达、盐城新一代SA天气雷达和江苏省7个自动雨量站(大丰、盐城、东台、射阳、建湖、兴化、阜宁)位置示意图。

图1 GLC-12A型移动X波段双线偏振雷达外场试验(a)雷达样机,(b)架设位置，(c)位置关系(黄点为雷达位置，红点为雨量站位置)Fig.1 Field Test of GLC-12A mobile X-band dual-polarization weather radar (a) radar hardware, (b) positon, (c) positional relation in Gaode Map(yellow dot: postion of radar, red dot: automatic rainfall station)

根据经纬度与雷达极坐标之间的转换公式(张林和杨洪平，2018；吴翀等，2014)，将雨量站、雷达站经纬度代入公式，计算得到每部雨量站相对于X波段双线偏振雷达、盐城S波段新一代天气雷达的方位和距离，如表2所示。

表2 雨量站与雷达距离Table 2 The distance between automatic rainfall station and radar

2 数据质量分析

2.1 降水过程的雨量站数据统计

阜宁雨量站距移动X波段双线偏振雷达为77.4 km，超出了该雷达在本次观测过程中的最大探测距离。本文采用6个雨量站(大丰、盐城、东台、射阳、建湖、兴化)和盐城新一代SA天气雷达数据作为对比，分析移动X波段双线偏振雷达数据质量。

2019年8月10—11日，在江苏省发生了大面积连续性降水过程。表3给出了6个雨量站在本次过程不同降水时段的最大小时降水量。

表3 2019年8月11日强降水观测过程Table 3 The heavy rainfall observed on 11 August 2019

2.2 雷达回波定性分析

选取2019年8月11日降水过程中03时观测数据，图2a为盐城新一代SA天气雷达观测回波(一圈代表50 km，最大探测距离230 km)，其中红色三角形和圆圈为移动X波段双线偏振雷达位置和观测范围，距盐城新一代SA天气雷达约25 km。图2b为移动X波段双线偏振雷达观测回波(第一圈代表50 km，第二圈代表90 km，回波最远探测距离是75 km)。从图中直观对比可知，移动X波段双线偏振雷达观测的回波结构和轮廓与新一代SA天气雷达基本一致，所不同的是移动X波段双线偏振雷达的距离分辨率(60 m)更高，有利于观测强对流、冰雹及龙卷等灾害天气的精细结构演变，但新一代SA天气雷达探测范围广，灵敏度更高，因此其弱回波的探测能力更强。

图2 2019年8月11日03时雷达观测回波(a)盐城新一代SA天气雷达0.5°仰角观测，(b)移动X波段双线偏振雷达1.5°仰角观测(红色三角和圆圈：X波段雷达架设位置和观测范围)Fig.2 The radar image at 03:00 BT 11 August 2019(a) CINRAD/SA of 0.5° elevation angle in Yancheng, (b) mobile X-band dual-polarization weather radar of 1.5° elevation angle(red triangle and circle: position and observation range of X-band dual-polarization weather radar respectively)

2.3 雷达和雨量站观测值定量分析

根据表2中雨量站在两部雷达上的极坐标(方位、距离)，取雨量站点在两部雷达坐标上的反射率因子观测数据，挑选与雷达观测时间最近时次的5 min 雨量站数据做对比(两部雷达均为6 min观测一次)。图3为6个雨量站数据随时间序列统计，其中横轴为时间序列，主纵轴(左)为反射率因子观测值，次纵轴(右)为5 min雨量值。

图3 两部雷达和6个雨量站观测值定量分析(a)大丰，(b)盐城，(c)东台，(d)射阳，(e)建湖，(f)兴化Fig.3 Quantitative analysis of radar and rainfall data (a) Dafeng, (b) Yancheng, (c) Dongtai, (d) Sheyang, (e) Jianhu, (f) Xinghua

求两部雷达在雨量站处平均反射率因子(dBzSA-dBzX)，如表4所示，平均偏差随距离增大也呈上升趋势。由此可知，X波段双线偏振雷达穿过雨区后受到的衰减与距离有关。针对同一次降水过程而言，距离越远，衰减越大。其中，盐城雨量站处(距离X波段雷达25.5 km)的平均反射率因子偏差较大，可能是因为X波段雷达观测的反射率因子高值区(超过40 dBz)的衰减严重或两者海拔高度的不一致等因素造成统计上存在较大的偏差。

表4 两部雷达在雨量站位置处平均反射率因子偏差Table 5 The average reflectivity factor deviation observed by two radars at rainfall station

3 数据偏差订正

雷达接收的回波功率Pr=Pr*+Nr，其中Pr*是气象目标回波功率，Nr是噪声功率，Pr*应等于Pr-Prmix，取dB后为10lgPr/Prmix=SNR，即信噪比(SNR)代表气象目标回波功率Pr*的dB值。由于接收机的噪声电平受温度影响会随时间变化，故某一时刻测得的回波功率Pr中所包含的噪声功率，不一定正好等于设定的最小可测功率Prmix，则由SNR决定的气象目标回波功率Pr*中，可能仍含有少量噪声或被过多地扣除了噪声，使SNR不能真正代表气象目标回波功率Pr*的dB值，即存在误差，由它反演出的差分反射率Zdr和零滞后相关系数ρhv(0)也会有误差。Zdr和ρhv(0)在SNR较小时的误差变大，这种现象会因发射功率随距离的增大而下降变得更加明显(张培昌等，2018)。而Zdr和ρhv(0)均为区分降水和非降水回波的有效参量，因此需要对Zdr和ρhv(0)数据进行噪声偏差订正，为后续资料质量控制、降水相态分类提供稳定可靠的数据。

3.1 信噪比计算

本文分析的GLC-12A型移动X波段双线偏振雷达没有直接给出SNR值，故SNR通过雷达反射率因子dBzh和探测距离R间接计算得到(张培昌等，2018)：

SNR=dBzh-20lg(R)+C

(1)

式中：SNR单位为dB，R单位为km，C是常数。

如图4所示，2019年8月11日03时为观测时刻，经计算得到1.5°仰角的反射率因子(图4a)和信噪比SNR测量值(图4b)。

3.2 零滞后相关系数ρhv(0)订正

零滞后相关系数ρhv(0)是指当水平和垂直偏振波同时发射与接收时，接收到回波水平偏振分量与垂直偏振分量之间的相关程度。

当信噪比大时，属于气象信息的相关部分占的比例大，故两种偏振的气象回波信号相关性好，相关系数就高；信噪比小时，气象信息相关部分占的比例小，相关系数就变小。从相关函数应反映气象目标特性而言，这个变小就是一种误差，必须通过信噪比的比值加以校正(张培昌等，2018)。ρhv(0)的订正公式如下所示：

ρcor(0)=ρhv(0)×(1+1/snr)

(2)

式中:snr=100.1SNR。图5为2019年8月11日03时1.5°仰角信噪比订正前后的相关系数。从图4a中可知，该时刻观测回波均为降水回波。图5a中50 km之外，受信噪比下降的影响，降水回波的相关系数降低到0.7左右，经过订正后，恢复到0.9以上。

图4 2019年8月11日03时1.5°仰角反射率因子(a)和信噪比计算值(b)Fig.4 Values of reflectivity (a) and SNR (b) at 1.5° elevation angle at 03:00 BT 11 August 2019

图5 2019年8月11日03时1.5°仰角信噪比订正前后的相关系数(a)订正前，(b)订正后Fig.5 The SNR correction of correlation coefficient of 1.5° elevation angle at 03:00 BT 11 August 2019(a) before correction, (b) after correction

根据S波段WSR-88D双线偏振雷达数据质量控制流程(张林和杨洪平，2018)中的非降水回波识别方法，将经过订正后的Zdr和ρhv(0)数据作为输入，识别和消除地物、生物等非降水回波。

3.3 差分传播相位φdp订正

差分传播相位：

φdp=φhh-φvv=δ+φdp

(3)

式中：φhh和φvv分别表示水平及垂直偏振发射波通过相同长度的一个降水区后，散射回天线处的相位值；δ为由于散射造成的双程差分后向散射相移；φdp为前向传播相位差，这里前向是指辐射源到接收体的方向，故粒子散射波射向天线也属于这种情况，它也是由于水平及垂直偏振波经过非球形粒子组成的降水区时造成的。

差分传播相位φdp是传播过程中非球形降水粒子对水平及垂直偏振波的传播相位变化不同而引起的，故φdp值的大小既与粒子形状、相态、取向有关，也与通过降水区的长度有关。φdp的大小能反映降水粒子的状况，其值随雨区增大而变大，当通过百千米以上的非球形雨区时，φdp值可达上百度(张培昌等，2018)。

φdp的准确与否，直接影响差分传播相移率Kdp的估计，在反射率因子Zh和差分反射率因子Zdr的衰减订正中，都要使用准确的φdp，因此需要对φdp数据进行质量控制，包括去折叠、抖动和滤波处理。

相位折叠是指当遇到大范围强降水时，降水回波远端的φdp可能会大于最大可测值而引起的折叠。肖艳姣等(2012)利用模糊逻辑区分地杂波和降水后，然后再看是否要对φdp值退折叠。本文选择将反射率因子和信噪比订正后的相关系数及差分反射率因子作为输入，进入双线偏振雷达数据质量控制流程消除地物、生物回波等杂波(张林和杨洪平，2018)，然后根据径向连续性检查对φdp退折叠处理(张培昌等，2018)。

设径向上相邻两个非杂波点的有效差分传播相位为φdp(i-1)、φdp(i)，若φdp(i-1)-φdp(i)≥140°，那么φdp(i)就被认为是折叠的，对其值加上180°(最大、最小可测量值之间的间隔)就完成了退折叠。

因为后向散射相移以及噪声、地物及旁瓣等因素，引起φdp的值显著抖动。对φdp滤波的基本思路是，对偏离φdp平均趋势较大的值，本文通过比较每个φdp点与其邻近方位和距离库的窗口(可选5×5或9×9)上的φdp值，偏差较大者进行滤波处理。对显著抖动的φdp值进行滤波方法是，首先沿径向对退折叠后的φdp进行4 km窗口的滑动平均，然后比较退折叠后的每个库中的φdp和其滑动平均值，如果二者的绝对差大于给定的阈值(本文取10°)，那么就用滑动平均值取代该距离库的φdp值。2019年8月11日03时观测时刻，0.5°仰角上φdp去抖动和滤波处理前后图像(图6)，红色圈内显示，滤波后的图像剔除了一些φdp抖动点。图7为该时刻239°方位上的φdp随距离的变化曲线，经过滤波后，该方位上消除了一些折叠或抖动的异常值。

图6 2019年8月11日03时1.5°仰角差分传播相位φdp去抖动和滤波处理(红色圈内存在一些抖动点)(a)滤波前，(b)滤波后Fig.6 The de-jitter and filter processing of differential propagation phase at 1.5° elevation angle at 03:00 BT 11 August 2019(a) before filtering, (b) after filtering(The red cicle indicates the presence of some wobble points)

图7 2019年8月11日03时1.5°仰角差分传播相位φdp(239°方位)随距离变化曲线(a)滤波前，(b)滤波后Fig.7 Variation curve of differential propagation phase (239° azimuth) at 1.5° elevation angle with distance at 03:00 BT 11 August 2019(a) before filtering, (b) after filtering

3.4 衰减订正

Bringi and Chandrasekar(2010)通过散射的数值模拟表明：衰减率Ah和差分衰减率Adp=Ah-Av(Ah、Av分别为水平和垂直偏振波在降水区中的衰减率，单位：dB·km-1)与差分传播相移率Kdp之间基本上为线性关系。因此，可利用Kdp或总差分传播相位φdp来表达在雨中距离雷达R处的总双向水平反射率因子衰减订正值ΔZh(R)和差分反射率因子的差分衰减订正ΔZdr(R)，其表达式如下：

ΔZh(R)=2α1·ΔR∑Kdp(i)=α1[φ(R)-φ(0)]

(4)

式中：i=0,1,2，…，N；∑为求i所有对应的Kdp之和。

ΔZdr(R)=2α2·ΔR∑Kdp(i)=α2[φ(R) -φ(0)]

(5)

式中：i=0,1,2，…，N；∑为求i所有对应的Kdp之和。

式(4)和式(5)中：N为雷达到距离R处的总探测库数；ΔR为库长(单位：km)，系数α1和α2依赖于雨滴椭球率。对于X波段双线偏振雷达而言，可取α1=0.25 dB·(°)-1，α2=0.033 dB·(°)-1(Matrosov et al，2002；张培昌等，2018)。φ(0)为沿径向的初始差分传播相位，由于离雷达2 km内差分传播相位受地物杂波和旁瓣等因素影响造成观测值不稳定，因此通过沿径向寻找2 km以外的一段距离，要求在该段距离上自第一个1 km开始连续向下，其φdp(经过去抖动和滤波处理后)在该段上的相关系数(经过订正后)都大于0.9，而且φdp在该段上的标准偏差小于10°(即取φdp波动小的一段)，当选到这样一段距离后，就取该段上的所有φdp的平均值，作为初始相位φ(0)。

采用式(4)和式(5)对反射率因子和差分反射率因子(经过信噪比订正后)进行衰减订正(图8和图9)。图8为反射率因子衰减订正效果，1.5°仰角反射率因子观测数据在方位300°附近存在缺测值，经过衰减订正后，采用2.5°仰角订正后的值补充1.5°仰角层上的缺测值。

图8 2019年8月11日03时衰减订正前后的1.5°仰角反射率因子(红色圆圈表示订正前后效果明显)(a)订正前，(b)订正后Fig.8 Reflectivity factor of before and after attenuation correction at 1.5° elevation angle at 03:00 BT 11 August 2019(a) before correction, (b) after correction(The red cicles indicate that the improvement effect is obvious before and after the revision)

图9 同图8，但为差分反射率因子Fig.9 Same as Fig.8, but for differential reflectivity factor

3.5 两部雷达观测定量对比

关于如何将两部不同地理位置、不同分辨率的雷达反射率因子进行匹配后定量对比，可采用雷达极坐标与经纬度坐标的转换公式(张林和杨洪平，2018；吴翀等，2014)，首先将新一代SA天气雷达极坐标转为经纬度坐标，再转换为X波段双线偏振雷达的极坐标，实现新一代SA天气雷达极坐标数据转为X波段雷达的极坐标数据，然后直接定量对比。

图10是0～30和30～60 km范围的两部雷达衰减订正前后观测反射率因子数据的定量对比散点图，横轴为盐城新一代SA天气雷达观测值，纵轴为X波段双线偏振雷达观测值。如图所示，在弱回波区间，X波段双线偏振雷达观测值略大于S波段新一代天气雷达观测值，随着回波值的增大，S波段新一代天气雷达观测值大于X波段双线偏振雷达观测值的比例逐渐变大，这种现象在30～60 km 区间比在0～30 km区间更明显。

图10 两部雷达衰减0～30 km(a，b)、30～60 km(c，d)范围订正前(a，c)和订正后(b，d)观测反射率因子定量对比Fig.10 Quantitative comparison of reflectivity factor between two radars before attenuation correction (a, c) and after attenuation correction (b, d) in the ranges of 0-30 km (a, b) and 30-60 km (c, d)

经过衰减订正后，由衰减造成的观测误差在整体上得到了改善。表5为统计S波段新一代天气雷达观测值大于X波段双线偏振雷达观测值所占的比例，0～30和30～60 km范围内在订正前的比例分别为48.18%和71.78%。经过衰减订正后，该比例分别降低为42.56%和63.41%，30～60 km范围内的改善更为明显。

表5 S波段新一代天气雷达观测值大于X波段双线偏振雷达观测值所占比例(单位：%)Table 5 Percentages of observation value by S-band weather radar larger than that by X-band dual-polarization weather radar (unit: %)