典型喀斯特山地城市不透水面時空格局

母浩江， 崔文剛， 魯 霞

(貴州師范大學(xué)地理與環(huán)境科學(xué)學(xué)院， 貴州貴陽550025)

近年來，隨著我國城市化進程不斷加快，大量農(nóng)村人口向城市集聚，在促進城市經(jīng)濟發(fā)展的同時也導(dǎo)致城市用地壓力日益增大，越來越多的森林、 耕地等自然地表被不透水地表所替代。不透水地表又被稱為不透水面(impervious surface，IS)，是指常見的道路、 廣場、 建筑等水不能通過其下滲到土壤中的城市人工景觀[1-2]。IS的大量增加會改變地表下墊面性質(zhì)[3]，使地表徑流和地下水的活動受到干擾，致使城市內(nèi)澇、 地下暗河斷流等現(xiàn)象產(chǎn)生[4-5]，進而給生態(tài)環(huán)境帶來負面影響。及時、 準確地掌握城市IS的分布情況和時空演變對于城市發(fā)展規(guī)劃及生態(tài)環(huán)境保護具有重要意義。

早期IS信息的獲取主要依靠專業(yè)人員進行實地采集， 這樣的獲取方式成本高， 周期長， 時效性差， 使得IS信息的采集和應(yīng)用十分困難[6]。 進入20世紀以來， 隨著遙感技術(shù)的飛速發(fā)展， 遙感數(shù)據(jù)獲取變得便利， 國內(nèi)外大量學(xué)者開始將遙感技術(shù)引入IS的研究。 隨著研究的深入， 多種基于遙感技術(shù)定量獲取IS信息的方法相繼被提出， 較典型的方法可概括為光譜混合分析(spectral mixture analysis，SMA)法[7-8]、 決策樹模型[9-10]、 回歸模型[11]、 指數(shù)法[12-13]和分類法[14]5種， 其中線性光譜混合分析(linear spectral mixture analysis，LSMA)法可以有效解決中低分辨率遙感數(shù)據(jù)的混合像元問題， 被廣泛用于提取IS信息的研究[15]。 Wu等[7]使用LSMA法成功提取了美國大都市區(qū)哥倫布的IS信息； 唐菲等[8]基于陸地衛(wèi)星7(Landsat 7)增強型專題制圖儀(enhanced thematic mapper plus, ETM+)影像， 利用LSMA法提取了北京等6個不同城市的IS信息， 效果良好； 李苗等[16]通過LSMA法提取區(qū)域多個時間節(jié)點的IS信息， 分析IS信息的時空變化規(guī)律， 進一步探討城市化進程的范圍和規(guī)模。 理論方法的日趨成熟也推動著IS的研究向不同的區(qū)域環(huán)境與應(yīng)用領(lǐng)域擴展， 但是通過查閱文獻發(fā)現(xiàn)， 相關(guān)研究主要集中于經(jīng)濟較為發(fā)達的沿海城市[17]和平原地區(qū)[18]， 對于城市化進程起步稍晚且地形條件更復(fù)雜的喀斯特地區(qū)， IS的研究極少。

貴陽市是貴州省省會， 同時也是世界上喀斯特地貌發(fā)育最典型的地區(qū)之一[19]， 境內(nèi)喀斯特地貌占行政區(qū)域面積的73.3%， 地形起伏較大， 景觀破碎化嚴重， 生態(tài)環(huán)境脆弱且具有明顯的空間異質(zhì)性。 近年來， 隨著經(jīng)濟實力和城鎮(zhèn)化水平不斷提高， 貴陽已成為我國西部陸海貿(mào)易的新通道、 西南地區(qū)的中心城市之一。 在城市化狀況快速發(fā)展的同時， IS急速增加所造成的環(huán)境問題也引起人們的重視， 然而目前尚未發(fā)現(xiàn)關(guān)于該地區(qū)IS信息及其空間變化的研究。

本文中以貴陽市主城區(qū)為研究區(qū)，采用基于植被覆蓋度提取法和LSMA法分別提取研究區(qū)的IS信息，利用混淆矩陣對提取結(jié)果進行精度驗證，選出最優(yōu)方法對研究區(qū)1996—2017年各時間節(jié)點的IS數(shù)據(jù)進行提取，依據(jù)結(jié)果定量分析21 a間研究區(qū)IS的時空分布和演化規(guī)律。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

