《智能计算导论课程设计-基于模糊贴近度的SAR图像变化检测.doc》由会员分享，可在线阅读，更多相关《智能计算导论课程设计-基于模糊贴近度的SAR图像变化检测.doc（18页珍藏版）》请在三一文库上搜索。
1、智能计算导论课程设计 题目：基于模糊贴近度的SAR图像变化检测 班级： 姓名： 学号：【摘要】：图像变化检测技术旨在检测不同时间段的同一场景图像之间发生的变化。图像变化检测技术主要依赖于辐射值或局部纹理的变化。这些变化可能是由于地表覆盖的真实变化引起的，或者是由照射角、大气条件、传感器精度。地面湿度等条件变化引起的。变化检测的基本前提是相对于由一些随机因素引起的变化由对象本身变化引起的辐射值或局部纹理的变化是可分的。合成孔径雷达具有全天候、全天时的特点，可方便地获得同一地区不同时间段的图像，于是SAR图像变化检测是它的一个重要应用领域。SAR图像变化检测是遥感技术的主要应用之一，它通过对不同时
2、期图像的比较分析，根据图像之间的差异来得到所需要的变化信息。SAR图像变化检测已经应用于很多方面，如对人工检测比较困难的热带雨林、沙漠等自然条件恶劣的地区进行监测，以了解生态环境变化的情况；对农田进行监测，分析农作物生长状况；对军事目标进行监测，了解兵力部署、军事调动等情报；对城区环境进行监控，合理规划城市布局，土地使用的管理和规范，非法占地和违章建筑的监控等。 随着研究环境和持续发展相关的许多地物变迁现象都是动态变化过程，随着遥感观测技术及SAR成像技术的发展，获取多波段、多极化、多时段的SAR图像的能力已大大提高，可为这些动态过程问题的研究提供丰富的信息源。但是由于地球系统的开放性、复杂性
3、和不确定性，如何有效地提取所需的信息，揭示出其变化等规律，是当前急需解决的问题，通过对不同时段的多通道(多波段或多极化等)的SAR图像之间的变化检测来获取地表覆盖的变化情况便成为一个日益重要的课题研究方向之一。 本文研究了基于模糊贴近度的图像进行变化检测，所提供的实验结果和验证算法的有效性，并显示了SAR图象在变化检测方面有良好的应用前景。【关键词】：SAR图象 模糊贴近度 变化检测 1 引言遥感变化检测是指通过对不同时期同一区域的遥感图像进行比较分析，根据图像之间的差异得到所需要的地物或目标的变化信息。现代遥感技术的飞速发展为变化检测提供了一种便捷的途径，遥感数据成为变化检测的主要数据源。与
4、可见光和红外遥感相比，微波遥感具有无可比拟的优点：微波能穿透云雾、雨雪，具有全天候、全天时的工作能力。二，微波对地物有一定穿透能力。三，采用侧视方式成像，覆盖面积大。正是这些优点，使得SAR 图像日益成为变化检测的重要数据源。SAR 变化检测技术的需求日益广泛。目前，全球坏境变化加剧，城市急速发展，洪水、地震等自然灾害时有发生，这些都需要及时掌握相关动态信息，为相关决策部门提供支持，而SAR 的种种优点为快速响应提供了技术支持和应急保障。2 变化检测的预处理2.1 图像配准图像的配准精度对变化检测的精度影响很大，如果图像没有较高精度的配准，则图像的错位会带来大量的虚警。因此，针对这一方面，Da
5、i 和Khorran1得出的结论为：好的变化检测结果需要高精度的配准算法。Jeannie Moulton 2等针对多时相SAR 图像的平移及偏转误差提出一种解决方法。意大利的Paolo 3针对高分辨SAR 影像配准的困难，提出基于象元与特征相结合的方法进行未精确配准的变化检测研究。2.2 噪声抑制基于 SAR 图像的变化检测对斑点的抑制要求很高：降班算法要既能有效降斑，还能保持图像的细节信息。因此，需要采用性能良好的自适应降斑方法来完成。相干斑抑制方法可以分为两类：一类是成像前的多视平滑处理，但同时也降低了图像的空间分辨率，所以目前常用的方法为成像后的去相干斑噪声的滤波技术。滤波方法主要有：（
6、1）传统方法，如均值滤波、中值滤波等；（2）模型方法，假定静态的噪声模型，采用相应的滤波器进行处理，如Kalman 滤波和Lee 滤波；（3）几何滤波方法，如Gamma MAP 滤波；（4）基于局域统计自适应滤波，这种方法能在平滑噪声的同时较好的保持边缘信息，是目前的SAR 变化检测常用的方法。本次实验采用的方法是frost滤波，Frost滤波器是特定大小窗口的像素值和指数脉冲响应m卷积的Wiener自适应滤波器：其中K是滤波器参数，t0代表中心像素的位置，|t|是距t0的距离。Frost滤波器采用的斑点噪声模型采用的形式如下：这里hij是系统响应函数，“*”为卷积算子。在通常的应用中，一般假
7、定hij为delta函数，最小均方滤波器形式如下：这里t对应于空间域中像素之间的距离。选择脉冲函数m(t)，使下式最小：按照频率域中Wiener滤波器的推导，就可以找到上式的解。 Frost滤波对于斑点噪声的抑制效果很好，同时对边缘、纹理和细节特征也有较好的保留作用，从实验结构可以证明这一点。3 常用的变化检测算法SAR 图像变化检测是研究不同时期的同一场景图像之间发生的变化。归结起来，近些年逐步形成的SAR 图像变化检测方法可以分为：(1)基于简单代数运算的变化检测，经典的变化检测方法包括图像差值法、图像比值法、对数比值法；(2)基于图像变换的变换检测，经典的变化检测方法包括主成分分析、变化
8、向量分析法、相关分析法图像变换法；(3)基于图像分类的变化检测方法。3.1 代数运算法早期变化检测的主要方法是利用图像差值法，主要利用两幅图像对应像素点的灰度差值作为变化判别依据，该方法简单、直接，但是容易受成像质量、噪声等的影响。图像差值法应用于海岸线环境、森林变化、沙漠化等。与差值法相比，比值差异图像与原始SAR 图像的平均强度值无关，对于SAR 图像的乘性噪声不敏感，而且减少了辐射定标误差的影响。图像比值法对于SAR 图像上的乘性噪声是不敏感的，被广泛用于城区变化检测。对数比值法，在得到对应像素的比值后，再取其对数，能压缩图像的变化范围，将乘性噪声转换为加性噪声。代数运算法的优点是直接、
9、简单，主要缺点在于(1)没有考虑多时相图像之间的相关性，简单计算可能出现虚警率；(2)变化阈值确定困难。3.2 图像变换法Lillesand提出的主分量法，又称K-L 变换，是一种经典的数学变换方法，把原来多个波段中的有用信息集中到互不相关的新成分图像中，达到冗余压缩和信息集中的目的。在进行主成分变换时，利用协方差矩阵和相关矩阵得到的主分量是不同的，由相关矩阵推导的主分量变换对于多时相分析是尤其有用的。Malila 提出了变化向量分析法(Change Vector Analysis,CVA),Johnson等详细阐述了这种方法的具体内容、优缺点。变化向量法是一种多变量的方法，描述从第一时间到第
10、二时间的光谱变化的强度和方向。如果变化向量的幅值超过给定的门限，则判定该像素发生变化，变化向量的方向包含变化类型信息。该方法可用在多通道极化SAR 图像的变化检测，或者用于多特征，如空间结构特征，纹理特征等分析。相关分析法，斑点的时相去相关提供了地表在结构或介电特性上可能的变化信息，并且这种变化检测不依赖于定标精度。实际应用中，通常使用一般的强度图像来计算斑点相关系数，不考虑相位信息。3.3 分类法分类后比较法，在70 年代末开始应用于Landsat 卫星影像的变化检测。分类后比较法是对每幅图像单独进行分类，然后对图像的分类结果图进行比较，以确定变化的类别和区域。如果对应像素的分类类别相同，则
11、认为该像素没有发生变化，否则认为该像素发生了变化。使用分类法的优点主要有：(1)可以给出区域变化的类型、位置、数量等有关地物性质信息；(2)受图像配准的影响要小。此方法的主要缺点是它受分类器的影响很大，不同的分类器能够决定变化检测的精度高低。同时分类法，将多时相图像构成一副复合影像进行分类，得到的分类结果图中的每一类代表一种变化类型。此方法能够减少分类时间，但是难于标记变化类别，而且分类更加复杂，对训练样本的选取相当高。随着计算机技术的发展，解决复杂问题的各种数学方法与专家知识被逐渐加入。人工神经网络、向量机、专家智能等方法成为热点。变化检测的各种方法各有优缺点，各有适用范围，实际应用中，哪种
12、好坏没有定论。区域不同、需要不同，方法的选择也不同。3.3 模糊差异图 使用的方法是基于模糊贴近度的差异图构造的方法.对已配准和校正的两时相SAR图像T1和T2，根据SAR图像的特点，采用如下的相似度公式来构造差异图（difference image，简记：DI）：其中N(x)表示位置x上的一个邻域的位置指示集，即取其空间上的一个邻域，把由邻域元素构成的集合看作是模糊向量，用上式度量两模糊向量间的相似度。若得到上式的值越大，则该位置上的像素点越接近于未变化类，反之，为变化类。利用此方法得到的差异图,在一定范围内解决了因图像辐射差异、配准误差和差异图像分类门限的选取造成的难以从检测信息中提取出关
13、键的变化的问题.4 变化信息图像后处理4.1双边滤波 根据以上变化检测方法对两个时相的遥感影像进行变化信息的提取，进行变化信息提取后，由实际野外观测提供目标区变化的真实信息，并借助相关资料和专题资料对变化信息进行后处理，归并小图斑，消除提取出的变化信息中的伪变化信息。采取的方法是对差异图进行双边滤波，以去除变化不明显的区域。由于双边滤波对于边缘有着较好的保留性，滤波后的图像边缘没有明显的模糊。4.2 FCM均值聚类在对差异图进行双边滤波后，得到较平滑的滤波图像，但对于来说，有用信息（两幅图像变化区域）还不够明显，所以在这里利用聚类分割图像的办法，将滤波图像进行特征突出。5实验流程及程序流程图：
14、原程序：（1）frost滤波程序%Frost filter for speckle noise reduction function ft=frost(I)% I is the noisy input imageticx y z=size(I);I=double(I);K=1;N=I;for i=1:x for j=1:y if (i1 & i1 & j=IM2(i,j) tmp1(i,j)=IM2(i,j); tmp2(i,j)=IM1(i,j); else tmp1(i,j)=IM1(i,j); tmp2(i,j)=IM2(i,j); end endendmintmp=zeros(A1-2
15、,B1-2);maxtmp=zeros(A1-2,B1-2); for i=1:1:A1 for j=1:1:B1 for m=-2:2 %3*3的滑窗 for n=-2:2 if (i+m)0 & (i+m)0 &(j+n) Image img % - winsize: spatial filter width % - sigma for intensity diff gaussain filter % - sigma for spatial filter = winsize/6 % Output - Filtered Image % Author: Amit Agrawal, 2004 %
16、 should be odd winsize = round( (winsize-1)/2 )*2 + 1; H,W = size(img); %Gaussian spatial filter g_filter = fspecial(gaussian,winsize,winsize/6); padnum = (winsize-1)/2; A = padarray(img,padnum padnum,replicate,both); img1 = zeros(size(img); for jj = padnum+1:(padnum+1+H-1) for kk = padnum+1:(padnum
17、+1+W-1) % Get a local neighborhood imgwin = A(jj-padnum:jj+padnum,kk-padnum:kk+padnum); % Find weights according to intensity diffs Wwin = exp(-abs(imgwin - imgwin(padnum+1,padnum+1)/sigma2); % Find composite filter newW = Wwin.*g_filter; t = sum(sum(newW); if(t0) newW = newW/t; end img1(jj-padnum,k
18、k-padnum) = sum(sum(imgwin.*newW); end end（4）聚类（fcm）程序function IX2 = fcm(IM);%IM是输入的原图像%IX2是分类结果maxX,maxY=size(IM);IMM=cat(3,IM,IM,IM);%初始化聚类中心（3类）cc1=8;cc2=100;cc3=200;ttFcm=0;while(ttFcm15) ttFcm=ttFcm+1 c1=repmat(cc1,maxX,maxY); c2=repmat(cc2,maxX,maxY); c3=repmat(cc3,maxX,maxY); c=cat(3,c1,c2,c3
19、); ree=repmat(0.000001,maxX,maxY); ree1=cat(3,ree,ree,ree); distance=IMM-c; distance=distance.*distance+ree1; daoShu=1./distance; daoShu2=daoShu(:,:,1)+daoShu(:,:,2)+daoShu(:,:,3);%计算隶属度u distance1=distance(:,:,1).*daoShu2; u1=1./distance1; distance2=distance(:,:,2).*daoShu2; u2=1./distance2; distan
20、ce3=distance(:,:,3).*daoShu2; u3=1./distance3;%计算聚类中心z ccc1=sum(sum(u1.*u1.*IM)/sum(sum(u1.*u1); ccc2=sum(sum(u2.*u2.*IM)/sum(sum(u2.*u2); ccc3=sum(sum(u3.*u3.*IM)/sum(sum(u3.*u3); tmpMatrix=abs(cc1-ccc1)/cc1,abs(cc2-ccc2)/cc2,abs(cc3-ccc3)/cc3; pp=cat(3,u1,u2,u3); for i=1:maxX for j=1:maxY if max(p
21、p(i,j,:)=u1(i,j) IX2(i,j)=1; elseif max(pp(i,j,:)=u2(i,j) IX2(i,j)=2; else IX2(i,j)=3; end end end%判结束条件 if max(tmpMatrix)0.0001 break; else cc1=ccc1; cc2=ccc2; cc3=ccc3; end for i=1:maxX for j=1:maxY if IX2(i,j)=3 IMMM(i,j)=240; elseif IX2(i,j)=2 IMMM(i,j)=130; else IMMM(i,j)=20; end end endendfor
22、i=1:maxX for j=1:maxY if IX2(i,j)=3 IMMM(i,j)=240; elseif IX2(i,j)=2 IMMM(i,j)=20; else IMMM(i,j)=20; end endend%显示最终分类结果IMMM=uint8(IMMM);IX2=IMMM;end（5）主程序image1_bmp= imread(1997.7.24.bmp); %读第一副图image1_uint8 = image1_bmp(:,:,1); %取出灰度图frost1_uint8=frost(image1_uint8); %进行frost滤波image1_double = dou
23、ble(frost1_uint8)+1; %uint8转换为double格式image1 = image1_double/255; %归一化image2_bmp= imread(1998.10.24.bmp); %读第二副图image2_uint8 = image2_bmp(:,:,1); %取出灰度图frost2_uint8=frost(image2_uint8); %进行frost滤波image2_double = double(frost2_uint8)+1; %uint8转换为double格式image2 = image2_double/255; %归一化di_image=di(ima
24、ge1,image2); %构造模糊差异图bilateral_di_image=bilateral(di_image,36,6); %双边滤波fcm_image=fcm(bilateral_di_image); %聚类（FCM）figure;subplot(3,3,1),imshow(image1_uint8); title(原图1); %显示原图1subplot(3,3,2),imshow(frost1_uint8); title(1的frost滤波图); %显示图1的frost滤波图subplot(3,3,3),imshow(image2_uint8); title(原图2); %显示原图
25、2subplot(3,3,4),imshow(frost2_uint8); title(2的frost滤波图); %显示图2的frost滤波图subplot(3,3,5),imshow(image1); title(第1幅归一化图); %显示第一幅归一化图subplot(3,3,6),imshow(image2); title(第2幅归一化图); %显示第二幅归一化图subplot(3,3,7),imshow(di_image); title(模糊差异图);%模糊差异图subplot(3,3,8),imshow(bilateral_di_image); title(双边滤波图);%双边滤波图subplot(3,3,9),imshow(fcm_image); title(聚类图);%聚类图6实验结果实验一：实验二 实验三7总结 本次实验是在大量研究各种图像处理办法后总结的一种较优的图像处理办法。主要流程就是在读入图像后，对图像进行预处理，如frost滤波、图像归一化等。然后是利用模糊贴近度算法获得模糊差异图，其优点是能较好的处理效果，如较低的出错率和漏检率等。接着利用双边滤波和均值聚类（fcm）进行处理，并得到较好的图像处理结果。然而由于时间仓促，并没有详细考虑阈值设定等问题。