多普勒天气雷达利用云雨目标物对电磁波的后向散射进行探测, 并测定其空间位置、强弱分布, 从而了解降水的生消、演变和移向移速等, 但其同样能探测到非降水目标物。当非气象目标物被处理进雷达基数据时就造成了杂波污染。由于全部产品和算法是建立在基数据基础上, 这就影响了导出产品, 特别敏感的是降水估计产品。Moszkowicz等^[1]、Joss等^[2]分别统计了2个月和60个月的WSR-88D降水产品, 仅由异常地物杂波造成的降水估计误差分别达15%和13%。因此杂波的去除问题也一直是雷达气象专家研究的重点之一。在新一代天气雷达中, 用RDASOT产生的旁路图来抑制正常的地物回波 (山、建筑物等) , 这种方法对雷达波束正常传播 (NP) 情况下的地物处理效果比较好, 而主要由雷达波束异常传播 (AP) 造成的瞬变杂波过滤则是通过定义一个杂波抑制区来实施, 对该区域选择抑制程度均匀作用在所有径向速度为零附近的杂波, 这种滤波方法显然会降低该区域内的反射率因子值, 从而对产品造成较大影响, 因此全国SA雷达基本上都没有采用这种抑制方法。而这些AP杂波的存在常常会影响雷达资料可信度, 为此, 人们研究了多种地物回波的识别方法, 国外在这方面展开的研究很多, 如在WSR-88D雷达基础上, 应用混沌理论对雷达回波进行分类及产生系列算法等进行滤波^[3-6], 应用雷达三维反射率因子特征进行杂波滤除^[7-8]及应用神经网络进行自动分类和识别^[9]等等。另外, 欧洲的一些常规天气雷达还结合卫星、地面降雨资料及数值预报产品等进行地物回波的识别和消除^[10-13]。国内也有专家在这方面作积极探索, 如刘黎平等开展了基于双线偏振雷达的地物回波识别方法的研究^[14]。所有这些方法都是应用非气象回波在雷达基数据产品上所表现的特征, 通过设计合理的算法, 从软件上再进行实时的补充处理过滤。本文统计分析了宁波2003—2006年多普勒雷达地物杂波的特点, 并结合当前国内外滤波的部分方法, 模拟人眼识别杂波过程中的模糊思维逻辑, 提出了本地化的人工智能的滤波方案, 并在实验中取得了较好的效果。

宁波雷达建在达蓬山顶, 海拔高度458 m, 四周没有遮挡角, 2003年正式投入业务运行。自运行以来, 在强对流的预警和监测中发挥了重要作用。但监测到的固定或非固定的强杂波也对预报分析造成了一定的困扰。统计2003 —2006年影响宁波雷达的非气象杂波, 概括起来主要有两类: ①地物杂波, 包括普通地物杂波 (NP杂波) 和异常传播导致的AP杂波; ②移动的物体 (船、飞机、鸟群、昆虫群等) 。最常见的对导出产品等影响较大的主要是地物杂波。

在宁波的西南山区50 km左右的小块区域常年出现相对固定杂波, 在0.5°仰角反射率因子产品上出现的等效回波强度有时达60 dBz左右 (以下对地物杂波的等效回波强度都用回波强度dBz表示) 。在南侧50 km左右同样还存在一窄带状的沿着山脉走向的东北—西南向的杂波带, 强度稍弱。这些回波总的特征是: ①通常出现在3 km高度以下, 在0.5°仰角的反射率因子产品图上回波强度强, 到1.5°仰角上则迅速减弱, 范围缩小, 即上下仰角面对应的回波差值较大; ②总体呈离散状分布, 强弱交杂, 虽然位置相对固定, 但前后时次还是略有位移, 强弱也有些许变化; ③在速度图上大都为接近零的离散状分布, 同时中间夹杂着较大风速的回波; ④在谱宽产品上0.5°仰角的特征不明显, 虽然大部分对应数值较小, 但中间有时也会出现大于10 m/s的谱宽值, 而到了1.5°仰角谱宽数据则基本上都是小于1 m/s的数据。另外在雷达站周围经常有强度低于30 dBz的大片点状回波, 这些回波在2.4°仰角下都有, 在0.5°仰角的速度图和谱宽产品上有时为不连续分布, 而当空中湿度较大或降雨前后, 速度图则较为连续。这类回波的出现很可能与雷达站周围大气中足够大的灰尘、液滴或昆虫等微粒有关, 在一定程度上能提供雷达站附近的风向, 因此这类杂波在文中没有考虑滤除。

