DMT机载云粒子图像形状识别及其应用

张 荣 李宏宇 周 旭 李 昊 胡向峰 夏 强

1)(中国气象局云雾物理环境重点开放实验室， 北京 100081)

2)(中国气象科学研究院灾害天气国家重点实验室， 北京 100081)

3)(河南省人工影响天气中心， 郑州 450003)

4)(河北省人工影响天气中心， 石家庄 050021)

5)(中国华云气象科技集团公司， 北京 100081)

引 言

飞机携带探测设备入云进行观测是目前获取云微物理特征最直接有效的办法。机载云粒子探测设备根据探测原理可分为3类：基于粒子挡光原理的二维光学阵列成像仪、基于粒子米散射理论的粒子谱仪和基于电荷耦合器件(charge coupled device,CCD)照相原理的粒子成像仪[1]。能够给出云粒子图像的二维光学阵列探测设备目前是国内最为常用的机载探测设备之一，最早由Knollenberg[2]研发，之后又有其他类似的探测设备陆续问世，如PMS(Particle Measuring System)公司的二维云粒子图像探头2D-C(two-dimensional cloud)和降水粒子图像探头2D-P(two-dimensional precipitation)、DMT(Droplet Measurement Technologies)公司的二维云粒子图像探头CIP(cloud imaging probe)和降水粒子图像探头PIP(precipitation imaging probe)及SPEC(Stratton Park Engineering Company)公司的二维立体成像仪2D-S(two-dimensional stereo)和广视野降水粒子成像仪HVPS(high-volume precipitation spectrometer)等。这些探测设备被广泛应用于云物理研究[3-6]，对云微物理结构和云降水过程机理的认识起重要推动作用。

我国自20世纪80年代引入PMS公司探测设备以来，在飞机探测云微物理结构和人工增雨试验等方面的研究取得了明显发展[7-13]。随着国家级人工影响天气作业飞机等航空平台的建设，越来越多的先进机载探测设备投入使用。目前，美国DMT公司生产的CIP与PIP是国内最为常用的两款云粒子成像探测仪器。尽管基于CIP和PIP探测数据开展的研究较多，但大多是基于其配套软件PADS(Particle Analysis and Display System)输出的有限信息开展的。由于PADS软件只输出逐秒的统计结果，并不输出每个粒子的详细信息，无法开展数据质量控制及粒子形状识别等相关工作，而且PADS软件很大程度上是个“黑盒子”，使用人员对软件内部处理算法的了解仅限于公开的说明文档，阻碍了国内机载探测数据的深入挖掘和分析。事实上，原始机载粒子图像探测数据中存在很多破碎粒子和虚假粒子，严重影响粒子探测结果的准确性。因此，开展机载云粒子图像原始数据解析并对数据进行质量控制有着迫切的现实需求。

能直接探测云粒子的二维形状是机载光学阵列成像探测设备相对雷达和卫星等遥感探测设备的一大优势。云粒子形状可影响其散射特性、生长率和下落末速度等，对云降水物理及地球辐射平衡有重要影响[14-16]。由于不同形状粒子具有不同的质量-尺度(m-D)关系，因此识别云粒子形状有助于提高云粒子水凝物含量计算的准确度。此外，识别云粒子形状还可用于验证和改进地基雷达[17-18]、天基卫星[19-20]等遥感探测设备的云相态反演算法及模拟云粒子相态的参数化方案[21-23]。虽然通过人眼识别粒子形状能达到较高的准确率，但观测到的粒子通常数以万计，人眼识别非常费时，且易带来人为主观性误差。因此，发展机载云粒子图像的自动识别技术十分必要。

黄敏松等[24-25]通过解析CIP和PIP粒子原始图像数据对破碎粒子及伪粒子的识别方法进行了研究，对提升机载粒子图像探测数据质量起重要作用。Holroyd[26]提出利用粒子图像几何特征对2D-C粒子图像进行识别的技术。王磊等[27]利用Holroyd的方法对灰度2D-C探头所测粒子形状识别进行研究。黄敏松等[28]也利用Holroyd提出的判别方法对CIP探测的粒子形状进行研究。近年，国内部分CIP已升级为灰度设备(即对粒子挡光幅度进行分档)，但对灰度CIP探测数据质量控制及粒子形状识别的研究鲜有报道。为此，中国气象局人工影响天气中心飞机运行团队经过近几年的努力，实现了对机载云粒子图像原始探测数据(普通设备及灰度设备)的读取和逐个粒子详细信息的提取，并基于粒子图像几何特征，对云粒子图像进行识别和分类。

