对高光谱遥感数据的分析与处理

对高光谱遥感数据的分析与处理
姓名：张俊飞
班级：021051
学号：02105058
E-mail：jeffei@
时间：2013年4月25日
对高光谱遥感数据的分析与处理
一、高光谱成像介绍
高光谱成像技术是近二十年来发展起来的基于非常多窄波段的影像数据技术，其最突出的应用是遥感探测领域，并在越来越多的民用领域有着更大的应用前景。

高光谱成像技术集中了光学、光电子学、电子学、信息处理、计算机科学等领域的先进技术，是传统的二维成像技术和光谱技术有机的结合在一起的一门新兴技术。

近几年年来，自然灾害频发，所以，及时、准确的灾情评估对决策部门制定科学和有效的救灾减灾方案具有关键性的作用。

遥感具有数据获取范围广、速度快等特点，应用在灾害评估中具有非常大的优势和潜力。

在我国近年来的多次重大自然灾害评估中，遥感技术都发挥了极其重要的作用。

遥感技术的应用不止于此。

下面列举了主要的应用方面：
1.气象：天气预报、全球气候演变研究；
2.农业：作物估产、作物长势及病虫害预报；
3.林业：调查森林资源、监测森林火灾和病虫害；
4.水文与海洋：水资源调查、水资源动态研究、冰雪监控、海洋渔业；
5.国土资源：国土资源调查、规划和政府决策；
6.环境监测：水污染、海洋油污染、大气污染、固体垃圾等及其预报；
7.测绘：航空摄影测量测绘地形图、编制各种类型的专题地图和影像地图；
8.地理信息系统：基础数据、更新数据。

虽然拥有诸多优点，但其本身带有很大的数据，对硬件和软件有很高的要求，本文中，先不对硬件进行讨论，就软件方面，对数据进行一系列处理，做到既不丢失其主要数据，又能降低其时空复杂度。

二、PCA理论基础
对测试数据库说明如下：
AVIRIS高光谱数据92AV3C：该场景由AVIRIS传感器于1992年6月获得，该数据为145*145大小，有220个波段。

该数据及真实标记图可以由因特网下载：http://www.ehu.es/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes。

该数据共包含16个类别。

该数据维数为200维，维数较高，我们希望找到一种简洁的算法来把它的维数降下来，这样处理数据的速度可以加快，节省人力物力。

2.1PCA简介
主成分分析(Principal Component Analysis,简称PCA)是一种常用的基于变量
协方差矩阵对信息进行处理、压缩和抽提的有效方法。

算法如下：
Step1分别求各维的平均值，然后对于所有的样例，都减去对应的均值，得到DataAdjust(m*n)
Step2求特征协方差矩阵。

Step3求协方差的特征值和特征向量，得到特征向量与特征值
Step4将特征值按照从大到小的顺序排序，选择其中最大的k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵。

Step5将样本点投影到选取的特征向量上。

假设样例数为m，特征数为n，减去均值后的样本矩阵为DataAdjust(m*n)，协方差矩阵是n*n，选取的k
个特征向量组成的矩阵为EigenVectors(n*k)。

那么投影后的数据FinalData
2.2算法合理性的求证----最大方差理论
要解释为什么协方差矩阵的特征向量就是k维理想特征，我看到的有三个理论：分别是最大方差理论、最小错误理论和坐标轴相关度理论。

这里简单探讨前两种，最后一种在讨论PCA意义时简单概述。

在信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。

如前面的图，样本在横轴上的投影方差较大，在纵轴上的投影方差较小，那么认为纵轴上的投影是由噪声引起的。

因此我们认为，最好的k维特征是将n维样本点转换为k维后，每一维上的样本方差都很大。

比如下图有5个样本点：（已经做过预处理，均值为0，特征方差归一）
下面将样本投影到某一维上，这里用一条过原点的直线表示（前处理的过程实质是将原点移到样本点的中心点）。

假设我们选择两条不同的直线做投影，那么左右两条中哪个好呢？根据我们之前的方差最大化理论，左边的好，因为投影后的样本点之间方差最大。

这里先解释一下投影的概念：
红色点表示样例，蓝色点表示在u上的投影，u是直线的斜率也是直线的方向向量，而且是单位向量。

蓝色点是在u上的投影点，离原点的距离是
（即或者）由于这些样本点（样例）的每一维特征均值都为0，因此投影到u上的样本点（只有一个到原点的距离值）的均值仍然是0。

回到上面左右图中的左图，我们要求的是最佳的u，使得投影后的样本点方差最大。

