SAR图像变换检测.docx
- 文档编号:10903054
- 上传时间:2023-05-28
- 格式:DOCX
- 页数:22
- 大小:458.47KB
SAR图像变换检测.docx
《SAR图像变换检测.docx》由会员分享,可在线阅读,更多相关《SAR图像变换检测.docx(22页珍藏版)》请在冰点文库上搜索。
SAR图像变换检测
题目:
图像的变化检测
学院:
电子工程学院
专业:
智能科学与技术
学生姓名:
****
学号:
********
图像变化检测
一、研究背景及方法
1.1研究背景
变化检测的实质是地表特征随时间变化发生的变化而引起的两个时期影像像元光谱响应的变化,所以利用不同时期的影像进行变化检测就能获得地物的变化信息。
SAR图像的变化检测专指利用多时相获取的同一地表区域的SAR图像来确定和分析地表变化,能提供地物的空间展布及其变化的定性与定量信息。
目前遥感图像变化检测已广泛应用于国民经济和国防建设的诸多领域。
民用方面,变化检测主要用于资源和环境监测中的土地利用和土地覆盖变化、森林和植被变化、湿地变化、城区变化、地形改变等变化信息获取;测绘中的地理空间数据更新;农业中的作物生长监测以及自然灾害中地震、洪水、泥石流和森林大火等灾情监测与评估。
在军事上,变化检测主要用于战场信息动态感知、军事目标和兵力部署监测等。
这些应用的需求促使了遥感图像变化检测技术的快速发展。
1.2现有的变化检测方法
变化检测技术从处理遥感数据的光谱波段数的多少可以分为多波段多时相遥感图像变化检测和单波段多时相遥感图像变化检测。
因为随着数据量的增大,多波段多时相图像数据冗余过量,其变化检测技术与比单波段遥感图像检测技术相比更为复杂,因此目前大多数变化检测技术集中在单波段多时相遥感图像变化检测方面。
对单波段多时相遥感图像已经提出的多种不同变化检测方法。
按照是否需要分类可分为直接变化检测和分类后变化检测;根据是否需要先验知识可分为监督变化检测和非监督变化检测;根据变化检测所选取的分析对象可分为面向像素、特征、对象的变化检测。
二、算法流程简介
2.1算法流程图
layer1layer2layer3layer4
2.2算法介绍
本文采用基于多尺度积和PCA的SAR图像变化检测方法的流程图如上所示。
首先对输入的已校正的两时相SAR图像I1、I2进行对数比运算,得到一幅对数比值图RI,对RI进行二维平稳小波4层分解,对每一分解层都采用SSNF的方法进行多尺度积去噪;其次将每一分解层多尺度积去噪后的小波系数进行逆二维平稳小波变换,得到对数比值图RI不同尺度图像layer1~layer4;然后对图像layer1~layer4采用PCA变换方法,并将PCA变换后的第一主成份图作为新的差异图;最后用模糊c均值聚类的方法对差异图NRI进行分割,从而得到变化检测结果。
2.2.1图像预处理
预处理包括图像配准和辐射校正,主要目的是尽可能的减少图像错位以及辐射度变化带来的影响。
图像配准将不同时相获得的图像进行几何对应,用于消除非地物变化的几何畸变。
辐射校正用于消除辐射失真,包括绝对辐射校正和相对辐射校正,其中绝对辐射校正需要通过己知的各种传感器和环境参数得到图像中地物的光谱反射值,相对辐射校正将某个时相遥感图像作为参照图,使相同地物在不同时相图像中具有一致的光谱反射值,该方法是一种克服辐射失真对变化检测影响的有效手段。
2.2.2对数比值法
普遍认为图像比值法比图像差值法对SAR图像中的乘性斑噪具有更好的抵抗能力。
对于SAR图像,一般采用对数比值法,可以先将SAR图像中乘性斑噪转化成加性噪声,便于对SAR图像进行有效地去噪。
2.2.3多尺度积
通过对数变换,将乘性噪声转化为加性噪声,易于多尺度积去噪。
根据小波分解的不同尺度间相关性,采用空间选择性的噪声滤波(SpatiallySelectiveNoiseFilter,SSNF)方法,该方法根据在相邻尺度下同一位置小波系数的相关性,在大的尺度下小波变换系数全部被保留,在较小尺度下只能保留高相关的系数,舍弃低相关系数,从而实现在有噪声信号中保留重要边缘和去噪的目的。
SSNF算法是一个迭代方法,令W(i,j)为图像小波变换系数,具体步骤如下:
(1)定义尺度间小波系数的相关性系数为相邻尺度的系数直接相乘。
其中m为当前层,d为水平、垂直、对角方向,L为小波分解层数
(2)计算相关性系数和小波系数的能量,即
(3)归一化相关性系数
(4)比较
和
图1尺度二维平稳小波变换后的图像
图2尺度积去噪后的图像
图3小波逆变换后的图像
2.2.4主成分分析(PCA)
对上述得到的对数比值图的每个尺度去噪之后的图像进行PCA变换以消除尺度间的相关性。
由于输入图像高度相关的结构成分在PCA变换后的第一主成份图中的占主导成分,而且变化信息将在第一主成份图中增强,因此本文用PCA变换后的第一主成份作为新的差异图,其具体步骤如下:
(1)将上述得到的对数比值图每个尺度去噪之后的图像layer列向量化,并组成矩阵A=(A1,A2,…,AI)
(2)计算矩阵A的相关系数矩阵R
其中矩阵元素
(3)计算相关系数矩阵尺的特征值和特征向量。
(4)图像序列经过PCA变换后得到新的图像Z=Au。
(5)将PCA变换后的第一主成份图ZI作为新的差异图NRI并输出。
图3小波逆变换后的图像
2.2.5模糊c均值聚类(FCM)
对上述得到的新差异图NRI用模糊C均值进行聚类,分为变化类和非变化类。
(1)随机初始化聚类中心ci,i=1,2。
(2)确定隶属度矩阵U
(3)计算代价函数,如果它小于某个确定的阈值,或其相对于上次价值函数的该变量小于某个阈值,则算法停止。
否则,转(4)
(4)修正聚类中心,返回
(2)。
三、实验结果分析
(1)墨西哥郊外
正确率
漏检数
虚警数
总错误数
97.5460%
5737
696
6433
(2)Ottawa地区水灾的RADARSATSAR图像
正确率
漏检数
虚警数
总错误数
89.9689%
6574
6401
173
四、实验代码
%---------------------------------主程序--------------------------%
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);
[mn]=size(I1);
fori=1:
m
forj=1:
n
if(I1(i,j)==0)
I1(i,j)=1;
end
if(I2(i,j)==0)
I2(i,j)=1;
end
end
end
%--------比值图的构造---------%
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);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)=HH(:
:
4);
%imageshow(layer1,layer2,layer3,layer4,2);
%--------尺度积去噪---------%
[tlayer1tlayer2tlayer3tlayer4]=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),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(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_numerror_numundetected_nummisdetected_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);
imshow(Standard_Image);title('±ê×¼²îÒìͼ');
subplot(2,2,4);
imshow(DI);title('±¾Ëã·¨µÄ²îÒìͼ');
disp('ÕýÈ·ÂÊΪ£º');
disp(effective_num/(m*n)*100);
disp('ÓÐЧ¸öÊýΪ£º');
disp(effective_num);
disp('×Ü´íÎóÊýΪ£º');
disp(error_num);
disp('©¼ìÊýΪ£º');
disp(undetected_num);
disp('Ð龯ÊýΪ£º');
disp(misdetected_num);
%-----------------------------------------------------------------%
%--------------------------比值图的构造-----------------------------%
functionRI=CreatRatioImage(I1,I2)
TI1=im2double(I1);
TI2=im2double(I2);
[mn]=size(I1);
RI=zeros(m,n);
fori=1:
m
forj=1:
n
RI(i,j)=log(TI1(i,j))-log(TI2(i,j));
end
end
end
%-----------------------------------------------------------------%
%----------------------------尺度积去噪声----------------------------%
function[tlayer1tlayer2tlayer3tlayer4]=ssnf(layer1,layer2,layer3,layer4)
%---¶¨Òå³ß¶È¼äС²¨ÏµÊýµÄÏà¹ØÐÔϵÊýΪÏàÁڳ߶ȵÄϵÊýÖ±½ÓÏà³Ë--5
corr1=layer1.*layer2;
corr2=layer2.*layer3;
corr3=layer3.*layer4;
corr4=layer4.*layer4;
%---¼ÆËãÏà¹ØÐÔϵÊýµÄÄÜÁ¿---%
[mnc]=size(layer1);
pw1=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pw1(1,d)=temp+layer1(i,j,d)*layer1(i,j,d);
end
end
end
[mnc]=size(layer2);
pw2=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pw2(1,d)=temp+layer2(i,j,d)*layer2(i,j,d);
end
end
end
[mnc]=size(layer3);
pw3=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pw3(1,d)=temp+layer3(i,j,d)*layer3(i,j,d);
end
end
end
[mnc]=size(layer4);
pw4=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pw4(1,d)=temp+layer4(i,j,d)*layer4(i,j,d);
end
end
end
%---¼ÆËãС²¨ÏµÊýµÄÄÜÁ¿---%
[mnc]=size(corr1);
pcorr1=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pcorr1(1,d)=temp+corr1(i,j,d)*corr1(i,j,d);
end
end
end
[mnc]=size(corr2);
pcorr2=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pcorr2(1,d)=temp+corr2(i,j,d)*corr2(i,j,d);
end
end
end
[mnc]=size(corr3);
pcorr3=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pcorr3(1,d)=temp+corr3(i,j,d)*corr3(i,j,d);
end
end
end
[mnc]=size(corr4);
pcorr4=zeros(1,c);
ford=1:
c
temp=0;
fori=1:
m
forj=1:
n
pcorr4(1,d)=temp+corr4(i,j,d)*corr4(i,j,d);
end
end
end
%---¹éÒ»»¯Ïà¹ØÐÔϵÊý---%
Ncorr1=corr1;
ford=1:
4
Ncorr1(:
:
d)=corr1(:
:
d)*sqrt((pw1(1,d))/(pcorr1(1,d)));
end
Ncorr2=corr2;
ford=1:
4
Ncorr2(:
:
d)=corr2(:
:
d)*sqrt((pw2(1,d))/(pcorr2(1,d)));
end
Ncorr3=corr3;
ford=1:
4
Ncorr3(:
:
d)=corr3(:
:
d)*sqrt((pw3(1,d))/(pcorr3(1,d)));
end
Ncorr4=corr4;
ford=1:
4
Ncorr4(:
:
d)=corr4(:
:
d)*sqrt((pw4(1,d))/(pcorr4(1,d)));
end
%---±È½Ï---%
[mnc]=size(layer1);
tlayer1=layer1;
ford=1:
c
fori=1:
m
forj=1:
n
if(abs(Ncorr1(i,j,d))<=abs(layer1(i,j,d)))
tlayer1(i,j,d)=0;
end
end
end
end
[mnc]=size(layer2);
tlayer2=layer2;
ford=1:
c
fori=1:
m
forj=1:
n
if(abs(Ncorr2(i,j,d))<=abs(layer2(i,j,d)))
tlayer2(i,j,d)=0;
end
end
end
end
[mnc]=size(layer3);
tlayer3=layer3;
ford=1:
c
fori=1:
m
forj=1:
n
if(abs(Ncorr3(i,j,d))<=abs(layer3(i,j,d)))
tlayer3(i,j,d)=0;
end
end
end
end
[mnc]=size(layer4);
tlayer4=layer4;
ford=1:
c
fori=1:
m
forj=1:
n
if(abs(Ncorr4(i,j,d))<=abs(layer4(i,j,d)))
tlayer4(i,j,d)=0;
end
end
end
end
end
%-----------------------------------------------------------------%
%---------------------------------PCA-----------------------------%
functionNRI=MPCA(layer1,layer2,layer3)
[rowcol]=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];
%¼ÆËãÏà¹ØϵÊý¾ØÕóR
mm=ones(row*col,1);
Ameans=A;
fork=1:
3
Ameans(:
k)=mean(A(:
k))*mm;
end
R=zeros(3,3);
fori=1:
3
forj=1:
3
sum1=0;sum2=0;sum3=0;
fork=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));
end
end
[VD]=eig(R);%¼ÆËã¾ØÕóRµÄÌØÕ÷ÏòÁ¿¾ØÕóVºÍÌØÕ÷ÖµD
Vsort=fliplr(V);
vsort=Vsort(:
1);
NRI=A*vsort;
NRI=reshape(NRI,row,col);
end
%-----------------------------------------------------------------%
%----------------------------阈值分割-------------------------------%
functionDI=CreatDifferenceImage(NRI)
[rowcol]=size(NRI);
I=reshape(NRI,row*col,1);
t_new=mean(I);
label=ones(row*col,1);
t_old=0;
while(t_new~=t_old)
t_old=t_new;
sum1=0;sum2=0;
n1=0;n2=0;
fori=1:
row*col
if(I(i,:
) label(i,: )=1; sum1=sum1+I(i,: ); n1=n1+1; else label(i,: )=2; sum2=sum2+I(i,: ); n2=n2+1; end end t1=sum1/n1; t2=sum2/n2; t_new=(t1+t2)/2; end fori=1: row*col if(label(i,: )==1) I(i,: )=0; else I(i,: )=255; end end DI=reshape(I,row,col); end %-----------------------------------------------------------------% %----------------------------阈值分割-------------------------------% function[effective_numerror_numundetected_nummisdetected_num]=CDperformence(Standard_Image,DI) [mn]=size(Standard_Image); undetected_num=0; misdetected_num=0; fori=1: m forj=1: n if(Standard_Image(i,j)==255&&DI(i,j)==0) undetected_num=undetected_num+1; elseif(Standard_Image(i,j)==0&&DI(i,j)==255) misdetected_num=misdete
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- SAR 图像 变换 检测