GRAPES-GEPS K-均值集合預報產品 開發及應用

齊倩倩 佟華* 陳靜

(1 中國氣象局數值預報中心，北京 100081; 2 國家氣象中心，北京 100081)

引言

受觀測資料、模式誤差及大氣混沌特性的影響，數值天氣預報存在著較大的不確定性。而集合預報是估計這種預報不確定的重要工具之一[1-3]。它通過對初始誤差、模式誤差的概率分布進行采樣，并多次積分數值天氣預報模式來獲取未來大氣可能狀態的概率分布，從而使單一確定性預報變為概率預報。如今，集合預報已成為數值預報，尤其是中期預報的重要組成部分之一。

集合預報包含了大量的信息，為更有效地為預報員們所利用，有必要從中識別和提取有用的預報信息。杜鈞[4]提出，集合預報產品大致分為3類：①集合平均或中值產品；②大氣可信度預報產品，通常用集合成員的離散度衡量；③概率預報產品。集合平均產品可濾掉預報樣本的隨機信息，比單一預報更準確，但不能包含預報的所有可能性。對于大氣可信度預報產品，集合預報成員間的離散度可量化天氣預報的可信度。概率預報產品能最大程度地包含實際大氣可能發生的全部情況。目前，國內對集合預報產品的應用技術已相當成熟，中央氣象臺業務內網上常用的集合預報產品有：集合平均和集合離散度產品、面條圖、郵票圖、概率分布平面圖、箱線圖等[5]。然而和其他一些國家如美國NCEP/NCAR和歐洲中心ECMWF數值預報產品相比，我國集合預報的產品種類還較為單一，不夠豐富[6-8]。因此有必要開發出新的集合預報產品，從而更充分地提取和利用集合預報樣本中的信息，以供預報員參考使用。

對于集合預報輸出的各個預報樣本，要在有限的時間內，找出不同預報結果的相似和差異特征是不現實的。聚類分析方法可對某一預報要素的不同集合成員進行分類，其中相似的成員歸為一類。通過執行聚類方法，可給出不同集合成員的相似性和差異性空間分布特征。

預報員可通過聚類產品大致了解所有集合樣本的預報分類情況，借此判斷天氣系統未來可能的發展方向，特別是對其中出現的異常情況，一目了然，起到提醒作用。目前，許多國家氣象業務部門已在使用聚類分析產品。ECMWF最初使用Ward聚類方法[9]，后改為管子法(Tubing clustering)對集合預報產品進行分類[10]，管子法給出了天氣預報的主要類別，并提供了距離平均最遠的極端天氣圖；NCEP采用距平相關系數分簇法對500 hPa高度場進行分組；法國氣象局采用Diday提出的動力模糊法[11-12]，初始劃分時的聚類中心由天氣類型確定，劃分用到的距離是位移和最大相關方法；瑞典氣象局使用神經元聚類法，該方法基于神經網絡原理；日本氣象廳采用中心聚類法(Central clustering)，即重心法，利用不同類重心間的距離作為標準進行劃分。然而，中國氣象局數值預報中心的GRAPES集合預報系統尚未開發和使用這些產品。K-均值聚類分析方法是由MacQueen提出的一種傳統的客觀分類方法[13-14]，即將一組個體按照相似性歸成若干類。該算法簡單有效，并能處理大量數據，尤其當樣本分布呈現類內距離小、類間距離大的聚類結構時，聚類效果較好。目前該方法已成功應用在多個行業[15-17]，如機器學習，數據挖掘，模式識別等多個領域。為此，本文將基于中國氣象局數值預報中心獨立自主研發的GRAPES全球集合預報系統(GRAPES-GEPS)[18-24]，開發K-均值聚類產品，并將該集合預報產品應用在一次具體的天氣過程中，給出初步的分析結果。

1 資料及天氣過程介紹

