CN112329790A - 一种城市不透水面信息快速提取方法 - Google Patents

Abstract

本发明公开一种城市不透水面(ISA)信息快速提取方法，包括：收集研究区Landsat‑5 TM、Landsat‑8 OLI或Sentinel‑2遥感影像中的任一种，采用简化的COST模型或sen2cor模型对影像进行大气校正，将原始影像灰度值转换为地表反射率值；根据影像类型，选择一组遥感分类特征，计算对应分类特征指数；基于多特征决策树模型，建立单规则分类或多规则分类指数的决策规则，进行ISA信息逐级分类；对提取的城市ISA信息进行精度评价，输出满足精度的城市ISA图。简化的COST模型解决了采用原始COST模型大气校正过程的复杂计算问题，提高了计算效率；提出的城市ISA信息快速提取方法可避免单纯利用不透水面指数单一分类特征的缺陷，极大提高了ISA信息提取精度，解决了城市ISA信息提取的关键技术问题。

Description

一种城市不透水面信息快速提取方法

技术领域

本发明属于图像智能处理技术领域，涉及一种遥感影像信息提取方法，特别涉及一种城市不透水面信息快速提取计算方法。

背景技术

城市不透水面(Impervious surface area，ISA)是一种硬质地表，由水泥、沥青、金属、玻璃等材质构成的建筑物屋顶、道路和停车场等。不透水面覆盖度(Impervioussurface area percentage，ISA％)是指单位面积地表中不透水面的面积的比例。随着城市化的快速推进，城市人口和城市空间规模急剧扩张，导致地表覆被和城市环境发生巨大变化，地表覆被方面主要表现为不透水面逐步取代了植被为主的地表自然景观，随之而来的对生态环境造成明显的负影响，如增强城市热岛、恶化城市水质和增加城市洪涝灾害风险等。不透水面被作为区域和城市环境的指示器被广泛应用，如城市生态水文过程模拟、城市热环境效应分析预测和城市建设规划等研究。因此在快速城市化进程中，如何快速准确提取精度较高的城市不透水面信息，进而分析其空间动态变化，对城市规划、城市生态环境及水文气候研究等方面有重要意义。

利用传统的观测方法对不透水面信息的提取和统计是较为困难的。而具有大范围、多尺度、可重复的对地观测优势的遥感技术的出现，为不透水面信息提取提供了新的技术手段。目前为止，从遥感卫星图像中提取ISA的算法主要有三大类，即机器学习方法(如人工神经网络)、光谱分解技术(如线性光谱混合分析和多端元光谱混合分析法)和光谱指数方法。光谱指数法，通常使用一个简单的线性组合公式来表示，公式中的多个参数可以通过原始波段、反射波段或其他特征指数波段获得。这种基于指数的方法不需要太多图像预处理的先验知识，与光谱分解技术相比，指数法无需复杂繁琐的光谱端元选择的过程；与机器学习方法相比，指数法具有实现简单、计算方便、易于理解与操作等优点。而且，一些优秀的光谱指数法(例如Deng和Wu提出的生物物理组成指数BCI方法)可以基于不同空间尺度的遥感影像中，绘制不透水面信息，其应用范围从低分辨率影像(如MODIS)、中等分辨率(如Landsat 5)到高分辨率影像(如Sentinel-2、IKONOS)。然而，由于遥感影像中，“同物异谱，同谱异物”的现象和受地物空间结构的复杂性影响，使得城市中的裸土和滩涂湿地等地表覆盖类型与不透水面相混淆，它们对不透水面提取的影响较大。因此单独的光谱指数法并不能很好地处理线性不可分的地物，其分类精度不高，往往难以达到分类要求。

大气校正是城市不透水面信息提取遥感信息提取的重要环节和前提，但常用的大气校正方法如大气层顶表观反射率模型TOA和基于影像本身的大气参数估计模型——COST(Chavez，1996)等，反射率计算过程相对比较复杂。因此如何简化原始的大气计算公式，提出快速大气校正的方法，同时综合多种分类方法的优点，探索利用中高空间分辨率影像解决图像的分类特征识别，快速准确提取城市不透水面信息，尚需深入研究。

