王 雷，劉小芳，趙良軍，石小仕

(1.四川輕化工大學(xué)自動化與信息工程學(xué)院，四川 宜賓 644000；2.四川輕化工大學(xué)計算機(jī)學(xué)院，四川 宜賓 644000)

1 引言

影像的多時相變化檢測是指對同一地點、同一角度下不同時間點影像進(jìn)行對比分析，從而得到地表變化信息的過程[1]。近年來，隨著遙感衛(wèi)星的商用程度不斷提高，可獲取的遙感影像數(shù)據(jù)越來越多，因此基于遙感影像的變化檢測方法在土地動態(tài)監(jiān)測、城市擴(kuò)展、災(zāi)害監(jiān)測、森林面積變化等領(lǐng)域發(fā)揮著重要作用，是遙感影像應(yīng)用的一個重要方向。但是，到目前為止還沒有出現(xiàn)一種穩(wěn)健的方法可以通過自動化的手段從多時相遙感影像中高準(zhǔn)確率地解析出變化信息，因此，發(fā)展可靠、高準(zhǔn)確率的算法仍然是研究的熱點[2]。目前遙感影像的變化檢測按照研究目標(biāo)的單位可以劃分為像素級[3-6]、特征級[7，8]、對象級[9-11]?；谙袼丶壍淖兓瘷z測因其具有直觀、可視化程度高的特點得到的廣泛的應(yīng)用，絕大多數(shù)變化檢測是在差異影像的基礎(chǔ)上來完成的，常用代數(shù)法，包括影像差分法、影像回歸法、影像比值法、相關(guān)系數(shù)法、變化向量分析法等；變化法，包括主成分分析法、K-T變換法、G-S法等。由于拍攝季節(jié)、拍攝角度姿態(tài)和陽光照射等客觀原因，預(yù)處理之后的影像往往還存在一定的輻射誤差和配準(zhǔn)誤差。因此傳統(tǒng)的基于差異影像的分割方法往往存在大量的椒鹽噪聲和像元斑點，導(dǎo)致變化檢測精度低。Celik[12]利用主成分分析(Principal component analysis，PCA)將局部信息投影到特征空間形成特征向量，然后利用K均值聚類算法(K-means clustering，Kmeans)進(jìn)行變化檢測(PCA-Kmeans)，該方法不能很好地去除噪聲，并且細(xì)節(jié)信息丟失較多。模糊C均值聚類算法(Fuzzy c-means clustering，F(xiàn)CM)[13，14]作為基于模糊理論產(chǎn)生的算法，在變化檢測中有很廣泛的應(yīng)用，但是FCM算法在聚類過程中易受噪聲和初試聚類中心的影響，進(jìn)而產(chǎn)生大量的虛警[15]。毛天祺等提出利用二進(jìn)小波增強(qiáng)與邊緣局部信息FCM算法的變化檢測方法，但是該方法在某種情況下會出現(xiàn)變化檢測精度低的問題[16]。羅文村[17]、郭臻[18]等人提出，在變化強(qiáng)度圖像基礎(chǔ)上，對變化強(qiáng)度設(shè)置合適的閾值，將變化強(qiáng)度圖二值化，進(jìn)而區(qū)分變化和未變化類別。目前已有大量關(guān)于自動閾值獲取方法的研究，并成功運(yùn)用于圖像分割，其中以O(shè)tsu法最為典型[19]，但基于閾值分割的方法無法充分利用多時相影像的變化信息，易受椒鹽噪聲影響導(dǎo)致精度較低。張道強(qiáng)等[20]提出使用基于核改進(jìn)的FCM算法做變化檢測的圖像分割(K-FCM)，實驗表明扔不能達(dá)到較高精度。分類結(jié)果衡量指標(biāo)方面，David L.Davies等人[21]提出的一種評估聚類算法優(yōu)劣程度的指標(biāo)(DB指標(biāo))。

針對上述變化檢測算法對遙感多時相影像變化檢測精度較低的問題，提出一種基于粒子群算法和改進(jìn)FCM算法的方法，將對數(shù)、差分、比值對應(yīng)的變化強(qiáng)度圖像各自的優(yōu)勢整合起來，形成3通道的變化強(qiáng)度圖，然后應(yīng)用PCA主成分分析法提取變化強(qiáng)度圖的主要特征，并將變化強(qiáng)度圖轉(zhuǎn)換為向量集投影到主成分空間形成特征空間，最后利用改進(jìn)的FCM算法進(jìn)行二分類聚類，劃分為變化和非變化兩個類別，最終生成結(jié)果圖。

2 改進(jìn)的FCM算法

