張 榮 李宏宇 周 旭 李 昊 胡向峰 夏 強

1)(中國氣象局云霧物理環(huán)境重點開放實驗室， 北京 100081)

2)(中國氣象科學研究院災害天氣國家重點實驗室， 北京 100081)

3)(河南省人工影響天氣中心， 鄭州 450003)

4)(河北省人工影響天氣中心， 石家莊 050021)

5)(中國華云氣象科技集團公司， 北京 100081)

引 言

飛機攜帶探測設備入云進行觀測是目前獲取云微物理特征最直接有效的辦法。機載云粒子探測設備根據探測原理可分為3類：基于粒子擋光原理的二維光學陣列成像儀、基于粒子米散射理論的粒子譜儀和基于電荷耦合器件(charge coupled device,CCD)照相原理的粒子成像儀[1]。能夠給出云粒子圖像的二維光學陣列探測設備目前是國內最為常用的機載探測設備之一，最早由Knollenberg[2]研發(fā)，之后又有其他類似的探測設備陸續(xù)問世，如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公司探測設備以來，在飛機探測云微物理結構和人工增雨試驗等方面的研究取得了明顯發(fā)展[7-13]。隨著國家級人工影響天氣作業(yè)飛機等航空平臺的建設，越來越多的先進機載探測設備投入使用。目前，美國DMT公司生產的CIP與PIP是國內最為常用的兩款云粒子成像探測儀器。盡管基于CIP和PIP探測數據開展的研究較多，但大多是基于其配套軟件PADS(Particle Analysis and Display System)輸出的有限信息開展的。由于PADS軟件只輸出逐秒的統計結果，并不輸出每個粒子的詳細信息，無法開展數據質量控制及粒子形狀識別等相關工作，而且PADS軟件很大程度上是個“黑盒子”，使用人員對軟件內部處理算法的了解僅限于公開的說明文檔，阻礙了國內機載探測數據的深入挖掘和分析。事實上，原始機載粒子圖像探測數據中存在很多破碎粒子和虛假粒子，嚴重影響粒子探測結果的準確性。因此，開展機載云粒子圖像原始數據解析并對數據進行質量控制有著迫切的現實需求。

能直接探測云粒子的二維形狀是機載光學陣列成像探測設備相對雷達和衛(wèi)星等遙感探測設備的一大優(yōu)勢。云粒子形狀可影響其散射特性、生長率和下落末速度等，對云降水物理及地球輻射平衡有重要影響[14-16]。由于不同形狀粒子具有不同的質量-尺度(m-D)關系，因此識別云粒子形狀有助于提高云粒子水凝物含量計算的準確度。此外，識別云粒子形狀還可用于驗證和改進地基雷達[17-18]、天基衛(wèi)星[19-20]等遙感探測設備的云相態(tài)反演算法及模擬云粒子相態(tài)的參數化方案[21-23]。雖然通過人眼識別粒子形狀能達到較高的準確率，但觀測到的粒子通常數以萬計，人眼識別非常費時，且易帶來人為主觀性誤差。因此，發(fā)展機載云粒子圖像的自動識別技術十分必要。

黃敏松等[24-25]通過解析CIP和PIP粒子原始圖像數據對破碎粒子及偽粒子的識別方法進行了研究，對提升機載粒子圖像探測數據質量起重要作用。Holroyd[26]提出利用粒子圖像幾何特征對2D-C粒子圖像進行識別的技術。王磊等[27]利用Holroyd的方法對灰度2D-C探頭所測粒子形狀識別進行研究。黃敏松等[28]也利用Holroyd提出的判別方法對CIP探測的粒子形狀進行研究。近年，國內部分CIP已升級為灰度設備(即對粒子擋光幅度進行分檔)，但對灰度CIP探測數據質量控制及粒子形狀識別的研究鮮有報道。為此，中國氣象局人工影響天氣中心飛機運行團隊經過近幾年的努力，實現了對機載云粒子圖像原始探測數據(普通設備及灰度設備)的讀取和逐個粒子詳細信息的提取，并基于粒子圖像幾何特征，對云粒子圖像進行識別和分類。

