利用夜光遥感影像进行城市化监测

王海军，孔祥冬，武克军，张 勃，卓 壮，任 琪

1.成都理工大学 工程技术学院，四川 乐山 614000

2乐山师范学院 四川旅游发展研究中心，四川 乐山 614000

3.西北师范大学 地理与环境学院，兰州 730072

1 前言

成渝城市圈是由成都和重庆两个大型城市为双核心，毗邻的中小城市为卫星城的重要城市群，也是中国西南地区城市化程度较高的地区。尤其是随着成渝、成绵乐高速铁路、公路网的完善，使得路网沿线和双核心辐射区城市化进程明显加快。城市化进程的提速，势必对周围的生态化环境、耕地资源等，产生胁迫和侵蚀。城市化扩展速度、方向不同，对于比邻区的生态胁迫程度也不同，因此对于成渝地区城市化进程的空间扩展情况的研究十分必要。

传统城市化监测方法主要有两种：第一，采用多期较高分辨率遥感数据进行自动分类、统计，从而定性定量地分析城市扩展情况。由于城市下垫面地物类型复杂，存在“同物异谱”和“同谱异物”的问题，很难得到较高精度的分类结果；第二，基于航拍数据，利用人工目视解译的方法，对城市地物分类，该种方法分类精度高，但是工作量巨大，时效性不好。鉴于上述两点，本文采用间接的方法来分析和监测城市扩展情况，利用城市夜晚灯光亮度（DMSP/OLS）来反映城市分布情况，长时间序列的亮度数据监测城市扩展情况。该方法的原理是：卫星传感器在夜晚开机成像，拍摄到城市夜晚的灯光亮度和分布情况，数据亮度值范围0～63。亮度值越高表明该区域是城市分布的概率就越高，反之亦然。该数据可以真实反映出区域城市分布情况。长时间序列的DMSP/OLS数据在国外已经得到了成熟的检验和应用[1-7]，近几年DMSP/OLS理论和应用研究也逐步受到国内学者的关注。

目前，对于DMSP/OLS图像上的城市区域提取方法主要采用的是灰度阈值法，通过设定阈值，大于该阈值的即为城市区，小于该阈值则为非城市区。该方法将图像的最小单位设定为像元（pixel），通过对于单个像元的DN值进行分析，来将像元进行归类统计。往往会对非城市区的高亮度像素归为城市，存在错分的现象。为解决此问题，本文采用梯度分割+灰度阈值的方法进行城市区域的识别与分类。该方法的主要原理是：首先基于图像的纹理、结构、灰度的聚类特点，将图像的像元（pixel）通过梯度分割法分割成大小不同的对象（object）。再通对图像的平均DN值进行统计，计算出城市区的阈值范围。从而将城市区从全景图像上识别出来。

2 数据获取与研究方法

2.1 数据获取

夜光遥感数据为美国国家地理中心（NGDC）在2013年发布的第5版的非福射定标夜间平均灯光强度数据（DMSP/OLS），该数据为全球夜间灯光影像，空间分辨率1 km，DN 值范围0～63，数据目前有F10、F12、F14、F15、F16、F18，6颗卫星探测到的1992—2013年共21年的数据。此夜光强度数据与对应年份的遥感影像数据相比，具有反映城市分布，更加直观、数据处理量小、数据时效性好。DMSP/OLS在夜间工作，能够探测到城市灯光、小规模居民点、交通道路、车流等发出的低强度灯光，可以综合反映人类活动信息；且与NOAA/AVHRR的空间分辨率、时间分辨率相当，适合动态监测大尺度城镇扩展；交通网、行政区划、地形等数据来源于国家基础地理信息中心，交通网数据用于与城市灯光分布点匹配，地形数据与行政区划数据用于夜光遥感数据分割与分类。

2.2 研究方法

2.2.1 图像定标

