施婷婷,徐涵秋,*,孙凤琴,陈善沐,杨绘婷

1 福州大学环境与资源学院，空间数据挖掘与信息共享教育部重点实验室，福州 350116 2 福州大学遥感信息工程研究所，福建省水土流失遥感监测评估重点实验室，福州 350116 3 福建省水土保持试验站，福州 350001

20世纪90年代以来,随着遥感和地理信息技术的迅速发展和广泛应用,生态环境质量监测与评价取得了大量成果,评价理论、方法及技术手段不断丰富和完善,尤其是遥感对地观测技术为快速、准确、实时地开展多尺度生态环境动态监测和评价提供了有效手段[1- 3]。如美国的全美自然保护区生态动态监测[4]、我国环保部开展的“全国生态环境十年变化(2000—2010)遥感调查与评估”,以及目前正在开展的国家生态保护红线划定工作,都是以遥感作为重要的支撑技术。

目前,应用遥感技术开展生态环境评价研究主要包括以下3个方面：(1) 基于景观分析的生态环境评价,研究者基于土地利用/覆盖的分布特征及变化趋势等,借助景观指数来开展生态效应研究。如Estoque等[5]利用遥感数据对1988、2002和2016年的La Mesa流域进行土地利用/覆盖分类,在此基础上定量化评价景观格局的变化,以探讨流域生态环境变化的影响因素;Jaafari等[6]基于遥感影像解译和景观生态学方法的综合运用定量分析景观格局的动态变化,从而实现城市化进程的环境影响评价;Zhang等[7]基于土地利用/覆盖变化数据,从景观格局角度开展艾比湖湿地保护区和甘家湖梭梭林自然保护区的生态环境效应研究。(2) 基于评价指数的生态环境评价,它通过构建不同的指标体系来综合评价生态环境效应。如Liang等[8]基于物理环境变量(土地利用/覆盖数据、地表温度和归一化植被指数等)和社会经济变量(人口密度、收入和受教育水平等)综合评价美国印第安纳州首府的城市环境质量;Wang等[9]基于压力-状态-响应(P-S-R)框架模型,建立了1995—2015年北京市生态安全评价指标体系;国家环境保护部的《生态环境状况评价技术规范》则采用生物丰度指数、植被覆盖指数、水网密度指数、土地胁迫指数和污染负荷指数等5个指标来综合反映区域生态环境状况[10]。(3) 基于遥感信息的空间化生态环境评价,通过遥感影像制作区域生态图件,实现生态环境变化空间可视化,而不仅仅用一个量化的数值来代表区域的生态环境。徐涵秋[11]选取了完全基于遥感信息的湿度、绿度、热度、干度4个指标来构建遥感生态指数(RSEI),实现区域生态环境质量快速准确的时空变化分析;Zhang等[12]选取了基于遥感数据的不透水面、归一化植被指数、地表温度和缨帽变换的亮度和绿度分量这5个指标构建生态评价模型,开展珠江三角洲地区生态环境变化的空间动态监测;芦颖等[13]选取了植被盖度、土壤指数、土壤湿度等遥感信息来构建综合评价模型,分析1990—2015年间贵州省乌江流域生态环境质量时空变化特征。相比于其他两类的评价方法,基于遥感信息的空间化生态环境评价方法以其指标易获取、空间可视化、快速简便等优点,已被广泛应用[14-17]。值得注意的是,目前这一类生态环境评价大多数是集中在中分辨率影像上(如Landsat系列数据),而在中高或高分辨率影像的应用则很少,主要是因为中高或高分辨率影像的波段数少,通常只有可见光和近红外的4个波段,缺乏以上生态环境遥感评价技术所需要的中红外波段。由于高分辨率影像可以更好地揭示地表生态环境细节[18],因此,如何利用只有可见光和近红外4个波段的中高或高分辨率遥感影像来进行区域生态质量监测与评价是当今基于遥感信息的空间化生态环境评价技术必须解决的一个关键问题。

本文以福建敖江流域的贵安开发区为重点研究对象,采用基于遥感信息的空间化生态环境评价技术中最具代表性的遥感生态指数(RSEI),构建了可用于只有4个波段的中高分辨率影像的生态质量指标,并将其应用于定量分析建筑用地的时空变化及其产生的生态效应,以期为细尺度的区域生态环境监测与评价提供一种新的技术。