正是基于以上分析，本案由此产生。

发明内容

本发明的目的，在于提供一种城市不透水面信息快速准确提取方法，其可避免单纯利用不透水面指数单一分类特征的缺陷，具有一定的通用性，极大提高了不透水面信息提取精度，解决了城市不透水面信息提取的关键技术问题。

本发明的另一目的，在于提供一种城市不透水面信息快速准确提取方法，其采用适合于Landsats系列影像的简化的大气校正方法，解决了Landsats系列影像大气校正过程复杂计算的问题。

为了达成上述目的，本发明的解决方案是：

一种城市不透水面信息快速提取方法，包括如下步骤：

步骤1，收集指定研究区的Landsat 5 TM、Landsat 8 OLI或Sentinel-2遥感影像中的任意一种，然后对收集的影像进行大气校正，将原始影像的灰度值转换为地表反射率值；

步骤2，根据遥感影像类型，选择一组遥感分类特征，分别计算对应的分类特征指数；

步骤3，基于遥感分类特征和城市土地覆盖类型，建立单规则分类或多规则分类指数的决策规则，并依据规则进行城市不透水面信息逐级分类；

步骤4，对提取的城市不透水面信息进行精度评价，输出满足精度的城市不透水面图。

上述步骤1中，对于Sentinel-2影像，采用sen2cor模型进行大气校正；对于Landsat 5 TM、Landsat 8 OLI影像，采用如下大气校正方法：

ρ_COST＝M_COST·(DN-DN_min)/cosθ+0.01

M_COST＝M/(TAU_z·TAU_v)

其中，DN为经过辐射校正的图像灰度值，DN_min为遥感器的最小光谱辐射值，θ为太阳天顶角，M是影响各波段的反射率调整参数，TAU_z为从地面到传感器的大气透射率；TAU_v为从太阳到地面的大气透射率；ρ_COST表示Landsats系列影像的地面反射率；M_COST表示大气校正参数。

上述步骤2中，若遥感影像类型为Landsat 5 TM或Landsat 8 OLI影像，选择的遥感分类特征包括不透水面指数BCI，修订的水体指数MNDWI，裸土指数BSI，湿度指数WI，亮度指数BI和绿度指数GI；若遥感影像类型为Sentinel-2遥感影像，选择的遥感分类特征包括不透水面指数BCI，水体指数NDWI，新裸土指数MFI，湿度指数WI，亮度指数BI和绿度指数GI。

上述步骤2中，水体指数MNDWI、修订的水体指数MNDWI的计算公式是：

MNDWI＝(ρ_Green-ρ_SWIR1)/(ρ_Green+ρ_SWIR1)

NDWI＝(ρ_Green-ρ_NIR)/(ρ_Green+ρ_NIR)

式中，ρ_Green、ρ_NIR和ρ_SWIR1分别为经过大气校正的绿光波段、近红外波段和短波红外波段1的反射率。

上述步骤2中，不透水面指数BCI的计算公式是：

BCI＝[(H+L)/2-V]/[(H+L)/2+V]

其中，H、L、V分别为归一化亮度分量、归一化湿度分量、归一化植被绿度分量。

上述步骤2中，H、L、V的计算公式为：

H＝(BI-BI_min)/(BI_max-BI_min)

V＝(GI-GI_min)/(GI_max-GI_min)

L＝(WI-WI_min)/(WI_max-WI_min)

其中，BI、GI、WI分别为亮度指数、绿度指数、湿度指数，BI_max、GI_max、WI_max分别为对应的亮度指数图像、绿度指数图像、湿度指数图像中的最大值，BI_min、GI_min、WI_min分别为对应的亮度指数图像、绿度指数图像、湿度指数图像中的最小值。

上述步骤2中，裸土指数BSI的计算公式是：

BSI＝[(ρ_Red+ρ_SWIR1)-(ρ_NIR+ρ_Blue)]/[(ρ_Red+ρ_SWIR1)+(ρ_NIR+ρ_Blue)]

