基于GEE的遥感生态指数时序计算

骆泓鉴， 明冬萍， 徐录

(中国地质大学(北京)信息工程学院，北京 100083)

0 引言

准确、良好的生态评价对地区的发展、规划起重要的支撑作用。通过地理信息的搜集、处理、分析，人们往往能避免主观判断的失误，对生态环境做出客观的评价。

徐涵秋提出了一种基于遥感技术的生态评价方法，以自然因子为主构建遥感生态指数(remote sensing ecological index, RSEI)，从而对城市的生态状况进行快速监测与评价[1]。不同的学者用RSEI或改进的RSEI进行生态评价或其他研究。贾浩巍等[2]以1995年、2005年、2015年为主要年份研究都兰县生态环境； 李红星等[3]以2005年、2010年、2015年为研究时间评估武汉市生态环境； 缪鑫辉等[4]选取2000年、2009年、2017 年的影像，分析甬江流域生态环境变化； 孙从建等[5]定量分析了 2002年、2009年、2017 年3个年份黄土高原的生态安全，指出治理区域。非时序的生态研究都是基于长时间间隔的影像开展的，是跳跃的，只能得出生态趋势，无法分析变化细节。

谷歌地球引擎(Google Earth Engine, GEE)是用于行星级地理空间分析的云平台[6]，能提供丰富的开源数据和强大的计算服务。GEE使得长时间序列的庞大计算成为可能，使得RSEI时序计算变得方便。已有学者基于GEE进行RSEI的连续分析，并取得相应的成果[7-9]。但这些研究在计算方式上仍采用传统的RSEI。传统的RSEI算法在分量的计算方式上会随传感器的变化而变化，加上不同传感器之间的成像差异，得到的RSEI计算结果差异较大，不利于长时间跨度的连续分析。综上所述，尚未有RSEI的时序计算研究。本文将探寻一种基于GEE的RSEI时序计算方式，在分量计算及其归一化方面对RESI算法进行改进，并进行分析验证。

1 研究区概况和数据预处理

1.1 研究区概况

新疆维吾尔自治区奎屯市坐落在天山北麓中段，准噶尔盆地西南部，高程在450～530 m之间，地势南高北低，西高东低，无明显地形起伏变化，不必考虑地形对研究的影响。奎屯市境内为大陆性气候，日照长，雨少干热，较易获取长时间无云或者少云的遥感影像。奎屯市近30 a来生态、发展都有巨大变化，利用遥感影像分析生态变化较为容易。

1.2 数据源与GEE数据预处理

研究数据采用Landsat的地表反射率产品和MODIS的地表温度数据产品(MOD11A1)。根据研究的需要，选取1989—2019年间每年6月1日—10月30日期间的影像进行计算分析。计算指数之前先要获取去云影像。GEE提供了针对Landsat 地表反射率产品的去云算法。但在实验中发现，上述算法在处理异物同谱物体时效果欠佳，容易导致像元丢失。本文提出了减云拼接的去云方法，其流程如图1所示，其中CLOUD-COVER代表云覆盖程度。完善之后的去云方法能补全缺失影像，如图2所示。

图1 减云拼接方法

(a) GEE提供的样例算法去云效果(b) 本文提出的减云拼接方法去云效果

2 RSEI时序计算

RSEI的核心思想是通过空间降维的方式，从一系列能在一定程度上反映生态环境质量的分量中，提取出最能表达当地生态环境的综合结果。RSEI的基本计算流程如图3所示。在众多的自然因素以及人文因素中，绿度、湿度、干度、温度是与人类活动密切相关的因素，常用于各类生态评价[1]。运用函数的方法对4大分量进行降维压缩，最终得到表征生态环境的RSEI。

图3 RSEI生产流程

2.1 RSEI分量计算

1)绿度计算。绿度的计算通常考虑植被指数。归一化差值植被指数(normalized difference vegetation index, NDVI)在生产中得到了最广泛的应用[10]，能在一定程度上表征生态环境。NDVI的计算公式为：

，

(1)

式中ρred和ρnir分别为红光、近红外波段地表反射率。

2)湿度计算。本文使用缨帽变化中的湿度分量进行湿度(Wetness)的计算。不同传感器之间的缨帽变化参数存在较大的差异[11-15]。在时序计算过程中，为了消除不同缨帽变化参数在计算过程中带来的结果差异，本文选择文献[15]中的缨帽变化湿度参数作为湿度计算的模型。其计算公式为：

