基于形態結構特征的對流單體自動分割方法研究

張軍，賀婷婷，侯謹毅，王萍

（天津大學電氣自動化與信息工程學院，天津 300072）

由于氣候變化，近年來對流風暴的頻率和強度都顯著增加[1-3].極端的對流風暴事件會造成嚴重的社會經濟損失.由于對流風暴高度動態的空間和時間過程，對流事件的研究仍然是一個具有挑戰性的問題.天氣雷達是監測強對流天氣（冰雹、大風、龍卷和暴洪）的主要工具之一，利用天氣雷達，可以更詳細地分析對流風暴的形成和運動過程.

在天氣雷達的反射率強度圖像上，對流風暴經常表現為多單體共存的結構.其中，有些單體交織在一起.雖然強對流災害可能由這些相互交織的單體群共同引發，但不同的對流單體所引發的災害類型、強度各有不同.因此，強對流天氣的自動識別及預警的前提是正確分割出每個對流單體.

自20 世紀50 年代以來，已經開發了許多基于天氣雷達的對流單體自動識別算法.由于對流單體的強度較高，在雷達反射率圖片上表現為局部極大值區域.最簡單的識別局部極大值區域的方法是基于閾值的方法，將雷達圖片上大于某個閾值的一片聯通區域識別為對流單體.目前兩種主流的對流單體識別方法，雷暴識別、跟蹤、分析和臨近預報方法[4]（Thunderstorm Identification，Tracking，Analysis and Nowcasting，TITAN）和SCIT[5]方法都是基于閾值的方法.其中TITAN 是一種固定閾值方法，而SCIT 是一種動態閾值方法.

一般來講，動態閾值方法是一種更加有效的方法，能夠有效的減少對流單體的分割錯誤.其他的一些基于動態閾值方法包括：侯正俊等[6]根據不同下墊面選擇不同閾值增強了TITAN 算法的地區適應能力；楊吉等[7]以SCIT 為基礎設計了一種基于三維雷達拼圖的單體分割算法；Crane[8]通過統計降水區域中的最大閾值獲得一個動態閾值對單體進行分割.

基于閾值的方法并不是唯一的能夠識別局部極大值區域的方法.Lakshmanan 等先后提出了基于K-均值聚類[9]和分水嶺算法[10]的單體識別算法，同樣通過檢測局部極大值區域來檢測對流單體，但是，這種方法容易受到噪音的干擾，容易產生過分割，需要對分割結果進一步合并處理.另外，Wang 等[11]提出了基于種子點生長和膨脹避讓的單體識別算法，該方法在閾值分割的基礎上，利用數學形態學方法獲得對流單體的更加完整的結構.Hou 等[12]使用區域樹結構描述方法，將天氣雷達圖片描述為一個區域樹形結構，然后通過剪枝法獲得圖片中的局部極大值區域，并以此作為對流單體的識別結果.

以上傳統的方法通過檢測天氣雷達反射率圖片中的局部極大值區域來識別對流單體.在實際應用中發現，簡單的通過檢測局部極大值區域來識別對流單體存在一些問題：首先，局部極大值區域不一定是對流單體，還可能是流場中的噪音和擾動；其次，傳統的對流單體識別算法只考慮對流單體的強度信息，而沒有考慮對流單體的形態結構信息，會得到錯誤或不完整等單體識別結果；最后，傳統的單體識別算法都是一個開環的結構，無反饋過程，換句話說，以往的單體識別算法不會檢驗分割結果是否正確，以及是否需要進一步對識別結果進行拆分等處理.

考慮到以上問題，本文設計了一種基于形態結構特征的對流單體動態識別算法.該算法在識別過程中加入了反饋環節，將單體識別過程設計為一個迭代識別過程，在每一步的迭代中，一個分類器通過分割結果的形態結構特征判斷每個識別結果是否為對流單體.若分割結果不是對流單體，那么進行進一步拆分.設計算法提取每個分割結果的形態結構特征也是本文工作的一個關鍵點.

