遥感数字图像处理实验报告

实验一 遥感图像统计特性

一、实验目的

掌握遥感图像常用的统计特性的意义和作用，能运用高级程序设计语言实现遥感图像统 计参数的计算。

二、实验内容

编程实现对遥感图像进行统计特性分析，均值、方差（均方差）、直方图、相关系数等。

三、实验原理

1.均值

像素值的算术平均值，反映图像中地物的平均反射强度。

11

00

(,)

N M j i f i j f MN

--===

∑∑

2.方差（或标准差）

像素值与平均值差异的平方和，反映了像素值的离散程度。也是衡量图像信息量大小的 重要参数。

11

2

00

2[(,)]

N M j i f i j f MN

σ--==-=

∑∑

3. 相关系数

反映了两个波段图像所包含信息的重叠程度。f , g 分别为两个波段的图像，它们之间的 相关系数计算公式为：

11

[((,))((,))]

(,)M N f

g f i j e

g i j e C f g ---?-=

∑∑

其中， e f , e g 分别为两个波段图像的均值。 四、实验步骤和内容 1.实验代码

clc clear all

I =imread ('m1.jpg'); whos I %显示图像信息 figure (1),imshow (I ); R =double (I (:,:,1)); G =double (I (:,:,2)); B =double (I (:,:,3));

%求图像的R,G,B 的均值,avg=mean(mean(I))

%求图像的R,G,B的均值

mean(R(:))

mean(G(:))

mean(B(:))

%求R,G,B的方差

varR=var(R(:));

varG=var(G(:))

varB=var(B(:))

%求RG,RB,GB的相关系数

corrcoef(R(:),G(:))

corrcoef(R(:),B(:))

corrcoef(B(:),G(:))

2.原始图像

Figure 1原始图像3.实验结果

R,G,B的均值

R =71.4789

G =99.6597

B =59.5149

R,G,B的方差

varR =1.0614e+03

varG =1.6791e+03

varB =192.6365

RG,RB,GB的相关系数

ans =1.0000 0.9536

0.9536 1.0000

ans =1.0000 0.3460

0.3460 1.0000

ans =1.0000 0.1401

0.1401 1.0000

4. 遥感图像的均值、方差、相关系数、信息熵等统计参数，以及图像直方图的意

义和作用。

①意义：

均值：像素值的算术平均值，反映的是图像中地物的平均反射强度，大小由图像中主体地物的光谱信息决定。

方差：像素值与平均值的平方和，表示像素值的离散程度。

相关系数：反映了两个波段图像所包含信息的重叠程度。

直方图：灰度级的函数，它表示图像中具有某种灰度级的像素的个数，反映了图像中某种灰度出现的频率。

②作用：

测量上的误差以及各种干扰因素的存在，图像的像素值变化具有随机性的特点。多数人认为，遥感图像中某一灰度级内像素出现的频率是服从高斯分布，即密度函数是正态的。一般来说，图像的概率分布难以用某一分析式来表示，但通过分析直方图，灰度级内的像素频数总是可以找出来的。从统计学角度来说，图像的数字特征可作为区分或识别图像中地物的依据。

把图像作为一个随机向量X，按照概率论可以有两种方法来表示，一种用密度函数来表示；另一种用统计参数来表示，这往往是密度函数不可知条件下的表示，如期望、方差、协方差等。使用统计特征可以用来对不同的图像或图像的处理效果进行比较。统计的图像范围根据需要确定，或者是整景图像，或者是指定的地图类型。

直方图：根据直方图的形态可以大致推断图像的反差，然后可通过有目的地改变直方图形态来改善图像的对比度。一般来说，如果图像的直方图形态接近正态分布，则这样的图像反差适中；如果直方图峰值位置偏向灰度值打的一遍，图像偏亮；如果峰值位置偏向灰度值小的一边，图像偏暗；峰值变化过陡、过窄，则说明图像的灰度值国语集中，反差小。

