遥感图像的光谱特征提取与比较方法

摘要

一种遥感图像的光谱特征提取与比较方法，用主分量变换完成多光谱图像谱间去相关以后，确定图像的各聚类中心矢量，按照“就近原则”将图像的各矢量元素分到各聚类中心代表的类中去，然后对聚类的结果进行合并与分裂操作，使得聚类的结果趋向平衡，得到特征提取的最终结果，最后在两幅图像完成特征提取后，按四个特征比较标准将矢量数组之间的比较进行标量化，使得图像间的相似和相异程度得到定性和定量的体现。本发明采用改进的迭代自组织数据分析技术ISODATA，在保证原算法效果的情况下，尤其在遥感图像波段较多时，能大幅度降低计算量，减少所需要的参数，使参数的控制更容易，并且使特征提取结果可以用于图像数据库检索。

说明书

技术领域：

本发明涉及一种遥感图像的光谱特征提取与比较方法，尤其涉及一种基于改进ISODATA(Iterative Self-Organizing Data Analysis Technique，迭代自组织数据分析技术A)算法的遥感图像的光谱特征提取与比较方法，在基于内容的图像数据库检索、图像分类和模式识别中均可有广泛运用。

背景技术：

基于内容的图像检索就是通过分析图像的内容(光谱、纹理、形状等)对图像进行检索的技术。其特点是：1、检索过程具有交互性，用户可以参与检索过程；2、引入了特征库和知识辅助的概念；3、它更注重的是信息的快速查询。

对遥感图像进行光谱特征提取以及特征比较是对遥感图像数据库进行检索的重要组成部分，由陈华等提出的基于主分量变换(K-L变换，Karhunen Loeve变换)和ISODATA聚类算法的多光谱遥感图像特征提取方法(陈华，安斌，陈书海，刘永昌.KL变换在多光谱图像聚类中的运用.红外与激光工程，2002，30(2)：79-82.)取得了不错的效果，但是，由于多光谱图像的数据量很大，使得该算法的计算量很大，另一方面，该算法中在ISODATA聚类中的参数过多，而且难以控制，同时该算法所提取的特征没有进行比较，无法运用于图像数据库检索。

发明内容：

本发明的目的在于针对现有技术存在的不足，提供一种基于改进ISODATA算法的遥感图像光谱特征提取与比较的方法，能够在保证原算法效果的情况下，(尤其在遥感图像波段较多时)大幅度的降低计算量，减少所需要的参数，使参数的控制更容易，并且提供特征比较的方法，使算法可以用于图像数据库检索。

为实现这样的目的，本发明的创新点在于保留原算法在初始聚类中得到的成果，把叠代的范围缩小，并且让叠代自然结束；将原算法的绝对性质的固定参数改为比例性质的动态参数；设计特征比较算法，使得图像间的相似和相异程度得到定性和定量的体现。用主分量变换完成多光谱图像谱间去相关以后，确定图像的各聚类中心矢量，按照“就近原则”将图像的各矢量元素分到各聚类中心代表的类中去。对聚类的结果进行合并与分裂操作，使得聚类的结果趋向平衡，得到特征提取的最终的结果。在两幅图像完成特征提取后，对它们的结果进行特征比较，得到它们的差异标准。

本发明实现方法主要包括图像矢量元素的初始聚类、类的合并与分裂、特征比较三个基本步骤：

1、图像矢量元素的初始聚类

先将原多光谱遥感图像进行主分量变换，根据实际需要的精度确定所要保留的主分量的个数，按照主分量对应的特征值由大到小的顺序保留主分量，得到主分量矩阵。令初始聚类中心数目等于主分量个数的平方，计算图像主分量矩阵的均值和标准差，以均值为中心，在左右2倍标准差之间等分，得到端点作为初始聚类中心。对主分量矩阵中的各矢量元素以“就近原则”进行初始聚类，即求图像的每个矢量元素与各聚类中心的二范数距离，若某矢量元素与某聚类中心的距离最小，就将该矢量元素划分到该类。最后得到每个类的元素集合与元素数目。

