基于高分辨率粒子成像仪的过冷液滴识别

王 楠，焦瑞莉，黄敏松

(1.北京信息科技大学 信息与通信工程学院，北京 100101；2.中国科学院大气物理研究所 云降水物理与强风暴重点实验室，北京 100029)

0 引言

飞机积冰是指飞机在云中飞行时，云中所含的过冷液滴撞击到飞机机身表面冷却结晶的现象[1]。飞机积冰严重危害飞行安全，轻则影响飞机的稳定性和操纵性，重则可能导致机毁人亡的恶性事故[2]。中国民用航空局的统计数据显示，从建国以来我国境内共发生过126起因飞机积冰造成的飞行事故。飞机积冰的问题应引起高度重视[3]。

飞机积冰强度与过冷液滴的含量、浓度、大小等密切相关，因此，分析过冷液滴的微物理特性是研究飞机积冰的有效且可靠的方式[4]。直接观测云内过冷液滴的含量、浓度、大小等是非常困难的，可通过飞机搭载仪器，例如美国 Stratton Park 工程公司的高分辨率粒子成像仪(cloud particle imager,CPI)对云中过冷液滴进行测量。然而，在使用CPI拍摄云中粒子时，除了观测到过冷液滴，还会记录其他形状的冰晶粒子等，因此，在对过冷液滴的微物理特性进行分析前，必须先将过冷液滴从不同形状类别的粒子中识别出来。

McFarquhar等[5]通过粒子圆形度的计算识别过冷液滴。李玉莲等[6]提出了一种基于多普勒谱的过冷液滴识别方法，通过对全局谱的谱峰识别，分离出不同类型粒子的局部谱，将具有慢峰局部谱特征的粒子识别为过冷液滴。本文提出一种基于CPI拍摄的云中过冷液滴识别方法，通过图像预处理操作将CPI记录的原始粒子图片转化为二值化粒子轮廓图像，提取粒子的归一化傅里叶描述子特征，计算与标准过冷液滴归一化描述子之间的欧氏距离，由此判定粒子是否为过冷液滴。

1 图像预处理

本文图像预处理主要包括粒子图像分割和粒子轮廓提取两部分。预处理操作可以加强粒子图像信息，为粒子特征提取过程做准备。

1.1 粒子图像分割

使用CPI拍摄云中过冷液滴时，CPI将多个时刻记录的不同类别粒子按照拍摄时间的先后顺序依次保存在一张图片中，每张原始粒子图片包含的粒子种类、粒子数量均不相同，原始粒子图片示例如图1所示。直接从原始粒子图片识别过冷液滴会导致识别准确率不高，需要将原始粒子图片中的所有粒子单独裁剪出来，根据裁剪后得到的粒子图像数据集对过冷液滴进行识别。本文应用基于连通域标记的图像分割方法将原始粒子图片中的所有粒子进行裁剪，该方法主要包括原始粒子图片连通域标记和粒子图像分割两个步骤。

图1 原始粒子图片

1.1.1 图像连通域标记

连通域标记是图像分割的必要步骤。连通域标记算法的原理是首先对二值图像进行一次完整的扫描，经过图像扫描，把所有目标像素点进行标记，获取等价连通域标记对(简称为等价对)，并进行记录。对于等价对的产生，开始时认为是两个不同的连通区域，随着像素扫描的不断深入，以及像素扫描顺序的变化，这两个区域又被发现是连通的。因此，需要对等价对进行记录，表明它们是属于同一个连通区域的，从而方便在首次遍历结束后对等价对进行后续的更新。需要注意的是，在对图像中的像素点进行连通域判断时，如果两个邻域的像素点像素值相同，但分别属于两个不同的连通域中，则出现连通域标记冲突，需要记录下该等价对，以便后续进行连通域标记冲突处理。