其中，ρ_Green和ρ_Red分别是绿光波段和红光波段的反射率，ρ_NIR和ρ_SWIR1分别为近红外波段和中红外波段1的反射率。

上述步骤3中，单规则分类或多规则分类指数的决策规则的具体内容是：

基于MNDWI或NDWI提取城市水体，区分水体与非水体；

基于亮度指数、绿度指数和湿度指数，计算不透水面指数；基于不透水面指数提取非植被，区分植被与非植被；

基于BSI或MFI提取裸土，区分裸土与非裸土；

基于湿度指数提取滩涂湿地，区分滩涂湿地与不透水面；

基于湿度指数变化检测，区分漏分的滩涂湿地或误分的高楼的影子。

上述步骤4的具体内容是：所述步骤4中，从地图下载器下载16级以上的GoogleEarth高分辨率影像数据，并进行重采样为1m；采用分层随机采样法为每种土地覆盖类型抽取N个样本点，建立混淆矩阵，计算总体分类精度与Kappa系数，如果分类精度不满足要求，则调整决策树中的关键参数。

上述步骤4后，还包括步骤5：根据步骤4得到的城市不透水面图，对不透水面指数BCI进行规一化，根据下式计算不透水面盖度ISA％：

其中，BCI_max、BCI_min分别为归一化不透水面指数的最大值和最小值；ISA％为不透水面盖度，表征不透水面的聚集密度，范围在0～1之间，ISA％越大，说明不透水面密度越高，0值说明该像元范围内无不透水面下垫面。

采用上述方案后，本发明的有益效果如下：

(1)针对现有遥感影像大气校正程序复杂的问题，本发明提出了适用于Landsat 5TM、Landsat 8 OLI两种遥感影像简化大气校正公式，把两种影像的大气校正模型的方程形式统一起来，大大简化了遥感影像的大气校正计算过程；

(2)本发明公开了一种中高分辨率的遥感数据城市不透水面提取方法，针对现有提取不透水面精度不高的问题，本发明采用中高空间分辨率的Landsat 5 TM、Landsat 8OLI和Sentinel-2遥感影像，构建了多特征决策树模型，相比现有BCI分类法，本发明的普适性更好，多特征模型决策规则参数比较稳定，避免了仅仅采用原始BCI方法过于单一，精度不高的缺陷；

(3)利用MNDWI或者NDWI指数法去除水体，然后利用BCI指数法对遥感影像数据进行植被区和非植被区的区分；

(4)对于不透水面与裸土易混淆的问题，提出使用BCI&BSI或BCI&MFI指数联合提取裸土信息，对二者进行区分；

(5)针对城市滩涂湿地与不透水面从光谱上难以区分的问题：本发明提出首先利用WI湿度指数方法区分出大部分滩涂湿地，然后提出使用两期影像的WI指数变化检测方法进一步区分，该方法能够快速准确区分误分和漏分的不透水面。该方法的能够区分细小滩涂湿地，明显提高了不透水面的分类精度高。

附图说明

图1是本发明的流程图；

图2是本发明中多特征决策树提取不透水面信息的流程图；

图3是实施区域不透水面信息提取结果图；

其中，(a)、(b)和(c)分别对应2009、2018和2019年。

具体实施方式

以下将结合附图，对本发明的技术方案及有益效果进行详细说明。

如图1所示，本发明提供一种城市不透水面信息快速提取方法，包括如下步骤：

步骤1，收集指定研究区的Landsat 5 TM、Landsat 8 OLI和Sentinel-2等遥感影像，分别对Landsat系列影像和Sentinel-2影像进行大气校正，将原始影像的灰度值转换为地表反射率值；

步骤2，在对遥感影像大气校正基础上，根据遥感影像类型和城市土地覆盖类型，选择对应的1组遥感分类特征进行计算，形成遥感分类特征指数影像的空间数据库；

步骤3，建立基于遥感分类特征指数影像数据的多特征决策树模型(见图2)，分别建立单规则分类或多规则分类指数的决策规则，并依据规则进行城市不透水面信息逐级分类，确定决策树的各个分支，包括水体、植被、裸土和不透水面等分支；

