多源多时相生态脆弱性遥感定量评价方法研究

杨明芬吴旭

（1.西藏自治区科技信息研究所，西藏 拉萨 850000；2.成都理工大学，四川 成都 610059）

生态环境是人类生存和发展的重要物质基础之一，对生态环境脆弱性变化展开研究，是制定区域可持续发展战略的重要前提。对生态环境进行检测和评价分析，有助于及时了解环境质量变化的规律并采取科学的应对措施。生态环境脆弱性评价是区域生态综合评价的重要内容和生态规划的重要依据，近年来在地理学、生态学和环境科学研究中广泛应用。

以RS和GIS技术支持的生态环境脆弱性研究被更多的学者采用。生态环境脆弱性评价可分为定性评价和定量评价，定量方法对环境的脆弱性、稳定性、敏感性和外部环境胁迫对系统造成的影响进行定量描述。［1］利用生态数学语言提出了生态脆弱带基础判定的数学方法。［2］用统计模型对生态环境脆弱性密度分布进行了初步测算。［3］通过多源遥感数据下植被覆盖变化研究对脆弱性评价进行了研究。由于生态环境脆弱性研究涉及地学、生态学等多个领域的研究，具有很强的综合性，因此传统评价方法使用的研究方法、评价指标体系和分级标准呈现多样化。各种研究方法的重点集中在确定生态脆弱性评价指标、权重设定、指定各指标稳定度的标准方面。这些方法需要较长的时间来处理前期数据资料，研究周期较长，无法快速定期对生态环境脆弱性现状和变化趋势进行评价监测，而且评价结果不够全面。

文章在前人建立的生态环境脆弱性定量评价方程基础上，研究一种新的区域态环境脆弱性定量评价模型。基于多源多时相的卫星遥感影像数据提取植被指数、土地覆盖分类、DEM坡度与坡向数据作为评价模型输入的原始数据。采用多元非线性回归方法提取生态环境脆弱性定量评价所需的模型参数。通过评价模型的计算实现研究区生态环境脆弱性的快速、定量评价，形成脆弱性评价分级图，并对模型评价的精度和稳定性进行了验证。

1 研究方法

1.1 基本思路

文章采用生态环境脆弱性的因子相关分析方法，参考［4］刘燕华和王经民的经验公式,构建脆弱性定量评价模型：

式（1）中F为生态环境脆弱性值，U为土地覆盖类型的变化梯度，f为变化区域的植被覆盖度，其中：

式（2）中U为变化的土地类型面积；Si为该类型区域前一年土地总面积为土地覆盖类型变化梯度；i为土地覆盖类型，参考前人研究基础并结合与植被有关的特性，划分为林地、草地、湿地、水体、人工用地、沙地等六种类型；

式（3）中f是估算的指标覆盖度，DNVI是归一化植被指数，NDVImax和NDVImin分别是研究区同一年内的最大值和最小值。

通过遥感监督分类提取NDVI、土地覆盖变化范围，并通过编写IDL程序计算土地变化的面积和总面积的最大值与最小值，求出变化梯度，根据NDVI求出植被覆盖度，代入式（4）求出由于变化梯度导致的生态环境脆弱性值。在IDL环境中，对研究区的数字地形进行高程、坡度与坡位的综合分析，选取与植被相关性最大的因素代入评价方程中，计算得出脆弱性评价结果。将数字地形模型与评价结果叠加，得出区域生态环境脆弱性综合评价等级。实验流程如图1。

图1 实验流程图

1.2 实验数据

研究区位于四川省若尔盖县，四川、甘肃、青海三省交界处，地处青藏高原东部边缘，海拔3400-3900m，东西与南北最大距离约150km，面积约10620km2，地理坐标为 102°08′E 至 103°39′E，32°56′至34°19′N之间，是长江和黄河流域上游地区重要的水源涵养地，也是中国残存面积最大的高原泥炭沼泽湿地，整个研究区形成一个相对完整的地理自然区，对于研究区域生态脆弱性变化极具代表性。

1.2.1 Landsat数据。Landsat1-8系列卫星是由美国航空航天局 NASA（National Aeronautics and Space Ad⁃ministration）发射的对地观测卫星。文章采用Land⁃sat5 TM数据提取土地覆盖分类信息，并计算土地类型变化率。TM的覆盖范围为之间的所有陆地区域，东西扫描范围约183km，南北扫描范围约为170km，数据更新周期为16天，空间分辨率30m。文章获取2001-2010年研究区内的Landsat5 TM遥感数据（每年4景，共40景数据）。考虑到研究每一景影像可能存在多云覆盖的问题，我们选择的单景影像平均云量小于10%的影像数据，成像时间集中在8-9月。选择出每一年内成像质量较好的影像数据，对各景数据进行镶嵌，根据研究区矢量地图裁剪得出研究区影像数据，对裁剪的研究区数据进行大气校正、正射校正、去除薄云处理。其中投影为UTM、坐标系为WGS-1984，如图2。

图2 研究区历年影像

