2.1 需水量与耗水量测定原理及数据标准化方法

2.1.1 水量平衡法

水量平衡法是一种间接的测定方法，通过计算代表范围内土壤水量的收入和支出之差来推求作物ET，即将作物耗水作为水量平衡方程的余项进行推求。其表达式为

式中：ET _c 为时段t内的作物需水量或耗水量；R为保存在土壤计划湿润层的有效雨量，一般取单次降雨量大于5mm的降雨为有效雨量；I为时段t内的灌水量；F为地表径流量，考虑到干旱区作物生育期很少存在地表径流，该项取为0；ΔW为时段t内土壤水平衡计算层内储水量的减少量，根据单次灌水量、土壤质地及制种玉米根系分布情况，土壤水平衡计算层取100cm；Q为时段t内耕作层下边界的地下水补给量或深层渗漏量，当地下水埋深较深时，地下水补给量可以忽略，且根据在干旱区甘肃河西走廊的土壤水分观测结果，灌水前后100cm处土壤含水量基本不变化，深层渗漏量可以忽略，故该项取值为0。式中各项均以mm计。

根据干旱区情况，式（2.1）可简化为

其中，ΔW可以用土壤含水量，由式（2.3）求得

式中：z _rt 为土壤水平衡计算深度，m；θ _z （t ₁ ）、θ _z （t ₂ ）分别为计算时段初和时段末的计划湿润层平均土壤含水量，m ³ /m ³ 。

该方法的优点是适用范围广，测量空间尺度可小至10m ² ，大至10km ² ，时间尺度从一周至一年 ^[1] 。非均匀下垫面和任何天气条件下都可以应用，不受微气象方法中许多条件的限制。因此，它最为适合下垫面不均一，土地利用状况复杂情况下的大面积耗水量测定，可用来对其他测定或估算方法进行检验或校核 ^[2] 。准确计算区域边界范围内外的水分交换量和取得足够精确的水量平衡各分量测定值，就可以得到较为准确的农田总耗水量。然而水量平衡法各分项，如地表径流项、深层渗漏项、地下水补给项准确值不易得到，降低了水量平衡法的精确性。在干旱气候环境下，地表径流项可忽略 ^[3，4] ，但是实际上地表径流项能否忽略依赖于降雨的特征（如雨量、持续时间和强度等），且仅对于特定的土壤种类（砂壤土）才能忽略 ^[1] 。很多研究表明，在干旱区可以忽略深层渗漏项，但是实际上此项是否能忽略依赖于土壤深度、渗透性及土壤储水性 ^[5] ，需要验证每个观测点土壤渗透性，确定土壤深层渗漏是否可以忽略。当地下水位较低时，地下水补给项可以忽略；当地下水位较高时，地下水补给项不可忽略。干旱区地表蒸发强烈，地下水补给项通常被忽略。

水量平衡法中，土壤含水量的测定影响较大。首先，由于土壤含水量具有空间变异性，导致取样点的土壤含水量情况不一定能够代表整个区域的土壤含水量特征。土壤含水量具有空间变异性有两层含义：其一局部灌溉模式条件下，土壤干湿分界明显，小尺度内不同测点土壤含水量差异较大，张宝忠在中国西北旱区应用水量平衡法时为了确保土壤含水量的代表性，分别测量沟里、垄上和遮阴三个位置的土壤含水量并取平均以增强其代表性 ^[6] ；其二由于土壤自身组成、结构不同引起土壤含水量具有空间变异性，为此应增加随机取样的数量，使所取得的样本在较高的置信水平下具有较高的精确性。Hupet等的研究表明，在6300m ² 的研究面积上，置信水平为 0.95，相对精确性为 10%条件下，需要17～31个样本点 ^[7] 。其次，土壤水分的测量方法影响测量值的准确性。常用的土壤水分测定方法主要有取土烘干法、中子仪测定法、时域反射仪测定法。取土烘干法是其他土壤水分测定方法的基础和依据，方法简单，精度高，但不能连续定位观测。中子仪测定法多用于科研，表层土壤含水量不易准确测定。时域反射仪测定法比较先进，能够连续、快速地对土壤水分进行定位观测，且无辐射，对土壤结构不会起破坏作用 ^[8] ，所测表层土壤含水率比中子仪测定法精度要高得多 ^[2] 。当水量平衡法各分项不能准确测定或是各分项对水量平衡所起作用较大而被忽略时，此方法产生的误差较大。最后，水量平衡法不能测定作物耗水的动态变化过程，也不能在短期内获得可靠的耗水量，而只能给出长时段（一般为一周以上）的总耗水量，若进行一周以下水量平衡的研究，需要增加取样点数量。

2.1.2 大型称重式蒸渗仪法

大型称重式蒸渗仪通过直接测定箱内土体重量的变化来测定植物蒸腾和土壤蒸发，测定精度较高，可自动记录小时，甚至分钟尺度数据，已经成为测定农田ET的标准仪器 ^[9，10] 。大型称重式蒸渗仪（简称蒸渗仪）测定作物ET的理论是水量平衡原理。对蒸渗仪内的土体使用水量平衡方程，其中，ΔW由称重系统直接称出，一般设置蒸渗仪称重时间间隔为0.5h或1.0h。称重系统在不同时间准确称取蒸渗仪内土体的重量，并由数据采集系统记录，相邻整点间重量的差值即该0.5h或1.0h内的ΔW，精度可达到0.02mm。由于作物种植在蒸渗仪内的土体中，地表径流量和地下水补给量为0。另外，深层渗漏量为土体排出的水量，可由排水系统准确获取。研究中蒸渗仪内土体深度可达2m，深层渗漏量可以忽略。在无降雨、无灌溉时段，式（2.1）可以简化为

该方法有以下优点：①测量精度高，例如中国农业大学石羊河实验站（简称石洋河实验站）新建的称重式蒸渗仪测量精度为0.02mm；②测定面积较大，深度较深，例如石羊河实验站的称重式蒸渗仪表面积为4.0m ² ，深度为2.0m；③同步测定多种指标，通过在蒸渗仪箱体内增加其他传感器，可同时测定土壤水分和温度的垂直剖面或提取溶液进行土壤养分分析。蒸渗仪测定农田ET时需要注意几个问题：

（1）绿洲效应。如果蒸渗仪箱体内作物长势和土壤水分都优于箱外大田作物，或蒸渗仪周围的保护区范围不够大，将会产生绿洲效应。这将会导致蒸渗仪测定值高于实际结果。因此要保证蒸渗仪周围有足够的风浪区，并且箱体内外的土壤水分尽量保持一致。

（2）边框效应。在作物没有完全覆盖地表前，蒸渗仪的箱体外框会反射太阳光，导致箱内的植物接受到更多的辐射，因而测定值偏高。因此，应尽量减少蒸渗仪边框的加热效应，使裸露的边框面积尽量少。

