典型喀斯特山地城市不透水面时空格局

母浩江， 崔文刚， 鲁 霞

(贵州师范大学地理与环境科学学院， 贵州贵阳550025)

近年来，随着我国城市化进程不断加快，大量农村人口向城市集聚，在促进城市经济发展的同时也导致城市用地压力日益增大，越来越多的森林、 耕地等自然地表被不透水地表所替代。不透水地表又被称为不透水面(impervious surface，IS)，是指常见的道路、 广场、 建筑等水不能通过其下渗到土壤中的城市人工景观[1-2]。IS的大量增加会改变地表下垫面性质[3]，使地表径流和地下水的活动受到干扰，致使城市内涝、 地下暗河断流等现象产生[4-5]，进而给生态环境带来负面影响。及时、 准确地掌握城市IS的分布情况和时空演变对于城市发展规划及生态环境保护具有重要意义。

早期IS信息的获取主要依靠专业人员进行实地采集， 这样的获取方式成本高， 周期长， 时效性差， 使得IS信息的采集和应用十分困难[6]。 进入20世纪以来， 随着遥感技术的飞速发展， 遥感数据获取变得便利， 国内外大量学者开始将遥感技术引入IS的研究。 随着研究的深入， 多种基于遥感技术定量获取IS信息的方法相继被提出， 较典型的方法可概括为光谱混合分析(spectral mixture analysis，SMA)法[7-8]、 决策树模型[9-10]、 回归模型[11]、 指数法[12-13]和分类法[14]5种， 其中线性光谱混合分析(linear spectral mixture analysis，LSMA)法可以有效解决中低分辨率遥感数据的混合像元问题， 被广泛用于提取IS信息的研究[15]。 Wu等[7]使用LSMA法成功提取了美国大都市区哥伦布的IS信息； 唐菲等[8]基于陆地卫星7(Landsat 7)增强型专题制图仪(enhanced thematic mapper plus, ETM+)影像， 利用LSMA法提取了北京等6个不同城市的IS信息， 效果良好； 李苗等[16]通过LSMA法提取区域多个时间节点的IS信息， 分析IS信息的时空变化规律， 进一步探讨城市化进程的范围和规模。 理论方法的日趋成熟也推动着IS的研究向不同的区域环境与应用领域扩展， 但是通过查阅文献发现， 相关研究主要集中于经济较为发达的沿海城市[17]和平原地区[18]， 对于城市化进程起步稍晚且地形条件更复杂的喀斯特地区， IS的研究极少。

贵阳市是贵州省省会， 同时也是世界上喀斯特地貌发育最典型的地区之一[19]， 境内喀斯特地貌占行政区域面积的73.3%， 地形起伏较大， 景观破碎化严重， 生态环境脆弱且具有明显的空间异质性。 近年来， 随着经济实力和城镇化水平不断提高， 贵阳已成为我国西部陆海贸易的新通道、 西南地区的中心城市之一。 在城市化状况快速发展的同时， IS急速增加所造成的环境问题也引起人们的重视， 然而目前尚未发现关于该地区IS信息及其空间变化的研究。

本文中以贵阳市主城区为研究区，采用基于植被覆盖度提取法和LSMA法分别提取研究区的IS信息，利用混淆矩阵对提取结果进行精度验证，选出最优方法对研究区1996—2017年各时间节点的IS数据进行提取，依据结果定量分析21 a间研究区IS的时空分布和演化规律。

1 研究区概况与数据源

1.1 研究区概况

贵阳市位于我国的西南部(东经106°07′—107°17′， 北纬26°11′—27°22′)，地处长江和珠江两大水系流域的分水岭地段，四面环山，地势西高东低，境内丘陵起伏，属典型的喀斯特地形地貌，气候类型为亚热带季风气候[20]。本文中选取该市主城区为研究区，包括南明区、 云岩区、 观山湖区、 白云区、 花溪区、 乌当区6个行政区，其中云岩区、 南明区大部分区域位于盆地，乌当区、 花溪区、 白云区和观山湖区大部分区域以山地为主。

1.2 数据源与预处理

