MLP在雷達定量降水估測中的應用

卓健，廖勝石，蘇傳程，鄧悅，張小瓊

(1.廣西壯族自治區氣象信息中心，廣西南寧 530022；2.廣西壯族自治區氣候中心，廣西南寧 530022)

1 引 言

由于天氣雷達探測資料高時空分辨率的特點，對于短臨預報、中小尺度天氣監測等氣象業務具有重要意義[1-3]。 雷達定量降水估測(Quantitative Precipitation Estimation，QPE)是天氣雷達探測資料應用的重點，在國內外有眾多研究和業務產品，其分辨率已達到1 km、分鐘級。由于Z-I關系不確定、地形遮擋、距離衰減、探測高度等原因造成雷達估測降水存在明顯的系統性偏差，研究人員在研制定量降水估測產品時采用反射率垂直廓線訂正、降水類型分類、雙偏振量質量控制、雙偏振量降水估計、雨量計偏差訂正等方法提升估測精度。針對國內雷達型號多、標定不均一等具體情況[4-7]，研究人員從波束阻擋、融化層及反射率垂直廓線訂正、降水類型分類Z-I關系、利用雙偏振量估測降水等技術方法進行深入研究，取得豐碩成果[8-26]。

近年來，大數據技術與人工智能技術的飛速發展為氣象領域充分挖掘海量歷史數據價值提供了生機，逐步實現以數據驅動的思路來解決氣象領域的問題。在雷達定量測量降水問題上，Xiao等[16]采用人工神經網絡用于雷達降水估測，通過將雷達觀測插值不同的高度層，輸入不同高度層的雷達反射率來估測降水，結果優于Z-I等傳統方法；Chen等[17]將地面雷達和TRMM 雷達觀測數據相結合作為深度神經網絡的輸入，證明了該創新的降雨估測算法的廣闊前景和通用性；Liu 等[18]研究發現基于人工神經網絡的定量降水估測比用Z-I關系更精確；邵月紅等[19]利用BP 神經網絡對臨沂地區的降水進行估測，取得效果比概率匹配法更好。傅德勝等[20]基于徑向基函數神經網絡，建立雷達定量估測降水模型，將其用于地面降水估測；殷志遠等[21]采用遺傳神經網絡進行了雷達定量估測和預報降水的實驗，整體上提高了洪水預報的精度。陳訓來等[22]采用梯度提升決策樹法，結果要優于固定Z-I關系與動態Z-I關系。張晨陽等[23]將加權隨機森林和反射率垂直廓線結合，并考慮地形特征，效果有了一定提升。

雷達定量估測降水研究取得很多成果，但有些研究結果的計算效率，估測精度仍需改進，在實際業務上推廣困難。由于雷達回波和降水關系存在很強的地域性，廣西地處亞熱帶地區，地形復雜，降雨天氣變化多端，“太陽雨”、“牛背雨”、“街區雨”、“東邊日出西邊雨”的事例比比皆是，目前針對廣西區域開展基于深度學習的定量降水估測算法的研究還比較少，國內外現有的經驗不能直接在廣西應用，通過采用深度學習等人工智能技術，提升廣西雷達資料估測精度是亟需解決的技術問題。

汪瑛等[24]指出一個好的定量降水估測算法關鍵在于能否實時刻畫降水系統的時間變化，由此提出了動態分級Z-I關系法。該方法不依賴氣候統計，改進了對短時強降水、尤其是大量級的短時極強降水低估的問題，計算速度快，便于移植。卓健等[25]進一步改進動態分級Z-I關系法，提出了快速動態分級法(Fast Dynamic Categorical method，FDC)，指出以偏差為0 且方差最小為條件時分級均值是唯一解，使得求解QPE 的計算速度進一步提升。動態分級Z-I關系法和快速動態分級法強調了時間相關性高的數據是更合理的選擇，部分解決估測精度和計算速度的問題。本文在此基礎上進一步提出一個假設：在同一時次，回波塊是由多個不同降水強度特征的小區域聯合組成。根據這個臨近時次、臨近空間的降水回波性質相似的假設，我們在設計廣西區域分鐘級雷達定量降水估測產品算法時，采用AI 結合快速動態分級法(FDC)為核心來設計降水估測產品模型；與眾多前人研究不同，本文還嘗試在同一時次單部雷達使用超過一個Z-I關系，以此進行估測精度提升試驗。

2 資料和方法

2.1 數據來源和質控方案