1 机载云粒子图像探测原理

目前，DMT公司生产的CIP和PIP是国内使用最为广泛的机载云粒子图像探测仪器。CIP和PIP属于光学阵列成像仪器，所采用的探测原理完全相同。设备内置的激光发生器产生的光强均匀分布的平行激光束通过一组由64个光电二极管组成的线性阵列接收。如图1所示，当粒子通过激光束时，部分激光束会被粒子遮挡，使得部分光电二极管接收到的光强降低。光强的变化以正比于粒子通过激光束速度的采集频率记录，设计的采集频率为粒子运动速度与仪器探测分辨率的比值[29]，如在100 m·s-1的真空速下，CIP的采集频率为4 MHz(CIP分辨率为25 μm)，PIP的采集频率为1 MHz(PIP分辨率为100 μm)。只有在采集频率正确(即真空速正确)的情况下，才能确保每个被粒子遮挡的阴影像素点准确代表一个边长等于仪器分辨率的正方形，否则会导致粒子图像变形(拉长或压扁)。当光电二极管所接收光强低于正常值的50%时，则认为该光电二极管被粒子遮挡，并记录为0，否则记录为1。每次记录的一组0和1的状态被称为一个切片，将连续的切片拼接起来即得到该粒子的图像。

图1 粒子通过光电二极管阵列时形成粒子图像的原理图

此外，DMT公司还开发了能表示粒子挡光强弱的CIP，被称为灰度(grayscale)CIP。它将每个阴影像素点的挡光强弱进一步细分为3档，如光强低于正常值的25%，50%和75%，显示粒子图像时将挡光幅度等级用不同颜色区分，因此灰度设备相比普通光学阵列设备能给出粒子的更多细节特征。

2 粒子图像数据质量控制

2.1 破碎粒子剔除

粒子在进入二维光学阵列探头采样区前，经常因碰到仪器的前端而破碎，从而导致数目众多的碎屑粒子被仪器记录。尽管本文使用的CIP已采用由Korolev等[30]设计的抗破碎锥面尖端，与早期CIP设备相比，进入采样区的破碎粒子已大幅减少，但依然有些破碎粒子不可避免地进入采样区并被仪器记录，进而对探测结果产生不可忽略的影响[31]。因此，在对粒子形状进行分类前，首先要剔除这些破碎粒子。本文采用Field等[32-33]提出的基于相邻粒子时间间隔阈值方法对破碎粒子进行剔除，即将相邻粒子时间间隔小于某阈值的粒子作为破碎粒子进行剔除。该方法基于的基本观测事实：破碎的粒子会形成一簇在空间上分布紧密的粒子群，它们通过仪器采样体积时，相邻粒子的时间间隔相对自然粒子明显较短。本文统计所用CIP资料的相邻粒子时间分布特征，这里取其中30 s的统计结果进行展示，如图2所示。图2a为该时段相邻粒子间隔时间及间隔距离(采用真空速100 m·s-1)散点图，图2b为该时段相邻粒子间隔时间分布图。由图2b可以看到，相邻粒子时间间隔最大值位于2×10-3s附近。综合上述统计结果，将2×10-5s作为CIP的粒子破碎时间间隔阈值，即在100 m·s-1的真空速下，将间隔距离小于2 mm的粒子作为破碎粒子进行剔除。利用该方法剔除破碎粒子的示例如图3所示，其中红色矩形框内粒子因与其前一相邻粒子的时间间隔小于2×10-5s从而被剔除。

图2 CIP相邻粒子间隔时间和间隔距离(a)及相邻粒子间隔时间分布(b)

此外，由于仪器记录的数据中有时会包含离散点的情况(如图3中所示粒子P)，这些离散点可能是自身或周围粒子破碎产生的碎屑，这些离散点恰好与当前粒子同时进入仪器采样区而被仪器记录。如果不将这些离散点剔除，会影响提取粒子几何特征参数的准确性，进而影响粒子形状分类结果。由于这些离散点被记录在当前粒子图像数据中且与当前粒子图像共用同一时间，因此无法采用时间间隔阈值方法对这些离散点进行剔除。本文通过对单个粒子所有阴影像素点周围的8个点(上、下、左、右及左上、左下、右上、右下)进行分析，根据阴影像素点周围8个点是否都为空挑选图像中的阴影像素连通区，并以像素点数最多的连通区作为该粒子有效图像进行粒子形状分类，而将像素点数较少的连通区(如图4中红圈所示)作为离散点剔除。