在宁波雷达站的西北面或北面80~200 km距离处, 在低层的反射率因子产品中还经常探测到超强的非降水回波, 等效回波强度甚至达到94 dBz。这些杂波通常出现在晴空的半夜到清晨, 或在降雨前后。其产生与大气折射很有关系, 通常出现在逆温或湿度随高度增加而迅速减少的气象条件下, 这就是通常所说的AP杂波^[15]。但在雷达站的南侧基本观测不到这种杂波, 从地形上来分析, 在宁波的西南面有天台山和四明山绵延起伏于境内, 而在北面, 雷达探测首先要经过杭州湾水面, 接着是杭嘉湖平原地区。产生AP的直接原因是雷达波束在传播过程中由于折射指数的不同导致波束弯曲, 从而观测到地面的目标物。而温度、湿度是造成大气密度变化乃至大气折射率变化的因素。在湖面、海面等大的水域, 大气中水汽成分较山区多, 对温度尤其敏感, 更容易形成上下不同高度大气密度不同。这也就解释了为什么AP杂波通常出现在北面。当宁波雷达监测到这种杂波时, 往往上海和舟山的雷达也能观测到, 只是强度不同。AP杂波的特点是: ①强度变化大, 从弱到强呈离散的无序分布, 每次出现形态范围大小不一, 其爆发和消失经常很迅速, 有时在同一天会反复出现; ②通常出现在0.5°仰角上, 其等效反射率因子值通常很高, 最高可达90 dBz以上, 在1.5°以上仰角则迅速减弱, 上下层的回波差值有时可达50 dB以上; ③平均径向速度图上的典型特征是仰角面上对应一片接近零的速度场, 在0.5°仰角面上中间有时也会孤立地镶嵌着一些非零值; ④0.5°仰角面谱宽产品中相对较乱, 大小谱宽值交杂, 而1.5°仰角以上其对应区域内的谱宽通常数值较低, 数值波动幅度很小。

非气象回波与降水回波最大的区别就是Z, V, W产品中表现为颜色的相对不连续性, 以及分布的相对离散性。而人眼对颜色的敏感性, 对经常出现的杂波形态及相对位置的模糊记忆能力及人类思维与生俱来模糊处理能力等, 使有经验的预报员能根据这些特点区分是降水回波还是杂波污染。因此, 在程序设计前, 首先要对杂波常出现的区域进行深入的统计分析, 并给出可能的同一个仰角面的水平回波梯度, 上下层仰角面的垂直梯度及认为是突变区域时的周边平均梯度等的宽泛阈值。有些杂波如一些弱的AP杂波, 当其混合在降水中时, 由于其本身也具有一定的连续性, 人眼也不能马上作出判断, 这时就需要预报员对这些可疑区域的回波作进一步的思考, 分析其是否是杂波。这就需要结合更高仰角的产品或其他的产品作辅助分析。如上层的回波在形态上是否与当前回波存在相似性, 回波强度是否迅速减小等等。从程序设计的角度来看, 就可以转化为统计可疑区域内的大部分数据是否在上层都存在着对应数据, 计算上下层对应的垂直变化等等。当上层也存在弱降水回波时, 这两个条件可能还不足以识别, 这时预报员还会结合径向速度和谱宽产品进行分析。这在程序设计上就可以转化为检测其风速和谱宽的连续性和数值大小等等。

程序具体设计时, 在充分考虑计算效率与计算效果之间的博弈关系的基础上, 需尽可能遵守以下的原则理念: ①保守性, 设置严格的条件保护有效数据, 防止真正的降水回波被误滤; ②独立性, 对杂波的每种特征进行相对独立的滤波处理; ③多元性, 在杂波的识别工作中尽可能多地利用多元数据层 (例如: Z, V, W) 和空间层 (不同仰角数据) 及其他的一些实况数据。

根据以上的设计理念, 参考Kessing等提出的模糊逻辑方法, 从回波强度、径向速度和速度谱宽等雷达资料中提取用于区分地物杂波的一些物理量, 根据2003—2006年的统计给出宽泛的阈值, 并设计识别和滤除程序。这些物理量包括:反映同一仰角面上回波强度的水平变化Z_H, 上下仰角面回波强度的垂直变化Z_V, 判断风速是否连续性的参数V_vave, 在模糊识别出可能是杂波的情况下, 判断在一个有限小区域内的风速是否具有离散性和低值性的参数V_ratio, 判断是否是孤立点的参数N_X, 另外还有小区域内的径向速度平均值v_MD, 速度谱宽平均值w_MD, 算背景场时需先计算平均回波强度Z_MD等。主要的物理量定义如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

