詹天明，宋 博，孫 樂，萬鳴華，楊國為

1.南京審計大學 信息工程學院，南京211815

2.南京審計大學 審計信息工程與技術協(xié)同創(chuàng)新中心，南京211815

3.南京信息工程大學 計算機學院，南京210094

變化檢測是地球科學以及遙感領域研究最多的問題之一。它的任務是利用同一場景的兩個或多個時間點的圖像，以檢測不同時間點內發(fā)生的變化。其應用范圍包括自然災害監(jiān)測、城市和森林的長期跟蹤、農作物的生長狀態(tài)觀測、土地資源開發(fā)利用審計等。變化檢測從任務角度主要分為：判斷是否發(fā)生變化、變化位置定位、變化發(fā)生的種類等。而變化檢測過程主要包括預處理、變化檢測算法和效果評估這三個步驟。其中，變化檢測算法是其核心步驟。目前，常被用于變化檢測的主要有合成孔徑雷達圖像、全色圖像、多光譜圖像以及高光譜圖像。具有高光譜分辨率的高光譜圖像可以提供比其他遙感圖像更豐富的光譜信息，因此高光譜圖像常被用于識別細微變化，反映出不同物體的詳細成分等問題。近年來，基于高光譜圖像的變化檢測被很多學者廣泛研究，但其較低的空間分辨率、光譜維冗余信息以及較少的數(shù)據(jù)集都給高光譜變化檢測方法的發(fā)展帶來了較大影響。

近年來，基于二分類的變化檢測方法較為常見，它的目的是將多時相遙感影像劃分為變化類和非變化類。Malila提出了變化向量分析（change vector analysis，CVA）方法，該方法計算兩個時相的光譜變化向量，再根據(jù)變化幅度和變化方向判定每個像素點是否發(fā)生變化。在此基礎上，Baisantry 等人提出了基于增強主成分分析（principal component analysis，PCA）和反三角函數(shù)閾值的變化矢量分析方法（PCA-CVA）。該方法核心思想是用兩時相圖像的主成分差分和反三角函數(shù)合理選擇閾值，進而判定當前像素點是否發(fā)生變化。Nielsen 等人利用典型性相關分析的降維屬性提出多元變化檢測（multivariate alteration detection，MAD）方法，該方法通過求解多元變化值以判斷各個像素是否在時間維上發(fā)生變化。此后，Nielsen 在多元變化檢測方法基礎上提出了迭代加權多元變化檢測方法（iteratively reweighted MAD，IR-MAD），該方法在迭代步驟中對每個像素賦予的權重進行更新，迭代結束后以權重作為判斷像素是否變化的標準。

基于監(jiān)督的遙感影像的變化檢測方法也被廣泛研究。其中一類方法的主要思想是通過監(jiān)督學習對各時相遙感圖像進行分類，然后根據(jù)每個像素點在不同時間點上的類別信息進行變化檢測。Bovolo 等人提出的變化檢測框架中先對兩時相的遙感影像分別進行逐像素分類，如果一個像素點在兩時相的影像中分類不同，就可以認為這個像素點產生了變化。另一類方法是對差異圖像進行監(jiān)督分類，將差異圖像分為變化和非變化兩類，從而完成變化檢測任務。Nemmour 等人提出了多支持向量機（multiple support vector machine，Multiple SVM）的變化檢測方法，該方法先對變化圖像進行二分類，再用模糊積分和吸引子動力學指標對分類結果進行優(yōu)化。隨著深度學習方法的發(fā)展，基于深度學習的高光譜變化檢測方法也被廣泛研究。Wang 等人提出端到端2D 卷積神經網(wǎng)絡框架的高光譜變化檢測方法。該方法利用兩時相的光譜信息以及解混結果為每個像素點構建一個混合矩陣，再用卷積神經網(wǎng)絡進行分類。

雖然遙感變化檢測方法被廣泛提出，但受限于高光譜圖像的低空間分辨率和高光譜分辨率，直接應用于高光譜變化檢測時其精度都不是太高。而現(xiàn)有基于深度學習方法的高光譜變化檢測方法需要大量的訓練樣本才能保證較高的檢測精度。此外，這些方法需要將三維高光譜圖像轉換為二維矩陣，一定程度上破壞了高光譜圖像空間關聯(lián)性。為解決該問題，將三維高光譜圖像表示為一個張量的形式進行處理，在高光譜圖像降維、解混、分類和融合中引起了廣泛關注。張量模型利用多線性代數(shù)的強大框架，可以聯(lián)合利用光譜和空間信息，有助于提升后續(xù)高光譜圖像處理的性能。

目前，張量建模在高光譜圖像融合以及分類領域已得到較為廣泛的研究，但在高光譜變化檢測問題中利用張量進行建模的方法較少。為充分利用張量模型中空譜聯(lián)合特征以提高高光譜變化檢測精度，本文提出一種基于協(xié)同稀疏與非局部低秩張量的高光譜變化檢測方法。該方法首先求得兩相位的高光譜差分圖像，再根據(jù)差分圖像中圖像塊非局部特點，提取不同的非局部相似張量簇。然后基于協(xié)同稀疏正則化和低秩正則化建立變化檢測模型，利用極小化方法求解得到表示系數(shù)。最后根據(jù)表示系數(shù)求出張量在變化和非變化兩類中的投影殘差，由投影殘差判斷該張量是否發(fā)生變化。該方法主要有如下兩個優(yōu)點：（1）將非局部圖像塊組合在一起形成一系列低階張量，可以深挖高光譜空譜聯(lián)合信息，同時降低維度；（2）協(xié)同稀疏表示分類可以在較少的訓練樣本下保證一定的分類精度，且在高光譜圖像稀疏分解方面比傳統(tǒng)的低階正則化分類器具有更好的性能，從而克服高光譜變化圖像數(shù)據(jù)集不足的問題。在仿真和真實數(shù)據(jù)集中進行對比實驗，結果表明本文方法可以取得較好的高光譜變化檢測精度。

1 高光譜圖像的張量表示

設χ ∈R為-階張量，各個I表示模式，張量中的每個元素表示為x。階張量的F 范數(shù)和范數(shù)定義為：

階張量沿模展開表示為∈R，將還原成張量則可表示為χ=fold()，其中fold是張量沿模展開的逆操作。設∈R表示大小為J×I的矩陣，和χ 的階積運算定義為：

設J∈R表示高光譜圖像，∈R是字典，S ∈R是系數(shù)張量，N∈R是噪聲。它們之間的關系可以表示為：

其中，、表示高光譜圖像的空間維的高和寬，表示高光譜圖像的波段數(shù)，表示字典中的原子個數(shù)。

2 本文方法

本文方法首先將兩時相圖像J和J做差分取絕對值得到J=|J-J|∈R，然后對差異絕對值圖建立字典={,}∈R，其中、分別表示變化和不變區(qū)域的子字典。建立基于協(xié)同稀疏與非局部低秩張量正則化模型，求出J 在字典中的表示系數(shù)，最后根據(jù)其在兩個子空間的投影殘差判斷每個像素是否發(fā)生變化?？傮w流程圖如圖1 所示。

圖1 本文方法總體流程圖Fig.1 General flow chart of proposed method

2.1 協(xié)同稀疏模型

由式（4）高光譜圖像在字典下的建模方式，兩時相圖像差異絕對值圖像可建模為：

其中，S ∈R是未知的表示系數(shù)的張量形式，N∈R是殘差。根據(jù)表示學習理論，通過求解如下優(yōu)化問題即可求得表示系數(shù)S ：

其中，(S)是S 的先驗函數(shù)。在表示學習理論中，將系數(shù)的張量形式乘以字典，必然導致系數(shù)的張量形式在空間維上僅包含少量非零值，但在最近鄰理論中，每個像素點與字典中的少數(shù)原子具有較強的相關性，因此每個像素點在字典下的表示系數(shù)應該具有稀疏性。為此，基于上述先驗，建立如下協(xié)同稀疏模型：

根據(jù)文獻[15]，協(xié)同稀疏正則項||S||定義如下：

2.2 協(xié)同稀疏與非局部低秩張量分解模型

將上述局部張量塊應用于整個圖像域中，表示系數(shù)張量形式S 的核范數(shù)可表示為：

將式（9）與式（6）相結合，得到如下的協(xié)同稀疏與非局部低秩張量模型：

其中，是非負常數(shù)，第一項是求得J 在字典中的一個投影系數(shù)，而第二項則體現(xiàn)了該系數(shù)的協(xié)同稀疏先驗，第三項S 的核范數(shù)即為圖像的非局部低秩約束，體現(xiàn)了該系數(shù)的空間非局部相關性先驗。對該模型極小化即可得到J在字典中最優(yōu)表示系數(shù)。

2.3 優(yōu)化算法

模型（10）無法直接求解，為此采用交替方向乘子法將目標函數(shù)轉化為如下多個子問題分別求解：

其中，U、U、是引入變量。式（11）的拉格朗日方程可表示為：

其中，是拉格朗日懲罰因子，V、V、V、V是拉格朗日乘子，是懲罰函數(shù)，當≥0 時，()=0，當＜0 時，()=+∞。

由交替方向乘子法，上述變量可通過交替迭代方式求解。求解變量S 的目標函數(shù)為：

其迭代解為：

其中，是單位矩陣，是的轉置。

求解變量U、U、U、U的目標函數(shù)為：

其中，函數(shù)表達式為：

函數(shù)表達式為：

(,)=diag(max{()-,0})(·)表示奇異值分解。

拉格朗日乘子V、V、V、V更新步驟為：

本文方法的總體算法步驟如下：

基于協(xié)同稀疏與非局部低秩張量的高光譜變化檢測算法

輸入：第一時相高光譜圖像J，第二時相高光譜圖像J。

3 實驗結果與分析

3.1 數(shù)據(jù)集

實驗采用了3 組數(shù)據(jù)，第一組是Farmland 數(shù)據(jù)集，第二組和第三組是UASFC（Urban area in San Francisco City）數(shù)據(jù)集。

Farmland數(shù)據(jù)集是由EO-1（Earth observing-1）Hyperion設備采集的2006年5月3 日 和2007 年4 月23 日江蘇省鹽城市一個農田的高光譜圖像。其光譜范圍為0.4～2.5 μm，光譜分辨率10 nm 左右，空間分辨率大約為30 m。該圖像共包含155 個波段，空間大小為450×140 像素。兩相位偽彩色圖及真實變化圖如圖2 所示。

圖2 Farmland 數(shù)據(jù)集Fig.2 Farmland dataset

UASFC 數(shù)據(jù)集是由無人機機載傳感器采集的2009年9月18日和2015年5月11日舊金山市市區(qū)的高光譜圖像。拍攝頻率為1.26 GHz，光譜空間分辨率為1.66 m，方位分辨率為1 m，入射角為25°～65°。共兩組數(shù)據(jù)：第一組數(shù)據(jù)（UASFC1）空間大小為200×200 像素，含138 波段；第二組數(shù)據(jù)空間（UASFC2）大小為100×100 像素，含138 波段。兩相位偽彩色圖及真實變化圖如圖3 和圖4 所示。

圖3 Urban area in San Francisco City 第一組數(shù)據(jù)集Fig.3 The first dataset of Urban area in San Francisco City

圖4 Urban area in San Francisco City 第二組數(shù)據(jù)集Fig.4 The second dataset of Urban area in San Francisco City

3.2 模型評價指標

本文使用總體準確度（overall accuracy，OA）和Kappa 系數(shù)兩個指標評估模型的變化檢測效果。為了計算OA 和Kappa 值，需要統(tǒng)計模型對每個像素的預測情況：模型正確分類了變化的像素數(shù)量記為真陽性（true positive，TP）；模型正確分類了未變化的樣本數(shù)量記為真陰性（true negative，TN）；模型錯誤地將未變化的像素分類為變化像素的數(shù)量記為假陽性（false positive，F(xiàn)P）；模型錯誤地將變化的像素分類為未變化像素的數(shù)量記為假陰性（false negative，F(xiàn)N）。記為像素總數(shù)，則有：

用于衡量模型正確的分類像素比例，計算公式為：

Kappa 系數(shù)衡量模型預測結果和地面真值圖的一致性，計算公式為：

其中，p的計算公式為：

3.3 參數(shù)設置及分析

本文的迭代終止條件為=10，若迭代次數(shù)超過30 次尚未達到該條件則終止迭代。本文模型中的兩個參數(shù)和分別控制聯(lián)合稀疏約束項和非局部低秩約束項，其值較低時約束較小，空譜聯(lián)合特征不明顯，檢測結果中噪聲過多會導致精度下降；而當其值較高時約束較大，檢測結果過多地受到空譜聯(lián)合約束的影響，數(shù)據(jù)本身對變化檢測的決策能力下降，變化檢測精度也不高。為此，根據(jù)各項的數(shù)值量級關系取∈[0.000 5,0.003 0],∈[0.001 0,0.005 0] 進行實驗。圖5 顯示了不同的和值下，本文方法的檢測精度。從圖中結果可以看出，本文方法在=0.001 0 和=0.002 5 時取得最優(yōu)的檢測精度，且在∈[0.001 0,0.002 5]，∈[0.001 0,0.005 0]的范圍內都保持一個較高的檢測精度。

圖5 選取不同參數(shù)時三組數(shù)據(jù)集的平均OA 值Fig.5 Average OA values of three datasets with different parameters

對于監(jiān)督學習的方法，訓練樣本個數(shù)對變化檢測精度的影響也是比較敏感的。為此，本文在每個數(shù)據(jù)集中的訓練樣本總數(shù)取1 000 到8 000 進行實驗，其中正負樣本個數(shù)各占50%。圖6 顯示了不同訓練樣本個數(shù)下本文方法的變化檢測精度。從圖中可以看出，本文方法在Farmland 數(shù)據(jù)集中選取2 000 個訓練樣本即可得到較高的檢測精度，在Urban area in San Francisco City 兩組數(shù)據(jù)集中選取4 000 個訓練樣本得到較高的檢測精度。

圖6 三組數(shù)據(jù)集取不同訓練樣本數(shù)量時的OA 值Fig.6 OA values of three datasets with different number of training samples

本文方法需要選擇圖像塊構造非局部張量，為選擇最優(yōu)尺度的圖像塊，分別選取3×3，5×5 和7×7 三種尺度在Farmland、UASFC1 和UASFC2 數(shù)據(jù)集上進行變化檢測。對應的OA 值如表1 所示。從表中可以看出，圖像塊尺度為3×3 大小時，本文方法在3 個數(shù)據(jù)集上都能取得較好的OA 值。這是因為在高光譜變化檢測問題中，變化區(qū)域包含有很多較小且不規(guī)則的圖像塊。若尺度較大，會對這些圖像塊進行平滑，從而影響變化檢測的精度。在本文實驗中，選取3×3 大小的圖像塊構造非局部張量進行變化檢測。

表1 不同尺度圖像塊在三組數(shù)據(jù)集上的OA 值Table 1 OA values of different scales image blocks on three datasets

3.4 對比實驗

為定量分析本文方法的變化檢測效果，選取基于SVM 的變化檢測方法、基于卷積神經網(wǎng)絡（convolutional neural network，CNN）的變化檢測方法、基于解混和CNN 的變化檢測方法（general end-to-end 2-D CNN，GETNET）、基于低秩約束的變化檢測方法（non-local low-rank sparse unmixing，NLLRSU）、基于張量回歸的變化檢測方法（generalized tensor regression，GTR）以及本文方法去除非局部張量約束的變化檢測方法（collaborative sparsity tensor，CST）與本文方法CSNLRT（collaborative sparsity non-local low-rank tensor）進行比較。圖7～圖9 分別顯示了這些方法在Farmland 數(shù)據(jù)集和Urban area in San Francisco City（UASFC）兩組數(shù)據(jù)集上的變化檢測結果。從圖中可以看出，基于SVM 的變化檢測方法和基于CNN 的變化檢測方法在空間上未進行任何約束，檢測結果中噪聲較多。而GETNET 方法、基于低秩空間約束方法NLLRSU 以及本文方法因使用了空譜聯(lián)合的約束，檢測結果具有較強的抗噪性。表2 給出了這些方法在各個數(shù)據(jù)集中的OA 值和Kappa 值，從結果中可以看出，本文方法在3 組數(shù)據(jù)集中都能取得最優(yōu)的OA 值和Kappa 值。從OA 值和Kappa 值中可以看出，本文所提的張量表示方法可以有效地從兩時相高光譜圖像中準確檢測出變化的部分。此外，本文所提CSNLRT 方法比CST 方法在OA 以及Kappa指標上都取得較好的結果，這說明本文所提變化檢測模型中的非局部低秩的約束項在抑制噪聲提高檢測精度上起到了非常好的效果。由于本文方法利用了非局部張量低秩建模以及迭代求解，本文方法的運行時間較長。

圖7 不同方法在Farmland 數(shù)據(jù)集上的檢測結果Fig.7 Detection results of different methods on Farmland dataset

圖8 不同方法在UASFC1 數(shù)據(jù)集上的檢測結果Fig.8 Detection results of different methods on UASFC1 dataset

圖9 不同方法在UASFC2 數(shù)據(jù)集上的檢測結果Fig.9 Detection results of different methods on UASFC2 dataset

表2 不同方法在3 組數(shù)據(jù)集上的實驗結果Table 2 Experimental results of different methods on 3 datasets

4 總結

本文針對高光譜變化檢測問題，提出一種基于協(xié)同稀疏與非局部低秩張量的高光譜圖像變化檢測方法。該方法首先求得前后時間點的高光譜差分圖像，再根據(jù)差分圖像中圖像塊的非局部分布特點，提取不同的非局部相似張量簇。然后基于協(xié)同稀疏正則化和低秩正則化建立變化檢測模型，利用極小化方法對模型求解得到表示系數(shù)。最后根據(jù)表示系數(shù)求得張量在不同類別中的投影殘差，進而判斷該張量是否發(fā)生變化。在Farmland 和UASFC 數(shù)據(jù)集中進行實驗，本文方法的OA 和Kappa 值都優(yōu)于其他方法所得結果，表明本文方法在兩時相高光譜數(shù)據(jù)變化檢測任務中可以取得較好的檢測精度。