基于PT-T和SBDART模型的桉树人工林蒸散量研究

杨静学，黄本胜 ，陈亮雄，洪昌红

(广东省水利水电科学研究院，广州510635)

桉树为世界三大速生树种之一，由于其生长迅速、适应性强、耐贫瘠、轮伐期短、经济价值高，成为我国南方重要的人工林经济林种. 广东省是我国最早引种桉树的省份之一，桉树引进已有一百多年历史，尤其在“九五”和“十五”时期在全省范围得到快速发展，现仍呈现增长趋势，数据显示，截止2008年，广东省桉树林已达到107.8 万公顷.

桉树为广东省农民带来巨大经济效益的同时，也引发许多争议. 一方面，水环境学家认为桉树耗水量大，大面积种植会破坏水量平衡，影响水分循环，从而导致林区无法蓄积洪水，致使洪峰期流量增大，枯水期径流减少［1］. 另一方面，林业学家尤其华南地区的林业学家认为，桉树在湿润地区并没有生长优势，也不会消耗过多的地下水，对区域水量平衡没有负面影响［2－4］. 同时认为桉树林大面积种植的确会导致水源缺乏，影响水文过程，但并非因桉树林本身耗水量大，而是因其种植方式不科学所导致的.如桉树人工林种植常采用劈杂、炼山、挖明穴回表土整地方式，破坏了土壤渗透率和水土保持能力，但此类现象存在于各类人工林，并不能以此说明桉树林的生长破坏了种植区的水资源和水环境［5－7］.

针对桉树林在广东省是否存在蒸散发量较大从而对当地水循环产生影响的问题，本文基于遥感影像，结合PT-T 蒸散模型［8－9］和SBDART 辐射传输模式［10］进行了区域(林分)尺度的探讨. 研究以天然混交林为参照，探讨该区域桉树林的蒸散耗水量状况.

1 数据来源

数据包括:1)2006年11月10日TM 影像数据;2)以及桉树林和天然混交林样区矢量边界图. 实验样区的选择应尽量使其地理条件一致，以保证两个样区的计算结果有可比性. 最终选择的实验样区如图1 所示:天然混交林样区位于高明市、桉树林样区位于鹤山市和开平市交界处. 2个实验样区地理条件基本一致:海拔都为100～300 m，且都为阳面，面积分别为1.0 km2和1.7 km2，仅相隔约26 km，气候条件基本一致. 通过高分辨率卫星影像和野外调查可以看出:桉树林有明显的人工耕作混迹;树冠呈有规律的带状分布特征;天然混交林纹理较凌乱(图2). 样区对应的高分辨率影像时间为2006年11月6日，该影像显示，此时期的样区桉树林处于生长期，为了防止用于模型运算的中分辨率TM 遥感影像的样区桉树林处于轮伐期，实验选择的TM影像获取时间为2006年11月10日.

图1 实验样区Figure 1 Experimental sample region

图2 桉树林和天然混交林Figure 2 Eucalyptus forest and natural mixed forests

2 研究方法

遥感技术被认为是可以定量估算区域尺度蒸散发的唯一可能途径［11］. 遥感蒸散发模型的假定:当忽略平流作用时，能量只沿竖直方向传输，若再忽略植被冠层储热、光合作用耗能等分量时，可得到地表能量平衡方程［12］

式中Rn为地表净辐射通量，H 为显热通量，λ ET 为潜热通量，G 为土壤热通量，各项单位均为W/m2.λ ET 为蒸腾和蒸散发作用所消耗的能量，其中λ 为潜热蒸发系数，通常取2.49 ×106W/(m2·mm)，ET是蒸散发量(mm).H 为加热地表而使用的能量.Rn为地表获得的辐射通量能量的总和，可基于辐射传输模式，通过输入影像拍摄瞬时的太阳天顶角、地表温度、反照率和地表反照率等计算.辐射传输模式选用SBDART 模式.

SBDART 模式是近年来国际上比较流行的一种计算辐射传输的模式，该模式基于可靠的物理模型旨在解决在卫星遥感和大气辐射能量平衡研究中遇到的各种辐射传输问题. 模式中应用了辐射传输的离散坐标法来计算辐射传输过程，它采用平面平行辐射，考虑了地表、晴空大气成分、云、气溶胶以及影响辐射过程所有重要光谱波段.

