CN109919875A - 一种高时频遥感图像特征辅助的居民地提取与分类方法 - Google Patents

Abstract

本发明公开一种高时频遥感图像特征辅助的居民地提取与分类方法，其步骤：首先收集覆盖研究区域的高时频遥感图像和相应区域的夜间灯光数据，其次对收集的高时频遥感图像进行辐射定标、几何校正和图像配准，对夜间灯光数据进行图像重采样，然后采用NDVI最大值合成阈值分割、NDWI最大值合成阈值分割、近红外波段辐亮度随太阳高度角的变化斜率阈值分割、夜间灯光数据约束等逐步剔除背景地类实现居民地提取，得到居民地提取结果图，最后通过居民地连续斑块的面积和夜间灯光的阈值分割实现居民地分类，得到居民地分类结果图。方法易行，操作简便，能够较好的将水体、裸地等与居民地区分开来，且能够实现居民地二级分类，居民地提取与分类精度较高。

Description

一种高时频遥感图像特征辅助的居民地提取与分类方法

技术领域

本发明属于遥感图像分类识别技术领域，更具体地涉及一种利用高时频遥感图像数据和夜间灯光数据进行居民地提取与分类的方法，适用于可以获取同一区域地物类型未发生变化的时间段内多个不同太阳高度角的多光谱遥感图像序列及该地区的夜间灯光数据的情况下的居民地调查与监测。

背景技术

居民地是人与自然相互作用的产物，是人类因生产和生活需要而集聚定居的场所。近年来，随着经济社会发展，我国居民地迅速扩张(王雷,李丛丛,应清,程晓,王晓昳,李雪艳,胡娈运,梁璐,俞乐,黄华兵,宫鹏.中国1990～2010年城市扩张卫星遥感制图[J].科学通报,2012,57(16):1388-1399.)，不仅改变了地球表面覆盖与形态，而且直接影响着局部、区域乃至全球的气候、生物化学、水文过程。快速、准确、全面、客观地获取与分析居民地的空间分布和时空变化，不仅能为人口、经济等统计数据的空间化提供基本的骨架，在居民地管理、耕地保护、生态环境动态监测与保护、地理国(世)情监测、灾害评估等方面也具有重要意义(陈军,陈利军,李然,廖安平,彭舒,鲁楠,张宇硕.基于GlobeLand30的全球城乡建设用地空间分布与变化统计分析[J].测绘学报,2015,11:1181-1188.)。新型城镇化及智慧城市建设的持续推进，给更精确的、更及时的监测居民地分布和变化提出了更高的要求。遥感技术具有快速、宏观、综合、准确、周期性、成本低等优势，已经成为居民地提取与监测的重要手段。

现阶段，利用遥感技术提取居民地信息的方法主要可以分为目视解译方法、遥感模型方法等。目视解译方法由于加入了人的知识，提取精度较高(王雷,李丛丛,应清,程晓,王晓昳,李雪艳,胡娈运,梁璐,俞乐,黄华兵,宫鹏.中国1990～2010年城市扩张卫星遥感制图[J].科学通报,2012,57(16):1388-1399.)，但工作较繁琐，费时费力，同时需要解译人员有丰富的遥感解译知识与经验。遥感模型方法是常用的居民地信息提取方法，该类方法就是根据一定的定量运算方法从影像上获得地物的空间分布。基于遥感模型的居民地信息提取方法大致可以分为四类：基于统计分类的模型、基于决策树的模型、基于知识发现的模型和面向对象的模型(顾娟,陈军,张宏伟,周启鸣.基于遥感影像的居民地提取研究综述与展望[J].计算机工程与应用,2007(30):1-4+10.)。

(1)基于统计的分类模型：通过提取各类训练样本确定判别函数和模型参数，从而把各像元划归到各个给定类(Gong P,Wang J,Le Y.Finer resolution observation andmonitoring of global land cover:first mapping results with Landsat TM and ETM+data[J].Int.J.Remote Sens.,2013,34(7):2607-2654.)。统计分类模型可以通过统计模型对居民地在任意维数特征空间中的分布进行建模，不同影像特征可以互为补充，因此被认为是一种稳定性高、鲁棒性好的方法。但目前大部分研究采用的是传统的高斯分布模型描述地物的概率分布，高斯分布主要描述对称的、单峰的数据集合。如果居民地在影像特征空间中的分布比较复杂，不符合高斯分布时，采用简单的单峰形式高斯分布密度函数来表征居民地在特征空间的分布，就可能造成与实际分布的较大偏差，导致分类精度下降。

(2)基于决策树的模型：通过一些判断条件对原始数据集逐步进行二分和细化(李振敏,王晓青,窦爱霞,杨海霞,黄树松,崔丽萍.基于Landsat-8与ZY3多源遥感影像的城镇居民地识别与再分类研究[J].地震,2015,03:123-135.)。决策树方法的优点是不需要依赖任何先验假设条件，而且可以方便地利用除亮度值以外的其它知识，所以在遥感影像分类和专题信息提取中有广泛的应用。当问题比较简单以及训练样本很少时，决策树方法十分有效。决策树方法的关键在于阈值的取定，常用的有样本观测、经验知识及基于信息熵等方法。

(3)基于知识发现的模型：利用地物的光谱特征、空间结构与形态、地物空间关系以及动态变化等知识所建立的遥感专题信息提取模型(杨智翔,何秀凤.基于改进的NDBI指数法的遥感影像城镇用地信息自动提取[J].河海大学学报(自然科学版),2010,38(2):181-184.)。如杨山(杨山.发达地区城乡聚落形态的信息提取与分形研究——以无锡市为例[J].地理学报,2000(06):671-678.)在深入研究NDVI(Normalized DifferenceVegetation Index)的基础上，发现利用TM4(近红外)和TM5(短波红外)两个波段间的运算，较好地提取城镇居民地。基于该思想查勇等(查勇,倪绍祥,杨山.一种利用TM图像自动提取城镇用地信息的有效方法[J].遥感学报,2003(01):37-40+82.)提出了归一化建筑指数(Normalized Difference Built-up Index，NDBI)，NDBI＝(TM5-TM4)/(TM5+TM4)，并利用NDBI自动提取了无锡市的城镇用地信息,取得了较好的结果。NDBI是居民地提取中常用的特征指数，当然基于NDBI的居民地提取结果中仍难以避免的混分了低密度植被覆盖区和裸地。同时，目前多数的卫星传感器没有设置短波红外波段，使得基于NDBI的居民地提取方法难以推广应用。基于知识发现的方法，在一定程度上可提高计算机解译精度，但远未达到实用阶段。

