基于GEE的遙感生態(tài)指數(shù)時序計算

駱泓鑒， 明冬萍， 徐錄

(中國地質(zhì)大學(xué)(北京)信息工程學(xué)院，北京 100083)

0 引言

準(zhǔn)確、良好的生態(tài)評價對地區(qū)的發(fā)展、規(guī)劃起重要的支撐作用。通過地理信息的搜集、處理、分析，人們往往能避免主觀判斷的失誤，對生態(tài)環(huán)境做出客觀的評價。

徐涵秋提出了一種基于遙感技術(shù)的生態(tài)評價方法，以自然因子為主構(gòu)建遙感生態(tài)指數(shù)(remote sensing ecological index, RSEI)，從而對城市的生態(tài)狀況進行快速監(jiān)測與評價[1]。不同的學(xué)者用RSEI或改進的RSEI進行生態(tài)評價或其他研究。賈浩巍等[2]以1995年、2005年、2015年為主要年份研究都蘭縣生態(tài)環(huán)境； 李紅星等[3]以2005年、2010年、2015年為研究時間評估武漢市生態(tài)環(huán)境； 繆鑫輝等[4]選取2000年、2009年、2017 年的影像，分析甬江流域生態(tài)環(huán)境變化； 孫從建等[5]定量分析了 2002年、2009年、2017 年3個年份黃土高原的生態(tài)安全，指出治理區(qū)域。非時序的生態(tài)研究都是基于長時間間隔的影像開展的，是跳躍的，只能得出生態(tài)趨勢，無法分析變化細(xì)節(jié)。

谷歌地球引擎(Google Earth Engine, GEE)是用于行星級地理空間分析的云平臺[6]，能提供豐富的開源數(shù)據(jù)和強大的計算服務(wù)。GEE使得長時間序列的龐大計算成為可能，使得RSEI時序計算變得方便。已有學(xué)者基于GEE進行RSEI的連續(xù)分析，并取得相應(yīng)的成果[7-9]。但這些研究在計算方式上仍采用傳統(tǒng)的RSEI。傳統(tǒng)的RSEI算法在分量的計算方式上會隨傳感器的變化而變化，加上不同傳感器之間的成像差異，得到的RSEI計算結(jié)果差異較大，不利于長時間跨度的連續(xù)分析。綜上所述，尚未有RSEI的時序計算研究。本文將探尋一種基于GEE的RSEI時序計算方式，在分量計算及其歸一化方面對RESI算法進行改進，并進行分析驗證。

1 研究區(qū)概況和數(shù)據(jù)預(yù)處理

1.1 研究區(qū)概況

新疆維吾爾自治區(qū)奎屯市坐落在天山北麓中段，準(zhǔn)噶爾盆地西南部，高程在450～530 m之間，地勢南高北低，西高東低，無明顯地形起伏變化，不必考慮地形對研究的影響。奎屯市境內(nèi)為大陸性氣候，日照長，雨少干熱，較易獲取長時間無云或者少云的遙感影像。奎屯市近30 a來生態(tài)、發(fā)展都有巨大變化，利用遙感影像分析生態(tài)變化較為容易。

1.2 數(shù)據(jù)源與GEE數(shù)據(jù)預(yù)處理

研究數(shù)據(jù)采用Landsat的地表反射率產(chǎn)品和MODIS的地表溫度數(shù)據(jù)產(chǎn)品(MOD11A1)。根據(jù)研究的需要，選取1989—2019年間每年6月1日—10月30日期間的影像進行計算分析。計算指數(shù)之前先要獲取去云影像。GEE提供了針對Landsat 地表反射率產(chǎn)品的去云算法。但在實驗中發(fā)現(xiàn)，上述算法在處理異物同譜物體時效果欠佳，容易導(dǎo)致像元丟失。本文提出了減云拼接的去云方法，其流程如圖1所示，其中CLOUD-COVER代表云覆蓋程度。完善之后的去云方法能補全缺失影像，如圖2所示。

