基于空间几何原理的四维遥感生态指数构建方法
技术领域
本发明涉及遥感生态指数构建领域，具体涉及一种基于空间几何原理的四维遥感生态指数构建方法。
背景技术
区域生态评价指标近几年来已经逐渐被人们所重视，遥感作为定量评估区域生态环境的重要手段，已经得到了广泛运用。遥感以一定的周期重复获取数据，同时可大面积同步观测，并真实地反映区域生态环境变化，快速地获取目标物属性特征，从而进行生态评价。
遥感数据在生态状况评价中最多的是利用遥感数据从单一的角度评价生态中的某一个方面，但是无法整体评价生态状况。有的研究人员通过主成分变换或者因子分析方法集成多个指标构建反映生态状况的综合遥感指数，如RSEI(徐涵秋，2013)，EEM(Zhang,2016)等，但这些指数物理意义不明确；VMTEI(周紫燕等,2020)利用的空间几何原理，但是只考虑了植被覆盖、地表湿度和地表温度3个指标。本发明基于多波段遥感数据，耦合了代表绿度的垂直植被指数PVI、代表干度的改进型垂直干旱指数MPDI(PVI和MPDI在空间中垂直)以及地表干化度NDSI和地表温度LST，引入空间几何原理构建空间几何生态指数SGEI，提高了生态指数SGEI的物理意义，可以综合反映区域的生态环境状况的空间分布。
发明内容
有鉴于此，本发明的目的在于提供一种基于空间几何原理的四维遥感生态指数构建方法，以克服现有技术中存在的缺陷。
为实现上述目的，本发明采用如下技术方案：
一种基于空间几何原理的四维遥感生态指数构建方法，包括以下步骤：
步骤S1：获取研究区包括蓝光、绿光、红光、近红外、短波红外、热红外波段的遥感数据并进行预处理，进一步得到据红光波段和近红外波段二维空间散点图，并根据该散点图拟合土壤线并求出土壤线方程；
步骤S2：计算垂直植被指数PVI；
步骤S3：计算垂直干旱指数PDI，并设混合像元两端极限是裸土和植被，用混合像元的PDI的值减去植被覆盖PDI的值得到改进型垂直干旱指数MPDI；
步骤S4：计算地表干化度NDSI；
步骤S5：计算地表温度LST；
步骤S6：以垂直植被指数PVI、改进型垂直干旱指数MPDI、地表干化度指数NDSI和地表温度指数LST为坐标轴，构建四维特征空间；
步骤S7：计算基于植被、地表干旱度、地表干化度和地表温度的遥感生态指数SGEI。
进一步的，所述步骤S1具体为：
步骤S11：获取研究区红光波段蓝光、绿光、红光、近红外、短波红外、热红外波段的研究区遥感数据，对遥感数据进行几何校正、辐射定标、大气校正等预处理，提取蓝光、绿光、红光、近红外、短波红外、热红外波段的反射率数据；
步骤S12：以红光波段反射率数据为横坐标，以近红外波段反射率数据为纵坐标，构建二维空间散点图；
步骤S13：基于构建出的空间二维散点图的分布，通过横坐标对应的最小纵坐标值构建(X1,Y1)点集，并通过最小二乘法拟合土壤线，得到土壤线的方程：
ρNIR＝M*ρRed+I (1)
式中，ρRed和ρNIR分别表示红光波段和近红外波段的反射率值(值域范围为[0,1])；M和I指的是土壤线方程的斜率和截距。
进一步的，所述步骤S2具体为：根据步骤S1确定的土壤线方程，垂直植被指数PVI的计算方法如下：
式中，ρNIR表示近红外波段的反射率，ρRed表示红光波段的反射率；M和I分别代表土壤线方程的斜率和截距。
进一步的，所述步骤S3具体为：
步骤S31：植被覆盖度指数fv通过归一化植被指数NDVI计算得出，用影像中全植被覆盖和无植被覆盖的裸土区域通过比例计算得出，植被覆盖度选择采用Baret等的半经验植被覆盖度计算方法，其表达式如下：
式中，NDVIMAX和NDVIMIN分别为100％植被覆盖下NDVI和裸土对应的NDVI；
步骤S32：根据步骤S1确定的土壤线方程，垂直干旱指数PDI的计算方法如下：
式中，ρNIR表示近红外波段的反射率，ρRed表示红光波段的反射率；M和I分别代表土壤线方程的斜率和截距；
步骤S33：利用改进型垂直干旱指数MPDI，根据步骤S1确定的土壤线方程，MPDI表达式为：
式中，ρS,Red表示裸土在红光波段的反射率，ρS,NIR表示裸土在近红外波段的反射率；
步骤S34：设混合像元两端极限是裸土和植被，则混合像元反射率是两者线性结合，以此推导出MPDI，其表达式为：
式中，ρRed和ρNIR分别为遥感影像中经过大气校正的红光、近红外波段反射率；ρV,Red和ρV,NIR分别为红光波段和近红外波段的植被反射率；M为土壤基线斜率。
进一步的，所述步骤S4具体为：
SI代表的是裸露土壤(沙土、裸土等)造成的地表干化，其计算方法为：
IBI是建筑物指数，代表的是城区的建筑用地造成的干化，其计算方法为：
式中ρSWIR1、ρblue、ρgreen分别为短波红外、蓝光和绿光波段的地表反射率，其他变量同上。
则地表干化度的计算方法如下：
进一步的，当所选择的数据为2000年前的Landsat5遥感数据时，所述步骤S5具体为：
(1)图像辐射定标，对Landsat5热红外波段进行辐射定标；
(2)植被覆盖度计算，利用混合像元分解法计算植被覆盖度Fv，将影像中的地类大致分为水体、植被和建筑，具体的计算公式如下：
Fv＝[(NDVI-NDVISoil)/(NDVIVeg-NDVISoil)] (10)
式中，NDVI为归一化植被指数，NDVISoil为无植被覆盖土壤区域的NDVI值，NDVIVeg则为完全被植被覆盖区域的NDVI值；
(3)地表比辐射率ε计算
水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率估算则分别根据下式进行计算：
εsurface＝0.9625+0.0614Fv-0.0461Fv 2 (11)
εbuilding＝0.9589+0.086Fv-0.0671Fv 2 (12)
式中，εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率；
(4)像元在传感器处的光谱辐射值Lλ计算
Lλ＝gain·Qλ+bias (13)
式中：Lλ为热红外波段的像元在传感器处的光谱辐射值，单位为W/(m2.ster.μm)；gain为对应波段的增益值；bias为对应波段的偏置值；Qλ为像元对应DN值；
(5)计算传感器处亮温值
Tλ＝K2/ln(K1/Lλ+1) (14)
式中，Tλ是传感器处亮温值，单位为K；
LST＝Tλ/[1+(λ·Tλ/ρ)lnε] (15)
式中，λ为热红外波段的中心波长或有效波长，单位为μm，；K1和K2分别为定标参数。
进一步的，当所选择的数据为2000年后的Landsat5遥感数据时，所述步骤S5具体为：
(1)图像辐射定标，对Landsat5热红外波段进行辐射定标；
(2)植被覆盖度计算，利用混合像元分解法计算植被覆盖度Fv，将影像中的地类大致分为水体、植被和建筑，具体的计算公式如下：
Fv＝[(NDVI-NDVISoil)/(NDVIVeg-NDVISoil)] (10)
式中，NDVI为归一化植被指数，NDVISoil为无植被覆盖区域的NDVI值，NDVIVeg则为完全被植被覆盖的区域的NDVI值；
(3)地表比辐射率ε计算
水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率估算则分别根据下式进行计算：
εsurface＝0.9625+0.0614Fv-0.0461Fv 2 (11)
εbuilding＝0.9589+0.086Fv-0.0671Fv 2 (12)
式中，εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率；
(4)计算相同温度下黑体的辐射亮度值
卫星传感器接收到的热红外辐射亮度值Lλ包括三部分：大气向上辐射亮度，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。其计算方法为：
Lλ＝[ε·B(TS)+(1-ε)L↓]·τ+L↑ (13)
这里，ε为地表辐射率，TS为地表真实温度，B(TS)为黑体在TS的热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS)＝[Lλ-L↑-τ·(1-ε)L↓]/τ·ε (14)
(5)反演地表温度
在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数，求得地表真实温度TS：
TS＝K2/ln(K1/B(TS)+1) (15)
式中，B(TS)为黑体辐射亮度。
进一步的，当所选择的数据为Landsat8遥感数据时，所述步骤S5具体为：
(1)图像辐射定标，对Landsat8热红外波段进行辐射定标；
(2)植被覆盖度计算，利用混合像元分解法计算植被覆盖度Fv，将影像中的地类大致分为水体、植被和建筑，具体的计算公式如下：
Fv＝[(NDVI-NDVISoil)/(NDVIVeg-NDVISoil)] (10)
式中，NDVI为归一化植被指数，NDVISoil为无植被覆盖区域的NDVI值，NDVIVeg则为完全被植被覆盖的区域的NDVI值；
(3)地表比辐射率ε计算，水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率估算则分别根据下式进行计算：
εsurface＝0.9625+0.0614Fv-0.0461Fv 2 (11)
εbuilding＝0.9589+0.086Fv-0.0671Fv 2 (12)
式中，εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率；
(4)计算相同温度下黑体的辐射亮度值
卫星传感器接收到的热红外辐射亮度值Lλ包括三部分：大气向上辐射亮度，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量；其计算方法为：
Lλ＝[ε·B(TS)+(1-ε)L↓]·τ+L↑ (13)
这里，ε为地表辐射率，TS为地表真实温度，B(TS)为黑体在TS的热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS)＝[Lλ-L↑-τ·(1-ε)L↓]/τ·ε (14)
在NASA官网(http://atmcorr.gsfc.nasa.gov/)中输入成影时间以及中心经纬度，则会提供上式中所需要的参数；
(5)反演地表温度
在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数，求得地表真实温度TS：
TS＝K2/ln(K1/B(TS)+1) (15)
式中，B(TS)为黑体辐射亮度；对于Landsat8 TIRS而言，K1＝774.89W/(m2*μm*sr)，K2＝1321.08K。
10.根据权利要求1所述的基于空间几何原理的四维遥感生态指数构建方法，其特征在于，所述步骤S7具体为：
步骤S71：将每个轴的值归一化至[0,1/2]，然后再计算SGEI，使SGEI的取值范围为[0,1]。其中，PVI进行正向归一化，LST、MPDI、NDSI则进行反向归一化
正向归一化公式：
反向归一化公式：
式中，Normalizedpi为归一化后的值；
步骤S72：在步骤6中所创建的三维特征空间中，构建表征生态质量的指标SGEI的表达式为：
SGEI即为基于绿度-土壤干度-地表温度-地表干化度的四维遥感生态指数，数值范围为[0,1]，SGEI越接近1生态状况越好，反之，生态状况越差。
本发明与现有技术相比具有以下有益效果：
本发明既考虑了与生态密切相关的多个指标，又通过引入空间几何原理提高了SGEI的物理意义。另外，SGEI不仅可以作为一个量化指标来描述空间生态质量，还可以将生态质量空间可视化，并进行时空分析。
附图说明
图1为本发明构建方法流程图；
图2为本发明一实施例中空间几何生态指数SGEI构建过程中坐标系的图解；
图3为本发明一实施例中拟合土壤线方程图；
图4为本发明一实施例中PVI、MPDI、NDSI和LST的空间分布图；
图5为本发明一实施例中SGEI的空间分布图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
请参照图1，本发明提供种基于空间几何原理的四维遥感生态指数(SGEI)构建方法，包括如下步骤：
步骤S1：获取研究区Landsat包括蓝光、绿光、红光、近红外、短波红外、热红外波段的遥感数据并进行预处理，根据红光波段和近红外二维空间散点图拟合土壤线并求出土壤线方程：
步骤S1的具体步骤如下：
步骤S11：获取研究区Landsat的蓝光、绿光、红光、近红外、短波红外、热红外波段的研究区遥感数据，对遥感数据进行几何校正、辐射定标、大气校正等预处理，提取蓝光、绿光、红光、近红外、短波红外、热红外波段的反射率数据。
步骤S12：根据红光波段反射率和近红外波段反射率的数据构建二维空间散点图。
步骤S13：根据构建出的空间二维散点图的分布，拟合土壤线方程：
ρNIR＝1.1373*ρRed-0.0825 (1)
式中ρNIR和ρRed分别表示近红外波段和红光波段的反射率值(均为0到1)
步骤S2：根据土壤线的斜率和截距，选取任意一点到土壤线的垂直距离表示区域内植被覆盖情况，不同类型植被反射率到土壤线的距离都不同，计算垂直植被指数PVI。
步骤S2的具体步骤如下：
根据步骤S13确定的土壤线方程，垂直植被指数PVI的计算方法如下：
式中，ρNIR表示近红外波段的反射率，ρRed表示红光波段的反射率；M和I分别代表土壤线方程的斜率和截距。根据方程(1)，此实例中，M＝1.1373，I＝0.0825，则PVI的具体表达式为：
步骤S3：用混合像元的PDI的值减去植被覆盖PDI的值得到MPDI。
所述步骤S3的具体步骤如下：
计算植被覆盖度，其表达式如下：
式中，NDVIMAX和NDVIMIN分别为100％植被覆盖下NDVI和裸土对应的NDVI。
根据步骤13确定的土壤线方程，MPDI表达式为：
式中，ρS,Red表示裸土在红光波段的反射率，ρS,NIR表示裸土在近红外波段的反射率。
改进的部分主要是考虑了植被覆盖的影响，MPDI表达式为：
式中，ρRed和ρNIR分别为遥感影像中经过大气校正的红光、近红外波段反射率，ρV,Red和ρV,NIR分别为红光波段和近红外波段的植被反射率。得到其数值表达式：
将ρV,Red和ρV,NIR分别采用0.05和0.5。M为土壤基线斜率，如果要计算多年的，可考虑利用土壤线基线平均斜率。
步骤S4：计算地表干化度NDSI，由不透水面建筑指数IBI(An Index-basedBuilt-up Index)和裸土指数SI(Soil Index)合成。
步骤S41：SI代表的是裸露土壤(沙土、裸土等)造成的地表干化，其计算方法为：
步骤S42：IBI代表的是建筑物指数，代表的是城区的建筑用地造成的干化，其计算方法为：
式中ρSWIR1、ρblue、ρgreen分别为短波红外、蓝光和绿光波段的地表反射率，其他变量同上。
步骤S43：地表干化度的计算方法如下：
步骤S5：计算地表温度LST。
本实例中，采用Landsat8 TIRS计算LST。
步骤S6：以垂直植被指数PVI、改进型垂直干旱指数MPDI、地表干化度指数NDSI和地表温度指数LST构建四维特征空间。
步骤S7：计算基于植被、地表干旱度、地表干化度和地表温度的遥感生态指数SGEI。
为了使指标更具合理性，便于进行比较，将每个轴的值归一化至[0,1/2]，然后再计算SGEI，使SGEI的取值范围为[0,1]。其中，PVI进行正向归一化，LST、MPDI、NDSI则进行反向归一化。
正向归一化公式：
反向归一化公式：
式中，Normalizedpi为归一化后的值。
在步骤7中所创建的四维特征空间中，构建表征生态质量的指标SGEI的表达式为：
SGEI即为基于绿度-土壤干度-地表温度-地表干化度的四维遥感生态指数，数值范围为[0,1]，SGEI越接近1生态状况越好，反之，生态状况越差。
本实施例中，利用福建省福州市内陆区域2013年10月的Landsat8数据来计算空间几何遥感生态指数SGEI。图4为PVI、MPDI、NDSI和LST的空间分布图，图5为本发明实施例中SGEI的空间分布图。
从图中可知，除主城区外，福州内陆区域植被覆盖度均较高；地表干度在闽江两岸主城区较高，在福州东南区域也较高；地表干化度与植被指数呈较高的负相关关系，地表温度和地表干度呈现出较明显的正相关关系；空间几何生态指数SGEI的分布符合垂直植被指数、改进型垂直干旱指数、地表干化度和地表温度的空间分布，整体来说，绿度越高、地表干旱度越低、地表干化度越低、地表温度越低的区域其表征的生态状况越好，反之生态状况越差。
以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。