(4)面向对象的模型：在分类时不仅利用影像的光谱特征，更多的是利用了纹理、形状及空间关系等特征(陈军,陈利军,李然,廖安平,彭舒,鲁楠,张宇硕.基于GlobeLand30的全球城乡建设用地空间分布与变化统计分析[J].测绘学报,2015,11:1181-1188.)。随着高分辨率卫星图像的深入应用，现有的居民地提取中图像的光谱信息以及尺寸、形状、纹理、邻近关系等空间细节信息得到了较为充分地利用(Bagan H,Yamagata Y.Landsatanalysis of urban growth:how Tokyo became the world's largest megacity duringthe last 40years[J].Remote Sensing of Environment,2012,127:210–222；金飞.基于纹理特征的遥感影像居民地提取技术研究,解放军信息工程大学,摄影测量与遥感,2013,博士；Sun C,Wu Z,Lv Z,Yao N,Wei J.Quantifying different types of urban growthand the change dynamic in Guangzhou using multi-temporal remote sensing data[J].International Journal of Applied Earth Observation and Geoinformation,2013,21:409-417；Wania A,Kemper T,Tiede D,Zeil P.Mapping recent built-up areachanges in the city of Harare with high resolution satellite imagery[J].Applied Geography,2014,46:35-44.)，从一定程度上提高了信息提取的精度。但是，居民地的信息提取精度随遥感影像空间分辨率的提高是有限的，原因主要在于空间分辨率的提高会放大居民地内部结构的细节，增加了一些干扰性的非目标噪声(汪权方,许纪承,陈媛媛,李家永,王新生.遥感影像空间分辨率对居民地信息提取的影响[J].资源科学,2012,34(1):159-165.)。

虽然目前已经形成了一系列基于遥感的居民地提取模型方法，但居民地提取精度问题在某些背景地物反射率较低或较高的地区仍然比较突出，如水体与低亮度居民地有相似的反射特征，裸地与高亮度居民地有相似的反射特征，目前的居民地提取方法容易将水体、裸地等与居民地混淆，从而导致居民地的提取精度下降，且常用的基于NDBI方法因波段限制难以适用于仅具有可见光-近红外波段的传感器影像。

另外，当前居民地提取使用的数据基本为中高分辨率的太阳同步轨道卫星图像，如Landsat系列卫星图像，这些卫星的重访周期较长，且通常在一天中的固定时间拍摄，难以获取短期(地物类型未发生变化的时间段内，如一日内)不同太阳高度角的图像，因此无法获取高时频信息。高分四号(GF-4)卫星是我国第一颗地球同步轨道遥感卫星，具备高时间分辨率、中等空间分辨率观测能力，配置地球同步轨道50m分辨率的凝视相机，可以完成分钟级别的遥感监测任务，能够获得同一区域短时内不同太阳高度角的图像，其时间维特征为基于遥感图像的地物识别和提取提供了新的思路。

夜间灯光数据作为人类活动的一种有效表征形式，在宏观尺度的城市研究中具有巨大的潜力和应用前景。基于夜间灯光数据通过阈值法可以实现城区提取并进行城市扩张监测(罗庆,李小建.基于VIIRS夜间灯光的中国城市中心的分异特征及其影响因素[J].地理研究,2019(01):155-166；刘佳,辛鑫,刘斌,邸凯昌,岳宗玉,王承安.基于DMSP/OLS夜间灯光影像的2000—2013年鄂尔多斯市城市扩张遥感制图与驱动因子分析[J].国土资源遥感,2018,30(01):166-172.)。目前的夜光遥感对地观测平台主要是美国国家海洋和大气管理局NOAA(National Oceanic and Atmospheric Administration)的DMSP/OLS(DefenseMeteorological Satellite Program-Operational Linescan System)(1992年-2013年)和NPP/VIIRS(National Polar-Orbiting Partnership’s Visible Infrared ImagingRadiometer Suite)(2012年-至今)。与1km分辨率的DMSP/OLS影像相比，NPP/VIIRS分辨率更高(500m)，较宽的辐射探测范围解决了DMSP/OLS数据的像元过饱和问题。但是由于“灯光溢出”效应的存在，夜间灯光数据探测到的反映人类活动的灯光分布区范围往往要大于实际情况，难以确定一个城市区域提取的最佳灯光强度阈值实现城市居民地的精确提取(赵敏,程维明.基于DMSP/OLS夜间灯光数据的城市空间扩展研究综述[J].测绘与空间地理信息,2015,38(03):64-68.)。然而，OLS与VIIRS能探测到城市灯光甚至小规模的居民地、车流等发出的低强度灯光，明显区别于黑暗的乡村背景，从而有利于城镇与农村居民点的区分。综合利用中高分辨率光学卫星遥感图像和夜间灯光数据的优势，可以提高居民地提取精度，并实现居民地二级分类。

针对目前基于遥感图像的居民地提取中存在的问题，结合GF-4卫星的高时频数据获取能力和夜间灯光数据的城镇与乡村区分能力，本发明提出了一种利用高时频遥感图像特征结合夜间灯光数据的居民地提取与分类方法，提高基于遥感的居民地信息提取与分类精度。

发明内容

为解决如何实现基于遥感图像的居民地信息准确提取的技术问题，本发明的目的是在于提供了一种利用高时频遥感图像数据和夜间灯光数据的居民地提取与分类方法，方法易行，操作简便，能够较好地将水体、裸地等易混淆地类与居民地区分开来，并实现了居民地二级分类。

为了实现上述的目的，本发明采用以下技术方案：

本发明技术构思是：首先是收集覆盖研究区域的不同太阳高度角的高时频多光谱遥感图像和夜间灯光数据；其次是对收集的高时频多光谱遥感图像进行几何校正、图像配准和辐射定标等预处理，对收集的夜间灯光数据进行图像重采样；再次是基于上述预处理后的高时频多光谱遥感图像，分别计算NDVI、NDWI最大值合成图，通过阈值分割剔除植被和大部分水体；基于上述预处理后的高时频多光谱遥感图像，采用趋势分析方法得到近红外波段辐亮度随太阳高度角的变化斜率图，通过阈值分割剔除易与居民地混淆的水体；基于上述预处理后的高时频多光谱遥感图像，计算图像的绿、红、近红外波段的辐亮度之和的最大值合成图，对其设置低亮阈值,剔除水体、阴影等低亮地表，获得宽约束下的居民地提取结果；在宽约束下的居民地提取结果基础上，对绿、红、近红外波段的辐亮度之和的最大值合成图设置高亮阈值,剔除云、雪、裸地等高亮地表(其中可能包含机场等亮度较高的居民地地表)，获得严约束下的居民地提取结果；基于宽约束下的居民地提取结果，对夜间灯光数据设置阈值,提取高亮居民地地表，并将其与上一步严约束下的居民地提取结果合并，获得完整的居民地提取结果；最后，利用上一步的居民地提取结果，对居民地连续斑块的面积和夜间灯光图像设置阈值，区分城镇与农村居民点，从而实现居民地的二级分类。

一种高时频遥感图像特征辅助的居民地提取与分类方法，其步骤是：

1.数据准备。收集覆盖研究区域的不同太阳高度角的高时频多光谱遥感图像和夜间灯光图像。

所述的研究区域高时频多光谱遥感图像收集的方法，其步骤是：

(1)每景遥感图像均覆盖研究区域，且每景图像的云覆盖区域均不大于10％；

(2)每景遥感图像均包含绿、红、近红外波段，最好是同一传感器拍摄的图像；

(3)收集不同时相遥感图像n景(n>＝2)，不同时相图像上的地物类型未发生变化，一般数据的时间跨度不超过10天；

(4)遥感图像拍摄时的太阳高度角在30－75度，且最大太阳高度角和最小太阳高度角之间的差异不小于10度；

所述的研究区域夜间灯光图像收集的方法，其步骤是：

(1)采用NPP/VIIRS DNB(Day/Night Band)年合成夜间灯光图像数据，数据采集时间与高时频多光谱遥感图像采集时间尽量接近；

2.图像预处理。对步骤1收集的高时频多光谱遥感图像进行辐射定标、几何校正和图像配准，对夜间灯光图像进行图像重采样处理。

所述的高时频遥感图像辐射定标、几何校正和图像配准的方法，其步骤是：

(1)辐射定标。通过从卫星图像的头文件中读取定标系数，使用以下计算公式将各时相影像的数字量化值DN(Digital Number)转换成辐亮度，实现图像的辐射定标：

L_t＝Gain×DN+Bias

式中，L_t为图像的辐亮度，Gain为增益，Bias为偏置，单位均为W﹒m^-2﹒sr^-1﹒μm^-1，DN为卫星载荷观测值，无量纲。

(2)几何校正。利用卫星图像自带的有理多项式系数(Rational PolynomialCoefficients，RPC)，通过有理函数模型(Rational Function Model，RFM)对图像进行几何校正。

(3)图像配准。选定一个时相的遥感图像为基准图像，其它(n-1)个时相的遥感图像作为待配准图像，通过选取基准图像与待配准图像上的同名特征点，然后利用这些同名特征点建立基准图像与待配准图像之间的几何畸变模型(如多项式模型)，再利用该几何畸变模型进行几何校正，从而实现不同时相遥感图像的图像配准。

所述的夜间灯光图像进行图像重采样的方法，其步骤是：

(1)图像重采样。通过图像灰度值重采样对夜间灯光图像进行空间分辨率变换，使得变换后的图像空间分辨率与高时频多光谱遥感图像的空间分辨率一致。采用数学上的双向线性内插法作为重采样方法计算分辨率变换后像元位置的灰度值。

双线性插值法是将与待求像元的四个邻近像元的灰度值在行和列两个方向上作线性内插。通常根据待求像元与邻近的四个像元的距离权重计算出灰度值，并赋给待求像元。双线性插值法公式为：

式中，g_x'y'为输出像元灰度值，g_i为邻近点i的灰度值，p_i为邻近点对投影点的权重(p_i＝1/d_i，d_i表示邻近点到投影点的距离)；g₁、g₂、g₃、g₄分别为四个邻近点的灰度值，p₁、p₂、p₃、p₄分别为四个邻近点对投影点的权重。

3.居民地提取。基于上述(步骤2)预处理后的高时频遥感图像和夜间灯光图像，通过逐步剔除背景地类实现居民地提取。

所述的逐步剔除背景地类实现居民地提取的方法，其步骤是：

(1)NDVI最大值合成阈值分割。基于上述(步骤2)图像预处理后的高时频多光谱遥感图像，根据NDVI计算公式计算每个时相的NDVI，并对各时相的NDVI进行最大值合成，得到NDVI最大值合成图。基于NDVI最大值合成图，通过阈值分割剔除植被覆盖地表。

根据NDVI计算公式计算每个时相的NDVI,包括以下步骤：

NDVI为归一化植被指数(Normalized Difference Vegetation Index)，其计算公式为：

其中，NDVI为归一化植被指数(无量纲)，NIR和R分别为图像的近红外波段和红波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1。

对各时相的NDVI进行最大值合成，得到NDVI最大值合成图，结果图中每个像元的值均为相应像元在各时相NDVI中的最大数值,包括以下步骤：

采用最大值合成法(Maximum Value Composites,MVC)进行NDVI最大值合成，即选取观测像元的在时间序列图像中的最大NDVI值作为该像元合成期的NDVI。MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数。

基于NDVI最大值合成图，通过阈值分割剔除植被覆盖地表，包括以下步骤：

1)对于NDVI最大值合成图，设定阈值T1，T1为正数；

