考虑植被盖度的城市蒸散发模拟研究

2022-12-19陈挚黄樱丁金山石喆邱国玉鄢春华

北京大学学报（自然科学版） 2022年6期

陈挚黄樱丁金山石喆邱国玉鄢春华

北京大学深圳研究生院环境与能源学院, 深圳 518055; †通信作者, E-mail: yanch@pkusz.edu.cn

随着全球气候变化, 以城市热岛为代表的一系列生态环境问题日益显著, 给实现联合国“可持续城市和社区”的可持续发展目标带来巨大挑战[1]。作为联系城市水分循环与能量平衡的纽带, 蒸散发是最经济有效的城市温度调节手段[2]。加强城市地区蒸散发的研究, 对缓解气候变化影响和优化城市绿化方案至关重要[3]。

目前对城市样点或斑块尺度植被蒸散发的研究较多, 由于观测方法的限制, 对人类直接生存的邻里空间(neighborhood scale)的城市蒸散发时空分布特征研究有限[4]。涡度相关法(eddy covariance, EC)根据湍流统计理论, 利用风、温、湿等气象要素的脉动数据, 计算相应脉动量的协方差, 进而得到地表显热通量和潜热通量(蒸散量)[4]。涡度相关法假设条件较少, 通过调整架设高度, 可以实现不同范围湍流通量的观测, 是最适合城市不同尺度蒸散发的直接观测系统[5]。1999年, Grimmond 等[6]通过涡度相关法对城市蒸散发进行测量。近年来, 利用涡度系统陆续开展不同城市区域蒸散发观测研究, 探究城乡区域蒸散发差异以及不同区域的蒸散发的时间变化特征, 定性分析气象因素对区域蒸散发的影响[7–9]。国内城市中也开展了涡度相关监测, 但主要用于碳通量的观测研究[10–14]。世界气象组织强调有必要进行更多城市的能量通量观测, 这些观测覆盖更多地点, 能够包括更全面的城市特征和受不同天气影响的城市通量的季节变化[15–16]。

由于设备故障、湍流的随机性、垂直风速与标量之间计算协方差的不确定性以及仪器误差, 会导致一系列随机误差[17], 测量数据会出现明显缺失。在城市复杂环境中, 数据缺失量能达到 30%, 甚至更高[18–19]。不完整的数据集会导致有偏差的不可靠的结果[20]。其次, 计算观测通量的年度或月度预算需要填充数据缺口来获得连续数据集。目前, 通常采用以下方法来填补通量以及气象驱动因素的数据缺口: 非线性回归法[21–23]、查表法[24](look-up tables, LUT)、人工神经网络[25–27](artificial neural networks, ANNs)、随机森林[28–29](random forest,RF)、边际分布抽样算法[30](marginal distribution sampling, MDS)等。这些插值方法基本上只考虑气象驱动因素, 在森林等均质研究区域可获得较好的插值效果[29]。由于城市具有空间异质性, 地表植被零散, 分布不均, 相似气象条件下不同区域的蒸散发有一定的差异, 在异质性城市使用现有方法进行 EC 空缺值插值的准确性有待进一步探究。边际分布抽样算法(MDS)是欧洲通量网(EUROFLUX)中标准的间隙填充方法之一[30], 具有与人工神经网络相近的性能, 比其他算法更加可靠。

本研究基于城市涡动观测系统 2017—2018 每30 分钟观测数据, 结合足迹源区范围及下垫面遥感数据, 获取源区下垫面信息–植被盖度, 构建引入植被覆盖率的城市蒸散发随机森林估算模型(RF_S 模型)。通过与未引入植被盖度的 RF 模型对比, 探究模型估算效果, 并与应用较为广泛的 MDS 插值模型对比插值效果, 同时通过敏感性分析, 探究影响城市蒸散发的潜在因子。深圳市 4—9月雨量充足,为典型湿季, 已有研究发现干、湿季植被蒸散发有显著的差异[31–32]。为了使随机模型可以更精准地学习蒸散发规律, 本研究将数据划分为湿季和干季。