貴陽市位于我國的西南部(東經(jīng)106°07′—107°17′， 北緯26°11′—27°22′)，地處長江和珠江兩大水系流域的分水嶺地段，四面環(huán)山，地勢西高東低，境內(nèi)丘陵起伏，屬典型的喀斯特地形地貌，氣候類型為亞熱帶季風氣候[20]。本文中選取該市主城區(qū)為研究區(qū)，包括南明區(qū)、 云巖區(qū)、 觀山湖區(qū)、 白云區(qū)、 花溪區(qū)、 烏當區(qū)6個行政區(qū)，其中云巖區(qū)、 南明區(qū)大部分區(qū)域位于盆地，烏當區(qū)、 花溪區(qū)、 白云區(qū)和觀山湖區(qū)大部分區(qū)域以山地為主。

1.2 數(shù)據(jù)源與預(yù)處理

為了避免由影像數(shù)據(jù)源不同、 拍攝季節(jié)不同、 云霧遮蔽等原因造成的誤差， 保證研究數(shù)據(jù)的可靠性， 本文中選取研究區(qū)內(nèi)無云遮蔽且拍攝月份相近的3期Landsat影像作為數(shù)據(jù)源， 影像數(shù)據(jù)包含4景陸地衛(wèi)星5號(Landsat 5)專題成像儀(TM)影像和2景陸地衛(wèi)星8號(Landsat 8)陸地成像儀(OLI)影像(見表1)。 影像數(shù)據(jù)來自中國科學(xué)院地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn)。前期需要對影像數(shù)據(jù)進行預(yù)處理， 包括輻射校正、 大氣校正、 拼接以及按研究區(qū)范圍裁剪。

表1 研究區(qū)遙感數(shù)據(jù)

2 研究方法

2.1 IS信息的提取

2.1.1 基于植被覆蓋度提取法

基于植被覆蓋度提取IS數(shù)據(jù)的方法最早由Carlson等[21]提出。該方法根據(jù)IS與植被覆蓋度的負相關(guān)關(guān)系提取IS信息，計算方法相對簡單，并且能夠快速獲取大范圍連續(xù)的IS信息，已經(jīng)成為現(xiàn)階段應(yīng)用最多的IS信息提取方法。該方法的表達式[13-14，21]為

I=1-F，

(1)

(2)

(3)

式中：I為IS指數(shù)，%；F為植被覆蓋度，%；N為歸一化植被指數(shù)；Nsoil為N在裸土或無植被區(qū)的值；Nveg為N在純植被區(qū)的值；NIR為近紅外波段的反射率，%；R為紅光波段反射率，%。由于影像存在不可避免的噪聲，因此在一般情況下，Nveg和Nsoil取一定置信度范圍內(nèi)的最大值與最小值。在參考前人經(jīng)驗[22]并多次試驗后，本文中選取累積分數(shù)95%～5%作為置信范圍，得到N的最大值和最小值。

在利用植被覆蓋度提取IS信息前，還需對研究區(qū)內(nèi)的水體信息進行掩膜。徐涵秋[23]提出的改進型歸一化差異水體指數(shù)(modified normalized difference water index，MNDWI)能夠清楚地分出陰影和水體，是現(xiàn)階段提取效果最好的水體提取法，因此本文中使用該指數(shù)對水體信息進行提取、 掩膜。該指數(shù)的表達式為

(4)

式中：M為歸一化水體指數(shù)；G為綠光波段的反射率，%。

2.1.2 LSMA法

LSMA法通過計算混合像元中各端元的組合比例，將混合像元的光譜分解為多種端元光譜的線性組合[8]。該方法假設(shè)一個像元的反射率由影像拍攝瞬間各地物的純凈反射率(端元)按一定比例組合而成，即

(5)

式中：Rx為波段x的光譜反射率，%；n為端元個數(shù)；fi為端元i在像元中所占的比例，%；Ri,x為端元i在波段x的反射率，%；ex為殘差。

理論上， 端元的數(shù)量在一定程度上決定了LSMA法的擬合程度， 但這并不意味著端元種類越多，模型效果越優(yōu)。Wu等[7]通過實際研究發(fā)現(xiàn)，選取4種端元時IS信息的反演效果最優(yōu)。本文中依據(jù)前人經(jīng)驗[7]對影像進行端元采集時選擇由玻璃、混凝土為代表的高反射率地物，傳統(tǒng)的舊式瓦片屋頂、瀝青路面為代表的低反射率地物，以及植被和土壤4種端元類型。除此之外，研究[2]發(fā)現(xiàn)，水體在LSMA法中的表現(xiàn)為低反射率。為了避免水體對研究結(jié)果的影響，本文中在提取純凈端元之前使用式(4)中的修正歸一化水體指數(shù)，對研究區(qū)內(nèi)的水體進行剔除，將最終得到的高反射率蓋度圖和低反射率蓋度圖相加得到研究區(qū)IS蓋度圖[24]。