1 數據和方法

1.1 數據和預處理

本文研究案例的雷達數據是由位于中國天津市塘沽區的S 波段多普勒雷達收集的.雷達站點位于中國北京東南100 km 處的（117°43′E，39°00′N）.雷達體掃模式為VCP 21，即每5 min 進行一次體積掃描.為了訓練和測試本文算法，我們一共收集了10個案例.這10 個案例中包含不同模式的對流風暴，每個案例的持續時間約為2 h.具體的案例相關信息見表1.其中，6 個案例作為訓練集，而剩余的4 個案例作為測試集.

表1 對流風暴案例列表Tab.1 List of convective storm cases

本方法設計用于處理0.5°仰角的雷達反射率圖像，原始的雷達數據被轉換到笛卡爾坐標系下，其中，笛卡爾網格的水平空間分辨率為0.01°0.01°（1 km 1 km）.為了濾除雷達圖片上的噪聲和地面雜波，對雷達圖片數據進行濾波操作，因此，利用中值濾波消除圖片中的噪音點.

1.2 對流單體識別方法

1.2.1 算法概述

圖1 為本文設計方法的流程圖.整體流程包括3個主要部分：初分割、迭代分割和判別器訓練.在初分割部分，首先利用區域樹方法檢測風暴圖像中的極大值區域，然后以極大值區域為核心，采用分水嶺方法分割風暴圖像；在迭代分割部分，首先提取分割結果的形態結構特征，然后通過分類器判斷分割結果是否為對流單體，若不是對流單體，則進行進一步拆分；迭代分割過程中用到的分類器通過人工標記的單體樣本圖片訓練得到.

圖1 算法流程圖Fig.1 Flow chart of the algorithm

1.2.2 對流單體的初分割

對流單體初分割方法是建立在采用區域樹結構描述對流風暴圖片的基礎上.區域樹結構描述對流風暴雷達圖片的方法見文獻[12].

圖2 給出了構造雷達反射率圖片區域樹結構的示意圖.圖2（a）為一個對流風暴的示意圖.圖2（b）中的P1到P2為不同閾值下的分割圖片，其中的聯通區域用標號標出.圖2（c）為描述對流風暴的抽象數據結構，其中每個節點對應一個區域，并且存儲了對應區域的相關信息.節點間的連線代表了區域的包含關系.

圖2 構建區域樹結構的示意圖Fig.2 Schematic diagram of building the region tree structure

當利用樹結構描述對流風暴的時候，需要引入一些樹結構相關的定義.樹中節點間的連線稱為邊，邊上強度值較低的節點稱為父節點，而強度值較高的節點稱為子節點.沒有子節點的節點稱為葉節點，存在多個子節點的節點為分支節點，沒有父節點的節點為樹結構中的根節點.

當對流風暴雷達反射率圖像被描述為區域樹結構之后，可以借助區域樹結構實現對風暴圖片濾波和檢測其中的局部極大值區域.在區域樹結構中，將面積小于5 的所有區域剔除，然后將剩下的區域疊加在一起，就實現了對風暴結構的濾波.如圖3（a）所示，中間為去除小塊區域之后的樹結構，右側為重新疊加區域得到的濾波結果.

在利用區域樹結構尋找局部極大值區域的時候，首先遍歷樹中的每一個葉節點，沿著其父節點向根部探索，當遇到一個分支節點時，探索過程結束，探索路徑上的節點集合對應著一個局部極大值區域.將探索路徑上的節點對應的區域疊加在一起，得到局部極大值圖像.圖3（b）給出了一個采用這種方式獲得局部極大值區域集合.其中間圖上的矩形框標記出所有探索得到的路徑，將這些路徑上的區域重疊在一起，得到了此案例中的局部極大值區域，見圖3（b）右側.

