CN102254191B - 一种基于图像处理的降水粒子相态识别方法 - Google Patents

Abstract

一种基于图像的降水粒子相态识别方法，属于数字图像识别领域，目的在自动检测降水粒子目标的基础上，识别并判定其相态，并进一步提取相关降水信息，以提高降水天气现象观测的自动化程度与精度。本发明用于降水天气现象观测中降水相态的自动判别与降水信息获取。顺序包括：(1)图像获取步骤，(2)降水粒子目标检测步骤，(3)特征提取步骤，(4)降水粒子相态识别步骤，(5)降水信息提取步骤。本发明结合图像处理和模式识别的方法能有效识别降水粒子相态，判别降水天气现象类型，从而提高降水天气现象观测的自动化程度与精度。

Description

一种基于图像处理的降水粒子相态识别方法

技术领域

本发明属于数字图像识别领域，具体涉及一种基于图像的降水粒子相态识别方法，用于自动识别降水粒子的相态，进而提取相关的降水信息。

背景技术

降水观测是天气现象观测中的重要组成部分，在国民经济建设服务中有着重要的作用，如防洪、抗旱、减灾、水利工程的设计、工农业生产等，同时对气候变化特别是季风研究等许多方面有着重要的影响。

在降水观测中很重要的一点就是实时得到降水粒子的相态、形状、尺度信息，估测降水强度和降水量，判别降水天气现象类型，提供准确及时的天气状况情报。

降水粒子相态可以分为两类，一类是液态粒子，如雨滴；一类是固态粒子，如雪花、冰粒等。考虑到实际降水过程中，液态粒子主要为雨滴，固态粒子主要为雪花的情况，本发明仅考虑识别降水粒子相态为雨滴或雪花。

目前国内外识别降水粒子相态的方法是获取具有双线偏振功能的雷达探测到的水平反射率因子(Z_H)、垂直反射率因子(Z_V)、差分反射率因子(Z_DR)、线性退偏振比(L_DR)、差传播相移(Φ_DP)、零滞后互相关系数(ρ_HV(0))等与降水粒子类型、形状和分布等要素密切相关的参数进行反演分析，推导降水粒子的相态。我国在利用偏振雷达识别降水粒子相态方面的研究工作不多，主要原因是我国偏振雷达技术的发展相对滞后。由于国外不同雷达的测量参数不尽相同，而且识别出的天气类型与我国的天气现象有较大的差异，因此这些方法在我国不能完全使用。目前国内关于降水粒子相态识别的研究还不成熟，缺乏实际可应用的用于识别降水粒子相态的系统。

1962年，美籍华人胡贵明首先提出了首先提出7个基于直角坐标系的几何矩(以下分别用φ_Hu₁、φ_Hu₂、φ_Hu₃、φ_Hu₄、φ_Hu₅、φ_Hu₆和φHu₇表示)，并应用于形状识别。

Chen提出了利用区域边界计算区域矩的快速算法(以下分别用φ_Chen₁、φ_Chen₂、φ_Chen₃、φ_Chen₄、φ_Chen₅、φ_Chen₆和φ_Chen₇表示)，该算法通过图形的边缘点来计算区域矩。

相对矩指的是利用矩之间的比值去掉与区域有关的面积比例因子，从而使不变矩(以下用R₁、R₂、R₃、R₄、R₅、R₆、R₇、R₈、R₉和R₁₀表示)与面积或结构的比例缩放没有关系，而仅与几何形状有关。

二维极坐标傅立叶描述子能够捕捉形状的径向特征和角向特征，可以以一定的精度描述目标的形状特征，可以对形状的相似性进行定量的判别，其干扰能力强，通用性高。

特征选择是模式识别中的重要过程，现有的特征选择方法有分支定界法、模拟退火算法、遗传算法、顺序前进法、顺序后退法等方法。

