基于夜间灯光数据与高分辨率影像的城市建成区提取方法

摘要

本发明涉及城市遥感信息技术应用领域，具体涉及一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，包括：采集同一区域的夜间灯光遥感图像和高分辨率全色遥感图像，并分别对这两种图像进行预处理；对高分辨率全色遥感图像样本内的建筑区域与其它区域进行分割标注；将同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率进行调整，使两种图像的分辨率相同；利用同一区域的两种图像中每行每列的像素及其权重构建交叉熵损失函数；利用交叉熵损失函数得到训练好的语义分割网络；将预处理后的两种图像输入到训练好的语义分割网络中，完成城市建成区的提取。上述方法用于城市建成区的提取，可提高建成区的提取准确性。

权利要求

1.一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，其特征在于，包括：S1：采集同一区域的NPP/VIIRS夜间灯光遥感图像和高分辨率全色遥感图像，并分别对这两种图像进行预处理；

S2：构建语义分割网络模型：

S201：对预处理后的高分辨率全色遥感图像样本内的建筑区域与其它区域进行分割标注；

S202：将同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率进行调整，使同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率相同；

S203：利用同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像中每行每列的像素及其权重构建交叉熵损失函数；

所述交叉熵损失函数的表达式为：

式中，表示交叉熵损失函数，M、N分别为图像的宽、高， 表示图像中第i行第j列像素的权重， 表示语义分割网络预测的第i行第j列像素的概率， 表示语义分割网络的第i行第j列像素的标签值；

所述图像中第i行第j列像素的权重的表达式为：=

式中， 表示图像中第i行第j列像素的权重， 分别表示三种不同情况得到的第i行第j列像素的初始权重；

所述 是按照以下方式获取：

将夜间灯光遥感图像分为三个区域，该三个区域包括：建成区、过渡区和非建成区；

对建成区和非建成区进行中值滤波去噪，然后在两个区域中各选取一个中心邻域，并对中心邻域内各个像素的灰度值按大小进行排序，选取排序序列中的中值作为该中心邻域中心像素的灰度值；

对建成区和非建成区进行形态重构，获取形态重构后的夜间灯光遥感图像；

获取形态重构后的夜间灯光遥感图像的梯度图像；

利用分水岭分割算法提取梯度图像中的建成区；

将梯度图像中的建成区和非建成区进行标记，得到城市建成区的二值图，将该城市建成区的二值图的像素值作为；

所述 的计算公式为：

式中，表示取最优阈值时所有城市斑块的周长之和， 表示每个城市斑块的平均灯光强度， 表示高斯核函数， 表示每个城市斑块中像素的坐标；

所述 的表达式为：

式中， 表示取最优阈值时第n个城市斑块的周长，Q表示取最优阈值时城市斑块的数量；

所述 的表达式为：

=

式中， 表示取最优阈值时第 个城市斑块的年平均夜间灯光总量；

所述 的表达式为：

式中，是带宽；

所述 是按照以下方式获取：

从高分辨率全色遥感图像样本中提取NDVI和NDBI，其中NDVI为归一化植被指数，NDBI为归一化建筑指数， 和 的表达式分别为：式中， 为近红外波段，为红外波段， 为中红外波段；

采用改进的城市夜间灯光指数法VBANUI提取城市建成区， 的计算公式如下：式中， 为夜间灯光遥感图像；

在完成 指标的计算后，令

式中， 为每个城市斑块的最大值；

S204：利用构建的交叉熵损失函数对语义分割网络模型进行监督训练，直至语义分割网络收敛为止，得到训练好的语义分割网络；

S3：将预处理后的同一区域的NPP/VIIRS夜间灯光遥感图像和高分辨率全色遥感图像输入到训练好的语义分割网络模型中，完成城市建成区的提取。

2.根据权利要求1所述的一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，其特征在于，所述最优阈值是按照以下方式获取：设置初始阈值，并将该初始阈值为基础，以间隔权重为间隔进行迭代增加；

利用迭代后的结果获得城市斑块的周长随阈值迭代的分布情况；

找出城市斑块的周长突变时所对应的辐射值，该辐射值减掉间隔权重即为最优阈值。

说明书

基于夜间灯光数据与高分辨率影像的城市建成区提取方法

技术领域

[0001] 本发明涉及城市遥感信息技术应用领域，具体涉及一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法。

