GRAPES-GEPS K-均值集合預(yù)報產(chǎn)品 開發(fā)及應(yīng)用

齊倩倩 佟華* 陳靜

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

引言

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

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

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

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

1 資料及天氣過程介紹

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

2 聚類分析產(chǎn)品

2.1 K-均值聚類方法介紹

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

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

(1)

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

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

(2)

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

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

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

2.2 500 hPa位勢高度場聚類產(chǎn)品

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

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

(5)

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

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

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

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

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

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

2.3 850 hPa溫度場聚類產(chǎn)品

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

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

圖5 GRAPES-GEPS的2020年2月10日12:00起報，96 h預(yù)報的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預(yù)報的850 hPa溫度場的集合預(yù)報離散度(b)

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

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

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

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

2.4 10 m風(fēng)速聚類產(chǎn)品

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

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

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

3 結(jié)論和討論

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

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

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

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

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