2.2 精度驗證

對影像分類提取結(jié)果驗證的最好方式是隨機提取采樣點進行實地驗證；但是，實地驗證成本較高，對過往年份展開實地驗證也不切實際，因此，本文中先利用地理信息系統(tǒng)(GIS)工具隨機生成驗證點，參照驗證點的空間位置對基準影像(2007、 2017年基準影像采用同期分辨率為1.07 m的Google影像，1996年的Google影像因效果不佳而采用原始影像)進行目視解譯，將解譯結(jié)果與二值化后的IS結(jié)果進行對比，對比結(jié)果通過混淆矩陣來反映提取精度。

2.3 IS的空間分布與演變

2.3.1 標準差橢圓分析法

為了探究研究區(qū)IS的空間變化情況，嘗試利用標準差橢圓分析法進行分析。通過分析標準差橢圓，可以從橢圓中心、 分布范圍、 方向和形狀等多個角度分析要素的空間分布特征及空間演變過程[25]。在地學(xué)領(lǐng)域，主要用橢圓的長軸半徑長度與扁率反映空間格局總體要素的集中密度，用長軸偏角反映分布格局的主導(dǎo)方向。

2.3.2 IS的擴展指標

IS擴展指標能夠全面反映研究區(qū)IS擴展的情況。 本文中選擇IS密度、 IS擴展速度這2個指標來反映各區(qū)域IS擴展趨勢[26]。 IS密度是指某區(qū)域中IS的面積占該區(qū)域總面積的比例； IS擴展速度是指在一段特定時間段內(nèi)， 以年為單位的IS面積增長情況。

IS密度計算公式為

(6)

式中：ρa為區(qū)域a的IS密度，%；Aa為區(qū)域a的IS的面積， km2；Aa,tot為區(qū)域a的總面積， km2。

IS擴展速度計算公式為

(7)

式中：v為區(qū)域a的IS擴展速度， km2/a；t為時間間隔，a；Apri和Aend分別為區(qū)域a研究初期和末期的IS面積， km2。

3 結(jié)果與分析

3.1 精度驗證結(jié)果分析

將2種IS信息提取方法獲取的研究區(qū)2017年的IS影像選擇合適閾值進行二值化處理。利用生成的驗證點，分別將2幅二值化IS影像與基準影像疊在一起進行驗證，結(jié)果見圖1。基于植被覆蓋度提取法提取的IS總體分類精度為0.89，Kappa系數(shù)為0.53，而LSMA法得到的總體分類精度為0.95，Kappa系數(shù)為0.74。此外，在目視判斷過程中發(fā)現(xiàn)基于植被覆蓋度提取IS信息的方法很容易將大片裸地誤分為IS，而用LSMA法時該現(xiàn)象很少出現(xiàn)。由此可知，采用LSMA法在喀斯特地區(qū)具有較好的可行性，并且提取效果優(yōu)于基于植被覆蓋度提取法。

在后續(xù)研究中選擇LSMA法繼續(xù)提取研究區(qū)1996、 2007年的IS信息， 并進行精度驗證， 結(jié)果如表2所示。 由表可知， 本文中提取的IS信息均能滿足研究要求， 可用于分析研究區(qū)IS空間格局演變。

表2 不同方法對貴陽市主城區(qū)不透水面的提取精度驗證

3.2 IS時空分布

利用GIS技術(shù)對得到的1996、 2007、 2017年的IS二值化圖(見圖2)進行面積統(tǒng)計， 由圖2可知， 研究區(qū)1996、 2007、 2017年的IS面積分別為50.82、 126.56、 234.08 km2， 21 a間貴陽市IS面積共增加183.26 km2，增長了近5倍，整體呈現(xiàn)出高速擴張趨勢，其中1996—2007年IS面積增加了75.76 km2，年平均增長量為6.89 km2/a；2008—2017年IS面積增加量達107.49 km2，年平均增長量為10.75 km2/a，前10 a的平均增長量明顯高于后10 a的。在2008—2017年中， 貴陽市內(nèi)、 外部交通快速發(fā)展， 內(nèi)部以環(huán)城高速和中環(huán)為代表的各等級道路竣工通車， 外部貴遵附線、 花安高速以及貴廣、 滬昆高速鐵路建成通車， 大力推動了貴陽市城市的發(fā)展。 除此之外， 觀山湖區(qū)作為緩解舊城區(qū)人口稠密、 交通擁堵的發(fā)展新區(qū)， 是城市擴展的重點方向， 在大量文化商貿(mào)區(qū)、 居住區(qū)、 公共基礎(chǔ)設(shè)施建設(shè)的同時， 也導(dǎo)致了IS的增長。

