1979-2020年深圳市不透水面信息提取及其时空格局变化分析

杨丽媛, 陈洋波

(中山大学 地理科学与规划学院, 广东 广州 510275)

不透水表面是指能够自然或人为隔离地表水渗入土壤的表面。不透水面会改变洪水径流和物质沉淀，它是评价城市进程和生态环境的重要指标[1]。研究表明，不透水表面的不断增长对水循环有诸多影响，地表的不透水性导致降雨难以渗入地下，从而降低了雨水的下渗率，下渗减少而地表径流增多，城市被大量不透水面覆盖，地表径流的流速和流量都大大增加，无疑加大了城市排水的压力和洪涝的风险[2]。不透水面作为研究城市的一个重要指标，其面积大小、增长变化趋势、所占总研究区域的比例等，可以对城市的发展变化适宜程度、人与自然和谐共处程度进行评估，对探究现阶段城市发展的问题、对城市未来发展的规划及海绵城市的建设具有重要的意义。

不透水面的提取方法是近年来的研究热点，它在城市规划设计、智慧城市建设、城市水文监测和水环境管理等方面都有所应用。近20 a中城市不透水面提取技术有了明显进步，针对各种情况探究出了多种方法[3]，大批量的卫星遥感影像解译主要采用遥感影像自动分类算法[4]，如决策树法[5]、人工神经网络法[6]、支持向量机法(SVM)[7]、AP聚类算法[8]等。以上分类方法虽然有着解译效率高、可以获取大范围内的数据的优点，但是在对像元进行分类时，仅将像元的地表覆盖类型分成一种类型，仅能判断一个像元是否为不透水面，并没有对像元的不透水率进行表示。近年来国内外开始研究计算像元不透水率的方法，Ridd提出了著名的V-I-S模型(vegetation-imperviousness-soil)的亚像元分解方法[9]；Wu等[10]基于Landsat ETM+对地表覆被的分类有了进一步的划分，分为了高反照率地物、低反照率地物、植被和土壤，实现了可操作的混合像元光谱分解。在LSMA的基础上，针对固定端元的问题，后续发展了多端元光谱混合分析法(MESMA)[11]。除此之外中外学者也从光谱波段入手，探寻一种直接的提取方法。Xu基于不透水面地类在热红外、中红外、近红外和可见光波段的反射率特征，构建了不透水面指数NDISI[12]。Deng等[13]基于不透水面、植被和土壤在缨帽变换分量上的差异，提出了BCI指数(biophysical composition index)，在实际应用中不透水面与土壤的区分度欠佳，同时该方法需要预先掩膜水体。然而，城市不透水面的精确提取仍然面临多种挑战，例如在实际应用中纯净像元较少并较难提取，单个像元往往包含多种地物类型成为混合像元，混合像元的增多会导致提取地物类型的复杂度上升，如何更精准高效的提取不透水面的问题亟需解决。

本文以深圳市作为研究区，基于30 m空间分辨率的Landsat影像采用完全约束最小二乘算法的混合像元分解方法提取不透水面，解译得到深圳市1979—2020年的不同不透水面比例数据集，并应用不透水面数据计算其平均中心、标准差椭圆，将不同比例的不透水面进行景观格局分析。并根据深圳市不透水面的空间格局变化，总结41 a来城市发展的规律并对深圳市水土保持与环境保护和绿色可持续发展提供参考意见。

1 研究区概况与数据来源

1.1 研究区概况

深圳，简称“深”，是广东省副省级市，位于东经113°43′至114°38′，北纬22°24′至22°52′之间，占地面积1 997.47 km2，地势东南高，西北低，大部分为低丘陵地，是粤港澳大湾区中心城市之一。1980年被国务院确定为经济特区、经过快速的城市化进程发展为全国性经济中心城市和国际化城市。

1.2 数据源及预处理

