计算等压面要素场的基本检验指标

August 26, 2020 (最后修改: August 21, 2021)

meteorology/verification

本文介绍计算等压面要素场的几种基本检验指标。

重要提示：本文介绍的部分代码参考自其他检验工具。

介绍

本文介绍的基本检验指标涉及以下三种数据：

预报场 (forecast_field)
分析场 (analysis_field)：当前模式的分析场或者使用 FNL 等其他分析场
气候场 (climate_field)

本文假定上述数据均被插值成 WMO GDPFS 手册中规定的标准网格，即 1.5 度 * 1.5 度。下面的指标计算不涉及数据插值问题。

计算指标还需要使用到网格点对应的纬度坐标值 (latitudes) 。

指标计算即可以针对全球范围，也可以针对特定的区域范围。本文假定上述四个参数都已按照选定的区域范围进行裁剪。

指标

Mean Error

Mean Error (ME)，也叫作 Bias，表示预报值与验证值之间的偏差的平均值。定义如下：

\begin{align*} \textrm{ME} &\equiv \left ( \sum_{i=1}^{n}w_iD_i \right ) \mathbin{/} \sum_{i=1}^{n}w_i \ \end{align*}

\begin{align*} D_i &= F_i - A_i \ \end{align*}

\begin{align*} w_i &= \frac{1}{n} (\text{or } \cos\phi_i \text{, and so on}) \end{align*}

其中

F_i 表示预报值
A_i 表示验证值 (verification value)
D_i 表示两者的差值
w_i 是权重系数，n 是样本的个数，phi_i 是纬度。

观测值，初始值或客观分析值通常被用作验证值。当预报完全正确时，ME 等于 0，被称为完美预报 (perfect forecast)。

计算北半球 (Northern Hemisphere) 等广阔区域的平均值时，需要考虑区域之间与纬度相关的差异，并使用加权系数进行评估。

例如，为了评估等角投影 (euirectangular projection) 中的客观分析，通常将加权系数 “wi = 1/n"替换为纬度的余弦 “cos phi_i”。本文中介绍的其他指标采用相同的计算方法。

回答如下问题：平均预报误差是多少？

范围：负无穷大到正无穷大

完美分数：0

特性：

简单，熟悉。也称为（加性）偏差 (additive bias)。不测量误差的大小。不测量预测值和观察值之间的对应关系，即，如果存在补偿误差，则可以为不良预测获得完美的分数。

Mean Absolute Error

\begin{align*} \textrm{MAE} &\equiv \left ( \sum_{i=1}^{n}w_i\left | F_i - A_i \right | \right ) \mathbin{/} \sum_{i=1}^{n}w_i \ \end{align*}

回答如下问题：预报误差的平均幅度是多少？

范围：0 到正无穷大

完美分数：0

特性：

简单，熟悉。不指示偏差的方向。

Standard Deviation

Standard Deviation (SD)

\begin{equation*} SD \equiv \sqrt{\left ( \sum_{i=1}^n w_i(D_i - \textrm{ME})^2 \right ) \mathbin{/} \sum_{i=1}^n w_i} \end{equation*}

Root Mean Square Error

Root Mean Square Error (RMSE) 通常用于表示预测准确性。定义如下：

\begin{align*} \textrm{RMSE} &\equiv \sqrt{\sum_{i=1}^{n}w_iD_i^2 \mathbin{/} \sum_{i=1}^{n}w_i} \end{align*}

\begin{align*} D_i &= F_i - A_i \ \end{align*}

\begin{align*} w_i &= \frac{1}{n} (\text{or } \cos\phi_i \text{, and so on}) \end{align*}

其中 D_i 表示预报和验证值的差值，w_i 表示权重系数，n 是样本个数。

RMSE接近零表示预测值更接近验证值。对于完美预测，RMSE等于零。将 ME 和随机误差分开，RMSE 可以表示如下：

\begin{equation*} \textrm{RMSE}^2 = \textrm{ME}^2 + \sigma_e^2 \end{equation*}

其中 \sigma_e^2 表示差值 D_i 的标准差 (standard deviation, SD)

\begin{equation*} \sigma_e^2 = \left ( \sum_{i=1}^n w_i(D_i - \textrm{ME})^2 \right ) \mathbin{/} \sum_{i=1}^n w_i \end{equation*}

回答如下问题：预报误差的平均幅度是多少？

范围：0 到正无穷大

完美分数：0

特性：

简单，熟悉。测量“平均”误差，并根据误差的平方加权。不指示偏差的方向。与较小的错误相比，RMSE 对较大的错误具有更大的影响，如果特别不希望出现较大的错误，则这可能是一件好事，但也可能鼓励保守的预测。

Anomaly Correlation Coefficient

异常相关系数 (Anomaly Correlation Coefficient, ACC) 是空间场验证中使用最广泛的量度之一 (Jolliffe and Stephenson 2003)，它代表预报异常与验证值与参考值（如气候数据）之间的相关性。 ACC 定义如下：