算法獲得的局部極大值區域與SCIT 方法獲得的對流單體檢測結果是相似的，其缺陷很明顯，當對流風暴內部有多個細小的反射率核的時候，對流單體的檢測結構不完整.為了獲得對流單體的完整結構，我們使用了基于距離變換的分水嶺分割算法[16-18]，利用距離變換方法計算出風暴區域內部所有點到局部極大值區域的最短歐式距離.如果一個點在局部極值內部，那么距離變換結果為0.風暴區域的距離變換示例如圖3（c）中部所示.對距離變換的結果采用分水嶺方法分割，得到對流單體的初分割結果，見圖3（c）右側.可見，初分割方法能夠獲得對流單體的完整結構.

圖3 基于區域樹結構的初分割方法示意圖Fig.3 Schematic diagram of the initial segmentation based on region tree structure

1.2.3 對流單體的迭代分割

對流單體的初分割過程只考慮了對流單體的強度信息，但是當對流單體發生分裂與合并的時候，容易造成錯誤的分割結果.因此需要對初分割的結果進行進一步處理，檢查每個分割結果的形態結構特征以判斷其是否為一個對流單體.當對流單體的形態結構表現為多個單體相互聯合的情況時，需要對其進一步拆分.因此，本文設計了一個迭代的再分割過程.

對流風暴的初分割結果送入迭代分割過程中，在迭代過程中的每一步，首先提取分割區域的形態結構特征，然后利用一個判別器決定該區域是否需要進一步進行分割，若該分割區域滿足對流單體的判據，則輸出該分割結果；否則，對分割結果進一步拆分，不斷重復以上迭代過程.在迭代過程中，重點是如何提取分割結果的幾何形態結構，以及如何對分割區域進行再分割.

在迭代分割過程中，計算每一個分割區域的以下3 個特征：面積衰減率、質心偏移率以及凹陷度.分割區域的特征計算需要借助區域樹結構描述方法.如圖4 所示，初分割結果相當于將原始的區域樹結構拆分為一組由單一路徑構成樹集合，每一個子樹對應一個初分割區域.對于每一個子樹，由于其只有一個單一的區域鏈，因此可以記錄為一個區域序列：T={R1，R2，…，Rn}.每個區域的面積和質心點同樣構成序列{S1，S2，…，Sn}，{x1，x2，…，xn}.下面給出三種特征的計算方法：

圖4 區域樹結構的初分割效果Fig.4 Initial segmentation effect of region tree structure

1）面積衰減率.首先計算相鄰兩個區域的面積比值：

這里i=2 ∶n.然后計算集合{ti}，i=2 ∶n 的方差dt，將dt作為面積衰減率的一個度量.面積衰減率度量了分割區域內部不同強度的面積的波動情況.對于一個對流單體，其內部各強度區域的面積變化較為平緩，而當分割區域內部包含多個單體時，其內部各強度區域的面積會存在一個跳變，導致dt較大，更趨向于被判定為需要進行再分割.

2）質心偏心率.與面積衰減率計算公式相似，首先計算相鄰兩個區域的質心偏移程度：

式中：d（xi-1，xi）為區域Ri-1和Ri的質心歐式距離；用于對距離進行歸一化.然后計算{ci}，i=2 ∶n的方差dc，將其作為質心偏心率的一個度量.通過分析，質心偏心率越大，相鄰兩反射率閾值下的連通域質心點距離越遠，分割區域內部更趨向于需要進行再分割.

3）凹陷度.具有分裂趨勢的單體通常會有凹陷的較高反射率的區域，而一個完整的對流單體則不會有凹陷較高反射率的區域，如圖5（a）所示.因此可以求取分割區域的凹陷度作為識別對流單體的一個結構特征.

圖5 給出了計算分割區域凹陷度的方法，首先計算每個區域{R1，R2，…，Rn} 的凸包多邊形的面積Sa1，Sa2，…，San，再計算每個區域面積與區域凸包多邊形的面積比值ei，然后計算不同閾值下最大的凹陷面積比em，作為凹陷度的度量.

