发明内容
针对现有技术的不足,本发明提供一种生态承载力的遥感估算方法,考虑空间异构性,构建生态承载力指标评价体系,快速、全面地进行生态承载力的评价。
一种生态承载力的遥感估算方法,包括以下步骤:
S1:构建生态承载力指标评价体系,分别从干扰强度、生态系统状态、人类响应三个方面建立对应的指标层;
S2:基础资料收集:获取位于待估算区域的卫星数据、数字地形数据、行政区划矢量图、集水单元矢量图、土壤图、土壤样品等;
S3:根据所述卫星数据和土壤图获取降水量分布图、土地利用分类图及植被覆盖度情况,构建植被覆盖度数据库;
S4:基于步骤S2-S3中所述基础数据、降水量分布图、土地利用分类图、植被覆盖度,分别以行政区划和集水单元为评价单元,进行生态承载力各单指标估算;
S5:对所述各单指标,利用层次分析法对各单指标进行权重因子确定,对待估算区域进行生态承载力综合评价。
进一步地,所述S1还基于“压力-状态-响应”(PSR)模型建立相应的生态承载力评价指标体系,确定生态承载力评价指标层。
进一步地,所述S2中的卫星数据具体包括:反射率数据、增强植被指数EVI和地表温度数据LST、降水卫星数据。
进一步地,所述所述S2中的土壤图为全国1:400万土壤类型分布图。
进一步地,所述S3还包括如下步骤:
S3.1:对所述遥感数据进行预处理;
S3.2:根据S3.1中预处理后的数据进行年降水量计算、土地覆盖精细分类及归一化植被指数NDVI计算;
S3.3:根据S3.2中得到的土地精细覆盖分类结果及归一化植被指数,结合土壤类型图,分别选取各种土地利用类型99.5%的置信区间作为该类型的NDVIV,选取各种土壤类型0.05%的置信区间作为NDVIS,获取流域NDVIV、NDVIS分布图,计算流域植被覆盖度图;
S3.4:根据S3.3得到的植被覆盖度结果,构建植被覆盖度数据库。
更进一步地,所述S3.1的预处理包括:格式转化、投影转换、数据重采样、裁剪、反射率数据降维、LST数据重建。
更进一步地,所述的LST数据重建采用光谱角距离加权算法,利用多时相数据的反射率光谱角作为加权系数进行LST数据重建;反射率数据当前时相和前后各特定个数时相的反射率数据用于计算参考像元和空值像元的光谱角距离;光谱角距离小于0.1,将参考像元放入参考像元集;通过下式对LST图像中空值像元进行重建:
其中,LSTFP是空值像元的重建结果,NRPs是参考像元集中的参考像元个数,LSTj是第j个参考像元的LST值,ωj为LSTj的加权系数,ST为相似度阈值,SADj是第j个参考像元和空值像元的光谱角距离。
进一步地,所述S4还包括如下步骤:
S4.1:基于S1中生态承载力指标评价体系及评价指标层,建立各个单指标评价模型;
S4.2:根据S4.2的单指标评价模型,以遥感像元为基本单元,计算8种单指标;
S4.3:分别以行政区划和集水单元为基本评价单元,计算8种单指标的估算量。
更进一步地,所述S4.2中8种单指标分别为人类干扰强度、灾害频率、污染排放强度、土壤侵蚀模数与等级、NPP、景观破碎度、综合弹性和成效指数。
进一步地,所述S5还包括如下步骤:
S5.1:建立生态承载力综合评价模型;
S5.2:利用层次分析法对各单指标进行权重因子确定;
S5.3:根据S5.1的生态承载力评价模型和S5.2中所述的各单指标进行权重因子,以行政区划和集水单元为基本评价单元,计算各评价单元的生态承载力。
与现有技术相比,本发明具有以下优点:
(1)本发明所述方法利用遥感数据获取方便优势,结合地面调查数据,构建一种生态承载力的遥感估算方法,能够快速、全面地实现待估算区域的生态承载力评价;
(2)本发明所述方法将行政区划单元与集水单元评价相结合,便于政府规划和管理部门准确定位生态承载力较差地区,并及时开展相应的控制和管理措施。
具体实施方式
下面结合附图和实施例,对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明,但不用来限制本发明的范围。
中分辨率成像光谱仪MODIS(Moderate Resolution ImagingSpectroradiometer)数据由于其光谱范围广,重返周期短,易获取等优点,被广泛应用与地表、生物圈、大气和海洋等的长期全球观测。MODIS数据产品算法成熟,类型多样,成为很多部门日常业务的重要数据源。本发明中所提到的用于土地覆盖精细分类的卫星优选为MODIS,卫星遥感数据尤指MODIS产品数据。热带测雨卫星TRMM(Tropical RainfallMeasuring Mission satellite)降水量估计产品的精度较高,现已有较多气象站点观测降水数据与TRMM数据的对比实验并验证了TRMM数据估算降水量的有效性,本发明中所提到的降水卫星数据尤指TRMM 3B43数据。
图1是本发明一种实施方式的生态承载力的遥感估算方法的流程图。参照图1,所述方法包括以下步骤:
S1:构建生态承载力指标评价体系,分别从干扰强度、生态系统状态、人类响应三个方面建立对应的指标层;
S2:基础资料收集:获取位于待估算区域的卫星数据、数字地形数据、行政区划矢量图、集水单元矢量图、土壤图、土壤样品等;
S3:根据所述卫星数据和土壤图获取降水量分布图、土地利用分类图及植被覆盖度情况,构建植被覆盖度数据库;
S4:基于步骤S2-S3中所述基础数据、降水量分布图、土地利用分类图、植被覆盖度,分别以行政区划和集水单元为评价单元,进行生态承载力各单指标估算;
S5:对所述各单指标,利用层次分析法对各单指标进行权重因子确定,对待估算区域进行生态承载力综合评价
所述S1具体包括:
基于“压力-状态-响应”(PSR)模型,建立相应的生态承载力评价指标体系,确定生态承载力评价指标层。
所述S2中卫星数据具体包括TRMM 3B43数据、FY-2卫星数据及4类MODIS产品数据:反射率数据、增强植被指数(EVI)、地表温度数据(LST)和净初级生产力数据(NPP)。
S3进一步包括如下步骤:
S3.1:对所述的TRMM/FY-2进行了格式转化、投影信息转换等处理。
对所述4类MODIS数据,以MODIS重投影工具MRT(MODIS Reprojection Tool)为平台,对卫星数据进行投影转换、拼接等处理,使得拼接后图像覆盖整个待估算区域。用专业的遥感处理软件ENVI对4类数据进行重采样,4类数据产品全部统一采样到500m空间分辨率,前3类数据产品采样到32天时间分辨率(一年12个时相),并计算每类产品一年的最大值、最小值和平均值3个波段,附加到每类数据末尾,其中反射率数据7个波段分别计算最大值、最小值和平均值,最终共计135个波段,作为数据集合参与地物精细分类。NPP数据则直接采用下载得到的时间分辨率为1年的合成数据。
反射率数据降维通过波段件相关性计算找到相关性高的波段组,比较其信噪比,去掉信噪比低的波段。
对有很多像元缺失的LST数据进行重建。LST数据重建采用光谱角距离加权算法,利用多时相数据的反射率光谱角作为加权系数进行LST数据重建。反射率数据当前时相和前后各特定个数时相的反射率数据用于计算参考像元和空值像元的光谱角距离。光谱角距离小于0.1,将参考像元放入参考像元集。通过下式对LST图像中空值像元进行重建:
其中,LSTFP是空值像元的重建结果,NRPs是参考像元集中的参考像元个数,LSTj是第j个参考像元的LST值,ωj为LSTj的加权系数,ST为相似度阈值,SADj是第j个参考像元和空值像元的光谱角距离。
S3.2:根据S3.1中预处理后的TRMM 3B43数据,在ENVI软件支持下,将TRMM 3B43降水速率层从HDF-EOS文件中提取出来,之后用降水速率乘以各月的总时间,生成月降水量格点数据,再进行年值统计。使用待估算区域矢量边界对其进行切割,提取出待估算区范围内的年降水量。
根据S3.1中预处理后的MODIS数据,基于中国科学院土地利用遥感监测分类系统对待估算区域进行分类,将土地利用类型分为6大类,20小类;根据S3.1预处理后的地表反射率数据进行归一化植被指数(NDVI)计算,公式为:
其中,αNIR是近红外波段的反射率,αRED是红外波段的反射率;
S3.3:根据S3.2中得到的土地精细覆盖分类结果及归一化植被指数,结合土壤类型图,分别选取各种土地利用类型99.5%的置信区间作为该类型的NDVIV,选取各种土壤类型0.05%的置信区间作为NDVIS,获取流域NDVIV、NDVIS分布图,计算流域植被覆盖度图。流域植被覆盖度Cov的计算公式如下:
S4进一步包括如下步骤:
S4.1:基于S1中生态承载力指标评价体系及评价指标层,建立各个单指标评价模型;
S4.2:根据S4.2的单指标评价模型,以遥感像元为基本单元,计算8种单指标;
S4.3:分别以行政区划和集水单元为基本评价单元,计算8种单指标的估算量。
人类干扰强度主要为人类干扰强度指数指标,建立人类干扰强度模型:
其中,PDI为人类干扰强度指数,BA为评价单元内建筑区的面积,AA为评价单元内农田的面积,A为评价单元的面积。
洪涝、干旱灾害用灾害频率指标表示,建立灾害频率模型:
其中,API为灾害频率,Y'为发生灾害的年数,Y为参与灾害统计的总年数。将将降水量小于年均降水量的25%的年份定位灾害年。
构建污染排放强度模型:
其中,PI为点源与非点源污染排放强度指数,P为评价单元内点源与非点源总污染排放,A为评价单元的面积,p为单位面积最大污染排放量。
构建土壤侵蚀强度模型:
I=R·K·LS·C·P
其中,SEM土壤侵蚀强度指数,I为评价单元的土壤侵蚀量,A为评价单元的面积,i为单位面积的最大土壤侵蚀量;R为降雨-径流侵蚀因子,K为土壤可侵蚀因子,LS为地形因子,C为覆盖-管理因子,P为水土保持措施因子。
生态活力用净初级生产力NPP来表征,构建生态活力模型:
RDI=(0.629+0.237·PER+0.00313·PER2)2
PER=(∑T/12)·58.93/r
其中T为月平均温度,PER为植被地区可能蒸散率,RDI为辐射干燥度,r为年降水量。
构建景观破碎度模型:
C=∑ni/A
其中,C为景观破碎度,∑ni为景观中所有景观类型的斑块总数,LA为景观的总面积,CI为归一化指数,Cmax为最大破碎度。
生态系统的恢复能力主要用综合弹性指标表示,构建综合弹性模型:
其中,E为综合弹性,c分别为评价单元内林地、水体、草地、农田、建设用地和未利用地的面积,A为评价单元的总面积,f分别为不同土地利用类型所占的权重,其中,
f林地、水体=1
f草地=0.8
f农田=0.6
f建筑用地=0.2
f未利用地=0.4
人类响应措施用成效指数(RI)表示,通过高分辨率遥感影像,和现场调查、文献调研,识别流域内维护生态系统的措施,比如水土保持措施、退耕还林还草措施、水资源净化设备等,结合专家打分,得到成效指数;
S5进一步包括如下步骤:
S5.1:建立生态承载力综合评价模型:
CA=∑F·f
其中,F分别为单指标评价因子,包括人类干扰强度指数(PDI)、灾害频率(API)、点源与非点源污染排放强度指数(PI)、土壤侵蚀模数(SEM)、生态活力(NPP)、景观破碎度(CI)、综合弹性(E)和成效指数(RI)指标;f分别为单指标所占权重;
S5.2:利用层次分析法对各单指标进行权重因子确定;
根据上述步骤PSR模型中确定的指标体系,将要素层其视为同层子目标,按总目标、各子目标、决策变量的顺序分解成不同的层次结构。将生成的集合按照层次结构分解后,应用专家评估法,构造两两比较的权重判断矩阵:
其中,Ann表示目标或情景之间的对比权值。
构造权重判断矩阵如下:
对构造出的权重判断矩阵,将判断矩阵:
A=(aij)n×n,n=1,2,……,m
元素按列作归一化处理,得
其中,
将矩阵的元素按行相加,得向量:
其中,
向量作归一化处理,得所求特征向量:
W=(ω1,ω2,…,ωn)Τ,
其中,
计算所述权重判断矩阵的特征向量,求出情景权值特征向量与对应的目标权值特征向量;求解所述判断矩阵的最大特征值,求出判断矩阵的最大特征值:
根据求得的最大特征值计算矩阵的一致性。只有当矩阵一致性小于0.1时,矩阵计算有效,若矩阵一致性大于或者等于0.1,则需要专家重新评估,生成新的判断矩阵并重新计算特征值来判定新矩阵是否有效。当矩阵有效时,则根据情景权值特征向量与对应的目标权值特征向量,确定出最优值,其中,
fPDI=-0.21
fAPI=-0.09
fPI=-0.2
fSEM=-0.1
fNPP=0.05
fCI=-0.05
fE=0.1
fRI=0.2
S5.3:根据S5.1的生态承载力评价模型和S5.2中所述的各单指标进行权重因子,以行政区划和集水单元为基本评价单元,计算各评价单元的生态承载力。
图2是本实施例得到的2010年海河流域生态承载力结果图(数值越大,生态承载力越好),分别以集水单元及行政区划为评价单元。图2表明,海河流域生态承载力较低的地方主要集中在城镇集中的中东和南部平原地区,北京、天津及西部地区次之,北部山区及太行山沿线相对比较好。行政区划与集水单元评价相结合便于政府规划和管理部门准确定位生态承载力情况较差地区,并及时开展相应的控制和管理措施。
以上实施方式仅用于说明本发明,而非对本发明的限制,有关技术领域的普通技术人员,在不脱离本发明的精神和范围的情况下,还可以做出各种变化和变型,因此所有等同的技术方案也属于本发明的范畴,本发明的专利保护范围应由权利要求限定。