也谈面向对象预报评分方法

持续关注学术前沿的读者应该会记得，很久之前我写过一篇题为《MODE预报评分方法》的介绍，讲解了MODE这种不...



持续关注学术前沿的读者应该会记得，很久之前我写过一篇题为《MODE预报评分方法》的介绍，讲解了MODE这种不同于传统TS、ETS或者FSS评分的预报检验算法。最近为了写代码看了一些文献，突然发现这种面向预报对象的评分方法正逐渐成为检验预报结果的主流算法，受到越来越多的科研工作者乃至业务工作者的青睐。究其原因，当然是因为面向对象的预报评分方法能更全面、客观地评估预报结果，通过模糊的判据得出接近人脑思维的结论而非简单地计算一个数值。（是不是感觉和现代天气预报从确定性预报走向概率预报有异曲同工的地方？）

而进一步的查阅资料我发现，国外气象界对面向对象评分方法的研究，已经不单单是基于最早的那几个算法（MODE、CRA、SAL），甚至已经萌生了各种各样的理论，以适用于不同的情况。

就连一开始只能应用于二维场的MODE，也开始向三维（时间轴）进化：

而反观国内，似乎传统的检验方法依然是业内主流，如此对比之下未免感受到巨大差距。所以个人觉得有必要对这些面向对象的预报检验方法进行介绍，至少是抛砖引玉吧。

不过呢，就像上面所说，面向对象的检验方法已经不少，这里不可能一一列举，我只能挑出其中个别较为有名的进行介绍。而且，虽然这些方法各有各的特点，但万变不离其宗，总体的设计思路都是类似的，希望看完后读者也能形成自己的一套面向对象的预报检验思路，进而开发、设计属于我们自己的算法。

一、MODE

MODE（Mothod for Object-Based Diagnostic Evaluation），这个名字相当点题：基于对象的诊断检验方法（真是非常直白的名字啊Orz），算法的详细介绍可在Davis的文章：《Object-Based Verification of Precipitation Forecasts》Part I-II中查看，文章已经全网免费。

MODE作为一个出现较早的面向对象预报检验方法在如今却依然受到广泛使用，经得起时间的考验：这与它严谨又巧妙的设计不无关系。该算法可分为如下几步：

1、获得观测场和预报场中的预报对象

对象对象，说到底就是预报员关注的天气系统、强度中心。对于这样的对象，人类可以一眼就识别出来，但机器却做不到，所以必须使用特别的算法。最简单的操作就是设定阈值，然后将不满足阈值的区域遮罩掉。但这样一来，可能会进入大量的孤立中心，从而干扰计算机的识别。为了解决这一点，MODE算法中使用了卷积对格点场进行平滑处理，将范围较小的孤立对象全部滤去，然后再用阈值法进行过滤：



如上图d，经过处理后格点场中的主要对象也被分离了出来。

光分离没有用，我们还需要让计算机对这些对象进行识别和标记。这里使用的就是图像处理领域的算法：“图像连通域标记算法”。（实际上绝大部分面向对象的检验方法都利用了图像处理领域的知识）算法思路很简单，就是根据像素点是否连接做连通域判断，然后进行标记。作为一种相当经典的图像处理方法，应该不难查阅资料，这里就不多赘述。

2、设计相似系数，得到相似度矩阵

得到了预报场与观测场中的对象，下一步就是计算相似度。MODE算法在此处的处理较为自由，使用者可以根据自己的需求或设想进行设计。但最后都应以如下方程的形式进行累加：F代表相似方程，它计算得到的是两个对象之间的相似系数（结余0-1之间），w则是该相似系数在总相似度中所占的权重，c代表了该相似系数的置信度，M为相似系数的总个数。

常用的相似系数有：面积比、重叠率、轴角差、几何中心距离、边界最短距离。每一种相似系数即其权重和置信度都可以由使用者按照自己的需求进行设计，例如文章《the method for object based diagnostic evaluation applied tonumerical forecasts from the 2005 NSSL/SPC spring program》中的相似系数设计：这里我一定要狠狠地吐个槽，绝大部分使用MODE检验方法的作者在设计相似系数和置信度的时候都是用图像的形式而非具体公式，搞得我在编码的时候还要看着图片的坐标来推算方程！

3、计算相似度矩阵，分析结果