实验二遥感图像增强处理

一、实验目的掌握常用遥感图像的增强方法，能运用高级程序设计语言实现遥感图像的增强处理。二、实验内容

编程实现对遥感图像的IHS 变换、IHS 逆变换、进行统计特性分析，均值、方差（均方差）、直方图、相关系数等。

三、实验原理

1. IHS 变换与IHS 逆变换

IHS 变换：

1

2

I R

V G

V B

????

????

=

????

????

????

?

??

其中，

2

1

arctan

V

H

V

??

= ?

??

S=

IHS 逆变换：

1

2

R I

G V

B V

?

?

?

????

????

=

????

????

????

?

?

?

2. 遥感图像的密度分割

密度分割：将图像的密度或亮度值分成若干等级的处理方法。对单波段遥感图像按灰度分级，对每级赋予不同的色彩，使之变为一幅彩色图像。

3. SPOT 图像模拟真彩色

模拟真彩色：通过某种形式的运算得到模拟的红、绿、蓝三个通道，然后通过彩色合成近似的产生真彩色图像。

（1）SPOT IMAGE 公司提供的方法该方法实际上是将原来的绿波段当作蓝波段，红波段（0.61-0.68 μm）仍采用原来的波

段，绿波段用绿波段、红波段、红外波段的算术平均值来代替。

（2）ERDAS IMAGING 软件中的方法 此法将原来的绿波段当作蓝波段，红波段仍采用原来的波段，绿波段用绿波段、红外波

段按 3：1 的加权算术平均值来代替。 （3）不确定参数法

引入了全色波段（用 P 表示），红色用(aP ＋(1－a)XS3)来表示，绿色用 2P ×XS2/(XS1+ XS2)表示，蓝色用 2P ×XS1/(XS1+XS2)表示。系数 a 根据遥感影像景观取值介于 0.1 到 0.5 之间。

4. 直方图规定化（直方图匹配）

直方图规定化是为了使一个波段图像的直方图变成规定形状的直方图而对图像进行转换的增强方法。

图 2-1 直方图规定化示意图

令 P(r) 为原始图像的灰度密度函数，P(z)是期望通过匹配的图像灰度密度函数。对 P(r)及 P(z) 作直方图均衡变换，通过直方图均衡为桥梁，实现 P(r) 与 P(z) 变换。 步骤： （1）由

()()r

s T r p r dr

==? 01

r ≤≤ 各点灰度由 r 映射成 s 。 （2）由

()()z

v G z p z dz

==?

01

z ≤≤

各点灰度由 z 映射成 v 。 （3）根据 v=G(z), z=G-1(v)

由于 v, s 有相同的分布,逐一取 v=s ，求出与 r 对应的 z=G-1(s)。 离散灰度级情况：

由（1）、（2）计算得两张表，从中选取一对 vk , sj ，使 vk ≈sj ，并从两张表中查得对 应的 rj ，zk 。于是，原始图像中灰度级为 rj 的所有象素均映射成灰度级 zk 。最终得到

所期望的图像。

四、实验代码及图像

1.HIS变换和逆变换

①实验代码

clc

clear all

T=imread('F:\桌面\实验\遥感数字图像处理\shiyan\L7.tif');

T=im2double(T);

R=T(:,:,1);

G=T(:,:,2);

B=T(:,:,3);

[row,col,q]=size(T);

x2=sqrt(2);

x3=sqrt(3);

x6=sqrt(6);

I=zeros(row,col);

V1=zeros(row,col);

V2=zeros(row,col);

H=zeros(row,col);

S=zeros(row,col);

%RGB到ISH ISH变换

for i=1:row

for j=1:col

x=(1/x3)*R(i,j)+(1/x3)*G(i,j)+(1/x3)*B(i,j);

I(i,j)=x;

y=(1/x6)*R(i,j)+(1/x6)*G(i,j)+(-2/x6)*B(i,j);