背景技术

[0002] NPP/VIIRS夜间灯光数据由于其较高的分辨率以及覆盖面广、成本低廉、效率高等特点，被广泛用来衡量人类活动的强度与广度，主要应用于城市建成区提取、区域经济、地缘政治等研究领域。在目前利用夜间灯光数据提取城市建成区的研究中，常用的最佳阈值分割法主要有四种：经验阈值法、突变检测法、统计数据验证法和较高分辨率影像比较法。关于高分辨率影像提取城市建成区的研究大致可以分为两类：一类是基于区域增长方法的半自动提取方法，另一类方法则主要是基于分类的思想，根据光谱和纹理特征进行分类并作进一步处理以达到提取城市建成区的目的。

[0003] 然而，利用NPP/VIIRS夜间灯光数据提取城市建成区仍然存在挑战：主要是NPP/VIIRS夜间灯光数据受到溢出效应的影响。由于非相干光从光源向各个方向辐射，比如，灯光亮度在周围区域的弥散，造成灯光亮度的溢出效应，导致对城市土地用地面积的高估，限制了NPP/VIIRS夜间灯光数据在城市土地用地精确提取的广泛应用。

发明内容

[0004] 为了克服上述现有技术的缺点，本发明提供了一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，以实现对城市建成区的精确提取。

[0005] 为达到上述目的，本发明采用以下技术方案，一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，包括：

[0006] S1：采集同一区域的NPP/VIIRS夜间灯光遥感图像和高分辨率全色遥感图像，并分别对这两种图像进行预处理；

[0007] S2：构建语义分割网络模型：

[0008] S201：对预处理后的高分辨率全色遥感图像样本内的建筑区域与其它区域进行分割标注；

[0009] S202：将同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率进行调整，使同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率相同；

[0010] S203：利用同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像中每行每列的像素及其权重构建交叉熵损失函数；

[0011] S204：利用构建的交叉熵损失函数对语义分割网络模型进行监督训练，直至语义分割网络收敛为止，得到训练好的语义分割网络；

[0012] S3：将预处理后的同一区域的NPP/VIIRS夜间灯光遥感图像和高分辨率全色遥感图像输入到训练好的语义分割网络模型中，完成城市建成区的提取。

[0013] 所述一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，所述交叉熵损失函数的表达式为：

[0014]

[0015] 式中，表示交叉熵损失函数，M、N分别为图像的宽、高， 表示图像中第i行第j列像素的权重， 表示语义分割网络预测的第i行第j列像素的概率， 表示语义分割网络的第i行第j列像素的标签值。

[0016] 所述一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，所述图像中第i行第j列像素的权重的表达式为：

[0017] =

[0018] 式中， 表示图像中第i行第j列像素的权重， 分别表示三种不同情况得到的第i行第j列像素的初始权重。

[0019] 所述一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，所述 是按照以下方式获取：

[0020] 将夜间灯光遥感图像分为三个区域，该三个区域包括：建成区、过渡区和非建成区；

[0021] 对建成区和非建成区进行中值滤波去噪，然后在两个区域中各选取一个中心邻域，并对中心邻域内各个像素的灰度值按大小进行排序，选取排序序列中的中值作为该中心邻域中心像素的灰度值；

[0022] 对建成区和非建成区进行形态重构，获取形态重构后的夜间灯光遥感图像；

[0023] 获取形态重构后的夜间灯光遥感图像的梯度图像；

[0024] 利用分水岭分割算法提取梯度图像中的建成区；

[0025] 将梯度图像中的建成区和非建成区进行标记，得到城市建成区的二值图，将该城市建成区的二值图的像素值作为 。

[0026] 所述一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，所述 的计算公式为：

[0027]

[0028] 式中，表示取最优阈值时所有城市斑块的周长之和， 表示每个城市斑块的平均灯光强度， 表示高斯核函数， 表示每个城市斑块中像素的坐标；

[0029] 所述 的表达式为：

[0030]

[0031] 式中， 表示取最优阈值时第n个城市斑块的周长，Q表示取最优阈值时城市斑块的数量；

[0032] 所述 的表达式为：

[0033] =

[0034] 式中： 表示取最优阈值时第 个城市斑块的年平均夜间灯光总量；

[0035] 所述 的表达式为：

[0036]

[0037] 式中，是带宽。

