程 巖,趙化啟,劉曉敏,王宇春,李國晶,袁東華,劉 琳,郭 浩
(1.佳木斯大學信息電子技術學院,佳木斯 154007;2.佳木斯大學材料科學與工程學院,佳木斯 154007)
如何提取同一地域不同時期信息的變化已經成為一項重要的研究課題,地球表面上的生態環境是非靜態的,不同時期觀察同一個對象或者某一種現象鑒別其狀態差別的過程稱為變化檢測[1]。從20世紀60年代開始,人類就開始使用遙感技術對地球表面進行觀測,統計地球表面物體變化,從而完成地表數據的更新。自20世紀80年代開始,美國、日本等國家相繼提出了對地觀測計劃(EOS),使人類對于地球表面的理解達到全新階段[2]。我國從“六五”時期開始也采取了衛星遙感對土地資源進行調查研究[3]。
國內外研究學者在地球表面的土地變化自動化提取方面進行了深入的研究。早在1998年,Ridd等[4]利用1986年和1990年時期的TM影像,使用圖像差分法、KT變換法等進行了土地的變化檢測,并進行了結果精度的評估。Ding等[5]使用差值法來進行探測土地覆蓋的變更情況。此后,Liu等[6]使用人工神經網絡(ANN)的方法來檢測衛星圖像中描述的城市中的變化情況,進而驗證了基于神經網絡的變化檢測,其檢測時間短、準確率高等特點。張華國等[7]經過幾何校正、彩色合成和數據融合等遙感數據處理方法,然后采用監督分類法、閾值法、植被指數法和人機交互法相結合的方法進行了土地覆蓋分類,獲得了南麂列島土地覆蓋的最新信息,其與傳統的土地檢測相比具有準確度高、時效性強等顯著優勢。
從技術流程上遙感變化檢測分為大致4個環節:數據選擇及預處理、差異圖生成、差異圖分析及精度評估。變化檢測的方法可分為以下4類:
(1)代數運算法[8-11]:是一種使用非常廣泛的變化檢測方法,主要由比值法、差值法、植被指數法構成。比值法和差值法是對圖像的灰度值相除或相減得到其差異圖,植被指數法差異圖的計算是由紅光和近紅外波段之間進行比值從而突出植被的信息,然后再對植被的指數差值處理而來。對于以上三種方法產生的差異圖,再通過閾值對生成的差異圖進行二值化,最終得到變化檢測圖。
(2)分類法[12-13]:包括直接分類法和分類后變化法兩類,最常用的是分類后變化法,其過程是對圖像進行分類,然后根據分類的結果找到對應的位置進行比較,最終得到圖像差異圖。直接分類法將多個圖像進行重合疊加,再根據所變化的類型來進行分類,不需要對每個圖像進行單獨的分類操作。
(3)變換法[14-15]:通過對圖像數據信息的變換對相關信息進行抑制對所變化的區域有更好的突出。最常用的變換方法有主成分分析法,其過程是將兩幅圖像數據都進行降維,得到主要的主成分,然后通過差值法來進行運算得到差異圖;另一種常用的方法是多元變化檢測法,其過程使用兩張圖像的線性組合的差值來進行變化檢測,得到差異圖,根據其差異圖再進行閾值的選定與分割,進而得到變化效果圖。
(4)深度學習法[16-19]:深度網絡學習架構相對于傳統的變化檢測方法能夠減少檢測的結果對變化檢測過程中所產生差異圖依賴的影響。其中涉及了卷積神經網絡[20]、深度信念網絡[21]、卷積自編碼器[22]等。
傳統的遙感土地變化檢測依賴于同質數據,由于光學數據的光照條件或濕度和降水等其他因素的影響,圖像可能不均勻。特別是當不同的傳感器出現時,這使得直接比較存在變化檢測不準確的問題,分類后比較法適用于同質圖像和異質圖像來進行變化檢測。這種方法依賴于在進行分類時分類算法的優劣。copula理論轉換[23],使兩幅圖達到相似的特征,需要進行合適的關聯且采用標記像素來進行訓練。核典型相關分析[24]、流形學習[25]等也嘗試應用到異質圖像變化檢測,但是這類方法需要進行人工分析,限制性很大。針對異質圖像的特點,將來源不同的圖像映射到同一個特征空間,這是一個亟待解決的問題。因此,本文提出了基于自編碼器網絡架構的梯度信息與差異圖融合變化檢測方法。
針對于傳統的變化檢測存在以下問題:
(1)在進行遙感圖像變化檢測時,來自于不同傳感器的遙感圖像灰度變化具有非線性關系,直接進行變化檢測存在不準確問題。
(2)原始圖像具有不明顯的輪廓和紋理,直接使用原始圖像的輪廓和紋理會導致變化檢測結果不準確。
(3)傳統的差分運算生成的差異圖存在圖像特征表達不準確,導致差異圖分析過程中檢測不準確問題。
本文對于以上傳統的變化檢測存在的問題提出了以下解決方案:
(1)針對異源圖像灰度變化具有非線性關系,本文基于多損失約束的自編碼器的方法來實現一個域的數據高度非線性轉換到另一個域。
(2)針對圖像具有不明顯的輪廓和紋理的問題,本文提出對原始圖像、形態學梯度圖像、方向梯度圖像進行融合用于自編碼器的輸入,提高了本文方法的性能。
(3)針對差異圖生成過程中傳統差分運算導致圖像特征表達不準確問題,本文提出區域特征融合與加權平均融合算法相結合的方法生成差異圖,有效增強了差異圖的特征表達能力。
針對來自于不同傳感器的遙感圖像在進行變化檢測時灰度變化是非線性關系,輪廓、紋理不明顯及差異圖特征表達不準確等問題,本文提出了一種基于自編碼器網絡架構的梯度信息與差異圖融合變化檢測方法,如圖1所示,給定兩幅圖像,使用自編碼器使圖像轉換到同一域下。其中,自編碼器分為編碼和解碼兩個過程,分別生成兩組編碼圖像,選取兩組最優編碼圖像進行圖像差分,再采取區域特征融合與加權平均融合的方法生成差異圖,最后使用大津法[26]進行差異圖分析。同時對原始圖像、形態學梯度圖像、方向梯度圖像進行融合,用于自編碼器的輸入,有效增強了圖像的輪廓和紋理,最終得到了更好的效果。

