CN114266968A - 一种城市不同土地覆盖类型遥感自动解译方法 - Google Patents

Abstract

本发明提出了一种城市不同土地覆盖类型遥感自动解译方法，其步骤为：首先，通过云计算平台获取X年的N幅Landsat卫星影像的地表反射率图像，并利用光谱采集软件采集城市不同土地覆盖类型的光谱信息；其次，分析不同地物之间的光谱差异，确定识别城市不同土地覆盖类型的有效光谱波段及其图像合成方案；然后，创建蓝色彩钢板房遥感指数、红色不透水面遥感指数；再基于数据处理后的各种图像和阈值，创建城市不同土地覆盖类型遥感解译模型，实现城市不同土地覆盖类型的快速分类；最后对分类结果进行孤立像元合并处理，得到城市不同土地覆盖类型图。本发明提出的城市不同土地覆盖类型遥感解译模型，可实现全球尺度上的城市用地类型的自动识别制图。

Description

一种城市不同土地覆盖类型遥感自动解译方法

技术领域

本发明涉及遥感目标识别技术领域，特别是指一种城市不同土地覆盖类型遥感自动解译方法。

背景技术

随着社会经济的快速发展和城镇化率的不断提高，城市用地规模呈现快速扩张的态势，进而造成城市用地与农业用地的矛盾日益突出。及时有效地监测城市用地规模和土地覆盖类型对研究评价城市生态环境、城市发展质量，以及城市规划等宏观决策具有重大意义。但是目前城市不同土地覆盖类型制图研究存在精度低、分类粗糙、更新时效慢等问题。

发明内容

针对上述技术背景中存在的不足，本发明提出了一种城市不同土地覆盖类型遥感自动解译方法，解决了城市不同土地覆盖类型制图研究存在自动化水平低、分类粗糙、更新时效慢等问题。

本发明的技术方案是这样实现的：

一种城市不同土地覆盖类型遥感自动解译方法，其步骤如下：

S1、利用遥感云平台采集X年的N幅Landsat卫星影像，并对N幅Landsat卫星影像进行数据预处理和云层剔除处理，得到N幅地表反射率图像；其中，Landsat卫星影像包括蓝色波段、绿色波段、红色波段、近红外波段、第一短波红外波段和第二短波红外波段这6个波段影像；

S2、基于地面调查数据确定城市包含的土地覆盖类型，基于GPS定位工具获取土地覆盖类型的地理坐标，并使用光谱采集软件采集分析地表反射率图像中的土地覆盖类型的像元光谱样本，得到用于识别土地覆盖类型的有效光谱波段；其中，土地覆盖类型包括蓝色彩钢板房、红色不透水面、水泥不透水面、其他不透水面、植被、水体、未识别类型这7种类型；

S3、基于水体遥感指数和步骤S1中得到的N幅地表反射率图像计算得到N幅水体遥感指数图像，采用均值图像重构方法对N幅水体遥感指数图像进行图像重构，得到X年水体遥感指数均值图像，记为WI图像；基于步骤S2中地面调查数据中包含的水体样本统计分析水体在WI图像中的分布区间，取其最小值作为阈值α；

S4、基于步骤S1中得到的N幅地表反射率图像中的第二短波红外波段图像，采用四分位数图像重构方法对N幅第二短波红外波段图像进行图像重构，取其重构的第零个四分位数图像为X年第二短波红外波段最小值图像，记为SWIRmin图像，取其重构的第二个四分位数图像为X年第二短波红外波段中值图像，记为SWIRM图像；基于步骤S2中地面调查数据中包含的水体样本统计分析水体在SWIRmin图像中的分布区间，取其最大值作为阈值β；基于步骤S2中地面调查数据中包含的红色不透水面样本统计分析红色不透水面在SWIRM图像中的分布区间，取其最小值作为阈值κ；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在SWIRM图像中的分布区间，取其最小值作为阈值ν；

S5、创建蓝色彩钢板房遥感指数，基于步骤S1中得到的N幅地表反射率图像计算得到N幅蓝色彩钢板房遥感指数图像，然后采用四分位数图像重构方法对N幅蓝色彩钢板房遥感指数图像进行图像重构，取其重构的第三个四分位数图像为X年蓝色彩钢板房遥感指数图像，记为BI图像；基于步骤S2中地面调查数据中包含的蓝色彩钢板房样本统计分析蓝色彩钢板房在BI图像中的分布区间，取其最小值作为阈值γ；

S6、基于步骤S1中得到的N幅地表反射率图像中的蓝色波段图像，采用四分位数图像重构方法对N幅蓝色波段图像进行图像重构，取其重构的二个四分位数图像为X年蓝色波段中值图像，记为BM图像；基于步骤S2中地面调查数据中包含的蓝色彩钢板房样本统计分析蓝色彩钢板房在BM图像中的分布区间，取其最小值作为阈值δ；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在BM图像中的分布区间，取其最小值作为阈值μ；

