0%

SpatialVx-package: Spatial Forecast Verification

过滤方法

  1. 领域方法:

    邻域方法通常将卷积核平滑器应用于验证集中的一个或两个字段,然后应用传统分数。 Ebert (2008, 2009)

  2. 尺度分离方法:

    带通滤波器:允许特定频段的波通过同时屏蔽其他频段的设备。

    尺度分离是指将带通滤波器(和/或进行多分辨率分析,MRA)应用于验证集的想法。通常,技能是按比例评估的。然而,也应用了其他技术。例如,在应用传统统计之前对字段进行去噪,使用变异函数,或应用基于变异函数的统计测试(最后这些与尺度分离思想的精神不太相似,但至少有些相关)。
    Briggs和Levine (1997)中提出的小波方法具有一定的功能性。特别是,要在应用传统的验证统计数据之前简单地对字段进行降噪,
    Casati等人(2004年)介绍的强度标度技术和Casati (2010年)提出的技术的新发展
    虽然不是严格意义上的“尺度分离”方法,但结构函数(变异函数是其特例)在分析不同分离距离的场的意义上是相同的,并且这些“尺度”是彼此分离的(即分数不一定随着尺度的增加而提高或降低)。
    也有对这些函数的微小修改(对fields函数的微小修改)来计算Harris等人(2001)的结构函数。

移位方法

在 Gilleland 等人 (2009) 中,该类别分为两种主要类型,即场变形和基于特征。前者将二进制图像测量/度量与场变形技术集中在一起,因为二值图像测量告知了两个场(跨整个场)的空间范围或模式之间的“相似性”(或不相似性)。在这里,它们被进一步分解为仅产生单个(或小向量)度量或度量(位置度量)的那些,以及具有移动网格点位置以更好地匹配字段的机制的那些空间上(场变形)。

  1. 基于距离和空间对齐的汇总措施:

Gilleland (2020) 引入了一种新的空间对齐汇总度量,该度量介于零和一之间,一个代表完美匹配,零代表不匹配。有一个用户可选择的参数/参数确定测量值向零下降的速率。另外两个汇总度量也包含强度信息。
除了上述新测度之外,还包括较早的知名测度,包括:Hausdorff 度量、partial-Hausdorff 度量、FQI (Venugopal et al., 2005)、Baddeley’s delta metric (Baddeley, 1992; Gilleland, 2011; Schwedler et al., 2011), metrV (Zhu et al., 2011),以及在 Baddeley, 1992 中描述的定位性能测量:平均误差距离、均方误差距离和 Pratt’s Figure of
Merit (FOM)。
图像矩可以提供有关位置误差的有用信息,并在基于特征的方法中使用,特别是 MODE,因为它们给出了图像(或特征)的质心,以及方向角,以及其他有用的属性。

  1. 场变形

这些函数执行 Marzban 和 Sandgathe (2010) 中描述的分析,并且基于 Lucas 和 Kanade (1981) 的工作。
可以使用rigider 函数来估计刚性变换。要使用指定的参数(x 和 y 平移和/或旋转)简单地对场(或特征)进行刚性变换,可以使用刚性变换函数。对于这些可能会导致变换不能完全落在网格点上的函数,函数 Fint2d 可用于从最近的网格点进行插值

  1. 基于特征的方法:

这些方法有时也称为基于对象的方法(此包中使用术语“特征”以区别于 R 对象),并且与基于对象的图像分析 (OBIA) 中使用的技术有许多相似之处,这是一种相对主要由于地球观测传感器和 GIS 科学的进步而出现的新研究领域(Blaschke 等人,2008 年)。它试图识别一个字段中的各个特征,然后逐个特征地分析这些字段。除了位置特定的错误信息之外,这可能还涉及强度错误信息。此外,可以使用命中、未命中和误报的新定义找到列联表验证统计信息(正确的否定更难以评估,但也可以完成)。

目前,有执行 Davis 等人介绍的分析的功能。 (2006,2009),包括 Gilleland 等人 (2008) 的合并/匹配算法,以及 Wernli 等人 (2008, 2009) 的 SAL 技术。复合分析的一些功能(Nachamkin,2004)是通过将单个特征放置到相对网格上来提供的,以便每个特征共享相同的质心。形状分析部分由识别边界点的函数支持(Micheas 等人 2007;Lack 等人 2010)。
增加了 Marzban 和 Sandgathe(2006;2008)的聚类分析方法。

  1. 几何特征的措施:
    也许这个子标题中的措施最好被描述为 2c 的一部分。它们当然在该领域很有用,但 AghaKouchak 等人也为整个领域提出了建议。 (2011);尽管在 MODE 等中已经应用了类似的措施。 AghaKouchak 等人介绍的措施。 (2011) 可在此处获得:连通性指数 (Cindex)、形状指数 (Sindex) 和面积指数 (Aindex):

空间(和/或时空)场的统计推断:

除了Gilleland et al.(2009)分类的方法外,还有用于两个空间域之间比较的函数。Elmore等人(2006)详细介绍的场显著性方法也需要时间维度,它涉及在每个网格点(或位置)单独使用圆形块引导(CBB)算法(通常针对平均误差),以确定网格点的显著性(零假设的平均误差为零),然后用半参数蒙特卡罗方法即Livezey和Chen(1983)确定场的意义
Hering和Genton(2011)提出的空间预测比较检验(spatial prediction comparison test, SPCT)。计算损失函数的支持函数包括:绝对误差(abserrloss)、平方误差(sqerrloss)和相关技能(corrskill),以及Gilleland(2013)中引入的距离图损失函数(distmaploss)。

其他

Mesinger(2008)中引入的偏置校正TS和ETS(或TS dHdA和ETS dHdA)
在Lakshmanan和Kain(2010)中引入的二维高斯混合模型(GMM)方法

通过S1和ACC函数可以得到S1评分和ACC异常相关性(anomaly correlation,简称ACC)。
有关这些统计数据的更多信息,请参阅Brown等人(2012) 和Thompson和Carter(1972)。
还包括Willmott等人(2007)的地理框图函数。

github链接

vgg

transforms 运行机制

常用的图像预处理方法,提高泛化能力

1
2
3
4
5
6
train_transform = transforms.Compose([
transforms.Resize((32, 32)), # 缩放为32*32
transforms.RandomCrop(32, padding=4), # 随机裁剪
transforms.ToTensor(), # 图片转张量,同时归一化0-255 --> 0-1
transforms.Normalize(norm_mean, norm_std), # 标准化均值为0标准差为1
])

torch.load_state_dict()函数

将预训练的参数权重加载到新的模型之中

model.eval()

参考链接

model.eval()与 model.train()的区别
  • model.train()是保证BN层能够用到每一批数据的均值和方差。

    • 对于Dropout,model.train()是随机取一部分网络连接来训练更新参数。
  • model.eval()是保证BN层能够用全部训练数据的均值和方差,即测试过程中要保证BN层的均值和方差不变。

    • 对于Dropout,model.eval()是利用到了所有网络连接,即不进行随机舍弃神经元。
model.eval()与 torch.no_grad()的区别
  • model.eval()足够得到正确的validation/test的结果;
  • torch.no_grad()则是更进一步加速和节省gpu空间(因为不用计算和存储梯度),从而可以更快计算,也可以跑更大的batch来测试。

tqdm

一个快速,可扩展的Python进度条,可以在 Python 长循环中添加一个进度提示信息,用户只需要封装任意的迭代器 tqdm(iterator)。

跨时空尺度的预测验证方法

问题

什么是预测验证?

  • 预测:对未来状态(天气、股票市场价格等)的预测
  • 预测验证评估预测质量的过程。

比较或验证 -> 预测 & 实际发生的相应观察结果(或对真实结果的某种良好估计)

验证结果

  • 定性的(“它看起来对吗?”)
  • 定量的(“它有多准确?”)。