本文使用的資料取自2019 年3 月10 日—10月10 日廣西雷達汛期運行模式覆蓋的全時段資料，包括10 部新一代多普勒雷達產品和約2 700個自動站(自動氣象站)降水量數據，所有資料均從CIMISS(全國綜合氣象信息共享系統)獲取。雷達產品使用各雷達站按業務考核要求上傳的PUP組合反射率產品，不做質控；自動站資料使用分鐘降水量和質控后的小時降水量數據，對降水數據進行簡單質控：對于某站的某小時數據，若該時次的分鐘降水量累計值不等于質控后的小時降水量則剔除該時次的分鐘和小時降水數據，若某時次的分鐘雨量有缺，但是累計值等于質控后的小時降水量，則將所缺分鐘降水量設置為0，其余分鐘降水量數據全部作為正常數據使用。分鐘雨量資料與雷達資料的匹配方式如下：用T-6 min 到T-1 min 累計雨量與T時刻雷達回波相匹配，將這一匹配得到的回波等級與降水的關系用于T-6 min 到T-1 min時間段的降水估測。

2.2 方 法

2.2.1 FDC方法及改進思路

FDC 方法將某一時次回波Z劃分為Lv1，Lv2，……，Lvn共n級，使用逐6 min 自動站降水數據與對應位置回波求Z-I關系，限制條件是逐時次分級的QPE 值應滿足偏差(Bias)為0 且方差(Var)最小，推導出同級回波的站點降水均值是該時次該等級回波Z-I關系的解。

方差用于描述離散程度，對于第n級回波，FDC方法估測降水的方差如式(1)，

圖1 FDC方法流程圖 表示各級回波的Z-I關系。

從AI 視角，FDC 方法屬于一種多分類算法，忽略求解Z-I關系過程細節，可以將圖1 簡化成感知器(圖2)，輸入的是雷達回波和自動站雨量數據，輸出的是Z-I關系。

圖2 FDC方法求Z-I關系感知器

第一次FDC 方法將回波和自動站分成了n類，每類對應一個強度范圍等級的雷達回波并包含一定數量的自動站點(圖3a)，每類回波對應一個QPE 值。要進一步提升FDC 方法的QPE 性能，應該在保持偏差為0 的前提，使方差盡量趨近0。FDC 方法按照自動站點所在位置的回波強度等級進行分類(圖3a)，在一個時次，每個強度范圍等級的雷達回波為一類，該類回波的QPE 是唯一的(圖3a 中綠色橫線)，但是站點觀測到的降水量可能與這個值并不一致，這就是FDC 方法方差產生的原因。

圖3 FDC分類示意圖 a.第一次分類；b.第二次分類。豎虛線表示回波分類間隔，圓點表示自動站降水量，橫線表示第一次分類小類回波的Z-I關系，橫虛線表示第二次分類小類回波的Z-I關系，圖b帶邊框圓點表示第二次分類后同顏色小類的強站點。

進一步細分可達到偏差為0，方差進一步縮小的目的，具體方法如下：首先將自動站降水值與該站點位置回波的QPE值進行比較，自動站降水值≥QPE 值的站點定義為強站點(紅色圓點，用AWSS表示)，其余則為弱站點(藍色圓點，用AWSW表示)；然后對第一次FDC 方法分出的所有回波和自動站類別進行類內二分類，即一個強度范圍等級的雷達回波，與該量級回波包含的強站點組成一個子類，與弱站點組成另一個子類，各子類分別重新求Z-I關系，這一方法將回波和自動站分成了2n類(圖3b)，各類分別使用類內站點均值作為新的該類回波Z-I關系，強站點的Z-I關系如圖中深紅色虛線表示，弱站點的Z-I關系如圖中淺綠色虛線表示。設分類前由自動站降水量順序序列｛H1，……，Hm，……，Hn｝組成，新類是在Hm處截斷形成的兩個新順序序列｛H1，……，Hm-1｝和｛Hm，……，Hn｝，由方差公式可知，兩個新序列的方差之和小于等于原序列，僅當一個新的序列為空集時兩個新序列的方差之和才等于原序列，這也可通過圖1直觀看出。理論上，可以根據誤差再繼續進行分類，從而不斷縮小方差，但是實際應用時過多分類容易發生過擬合現象：建模的訓練集得到一個“完美”的公式，在獨立檢驗的驗證集里表現不佳。

2.2.2 個 例

隨意抽取南寧SA 雷達一個時次(2019 年8 月1日03:48(世界時，下同))的雷達回波(圖4a)。該時次回波位置對應自動站雨量838 個樣本值空間分布（圖4b），傳統方法[26]，如固定關系法(Z=300I1.4)和最優法(Z=3I2.3)，擬合得到的Z-I曲線分別對應圖4c的灰、黃點線。若站點觀測值落在曲線上，則估測值與觀測值為0誤差，很多站點的觀測值并未落在擬合曲線上，這就是估測值與觀測值的均方根誤差來源。

