一种基于MODIS数据的区域生态环境质量综合评价方法

技术领域

本发明涉及一种基于MODIS数据的区域生态环境质量综合评价方法，属于生态环境监测和遥感地学应用技术领域。

背景技术

气候变化和人类活动的双重影响加剧了对生态环境的干扰和破坏，导致全球生态环境问题日益突出。客观高效地监测评价区域生态环境质量，科学认识区域生态环境状况及变化趋势，对于促进人类社会与自然环境和谐共处，加强区域生态环境治理和国土空间开发利用规划，实现社会经济绿色高质量发展具有重要意义。

生态环境是指生物和非生物构成的环境的总称，生态环境质量评价是指按照特定评价标准，定性分析或定量评判生态环境质量的好坏优劣。当前，生态环境质量评价研究多采用指标体系法，基于压力-状态-响应(PSR)模型筛选构建指标体系，通过层次分析法、模糊数学法等确定权重从而评价生态环境质量状况，但这些方法易受人为主观因素的影响，数据多来源于社会经济统计方面，一定程度上限制了其使用。

遥感技术具有客观准确、范围广、精度高、周期性强等特点，在生态环境监测领域应用潜力巨大。目前，大多采用单一遥感监测指标对生态系统某一方面的特征进行评价，如利用不透水地表覆盖度评价城市土地利用，地表温度评测城市热环境等。《生态环境状况评价技术规范》中提出了主要基于遥感技术的生态环境状况指数EI，由植被覆盖指数、水网密度指数、生物丰度指数、土地退化指数和污染负荷指数构成，前3个指标可由遥感获取，但其他指标需结合地面观测和各地区年度统计数据，受地域、时间、尺度的限制较大，同时权重设置受人为因素影响较大，且EI指数只能笼统说明一个地区的生态状况，无法实现空间的可视化。徐涵秋基于Landsat遥感影像提取了绿度、热度、湿度、干度4个指标，采用主成分分析法构建了传统遥感生态指数RSEI，但由于Landsat数据源的限制，RSEI研究范围较小，适用性不高，多应用于市县、流域、矿区等小范围的生态环境质量评价。

Landsat数据时间分辨率较低，受天气、地形等条件的影响难以获取区域大尺度范围内同一时期的高质量影像。同时RSEI对区域生态环境质量的评价不够综合全面，人类社会活动如农业发展、城市扩张等都会对生态系统生产力和大气环境状况产生重要影响，且RSEI未能充分利用遥感影像的有用信息，使得RSEI在更广域尺度上的应用具有一定的局限性和不确定性。

发明内容

针对现有技术存在的不足，本发明提出了一种基于MODIS数据的区域生态环境质量综合评价方法，构建了基于MODIS的新型遥感生态指数RSEI

为解决上述问题，本发明所采取的技术方案是：

一种基于MODIS数据的区域生态环境质量综合评价方法，包括以下步骤：

S1，基于获取的待研究区任意年份区间的MODIS数据进行预处理并根据IGBP分类体系提取水体得到水体信息二值图；

S2，基于预处理后的MODIS数据提取生态因子，然后将生态因子进行标准化处理，使值域范围统一在0～1内，其中生态因子包括绿度因子、湿度因子、热度因子、干度因子、灰度因子和生产力因子；

S3，基于标准化处理后的生态因子进行主成分分析，以累计方差贡献率大于90％为阈值选取主成分波段，依照所选各波段的方差贡献率对相应主成分进行加权求和，归一化后得到基于MODIS的遥感生态指数RSEI

S4，基于遥感生态指数RSEI

进一步地，所述步骤S1中的MODIS数据包括时间序列的地表反射率数据、地表温度数据、总初级生产力数据和年度土地覆被数据；

其中，地表反射率数据、地表温度数据和总初级生产力数据用于计算生态因子，年度土地覆被数据用于水体信息掩膜。

进一步地，所述步骤S2中提取绿度因子包括以下步骤：

