王海軍，孔祥冬，武克軍，張 勃，卓 壯，任 琪

1.成都理工大學 工程技術學院，四川 樂山 614000

2樂山師范學院 四川旅游發(fā)展研究中心，四川 樂山 614000

3.西北師范大學 地理與環(huán)境學院，蘭州 730072

1 前言

成渝城市圈是由成都和重慶兩個大型城市為雙核心，毗鄰的中小城市為衛(wèi)星城的重要城市群，也是中國西南地區(qū)城市化程度較高的地區(qū)。尤其是隨著成渝、成綿樂高速鐵路、公路網(wǎng)的完善，使得路網(wǎng)沿線和雙核心輻射區(qū)城市化進程明顯加快。城市化進程的提速，勢必對周圍的生態(tài)化環(huán)境、耕地資源等，產(chǎn)生脅迫和侵蝕。城市化擴展速度、方向不同，對于比鄰區(qū)的生態(tài)脅迫程度也不同，因此對于成渝地區(qū)城市化進程的空間擴展情況的研究十分必要。

傳統(tǒng)城市化監(jiān)測方法主要有兩種：第一，采用多期較高分辨率遙感數(shù)據(jù)進行自動分類、統(tǒng)計，從而定性定量地分析城市擴展情況。由于城市下墊面地物類型復雜，存在“同物異譜”和“同譜異物”的問題，很難得到較高精度的分類結(jié)果；第二，基于航拍數(shù)據(jù)，利用人工目視解譯的方法，對城市地物分類，該種方法分類精度高，但是工作量巨大，時效性不好。鑒于上述兩點，本文采用間接的方法來分析和監(jiān)測城市擴展情況，利用城市夜晚燈光亮度（DMSP/OLS）來反映城市分布情況，長時間序列的亮度數(shù)據(jù)監(jiān)測城市擴展情況。該方法的原理是：衛(wèi)星傳感器在夜晚開機成像，拍攝到城市夜晚的燈光亮度和分布情況，數(shù)據(jù)亮度值范圍0～63。亮度值越高表明該區(qū)域是城市分布的概率就越高，反之亦然。該數(shù)據(jù)可以真實反映出區(qū)域城市分布情況。長時間序列的DMSP/OLS數(shù)據(jù)在國外已經(jīng)得到了成熟的檢驗和應用[1-7]，近幾年DMSP/OLS理論和應用研究也逐步受到國內(nèi)學者的關注。

目前，對于DMSP/OLS圖像上的城市區(qū)域提取方法主要采用的是灰度閾值法，通過設定閾值，大于該閾值的即為城市區(qū)，小于該閾值則為非城市區(qū)。該方法將圖像的最小單位設定為像元（pixel），通過對于單個像元的DN值進行分析，來將像元進行歸類統(tǒng)計。往往會對非城市區(qū)的高亮度像素歸為城市，存在錯分的現(xiàn)象。為解決此問題，本文采用梯度分割+灰度閾值的方法進行城市區(qū)域的識別與分類。該方法的主要原理是：首先基于圖像的紋理、結(jié)構(gòu)、灰度的聚類特點，將圖像的像元（pixel）通過梯度分割法分割成大小不同的對象（object）。再通對圖像的平均DN值進行統(tǒng)計，計算出城市區(qū)的閾值范圍。從而將城市區(qū)從全景圖像上識別出來。

2 數(shù)據(jù)獲取與研究方法

2.1 數(shù)據(jù)獲取

