1 引言

近几十年来，全球气候变暖越来越受到学者们的关注^{[1, 2, 3, 4, 5]}。对过去温度的重建是研究不同时间尺度气候变化特征的重要手段^{[6, 7, 8]}。由于器测数据大多不超过150年^[9]，所以目前主要依靠气候代用资料(树轮、 冰芯、 历史文献、 石笋、 孢粉等)来重建历史时期的温度变化^{[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]}。基本方法是: 首先得到校准时期(即器测数据和代用数据重叠时段)代用数据和器测数据间的线性函数关系，然后利用所得关系和代用数据重建过去的温度变化^{[22, 23, 24, 25]}。

区域优化平均法(Optimal Regional Averaging Method)是一种计算区域平均温度序列的方法，最早由Soviet提出，并由Gandin^[26]和Kagan^[27]进一步发展而来。与简单的面积加权相比，该方法可在使用较少数据的情况下得到较精确的重建结果^{[28, 29]}。Shen等^{[30, 31]}利用中国北方23个观测站月平均地表温度研究1961~1990年中国东北地区的平均温度，并得到较好的结果。但目前区域化平均法计算协方差模式时，所使用的经验公式只局限于局部区域，在较大空间范围使用时往往存在问题。王聪等^[29]提出一种适用于较大区域的协方差模型(Covariance Pattern)，但该模型要求代用数据在时间上必须是连续的，这样会浪费较多数据，从而降低计算精度。因此，找到一种适用于较大区域，且可以使用非连续性数据的协方差模式是必要的。本文拟通过引入球谐分析法(Spherical Harmonic Analysis)改进现有方法来解决这一问题。

在面积为S的区域Ω中，T(r，t)表示区域内点r，时刻t的温度。在理想状况下，区域内所有点的温度均已知，t时刻区域的平均温度为:

但在实际计算中，通常只能得到区域内部分点的温度，因此不能直接用公式(1)来计算平均温度，但可以用对已知的各点温度进行加权平均法来计算平均温度。假设区域Ω内已知N个点的温度，w(i)为i点的权值，则区域Ω的近似平均温度可表示为:

权值满足的条件为:

因此，在t时刻真实平均温度与近似平均温度的均方误差可表示为:

其中，ρ(i，j)=T(i，t)*T( j，t)。

当ε²最小时，与

的差值最小，所计算的平均温度最为精确，此时的权值即为最优权值。选取拉格朗日乘子法计算ε²的极值，从而确定最优权值。选取-2λ为拉格朗日乘子，集合式(3)和(4)可构造新函数F:

将函数F对各自变量求导并令导数为零，可得如下方程组:

其中为协方差模式，表示点i位置的温度与整个研究区域的平均温度的相关性，其计算方法将在下节中介绍。

通过求解公式(6)即可得到最优权值w(i)，再将w(i)代入公式(2)即可求得区域内的平均温度。

通过协方差模式的计算式可以看出，只要计算出研究区域内的各站点的温度，即可直接计算出，而通过球谐分析法可以较为容易地计算出温度分布。

球谐分析是一种模拟球面上物理量随位置变化的方法，已广泛应用于地球科学和天文学研究中^{[32, 33]}。

球谐函数是球坐标系下拉普拉斯方程角度部分的一组解，其正则化形式可以表示为:

其中P^m_n为伴随勒让得多项式，θ为π/2减纬度，φ为经度。

可以证明正则化的球谐函数是完备正交的。因此，可以将全球各站点温度用球谐函数展开。

其中，a^m_n(t)和b^m_n(t)为球谐系数。

为了实现球谐函数的求解，假设全球温度场可以用有限项球谐函数得到，取公式(8)中(N＋1)²项近似代替:

利用Tikhonov和Arsenin^[34]的正则化方法(Regularization Method)求解公式(9)即可求得球谐系数，再将球谐系数代回公式(9)即可计算出全球温度场。藉此就可计算出协方差模式。

本文首先采用Climatic Research Unit(CRU)的数据集CRUTEM4(2014)计算全球1951~2013年的年平均温度，以检验本文改进的区域优化平均法的性能。

为了模拟站点在时间上不连续的情况，采取以下步骤计算温度序列(以30~40个站点组为例)，并在图 1中给出3组站点各年站点选取结果:1)选定年份； 2)在30~40之间随机选取一个数X； 3)在选定年份的全球站点中随机选取Ⅹ个站点； 4)使用球谐分析法利用这些站点数据计算全球各站点温度； 5)利用全球各站点温度计算协方差模式； 6)使用区域优化平均法计算该年平均温度； 7)重复以上步骤(1)~(6)从而得到整个时间段内的年平均温度序列。

