格陵兰近千年气候变化的多尺度分析

格陵兰近千年气候变化的多尺度分析1

杜建丽，林振山，俞鸣同，张真真，刘会玉

南京师范大学地理科学学院，南京 (210046)

E-mail: linzhenshan@

摘要：用经验模态分解（EMD）方法对格陵兰冰盖GISP2冰芯古气候代用指标δ18O序列进行分析，结果表明：北半球高纬地区气候波动具有3、6.5、12、24、49、96、213、468年准周期性波动，既有NESO的影响，也受制于太阳活动的周期。近1000年来气候变化总趋势出现历时约490年的中世纪暖期和历时约570年的小冰期，其间还存在次级的冷暖期变化；EMD第7，第8本征模函数IMF7和IMF8波动振幅以及EMD趋势项在1350年出现明显的转型，表明1350年为中世纪暖期和小冰期的分界。研究结果对北半球高纬地区近千年气候变化的波动规律和变化分期作了较详细的划分。

关键词：格陵兰冰芯；EMD；波动周期；小冰期

0. 前言

过去2ka的气候变化是“过去全球变化”（PAGES）和“气候变率与可预测性”（CLIVAR）两大国际研究计划的重要研究内容[1,2,3]，其中1 ka来的气候变化是现代器测资料与过去代用指标序列衔接的关键时段,也是人类活动影响不断加剧的时期,气候冷暖变化阶段和周期以及现代气温的快速增暖成为各国科学家关注的热题。

目前，对2ka来气候研究主要关注多种代用指标所揭示的中世纪暖期、小冰期在不同地区的反应，中世纪暖期、小冰期是否具有全球意义，暖期和小冰期内的次级气候波动及其峰谷值点，中世纪暖期向小冰期转型的分界，以及现代增暖期的变幅等热点问题。中世纪暖期以及随后的小冰期在欧洲、北美中高纬地区的多种代用指标里有所反映。Briffa K R, Lipp J, 等研究西欧多处树轮发现，大约950～1200AD气候暖期对树木生长和稳定同位素变化具有明显的影响[4]；Dansgaard W，Grove J M 和Stuiver M,等研究冰岛和格陵兰冰芯δ18O序列时指出其中具有中世纪暖期的信号[5，6，7]；Bodri 和Cermak从捷克深钻孔反演的温度变化也发现了中世纪暖期的存在[8]。中国东部和西部对中世纪暖期的反应差别较大。王绍武等根据史料、树木年轮、冰芯等代用资料建立了1000～2000 AD中国东部、西部及全国平均气温距平序列, 认为中国东部在1040～1130AD；1170～1260AD存在中世纪暖期,而西部无明显的中世纪暖期[9]。北半球小冰期研究已有较多的报道。捷克钻孔对古地温研究表明小冰期出现在1400～1500AD，1600～1700AD，最冷时期为1650AD；Stuiver M,对格陵兰冰芯δ18O 研究表明小冰期出现在1300～1900AD；王绍武研究中国近千年气候变化指出，中国东部小冰期出现1600～1690AD，1780～1870AD，西部的反应可提前至11世纪。

这些研究表明千年来气候变化出现中世纪暖期和小冰期现象，但是在各地区的反应有所差异，表现在暖期和冰期的起止年代不同，气候期内的次级波动次数和峰谷极值出现的年代也不同；中世纪暖期向小冰期气候转型的分界年代也不一致，以及千年来对气候变化影响显著的小尺度高频波动的周期规律尚不清楚。这些研究多数从时间序列数据的直观变化中总结规律，发现事件的突变，或作简单的数据平滑处理，缺少严密的科学方法推导非平稳序列的数学模式，蕴涵在非平稳序列中的多尺度变化规律难以发现，阈值和临界点的确定带有较大的主观臆断，因此众说纷纭，难以达成共识。

1本课题得到教育部博士点基金项目(20060139010)“1ka-10ka全新世东亚季风的多尺度诊断和预测”的资助。

经验模态分解（EMD ）方法是目前最好、最有效地进行非平稳序列多尺度分析的方法[10]。1998年Huang [11]提出了经验模态分解，1999年又做了一些改进[12]。EMD 可以对一个信号序列进行不同尺度（频率）的波动或趋势逐级分解，产生一系列具有不同特征尺度的数据序列，所表达的物理意义明确。EMD 方法具有自适应的特性,适宜于非平稳信号的分解，而被美国NASA 宇航中心采用[10]。

1. 研究方法

1.1 经验模态分解(EMD)简介

自然界的波谱信号往往是多种不同尺度信号的迭代结果。被检测的信号综合了不同尺度信号的因素，使得信号波谱杂乱无章，难以分辨信号中隐含的不同尺度信号的物理意义。过去对第四纪样品代用指标检测数据的分析，多数以目视方法、概率统计方法做简单解释；或用简单的数据平滑处理后，判读数据的地学意义。这些处理没有把数据的原本真实的物理意义反映出来，或扭曲了数据的真实物理意义，使解释意义失真。

经验模态分解（EMD ）方法是一种处理非平稳数据序列的全新方法，其本质是对一个信号进行平稳化处理，将信号中不同尺度的波动或趋势逐级分解开来，形成一系列具有不同尺度的数据序列，每一个序列称为一个本征模函数（intrinsic mode function, IMF ）分量，其中不同尺度数据的物理意义得以保留，并很好地被提取和表达。最低频率的分量（IMF ，res ）代表原始信号的总趋势或均值的时间序列。

EMD 方法的处理过程相对简单，适合于对大容量数据的处理。具体方法[10]是：找出原始数据序列X(t)所有的极大值点，其中t 是时刻，t=1,…,T 。用三次样条函数拟合，形成原数据序列的上包络线；相应地找出X(t)所有的极小值点，同样用三次样条函数拟合，形成原数据序列的下包络线；上下包络线的均值为原数据序列的平均包络线m 1(t)；将原始数据序列X(t)减去平均包络，得到一个去掉低频的新数据序列h 1(t)。

h 1(t)=X(t)－m 1(t) （1）

通常情况下，一次处理后h 1(t)仍然不是一个平稳数据序列，因此，将其重复上述处理过程k 次。直到所得到的平均包络趋于零为止。

h k (t) = h (k-1)(t) －m k (t) (2)

由此，得到第一本征模函数（IMF1）分量C 1（t ）：

C 1(t)= h k (t) (3)

IMF1分量代表原始数据序列中最高频的组分，保留了原始数据中最高频信号的物理特征。原始数据中其他尺度的信号可以用同样的方法继续提取。将原始数据序列X(t)减去第一个IMF 分量C 1(t)，可以得到一个去掉高频组分的差值数据序列r 1(t)。对r 1(t)进行上述平稳

化处理过程可以得到第2个IMF 分量C 2(t)，

如此重复下去直到最后一个差值序列r n (t)不可再分解为止。此时，r n (t)代表原始数据序列的趋势或均值的时间序列：

r 2(t) =r 1(t)- C 2(t),···, r n (t) = r n-1(t)- C n (t) (4)

Huang 将这样的处理过程形象地比喻为“筛”过程。最后，原始的数据序列即可由这些IMF 分量以及一个均值或趋势项表示：

()()()1

.n

j n i X t C t r t ==+∑ (5)