基于MDL準(zhǔn)則和EM算法的多光譜遙感影像變化檢測(cè)

石愛業(yè)，吳國(guó)寶，譚德寶，申邵洪，黃鳳辰，馬貞立，王 超

(1.河海大學(xué)計(jì)算機(jī)與信息學(xué)院，江蘇 南京 210098;2.長(zhǎng)江科學(xué)院空間信息技術(shù)應(yīng)用研究所，湖北武漢 430010)

基于遙感影像的變化檢測(cè)是根據(jù)同一區(qū)域不同時(shí)相的遙感影像辨識(shí)、觀測(cè)對(duì)象或者現(xiàn)象的狀態(tài)變化過程［1-2］。它已廣泛應(yīng)用于資源管理與規(guī)劃、環(huán)境保護(hù)等諸多領(lǐng)域［2-4］，并為相關(guān)部門提供科學(xué)決策的依據(jù)。當(dāng)前，遙感影像的變化檢測(cè)已成為一個(gè)非常活躍的研究方向［2-3，5-8］。遙感影像變化檢測(cè)中的代數(shù)法是當(dāng)前使用最廣泛的方法之一，但其存在一些不足:構(gòu)造差異影像的方法單一，不能很好地反映復(fù)雜地物的真實(shí)變化情況［9-10］;對(duì)于傳統(tǒng)的EM(expectation-maximization)算法［11］，當(dāng)差異影像含有多種地物變化類型、光譜信息較復(fù)雜時(shí)，其變換檢測(cè)精度較低［12］。

針對(duì)上述問題，筆者首先采用主成分(principle component analysis，PCA)變換結(jié)合相關(guān)系數(shù)融合法［9］來構(gòu)造差異影像;其次，利用分支數(shù)為k的高斯分布混合模型(gaussian mixture model，GMM)對(duì)差異影像進(jìn)行建模，并引入MDL(minmum description length)-EM算法來估計(jì)GMM的各參數(shù);最后，利用基于統(tǒng)計(jì)最小錯(cuò)誤率的Bayes判別準(zhǔn)則來確定變化檢測(cè)的閾值，得到變化檢測(cè)結(jié)果。通過實(shí)驗(yàn)數(shù)據(jù)比較分析，驗(yàn)證所提算法的可行性和有效性。

1 原理與步驟

1.1 基于PCA變換結(jié)合相關(guān)系數(shù)融合法的差異影像

多光譜遙感影像各波段上包含的信息存在一定程度的相關(guān)和冗余，因?yàn)镻CA變換能很好地把原影像中多個(gè)波段的有用信息壓縮到數(shù)目盡可能少的新成分影像中，并能使這些成分之間互不相關(guān)［12］，因此筆者首先對(duì)遙感影像數(shù)據(jù)進(jìn)行PCA變換提取第一主成分，然后采用文獻(xiàn)［9］所提的相關(guān)系數(shù)融合法來構(gòu)造不同時(shí)相間的差異影像，以實(shí)現(xiàn)最大限度地消除各光譜波段間的相關(guān)性，從而有效地集中和突出不同時(shí)相間的變化信息，并最大限度地保留變化的細(xì)節(jié)信息，減少偽變化信息及背景噪聲的干擾，為下一步的變化檢測(cè)創(chuàng)造有利條件［10］。

1.2 結(jié)合MDL準(zhǔn)則的EM算法

直方圖擬合是直方圖分析法的一種［13］，利用模型對(duì)目標(biāo)影像直方圖進(jìn)行擬合，這與進(jìn)行遙感影像變化檢測(cè)時(shí)將差異影像多類地物的光譜信息視作k類高斯分布的疊加相吻合。基于此，提出采用GMM來擬合差異影像直方圖。設(shè)L為差異影像灰度值集合(一般為0，1，…，255)，h(fi)為直方圖，其中fi為其第i個(gè)像元灰度值，利用EM算法估計(jì) GMM各參數(shù)［11］的迭代公式為

式中:p(j)，μj，——GMM中分支j的權(quán)重、均值和方差;t——迭代次數(shù)。GMM最優(yōu)模型的選擇由作為適應(yīng)度函數(shù)的MDL準(zhǔn)則來確定，其表達(dá)式為

其中 Θ =(p(1)，p(2)，…，p(k);μ1，μ2，…，μk;，，…)式中:k——GMM的分支數(shù);N——GMM的像元總數(shù);M——GMM中像元的維數(shù);ln P完全數(shù)據(jù)集的對(duì)數(shù)似然函數(shù);——引入的針對(duì)GMM分支數(shù)k的懲罰函數(shù)。