无论哪种情况,它都应该为您提供有关预测误差性质的信息。

为什么要验证?

需要验证预测的三个最重要的原因是:

  1. 监控预测质量 - 预测的准确性如何?随着时间的推移,它们是否会有所改善?
  2. 提高预测质量 - 变得更好的第一步是发现你做错了什么。
  3. 比较不同预测系统的质量 - 一个预测系统在多大程度上比另一个预测系统提供更好的预测,以及该系统在哪些方面更好?

预测和验证的种类

有许多类型的预测,每一种对应略有不同的验证方法
下表列出了区分预测的一种方法,以及适用于该类型预测的验证方法。 David Stephenson 提出了一种预测分类方案。 通常可以通过重新排列分类对数据进行阈值处理来将一种类型的预测转换为另一种类型的预测。

预测性质 例子 验证方法
确定的 定量降水预报 视觉、二分类、多分类、连续、空间
概率的 降水概率、集合预报 视觉、概率、集合
定性 5天展望 视觉、二法类、多分类
时空域 例子 验证方法
事件序列 一个城市的每日最高气温预报 视觉、二分类、连续、概率
空间分布 位势高度图、降雨量图 视觉、二分类、多分类、连续、概率、空间、集合
汇集空间和时间 全球月平均异常气温 二分类、多分类、连续、概率、集合
预测的特异性 例子 验证方法
二分类(是/否) 雾的形成 视觉、二分类、概率、空间、集合
多分类 寒冷、正常或温暖的条件 视觉、多分类、概率、空间、集合
连续 最高温度 视觉、连续、概率、空间、集合
面向对象或面向事件 热带气旋运动和强度 视觉、二分类、多分类、连续、概率、空间

特异性(specificity)

本质含义是“Connected with one particular thing only”即只与唯一的特定事物相关,具有专一性。在免疫学上我们会说“抗体具有特异性”,指的就是抗体具有专一性,某一特定的抗体只能与唯一一种抗原结合。其实准确来讲是一类具有特定抗原表位的抗原,因为某些不同的抗原具有相同的抗原表位(称为交叉抗原),也可与同一种抗体特异性结合。

怎样使预测变好?

预测验证领域的先驱 Allan Murphy 写了一篇关于什么是预测“好”的文章(What Is a Good Forecast An Essay on the Nature of Goodness in Weather Forecasting)。他区分了三种类型的“好”:

  • 一致性 - 基于预测者的知识库,预测的程度对应预测者对情况的最佳判断
  • 质量     - 预测与实际发生的情况相对应的程度
  • 价值     - 预测帮助决策者实现一些经济 和/或 其他利益增量的程度

由于我们对预测验证感兴趣,让我们更深入地了解预测质量。Murphy 描述了有助于预测质量的九个方面(称为“属性”)。 这些是:

  1. 偏差     - 平均预测平均观测之间的对应关系。
  2. 关联     - 预测和观测之间的线性关系的强度(例如,相关系数衡量这种线性关系)
  3. 准确性 - 预测与事实之间的一致性水平(以观测为代表)。预测和观测之间的差异就是误差。误差越小,准确度越高。
  4. 技能     - 预测相对于某些参考预测的相对准确性参考预报通常是没有技能的预报,例如随机机会持续性(定义为最近的一组观察结果,“持续性”意味着条件没有变化)或气候学。技能是指纯粹由于预测系统的“智能”而提高的准确性。天气预报可能更准确仅仅因为天气更容易预测 —— 技能会考虑到这一点。
  5. 可靠性 - 预测值和观测值之间的平均一致性。如果将所有预测一起考虑,则整体可靠性偏差相同。如果将预测分为不同的范围或类别,则可靠性与条件偏差相同,即每个类别具有不同的值。
  6. 分辨率 - 预测将事件集分类或分解为具有不同频率分布子集的能力。这意味着预测“A”时的结果分布与预测“B”时的结果分布不同。即使预测是错误的,如果预测系统能够成功地将一种结果与另一种结果区分开来,它也有分辨率。
  7. 锐度     - 预测极值的趋势。举个反例,“气候学”的预测没有锐度。锐度只是预测的属性,与分辨率一样,预测即使是错误的也可以具有此属性(在这种情况下,它的可靠性会很差)。
  8. 辨别力 - 预测区分观测结果的能力,即每当结果发生时,对结果具有更高的预测频率。
  9. 不确定性 - 观察的可变性。不确定性越大,预测就越困难。

传统上,预测验证强调准确性技能。 需要注意的是,其他预测属性的效果也对预测结果有很大影响。

预测质量 vs 预测值

预测质量预测值不同。如果预测根据某些客观或主观标准很好地预测了观测条件,则预测具有高质量。如果它能帮助用户做出更好的决定,它就有价值。

想象这样一种情况,高分辨率数值天气预报模型 预测特定地区孤立雷暴的形成,并且在该地区确实观测到了雷暴,但在模型建议的特定地点没有观察到。根据大多数标准验证措施,该预报质量较差,但对于预报员发布公共天气预报可能非常有价值。

一个高质量但价值不大的预测示例是对旱季撒哈拉沙漠上空晴朗天空的预测。

错过事件的成本很高时,故意过度预测罕见事件可能是合理的,即使也可能导致大量错误警报。

这种情况的一个例子是机场出现雾。在这种情况下,二次评分规则(那些涉及平方误差的规则)将倾向于严厉惩罚此类预测,并且诸如“命中率”之类的正向评分可能更有用。

Katz 和 Murphy (1997)Thornes 和 Stephenson (2001) 以及 Wilks (2001) 描述了评估天气预报价值的方法。相对值图有时用作验证诊断。

什么是“事实”

我们用来验证预测的“真实”数据通常来自观测数据。这些可能是雨量计测量、温度观测、卫星衍生的云量、位势高度分析等。

在许多情况下,很难知道确切的真相,因为观测存在错误不确定性的来源包括测量本身的随机误差偏差误差抽样误差和其他代表性误差,以及在分析或以其他方式更改观测数据以匹配预测规模时的分析误差

无论对错,大多数时候我们都忽略观测数据中的错误。如果观测中的误差 远小于 预测中的预期误差高信噪比),我们可以避免这种情况。在比较不同的预测方法时,即使是偏斜或抽样不足的验证数据也可以让我们很好地了解哪些预测产品比其他产品更好。解释当前正在研究的验证数据中的错误的方法。

验证结果的有效性

验证数据的数量和质量越高,验证结果自然越可信。在验证结果本身设置一些误差范围总是一个好主意。尤其重要的是(a)对于样本容量通常很小的罕见事件,(b)当数据显示出很大的可变性,以及(c)当你想知道一种预测产品是否比另一种得多(在统计学意义上)。

通常的方法是使用分析近似引导方法(取决于分数)确定验证分数的置信区间。这方面的一些很好的气象参考资料有Seaman et al.(1996)、Wilks(2011,第5章)、Hamill(1999)以及Kane和Brown(2000)。

汇集和分层的结果

为了获得可靠的验证统计数据,可以将大量的预测/观测对pairs(样本)按时间和/或空间进行汇总。样本数越大,验证结果越可靠。您还可以通过聚合较长一段时间内的验证统计信息来获得汇总结果,但要小心正确地处理非线性分数

然而,集中样本的危险在于,当数据不均匀时,它可能会掩盖预测性能的变化。它可以将结果偏向常见的采样情况(例如,站密度较高的地区,或没有恶劣天气的日子)。非同质样本可能导致使用一些常用指标高估预测技能——Hamill和Juras(2005)提供了一些明确的例子,说明这是如何发生的。

将样本分层为准均匀子集(按季节、地理区域、观测强度等)有助于梳理出特定地区的预测行为。当这样做时,请确保子集包含足够的样本,以给出值得信赖的验证结果。