1 研究方法

1.1 研究区概况

敖江位于福建省东北部,为福建省第六大河流,是闽东独立水系,其中游修建的山仔水库被定为福州市的第二水源。敖江流经的连江段,近年来随着经济的快速发展,沿岸大片的土地被开发,大量的泥沙排入江中,特别是近年来新发展的贵安开发区,其大规模的开发建设造成水土流失问题严重,对水源地的水质形成了潜在的威胁。

本文以建设项目引发的生态变化区域为重点研究对象,最终确定了涵盖山仔水库至贵安开发区的敖江段及其沿岸地区为主要研究区(图1左,面积451.98 km2),其中大规模开发的贵安新区是本研究的重点区域(图1右)。

图1 研究区范围Fig.1 Study area

1.2 遥感数据源及其预处理

本文选用2幅中高分辨率遥感影像作为主要数据源,分别为2010年8月3日的ALOS影像和2016年9月3日的GF- 1影像(表1)。两幅影像季相相近,均为夏季,因此,植被具有相近的物候和生长状态,保证了实验结果的可比性。

表1 主要数据源

由于ALOS和GF- 1影像缺少计算RSEI生态指数所需的热红外波段,因此选用Landsat影像的热红外波段来补充。但由于缺乏与所用ALOS和GF- 1影像同一年份的Landsat影像,本文选用了它们前一年份的Landsat影像,时间分别为2009年6月6日的Landsat 5 TM 6波段影像和2015年9月27日的Landsat 8 TIRS 10波段影像。

构建RSEI的湿度分量需要选取多幅ALOS和GF- 1影像及其同步或近同步的Landsat影像(表2),因此本文还选用了江苏连云港和甘肃酒泉市的影像来补充。采用多个地理位置的影像有利于所推导的湿度分量有较好的通用性。

表2 辅助数据源