为了避免由影像数据源不同、 拍摄季节不同、 云雾遮蔽等原因造成的误差， 保证研究数据的可靠性， 本文中选取研究区内无云遮蔽且拍摄月份相近的3期Landsat影像作为数据源， 影像数据包含4景陆地卫星5号(Landsat 5)专题成像仪(TM)影像和2景陆地卫星8号(Landsat 8)陆地成像仪(OLI)影像(见表1)。 影像数据来自中国科学院地理空间数据云平台(http://www.gscloud.cn)。前期需要对影像数据进行预处理， 包括辐射校正、 大气校正、 拼接以及按研究区范围裁剪。

表1 研究区遥感数据

2 研究方法

2.1 IS信息的提取

2.1.1 基于植被覆盖度提取法

基于植被覆盖度提取IS数据的方法最早由Carlson等[21]提出。该方法根据IS与植被覆盖度的负相关关系提取IS信息，计算方法相对简单，并且能够快速获取大范围连续的IS信息，已经成为现阶段应用最多的IS信息提取方法。该方法的表达式[13-14，21]为

I=1-F，

(1)

(2)

(3)

式中：I为IS指数，%；F为植被覆盖度，%；N为归一化植被指数；Nsoil为N在裸土或无植被区的值；Nveg为N在纯植被区的值；NIR为近红外波段的反射率，%；R为红光波段反射率，%。由于影像存在不可避免的噪声，因此在一般情况下，Nveg和Nsoil取一定置信度范围内的最大值与最小值。在参考前人经验[22]并多次试验后，本文中选取累积分数95%～5%作为置信范围，得到N的最大值和最小值。

在利用植被覆盖度提取IS信息前，还需对研究区内的水体信息进行掩膜。徐涵秋[23]提出的改进型归一化差异水体指数(modified normalized difference water index，MNDWI)能够清楚地分出阴影和水体，是现阶段提取效果最好的水体提取法，因此本文中使用该指数对水体信息进行提取、 掩膜。该指数的表达式为

(4)

式中：M为归一化水体指数；G为绿光波段的反射率，%。

2.1.2 LSMA法

LSMA法通过计算混合像元中各端元的组合比例，将混合像元的光谱分解为多种端元光谱的线性组合[8]。该方法假设一个像元的反射率由影像拍摄瞬间各地物的纯净反射率(端元)按一定比例组合而成，即

(5)

式中：Rx为波段x的光谱反射率，%；n为端元个数；fi为端元i在像元中所占的比例，%；Ri,x为端元i在波段x的反射率，%；ex为残差。

理论上， 端元的数量在一定程度上决定了LSMA法的拟合程度， 但这并不意味着端元种类越多，模型效果越优。Wu等[7]通过实际研究发现，选取4种端元时IS信息的反演效果最优。本文中依据前人经验[7]对影像进行端元采集时选择由玻璃、混凝土为代表的高反射率地物，传统的旧式瓦片屋顶、沥青路面为代表的低反射率地物，以及植被和土壤4种端元类型。除此之外，研究[2]发现，水体在LSMA法中的表现为低反射率。为了避免水体对研究结果的影响，本文中在提取纯净端元之前使用式(4)中的修正归一化水体指数，对研究区内的水体进行剔除，将最终得到的高反射率盖度图和低反射率盖度图相加得到研究区IS盖度图[24]。

2.2 精度验证

对影像分类提取结果验证的最好方式是随机提取采样点进行实地验证；但是，实地验证成本较高，对过往年份展开实地验证也不切实际，因此，本文中先利用地理信息系统(GIS)工具随机生成验证点，参照验证点的空间位置对基准影像(2007、 2017年基准影像采用同期分辨率为1.07 m的Google影像，1996年的Google影像因效果不佳而采用原始影像)进行目视解译，将解译结果与二值化后的IS结果进行对比，对比结果通过混淆矩阵来反映提取精度。

2.3 IS的空间分布与演变

2.3.1 标准差椭圆分析法

为了探究研究区IS的空间变化情况，尝试利用标准差椭圆分析法进行分析。通过分析标准差椭圆，可以从椭圆中心、 分布范围、 方向和形状等多个角度分析要素的空间分布特征及空间演变过程[25]。在地学领域，主要用椭圆的长轴半径长度与扁率反映空间格局总体要素的集中密度，用长轴偏角反映分布格局的主导方向。

2.3.2 IS的扩展指标

IS扩展指标能够全面反映研究区IS扩展的情况。 本文中选择IS密度、 IS扩展速度这2个指标来反映各区域IS扩展趋势[26]。 IS密度是指某区域中IS的面积占该区域总面积的比例； IS扩展速度是指在一段特定时间段内， 以年为单位的IS面积增长情况。

IS密度计算公式为

(6)

式中：ρa为区域a的IS密度，%；Aa为区域a的IS的面积， km2；Aa,tot为区域a的总面积， km2。

IS扩展速度计算公式为

(7)

式中：v为区域a的IS扩展速度， km2/a；t为时间间隔，a；Apri和Aend分别为区域a研究初期和末期的IS面积， km2。

3 结果与分析

3.1 精度验证结果分析

将2种IS信息提取方法获取的研究区2017年的IS影像选择合适阈值进行二值化处理。利用生成的验证点，分别将2幅二值化IS影像与基准影像叠在一起进行验证，结果见图1。基于植被覆盖度提取法提取的IS总体分类精度为0.89，Kappa系数为0.53，而LSMA法得到的总体分类精度为0.95，Kappa系数为0.74。此外，在目视判断过程中发现基于植被覆盖度提取IS信息的方法很容易将大片裸地误分为IS，而用LSMA法时该现象很少出现。由此可知，采用LSMA法在喀斯特地区具有较好的可行性，并且提取效果优于基于植被覆盖度提取法。

在后续研究中选择LSMA法继续提取研究区1996、 2007年的IS信息， 并进行精度验证， 结果如表2所示。 由表可知， 本文中提取的IS信息均能满足研究要求， 可用于分析研究区IS空间格局演变。

表2 不同方法对贵阳市主城区不透水面的提取精度验证

3.2 IS时空分布

利用GIS技术对得到的1996、 2007、 2017年的IS二值化图(见图2)进行面积统计， 由图2可知， 研究区1996、 2007、 2017年的IS面积分别为50.82、 126.56、 234.08 km2， 21 a间贵阳市IS面积共增加183.26 km2，增长了近5倍，整体呈现出高速扩张趋势，其中1996—2007年IS面积增加了75.76 km2，年平均增长量为6.89 km2/a；2008—2017年IS面积增加量达107.49 km2，年平均增长量为10.75 km2/a，前10 a的平均增长量明显高于后10 a的。在2008—2017年中， 贵阳市内、 外部交通快速发展， 内部以环城高速和中环为代表的各等级道路竣工通车， 外部贵遵附线、 花安高速以及贵广、 沪昆高速铁路建成通车， 大力推动了贵阳市城市的发展。 除此之外， 观山湖区作为缓解旧城区人口稠密、 交通拥堵的发展新区， 是城市扩展的重点方向， 在大量文化商贸区、 居住区、 公共基础设施建设的同时， 也导致了IS的增长。

将提取的IS信息按行政区进行分区统计， 结果见表3。 由表可看出， 1996—2007年IS面积增长量最大的是花溪区， 增长了25.72 km2，增长量最小的是云岩区，仅为8.59 km2。通过观察各区IS面积增长率发现，观山湖区的增长率最高，达到了610.67%， 其次是花溪区、 乌当区， 增长率分别为377.13%、 369.85%，而南明区的增长率最低，仅为45.91%。2008—2017年IS面积增长最多的是观山湖区，其次是花溪区，增长最少的是南明区。对比1996—2007年的数据，各区IS面积增长率明显降低且区域间降低幅度差异较大，其中观山湖区增长率下降超400%，但其增长率仍居各行政区值之首，南明区21 a间IS面积增长率都最低，但其一直处于平稳增长状态。

表3 1996—2007年贵阳市主城区各区域不透水面面积的变化情况

从贵阳市各行政区IS面积的总体变化趋势可以看出： 1996—2007年观山湖区、 花溪区、 乌当区处于城市化扩张爆发期， 大量的IS井喷式出现； 南明区处于城市稳定发展期， IS格局基本稳定， 增长相对缓慢； 云岩区、 白云区则属于典型的中速发展期， IS平稳扩展。 2008—2017年观山湖区、 乌当区仍处于高速发展阶段， 强度明显减弱； 花溪区与云岩区、 白云区一样进入城市中速发展期， IS继续平稳扩张。

为了更直观地展现各阶段研究区IS面积变化的空间分布情况， 利用GIS工具计算研究区不同阶段之间IS面积差值， 由此得到不同阶段IS面积变化的图斑， 结果见图3。 由图可以看出， 1996年以前研究区IS分布较集中， 主要分布在以南明区、 云岩区为主的市中心； 1996—2007年IS面积变化主要分布在各行政区中心周围， 同时城市中心的IS面积增加也很显著； 2008—2017年， IS的增加使得之前许多破碎的斑块连接起来， IS开始出现成片聚集现象， 还有很多新增IS以“飞地”形式出现在各行政区中。

3.3 IS格局演变分析

3.3.1 IS空间变化特征

通过分析研究区1996、 2007、 2017年IS的标准差椭圆，探讨该区域IS的整体变化趋势， 如图4所示。由图可以看出：1996—2007年，标准差椭圆面积明显增大且整体向南移动，表明此阶段IS扩展重心正在向南倾斜；主轴方向偏角从9.87°减小到2.30°， 意味着该区域IS的分布趋势由东北—西南方向往正北—正南方向发展； 标准差椭圆的长半轴由19.11 km增大到22.58 km， 短半轴由10.46 km增大到13.51 km，表明IS的空间分布在主轴方向(南北方向)和辅轴方向(东西方向)都呈现较弱的离散化现象。

图3 贵阳市主城区不同时期不透水面的变化

2008—2017年研究区IS标准差椭圆范围逐渐缩小且向东北方向移动,扩展重心逐渐向西北方向移动； 主轴方向的偏角从2007年的2.30°增至17.02°， 该阶段贵阳市主城区IS的分布方向转变为东北—西南方向， 与1996年相同且趋势更加显著； 长半轴由22.58 km减至22.2 km， 短半轴由13.51 km增至14.67 km， 可知IS的空间分布在主轴方向(南北方向)的离散化趋势有所缓解，在辅轴方向(东西方向)离散化程度进一步加强。

(a)1996年 (b)1996—2007年 (c)1996—2017年

纵观研究区IS的变化情况， 与研究区地形密切相关。 贵阳市属于典型的喀斯特山地城市， 地表景观破碎化严重。 南明区、 云岩区是贵阳市的老城区， 大部分区域位于贵阳—中曹司盆地内， 四面环山， 受地形严重限制， 迫使其城市范围不断寻求向外扩展， 进而致使IS向其他行政区扩张， 但内部填充式扩张仍是其扩展的重要方式， 这种状况一直保持到2000年。 2000年后， 大量的交通干线开始建设， 以垭口和河谷为走廊逐步向地势平坦的花溪区、 观山湖区延伸[27]。 2001—2006年间，花溪区成为贵阳市交通的主要扩展地区， 交通建设量占同期工程建设总量的45%。 从2006 年开始， 贵阳市城市建设发展的主方向转移至观山湖区， 其次是花溪区， 这2个方向扩展的比例达到同期工程建设总量的80%[27]， 这些变化情况都与研究区IS标准差椭圆的变化相呼应。

3.3.2 IS扩展趋势

1996—2017年贵阳市主城区各区域IS的扩展指标如图5所示。 从图5(a)中可以看出， 1996、 2007、 2017年研究区的IS密度分别为2.14%、 5.32%、 9.84%， 22 a间增加了7.7%。各行政区的IS密度变化存在明显差异，其中云岩区的IS密度在2个研究阶段都是最高的，1996、 2007、 2017年的IS密度分别为12.23%、 21.62%、 36.33%，而乌当区在整个研究阶段IS密度最低， 增长最缓慢。从图5(b)中可以看出：1996—2007年研究区IS的扩展速度为6.89 km2/a，各行政区IS扩展极不均衡，其中花溪区的扩展最快，扩展速度为2.34 km2/a，是各行政区扩展速度平均值的2倍，乌当区的IS扩展最慢，扩展速度仅为0.67 km2/a；2008—2017年，研究区IS扩展速度整体上升至10.75 km2/a，花溪区虽然保持较高的扩展速度，但是扩展速度呈倒退趋势，其余各区扩展速度明显提升，特别是观山湖区扩展速度上升最多，达到了2.65 km2/a，扩展速度最低的南明区也较前一阶段提升了0.4 km2/a。

(a)IS密度

4 结论

1)通过对LSMA法和基于植被覆盖度提取法得到的贵阳市1996—2017年IS信息进行精度验证可知，LSMA法的提取结果明显优于基于植被覆盖度提取法的，说明在喀斯特地貌环境下使用LSMA法提取IS信息是适用的，并且效果更优。

2)1996—2017年贵阳市IS面积从50.82 km2增加到210.12 km2，该市的IS扩展可分为3个阶段：前期(1996年)IS分布比较集中，主要分布在城市中心区域；中期(1996—2007年)主要呈现出以各行政区为发散中心向四周扩张，同时伴随着城市中心区域的内部填充；后期(2008—2017年)IS的扩充得到进一步发展，扩展方式主要呈现出多中心发散的特点，并出现IS“飞地”，其中观山湖区扩张尤为突出。

3)研究区IS面积的标准差椭圆偏角从1996年的9.87°增加到2017年的17.02°，表明研究区IS的分布趋势表现为东北—西南方向，并且随着城市的发展，该分布态势得到进一步加强。

4)本文中将LSMA法应用于喀斯特地区IS信息的提取研究，与基于植被覆盖度提取法提取的IS信息相比，LSMA法更适用于喀斯特地区。本文中还提取了研究区多个时间节点的IS信息，对IS信息的时空演变进行了进一步探讨，研究结果为喀斯特地区IS研究和贵阳市后期城市发展规划、 生态文明城市建设提供了参考。对提取到的IS信息开展相关应用研究将是后期研究的重点。