基于變權(quán)馬爾可夫隨機(jī)場(chǎng)的遙感圖像變化檢測(cè)

高 鵬 翔

(中國(guó)工程物理研究院計(jì)算機(jī)應(yīng)用研究所 四川 綿陽(yáng) 621900)

0 引 言

隨著人類科技水平的迅猛發(fā)展，人類在地球上的活動(dòng)范圍愈發(fā)廣泛。因此及時(shí)、精確地開展對(duì)地觀測(cè)來(lái)實(shí)時(shí)監(jiān)測(cè)地表的覆蓋變化情況，對(duì)人類的發(fā)展具有重要意義[1]。當(dāng)代航空、航天和傳感器等技術(shù)的進(jìn)步，促使了遙感圖像獲取手段的日益成熟[2]。當(dāng)下采用遙感數(shù)據(jù)進(jìn)行變化檢測(cè)已是監(jiān)測(cè)地表覆蓋變化最精確、高效的技術(shù)手段，亦是遙感研究的熱點(diǎn)。遙感圖像變化檢測(cè)是指通過(guò)分析一對(duì)或者一系列在相同地區(qū)不同時(shí)間點(diǎn)獲取的遙感圖像，來(lái)確定地表是否發(fā)生改變的技術(shù)[3]。

變化檢測(cè)已在城市規(guī)劃與建設(shè)農(nóng)業(yè)、牧業(yè)、林業(yè)、工業(yè)、軍事等領(lǐng)域內(nèi)得到了普遍地運(yùn)用，如土地利用/土地覆蓋情況監(jiān)測(cè)、地理數(shù)據(jù)庫(kù)更新、災(zāi)害評(píng)估、森林植被覆蓋監(jiān)測(cè)、土地沙漠化監(jiān)測(cè)、水系變化監(jiān)測(cè)、城市建設(shè)規(guī)劃、軍事目標(biāo)監(jiān)測(cè)等[4-8]。

馬爾可夫隨機(jī)場(chǎng)模型(MRF)由于其完備的數(shù)學(xué)理論，以及能夠充分利用圖像中的空間上下文信息[9]，在遙感圖像變化檢測(cè)上取得了可喜的成就。Bruzzone等[10]首先使用MRF進(jìn)行遙感圖像變化檢測(cè)，對(duì)差值圖像使用EM算法，分別估計(jì)不變類和變化類概率密度，之后引入MRF，通過(guò)條件迭代式(ICM)算法[11]求解。Benedek等[12]采用條件混合馬爾可夫模型來(lái)進(jìn)行變化檢測(cè)，這是一種多層MRF模型，結(jié)合了混合馬爾可夫模型和信號(hào)的條件獨(dú)立隨機(jī)場(chǎng)。Chen等[13]提出了一種上下文敏感的弱監(jiān)督變化檢測(cè)技術(shù)，該方法通過(guò)MRF模型分析概率高斯過(guò)程分類器的后驗(yàn)概率實(shí)現(xiàn)。Chen等[14]針對(duì)傳統(tǒng)的MRF方法在不連續(xù)區(qū)域(如邊界、山脊、山谷等)經(jīng)常得到錯(cuò)誤的變化檢測(cè)結(jié)果，提出了一種無(wú)監(jiān)督變化檢測(cè)算法。Gong等[9]提出了一種基于模糊聚類和改進(jìn)MRF能量函數(shù)的SAR變化檢測(cè)算法，該方法在計(jì)算能量函數(shù)時(shí)同時(shí)考慮了像素點(diǎn)的隸屬度和同種類別的鄰域像素?cái)?shù)目改變。Subudhi等[15]提出了一種基于模糊聚類和MRF的變化檢測(cè)算法，使用模糊吉布斯馬爾可夫隨機(jī)場(chǎng)(GMRF)來(lái)對(duì)差值圖像的空間灰度特性建模，采用基于可變鄰域搜索的全局收斂準(zhǔn)則來(lái)迭代估計(jì)模糊GMRF模型的參數(shù)，避免了傳統(tǒng)算法容易得到局部最優(yōu)值問(wèn)題。范奎奎等[16]將雙樹復(fù)小波變換和MRF相結(jié)合，解決了變化檢測(cè)中多尺度圖像降噪時(shí)丟失高頻信息與單像素孤立問(wèn)題。

然而，傳統(tǒng)的基于馬爾可夫隨機(jī)場(chǎng)的變化檢測(cè)算法，對(duì)于中心像素和鄰域像素使用固定的權(quán)重，會(huì)造成空間鄰域信息的過(guò)度使用[14-15]，使得傳統(tǒng)的標(biāo)記場(chǎng)不能精確地確定鄰域像素間的空間關(guān)系，導(dǎo)致變化檢測(cè)結(jié)果過(guò)度平滑，影響變化檢測(cè)的精度。