2)对于NDVI最大值合成图中像元值大于T1的像元，将其归入植被覆盖地表类别(非居民地)，赋值为0，其他像元归入备选居民地类别，赋值为1。由于植被覆盖地表的NDVI通常大于居民地，通过设置合适的阈值，可以将植被覆盖地表从图像中剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域。

(2)NDWI最大值合成阈值分割。基于上述(步骤2)图像预处理后的高时频多光谱遥感图像，根据NDWI计算公式计算每个时相的NDWI，并对各时相的NDWI进行最大值合成，得到NDWI最大值合成图。基于NDWI最大值合成图，通过阈值分割剔除水体。

根据NDWI计算公式计算每个时相的NDWI,包括以下步骤：

NDWI为归一化水体指数(Normalized Difference Water Index)，其计算公式为：

其中，NDWI为归一化水体指数(无量纲)，NIR和G分别为图像的近红外波段和绿波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1。

对各时相的NDWI进行最大值合成，得到NDWI最大值合成图，结果图中每个像元的值均为相应像元在各时相NDWI中的最大数值,包括以下步骤：

采用最大值合成法(Maximum Value Composites,MVC)进行NDWI最大值合成，即选取观测像元的在时间序列图像中的最大NDWI值作为该像元合成期的NDWI。MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数。

基于NDWI最大值合成图，通过阈值分割剔除水体，包括以下步骤：

1)对于NDWI最大值合成图，设定阈值T2，T2为正数；

2)针对NDVI最大值合成阈值分割结果图中像元值为1的区域(即备选居民地区域)，将NDWI最大值合成图中像元值小于T2的像元归入备选居民地区域，赋值为1，其他像元归入水体(非居民地)，赋值为0。由于水体的NDWI通常大于居民地，通过设置合适的阈值，可以将水体从图像中剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域。

(3)近红外波段辐亮度变化斜率阈值分割。基于上述(步骤2)图像预处理后的高时频多光谱遥感图像，采用趋势分析方法分析近红外波段辐亮度随太阳高度角变化的趋势特征，得到近红外波段辐亮度随太阳高度角的变化斜率图。对近红外波段辐亮度随太阳高度角的变化斜率图进行阈值分割剔除易与居民地混淆的水体。

采用趋势分析方法获取近红外波段辐亮度随太阳高度角的变化斜率图,包括以下步骤：

采用Theil-Sen Median趋势分析方法，计算n(n-1)/2个数据组合的斜率的中位数。该方法有由Henri Theil于1950年、Pranab K.Sen于1968年分别提出的研究长时序变化的方法(Sen,Pranab Kumar.Estimates of the regression coefficient based onKendall's tau.Journal of the American Statistical Association,1968,63(324):1379–1389.),其公式为：

其中，slope为图像像元近红外波段辐亮度值随太阳高度角的变化斜率(无量纲)，NIR_i和NIR_j分别为i时相和j时相的近红外波段辐亮度值(无量纲)，Sun_i和Sun_j分别为i时相和j时相的太阳高度角值(无量纲)。median是一种计算机函数，返回一组数值的中值，即将这组数值按照从小到大的顺序排列，返回居于中间的数值。如果这组数值集合中包含偶数个数值，则返回位于中间的两个数的平均值。

按照以上公式计算各像元近红外波段辐射亮度值随太阳高度角的变化斜率，生成近红外波段辐亮度变化斜率图。

对近红外波段辐亮度变化斜率图设定阈值进行水体剔除的方法，其步骤是：

1)对于生成的近红外波段辐亮度变化斜率图，设定斜率阈值T3，T3为正数；

2)针对NDWI最大值合成阈值分割结果图中像元值为1的区域(即备选居民地区域)，将斜率slope大于T3的像元归入备选居民地区域，赋值为1，其他像元归入水体(非居民地)，赋值为0。由于水体像元的斜率slope通常小于居民地，通过设置合适的阈值，可以将低亮居民地与水体区分。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域。

(4)绿、红、近红外波段辐亮度之和最大值合成。基于上述(步骤2)图像预处理后的高时频多光谱遥感图像，计算每个时相图像的绿、红、近红外波段辐亮度之和，并对各时相的辐亮度和进行最大值合成，得到辐亮度和的最大值合成图。

计算每个时相图像的绿、红、近红外波段辐亮度之和包括以下步骤：

绿、红、近红外波段辐亮度之和计算公式为：

MSS＝G+R+NIR

其中，MSS为绿、红、近红外波段辐亮度之和，G、R和NIR分别为图像的绿、红、近红外波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1。