用比值法、自然對數(shù)法、差值法獲取三種變化強(qiáng)度圖像，并合成三通道的變化強(qiáng)度圖，然后用PCA做主成分提取，并將變化強(qiáng)度圖映射到主成分空間中，再用傳統(tǒng)FCM算法進(jìn)行二分類聚類獲得聚類中心用于初始化粒子群，粒子群以DB指標(biāo)為適應(yīng)度，在全局空間內(nèi)搜索最優(yōu)值，獲得最優(yōu)劃分聚類中心，最后輸出變化檢測結(jié)果圖。

2.1 獲取三通道變化強(qiáng)度圖像

假設(shè)時相1影像的矩陣表達(dá)為XT1，時相2影像表達(dá)為XT2。

2.1.1 對數(shù)比率法

對數(shù)比率圖像Id1通過對不同時相影像對應(yīng)像素值做比的結(jié)果取對數(shù)計算得到，使用式(1)進(jìn)行計算

(1)

其中e代表自然對數(shù)下標(biāo)，C代表小常數(shù)，對數(shù)比率法的特點是能夠加強(qiáng)比值圖像強(qiáng)度較低的部分。

2.1.2 絕對差值法

絕對差值圖像Id2是由兩時相影像做差取絕對值獲得的，在兩時相影像中變化的區(qū)域在差值圖像中擁有較大的絕對值，未變化的區(qū)域則擁有更小的絕對值。絕對差值變化圖像使用式(2)進(jìn)行計算

Id2=|XT1-XT2|

(2)

2.1.3 簡單比值法

在簡單比值法圖像Id3中，將對應(yīng)多時相影像對應(yīng)位置上的像素值做除，未變化的區(qū)域比值接近于1，而變化區(qū)域的比值與1的絕對差值較大。使用式(3)進(jìn)行計算

(3)

簡單比值法生成的圖像能夠有效減少太陽仰角、陰影和地形的影響。

2.1.4 3通道變化強(qiáng)度圖像

將由式(1)-(3)得到的三種不同表達(dá)的變化強(qiáng)度圖像表示成一個三通道的變化強(qiáng)度圖像，其中一個像素用Id(x，y)來表達(dá)，其中x代表該像素點在圖像中的行坐標(biāo)，y代表在圖像中的列坐標(biāo)，因此對于每一個像素點可以用式(4)描述為

(4)

假設(shè)最終的變化圖像的寬度為W，高度為H，那么Id(x，y)這樣的向量就有H×W個。為了方便后續(xù)表達(dá)，將Id(x，y)表示為Id，k，其中k表示的是圖像中的第k個像素，1≤k≤N，N=W×H。計算三通道像素的平均值，得到平均值向量χ，χ根據(jù)式(5)進(jìn)行計算

(5)

對變化強(qiáng)度圖像做均值化處理得到R，R由式(6)計算

Rdi=Idi(x，y)-χdi

(6)

2.2 用PCA提取變化強(qiáng)度圖像主要特征

PCA主成分提取技術(shù)已經(jīng)廣泛應(yīng)用于特征的提取和降維計算[22]，于是根據(jù)式(6)計算的到的變化強(qiáng)度圖像各位置上的Rdi，其協(xié)方差矩陣CI可以由式(7)表示出來

(7)

由式(8)計算協(xié)方差矩陣CI的特征值

CIvi=λivi

(8)

其中，λi表示協(xié)方差矩陣的第i個特征值，與其對應(yīng)的特征向量是vi。對特征向量進(jìn)行降序排列，為

D=desc(λ1，λ2，…，λn)

(9)

其中，λ1，λ2，…，λn是協(xié)方差矩陣的特征值的降序排列，對應(yīng)的特征向量a1，a2…，an所構(gòu)成的正交向量為

A=(a1，a2…，an)

(10)

協(xié)方差的變換矩陣為

(11)

其中，(i1，i2，…in)T為原始圖像的n維向量表達(dá)，a11，a12，…ann為特征值對應(yīng)特征向量鎖構(gòu)成的協(xié)方差矩陣對應(yīng)位置上的值，(y1，y2，…yn)T為原始圖像的n維向量表達(dá)在PCA主成分空間上的投影，矩陣Y為變換后的主成分，前幾個行向量為所求原始影像的主要特征信息，每個特征的權(quán)重由其特征值在所有特征值中所占的比重來決定，可以由式(12)計算得到，本實驗選取特征權(quán)重是比例是0.9，即取累計權(quán)重超過0.9的前幾維特征。

(12)

2.3 基于粒子群和DB指標(biāo)改進(jìn)的FCM算法

