一种用于气象学中滤除超折射的方法

摘要

本发明公开了一种用于气象学中滤除超折射的方法，滤除多普勒气象雷达0.5°仰角反射率图像中小于35dBz的像素点，获取第一滤除图像；滤除面积小于50个单位的连通域，获取第二滤除图像；采用半径为1的圆形结构元素进行形态学开运算后再滤除面积小于50个单位的连通域，获取第三滤除图像；通过灰度共生矩阵得到第三滤除图像中每个连通域的惯性矩特征；滤除面积小于100且惯性矩大于1.9和面积大于100且惯性矩大于1.5的连通域，获得第四滤除图像；滤除径向速度V∈[-1，1]的像素点，获取第五滤除图像；采用半径为1的圆形结构元素进行形态学闭运算。在滤净超折射的同时，保留了有效降水云团，使对超折射的漏识率和误判率更小。

说明书

技术领域

本发明涉及气象领域，特别涉及一种用于气象学中滤除超折射的方法。

背景技术

多普勒天气雷达是利用云雨目标物对电磁波的后散射来发现气象目标，并测 定气象目标物的空间位置和强弱分布，但多普勒天气雷达同样能探测到非气象目 标，这些非气象目标回波会污染多普勒天气雷达数据，使天气雷达三个基本的物 理参数平均功率、平均多普勒速度和谱宽的估计出现较大的偏差，所以非气象目 标回波抑制成为雷达数据质量控制的重要问题。

非气象目标回波包括雷达波束正常传播(NP)情况下探测到的地物回波和异 常传播(AP)情况下探测到的超折射回波。特殊气象条件如低空出现逆温大气层会 引发雷达波速异常传播，而特殊的气象条件每天、每小时、甚至每个探测时刻都 在变化，因此，超折射是多年来研究人员致力于消除的主要对象。

在现有的技术中，超折射滤除效果较好的是由江源^[1]改进的Kessinger^[2]方法， 该方法首先构建反映超折射和降水回波差异的6个参量，即回波强度的水平变化 T_DBZ、垂直变化G_DBZ等；然后根据超折射的分布参量构建6个折线型隶属函数， 并借此确定6个判断阈值，最后进行组合判断。

发明人在实现本发明的过程中，发现现有技术中至少存在以下不足：

现有技术中选用的6个参量的类间混叠比较严重，虽然使用了隶属函数及多 参量的组合判断，依然在没有滤净超折射的同时，过多地滤除了有效降水云团； 同时该方法特别要求两类样本的量要足够的多，且具有足够的代表性，在这些不 满足时，4％的漏识率则偏于乐观，而4％的误判率仍然较高。

发明内容

本发明提供了一种用于气象学中滤除超折射的方法，本发明避免了过多滤除 有效降水云团，降低了漏识率，详见下文描述：

一种用于气象学中滤除超折射的方法，所述方法包括以下步骤：

一种用于气象学中滤除超折射的方法，所述方法包括以下步骤：

(1)滤除多普勒气象雷达0.5°仰角的反射率图像中小于35dBz的像素点，获取 第一滤除图像；

(2)从所述第一滤除图像中滤除面积小于50个单位的连通域，获取第二滤除 图像；

(3)对所述第二滤除图像采用半径为1的圆形结构元素进行形态学开运算， 获取运算后图像；

(4)从所述运算后图像中滤除面积小于50个单位的连通域，获取第三滤除图 像；

(5)将所述第三滤除图像中的每一个连通域作为生成灰度共生矩阵的对象， 计算其惯性矩特征；滤除面积小于100且惯性矩大于1.9和面积大于100且惯性 矩大于1.5的连通域，获取第四滤除图像。

(6)从所述第四滤除图像中滤除径向速度V∈[-1，1]的像素点，获取第五滤除图 像；

(7)对所述第五滤除图像采用半径为1的圆形结构元素进行形态学闭运算， 获取不包含超折射的反射率图像。

所述将所述第三滤除图像中的每一个连通域作为生成灰度共生矩阵的对象， 计算其惯性矩特征；滤除面积小于100且惯性矩大于1.9和面积大于100且惯性 矩大于1.5的连通域，获取第四滤除图像具体为：