对各时相的辐亮度和进行最大值合成，得到辐亮度和的最大值合成图,结果图中每个像元的值均为相应像元在各时相辐亮度和中的最大数值,包括以下步骤：

采用最大值合成法(Maximum Value Composites,MVC)进行辐亮度和最大值合成，即选取观测像元在时间序列图像中的最大辐亮度和(即像元在各时相图像中的绿、红、近红外波段辐亮度之和的最大值)作为该像元合成期的辐亮度和。MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数。

(5)宽约束下的居民地初提取，设置阈值在保证居民地完整的前提下(即尽量减小漏检误差)，尽量剔除背景地类。基于绿、红、近红外波段辐亮度之和最大值合成结果图，通过设置阈值剔除辐亮度较低的阴影、水体，得到宽约束下的居民地初提取结果。包括以下步骤：

1)对于绿、红、近红外波段辐亮度之和最大值合成结果图，设定阈值T4，T4为正数，T4的设置在保证居民地完整的前提下，尽量剔除阴影、水体等混淆地物；

2)针对近红外波段辐亮度变化斜率阈值分割结果图中像元值为1的区域(即备选居民地区域)，将辐亮度和最大值合成图中像元值大于T4的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0。由于水体、阴影的辐亮度和通常小于居民地，通过设置合适的阈值，可以在保证居民地完整的前提下，将这些非居民地区域从图像中尽量剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为宽约束下的居民地区域。宽约束下的居民地初提取结果中包含了较完整的居民地，但仍包含部分错提的辐亮度较高的云、雪、裸地等像元。

(6)严约束下的居民地初提取，设置阈值在保证较彻底剔除背景地类的前提下(即尽量减小错检误差)，尽量保留较完整的居民地。基于绿、红、近红外波段辐亮度之和最大值合成结果图，通过设置阈值剔除辐亮度较高的云、雪、裸地等，得到严约束下的居民地初提取结果。包括以下步骤：

1)对于绿、红、近红外波段辐亮度之和最大值合成结果图，设定阈值T5，T5为正数，T5的设置在保证较彻底剔除云、雪、裸地等非居民地的前提下，尽量保留较完整的居民地；

2)针对宽约束下的居民地初提取结果图中像元值为1的区域(即备选居民地区域)，将辐亮度和最大值合成图中像元值小于T5的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0。由于云、雪、裸地的辐亮度和通常大于居民地，通过设置辐亮度和大于居民地而小于云、雪、裸地的阈值，可以将这些非居民地区域从图像中较彻底的剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为严约束下的居民地区域。严约束下的居民地初提取结果中较彻底的剔除了云、雪、裸地等非居民地噪声，但居民地区域提取不够完整，一些辐亮度较高的城镇区域容易漏提取。

(7)夜间灯光严约束下的城镇提取。基于上述(步骤2)图像预处理后的夜间灯光数据，通过设置较严格的阈值得到基于灯光指数的城镇区域提取结果，并将其与宽约束下的居民地初提取结果结合，得到夜间灯光严约束下的城镇提取结果。包括以下步骤：

1)对于图像预处理后的夜间灯光数据，设定阈值T6，T6为正数，T6的设置在保证较彻底剔除非城镇区域的前提下，尽量保证城镇的完整；

2)针对宽约束下的居民地初提取结果图中像元值为1的区域(宽约束下的居民地区域)，将夜间灯光指数大于T6的像元归入城镇区域，赋值为1，其他像元归入非居民地，赋值为0。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为夜间灯光严约束下的城镇区域。

(8)结果合并。将严约束下的居民地初提取结果与夜间灯光严约束下的城镇提取结果进行合并，得到完整且噪声尽量小的最终的居民地提取结果。包括以下步骤：

针对严约束下的居民地初提取结果图中像元值为1的区域(严约束下的居民地区域)和夜间灯光严约束下的城镇提取结果图中像元值为1的区域(夜间灯光严约束下的城镇区域)归入居民地区域，赋值为1，其他像元归入非居民地，赋值为0。

居民地提取结果图中像元值为0的区域为居民地区域，像元值为1的区域为居民地区域。

4.居民地分类。基于上述(步骤3)居民地提取结果，通过居民地连续斑块的面积和夜间灯光的阈值分割实现居民地分类。

(1)连续斑块面积、夜间灯光阈值分割。基于上述居民地提取(步骤3)结果图，通过设置居民地图斑连续斑块的面积阈值和夜间灯光阈值，将居民地中的城镇与农村居民点区分开来。包括以下步骤：

1)对于居民地提取结果图，设定连续斑块面积阈值T7，T7为正数，T7设置为保证小于最小的城镇图斑面积的前提下，尽量大于农村居民点图斑面积；

2)对于图像预处理后的夜间灯光数据，设定阈值T8，T8为正数，T8的设置在保证城镇完成的前提下，尽量剔除非城镇区域；

3)对居民地提取结果图中连续斑块的面积大于阈值T7的居民地图斑，将平均夜间灯光指数大于T8的图斑归入城镇，图斑内的像元赋值为2；对居民地提取结果图中连续斑块的面积大于阈值T7的居民地图斑，将平均夜间灯光指数小于T8的图斑归入农村居民点，图斑内的像元赋值为1；对居民地提取结果图中连续斑块的面积小于等于阈值T7的居民地图斑，归入农村居民点，图斑内的像元赋值为1；其他像元归入非居民地，赋值为0。

(2)居民地分类结果居民地分类结果图中像元值为0的区域为非居民地区域，像元值为1的区域为农村居民点，像元值为2的区域为城镇。

上述技术措施，填补了目前基于高时频遥感图像进行居民地提取的空白，通过利用居民地与其他地类在光谱与时序变化特征上的差异实现了基于高时频遥感图像的居民地高精度提取，同时结合夜间灯光数据区分城镇与乡村背景的优势，实现了居民地的二级分类。

关键步骤单元居民地提取是基于高时频遥感图像和夜间灯光图像，通过逐步剔除背景地类实现居民地提取。

解决的技术难题：解决了居民地的高精度自动提取和二级分类问题，尤其是居民地提取中容易混入低反射的水体和高反射的裸地的问题。

主要解决的技术问题：如何将背景地类(尤其是易与居民地混淆的地类)逐步剔除，提高居民地提取的精度。

达到的技术效果：能够实现居民地的高精度提取，对于基于光谱的居民地提取方法中容易与居民地混淆的水体、裸地能够较好的区分，同时能够实现居民地的二级分类。

发明相对于现有技术的进步：在居民地提取中引入了地物反射光谱随时间(太阳高度角)的变化特征信息，并引入了夜间灯光数据的约束，提高了居民地提取的精度，尤其是减少了由水体、裸地造成的错检误差。

技术方案与现有技术的主要区别：现有的遥感图像居民地提取技术一般是基于单时相的遥感图像，仅利用了遥感图像的光谱及纹理等信息，基本没有引入时序特征信息，本发明通过引入地物反射光谱随时间(太阳高度角)的变化特征信息实现居民地的遥感图像提取，为居民地提取开辟了新的途径。

本发明与现有技术相比，具有以下有益效果及优点：

1、填补了目前基于高时频遥感图像进行居民地信息提取的空白，通过利用居民地与其他地类在近红外波段的反射光谱随太阳高度角的变化特征差异实现了基于高时频遥感图像的居民地提取。

2、与目前常用的监督分类法等相比，本发明通过引入地物时间维的变化特征实现居民地的提取，并结合夜间灯光数据的优势，能够得到更高的提取精度，能够较好的将水体、裸地与居民地区分开来，且能够很好的实现居民地二级分类。