1 機載云粒子圖像探測原理

目前，DMT公司生產的CIP和PIP是國內使用最為廣泛的機載云粒子圖像探測儀器。CIP和PIP屬于光學陣列成像儀器，所采用的探測原理完全相同。設備內置的激光發(fā)生器產生的光強均勻分布的平行激光束通過一組由64個光電二極管組成的線性陣列接收。如圖1所示，當粒子通過激光束時，部分激光束會被粒子遮擋，使得部分光電二極管接收到的光強降低。光強的變化以正比于粒子通過激光束速度的采集頻率記錄，設計的采集頻率為粒子運動速度與儀器探測分辨率的比值[29]，如在100 m·s-1的真空速下，CIP的采集頻率為4 MHz(CIP分辨率為25 μm)，PIP的采集頻率為1 MHz(PIP分辨率為100 μm)。只有在采集頻率正確(即真空速正確)的情況下，才能確保每個被粒子遮擋的陰影像素點準確代表一個邊長等于儀器分辨率的正方形，否則會導致粒子圖像變形(拉長或壓扁)。當光電二極管所接收光強低于正常值的50%時，則認為該光電二極管被粒子遮擋，并記錄為0，否則記錄為1。每次記錄的一組0和1的狀態(tài)被稱為一個切片，將連續(xù)的切片拼接起來即得到該粒子的圖像。

圖1 粒子通過光電二極管陣列時形成粒子圖像的原理圖

此外，DMT公司還開發(fā)了能表示粒子擋光強弱的CIP，被稱為灰度(grayscale)CIP。它將每個陰影像素點的擋光強弱進一步細分為3檔，如光強低于正常值的25%，50%和75%，顯示粒子圖像時將擋光幅度等級用不同顏色區(qū)分，因此灰度設備相比普通光學陣列設備能給出粒子的更多細節(jié)特征。

2 粒子圖像數據質量控制

2.1 破碎粒子剔除

粒子在進入二維光學陣列探頭采樣區(qū)前，經常因碰到儀器的前端而破碎，從而導致數目眾多的碎屑粒子被儀器記錄。盡管本文使用的CIP已采用由Korolev等[30]設計的抗破碎錐面尖端，與早期CIP設備相比，進入采樣區(qū)的破碎粒子已大幅減少，但依然有些破碎粒子不可避免地進入采樣區(qū)并被儀器記錄，進而對探測結果產生不可忽略的影響[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)，這些離散點可能是自身或周圍粒子破碎產生的碎屑，這些離散點恰好與當前粒子同時進入儀器采樣區(qū)而被儀器記錄。如果不將這些離散點剔除，會影響提取粒子幾何特征參數的準確性，進而影響粒子形狀分類結果。由于這些離散點被記錄在當前粒子圖像數據中且與當前粒子圖像共用同一時間，因此無法采用時間間隔閾值方法對這些離散點進行剔除。本文通過對單個粒子所有陰影像素點周圍的8個點(上、下、左、右及左上、左下、右上、右下)進行分析，根據陰影像素點周圍8個點是否都為空挑選圖像中的陰影像素連通區(qū)，并以像素點數最多的連通區(qū)作為該粒子有效圖像進行粒子形狀分類，而將像素點數較少的連通區(qū)(如圖4中紅圈所示)作為離散點剔除。

圖3 利用相鄰粒子時間間隔閾值法剔除破碎粒子示例(紅色矩形所標記粒子為破碎粒子，灰色豎線用以分割不同粒子)

圖4 粒子本身存在離散點(紅圈所示)的原始圖像(a)及剔除離散點后的圖像(b)

2.2 偽粒子剔除

這里的偽粒子指在二極管陣列方向上或飛行方向上只有1個像素的噪點及條狀粒子(圖5)。偽粒子主要由電信號的干擾或儀器鏡頭霧化造成，在實際探測過程中一般無法避免，它們的存在會造成小粒子數量的異常增高及線狀粒子數量的偏多。本文通過判斷粒子在二極管陣列方向或飛行方向上是否只有1個像素對這類粒子進行識別并剔除。

