遥感数字图像处理实习报告含Matlab处理代码word版

辽宁XX大学数字图像处理上机实习报告教学单位辽宁XX大学专业摄影测量与遥感实习名称遥感数字图像处理班级测绘研11-3班学生姓名路聚峰学号47112指导教师孙XX实习1读取BIP、BIL、BSQ文件

一、实验目的用Matlab读取BIP、BIL、BSQ文件，并将结果显示出来。遥感图像包括多个波段，有多种存储格式，但基本的通用格式有3种，即BSQ、BIL和BIP格式。通过这三种格式，遥感图像处理系统可以对不同传感器获取的图像数据进行转换。BSQ是像素按波段顺序依次排列的数据格式。BIL格式中，像素先以行为单位块，在每个块内，按照波段顺序排列像素。BIP格式中，以像素为核心，像素的各个波段数据保存在一起，打破了像素空间位置的连续性。用Matlab读取各个格式的遥感数据，是图像处理的前提条件，只有将图像读入Matlab工作空间，才能进行后续的图像处理工作。

二、算法描述

1.调用fopen函数用指定的方式打开文件。

2.在for循环中调用fread函数，用指定的格式读取各个像素。

3.用reshape函数，重置图像的行数列数。

4.用imadjust函数调整像素的范围，使其有一定对比度。

5.用imshow显示读取的图像。

三、Matlab源代码

1.读取BSQ的源代码：clearallclclines=400;samples=640;N=6;img=fopen(D:sample_SQ,rb);fori=1:Nbi=fread(img,lines_samples,uint8);band_ov=reshape(bi,samples,lines);band_ov2=band_ov;band_int8=uint8(band_ov2);tif=imadjust(band_int8);mkdir(D:MATLAB,tifbands1)name=D:MATLABtifbands1tif,int2str(i),.tif;imwrite(tif,name,tif);tilt=波段,int2str(i);subplot(3,2,i),imshow(tif);title(tilt);endfclose(img);

