SAR圖像變換檢測(cè)

1、題 目： 圖像的變化檢測(cè) 學(xué) 院： 電子工程學(xué)院 專 業(yè)： 智能科學(xué)與技術(shù) 學(xué)生姓名： * 學(xué) 號(hào)： 02115* 圖像變化檢測(cè)一、 研究背景及方法1.1 研究背景變化檢測(cè)的實(shí)質(zhì)是地表特征隨時(shí)間變化發(fā)生的變化而引起的兩個(gè)時(shí)期影像像元光譜響應(yīng)的變化，所以利用不同時(shí)期的影像進(jìn)行變化檢測(cè)就能獲得地物的變化信息。SAR圖像的變化檢測(cè)專指利用多時(shí)相獲取的同一地表區(qū)域的SAR圖像來(lái)確定和分析地表變化，能提供地物的空間展布及其變化的定性與定量信息。目前遙感圖像變化檢測(cè)已廣泛應(yīng)用于國(guó)民經(jīng)濟(jì)和國(guó)防建設(shè)的諸多領(lǐng)域。民用方面，變化檢測(cè)主要用于資源和環(huán)境監(jiān)測(cè)中的土地利用和土地覆蓋變化、森林和植被變化、濕地變化、城區(qū)變

2、化、地形改變等變化信息獲取；測(cè)繪中的地理空間數(shù)據(jù)更新；農(nóng)業(yè)中的作物生長(zhǎng)監(jiān)測(cè)以及自然災(zāi)害中地震、洪水、泥石流和森林大火等災(zāi)情監(jiān)測(cè)與評(píng)估。在軍事上，變化檢測(cè)主要用于戰(zhàn)場(chǎng)信息動(dòng)態(tài)感知、軍事目標(biāo)和兵力部署監(jiān)測(cè)等。這些應(yīng)用的需求促使了遙感圖像變化檢測(cè)技術(shù)的快速發(fā)展。1.2 現(xiàn)有的變化檢測(cè)方法 變化檢測(cè)技術(shù)從處理遙感數(shù)據(jù)的光譜波段數(shù)的多少可以分為多波段多時(shí)相遙感圖像變化檢測(cè)和單波段多時(shí)相遙感圖像變化檢測(cè)。因?yàn)殡S著數(shù)據(jù)量的增大，多波段多時(shí)相圖像數(shù)據(jù)冗余過量，其變化檢測(cè)技術(shù)與比單波段遙感圖像檢測(cè)技術(shù)相比更為復(fù)雜，因此目前大多數(shù)變化檢測(cè)技術(shù)集中在單波段多時(shí)相遙感圖像變化檢測(cè)方面。對(duì)單波段多時(shí)相遙感圖像已經(jīng)提出的

3、多種不同變化檢測(cè)方法。按照是否需要分類可分為直接變化檢測(cè)和分類后變化檢測(cè)；根據(jù)是否需要先驗(yàn)知識(shí)可分為監(jiān)督變化檢測(cè)和非監(jiān)督變化檢測(cè)；根據(jù)變化檢測(cè)所選取的分析對(duì)象可分為面向像素、特征、對(duì)象的變化檢測(cè)。二、 算法流程簡(jiǎn)介2.1 算法流程圖預(yù)處理（配準(zhǔn)和輻射校正）時(shí)相1圖像I1時(shí)相2圖像I2圖像比較（對(duì)數(shù)比值法）RI多尺度二維平穩(wěn)小波變換尺度積去噪尺度積去噪尺度積去噪小波逆變換小波逆變換小波逆變換PCA變換后第一主成分圖作為新差異圖閾值分割新差異圖（FCM） layer1 layer2 layer3 layer42.2 算法介紹 本文采用基于多尺度積和PCA的SAR圖像變化檢測(cè)方法的流程圖如上所示。首

