欽州市城市擴張的生態響應*

劉文良，李 偉，郝連成

(中國地質調查局煙臺海岸帶地質調查中心，山東煙臺 264004)

0 引言

城市化是指人口向城市聚集，并由此引起一系列社會變化的過程[1]。城市向郊區擴張是這一過程最為直觀的特征。隨著城市的擴張出現了一些無法避免的生態環境問題，例如城市熱島、水污染等。城市建設和生態保護之間的矛盾一直是科學研究的重要內容。客觀、準確地評價城市生態環境的變化并發現其存在的問題，為城市的規劃和建設提供建議，有利于區域的持續發展。

當前遙感技術以大面積、快速、周期性觀測的優點，在生態環境監測領域[2,3]得到了廣泛的應用。但目前以單一指標評價居多，如使用植被指數監測地表覆被及其對氣候的響應[4]，利用地表溫度估測城市熱島效應[5]，通過反演不透水面分析城市擴張變化[6]，等等。這種單一指標的評價相對片面，不能完整地反映城市生態的變化。綜合多個指標的生態環境評價方法有層次分析法(AHP)[7]、模糊數學評判法[8]、壓力-狀態-響應(PSR)模型法[9]等，這些方法存在人為確定權重的問題，且要使用大量的統計數據，模型構建相對復雜。徐涵秋[10]提出的遙感生態指數(Remote Sensing Ecological Index，RSEI)基于遙感技術，用主成分分析的方法集成了多個指標，能夠簡單、快速、客觀、定量地評價區域整體的生態環境狀況，還可以進行長時間序列的可視化空間表達，目前已經在許多城市[11,12]和地區[13,14]得到了應用。

欽州市是我國大西南的一個典型港城分離型城市。自20世紀90年代以來，欽州市城區和欽州港港區規模有了大幅度的增長，且有進一步擴大的趨勢，快速的城市擴張是否影響了該市的生態環境是一個值得關注的問題。對欽州市進行整體的生態評價，研究制約生態環境變化的驅動因子是非常有必要的。本文基于RSEI，對欽州市2000-2017年的生態環境變化進行量化分析評價，發掘影響該市生態環境狀況的關鍵因子，為城市建設和生態保護工作提供理論依據。

1 材料與方法

1.1 研究區概況

欽州市位于廣西南部沿海，在北部灣經濟區居中心位置，是面向東盟開放的重要窗口之一。該地區屬海洋季風氣候，夏季高溫多雨，冬季溫和少雨。境內地貌以低山丘陵為主，地勢自北向南傾斜。1996年，經廣西壯族自治區人民政府批準，欽州市成立欽州港經濟技術開發區，2010年升級為國家級開發區。2019年8月，中國(廣西)自由貿易試驗區欽州港片區獲批。自欽州港開發區成立以來，港口的建設和發展推動了欽州市城市化的進程，現已形成以老城區為核心的生活中心和以欽州港為核心的工業中心。本研究選取老城區和欽州港所在的欽南區作為研究區，該區域是近20年來欽州城市發展最快的區域，能集中反映城市生態變化的狀況，其地理跨度為21°36′—22°08′N，108°24′—109°09′E，總面積約為2 383 km2。

1.2 數據來源

本研究選取Landsat系列影像作為數據源，包括2000年11月6日的ETM+影像、2009年10月6日的TM影像和2017年10月28日的OLI/TIRS影像。影像獲取季節一致，避免季相差異所帶來的影響。為保證各指數計算的準確性，首先分別對3期影像進行輻射定標、大氣校正、幾何校正等預處理，再在此基礎上進行研究區的裁剪和各指標分量的計算。

1.3 方法

遙感生態指數(RSEI)基于遙感技術，利用主成分分析的方法集成綠度、濕度、干度、熱度4個指標，以綜合客觀反映區域的生態環境狀況。

1.3.1 指標提取

1.3.1.1 綠度指標

綠度指標用歸一化植被指數(Normalized Difference Vegetation Index,NDVI)來表示，該指數是最常用的反映植被覆蓋和生長狀況的參量，其計算公式為

NDVI=(ρNIR-ρR)/(ρNIR+ρR)，

(1)

式(1)中，ρR、ρNIR分別為影像紅波段和近紅波段的反射率，分別對應TM/ETM+影像的第3，4波段，OLI影像的第4，5波段。

1.3.1.2 濕度指標

纓帽變換技術獲取的濕度分量能夠反映植被和土壤的濕度狀況，廣泛應用于生態監測[15]，故用來表示濕度指標，其計算公式為

Wet=c1ρB+c2ρG+c3ρR+c4ρNIR+c5ρSWIR1+c6ρSWIR2，

(2)