中国专利文献201010145822.0公开了一种基于梯度的雨滴识别方法，该方法通过将待识别雨滴图像及背景图像分别进行各向异性扩散平滑滤波，结果图像进行差分操作，获取差分图像，再对差分图像提取边缘梯度并进行二值化，经过形态学操作处理，即得到雨滴目标检测结果。

上述公开的方法中，其只能识别出雨滴，没有对降水粒子相态进行区分，对于其他相态的降水粒子如雪花，也只能识别为雨滴，即该方法存在一定局限性。

发明内容

为克服现有方法中所存在的问题，本发明提供一种基于图像的降水粒子相态识别方法，能够自动识别降水图像中的降水粒子目标的相态并提取相关的粒子信息，提高了降雨天气现象观测的自动化程度与精度。

本发明的一种基于图像处理的降水相态识别方法，包括如下步骤：

(1)图像获取步骤，分别从成像设备获取大小一致的待识别的降水图像及其所对应的背景图像；

(2)对获取的图像进行降水目标检测，获得降水目标检测结果图像。

将待识别降水图像及背景图像分别进行各向异性扩散平滑滤波，得到的两幅结果图像进行差分操作，获取差分图像。对差分图像提取边缘梯度并进行二值化，经过形态学操作处理，得到降水目标检测结果图像。

(3)对降水目标检测结果图像中每个目标，提取其图像特征，包括对比度、凹凸度、致密度、圆形度、Hu矩、Chen矩、相对矩、二维极坐标傅立叶描述子。

(4)降水目标识别步骤，

首先，训练SVM分类器，即先提取训练用的目标样本的图像特征值，并根据该目标样本对应的相态类别，训练SVM分类器。

然后，提取待识别降水目标的图像特征值，作为训练好的SVM分类器的输入，得到分类结果，即为待识别降水目标相态的识别结果，包括雨滴或雪花。

本方法中，还可以包括对识别出来的结果进行进一步的信息提取，具体包括：提取识别结果中的雨滴目标数目、雨滴目标平面平均直径、雨滴目标立体平均直径分布谱、雪花目标数目、雪花平面面积。具体为：

(5.1)雨滴目标的数目Num_Raindrop，即R(x，y)中灰度值为255的连通域的数目；

(5.2)设雨滴目标Ri最小外接椭圆的长轴长度为MajorL_i，短轴长度为MinL_i，则该雨滴目标的平均平面直径MeanFlatDia_i定义为：

${\mathrm{MeanFlatDia}}_i=\sqrt{{\mathrm{MajorL}}_i\×{\mathrm{MinL}}_i}$

统计每个雨滴目标的最小外接椭圆的长轴长度MajorL_i和短轴长度为MinL_i，即可得出雨滴目标的平均平面直径MeanFlatDia_i，通过实验获取平面直径与立体直径的函数关系，得到雨滴目标的平均立体直径MeanSolidDia_i；

(5.3)设降水图像中雨滴目标平均立体直径的值域为RN＝[MinDia，MaxDia]，将RN平均分成Num1个子值域{RN_j，j＝1，2，…Num1-2，Num1-1，Num1}，雨滴目标平均立体直径分布谱DiaSpectrum描述的是雨滴平均立体直径MeanSolidDia_i分布在RN中各个子值域上的数目。具体实施中Num1取[20，30]间的整数；

(5.4)雪花目标的数目Num_Snowflake，即R(x，y)中灰度值为128的连通域的数目；

(5.5)雪花目标平面面积Area_Snowflake_j描述的是雪花目标S_j实际所占平面面积；

(5.6)统计雨滴目标总体积和雪花目标平面总面积，结合成像载玻片的面积和降水采集时间，得出降雨量、降雪量、降水量、降雨强度和降雪强度，根据国家气象局颁布的降水强度等级划分标准得出降水类型为小雨、中雨、大雨、暴雨、大暴雨、特大暴雨、雨夹雪、小雪、中雪、大雪、暴雪、大暴雪或特大暴雪。