，

(2)

式中ρblue，ρgreen，ρswir1，ρswir2分别为蓝光、绿光、第1短波红外、第2短波红外波段的地表反射率。

3)干度计算。徐涵秋用建筑物指数(index-based built-up index, IBI)和裸土指数(bare soil index, BSI)计算干度指数(normalized difference building and soil index, NDBSI)[16-18]。干度的计算公式为：

NDBSI=γIBI+(1-γ)BSI

，

(3)

(4)

，

(5)

式中:γ为参数系数，通常取值0～1.0，这里根据前人经验取γ= 0.5[1]；ρmid为中红外波段地表反射率，但在Landsat中无法找到对应的波段，通常采用第1短波红外波段反射率(ρswir1)来代替。

4)温度计算。温度指的是地表温度(land surface temperature, LST)。它可以通过遥感的方式，用获取到的热红外波段反演得到。其原理是基于普朗克定律量化所构成的热辐射传输方程。

前人提出不同的温度反演方法，有单窗算法、单通道算法和劈窗算法等。本文的实验数据是处理之后的拼接影像，不同像元来自不同时间段，而辐射传输模型、单窗算法和单通道算法都需要知道具体时间下的大气参数数据，不适合本次实验。

本文采用亮温转换的方式来进行温度的反演，其公式为：

，

(6)

式中： C2为普朗克函数常量，14 387.79 μm·K；B为DN值经过辐射定标和转换后的亮度温度；ε为计算得到的比辐射率；λ为波段的有效波长，不同传感器的有效波长不同[19]。

2.2 整体归一化

本次实验是基于不同传感器、不同拍摄时间的拼接影像开展的。不同传感器结果之间存在较大差异，如图4所示。常见的归一化方法容易造成标准的不统一，不利于后续计算分析。

(a) 2016年Landsat7影像湿度归一化结果 (b) 2016年Landsat8影像湿度归一化结果

为解决上述问题，本文用不同年份计算结果的整体最值进行归一化，为不同传感器的计算结果制定相同的标准，方便进行对比分析。整体最值归一化的公式为：

，

(7)

式中：Ni为归一化后结果值；Vi为原始像元值；VW-max和VW-min分别为所有计算结果中的最大值和最小值。

2种不同的归一化方式的差别如图5所示。

(a) 单一影像归一化 (b) 整体最值归一化

2.3 RSEI最终计算

绿度、湿度、干度、温度均能在一定程度上表征生态优劣。主成分分析是通过在多维空间中变换坐标方向的一种识别数据模式的方法，是一种突出数据相似性和差异性的方式来表达数据的方式[20]。通过主成分分析，能够综合考虑4大分量的作用，构建一个表征生态的指数，用于生态评价。RSEI计算公式为：

RSEI0=f(NDVI,NDBSI,Wetness,LST) ，

(8)

。

(9)

式(8)表示对归一化所得的4项指数进行主成分分析，并取第一主成分分量作为RSEI0。式(9)表示将主成分分析所得的RSEI0进行归一化处理，得到最终的RSEI。计算结果显示，第一主成分分量的平均信息包含量超过了85.37%，能够较好地代替4大分量，实现数据压缩。

3 RSEI时序分析实验

不同于传统的RSEI计算方式，RSEI时序计算在温度、湿度分量计算以及归一化方式上都有变化。本节将通过对比实验进行方法验证。

3.1 湿度分量验证

由于计算湿度分量采用的参数不统一，所以本文对文献[12—15]的4种缨帽参数的计算方式进行了比较，使用2001年和2002年的Landsat7影像，2010年和2011年的Landsat5影像，以及2017年和2018年的Landsat8影像进行测试。

对比表1发现，文献[13]和文献[14]中的湿度计算结果在上下兼容性上略显不足，不适用于长时间序列分析。根据地物的反射波谱可知，含水地域的地表反射率普遍较低，在真彩色影像上表现为深色、暗色调。以此为依据对影像进行目视判读，观察真彩色影像中的明亮区域，并与湿度影像中的对应区域进行比较。发现文献[12]在一些地区的结果与目视判读不相符，而本文选择的参数的计算结果更符合目视判读。综上，本文所选的文献[15]中的湿度参数在传感器兼容性以及合理性上的表现更据优势，更适合时序分析。

