NWPC消息平台：标准时间统计算法

本文属于介绍 NWPC 消息平台 系列文章。

简介

数值天气预报是基于数学物理方法客观定量计算未来天气演变的科学，是一个跨学科的复杂而严格的系统性工程 [1]。 数值天气预报业务系统通常由多个步骤构成，一般分为资料前处理、模式计算和产品后处理等三个步骤 [2]。 面向用户的最终产品由产品后处理步骤生成，产品生成的时效受前面步骤的影响。

NWPC 数值预报业务系统的前处理步骤需要从多种来源获取观测数据，包括 HPC 共享存储、CIMISS、FTP 等，并对这些资料进行处理，转为模式系统可以识别的输入数据。 受 HPC 文件系统 I/O 速度、网络传输速率、数据接口响应时间等因素影响，每个时次的资料前处理步骤运行时长会有一定的浮动。 模式同化和模式积分等计算任务需要使用大量计算节点，不同时次模式计算步骤的运行时长有一定的浮动。 产品后处理涉及大量文件系统 I/O 操作，运行时间也会有一定的浮动。 综上所述，数值预报业务系统各个时次的运行时间段常常不够稳定，模式产品的生成时间通常会在一定的区间范围内浮动。 所以，需要一种有效的方法计算模式产品的标准生成时间段。

本文介绍 NWPC 消息平台目前使用的一种计算标准时间段的方法。

资料

使用 NWPC 消息平台产品事件消息类型中的原始分辨率 GRIB 2 产品完成上传二级存储的消息作为产品生成的时间数据。 统计 GRAPES GFS、GRAPES MESO 10KM、GRAPES MESO 3KM 和 GRAPES TYM 四个模式所有时次的产品生成标准时间段。

GRAPES GFS 模式 10 月 3 日 06 时次 120 小时 GRIB 2 产品生成的事件消息数据如下所示。

{
  "app": "nwpc-message-client",
  "type": "production",
  "time": "2020-10-03T10:47:20.395949667Z",
  "data": {
    "event": "storage",
    "forecast_time": "120h",
    "name": "orig",
    "start_time": "2020-10-03T06:00:00Z",
    "status": 1,
    "stream": "oper",
    "system": "grapes_gfs_gmf",
    "type": "grib2"
  }
}

产品生成消息的详细介绍请参看文章《NWPC消息平台：产品事件消息》。

方法

因为产品生成时间的波动性，使用单一时间（例如平均值）不能很好地代表产品生成的普遍情况。 本文将产品标准时间表示为一个时间段，使用置信区间的上下界作为时间段的两个端点。

对于某个特定时效的产品，一天只有一个数据，计算置信区间的数据量太少，容易受到离群值的影响。 本文使用自助法计算置信区间，对数据进行多次重采样并求均值，得到均值数据，再计算均值的分位数，得到置信区间的上下界。

结果与分析

本文使用上述方法统计 GRAPES GFS，GRAPES MESO 10KM，GRAPES MESO 3KM 和 GRAPES TYM 四个模式所有时次的产品生成情况。 使用 2020 年 9 月 1 日至 11 月 30 日共 61 天的数据，重采样 10000 次，每次选择 20 个样本，计算 95% 置信区间。

标准时间计算结果

GRAPES GFS

GRAPES GFS 全球预报系统每天运行 4 次，其中 00 和 12 时次生成 240 小时预报产品，06 和 18 时次生成 120 小时预报产品。 下图是 2020 年 9 月 1-10 日 00 时次所有产品生成时间段，可以看到产品生成的起止时间有一定的波动性。

使用自助法计算置信区间的结果如下图所示。

四个时次标准时间的区间都在 20-30 分钟之间。

GRAPES MESO 10KM

GRAPES MESO 10KM 区域预报系统每天运行 2 次，每次生成 84 小时预报产品。

使用自助法计算置信区间的结果如下图所示。

GRAPES MESO 10KM 产品标准时间的区间在 10-20 分钟之间，12 时次比 00 时次更稳定。

GRAPES MESO 3KM

GRAPES MESO 3KM 区域高分辨率预报系统每天运行 8 次，每次生成 36 小时预报产品。

使用自助法计算置信区间的结果如下图所示。

标准时间的区间在 20-60 分钟之间。 03 和 15 时次中间时效有明显的突起现象，具体原因有待进一步检查。

