研究区位于内蒙古自治区通辽市科尔沁左翼后旗阿古拉镇,地处科尔沁沙地南缘,地理坐标为122°33'00″—122°41'00″E,43°18'48″—43°42'24″N,面积约55 km²。总体上研究区地形自西向东、由南向北缓慢倾斜,西高东低。境内海拔最高为232 m,最低为186 m。研究区为固定沙丘、半固定沙丘、流动沙丘相结合,坨甸相间分布,中间地带的草甸,农田和牧场镶嵌交错分布,沙丘、草甸、湖泊、农田和村庄分别占研究区总面积的54.5%,26.6%,5.2%,10.4%,3.3%。属于典型的半干旱荒漠化地区,温带半干旱大陆性季风气候,Ф20 cm蒸发皿多年平均水面蒸发量为1 400.3 mm,多年平均降水量为379.8 mm。研究区植被种类丰富,主要有差巴嘎蒿、沙蓬、少花蒺藜草、草麻黄、虫实、冷蒿、芦苇、羊草、蒲公英等。土壤类型有:砂土、壤砂土、壤土、黑钙土、栗钙土等。土壤中砂粒占80%以上,有的甚至达到100%。

沙丘、草甸试验站点分别设有高10 m的微气象观测塔,在塔的不同高度分别布设了风、温、湿梯度传感器,并布设了辐射及涡度相关系统。两个站点所安装的均为开路式涡度通量观测系统,附加了温湿度传感器用于修正涡动相关系统所测定的空气温湿度。站点的涡动相关系统架高分别为2.62 m,4.61 m,采集频率均为10 Hz。具体观测项目见表1。

表1 涡度相关塔测量项目及仪器

本文利用的遥感数据为Landsat-8/OLI影像(2016年5—10月),获取自USGS(https:∥earthexplorer.usgs.gov/)。取清晰度高、无云或云量较少的6期影像(5—9月),每月一幅,之后对影像进行几何校正、辐射定标、大气校正等预处理工作。

气象数据(气温、标准等压面的气压、露点温度)源自中国气象数据网(http:∥data.cma.cn/),使用通辽国家气象站实测数据(2016年5—10月)。其他气象数据(气温、风速、降水、土壤温度等)分别源自研究区试验站2016年(5—10月)沙丘与草甸两个观测站的实测数据。

涡度相关系统安装的理想条件是下垫面水平均一,气流稳定等,而在实际测量中,很难达到理想状态,因此涡度数据在使用前需要预处理,处理步骤主要有:异常值及野点剔除、坐标旋转修正、数据插补和求日平均等。

DEM数据为SRTM(ShuttleRadarTopographyMission)数据,经过简单裁剪拼接等处理工作后即可使用,其空间分辨率为30 m,由于空气温度随高程变化,因此需要DEM数据进行修正,数据来源于地理空间数据云(http:∥www.gscloud.cn/)。

Wang和Bras在估算蒸散发的过程中,引入最大熵增原理与地表能量余项法,这一方法使得估算蒸散发更加简便^[19-28]。本次研究使用他们的研究成果,并给出基于最大熵增原理的蒸散发模型的推理过程和公式。

对于无植被覆盖的地表(裸地),熵增函数D(H,E,G)的计算公式如下:

D(H,E,G)≡(2H²)/(I_s)+(2E²)/(I_a)+(2G²)/(I_e)(1)

式中:H是显热通量(W/m²); I_s为H的热惯性参数(W·m²·K·s^1/2); E是潜热通量(W/m²); I_a为E的热惯性参数(W·m²·K·s^1/2); G为土壤热通量(W/m²); I_e为G的热惯性参数(W·m²·K·s^1/2)。H,E和G是利用最大熵增原理,在R_n确定的前提下,运用拉格朗日乘数法极值化D(H,E,G),得到的计算公式:

式中:R_n是净辐射(W/m²),计算方法见公式(14); q_s为表面比湿(g/kg); T_s为地表温度(℃); σ为无量纲参数,可通过q_s与T_s得出; I₀为与显热通量独立的简洁表达式,在计算H的过程中需要利用二分法或牛顿迭代法,其计算精度需要提前确定并作为约束条件。

在MEP模型中认为,植被覆盖的地表为单一封闭冠层,这时G相对于H和E可忽略不计,即I_s=0,冠层上的能量平衡可以表示为:R_nE_v+H,E_v和H的计算公式如下:

式中:E_v为植被蒸腾作用的潜热通量(W/m²),此时的T_s为叶面温度; q_s为叶面的表面比湿。