圖1 減云拼接方法

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

2 RSEI時序計算

RSEI的核心思想是通過空間降維的方式，從一系列能在一定程度上反映生態(tài)環(huán)境質(zhì)量的分量中，提取出最能表達(dá)當(dāng)?shù)厣鷳B(tài)環(huán)境的綜合結(jié)果。RSEI的基本計算流程如圖3所示。在眾多的自然因素以及人文因素中，綠度、濕度、干度、溫度是與人類活動密切相關(guān)的因素，常用于各類生態(tài)評價[1]。運用函數(shù)的方法對4大分量進行降維壓縮，最終得到表征生態(tài)環(huán)境的RSEI。

圖3 RSEI生產(chǎn)流程

2.1 RSEI分量計算

1)綠度計算。綠度的計算通常考慮植被指數(shù)。歸一化差值植被指數(shù)(normalized difference vegetation index, NDVI)在生產(chǎn)中得到了最廣泛的應(yīng)用[10]，能在一定程度上表征生態(tài)環(huán)境。NDVI的計算公式為：

，

(1)

式中ρred和ρnir分別為紅光、近紅外波段地表反射率。

2)濕度計算。本文使用纓帽變化中的濕度分量進行濕度(Wetness)的計算。不同傳感器之間的纓帽變化參數(shù)存在較大的差異[11-15]。在時序計算過程中，為了消除不同纓帽變化參數(shù)在計算過程中帶來的結(jié)果差異，本文選擇文獻[15]中的纓帽變化濕度參數(shù)作為濕度計算的模型。其計算公式為：

，

(2)

式中ρblue，ρgreen，ρswir1，ρswir2分別為藍(lán)光、綠光、第1短波紅外、第2短波紅外波段的地表反射率。

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

NDBSI=γIBI+(1-γ)BSI

，

(3)

(4)

，

(5)

式中:γ為參數(shù)系數(shù)，通常取值0～1.0，這里根據(jù)前人經(jīng)驗取γ= 0.5[1]；ρmid為中紅外波段地表反射率，但在Landsat中無法找到對應(yīng)的波段，通常采用第1短波紅外波段反射率(ρswir1)來代替。

4)溫度計算。溫度指的是地表溫度(land surface temperature, LST)。它可以通過遙感的方式，用獲取到的熱紅外波段反演得到。其原理是基于普朗克定律量化所構(gòu)成的熱輻射傳輸方程。

前人提出不同的溫度反演方法，有單窗算法、單通道算法和劈窗算法等。本文的實驗數(shù)據(jù)是處理之后的拼接影像，不同像元來自不同時間段，而輻射傳輸模型、單窗算法和單通道算法都需要知道具體時間下的大氣參數(shù)數(shù)據(jù)，不適合本次實驗。

本文采用亮溫轉(zhuǎn)換的方式來進行溫度的反演，其公式為：

，

(6)

式中： C2為普朗克函數(shù)常量，14 387.79 μm·K；B為DN值經(jīng)過輻射定標(biāo)和轉(zhuǎn)換后的亮度溫度；ε為計算得到的比輻射率；λ為波段的有效波長，不同傳感器的有效波長不同[19]。

2.2 整體歸一化

本次實驗是基于不同傳感器、不同拍攝時間的拼接影像開展的。不同傳感器結(jié)果之間存在較大差異，如圖4所示。常見的歸一化方法容易造成標(biāo)準(zhǔn)的不統(tǒng)一，不利于后續(xù)計算分析。

(a) 2016年Landsat7影像濕度歸一化結(jié)果 (b) 2016年Landsat8影像濕度歸一化結(jié)果

為解決上述問題，本文用不同年份計算結(jié)果的整體最值進行歸一化，為不同傳感器的計算結(jié)果制定相同的標(biāo)準(zhǔn)，方便進行對比分析。整體最值歸一化的公式為：

，

(7)