DMSPOLS夜光图像据是由多颗卫星获取，由于OLS传感器未做星上定标，因此F10（1992—1994年）、F12（1994—1999年）、F14（1997—2003年）、F15（2000—2007年）、F16（2004—2009年）、F18（2010—2013年）多颗卫星获取的图像缺乏可比性[8-9]。因此本文采用图像间相互校正法[8]对1992—2013年21期夜光图像进行校正处理。F16（2007年）图像DN累计值最高，选取其作为参考数据。结合本文实验区内绵阳市城市发展较为平稳，DN值变化幅度较小的特点，选取绵阳市地区作为定标区。构建定标区内参考图像与待矫正图像DN值一元二次回归模型并计算回归系数。以此回归系数采用式（1）对1992—2013年的数据进行校正（校正采用Arcgis10空间分析模块实现）。

其中IMG与IMG0（未定标）分别是校正前后的图像，a、b、c为回归系数。

2.2.2 图像识别

（1）图像分割

图像分割是图像识别的基础，图像分割目的在于将图像上不同的目标分割成对象（object）。图像分割依据图像本身的颜色、纹理、灰度、结构的差异而进行的。本文通过对目前比较成熟的图像处理平台（Matlab、Envi、e-Cognition）进行分析，选取从分割数据纹理、结果最适合的两种方法：多尺度分割（Multi-resolution segmentation）和梯度分割（Contrast filter segmentation）[10-11]。多尺度分割技术源于医学图像识别，是一种自下而上逐级合并技术。图像分割最大尺度（scale）即为图像本身（1行×1列），分割最小尺度（scale）为一个像素并且图像被分割成（n行×m列），实际应用中往往是先将图像分成像素级，在此基础上逐级合并，直至合并后对象内部异质性最小，合并即停止。多尺度分割多用于多光谱和高光谱图像的分割。梯度分割技术是对图像灰度逐行检测，检测像素DN值变化梯度，如果发生明显改变则视为图像上不同斑块的边缘区，全景图像检测后，将全景图像再次按照像素级别分割。梯度分割多用于全色波段图像的分割，是一种自上而下的图像分割技术。

（2）图像分类

图像分类是将分割后的图像对象（object）定义到指定的分类系统中的过程。本文分类的目的在于将夜光图像上的城市区和非城市区分开，空间显示出1992—2013年期间成渝地区城市化的进程，因此本文建立的分类系统为Urban和Un-urban。比较了最邻近分类和灰度阈值[12-13]分类两种分类算法。阈值分类更适合全色波段图像分类。通过对分割后的图像灰度值的数据挖掘和训练，确定了不同年份图像城市区阈值范围，进而完成图像分类。

（3）精度验证

图像上单个对象的分类精度，采用该对象（城市区域）的分类统计面积与对象实际面积只差再与实际面积进行比值计算。对于全景图像分类精度采用每个对象的面积差累计求和后与实际面积累计求和进行比值[14-15]，计算原理如下公式所示：

其中，Ai表示图像上单个对象分类精度；An表示全景图像上目标分类总精度；si表示城市区实际的面积；soi表示图像上城市区分类统计的面积；i→m表示图像上分类区的目标个数。

3 成渝地区城市化进程检测

3.1 夜光图像定标

根据2.2.1小节所阐述的方法，构建定标区内参考图像与待矫正图像DN值一元二次回归模型并计算回归系数。通过构建参考图像与待矫正图像的回归模型（图1），统计分析出参考图像F16-2007，与带矫正图像1992—2013年（除2007年）回归模型的数学方程，并对其a、b、c三个参数进行计算，建立起回归系数表。利用公式（1）在Arcgis10空间分析模块中，采用所建立的函数关系对各期图像进行运算处理，得到的图像即为定标后的图像。因为该数据均是以2007年图像为参考数据，因此各期图像具备可比性。可以用于后期图像分类后的城市化进程的叠加分析。

3.2 图像城市区识别