式(2)中，Wet表示濕度分量，ρB、ρG、ρR、ρNIR、ρSWIR1、ρSWIR2在TM/ETM+影像中分別對應第1，2，3，4，5，7波段的反射率，在OLI影像中分別對應第2，3，4，5，6，7波段的地表反射率。c1、c2、c3、c4、c5、c6為纓帽變換中各波段的系數，不同傳感器的系數不同；對于TM影像，Crist[16]推導的系數依次為0.031 5，0.202 1，0.310 2，0.159 4，-0.680 6 和0.610 9；對于ETM+影像，Huang等[17]推導的系數為0.262 6，0.2141，0.092 6，0.065 6，-0.762 9，-0.538 8；對于OLI影像，李博倫等[18]推導的系數依次為0.265 1，0.236 7，0.129 6，0.059 0，-0.750 6，-0.538 6。

1.3.1.3 干度指標

自然環境中，造成地表干化的地物主要是建筑物和裸土，歸一化建筑-土壤指數(Normalized Difference Building-Soil Index,NDBSI)能表征建筑物和裸土的信息，由IBI和SI組合而來，本文用兩者的均值來表示干度指標，其計算公式為

NDBSI=(IBI+SI)/2，

(3)

IBI=[2ρSWIR1/(ρSWIR1+ρNIR)-ρNIR/(ρNIR+

ρR)-ρG/(ρG+ρSWIR1)]/[2ρSWIR1/(ρSWIR1+ρNIR)+ρNIR/(ρNIR+ρR)+ρG/(ρG+ρSWIR1)]，

(4)

(5)

1.3.1.4 熱度指標

地表溫度(Land Surface Temperature,LST)在城市熱環境研究、生態環境監測以及災害預警等領域具有廣泛的應用，用來表征熱度指標。Landsat系列衛星一直是地表溫度反演的重要數據來源之一，其中，Landsat 5和Landsat 7各有一個熱紅外波段(TM 6、ETM+ 6)，Landsat 8 TIRS傳感器有兩個熱紅外波段(TIRS 10，11)，由于TIRS 11波段的定標參數存在誤差，美國地質調查局(United States Geological Survey，USGS)不鼓勵使用該波段進行地表溫度反演[19]。綜合考慮，本文使用普適性強，需求參數少的單通道算法對地表溫度進行反演。單通道算法由Jiménez-Munoz等[20]于2003年提出，經過多次改進[21,22]，適用于Landsat系列衛星傳感器，并且可以取得1.5 K的精度[23,24]。該算法主要公式為

LST=γ[ε-1(φ1L+φ2)+φ3]+δ，

(6)

γ=T2/(bγL)，

(7)

δ=T-T2/bγ，

(8)

公式(6)—(8)中，LST為地表溫度，單位為開爾文(K)；對于TM 6、ETM+ 6、TIRS 10波段，參數bγ分別為1 256，1 277和1 324；L為TM 6、ETM+ 6、TIRS 10波段輻射亮度值，可根據影像頭文件中的gain和offset參數定標得到。φ1、φ2、φ3用下列公式計算：

φ1=1/τ，

(9)

φ2=-L↓-L↑，

(10)

φ3=L↓，

(11)

式(9)—(11)中，τ、L↓、L↑分別為大氣透過率、大氣下行輻射和大氣上行輻射，這些參數可以由網站http://atmorr.gsfc.nasa.gov查詢而來；T為熱紅外波段在傳感器處的亮度溫度，其公式為

T=K2/ln(K1/L+1)，

(12)

式(12)中，K2、K1為常數，對于TM 6、ETM+ 6、TIRS 10波段，K1的數值分別為607.76,666.09,774.89，單位為W·m-2·sr-1·μm-1，K2數值分別為1 260.56,1 282.71,1 321.08，單位為K；ε為地表比輻射率，本文使用覃志豪等[25]提出的方法進行反演。

1.3.2 遙感生態指數構建

1.3.2.1 指標分量標準化處理

上述公式計算的4個指標量綱是不一致的，因此，在計算遙感生態指數之前需要進行標準化處理，將各指標的量綱都統一到[0,1]之間。其公式為

NIi=(Ii-Imax)/(Imax-Imin)，

(13)

式(13)中，NIi為經過標準化處理之后的像元值，Ii為指標分量在i像元處的值，Imax、Imin為該指標分量的最大和最小像元值。

1.3.2.2 指標分量的水體掩膜處理

遙感生態指數中濕度分量(Wet)研究的是與植被和土壤濕度有關的參量，因此在計算時需要將大面積的水體提取出來并掩膜掉,以避免水域對主成分分析的載荷造成影響。本文對水體的提取采用徐涵秋[26]提出的歸一化差異水體指數(Modified Normalized Difference Water Index，MNDWI)，其公式為

MNDWI=(ρG-ρSWIR1)/(ρG+ρSWIR1)。

(14)

1.3.2.3 計算遙感生態指數

