閆星光, 李 晶, 閆蕭蕭, 馬天躍, 蘇怡婷, 邵嘉豪, 張 瑞

中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院, 北京 100083

引 言

近幾十年, 衛(wèi)星遙感影像成為自然資源調(diào)查、 植被生態(tài)反演、 生態(tài)環(huán)境調(diào)查的主要數(shù)據(jù)來(lái)源[1-3]。 高質(zhì)量的遙感影像應(yīng)具備鑲嵌精度高、 信息豐富和色調(diào)和諧3個(gè)條件; 但由于影像在拍攝過(guò)程中的受到氣候變化、 季節(jié)更替等外部因素以及傳感器自身和衛(wèi)星重返周期的內(nèi)部因素影響, 不可避免地導(dǎo)致影像產(chǎn)生光照分布不均勻的情況, 使得各景影像間的灰度值存在不同程度的差異[4-5]。 特別是在長(zhǎng)時(shí)間序列的中、 大尺度遙感監(jiān)測(cè)中, 通常需要對(duì)比多期遙感影像來(lái)判別地物的變化情況。 以月、 季和年等不同時(shí)間尺度鑲嵌而成的影像數(shù)據(jù), 多會(huì)產(chǎn)生影像灰度值不同的問(wèn)題和相鄰影像條帶的色調(diào)不均勻的現(xiàn)象。

遙感影像勻光處理技術(shù)方法歸結(jié)起來(lái)可以分為三類: 基于照度與反射模型算法、 基于加性噪聲模型算法和基于統(tǒng)計(jì)型的方法[6]。 基于照度與反射模型算法和加性噪聲模型算法方法都是按照數(shù)學(xué)模型來(lái)模擬影像中的亮度分布, 對(duì)單景影像內(nèi)的色調(diào)進(jìn)行不同程度的補(bǔ)償。 例如Mask勻光算法[7]、 Retinex算法[8]和Wallis算法[9]等, 雖然使用Mask勻光算法和Retinex算法能得到較好的效果, 但影像內(nèi)部較暗的區(qū)域容易出現(xiàn)色差和影像邊緣偏色的現(xiàn)象, 從而導(dǎo)致影像灰度值失真的情況; Wallis算法依靠線性關(guān)系進(jìn)行傳遞, 能使原有影像的整體亮度和色調(diào)不發(fā)生較大變化, 但因線性累積誤差的傳遞, 使得離基準(zhǔn)影像較遠(yuǎn)的影像容易出現(xiàn)偏色現(xiàn)象。 基于統(tǒng)計(jì)型的方法是統(tǒng)計(jì)影像的特征和灰度值, 按照一定的規(guī)則使影像間不同位置的色調(diào)和亮度達(dá)到一致, 多以線性變換法、 直方圖匹配法和基于均值和標(biāo)準(zhǔn)偏統(tǒng)計(jì)的方法為主, 其中, 線性變換法處理后的影像雖能同時(shí)考慮區(qū)域范圍內(nèi)多景影像的色調(diào)一致性, 但無(wú)法保持影像非線性的特點(diǎn), 會(huì)導(dǎo)致影像產(chǎn)生局部色調(diào)偏差; 基于直方圖方法的勻色處理, 能很好的解決多景影像灰度值分布不連續(xù)的現(xiàn)象, 適用于影像內(nèi)部特征簡(jiǎn)單的影像; 基于均值和標(biāo)準(zhǔn)偏統(tǒng)計(jì)的影像修復(fù)方法, 僅能通過(guò)單一的統(tǒng)計(jì)值來(lái)進(jìn)行影像修復(fù), 缺乏影像不同特征區(qū)域的精確表達(dá)。 有多景影像灰度值的遙感影像的余光處理過(guò)程中, 主要存在以下三個(gè)方面問(wèn)題: 其一, 中、 大尺度多景影像拼接所造成的條帶色差研究中, 現(xiàn)有的勻光算法在處理后很難實(shí)現(xiàn)影像條帶間亮度和色調(diào)的基本一致; 其二, 在長(zhǎng)時(shí)序遙感監(jiān)測(cè)中, 基于多源遙感的影像修復(fù)規(guī)則不一, 無(wú)法有效保證當(dāng)期影像質(zhì)量的客觀性和長(zhǎng)時(shí)序影像反演的精度; 其三, 現(xiàn)有的研究方法多基于本地ERDAS和GeoDodging等軟件, 在長(zhǎng)時(shí)間序列的中、 大尺度遙感影像勻光處理中, 難以高效和快速的實(shí)現(xiàn)影像修復(fù)的批量化處理[10-12]。 在長(zhǎng)時(shí)間序列和中、 大尺度范圍的遙感影像勻光處理中, 如何高效快速地實(shí)現(xiàn)基于當(dāng)期同源遙感影像的修復(fù)問(wèn)題亟需解決。 因此, 提出了一種基于隨機(jī)森林方法的直方圖匹配方法, 減少影像直方圖形態(tài)特征差異, 從而降低影像失真的現(xiàn)象。