針對(duì)上述問(wèn)題，本文提出一種基于變權(quán)馬爾可夫隨機(jī)場(chǎng)的遙感圖像變化檢測(cè)算法。首先通過(guò)兩時(shí)圖像生成差值圖像；隨后使用模糊C均值聚類算法對(duì)差值圖像上的像素點(diǎn)進(jìn)行分類；接著使用一個(gè)改進(jìn)的MRF模型來(lái)對(duì)差值圖像建模；最后使用ICM算法求解能量函數(shù)最小化得到變化檢測(cè)結(jié)果。

1 傳統(tǒng)的基于MRF的變化檢測(cè)

定義大小為M×N的兩時(shí)差值圖像為XD，其滿足XD={x1,x2，…，xn}，n=M×N。變化類別與未變化類別為ω={ωch,ωun}，Y為類別標(biāo)記場(chǎng)，對(duì)于二值變化檢測(cè)問(wèn)題，將圖像上所有像素點(diǎn)分為變化與未變化兩個(gè)類別。根據(jù)MRF理論有：

Y=arg max{P(Y)P(XD|Y)}

(1)

式中：P(Y)為類別標(biāo)記的先驗(yàn)概率；P(XD|Y)是數(shù)據(jù)集中像素點(diǎn)的條件概率密度。

根據(jù)Hammersley-Clifford定理，可知吉布斯隨機(jī)場(chǎng)(GRF)與馬爾可夫隨機(jī)場(chǎng)(MRF)等價(jià)，則式(1)值最大時(shí)等同于MRF能量函數(shù)最小。對(duì)于每個(gè)像素點(diǎn)xi，其能量函數(shù)UMRF(xi)定義為兩部分組成：

UMRF(xi)=Uspectral(xi)+Uspatial(xi)

(2)

式中：Uspectral(xi)是光譜能量函數(shù)，描述了像素點(diǎn)xi的光譜特性，是該像素點(diǎn)的灰度統(tǒng)計(jì)描述；Uspatial(xi)是空間能量函數(shù)，描述了像素點(diǎn)xi與其所在鄰域內(nèi)其余像素點(diǎn)之間的類別依賴。

通常假定兩時(shí)差值圖像的灰度統(tǒng)計(jì)特性滿足高斯分布，于是采用高斯分布來(lái)描述光譜能量函數(shù)Uspectral(xi)：

(3)

對(duì)于空間能量函數(shù)Uspatial(xi)，其根據(jù)中心像素點(diǎn)xi的標(biāo)記與xi的鄰域像素的標(biāo)記情況來(lái)確定。且Uspatial(xi)滿足：

(4)

(5)

式中：β為懲罰系數(shù)，通常β>0，其值由用戶設(shè)定，通過(guò)取不同的β值可以控制鄰域像素對(duì)中心像素的影響，調(diào)整空間上下文信息對(duì)變化檢測(cè)過(guò)程的影響；Ni是像素點(diǎn)xi的鄰域，通常采用二階鄰域，即對(duì)于每一個(gè)中心像素點(diǎn)周圍有8個(gè)像素點(diǎn)；l(xi)和l(xj)是像素點(diǎn)類別標(biāo)記，采用Potts模型來(lái)定義中心像素點(diǎn)與其鄰域像素點(diǎn)類別標(biāo)記關(guān)系I(l(xi),l(xj))；對(duì)于二值指示函數(shù)I(l(xi),l(xj))，當(dāng)中心像素點(diǎn)標(biāo)記與其鄰域像素的標(biāo)記相同時(shí)其值取為-1，否則取為0。

變化檢測(cè)可視為典型的組合優(yōu)化問(wèn)題，其結(jié)果可通過(guò)迭代搜索算法所獲取的全局或局部最優(yōu)解來(lái)得到。使用ICM算法便可以獲得式(2)的最優(yōu)值。

2 基于變權(quán)MRF的變化檢測(cè)

2.1 變權(quán)懲罰系數(shù)

在傳統(tǒng)的基于MRF模型的變化檢測(cè)算法中，Potts模型使用一個(gè)固定的懲罰系數(shù)β來(lái)計(jì)算空間能量函數(shù)。對(duì)于差值圖像中所有的像素點(diǎn)來(lái)說(shuō)，β的取值是相同的，這沒有考慮差值圖像中像素灰度值的分布特性，易造成差值圖像中灰度值極大或者極小的像素點(diǎn)對(duì)于圖像鄰域信息的使用過(guò)度，從而影響變化檢測(cè)精度。

對(duì)于差值圖像首先使用FCM算法，將圖像上的像素點(diǎn)初步分為兩個(gè)類別，這會(huì)得到兩個(gè)聚類中心。基于此結(jié)果，可以得到兩個(gè)閾值H1、H2，因此可以將差值圖像分為三部分。其中H1左側(cè)為未變化類別，H1與H2之間為不確定類別，H2右側(cè)為變化類別。H1與H2的計(jì)算方式如下：