2、类的合并与分裂

计算元素数目的百分比，将百分比小于某一给定阈值的类去掉。计算各类的均值与标准差以及总体的标准差。如果某一类的标准差与总体标准差的比值大于某一给定阈值，将该类分解，即以该类的聚类中心与标准差的和与差形成一对新的聚类中心，然后将该类元素按照“就近原则”划分到新类中。如果某一类的标准差与总体标准差的比值小于某一给定阈值，将它与下一个和它相邻的类合并，计算合并以后的类的均值和标准差作为新类的聚类中心和标准差，但如果是最后一个类出现此种情况，就不用作任何处理。重复计算各类均值与标准差以及总体的标准差并进行以上判断，直到没有类的合并与分裂发生为止。最后形成的聚类中心、各类元素数目、各类标准差即为特征提取结果。

3、特征比较

将两幅图像的聚类中心和标准差分别标量化，各类元素数目进行归一化，然后按以下四个标准进行比较。标准1：以聚类中心为横坐标，各类元素百分比为纵坐标在同一个坐标系下作曲线，定性比较它们的相似程度；标准2：计算曲线的积分面积以及它们差的绝对值，得到面积差与面积之间的比值；标准3：将两幅图像的聚类中心以“就近原则”配对，并计算它们对应的原矢量之间二范数的平均值；标准4：计算由标准3配对的类的标准差对应的原矢量之间的二范数，计算两聚类中心标准差之间的距离与它们的几何平均值的比值。

本发明的光谱特征提取与比较方法具有如下有益效果：

通过保留ISODATA聚类算法在初始聚类中得到的成果，把叠代的范围缩小为类的合并与分裂，并让叠代自然结束，使得在保证原算法效果的情况下，尤其在遥感图像波段较多时，大幅度的降低了计算量；通过将原算法的绝对性质的固定参数改为比例性质的动态参数，使得在保证原算法效果的情况下，所需要的参数减少，参数控制更容易，参数意义更加明确；通过设计四个特征比较标准，将矢量数组之间的比较进行标量化，使得图像间的相似和相异程度得到定性和定量的体现，并使得本发明方法的特征提取结果可以用于图像数据库检索。

附图说明：

图1为本发明基于改进ISODATA算法的遥感图像光谱特征提取与比较方法的流程图。

如图所示：对原图像进行主分量变换完成以后，确定各聚类中心矢量；按照“就近原则”将图像的各矢量元素分到各类中去；对聚类的结果进行合并与分裂，使聚类结果趋向平衡；对聚类结果进行特征比较，得到四个差异标准；给出两幅图像相似程度的定性和定量说明。

图2为本发明实施例采用的遥感图像。

图2所示为“数字上海第一次飞行”遥感图像数据库中三幅200×200×128的图像，图中显示均是第120波段，背景均是黄浦江边陆地。

图3为图2中三幅图比较的标准1曲线，图3(a)为图2(a)与(b)比较曲线，图3(b)为图2(a)与(c)比较曲线，图3(c)为图2(b)与(c)比较曲线。

具体实施方式：

为了更好地理解本发明的技术方案，以下结合附图和实施例对本发明的实施方式作进一步描述。

本发明的一个实施例采用的原始遥感图像为“数字上海第一次飞行”遥感图像数据库中三幅200×200×128的图像，如图2所示，图中显示均是第120波段，背景为黄浦江边陆地。本发明遥感图像的光谱特征提取与比较方法流程如图1所示，具体实施方法按如下步骤进行：

1、图像矢量元素的初始聚类

先对原图像进行主分量变换，去除谱间相关性。将多光谱图像每个波段记为灰度矩阵向量形式：f_i＝[α₁₁...α_1nα₂₁...α_2nα₃₁......α_m1...α_mn]，则B＝[f₁f₂......f_p]¹表达一幅图像，B的协方差矩阵C的特征向量ν_k按特征值λ_k从大到小对应排列，组成变换阵A，C_D＝A′CA，C_D为对角阵。线性变换D＝A′B，D的行称为主分量，λ_k是主分量方差。若主分量较多，去掉相对较小的特征值对应的主分量，对结果没有太大影响，却能减小计算量。本实施例中保留了前8个主分量，即主分量为200×200×8的3维矩阵。