本文应用连通域标记法对原始粒子图片进行图像标记的具体实现流程为：首先，对原始灰度化粒子图片应用最大类间方差法进行二值化处理；接着，对二值化粒子图片进行图像开运算操作以消除图片中无用的粒子时间信息，并对开运算操作后的粒子图片进行粒子孔洞填充以提高连通域标记的准确性；然后，应用连通域标记法将粒子图片中的每个粒子视为一个连通域进行标记。

本文分别采用了四连通域标记和八连通域标记法对粒子图片进行标记，并将标记结果进行对比。结果表明，经过孔洞填充后，应用两种方法标记的粒子结果相同，原始粒子图片中的所有粒子均被成功标记。对图1应用八连通域标记法程序运行时间是0.536 s，应用四连通域标记法程序的运行时间是0.489 s，标记结果如图2所示。因此，在结果相同的情况下，考虑到算法复杂度和运行时间问题，本文选择四连通域标记法对原始粒子图片进行标记。

图2 四连通域标记法标记后的粒子图像

1.1.2 图像分割

对原始粒子图片进行连通域标记后，计算出图片中每个连通域的外接矩形，根据外接矩形坐标对图片进行裁剪，将图片中的所有粒子单独裁剪出来并进行保存。对图2进行图像分割后获得的单个粒子图像如图3所示。

图3 图像分割后的单个粒子图像

1.2 粒子轮廓提取

对分割后得到的单个粒子图像进行Canny算子边缘检测及轮廓跟踪等程序，将二值化粒子图像转化为粒子轮廓图。粒子轮廓提取结果如图4所示。

图4 粒子轮廓提取结果

2 粒子特征提取

经过粒子图像预处理后，得到粒子的轮廓。为了实现过冷液滴的识别，必须对粒子的形状进行描述，即提取粒子形状的某些特征表示该目标。本文应用傅里叶描述子特征对过冷液滴进行识别，并与McFarquhar等[5]提出的粒子圆形度特征(下文简称为“粒子圆形度特征”)对过冷液滴的识别效果进行比较。

2.1 粒子圆形度特征

粒子圆形度的计算是最直观判别过冷液滴的方法。粒子圆形度为

(1)

式中：A为粒子的面积，即粒子区域中像素为1的个数；L为粒子的周长。

2.2 傅里叶描述子特征

傅里叶描述子是物体形状边界曲线的傅里叶变换系数，它是物体边界曲线信号的频域分析的结果，是一种描述物体形状的有效特征。基本原理是假定粒子的闭合边界曲线表示为一个总长度为N的离散的坐标序列{x(n),y(n)|n=0,1,…,N-1}，则它的复数形式表示为

z(n)=x(n)+jy(n)n=0,1,…,N-1

(2)

这样这条闭合曲线的边界就可以在一维空间上表示。由于闭合曲线的周期为N，这一边界轮廓曲线的坐标序列也是周期的，周期长度为N。

一维边界轮廓曲线坐标序列的傅里叶变换如下：

(3)

坐标序列的逆变换为

(4)

由于傅里叶系数具有能量向低频集中的特性，所以，仅根据少数的参数值就能够对粒子的边缘线进行确定，以完成粒子形状特征的识别与判断。直接得到的傅里叶系数是对位移、旋转、尺度以及起点具有依赖性的，可以对傅里叶系数进行归一化，仅利用傅里叶系数的幅度信息。令系数Z(0)为0，从而使‖Z(k)‖与旋转、平移和曲线起点的选择无关。把每一个系数的幅值‖Z(k)‖除以‖Z(1)‖,得到归一化傅里叶描述子D(k)：

(5)

D(k)是不随尺度变化而变化的，所以D(k)同时具有旋转、平移、尺度不变性，并且与曲线的起点位置的选择无关，把它作为傅里叶描述子准确描述粒子的形状。

将图4粒子轮廓坐标以复数形式表示，按照式(3)逐点进行傅里叶变换，得到复数形式的傅里叶变换系数，通过计算D(k)，得到归一化的傅里叶系数，表1是计算得到的归一化描述子特征向量值。