圖1 梯度信息與差異圖融合變化檢測框架圖
在數據的預處理階段,本文采用網絡架構以自卷積自編碼器為主體。卷積自編碼器不僅可以提取圖像的有效信息和去除部分的冗余信息,更重要的是可以將兩幅圖像映射到相同的特征空間,進行差異圖的生成,有效地解決異源遙感圖像變化檢測問題,其本質是用輸出圖像來高度還原輸入圖像的一個系統。
圖2中組成的系統是自編碼器網絡的結構圖,來實現圖像域到域之間的映射,自編碼器網絡的訓練過程包含編碼和解碼兩個部分,在進行解碼和編碼的過程中,本文引入了文獻[22]中的重建損失、周期一致性損失、加權翻譯損失、編碼相關損失來對自編碼器進行約束,使圖像在閾間轉換過程中更好地對齊。在編碼階段,對輸入的圖像進行編碼,得到編碼層M;在解碼階段,對編碼層M進行解碼,得到輸入圖像的重構,然后通過多損失約束函數使重構誤差達到最小。

圖2 自編碼器網絡結構圖
自編碼器網絡的編碼和解碼函數如公式(1)所示:
sf為編碼器的激活函數,如公式(2)所示:
基于多損失約束的卷積自編碼器網絡偽代碼如圖3所示。