4、先對(duì)輸入的已校正的兩時(shí)相SAR圖像I1、I2進(jìn)行對(duì)數(shù)比運(yùn)算，得到一幅對(duì)數(shù)比值圖RI，對(duì)RI進(jìn)行二維平穩(wěn)小波4層分解，對(duì)每一分解層都采用SSNF的方法進(jìn)行多尺度積去噪；其次將每一分解層多尺度積去噪后的小波系數(shù)進(jìn)行逆二維平穩(wěn)小波變換，得到對(duì)數(shù)比值圖RI不同尺度圖像layer1layer4；然后對(duì)圖像layer1layer4采用PCA變換方法，并將PCA變換后的第一主成份圖作為新的差異圖；最后用模糊c均值聚類的方法對(duì)差異圖NRI進(jìn)行分割，從而得到變化檢測(cè)結(jié)果。2.2.1圖像預(yù)處理預(yù)處理包括圖像配準(zhǔn)和輻射校正，主要目的是盡可能的減少圖像錯(cuò)位以及輻射度變化帶來(lái)的影響。圖像配準(zhǔn)將不同時(shí)相獲得的圖像進(jìn)行幾何

5、對(duì)應(yīng)，用于消除非地物變化的幾何畸變。輻射校正用于消除輻射失真，包括絕對(duì)輻射校正和相對(duì)輻射校正，其中絕對(duì)輻射校正需要通過己知的各種傳感器和環(huán)境參數(shù)得到圖像中地物的光譜反射值，相對(duì)輻射校正將某個(gè)時(shí)相遙感圖像作為參照?qǐng)D，使相同地物在不同時(shí)相圖像中具有一致的光譜反射值，該方法是一種克服輻射失真對(duì)變化檢測(cè)影響的有效手段。2.2.2對(duì)數(shù)比值法普遍認(rèn)為圖像比值法比圖像差值法對(duì)SAR圖像中的乘性斑噪具有更好的抵抗能力。對(duì)于SAR圖像，一般采用對(duì)數(shù)比值法，可以先將SAR圖像中乘性斑噪轉(zhuǎn)化成加性噪聲，便于對(duì)SAR圖像進(jìn)行有效地去噪。RI=logI1I2=log I1-logI22.2.3多尺度積通過對(duì)數(shù)變換，將乘

6、性噪聲轉(zhuǎn)化為加性噪聲，易于多尺度積去噪。根據(jù)小波分解的不同尺度間相關(guān)性，采用空間選擇性的噪聲濾波(Spatially Selective Noise Filter，SSNF)方法，該方法根據(jù)在相鄰尺度下同一位置小波系數(shù)的相關(guān)性，在大的尺度下小波變換系數(shù)全部被保留，在較小尺度下只能保留高相關(guān)的系數(shù)，舍棄低相關(guān)系數(shù)，從而實(shí)現(xiàn)在有噪聲信號(hào)中保留重要邊緣和去噪的目的。SSNF算法是一個(gè)迭代方法，令W(i，j)為圖像小波變換系數(shù)，具體步驟如下：(1) 定義尺度間小波系數(shù)的相關(guān)性系數(shù)為相鄰尺度的系數(shù)直接相乘。corrmdi,j=mm+1wmdi,j m=1,2,L 其中m為當(dāng)前層，d為水平、垂直、對(duì)角方向

7、，L為小波分解層數(shù)(2) 計(jì)算相關(guān)性系數(shù)和小波系數(shù)的能量，即pwmdi,j=i,jwmd2(i,j)pcorrmdi,j=i,jcorrmd2(i,j)(3) 歸一化相關(guān)性系數(shù)Ncorrmdi,j=corrmdi,j*(pwmd/pcorrmd)(4) 比較Ncorrmdi,j和wmdi,jwmdi,j=wmdi,j，ifNcorrmdi,j>|wmdi,j|0， else圖1 尺度二維平穩(wěn)小波變換后的圖像圖2 尺度積去噪后的圖像圖3 小波逆變換后的圖像2.2.4主成分分析（PCA）對(duì)上述得到的對(duì)數(shù)比值圖的每個(gè)尺度去噪之后的圖像進(jìn)行PCA變換以消除尺度間的相關(guān)性。由于輸入圖像高度相關(guān)的結(jié)