隨著云計(jì)算的興起, 大量的遙感監(jiān)測(cè)和光譜指數(shù)分析都在云平臺(tái)上實(shí)現(xiàn)。 Google Earth Engine(GEE)是衛(wèi)星遙感影像可視化計(jì)算和分析處理的云平臺(tái), 包含了Landsat、 MODIS和Sentinel系列等眾多衛(wèi)星影像數(shù)據(jù), 現(xiàn)有超600個(gè)公共數(shù)據(jù)集[13]。 特別是在長(zhǎng)時(shí)間序列、 大范圍的遙感監(jiān)測(cè)研究中, GEE平臺(tái)憑借其強(qiáng)大的云計(jì)算能力, 極大的提升影像處理時(shí)間和工作效率。 已有諸多學(xué)者在GEE平臺(tái)完成眾多的遙感監(jiān)測(cè)和生態(tài)反演[14-15], 但現(xiàn)有多數(shù)研究中, 都只是將影像進(jìn)行輻射修復(fù)和大氣修復(fù)后直接進(jìn)行植被反演和分析, 忽視由影像拼接造成的條帶色差的問(wèn)題。 針對(duì)上述問(wèn)題, 我們以山西省為研究區(qū), 基于GEE云計(jì)算平臺(tái)采用隨機(jī)森林的直方圖匹配方法, 消除鑲嵌后NDVI影像產(chǎn)生的色調(diào)不均衡現(xiàn)象, 并將1986年—2020年修復(fù)前后的NDVI影像進(jìn)行對(duì)比分析, 旨在為長(zhǎng)時(shí)間序列的遙感影像自動(dòng)化勻色處理尋求一種簡(jiǎn)單、 高效、 適應(yīng)性強(qiáng)的方法。

1 研究區(qū)域及數(shù)據(jù)預(yù)處理

1.1 研究區(qū)域

如圖1所示, 山西省位于中國(guó)華北地區(qū)(圖1), 介于北緯34°34′ N—40°44′ N, 東經(jīng)110°14′ E—114°33′ E之間, 總面積15.67 km2。 山西省呈現(xiàn)“兩山夾一川”的態(tài)勢(shì), 內(nèi)部起伏不平, 山區(qū)面積占總面積的80.1%, 是典型的黃土覆蓋的山地高原, 地勢(shì)東北高西南低。 山西省地跨黃河、 海河兩大水系, 河流屬于自產(chǎn)外流型水系。 黃河流域生態(tài)保護(hù)和高質(zhì)量發(fā)展成為國(guó)家重要戰(zhàn)略, 山西省是黃河流域的重要組成部分, 2021年10月08日中共中央、 國(guó)務(wù)院印發(fā)《黃河流域生態(tài)保護(hù)和高質(zhì)量發(fā)展規(guī)劃綱要》中新的戰(zhàn)略布局當(dāng)中, 在構(gòu)建黃河流域“一軸兩區(qū)五極”的發(fā)展動(dòng)力格局中, 山西省處于“兩區(qū)”即糧食主產(chǎn)區(qū)和能源富集區(qū)的重要地位。

圖1 研究區(qū)示意(a), 山西省位置 (b)

1.2 數(shù)據(jù)處理

