若尔盖湿地参考作物蒸散量时空演变特征及成因分析

梁 霄, 靳晓言, 强皓凡

(1.四川大学水力学与山区河流开发保护国家重点实验室 水利水电学院,成都 610065;2.南方丘区节水农业研究四川省重点实验室,成都 610066)

蒸散是水文循环中十分重要的环节，是影响地区水热平衡的重要气候因子和参数[1]。研究表明，每年蒸散进入大气的水分约占降水量的2/3，蒸散是当前研究陆面—大气水文循环过程的重要内容之一[2]。参考作物蒸散(Reference Crop Evapotranspiration，ET0)既是水分循环的重要部分，也是能量平衡的重要部分，表示在一定气象条件下水分供应不受限制时，某一固定下垫面能达到的最大蒸散量，也称潜在蒸散量，是评价某一地区干旱程度、研究作物需水及指导灌溉的重要因子[3]。目前估算ET0的方法大体分为：Penman-Monteith法、温度估计法、辐射能量法、质量传导法和水平衡法等[4]。其中，1998年由联合国粮农组织(FAO)推荐的Penman-Monteith模型综合了空气动力学和太阳辐射等关键因素的影响，具有良好的水文气象物理基础，在全球范围应用广泛[5]。

黄河源区地处青藏高原东北部，近年来该区域主河道断流、水土流失、草地生态系统持续退化、土地沙化及荒漠化日趋严重，该区域的蒸散一直是国内外研究热点[6]。其中，位于源区东缘的若尔盖湿地是黄河源区重要的水源涵养地，素有黄河“蓄水池”之称，每年供给黄河源超过30%的径流量[7]，同时是高原面积最大、最典型的高寒沼泽湿地，其主导功能是水源涵养，并具有径流调节、生物多样性保护、水土保持、沙化控制、调节局部区域小气候、环境自净及固碳等辅助生态功能[8]。近年来，由于人类活动和自然因素的双重作用，湿地逐渐退化，沙漠化面积与强度明显增加，环境不断恶化[9]。

目前专门针对若尔盖湿地参考作物蒸散的研究较少，且主要集中于ET0总体变化趋势与气象因素单因子分析[1, 10]，未对该区ET0时空分布与季节变化进行较为系统的研究。本文基于辐射矫正的Penman-Monteith模型[11]，采用累积距平、Mann-Kendall检验、Pettitt检验、Theil-Sen趋势度、EOF分解和通径分析等对1960-2015年若尔盖湿地ET0时空变化及影响因子进行分析。研究对湿地水资源的科学管理和脆弱生态环境恢复具有积极意义。

1 研究区域与数据

1.1 研究区概况

若尔盖湿地(32°20′～34°00′N，101°36′～103°30′E)地处青藏高原东缘，隶属四川省若尔盖县、红原县和阿坝县，以及甘肃省的玛曲县和碌曲县等，是我国特有的沼泽分布区，还是全球面积最大的高原湿地，生态研究意义重大。该区域属大陆性高原气候，寒冷湿润，霜冻期长，日温差大，平均海拔3 500 m，年均降水量750 mm，年平均气温1 ℃左右(数据由区域内各气象站点算术平均得来)[12]。

图1 若尔盖湿地及其周边气象站分布Fig.1 Distribution of meteorological stations of the Zoige Wetland and its surrounding area

1.2 数据来源

本文数据来源于中国气象数据共享服务网(http:∥cdc.cma.gov.cn/home.do)，包括若尔盖湿地及周边地区(图1)共19个气象站1960-2015年逐日气象资料，即降水量(P)、最高气温(Tmax)、最低气温(Tmin)、平均气温(Tmean)、相对湿度(RH)、10 m风速(U10)、大气压强(p)和日照时数(n)。

2 研究方法

2.1 Penman-Monteith模型

前人研究表明，P-M模型计算结果与黄河上游地区蒸发皿实测值之间存在良好的相关性，采用P-M模型计算ET0是可行的[13]。计算公式如下[14]：

