基于Python的Landsat8OLI影像建设用地自动识别与提取

吴亚玲 黄义忠 彭秋志

摘要：为了提高从海量遥感影像中解译建设用地的效率，设计一种基于Python从Landsat8 OLI遥感影像中自动识别与提取建设用地的软件系统。利用时域中值濾波去除云噪声，并通过z值标准化获得可适用统一阈值划分的稳定像元；结合多种地物指数开展阈值划分以初步识别可能的建设用地像元；结合像元邻域组合关系进一步精化提取结果，并自动计算对应的Kappa系数。研究表明，与传统基于大型专业遥感影像解译平台相比，该系统不仅能保持较高的解译精度，而且具有体积小、速度快、自动化程度高、可迁移性和可扩展性强等特点。

关键词：Python；遥感影像；自动识别；建设用地

DOIDOI：10.11907/rjdk.181322

中图分类号：TP319

文献标识码：A 文章编号：1672-7800（2018）010-0161-04

英文摘要Abstract：In order to improve the efficiency of interpreting construction land from massive remote sensing images， the software system based on Python which can automatically identify and extract construction land from Landsat 8 OLI remote sensing images was designed. Firstly， median filtration was used to remove cloud interference， the stable pixels of uniform threshold value partitioning that was obtained through normalization of z value. Then，the threshold partitioning combined with various features indices are employed to identify possible pixels of construction land. Finally， the meta-neighborhood combination relationship is further refined to extract the result， and the corresponding Kappa coefficient is automatically calculated. The results show that the proposed system can not only maintain high accuracy of interpretation， but also has the advantages of small size， high speed and high degree of automation compared with the traditional method of recognition and extraction of construction land based on large professional remote sensing images interpretation system.

英文关键词Key Words：Python；remote sensing image；automatic identification；construction land

0 引言

社会经济的持续发展以及城市化进程中伴随的城市人口、社会发展要素重组等，对城市建设用地利用方式产生了巨大影响[1-2]。城市建设用地动态监测空间基础数据中遥感影像分辨率的提高、分类方法的改进，对实现城市建设用地的高效利用具有重要意义[3]。

随着遥感影像大范围普及，高效、快捷、自动化地获取建设用地信息已经成为一种趋势[3-4]。近年来，为了获得更高的分类精度，众多从遥感影像中提取建设用地信息的方法相继被提出。这些方法就分类单元而言可分为基于像元的分类和基于对象的分类[5]，就分类途径而言可分为监督分类和非监督分类[6]，就分类结果表现形式而言可分为确定性分类和模糊分类[7-8]，就分类所用信息类别而言可分为仅基于光谱信息的分类和以光谱信息为主的多信息组合分类等[9-10]。然而，日益多样的方法虽然带来更多选择可能，但同时也增加了学习成本、增大了选择困难。在实际工程应用中，不仅需要考虑解译精度，而且需要兼顾简单易学、便捷高效、成本低廉等要求。因此，有必要针对特定的专题解译需求，开发出一系列自动化软件系统，以减轻对大型专业遥感影像处理平台的依赖，大幅降低成本投入。

本文借助面向对象的解释型计算机程序设计语言——Python，综合运用多波段指数法[11]和归一化阈值控制法[12 -14]，提出一种节约时间、更新速度快的建设用地识别与提取技术，以期自动化和批量化从海量遥感影像中解译建设用地。

1 系统原理与流程设计

建设用地自动识别与提取系统流程如图1所示。基本原理是，利用自动化文件遍历和迭代方法串联起每一个解译步骤，实现“进去的是原始数据，出来的是最终产品”之目的。该系统主要包括6个关键解译步骤：①对下载的原始遥感影像进行数据预处理，以便实现数据输入自动化，主要模块包括解压缩文件、重命名文件、提取关键波段文件和获取影像最小重叠区；②叠合同一时期多期影像，利用时域中值滤波法实现去云降噪；③通过z值标准化获得可适用统一阈值划分的稳定像元；④结合多种已被证明效果良好且稳定的地物指数开展阈值划分，以初步识别可能的建设用地像元；⑤结合像元邻域组合关系进一步精化提取结果；⑥自动计算Kappa系数，并根据该系数对前两步进行迭代式调整，保证解译结果的质量。