\begin{equation*} \textrm{ACC} \equiv \frac { \sum\limits_{i=1}^{n}w_i(f_i - \bar{f})(a_i - \bar{a}) } { \sqrt{\sum\limits_{i=1}^{n} w_i(f_i - \bar{f})^2 \sum\limits_{i=1}^{n} w_i(a_i - \bar{a})^2} } \end{equation*}

其中

\begin{equation*} f_i = F_i - C_i \end{equation*}

\begin{equation*} \bar{f} = \left ( \sum_{i=1}^{n}w_if_i \right ) \mathbin{/} \sum_{i=1}^{n}w_i \end{equation*}

\begin{equation*} a_i = A_i - C_i \end{equation*}

\begin{equation*} \bar{a} = \left ( \sum_{i=1}^{n}w_ia_i \right ) \mathbin{/} \sum_{i=1}^{n}w_i \end{equation*}

F_i 表示预报值
A_i 表示验证值
C_i 参考值，例如气候值
\bar{f} 是 f_i 的均值
\bar{a} 是 a_i 的均值
w_i 表示权重系数

如果预测异常的变化模式与验证异常的变化模式完全一致，则 ACC 将采用最大值 1。相反，如果变化模式完全反转，则将采用最小值 -1。

回答如下问题：预报异常与观测异常的对应程度如何？

范围：-1 到 1

完美分数：1

特性：

测量预报和观测值之间的对应关系或相位差，减去每个点的气候平均值 C，而不是样本平均值。经常使用异常相关性来验证数值天气预报（NWP）模式的输出。 ACC 对预报偏差不敏感，因此良好的异常相关性不能保证准确的预报。

计算

使用 numpy 或 xarray 库保存数组

import numpy as np
import xarray as xr

函数参数说明：

forecast_field：预报场
analysis_field：分析场
climate_field：气候场
latitudes：纬度坐标数组

ME

def bias(
    forecast_field: np.ndarray or xr.DataArray,
    analysis_field: np.ndarray or xr.DataArray,
    latitudes: np.ndarray or xr.DataArray
) -> np.ndarray or xr.DataArray:
    result = np.sum(
        (forecast_field - analysis_field) * np.cos(latitudes * np.pi / 180.0)
    ) / np.sum(np.cos(latitudes * np.pi / 180.0))
    return result

MAE

def absolute_bias(
    forecast_field: np.ndarray or xr.DataArray,
    analysis_field: np.ndarray or xr.DataArray,
    latitudes: np.ndarray or xr.DataArray
) -> np.ndarray or xr.DataArray:
    result = np.sum(
        np.abs(forecast_field - analysis_field) * np.cos(latitudes * np.pi / 180.0)
    ) / np.sum(np.cos(latitudes * np.pi / 180.0))
    return result

SD

def std(
    forecast_field: np.ndarray or xr.DataArray,
    analysis_field: np.ndarray or xr.DataArray,
    latitudes: np.ndarray or xr.DataArray
) -> np.ndarray or xr.DataArray:
    bias_result = bias(forecast_field, analysis_field, latitudes)
    result = np.sqrt(
        np.sum(
            np.power(forecast_field - analysis_field - bias_result, 2) * np.cos(latitudes * np.pi / 180.)
        ) / np.sum(np.cos(latitudes * np.pi / 180.))
    )
    return result

MSE

def mse(
    forecast_field: np.ndarray or xr.DataArray,
    analysis_field: np.ndarray or xr.DataArray,
    latitudes: np.ndarray or xr.DataArray
) -> np.ndarray or xr.DataArray:
    result = np.sum(
        np.power(forecast_field - analysis_field, 2) * np.cos(latitudes * np.pi / 180.0)
    ) / np.sum(np.cos(latitudes * np.pi / 180.0))
    return result

ACC

def acc(
    forecast_field: np.ndarray or xr.DataArray,
    analysis_field: np.ndarray or xr.DataArray,
    climate_field: np.ndarray or xr.DataArray,
    latitudes: np.ndarray or xr.DataArray
) -> np.ndarray or xr.DataArray:
    forecast_climate = bias(forecast_field, climate_field, latitudes)
    obs_climate = bias(analysis_field, climate_field, latitudes)
    acc1 = np.sum(
        (forecast_field - climate_field - forecast_climate) * (analysis_field - climate_field - obs_climate) * np.cos(latitudes * np.pi / 180.)
    )
    acc2 = np.sum(
        np.power(forecast_field - climate_field - forecast_climate, 2) * np.cos(latitudes * np.pi / 180.)
    )
    acc3 = np.sum(
        np.power(analysis_field - climate_field - obs_climate, 2) * np.cos(latitudes * np.pi / 180.)
    )
    result = acc1 / np.sqrt(acc2 * acc3)
    return result

参考

https://www.cawcr.gov.au/projects/verification/

https://www.jma.go.jp/jma/jma-eng/jma-center/nwp/outline2019-nwp/pdf/outline2019_Appendix_A.pdf