SAR图像变换检测

1、题 目： 图像的变化检测 学 院： 电子工程学院 专 业： 智能科学与技术 学生姓名： * 学 号： 02115* 图像变化检测一、 研究背景及方法1.1 研究背景变化检测的实质是地表特征随时间变化发生的变化而引起的两个时期影像像元光谱响应的变化，所以利用不同时期的影像进行变化检测就能获得地物的变化信息。SAR图像的变化检测专指利用多时相获取的同一地表区域的SAR图像来确定和分析地表变化，能提供地物的空间展布及其变化的定性与定量信息。目前遥感图像变化检测已广泛应用于国民经济和国防建设的诸多领域。民用方面，变化检测主要用于资源和环境监测中的土地利用和土地覆盖变化、森林和植被变化、湿地变化、城区变

2、化、地形改变等变化信息获取；测绘中的地理空间数据更新；农业中的作物生长监测以及自然灾害中地震、洪水、泥石流和森林大火等灾情监测与评估。在军事上，变化检测主要用于战场信息动态感知、军事目标和兵力部署监测等。这些应用的需求促使了遥感图像变化检测技术的快速发展。1.2 现有的变化检测方法 变化检测技术从处理遥感数据的光谱波段数的多少可以分为多波段多时相遥感图像变化检测和单波段多时相遥感图像变化检测。因为随着数据量的增大，多波段多时相图像数据冗余过量，其变化检测技术与比单波段遥感图像检测技术相比更为复杂，因此目前大多数变化检测技术集中在单波段多时相遥感图像变化检测方面。对单波段多时相遥感图像已经提出的

3、多种不同变化检测方法。按照是否需要分类可分为直接变化检测和分类后变化检测；根据是否需要先验知识可分为监督变化检测和非监督变化检测；根据变化检测所选取的分析对象可分为面向像素、特征、对象的变化检测。二、 算法流程简介2.1 算法流程图预处理（配准和辐射校正）时相1图像I1时相2图像I2图像比较（对数比值法）RI多尺度二维平稳小波变换尺度积去噪尺度积去噪尺度积去噪小波逆变换小波逆变换小波逆变换PCA变换后第一主成分图作为新差异图阈值分割新差异图（FCM） layer1 layer2 layer3 layer42.2 算法介绍 本文采用基于多尺度积和PCA的SAR图像变化检测方法的流程图如上所示。首

4、先对输入的已校正的两时相SAR图像I1、I2进行对数比运算，得到一幅对数比值图RI，对RI进行二维平稳小波4层分解，对每一分解层都采用SSNF的方法进行多尺度积去噪；其次将每一分解层多尺度积去噪后的小波系数进行逆二维平稳小波变换，得到对数比值图RI不同尺度图像layer1layer4；然后对图像layer1layer4采用PCA变换方法，并将PCA变换后的第一主成份图作为新的差异图；最后用模糊c均值聚类的方法对差异图NRI进行分割，从而得到变化检测结果。2.2.1图像预处理预处理包括图像配准和辐射校正，主要目的是尽可能的减少图像错位以及辐射度变化带来的影响。图像配准将不同时相获得的图像进行几何

5、对应，用于消除非地物变化的几何畸变。辐射校正用于消除辐射失真，包括绝对辐射校正和相对辐射校正，其中绝对辐射校正需要通过己知的各种传感器和环境参数得到图像中地物的光谱反射值，相对辐射校正将某个时相遥感图像作为参照图，使相同地物在不同时相图像中具有一致的光谱反射值，该方法是一种克服辐射失真对变化检测影响的有效手段。2.2.2对数比值法普遍认为图像比值法比图像差值法对SAR图像中的乘性斑噪具有更好的抵抗能力。对于SAR图像，一般采用对数比值法，可以先将SAR图像中乘性斑噪转化成加性噪声，便于对SAR图像进行有效地去噪。RI=logI1I2=log I1-logI22.2.3多尺度积通过对数变换，将乘

6、性噪声转化为加性噪声，易于多尺度积去噪。根据小波分解的不同尺度间相关性，采用空间选择性的噪声滤波(Spatially Selective Noise Filter，SSNF)方法，该方法根据在相邻尺度下同一位置小波系数的相关性，在大的尺度下小波变换系数全部被保留，在较小尺度下只能保留高相关的系数，舍弃低相关系数，从而实现在有噪声信号中保留重要边缘和去噪的目的。SSNF算法是一个迭代方法，令W(i，j)为图像小波变换系数，具体步骤如下：(1) 定义尺度间小波系数的相关性系数为相邻尺度的系数直接相乘。corrmdi,j=mm+1wmdi,j m=1,2,L 其中m为当前层，d为水平、垂直、对角方向

7、，L为小波分解层数(2) 计算相关性系数和小波系数的能量，即pwmdi,j=i,jwmd2(i,j)pcorrmdi,j=i,jcorrmd2(i,j)(3) 归一化相关性系数Ncorrmdi,j=corrmdi,j*(pwmd/pcorrmd)(4) 比较Ncorrmdi,j和wmdi,jwmdi,j=wmdi,j，ifNcorrmdi,j>|wmdi,j|0， else图1 尺度二维平稳小波变换后的图像图2 尺度积去噪后的图像图3 小波逆变换后的图像2.2.4主成分分析（PCA）对上述得到的对数比值图的每个尺度去噪之后的图像进行PCA变换以消除尺度间的相关性。由于输入图像高度相关的结

8、构成分在PCA变换后的第一主成份图中的占主导成分，而且变化信息将在第一主成份图中增强，因此本文用PCA变换后的第一主成份作为新的差异图，其具体步骤如下：（1） 将上述得到的对数比值图每个尺度去噪之后的图像layer列向量化，并组成矩阵A=(A1，A2，AI)（2） 计算矩阵A的相关系数矩阵RR=r11r12r21r22r13r14r23r24r31r32r41r42r33r34r43r44其中矩阵元素rij=k=1nAk-Ai(Ak-Aj)k=1nAk-Ai2k=1nAk-Aj2（3） 计算相关系数矩阵尺的特征值和特征向量。（4） 图像序列经过PCA变换后得到新的图像Z=Au。（5） 将PCA

9、变换后的第一主成份图ZI作为新的差异图NRI并输出。图3 小波逆变换后的图像2.2.5模糊c均值聚类（FCM）对上述得到的新差异图NRI用模糊C均值进行聚类，分为变化类和非变化类。（1） 随机初始化聚类中心ci,i=1,2。（2） 确定隶属度矩阵Uuij=1k=1c(dijdkj)2（3） 计算代价函数，如果它小于某个确定的阈值，或其相对于上次价值函数的该变量小于某个阈值，则算法停止。否则，转（4）（4） 修正聚类中心,返回（2）。ci=j=1nuijmxjk=1c(dijdkj)2(m-1)三、 实验结果分析（1）墨西哥郊外正确率漏检数虚警数总错误数97.5460%57376966433 （

10、2）Ottawa地区水灾的RADARSAT SAR图像正确率漏检数虚警数总错误数89.9689%65746401173四、 实验代码%-主程序-%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%-比值图的构造-%RI=CreatRatioImage(I1,I2);%figure(1)%imshow(RI);title('比值图像');%-多尺度二维平稳小波变换-%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);%-变化检测的性能评估-%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);%-%-比值图的构造-%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ºÍ