MEP模型的优点在于输入参数少,只需要地表温度(T_s)、净辐射(R_n)以及表面比湿(q_s),地表温度的计算使用Ren等^[29-30]应用通用劈窗算法针对Landsat8数据开发的模型,公式如下:

式中:T_i是第10波段(Band 10:10.6～11.19 μm)的亮度温度; T_j是第11波段(Band 11:11.5～12.51 μm)的亮度温度; ε为第10波段与第11波段的比辐射率平均值; Δε为第10波段与第11波段比辐射率的差值; b_k(k=0,1,…,7)为劈窗算法的系数; b_k的取值随着大气水汽含量的变化而有所差异。使用劈窗协方差表示方差比法(Split Window Covariance-Variance Ratio,SWCVR)可以计算大气水汽含量,计算公式如下:

W=a+bτ_j/τ_i(9)

式中:W为大气水汽含量(g/cm²); a和b为系数从模拟数据中获得; τ_i和τ_j分别为i和j波段的大气透过率。

为了方便后续模型的计算以及数据验证,将表面比湿换算为相对湿度(RH),本文将Peng等^[31]使用MODIS数据估算大气水汽含量方法移植到Landsat8数据,公式如下:

RH=e/e_s(10)

e_s=611exp((17.27T_a)/(237.3+T_a))(11)

e=q_s×(p_a)/(0.622)(12)

式中:e为地面水汽压(hpa); e_s为饱和水汽压(hpa); T_a为大气温度(℃); P_a为大气压(hpa); W与q_s之间有着良好的线性关系,使用通辽站的实测探空数据,得到W与q_s的回归方程,如式(13):

q_s=0.001(-0.682W²+6.677W-1.1123)(13)

SEBAL模型^[11-13]在世界上得到了广泛应用,在不同地区、气候类型以及不同土地覆被类型下,有大量学者对SEBAL模型进行了验证,相比其他遥感模型其精度相对较高,因此本文使用SEBAL模型进行对比验证,计算过程如下:

R_n的计算公式如式(14):

R_n=(1-α)R_S↓+R_L↓-R_L↑-(1-ε₀)R_L↓(14)

式中:R_s↓为入射短波辐射(W/m²); α为地表反照率; R_L↓为入射长波辐射(W/m²); R_L↑为出射长波辐射(W/m²); ε_o为宽波段表面热发射。

土壤热通量(G)计算公式如下:

G=(R_n(T_s-273.15)(0.0038α+0.007α²)(1-0.98NDVI⁴))/α(15)

显热通量(H)计算公式如下:

H=ρ_airC_p(T_s-T_a)/(r_ah)=ρ_airC_p(dT)/(r_ah)(16)

式中:ρ_air为空气密度(kg/m³); C_p为空气定压比热(取1 004.07 J/kg/K); T_s为地表温度或冠层的表面温度(K); T_a是2 m处的气温(K); dT为Z₁与Z₂两个高度之间的温差(K); r_ah是Z₁与Z₂之间空气动力学阻抗; Z₁,Z₂分别是指0.01 m和2 m高度处。

遥感模型中计算出的ET是卫星过境时刻的瞬时ET,在实际应用中,往往需要的是日ET,在本次研究中通过式(17)进行日尺度扩展^[29]。

ET₂₄=(86400ET_rF(R_n24-G₂₄))/λ(17)

式中:ET₂₄为日蒸散量(mm); λ为水分汽化潜热指数; ET_rF为卫星过境瞬时参考蒸散比(mm); R_{n_24}为卫星过境日的净辐射通量(W/m²); G₂₄为全天的土壤热通量(W/m²)。

利用2016年5—10月沙丘和草甸两个观测站的涡度相关数据以及气象观测数据与MEP模型估算结果、SEBAL模型估算结果、地表参数进行对比。为了更好地评价两种模型,本文选用统计学中均方根误差(RMSE)、相对误差(MRE)与拟合系数(R²)作为评价指标,对比结果见表2与表3,图1。

表2 地表参数及MEP模型结果估算误差

表3 SEBAL模型估算误差对比

由表2与表3得到,MEP模型反演G值的均方根误差与相对误差最大分别为61.50 W/m²,34%,G的拟合系数最小为0.57; SEBAL模型反演E值的均方根误差最大为47.93 W/m²,G的相对误差最大为15%,拟合系数最低为0.32; 地表参数验证结果相对较好,R_n的模拟结果误差较小,其拟合系数最高为0.98,T_s的平均相对误差最小为2%。由图1看出:SEBAL模型与MEP模型之间反演结果整体误差较小,其中MEP模型的H值偏高,G的误差相对较小,个别日期的E值有着较大的差异。