从表1可看出边缘描述的形状的能量大多集中在低频部分，傅里叶变换的高频分量往往很小且容易受到高频噪声的干扰，一般只使用归一化傅里叶描述子的低频分量计算粒子形状的相似差异。当k>20时，归一化系数已经接近0，所以k取20。当k=20时，对图5应用傅里叶描述子恢复的粒子形状如图5所示。可见，傅里叶描述子可以很好地描述粒子形状，对过冷液滴的识别具有重要应用价值。

表1 归一化描述子特征向量值

图5 利用傅里叶描述子恢复的粒子形状

傅立叶变换的各频率分量是相互正交的，应用欧氏距离计算归一化傅里叶描述子间的形状差异，欧氏距离为

(6)

3 过冷液滴识别实验

本文分别提取粒子的归一化傅里叶描述子特征和粒子圆形度特征进行过冷液滴识别实验，以比较两种特征对于过冷液滴识别的有效性。

3.1 测试集图像选取

从原始粒子图片经过图像预处理后得到的粒子轮廓图像集中选取200幅过冷液滴轮廓图和400幅其他形状的的粒子轮廓图(包括霰状、柱状、不规则小粒子和板状等)共600幅粒子图像作为测试集图像。

3.2 应用粒子圆形度识别过冷液滴

根据式(1)计算出600幅测试集图像的粒子圆形度值，按照降序原则对粒子圆形度值进行排序，将圆形度值最高的200幅粒子图像判别为过冷液滴图像。

应用圆形度值计算对测试集图像进行过冷液滴识别的程序运行时间为1.418 s，满足实时过冷液滴识别的要求。测试集图像的识别结果中，有13个其他形状的粒子被误判别为过冷液滴，过冷液滴识别的准确率为93.5%(识别准确率=正确识别过冷液滴数/过冷液滴总数)，被误判别的13个粒子如图6所示。表2是被误判别粒子在圆形度降序序列中对应的排列序号和圆形度值。

表2 被误判别粒子对应的排列序号和圆形度值

图6 被误判别的粒子图像

3.3 应用傅里叶描述子识别过冷液滴

3.3.1 标准过冷液滴选取

实验中选取标准过冷液滴用于计算和测试集图像归一化描述子之间的欧氏距离。首先从原始粒子图片经图像预处理后得到的图像数据集中选取10张过冷液滴图像，提取这10个过冷液滴的前20维归一化描述子并计算归一化描述子之间的欧氏距离，欧氏距离值均为0。接着，计算这10个过冷液滴的圆形度值，选择圆形度值最接近1的过冷液滴作为标准过冷液滴图像。

3.3.2 傅里叶描述子特征对过冷液滴识别

提取测试集图像的前20维归一化描述子特征，计算测试集图像与标准过冷液滴傅里叶描述子之间的欧氏距离，按照升序原则对距离值进行排序，将欧氏距离值最小的200个粒子判别为过冷液滴。

应用傅里叶描述子特征识别过冷液滴的程序运行时间为2.105 s,满足实时过冷液滴识别的要求。测试集图像的识别结果中，有8个其他形状的粒子被误判别为过冷液滴，过冷液滴识别的准确率为96%，误识别的8个粒子如图7所示。表3是被误判别粒子在欧氏距离升序序列中对应的排列序号和欧氏距离值。

图7 误识别的粒子图像

表3 被误判别粒子对应的排列序号和欧氏距离值

3.4 实验结果分析

应用粒子圆形度特征对测试集图像的过冷液滴识别准确率为93.5%，应用傅里叶描述子特征识别过冷液滴的准确率为96%，由此可证实与粒子圆形度特征相比，应用傅里叶描述子特征识别过冷液滴具有更高的准确率。在时间开销方面，应用粒子圆形度识别过冷液滴的程序运行时间是1.418 s，而应用傅里叶描述子识别的程序运行时间为2.105 s，两种方法均满足实时识别过冷液滴的要求。