同時記錄最大凹陷面積比所屬的強度Im.凹陷度值越大，說明分割區域的凹陷程度越厲害，則越傾向于被判定為需要進一步分割.在圖5 中，圖5（b）是分割結果，其中具有45 dBZ 的反射率區域（見圖5（c））具有最大的區域凹陷面積比，該區域的凸包多邊形如圖5（d）所示.

圖5 凹陷度計算示意圖Fig.5 Schematic diagram of concavity calculation

對于每一個分割結果，計算其3 個形態結構特征，得到一個三維特征向量，在迭代過程中，將該向量輸入到二分類SVM 中，對于被判定為需要進一步分割的多單體風暴，采用數學形態學開運算對其再分割.開運算是指對二值圖片先腐蝕后膨脹，開運算操作可以平滑對象的輪廓，消除細小的尖刺，斷開窄小的連接.詳細的再分割流程圖見圖6，該流程的說明如下：

圖6 迭代分割方法流程圖Fig.6 Flow chart of the iterative segmentation method

①利用最大凹陷面積比強度Im對原始單體區域進行二值化，得到區域核.

②對強度Im下的區域核進行一次開運算處理，對開運算處理過的區域核進行中值濾波，濾除斷開的細小的毛刺區域.

③識別所有開運算及濾波之后的聯通區域，并且根據面積進行排序.找出面積最大的聯通區域的面積，計算其他區域與最大面積的比值.

④若所有的面積比值都小于0.5，說明沒有完成區域拆分，轉步驟②.否則，說明開運算將該區域拆分為兩個及兩個以上的子區域.

⑤將面積最大的兩個連通域作為符合要求的核，基于原本單體的掩膜進行分水嶺分割操作，得到兩個分割區域.

圖7 展示了一個案例來說明如何使用數學形態學方法對單體區域進行再分割.首先，在圖7（a）中的所有反射率區域中，利用Im提取圖7（a）中的連通區域，得到如圖7（b）所示的圖像.對圖7（b）圖像進行多次的開運算和中值濾波處理，直到該區域被拆分為兩個核區，如圖7（c）所示，然后利用圖7（c）中的多個核區對原始的分割區域進行距離變換和分水嶺分割操作.得到兩個再分割的區域，如圖7（d）所示.此時，再分割過程結束.

圖7 區域再分割示例Fig.7 An example of region re-segmentation

1.2.4 分割結果判別器訓練

在迭代分割過程中，分割結果的判別和進一步分割不斷交替進行，判別過程利用了一個二分類的SVM 分類器判斷是否需要對單體進行進一步的分割.SVM 分類器可以解決小樣本下機器學習的問題，提高泛化性能[19]，考慮到訓練集樣本的數據量及三特征的空間分布，本文選擇SVM 分類器判斷分割區域是否為一個對流單體.

為了獲得SVM 模型，需要通過人工標記樣本，樣本包括需要再進一步拆分的分割結果和不需要進行拆分的對流單體，人工獲取樣本的過程相當于在圖1 所示的流程圖中用人工取代SVM 模型.在每一次分割的迭代過程中，通過人工判斷分割的結果是否為對流單體，如果是對流單體，那么標記為正樣本，否則，將其標記為負樣本.該標記過程不斷重復，直到所有的單體都被識別出來.

對所有人工標記的正負樣本，計算其3 個幾何特征（如前所述）.SVM 的訓練過程其實是為了找到正負樣本在高維特征空間的一個超平面，正樣本和負樣本分別分布在超平面的兩側.

2 實驗與結果

本文算法與兩種傳統的對流單體識別算法進行對比來驗證算法的有效性，這兩種對流單體識別算法分別是SCIT 方法和單閾值算法.SCIT 方法是一個基于動態閾值的極大值區域檢測算法，而單閾值方法是一個基于固定閾值的方法.評估系統平臺為個人PC 機，處理器為英特酷睿5，系統內存為8 Gb，評估實驗軟件及評估算法通過C++語言編寫.

