厄尔尼诺介绍与指数计算

发表于 2019-06-29 分类于闭门造车阅读次数： Disqus：

引言

之前我们简单介绍了NetCDF数据格式以及如何用Python进行简单的处理和分析。近年来，随着全球气候的变化，各类气象极端事件也大有愈演愈烈之势，例如2014-2016年就发生了记载以来最强的厄尔尼诺事件。就本文将在此基础上，介绍厄尔尼诺指数的计算和分析过程。

厄尔尼诺介绍

基本概念

第一次接触“厄尔尼诺”的概念，还是上小学的时候，看我爸的工程师英语教材（上面有很多用于阅读理解的科普文章，很有意思）。

厄尔尼诺現象（西班牙語：El Niño），指东太平洋海水每隔数年就会异常升温的现象，是厄尔尼诺-南方振荡现象（El Niño-Southern Oscillation，简称ENSO）中，东太平洋升温的阶段。厄尔尼诺是西班牙语，意指“小男孩”，指的是耶稣，因为南美太平洋的变暖时期通常都在圣诞节附近，所以厄尔尼诺现象也被叫做“圣婴”现象。它与中太平洋和东太平洋（约在国际换日线及西经120度）赤道位置产生的温暖海流有关，厄尔尼诺-南方振荡现象是指中太平洋和东太平洋赤道位置海面温度的高低温循环。

此外，厄尔尼诺-南方振荡现象中的低温阶段称为拉尼娜现象（也称为反圣婴现象），是指东太平洋的海面温度高于平均值，以及西太平洋的气压较低及东太平洋的气压较高所造成。

厄尔尼诺影响

厄尔尼诺/拉尼娜是指赤道中、东太平洋海表大范围持续异常偏暖/冷的现象，是气候系统年际气候变化中的最强信号。厄尔尼诺/拉尼娜事件的发生，不仅会直接造成热带太平洋及其附近地区的干旱、暴雨等灾害性极端天气气候事件，还会以遥相关的形式间接地影响到全球其它地区天气气候并引发气象灾害。

2014-2016年发生的超强厄尔尼诺事件，导致全球气候异常变化，包括我国在内的多个国家都经历了严重的洪涝灾害或旱灾，导致多国粮食严重减产，人民生命财产安全受损。1998年的极强厄尔尼诺事件会造成我国长江流域的严重洪涝灾害，给人民生命财产安全和我国经济发展带来巨大影响。我国历史上长江流域的大水年（1969年、1983年、1987年、1991年和1998年）都发生在厄尔尼诺的次年。

enso-warm-episode-djf

厄尔尼诺指数

指数简介

厄尔尼诺现象的影响如此巨大，显然我们迫切需要了解目前甚至预测未来厄尔尼诺现象的状态。对人来说，最理想的情况莫过于有一个直观的指标来表征这一切。从历史上来看，我们也正是这么多的：将描述厄尔尼诺现象所需要的复杂因素归结为一个简单的数字，然后跟踪这个数字的变化。

目前，国际上流行的厄尔尼诺指数简单介绍如下：

Oceanic Niño Index (ONI)，这是NOAA的官方指标，基于中东部热带太平洋（Nino 3.4区域）的三个月海表温度异常平均值（sea surface temperature anomalies, SSTA）计算；
Niño 1~4 指数，将赤道中东太平洋划分为几个区域包括Nino 1,2,3,4和3.4区，对每个区域计算SSTA；
Souther OScillation Index (SOI)，塔希提（Tahiti）与达尔文（Darwin）站月平均海平面气压差序列的标准化值
Multivariate ENSO Index (MEI)，基于赤道太平洋海6个气象变量的加权线性组合计算得到；
其他指数，参见 https://mrcc.illinois.edu/mw_climate/elNino/climatology.jsp

enso_indices

我们可能会有疑问：一方面我们希望将ENSO简化为一个简单的指数，另一方面却又提出了各种各样的指数将问题复杂化，为什么不用一个最好的指数呢？

首先，ENSO作为一种气象现象，他本身就是人为定义出来的，我们可以认为有定性描述的现象是客观的，但是定量的界定却是主观的，因此在没有客观标准的前提下，我们并不能说哪一个指数就真正表征ENSO。

其次，ENSO本身作为一个极其复杂的气象事件，目前看来任何指数都是对这一事件的简化，就像盲人摸象一样，一个指数可能是ENSO在某一方面的投射，所以如果像尽可能全地描述ENSO就需要多个指数。