1 材料与方法

1.1 研究区域

研究区域位于广东省深圳市北京大学深圳研究生院内(简称北大园区)。观测仪器安置于北大园区环境与能源学院 E 栋楼顶, 具体位置见图1。以观测点为圆心, 250 m 为半径, 观测区域下垫面如图1(c)所示。一条公路横贯研究区域, 公路以北下垫面以草坪林地为主, 人类活动较少; 公路以南为北大园区, 建筑物和不透水面增加, 人类活动频繁。

图1 研究区域概况Fig.1 Study area

1.2 研究方法

1.2.1 涡动相关数据处理

通过涡动相关系统, 对研究区域的水热通量、净辐射及常规气象要素进行长期观测。该系统主要由三维超声风速温度计、快速响应红外线气体分析仪及净辐射测量系统组成, 以 10 Hz 的采样频率采集传感器高度上的水平风速、垂直风速、温度、水汽密度和二氧化碳密度等。主要观测要素及仪器如表1 所示。通过涡动相关系统内部软件进行数据前处理, 获得每 30 分钟平均湍流数据。

表1 涡动相关系统仪器说明Table 1 Detail information about eddy covariance system

数据处理过程遵循标准程序, 包括野点去除、坐标旋转、频率损失修正、感热通量的超声虚温修正和 WPL 修正。涡动系统是连续测量, 但仪器故障、雨雪条件影响及异常值等会造成数据缺失, 同时非湍流运动会给涡动相关通量计算带来偏差。通过摩擦速度u*的阈值筛选, 剔除低湍流状态数据。输出数据质量标记为“0”的是高质量原始数据; 标记为“1”的较为可靠, 可用于长期通量统计; 标记为“2”和“3”的次之[33–34]。本研究中将经筛选后标记为“2”和“3”的无效数据剔除, 记为空缺值。

边际分布抽样算法(MDS)综合了平均昼夜变化法和查表法两种算法[28]。本研究将该算法与构建模型进行对比, 对每 30 分钟数据空缺值进行插值,分析模型空缺值插值效果。

该方法中的插值算法取决于气象数据太阳辐射(Rs)、气温(Ta)以及饱和水汽压差(VPD)的可利用性。1) 当气象数据中 Rs, Ta 和 VPD 均可获得时, 用相邻7天中相似气象条件下的有效观测值的平均值来填补缺失值(LUT 算法)。相似气象条件定义为 Rs, Ta 和 VPD 差异分别不超过 50 W/m2, 2.5℃和 5.0 hPa。当 7 天窗口不存在相似气象条件时, 查询窗口延长为 14 天。2) 当 Rs, Ta 和 VPD 缺失, 但Rs 观测正常时, 仍然采用上述 LUT 算法, 但相似条件仅定义为 Rs 差异不超过 50 W/m2。3) 当 Rs, Ta 和VPD 同时缺失, 采用平均昼夜变化法算法进行插值。首先采用同一天相邻 1 小时的有效观测数据进行线性插值; 若同一天相邻 1 小时的观测数据缺失,将窗口扩大到相邻 1 天或 2 天。如果执行完以上 3个插值步骤后, 缺失值仍未被完全插值, 则继续增加观察窗口天数, 做进一步的插值。

u*过滤及 MDS 空缺值填充均使用 R 语言的通量数据处理包 REddyProc 实现[34]。

1.2.2 蒸散发源区及植被盖度获取

EC 观测值代表其迎风向一定区域内观测要素的均值, 其观测区域(即代表的下垫面源区范围, 也称足迹范围)随风向不断变化。本研究使用 Kljun等[35]2015年提出的二维足迹模型进行每半小时周期足迹计算。该模型是基于拉格朗日随机粒子弥散模型 LPDM-B 的一种二维通量足迹参数化方法(flux footprint parameterization, FFP), 适用于大范围的边界层条件和不同粗糙度表面[36–37]。