評估過程包括定性評估和定量評估，在定性評估中，首先分析了本文設計的對流單體幾何特征對于分割對流單體的有效性，然后通過幾個典型的案例展示，在對流風暴演化的過程中（如生長、消失、分裂、合并等過程），本文算法如何分割不同風暴類型（如線性風暴、團狀風暴和孤立單體風暴）中的對流單體，以及本文方法與傳統方法的區別.在定量評估過程階段，通過人工評估對流單體的識別結果，給出不同單體分割方法的主觀定量評估結果.

2.1 定性評估

1）特征有效性.圖8（a）（b）（c）分別是3 個特征的可分性對比圖，左側分布圖表示無需再分割的單體特征，右側分布圖表示需要再分割的分割結果的特征.

可以看出：無需再分割的單體，其面積衰減率、質心偏心率和凹陷度集中分布在低值區域，而內含多單體區域的這三個特征較分散的分布于高值區域.因此可以說三特征描述的單核單體和多核“單體”分別來自兩個均值、方差均具顯著性差異的總體.簡單地說：三特征各自的分類能力是顯著的；假設令每個特征獨立承擔辯證上有待再做分割處理的任務，并取圖8 中低值類上限值做閾值，面積衰減率的漏分辨率僅為9.09%，質心偏心率和凹陷度的漏分割率也分別低于9.1%和21.05%.聯合建立三特征空間樣本分布如圖8（d）所示，顯然，這兩類樣本在三維空間的分布狀態為支持向量機將其變換到高維空間后獲得優質分類能力奠定了重要基礎，詳見后面的分析.

圖8 正負樣本的形態特征分布Fig.8 Distributions of morphological characteristics of positive and negative samples

2）案例分析.在本節中，采用一個線性風暴展示本文算法分割對流單體的各步中間結果，然后使用不同類型的風暴說明本文方法與傳統方法的區別.

線性對流風暴.圖9 展示了一個線性風暴在天氣雷達0.5°仰角上的觀測結果.雷達反射率圖片的大小為512×512，圖片中心為雷達所在的位置，3 個距離圈分別標記出離雷達中心50 km、150 km、230 km 的范圍，各條斜線的角度分別為30°、60°、…、360°.圖9 中虛線框中的部分存在多個聚集的單體，而且其上部的單體正在發生分裂，我們使用線性風暴的這一部分展示本文算法是如何工作的.

圖9 線狀風暴案例Fig.9 Linear storm case

圖10 展示了初分割方法各個步驟的輸出結果，以及最終的初分割結果.天氣雷達的反射率圖片首先需要從極坐標系轉換到笛卡爾坐標系，如圖10（a）所示，然后對雷達圖片進行中值濾波，去除噪點，得到圖10（b）.利用圖10（b）的區域樹結構，可以得到圖中的局部極大值區域，如圖10（c）所示，圖中陰影部分為局部極大值區域.同時，圖10（c）中還展示了30 dBZ 的風暴區域（圖中白色部分）.在風暴區域內部以局部極大值區域為核進行分水嶺分割，可以得到對流風暴中單體的初分割結果，如圖10（d）和圖10（e）所示.

從圖10（e）展示的對流單體初步分割的結果可以看出，正在發生分裂或者合并的相鄰單體，由于這些單體共用一個內核，難以通過局部極大值區域檢測分開.為了能夠區分并分割這些單體，需要結合對流單體的形態結構信息，對初分割結果進一步拆分.

在迭代分割過程中，第一步的輸入為圖10（e）中的所有初分割區域結果.首先計算圖10（e）中每一個分割結果的向量形態結構特征，將特征向量輸入到SVM 分類器中，根據SVM 分類器的判斷，圖10（e）中標號為14 號的分割結果不符合對流單體的判據，而其余的3 個分割結果都滿足對流單體的判定條件，因此需要對14 號分割結果進一步拆分.