式中n_a, n_r表示在方位和距离方向定义的计算范围, 一般取值为7~11, Z_i,j, V_i,j为任意点的回波强度和速度, Z_up, Z_low分别表示上下层的PPI回波强度, Y(r), W(r) 表示与距离有关的权重, V_v表示单位距离 (默认为1 km) 风速的平均变化, S_ign(v), V_alid(v), F_og(v) 分别判断风速是正是负, 是否有效及是否为退模糊失败数据, V_vave则是判断在PPI上风速值是否连续变化, 其中A的取值为10~15 m/s, V_ratio则是对小区域内疑似杂波进一步分析其风速在给定区域内风速值在零风速附近的比例大小, V_sign则是用来判断风速是否是接近零的值, 其中V_thresh的取值在1~3 m/s。N_X是判断孤立点的参数, 对于给定的反射率库X, 以它为中心, 取M_total为5×5数据, 得到非空反射率探测值个数M, 两者相除即得非空值的百分比。

据统计分析, 层状云降水的回波特征表现为轻到中等的回波强度 (大多数情况下小于40 dBz), 弱的水平、垂直反射率梯度及范围较大 (通常远大于500 km²), 这有别于地物杂波的水平、垂直反射率梯度大, 面积较小等特征。而对于强对流回波, 相对来说在水平和垂直方向上回波强度的梯度变化较大, 但其回波顶高一般远大于5 km, 且其风速值相对较连续分布且值较大, 而且谱宽值也相对较大, 这也是有别于地物杂波接近零风速的分布及1.5°仰角PPI上谱宽值小于2 m/s的特征。因此根据降水回波及地物杂波回波特征的不同, 在统计的基础上, 设定一系列的参数阈值以滤除杂波。下面介绍杂波识别流程。

①背景场的构建

资料采用的是极坐标格式的基数据, 为避免降水回波误滤及能针对不同背景场里的杂波进行不同方式的滤波, 以得到更好的滤波效果。首先对资料进行一个简单的分类, 按回波强度的大小和垂直梯度分成未知类型, AP杂波, 杂波点, 弱降水或者晴空, 普通降水, 强回波边缘, 强回波等类型, 并进行标记, 统计给定区域里出现概率最大的回波类型, 从而得出该区域内的背景数据场。这里要说明的是只要0.5°和1.5°两仰角面对应的垂直变化在50 km内小于5 dB和50 km外小于10 dB时就认为是降水回波并予以一定程度上的保护。

②对孤立点的滤除

这主要是针对移动的物体 (船、飞机、鸟群、昆虫群等) 造成的孤立点及斑点的噪声回波。判断孤立点系数N_X, 如果N_X小于给定的阈值 (可取0.75) 则有条件滤除, 这个参数对于识别斑点的噪声回波是非常有效的, 但有时对降水回波边界上折角处的回波也会产生一定的平滑作用。

③联合各物理量参数进行识别过滤

根据强杂波的水平变化大及在垂直方向上迅速递减这一规律, 设置0.5°仰角面的Z_H及0.5°和1.5°仰角面间的Z_V的阈值, 如Z_H随距离的不同阈值设为25~40 dB, Z_V设为20~40 dB, 为避免小尺度强对流降水被误滤, 又增加了速度的连续性检测以及是否是零风速的判断, V_vave的值在大于阈值 (可取2) 附近时表示风速为不连续分布, V_ratio < 10%表示小区域内风速值较小且不连续等。对于超强的AP杂波, 在1.5°仰角上有时还会有中等强度的杂波存在, 这些杂波对应的风速谱宽值都是接近零的数值, 因此又可增加判断1.5°仰角上对应的w_MD, v_MD的绝对值是否是小值, 给其设定的值分别为小于2 m/s和1 m/s。

另外, 为了更快地发现可疑区域, 分析是否是杂波, 提出了突变区域的过滤方法。这种算法首先需要确定边界, 对于大于30 dBz回波点, 计算其上层对应点周围7×7网格的平均回波强度Z_MD, 当该点与Z_MD的Z_V大于10 dB时, 如果该点与周围8个点的任意Z_H大于给定的阈值M (15~25 dB), 则搜索下一个点, 直到另一个点与周边相邻的任一点Z_H大于M时结束。最后对这个回波团进行突变区域的检测, 如果回波团的点数小于阈值 (可取50~100), 且回波团的突变边界的比例系数 (指回波团边上的点与周边点的反射率差值超过阈值N, 则认为该点存在突变, 统计回波团边界上这种点所占的比例即为突变边界系数), 达到一定阈值, 则认为该区域可能为可疑杂波团, 经进一步分析, 满足一定的严格条件后, 如这团回波w_MD, v_MD的绝对值是否是小值, 或者其背景场的回波强度是否在15到30 dBz的范围内 (这主要是考虑到小尺度强对流回波团与周边的回波强度也可能存在梯度较大的情况), 才考虑对该区域进行滤波。

