何佳 惠建忠* 何險峰 王曙東 高金兵

(1 中國氣象局公共氣象服務(wù)中心, 北京 100081; 2 華風(fēng)氣象傳媒集團有限責(zé)任公司, 北京 100081)

引言

氣象服務(wù)的核心競爭力是不斷提升的預(yù)報質(zhì)量。通過檢驗評估給出客觀結(jié)果，一方面便于管理人員進行考核，另一方面便于研發(fā)人員對比不同方法帶來的預(yù)報性能改變，促進產(chǎn)品不斷改進完善。

國家氣象中心建設(shè)開發(fā)的基于Web的國家級天氣預(yù)報檢驗分析系統(tǒng)[1]，實現(xiàn)了涵蓋國家級省級智能網(wǎng)格預(yù)報、全國城鎮(zhèn)天氣預(yù)報等數(shù)十項檢驗業(yè)務(wù)產(chǎn)品的檢驗評估。用戶可以通過瀏覽器查看相應(yīng)產(chǎn)品的檢驗反饋信息，為各省以及國家級預(yù)報業(yè)務(wù)考核提供了信息支撐。萬夫敬[2]等分析了ECMWF模式氣溫預(yù)報在青島地區(qū)的預(yù)報性能，根據(jù)模式誤差特點，給出主觀訂正參考值，提高了氣溫預(yù)報正確率。劉麗敏[3]等通過對中央臺精細化溫度預(yù)報在伊春市的檢驗和誤差分析，掌握其預(yù)報性能，在實際工作中有針對性加以訂正，為本地釋用提供依據(jù)。趙平偉[4]等開展了GPM IMAGE、ERA5降水?dāng)?shù)據(jù)在云南的適用性評估和對比。蔡曉杰[5]等利用上海沿岸海域19個站點風(fēng)場資料使用平均誤差、預(yù)報準確率等指標對風(fēng)速、風(fēng)向、大風(fēng)過程進行檢驗。吳瑞姣[6]等運用兩分類TS評分、鄰域法FSS評分方法對安徽省WRF模式短時強降水開展檢驗。鄭琳琳[7]等采用均方誤差、TS評分等方法對改進的降水臨近外推預(yù)報技術(shù)方法進行效果檢驗，得出新方法更具優(yōu)勢的結(jié)論。劉湊華[8]等研究了基于目標的降水檢驗方法并對其改進后進行業(yè)務(wù)應(yīng)用。可見檢驗在業(yè)務(wù)中發(fā)揮著重要作用。目前檢驗針對的產(chǎn)品在時次、時效上主要是天或小時量級，少有針對分鐘級這樣高頻次發(fā)布，高時間分辨率產(chǎn)品的檢驗。產(chǎn)品在時間維度精細化的提升意味著檢驗計算時間復(fù)雜度的增加。

Python是近幾年越來越受到開發(fā)人員關(guān)注和使用的編程語言。其簡潔、易于理解、功能豐富的特點無論是專業(yè)編程人員還是普通用戶都能夠快速上手。在氣象領(lǐng)域，推出了很多易用工具包。如標簽化多維數(shù)組和數(shù)據(jù)集xarray[9]——在吸收用于實現(xiàn)科學(xué)數(shù)據(jù)自我描述的通用數(shù)據(jù)模型(Common Data Model)基礎(chǔ)上，建立類似pandas的工具，使分析多維數(shù)組變得簡單、高效、不易出錯。國家氣象中心預(yù)報技術(shù)研發(fā)室使用Python研發(fā)了全流程檢驗程序庫[10]，在數(shù)據(jù)層提供了數(shù)據(jù)讀寫、轉(zhuǎn)換等基礎(chǔ)函數(shù)；在檢驗算法層提供有無預(yù)報、多分類、概率預(yù)報等檢驗算法；在檢驗產(chǎn)品制作層提供了數(shù)值型、表格型、圖片型、誤差序列分析等函數(shù)，成為國內(nèi)第一款專門用于氣象預(yù)報檢驗的python程序庫。本文在檢驗結(jié)果可視化階段使用了該庫提供的功能。