式中：Δ为饱和水汽压曲线斜率，kPa/℃；γ为干湿计常数，kPa/℃；U2为2m高处风速，m/s；Rn为净辐射，MJ/m2；G为土壤热通量，MJ/m2；T为平均气温，℃；es和ea分别为饱和水汽压和实际水气压，kPa。除净辐射Rn应进行地区校正外[11]，各变量根据FAO方法计算[14]。本文采用Yin等辐射矫正的经验系数计算Rn[11]，公式如下：

式中：σ为Stefan-Boltzmann常量[4.903×10-9MJ/(K4·m2)]；n为实际日照时数，h；N为最大日照时数，h；Rso为晴天辐射，MJ/m2；Tmax,k、Tmin,k分别为最高和最低气温，K。

2.2 数据分析方法

本文采用Theil-Sen趋势度、累积距平和Mann-Kendall检验法分析ET0趋势变化，采用Pettitt检验[15]确定突变时间，对ET0变化的强度、置信度和趋势度进行反距离权重插值(IDW)并对年ET0进行EOF分解分析空间演变特征，最后采用通径分析研究ET0变化成因。

2.2.1Mann-Kendall趋势检验法

Mann-Kendall(M-K)趋势检验法是一种非参数统计检验，不需要样本服从一定的分布，同时也不受少数异常值的干扰，更适用于类型变量和顺序变量，适用于水文、气象等非正态分布的数据。本文时间序列数据( x1，x2，…，xn) 是样本容量n=56的随机变量，其正态分布的M-K统计值定义如下[16]：

2.2.2Theil-Sen趋势度

ET0的年际变化率由Theil-Sen趋势度来表示，该方法可以减少数据异常值的影响，是一种稳健的非参数统计趋势计算方法，其计算公式为[17]：

式中：1

2.2.3 经验正交函数(EOF)分解

经验正交函数(EOF)分解是气候变化领域常用的时空分解方法，由统计学家Pearson于1901年提出，并由Lorenz在20世纪50年代首次引入气象研究领域，具有展开和收敛速度快、资料信息可集中量大、能在有限区域内对不规则分布站点进行分解且所分解空间结构物理意义明确等优点[18,19]。本文建立湿地年ET0标准化矩阵，在MATLAB中计算得到相应特征值及贡献率，选取贡献率较大的前3个特征向量，利用IDW法进行空间插值进而分析湿地ET0的空间分布特征。

2.2.4 通径分析

通径分析是数量遗传学家SewallWright于1921年提出的一种多元统计技术，通过将自变量与因变量之间的相关系数分解为直接作用和间接作用的代数和，来研究变量间的相互关系、自变量对因变量的作用方式和影响程度[20]。该方法不受自变量间度量单位和变异程度的影响，从而能为研究结果提供可靠依据[21]。由于ET0受各气象因子综合影响，且各因子间存在较强的相关性，很难确定单个因子对ET0作用的程度[21]，通径分析方法能很好地解决这一问题。

3 结果与分析

3.1 若尔盖湿地ET0时间变化趋势

以各气象站点经辐射矫正的Penman-Monteith模型计算的ET0取平均作为区域代表值研究时间变化趋势，以准确、全面地表述若尔盖湿地的整体气候情况。由图2(a)，1960-2015年若尔盖湿地年ET0均值为625.3mm，这低于王建兵等[1](761.3mm)的研究结果，一方面是由计算年限及站点差异所致，另一方面是因为本文采用了辐射校正的P-M模型，而未经校正的计算结果普遍偏高[11]。其中最小值出现在1989年，为576.3mm，最大值出现在2006年，为674.8mm，极差为98.5mm；Theil-Sen趋势度表明，年ET0以4.89mm/10a的速率显著上升(MK检验Z=2.39，p<0.01)，这与杜家强等[13]研究结果一致。Pettitt突变检验表明，年ET0于1968年发生突变(p<0.01)，突变前(1960-1967年)多年平均ET0为598.1mm/a，突变后(1968-2015年)增至629.8mm/a，增幅为5.3%。气象数据分析表明，湿地平均气温和降水量分别于1967和1968年发生突变上升，可能由此导致了ET0的突变。另外，气象突变是某个区域上同步发生的事件，未来研究还需附近区域相比较来进一步印证这一结论。