（3）裂缝效应。蒸渗仪箱体内的土壤由于干燥可能会出现裂缝，当采用低频地面灌溉时，灌溉后水分沿着裂缝迅速下渗到底层，导致灌溉不均匀。为了减小裂缝对灌溉不均匀的影响，可以采用少量多次灌溉的方式实现。

2.1.3 涡度相关法测定原理及数据处理标准化方法

2.1.3.1 涡度相关法测定原理及数据处理

Swinbank提出利用涡度相关系统测量温度、湿度、风速的脉动值，从而计算近地层潜热通量（可转化为ET）和感热通量，其公式为

式中：ρ _a 为空气密度，kg/m ³ ；w'为垂直风速脉动值，m/s；q'为比湿的脉动值，g/g； 为垂直风速与比湿脉动的协方差；λET为水汽垂直通量，W/m ² ；ET为作物需水量或耗水量，kg/（m ² /s）；λ为水的汽化潜热，J/kg；H _s 为感热通量，W/m ² ；C _p 为空气的定压比热，J/（kg/K）；T'为虚温的脉动值，K； 为垂直风速与温度脉动的协方差。式（2.5）和式（2.6）表明，只需测量垂直风速与比湿、温度脉动的协方差，便可求出相应的通量。水汽通量除以汽化潜热便得到该时段内的作物ET。

涡度相关法的原理较为简单，但如果直接将观测的数据代入公式进行计算而不对观测数据进行必要的校正和插补，在某些情况下会产生很大的误差。因此，有必要对观测结果进行必要的校正和插补，保证数据的质量。在研究中采用的涡度相关系统为开路式系统。根据前人的研究，数据校正步骤包括：异常值剔除、传感器的频率响应校正、坐标轴旋转、潜热通量的空气密度校正、缺失通量数据的插补。

（1）异常值剔除。异常值通常由随机电信号异常或超声传感器障碍（如降雨）引起。在实际通量计算中，异常值和超出临界值的数据需要剔除，否则会导致通量数据异常。参照通量资料剔除的普适标准，具体的剔除标准为：①异常数据，当某一时刻通量观测数据与前5个时刻通量观测数据的平均值之差的绝对值大于这5个数据方差的5倍时，视为异常数据，予以剔除；②降雨时段数据，降雨时，KH ₂ O水汽计无法正常工作，潜热数据通常显示为“NAN”，这种情况下潜热数据予以剔除；③超出合理范围，制种玉米夜间凝结放热一般不会超过100W/m ² ，因此低于-100W/m ² 的潜热数据予以剔除。

（2）传感器的频率响应校正。传感器对信号波动响应的能力被称为传感器频率响应特征，对传感器的这些局限性的校正称为频率响应校正。频率响应局限可能来源于传感器布置、信号采集和处理。信号频率响应问题在涡度相关系统发展之初就已被意识到，但直到20世纪80年代后期Moore给出了进行频率响应校正的简单算法后才普遍使用 ^[11] 。通过应用频率响应校正可以消除传感器以及物理分离引起的通量误差。测定通量的不准确可能来源于传感器不能对较快和较慢信号波动的准确响应。这种频率响应误差可以通过测定信号谱与实际信号谱的比值来确定。依赖频率的比值函数称为传递函数，某值为0～1.0。频率响应传递函数受传感器时间常数，模拟和数字滤波器的频率响应，信号路径，传感器分离，闭路系统的管道衰减等因素影响。校正频率响应衰减的通量需要知道实际和衰减协谱。

为了进行频率响应校正，经常把理想协谱与频率响应传递函数相乘得到衰减协谱。通过实际模型协谱积分值与衰减模型协谱积分值的比值可以得到通量校正因子。湍流通量的误差ΔF/F，可由式（2.7）得到

式中：T _ws 为传感器垂直风速w和标量s的传递函数，是多个频率响应函数的乘积；S _ws 为w和s关于实际频率f（单位：Hz）的理想协谱。

涡度相关方法的标量s垂直通量的传递函数T _ws 可以表示为实际频率f的函数：

式中：G _w 、G _s 分别为传感器的高、低频损失；T _ss 为由于传感器分离引起的高频损失，例如三维风速仪与水汽计的垂直分离；T _pw 为风速向量线性平均引起的高频损失；T _ps 为标量的线性平均引起的高频损失。

G _w（s） 可以表示为

式中：τ _w（s） 为时间常数，由传感器确定。

根据经验拟合给出了T _ss 一个简单表达式：

式中：n为标准化频率； 为平均水平风速；d _ss 为传感器间的有效侧向距离， d _sa 为传感器实际距离；sin（β _d ）为传感器连线与风向夹角。

风向量线性平均的传递函数T _pw ：

式中：d _pl 为三维超声风速仪的路径长度。

标量线性平均的传递函数T _ps ：

式中：d _s 为三维超声风速仪的路径长度或KH _{2 O} 水汽计的路径长度。

理论计算需要谱模型。已有很多著作给出了湍流标量和矢量的谱和协谱的计算方法 ^[12，13] 。开路涡度相关系统对风速依赖性不太强，主要依赖于传感器路径长度和分离距离。在不稳定和中性大气条件下，频率响应的校正因子与稳定度关系不大。然而，它们在稳定条件下关系很明显的原因可能是协谱变化显著 ^[14] 。

图2.1是典型晴天中午时段（2009-07-14，12：00—12：30）的水平风速（u）、垂直风速（w）、气温（T）和水汽密度（H ₂ O）的功率谱与实际频率的关系图。从图中可以看出，在惯性子区间内（能量传递区），4个属性的功率基本满足能量从低频区按一定的斜率传递到高频区规律。这表明三维风速仪和KH ₂ O水汽计对高频信号的响应能力基本可以满足实际需要。但是，当频率高于2.0Hz时，气温和水汽密度的功率谱曲线呈轻微向上凸起趋势[图2.1（c）和图2.1（d）]，这可能是由于仪器高频损失引起，因此需要进行频率响应校正。

图2.2是垂直风速与气温和水汽密度的协谱关系。从图中可以看出，在惯性子区间内，气度和水汽密度的协谱满足能量从低频区按一定的斜率传递到高频区的规律。类似于图2.1的功率谱，当频率高于2.0Hz时，协谱值也有轻微向上凸起趋势，这表明高频采集区可能有通量损失，需要进行频率响应校正 ^[11] 。