夜光遙感數(shù)據(jù)為美國國家地理中心（NGDC）在2013年發(fā)布的第5版的非福射定標夜間平均燈光強度數(shù)據(jù)（DMSP/OLS），該數(shù)據(jù)為全球夜間燈光影像，空間分辨率1 km，DN 值范圍0～63，數(shù)據(jù)目前有F10、F12、F14、F15、F16、F18，6顆衛(wèi)星探測到的1992—2013年共21年的數(shù)據(jù)。此夜光強度數(shù)據(jù)與對應年份的遙感影像數(shù)據(jù)相比，具有反映城市分布，更加直觀、數(shù)據(jù)處理量小、數(shù)據(jù)時效性好。DMSP/OLS在夜間工作，能夠探測到城市燈光、小規(guī)模居民點、交通道路、車流等發(fā)出的低強度燈光，可以綜合反映人類活動信息；且與NOAA/AVHRR的空間分辨率、時間分辨率相當，適合動態(tài)監(jiān)測大尺度城鎮(zhèn)擴展；交通網(wǎng)、行政區(qū)劃、地形等數(shù)據(jù)來源于國家基礎地理信息中心，交通網(wǎng)數(shù)據(jù)用于與城市燈光分布點匹配，地形數(shù)據(jù)與行政區(qū)劃數(shù)據(jù)用于夜光遙感數(shù)據(jù)分割與分類。

2.2 研究方法

2.2.1 圖像定標

DMSPOLS夜光圖像據(jù)是由多顆衛(wèi)星獲取，由于OLS傳感器未做星上定標，因此F10（1992—1994年）、F12（1994—1999年）、F14（1997—2003年）、F15（2000—2007年）、F16（2004—2009年）、F18（2010—2013年）多顆衛(wèi)星獲取的圖像缺乏可比性[8-9]。因此本文采用圖像間相互校正法[8]對1992—2013年21期夜光圖像進行校正處理。F16（2007年）圖像DN累計值最高，選取其作為參考數(shù)據(jù)。結(jié)合本文實驗區(qū)內(nèi)綿陽市城市發(fā)展較為平穩(wěn)，DN值變化幅度較小的特點，選取綿陽市地區(qū)作為定標區(qū)。構(gòu)建定標區(qū)內(nèi)參考圖像與待矯正圖像DN值一元二次回歸模型并計算回歸系數(shù)。以此回歸系數(shù)采用式（1）對1992—2013年的數(shù)據(jù)進行校正（校正采用Arcgis10空間分析模塊實現(xiàn)）。

其中IMG與IMG0（未定標）分別是校正前后的圖像，a、b、c為回歸系數(shù)。

2.2.2 圖像識別

（1）圖像分割

圖像分割是圖像識別的基礎，圖像分割目的在于將圖像上不同的目標分割成對象（object）。圖像分割依據(jù)圖像本身的顏色、紋理、灰度、結(jié)構(gòu)的差異而進行的。本文通過對目前比較成熟的圖像處理平臺（Matlab、Envi、e-Cognition）進行分析，選取從分割數(shù)據(jù)紋理、結(jié)果最適合的兩種方法：多尺度分割（Multi-resolution segmentation）和梯度分割（Contrast filter segmentation）[10-11]。多尺度分割技術源于醫(yī)學圖像識別，是一種自下而上逐級合并技術。圖像分割最大尺度（scale）即為圖像本身（1行×1列），分割最小尺度（scale）為一個像素并且圖像被分割成（n行×m列），實際應用中往往是先將圖像分成像素級，在此基礎上逐級合并，直至合并后對象內(nèi)部異質(zhì)性最小，合并即停止。多尺度分割多用于多光譜和高光譜圖像的分割。梯度分割技術是對圖像灰度逐行檢測，檢測像素DN值變化梯度，如果發(fā)生明顯改變則視為圖像上不同斑塊的邊緣區(qū)，全景圖像檢測后，將全景圖像再次按照像素級別分割。梯度分割多用于全色波段圖像的分割，是一種自上而下的圖像分割技術。

（2）圖像分類

圖像分類是將分割后的圖像對象（object）定義到指定的分類系統(tǒng)中的過程。本文分類的目的在于將夜光圖像上的城市區(qū)和非城市區(qū)分開，空間顯示出1992—2013年期間成渝地區(qū)城市化的進程，因此本文建立的分類系統(tǒng)為Urban和Un-urban。比較了最鄰近分類和灰度閾值[12-13]分類兩種分類算法。閾值分類更適合全色波段圖像分類。通過對分割后的圖像灰度值的數(shù)據(jù)挖掘和訓練，確定了不同年份圖像城市區(qū)閾值范圍，進而完成圖像分類。

