[发明专利]一种基于MODIS数据的区域生态环境质量综合评价方法

权利要求书

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

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

MODIS数据包括时间序列的地表反射率数据、地表温度数据、总初级生产力数据和年度土地覆被数据；

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

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

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

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-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-31，根据步骤S1所得的地表温度数据计算时间序列的地表温度LST，得到原始时序LST数据；

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

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

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

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

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

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

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

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

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

灰度因子为差值指数DI的计算公式为：

DI＝ρ_red-ρ_nir；

其中，ρ_red、ρ_nir为MODIS数据红波段及近红外波段的反射率值，分别对应地表反射率数据的第1波段和第2波段；

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的计算公式为：

GPP＝0.1*DN；

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

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

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

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

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

RSEI_M为基于MODIS的新型遥感生态指数，通过集成表征生态系统绿度、湿度、热度、干度、灰度、生产力的NDVI数据、SWCI数据、LST数据、NDSI_M数据、DI数据和GPP数据六个指标，基于主成分分析结果，以累计方差贡献率超过一定阈值选取主成分，加权计算得到；

基于MODIS的遥感生态指数RSEI_M的计算过程为：

PC_i＝f(NDVI,SWCI,LST,NDSI_M,DI,GPP)，

其中，f(·)表示主成分运算，PC_i表示获取的主成分分量，NDVI代表绿度因子,SWCI代表湿度因子,LST代表热度因子,NDSI_M代表干度因子,DI代表灰度因子,GPP代表生产力因子，i表示第i个主成分波段，e_i为第i个主成分波段的方差贡献率，k为选取的主成分波段个数,RSEI₀为初始的遥感生态指数值，RSEI_0min为初始遥感生态指数RSEI₀的最小值，RSEI_0max为初始遥感生态指数RSEI₀的最大值，RSEI_M为基于MODIS的遥感生态指数，RSEI_M值越接近1，表示生态环境质量越好，越接近0，表示生态环境质量越差；

S4，基于遥感生态指数RSEI_M数据，按照分级赋值方法对研究区的生态环境质量做出综合评价，并利用差值法得出年份区间的生态环境质量变化。