1)反射率图上连通域区域的灰度共生矩阵的构建

灰度共生矩阵P_d，θ反映数字图像f中θ方向上相距d的两像素间灰度变化的 情况，令d＝1，θ为0°或45°或90°或135°，针对一个反射率图像中的一个连通域Ω， 得到4个灰度共生矩阵，记为P_1，0°、P_1，90°、P_1，45°和P_1，135°，各矩阵元素的计算公式如 下：

q_1，0°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，k＝m，|l-n|＝1，(k，l)∈Ω，(m，n)∈Ω}      (1)

q_1，90°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，|k-m|＝1，l＝n，(k，l)∈Ω，(m，n)∈Ω}     (2)

q_1，45°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，k-m＝±1，l-n＝m1，(k，l)∈Ω，(m，n)∈Ω}  (3)

q_1，135°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，k-m＝ml，l-n＝±1，(k，l)∈Ω，(m，n)∈Ω} (4)

其中，灰度共生矩阵的阶数就是对象的灰度等级，记为N，i，j＝0，1，L，N-1， &表示括号内条件成立的像素对数目，多普勒雷达反射率图像中的各像素点取值 (反射率强度)使用16个灰度等级，从0～15灰度级表示的反射率强度依次为背景 色、-5dBz、0、5dBz、10dBz、15dBz、20dBz、25dBz、30dBz、35dBz、40dBz、 45dBz、50dBz、55dBz、60dBz、65dBz及其以上；

2)反射率图上连通域区域的灰度共生矩阵

3)计算惯性矩特征T_I

T_I＝(∑_{θ＝0°，90°，45°，135°}I_θ)/4

$I_{\theta }={\Sigma }_0^{N-1}{\Sigma }_0^{N-1}{(i-j)}^2{p}_{d,\theta }(i,j)$

$p_{d,\theta }(i,j)=\frac{q_{d,\theta }(i,j)-\min \left\{q_{d,\theta }(i,j)\right\}}{\max \left\{q_{d,\theta }(i,j)\right\}-\min \left\{q_{d,\theta }(i,j)\right\}}$

4)滤除面积小于100且T_I＞1.9和面积大于100且T_I＞1.5的连通域，获取所述第 四滤除图像。

本发明提供的技术方案的有益效果是：本方法首先找到基于灰度共生矩阵的 与旋转无关的高类间区分度的主力特征，配合使用数学形态学的方法及来自径向 速度的辅助特征，经多层滤波及基于有效特征的识别，完成超折射的滤除。使得 在保留了有效降水云团的前提下，几乎滤净了超折射，使得误判率从大于4％降 低到1.01％。

附图说明

图1为本发明提供的一种用于气象学中滤除超折射的方法的流程图；

图2为本发明提供的生成灰度共生矩阵的有效区域示意图；

图3为本发明提供的超折射与降水回波混合样本集关于3个纹理特征的分 布图；

图4为各类样本在特征惯性矩T_I上的分布图；

图5-1、图5-2、图5-3、图5-4、图5-5、图5-6和图5-7为本发明提供的超 折射滤除结果示例图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实 施方式作进一步地详细描述。

为了在尽量滤除超折射的同时，避免过多滤除有效降水云团，降低漏识率和 误判率，本发明实施例提供了一种用于气象学中滤除超折射的方法，参见图1， 详见下文描述：

就反射率强度及强度分布特点而言，对暴雨、冰雹等灾害性降水云团回波与 超折射回波进行比较，不难发现以下特点：

1、整体上，超折射通常呈现为具有多个高反射率中心的间歇性小团状杂波 群，各小团或放射性分布、或无规则分布；降水云团则多为具有一个或多个高反 射率中心的单核个体或多核个体，其中，部分多核个体与超折射小团在其分布及 区域面积上具有相似性；

2、超折射小团内及边缘处的反射率强度跳变剧烈；降水云团的反射率强度 却逐次递增或递减，几乎不发生跳变；