S7、基于植被遥感指数和步骤S1中得到的N幅地表反射率图像计算得到N幅植被遥感指数图像，采用四分位数图像重构方法对N幅植被遥感指数图像进行图像重构，取其重构的第四个四分位数图像为X年植被遥感指数最大值图像，记为VI图像；基于步骤S2中地面调查数据中包含的植被样本统计分析纯净像元的植被在VI图像中的分布区间，取其最小值作为阈值ε；统计分析混合像元的植被在VI图像中的分布区间，取其最小值作为阈值ζ；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析纯净像元的其他不透水面在VI图像中的分布区间，取其最大值作为阈值π；

S8、基于步骤S1中得到的N幅地表反射率图像计算得到N幅蓝色波段与红色波段的和值图像，采用四分位数图像重构方法对N幅差值图像进行图像重构，取其重构的第二个四分位数图像为X年蓝色波段与红色波段的和值图像，记为RBsum图像；基于步骤S2中地面调查数据中包含的植被样本统计分析植被在RBsum图像中的分布区间，取其最大值作为阈值η；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析其他不透水面在RBsum图像中的分布区间，取其最小值作为阈值σ；

S9、基于步骤S1中得到的N幅地表反射率图像中的近红外波段图像，采用四分位数图像重构方法对N幅近红外波段图像进行图像重构，取其重构的四个四分位数图像为X年近红外波段最大值图像，记为NIRmax图像；基于步骤S2中地面调查数据中包含的植被样本统计分析植被在NIRmax图像中的分布区间，取其最小值作为阈值θ；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析其他不透水面在NIRmax图像中的分布区间，取其最大值作为阈值ο；

S10、创建红色不透水面遥感指数，基于步骤S1中得到的N幅地表反射率图像计算得到N幅红色不透水面遥感指数图像，然后采用四分位数图像重构方法对N幅红色不透水面遥感指数图像进行图像重构，取其重构的第三个四分位数图像为X年红色不透水面遥感指数图像，记为RI图像；基于步骤S2中地面调查数据中包含的红色不透水面样本统计分析红色不透水面在RI图像中的分布区间，取其最小值作为阈值ι；

S11、基于步骤S1中得到的N幅地表反射率图像中的红色波段图像，采用四分位数图像重构方法对N幅红色波段图像进行图像重构，取其重构的二个四分位数图像为X年红色波段中值图像，记为RM图像；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在RM图像中的分布区间，取其最小值作为阈值ξ；

S12、判断像元i在步骤S3中得到的水体遥感指数均值图像中的像元值是否大于阈值α，若不是，执行步骤S13，若是，判断像元i在步骤S4中得到的第二短波红外波段最小值图像中的像元值是否小于阈值β，若是，像元i为水体，若不是，执行步骤S13；

S13、判断像元i在步骤S5中得到的蓝色彩钢板房遥感指数图像中的像元值是否大于阈值γ，若不是，执行步骤S14，若是，判断像元i在步骤S6中得到的蓝色波段中值图像中的像元值是否大于阈值δ，若是，像元i为蓝色彩钢板房，若不是，执行步骤S14；

S14、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否大于阈值ε，若不是，执行步骤S15，若是，像元i为植被；

S15、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否属于区间范围[ζ，ε]，若不是，执行步骤S16，若是，判断像元i在步骤S8中得到的蓝色波段与红色波段的和值图像中的像元值是否小于阈值η，若不是，执行步骤S16，若是，判断像元i在步骤S9中得到的近红外波段最大值图像中的像元值是否大于阈值θ，若不是，执行步骤S16，若是，像元i为植被；

S16、判断像元i在步骤S10中得到的红色不透水面遥感指数图像中的像元值是否大于阈值ι，若不是，执行步骤S17，若是，判断像元i在步骤S4中得到的第二短波红外波段中值图像中的像元值是否大于阈值κ，若不是，执行步骤S17，若是，像元i为红色不透水面；

S17、判断像元i在步骤S6中得到的蓝色波段中值图像中的像元值是否大于阈值μ，若不是，执行步骤S18，若是，判断像元i在步骤S4中得到的第二短波红外波段中值图像中的像元值是否大于阈值ν，若不是，执行步骤S18，若是，判断像元i在步骤S11中得到的红色波段中值图像中的像元值是否大于阈值ξ，若不是，执行步骤S18，若是，像元i为水泥不透水面；

S18、判断像元i在步骤S9中得到的近红外波段最大值图像中的像元值是否小于阈值ο，若不是，像元i为未识别类型，若是，执行步骤S19；

S19、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否小于阈值π，若是，像元i为其他不透水面，若不是，执行步骤S20；

S20、判断像元i在步骤S8中得到的蓝色波段与红色波段的和值图像中的像元值是否大于阈值σ，若是，像元i为其他不透水面，若不是，像元i为未识别类型；

S21、循环执行步骤S12至步骤S20，直至遍历完图像中的所有像元位置，完成研究区域的城市不同土地覆盖类型的遥感识别分类；

