基于Landsat8辐射模块耦合SEBAL模型蒸散发估算

翟劭燚，王文种，刘九夫，2，王 欢，2，陆之昂，2

(1.南京水利科学研究院，江苏南京210029；2.水利部应对气候变化研究中心，江苏南京210029)

0 引 言

蒸散发是地表和大气间水分、能量交换系统中的重要过程，是水文循环和能量循环的重要组成部分[1]；同时也是陆地表层水循环最大、最难估算的部分。传统的估算蒸散发方法是基于点源的计算，研究方法包括：蒸渗仪法[2]、波文比能量平衡法[3]、涡度相关法[4]、蒸散速率与其他气象要素的关系建立经验公式推算法和水量平衡法[5]。传统的气候水文学方法对地表蒸散量的计算，仅限于离散点的观测与小区域范围内估算。由于自然下垫面的异质性，基于点尺度观测数据在区域尺度上代表性较差，进行空间尺度扩展的准确性未必能满足需求。

随着遥感技术的不断发展，连续、动态地获取大尺度区域地表信息的能力不断提高，利用遥感手段对区域尺度与蒸散发相关的地表信息和遥感参数进行监测，利用遥感技术测量和估算区域尺度的蒸散发算法也成为新的技术手段，近年正蓬勃发展。遥感数据被引入各种估算模型，与气候数据相结合，能够实现大空间尺度的蒸散量计算，对水文水资源的综合管理以及气候变化、农业发展等都有重要意义。对于干旱、半干旱地区生态系统合理利用和分配水资源的研究，迫切需要深入了解不同植被覆盖和土地利用条件下的蒸散发情况[6]。在水文、气候、农业和生态学中研究，对地表蒸散量的准确估算将会产生巨大的经济、社会、生态和科学效益。

SEBAL模型是基于能量平衡原理通过遥感技术反演蒸散发的典型方法。模型假设研究区域存在冷热极限像元，很好地考虑了干旱地区空间异质；因此，SEBAL模型[7- 8](包括METRIC模型[9]和M-SEBAL模型[10])在干旱区遥感蒸散发估算中得到了广泛的应用[11]。然而，SEBAL模型计算中采用了许多经验公式，特别是净辐射和土壤热通量等输入参数的估算，具有区域适用性；这样就会由于SEBAL模型实际应用时的区域、时间、大气条件、影像质量等的不同使参数的普适性受到影响，最终会影响反演精度[12]。这些公式的地表参数在应用到其他区域时，需要地面观测数据核实验证。此外，为减少净辐射估算的不准确性，大多数以前的研究选用有限的遥感影像，这限制了时间尺度扩展。因此，为使SEBAL模型更具有普遍适用性，需要研究机理性方法。

龙笛[13]研究表明，SEBAL模型中显热通量计算对于净辐射(RN)不敏感，根据地表能量平衡方程的余项法计算得到的潜热通量，应该是敏感的。RN作为主计算潜热通量的主要驱动力决定蒸散量(ET)的空间分布[17]。RN受许多因素影响，包括时间、地点、土地覆盖的影响，大气条件等，通常是额外地面太阳辐射、大气辐射、大气透射率函数、地表反照率、表面温度和表面发射率等参数的函数[14-15]，SEBAL模型采用经验公式简化计算这些参数。例如，大气透过率，这是由如气溶胶因素约束的水汽、臭氧深度，潜热通量对于大气透过率敏感性较高，SEBAL模型中简化为只有高度的函数[16]，这当然会导致高估，如果不仔细考虑各种大气因素操作。Bisht等人[14,17]提出的机械方案可能是一个合适的选择，因为它能够捕捉到RN空间分布且无需地面数据作为模型的输入，也清楚地认识到在空间上变化的输入参数的要求。

Landsat8数据空间分辨率高、光谱信息丰富、定标精确，通过适当的方法就可以反演丰富的地面信息，是反演区域蒸散发的理想数据源，为区域蒸散量的监测提供支持。本文在总结文献的基础上，分析采用Landsat8数据反演区域蒸散发基本原理和过程并介绍SEBAL模型中所需参数的估算方法，耦合新的净辐射能量估计方案改进SEBAL模型，并在半干旱地区进行了对比验证。对利用SEBAL模型采用Landsat8遥感影像估算区域蒸散发的研究及应用有重要的指导意义。

1 模型基本原理

改进的SEBAL模型与SEBAL模型都是基于能量平衡方程，对于研究区域，网格内陆面能量平衡方程可以表达为

RN=G+H+λET

(1)

式中，RN是净辐射量，W·m-2；G是土壤热通量，W/m2，即地面与土壤间的热量交换；H为感热通量，W/m2；λ是水的汽化潜热，W/(m2·mm2)。