图3 利用相邻粒子时间间隔阈值法剔除破碎粒子示例(红色矩形所标记粒子为破碎粒子，灰色竖线用以分割不同粒子)

图4 粒子本身存在离散点(红圈所示)的原始图像(a)及剔除离散点后的图像(b)

2.2 伪粒子剔除

这里的伪粒子指在二极管阵列方向上或飞行方向上只有1个像素的噪点及条状粒子(图5)。伪粒子主要由电信号的干扰或仪器镜头雾化造成，在实际探测过程中一般无法避免，它们的存在会造成小粒子数量的异常增高及线状粒子数量的偏多。本文通过判断粒子在二极管阵列方向或飞行方向上是否只有1个像素对这类粒子进行识别并剔除。

图5 在二极管阵列方向或飞行方向上只有1个像素的噪点(a)及条状粒子(b)

3 粒子形状识别

本文通过提取粒子形状几何特征参量，基于Holroyd[26]的方法，将云粒子形状分为8类：微小、线状、聚合状、霰状、球状、板状、枝状和不规则状。下面详细介绍分类的具体方法和步骤。

3.1 粒子几何参量定义

粒子形状分类所用的几何特征参量如图6所示。Dx为粒子在飞行方向(x方向)上的尺度；Dy为粒子在二极管阵列方向(y方向)上的尺度；d为粒子最大弦长，即粒子在任意方向上的最大尺度；w为垂直于最大弦长d方向的最大尺度；p为粒子图像周长，即粒子边缘像素点连线(图6绿线所示)的长度。本文所有长度单位均为像素点数量，像素点为所有挡光幅度超过30%(本文所用灰度CIP挡光幅度阈值为30%, 50%和70%)的点。

图6 粒子形状识别所用几何参量

r为x方向和y方向上像素点的线性相关系数，可衡量粒子图像的线性程度，其计算公式如下：

(1)

3.2 分类方法

当破碎粒子、单个粒子中的离散点及伪粒子被剔除后，采用顺序判断方式对不同形状粒子进行判别(如表1所示)。该判别方法遵循Holroyd[26]的原理和步骤，即只有当粒子不满足当前判别标准时，才进入下一判别标准。由于Holroyd判别标准是基于早期的2D-C设备(光学阵列由32个光电二极管构成，分辨率范围为0.025～0.040 mm)得出，因此本文在使用过程中对部分阈值指标进行适当调整(经过大量调整阈值的试验，使大多数判定的粒子形状与人眼识别结果一致)，调整前后的判别条件见表1。根据该方法对粒子形状分类的结果示例见图7。

表1 粒子形状判别流程

图7 粒子形状识别示例

续图7

3.3 分类结果说明

考虑到CIP粒子图像仅为真实粒子的二维投影，以及其较低的像素分辨率(25 μm)和有限的探测尺度范围(约25～1550 μm)，本文对粒子形状的分类结果可能与云微物理学定义的粒子分类结果不完全一致。以下是对本文识别粒子形状的一些说明：

①本文所用的粒子图像探测设备CIP的像素分辨率为25 μm，且在对粒子形状进行分类时，将像素数少于23个的粒子归为微小粒子而不做具体形状分类。因此，本文识别出的球状粒子尺度大于50 μm，主要是冻滴或者是霰(如果在暖区，则是雨滴)。

②线状粒子主要是柱状和针状冰晶。

③在区分霰状粒子和聚合体时，主要考虑两个因素：粒子的密实程度(S是否大于0.7)和接近圆的程度(F值是否小于9)。从分类结果看，绝大部分霰状粒子就是通常认为的霰粒子，也有少量是聚合体。

④板状粒子主要以片状冰晶为主。

⑤聚合状粒子主要是聚合体。

⑥枝状粒子主要是雪晶，也可能包含一些不同线状粒子的聚合体。

4 应 用