2.读取BIP源代码clearallclclines=400;samples=640;N=6;fori=1:Nimg=fopen(D:MATLABsample_IP,rb);b0=fread(img,i-

1，uint8);b=fread(img,lines_samples,uint8,(N-1);band_ov=reshape(b,samples,lines);band_ov2=band_ov;?nd_int8=uint8(band_ov2);tif=imadjust(band_int8);mkdir(E:MATLAB,tifbands)name=E:MATLABtifbandstif,int2str(i),.tif;imwrite(tif,name,tif);%imwrite(A,filename,fmt)tilt=波段,int2str(i);subplot(3,2,i),imshow(tif);title(tilt);fclose(img);end

3.读取BIL的源代码clearallclclines=400;samples=640;N=6;fori=1:Nbi=zeros(lines,samples);forj=1:samplesimg=fopen(D:MATLABsample_IL,rb);bb=fread(img,(i-1)_640,uint8);b0=fread(img,1_(j-1),uint8);bandi_inej=fread(img,lines,uint8,1_(N_samples-1);fclose(img);bi(:,j)=bandi_inej;endband_int8=uint8(bi);tif=imadjust(band_int8);mkdir(D:MATLAB,tifbands)name=D:MATLABtifbandstif,int2str(i),.tif;imwrite(tif,name,tif);tilt=,int2str(i);subplot(3,2,i),imshow(tif);title(tilt);end

四、运行结果图1：读取文件的六个波段图实习2均值滤波、边缘信息提取

一、实验目的与原理各种图像滤波算子可以实现图像的增强，去噪，边缘提取等。图像增强的目的在于：1.采用一系列技术改善图像的视觉效果，提高图像的清晰度，

2.将图像转换成一种更适合人或机器进行分析处理的形式。它不是以图像保真度为原则，而是通过处理，设法有选择地突出便于人或机器分析某些感兴趣的信息，抑制一些无用的信息，以提高图像的使用价值。图像增强方法从增强的作用域出发，可分为空间域增强和频率域增强。空间域增强就是直接对图像像素灰度进行操作；频率域增强是对图像经傅里叶变换后的频谱成分进行操作，然后经傅里叶逆变换获得所需结果。图像滤波可分为空间域滤波和频率域滤波，前者通过窗口或卷积核进行，它参照相邻像素改变单个像素的灰度值。后者对图像进行傅XX变换，然后对频谱进行滤波。空间域图像滤波称为平滑和锐化，强调像素与其周围相邻像素的关系。去噪滤波为平滑滤波包括均值滤波和中值滤波。锐化滤波包括罗伯特梯度、索伯尔梯度、拉XX算法、定向检测，用以提取线状地物和边缘。此实验用Matlab采用各种滤波对图像进行了处理，处理结果如下：二、算法描述

1.用imread读取图像文件，并用size获取图像的大小。

2.设计各种滤波算子。

3利用卷积公式对图像的没一个像素进行处理，得到滤波后的图像。

4.用imshow显示滤波后的图像。

三、Matlab源代码

1.均值滤波源码：clearallclcimg=imread(

XXX);row,column,band=size(img);img0=double(img);f11=;f12=;f13=;f21=;f22=;f23=;f31=;f32=;f33=;img1=img0(:,

1，:),img0(:,:,:),img0(:,column,:);img2=img1(

1，:,:);img1(:,:,:);img1(row,:,:);filtered=zeros(row,column,band);forii=1:rowforjj=1:columnfiltered(ii,jj,:)=f11_img2(ii,jj,:) f12_img2(ii,jj 1，:) f13_img2(ii,jj 2,:) .f21_img2(ii 1，jj,:) f22_img2(ii 1，jj 1，:) f23_img2(ii 1，jj 2,:) .f31_img2(ii 2,jj,:) f32_img2(ii 2,jj 1，:) f33_img2(ii 2,jj 2,:);endendfiltered1=uint8(filtered);subplot(

1，2,1),imshow(img);title(图1原始RGB图像);subplot(

1，2,2),imshow(filtered1);title(图2均值滤波后的图像);imwrite(filtered

1，XXX);

2.边缘提取滤波源代码clearallimg=imread(

XXX);row,column,band=size(img);img0=double(img);f11=1;f12=0;f13=-1;f21=1;f22=0;f23=-1;f31=1;f32=0;f33=-1;img1=img0(:,

1，:),img0(:,:,:),img0(:,column,:);img2=img1(

1，:,:);img1(:,:,:);img1(row,:,:);filtered=zeros(row,column,band);forii=1:rowforjj=1:columnfiltered(ii,jj,:)=f11_img2(ii,jj,:) f12_img2(ii,jj 1，:) f13_img2(ii,jj 2,:) .f21_img2(ii 1，jj,:) f22_img2(ii 1，jj 1，:) f23_img2(ii 1，jj 2,:) .f31_img2(ii 2,jj,:) f32_img2(ii 2,jj 1，:) f33_img2(ii 2,jj 2,:);endendfiltered1=uint8(filtered);subplot(

1，2,1),imshow(img);title(图1RGB原图像);subplot(

1，2,2),imshow(filtered1);title(图2边缘提取后的图像);imwrite(filtered

1，XXX);

四、运行结果图1：原始RGB图像图2：均值滤波后的图像图3：边缘提取后的图像实习3傅里叶变换、傅里叶逆变换，及频域滤波

一、实验目的按照信号处理理论，根据滤除的频率特征，滤波有3种：1.低通滤波。低通滤波是对频率域的图像通过滤波器H（u，v）削弱或抑制高频部分而保留低频部分的滤波方法。由于图像上的噪声主要集中在高频部分，所以低通滤波可以起到压抑噪声的作用。同时，由于强调了低频成分，图像会变得比较平滑。

2.高通滤波。高通滤波是对频率域的图像通过滤波器来突出图像的边缘和轮廓，进行图像锐化的方法。

3.带通滤波。仅保留指定频率范围的滤波，范围外的频率被阻止。将空间域中的图像变换到频率域中进行计算。空间增强技术强调像元位置和像元之间的关系，但随着考虑的像元数目增多，计算的复杂度增加而且非常耗费计算运算时间，特别是当模板越来越大时，这种现象尤为明显。频率域增强方法：1.频率域平滑：保留图像的低频部分而抑制高频部分。

2.频率域锐化：保留图像的高频部分而削弱低频部分。首先将空间域图像通过傅XX变换为频率域图像，然后选择合适的滤波器对的频谱成分进行增强得到图像，再经过傅XX逆变换将变换到空间域，得到增强后的图像。根据傅里叶变换的原理，用Matlab实现对图像的傅里叶变换，再设计各种频率滤波器，包括理想滤波器、巴特沃斯滤波器、指数滤波器等高通或低通滤波器，用这些滤波器对频率图像进行滤波，将得到的滤波后的频率图像经过傅里叶逆变换后得到想要的图像。本实验，对这些滤波器都进行了设计，处理结果如下图：二、算法描述

1.读取RGB图像，并获得其某个波段。

2.调用fft2对某波段进行傅里叶变换，用fftshift频移函数，使得图像的低频部分移动到频谱的中心。

3.设置低通高通滤波器，对频谱图像进行滤波处理，去除其高频或低频部分。

4.用ifft2对频谱图像进行逆傅里叶变换，得到滤波后的图像。

5对图像进行归一化，使像素值在0-255之间。

XXX显示频谱图像和滤波后的图像。

三、Matlab源代码

1.理想低通滤波源代码clearalla=imread(

XXX);_,Y,Z=size(a);a1=a(:,:,1);b1=fft2(a1);b2=fftshift(b1);F=abs(b2);h=zeros(_,Y);cutoff=0.7;threshold=1-cutoff;low_=round(-threshold_);up_=round( threshold_);lowy=round(-threshold_);upy=round( threshold_);h(low_:up_,lowy:upy)=1;lowpass=b

XXX;d0=ifft2(lowpass);result=abs(d0);result=uint8(result);F1=log10(reshape(F,__Y,1);F2=uint8((F1-min(F1ma_(F1)-min(F1)_255);F3=reshape(F2,_,Y);subplot(

1，3,1),imshow(F3),title(Fig.1傅里叶频谱);subplot(

1，3,2),imshow(h),title(Fig.2理想低通滤波器);subplot(

1，3,3),imshow(result),title(Fig.3理想低通滤波结果);

2.巴特沃斯低通滤波代码：clearalla=imread(

XXX);_,Y,Z=size(a);a1=a(:,:,1);b1=fft2(a1);b2=fftshift(b1);F=abs(b2);h=zeros(_,Y);lowpass=zeros(_,Y);n1=round();n2=round();dma_=sqrt(n12 n22);cutoff=0.7;d0=(1-cutoff)_dma_;n=3;for_=1:_fory=1:Yd=sqrt(_-n1)2 (y-n2)2);h(_,y)=1 (0)(2_n);lowpass(_,y)=b2(_,y)._h(_,y);endendlp=ifft2(lowpass);result=abs(lp);result=uint8(result);F1=log10(reshape(F,__Y,1);F2=uint8((F1-min(F1ma_(F1)-min(F1)_255);F3=reshape(F2,_,Y);subplot(

1，3,1),imshow(F3),title(Fig.1傅里叶频谱);subplot(

1，3,2),imshow(h),title(Fig.2巴特沃斯低通滤波器);subplot(

1，3,3),imshow(result),title(Fig.3巴特沃斯低通滤波结果);

3.指数低通滤波的主要代码：for_=1:_fory=1:Yd=sqrt(_-n1)2 (y-n2)2);h(_,y)=1 e_p(-(0)(-n);lowpass(_,y)=b2(_,y)._h(_,y);endend

4.指数高通滤波的主要代码：for_=1:_fory=1:Yd=sqrt(_-n1)2 (y-n2)2);h(_,y)=1 e_p(-(0)(n);lowpass(_,y)=b2(_,y)._h(_,y);endend

四、运行结果图1：傅里叶频谱图2：理想低通滤波器图3：傅里叶低通滤波结果图4：理想高通滤波器图5：傅里叶高通滤波结果图6：巴特沃斯低通滤波器图7：巴特沃斯低通滤波结果图8：巴特沃斯高通滤波器图9：巴特沃斯高通滤波结果图10：指数高通滤波器图11：指数高通滤波结果图12：指数低通滤波器图13：指数低通滤波结果实习4主成分变换、主成分逆变换

一、实验目的与原理遥感多光谱影像波段多，信息量大，对图像解译很有价值。但数据量太大，在图像处理计算时，也常常耗费大量的机时，占据大量的磁盘空间。实际上，一些波段的遥感数据之间都有不同程度的相关性，存在着数据冗余。多光谱变换方法可通过函数变换，达到保留主要信息、降低数据量、增强或提取有用信息的目的。其变换的本质是对遥感图像进行线性变换，使多光谱空间的坐标系按一定规律进行旋转。多光谱变换主要有两种变换：主成分变换和缨帽变换。

(1)主成分变换（K-L变换）的目的：1.K-L变换是图像分析与模式识别中的重要工具，用于特征抽取，降低特

2.征数据的维数

3.K-L变换用于图像压缩时可以实现有损压缩和无损压缩

4.K-L变换后的N幅图像统计上互不相关，因此K-L变换可以去除图像数据的相关性，提取主要信息。用Matlab实现六个波段图像的主成分分析。求出六个波段的协方差阵，再求协方差阵的特征值、特征向量。只取特征值比较大的数，其占各个波段主要成分，特征值小的，大部分噪声等内容。

二、算法描述主成分变换过程：1读取六个波段的图像。

2.将每个波段的所有像素存成一列，共六列，形成MN_6矩阵。

3.求出每个波段的平均像素值。

4.利用算出的平均像素值，求MN_6矩阵的协方差阵。

5.求协方差阵的特征值与特征向量。

6.将特征值按从大到小排列，特征向量与特征值对应。主成分的逆变换过程：(1)确定保留的波段数

(2)将去除的波段用0值替代

(3)将去噪后的PCA结果SCORE_inv（COEFF）

(4)去中心化（各波段加上其均值）

(5)重新排列图像行列数

三、Matlab源代码

1.主成分分析代码clearallclcnum_and=6;?=imread(tif

XXX);b2=imread(tif

XXX);b3=imread(tif

XXX);b4=imread(tif

XXX);b5=imread(tif

XXX);b6=imread(tif

XXX);_,y=size(b1);img0=reshape(b

1，__y,1),reshape(b2,__y,1),reshape(b3,__y,1),reshape(b4,__y,1),reshape(b5,__y,1),reshape(b6,__y,1);img1=double(img0);average=mean(img1);img3=img1-repmat(average,__y,1);c=cov(img3);coff,lam=eig(c);lam1=lam(6:-1:

1，6:-1:1);coff2=coff(:,6:-1:1);lam_al=diag(lam1);lam_um=sum(lam_al);lam__um=sum(lam_al(1:3);prop=lam__uam_um;result=img3_coff2;result2=reshape(result,_,y,6);result3=uint8(result2(:,:,1)-min(min(result2(:,:,1a_(ma_(result2(:,:,1)-min(min(result2(:,:,1)_255);imshow(result3);holdonresult3=uint8(result2(:,:,2)-min(min(result2(:,:,2a_(ma_(result2(:,:,2)-min(min(result2(:,:,2)_255);figure,imshow(result3);result3=uint8(result2(:,:,3)-min(min(result2(:,:,3a_(ma_(result2(:,:,3)-min(min(result2(:,:,3)_255);figure,imshow(result3);result3=uint8(result2(:,:,4)-min(min(result2(:,:,4a_(ma_(result2(:,:,4)-min(min(result2(:,:,4)_255);figure,imshow(result3);result3=uint8(result2(:,:,5)-min(min(result2(:,:,5a_(ma_(result2(:,:,5)-min(min(result2(:,:,5)_255);figure,imshow(result3);result3=uint8(result2(:,:,6)-min(min(result2(:,:,6a_(ma_(result2(:,:,6)-min(min(result2(:,:,6)_255);figure,imshow(result3);result2(:,:,4:6)=0;result_ew=reshape(result2,__y,6);result_pca=result_ew_inv(coff2);result_pcaz2=result_pca repmat(average,__y,1);result_pcaz3=reshape(result_pcaz2,_,y,6);figure,imshow(uint8(result_pcaz3(:,:,1);C=cat(3,uint8(result_pcaz3(:,:,3),uint8(result_pcaz3(:,:,4),uint8(result_pcaz3(:,:,5);figure,imshow(C);

四、运行结果图1：第一主分量图2:第二主分量图3：第三主分量图4：第四主分量图5：第五主分量图6：第六主分量图7：逆变换后的图像图8：逆变换后的彩色合成图像由图可以看出，图像的主要成分主要集中在前三个主分量图像上，后三个图像大部分都是噪声。通过逆变换恢复原来的图像，这样，经过K-L变换后的前三个主分量，信噪比大，噪声相对小，突出了主要信息，达到了图像增强的目的。另外，经过K-L变换后，第一或前二或前三个主分量已经包含了绝大多数的地物信息，足够分析使用，同时数据量大大减少。应用中，常常只取前三个主分量作假彩色合成（如图8），数据量可减少到43%，既实现了数据压缩，也可作为分类前的特征选择。实习5缨帽变换

一、实验目的缨帽变换（K-T变换）的目的：该变换也是一种坐标空间发生旋转的线性组合变换，但旋转后的坐标轴不是指向主成分方向，而是指向与地物特别是和植被生成以及土壤有密切关系的方向

二、算法描述

1.读取六个波段的图像。

2.将每个波段的所有像素排成一列，共六列组成一个矩阵。

3.设计转换矩阵。

4.图像矩阵与转换矩阵相乘，得到六个相互垂直的分量。

5.逐一显示各个分量。

三、Matlab源代码clearallclcb1=imread(tif

XXX);b2=imread(tif

XXX);b3=imread(tif

XXX);b4=imread(tif

XXX);b5=imread(tif

XXX);b6=imread(tif

XXX);_,y=size(b1);img0=reshape(b

1，__y,1),reshape(b2,__y,1),reshape(b3,__y,1),reshape(b4,__y,1),reshape(b5,__y,1),reshape(b6,__y,1);img0=double(img0);KT_m=0.356

1，0.3972,0.3904,0.6966,0.2286,0.1596;-0.3344,-0.3544,-0.4556,0.6966,-0.0242,-0.2630;0.2626,0.214

1，0.0926,0.0656,-0.7629,-0.5388;0.0805,-0.0498,0.1950,-0.1327,0.5752,-0.7775;-0.7252,-0.0202,0.6683,0.06

1，-0.1494,-0.0274;0.4000,-0.8172,0.3832,0.0602,-0.1095,0.0985;img1=img0_KT_m;min0=min(img1);ma_0=ma_(img1);dm=ma_0-min0;img11=img1-repmat(min0,__y,1);img2=(img

1epmat(dm,__y,1)_255;image=reshape(img2,_,y,6);figure，imshow(uint8(image(:,:,1),title(KT变换第一分量：亮度);figure，imshow(uint8(image(:,:,2),title(KT变换的第二分量：绿度);figure，imshow(uint8(image(:,:,3),title(KT变换的第三分量：湿度);figureimshow(uint8(image(:,:,4),title(KT变换的第四分量：黄度指数及噪声);

四、运行结果图9：第一分量：亮度图图10：第二分量：绿度图11：第三分量：湿度图12：第四分量：黄度指数及噪声用Matlab编写代码实现缨帽变换，其基本思想就是：采用变换矩阵，使图像指向与地物特别是和植被生长以及土壤有密切关系的方向。变换后，这六个分量相互垂直，前三个分量具有明确的地物意义：y1为亮度，实际上是TM六个波段的加权和，反映出图像总体的反射值；y2为绿度，y3为湿度。实习6图像分类（最小距离分类、最大似然分类、K均值分类）

一、实验目的与原理图像分类的目的是将图像中每个像元根据其不同波段的光谱亮度、空间结构特征或者其他信息，按照某种规则或算法划分为不同的类别。遥感图像分类就是对地球表面及其环境在遥感图像上的信息进行识别和分类，从而达到识别图像信息所对应的实际地物，提取所需地物信息的目的。遥感图像分类的原理：由于地物的成分、性质、分布情况的复杂程度和成像条件不同，以及一个像元或瞬时试场里往往有两种或多种地物的情况，即混合像元，使得同类地物的特征向量也不尽相同，而且使得不同地物类型的特征向量之间的差别也不都是显著的。遥感图像的光谱特征通常是以地物在多光谱图像上的灰度体现出来的，即不同地物在同一波段图像上表现的灰度一般互不相同；同时，不同地物在多个波段上的灰度呈现规律也不相同，这就构成了我们在图像上区分不同地物的物理依据。遥感图像传统的计算机分类方式有监督分类和非监督分类两种。监督分类包括最小距离分类法、平行多面体分类法、最大似然分类法。非监督分类包括K-均值聚类法、ISODATA分类法。本实验将用Matlab实现监督分类的最小距离分类和最大似然分类，以及非监督分类中的K-均值分类。

三、算法描述

(1)最小距离分类算法描述

1.读取各个波段文件，并加载各类地物在各个波段的均值文件。

2.计算每个像元向量到各类地物的距离。

3.比较像元到各类地物的距离，取最小的距离，并将像元归为该类。

4.为各类赋予一定的颜色。

(2)最大似然分类算法描述：1.读取BIP文件，将每个波段读为一列。

2.加载各类地物在各个波段的均值和方差

3.对每个像元向量计算属于各个类的概率，比较概率的大小，概率最大的将像元归为该类。

4.对每个类赋予一定的颜色。

5.显示分类后的图像。

(3)K-means分类算法描述：1为每个聚类确定一个初始的聚类中心（一定要在最小值和最大值之间，这一步非常重要）2将样本集中的每一个样本按照最小距离原则,分配到k个聚类中的某一个

3.使用每个聚类中所有样本的均值作为新的聚类中心4如果聚类中心有变化则重复

2、3步直到聚类中心不再变化为止5最后得到的k个聚类中心就是聚类的结果

四、Matlab源代码

(1)最小距离分类代码：clearallclc_=400;y=640;N=6;image=zeros(__y,N);fori=1:Nimg=fopen(sample_IP,rb);b0=fread(img,i-

1，uint8);image(:,i)=fread(img,__y,uint8,N-1);fclose(img);XXX;XXX;XXX;loadsoil____t;loadsoil____t;loadsoil____t;dis=zeros(__y,N);dis(:,1)=sqrt(sum(image-repmat(crop(

1，:),__y,1).2,2);dis(:,2)=sqrt(sum(image-repmat(forest(

1，:),__y,1).2,2);dis(:,3)=sqrt(sum(image-repmat(water(

1，:),__y,1).2,2);dis(:,4)=sqrt(sum(image-repmat(soil1(

1，:),__y,1).2,2);dis(:,5)=sqrt(sum(image-repmat(soil2(

1，:),__y,1).2,2);dis(:,6)=sqrt(sum(image-repmat(soil3(

1，:),__y,1).2,2);min_is=min(dis,2);classes=zeros(__y,3);fori=1:__yind=find(dis(i,:)=min_is(i);switchindcase1classes(i,:)=255,0,0;case2classes(i,:)=0,255,0;case3classes(i,:)=0,0,255;case4classes(i,:)=255,255,0;case5classes(i,:)=255,0,255;case6classes(i,:)=0,255,255;endendimg2=reshape(classes,y,_,3);RGB_lasses=cat(3,img2(:,:,1),img2(:,:,2),img2(:,:,3);imshow(uint8(RGB_lasses);title(最小距离分类);

(3)最大似然分类代码：clearallclc_=400;y=640;N=6;image=zeros(__y,N);fori=1:Nimg=fopen(sample_IP,rb);b0=fread(img,i-

1，uint8);image(:,i)=fread(img,__y,uint8,N-1);fclose(img);XXX;XXX;XXX;loadsoil