实验十三 生态安全评价

1 背景及目的

        生态安全是生态系统完整性和健康性的整体反映，完整健康的生态系统具有调节气候净化污染、涵养水源、保持水土、防风固沙、减轻灾害、保护生物多样性等功能。维护生态安全对于人类生产、生活、健康及可持续发展至关重要。随着城市化进程的不断推进，生态系统受到来自人类社会日益严重的威胁，保障生态安全作为迫切的社会需求已经成为人类社会可持续发展面临的新主题。因此，对区域生态安全水平进行评价，为改善区域生态问题的对策与措施提供参考具有重要的现实意义。

       通过本实验希望能够达到以下目的：

       ①掌握生态安全评价指标体系的构建与计算

       ②掌握景观指数移动窗口计算方法

       ③掌握突变模型在指标集成中的应用

2 实验步骤

2.1 数据的预处理

2.1.1 设置工作空间

添加数据：

地理处理-->环境设置：

对工作空间进行设置主要为方便管理过程中产生的数据，也方便自己找得到想要用的数据，更改到当前实验文件夹：

2.1.2 栅格数据的统一处理：栅格重采样

数据管理工具-->栅格--栅格处理-->重采样

这一步骤的作用在于本实验的统一栅格数据为10000m精度，但原数据土地利用数据为30m，因此采用重采样工具对土地利用数据进行处理，将分辨率更改为1000×1000：

2.2 生态服务能力指标计算

不同土地利用类型对应不同生态服务价值系数，利用重分类工具更改栅格中的值：

Spatial Analyst 工具-->重分类-->重分类

更改后如图：

2.3 生态敏感性指标计算

2.3.1 距道路距离的计算

Spatial Analyst 工具-->距离-->欧氏距离

在环境中对输出范围和坐标系进行设置：

结果(我也不懂为什么图是歪的)：

以区域面数据为范围对得到的距道路距离栅格进行提取：

Spatial Analyst 工具-->提取分析-->按掩膜提取：

2.3.2 距城市距离

这一步骤可以得到土地利用的矢量数据，转换工具-->由栅格转出-->栅格转面：

右键打开土地利用类型的属性表，按属性筛选出对应值为80的数据：

将刚才选中的数据导出：

对图层进行编辑：

 选择面积较大的部分作为城市区域，再将其他区域进行删除，最后保存编辑结果：

打开建设用地属性表，刚才选择的地区都会进行高亮显示，然会选择结果进行反选删除不要的数据：

继续利用欧氏距离和按掩膜提取将处理过的建设用地设为输入栅格数据，得出距城市距离栅格数据结果：

结果：

右键属性调整一下色带：

结果：

2.3.3 植被覆盖

对土地利用数据进行重分类：Spatial Analyst 工具-->重分类-->重分类

结果：

 2.3.4 坡度

Spatial Analyst 工具-->表面分析-->坡度

结果：

右键土地利用，导出数据：

2.4 生态组织架构

打开fragstats，把上一步生成的文件导入：

 

创建一个记事本文件，文件后缀为fcd，输入如下内容，老实英文状态下手敲，不要文字识别，容易后面运行报错：

导入：

 设置分析参数：

勾选景观指数：

 

 

运行：

 报错了，查看下日志：

Error: No level selected for the analysis, please select at least one level.

提示没有为分析选择级别，倒回去看步骤，发现确实少勾选每一个条件：

再次运行，成功了，系统字自动生成一个文件夹：

文件位置： 

由于在Fragstats中空值会显示为-999，为不影响归一化结果，需要用“栅格计算器将空值去除

先把Fragstats中生成的6个景观指数加载进ArcMap：

利用栅格计算器去除计算出现的异常值：

Spatial Analyst 工具-->地图代数-->栅格计算器

 另外五个都是同样的操作：

2.5 突变模型

利用突变模型对指标体系进行架构

2.5.1 剔除水体

打开土地利用类型属性表，筛选出水体对应值为60的：

右键图层将水体数据导出：

取消上一步的选择：

使用分析工具中的擦除，将水体数据从土地利用数据中擦除：

分析工具-->叠加分析-->擦除

利用区域面的范围对得到的土地利用类型进行裁剪：

分析工具-->提取分析-->裁剪

