一種基于分層策略的時空融合模型

張愛竹，王 偉，鄭雄偉，姚延娟，孫根云，辛 蕾，王 寧，胡 光

(1.中國石油大學(華東)海洋與空間信息學院,青島 266580；2.青島海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室,青島 266237；3.中國自然資源航空物探遙感中心，北京 100083；4.環境保護部衛星環境應用中心,北京 100094；5.國家海洋局北海預報中心，青島 266061；6.武漢中測晟圖遙感技術有限公司，武漢 430223)

0 引言

長時間序列的高空間分辨率遙感數據是高精度地表快速監測的基礎，廣泛應用于土地利用制圖和變化檢測、農作物的生長監測、地表溫度監控等領域[1]。然而，遙感數據在空間分辨率和時間分辨率上的相互制約，導致單一傳感器無法同時獲取高空間和高時間分辨率的遙感數據。時空融合技術可較好地解決這一矛盾，在植被動態監測[2]、葉面積指數[3]、城市熱島監測[4]、地表溫度生成[5]、土地覆蓋分類[6]、病毒傳播監測[7]等多個領域得到應用。

在眾多已有的遙感時空融合算法中，簡單高效的線性模型得到了最為廣泛的應用[8-13]。其本質思想是通過構建高低空間分辨率圖像中對應像元反射率之間的線性關系，將低空間分辨率圖像的時間變化信息融入到已有的高空間分辨率圖像中[14]。Gao等[15]最早提出了一種時空自適應反射率融合模型(spatial and temporal adaptive reflectance fusion model,STARFM)，該模型通過融合MODIS圖像和Landsat圖像生成了具有MODIS時間分辨率和Landsat空間分辨率的地表反射率數據。雖然STARFM模型能夠獲取預測數據，但是不適用于反射率變化較大的區域和異質性較高的區域[16]。針對前者，Hilker等[17]提出了一種反射率變化映射的時空自適應算法，該算法首先確定不同時相Landsat的空間變化和MODIS的時間變化，然后選擇地表覆蓋變化與MODIS最接近的Landsat圖像作為數據預測的基準數據，從而提高了Landsat圖像的預測精度。針對高異質性問題，Zhu等[18]提出了增強型時空自適應反射率融合模型(enhanced spatial and temporal adaptive reflectance fusion model,ESTARFM)，假設地物端元在短時間內只發生線性變換，通過引入轉換系數計算出不同端元的反射率變化，提高遙感圖像異質性區域的預測精度。雖然這些改進模型取得了一定的應用效果，但是它們的普遍問題是不能適用于土地覆蓋類型發生突變區域的預測[19]。

近些年來，為了進一步提高時空數據融合的預測精度，研究人員嘗試將短時間內同一研究區的數據變化關系建模為線性變化模型[20]。其本質思想是將基準和預測時刻的兩幅低空間分辨率圖像的線性方程的系數直接應用于高空間分辨率圖像的預測。例如，Cheng等[21]提出了一種基于非局部濾波的時空融合模型，利用兩個回歸系數更準確地描述地表覆蓋變化信息，提高了對復雜變化地物的預測能力；Wang等[22]提出了一種名為Fit-FC的新型時空融合算法，綜合利用回歸模型、空間濾波和殘差補償提高預測精度。但是，如何對土地覆蓋類型發生突變的區域進行數據融合與預測仍然是個難題。

針對時空融合模型在突變區域預測精度不高的問題，本文提出了一種基于分層策略的時空融合模型(hierarchical spatial-temporal fusion model,H-STFM)。通過分別捕捉漸變和突變的土地覆蓋類型變化區域，并設計不同的融合方法，有效提高預測圖像異質性區域的預測精度。該模型由線性回歸模塊和加權濾波模塊組成，分別用于預測物候變化和突變區域。首先判別目標像素(待預測像素)的變化性質，然后對物候變化像素進行線性回歸預測；對突變像素進行加權濾波處理，在超像素鄰域中選擇變化同質性像素作為相似像素進行加權濾波預測；最后將物候變化和突變區域的預測結果利用優化的時間加權函數融合生成最后預測圖像。旨在充分考慮地表覆蓋的突出變化，解決時空融合模型在地表突出變化上預測效果不好的問題，使得預測圖像更加接近真實地表數據。

