1 引言

气温是天气及气候事件预报的重要要素之一，同时也是气候变化研究重要的指标之一，随着气候变化研究的不断深入，研究人员越来越迫切需要长序列、高质量的气候资料。由于现代气象观测数据的时间比较短，资料不太齐全，对更长期的区域性气候变化的趋势、时空特征及其成因的认识还远远不足。当前气候变化研究用的长序列气温序列资料源主要有3种：一是器测气温，它相对来说，是能够最为精确地代表观测点的气温资料，全球最为古老的气温仪器观测始于17世纪的中英格兰地区(Manley，1974)，但对大部分地区来说，其观测序列相对较短，尤其在中国西部地区；二是树轮密度年表(吴祥定等，1988；刘洪滨等，2003)、冰芯(姚檀栋等，2006)、沉积(何尧启等，2005；王君波等，2005)、历史文献(郑景云等，2003；Ge，et al，2005)等重建的气温代用资料序列，但由于代用资料毕竟不是直接观测的气温序列，其与器测资料的误差可能随着地域、时间及资料种类等有明显的不同。以树轮为例，由于存在季节差异，其取样间隔通常也不是常规意义上的月或年序列(袁玉江等，1999；张同文等，2008)；三是气候模式分析模拟得到的历史气温模拟序列，模式数据有着时空连续、不同气象要素之间协调性好的优势，但由于模式的初值误差、模式误差和大气系统的混沌性(Lorenz，1965)，使得数值模式存在不确定性。国际耦合模式比较计划第5阶段(Coupled Model Intercomparison Project Phase 5，CMIP5)(Taylor，et al，2012)由世界气候研究计划(WCRP)组织实施，全球近30个模式研发机构参与了该计划。CMIP5模式为气候变化机理研究和未来预估提供重要的数值模拟数据，与前几个阶段的CMIP模式相比，大部分CMIP5模式在物理过程、耦合碳循环等方面有明显的改进。目前，基于CMIP5模式历史模拟结果开展了很多分析研究(Xu，et al, 2012；姚遥等，2012；张艳武等，2016；Li，et al，2017)。因此，如何有效结合不同来源气温资料各自的优势特点，发展长时期的气温融合序列，是解决当前气候与气候变化业务和科研迫切需求的一种可能有效途径。

与经典统计学相比，贝叶斯统计学对于概率的理解更具有主观性，在总体信息和样本信息的基础上还利用到了先验信息，而先验信息既可以是客观的，也可以是主观的(黄嘉佑等，2015)。贝叶斯融合理论是一种概率方法，能利用动力模型集合平均预报更新气候(先验)分布，从而获得条件后验分布。该方法具有两个优点：订正和融合，订正指通过似然估计来订正耦合气候模型的模拟值；融合是指模型预报技能较低时，可以从历史气候态中获得与之相似的先验信息(Coelho，et al，2004；Luo，et al，2007)。Epstein(1962, 1985)讨论了贝叶斯方法应用于气象学的可行性并且已经应用于气候的统计推断和预测。Rajagopalan等(2002)使用贝叶斯方法在两个不同季节最佳地融合了全球季节降水和气温预测，并发现对单个模型的预测技能有普遍的提高。中国的一些学者也利用贝叶斯方法做了很多研究，潘旸等(2015)基于贝叶斯方法研究了实时观测的地面、卫星、雷达3种逐时降水资料的融合，结果表明该方法生成的产品精度优于任一单一来源的降水资料。郎杨(2015)利用贝叶斯融合方法对NCEP CFSv2在中国区域的季节气温预报进行订正，表明在每个季节和区域，预报技能都得到明显改善。

全球气候变暖是当前广泛关注的一个热门问题。为了确定中国气候变暖的趋势，就要建立一个有代表性的中国气温序列(王绍武等，1998)。中国最早的气象观测始于19世纪40年代，但直至20世纪初，观测站点仍然十分有限，特别是在西部地区。目前中国国家级业务单位馆藏的气温观测资料最早始于20世纪30年代，因此，在计算百年尺度中国气候变化序列或者估计近百年气候变化平均趋势时，往往存在着较大的不确定性(Li，et al, 2010, 2017；翟盘茂等，2016)。确保气候序列的均一性、降低气候资料中的系统性偏差(如城市化/土地利用等导致的偏差等)和尽可能保证资料分布的覆盖度等是降低气候变化趋势检测的不确定性的主要手段。本研究基于均一化的器测气温(李庆祥，2011)、气温代用资料与CMIP5模式模拟气温资料，利用统计后处理方法—贝叶斯融合模型开展中国北疆地区近百年气温序列融合试验，同时也利用观测资料进行了融合效果检验，探讨该方法对不同来源的多源气候序列进行融合的适用性。