Landsat和MODIS系列影像數(shù)據(jù)獲取均來(lái)自于GEE云平臺(tái)公共數(shù)據(jù)集。 圖2為技術(shù)流程, 主要分為數(shù)據(jù)預(yù)處理、 影像匹配和驗(yàn)證分析三個(gè)部分。 數(shù)據(jù)預(yù)處理過(guò)程主要包括: 上傳山西省矢量邊界; 逐年影像時(shí)間篩選(每年1月1日至12月31日)和波段選擇(“Red”、 “NIR”和“pixel_qa”); “pixel_qa”波段去云和NDVI波段計(jì)算; 最后, 按照qualityMosaic函數(shù)進(jìn)行影像拼接和鑲嵌。 影像匹配部分, 首先, 通過(guò)目視解譯的方法, 判斷逐年NDVI影像是否需要進(jìn)行影像勻光處理, 將需影像修復(fù)的部分作為目標(biāo)影像, 參考影像為目標(biāo)影像相鄰的區(qū)域, 目的是讓目標(biāo)影像獲取和相鄰影像的一致的色調(diào); 其次, 分別統(tǒng)計(jì)參考影像和目標(biāo)影像NDVI的DN (digital number)值, 進(jìn)而計(jì)算概率密度函數(shù)和累積分布函數(shù), 按照隨機(jī)森林方法進(jìn)行直方圖匹配, 從而獲得勻光處理目標(biāo)影像。 驗(yàn)證分析分為三個(gè)部分: 第一, 山西省逐年NDVI影像修復(fù)前后對(duì)比; 第二, 1986年—2020年勻光處理前后NDVI影像的時(shí)序分析; 第三, 勻光處理后的NDVI影像和2000年后MODIS影像對(duì)比分析。

圖2 衛(wèi)星遙感影像數(shù)據(jù)處理流程和分析圖

山西省邊界遠(yuǎn)超出單景影像的覆蓋范圍, 因此在影像拼接中將由多幅不同軌跡的影像組成, 目標(biāo)影像選取的原則是按照小于研究區(qū)總面積的50%進(jìn)行修復(fù), 以最大程度的減小勻光處理后影像對(duì)原始數(shù)據(jù)的影響。 經(jīng)統(tǒng)計(jì), 所有年份的目標(biāo)影像面積占比均小于總研究區(qū)的30%, 參考影像色差所選取的區(qū)域均為目標(biāo)影像的銜接條帶。

在影像數(shù)據(jù)的選取上, 根據(jù)現(xiàn)有長(zhǎng)時(shí)間序列的研究, 所使用的影像數(shù)據(jù)多為L(zhǎng)andsat TOA和SR影像數(shù)據(jù), 最終選取Landsat 5 TM、 Landsat 7 ETM+和Landsat 8OLI作為研究數(shù)據(jù), MODIS(MOD13Q1、 MOD13A1和 MOD13A2)數(shù)據(jù)作為2000年修復(fù)后Landsat影像的驗(yàn)證對(duì)比, 具體影像數(shù)據(jù)信息如表1所示。

表1 Landsat和MODIS影像數(shù)據(jù)列表

2 研究方法和原理

2.1 直方圖匹配算法

直方圖匹配算法, 又稱直方圖規(guī)定化匹配方法, 是一種通過(guò)直方圖變化對(duì)圖像進(jìn)行視覺(jué)增強(qiáng)的方法, 多用于圖像增強(qiáng)。 直方圖匹配用于衛(wèi)星遙感影像之間的均衡勻色處理時(shí), 是以一幅影像為參考影像, 將其通過(guò)直方圖映射的方式修復(fù)目標(biāo)影像, 從而達(dá)到多幅影像的色彩一致性[16]。

直方圖匹配主要有3個(gè)步驟: (1)統(tǒng)計(jì)參考影像和目標(biāo)影像逐像素DN值個(gè)數(shù); (2)建立兩者的灰度累積直方圖; (3)通過(guò)一定的映射規(guī)則, 從而使目標(biāo)影像獲取和參考影像相近的累積分布概率[17]。 設(shè)參考影像和目標(biāo)影響的概率密度函數(shù)與對(duì)應(yīng)的累積分布函數(shù)分別為Pr(gr),Pt(gt),Hr(gr),Ht(gt), 首先, 要使目標(biāo)影像的灰度值等于參考影像的灰度值首先滿足式(1)—式(3)

Hr(gr)=Ht(gt)

(1)

(2)

(3)

2.2 隨機(jī)森林方法

隨機(jī)森林是通過(guò)建立多個(gè)獨(dú)立的決策樹, 建立分類樹的集合, 其目的就是通過(guò)大量的基礎(chǔ)樹模型找到最穩(wěn)定可靠的結(jié)果[18]。 無(wú)論所需處理的是離散型數(shù)據(jù), 還是連續(xù)性數(shù)據(jù), 隨機(jī)森林方法都可以有效的處理。 本文所需勻光處理的逐年參考影像位置都是不定的, 影像波段DN值數(shù)據(jù)皆無(wú)固定規(guī)律。 因此, 調(diào)用GEE云平臺(tái)中隨機(jī)森林函數(shù)ee.Classifier.smileRandomForest()參與直方圖匹配, 分別將參考影像和目標(biāo)影像DN值及其概率分布作為訓(xùn)練樣本, 將分類后的目標(biāo)影像的值按照參考影像DN值的概率分布作為直方圖匹配的映射規(guī)則進(jìn)行直方圖匹配。