8、構(gòu)成分在PCA變換后的第一主成份圖中的占主導(dǎo)成分，而且變化信息將在第一主成份圖中增強(qiáng)，因此本文用PCA變換后的第一主成份作為新的差異圖，其具體步驟如下：（1） 將上述得到的對(duì)數(shù)比值圖每個(gè)尺度去噪之后的圖像layer列向量化，并組成矩陣A=(A1，A2，AI)（2） 計(jì)算矩陣A的相關(guān)系數(shù)矩陣RR=r11r12r21r22r13r14r23r24r31r32r41r42r33r34r43r44其中矩陣元素rij=k=1nAk-Ai(Ak-Aj)k=1nAk-Ai2k=1nAk-Aj2（3） 計(jì)算相關(guān)系數(shù)矩陣尺的特征值和特征向量。（4） 圖像序列經(jīng)過PCA變換后得到新的圖像Z=Au。（5） 將PCA

9、變換后的第一主成份圖ZI作為新的差異圖NRI并輸出。圖3 小波逆變換后的圖像2.2.5模糊c均值聚類（FCM）對(duì)上述得到的新差異圖NRI用模糊C均值進(jìn)行聚類，分為變化類和非變化類。（1） 隨機(jī)初始化聚類中心ci,i=1,2。（2） 確定隸屬度矩陣Uuij=1k=1c(dijdkj)2（3） 計(jì)算代價(jià)函數(shù)，如果它小于某個(gè)確定的閾值，或其相對(duì)于上次價(jià)值函數(shù)的該變量小于某個(gè)閾值，則算法停止。否則，轉(zhuǎn)（4）（4） 修正聚類中心,返回（2）。ci=j=1nuijmxjk=1c(dijdkj)2(m-1)三、 實(shí)驗(yàn)結(jié)果分析（1）墨西哥郊外正確率漏檢數(shù)虛警數(shù)總錯(cuò)誤數(shù)97.5460%57376966433 （

10、2）Ottawa地區(qū)水災(zāi)的RADARSAT SAR圖像正確率漏檢數(shù)虛警數(shù)總錯(cuò)誤數(shù)89.9689%65746401173四、 實(shí)驗(yàn)代碼%-主程序-%clear;clc;I1=imread('img1.bmp');I2=imread('img2.bmp');Standard_Image=imread('cankao.bmp');I1=rgb2gray(I1);I2=rgb2gray(I2);Standard_Image=rgb2gray(Standard_Image);m n=size(I1);for i=1:m for j=1:n if(I1(i,

11、j)=0) I1(i,j)=1; end if(I2(i,j)=0) I2(i,j)=1; end endend%-比值圖的構(gòu)造-%RI=CreatRatioImage(I1,I2);%figure(1)%imshow(RI);title('比值圖像');%-多尺度二維平穩(wěn)小波變換-%LL,LH,HL,HH=swt2(RI,4,'db4');layer1=LL;layer2=LL;layer3=LL;layer4=LL;layer1(:,:,1)=LL(:,:,1);layer1(:,:,2)=LH(:,:,1);layer1(:,:,3)=HL(:,:,1);