图2.3显示了经过频率响应校正前后的潜热通量和感热通量的比较，下标0表示校正前数据，下标1表示校正后数据。经过频率响应后，潜热和感热通量都有所增加，潜热通量增加6%，感热通量增加2%。这表明频率响应校正对于湍流通量的影响较大，是通量数据后处理过程的必须环节。Wolf等的研究表明，频率响应校正之后，潜热通量增加15%～18%，感热通量增加2% ^[15] ，潜热通量增加量大于感热通量，可能是不同作物以及传感器的原因。Moore等研究对感热通量和潜热通量的校正量为5%～10%，对于某些特定的条件影响更大 ^[11] ，与已有结果类似。

（3）坐标轴旋转。坐标轴旋转的目的是剔除仪器倾斜误差和侧风影响湍流通量矢量成分，以消除平均垂直风速不为0的影响，也就是说使速度和标量浓度梯度仅存在于垂直方向，因而不存在水平平流，也没有风向切变导致的侧风向动量通量。用超声风速计观测风速时面临着怎样安装风速计和确定安装角度等问题。即使地形水平均一，仪器进行严格的水平设置调试，平均气流也不能保证水平。Lee等认为对流效应及中尺度环流、地表不均一和倾斜地表引起的局地热环流以及大气湍流尺度降低阻碍边界层发展是造成垂直风速非零的主要原因 ^[16] 。因此，需要采用适当的坐标轴旋转方法，以消除其影响。目前采用较多的坐标轴旋转方法有：二次坐标轴旋转法、三次坐标轴旋转法以及平面拟合法。

二次坐标轴旋转法的做法为，首先沿着垂直Z轴旋转XOY平面，使得X轴平行于水平风速方向，然后沿着Y轴旋转XOZ平面，使得垂直风速在新的XOY平面上为0。三次坐标轴旋转法在二次坐标轴旋转的基础上，通过进一步旋转使 为0，相应的平均侧风应力为0。二次或三次坐标轴旋转方法在斜坡地形下的通量计算中得到了广泛的应用 ^[12，17] ，然而，二次或者三次坐标轴旋转在长期通量研究中存在着通量平均期间内实际的平均垂直风速可能不为0的不足。对于每半小时的通量数据，采用二次或三次坐标轴旋转法可能会产生显著的偏差或长期通量平衡中的系统低估。

针对此问题，Paw等提出了可以用来估计平均垂直风速的平面拟合（PF）法，表明平面拟合法在比较规则或均匀的平地或坡地能取得较好的校正效果，但在非平坦的下垫面条件下并不适用 ^[18] 。李思恩的研究也表明，应用平面拟合法校正下垫面均一、平坦的农田生态系统通量效果比较理想 ^[19] 。但是其不能用于单个通量的计算，必须有多组通量数据的平均周期才能使用。因此，我们在研究中选用平面拟合法对30min通量数据进行坐标轴旋转校正。具体做法为：通过数据和统计的方法拟合出一个新平面，在该平面上平均垂直风速可以表达为两个方向水平风速的函数且保证垂直风速的平均值为0。图2.4为在坡地上利用平面拟合法进行坐标旋转的示意图。

OZ与水平面XOY垂直，风向如箭头指示。X ₁ O ₁ Y ₁ 为拟合得到的新平面，其中O ₁ Z ₁ 垂直于X ₁ O ₁ Y ₁ ，O ₁ X ₁ 方向是新坐标系的水平方向。于是，平均垂直风速w可以表示为水平风速的线性关系函数，即：

式中：b ₀ 、b ₁ 和b ₂ 为回归参数。

在各生育期选取共10d的30min三维超声风速仪测定的风速资料，采用SPSS软件进行回归拟合，使得函数S的值达到最小：

式中： 、 和 为选取的30min风速资料。

一旦 b ₀ 、b ₁ 和b ₂ 被确定，用平面拟合法校正后的潜热和感热通量值可以按式（2.18）和式（2.19）计算：

式中：w _p 为在新平面上的垂直风速。

P ₃₁ 、P ₃₂ 和P ₃₃ 计算公式如下：

利用平面拟合校正法对中国农业大学石羊河实验站2011年制种玉米生育期内感热和潜热数据进行校正，校正前后的比较如图2.5所示。平面拟合校正后潜热和感热数据均略有降低，但变化不明显，基本可以忽略。这是因为研究区域较为平坦且仪器经过严格的水平设置调试，这与张鑫在西北干旱区葡萄园的研究结果相似 ^[20] 。因此，我们可以得出在西北干旱区平坦农田使用涡度相关系统进行耗水监测不需要进行坐标轴旋转校正。

（4）潜热通量的空气密度校正。当用涡度相关系统观测某气体成分的湍流通量时，需要考虑因热量或水汽通量的传输引起的气体的密度变化。如果测量的是该气体相对于干空气混合比的脉动或混合比的平均梯度变化，不需要进行校正。但是如果测量的是相对于湿空气的混合比，该项校正是十分必要的。因此，对于潜热通量需要进行空气密度校正。Webb等详细阐述了空气密度校正的方法和步骤，对于潜热通量，以协方差形式给出了其表达式[式（2.23）]，通过此式可以校正感热和潜热通量对水汽通量的影响 ^[21] 。

式中：M _d 、M _v 分别为干空气和水汽分子量；ρ _d 、ρ _v 分别为干空气密度和水汽密度，kg/m ³ 。

我们在西北干旱区利用空气密度校正法对2011年制种玉米生育期内潜热数据进行了校正，校正前后的潜热通量相关关系如图2.6所示。从图中可以看出，经过空气密度校正，潜热通量提高大约2%。

（5）缺失通量数据的插补。对于断电、系统故障等因素造成的数据缺失以及异常值剔除后的数据序列，需要建立一套完整的数据插补方法，形成完整的数据序列。目前在通量界常用的数据插补方法有平均昼夜变化法、人工神经网络法和半经验法。我们在研究中数据插补的具体方法为：①缺失的降雨时段潜热通量按0处理；②对于小于2h的缺失数据，采用线性内插的方法进行插补；③对于2h以上缺失的潜热数据，采用P-T系数平均昼夜法插补，对于2h以上缺失的感热数据，采用平均昼夜法插补 ^[22] ，插补周期为7d。

2.1.3.2 涡度相关法测定数据的能量闭合度评价

评价通量观测数据质量的方法很多，包括原始数据分析、谱分析、稳态测试、大气湍流统计特征和能量平衡闭合等，其中能量平衡闭合程度的评价是数据质量评价的重要参照方法之一。其表达式为

式中：R _n 为净辐射，W/m ² ；G为土壤热通量，W/m ² ；λET为潜热通量，W/m ² ；H _s 为感热通量，W/m ² 。