本发明中，还可以在步骤(3)所提取的图像特征中，先进行特征选择以选择出识别最有效的若干特征，再进行降水目标识别步骤，以提高识别准确率。选择出识别最有效的若干特征可以使用现有常规的特征选择算法。

本发明将数字图像自动识别技术引入到了降雨天气现象观测中，可以自动识别降水图像中的降水粒子相态并提取相关的信息，为降水粒子相态的判别提供了一种新的方法。结合时间信息，本发明可进一步应用于降水量统计、实时降水强度监测以及降水相态与微结构分析等，在国民经济建设服务中有着重要的作用。

附图说明

图1是本发明流程示意图；

图2是一幅待识别的降水图像；

图3是待识别降水图像所对应的背景图像；

图4是降水图像中降水目标检测结果图像，图像中黑色部分表示背景，每一块白色区域表示一个降水目标；

图5是降水图像中降水目标粒子相态识别结果，图像中黑色部分表示背景，每一块白色区域表示一个雨滴目标，每一块灰色区域表示一个雪花目标；

图6示意了一幅雨滴图像的灰度分布；

图7示意了一幅雪花图像的灰度分布；

图8示意了若干规则和不规则的二值形状图像；

图9示意了点集的凸包。

具体实施方式

本发明以灰度和形状信息作为识别降水粒子相态的显著特征，通过结合图像处理和模式识别的方法能有效识别降水粒子相态并提取相关的降水信息，判别降水天气现象类型，处理流程图如图1所示。

(1)图像获取步骤，分别从成像设备获取尺寸大小一致的待识别的降水图像及其所对应的背景图像。背景图像是指按一定条件清除了观测系统中的所有降水粒子后，成像设备所拍摄的图像。图2所示的是待识别降水图像，图3所示的是图2所对应的背景图像；

(2)对获取的图像进行降水目标检测，获得降水目标检测结果图像。

201010145822.0公开的技术方案中已经给出了对雨滴图像进行目标检测，获得雨滴目标检测结果的方法。本实施例中，采用上述方案对待识别降水图像进行检测，得到降水目标检测结果图像。该结果图像中一般会含有多个检测出的目标，每个目标可以为雨滴或雪花。

具体过程为：将待识别降水图像及背景图像分别进行各向异性扩散平滑滤波，结果图像进行差分操作，获取差分图像。对差分图像提取边缘梯度并进行二值化，经过形态学操作处理，得到降水目标检测结果图像。

(3)特征提取步骤，即对降水目标检测结果图像中每个目标，提取其图像特征。

(3.1)对比度。

对比度定义为：将降水目标所在区域的灰度图像用基于OTSU阈值分割的方法分为两类，具有较高灰度值一类的灰度均值与较低灰度值一类的灰度均值的比值作为对比度。

图6所示的是一幅降水目标为雨滴的图像经各向异性扩散滤波后图像的灰度分布。在图6中可以看到，直方图有双峰，灰度值高的峰对应雨滴的边缘较亮的部分，灰度值低的峰对应雨滴内部较暗的部分。灰度较高一类均值和灰度较低一类均值相差很大，对比度为3.56。

图7所示的是一幅雪花图像经各向异性扩散滤波后图像的灰度分布。在图7中可以看到，与图6相比，直方图没有明显的双峰，仅有单峰，对应的是雪花区域均匀分布的灰度值。灰度较高一类均值和灰度较低一类均值相差较小，对比度为1.40。

(3.2)凹凸度。

几何形状的凹凸度与形状的凸包有关。点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。图9中由黑色线段表示的多边形就是点集Q＝{p0，p1，...p12}的凸包。可以形象地想成这样：在地上放置一些不可移动的木桩，用一根橡皮筋从最外围套住这些木桩，橡皮筋最后收缩成的形状就是凸包了。