得到相似系数和置信度，接下来就能一一计算出预报场和观测场中对象的总相似度，注意是一一计算。预报场中三个对象，观测场两个对象，那么就应该计算出6个总相似度，从而得到这样的矩阵：有了这个矩阵，不说预报员，就连机器也能对预报结果做出较为清晰的判断。下面这个图是笔者之前绘制的MODE算法检验结果，可以发现，预报场中有四个对象（四个预报系统），观测场中只有三个，而根据最右得出的相似度矩阵可知，预报对象2与观测对象1的相似度较高，预报对象1与观测对象2相似度较高，预报对象4与观测对象3相似度较高，但不如前两组对象。这说明预报结果中对观测到的三个系统模拟都较好，但其中对观测场中第三个系统的预报有较大误差。同时还可以发现，预报对象3跟所有观测对象都不能进行较好的匹配（相似度较低），所以它是一个空报。



从上面这段分析读者应该能发现，面向对象的检验方法并不是为了得到一个评估整体预报效果的数值，而是方便预报员进一步从细节入手分析预报结果，甚至是方便计算机智能地自动对预报结果进行分析，从而极大地节省劳动力。

二、CRA

CRA（Contiguous Rain Area）连续雨区法。这个算法最早是由澳大利亚学者Ebert设计，具体参考文献可见：《Verificationof precipitation in weather systems: determination of systematic errors》和《toward better understanding ofthe contiguous rain area(CRA) method for spatial forecast verification》，均已全网免费。CRA是一种比MODE还要老牌的面向对象预报检验方法，其主要算法流程如下：

1、确定连续雨区

什么叫连续雨区？就是大于某个阈值的降水区域。不过CRA中的连续雨区并不是基于观测场或者预报场，而是基于二者的融合场：

CRA=MAX（CRA_obs+CRA_fcst）

即是说，将观测和预报场按最大值进行叠加，在叠加后的二维场中获得满足阈值条件的连续雨区。而CRA中的对象，就是这些连续雨区。识别连续雨区的方法与MODE相同，只不过CRA不适用卷积进行平滑，而是根据使用者设定，舍弃面积（像素点数量）较小的连通域，只保留面积较大的对象。

2、移动预报场

虽然是基于预报、观测融合场得到连续雨区，但连续雨区中各部分的隶属依然要明确，如下图中，预报是横线阴影，观测是竖线阴影：



明确了隶属关系，下一步就是将预报场向观测场移动。简单来说就是把上图中横线阴影的区域整体地平移到新的区域。移到什么位置停止呢？Ebert认为最小误差和最大相似系数都可以作为停止判据。我个人比较喜欢最大相似系数判据。具体操作就是：预报雨区每移动一次，得到一个新的连续雨区，从而能计算出这个新连续雨区中的相似系数，然后再移动，再计算，直到遍历全场。

当然，这样算法效率很慢，但奈何官方也是如此计算，同时官方代码中也对此有所吐槽Orz：



需要注意的是：移动的并非整个预报场，而是以某一连续雨区对象为中心的小区域。这样能极大地减少计算量（使用者还可以设计很多其他判据优化代码）。

3、计算误差

CRA方法中将误差拆分为三个部分：位移误差（displacement）、雨量/强度误差（volume）、形势误差（pattern）：



其中，位移误差为位移前后连续雨区中预报场同观测场的均方差（不开根！）：



雨量误差为移动后的预报场在连续雨区中的平均值同观测场平均值的误差平方：



形势误差则是用位移后连续雨区内预报误差减去雨量误差：



最后CRA输出的并非一个整体误差，而是以上三种误差各自所占比例，由此可看出预报结果到底哪里出了问题，方便预报员对结论进行分析和订正。另外，也可根据雨区面积大小以及对象匹配判断是否出现空报、漏报。

然而，这里我很想吐个槽：

实际上与CRA有关的文献中对这几种误差的计算都描述的相当模糊不清，甚至有前后矛盾的情况出现。包括官方CRA代码，实际上存在很多问题与不合理的地方。这也是我不喜欢CRA检验方法的原因：这个算法本身设计不够严谨，规则制定并不完善。

更不用说那个慢得要死的位移方法了Orz。

三、SAL

SAL（Structure，Amplitude and Location）检验方法，详细内容可参见文章：《SAL-ANovel Quality Measure for the Verification of Quantitative PrecipitationForecasts》。

说道SAL检验，我感动的快要哭了。在饱受MODE和CRA折磨虐待后，SAL是算法流程、算法思路、公式写得最清楚的一个面向对象检验方法！几乎详细到每一步怎么做，简直业界良心！

与CRA类似，SAL最终目的也是计算三个误差：structure（S），amplitude（A）和location（L）。其算法流程如下（注意下文中的R均代表对象区域，D代表整个网格区域）：

1、划分目标对象

与之前相同，第一步依然是划分目标对象。文献中所用为简单的阈值法：



其中上标max代表区域中降水最大值。其实用多少阈值都是用户自己决定，识别预报对象、观测对象的算法在所有的面向对象检验方法中都差不多。

2、计算S、A、L

我们按照文中顺序来。首先计算的是A：