2 资料和预处理方法

在融合试验中，使用中国新疆及周边区域气温观测序列、新疆周边区域气温代用资料序列以及CMIP5模式模拟序列，融合时选取1900年后气温重建序列与CMIP5模式序列，以1951年后资料序列进行效果检验。树轮气温重建序列来自中国气象局乌鲁木齐沙漠气象研究所，一共8条，分别为乌鲁木齐河(大西沟站)上年12月—当年3月平均最低气温重建序列(袁玉江等，1999)、乌鲁木齐河(大西沟站)7月平均最高气温重建序列(陈峰等，2011)、天山中段巩乃斯地区5—8月平均最高气温重建序列(张同文等，2011)、新源站4—8月平均气温重建序列(喻树龙，2010)、阿勒泰西部5—9月平均气温重建序列(张同文等，2008)、阿勒泰西部7月平均气温重建序列(张瑞波，2014)、新疆阿勒泰站6月平均气温重建序列(尚华明等，2010)、阿勒泰东部6—7月平均气温重建序列(胡义成，2012)。器测站点资料来自中国国家气象信息中心，所有资料均已进行均一化处理(Li，et al，2010；Xu，et al, 2017)，数据集的时段虽然是从1900年开始，但实际上，中国西部的器测站点气温资料最早是自20世纪30年代起。在计算过程中，均使用1961—1990年的气候平均值作为标准气候值计算距平序列。首先采用双线性插值将经过均一化处理的气温观测序列插值到图 1所示红色矩形内斜线填充的网格。

以上述气温代用资料序列为自变量，器测资料的网格序列为因变量，利用线性回归分析方法建立各网格的重建年气温序列，在文中命名为“网格重建序列”，并利用留一验证法对网格重建序列进行交叉检验。留一验证法的具体步骤为：在校准期内，剔除某一年的值，建立回归方程，然后将所剔除年代的重建值代入回归方程，得到对应年份气候要素的估计值。重复上述过程，直到得到整个重叠时段气候要素的估计序列，并同器测值进行比较，从而来检验回归方程的稳定性。图 1中红色矩形内为重建区域，包含15个网格(文用该网格中心点经纬度表示)，由于器测资料中网格(48.75°N，83.75°E)资料缺失，一共重建了14个网格的序列，其中，网格(48.75°N，81.25°E)、(43.75°N，83.75°E)、(46.25°N，83.75°E)、(43.75°N，86.25°E)、(46.25°N，86.25°E)、(48.75°N，88.75°E)、(48.75°N，91.25°E)未通过上述留一交叉检验，予以剔除，选取通过检验的7个网格参与融合试验，如图 1中反斜线所填充的网格。图 2为网格(43.75°N，81.25°E)经过线性回归所得的年重建序列与观测序列的对比，可见网格重建序列与原始观测序列的变化趋势基本一致，只是重建序列的趋势略偏缓，两者的线性趋势均通过了α=0.05显著性水平的t检验，两条序列的相关系数达到0.71，重建序列解释了观测序列50%左右的方差。

选取CMIP5模式41个成员1900年以来的模拟结果，将各个成员模拟资料采用双线性插值到如图 1所示的网格，为了客观衡量CMIP5模式各成员对气温的模拟能力，参考姜大膀等(2013)，绘制出1951—2005年各个模式成员与网格观测序列的泰勒图(图 3)，可以看出41个气候模式在模拟1951—2005年气温序列时的性能差异，模式序列与参考序列之间的距离越接近，即模式序列与观测序列的均方根误差越小，表明相应模式的模拟性能越好。可见离参考序列最近的8个模式序列分别为ACCESS1.0、BCC-CSM1.1、BCC-CSM1.1(m)、CNRM-CM5、CESM1(CAM5.1, FV2)、GFDL-ESM2M、MPI-ESM-MR和NorESM1-M。值得指出的是，尽管如此，单个模式对中国西部地区的模拟总体效果并不理想，相关系数基本在0.5以下，因此，本研究将这8个模式的平均序列进行融合，并深入分析这8个模式的平均序列与观测序列的泰勒图参数，结果发现其对1951—2005年气温序列的模拟性能有了较大的提升，相关系数达到了0.67，均方根误差也减小至0.71℃。选取上述8个模式进行融合试验(表 1)，相关系数均通过了α=0.05的显著性t检验。

3 贝叶斯融合方法

