1 引言

蒸散发量(Evapotranspiration)是植被及地面向大气输送的水汽总通量，区域蒸散发量包括地表土壤蒸发和植被蒸腾两部分，其中土壤蒸发是土壤中的水分通过上升和汽化从土壤表面进入大气的过程，受太阳辐射、湿度、温度、风速等一系列气候因素的影响，而植被蒸腾是植被体内水分散失的主要原因，与植被的光合作用和呼吸作用有关(Chen et al., 2003)。目前，借助地面仪器观测和遥感估算是获取蒸散发量的有效手段(Wang and Dickinson, 2012；Liu et at.，2013)。Anderson et al.(2012)利用遥感获取的地表温度和能量平衡模型，估算了区域面上的蒸散发量。Liu et al.(2013)利用涡动相关系统和大孔径闪烁仪，分析了海河流域蒸散发量的季节和年变化趋势。此外，模式和再分析资料也可以提供单点及区域面上蒸散发的数据集(Roads and Chen, 2001；Ma and Ma, 2015)。

黄河源区地处青藏高原的东北部，近年来该区域主河道断流、水土流失、草地生态系统持续退化、土地沙化及荒漠化日趋严重，该区域的蒸散发量一直是国内外关注的热点。Li et al.(2013)基于能量平衡模型和遥感资料，估算了黄河源区湿地下垫面的蒸散发量。Liu et al.(2015)利用静止气象卫星资料和遥感模型估算了黄河源区蒸散发量的日变化过程。以上模型大都针对日尺度进行评估，而分析长时间序列蒸散发量的变化情况，对了解区域内水分和能量的空间分配，区域水资源的可持续开发利用和生态保护都具有重要的意义。周德刚和黄荣辉(2006)通过地面台站的观测资料研究表明，近50年来黄河源区气温呈持续上升的趋势，降水和水文断面径流量呈减少的趋势，但由于观测手段和时间的限制，作为水分循环中重要因子——蒸散发量的变化，还有待商榷。再分析资料可以提供较长时间尺度的区域蒸散发产品，有潜力为黄河源区水文研究提供基础的数据支持(王霁吟等，2015)。目前，利用再分析资料提供的蒸散发数据集进行黄河源区蒸散发量的评估工作并不常见。

本文选取的再分析资料为：欧洲中期天气预报中心(ECMWF)提供的ERA-interim和美国国家环境预测中心(NCEP)和大气研究中心(NCAR)提供的再分析资料数据集，但两种再分析资料地处青藏高原东北部的黄河源区的可信度和质量亟待评估。因此，本文首先基于涡动相关系统观测的2010年全年的蒸散发资料，评估了两种再分析资料在黄河源区的适用情况。其次，基于统计方法分析了源区蒸散发量的变化趋势。最后分析了黄河源区蒸散发量突变的年份，以期对黄河源区水循环和水量平衡提供可信的蒸散发资料。

2 研究方法

利用两种蒸散发再分析资料，基于2010年涡动相关系统观测的蒸散发结果进行分析，得出两种蒸散发再分析资料在黄河源区的适用情况；利用MK(Mann Kendall)方法和Sen斜率检测法，分析黄河源区长时间序列蒸散发量的变化趋势和速率；利用滑动t检验和SQMK(Sequential Mann Kendall)方法，检测黄河源区蒸散发量发生突变的时间和区域。

2.1 趋势分析方法

MK检验和Sen斜率检测法是基于秩和变化斜率的检验方法，其优点在于不需要预先假定样本的分布，不受异常值的干扰，能很好地处理非正态分布的数据。因此，这2种方法被广泛地用于水文、气象数据序列变化趋势的研究。

2.1.1 MK方法

MK方法是一种非参数统计检验方法，变量可以不具有正态分布特征，因此适用于水文变量的趋势检验(张建云等，2007)。假设X₁, X₂, …, X_n为时间序列变量，n为时间序列的长度，MK方法定义统计量S为

$S = \sum\limits_{j = 1}^{n - 1} {\sum\limits_{k = j + 1}^n {s({X_k} - {X_j})} } $

(1)

其中，