將提取的IS信息按行政區(qū)進行分區(qū)統(tǒng)計， 結(jié)果見表3。 由表可看出， 1996—2007年IS面積增長量最大的是花溪區(qū)， 增長了25.72 km2，增長量最小的是云巖區(qū)，僅為8.59 km2。通過觀察各區(qū)IS面積增長率發(fā)現(xiàn)，觀山湖區(qū)的增長率最高，達到了610.67%， 其次是花溪區(qū)、 烏當區(qū)， 增長率分別為377.13%、 369.85%，而南明區(qū)的增長率最低，僅為45.91%。2008—2017年IS面積增長最多的是觀山湖區(qū)，其次是花溪區(qū)，增長最少的是南明區(qū)。對比1996—2007年的數(shù)據(jù)，各區(qū)IS面積增長率明顯降低且區(qū)域間降低幅度差異較大，其中觀山湖區(qū)增長率下降超400%，但其增長率仍居各行政區(qū)值之首，南明區(qū)21 a間IS面積增長率都最低，但其一直處于平穩(wěn)增長狀態(tài)。

表3 1996—2007年貴陽市主城區(qū)各區(qū)域不透水面面積的變化情況

從貴陽市各行政區(qū)IS面積的總體變化趨勢可以看出： 1996—2007年觀山湖區(qū)、 花溪區(qū)、 烏當區(qū)處于城市化擴張爆發(fā)期， 大量的IS井噴式出現(xiàn)； 南明區(qū)處于城市穩(wěn)定發(fā)展期， IS格局基本穩(wěn)定， 增長相對緩慢； 云巖區(qū)、 白云區(qū)則屬于典型的中速發(fā)展期， IS平穩(wěn)擴展。 2008—2017年觀山湖區(qū)、 烏當區(qū)仍處于高速發(fā)展階段， 強度明顯減弱； 花溪區(qū)與云巖區(qū)、 白云區(qū)一樣進入城市中速發(fā)展期， IS繼續(xù)平穩(wěn)擴張。

為了更直觀地展現(xiàn)各階段研究區(qū)IS面積變化的空間分布情況， 利用GIS工具計算研究區(qū)不同階段之間IS面積差值， 由此得到不同階段IS面積變化的圖斑， 結(jié)果見圖3。 由圖可以看出， 1996年以前研究區(qū)IS分布較集中， 主要分布在以南明區(qū)、 云巖區(qū)為主的市中心； 1996—2007年IS面積變化主要分布在各行政區(qū)中心周圍， 同時城市中心的IS面積增加也很顯著； 2008—2017年， IS的增加使得之前許多破碎的斑塊連接起來， IS開始出現(xiàn)成片聚集現(xiàn)象， 還有很多新增IS以“飛地”形式出現(xiàn)在各行政區(qū)中。

3.3 IS格局演變分析

3.3.1 IS空間變化特征

通過分析研究區(qū)1996、 2007、 2017年IS的標準差橢圓，探討該區(qū)域IS的整體變化趨勢， 如圖4所示。由圖可以看出：1996—2007年，標準差橢圓面積明顯增大且整體向南移動，表明此階段IS擴展重心正在向南傾斜；主軸方向偏角從9.87°減小到2.30°， 意味著該區(qū)域IS的分布趨勢由東北—西南方向往正北—正南方向發(fā)展； 標準差橢圓的長半軸由19.11 km增大到22.58 km， 短半軸由10.46 km增大到13.51 km，表明IS的空間分布在主軸方向(南北方向)和輔軸方向(東西方向)都呈現(xiàn)較弱的離散化現(xiàn)象。

圖3 貴陽市主城區(qū)不同時期不透水面的變化

2008—2017年研究區(qū)IS標準差橢圓范圍逐漸縮小且向東北方向移動,擴展重心逐漸向西北方向移動； 主軸方向的偏角從2007年的2.30°增至17.02°， 該階段貴陽市主城區(qū)IS的分布方向轉(zhuǎn)變?yōu)闁|北—西南方向， 與1996年相同且趨勢更加顯著； 長半軸由22.58 km減至22.2 km， 短半軸由13.51 km增至14.67 km， 可知IS的空間分布在主軸方向(南北方向)的離散化趨勢有所緩解，在輔軸方向(東西方向)離散化程度進一步加強。

(a)1996年 (b)1996—2007年 (c)1996—2017年