利用晴朗天气下的遥感数据反演相关地表参数(地表反照率、植被指数、地表发射率、地表温度等)，结合相应气象及地表观测资料(气温、风速、平均高度等)，获得区域范围的净辐射通量、土壤热量，并根据遥感图像中冷、热像元点的选取确定地表温度与温度梯度差的线性关系，通过Monin-Obukhov相似理论迭代计算求得显热通量，从而根据余项法得到区域潜热通量和瞬时蒸散发λET值，最后通过蒸发比不变原理[18]进行时间尺度扩展，求得区域日蒸散量或更长时间的蒸散总量。

图1 改进SEBAL模型流程

1.1 净辐射计算

RN=(1-α)RS↓+RL↓-RL↑-(1-εs)RL↓

(2)

式中，α是地面反照率，W/m2；RS↓是下行的太阳短波辐射；RL↓是下行的长波辐射；RL↑是上行的长波辐射；是地面比辐射率。

太阳短波辐射能，它与太阳高度角、太阳辐射强度、大气条件密切相关。Niemela[19]对比分析了6种不同晴空条件下的太阳短波辐射，认为Zillman[20]提供的晴空条件下的短波辐射计算方法简单，所需参数较少。公式为

(3)

式中，I0为太阳常数，取值1 367；θ为太阳天顶角，rad；1/R2为日地距离订正因子，无量纲；τsw为短波大气透射率，某特定地区上空大气透过率受气象因素影响很大，在较短时间内，大气的压强、湿度、气体密度可发生明显的变化，透过率会因此而发生较大程度的改变。大气单向透射率的值约为0.55～0.85，如果研究区的面积很小、海拔较低的时候可以用一个值代替，约取0.75，当研究区面积较大时，必须将τsw转化为研究区面上的计算，本文采用ASCE-EWRI推荐的公式[9]

(4)

式中，P为大气压强，kPa；K为无量纲校正系数，晴朗天气下K=1.0；浑浊天气下取K=0.5；w为大气水汽含量，mm。王猛猛等[21]采用SWCVR算法应用于Landsat8数据的水汽含量反演结果精度较高，得到Landsat8热红外大气透过率比值和水汽含量的关系。即

(5)

(6)

式中，τ10、211为第10、11波段的大气透过率；z为地面高程(m)。到达地面的长波辐射能量可由斯蒂芬-玻耳兹曼(Stefan-Boltzmann)方程计算。即

(7)

(8)

式中，σ为常数；εa为空气比辐射率，Ta为空气温度；Ts为地表温度。

在地表反演过程中，水汽是大气透过率估算所必需考虑的主要因素。通常的做法是通过MODT-RAN、6S等大气模型软件模拟大气透过率与水汽含量的关系。对于Landsat 8数据而言，大气水汽含量很难从影像反演得到。宋挺等[22]使用与Landsat8影像同期过境时刻最为接近的采用多通道反演的MOD05-L2产品[23]，采用3次卷积内插法对MOD05-L2水汽产品数据进行重采样处理成与Landsat8数据分辨率一致，并用MODIS MOD11-L2地表温度产品与基于Landsat8反演的地表温度结果进行对比验证及精度分析，表明Juan C.Jiménez-Muoz劈窗算法(SW2)对输入参数敏感性最低。

Ts=T10+1.378(T10-T11)+0.183(T10-T11)2-

0.268+(54.30-2.238ω)(1-ε)-

(129.20-16.40ω)ε

(9)

式中，T10、T11为第10、11波段的亮度温度，K；ε为平均地表比辐射率。

1.2 土壤热通量(G)

土壤热通量是指进入土壤和植被内部的热交换能量，表征表层土壤与下层土壤之间的热交量换，是地表能量平衡中能量再分配的重要因素，是所有能量平衡方程分析中要涉及的关键因素。土壤热通量与净辐射热通量、显热能量相比，土壤热通量的值较小，其变化趋势与净辐射通量大致相同，采用参数化公式计算。即

G=0.3Rn(1-0.98NDVI4)

(10)

式中，NDVI为归一化植被指数。

1.3 感热通量(H)

感热通量(Sensible Heat Flux)是指以对流和传导形式进入空气中的那部分用于加热空气的热能量，是由于温度的梯度差异作用而造成的，其计算公式为

(11)

式中，ρ是空气密度，1.21 kg/m3；Cp表示空气定压比热，取值1004 J/(kg·K)；dT是零平面位移以上不同高度间的温差，K；rah为空气动力学的阻抗，s/m。rah是未知量，它们与不同梯度的温度、地表粗糙度和风速相关。

空气动力阻抗是摩擦风速、大气稳定度和地表动量粗糙度的函数。运用风速廓线理论计算

(12)

摩擦风速

(13)