本文选用Landsat影像为主要数据源。Landsat是20世纪70年代早期唯一可用的卫星图像[14]，此前的历史遥感资料大多为航空影像和侦察卫星。深圳土地利用/覆被遥感监测选用1979—2020年的逐年高质量无云Landsat系列影像作为分类提取底图数据。深圳处于常年多云区，因此尽量选择冬季影像以避免云量和植被生长对土地利用/覆被分类精度的影响，并且2景影像选择相近成像时间。收集到的各期影像都已经进行了几何配准和大气校正等预处理，精度控制在2～3个像元以内，基本满足本研究的精度要求。

2 研究方法

2.1 完全约束最小二乘算法的线性光谱混合分解方法

线性光谱混合分解方法适合中分辨率的遥感影像，但是由于无约束最小二乘混合像元分解算法即线性光谱混合模型(linear spectral unmixing)只有一个约束条件，即解译所得不同地表覆被的盖度图结果的和设为1，所以得到的不透水面盖度图经常会出现负值的情况，很显然这是不合理的。针对此算法提取地物端元丰度出现的局限性，通过完全约束的最小二乘混合像元分解算法(FCLS spectral unmixing)完成对研究区典型地物端元丰度的提取。该方法较好地解决了提取的端元丰度信息出现负值的情况，直接获得每个像元内各种端元的盖度图，DN值在0～1之间，很好的表示了单元间不透水面的程度，并且和为1。公式如下：

(1)

式中：Riλ为第i个像元在波段λ上的反射率;fki为第k个端元类型在第i个像元中所占的比例;Ckλ为第k个端元类型在波段λ上的反射率;ε为残差项。

2.2 精度评估

精度验证使用Google高分辨率影像作为参考数据。在Arcmap中随机生成100个以1 km为直径的圆，为样本单元，分别统计解译的每个样本的不透水面比例;将样本进行矢量化并加载到谷歌高分辨率影像上，采用目视解译的方法人工提取不透水面并计算不透水面比例。通过均方根误差(RMSE)来计算提取精度，计算公式如下：

(2)

式中：Yi为在高分辨率影像上提取的不透水面比例;Xi为混合像元分解提取的不透水面比例,n为样本数。RMSE可以计算混合像元分解不透水面比例和原始影像不透水面比例之间的差异,并且随着两者的差异减小RMSE值下降,所以RMSE值越小意味着混合像元分解的精度越高[15]。

2.3 平均中心和标准差椭圆

不透水面平均中心的计算可以揭示城市迁移和扩张的方向。计算公式为:

(3)

(4)

式中：x和y分别为不透水面平均中心的经纬度;X和Y为算术平均中心;n为不透水面斑块的个数。计算出每年的深圳市不透水面平均中心,分析深圳市不透水面整体变化方向。

本文采用标准差椭圆(standard deviation ellipse， SDE)衡量城市扩张的方向和时空发展趋势[16]。其方位角计算公式为:

(5)

(6)

(7)

(8)

(9)

(10)

方位角可以表示不同时期不透水面的变化方向，长轴与短轴之间的比例表示不透水面的聚集或破碎程度。比值大于1时椭圆扁率较大，此时的不透水面分布具有明显的方向性，分布较为分散;等于1时椭圆扁率小，此时的不透水面分布几乎没有方向性，分布较为聚集。

2.4 景观格局指数

本文选择景观格局指数来定量分析不同密度等级的不透水面的变化特征[17]。其中蔓延度指数(CONTAG)和香农多样性指数(SHDI)可以用来探究景观的扩散程度和异质性，香农多样性指数(SHEI)表示各斑块组成类型均等分布的程度即斑块多样性程度。同时应用斑块密度(PD)和景观形状指数(LSI)来描述不同密度水平下不透水地表的空间破碎度和形状复杂度，应用斑块类型所占景观面积的比例(PLAND)来确定不同时期占优势的不透水面类型，应用斑块聚集指数(AI)来描述不同密度不透水表面的聚集程度。

3 结果与分析

3.1 不透水面时空分布提取结果及其变化分析

