写在前面

当前对于高温热浪事件指标的定义缺乏统一标准，不同学者有不同的算法。Russo（2015）提出了一种每日热浪幅度指数（HWMId，the heat wave magnitude index daily）计算方法，适合对全球的热浪强度进行计算和对比分析。该方法也被wiki百科对热浪指数的介绍中采用（wiki）。
该方法作者提供了相关的R语言包和函数，本文对此进行原理和代码的介绍。

放上wiki上的一张图（Zampieri, M., 2016）

1. 原理介绍

HWMId定义为每年最大热浪事件的幅度，包括热浪事件强度、事件开始时间和持续时间。该指数将高温热浪事件定义为连续3天日高温大于热浪事件温度阈值的一段时期。高温热浪事件阈值根据基准期31天窗口内90%分位数的日高温进行计算。
若将1981-2010年作为基准期（基准期最少为30年），若判断某日d是否发生高温热浪，则高温热浪事件阈值计算公式为：

U表示数据集的联合，是y年中第i天的日高温。该日的高温热浪幅度（Magnitude）为：

是高温热浪事件发生日d的最高气温，和分别是基准期31天窗口内日高温时间序列的第25%分位数和第75%分位数。高温热浪事件往往发生数天，则整个事件的幅度为每日幅度之和。

更详细的介绍请参见：Russo（2015）。

2. 实操

HWMId（RDocumentation）计算可使用R包extRemes（cran，RDocumentation）实现。

2.1 关键函数

该指标计算主要通过以下函数实现：
hwmid(yTref, Tref, yTemp, Temp)
输入数据分别为：参考数据的开始年份，参考数据，待计算数据的开始年份，待计算数据；
输出的结果分别为：相应年份的热浪指数，持续天数，开始时间；每天的参考温度；75th和25th分位数数值

2.2 输入数据：

yTref和yTemp均是具体的年份值，如1970、2000。Tref和Temp均是一维矩阵数据。

• yTref代表Tref的开始年份。函数要求基准的参考数据最少为30年，因此Tref至少是包含30年数据的矩阵，时间步长为1天，单个站点的一维数据。

• yTemp代表Temp的所在年份。Temp按照函数要求为1整年的数据，时间步长为1天，单个站点的一维数据。但具体时间长度可以按需修改。Temp既可以位于基准期内，也可以位于基准期外。

2.3 输出数据：

仅介绍“相应年份的热浪指数，持续天数，开始时间”三个常用结果，其他可输出的值可在参考文档中自行查看。
通过原理介绍中列出的公式大家也可看出该指标计算的热浪指数强度表示的是气温相对于一定阈值的相对大小，因此没有单位。持续天数和开始时间表示的均是年内最强一次热浪的持续天数和开始时间。

2.4 计算

hwmidFr <- hwmid(1980, Tref, year, Temp)

hwmid1 <- hwmidFr$hwmid[1] #相应年份的热浪指数
hwmid2 <- hwmidFr$hwmid[2] #持续天数
hwmid3 <- hwmidFr$hwmid[3] #开始时间

以上代码可分别进行多年多个站点或多像元的计算，以满足自己的计算需求。
对于基于栅格数据进行整个区域的运算，由于R语言不方便进行像元尺度的计算，之后我将介绍如何在matlab上直接实现从.nc数据直接算出每年的区域上的热浪强度，实现文章开头图片的效果。

References

Russo, S., Sillmann, J., Fischer, E.M., 2015. Top ten European heatwaves since 1950 and their occurrence in the coming decades. Environ. Res. Lett. 10. https://doi.org/10.1088/1748-9326/10/12/124003
Zampieri, M., Russo, S., di Sabatino, S., Michetti, M., Scoccimarro, E., Gualdi, S., 2016. Global assessment of heat wave magnitudes from 1901 to 2010 and implications for the river discharge of the Alps. Sci. Total Environ. 571, 1330–1339. https://doi.org/10.1016/j.scitotenv.2016.07.008

​