2018年12月—2019年3月，中国气象局新舟60人工增雨飞机(编号：B-3435)在河南省进行人工增雨(雪)作业及大气探测任务，本文以其中3次飞机探测结果为例，展示图像形状识别技术在云物理研究方面的拓展应用。

4.1 探测个例简介

2018年12月10日、2019年1月8日和2月26日在河南省境内发生了3次由冷锋过境引起的大范围层云降雪过程。中国气象局新舟60人工增雨飞机(编号：B-3435)对这3次过程进行了探测。该架飞机上搭载有CIP(灰度型)、PIP、气象综合探测系统(AIMMS)等多种云物理探测设备[34]。本文选取这3次探测个例中飞机在云内飞行且探测到丰富云粒子图像的时段进行分析，分析时段的个例情况见表2(本文所用时间均为北京时)。

表2 分析时段的个例情况

4.2 粒子形状统计特征

为确保统计结果的准确性，本文剔除了遮挡光电二极管阵列末端(第1个或第64个二极管)而没有被完整探测到的粒子，只统计被完整探测到的粒子。3次个例中不同形状粒子出现频率及不同形状粒子平均面积的统计结果如图8及表3所示，图8还给出不同个例中的典型粒子图像。可以看到，个例20181210中球状粒子出现频率最高，占比高达65.76%，其次为微小粒子(21.17%)和板状粒子(9.84%)。该个例中的典型粒子如图8所示，该个例中大部分粒子为球形，说明粒子主要以过冷水滴形式存在。个例20190108中，板状粒子出现频率最高，占比为46.84%，符合Bailey等[35]从实验室观测的及超过100万张云粒子成像仪CPI(cloud particle imager)直接观测的冰晶图像中总结的冰晶形状与温度的关系，即温度为-22～-8℃有利于板状冰晶生长。其次为球状粒子(22.72%)，说明存在少量过冷水滴。同时，线状粒子出现频率也较高(14.33%)，说明水汽凝华增长过程在该个例中对云粒子长大形成起重要作用。个例20190226中也是板状粒子出现频率最高，占比为43.28%，符合Bailey等[35]的结果。其次为霰状粒子(15.12%)和球状粒子(13.04%)，不规则状粒子占14.22%，聚合状粒子占2.73%，说明此次过程中除了水汽凝华增长外，凇附增长和聚并增长也对云粒子形成起重要作用。3次个例中，不同形状粒子的平均面积表现出基本一致的特点：聚合状粒子面积最大，其次为霰状粒子和枝状粒子，球状、板状和线状粒子面积相对较小，微小粒子面积最小。需要指出的是，由于本文剔除了遮挡光电二极管阵列末端而未被完整探测到的粒子，但实际上这类未被完整探测到的粒子大多数是聚合状或霰状等尺度较大的粒子，因此，聚合状和霰状粒子的实际占比和面积很可能高于本文结果。

表3 统计结果

图8 各形状粒子出现频率、平均面积及不同个例典型粒子图像

4.3 水凝物含量计算

DMT配套采集软件PADS未对粒子形状进行分类，它将所有粒子都当成球形液滴计算得到水凝物含量，其计算结果对冰云或冰水混合云显然不合理。对于冰粒子，粒子质量和尺度之间的关系依赖于粒子形状。本文采用Holroyd[26]总结的不同形状冰粒子质量的经验公式(详见文献26中表2)得到各粒子质量。

要计算粒子水凝物含量，还需求出对应时段的采样体积。对于光学阵列类探测设备，单位时间内的采样体积计算公式为

V=J×K×T。

(2)

式(2)中，V为单位时间(1 s)内的采样体积，J为仪器景深，K为有效阵列宽度，T为飞机真空速。对于光学阵列类探测设备，景深随粒子尺度的增大而增大，但受限于探测设备两臂(发射臂和接收臂)之间的距离，其计算公式[2]为

(3)