对主分量矩阵进行初始聚类。ISODATA算法的初始聚类中心数目K是根据经验确定的，一般取主分量数p的平方。本发明中K的初值也这样确定，但是，本发明方法中步骤2的迭代是自然结束的，迭代结果是个平衡点，则K对结果的影响大大降低，算法本身的鲁棒性得到加强，即K在p²一定范围内取值，如p²＝9，在取K＝8，9，10时，结果变化很小。计算D的均值向量u和标准差向量σ：

矩阵D的均值向量：

$>>u>=>>1>mn>>>\Sigma >>l>=>1>>mn>>>x>l>>->->->>(>1>)>>>s>$

矩阵D的标准差向量：

$>>\sigma >=>>>1>>mn>->1>>>>\Sigma >>l>=>1>>mn>>>>(>>x>l>>->u>)>>2>\; >->->->>(>2>)>>>s>$

在(u-2σ)和(u+2σ)间等分(K-1)份，得K个初始聚类中心center(i)，将矩阵D所有的元素以“就近原则”聚类。设x为一向量元素，min是求最小值，norm2是求二范数，Index是求序号。如下所示，x划分到第r类：

$>>center>>(>i>)>>=>>(>u>->2>\sigma >)>>+>>>4>\sigma >>>K>->1>>>>(>i>->1>)>>,>i>=>1,2>,>.>.>.>.>.>.>K>->->->>(>3>)>>>s>$