[0038] 所述一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，所述最优阈值是按照以下方式获取：

[0039] 设置初始阈值，并将该初始阈值为基础，以间隔权重为间隔进行迭代增加；

[0040] 利用迭代后的结果获得城市斑块的周长随阈值迭代的分布情况；

[0041] 找出城市斑块的周长突变时所对应的辐射值，该辐射值减掉间隔权重即为最优阈值。

[0042] 所述一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，所述 是按照以下方式获取：

[0043] 从高分辨率全色遥感图像样本中提取NDVI和NDBI，其中NDVI为归一化植被指数，NDBI为归一化建筑指数， 和 的表达式分别为：

[0044]

[0045]

[0046] 式中， 为近红外波段，为红外波段， 为中红外波段；

[0047] 采用改进的城市夜间灯光指数法VBANUI提取城市建成区， 的计算公式如下：

[0048]

[0049] 式中， 为夜间灯光遥感图像；

[0050] 在完成 指标的计算后，令

[0051]

[0052] 式中， 为每个城市斑块的 最大值。

[0053] 本发明的有益效果是：通过结合NPP/VIIRS夜间灯光数据和高分辨率影像提供了一种城市建成区的提取方法，该提取方法可以有效减少NPP/VIIRS夜间灯光数据溢出效应的影响，尤其在水体、建筑区内的植被区域等的夜间灯光亮度溢出现象，弥补当前已有方法的不足。同时本发明构建一种语义分割网络模型：利用夜间灯光数据提取不同形态城市建成区的最佳阈值，获取加权交叉熵。

附图说明

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

[0055] 图1是本发明的城市建成区提取方法的流程示意图；

[0056] 图2是本发明的语义分割网络模型的构建过程示意图；

[0057] 图3是本发明的城市建成区的最优阈值获取过程示意图。

具体实施方式

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

[0059] 本实施例提供了一种基于夜间灯光数据与高分辨率影像的城市建成区提取方法，如图1所示，包括：

[0060] S1：采集同一区域的NPP/VIIRS夜间灯光遥感图像和高分辨率全色遥感图像，并分别对这两种图像进行预处理。

[0061] 1.对NPP/VIIRS夜间灯光遥感图像进行预处理。

[0062] 首先用ArcMap完成NPP/VIIRS夜间灯光遥感图像的裁剪、重采样和重投影，裁剪、重采样和重投影属于常见的遥感图像预处理，过程细节不再赘述。NPP/VIIRS夜间灯光遥感图像的原始数据采用WGS1984坐标系，将其转换为Albers等积投影坐标系。然后将其重采样到0.5km×0.5km的网格中，本步骤主要是为了避免由坐标系引起的图像失真并简化图像区域的计算。

[0063] 对NPP/VIIRS夜间灯光遥感图像进行年度均值图像的合成。在下载时月合成数据出现5、6月份缺失的情况，这可能是因为极光等污染光源的影响，且这种现象在每年的月合成数据中基本都会出现，因此采取了剔除这两个月的数据然后把剩余10个月数据求平均的方法来合成年度均值图像。合成公式如下所示：

[0064]

[0065] 式中： 表示第 个月的夜间灯光总量， 表示年平均夜间灯光总量。

[0066] 对NPP/VIIRS夜间灯光遥感图像进行不稳定光源和背景噪声消除、极值消除，把小于零的值设为0.001（近似为0，对统计分析没影响），把背景值设为0；把大于235的数值设为235（根据推理经验值），以删除一些异常值。

[0067] 2.对高分辨率全色遥感图像进行预处理。

[0068] 首先对高分辨率全色遥感图像进行辐射校正、大气校正，然后利用小波变换算法对校正后的图像进行去云处理，对图像进行直方图均衡化拉伸，其作用是通过拉伸像素强度分布范围来增强图像对比度。辐射校正、大气校正和直方图均衡化拉伸属于常见的遥感图像预处理，过程细节不再赘述。

[0069] S2：构建语义分割网络模型。

[0070] 构建的网络模型采用基于深度学习的语义分割方法，具体有Unet、SegNet和FCN等模型，采用编码器‑解码器结构。

[0071] 如图2所示，本实施例所述语义分割网络模型的构建过程具体为：

[0072] S201：对预处理后的高分辨率全色遥感图像样本内的建筑区域与其它区域进行分割标注。