累积距平法是常用的用曲线直观判断因素变化趋势的方法，同时可根据曲线变化情况判断因素变化的阶段性特征[18]。由图2(b)累积距平检验表明，湿地年ET0总体呈“低—高—低—高”变化趋势，4个阶段分别为1960-1968年、1969-1980年、1981-2000年和2001-2015年，各阶段年平均ET0为598.9、641.5、616.7、639.5mm。其中1968年既是转折点也是突变点(通过0.01信度检验)，1980年和2000年只是转折点而非突变点(未通过显著性检验)。

图2 若尔盖湿地年际ET0变化与累积距平曲线Fig.2 Annual ET0 and its accumulated departure of the Zoige Wetland

季节变化上(图3)，夏季ET0均值最大，为248.2mm，其次为春季188.5mm，秋季125.5mm和冬季63.2mm，依次约占全年ET0的39.69%、30.14%、20.07%和10.10%。趋势变化上，各季ET0均呈上升趋势，其中秋季上升最明显，速率为1.84mm/10a(Z=3.57，p<0.01)，其次依次为：夏季1.73mm/10a(Z=1.25，不显著)，春季0.96mm/10a(Z=1.16，不显著)，冬季0.73mm/10a(Z=1.45，p<0.1)。Pettitt突变检验表明，秋季ET0于1997年发生突变(p<0.01)，突变前年均ET0为123.3mm/a，突变后增至129.8mm/a，增幅为5.3%。冬季ET0于2003年发生突变(p<0.1)，突变前年均ET0为62.3mm/a，突变后增至66.1mm/a，增幅为6.2%。春、夏两季未出现突变。

图3 若尔盖湿地各季ET0及Pettitt突变检验Fig.3 Seasonal ET0 and its Pettitt test of the Zoige Wetland

3.2 若尔盖湿地ET0空间分布特征

由于自然地理环境的空间异质性，若尔盖湿地不同空间位置的ET0表现出不同大小与变化趋势。由于若尔盖湿地地处青藏高原边缘，周边气象站点分布稀疏且数据的时间序列长度往往较短，前人文章大多仅选取3～5个站点进行分析[1,12]。为增强本文研究对象的数据代表性，采用湿地及其周边气象数据连续完整的共19个站点，且这些站点较为均匀地分布在湿地的各个方位上，因此空间插值分析具有良好的自然地理代表性。由图4(a)，1960-2015年若尔盖湿地年均ET0空间分布在580.9～672.2mm之间变化，高值主要分布在湿地南部及东部边缘，低值主要分布在湿地中部的若尔盖(596.2mm)与红原(591.3mm)、西部的久治(580.9mm)和东南部的松潘(583.9mm)等地，整体呈南部、东部边缘高、西北—东南一线较低的空间分布特征。

利用MK检验分析ET0年际变化趋势[图4(b)]可知，湿地ET0变化总体呈东北高、西南低的趋势。阿坝—红原以北地区变化最为显著，呈明显增高趋势；松潘一带ET0变化并不明显，呈不显著的缓慢上升趋势；西部班玛以北及南部部分地区ET0呈减少趋势且变化不显著。

图4(c)为利用Theil-Sen趋势度方法获取的湿地ET0时空变化趋势，所得结果与MK方法图4(b)相似，均呈由东北向西南递减趋势。趋势度β在-0.35～3.14mm/a之间变化，区域内迭部站ET0上升速率最大，为3.14mm/a；松潘站ET0上升速率最小，为0.34mm/a。东北部迭部—玛曲—若尔盖地区ET0明显增高(1.0～3.0mm/a)，是全区变化最明显的区域；中部久治—阿坝—红原—松潘一带及碌曲周边ET0缓慢上升(0～1.0mm/a)；西部班玛以北及南部马尔康、黑水之间地区ET0呈缓慢下降趋势。

