結(jié)合鄰域信息和結(jié)構(gòu)特征的遙感影像變化檢測

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

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

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

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

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

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

1 本文方法

本文方法主要包括以下3個部分：①獲取鄰域信息：通過提取雙時相影像間的NCI特征，獲得局部鄰域的上下文信息，同時在影像間利用鄰域像素的互相關(guān)性進(jìn)行模板匹配獲得匹配誤差。②結(jié)構(gòu)特征提取：利用方向梯度信息構(gòu)建結(jié)構(gòu)特征描述符，提取雙時相影像的結(jié)構(gòu)特征。③生成檢測結(jié)果：將NCI、匹配誤差、雙時相影像結(jié)構(gòu)特征作為屬性信息進(jìn)行決策樹分類，獲得初始變化檢測結(jié)果，然后以雙時相影像波段差異圖為特征場，初始變化檢測結(jié)果為標(biāo)記場進(jìn)行MRF優(yōu)化，獲得變化檢測二值圖。方法的總體流程如圖1所示。

1.1 獲取鄰域信息

1.1.1 鄰域相關(guān)影像

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

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

(1)

(2)

(3)

(4)

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

1.1.2 匹配誤差

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

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

1.2 結(jié)構(gòu)特征提取

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

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

(5)

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

1.3 生成檢測結(jié)果

1.3.1 決策樹分類

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

1.3.2 MRF優(yōu)化

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

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

2 試驗(yàn)與分析

2.1 試驗(yàn)數(shù)據(jù)

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

第1組數(shù)據(jù)集為墨西哥數(shù)據(jù)集(圖4)，其中，圖4(a)、(b)分別為在2000年4月、2002年5月獲得的墨西哥郊外的兩幅Landsat 7 ETM+4遙感影像，影像大小均為512×512像素，分辨率為30 m。前時相影像顯示了火災(zāi)尚未發(fā)生時的情形，后時相影像可以清楚地看出火災(zāi)發(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)中矩形框所示。參考變化圖通過目視解譯獲得。

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

圖3 基于CFOG的結(jié)構(gòu)特征構(gòu)建過程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àn)設(shè)計(jì)

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

2.3 試驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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