V1(i,j)=y;

z=(1/x2)*R(i,j)+(-1/x2)*G(i,j)+0*B(i,j);

V2(i,j)=z;

H(i,j)=atan(V1(i,j)/V2(i,j));

S(i,j)=sqrt(V1(i,j)^2+V2(i,j)^2);

end

end

ish=cat(3,I,H,S);

figure,subplot(3,2,1),imshow(T),title('原始图像');

subplot(3,2,2),imshow(ish),title('ISH图像');

subplot(3,2,3),imshow(I),title('亮度I');

subplot(3,2,4),imshow(S),title('饱和度S');

subplot(3,2,5),imshow(H),title('色调H');

%ISH到RGB ISH逆变换

R1=zeros(row,col);

G1=zeros(row,col);

B1=zeros(row,col);

for s=1:row

for t=1:col

x1=(1/x3)*I(s,t)+(1/x6)*V1(s,t)+(1/x2)*V2(s,t);

y1=(1/x3)*I(s,t)+(1/x6)*V1(s,t)+(-1/x2)*V2(s,t);

z1=(1/x3)*I(s,t)+(-2/x6)*V1(s,t)+0*V2(s,t);

R1(s,t)=x1;

G1(s,t)=y1;

B1(s,t)=z1;

end

end

RGB=cat(3,R1,G1,B1);

figure,subplot(3,2,1),imshow(ish),title('ISH图像');

subplot(3,2,2),imshow(RGB),title('RGB图像');

subplot(3,2,3),imshow(R),title('红色波段R');

subplot(3,2,4),imshow(G),title('绿色波段G');

subplot(3,2,5),imshow(B),title('蓝色波段B');

②实验结果

图1 原始图像

图2 HIS变化后的图像

图3-1 正变换H分量图3-2 正变换S分量

图3 正变换I分量

图4 HIS逆变换及其R,G,B分量

2.SPOT图像模拟真彩色

①实验代码

clc

clear all

T=imread('F:\桌面\实验\遥感数字图像处理\shiyan\test23.tif');

T=im2double(T);

G=T(:,:,1);%原图的绿波段

R=T(:,:,2);%原图的红波段

NR=T(:,:,3);%原图的近红外波段

[row,col,q]=size(T);

B1=G;

R1=R;

%SPOT IMAGE公司提供的拟彩色方法

G1=zeros(row,col);

for i=1:row

for j=1:col

G1(i,j)=(G(i,j)+R(i,j)+NR(i,j))/3;

end

end

ncs1=cat(3,R1,G1,B1);

figure,subplot(1,2,1),imshow(T),title('原始图像'); subplot(1,2,2),imshow(ncs1),title('拟彩色图像');

%ERDAS IMAGING软件中提供的拟彩色方法

G2=zeros(row,col);

for i=1:row

for j=1:col

G2(i,j)=(3*G(i,j)+NR(i,j))/4;

end

end

ncs2=cat(3,R1,G2,B1);

figure,subplot(1,2,1),imshow(T),title('原始图像'); subplot(1,2,2),imshow(ncs2),title('拟彩色图像');

五、分析与讨论

实验三遥感图像融合

一、实验目的

掌握多源遥感图像融合的原理与方法，能运用高级程序设计语言实现遥感图像的融合。

二、实验要求

选择IHS 变换、PCA 变换和Brovey 变换三种方法中的一种，编程实现多源遥感图像融合，即将低空间分辨率的多光谱图像与高空间分辨率的全色图像实现融合。原始图像多光谱图像为ETM 或TM 图像，全色图像为SPOT 图像。

三、实验原理

图像融合就是将不同来源的同一对象的图像数据进行空间配准，然后采用一定的算法将各图像数据中所含的信息优势或互补性有机地结合起来产生新图像数据的技术。这种新数据具有描述所研究对象的较多优化的信息表征，同单一信息源相比，能减少或抑制对被感知对象或环境解释中可能存在的多义性、不完整性、不确定性和误差，最大限度地利用各种信息源提供的信息。