縱觀研究區(qū)IS的變化情況， 與研究區(qū)地形密切相關(guān)。 貴陽市屬于典型的喀斯特山地城市， 地表景觀破碎化嚴重。 南明區(qū)、 云巖區(qū)是貴陽市的老城區(qū)， 大部分區(qū)域位于貴陽—中曹司盆地內(nèi)， 四面環(huán)山， 受地形嚴重限制， 迫使其城市范圍不斷尋求向外擴展， 進而致使IS向其他行政區(qū)擴張， 但內(nèi)部填充式擴張仍是其擴展的重要方式， 這種狀況一直保持到2000年。 2000年后， 大量的交通干線開始建設(shè)， 以埡口和河谷為走廊逐步向地勢平坦的花溪區(qū)、 觀山湖區(qū)延伸[27]。 2001—2006年間，花溪區(qū)成為貴陽市交通的主要擴展地區(qū)， 交通建設(shè)量占同期工程建設(shè)總量的45%。 從2006 年開始， 貴陽市城市建設(shè)發(fā)展的主方向轉(zhuǎn)移至觀山湖區(qū)， 其次是花溪區(qū)， 這2個方向擴展的比例達到同期工程建設(shè)總量的80%[27]， 這些變化情況都與研究區(qū)IS標準差橢圓的變化相呼應(yīng)。

3.3.2 IS擴展趨勢

1996—2017年貴陽市主城區(qū)各區(qū)域IS的擴展指標如圖5所示。 從圖5(a)中可以看出， 1996、 2007、 2017年研究區(qū)的IS密度分別為2.14%、 5.32%、 9.84%， 22 a間增加了7.7%。各行政區(qū)的IS密度變化存在明顯差異，其中云巖區(qū)的IS密度在2個研究階段都是最高的，1996、 2007、 2017年的IS密度分別為12.23%、 21.62%、 36.33%，而烏當區(qū)在整個研究階段IS密度最低， 增長最緩慢。從圖5(b)中可以看出：1996—2007年研究區(qū)IS的擴展速度為6.89 km2/a，各行政區(qū)IS擴展極不均衡，其中花溪區(qū)的擴展最快，擴展速度為2.34 km2/a，是各行政區(qū)擴展速度平均值的2倍，烏當區(qū)的IS擴展最慢，擴展速度僅為0.67 km2/a；2008—2017年，研究區(qū)IS擴展速度整體上升至10.75 km2/a，花溪區(qū)雖然保持較高的擴展速度，但是擴展速度呈倒退趨勢，其余各區(qū)擴展速度明顯提升，特別是觀山湖區(qū)擴展速度上升最多，達到了2.65 km2/a，擴展速度最低的南明區(qū)也較前一階段提升了0.4 km2/a。

(a)IS密度

4 結(jié)論

1)通過對LSMA法和基于植被覆蓋度提取法得到的貴陽市1996—2017年IS信息進行精度驗證可知，LSMA法的提取結(jié)果明顯優(yōu)于基于植被覆蓋度提取法的，說明在喀斯特地貌環(huán)境下使用LSMA法提取IS信息是適用的，并且效果更優(yōu)。

2)1996—2017年貴陽市IS面積從50.82 km2增加到210.12 km2，該市的IS擴展可分為3個階段：前期(1996年)IS分布比較集中，主要分布在城市中心區(qū)域；中期(1996—2007年)主要呈現(xiàn)出以各行政區(qū)為發(fā)散中心向四周擴張，同時伴隨著城市中心區(qū)域的內(nèi)部填充；后期(2008—2017年)IS的擴充得到進一步發(fā)展，擴展方式主要呈現(xiàn)出多中心發(fā)散的特點，并出現(xiàn)IS“飛地”，其中觀山湖區(qū)擴張尤為突出。

3)研究區(qū)IS面積的標準差橢圓偏角從1996年的9.87°增加到2017年的17.02°，表明研究區(qū)IS的分布趨勢表現(xiàn)為東北—西南方向，并且隨著城市的發(fā)展，該分布態(tài)勢得到進一步加強。

4)本文中將LSMA法應(yīng)用于喀斯特地區(qū)IS信息的提取研究，與基于植被覆蓋度提取法提取的IS信息相比，LSMA法更適用于喀斯特地區(qū)。本文中還提取了研究區(qū)多個時間節(jié)點的IS信息，對IS信息的時空演變進行了進一步探討，研究結(jié)果為喀斯特地區(qū)IS研究和貴陽市后期城市發(fā)展規(guī)劃、 生態(tài)文明城市建設(shè)提供了參考。對提取到的IS信息開展相關(guān)應(yīng)用研究將是后期研究的重點。