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

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

引 言

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

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

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

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

1.1 研究區(qū)域

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

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

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

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

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

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

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

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

2 研究方法和原理

2.1 直方圖匹配算法

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

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

Hr(gr)=Ht(gt)

(1)

(2)

(3)

2.2 隨機森林方法

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

3 結(jié)果與討論

3.1 逐年影像分析

通過對1986年—2020年逐年NDVI影像進行篩選, 最終共找出1986、 1987、 1988、 1989、 1991、 1993、 1994、 1995、 1996、 1997、 1999、 2000、 2003、 2005、 2006、 2008、 2009、 2010、 2011、 2017年共計20年的影像存在條帶色差問題。 分別對逐年的Landsat TOA和SR影像數(shù)據(jù)進行直方圖匹配修復(fù), 并對比修復(fù)前后的NDVI均值和標(biāo)準(zhǔ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ù)前后影像的進行分析為例, 1994年Landsat SR和TOA影像修復(fù)前后的NDVI平均值分析結(jié)果見表3。 Landsat TOA和Landsat SR圖像修復(fù)前相比, 影像修復(fù)區(qū)的NDVI平均值分別提高了0.111 2和0.121 9, 提高比例的分別高達32.6%和29.03%。 因目標(biāo)影像面積只占總面積的6.76%, 山西省整個研究區(qū)的NDVI值變化并不明顯, Landsat TOA和Landsat SR圖像NDVI值僅分別增加了0.007 7和0.008 5。 總體來看, 修復(fù)后的目標(biāo)影像更加接近于山西省整體的NDVI值。 圖3為修復(fù)前后的影像, 從圖像的平滑度來看修復(fù)后的圖像整體過度更加平滑, 能更好地反映山西省歸一化植被指數(shù)的真實情況。

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

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

1994年影像修復(fù)前后NDVI剖面的對比分析如圖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ù)前更加符合影像的整體過度。

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

3.2 影像時序分析

為了驗證影像修復(fù)前后對于山西省34年間NDVI的影響, 分別對影像修復(fù)前后的NDVI值進行時序分析, 如圖5所示。 一元線性回歸的趨勢性分析結(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é)果在長時間序列的變化波動性更小, 趨勢更加平滑。 修復(fù)后的Landsat SR影像比Laodsat TOA影像提升幅度更明顯。

圖5 1986年—2020年NDVI的時間序列分析

3.3 多源影像對比分析

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

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

圖7 2003年不同分辨率的MODIS驗證圖像

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

圖8 2000年至2020年的NDVI時間序列分析

4 結(jié) 論

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

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