方法

标准验证方法

“眼球”验证

最古老和最好的验证方法之一是老式的视觉或“眼球”方法:

将预测和观测放在一起,用人类的判断来辨别预测的错误。表示数据的常用方法是时间序列地图

http://shuwenlovestudy.top/Forecast_Verification_methods_Across_Time_and_Space_Scales/timeseries.gif
http://shuwenlovestudy.top/Forecast_Verification_methods_Across_Time_and_Space_Scales/DWDmaps.gif

如果你只有几个预测,或者你有很多时间,或者你对定量验证统计不感兴趣,那么眼球法是很好的。即使当您需要统计数据时,不时地查看数据也是一个非常好的主意!

然而,眼球法不定量的,它很容易产生解释的个体、主观偏见。因此,在任何正式的核查程序中都必须谨慎使用。

下面几节相当简短地描述了二分类、多分类、连续和概率预测的标准验证方法和评分。有关标准方法的更多细节和讨论,请参阅Stanski等人(1989)或一本关于预测验证和统计的优秀书籍。

二分(是/否)预测法

二分预测说,“是的,一个事件会发生”,或者“不,这个事件不会发生”。雨和雾预报是不是预报的常见例子。对于某些应用,可以指定一个阈值来区分“是”和“否”,例如,风速大于50节。

为了验证这种类型的预测,我们从一个列联表开始,该表显示了“是”和“否”预测和发生的频率。预测(是或否)和观测(是或否)的四种组合称为联合分布:

  • hit - 预测要发生的事,实际上也发生了
  • miss - 事件预测不会发生,但实际上发生了
  • false alarm - 事件预测会发生,但实际上没有发生
  • correct negative - 事件预报不发生,而实际上也没有发生

在列联表的下半部分和右半部分给出了观测和预测的发生和不发生的总数,称为边际分布。

http://shuwenlovestudy.top/Forecast_Verification_methods_Across_Time_and_Space_Scales/%E5%88%97%E8%81%94%E8%A1%A8.jpg

列联表是查看所犯错误类型的有效方法。一个完美的预测系统只会产生hit和correct negative,不会出现漏报或误报。

从列联表中的元素计算出大量种类繁多的分类统计数据,以描述预测性能的特定方面。我们将用一个(编造的)例子来说明这些统计数据。假设一年的官方每日降雨预报和观测产生了以下列联表:

http://shuwenlovestudy.top/Forecast_Verification_methods_Across_Time_and_Space_Scales/%E5%88%97%E8%81%94%E8%A1%A8_%E6%95%B0%E6%8D%AE.jpg

可以从yes/no列联表中计算的分类统计如下所示。有时这些分数以括号中显示的替代名称来表示。


精度(分数正确) - $Accuracy = \frac{hits \ + \ correct \ negatives}{total}$

  • 回答了以下问题:总体而言,这些预测有多少是正确的?
  • 取值范围:0 ~ 1。完美的分数:1。
  • 特点:简单、直观。

可能会引起误解,因为它严重受最常见类别影响,在罕见天气的情况下,通常是“无事件”。

在上面的例子中,Accuracy =(82 + 222) / 365 = 0.83,表明83%的预测是正确的。


偏差得分(频率偏差) - $BIAS = \frac{hits \ + \ false \ alarms} {hits \ + \ misses}$

  • 回答了以下问题:“是”事件的预测频率与观察到的“是”事件的观测频率相比如何?
  • 范围:0到∞。完美的分数:1。
  • 特征:测量预测事件的频率与观测事件的频率的比率。
    • 表示预测系统是否倾向于低预测(BIAS<1)或高预测(BIAS>1)事件。
    • 不测量预报与观测结果的对应程度,只测量相对频率

在上例中,BIAS =(82+38) / (82+23) = 1.14,表明雨频有轻微超预报。


检测概率(命中率) - $POD = \frac{hits}{hits \ + \ misses}$ (也记作H)

  • 回答了以下问题:观测到的“是”事件中有多少是被正确预测的?
  • 取值范围:0 ~ 1。完美的分数:1
  • 特征:
    • 命中敏感,但忽略错误警报。
    • 对气候频率非常敏感的事件。对罕见事件很有帮助。
    • 可以通过发布更多“是”的预测来人为地提高点击率。
    • 应与误报率(FAR)结合使用。
    • 相对工作特性(Relative Operating Characteristic, ROC)被广泛应用于概率预测中,POD也是其重要组成部分。

在上例中,POD = 82 /(82+23) = 0.78,表明观测到的降雨事件中约有3/4的yes是正确预测的。


误报率 - $FAR = \frac{false alarms}{hits \ + \ false alarms}$

  • 回答了以下问题: 预测“是”的比例是多少,但事件实际上没有发生(即,是误报)?
  • 取值范围:0 ~ 1。完美的分数:1
  • 特点:
    • 对误报敏感,但忽略漏报。
    • 对事件的气候频率非常敏感。
    • 应该与检测概率结合使用(POD)。

在上面的例子中, FAR = 38 / (82+38) = 0.32,表明在大约 1/3 的预报降雨事件中, 没有观察到下雨。


误检概率(false alarm rate) - $POFD = \frac{false alarms}{correct negatives \ + \ false alarms}$ (也表示为 F )

  • 回答以下问题: 观察到的“否”占多大比例,但事件被错误地预测为“是”?
  • 取值范围:0 ~ 1。完美的分数:1
  • 特点:
    • 对误报敏感,但忽略漏报。
    • 可以通过发布更少的“是”预测来人为地改进以减少误报的数量。
    • 不经常报告确定性预测,但是是相对操作特征 (ROC) 广泛用于概率预测。

在上面的例子中, POFD = 38 / (222+38) = 0.15,表明对于观测到的“无雨”事件的 15% 预测不正确。


待更新。。。 太多了,而且不太重要。。。

多类别预测法

验证多类别预测的方法也从列联表开始,列联表显示各种箱子中预测和观测的频率。它类似于分类的散点图。

连续变量预测法

概率预测法

科学或诊断验证方法

空间预测法

概率预测法,包括集合预测系统

针对罕见事件的方法

其他方法

样本预测数据集

芬利龙卷风预报

降水预报概率

链接:https://cawcr.gov.au/projects/verification

未完待续。。。

翻译着翻译着,越来越看不懂,等再积累一些知识再回来翻译

验证基础

基础验证概念

什么是验证?

  1. 验证是将预测与相关观测结果进行比较的过程
  2. 验证是衡量预测的质量(相对其价值)
  3. (验证)对于许多方面,一个更合适的术语是“评价

为什么要验证?

  1. 监控模型性能
  2. 识别模型缺陷 -> 帮助操作预测员了解模型的偏差
  3. 改进决策、改善预测 -> 在不同条件下使用对应合适的模型
  4. 选择模型或模型配置(模型是否有所改进?)
  5. 纠正模型缺陷
  6. 确定预测的劣势、优势、差异 -> 帮助用户理解(解释)预测

确定验证目标

思考如下问题:

  1. 该模型在哪些位置的性能最好
  2. 是否可以调节使得预测变得更好/更坏
  3. 概率预测是否得到了良好的校准(即是否可靠)?
  4. 天气预报是否正确地捕捉到了天气的自然变化

应该衡量哪些预测性能属性?

预测”好”

取决于两个方面:

  1. 预测的质量

  2. 用户及其对预测信息的应用

    例子:

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/F_O.png

    许多验证方法结论表示,这种预测没有任何技能,而且非常不准确。

    如果我是这个流域的水务经理,这是一个相当糟糕的预测。

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/F_O_%E6%B0%B4%E5%9F%9F.png

    但如果我是一名航空交通战略规划者。。这可能是一个很好的预测

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/F_O_%E8%88%AA%E7%A9%BA.png

    1. 不同的用户对预测的好坏有不同的看法

    2. 不同的验证方法可以衡量不同类型的“好”

  3. 预测质量只是预测“好坏”的一个方面

  4. 预测值与预测质量相关(但预测质量提高,某些方面的预测值可能下降)