GRAPES全球集合預報系統(GRAPES-GEPS)是在GRAPES全球數值預報模式的基礎上，由中國氣象局數值預報中心獨立自主研發的預報系統。該系統采用奇異向量法產生初值擾動[21-22]，使用SPPT方法和SKEB擾動方案產生模式隨機擾動[23-24]，共生成31個集合預報成員，包括30個擾動預報和1個控制預報。集合成員水平分辨率為50 km，即0.5°×0.5°，垂直60層，輸出范圍為全球。該業務系統每天00:00UTC和12:00 UTC起報，預報時效為360 h，其中，84 h以內預報模式輸出的時間間隔為3 h，84～192 h模式輸出間隔為6 h，192～360 h模式輸出間隔為24 h。

2 聚類分析產品

2.1 K-均值聚類方法介紹

設聚類的樣本集為：X={x1,x2,…,xn},n為樣本總數，令Aj(j=1,2,…,K)表示聚類的K個類別，nj(j=1,2,…,K)表示第j類的樣本數，且滿足n1+n2+n3+…+nK=n。

記K個聚類中心為{z1,z2,…,zK}，則有

(1)

K-均值聚類算法通過不斷地迭代，最終確定穩定的聚類中心和聚類中心附近的樣本。具體地，其基本步驟如下[26-27]：①從所有數據中隨機選擇K個樣本z1(1),z2(1),…,zK(1)作為初始聚類中心；②計算每個樣本xi(i=1,2,…,n)到K個種子聚類中心的歐式距離，按照距離最近的原則，將每個樣本分配到歐氏距離最小的聚類中心的類中，從而對樣本集合聚類，確定每個樣本的類屬關系；③按照式(1)重新計算新的聚類中心z1(m),z2(m),…,zK(m),m表示迭代次數；④重復執行步驟2～4，直到聚類中心穩定或達到指定的迭代次數。

K-均值算法思想較為簡單，但確定一個合適的分類數量K對聚類效果的好壞有很大的影響。關于聚類數量的確定有以下方法：經驗法、密度法、逐個歸類法、遞推法和爬山法[28-30]。爬山法，又稱為肘方法，是確定最優聚類數的邏輯判定方法。該方法簡單明了，易于實現，已被廣泛應用在K-均值聚類算法中，以確定最佳的聚類數量[31-32]。為此，本文將采用爬山法來確定分類數量K。具體地，采用誤差平方和(Sum of the Squared Errors, SSE)SSE作為約束條件，該聚類判定函數定義如下:

(2)

其中，K是聚類數量，p為類別Ak的數據元素，zk是第k個聚類的中心點，即為類別Ai的均值(中心)，SSE為每個樣本數據點到相應聚類中心的距離平方和，即分為K類后，數據與各自所屬中心的誤差平方和。

K越大，每個類別中數量越少，類內部距離越近，此時SSE越小，說明樣本聚合程度越高。當K小于“真實聚類數”時，隨著K的增大，每個簇的聚類程度將大幅增加，故SSE下降幅度很大。而當K到達“真實聚類數”時，再增加K得到的聚合程度回報會迅速變小，此時，SSE的下降幅度逐漸變小，并最終趨于平穩，即認為再增加聚類數量，不能增強聚類效果。因此，最合適的分類數量K應選為SSE下降幅度變慢，并最先趨于平穩的點。

本文將基于GRAPES-GEPS預報資料的500 hPa位勢高度場、850 hPa溫度場和地面要素10 m風場，首先利用樣本聚類判定函數SSE(K)確定最合適的分類數量K，然后分別對其執行K-均值聚類，進一步地，結合該時段的天氣過程對該聚類產品給出分析。

2.2 500 hPa位勢高度場聚類產品