常用的遥感图像融合方法有：IHS 变换融合法、主成分（PCA）融合法、Brovey 变换融合法、小波变换融合法等。本实验只给出IHS 变换融合法、主成分（PCA）融合法、Brovey 变换融合法三种方法。

1. 多源遥感图像IHS 变换融合法

RGB 空间和IHS 空间相互转化的处理过程称之为IHS 变换和反变换。IHS 变换可用于两个方面：第一，用于颜色增强；第二，用于多源配准图像合成。关于IHS 用于图像增强的内容本文将不再详述。

IHS 变换是多源遥感图像数据融合常用的方法之一，可以实现不同空间分辨率的遥感图像之间的几何信息的叠加。它首先将RGB 颜色空间的3 个波段的多光谱图像转化到IHS 空间的3 个分量，即色调（H）、饱和度（S）、亮度（I）。然后，将高空间分辨率图像进行对比度拉伸，使它和亮度分量（I）有相同的均值和方差；最后，用拉伸后的高空间分辨

率图像代替亮度分量（I），把它同色调（H）、饱和度（S）进行IHS 逆变换得到融合图像。

IHS 变换：

1

2

I R

V G

V B

????

????

=

????

????

????

?

??

其中，2

1

arctan

V

H

V

??

= ?

??

S=

IHS 逆变换：

1

2

R I

G V

B V

?

?

?

????

????

=

????

????

????

?

?

?

2. 多源遥感图像的PCA 变换融合法

这一基于PCA 的遥感图像融合算法是以3.1 中所给出的PCA 算法为基础的，下面将详细介绍。

（1）PCA 分析

给定M 个波段的相关图像，主成分分析可以用一个组合矩阵S 进行，它是一个M ×M 维的对称矩阵，其中，S(i,i)是第i波段的方差，S(i, j)是波段i与波段j之间的协方差。方差可以通过下式计算

22

11111

111

()

1(1)

N N N N N

i i i i i

i i i i i

VAR x x x N x x x

N N N N

=====

??

????

=-=-

??

???

--

????

??

∑∑∑∑∑

其中，N 为像素的总数，x代表一个波段的图像，y 代表另一个波段的图像。

设一幅多光谱图像由7 个波段组成，需要计算的方差、协方差矩阵S ，其大小就是7×7。这个矩阵的对角元素，由7 个波段的方差组成。另外，有关于主对角线上下对称的21 个元素，代表了由7 个波段得到的任意两个波段之间的协方差。

（2）图像变换

一旦上述变换中矩阵S 的特征向量求得以后，这些特征向量就可用以实现多种图像变换。令X 表示一个矩阵，其各列为求得的这M 个特征向量，因此它是一个M ×M 的矩阵。又令矩阵A表示多光谱或高光谱图像的集合。在这种情况下，A为m×n矩阵，其中N 为每幅图像的像素总数，M 为图像A 的波段数。于是有矩阵变换

'

A TA

='

T

A T A

=

其中，矩阵A′的第一行就是第一主成分图像，第二行就是第二主成分图像，等等。

设M 个特征向量按特征值递减顺序排列，分别为1,,M u u u …, ，

[]1,,T

M T u u u =…, 又令[][]'

1212,,,,,,T

T

M M A Y Y Y A A A A ==…，…

（3）主成分替换与图像重建

将全色（高分辨率）图像与第一主成分图像1 Y 进行直方图匹配后，代替第一主成分然

后进行主成分反变换，即可得到融合后的图像。

（4）算法流程

基于PCA 的遥感图像融合算法的流程如图3-1 所示。

3. Brovey 图像融合法

Brovey 图像融合也称为色彩标准化（Color Normalized ）变换融合，由美国学者Brovey 建立模型并推广而得名。其思想是将多光谱图像的像方空间分解为色彩和亮度成分并进行计算。它是一种常用于多光谱图像增强的比值变换融合方法，比值变换能消除空间或时间变换后产生的增益和偏置因子。其计算公式为：

