城市化进程所引发的热环境问题是当前“城市气候与环境”的研究热点。在我国西南地区的广西、贵州、云南等存在大面积喀斯特地貌, 构成世界面积最大的喀斯特地区之一^[1]。桂林自古享有“山水甲天下”的美誉, 是该地区典型的喀斯特城市, 特点是在城区中镶嵌着一些喀斯特自然植被覆盖的石灰岩孤峰或峰丛, 岩溶山峰接收的太阳辐射受复杂地形因素的影响。近20多年来桂林城市随着向郊区、向高空的快速扩展呈现越来越明显的热环境问题, 中高层建筑逐渐增多, 建筑物、道路等不透水地表逐渐代替了近郊有植被覆盖的农田、林地以及水塘等自然地表, 一些喀斯特山峰也逐渐进入城区, 但少数城市周边的喀斯特山峰被开辟为采石场, 成为获取城市建设工程石料的便捷途径, 导致喀斯特植被遭到毁灭性破坏, 山体的基岩大面积完全裸露。2014年桂林获准列入“中国南方喀斯特第二期”世界遗产名录, 如何应对城市化进程中的生态环境问题、实现可持续发展成为其城市建设急需研究的内容。

地表水热通量(潜热与显热通量)是描述近地层大气和下垫面间水分、能量交换的参数, 城市地物覆盖变化使得地表物理特性参数发生变化, 影响地表能量平衡和水热通量, 从而影响城市热环境^[2]。Ando和Ueyama^[3]利用涡度相关系统对日本Sakai城市建筑区连续两年的观测数据分析潜热、显热等能量通量的变化规律。中科院地理所匡文慧^[4]研发了城市地表结构组分与热环境生态调控模型(EcoCity), 将潜热通量和显热通量列为模型的10个输入参数之中。研究城市地表水热通量的时空变化, 可从能量角度分析热环境问题形成的物理机理, 为减缓城市热岛提供科学数据。

水热通量的地面观测方法主要有空气动力学法、波文比能量平衡法、涡度相关法^[5-6]和大孔径闪烁仪法^[6-7], 然而其数值的尺度为局地观测点上风向一定范围内, 或大孔径闪烁仪几百米到数公里的较大尺度。对于城市生态系统, 水热通量受地表覆盖类型和环境因素等的影响而高度时空变化^{[3, 8-9]}, 以有限的观测点的局地实验值来描述整个城市将带来较大误差, 而遥感反演方法提供了准确获取整个城市区域水热通量信息的唯一经济可行方法。

遥感反演水热通量的方法有经验方法(包括Penman-Monteith法^[10]、波文比法、水量平衡法、Dalton法等)和物理方法(梯度法、热惯量法等)。经验方法通过拟合通量地面观测值与相关遥感参数的数学方程, 较大程度受局地环境因素的影响。热惯量法不再依赖非遥感数据, 热通量的计算利用地表土壤对吸收辐射的响应方程, 但目前只适用于裸土或沙地情形。梯度法近年来常用有单层模型(SEBAL模型、SEBS模型、METRIC模型等)和双层模型(TSEB模型等), 其首先遥感反演地表温度, 利用阻抗公式结合气温来计算显热通量, 然后潜热通量的计算利用地表能量平衡方程。双层模型的计算过程较复杂, 其将地表分为植被和土壤, 增加了模型相应的输入参数, 需计算大量的阻抗。SEBS模型提出参数化计算kB^-1方法, 当在没有研究区先验了解时存在kB^-1不确定性, 需更多下垫面信息才能合理估计kB^-1。SEBAL模型的物理概念清晰, 运行时只需输入少量气象数据, 研究对象主要为自然生态系统^[10], 近年来模型被尝试应用于中国南京、印度德里等地形平坦的城市区域^[11-12], 由于模型的输入参数未考虑地形因素, 而喀斯特城市中存在喀斯特山峰使得研究区地形复杂, 使用SEBAL模型将导致较大误差。针对SEBAL模型存在的问题, METRIC模型综合考虑了高程、坡度和坡向对辐射的影响^[13], 已成功应用于绿洲、农田、流域等^[14-21], 但未见研究城市区域, 原因主要是我国西南喀斯特位于经济欠发达的西部地区, 缺少项目资金支持对喀斯特城市复杂地形地表水热通量的研究。因此, 针对我国西南地区桂林喀斯特城市近20多年来快速扩展所引发的热环境问题, 改进METRIC模型使其适用于喀斯特城市实际状况, 利用模型和Landsat遥感时序图像反演地表水热通量, 定量分析水热通量的时空变化规律。