基于機器學(xué)習(xí)方法的短臨多要素氣象預(yù)報系統(tǒng)(Weather Elements Nowcasting based on machine learning, 簡稱WEN)，每10 min滾動更新發(fā)布未來2 h多氣象要素，覆蓋中國區(qū)域的短時臨近預(yù)報,具有高更新頻次、高時空分辨率特點。系統(tǒng)在分布式計算環(huán)境下，建立72候(1年)×12時辰(1天每2 h為1個時辰)多預(yù)報時段復(fù)雜預(yù)報模型。如何快速獲得檢驗結(jié)果是本項工作面臨的挑戰(zhàn)。

如圖1所示，WEN檢驗子系統(tǒng)主要分為數(shù)據(jù)準備、檢驗計算和檢驗分析3個階段。第1階段進行數(shù)據(jù)讀取、合并，目標是設(shè)計通用的數(shù)據(jù)模型。第2階段進行檢驗計算，包括連續(xù)型數(shù)值預(yù)報檢驗、晴雨檢驗和降水分級檢驗3類檢驗方法，全部要素共計10個檢驗指標。第3階段分要素、區(qū)域?qū)EN模型進行性能評估。

圖1 WEN檢驗業(yè)務(wù)流程

1 數(shù)據(jù)模型

1.1 檢驗數(shù)據(jù)

多維數(shù)組(又稱n維，ND ，有時稱為張量)是計算科學(xué)的一個重要部分，在很多領(lǐng)域都有應(yīng)用，包括物理、天文學(xué)、地球科學(xué)、生物信息學(xué)、工程、金融和深度學(xué)習(xí)。在Python中，NumPy為處理原始ND數(shù)組提供了基本的數(shù)據(jù)結(jié)構(gòu)和API。然而，真實的數(shù)據(jù)集通常不僅僅是原始的數(shù)字，它們有標簽來編碼關(guān)于數(shù)組值如何映射到空間、時間等位置的信息。Xarray在原始的與Numpy類似的多維數(shù)組上引入了維度、坐標和屬性形式的標簽，使得開發(fā)更直觀、簡潔、更少出錯。其核心數(shù)據(jù)結(jié)構(gòu)Dataset被設(shè)計為NetCDF文件格式的數(shù)據(jù)模型的內(nèi)存表示。

WEN系統(tǒng)采用網(wǎng)絡(luò)通用數(shù)據(jù)格式(Network Common Data Form，NetCDF)，使用TDS工具(THREDDS Data Server，TDS)管理文件并提供數(shù)據(jù)訪問服務(wù)。預(yù)報和觀測數(shù)據(jù)(以下統(tǒng)稱檢驗數(shù)據(jù))均為每10 min發(fā)布一次。每個時次發(fā)布的站點預(yù)報數(shù)據(jù)包括預(yù)報時刻、站點、氣象要素3個維度，觀測數(shù)據(jù)包括站點和氣象要素2個維度。

為避免每次在線獲取耗時缺點，同時提高處理和共享能力，采用地球科學(xué)科研領(lǐng)域元數(shù)據(jù)標準CF(Climate and Forecast Metadata Conventions)[11]，將檢驗數(shù)據(jù)按天合并后保存到本地。數(shù)據(jù)模型構(gòu)建分析如下：①檢驗用的站點數(shù)據(jù)屬于離散采樣幾何類型(Discrete Sampling Geometries)；②一天多個時次更新發(fā)布，構(gòu)成了時間序列(Time Series)?；诖?，檢驗數(shù)據(jù)采用正交多維數(shù)組的時間序列數(shù)據(jù)表示方法，站點預(yù)報數(shù)據(jù)模型如圖2所示。預(yù)測數(shù)據(jù)模型去掉“預(yù)報時效”維度即為觀測數(shù)據(jù)模型。在計算機內(nèi)存中對應(yīng)Xarray的Dataset。在預(yù)測和觀測數(shù)據(jù)匹配時，通過表示位置的“station”和表示時間的“time”兩個“標簽”即可實現(xiàn)檢驗數(shù)據(jù)配對，體現(xiàn)了xarray直觀、簡潔的優(yōu)點。