3 結(jié)果與討論

3.1 逐年影像分析

通過(guò)對(duì)1986年—2020年逐年NDVI影像進(jìn)行篩選, 最終共找出1986、 1987、 1988、 1989、 1991、 1993、 1994、 1995、 1996、 1997、 1999、 2000、 2003、 2005、 2006、 2008、 2009、 2010、 2011、 2017年共計(jì)20年的影像存在條帶色差問(wèn)題。 分別對(duì)逐年的Landsat TOA和SR影像數(shù)據(jù)進(jìn)行直方圖匹配修復(fù), 并對(duì)比修復(fù)前后的NDVI均值和標(biāo)準(zhǔn)差進(jìn)行了分析, 結(jié)果如表2所示。

表2 山西省陸地衛(wèi)星圖像修復(fù)前后的NDVI值和標(biāo)準(zhǔn)偏差

1986年—2020年, 經(jīng)修復(fù)后的影像精度有不同程度的提高, 其中修復(fù)后的Landsat TOA影像NDVI平均值比修復(fù)前提高了0.008 1, 影像的平均標(biāo)準(zhǔn)差下降了0.005 8; 修復(fù)后的Landsat SR影像NDVI平均值比修復(fù)前提高了0.009 2, 影像的平均標(biāo)準(zhǔn)差下降了0.004 4。

以1994年修復(fù)前后影像的進(jìn)行分析為例, 1994年Landsat SR和TOA影像修復(fù)前后的NDVI平均值分析結(jié)果見(jiàn)表3。 Landsat TOA和Landsat SR圖像修復(fù)前相比, 影像修復(fù)區(qū)的NDVI平均值分別提高了0.111 2和0.121 9, 提高比例的分別高達(dá)32.6%和29.03%。 因目標(biāo)影像面積只占總面積的6.76%, 山西省整個(gè)研究區(qū)的NDVI值變化并不明顯, Landsat TOA和Landsat SR圖像NDVI值僅分別增加了0.007 7和0.008 5。 總體來(lái)看, 修復(fù)后的目標(biāo)影像更加接近于山西省整體的NDVI值。 圖3為修復(fù)前后的影像, 從圖像的平滑度來(lái)看修復(fù)后的圖像整體過(guò)度更加平滑, 能更好地反映山西省歸一化植被指數(shù)的真實(shí)情況。

表3 1994年影像修復(fù)前后NDVI值

圖3 1994年修復(fù)前后的NDVI圖像

1994年影像修復(fù)前后NDVI剖面的對(duì)比分析如圖4所示, 修復(fù)后的NDVI影像能平穩(wěn)的反映該區(qū)域的NDVI值。 修復(fù)前Landsat TOA影像的擬合度為0.008 4, 修復(fù)后為0.170 7; 修復(fù)前的Landsat SR影像的擬合僅有0.002 4, 修復(fù)后為0.120 4。 結(jié)果表明, 修復(fù)后的影像擬合度更高, 表明修復(fù)后的影像比修復(fù)前更加符合影像的整體過(guò)度。

圖4 NDVI指數(shù)剖面前后的1994年圖像比較

3.2 影像時(shí)序分析

為了驗(yàn)證影像修復(fù)前后對(duì)于山西省34年間NDVI的影響, 分別對(duì)影像修復(fù)前后的NDVI值進(jìn)行時(shí)序分析, 如圖5所示。 一元線性回歸的趨勢(shì)性分析結(jié)果表明: Landsat TOA影像修復(fù)后的NDVI值斜率為0.006 2小于影像修復(fù)前, 擬合度R2為0.810 6高于修復(fù)前; Landsat SR影像修復(fù)前的NDVI斜率為0.007 1, 而修復(fù)后的斜率為0.006 7, 且擬合優(yōu)度R2值0.836 3大于修復(fù)前的0.829 3, 表明影像修復(fù)后的擬合度更高。 整體上看Landsat影像修復(fù)后的結(jié)果在長(zhǎng)時(shí)間序列的變化波動(dòng)性更小, 趨勢(shì)更加平滑。 修復(fù)后的Landsat SR影像比Laodsat TOA影像提升幅度更明顯。