批量处理14个指标数据--->按掩膜提取

Spatial Analyst 工具-->提取分析-->按掩膜提取

2.5.2 指标归一化

打开python：

写入代码前需要知道哪些是正向指标，哪些是负向指标，二者代码存在区别：

 

我按上图顺序将代码附上：

注意正向指标与负向指标的公式语句存在差别，代码我也分开放了：

//1 NDVI1  正向
import arcpy
from arcpy import env
from arcpy.sa import *rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_NDVI1.tif"//需要归一化的数据所在位置，一定是上一步进行掩膜提取后的数据
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\NDVI.tif"//归一化后的数据存放位置maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：0.89999997615814    原栅格最小值：0.11599999666214//2 NPP2  正向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_NPP2.tif"//NNP2
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\NPP2.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：2049.6000976563  原栅格最小值：0//3 生态服务价值   正向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_Land_se1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\生态服务价值.tif" maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：3149   原栅格最小值：0//4 距城市距离   正向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_城市距离掩膜提1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\距城市距离.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：60373.8359375   原栅格最小值：0//5 距道路距离   正向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_Extract1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\距道路距离.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：20396.078125   原栅格最小值：0//6 植被覆盖   正向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_植被覆盖1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\植被覆盖.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：4   原栅格最小值：0//14 斑块连通度 正向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_斑块连接度co1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\斑块连通度.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(Raster(rawRaster)-float(minValue))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：100    原栅格最小值：20.166975021362

//7 降雨量   负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_降雨量1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\降雨量.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：13948.84375    原栅格最小值：11076.044921875//8 坡度   负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_Slope_D1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\坡度.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：30.61713218689    原栅格最小值：0//9 斑块密度   负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_板块破碎度1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\斑块密度.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!"
//我的结果：原栅格最大值：0.77777779102325    原栅格最小值：0.11111111193895//10 边缘密度   负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_边缘密度ED1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\边缘密度.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!" 
//我的结果：原栅格最大值：13.333333015442   原栅格最小值：0//11 面积加权平均形状指数  负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_面积平均加权形1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\面积加权平均形状指数.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!" 
//我的结果：原栅格最大值：2   原栅格最小值：1//12 香农多样性指数   负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_香农多样性SH1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\香农多样性.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!" 
//我的结果：原栅格最大值：1.6769877672195   原栅格最小值：0//13 斑块平均距离 负向
import arcpy
from arcpy import env
from arcpy.sa import *
rawRaster="E:\\DengSL zuoye\\exp13\\exp13\\data_exp13\\Extract_斑块平均距离E1.tif"
NormalizedRasterpath="E:\\DengSL zuoye\\exp13\\guiyihua\\斑块平均距离.tif"
maxValueDS=arcpy.GetRasterProperties_management(rawRaster,"MAXIMUM")
maxValue=maxValueDS.getOutput(0)
print"原栅格最大值："+str(maxValue)
minValueDS=arcpy.GetRasterProperties_management(rawRaster,"MINIMUM")
minValue=minValueDS.getOutput(0)
print"原栅格最小值："+str(minValue)
NormalizedRaster=(float(maxValue)-Raster(rawRaster))/(float(maxValue)-float(minValue))
NormalizedRaster.save(NormalizedRasterpath)
print"归一化成功!" 
//我的结果：原栅格最大值：2828.4272460938    原栅格最小值：0

2.5.3 突变模型

突变模型对于不同数量变量有不同的归一化公式

观察上面两个图，先算B级指数，即活力B1、经济压力B3、自然生态B4、景观破碎度B5、景观连通度B6。

如何去计算这几个指数？则根据指标数量（控制变量数量）来带入公式，此处以景观破碎都B5为例，由图可知景观破碎度B5具有斑块密度C9、边缘密度C10、面颊加权平均形状指数C11、香浓多样性指数C12这四个控制变量，那么带入的公式为：

              

接下来打开栅格计算器：

Spatial Analyst 工具-->地图代数-->栅格计算器

结果：

以此类推得活力B1：

社会经济压力B3：

自然生态B4：

景观连通度B6：

生态服务能力：

 生态敏感性A2：

生态安全指数：

所有数据计算完毕，调整细节出图：