本发明在我国三个100km*100km范围的不同地区的GF-4卫星高时频影像上分别进行了居民地提取实验，居民地的漏检误差均不超过10％，错检误差均不超过6％，证明了该发明的有效性。

附图说明

图1为一种高时频遥感图像特征辅助的居民地提取与分类方法的流程图；

图2为一种所采用的包含居民地的GF-4高时频遥感图像；

a)获取时间为2016年8月28日9:09:00，太阳高度角39.75度；b)获取时间为2016年8月27日9:43:38，太阳高度角43.56度；c)获取时间为2016年8月27日10:54:12，太阳高度角55.64度。

图3为一种所采用的与GF-4高时频遥感图像区域相同的夜间灯光数据图像；

图4为一种高时频遥感图像特征辅助的居民地提取与分类方法的居民地提取结果图；

图5为一种高时频遥感图像特征辅助的居民地提取与分类方法的居民地分类结果图；

图6为一种居民地提取结果与常规方法居民地提取结果的比较；

a)SVM监督分类方法居民地提取结果；b)本发明方法居民地提取结果。

SVM监督分类方法的结果对训练区的选择具有较强的依赖性，且分类速度较慢，耗时长，易将水体、裸地与居民地混淆。本发明方法速度快，能够较好的将水体、裸地等与居民地区分开来，有效减少错检误差。

图7为一种居民地提取与常规方法的精度定量比较。

相比于基于光谱的SVM监督分类居民地提取方法，本发明方法具有更高的居民地提取精度，其漏检误差和错检误差均有所降低，尤其错检误差显著降低。

具体实施方式

实施例1：

下面结合附图通过实施例对本发明做进一步的详细说明。

根据图1、图2、图3、图4、图5、图6、图7可知，一种高时频遥感图像特征辅助的居民地提取与分类方法，具体步骤包括：

A、单元100数据准备。收集覆盖研究区域的不同太阳高度角的高时频多光谱遥感图像(如图2)和夜间灯光图像(如图3)。

研究区域高时频多光谱遥感图像收集的方法，其步骤是：

(1)每景遥感图像均覆盖研究区域，且每景图像的云覆盖区域均不大于10％；

(2)每景遥感图像均包含绿、红、近红外波段，最好是同一传感器拍摄的图像，本例中所用的遥感图像均为GF-4多光谱传感器图像；；

(3)收集不同时相遥感图像n景(n>＝2)，不同时相图像上的地物类型未发生变化，一般数据的时间跨度不超过10天，本例中n＝3，所用图像均为同一天拍摄；

(4)遥感图像拍摄时的太阳高度角在30-75度，且最大太阳高度角和最小太阳高度角之间的差异不小于10度，本例中遥感图像拍摄时的太阳高度角分别是39.75度、43.56度和55.64度。

研究区域夜间灯光图像收集的方法，其步骤是：

(1)采用NPP/VIIRS DNB(Day/Night Band)年合成夜间灯光图像数据，数据采集时间与高时频多光谱遥感图像采集时间尽量接近，本例采用2016年NPP/VIIRS DNB年合成夜间灯光数据。夜间灯光数据来源于美国国家海洋和大气管理局NOAA(National Oceanicand Atmospheric Administration)下属的国家地球物理数据中心NGDC(NationalGeophysical Data Center,https://ngdc.noaa.gov/eog/viirs/download_dnb_composites.html)。

B、单元101图像预处理。对收集的高时频多光谱遥感图像进行辐射定标、几何校正和图像配准，对夜间灯光图像进行图像重采样处理。步骤如下：

高时频遥感图像辐射定标、几何校正和图像配准包括以下步骤：

(1)辐射定标1011。通过从卫星图像的头文件中读取定标系数，使用以下计算公式将各时相影像的数字量化值DN(Digital Number)转换成辐亮度，实现图像的辐射定标：

L_t＝Gain×DN+Bias

式中，L_t为图像的辐亮度，Gain为增益，Bias为偏置，单位均为W﹒m^-2﹒sr^-1﹒μm^-1，DN为卫星载荷观测值，无量纲。

(2)几何校正1012。利用卫星图像自带的有理多项式系数(Rational PolynomialCoefficients，RPC)，通过有理函数模型(Rational Function Model，RFM)对图像进行几何校正。

(3)图像配准1013。选定一个时相的遥感图像为基准图像，其它(n-1)个时相的遥感图像作为待配准图像，通过选取基准图像与待配准图像上的同名特征点，然后利用这些同名特征点建立基准图像与待配准图像之间的几何畸变模型，本例中选择二次多项式模型，再利用该几何畸变模型进行几何校正，从而实现不同时相遥感图像的图像配准。

夜间灯光图像的图像重采样的方法，其步骤是：

(1)图像重采样1014。通过图像灰度值重采样对夜间灯光图像进行空间分辨率变换，使得变换后的图像空间分辨率与高时频多光谱遥感图像的空间分辨率一致。采用数学上的双向线性内插法作为重采样方法计算分辨率变换后像元位置的灰度值。

双线性插值法是将与待求像元的四个邻近像元的灰度值在行和列两个方向上作线性内插。通常根据待求像元与邻近的四个像元的距离权重计算出灰度值，并赋给待求像元。双线性插值法公式为：

式中，g_x'y'为输出像元灰度值，g_i为邻近点i的灰度值，p_i为邻近点对投影点的权重(p_i＝1/d_i，d_i表示邻近点到投影点的距离)。g₁、g₂、g₃、g₄分别为四个邻近点的灰度值，p₁、p₂、p₃、p₄分别为四个邻近点对投影点的权重。

C、单元102居民地提取。基于上述预处理101后的高时频遥感图像和夜间灯光图像，通过逐步剔除背景地类实现居民地提取。步骤如下：

(1)NDVI最大值合成阈值分割1021。基于上述图像预处理101后的高时频多光谱遥感图像，根据NDVI计算公式计算每个时相的NDVI，并对各时相的NDVI进行最大值合成，得到NDVI最大值合成图。基于NDVI最大值合成图，通过阈值分割剔除植被覆盖地表。

根据NDVI计算公式计算每个时相的NDVI，包括以下步骤：

NDVI为归一化植被指数(Normalized Difference Vegetation Index)，其计算公式为：

其中，NDVI为归一化植被指数(无量纲)，NIR和R分别为图像的近红外波段和红波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1。

对各时相的NDVI进行最大值合成，得到NDVI最大值合成图，包括以下步骤：

采用最大值合成法(Maximum Value Composites,MVC)进行NDVI最大值合成，即选取观测像元的在时间序列图像中的最大NDVI值作为该像元合成期的NDVI。MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数。

基于NDVI最大值合成图，通过阈值分割剔除植被覆盖地表，包括以下步骤：

1)对于NDVI最大值合成图，设定阈值T1，T1为正数,本例中T1＝0.25；

2)对于NDVI最大值合成图中像元值大于T1的像元，将其归入植被覆盖地表类别(非居民地)，赋值为0，其他像元归入备选居民地类别，赋值为1。由于植被覆盖地表的NDVI通常大于居民地，通过设置合适的阈值，可以将植被覆盖地表从图像中剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域。

(2)NDWI最大值合成阈值分割1022。基于上述图像预处理101后的高时频多光谱遥感图像，根据NDWI计算公式计算每个时相的NDWI，并对各时相的NDWI进行最大值合成，得到NDWI最大值合成图。基于NDWI最大值合成图，通过阈值分割剔除水体。

根据NDWI计算公式计算每个时相的NDWI，包括以下步骤：