开展验证研究的基本指南

  1. 考虑用户(谁对验证结果感兴趣)
  2. 用户最关心的是预测质量的哪些方面
  3. 确定代表正在被预测的事件的观测结果
  4. 确定可以提供感兴趣的问题的答案的多个验证属性
  5. 选择适当度量和表示感兴趣的属性的度量图形
  6. 确定一个提供技能参考水平的比较标准(例如:持久性、气候学、旧模型)

预测和观测

预测、观测的类型

连续的

  1. 温度
  2. 雨量
  3. 500mb高度

分类

  1. 二分类
    1. 雨与无雨
      1. 强风与无强风
      2. 夜间霜冻与无霜冻
  2. 多分类
    1. 云量类别
    2. 降水类型(小雨、中雨、暴雨等)

匹配预测和观测结果

  1. 点对点网格

    1. 匹配obs到最近的网格点

    2. http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E7%82%B9%E5%AF%B9%E7%82%B9.png

      http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E7%82%B9%E5%AF%B9%E7%82%B9%E4%BE%8B%E5%AD%90.png

      匹配雨量为最近的网格点的值,最近的点(右上角)值为0,所以fcst = 0

  2. 网格点对点

    1. 插值?

    2. 取最大的价值吗?

    3. http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E7%BD%91%E6%A0%BC%E5%AF%B9%E7%82%B9.png

      http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E7%BD%91%E6%A0%BC%E5%AF%B9%E7%82%B9%E4%BE%8B%E5%AD%90.png

      将网格值插值到雨量的位置(粗略假设:每个网格点的权重相等)

      20 * 0.25 *3 = 15

不建议使用模型分析作为验证的“观测”

因为:不独立!!

非独立会影响什么?“更好的”的分数……(不具有代表性)

观测特征及其影响

观测并不完美!

观测误差 vs 可预测性和预测误差/不确定性

相同参数的不同观测类型(手动或自动)可能会影响结果

典型的仪表错误有:

  • 温度:+/- 0.1℃
  • 风速:速度相关误差,+/- 0.5m/s
  • 降水(仪表):+/- 0.1mm(half tip),但高达50%

其他问题:定位问题(例如,屏蔽/暴露)

在某些情况下,“预测”误差与仪器的限制非常相似

观察误差的影响

观测误差增加了验证结果的不确定性

对验证结果的影响

  1. RMSE - 高估
  2. 传播 - 更多 obs 异常值使整体看起来分散不足
  3. 可靠性 - 较差
  4. 分辨率 - BS 分解更大,但 ROC 区域更差
  5. CRPS - 较差的平均值

更多的样本可以有所帮助(结果的可靠性)

供验证的统计依据

验证的统计依据

任何验证活动都应从彻底检查预测和观测结果的统计特性开始。

  1. 例如,许多工具都是基于正态性(高斯分布)的假设。这是否适用于有问题的数据集?
  2. 预测是否捕捉到了观测到的范围
  3. 预测和观察到的分布是否匹配/一致?
  4. 他们有相同的平均值,变化特征等吗?

除了需要评估数据的特征之外

联合分布边际分布条件分布有助于理解预测验证的统计基础

  1. 这些分布可以与验证中使用的具体总结性能度量有关
  2. 对验证感兴趣的特定属性是由这些分布来衡量的

基本概率

又称(边际概率)

$p_x = p(X=x)$

一个随机变量X将取值x的概率

联合概率

$p_{x,y} = p(X=x,Y=y)$

事件x和y同时发生的概率

条件概率

$p_{x,y} = p(X=x|Y=y)$

给定事件y为真(或发生)时,事件x为真(或发生)的概率


验证可以表示为评价预测和观测的联合分布的过程

(验证就是预测和观测场一样,也就是两者同时发生的情况,也就是联合分布)

  1. 所有关于预测、观测及其关系的信息都用这种分布来表示
  2. 联合分布可以被分解成两对条件分布和边际分布
    • $p(f,x) = p(F=f|X=x)p(X=x)$
    • $p(f,x) = p(X=x|F=f)p(F=f)$

分布的图形化表示

联合分布

  1. 散点图
  2. 密度图
  3. 三维直方图
  4. 等高线图

边际分布

  1. 茎叶图 http://blog.sina.com.cn/s/blog_4da7fafa0100y9i3.html

  2. 直方图/柱状图

  3. 方框图

  4. 累积分布图

  5. 分位数-分位数图/Q-Q图

  6. 密度函数图(函数做y,取值做x)

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E5%AF%86%E5%BA%A6%E5%87%BD%E6%95%B0%E5%9B%BE.png

  7. 累积分布图(累计值做y,取值做x)

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E7%B4%AF%E7%A7%AF%E5%88%86%E5%B8%83%E5%9B%BE.png

条件分布

  1. 条件分位数图

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E6%9D%A1%E4%BB%B6%E5%88%86%E4%BD%8D%E6%95%B0%E5%9B%BE.png

  2. 条件箱形图

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E6%9D%A1%E4%BB%B6%E7%AE%B1%E5%BD%A2%E5%9B%BE.png

  3. 茎叶图

    http://shuwenlovestudy.top/7th_International_Verification_Methods_Workshop_Tutorial_Talks/%E8%8C%8E%E5%8F%B6%E5%9B%BE.png

比较和推理

技能分数(Skill scores)

  1. 技能分数是对相对表现的一种衡量标准
  2. 测量超过标准的改进百分比
  3. 正向导向(越大越好)
  4. 标准的选择很重要

通用技能分数定义:

$\frac {M - M_{ref}}{M_{perf} - M_{ref}}$

M是对预测的验证度量

$M_{ref}$是衡量参考预测的指标

$M_{perf}$是对完美预测的衡量标准

参考类型

类型 例子 特性
Random
随机
Equitable Threat Score
(公平的威胁得分)
很好理解的统计基准
没有物理意义
Persistence
持续
Constructed skill score
构建的技能分数
可预测性的度量(当持久性是一个糟糕的预测时,可预测性很低)
显示运行 NWP 模型的附加值
Sample climate
样品气候
Constructed skill score
构建的技能分数
比持久性更进一步,即平滑
由于政权依赖性,保留了可预测性元素
Long-term climatology
长期气候学
Constructed skill score,extremes
构建的技能得分,极端值
最简单的节拍参考,最流畅
对代表性、汇集问题、气候变化趋势的关注

应该尽可能地估计分数和测量方法的不确定性!

不确定性来自于

  • 抽样变化性
  • 观察误差
  • 代表性差异
  • 其他?

置信区间和假设检验的方法

  • 参数化(即,取决于一个统计模型)

  • 非参数(例如,源自重采样过程,通常称为“引导”)

验证属性

验证属性度量预测质量的不同方面

  • 表示应该考虑的一系列特性
  • 许多可能与预测和观测的联合、条件和边际分布有关

例子

  1. 偏见
    (边际分布)
  2. 相关性
    整体关联(联合分布)
  3. 准确度
    差异(联合分布)
  4. 校准
    测量条件偏差(条件分布)
  5. 区别
    预测区分不同观测值的程度(条件分布)

验证措施的理想特点

统计有效性

特性(概率预测)

  1. 当预测与预测者的最佳判断相一致时,就会达到“最佳”得分
  2. “对冲”是惩罚
  3. 示例:Brier分数

公平性

  1. 恒定和随机预测应获得相同的分数
    示例:Gilbert 技能得分(2x2 案例); Gerrity 评分
  2. 没有分数达到更严格意义上的这一点
    例如:大多数分数对偏差、事件频率很敏感

