X波段双偏振天气雷达组网反射率因子订正研究

高 霞, 周 晔, 周淑玥, 杜 捷, 王海江

(1.中国航空工业集团公司雷华电子技术研究所,江苏 无锡 214000;2.成都信息工程大学电子工程学院,四川 成都 610225)

0 引言

雷达电磁波在传输过程中,能量会被沿途接触到的各类粒子吸收转变成机械能、内能、热能等[1],电磁波的衰减与波长成反比。X波段双偏振天气雷达由于其波长较短,在探测降雨时能量会被降水粒子吸收等,导致能量衰减,在远距离处的弱雷达回波返回时很难被雷达接收机接收,导致远距离处信息丢失,严重影响雷达的探测精度以及雷达探测资料的应用。

X波段双偏振天气雷达存储有差分传播相移ΦDP,该参数为距离累积量不受衰减影响[2],因此被广泛用于反射率因子的衰减订正。目前较为广泛的衰减订正方法是通过计算衰减率来完成反射率因子的订正,但衰减率大都是给定一个固定的区间范围[3],存在诸多弊端,如不能灵活适应不同降水粒子的衰减。因此可以根据差分传播相移ΦDP不受衰减的影响来设计自适应算法完成衰减订正工作[4-5]。

首先对单部雷达进行反射率因子衰减订正处理,其次将衰减订正处理后的3部X波段双偏振天气雷达反射率因子插值到统一的CAPPI网格上,再对组网反射率因子CAPPI进行缺空补值和斑点滤波处理,最后设计区域检测算法,对组网反射率因子进行订正。

1 资料选取

3部X波双偏振天气雷达隶属于北京市气象局各区分局的雷达站点,分别是北京房山站雷达、北京昌平站雷达、北京顺义站雷达,另外用于验证的1部SA波段雷达为北京大兴站雷达。各站点雷达的详细参数如表1所示。

表1 雷达站点及详细参数

研究数据来源于2020年夏季(8-9月)雷达探测数据(表1),由于夏季降水丰沛,雷达联合组网后效果更佳。

2 单部雷达反射率因子衰减订正

因为衰减是双程的累积量,故衰减订正为

2.1 AH-KDP法

对于衰减率的计算,Bringi等[5]通过散射的模拟实验发现反射率因子的衰减率AH与差分传播相移率KDP基本上呈线性关系,并在实验的基础上提出了衰减订正公式:

其中,α为衰减系数,Testud指出对于X波段雷达,衰减系数α在0.139～0.335[3],本文取其均值0.237。

2.2 自适应算法

为更好地适应不同降雨类型,Bringi等[5]对上节算法进行改进,命名为自适应约束算法。自适应约束算法在计算中使用了差分传播相移ΦDP和水平反射率因子ZH,并根据误差大小不断调整衰减系数α,将其调整为最佳值。自适应算法流程多,迭代次数大,算法涉及大量区间内积分运算,对差分传播相移有较高的质量要求,因此需要对其进行质控,该算法流程如图1所示。

图1 自适应算法流程

自适应算法流程可以概括为:对每个体扫数据的所有仰角层逐层订正,每层进行逐径向订正,将每个待订正的径向分为若干区间,分别计算这些区间内对应的衰减系数α,再根据径向内所有区间的衰减系数对该径向进行订正,然后完成整层数据的订正,最后再完成整个体扫数据的订正。算法具体步骤如下:

(1)给定衰减系数α的初始值为0.01。

(2)利用ZH、降水量I和差分传播相移ΦDP计算衰减率AH。

其中,I为降水量,b为衰减指数常数,取值0.76～0.84,本文取均值0.8,α为区间待求衰减系数,ZH为订正前水平反射率因子(雷达存储值)。

(3)根据KDP和AH的关系、KDP和ΦDP的关系估算

(5)衰减系数α以0.05为步长,递增,重复步骤(4),当两者的差值最小时,即得到最佳衰减系数。

(6)应用最佳衰减系数α进行订正。

2.3 两种衰减订正方法效果对比分析

