王振興,劉 東*,王 敏,于依萍
(1. 中國科學院大學資源與環境學院,北京 100049)
沿海城市是我國經濟發展的重要引擎,具有獨特的擴張特征與生態優勢。目前針對沿海城市擴張研究相對較少,探討沿海城市擴張及其生態環境效應有助于理解沿海城市擴張規律及其與生態要素之間的關系。青島市是我國東部濱海度假旅游城市,近幾十年來經濟社會快速發展,成為區域的重要發展極。本文以青島市為研究區,融合多源夜間燈光影像(DMSP/OLS 和NPP/VIIRS)提取城市建成區,分析近20 a(2000—2019年)青島城市擴張的時空特征;利用MODIS 數據,提取植被覆蓋、熱島強度,通過相關性分析,揭示城市擴張過程對生態環境的影響,以期為城市的規劃與生態環境的保護提供政策建議和參考。
青島市(119°30’~121°00’E,35°35’~7°09’N)地處我國山東半島東南部,瀕臨黃海,是山東省重要中心城市和國家計劃單列城市,也是“一帶一路”新亞歐大陸橋經濟走廊建設和海上合作中主要節點城市。青島市包含7 個市轄區和3 個縣級市,總面積約11293 km2。地形以海濱丘陵為主,屬于溫帶季風氣候,夏季無酷暑,冬季少嚴寒,降水充沛,四季分明。2019年底,全市常住總人口為949.98萬,生產總值達11741.31億元[1]。
研究數據包括夜光遙感影像(DMSP/OLS 和NPP/VIIRS)、Landsat TM影像、MODIS數據(MOD11A2和MOD13A2)、行政區劃數據及統計年鑒數據。DMSP/OLS 來源于美國國家地球物理數據中心(https://www.ngdc.noaa.gov/)發布的第四版穩定的年合成影像數據,分別選取F14 衛星(2000年),F15 衛星(2000年、2004年),F16 衛星(2004、2008年) 和F18 衛星(2012年)等6幅影像。DMSP/OLS影像空間分辨率約1 km×1 km,圖像灰度范圍在0~63。因缺乏輻射定標,數據間不具備一致性。NPP/VIIRS來源于科羅拉多礦業大學(https://eogdata.mines.edu/products/vnl/)提供的NPP月合成數據,空間分辨率約500 m×500 m,已完成星上輻射定標。本文選取2012、2016、2019年等多幅月產品,進行平均值運算得到年合成數據。
Landsat TM 來源于美國地質勘探局(https://www.usgs.gov/),空間分辨率30 m×30 m,選取與夜光影像相同年份影像(2000、2004、2008、2012、2016年和2019年),用于各年份城市建成區提取的空間驗證。
MOD11A2為8日合成陸地溫度產品,MOD13A2為16 d 合成植被指數產品,分辨率為1 km×1 km,來源于美國航空航天局官方網站(https://www.nasa.gov)。6-8 月份是植被生長旺盛和氣溫最高時期,但影像易受云、雨等天氣影響。通過質量控制波段剔除受云等因素影響的低質量像元,統計研究年份6-8月各日期青島市內可用像元數和地表溫度(LST)與植被指數(NDVI)的均值(圖1)。可用像元數表現為先下降后升高趨勢,7 月份達到谷值,8 月份逐漸恢復,可見6、7 月份受云雨天氣影響比較嚴重。LST 和NDVI 則均先升高后下降,分別在7 月3 日和8 月12 日達到峰值,隨后逐漸減小。綜合來看,MOD11A2數據的6月份影像可用像元數可基本覆蓋全域,同時LST達到峰值,利于分析熱島效應;MOD13A2 數據的8 月份影像可用像元數達到最大,且NDVI 在峰值附近,利于分析植被覆被變化。因此,本文分別選取2000、2004、 2008、 2012、 2016年和2019年6 月份MOD11A2 數據和8 月份的MOD13A2 數據分析熱島強度和植被覆蓋度。

圖1 6-8月MODIS可用像元數和NDVI、LST均值
行政區劃邊界來源于國家基礎地理信息中心(https://www.webmap.cn/)。建成區統計數據來源于《中國城市建設統計年鑒》。所有遙感影像經過預處理,將影像鑲嵌、裁剪為研究區域;坐標系轉換為WGS-84 UTM投影坐標系,像元分辨率通過最近鄰法重采樣至1 km×1 km。
DMSP/OLS 夜光影像由多個DMSP 衛星搭載著不同的OLS傳感器獲得,受傳感器輻射探測差異、器件老化等影響,存在相同位置下同一傳感器獲取的連續不同年份像元值之間存在異常,不同傳感器獲取的同一年份的數據之間差異明顯等數據未定標的問題[2-3]。本文選用不變區域目標法對影像進行輻射校正,以城市發展相對平穩的黑龍江雞西市為不變區域,選擇F15-2004年影像為參考影像,利用二次多項式模型進行相對輻射校正(公式1)。
式中,DNcal為校正后影像DN值;DN為校正前影像DN值;a、b 為回歸系數;c 為常數項。每幅待校正影像與參考影像的校正系數與擬合R2見表1。每幅影像校正完成后,將同一年份不同傳感器的影像取平均值,獲得該年份的合成產品,得到一致的DMSP/OLS序列數據。

表1 NPP/VIIRS與DMSP/OLS一致性校正系數與精度
對NPP/VIIRS進行降噪,以0.3為下限閾值,去除影像中的微弱燈光,同時將影像中的負值修改為0;以同年份北京、上海、廣州3 個城市中亮度均值與一倍標準差之和的最大值為上限閾值,剔除研究區內極高亮度的異常像元,一定程度上縮小NPP/VIIRS 影像的值域范圍,實現DMSP/OLS 高亮度飽和的效果[4]。NPP/VIIRS 與DMSP/OLS 數據量綱并不相同,為統一2 種類型的數據,利用兩者共有的2012年影像,利用公式(1)得到擬合方程與相關系數(表1)。參考Zhao[5]等研究,基于擬合方程生成2016、2019年類DMSP/OLS 影像,最終得到2000—2019年一致穩定的長時間序列夜光遙感影像。
2.2.1 建成區提取與精度驗證
建成區指行政區域內用于建設發展的非農業土地。本文利用動態閾值法提取建成區,根據統計年鑒中青島市各年份的建成區面積,動態地調整夜光影像中非建成區與建成區的閾值,當提取建成區面積與統計年鑒中數據的絕對差值最小時,即為最佳閾值。大于等于該閾值的像元是建成區,小于該閾值的像元被認為是非建成區。
建成區精度驗證包括數量精度驗證與空間精度驗證,數量精度在動態閾值選擇過程中已經確定,空間精度由Landsat TM/OLI 影像隨機選擇樣本點計算總體分類精度加以評判,具體結果見表2。

表2 2000-2019年青島市建成區提取精度評價
2.2.2 標準差橢圓
標準差橢圓是一種能夠直觀反映出點要素密集程度、方向性以及空間分布特征的空間統計方法[6]。標準差橢圓由3個要素確定:重心、方位角和X、Y軸長度。其中,重心是指城市建成區的重心位置,可反映城市建設重心的移動趨勢;方位角是橢圓長軸按順時針旋轉與正北方向的夾角,表征城市分布的主趨勢方向。X、Y 軸由所在軸的標準差決定,長軸和短軸的長度分別表征城市分布方向和分布范圍。長短軸與半軸的差距越大(扁率越大),表示城市布局的方向性越明顯。
2.2.3 燈光指數
綜合夜間燈光指數(CNLI)與中國城市化水平復合指標之間顯著相關,可以概括出地區的城市化水平[7]。CNLI 為某地區燈光面積比例(LAP)與平均燈光強度(MLI)的乘積。
式中,Arealight為燈光區域面積;Area 為地區的總面積;DNi為燈光像元的亮度值;Ci為像元值是DNi的像元數量;DNmax為像元亮度最大值。
歸一化植被指數(NDVI)是監測區域內植被狀態的有效指示因子。本文利用NDVI 數據計算建成區、非建成區以及全市NDVI 平均值,觀察其變化趨勢,反映青島市植被覆蓋的時空變化[8]。變化趨勢分析通過計算趨勢斜率slope,得出NDVI 隨時間的變化趨勢。正值表示植被覆蓋呈增加趨勢,負值表示減少趨勢(公式5)。
式中,n為年數,NDVIi、NDVI 為第i年NDVI和多年平均值。
城市熱島效應是指城區地表溫度高于周圍郊區的現象。本文通過構建熱島強度指數來描繪青島市熱島效應的時空變化[9]。借助ArcGIS 在建成區創建10 km緩沖區代表郊區,用建成區溫度與郊區平均溫度的差值作為熱島強度指數。
式中,Ii為影像中第i個像元的熱島強度指數;Ti為第i個建成區像元的溫度值;n為郊區有效像元數量;Tj為郊區第j個像元的溫度值。為反映城市熱島效應等級分布,根據熱島強度指數的范圍,本文將青島市劃分為強負熱島([-∞,-5.0))、弱負熱島([-5.0,-1.0))、無熱島([-1.0,1.0] )、弱熱島((1.0,5.0] )和強熱島([5.0,+∞))等5個等級。
利用Pearson 線性相關的方法,計算燈光指數與植被指數和熱島強度指數之間的相關性,反映城市擴張過程對生態環境的影響(公式7)。顯著性由t 檢驗決定,采用p<0.05 和p<0.01分別定義統計顯著和極顯著水平。
式中,r為相關系數;x和y分別為統計變量;xˉ和yˉ分別為統計變量平均值;n為樣本數。
表2 為青島市建成區提取數量精度和空間精度。各年份建成區提取的數量精度均超過91%,空間精度均大于82%(2016年稍低,為79%),結果滿足研究需求。
近20 a青島市建成區范圍逐步擴大,從2000年的199 km2增長到2019年的934 km2,面積增加約735 km2,年均增長38.7 km2。2008—2012年增長最慢,年均增長率僅19.0 km2;2012—2016年增長最快,年均增長達54.0 km2;2004—2008年和2016—2019年增長相對較快,年均增長接近50 km2。大體上,青島市建成區擴張速率經歷了“上升-下降-再上升-再下降”的波動式上升趨勢(表2)。空間上,東南部沿海區變化最為明顯,形成了以市南區、市北區、李滄區、城陽區和黃島區為主體的環膠州灣市轄區;部分內陸地區(如平度市、膠州市)等則以縣城為中心向外圍發散擴張。總體上看,2012年以來青島市建成區擴張明顯加快,在原城區周圍出現了大量零散城鎮,城市建設整體呈現出由沿海向內陸擴展,表現為“核心市區不斷擴大,周邊鄉鎮快速興起”的點面結合的城市發展特征。
從標準差橢圓分布看,橢圓的面積、方位角、重心等參數都發生了明顯變化(圖2)。標準差橢圓面積連年增大,說明城市建成區分布范圍持續向外擴張,以2000—2008年變化最大,隨后逐漸平穩(圖2a、2b)。橢圓方位角從2000年由5°增長至2019年的17°,總體略微向東北方向偏移(圖2d)。X、Y 軸之比2000—2004年驟然下降,隨后穩定在1.6 左右,表明城市前期擴張的方向性明顯,后期基本穩定不變(圖2c)。橢圓重心表現為由東南向西北移動的趨勢,反映了建成區逐漸向西、北方向擴展。重心移動過程以2000—2008年最為明顯,移動了7.5 km,2012年后移動變化趨于緩和(圖2a)。

圖2 2000-2019年青島市標準差橢圓及參數
近20 a青島市平均燈光指數呈平穩上升趨勢,指數由0.14 增長到0.26,年均增長率約0.006(圖3)。指數幅度和增長率較低,表明青島市城市建設進程持續向前,但整體程度不高。不同縣區燈光指數均有不同程度增長,但差異明顯。其中,市南區、市北區和李滄區始終最高(大于0.8),市北區接近1;萊西市和平度市始終最低(不足0.2)。

圖3 2000-2019年青島市不同區縣燈光指數
近20 a 青島市建成區平均NDVI 呈持續上升趨勢,年均增長率約0.0022。非建成區平均NDVI 經歷了先上升(2000—2008年)再下降(2008—2019年)的變化過程;全市平均NDVI 整體變化與非建成區一致(圖4a)。近年來青島市整體植被覆蓋表現出建成區內部增加,外部減少的空間態勢,反映了建成區在向外擴張過程中其他土地類型(如林地、耕地等)向建設用地轉化,降低了植被覆蓋程度;城市建成區植被覆蓋指數連年增加(圖4b),增加區域集中在建成區內,下降區域分布在市轄區邊緣及新興城鎮附近(圖4c),建成區植被覆蓋上升主要由老城區內部植被覆蓋增加所致,一方面說明城市綠化水平不斷提高,生態呈向好趨勢;另一方面說明城市邊緣及新興城鎮面臨著植被退化的問題。

圖4 2000-2019年青島市和城市建成區植被變化
從平均熱島強度看,近20 a青島市始終處在弱熱島區,但呈增加趨勢。2000—2012年增幅較大,2012—2019年增幅變緩,期間總增加約1.5℃,年增長率約0.0692℃(圖5a)。無熱島類型呈下降趨勢,弱熱島和強負熱島波動變化;強熱島類型表現為先增加后下降,2000—2012年間急劇增加,面積從1.7%增至27.5%,2012年后開始逐漸下降。從空間分布來看(圖5b),青島市以弱熱島類型為主,分布范圍最廣泛,強熱島、無熱島和弱負熱島所比例較少,缺少強負熱島。強熱島主要分布在市轄區,無熱島和弱負熱島區集中在靠近海岸或嶗山區域。綜合來看,隨著建成區生產、生活強度日益增加,市區與郊區的溫度差異擴大,熱島效應日益明顯。

圖5 2000-2019年青島市熱島強度指數和空間分布
全市和市轄區燈光指數與熱島強度指數呈顯著正相關,全市相關性最強(r=0.8888,p<0.05),建成區稍弱(r=0.8645,p<0.05)。縣級市燈光指數與熱島強度指數相關性極弱且不顯著(r=0.1897,p>0.05)(表3)。總體看,城市燈光指數與熱島效應具有明顯正相關性,相對于縣級市,城市建成區熱島效應更加明顯。

表3 2000-2019年青島市燈光指數與植被指數和熱島強度相關性
本文以快速發展的沿海城市—青島市為例,基于多源數據(遙感影像、統計數據等),分析了2000年以來青島城市擴張過程及時空演化特征,揭示了城市擴張對生態環境的影響。主要結論有:
1)利用燈光影像可實現城市擴張的連續監測。青島市建成區由沿海向內陸擴展不斷擴張。
2)全市平均NDVI呈先上升后下降趨勢,建成區NDVI 一直增長,升高區集中于市轄區,下降區多分布在城市邊緣和新興城鎮附近。熱島強度呈增強趨勢,但整體始終處于弱熱島態,強熱島分布在建成區內部,空間范圍逐漸增大。
3)不同區域的燈光指數與植被指數和熱島強度指數的相關性不同。城市燈光指數與熱島效應具有明顯正相關性,城市建成區熱島效應更加明顯;而城市燈光指數與植被指數沒有明顯相關性。隨著城市建成區的擴張,熱島效應逐漸增強。