2 系统实现

系统实现基于Python语言进行编码，数据预处理和质量控制部分通过自行编码实现，主要空间分析函数调用自ArcGIS中的空間分析模块。

2.1 数据预处理

（1）对下载的原始数据进行批量解压。对海量原始影像压缩文件夹进行遍历，识别并逐次解压。文件解压关键代码如下：

targzfile.extract（file， unzipedfolder） #从压缩文件targzfile中把解压出的文件file放到文件夹unzipedfolder里

（2）重命名。在系统中，为了后续流程的完整性，对已解压文件名不统一的波段重新命名。文件重命名关键代码如下：

arcpy.Rename_management（name_fr， name_to） #从将文件name_fr重命名为name_to

（3）提取关键波段。海量影像的所有波段中，依据建设用地识别方法，对波段进行筛选，单独提取出关键的波段数，并进行另存。

（4）获取最小重叠区。去除空值和影像边界，其留下的最小边界作为后面研究范围。以关键波段中某一波段为代表，利用Python脚本对其影像进行批量化掩摸，计算最小重叠区。关键代码如下：

arcpy.gp.CellStatistics_sa（in， out， "MINIMUM"； "DATA"） #获取所有输入图层的最小值

剔除外围以0值表示的空值黑边，只保留有值区，裁剪出最小重叠区。关键代码如下：

arcpy.gp.Con_sa（in ， 1， out， ""， "Value > 0"） #剔除外围黑边

2.2 云噪声去除

为了消除时相干扰，去除云噪声，需进行时域中值滤波。其原理是云的反射率在各波段都位于异常极大值区，而云影大多位于异常极小值区，取中值具有很好的去云和去云影效果。关键代码如下：

arcpy.gp.CellStatistics_sa（in， out， "MEDIAN"， "DATA"） #获取所有输入图层的中值

2.3 稳定像元获取

不同影像中波段值存在时相差异，需对生成的波段中值进行z值标准化，并计算影像的z值标准差，再以0.5倍标准差为阈值，替换不稳定像元，最后反z值标准化，即基于统一的平均值和标准差系数反算回去，获取每个单波段的稳定像元图层。该步骤能保证在后续步骤中影像均采用相同阈值进行分割。关键代码如下：

MEAN = arcpy.GetRasterProperties_management（in， "MEAN"， "Band_1"） #计算多期影像平均值

STD = arcpy.GetRasterProperties_management（in， "STD"， "Band_1"） #计算多期影像标准差

arcpy.gp.Minus_sa（in， MEAN， numerator） #获得每期影像与多期影像平均值之差，以作为标准化的分子numerator

arcpy.gp.Divide_sa（numerator， STD， out） #获得每期影像与numerator之比，得到该期影像的标准化值

arcpy.gp.Con_sa（in ，1， out， ""， "Value < 0.5"） #获得均值两侧0.5倍标准差范围内的像元，即稳定像元

2.4 建设用地初步识别

建设用地提取过程实质是遥感影像信息的处理过程[15] ，遥感技术手段具有客观、宏观、快速、准确、动态等特点[16]。不同类型的地物反射光谱强度在不同波段上具有不一致性，分析和找出建设用地与其它背景地物间的差异，充分利用遥感影像光谱信息、时相信息、地形信息和邻域信息，逐步排除确定性无关地物。运用波段差异提取出其特征信息的NDVI（归一化植被指数）[17]、MNDWI（归一化水体指数）[18]，其值介于（-1，1）之间，经过其处理后比值运算用于消除地形差异的影响。

通过地物指数计算阈值分割方法，在已选用的地物指数参数中，构建各参数的阈值范围，如：NDVI、MNDWI。同时在提取城市建设用地等地物中，仅依据指数去除植被和水体是完全不够的，在系统中根据每两个波段间的归一化差值取阈值去除高亮地物（旱裸土、部分滩涂、部分高亮人工地面、部分高亮屋顶）、水体及阴坡植被、裸地及植被、海岛湖泊及滩涂等其它不属于建设用地的地物类型。其中地物参数选取和阈值分割均在Python中批量化处理，用到的波段组合原理为：

