本发明属于夜光遥感在城市发展研究领域的应用，提出一种全新的基于夜光遥感影像的城市范围提取方法。
背景技术：
中国城市化进程在20世纪80年代开始加快。城市化是衡量一个国家发展状况的主要标志之一。但是这种快速的城市化进程同时也为人们的生活环境和生态系统带来很多的影响。城市化与工业化的关系、城市化是否滞后、未来20年城市化的速度应该为多少、应该发展什么规模的城市这四个问题是我国目前城市化过程中遇到的主要问题。为解决这些问题，我们必须能够及时地获得城市确切的范围和城市变化的时空趋势。
航空航天遥感技术的发展为城市范围提取提供了一种高效的手段，相比于传统的地面实际测量的方法，遥感技术有着工作量小、成本低、周期短、效率高等特点，并且能够满足当前研究城市化的需求。传统的遥感技术提取城市范围的方式一般是利用一年内较高分辨率的多光谱的遥感影像进行提取，其提取过程主要包括影像的几何校正、大气校正、镶嵌、裁剪、分类等步骤。由于多光谱的遥感影像获取受到天气的影响较大，因此所得到的影像在进行几何校正、镶嵌等工作时由于成像时间的不同导致其操作比较困难。并且一般使用的多光谱影像的空间分辨率较高，这样会使得数据量变大，对于数据处理、存储、管理、发布等都造成一定的困难。
夜光遥感始于20世纪70年代美国国防气象卫星计划的线性扫描业务系统(defensemeteorologicalsatelliteprogram/operationallinescansystem,dmsp/ols)，该系统起初是为探测夜间的云层反射的微弱夜光来获得夜间云层分布的信息[1]，但是相关科学家发现该传感器能在无云的条件下获得地球表面上发射的可见光，比如城镇照明设施的发光、渔船发光、油气井的燃烧发光等。相比于传统的日间遥感，夜光遥感影像成像于夜间，而夜间唯一能直接反映人类活动的是照明设施发光，鉴于夜光遥感的这种独特的特性，因此它更多的是反映地球表面城市的时空分布，这为提取城市范围提供了一种新的数据源。目前所能获得的夜光影像有dmsp/ols的年产品和美国极轨卫星(suominationalpoalorbitingpartnership，suominpp)搭载的可见光近红外成像辐射传感器(visibleinfraredimagingradiometersuite,viirs)所得到的月产品和年产品。目前发布的这两种夜光遥感影像全都为全球合成产品，不需要进行大气校正、镶嵌等处理，并且由于数据量小，因此利用夜光遥感影像进行城市范围提取是研究城市发展的一种新的高效的手段。
目前已有很多学者研究了各种利用夜光影像提取城市范围的方法，其主要以阈值分割方法为主。其阈值分割方法又包括单一阈值法、分区域阈值法和基于对象阈值法。
单一阈值方法最初由marcl.imhoff等学者[2]提出，该方法则是在研究区域内选择一个灯光亮度值作为城市分割阈值，即亮度值大于该阈值的像元为城市像元，反之则为非城市像元。该方法简单易于实现，但是随着研究区域范围的变大，研究区域内的发展越不平衡，导致夜间灯光分布的差异性变大，则单一阈值法无法再满足大范围研究要求[3]。
针对单一阈值方法的缺陷，zhang等学者[4]最初提出分区域阈值的方法。该方法则是根据研究区域内的经济、人口、地理位置等，将研究区域分为多个较小的研究区域，然后在每个分区中独立选择阈值，最后根据区域阈值得到整个研究区域的城市范围。虽然该方法较单一阈值方法更加合理，但是其所得城市范围会随着分区标准的不同而发生变化，并且即使在城市内部灯光的分布任然存在着很大的差异。
为了弥补单一阈值法和分区域阈值法的不足，已有文献[5-6]提出基于对象的阈值方法。即是先对夜光影像根据灯光亮度值进行分割，然后再计算每个对象内的最佳阈值。这种基于对象的方法可以解决由于灯光的空间和亮度分布差异造成的城市范围精度问题。文献[5]同时也表明灯光阈值的选择不仅与灯光亮度值之间存在关系，还与分割对象的尺寸、对象内均值等因素相关，并且符合一定的函数分布关系。但是函数模型随着数据年份的不同而不同，需要根据不同的年份来设计函数的参数，并且通过散点图和统计数据，很难较准确的模拟阈值和灯光数据之间的关系。
人工神经网络(artificialneuralnetwork,ann)是生物神经网络研究的基础上所发展出的一种计算方法，该方法基于学习训练样例以完成分类等特定任务，并能够在训练的过程中逐渐改进性能[7]。人工神经网络是一种试图模拟生物神经网络的结构和功能的数学模型。其基本组成是人工神经元，即简单的数学函数。人工神经元包含三组运算规则：乘法、求和和激活。在人工神经元的入口处，输入被加权，即每个输入值乘以单独的权重；在神经元的中间部分，对所有输入的加权值和其偏差进行求和。在出口处，先前加权的输入和偏差的总和将通过激活函数进行最后的运算。
为了充分发挥人工神经元的互连而得到数学复杂性的优势，又同时避免提高系统的复杂性导致其难以管理，人工神经元的连接通常依照一定的规律和范式。一些“标准化”的人工神经网络拓扑图已被提出，这些预定义的形状可以帮助使用者更简单快速有效地解决问题[8]。不同类型的人工神经网络拓扑适用于解决不同类型的问题。在确定给定问题的类型后，决定将采用的人工神经网络的拓扑结构，同时为了发挥神经网络的最优效用，拓扑及其参数均需要进行微调。完成人工神经网络的拓扑构建只是完成了一半的任务。正如生物神经网络需要学习如何对环境中给定输入得出适当的输出，人工神经网络也需要通过监督、非监督或强化学习来得出如何根据输入得出反馈，最后得到基于样本数据的最优神经网络模型。
近年来，已有多种人工神经网络模型被提出并得以深入研究。其中80％～90％的人工神经网络模型是采用前馈反向传播网络(back-propagationnetwork简称bp网络)或它的改进形式[9]。bp神经网络它有一个输入层、一个输出层和一个或多个隐含层,同层神经元间无关联,异层神经元间向前连接。根据对象的复杂程度,选择适当的网络结构,就可以实现从输入空间到输出空间的任意非线性函数的映射。bp神经网络主要用于：函数逼近、系统辨识和预测、分类、数据压缩等领域[10]。
综合前人的研究，目前利用夜光遥感影像提取城市范围还没有比较成熟的方法，其中绝大多数方法都需要人工设置阈值，这将会使研究结果存在较大的人为因素造成的误差。同时基于对象提取城市范围方法中表明对象阈值与对象参数之间存在一定的函数关系。本专利则通过使用bp神经网络模型，以分割对象内的最大值、最小值、均值、方差、像元数和最佳阈值为样本数据对模型进行训练，利用训练好的bp神经网络来模拟输入参数和输出最佳阈值之间的某种函数关系。该方法能够自动设置阈值和提取城市范围，使得利用夜光影像提取城市范围摆脱人为因素的干扰，得到更加准确和可靠的结果。
相关文献如下：
[1]croftta.nighttimeimagesoftheearthfromspace[j].scientificamerican,1978,239(1):86-98.
[2]imhoffml,lawrencewt,stutzerdc,etal.atechniqueforusingcompositedmsp/ols“citylights”satellitedatatomapurbanarea[j].remotesensingofenvironment,1997,61(3):361-370.
[3]zhangq,setokc.cannight-timelightdataidentifytypologiesofurbanization？aglobalassessmentofsuccessesandfailures[j].remotesensing,2013,5(5):3476-3494.
[4]liuz,hec,zhangq,etal.extractingthedynamicsofurbanexpansioninchinausingdmsp-olsnighttimelightdatafrom1992to2008[j].landscape&urbanplanning,2012,106(106):62-72.
[5]zhouy,smithsj,elvidgecd,etal.acluster-basedmethodtomapurbanareafromdmsp/olsnightlights[j].remotesensingofenvironment,2014,147(18):173–185.
[6]xiey,wengq.updatingurbanextentswithnighttimelightimagerybyusinganobject-basedthresholdingmethod[j].remotesensingofenvironment,2016,187:1-13.
[7]whitleyd,starkweathert,bogartc.geneticalgorithmsandneuralnetworks:optimizingconnectionsandconnectivity[j].parallelcomputing,1990,14(3):347-361.
[8]李双成,郑度.人工神经网络模型在地学研究中的应用进展[j].地球科学进展,2003,18(1):68-76.
[9]杨兆升,朱中.基于bp神经网络的路径行程时间实时预测模型[j].系统工程理论与实践,1999,19(8):59-64.
[10]李萍,曾令可,税安泽,金雪莉,刘艳春,王慧.基于matlab的bp神经网络预测系统的设计[j].计算机应用与软件,2008(04):149-150+184.
技术实现要素：
综合前人研究并且针对目前利用夜光遥感影像提取城市范围的不足，本发明提出一种新的基于夜光遥感影像的城市范围提取方法。
为了解决上述技术问题，本发明采用如下的技术方案：
一种基于夜光遥感影像的城市范围提取方法，包括以下步骤：
步骤1，获取研究区域的夜间灯光dmsp/ols数据和部分研究区域土地利用类型数据，以及modis水体分布图；
步骤2，利用modis水体分布图建立二值图像对dmsp/ols夜光影像进行掩膜处理；
步骤3，对步骤2所得掩膜处理后的dmsp/ols夜光影像进行分割；
步骤4，计算步骤3所得每个分割对象内的像元数、夜光均值、最大值、最小值和方差；
步骤5，选择已有相应土地利用类型数据的分割对象，根据土地利用类型数据中的城市像元个数计算各分割对象的最佳灯光分割阈值；
步骤6，以步骤5所得分割对象内像元数、均值、最大值、最小值和方差为输入数据，相应最佳灯光分割阈值为输出值，对bp神经网络模型进行训练；
步骤7，利用步骤6训练好的bp神经网络模型计算所有分割对象的最佳阈值，利用阈值得到研究区域的最终城市范围。
而且，步骤5中，对每个分割对象分别执行以下子步骤，
步骤5.1，以夜光影像中对象内的像元最小值为阈值初值threshold_urban，大于阈值的灯光像元则为城市像元，反之则为非城市像元；
步骤5.2，分别计算对象内利用阈值得到的城市像元个数urban_dmsp和土地利用类型数据中的城市像元个数urban_nlcd；
步骤5.3，比较利用阈值提取的城市像元数和土地利用类型数据中的已知的城市像元数，计算两者的绝对误差值t＝|urban_dmsp-urban_nlcd|；
步骤5.4，从小到大依次选择对象内的像元值为新的阈值初值threshold_urban，重复步骤5.1-5.3，直到对象内所有像元的灯光亮度值全部经过计算；
步骤5.5，选择t最小时所对应的灯光亮亮度值为最佳灯光阈值threshold_urban。
而且，步骤5中，bp神经网络、的层数设置为输入层、2个隐含层和1个输出层，其中输入层有5个神经元，2个隐含层都含有10个神经元，最后输出层为1个神经元。
而且，步骤3中，采用标记控制分水岭分割算法对掩膜处理后的dmsp/ols影像进行分割。
与现有技术相比，本发明具有以下优点和有益效果：
1、本发明利用基于对象和神经网络的方法，能够使得计算结果更加合理和准确。
2、本发明使用神经网络代替具体的函数模型，避免了因函数模型的不准确而造成的计算结果有偏差，并且本发明方法数据处理和计算更加简单。
附图说明
图1为本发明得提取流程示意图。
图2为本发明实施例中bp神经网络模型的构造示意图以及输入输出变量。
具体实施方式
以下结合附图和实施例对本发明技术方案进行说明。
参见图1，本发明实施例包括以下步骤：
步骤1，获取研究区域的夜间灯光dmsp/ols影像和部分研究区域土地利用类型数据，以及modis水体分布图；
优选采用同一年份的数据，本步骤的具体操作如下：
分别下载研究区域同一年份的dmsp/ols夜光影像和研究区域中部分区域的土地利用类型数据(只需要部分区域的数据用于训练网络)以及modis水体分布图。
步骤2，利用modis水体分布图建立二值图像对dmsp/ols夜光影像进行掩膜处理；
本步骤的具体操作如下：
本步骤主要考虑到水体的反射光会增加地表的灯光量，通过掩膜处理则可以消除该影响，其计算公式如下：
式(3)中dmsp为夜光影像，本发明采用dmsp/ols夜光影像，lc_water为水体二值图像，i,j为影像的行列号，这样dmsp(i,j)表示夜光影像中像元(i,j)的取值，lc_water(i,j)表示水体二值图像中像元(i,j)的取值。
步骤3，利用标记控制分水岭分割算法对步骤2所得掩膜处理后的dmsp/ols夜光影像进行分割；
本步骤具体操作如下：
采用标记控制分水岭分割算法对掩膜处理后的dmsp/ols影像进行分割，分割窗口为3×3。标记控制分水岭分割算法可参见parvatik,raobsp,dasmm.imagesegmentationusinggray-scalemorphologyandmarker-controlledwatershedtransformation[j].discretedynamicsinnature&society,2013,2008(1026-0226):307-318.
步骤4，计算步骤3所得每个分割对象内的像元数、夜光均值、最大值、最小值和方差；
本步骤具体操作如下：
分别统计分割后每个对象的像元数、灯光均值、最大值、最小值、方差，其计算公式如下：
obj_size＝n.……(4)
obj_max＝max(dmsp(i)),i＝1,2,…,n.……(6)
obj_min＝min(dmsp(i)),i＝1,2,…,n.……(7)
式中obj_size、obj_mean、obj_max、obj_min、obj_var分别为对象的像元数、对象内像元均值、最大值、最小值、方差，其中n为该对象内的像元数，dmsp(i)为对象内的第i个像元在步骤2所得掩膜处理后的dmsp/ols夜光影像中的相应取值。
步骤5，选择已有相应土地利类型数据的分割对象，根据这些对象在土地利用类型数据中的城市像元个数计算最佳灯光阈值；
上述步骤5进一步包括对各对象分别执行以下子步骤：
步骤5.1，以夜光影像中对象内的像元最小值为阈值初值threshold_urban，大于阈值的灯光像元则为城市像元，反之则为非城市像元；
步骤5.2，分别计算对象内利用阈值得到的城市像元个数urban_dmsp和土地利用类型数据中的城市像元个数urban_nlcd；
步骤5.3，比较利用阈值提取的城市像元数和土地利用类型数据中的已知的城市像元数，计算两者的绝对误差值：t＝|urban_dmsp-urban_nlcd|；
步骤5.4，从小到大依次选择对象内的像元值为新的阈值初值threshold_urban，重复步骤5.1、5.2、5.3，直到对象内所有像元的灯光亮度值全部经过计算；
步骤5.5，选择t最小时所对应的灯光亮度值为最佳灯光阈值threshold_urban。
实施例中，本步骤的具体操作如下：
随机选择部分的对象，计算每个对象的最佳阈值，当|lc_sum-obj_sum|取最小值时所对应的阈值为最佳阈值，其计算公式如下：
obj_threshold(j)∈[min(dmsp),max(dmsp)].……(9)
其中dmsp为当前处理的分割对象，min(dmsp)为对象内的像元最大值，max(dmsp)为对象内的像元最小值，obj_threshold(j)为当前对象选择的第j个阈值初值threshold_urban，dmsp_urban(i)为当前对象第i个像元的值dmsp(i)经过阈值分割后的结果，1代表城市像元，0代表非城市像元，lc_urban(i)为当前对象第i个像元的土地利用类型数据，1代表城市像元，0代表非城市像元，lc_sum为土地利用类型数据中的城市像元个数，obj_sum为在某一阈值下dmsp/ols夜光影像对象中的城市像元个数，n为对象内像元总数。从小到大遍历所有的obj_threshold值，当|lc_sum-obj_sum|为最小值时threshold_urban＝obj_threshold(j)。
步骤6，以步骤5所得分割对象内像元数、夜光均值、最大值、最小值和方差为输入数据，相应最佳灯光分割阈值为输出值，作为网络训练样本对构建好的bp神经网络模型进行训练；
上述步骤6进一步包括以下子步骤：
步骤6.1，首先设计bp神经网络的层数和每个层的神经元激活函数，本专利中bp神经网络的层数设置为输入层、两个隐含层和一个输出层。其中输入层有五个神经元，两个隐含层都含有十个神经元，最后输出层为一个神经元。其中隐含层和输出层的激活函数为“tansig”和“purelin”函数。
步骤6.2，选择部分对象，分别以这些对象的最大值、最小值、均值、方差和对象尺寸为输入，其中选择的这部分对象的最佳阈值作为输出值，对模型进行训练。
步骤6.3，利用选择的部分对象训练的bp神经网络模型，对所有的对象进行处理，计算得到每个对象的最佳阈值。
实施例中，本步骤的具体操作如下：
如图2中所示，首先构造一个4层的bp神经网络初始模型，其中第一层为输入层，其第二、三层为隐含层，网络结点都为10个，第四层为输出层。
随机选择部分的对象进行bp神经网络模型的训练，其中样本的数据组成为输入数据(obj_size、obj_mean、obj_max、obj_min、obj_var)和输出数据(obj_threshold)。
步骤7，利用步骤6训练好的bp神经网络模型计算所有分割对象的最佳阈值，利用阈值得到研究区域的最终城市范围。
本步骤的具体操作如下：
7.1，利用训练好的网络模型对全部的分割对象进行处理，得到每个对象的最佳阈值。
7.2，根据bp神经网络模型计算得到的最佳阈值计算城市范围，包括对每个对象分别采用相应最佳阈值进行提取，大于阈值的灯光像元则为城市像元，反之则为非城市像元，综合所有对象提取结果得到城市范围。
下面将结合具体实例应用进一步说明本发明的技术方案及有益效果。
选择2005年中国区域的dmsp/ols夜光影像和全国部分城市土地利用类型数据以及modis水体分布图。首先利用利用modis水体影像对dmsp/ols夜光影像进行掩膜处理，然后对掩膜处理后的dmsp/ols夜光影像利用图像分割算法进行分割。随机选择部分对bp神经网络模型进行训练，最后将训练好的bp神经网络模型运用到整个影像中，得到2005年城市范围。对实验结果和真实土地利用数据在省级尺度上进行回归分析，结果显示回归所得r2＝0.85,rmse＝421.8。
具体实施时，本发明提供的流程可采用软件技术实现自动运行。
以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。