結合2020年2月13—16日的寒潮天氣過程，本文將選取GRAPES-GEPS的2月10日起報，超前96 h預報的500 hPa位勢高度場資料，針對中國及中國以北的部分區域(15°～75°N，70°～140°E)，做K-均值聚類分析。首先，對于固定的聚類數量K, 通過K-均值迭代算法(步驟1～4)確定穩定的聚類中心及聚類中心附近的樣本；其次，基于迭代的聚類中心及聚類中心周圍的樣本，利用式(2)，計算誤差平方和SSE函數。當聚類數量變化時，SSE隨之變化，進一步給出誤差平方和SSE隨聚類數量K的變化，如圖1所示。然后，采用爬山法確定最佳聚類數量。隨著K的增大，SSE逐漸降低，當K>3時，SSE降低幅度變小并逐漸趨于穩定，因此，最佳聚類數量K應為3。對該500 hPa位勢高度場預報資料做K-均值聚類，并分為3類，同時，我們也計算了每一類別的概率，按照概率從高到低排列，對應的每一類別的聚類中心如圖2所示。為了便于分析聚類產品，給出了2月10日12:00UTC時次起報，96 h預報的集合離散度及2月14日12:00UTC時的500 hPa位勢高度場實況，如圖3所示。進一步地，為了更客觀地衡量不同類別的聚類產品與實況的相似程度，我們采用相似系數來度量。相似系數的計算公式如下[2]：

圖1 500 hPa位勢高度誤差平方和(SSE)隨聚類數目K的變化

(5)

其中，ei表示第i個空間場，是不同空間場ei和ej的內積，且‖ei‖=,i=1,2。

李彬等[28]研究膜-生物反應器處理高鹽廢水結果表明，污泥中無機成分含量增加，絮體更為緊密，沉降性能變好，膜面污染物的成分為蛋白質、糖類和腐殖酸等。

2月14日12:00實況結果表明，亞洲中高緯度地區為“一槽一脊”型，內蒙古上空的后部高壓脊在貝加爾湖附近形成阻塞高壓，其中心強度為556 dagpm，鄂霍次克海地區有低壓冷渦系統，其中心強度低于512 dagpm。低壓系統后部有非常強的冷平流，有利于引導高緯度冷空氣向東北地區堆積。對比圖2和圖3，結果表明，通過對所有集合預報成員96 h預報時效的500 hPa位勢高度場進行3類的K-均值聚類，基本上所有類別都能呈現出中高緯500 hPa位勢高度Ω形的環流形勢及低壓系統后部冷平流的環流走向，這與實況相符。另外，表1的結果也表明，所有類別的聚類中心與實況的空間相似系數均在0.90以上。實況中呈現的貝加爾湖附近高壓脊形成的阻塞高壓形勢，在發生概率最高的第1類聚類類別中有所體現，同時該類別中高壓脊的強度、范圍及背后的低壓系統走向均與實況更為吻合，其空間相似系數達到0.978。發生概率最低的第3類聚類類別，其高壓脊范圍與實況比偏北、偏東，強度偏強，其脊后的低壓系統比實況偏西、偏強。發生概率次高的第2類聚類類別，其環流形勢與實況的相似度介于第1類和第3類之間。基于GRAPES-GEPS采用的K-均值聚類法的聚類產品，其最高發生概率對應的環流形勢與實況非常接近，最低發生概率對應的環流形勢最偏離實況，這表明，該方法能有效地劃分出最有可能發生的環流形勢，給預報員提供有價值的預報信息，同時也反映出采用K-均值聚類分析方法開發出的聚類產品是可行且合理的。

圖2 GRAPES-GEPS的2020年2月10日12:00起報， 96 h預報的500 hPa位勢高度場(dagpm)的K-均值 聚類：(a)第1類聚類中心(概率:47%),(b)第2類聚類 中心(概率：44%),(c)第3類聚類中心(概率:9%) (時間為世界時，下同)

圖3 GRAPES-GEPS的2020年2月10日12:00時起報， 96 h預報的500 hPa位勢高度場的集合離散度(陰影) 及2月14日12:00時500 hPa位勢高度場實況 (黑色實線,單位：dagpm)

表1 GRAPES-GEPS 500 hPa位勢高度不同類別的聚類中心

2.3 850 hPa溫度場聚類產品

