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

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

0 引言

多光譜圖像的變化檢測是對不同時間獲取的多光譜遙感圖像進行定量分析，并確定圖像區(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)性等)，丹麥學者Nielsen等[9-11]率先提出了多元變化檢測(multivariate alteration detection，MAD)的概念和方法。廖明生等[12]和陳壘等[13]將這種方法應(yīng)用于多時相多波段遙感圖像變化檢測，實驗結(jié)果也顯示出了該方法的明顯優(yōu)勢和應(yīng)用潛力?；贛AD方法，Nielsen[14]又提出了迭代加權(quán)多元變化檢測(iterative weighted multivariate alteration detection，IRMAD)方法; Canty等[15]通過實驗證實，該方法在干旱或沙漠地區(qū)以及場景中變化部分所占比例相對較大的多時相多波段遙感圖像變化檢測中，效果優(yōu)于MAD方法； 但在變化部分所占比例相對較小的非干旱或沙漠地區(qū)，并不能達到很好的效果。為此，本文提出了一種基于選權(quán)迭代估計(iterative estimation with weight selection，IEWS)與非監(jiān)督分類(unsupervised classification，UC)的多光譜圖像變化檢測方法。通過借鑒IEWS的思想[16]，并以類似于IRMAD的迭代模式進行回歸估計，剔除掉大部分的非變化信息，得到初始變化信息； 對于一個非單一場景，由于地物的多樣性，2個時相圖像中不同地物的非變化信息之間存在不同的線性關(guān)系，因而通過基于初始變化信息的分類處理，以及對不同類別的IEWS，剔除殘存的非變化信息，得到最終的變化信息。利用真實TM圖像處理的結(jié)果驗證了本文方法的有效性，所得到的變化信息在空間位置上同該區(qū)域相應(yīng)時間段內(nèi)土地覆蓋類型變化情況具有很好的一致性； 同時與采用MAD，IRMAD方法的變化檢測結(jié)果相比較，表明本文方法對相對較小的變化信息具有更好的變化檢測能力。

1 算法的基本原理

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

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

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

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

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

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

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

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

,

(1)

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

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

,

(2)

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

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

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

4 實驗結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

前已述及，IRMAD方法在干旱或沙漠地區(qū)以及場景中變化部分所占比例相對較大的多時相多波段遙感圖像變化檢測時，其性能優(yōu)于MAD方法，但在變化部分所占比例相對較小的非干旱或沙漠地區(qū)，卻不一定能達到很好的效果。采用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)本文方法利用第一次估計獲得的初步變化結(jié)果進行分類，并對分類后的圖像再次進行估計，能更為精確地剔除變化檢測結(jié)果中的非變化信息，保留變化信息。通過計算2次估計結(jié)果的JM距離，獲得的最終估計結(jié)果對變化信息和非變化信息具有更好的可分離性。

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

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

參考文獻(References)：

[1] 趙英時.遙感應(yīng)用分析原理與方法[M].北京:科學出版社,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ù)據(jù)土地利用與覆蓋變化自動檢測方法[J].遙感學報,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].電子與信息學報,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].遙感學報,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ì)通報,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].武漢:武漢大學出版社,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ù)在光學遙感器成像仿真中的應(yīng)用[J].系統(tǒng)仿真學報,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.