基于能量平衡算法的精河流域绿洲蒸散发时空变化模拟

于 辉 ，代鹏超，张金燕，毋兆鹏,2

(1.新疆师范大学地理科学与旅游学院，新疆 乌鲁木齐 830054；2.新疆干旱区湖泊环境与资源实验室，新疆 乌鲁木齐 830054)

水资源利用是干旱区资源环境和生态系统的核心问题，其中蒸散是干旱内陆水循环中水分消耗的最终途径，也是区域水量平衡和能量平衡的最活跃因子。精确估算地表蒸散发的时空变化，对于了解水文循环和能量平衡的整个过程，评价区域水循环和水平衡的功能，实现水资源、生态、社会和谐发展具有十分重要意义[1-2]。

1802年，Dalton提出关于蒸散发的计算公式，开启了国内外研究蒸散发的历史性时代。随后发展了一系列蒸散发观测和计算方法[3]。应用遥感技术是蒸散发估算领域新的一大亮点，为区域空间尺度的蒸散发扩展计算找到了解决方法，许多适用不同区域的和改进的蒸散发计算模型就是由此而生，如陆地表面能量平衡算法(Surface energy balance algorithm for land ,SEBAL)[4]、地表能量平衡指数(Surface energy balance index, SEBI)[5]、简化的地表能量平衡指数(S-SEBAL)[6]、地表能量平衡系统(Surface energy balance system, SEBS)[7]、MODIS蒸发比模型[8]和双层蒸散模型TTME[9]等。

研究区属天山北坡地区，位于全国“两横三纵”城市化战略格局中陆桥通道的西端，是全国主体功能区规划[10]确定的国家层面重点开发区域。研究区是北疆重要的优质棉基地之一，但随着近几十年区域人口增加及城市化进程加剧，加之干旱气候条件、水资源分布空间差异的综合作用，在绿洲开发过程中出现了一系列生态危机。2015年3月，新疆被确定为“丝绸之路”经济带核心区，使研究区从对外开放的“末梢”变成了开放前沿,也更使得这里成为了新时期我国开展生态环境国际合作和生态文明建设的重点区域。目前针对该区域蒸散发的研究虽已有部分成果，但数量较少且难以反映在“山地-绿洲-荒漠”系统背景下蒸散发的空间异质性，对全区未来蒸散发变化趋势的定量研究还相对比较薄弱。因此，本文利用基于能量平衡开发的Surfance energy balance alogrithm for land (SEBAL)模型，利用Landsat_5 TM遥感数据和地面气象数据，对精河流域绿洲地表通量及日蒸散发进行了估算，同时分析了蒸散发的时空变化规律和未来的变化趋势，以期为该区域水资源管理规划及其可持续发展提供科学依据。

1 研究区概况

研究区精河流域绿洲，位于新疆维吾尔自治区西北部准噶尔盆地西南边缘，天山支脉婆罗科努山北麓，介于82°32′50″E-83°17′49″E、44°36′30″N-45°08′23″N之间，包括整个艾比湖湖区及精河流域平原绿洲，总面积3 503.41 km2(图1)。该地区地处亚欧大陆腹地，距海洋较远，地势南高北低，属于北温带干旱荒漠型大陆性气候，日照充足，冬冷夏热，干燥少雨，蒸发量大，冬夏长，春秋时间较短。年平均太阳总辐射量冬季和夏季相差较大，最大可达到36.7 kcal·cm-2。降水量由南部山区向中、北部平原逐渐减少，年均降水102 mm左右，其中山区降水较为丰富，最大时可达约700 mm。研究区全年的蒸发主要集中在4～9月，多年平均蒸发量约为1 625 mm，且蒸发量大于降水量。

图1 研究区位置图Fig.1 Location of the study area

2 数据来源与研究方法

2.1 数据来源

研究所使用的1998、2007年和2011年遥感影像数据为Landsat-5 TM影像，拍摄日期分别为9月25日、9月18日、9月13日，三期数据含云量均低于5%，影像清晰度良好，经过几何精校正、辐射定标、Flaash大气校正后适合本次建模的基础条件。

土地利用数据参考中国科学院土地资源遥感调查与监测技术规程[11]，结合研究区的实际情况，将研究区土地利用类型划分为耕地、林地、草地、水域、建设用地、未利用土地6个类型，基于TM影像结合地形图、专题图等相关资料自行解译完成，并通过调查验证，误差控制在1个像元内，最终在ARC GIS平台下生成研究区土地利用数据库。