圖3 多損失約束的卷積自編碼器網絡偽代碼
異源圖像之間各個方面存在著比較大的差異,這些差異使圖像在進行變化檢測過程中存在檢測不準確問題,為使異源圖像的顯著信息在變化檢測中更好地體現,本文在數據預處理階段做出了一些改變。本文將方向梯度與形態學梯度圖像與原圖像放到圖像不同的通道中,來進行圖像的變化檢測。
圖像的方向梯度反映了圖像的紋理、邊緣等顯著的信息。梯度值的灰度值變大,圖像對比度增強,此時圖像的紋理更加清晰,灰度變化更加強烈。圖像函數f(x,y)在Q處的梯度,如公式(3)所示:
圖像形態學梯度[27]在去除細節和噪聲的同時也保留了重要區域的輪廓,圖像的輪廓是對物體形態的有力描述,基于形態學梯度算法如公式(4)所示。其中⊕表示膨脹,Θ表示腐蝕,g(x,y)表示輸入圖像,B(i,j)表示結構元素。
腐蝕操作:是求局部最小值算子的操作,使用結構元素B(i,j)對輸入圖像進行腐蝕,將結構體元素B(i,j)占用面積最小的像素點分配給參考的指定像素,使像素圖像中整體灰度值下降,高亮區域逐漸減少。腐蝕公式如公式(5)所示,結構元素對輸入圖像進行腐蝕,是以結構元素B(i,j)的原點為中心,在輸入圖像的區域進行平移,與之并集的就是腐蝕后的圖像。
膨脹操作:是求局部最大值算子的操作,使用結構元素B(i,j)對輸入圖像進行膨脹,將結構體元素B(i,j)占用面積最大的像素點分配給參考的指定像素,使像素圖像中整體灰度值上升,高亮區域逐漸增加。在輸入圖像的區域進行平移,與之交集的就是膨脹后的圖像。膨脹公式如公式(6)所示:
無論膨脹或者腐蝕操作,其過程使結構體元素B(i,j)在圖像上進行平移計算,從而得出膨脹和腐蝕圖像,本文采用3×3的結構元素,過大結構元素會導致圖像邊緣間產生影響,經過計算得到形態學梯度圖像。
對于差異圖的生成,如圖4所示,SAR圖像與光學圖像經過自編碼器網絡后生成了兩組差分圖像,然后以圖像間的相關系數為閾值進行圖像的融合,相關系數表現了圖像之間的相關程度,如果相關性一般,則選取加權平均融合,如果相關性很大,則選取區域特征融合,更加凸顯了差異圖像的特征表達[28]。

圖4 差異圖生成算法結構圖
相關系數計算如公式(7)所示:
具體融合規則如下,利用公式(7)中的相關性數值K,以融合圖像像素點(m,n)為中心,窗口大小為3×3進行相關性數值K計算,一般情況下K≥0.8時,認為兩個變量之間具有很強的關聯性,因此閾值選取為0.8。
如果相關性數值K<0.8,使用加權平均融合算法,圖像方差越大,表示圖像對應的像素值變化越大,以融合圖像像素點(m,n)為中心,計算窗口大小為3×3的像素塊的方差,根據每個窗口的方差除以兩個圖像窗口方差之和得到圖像融合時的權重,計算得到兩張圖像上像素值的加權平均,循環迭代,進行整張圖的融合。
如果相關性數值K≥0.8,使用區域特征融合算法,以融合圖像像素點(m,n)為中心,窗口大小為3×3計算區域能量,采用窗口能量取小法,取圖像能量較小的一方可以極大程度上減小噪聲及誤判區域對差異圖分析的影響,其定義如公式(10)所示:
比較像素點(i,j)為中心的窗口能量,采用窗口能量取小法,如公式(11)所示:
區域特征與加權平均融合的差異圖生成偽代碼如圖5所示。

圖5 區域特征與加權平均融合的差異圖生成偽代碼
此實驗的實驗環境為Windows10系統;處理器型號Intel i5;處理器主頻2.53 GHz和16 GB的運行內存(RAM);PyCharm Community Edition 2020.3.2 x64;Python3.8.0;Matlab2016。
該數據集為中國曙光村土地利用現狀情況,如圖6所示,其中SAR圖像尺寸大小為593×921×1,其圖像來自Radarsat-2的C波段,拍攝時間于2008年6月;光學圖像尺寸大小為593×921×3,其圖像來自Google Earth,拍攝時間于2012年9月,光學圖像包括RGB三個波段。