1 数据与方法 1.1 研究区概况

研究区位于中国西南岩溶地区的广西东北部桂林主城区, 石灰岩孤峰、峰林广布地面, 是典型的喀斯特岩溶地貌城市。属亚热带湿润季风气候, 气候温和, 年平均气温18—19℃, 最热的7、8月份平均气温约28℃, 最冷的1、2月份平均气温约9℃, 最低气温偶尔降到0℃以下, 年平均相对湿度73%—79%。四季分明, 夏长冬短, 无霜期长, 年平均无霜期309 d。光照充足, 年日照时数1447.1 h。雨量充足, 年平均降雨量1887.6 mm, 降雨量主要集中于4—7月, 年分配不均, 雨热基本同季。全年风向以偏北风为主, 平均风速2.2—2.7 m/s。

1.2 数据

由于桂林位于亚热带气候区, 春季、夏季和冬季多云雾和降雨, 遥感图像中云量太大, 无法准确提供地物信息, 因此遥感数据选取覆盖桂林主城区1994年10月22日TM、2000年10月30日ETM+、2006年9月21日TM、2010年11月11日TM和2015年10月16日OLI/TIRS秋季5景Landsat卫星图像, 经辐射校正、配准、几何精纠正和大气校正。数字高程DEM数据为ASTER GDEM, 空间分辨率为30 m。遥感模型运行所需的卫星过境气象数据来自桂林国家基本气象站。

1.3 METRIC模型原理

地表能量平衡方程是模型的原理^[11]：

(1)

其中LE为潜热通量, 包括植被蒸腾和地面蒸发的能量, LE=λET, ET为蒸散量, λ为水的汽化潜热；R_n为净辐射, A为人为热源释放的能量, H_S为显热通量；G为土壤热通量, ΔS为下垫面热能储存。

(2)

(3)

其中S_in为太阳入射短波辐射, L_in为大气入射长波辐射, L_out为地表发射长波辐射, ε、α、T_S分别为地表的比辐射率、反照率、温度, Stefan-Boltzman常数σ=5.67×10^-8 W/m² K⁴。

(4)

其中T₁、T₂分别为高度z₁、z₂空气温度, ρ_a、r_ah分别为空气的密度、动力学阻抗, 空气定压比热C_P＝1004 J/kg K。

METRIC模型的特点为：(1)假设地-气温差与地表温度T_S呈线性相关方程, 通过选取极端干湿“热点”和“冷点”来计算回归系数, 避免像元空气温度的确定。(2)利用迭代法计算像元显热, 确保热量传输粗糙度、温度梯度和显热通量之间正确的物理耦合关系。最后瞬时潜热通量LE由能量平衡方程计算。

1.4 模型的改进与输入参数的确定

研究区“热点”选取桂林城区中心广场, 为较大面积的裸露水泥地表, 潜热近似为零；研究区“冷点”选取漓江水体中心位置, 显热近似为零, 所有能量都转变为潜热。利用“热点”和“冷点”的地表参数, 经多次循环迭代运算获得地-气温差和地表温度T_S之间的线性相关方程系数^[13]。

模型输入参数改动部分如下, 其他输入参数按Allen等^[13]文献计算。

1.4.1 地表反照率的计算

地表反照率的计算采用Liang算法^[22]：