贝叶斯统计学认为，参数取值是随机的、经验的或者数据的信息形成参数的先验概率分布。在得到样本之前，从先验分布中产生了某个样本，即参数值。利用随机试验中不断获得新的信息修正先验分布，得到更为接近实际(更准确)的后验概率分布，进而降低期望损失。文中采用贝叶斯理论构建基于概率统计方法来融合各种不同来源信息的思路，可利用动力模式集合预报更新气候信息(先验信息)以获得条件后验分布。所获得的新的概率分布反映变量新的可信度水平。假设一随机变量气温θ，在抽样之前，关于θ的认知仅为历史资料的概率分布，即概率密度函数p(θ)，如果没有任何其他新的数据来源，例如气候模式，p(θ)即为从θ历史记录中得到的气候分布，现假设有一个气候模式，能产生该变量的预测值y，当给定新预测值y时，条件分布p(θ|y)则反映θ新的可信度，按照贝叶斯理论，则

(1)

式中，θ_t为时间t时的观测变量，y_t为同一时间动力模式成员集合平均值。p(θ)为(相对于在获取新信息y之前)θ的先验分布，p(θ|y)则称为θ的后验分布，p(y|θ)称为似然函数。贝叶斯方法有3个主要的步骤：先验分布的选择、似然函数的建模以及后验分布的确定。

3.1 先验分布选择

如前所述，先验分布一般是经验的或者数据的，但也可以有多种选择，如气候预测值或统计模型预测值(黄嘉佑等，2015；Coelho，et al, 2004；Luo，et al, 2007)。文中，将与观测数据最接近的数据作为先验分布的原则，以通过线性回归计算的网格重建序列作为先验分布均值，而利用留一验证法计算重建序列的方差作为先验分布方差。

3.2 似然函数

似然函数可以通过重建序列与模式模拟序列的关系估算。文中，对同一个网格，CMIP5模式有8个模式成员，采用模式成员平均值的条件分布p(y|θ)估算似然函数p(y|θ)，当模式成员平均值为y时，模式集合的条件分布为

(2)

试验中，采取线性回归模型构建网格重建值与模式平均值的似然函数

(3)

式中，α和β分别是截距和斜率。变量ε为回归残差(均值为0)并且假定为正态分布，其方差φ_ε反映回归的有效性。在线性模型下，y遵循正态分布

(4)

式中，α+βθ为平均值，φ_ε为残差的方差。假设模式模拟为正态分布，方差为φ_y，当给定模式模拟平均值时，y的条件分布为

(5)

假定模式模拟方差φ_y独立于线性回归误差φ_ε，则似然函数为

(6)

3.3 后验分布

从贝叶斯理论可知，当先验分布和似然函数均遵循正态分布时，后验分布p(θ|y)也遵循正态分布

(7)

式中，θ_p和φ_p分别为后验分布的平均值和方差，且

(8)

式中，θ₀和φ₀分别为先验分布均值和先验分布方差，可见与先验分布相比，后验分布的方差更小，这表明当先验信息与模式信息融合后，变量的不确定性将降低。

3.4 融合试验及与器测序列的对比

为了体现两种资料融合效果，这里以器测资料作为对比标准，选取20世纪50年代以来的时段开展融合试验。图 4为网格(48.75°N，86.25°E)1951年后观测、融合、重建和模式平均气温距平序列变化趋势，可见重建序列和模式序列的距平序列虽均为上升趋势，但二者的变化趋势值均较小，融合气温距平的上升趋势明显高于上述二者，与观测气温距平的趋势最接近且通过了α=0.05显著性水平的t检验。并进一步可以看出，在前期(20世纪80年代之前)，重建与器测资料吻合较好，而模式模拟资料效果差，φ_ε偏大，融合序列中模式资料的权重偏小，融合与重建资料的符合度也较好；此后重建相对于器测资料出现了明显的低估，但而后的模式模拟效果往往较好(与气候场相比)，φ_ε偏小，此时融合资料中模式权重增加，则较好地避免了重建资料的这种低估，这显然是一般简单的算术平均或权重平均所无法达到的。这也进一步说明，经过贝叶斯融合处理以后，与重建(先验分布)或模式平均气温距平相比，后验分布(即融合序列)与观测气温趋势更一致。

3.5 误差评估

融合过程中，遵循留一验证原则，所有后验分布参数均使用除本年外的剩余序列计算而得。融合结果使用均方根误差(RMSE)和平均绝对误差(MAD)进行评估

(9)

(10)

式中，N为样本的总数，T_i为网格观测气温，M_i为网格重建气温、网格融合气温或模式平均气温。

4 融合结果 4.1 融合精度评估