S22、采用3×3运算窗口对步骤S21中得到的分类结果进行滤波处理，消除分类结果中的孤立像元，完成最终的城市不同土地覆盖类型的遥感识别分类。

步骤S1中对Landsat卫星影像进行数据预处理和云层剔除处理的方法为：对Landsat卫星影像中所有像元值乘以修正系数0.0000275，再减去0.2，得到Landsat卫星影像的地表反射率图像；在地表反射率图像中像元i的位置上，如果蓝色波段、绿色波段、红色波段上像元值的极差的绝对值小于0.08，并且绿色波段上的像元值大于0.25，那么像元i的位置上的地物类型为云层，则将像元i的位置上所有波段的像元值改写为空值，实现云层剔除处理。

所述采用均值图像重构方法对N幅水体遥感指数图像进行图像重构的方法为：首先创建1幅像元行列数与水体遥感指数图像的像元行列数一致的空值图像；然后将N幅水体遥感指数图像进行图层堆叠，像元i的位置上存在n个水体遥感指数值，计算得到这n个数值的平均值，并写入空值图像的相应位置，依次遍历所有像元位置，得到水体遥感指数均值图像。

所述采用四分位数图像重构方法对N幅第二短波红外波段图像进行图像重构的方法为：首先创建5幅像元行列数与第二短波红外波段图像的像元行列数一致的空值图像；然后将N幅第二短波红外波段图像进行图层堆叠，像元i的位置上存在n个第二短波红外波段像元值，对这n个值按照从小到大的顺序排列成数列，分别将数列中位于第零、一、二、三、四的四分位点的值写入空值图像中的相应像元位置，依次遍历所有像元位置，即分别得到第二短波红外波段图像的第零、一、二、三、四的四分位数图像；其中第零、二、四的四分位点分别表示数组的最小值、中值、最大值。

所述步骤S5中创建蓝色彩钢板房遥感指数的方法为：

BPI_i＝(ξ_blue,i-ξ_red,i)/(ξ_blue,i+ξ_red,i)；

其中，BPI_i表示蓝色彩钢板房遥感指数图像中像元i上的蓝色彩钢板房遥感指数值，ξ_blue,i表示地表反射率图像中像元i上的蓝色波段的地表反射率，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

所述步骤S8中创建红色不透水面遥感指数的方法为：

RRI_i＝(ξ_red,i-ξ_green,i)/(ξ_red,i+ξ_green,i)；

其中，RRI_i表示红色不透水面遥感指数图像中像元i上的红色不透水面遥感指数值，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，ξ_green,i表示地表反射率图像中像元i上的绿色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

所述步骤S22中消除分类结果中孤立像元的方法为：对于像元(x，y)，图像中最大行数为X，最大列数为Y，其中x＝2，3，……，X-1；y＝2，3，……，Y-1；统计以像元(x，y)为中心的3×3窗口范围内不同类别的像元个数，设置像元(x，y)的数量权重为5，与像元(x，y)类别一致的像元个数为5+m，m为3×3窗口中与像元(x，y)类别一致的像元个数，令J表示为3×3窗口中与像元(x，y)类别不一致的类别中的最大像元个数，如果5+m小于J，则像元(x，y)的类别属性改写为与J类别一致的类别属性，否则，像元(x，y)的类别属性不变；依次遍历图像中所有可用像元位置，迭代3次，完成分类结果中孤立像元的消除。

与现有技术相比，本发明产生的有益效果为：

(1)本发明提出了基于光学卫星影像的蓝色彩钢板房遥感指数Blue Plate Index(BPI)、红色不透水面遥感指数Red Surface Index(RSI)，为城市不透水面遥感识别提供了理论基础；

(2)本发明创建的城市不同土地覆盖类型遥感自动解译模型算法，充分利用了不同土地覆盖类型的光谱特征，能够适用于全球尺度上城市不同土地覆盖类型的快速识别，为相关部门、行业提供精准的城市用地数据。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图；

图2为本发明的开封市不同土地覆盖类型识别结果。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明实施例提供了一种城市不同土地覆盖类型遥感自动解译方法，其步骤如下：

S1、利用遥感云平台采集2021年5月至9月份的N(N＝15)幅Landsat卫星影像，并对N幅Landsat卫星影像进行数据预处理和云层剔除处理，得到N幅地表反射率图像；其中，Landsat卫星影像包括蓝色波段、绿色波段、红色波段、近红外波段、第一短波红外波段和第二短波红外波段这6个波段影像；

对Landsat卫星影像进行数据预处理的方法为：对Landsat卫星影像中所有像元值乘以修正系数0.0000275，再减去0.2，得到Landsat卫星影像的地表反射率图像；实现数据预处理。