总结

  1. 所有的预测都应该经过验证——如果有什么东西值得预测,那就值得验证

  2. 分层和聚合

    1. 聚合可以帮助增加样本大小和统计健壮性,但也可以隐藏性能的重要方面
      1. 最常见的制度可能主导结果,掩盖了性能的变化
      2. 因此,将结果分层为有意义的同质子组是非常重要的
  3. 观测场

    1. 没有所谓的“真相”!!

    2. 观察结果通常比模型分析更“真实”(至少它们相对更独立)

    3. 无论以任何可能的方式,都应该考虑到观测结果的不确定性

      例如,相邻的观测结果彼此匹配得如何?

  4. 4w+h

    1. who:谁想知道
    2. what
      1. 用户关心吗
      2. 我们在评估什么参数吗?它的特征是什么(例如,连续的、概率性的)?
      3. 阈值很重要(如果有的话)?
      4. 预测分辨率是否相关(例如,特定地点、区域平均)?
      5. obs的特征(例如质量、不确定性)?
      6. 是否有适当的方法?
    3. why:我们需要验证吗?
    4. how:您是否需要/显示相关结果(例如,分层/聚合)?
    5. which
      1. 方法和度量标准是否合适吗?
      2. 需要采用方法(例如,偏差、事件频率、样本量)
  5. 茎状图和叶图:边缘分布和条件分布

数据准备

观测数据来源

  1. 如果我们对预报有效期内的每个地点和每个时间点进行观测,这不是很好吗?

    这样我们就可以对任何预测进行完全的验证

  2. 观测结果代表了大气在空间和时间上的真实状态的一个“样本”

  3. 观测也可能在某个点一个区域上有效

    实地观测或遥感的

  4. 实地观测 - 地面高空大气

    1. 在现场,点有效
    2. 高分辨率,但在空间中采样严重不足
    3. 的仪器几乎可以连续地取样
    4. 唯一重要的误差是仪器误差,通常很

遥感观察

卫星雷达是最常见的

  1. 雷达
    1. 测量地表上方体积内水凝物的反向散射
    2. 与感测体积中的降雨率的关系是一个复杂的函数,但已知
    3. 感知到的平均雨率与雨率(或地表的总降雨量)之间的联系要脆弱得多
    4. 误差的几种来源:衰减异常传播接近冰点的明亮频带等。
  2. 卫星
    1. 根据仪器测量一个或多个频带内的后向散射辐射。
    2. 通常低垂直分辨率 - 可以测量总柱水分
    3. 传递函数需要将返回值转换为感兴趣的变量的估计数
    4. 最有用,特别是与表面观测相结合

遥感数据

  1. 大数据
  2. 检测到的变量通常不是要验证的变量 – 需要传递函数误差源之一
  3. 分辨率取决于仪器,雷达为几米,卫星数据为1公里左右。
  4. 高覆盖率,在时间上可能是零星的
  5. 注意由于外部影响信号而造成的错误

数据特征总结

实地考察 雷达 卫星
分辨率 - 空间 高 - 点 相当高 - 雷达量平均数 取决于足迹在1公里左右
分辨率 - 时间
空间采样频率 这是很低的,除了特殊的网络 高 - 基本上是连续的 其域内的地理位置高 极地轨道变量
时间采样频率 可以是高 高,通常为10min左右 中等用于地球轨道。
低,适用于极地轨道飞行

分辨率:定义观测值的时间或空间距离

采样频率(粒度):在时间或空间上的观测频率

误差和不确定性的来源

  1. 在频率或值上的偏差

  2. 仪器偏差

  3. 随机误差或噪声

  4. 报告错误

  5. 主观观察

    例如云覆盖

  6. 精度误差

  7. 传输函数错误

  8. 分析误差

    当使用分析时

观测的质量控制

做这件事是绝对必要的,即使是“好的” 观测站点

基本方法

 1. **伙伴检查**(空间和时间)
 2. 简单的**范围检查**
 3. **趋势检查**(与附近的独立时空检查)
 4. 绝对值检查。
 5. 在不消除太多“好”数据的情况下删除“坏”数据
 6. 但不做预测 - 观测的差异检查

使用模型作为观察的比较标准不是一个好主意,作为一个过滤器来消除模型无法解决的极端情况

​ 1. 使观测数据依赖于模型

 2. 在qc中使用的模型得到了更好的验证结果

了解有关仪器及其错误的细节也很重要。

预测有效性类型

用于客观验证

“预测必须被陈述,以便它们可核实”

一个预测的意义是什么?精确度吗?

  1. 需要进行客观验证
  2. 如果验证是面向用户的,那么用户的理解是很重要的
  3. 所有预测对空间点或区域有效
  4. 在这个地区的所有点上吗?

同样对于时间:一个预测可能是

  1. 时间的一瞬间
  2. 一个瞬间在时间上,但“某个时间”在一个范围内
  3. 一段时间,例如24小时折旧
  4. 在一段时间内的极端情况吗?

预测数据来源,以供验证

所有类型的NWP模型

  1. 主要变量(P或Z、T、U、V、RH或Td)的确定性预测,通常是在模型的三维域上的网格点
  2. 其他衍生变量:由模型计算出的折旧率、折旧总额、云量和高度等,可能无法观测到
  3. 空间和时间表示被认为是连续的,但有限的尺度集可以被解决。

后处理模型输出

  1. 统计方法,如MOS
  2. 动态的或经验性的方法,例如折旧类型
  3. 相互依赖的模型,如海浪

操作预测

  1. 格式取决于用户的需要
  2. 可以是点,可以是一个区域或一段时间内的最大平均值或次要平均值

“一切都应该得到核实”

变量类型

  1. 连续的

    1. 可以在其范围内承担任何值(接近)
    2. 例如温度、风
    3. 预测是针对特定的值进行的
  2. 分类的

    1. 只能接受一小部分特定的值
    2. 可以这样观察到,如降水、降水类型、视觉障碍
    3. 可以从一个连续的变量中“分类”,例如降水量、上限、vis、云量
    4. 如果有,验证为分类发生概率
  3. 概率分布

    1. 验证为概率分布函数累积分布函数
  4. 转换变量

    1. 数值已从原始观测值中有所改变

    2. 示例:

      1. 准连续变量的分类,例如云量

      2. 要根据用户需求进行评估:

        1. “升级”到模型网格盒
        2. 插补
      3. 转换观测值的分布:

        例如,通过子设置来选择极端情况

数据匹配问题

例如,预测可以在空间上定义为“威胁区域”,或在网格(模型)上表示

  1. 有限尺度集
  2. 在空间和时间上相关

观测结果是分散的点值

  1. 代表所有刻度,但仅在车站有效
  2. 采样不足作为场

预测至观测技术:

  1. 问:在验证地点的预测是什么?
  2. 推荐的验证方法-不要考虑观测值。
  3. 插值到观测位置的插值-为平滑变量
  4. 最近的网格点-用于“情景性”或空间分类变量
  5. 除QC外,观察结果保持不变
  6. 有时,通过将模型预测转换为“如果预测正确,卫星将看到什么”,可以对遥感数据进行验证

观测预测技术(适用针对建模人员):

  1. 放大——在网格上求平均值——仅当这是预测(模型)的真正定义时,例如Cherubini等人2002
    1. 本地验证
  2. 只验证那些有数据的地方!

例子

降水验证项目:方法论 - 欧洲

升尺度:

1x1个网格框,模型分辨率的限制

网格箱上的平均磅,每个网格箱至少9吨(欧洲数据)

网格盒上的平均 obs,每个网格盒至少 9 个 stns(欧洲数据)

仅在有足够的数据处进行验证

回答在模型的能力范围内关于预测的质量的问题

最有可能的用户是建模者