[0073] 对高分辨率全色遥感图像样本进行标注，像素位置为建筑区域的标注为0，其它标注为1，最终得到一张像素值包含0,1的标注图像。

[0074] S202：将同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率进行调整，使同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像的分辨率相同。

[0075] 对同一区域的夜间灯光遥感图像和高分辨率全色遥感图像样本进行分块处理并调整分辨率。由于遥感图像数据分辨率较大，输入网络时需要进行分块处理，例如一个10000*10000的图像，可以分块成400块500*500大小的。同时由于两种图像分辨率不同，需要采用图像插值的方法把两种图像调到相同分辨率，常用的图像插值方法有：最邻近元法、双线性内插法和三次内插法。本实施例图像插值方法采用三次内插法，该方法利用三次多项式求逼近理论上最佳插值函数。

[0076] S203：利用同一区域的高分辨率全色遥感图像样本和夜间灯光遥感图像中每行每列的像素及其权重构建交叉熵损失函数。

[0077] 语义分割网络模型采用交叉熵，它是对给定随机变量或事件集的两个概率分布之间的差异的度量，分类任务采用softmax激活函数+交叉熵损失函数，softmax激活函数将一个向量进行“归一化”成概率分布的形式，再采用交叉熵损失函数计算Loss。

[0078] 交叉熵损失函数 的具体表达为：

[0079]

[0080] 其中M、N分别为图像的宽、高， 表示第i行第j列像素的权重，加权交叉熵基于图像位置进行加权，可以使得网络更注意VBANUI大的区域，且可以学习到DN与建成区的关系，提高提取城市建成区的准确性。 表示语义分割网络预测的第i行第j列像素的概率。 表示语义分割网络的第i行第j列像素的标签值。

[0081] 上述加权交叉熵是基于图像位置进行加权，可以使得网络更注意VBANUI大的区域，且可以学习到DN与建成区的关系，提高城市建成区的提取准确性，每个像素的权重的获取方法具体如下：

[0082] （1）将夜间灯光遥感图像分为三个区域：建成区、非建成区和过渡区。将亮度辐射值大于c的像元作为建成区，像元亮度辐射值介于a和c的作为过渡区，亮度辐射值低于a的像元作为非建成区。将整幅夜间灯光遥感图像输入做如下处理：

[0083] ①对建成区和非建成区进行中值滤波去噪，各选取一个中心邻域，然后对中心邻域内各个像素的灰度值按大小进行排序，选取排序序列中的中值作为该中心邻域中心像素点的灰度值；

[0084] ②对建成区和非建成区进行形态重构，选择3维矩阵 作为结构元，对得到的图像进行闭运算，此操作可以较好地提高图像的轮廓清晰度，进而得到三个分区的轮廓信息；

[0085] ③利用Sobel算子计算梯度数据，该算子是一种1阶微分算子，它利用像素邻近区域的梯度值来计算1个像素的梯度，然后根据一定的绝对值来取舍。

[0086] ④利用分水岭分割算法提取建成区，通常把梯度图像作为输入图像。

[0087] 通过形态重构和标记得到二值图，其像素值记为 ，有效地抑制了提取建成区的过程中存在欠分割的问题。

[0088] （2）对预处理后的夜间遥感图像采用周长突变检测法，由于城市建成区是聚合人类活动的整体，且NPP/VIIRS夜间灯光遥感图像分辨率较低，所以城市建成区内部应该很少有破碎斑块，且其夜间灯光数值应该比非城市建成区和城郊都更高。因此，随着设置的阈值不断增大，划归城市建成区的像元不断减少，城市建成区的周长相应也变小。但是，当阈值增大到一定程度时，城市建成区内部开始破碎，进而导致城市建成区周长没有变小，反而骤然增大。正确的阈值应该能够保证提取出来的城市结构的整体性，使之鲜少存在破碎区，所以将突变点之前的灯光值识别为最优阈值。

[0089] 如图3所示，获取最优阈值的具体实现步骤如下：

[0090] 1.设置初始阈值为 ，将辐射值不小于 的像元划分为城市建成区像元；

[0091] 2.将连续的城市像元聚合成城市斑块，计算城市斑块的总数量，记为Q；

[0092] 3.从1到Q标识城市斑块，计算每一个城市斑块的周长，并求和；