(6)

式中：C1為未變化類別聚類中心；C2為變化類別聚類中心；Mmid為模糊C均值聚類后具有相同變化與未變化的隸屬度的像素點(diǎn)所具有的灰度值，即該像素點(diǎn)屬于變化類別與未變化類別的可能性相同，兩種隸屬度值均為0.5；α1和α2是兩個(gè)常數(shù)，可以用來(lái)調(diào)節(jié)三個(gè)部分的面積。該劃分如圖1所示。

圖1 變化、不確定、未變化類別劃分示意圖

在此基礎(chǔ)上，根據(jù)差值圖像中像素點(diǎn)的灰度值所在的區(qū)間，采用三個(gè)不同的表達(dá)式來(lái)計(jì)算所在區(qū)間的懲罰系數(shù)。于是修正式(4)中的固定懲罰系數(shù)β。其計(jì)算公式為：

(7)

式中：βm(X(i,j))是像素點(diǎn)位置為(i,j)的懲罰系數(shù)；Xmin和Xmax分別是差值圖像中灰度最小和最大的值。

當(dāng)圖像中灰度值小于H1時(shí)，隨著灰度值從H1減少到Xmin時(shí)，懲罰系數(shù)線性減小，這減小了具有相對(duì)較小的灰度值的像素對(duì)于空間鄰域信息的影響，這是因?yàn)檫@些像素點(diǎn)有極大的概率保持不變。同理，當(dāng)灰度值大于H2時(shí)，隨著灰度值由H2增加到Xmax，懲罰系數(shù)線性減小，這減小了具有相對(duì)較大的灰度值的像素對(duì)于空間鄰域信息的影響，這是因?yàn)檫@些像素點(diǎn)有極大的概率保持變化。對(duì)于灰度值處于H1與H2中間時(shí)，其屬于不確定的類別，對(duì)于這類像素點(diǎn)采用一個(gè)固定的懲罰系數(shù)。以此重新定義了使用的空間信息的權(quán)重。

2.2 變權(quán)指示函數(shù)

在傳統(tǒng)的基于MRF模型的變化檢測(cè)算法中假設(shè)某個(gè)像素點(diǎn)常常與其所在鄰域中的像素點(diǎn)具有相同的變化類別標(biāo)記。然而由于外部條件所限，獲得的遙感數(shù)據(jù)中包含大量混合像素。此外由于圖像數(shù)據(jù)和算法本身的固有不確定性，差異圖像的初始標(biāo)記類別并不是確定的，這尤其體現(xiàn)在圖像細(xì)節(jié)豐富的區(qū)域中。所以純粹地將指示函數(shù)的值取為-1和0是不合適的，該操作會(huì)顯得太絕對(duì)和不精確。本文對(duì)Potts模型進(jìn)行修正，來(lái)重新定義鄰域像素間的空間關(guān)系。

Foody[17]提出了一種基于熵的方法來(lái)度量鄰域空間的像素標(biāo)記的不確定性。采用該方法在模糊C均值聚類后，對(duì)于任意像素點(diǎn)可以得到其屬于變化與未變化兩個(gè)類別的隸屬度。故可得像素點(diǎn)xj的熵，定義如下：

Ej=-uj,ch×log2uj,ch-uj,un×log2uj,un

(8)

式中：uj,ch和uj,un分別是該像素點(diǎn)屬于變化與未變化類別的隸屬度值。于是可以將式(5)修正為：

(9)

(10)

這樣在計(jì)算空間能量函數(shù)時(shí)，可以充分考慮到鄰域像素的類別標(biāo)記的不確定性，使變化檢測(cè)結(jié)果的過(guò)度平滑問(wèn)題得到了很好的控制。故式(4)可以重新改寫為：

(11)

2.3 算法描述

該變化檢測(cè)算法的主要步驟概括如下：

1) 對(duì)配準(zhǔn)、校正后的兩時(shí)圖像做差獲得差值圖像。

3) 使用變權(quán)MRF模型來(lái)修正空間能量函數(shù)Uspatial(xi)的計(jì)算，使用改進(jìn)的變權(quán)懲罰系數(shù)和變權(quán)指示函數(shù)替換原始計(jì)算方式，重新定義空間上下文信息。

4) 使用ICM算法求解能量函數(shù)最小值，得到最終的變化檢測(cè)結(jié)果。

該算法的示意圖如圖2所示。

圖2 基于變權(quán)MRF變化檢測(cè)算法示意圖

3 實(shí) 驗(yàn)

3.1 實(shí)驗(yàn)數(shù)據(jù)及環(huán)境

