閆 軍, 王新舒, 韓旭日
(內(nèi)蒙古自治區(qū)氣象科學(xué)研究所, 內(nèi)蒙古 呼和浩特 010051)
雷暴等強(qiáng)對流天氣是世界公認(rèn)的嚴(yán)重威脅航空器飛行安全的天氣系統(tǒng),該類系統(tǒng)所產(chǎn)生的風(fēng)切變、強(qiáng)降水、冰雹等天氣現(xiàn)象是造成飛行事故的直接原因。因此,一直以來強(qiáng)對流都是航空氣象的研究重點(diǎn),同時(shí)也是難點(diǎn)。一方面,針對高頻次發(fā)生的強(qiáng)對流天氣,國內(nèi)外學(xué)者雖然已經(jīng)進(jìn)行了長期的研究,但是對其形成演變的機(jī)理認(rèn)知十分有限,系統(tǒng)的復(fù)雜性使得諸多問題依然沒有被解決;另一方面,社會及行業(yè)重點(diǎn)領(lǐng)域?qū)υ擁?xiàng)業(yè)務(wù)服務(wù)的要求越來越高,預(yù)報(bào)預(yù)警的準(zhǔn)確率、時(shí)空分辨率、時(shí)效性需求極高,航空安全在其中尤為突出。因此,如何利用新的技術(shù)方法快速、準(zhǔn)確地獲取強(qiáng)對流天氣影響范圍、精度及程度是臨近預(yù)報(bào)業(yè)務(wù)的研究熱點(diǎn),得到的研究成果在航空氣象勤務(wù)中能發(fā)揮有效作用[1]。
層狀云與對流云自動(dòng)識別算法的研究一直是天氣雷達(dá)應(yīng)用研究領(lǐng)域的熱點(diǎn),并且已經(jīng)產(chǎn)生了一批切實(shí)可行的成果。識別算法主要分為背景場閾值技術(shù)(BET)和概率匹配方法兩大類。D.D.Churchill 等人在雷達(dá)反射率數(shù)據(jù)基礎(chǔ)上使用BET 方法,該方法集中在核的識別上,為每一個(gè)已識別核施加一個(gè)固定的影響半徑,以此確定對流云區(qū)域[2]。M.Steiner 等人認(rèn)為文獻(xiàn)[2]所使用的固定影響半徑依據(jù)是不充分的,他們將影響半徑作為變量,它是背景場平均反射率的函數(shù),設(shè)定的門限閾值取決于背景場平均反射率,該算法已被TRMM(the Tropical Rainfall Measuring Mission)項(xiàng)目采納為降水分類的業(yè)務(wù)運(yùn)行算法[3]。DeMot 等人將上述方法擴(kuò)展到三維,通過使用采自雙極化雷達(dá)的垂直速度,證明了三維分析改進(jìn)了降水分類的精度[4]。Houze 等人也將他們的研究建立在三維分析的基礎(chǔ)上,但不同于其他的BET 方法,他們將重點(diǎn)放在通過窗區(qū)概率匹配方法來改進(jìn)降水估測,將層狀云和對流云不同降水增長機(jī)制的影響與雷達(dá)數(shù)據(jù)固有的限制相結(jié)合,由3 個(gè)參數(shù)代表降水分割結(jié)果:徑向反射率梯度、亮帶小數(shù)、有效性(測量云的深度)[5]。國內(nèi)從事這一領(lǐng)域研究的人員相對較少。王改利在Lakshmanan 的對雷暴識別聚類分析方法基礎(chǔ)上提出一種K-means 方法,該方法基于K-means 的分層聚類能夠?qū)崿F(xiàn)暴雨回波的多尺度識別。肖艷姣選取了組合反射率及其水平梯度、反射率因子等于35 dBZ 的回波頂高及其水平梯度、垂直累積液態(tài)水含量及其密度等6 個(gè)候選參數(shù),并統(tǒng)計(jì)這6 個(gè)候選識別參數(shù)分布的概率密度特征,利用模糊邏輯法進(jìn)行層狀云和對流云的識別[6-7]。
本研究是在充分汲取前人成果的基礎(chǔ)上,嘗試使用深度學(xué)習(xí)方法來試驗(yàn)對流云區(qū)和層狀云區(qū)的識別[8-10]。深度學(xué)習(xí)的技術(shù)理論已經(jīng)在多個(gè)領(lǐng)域展示出較高的應(yīng)用價(jià)值,并且不斷在計(jì)算機(jī)視覺、自然語言處理、自動(dòng)駕駛等多個(gè)方面取得重大突破。其算法模型的建立依賴大量的樣本數(shù)據(jù),天氣雷達(dá)具備數(shù)量龐大的高時(shí)空分辨率觀測數(shù)據(jù),這就為將深度學(xué)習(xí)理論引入天氣系統(tǒng)識別和預(yù)測提供了可能。在眾多的機(jī)器學(xué)習(xí)研究方向中,表征學(xué)習(xí)的重點(diǎn)是如何自動(dòng)找出表示數(shù)據(jù)的合適方式,以便能更好地將輸入變換為輸出;而本研究所采用的深度學(xué)習(xí)理論是具有多級表示的表征學(xué)方法,它可以逐級表示越來越抽象的概念或模型。在該模型的框架下,構(gòu)建數(shù)據(jù)集合和特征標(biāo)簽,以此來訓(xùn)練端到端的多層神經(jīng)網(wǎng)絡(luò),雷達(dá)數(shù)據(jù)中的對流云區(qū)可以逐級表征為特定的位置和角度的邊緣,即通過具備一定模型深度的多個(gè)隱藏層的神經(jīng)網(wǎng)絡(luò),逐層組合特征非線性變換,將原始特征映射到高維度的特征空間,從而語義分割出對流云區(qū)。
訓(xùn)練、驗(yàn)證和測試數(shù)據(jù)集的建立采用呼和浩特新一代天氣雷達(dá)2015—2022 年(呼和浩特C 波段HHHTRC)基數(shù)據(jù),該雷達(dá)饋源海拔高度為2 062 m,預(yù)警業(yè)務(wù)觀測使用VCP21 模式(9 層立體掃描),數(shù)據(jù)流間隔6 min,數(shù)據(jù)庫長為250 m。雷達(dá)站所處位置凈空條件良好,無明顯遮擋及干擾源,數(shù)據(jù)質(zhì)量較高。
1) 依據(jù)雷達(dá)觀測記錄挑選出較為顯著的對流性天氣過程的基數(shù)據(jù);
2) 由于雷達(dá)基數(shù)據(jù)獲取過程中使用的是徑向極坐標(biāo)采集和存儲的方式,加之地球曲面的變化,需要對數(shù)據(jù)做等距保角投影轉(zhuǎn)換,使之投影到與地球相平行的平面上[6]。因此,需要將極坐標(biāo)系數(shù)據(jù)通過坐標(biāo)變換的方式轉(zhuǎn)換成地理經(jīng)緯坐標(biāo)系三維數(shù)據(jù)場,坐標(biāo)轉(zhuǎn)換公式為:
式中:λ0、φ0為雷達(dá)天線所在位置的經(jīng)緯度;λ、φ為目標(biāo)物的經(jīng)緯度;L為目標(biāo)物到天線之間的距離;R為等效地球半徑;h為天線到達(dá)地面的高度;H為回波所在高度。
3) 原始基數(shù)據(jù)是沿雷達(dá)掃描線徑向分布的,為保持探測數(shù)據(jù)的連續(xù)性,相鄰掃描線之間和相鄰掃描層之間存在數(shù)據(jù)點(diǎn),需要填充以免造成信息的缺失。本文采用四線插值算法,該算法在確定掃描盲區(qū)的一個(gè)空間點(diǎn)物理量數(shù)值時(shí)需4 條雷達(dá)掃描線上的8 個(gè)極坐標(biāo)點(diǎn)。其算法為:將坐標(biāo)系中的空間某一數(shù)據(jù)點(diǎn)(xn,yn,zn),通過公式轉(zhuǎn)換到極坐標(biāo)系中去,設(shè)它在極坐標(biāo)系中的相應(yīng)位置為(Ri,θj,φk),對其取整數(shù),從而獲得與它靠近的8 個(gè)極坐標(biāo)數(shù)據(jù)點(diǎn)(Ri1,θj1,φk1)、(Ri2,θj1,φk1)、(Ri1,θj2,φk1)、(Ri2,θj2,φk1)、(Ri1,θj1,φk2)、(Ri2,θj1,φk2)、(Ri1,θj2,φk2)、(Ri2,θj2,φk2)。所對應(yīng)的回波數(shù)據(jù)如圖1 所示。
圖1 本研究邏輯結(jié)構(gòu)
根據(jù)(Ri,θj,φk)與相鄰8 個(gè)數(shù)據(jù)點(diǎn)的距離遠(yuǎn)近,對回波數(shù)據(jù)進(jìn)行四線性插值,得出(Ri,θj,φk)點(diǎn)上的數(shù)據(jù),這也就是(xn,yn,zn)上的數(shù)據(jù);最后,綜合雷達(dá)數(shù)據(jù)和深度學(xué)習(xí)算法的可行性,將經(jīng)上述處理后雷達(dá)數(shù)據(jù)場數(shù)據(jù)轉(zhuǎn)化為垂直最大回波強(qiáng)度顯示分析產(chǎn)品(CR)建立模型數(shù)據(jù)集,CR 產(chǎn)品利用數(shù)據(jù)場強(qiáng)度數(shù)據(jù),在以1 km×1 km(或2 km×2 km)為底面積的數(shù)據(jù)點(diǎn)垂直向上到回波頂?shù)拇怪敝w中,對所有位于該柱體中的回波強(qiáng)度數(shù)據(jù)點(diǎn)進(jìn)行比較,挑選出最大回波強(qiáng)度。再用測高公式計(jì)算最大回波強(qiáng)度的所有高度,從而得到最大回波強(qiáng)度分布圖像。這種產(chǎn)品有助于用戶快速查看最大回波強(qiáng)度的分布,對有些突發(fā)性的強(qiáng)對流天氣,其初始回波往往出現(xiàn)在中空,用戶使用這種產(chǎn)品可以有效監(jiān)測這一類的強(qiáng)對流天氣,在穩(wěn)定性降水條件下,還有助于用戶識別0 ℃層亮帶。由于在降雹區(qū)域,相應(yīng)的中空可能存在水分累積區(qū),所以CR 產(chǎn)品還可以作為監(jiān)測冰雹的發(fā)生發(fā)展的工具之一。
主要選取呼和浩特雷達(dá)站汛期(5 月—10 月)預(yù)警區(qū)域內(nèi)所發(fā)生的強(qiáng)對流天氣過程的探測基數(shù)據(jù),經(jīng)預(yù)處理后組成算法數(shù)據(jù)集,其中231 個(gè)數(shù)據(jù)樣本作為訓(xùn)練集,69 個(gè)數(shù)據(jù)樣本作為驗(yàn)證集,85 個(gè)數(shù)據(jù)樣本作為測試集,比例設(shè)置為3∶1∶1。數(shù)據(jù)集的時(shí)間跨度為2015—2022 年,其中2022 年的數(shù)據(jù)作為測試集樣本,在實(shí)驗(yàn)過程中2022 年的樣本不參與模型訓(xùn)練和調(diào)參,這有利于客觀地驗(yàn)證模型的泛化和遷移能力以及識別效果,如表1 所示。
表1 數(shù)據(jù)集詳細(xì)信息
基于雷達(dá)強(qiáng)度數(shù)據(jù)對流云區(qū)識別算法模型的建立可以充分借鑒計(jì)算機(jī)視覺中圖像語義分割的理論,后者著重像素級的圖像目標(biāo)識別,即通過算法標(biāo)注出每個(gè)像素所屬的對象類別,而對流云區(qū)的識別實(shí)質(zhì)上是通過強(qiáng)度數(shù)據(jù)點(diǎn)的識別劃分出不同目標(biāo)區(qū)域。
鑒于雷達(dá)數(shù)據(jù)量的原因,模型設(shè)計(jì)之初應(yīng)該考慮到在此種分割場景下算法要足夠輕量。一般的分割模型網(wǎng)絡(luò)結(jié)構(gòu)通常:輸出是一個(gè)N×H×W的概率圖,其中N是一個(gè)向量,向量的維數(shù)代表類別的數(shù)量,即向量的每個(gè)元素對應(yīng)一個(gè)類別,而這個(gè)元素的數(shù)值是概率值,表示對應(yīng)識別類別的概率,給出每一個(gè)數(shù)值點(diǎn)類別,所以H×W大小要與原數(shù)據(jù)場保持一致。分割網(wǎng)絡(luò)實(shí)質(zhì)上是一個(gè)沒有全連接層的全卷積網(wǎng)絡(luò),該類型網(wǎng)絡(luò)結(jié)構(gòu)可以分為Encoder 和Decoder 兩個(gè)部分。Encoder 部分特征圖的長寬會逐步縮小,需要下采樣;Decoder 部分特征圖的長寬會逐步增大,需要上采樣。
一個(gè)切實(shí)可行的分割網(wǎng)絡(luò)需要解決如下問題:
1) 通過盡可能少的Encoder 層數(shù)獲得較大的感受野。所謂感受野是指輸出的二維數(shù)據(jù)場上每個(gè)點(diǎn)的數(shù)值,是由輸入二維數(shù)據(jù)場上大小為kh×kw的區(qū)域的元素與卷積核每個(gè)元素相乘再相加得到的,所以輸入數(shù)據(jù)場區(qū)域內(nèi)每個(gè)元素?cái)?shù)值的改變都會影響輸出點(diǎn)的像素值,將這個(gè)區(qū)域叫作輸出數(shù)據(jù)場上對應(yīng)點(diǎn)的感受野。感受野的大小是評價(jià)分割網(wǎng)絡(luò)的特征表示能力的重要指標(biāo)之一。
2) 下采樣過程中雖然得到了高階特征,但其長和寬在不斷縮小,如果直接作為上采樣輸入,則會丟失信息。解決這一問題的方法就是通過跳躍連接,使得上采樣過程不僅使用高階特征,也可以使用同維度的低階特征。
3) 多尺度問題。在雷達(dá)探測資料中,對流性云區(qū)目標(biāo)有大有小,如何保證不同尺度的目標(biāo)均能被很好的處理呢?解決的方式是不同尺度的目標(biāo)交給不同感受野的卷積層來處理。
本研究借鑒分割理論在生物醫(yī)學(xué)領(lǐng)域的應(yīng)用成果,其中U-Net 網(wǎng)絡(luò)是采用編碼器(下采樣)-解碼器(上采樣)結(jié)構(gòu)和跳躍連接的一種設(shè)計(jì)方法,也是比較基礎(chǔ)的方法,后續(xù)的許多卷積神經(jīng)網(wǎng)絡(luò)仍延續(xù)了該算法的核心思想。U-Net 網(wǎng)絡(luò)如圖2 所示,其基本特點(diǎn)是使用全卷積網(wǎng)絡(luò)(FCN)。它與卷積神經(jīng)網(wǎng)絡(luò)(CNN)的區(qū)別在于將CNN 最后的全連接層(FC)替換為卷積層(Conv),因此該網(wǎng)絡(luò)是一個(gè)端到端的網(wǎng)絡(luò);網(wǎng)絡(luò)的Contracting path(如 圖2 左 側(cè))使 用Conv 和MaxPooling;網(wǎng) 絡(luò) 的Expansive path(如圖2 右側(cè))使用上采樣(upsampling)與左側(cè)的Contracting path、poolingde featuremap 相結(jié)合,然后逐層上采樣到heatmap;最后經(jīng)過兩次Conv,達(dá)到最終的heatmap,再用1×1 的Conv 做分類。
圖2 U-Net 網(wǎng)絡(luò)結(jié)構(gòu)圖
對流云和層狀云無論從外部形態(tài)還是內(nèi)部機(jī)理都存在很大的區(qū)別:對流云是由大氣層結(jié)構(gòu)不穩(wěn)定引起的氣流垂直運(yùn)動(dòng)所致,在雷達(dá)探測分析產(chǎn)品中對流性降水表現(xiàn)為較高的反射率數(shù)值,云體水平反射率梯度較大,云區(qū)形態(tài)分布面積較小,較大的垂直厚度,云頂部不平整,云頂高度梯度大;層狀云在雷達(dá)探測分析產(chǎn)品中所反映的反射率強(qiáng)度相對較弱,分布較為均勻,水平反射率梯度較小,其外部形態(tài)表現(xiàn)為云區(qū)分布面積較大,比較薄的垂直厚度, 頂部比較平整, 云頂高度梯度不大,常出現(xiàn)零度層亮帶等特征。根據(jù)上述云體特征與模型技術(shù)要求進(jìn)行算法實(shí)驗(yàn)。
整個(gè)模型網(wǎng)絡(luò)有23 層卷積層,其中左側(cè)為收縮網(wǎng)絡(luò),右側(cè)是擴(kuò)張網(wǎng)絡(luò)。收縮網(wǎng)絡(luò)由多組卷積神經(jīng)網(wǎng)絡(luò)組成,每組由兩層 3×3 unpadded 卷積層組成,每個(gè)卷積層后使用ReLU。每組之間使用size=2、stride=2 的最大池化進(jìn)行下采樣。每組第一次卷積時(shí),將特征圖通道數(shù)增加1 倍。擴(kuò)張網(wǎng)絡(luò)也是由多組卷積層組成,每組先使用2×2 的up-convolution(轉(zhuǎn)置卷積)將分辨率翻倍且通道數(shù)減半,然后在與收縮網(wǎng)絡(luò)中對應(yīng)的層,輸出特征結(jié)果concatenation。但是因?yàn)閁-Net 全部使用unpadded 卷積,對應(yīng)的收縮網(wǎng)絡(luò)的特征圖的分辨率更高。由于算法實(shí)驗(yàn)中采用了Overlap-tile 策略,可以將收縮網(wǎng)絡(luò)對應(yīng)的特征結(jié)果四周裁剪后,再拼接到擴(kuò)張網(wǎng)絡(luò)特征圖中,拼接后接兩個(gè)unpadde 的3×3 卷積層,每個(gè)卷積后接ReLU,擴(kuò)張網(wǎng)絡(luò)的最后一層是1×1 卷積層,將通道數(shù)轉(zhuǎn)換為類別數(shù)。
為了減小開銷,U-Net 網(wǎng)絡(luò)在訓(xùn)練參數(shù)上將batch size 設(shè)置為1,所使用的SGD 的動(dòng)量系數(shù)為0.99。加權(quán)損失部分最后一層使用逐強(qiáng)度點(diǎn)的softmax 和交叉熵?fù)p失,公式如下:
式中:ak(X)是二維數(shù)據(jù)場中數(shù)據(jù)點(diǎn)的位置X處通道k的激活值;K是通道總數(shù)(也就是類別數(shù));?∈{1 ,2,…,K},表示每一個(gè)數(shù)據(jù)點(diǎn)的類別;ω(X)是X數(shù)據(jù)點(diǎn)的損失權(quán)重。
實(shí)驗(yàn)過程中還涉及數(shù)據(jù)增強(qiáng),這是由于樣本數(shù)量相對較少,采用的方法是使用3×3 組網(wǎng)格上的隨機(jī)唯一向量生成平滑形變。
算法模型把對流云區(qū)識別分割后,需要將模型預(yù)測的結(jié)果與人工標(biāo)注的驗(yàn)證標(biāo)簽進(jìn)行比對,來衡量分割算法的優(yōu)劣。本研究采用該領(lǐng)域通用的評價(jià)指標(biāo):準(zhǔn)確率(Accuracy)、召回率(Recall)、Dice 系數(shù)。
準(zhǔn)確率(Accuracy)為在雷達(dá)分析產(chǎn)品強(qiáng)度數(shù)據(jù)總數(shù)中被正確分類的數(shù)據(jù)點(diǎn)數(shù),該指標(biāo)通常被用來衡量圖像分類網(wǎng)絡(luò)模型對像素分類的性能;召回率(Recall)是真實(shí)類別為對流云區(qū)的數(shù)據(jù)點(diǎn)被網(wǎng)絡(luò)模型正確分類的概率;Dice 系數(shù)通常在醫(yī)學(xué)圖像分割領(lǐng)域用作評價(jià)指標(biāo),該指標(biāo)用于對網(wǎng)絡(luò)模型預(yù)測值與真實(shí)標(biāo)注的相似度進(jìn)行評估。Dice 系數(shù)取值范圍為[0,1],其評估的分值越高,說明該模型的分割效果越好,用于本研究則表示算法識別結(jié)果與真實(shí)標(biāo)注越接近。準(zhǔn)確率(Accuracy)、召回率(Recall)、Dice 系數(shù)計(jì)算公式分別為:
式中:TP 為測試集中對流云區(qū)真實(shí)的數(shù)據(jù)點(diǎn)總數(shù)(被標(biāo)記為對流云區(qū)的數(shù)據(jù)點(diǎn),并被正確預(yù)測為對流云區(qū)的數(shù)據(jù)點(diǎn));TN 為測試集中非對流云區(qū)的數(shù)據(jù)點(diǎn)總數(shù)(被標(biāo)記為非對流云區(qū)的數(shù)據(jù)點(diǎn),并被正確預(yù)測為非對流云區(qū)的數(shù)據(jù)點(diǎn));FP 為測試集中對流云區(qū)真實(shí)的數(shù)據(jù)點(diǎn)總數(shù)(被標(biāo)記為非對流云區(qū)的數(shù)據(jù)點(diǎn),但被預(yù)測為對流云區(qū)的數(shù)據(jù)點(diǎn));FN 為測試集中對流云區(qū)真實(shí)的數(shù)據(jù)點(diǎn)總數(shù)(被標(biāo)記為對流云區(qū)的數(shù)據(jù)點(diǎn),但被預(yù)測為非對流云區(qū)的數(shù)據(jù)點(diǎn))[10-11]。
對實(shí)驗(yàn)的測試結(jié)果采用上述評價(jià)指標(biāo)進(jìn)行計(jì)算,最終統(tǒng)計(jì)評價(jià)指標(biāo),如表2 所示。由表2 可知:準(zhǔn)確率平均值為76.63%,最高值為82.02%;召回率與準(zhǔn)確率相近,分別為76.34%和81.89%。
表2 實(shí)驗(yàn)結(jié)果統(tǒng)計(jì) %
Dice 系數(shù)則反映出一個(gè)分割效果較好的數(shù)值,從圖中總體識別效果看:
1) 算法能夠準(zhǔn)確地識別并分割出對流云區(qū)域,特別是相對面積較大的區(qū)域(如圖3a)所示);
圖3 識別結(jié)果
2) 對于正在發(fā)展,強(qiáng)度不是很強(qiáng)的對流性云區(qū)域識別率不高,未能識別出整個(gè)大的區(qū)域,將整片區(qū)域分割為若干區(qū)域(如圖3b)所示);
3) 盡管大的區(qū)域捕捉較為準(zhǔn)確,但區(qū)域的邊緣識別不夠精準(zhǔn),邊緣識別結(jié)果較為粗糙;
4) 面積小的對流性云區(qū)域識別度較低,這與雷達(dá)數(shù)據(jù)本身存在噪聲點(diǎn)、雜波及地物干擾項(xiàng)有關(guān)(如圖3c)、圖3d)方框區(qū)域所示)。圖中,白色區(qū)域?yàn)樽R別出新一代天氣雷達(dá)對流云區(qū)大致位置。
本文嘗試將卷積神經(jīng)網(wǎng)絡(luò)模型用于天氣雷達(dá)探測特征天氣系統(tǒng)領(lǐng)域,解決對流云區(qū)域的識別及分割問題。在數(shù)據(jù)集的建立過程中用到了坐標(biāo)轉(zhuǎn)換、四線插值和組合反射率等算法,并建立了U-Net 網(wǎng)絡(luò)架構(gòu)的算法模型。經(jīng)過訓(xùn)練和調(diào)參的實(shí)驗(yàn)過程,從而找到識別效果最佳的模型網(wǎng)絡(luò)組成。最后利用準(zhǔn)確率、召回率和Dice系數(shù)等評價(jià)指標(biāo)建立識別分割效果的統(tǒng)計(jì)數(shù)值支撐。從研究成果看,總體達(dá)到了設(shè)計(jì)之初的實(shí)驗(yàn)效果;但從識別結(jié)果中也發(fā)現(xiàn)了諸多問題,在后續(xù)的研究中仍有較大的提升空間。