模型技术观察:

  1. 在模型网格上的观测数据分析
    1. 经常做,但除了一些模型研究外,不是一个验证的好主意
    2. 使用模型独立的方法进行分析,如巴恩斯
    3. 使用依赖模型的方法进行分析-数据同化(验证错误!)例如,Park等人,2008年

不同的“真理”的影响

将匹配点obs与区域延迟的预测相匹配:事件是什么?

对于分类预测,我们必须清楚正在预测的“事件”

  1. 预测有效的位置或区域
  2. 对它有效的时间范围
  3. 类别定义

现在,什么被定义为正确的预测呢?

  1. 该事件被预测出来,并在该地区的任何地方被观察到吗?超过一定比例的面积?
  2. 规模考虑因素

收集数据进行验证

存档预测和观测

  1. 你自己的:观测站观测和相应的预报
  2. 大多数NWP中心将他们的预测和观察结果存档;如果你使用他们的模型,你可能会让他们给你相关的数据进行验证

目标:生成一套相匹配的预测和观测结果


emm, ppt看的有点迷呀,,,唉

文件下载

链接:https://www.7thverificationworkshop.de/Tutorial/Tutorial-Talks/index.html

  1. 验证基础
  2. 数据准备

思想

问题

http://shuwenlovestudy.top/FSS_FractionScoreSkill%20/%E5%BC%95%E5%85%A5%E9%97%AE%E9%A2%98.png

对于如上预测场和观测场,虽然降水都分布在7个降水格点上,但位置均未能一一对应,所以其TS评分等评分为0(没有降水预报技巧),但从概率角度看,观测与预报的降水发生概率是一样的,7/49,具有相同的降水预报面积

FSS

作者:Robert、Lean(2008)

作用:考察不同尺度内降水发生概率

具体公式

按步骤

  1. Brier评分的变形:比较预报与观测的降水频率 FBS (Fraction Brier Score)

    公式:$FBS = \frac{1}{N} \displaystyle \sum_{N}(P_{fcst} - P_{obs})^2$

    $P_{fcst}$、$P_{obs}$:每个领域尺度内 预报与预测降水发生概率(0 - 1)

    N:网格数量

  2. 方差技巧评分:获取正定的(特征值均为正)评分

    公式:$FSS = 1 - \frac{FBS} {\frac{1}{N}(\displaystyle \sum_N P_{fcst}^2 + \displaystyle \sum_N P_{obs}^2)}$

    FSS

    1. 范围(0 - 1)

      0 :完全不匹配

      1:完美匹配

    2. 特点:
      领域尺度增加 -> FSS评分增大

      当$P_{fcst} = bP_{obs}$,FSS评分向 $2b/(b^2+1)$ 渐近

简写

$FSS = 1 - \frac{\frac{1}{N} \displaystyle \sum_{N}(P_{fcst} - P_{obs})^2} {\frac{1}{N}(\displaystyle \sum_N P_{fcst}^2 + \displaystyle \sum_N P_{obs}^2)}$

代码

R语言

代码

链接:https://rdrr.io/cran/verification/src/R/fss.R#heading-0

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
fss <- function(obs, pred, w = 0, FUN = mean, ...){
### compare matrixes of forecast of observed values and forecast.
### values can be calcuated using different windows.

### with a window size of 0, obs is returned.
obs.matrix <- matrix.func(DAT = obs, w = w, FUN = FUN)

### with a window size of 0, obs is returned.
frcs.matrix <- matrix.func(DAT = pred, w = w, FUN = FUN)

if(nrow(obs)!= nrow(pred) & ncol(obs)!= nrow(obs) ) stop("Observation matrix and forecast matrix different sizes")

n <- prod(dim(obs.matrix)) ### number of gridpoints

N <- sum((obs.matrix-frcs.matrix)^2, na.rm = TRUE)/n ### numerator
D <- (sum(obs.matrix^2, na.rm = TRUE) +sum(frcs.matrix^2, na.rm = TRUE))/n ### denominator

FSS <- 1 - N/D
return(FSS)
}

matrix.func <- function(DAT, w = 0, FUN = mean, ...){

### w is the '' radius'' of window. eg. w = 2, defines a 5 by 5 square

### define function
FUN <- match.fun(FUN)

### define output dimension
II <- nrow(DAT) - 2*w ## output row dimension
JJ <- ncol(DAT) - 2*w

if(JJ<=0|II <= 0) {stop("The window exceeds the size of the observation" ) }
OUT <- matrix(NA, nrow= II, ncol = JJ)

for(i in 1:II){
for(j in 1:JJ){
sub <- DAT[ i :(i + 2*w ),
j :(j + 2*w ) ] # subset data

OUT[i,j] <- FUN(sub,...)
} ## close J
} ## close I

return(OUT)
} ## close functionr

测试

链接:https://rdrr.io/cran/verification/man/fss.html

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
library(verification)

grid<- list( x= seq( 0,5,,100), y= seq(0,5,,100))
obj<-Exp.image.cov( grid=grid, theta=.5, setup=TRUE)
look<- sim.rf( obj)
look[look < 0] <- 0

look2 <- sim.rf( obj)
look2[look2<0] <- 0

fss(look, look2, w=5)


## Not run:
# The following example replicates Figure 4 in Roberts and Lean (2008).
#### examples

LAG <- c(1,3,11,21)
box.radius <- seq(0,24,2)

OUT <- matrix(NA, nrow = length(box.radius), ncol = length(LAG) )

for(w in 1:4){

FRCS <- OBS <- matrix(0, nrow = 100, ncol = 100)

obs.id <- 50
OBS[,obs.id] <- 1

FRCS[, obs.id + LAG[w]] <- 1

for(i in 1:length(box.radius)){

OUT[i, w] <- fss(obs = OBS, pred = FRCS, w = box.radius[i] )

} ### close i
} ### close w

b <- mean(OBS) ## base rate

fss.uniform <- 0.5 + b/2
fss.random <- b

matplot(OUT, ylim = c(0,1), ylab = "FSS", xlab = "grid squares",
type = "b", lty = 1, axes = FALSE , lwd = 2)

abline(h = c(fss.uniform, fss.random), lty = 2)
axis(2)
box()
axis(1, at = 1:length(box.radius), lab = 2*box.radius + 1)
grid()

legend("bottomright", legend = LAG, col = 1:4, pch = as.character(1:4),
title = "Spatial Lag", inset = 0.02, lwd = 2 )

## End(Not run)

结果截图

http://shuwenlovestudy.top/FSS_FractionScoreSkill%20/r%E8%AF%AD%E8%A8%80%E7%BB%93%E6%9E%9C%E5%9B%BE.png

Python

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import numpy as np 

def func(data, w = 0):
row = data.shape[0]
col = data.shape[1]

ii = row - 2 * w
jj = col - 2 * w

if(ii <= 0 or jj <=0):
print("error!")

out = np.zeros((ii, jj))

for i in range(1, ii):
for j in range(1, jj):

sub = data[i : (i + 2 * w + 1), j : (j + 2 * w + 1)]

out[i, j] = np.mean(sub)
return out

def fss(obs, fc, w = 0):
obs_m = func(obs, w)
# print(obs_m[0][16])
# print(obs_m[0][17])
# print(obs_m[0][18])
fc_m = func(fc, w)
# print(obs_m)
# print(fc_m)
n = obs_m.shape[0] * obs_m.shape[1]
N = ((obs_m - fc_m) ** 2).sum() / n
# print(N.shape)
D = ((obs_m ** 2).sum() + (fc_m ** 2).sum()) / n
# print(n)
# print(N)
# print(D)
return 1 - N / D

测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import matplotlib.pyplot as plt

lag = [1, 3, 11, 21]
radius = np.arange(0, 24 + 2, 2)

out = np.zeros((len(radius), len(lag)))