式(1)中，G 相对于H 和λ ET 而言较小可忽略，因此设定G=0，式(1)可变为:

以蒸散发量ET 为目标参数:

式(3)中，由于λ 是常数，求ET 的关键在于求潜热通量λ ET. 基于遥感模型反演λ ET 参数通常通过计算蒸散发比EF 值和净辐射通量Rn来反演:

其中，EF 是指潜热通量与有效能量的比值，当G 为0时，它能非常直观地体现区域的水分和光能消耗效率.

本实验中EF 值的反演采用的遥感模型为PT-T模型［8－9］，它属于半经验模型. PT-T 模型最早由Wang 等［8］基于Jiang 等［13－14］的蒸散发模型提出的，它属于Priestley-Taylor 模型［15］的扩展. Priestley-Taylor 方程为:

其中，Δ/(Δ + γ)的值主要取决于空气温度与高程［16］，Δ 的计算公式为:

其中，

式(5)中，γ 可以由下式求得:

式(10)中z 是海拔高程(m)，p0是海平面的标准大气压(101.315 kPa).

图3 Δ/(Δ+γ)与高程(气温为294 K)和气温(高程为200 m)的相关关系Figure 3 Relationship between Δ/(Δ+γ)and elevation (T=294 K)and Relationship between Δ/(Δ+γ)and temperature (elevation=200 m)

对空气温度和海拔做敏感性分析(图3):当气温不变时(设为299 K)，高程每增加35 m，Δ/(Δ +γ)仅增加约0.001，而当高程不变时(设为200 m)，气温每增加1 K，Δ/(Δ+γ)增加约0.01. 根据以上分析，Δ/(Δ+γ)对高程的敏感程度要远小于气温，由于本实验样区高程为100～200 m 范围内，高程可以视为常数;另外，虽然2 实验样区之间气温2 K 左右的差异使得2 地Δ/(Δ +γ)有约0.02 的差值，但由于2个实验样区面积为1～2 km，在各实验样区内，气温差异不明显，Δ/(Δ+γ)为常数.

图4 PT-T 模型计算α 的散点示意图Figure 4 Calculation of α based on PT-T model and scatter diagram

PT-T 模型计算EF 的关键在于计算α 值. 图4为α 计算示意图. X 轴为NDVI 值(植被归一化指数)，Y 轴为昼夜地表温差值或为地表温度. 该模型遵循3个假设:(1)温度(或温差)最小的CD 边为湿边，该边上所有点的α 值为1.26;(2)温度(或温差)最大的AB 边为干边，该边上点的α 值为0. (3)α 与T 的空间变化关系是线性关系. 因为△HⅠG 与△AGC 为相似三角形，可导出以下方程式:

由于αmax=1.26，αmin=0，式(10)可写为:

由式(14)可知，PT-T 模型中α 的大小主要与地温Ts有关，与NDVI 无关. 地表温度Ts的反演采用单窗算法［17］，精度为0.4 K. 需要注意的是，由于研究区在同一影像上，实验样区各点计算EF 时的Tmax和Tmin为同一数值，因此不具有空间奇异性.

Rn净辐射通量的计算采用SBDART 模型计算，模型输入参数如表1 所示. 其中实验样区气溶胶光学厚度值来源于MODIS 的气溶胶产品MOD04_L2，时间为2006年11月10日. 数据显示，当日天然混交林样区气溶胶光学厚度值为0.223，人工桉树林样区气溶胶光学厚度值为0.244.

表1 SBDART 模型输入参数Table 1 Input parameters of SBDART model

实验计算流程如图5 所示. 分别计算EF 和Rn值，再计算λ ET:

最终通过式(3)计算蒸散发量ET，从而进行实验样区间的蒸散发量对比分析.

图5 实验计算流程图Figure 5 Experimental calculation flow chart

3 结果与讨论