圖2 WEN檢驗數(shù)據(jù)模型(按天合并的站點預(yù)報)

1.2 結(jié)果數(shù)據(jù)

檢驗方法不同對應(yīng)指標、要素不同。如晴雨檢驗和降水分級檢驗僅針對降水要素，包括TS評分、空報率等指標。連續(xù)型預(yù)報數(shù)值型檢驗包括平均誤差、均方根誤差等指標。WEN規(guī)定按候輸出檢驗結(jié)果?；诖?，檢驗結(jié)果按檢驗方法分類生成，在時間維度包括時效、時辰等。

1.3 站點數(shù)據(jù)

傳統(tǒng)檢驗方法基于點對點的對比，即在一個時間點上參與檢驗的樣本數(shù)為該區(qū)域的站點數(shù)。WEN系統(tǒng)包括105567個站點?？紤]觀測數(shù)據(jù)質(zhì)量，選取包括國家級氣象站(基準站、基本站、一般站)和質(zhì)控后的部分區(qū)域自動站，共計10656個站點。在此基礎(chǔ)上將全國分成8個區(qū)域，分區(qū)范圍見圖3，各分區(qū)站點數(shù)見表1。

圖3 檢驗分區(qū)示意

表1 分區(qū)站點數(shù)

2 檢驗計算

WEN采用連續(xù)型預(yù)報數(shù)值檢驗、有無預(yù)報檢驗(晴雨預(yù)報檢驗)和多分類檢驗(降雨分級檢驗)等3類檢驗算法。后兩類屬于分類問題，前者是常見的二分類，后者是類別數(shù)大于2的分類。在分類算法的實現(xiàn)中應(yīng)用混淆矩陣獲得分類樣本統(tǒng)計。針對耗時任務(wù)給出了并行解決方案。

2.1 檢驗方法

連續(xù)型預(yù)報檢驗采用平均誤差、平均絕對誤差、均方根誤差、偏差4種統(tǒng)計方法，檢驗變量包括氣壓、溫度、相對濕度、風(fēng)速等近地面要素。降水采用技巧評分(TS)、相當(dāng)技巧評分(ETS)、空報率、漏報率、預(yù)報偏差、準確率6種檢驗方法。上述方法已在業(yè)務(wù)中廣泛應(yīng)用，在此不做詳細介紹。計算公式主要參考由陳昊明等[12]的高分辨率區(qū)域數(shù)值預(yù)報業(yè)務(wù)檢驗評估指標及算法。降水分級檢驗主要參考關(guān)于下發(fā)中短期天氣預(yù)報質(zhì)量檢驗辦法(試行)的通知(氣發(fā)〔2005〕109號)中中短期天氣預(yù)報質(zhì)量檢驗辦法(試行)解釋說明(表2)和短時氣象服務(wù)降雨量等級劃分[13](表3)。

表2 降水分級檢驗評定簡表

表3 短時氣象服務(wù)降雨量劃分 mm

2.2 數(shù)據(jù)缺失處理

為保證檢驗結(jié)果的正確性及連續(xù)性，缺失數(shù)據(jù)的處理是必不可少的關(guān)鍵步驟。數(shù)據(jù)缺失有兩種表現(xiàn)，一是數(shù)據(jù)文件缺失，二是數(shù)據(jù)文件中存在空數(shù)據(jù)。第1種情況表明某時次沒有可用于檢驗的數(shù)據(jù)，可將檢驗結(jié)果置為一個約定值，例如-999。第2種情況則需要將空數(shù)據(jù)去除，從而保證統(tǒng)計樣本的有效性。需要注意的是空數(shù)據(jù)并非一一對應(yīng)，可應(yīng)用公共掩碼矩陣對檢驗數(shù)據(jù)進行剔除。技術(shù)實現(xiàn)的關(guān)鍵代碼見圖4。