3、有部分超折射区域的径向速度值分布在[-1，1]m/s之间；降水云团的径向 速度值往往超出此范围。

根据以上降水回波区域与超折射区域之间不同的图像特点，为了更加有效的 去除超折射，该方法包括以下步骤：

101：滤除多普勒气象雷达0.5°仰角的反射率图像中小于35dBz的像素点， 获取第一滤除图像；

考虑到多普勒天气雷达回波数据的可视化软件通常把反射率数据离散化为 15个等级，即反射率值R＝{-5，0，5，10，15，20，25，30，35，40，45，50，55，60，≥65}dBz，一般来讲形 成灾害性天气的云团需达到一定的反射率强度，反射率强度小于35dBz的地方不 会形成灾害性天气。经验数据为风暴云团一定存在反射率值R≥35dBz的区域，另 外，超折射多出现在低空，多见于0.5°仰角的反射率图像之中，所以滤除多普 勒气象雷达0.5°仰角的反射率图像中小于35dBz的像素点，获取第一滤除图像， 以减小计算区域。

102：从第一滤除图像中滤除面积小于50个单位的连通域，获取第二滤除图 像；

由于形成灾害性天气的云团在低仰角面积的经验数据是大于等于50km²， 为此从第一滤除图像中滤除面积小于50个单位(1个单位是1个像素点，1个像 素点代表1km)的连通域，获取第二滤除图像。

103：对第二滤除图像使用半径为1的圆形结构元素进行形态学开运算，获 取运算后图像；

其中，形态学开运算具体为先腐蚀后膨胀的过程，它具有消除细小物体、在 纤细点处分离物体和在平滑较大物体的边界时又不明显改变其面积的作用。在滤 除面积小于50个单位的连通域之后，有可能出现在被保留下来的区域之间出现轻 微粘连，所以要用半径为1的圆形结构元素进行形态学开运算，将其分开。

104：从形态学开运算后图像中滤除面积小于50个单位的连通域，获取第三 滤除图像；

其中，在对第二滤除图像进行形态学开运算之后，有的区域面积可能会减小， 所以要再次滤除面积小于50个单位的连通域，获取第三滤除图像。

105：将第三滤除图像中的每一个连通域作为生成灰度共生矩阵的对象，计 算基于灰度共生矩阵的惯性矩特征T_I；利用惯性矩特征，滤除面积小于100且 T_I＞1.9和面积大于100且T_I＞1.5的连通域，获取第四滤除图像；

其中，该步骤具体为：

1)多普勒反射率图像中连通区域的灰度共生矩阵的构建；

灰度共生矩阵P_d，θ反映数字图像f中θ方向上相距d的两像素间灰度变化的 情况，在d和θ确定之后，共生矩阵的元素值q_d，θ(i，j)便等于满足条件f(k，l)＝i、 f(m，n)＝j的像素对数目。设待分析的连通区域为Ω。为充分反映超折射区域内反 射率强度跳变的情况，令d＝1，θ分别取0°、90°、45°、135°，得到4个灰度共生矩阵P_1，0°、 P_1，90°、P_1，45°和P_1，135°，各矩阵元素的计算公式如下：

q_1，0°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，k＝m，|l-n|＝1，(k，l)∈Ω，(m，n)∈Ω}      (1)

q_1，90°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，|k-m|＝1，l＝n，(k，l)∈Ω，(m，n)∈Ω}     (2)

q_1，45°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，k-m＝±1，l-n＝m1，(k，l)∈Ω，(m，n)∈Ω}  (3)

q_1，135°(i，j)＝&{f(k，l)＝i，f(m，n)＝j，k-m＝ml，l-n＝±1，(k，l)∈Ω，(m，n)∈Ω} (4)

其中，灰度共生矩阵的阶数就是对象的灰度等级，记为N，i，j＝0，1，L，N-1， &表示括号内条件成立的像素对数目。多普勒雷达反射率图像中的各像素点取值 (反射率强度)使用16个灰度等级，从0～15灰度级表示的反射率强度依次为背景 色、-5dBz、0、5dBz、10dBz、15dBz、20dBz、25dBz、30dBz、35dBz、40dBz、 45dBz、50dBz、55dBz、60dBz、65dBz及其以上。因此由式(1)～(4)定义的元素组 成的四个灰度共生矩阵为：