NDWI为归一化水体指数(Normalized Difference Water Index)，其计算公式为：

其中，NDWI为归一化水体指数(无量纲)，NIR和G分别为图像的近红外波段和绿波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1。

对各时相的NDWI进行最大值合成，得到NDWI最大值合成图，包括以下步骤：

采用最大值合成法(Maximum Value Composites,MVC)进行NDWI最大值合成，即选取观测像元的在时间序列图像中的最大NDWI值作为该像元合成期的NDWI。MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数。

基于NDWI最大值合成图，通过阈值分割剔除水体，包括以下步骤：

1)对于NDWI最大值合成图，设定阈值T2，T2为正数，本例中T2＝0.30；

2)针对NDVI最大值合成阈值分割1021结果图中像元值为1的区域(即备选居民地区域)，将NDWI最大值合成图中像元值小于T2的像元归入备选居民地区域，赋值为1，其他像元归入水体(非居民地)，赋值为0。由于大部分水体的NDWI通常大于居民地，通过设置合适的阈值，可以将水体从图像中剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域。

(3)近红外波段辐亮度变化斜率阈值分割1023。基于上述图像预处理101后的高时频多光谱遥感图像，采用趋势分析方法分析近红外波段辐亮度随太阳高度角变化的趋势特征，得到近红外波段辐亮度随太阳高度角的变化斜率图。基于近红外波段辐亮度随太阳高度角的变化斜率图，通过阈值分割剔除易与居民地混淆的水体。

采用趋势分析方法获取近红外波段辐亮度随太阳高度角的变化斜率图，包括以下步骤：

基于图像预处理101后的高时频影像，采用Theil-Sen Median趋势分析法分析近红外波段辐亮度随太阳高度角的变化趋势特征，具体是按照以下公式计算各像元近红外波段辐射亮度值随太阳高度角的变化斜率，生成近红外波段辐亮度变化斜率图。

其中，slope为图像像元近红外波段辐亮度值随太阳高度角的变化斜率(无量纲)，NIR_i和NIR_j分别为i时相和j时相的近红外波段辐亮度值(无量纲)，Sun_i和Sun_j分别为i时相和j时相的太阳高度角值(无量纲)。median是一种计算机函数，返回一组数值的中值，即将这组数值按照从小到大的顺序排列，返回居于中间的数值。如果这组数值集合中包含偶数个数值，则返回位于中间的两个数的平均值。

对近红外波段辐亮度变化斜率图设定阈值进行水体剔除，包括以下步骤：

1)对于生成的近红外波段辐亮度变化斜率图，设定斜率阈值T3，T3为正数,本例中T3＝0.35；

2)针对NDWI最大值合成阈值分割1022结果图中像元值为1的区域(即备选居民地区域)，将斜率slope大于T3的像元归入备选居民地区域，赋值为1，其他像元归入水体(非居民地)，赋值为0。由于水体像元的斜率slope通常小于居民地，通过设置合适的阈值，可以将低亮居民地与水体区分。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域。

(4)绿、红、近红外波段辐亮度之和最大值合成1024。基于上述图像预处理(101)后的高时频多光谱遥感图像，计算每个时相图像的绿、红、近红外波段辐亮度之和，并对各时相的辐亮度和进行最大值合成，得到辐亮度和的最大值合成图。

计算每个时相图像的绿、红、近红外波段辐亮度之和包括以下步骤：

绿、红、近红外波段辐亮度之和计算公式为：

MSS＝G+R+NIR

其中，MSS为绿、红、近红外波段辐亮度之和，G、R和NIR分别为图像的绿、红、近红外波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1。

对各时相的辐亮度和进行最大值合成，得到辐亮度和的最大值合成图，包括以下步骤：

采用最大值合成法(Maximum Value Composites,MVC)进行辐亮度和最大值合成，即选取观测像元的在时间序列图像中的最大辐亮度和值作为该像元合成期的辐亮度和。MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数。

(5)宽约束下的居民地初提取1025。基于绿、红、近红外波段辐亮度之和最大值合成1024结果图，通过设置阈值剔除辐亮度较低的阴影、水体，得到宽约束下的居民地初提取结果。包括以下步骤：

1)对于绿、红、近红外波段辐亮度之和最大值合成1024结果图，设定阈值T4，T4为正数，T4的设置在保证居民地完整的前提下，尽量剔除阴影、水体等非居民地,本例中T4＝110W﹒m^-2﹒sr^-1﹒μm^-1；

2)针对近红外波段辐亮度变化斜率阈值分割1023结果图中像元值为1的区域(即备选居民地区域)，将辐亮度和最大值合成图中像元值大于T4的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0。由于水体、阴影的辐亮度和通常小于居民地，通过设置合适的阈值，可以在保证居民地完整的前提下，将这些非居民地区域从图像中尽量剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为宽约束下的居民地区域。宽约束下的居民地初提取结果中包含了完整的居民地，但仍包含部分错提的辐亮度较高的云、雪、裸地等像元。

(6)严约束下的居民地初提取1026。基于绿、红、近红外波段辐亮度之和最大值合成1024结果图，通过设置阈值剔除辐亮度较高的云、雪、裸地等，得到严约束下的居民地初提取结果。包括以下步骤：

1)对于绿、红、近红外波段辐亮度之和最大值合成1024结果图，设定阈值T5，T5为正数，T5的设置在保证较彻底剔除云、雪、裸地等非居民地的前提下，尽量保留较完整的居民地,本例中T5＝230W﹒m^-2﹒sr^-1﹒μm^-1；

2)针对宽约束下的居民地初提取1025结果图中像元值为1的区域(即备选居民地区域)，将辐亮度和最大值合成图中像元值小于T5的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0。由于云、雪、裸地的辐亮度和通常大于居民地，通过设置合适的阈值，可以将这些非居民地区域从图像中较彻底的剔除。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为严约束下的居民地区域。严约束下的居民地初提取结果中较彻底的剔除了云、雪、裸地等非居民地噪声，但居民地区域提取不够完整，一些辐亮度较高的城镇区域容易漏提取。

(7)夜间灯光严约束下的城镇提取1027。基于图像预处理101后的夜间灯光数据，通过设置较严格的阈值得到基于灯光指数的城镇区域提取结果，并将其与宽约束下的居民地初提取(1025)结果结合，得到夜间灯光严约束下的城镇提取结果。包括以下步骤：

1)对于图像预处理101后的夜间灯光数据，设定阈值T6，T6为正数，T6的设置在保证较彻底剔除非城镇区域的前提下，尽量保证城镇的完整,本例中T6＝240；

2)针对宽约束下的居民地初提取1025结果图中像元值为1的区域(宽约束下的居民地区域)，将夜间灯光指数大于T6的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0。

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为夜间灯光严约束下的城镇区域。

(8)结果合并1028将严约束下的居民地初提取1026。结果与夜间灯光严约束下的城镇提取1027结果进行合并，得到完整且噪声尽量小的最终的居民地提取结果。包括以下步骤：

针对严约束下的居民地初提取1026结果图中像元值为1的区域(严约束下的居民地区域)和夜间灯光严约束下的城镇提取1027结果图中像元值为1的区域(夜间灯光严约束下的城镇区域)归入居民地区域，赋值为1，其他像元归入非居民地，赋值为0。

居民地提取结果图(图4)中像元值为0的区域为居民地区域，像元值为1的区域为居民地区域。

D、单元103居民地分类。基于上述居民地提取102结果，通过居民地连续斑块的面积和夜间灯光的阈值分割实现居民地分类。