以1979—2020年Landsat遥感影像为数据，基于数据预处理，最小噪声分离变换，纯净像元提取及端元选择的研究方法，根据完全约束的最小二乘混合像元分解算法对深圳市长时间序列不透水面进行了提取，随后利用高分辨率Google影像对提取得到的不透水面进行精度评估，精度验证结果表明，估算结果均方根误差(root mean squar eerror，RMSE)为8.94%～13.26%，表明混合像元分解解译具有较好的精度，能够较准确地反映不透水面的空间分布。图1是提取结果及其时空分布。

1979—2020年深圳市不透水面比例快速增长，表1和表2是研究区1979—2020年不透水面面积、面积占比、变化速率Ⅴ及变化强度AGR的计算结果。由表1—2和图1—2可知研究区近41 a不透水面的面积不断增加，占总研究区的比例逐渐增长，但是其增长速度不是均匀的。在1990—2010年间增长速率和增长强度分别在2005年和1995年达到最高，后逐渐变缓，但还是保持着增长的趋势。最终不透水面面积由1979年的18.64 km2增长至2020年的949 km2，不透水面占整个深圳市的比例也由0.94%增至47.51%。

图1 1979-2020年深圳市不透水面信息提取结果

表1 研究区不透水面面积及其变化统计情况

图2 深圳市不透水面面积变化和不透水面比例变化

表2 深圳市各区不透水面变化速率和变化强度

根据不透水面扩张强度和变化速率的公式计算得到相应的数据(表1)，按陈龙乾等[18]研究将年均扩张变化速率划分为：高速扩张(>10 km2/a)，快速扩张(6～10 km2/a)，中速扩张(2～6 km2/a)，低速扩张(<2 km2/a)。由表1可以看出，研究区在1979—1990年期间扩张速率较低，为中速、快速扩张，在1990—2020年期间均为高速扩张，其中1990—2010年期间增长速度最快，到了2010—2020年增长速度放缓。同时由表1可以看出，扩张强度指数前中后时间段强度差异较大，前期扩张强度指数逐渐升高，至1995年达到最高后逐渐下降，并整体呈现逐渐降低的趋势。整体来看扩张强度较后期扩张强度大，前后相差3倍有余，这与不透水面扩张速率的变化相符。

研究区内各区的增长有所不同，为了方便研究，根据研究区的变化速率，将8个时间段划分为4个时间段，分别计算在这4个时间段内各个区的城镇用地扩张相关指数，从表2可以看出各区城镇用地的扩张速率、程度等变化不尽相同，分异明显。整体来看各区扩张变化也存在明显差距：经济发展较早的地区如福田区、罗湖区，整体的扩张强度和变化速率前后差距较小，盐田区和南山区因为地形限制则前段时间为快速扩张，很快速度转为慢速，其他区在前期也是快速的扩张，2000年以后还保持着较大的增长趋势，2010年以后速度下降。

由上分析可以发现，不透水面的变化速率与扩张强度变化总趋势是一致的，均呈现逐渐加快后又放缓的趋势，表明了研究区不透水面面积、变化速率和扩张程度呈明显的正相关关系，即面积增长快的同时变化速率和扩张强度相应较高。通过分析深圳市的城市发展影响因素也不难发现这一变化规律：在1979—1986年期间，由于经济特区刚刚建立，经济条件落后的深圳还没有迅速崛起的能力，土地变化不是短时间内可以剧烈变化的，是一个逐渐增长的过程，所以扩张刚刚起步速率较慢；随着国家政府的政策支持下开始了招商引资、工业区大量建设、房地产投资的快速增长期，大大增加了土地的开发与应用，1986年后城镇化进程加快，直至1990—1995年变化速率达到最大，扩张强度也较前期有了巨大的增长；2000—2010年，深圳市制定和修订了一系列优惠政策来吸引外国投资，第二、三产业的增长率大大提高，该阶段的不透水面变化速率达到了最高。而随着后期适宜开发的土地基本开发完毕，基础设施也逐渐完备，盲目的扩张城市用地导致了土地供需矛盾的升级以及环境的恶化，在2010—2020年时期，深圳新一轮的城市总体规划颁布，放慢城市发展速度，重视土地资源的浪费问题，使得此后的扩张强度在逐渐下降，变化速度放缓。