记降水目标所在图像的连通域为S，其边缘上各点组成的点集Q的凸包为C，即C是连通域S的凸包。

连通域的面积凹凸度(Convexity_Area)定义为：

$\mathrm{CA}=\frac{\mathrm{Area}\left(S\right)}{\mathrm{Area}\left(C\right)}---\left(1\right)$

连通域的周长凹凸度(Convexity_Perimeter)定义为：

$\mathrm{CP}=\frac{\mathrm{Perimeter}\left(C\right)}{\mathrm{Perimeter}\left(S\right)}---\left(2\right)$

其中Area、Perimeter函数分别表示求连通区域的面积和周长。连通域的面积Area(S)为该连通域的像素数。连通域的周长计算方法如下：通过边缘跟踪得到连通域边缘像素点的坐标P₁(x₁，y₁)、P₂(x₂，y₂)......P_N(x_N，y_N)。则连通域的周长为：

$\mathrm{Perimeter}\left(S\right)=\underset{n=1}{\overset{N-1}{\mathrm{\Σ}}}\sqrt{{(x_n-x_{n+1})}^2+{(y_n-y_{n+1})}^2}---\left(3\right)$

连通域的凹凸度(Convexity)定义为：

C＝CA×CP    (4)

由于Area(S)≤Area(C)，Perimeter(S)≥Perimeter(C)，故CA、CP的取值范围分别为0＜CA≤1、0＜CP≤1、0≤C≤1。

下表(表1)中给出了图8中8幅图像对应的凹凸度的特征参数值。

表1图8中8幅图像的凹凸度特征参数值

图像	凹凸度
		A1	0.9890
A2	0.9853
		B1	0.9734
B2	0.9571
		B3	0.9734
C1	0.5174

C2	0.5250
		C3	0.5173

(3.3)致密度。

致密度的定义如下：

$\mathrm{Compactness}=\frac{P^2}{A}---\left(5\right)$

式中P是目标形状的周长，A是目标形状的面积。致密度即目标周长的平方与目标面积之比。致密度描述了目标的紧凑程度。目标形状越复杂，致密度越大。

(3.4)圆形度。

圆形度的定义如下：

$\mathrm{Roundness}=\frac{A}{\mathrm{\π}{\stackrel{\‾}{r}}^2}---\left(6\right)$

式中A是目标形状的面积，

是目标边界点到图像重心的平均距离。圆形度描述了目标与标准圆的类似程度。

下表(表2)中给出了图8中8幅图像对应的致密度和圆形度的特征参数值。

表2图8中8幅图像的致密度、圆形度特征参数值

图像	致密度	圆形度
			A1	13.6346	1.0220
A2	13.1577	1.0450
			B1	17.4666	0.9593
B2	16.8991	0.9793
			B3	17.4666	0.9593
C1	40.2410	1.0885
			C2	38.8028	1.1034
C3	40.2410	1.0885

分析上表(表2)数据可以看出，类圆形形状的致密度数值较小，不规则形状的致密度数值较大。两类形状的致密度区别较大，而圆形度区别较小。

(3.5)Hu矩。

数字图像f(x，y)在点(x，y)处的(p+q)阶矩定义为(p、q为非负整数)：

$m_{p,q}=\underset{x}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}{x}^p{y}^{q}f(x,y)---\left(7\right)$

点(x，y)的(p+q)阶中心矩定义为：

${\mathrm{\μ}}_{p,q}=\underset{x}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}{(x-\stackrel{\‾}{x})}^p{(y-\stackrel{\‾}{y})}^q(x,y)---\left(8\right)$

表示重心点的坐标，其中

用零阶中心矩对其余各阶中心矩进行规格化，可以得到规格化的中心矩：

${\mathrm{\η}}_{p,q}=\frac{{\mathrm{\μ}}_{p,q}}{{{\mathrm{\μ}}_{\mathrm{0,0}}}^{\mathrm{\γ}}}$ $(\mathrm{\γ}=\frac{p+q+2}{2})---\left(9\right)$