④杂波处理

在杂波被识别出后, 依据前面所提出的原则理念, 对于杂波的处理应该结合直接全部剔除和近似数值填充机制等多种处理方式。采用这样的杂波处理策略的好处是可以把误滤所造成的损失降到最低, 这一方面对数据起到了保护作用, 另一方面也在一定程度上还原了真实的回波场。例如对于混在降水中的杂波, 本文是综合相邻仰角层以及周边数据经过线性插值后再进行保守性的填充, 而对于无降水区域中的杂波, 识别后则采用直接去除的方法。

本文抽样测试了2004—2006年共190组基数据, 分类分距离按级别进行效果检验, 结果如表 1所示。表中A表示无降水回波, B表示杂波与普通降水混合或相邻的个例, C表示杂波与对流性降水混合或相邻的个例, D则表示真正的降水回波没有被误滤的情况, 字母后的1, 2, 3, 4分别表示距离R的范围: 0 < R≤50 km, 50 < R≤100 km, 100 < R≤150 km, R >150 km等4种情况。其效果检测为“优秀”表示在该范围内的杂波都被明显滤除, “良好”指部分杂波没有清除干净, “失败”则指该范围内的个别杂波明显没有被滤除。而针对D组, “优秀”表示降水回波明显没误滤, “良好”指30 dBz以下的弱回波边缘有部分被滤掉, “失败”则指个别小尺度回波团或40 dBz以下的回波边缘被明显误滤。从表中可看出, 在无降水回波的情况下, 150 km内杂波的优秀识别率达到94 %以上, 失败的个例很少。而在大于150 km之外, 则滤波效果很差。这是因为在识别过程中始终结合速度和谱宽等进行处理, 而在接近150 km和150 km之外, 速度和谱宽产品的数据很不可靠且速度模糊的现象明显增多, 这也同样导致了B, C类型在远距离杂波识别中失败率的增加。在杂波与普通降水混合或相邻的情况下, 在100 km范围内没有明显失败的个例, 优秀率也较高, 而在100 km之外, 失败的个例明显增加。在杂波与对流性降水混合的情况下, 杂波的识别率还要优于前者, 这主要是因为对流性降水回波高度较高, 远距离处在1.5°仰角上依然还有回波存在, 另外, 从现有的个例来看, 出现的AP杂波强度也相对更强, 而普通降水情况下对应的上层不一定有回波, 另外, 混合的AP杂波强度相对较弱。滤波方案在综合参考上层的回波情况及相关背景场等情况下, 强对流回波滤波效果更好。分析降水回波没有被误滤的情况, 由表 1可看出, 在100 km之内, 基本不存在降水被误滤的情况, 而超过100 km, 从个例分析来看, 部分强对流回波边缘突起部分有小于30 dBz的回波被滤掉, 从算法上来分析, 主要是被孤立点算法所滤除, 另外有个别强度小于40 dBz小尺度的回波团由于在远距离处风速质量较差, 有回波处对应无效数据, 从而引起部分误滤。但总的来说, 在150 km内, 滤波工作还是很有效的, 下面给出相关的个例分析。

表 1

表 1 3 种情况下不同距离范围内的滤波效果对比分析与降水误滤分析 Table 1 The contrast analysis on filtering impact in different ranges under the three instances and the percentage of reserved real precipitation

表 1 3 种情况下不同距离范围内的滤波效果对比分析与降水误滤分析 Table 1 The contrast analysis on filtering impact in different ranges under the three instances and the percentage of reserved real precipitation

彩图 1b中显示的都为晴空条件下出现的杂波, 西北面大片都为AP杂波, 在对应的速度图上, 在150 km范围内对应的都是零风速区, 谱宽值也相对较小。而在150 km之外则是由于距离折叠造成的退模糊失败的区域, 部分速度值超过60 m/s, 数据质量显然不可靠。在1.5°仰角上杂波强度和范围迅速减小。雷达站的南面50 km附近的地物杂波强度超过50 dBz, 对应的速度图上表现为离散的零风速分布, 谱宽值相应较小。在1.5°仰角上这种杂波基本上没有。在雷达站50 km内则是相对较弱杂波, 回波强度大多小于15 dBz, 对导出产品的影响不大, 在这里都不作滤除。从彩图 1a, e可看出, 在150 km之内, AP杂波和地物强杂波都得到了有效的滤除。而在150 km之外, 由于速度质量的不可靠而不在质量控制流程之内而得以保留下来。

