在全球气候变化和快速城市化背景下, 城市洪涝灾害风险逐渐加剧^[1-2]。研究表明, 全球超过600个城市正面临着严重的洪涝灾害风险。其中, 中国有超过60%的城市发生过内涝^[3-5]。城市化过程中土地利用/覆盖的剧烈变化尤其是不透水面的迅速扩张, 导致城市区域自然水文过程明显改变, 削弱了城市生态系统的雨洪调节能力^[6-8]。

城市生态系统雨洪调节服务是人类从生态系统获取的重要惠益之一, 反映了城市生态系统通过植被的冠层截留、土壤下渗等方式影响地表径流产生过程和径流量, 从而减小洪涝灾害风险与造成社会经济损失的能力^[9-11]。目前, 较多研究集中于区域尺度的生态系统洪涝调节服务空间评价与制图^[12-14], 但由于城市内部下垫面结构具有高度的空间异质性, 城市雨洪调节服务的量化评估存在一定困难。在城市流域的研究中, SWAT、SWMM、HEC-HMS等分布式水文模型被广泛应用于模拟不同降雨情景下地表径流的变化特征^[15]。为了评估城市雨洪调节服务, 前人基于地表径流构建了一些表征雨洪调节服务的参数指标, 包括城市绿地径流调蓄量和调节率、地表径流下渗率、地表径流系数等^[16-19]。相关研究表明自然植被较其他土地利用类型能够提供更多的雨洪调节服务, 城市公园绿地区较工业区、商业区等功能区具有更高的雨洪调节能力^{[13, 16]}。同时, 由于城市不透水面比例的上升使得地表产流量增加和洪峰提前, 绿地空间的减少导致城市对雨水的滞留能力下降, 这一过程加剧了城市雨洪调节服务的退化程度^[19-20]。因此, 揭示北京城市土地利用/覆盖变化对雨洪调节服务的影响对优化城市生态结构与提升城市生态系统可持续性尤为重要。

近年来, 快速城市化导致北京城市不透水面迅速增长, 进一步加剧了城市雨洪产流和内涝风险。例如, 2012年“7.21”的特大暴雨引发的城市内涝灾害, 造成了严重的社会经济损失^[21-22]。虽然已有研究评估了北京城市绿地空间的滞蓄雨水径流调节能力与雨洪产流风险^{[3, 15]}, 但对长时间序列的土地利用/覆盖变化特征及其对雨洪调节服务影响的研究仍然较少。因此, 本研究通过分析1991—2014年北京城市土地利用/覆盖和雨洪调节服务的变化, 定量评估城市土地利用/覆盖变化对雨洪调节服务的影响, 以期为城市生态结构优化、雨洪调节服务提升及海绵城市规划建设提供科学依据。

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

北京市地理位置为39.4°N—41.6°N、115.7°E—117.4°E, 市域面积约1.64万km², 其中城市核心区(东城、西城、朝阳、丰台、石景山、海淀)的面积为1369 km²。过去30年间, 北京市经历了快速的城市化过程。到2018年, 全市常住人口为2154.2万, 人口城镇化水平达到86.5%^[23]。北京市的气候为温带大陆性季风气候, 雨热同期, 多年平均降水量约590 mm, 其中夏季降雨量占全年的70%以上^[24]。2016年, 住建部将北京市列入第二批“海绵城市”试点城市, 旨在通过提升排水能力和建设低影响开发设施以提升城市雨洪调节服务, 达到“小雨不积水, 大雨不内涝”的成效^[25]。

1.2 数据源