本文仍選取GRAPES-GEPS的2020年2月10日起報，96 h預報的850 hPa溫度場資料，針對包含中國的區域(15°～55°N，70°～140°E)，做K-均值聚類分析。首先，采用樣本聚類誤差平方和的方法，計算得出最合適的聚類數量K應為5，如圖4所示。對該850 hPa溫度場資料做K-均值聚類，并分為5類，其中這5類的概率分別為：37%、32%、17%、13%和1%，結果如圖5所示。同樣地，我們給出2月10日12:00UTC起報，96 h預報的集合離散度及2月14日12:00UTC時次的850 hPa溫度場實況，如圖6所示。為了便于進一步地對比該聚類產品和實況的差異，圖7給出了K-均值聚類結果與對應時刻850 hPa溫度場實況的差的空間分布。

圖4 850 hPa溫度誤差平方和(SSE)隨聚類數目K的變化

圖5 GRAPES-GEPS的2020年2月10日12:00起報，96 h預報的850 hPa溫度場的K-均值聚類：(a)第1類聚類中心(概率：37%),(b)第2類聚類中心(概率:32%),(c)第3類聚類中心(概率：17%),(d)第4類聚類中心(概率:13%),(e)第5類聚類中心(概率:1%)

圖6 2020年2月14日12:00UTC 850 hPa溫度場實況(a)及GRAPES-GEPS的2月10日 12:00UTC起報，96 h預報的850 hPa溫度場的集合預報離散度(b)

圖7 GRAPES-GEPS的2020年2月10日12:00起報，96 h預報的850 hPa溫度場的K-均值聚類結果與2月14日12:00 850 hPa溫度場實況的差：(a)第1類,(b)第2類,(c)第3類,(d)第4類,(e)第5類

通過對所有集合預報成員96 h預報的850 hPa溫度場進行5類的K-均值聚類，基本上所有類別都能呈現出全國2 m溫度緯向梯度的空間分布特征，即全國溫度從北到南呈帶狀逐漸增加，東北及內蒙古大部分地區溫度在-10 ℃以下，與實況相符。從96 h預報的集合離散度可看出，2月14日12:00，甘肅中東部、陜西、內蒙古、山西、河北、河南、山東、東北等大部分地區預報離散度較高，而5類的K-均值聚類分析產品，其溫度分布在上述相應的地方也有較大差異。

具體地，通過分析K-均值聚類結果與對應時刻850 hPa溫度場實況的差的空間分布，發現，發生概率最高及次高的第1類和第2類聚類產品，和實況相比，其誤差均較小，且第1類產品最接近實況。發生概率最低及次低的第5類和第4類聚類產品，在東北、華東、甘肅等大部分地區，與實況相比，存在大的冷偏差。發生概率位于中間的第3類聚類產品，其在華中大部分地區存在暖偏差。因此，通過整體權衡不同類別的產品，發生概率最高的第1類聚類產品最能反映實況中850 hPa溫度的空間分布，而發生概率最低的第5類聚類產品最偏離實況，誤差最大。這表明，K-均值聚類分析方法開發出的該時次溫度預報的聚類產品是可行且合理的，預報員可通過某一時次預報的發生概率最高的聚類產品，并結合其他不同類別的產品進行綜合分析，從而判斷某一時次預報的溫度。

和實況相比，即使是發生概率最高的第1類聚類產品，其在東北、華東等地區仍存在一定程度的冷偏差，在華北等部分地區存在一定程度的暖偏差。一方面可能是因為，集合預報系統在某些位置處本身存在系統偏差，基本上所有的預報樣本都不能準確描述某些位置的溫度，所以即便聚類算法之后，這些系統偏差依然會存在；另一方面，可能是因為聚類產品是類別內集合成員的算術平均，這就有可能損失部分有用的信息，個別成員預報有可能會優于聚類分析產品。

2.4 10 m風速聚類產品