4 验证实验与分析

为了进一步验证傅里叶描述子对过冷液滴的识别能力，设置验证实验比较傅里叶描述子和粒子圆形度特征对过冷液滴的识别效果。

4.1 验证实验图像数据集选取

任意选取30张原始粒子图片，将图片平均分为3组进行图像预处理操作，将预处理后得到的粒子轮廓图像集作为验证实验数据集。

4.2 验证实验阈值设定

根据过冷液滴识别实验结果设定一系列阈值对验证集粒子进行过冷液滴识别，以比较归一化描述子特征和粒子圆形度特征对过冷液滴的识别效果。

从表2可看出，应用粒子圆形度特征对过冷液滴识别得到的圆形度降序序列中，第一个被误识别粒子的排列序号为142，其他被误识别粒子普遍集中在序号170以后，当粒子圆形度阈值设定在圆形度序列序号140～170之间对应的值，能具有较好的识别效果。而从表3可看出，应用归一化描述子对过冷液滴识别得到的欧氏距离升序序列中，第一个被误识别粒子的排列序号为155，其他被误识别粒子普遍集中在序号180以后，当归一化描述子阈值设定在序号150～180之间对应的值，能具有较好的识别效果并保证识别的准确率。

所以，为了保证阈值选取的合理性并有效比较两种特征用于过冷液滴识别的识别效果，分别将粒子圆形度降序序列和欧氏距离升序序列中序号140、150、160、170、180对应的值设定为阈值，对验证实验数据集进行过冷液滴的识别，表4是序列序号及相对应的阈值。

表4 排序序列序号及对应阈值

4.3 验证实验

4.3.1 应用傅里叶描述子特征识别过冷液滴

分别提取3组验证集粒子图像的前20维傅里叶描述子特征，计算验证集粒子与标准过冷液滴傅里叶描述子之间的欧氏距离。将计算得到的欧氏距离值与设定的阈值进行比较，如果验证集粒子与标准过冷液滴之间的欧氏距离值小于阈值，将其判定为过冷液滴。表5是应用傅里叶描述子特征对3组验证集粒子进行过冷液滴识别的结果。

表5 应用傅里叶描述子特征对3组验证集粒子的识别结果 个

4.3.2 应用粒子圆形度特征识别过冷液滴

分别计算3组验证集粒子的圆形度值，将计算得到的圆形度值与设定的阈值进行比较，如果验证集粒子的圆形度值大于阈值，将其判定为过冷液滴。表6是应用粒子圆形度特征对3组验证集粒子进行过冷液滴识别的结果。

表6 应用粒子圆形度特征对3组验证集粒子的识别结果 个

4.4 实验结果分析

将3组实验结果绘成柱状图如图8所示。从图8可看出，与粒子圆形度特征相比，应用归一化描述子特征能够识别出更多的过冷液滴，并且具有更高的识别准确率。由此可证实，归一化傅里叶描述子特征对于过冷液滴的识别具有重要的应用价值。

图8 三组实验结果分析

5 结束语

本文提出了一种基于CPI拍摄的云中过冷液滴识别方法。首先，对CPI记录的原始粒子图片进行图像预处理操作，将原始粒子图片转化为二值化粒子轮廓图像。接着，提取粒子的归一化傅里叶描述子特征，计算与标准过冷液滴归一化描述子之间的欧氏距离，根据计算得到的欧氏距离值判定粒子是否为过冷液滴。本文通过实验比较了归一化描述子特征与McFarquhar等提出的粒子圆形度特征对过冷液滴的识别效果，实验结果表明，具有旋转、平移、尺度不变性并且和曲线起始位置无关的归一化傅里叶描述子对过冷液滴的识别具有更高的准确率，说明傅里叶描述子对过冷液滴的识别具有重要的应用价值。