由于TM 热红外波段分辨率为120 m，可见近红外波段为30 m，因此反演的Albedo 和NDVI 分辨率为30 m，地表温度Ts、EF、λ ET 以及ET 分辨率为120 m. 反演的各样区地表温度如图6 所示. 每个样区内的地表温度图都有7 种灰度值，分别代表7个不同的地表温度，从而对应7个不同的EF 值. 将这7个不同的EF 值分别乘以对应区域的Rn均值可获得各区域的λ ET，最终计算出蒸散发量ET. 需要注意的是，此时求得的ET 为瞬时值，即每秒的蒸散发量. 若设其为瞬时均值，将其乘以86 400 可获得一天的蒸散发量. 由于该瞬时值采集时间为当地时间上午10:35，将其设为瞬时均值是偏大的，因为夜晚没有阳光，蒸散发量会相对较弱，因此，以此估算的日均蒸散发量也会相对偏大，但该累积算法不会改变异温像元的蒸散发量之间的相关关系，也不影响本文的对比研究分析.

图6 样区地表温度反演结果图Figure 6 Land surface temperature inversion results of experimental sample regions

为了辅助下一步的分析，计算出7个Ts值区域对应的NDVI 均值、EF、Rn、λ ET 以及ET 值(表2)，表中ETs、ETd分别表示每秒和每天蒸散发量，并通过回归分析，作出各实验样区的NDVI 与EF 以及24 h 的ET 相关关系(图7).

在选择2个样区时充分考虑了地理环境条件一致性的问题，以保证反演参数的可比性. 由表2 可知:(1)样区桉树林的NDVI 值为0.293～0.399，天然混交林NDVI 为0.378～0.527，这说明桉树林样区的植被覆盖率相对天然混交林样区略低;(2)样区桉树林的EF 为0.370～0.558，天然混交林EF 为0.277～0.495，这说明桉树林蒸散发比相对天然混交林偏高;(3)样区桉树林的日均ET 为5.31～8.02，天然混交林为3.98～7.11，这说明桉树林蒸散发量相对天然混交林偏高;(4)随着EF 降低，相应的ET 也越低，这说明蒸散比越低，蒸散量也随之降低.

表2 实验样区统计分析表Table 2 Statistical analysis table of experiment sample regions

从EF-NDVI 和ET-NDVI 的二维散点图(图7)可以看出:天然混交林/桉树林的NDVI 和EF/ET都呈正相关关系，也就是说植被覆盖率越高，蒸散比和蒸散量都越高;从散点图中看出2 类林区的EFNDVI 和LE-NDVI 回归线都呈现平行趋势，而相同林区的EF-NDVI 和ET-NDVI 相关系数差异不大，说明两类林区的区域尺度的蒸散耗水量增加率与植被指数的增加率相当;另外，相同植被指数条件下桉树林区EF 和ET 都比天然混交林高，说明地理环境条件相似时，相同植被覆盖率的情况下，桉树林蒸散发量要大于天然混交林.

在广东省地区，当充分考虑天然混交林和桉树林样区各条件的一致性时，桉树林耗水量大于天然混交林区域. 说明在桉树林在广东省这类湿润气候地区可能依然存在生长优势，也可能进一步影响该区域的水量平衡.

图7 蒸散比EF(A)、24 h 蒸散量ET(B)分别与样区NDVI 的关系Figure 7 Relationships between NDVI and EF (A)，NDVI and 24-h-ET (B)

4 结论

利用遥感数据，基于PT-T 蒸散发模型和SBDART 模型，对广东省桉树林相对于天然混交林区域的蒸散量进行了计算分析. 两类林区的实验样地的选择考虑了各地理环境条件的一致性，从而保证了计算结果的客观可比性. 分析结果显示:区域尺度的植被蒸散发量与植被覆盖呈正相关关系;两类林区的区域尺度的蒸散发量增加率与植被指数的增加率相当;相同植被覆盖条件下，桉树林区蒸散发量要高于天然混交林，当植被指数为0.4 时，24 h 的ET 相差约0.18 mm. 结果表明，桉树林在广东省这类湿润气候地区可能依然存在生长优势，相对于天然混交林而言能消耗更多的水分进行蒸腾蒸散作用，从而对该区域的水量平衡造成影响.

由于迄今未见有基于遥感技术对桉树林区ET进行估算的研究，本研究结果充分说明了利用遥感数据和模型可对桉树林区的蒸散发量进行有效监测. 遥感数据的获取方式和成本相对于实测台站数据更为便捷和廉价，能够对更多的桉树林区进行长时间序列的蒸散发量的监测，这些数据也可以作为缺乏台站数据的区域的有效补充.