图 5和6分别为网格(43.75°N，81.25°E)、(48.75°N，86.25°E)观测、融合、重建和模式平均气温距平1900—1995年和1900—2003年的时间序列，可以看出，1951年以前，网格(43.75°N，81.25°E)的融合与采用线性回归的重建气温距平序列较一致，有较明显的升温趋势，模式平均气温距平序列比较平缓，升温不明显；1951年以后，观测、融合、采用线性回归的重建、模式平均气温距平序列的变化趋势基本一致，模式平均气温距平序列在1975年之前与其他3种序列相差较大，1985年以后与其他3种序列较一致，在整个期间，模式气温距平的变化比较平缓，网格(48.75°N，86.25°E)也有相似的情况。

进一步，对成功重建的7个网格序列和8个气候模式模拟的平均序列进行了融合。从各网格融合、重建、模式平均气温距平序列与观测气温距平序列的均方根误差和平均绝对误差亦能表达同样的结果。从表 2可知，融合后的距平序列与重建距平序列相比，各个网格的均方根误差和平均绝对误差均有所减小，不确定性降低。两者相差较大的网格为(46.25°N，88.75°E)，均方根误差和平均绝对误差分别相差0.19和0.15℃；与模式平均距平序列相比，融合后的距平序列的均方根误差和平均绝对误差也均小于模式平均距平的均方根误差和平均绝对误差，相差最大的为网格(43.75°N，81.25°E)，分别为0.18和0.11℃。这说明，通过贝叶斯融合以后的序列误差更小，精度更高。

4.2 融合序列长期趋势

鉴于CMIP各模式初始化时存在的差异，考察模式(特别是单个模式)序列对各个具体年份的模拟往往没有实际意义，而通过对多个模式的集合后，多模式序列在年代际尺度或者长期趋势才具有较好的代表性。因此，文中也比较了融合序列在长期趋势上的改善(即不同时期的相对差异)。如表 3所示，参与融合的7个网格，重建气温和模式气温的距平序列虽均为上升趋势，但二者的趋势值均较小，而从表中融合序列可以看出，融合气温距平的上升趋势均比上述二者明显，与观测气温距平的趋势最接近且均通过了α=0.05显著性水平的t检验；表 4为1900年以来各个网格融合气温、重建气温以及模式平均气温距平线性变化，可见各个网格3种气温距平均为上升趋势，除网格(46.25°N，81.25°E)和(46.25°N，91.25°E)融合气温和模式平均气温距平一致外，其他5个网格，融合气温距平比其他两种气温距平的趋势上升明显且均通过了α=0.05显著性水平的t检验。

5 结论与讨论

本研究尝试将贝叶斯方法应用于新疆北部地区若干个网格的树轮气温重建序列、CIMP5模式序列的融合，并从序列线性趋势及统计分析的角度检验了该方法的有效性。结果表明：

(1) 简单的线性回归所产生的重建序列与模式模拟序列均包含气候趋势变化的信息，在贝叶斯融合中，通过似然函数有效地去除了模型预测的偏差，同时分别提取先验分布和模式模拟的有用信息进行融合以产生最小方差的后验分布，与任一单一序列相比，融合产生的后验序列误差更小，精度更高，线性趋势也更接近于观测序列。

(2) 贝叶斯融合能够容易地应用于先验信息与气候模式序列的融合，其巧妙地结合了统计模型和动力模型各自的优点，产生更具代表性的融合结果，融合结果与先验分布及模式模拟性能有着密切的关系，融合结果能够纠正先验分布及气候模拟数据的偏差。当存在气候模式序列或其他气候序列如树轮重建序列而缺乏观测序列时，使用贝叶斯融合方法创建气候替代序列进行气候分析不失为一个可行的途径。

需要说明的是，本研究旨在探讨贝叶斯融合方法在多来源气候序列融合中应用的可行性，寻求当无仪器观测序列时，是否可找到尽可能精确的替代序列的方法。在融合时选取的信息源相对较少(仅仅是多模式平均和古气候代用资料两种)；此外，树轮重建资料及模式资料本身的质量水平和精度也直接影响最后融合资料的质量，毕竟如果没有高质量的先验分布，不管融合手段如何先进，仍然难以得到很好的融合结果；CMIP5模式的模拟能力对融合结果也至关重要，特别在20世纪50年代之前，CMIP5模式对北疆地区气温序列模拟效果本身较差，融合序列基本与先验分布(即重建序列)重合，而此后，模式模拟效果逐渐转好，才进一步体现出融合序列的优越性。因此，融合分量的选择也非常重要，本研究仅是一个初步尝试。在实际应用中，还需要进一步根据实际数据特点，选择最好的融合分量，设计最佳融合方案，以获取更佳融合效果。