r＝Index(min(norm2(x-center(1)，norm2(x-center(2)，......，norm2(x-center(K))))    (4)

对每个元素重复以上操作，得到每个类的元素集合M_i和数目N_i。本实施例形成了64个类，每个类中有数目不等的矢量元素。

2、类的合并与分裂

如果某类满足

均值：

$>>aver>>(>i>)>>=>>1>>N>i>>>>\Sigma >>l>=>1>>>N>i>>>>M>i>>>(>l>)>>->->->>(>6>)>>>s>$

标准差极大值：

$>>stda>>(>i>)>>=>max>>(>>>1>>>N>i>>->1>>>>\Sigma >>l>=>1>>>N>i>>>>>(>>M>i>>>(>l>)>>->aver>>(>i>)>>)>>2>\; >)>>->->->>(>7>)>>>s>$

标准差极大值均值：

$>>av>\_>stda>=>>1>W>>>\Sigma >>i>=>1>>W>>stda>>(>i>)>>,>i>=>1,2>,>.>.>.>.>.>.>W>->->->>(>8>)>>>s>$

设std为标准差运算，类的合并如下所示：

若stda(i)＜γ_c(av_stda)但是i＝W，不做处理；

若stda(i)＜γ_c(av_stda)并且i≠W，则M_i+1＝M_i+M_i+1，N_i+1＝N_i+N_i+1，

$>>aver>>(>i>+>1>)>>=>>1>>>N>i>>+>>N>>i>+>1>>>>>>(>>\Sigma >>l>=>1>>>N>i>>>>M>i>>>(>l>)>>+>>\Sigma >>l>=>1>>>N>>i>+>1>>>>>M>>i>+>1>>>>(>l>)>>)>>,>stda>>(>>M>>i>+>1>>>)>>=>std>>(>>M>>i>+>1>>>)>>->->->>(>9>)>>>s>$

抛弃aver(i)、stda(i)、M_i和N_i，γ_c＜1，γ_c越接近1就越倾向合并。

类的分裂运算如下所示：

若stda(i)＞γ_s(av_stda)

则aver(W+1)＝aver(i)-stda(i)，aver(W+2)＝aver(i)+stda(i)        (10)

以“就近原则”将M_i中的元素分到M_W+1与M_W+2中，并计算N_W+1、N_W+2和stda(W+1)、stda(W+2)，抛弃aver(i)、stda(i)、M₁和N₁。γ_s＞1，γ_s越接近1就越倾向分裂。将新的聚类中心以“就近原则”插入到原聚类中心数组里。

根据新的聚类中心，重复步骤2，直到没有类的合并与分裂为止，最终形成三个特征量：聚类中心矢量组aver(i)、各类标准差矢量组stda(i)、各类元素数目N₁。本实施例的三幅图各形成34、36和34个类，以及各类对应的聚类中心和标准差。

3、特征比较

将aver1(i)，aver2(j)中维数大的后几维去掉，使得维数相同，N1_i，N2_j归一化。

令r＝norm2(aver1(1)-aver2(1))，s＝norm2(aver1(1)-aver2(2))，t＝norm2(aver2(1)-aver2(2))，

  u＝norm2(aver1(m)-aver2(n))，v＝norm2(aver1(m)-aver2(n-1))，

  w＝norm2(aver2(n)-aver2(n-1))

若r＜s并且t＜s，min＝aver1(1) 否则min＝aver2(1)，

若u＜v并且w＜v，max＝aver1(m) 否则max＝aver2(n)。

aver1(i)＝norm2(aver1(i)-min)/norm2(max-min)，

aver2(i)＝norm2(aver2(i)-min)/norm2(max-min)

stda1(i)和stda2(j)同理处理。

标准1：以aver1(i)，aver2(j)为横坐标，N_i，N2_j为纵坐标作图g₁，g₂定性比较；

标准2：

$>>q>=>abs>>(>S>)>>/>>>S>1>>>S>2>\; >->->->>(>11>)>>>s>$

S₁、S₂和S分别是g₁、g₂和g₁-g₂的积分值，q越小，图像越相似；

标准3：Index(i，j)，两幅图像的aver1(i)、aver2(j)以“就近原则”配对，并计算它们对应的原矢量之间距离的平均值，平均值越小图像越相似；

标准4：Qstda(i，j)＝mean(stda1(i)-stda2(j))                       (12)

其中i，j来自Index(i，j)，mean表示求平均，Qstda(i，j)越小图像越相似。

对图2所示的3幅遥感图像进行特征比较如下所示：

(1)当参数α＝0.01，γ_s＝0.5，γ_c＝2时，图2特征比较标准1见图3。(图3(a)中的实线表示图2(a)，虚线表示图2(b)；图3(b)中实线表示图2(a)，虚线表示图2(c)；图3(c)中实线表示图2(b)，虚线表示图2(c))。

       表1：当参数α＝0.01，γ_s＝0.5，γ_c＝2时，标准2、3、4结果

图2(a)与图2(b)图2(a)与图2(c)图2(b)与图2(c)  面积比0.660480.851340.87367  聚类中心平均距离5.22229.77277.6111  标准差平均距离9.121510.9099.0091

根据4个标准的定义，从图3和表1得出判断：图2(a)和图2(b)的灰度分布范围和分布规律都比较相似，图2(a)和图2(c)的以及图2(b)和图2(c)的相差较大。(2)采用原ISODATA算法以及上文所用的比较标准，使用参数：类元素数目最小值参数θ_N＝400、类分裂参数θ_S＝56.9、类合并参数θ_C＝151.4、L＝3，θ_N是原图像总象素的1％，θ_S是原图像初始聚类中心之间最大距离的5％，θ_C是首次初始聚类后标准差平均值中的最大值。

表2：原算法和本发明的改进算法在Matlab6.5下完成聚类所需时间(CPU主频：866MHZ，单位：秒)

图2(a)与图2(b)图2(a)与图2(c)图2(b)与图2(c)改进算法仿真(1)聚类时间65.33466.37566.967原算法仿真(2)聚类时间254.33337.69245.34

从表2可以发现本发明的改进算法与原算法所用时间是大相径庭，原算法所用时间是改进算法所用时间的3倍左右，可见本发明的方法在计算量上确实比原算法优越。