一般来说，涡度相关系统观测的潜热与感热通量之和λET+H _s 小于可供能量R _n G，即普遍存在能量不闭合现象。根据前人的研究，能量不闭合的主要原因有：①通量观测中的采样误差，即涡度相关系统的通量贡献区面积与能量分量测定仪器的测量面积不同造成的湍流能量与有效能量之间的误差；②测量仪器不准确标定带来的仪器系统偏差，以及数据处理不规范对能量平衡的影响；③其他能量吸收项如冠层热储量、植物光合耗能等的忽略；④高频与低频湍流通量损失，理论和实践经验说明应当考虑高频和低频湍流通量的损失，但目前仍没有一种方法能对频率的响应进行准确校正；⑤平流效应的影响，即使在较为平坦的下垫面，当大气层结稳定性很强时由于夜间泄流和平流现象的发生，也会影响能量平衡的闭合程度。

Wilson等详细分析了FLUXNET各站点的能量平衡闭合状况，闭合度分布为0.56～1.20 ^[23] 。Li Z Q等详细评价了ChinaFLUX各站点的能量平衡闭合状况，其平均闭合度为0.73 ^[24] 。本书的研究能量平衡闭合情况如图2.7所示，制种玉米生育期内半小时尺度通量的能量闭合度为0.854，这表明经过数据校正后，涡度相关系统观测的通量数据较为可靠。

另外，在制种玉米全生育期共选取典型晴天和阴天各10d（5—9月每月各2d），分析不同天气情况下的能量闭合状况，如图2.8所示。晴天情况下：λET+H _s =0.854（R _n -G），R ² =0.893；阴天情况下：λET+H _s =0.920（R _n -G），R ² =0.907。晴天的能量闭合率低于阴天。从散点的离散程度还可以看出，晴天情况下不同时刻的能量闭合率差异较阴天情况大。

能量强制闭合被认为是一种提升涡度相关系统测定通量数据质量的可靠方法，通过能量强制闭合可以有效地修正潜热的低估现象 ^[25，26] 。我们在研究中对于白天（R _n ＞0）的通量数据，采用波文比强制闭合法进行能量强制闭合。其具体做法为：①计算出可利用能量与潜热、感热通量和之间的差值D；②将该差值按照每日10：00—15：00平均的波文比β分配给潜热和感热。具体公式如下：

根据Tolk等的研究，夜间的潜热耗能占全天的3%～10%，且其只与气象因子（主要是VPD）有关 ^[27] 。因此，研究中采用Ding介绍的过滤/插值法对夜间（R _n ＜0）潜热通量数据进行处理 ^[28] 。其具体做法为：①当夜间潜热通量小于0或摩擦风速小于0.15m/s时，将对应的潜热数据删除；②在每个生育期分别在剩余的夜间潜热通量和VPD之间建立经验公式；③根据建立起的经验公式对删除的潜热数据进行插补。

能量强制闭合前后制种玉米半小时尺度潜热通量的相关性如图2.9所示，变化过程的比较如图2.10所示。从图2.9可以看出，经过能量强制闭合，半小时尺度潜热通量提升大约8.6%，有效地修正了涡度相关法潜热通量的低估现象。从图2.10可以看出，能量强制闭合对制种玉米不同生育阶段潜热能量的提升幅度不同，在生育初期和生育末期，能量强制闭合对潜热通量的提升效应不明显；在生育中期，潜热通量有显著提升。这说明，在制种玉米生长旺盛期，涡度相关法监测制种玉米潜热通量能量低估现象最为显著。以后在论述中未加特殊说明涡度相关法数据均采用能量强制闭合后数据 ^[29] 。

图2.9 利用能量强制闭合前后制种玉米半小时尺度潜热通量的相关性（2011年）

2.1.3.3 涡度相关法测定数据的空间代表性分析

利用涡度相关法测定通量时，传感器被安置在地表以上的一定高度处，所测得的下垫面与大气之间的物质、能量和动量交换的信息只能代表发生在特定下垫面或某一部分下垫面的物理过程 ^[30] 。通量观测的空间代表性，也就是由点到面的代表性，是指仪器所在点的测量值能在多大程度上反映实际下垫面的平均或累积状况 ^[30] 。与此相关的研究可以应用在通量观测塔的选址，通量观测数据的质量控制和通量观测的尺度扩展等方面 ^[31] 。对于代表性问题的研究现在还没有绝对可靠的定量方法，只能提出一个衡量标准，即实际值与观测值的差在多大范围内变化是可以接受的，对可接受值作统计分析即可做出代表性的大致评价。对于通量观测来说，真实的通量值是非常难以获得的，因此对空间代表性的研究还存在着一定的工作难度和局限性。

通量贡献区（footprint）或源面积（source area）是指对近地面层某一点所观测到的湍流交换过程有贡献的有效源汇区域或者说对观测点的通量大小产生主要影响的下垫面区域 ^[32，33] 。通量贡献区分析可以有效地将对测量值产生重要影响的下垫面区域与源权重很小的不重要区域分开，评价下垫面特征对通量数据产生的影响。通量贡献区的大小与位置随着风向、仪器观测高度、下垫面粗糙度和边界层特征（如大气稳定度等）的变化会发生瞬时的改变 ^{[30，34，35]} 。与通量贡献区相联系的另一个概念是通量贡献函数或称源权重函数，它是描述近地面层表面源汇空间分布和仪器观测值之间关系的函数。通量贡献函数的值可以理解为表面区域某一点源对通量观测值贡献程度的大小。目前在通量界使用的“footprint”既可指通量贡献区，也可指通量贡献函数。通量贡献分析是指对通量贡献函数的求解，可靠的通量贡献函数解可以用来对观测数据进行质量评价，以保证通量来源于所关注的研究区域。通量贡献分析是评价通量观测空间代表性的重要工具。

通量贡献区的概念自提出以来，大量的实验学家致力于发展通量贡献模型并将其应用到实验研究中 ^[36] 。这使得对通量贡献的评价，用基于物理方法的标准代替经验方法成为可能。经典的风浪区长度与测量高度之间的100∶1经验性法则，虽可以作为通量观测中通量贡献区大小的粗略估计，但正逐渐被基于物理方法的标准所替代。但是求解通量贡献函数不是很简单，有很多估算通量贡献函数的模型，包括解析模型、Lagrangian随机质点弥散模型、大涡模拟模型和总体平均闭合模型。另外，这些方法的一部分已经被参数化，为了实际应用的目的简化了原来的算法 ^[37-39] 。关于footprint概念的发展的详细介绍，可以参阅有关文献 ^[34，40] 。我们选用物理意义明确的解析模型FSAM进行分析。

有很多估算footprint函数的模型，因为二维欧拉解析通量源区模型（FSAM） ^[41] 具有数学简单、计算成本合理和精度满意的特征，所以被用来确定我们通量测定的贡献区。FSAM算法是基于平均风速平行并与x轴方向相反的烟羽假定。点（0，0，z _m ）的通量值η可由下垫面上的源强分布Q _η 与footprint函数得到