12、layer1(:,:,4)=HH(:,:,1);layer2(:,:,1)=LL(:,:,2);layer2(:,:,2)=LH(:,:,2);layer2(:,:,3)=HL(:,:,3);layer2(:,:,4)=HH(:,:,2);layer3(:,:,1)=LL(:,:,3);layer3(:,:,2)=LH(:,:,3);layer3(:,:,3)=HL(:,:,3);layer3(:,:,4)=HH(:,:,3);layer4(:,:,1)=LL(:,:,4);layer4(:,:,2)=LH(:,:,4);layer4(:,:,3)=HL(:,:,4);layer4(:,:,4

13、)=HH(:,:,4);%imageshow(layer1,layer2,layer3,layer4,2);%-尺度積去噪-%tlayer1 tlayer2 tlayer3 tlayer4=ssnf(layer1,layer2,layer3,layer4);%imageshow(tlayer1,tlayer2,tlayer3,tlayer4,3);%-小波逆變換-%layer1=iswt2(tlayer1(:,:,1),tlayer1(:,:,2),tlayer1(:,:,3),tlayer1(:,:,4),'db4');layer2=iswt2(tlayer2(:,:,1),

14、tlayer2(:,:,2),tlayer2(:,:,3),tlayer2(:,:,4),'db4');layer3=iswt2(tlayer3(:,:,1),tlayer3(:,:,2),tlayer3(:,:,3),tlayer3(:,:,4),'db4');layer4=iswt2(tlayer4(:,:,1),tlayer4(:,:,2),tlayer4(:,:,3),tlayer4(:,:,4),'db4');%figure(4)%subplot(1,4,1);imshow(layer1);%subplot(1,4,2);imshow(

15、layer2);%subplot(1,4,3);imshow(layer3);%subplot(1,4,4);imshow(layer4);%-PCA變換-%NRI=MPCA(layer1,layer2,layer3);%figure(5)%imshow(NRI);title('PCA±ä»¯ºóµÄͼÏñ');%-閾值分割-%DI=CreatDifferenceImage(NRI);%DI=FCM(NRI,2);%-變化檢測(cè)的性能評(píng)估-%effective

