格陵兰近千年气候变化的多尺度分析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.nj n i X t C t r t ==+∑ (5)1.2 代用指标的含义和来源研究古气候变化规律和机制最主要的参数是温度变化。
在古气候变化多种代用指标中，冰盖降水δ18O变化对气温变化的敏感性和时标的测定准确性相对较高。
其基本原理是水汽海陆循环中，δ18O是地表温度的函数。
温度升高时，降水中δ18O较高；温度降低时，降水中δ18O 较低。
因此，用δ18O变化的时间序列直接表达气温的变化规律是切实可行的（请引用前人文献）。
格陵兰冰盖钻孔GISP2冰芯δ18O变化代表了北半球高纬地区近千年气候变化的规律。
对这一规律的深入研究和拓展，有助于对北半球未来气候变化的预测。
90年代，格陵兰冰盖研究项目取得了大量有关第四纪气候变化代用指标的数据。
美国华盛顿大学第四纪同位素实验室完成的格陵兰冰盖GISP2冰芯δ18O数据序列，时间跨度从公元818～1987年,具有每年间隔，信号分辨率高，数据连续性好，间隔均匀，数据量大，内含气候变化信息丰富等优点，适合于做多尺度变化研究（注明数据来源：如世界数据中心WDC 的古环境数据库，或原文出处）。
2. 格陵兰冰盖GISP2冰芯δ18O序列的EMD分析从数据的直角坐标显示（图1），几乎找不到什么规律,既看不出δ18O的多时间尺度特征性,也看不出周期变化在时间域中的分布情况，因此，必须对数据做必要的处理。
图1 格陵兰冰芯δ18O近1000年来的趋势利用EMD运算软件，对格陵兰冰盖GISP2冰芯公元818a～1987a近千个数据的δ18O 序列进行分解，得到的8个本征模函数（intrinsic mode function, IMF）分量及其数据趋势分量res（图2，表1）。
每个IMF分量图形表示不同尺度或一个窄波段的δ18O变化特征信号。
I M F 1I M F 2I M F 3I M F 4I M F 5I M F 6I M F 7I M F 8R图2 近千年格陵兰冰芯的δ18O 变化的IMF 分量及其趋势量res表1 不同IMF 分量的方差贡献率冰芯IMF1 IMF 2 IMF 3IMF 4IMF 5IMF 6IMF 7 IMF8 Res 方差贡献率（%） 44.17 20.11 14.5210.61 6.39 1.08 1.06 1.1 0.97周期 3 6.5 12 24 49 96 213 468 图2中每一个IMF 分量信号特征的还原性很强，变化规律符合自然信号非线性变化特点。
其波动周期虽然不是严格的函数周期，但是具有相对稳定的变化准周期，可以将其求取平均周期，并将每种尺度信号波动频率和振幅对原数据总体特征影响程度用方差贡献率表示出来（表1）。
3. 结论3.1 多尺度变化周期IMF 的各个分量具有三个不同数量级相对差异的特征，表现在频率和振幅的方差贡献率＞20%，20%～10%，＜10%三个数量级。
对应的分量分别是IMF 1，IMF 2；IMF 3，IMF 4；IMF 5，IMF 6，IMF 7，IMF 8。
最大数量级分量中，IMF 1表示的是一个3年的准周期性波动，IMF 2分量表示6～7年准周期波动。
第二数量级分量中，IMF 3分量表示12年准周期的波动，IMF 4表示24年准周期波动。
较次要数量级的IMF 5、IMF 6、IMF 7、IMF 8分量表示的分别是准49、96、213、468年准周期的波动。
趋势项（res ）表示近1000年来的气温变化的总体趋势。
从IMF 分量的方差贡献率来看，IMF 1，IMF 2的波动频率和振幅的方差贡献率之和达到64.28%，是影响千年气候变化的主导周期，分别是3年和6～7年的准周期。
IMF 3，IMF 4的方差贡献率之和为25.13%，在千年气候变化中次于3年和6～7年的周期的影响。