表1 部分年份湿度分量结果对比

3.2 温度分量验证

选择2007—2011年Landsat5和Landsat7影像及2013—2018年Landsat8影像进行测试，在温度计算上比较文献[21]中的模型和亮温转换模型，在比辐射率计算上比较NDVI阈值法、文献[22]和文献[23]中的方法[22-25]。计算结果通过散点图的方式和相同时期下的MOD11A1数据进行对比。

表2 基于2种LST模型的3种比辐射率算法测试结果

对比表2结果，优先考虑方差较小的情况，并协同比较平均差值和中位差值，验证得到本文选择的亮温转换模型和NDVI阈值法作为LST的反演方法是最为合适的。

3.3 时序分析验证

3.3.1 贡献率验证

用主成分分析的方法进行数据压缩，压缩后的第一主成分分量能在最大程度上代表原始数据的主要信息，从而代替生态分量表征生态环境。贡献率越高，则第一主成分分量的代表能力越强。表3是2种不同RSEI计算过程中第一主成分分量的贡献率对比。

表3 部分贡献率对比

从表3中可以看出，传统计算方式下的贡献率同时序计算方式下的贡献率相近，在个别年份上表现为略高于时序计算下的贡献率，但整体平均贡献率不足于时序计算下的贡献率。运用时序计算方式所得到的贡献率在不同年份上存在差别，大部分结果的贡献率能达到85%以上，部分结果的贡献率不足80%，整体上贡献率均值为85.37%。这说明时序计算方式所得的结果能较好地表征生态环境。

3.3.2 多项式拟合

通过对研究区进行取均值的方法对生态进行评价。用像元均值来代表该年的生态情况，并通过折线图的形式表示连续年份的生态变化情况。

分别采用传统的RSEI计算方式以及RSEI时序计算方式进行对照实验。结果如图6所示，RSEI1代表时序计算方式的结果，RSEI2代表传统计算方式的结果。两者都表现出类似的生态变化走向。通过分析可以得出： 奎屯市的总体生态情况可分为3个阶段。第一阶段是1989—2008年的起伏增长阶段； 第二阶段是2008—2012年的高速增长阶段; 第三阶段是从2012年至今的平缓稳定阶段。奎屯市位于内陆，地势平坦，受自然灾害的影响有限，不易发生生态环境的突变。利用光滑函数拟合的方式对生态进行评价有一定合理性。

图6 RSEI均值时序

通过对折线进行趋势线拟合并比较R2差别，可以明显看出2种计算结果之间的差别。由RSEI折线本身的变化趋势得到，不宜采用对数、指数等函数进行拟合。本文主要采用线性、二次多项式和三次多项式进行拟合对比，结果如表4所示。

表4 趋势线R2比较

从表4中可以看出，2种不同的计算方式在多项式拟合的结果上存在差异。其中，RSEI时序计算方式取得了更好的拟合效果，RSEI传统计算方式则存在更大的扰动。在3种函数的拟合结果上，2种方法的R2都存在高于1%、接近2%的差异。

4 结论

本文以奎屯市为研究区，使用GEE作为主要研究平台，配合以QGIS等工具，进行了RSEI在长时间序列分析中的探索。本文改进了RSEI的计算方式，使之更适用于GEE平台的长时间序列分析。该方式是以增强信息之间可对比性为目的提出的。主要结论如下：

1)对RSEI各个分量的计算方式进行了调整。在湿度指数的计算上，对比了以往研究中出现的缨帽变换参数，并根据目视判读等原则选取更合适的缨帽变换参数； 在温度指数计算方面，对比了2类传统的地表温度反演方法，并与相同时间尺度下的MODIS影像结果进行对比，采用统计精度最高的计算方法。

2)在指数的归一化上采用整体最值的归一化方式，规范了数据标准，增加数据之间可比性。对2种计算方式的最终计算结果进行了第一主成分分量的贡献率对比，得出RSEI时序方法在贡献率方面的优势。并分别采用线性、二次多项式和三次多项式对RSEI长时间序列结果进行函数拟合，验证了RSEI时序计算的优化性。