f（-x，-y，z _m -z ₀ ）将点（0，0，z _m ）的通量η与下垫面上的源强分布Q _η 联系起来，因此把它称为通量贡献函数。该模型不考虑平均风向的湍流扩散，并假定横向风速为高斯分布。源权重函数值可以理解为某一给定点源如（x _s ，y _s ，z ₀ ）对观测点通量值δ贡献的相对权重。因此，源权重函数的大小依赖于点源与某观测点之间的距离。可以认为在点（x _s ，y _s ，z ₀ ）处有一个单位强度的点源，则

式中：Q _η，u 为下垫面上连续的单位源强常数；δ为Dirac-delta分布函数。

因而，把式（2.31）代入式（2.30）卷积中，得到η（0，0，z _m ）的值与源权重函数f成比例的式（2.32）：

源权重函数提供了关于单个点源相对权重的信息。然而，实际中经常需要得到表面的什么区域对高度z _m 处的η值最有影响。换句话说，对η值产生贡献的P水平的最小区域是什么？这个最小区域Ω _p 被定义为水平P的源区。

运用K理论，并参阅参考文献 ^[42] ，浓度扩散的垂直通量F可以表达为

式中：K _C （z）为涡度扩散率； 为侧风向的积分通量；D _y （x，y）为横向的浓度分布函数。

通过下面的二维平流扩散方程，把横向的积分通量 与侧风向的积分浓度 和平均风速廓线 相联系：

积分式（2.34），在z _m 高度的垂直通量表达式为

应用地表横向风积分单位点源 的边界条件： ，式（2.35）被简化为

应用式（2.32）、式（2.33）、式（2.35）和式（2.36），二维通量源权重函数为

式中：F为垂直湍流通量，对应式（2.32）的η；F _u 为表面单位点源（汇）强度，对应于式（2.32）的Q _η，u 。

P水平的通量源区Ω _p 可以用F（x，y，z _m -z ₀ ）=F _P 的等值线限定的区域来表示，即P是footprint函数整体积分φ _tot 的一部分。因此，P水平源区被定义为能达到P水平通量贡献率的最小区域上的通量贡献函数的积分。

式中：φ _P 为通量贡献函数在Ω _p 上的积分。

因为φ _tot =1.0 ^[41] ，所以式（2.38）简化为

对于任何的P _F ，除了P _F =1.0，Ω _pF 的边界被限制在x＞0，因而，应用式（2.36）和式（2.37），式（2.39）被写为

源区Ω _p 依靠有效测量高度z _m （传感器安装高度与零平面位移之间距离）、表面粗糙长度z ₀ 、大气稳定度z _m /L（L为Obukhov长度）、横向风速波动强度σ _v /u*（σ _v 为横向风速的标准差，u*为摩擦风速）。大气稳定度包括不稳定状态（z _m /L＜0），稳定状态（z _m /L＞0）和中性状态（z _m /L≈0）。

式中：T为大气平均温度，K；k为冯·卡门常数，取 0.4；g为重力加速度，取 9.8 m/s ² ； 为感热通量。

输入满足条件的参数，模型输出不同贡献率水平P的通量源区位置和尺度参数值，进而可以绘出P水平的源区范围。

图2.11显示了大气稳定条件下（z _m /L=0.046，L为Obukhov长度）90%通量贡献源区的权重函数分布情况。从图2.11中可知，涡度相关系统测定的通量贡献区在上风向是非正态分布，主要贡献区分布在150m以内。最大通量贡献值出现在上风区大约50m，在主要贡献区的1/3处。测定通量贡献源区在侧向区满足标准正态分布。图2.12显示了玉米全生育期的不稳定和稳定大气层结条件下通量贡献源区的空间分布，稳定大气层结通量贡献区大约是不稳定大气层结的2倍 ^[43] 。

2.1.4 波文比-能量平衡法原理及数据处理

利用波文比-能量平衡法测定作物耗水的理论基础是能量平衡原理与近地层梯度扩散理论。根据能量平衡原理，对于给定下垫面，能量平衡方程可表达为

式中：AD为能量水平交换量，W/m ² ；PH为光合作用耗能；M为冠层储能及作物新陈代谢等耗能，W/m ² 。

对于大面积的均一下垫面，能量垂直交换量远大于水平交换量，AD项可略去不计。此外，PH项和M项远小于其他各项，一般情况下也可以忽略，式（2.42）可简化为

根据近地层梯度扩散理论，单位时间内因湍流交换产生的潜热和感热通量可表达为

式中：K _v 为潜热湍流交换系数，m ² /s；∂q/∂z为近地层空气的比湿梯度，1/m；K _h 为感热湍流交换系数，m ² /s；∂T/∂z为近地层的温度梯度，℃/m。

1926年，Bowen提出用一个比值，即波文比（β）来反映能量平衡中感热通量和潜热通量的比例关系，其表达式为

式中：P _a 为大气压，kPa；λ为水的汽化潜热，J/kg；ΔT为两个高度的温度差，℃；Δe为两个高度的水汽压差，kPa。

根据雷诺相似理论，感热和潜热的湍流交换系数相等，其余项为常数。式（2.46）可化简为

其中，γ为湿度计常数，其表达式为

将式（2.43）和式（2.46）联立可以得出

根据波文比系统测得R _n 、G两个高度的T和e后，计算出两个高度间的ΔT和Δe，将其代入上述各式，便可计算出β、λET和H _s ，并由潜热通量λET除以汽化潜热λ，换算得到作物需水量或耗水量。

波文比-能量平衡法具有较多的理论假设和限制条件，当气象、作物等因素及仪器精度满足不了其条件时，会产生较大的误差，因此需要对其数据进行严格的筛选，选出满足条件的数据，并对不满足条件的数据剔除，选用适当的方法进行插补。我们在研究中按以下步骤进行前期数据处理。

（1）利用方向性判断筛选数据。Perez等通过数学和物理分析，确定了Δe、ΔT与β以及（R _n -G）之间相互制约的定性关系，如图2.13所示 ^[44] 。当λET和H符号不满足要求时，能量传输方向出现错误，数据需要剔除。

我们在中国农业大学石羊河实验站研究中波文比系统在制种玉米全生育期采集数据方向性判断结果见表2.1。从表中可以看出，2011年全生育期共6787组数据，合格数据为 4228 组，合格数据占总数据的62.30%，合格数据比例较低。进一步分析，白天数据（7：30—19：00）共3553组，其中合格数据3324组，合格数据占白天总数据的93.55%，合格比例很高；夜间数据（0：00—7：00 和 19：30—24：00）共3234组，合格数据904组，合格数据仅占夜间总数据的27.95%。夜间不合格数据占不合格数据总数的91.05%，白天不合格数据仅占8.95%。夜间不合格数据较多是因为夜间R _n 较小，（R _n -G）的符号正负交错出现，容易发生错误；且夜间大气多处于层结稳定状态，潜热和感热的湍流交换系数不满足雷诺相似理论，不符合波文比法的应用条件。方向性判断结果表明，白天数据质量较好，夜间数据质量较差。