3.3 若尔盖湿地ET0的EOF分析

为进一步分析若尔盖湿地参考作物蒸散量变化特征，对19个气象站点1961-2016年的年ET0距平矩阵进行EOF分解，得到代表湿地ET0空间分布类型的相互正交特征向量。贡献率越大的特征向量，其模态对应的ET0分布形式越典型。每一模态的极大值中心即ET0变化的敏感中心。EOF分解结果表明，特征值最大的前3个特征向量累计方差贡献率达75.91%，已能反映湿地ET0变化的主要空间分布特征。

第一特征向量方差贡献率为48.08%，代表了湿地蒸散变化的最主要特征。由图5(a)，ET0第一特征向量均为正值，说明若尔盖湿地ET0变化保持高度区域一致性，大尺度气候系统的作用使整个区域ET0有相同的变化趋势。第一特征向量高值位于若尔盖湿地东北部，说明迭部—碌曲—玛曲—若尔盖为年尺度ET0浮动变化强烈区域，对气候变化反应敏感；低值位于湿地西部及南部边缘，说明这些区域ET0浮动变化相对较缓。整个区域第一特征向量大体呈自东北向西南递减的空间分布趋势，说明从整个湿地尺度来看，ET0的浮动变化强度自东北向西南递减，这与图4中MK检验下和Theil-Sen趋势度检验下ET0变化趋势结果相一致。第二特征向量方差贡献率为18.02%，也是若尔盖湿地ET0空间分布的一个较重要形式。由图5(b)，第二特征向量在-0.51-0.33间变化，空间系数数值北负南正的纬向分布特征表明湿地南、北地区ET0具有反相位变化的空间特征，即北部ET0偏高(偏干)时，南部ET0偏低(偏湿)，反之亦然。第三特征向量特征贡献率为9.81%，也能对若尔盖湿地ET0空间分布做出一定解释。由图5(c)，第三特征向量在-0.44～0.25间变化，呈东西两端高而中部低的分布形式，表明东西两端与中部具有ET0反相位变化的经向分布空间特征。

图4 若尔盖湿地ET0空间分布特征Fig.4 Spatial distribution characteristics of ET0 in the Zoige Wetland

图5 若尔盖湿地ET0前三特征向量空间分布Fig.5 Spatial distribution of the first three feature vectors of ET0 in the Zoige Wetland

EOF分析得到的特征向量反映了若尔盖湿地ET0的时空演变主要特征。从年尺度ET0的空间分布特征来看，其在整体上受大尺度气候系统的影响，表现为全区蒸散变动的一致性，南、北反向差异和东西、中部反向差异分别为第二和第三空间结构特征。

3.4 若尔盖湿地ET0变化成因分析

根据通径分析方法对ET0和气象因子(降水量、最高气温、最低气温、平均气温、相对湿度、风速、日照时数和净辐射)进行逐步回归，找出彼此间相关性强而引起多重共线性的自变量。相对湿度、净辐射、风速和平均气温对ET0的影响显著(p<0.05)，其他气象因子影响不显著。将不显著因子排除，建立最优的回归方程。对ET0与影响显著因子进行通径分析，计算各因子对ET0的直接作用(通径系数)、间接作用、相关系数，最后分析各气象因子对回归方程估测可靠程度E的总贡献，结果见表1。其中分析确定ET0与气象因子关系时，计算ET0和气象因子均采用以年作基准时间尺度的1960-2015年气象数据，且仍采用19个站点数据均值作为区域的整体代表值，通径分析的具体计算方法参见文献[20]。