S2-11，根据步骤S1中所得的地表反射率数据计算时间序列的植被指数NDVI，得到原始时序NDVI数据；

S2-12针对步骤S2-11所得的原始时序NDVI数据，根据数据质量控制波段检测云和阴影等噪声，剔除异常值，利用线性插值对缺失值进行补齐，得到插值时序NDVI数据；

S2-13，针对步骤S2-12所得的插值时序NDVI数据，采用S-G滤波进行重构处理，得到高质量的重构时序NDVI数据，对每年7-9月的重构时序NDVI数据取平均值，得到NDVI均值影像；

S2-14，针对步骤S2-13所得的NDVI均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区绿度因子专题图。

进一步地，所述步骤S2中提取湿度因子包括以下步骤：

S2-21，根据步骤S1所得的地表反射率数据计算时间序列的地表含水量指数SWCI，得到原始时序SWCI数据；

S2-22，针对步骤S2-21所得的原始时序SWCI数据，根据产品质量控制波段检测云和阴影等噪声，剔除异常值，利用线性插值对缺失值进行补齐，得到插值时序SWCI数据；

S2-23，针对步骤S2-22所得的插值时序SWCI数据，采用S-G滤波进行重构处理，得到高质量的重构时序SWCI数据，对每年7-9月的重构时序SWCI数据取平均值，得到SWCI均值影像；

S2-24，针对步骤S2-23所得的SWCI均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区湿度因子专题图。

进一步地，所述步骤S2中提取热度因子包括以下步骤：

S2-31，根据步骤S1所得的地表温度数据计算时间序列的地表温度LST，得到原始时序LST数据；

S2-32，针对步骤S2-31所得的原始时序LST数据，根据产品质量控制波段剔除低质量无效观测及异常值，得到高质量时序LST数据，对每年7-9月的高质量时序LST数据取平均值，并得出LST均值影像；

S2-33，针对步骤S2-32所得的LST均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影以及重采样，再利用研究区矢量文件对影像进行裁剪，得到研究区热度因子专题图。

进一步地，所述步骤S2中提取干度因子包括以下步骤：

S2-41，根据步骤S1所得的地表反射率数据计算时间序列的建筑-裸土指数NDSI

S2-42，针对步骤S2-41所得的原始时序NDSI

S2-43，针对步骤S2-42所得的插值时序NDSI

S2-44，针对步骤S2-43所得的NDSI

进一步地，所述步骤S2中提取灰度因子包括以下步骤：

S2-51，根据步骤S1所得的地表反射率数据计算时间序列的差值指数DI，得到原始时序DI数据；

S2-52，针对步骤S2-51所得的原始时序DI数据，根据产品质量控制波段检测云和阴影等噪声，剔除异常值，利用线性插值对缺失值进行补齐，得到插值时序DI数据；

S2-53，针对步骤S2-52所得的插值时序DI数据，采用S-G滤波进行重构处理，得到高质量的重构时序DI数据，对每年7-9月的重构时序DI数据取平均值，并得出DI均值影像；

S2-54，针对步骤S2-53所得的DI均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区灰度因子专题图。

进一步地，所述步骤S2中提取生产力因子包括以下步骤：

S2-61，根据步骤S1所得的总初级生产力数据计算时间序列的GPP，得到原始时序GPP数据；

S2-62，针对步骤S2-61所得的原始时序GPP数据，根据产品质量控制波段剔除低质量无效观测及异常值，利用线性插值对缺失值进行补齐，得到插值时序GPP数据；

S2-63，针对步骤S2-62所得的插值时序GPP数据，采用S-G滤波进行重构处理，得到高质量的重构时序GPP数据，对每年7-9月的重构时序GPP数据取平均值；

S2-64，针对步骤S2-63所得的GPP均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区生产力因子专题图。

进一步地，所述步骤S2-11中绿度因子为归一化差值植被指数NDVI，其计算公式为：

其中，ρ

进一步地，所述步骤S2-21中湿度因子为地表含水量指数SWCI的计算公式为：