圖10 初分割方法的輸出結果Fig.10 Results of the initial segmentation procedures

圖11 展示了對圖10（e）中14 號分割結果進一步拆分的各個中間步驟.在14 號分割結果的特征提取過程中，記錄了具有最大面積凹陷比的區域強度Im以反射率強度為Im的區域作為再分割目標區域，如圖11（a）所示.對圖11（a）中的區域不斷進行開運算和濾波處理，直到將該區域拆分為多個區域并且滿足拆分的判據，拆分結果見圖11（b）.在圖11（a）所示的區域中，以拆分結果為核，進行距離變換，得到圖11（c），對距離變換的結果圖進行分水嶺分割（圖11（d）是距離變換的局部放大示意圖）可以得到圖11（e）的再分割結果.

圖11 再分割各步驟結果Fig.11 Results of the re-segmentation procedures

接下來，對再分割得到的單體進行特征提取、SVM 分類，判斷是否需要進一步迭代分割.在此案例中，圖11（e）中的兩個分割結果的集合結構特征全部滿足對流單體的判據，因此無需進行再次分割，對流單體識別結束.最終利用本文方法得到的分割結果如圖12（a）所示.

圖12 中同時列出了利用單閾值分割算法及SCIT 算法的單體分割結果.為了詳細對比3 種算法的分割結果，不同算法的單體分割結果采用兩種方式標記出來，圖12（a）（c）（e）采用彩色掩膜的方式標記出本文算法、單閾值方法及SCIT 方法的單體分割結果；而圖12（b）（d）（f）再采用矩形框的形式標記出3 種算法的單體分割結果.

圖12 線性風暴分割結果Fig.12 Segmentation results of line-typed storms

由圖12 可見，單閾值方法將風暴上部的5 個相互黏連的、具有分裂趨勢的對流單體識別為了一個結構，單體分割錯誤非常明顯.相比之下，SCIT 方法能夠將線性風暴中相互臨近的對流單體分割出來，但是SCIT 方法將一對正在分裂的對流單體識別為了一個單體.另外，SCIT 方法識別出的單體結構明顯不夠完整，這些識別結果中不包含30 dBZ 及35 dBZ的反射率區域.在此案例中，本文算法獲得了最優的單體分割結果，能夠將所有對流單體正確分割，而且相比于SCIT 方法，本文算法能夠獲得更加完整的對流單體結構.

團狀風暴與孤立單體風暴的分割過程與線狀單體相同，采用的雷達基數據格式相同，且雷達反射率圖像與線狀單體所示的分辨率相同，因此，在本小節中，案例分析過程參照線狀單體，接下來，展示本文方法與單閾值方法及SCIT 算法在團狀風暴及孤立單體風暴情況下的單體分割結果.

圖13 和圖14 分別為3 種單體分割算法對于團狀風暴和孤立單體風暴中的單體的識別效果.圖13和圖14 采用了與圖12 相同的展示方法，使用色標和矩形框標記對流單體的分割結果.

圖13 團狀風暴分割結果Fig.13 Segmentation results of clustered storms

圖14 孤立單體風暴分割結果Fig.14 Segmentation results of single-cell storms

從圖13 中可以看出，在團狀風暴的案例中，多個單體相互臨近，但是沒有單體分裂/合并的現象出現.在這種情況下，單閾值分割方法同樣將一些相鄰的對流單體識別為一個單體.SCIT 方法與本文算法獲得了相似的識別結果，相對比本文方法，SCIT 方法的缺陷主要在于單體結構不夠完整，忽視了單體結構的完整性，丟失了部分信息.

從圖14 中可以看出，對于孤立單體風暴的案例，由于各個孤立單體風暴之間的距離較遠，3 種單體分割算法獲得了幾乎相同的結果，而且各種方法都能夠獲得對流單體的完整結構.

2.2 定量評估