Hu矩的定义如下：

φ₁＝η_2，0+η_0，2

φ₂＝(η_2，0-η_0，2)²+4η_1，1 ²

φ₃＝(η_3，0-3η_1，2)²+(3η_2，1-η_0，3)²

φ₄＝(η_3，0+η_1，2)²+(η_2，1+η_0，3)²

φ₅＝(η_3，0-3η_1，2)(η_3，0+η_1，2)[(η_3，0+η_1，2)²-3(η_2，1+η_0，3)²](10)

+(3η_2，1-η_0，3)(η_2，1+η_0，3)[3(η_3，0+η_1，2)²-(η_2，1+η_0，3)²]

φ₆＝(η_2，0-η_0，2)[(η_3，0+η_1，2)²-(η_2，1+η_0，3)²]+4η_1，1(η_3，0+η_1，2)(η_2，1+η_0，3)

φ₇＝(3η_2，1-η_0，3)(η_3，0+η_1，2)[(η_3，0+η_1，2)²-3(η_2，1+η_0，3)²]

+(η_3，0-3η_1，2)(η_2，1+η_0，3)[3(η_3，0+η_1，2)²-(η_2，1+η_0，3)²]

下表(表3)中给出了图8中8幅图像对应的Hu矩参数值。

表3图8中8幅图像的Hu矩参数值

注：Inf表示值极大，超过10⁶。

分析上表(表3)中数据可知，对于规则形状，Hu矩可能会出现无穷大的情况。对比图像C1、C2、C3的Hu矩数据可知，Hu矩具有缩放、旋转不变性。对比图A2、B2、C2的Hu矩数据可知，不同形状的Hu矩数值有较大区别，说明Hu矩具有区分不同形状的能力。

(3.6)Chen矩。

通过对数字图像分割、边界提取及跟踪等步骤，得到目标的边界每个像素点的坐标值，定义点(x，y)处的(p+q)阶矩为：

$m_{p,q}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{x_i}^p{y_i}^q---\left(11\right)$

点(x，y)的(p+q)阶中心矩定义为：