圖4 2019年8月1日03:48南寧SA雷達組合反射率因子產品(a)、回波區6 min自動站降水量(b)、傳統方法數據擬合曲線(c)、兩次分類FDC估測結果(d) b中縱橫坐標為距離雷達中心經緯度度數，c中藍點表示站點觀測值(Z表示雷達回波強度，I表示降水強度)，d中紅色表示強站點，藍色表示弱站點，用圓圈、方框和三角符號表示該時次不同等級回波對應雨強的樣本數量。

使用FDC方法進行二次分類后，擬合得到強、弱站點的兩類Z-I關系(圖4d)，由于降水雨強為韋伯分布(Weibull distribution)，擬合得到的曲線更接近占多數的小量值，即曲線偏向該類數據的下方。值得注意的是這兩條線不同于傳統方法得到的指數形式曲線，并未呈現單調遞增的特點，體現了動態方法是用實時數據來擬合Z-I關系，每個時次參與擬合的數據不同，得到的Z-I關系是動態變化的。

通過2.2.1節的討論可知，經過兩次分類后，估測值與觀測值的均方根誤差縮小了。如-5～15 dBZ 對應的286 個站點觀測值為0 mm，兩次分類FDC 的估測值也為0 mm，這部分站點估測誤差為0，圖4d 中的藍實線和橫坐標基本重合；而圖4c 中可看到，固定關系法和最優法的Z-I曲線均位于橫坐標軸上方，即估測值＞0 mm，存在不同程度的估測誤差，其中最優法5～10 dBZ 的6 min 估測值偏差是最大的。可以預見，將分好的類再進一步分類，會有更多估測值落在Z-I關系曲線上，從而進一步縮小均方根誤差。將算法在實際系統中實現時，需要考慮可用站點數，避免過擬合現象發生，為研究方便，本文僅進行兩次分類。

2.2.3 人工智能技術選擇

機器學習算法可分為回歸(Regression)、分類(Classification)和聚類(Clustering)，分類或者聚類都比較適合用于區分任務，聚類和分類的差異在于是否使用標簽(label)。本研究將第二次分類結果作為同級回波分類的標簽，因此從分類算法中選擇算法。

兩次分類的FDC方法從結構上看是一個多層感知器(Multilayer Perceptron，MLP)，FDC 的MLP模型包含兩層感知器(圖5)，第一層感知器使用雷達回波與全部自動站觀測數據建立第一次Z-I關系G1，根據偏差區分出強站點和弱站點后，再次使用FDC 方法，第二層感知器將強站點和弱站點分別與雷達回波建立新的Z-I關系G2S和G2W。

圖5 MLP模型求Z-I關系 AWSS表示強站點，AWSW表示弱站點。

2.2.4 多部雷達的QPE拼圖

多部雷達組合反射率拼圖有平均值法和最大值法，試驗本研究多部雷達定量降水估測拼圖兩種方法小樣本測試后，發現最大值拼圖整體偏強，偏強幅度超過40%，平均值拼圖整體偏弱，偏弱幅度很小。分析平均值拼圖整體偏弱的原因：由于有一定比例的樣本處于無回波區域，這部分樣本的加入，導致整體評估時誤差呈負值，這是可接受的由樣本因素引起的正常偏差，因此取平均值拼圖法做為本研究的拼圖方式。

將各雷達降水估測分別作為一個分支，通過連結方式得到聯合雷達降水定量估測模型(圖6)，本模型的特點。

圖6 多部雷達QPE簡圖 輸出層的QPE表示最終計算結果。

(1) 模型的核心計算為FDC是可解釋的，所以這是一種沒有“黑盒子”的AI模型。

(2) 模型將所有雷達等權重看待，即對于任意位置，可觀測到該位置雷達回波的這部雷達對該位置的降水估測只由其本身性能決定，不需要考慮其他雷達從另外觀測點探測結果，也就規避了雷達探測性能均一性差異問題；從觀測角度，對于同一個氣象回波，不同雷達從不同方位對其進行探測，因為雷達的不一致，觀測角度不同，這就如同測量數據一樣多次測量，取其平均值是一個合理的做法。

(3) 隱藏層2 計算輸出層的QPE 見公式(2)，ki取值1 或0，由站點是否屬于對應的強或弱站點確定，n為kiG2i不為0 值的數量，由于降水的非負性，若某區域某部雷達被遮擋導致無法觀測，則該雷達這一位置的站點不進行估測也不納入最后計算，不影響QPE值。

