蒸散发模型计算过程

一、蒸散发模型计算过程

（一）Penman-Monteith模型计算过程

Penman-Monteith模型是典型的模型计算方法，由于可以在较长时间尺度上获得其所需的气温、风速、湿度和太阳辐射4种参数而具有较好的适用性，被联合国粮食及农业组织（FAO）推荐为参考作物蒸散发量的计算方法（Richard，1998），因此其也被称为FAO Penman-Monteith模型，并在全世界范围内得到了大量的验证与应用（Fitzmaurice，2005），已取得良好的效果。

在FAO Penman-Monteith模型中，除了需要以上4种可以在常规气象站观测中获得的基本参数外（图7-1），还需要借助一系列的地表通量、日地轨道数据等过程参数（表7-1）在理论公式中进一步推导计算，最终得到参考作物蒸散发量（ET0），如图7-2所示。

尽管FAO Penman-Monteith公式具有坚实的物理理论基础，但是此方法也存在一定的缺点，即公式中的参数辐射和土壤热通量项是不能直接测量的，只能通过公式计算。当地表有植被覆盖或供水状况不充分时，蒸发、蒸腾量的估算就会出现偏差，此时，就需要对该公式进行修正。因此，考虑研究区当地农业发展和作物种植情况，本研究在上述ET0计算的基础上，借助地表作物在不同生长时期中需水量的生物学特性的作物系数Kc，进一步计算实际蒸散发量ET：

由于研究区支柱产业为棉花种植，因此，本研究根据新疆当地棉花的生长周期并考虑不同生育阶段，逐月选取相应的作物系数Kc值（王梅等，2016）进行反演（表7-2）。而在冬季非作物生长期内，精河县域内地表基本被积雪覆盖，地表实际蒸散发量趋近于零，且近似于潜在蒸散发量。

修正后的FAO Penman-Monteith模型既考虑了影响蒸散发量的大气物理特性，又结合了植被在不同生长周期里的生理特性，具有很好的物理依据，能比较合理、有效地反演蒸散发量的变化过程，也为非饱和下垫面蒸散发量的反演提供了一种有效且简单易操作的方案（图7-3）。

然而，由于FAO Penman-Monteith模型在理论基础上依然为单层模型，并没有将植物蒸腾和土壤蒸发分开计算，只是将地表土壤和植被看作一个整体进行计算，完全忽视了土壤蒸发与植被蒸腾差异性，并且也没有考虑下垫面的粗糙度对植被阻抗和空气动力学阻抗的影响，所以该模型只适用于完全覆盖低矮植被地表的条件下的蒸散发量估算，而在非均质地表面的应用中往往受到限制。

（二）遥感蒸散发模型模拟过程

20世纪90年代以来，随着遥感技术的不断发展，遥感获取地表通量信息的能力越来越强，出现了不少以SVAT模型为基础的遥感蒸散发模型。利用遥感方法计算地表蒸散发量的难点是，遥感反演的地面辐射温度与空气动力学温度不相等，使得感热通量（H）的计算十分复杂。遥感蒸散发模型主要解决H的精确反演，或者避开H的求解，并尽量避免使用遥感技术无法获取的部分气象参数。遥感参数化的模型主要有SEBAL模型、MODIS蒸散发比模型、SEBS模型和基于分层能量切割算法的双层蒸散发模型等。本研究主要基于SEBAL模型和SEBS模型来估算研究区蒸散发量，并采用蒸散发比不变法推算研究区的日蒸散发量，最终根据土地覆被类型，对比分析蒸散发量的时空格局变化（图7-4）。

1.SEBAL模型

地表能量平衡算法中SEBAL出现于1989年，是从对埃及西部沙漠的浅层地下水的蒸散发量的估计发展而来的。SEBAL最基本的特征是它不同于其他ET的估算——热红外算法，它是采用能量平衡估算水循环的过程，反演蒸散发量、生物量、水分耗散和土壤水分，能量平衡可以通过卫星数据量化，即从卫星影像中得到地表反照率、叶面积指数、植被指数和地表温度等地表特征参数。除了卫星图像外，SEBAL模型还需要一些气象数据，如风速、湿度、太阳辐射和空气温度。由于SEBAL模型是基于能量平衡，并不是水平衡，所以土地覆被数据、土壤类型、水文条件并不是必需的。无论是从全球范围到区域范围，还是从区域到农场农田，能量平衡适用于各级空间和时间分辨率的卫星图像。该模型的技术流程如图7-5所示。

SEBAL模型有几个不同于其他遥感通量算法的特点：

（1）感热通量取决于所谓的“热点”和“冷点”像元。“热点”指非常干燥的区域（潜热通量为零）。“冷点”指非常潮湿的区域（感热通量为零）。通过这两个极端点可以确定感热通量或蒸散发比。