${\text{s}}({X_k} - {X_j}) = \left\{ \begin{gathered} + 1, {\text{ }}{X_k} > {X_j} \hfill \\ 0, {\text{ }}{X_k} = {X_j} \hfill \\ - 1, {\text{ }}{X_k} < {X_j} \hfill \\ \end{gathered} \right.$

(2)

式中X_j、X_k分别为j、k年的相应测量值，且k > j。

$Z = \left\{ \begin{gathered} \frac{{S - 1}}{{\sqrt {V(S)} }}, {\text{ }}S > 0 \hfill \\ 0, {\text{ }}S = 0 \hfill \\ \frac{{S + 1}}{{\sqrt {V(S)} }}, {\text{ }}S < 0 \hfill \\ \end{gathered} \right.$

(3)

其中，Z为一个正态分布的统计量，V(S)为方差。在给定的α置信水平上，如果$\left| Z \right| \geqslant {Z_{1 -a/2}}$，即在α置信水平上，时间序列数据存在明显的上升或下降趋势。$\left| Z \right| \geqslant {\text{1}}{\text{.96}}$和$\left| Z \right| \geqslant 2.58$分别表示其通过了置信度为95%和99%的显著性检验。

2.1.2 Sen斜率检测法

Sen斜率检测法能降低或避免数据缺失及异常对统计结果的影响，以样本在不同长度的变化率构造秩序列，基于一定显著性水平α进行统计量检验，得出变化率数值区间，并以中值大小判断时间序列变化趋势及程度(Gocic and Trajkovic, 2013)。Sen斜率公式为

${S_{ij}} = {M_{}}\left[ {\frac{{{X_j} - {X_i}}}{{j - i}}} \right]$

(4)

其中，X_i和X_j分别为第i和第j时刻的序列值，1 < i < j < n，n为序列长度。Sen斜率即斜率中值，由序列长度n决定的S_ij总数N的奇偶性决定，k为与序列长度有关的整数，$N=n(n -1)$：

$S = \left\{ \begin{gathered} {S_{k + 1}}, {\text{ }}N = 2k + 1 \hfill \\ ({S_k} + {S_{k + 1}})/2, {\text{ }}N = 2k \hfill \\ \end{gathered} \right.$

(5)

2.2 突变检测方法 2.2.1 SQMK方法

用SQMK方法检验气候序列突变时，突变检验统计量计算为：由X₁, X₂, …, X_n构造一秩序列m_i (m₁, m₂, …, m_n)，m_i为X_i＞X_j(1≤j≤i)的样本累计数(祁威等，2013)。定义变量d_k为

${d_k} = \sum\limits_i^k {{m_i}, {\text{ }}2 \leqslant k \leqslant N} $

(6)

d_k平均值以及方差定义如下：

$E({d_k}) = \frac{{k(k - 1)}}{4}$

(7)

$V({d_k}) = \frac{{k(k - 1)(2k + 5)}}{{72}}$

(8)

在时间序列随机独立假定下，定义U_k(统计量)为

${U_k} = \frac{{{d_k} - E({d_k})}}{{\sqrt {V({d_k})} }}$

(9)

SQMK方法突变检验中，在正序列曲线(U_k)超过临界值信度线的前提下，若正序列和反序列两个统计量序列仅有一个明显的交叉点，且位于信度线之间，则表明该交叉点为突变点；若交叉点位于信度线之外，或者是存在多个明显的交叉点，则不确定是否为突变点。

2.2.2 滑动t检验方法

在滑动t检验方法中，设变异点前后两序列总体分布函数各为F₁(x)和F₂(x)，从总体F₁(x)和F₂(x)中分别抽取n₁和n₂两个样本，X₁和X₂为样本均值，${S_1}^2$和${S_2}^2$分别为方差(余予等，2013)。突变检验统计量T为

$T=\frac{\overline{{{X}_{1}}}-\overline{{{X}_{2}}}}{\sqrt{{\frac{({{n}_{1}}-1)s_{1}^{2}+({{n}_{2}}-1)s_{2}^{2}}{{{n}_{1}}+{{n}_{2}}-2}}/{(1/{{n}_{1}}+1/{{n}_{2}})}\;}}$

(10)

T遵从自由度$v={n_1} + {n_2} -2$的t分布。利用滑动t检验法，连续设置基准点，可以得到统计量序列T_i，给定显著性水平(α)，查t分布表得到临界值(q_α)，若$\left| {{q_i}} \right| > \left| {{q_\alpha }} \right|$，则认为在基准点时刻出现了突变。

3 研究区域和数据来源

黄河源区(31.5°N～36.5°N，95.5°E～103.5°E)位于青藏高原东北部，面积约12.2×10⁴km²，横跨青海、四川、甘肃三省，绝大部分在海拔3000～5000 m以上，该地区陆面蒸散发量的大小与黄河水资源变化有密切的关系(康悦等，2011)。图 1为黄河源区在中国的位置及黄河流经黄河源区的主要站点分布情况，其中五角星为架设涡动相关系统的站点。黄河源区内植被类型以高原灌丛草甸为主，间以零星灌木丛。左上角两块蓝色区域为扎陵湖和鄂陵湖，布尔汗不达山脉以北为柴达木盆地，为灌木丛覆盖。

本文以中国科学院寒区旱区环境与工程研究所若尔盖高原湿地试验站以及源头麻多乡陆面过程试验站为依托，收集了该地区2010年1～12月逐日地面气象观测资料和涡动相关系统观测资料。涡动相关系统观测资料包括辐射通量、感热通量、潜热通量、土壤热通量等。观测点所在位置较为平坦，下垫面土壤类型为高山草甸土，植被类型单一，以10 cm高的草甸覆盖为主。涡动相关系统架设高度3.02 m，观测资料包括短波辐射和长波辐射通量(向上和向下)、感热通量、潜热通量、土壤热通量(2层：-10 cm和-20 cm)等，采用频率10 Hz，参照《涡动相关通量观测指导手册》(王介民，2012)进行处理，每10 min输出一次该时段平均通量，最后将结果累加，得到日尺度的平均通量值。

ERA-interim所采用的同化系统是ECMWF集成预报系统，其中包括四维空间变量分析的运用、对卫星资料变量偏差的修正以及对观测系统中其他数据的处理等。本研究所选取的ERA-interim蒸散发再分析资料空间分辨率为0.125°(纬度)×0.125°(经度)，变量为日蒸散总量，获取网址为http://apps.ecmwf.int/datasets/[2015-02-01]。

NCEP-DOE 2再分析资料校正了NCEP/NCAR中存在的一些已知误差问题，增加了对陆地降水的简单同化，而且修改了数值模式中的一些物理过程和参数化方案，因而被看作是NCEP/NCAR再分析资料的更新和订正版本。蒸散发再分析资料空间分辨率约为2.0°(纬度)×2.0°(经度)，变量为日蒸散总量，获取网址为http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.gaussian.html[2015-02-01]。

4 黄河源区蒸散发量时空变化趋势及突变分析 4.1 两种蒸散发再分析资料的相关性分析

再分析资料是利用资料同化技术把数值预报产品与观测资料融合起来的产物，包含有数值模式、同化方案和观测系统变更等引入的误差。因此，再分析资料的质量检验和评估不仅有助于提高气候变化研究结果的可靠性，而且对于气候模式评估和资料同化技术的改进等也有十分重要的意义。本研究首先将ERA-interim和NCEP-DOE 2提供的蒸散发再分析资料处理为空间分辨率为2.0°(纬度)×2.0°(经度)的格式，并利用平均汽化潜热系数(2.49×10⁶W m^-2mm^-1)求解出蒸散发量，使再分析资料具有统一的空间分辨率和单位；其次，由于ERA-interim和NCEP-DOE 2提供的蒸散发再分析资料为日尺度，需要将涡动相关系统观测的地面蒸散发资料进行整理，将每10 min一次的观测资料累加处理成日资料，缺失数据和有降水过程发生时刻的数据进行剔除，所以再分析资料和观测结果在时间尺度上的误差是造成对比结果差异的最可能原因。Liu et al.(2011)分析了不同风向下涡动相关系统通量源区的情况，探讨了单点尺度与区域尺度通量观测尺度的关系。由于本研究观测点所在位置较为平坦，下垫面比较单一，根据湍流观测中的印痕(footprint)效应，涡动相关系统观测的量不能认为是个点上的值，而具有几十米到几公里的空间代表性。

图 2是利用2010年1～12月黄河源涡动相关系统观测的蒸散发量与两种蒸散发再分析产品的对比情况，左图为时间序列的对比图，右图为散点图。从左图可以得出，ERA-interim和NCEP-DOE 2两种蒸散发再分析产品都呈现较明显的月变化分布特征，即1～3月值偏低，4月开始逐渐升高，7月呈现全年蒸散发的最大值，6月中旬到7月中旬有轻微的减少，7月中旬到8月中旬增大，8月中旬到12月缓慢减少。其中，NCEP-DOE 2提供的蒸散发再分析资料在4～6月以及10～12月较地面观测资料有较明显的偏大。由右图的散点图也可以得出，NCEP-DOE 2提供的蒸散发再分析资料在1:1等值线附近更为离散。表 1为两种再分析资料与地面观测资料的误差分析情况，其中回归方程的斜率越接近1，误差平方和越小，相关系数越接近1，表示测量值(再分析资料)与真实值(地面观测值)越逼近，均方根误差表示数据本身的离散程度。从表 1中可以得出，利用ERA-interim提供的蒸散发再分析资料与地面观测资料对比，回归方程的斜率为0.8827，误差平方和为147.8，相关系数为0.8667，均方根误差为0.6381，ERA-interim较NCEP-DOE 2提供的蒸散发再分析资料更适用于黄河源区以高原灌丛草甸为主的下垫面。

4.2 黄河源区蒸散发量时空变化趋势

图 3和图 4分别为利用MK方法和Sen斜率检测法获取的黄河源蒸散发时空变化趋势，其中图 3a、图 4a为年变化趋势，另外还选取了全年中蒸散发变化最为剧烈的3个时段进行分析，分别为5～6月、7～8月及9～10月。

利用MK方法统计1979～2014期间黄河源区蒸散发量的年变化趋势可知(图 3a)，黄河源区蒸散发量总体呈现北高、西南低、东南变化不明显的趋势，北部兴海—共和—贵德地区变化最为显著，呈现明显增高的趋势，是全区蒸发量变化最明显的区域，大部分区域超过了99%的信度检验；西南部曲麻莱—治多—玉树地区蒸散发量明显减少，介于95%～99%的置信区间；东南部玛沁—玛曲—久治地区蒸散发量的变化并不明显，呈现缓慢上升的趋势，大部分区域没有超过95%的信度检验。图 3b、图 3c、图 3d分别显示了全年中蒸散发变化最为剧烈的3个时段的情况，其中5～6月的变化趋势和年变化趋势相近，呈现出现北高、西南低、东部变化不大的形态，但是大部分变化区域都没有超过99%的信度检验；西南部曲麻莱—治多—玉树地区在7～8月的变化趋势与年变化趋势不同，并没有呈现出明显减少的趋势，而是和东南部玛沁—玛曲—久治地区类似，呈现出变化不明显的特征；东南部的蒸散发量在年变化、5～6月以及7～8月都变化不明显，但在9～10月呈现出明显增加的趋势，可能原因是由于东南部玛沁—玛曲—久治地区大部分下垫面为季节性湿地覆盖，该区域季节性湿地在9～10月比较活跃，影响了区域蒸散发量的大小。

图 4为利用Sen斜率检测法得到的1979～2014年黄河源区蒸散发变化情况，利用该方法计算的年变化结果(图 4a)与MK方法(图 3a)相似，都呈现北高、西南低、东南变化不明显的趋势，北部兴海—共和—贵德地区呈现明显增高的趋势，在1.5～2.5 mm/a，是全区蒸发量变化最明显的区域；西南部曲麻莱—治多—玉树地区蒸散发量明显减少，在－1.0～－0.5 mm/a；东南部玛沁—玛曲—久治地区蒸散发量的变化并不明显，呈现缓慢上升的趋势，蒸散发量的变化在0.5～1.0 mm/a。图 4b、图 4c、图 4d分别显示了全年中蒸散发变化最为剧烈的3个时段的情况，其中5～6月的变化趋势和年变化趋势相近，呈现北高、西南低、东部变化不大的形态；北部兴海—共和—贵德地区在7～8月呈现明显增高的趋势，大于1.0 mm/a，是季节变化最明显的区域；东南部的蒸散发量在年变化、5～6月以及7～8月变化都不明显，但在9～10月呈现出明显增加的趋势，大于1.0 mm/a。

4.3 黄河源区蒸散发量突变情况

图 5和图 6分别为利用SQMK和滑动t检验的方法检测若尔盖高原湿地试验站蒸散发变化的突变年，其中图 5a、图 6a为年变化趋势，另外还选取了全年中蒸散发变化最为剧烈的3个时段进行分析，分别为5～6月、7～8月及9～10月。

SQMK方法中正序和逆序的交点且位于信度线之间，就是蒸发突变年，若交叉点位于信度线之外，则不确定是否为突变年。由图 5a可知，黄河源区若尔盖高原湿地试验站的蒸散发量在1979～2014年期间有3个突变年，突变在1989年，后持上升趋势，1995年后上升趋势显著，正序线和逆序线在1984和1989年有交点且趋势上升，说明1984年和1989年的蒸发量有明显增加；图 5b显示突变年份又较图 5d晚，且上升趋势并不显著；图 5d和全年最为相似，但突变年份为1990年，较全年蒸散稍晚，且1995年上升趋势显著。蒸散发变化最剧烈的3个时段中，5～6月和9～10月的突变年较少，都集中在20世纪80年代和90年代，其中5～6月蒸发突变年为1990年和1993年，9～10月蒸发突变年为1984、1987和1990年，而7～8月的交点比较多，依照MK方法规则，由于其存在多个明显的交叉点，则不确定是否为突变点。

SQMK方法在气象水文领域统计突变点应用广泛，但是由于该方法检测突变点存在局限性，因此辅以滑动t检验法来检测黄河源若尔盖地区蒸散发序列的突变情况。图 6为利用滑动t方法检测黄河源区1979～2014年蒸散发突变情况，上下两条黑线内为置信区间，超过绝对值2.31表示通过95%的信度检验，年变化情况如6a所示，利用滑动t检验法1989年检测到显著的蒸散发突变，即1989年为若尔盖地区的蒸发突变年，并且趋势向上，且1990年9～10月的蒸散发量同时发生了突变，时间上略晚于全年蒸散发量突变时间，这与根据SQMK方法得出结论是相同的，不同的是5～6月蒸散发在1985年发生了明显突变，而SQMK并无明显趋势与之对应，通过上述两种方法的分析，基本上可以确定在1990年左右，黄河源区的蒸散发量发生了明显的突变，趋势为增大，而其中秋季蒸散趋势的响应全年蒸散增大最为显著，而7～8月的蒸散发变化并不显著，而5～6月蒸散发突变年份可能在1985年，趋势与全年变化相同，并且时间大幅度提前于全年变化趋势，但变化量并不显著。利用该方法检测蒸散发变化最为剧烈的3个时段的结果显示，5～6月的突变年为1985年，7～8月期间没有检测出蒸发突变年，而1990年为9～10月的蒸发突变年。

为了更深入的了解整个黄河源区蒸散发突变的情况，表 2统计了沿黄河一线中不同站点发生突变的情况，其中，麻多X和麻多代表黄河源头区，达日、久治、玛曲代表黄河源的主要产流区，唐乃亥代表黄河在黄河源区的出口区。6个站点沿黄河一线依次排列，海拔的分布由高(麻多X，4450 m)到低(唐乃亥，297 0m)；年平均蒸散发量除了唐乃亥之外呈现由低到高的次序排列，其中麻多X和玛曲年平均蒸散发量分别为364.45 mm和493.32 mm，是黄河源蒸散发量最低和最高的区域；变差系数(Coefficient of Variation，CV)是一个表示标准差相对于平均数大小的相对量，统计年度最大蒸散月的蒸散量在1979～2014年期间的变率(CV of max monthly ET)可知，久治(0.087)和玛曲(0.083)的蒸散发变化最剧烈；此外，表 2还统计了6个站点蒸散发在年尺度和3个不同时段的突变情况。以年变化而言，除了久治在1998年检测到突变之外，其他各站点的突变都集中在1984～1989年间，特别是1989年发生了明显的突变。针对久治的情况，5～10月蒸散发量在1998年并未发生明显突变，由此可以推断可能是由于非植被生长期蒸散发量发生了突变。而其他各站点在20世纪80年代末期至90年代初期变化较为活跃，而源头区除此时期外，进入20世纪后，也有明显的蒸散发量减小的趋势。以蒸散发变化最为剧烈的3个时段而言，除了5～6月在麻多X、玛曲和达日在2000年之后有突变之外，各站点的突变都集中在20世纪80世纪和90年代。

5 结论

本研究基于ERA-interim和NCEP-DOE 2提供的蒸散发再分析资料，与若尔盖高原湿地试验站地面实测资料进行对比，通过分析相关系数和均方根误差，得出ERA-interim较NCEP-DOE 2提供的蒸散发再分析资料更适用于黄河源区以高原灌丛草甸为主的下垫面。

利用ERA-interim提供的蒸散发再分析资料，基于MK和Sen斜率检测法评估了黄河源区蒸散发量的变化趋势。1979～2014年黄河源区蒸散发量总体上呈现北高南低的年变化趋势，北部兴海—共和—贵德地区增加最为迅速，在1.5～2.5 mm/a，西南部曲麻莱—治多—玉树地区减少最为明显，在－1.0～－0.5 mm/a，东南部玛沁—玛曲—久治地区蒸散发量的变化在0. 5～1.0 mm/a。5～6月的变化趋势和年变化趋势相近，呈现北高、西南低、东部变化不大的形态。

利用ERA-interim提供的蒸散发再分析资料，基于SQMK和滑动t检验方法评估了黄河源区蒸散发量的突变情况，得出1990年左右，黄河源若尔盖地区的蒸散发量发生了明显的突变，趋势为增大。此外，还统计了黄河源区其他6个站点蒸散发的突变情况。除了久治在1998年检测到突变之外，其他各站点的突变都集中在1984～1989年间，特别是1989年发生了明显的突变。