由于投影后均值为0，因此方差为：
中间那部分很熟悉啊，不就是样本特征的协方差矩阵么（的均值为0，一般协方差矩阵都除以m-1，这里用m）。

用来表示，表示，那么上式写作
由于u是单位向量，即，上式两边都左乘u得，
即
就是的特征值，u是特征向量。

最佳的投影直线是特征值最大时对应的特征向量，其次是第二大对应的特征向量，依次类推。

因此，我们只需要对协方差矩阵进行特征值分解，得到的前k大特征值对应的特征向量就是最佳的k维新特征，而且这k维新特征是正交的。

得到前k个u 以后，样例通过以下变换可以得到新的样本。

其中的第j维就是在上的投影。

通过选取最大的k个u，使得方差较小的特征（如噪声）被丢弃。

三、算法性能评估
首先我们用k近邻算法对原始数据进行分类，得到识别率为0.7147.
其次，我们使用PCA算法对原始数据进行相应的降维处理，把200维的数据依次将至1——200维，观察其识别率，得到下面的识别率统计图。

图3.1将200维降至1——200维时的识别率统计
由上图可以看出，将数据维数降至40——80维时识别率已经趋于稳定，与数据不作处理时的识别率高度接近。

由此也证明了PCA算法的准确性。

下图给出了上图在1——20维时的放大图
由上图可以看出，将数据压缩的维数过低的话，识别率会出现较大波动。

但随着维数增加，相应的识别率也会趋于稳定。

下图反映了维数降至20--70维时的识别率波动图。

图3.3将200维降至20—62维时的识别率统计
由上图可以看出，当维数降至40-70维时，识别率波动不超过0.02，这是可以接受的一个数据。

再次说明了，200维数据的地物描述可以用更低维的数据来描述，既可以减少计算量，又可以保证较好的识别率。

四、用分类图证明算法的合理性
20406080100120140
图4.1数据的标准RGB图
上图给出标准地物图的RGB图，而下图则是未经PCA处理，直接得出的分
类图，数据维数为200维，我们把它记为对照图。

20
40
60
80
100
120
140
20406080100120140
图4.2未经PCA 处理的分类图（对照图）
20
40
60
80
100
120
140
由图4.3可以看出，将数据做PCA处理，将200维压缩成200维，其识别效果
20
40
60
80
100
120
140
20
40
60
80
20406080100120140
图4.5降至60维的分类图
20
40
60
80
20
40
60
80
由以上图4.4—4.7，我们不难看出将位数降至40-80维时，识别效果较之对照
图还是相当可观的，可接受的。

如果将维数降的太低，则会出现识别率过低的情况。

下图就展示降至低维时的情形。

20
40
60
80
100
120
140
20
40
60
80
100
120
140
五、结论
通过对降至各维时的识别率的统计与相应的识别后的RGB图，可以知将200维的数据降至45—80维时就能达到较好的识别率，这样既有效的节约了资源，同时也保证了识别率。

在近几年来，数据迅速膨胀并变大，它决定着企业的未来发展，虽然现在企业可能并没有意识到数据爆炸性增长带来问题的隐患，但是随着时间的推移，人们将越来越多的意识到数据对企业的重要性。

大数据时代对人类的数据驾驭能力提出了新的挑战，也为人们获得更为深刻、全面的洞察能力提供了前所未有的空间与潜力。

在大数据时代，数据的挖掘技术显得尤为重要。

这就需要大家潜心研究，找出更好更优的数据挖掘技术。

