基于選權(quán)迭代估計(jì)與非監(jiān)督分類的多光譜圖像變化檢測

李莎， 倪維平， 嚴(yán)衛(wèi)東， 吳俊政， 張晗

(西北核技術(shù)研究所，西安 710024)

0 引言

多光譜圖像的變化檢測是對不同時(shí)間獲取的多光譜遙感圖像進(jìn)行定量分析，并確定圖像區(qū)域內(nèi)地表變化特征和過程的技術(shù)[1]。從技術(shù)處理流程來看，一般有圖像預(yù)處理、變化信息發(fā)現(xiàn)、變化區(qū)域提取與變化類型確定等過程[2]，其中最為關(guān)鍵的是變化信息發(fā)現(xiàn)，大多數(shù)研究都是圍繞這一問題開展的。近年來，國內(nèi)外相繼發(fā)展了許多適用于遙感圖像變化檢測的方法[3-6]，但沒有一種方法是最優(yōu)的，每種方法都只具有一定的適用范圍。為了解決傳統(tǒng)方法(圖像比值法、變化向量分析法、主成分分析法[7-8]等)應(yīng)用于多波段遙感圖像變化檢測中的不足(抑制噪聲、消除相關(guān)性等)，丹麥學(xué)者Nielsen等[9-11]率先提出了多元變化檢測(multivariate alteration detection，MAD)的概念和方法。廖明生等[12]和陳壘等[13]將這種方法應(yīng)用于多時(shí)相多波段遙感圖像變化檢測，實(shí)驗(yàn)結(jié)果也顯示出了該方法的明顯優(yōu)勢和應(yīng)用潛力。基于MAD方法，Nielsen[14]又提出了迭代加權(quán)多元變化檢測(iterative weighted multivariate alteration detection，IRMAD)方法; Canty等[15]通過實(shí)驗(yàn)證實(shí)，該方法在干旱或沙漠地區(qū)以及場景中變化部分所占比例相對較大的多時(shí)相多波段遙感圖像變化檢測中，效果優(yōu)于MAD方法； 但在變化部分所占比例相對較小的非干旱或沙漠地區(qū)，并不能達(dá)到很好的效果。為此，本文提出了一種基于選權(quán)迭代估計(jì)(iterative estimation with weight selection，IEWS)與非監(jiān)督分類(unsupervised classification，UC)的多光譜圖像變化檢測方法。通過借鑒IEWS的思想[16]，并以類似于IRMAD的迭代模式進(jìn)行回歸估計(jì)，剔除掉大部分的非變化信息，得到初始變化信息； 對于一個(gè)非單一場景，由于地物的多樣性，2個(gè)時(shí)相圖像中不同地物的非變化信息之間存在不同的線性關(guān)系，因而通過基于初始變化信息的分類處理，以及對不同類別的IEWS，剔除殘存的非變化信息，得到最終的變化信息。利用真實(shí)TM圖像處理的結(jié)果驗(yàn)證了本文方法的有效性，所得到的變化信息在空間位置上同該區(qū)域相應(yīng)時(shí)間段內(nèi)土地覆蓋類型變化情況具有很好的一致性； 同時(shí)與采用MAD，IRMAD方法的變化檢測結(jié)果相比較，表明本文方法對相對較小的變化信息具有更好的變化檢測能力。

1 算法的基本原理

借鑒IEWS思想，本文先用最小二乘法進(jìn)行參數(shù)估計(jì)； 每次估計(jì)后，用差值D計(jì)算下一次迭代中圖像每一個(gè)像元X的權(quán)值P，使非線性變化像元的權(quán)值越來越小，直到趨近于0，這樣發(fā)生非線性變化的像元就被排除在估計(jì)樣本之外； 最后獲得的估計(jì)參數(shù)將不受變化信息的影響。IEWS的計(jì)算過程如圖1所示。

圖1 選權(quán)迭代估計(jì)的計(jì)算步驟

通過一次IEWS獲得的變化信息將是圖像中可能發(fā)生變化的所有像元值。由于地物類型的多樣性，并不能一次將非變化信息完全剔除。為了從初始變化信息中找出更為精確的土地覆蓋變化信息，需要對差值D做進(jìn)一步處理。對于同類地物，其所在的前后2個(gè)時(shí)相的遙感圖像應(yīng)該具有很好的線性相關(guān)性，因而通過基于變化信息進(jìn)行的分類，并對每類像元進(jìn)行IEWS，則可以剔除掉每類像元中的線性變化部分。所有N種類型地物非線性變化像元合集即為2個(gè)時(shí)相圖像之間真正的變化信息。變化信息獲取過程如圖2所示。

圖2 基于選權(quán)迭代估計(jì)與非監(jiān)督分類的變化檢測流程