（2）利用拒绝域判断筛选数据。在波文比数据满足方向性判断的基础上，潜热通量能否准确测定取决于β能否准确测量。尤其是当β值接近-1时，β测量的偏差将会导致潜热通量和感热通量出现严重偏离，而这些异常值是方向性判断无法剔除的。根据式（2.46），β的精度与温湿度梯度的准确测定密不可分，但是以往的研究只是简单地直接给出β的取舍区间，再利用相邻数据进行插补，没有建立β取值范围与传感器精度的关系，缺乏理论基础，随意性较强。针对此问题，Perez等综合考虑水汽压梯度和温湿度探头的辨别率，提出了拒绝域理论。利用该理论，可以确定β取值的动态拒绝区间。在我们的研究中，温湿度探头精度δΔT和δΔe分别为0.1℃和0.04kPa，β的拒绝域是以β=-1±ε曲线围成的一个动态拒绝域，其中：

具体拒绝域类型见表2.2。拒绝域判断结果显示，不合理的数据也主要集中出现在夜间，与方向性判断结果一致。

（3）数据插补。对于不满足方向性判断和拒绝域判断剔除的数据和由于仪器故障导致缺失的数据予以插补。白天缺失和舍去数据时段较短，可利用相邻数据，采用线性内插的方法进行插补。夜间大气多处于层结稳定状态，温湿度廓线非相似性，不符合波文比法的应用条件，不合格数据一般占夜间总数据量的一半以上。随着仪器的使用年限增加，仪器精度降低，夜间不合格数据进一步增多，给数据插补造成很大困难。夜间ET占总ET的3%～10% ^[27] ，夜间潜热数据的缺失将导致ET的低估。

针对此问题，假定对于某一具体区域指定作物夜间累积ET占全天ET的比例是相对固定的，尝试运用波文比法测定的白天ET推求制种玉米全天ET。从表2.3中可以看出，虽然蒸渗仪法和涡度相关法测定的制种玉米全生育期ET分别为 564.69mm和474.17mm，全生育期白天ET分别为514.43mm和433.45mm，两种方法测定的全生育期ET和白天ET结果不同，但是，白天ET占全天ET的比例在不同生育阶段和全生育期是十分接近的：蒸渗仪法在苗期、拔节期、抽穗期、灌浆期、成熟期和全生育期该比例分别为82.35%、90.42%、93.82%、92.14%、93.12%和91.10%，涡度相关法该比例对应为86.19%、91.27%、94.30%、92.52%、90.48%和91.41%。这表明假定是成立的，对于某一具体区域指定作物夜间ET占全天ET的比例相对固定，可以利用其他测定方法（如涡度相关法或蒸渗仪法）得到的白天ET占全天ET的比例和波文比法测定的白天ET值推求波文比法全天ET，结果见表2.4，波文比法测定苗期、拔节期、抽穗期、灌浆期、成熟期和全生育期ET分别为 109.61mm、106.30mm、106.24mm、97.60mm、78.18mm和497.93mm。

2.1.5 茎流计法原理及数据处理

茎流计法可根据原理不同分为：热脉冲法、热扩散法和热平衡法。热脉冲法（heat pulse velocity，HPV）是1932年德国植物生理学家Huber首次提出，并率先运用到实际研究中的。该方法只适用于直径大于 30mm的木本，可以用于分析树干剖面液流的特征 ^[45] ，其准确性在树木研究中得到检验，但树干液流较低时热脉冲技术测定值不准确 ^[46] 。

热脉冲法是在植物茎秆部安装热脉冲发射器（热源），定时发射短时热脉冲，加热汁液，随着植物向上液流热脉冲向上运动，由在热源的上方一定距离处安装的热敏探针T ₁ 探测其温度峰值，确定热脉冲到达时间，测定植物液流速度的。目前为了消除环境温度变化的影响，在植物茎秆下方不受热源影响的地方，安装另一个热敏探针T ₂ ，通过探测温差（T ₁ -T ₂ ）峰值，确定热脉冲到达时间来测定植物液流速度的，假设从加热到两个热敏探针的温差出现峰值的时差为t _e 可用式（2.52）计算热脉冲速率：

式中：V为热脉冲速率，mm/s；X _d 为加热针上方的热敏探针与加热探针的距离，mm，一般为10mm；X _u 为加热针下方的热敏探针与加热探针的距离，mm，一般为5mm；t _e 为加热到两个热敏探针的温差出现峰值的时差，s。

探针对植物茎秆液流有干扰，同时在钻洞和插入过程中对周围组织也有损伤。上述测量方法获得的热脉冲速率V并非未受干扰植物茎秆中的真实的热脉冲速率，而是明显地对植物茎秆液流速率“真值”的低估。为植物茎秆导管中的液流速率HPV可用二次抛物线校正方程：

式中：a、b、c是常数，决定于受影响组织面积的大小，其值则与探针的材料和探针间隔有关。

植物茎秆液流速率SFD可用式（2.54）计算：

式中：ρ、ρ _s 分别为植物茎秆木质部和树液的密度；c、c _s 分别为植物茎秆木质部和树液比热。

ρ _c 的大小可通过对边材样品的简单测量得到，可使用以下几种方法。

（1）由植物茎秆密度与湿度比确定：

式中：ρ _b 为植物茎秆的密度（干重/新鲜时体积）；m为茎秆的含水率；c _w 为植物茎秆的比热，可通过温度T（℃）估算。

（2）由植物茎秆中液体与木质的体积比确定：

式中：F ₁ 、F _m 分别为由植物茎秆中液体与木质占整个茎秆的体积比。

我们采用该方法计算，苹果树分别取0.45和0.5。

植物茎秆液流速率在茎秆的径向不同位置处是不同的，因此计算蒸腾流时要测定茎秆不同深度处的液流速率值，并用式（2.58）计算：

式中：Q _m 为液流量，mm ³ /s；R ₁ 为植物茎秆木质部外边界半径，mm；R ₂ 为植物茎秆木质部内边界半径，mm；r为测定点半径，mm；U' _v （r）为测定点r处的液流速率函数，mm/s，由植物茎秆木质部不同径向测定点的液流速率SFD值拟合得出。

由此即可推求液流通量。实际工作中，对于较细茎秆或枝条的植物，可用标定法确定。对较粗茎秆或枝条的植物有以下方法：①简单平均法，用不同点上测出的SFD的算术平均值乘以木质部切面积（SA）即可得，但精度较差；②抛物线积分法，即假设SFD随植物茎秆径向按二阶多项式变化，求得多项式系数，用上式积分获得树液流量。