在選定初始最大分支數(shù)kmax的情況下，通過引入合并最小距離分支［14］的辦法自適應(yīng)地依次減少分支數(shù)至最小分支數(shù)kmin，并依據(jù)式(5)確定最優(yōu)的分支數(shù)k*和參數(shù)Θ*。

利用上述MDL-EM算法能自適應(yīng)地得到GMM的k*和Θ*，進(jìn)而得到對(duì)差異影像直方圖的最佳擬合，再利用基于最小錯(cuò)誤率的Bayes判別準(zhǔn)則確定變化檢測(cè)的閾值［11］。

1.3 新算法實(shí)現(xiàn)步驟

步驟1 對(duì)選取的2幅遙感影像TM1，TM2進(jìn)行幾何校正及輻射歸一化等預(yù)處理。

步驟2 分別對(duì)經(jīng)過預(yù)處理的TM1，TM2進(jìn)行PCA變換，并提取第一主成分X1pc1，X2pc1。

步驟3 求X1pc1，X2pc1的差值、比值影像，并通過相關(guān)系數(shù)融合法構(gòu)造差異影像F。

步驟4 設(shè)定初始參數(shù)kmax，kmin及Θ0，運(yùn)用MDL-EM算法估計(jì)k*及Θ*。

步驟5 使用GMM對(duì)差異影像F的直方圖進(jìn)行擬合。

步驟6 采用基于最小錯(cuò)誤率的Bayes判別準(zhǔn)則確定變化檢測(cè)的閾值，得到變化檢測(cè)結(jié)果。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)1

采用巴西朗多尼亞州亞馬遜流域的2幅Landsat TM(Thematic Mapper)遙感影像為實(shí)驗(yàn)對(duì)象，對(duì)其森林砍伐面積進(jìn)行變化檢測(cè)。2幅遙感影像的獲取時(shí)間分別為2000年7月和2006年7月，實(shí)驗(yàn)區(qū)大小為512像元 ×512像元，如圖1所示。該區(qū)域主要由河流、森林及砍伐后的裸地3部分組成。

在幾何校正及輻射歸一化的基礎(chǔ)上，采用基于PCA變換結(jié)合相關(guān)系數(shù)融合法對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行構(gòu)造，得到的差異影像如圖2所示，其像元可視為由分支數(shù)為k的高斯混合分布組成，對(duì)應(yīng)的直方圖如圖3所示。

圖1 Landsat TM數(shù)據(jù)Fig.1 Landsat TM data(true color)

圖2 差異影像Fig.2 Difference image

圖3 差異影像直方圖Fig.3 Histogram of difference image

在MDL-EM算法中，設(shè)置GMM的初始值kmax=30，kmin=2。經(jīng)過計(jì)算，各分支數(shù)k對(duì)應(yīng)的R如圖4所示。由圖4可以看到k=9時(shí)對(duì)應(yīng)的R最小，因此由MDL準(zhǔn)則可知，當(dāng)k=9時(shí)的模型是最佳的擬合模型，此時(shí)對(duì)差異影像直方圖擬合的效果最好，對(duì)應(yīng)的GMM各分支參數(shù)如表1所示。由于本實(shí)驗(yàn)的主要目的是檢測(cè)森林砍伐區(qū)變化(只檢測(cè)森林是否被砍伐)，暫不考慮其他的變化類型，所以將砍伐區(qū)統(tǒng)一歸為“變化類”，未砍伐區(qū)統(tǒng)一歸為“非變化類”。

根據(jù)表1，利用GMM對(duì)圖3所示的差異影像直方圖進(jìn)行擬合。由圖3可知，當(dāng)灰度值大于100時(shí)，其像素?cái)?shù)很少，直方圖與橫坐標(biāo)幾乎重合，此時(shí)進(jìn)行擬合結(jié)果比較意義不大。為了便于對(duì)擬合結(jié)果進(jìn)行較清晰的比較，僅將像元灰度值在0～100之間的擬合部分進(jìn)行放大顯示，k=9時(shí)的直方圖擬合結(jié)果如圖5(a)所示。為了比較MDL-EM算法求得的不同分支數(shù)對(duì)差異影像直方圖的擬合效果，選取具有代表性的k=2時(shí)的直方圖擬合結(jié)果作為比較，如圖5(b)所示。由圖5可以看出，在整體的擬合效果上，k=9的擬合效果要好于k=2的擬合效果。依據(jù)MDL-EM算法參數(shù)估計(jì)的結(jié)果和最小錯(cuò)誤率的Bayes判別準(zhǔn)則，確定變化檢測(cè)的閾值，進(jìn)而對(duì)圖2進(jìn)行二值化處理，得到最終變化檢測(cè)結(jié)果如圖6(a)所示。