2 權(quán)函數(shù)的選取

如何確定權(quán)值是上述估計(jì)過程中需要解決的關(guān)鍵問題。本文利用類似于IRMAD的迭代模式[17]來確定每次迭代過程中的權(quán)值P。

2個(gè)時(shí)相圖像X,Y都是多波段圖像，其差值圖像D包含有多個(gè)圖像層。通常假設(shè)差值圖像中每一層圖像的像元值近似服從高斯分布，并且原始圖像的每個(gè)波段之間相互獨(dú)立，從而每個(gè)像元的D值的平方和歸一化方差將近似服從自由度與波段數(shù)n相等的Chi-square分布(即“卡方分布”)。在曲線擬合中，可使用Chi-square作為模擬擬合優(yōu)劣的一個(gè)指標(biāo)，Chi-square越小越好。

對于第j個(gè)像元，cj為Chi-square分布中的一個(gè)樣本，有

,

(1)

式中：n為波段數(shù)；σDi為第i個(gè)波段差值圖像Di的標(biāo)準(zhǔn)差。進(jìn)而可將像元不變化的概率Pr[5]表示為

Pr=1-Pχ2;n(c)

,

(2)

式中：χ2為自由度為n的Chi-square分布；c為變化的像元數(shù)，c值越小，不變化概率越大。基于此，迭代過程中將像元不變化概率Pr作為觀測值的權(quán)值，進(jìn)行加權(quán)迭代計(jì)算。這樣，變化信息的權(quán)值將逐漸趨近于零，含有變化信息的觀測值將被剔除，則最后獲得觀測值的權(quán)等于不包含變化信息的觀測值的權(quán)。通過IEWS后獲得的差值圖像即為2個(gè)時(shí)相圖像之間的初始變化信息。

3 分類方法的選擇與實(shí)現(xiàn)

對于分類方法的選擇，目前常用的遙感分類方法主要包括監(jiān)督分類與非監(jiān)督分類。由于監(jiān)督分類法需要事先確定訓(xùn)練樣本數(shù)據(jù)，而在很多情況下訓(xùn)練樣本難以獲取； 而非監(jiān)督分類直接根據(jù)待分類數(shù)據(jù)個(gè)體間的相似性測度進(jìn)行類別劃分，無需訓(xùn)練樣本數(shù)據(jù)。為了實(shí)現(xiàn)較高的自動(dòng)化程度，本文選用非監(jiān)督分類ISODATA算法進(jìn)行分類。ISODATA算法是在沒有先驗(yàn)知識(shí)的情況下對數(shù)據(jù)進(jìn)行分類，能夠吸取中間結(jié)果所得到的經(jīng)驗(yàn)，具有自組織性，能夠自動(dòng)地進(jìn)行類別的“合并”和“分裂”，其各個(gè)參數(shù)可在聚類調(diào)整中逐漸確定，并最終構(gòu)建所需要的判別函數(shù)，且運(yùn)算速度與精度均較好[18]。為此，本文選用ISODATA法作為分類算法。

4 實(shí)驗(yàn)結(jié)果與分析

本文選取香港維多利亞港周圍區(qū)域(覆蓋了九龍半島及香港島和青衣島的一部分)進(jìn)行研究。從20世紀(jì)90年代開始的香港填海造地等工程使該區(qū)域內(nèi)的土地覆蓋和利用均發(fā)生了很大的變化。為了盡量避免植被覆蓋、氣候變化等因素對真實(shí)變化的影響，用作變化檢測的2景圖像的成像時(shí)間處于不同年份的同一月份。2景圖像均為Landsat5 TM圖像，空間分辨率均為30 m，成像時(shí)間分別為1989年11月和2005年11月。選用除熱紅外波段(TM6)以外的6個(gè)可見光和近紅外波段，對原始圖像進(jìn)行幾何配準(zhǔn)后截取出360像元×380像元的區(qū)域作為研究對象(圖3)。

(a) 1989年11月20日(b) 2005年11月23日

對2個(gè)時(shí)相圖像運(yùn)用本文方法進(jìn)行變化檢測，通過一次IEWS得到的初始變化信息圖像如圖4(a)所示，顏色很暗或很亮的地方為變化區(qū)域。對初始變化信息進(jìn)行非監(jiān)督的ISODATA分類，并對每類像元單獨(dú)進(jìn)行估計(jì)，得到每類地物的變化信息，如圖4(b)所示。經(jīng)過分類后的再次IEWS，進(jìn)一步清楚地顯示出感興趣的變化信息，使得變化區(qū)域與非變化區(qū)域的對比更為明顯。

(a) 第一次估計(jì)結(jié)果(b) 最后估計(jì)結(jié)果