MS F

P MS MS MS

MS F

P MS MS MS MS

F P MS MS MS R R I R

G B G G I R G B B B I R G B ?=??++??=??++??=

??++?

式中， P I 为源高分辨率全色图像； R MS 、G MS 、B MS 分别为源低分辨率多光谱图像的红、绿、蓝通道图像； R F 、G F 、B F 分别为融合图像的红、绿、蓝通道图像。得到融

合图的各个波段数据后，还必须对其进行灰度拉伸，使其灰度值范围与源多光谱图像各波段的灰度值范围一致，最终得到彩色的融合图像。 四、实验内容 1.实验代码

clear all;

%%%读取高分辨率图像

f1=imread('pan00.jpg ');

subplot(2,2,1),imshow(f1)

%%%利用插值将多光谱图像放大到与高分辨率图像一样大小[M,N]=size(f1);

f2=imread('mul00.jpg ');

f2=imresize(f2,[M,N],'bilinear');

subplot(2,2,2);imshow(f2);

%%%将RGB空间转换为HIS

f1=double(f1);

f2_hsi=rgb2hsv(f2);

f2_h=f2(:,:,1);

f2_s=f2(:,:,2);

f2_i =f2(:,:,3);

%%%进行小波分解

[c1 s1]=wavedec2(f1,1,'sym4');

f1=im2double(f1);

[c2_h s2_h]=wavedec2(f2_h,1,'sym4');

[c2_s s2_s]=wavedec2(f2_s,1,'sym4');

[c2_i s2_i]=wavedec2(f2_i,1,'sym4');

%%%对系数进行融合

c_h=0.5*(c2_h+c1);

c_s=0.5*(c2_s+c1);

c_i =c1;

%%%分别对H分量、I分量、S分量进行直方图均衡化

f_h= waverec2(c_h,s1,'sym4');

f_h=histeq(f_h);

f_s= waverec2(c_s,s1,'sym4');

f_s=histeq(f_s);

f_i = waverec2(c_i,s1,'sym4');

f_i=histeq(f_i);

%%%显示融合后图像

g=cat(3,f_h,f_s,f_i);

subplot(2,2,3);imshow(g);

2.实验结果

图1 多光谱图像

图2 高分辨率全射图像

图3 融合后的图像

五、分析与思考

1.图像融合与色彩合成之间的区别

实验四图像变换及信息提取

一、实验目的

1.掌握遥感图像变换方法傅里叶变换和代数运算；

2.会应用这些变换提取遥感图像的相关信息，如植被指数等；

3.能运用高级程序设计语言实现遥感图像的变换处理。

二、实验要求

1.编程实现图像的Fourier变换及其反变换，输出图像的频谱图像；

2.编程计算遥感图像的植被指数。

三、实验原理

1．图像的傅里叶变换

(1)应用傅里叶变换进进行图像处理理

傅里叶变换换是线性系统统分析的一个有力工具，它它能够定量地地分析诸如数数

字化系统、采样点、电子放大器器、卷积滤波器器、噪音和显显示点等的作作用。通

过实实验培养这项项技能，将有助于解决决大多数图像像处理问题。对对任何想在工工作中有效应应用数字图像像处理技术的的人来说，把把时间用在在学习和掌握握博里叶变换换上是很有必要的。 (2)傅里叶（Fourier ）变换的定义

对于二维信号，二维Fourier 变换定义为：

2()2()

(,)(,)(,)(,)cos sin j ux vy j ux vy j F u v f x y e dxdy

f x y F u v e dudv e j ππθθθ∞∞-+-∞-∞∞∞+-∞-∞?=?

?=??

=+?

?

????

二维离散傅里叶变换为：

11