式中：Ni為歸一化后結(jié)果值；Vi為原始像元值；VW-max和VW-min分別為所有計算結(jié)果中的最大值和最小值。

2種不同的歸一化方式的差別如圖5所示。

(a) 單一影像歸一化 (b) 整體最值歸一化

2.3 RSEI最終計算

綠度、濕度、干度、溫度均能在一定程度上表征生態(tài)優(yōu)劣。主成分分析是通過在多維空間中變換坐標(biāo)方向的一種識別數(shù)據(jù)模式的方法，是一種突出數(shù)據(jù)相似性和差異性的方式來表達(dá)數(shù)據(jù)的方式[20]。通過主成分分析，能夠綜合考慮4大分量的作用，構(gòu)建一個表征生態(tài)的指數(shù)，用于生態(tài)評價。RSEI計算公式為：

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

(8)

。

(9)

式(8)表示對歸一化所得的4項指數(shù)進行主成分分析，并取第一主成分分量作為RSEI0。式(9)表示將主成分分析所得的RSEI0進行歸一化處理，得到最終的RSEI。計算結(jié)果顯示，第一主成分分量的平均信息包含量超過了85.37%，能夠較好地代替4大分量，實現(xiàn)數(shù)據(jù)壓縮。

3 RSEI時序分析實驗

不同于傳統(tǒng)的RSEI計算方式，RSEI時序計算在溫度、濕度分量計算以及歸一化方式上都有變化。本節(jié)將通過對比實驗進行方法驗證。

3.1 濕度分量驗證

由于計算濕度分量采用的參數(shù)不統(tǒng)一，所以本文對文獻[12—15]的4種纓帽參數(shù)的計算方式進行了比較，使用2001年和2002年的Landsat7影像，2010年和2011年的Landsat5影像，以及2017年和2018年的Landsat8影像進行測試。

對比表1發(fā)現(xiàn)，文獻[13]和文獻[14]中的濕度計算結(jié)果在上下兼容性上略顯不足，不適用于長時間序列分析。根據(jù)地物的反射波譜可知，含水地域的地表反射率普遍較低，在真彩色影像上表現(xiàn)為深色、暗色調(diào)。以此為依據(jù)對影像進行目視判讀，觀察真彩色影像中的明亮區(qū)域，并與濕度影像中的對應(yīng)區(qū)域進行比較。發(fā)現(xiàn)文獻[12]在一些地區(qū)的結(jié)果與目視判讀不相符，而本文選擇的參數(shù)的計算結(jié)果更符合目視判讀。綜上，本文所選的文獻[15]中的濕度參數(shù)在傳感器兼容性以及合理性上的表現(xiàn)更據(jù)優(yōu)勢，更適合時序分析。

表1 部分年份濕度分量結(jié)果對比

3.2 溫度分量驗證

選擇2007—2011年Landsat5和Landsat7影像及2013—2018年Landsat8影像進行測試，在溫度計算上比較文獻[21]中的模型和亮溫轉(zhuǎn)換模型，在比輻射率計算上比較NDVI閾值法、文獻[22]和文獻[23]中的方法[22-25]。計算結(jié)果通過散點圖的方式和相同時期下的MOD11A1數(shù)據(jù)進行對比。

表2 基于2種LST模型的3種比輻射率算法測試結(jié)果

對比表2結(jié)果，優(yōu)先考慮方差較小的情況，并協(xié)同比較平均差值和中位差值，驗證得到本文選擇的亮溫轉(zhuǎn)換模型和NDVI閾值法作為LST的反演方法是最為合適的。

3.3 時序分析驗證

3.3.1 貢獻率驗證

用主成分分析的方法進行數(shù)據(jù)壓縮，壓縮后的第一主成分分量能在最大程度上代表原始數(shù)據(jù)的主要信息，從而代替生態(tài)分量表征生態(tài)環(huán)境。貢獻率越高，則第一主成分分量的代表能力越強。表3是2種不同RSEI計算過程中第一主成分分量的貢獻率對比。