本文選用可分離性測度Jeffries-Matusita distance,即JM距離來評價(jià)圖4中2次估計(jì)結(jié)果的有效性。根據(jù)先驗(yàn)知識(shí)，選取了2個(gè)典型樣本區(qū)(圖4)，典型樣本區(qū)1在2個(gè)時(shí)相期間內(nèi)由水體變?yōu)殛懙兀瑢⑵渥鳛樽兓瘶颖荆?典型樣本區(qū)2在2個(gè)時(shí)相都為林地，故將其作為非變化樣本。2次估計(jì)結(jié)果在各波段上2類樣本間的JM距離如表1所示，除波段6外，最后估計(jì)結(jié)果的JM距離都大于第1次估計(jì)結(jié)果，這說明最后估計(jì)結(jié)果具有更好的可分離性。

表1 2次估計(jì)結(jié)果在各波段上2類樣本間的JM距離

將最后變化圖像分為2部分，圖5(a)為像元灰度值相對變大的情況，圖5(b)為像元灰度值相對變小的區(qū)域，分別對應(yīng)于第2時(shí)相遙感圖像相對第1時(shí)相遙感圖像變亮和變暗的區(qū)域。

(a) 灰度值相對變大區(qū)域(b) 灰度值相對變小區(qū)域

圖5(a)中的變化像元顯示出了1989—2005年間該區(qū)域因填海造地工程而新出現(xiàn)的填海土地、碼頭及橋梁等； 而在圖5(b)中，紅色圓圈標(biāo)記出的地塊，則顯示在該期間土地覆蓋情況的變化導(dǎo)致了遙感圖像亮度變暗。這表明本文方法可以分辨出每個(gè)波段上像元灰度值是相對變大還是變小了，有助于對變化情況做進(jìn)一步的分析。

對2個(gè)時(shí)相遙感圖像使用MAD方法得到6個(gè)檢測圖像(圖6)。可以看出，僅MAD5圖像集中反映出區(qū)域內(nèi)填海造地工程帶來的變化，其他圖像反映的則是少數(shù)次要變化和大量噪聲成分。

(a) MAD1(b) MAD2(c) MAD3

(d) MAD4(e) MAD5(f) MAD6

采用本文方法最后得到的變化像元單波段圖像與MAD5圖像使用相同方法進(jìn)行變化區(qū)域提取的變化圖像的對比如圖7所示。本文方法所得到的結(jié)果在表達(dá)影像細(xì)節(jié)方面(圖7中綠色圓圈所標(biāo)注的橋梁和碼頭等處)有更好的效果。

為進(jìn)一步說明本文方法對細(xì)微目標(biāo)的檢測能力，選取圖7(a)中間和右側(cè)綠色標(biāo)注出的2個(gè)區(qū)域作為典型區(qū)域，分別采用本文方法和MAD5圖像進(jìn)行比較。先對這2個(gè)典型區(qū)域進(jìn)行目視解譯獲取其真實(shí)變化情況，再將本文方法及MAD5圖像進(jìn)行閾值提取后的變化結(jié)果與真實(shí)變化情況進(jìn)行對比，結(jié)果如表2所示。表中MAD5的提取結(jié)果2是該方法誤檢率和漏檢率均為最低情況下的統(tǒng)計(jì)結(jié)果。將其與本文方法提取結(jié)果比較可知，本文方法有更低的誤檢率和漏檢率，能更好地區(qū)分變化與非變化信息。

(a) 變化區(qū)域提取圖像 (b) MAD5變化區(qū)域提取圖像

前已述及，IRMAD方法在干旱或沙漠地區(qū)以及場景中變化部分所占比例相對較大的多時(shí)相多波段遙感圖像變化檢測時(shí)，其性能優(yōu)于MAD方法，但在變化部分所占比例相對較小的非干旱或沙漠地區(qū)，卻不一定能達(dá)到很好的效果。采用IRMAD方法對

表2 典型區(qū)域提取結(jié)果對比

本文研究區(qū)檢測的變量圖像如圖8所示。目測便可發(fā)現(xiàn)，其變化檢測效果均遜于MAD方法和本文方法。

(a) IRMAD1(b) IRMAD2(c) IRMAD3

(d) IRMAD4(e) IRMAD5(f) IRMAD6

5 結(jié)論

1)本文方法利用第一次估計(jì)獲得的初步變化結(jié)果進(jìn)行分類，并對分類后的圖像再次進(jìn)行估計(jì)，能更為精確地剔除變化檢測結(jié)果中的非變化信息，保留變化信息。通過計(jì)算2次估計(jì)結(jié)果的JM距離，獲得的最終估計(jì)結(jié)果對變化信息和非變化信息具有更好的可分離性。

2)本文通過與MAD及IRMAD方法的實(shí)驗(yàn)結(jié)果相比較，證實(shí)了本文提出方法的有效性和優(yōu)越性。