2.3.1 標(biāo)準(zhǔn)FCM算法

FCM算法是一種對數(shù)據(jù)樣本點做無監(jiān)督自動分類的方法，在數(shù)據(jù)挖掘領(lǐng)域得到了廣泛的運(yùn)用，該算法通過計算數(shù)據(jù)點屬于各個樣本中心的隸屬度，對各個樣本點進(jìn)行分類。目標(biāo)函數(shù)Jm(U，V)如式(13)所示，聚類的最終目的是尋找一組聚類中心點，使得目標(biāo)函數(shù)取得極小值，具體實現(xiàn)步驟如下

(13)

step1：其中c(2≤c≤n)是聚類的種數(shù)，n是樣本的數(shù)量，dij表示第j個樣本點到第i個中心點的歐氏距離，uij代表第j個樣本點屬于第i個聚類中心的隸屬度，設(shè)定迭代的初始聚類中心V0，和初始隸屬度矩陣U0，設(shè)定模糊權(quán)重指數(shù)m(m>1)，最大迭代次數(shù)Cmax和兩次聚類中心變化差值最小的停止閾值ε。

(14)

step3：根據(jù)新一輪的隸屬度矩陣Ur+1通過式(15)計算新的聚類中心Vr+1，其中xk表示第k個數(shù)據(jù)樣本點的向量表達(dá)。

(15)

2.3.2 DB指標(biāo)

用于衡量分類結(jié)果的好壞，該指標(biāo)不僅考慮類間方差對聚類結(jié)果的影響，還將聚類中心之間的相似程度作為考慮因素，為聚類操作提供指導(dǎo)，這個指標(biāo)越小代表聚類結(jié)果越好，其具體的計算步驟如下

step1：根據(jù)式(16)計算Si，q

(16)

其中xj代表第j個數(shù)據(jù)，zi指第i個聚類中心，Ni是類別劃分為聚類中心zi的數(shù)據(jù)個數(shù)，q是常數(shù)，當(dāng)q=2時，Si，q等價于所有樣本點到對應(yīng)聚類中心距離的標(biāo)準(zhǔn)差。

step2：根據(jù)式(17)計算Mij

(17)

ak，i代表中心i的第k個值，Mij定義為聚類中心i和聚類中心j的質(zhì)心距離。

step3：根據(jù)式(18)、(19)計算Ri，j和Ri，qt

(18)

(19)

step4：根據(jù)式(20)計算DB指標(biāo)

(20)

式中K為聚類中心的數(shù)量。

2.3.3 粒子群算法PSO

PSO粒子群算法由Kennedy等人于1995年提出[23]，它同時具備進(jìn)化和群智能計算的優(yōu)點，是一種啟發(fā)式的全局優(yōu)化算法。與其它群算法的思路類似，PSO算法是通過個體之間的協(xié)作與競爭實現(xiàn)復(fù)雜問題最優(yōu)解的全局搜索。

PSO算法可以描述為，設(shè)粒子在D維的解空間進(jìn)行最優(yōu)解搜索，粒子群包含的粒子數(shù)量為N，第k個粒子在解空間中的位置描述為Xk=(xk1，xk2，…，xkD)，粒子在解空間中的飛行速度可以描述為Vk=(vk1，vk2，…，vkD)，每一個粒子代表空間中的一個解，粒子通過不斷在解空間中搜索以尋找最優(yōu)解。第k個粒子到目前為止找到的最優(yōu)位置記錄為Pk=(pk1，pk2，…，pkD)，整個粒子群所搜索到的歷史最優(yōu)位置記錄為，每個粒子在一次迭代中的位置和速度按照式(21)和(22)進(jìn)行變化

vkd(t+1)=wvkd+c1r1(pkd(t)-xkd(t))+

c2rr(pgd(t)-xkd(t))

(21)

xkd(t+1)=xkd(t)+vkd(t+1)

(22)

式中，r1和r2分別為范圍在[0,1]之間的隨機(jī)數(shù)；c1和c2分別代表粒子的自我認(rèn)知系數(shù)和社會認(rèn)知系數(shù)；w是粒子搜索的慣性系數(shù)，可以控制粒子更新速度受上一次速度的影響，每個粒子第d維的速度和位置變化范圍為[-vd,max,vd,max]和[-xd,max,xd,max]，粒子速度太大可能會使粒子飛過最好解，速度太小可能使粒子群陷入局部最優(yōu)解。

2.3.4 基于粒子群和DB指標(biāo)的優(yōu)化FCM算法的實現(xiàn)步驟

