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 km²。

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)来表示，该指数是最常用的反映植被覆盖和生长状况的参量，其计算公式为

$ \mathrm{NDVI}=\left(\rho_{\mathrm{NIR}}-\rho_{\mathrm{R}}\right) /\left(\rho_{\mathrm{NIR}}+\rho_{\mathrm{R}}\right), $

(1)

式(1)中，ρ_R、ρ_NIR分别为影像红波段和近红波段的反射率，分别对应TM/ETM+影像的第3，4波段，OLI影像的第4，5波段。

1.3.1.2 湿度指标

缨帽变换技术获取的湿度分量能够反映植被和土壤的湿度状况，广泛应用于生态监测^[15]，故用来表示湿度指标，其计算公式为

$ \mathrm{Wet}=c 1 \rho_{\mathrm{B}}+c 2 \rho_{\mathrm{G}}+c 3 \rho_{\mathrm{R}}+c 4 \rho_{\mathrm{NIR}}+c 5 \rho_{\mathrm{SWIR} 1}+ \\ c 6 \rho_{\text {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组合而来，本文用两者的均值来表示干度指标，其计算公式为

$ \mathrm{NDBSI}=(\mathrm{IBI}+\mathrm{SI}) / 2, $

(3)

$ \mathrm{IBI}=\left[2 \rho_{\mathrm{SWIR} 1} /\left(\rho_{\mathrm{SWIR} 1}+\rho_{\mathrm{NIR}}\right)-\rho_{\mathrm{NIR}} /\left(\rho_{\mathrm{NIR}}+\right.\right. \\ \left.\left.\rho_{\mathrm{R}}\right)-\rho_{\mathrm{G}} /\left(\rho_{\mathrm{G}}+\rho_{\mathrm{SWIR} 1}\right)\right] /\left[2 \rho_{\mathrm{SWIR} 1} /\left(\rho_{\mathrm{SWIR} 1}+\rho_{\mathrm{NIR}}\right)+\right. \\ \left.\rho_{\mathrm{NIR}} /\left(\rho_{\mathrm{NIR}}+\rho_{\mathrm{R}}\right)+\rho_{\mathrm{G}} /\left(\rho_{\mathrm{G}}+\rho_{\mathrm{SWIR} 1}\right)\right], $

(4)

$ \mathrm{SI}=\frac{\left(\rho_{\mathrm{SWIR} 1}+\rho_{\mathrm{R}}\right)-\left(\rho_{\mathrm{NIR}}+\rho_{\mathrm{B}}\right)}{\left(\rho_{\mathrm{SWIR} 1}+\rho_{\mathrm{R}}\right)+\left(\rho_{\mathrm{NIR}}+\rho_{\mathrm{B}}\right)}。$

(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]}。该算法主要公式为

$ \mathrm{LST}=\gamma\left[\varepsilon^{-1}\left(\varphi_{1} L+\varphi_{2}\right)+\varphi_{3}\right]+\delta, $

(6)

$ \gamma=T^{2} /\left(b_{\gamma} L\right), $

(7)

$ \delta=T-T^{2} / b_{\gamma}, $

(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参数定标得到。φ₁、φ₂、φ₃用下列公式计算：

$ \varphi_{1}=1 / \tau, $

(9)

$ \varphi_{2}=-L \downarrow-L \uparrow, $

(10)

$ \varphi_{3}=L \downarrow, $

(11)

式(9)—(11)中，τ、L↓、L↑分别为大气透过率、大气下行辐射和大气上行辐射，这些参数可以由网站http://atmorr.gsfc.nasa.gov查询而来；T为热红外波段在传感器处的亮度温度，其公式为

$ T=K_{2} / \ln \left(K_{1} / L+1\right), $

(12)

式(12)中，K₂、K₁为常数，对于TM 6、ETM+ 6、TIRS 10波段，K₁的数值分别为607.76, 666.09, 774.89，单位为W·m^-2·sr^-1·μm^-1，K₂数值分别为1 260.56, 1 282.71, 1 321.08，单位为K；ε为地表比辐射率，本文使用覃志豪等^[25]提出的方法进行反演。

1.3.2 遥感生态指数构建 1.3.2.1 指标分量标准化处理

上述公式计算的4个指标量纲是不一致的，因此，在计算遥感生态指数之前需要进行标准化处理，将各指标的量纲都统一到[0, 1]之间。其公式为

$ N I_{i}=\left(I_{i}-I_{\max }\right) /\left(I_{\max }-I_{\min }\right), $

(13)

式(13)中，NI_i为经过标准化处理之后的像元值，I_i为指标分量在i像元处的值，I_max、I_min为该指标分量的最大和最小像元值。

1.3.2.2 指标分量的水体掩膜处理

遥感生态指数中湿度分量(Wet)研究的是与植被和土壤湿度有关的参量，因此在计算时需要将大面积的水体提取出来并掩膜掉, 以避免水域对主成分分析的载荷造成影响。本文对水体的提取采用徐涵秋^[26]提出的归一化差异水体指数(Modified Normalized Difference Water Index，MNDWI)，其公式为

$ \mathrm{MNDWI}=\left(\rho_{\mathrm{G}}-\rho_{\mathrm{SWIR} 1}\right) /\left(\rho_{\mathrm{G}}+\rho_{\mathrm{SWIR} 1}\right)。$

(14)

1.3.2.3 计算遥感生态指数

遥感生态指数(RSEI)是可以综合以上4个指标信息的参量。主成分分析方法(Principal Component Analysis, PCA)可以去除波段之间多余信息，将多波段的图像信息压缩到少数几个波段，其中第一主成分(PC1)包含了绝大多数的信息。该方法的最大优点是根据各指标对主分量的贡献值客观地确定各指标的权重，避免人为误差，将经过标准化和水体掩膜处理的指标分量重新合成一幅新的影像。使用ENVI的主成分分析模块进行计算，得出PC1和相关统计数据。然后依次将PC1进行正负值转换、标准化处理，得到遥感生态指数(RSEI)，主要公式为

$ \mathrm{RSEI}_{0}=1-\mathrm{PC} 1, $

(15)

$ \mathrm{RSEI}=\left(\mathrm{RSEI}_{0}-\mathrm{RSEI}_{0_{-} \min }\right) /\left(\mathrm{RSEI}_{0_{-} \max }-\right. \\ \left.\operatorname{RSEI}_{0_{-} \min }\right), $

(16)

式(15)中，RSEI₀为初始生态指数，PC1为ENVI计算得到的第一主成分PC1。式(16)中，RSEI数值介于[0, 1]之间，值越大，生态越好，反之越差；RSEI_{0_min}、RSEI_{0_max}为初始生态指数RSEI₀的最小值和最大值。

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年以后趋于稳定。

为定量化分析生态状况，将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 km², 同比增长14.07%；2)生态质量等级较差和差的区域面积和比例先增加后降低。2000年，二者所占面积之和为45.33 km²，约占研究区当期总面积的1.94%；2009年，面积之和为59.38 km²，约占研究区当期总面积的2.52%；到2017年，面积为53.82 km²，约占当期总面积的2.26%。

2.2 生态环境变化监测

为分析研究区域18年来生态质量的差异变化，在RSEI分级的基础上，对研究区2000年和2017年的RSEI指数做差值运算，并将结果划分为变差、不变、变好三大类，结果如图 2所示。经统计，钦州市的总面积由2000年的2 342.19 km²增长至2017年的2 382.31 km²，面积增长了约1.7%，主要是围填海(图 2紫色图斑)的原因，且95%以上的围填海区域都集中在南侧的钦州港。由图 2可知，生态环境恶化的区域主要集中在城区和钦州港东北部，城区的扩张和钦州港工业园的建设是生态恶化的主要原因。生态环境有所改善的区域分布较广，其中以城区东北部为典型，这些区域的地貌以低山丘陵为主，地物类型主要是林地和园地。

经统计(不含围填海区)，大部分区域生态质量无明显变化(图 2淡黄色图斑)，占总面积的58.21%；生态质量变差的区域(图 2橙色图斑)面积为282.75 km²，占总面积的12.72%；生态环境有所改善的区域(图 2浅绿色图斑)面积为646.44 km²，占占面积的29.07%。

2.3 城市建成区扩张对生态的影响

分别对2000年、2009年、2017年的钦州市城区和钦州港建成区进行解译，其建筑面积由2000年的17.07 km²增长到2017年的107.27 km²，增加了5倍多。为研究城市的扩张对生态质量等级重心的影响，用ArcGIS分别将3期的RSEI栅格影像转换为矢量要素点。由前文可知，建筑用地的生态质量等级大都为差和较差，因此分别提取城区和钦州港生态等级差和较差的要素点，并利用ArcGIS分别计算这些要素点的加权平均中心(紫色点为城区的，红色点为钦州港的)，如图 3所示。图 3反映了城市建筑面积扩张和低生态等级区域中心的迁移关系。平均中心迁移的方向和距离反映了低生态等级区域的分布状况，平均中心基本上随建筑用地扩张的方向迁移，以上情况进一步表明，城市的扩张对生态有一定的负面影响。从迁移距离来看，钦州港平均中心大于城区，说明钦州港区域的建设力度比城区大。

2.4 建模分析

建模结果如下：

RSEI=0.5167+0.5069NDVI+0.8191Wet－0.1019NDBSI－0.0033LST，R²=1 (2000年);

RSEI=0.9410+0.5079NDVI+0.2583Wet－0.1339NDBSI－0.0213LST，R²=1 (2009年);

RSEI=0.5413+0.5883NDVI+0.0438Wet－0.1277NDBSI－0.0080LST，R²=1 (2017年)。

从建模情况来看，3个方程的R²均为1，模型拟合度高，建模过程中，在4个变量都保留的情况下，标准误差RMSE的值最低，说明这4个因子都是必要的。其中，NDVI和Wet的系数为正，与RSEI成正相关关系；NDBSI和LST系数为负，与RSEI成负相关关系，与主成分分析的结果相一致。

散点图能直观地反映自变量与因变量的相关性，以2017年为例(图 4)，NDVI、NDBSI与RSEI之间的线性关系较强，Wet、LST对RSEI的线性趋势较弱。从2017年的建模方程和图 4a来看，NDVI是影响研究区生态状况的关键因子，建模系数的绝对值远大于其他因子，散点图趋势线的斜率达到0.68，这表明当地植被一旦遭到规模性破坏，其生态环境将迅速恶化。

3 结论

本研究基于RSEI，对钦州市2000－2017年的生态环境状况进行定量分析评价，从时间和空间两个角度直观地表达了近20年来钦州城市扩张的生态效应。结果表明：1)2000—2017年间，钦州市RSEI均值均在0.7左右，生态等级以良(RSEI 0.6-0.8)为主，且2017年的遥感生态指数均值最大，表明钦州市的生态环境质量在改善，城市的扩张给钦州市整体生态状况带来的负面影响不大。从这方面来看，钦州市仍有很大的发展潜力。2)从主成分分析和建模方程的系数来看，绿度对钦州市的RSEI贡献程度最大，是影响生态的关键因子。3)从城市建成区扩张对低生态等级平均中心的影响结果来看，平均中心的迁移方向随城市扩张方向移动，从迁移距离来看，钦州港的建设力度大于城区。4)生态等级的分布和变化监测结果表明，部分城市建设区域的生态等级由差变好，建设初期破坏的自然生态表面后期得到了修复，表明恢复植被覆盖的园林绿化工程能有效地改善生态环境状况。