[0093] 4.城市斑块的周长计算公式如下：

[0094]

[0095] 式中， 表示取最优阈值时第n个城市斑块的周长，Q表示取最优阈值时城市斑块的数量。求周长是斑块连通域的一个属性，通过OpenCV直接调用，这里不再赘述。

[0096] 5.回到第一步，将阈值 以0.01（根据推理经验值）为间隔，迭代增加；观察并统计随着阈值的增加，城市斑块的周长的分布情况。

[0097] 6.找出最优阈值。找到城市斑块的周长突然剧烈增加时对应的辐射值，这个辐射值减掉0.01就是最优阈值；

[0098] 7.基于城市斑块的高斯赋权，采用高斯核函数，得到 ，高斯核函数可以将有限维数据映射到高维空间，高斯核函数圆心为每个斑块的质心，其定义为：

[0099]

[0100] 其中，是带宽（根据经验值，在二维情况下 =1），控制径向作用范围，换句话说，控制高斯核函数的局部作用范围。在图形上，正态分布是一种钟形曲线，越接近中心，取值越大，越远离中心，取值越小。 表示每个城市斑块中像素的坐标，当 和 的欧式距离处于某一个区间范围内的时候，假设固定 ， 随 的变化而变化，且变化程度相当显著。

[0101] 8.提取每一个城市斑块内的年平均夜间灯光总量，最终每个城市斑块都得到一个权重，则每个城市斑块的平均灯光强度计算公式如下：

[0102] =

[0103] 表示取最优阈值时第 个城市斑块的年平均夜间灯光总量。

[0104]

[0105] 表示取最优阈值时所有城市斑块的周长之和， 表示每个城市斑块的平均灯光强度， 表示高斯核函数。

[0106] （3）从高分辨率全色遥感图像样本中提取NDVI和NDBI，其中NDVI为归一化植被指数，NDBI为归一化建筑指数，NDVI和NDBI的表达式如下：

[0107]

[0108]

[0109] 其中， 为近红外波段，为红外波段， 为中红外波段。

[0110] 采用改进的城市夜间灯光指数法VBANUI提取城市建成区， 的计算公式如下：

[0111]

[0112] 其中， 为夜间灯光遥感图像，  为归一化植被指数， 为归一化建筑指数。

[0113] NDVI取值为[‑1,1]，NDVI>0为植被覆盖区，NDVI<0为非植被覆盖区。表达式(1‑NDVI)则表征了较大的非植被覆盖区权重在城市核心地区，(1‑NDVI)与NTL的结合可减少城市核心区夜间灯光的饱和现象，并增加城市内部核心区变化特征的快速识别。

[0114] 其次，NDBI取值范围也在[‑1，1]之间。研究表明NDBI正值为城市用地区域，NDBI负值为非城市用地区域。融入NDBI后，在水体对城市建成区混淆或者NDVI和NTL有溢出效应时，都可以有效控制并减弱其影响。

[0115] 在完成 指标的计算后，令

[0116]

[0117] 其中， 为每个城市斑块的 最大值。

[0118] 至此，对于图像中的每个像素都有三个初始权重，分别为 。

[0119] （4）由上述公式计算可知：

[0120]

[0121]

[0122]

[0123] 最终得到每个像素的权重的表达式为：

[0124] =

[0125] 将 代入上述交叉熵损失函数中进行计算。

[0126] S204：利用构建的交叉熵损失函数对语义分割网络模型进行监督训练，直至语义分割网络收敛为止，得到训练好的语义分割网络。

[0127] S3：将预处理后的同一区域的NPP/VIIRS夜间灯光遥感图像和高分辨率全色遥感图像输入到训练好的语义分割网络模型中，完成城市建成区的提取。

[0128] 将上述两种图像输入编码器，经过编码器特征提取，输出为特征图，解码器输入为特征图，经过解码器的上采样和拟合，得到目标分割图。

[0129] 对语义分割网络进行优化，常用优化器有：SGD、Adam和Lookahead优化器，本实施例优选使用Lookahead优化器。

[0130] 综上所述，本实施例提出的方法能够结合夜间灯光数据和高分辨率数据，提出加权交叉熵的方法，再引入NDVI和NDBI指标，有效减少了灯光亮度与NDVI溢出效应的影响，达到提高城市建成区提取准确性的目的。

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