圖5 1986年—2020年NDVI的時(shí)間序列分析

3.3 多源影像對(duì)比分析

為了驗(yàn)證影像修復(fù)后的準(zhǔn)確性, 將2000年后修復(fù)的Landsat圖像分別與MODIS系列數(shù)據(jù)(250 m、 500 m和1 km)進(jìn)行相關(guān)性分析。 結(jié)果顯示, 2000年、 2003年、 2005年、 2011年和2017年的影像都有明顯改善。 圖6和圖7分別顯示了2003年影像修復(fù)前后的Landsat SR和TOA影像與MODIS(250 m、 500 m和1 km)影像。 直方圖方法修復(fù)后的Landsat影像能更好地反映色彩平衡, 整體視覺(jué)效果更加符合MODIS影像NDVI的分布特征。

圖6 2003年修復(fù)前后的NDVI圖像

圖7 2003年不同分辨率的MODIS驗(yàn)證圖像

2003年修復(fù)后的Landsat影像和MODIS影像的Pearson相關(guān)性分析如表4所示, 修復(fù)后的Landsat SR和TOA與MODIS影像的相關(guān)系數(shù)分別提高了0.049和0.061(p<0.05); 修復(fù)后的Landsat SR與MOD13Q1、 MOD13A1、 MOD13A2影像相關(guān)系數(shù)分別提高了0.050、 0.047和0.049; 修復(fù)后的Landsat TOA與MOD13Q1、 MOD13A1、 MOD13A2影像的相關(guān)系數(shù)分別提高了0.066、 0.060和0.059。 結(jié)果表明, 經(jīng)圖像修復(fù)后的Landsat數(shù)據(jù)與MODIS影像更加一致。

表4 Landsat圖像和MODIS圖像的Pearson相關(guān)分析

修復(fù)前后Landsat和MODIS影像NDVI平均時(shí)間序列變化的對(duì)比結(jié)果見(jiàn)圖8。 與修復(fù)前相比, 修復(fù)后的Landsat影像與MODIS影像的整體趨勢(shì)更加一致, 其中修復(fù)后的Landsat TOA影像的擬合度由0.812 5提高到0.871 1; 修復(fù)后的Landsat SR圖像的擬合度由0.821 8提高到0.853 7。

圖8 2000年至2020年的NDVI時(shí)間序列分析

4 結(jié) 論

目前, 在影像修復(fù)過(guò)程中, 現(xiàn)有的研究多是在研究區(qū)的影像鑲嵌和NDVI計(jì)算之前完成, 這大大限制了影像的處理速度。 本工作針對(duì)歸一化植被指數(shù)NDVI影像拼接后存在的影像帶狀斑塊效應(yīng)和色彩不均勻問(wèn)題, 以山西省作為研究區(qū), 利用GEE平臺(tái)調(diào)用隨機(jī)森林函數(shù)提出一種基于云端快速進(jìn)行直方圖影像勻光處理的方法, 極大的提高了影像修復(fù)的效率。

基于同源影像的直方圖匹配能最大程度地保留當(dāng)期原始影像的DN值和色彩亮度, 且在長(zhǎng)時(shí)序研究中使當(dāng)期影像的NDVI反演結(jié)果更加準(zhǔn)確, 經(jīng)過(guò)該方法處理后影像的色彩一致性較好, 同時(shí)無(wú)需考慮研究區(qū)的地理差異和空間異質(zhì)性。 此外, 通過(guò)對(duì)比1986年—2020年逐年影像修復(fù)前后的結(jié)果, 經(jīng)過(guò)本方法修復(fù)后的影像在長(zhǎng)時(shí)間序列的植被監(jiān)測(cè)過(guò)程中能更精確、 可靠的得出影像的修復(fù)結(jié)果, 有效減少NDVI值在長(zhǎng)時(shí)間序列的突變, 提高長(zhǎng)時(shí)間序列分析的準(zhǔn)確性和穩(wěn)定性。 本方法能有效改善影像條帶色差較大的區(qū)域, 但對(duì)于影像條帶邊界不明顯的區(qū)域識(shí)別仍需提升, 后續(xù)研究的重點(diǎn)將圍繞影像色差邊界的自動(dòng)識(shí)別和修復(fù)展開。