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的有效性，还需采用更多的天气过程进行检验。另外，受集合预报样本数的限制及实际天气变化复杂性的影响，在更长时间的预报时效上，集合预报本身预报能力也是有限的，因此相应聚类产品的展示效果也会受到影响。其次，针对不同的天气形势，需采用不同的聚类方法产生聚类产品。在未来，有必要将更多的聚类方法应用在各类天气形势中，从而找到理想的聚类方法，以满足对各类天气形势的预报效果。