（2）dT（空气温度的垂直差异）是反演计算感热通量的关键点。这意味着dT的计算不需要测量辐射表面温度和空气温度。

（3）dT和辐射表面温度呈线性关系。这种关系取决于选择的卫星影像（面积、气候、卫星过境时间），这种方法通常被称为“自校准”。

（4）表面热通量分数在一天中是恒定的。但是，对蒸散发比和相对蒸散发量而言，一个简单的修正都可以对其产生影响。

2.SEBS模型

利用遥感信息来估算陆地表面和大气之间热交换的方法可大致分为两种：先计算感热通量，然后求得潜热通量作为能量平衡方程的残差；或者先用一个组合方程计算一个指数（如作物水分胁迫指数），进而估算相对蒸散发量。感热通量的估算虽然在小尺度均质表面可以实现，但其在一些结构和热量非均质的冠层的应用仍存在困难。经典的遥感通量算法基于将地表温度测量和空间恒定的气象参数相结合，可用于小尺度地表通量的估算，但在较大的尺度上，地表气象参数不再是恒定的，地表结构和热量情况既不均匀，也不恒定。因此，需要设计更先进的，可应用于更大的尺度范围、具有非均质表面的地表的算法模型。在此过程中，由荷兰学者苏中波提出的蒸散发量计算方法——表面能量平衡系统（SEBS）模型，在反演各种尺度下的湍流热通量和蒸散发数据时具有较高的精度。

SEBS模型理论包括：一套用于地表物理参数反演的工具参数，如地表反照率、地表比辐射率、地表温度、植被覆盖等光谱反射和辐射参数；一个计算热传导粗糙度长度的扩展模型；一个基于能量平衡在极端状态下计算蒸散发比的新方法，该模型的技术流程如图7-6所示。

SEBS模型需要输入三组参数信息。第一组参数，包括由遥感反演得到的地表反照率、地表比辐射率、地表温度、植被盖度、叶面积指数以及植被高度（或地表粗糙度）。其中，当植被没有明确可用信息时，可用归一化植被指数替代。第二组参数，包括参考高度下的气压、温度、湿度和风速。参考高度为点应用的测量高度和行星边界层的高度。这组数据也可以由大尺度气象模型估算。第三组参数，包括向下的太阳辐射和长波辐射，可以是直接测量数据，也可以是模型输出或参数设置。

为了得到感热通量和潜热通量，SEBS模型虽采用与SEBAL模型相似的理论进行计算，但与SEBAL模型不同的是，SEBS模型区分了大气边界层和大气近地面层。也就是说，SEBS模型基于总体大气边界层相似理论和莫宁-奥布霍夫相似理论，对大气近地面层剖面动量交换和感热交换进行稳定度修正，将地表参数与大气边界层平均风速、平均温度及地表位温差等相结合，通过整体参数化求解，最终获得摩擦风速、奥布霍夫稳定度长度及感热通量。

3.遥感蒸散发模型主要参数反演结果

（1）归一化植被指数与地表植被覆盖度。

归一化植被指数（NDVI）是目前反映效果最好，且应用范围最广的一种表征地表植被覆盖情况的遥感指标。其取值范围在-1～1，-1～0表示地表覆盖为雪、云或水体等，大于0表示陆地表面有植被覆盖，其中也包含低度植被覆盖的裸地、戈壁等。

归一化植被指数的遥感反演结果（图7-7）显示，研究区西北部和东南部植被指数较高，而西南部和东北部较低。这主要是由于研究区西北部为精河下游绿洲，分布着广袤耕地；东南部为山地，多分布着山前草甸与天山云杉等；而东北部为艾比湖及湖周盐碱地，植物无法生长；西南部是戈壁、荒漠。通过多年变化对比发现，精河县的高植被区域不断扩大，区域平均植被指数已由1990年的不足0.01增长到2016年的0.12。其中，高度植被覆盖区域（0.7＜NDVI≤1）由0.588％增长到4.720％，低度植被覆盖的裸地（0＜NDVI≤0.2）区域由29.449％减少到25.843％。图7-8显示，第一个波峰在不断从趋近0的位置向右侧移动，第二个波峰即高度植被覆盖区域，也在逐步显现、增高。

地表植被覆盖度是植被在地面的垂直投影面积占总统计区面积的百分比，主要基于像元二分模型计算，计算方法如公式（7-2）、公式（7-3）所示。其中，NDVIsoil表示低度植被覆盖或完全是裸土区域的NDVI值，NDVIveg则表示完全被植被所覆盖区域（即纯植被）像元的NDVI值。为避免噪声影响，NDVI最大值和最小值一般根据研究区实际情况在一定置信度范围内取值。