遙感生態指數(RSEI)是可以綜合以上4個指標信息的參量。主成分分析方法(Principal Component Analysis,PCA)可以去除波段之間多余信息，將多波段的圖像信息壓縮到少數幾個波段，其中第一主成分(PC1)包含了絕大多數的信息。該方法的最大優點是根據各指標對主分量的貢獻值客觀地確定各指標的權重，避免人為誤差，將經過標準化和水體掩膜處理的指標分量重新合成一幅新的影像。使用ENVI的主成分分析模塊進行計算，得出PC1和相關統計數據。然后依次將PC1進行正負值轉換、標準化處理，得到遙感生態指數(RSEI)，主要公式為

RSEI0=1-PC1，

(15)

RSEI=(RSEI0-RSEI0_min)/(RSEI0_max-RSEI0_min)，

(16)

式(15)中，RSEI0為初始生態指數，PC1為ENVI計算得到的第一主成分PC1。式(16)中，RSEI數值介于[0，1]之間，值越大，生態越好，反之越差；RSEI0_min、RSEI0_max為初始生態指數RSEI0的最小值和最大值。

1.3.3 建模分析

分別在研究區3個時期的NDVI、Wet、NDBSI、LST和RSEI影像隨機選取了3 000個樣點，基于樣點數據進行逐步多元回歸分析，并繪制散點圖，分析4個變量因子及其與生態環境之間的關系。

2 結果與分析

2.1 生態變化整體分析

表1統計了研究區4個指標分量和RSEI的均值,以及各指標對第一主成分(PC1)的載荷。由表1可以看出，NDVI和Wet對PC1的載荷均為正值，說明兩者對生態起正面作用;NDBSI和LST對PC1的載荷為負值，兩者對生態起負面作用，與客觀實際一致。2000年、2009年、2017年各分量對PC1的貢獻率分別為73.94%、74.59%、87.66%，表明PC1已經集成了4個分量的絕大多數信息。統計表明，2000年、2009年和2017年研究區的RSEI均值分別為0.676，0.698和0.725，整體增長了7.2%，生態環境有變好的趨勢。分析表1中4個指標對PC1的載荷可知， NDVI和Wet對PC1載荷的絕對值之和遠大于NDBSI和LST的絕對值之和，說明在該地區植被覆蓋和土壤濕度對生態環境的優化作用大于地面建筑和地表溫度對生態環境的破壞作用。4個指標中，又以NDVI對PC1載荷的絕對值最大，說明植被的覆蓋是影響該地區生態的主要因素。2000年、2009年、2017年NDBSI對PC1載荷的絕對值呈先增長后穩定不變的趨勢，表明2009年以前研究區地表的“干化”對生態環境的負面作用在增加，2009年以后趨于穩定。

表1 研究區4個指標和RSEI統計表

為定量化分析生態狀況，將RSEI以0.2為間隔分成差(0—0.2)、較差(0.2—0.4)、中(0.4—0.6)、良(0.6—0.8)和優(0.8—1.0)5個等級(界限值屬于前一級)。圖1顯示了研究區2000-2017年地表景觀的變化和RSEI等級的空間變化情況。整體來看欽州市植被覆蓋率高，生態等級以良為主。2000年為欽州港建設的初期，港口處出現小規模的填海，城區規模較小，生態質量差和較差的區域集中在港口的北部，其他區域有零星分布，主要地物類型為建筑用地、裸地等，城區的生態等級以中為主；2000-2009年，城區和港口的范圍都有了一定規模的擴張，生態等級差和較差的區域大都集中在這兩個區域，城區的生態等級轉變為以較差為主；2009-2017年，城區和港口的范圍雖然出現較大規模的擴張，但是部分區域建設初期破壞的自然表面得到了修復，如欽州港西部和老城區的一些區域，生態等級差和較差的區域得到稀釋，港口的生態等級以較差為主，城區生態等級以中為主，說明生態保護的理念逐步深入人心。

表2統計了3個時期RSEI各等級的面積及其所占比例，3個時期的總面積各不相同且比欽州市總面積稍小，主要是因為計算RSEI時對水體進行了剔除。結果表明：1)2000年、2009年、2017年研究區的生態質量等級以優、良為主，二者之和均占總面積的70%以上，且二者所占的比例之和有增長的趨勢，由2000年的74.97%增長到2017年的84.08%，二者面積之和增長247.1 km2,同比增長14.07%；2)生態質量等級較差和差的區域面積和比例先增加后降低。2000年，二者所占面積之和為45.33 km2，約占研究區當期總面積的1.94%；2009年，面積之和為59.38 km2，約占研究區當期總面積的2.52%；到2017年，面積為53.82 km2，約占當期總面積的2.26%。

2.2 生態環境變化監測