对Landsat卫星影像进行云层剔除处理的方法为：在地表反射率图像中像元i的位置上，如果蓝色波段、绿色波段、红色波段上像元值的极差的绝对值小于0.08，并且绿色波段上的像元值大于0.25，那么像元i的位置上的地物类型为云层，则将像元i的位置上所有波段的像元值改写为空值，实现云层剔除处理。

S2、基于地面调查数据确定城市包含的土地覆盖类型，基于GPS定位工具获取土地覆盖类型的地理坐标，并使用光谱采集软件采集分析地表反射率图像中的土地覆盖类型的像元光谱样本，得到用于识别土地覆盖类型的有效光谱波段；其中，土地覆盖类型包括蓝色彩钢板房、红色不透水面、水泥不透水面、其他不透水面、植被、水体、未识别类型这7种类型；

S3、基于水体遥感指数和步骤S1中得到的N幅地表反射率图像计算得到N幅水体遥感指数图像，采用均值图像重构方法对N幅水体遥感指数图像进行图像重构，得到2021年水体遥感指数均值图像，记为WI图像；基于步骤S2中地面调查数据中包含的水体样本统计分析水体在WI图像中的分布区间，取其最小值作为阈值α，α＝0。

所述水体遥感指数的计算方法为：

WI_i＝(ξ_green,i-ξ_nir,i)/(ξ_green,i+ξ_nir,i)；

其中，WI_i表示蓝色彩钢板房遥感指数图像中像元i上的蓝色彩钢板房遥感指数值，ξ_green,i表示地表反射率图像中像元i上的蓝色波段的地表反射率，ξ_nir,i表示地表反射率图像中像元i上的红色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

所述采用均值图像重构方法对N幅水体遥感指数图像进行图像重构的方法为：首先创建1幅像元行列数与水体遥感指数图像的像元行列数一致的空值图像；然后将15幅水体遥感指数图像进行图层堆叠，那么像元位置i上存在n个水体遥感指数值，计算得到这n个数值的平均值，并写入空值图像的相应位置，依次遍历所有像元位置，即得到水体遥感指数均值图像。

S4、基于步骤S1中得到的N幅地表反射率图像中的第二短波红外波段图像，采用四分位数图像重构方法对N幅第二短波红外波段图像进行图像重构，取其重构的第零个四分位数图像为2021年第二短波红外波段最小值图像，记为SWIRmin图像，取其重构的第二个四分位数图像为2021年第二短波红外波段中值图像，记为SWIRM图像；基于步骤S2中地面调查数据中包含的水体样本统计分析水体在SWIRmin图像中的分布区间，取其最大值作为阈值β，β＝0.1；基于步骤S2中地面调查数据中包含的红色不透水面样本统计分析红色不透水面在SWIRM图像中的分布区间，取其最小值作为阈值κ，κ＝0.17；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在SWIRM图像中的分布区间，取其最小值作为阈值ν，ν＝0.17。

所述采用四分位数图像重构方法对N幅第二短波红外波段图像进行图像重构的方法为：首先创建5幅像元行列数与第二短波红外波段图像的像元行列数一致的空值图像；然后将N幅第二短波红外波段图像进行图层堆叠，像元i的位置上存在n个第二短波红外波段像元值，对这n个值按照从小到大的顺序排列成数列，分别将数列中位于第零、一、二、三、四的四分位点的值写入空值图像中的相应像元位置，依次遍历所有像元位置，即分别得到第二短波红外波段图像的第零、一、二、三、四的四分位数图像；其中第零、二、四的四分位点分别表示数组的最小值、中值、最大值。

S5、创建蓝色彩钢板房遥感指数，基于步骤S1中得到的N幅地表反射率图像计算得到N幅蓝色彩钢板房遥感指数图像，然后采用四分位数图像重构方法对N幅蓝色彩钢板房遥感指数图像进行图像重构，取其重构的第三个四分位数图像为2021年蓝色彩钢板房遥感指数图像，记为BI图像；基于步骤S2中地面调查数据中包含的蓝色彩钢板房样本统计分析蓝色彩钢板房在BI图像中的分布区间，取其最小值作为阈值γ，γ＝0。

所述创建蓝色彩钢板房遥感指数的方法为：

BPI_i＝(ξ_blue,i-ξ_red,i-0.01)/(ξ_blue,i+ξ_red,i)；

其中，BPI_i表示蓝色彩钢板房遥感指数图像中像元i上的蓝色彩钢板房遥感指数值，ξ_blue,i表示地表反射率图像中像元i上的蓝色波段的地表反射率，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

S6、基于步骤S1中得到的N幅地表反射率图像中的蓝色波段图像，采用四分位数图像重构方法对N幅蓝色波段图像进行图像重构，取其重构的二个四分位数图像为2021年蓝色波段中值图像，记为BM图像；基于步骤S2中地面调查数据中包含的蓝色彩钢板房样本统计分析蓝色彩钢板房在BM图像中的分布区间，取其最小值作为阈值δ，δ＝0.1；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在BM图像中的分布区间，取其最小值作为阈值μ，μ＝0.1。