图 1(Fig. 1)

图 1 3个不同站点组计算过程中各年站点选取结果 Fig. 1 Variations of station number over time with respect to three different groups of stations periods in each year

本文应用上述方法计算了30~40、 40~50和50~60个站点3种情况下的全球年平均温度序列并与CRU结果(使用面积加权平均法)进行对比(图 2)，并计算了3种站点段相对于CRU结果的累积误差(图 3)。

图 2(Fig. 2)

图 2 利用改进的区域优化平均法重建的全球年平均 温度变化序列及其与CRU的对比 Fig. 2 Comparison between global average temperature obtained by improved Optimal Regional Average Method with different stations and global average temperature of CRU

图 3(Fig. 3)

图 3 3种站点段结果的累积误差的比较 Fig. 3 Comparison of cumulative error between the results of three stations periods

从图 2和3中可以看出，应用改进的区域优化平均法计算得到的全球年温度序列与CRU的结果相比，在使用时间上不连续的站点的情况下仍可以很好地描述温度序列的变化，且计算精度随站点的增多而增加。图 3中的累积误差为相对于CRU结果的误差(重建结果减去CRU的结果)绝对值的累计值和。

为了进一步验证结果的可靠性，本文选用简单的算术平均法进行对比(图 4和5)。从图 5可以看出，使用相同站点的情况下，改进的区域优化平均法的累积误差小于算术平均法结果。

从以上结果可以看出，改进的区域优化平均法可以较好地利用非连续的站点数据计算温度序列，且改进的区域优化平均法的计算精度优于简单的算术平均法。由于器测资料太短，在重建过去温度序列时只能依赖于代用资料。代用资料有可能在时间上不连续，也可能在记录长度上不一致。本文的方法无疑对于这类数据的使用有重要的价值。

图 4(Fig. 4)

图 4 使用相同站点数据的情况下改进的区域优化平均法 计算全球年平均温度的结果与简单算术平均结果的对比 Fig. 4 Comparison of global average temperature obtained by improved Optimal Regional Average Method with that derived from simple arithmetic average using the same dataset

图 5(Fig. 5)

图 5 使用30~40站点的情况下改进的区域优化平均法 与简单算术平均法的累积误差比较 Fig. 5 Comparison of cumulative error between improved Optimal Regional Average Method and simple arithmetic average method using 30~40 stations

本文选取Schneider等^[35]文章中的15组树轮最大晚材密度(Maximum-latewood-density，简称MXD)序列(表 1)重建过去1400年(公元600~2002年)北半球夏季平均温度序列，并与Schneider等^[35]原有的重建结果进行对比。表 1中的r为1901~1976年MXD序列与邻近的CRUTEM4v网格点的相关系数。

表 1(Table 1)

表 1 15组MXD序列的站点信息Table 1 Site information of the 15 MXD-datasets 站点名称 国家 经纬度 研究树种 年表长度/A.D. 样本量 r 参考文献 *表示90％的置信水平 **表示95％的置信水平 Alaska US 68.8°N，142.4°W Picea glauca 1073~2002 246 0.31** Anchukaitis等，2013[36] Altai Russia 50.0°N，88.°E Larix sibirica 462~2007 59 0.60** Myglan等，2012[37] Athabasca Canada 52.3°N，117.3°W Picea engelmannii 1072~1994 102 0.56** Luckman 等，2005[38] Campbell Canada 68.3N，133.3°W Picea glauca 1175~1992 58 0.47** ITRDB[39] Jmtland Sweden 63.5°N，15.5°E Pinus silvestris，Picea abies 1107~1978 158 0.63** Schweingruber 等，1988[40] Lauenen Switzerland 46.4°N，7.3°E Picea abies 982~1976 206 0.33** Schweingruber等，1988[40] Ltschental Switzerland 47.5°N，7.5°E Larix decidua 735~2004 180 0.62** Büntgen 等，2006[41] Mangazeja Russia 66.7°N，82.3°E Larix sibirica，Picea obovata 1246~1990 143 0.70** ITRDB[39] Nscan Finland 67.5°N，22.5°E Pinus silvestris -181~2006 587 0.82** Esper等，2012[42] Polarural Russia 66.9°N，65.6°E Larix sibirica，Picea obovata 778~2006 141 0.76** Briffa等，2013[43] Pyrenees Spain 42.5°N，2.5°E Pinus uncinata 1044~2005 203 0.40** Büntgen 等，2008[44] Quebec Canada 57.5°N，76.0°W Picea mariana 1352~1989 83 0.59** ITRDB[39] Tornetrsk Sweden 68.2°N，19.7°E Pinus silvestris 441~2010 124 0.79** Melvin 等，2013[45] Tyrol Austria 47.5°N，12.5°E Picea abies 1047~2003 227 0.42** Esper 等，2007[46] Zhaschiviersk Russia 67.5°N，142.6°E Larix sibirica 1311~1991 55 0.30* ITRDB[39]