基于灰度共生矩阵的纹理特征有14种之多^[3]，其中最常用的是共生矩阵的能 量E、熵H和惯性矩I，计算方法如下：

$E={\Sigma }_0^{N-1}{\Sigma }_0^{N-1}{p_{d,\theta }}^2(i,j)---\left(5\right)$

$H=-{\Sigma }_0^{N-1}{\Sigma }_0^{N-1}{p}_{d,\theta }(i,j)\times \log p_{d,\theta }(i,j)---\left(6\right)$

$I={\Sigma }_0^{N-1}{\Sigma }_0^{N-1}{(i-j)}^2{p}_{d,\theta }(i,j)---\left(7\right)$

式中，P_d，θ(i，j)是被归一化了的q_d，θ(i，j)：

$p_{d,\theta }(i,j)=\frac{q_{d,\theta }(i,j)-\min \left\{q_{d,\theta }(i,j)\right\}}{\max \left\{q_{d,\theta }(i,j)\right\}-\min \left\{q_{d,\theta }(i,j)\right\}}---\left(8\right)$

设图2是一个待识别的连通域，它的边界一般为35dBz，灰度值取9，它与相 邻的背景0值出现很大的跳变，这种跳变是因以35dBz做分割阈值造成的，因此是 不真实的，为防止这种不真实的跳变混淆超折射区域中灰度值之间真实的跳变， 同时，不遗漏超折射区域35dBz向下的剧烈跳变，沿区域边界外扩一个单位，如 图2所示。在计算灰度共生矩阵中各元素的取值时，对像素点的区域约束从 (m，n)∈Ω改成(m，n)∈Ω′。

2)基于灰度共生矩阵的能量E、熵H和惯性矩I惯性矩特征T_I分类有效性分 析。

从图3来看，超折射与降水回波混合样本集在特征T_E上没有形成波谷，说明 此特征的类间区分度极差；在特征T_H上虽出现波谷，但因其两侧的样本数悬殊较 大，说明样本群在该特征上表现出较大的混叠，类间区分度不高；与特征T_I的谷 值1.5相对应，左右两个子样本群的样本数目与表1数据基本匹配，预示着特征T_I可 能具有较强的类间区分度。于是，分别用4593个超折射样本和4155个降水样本(面 积超过100单位的2326个、小于100单位的1829个)统计各自在特征T_I上的分布， 如图4所示。

从图4不难看出，就统计样本而言，降水区域和超折射区域这两类样本在 T_I＝1.5～1.9的区域上发生了混叠，而这部分样本数为709个，占总样本数的8.1％， 再观察这8.1％的样本，发现有相当的数量属于超折射区域与降水区域发生轻微粘 连而被作为一个连通域样本的情形，如果对这种“样本”进行有效地拆解，混叠 的样本数会低于8.1％。将图4与图3进行对比，会发现判别阈值1.5正好与图3中T_I分布的波谷值相对应。如果以面积小于100且惯性矩大于1.5和面积大于100且惯 性矩大于1.9为条件来进行区域滤除，会有低于8％的超折射区域被遗漏，而有效 降水区域不会被误删除。

“超折射小团内及边缘处的反射率强度跳变剧烈，降水云团的反射率值却逐 次递增或递减，几乎不发生跳变”，结合这一特点分析由式(5)～(7)定义的3个纹 理特征，可以解释如上实验结果。具体的，特征T_E和T_H反映的是区域内相邻反 射率强度变化的次数而不是变化的程度，因此，它们无法体现相邻反射率强度的 跳变情况；特征T_I将变化次数乘以因子(i-j)²，对于渐进变化的降水样本来讲， 此因子约为1，而对于超折射样本，此因子则几乎总是大于1，且跳变程度越高， 大于1的程度越大，例如，某超折射区域内的一对相邻像素的反射率值分别为 55dBz和40dBz，则(i-j)²＝(13-10)²＝9＞＞1。