为了减少不同时相影像在地形、光照等方面的差异以及保证影像间空间叠加分析的准确性,所有的ALOS和GF- 1影像需要经过正射校正和辐射校正。首先进行影像的正射校正。由于购买的GF- 1影像是Level 1A级别,因此可基于影像自带的RPC文件,利用ENVI软件进行RPC正射校正。而ALOS影像缺乏自带的RPC文件,因此需要自定义RPC进行正射校正。辐射校正采用的是中国资源卫星应用中心官方网站(http://218.247.138.119/CN/index.shtml)提供的GF- 1辐射定标参数,以及Bouvet等[19]提供的ALOS辐射定标参数,并采用Chavez的COST大气校正模型[20]。对GF- 1和ALOS影像进行辐射校正,将影像的亮度值转换为传感器处反射率,其公式为：

L= gain ·Q+ bias

(1)

ρ= [π·(L-Lh) ·d2] /(Es·cosθs)

(2)

式中：L是传感器接收到的辐射亮度;gain和bias8分别为增益值与偏置值;Q为亮度值(DN);Lh为大气影响修正值,可以通过最暗像元法获得[20];d为日-地天文单位距离;Es为大气顶部平均太阳辐照度;θs为太阳天顶角。以上参数可从影像头文件中获得。

由于从USGS网站下载的Landsat 5 TM和Landsat 8 OLI影像均是L1T数据产品,已经过系统辐射校正和地面控制点的几何校正,并且使用DEM数据进行地形校正,因此可直接用于叠加分析。采用Landsat用户手册的模型[21]和Chavez的COST大气校正模型[20]将各个Landsat影像的亮度值均转换为传感器处反射率。

1.3 建筑用地信息提取

由于ALOS、GF- 1影像缺乏中红外波段,因此无法将传统的适用于中等空间分辨率影像的各种建筑用地指数应用于该类影像。为此,本文参考Xu[22]基于规则提取不透水面的思路,通过分析各种地类在ALOS、GF- 1影像中的光谱特征,最终建立了基于ALOS、GF- 1影像的建筑用地提取规则和模型(图2)。其中选取了归一化水体指数NDWI[23]和归一化植被指数NDVI参与模型构建,计算公式如下：

NDWI = (ρ2-ρ4)/(ρ2+ρ4)

(3)

NDVI = (ρ4-ρ3)/(ρ4+ρ3)

(4)

式中：ρi(i=1, 2, 3, 4)分别代表ALOS或GF- 1影像的第1、2、3和4波段的反射率。

图2 ALOS和GF- 1影像的建筑用地提取规则和模型Fig.2 Rules and models used for extracting built-up land feature from ALOS and GF- 1 images

1.4 RSEI指数计算

遥感生态指数RSEI是近年提出的完全基于遥感信息的生态评价指数[11],已被广泛应用于不同生态条件的区域和城市[14],[24-25]。该指数耦合了与生态密切相关,且人类能够直接感知的湿度、绿度、热度和干度4大指标。RSEI指数各指标可通过以下遥感指数或参量求取：

(1) 湿度指标：通常用缨帽变换中的湿度分量来代表。但目前尚无针对ALOS和GF- 1数据的缨帽变换算法及其湿度分量。因此,本文基于Landsat系列数据缨帽变换的湿度分量,及其与ALOS、GF- 1数据的拟合关系来求出它们的湿度分量。

图3为求取GF- 1影像湿度分量系数的流程图。首先,从表2江苏连云港市和甘肃酒泉市这2幅GF- 1与Landsat 8的同步影像对中合计选取了一一对应的40384个随机样本点,利用Baig等[26]的算法先求出Landsat 8影像缨帽变换的湿度分量,并将其作为因变量,然后以GF- 1各波段作为自变量,采用多元线性回归对二者进行拟合,求出GF- 1湿度分量的经验系数,并将其单位化,即获得GF- 1湿度分量的系数(方程通过P<0.001的显著性检验)：

WetnessGF- 1= -0.1591ρ1+0.7574ρ2-0.6311ρ3-0.0523ρ4-0.003 (R=0.829)

(5)

图3 GF- 1影像求取湿度分量系数的流程图 Fig.3 Flowchart of the derivation of the wetness-component coefficients for GF- 1 imagery

同理,利用Crist[27]的Landsat 5 TM缨帽变换湿度分量与ALOS数据间的拟合关系求取ALOS影像缨帽变换湿度分量(方程通过P<0.001的显著性检验)：

WetnessALOS=-0.1912ρ1+0.6081ρ2-0.7381ρ3-0.2210ρ4+0.003 (R=0.807)

(6)

(2) 绿度指标：选用与植物生物量、叶面积指数和植被覆盖度都密切相关的归一化植被指数(NDVI)来表示,见公式(4)。

(3) 热度指标：由于ALOS和GF- 1影像缺乏热红外波段,采用TM 6波段和TIRS 10波段反演的地表温度(LST)代表热度指标。本文采用Jiménez-Muoz等[28-29]提出的单通道算法来反演地表温度。其主要公式为：

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

(7)

(8)

式中：LST为地表温度;对于TM 6波段和TIRS 10波段,参数bγ分别为1256和1324;γ和δ是基于Planck函数的2个参数;ε为地表比辐射率,对于TM 6波段采用Sobrino等提出的基于植被覆盖度的算法[30],对于TIRS 10波段,其取值可见参考文献[31-32];ψ1、ψ2、ψ3可用下式计算：

ψ1= 1/τ,ψ2 =-L↓-L↑/τ,ψ3=L↓

(9)

式中：τ为大气透过率,L↑和L↓为大气上行和下行辐射强度,可根据影像的时间以及所处地理位置和季相,通过MODTRAN模拟得到;T为传感器处的亮温,公式为：

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

(10)

式中：K1和K2分别为定标参数,对于TM 6波段,K1=607.76 W m-2ster-1μm-1、K2=1260.56 K,对于TIRS 10波段,K1=774.89 W m-2ster-1μm-1、K2=1321.08 K。L为TM 6波段或TIRS 10波段的传感器处光谱辐射值,对于TM 6波段,L可以从公式(1)求取,对于TIRS 10波段,L可以用下式求取：

L=ML·Q+AL

(11)

式中：ML和AL分别为TIRS 10波段的乘法调整系数和加法调整系数(可从影像的头文件中获得)。

(4) 干度指标：由于缺乏中红外波段,RSEI中采用的代表干度的IBI建筑指数[33]、SI裸土指数[34]及二者的合成指数NDBSI[11]等3个干度指标都无法应用,因此,通过考察类似的指标,考虑采用不需要中红外波段的土壤盐渍度指数(Soil Salinity Index)[35- 38]来表示。该指数不仅能够增强土壤盐渍化状况,也能够有效地增强建筑和裸土等地表信息。表3为所用的2幅Landsat影像的各个土壤盐渍度指数与其3个干度指数之间的相关系数R。研究表明：SI-D、SI1、SI3和S5这4个土壤盐渍度指数与干度指数的相关性较好,均为正相关,2幅影像的R的均值均高于0.84。为了进一步确定上述4个土壤盐渍度指数在ALOS和GF- 1影像的适用性,分别计算出ALOS和GF- 1影像这4个土壤盐渍度指数,然后将它们与近同步的Landsat影像的3个干度指数进行相关性分析(表4),发现对于ALOS和GF- 1影像,Douaoui等[37]提出的盐渍度指数(SI-D)与各个干度指标的相关系数R的平均值最大,因此本文选其作为干度指标,公式为：

SI-D = (ρ1×ρ3)1/2

(12)

RSEI采用第一主成分(PC1)来构建指数[39]。主成分变换通过对特征光谱空间坐标轴的旋转来去掉各指标间的相关性,从而将最主要的信息集中到第一主成分上。其最大优点就是集成各指标的权重不是人为确定,而是根据各个指标对各主分量的贡献度来自动、客观地确定,从而避免人为设定权重造成的结果偏差。

表3 Landsat影像的土壤盐渍度指数与干度指数之间的相关系数R

表4 ALOS/GF- 1影像的土壤盐渍度指数与近同步的Landsat影像的干度指数之间的相关系数R

Table 4 The correlation coefficients(R)between the soil salinity indices of ALOS or GF- 1 images and the dryness indices of near synchronous Landsat images

文献Literature指数IndexALOS-TM GF-1-OLI IBISINDBSIIBISINDBSI平均AverageDouaoui等[37]Salinity Index (SI-D)0.681 0.765 0.739 0.726 0.834 0.799 0.757 Salinity Index 1 (SI1)0.680 0.772 0.742 0.714 0.824 0.788 0.753 Salinity Index 3 (SI3)0.677 0.767 0.738 0.710 0.820 0.784 0.749 Abbas和Khan[38]Salinity Index (S5)0.663 0.729 0.710 0.742 0.834 0.807 0.748

由于以上4个指标量纲不统一,需将它们进行正规化,将它们的值都统一到[0, 1]之间。 经过正规化后的4个指标就可以用来计算PC1。为使PC1大的数值代表生态条件好,进一步用1-PC1获得RSEI：

RSEI=1- {PC1[f(Wetness, NDVI, LST, SI-D)]}

(13)

为了便于比较,进一步将RSEI进行归一化,使其值介于[0, 1]之间。其值越接近1,表示生态质量越好。

2 结果与分析

2.1 建筑用地的时空变化

采用以上ALOS和GF- 1影像的建筑用地提取模型获得研究区2010、2016年的建筑用地信息, 并对提取结果再进行适当的人工修改(图4)。参照ALOS全色影像(2.5 m)和GF- 1全色影像(2 m),分别采用300个随机样点对提取结果进行精度验证。结果表明,2个研究年份的建筑用地提取的总精度及各类用地的使用者精度和生产者精度均在88%以上(表5),满足精度要求。

为了更好地了解研究区的建筑用地变化,对研究区和贵安开发区所提取的建筑用地信息进行统计分析(表6)。结果表明,研究区建筑用地面积在建设前后快速增长,从2010年的7.43 km2增长到2016年的17.79 km2,共增加了10.36 km2,增幅达139%。其中2010—2016年间贵安开发区的建筑用地变化面积占整个研究区的86%,可见该区域建筑用地扩展主要是由大规模的建设项目引起的。

图4 2010和2016年研究区遥感影像及其对应的建筑用地提取影像Fig.4 Remote sensing images and the corresponding built-up land maps of the study area in 2010 and 2016

年份Year数据类型Data type类别Class验证数据Reference data建筑用地Built-up land非建筑用地Non-built-up land总计Row totals使用者精度/%User′s accuracy 2010分类数据建筑用地4544991.84非建筑用地624525197.61列合计51249300生产者精度/%88.2498.39总精度Overall accuracy =96.67%;Kappa=0.88002016分类数据建筑用地5265889.66非建筑用地723524297.11列合计59241300生产者精度/%88.1497.51总精度Overall accuracy=95.67%;Kappa=0.8619

图5为2010—2016年引发建筑用地发生变化的大型建设项目分布位置图。近6年来,贵安开发区呈现爆发式发展,特别是在贵安新天地园区中,一系列大型配套项目如国际创新中心贵安基地、贵安花海拓展公园、贵安水世界、贵安欢乐世界、奥特莱斯购物广场、福州一中贵安学校、福州百姓文化长廊、世纪金源大饭店等以及一批低密度住宅、花园洋房、高层公寓等住宅区相继建成,使该区一跃成为“海西首席旅游商住综合体”。此外,贵安通往福州市区的主要道路沿线也有大面积住宅区建成。

表6 2010—2016年建筑用地变化

图5 2010—2016年建筑用地主要变化区的大型建设项目位置示意图Fig.5 Map showing the locations of large construction projects developed in the years 2010—2016

2.2 生态质量的时空变化

首先根据推导出的湿度分量系数反演出研究区2010年ALOS和2016年GF- 1影像的湿度分量图(图6)。图中可以看出,水体具有最高的湿度(红色),植被次之(绿色),而建筑用地最低(蓝色),说明本文建立的湿度分量可以客观反映地表湿度的实际情况。

图6 2010和2016年湿度分量分级图Fig.6 3-leveled wetness component maps of the study area in 2010 and 2016

表7为遥感生态指数RSEI的主成分分析数据。从表中可知,在这2个年份中,PC1的特征值所占的比例都大于80%,说明用PC1可以较好地代表4个指标变量。在这4个指标中,湿度和绿度均为正值,说明它们对生态的贡献起正向作用;而热度和干度均为负值,说明它们对生态的贡献起负向作用,这与实际情况相符。

表8为研究区RSEI的相关统计数据。从表中可知,整个研究区的RSEI均值从2010年的0.813下降到2015年的0.765,下降了6%,说明生态质量总体下降。而就贵安开发区而言,其生态质量下降得更明显,RSEI均值从建设前的0.787下降到建设后的0.689,下降了12.5%。从4个分指标变化情况来看,对于整个研究区,对生态有利的湿度和绿度指标的均值在2010—2016年间都表现为下降,而对生态起负面作用的热度和干度指标的均值都表现为上升。这4个指标的各自表现也说明了整个研究区的生态质量呈下降趋势。而对于贵安开发区,其变化趋势与整个研究区类似,但总体变化幅度更大,这也使得贵安开发区的生态质量下降得更为明显,其下降幅度是整个研究区的2倍。总的看来,所构建的生态指数RSEI可以很好地综合反映4个指标的变化情况。

表7 主成分分析数据

表8 4个指标分量和RSEI指数的变化

RSEI自2013年提出以来,已被广泛应用,从数十个地区的应用情况来看,该指数具有较好的鲁棒性,其结果也都有很好的可比性。如沙漠地区一般<0.35[40],城市一般在0.4—0.6[15,41],林区一般>0.65[42- 43]。本次研究区以林区为主,2个年份的RSEI值都在0.65以上,与现有的研究(如南瓮河自然保护区[42]、长白山自然保护区[43])基本一致。

图7为基于RSEI生态指数反演的研究区2010、2016年生态状况图,图中颜色从红到绿代表生态从差到优。建筑用地主要呈红色调,植被呈绿色调,水体根据浑浊程度呈黄到红色调。

图7 研究区2010、2016年RSEI生态状况图Fig.7 RSEI maps of the study area in 2010 and 2016

进一步将各年份的生态指数以0.2为间隔分为5级,分别代表差、较差、中等、良、优5个生态等级,然后据此对建设前后2个年份的生态指数影像进行分级(表9)。统计表明：2010—2016年,研究区的生态级别为差到中级(1—3级)所占的面积比例从5.51%上升到8.31%,而优良等级(4—5级)所占的面积从94.49%下降91.69%,表明整个研究区的生态质量虽有下降但幅度较小。然而,贵安开发区则明显不同,其生态级别为差到中级(1—3级)所占的面积比例从9.29%上升到20.9%,而优良等级(4—5级)所占的面积从90.71%下降79.1%,表明开发区的生态质量明显下滑。

表9 2010—2016年各RSEI级别的面积变化/km2

图8 基于红绿法的2010—2016年RSEI变化检测 Fig.8 Change detection of RSEI between 2010 and 2016 based on red-green differencing method

RSEI除了能够作为刻画区域生态状况的一个定量指标外,还可用于区域生态时空变化的可视化,有利于不同时期生态状况的时空对比。利用红绿法对研究区开发前后生态质量的空间变化进行差值变化检测(图8)。从图8看出,代表生态条件变差的大片区域主要集中分布在贵安开发区,特别是在贵安新天地园区,这也与建筑用地空间扩展情况一致。从变化检测结果来看(表10),2010—2016年间,研究区等级上升(变好)的地区面积为28.52 km2,占总面积的6.3%,而生态等级下降(变差)的地区面积为153.47 km2,占33.9%,大大超过前者。

2.3 建筑用地空间变化与生态变化的关系

为了揭示建筑用地与生态质量二者时空变化的定量关系,利用RSEI可刻画空间变化的优势将二者进行回归分析。为避免不同空间单元尺度对研究结果造成影响,本文按照50 m×50 m、100 m×100 m、300 m×300 m和500 m×500 m这4种单元网格分别对研究区2010、2016年的RSEI和建筑用地影像进行系统采样,求取每个网格在两年份间生态质量的差值。同时,计算出两年份间各单元网格中的建筑用地面积占比变化情况,然后对二者进行回归分析,以分析区域建筑用地空间变化与生态变化之间的定量关系。

回归结果表明(图9),RSEI与建筑用地面积占比在4种单元网格中均呈现出显著的线性负相关关系(R2均大于0.8),表明建筑用地占比越大,其生态质量下降越显著。以R2最大的100 m×100 m网格的回归方程为例,若建筑用地面积占比每增加10%,其RSEI值将下降0.041。也就是说,如果要将贵安开发区生态质量控制在不低于良好等级,即RSEI值不小于0.6,则建筑用地面积的上升空间只有10.33 km2,一旦超出这一界限,则生态质量就会落入中等级别。

表10 2010—2016年研究区生态质量变化检测

图9 不同采样网格的建筑用地面积占比变化与生态质量变化的回归关系Fig.9 Relationship between built-up land and ecological quality tested with different sampling grids

3 结论

本研究通过多元回归分析构建湿度分量和引入盐渍度指数,首次将遥感生态指数(RSEI)运用在只有可见光和近红外4个波段的中高分辨率遥感影像上,解决了开展生态环境遥感评价技术依赖中红外波段的问题,实现了区域生态质量在空间细尺度上的定量综合评价。

总的看来,贵安开发区的建设是引起研究区生态质量总体下降的关键原因。贵安开发区的建筑用地面积在2010—2016年间共增加了8.94 km2,增幅达266%,占整个研究区建筑用地变化面积的86%。建筑用地面积的增长导致了开发区范围生态质量的显著下降,表现在其生态指数RSEI均值在建设前后下降了12.5%。建筑用地面积占比与生态质量呈显著的负相关关系,建筑用地面积占比每增加10%,其RSEI值将下降0.041。显然,贵安开发区的建设,破坏了该区原有以植被为主的自然生态,造成了区域生态质量的显著下降。因此,在区域经济发展的决策规划中,要融入生态的理念,严格控制建筑用地面积,有效地保护青山绿水,实现区域的健康可持续发展。

对一个区域的生态研究,一般应采用同一数据源的遥感影像,但本文因受到可用影像数据源的限制,不得不采用ALOS AVNIR- 2和GF- 1 PMS1这2种不同的传感器影像,因此,可能会对研究结果会带来一定的影响,如不同传感器之间的系统差异等。