连续斑块面积、夜间灯光阈值分割1031。基于上述居民地提取102结果图，通过设置居民地图斑连续斑块的面积阈值和夜间灯光阈值，将居民地中的城镇与农村居民点区分开来。包括以下步骤：

1)对于居民地提取102结果图，设定连续斑块面积阈值T7，T7为正数，T7设置为保证小于最小的城镇图斑面积的前提下，尽量大于农村居民点图斑面积，本例中T7＝4km²；

2)对于图像预处理101后的夜间灯光数据，设定阈值T8，T8为正数，T8的设置在保证城镇完成的前提下，尽量剔除非城镇区域，本例中T8＝65；

3)对居民地提取102结果图中连续斑块的面积大于阈值T7的居民地图斑，将平均夜间灯光指数大于T8的图斑归入城镇，图斑内的像元赋值为2；对居民地提取102结果图中连续斑块的面积大于阈值T7的居民地图斑，将平均夜间灯光指数小于T8的图斑归入农村居民点，图斑内的像元赋值为1；对居民地提取102结果图中连续斑块的面积小于等于阈值T7的居民地图斑，归入农村居民点，图斑内的像元赋值为1；其他像元归入非居民地，赋值为0。

居民地分类结果图(图5)中像元值为0的区域为非居民地区域，像元值为1的区域为农村居民点，像元值为2的区域为城镇。

E、对采用本发明方法进行高时频遥感图像特征辅助的居民地提取与采用SVM监督分类方法(Cortes和Vapnik于1995年提出，是机器学习领域一种有监督的学习模型，通常用来进行模式识别、分类以及回归分析)居民地提取的结果进行比较(图6)。SVM监督分类方法的结果对训练区的选择具有较强的依赖性，且分类速度较慢，耗时长，易将水体、裸地与居民地混淆。本发明方法速度快，能够较好的将水体、裸地等与居民地区分开来，提高检测精度。

F、对采用本发明方法进行高时频遥感图像特征辅助的居民地提取与采用SVM监督分类方法居民地提取的精度定量比较(图7)。图中的居民地提取精度的评价方法如下：

首先对需要进行精度评价的区域进行基于人工目视解译的居民地提取，然后将目视解译的居民地提取结果作为参考依据对自动方法的居民地提取结果进行基于混淆矩阵的精度评价，混淆矩阵格式如下：

居民地提取混淆矩阵

通过计算居民地提取的漏检误差A_Omission和错检误差A_Commision评价居民地提取的精度，它们的计算公式分别是：

其中，N₁₁表示在自动方法居民地提取结果中为居民地而参考数据中也为居民地的像素数，N₁₂表示在自动方法居民地提取结果中为居民地而参考数据中为非居民地的像素数，N₂₁表示在自动方法居民地提取结果中为非居民地而参考数据中为居民地的像素数，N₂₂表示在自动方法居民地提取结果中为非居民地而参考数据中也为非居民地的像素数。

本实施例中相比于传统的监督分类法等基于光谱的居民地提取方法，本发明方法具有更高的居民地提取精度，其漏检误差和错检误差均降低，能够更好的消除水体、裸地的影响，提高居民地检测精度。

通过上述具体实施例，在居民地提取中引入地物反射光谱随时间(太阳高度角)的变化特征信息，实现了实验区居民地的高精度提取，对于基于光谱的居民地提取方法中容易与居民地混淆的水体和裸地能够较好的区分，同时结合夜间灯光数据能够很好的实现居民地二级分类。同是也一定程度上提高了居民地提取的速度和自动化程度。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。

Claims (1)

1.一种高时频遥感图像特征辅助的居民地提取与分类方法，其步骤是：

(A)收集覆盖研究区域的不同太阳高度角的高时频多光谱遥感图像和夜间灯光图像，研究区域高时频多光谱遥感图像收集的包括：

1)每景遥感图像均覆盖研究区域，每景图像的云覆盖区域不大于10％；

2)每景遥感图像均包含绿、红、近红外波段，同一传感器拍摄的图像；

3)收集不同时相遥感图像n景(n>＝2)，不同时相图像上的地物类型未发生变化，数据的时间跨度不超过10天；

4)遥感图像拍摄时的太阳高度角在30－75度，太阳高度角和最小太阳高度角之间的差异不小于10度；

研究区域夜间灯光图像收集的包括：采用NPP/VIIRS DNB年合成夜间灯光图像数据，数据采集时间与高时频多光谱遥感图像采集时间接近；

(B)对步骤(A)收集的高时频多光谱遥感图像进行辐射定标、几何校正和图像配准，对夜间灯光图像进行图像重采样处理；

1)通过从卫星图像的头文件中读取定标系数，使用以下计算公式将各时相多光谱遥感图像的数字量化值DN转换成辐亮度，实现图像的辐射定标：

L_t＝Gain×DN+Bias

式中，L_t为图像的辐亮度，Gain为增益，Bias为偏置，单位均为W﹒m^-2﹒sr^-1﹒μm^-1，DN为卫星载荷观测值，无量纲；

2)利用卫星图像自带的有理多项式系数，通过有理函数模型对各时相多光谱遥感图像进行几何校正；

3)选定多时相遥感图像中一个时相的遥感图像为基准图像，其它(n-1)个时相的遥感图像作为待配准图像，通过选取基准图像与待配准图像上的同名点，然后利用这些同名点建立基准图像与待配准图像之间的几何畸变模型，再利用该几何畸变模型进行几何校正；

4)通过图像灰度值重采样对夜间灯光图像进行空间分辨率变换，使得变换后的图像空间分辨率与高时频多光谱遥感图像的空间分辨率一致，采用数学上的双向线性内插法为重采样方法计算分辨率变换后像元位置的灰度值；

双线性插值法是将与待求像元的四个邻近像元的灰度值在行和列两个方向上作线性内插，通常根据待求像元与邻近的四个像元的距离权重计算出灰度值，并赋给待求像元，双线性插值法公式为：

式中，g_x'_y'为输出像元灰度值，g_i为邻近点i的灰度值，p_i为邻近点对投影点的权重，p_i＝1/d_i，d_i表示邻近点到投影点的距离；g₁、g₂、g₃、g₄分别为四个邻近点的灰度值，p₁、p₂、p₃、p₄分别为四个邻近点对投影点的权重；

(C)基于(步骤B)预处理后的高时频遥感图像和夜间灯光图像，通过逐步剔除背景地类实现居民地提取，其步骤包括：

1)基于(步骤B)图像预处理后的高时频多光谱遥感图像，根据NDVI计算公式计算每个时相的NDVI，并对各时相的NDVI进行最大值合成，得到NDVI最大值合成图，基于NDVI最大值合成图，通过阈值分割剔除植被覆盖地表；

NDVI为归一化植被指数，其计算公式为：

其中，NDVI为归一化植被指数，NIR和R分别为图像的近红外波段和红波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1；

采用最大值合成法对各时相的NDVI进行最大值合成，即选取观测像元的在时间序列图像中的最大NDVI值作为该像元合成期的NDVI，MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数；

基于NDVI最大值合成图，通过阈值分割剔除植被覆盖地表，包括以下步骤：

①对于NDVI最大值合成图，设定阈值T1，T1为正数；

②对于NDVI最大值合成图中像元值大于T1的像元，将其归入植被覆盖地表类别，赋值为0，其他像元归入备选居民地类别，赋值为1，植被覆盖地表的NDVI通常大于居民地，通过设置阈值，将植被覆盖地表从图像中剔除；