采用多尺度与梯度分割方法，对成渝地区定标后的夜光图像进行分割实验。首先采用多尺度分割方法进行分割，并且对分割参数：尺度（scale）、形状（Shape）、紧致（Compactness）进行设定，原则是自下而上，然后逐级合并。三项参数分别设置为：30、0.3、0.3，分割效果如图2（a）、图3（a）所示，对于小目标（small object）分割效果较好，可以完全从背景图像中识别出来，但是大的对象（big object）分割效果则不理想，存在碎化并且识别不完全的现象。在此基础上，对分割结果进行逐级合并，直至参数值为60、0.5、0.6时，大目标分割效果较为理想，但是小目标存在漏分的现象（图2（b）、图3（b））。对同一幅图像，采用梯度分割方法进行分割实验，分割参数 object size=10 ，分割结果如图2（c）、图3（c）所示，对于成都和重庆两个城市及毗邻地区的中小城市目标分割都较为理想，不存在错分和漏分的问题，表明梯度分割方法适合该图像的分割，并且效果理想，可用于该区域的图像分割。通过梯度滤波的分割方法，将全景图像上的目标划分成不同的对象。采用阈值分类的方法，通过阈值实验，结果阈值在15～63范围内区域为城市区。为了对分类范围的验证，对2013年图像上眉山、乐山、资阳、自贡、泸州、南充、广安、万州区、涪陵区，单个图像对象的分类精度采用公式（2）进行精度评价，结果Ai值均小于0.071，同时对全景图像目标的总分类精度进行计算，结果An为0.052，表明对于图像分类的面积与真实面积偏差很小，说明分类结果是可信的。同时采用百度地图数据与分类结果数据进行叠加显示，城区边界范围吻合较好，说明夜光图像分类结果精度较理想，可以用于城市变化监测分析。

图1 1992—2013年图像校正模型及参数（部分回归模型散点图）

图2 成都地区夜光图像分割结果（2013年）

图3 重庆地区夜光图像分割结果（2013年）

3.3 成渝地区城市化进程分析

通过图像分割、分类提取出城市区的1992—2013年每年矢量边界。对每年的矢量边界进行叠加分析，显示出城市区在不同年份扩展情况（如图4所示）。从图上可以看出，成都和重庆以原有的城市中心区向外扩展，1992—2003年来城市扩展速度较慢，而从2004年以来扩展速度很快，扩展的范围较大。成都及其比邻区到2013年时，形成了以成都为核心，以绵阳-德阳-成都-眉山-乐山为轴线的发展趋势。扩展方向主要向东北和西南。结合成都及其比邻区的重要交通网（图5（a））发现，该区域已形成明显的以成都为中心城和以绵阳、遂宁、内江、自贡、乐山、雅安为卫星城的环状城市发展结构。同时中心城与卫星城都有重要的交通网络连接，形成了以成都市外环、S106、卫星城间的交通线等，三环网络结构。重庆地区在2004—2013年城市扩展速度较快，但是扩展区域小于成都地区。重庆向四周扩展的差异性较小，以向西为主。结合重庆周围的交通网（图5（b））可见，重庆城市发展结构特点呈现星型形状，以重庆市为中心向四周辐射。而辐射区域主要是面向达州、南充、涪陵、泸州。与成都相比，没有形成明显的环状城市发展与交通网络结构。

图4 成渝地区1992—2013年城市化进程空间分布

图5 成渝地区1992—2013年城市化结构空间分布

4 结论