S7、基于植被遥感指数和步骤S1中得到的N幅地表反射率图像计算得到N幅植被遥感指数图像，采用四分位数图像重构方法对N幅植被遥感指数图像进行图像重构，取其重构的第四个四分位数图像为2021年植被遥感指数最大值图像，记为VI图像；基于步骤S2中地面调查数据中包含的植被样本统计分析纯净像元的植被在VI图像中的分布区间，取其最小值作为阈值ε，ε＝0.6；统计分析混合像元的植被在VI图像中的分布区间，取其最小值作为阈值ζ，ζ＝0.5；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析纯净像元的其他不透水面在VI图像中的分布区间，取其最大值作为阈值π，π＝0.4。

所述植被遥感指数的计算方法为：

VI_i＝(ξ_nir,i-ξ_red,i)/(ξ_red,i+ξ_nir,i)；

其中，VI_i表示蓝色彩钢板房遥感指数图像中像元i上的蓝色彩钢板房遥感指数值，ξ_nir,i表示地表反射率图像中像元i上的蓝色波段的地表反射率，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

S8、基于步骤S1中得到的N幅地表反射率图像计算得到N幅蓝色波段与红色波段的和值图像，采用四分位数图像重构方法对N幅差值图像进行图像重构，取其重构的第二个四分位数图像为2021年蓝色波段与红色波段的和值图像，记为RBsum图像；基于步骤S2中地面调查数据中包含的植被样本统计分析植被在RBsum图像中的分布区间，取其最大值作为阈值η，η＝0.18；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析其他不透水面在RBsum图像中的分布区间，取其最小值作为阈值σ，σ＝0.16。

S9、基于步骤S1中得到的N幅地表反射率图像中的近红外波段图像，采用四分位数图像重构方法对N幅近红外波段图像进行图像重构，取其重构的四个四分位数图像为2021年近红外波段最大值图像，记为NIRmax图像；基于步骤S2中地面调查数据中包含的植被样本统计分析植被在NIRmax图像中的分布区间，取其最小值作为阈值θ，θ＝0.3；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析其他不透水面在NIRmax图像中的分布区间，取其最大值作为阈值ο，ο＝0.45。

S10、创建红色不透水面遥感指数，基于步骤S1中得到的N幅地表反射率图像计算得到N幅红色不透水面遥感指数图像，然后采用四分位数图像重构方法对N幅红色不透水面遥感指数图像进行图像重构，取其重构的第三个四分位数图像为2021年红色不透水面遥感指数图像，记为RI图像；基于步骤S2中地面调查数据中包含的红色不透水面样本统计分析红色不透水面在RI图像中的分布区间，取其最小值作为阈值ι，ι＝0.09。

所述创建红色不透水面遥感指数的方法为：

RSI_i＝(ξ_red,i-ξ_green,i)/(ξ_red,i+ξ_green,i)；

其中，RSI_i表示红色不透水面遥感指数图像中像元i上的红色不透水面遥感指数值，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，ξ_green,i表示地表反射率图像中像元i上的绿色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

S11、基于步骤S1中得到的N幅地表反射率图像中的红色波段图像，采用四分位数图像重构方法对N幅红色波段图像进行图像重构，取其重构的二个四分位数图像为2021年红色波段中值图像，记为RM图像；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在RM图像中的分布区间，取其最小值作为阈值ξ，ξ＝0.13。

S12、判断像元i在步骤S3中得到的水体遥感指数均值图像中的像元值是否大于阈值α，若不是，执行步骤S13，若是，判断像元i在步骤S4中得到的第二短波红外波段最小值图像中的像元值是否小于阈值β，若是，像元i为水体，若不是，执行步骤S13。

S13、判断像元i在步骤S5中得到的蓝色彩钢板房遥感指数图像中的像元值是否大于阈值γ，若不是，执行步骤S14，若是，判断像元i在步骤S6中得到的蓝色波段中值图像中的像元值是否大于阈值δ，若是，像元i为蓝色彩钢板房，若不是，执行步骤S14。

S14、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否大于阈值ε，若不是，执行步骤S15，若是，像元i为植被。

S15、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否属于区间范围[ζ，ε]，若不是，执行步骤S16，若是，判断像元i在步骤S8中得到的蓝色波段与红色波段的和值图像中的像元值是否小于阈值η，若不是，执行步骤S16，若是，判断像元i在步骤S9中得到的近红外波段最大值图像中的像元值是否大于阈值θ，若不是，执行步骤S16，若是，像元i为植被。

S16、判断像元i在步骤S10中得到的红色不透水面遥感指数图像中的像元值是否大于阈值ι，若不是，执行步骤S17，若是，判断像元i在步骤S4中得到的第二短波红外波段中值图像中的像元值是否大于阈值κ，若不是，执行步骤S17，若是，像元i为红色不透水面。