（3）精度驗證

圖像上單個對象的分類精度，采用該對象（城市區(qū)域）的分類統(tǒng)計面積與對象實際面積只差再與實際面積進行比值計算。對于全景圖像分類精度采用每個對象的面積差累計求和后與實際面積累計求和進行比值[14-15]，計算原理如下公式所示：

其中，Ai表示圖像上單個對象分類精度；An表示全景圖像上目標分類總精度；si表示城市區(qū)實際的面積；soi表示圖像上城市區(qū)分類統(tǒng)計的面積；i→m表示圖像上分類區(qū)的目標個數(shù)。

3 成渝地區(qū)城市化進程檢測

3.1 夜光圖像定標

根據(jù)2.2.1小節(jié)所闡述的方法，構(gòu)建定標區(qū)內(nèi)參考圖像與待矯正圖像DN值一元二次回歸模型并計算回歸系數(shù)。通過構(gòu)建參考圖像與待矯正圖像的回歸模型（圖1），統(tǒng)計分析出參考圖像F16-2007，與帶矯正圖像1992—2013年（除2007年）回歸模型的數(shù)學方程，并對其a、b、c三個參數(shù)進行計算，建立起回歸系數(shù)表。利用公式（1）在Arcgis10空間分析模塊中，采用所建立的函數(shù)關系對各期圖像進行運算處理，得到的圖像即為定標后的圖像。因為該數(shù)據(jù)均是以2007年圖像為參考數(shù)據(jù)，因此各期圖像具備可比性?？梢杂糜诤笃趫D像分類后的城市化進程的疊加分析。

3.2 圖像城市區(qū)識別

采用多尺度與梯度分割方法，對成渝地區(qū)定標后的夜光圖像進行分割實驗。首先采用多尺度分割方法進行分割，并且對分割參數(shù)：尺度（scale）、形狀（Shape）、緊致（Compactness）進行設定，原則是自下而上，然后逐級合并。三項參數(shù)分別設置為：30、0.3、0.3，分割效果如圖2（a）、圖3（a）所示，對于小目標（small object）分割效果較好，可以完全從背景圖像中識別出來，但是大的對象（big object）分割效果則不理想，存在碎化并且識別不完全的現(xiàn)象。在此基礎上，對分割結(jié)果進行逐級合并，直至參數(shù)值為60、0.5、0.6時，大目標分割效果較為理想，但是小目標存在漏分的現(xiàn)象（圖2（b）、圖3（b））。對同一幅圖像，采用梯度分割方法進行分割實驗，分割參數(shù) object size=10 ，分割結(jié)果如圖2（c）、圖3（c）所示，對于成都和重慶兩個城市及毗鄰地區(qū)的中小城市目標分割都較為理想，不存在錯分和漏分的問題，表明梯度分割方法適合該圖像的分割，并且效果理想，可用于該區(qū)域的圖像分割。通過梯度濾波的分割方法，將全景圖像上的目標劃分成不同的對象。采用閾值分類的方法，通過閾值實驗，結(jié)果閾值在15～63范圍內(nèi)區(qū)域為城市區(qū)。為了對分類范圍的驗證，對2013年圖像上眉山、樂山、資陽、自貢、瀘州、南充、廣安、萬州區(qū)、涪陵區(qū)，單個圖像對象的分類精度采用公式（2）進行精度評價，結(jié)果Ai值均小于0.071，同時對全景圖像目標的總分類精度進行計算，結(jié)果An為0.052，表明對于圖像分類的面積與真實面積偏差很小，說明分類結(jié)果是可信的。同時采用百度地圖數(shù)據(jù)與分類結(jié)果數(shù)據(jù)進行疊加顯示，城區(qū)邊界范圍吻合較好，說明夜光圖像分類結(jié)果精度較理想，可以用于城市變化監(jiān)測分析。

圖1 1992—2013年圖像校正模型及參數(shù)（部分回歸模型散點圖）

圖2 成都地區(qū)夜光圖像分割結(jié)果（2013年）