1.2.2 MODIS数据。文章采用MOD13A1 NDVI数据产品，共获取2001、2003、2006、2008年4个年份，16天间隔的产品数据（每年23景，共计92景影像数据，空间分辨率500m）。在ENVI软件环境中对每年数据及十年数据进行HANTS滤波处理，生成长时空序列数据集，采用波段运算工具，对NDVI的月均、年均值进行计算。对数据集进行Hants滤波处理，然后利用波段运算工具计算2001-2010年均NDVI的数据集，如图3。

图3 NDVI年均数据集

1.2.3 地形数据。本文选取的DEM数据为SRTMDE⁃MUTM（90m分辨率），提取高程数据；坡度数据为SRTMSLOPE（90m分辨率）；坡向数据为SRTMAS⁃PECT（90m分辨率），如图4所示。这三种数据用于分析其与NDVI的相关性，估算未知区域NDVI的取值，为定量评价模型提供所需参数。

图4 地形数据

1.3 实验步骤

1.3.1 计算土地覆盖类型的变化梯度。研究区典型地物类型包括：河流、湖泊、沼泽、沼泽化草甸、湿草甸、灌丛、灌木林、森林、裸地、裸岩、沙地和建筑用地等［5］。其中河流、湖泊、泥炭沼泽、沼泽化草甸和湿草甸是研究区的主要湿地类型［7］。文章参考《全国遥感监测土地利用/覆盖分类体系》、前人研究成果，并考虑模型计算的简便性，把土地覆盖类型分为：草地、林地、裸地（裸土、沙地、裸岩、积雪）、水体（河流、湖泊）、沼泽、人工用地（居民地、耕地、道路）六类［6］。为保证分类精度和信息提取的快捷，文章采用ENVI软件下的支持向量机分类（Support Vector Machine，SVM）方法对遥感图像进行监督分类，再通过人工操作修改错误分类。分类过程中结合野外调查的先验知识选择训练样本。影像中难免存在不可去除的云层干扰等因素，由此产生错误分类，文章通过与影像观测时间最接近的数据进行比较人工修改分类。把相同年份的遥感预处理数据与分类结果数据叠加，在遥感预处理影像中随机选取各类型纯净像元，采用ENVI的分类后处理工具中的误差混合矩阵进行分类精度评价，分类结果精度如表1所示（以2001年为例）。

表1 土地覆盖类型分类精度表

从各时相分类精度可以看出，总体精度在85%以上，可以满足应用需求。其中沼泽与水体、沼泽与草地、裸地与草地误分的现象明显，这主要是由于高植被覆盖沼泽与草地的光谱信息比较接近，地覆盖草地与裸露地表、沙地等裸地在图像上比较相似，从而导致误分；同时少量区域由于厚云根本上无法提出，对林地、草地也存在误分。

文章在获取2001、2003、2006、2008年土地覆盖类型的分类数据基础上，采用的分类统计工具计算每种覆盖类型相邻年份的变化值，计算出土地覆盖类型的变化梯度，对土地覆盖类型的动态变化进行分析。通过计算结果获得变化的类型、面积、像元数量，并统计出土地覆盖类型变化的百分比。

1.3.2 植被覆盖度估算。然后基于多源遥感数据的生态环境脆弱性评价，评价模型参数组要依赖于植被覆盖度和土地覆盖类型的变化梯度。但一个地区的生态环境脆弱性往往是多种因素共同的结果。研究植被指数、地形高度、坡度以及坡向变化的规律有利于区域生态环境脆弱性定量评价的准确性。文章对获取的2001-2010年NDVI数据集与研究区地形数据的相关性分析，步骤如下：

将NDVI数据集与研究区高程、坡度、坡向数据叠加，并配准；

编写程序，在研究区内指定的矩形窗口范围内生成随机采样点，获取每一年份同一经纬度的NDVI值、海拔高程、坡度值、坡向值；

通过程序自动提取对应坐标点的植被指数、值高程值、坡度和坡向值，生成NDVI-高程、NDVI-坡度、NDVI-坡向的数据记录文件，以便进一步分析它们之间的关系。

1.3.3 NDVI与地形数据的相关性分析。将采集到的10个时相NDVI数据与地形高程进行研究：利用IDL编程绘制出对应坐标点的NDVI随地形高度变化的散点图，对每一时相2000多个采样点进行统计分析。从绘制出的散点图可以观察出：NDVI值随着高度的增加，是先增加后减少，其变化规律和植被分布的垂直地带性是相符合。坡度为0的时候，ND⁃VI的值在坡度0.4-0.5区间分布均匀，这表明在草地这样的平坦地形，NDVI的值是没有反映坡度的变化规律。因此NDVI随坡度的变化表现平直，没有明显的升降变化。同时，NDVI与坡向的变化分布均匀，没有明显的规律性变化，数据与各种拟合函数相关性极低。NDVI与坡度和坡向的相关系数大都小于0.2，因此跟研究中只针对NDVI与高度值进行回归分析和函数拟合。根据散点图显示，NDVI与高度之间有明显的规律性变化，但无明显的线性或非线性函数关系，考虑对NDVI值进行某种转换，研究分析转换后的值与高程之间的相关关系，如表2所示。