frcs = np.zeros((100, 100))
obs = np.zeros((100, 100))

id = 50 - 1
obs[ : , id] = 1

for l in range(0, len(lag)):
for r in range(0, len(radius)):
frcs[ : , id + lag[l - 1]] = 0
frcs[ : , id + lag[l]] = 1
out[r][l] = fss(obs, frcs, radius[r])
b = np.mean(obs)
uniform = 0.5 + b / 2
random = b

plt.plot(out)
plt.ylabel("FSS")
plt.xlabel("grid squares")

结果截图

http://shuwenlovestudy.top/FSS_FractionScoreSkill%20/python%E7%BB%93%E6%9E%9C%E5%9B%BE.png

坐标可以自定义修改,值是一样的

http://shuwenlovestudy.top/FSS_FractionScoreSkill%20/%E7%BB%93%E6%9E%9C%E7%9F%A9%E9%98%B5.png

下载

Roberts2008原文

思想

SAL方法 通过三个分量,分别从 场结构、降水强度和场位置偏差 对预测结果进行评价

threshold

确定方法是经验法

阈值的设定,效果是大于阈值的格点才会被定为降水主体成员降水体:降水主体内(各不连续的小降水区域)每个连续的小区域),只有降水主体成员才参与上述S、L2的运算

three component

  1. 计算总雨带特征:A、L1 (与阈值无关)

  2. 计算内部结构特征:S、L2

Structure

S表示Structure 结构,提供 预测场和观测场 有关其大小形状的信息。

  1. 正值 -> 表示模拟的降水对象太大或太平
  2. 负值 -> 表示对象太尖锐或太小

http://shuwenlovestudy.top/SAL_StructureAmplitudeLocation/S%E8%AF%84%E5%88%86%E7%AE%80%E4%BB%8B.png

Amplitude

A表示Amplitude (振幅)降水强度

  1. 正值 -> 表示对总降水量的估计过高

  2. 负值 -> 表示对总降水量的估计过低

Location

L表示Location位置,提供预测场和观测场 位置偏差信息

http://shuwenlovestudy.top/SAL_StructureAmplitudeLocation/L%E8%AF%84%E5%88%86%E7%AE%80%E4%BB%8B.png

分为L1、L2两部分

  1. L1表示预报场与观测场质心之间的归一化距离。(两个场质心之间的距离)

  2. L2是总降水场的质心与观测和预报的个别降水对象的归一化距离之差。(每个降水对象到质心的距离)

    L2是为了预防下图情况,由于观测场分为两块,质心在中间,与预测场重合

    http://shuwenlovestudy.top/SAL_StructureAmplitudeLocation/L2%E8%AF%84%E5%88%86%E7%AE%80%E4%BB%8B.png

    some example

http://shuwenlovestudy.top/SAL_StructureAmplitudeLocation/SAL%E4%BE%8B%E5%AD%90.png

( S = 0、A=0、L=0表示完美匹配)

具体公式

S - 结构

公式:$S = \frac{V(R_{mod})-V(R_{obs})}{0.5[V(R_{mod})+V(R_{obs})]}$

  1. $V(R_{mod})$:预测场中所有降水体以体内总降水量为权重的$V_n$的加权平均

    1. 公式:$V(R) = \frac{\sum^{m}_{n=1}R_nV_n}{\sum^m_{n=1}R_n}$
      1. $R_n$:降水体的总降水量
  2. $V_n$:第n个降水体内的总降水量最大降水量

    😃降水量就是某一个降水体每个降水值的总和

    1. 公式:$V_n = {\displaystyle \sum_{(i,j)\in R_n} {R_{i,j}/R^{max}_n} = R_n/R^{max}_n}$

      1. $R_{i,j}$:格点(i,j)处的降水

      2. $R^{max}_{n}$:降水体n内的最大降水值

        😃降水体相当于是面积,最大值表示这块面积里面的某一个点的最大值

A - 强度

公式:$A = \frac {D(R_{mod}) - D(R_{obs})} {0.5[D(R_{mod}) + D(R_{obs})]}$

  1. $D(R) = \frac{1}{N} \displaystyle \sum_{(i,j)\in D}{R_{i,j}}$ → $D(R)$为所取区域内非缺省格点降水的平均值

L - 位置

公式

  1. $L = L_1 + L_2$

    1. $L_1 = \frac{|x(R_{mod})-x(R_{obs})|}{d}$ → L1为区域D实况与预报降水主体中心之间的距离

      1. $x(R)$:表示所取区域D内降水主体的重心位置
      2. $d$:区域D内最大距离
    2. $L_2 = 2\cdot{\frac{r(R_{mod})-r(R_{obs})}{d}}$ → L2为降水主体重心与降水场每个降水体重心之间的平均距离

      1. $r = \frac{\sum ^{m}_{n=1}{R_n \cdot{|x-x_n|}}}{\sum^{m}_{n=1}R_n}$

        1. $m$ 表示降水体个数
        2. $r(R)$ 表示各降水体重心与降水主体中心之间距离加权平均(权重是雨量)的权重
        3. 降水体雨量越大 → 离降水主体重心越远 → r值越大 → 最大值为(1/2)d
      2. $R_n = \displaystyle \sum_{(i,j)\in R_n}{R_{i,j}}$ → 区域内第n个降水体的降水量

        😃降水体是一个小连续区域,就是它的总降水量

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
import numpy as np
from scipy import ndimage
#numpy.linalg模块包含线性代数的函数。使用这个模块,可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等。
#norm则表示范数
from numpy.linalg import norm
import cv2
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(threshold=np.inf)

#obs:观测场;fc:预测场
#返回值前三个是S、A、L,还有其他信息
def compute_SAL(field_obs, field_fc, user_threshold=None):
"""SAL-score (Structure, Amplitude, Location)

Parameters
----------
field_obs : numpy.ndarray
2D observation (or reference) field
field_fc : numpy.ndarray
2D forecast field (same shape as `field_obs`)
user_threshold : float, optional
If set, use this threshold for object separation. 如果已设置,请使用此阈值进行对象分离。
If not set, use the default method of Wernli is used: ``threshold = field_obs.max() * 1/15``

Returns
-------
S : float
Structure score
A : float
Amplitude score
L : float
Location score
L1, L2 : float
Location score components `L1` and `L2`
NObjects_obs : float
Number of thresholded objects in observation 观察中的阈值对象数
NObjects_mod : float
Number of objects in model (forecast `fc`) 模型中的对象数(预测“fc”)
R_max : float
Maximum Rain-value of observations `field_obs` 观测值的最大降雨值`
user_treshold : float
Threshold that was actually used 实际使用的阈值
"""
if (field_obs.shape != field_fc.shape):
raise Exception("field_obs and field_fc need to have the same shape" +\
" (this function also assumes that they cover the same domain!)")

if (user_threshold is not None) and (type(user_threshold) != float):#user_threshold是给定的阈值
raise Exception("user_threshold must be a float")

R_mod = field_fc # "mod" is for "model" 预测场的数据
R_obs = field_obs # "obs" is for "observation" 观测场的数据

#不能出现有数据小于0的,出现即报错
if np.any(R_mod < 0.):
raise Exception("negative values in forecast." +\
"use different field or apply SAL.threshold_to_zero(field)")

if np.any(R_obs < 0.):
raise Exception("negative values in forecast." +\
"use different field or apply SAL.threshold_to_zero(field)")

f_threshold = 1/15.

#nanpercentile:不算nan值找到第95个百分数
#https://www.cjavapy.com/article/1088/
R_max_obs = np.nanpercentile(np.where(R_obs>0.1, R_obs, np.nan), 95)#最大*0.95
threshold_obs = R_max_obs * f_threshold#1/15

R_max_mod = np.nanpercentile(np.where(R_mod>0.1, R_mod,np.nan), 95)
threshold_mod = R_max_mod * f_threshold