实际蒸散量主要根据以下5个气象站点：阿拉山口气象站点(站点编号：51232)、托里气象站点(站点编号：51241)、温泉气象站点(站点编号：51330)、精河气象站点(站点编号：51334)以及乌苏气象站点(站点编号：51346)的日值数据，利用联合国粮农组织推荐的Penman-Monteith公式，并借助表征作物在不同生长时期中需水量的生物学特性的作物系数Kc计算获得。

2.2 基于能量平衡算法的SEBAL模型

Bastiaanssen博士在1998年提出SEBAL模型，已经在多个国家和地区得到了广泛地应用。该模型依据陆地表面能量平衡原理，即地球表面所获得的净辐射量等于土壤热通量、感热通量和潜热通量之和。SEBAL模型的计算公式为：

Rn=H+LE+G

(1)

其中,Rn为地表净辐射通量(W·m-2)，H为感热通量(W·m-2),LE为潜热通量(W·m-2),G为土壤热通量(W·m-2)。

2.2.1 地面净辐射(Rn)的计算 地面净辐射是地表的主要能量来源，是地表通过太阳辐射的短波辐射、长波辐射过程得到的净能量(图2a)，计算公式如下：

Rn=(1-α)Rs+(Lin-Lout)-(1-ε)Lin

(2)

其中：Rn为地表净辐射(W·m-2)；Rs为太阳总辐射(W·m-2)；α为地面反射率；ε为地面比辐射率；Lin为大气长波辐射；Lout为地面长波辐射。

图2 研究区模型关键参数运算结果Fig.2 Diagram of operational results of key parameter of study area model

2.2.2 土壤热通量(G)的计算 土壤热通量指存储在土壤或植被中的能量，是一个相对较小的量，直接计算较困难，一般通过G与地表温度(Ts)、Rn、α、植被指数(NDVI)的统计关系求得，同时根据卫星过境时间进行适当的校正(图2b)，计算公式如下：

(1-0.978NDVI4)×Rn

(3)

其中：c11为卫星过境时间，在地方时12∶00之前取0.9，12∶00-14∶00取1.0，14∶00-16∶00取1.1。

2.2.3 感热通量(H)的计算 感热通量是大气稳定度、风速和表面粗糙度的函数(图2c)，计算公式如下：

(4)

式中,ρair是空气密度(kg·m-2)；Cp为空气定压比热(≈1004J·kg-1·K-1)；dT(K)为零平面以上高度z1和z2处的温差(T1-T2)；rah为热量传输的空气动力学阻抗(s·m-1)。

上式中，dT、rah均为未知量且彼此相关，因此借用Monin-Obukhov理论，通过循环递归算法进行求解(图3)。图中U为风速，U*为摩擦风速，Zom为动力粗糙度；L是Monin-Obukhov长度，反映地面层湍流特性的关键参量；Ψh、Ψm是感热通量稳定度修订函数，可以根据Paulson[12]提出的大气稳定度修正公式计算得到。1998年rah值在33.698～35.794之间震荡，取第14次值；2007年rah值在32.567～32.872之间震荡，取第12次值；2011年rah值在34.152～34.083之间震荡，取第13次值。

图3 Monin-Obukhov循环递归流程图Fig.3 The flow chart of Monin-Obukhov circular recursion

2.2.4 潜热通量(LE)的计算 通过模型分别反演出了净辐射、土壤热通量和感热通量后，运用公式5即可反演出卫星过镜时刻的潜热通量。

LE=Rn-H-G

(5)

公式中通量单位均为W·m-2，其中Rn、G、H都是基于像元得到的栅格数据。

2.2.5 全天蒸散发估算 运用瞬时潜热通量计算出蒸发比，将瞬时值延伸至全天的蒸散发值。SEBAL模型中，假设24 h之内蒸发比EF是相对稳定的，即

(6)

卫星过境时刻的瞬时蒸散发量计算公式为：

ETinst=3 600×(Rn-H-G)/I

(7)

其中，I为汽化潜热(J·kg-1),3 600是秒到小时的转换系数，所以单位是mm·h-1。则24 h蒸散发ET24为：

(8)

式中,Rn24、G24为日平均净辐射通量和日平均土壤热通量。

2.3 蒸散发周期预测

小波是一种特殊的长度有限、平均值为0的波形，倾向于不规则和不对称，能分析出时间序列周期性变化的局部特性，同时可以更清楚地看出各周期随时间的变化情况[14]。本文采用Morlet小波研究ET序列的特征尺度和周期性，其小波函数为：

(9)

在此基础上，结合Mann-Kendall(M-K)突变检验进行分析，序列大于0表明参数呈上升趋势，反之呈下降趋势，当其超过0.05显著性水平检验临界值(1.96)时，则上升或下降趋势显著[15]。