3.2 城市迁移、扩张方向与空间变化分析

本文利用9期不透水面比例数据计算不透水面平均中心(图3)，探究城市迁移和扩张的方向。1979—2000年，深圳市的不透水面变化具有明显方向性，持续向着东北方向移动，原因是在此期间光明、龙岗、龙华、坪山区以及宝安区的北部城镇化迅速发展，其中坪山区的发展强度非常高，使得不透水面重心不断向北偏东偏移。2000—2010年，西部沿海地区填海造陆扩大城市面积，同时政策“三轴两带”沟通了深圳市东西地区的联系，全市城镇用地重心向西方向迁移。2010年后，随着大鹏新区的建立，惠深沿海高速连接西东，大大促进了深圳北部的发展，发展速度快于中部和西部，使得不透水面重心在此时期向东偏移。但是由于东部的地形原因，城市扩张到一定程度后就变化很小，而西部已经发展成熟的宝安、南山区等填海造陆扩张了城市面积，并且在原有基础上发展为更为高密度的不透水面，所以在2015年后不透水面重心又向西偏移。

图3 深圳市不透水面平均中心点和标准差椭圆

本文利用9期的不透水面数据计算标准差椭圆，分析城市的空间变化。1979—1990年，深圳市的不透水面发展方向变化明显，由表3和图3可知，标准差椭圆的方位角在1979—1990年期间由116.86°下降到83.26°，1990年之后开始上升，到2015年上升至95.74°，表明1979—1990年深圳市不透水面长轴先由西北—东南指向向东—西指向移动，1990年之后不透水面长轴主要沿着东—西方向延伸。

表3 不同时期标准差椭圆参数

整体来看东—西剖面是深圳市不透水面发展的主要方向。1979—2000年，标准差椭圆长短轴差距逐渐变大，长轴从23 237.47 m增长到31 608.38 m，短轴从11 543.35 m增长到12 601.71 m，比率变化为总体上升，椭圆扁率逐渐增大，表明不透水面分布的方向性逐渐明显，且分布由小面积的聚集变为大范围的分散。2000—2020年，除2010—2015年深圳东部大鹏新区的发展导致标准差椭圆扩大外，总体呈现标准偏差椭圆的长轴持续缩小，短轴逐渐增大，长轴从31 608.38 m下降至28 912.37 m，短轴从12 601.71 m增长到13 193.58 m，长轴与短轴的比率变化为波动下降，表明深圳各地不透水面发展的聚集趋势慢慢增加。

3.3 不透水面景观格局分析

表4为深圳市1979—2020年期间的不透水面景观格局指数。由表4可以看出，深圳市的不透水面斑块个数(NP)呈不断增加的趋势，在1979—1990年期间变化较慢，在1990—2010年期间变化迅速。整体不透水面斑块密度(PD)由1979—1990年期间快速增长后逐渐下降，说明前期发展零散而迅速，不透水面景观破碎程度高，后期随着斑块数量不断增长，城市发展不断扩大，斑块相互连接，不透水面景观破碎程度逐渐降低，团聚程度提高。同时斑块聚集指数(AI)显示除1979年团聚程度较高外，此后聚合度呈不断上升趋势最终聚合度达到了79.29，相比一开始的86.61还是有所降低的。1979年的不透水面范围小而聚集，聚合度高，随着全市范围的不断发展，扩大到现在的面积，聚合度较高，结构连通性较好。香农多样性指数(SHDI)增大，说明各拼块类型在景观中呈均衡化趋势分布且景观趋于丰富多元化；香农均匀度指数(SHEI)其值趋近于1时说明景观中没有明显的优势类型且各斑块类型在景观中均匀分布，由表4可以看出1979—2015年香农均匀度指数逐渐上升趋近1，各斑块类型在景观中的均衡化趋势增加。由表4及图4可知，深圳市不透水面景观的蔓延度CONTAG指数总体呈现波动下降，不透水面SHDI与SHEI总体呈现波动上升，说明各等级的不透水面景观整体仍呈增长态势，不同不透水面类型的多样性呈增加趋势，异质化程度增加，逐渐趋于均衡化。