由表1可以看出，相对湿度对ET0的通径系数为-0.57，对回归方程估测可靠程度E的总贡献为0.40，在各项指标中最大，说明相对湿度是影响若尔盖湿地ET0的最重要气象因子；净辐射、平均气温和风速对ET0的通径系数分别为0.55、0.53、0.35，对E的贡献分别为0.38、0.24和0.13，对ET0的影响依次减弱。其中，相对湿度对ET0的变化起负向直接作用(通径系数为负)，净辐射、平均气温和风速则起正向直接作用。这是因为辐射越强，则地温升高，促进水热交换。而气温上升，将引起下垫面水分运动加快，并促进植被蒸腾，说明在全球变暖的大背景下，气温的升高已经开始对区域蒸散产生影响。风速增大则利于大气中水汽扩散与热能传递，促进蒸散。结合各气象因子变化的趋势度及置信度和通径分析可知，相对湿度以-0.40%/10a趋势显著降低(p<0.01)，而净辐射、平均气温和风速分别以0.01MJ/(m-2·10a)趋势(p<0.1)、0.36 ℃/10a趋势(p<0.01)和0.04m/s趋势(p<0.01)显著上升。故近年来若尔盖湿地ET0的上升主要是因为相对湿度的显著下降，其次是净辐射、平均气温和风速的显著上升。若尔盖湿地位于青藏高原黄河源区，本文研究结论与王朕等相对湿度下降和平均气温上升是当地干旱化主要原因的结论相一致[22]，亦与王建兵等气温上升、相对湿度下降和降水减少是导致若尔盖湿地ET0上升的主要因子的结论相似[1]。然而国内大部分地区蒸散主要受风速和辐射、日照时数影响[23]，这与本文相对湿度影响最为显著、其次才为净辐射的结论存在差异。这主要是由当地独特的地形地貌及气候条件导致的。据尹云鹤等的结论，我国ET0对相对湿度的敏感性最高, 但由于相对湿度变化趋势往往不明显, 因此在我国大部分地区其并不是ET0的主导气候因子[23]。然而若尔盖湿地近年来相对湿度变化显著，其Z值达到-3.49，故其导致的ET0的变化更加明显。另外，近年来湿地低云量增幅显著，由此导致日照时数减少，且低云对太阳总辐射的削弱作用很强，这也在一定程度上减弱了净辐射的影响程度。

由各因子的间接作用可知，净辐射间接和为0.66，相关系数为0.72，且主要是通过相对湿度、平均气温对ET0产生正向影响，作用系数分别为0.42、0.19。这表明，净辐射与ET0的间接和及相关系数较大，主要是通过相对湿度和平均气温实现的，且净辐射、平均气温与风速的间接作用分析中，相对湿度的作用系数明显高于其他因子，依次为0.42、0.18和0.15，进一步反映出相对湿度为若尔盖湿地ET0的主要气象影响因子。

可见，各气象因子之间相互制约、相互影响。影响若尔盖湿地ET0的主要气象因子为相对湿度，且能够综合其他因子对ET0产生作用。

表1 气象因子对参考作物蒸散量的通径分析Tab.1 Path analysis of reference crop evapotranspiration by meteorological factors

4 结 论

(1)1960-2015年若尔盖湿地年ET0均值为625.3mm，并以4.89mm/10a的速率显著上升(p<0.01)，四季ET0表现为夏季>春季>秋季>冬季。年、秋、冬ET0分别在1968年(p<0.01)、1997年(p<0.01)、2003年(p<0.1)发生突变上升，春、夏两季未出现突变。

(2)湿地年均ET0呈南部、东部边缘高、西北—东南一线较低的空间分布特征，且变化速率呈由东北向西南递减趋势，其中西部班玛以北及南部马尔康、黑水之间地区ET0呈缓慢下降趋势。

(3)ET0第一特征向量均为正值，说明若尔盖湿地ET0变化保持高度区域一致性，南、北反向差异和东西、中部反向差异分别为第二和第三空间结构特征。

(4)影响若尔盖湿地ET0的主要气象因子为相对湿度，且能够综合其它因子对ET0产生作用。近年来若尔盖湿地ET0的上升主要是因为相对湿度的显著下降，其次是净辐射、平均气温和风速的显著上升。