因S波段天气雷达波长较长,电磁波受雨区衰减影响较小,因此验证数据集为北京大兴站同时刻(两雷达探测误差在1分钟以内)同区域内的S波段单偏振天气雷达数据。通过两种方法对不同地区雷达反射率因子数据进行订正发现,固定衰减系数的AH-KDP法在近雷达处能够对反射率因子有较弱订正,在远离雷达处效果不佳。AH-KDP法的反射率因子订正前后距离廓线对比如图2所示。反射率因子订正前后对比PPI图如图3所示。而自适应算法在近、远离雷达处效果都较明显,能较好地恢复真实降雨区数据。自适应算法的反射率因子订正前后距离廓线对比如图4所示。反射率因子订正前后对比PPI图如图5所示。

图2 AH-KDP法订正值与探测值及验证值对比图

图3 AH-KDP法衰减订正效果对比图

图4 自适应算法订正值与探测值及验证值对比图

图5 自适应算法衰减订正效果对比图

为更直观反映两种算法的优劣,对同种X波段雷达探测数据分别用两种算法进行订正,并以同时刻(误差在1分钟左右)S波段雷达探测数据作验证,结果表明自适应算法明显优于AH-KDP法,对比图见图6。

图6 AH-KDP法与自适应算法效果对比图

根据自适应算法和AH-KDP法的效果对比,本文最终选择自适应算法对组网前的3部X波段双偏振天气雷达的反射率因子数据进行第一步衰减订正处理。

3 网格化处理

3.1 插值处理

常用的线性插值算法包括:垂直线性插值、垂直水平线性插值和八点线性插值算法。八点线性内插(EPI)选取网格点周围8个相关联的距离库作为参考点,将参考维度扩展至仰角、方位角、径向上,能够提升内插的可靠性。该算法解释为选取该网格点对应在极坐标系下最近的两个仰角层和方位角(或径向)上六面体的8个顶点作为参考点[6]。

如图7所示,p1、p2、p3、p4为位于该网格点p上方仰角层内最近的4个距离库参考点,p5、p6、p7、p8为位于该网格点p下方仰角层内最近的4个距离库参考点;p1、p3、p5、p7 位于同一方位角,p2、p4、p6、p8 位于同一方位角;p1、p2、p5、p6拥有相同的径向距离库,p3、p4、p7、p8拥有相同的径向距离库。p1～p8对应坐标为 p1(e1,a1,r1)、p2(e1,a2,r1)、p3(e1,a1,r2)、p4(e1,a2,r2)、p5(e2,a1,r1)、p6(e2,a2,r1)、p7(e2,a1,r2)、p8(e2,a2,r2)。p1～p8对应值为V1～V8,待插值点p在空间笛卡儿坐标系下的坐标为p(x,y,z),需要将其转换至极坐标系下的坐标p(e,a,r),转换公式如式(7),最后用三个维度的线性插值对待插值点p的值Vp进行计算,计算公式如式(8)。

图7 八点线性插值法示意图

其中的we、wr、wa分别为各参考点在仰角、径向距离和方位角上的权重系数,其计算公式如式(9)～(11)所示。

3.2 组网拼图处理

在得到多部雷达网格化处理后的雷达资料后,常见的拼图算有以下几种:最临近法、最大值法、反距离权重函数法[7]。反距离权重函数法相较于最大值法和最临近法,综合考虑了组网所用的多部雷达反射率因子在同一网格点的贡献,能够明显地减少其他两种方法带来的色标跳变现象,组网图像更连续。因此,组网处理部分最终选择反距离权重函数法,其示意图如图8所示。

图8 反距离权重法示意图

由图8可知,3部雷达形成一个三角形区域,网格点p的值由该网格点到3部雷达R1、R2、R3距离决定,最后根据网格点p到3部雷达在三维x-y-z空间上的直线距离r1、r2、r3进行反距离加权,设 r1、r2、r3大小递减。3部雷达的坐标为R1(x1,y1,z1)、R2(x2,y2,z2)、R3(x3,y3,z3),p点的坐标为p(x,y,z)。3部雷达在p点的值分别为V1、V2、V3,p点拼图后的值应为Vp。拼图公式如下:

3.3 网格化处理结果

对2020年8月12日10:00:00北京房山站、北京昌平站、北京顺义站X波段双偏振天气雷达的水平反射率因子PPI以八点线性插值方法进行3 km高的CAPPI插值与组网处理,其插值结果、组网拼图结果分别如图9所示。

由图9可知,EPI法在3个维度(仰角、方位角、距离)上选择8个参考点,能够对3部X波段双偏振天气雷达极坐标系下的反射率因子资料合适地插值到统一的空间笛卡尔坐标系下(CAPPI)。由图10可知,反距离权重函数法基本能将3部雷达的插值结果较完整的组网到同一个网格上。至此,完成前期的3部雷达的反射率因子组网处理。