表4 1979-2020年深圳市不透水面景观格局指数

图4 1979-2020年深圳市不透水面CONTAG指数与SHDI指数变化

组分空间构型分析(图5)表明，1979年不透水面的AI指数较高，说明其聚集程度高，但是同时PD和LSI指数可以看出不透水面密度与形状复杂度较低，此时不透水面少而聚集。低、中低密度不透水面景观在1979—1986年有了较大幅度的增长，特别是低密度不透水面PD指数突然增大，斑块密度达到最大，同时AI显示1979—1986年聚集度剧烈下降，说明此时低、中低密度不透水面成为密度等级中最为破碎的景观，之后PD开始逐渐下降，AI波动上升，表明不透水面景观不断发展，聚集程度也在不断提高。中高、高密度不透水面景观的PD保持稳定，小幅度增加，破碎化程度变化小；除高密度不透水面景观外其他等级不透水面的LSI保持增长的趋势，表明这4种密度不透水面景观形状的复杂性逐渐增长，但是高不透水面景观的LSI以2005年为转折点前期增长后期下降，表明高密度不透水面景观的形状逐渐变得复杂后又转为简单。PLAND指数显示在1979—1986年中低密度不透水面比较大幅度增长，在不透水面总体中占到51.07%，随后逐渐下降，高密度不透水面占比在波动中逐渐上升，在2020年达到39.75%，其他密度不透水面的占比在小范围内浮动，总体变化并不大。从斑块连接情况来看，高密度不透水面AI保持较高水平且缓慢增长并在2020年达到88.86，其他盖度不透水面的AI随有阶段性的小幅度上升，但总体指数不高，表明高密度不透水面景观聚集度在不断增强，其他密度不透水面景观的聚集度变化不大，只有小幅度的波动，聚集度以及连通性较高密度不透水地表景观弱。

图5 1979-2020年深圳市不透水面景观格局变化

4 结 论

本文运用完全约束最小二乘算法混合像元分解得到深圳市不透水面比例数据集，并应用不透水面数据计算其变化速率与强度，利用平均中心、标准差椭圆景观格局指数的方法进行了深圳市不透水面发展变化分析。

(1) 1979—2020年深圳市不透水面不断扩张，但整体扩张呈逐渐加快又变缓的趋势，1990—2010年增长速率和增长强度最高，后逐渐变缓。

(2) 1979—2000年，深圳市不透水面发展具有明显方向性，持续向着东北方向移动，随后2000—2020年，全市城镇用地重心延东西方向移动，变化较前期较小。不透水面分布由小面积的聚集变为大范围的分散，后随着不透水面密度逐渐增加而聚集度上升。

(3) 1979—2020年，深圳市不同不透水面类型的多样性呈增加趋势，异质化程度增加，逐渐趋于均衡化；低、中低密度不透水面景观逐渐从最为破碎的且占较大比例的景观开始逐渐变得较为团聚，并且面积逐渐不占优势；中、中高密度不透水面景观形状变得越来越复杂，聚集度不高；高密度不透水面景观面积不断增长，连接性在逐渐增强。

研究结果表明，高密度不透水面聚集度逐渐增长，但高聚集度不透水面与地表径流有显著的相关性[19]，中高密度区和高密度区不透水面对城市热环境影响最大[20]，故未来在城市发展中要。本文从不透水面比例的时空变化分析了深圳市的发展变化规律，面向实际的应用还存在一些不足，将来需要对不透水面与其他生态环境的关联性方面展开深入研究。