圖4 缺失值處理關(guān)鍵代碼

2.3 混淆矩陣

混淆矩陣是機器學(xué)習(xí)中總結(jié)分類模型預(yù)測結(jié)果的情形分析表，以矩陣形式將數(shù)據(jù)集中的記錄按照真實的類別與分類模型預(yù)測的類別進行判斷匯總。其中矩陣的行表示真實值，矩陣的列表示預(yù)測值?；赑ython語言的機器學(xué)習(xí)工具scikit-learn的metrics模塊的混淆矩陣函數(shù)confusion_matrix實現(xiàn)了混淆矩陣的計算。通過傳遞真實類別和預(yù)測類別樣本，返回一個n階矩陣C(n表示類別數(shù)，如二分類n=2)，Ci,j表示真實類別為i預(yù)測類別為j的樣本數(shù)量。以二分類晴雨檢驗為例，C0,0表示真反例即命中否的樣本數(shù)，對應(yīng)表2中ND;C1,0表示假反例即漏報的樣本數(shù)，對應(yīng)表2中NC;C0,1表示假正例即空報的樣本數(shù)，對應(yīng)表2中NB;C1,1表示真正例即命中的樣本數(shù)，對應(yīng)表2中NA。二分類混淆矩陣與晴雨判定對應(yīng)關(guān)系見表4。通過將氣象領(lǐng)域分類檢驗的判定表與機器學(xué)習(xí)領(lǐng)域分類模型評估分析表即混淆矩陣建立對應(yīng)關(guān)系，應(yīng)用混淆矩陣函數(shù)實現(xiàn)分類檢驗計算公式中各參量，再根據(jù)公式最終獲得檢驗計算結(jié)果。晴雨檢驗計算實現(xiàn)的關(guān)鍵代碼見圖5。

表4 混淆矩陣晴雨檢驗對應(yīng)關(guān)系

圖5 晴雨檢驗計算實現(xiàn)代碼

2.4 并行計算

2.4.1 并行計算概念、方式、環(huán)境

從計算復(fù)雜性的角度來看，一個算法的復(fù)雜性可表示為空間復(fù)雜性和時間復(fù)雜性兩個方面。并行計算(Parallel Computing)的基本思想是用多個處理器來協(xié)同求解同一個問題，即將被求解的問題分解成若干個部分，各部分均由一個獨立的處理器或處理機來并行計算。并行計算的目標就是盡量減少時間復(fù)雜性，通常這是通過增加空間復(fù)雜性來實現(xiàn)的。并行計算可分為時間上的并行和空間上并行。時間上的并行就是指流水線技術(shù)。空間上的并行則是指多個處理器并發(fā)的執(zhí)行計算。并行環(huán)境多為MPP(Massively Parallel Processing，大規(guī)模并行處理)并行機。不同于普通PC機，通常并行機有特定的系統(tǒng)結(jié)構(gòu)，支持指定的并行環(huán)境(常見的有MPI(Message Passing Interface))，硬件價格昂貴，同時并行計算的實現(xiàn)也是一項復(fù)雜而艱巨的工作。