S17、判断像元i在步骤S6中得到的蓝色波段中值图像中的像元值是否大于阈值μ，若不是，执行步骤S18，若是，判断像元i在步骤S4中得到的第二短波红外波段中值图像中的像元值是否大于阈值ν，若不是，执行步骤S18，若是，判断像元i在步骤S11中得到的红色波段中值图像中的像元值是否大于阈值ξ，若不是，执行步骤S18，若是，像元i为水泥不透水面。

S18、判断像元i在步骤S9中得到的近红外波段最大值图像中的像元值是否小于阈值ο，若不是，像元i为未识别类型，若是，执行步骤S19。

S19、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否小于阈值π，若是，像元i为其他不透水面，若不是，执行步骤S20。

S20、判断像元i在步骤S8中得到的蓝色波段与红色波段的和值图像中的像元值是否大于阈值σ，若是，像元i为其他不透水面，若不是，像元i为未识别类型。

S21、循环执行步骤S12至步骤S20，直至遍历完图像中的所有像元位置，完成研究区域的城市不同土地覆盖类型的遥感识别分类。

S22、采用3×3运算窗口对步骤S21中得到的分类结果进行滤波处理，消除分类结果中的孤立像元，完成最终的城市不同土地覆盖类型的遥感识别分类。

所述消除分类结果中孤立像元的方法为：对于像元(x，y)，图像中最大行数为X，最大列数为Y，其中x＝2，3，……，X-1；y＝2，3，……，Y-1；统计以像元(x，y)为中心的3×3窗口范围内不同类别的像元个数，设置像元(x，y)的数量权重为5，与像元(x，y)类别一致的像元个数为5+m，m为3×3窗口中与像元(x，y)类别一致的像元个数，令J表示为3×3窗口中与像元(x，y)类别不一致的类别中的最大像元个数，如果5+m小于J，则像元(x，y)的类别属性改写为与J类别一致的类别属性，否则，像元(x，y)的类别属性不变；依次遍历图像中所有可用像元位置，迭代3次，完成分类结果中孤立像元的消除。

为了验证本发明的效果，本发明以2021年开封市不同城市土地覆盖类型为实验对象，卫星影像为Landsat-8光学卫星影像，云平台为Google Earth Engine遥感云计算平台，识别结果如图2所示。通过图2可以看出城市土地覆盖类型的边界等纹理信息完整，说明了本发明对冰川识别的可靠性、准确性。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种城市不同土地覆盖类型遥感自动解译方法，其特征在于，其步骤如下：

S1、利用遥感云平台采集X年的N幅Landsat卫星影像，并对N幅Landsat卫星影像进行数据预处理和云层剔除处理，得到N幅地表反射率图像；其中，Landsat卫星影像包括蓝色波段、绿色波段、红色波段、近红外波段、第一短波红外波段和第二短波红外波段这6个波段影像；

S2、基于地面调查数据确定城市包含的土地覆盖类型，基于GPS定位工具获取土地覆盖类型的地理坐标，并使用光谱采集软件采集分析地表反射率图像中的土地覆盖类型的像元光谱样本，得到用于识别土地覆盖类型的有效光谱波段；其中，土地覆盖类型包括蓝色彩钢板房、红色不透水面、水泥不透水面、其他不透水面、植被、水体、未识别类型这7种类型；

S3、基于水体遥感指数和步骤S1中得到的N幅地表反射率图像计算得到N幅水体遥感指数图像，采用均值图像重构方法对N幅水体遥感指数图像进行图像重构，得到X年水体遥感指数均值图像，记为WI图像；基于步骤S2中地面调查数据中包含的水体样本统计分析水体在WI图像中的分布区间，取其最小值作为阈值α；

S4、基于步骤S1中得到的N幅地表反射率图像中的第二短波红外波段图像，采用四分位数图像重构方法对N幅第二短波红外波段图像进行图像重构，取其重构的第零个四分位数图像为X年第二短波红外波段最小值图像，记为SWIRmin图像，取其重构的第二个四分位数图像为X年第二短波红外波段中值图像，记为SWIRM图像；基于步骤S2中地面调查数据中包含的水体样本统计分析水体在SWIRmin图像中的分布区间，取其最大值作为阈值β；基于步骤S2中地面调查数据中包含的红色不透水面样本统计分析红色不透水面在SWIRM图像中的分布区间，取其最小值作为阈值κ；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在SWIRM图像中的分布区间，取其最小值作为阈值ν；

S5、创建蓝色彩钢板房遥感指数，基于步骤S1中得到的N幅地表反射率图像计算得到N幅蓝色彩钢板房遥感指数图像，然后采用四分位数图像重构方法对N幅蓝色彩钢板房遥感指数图像进行图像重构，取其重构的第三个四分位数图像为X年蓝色彩钢板房遥感指数图像，记为BI图像；基于步骤S2中地面调查数据中包含的蓝色彩钢板房样本统计分析蓝色彩钢板房在BI图像中的分布区间，取其最小值作为阈值γ；