式(3)中，B为探测仪器两臂间的距离(本文所用CIP为70 mm)，C是表示粒子离设备焦平面(在仪器两臂正中间位置)距离的无量纲常数，通常取值范围在3～8之间[36]。当采用50%挡光幅度作为粒子探测阈值时，DMT采用C值为3[29]。由于本文采用挡光幅度30%作为粒子探测阈值，故C值应适当调大，参考O’Shea等[37]最新结果采用C值为4.6。D为粒子尺度(为与PADS保持一致，本文采用Dy作为粒子尺度)，λ为激光波长(658 nm)。式(2)中K被称为有效阵列宽度，是因为其随粒子尺度的变化而变化。对于有效阵列宽度的计算主要有整体在内(all-in)和中心在内(center-in)两种技术[38]。整体在内技术只考虑被完整探测的粒子而排除那些一端或两端遮挡光电二极管阵列末端的粒子，中心在内技术不仅考虑被完整探测的粒子还考虑粒子中心也被探测到的不完整粒子。尽管中心在内技术能得到较大的采样体积从而能接受更多的粒子，但它一般只用于圆形或准圆形的粒子，将其应用于非球形的冰粒子时会带来较大的不确定性[39]。由于本文所用探测数据均在零度层以上获取，非球形粒子占比很大，因此在计算有效阵列宽度K时采用整体在内技术，其计算公式为

K=δ×(N-X-1)。

(4)

式(4)中，δ为仪器分辨率(25 μm)，N为阵列中光电二极管数量(64)，X为在光电二极管阵列方向上被粒子遮挡的光电二极管数量(图6中的Dy)。由式(4)可知，采用整体在内技术得到的有效阵列宽度随粒子尺度增大而减小，原因是随粒子尺度增大，粒子遮挡光电二极管阵列末端从而被排除的概率增大。

基于以上方法计算3次个例的水凝物含量逐秒分布如图9蓝线所示，其平均值分别为0.082 g·m-3, 0.265 g·m-3和0.022 g·m-3。图9红线所示为将所有粒子当作球形液滴(即PADS输出结果所用算法)求得的水凝物含量，其平均值分别为0.577 g·m-3, 3.199 g·m-3和0.259 g·m-3。虽然二者变化趋势基本一致，但3次个例中将所有粒子当作球形液滴求得的水凝物含量平均值比按照不同形状的质量-尺度关系求得的水凝物含量平均值分别偏大7.0倍、12.1倍和11.7倍。由此证实，DMT配套软件PADS在处理CIP数据时所用算法严重高估了冰云或冰水混合云中的水凝物含量，应慎用。

图9 根据不同形状粒子的质量公式所得水凝物含量及将所有粒子视作球形液滴所得水凝物含量

5 结论与讨论

本文介绍机载云粒子图像数据的质量控制方法，在此基础上实现了云粒子形状自动分类。通过在河南省3次飞行探测个例中的灰度CIP数据，展示这些技术在云物理分析中的拓展应用，得到以下结论：

1) 通过解析云粒子图像原始数据，对其中的破碎粒子和伪粒子进行识别并剔除。在此基础上，根据Holroyd方法，实现对云粒子形状的8种分类：微小、线状、聚合状、霰状、球状、板状、枝状和不规则状。该技术解决了通过人眼对粒子形状进行识别带来的主观性误差及效率低下的问题。

2) 对比分析基于不同粒子形状的质量-尺度关系和将所有粒子视作球形液滴计算所得的水凝物含量，发现后者平均值超过前者约1个量级，说明云粒子形状识别技术可提高冰云或冰水混合云中水凝物含量计算的准确性。

对于光学阵列类探测设备，采样体积与粒子尺度密切相关，采用不同粒子尺度定义得到的采样体积不同，进而会导致基于采样体积得到的参量(如粒子数浓度、水凝物含量等)有所不同。因此，相关研究在给出需要用到采样体积的参量时，应说明所采用粒子尺度的定义，同时应给出对于部分形状粒子所采用的处理技术(如整体在内、中心在内等)，如果是灰度设备，还应给出所选用的挡光幅度阈值，以便对不同研究结果进行客观合理的比对。此外，本文在基于Holroyd方法识别粒子形状时采用的阈值根据河南省3次冬季探测个例得出的，将其应用于其他地区时可能需适当调整部分阈值以适应当地的云粒子形状特征。

近年来，机器学习等新技术迅猛发展，在图像识别等领域展示了相比传统阈值方法的明显优越性，特别是卷积神经网络(CNN)作为深度学习中使用的代表性算法之一，已广泛应用于图像识别领域。但目前将机器学习方法用于机载云粒子图像形状识别的工作还很少，因此，采用机器学习等新技术手段丰富云粒子形状识别种类并提升粒子形状识别准确率是今后的发展方向之一。