圖3 重慶地區(qū)夜光圖像分割結(jié)果（2013年）

3.3 成渝地區(qū)城市化進程分析

通過圖像分割、分類提取出城市區(qū)的1992—2013年每年矢量邊界。對每年的矢量邊界進行疊加分析，顯示出城市區(qū)在不同年份擴展情況（如圖4所示）。從圖上可以看出，成都和重慶以原有的城市中心區(qū)向外擴展，1992—2003年來城市擴展速度較慢，而從2004年以來擴展速度很快，擴展的范圍較大。成都及其比鄰區(qū)到2013年時，形成了以成都為核心，以綿陽-德陽-成都-眉山-樂山為軸線的發(fā)展趨勢。擴展方向主要向東北和西南。結(jié)合成都及其比鄰區(qū)的重要交通網(wǎng)（圖5（a））發(fā)現(xiàn)，該區(qū)域已形成明顯的以成都為中心城和以綿陽、遂寧、內(nèi)江、自貢、樂山、雅安為衛(wèi)星城的環(huán)狀城市發(fā)展結(jié)構(gòu)。同時中心城與衛(wèi)星城都有重要的交通網(wǎng)絡連接，形成了以成都市外環(huán)、S106、衛(wèi)星城間的交通線等，三環(huán)網(wǎng)絡結(jié)構(gòu)。重慶地區(qū)在2004—2013年城市擴展速度較快，但是擴展區(qū)域小于成都地區(qū)。重慶向四周擴展的差異性較小，以向西為主。結(jié)合重慶周圍的交通網(wǎng)（圖5（b））可見，重慶城市發(fā)展結(jié)構(gòu)特點呈現(xiàn)星型形狀，以重慶市為中心向四周輻射。而輻射區(qū)域主要是面向達州、南充、涪陵、瀘州。與成都相比，沒有形成明顯的環(huán)狀城市發(fā)展與交通網(wǎng)絡結(jié)構(gòu)。

圖4 成渝地區(qū)1992—2013年城市化進程空間分布

圖5 成渝地區(qū)1992—2013年城市化結(jié)構(gòu)空間分布

4 結(jié)論

采用夜光遙感圖像對成渝城市群近21年的城市化進行監(jiān)測分析。改變了傳統(tǒng)基于多光譜遙感數(shù)據(jù)直接進行監(jiān)測的方法，而是利用城市燈光強度在長時間序列內(nèi)的輻射變化，反映城市擴展情況，為今后的城市發(fā)展與監(jiān)測提供一種新的思路與理念。同時在對于夜光圖像城市區(qū)識別方法上，本文對傳統(tǒng)方法進行了創(chuàng)新，采用了梯度分割和灰度閾值分類的組合方式，提高了對圖像城市區(qū)的識別精度。研究結(jié)果表明，對于定標處理后圖像的識別精度An達到了0.052，提取的城市區(qū)面積范圍與百度地區(qū)城市區(qū)面積范圍很好的吻合，該方法大大提高了城市化監(jiān)測的時效性；成渝地區(qū)城市化擴展呈現(xiàn)不同特點，成都以及毗鄰區(qū)主要呈現(xiàn)以綿陽-德陽-成都-眉山-樂山為主軸的東北西南向擴展特點，而城市擴展提速是從2003年以后開始，近10年該地區(qū)城市擴展范圍增加明顯。而重慶地區(qū)則是重慶為核心向四周擴展，沒有明顯的擴展較快的方向；近21年成都和重慶的城市擴展呈現(xiàn)了獨特的擴展結(jié)構(gòu)，成都及其毗鄰區(qū)呈現(xiàn)明顯的三層環(huán)狀結(jié)構(gòu)，第一層為成都外環(huán)、第二層為s106省道，將彭州、樂至、中江、眉山、都江堰等城市進行了串聯(lián)、第三層為綿陽、遂寧、內(nèi)江、自貢、樂山。而重慶的擴展結(jié)構(gòu)呈現(xiàn)星型，以重慶為中心向外輻射，輻射方向主要是西北和西南，尤其是向成都方向。

[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.