1 H-STFM模型

H-STFM模型主要包括4個步驟：變化像素的判斷、加權濾波預測、線性回歸預測、時間加權融合，具體流程如圖1所示。為了便于討論，本文將低空間分辨率且重訪周期短的遙感圖像稱為“低分辨率圖像”，將高空間分辨率但重訪周期長的遙感圖像稱為“高分辨率圖像”。

圖1 算法流程圖Fig.1 Algorithm flowchart

1.1 判斷變化像素

H-STFM模型對變化像素進行判斷，需要兩個基準時刻(圖1所示t1和t3時刻)的兩幅高、低分辨率圖像對和預測時刻(圖1所示t2時刻)的低分辨率圖像。首先，對上述5幅遙感圖像進行數據預處理，將3幅低分辨率圖像進行重采樣和剪切處理，使其與高分辨率圖像具有相同的空間分辨率和分布范圍。然后，計算t1和t3時刻低分辨率圖像分別與t2時刻低分辨率圖像的反射率差值并進行閾值判斷，將大于閾值的數據選出作為突變數據，其他作為物候變化數據。以t1與t2時刻低分辨率數據為例，判斷方法如式(1)所示：

(1)

式中：M為已知的低分辨率圖像；(x,y)為目標像素的坐標；B為圖像的波段；σ(M2-M1)為兩期低分辨率圖像反射率差值的標準差；N為土地覆蓋類型的數量。即閾值由反射率差值的標準差與土地覆蓋類型數量決定。基于以上方法，可以分別得到t1與t2時刻以及t3與t2時刻低分辨率圖像的突變數據和物候變化數據。

1.2 基于線性回歸的預測

本文采取兩層處理的策略分別對物候變化數據和突變數據進行預測，即對整體物候變化和局部突變分別進行預測。對于物候變化數據，本文采用全局線性回歸模型，即對兩幅低分辨率圖像使用全局窗口來對每個波段中像素進行線性回歸擬合，如式(2)所示：

M(x,y,t2,B)=a×M(x,y,t1,B)+b+R，

(2)

式中：a和b為線性回歸的兩個系數，可用最小二乘估計法得到[10]；R為系統殘差。

線性回歸模型的理論基礎是在同一時刻的遙感圖像中，低分辨率圖像和高分辨率圖像在相同時間范圍內的變化是一致的。因此，可以直接將基于t1時刻低分辨率圖像構建的回歸模型，應用于t2時刻高分辨率圖像的預測中。進而得到物候變化區域t2時刻的高分辨率圖像，如式(3)所示：

LR(x,y,t2-1,B)=a×L(x,y,t1,B)+b+R，

(3)

式中：LR(x,y,t2-1,B)為基于回歸模型基于t1時刻預測的t2時刻高分辨率圖像；L(x,y,t1,B)為t1時刻的高分辨率圖像。同理，可以基于t3時刻圖像預測得到t2時刻的高分辨率圖像LR(x,y,t2-3,B)。

1.3 基于加權濾波的預測

對于突變數據采用加權濾波的方法進行預測。經典時空融合算法僅基于光譜信息在矩形窗口中選擇相似性像素，不能充分利用高分辨率圖像中豐富的空間特征[23]，不適用于突變數據選擇相似性像素。為充分利用圖像的空-譜信息，本文先對t1時刻的高分辨率圖像利用分形網絡分割算法進行超像素分割[24]，然后將超像素作為鄰域窗口進一步約束相似像素的選擇范圍。此外，由于目標像素是突變數據，即在該時間段內地表覆蓋發生了較大的變化。為了更準確地選擇相似像素，本文提出基于混合像元分解的相似像素選擇方法：先計算t1與t2時刻低分辨率圖像目標像素的光譜反射率差值D，根據高分辨率圖像的標準差和土地覆蓋類別設置閾值Tp，得到閾值范圍為[D-Tp,D+Tp]，如式(4)所示：

Tp=σ(L1)/2×N，

(4)

式中：Tp為判別異常像元的閾值;σ(L1)為t1時刻高分辨率圖像反射率的標準差。