3 结果与分析

3.1 SEBAL模型精度验证

研究区内缺乏实际地表通量观测数据，因此，参考朱明承等[16]、刘文娟[17]和姜红[18]的方法，先基于研究区及周边站点气象数据，采用联合国粮农组织(FAO)推荐的P-M法计算参考作物蒸散发，再结合作物系数得出主要土地利用类型的实际蒸散发，并将此结果与遥感估算值相对比。结果表明，SEBAL模型计算结果相对误差均在18%之内，能够满足对反演结果的精度要求(表1)。

表1 SEBAL模型精度验证

Table 1 Precision validation of the SEBAL model

土地利用类型LandusetapeKC系数KCcoefficientFAOPM蒸散发均值/mmFAOPMevapotranspirationaverage遥感反演蒸散发均值/mmEvapotranspirationaveragebyremotesensing绝对误差Absoluteerror相对误差/%Relativeerror耕地Cultivatedland1．155．926．150．233．89林地Woodland1．085．566．430．8715．65草地Grassland0．854．385．150．7717．57

3.2 日蒸散量时空变化分析

3.2.1 日蒸散量空间分异特征 就空间尺度而言，研究区位于中国西北部典型干旱区，海拔落差较小，土地利用类型结构比较单一。在气候、降水、温度和土壤含水量等因素的综合影响下，该区域的日蒸散发模型反演结果有显著的空间差异性。图4反映了研究区不同年份夏季的日蒸散发分布情况，结果表明中部艾比湖湖面为极高值范围，这是由于湖面上方水分相对于其他地类比较充足，湖水比热容大，湖区温度上升和下降的幅度小。西南区除托里乡旦达盖沙漠外为蒸散发次高值区域，因为这部分区域地形较为平坦、海拔低，研究区耕地主要分布于此，基础灌溉设施完善，棉花此时也正处于花铃期末、蒸腾作用强烈，蒸散发值高。研究区的西北部是著名的风区阿拉山口，导致植被遭到破坏或覆盖度低，土地沙化现象较为严重，植物的蒸发、蒸腾量较少。艾比湖的北部、东部梭梭、红柳等灌木由于数年前的连年樵采而受到不同程度的破坏，导致地表覆盖较差出现了蒸散量低值区。研究区东南部开始进入天山山区，地形较高造成小范围水分的重新分配，林地面积占比高，故也造成该区域日蒸散发量较高。

从时间上看，1998-2007年，日蒸散发艾比湖湖区的高值范围变化不大，但西南区的次高值区域明显减少，取而代之的是低值区域的增加。2007-2011年湖区的日蒸散高值范围明显减少，同时西南区的次高值区域则明显扩大，低值区域明显减少。这一方面是由于该区域的耕地面积呈增加趋势，且增加趋势在加强；另一方面则是由于艾比湖湿地自然保护区在2007年晋升为国家级自然保护区，周围生态环境得到保护而导致这一结果。

图4 研究区日蒸散发空间分布Fig.4 Spatial distribution diagram of daily evapotranspiration in the study area

3.2.2 日蒸散量时间分异特征 利用精河气象站近40年数据和P-M公式，结合Mann-Kendall 检验对1970-2015年研究区年实际蒸散量序列进行分析(图5)，结果表明，1970-2015年正序列曲线(UF)总体为持续下降趋势，UF和反序列曲线(UB)于1973年在临界线区间内有1个明显的交点，表明研究区蒸散发在1973年发生突变，突变后精河流域平均年蒸散量较突变前减少218.318 mm，降低幅度为22.019%。进一步根据1970-2015年精河流域实际蒸散量数据，绘制Morlet小波变换系数的实部等值线图(图6)。结果表明，研究区蒸散量主要存在26～30 a的周期变化。

为探究研究区未来的蒸散量变化情况，本文基于实际蒸散量小波方差(图7)，对28 a处的小波系数值[f(t)]与年份t建立回归方程：f(t)=357sin(0.337 6t-115.7)，R2=0.98，P<0.05，并选择2018-2030年为预测期对研究区蒸散量变化趋势进行预测(图8)。结果表明，在现有条件不发生大的改变前提下，预测其整体蒸散量呈波动变化。其中，自当前到2022年蒸散量不断增加，并将于2022年发生突变，进入下降期。预测到2030年，蒸散量将再次进入上升周期。

图5 研究区蒸散量Mann-Kendall突变检验Fig.5 Mann-Kendall tests of the evapotranspiration rate in the study area for sudden changes

图6 研究区蒸散量Morlet小波系数实部等值线图Fig.6 Contour map isogram of real part of Morlet wavelet coefficients of ET in the study area