采用夜光遥感图像对成渝城市群近21年的城市化进行监测分析。改变了传统基于多光谱遥感数据直接进行监测的方法，而是利用城市灯光强度在长时间序列内的辐射变化，反映城市扩展情况，为今后的城市发展与监测提供一种新的思路与理念。同时在对于夜光图像城市区识别方法上，本文对传统方法进行了创新，采用了梯度分割和灰度阈值分类的组合方式，提高了对图像城市区的识别精度。研究结果表明，对于定标处理后图像的识别精度An达到了0.052，提取的城市区面积范围与百度地区城市区面积范围很好的吻合，该方法大大提高了城市化监测的时效性；成渝地区城市化扩展呈现不同特点，成都以及毗邻区主要呈现以绵阳-德阳-成都-眉山-乐山为主轴的东北西南向扩展特点，而城市扩展提速是从2003年以后开始，近10年该地区城市扩展范围增加明显。而重庆地区则是重庆为核心向四周扩展，没有明显的扩展较快的方向；近21年成都和重庆的城市扩展呈现了独特的扩展结构，成都及其毗邻区呈现明显的三层环状结构，第一层为成都外环、第二层为s106省道，将彭州、乐至、中江、眉山、都江堰等城市进行了串联、第三层为绵阳、遂宁、内江、自贡、乐山。而重庆的扩展结构呈现星型，以重庆为中心向外辐射，辐射方向主要是西北和西南，尤其是向成都方向。

[1]Liu Z F，He C Y，Zhang Q F，et al.Extracting the dynamics of urban expansion in China using DMSP-OLS nighttime light data from 1992 to 2008[J].Landscape and Urban Planning，2012，106（1）：62-72.

[2]Wu J，He S，PENG J，et al.Inter calibration of DMSPOLS nighttime light data by the invariant region method[J].International Journal of Remote Sensing，2013，34（20）：7356-7368.

[3]Pandey B，Joshi P K，Seto K C.Monitoring urbanization dynamics in India using DMSP/OLS night time lights and SPOT-VGT data[J].International Journal of Applied Earth Observation and Geo-information，2013，23：49-61.

[4]Li X，Xu H M，Chen X L，et al.Potential of NPP-VIIRS nighttime light imagery for modeling the regional economy of China[J].Remote Sensing，2013，5（6）：3057-3081.

[5]Forbes D J.Multi-scale analysis of the relationship between economic statistics and DMSP-OLS night Light images[J].GIs Science and Remote Sensing，2013，50（5）：483-499.

[6]Keola S，Andersson M，Hall O.Monitoring economic development from space：Using night light and land cover data to measure economic growth[J].World Development，2015，66：322-334.

[7]Min B，Gaba K M，Sarr O F，et al.Detection of rural electrification in Africa using DMSP-OLS night lights imagery[J].International Journal of Remote Sensing，2013，34（22）：8118-8141.

[8]Elvidge C D，Ziskin D，Baugh K E，et al.A fifteen year record of global natural gas flaring derived from satellite data[J].Energies，2009，2（3）：595-622.

[9]Baugh K，Elvidge C，Ghosh T，et al.Development of a 2009 stable lights product using DMSP-OLS data[C]//Proceedings of the Asia Pacific Advanced Network，2010.

[10]Benz U C，Hofmann P，Willhauck G，et al.Multi-resolution，object-oriented fuzzy analysis of remote sensing data for GIS-Ready information[J].SPRS Journal of Photogrammetry&Remote Sensing，2004，58：239-258.

[11]Sreejini K S，Govindan V K.Improved multi-scale matched filter for retina vessel segmentation using PSO algorithm[J].Egyptian Informatics Journal，2015，16：253-260.

[12]Baraldi P，Cannarile F，Maio F D，et al.Hierarchical k-nearest neighbors classification and binary differential evolution for fault diagnostics of automotive bearings operating under variable conditions[J].Engineering Applications of Artificial Intelligence，2016，56：1-13.

[13]Issac A，Sarathi M P，Dutta M K.An adaptive threshold based image processing technique forimproved glaucoma detection and classification[J].Computer Methods and Programs in Biomedicine，2015，122（2）：229-244.

[14]Wang Haijun，Kong Xiangdong，Zhang Bo.The simulation of lucc based on logistic-CA-Markov model in Qilian Mountain area，China[J].Sciences in Cold and Arid Regions，2016，8（4）：350-358.

[15]Wang Haijun，Dai Shengpei，Huang Xiaobin.The remote sensing monitoring analysis based on object-oriented classification method[C]//ChineseConferenceon Image and Graphics Technologies，2013：92-101.