NDVI=NIR-RedNIR+Red（1）

MNDWI=Green-SWIR1Green+SWIR1（2）

D73=SWIR2×Red（3）

ND75=SWIR2-SWIR1SWIR2-SWIR1（4）

ND42=NIR-GreenNIR+Green（5）

ND75ND42=ND75-ND42（6）

式（1）-式（6）均表示在Python中计算用到的指数。针对实验区现状，反复试验后，优先构建出一套完整的建设用地提取流程，用来提取建设用地信息。其中，部分关键代码如下：

（1）根据波段计算出NDVI、MNDWI指数。

（2）采用归一化阈值法确定植被的阈值范围为0.25，粗略去除植被。

arcpy.gp.Con_sa（MNDWI ，1， Pout， 0， 'Value > 0.25'） #获取像元大于0.25的NDVI，用于识别植被

（3）采用步骤（2）中的归一化阈值法，确定水体的阈值范围为0.1，去除水体。

arcpy.gp.Con_sa（MNDWI ，1， Pout， 0， 'Value > 0.10'） #获取像元大于0.1的MNDWI，用于识别水体

（4）在前3步的基础上，通过光谱和纹理分析，去除试验区内比较明显的植被、水体以后，通过光谱差异和借用实验区范围内的Google Earth影像对比，进一步识别出需要识别去除的高亮地物（含裸土、部分滩涂、部分高亮人工地面、部分高亮屋顶）、水体及偏阴坡植被、裸地及植被。

arcpy.gp.Con_sa（D73 ， 1， Pout， 0， 'Value > 0.06'） #获取像元大于0.06的D73，用于识别高亮地物

arcpy.gp.Con_sa（D73 ，1， Pout，0， 'Value < 0.01'） #获取像元小于0.01的D73，用于识别水体及偏阴坡植被

arcpy.gp.Con_sa（ND75ND42 ，1， Pout， 0， 'Value < -0.40'） #获取像元大于-0.04的ND75ND42，识别裸地及植被

2.5 建设用地精化提取

初步识别可能的建设用地后，系统用于精化提取建设用地的方法主要是采用空间邻域关系计算。分别对水体、植被的边界像元通过邻域均值计算后，采用条件函数逐一判别消除其它地物。其关键代码如下：

arcpy.gp.Con_sa（Vm ， 1， Vm01 ， 0， "VALUE > 0.25"） #获取邻域内均值Vm超过0.25的植被

arcpy.gp.Con_sa（Wm ， 1， Wm01 ， 0， "VALUE > 0.25"） ##获取邻域内均值Wm超过0.25的水体

arcpy.gp.CellStatistics_sa（Vm01&Wm01;， Vm01sumWm01， "SUM"， "DATA"） #获取邻域内水体和植被Vm01& Wm01均超过0.25的像元的和

arcpy.gp.Con_sa（Vm01sumWm01 ，1 ， Pout， 0， "VALUE = 2"） #获取上一步中VALUE = 2的值

2.6 质量控制

对遥感分类精度进行评价，需要参考数据作为评价基准。系统选择Google Earth影像作为参考数据，在针对完整流程提取建设用地的基础上，于所选实验区内随机生成1 000个参考点，对建设用地提取结果进行精度评价。采用的评价指标为Kappa系数，其计算在Python中进行。

3 实例验证

3.1 实验区选择与数据下载

实验选取昆明中心城区（五华区、盘龙区、西山区、官渡区、呈贡新区和空港经济区）作为实验区，如图2所示。采用2015年的Landsat 8 OLI遥感影像数据，数据下载自https：//glovis.usgs.gov/。为保障影像质量，避免云量干扰，下载的影像数据云量均控制在10%以内。对地面同一地点，下载同时区影像数量确保大于或等于10，本文实验区所在的影像行列号为129043。以上条件限制主要为了在后续自动提取建设用地过程中，能利用尽量多的时相信息消除云干扰影响，减少季相差异，从而提高分类精度。其Landsat 8 OLI遥感影像数据均经过Level1标准地形校正，即经过系列辐射校正、地面控制点几何校正和DEM地形校正，并以单波段输出格式为GeoTIFF的形式封装在压缩包文件内，取样方式均为三次卷积，地图投影均采用UTM-WGS84投影。本文中运用到的波段为：Band3（Green）、Band4（Red）、Band5（NIR）、Band6（SWIR1）、Band7（SWIR2）5个波段，分辨率均为30m。