而且，一个非常重要的原因是，不用指数反映了ENSO的不同侧面，而不同地区受ENSO的影响也是不同的。说白了，ENSO本来也是为了解释全球范围内的气象变化，那么不同国家的气象机构关心的是其受EMSO影响最大的方面，因此在指数的选取上不同国家有着不同的偏好。例如，热带太平洋岛屿上的人们可能更关心SST，而在远离热带太平洋的国家，他们可能更关心海平面气压的大规模变化。

当然，在这些指数里面，最常用的是ONI和NINO3.4，两者都是基于nino 3.4区的SSTA计算，区别在于计算SSTA均值的滑动窗口和判别标准，例如在判定厄尔尼诺时，ONI通常采用三个月滑动平均SSTA连续5个月超过0.5摄氏度，而NINO3.4为5个月滑动平均连续6个月超过0.4摄氏度。

我国2017年新制定的国标实际上选用的就是ONI（虽然写的是NINO3.4）。

国标规定

2014年至2016年的超强厄尔尼诺事件暴露出厄尔尼诺/拉尼娜事件认定存在分歧的问题。针对国内尚缺乏统一的厄尔尼诺/拉尼娜事件判别标准的现状，为了规范判别标准，2017年，国家气候中心牵头制定《厄尔尼诺/拉尼娜事件判别方法》国家标准。

标准对主要海温监测关键区及指数、判别方法等做了明确规定：

标准采用了国际通用的NINO3.4来判别事件的发生、强度、持续时间等，用东部型指数（ $I_{EP}$ ）和中部型指数（ $I_{CP}$ ）来判定事件的类型。
NINO3.4的3个月滑动平均绝对值达到或超过0.5℃、持续至少5个月，判定为一次厄尔尼诺/拉尼娜事件（指数≥0.5℃为厄尔尼诺事件；指数≤-0.5℃为拉尼娜事件）；以NINO3.4满足事件判定的时间为持续时间。
在事件过程中，NINO3.4的3个月滑动平均绝对值达到最大的时间和数值分别定义为事件的峰值时间和峰值强度。然后，用 $I_{EP}$ 和 $I_{CP}$ 来判定事件的类型：事件过程中 $I_{EP}$ 的绝对值达到或超过0.5℃且持续至少3个月的类型判定为东部型事件；事件过程中 $I_{CP}$ 的绝对值达到或超过0.5℃且持续至少3个月的类型判定为中部型事件。
若一次事件中同时包含上述两种情况、存在两种类型间的转换，则将事件峰值所在类型定义为事件主体类型，另一种为非主体类型，整个事件的类型以事件主体类型为准。

用 $I_{EP}$ 和 $I_{CP}$ 可以衡量出两类事件对中国气候不同的影响关系。就中国东部地区夏季降水而言， $I_{CP}$ 所反映的中部型影响与 $I_{EP}$ 所反映的东部型影响都很明显，且基本上影响相反。例如，中部型会导致华北-西北东部地区多雨、长江流域少雨，东部型则使得华北-西北地区东部少雨、长江流域多雨。相比之下，NINO3.4区介于东部型和中部型之间，其海温异常的影响也是介于两种类型指数影响之间，即无论是正负影响区的振幅都明显弱于 $I_{EP}$ 和 $I_{CP}$ 的影响。这充分说明，从气候影响来看，采用IEP和ICP作为监测指标来判别事件，要比单纯用ININO3.4对我国气候异常更具有指示意义。

（以上来源：《中国气象报》2017年6月20日三版）

计算公式

厄尔尼诺/拉尼娜事件的主要监测关键区，包括NINO1+2区（90°W-80°W，10°S-0°）、NINO3区（150°W-90°W，5°S-5°N）、NINO4区（160°E-150°W，5°S-5°N）和NINO3.4区（170°W-120°W，5°S-5°N）。根据我国国标，需要的指数分别计算如下（事件判别标准请查阅国标）：

ENSOindices_locations

NINO3指数

NINO3区SSTA 的平均值。

NINO4指数

NINO4区SSTA 的平均值。

NINO3.4指数

NINO3.4区SSTA 的平均值。

东部型厄尔尼诺/拉尼娜指数

I_{EP} = I_{NINO3}-\alpha\times{I_{NINO4}}

其中，