然后計算t2時刻高分辨率圖像目標像素與超像素內其他像素的反射率差值，如果某一像素對應的差值在閾值范圍內，則被選定為備選相似像素；接著利用基于純凈像元指數的端元提取方法[25]與完全約束最小二乘法[26]完成基準時刻高分辨率圖像的混合像元分解，得到其端元豐度圖；最后，根據實驗經驗設置端元的變化閾值，將備選相似像素與目標像素反射率差值大于閾值的像素作為最終的相似像素，相似像素個數標記為N(x,v)。該方法基于地表覆蓋變化信息選擇相似像素，即選擇該像素經過土地覆蓋變化后的光譜相似像素，而非原始像素的相似像素，與傳統方法差異較大。為了避免混淆，本文將相似像素定義為變化同質性像素。

圖2以兩個時刻Landsat與MODIS圖像中某一超像素為例，展示了本文變化同質性像素的選擇方法。如圖2(a)所示，t1到t2時刻的MODIS圖像目標像素的反射率變化為5(反射率從5變化到10)。閾值Tp根據式(4)計算得到，數值為0.5，因此本文選擇變化同質性像素的閾值范圍為[4.5,5.5]。即在t1時刻Landsat的超像素鄰域中，選擇與該目標像素反射率之差在閾值范圍的像素作為變化同質性像素的備選像素，如圖2(a)中標記為紅色的反射率為9.6,9.7,10.4和10.5的像素。圖2(b)為t1時刻的Landsat圖像混合像元分解之后得到的3幅端元豐度圖。由于本小節處理的數據為突變數據，假設目標像素已發生較大地物變化，所以應該選擇在t1時刻與目標像素分度值差異較大的像素。本文根據實驗經驗設置端元變化閾值為0.5，選擇與目標像素差值大于端元變化閾值的備選像素，最后將3幅豐度圖的備選像素取交集作為最終的變化同質性像素，如圖2(b)所示。

(a)基于閾值的變化同質像素初步選擇 (b)基于混合像元解混的變化同質像素約束圖2 基于超像素的變化同質性像素篩選Fig.2 Superpixel-based change pixel filtering

獲取變化同質性像素后，需要據其構建加權濾波的權重函數。為充分利用圖像的空間與光譜信息，本文綜合利用目標像素與變化同質性像素的光譜差異和距離差異構建權重函數。對于坐標為(x,y)的目標像素，其變化同質性像素的坐標可以定義為(xi,yi)，其中i∈[1,2，…，Nx,y]，j∈[1,2，…，Nx,y]。目標像素(x,y)與變化同質性像素(xi,yi)之間的光譜差異Sij表示為：

Sij=|[L(x,y,t1,B)-L(xi,yi,t1,B)]-[M(x,y,t2,B)-M(x,y,t1,B)]|。

(5)

對應的距離權重Dij由變化同質性像素距離目標像素的空間距離定義，公式為：

(6)

式中A為比例系數，調整距離權重與綜合權重的限制常數。

變化同質性像素與目標像素光譜差異越小、空間距離越接近，表示與目標像素的相似程度越高，因此應當被賦予更高的權值。基于此，本文將光譜差異和距離差異進行歸一化，得到每個變化同質性像素對目標像素的貢獻權重Wij，公式為：

(7)

求得歸一化權重后，對于突變數據，利用基準時刻t1的低分辨率圖像計算預測時刻t2的高分辨率圖像，公式為：

(8)

式中LF(x,y,t2-1,B)為基于濾波模型利用t1時刻預測的t2時刻高分辨率圖像。同樣，可以利用上述方法得到基于t3時刻的基準圖像預測得到t2時刻的高分辨率圖像LF(x,y,t2-3,B)。

基于1.2和1.3小節的方法，可以分別針對物候變化區域和突變區域得到基于t1時刻圖像預測的t2時刻圖像，如式(3)與式(8)所示。兩者綜合得到整體的t2時刻預測圖像，公式為：

L(x,y,t2-1,B)=LR(x,y,t2-1,B)+LF(x,y,t2-1,B)。

(9)