以上理论分析及统计实验均说明特征T_I最适合于担当鉴别超折射的工作。

106：从第四滤除图像中滤除径向速度V∈[-1，1]的像素点，获取第五滤除图像；

由于第四滤除图像中还会有低于8％的残留超折射，利用降水回波和超折射 在径向速度值上的差异，即有部分超折射区域的径向速度值分布在[-1，1]m/s之间 而降水云团的径向速度值往往超出此范围，所以再滤除径向速度V∈[-1，1]的像素 点，获取第五滤除图像。

107：对第五滤除图像采用半径为1的圆形结构元素进行形态学闭运算，获 取不包含超折射的反射率图像。

由于进行上述处理之后，连通域中的元素可能会出现细小空洞，用半径为1 的圆形结构元素进行形态学闭运算以弥补因上步操作造成的小空洞，得到系列处 理后的反射率图像。

下面以一个具体的试验来验证本发明实施例提供的一种用于气象学中滤除 超折射的方法的可行性，详见下文描述：

历史样本中，2006年6月27-28日天津地区的多普勒雷达回波图像中连续出现 了比较严重的超折射，表1给出了这组样本的细节信息。

表1 统计实验的样本信息

参见图5-1、图5-2、图5-3、图5-4、图5-5、图5-6和图5-7，图5-1描述了 多普勒气象雷达0.5°仰角的反射率图像中大于35dBz的像素点的反射率图；图5-2 中描述了反射率大于35dBz和连通域面积大于50个单位的反射率图；图5-3对 图5-2中的反射率图用半径为1的圆形结构元素进行形态学开运算之后的图；图 5-4进行完上述的形态学开运算以后再滤除面积小于50个单位的连通域；图5-5 滤除面积小于100且惯性矩大于1.9和面积大于100且惯性矩大于1.5的连通域 之后的反射率图；图5-6在图5-5的基础上再滤除径向速度V∈[-1，1]的像素点之后 的反射率图像；图5-7在图5-6的基础上再用半径为1的圆形结构元素进行闭运 算(补空洞)，得到最终的反射率图像。从图5-7中可以得出在有效降水区域全部 保留的前提下，残留了极少量的超折射区域。

测试图像选自天津地区2006.6.27-28、2007.7.11、2007.8.1和8.26日的同时含 超折射和降水回波的反射率图1193张，运用本文方法处理后，得到测试结果如表 2所示。

表2 识别率统计

理论分析及统计实验说明了选用本方法的合理性和有效性。与使用6特征协 同隶属度鉴别的方法^[1]相比，错误率从大于4％降低到不足1.01％，同时，文献 ^[1]的接近4％的错误率是以像素点做样本量的统计(共计42372个超折射像素点)、 且包含了对有效降水的漏识，于是形成滤除超折射的较大损失，而本方法做到的 是低错误率下的1.01％损失，它以连通域做样本，超折射像素点数的规模超过文 献^[1]的14982×(100～10000)/42372＝35～3500倍，使测试结果的置信度更高。

综上所述，本发明实施例提供了一种用于气象学中滤除超折射的方法，本方 法首先找到基于灰度共生矩阵的与旋转无关的高类间区分度的主特征，配合使用 数学形态学的方法及来自径向速度的辅助特征，完成对超折射的滤除。使得在保 留了有效降水云团的前提下，几乎滤净了超折射，使对超折射的误判率和漏识率 更小。

参考文献

[1]江源，刘黎平，庄薇.多普勒天气雷达地物回波特征及其识别方法改进[J].应用气 象学报，2009，20(2)：203-212；

[2]Kessinger C，Ellis S，Vanandel J，et al.The AP Clutter Mitigation Scheme f or the  WSR-88D//Preprints，31st Conference on Radar Meteorology，Amer Meteor Soc， 2003：526～529.

[3]RMHaralick，K Shangmugam，etc.Texture feature for image classification[J]. IEEE Transaction on Systems，1973，SMC-3(6)：768-780.

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实 施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精 神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护 范围之内。