對于地面要素10 m風速，同樣地，我們計算了SSE隨聚類數量K的變化，進而分析出最合適的聚類數量K為3，圖略。基于GRAPES-GEPS的2月10日12:00UTC起報，96 h預報的10 m風速預報資料，針對包含中國的區域(15°～55°N，70°～140°E)，執行K-均值聚類算法，并分為3類，其中這3類聚類產品發生的概率相當，分別為：38%、33%和29%，結果如圖8所示。同樣地，作為參考，圖9給出2月10日12:00起報，96 h預報的集合離散度及2月14日12:00 UTC時次的10 m風速實況。為了便于進一步地對比10 m風速聚類產品和實況的差異，圖10給出了K-均值聚類結果與對應時刻10 m風速實況的差的空間分布。

K-均值聚類產品及實況和集合離散度的結果表明，基于集合預報樣本所分的3個類別都能呈現出全國10 m大風速的空間分布特征(圖8)，即全國10 m風速大值區主要集中在西藏、青海、華北平原等大部分地區，其中天津10 m風速可達13 m/s以上，與實況分布情況相符(圖9a)。在風速較大處，集合樣本離散度也較大，尤其是東部渤海、黃海、東海等海上的10 m風速離散度大于陸地 (圖9b)，因此，不同類別的風速在東部海區和天津等地，其大小差異較為顯著。具體地， 通過對比K-均值聚類產品對應時刻10 m風速實況的差的空間分布(圖10)，發現，這3類聚類產品與實況相比，在華北平原等大部分地區呈現負偏差，而在黃海等地區呈現強的正偏差。3類產品正偏差程度較為一致，而發生概率相對較高的第1類和第2類聚類產品，其負偏差程度相對較低。尤其是第1類聚類產品，其對天津及周邊10 m風速的空間分布及強度均描述非常準確。這表明，基于K-均值聚類法的GRAPES-GEPS聚類產品，能夠方便有效地實現集合預報10 m風速樣本的濃縮信息，并能提供有價值的預報信息。

圖10 GRAPES-GEPS的2020年2月10日12:00起報，96 h 預報的10 m風速的K-均值聚類結果與2月14日12:00 10 m風速實況的差：(a)第1類,(b)第2類,(c)第3類

3 結論和討論

本文基于GRAPES全球集合預報系統(GRAPES-GEPS)，采用爬山法確定最合適的聚類數量，并采用K-均值聚類算法產生聚類產品，進一步分析了該集合預報產品在2020年2月13日至2月16日全國寒潮天氣過程中的實際預報價值和應用效果。

(1)該方法的500 hPa位勢高度場聚類產品的所有類別都能呈現出中高緯Ω形的環流形勢及低壓系統后部冷平流的環流走向。另外，該方法能有效劃分出發生概率最高的環流形勢的分布和走向，為預報員提供有價值的預報信息。

(2)基于K-均值聚類方法的850 hPa溫度場聚類產品，所有類別都能呈現出全國2 m溫度從北到南呈帶狀逐漸增加的空間分布特征；在離散度較大處，每一類別差別較大；發生概率最高的第1類聚類產品最能反映實況中850 hPa溫度的空間分布。

(3)基于K-均值聚類方法的風速聚類產品，在風速較大處，集合樣本離散度也較大，不同類別的風速大小差異也較為顯著；發生概率較高的第1類聚類產品，其對天津及周邊10 m風速的分布及強度描述均較準確。

本文基于GRAPES-GEPS開發了K-均值聚類分析產品，并將其應用在2020年2月13日至16日的一次全國寒潮天氣過程中。結果表明，該集合預報產品對寒潮天氣預報的指導是有效且合理的，且可為預報員提供有價值的預報信息。然而，本文僅針對一次重大天氣過程的個例進行分析，全面評價集合預報K-均值聚類產品在GRAPES-GEPS的有效性，還需采用更多的天氣過程進行檢驗。另外，受集合預報樣本數的限制及實際天氣變化復雜性的影響，在更長時間的預報時效上，集合預報本身預報能力也是有限的，因此相應聚類產品的展示效果也會受到影響。其次，針對不同的天氣形勢，需采用不同的聚類方法產生聚類產品。在未來，有必要將更多的聚類方法應用在各類天氣形勢中，從而找到理想的聚類方法，以滿足對各類天氣形勢的預報效果。