步骤4，对提取的城市不透水面信息进行精度评价，输出满足精度的城市不透水面图；

步骤5，基于城市不透水面图，计算不透水面盖度。

所述步骤1中，收集Landsat系列影像和Sentinel-2等遥感影像，分别对这些影像进行大气校正，将原始影像的灰度值转换为地表反射率值。其中，对于Sentinel-2影像，可以采用sen2cor模型进行大气校正；对于Landsat系列影像，采用本发明提出的简化大气校正方法，其主要的推导过程如下：

(1)基于Landsat-5 TM的简化TOA大气校正模型

原始的大气层顶表观反射率TOA模型，其计算公式如式(1)-(2)：

ρ_TOA＝L_sat·π·d²/(E₀·cosθ) (1)

L_sat＝Gain·DN+Bias (2)

式中，DN(或QCAL)为经过辐射校正的图像灰度值，π为常量；E₀为大气顶层太阳辐照度；θ为太阳天顶角(与太阳高度角互余)；d为日地天文距离，L_sat为传感器所接受的辐射能量。Bias和Gain分别为图像辐亮度值的偏移参数和增益参数。各波段光谱通道的Bias和Gain可以从遥感头文件中获取，d值和θ也在影像头文件获得。

把式(2)代入式(1)，并设K_TOA＝π·d²/E₀，可以得到式(3)：

ρ_TOA＝K_TOA(Gain·DN+Bias)/cosθ (3)

式(3)可简化为式(4)：

ρ_TOA＝(M_TOA·DN+A_TOA)/cosθ (4)

式(4)中，

在2016年后的Landsat 5 TM头文件中，不仅给出了增益Gain和偏置Bias的参数信息，同时增加了反射率定标参数Reflectance_mult_band和Reflectance_add_band，这两个反射率定标参数分别与公式(4)中的M_TOA、A_TOA对应。因此从头文件中可以直接读取M_TOA、A_TOA文件，公式(4)简化了计算过程。Landsat 8 OLI卫星影像提供了新的表观反射率模型，如式(5)：

ρ_L8TOA＝(M_p·DN+A_p)/cosθ (5)

式(5)中，M_p和A_p分别表示各波段反射率的增益参数和偏移参数，在遥感影像提供的头文件中读取。可以看出公式(4)和(5)是一致的，本发明把Landsat5 TM和Landsat 8两种影像的TOA大气校正模型进行统一。

(2)基于Landsat-5 TM和Landsat 8 OLI的简化COST大气校正模型

COST模型考虑了大气散射和吸收，校正了大气程辐射和部分大气透过率造成的影响，非常适合于历史遥感数据的处理。但原始的COST大气校正模型由公式(6)-(9)组成，计算过程复杂。

Lhaze＝L_min-L_1％ (7)

L_min＝Gain·DN_min+Bias (8)

公式(6)-(9)中，Lhaze为大气层辐射值；DN_min为遥感器的最小光谱辐射值，L_min为影像直方图频数累加和等于总像元个数0.01％的像元DN值所对应的辐射值；L_1％为假设黑体反射率为1％各波段的黑体辐射值；TAU_z为从地面到传感器的大气透射率(设置为常数1或cosθ)；TAU_v为从太阳到地面的大气透射率(设置为常数1)。

设

则式(6)和(9)可分别写为式(10)和(11)，把式(8)、(9)和(11)代入(7)得式(12)：

ρ_COST＝(L_sat-L_haze)·K_COST/cosθ (10)

因此把式(12)代入式(10)，进一步计算可得式(13)：

式(13)可以简化为：

ρ_COST＝M_COST·(DN-DN_min)/cosθ+0.01 (14)

式(14)中，M_COST＝K_COST·Gain＝M_TOA/(TAU_z·TAU_v)。

同样，基于Landsat 8 OLI影像数据的COST模型，即公式(15)

ρ_L8COST＝M_COST·(DN-DN_min)/cosθ+0.01 (15)

式(15)中，M_COST＝M_p/(TAU_z·TAU_v)。