GRAPES TYM

GRAPES TYM 区域台风预报系统每天运行 4 次，每次生成 120 小时预报产品。

使用自助法计算置信区间的结果如下图所示。

标准时间的区间在 15-30 分钟之间。 00 和 18 时次中间时效有明显的突起现象，具体原因有待进一步检查。

监控报警时间阈值

NMC 统一监控平台已对 GRAPES 系统模式产品进行监控，当产品生成延迟时，会发送报警短信。 根据以上标准时间统计结果，本文为 GRAPES GFS 和 GRAPES MESO 3KM 系统设置产品延迟的报警时间阈值，如下表所示。

结论

本文基于 NWPC 消息平台中的 GRIB 2 产品完成上传二级存储的消息数据，设计并实现一种计算产品生成标准时间段的方法。 对 GRAPES GFS，GRAPES MESO 10KM，GRAPES MESO 3KM 和 GRAPES TYM 四种确定性模式业务系统在 2020 年 9 月 1 日至 11 月 30 日的消息数据进行分析，制作所有时效产品的标准时间段折线图， 并根据标准时间段的计算结果为 NMC 统一监控平台设置 GRAPES GFS 和 GRAPES MESO 3KM 两个模式产品延迟报警时间阈值。

本文结论如下：

1. 使用自助法计算置信区间作为产品生成的标准时间，能有效反映模式系统产品的普遍生成情况，从一定程度上也能体现产品生成时间的不稳定性。
2. 当前的计算结果尚有不合理的地方，计算方法有改进的空间，尤其需要评估离群点对整个计算结果的影响。

实现

使用产品时间消息 df_start_hour 作为原始输入数据，为每个时效 forecast_hour 单独计算标准时间段 standard_time。

for hour_index in forecast_hours_range:
    forecast_hour = forecast_hours[hour_index]
    clock_df_hour = df_start_hour[df_start_hour["forecast_hour"] == forecast_hour]

    means = self._get_means(
        clock_df_hour,
        bootstrap_sample=self.bootstrap_sample,
        bootstrap_count=self.bootstrap_count,
    )

    standard_time = self._get_bound(
        means,
        forecast_hour=forecast_hour,
        quantile=self.quantile,
    )

    standard_times.append(standard_time)

_get_means() 函数返回多次重采样数据的均值列表，其中 bootstrap_sample 是每次采样的个数，bootstrap_count 是重采样的次数。

def _get_means(
        self,
        clock_df,
        bootstrap_sample,
        bootstrap_count
) -> typing.List[pd.Timedelta]:
    means = []
    for i in range(bootstrap_count):
        mean = self._get_mean(clock_df, bootstrap_sample)
        means.append(mean)
    return means

_get_mean() 函数完成单次重采样，并返回均值

def _get_mean(
        self,
        clock_df: pd.DataFrame,
        bootstrap_sample: int,
) -> pd.Timedelta:
    sampled_data = clock_df["clock"].sample(
        n=bootstrap_sample,
        replace=True,
    )
    return sampled_data.mean()

_get_bound() 函数返回置信区间的上下界。 其中 quantile 是置信区间的宽度，一般为 0.95 或 0.99。

def _get_bound(
        self,
        means: typing.List[pd.Timedelta],
        forecast_hour: int,
        quantile: float,
) -> typing.Dict:
    bdf = pd.DataFrame(means).applymap(lambda x: x.ceil("s"))
    upper_bound = bdf.quantile(
        quantile + (1 - quantile) / 2,
        numeric_only=False,
        interpolation="nearest"
    )
    lower_bound = bdf.quantile(
        (1 - quantile) / 2,
        numeric_only=False,
        interpolation="nearest",
    )
    return {
        "forecast_hour": forecast_hour,
        "upper_bound": upper_bound[0],
        "lower_bound": lower_bound[0],
    }

参考

参考文献

[1] 沈学顺，王建捷，李泽椿，陈德辉，龚建东. 2020. 中国数值天气预报的自主创新发展. 气象学报，78(3)：451-476

[2] 佟华,胡江林,张玉涛.GRAPES模式后处理技术改进应用研究[J].气象科技,2020,48(4):511~517

NWPC 消息平台项目

nwpc-oper/nwpc-message-client

nwpc-oper/nwpc-message-tool

产品消息

《NWPC消息平台：产品事件消息》