基于上述原理，提出一種基于DB指標(biāo)改進(jìn)的FCM算法，并用PSO算法克服FCM算法易受噪聲影響和對初始聚類中心敏感的缺陷，用于遙感滑坡的變化檢測，具體計算步驟如下：

step1：根據(jù)多時相遙感影像，采用式(1)-(3)從不同角度獲取變化強(qiáng)度圖。

step2：利用式(4)-(6)獲取3通道變化圖像，并做均值化處理。

step3：根據(jù)式(7)-(11)對變化強(qiáng)度圖做主成分特征提取，然后將其對應(yīng)的協(xié)方差矩陣的特征值按大小降序排列，用式(12)計算各個特征值的權(quán)重，最后根據(jù)特征值權(quán)重累加閾值比例，提取前面數(shù)維特征為主要特征信息。

step4：根據(jù)式(13)-(15)描述的傳統(tǒng)FCM算法在step3中形成的向量空間中進(jìn)行聚類，得到一個聚類中心。

step5：初始化粒子群，將step4中獲得的聚類中心作為一個粒子初始化到粒子群中，起加速收斂的作用。

step6：根據(jù)式(16)-(20)，計算各個粒子的適應(yīng)度DB指標(biāo)。

step7：根據(jù)適應(yīng)度，更新粒子群歷史最優(yōu)位置Pg和每一個粒子對應(yīng)的歷史最優(yōu)位置Pk。

step8：根據(jù)式(21)-(22)，更新粒子的速度Vk(t+1)和位置Xk(t+1)。

step9：重復(fù)step 6)至step 8)，直到達(dá)到最大迭代次數(shù)或滿足收斂條件。

step10：輸出最優(yōu)聚類中心。

3 實驗與仿真

3.1 實驗數(shù)據(jù)

實驗采用截取部分汶川縣綿虒鎮(zhèn)2005年和2011年兩個時相的Google Earth的19級遙感影像作為數(shù)據(jù)，其空間分辨率大小為2m，圖像為RGB3通道彩色影像，通過對比2008年汶川地震前后影響差異提取變化信息。在進(jìn)行試驗之前先對兩時相影響進(jìn)行了配準(zhǔn)，配準(zhǔn)誤差小于0.5像素，兩時相影像如圖1(a)、(b)所示，圖片大小為2261×2215像素，從圖1中可以看出，經(jīng)過汶川地震之后，兩時相數(shù)據(jù)發(fā)生了明顯的變化，主要變化包含新生滑坡信息和耕地變化。

圖1 汶川綿虒鎮(zhèn)2005和2011兩時相影像

為了驗證本文算法的有效性，對上述數(shù)據(jù)進(jìn)行變化檢測實驗，實驗硬件如下。處理器為Intel Core i7-6700HQ，內(nèi)存為16GB。實驗軟件如下：操作系統(tǒng)Windows 10專業(yè)版，仿真平臺為Matlab 2016b。

3.2 實驗步驟

實驗可以分為三大步驟，首先是準(zhǔn)備數(shù)據(jù)，計算獲取變化強(qiáng)度圖，再利用主成分分析法提取變化強(qiáng)度圖像的主要特征，并將樣本集投影到特征空間得到樣本空間，最后利用改進(jìn)的FCM算法做聚類得到最優(yōu)中心，獲取變化結(jié)果圖。操作流程圖如圖2所示。

圖2 滑坡區(qū)域變化檢測操作流程圖

3.3 實驗結(jié)果分析

首先利用對數(shù)法、差值法、比值法獲取不同的變化強(qiáng)度圖像如圖3(a)、(b)、(c)所示，然后將其組成一個3通道圖像，將變化強(qiáng)度圖像按像素展開成1×3的向量，構(gòu)成向量空間，使用PCA主成分分析法提取第一主分量，提取到的特征如圖4所示。

圖3 變化強(qiáng)度圖

圖4 三通道變化強(qiáng)度影像主成分特征圖

為了達(dá)到在空間中盡可能分散數(shù)據(jù)點以提升聚類效果的目的，將變化圖像主要特征圖按照4×4的像素塊展開為以1×16向量構(gòu)成的向量空間，使用PCA主成分分析法提取相應(yīng)協(xié)方差矩陣的特征值和特征向量，然后以0.9為特征值累加閾值選取對應(yīng)特征向量前幾維分量作為新特征空間的正交基，最后將變化強(qiáng)度圖像像素對應(yīng)的向量空間投影到新特征空間中，最后使用本文提出的算法進(jìn)行聚類，輸出差異二值圖。