足迹函数f(x,y,zm)用测风累积函数和高斯测风色散函数Dy表示(假设侧风湍流弥散独立于垂直或流向输运):

其中,σy为测风距离标准偏差,X*为缩放的无量纲迎风距离,zm为仪器观测高度(实际观测高度减去零平面位移高度),h为边界c 层高度,u*为摩擦速度,为观测高度平均风速,k为von Karman常数。

足迹计算通过 www.footprint.kljun.net 网站获取标准 FFP 程序包, 输入函数所需参数, 在 Matlab 中实现每半小时足迹函数计算, 空间分辨率为 1 m。模型主要输入参数如表2 所示。

表2 通量足迹模型输入参数详细信息Table 2 Detailed information of input variables for FFP model

参数 umean,u*和σv=可由涡动相关数据处理获得。参数zm和z0根据下垫面元素平均高度zH计算,zm=z-zd=z-fd·zH,z0=f0·zH(fd=0.7,f0=0.1)[37–38]。奥布霍夫长度边界层高度计算通常参考 Kljun 等[35]建议的方法。以上计算所需参数均可由涡度及无人机遥感数据获取。

使用无人机获取覆盖整个研究区域的正射影像, 将图片导入 ENVI 遥感数据处理软件, 通过最大似然法对图片进行监督分类, 将研究区域分为植被覆盖区及非植被覆盖区。图像于 2019年拍摄, 由于深圳市地处亚热带季风气候区, 且研究区域内植被以常绿阔叶植物为主, 该区域内植被年内无明显变化, 同时下垫面构成相对稳定。因此, 基于涡度每 30 分钟周期平均湍流数据, 通过足迹分析获取相应时刻的源区范围。图2 为基于 2018年7月29日夜间 3:00 和日间 14:30 时的数据，计算得到相应的源区范围, 结合下垫面分类图, 即可获得该时刻源区范围内的植被盖度。

图2 蒸散发足迹源区范围及植被覆盖率(S )Fig.2 Source area and vegetation coverage (S)of evapotranspiration

1.2.3 随机森林(RF)模型