$I_{EP}$ —东部型厄尔尼诺/拉尼娜指数，单位为摄氏度（℃）；
$I_{NINO3}$ —NINO3指数，单位为摄氏度（℃）；
$I_{NINO4}$ —NINO4指数，单位为摄氏度（℃）；
$\alpha$ —当 $I_{NINO3}\times{I_{NINO4}>0.4}$ 时， $\alpha=0.4$ ；否则 $\alpha=0$ 。

中部型厄尔尼诺/拉尼娜指数

I_{CP} = I_{NINO4}-\alpha\times{I_{NINO3}}

其中，

$I_{CP}$ —中部型厄尔尼诺/拉尼娜指数，单位为摄氏度（℃）。

ONI时间序列计算

导入数据

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import datetime
import xarray

mpl.rcParams['figure.figsize']=(20,9)
mpl.style.use('ggplot')

我们注意到计算ONI只需要海洋表面温度的历史数据，可以从NOAA Earth System Research Laboratory Physical Sciences Division (https://www.esrl.noaa.gov/psd/) 下载该数据。

1
2
3

sst_file = 'data/skt.mon.mean.nc'
sst_data = xarray.open_dataset(sst_file)
sst_data

<xarray.Dataset>
Dimensions:  (lat: 94, lon: 192, time: 857)
Coordinates:
  * lat      (lat) float32 88.542 86.6531 84.7532 82.8508 80.9473 79.0435 ...
  * lon      (lon) float32 0.0 1.875 3.75 5.625 7.5 9.375 11.25 13.125 15.0 ...
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 1948-03-01 ...
Data variables:
    skt      (time, lat, lon) float32 ...
Attributes:
    description:  Data is from NMC initialized reanalysis\n(4x/day).  It cons...
    platform:     Model
    Conventions:  COARDS
    NCO:          20121013
    history:      Tue Jul  6 00:05:45 1999: ncrcat skt.mon.mean.nc /Datasets/...
    title:        monthly mean skt.sfc from the NCEP Reanalysis
    References:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...

lon = sst_data.lon.values
lat = sst_data.lat.values
nctime = sst_data.time.values
nctime = nctime.astype('M8[ms]').astype('O')
t_unit = sst_data.skt.units
skt = sst_data.skt.values

获取Nino3.4区的经纬度范围和需要计算的时间范围：

# subregion for nino3 area
idx_lat = (lat>=-5.0) & (lat<=5.0)
idx_lon = (lon>=210.0) & (lon<=270.0)

years = np.array([idx.year for idx in nctime])
idx_time = (years >= 1981) & (years <= 2010)

提取给定范围内的数据：

lat_oni = lat[idx_lat]
lon_oni = lon[idx_lon]
dates_oni = nctime[idx_time]

#采用Mask数组索引多维数组与普通方法会有些不同
skt_oni = skt[idx_time][:,idx_lat][:,:,idx_lon]

指数计算

看我们前面的做法，基本思路是将数据放到ndarray里面计算，但是！这些变量和我们要计算的指数本质上都是时间序列啊，如果能用Pandas是不是会简单很多呢？

1
2
3

sst_df = sst_data.to_dataframe()
print(sst_df.index.names)
sst_df.head()

['lat', 'lon', 'time']

			skt
lat	lon	time
88.542	0.0	1948-01-01	-39.854198
		1948-02-01	-37.212414
		1948-03-01	-32.056770
		1948-04-01	-26.997002
		1948-05-01	-6.077735

为了看上去更加直观，我们只保留time作为索引，而把lat和lon恢复为普通的列，其实如果把普通列进行排序后进行查询，可能比用索引速度会更快一些。

sst_df.reset_index(inplace=True)
sst_df.sort_values(['time', 'lat','lon'], inplace=True)
sst_df.set_index('time',inplace=True)

sst_selected = sst_df[(sst_df.lat<=5)&(sst_df.lat>=-5) & (sst_df.lon>=210.0) & (sst_df.lon<=270.0)]
sst_selected.head()

	lat	lon	skt
time
1948-01-01	-4.76184	210.000	27.212244
1948-01-01	-4.76184	211.875	27.117422
1948-01-01	-4.76184	213.750	26.931597
1948-01-01	-4.76184	215.625	26.635796
1948-01-01	-4.76184	217.500	26.255482

接下来以月为单位，计算每个月的历史均值（需要注意的是，由于存在季节性，不应该计算年度的均值）。

此外，根据WMO的一些规定，在计算标准气候平均值（standard climatological normals）时，根据WMO的规定，以历史上最近整30年的30个样本平均值作为相应气候要素常年值的参考点，并10年更新一次。也就是说，如果我们计算近年来的SSTA，应当使用1981~2010年的其后资料统计值。

让我哭笑不得的是，虽然DataFrame看上去很美好，但是在后面的计算中速度却不太让人满意，提高速度的方法是将其向量化（Vectorization），其中一个方法就是使用numpy——所以，我。。。

# 按行填充，依次为 年 月 维度 经度
skt_oni = np.reshape(skt_oni, (30, 12, 6, 33))
sst_mean = np.mean(skt_oni, axis=0)
sst_mean.shape

(12, 6, 33)

计算海水温度距平值（SSTA）

num_repeats = 30 # 30 years
sst_mean = np.vstack([sst_mean] * num_repeats)
sst_mean = np.reshape(sst_mean, (num_repeats,12,6,33))
# ssta over the area
ssta = skt_oni - sst_mean

# calcuate overall average
oni_ts = np.mean(ssta, axis=(2,3))
oni_ts = np.hstack(oni_ts)

结果可视化

dates = pd.date_range('1981','2011', freq='MS', closed='left')
oni_df = pd.DataFrame({'date':dates,'oni':oni_ts})
oni_df.set_index('date', inplace=True)
plt.figure(figsize=(20,9))
plt.plot(oni_df.index, oni_df.oni, 'black', alpha=1, linewidth=2)
plt.fill_between(oni_df.index, 0, oni_df.oni, oni_df.oni>0,
                color='red', alpha=0.75)
plt.fill_between(oni_df.index, 0, oni_df.oni, oni_df.oni<0,
                color='blue', alpha=0.75)
plt.xlabel('Years')
plt.ylabel('[$^oC$]')
plt.title('ONI SSTA 30-year (1981-2010)', fontsize=12)

Text(0.5, 1.0, 'ONI SSTA 30-year (1981-2010)')

output_34_1

1 2	print(oni_df[oni_df.oni>3].index) print(oni_df[oni_df.oni<-1.5].index)

DatetimeIndex(['1982-12-01', '1983-01-01', '1997-10-01', '1997-11-01',
               '1997-12-01', '1998-01-01'],
              dtype='datetime64[ns]', name='date', freq=None)
DatetimeIndex(['1988-05-01', '1988-06-01', '1988-07-01', '1988-08-01',
               '1988-10-01', '1988-11-01', '1988-12-01', '1999-12-01',
               '2000-01-01', '2007-10-01', '2007-11-01', '2007-12-01',
               '2008-01-01', '2010-10-01', '2010-11-01', '2010-12-01'],
              dtype='datetime64[ns]', name='date', freq=None)

基于上面的数据，我们大致可以判断82-83年和97-98年的厄尔尼诺较为严重，而88年和11年的拉尼娜事件较为严重。在两次强厄尔尼诺之后，我国发生了严重的洪灾；而在两次拉尼娜之后，我国则出现了相当程度的旱灾，厄尔尼诺/拉尼娜似乎确实与我国的自然灾害有某种关系。当然，如果希望得到更可靠的结论，一方面需要更长的时间尺度，另一方面最好能够有机理上的解释。

小结

我们简单介绍了厄尔尼诺的基本概念和指数，并且演示了基于netCDF4格式数据，一步一步计算厄尔尼诺系数的过程。本文现在做的还是较为简单的，其实可以做更加深入的分析，例如计算厄尔尼诺指数与其他气象变量，例如我国某些区域的降雨量等之间的关系。

参考文献

https://www.climate.gov/news-features/blogs/enso/why-are-there-so-many-enso-indexes-instead-just-one
“The Climate Data Guide: Nino SST Indices (Nino 1+2, 3, 3.4, 4; ONI and TNI).” Retrieved from https://climatedataguide.ucar.edu/climate-data/nino-sst-indices-nino-12-3-34-4-oni-and-tni.
https://www.esrl.noaa.gov/psd/enso/
https://www.weather.gov/fwd/
https://en.wikipedia.org/wiki/El_Niño
http://qh323.com/portal/getNewsDetails?name=encyclopedias&id=055a483584e04596ad447873e13fa29d
厄尔尼诺/拉尼娜事件判别方法 http://www.cma.gov.cn/root7/auto13139/201703/t20170306_397953.html
https://climatedataguide.ucar.edu/climate-data/sst-data-sets-overview-comparison-table
https://www.linkedin.com/pulse/create-nino3-time-series-chonghua-yin/
http://public.wmo.int/en/resources/meteoterm