2)基于(步骤B)图像预处理后的高时频多光谱遥感图像，根据NDWI计算公式计算每个时相的NDWI，并对各时相的NDWI进行最大值合成，得到NDWI最大值合成图，基于NDWI最大值合成图，通过阈值分割剔除水体；

NDWI为归一化水体指数，其计算公式为：

其中，NDWI为归一化水体指数，NIR和G分别为图像的近红外波段和绿波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1；

采用最大值合成法进行NDWI最大值合成，选取观测像元的在时间序列图像中的最大NDWI值作为该像元合成期的NDWI，MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数；

基于NDWI最大值合成图，通过阈值分割剔除水体，包括以下步骤：

①对于NDWI最大值合成图，设定阈值T2，T2为正数；

②针对NDVI最大值合成阈值分割结果图中像元值为1的区域，将NDWI最大值合成图中像元值小于T2的像元归入备选居民地区域，赋值为1，其他像元归入水体，赋值为0，水体的NDWI通常大于居民地，通过设置合适的阈值，将水体从图像中剔除；

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域；

3)基于(步骤B)图像预处理后的高时频多光谱遥感图像，采用趋势分析方法分析近红外波段辐亮度随太阳高度角变化的趋势特征，得到近红外波段辐亮度随太阳高度角的变化斜率图，基于近红外波段辐亮度随太阳高度角的变化斜率图，通过阈值分割剔除易与居民地混淆的水体；

采用趋势分析方法，计算n(n-1)/2个数据组合的斜率的中位数，获得近红外波段辐亮度随太阳高度角的变化斜率图，趋势分析方法的公式为：

其中，slope为图像像元近红外波段辐亮度值随太阳高度角的变化斜率，NIR_i和NIR_j分别为i时相和j时相的近红外波段辐亮度值，Sun_i和Sun_j分别为i时相和j时相的太阳高度角值；median是一种计算机函数，返回一组数值的中值，将这组数值按照从小到大的顺序排列，返回居于中间的数值，这组数值集合中包含偶数个数值，返回位于中间的两个数的平均值；

对近红外波段辐亮度变化斜率图设定阈值进行水体剔除，包括以下步骤：

①对于生成的近红外波段辐亮度变化斜率图，设定斜率阈值T3，T3为正数；

②针对NDWI最大值合成阈值分割结果图中像元值为1的区域，将斜率slope大于T3的像元归入备选居民地区域，赋值为1，其他像元归入水体，赋值为0，水体像元的斜率slope通常小于居民地，通过设置合适的阈值，将低亮居民地与水体区分；

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为备选居民地区域；

4)基于(步骤B)图像预处理后的高时频多光谱遥感图像，计算每个时相图像的绿、红、近红外波段辐亮度之和，并对各时相的辐亮度和进行最大值合成，得到辐亮度和的最大值合成图；

绿、红、近红外波段辐亮度之和计算公式为：

MSS＝G+R+NIR

其中，MSS为绿、红、近红外波段辐亮度之和，G、R和NIR分别为图像的绿、红、近红外波段辐亮度值，单位均为W﹒m^-2﹒sr^-1﹒μm^-1；

采用最大值合成法进行辐亮度和最大值合成，选取观测像元的在时间序列图像中的最大辐亮度和值作为该像元合成期的辐亮度和，得到辐亮度和的最大值合成图，MVC法的公式为：

X＝Max(X₁，...，X_i，...，X_n)，1≤i≤n

其中，X为某像元的最大值合成结果，X_i为该像元第i时相的值，n为图像时相数；Max为求取一组数中的最大值的函数；

5)基于绿、红、近红外波段辐亮度之和最大值合成结果图，通过设置阈值剔除辐亮度的阴影、水体，得到宽约束下的居民地初提取结果，包括以下步骤：

①对于绿、红、近红外波段辐亮度之和最大值合成结果图，设定阈值T4，T4为正数，T4的设置在居民地完整的前提下，剔除阴影、水体非居民地；

②针对近红外波段辐亮度变化斜率阈值分割结果图中像元值为1的区域，将辐亮度和最大值合成图中像元值大于T4的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0，水体、阴影的辐亮度和通常小于居民地，通过设置阈值，在居民地完整的前提下，将这些非居民地区域从图像中尽量剔除；

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为宽约束下的居民地区域，宽约束下的居民地初提取结果中包含了完整的居民地，包含错提的辐亮度的云、雪、裸地像元；

6)基于绿、红、近红外波段辐亮度之和最大值合成结果图，通过设置阈值剔除辐亮度的云、雪、裸地，得到居民地初提取结果，包括以下步骤：

①对于绿、红、近红外波段辐亮度之和最大值合成结果图，设定阈值T5，T5为正数，T5的设置在剔除云、雪、裸地非居民地的前提下，保留较完整的居民地；

②针对宽约束下的居民地初提取结果图中像元值为1的区域，将辐亮度和最大值合成图中像元值小于T5的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0，云、雪、裸地的辐亮度和通常大于居民地，通过设置阈值，将非居民地区域从图像中彻底的剔除；

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为严约束下的居民地区域，严约束下的居民地初提取结果中彻底的剔除了云、雪、裸地非居民地噪声，居民地区域提取不够完整，辐亮度高的城镇区域容易漏提取；

7)基于(步骤B)图像预处理后的夜间灯光数据，通过设置的阈值得到基于灯光指数的城镇区域提取结果，并将其与宽约束下的居民地初提取结果进行结合，得到夜间灯光严约束下的城镇提取结果，包括以下步骤：

①对于图像预处理后的夜间灯光数据，设定阈值T6，T6为正数，T6的设置在彻底剔除非城镇区域的前提下，城镇的完整；

②针对宽约束下的居民地初提取结果图中像元值为1的区域，将夜间灯光指数大于T6的像元归入备选居民地区域，赋值为1，其他像元归入非居民地，赋值为0；

结果图中像元值为0的区域为非居民地区域，像元值为1的区域为夜间灯光严约束下的城镇区域；

8)将严约束下的居民地初提取结果与夜间灯光严约束下的城镇提取结果进行合并，得到完整且噪声的最终的居民地提取结果，包括以下步骤：

针对严约束下的居民地初提取结果图中像元值为1的区域和夜间灯光严约束下的城镇提取结果图中像元值为1的区域归入居民地区域，赋值为1，其他像元归入非居民地，赋值为0，居民地提取结果图中像元值为0的区域为居民地区域，像元值为1的区域为居民地区域；

(D)基于(步骤C)居民地提取结果，通过居民地连续斑块的面积和夜间灯光的阈值分割，将居民地中的城镇与农村居民点区分开来，实现居民地分类，得到居民地分类图，包括以下步骤：

①对于居民地提取结果图，设定连续斑块面积阈值T7，T7为正数，T7设置为小于最小的城镇图斑面积的前提下，大于农村居民点图斑面积；

②对于图像预处理后的夜间灯光数据，设定阈值T8，T8为正数，T8的设置在城镇完成的前提下，剔除非城镇区域；

③对居民地提取结果图中连续斑块的面积大于阈值T7的居民地图斑，将平均夜间灯光指数大于T8的图斑归入城镇，图斑内的像元赋值为2；对居民地提取结果图中连续斑块的面积大于阈值T7的居民地图斑，将平均夜间灯光指数小于T8的图斑归入农村居民点，图斑内的像元赋值为1；对居民地提取结果图中连续斑块的面积小于等于阈值T7的居民地图斑，归入农村居民点，图斑内的像元赋值为1；其他像元归入非居民地，赋值为0，得到居民地分类结果图；

居民地分类结果图中像元值为0的区域为非居民地区域，像元值为1的区域为农村居民点，像元值为2的区域为城镇。