if (user_threshold is not None):
threshold_obs = user_threshold
threshold_mod = user_threshold

R_mod_thr = np.where(R_mod > threshold_mod, R_mod, 0.) # 所有的降水体
R_obs_thr = np.where(R_obs > threshold_obs, R_obs, 0.) # maskiert

D_R_mod = np.mean(R_mod)#平均值
D_R_obs = np.mean(R_obs)

A = (D_R_mod - D_R_obs) / (0.5 * (D_R_mod + D_R_obs))#与阈值无关

# x_R_mod: center of mass of mod
x_R_mod = np.array(ndimage.measurements.center_of_mass(R_mod))#nb,直接调用函数找出质点
# x_R_obs: center of mass of obs
x_R_obs = np.array(ndimage.measurements.center_of_mass(R_obs))

#norm默认是2范式,每个元素先平方,再根号,也就是对角线距离
d_diagonal = norm(#field_obs.shape)

L1 = norm(x_R_mod - x_R_obs) / d_diagonal#与阈值无关

##label标记连接成分
#https://blog.csdn.net/Monkey_dada/article/details/119353323?spm=1001.2101.3001.6661.1&utm_medium=distribute.pc_relevant_t0.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1.no_search_link&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1.no_search_link
labels_mod, NObjects_mod = ndimage.measurements.label(R_mod_thr)#label标记连接成分,在找降水体
labels_obs, NObjects_obs = ndimage.measurements.label(R_obs_thr)

labels_list = [labels_mod, labels_obs]
NObjects_list = [NObjects_mod, NObjects_obs]

R_list = [R_mod, R_obs]
x_R_list = [x_R_mod, x_R_obs]

r_list = [0,0] # fuer L2
V_list = [0,0] # fuer S

for iList, NObjects in enumerate(NObjects_list):
R_n_list = np.zeros(NObjects)
V_n_list = np.zeros(NObjects)
distance_n_list = np.zeros(NObjects)

for n in range(NObjects):

label = n + 1

object_n = np.where(labels_list[iList] == label, R_list[iList], 0)
R_n_list[n] = object_n.sum()
V_n_list[n] = R_n_list[n] / object_n.max()
x_n = np.array(ndimage.measurements.center_of_mass(object_n))
distance_n_list[n] = norm(x_R_list[iList] - x_n)

R_n_list = np.array(R_n_list)
distance_n_list = np.array(distance_n_list)
r_list[iList] = (R_n_list*distance_n_list).sum() / R_n_list.sum()
V_list[iList] = (R_n_list*V_n_list).sum() / R_n_list.sum()

L2 = 2.*norm(r_list[0] - r_list[1]) / d_diagonal
L = L1 + L2
S = (V_list[0] - V_list[1]) / (0.5 * (V_list[0] + V_list[1])) # 0: mod, 1:obs

return S, A, L, L1, L2, NObjects_obs, NObjects_mod, R_max_obs, R_max_mod, threshold_obs, threshold_mod

#打印结果
def PrintResultSAL(field_obs, field_fc, user_threshold=None):
S, A, L, L1, L2, NObjects_obs, NObjects_mod, R_max_obs, R_max_mod, threshold_obs, threshold_mod = compute_SAL(field_obs,field_fc)
print("S = ", S)
print("A = ", A)
print("L = ", L)
# print("L1 = ", L1)
# print("L2 = ", L2)
# print("NObjects_obs = ", NObjects_obs)
# print("NObjects_mod = ", NObjects_mod)
# print("R_max_obs = ", R_max_obs)
# print("R_max_mod = ", R_max_mod)
# print("threshold_obs = ", threshold_obs)
# print("threshold_mod = ", threshold_mod)
# print()

def matrix_SAL(filename1, filename2):
#预测
EC_fst = np.load(filename1)
file_fst = EC_fst['f']
result_fst = file_fst[0].reshape(210,211)
#真值

EC_obs = np.load(filename2)
file_obs = EC_obs['l']
result_obs = file_obs[0].reshape(210,211)
PrintResultSAL(result_fst,result_obs)

#之前强降水数据集使用例子
#matrix_SAL('testset_ECrainPredi_03h_202008.npz','../CMPA/08/testset_CMPA_03h_202008.npz')

下载

  1. wernli2008原文
  2. 一篇较好的总结资料

2021年推免结束了,最终去了国防科技大学 气象海洋学院,忙里偷闲,记录一下自己的一些感受。

image-20211014221958513

个人情况

本科双非,坐标长春理工大学

成绩:2/254,前1%

英语:四六级通过(六级飘过)

比赛:网安选手,一两个国赛、n个省赛、n个校赛

项目:一个导师的在研项目

夏令营

特别后悔开始投了很多学硕,双非孩子不要想了,全投专硕!!

一定要从夏令营开始重视,不要想着开始先玩一玩,你永远不知道每年的保研形势是什么样子!!

双非本来就非常不占优势了,一定要做好每一步,你才有机会!!

投递了很多学校,但是效果甚微,一方面文书没有准备好,一方面投递学硕策略不正确。


入营:山东大学、西北工业大学、吉林大学、东南网安(无锡)、浙大软院(报名就给)、国防科大 气象海洋

offer:吉林大学、国防科大 气象海洋


emmm、山大和西工大因为地理不喜欢,也没有好好准备,没有拿到offer

吉林大学可能因为地理原因,基本都给了offer

东南今年初审直接分配了无锡,直接无了

浙大软院因为在宁波,和三本学校共用校区也无了

国防科大其实偶然机会报名了,但是比较幸运进去了,也是最终去向。

这里还打了信安国赛,参加了哈工大的暑期学校,挺好的,哈哈哈

预推免

预推免并没有比夏令营简单,哭哭,

这时候简历什么的文书已经相对完善了,但是机会不多了


入营:西工大(我又来啦)、电科-深圳、华南理工大学

offer:西工大、电科深、华工


差不多入营了都拿到offer了,感觉夏令营和预推免难度差不多,,,电科和华工投的都是专硕

结果

image-20211014224904410

  1. 电科

其实电科深的研究发现我很喜欢,导师是二进制的,太难得了,但是导师没回过我,而且电科深据说都放羊,就放弃了

电科发消息很快,感受很好,我和朋友都很快收到了录取通知

  1. 华工

简历都是网安,问的我全是计算机,网安哭哭,但是还是过了,

感觉被养🐟,和同学差了一分多,他很早收到通知了,我一直没收到,,,

  1. 西工大

方向不感兴趣,地理也不好,当然不会去了

  1. 国防科大

其实电科和华工没出来之前,难受了很久,难道自己就只能去西工大或者吉大了么(因为非常想回南方)。

然后电科深和华工出来后很开心,不过因为和国防科大老师说好了后,就决定去那里了。

国防科大,我是没有想到的,军校,曾经自己也有一段军旅梦呀,未来好好加油!

我的推免,满意也算满意,不满意也有些不满意,但总的来说,还是不错的。

如果从来一次,自己肯定能做的更好,冲个华5 网安呀,梦校

但是国防科大一方面学校title挺好的,导师也很好,push型是我爱的,学校补贴,小组氛围也很好,唯一遗憾的是专业变成大气科学了,有机会读个博士变回来哈哈哈哈~

勉励

唉,双非学校不太看重科研论文,竞赛权重更重要。

如果认准了保研,首先绩点必须稳、最好rank1;其次,英语四六级550+,保研的时候,你会发现大家都500+的,,,;然后,最好参加ACM竞赛,拿几个国家级、省级的竞赛;最后,大二左右跟随学校某实验室老师准备科研/有钱可以考虑花钱跟机构做(产出保证)。高质量论文多多益善(1-2篇)、其他的项目(最好是国家级大创),再来几个就行了。

感觉,能做到这几点,出来无敌了。