本研究所用的数据包括遥感影像数据、土地利用/覆盖数据和土壤类型数据。(1)Landsat遥感影像(行列号：123/032), 包括1991、2000、2009和2014年4期, 空间分辨率为30 m, 数据下载自美国地质调查局(https://glovis.usgs.gov/), 用于提取城市土地利用/覆盖数据；(2)高分辨率遥感影像, 用于验证城市土地利用/覆盖数据, 时间为2014年, 空间分辨率为0.8 m, 数据获取自地理监测云平台(http://www.dsac.cn/)；(3)土地利用/覆盖数据包括北京市市域1990、2000、2010和2015年的1:10万土地利用/覆盖数据, 该数据城市用地的分类精度能够达到90%以上^[26-27], 将用于辅助提取城市用地边界, 空间分辨率为30 m；(4)土壤类型数据则作为水文模型的输入数据, 空间分辨率为1 km。土地利用/覆盖数据和土壤类型数据均获取自中国科学院资源与环境科学数据中心(http://www.resdc.cn/)。

1.3 城市土地利用/覆盖数据提取 1.3.1 城市用地提取

基于北京市1990、2000、2010和2015年土地利用/覆盖数据的城市用地图层, 利用1991、2000、2009和2014年的Landsat影像和人工目视解译方式在ArcGIS 10.2软件环境下对4期城市用地边界进行修订, 获取对应年份的城市用地数据, 即城市建成区边界数据。为了评估城市土地利用/覆盖变化对雨洪调节服务的影响, 本研究将各个时段的城市用地划分为老城区和新扩张区, 具体见图 1。

1.3.2 城市地表覆盖组分数据提取

本研究将城市土地覆盖划分为不透水面、绿地空间、水域和裸土四种类型, 并利用线性混合像元分解方法提取各类的组分信息, 该方法目前已经广泛应用于城市不透水面信息提取研究中^[28-31]。线性混合像元分解是假设像元在某一波段的反射率等于各个组分的反射率与其所占像元面积比例的加权和, 具体可表示为：

(1)

式中, R为地表反射率；f_i为端元e_i组分信息；n为选取的端元数目。本研究中选取高反照、低反照、植被和裸土4个端元进行像元分解, 并用于后续不透水面、植被、水域和裸土信息的提取^[30]。为验证城市土地覆盖组分的提取精度, 在2014年建成区内随机采样了187个90 m×90 m的网格, 利用2014年的高分辨率影像数据和人工目视解译方式获取各个网格内不透水面和绿地空间比例。通过对比基于Landsat提取的不透水面和绿地空间比例与高分辨率遥感影像提取的地表覆盖比例, 均方根误差小于10%, 表明获取的数据精度能够满足后续研究需求(图 2)。

1.4 城市雨洪调节服务评估

本研究选取城市绿地空间径流调节率^[9]表征地表径流调节能力, 用于刻画城市雨洪调节服务。具体计算公式如下：

(2)

式中, ROR为各区域的地表径流调节率；ΔV为各区域绿地空间的地表径流调蓄量；P为降水量；A为各区域的面积。径流调节率越高, 则表明绿地空间的提供的雨洪调节服务越高。

对于绿地空间的地表径流调蓄量^[9](ΔV), 具体计算公式如下：

(3)

式中, ΔV是像元内绿地空间的地表径流调蓄量；Q_i是比例为100%的不透水面产生的地表径流量；Q_a是实际产生的地表径流量；A_i是每个像元的面积。

本研究利用SCS—CN水文模型^[32]模拟不同降雨情景下地表产流(Q), 具体计算公式如下：

(4)

(5)

(6)

式中, Q是地表产流量；P是降水量；S是降水潜在最大滞留率；CN是表征地表产流能力的参数, 取值范围是0—100；I_a是降水的初始下渗量；λ是土壤的下渗系数, 一般取值为0.2^[15-16]。

CN是SCS—CN模型最关键的参数之一, 本研究参考Fan等^[33]提出的方法进行计算, 即利用城市地表覆盖比例及其各自CN值进行加权计算得到, 具体计算公式如下：

(7)

式中, f_isa, f_veg和f_soil分别代表不透水面、植被和裸土在像元内的百分比, CN, CN_isa, CN_veg, CN_soil分别代表计算得到的综合CN值及不透水面、植被和裸土的CN值, 不透水面和裸土的初始CN值分别为98和91, 植被的初始CN值的计算方法见参考文献^[32-34]。此外, 土壤湿度设置为中等情景, 降雨情景设置1、10、25a和100a一遇4种情景, 北京市对应的降雨强度分别为39.7、72.1、84.9 mm和104.4 mm^[15]。

1.5 城市雨洪调节服务空间集聚特征分析 1.5.1 Moran′s I空间自相关分析

全局Moran′s I^[35]可以测度区域整体的自相关性, 本研究选取该指数评估城市雨洪调节服务空间分布的集聚特征, 其计算公式如下：

(8)

式中, x为k个区域地表径流调节率的均值, W_ij为空间权重矩阵, x_i为第i个区域的地表径流调节率, 且i ≠ j。

1.5.2 冷热点分析

冷热点分析可用于评估雨洪调节服务高低的聚集状况, 具有显著统计意义的热点区, 其本身属性有高值, 其周边也为高值要素。Getis-Ord Gi^*指数^[36]是刻画区域冷热点的常用指标, 其计算公式为：

(9)

式中, G_i(d)是第i个街道的集聚特征统计量；x是第i或j个街道的地表径流调节率；w_{i, j}是基于距离d范围内的空间相邻权重矩阵；S_i²是地表径流调节率的方差；是所有街道径流调节率的均值。

1.6 城市土地利用/覆盖变化对雨洪调节服务的影响评估

为定量评估城市土地利用/覆盖变化对雨洪调节服务的影响, 本研究首先计算t1到t2时段老城区和新扩张区绿地空间径流调蓄量变化, 然后计算其占城市用地总径流调蓄量变化的比例, 以表征老城区和新扩张区土地利用/覆盖变化对雨洪调节服务的变化贡献。具体计算公式如下：

(10)

(11)

式中, C_old和C_new分别为t1到t2时段老城区和新扩张区土地利用/覆盖变化对雨洪调节服务影响的贡献；V_t1和V_t2分别为t1和t2时间的径流调蓄量；V_{t1_old}和V_{t2_old}分别为老城区在t1和t2时间的径流调蓄量；V_{t1_new}和V_{t2_new}分别为新扩张区在t1和t2时间的径流调蓄量。

2 结果与分析 2.1 城市土地利用/覆盖变化特征分析

在人口和社会经济增长的驱动下, 北京城市用地呈快速扩张趋势, 但增长趋势逐渐放缓。1991—2014年, 城市用地面积从380.44 km²增长1476.96 km², 扩张速度由1991—2000年的54.08 km²/a下降到2009—2014年的39.22 km²/a(图 3)。空间上来看, 城市用地主要向主城区的东部和南部扩张, 新增城市用地多分布在朝阳区、海淀区和丰台区, 并呈圈层式向外扩展的模式(图 4)。

同时, 城市内部土地覆盖结构也发生显著变化。其中, 城市建成区内不透水面的面积由1991年的204.22 km²增加到2014年的937.96 km², 比例则由53.68%上升到63.51%(图 3)。特别是2000年以来, 在新一轮国土空间开发的驱动下, 北京城市内部不透水面也呈快速扩张态势, 2000—2014年间增长了308.35 km², 比例上升了12.67%。虽然城市公园绿化建设取得一系列成效, 使得公共绿地空间的面积由139.93 km²增长到483.74 km², 但在整个建成区内绿地空间比例却下降到32.75%(图 3)。

2.2 城市雨洪调节服务变化特征分析

总的来看, 1991—2014年北京城市的雨洪调节服务呈先上升后下降趋势。在1年一遇降雨情景下, 2014年北京城市建成区的平均地表径流调节率为34.05%, 表明城市绿地空间能够有效滞留34.05%的雨水, 一定程度上降低城市的雨洪产流风险。随降雨强度增大, 城市绿地的径流调节能力逐渐下降, 10a、25a和100a一遇降雨情景下平均径流调节率分别下降到26.66%、24.54%和21.89%(图 5)。与1991年相比, 城市建成区内平均地表径流调节率在4种降雨情景的下降幅度分别为5.37%、2.40%、1.75%和1.04%, 表明城市绿地空间对雨水的调节能力有所下降, 一定程度上加剧了城市雨洪产流风险。其中, 2000年的城市地表径流调节率明显高于其他年份, 主要由于该时段内城市绿地空间比例较高, 这也表明城市绿地空间的合理配置一定程度上能够提升雨洪调节服务。

从空间上来看, 由于老城区和新扩张区内不透水面和绿地空间结构和配置的差异, 使城市地表径流调节率呈“中间低-四周高”的空间分布格局, 并具有明显的空间集聚特征(图 6)。在10a一遇降雨情景下, 在99%的置信区间内1991年、2000年、2009年、2014年径流调节率的Moran′s I指数分别为0.55、0.61、0.53和0.49, 表明街道尺度的地表径流调节率分布具有高度集聚性, 整体呈现以东城和西城为核心的冷点区以及以西部和东南部为核心的热点区。同时, 2014年在95%的置信区间内雨洪调节服务低值聚集区较1991年也有所扩张, 高值聚集区距城市中心更远, 表明城市内部的雨洪产流风险更高(图 6)。

2.3 城市土地利用/覆盖变化对雨洪调节服务的影响分析

由于城市老城区和新扩张区不同时段土地利用/覆盖变化存在差异, 其对城市雨洪调节服务变化的影响程度也有所不同。总的来看, 1991—2014年城市新扩张区土地利用/覆盖变化对雨洪调节服务的变化起主导作用, 贡献度超过84%, 老城区的贡献度不足16%(图 7)。这主要是由于老城区内部土地覆盖结构相对稳定, 变化幅度较小；而新扩张区由于初始时段的自然植被比例相对较高, 变化程度更剧烈。结果也表明不同的降雨强度情景下老城区和新扩张区土地利用/覆盖变化对雨洪调节服务的影响并无明显差异。

同时, 由于城市新扩张区土地利用/覆盖变化的速度和强度逐渐下降, 对雨洪调节服务的影响程度随之逐渐下降, 而老城区对雨洪调节服务的影响程度逐渐增强。从不同时段来看, 1991—2000年由于城市扩张速度较快, 新扩张区对雨洪调节服务变化的贡献度超过69%；2000—2009年, 在城市“外延式”扩张的同时老城区内不透水面也呈“填充式”增长, 导致其不透水面比例不断上升, 使得老城区对雨洪调节服务的影响有所增强, 但新扩张区的贡献率仍然高于53%；2009—2014年, 城市扩张速度明显放缓, 建成区内不透水面以“填充式”模式发展, 老城区对雨洪调节服务的影响逐渐高于新扩张区(图 7)。

3 结论与讨论

本文以北京为研究区域, 基于多期Landsat遥感影像提取了城市土地利用/覆盖数据, 量化了不同降雨情景下的雨洪调节服务, 并定量评估了城市土地利用/覆盖变化对建成区内雨洪调节服务变化的影响。研究主要结论如下：1991—2014年北京城市用地呈快速扩张态势, 城市内部不透水面比例明显上升, 导致城市绿地空间比例下降幅度达到9.90%。地表径流调节率对土地利用/覆盖变化响应剧烈, 城市不透水面的快速增长导致城市雨洪调节服务呈现退化态势。城市新扩张区较老城区相比, 雨洪调节能力更高, 也主导着城市建成区内雨洪调节服务的变化, 但老城区的贡献逐渐增加。在未来城市建设规划中应合理布局绿地空间以提升城市雨洪调节服务供给能力, 减小城市雨洪和内涝风险。

研究发现, 城市地表径流调节能力随降雨强度的增大逐渐下降, 但下降趋势有所减缓, 这可能是由于城市绿地植被和地表对降雨的滞留和下渗能力有限所致, 与已有研究中不同功能区的地表径流调节率随降雨强度下降结论基本一致^{[17, 30]}。同时, 1991—2014年城市建成区内地表径流调节率下降的幅度介于1.04%—5.37%, 但在2000—2010年的下降幅度要高于4.33%, 较已有研究中北京六环内土地利用的变化导致地表径流调节率下降6%略低^[6]。这可能是由于北京六环的范围远大于城市建成区, 土地利用/覆盖变化程度更加剧烈, 其雨洪调节服务退化程度更高。研究结果也表明城市化过程中新扩张区的地表覆盖变化是主导着城市雨洪调节服务的变化。因此, 优化配置城市新扩张区的地表覆盖结构和绿地空间将有助于提升城市整体的雨洪调节能力。对于老城区而言, 其不透水面比例较高, 地表覆盖类型改造难度较大, 可以通过地下排水管网改造以提升城市排水能力或建设适宜性较高的绿色基础措施以提升其雨洪调节服务。此外, 城市地表产流虽主要受土地覆盖比例影响, 但也与景观空间结构和配置有一定关系^{[6, 37]}。本研究未考虑地表覆盖结构配置的影响, 在后续研究中可应用高分辨率遥感影像提取更加精细的城市土地覆盖数据以开展城市雨洪调节服务评估。