同理，可以基于t3時刻的高低分辨率圖像，預測得到t2時刻的高分辨圖像L(x,y,t2-3,B)。

1.4 時間加權

獲得兩幅t2時刻的預測圖像后，對其進行時間加權融合得到最終的預測結果。本文選用余弦相似度作為權重的賦值標準：將2期低分辨率圖像看作向量，計算2個向量的夾角余弦值來評估它們的相似度，即余弦值越大相似性越高，設置的權重越大。余弦相似度不僅可以體現光譜值差異，還可以體現光譜曲線形狀上的差異，計算的相似度量更為可靠。具體的時間權重函數Thk表達式為：

(10)

其中，余弦相似度表達式為：

(11)

式中：k為不同時刻的基準圖像，值為1或3；cos為夾角余弦值；θ1為t1和t2時刻低分辨率像素向量夾角；θ3為t3和t2時刻低分辨率像素向量夾角。基于以上權重，可綜合求得t2時刻的高分辨率圖像L(x,y,t2,B)，公式為：

L(x,y,t2,B)=L(x,y,t2-1,B)×T1+L(x,y,t2-3,B)×T3。

(12)

2 實驗

2.1 實驗數據

本文采用兩個Landsat7 ETM+和MODIS數據集對算法性能進行測試，每個數據集包含3幅Landsat圖像和3幅MODIS圖像。經過投影校正、大氣校正、重采樣和裁剪等預處理，Landsat圖像分辨率為30 m，MODIS圖像分辨率為500 m。本文采用綠、紅、近紅外3個波段進行預測實驗，即Landsat的4，3，2波段和MODIS的3，1，2波段。