表 1 15组MXD序列的站点信息 Table 1 Site information of the 15 MXD-datasets

在1850~2002年时段，本文计算了改进的区域优化平均法和Schneider等^[35]的重建结果相对于CRU结果的累积误差值。从图 6可以看出，与Schneider等^[35]结果比较，改进的区域优化平均法误差稍大，可能是因为使用改进的区域优化平均法重建温度序列时需要先计算重建区域内的温度分布，而本研究选取的15组站点数据均集中分布于北半球中纬度地区，使用这些数据计算的北半球温度分布可能会产生较大的误差，从而造成最终的结果误差偏大。

图 6(Fig. 6)

图 6 在1850~2002年时段改进的区域优化平均法的重建结果和Schneider等(2015)[35]的重建结果的累积误差的比较 Fig. 6 Comparison of cumulative error of reconstruction results between improved Optimal Regional Average Method and Schneider et al.(2015)[35] during 1850~2002A.D.

本文将重建的过去1400年北半球夏季平均温度序列与Schneider等^[35]的结果进行对比(图 7)，并计算了公元800~1100年、 1300~1500年、 1600~1800年以及1900年以后几个特征时期的相对于20世纪器测夏季温度的距平(表 2)。从图 7和表 2中可以看出，相较于Schneider等^[35]的结果，本文重建的结果显示中世纪(公元800~1100年)的平均温度比14和15世纪的平均温度高，即中世纪暖期更为显著，而Schneider等^[35]的结果恰恰相反，没有显示更为明显的中世纪暖期。而17~19世纪的小冰期和20世纪气候增暖事件在本研究和Schneider等^[35]的重建中均有明显表现。两种方法重建结果在中世纪及14、 15世纪的反差可能是因为两种方法所选用的权值不同，Schneider等^[35]选择的是简单的相关系数加权法； 而本文的方法是通过球谐分析法计算校准时段逐年北半球的温度分布，并通过温度分布使用区域优化平均法逐年计算权值。相比而言，尤其当样本量充足(如大于30时，如图 4所示)本文的方法可能会更有效地重建出温度变化的真实特征。

图 7(Fig. 7)

图 7 改进的区域优化平均法的结果与Schneider等[35]的结果相对于1900~1999年CRU平均温度的距平 Fig. 7 Temperature anomalies relative to 1900~1999A.D. average of improved Optimal Regional Average Method and Schneider et al.[35]

表 2(Table 2)

表 2 几个特征时期相对于1900~1999年CRU平均温度的距平Table 2 Average temperature anomalies during several typical periods of temperature change relative to 1900~1999A.D. average 时期／A.D. 方法 Schneider等[35] 改进的区域优化平均法 800~1100 -0.149 -0.143 1300~1500 -0.024 -0.171 1600~1800 -0.272 -0.323 1900~2002 -0.010 -0.042

表 2 几个特征时期相对于1900~1999年CRU平均温度的距平 Table 2 Average temperature anomalies during several typical periods of temperature change relative to 1900~1999A.D. average

依据本文新建的温度序列，计算了中世纪最暖100年和最冷100年相对于20世纪CRU平均温度的距平。其中，最暖期出现在9世纪中叶到10世纪中叶(866~965年)，比20世纪平均值低0.064℃，最冷期出现在16世纪末期到17世纪(1578~1677年)，比20世纪平均值低0.342℃。

本文应用球谐分析法对原有的区域优化平均法进行改进，使区域优化平均法适用范围推广到北半球。利用CRU数据计算1951~2013年全球年温度序列，得到较好的结果: 在使用非连续性站点仍可以较好的描述温度序列的变化趋势； 重建结果的精度会随使用数据点数的增多而增大； 使用相同站点的情况下结果优于简单算术平均的结果。因此，利用改进的区域优化平均法可以较好地使用非连续性数据重建过去温度，扩大了代用数据的应用范围。本文利用Schneider等^[35]的数据对该方法进行了验证，重建了北半球过去1400年夏季平均温度，结果显示该方法是可行的，重建结果显示显著的中世纪暖期，小冰期及20世纪的暖期。

致谢 感谢审稿专家和编辑部老师建设性的修改意见，使文章得以完善。