因此对于Landsat 5 TM和Landsat 8 OLI可以采用统一的COST大气校正模型

ρ_COST＝M_COST·(DN-DN_min)/cosθ+0.01 (16)

M_COST＝M/(TAU_z·TAU_v) (17)

式(16)-(17)中，M是影响各波段的反射率调整参数，表示各波段反射率的增益参数，可以从头文件中可以直接读取。

原始的COST模型中分3步计算大气校正，其过程包括：(a)根据遥感器的增益与偏移进行遥感器定标；(b)遥感器的光谱辐射值转换为遥感器的相对反射值；(c)大气纠正，消除因大气吸收和散射造成的大气影响，计算地球表面像元相对反射率。但本发明修订的COST模型，针对Landsat 5 TM和Landsat 8 OLI影像，把3个计算过程合并到一个模型中，该方法也把Landsat 5 TM和Landsat 8 OLI两种影像的大气校正方法统一到同一个模型中。

所述步骤2中，遥感分类指数集包括不透水面指数BCI，水体指数NDWI，修订的水体指数MNDWI，裸土指数BSI，新裸土指数MFI，湿度指数WI，亮度指数BI，绿度指数GI等；

MNDWI＝(ρ_Green-ρ_SWIR1)/(ρ_Green+ρ_SWIR1) (18)

NDWI＝(ρ_Green-ρ_NIR)/(ρ_Green+ρ_NIR) (19)

式中，ρ_Green、ρ_NIR和ρ_SWIR1分别为经过大气校正的绿光波段、近红外波段和短波红外波段1的反射率。

BCI指数的计算公式为：

BCI＝[(H+L)/2-V]/[(H+L)/2+V] (20)

式中，H、L、V分别为归一化亮度分量(高反射率)，归一化湿度分量(低反射率)，归一化植被绿度分量。这三个分量的计算公式为：

H＝(BI-BI_min)/(BI_max-BI_min) (21)

V＝(GI-GI_min)/(GI_max-GI_min) (22)

L＝(WI-WI_min)/(WI_max-WI_min) (23)

对于Landsat 5 TM影像，采用式(24)-(26)计算H、V和L三个分量：

对于Landsat 8 OLI影像，采用式(27)-(29)计算H、V和L三个分量：

对于Sentinel-2影像，采用式(30)-(32)计算H、V和L三个分量：

式(24-32)中，BI_TM，GI_TM，WI_TM，BI_OLI，GI_OLI，WI_OLI，BI_{Sen_2}，GI_{Sen_2}，WI_{Sen_2}分别表示Landat-5 TM，Landsat 8 OLI和Sentinel-2遥感影像的缨帽变换亮度指数，绿度指数和湿度指数；ρ_Coastal、ρ_Blue、ρ_Green、ρ_Red、ρ_VRE1、ρ_VRE2、ρ_VRE3、ρ_NIR、ρ_VRE4、ρ_Water、ρ_Cirrus、ρ_SWIR1和ρ_SWIR2分别为经过大气校正的对应遥感影像的气溶胶波段、蓝光波段、绿光波段、红光波段、红边1波段、红边2波段、红边3波段、近红外波段、红边4波段、水气波段、卷云波段、短波红外波段1和短波红外波段2的反射率。

式(33)的裸土指数(bare soil index,BSI)中，ρ_Blue和ρ_Red分别是蓝光波段和红光波段的反射率，ρ_NIR和ρ_SWIR1分别为近红外波段和中红外波段1的反射率。

式(34)-(35)中，ρ_VRGEi表示Sentinel-2影像中四个植被红边波段的反射率，λ_i表示各个植被红边波段的中心波长(705,740,783和865nm)，ρ_Bλi指四个植被红边波段的基准反射率。MFI指数曾被提出用于提取沿海红树林的提取，本发明发现该指数可很好地提取Sentinel-2影像中的裸土信息。

所述步骤3的具体实施过程如下：

(31)基于水体指数提取城市水体，区分水体与非水体；

(32)基于不透水面指数提取非植被，区分植被与非植被；

(33)基于裸土指数提取裸土，区分裸土与非裸土；