圖4 不同k對(duì)應(yīng)的RFig.4 R with different k values

表1 k*=9時(shí)GMM各分支參數(shù)Table 1 GMM branch parameters by k*=9

圖5 差異影像直方圖擬合結(jié)果Fig.5 Fitted histogram of difference image

為了比較新算法的優(yōu)劣，選取兩類變化檢測(cè)算法做進(jìn)一步的對(duì)比分析:(a)Bruzzone所提的CVA-EM算法［11］;(b)PCA融合－EM算法。所有對(duì)比實(shí)驗(yàn)的EM算法均設(shè)定 GMM的分支數(shù)k=2，EM算法初始值設(shè)置方法和文獻(xiàn)［11］相同。2種算法的檢測(cè)結(jié)果分別如圖6(b)(c)所示;目視解譯得到的理想變化檢測(cè)結(jié)果如圖6(d)所示。

圖6 實(shí)驗(yàn)數(shù)據(jù)1的3種算法及目視解譯變化檢測(cè)結(jié)果Fig.6 Change detection map by three algorithms and visual interpretation in experimental group 1

采用誤檢率、漏檢率、總體精度、Kappa系數(shù)和計(jì)算時(shí)間5個(gè)指標(biāo)來評(píng)價(jià)新算法及其他2種變化檢測(cè)算法的性能，結(jié)果如表2所示。由表2可知，較于其他算法，新算法除計(jì)算代價(jià)稍大外，在總體檢測(cè)性能上是最好的，因此新算法提高了多光譜遙感影像的變化檢測(cè)精度。

表2 實(shí)驗(yàn)數(shù)據(jù)1的各變化檢測(cè)算法性能評(píng)價(jià)Table 2 Performance evaluation of three algorithms in experimental group 1

2.2 實(shí)驗(yàn)數(shù)據(jù)2

采用美國(guó)加利福尼亞州索爾頓湖卡利帕特里亞地區(qū)的2幅遙感影像，獲取時(shí)間分別為2008年8月和2009年8月，實(shí)驗(yàn)區(qū)大小為400像元×400像元，如圖7所示。該區(qū)域地物類型較為豐富，包含有水域、山地、城市、農(nóng)用耕地、裸地等，同時(shí)通過目視解譯發(fā)現(xiàn)該區(qū)域地物變化信息也較為復(fù)雜，如水域變成裸地、山地變成農(nóng)用耕地、農(nóng)用耕地變成裸地等。

利用MDL-EM算法對(duì)GMM的各參數(shù)進(jìn)行估計(jì)，得到k*=10，最終的變化檢測(cè)結(jié)果如圖8(a)所示。采用與第1個(gè)實(shí)驗(yàn)相同的比較算法，且各算法的初始參數(shù)設(shè)置亦相同，得到的變化檢測(cè)結(jié)果如圖8(b)(c)所示;目視解譯得到的理想變化檢測(cè)結(jié)果如圖8(d)所示。各算法的性能評(píng)價(jià)見表3。

通過比較圖8的變化檢測(cè)結(jié)果及由表3的各變化檢測(cè)算法性能評(píng)價(jià)可知，相比于其他的算法，本文所提新算法在總體檢測(cè)性能上是最好的，這也表明所提新算法用于多光譜遙影像的變化檢測(cè)是可行的。

圖7 Landsat TM數(shù)據(jù)Fig.7 Landsat TM data(true color)

圖8 實(shí)驗(yàn)數(shù)據(jù)2的3種算法及目視解譯的變化檢測(cè)結(jié)果Fig.8 Change detection map by three algorithms and visual interpretation in experimental group 2

表3 實(shí)驗(yàn)數(shù)據(jù)2的各變化檢測(cè)算法性能評(píng)價(jià)Table 3 Performance evaluation of three algorithms in experimental group 2

3 結(jié) 語