S6、基于步骤S1中得到的N幅地表反射率图像中的蓝色波段图像，采用四分位数图像重构方法对N幅蓝色波段图像进行图像重构，取其重构的二个四分位数图像为X年蓝色波段中值图像，记为BM图像；基于步骤S2中地面调查数据中包含的蓝色彩钢板房样本统计分析蓝色彩钢板房在BM图像中的分布区间，取其最小值作为阈值δ；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在BM图像中的分布区间，取其最小值作为阈值μ；

S7、基于植被遥感指数和步骤S1中得到的N幅地表反射率图像计算得到N幅植被遥感指数图像，采用四分位数图像重构方法对N幅植被遥感指数图像进行图像重构，取其重构的第四个四分位数图像为X年植被遥感指数最大值图像，记为VI图像；基于步骤S2中地面调查数据中包含的植被样本统计分析纯净像元的植被在VI图像中的分布区间，取其最小值作为阈值ε；统计分析混合像元的植被在VI图像中的分布区间，取其最小值作为阈值ζ；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析纯净像元的其他不透水面在VI图像中的分布区间，取其最大值作为阈值π；

S8、基于步骤S1中得到的N幅地表反射率图像计算得到N幅蓝色波段与红色波段的和值图像，采用四分位数图像重构方法对N幅差值图像进行图像重构，取其重构的第二个四分位数图像为X年蓝色波段与红色波段的和值图像，记为RBsum图像；基于步骤S2中地面调查数据中包含的植被样本统计分析植被在RBsum图像中的分布区间，取其最大值作为阈值η；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析其他不透水面在RBsum图像中的分布区间，取其最小值作为阈值σ；

S9、基于步骤S1中得到的N幅地表反射率图像中的近红外波段图像，采用四分位数图像重构方法对N幅近红外波段图像进行图像重构，取其重构的四个四分位数图像为X年近红外波段最大值图像，记为NIRmax图像；基于步骤S2中地面调查数据中包含的植被样本统计分析植被在NIRmax图像中的分布区间，取其最小值作为阈值θ；基于步骤S2中地面调查数据中包含的其他不透水面样本统计分析其他不透水面在NIRmax图像中的分布区间，取其最大值作为阈值ο；

S10、创建红色不透水面遥感指数，基于步骤S1中得到的N幅地表反射率图像计算得到N幅红色不透水面遥感指数图像，然后采用四分位数图像重构方法对N幅红色不透水面遥感指数图像进行图像重构，取其重构的第三个四分位数图像为X年红色不透水面遥感指数图像，记为RI图像；基于步骤S2中地面调查数据中包含的红色不透水面样本统计分析红色不透水面在RI图像中的分布区间，取其最小值作为阈值ι；

S11、基于步骤S1中得到的N幅地表反射率图像中的红色波段图像，采用四分位数图像重构方法对N幅红色波段图像进行图像重构，取其重构的二个四分位数图像为X年红色波段中值图像，记为RM图像；基于步骤S2中地面调查数据中包含的水泥不透水面样本统计分析水泥不透水面在RM图像中的分布区间，取其最小值作为阈值ξ；

S12、判断像元i在步骤S3中得到的水体遥感指数均值图像中的像元值是否大于阈值α，若不是，执行步骤S13，若是，判断像元i在步骤S4中得到的第二短波红外波段最小值图像中的像元值是否小于阈值β，若是，像元i为水体，若不是，执行步骤S13；

S13、判断像元i在步骤S5中得到的蓝色彩钢板房遥感指数图像中的像元值是否大于阈值γ，若不是，执行步骤S14，若是，判断像元i在步骤S6中得到的蓝色波段中值图像中的像元值是否大于阈值δ，若是，像元i为蓝色彩钢板房，若不是，执行步骤S14；

S14、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否大于阈值ε，若不是，执行步骤S15，若是，像元i为植被；

S15、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否属于区间范围[ζ，ε]，若不是，执行步骤S16，若是，判断像元i在步骤S8中得到的蓝色波段与红色波段的和值图像中的像元值是否小于阈值η，若不是，执行步骤S16，若是，判断像元i在步骤S9中得到的近红外波段最大值图像中的像元值是否大于阈值θ，若不是，执行步骤S16，若是，像元i为植被；

S16、判断像元i在步骤S10中得到的红色不透水面遥感指数图像中的像元值是否大于阈值ι，若不是，执行步骤S17，若是，判断像元i在步骤S4中得到的第二短波红外波段中值图像中的像元值是否大于阈值κ，若不是，执行步骤S17，若是，像元i为红色不透水面；

S17、判断像元i在步骤S6中得到的蓝色波段中值图像中的像元值是否大于阈值μ，若不是，执行步骤S18，若是，判断像元i在步骤S4中得到的第二短波红外波段中值图像中的像元值是否大于阈值ν，若不是，执行步骤S18，若是，判断像元i在步骤S11中得到的红色波段中值图像中的像元值是否大于阈值ξ，若不是，执行步骤S18，若是，像元i为水泥不透水面；