(5)

(6)

其中, ρ_B、ρ_R、ρ_NIR、ρ_SWIR1、ρ_SWIR2分别为蓝、红、近红外、短波红外1和短波红外2等波段反射率, α为地表反照率, α_toa为大气外反照率, α_path为程辐射(取值0.025—0.04之间), τ为大气单向透射率。

1.4.2 地表比辐射率计算的改进

根据喀斯特城市的实际情况, 提出估算喀斯特山峰混合像元比辐射率的方法。首先从ASTER光谱实验数据库(http://speclib.jpl.nasa.gov)确定典型地物的比辐射率。对TIRS波段10数据, 比辐射率水体ε_W取0.997, 植被ε_V取0.987, 水泥建筑ε_B取0.965, 干燥土壤ε_S取0.968, 石灰岩ε_R取0.958；对TM波段6数据, 比辐射率水体ε_W取0.995, 植被ε_V取0.985, 水泥建筑ε_B取0.968, 干燥土壤ε_S取0.973, 石灰岩ε_R取0.960。Landsat热红外波段的地面空间分辨率TIRS为100 m, TM为120 m, 城市的遥感像元实际上是多种地物构成的混合像元。由于亚热带地区较好的水热条件利于植被生长, 植被基本完全覆盖裸岩间的土壤, 喀斯特山峰在遥感图像中是由植被与裸岩构成混合像元, 比辐射率计算公式为：

(7)

(8)

建筑、道路与绿化植被构成的混合像元的比辐射率计算公式为：

(9)

绿化植被与土壤构成的混合像元的比辐射率计算公式为：

(10)

其中, P_V为植被覆盖度；R_V、R_B、R_S、R_R分别为植被、建筑/道路、土壤、裸岩的温度比率, 可由P_V估算^[23]；dε为混合像元中由于地形起伏引起的比辐射率修正项；地形因子F根据不同几何分布取值^[24]。

1.4.3 地表温度反演的改进

在经典METRIC模型中使用热红外波段通过对比辐射率的校正计算地表温度, 而本文地表温度的反演利用包含大气与地表影响的单窗算法^[25]。对TM/ETM+数据使用热红外波段6, 而对TIRS数据有两个热红外波段, 鉴于波段11的定标存在较大的不确定性, 故使用波段10^[26]。

(11)

(12)

(13)

其中, L_λ、T_B分别为热红外波段的辐射亮度、亮度温度；k₁、k₂为常数；a、b为常量^[25], T_S为地表温度, T_a为大气平均作用温度, 可由地面附近气温计算。

1.4.4 土壤热通量和热能储存、地表粗糙度的计算

参考已有的研究, 利用经验公式对土壤热通量G进行估算^[27]：

对植被覆盖下垫面

(14)

对非植被下垫面

(15)

其中, NDVI为归一化植被指数。系数C_g取值, 建筑和道路、水体取0.40, 裸土、裸岩取0.30。城市下垫面建筑和道路等的热能储存, 参考以往的研究用其占净辐射比例合并到土壤热通量中进行估算。

地表粗糙度参考已有文献的研究结果^[27]。建筑取1.0 m, 以林地为主的植被取0.8 m, 水体取0.000034 m, 裸土取0.001 m。

2 结果与分析

由于桂林是我国西南地区的旅游城市, 主城区基本没有大型工业厂房；又属于我国西部地区经济欠发达的三线城市, 汽车数量较少, 汽车交通流密度相对不太高；遥感图像获取时间在气温舒适的秋季, 一般不使用空调, 因此人为热源释放的热量较少, 远小于太阳辐射, 在能量平衡方程中可以近似忽略。

使用机器学习方法支持向量机法对遥感图像进行地物分类, 然后根据分类结果估算地表比辐射率^[26]。在遥感图像中根据地面GPS定位数据提取典型地物样本, 分类特征向量的元素由归一化植被指数NDVI、归一化建筑指数NDBI、水体指数MNDWI、光谱特征和数字高程DEM组成。支持向量机SVM的内积函数取径向基RBF函数, 使用交叉验证法确定最优惩罚系数和间隔松弛因子。SVM的学习使用样本训练集, 分类精度使用测试集进行检测, 将地物分为建筑(道路)、植被、水体、裸土、裸岩。2015年、2010年、2006年、2000年、1994年的总体分类精度分别为89.6%、88.7%、89.3%、87.5%、88.2%, 接近90%, Kappa系数分别为0.881、0.856、0.877、0.827、0.843。

使用改进的METRIC模型反演桂林城区地表显热和潜热通量, 结果见图 1、图 2。反演结果的验证, 利用BREB法对2015年卫星过境时广西师范大学校园内草地上架设的HOBO小型自动气象站、四分量净辐射传感器CNR4的空气温湿度梯度、辐射等气象数据估算地表显热与潜热通量, 根据GPS数据定位, 模型值与观测值的相对误差见表 1, 模型改进前显热、潜热通量反演值与观测值的相对误差分别为16.5%、17.8%, 模型改进后显热、潜热通量反演值与观测值的相对误差分别为14.4%、13.7%, 比模型改进前精度有所提高, 反演误差与已有的研究相近^{[14, 17]}, 在该模型正常误差范围之内, 说明模型估算结果是合理的。

2.1 水热通量的空间变化

统计研究区典型地物水热通量均值和波文比(显热通量与潜热通量比值)均值, 结果见表 2、表 3。

由图 1、图 2、表 2、表 3明显看出, 水体和植被覆盖地表的潜热远高于显热, 水体的潜热高于植被覆盖地表的潜热而显热低于植被覆盖地表的显热, 水体的显热低于100 W/m², 波文比数值范围：水体在0.10—0.25, 地面植被(以林地为主)和喀斯特山峰阳坡植被在0.35—0.50, 喀斯特山峰阴坡植被在0.30—0.45；在太阳辐射下受地形影响, 潜热和显热通量均值以喀斯特山峰阳坡植被、地面植被明显高于喀斯特山峰阴坡植被。喀斯特山峰植被像元实际上是包含少量比例裸岩的混合像元。水体(穿城而过的漓江和城区的桂湖、榕湖、杉湖等)的潜热最高而显热最低；而城区茂密植被覆盖的较大型喀斯特山峰(普陀山、叠彩山、西山、穿山、骝马山、老人山、南溪山等)的潜热也相对较高而显热相对较低。

建筑/道路、地面裸土和喀斯特山峰裸岩的显热高于潜热, 波文比数值范围：建筑/道路在1.10—1.65, 地面裸土在1.10—1.80, 喀斯特山峰裸岩在1.40—2.50。显热通量均值以喀斯特山峰裸岩和建筑/道路高于地面裸土, 潜热通量均值以建筑/道路和地面裸土高于喀斯特山峰裸岩, 原因是建筑/道路周边常种植树木等绿化植被, 建筑/道路实际上是包含少量比例绿化植被的混合像元, 植被蒸腾产生潜热；研究区中地面裸土面积不多, 主要是建筑工地开挖土方, 常用浇水降尘措施减少扬尘对周边环境的影响, 裸土中具有一定的含水量, 土壤蒸发产生潜热；喀斯特山峰裸岩间土壤土层瘠薄、含水量低, 植被的生长速度十分缓慢, 故实际上是包含少量比例草本植被、低矮灌木的混合像元, 植被蒸腾和土壤蒸发量较少。

2.2 水热通量的时间变化

研究区水热通量随时间的变化受政府城市新开发区建设和绿化建设导致的地表覆盖类型变化的综合影响。由表 2可以看出研究区波文比在1994年最高, 达到1.62, 原因是1991年国务院批准建设桂林国家级高新技术产业开发区, 1994年广西政府批准建设桂林西城经济开发区, 带动了桂林城区的快速扩展。

然而桂林作为旅游城市, 新开发区的绿化建设也同时受到政府的重视, 常使用大树移植方式快速提高新开发区的植被覆盖率, 加上水热条件较好, 樟树、榕树等本土树种生长速度较快, 1999年桂林被命名“全国园林绿化先进城市”, 2000年研究区波文比下降到20年来的最低点, 为1.24。

2003年建设高新区铁山工业园和信息产业园；2005年建设高新区英才科技园；为保护漓江水环境, 2007年8月临桂新区获批建设计划到2015年完成13.96 km²建城区目标；2010年开始建设临桂新区秧塘工业园(占地10 km²), 新开发区的显热较高, 而潜热较低, 2000—2015年研究区波文比呈现逐渐升高的趋势, 2015年上升到1.51, 接近1994年数值。

分别将显热、潜热通量的数值范围从低到高的40%、30%、40%定义为低值区、中值区、高值区, 统计各区间像元比例, 结果见表 4、表 5。研究区显热高值区和潜热低值区主要为植被覆盖率非常低的较大面积的新建高密度建筑区和喀斯特山峰裸岩, 比例低于10.0%；显热低值区和潜热高值区主要为水体(漓江、桂湖、榕湖、杉湖等)、茂密植被覆盖的喀斯特山峰、城市公园林地等, 受到人为保护, 比例较高, 显热低值区比例在35.0%—50.0%, 而潜热高值区比例在35.0%—55.0%。喀斯特城市扩展过程出现的显热高值区比例以1994年最高, 为10.0%, 2000年下降到5.4%, 之后呈逐渐上升到2010年的9.4%, 但2015年下降到7.1%, 而潜热低值区比例的变化趋势与显热高值区比例基本相同。显热高值区和潜热低值区比例的变化导致显热中、低值区比例和潜热中、高值区比例的变化, 城市扩展将显热和潜热中值区的低密度建筑区, 显热低值区和潜热高值区的农田、林地、池塘等, 转变为显热高值区和潜热低值区的新建高密度建筑区, 部分新建高密度建筑区经过若干年的绿化建设后植被覆盖率增加可以转变为显热和潜热中值区。

2.3 植被覆盖度变化对水热通量的影响

城市地表水热通量受地物覆盖变化、气候条件等影响, 植被覆盖度变化是最重要的因素之一, 将其划分为10等份, 统计各等份水热通量均值, 结果见图 3。统计植被覆盖度增加0.1时水热通量变化的均值, 结果见表 6、表 7。植被覆盖度与显热通量呈明显负相关关系, 而与潜热通量呈明显正相关关系。在不同植被覆盖度范围, 水热通量随植被覆盖度的变化程度存在差异, 植被覆盖度在0.0—0.1范围极稀疏和在0.8—1.0范围极茂密时其变化对水热通量的影响相对较弱, 植被覆盖度增加0.1, 显热降低、潜热升高4—10 W/m²；而植被覆盖度在0.1—0.8范围时其变化对水热通量的影响相对较强, 水热通量与植被覆盖度近似呈线性关系, 植被覆盖度增加0.1, 显热降低8—27 W/m², 而潜热升高8—24 W/m²。

3 讨论

使用METRIC模型反演桂林喀斯特城市地表潜热、显热通量时, 由于实验条件限制, 未能对不同类型下垫面的土壤热通量、热能储存、地表粗糙度、零平面位移等进行实验观测, 而参考已有的研究进行估算, 将导致一定的反演误差；模型验证时尽管潜热、显热通量遥感模型像元反演值与点上实验观测值由于所代表的空间尺度不同, 存在验证误差。但由于反演误差主要影响地表通量数值, 对地表通量空间整体分布趋势、显热通量与潜热通量比值的波文比、水热通量与植被覆盖度关系的影响较小。

桂林喀斯特城市城区中镶嵌着一些拔地而起的喀斯特孤峰, 常被设置为公园, 喀斯特山峰上的植被受人为保护而生长状况较好, 植被覆盖度普遍相对较高, 植被覆盖阻挡、反射太阳对石灰岩山体的部分直接辐射, 茂密的喀斯特植被叶面整体的蒸腾作用较大而吸收较多热量, 潜热较高而显热较低, 影响其上方及周边空气, 形成局地小气候环境, 对调节喀斯特城市微气候环境具有重要作用。

喀斯特山峰裸岩一般是采石场, 其上以灌丛和阔叶林为主的喀斯特植被遭受人为破坏, 形成大面积的喀斯特山峰基岩裸露。由于裸岩间的土壤土层瘠薄, 植被生长缓慢, 很难在短时间内恢复喀斯特山峰植被覆盖, 容易造成土壤侵蚀, 形成石漠化现象。而热容量较大的大面积裸露石灰岩在太阳直接辐射下吸收、贮存大量的热量, 显热较高而潜热较低, 对城市热环境产生不透水地表相似的影响效应。因此镶嵌于城区中喀斯特山峰的植被保护对喀斯特城市热环境的改善至关重要。

植被覆盖率非常低的较大面积的新建高密度建筑区种植幼树等少量的植被, 植被覆盖率非常高的较大面积的林地继续增加植被覆盖, 对降低显热和提高潜热的效果不显著, 对环境状况的改变影响不大。

在太阳辐射下喀斯特山峰的水热通量受地形因素的影响, 但桂林主城区的喀斯特山峰距地面的相对高度一般不超过200 m, 高程的影响较小。

论文的创新点在于根据我国西南地区桂林喀斯特城市城区存在喀斯特山峰复杂地形的实际情况改进METRIC模型, 使模型适用于我国西南地区喀斯特城市地表水热通量的反演, 使用地面观测数据进行验证, 分析喀斯特城市地表水热通量的时空变化规律, 迄今未见相关研究论文发表。

4 结论

(1) 水体和植被覆盖地表的潜热通量远高于显热通量, 而建筑/道路、裸土和喀斯特山峰裸岩的显热高于潜热；潜热通量从高到低依次为水体、喀斯特山峰阳坡植被、地面植被、喀斯特山峰阴坡植被、建筑/道路和裸土、喀斯特山峰裸岩, 显热通量从高到低依次为喀斯特山峰裸岩和建筑/道路、裸土、喀斯特山峰阳坡植被、地面植被、喀斯特山峰阴坡植被、水体。

(2) 喀斯特城市水热通量随时间的变化受政府城市新开发区建设和绿化建设导致的地表覆盖类型变化的影响。研究区波文比在1994年最高, 达到1.62, 2000年下降到20年来的最低点, 为1.24, 之后逐渐升高到2015年的1.51, 接近1994年数值。

(3) 喀斯特城市扩展过程出现的显热高值区和潜热低值区比例低于10.0%, 其变化引发显热中低值区和潜热中高值区比例的变化。显热高值区比例在1994年最高, 为10.0%, 2000年下降到5.4%, 之后至2010年逐渐上升到9.4%, 但2015年下降到7.1%, 潜热低值区比例的变化趋势与显热高值区比例基本相同。

(4) 植被覆盖度在不同范围对水热通量的影响程度存在差异, 植被覆盖度在0.0—0.1范围很低和在0.8—1.0范围很高时其增加0.1, 显热通量降低、潜热通量升高4—10 W/m², 影响相对较弱；而植被覆盖度在0.1—0.8范围时其增加0.1, 显热通量降低8—27 W/m², 而潜热通量升高8—24 W/m², 影响相对更显著。

(5) METRIC遥感模型提供了定量分析喀斯特城市水热通量时空变化的快速、经济可行的方法, 并且取得较好的结果。

致谢: 本研究得到了“广西八桂学者创新研究团队”、“广西区域多源信息集成与智能处理协同创新中心”的支持, 特此致谢。