本文使用随机森林, 基于半小时周期数据估算城市度蒸散发。随机森林模型是一种基于集成决策树的算法, 可用于分类、回归和聚类[39]。其原理如下。1) 通过重采样, 从训练集N中有放回地随机抽取n(n

模型通过 MATLAB 中 TreeBagger 函数实现。通过均方根误差, 对超参数 mtryMinLeafSize 进行遍历寻优, 最终确定为 5。当决策树的数量为 400~500 时, 模型 OOB 误差不再有明显的波动, 本研究将决策树的数量设置为 500。以往研究也表明, 在该数值下模型即可表现较为优异[40]。

模型所有输入参数信息见表3。为了使随机森林能够更好地学习不同季节和昼夜的蒸散规律, 本研究将湿季和干季区分研究。湿季和干季分别标记为 1 和 2, 白天和夜间分别标记为 1 和 2。基于 2017—2018年每半小时周期数据集, 剔除输入变量存在空缺值以及降雨时数据, 将输出标记为“0”的高质量数据作为原始数据集, 从中随机选取 70%的数据作为训练集(数据量N=5091), 30%的数据作为验证集, 验证模型效果(数据量N=2181)。为了对比植被盖度对城市区域蒸散发的影响, 对包含植被盖度(输入变量M=10)和剔除植被盖度(输入变量M=9)的数据集分别构建模型, 模拟每半小时蒸散发值。

表3 随机森林模型输入参数Table 3 Detailed information of input variables for random forest

本研究使用决定系数R2、均方根误差 RMSE、平均绝对误差 MAE 和偏差百分比 pbias 来评估蒸发模型:

其中,N为样本数量 ;Oi为样本i的实测值 ;pi为样本i的模型估计值; pbias 表示模型倾向于高估或低估了蒸发量, 并表示为观测值的百分比。

2 结果与讨论

2.1 数据质量分析

经标准程序处理及u*过滤[34]后, 2017 和 2018年数据空缺值占比分别为 48.6%和 39.3%。图3 为2017—2018年蒸散发时间分布示意图, 其中夜间数据空缺值比例显著高于白天。2017 和 2018年夜间数据空缺值占比分别为 61.6%和 53.9%。由于仪器故障, 2017–08–16 至2017–09–20 没有观测数据, 2017年数据缺失率高于 2018年。蒸散发年度或月度预算需基于连续的涡动数据, 因此空缺值的填充十分必要。由于夜间湍流状态较为稳定, 夜间涡动系统通量观测值整体偏低[41]。因此, 夜间数据经摩擦速度筛选后, 可用于模型建立与验证的高质量数据集较少[34]。

图3 基于质量控制筛选后的 2017—2018年蒸散发的时间变化Fig.3 Annual and diurnal course of observed evapotranspiration (2017–2018)

2.2 RF_S 模型与 RF 估算模型对比分析

表4 列出 RF 模型输入参数主要气象因子及植被盖度与蒸散发的 Pearson 相关系数。可以看出,植被盖度与 LE 显著正相关, 且相关性呈现明显的干、湿季差异。其中, 湿季植被盖度与蒸散发的 Pearson相关系数为 0.21, 高于风速和土壤含水量。这一结果表明, 将植被覆盖度作为模型输入变量的蒸散发模拟及插值十分必要。

表4 输入参数与蒸散发的Pearson相关性Table 4 Pearson correlation between input parameters and evapotranspiration

图4 为 RF_S 模型及 RF 模型的模拟效果, 可以看出, 无论训练集还是测试集, 考虑植被盖度的RF_S 模型拟合效果均好于传统 RF 模型。对于测试集 , RF 模型的决定系数R2为 0.66, RMSE 为 23.7 W/m2, MAE 为 15.2 W/m2。 RF_S 模型的R2可达0.73, RMSE 为 20.5 W/m2, MAE 为 13.3 W/m2。上述结果表明, 增加植被盖度作为输入参数可以进一步提高随机森林模型对研究区域蒸散发的估算精度。RF_S 测试集中回归拟合线的斜率为 1.06 (接近 1),由图4(b)可知, 散点均匀地分布在 1:1 线两侧, 表明RF_S 的估算结果非常接近涡度实测数据。

图4 基于随机森林模拟的城市蒸散发模拟与实测值比较Fig.4 Performance of the evapotranspiration estimation by random forest

2.3 RF_S 模型与 MDS 模型插值效果对比分析

图5 为 RF_S 模型与 MDS 的估算效果对比。涡动观测系统不仅会在不同时间点出现零散空缺值,受观测条件影响, 还可能产生时间较长的连续缺口。因此将随机选取的测试集 LE 设为空, 对比 RF_S 模型与 MDS 的估算效果(图5(a))。同时, 选取数据质量较高的一段连续时间数据集(2018–05–14 至2018–05–19), 对比长时间序列蒸散发空缺值填充效果(图5(b))。对于随机空缺值填充, MDS 算法的R2为0.57, 明显小于 RF_S 模型(R2=0.73); RMSE 为 26.2 W/m2, MAE 为 17.2 W/m2, 均大于 RF_S 相应统计指标。蒸散发出现连续空缺时, MDS 算法的精度有所下降, 且拟合效果依旧劣于 RF_S。因此, 对于城市蒸散发空缺值填充, 考虑植被盖度的随机森林估算模型明显优于边际分布抽样(MDS)算法。

图5 基于RF_S随机森林模型与MDS算法的城市蒸散发估算值与真实值比较Fig.5 Comparison of estimation effect between RF_S and MDS algorithm

深圳地区呈现典型干湿季特征, 随机森林及MDS 算法在不同时期日间及夜间对蒸散发估算效果如图6 所示。可见, 不同时期日间及夜间随机森林模型整体拟合效果均优于 MDS 算法。对于随机森林模型, 干季拟合效果 (日间,R2=0.70, pbias=0.05%; 夜间,R2=0.41, pbias=–1.0%)优于湿季(日间,R2=0.66, pbias=2.3%; 夜间,R2=0.34, pbias=–14.5%)。日间, 蒸散发估算存在一定的高估现象, 而夜间,存在一定的低估现象。基于实测值和随机森林模型估算值, 获取干湿季日内不同时刻蒸散发均值, 两者保持较好的一致性。对于湿季白天较高蒸散发区间(10:00—17:00), LE 真实值的均值为 70.4 W/m2,RF_S 模型估算值为 67.1 W/m2, MDS 模型估算值为 61.7 W/m2。与真实值相比, RF_S 及 MDS 模型分别低估 4.7%和 12.4%。对于干季白天较高蒸散发区间(10:00—14:00), LE 真实值的均值为 59.3 W/m2,RF_S 模型估算值为 58.4 W/m2, MDS 模型估算值为 59.1 W/m2。与真实值相比, RF_S 及 MDS 模型分别低估 1.5%和0.3%。

图6 干湿季基于MDS、RF_S模型的不同时刻蒸散发均值分析及拟合效果分析Fig.6 Estimation effect of evapotranspiration in dry and wet seasons

综上所述, RF_S 估算模型能获得较高的蒸散发估算值, MDS 模型存在明显的低估现象。湿季,日间 MDS 模型比 RF_S 模型存在更显著低估现象。由于 MDS 仅考虑相似气象条件, 未考虑下垫面类型,因此, 相比 MDS 算法, 通过 RF_S 模型对城市蒸散发进行插值能获得更准确的年度及月度蒸散发量。

2.4 城市蒸散发影响因子敏感性分析

本研究纳入影响蒸散发的因子(如植被盖度),选择表现更优的模型, 直接建立各因子与蒸散发的关系。通过控制其他气象因子不变, 对某一因子变量增加扰动, 获得在多因子共同控制下的某一因子变化对估算值的影响, 进一步分析各因子对蒸散发的相对控制作用, 探讨控制蒸散发的潜在因素。给每个输入变量增加 10%~90%标准差的扰动后, 估算的 LE 值与实际值的R2和 RMSE 如图7 所示。R2和 RMSE 变化趋势表明, 从全年整体上看, 植被盖度(S)、气温(T)和净辐射(Rn)是城市蒸散发估算中最主要的环境因子。扰动前R2和 RMSE 分别为0.73 和 20.5 W/m2, 扰动后植被盖度(S)、气温(T)和净辐射(Rn)的平均R2分别为 0.70, 0.70 和 0.71, 平均 RMSE 分别为 21.6, 21.8 和 21.9 W/m2。大气压(P)和风速(V)影响相对最小, 扰动后平均R2和 RMSE依次为 0.72~0.73 和 20.7~20.8 W/m2。干湿季主控因子存在明显差异。湿季植被盖度(S)是主要控制因素, 增加扰动后R2由 0.74 变为 0.70, RMSE 由 23.3 W/m2变为 24.8 W/m2。其次是大气压(P)、温度(T)和净辐射(Rn), 相对湿度(RH)和风速(V)的影响最弱。湿季深圳雨量充足且净辐射较大, 因此植被盖度是蒸散发主控因子。干季深圳降雨少, 净辐射小,植被盖度(S)的影响明显降低。温度(T)、净辐射(Rn)和饱和蒸气压差(VPD)是主要控制因子, 土壤含水量(SWC)的影响大于干季, 相对湿度(RH)的影响最弱。

图7 随机森林估算模型各输入变量的敏感性分析Fig.7 Sensitivity analysis for evapotranspiration estimations by deep neural network for each input variable over the test dataset