(34)基于湿度指数提取滩涂湿地，区分滩涂湿地与不透水面；

(35)基于湿度指数变化检测，进一步区分漏分的滩涂湿地或误分的高楼影子。

步骤(31)中，以基于简化的COST大气校正模型计算的遥感影像为根节点，将MNDWI≥P₁或NDWI≥P₁₁遥感指数作为约束条件，P₁或P₁₁为阈值，完成1级分类，城市土地覆盖类型中满足约束条件的区域为水体，否则为非水体；

步骤(32)中，以1级别判断中的非水体为分支节点，将BCI≥P₂作为约束条件，P₂为阈值，完成2级分类，城市土地覆盖类型中满足约束条件的区域为非植被，否则为植被；

步骤(33)中，以2级别判断中的非植被为分支节点，将BSI≥P₃或MFI≥P₃₁作为约束条件，P₃或P₃₁为阈值，完成3级分类，城市土地覆盖类型中满足约束条件的区域为非裸土，否则为裸土；

步骤(34)中，以3级别判断中的非裸土为分支节点，将WI≥P₄作为约束条件，P₄为阈值，完成4级分类，城市土地覆盖类型中满足约束条件的区域为滩涂湿地，否则为不透水面；

步骤(35)中，提取漏分的滩涂湿地和高楼影子，利用两期影像A(前期影像)和B(分类影像)，在提取影像A的水体基础上，把两期WI指数影像掩膜，做两期WI影像的差值，进行WI指数的变化检测，判断漏分的滩涂湿地和误分的高楼影子。

所述步骤4中，Google Earth高清晰影像中采集地面验证点，其做法是从地图下载器下载16级以上的Google Earth高分辨率影像数据(避免了影像配准步骤)，并进行重采样为1m。采用分层随机采样法为每种土地覆盖类型抽取至少60个样本点，建立混淆矩阵，计算总体分类精度与Kappa系数等指标。如果分类精度不满足要求(如设定总体精度低于90％)，需要调整决策树中的关键参数。

所述步骤5中，对BCI指数进行规一化，计算不透水面盖度ISA％，公式如下：

式(36)中，BCI_max、BCI_min分别为归一化不透水面指数的最大值和最小值；ISA％为不透水面盖度，表征不透水面的聚集密度，范围在0～1之间，ISA％越大，说明不透水面密度越高，0值说明该像元范围内无不透水面下垫面。实际应用中，高不透水面区的BCI_max取值如大型建筑物屋顶，大型停车场或港口码头等，低不透水面区的BCI_min取值如浓密的森林。

以下选取厦门市为实施案例，采用2个时相Landsat 5 TM和Landat 8 OLs遥感影像和1个Sentinel-2影像为数据源，30m空间分辨率的Landsats影像成像时间分别为2009/06/06和2018/03/11(范围为整个厦门市)，10m空间分辨率的Sentinel-2影像时间为2019/12/17(范围为厦门市主城区)。以此数据验证沿海城市时间不透水面信息提取方法的有效性。

案例实施过程中将根据研究区厦门市土地覆盖类型，将其划分为四大类：水体(含滩涂湿地)、植被、不透水面和裸地。应用本发明所述方法对厦门市3个时相的遥感影像数据进行了不透水面信息提取，并将其结果原始的BCI指数提取结果进行比较。需要说明的是原始的BCI指数法无法区分裸土、滩涂及细小的湿地类型。

具体实施步骤如下：

步骤一、对原始影像进行数据预处理，本发明技术方案中公式(16-17)，简化的大气校正COST模型进行预处理；

步骤二、对经过大气校正的遥感影像，按本发明技术方案中公式(18-35)，分别计算出对应的分类特征指数。对应Landsats系列影像，采用MNDWI、BCI、BSI、WI的组合；对应Sentinel-2影像，采用NDWI、BCI、MFI和WI指数的组合，从而构建遥感分类特征指数数据库；

步骤三、建立基于多特征的决策树分类模型(见图2)，模型中涉及的阈值P1-P5的选取如表1所示。以2018年厦门市Landsat 8影像提取不透水面模型的建立过程详细说明：