图9 三地区插值处理CAPPI图

图10 组网拼图处理CAPPI图

4 组网反射率因子订正

4.1 预处理

预处理部分主要分为缺空补值和斑点滤波[8-9],关于反射率因子CAPPI的缺空补值处理,设计NA×NR窗口滤波函数[10]对所有网格点上的值进行滤波检测。该窗口滤波函数可描述为:在每格网格点的第i个经度和第j个纬度处的有效值,其经度方向上的窗口尺度为NA,其纬度方向上的窗口尺度为NR,当窗口内的有效数据总数Ti,j在窗口内总库数占比Pi,j小于阈值Thr1时,则将该库视为孤立回波[11],并将该库处的值设置为无效值NaN,如公式(15)。当窗口内的有效数据总数Ti,j在窗口内总库数占比Pi,j大于阈值Thr2时,则将该库视为缺测值库,并将该库处的值从NaN设置为一个有效值,该有效值F为该窗口内所有有效值的均值,如式(16)。

其中阈值Thr1设置为0.45,阈值Thr2设置为0.3,对北京房山、昌平、顺义地区经过EPI法得到的反射率因子CAPPI进行缺空补值、斑点滤波处理前后的对比如图11所示。

图11 预处理前后CAPPI图

由图11预处理前后对比,可知预处理能够对边缘空缺值进行有效回填,并对组网区域内的斑点杂波进行相应的去除。

4.2 区域检测算法

关于组网反射率因子CAPPI的订正处理,由于在应用反射率因子进行CAPPI组网拼图前,已经对3部雷达的反射率因子以自适应算法进行订正处理,但这种处理不够彻底。因为订正是基于雷达电磁波在发射方向上不断衰减能量的前提下,运用最佳衰减系数和差分传播相移使探测值尽量靠近真实值。因此在反射率因子CAPPI组网之后,也要对网格点上的反射率因子值进行适当订正。本文设计使用区域检测算法对CAPPI每个网格点进行订正。该算法整体思路为以S波段雷达为参考,对拼图CAPPI和S波段雷达CAPPI的每个网格点建立维度为N的搜索框,并分别计算其均值Mean_S和Mean_X,当前者与后者的差值大于阈值Thr_C时,则判定为偏差严重,该偏差包括系统误差和衰减,需要进行订正。为防止订正后的值大于真实值,以尽量靠近真实值的为原则,订正后的值Vp为当前拼图CAPPI值VX加上Mean_S与Mean_X的差值,其中N取5,阈值Thr_C取2 dBZ。建立的搜索框示意图如图12。

图12 搜索框示意图

订正公式如式(17),组网反射率因子CAPPI订正前后对比效果如图13。

图13 订正处理前后对比图

由图13可知,订正部分主要在39°N～40.3°N、115.5°E～118°E附近,将订正后的组网反射率因子CAPPI与同时刻同地点的S波段雷达反射率因子CAPPI进行对比,如图14。

图14 订正CAPPI与验证CAPPI对比图

4.3 订正结果分析

由图 13 可 知,115.5°E～118°E,39°N～40.3°N为主要订正区域,区域检测算法以待订区域内同时刻下S波段天气雷达探测资料为参考,对每个待订正的网格点反射率因子建立搜索框,并计算两个搜索框内的区域均值,当两均值差异较大时,则判定为偏差过大,利用S波段雷达资料对其进行订正处理。由图14可知,主要订正区域内的订正结果基本与验证用S波段雷达反射率因子CAPPI相吻合,该算法订正效果表现良好。

5 结束语

在前人研究的基础上,首先通过分析2种常规衰减订正算法的优劣,选择较优算法对单部雷达反射率因子订正处理;其次运用经典八点线性插值算法和反距离权重函数法对订正后的3部X波段双偏振天气雷达反射率因子进行网格化处理,最后设计了区域检测算法对组网反射率因子CAPPI进行适当的订正处理,实验结果表明,本文设计的区域检测算法在CAPPI网格上表现良好的效果,具有一定的参考价值,但该方法局限于组网区域内需要存在衰减较小的长波长天气雷达用以验证算法效果,极大限制了该方法在业务上的应用,未来需设计其他数值订正方案来优化该算法。