曹曉鐘等[14]針對大量神經(jīng)網(wǎng)絡(luò)訓(xùn)練、預(yù)報的需求，在SP2并行機上實現(xiàn)了基于主從計算模式，采用粗粒度任務(wù)劃分的并行方案。趙軍等[15]在多核高性能計算機上使用MPI/OpenMP混合并行計算編程，以增量方式和較小的并行化代價實現(xiàn)了對AREM模式并行化改進，取得了較好的并行計算效果。二者均是在特定并行機環(huán)境下實現(xiàn)的并行化。周欽強等[16]采用Java語言基于Eclispes平臺開發(fā)實現(xiàn)了基于TCP多連接通信實時并發(fā)數(shù)據(jù)處理系統(tǒng)，滿足對周期性大批量浪涌數(shù)據(jù)進行實時快速處理的要求。系統(tǒng)運行于單處理機上。其主要利用線程池化技術(shù)和并發(fā)管理機制，實現(xiàn)過程不僅需要對不同線程進行詳細的程序流程分析，還要處理復(fù)雜的線程安全問題，開發(fā)成本相對較大。韓帥[17]等利用多節(jié)點分塊并行與模式并行結(jié)合的計算方案，解決了系統(tǒng)文件大小受限、計算資源和計算時效等問題。是一種時間和空間二維并行的方案。

2.4.2 DASK并行計算框架

DASK是一款開源免費、與現(xiàn)有科學(xué)計算項目Numpy、Pandas、Scikit-Learn高度協(xié)同的基于python語言的并行計算框架[18]。它以一種更方便簡潔的方式處理大數(shù)據(jù)量，能夠在普通單機或集群中進行分布式并行計算。DASK由兩個部分組成，動態(tài)任務(wù)調(diào)度和“大數(shù)據(jù)”集合。動態(tài)任務(wù)調(diào)度包括單機調(diào)度器和分布式調(diào)度器，分布式調(diào)度器可以在單機上運行，提供了比單機調(diào)度器更豐富的功能，如診斷面板。在單機中通過scheduler參數(shù)設(shè)置線程或者進程來并行處理，也稱為偽分布式。“大數(shù)據(jù)”集合提供的3種基本數(shù)據(jù)結(jié)構(gòu)分別是arrays、dataframes、bags。Arrays是對Numpy中的ndarray的部分接口進行了改進，從而方便處理大數(shù)據(jù)量。Dataframes是基于Pandas DataFrame改進的一個可以并行處理大數(shù)據(jù)量的數(shù)據(jù)結(jié)構(gòu)。以上兩者通過對大數(shù)據(jù)量進行分塊處理，提供了大于內(nèi)存數(shù)據(jù)的解決方案。Bags主要用于半結(jié)構(gòu)化的數(shù)據(jù)集，比如日志或者博客等。DASK提出的延遲計算是在一個普通的Python函數(shù)上面通過dask.delayed函數(shù)進行封裝，得到一個delayed對象。該對象只是構(gòu)建一個任務(wù)圖(Task Graph)，并沒有進行實際的計算，只有調(diào)用cumpute的時候才開始進行計算。延遲計算非常適合當(dāng)面對的問題是可并行的但是不適合高度抽象如Dask Array或Dask Dataframe的場景。下面給出一段在單機上使用分布式調(diào)度器、延遲計算、多進程實現(xiàn)并進行化的代碼示例，見圖6。

圖6 一種DASK并行化示例代碼

2.4.3 WEN檢驗并行化實現(xiàn)

并行計算首先是一種解決計算問題的思想。分析WEN檢驗，以數(shù)據(jù)準備環(huán)節(jié)為例，處理過程分為兩個步驟：①讀取數(shù)據(jù)，這個操作對每個文件是相同的；②所有文件讀取結(jié)束后進行合并操作。讀取一個文件是一個獨立動作，互相之間無依耐關(guān)系，可同時執(zhí)行。利用dask.delayed將讀取文件進行并行化改造。以處理156個觀測文件(1 d+2 h)為例，程序分別在兩臺具備不同核數(shù)的普通PC機上運行?？紤]網(wǎng)絡(luò)帶寬時間維度不均勻性，為使測試結(jié)果更客觀，選取2020年9月1—5日，每個日期在不同時間分3次運行，分別記錄每次的運行時間，再計算平均用時。并行前后程序執(zhí)行時間見表5。以平均數(shù)衡量，56核處理機效率提升11倍，96核處理機效率提升9倍。進行并行化改造增加的代碼不超過10行。

