基于改進FCM算法的遙感影像滑坡變化檢測方法

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

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

1 引言

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

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

2 改進的FCM算法

用比值法、自然對數法、差值法獲取三種變化強度圖像，并合成三通道的變化強度圖，然后用PCA做主成分提取，并將變化強度圖映射到主成分空間中，再用傳統FCM算法進行二分類聚類獲得聚類中心用于初始化粒子群，粒子群以DB指標為適應度，在全局空間內搜索最優值，獲得最優劃分聚類中心，最后輸出變化檢測結果圖。

2.1 獲取三通道變化強度圖像

假設時相1影像的矩陣表達為XT1，時相2影像表達為XT2。

2.1.1 對數比率法

對數比率圖像Id1通過對不同時相影像對應像素值做比的結果取對數計算得到，使用式(1)進行計算

(1)

其中e代表自然對數下標，C代表小常數，對數比率法的特點是能夠加強比值圖像強度較低的部分。

2.1.2 絕對差值法

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

Id2=|XT1-XT2|

(2)

2.1.3 簡單比值法

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

(3)

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

2.1.4 3通道變化強度圖像

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

(4)

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

(5)

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

Rdi=Idi(x，y)-χdi

(6)

2.2 用PCA提取變化強度圖像主要特征

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

(7)

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

CIvi=λivi

(8)

其中，λi表示協方差矩陣的第i個特征值，與其對應的特征向量是vi。對特征向量進行降序排列，為

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

(9)

其中，λ1，λ2，…，λn是協方差矩陣的特征值的降序排列，對應的特征向量a1，a2…，an所構成的正交向量為

A=(a1，a2…，an)

(10)

協方差的變換矩陣為

(11)

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

(12)

2.3 基于粒子群和DB指標改進的FCM算法

2.3.1 標準FCM算法

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

(13)

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

(14)

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

(15)

2.3.2 DB指標

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

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

(16)

其中xj代表第j個數據，zi指第i個聚類中心，Ni是類別劃分為聚類中心zi的數據個數，q是常數，當q=2時，Si，q等價于所有樣本點到對應聚類中心距離的標準差。

step2：根據式(17)計算Mij

(17)

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

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

(18)

(19)

step4：根據式(20)計算DB指標

(20)

式中K為聚類中心的數量。

2.3.3 粒子群算法PSO

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

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

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]之間的隨機數；c1和c2分別代表粒子的自我認知系數和社會認知系數；w是粒子搜索的慣性系數，可以控制粒子更新速度受上一次速度的影響，每個粒子第d維的速度和位置變化范圍為[-vd,max,vd,max]和[-xd,max,xd,max]，粒子速度太大可能會使粒子飛過最好解，速度太小可能使粒子群陷入局部最優解。

2.3.4 基于粒子群和DB指標的優化FCM算法的實現步驟

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

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

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

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

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

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

step6：根據式(16)-(20)，計算各個粒子的適應度DB指標。

step7：根據適應度，更新粒子群歷史最優位置Pg和每一個粒子對應的歷史最優位置Pk。

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

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

step10：輸出最優聚類中心。

3 實驗與仿真

3.1 實驗數據

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

圖1 汶川綿虒鎮2005和2011兩時相影像

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

3.2 實驗步驟

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

圖2 滑坡區域變化檢測操作流程圖

3.3 實驗結果分析

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

圖3 變化強度圖

圖4 三通道變化強度影像主成分特征圖

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

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

本方法和各方法執行結果以及目視解譯標準結果如圖5所示。

圖5 各方法結果對比

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

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

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

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

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

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

4 結論

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