式中，ρNIR为近红外波段波谱反射率，ρRED为红色可见光段波谱反射率，Fc为植被覆盖度。

研究区植被覆盖度在空间分布状况与归一化植被指数基本一致，多年整体变化状况比较明显（图7-9）。1990年研究区平均植被覆盖度为0.076，1994年为0.202，1998年为0.232，2002年为0.235，2007年为0.257，2011年为0.289，2013年为0.262，2016年为0.319。近20年，植被覆盖度以年均15.986％的幅度不断增长。其中，耕地面积的不断扩大是主要贡献因素。其次，山区林地由于多年保护恢复，植被冠层得到良好生长，植被覆盖度也显著提升。

根据图7-10分析，整体像元统计曲线呈非标准正态单侧分布，年平均值不断增大，曲线的中心位置向右移动，1990年标准差为0.156，2016年为0.318。其中，1990的统计曲线不平滑且呈波动性明显，这与早期精河流域绿洲开发导致土壤荒漠化、盐渍化有较大关系。2013年和2016年的曲线平滑稳定，基本无波动变化，说明在进入21世纪后，精河的破碎化地表状况得到了较高的抑制，保护恢复成效显著。

（2）地表比辐射率和地表反照率。

地表比辐射率是地表辐射能量与同温度下黑体（理想的辐射体）辐射能的比值，它的大小与研究区植被覆盖度有很大的关系，也能从侧面反映研究区对辐射能量的吸收状况，是精确计算地表温度的基础。往往地面比辐射率的计算是要区分地表覆被情况的，对典型地物分类进行计算，但由于上文中对植被覆盖度作像元二分法，相当于已进行分量计算，因此地表比辐射率计算公式如下：

式中，ε为地表比辐射率，Fc为植被覆盖度。

在所有到达地面的总辐射中，被地面反射回大气的辐射被称为地表反射辐射，用向上反射的辐射能量与入射辐射的总能量的比值表示地面反射能力的大小，称为地面反照率。本研究主要通过Liang等建立的Landsat TM数据反演，计算公式如下：

式中，Aldobe为地表反照率，ρTMn为第n波段波谱反射率。

由于Landsat-8与Landsat-7卫星的传感器波段内容和波段号略有区别，因此需对Landsat-8宽波段反照率反演的方程进行修改，计算公式如下：

根据图7-11分析，地表比辐射率的高值区和低值区与植被覆盖度空间分布完全一致，变化范围在0.986～0.990，年平均值的变化幅度为0.0001，1990—2016年间多年标准差分别为0.000484、0.000932、0.001036、0.00093、0.001183、0.001268、0.001292和0.001274。

具体而言，1990—2016年研究区地表比辐射率平均值分别为0.250、0.251、0.288、0.274、0.176、0.166、0.191和0.193，其中98.204％面积区域主要集中于0～0.4的范围内。高值大部分出现在农田灌溉渠下游和干涸湖盆外围的未利用地，其成因与农田灌溉、“洗盐”排水具有较大关系。低值大部分出现在艾比湖湖区及山区林地，主要是因为水体和林木树冠层可以吸收大量太阳辐射（图7-12）。

图7-13表明，研究区整体地表反照率逐年降低。在不同的土地利用类型间对比分析发现，1990—2016年未利用地的地表反照率平均下降了0.067，而其他土地利用类型变化程度都较小，这充分说明多年来，荒漠植被已经在恢复生长。

（3）地表温度。

地表温度（Ts）是太阳辐射到达地面后，被地面吸收产生的温度，是遥感反演蒸散发量过程中的一个重要参数。地表温度直接决定了大气长波辐射和地表长波辐射的分布状况，同时也是计算土壤热通量、感热通量和潜热通量几个地表通量重要参数的依据。

本研究利用遥感传感器TM/ETM+/TIRS热红外波段，依据大气校正法计算地表温度：

式中，K1、K2为热红外波段信息参数，对于不同传感器，其参数值也不同，其中：传感器TM中，K1=607.76，K2=1260.56；传感器ETM+中，K1=666.09，K2=1282.71；传感器TIRS中，K1=774.89，K2=1321.08。

黑体热辐射亮度（BTS）利用辐射传输方程计算获取：

式中，BTS为黑体热辐射亮度，Lλ为热红外波段的辐射亮度，τ为大气中热红外波段的透过率。

上述计算中所需要的近地表相关参数，主要来自NASA公开的大气剖面反演模型网站（http：//atmcorr.gsfc.nasa.gov）（表7-3）。

地表温度的反演结果（图7-14）显示，研究区内地表温度的范围在250～330 K，各年份均值分别293.86 K、292.83 K、297.17 K、293.24 K、302.00 K、299.70 K、306.21 K、310.52 K，上升幅度达到16.66 K。最大值出现在艾比湖以东荒漠，其地表主要是沙石砾漠，无植被，最小值则出现在精河县南部平均海拔在3500米以上的高山。艾比湖作为咸水湖泊，1990年水体表面的平均温度为286.30 K，2016年为293.54 K。