S18、判断像元i在步骤S9中得到的近红外波段最大值图像中的像元值是否小于阈值ο，若不是，像元i为未识别类型，若是，执行步骤S19；

S19、判断像元i在步骤S7中得到的植被遥感指数最大值图像中的像元值是否小于阈值π，若是，像元i为其他不透水面，若不是，执行步骤S20；

S20、判断像元i在步骤S8中得到的蓝色波段与红色波段的和值图像中的像元值是否大于阈值σ，若是，像元i为其他不透水面，若不是，像元i为未识别类型；

S21、循环执行步骤S12至步骤S20，直至遍历完图像中的所有像元位置，完成研究区域的城市不同土地覆盖类型的遥感识别分类；

S22、采用3×3运算窗口对步骤S21中得到的分类结果进行滤波处理，消除分类结果中的孤立像元，完成最终的城市不同土地覆盖类型的遥感识别分类。

2.根据权利要求1所述的城市不同土地覆盖类型遥感自动解译方法，其特征在于，步骤S1中对Landsat卫星影像进行数据预处理和云层剔除处理的方法为：对Landsat卫星影像中所有像元值乘以修正系数0.0000275，再减去0.2，得到Landsat卫星影像的地表反射率图像；在地表反射率图像中像元i的位置上，如果蓝色波段、绿色波段、红色波段上像元值的极差的绝对值小于0.08，并且绿色波段上的像元值大于0.25，那么像元i的位置上的地物类型为云层，则将像元i的位置上所有波段的像元值改写为空值，实现云层剔除处理。

3.根据权利要求1或2所述的城市不同土地覆盖类型遥感自动解译方法，其特征在于，所述采用均值图像重构方法对N幅水体遥感指数图像进行图像重构的方法为：首先创建1幅像元行列数与水体遥感指数图像的像元行列数一致的空值图像；然后将N幅水体遥感指数图像进行图层堆叠，像元i的位置上存在n个水体遥感指数值，计算得到这n个数值的平均值，并写入空值图像的相应位置，依次遍历所有像元位置，得到水体遥感指数均值图像。

4.根据权利要求1所述的城市不同土地覆盖类型遥感自动解译方法，其特征在于，所述采用四分位数图像重构方法对N幅第二短波红外波段图像进行图像重构的方法为：首先创建5幅像元行列数与第二短波红外波段图像的像元行列数一致的空值图像；然后将N幅第二短波红外波段图像进行图层堆叠，像元i的位置上存在n个第二短波红外波段像元值，对这n个值按照从小到大的顺序排列成数列，分别将数列中位于第零、一、二、三、四的四分位点的值写入空值图像中的相应像元位置，依次遍历所有像元位置，即分别得到第二短波红外波段图像的第零、一、二、三、四的四分位数图像；其中第零、二、四的四分位点分别表示数组的最小值、中值、最大值。

5.根据权利要求1所述的城市不同土地覆盖类型遥感自动解译方法，其特征在于，所述步骤S5中创建蓝色彩钢板房遥感指数的方法为：

BPI_i＝(ξ_blue,i-ξ_red,i)/(ξ_blue,i+ξ_red,i)；

其中，BPI_i表示蓝色彩钢板房遥感指数图像中像元i上的蓝色彩钢板房遥感指数值，ξ_blue,i表示地表反射率图像中像元i上的蓝色波段的地表反射率，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

6.根据权利要求1或5所述的城市不同土地覆盖类型遥感自动解译方法，其特征在于，所述步骤S8中创建红色不透水面遥感指数的方法为：

RRI_i＝(ξ_red,i-ξ_green,i)/(ξ_red,i+ξ_green,i)；

其中，RRI_i表示红色不透水面遥感指数图像中像元i上的红色不透水面遥感指数值，ξ_red,i表示地表反射率图像中像元i上的红色波段的地表反射率，ξ_green,i表示地表反射率图像中像元i上的绿色波段的地表反射率，i＝1,2,…,n，n为地表反射率图像中像元的总数量。

7.根据权利要求1所述的城市不同土地覆盖类型遥感自动解译方法，其特征在于，所述步骤S22中消除分类结果中孤立像元的方法为：对于像元(x，y)，图像中最大行数为X，最大列数为Y，其中x＝2，3，……，X-1；y＝2，3，……，Y-1；统计以像元(x，y)为中心的3×3窗口范围内不同类别的像元个数，设置像元(x，y)的数量权重为5，与像元(x，y)类别一致的像元个数为5+m，m为3×3窗口中与像元(x，y)类别一致的像元个数，令J表示为3×3窗口中与像元(x，y)类别不一致的类别中的最大像元个数，如果5+m小于J，则像元(x，y)的类别属性改写为与J类别一致的类别属性，否则，像元(x，y)的类别属性不变；依次遍历图像中所有可用像元位置，迭代3次，完成分类结果中孤立像元的消除。

Images