图7 研究区实际蒸散发量小波方差Fig.7 Wavelet variances of actual ET in the study area

图8 研究区实际蒸散发量小波系数(28 a)Fig.8 Wavelet coefficients of actual ET in the study area (28 a)

3.3 不同土地利用类型日蒸散发变化

不同土地利用类型使地表土壤的湿度和地表温度状况发生改变，影响陆面能量平衡，因此对区域日蒸散发的影响非常显著。利用ARC GIS空间分析技术对各土地利用类型平均日蒸散发量的研究表明(图9)，三期的土地利用类型的日平均蒸散发量均为：水域﹥林地﹥耕地﹥草地﹥建设用地﹥未利用地。

研究区的水域主要是艾比湖和精河绿洲的河流，由于水体位于北部和东北部地表比较平坦的区域，接受的太阳能辐射充足，因此实际蒸散发大于等于该区域的潜在蒸散发，也成为整个研究区日蒸散发最大区域，达到7.3～9.32 mm之间。研究区的林地主要分布在耕地边缘和河流附近等海拔相对较低但水分较为充足的地带，由于对地下水源的涵养能力比人工灌溉的耕地强，因此其日蒸散发量相比其他土地利用类型较高，均值达到6.43 mm。研究区耕地是主要土地利用类型且区域比较集中，由于主要作物棉花在9月初即将步入采摘期，棉农会停止浇水、施肥等田间管理活动，因此实际蒸散发情况相对于生长旺盛期而言较低，均值为6.15 mm。研究区草地自1998年到2011年呈现缩减的趋势，多数转变为耕地(表2)，加之草地的水源主要来自降水，因此由于本区域降水的缺乏，限制了该土地利用类型的蒸散发活动过程，使得该地区的实际蒸散发量不是很大，均值为5.37 mm。研究区建设用地在精河绿洲的占地极少，呈零星分布，主要由交通建设用地和建设用地等组成，不透水面积大，不利于水分扩散，导致潜热通量变小，使得整体蒸散发较低，日蒸散发值主要在4.49～5.34 mm之间，在研究区土地利用类型的蒸散发量图上表现出较低的区域。未利用土地在研究区为裸露岩石、戈壁、荒漠和沙地，是整个研究区分布最多的地类，由于地表覆盖度极低致使该类型的日蒸散发量成为了整个研究区的极低值区域。

图9 不同土地利用类型日平均蒸散量Fig.9 The daily mean ET of different land types in the study area

表2 1998-2011年研究区土地利用类型转移矩阵/km2Table 2 Land use transit matrix of study area from 1998 to 2011

4 结论与讨论

本文在蒸散发遥感反演的基础上，综合分析了研究区日蒸散发特征，得出以下结论：

1)研究区实际蒸散量近20 a间总体呈下降态势，虽然1995-2004年出现小幅增长趋势，但研究区蒸散量存在26～30 a的周期变化规律，故可预见2022年将再次转入下降，而2030年蒸散量则再次进入上升周期。

2)研究区中部艾比湖湖面为极高值范围，虽然东南部开始进入天山山区日蒸散发量较高，但西南区除托里乡旦达盖沙漠外的绿洲耕作区，实为较高值蒸散发的主体分布区域，西北部、北部、东部则为蒸散发低值区。在不同土地利用类型中，未利用地日蒸散发最小，其次是建设用地，除水体外林地和耕地蒸散发最高。

3)在自然界中，潜热(LE)和感热通量(H)都是伴随着地面净辐射通量(Rn)的变化而变化，同时LE和H的相互比例决定了能量的转换耗散方式。研究区三期潜热通量均值为450.6 W·m-2，感热通量均值为392.6 W·m-2，前者大于后者说明研究区内植物正处于生长旺季，较容易的水分来源与多样化的水分利用方式决定了其更高的LE过程，并主要以蒸散发的形式消耗能量。随着绿洲感热通量变小，潜热通量变大，湿度增强，对绿洲的保护范围加大，所以增加地面植被的覆盖程度，对绿洲应对灾害天气有一定的减轻作用。

4)研究区由于降水量较少，生态系统简单而脆弱，实际蒸散发波动一定会引起区域水热资源的多层次变化。一方面，实际蒸散发可减小辐射向感热的转化，增加空气湿度，提高最低气温及降低最高气温，起到调节气候的作用；但另一方面，则会增加流域农业需水量，促使农垦区增加灌溉，进而造成精河径流减少，加速湖泊面积萎缩，加剧流域的干旱和荒漠化程度，使原本脆弱的生态环境更加恶化。因此，有必要开展保护生态环境，加强生态管理，促进生态系统良性发展。