${\mathrm{\μ}}_{p,q}=\underset{i=1}{\stackrel{N}{\mathrm{\Σ}}}\underset{j=1}{\stackrel{N}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^p{(y_i-\stackrel{\‾}{y})}^q---\left(12\right)$

表示重心点的坐标，其中

用零阶中心矩对其余各阶中心矩进行规格化，可以得到规格化的中心矩：

${\mathrm{\η}}_{p,q}=\frac{{\mathrm{\μ}}_{p,q}}{{{\mathrm{\μ}}_{\mathrm{0,0}}}^{\mathrm{\γ}}}$ $(\mathrm{\γ}=\frac{p+q+2}{2})---\left(13\right)$

Chen矩定义如下：

φ₁＝η_2，0+η_0，2

φ₂＝(η_2，0-η_0，2)²+4η_1，1 ²

φ₃＝(η_3，0-3η_1，2)²+(3η_2，1-η_0，3)²

φ₄＝(η_3，0+η_1，2)²+(η_2，1+η_0，3)²

φ₅＝(η_3，0-3η_1，2)(η_3，0+η_1，2)[(η_3，0+η_1，2)²-3(η_2，1+η_0，3)²](14)

+(3η_2，1-η_0，3)(η_2，1+η_0，3)[3(η_3，0+η_1，2)²-(η_2，1+η_0，3)²]

φ₆＝(η_2，0-η_0，2)[(η_3，0+η_1，2)²-(η_2，1+η_0，3)²]+4η_1，1(η_3，0+η_1，2)(η_2，1+η_0，3)

φ₇＝(3η_2，1-η_0，3)(η_3，0+η_1，2)[(η_3，0+η_1，2)²-3(η_2，1+η_0，3)²]

+(η_3，0-3η_1，2)(η_2，1+η_0，3)[3(η_3，0+η_1，2)²-(η_2，1+η_0，3)²]

下表(表4)中给出了图8中8幅图像对应的Chen矩参数值。

表4图8中8幅图像的Chen矩参数值

注：Inf表示值极大，超过10⁶。

分析上表(表4)中数据可知，对于规则形状，Chen矩同样可能会出现无穷大的情况。对比图C1、C2、C3的Chen矩数据可知，Chen矩抗缩放、旋转不变性较差。对比图A2、B2、C2的Chen矩数据可知，不同形状的Chen矩数值有较大区别，说明Chen矩具有一定的区分不同形状的能力。

(3.7)相对矩。

相对矩的定义为：

$R_1=\frac{\sqrt{{\mathrm{\φ}}_2}}{{\mathrm{\φ}}_1}$ $R_2=\frac{{\mathrm{\φ}}_1+\sqrt{{\mathrm{\φ}}_2}}{{\mathrm{\φ}}_1-\sqrt{{\mathrm{\φ}}_2}}$

$R_3=\frac{\sqrt{{\mathrm{\φ}}_3}}{{\mathrm{\φ}}_4}$ $R_4=\frac{\sqrt{{\mathrm{\φ}}_3}}{\sqrt[4]{\left|{\mathrm{\φ}}_5\right|}}$

$R_5=\frac{\sqrt{{\mathrm{\φ}}_4}}{\sqrt[4]{\left|{\mathrm{\φ}}_5\right|}}$ $R_6=\frac{\left|{\mathrm{\φ}}_6\right|}{{\mathrm{\φ}}_1\·{\mathrm{\φ}}_3}---\left(15\right)$

$R_7=\frac{\left|{\mathrm{\φ}}_6\right|}{{\mathrm{\φ}}_1\·\sqrt{\left|{\mathrm{\φ}}_5\right|}}$ $R_8=\frac{\left|{\mathrm{\φ}}_6\right|}{{\mathrm{\φ}}_3\·\sqrt{\left|{\mathrm{\φ}}_2\right|}}$

$R_9=\frac{\left|{\mathrm{\φ}}_6\right|}{\sqrt{{\mathrm{\φ}}_2\·\left|{\mathrm{\φ}}_5\right|}}$ $R_{10}=\frac{\left|{\mathrm{\φ}}_5\right|}{{\mathrm{\φ}}_3\·{\mathrm{\φ}}_4}$

其中φ₁～φ₆的定义如式(10)所示。

下表(表5)中给出了图8中8幅图像对应的相对矩参数值。

表5图8中8幅图像的相对矩参数值

φ注：Inf表示值极大，超过10⁶。

分析上表(表5)中数据可知，对于规则形状，相对矩可能会出现异常情况(-1)。对比图C1、C2、C3的相对矩数据可知，相对矩抗旋转能力强，抗缩放能力也较强。对比图B2、C2的相对矩数据可知，不同形状的相对矩数值有较大区别，说明相对矩具有较强的区分不同形状的能力。

(3.8)二维极坐标傅立叶描述子。

给定一幅图像I＝{f(x，y)；0≤x＜M，0≤y＜N}，将I从笛卡尔坐标空间变换到极坐标空间，得到I_p＝{f(r，θ)；0≤r＜R，0≤θ＜2π}，式中R为径向距离最大值，r＝[(x-x_c)²+(y-y_c)²]^1/2，

将极坐标空间原点定为目标形状的质心(x_c，y_c)， $x_c=\frac{1}{M}\underset{t=0}{\overset{M-1}{\mathrm{\Σ}}}x\left(t\right),$ $y_c=\frac{1}{N}\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}y\left(t\right).$

接着对I_p进行二维傅立叶变换，得到初始变换系数数列PF：