表5 觀測數(shù)據(jù)采集合并運行時間測試結(jié)果 s

運用相同思路，對檢驗計算環(huán)節(jié)代碼進行并行化改造。以計算1天數(shù)據(jù)每個時次每個時效數(shù)值型預(yù)報檢驗為例(僅9月1日數(shù)據(jù)缺失36個時次，其它日期完整)，時次和時效兩個維度同時并行，即并行任務(wù)有144×12=1728個，測試結(jié)果見表6。并行后，56核環(huán)境效率提升20倍，96核環(huán)境提升51倍。改造后的關(guān)鍵代碼見圖7。

圖7 WEN數(shù)值型檢驗方法并行化實現(xiàn)關(guān)鍵代碼

表6 連續(xù)型數(shù)值型預(yù)報檢驗(每時次每時效)運行時間 s

從測試結(jié)果來看，效率提升并非與CPU核數(shù)成線性遞增。首先，該測試結(jié)果采用單機多進程并行方案，且并行任務(wù)數(shù)大于核數(shù)，即并行任務(wù)與CPU并非一次分配到位，而是進行了多次分配。其次，當(dāng)進程間傳輸數(shù)據(jù)較大時，會產(chǎn)生性能損失。從表5和表6并行方案測試結(jié)果來看，96核環(huán)境耗時均高于56核環(huán)境。但總體上計算效率并行方案都遠高于串行方案。

對于串行方案，表5測試結(jié)果顯示，56核環(huán)境略高于96核，但相差不大。主要原因是表5對應(yīng)的業(yè)務(wù)屬于IO密集型，主要受實時網(wǎng)絡(luò)流量影響，受CPU型號差異帶來的性能影響較小。表6測試結(jié)果，96核環(huán)境耗時是56核環(huán)境的3倍左右。表6對應(yīng)的業(yè)務(wù)屬于計算密集型，主要受CPU性能影響。56核環(huán)境CPU型號高于96核環(huán)境。兩者的型號分別是56 Intel(R) Xeon(R) Gold 5120 CPU @ 2.20 GHz和96 Intel(R) Xeon(R) CPU E7-8850 v2 @ 2.30 GHz。

3 檢驗分析

圖8是2021年3月1～6候模型的全國區(qū)域氣溫要素檢驗結(jié)果(篇幅有限僅以氣溫要素均方根誤差指標為例進行說明)。圖中橫軸表示一個候的預(yù)報模型，每個刻度代表一個模型(即每天12個時辰之一，世界時)。從圖中可以看出，①3月WEN模型氣溫均方根誤差最大4.77 ℃(第1候第11時辰)，最小1.13 ℃(第4候第5時辰)；②每候總的趨勢表現(xiàn)為：后5個時辰較前7個時辰大，即夜間(北京時晚10:00至次日早08:00)誤差較白天大。對比了21年3—5月的檢驗結(jié)果，發(fā)現(xiàn)不同月份表現(xiàn)差的模型均集中在夜間。

圖8 2021年3月第1～6候WEN模型全國氣溫均方根誤差

4 結(jié)論

本文在分析業(yè)務(wù)需求的基礎(chǔ)上，設(shè)定快速計算為目標，重點圍繞氣象領(lǐng)域檢驗分類問題和應(yīng)用并行計算改善時效性兩個方面，設(shè)計并實現(xiàn)了針對高頻次、高時間分辨率短臨氣象預(yù)報產(chǎn)品的檢驗系統(tǒng)。今后將圍繞檢驗報告自動生成、容器化部署等開展工作，推動公共氣象服務(wù)中心產(chǎn)品檢驗業(yè)務(wù)常態(tài)運行，為實現(xiàn)氣象服務(wù)高質(zhì)量發(fā)展貢獻力量。