1）對流單體分割臨界成功指數.為了定量評估各類對流單體分割算法的性能，需要人工檢驗各類分割算法的識別結果.首先人工標記出每個算法正確分割的單體，錯誤分割的單體，以及漏識的單體；然后統計每種單體出現的數目，根據統計結果計算出單體分割算法的量化評分值.

在人工檢驗過程中，算法正確分割的對流單體為具有完整結構的，存活時間超過兩個連續體掃的局部極值區域.一個對流單體的空間結構在演化過程中應該保持相對穩定.算法正確分割的單體個數記為X.錯誤的對流單體分割結果是指算法將由于擾動或者噪音造成的局部極大值區域識別為對流單體.所有錯誤的對流單體分割結果個數記為Y.漏識的單體是指由于面積較小，或強度過低等原因造成的算法不能夠分割的對流單體.所有漏識的對流單體個數記為Z.另外，當對流單體發生分裂或合并時，如果n 相鄰的對流單體被識別為一個對流單體，那么該對流單體識別結果被認為是錯誤的分割結果，同時認為算法漏識了n 個對流單體.由于分裂或合并造成對流單體分割錯誤則對分割的對流單體的數目進行單獨的統計，其數目分別記為Y′和Z′.

當通過人工標記的方式對算法的分割結果進行檢驗和統計之后，計算以下三項評分值：對流單體的識別率：

誤報率：

和臨界成功指數：

在定量評估過程中，以CSI 作為主要的算法量化評價指標，該值的范圍從0～1，CSI 評分值越高，算法性能越好.

2）形態學特征有效性驗證.本文提取了面積衰減率、偏心率及凹陷面積比3 個形態學特征以判斷單體是否需要進行迭代分割.為了驗證形態學特征的有效性，本文設計了關于特征有效性的測試，利用CSI 作為算法的評價指標，在所有流程不變的基礎上，使用最為基本的、在氣象領域最為常見的特征：如面積、最大強度、最小強度等.然后依次增加不同的形態特征，判斷選擇不同的特征對于算法有何種影響.

在2.1 節特征有效性分析中，分析了3 種形態特征對于正負樣本具有較好的可分性，且計算出面積衰減率、偏心率及凹陷面積比的漏分割率分別為9.09%、9.1%和21.05%.圖15 為正負樣本中氣象常用特征（面積、最大強度、最小強度）的分布情況.只有面積特征可以用于區分單體正負樣本，但是面積特征的漏分辨率達到43.40%.通過對比可知，本文提出的3 種形態特征更加適合于單體分割判別.

圖15 正負樣本的常用特征分布Fig.15 The distribution of common characteristics of positive and negative samples

此外，為了進一步驗證本文形態特征的有效性，我們首先使用氣象基本特征：面積、最大強度、最小強度對單體進行分割，得到CSI 評分值，然后，在使用基本特征訓練分類器的基礎上依次加入面積衰減率、偏心率及凹陷面積比三個形態特征，并且計算對應的CSI 評分值，結果如表2 所示.由表2 可知，當在訓練分類器的過程中加入形態結構特征的時候，分類器的CSI 評分值會得到改善，當加入3 種形態特征時，算法的CSI 評分值達到了最大（0.82）.該實驗結果顯示，本文選取的形態特征有助于得到更好的單體分割效果.

表2 特征有效性分析驗證結果Tab.2 Validation results of feature validity analysis

3）迭代分割的有效性驗證.在進行迭代分割的有效性評估之前，首先對算法的迭代次數進行了評估.由于迭代分割是在初分割的基礎上進行的，往往迭代分割的次數不會太多，大部分單體的初分割結果在經歷1～2 次迭代分割后即可滿足終止條件.因此，在單體分割的過程中設計了4 種情況的對比：僅僅包含初分割（不含迭代分割），加入1 次、2 次迭代分割和自適應迭代分割.其中，加入1 次、2 次迭代分割和自適應迭代分割3 種情況是在加入SVM 分類器的條件下進行的，即通過SVM 判斷分割結果是否需要進一步迭代分割.這4 種方法得到的單體分割評分值如表3 所示.當加入迭代分割之后，分割評分值會有明顯提升，這證明了迭代分割的有效性.