為了證明本文算法的有效性，將算法結(jié)果與四種現(xiàn)有的變化檢測方法進(jìn)行比較，這四種方法包括：Otsu自動確定閾值法、PCA-Kmeans法、PCA-FCM算法和核聚類方法(KFCM)。

本方法和各方法執(zhí)行結(jié)果以及目視解譯標(biāo)準(zhǔn)結(jié)果如圖5所示。

圖5 各方法結(jié)果對比

為了更清楚地觀察各個算法的變化檢測結(jié)果，在每個結(jié)果圖中選擇了兩個比較具有代表性的區(qū)域，用綠色實線矩形將展示區(qū)域從結(jié)果圖中框出，如上圖(a)-(f)所示，其中每幅圖右邊放大圖像即結(jié)果圖中對應(yīng)位置上的細(xì)節(jié)圖，細(xì)節(jié)圖和對應(yīng)位置由紅色箭頭連接。

對比圖5中的(a)、(b)、(d)三幅結(jié)果圖像，可以看出，與傳統(tǒng)的K-FCM算法相比，基于主成分分析的PCA-Kmeans和PCA-FCM算法具有較小的變化檢測區(qū)域，椒鹽噪聲明顯少于KFCM算法，但是相對于本文算法結(jié)果(e)，丟失細(xì)節(jié)信息較多；KFCM算法結(jié)果中變化區(qū)域較豐富，但虛檢現(xiàn)象較嚴(yán)重；Otsu算法結(jié)果圖(c)與本文算法較接近，對比本文結(jié)果(e)的細(xì)節(jié)圖可以明顯分辨出其受椒鹽噪聲影響更嚴(yán)重。

綜上所述，相較于對比的4個算法，本文算法可以在較大程度保留變化細(xì)節(jié)特征的同時，有效減少椒鹽噪聲，獲得更好的可視化結(jié)果。

為了定量評估本文算法和各種對比算法獲得的結(jié)果，引入四個量化計算指標(biāo)，1)漏檢率，即事實上變化了的像素未被檢測出的像素占總像素數(shù)量的百分比；2)虛檢率，即原本未發(fā)生變化的像素被檢測為變化像素數(shù)量所占總像素的百分比；3)總體正確率，即正確檢測的變化像素和非變化像素數(shù)量占總像素的百分比；4)Kappa系數(shù)，Kappa系數(shù)用來在變化檢測中檢測實驗結(jié)果與地面真實情況之間的相似程度，Kappa系數(shù)越大，說明檢測結(jié)果與真實變化情況越相似，當(dāng)Kappa系數(shù)為1的時候，說明兩幅圖完全一致。各個變化檢測算法定量分析結(jié)果對照如表1所示

表1 各個變化檢測算法定量分析結(jié)果

由表1數(shù)據(jù)可知，本文算法的漏檢率為2.99%稍高于K-FCM算法的1.23%，這是因為K-FCM算法的虛檢率高達(dá)14.42%，進(jìn)而導(dǎo)致其漏檢率偏低；本文算法虛檢率僅為1%，遠(yuǎn)低于其余四種算法；正確率94.03%相較于其它算法也有明顯優(yōu)勢；Kappa系數(shù)方面，本文算法結(jié)果的Kappa系數(shù)為0.86明顯高于其它幾種算法，說明本文的變化檢測二分類結(jié)果與地面實際變化情況相似度較高，因此，在遙感圖像變化檢測領(lǐng)域，本文算法具有較好的實驗結(jié)果。

4 結(jié)論

利用三種不同的多時相影像變化強(qiáng)度計算方式獲取三種不同的變化強(qiáng)度圖，組成一個三通道的變化強(qiáng)度圖像，解決了多時相影像變化信息利用率低的問題；利用主成分分析法獲取變化的第一主分量信息，再將主要變化信息投影到主成分空間；引入DB指標(biāo)聚類評價指標(biāo)優(yōu)化傳統(tǒng)FCM算法過程，使用PSO解決傳統(tǒng)FCM算法易受孤立點和初始聚類中心影響而陷入局部最優(yōu)的問題；使用改進(jìn)的FCM算法對投影到主成分空間中的主要變化信息進(jìn)行聚類。實驗結(jié)果表明，相對于其它4種常見變化檢測算法，本文算法能夠一定程度上減小椒鹽噪聲對變化檢測結(jié)果的影響，并且可以很好地保留圖像的特征信息，提升了變化檢測精度。但是，本文算法扔無法完全解決基于像素的變化檢測方法易受椒鹽噪聲影響的問題，因此本算法如何在保持細(xì)節(jié)特征不丟失的前提下，進(jìn)一步減少椒鹽噪聲的問題是進(jìn)一步的研究方向。