3)今后還有待于針對不同的數(shù)據(jù)源和不同的地形地貌情況做進(jìn)一步的研究，以驗(yàn)證該方法在不同情況下的適用性。

參考文獻(xiàn)(References)：

[1] 趙英時(shí).遙感應(yīng)用分析原理與方法[M].北京:科學(xué)出版社,2003:241-242.

Zhao Y S.The Theory and Method of Remote Sensing Application and Analysis[M].Beijing:Science Press,2003:241-242.

[2] 張繼賢,楊貴軍.單一時(shí)相遙感數(shù)據(jù)土地利用與覆蓋變化自動(dòng)檢測方法[J].遙感學(xué)報(bào),2005(3):294-299.

Zhang J X,Yang G J.Automatic land use and land cover change detection with one temporary remote sensing image[J].Journal of Remote Sensing,2005(3):294-299.

[3] Sing A.Digital change detection techniques using remotely-sensed data[J].International Journal of Remote Sensing,1989,10:989-1003.

[4] Mouat D A,Mahin G C,Lancaster J.Remote sensing techniques in the analysis of change detection[J].Geocarto International,1993,8(2):39-50.

[5] Arzandeh S,Wang J.Monitoring the change of phragmites and distribution using satellite data[J].Canadian Journal of Remote Sensing,2003,29(1):24-35.

[6] Li D R.Remotely sensed images and GIS data fusion for automatic change detection[J].International Journal of Image and Data Fusion,2010,1(1):99-108.

[7] Kwarteng A Y,Chavez P S.Change detection study of Kuwait City and environs using multi-temporal Landsat Thematic Mapper data[J].International Journal of Remote Sensing,1998,19(3):1651-1661.

[8] 張輝,王建國.一種基于主分量分析的SAR圖像變化檢測算法[J].電子與信息學(xué)報(bào),2008,30(7):1728-1730.

Zhang H,Wang J G.A SAR image change detection algorithm based on principal component analysis[J].Journal of Electronics and Information Technology,2008,30(7):1728-1730.

[9] Nielsen A A,Conradsen K,Simpson J J.Multivariate alteration detection(MAD)and MAF postprocessing in multispectral,bitemporal image data:New approaches to change detection studies[J].Remote Sensing of Environment,1998,64(1):1-19.

[10]Canty M J,Nielsenb A A,Schmidtc M.Automatic radiometric normalization of multi-temporal satellite imagery[J].Remote Sensing of Environment,2004,91(3/4):441-451.

[11]Canty M J,Nielsen A A.Visualization and unsupervised classification of changes in multispectral satellite imagery[J].International Journal of Remote Sensing,2006,27(18):3961-3975.

[12]廖明生,朱攀,龔健雅.基于典型相關(guān)分析的多元變化檢測[J].遙感學(xué)報(bào),2000,4(3):197-201.

Liao M S,Zhu P,Gong J Y.Multivariate change detection based on canonical transformation[J].Journal of Remote Sensing,2000,4(3):197-201.

[13]陳壘,馬潤賡,申維.基于典型相關(guān)分析的遙感影像變化檢測[J].地質(zhì)通報(bào),2007,26(7):916-920.

Chen L,Ma R G,Shen W.Detection of remote sensing image alteration based on canonical correlation analysis[J].Geological Bulletin of China,2007,26(7):916-920.

[14]Nielsen A A.The regularized iteratively reweighted MAD method for change detection in multi- and hyperspectral data[J].IEEE Transactions on Image Processing,2007,16(2):463-478.

[15]Canty M J,Nielsen A A.Automatic radiometric normalization of multi-temporal satellite imagery with the iteratively re-weighted MAD transformation[J].Remote Sensing of Environment,2008,112(3):1025-1036.

[16]邱衛(wèi)寧,陶本藻,姚宜斌,等.測量數(shù)據(jù)處理理論與方法[M].武漢:武漢大學(xué)出版社,2008:123-145.

Qiu W N,Tao B Z,Yao Y B,et al.The Theory and Method of Surveying Data Processing[M].Wuhan:Wuhan University Press,2008:123-145.

[17]Marpu P R,Gamba P,Canty M J.Change detection using iteratively reweighted regression with neural networks[C]//2010 IEEE International Geoscience and Remote Sensing Symposium.Honolulu,HI:IEEE,2010:2563-2566.

[18]顧有林,韓幫春.快鳥數(shù)據(jù)在光學(xué)遙感器成像仿真中的應(yīng)用[J].系統(tǒng)仿真學(xué)報(bào),2008,20(22):6265-6267.

Gu Y L,Han B C.Application of high resolution quickbird satellite data in imaging simulation of optical remote sensor[J].Journal of System Simulation,2008,20(22):6265-6267.