随着热技术的发展，先后有两类方法发展起来，即热扩散法与热平衡法。热扩散式探针法（thermal dissipation probe，TDP）是Granier ^[47] 在热脉冲法的基础上进行改造形成的方法，又称恒定热流传感器法（constant-heat flow sensors）或Granier方法。该方法是将一对内置热电偶的探针平行插入树干木质部，上方探针加热，下方探针不加热，作为参考，监测上下两个探针之间的温差。根据温差与液流速率的经验关系确定液流速率 ^[48] 。

热平衡法是 于1973年 ^[49] 提出并用于测定树木的液流，又称为树干热平衡法（trunk heat balance，THB）或 方法。该方法通常是在树干平行插入5个加热片，并在相应位置插入4对热电偶，然后给定功率加热，在忽略树干热容影响的情况下，用热平衡公式计算，直接得到液流通量，无需校正。基于组织热平衡原理，Sakuratani 1981年 ^[50] 发明了一种不用插入加热探针而是以加热片加热的包裹式液流计来测定树干液流，即茎部热平衡法（stem heat balance，SHB），适用于直径较小的植株。该方法经Steinberg和Baker等人发展应用于直径较小的树木或草本，无伤测定，精度较高。该方法在灌木及农作物的液流测定中推广应用 ^[[51-53]] ，都收到了良好的效果。在了解该方法原理的基础上，正确使用该方法所得结果比较准确，在许多研究中已经被证实 ^[53] 。目前已有的包裹式液流计探头直径型号为2～120mm，其中直径60mm以下的探头应用时误差较小 ^[54] 。

包裹式茎流计采用热平衡原理，通过在植物茎秆部施以恒定电压的热量，并通过上下两个热电偶来测定茎流携带走的热量，从而对植物茎秆的液流量与液流速率进行推算并连续监测。其热平衡方程表达式为

式中：P _in 为输入的热量，根据欧姆定律求得，P _in =V ² /R；Q _f 为植物茎秆液流带走的热量；Q _v 为沿树干方向散失的热量。其中，

式中：q _u 为向上散失的热量；q _d 为向下散失的热量。

其计算方程分别为

式中：Kst为茎秆的热传导率，W/（m·K），木本植物一般为0.42；A为茎秆横截面积，m ² ；dT _u /dX和dT _d /dX为温度梯度，K/m；dX是热电偶节点间距，m；dT _u 、dT _d 分别为探头上方和下方的两个热电偶节点间温差。

式中：A、B、H _b 和H _a 分别为探头测定的电压，mV；0.04为热电偶信号转化系数，mV/℃。

从式（2.60）～式（2.64）可得到下列方程：

Q _r 为沿水平方向散失的热量，计算公式为

式中：C-H _c 为探头水平方向两个热电偶的电压差；Ksh为鞘传导率，即水平方向的热传导速率。

Ksh在每天液流为0值时（通常认为在凌晨4：00—5：00）通过热平衡原理进行计算，计算公式如下：

最后，茎秆液流量Q _f （g/s）为

式中：C _p 为水的比热，J/（g/℃）；dT为上下两组热电偶的温度差，其值为：dT=[（A+B）-（H _a +H _b ）]/2。

为简化参数，仪器中取AH=A-H _a ，BH=B-H _b ，CH=C-H _c ，则

上述方法已被广泛地应用于多个树种液流量的测定且准确性已得到证明 ^[55，56] 。许多学者对热脉冲法、热平衡法和热扩散法及其他方法进行对比，均得到较为一致的结果 ^[57-60] 。茎流计法能在保持植株完整、不破坏周围生长环境的情况下测量单株的液流速度，从而较准确地确定单株的蒸腾量。该方法可以在较短时间尺度（小时）上揭示生理和环境因子对单株尺度的蒸腾响应。它的主要优点是不受复杂地形和空间异质性的影响，但由单株尺度提升到农田尺度有一定的难度，需要选取适当的方法进行尺度提升 ^[61，62] 。

2.1.6 遥感结合SEBAL方法测定作物ET原理

遥感传感器可以快捷、周期地获取大范围的二维甚至三维分布的地表电磁波信息，它已越来越广泛地应用在农业、地理、地质、海洋、水文、气象环境监测、地球资源勘探、军事侦察等多个方面。在农业水土工程领域中，遥感数据在水土资源动态变化监测、作物种植面积提取、作物长势监测和估产及作物水分亏缺估计等方面已有广泛的应用。

可见光和红外遥感获取土壤水分状况或通过建立光谱反射率与土壤湿度经验关系来实现，或结合光谱植被指数和热红外遥感获取的陆地表面温度提取土壤水分含量信息，从而进行农田水热时空动态估计，监测农作物亏水程度以指导灌溉。也可以利用遥感信息建立区域作物耗水的遥感反演方法，并根据气象数据计算出潜在耗水，从而对作物受水分胁迫状况出评价。

SEBAL（surface energy balance algorithm for land）是由Bastiaansen等提出的多步骤求取地面特征参数进而得到区域作物耗水量 ^[63] 。利用遥感数据反演得到的地表温度T _r ，半球反照率r ₀ 和归一化植被指数NDVI及其相互关系得出不同地表类型的宽带地表通量（包括感热通量和土壤热通量）后，用余项法逐像元地计算区域ET的分布 ^[64] 。

求解区域的H，需要区域分布的气温、地表粗糙长度及风速等地面资料，一般较难获得。为了解决这个问题，SEBAL方法采用了冷点和热点作为边界条件，假设地表与空气温差δT _a 和地表温度呈线性关系，并利用Monin-Obukhov相似假设对方程进行迭代求解。冷（湿）点，指影像中水分供应充足，植被较为密集的那个像元，在这一点，地表可利用能量完全用作作物耗水，λET≅R _n -G，δT _a ≅0；热（干）点是指干燥且无植被覆盖的像元，此点，作物耗水约为零，H≅R _n -G，δT _a ≅（R _n -G）r _a /ρ _a C _p 。

SEBAL无论作为计算方法还是验证方法已在许多研究中得到了成功地应用。其主要优点是物理基础较为坚实，适合于不同的气候条件区；另外，方法可以利用各种具有可见光、近红外和热红外波段的卫星遥感数据，并结合常规地面资料（如风速、气温、太阳净辐射等）计算能量平衡的各分量，因而可得出不同时空分辨率的ET分布图。运用空气和地表温度的实验方程，消去难以获取的空气温度参数。

通过蒸发比W，并假设蒸发比在24h内大致保持不变，从而将计算得到的瞬时ET扩展为日ET值：

SEBAL方法计算中所需的遥感参数包括：地表反照率α，归一化差值植被指数ND VI，地表比辐射率ε和地表温度T _s ，对于LANDSAT7 ETM+遥感数据，利用SEBAL方法计算区域日ET。