圖6 曙光村數據集圖像
變化檢測評估指標,總精確度(Overall Accuracy,OA),總精確度代表正確像素的比例,計算總精確度引入以下四個指標:真陽性(True Positive,TP),假陰性(False Negative,FN),假陽性(False Positive,FP),真陰性(True Negative,TN)。O A計算方式如公式(12)所示:
另外一個評估指標是Ka ppa系數,它反映了變化檢測結果圖與參考圖的一致性,Kappa系數計算方式如公式(13)和公式(14)所示:
其中Nc=TP+FN,Nu=F P+TN。
經過方向梯度公式計算后的各個方向梯度圖像,如圖7所示,(a)是SAR水平方向梯度圖像,(b)是光學水平方向梯度圖像,(c)是SAR垂直方向梯度圖像,(d)是光學垂直方向梯度圖像,各個方向的梯度圖像使圖像在細節處紋理更加顯著。

圖7 方向梯度圖像
經過形態學梯度計算公式得到形態學梯度圖像,如圖8所示,(a)是SAR形態學梯度圖像,(b)是光學形態學梯度圖像,形態學梯度圖像很好地描述出圖像的輪廓信息,使得圖像邊緣輪廓處更加明顯。

圖8 形態學梯度圖像
采用不同方法進行圖像變化檢測,在曙光數據集上檢測結果效果圖如圖9所示,(a)是CVA變化向量檢測方法,是通過兩個時相間圖像的向量大小顯示變化的程度,來確定分隔的閾值,從而確定兩個時相之間變化信息;(b)是DPCA采用主成分分析和混合分類方法來對圖像進行變化檢測;(c)是SCCN采用無監督的深卷積耦合網絡的方式來進行變化檢測;(d)是采用文獻[22]中的方法來進行變化檢測,記為方法HR;(e)是采用本文所提出的只含有梯度信息融合的方法所得到的效果圖,記為方法HRA;(f)是本文提出基于梯度信息與差異圖融合的方法得到的效果圖,記為方法HRB。

圖9 所有方法在曙光數據集上的變化檢測圖
使用以下幾種變化檢測方法與本文提出方法進行實驗對比,HR在圖像邊緣處及更多細節處有多處的噪聲,HRA方法極大地減少了噪聲的影響,但是HRA方法使得圖像上的某些區域亮度太高,造成變化檢測不準確。從表1中分析可以得出,HRA與HR相比,Kap pa系數提高了0.07左右,總體精度OA提高了0.84左右,從數據上可以看出HRA在原有方法HR之上性能有了一定程度上的提高,實驗數據表明方法HRA性能良好。

表1 各種變化檢測方法客觀性能數據比較
CVA變化向量檢測方法效果圖噪聲較多,數據表明實驗效果一般;DPCA的實驗數據表明其總體精度OA良好,但Kappa系數低于HRA的值;SCCN的Ka ppa系數略低于方法HRA,總體精度OA略高于方法HRA,實驗表明其性能良好。
HRA經過區域特征融合與加權平均融合算法得到方法HRB,HRB方法更好地抑制了噪聲的產生,使得變化檢測圖輪廓更加清晰,HRB相對于HRA在Kap pa系數上提高了將近0.05,在精確度OA上提高了2.08,說明融合后差異圖的特征表達更為準確,在進行差異圖分析過程中變化檢測精確度有了一定幅度的提高,如表2所示。

表2 基于區域特征與加權平均融合性能評價數據
針對來自異源圖像輪廓、紋理不明顯及差異圖特征表達不準確導致變化檢測不準確的問題,本文提出了一種基于自編碼器網絡架構的梯度信息與差異圖融合的方法,通過自編碼網絡將異源圖像映射到同一空間來進行圖像差分,使用區域特征融合與加權平均融合算法來進行差分圖像的像素級融合,從而生成差異圖,其中,自編碼器網絡的輸入使用了形態學梯度、方向梯度的融合圖像,差異圖分析采用了大津法。從大量實驗結果平均值可以得出,本文從數值指標上取得了不錯的結果。遙感變化檢測技術在推動著社會的不斷發展,在實際生活中,利用變化監測對土地進行監管對于社會的發展具有深遠的影響。