其中，ρ

进一步地，所述步骤S2-31中热度因子为地表温度LST的计算公式为：

LST＝0.02*DN-273.15；

其中，DN表示地表温度数据影像的灰度值。

进一步地，所述步骤S2-41中干度因子为建筑-裸土指数NDSI

其中，NDBI

进一步地，所述步骤S2-51中灰度因子为差值指数DI的计算公式为：

DI＝ρ

其中，ρ

进一步地，所述步骤S2-61中生产力因子为地表总初级生产力GPP的计算公式为：

GPP＝0.1*DN；

其中，DN表示总初级生产力数据影像的灰度值。

进一步地，所述线性插值是指对云或阴影覆盖区域的植被指数，利用前后晴空时相的植被指数值进行替补，线性插值为运用一次多项式插值方式补齐时序数据的缺失，公式如下：

其中，

进一步地，所述S-G滤波是一种重建时间序列植被指数的迭代平滑方法，S-G滤波实质上是一种移动窗口的加权平均算法，加权系数通过在滑动窗口内对于给定高阶多项式的最小二乘拟合得出，其过程可用下式表示：

其中，VI是原始植被指数，VI

进一步地，所述步骤S2中生态因子标准化计算公式为：

其中，Y

进一步地，所述步骤S3中RSEI

进一步地，所述步骤S3中基于MODIS的遥感生态指数RSEI

PC

其中，f(·)表示主成分运算，PC

进一步地，所述步骤S4中的分级赋值是指将特定年份的RSEI

所述差值法是指将最近年份和起始年份的RSEI

本发明的有益效果在于：

1)本发明所提出的区域生态环境质量综合评价方法，构建了一种基于MODIS的新型遥感生态指数RSEI

2)本发明所提出的方法构建了一种基于MODIS的新型遥感生态指数RSEI

3)本发明所提出的方法确定指标权重时，相比传统RSEI，不仅仅利用第一主成分，而是考虑主成分累计方差贡献率，以累计方差贡献率达到合理阈值选取主成分，充分利用了遥感影像的有用信息，避免了人为主观因素的影响。评价指标计算过程清晰，赋权科学合理，评价流程科学，提高了区域生态质量评价的准确性和客观性。

附图说明

下面结合附图对本发明作进一步的说明。

图1为本发明的流程示意图；

图2为时序NDVI重构前后的对比图；

图3是太行山-燕山东南部平原及沿海地带的区域2001年的生态环境质量等级空间分布图；

图4是太行山-燕山东南部平原及沿海地带的区域2010年的生态环境质量等级空间分布图；

图5是太行山-燕山东南部平原及沿海地带的区域2019年的生态环境质量等级空间分布图；

图6是太行山-燕山东南部平原及沿海地带的区域2001-2010年的生态环境质量变化空间分布图；

图7是太行山-燕山东南部平原及沿海地带的区域2010-2019年的生态环境质量变化空间分布图；

图8是太行山-燕山东南部平原及沿海地带的区域2001-2019年的生态环境质量变化空间分布图；

图9为实例局部结果与Landsat高精度影像的对比结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本实施例以太行山-燕山东南部平原及沿海地带为研究区域，监测时段为2001-2019年进行说明。

参照图1，图1为本发明一种基于MODIS数据的区域生态环境质量综合评价方法的流程示意图，包括以下步骤：

S1，基于获取的待研究区任意年份区间的MODIS数据进行预处理并根据IGBP分类体系提取水体得到水体信息二值图；

具体的，下载研究区范围2001年、2010年、2019年的MODIS数据，MODIS数据包括时间序列的地表反射率数据MOD09A1、地表温度数据MOD11A2、总初级生产力数据MOD17A2H和年度土地覆被数据MCD12Q1，地表反射率数据MOD09A1、地表温度数据MOD11A2、总初级生产力数据MOD17A2H用于计算生态指标因子，时间分辨率均为8天，空间分辨率MOD11A2为1000米，其余均为500米；年度土地覆被数据MCD12Q1用于水体信息掩膜，根据IGBP分类体系提取水体，将水体像元赋值为0，其他像元赋值为1，得到水体信息二值图。