$\mathrm{PF}(\mathrm{\ρ},\mathrm{\ψ})=\underset{r}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}f(r,{\mathrm{\θ}}_i)\×\exp [-j2\mathrm{\π}(\frac{r}{R}\mathrm{\ρ}+\frac{2\mathrm{\πi}}{T}\mathrm{\ψ})]---\left(16\right)$

式中：0≤r＜R，θ_i＝i(2π/T)，0≤ρ＜R，0≤ψ＜T；R，T分别为径向频率和角频率的分辨率。

对PF进行如下处理，得到变换系数数列PF′：

PF′(ρ，ψ)＝|PF(ρ，ψ)|    (17)

最后，对PF′进行尺度归一化，取归一化之后的向量FD作为特征向量：

$\mathrm{FD}=\{\frac{\left|\mathrm{PF}\left(\mathrm{0,0}\right)\right|}{A},\frac{\left|\mathrm{PF}\left(\mathrm{0,1}\right)\right|}{\mathrm{PF}\left(\mathrm{0,0}\right)},...,\frac{\left|\mathrm{PF}(0,n)\right|}{\left|\mathrm{PF}\left(\mathrm{0,0}\right)\right|},...,\frac{\left|\mathrm{PF}(m,0)\right|}{\left|\mathrm{PF}\left(\mathrm{0,0}\right)\right|},...,\frac{\left|\mathrm{PF}(m,n)\right|}{\left|\mathrm{PF}\left(\mathrm{0,0}\right)\right|}\}---\left(18\right)$

式中：A为目标区域面积，m，n分别为径向频率分辨率和角频率分辨率。本发明中径向频率分辨率m的取值范围为[1，5]，角频率分辨率n的取值范围为[4，36]。

(4)降水目标识别步骤，采用SVM分类器进行识别。SVM分类器的参数如下：核函数选择为RBF核函数，即

σ取值为[0.1，2]，惩罚因子取值为[2，100]。先提取训练用的目标样本的图像特征值，并根据该目标样本对应的相态类别，训练SVM分类器。然后，提取待识别降水目标的图像特征值，作为训练好的SVM分类器的输入，得到分类结果，即为待识别降水目标相态的识别结果，包括雨滴或雪花。

本发明中，在步骤(4)识别出降水目标相态的识别结果后，还可以对识别出来的结果R(x，y)进行进一步的信息提取，具体包括：降水信息提取步骤，提取雨滴目标数目、雨滴目标平面平均直径、雨滴目标立体平均直径分布谱、雪花目标数目、雪花平面面积。具体提取过程包括：

(5.1)雨滴目标的数目Num_Raindrop，即R(x，y)中灰度值为255的连通域的数目；

(5.2)设雨滴目标Ri最小外接椭圆的长轴长度为MajorL_i，短轴长度为MinL_i，则该雨滴目标的平均平面直径MeanFlatDia_i定义为：

${\mathrm{MeanFlatDia}}_i=\sqrt{{\mathrm{MajorL}}_i\×{\mathrm{MinL}}_i}$

统计每个雨滴目标的最小外接椭圆的长轴长度MajorL_i和短轴长度为MinL_i，即可得出雨滴目标的平均平面直径MeanFlatDia_i，通过实验获取平面直径与立体直径的函数关系，得到雨滴目标的平均立体直径MeanSolidDia_i；

(5.3)雨滴目标平均直径分布谱DiaSpectrum描述的是雨滴平均立体直径MeanSolidDia_i分布在雨滴目标平均立体直径值域RN中各个子值域上的数目。其中，降水图像中雨滴目标平均立体直径的值域为RN＝[MinDia，MaxDia]，MinDia为雨滴目标平均立体直径的最小值，MaxDia为雨滴目标平均立体直径的最大值，将RN平均分成Num1个子值域{RN_j，j＝1，2，…Num1-2，Num1-1，Num1}，Num1为自然数。根据雨滴目标的平均立体直径MeanDiai及上述步骤得出的各子值域即可确定出雨滴目标平均立体直径分布谱DiaSpectrum。具体实施中Num1取[20，30]间的整数；