图1 MEP模型与SEBAL模型对比验证结果

利用沙丘和草甸观测数据,通过MEP模型分别反演了两点的日ET,使用实测ET数据、遥感估算的ET数据以及Penmam-Monteith模拟ET之间进行对比,对比结果见图2。

图2 沙丘地、草甸地降雨,模型估算与观测ET值

从图中看出MEP模型估算ET和实测ET变化趋势基本一致,个别日期有所差异,遥感估算ET与实测ET以遥感反演ET值与观测ET以及Penmam-Monteith模拟ET之间误差较小,结果表明MEP模型能在不同的土地覆被类型模拟出精度较高的ET值。研究区ET总体波动较为平缓,表现为5—6月呈上升趋势,7—8月相对较高,9—10月处于下降趋势。沙丘的ET明显低于草甸,随着降水量的变化,沙丘的ET波动较为剧烈。降水对ET有着较大的影响,当有降雨事件发生的日期ET会随之降低,出现一个波谷,降雨事件结束之后ET会有明显升高,出现一个波峰。总体来说:在不同土地覆被类型下,MEP模型可以估算出相对精确的ET以及较为合理的变化趋势。

本文通过MEP模型与SEBAL模型分别反演出研究区瞬时潜热通量,通过日ET尺度扩展方法的到日ET值,日ET时空分布见附图13,附图14。

研究区日ET最大可以达到7 mm以上,主要集中在湖泊,湖泊相对较为稳定,日ET值在4～7 mm; 农田与草甸由于植被与农作物的影响日ET变化较为明显,尤其在6月、7月、8月,个别区域会有最大值出现,日ET保持在3～7 mm; 沙丘日ET一直处于较低水平,保持在1～3 mm。

SEBAL模型与MEP模型得到的ET空间分布图在5—7月差异较小,8—10月差异较大,二者差异较大的区域在半固定沙丘,误差在0.5～1 mm; 农田与草甸差异较小,误差在0.2～0.6 mm。主要原因有:(1)MEP模型以能量闭合为边界条件,而SEBAL模型会出现能量不闭合的情况,(2)输入参数不同,MEP模型所需参数少,只需要3个地表参数,而SEBAL模型输入参数多,参数的不同引起结果差异,(3)对于植被处理方法的不同:对于植被覆盖的地表MEP模型认为土壤热通量为零,直接计算显热通量,而SEBAL模型通过选取冷热点的方式估算显热通量,5—7月研究区处于春季,植被较少,8—10月研究区处于夏季植被生长茂盛,因此导致8—10月差异较大。

为进一步分析MEP模型与SEBAL模型的差异以及ET时空变化,按照不同土地覆被类型对ET时空分布图的ET平均值进行统计,得到各类土地覆被类型ET平均值变化曲线,见图3。

图3 2016年生长季不同土地覆被类型平均ET统计对比

随着时间的变化不同土地覆被类型下的ET值波动情况有所差异,但是整体呈现先增长后降低的变化趋势。5月所有土地覆被类型ET值均保持在较低水平,ET在2.5～4 mm; 6月—8月大部分区域的ET值均有较为明显的升高,而沙丘地由于水分补给不充足反而降低; 在9—10月ET呈现明显的下降趋势,10月除湖泊仍保持在相对较高的水平,其他区域ET值较低,ET空间分布情况与5月相似。湖泊、草甸、农田整体呈较高水平,受月份影响较为明显; 沙丘波动较为明显,受月份影响相对较小,沙丘类型不同,ET分布有较大差异,其中固定沙丘和半固定沙丘相似,而流动沙丘ET处于较低水平,尤其在6—8月差异较为明显; 村庄受人类活动影响较大,ET值分布相对较为集中并处于较低水平。

影响研究区ET时空分布的主要因素有:土壤类型不同,草甸地主要为黑钙土、栗钙土,沙丘地以砂土、壤砂土、壤土为主要的土壤类型,不同土壤类型下地表反照率有所差异,对太阳辐射的吸收也有差别^[32],导致ET的不同; 植被类型不同,随着生长季的到来植被的蒸腾作用所占蒸散发总量的比重越来越大,不同植被类型蒸腾作用有着较大差异,导致ET的不同^[33]; 土壤含水量以及供水条件不同,农田与草甸地的土壤含水量较高且水分供应充足,而沙丘地土壤含水量低,地下水埋深较深,除天然降水外,没有其他水源供给,因此沙丘ET要低于农田与草甸^[34]。