［1］张雷，邹进，胡吉敏. 浅析桉树种植对林区水文特性的影响［J］. 林业实用技术，2011(8):11－13.

［2］花忠利. 雷州半岛桉树人工林3PG 模型的研究［D］.福建:福建农林大学，2004:31－35.Hua Z L.Study on 3PG model for euclyptus plantations in Leizhou peninsula［D］. Fujian:Fujian Agriculture and Forestry University，2004:31－35.

［3］李春阳，Tuomela K. 桉树的抗旱性研究进展［J］. 世界林业研究，1998，11(3):22－27.Li C Y，Tuomela K. Advances in research on drought resistance of eucalyptus［J］. World Forestry Research，1998，11(3):22－27.

［4］张宁南. 桉树人工林耗水量研究［D］.北京:中国林业科学研究院，2010:21－78.Zhang N N. Studies on water use of eucalyptus plantations in guangdong［D］. Beijing:Chinese Academy of Forestry，2010:21－78.

［5］韩艺师，魏彦昌，欧阳志云，曹云. 连栽措施对桉树人工林结构及持水性能的影响［J］. 生态学报，2008，28(9):4609－4617.Han Y S，Wei Y C，Ouyang Z Y，Cao Y. Effects of continuous planting rotation on forest structural characteristics and water holding capacity of eucalyptus plantations［J］.Acta Ecologica Sinica，2008，28(9):4609－4617.

［6］罗素梅. 闽西北桉树人工林可持续经营技术研究［D］. 福建:福建农林大学，2010:53－55.Luo S M. Study on sustainable management techniques of eucalyptus plantations in the northwest of Fujian province［D］.Fujian:Agriculture and Forestry University，2010:53－55.

［7］王纪杰. 桉树人工林土壤质量变化特征［D］. 南京:南京林业大学，2011:110－113.Wang J J. Veriation characteristics of soil quality in eucalyptus plantations［D］. Nanjing:Nanjing Forestry University，2011:110－113.

［8］Wang K，Li Z，Cribb M. Estimation of evaporative fraction from a combination of day and night land surface temperatures and NDVI:A new method to determine the Priestley-Taylor parameter［J］. Remote Sensing of Environment，2006，102:293－305.

［9］Yang J，Wang Y. Estimating evapotranspiration fraction by modeling two-dimensional space of NDVI/albedo and day-night land surface temperature difference:A comparative study［J］. Advances in Water Resources，2011，34(4):512－518.

［10］Ricchiazzi P，Yang S，Gautier C. SBDART:A research and teaching software tool for Plane-Parallel radiative transfer in the Earth's atmosphere［J］. Bulletin of the American Meteorological Society，1998，79(10):2101－2114.

［11］冯景泽，王忠静. 遥感蒸散发模型研究进展综述［J］.水利学报，2012，43(8):914－925.

［12］Carlson T N，Capehart W J，Gillies R R. A new look at the simplified method for remote sensing of daily evapotranspiration［J］. Remote Sensing of Environment，1995，54:162－167.

［13］Jiang L，Islam S. Estimation of surface evaporation map over Southern Great Plains using remote sensing data［J］.Water Resources Research，2001，37(2):329－340.

［14］Jiang L，Islam S. An intercomparison of regional latent heat flux estimation using remote sensing data［J］. International Journal of Remote Sensing，2003，24(11):2221－2236.

［15］Priestley C H B，Taylor R J. On the assessment of surface heat flux and evaporation using large-scale parameters［J］. Monthly Weather Review，1972，100:81－92.

［16］Richards J M. Simple expression for the saturation vapor pressure of water in the range－ 50℃ to 140℃［J］.British Journal of Applied Physics，1971，4:115－118.

［17］覃志豪，Zhang M，Arnon K. 用陆地卫星TM6 数据演算地表温度的单窗算法［J］. 地理学报，2001，56(4):456－466.Qin Z H，Zhang M，Arnon K. Mono-window algorithm for retrieving land surface temperature from Landsat TM 6 data［J］. Acta Geographica Sinica，2001，56(4):456－466.