第一步，提取水体。应用MNDWI指数区分出水体与非水体，MNDWI≥0为水体；需要说明的是:对于Landsats影像采用MNDWI指数，而对于Sentinel-2影像采用NDWI指数(由于其空间分辨率比较高)；单独采用MNDWI或NDWI指数，存在滩涂湿地漏提的现象，特别是随着城市化的快速推进，滩涂湿地变得更加破碎化，导致不易辨识；因此在非水体区域中，存在着一些滩涂湿地，它们与不透水面信息相混淆。

第二步，在非水体中，应用BCI指数区分植被与非植被，BCI≥-0.07为非植被区；

第三步，在非植被区中，联合BSI指数和BCI指数区分裸地与不透水面，设置BSI＜0.15且BCI＜0.4为裸地。需要特别指出的是，单独的BSI指数无法有效提取裸土，因为部分不透水面区与裸土呈现相似的光谱；单独的BCI指数也无法区分裸土，尽管通常的裸土具有较低的BCI值，但并非都如此。

第四步，在非裸土中，区分不透水面与滩涂湿地。滩涂湿地是城市特别是沿海城市的特点之一，但它们不易区分。在WI湿度指数图中，相比不透水面，滩涂湿地具有较高的湿度指数值，设置WI≥-0.03为滩涂湿地，否则为不透水面。

第五步，单独的WI指数往往效果不理想，仍然有少量的细小滩涂湿地混淆在不透水面中，这时可以采用两期影像，首先获得两期影像的水体区域的WI指数图，然后进行WI图的变化监测。本例中，获得2009年和2018年的水体区域的WI图分别取名为WI09和WI18，然后进行WI指数图的变化监测，从而将漏分的滩涂湿地进行正确归类。

表1 3个时相厦门市影像在多特征决策树分类模型中的阈值表

步骤四、以具有高空间分辨率的同期的谷歌影像数据，提取800个采样点，作为参考影像数据，分别对案例实施过程中多个时相的影像提取的不透水面分类结果进行精度评价，并将其结果与原始的BCI方法提取结果进行比较；输出满足精度的城市不透水面图(见图3)。

测试结果表明，运用本发明进行厦门市不透水面信息提取，提取结果的精度明显高于原始的BCI方法，可参见表2。基于多特征决策树模型三个时相的不透水面分类结果中，平均的总体分类精度和Kappa系数分别为92.52％和0.89，相比原始BCI指数法，总体分类精度提高了7.17％，Kappa系数提高了0.11。

表2 多特征决策树方法与原始BCI法提取不透水面精度评价对比

步骤五、计算研究区不同时相的不透水面盖度。

综上，本发明首先提出一种简化的大气校正方法，避免了原始气校正方法的繁琐复杂的计算过程；该方法采用多特征决策树模型分层提取土地利用覆被信息，避免了单纯利用不透水面指数单一分类特征的缺陷，其具有一定的通用性，极大提高了不透水面信息提取精度，解决了城市不透水面信息提取的关键技术问题。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (10)

1.一种城市不透水面信息快速提取方法，其特征在于包括如下步骤：

步骤1，收集指定研究区的Landsat 5 TM、Landsat 8OLI或Sentinel-2遥感影像中的任意一种，然后对收集的影像进行大气校正，将原始影像的灰度值转换为地表反射率值；

步骤2，根据遥感影像类型，选择一组遥感分类特征，分别计算对应的分类特征指数；

步骤3，基于遥感分类特征和城市土地覆盖类型，建立单规则分类或多规则分类指数的决策规则，并依据规则进行城市不透水面信息逐级分类；

步骤4，对提取的城市不透水面信息进行精度评价，输出满足精度的城市不透水面图。

2.如权利要求1所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤1中，对于Sentinel-2影像，采用sen2cor模型进行大气校正；对于Landsat 5 TM、Landsat 8OLI影像，采用如下大气校正方法：

ρ_COST＝M_COST·(DN-DN_min)/cosθ+0.01

M_COST＝M/(TAU_z·TAU_v)