16、_num error_num undetected_num misdetected_num=CDperformence(Standard_Image,DI);%-½á¹ûÏÔʾ-%figure(6);subplot(2,2,1);imshow(I1);title('2000Äê4ÔÂ');subplot(2,2,2);imshow(I2);title('2002Äê5ÔÂ');subplot(2,2,3

17、);imshow(Standard_Image);title('±ê×¼²îÒìͼ');subplot(2,2,4);imshow(DI);title('±¾Ëã·¨µÄ²îÒìͼ'); disp('ÕýÈ·ÂÊΪ£

18、6;');disp(effective_num/(m*n)*100);disp('ÓÐЧ¸öÊýΪ£º');disp(effective_num);disp('×Ü´íÎóÊýΪ£º');disp(error_num);disp('©¼ìÊý&

19、#206;ª£º');disp(undetected_num);disp('Ð龯ÊýΪ£º');disp(misdetected_num);%-%-比值圖的構(gòu)造-%function RI=CreatRatioImage(I1,I2)TI1=im2double(I1);TI2=im2double(I2);m n=size(I1);RI=zeros(m,n);for i=1:m for j=1:n RI(i,j)=log(TI1(i,j)-

20、log(TI2(i,j); endendend%-%-尺度積去噪聲-%function tlayer1 tlayer2 tlayer3 tlayer4=ssnf(layer1,layer2,layer3,layer4)%-¶¨Òå³ß¶È¼äС²¨ÏµÊýµÄÏà¹ØÐÔϵÊýÎ

21、70;ÏàÁڳ߶ȵÄϵÊýÖ±½ÓÏà³Ë-5corr1=layer1.*layer2;corr2=layer2.*layer3;corr3=layer3.*layer4;corr4=layer4.*layer4;%-¼ÆËãÏà¹ØÐÔϵÊ

22、;ýµÄÄÜÁ¿-%m n c=size(layer1);pw1=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pw1(1,d)=temp+layer1(i,j,d)*layer1(i,j,d); end endendm n c=size(layer2);pw2=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pw2(1,d)=temp+layer2(i,j,d)*layer2(i,j,d); end endendm n c

23、=size(layer3);pw3=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pw3(1,d)=temp+layer3(i,j,d)*layer3(i,j,d); end endendm n c=size(layer4);pw4=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pw4(1,d)=temp+layer4(i,j,d)*layer4(i,j,d); end endend%-¼ÆËãС²¨Ï

24、;µÊýµÄÄÜÁ¿-%m n c=size(corr1);pcorr1=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pcorr1(1,d)=temp+corr1(i,j,d)*corr1(i,j,d); end endendm n c=size(corr2);pcorr2=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pcorr2(1,d)=temp+corr2(i,j,d)*corr2(i,j,d

25、); end endendm n c=size(corr3);pcorr3=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pcorr3(1,d)=temp+corr3(i,j,d)*corr3(i,j,d); end endendm n c=size(corr4);pcorr4=zeros(1,c);for d=1:c temp=0; for i=1:m for j=1:n pcorr4(1,d)=temp+corr4(i,j,d)*corr4(i,j,d); end endend%-¹éÒ»»

26、;¯Ïà¹ØÐÔϵÊý-%Ncorr1=corr1;for d=1:4 Ncorr1(:,:,d)=corr1(:,:,d)*sqrt(pw1(1,d)/(pcorr1(1,d);endNcorr2=corr2;for d=1:4 Ncorr2(:,:,d)=corr2(:,:,d)*sqrt(pw2(1,d)/(pcorr2(1,d);endNcorr3=corr3;for d=1:4 Ncorr3(:,:,d)=corr3(:,:,d)*sqrt(pw3(1,d)/(pcorr3

27、(1,d);endNcorr4=corr4;for d=1:4 Ncorr4(:,:,d)=corr4(:,:,d)*sqrt(pw4(1,d)/(pcorr4(1,d);end%-±È½Ï-%m n c=size(layer1);tlayer1=layer1;for d=1:c for i=1:m for j=1:n if(abs(Ncorr1(i,j,d)<=abs(layer1(i,j,d) tlayer1(i,j,d)=0; end end endendm n c=size(layer2);tlayer2=layer2;for d=1:c

28、for i=1:m for j=1:n if(abs(Ncorr2(i,j,d)<=abs(layer2(i,j,d) tlayer2(i,j,d)=0; end end endendm n c=size(layer3);tlayer3=layer3;for d=1:c for i=1:m for j=1:n if(abs(Ncorr3(i,j,d)<=abs(layer3(i,j,d) tlayer3(i,j,d)=0; end end endendm n c=size(layer4);tlayer4=layer4;for d=1:c for i=1:m for j=1:n if

29、(abs(Ncorr4(i,j,d)<=abs(layer4(i,j,d) tlayer4(i,j,d)=0; end end endendend%-%-PCA-%function NRI=MPCA(layer1,layer2,layer3)row col=size(layer1);A1=reshape(layer1,row*col,1);A2=reshape(layer2,row*col,1);A3=reshape(layer3,row*col,1);%A4=reshape(layer4,row*col,1);%A=A1,A2,A3,A4;A=A1,A2,A3;%¼Æ

30、;ËãÏà¹ØϵÊý¾ØÕóRmm=ones(row*col,1);Ameans=A;for k=1:3 Ameans(:,k)=mean(A(:,k)*mm;endR=zeros(3,3);for i=1:3 for j=1:3 sum1=0;sum2=0;sum3=0; for k=1:3 sum1=sum1+(A(:,k)-Ameans(:,i)'*(A(:,k)-Ameans(:,j); sum2=sum2+(A(:,k)-Ameans(:,i)'*(A(:,k)-Ameans(:,i); sum3=sum3+(A(:,k)-Ameans(:,j)'*(A(:,k)-Ameans(:,j); end R(i,j)=sum1/(sqrt(sum2*sum3); endendV D=eig(R);%¼ÆËã¾ØÕóRµÄÌØÕ÷ÏòÁ¿¾ØÕóVºÍ