图 1

图 1. 2016年6月13日宁波雷达05:50（北京时，下同）径向速度图及滤波前后反射率因子产品对比图 （距离圈相距50 km） Fig 1. The Nigbo radar products of radia velocity and refelctivity before and after the filtering at 05:50 on June 13, 2006. (the distance between adjaceat circles is 50 km)

如彩图 2b所示, 在雷达的西北面, 东北—西南向的带状区域为AP杂波, 其西南面还有降水回波相邻。在雷达的南面零星分布着点状或小团块状的则为相对静止的地物杂波。质量控制算法在这个个例中成功删除了地物强杂波, 只是在150 km之外在发生速度模糊的区域还残留着些许杂波。在滤波方案中, 在远距离处, 虽然回波的Z_H, Z_V较大, 但为了避免将小尺度的强对流回波团误滤, 在做速度检测时, 碰到有速度模糊的区域会把它保留下来, 这也是残留杂波没有滤除的原因, 仍需改进。

图 2

图 2. 同图 1，但是2006年4月10日06:55 Fig 2. Same as in Fig.1, but for 06:55 on April 10, 2006

如彩图 3所示, 在雷达西北面100 km附近的区域为弱的AP杂波, 而西南面50 km左右的NP杂波混在20~40 dBz的降水回波里, 南面还有离散的点状或小团块状的强杂波镶嵌在降水或无降水区域, 在1.5°仰角的反射率因子产品上看不到有杂波分布, 经过质量控制, 杂波在0.5°上得到了有效的抑制。

图 3

图 3. 同图 1，但是2006年4月10日11:54 Fig 3. Same as in Fig.1, but for 11:54 on April 10, 2006

在彩图 4b中, 在雷达站北面100 km之外有中等强度的AP杂波混合在降水回波中, 从对应的速度图上看, 降水的零风速带和杂波的零风速相交。应用突变区域加零风速等滤波方法并没有产生误滤的现象。从彩图 4a中可看出滤波效果相当好。

图 4

图 4. 同图 1，但是2006年7月20日05:55 Fig 4. Same as in Fig.1, but for 05:55 on July 20, 2006

在彩图 5中, 在雷达的西北面为超强的AP杂波, 回波强度大于90 dBz, 左侧相邻的是强对流回波。对应的速度产品上表现为:在150 km内, 0.5°仰角上对应为大片的零风速区, 在150 km之外则强弱相间, 杂乱无章。谱宽产品在0.5°仰角上表现得有些离散, 在150 km之内大部分小于2 m/s, 在150 km之外则由于速度模糊, 表现得相当杂乱。而在1.5°仰角的产品中, 速度图上杂波对应的地方基本为零, 谱宽对应的也都是小值 (谱宽图略)。因此在滤波过程中, 更好地运用1.5°仰角的速度和谱宽的杂波特征, 有利于低层杂波的滤除。经过质量控制后, 杂波基本上得到了滤除。小尺度的强对流回波团在滤波过程中也没有发生误滤的现象, 说明滤波方法是安全可靠的。

图 5

图 5. 同图 1，但是2006年6月8日23:00 Fig 5. Same as in Fig.1, but for 23:00 on July 8, 2006

1) 宁波雷达数据造成污染的主要原因是普通地物杂波和AP杂波。这些杂波总体表现为在低仰角的回波强度总体呈离散状分布, 强弱交杂, 垂直变化大, 0.5°仰角面的风速是间杂的接近零的数值, 谱宽是无序的强弱分布的离散值。而在1.5°仰角面上, 风速对应的值基本上等于零, 谱宽较小。因此在程序设计时, 可以在分析反射率因子数据基础上, 依次参考速度和谱宽数据, 对其进行滤波。

2) 实际设计中, 可根据经验总结出特征规律, 并考虑多种情况, 构建一些简化之后的线性关系并设置必需的阈值, 从而达到人工智能的效果。

3) 人工智能识别滤波方法经过190个例分类、分距离、分级别进行验证, 总体效果优异, 特别适用于150 km之内, 在确保数据的完整、准确的基础上能够有效地过滤非气象杂波。但在150 km之外, 由于小尺度的对流回波的Z_H, Z_v较大以及远距离的速度质量较差等客观因素, 导致对部分AP杂波过滤效果不理想。

在算法设计过程中, 坚持杂波每种具体特征的滤波算法相对独立, 因此其具有非常好的可拓展性。在未来如果可以融入卫星图像数据和地面实况数据, 增加反射率因子产品在时间上的连续性检测等, 理论上可以把滤波的有效范围拓展到200 km甚至更远。