本发明公开了一种高光谱图像监督分类方法，主要解决现有方法对高光谱图像像素之间的空间信息利用不充足的问题，其分类过程为：对于待分类的高光谱图像，依据训练样本集建立二叉树模型，得到内部节点索引表；依据该索引表，依序找出本次待分类的内部节点标号和标号对应的高光谱图像像素点的位置；逐个计算相应位置处的每个像素分别属于本次内部节点的左、右子节点的似然概率值，得到高光谱图像的初始类标记；依据交替优化算法，更新初始类标记得到最终类标记；分别记录最终类标记中的两个子节点包含像素点的位置，存储到相应的子节点中；依次处理每个内部节点对应的高光谱图像数据，直至处理完所有内部节点。本发明可用于遥感图像分类与识别。
(1)对于待分类的高光谱图像，依据训练样本集，建立二叉树模型，得到节点标号，将内部节点标号及其包含的类别组存储于索引表中：
1a)计算训练样本集中各个类别的均值向量，将其作为各类别特征；
1b)使用欧式距离作为衡量各类别之间差异的度量，计算出每两个类别i1，i2之间的空间距离diff(i1，i2)，所有类别之间的空间距离构成一个大小为K×K的可分性判别矩阵D，K为类别总数；
1c)找出可分性判别矩阵中最大值max(D)对应的两个类别i，j，将其作为两个类别组C1和C2的初始值，即C1＝{i}，C2＝{j}；逐个比较每个剩余类别i′与类别i和类别j的空间距离，若i′与类别i的空间距离diff(i′，i)小于i′与类别j的空间距离diff(i′，j)，则将类别i′划到第一类别组C1中；否则将i′划到第二类别组C2中，即下式：其中i′∈{1，2，...，K}，i′≠i，j，U表示所有i′的集合；
1d)按照1c)中的方法，不断将每个类别组细分为更小的类别组，直至所有类别组都只包含一个最终类别，树模型即建立成功；树模型建立过程中的每个类别组以及最终类别统称为树的节点，其中类别组称为内部节点，最终类别称为叶子节点；给所有节点分配节点标号，树的根节点标号为1，该根节点标号1对应整个原始图像，包含图像中所有像素点的位置；对树中任一个节点t，其左子节点l(t)的节点标号为2t，右子节点r(t)的节点标号为
2t+1，将内部节点标号及其包含的类别组存储于索引表中；
(2)将高光谱图像数据和树节点对应，依据数模型中的内部节点索引表，依序找出本次待分类的内部节点标号和标号对应的高光谱图像像素点的位置；
(3)根据得到的高光谱图像像素点位置，逐个计算相应位置处的每个像素分别属于本次内部节点的左、右子节点的似然概率值，将两个似然概率值中最大值对应的子节点的类别作为该像素类别，由此得到高光谱图像的初始类标记；
3a)对于节点t，它所包含的像素点在空间上独立，即每个像素点的分布与其它点不相关，则节点的似然概率值 表示成各个像素似然概率 的乘积，如下式所示：
其中y 和x 分别为节点t对应的像素集合和相应的类标记集合， 和 分别为y 和x中位于位置s处的像素和类标记，w(t)为节点t的祖先节点集，指树中先于节点t产生的所有节点的集合， 为祖先节点集w(t)的类标记估计， 为 中位于位置s处的类标记估计， l(t)和r(t)分别为节点t的左、右子节点；
3b)将节点的似然概率值分解为各个像素似然概率值的乘积之后，对于节点t中每一个像素点 和 分别表示该像素属于节点t的左子节点l(t)和右子节点r(t)的似然概率，在服从高斯分布的条件下，似然概率的计算式如下：和 分别为节点t中的像素 所属的类别 对应的图像数据的均值向量和协方差矩阵，B为高光谱图像的维数，T为取转置操作，(·)-1为取逆操作，|·|为取行列式操作；
3c)由<3>式计算出似然概率后，对于节点t中每一个像素点 若该像素属于节点t的左子节点l(t)的似然概率 大于属于右子节点r(t)的似然概率则该像素的类标记为左子节点l(t)的类别；否则，像素的类标记为右子节点r(t)的类别；
(4)依据交替优化算法，将得到的初始类标记进行类标记更新，得到本次的最终类标记：t t
4a)对节点t的类标记x 的更新和参数β 的估计交替进行，令k＝0， 由下式计算 和t t
其中 和 分别为节点t对应数据的类标记x 和边缘罚函数β 在交替优化算法的第k次估计中的结果， 为节点t中，给定祖先节点集合的类标记估计先验分布在交替优化算法的第k次估计中的结果； 为节点t中，给定类标t记x 以及祖先节点集合的类标记估计 似然分布在交替优化算法的第k次估计中的结果；
4b)令k＝k+1，已知第k-1次类标记估计结果 则第k次的边缘罚函数的估计用极大伪似然估计算法MPL计算，表达式如下：令 在设定的经验值范围[a，b]内，以固定步长λ增加 则使得<5>式等t
号右边的表达式值最大的 即MPL估计结果；其中t为当前节点标号，S 为节点it对应像素点的位置集，s为S 中的一个像素点的位置，η(s)为位置s的邻域，r为邻域η(s)内任一点， 为位置s的邻域中与其具有相同类标记的像素个数，L为标记空间；
4c)将估计出的参数 带入<4>式，应用迭代条件模式ICM算法逐点更新图像中像素点的类标记 具体计算公式如下：和 分别为节点t中的像素 所属的类别 对应的图像数据的均值向量和协方差矩阵， 为位置s的邻域中与其具有相同类标记的像素个数；
4d)将本次优化后的类标记 与上一次的类标记 进行比较，若两者相比变化不大，即 则得到的类标记 就是最终类标记优化结果 即 结束交替优化算法；否则重复步骤4b)-4c)；
(5)分别记录本次分类后的两个子节点包含像素点的位置，存储到相应的子节点中；
(6)重复步骤(2)到步骤(5)，依次处理每个内部节点对应的高光谱图像数据，直至处理完所有内部节点对应的数据，完成对高光谱图像的分类。
高光谱图像监督分类方法
技术领域
[0001] 本发明属于图像处理领域，涉及一种高光谱图像监督分类方法，可用于遥感图像分类与识别。
背景技术
[0002] 高光谱图像的分类与识别是遥感图像处理领域的一个新的重要分支，近年来在植被研究、矿物识别以及水文和海洋勘测等领域有着广泛的应用。高光谱遥感，是指具有很高的光谱分辨率的遥感技术，它能够在很窄的连续波段上记录地物，如地表、水域、大气的光谱反射特性，使得获取待测地物类别的细节成为现实。成像光谱仪能够为每个像素提供数十个至数百个窄波段，高光谱图像的每一帧都是对同一场景在不同波长上的成像，从而产生一条完整而连续的光谱曲线，提供所观测地物丰富的光谱信息。
[0003] 针对高光谱遥感技术的上述特点，常用的高光谱图像的分类方法有：
[0004] 从光谱角度出发，有SAM算法和混合像素分解技术。SAM方法通过比较目标光谱曲线和已有的样本库中的各类平均光谱曲线之间的相似程度来判定待分类光谱的类别归属；而混合像素光谱分解技术是通过将每个像素点分解为不同的基本单元即端元，来实现类别的划分。该类方法通常需要预先建立待分类图像所包含地物的光谱库，因而在实现上比较困难。
[0005] 从统计模型角度出发，有基于统计距离的方法和基于概率模型的方法，最小距离分类法是基于统计距离的方法中最简单的一种，该方法通过比较待分类像素和利用训练样本计算出的各类中心点之间的距离，来实现像素点的类别划分。而基于概率模型的方法中最典型的就是高斯极大似然分类方法，该方法将像素向量作为特征向量，在假定各类服从正态分布的基础上，根据有关概率判决函数的贝叶斯准则进行分类。该类方法应用广泛，但由于算法的复杂度和数据维数的平方成正比，且类别参数由于训练样本的相对较少无法可靠估计。
[0006] 还有基于决策树的分类方法，如二元层级分类器BHC，该方法假设高光谱数据中各类服从一定分布，将贝叶斯规则用在每层的决策中，不断的将一个K类分类问题转化为K-1个的二元分类问题，从而解决分类问题。该方法相比之前的方法，计算简单且易于实现，但由于没有很好的加入先验信息，忽略了图像中像素之间的相互关系，导致分类过程中对空间信息的利用不足，限制了分类精度的提高。
发明内容
[0007] 本发明的目的在于克服上述已有技术的不足，提供一种树状的马尔可夫随机场模型TS-MRF的高光谱图像监督分类方法，以有效地利用图像中像素之间的空间关系，提高分类精度。
[0008] 实现本发明目的的技术方案是：使用简洁高效的二叉树结构，结合马尔可夫随机场对图像像素之间的空间信息的利用，对高光谱图像进行监督分类，具体实现步骤如下：
[0009] (1)对于待分类的高光谱图像，依据训练样本集，建立二叉树模型，得到节点标号，将内部节点标号及其包含的类别组存储于索引表中：
[0010] 1a)计算训练样本集中各个类别的均值向量，将其作为各类特征；
[0011] 1b)使用欧式距离作为衡量各类别之间差异的度量，计算出每两个类别i1，i2之间的空间距离diff(i1，i2)，所有类别之间的空间距离构成一个大小为K×K的可分性判别矩阵D，K为类别总数；
[0012] 1c)找出可分性判别矩阵中最大值max(D)对应的两个类别i，j，将其作为两个类别组C1和C2的初始值，即C1＝{i}，C2＝{j}；逐个比较每个剩余类别i′与类别i和类别j的空间距离，若i′与类别i的空间距离diff(i′，i)小于i′与类别j的空间距离diff(i′，j)，则将类别i′划到第一类别组C1中；否则将i′划到第二类别组C2中，即下式：
[0013]
[0014] 其中i′∈{1，2，...，K}，i′≠i，j，U表示所有i′的集合；
[0015] 1d)按照1c)中的方法，不断将每个类别组细分为更小的类别组，直至所有类别组都只包含一个最终类别，树模型即建立成功；树模型建立过程中的每个类别组以及最终类别统称为树的节点，其中类别组称为内部节点，最终类别称为叶子节点；给所有节点分配节点标号，树的根节点标号为1，该根节点标号1对应整个原始图像，包含图像中所有像素点的位置；对树中任一个节点t，其左子节点l(t)的节点标号为2t，右子节点r(t)的节点标号为2t+1，将内部节点标号及其包含的类别组存储于索引表中；
[0016] (2)将高光谱图像数据和树节点对应，依据树模型中的内部节点索引表，依序找出本次待分类的内部节点标号和标号对应的图像像素点的位置；
[0017] (3)根据得到的图像像素点位置，逐个计算相应位置处的每个像素分别属于本次内部节点的左、右子节点的似然概率值，将两个似然概率值中最大值对应的子节点的类别作为该像素类别，由此得到图像的初始类标记；
[0018] (4)依据交替优化算法，将得到的初始类标记进行类标记更新，得到本次的最终类标记；
[0019] (5)得到最终类标记后，分别记录本次分类后的两个子节点包含像素点的位置，存储到相应的子节点中；
[0020] (6)重复步骤(2)到步骤(5)，依次处理每个内部节点对应的图像数据，直至处理完所有内部节点对应的数据，完成对高光谱图像的分类。
[0021] 本发明与现有的技术相比有以下优点：
[0022] 1)本发明由于采用了马尔可夫随机场，更好地利用了像素之间的空间依赖关系，使得图像的空间信息得到充分利用，提高了分类精度；
[0023] 2)本发明由于采用了二叉树状的结构，使得分类的过程可以递归地进行，简洁有效；
[0024] 3)本发明由于采用了一种可移植的树模型建立方法，使得基于TS-MRF模型的方法的树模型建立过程有章可依，给予树模型建立以有效的理论基础，增强该算法的适用性。
附图说明
[0025] 图1是本发明的流程图；
[0026] 图2是本发明所采用的高光谱图像三个波段合成的示意图；
[0027] 图3是本发明所建立的二叉树模型图；
[0028] 图4是本发明对图2所示高光谱图像的仿真结果图。
具体实施方式
[0029] 参照图1，本发明的具体实现步骤如下：
[0030] 步骤一，对于图2所示的待分类高光谱图像Y＝{y1，y2，...，yn}，n为图像的总像素数，依据从高光谱图像中得到的带有样本标记的训练样本集，按照逐次二元划分的方式建立如图3所示的二叉树模型，得到树模型中各个节点的节点标号，将其中的内部节点标号及其包含的类别组存储于索引表中，具体步骤如下：
[0031] 1a)计算训练样本集中各类别所包含的高光谱图像数据的均值向量，将其作为各类特征；
[0032] 1b)使用欧式距离作为衡量各类别之间差异的度量，计算出每两个类别i1，i2之间的空间距离diff(i1，i2)，所有类别之间的空间距离构成一个大小为K×K的可分性判别矩阵D，K为类别总数；
[0033] 1c)找出可分性判别矩阵D中最大值max(D)对应的两个类别i，j，将其作为两个类别组C1和C2的初始值，即C1＝{i}，C2＝{j}；逐个比较每个剩余类别i′与类别i和类别j的空间距离，若i′与类别i的空间距离diff(i′，i)小于i′与类别j的空间距离diff(i′，j)，则将类别i′划到第一类别组C1中；否则将i′划到第二类别组C2中，即下式：
[0034]
[0035] 其中i′∈{1，2，...，K}，i′≠i，j，U表示所有i′的集合；
[0036] 1d)按照1c)中的方法，不断将每个类别组细分为更小的类别组，直至所有类别组都只包含一个最终类别，树模型即建立成功；树模型建立过程中的每个类别组以及最终类别统称为树的节点，其中类别组称为内部节点，最终类别称为叶子节点；给所有节点分配节点标号，树的根节点标号为1，该根节点标号1对应整个原始图像，包含高光谱图像中所有像素点的位置S＝{1，2，..，n}；对树中任一个节点t，其左子节点l(t)的节点标号为2t，右子节点r(t)的节点标号为2t+1，将内部节点标号及其包含的类别组存储于索引表ID中。
[0037] 步骤二，将高光谱图像数据yt和树中的节点t对应，依据树模型中的内部节点索引表ID，依序找出本次待分类的内部节点标号ID(t)和标号对应的高光谱图像像素点的位t置S，除根节点之外的其它内部节点对应的高光谱图像像素点的位置均为原始图像位置集t
的一部分，即S ∈S。
[0038] 步骤三，根据得到的高光谱图像像素点位置St，逐个计算其中的每个位置s处的像素 分别属于本次内部节点t的左子节点l(t)和右子节点r(t)的似然概率值，将两个似然概率值中最大值对应的子节点的类别作为该像素的类别，由此得到高光谱图像的初始类标记，具体步骤如下：
[0039] 3a)对于节点t，它所包含的高光谱图像数据yt中的像素点 在空间上独立，即每个像素点的分布与其它点不相关，则节点的似然概率值 表示成各个像素似然概率 的乘积，如下式所示：
[0040]
[0041] 其中yt和xt分别为节点t对应的像素集合和相应的类标记集合， 和 分别为ytt和x 中位于位置s处的像素和类标记，w(t)为节点t的祖先节点集，指树中先于节点t产生的所有节点的集合， 为祖先节点集w(t)的类标记估计， 为 中位于位置s处的类标记估计， l(t)和r(t)分别为节点t的左、右子节点；
[0042] 3b)将节点的似然概率值分解为各个像素似然概率值的乘积之后，对于节点t中每一个像素点 和 分别表示该像素属于节点t的左子节点l(t)和右子节点r(t)的似然概率，在服从高斯分布的条件下，似然概率的计算式如下：
[0043]
[0044] 和 分别为节点t中的像素 所属的类别 对应的高光谱图像数据的均值向-1量和协方差矩阵，B为高光谱图像的维数，T为取转置操作，(·) 为取逆操作，|·|为取行列式操作；
[0045] 3c)由<3>式计算出似然概率后，对于节点t中每一个像素点 ，若该像素属于节点t的左子节点l(t)的似然概率 大于属于右子节点r(t)的似然概率则该像素的类标记为左子节点l(t)的类别；否则，像素的类标记为
右子节点r(t)的类别。
t
[0046] 步骤四，依据交替优化算法，将得到的初始类标记x 进行类标记更新，得到本次的最终类标记
[0047] 4a)对节点t的类标记xt的更新和参数βt的估计交替进行，令k＝0，由下式计算 和
[0048]t t
[0049] 其中 和 分别为节点t对应数据的类标记x 和边缘罚函数β 在交替优化算法的第k次估计中的结果， 为节点t中，给定祖先节点集合的类标记估计，先验分布在交替优化算法的第k次估计中的结果； 为节点t中，给定类t标记x 以及祖先节点集合的类标记估计 似然分布在交替优化算法的第k次估计中的结果；
[0050] 4b)令k＝k+1，已知第k-1次类标记估计结果 则第k次的边缘罚函数的估计 用极大伪似然估计算法MPL计算，表达式如下：
[0051]
[0052] 令 在设定的经验值范围[a，b]内，以固定步长λ增加 则使得<5>t式等号右边的表达式值最大的 即MPL估计结果；其中t为当前节点标号，S 为节点t对应像素点的位置集，s为St中的一个像素点的位置，η(S)为位置s的邻域，r为邻域η(s)内任一点， 为位置s的邻域中与其具有相同类标记的像素个数，L为标记空间；
[0053] 4c)将估计出的参数 带入<4>式，应用迭代条件模式ICM算法，逐点更新高光谱图像中像素点的类标记 具体步骤如下：
[0054] 4c1)针对每个像素点，利用初始分类得到的两类高光谱图像数据计算各自的均值向量 和协方差矩阵 和由<5>式估计出的边缘罚函数 一起，代入下面的<6>式进行优化：
[0055]
[0056] 使得该式最小的类标记即该像素点的类标记，由此将高光谱图像数据划分为两类，得到初始类标记；其中 和 分别为节点t中的像素 所属的类别 对应的高光谱图像数据的均值向量和协方差矩阵， 为位置s的邻域中与其具有相同类标记的像素个数；
[0057] 4c2)用步骤4c1)得到的新的两类高光谱图像数据计算各自的均值向量 和协方差矩阵 代入<6>式进行优化，得到新的类标记；
[0058] 4c3)检验是否满足迭代停止条件，即类标记发生变化的像素点总数小于预先设定的阈值，若满足停止条件，则结束ICM算法；否则重复4c1)-4c3)的步骤；
[0059] 4d)将本次优化后的类标记 与上一次的类标记 进行比较，若两者相比变化不大，即 则得到的类标记 就是最终类标记优化结果 即 结束交替优化算法；否则重复步骤4b)-4c)。
[0060] 步骤五，得到最终类标记 之后，分别记录本次分类后节点t的两个子节点l(t)l(t) r(t)和r(t)包含像素点的位置S 和S ，将其存储到相应的子节点中。
[0061] 步骤六，重复步骤二到步骤五，按照内部节点索引表ID中的节点存储顺序，依次t处理每个内部节点t对应的高光谱图像数据y，得到相应的节点标号 ，直至处理完所有内部节点对应的数据，完成对高光谱图像的分类。
[0062] 本发明的效果可通过以下仿真进一步说明：
[0063] 1.仿真内容
[0064] 应用本发明的方法和现有的高斯极大似然分类方法，对一幅高光谱图像进行分类实验，主要从各个类别的内部一致性、类别边缘的完整性和边缘定位准确度来评价这两种方法的性能。
[0065] 2仿真实验结果
[0066] 用本发明方法以及高斯极大似然分类方法对图2所示的高光谱图像进行分类，分类结果如图4所示。其中图4(a)为高斯极大似然分类的分类结果，图4(b)为本发明的分类结果，图4(c)为高斯极大似然分类的标记区域分类结果，4(d)为本发明的标记区域分类结果。由图4可见，本发明中，各个类别的一致性较为良好，类别边缘完整，且边缘定位准确，分类效果良好。相比之下，高斯极大似然分类方法各个类别的一致性也较为良好，但是类别边缘的定位很不准确，出现很多错误边缘，导致分类结果不佳。