式中，k为卡尔曼常数；z1定义为地表面植被覆盖处辐射能力转变为感热通量的高度，一般取值0.01 m；z2为不再受地表面粗糙影响的边界层高度，一般为2.00 m；u200为200 m高度处的风速，m/s；ψm为动量传输修正度；ψh为热量传输修正度；zzom为地表动量粗糙度，m。

由于表面加热导致低层大气内部的不稳定性，SEBAL根据Monin-Obukhov相似理论计算Monin-Obukhov相似长度及相应大气稳定度修正因子ψm和ψh，对rah进行校正，根据校正后的rah可重新求得dT，如此经过若干次迭代，直至H达到稳定，从而得到dT、Ts的经验回归系数a和b的最终值。根据区域内TS的空间分布求得dT的空间分布，进而求得每个像元的H，再根据余项法求解LE，通过时间尺度扩展便得到一天的蒸散量。

2 研究应用

Walnut Gulch流域(WG)位于美国东南部亚利桑那州的东南部((31°43′N， 110°41′W) ，是美国农业研究中心(USDA-ARS)和陆地水文研究(NASA)的长期研究区域，面积为150 km2。研究区属于圣佩德罗河上游流域(跨越索诺拉州、墨西哥州和亚利桑那州)，位于索诺兰沙漠和齐瓦瓦沙漠之间的过渡地带，属于半干旱地区，流域的海拔高度为 1 250～1 885 m之间(见图2)。

图2 Walnut Gulch流域示意

WG被认为是世界上观测数据最为密集的半干旱区域，提供丰富的数据资料可供下载(详见网址http:∥www.tucson.ars.ag.gov/dap/)。其包含丰富气象水文数据[26]、地理信息[27-28]、通量[29]、遥感数据[28]等数据资料。这些资料为基于遥感卫星的反演区域蒸发模型研究提供了先验知识和验证资料。该实验流域共设有三个自动观测气象站点，气象数据主要有：气压、气温、风速、相对湿度。区域内的两个地表通量观测站点记录的通量数据用于对比验证。

Landsat系列卫星为世界上应用范围最广应用次数最多的卫星之一，遥感卫星影像的空间分辨率和光谱分辨率等也在不断的提高。Landsat 8上携带有两个主要载荷：OLI(Operational Land Imager，陆地成像仪)和TIRS(Thermal Infrared Sensor，热红外传感器)。OLI包括了ETM+传感器所有的波段，对波段进行了重新调整，排除水汽吸收特征；OLI全色波段Band8波段范围较窄，这种方式可以在全色图像上更好区分植被和无植被特征；Landsat8上携带的TIRS载荷，将是有史以来最先进，性能最好的TIRS。TIRS将收集地球两个热区地带的热量流失，目标是了解所观测地带水分消耗，特别是美国西部干旱地区。数据可以通过http:∥glovis.usgs.gov/下载，本研究以2015337影像为例反演该区域地表蒸散发。

由计算结果提取草地和稀疏灌木站点的地表温度和通量数据，稀疏灌木和草地两个站点的landsat8估算值与地而观测值的对比见表1。从表1中可以看出，草地和稀疏灌木两个站点利用SEBAL模型估算的净辐射都存在高估的问题，与观测值相比偏差为175.7 W/m2和210.6 W/m2，而本文采用的计算结果偏差为7.4 W/m2和33.2 W/m2，计算的结果总体上要比SEBAL模型估算的结果更靠近1∶1线(见图3)。图3显示，净辐射计算值比SEBAL模型采用经验公式估算的结果更接近真实值。

表1 反演结果对比

图3 两种模型估算的净辐射能量反演结果对比

草地和稀疏灌木两个站点的观测潜热通量分别为24.0 W/m2和30.4 W/m2，残余法计算的潜热通量分别为109.2 W/m2和0.35 W/m2，SEBAL模型相对应于观测的计算结果偏差为30.5 W/m2和18.6 W/m2，相对应于残余法的计算结果偏差为-54.7 W/m2和48.6 W/m2；而本文采用的耦合模型计算结果对应于观测的偏差为15.6 W/m2和-26.3 W/m2，相对应于残余法的计算结果偏差为-69.6 W/m2和3.7 W/m2，模型算法的可靠性还需要进一步验证。2015337期反演区域地表蒸散发如图4所示。

图4 反演潜热通量结果(单位：W/m2)

3 结 语

SEBAL模型利用很少的输入参数就能反演区域蒸散量，反演的过程中， 模型中所需参数的估算是影响反演精度的主要因素 ，本文尝试利用Landsat8数据计算水汽含量及大气透射率代替经验公式来提高净辐射能量的估算精度及空间分布的准确性、可靠性，从能量平衡的角度阐述了计算净辐射耦合SEBAL模型的计算过程、步骤所需参数的估算方法，为模型在蒸散发反演中的应用提供帮助。