表3 部分貢獻率對比

從表3中可以看出，傳統(tǒng)計算方式下的貢獻率同時序計算方式下的貢獻率相近，在個別年份上表現(xiàn)為略高于時序計算下的貢獻率，但整體平均貢獻率不足于時序計算下的貢獻率。運用時序計算方式所得到的貢獻率在不同年份上存在差別，大部分結(jié)果的貢獻率能達(dá)到85%以上，部分結(jié)果的貢獻率不足80%，整體上貢獻率均值為85.37%。這說明時序計算方式所得的結(jié)果能較好地表征生態(tài)環(huán)境。

3.3.2 多項式擬合

通過對研究區(qū)進行取均值的方法對生態(tài)進行評價。用像元均值來代表該年的生態(tài)情況，并通過折線圖的形式表示連續(xù)年份的生態(tài)變化情況。

分別采用傳統(tǒng)的RSEI計算方式以及RSEI時序計算方式進行對照實驗。結(jié)果如圖6所示，RSEI1代表時序計算方式的結(jié)果，RSEI2代表傳統(tǒng)計算方式的結(jié)果。兩者都表現(xiàn)出類似的生態(tài)變化走向。通過分析可以得出： 奎屯市的總體生態(tài)情況可分為3個階段。第一階段是1989—2008年的起伏增長階段； 第二階段是2008—2012年的高速增長階段; 第三階段是從2012年至今的平緩穩(wěn)定階段。奎屯市位于內(nèi)陸，地勢平坦，受自然災(zāi)害的影響有限，不易發(fā)生生態(tài)環(huán)境的突變。利用光滑函數(shù)擬合的方式對生態(tài)進行評價有一定合理性。

圖6 RSEI均值時序

通過對折線進行趨勢線擬合并比較R2差別，可以明顯看出2種計算結(jié)果之間的差別。由RSEI折線本身的變化趨勢得到，不宜采用對數(shù)、指數(shù)等函數(shù)進行擬合。本文主要采用線性、二次多項式和三次多項式進行擬合對比，結(jié)果如表4所示。

表4 趨勢線R2比較

從表4中可以看出，2種不同的計算方式在多項式擬合的結(jié)果上存在差異。其中，RSEI時序計算方式取得了更好的擬合效果，RSEI傳統(tǒng)計算方式則存在更大的擾動。在3種函數(shù)的擬合結(jié)果上，2種方法的R2都存在高于1%、接近2%的差異。

4 結(jié)論

本文以奎屯市為研究區(qū)，使用GEE作為主要研究平臺，配合以QGIS等工具，進行了RSEI在長時間序列分析中的探索。本文改進了RSEI的計算方式，使之更適用于GEE平臺的長時間序列分析。該方式是以增強信息之間可對比性為目的提出的。主要結(jié)論如下：

1)對RSEI各個分量的計算方式進行了調(diào)整。在濕度指數(shù)的計算上，對比了以往研究中出現(xiàn)的纓帽變換參數(shù)，并根據(jù)目視判讀等原則選取更合適的纓帽變換參數(shù)； 在溫度指數(shù)計算方面，對比了2類傳統(tǒng)的地表溫度反演方法，并與相同時間尺度下的MODIS影像結(jié)果進行對比，采用統(tǒng)計精度最高的計算方法。

2)在指數(shù)的歸一化上采用整體最值的歸一化方式，規(guī)范了數(shù)據(jù)標(biāo)準(zhǔn)，增加數(shù)據(jù)之間可比性。對2種計算方式的最終計算結(jié)果進行了第一主成分分量的貢獻率對比，得出RSEI時序方法在貢獻率方面的優(yōu)勢。并分別采用線性、二次多項式和三次多項式對RSEI長時間序列結(jié)果進行函數(shù)擬合，驗證了RSEI時序計算的優(yōu)化性。