地表温度逐年分量统计结果（图7-15）显示，整体地表温度分布均呈现向高值区转移的趋势，即像元累积频率在显著提升到某一温度值时会趋于平缓。这一温度值在1990—2016年间，分别为302.8 K、302.775 K、307.314 K、305.562 K、314.515 K、313.804 K、317.601 K、320.743 K，平均上升幅度为19.943 K。

不同的地表覆被对地表温度的影响是不同的（表7-4）。我们结合土地利用和覆被分类图和气象站测量的当日温度，分析不同的下垫面类型的地表温度。最终发现，未利用地的地表温度最高，林地的最低，且各项数据都与实际状况相符，不存在异常，较好地反映了地表能量分布状况。

（4）净辐射通量。

地表净辐射通量称为辐射平衡或辐射差额，即可用的实际辐射能量。研究区净辐射通量的计算公式如下：

式中，Rn为地表净辐射通量（W/m2），RS为入射的短波辐射量（W/m2），α为地表反照率，R为入射的长波辐射量（W/m2），R为向外的长波辐射量（W/m2），ε为地表比辐射率。

入射长波辐射量（R）根据史蒂夫·波尔兹曼定律估算。

式中，R为入射的长波辐射量（W/m2），τsw为大气单向透过率，σ为史蒂夫·波尔兹曼常数（5.67×10-8W/（m2·K4）），Tα为空气温度（K）。其中，空间空气温度数据可用地表温度的中纬度修正函数获得。

出射长波辐射量（R）的计算公式如下：

式中，R为出射的长波辐射量（W/m2），ε是地表比辐射率，σ为史蒂夫·波尔兹曼常数，Ts为地表温度（K）。

其余涉及的相关辐射参数计算公式如下：

式中，RS为入射的短波辐射量（W/m2），Gsc为太阳常数（1367 W/m2），dr为日地距离修正系数，τsw为大气单向透过率，θ为太阳天顶角（rad）。

精河县净辐射通量的反演结果（图7-16）表明，高值区主要出现在山区林地，平均值达743 W/m2，且研究区内的水体的净辐射通量较大，为630 W/m2左右，低值通常出现在农田灌区下游和精河县中部戈壁地带，平均值为309.433 W/m2。通过多年净辐射通量的比较，1990—2016年如图所示年份的平均值分别是339.464 W/m2、348.543 W/m2、351.217 W/m2、369.225 W/m2、470.708 W/m2、519.519 W/m2、471.558 W/m2、480.372 W/m2，其最小值出现在1990年，最大值出现在2011年。

各年份的标准差分别为86.249、83.805、92.338、78.110、101.551、100.006、96.642和99.898，表明各个年份的净辐射通量分布较为集中，且相互之间差异不大。像元统计同时表明（图7-17），虽然随着时间变化，与之相关的地表参数都发生了变化，但净辐射通量的波动仍然保持大体一致。

（5）土壤热通量。

土壤热通量是指由于传导作用而存储在土壤和植被中的那部分能量，是土壤向下传输的热量，也是地表与深层土壤温度的温差导致的土壤热量交换，它对蒸散发和地表能量交换都有影响，其计算公式如下：

公式（7-15）在植被覆盖区应用效果较好，但在非植被的非均质下垫面中，仍需要参照土壤热通量与地表净辐射通量比值表进行计算（表7-5）。

图7-18显示，精河县的土壤热通量的空间分布存在明显差别。与净辐射通量相似，水体和南部山区的高海拔区域，土壤热通量值较高，而绿洲农业区的土壤热通量值较低。土壤热通量像元统计分析表明（图7-19），1990年到2016年间，最大值分别为216.709 W/m2、217.751 W/m2、231.614 W/m2、305.995 W/m2、266.442 W/m2、279.694 W/m2、268.821 W/m2、274.225 W/m2，最小值分别为6.642 W/m2、-0.773 W/m2、-2.994 W/m2、-13.362 W/m2、5.216 W/m2、21.601 W/m2、0.151 W/m2、8.932 W/m2，平均值分别为98.989 W/m2、88.738 W/m2、85.715 W/m2、91.542 W/m2、112.774 W/m2、120.017 W/m2、112.533 W/m2、108.390 W/m2。

研究区土壤热通量与净辐射的比值（G/Rn）结果为：耕地0.050、林地0.041、草地0.052、未利用地0.278、水体0.315。与参考数值对比并做误差分析，相对误差均小于15％，说明所反演的土壤热通量符合实际情况，同时也符合蒸散发量估算条件。