表2 NDVI与地形数据的相关性分析

文章尝试对NDVI与NDVI均值差的绝对值取自然对数，然后分析转换后的值与高程的关系、转换值，通过编写程序生成与高程值的散点图，如图5所示。根据散点图的分布情况进一步分析拟合函数的形态。

图5 NDVI与高程的非线性分析

1.3.4 基于多元非线性回归的NDVI估值。对多年份数据的转换值Y与高度值的散点分布进行分析，发现散点的分布特征十分接近e的幂函数类型，因此设计一种多元非线性回归方程对数据点进行拟合，并通过IDL编程计算回归参数：Y=C*eax+b，其中x为高度值，a、b、c为待定参数。通过以上计算阿对未知点的NDVI进行估值，从而估算出研究区的植被覆盖度。图6描述了计算的待定参数形成的拟合函数与Y的关系。

拟合函数对NDVI估计值的精度检验如表2所示。

表2 NDVI估值精度检验

图6 NDVI非线性函数拟合

1.4 定量评价模型

非线性回归模型表示为：

其中f为期望函数，xn为第n个响应的回归向量或自变量向量。模型的期望响应是关于参数的非线性函数。θ表示未知参数，当分析一个特定的数据集时，假定向量xn，，n=1，…N，N是固定的，并几种考虑期望响应关于θ的相依关系。构造N维向量η（θ），其第n个向量为：ηn（θ）=f（xn，θ），n=1，2…，N；记非线性回归模型为：

对于非线性的情形，利用迭代的方法求解最小二乘估计量。利用期望函数的线性近似，通过迭代改善的初始值θ0，直到迭代无变化为止。

将期望函数f（xn，θ）在θ0处进行一阶泰勒展开：

与η（θ0）相比，现在点y。因此取更好的参数值θ1=θ0+δ0，开始另一次迭代，继续计算新的残差z1=y-η（θ），新的导数矩阵V1和新的增量。重复上述过程，直至达到收敛为止。

2 实验结果与分析

在生态环境评价中，由于评价单元不同，评价结果会存在差异,因此评价单元的选择与研究区范围密切相关［7］。文章以遥感影像中的栅格数据为载体，以每一格栅格单元为评价的最小单位，在模型计算的基础上，按照指数的数值大小进行等级划分，确定其优劣程度。把脆弱性指标值划分为5个等级：无脆弱性、轻微脆弱性、一般脆弱性、中度脆弱性和严重脆弱性。文章采用此方法对2003、2006、2008三年的评价结果进行分级，以此分级阈值作为评价，最终得到研究区的生态环境脆弱性定量评价结果，如图7所示。

图7 区域生态环境脆弱性定量评价计算结果

由计算结果可知研究区评价值主要集中在0.2-0.8之间，生态环境脆弱性严重的区域主要集中在北部和中部部分地方。这里主要是裸露地表、沼泽边缘和建筑用地、耕地等受人类活动影响较为频繁的区域。从图7可以看出，各评价等级相间分布，各等级间没有明显的界线。

3 结论

文章从少量因子出发，先抽取和植被生态方面有关的因子（土地覆盖类型变化率、植被覆盖度），研究生态环境脆弱性和它们之间的关系。再研究植被指数和地形变化的规律，逐个因子添加不仅可以使研究的复杂问题变的简单而且有助于揭示它们之间的联系。主要研究成果包括：实现了生态环境脆弱性的遥感定量研究。解决了如何从遥感影像中提取出反映生态环境脆弱性的因子，并把这种参数带入到生态环境脆弱性的方程中，从而实现生态环境脆弱性的遥感定量研究是问题的关键。从植被的角度出发把造成生态环境脆弱性的因子划分为：植被覆盖度、土地覆盖类型变化梯度、地形高度、坡度以及坡向。从遥感影像中提取出土地覆盖类型变化的梯度，利用NDVI和植被覆盖度的关系，定量表达植被覆盖度。通过模型求出生态环境脆弱度，实现生态环境脆弱性的遥感定量研究，研究具有创新意义。揭示了研究区植被指数随地形高度、坡度、坡向变化的规律。研究中对每一时相采集到2000多个具有典型意义的点，数据丰富，资料可靠。通过研究发现，NDVI和高程存在着明显的相关关系，通过IDL编程求出了具体的定量表达关系式。实验表明植被指数随地形坡度和坡向的变化规律不明显，只简单地分布在均值周围。基于多元非线性回归方法，建立了生态环境脆弱性植被遥感定量研究的模型，通过改造传统意义上的脆弱度方程，实现了参数的替代。该模型中脆弱度被定义为土地覆盖类型变化的梯度和植被覆盖度的商。通过反映生态环境脆弱性本底值的方法来评价划分生态环境脆弱性等级。