(4) 模型可適應自動站和雷達站數量的動態變化。

2.2.5 KNN

更多的回波區域沒有自動站點，只計算1次ZI關系時一級回波只對應有一個QPE 值，使用MLP 模型對回波強度進行二次分類后，一級回波對應的Z-I關系存在雙值。以圖4 相同時次為例，這些站點強站點和弱站點空間分布具有強站點和弱站點有相對聚攏的趨勢(圖7a)，即同類強弱站點的具有水平空間歐氏距離小，不同類強弱站點水平空間歐氏距離大的特點，由此提出假設：“在同一時次，回波塊是由多個不同降水強度特征的小區域聯合組成”，這種相似性質站點空間聚攏特點是解決一級回波Z-I關系對應雙值或者多值的依據。

常見的機器學習分類算法有K 近鄰(KNearest Neighbor，KNN)，決策樹(Decision Tree)，樸 素 貝 葉 斯(Naive Bayesian，NB)，邏 輯 回 歸(Logistic Regression)，支持向量機(Support Vector Machine，SVM)，隨機森林(Random Forest，RF)等。根據臨近時次、臨近空間的降水回波性質相似假設，選擇KNN 對無站點區域回波進行分類，以此確定某一無站點回波區域是強區(類似強站點，圖7b 紅色區域)或者弱區(類似弱站點，圖7b 藍色區域)，對應該區域QPE值應該取值G2S還是G2W。

使用K=1 選擇2019 年8 月整月數據進行小樣本實驗，發現改進后的QPE效果存在明顯改進，對降水細節刻畫得更清晰。FDC方法QPE結果如圖7c 所示，改進后的QPE 結果如圖7d 所示，可看出兩方法的整體形狀基本一致，但是QPE 結果比圖8a更能體現較強降水區的細節，變化的區域很多，例如圖7d的淺藍色區域遠大于圖7c的，而圖7d深藍區域則小于圖7c的。用箭頭簡單標示出一些明顯改動的區域，紅色箭頭標示雷達回波較強，在第一次QPE時估測值較高，通過第二次分類，QPE降低，紫色箭頭則相反，第二次分類后QPE 提升，這些變化都更符合自動站實際觀測值，并且部分虛假回波造成的錯誤降水被正確訂正為無降水(圖7c、圖7d黑色空心箭頭)。

圖8 2019年3—10月估測值與觀測值逐站累計值散點圖 a.Train組(2 634站)；b.Test組(91站)。

3 統計檢驗方法及結果分析

3.1 檢驗方法評估指標

評估指標為常用的相關系數(r)、偏差(ME)、均方根誤差(RMSE)，

式中Ra(i)為第i點上的估測值，Rg(i)為第i點上的觀測值，Rˉa、Rˉg分別為估測值和觀測值的均值。

3.2 個例檢驗

前文從理論角度探討了本文研究方法QPE估測精度相對之前方法的改進思路，表1 給出圖4 例子的統計指標，比較傳統固定關系法[26]、最優法[26]、FDC方法和本文介紹的兩次分類FDC方法在單站雷達定量降水估測時的估測精度。

表1 幾種單雷達QPE方法，單時次統計指標

可看出兩次分類FDC得到的幾個統計指標在四種方法比較中都是最優。FDC 方法強調整體無偏，所以不管多少次分類，在樣本沒有錯誤和缺漏的情況下，ME 都接近0，更多次分類的RMSE 會減小。

3.3 獨立檢驗

對本文介紹的MLP 模型+KNN 方法，生成組網雷達定量降水估測的精度進行檢驗，檢驗時間段為2019 年3 月10 日—10 月10 日。獨立樣本檢驗方法更能說明產品性能指標，具體方法為：將質控后的自動站小時降水量作為“真值”數據(簡稱為觀測值，下同)；把自動站資料分為訓練集(Train組)和驗證集(Test 組)兩組，分組原則等采用潘旸等[27]的方法，將所有區域氣象站作為Train 組(2 634 站)，所有國家氣象站作為Test 組(91 站)，Test組分鐘數據不參與QPE研制。

首先按本研究的方法使用Train 組自動站降水分鐘數據和全時段雷達產品制作出逐6 min 雷達定量降水估測產品，再將每一時次的10 個逐6 min 產品累加得到小時產品；采用K 近鄰算法，提取Train 和Test 兩組自動站點(共2 725 站)在小時產品里對應格點的估測值(簡稱為估測值，下同)，使用估測值和觀測值進行檢驗評估。表2 給出估測產品分組的統計指標，整體看各項指標Train 組均高于Test組。

