CN115795237B - 一种基于四元Copula的综合遥感生态指数并行计算方法 - Google Patents

一种基于四元Copula的综合遥感生态指数并行计算方法 Download PDF

Info

Publication number
CN115795237B
CN115795237B CN202211514932.9A CN202211514932A CN115795237B CN 115795237 B CN115795237 B CN 115795237B CN 202211514932 A CN202211514932 A CN 202211514932A CN 115795237 B CN115795237 B CN 115795237B
Authority
CN
China
Prior art keywords
index
remote sensing
ecological
copula
quaternary
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211514932.9A
Other languages
English (en)
Other versions
CN115795237A (zh
Inventor
杨海波
李其钊
王宗敏
梁黄陶钧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhengzhou University
Original Assignee
Zhengzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhengzhou University filed Critical Zhengzhou University
Priority to CN202211514932.9A priority Critical patent/CN115795237B/zh
Publication of CN115795237A publication Critical patent/CN115795237A/zh
Application granted granted Critical
Publication of CN115795237B publication Critical patent/CN115795237B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A30/00Adapting or protecting infrastructure or their operation
    • Y02A30/60Planning or developing urban green infrastructure

Landscapes

  • Testing Or Calibration Of Command Recording Devices (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及遥感技术领域,应用遥感生态指标反映区域生态环境状态。选取绿度、湿度、干度以及热度四个遥感生态指标,使用多线程并行计算以及基于多核CPU的并行计算相结合的方式,创新性地使用四元Copula函数作为生态指标的计算工具,利用Gringorten公式、Kendall相关系数和最大似然估计法构建出四元Copula函数模型,并最终计算出基于Copula的综合遥感生态指数(Copula‑base remote sensing ecological index,CRSEI)。本发明构建的综合遥感生态指数,能够最大程度囊括四个指标所反映的地表类型状况信息以及四个指标间的相关性信息,更加综合、全面地通过栅格可视化结果从时间以及空间上反映生态环境质量,为生态环境的监测、评价以及治理提供可靠、有力的技术支持。

Description

一种基于四元Copula的综合遥感生态指数并行计算方法
技术领域
本发明涉及遥感技术领域,更具体涉及一种基于四元Copula的综合遥感生态指数并行计算方法。
背景技术
综合遥感生态指数是反映区域生态环境质量状况的一系列遥感生态指标的综合。目前常用的生态指标包括绿度(NDVI)、湿度(WET)、干度(NDBSI)、热度(LST)等,单一的遥感生态指标评价环境的方法只能监测某一种环境特征,如用植被指数监测森林覆盖;用地表温度监测城市热岛效应;用建筑指数监测城市扩张速率等,但不能对复杂的生态环境进行综合全面的评价。同时常用的生态评价方法基于权重组合的生态环境状况指数EI和基于主成分分析法(PCA)计算的遥感生态指数(Remote sensing ecological index,RSEI);其中EI不能可视化的反映区域内生态环境状况的分布情况,RSEI仅使用第一主要成分计算遥感生态指数,第一主成分的贡献率有时会比较低,无法保证生态指标的大部分信息被合理利用。
本发明基于Copula函数对遥感栅格数据使用并行计算的方式计算出基于Copula的综合遥感生态指数(Copula-base remote sensing ecological index,CRSEI),能够高效计算Copula函数模型并且将结果可视化,同时Copula函数囊括遥感生态指标的更多信息外,还能够考虑指标变量的相关性,从时间变化和空间演变上更合理准确地反映生态环境质量,为生态环境的监测、评价以及治理提供可靠、有力的技术支持。
发明内容
为了更好的解决上述问题,本发明提供一种基于四元Copula的综合遥感生态指数并行计算方法,包括,
S1.收集多研究区域十年的光谱遥感栅格影像数据、绿度指标月合成数据、热度指标月合成数据及年尺度的土地利用数据;用所述多光谱遥感栅格影像数据分别计算月尺度下的湿度指标和干度指标;并对四个遥感生态指标数据进行归一化处理,所述四个遥感生态指标包括绿度指标、湿度指标、干度指标和热度指标;
S2.将四个遥感生态指标的栅格影像数据以像元为单位逐个读取成四个遥感生态指标时间序列,构成十年120个月尺度的四个遥感生态指标时间序列,所述四个遥感生态指标时间序列包括绿度指标时间序列、湿度指标时间序列、干度指标时间序列及热度指标时间序列;
S3.对每个像元的四个遥感生态指标时间序列分别计算四个遥感生态指标边缘分布函数,所述四个遥感生态指标边缘分布函数包括绿度指标边缘分布函数、湿度指标边缘分布函数、干度指标边缘分布函数和热度指标边缘分布函数;
S4.利用每个像元的四个遥感生态指标边缘分布函数确定四元Copula函数的模型参数,构建月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数模型,将所述四个指标时间序列重新代入确定参数后的四元Copula函数模型,计算出综合遥感生态指数;
S5.将所述综合遥感生态指数结果重新生成可视化栅格影像数据,并将所述可视化栅格影像数据归一化;
S6.设定综合遥感生态指数的阈值,将所述综合遥感生态指数根据所述阈值划分成不同生态环境质量等级。
作为本发明的一种优选技术方案,所述S1中,可以方便地获取多光谱遥感影像数据用于指标的计算,所述遥感影像数据包括含红光、绿光、蓝光、第一近红外光、第二近红外光、第一短波红外、第二短波红外的7波段的8天合成数据。
作为本发明的一种优选技术方案,所述S1中,采用多线程并行的方式高效地计算遥感生态指标并归一化,所述月尺度下的湿度指标和所述月尺度下的干度指标采用python开启线程池多线程并行计算方式如下:
S11.根据所述遥感数据计算湿度指标,计算公式为
WETNESS=0.1147×Red+0.2408×NIR1+0.3283×Blue+0.3132×Green-0.3122×NIR2-0.6416×SWIR1-0.5087×SWIR2
其中,SWIR1、SWIR2、NIR1、NIR2、Red、Blue、Green分别为第一短波红外、短二短波红外、第一近红外光、第二近红外光、红光、蓝光、绿光波段的反射率;
S12.根据所述遥感影像数据计算干度指标,计算公式为:
BSI=[(SWIR1+Red)-(NIR1+Blue)]/[(SWIR1+Red)+(NIR1+Blue)]
NDBSI=(BSI+IBI)/2
其中,SWIR1、NIR1、Red、Green、Blue分别为第一短波红外、第一近红外光、红光、绿光、蓝光波段的反射率;
S13.将上述得到的湿度指标和干度指标去除异常值,并利用最大值合成法将数据合成月尺度下的湿度指标和月尺度下的干度指标;
S14.将所述月尺度下的湿度指标和月尺度下的干度指标进行归一化处理,归一化处理计算公式为:
NIi=(Ii-Imin)/(Imax-Imin)
其中,NIi为某一正规化后的指标值,Ii为该指标在像元i处的值;Imax为该指标的最大值;Imin为该指标的最小值。
作为本发明的一种优选技术方案,所述S2中,对四元Copula函数实现基于遥感栅格数据的读取计算,使用Matlab将归一化后的四个遥感生态指标的栅格影像数据读取为矩阵形式,按照同一位置的像元点逐个读取十年120个月的月尺度的四个遥感生态指标时间序列。
作为本发明的一种优选技术方案,所述S3中,按照逐像元的方式对四个遥感生态指标时间序列计算边缘分布函数,使用Gringorten公式,逐个像元点分别计算所述四个遥感生态指标边缘分布函数,当时间序列共有N项,所述时间序列x1≤x2≤…≤xi≤…≤xN中小于等于xi的出现次数为i时,其公式如下:
作为本发明的一种优选技术方案,所述S4中,使用四元Copula函数作为边缘分布函数的连接工具,栅格影像中每个像元采用并行的方式计算出各月月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数参数,完成Copula函数模型的构建;所述四元Copula函数选择四元t-Copula函数,将所述月尺度下的绿度指标、湿度指标、干度指标及热度指标代入上述确定参数后的四元t-Copula函数,计算出对应的综合遥感生态指数。
作为本发明的一种优选技术方案,根据Sklar定理,逐像元地计算四个遥感生态指标时间序列的Copula函数模型,所述四个遥感生态指标时间序列的联合分布函数P计算公式如下:
P(x1≤X1,x2≤X2,x3≤X3,x4≤X4)=C[U(X1),V(X2),G(X3),W(X4)]=p
其中,x1为绿度指标时间序列对应变量、x2为湿度指标时间序列对应变量、x3为干度指标时间序列对应变量、x4为热度指标时间序列对应变量;U(x1)为绿度指标的边缘分布函数、V(x2)为湿度指标的边缘分布函数、G(x3)为干度指标的边缘分布函数、W(x4)为热度指标的边缘分布函数;P为绿度指标、湿度指标、干度指标及热度指标的联合分布函数;C为四元t-Copula函数,表达式为:
密度函数的表达式为:
其中,ρ为对角线上的元素全为1的4阶对称正定矩阵,|ρ|表示方阵ρ的行列式,u1为绿度指标边缘分布函数U(x1),u2为湿度指标边缘分布函数V(x2),u3为干度指标边缘分布函数G(x3),u4为热度指标边缘分布函数W(x4);tρ,k表示相关系数矩阵为ρ、自由度为k的标准四元t分布的分布函数,表示自由度为k的标准一元t分布的分布函数的逆函数;/>
作为本发明的一种优选技术方案,为了确定四元t-Copula函数模型的参数,根据Kendall相关系数估计相关系数、最大似然函数法估计自由度,算法过程如下:
S41.对四个遥感生态指标运用经验分布或边缘分布函数进行概率转化,所述边缘分布函数包括绿度指标边缘分布函数U(X1)、湿度指标边缘分布函数V(X2)、干度指数边缘分布函数G(X3)和热度指标边缘分布函数W(X4);
S42.根据利用Kendall相关系数的τ来估计四个遥感生态指标的相关系数矩阵/>
S43.运用极大似然估计方法,估计自由度参数k,计算公式如下:
作为本发明的一种优选技术方案,所述S3、S4、S5中,遥感栅格影像有数以万计的像元点,每个像元点都要计算一个与之对应的四元t-Copula函数模型,为了实现多个像元的并行计算,提升计算效率,使用Matlab的主从结构的分布式计算方式实现基于多核CPU的并行计算。初始化Matlab并行计算环境时,最初的Matlab进程自动成为主节点,同时初始化多个Matlab计算子节点,所述子节点的个数可根据计算机CPU核心数量设定,利用Parfor将子节点调用用于计算。Parfor运行之初,主节点会将Parfor循环程序之外变量传递给计算子节点。子节点运算过程时互不干扰,运算完毕,将各子节点得到的结果组合到同一个矩阵变量中,并返回到Matlab主节点,最终计算完毕则关闭计算子节点。
作为本发明的一种优选技术方案,所述S5中,将综合遥感生态指数值根据每个像元点的位置、时间信息重新生成可视化栅格数据,并将综合遥感生态指数值归一化;根据综合遥感生态指数划分生态质量等级,将综合遥感生态指数从0到1,以0.2为间隔分成5个等级,所述综合生态质量等级根据综合遥感生态指数由低到高分别是差、较差、中等、良和优5个等级,并根据综合遥感生态指数等级分析综合遥感生态指数时间变化和空间分布规律。
有益效果
1、本发明的技术方案生成的基于Copula的综合遥感生态指数(Copula-baseremote sensing ecological index,CRSEI),创新性地采用四元Copula函数作为四个遥感生态指标地连接工具计算出基于Copula的综合遥感生态指数,能够最大程度反映绿度指标、湿度指标、干度指标及热度指标四个指标间的相关信息,基于Copula的综合遥感生态指数实现了对栅格影像数据的Copula拟合计算,根据每个像元的位置、时间信息生成可视化栅格影像数据,从时间以及空间上反映生态环境质量,为生态环境的监测、评价以及治理提供可靠、有力的技术支持。
2、本发明的技术方案通过绿度指标、湿度指标、干度指标及热度指标的栅格影像数据组成的时间序列构建四元Copula函数模型,采用主从结构的分布式计算方式实现基于多核CPU的并行计算从而更高效地计算出基于Copula的综合遥感生态指数(Copula-baseremote sensing ecological index,CRSEI)并生成栅格化可视数据,更加准确、综合全面的反映区域内生态环境质量,实现对复杂的生态环境进行监测。
3、本发明通过设定基于Copula的综合遥感生态指数的阈值,将基于Copula的综合遥感生态指数划分不同生态环境质量等级,可以直接通过生态环境质量等级直观反映出区域生态环境的质量。
附图说明
图1为本发明的流程示意图;
图2为四元t-Copula计算流程图;
图3为本发明月尺度下2001年10月的综合遥感生态指数CRSEI与遥感生态指数RSEI空间分布对比图;
图4为本发明月尺度下2001年1月、4月、7月以及10月的CRSEI指数分级图;
图5为月尺度下2001-2010年的NDVI、WET、NDBSI、LST、CRSEI、RSEI的影像均值变化对比图;
图6为月尺度下2001-2010年的NDVI、WET、NDBSI、LST、CRSEI、RSEI的影像抽样点指数值对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
本发明提供一种基于四元Copula的综合遥感生态指数并行计算方法,具体的,本实施例研究区域选用河南省的一小块区域,遥感数据采用MODIS数据产品编号是MOD09A1数据、中国500mNDVI月合成数据、以及中国1kmLST月合成数据,MOD09A1数据从美国国家航空航天局(National Aeronautics and Space Administration)戈达德航天中心的数据网站获取,而NDVI和LST数据则通过地理空间数据云获取下载,使用三次卷积方法重采样方法将数据分辨率统一为500m,经过最大值合成计算调整时间分辨率为一个月,选取的时段为2001-2010年120个月。数据结合GIS、ENVI软件和编程软件进行各个步骤的处理。
如图1所示,所述基于Copula的综合遥感生态指数并行计算方法,包括S1.收集研究区域十年的多光谱遥感栅格影像数据、绿度指标月合成数据、热度指标月合成数据及年尺度的土地利用数据;用所述多光谱遥感栅格影像数据分别计算月尺度下的湿度指标和干度指标;并对四个遥感生态指标数据进行归一化处理,所述四个遥感生态指标包括绿度指标、湿度指标、干度指标数据和热度指标;S2.将四个遥感生态指标的栅格影像数据以像元为单位逐个读取成四个遥感生态指标时间序列,构成十年120个月尺度的四个遥感生态指标时间序列,所述四个遥感生态指标时间序列包括绿度指标时间序列、湿度指标时间序列、干度指标时间序列及热度指标时间序列;S3.对每个像元的四个遥感生态指标时间序列分别计算四个遥感生态指标边缘分布函数,所述四个遥感生态指标边缘分布函数包括绿度指标边缘分布函数、湿度指标边缘分布函数、干度指标边缘分布函数和热度指标边缘分布函数;
S4.利用每个像元的四个遥感生态指标边缘分布函数确定四元Copula函数的模型参数,构建月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数模型,将所述四个指标时间序列重新代入确定参数后的四元Copula函数模型,计算出综合遥感生态指数;
S5.将所述综合遥感生态指数结果重新生成可视化栅格影像数据,并将所述可视化栅格影像数据归一化;S6.设定综合遥感生态指数的阈值,将所述综合遥感生态指数根据所述阈值划分成不同生态环境质量等级。
具体的,本实施例收集研究区域内2001-2010年十年的多光谱遥感影像栅格数据、绿度指标月合成数据、热度指标月合成数据及年尺度的土地利用数据,用多光谱遥感栅格影像数据分别计算月尺度下的湿度指标和干度指标,并对四个遥感生态指标数据进行归一化处理;将四个遥感生态指标的栅格影像数据以像元为单位逐个读取成时间序列数据,构成十年120个月尺度的四个遥感生态指标时间序列数据;对每个像元的四个遥感生态指标时间序列分别计算其边缘分布函数;利用每个像元的四个遥感生态指标的边缘分布函数确定四元Copula函数的模型参数,构建月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数模型,将所述四个遥感生态指标时间序列数据重新代入确定参数后的四元Copula函数模型,计算出综合遥感生态指数;将综合遥感生态指数生成可视化栅格影像数据,更加准确、综合全面的反映区域内生态环境质量,而且能从时间以及空间上反映生态环境质量,同时还可以对生态环境比较复杂的区域进行监测评估。
进一步地,所述S1中,能够方便地获取多光谱遥感影像数据用于指标的计算,所述遥感影像数据包括含红光、绿光、蓝光、第一近红外光、第二近红外光、第一短波红外光、第二短波红外光的7波段的8天合成数。
进一步地,所述S1中,采用多线程并行的方式高效地计算遥感生态指标并归一化,所述月尺度下的湿度指标和所述月尺度下的干度指标采用python开启线程池多线程并行计算方式如下:
S11.根据所述遥感数据计算湿度指标,计算公式为
WETNESS=0.1147×Red+0.2408×NIR1+0.3283×Blue+0.3132×Green-0.3122
×NIR2-0.6416×SWIR1-0.5087×SWIR2
其中,SWIR1、SWIR2、NIR1、NIR2、Red、Blue、Green分别为第一短波红外、短二短波红外、第一近红外光、第二近红外光、红光、蓝光、绿光波段的反射率;
S12.根据所述遥感影像数据计算干度指标,计算公式为:
BSI=[(SWIR1+Red)-(NIR1+Blue)]/[(SWIR1+Red)+(NIR1+Blue)]
NDBSI=(BSI+IBI)/2
其中,SWIR1、NIR1、Red、Green、Blue分别为第一短波红外、第一近红外光、红光、绿光、蓝光波段的反射率;
S13.将上述得到的湿度指标和干度指标去除异常值,并利用最大值合成法将数据合成为月尺度下的湿度指标和月尺度下的干度指标;
S14.将所述月尺度下四个遥感生态指标进行归一化处理,归一化处理计算公式为:
NIi=(Ii-Imin)/(Imax-Imin)
其中,NIi为某一正规化后的指标值,Ii为该指标在像元i处的值;Imax为该指标的最大值;Imin为该指标的最小值。
进一步地,所述S2中,对四元Copula函数实现基于遥感栅格数据的读取计算,使用Matlab将归一化后的四个遥感生态指标的栅格影像数据读取为矩阵形式,按照同一位置的像元点逐个读取十年120个月的月尺度的四个遥感生态指标时间序列。
具体的,时间序列数据的读取过程为:从第一个位置像元点N1可以按照时间顺序读取为绿度指标时间序列数据湿度指标时间序列数据干度指标时间序列/>热度指标时间序列依次读取完所有像元点。其中如图2所示,F1(x)为热度指标时间序列F2(x)为干度指标时间序列/>F3(x)为绿度指标时间序列/>F4(x)为湿度指标时间序列/>
进一步地,如图2所示,所述S3中,按照逐像元的方式对四个遥感生态指标时间序列计算边缘分布函数,使用Gringorten公式,逐个像元点分别计算所述四个遥感生态指标边缘分布函数,当时间序列共有N项,所述时间序列x1≤x2≤…≤xi≤…≤xN中小于等于xi的出现次数为i时,其公式如下:
具体的,所述序列边缘分布函数为含有120项的序列,所述120项的序列来自十年12个月的生态指标数据分别计算得到。
进一步地,所述S4中,使用四元Copula函数作为边缘分布函数的连接工具,栅格影像中每个像元采用并行的方式计算出各月月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数参数,完成Copula函数模型的构建;所述四元Copula函数选择四元t-Copula函数,将所述月尺度下的绿度指标、湿度指标、干度指标及热度指标代入上述确定参数后的四元t-Copula函数,计算出对应的综合遥感生态指数。
所述四元t-Copula函数的参数囊括了绿度指标、湿度指标、干度指标及热度指标的数据,因此构建的综合遥感生态指数能够更加全面综合的反映研究区域的生态环境质量。
进一步地,根据Sklar定理,逐像元地计算四个遥感生态指标时间序列的Copula函数模型,所述四个遥感生态指标时间序列的联合分布函数P计算公式如下:
P(x1≤X1,x2≤X2,x3≤X3,x4≤X4)=C[U(X1),V(X2),G(X3),W(X4)]=p
其中,x1为绿度指标时间序列对应变量、x2为湿度指标时间序列对应变量、x3为干度指标时间序列对应变量、x4为热度指标时间序列对应变量;U(x1)为绿度指标的边缘分布函数、V(x2)为湿度指标的边缘分布函数、G(x3)为干度指标的边缘分布函数、W(x4)为热度指标的边缘分布函数;P为绿度指标、湿度指标、干度指标及热度指标的联合分布函数;C为四元t-Copula函数,表达式为:
密度函数的表达式为:
其中,ρ为对角线上的元素全为1的4阶对称正定矩阵,|ρ|表示方阵ρ的行列式,u1为绿度指标边缘分布函数U(x1),u2为湿度指标边缘分布函数V(x2),u3为干度指数边缘分布函数G(x3),u4为热度指标边缘分布函数W(x4);tρ,k表示相关系数矩阵为ρ、自由度为k的标准四元t分布的分布函数,表示自由度为k的标准一元t分布的分布函数的逆函数;/>
进一步地,为了确定四元t-Copula函数模型的参数,根据Kendall相关系数估计相关系数、最大似然函数法估计自由度,算法过程如下:
S41.对四个遥感生态指标运用经验分布或边缘分布函数进行概率转化,所述边缘分布函数包括绿度指标边缘分布函数U(X1)、湿度指标边缘分布函数V(X2)、干度指数边缘分布函数G(X3)和热度指标边缘分布函数W(X4);
S42.根据利用Kendall相关系数的τ来估计四个遥感生态指标的相关系数矩阵/>
S43.运用极大似然估计方法,估计自由度参数k,计算公式如下:
进一步地,所述S3、S4、S5中,遥感栅格影像有数以万计的像元点,每个像元点都要计算一个与之对应的四元t-Copula函数模型,为了实现多个像元的并行计算,提升计算效率,使用Matlab的主从结构的分布式计算方式实现基于多核CPU的并行计算。初始化Matlab并行计算环境时,最初的Matlab进程自动成为主节点,同时初始化多个Matlab计算子节点,所述子节点的个数可根据计算机CPU核心数量设定,利用Parfor(Matlab中关键字)将子节点调用用于计算。Parfor运行之初,主节点会将Parfor循环程序之外变量传递给计算子节点。子节点运算过程时互不干扰,运算完毕,将各子节点得到的结果组合到同一个矩阵变量中,并返回到Matlab主节点,最终计算完毕则关闭计算子节点。
进一步地,所述S5中,将综合遥感生态指数值根据每个像元点的位置、时间信息重新生成可视化栅格数据,并将综合遥感生态指数值归一化;根据综合遥感生态指数划分生态质量等级,将综合遥感生态指数从0到1,以0.2为间隔分成5个等级,所述综合生态质量等级根据综合遥感生态指数由低到高分别是差、较差、中等、良和优5个等级,并根据综合遥感生态指数等级分析综合遥感生态指数时间变化和空间分布规律。
具体的,参考绿都指标、湿度指标、干度指标以及热度指标等标准指数的等级划分对综合遥感生态指数进行了等级划分,其中基于Copula的综合遥感生态指数也称为CRSEI,结果如表1所示:
表1CRSEI生态等级划分
其中Index为CRSEI值。
通过步骤S1-S6计算出CRSEI值后,对照表1可对生态环境质量进行判定。当0.8≤CRSEI≤1时,可认为该地区的生态环境质量为优;当0.6≤CRSEI<0.8时,可认为该地区的生态环境质量为良;当0.4≤CRSEI<0.6时,可认为该地区的生态环境质量为中等;当0.2≤CRSEI<0.4时,可认为该地区的生态环境质量为较差;当0≤CRSEI<0.2时,可认为该地区的生态环境质量为差。
对比例:
分析研究区域的生态环境质量时空演变特征,采用月尺度下的RSEI和CRSEI,RSEI为遥感生态指数。
图3是根据上述方法,计算的月尺度下研究区域2001年10月份的CRSEI、RSEI分级图,并结合对应年份的土地利用类型数据与CRSEI进行比较,以检查计算的CRSEI的性能,图3中,a为2001年10月CRSEI分级图,b为2001年土地利用分布图,c为2001年10月RSEI分级图。总体而言,可以发现研究区域的CRSEI空间分布规律与RSEI指数基本一致,间接验证计算出的CRSEI指数的有效性,并且结合土地利用类型数据分析,林地、湿地以及植被分布的区域CRSEI值较高,反映出这些区域的生态环境较好,而裸土、城市和建筑分布的区域CRSEI值较低,反映出这些区域的生态环境较差,这符合实际的认知情况,直接验证计算的CRSEI能够合理地反映生态环境质量的空间分布真实状况。再由图4中对应影像的黑框区域A可以看出此区域的分布着大量的植被,此时的CRSEI反映该区域的生态环境质量为良或优,而RSEI则反映为中等或较差;再从黑框区域B可以看出此区域的土地利用类型大多为不透水层,此时的CRSEI反映该区域的生态环境质量为中等或较差,而RSEI指数则反映为良或优,这说明CRSEI在反映生态环境质量的空间变化规律上比RSEI指数更为合理准确。
图4是选取上述方法计算得到的月尺度下研究区域2001年1月份、4月份、7月份以及10月份的CRSEI分级图,其中图4中,a为2001年1月CRSEI分级图,b为2001年4月CRSEI分级图,c为2001年7月CRSEI分级图,d为2001年10月CRSEI分级图;从图4可以看出夏季7月份由CRSEI反映出来的遥感生态质量整体较高,春季4月份以及秋季10月份由CRSEI反映出来的遥感生态质量整体上比夏季7月份差一些,而冬季1月份由CRSEI反映出来的遥感生态质量整体上最差,这也符合夏季植物等生长旺盛、雨水较多从而生态环境质量较好,春秋季节次之,冬季较差的季节规律,表明CRSEI在反映生态环境质量的时间变化规律上的合理性。
根据已有研究可知,绿度指标和湿度指标对生态环境质量起积极作用,即绿度指标和湿度指标值较高的区域,生态环境质量好,对应的CRSEI与RSEI值应该较高,绿度指标(NDVI)和湿度指标(WET)值较低的区域,生态环境质量差,对应的CRSEI与RSEI值应该较低;而干度指标(NDBSI)和热度指标(LST)对生态环境质量起消极作用,即干度指标(NDBSI)和热度指标(LST)值较高的区域,生态环境质量差,对应的CRSEI与RSEI值应该较低,NDBSI和LST值较低的区域,生态环境质量好,对应的CRSEI与RSEI值应该较高。根据图5的均值对比图中的矩形框a、b、c以及根据图6的采样点对比图中的矩形框a、b、c、d所示,CRSEI符合上述的研究规律,这说明CRSEI无论是整体月尺度的均值变化还是采样点值变化都具有一定的合理性;其中图5中,右上侧的均值参数由上至下分别为NDVI-MEAN(绿度指标均值)、WET-MEAN(湿度指标均值)、NDBSI-MEAN(干度指标均值)、LST-MEAN(热度指标均值)、CRSEI-MEAN(综合遥感生态指数均值)、RSEI-MEAN(遥感生态指数均值);图6中,右上侧的采样参数由上至下分别为NDVI、WET、NDBSI、LST、CRSEI、RSEI。
以上这些结果证明了本发明基于栅格化的NDVI、WET、NDBSI、LST四种遥感生态指标与四元Copula函数构建的综合遥感生态指数CRSEI可以有效综合绿度指标(NDVI)、湿度指标(WET)、干度指标(NDBSI)、热度指标(LST)四个遥感生态指标的相关性信息,能够表征这四个指标的联合特征,四个生态指标与CRSEI指标都基于栅格可视化,从时间变化和空间演变上更合理准确地反映生态环境质量,为监测和评价生态环境质量提供可靠和有效的工具。
应该理解的是,虽然本发明各实施例的流程图中的各个步骤按照箭头的指示依次显示,但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明,这些步骤的执行并没有严格的顺序限制,这些步骤可以以其它的顺序执行。而且,各实施例中的至少一部分步骤可以包括多个子步骤或者多个阶段,这些子步骤或者阶段并不必然是在同一时刻执行完成,而是可以在不同的时刻执行,这些子步骤或者阶段的执行顺序也不必然是依次进行,而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程,是可以通过计算机程序来指令相关的硬件来完成,上述的程序可存储于一个非易失性计算机可读取存储介质中,该程序在执行时,可包括如上述各方法的实施例的流程。其中,本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用,均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限,RAM以多种形式可得,诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
以上上述的实施例的各技术特征可以进行任意的组合,为使描述简洁,未对上述实施例中的各个技术特征所有可能的组合都进行描述,然而,只要这些技术特征的组合不存在矛盾,都应当认为是本说明书记载的范围。
以上上述的实施例仅表达了本发明的几种实施方式,其描述较为具体和详细,但并不能因此而理解为对本发明专利范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变形和改进,这些都属于本发明的保护范围。因此,本发明专利的保护范围应以所附权利要求为准。
以上上述的仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,将Copula函数基于遥感栅格数据进行并行计算,提高计算效率并且将计算结果可视化,同时Copula函数作为一种新的多变量分析工具对遥感生态指标进行拟合分析,除了囊括遥感生态指标的更多信息外,还能够考虑指标变量的相关性,从时间变化和空间演变上更合理准确地反映生态环境质量,为生态环境的监测、评价以及治理提供可靠、有力的技术支持,包括如下步骤:
S1.收集研究区域内十年的多光谱遥感栅格影像数据、绿度指标月合成数据、热度指标月合成数据及年尺度的土地利用数据;用所述多光谱遥感栅格影像数据分别计算月尺度下的湿度指标和干度指标;并对四个遥感生态指标数据进行归一化处理,所述四个遥感生态指标包括绿度指标、湿度指标、干度指标和热度指标;
S2.将四个遥感生态指标的栅格影像数据以像元为单位逐个读取成四个遥感生态指标时间序列,构成十年120个月尺度的四个遥感生态指标时间序列,所述四个遥感生态指标时间序列包括绿度指标时间序列、湿度指标时间序列、干度指标时间序列及热度指标时间序列;
S3.对每个像元的四个遥感生态指标时间序列分别计算其遥感生态指标边缘分布函数,所述四个遥感生态指标边缘分布函数包括绿度指标边缘分布函数、湿度指标边缘分布函数、干度指标边缘分布函数和热度指标边缘分布函数;
S4.利用每个像元的四个遥感生态指标边缘分布函数确定四元Copula函数的模型参数,构建月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数模型,将所述四个指标时间序列重新代入确定参数后的四元Copula函数模型,计算出综合遥感生态指数;
S5.将所述综合遥感生态指数结果重新生成可视化栅格影像数据,并将所述可视化栅格影像数据归一化;
S6.设定综合遥感生态指数的阈值,将所述综合遥感生态指数根据所述阈值划分成不同生态环境质量等级;
所述S1中,能够方便地获取多光谱遥感影像数据用于指标的计算,所述遥感影像数据包括含红光、绿光、蓝光、第一近红外光、第二近红外光、第一短波红外、第二短波红外的7波段的8天合成数据;
所述S2中,对四元Copula函数实现基于遥感栅格数据的读取计算,使用Matlab将归一化后的四个遥感生态指标的栅格影像数据读取为矩阵形式,按照同一位置的像元点逐个读取十年120个月的月尺度的四个遥感生态指标时间序列。
2.根据权利要求1所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,所述S1中,采用多线程并行的方式高效地计算遥感生态指标并归一化,所述月尺度下的湿度指标和所述月尺度下的干度指标采用python开启线程池多线程并行计算方式如下:
S11.根据遥感数据计算湿度指标,计算公式为
WETNESS=0.1147×Red+0.2408×NIR1+0.3283×Blue+0.3132×Gree-0.3122×NIR2-0.6416×SWIR1-0.5087×SWIR2
其中,SWIR1、SWIR2、NIR1、NIR2、Red、Blue、Green分别为第一短波红外、短二短波红外、第一近红外光、第二近红外光、红光、蓝光、绿光波段的反射率;
S12.根据遥感影像数据计算干度指标,计算公式为:
BSI=[(SWIR1+Red)-(NIR1+Blue)]/[(SWIR1+Red)+(NIR1+Blue)]
NDBSI=(BSI+IBI)/2
其中,SWIR1、NIR1、Red、Green、Blue分别为第一短波红外、第一近红外光、红光、绿光、蓝光波段的反射率;
S13.将所述湿度指标和干度指标去除异常值,并利用最大值合成法将数据合成月尺度下的湿度指标和月尺度下的干度指标;
S14.将所述月尺度下的湿度指标和月尺度下的干度指标进行归一化处理,归一化处理计算公式为:
NIi=(Ii-Imin)/(Imax-Imin)
其中,NIi为某一正规化后的指标值,ii为该指标在像元i处的值;Imax为该指标的最大值;Imin为该指标的最小值。
3.根据权利要求1所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,所述S3中,按照逐像元的方式对四个遥感生态指标组成的四组时间序列计算边缘分布函数,使用Gringorten公式,逐个像元点分别计算所述四个遥感生态指标边缘分布函数,当时间序列共有N项,所述时间序列x1≤x2≤…≤xi≤…≤xN中小于等于xi的出现次数为i时,其公式如下:
4.根据权利要求1所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,所述S4中,使用四元Copula函数作为边缘分布函数的连接工具,栅格影像中每个像元采用并行的方式计算出各月月尺度下的绿度指标、湿度指标、干度指标及热度指标的四元Copula函数参数,完成Copula函数模型的构建;所述四元Copula函数选择四元t-Copula函数,将所述月尺度下的绿度指标、湿度指标、干度指标及热度指标代入上述确定参数后的四元t-Copula函数,计算出对应的综合遥感生态指数。
5.根据权利要求4所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,根据Sklar定理,逐像元地计算四个遥感生态指标时间序列的Copula函数模型,所述四个遥感生态指标时间序列的联合分布函数P计算公式如下:
P(x1≤X1,x2≤X2,x3≤X3,x4≤X4)=C[U(X1),V(X2),G(X2),W(X4)]=p
其中,x1为绿度指标时间序列对应变量、x2为湿度指标时间序列对应变量、x3为干度指标时间序列对应变量、x4为热度指标时间序列对应变量;U(x1)为绿度指标的边缘分布函数、V(x2)为湿度指标的边缘分布函数、G(x)为干度指标的边缘分布函数、W(x4)为热度指标的边缘分布函数;P为绿度指标、湿度指标、干度指标及热度指标的联合分布函数;C为四元t-Copula函数,表达式为:
密度函数的表达式为:
其中,ρ为对角线上的元素全为1的4阶对称正定矩阵,|ρ|表示方阵ρ的行列式,u1
为绿度指标边缘分布函数U(x1),u2为湿度指标边缘分布函数V(x2),u3为干度指标边缘分布函数G(x3),u4为热度指标边缘分布函数W(x4);tρ,k表示相关系数矩阵为ρ、自由度为k的标准四元t分布的分布函数,表示自由度为k的标准一元t分布的分布函数的逆函数;
6.根据权利要求5所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,为了确定四元t-Copula函数模型的参数,根据Kendall相关系数估计相关系数、最大似然函数法估计自由度,方法如下:
S41.对四个遥感生态指标运用经验分布或边缘分布函数进行概率转化,所述边缘分布函数包括绿度指标边缘分布函数U(X1)、湿度指标边缘分布函数V(X2)、干度指数边缘分布函数G(X3)和热度指标边缘分布函数W(X4);
S42.根据利用Kendall相关系数的τ来估计四个遥感生态指标的相关系数矩阵/>i=1,2,3,4;j=1,2,3,4;
S43.运用极大似然估计方法,估计自由度参数k,计算公式如下:
7.根据权利要求4所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,所述S3、S4、S5中,遥感栅格影像有数以万计的像元点,每个像元点都要计算一个与之对应的四元t-Copula函数模型,为了实现多个像元的并行计算,提升计算效率,使用Matlab的主从结构的分布式计算方式实现基于多核CPU的并行计算;初始化Matlab并行计算环境时,最初的Matlab进程自动成为主节点,同时初始化多个Matlab计算子节点,所述子节点的个数可根据计算机CPU核心数量设定,利用Parfor将子节点调用用于计算;Parfor运行之初,主节点会将Parfor循环程序之外变量传递给计算子节点;子节点运算过程时互不干扰,运算完毕,将各子节点得到的结果组合到同一个矩阵变量中,并返回到Matlab主节点,最终计算完毕则关闭计算子节点。
8.根据权利要求1所述一种基于四元Copula的综合遥感生态指数并行计算方法,其特征在于,所述S5中,将综合遥感生态指数值根据每个像元点的位置、时间信息重新生成可视化栅格数据,并将综合遥感生态指数值归一化;根据综合遥感生态指数划分生态质量等级,将综合遥感生态指数从0到1,以0.2为间隔分成5个等级,所述生态质量等级根据综合遥感生态指数由低到高分别是差、较差、中等、良和优5个等级,并根据综合遥感生态指数等级分析综合遥感生态指数时间变化和空间分布规律。
CN202211514932.9A 2022-11-25 2022-11-25 一种基于四元Copula的综合遥感生态指数并行计算方法 Active CN115795237B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211514932.9A CN115795237B (zh) 2022-11-25 2022-11-25 一种基于四元Copula的综合遥感生态指数并行计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211514932.9A CN115795237B (zh) 2022-11-25 2022-11-25 一种基于四元Copula的综合遥感生态指数并行计算方法

Publications (2)

Publication Number Publication Date
CN115795237A CN115795237A (zh) 2023-03-14
CN115795237B true CN115795237B (zh) 2023-11-03

Family

ID=85443350

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211514932.9A Active CN115795237B (zh) 2022-11-25 2022-11-25 一种基于四元Copula的综合遥感生态指数并行计算方法

Country Status (1)

Country Link
CN (1) CN115795237B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117274082A (zh) * 2023-09-18 2023-12-22 深圳市中科云驰环境科技有限公司 一种基于遥感生态指数的水体生态环境质量分析方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115099453A (zh) * 2022-05-06 2022-09-23 河海大学 一种多变量栅格化卫星遥感综合干旱风险评估方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7260766B2 (ja) * 2019-03-29 2023-04-19 日本電信電話株式会社 シミュレーション値算出方法及びシミュレーション値算出装置

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115099453A (zh) * 2022-05-06 2022-09-23 河海大学 一种多变量栅格化卫星遥感综合干旱风险评估方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Analysis of Vegetation Vulnerability Dynamics and Driving Forces to Multiple Drought Stresses in a Changing Environment;Xiaoting Wei等;《remote sensing》;全文 *
Copula-Based Drought Analysis Using Standardized Precipitation Evapotranspiration Index: A Case Study in the Yellow River Basin, China;Fei Wang等;《water》;全文 *
区域生态环境变化的遥感评价指数;徐涵秋;《中国环境科学》;参见2-5页 *
融合气象和遥感植被信息的综合干旱指数构建及精度评价;王丹钰等;武汉大学学报(信息科学版);1-7页 *

Also Published As

Publication number Publication date
CN115795237A (zh) 2023-03-14

Similar Documents

Publication Publication Date Title
Jaseena et al. Decomposition-based hybrid wind speed forecasting model using deep bidirectional LSTM networks
Guermoui et al. A comprehensive review of hybrid models for solar radiation forecasting
Khodayar et al. Spatio-temporal graph deep neural network for short-term wind speed forecasting
Papagiannopoulou et al. A non-linear Granger-causality framework to investigate climate–vegetation dynamics
Zhang et al. Regionalization and spatial changing properties of droughts across the Pearl River basin, China
Fan et al. Prediction of crop yield using big data
Hu et al. Research and application of a hybrid model based on Meta learning strategy for wind power deterministic and probabilistic forecasting
CN115795237B (zh) 一种基于四元Copula的综合遥感生态指数并行计算方法
US11264799B2 (en) Submodular load clustering with robust principal component analysis
Shahrin et al. Agricultural analysis and crop yield prediction of habiganj using multispectral bands of satellite imagery with machine learning
Ryu et al. Residential load profile clustering via deep convolutional autoencoder
CN116128141A (zh) 风暴潮预测方法、装置、存储介质及电子设备
Balti et al. Big data based architecture for drought forecasting using LSTM, ARIMA, and Prophet: Case study of the Jiangsu Province, China
Karevan et al. Black-box modeling for temperature prediction in weather forecasting
Akrami et al. Graph-based local climate classification in Iran
Fan et al. Spatio-temporal denoising graph autoencoders with data augmentation for photovoltaic timeseries data imputation
CN104699979B (zh) 基于复杂网络的城市湖库藻类水华混沌时间序列预测方法
Xie et al. Investigating long-term trends of climate change and their spatial variations caused by regional and local environments through data mining
Kumawat et al. Time-Variant Satellite Vegetation Classification Enabled by Hybrid Metaheuristic-Based Adaptive Time-Weighted Dynamic Time Warping
Ayoub Self-organizing profiles to characterize representative temporal settings for daylight simulations
CN112819208A (zh) 一种基于特征子集耦合模型的空间相似性地质灾害预测方法
CN116937559A (zh) 基于循环神经网络和张量分解的电力系统负荷预测系统和方法
Wijayanto et al. Estimating Rice production using machine learning models on multitemporal Landsat-8 satellite images (case study: Ngawi regency, East Java, Indonesia)
CN115564989A (zh) 面向土地利用分类的随机森林算法
CN114782835A (zh) 作物倒伏面积比例检测方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant