葉沅鑫，孫苗苗，王蒙蒙，譚 鑫

1. 西南交通大學地球科學與環(huán)境工程學院，四川 成都611756； 2. 高速鐵路安全運營空間信息技術國家地方聯(lián)合工程實驗室, 四川 成都611756

隨著我國城市化與工業(yè)化進程不斷加快，土地作為社會經(jīng)濟發(fā)展的重要支撐，為人類活動和社會發(fā)展提供了基礎的生產(chǎn)資料。準確快速地獲取土地利用變化信息，對社會發(fā)展、環(huán)境保護、自然資源管理等具有重要意義[1-3]。而常規(guī)的地面調(diào)查手段不僅耗時長、成本高，且部分區(qū)域難以到達，無法人為判斷。近年來，遙感成像技術的快速發(fā)展為人類獲取地面信息提供了新的手段，如利用遙感影像對土地利用變化信息進行檢測，此方法已成為變化檢測領域的研究熱點。

目前，遙感影像變化檢測方法，主要可以分為面向?qū)ο蠓ê突谙裨?。面向?qū)ο蠓╗4-6]是將空間鄰近、光譜相似的對象作為基本處理單元，該方法雖然能夠利用高分辨率影像豐富的空間信息，但在生成對象時的最優(yōu)分割尺度卻難以確定[7]。基于像元法以單個像元為處理單元，主要包括代數(shù)法[8-9]、基于深度學習的方法[10-13]和結合多特征的方法。其中，代數(shù)法雖簡單易行但檢測結果椒鹽現(xiàn)象嚴重。基于深度學習的方法雖然能夠通過深層網(wǎng)絡結構自動、多層次地提取復雜地物的抽象特征，但由于目前公開的變化檢測數(shù)據(jù)集較少，難以獲取大量的訓練樣本，在現(xiàn)階段具有一定的局限性[7]。結合多特征的方法[14]可以綜合利用影像的光譜、紋理等多種信息，便于構建穩(wěn)健和適用的變化檢測模型，因此得到了學者們的廣泛關注。

在結合多特征的變化檢測方法中，鄰域信息是最常用的特征之一，可以將其與灰度信息、空間信息等影像特征進行結合，優(yōu)化檢測結果[15]，還可以與紋理信息結合，提高檢測精度[16]。其中鄰域相關影像(neighborhood correlation image，NCI)[15]是一種能夠反映影像局部鄰域相關性的上下文信息，對影像是否發(fā)生變化十分敏銳。但NCI是通過固定窗口分析獲得的，因此窗口大小會直接影響檢測結果，窗口較大時會削弱影像細節(jié)處的變化，而窗口較小會產(chǎn)生嚴重的椒鹽現(xiàn)象。文獻[17—19]提出了基于自適應延伸技術的變化檢測方法，該類方法依據(jù)同類地物像元間的光譜相似性獲得以像元為中心的自適應區(qū)域，然后基于自適應區(qū)域提取變化信息。雖然能夠有效削弱檢測結果中的椒鹽現(xiàn)象，但由于該類方法通常將自適應區(qū)域的平均灰度值作為中心像元值，因此會在一定程度上影響變化地物邊界的準確性。故僅利用單一的空間鄰域信息難以獲得理想的檢測結果。基于此，本文引入一種反映鄰域相關性的測度(匹配誤差)作為NCI的補充，其基本思想是：在完全配準的雙時相影像中，未變化區(qū)域?qū)木植苦徲蚓哂泻芨叩南嚓P性，且相關性最強的鄰域位置幾乎沒有偏移，或者偏移很小，即匹配誤差很?。欢鴮τ谧兓瘏^(qū)域，由于地物類型發(fā)生了變化，雙時相影像中相關性最強的鄰域會遠離對應鄰域位置，即匹配誤差較大。

此外，考慮到多時相遙感影像中光譜差異通常較大，而幾何結構特征相對于灰度信息更加穩(wěn)定，受光譜差異影響更小。故本文將幾何結構特征屬性融入變化檢測方法中，提出了一種鄰域信息與結構特征相結合的遙感影像變化檢測方法。該方法通過將NCI和匹配誤差相結合獲得更加穩(wěn)健的鄰域信息，利用方向梯度通道信息[20]提取影像結構特征，然后將上述3種屬性特征作為決策樹的分類屬性，進行二值分類，獲得初始變化檢測結果，最后借鑒馬爾可夫隨機場(Markov random field，MRF)模型在表征圖像像素空間關系方面的成功經(jīng)驗[21]，對初始檢測結果進行MRF優(yōu)化，得到具有空間一致性的變化檢測結果。

1 本文方法

本文方法主要包括以下3個部分：①獲取鄰域信息：通過提取雙時相影像間的NCI特征，獲得局部鄰域的上下文信息，同時在影像間利用鄰域像素的互相關性進行模板匹配獲得匹配誤差。②結構特征提?。豪梅较蛱荻刃畔嫿ńY構特征描述符，提取雙時相影像的結構特征。③生成檢測結果：將NCI、匹配誤差、雙時相影像結構特征作為屬性信息進行決策樹分類，獲得初始變化檢測結果，然后以雙時相影像波段差異圖為特征場，初始變化檢測結果為標記場進行MRF優(yōu)化，獲得變化檢測二值圖。方法的總體流程如圖1所示。

1.1 獲取鄰域信息

1.1.1 鄰域相關影像

鄰域相關影像是利用雙時相影像間局部鄰域的光譜屬性得到3個表示上下文信息的特征影像，即相關系數(shù)影像，斜率影像，截距影像。它們用于變化檢測的基本原理為：若雙時相影像對應鄰域范圍內(nèi)的地物沒有發(fā)生變化，則該區(qū)域的光譜值傾向于線性相關，即相關系數(shù)較大、趨近于1，斜率趨近于1，截距趨近于0；若對應鄰域范圍內(nèi)的地物發(fā)生了變化，則該區(qū)域的光譜值傾向于線性無關，相關系數(shù)較小、偏離1，斜率偏離1，截距偏離0。為獲取鄰域相關影像，在雙時相影像對應像元位置開辟一個3×3的動態(tài)窗口，然后采用式(1)—式(4)分別計算對應窗口鄰域的相關系數(shù)r，斜率a，截距b，并將值賦給窗口中心像元位置，構建相關系數(shù)影像，斜率影像，截距影像，統(tǒng)稱NCI[15]。

圖1 本文方法總體流程Fig.1 Flowchart of proposed method

(1)

(2)

(3)

(4)

式中，cov12為雙時相影像對應鄰域窗口內(nèi)所有波段亮度值的協(xié)方差；BVi1、BVi2分別為前、后時相對應鄰域窗口內(nèi)所有波段第i個像元的亮度值；μ1、μ2分別為前、后時相對應鄰域窗口內(nèi)所有波段亮度值的平均值；n為鄰域窗口內(nèi)所有波段的像元總數(shù)；s1、s2分別為前、后時相對應鄰域窗口內(nèi)所有波段亮度值的標準差。

1.1.2 匹配誤差

雖然NCI能度量影像間的變化信息，但檢測結果中易存在大量椒鹽。為充分利用影像的空間信息，本文在計算雙時相影像的相關性時，不僅計算對應鄰域位置(記為中心位置)的相關系數(shù)，還通過模板匹配的方式，利用相關系數(shù)計算局部鄰域的匹配誤差，以此來獲取更多的變化信息。其中，匹配誤差用于變化檢測的基本原理為：對于已經(jīng)配準過的雙時相影像，若對應鄰域未發(fā)生變化，則相關性最強的位置幾乎沒有偏移，匹配誤差較小。若對應鄰域發(fā)生了變化，相關性最強的位置會存在較大的偏移，匹配誤差較大。

為獲得匹配誤差，將前時相影像作為基準影像，在后時相影像上對應像元位置開辟一個9×9大小的搜索窗口，然后在搜索窗口內(nèi)進行模板滑動，計算模板與前時相鄰域的相關系數(shù)，獲得相關系數(shù)矩陣，最后通過計算相關系數(shù)最大的點到矩陣中心的歐氏距離d得到匹配誤差(匹配誤差=d)。圖2以黃色點(黃色十字標識符標識的點)和粉色點(粉色十字標識符標識的點)為例對未變化點和變化點的匹配誤差進行可視化展示。其中，圖2(a)和圖2(b)分別展示了兩點在前、后時相影像中的位置，圖2(c)、圖2(d)分別為黃色點和粉色點的相關系數(shù)矩陣可視化柱狀圖，其中采用紅色十字標識符標識相關系數(shù)矩陣中心位置，紅色柱子標識相關系數(shù)最大的點所在的位置。由圖2可以看出，黃色點為未變化點，矩陣中心的相關系數(shù)較小，難以給出正確的變化信息，但該點的匹配誤差較小(匹配誤差為1)，可以給出正確的變化信息。粉色點為變化點，矩陣中心的相關系數(shù)較大，難以給出正確的變化信息，但該點的匹配誤差較大(匹配誤差約為5.66)，可以給出正確的變化信息。以上初步說明了，相關系數(shù)的大小難以精確地反映影像的變化信息，將匹配誤差引入到變化檢測中，可作為NCI的有效補充。

1.2 結構特征提取

考慮到通常情況下，雙時相影像間光譜和灰度信息差異較大，本文基于方向梯度通道特征(channel features of orientated gradients,CFOG)[20]構建結構特征，以抵抗影像間的輻射差異和灰度差異對檢測結果的影響。整個構建過程主要包括以下3個步驟：①利用多方向梯度信息構建方向梯度通道；②對方向梯度通道進行特征通道卷積形成CFOG；③對CFOG進行通道疊加形成結構特征。其中，方向梯度通道由m層構成，每層的方向梯度通道gθ(θ為該層通道所在的劃分方向)在點(h,l)處的梯度幅值gθ(h,l)是通過點(h,l)處的水平方向梯度幅值gx(h,l)和垂直方向的梯度幅值gy(h,l)根據(jù)式(5)[20]計算得到

gθ(h,l)=abs(cosθ×gx(h,l)+sinθ×gy(h,l))

(5)

在形成方向梯度通道之后，在X和Y方向上采用二維高斯核、梯度方向上采用[1,2,1]T核進行特征通道卷積和歸一化處理形成CFOG，最后將CFOG所有通道對應位置求和取平均構建單通道結構特征，以減少特征維數(shù)。圖3顯示了雙時相影像結構特征的構建過程，從中可以清晰地看出，由于CFOG不受大氣折光和灰度差異的影響，因此基于CFOG的結構特征能夠較好地反映地物的幾何結構信息。

1.3 生成檢測結果

1.3.1 決策樹分類

在結合鄰域信息和結構特征的多特征變化檢測問題中，涉及多個屬性特征。決策樹作為一種非參數(shù)化以及能夠在樣本學習中自動估計屬性信息重要程度的分類器，能較好地適用于結合多特征的變化檢測方法，故本文采用決策樹作為分類器獲得初始變化檢測結果。在進行決策樹分類時，首先依據(jù)參考變化圖中變化像元和未變化像元的比例選取訓練樣本。然后以訓練樣本中像元的屬性特征(即NCI、匹配誤差、前時相結構特征、后時相結構特征)作為輸入，其對應的檢測結果(變化像元記為255，未變化像元記為0)作為輸出，采用自頂向下的遞歸方法，以基尼系數(shù)為標準確定樹節(jié)點的最優(yōu)屬性及最佳分割點，并在結點上進行屬性值的比較，依據(jù)像元對應的不同屬性值判斷從該結點向下的分支，直到某個屬性的子集中所有樣本都屬于同一類別，完成決策樹訓練。最后依據(jù)訓練好的決策樹分類規(guī)則對全部像元進行判斷，獲得初始變化檢測結果。

1.3.2 MRF優(yōu)化

考慮到?jīng)Q策樹分類方法是基于像元的分類方法，其檢測結果過于破碎，存在一定的椒鹽噪聲現(xiàn)象。MRF作為概率論的一個分支理論，可以很好地刻畫影像中鄰域像元屬性間的依賴關系，具有較強的抗噪性，其在圖像分割[22]、變化檢測領域取得了廣泛應用[23-24]，故可將MRF模型用于優(yōu)化決策樹檢測結果。即以雙時相影像波段差異圖為觀測場、決策樹檢測結果為初始標記場，依據(jù)MRF與Gibbs隨機場的等價性將標記場的全局最優(yōu)估計問題轉(zhuǎn)化為能量之和的最小化問題，并通過迭代獲得最終二值變化圖。

具體實現(xiàn)思路為：①參數(shù)設置，即初始化標記場，初始化迭代次數(shù)為1，總迭代次數(shù)mtotal=5(考慮到基于MRF模型的變化檢測方法中存在過度平滑的問題[25]，本文采用較小的迭代次數(shù)對初始檢測結果進行優(yōu)化)；②依據(jù)當前標記場分類情況，計算特征場參數(shù)和標準差；③計算特征場能量和標記場能量，并依據(jù)能量最小化原則更新當前分割結果；④重復步驟②③④直到迭代輪數(shù)達到mtotal，獲得最終二值變化圖。

2 試驗與分析

2.1 試驗數(shù)據(jù)

本文通過選取兩組不同傳感器的數(shù)據(jù)集，檢驗算法的有效性，其中每組試驗區(qū)域均包含兩幅已配準的不同時相的遙感影像和一幅參考變化圖，在參考變化圖中白色部分表示變化區(qū)域，黑色部分表示未變化區(qū)域。

第1組數(shù)據(jù)集為墨西哥數(shù)據(jù)集(圖4)，其中，圖4(a)、(b)分別為在2000年4月、2002年5月獲得的墨西哥郊外的兩幅Landsat 7 ETM+4遙感影像，影像大小均為512×512像素，分辨率為30 m。前時相影像顯示了火災尚未發(fā)生時的情形，后時相影像可以清楚地看出火災發(fā)生的位置及范圍(圖4(b)中顏色較暗的區(qū)域)；參考變化圖通過目視解譯獲得。

第2組數(shù)據(jù)集為印度尼西亞數(shù)據(jù)集(圖5)，該數(shù)據(jù)集是由Quick Bird衛(wèi)星于2004年4月和2005年1月拍攝，影像大小為1500×1500像素，分辨率為2.5 m。該數(shù)據(jù)集反映的是印度尼西亞地區(qū)受海嘯侵襲影響的地表變化情況，此時間段正直海嘯過后，植被、水體等土地覆蓋類型發(fā)生顯著變化，且兩幅影像之間存在大氣折光差異(前時相影像中存在大氣折光，如圖5(a)中矩形框所示。參考變化圖通過目視解譯獲得。

注：黃色點為未變化點，粉色點為變化點圖2 匹配誤差示意Fig.2 Schematic diagram of matching errors

圖3 基于CFOG的結構特征構建過程Fig.3 Construction process of structural features based on CFOG

圖4 墨西哥數(shù)據(jù)集Fig.4 Mexico data set

圖5 印度尼西亞數(shù)據(jù)集Fig.5 Indonesia data set

2.2 試驗設計

為了驗證本文方法的優(yōu)越性，將其與4種基于像元的變化檢測方法進行比較，并在表1中對本文方法和其他4種方法所使用的屬性特征進行了說明。其中，方法I為基于灰度信息的變化向量分析(change vector analysis，CVA)法；方法Ⅱ為將影像灰度信息和NCI作為分類屬性的NCI法[15]；方法Ⅲ為基于自適應上下文信息(adaptive contextual information)的方法[19]，記為ACI法；方法Ⅳ為CVA與局部自相似紋理差分度量(local-similarity-based texture difference measure，LSTDM)相結合的方法[16]，記為CVA+LSTDM法；方法Ⅴ為結合鄰域信息(NCI、匹配誤差)和結構特征的本文方法。

2.3 試驗結果與分析

采用5種變化檢測方法分別對兩組數(shù)據(jù)集進行檢測，檢測結果如圖6、圖7所示，其中所有的圖(a)為參考變化圖，圖(b)為CVA法檢測結果，圖(c)為NCI法檢測結果，圖(d)為ACI法檢測結果，圖(e)為CVA+LSTDM法檢測結果，圖(f)為本文方法檢測結果。其中白色表示變化，黑色表示未變化。檢測精度分別見表2、表3。

表1 本文對比的變化檢測方法

由圖6(b)、圖7(b)可以看出，CVA法可以識別墨西哥數(shù)據(jù)集中的主要變化區(qū)域，但椒鹽現(xiàn)象明顯，虛警率較高。對于印度尼西亞數(shù)據(jù)集，CVA法無法正確識別植被區(qū)域內(nèi)像元間的差異，虛警率也較高。這是因為CVA利用的是像元間的灰度差來探測變化像元，對光譜差異敏感。

由圖6(c)可以看出，NCI對存在光譜差異的地區(qū)靈敏度較高，但存在大量虛檢，無法正確識別變化區(qū)域內(nèi)部的小塊未變化區(qū)域(如圖6(c)中灰色框所示)。由圖7(c)可以看出，NCI可以通過局部相關性分析，較好地克服單一像元的孤立性，正確識別植被區(qū)域內(nèi)單個像元間的光譜差異，但難以正確識別存在大氣折光和水質(zhì)差異的未變化區(qū)域(如圖7(c)中灰色箭頭、灰色框所示)，虛警率仍較高。

圖6 墨西哥數(shù)據(jù)集檢測結果Fig.6 The detection results of Mexico data set

圖7 印度尼西亞數(shù)據(jù)集檢測結果Fig.7 The detection results of Indonesia data set

表2 墨西哥數(shù)據(jù)集5種方法檢測精度

表3 印度尼西亞數(shù)據(jù)集5種方法檢測精度

從ACI法對兩個數(shù)據(jù)集的檢測結果(圖6(d)和圖7(d))可以看出，基于自適應延伸技術的鄰域信息可以有效減弱椒鹽現(xiàn)象，削弱大氣折光和水質(zhì)差異對印度尼西亞數(shù)據(jù)集檢測結果的影響(如圖7(d)中灰色箭頭、灰色框所示)，總分類精度和Kappa系數(shù)較高、虛警率較低，但由于ACI法將延伸區(qū)域的灰度均值作為中心像元值，會影響變化區(qū)域邊界的準確性(如圖6(d)中灰色框所示)，且存在漏檢現(xiàn)象(如圖6(d)和圖7(d)中灰色橢圓框所示)。

圖6(e)、圖7(e)檢測結果表明，CVA和LSTDM相結合的方法受大氣折光、水質(zhì)差異影響較小、椒鹽現(xiàn)象較弱，但該方法對變化信息的靈敏性較弱，無法正確識別光譜差異較弱的變化區(qū)域(如圖6(e)和圖7(e)中灰色橢圓框所示)，漏檢率較高。

圖6(f)、圖7(f)的檢測結果表明，本文方法受大氣折光和水質(zhì)差異影響較小，椒鹽現(xiàn)象不明顯，相比于以上4種方法，具有更高的總分類精度和Kappa系數(shù)。這是因為匹配誤差通過刻畫局部鄰域?qū)袼亻g的偏移距離，提供了一種反映鄰域相關性的測度，對NCI進行了有效補充，并通過引入結構特征，增強算法對影像間光譜差異的穩(wěn)健性，獲得較好的檢測結果。

綜合上述試驗結果可以看出，結合鄰域信息和結構特征的本文方法可以有效檢測出地物的變化情況，降低檢測結果中的椒鹽現(xiàn)象。獲得5種方法中總分類精度最高、Kappa系數(shù)最高，漏檢率與虛警率都較低的檢測結果。

3 結 語

針對利用傳統(tǒng)鄰域信息進行變化檢測時存在椒鹽現(xiàn)象嚴重，漏檢率和虛警率較高等問題，本文提出了一種結合鄰域信息和結構特征的變化檢測方法。該方法構建了一種新的反映鄰域相關性的度量測度(即匹配誤差)，并將其與NCI結合來獲取更加穩(wěn)健的鄰域信息，與此同時，引入結構特征增強算法對影像間光譜差異的穩(wěn)健性，然后通過決策樹分類和MRF優(yōu)化獲得具有空間一致性的檢測結果。通過采用兩組試驗數(shù)據(jù)，與CVA和NCI等4種方法進行對比，證明了本文方法能夠獲得精度更高、虛警率和漏檢率均較低的變化檢測結果。由于在本文方法中匹配誤差的搜索窗口是一個重要的參數(shù)，針對不同復雜度的雙時相影像需要對搜索窗口的大小進行調(diào)整，因此搜索窗口大小的自適應性是未來的研究方向之一。