實(shí)驗(yàn)使用的數(shù)據(jù)集為AICDDataset數(shù)據(jù)集和Szada數(shù)據(jù)集。AICDDataset為光學(xué)遙感模擬數(shù)據(jù)集，該數(shù)據(jù)集包含1 000組大小為800×600的兩時(shí)圖像及對(duì)應(yīng)的變化檢測(cè)參考圖像，空間分辨率為0.5 m×0.5 m。Szada為真實(shí)遙感影像數(shù)據(jù)集，包含7組大小為952×640的兩時(shí)圖像以及對(duì)應(yīng)的變化檢測(cè)參考圖像，空間分辨率為1.5 m×1.5 m。

分別進(jìn)行如下實(shí)驗(yàn)：(1) 從2個(gè)數(shù)據(jù)集中隨機(jī)選取2組數(shù)據(jù)，分別進(jìn)行2組模擬實(shí)驗(yàn)、2組真實(shí)實(shí)驗(yàn)進(jìn)行算法效果展示，所用的數(shù)據(jù)如圖3所示。(2) 從AICDDataset數(shù)據(jù)集中隨機(jī)抽取100組數(shù)據(jù)，比較本文算法與對(duì)比實(shí)驗(yàn)算法的效果。變化參考圖中白色區(qū)域代表兩時(shí)圖像發(fā)生變化的位置，黑色區(qū)域代表未發(fā)生變化的位置。

圖3 部分實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)所使用的數(shù)據(jù)均為經(jīng)過(guò)精確配準(zhǔn)、幾何校正、輻射校正的數(shù)據(jù)，消除了由于外部不良因素導(dǎo)致的偽變化，通過(guò)MATLAB R2015b編制程序求解。

3.2 實(shí)驗(yàn)結(jié)果及分析

將本文算法與文獻(xiàn)[10]和文獻(xiàn)[15]的變化檢測(cè)算法進(jìn)行對(duì)比，以上4組實(shí)驗(yàn)結(jié)果分別如圖4-圖7所示。本文采用Kappa系數(shù)來(lái)評(píng)價(jià)變化檢測(cè)的結(jié)果，其取值范圍在[-1,1]之間，Kappa系數(shù)的大小反映了變化檢測(cè)結(jié)果圖與變化參考圖之間一致性的程度，其值越接近于1，表示變化檢測(cè)結(jié)果圖與變化參考圖的匹配程度越高。

圖4 實(shí)驗(yàn)(1)結(jié)果圖

圖5 實(shí)驗(yàn)(2)結(jié)果圖

圖6 實(shí)驗(yàn)(3)結(jié)果圖

圖7 實(shí)驗(yàn)(4)結(jié)果圖

上述實(shí)驗(yàn)的Kappa系數(shù)統(tǒng)計(jì)結(jié)果如表1所示。

表1 Kappa系數(shù)統(tǒng)計(jì)表

從AICDDataset數(shù)據(jù)集中隨機(jī)抽取100組數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，分別統(tǒng)計(jì)本文算法與對(duì)比實(shí)驗(yàn)算法的Kappa系數(shù)，結(jié)果如圖8所示。

圖8 Kappa系數(shù)對(duì)比圖

上述兩次實(shí)驗(yàn)結(jié)果表明，在絕大部分情況下，本文提出的變化檢測(cè)算法與變化參考圖具有最高的相似度，取得了比對(duì)比算法更好的變化檢測(cè)結(jié)果。本文算法實(shí)現(xiàn)了對(duì)變化信息的有效提取，其含有較少的漏檢和誤檢像素，對(duì)噪聲等偽變化信息實(shí)現(xiàn)了良好的甄別，噪聲點(diǎn)更少，實(shí)驗(yàn)結(jié)果也更平滑。

4 結(jié) 語(yǔ)

針對(duì)傳統(tǒng)基于MRF模型變化檢測(cè)算法在能量函數(shù)的計(jì)算中使用固定的懲罰系數(shù)和固定指示函數(shù)的問(wèn)題，提出了一種基于變權(quán)MRF模型的無(wú)監(jiān)督的遙感變化檢測(cè)算法。該算法通過(guò)使用變權(quán)懲罰系數(shù)和基于熵的局部不確定鄰域標(biāo)記的變權(quán)指示函數(shù)，重新定義了空間信息權(quán)重和MRF模型中的鄰域像素的空間關(guān)系，進(jìn)而實(shí)現(xiàn)對(duì)空間能量函數(shù)的改造，優(yōu)化了傳統(tǒng)MRF模型。實(shí)驗(yàn)結(jié)果表明本文算法提供了比傳統(tǒng)MRF模型更精確的變化檢測(cè)，避免了變化檢測(cè)的過(guò)度平滑問(wèn)題。因此本文算法是一種可行、有效的變化檢測(cè)算法。