圖5 在二極管陣列方向或飛行方向上只有1個像素的噪點(a)及條狀粒子(b)

3 粒子形狀識別

本文通過提取粒子形狀幾何特征參量，基于Holroyd[26]的方法，將云粒子形狀分為8類：微小、線狀、聚合狀、霰狀、球狀、板狀、枝狀和不規(guī)則狀。下面詳細介紹分類的具體方法和步驟。

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 粒子形狀識別示例

續(xù)圖7

3.3 分類結果說明

考慮到CIP粒子圖像僅為真實粒子的二維投影，以及其較低的像素分辨率(25 μm)和有限的探測尺度范圍(約25～1550 μm)，本文對粒子形狀的分類結果可能與云微物理學定義的粒子分類結果不完全一致。以下是對本文識別粒子形狀的一些說明：

①本文所用的粒子圖像探測設備CIP的像素分辨率為25 μm，且在對粒子形狀進行分類時，將像素數少于23個的粒子歸為微小粒子而不做具體形狀分類。因此，本文識別出的球狀粒子尺度大于50 μm，主要是凍滴或者是霰(如果在暖區(qū)，則是雨滴)。

②線狀粒子主要是柱狀和針狀冰晶。

③在區(qū)分霰狀粒子和聚合體時，主要考慮兩個因素：粒子的密實程度(S是否大于0.7)和接近圓的程度(F值是否小于9)。從分類結果看，絕大部分霰狀粒子就是通常認為的霰粒子，也有少量是聚合體。

④板狀粒子主要以片狀冰晶為主。

⑤聚合狀粒子主要是聚合體。

⑥枝狀粒子主要是雪晶，也可能包含一些不同線狀粒子的聚合體。

4 應 用

2018年12月—2019年3月，中國氣象局新舟60人工增雨飛機(編號：B-3435)在河南省進行人工增雨(雪)作業(yè)及大氣探測任務，本文以其中3次飛機探測結果為例，展示圖像形狀識別技術在云物理研究方面的拓展應用。

4.1 探測個例簡介

2018年12月10日、2019年1月8日和2月26日在河南省境內發(fā)生了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%)，不規(guī)則狀粒子占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為飛機真空速。對于光學陣列類探測設備，景深隨粒子尺度的增大而增大，但受限于探測設備兩臂(發(fā)射臂和接收臂)之間的距離，其計算公式[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種分類：微小、線狀、聚合狀、霰狀、球狀、板狀、枝狀和不規(guī)則狀。該技術解決了通過人眼對粒子形狀進行識別帶來的主觀性誤差及效率低下的問題。

2) 對比分析基于不同粒子形狀的質量-尺度關系和將所有粒子視作球形液滴計算所得的水凝物含量，發(fā)現后者平均值超過前者約1個量級，說明云粒子形狀識別技術可提高冰云或冰水混合云中水凝物含量計算的準確性。

對于光學陣列類探測設備，采樣體積與粒子尺度密切相關，采用不同粒子尺度定義得到的采樣體積不同，進而會導致基于采樣體積得到的參量(如粒子數濃度、水凝物含量等)有所不同。因此，相關研究在給出需要用到采樣體積的參量時，應說明所采用粒子尺度的定義，同時應給出對于部分形狀粒子所采用的處理技術(如整體在內、中心在內等)，如果是灰度設備，還應給出所選用的擋光幅度閾值，以便對不同研究結果進行客觀合理的比對。此外，本文在基于Holroyd方法識別粒子形狀時采用的閾值根據河南省3次冬季探測個例得出的，將其應用于其他地區(qū)時可能需適當調整部分閾值以適應當地的云粒子形狀特征。

近年來，機器學習等新技術迅猛發(fā)展，在圖像識別等領域展示了相比傳統閾值方法的明顯優(yōu)越性，特別是卷積神經網絡(CNN)作為深度學習中使用的代表性算法之一，已廣泛應用于圖像識別領域。但目前將機器學習方法用于機載云粒子圖像形狀識別的工作還很少，因此，采用機器學習等新技術手段豐富云粒子形狀識別種類并提升粒子形狀識別準確率是今后的發(fā)展方向之一。