2(//)001(,)(,)(0,1,2......1;0,1,2.......1)

N M j ux N vy M x y F u v f x y e MN u N v M π---+====-=-∑∑

11

2(//)

00(,)(,)(0,1,2.......1;0,1,2.......1)

N M j ux N vy M u v f x y F u v e x N y M π--+====-=-∑∑

二维离散傅里叶变换的频谱和相位角：

(,)(,)

(,)arctan

(,)

F u v I u v u v R u v ?==

（3）图像的快速傅立叶变换

图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法。实际上，现在有实现傅立叶变换的芯片，可以实时实现傅立叶变换。

快速傅立叶变换（FFT ）的计算公式:

11

2(//)

00

(,)(,)N M j ux N vy M u v F u v F u v e π--+==←∑∑ 快速傅立叶逆变换（IFFT ）的计算公式:

11

2(//)

00

(,)(,)N M j ux N vy M u v f u v f u v e π--+==←∑∑

2.遥感图像的植被指数 （1）比值植被指数RVI ：IR RVI R

=

① 绿色健康植被覆盖地区的 RVI 远大于1，而无植被覆盖的地面（裸土、人工建筑、水体、植被枯死或严重虫害）的RVI 在1 附近。植被的RVI 通常大于2；

②RVI 是绿色植物的灵敏指示参数，与LAI 、叶干生物量(DM)、叶绿素含量相关性高，

可用于检测和估算植物生物量；

③植被覆盖度影响RVI，当植被覆盖度较高时，RVI 对植被十分敏感；当植被覆盖度<50%时，这种敏感性显著降低；

④RVI 受大气条件影响，大气效应大大降低对植被检测的灵敏度，所以在计算前需要进行大气校正，或用反射率计算RVI。

（2）归一化植被指数NDVI：

IR R NDVI

IR R

-

=

+

①NDVI 的应用：检测植被生长状态、植被覆盖度和消除部分辐射误差等；

②-1≤NDVI≤1，负值表示地面覆盖为云、水、雪等，对可见光高反射；0 表示有岩石或裸土等，NIR 和R 近似相等；正值，表示有植被覆盖，且随覆盖度增大而增大；

③NDVI 的局限性表现在，用非线性拉伸的方式增强了NIR 和R 的反射率的对比度。对于同一幅图象，分别求RVI 和NDVI 时会发现，RVI 值增加的速度高于NDVI 增加速度，即NDVI 对高植被区具有较低的灵敏度；

④NDVI 能反映出植物冠层的背景影响，如土壤、潮湿地面、学、枯叶、粗超度等，且与植被覆盖有关；

⑤NDVI 对绿色植物的生长状态和空间分布密度反映敏感，但受土壤背景的影响较大，适用于作物生长早期或植被覆盖度较低的区域。

（3）正交植被指数PVI

PVI =1.6225(IR) ?2.2978(R) +11.0656

PVI = 0.939(IR) ?0.344(R) + 0.09

其中，IR 为遥感多波段图像中的近红外（Infrared）波段；R 为红波段。前式适用于NOAA 卫星的AVHRR 数据，后式适用于NOAA 卫星的AVHRR 数据。

四、实验内容

1.Fourier变换及其反变换

①实验代码

clear all

clc

RGB=imread('image1.jpg');

figure(1),imshow(RGB);

I=rgb2gray(RGB);

figure(2),imshow(I);

B=fftshift(fft2(I));% 将频谱的坐标原点由左上角移至屏幕中央

A1=abs(B);

R=real(B);I=imag(B);% 取Fourier变换结果的实部R 和虚部I A2=sqrt(R.^2+I.^2);% 频谱

P1=log(abs(B))+1;% 方法1：使用对数变换，获取更好的效果图

figure(3),imshow(P1,[])

m=min(min(min(A2)));

n=max(max(max(A2)));

P2=round((A2-m)/(n-m)*255);%方法2：归一化,获取更好的效果图

figure(4),imshow(P2)