通過引入PCA變換并結(jié)合相關(guān)系數(shù)融合法來構(gòu)造差異影像，將MDL準(zhǔn)則引入EM算法中估計(jì)高斯混合模型的分支數(shù)及各分支參數(shù)。實(shí)驗(yàn)分析表明，所提新的變化檢測(cè)算法是有效的。新算法中，GMM的初始分支數(shù)通過試錯(cuò)法來選擇。另外，新算法是在假定差異影像像元間相互獨(dú)立的基礎(chǔ)上進(jìn)行的，并未考慮像元的空間關(guān)系，如何將這種空間關(guān)系嵌入到變化檢測(cè)算法中以提高檢測(cè)性能也是今后需要進(jìn)一步研究的問題。

［1］SING A.Change detection techniques using remotely sensed data［J］.International Journal of Remote Sensing，1989，10(6):989-1003.

［2］LU D，MAUSEL P，MORAN E.Change detection techniques［J］.International Journal of Remote Sensing，2004，25(12):2365-2407.

［3］BOVOLO F，BRUZZONE L.A split-based approach to unsupervised change detection in large-size multi-temporal images:application to tsunami-damage assessment［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45(6):1658-1670.

［4］CHINI M，PACIFICI F，EMERY W.Comparing statistical and neural network methods applied to very high resolution satellite images showing changes in man-made structures at rocky flats［J］.IEEE Transactions on Geoscience and Remote Sensing，2008，46(6):1812-1821.

［5］GABRIELE M，ELENA A，SEBASTIANO B S.Multiscale unsupervised change detection on optical images by Markov random fields and wavelets［J］.IEEE Geoscience and Remote Sensing Letters，2011，8(4):725-729.

［6］MERCIER G，SERPICO S.Conditional copulas for change detection in heterogeneous remote sensing images［J］.IEEETransactions on Geoscience and Remote Sensing，2008，46(5):1428-1441.

［7］LUO W，LI H L.PCA based unsupervised change detection for color satellite images under the quaternion model［C］//International Symposium on Intelligent Signal Processing and Communication Systems.Chengdu:Journal of University of Electronic Science and Technology of China，2010:1-4.

［8］GABRIELE M，ELENA A，SEBASTIANO B S.A contextual multiscale unsupervised method for change detection with multitemporal remote-sensing images［C］//2009 Ninth International Conference on Intelligent Systems Design and Application.Pisa:Journal of University of Pisa，2009:572-577.

［9］魏立飛，鐘燕飛，張良培，等.遙感影像融合的自適應(yīng)變化檢測(cè)［J］.遙感學(xué)報(bào)，2010，14(6):1204-1211.(WEI Lifei，ZHONG Yanfei，ZHANG Liangpei，et al.Adaptive change method of remote sensing image fusion［J］.Journal of Remote Sensing，2010，14(6):1204-1211.(in Chinese))

［10］柳思聰，杜培軍，陳紹杰.決策級(jí)融合的多分辨率遙感影像變化檢測(cè)［J］.遙感學(xué)報(bào)，2011，15(4):846-862.(LIU Sicong，DU Peijun，CHEN Shaojie.A novel change detection method of multi-resolution remotely sensed images based on the decision level fusion［J］.Journal of Remote Sensing，2011，15(4):846-862.(in Chinese))

［11］BRUZZONE L，PRIETOD F.Automatic analysis of difference image for unsupervised change detection［J］.IEEE Transactions on Geoscience and Remote Sensing，2000，38(3):1171-1182.

［12］吳柯，牛瑞卿，王毅，等.基于PCA與EM算法的多光譜遙感影像變化檢測(cè)研究［J］.計(jì)算機(jī)科學(xué)，2010，37(3):282-284.(WU Ke，NIU Ruiqing，WANG Yi，et al.Change detection of multispectral remote sensed images based on PCA and EM algorithm［J］.Computer Science，2010，37(3):282-284.(in Chinese))

［13］鄒丹平，馮濤，李咸偉，等.基于EM的直方圖逼近及其應(yīng)用［J］.中國(guó)圖像圖形學(xué)報(bào)，2005，10(11):1458-1462.(ZOU Danping，F(xiàn)ENG Tao，LI Xianwei，et al.Histogram approximation based on expectation maximization algorithm and its application［J］.Journal of Image and Graphics，2005，10(11):1458-1462.(in Chinese))

［14］XIE W B，WANG X H.EM algorithm based MDL application to estate the mixture model clustering parameters［C］//The 9th International Conference on Signal Processing.Beijing:International Journal of Scientific＆ Technology Research，2008:1637-1640.