S2，基于预处理后的MODIS数据提取生态因子，然后将生态因子进行标准化处理，使值域范围统一在0～1内，其中生态因子包括绿度因子、湿度因子、热度因子、干度因子、灰度因子和生产力因子；

具体的，提取绿度因子、湿度因子、热度因子、干度因子、灰度因子和生产力因子的步骤如下：

提取绿度因子包括以下步骤：

S2-11，根据步骤S1中所得的地表反射率数据计算时间序列的植被指数NDVI，得到原始时序NDVI数据；

NDVI与植被生物量、叶面积指数、植被覆盖度等密切相关，是表征地表植被覆盖与生长状况的有效指标，计算公式为：

其中，ρ

植被指数基于遥感光谱信息计算得到，但遥感光学数据获取过程中，受云层、天气、传感器等因素影响，得到的时序植被指数存在大量噪声，因此需要对时序数据去除噪声，重构后得到高质量的时序数据集。NDVI、SWCI、NDSI

S2-12，针对步骤S2-11所得的原始时序NDVI数据，根据产品质量控制波段检测云和阴影等噪声，剔除超过NDVI取值范围(-1,1)的异常值，利用线性插值对缺失值进行补齐；基于MOD09A1地表反射率计算得到的原始时序NDVI，根据产品自带的质量标识识别云和阴影等噪声，先采用线性插值方式，利用前后晴空时相的植被指数值对噪声和缺失进行替补，得到插值时序NDVI数据；

S2-13，针对步骤S2-12所得的插值时序NDVI数据，采用S-G滤波进行重构处理，得到高质量的重构时序NDVI数据，对每年7-9月的重构时序NDVI数据取平均值，并得到NDVI均值影像；

S2-14，针对步骤S2-13所得的NDVI均值影像，根据步骤S1中的水体信息二值图进行掩膜，拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区绿度因子专题图；

提取湿度因子包括以下步骤：

S2-21，根据步骤S1所得的地表反射率数据计算时间序列的地表含水量指数SWCI，得到原始时序SWCI数据；

SWCI能有效提取植被冠层及地表的水分含量，计算公式为：

其中，ρ

S2-22，针对步骤S2-21所得的原始时序SWCI数据，根据产品质量控制波段检测云和阴影等噪声，剔除超过SWCI取值范围(-1,1)的异常值，利用线性插值对缺失值进行补齐，得到插值时序SWCI数据；

S2-23，针对步骤S2-22所得的插值时序SWCI数据，采用S-G滤波进行重构处理，得到高质量的重构时序SWCI数据，对每年7-9月的重构时序SWCI数据取平均值，得到SWCI均值影像；

S2-24，针对步骤S2-23所得的SWCI均值影像，根据步骤S1中的水体信息二值图进行掩膜，拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区湿度因子专题图。

提取热度因子包括以下步骤：

S2-31，根据步骤S1所得的地表温度数据计算时间序列的地表温度LST，得到原始时序LST数据；

LST由MOD11A2地表温度数据直接计算得到，单位为℃，计算公式为：

LST＝0.02*DN-273.15；

其中，DN表示地表温度数据MOD11A2影像的灰度值。MOD11A2地表温度数据自带质量标识波段，用质量标识去除低质量无效观测，能有效降低噪声影响，得到高质量时序LST数据。

S2-32，针对步骤S2-31所得的原始时序LST数据，根据产品质量控制波段剔除低质量无效观测及异常值，得到高质量时序LST数据，对每年7-9月的高质量时序LST数据取平均值，并得出LST均值影像；

S2-33，针对步骤S2-32所得的LST均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影以及重采样，再利用研究区矢量文件对影像进行裁剪，得到研究区热度因子专题图。

提取干度因子包括以下步骤：