净辐射是地表能量、动量、水分输送与交换过程中的主要来源，可根据地表辐射平衡方程由入射能量减去出射能量求得。通常，净辐射通量R _n 白天为正，夜晚为负。

式中：K _in 为入射短波辐射，W/m ² ；K _out 为出射短波辐射，W/m ² ；L _in 为入射长波辐射，W/m ² ；L _out 为出射长波辐射，W/m ² ；σ为Stefan-Boltzman常数[5.67×10 _{-8 W} /（m ² K ⁴ ）]；ε _a 为大气比辐射率；T _0ref 为参考高度的空气温度，K，可以取灌水充足的植被表面温度；T _s 为地表温度，K。

（1）地表反照率α。计算地表反照率需要各波段的辐射亮度、大气外光谱反射率、各波段的权重系数、大气外反照率和大气单向透射率等参数，其求解过程如下：

式中：L _λ 为λ波段的辐射亮度，W/（m ² ·μm·sr）；Q _cal 为像元灰度值；Q _calmax 为最大DN值，即Q _calmax =255；LMAX _λ 和LMIN _λ 分别为遥感器所接收到的λ波段的最大和最小辐射亮度，即相对应于Q _calmax =255和Q _calmin =0（NLAPS产品）或Q _calmin =1（LPGS产品）时的最大和最小辐射亮度。

式中：α _toa 为大气外反照率；c _λ 为λ波段的权重系数，由λ波段的α _toa 求得，∑c _λ =1。

式中：ρ _λ 为λ波段的大气外光谱反射率；d _r 为日地天文单位距离；ESUN _λ 为大气外光谱辐照度，W/（m ² ·μm）；θ为太阳天顶角，rad。

式中：α _{path-radiance} 为程辐射率，W/（m ² ·μm·sr），其值通常为0.025～0.04；τ _sw 为大气单向透射率。

在晴空且较为干燥的大气条件下：

式中：z为海拔高度，m，可以从DEM数据中获得。

（2）归一化差值植被指数NDVI。

式中：ρ ₄ 、ρ ₃ 分别为近红外、红光波段的反射率。

（3）地表比辐射率ε。比辐射率为物体的辐射出射度与同温度、同波长下的黑体辐射出射度的比值。地表比辐射率通过与NDVI的经验关系求取 ^[65] ：

式中：NDVI＞0。否则，假定ε等于1（如水体）。

（4）地表温度T _s 。地表温度可以由热红外波段通过大气校正法、单窗算法或单通道法反演。SEBAL方法在利用Plank公式由LANDSAT7 ETM+波段6计算出地面物体的亮度温度后，将结果经过地表比辐射率简单校正后获得地表温度：

式中：L ₆ 为LANDSAT7 ETM+第6波段的大气外光谱辐射强度；K ₁ 、K ₂ 为热红外波段校正参数，W/（m ² ·μm·sr），K ₁ =666.09、K ₂ =1282.71（K）。

（5）大气比辐射率ε _a 。

（6）入射波辐射K _in 。

式中：G _sc 为太阳常数，G _sc =1367W/m ² 。

土壤热通量指的是由于传导作用而存储在土壤和植被中的那部分能量，与热流方向的土温梯度、土壤热容量、热传导率成正比，在热量平衡中是一个相对较小的量，直接计算较为困难，一般通过G与T _s 、R _n 、α、NDVI的统计关系求得，根据卫星过境时间进行适当修正：

式中：c ₁₁ 为卫星过境时间对G的影响，过境时间在地方时12点以前c ₁₁ 取0.9，在12点到14点之间取1.0；在14点到16点之间取1.1。

感热通量是指由于传导和对流作用而散失到大气中的那部分能量，是关于大气稳定度、风速和表面粗糙度的函数。其计算公式为

式中：T为空气温度，K；Z为高程。

计算感热通量的公式中，H、dT和r _a 均是未知量，且彼此直接相关，SEBAL中采用迭代方法计算。步骤如下：

（1）根据稳定表面风廓线关系估算摩擦速度的空间分布。稳定表面的风廓线关系为

式中：u _x 为高度z _x 处的风速，m/s；u*为摩擦速度；k为冯·卡门（Karman）常数，k=0.41；z _om 为动量传输表面粗糙度，m，可采用Su提出的经验公式 ^[66] ：

（2）根据气温在一定范围内随高程增加而降低，在研究区内选定某一参考高度z _ref ，利用DEM数据对T _s 进行校正：

式中：z为DEM数据中的高程值，单位应与z _ref 一致。

（3）假设dT与T _s *呈线性关系，建立经验公式：

参数c、d通过在遥感影像上选定两个极端“指示”像元——“干点”“湿点”来确定。“干点”指干燥的天然裸地或没有植被覆盖的闲置农田，假设在干点满足H≈R _n -G ₀ （即，可利用能量完全用于表面加热）；“湿点”指影像中水分供应充足、处于潜在作物耗水水平的像元，一般出现在刚灌水后的农田，在湿点H≈0（即，可利用能量完全用于作物耗水），基于此假设，得到dT值。

（4）计算空气动力学阻力r _a ：

通常，z ₁ 取值略高于植被冠层的平均高度，z ₂ 取值略低于边界层高度，实际应用中，一般取z ₁ =0.01m，z ₂ =2m。

（5）将r _a 代入式（2.86），得到感热通量H。

（6）由于表面加热导致近地层大气处于不稳定状态，SEBAL方法应用Monin-Obuk hov相似理论，引入大气热量传输与动量传输的稳定度订正因子Ψ _h 和Ψ _m ，并计算Monin-Obukhov长度L，对空气动力学阻力r _a 进行校正后，迭代求解H。Ψ _h 和Ψ _m 的具体求解方法如下：

1）L＞0，稳定状态：

2）L＜0，非稳定状态：

3）L=0，中性状态：

式中：g为重力加速度，g=9.81m/s ² ；x（z）为z高度的参数。

（7）将Ψ _{h（z 1 ）} 、Ψ _{h（z 2 ）} 和Ψ _{m（z x ）} 代入下列公式中，对r _a 进行校正：

一般地，重复运行步骤（5）～（7），直到得到稳定的感热通量H。

最后计算得到作物ET。

（1）计算日净辐射量R _n24 。

式中：φ为像元的地理纬度，rad；δ为太阳赤纬。

（2）引入蒸发比，日作物ET ₂₄ 可表达为

式中：R _n24 为日净辐射通量；G ₂₄ 为日土壤热通量；λ为水的汽化潜热。

将ET ₂₄ 的单位换算为mm/d： bwhtNo0cy6gUIe2fSFTRDgPFjEmsbpVuX/mxR6qaYCZsXCO4EzuHAdktrKql3t/k