數據集1實驗區位于加拿大北部區域(http://ledaps.nascom.nasa.gov/ledaps/tools/StarFM.htm)，如圖3所示。研究區內土地覆蓋類型以森林(云杉、松樹、白楊)為主，輔之以沼澤和稀疏植被(土壤、巖石、燒傷痕)，植被生長季節短，物候變化迅速。

(a)t1時刻MODIS圖像 (b)t2時刻MODIS圖像 (c)t3時刻MODIS圖像

數據集2位于中國內蒙古區域(https://drive.google.com/open?id=1yzw-4TaY6GcLPIRNF BpchETrFKno30he)。研究區內地物目標豐富，包含農田、牧場、山區、城鎮、城市等，如圖4所示。可以看出，數據集2相比于數據集1具有更高的光譜異質性和空間異質性，預測難度更高。

(a)t1時刻MODIS圖像 (b)t2時刻MODIS圖像 (c)t3時刻MODIS圖像

2.2 融合評價指標

本文采用的時空融合評價指標包括結構相似性指數(structural similarity,SSIM)、平均絕對差(average absolute deviation,AAD)、方差誤差(variance of error,VOE)和相對無量綱全局誤差(erreur relative globale adimensionnelle de synthèse,ERGAS)。SSIM表示預測圖像與真實圖像之間的相似度;AAD代表預測圖像與真實圖像之間的偏差;VOE表示預測圖像的誤差波動;ERGAS評估光譜范圍內所有波段的光譜質量。SSIM值越大表示算法性能越好，最佳數值為1；而AAD,VOE,ERGAS值越低表示圖像融合質量越好。對比算法包括STARFM和ESTARFM兩種。

2.3 實驗結果和討論

圖5展示了數據集1的預測圖像。其中，圖5(a)為真實Landsat圖像，圖5(b)—(d)為3種時空融合算法生成的預測圖像。可以看出，3個模型的預測結果整體相似，H-STFM模型獲取的預測圖像與真實圖像整體色調上更接近。而從局部放大圖可以看出，雖然STARFM和ESTARFM模型也能夠捕捉到圖像的整體變化，但是對很多空間細節信息的預測不夠準確，且有較嚴重的顏色失真和色斑，H-STFM模型獲取的預測圖像更好地保留了地物的顏色和空間細節信息。

(a)真實圖像 (b)STARFM (c)ESTARFM (d)H-STFM

圖6是以紅光波段為例展示了數據集1 的圖像反射率預測值與真實值的散點圖，其他數據結果展示如表1，其中加粗體為各指標最佳結果。如圖6所示，相比于兩種對比算法，H-STFM模型得到的散點圖中所有的數據都接近1∶1線，即其能夠較好地捕捉到物候引起的反射率變化，擬合的效果更好。

(a)STARFM (b)ESTARFM (c)H-STFM

表1 數據集1的定量評估結果Tab.1 Quantitative evaluation results of data set 1

從表1中可以看出，在綠波段和紅波段，H-STFM模型得到的預測圖像的偏差AAD和誤差波動VOE更小，預測效果更好。雖然在近紅外波段上H-STFM模型獲得的AAD略高于ESTARFM模型，但是H-STFM模型在各個波段上的SSIM和VOE都是最優的，表明H-STFM模型得到的預測圖像與真實圖像之間整體的相似度最高并且誤差波動最小。H-STFM模型得到的ERGAS指標為0.455 5，明顯優于兩種對比算法，表明H-STFM模型在相關光譜范圍內所有波段的光譜質量最高。

圖7是數據集2的預測結果圖。其中，圖7(a)為真實Landsat圖像，圖7(b)—(d)為3種時空融合算法生成的預測圖像。由圖7可以看出，ESTARFM和H-STFM模型可以得到較STARFM模型質量更高的預測效果；而在黃色圓圈區域，H-STFM模型的預測效果最好，與真實圖像最為接近。在紅框標示的局部區域可以看出，H-STFM模型的預測圖像更好地保留了真實圖像的整體色調，且更好地預測了圖像空間信息，而STARFM的預測圖像出現了較為明顯的斑塊效應，ESTARFM模型在色彩飽和度方面的預測效果較差。結果表明，H-STFM模型在異質性區域同樣能夠獲取較好的預測結果，散點圖和定量評價結果也從定量的角度驗證了這一點。

(a)真實圖像 (b)STARFM (c)ESTARFM (d)H-STFM

圖8同樣以紅光波段為例，展示了數據集2的Landsat圖像反射率預測值與真實值的散點圖，其他數據見表2。與圖6類似，H-STFM模型得到的散點圖中所有的數據較為聚集且接近1∶1線，數據更為集中，即該模型能夠較好地捕捉到突變引起的反射率變化，散點數據擬合的效果最好。

(a)STARFM (b)ESTARFM (c)H-STFM

表2 數據集2的定量評估結果Tab.2 Quantitative evaluation results of data set 2

如表2所示，在紅波段，H-STFM模型得到的預測圖像的AAD和SSIM更小，表示預測效果最好。H-STFM模型在綠波段和近紅外波段上AAD與SSIM表現不穩定，但是在各個波段上的VOE表現穩定且都為最優值，表明通過H-STFM模型獲取的預測圖像與真實圖像之間的相似度最高。與數據集1實驗結果類似，H-STFM的ERGAS指標也優于兩種對比算法，預測圖像的光譜質量最高。

3 結論

針對地表覆蓋的正常物候變化和異常突出變化，設計了一種分層預測的策略，提出H-STFM模型。根據在時間段內地表類型的變化情況，將待預測像素分為物候變化像素和突變像素，對兩種類型的像素分別進行線性回歸和加權濾波預測。同時，對相似像素的選擇方法進行改進，本文利用超像素鄰域窗口和混合像元分解豐度圖進行約束，篩選出更準確的相似像素。最后將時間加權函數進行優化改進，以提高圖像融合精度。在兩個實驗區的實驗結果表明，H-STFM模型可以得到更高的預測精度，更有利于預測土地覆蓋異常變化。主要研究結論如下：

1)利用分層思想，根據像素變化有針對性構建了不同變化數據的預測方法，預測精度明顯提高。

2)引入超像素作為備選鄰域篩選變化同質性像素，并利用豐度圖約束像素選擇，選擇相似像素的準確度明顯提高。

3)選用余弦相似度作為時間加權函數的度量準則，有效提高了預測精度。

本文算法還存在一定的局限性，如將低分辨率圖像建立的回歸模型直接應用到了低分辨率圖像的求解，存在一定的系統誤差。后期研究將進一步考慮不同空間分辨率數據成像系統差異，提高圖像預測的精度。