表3 迭代分割的有效性驗證結果Tab.3 Validation results of iterative segmentation

4）SVM 分類器的選擇及有效性能測試.本文中SVM 分類器的訓練用到了60 個樣本，其中正負樣本各30 例.由于我們關注的重點不是分類器的性能，而是利用一個分類器能夠終止分割過程.因此，本文中分類器本身是可以替換的，選擇哪一種分類器，相當于對于整體算法的一個調參過程.在測試中，我們嘗試了不同類型的分類器，檢測了各種分類器對本文算法性能的評分效果.

表4 為3 種分類器進行單體分割實驗結果對比，利用CSI 作為評價指標，得到3 種分類器對應的CSI 評分值分別為0.80、0.83 及0.79.可見，采用哪一種分類器對算法的分割效果影響不大，但使用SVM的分割結果略高，因此最終選取SVM 作為本文算法使用的分類器.

表4 3 種分類器進行單體分割實驗結果對比Tab.4 Comparison of experimental results of cells segmentation of three classifiers

為了測試SVM 的作用，在迭代分割中去掉SVM模型，也就是說，在每一個迭代分割過程中，所有分割結果不經過SVM 模型，直接進入下一步迭代分割，用這種分割方法得到的單體分割評分值見表5.

表5 不加SVM 控制迭代次數的實驗結果Tab.5 Experimental results without SVM to control the number of iterations

不加SVM 分類得到的一次迭代與二次迭代分割結果分別為0.76 和0.71，低于加入SVM 的評分值，這是因為原本通過SVM 分類器可能被判定為單體的初分割結果也被再次分割了，破壞了原本單體的結構，影響了分割結果.可見，在迭代分割過程中加入SVM 能夠提高最終分割結果的評分值，這驗證了SVM 的有效性.

5）單體定量分割結果.由于人工標記過程工作量較大，在實驗中只使用了6 個案例進行評估，這6個案例包含了3 種風暴類型，每個案例的時長約為2 h.表6 給出了人工檢驗3 種對流單體識別算法的統計評分結果.從表6 中可見，本文算法獲得了最高的CSI 評分值，達到了0.84，遠高于單閾值方法和SCIT 算法.具體分析各種算法的統計結果，可以看出本文算法性能獲得提升的主要原因是算法在對流風暴發生分裂或合并的時候，有效的降低了漏報和誤報的對流單體數目.

表6 單體定量分割結果Tab.6 Quantitative segmentation results of cells

6）算法運行效率.統計了3 種算法在10 個案例中的平均運行時間，本文算法、單閾值方法及SCIT方法處理一張雷達反射率圖片的平均時間分別為1.858 s、0.051 s 和0.143 s.對于3 種算法而言，本文算法的運行時間最長，這是由于本文方法是一種迭代的分割方法，而且在迭代的過程中需要大量的數學形態學運算.雖然本文方法運行時間較長，但相對于雷達的采樣周期（5～6 min），本文算法的運行時間可以忽略不計，完全能夠滿足雷達數據在線處理的要求.

3 結論

針對對流風暴內部單體分割困難的問題，本文提出了一種基于形態結構特征的對流單體自動分割方法.與傳統方法相比，本文方法的特點是：1）在對流單體的識別過程中考慮到了對流單體的形態結構信息，而不僅僅是反射率強度信息；2）對流單體的分割過程設計為一個迭代過程，用于消除錯誤的單體分割結果；3）該方法是一種基于數據驅動的方法，而傳統方法是基于規則的目標識別方法.實驗結果顯示，本文算法在對流單體聚集或發生分裂與合并的時候，能夠更加有效的分割出對流單體的完整結構，而且本文算法的量化評分值顯著高于傳統的SCIT方法和單閾值方法.