六、部分代码展示
6.1降至各维时的总识别率统计
function my_tu=zhixiantu(from,to)
load('C:\IndianaPines.mat')
Label=double(Label);
Class1=[pixels;Label'];
zer=0;
for i=1:size(Class1,2)
if Class1(201,i)~=0
zer=zer+1;
a(:,zer)=Class1(:,i);
end
end
%%归一化
for i=1:size(a,1)-1
a(i,:)=(a(i,:)-min(a(i,:)))/(max(a(i,:))-min(a(i,:)));
end
class=a';
%%求DataAjust
for j=1:size(class,2)-1
u=mean(class(:,j));
DataAjust(:,j)=class(:,j)-u;
end
d=cov(class(:,1:size(class,2)-1));
[b2,b1]=eig(d);
eigvalue=eig(d);
vector=b2;
[e,f]=sort(eigvalue);
%k=4;%说明压缩维数
for k=from:to
for i=1:k
Eigvector(:,i)=vector(:,f(201-i));
end
final=DataAjust*Eigvector;
final(:,k+1)=class(:,201);
w(:,1)=final(:,k+1);
for Tencount=1:3
c=randperm(size(final,1));
m=1;
p=1;
for category=1:16
n=0;
for j=1:size(c,2)
if w(c(j))==category
n=n+1;
b(n)=c(j);
end
end
d=floor(size(b,2)*0.1);
train(m:m+d-1,:)=final(b(1:d),:);
m=m+d;
%train(201,:)=k;
%test(:,p:p+size(b,2)-d-1)=a(:,b(d+1:size(b,2)));
%p=p+size(b,2)-d;
%test(201,:)=k;
b=[];
end
%%分类
test=final;
kmeans=3;%指定近邻方式
for x=1:size(test,1)
for y=1:size(train,1)
%%
distance=0;
for k1=1:k
distance=distance+(test(x,k1)-train(y,k1)).^2;
end
dis(y)=distance;
end
[z,v]=sort(dis);
p1=zeros(1,16);
for n=1:kmeans
for q1=1:16
if train(v(n),k+1)==q1
p1(q1)=p1(q1)+1;
end
end
end
[z1,v1]=max(p1);
test(x,k+2)=v1;
dis=[];
end
%%统计识别率
count=0;
for i=1:size(test,1)
if test(i,k+1)==test(i,k+2)
count=count+1;
end
end
lv(Tencount)=count/10366;
train=[];
end
lv1(k-19)=mean(lv);
mean(lv)
k
%train=[];
end
k=20:62;
lv2=lv1;
plot(k,lv2,'b-','Marker','o','LineWidth',2); hold on;
m2p=linspace(0,20,28);
plot(m2p,spline(k,lv2,m2p),'r-','LineWidth',2); xlabel('把200维降至1-20维');ylabel('相应的识别率'); grid on;
title('识别率统计图（续）');
end
6.2关于判定类别
function h=RGB(k)
load('C:\IndianaPines.mat')
Label=double(Label);
Class1=[pixels;Label'];
a=Class1;
w=Label';%类型不一致
%%去除类标为零的，但保存类别号。

画图用。

zer=0;zer0=0;
for i=1:size(a,2)
if a(201,i)~=0
zer=zer+1;
pa(:,zer)=a(:,i);
zn0(zer)=i;
else
zer0=zer0+1;
z0(zer0)=i;
end
end
%%归一化
for i=1:size(a,1)-1
pa(i,:)=(pa(i,:)-min(pa(i,:)))/(max(pa(i,:))-min(pa(i,:))); end
%%PCA处理
class=pa';
zer=0;zer0=0;
for i=1:size(a,2)
if a(201,i)~=0
zer=zer+1;
pa(:,zer)=a(:,i);
zn0(zer)=i;
else
zer0=zer0+1;
z0(zer0)=i;
end
end
for j=1:size(class,2)-1
u=mean(class(:,j));
DataAjust(:,j)=class(:,j)-u;
end
d=cov(class(:,1:size(class,2)-1));
[b2,b1]=eig(d);
eigvalue=eig(d);
vector=b2;
[e,f]=sort(eigvalue);
%k=200;%说明压缩维数
for i=1:k
Eigvector(:,i)=vector(:,f(201-i));
end
final=DataAjust*Eigvector;
final(:,k+1)=class(:,201);
w=final(:,k+1);
%%选择训练、测试样本
%Avq=zeros(1,16);
test=final;
result=zeros(size(test,1),16);
for tenCount=1:5
disp(sprintf('第%d次识别\n',tenCount));
c=randperm(10366);
m=1;
p=1;
for category=1:16
n=0;
for j=1:size(c,2)
if w(c(j))==category
n=n+1;
b(n)=c(j);
end
end
d=floor(size(b,2)*0.1);
train(m:m+d-1,:)=final(b(1:d),:);
m=m+d;
b=[];
end
%%分类开始
kmean=3;%指定近邻方式
for x=1:size(test,1)
%for x=1:200
for y=1:size(train,1)
%%
distance=0;
for k1=1:k
distance=distance+(test(x,k1)-train(y,k1)).^2;
end
dis(y)=distance;
end
[z,v]=sort(dis);
p1=zeros(1,16);
for n=1:kmean
for q1=1:16
if train(v(n),k+1)==q1
p1(q1)=p1(q1)+1;
end
end
end
[z1,v1]=max(p1);
test(x,k+1+tenCount)=v1;
result(x,v1)=result(x,v1)+1;
end
train=[];
end
%%十次分类中次数较多者为最终分类
for x=1:size(test,1)
[z2,v2]=max(result(x,:));
test(x,k+1+tenCount+1)=v2;
end
%%导出分类标签
mylabel(zn0,1)=test(:,k+1+tenCount+1); mylabel(z0,1)=0;
h=mylabel;
end。