S2-41，根据步骤S1所得的地表反射率数据计算时间序列的建筑-裸土指数NDSI

NDSI

其中，NDBI

S2-42，针对步骤S2-41所得的原始时序NDSI

S2-43，针对步骤S2-42所得的插值时序NDSI

S2-44，针对步骤S2-43所得的NDSI

提取灰度因子包括以下步骤：

S2-51，根据步骤S1所得的地表反射率数据计算时间序列的差值指数DI，得到原始时序DI数据；

DI能反映大气中颗粒物PM2.5浓度，表征大气环境质量状况，计算公式为：

DI＝ρ

其中，ρ

S2-52，针对步骤S2-51所得的原始时序DI数据，根据产品质量控制波段检测云和阴影噪声，剔除异常值，利用线性插值对缺失值进行补齐，得到插值时序DI数据；

S2-53，针对步骤S2-52所得的插值时序DI数据，采用S-G滤波进行重构处理，得到高质量的重构时序DI数据，对每年7-9月的高质量时序DI数据取平均值，并得出DI均值影像；

S2-54，针对步骤S2-53所得的DI均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区灰度因子专题图。

提取生产力因子包括以下步骤：

S2-61，根据步骤S1所得的总初级生产力数据计算时间序列的GPP，并得出原始时序GPP数据；

GPP表征生态系统总初级生产力，单位为gC/m

GPP＝0.1*DN；

其中，DN表示总初级生产力数据MOD17A2H影像的灰度值。

S2-62，针对步骤S2-61所得的原始时序GPP数据，根据产品质量控制波段剔除低质量无效观测及异常值，利用线性插值对缺失值进行补齐，得到插值时序GPP数据；

S2-63，针对步骤S2-62所得的插值时序GPP数据，采用S-G滤波进行重构处理，得到高质量的重构时序GPP数据，对每年7-9月的重构时序GPP数据取平均值；

S2-64，针对步骤S2-63所得的GPP均值影像，根据步骤S1中的水体信息二值图进行掩膜，并拼接转投影，再利用研究区矢量文件对影像进行裁剪，得到研究区生产力因子专题图。

上述提取绿度因子、湿度因子、热度因子、干度因子、灰度因子和生产力因子中，线性插值是指对云或阴影覆盖区域的植被指数，利用前后晴空时相的植被指数值进行替补，线性插值为运用一次多项式插值方式补齐时序数据的缺失，线性插值的计算公式为：

其中，

上述提取绿度因子、湿度因子、热度因子、干度因子、灰度因子和生产力因子中，其中，S-G滤波是一种常用的重建时间序列植被指数的迭代平滑方法，能够清晰地描述时间序列的季节性变化，不受数据时间空间尺度和传感器的限制，S-G滤波计算过程可用下式表示：

其中，VI是原始植被指数，VI

S2-7，基于上述所得的各项生态因子进行标准化处理，消除量纲，使值域范围统一在(0,1)内，标准化计算公式为：

其中，Y

S3，基于标准化处理后的生态因子进行主成分分析，以累计方差贡献率大于90％为阈值选取主成分波段，依照所选各波段的方差贡献率对相应主成分进行加权求，归一化后得到基于MODIS的遥感生态指数RSEI

具体的，RSEI

针对步骤S2所得的各生态因子标准化数据，利用ENVI软件进行主成分分析(PCA)，以主成分累计方差贡献率大于90％选取主成分波段，依照方差贡献率对相应主成分进行加权求和，归一化后得到基于MODIS的遥感生态指数RSEI

PC

式中，f(·)表示主成分运算，PC

S4，基于遥感生态指数RSEI

具体的，对研究区2001年、2010年、2019年的的遥感生态指数RSEI

尽管这里参照本发明的多个解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。更具体地说，在本申请公开、附图和权利要求的范围内，可以对主题组合布局的组成部件和/或布局进行多种变型和改进。除了对组成部件和/或布局进行的变形和改进外，对于本领域技术人员来说，其他的用途也将是明显的。