3.2 建设用地提取结果

图2显示了研究区Landsat 8 OLI影像的RGB波段合成效果，图3是利用本文所述方法自动提取的建设用地效果。直观对比，本文所述方法能够很好地将建设用地与其它地物区别开来，特别是能很好地剔除耕地和滩涂；就精度定量评估而言，建设用地提取结果的Kappa系数为0.812 5，与大多数现有方法提取建设用地结果的精度水平相当，可以满足大部分相关研究和应用的精度要求。

4 结语

本文在回顾和总结遥感影像建设用地提取研究的基础上，介绍了一种基于Python从Landsat8 OLI遥感影像中自动识别与提取建设用地的技术方法。与其它建设用地提取方法相比，其可以在大范围遥感影像覆盖区域实时、高效、准确地自动提取建设用地信息，为快速监测和评价建设用地变化提供新的途径。在后续研究中，可将本文研究方法延伸到其它地物的获取、监测等领域。

参考文献：

[1] 袁丽丽.城市化进程中城市土地可持续利用研究[D].武汉：华中农业大学， 2005.

[2] 马兵.扬州市城市化进程中城乡建设用地动态变化研究[D].南京：南京农业大学， 2009.

[3] 王宏胜，李永树，吴玺，等.结合空间分析的面向對象无人机影像土地利用分类[J]. 测绘工程， 2018（2）：57-61.

[4] 何少林，徐京华，张帅毅.面向对象的多尺度无人机影像土地利用信息提取[J].国土资源遥感，2013，25（2）：107-112.

[5] 张洋华，刘慧平，杨晓彤.基于像元与对象分类的景观指数差异性分析[J].遥感技术与应用，2016，31（1）：119-125.

[6] 赵春霞，钱乐祥.遥感影像监督分类与非监督分类的比较[J].河南大学学报：自然科学版，2004，34（3）：90-93.

[7] 张卉，张超，赵冬玲.特征权重优化高分辨率遥感影像模糊分类研究[J].遥感信息，2010（1）：94-98.

[8] 刘艳芳，兰泽英，刘洋，等.基于混合熵模型的遥感分类不确定性的多尺度评价方法研究[J].测绘学报，2009，38（1）：82-87.

[9] 陈晨，张友静.基于多尺度纹理和光谱信息的SVM分类研究[J].测绘科学，2009，34（1）：29-31.

[10] 厉小润，朱洁尔，王晶，等.组合核支持向量机高光谱图像分类[J].浙江大学学报：工学版，2013，47（8）：1403-1410.

[11] 樊彦国，李翔宇，张磊，等.基于多波段分析的无阈值自动光谱角制图分类法[J].地理与地理信息科学，2010， 26（2）：38-41.

[12] VIOVY N， ARINO O， BELWARD A S. The best index slope ex-traction（BISE）：a method for reducing noise in NDVI time series[J].International Journal of Remote Sensing，1992，13：1585-1590.

[13] 徐影，丁一汇，赵宗慈.人类活动引起的我国西北地区21世纪温度和降水变化情景分析[J].冰川冻土，2003，25（3）：327-330.

[14] 付忠良.一些新的图像阈值选取方法[J].计算机应用，2000，20（10）：15-17.

[15] 李爱民.基于遥感影像的城市建成区扩张与用地规模研究[D].郑州：解放军信息工程大学，2009.

[16] 邢雅娟，刘东升，王鹏新.遥感信息与作物生长模型的耦合应用研究进展[J].地球科学进展，2009，24（4）：444-451.

[17] DEFRIES R， TOWNSHEND J. NDVI-derived land cover classification at global scales[J]. International Journal of Remote Sensing， 1994， 15（17）：3567-3586.

[18] 徐涵秋.福清市城鎮空间扩展规律及其驱动机制分析[J].遥感技术与应用， 2002， 17（2）：86-92.

（责任编辑：何 丽）