这里的D代表计算整个网格空间的数值平均（注意不是某一对象）：



接下来计算L：



其中：



代表标准化后的质心距离。X代表计算降水质心坐标（整个预报网格的质心和整个观测网格的质心，而非某一对象质心），d表示整个网格空间的最大可能距离（网格对角线之类的）。

接下来又计算：



其中



这里Rn表示对象区域中的降水总量，r计算的是某一网格内所有对象质心到整个网格空间质心距离的加权平均。这样就能计算得到观测空间的r和预报空间的r两个数值，再利用：



进行标准化，得到L2

最后我们来计算S：

首先计算：



这里
代表对象区域内的最大值。接下来计算：



即用刚刚得到的Vn为权，计算整个网格区域内的对象区域总降水量加权平均。最后：



到分别将上述算式应用到观测网格和预报网格，计算出参数S。

3、分析结果

从上可知，SAL检验方法与CRA和MODE不同，它基于网格中的对象，最后得到了三个能够对整个网格预报结果进行描述的数值，而非将对象拿出来做一对一对比。参数A描述了预报场的平均降水预报准确度，参数L描述了预报场的平均落区准确度，参数S表示了预报场对降水区域结构的预报是否准确。三个参数的完美值都是0。

哎，SAL真是简单又方便

四、小结

上面介绍了三种面向对象的预报检验方法，虽然各有各的特点，但同时也可发现，这三种方法的第一步都是相同的：获得网格中的对象。针对对象进行预报检验是面向对象检验方法的精髓。

确认了对象后，接下来的操作实际上有着很大自由度。用户可以自行设计误差判据，可以计算观测和预报场中对象一一对应的关系、也可以计算整体特征。甚至可以把上面这些方法进行融合使用：例如在使用CRA进行检验时，位移停止判据可以不用相关系数，而采用计算MODE中的相似度作为判据。

个人认为，预报检验方法的相互融合是今后发展的大趋势，就像资料同化的hybrid一样，取长补短才能进步。例如我可以把CRA同DAS做个结合：将DAS算法中的位移方法应用到CRA中。为了方便说明，这里就简单提一下DAS。

DAS（Displacement and Amplitude Score），详细算法流程可见Keil等的文章《A displacement-based error measure applied in a regional ensembleforecasting system》和《a displacement and amplitude score employing an optical flowtechnique》。

算法的一大核心就在预报场位移上。在CRA中，位移是简单的逐格点移动，非常的费时费力；而DAS则首先将原始像素网格粗化。具体做法就是每隔特定距离保留网格点，并计算该像素点为中心的特定范围网格平均值作为当前格点值。就像下面这个梗一样：



经过粗化后的网格数量得到极大减少，然后在这样的网格空间内计算位移向量。其实与其说是位移，不如说是形变。具体算法为：确定某一像素（网格点），并将整个网格做正负2格范围内的移动，并计算该像素点为中心的5x5区域与观测场的方差。当方差最小时停止移动，此时的位移向量为该像素点的位移向量。以此类推可得到所有像素点的位移向量。然后将最终得到的位移向量场应用到原始（未粗化）网格。接下来按较小比例再次粗化，重复上述操作（即进行微调）。最后把每一次得到的位移向量加起来，计算出最终的位移向量。

我们可以思考一下，如果把这个位移方法的思路应用到CRA里，便可这样操作：移动区域限定在连续雨区，然后进行粗化，接着移动整个连续雨区网格，当满足条件（如方差最小、相关系数最大等）后停止位移，把位移向量按粗纹理化的比例放大，应用到原始网格空间；然后按较小比例再次粗化，重复上述操作（即进行微调）。最后把每一次得到的位移向量加起来，计算出最终的位移向量。这样一来计算量将远小于逐网格点位移。

而DAS最后计算的AMP与DIS误差参数的算法，也可以同SAL进行结合。CRA如果想要计算整个预报场的误差状况，也可以使用SAL的加权平均的思路。

类似的组合方式是无穷无尽的，就怕你想不到。这也是面向对象检验方法的一个特色：没有绝对固定的套路，规矩都是模糊的。模糊的规则才更接近人脑的思维模式，充满可变通性。

当然了，并不是说面向对象的检验方法会取代传统的检验算法。上述的检验方法都并非完美，具体使用还应根据情况而选择。在这一点上，传统检验方法便略胜一筹：简单粗暴的特点，使其在各种情况下都能发挥作用（至少能计算出一个相对合理的数值）；另外，各种各样的判据和对象变化也使得面向对象检验算法变得非常繁琐，编码极度复杂。

不过呢，笔者依然相信在未来的气象领域，面向对象的检验算法一定会得到不断发展，被越来越多的人所使用。

    关注 格点气象