(5.4)雪花目标的数目Num_Snowflake，即R(x，y)中灰度值为128的连通域的数目；

(5.5)雪花目标平面面积Area_Snowflake_j描述的是雪花目标S_j实际所占平面面积；

(5.6)统计雨滴目标总体积和雪花目标平面总面积，结合成像载玻片的面积和降水采集时间，得出降雨量、降雪量、降水量、降雨强度和降雪强度，根据国家气象局颁布的降水强度等级划分标准得出降水类型为小雨、中雨、大雨、暴雨、大暴雨、特大暴雨、雨夹雪、小雪、中雪、大雪、暴雪、大暴雪或特大暴雪。

本发明中，在步骤(3)提取出目标的图像特征后，还可以先对这些图像特征进行优化选择的步骤，以选择出识别最有效的若干特征，再进行降水目标识别步骤，以提高识别准确率。选择出识别最有效的若干特征可以使用现有常规的特征选择算法。

Claims (5)

1.一种基于图像处理的降水粒子相态识别方法，包括如下步骤：

(1)图像获取步骤，分别从成像设备中获取待识别的降水图像及其所对应的背景图像；

(2)对获取的图像进行降水目标检测，获得降水目标检测结果图像；

(3)提取出所述降水目标检测结果图像中的每个降水目标的多个图像特征，其中，所述图像特征包括对比度、凹凸度、致密度、圆形度、Hu矩、Chen矩、相对矩和二维极坐标傅立叶描述子；

(4)降水目标识别步骤，

首先，训练SVM分类器，即先提取训练用的目标样本的各图像特征值，并根据该目标样本对应的相态类别，训练SVM分类器；

然后，将所述每个降水目标的图像特征值，作为训练好的SVM分类器的输入，得到分类结果，即为待识别降水目标相态的识别结果，该识别结果包括雨滴和雪花。

2.根据权利要求1所述的方法，其特征在于，所述步骤（2）具体为：将待识别降水图像及背景图像分别进行各向异性扩散平滑滤波，再进行差分操作，获取差分图像，然后对差分图像提取边缘梯度并进行二值化，经过形态学操作处理，即得到降水目标检测结果图像。

3.根据权利要求1或2所述的方法，其特征在于，该方法中还包括对识别结果进行信息提取的步骤，包括：雨滴或雪花目标数目、雨滴目标平面平均直径、雨滴目标立体平均直径分布谱及雪花平面面积，其中：

（5.1）识别结果R(x,y)中灰度值为255的连通域的数目即为雨滴目标的数目，灰度值为128的连通域的数目即为雪花目标的数目；

（5.2）该雨滴目标的平均平面直径MeanFlatDia_i定义为：MajorL_i为任一雨滴目标R_i最小外接椭圆的长轴长度，MinL_i为短轴长度；

（5.3）雨滴目标的平均立体直径通过平面直径与立体直径的函数关系得到；

（5.4）设雨滴目标平均立体直径的值域为RN=[MinDia,MaxDia]，将RN平均分成Num1个子值域{RN_j,j=1,2,…Num1-2,Num1-1,Num1}，雨滴目标平均立体直径分布谱DiaSpectrum即指雨滴平均立体直径MeanSolidDia_i分布在RN中各个子值域上的数目；

(5.5)雪花目标平面面积指雪花目标实际所占的平面面积。

4.根据权利要求1或2所述的方法，其特征在于，在步骤（3）所提取的各图像特征中，先进行特征选择以选择出识别最有效的特征，再进行降水目标识别步骤。

5.根据权利要求4所述的方法，其特征在于，所述进行特征选择的方法为分支定界法、模拟退火算法、遗传算法、顺序前进法或顺序后退法。

Images