為分析研究區域18年來生態質量的差異變化，在RSEI分級的基礎上，對研究區2000年和2017年的RSEI指數做差值運算，并將結果劃分為變差、不變、變好三大類，結果如圖2所示。經統計，欽州市的總面積由2000年的2 342.19 km2增長至2017年的2 382.31 km2，面積增長了約1.7%，主要是圍填海(圖2紫色圖斑)的原因，且95%以上的圍填海區域都集中在南側的欽州港。由圖2可知，生態環境惡化的區域主要集中在城區和欽州港東北部，城區的擴張和欽州港工業園的建設是生態惡化的主要原因。生態環境有所改善的區域分布較廣，其中以城區東北部為典型，這些區域的地貌以低山丘陵為主，地物類型主要是林地和園地。

經統計(不含圍填海區)，大部分區域生態質量無明顯變化(圖2淡黃色圖斑)，占總面積的58.21%；生態質量變差的區域(圖2橙色圖斑)面積為282.75 km2，占總面積的12.72%；生態環境有所改善的區域(圖2淺綠色圖斑)面積為646.44 km2，占占面積的29.07%。

圖1 研究區2000年、2009年、2017年遙感影像和RSEI等級分布圖

表2 研究區生態環境評價等級統計表

圖2 研究區2000-2017年RSEI等級變化空間分布

2.3 城市建成區擴張對生態的影響

分別對2000年、2009年、2017年的欽州市城區和欽州港建成區進行解譯，其建筑面積由2000年的17.07 km2增長到2017年的107.27 km2，增加了5倍多。為研究城市的擴張對生態質量等級重心的影響，用ArcGIS分別將3期的RSEI柵格影像轉換為矢量要素點。由前文可知，建筑用地的生態質量等級大都為差和較差，因此分別提取城區和欽州港生態等級差和較差的要素點，并利用ArcGIS分別計算這些要素點的加權平均中心(紫色點為城區的，紅色點為欽州港的)，如圖3所示。圖3反映了城市建筑面積擴張和低生態等級區域中心的遷移關系。平均中心遷移的方向和距離反映了低生態等級區域的分布狀況，平均中心基本上隨建筑用地擴張的方向遷移，以上情況進一步表明，城市的擴張對生態有一定的負面影響。從遷移距離來看，欽州港平均中心大于城區，說明欽州港區域的建設力度比城區大。

2.4 建模分析

建模結果如下：

RSEI=0.5167+0.5069NDVI+0.8191Wet-0.1019NDBSI-0.0033LST，R2=1 (2000年);

圖3 城市擴張及低生態等級平均中心遷移圖

RSEI=0.9410+0.5079NDVI+0.2583Wet-0.1339NDBSI-0.0213LST，R2=1 (2009年);

RSEI=0.5413+0.5883NDVI+0.0438Wet-0.1277NDBSI-0.0080LST，R2=1 (2017年)。

從建模情況來看，3個方程的R2均為1，模型擬合度高，建模過程中，在4個變量都保留的情況下，標準誤差RMSE的值最低，說明這4個因子都是必要的。其中，NDVI和Wet的系數為正，與RSEI成正相關關系；NDBSI和LST系數為負，與RSEI成負相關關系，與主成分分析的結果相一致。

散點圖能直觀地反映自變量與因變量的相關性，以2017年為例(圖4)，NDVI、NDBSI與RSEI之間的線性關系較強，Wet、LST對RSEI的線性趨勢較弱。從2017年的建模方程和圖4a來看，NDVI是影響研究區生態狀況的關鍵因子，建模系數的絕對值遠大于其他因子，散點圖趨勢線的斜率達到0.68，這表明當地植被一旦遭到規模性破壞，其生態環境將迅速惡化。

圖4 2017年4個因子與RSEI的相關性散點圖

3 結論

本研究基于RSEI，對欽州市2000-2017年的生態環境狀況進行定量分析評價，從時間和空間兩個角度直觀地表達了近20年來欽州城市擴張的生態效應。結果表明：1)2000—2017年間，欽州市RSEI均值均在0.7左右，生態等級以良(RSEI 0.6-0.8)為主，且2017年的遙感生態指數均值最大，表明欽州市的生態環境質量在改善，城市的擴張給欽州市整體生態狀況帶來的負面影響不大。從這方面來看，欽州市仍有很大的發展潛力。2)從主成分分析和建模方程的系數來看，綠度對欽州市的RSEI貢獻程度最大，是影響生態的關鍵因子。3)從城市建成區擴張對低生態等級平均中心的影響結果來看，平均中心的遷移方向隨城市擴張方向移動，從遷移距離來看，欽州港的建設力度大于城區。4)生態等級的分布和變化監測結果表明，部分城市建設區域的生態等級由差變好，建設初期破壞的自然生態表面后期得到了修復，表明恢復植被覆蓋的園林綠化工程能有效地改善生態環境狀況。