表2 Train/Test組統計指標

表3 2019年Train/Test組r站數分布

對相關系數按等級分別進行劃分(表2)。Train組98.75%的站點r＞0.9，該組僅有1站相關系數小于0.5，分析該站逐時降水概率僅為1.7%，觀測累計降水量僅為23.8 mm，與同期同組站點平均降水概率11.2%，站點累計降水量平均值1 065.2 mm差距極大，分析該站因各種因素導致長時間雨量數據為0，同樣情況出現在r∈(0.6，0.7]區間的1站；Test 組r主要分布在(0.8，0.9]段，沒有r＜0.5 的站點，作為獨立檢驗組，相關系數屬于較高水準，組成員均為國家站，與一般認為國家級站點可靠性強于區域站的觀點吻合。站點累計估測值和觀測值散點圖有較好聚攏性(圖8)，特別是Train組。

對同時次所有站點估測值與觀測值分別進行累加，Train 組各時次的估測合計值與觀測合計值基本相等(圖9a)，表現在散點圖上就是基本聚攏為一條直線，體現本文方法中作為Z-I關系的FDC方法強調整體無偏的特點，Test 組則表現稍差(圖9b)，說明存在一定量的空估和漏估，這是需要進一步改進的地方。對Test 組進行不同閾值的降水分級檢驗(表4)，發現除了在[0.1，2)級別，ME 為正值，其余組ME 均為負值，說明本方法在低級別高估多于低估，高級別低估多于高估，降雨量級越大，低估越明顯，對于小于10 mm/h 量級降水，估測值和觀測值還是比較接近的。說明本方法在小量級降水的估測精度比較高，對于強降水則有一定的低估偏差。

表4 Test組分級檢驗統計指標(分級閾值)

圖9 2019年3—10月逐時估測值與觀測值合計值散點圖 a.Train組(2 634站)；b.Test組(91站)。

3.4 獨立檢驗組相關系數差異最大兩個站點簡單分析

在Test組中，相關系數最差的站(全州)相關系數r=0.539 1，估測值和觀測值散點分布如圖10a所示，關系數最好的站(德保)相關系數r=0.953 4，散點分布如圖10b 所示。德保站的散點更聚攏在估測值=觀測值的紅線附近，全州站的散點則呈離散分布狀態。初步分析這兩個站估測性能差異明顯的原因，全州站位于廣西東北角，目前使用的十部組網雷達僅桂林雷達可探測到該區域，全州站至桂林雷達直線距離約110 km，該區域山體較高，遮擋影響大，沒有第二部雷達參與估測，從而導致估測精度不足，短時強降水(超20 mm/h)的估測準確率，TS=0%；德保站位于廣西西南角，為百色、崇左兩部SA 波段雷達探測覆蓋區域，短時強降水估測準確率，TS=100%。估測精度不足的部分區域，可通過引進周邊省份雷達觀測數據，如永州雷達至全州站直線距離約80 km 且遮擋更少，或者期待廣西用于填補雷達弱區、盲區的新增雷達建設成果。

4 小 結

本文提出了一種使用MLP 為框架的QPE 方法，在隱藏層各雷達分別對探測區域進行估測，各時刻同等級回波使用超過一個Z-I關系，從而達到在輸出層得到精度較高的估測結果，算法關鍵點如下。

(1) 以FDC方法為核心，根據單站雷達各級回波的降水估測結果誤差將回波區內的站點分為強站點和弱站點兩類，比只使用一次FDC 方法對雷達降水回波的估測精度更優，分好的類可進一步繼續分類。

(2) 使用MLP 技術，將多部雷達的FDC 兩次分類方法組成一個深度學習QPE框架。隱藏層到輸出層，神經元的值為0 時，參與輸出計算的權重為0，有效避免被遮擋區域估測值為0，導致拉低QPE值。

(3) 提出一個假設：“同一時次時，臨近空間的降水回波性質相似”。依據這一假設，采用KNN方法將無自動站點區域劃歸強區或弱區，確定任意位置回波對應的Z-I關系。

本文的MLP 模型和FDC 方法都是可解釋的，是一個沒有“黑盒子”的深度學習模型，模型能動態適應雷達和自動站的動態增減。對基于該模型的2019 年3—10 月的廣西新一代天氣雷達定量估測降水產品進行檢驗分析，Train 組相關系數達0.973 7，Test組相關系數達0.825 6，結果表明該方法有較好的雷達定量估測降水能力。

根據本文方法研制的廣西地面1 km 分鐘級實況降水分析產品已于2022 年1 月起通過天擎向廣西區內用戶提供，并在廣西氣象業務內網展示。