其中，DN为经过辐射校正的图像灰度值，DN_min为遥感器的最小光谱辐射值，θ为太阳天顶角，M是影响各波段的反射率调整参数，TAU_z为从地面到传感器的大气透射率；TAU_v为从太阳到地面的大气透射率；ρ_COST表示Landsats系列影像的地面反射率；M_COST表示大气校正参数。

3.如权利要求1所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤2中，若遥感影像类型为Landsat 5 TM或Landsat 8 OLI影像，选择的遥感分类特征包括不透水面指数BCI，修订的水体指数MNDWI，裸土指数BSI，湿度指数WI，亮度指数BI和绿度指数GI；若遥感影像类型为Sentinel-2遥感影像，选择的遥感分类特征包括不透水面指数BCI，水体指数NDWI，新裸土指数MFI，湿度指数WI，亮度指数BI和绿度指数GI。

4.如权利要求3所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤2中，水体指数MNDWI、修订的水体指数MNDWI的计算公式是：

MNDWI＝(ρ_Green-ρ_SWIR1)/(ρ_Green+ρ_SWIR1)

NDWI＝(ρ_Green-ρ_NIR)/(ρ_Green+ρ_NIR)

式中，ρ_Green、ρ_NIR和ρ_SWIR1分别为经过大气校正的绿光波段、近红外波段和短波红外波段1的反射率。

5.如权利要求3所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤2中，不透水面指数BCI的计算公式是：

BCI＝[(H+L)/2-V]/[(H+L)/2+V]

其中，H、L、V分别为归一化亮度分量、归一化湿度分量、归一化植被绿度分量。

6.如权利要求5所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤2中，H、L、V的计算公式为：

H＝(BI-BI_min)/(BI_max-BI_min)

V＝(GI-GI_min)/(GI_max-GI_min)

L＝(WI-WI_min)/(WI_max-WI_min)

其中，BI、GI、WI分别为亮度指数、绿度指数、湿度指数，BI_max、GI_max、WI_max分别为对应的亮度指数图像、绿度指数图像、湿度指数图像中的最大值，BI_min、GI_min、WI_min分别为对应的亮度指数图像、绿度指数图像、湿度指数图像中的最小值。

7.如权利要求3所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤2中，裸土指数BSI的计算公式是：

BSI＝[(ρ_Red+ρ_SWIR1)-(ρ_NIR+ρ_Blue)]/[(ρ_Red+ρ_SWIR1)+(ρ_NIR+ρ_Blue)]

其中，ρ_Green和ρ_Red分别是绿光波段和红光波段的反射率，ρ_NIR和ρ_SWIR1分别为近红外波段和中红外波段1的反射率。

8.如权利要求3所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤3中，单规则分类或多规则分类指数的决策规则的具体内容是：

基于MNDWI或NDWI提取城市水体，区分水体与非水体；

基于亮度指数、绿度指数和湿度指数，计算不透水面指数；基于不透水面指数提取非植被，区分植被与非植被；

基于BSI或MFI提取裸土，区分裸土与非裸土；

基于湿度指数提取滩涂湿地，区分滩涂湿地与不透水面；

基于湿度指数变化检测，区分漏分的滩涂湿地或误分的高楼的影子。

9.如权利要求1所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤4的具体内容是：所述步骤4中，从地图下载器下载16级以上的Google Earth高分辨率影像数据，并进行重采样为1m；采用分层随机采样法为每种土地覆盖类型抽取N个样本点，建立混淆矩阵，计算总体分类精度与Kappa系数，如果分类精度不满足要求，则调整决策树中的关键参数。

10.如权利要求3所述的一种城市不透水面信息快速提取方法，其特征在于：所述步骤4后，还包括步骤5：根据步骤4得到的城市不透水面图，对不透水面指数BCI进行规一化，根据下式计算不透水面盖度ISA％：

其中，BCI_max、BCI_min分别为归一化不透水面指数的最大值和最小值；ISA％为不透水面盖度，表征不透水面的聚集密度，范围在0～1之间，ISA％越大，说明不透水面密度越高，0值说明该像元范围内无不透水面下垫面。

Images