基于半监督学习的海岸海洋遥感影像特征地类的识别方法

摘要

本发明公开了一种基于半监督学习的海岸海洋遥感影像特征地类的识别方法，属于半自动遥感影像识别领域。其步骤为：为每一类特征地物选取标记样本；构建面向对象的遥感影像的分割结果；计算出所有样本像元隶属于各特征地类的初估概率值，计算出样本数据在归为各个特征地类分量的概率；使用特征空间规则对概率图像进行修正；判定其所属特征地类，实现特征地类的识别，并输出识别结果图。本发明结合了先验知识与数据的统计特性，能够用地学先验知识引导数据挖掘过程，实践证明，该算法够能有效地进行遥感影像分类，得到比较满意的结果，并具有高效率、高精度的特点，能够直接应用于国家各级基础地理信息数据库遥感专题信息的维护与更新。

说明书

技术领域

本发明涉及一种遥感影像数据的特征地类识别方法，特别是涉及一种基于半监督学习的海岸海洋遥感影像特征地类的识别方法。

背景技术

遥感作为一种非接触式的观测技术具有悠久的历史，其所提供的遥感影像数据具有时效性高、覆盖范围广、信息丰富等优点，已经被广泛地运用于土地利用、资源勘探、生态环境监测，以及海岸海洋特征地类识别等诸多领域，对社会经济的发展发挥着越来越重要的作用。海岸带区域作为社会和经济发展的核心地带，随着海洋经济开发力度的加强，海岸带区域的土地利用、海域使用及环境变化也日新月异。为了更好地保护和开发利用海洋资源，利用遥感技术的快速动态监测优势，能及时获取海岸带环境变化信息，并为海岸带区域的特征地类识别及资源环境监测提供了技术支撑。对一些敏感地区进行的周期性业务化遥感监测，无论在军事上，还是在海岸带环境资源开发利用管理上，都具有十分重要的意义。

随着对地观测技术的飞速发展，一个多层次、全天候、多分辨率互补的全球观测网络正在逐步形成，遥感影像数据源的种类与数量都呈现出几何级的增长，面对如此海量的遥感影像数据，现有的针对单一遥感数据源的人工目视解译与计算机辅助识别相结合的方法已远远不能满足多源、大范围、长时间序列遥感特征地类识别的需求。海岸海洋特征地类应用研究作为遥感应用的一个重要分支，也亟需发展面向多源遥感影像数据的、自动化程度较高的遥感信息识别方法。

遥感影像特征地类识别的主要目的是对基元进行特征分析，结合基元对象所属大区域的环境知识，构造相应的特征地类识别规则，把未知类别的样本映射到所构造的类别中，从而获得遥感影像与实际地物的对应信息。目前，遥感影像的特征地类识别以数据挖掘为技术手段。数据挖掘用于发现大量数据中所蕴含的知识或者是规则，在建立了正确的数据表示以后，利用适用于当前数据来对数据潜在的特征进行挖掘。数据挖掘中常用的模型有概率模型、决策树模型、神经网络模型、聚类分析模型、模糊集模型等。1997年，Huang等在期刊《Photogrammetry Engineering andRemote Sensing》63卷第10期中发表“A Machine-Learning Approach to Automated Knowledge Based Building for Remote Sensing Image Analysis with GIS Data”一文，用决策树方法从GIS数据和SPOT遥感影像中提取知识进行湿地分类；1998年，李祚泳在期刊《红外与毫米波学报》17卷第2期发表“用BP神经网络实现多波段遥感图像的监督分类”一文，研究了基于神经网络的多波段遥感影像监督分类方法；2006年，郭小卫等在期刊《遥感信息》第6期中发表“一种多尺度无监督遥感影像分割方法”一文，提出了一种基于聚类的多尺度无监督遥感影像分类方法，并利用该方法对SAR图像进行分类研究。然而在很多时候，单纯利用数据挖掘技术进行特征地类的识别，往往难以取得较好的效果。这是由于具有一定规模的自然区域，物质组成多样，地物在时空上相互交错，相互转换形成结构和形式各具特点的复杂混合体，导致同类地物光谱响应的多样化，使得“同物异谱，异物同谱”现象发生。

在遥感影像特征地物识别过程中，根据是否引入训练样本可分为监督学习、非监督学习和半监督学习方法。在遥感影像监督学习中，学习器通过对大量有标记的训练样本进行学习，从而建立模型用于预测未标记的图像单元地物类别。然而，从遥感影像中采集大量有标记的训练样本非常困难，并且需要耗费大量的人力物力；如果只使用少量的有标记示例，那么利用它们所训练出的学习系统往往很难具有强泛化能力。与之相对应的非监督学习中，学习器不依赖任何有标记的训练样本，直接对全部未标记的图像单元进行统计分析，建立相应的规则模型，并对这些图像单元进行地类识别。然而，由于完全没有先验知识的输入，非监督方法的结果很难直接作为成果数据产品。半监督学习是介于前两者之间的学习方式，，在少量已标记类别样本提供的监督信息的引导下，学习全部样本或只学习未标记类别样本。该方法只需提供少量的标记样本，而通过全部样本的学习又可以获得相对于非监督学习更好的学习效果。

半监督学习的基本设置是给定一个来自某未知分布的有标记示例集L＝{(x₁，y₁)，(x₂，y₂)，…，(x_|L|，y_|L|)}以及一个未标记示例集U＝{x₁′，x₂′，…，x_|U|′}，期望学得函数f：X→Y可以准确地对示例x预测其标记y。这里x_i，x_j′∈X均为d维向量，y_i∈Y为示例x_i的标记，|L|和|U|分别为L和U的大小，即它们所包含的示例数。

发明内容

1.发明要解决的技术问题

如前所述，单纯利用数据挖掘技术进行特征地类识别，往往难以取得较好的识别效果；采用监督学习方法根据已标记样本进行建模，则会造成大量未标记样本的浪费，且标记样本少则精度下降，标记样本多则成本上升；而运用非监督学习方法根据未标记样本建模，因为其未使用任何先验知识，故精度偏低，很难直接作为成果数据产品。

为此，本发明将遥感信息地学知识和数据挖掘方法结合起来，基于半监督学习的海岸海洋遥感影像特征地类的识别方法。该基于半监督学习的遥感影像地类识别方法选取少量已标记的样本来引导大量的未标记数据，并构建相应的模型模拟全体样本数据的联合分布概率，最终实现海岸海洋遥感影像特征地类识别，在精度和效率上均得到比较满意的特征地类识别结果，并具备良好的通用性。

2.本发明的技术方案如下：

原理：遥感影像识别是根据遥感影像中像元集的特征构造一个地类特征规则模型把未知类别的样本映射到相应的特征地类中，从而获得遥感影像中与实际地物的对应信息。本发明中采用高斯混合模型(Gaussian Mixture Model，简称GMM)捕捉遥感影像数据的特征，模拟出样本数据所服从的联合概率密度函数。其中，GMM的混合系数π_k、μ_k和∑_k是通过运用期望最大化方法(Expectation Maximization，简称EM)求解的极大似然估计值(Maximum Likelihood Estimate，简称MLE)。

基于半监督学习的海岸海洋遥感影像特征地类的识别方法，包括以下步骤(具体的步骤流程图见图1)：

步骤1：通过卫星遥感技术手段，获取描述地表物体对电磁波的反射和其自身发射的电磁波的遥感影像，根据不同类型地物在遥感影像中具有特定的光谱显示特征，为每一类特征地物在遥感影像中选取标记样本，构建一个长度为m的一维数组(m＝影像行数×影像列数×波段数)，依次读取各波段遥感影像数据并存入该数组；

步骤2：利用种子生长方法对经过预处理的遥感影像进行分割处理，得到一系列在光谱信息上表现相近似的像元，通过赋以相同属性值来合并这些表现相近的像元，并求取该对象中各像素光谱均值作为此对象的光谱信息，最终实现遥感影像的分割；

步骤3：初始化各类别密度分布待估计的参数Θ，包括各类别的比例、均值向量μ和协方差矩阵∑；根据步骤1中所选取的已标记样本进行最近邻分类，得到特征空间内的高维距离矩阵，进而计算出所有样本像元隶属于各特征地类的初估概率值；

步骤4：通过EM迭代算法对所有样本像元进行迭代运算，直到算法收敛，估计出高斯混合模型的参数θ_k＝{π_k，μ_k，∑_k}的值，并计算出样本数据在归为各个特征地类分量的概率；假设在迭代过程中，当前估计的参数为Θ^old，则通过一次E步和M步的计算后估计出的参数为Θ^new，引入辅助函数Q(Θ，Θ^old)表示每一次迭代后完整数据集的对数似然的数学期望，则E步和M步的实质分别为：

$Q(\Theta ,{\Theta }^{\mathrm{old}})=E\left[\ln L\left(\Theta \right|Y)\right]=E\left[\ln \underset{n=1}{\overset{N}{\Pi }}p(x_n,z_n|\Theta )\right]$①

Θ^new＝arg maxQ(Θ，Θ^old)    ②

具体E步和M步计算方法为：

1)期望步运算(Expectation Step，简称E步)：利用当前的参数估计结果计算所有样本的后验概率；具体通过公式③和④计算出样本数据向量是由第k个混合分量产生的概率，从数据挖掘的角度来看，即为样本x属于类别k的概率；

$p(z=k|x)=\frac{p(z=k)p\left(x\right|z=k)}{\underset{k=1}{\overset{K}{\Sigma }}p(z=k)p\left(x\right|z=k)}=\frac{{\pi }_{k}N\left(x\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x\right|{\mu }_k,{\Sigma }_k)}$③

$N\left(x\right|{\mu }_k,{\Sigma }_k)=\frac{1}{{\left(2\pi \right)}^{D/2}\sqrt{\left|{\Sigma }_k\right|}}\exp (-\frac{1}{2}{(x-{\mu }_k)}^T{\Sigma }^{-1}(x-{\mu }_k))$④

其中，D为样本数据向量x的特征维数(本发明中D＝1)，k为构成GMM的高斯模型的个数，而π_k即为先验概率p(z＝k)，N(x|μ_k，∑_k)即为条件概率p(x|z＝k)，概率密度函数为

2)极大似然步运算(Maximum Likelihood Step，简称M步)：利用当前所有样本的后验概率重新进行参数估计，并根据标记样本对参数进行加权叠加修正；具体通过公式⑤、⑥、⑦计算出高斯混合模型的参数θ_k＝{π_k，μ_k，∑_k}的值，其中，π_k为混合系数，表示一个随机抽取的样本是由第k个混合分量产生的概率，均值向量μ_k和协方差矩阵∑_k为每个独立的高斯模型N(x|μ_k，∑_k)自身的参数(见公式④)；

${\pi }_k=\frac{\underset{n=1}{\overset{N}{\Sigma }}\frac{{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}}{N}=\frac{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)}{N}$⑤

${\mu }_k=\frac{\underset{n=1}{\overset{N}{\Sigma }}\frac{{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}{x}_n}{\underset{n=1}{\overset{N}{\Sigma }}\frac{{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}}=\frac{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)x_n}{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)}$⑥

${\Sigma }_k=\frac{\underset{n=1}{\overset{N}{\Sigma }}\frac{{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}(x_n-{\mu }_k){(x_n-{\mu }_k)}^T}{\underset{n=1}{\overset{N}{\Sigma }}\frac{{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x_n\right|{\mu }_k,{\Sigma }_k)}}=\frac{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)(x_n-{\mu }_k){(x_n-{\mu }_k)}^T}{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)}$⑦

3)对E步和M步进行迭代计算，直到算法收敛；

步骤5：通过不同地物在遥感影像中的光谱特征，分别计算可见光指数(VisibleLight Index，简称VLI，公式⑧)、归一化植被指数(Normorlized Vegetation Index，简称NDVI，公式⑨)、归一化水体指数(Normorlized Water Index，简称NDWI，公式⑩)来的构建新的光谱空间规律，利用所构建的光谱空间规律对概率图像进行修正；

VLI＝BLUE+GREEN+RED              ⑧

NDVI＝(NIR-RED)/(NIR+RED)        ⑨

NDWI＝(GREEN-NIR)/(GREEN+NIR)    ⑩

式中，BLUE表示遥感影像的蓝光波段，GREEN表示遥感影像的绿光波段，RED表示遥感影像的红光波段，NIR表示遥感影像的近红外波段；

步骤6：比对所计算出的每个样本所隶属于各自类别的概率，判定所属概率最大的类别即为其所属特征地类，进行标记和赋值，实现特征地类的识别，输出识别结果图。

步骤1中改进目标识别前的样本选取，与传统的遥感影像目标识别方法(如最大似然法、神经网络法、决策树法等)需要大量样本进行统计分析不同，本算法仅要求为每个地类选取一个样本作为初始聚类中心，来替代标准EM算法中用随机生成法产生的初始聚类中心。

步骤4中改进目标识别中的迭代优化。在每轮迭代的M步中，用各地类的标记样本同模型估计的各聚类中心进行加权叠加，使聚类中心向标记样本漂移。该步骤确保了聚类中心不会与实际情况(依靠先验知识精选的样本)偏离太远，同时在一定程度上增加了算法的随机性，使其陷入局部最优的概率降低。

步骤5中改进目标识别后的结果修正；计算出各地类的概率分布后，通过特征空间中隐含的遥感信息地学知识对结果进行修正，这一步主要依赖于遥感影像特征空间的构建。

3.有益效果

相比现有技术，本发明将遥感信息地学知识和数据挖掘方法结合起来，并合理建立有标记样本和未标记样本之间的联系，提高了目标识别的性能。具体有益效果如下：

第一，改进目标识别前的样本选取。与传统的遥感影像地类识别方法(如最大似然法、神经网络法、决策树法等)需要大量样本进行统计分析不同，本方法仅要求为每个地类选取一个样本作为初始聚类中心，来替代标准EM方法中用随机生成法产生的初始聚类中心，提高识别性能。

第二，改进目标识别中的迭代优化。在每轮迭代的M步中，用各地类的标记样本同模型估计的各聚类中心进行加权叠加，使聚类中心向标记样本漂移。该步骤确保了聚类中心不会与实际情况(依靠先验知识精选的样本)偏离太远，同时在一定程度上增加了方法的随机性，使其陷入局部最优的概率降低。

第三，改进目标识别后的结果修正。计算出各地类的概率分布后，通过特征空间中隐含的遥感信息地学知识对结果进行修正，这一步主要依赖于遥感影像特征空间的构建。

综上，本发明结合了先验知识与数据的统计特性，能够用地学先验知识引导数据挖掘过程，实践证明，该方法够能有效地进行遥感影像目标识别，得到比较满意的结果，并具有高效率、高精度的特点。

附图说明

附图1为半监督遥感图像特征地类识别流程图；

附图2为实施例中所用卫星遥感影像的分割影像；

附图3为水体地类概率图像(白色到黑色渐变，代表概率由高到低)。

附图4为植被地类概率图像(白色到黑色渐变，代表概率由高到低)。

附图5为建筑物地类概率图像(白色到黑色渐变，代表概率由高到低)。

附图6为半监督学习特征地类识别结果图像。

具体实施方式

本实施例采用附图2所示遥感影像作为海岸海洋特征地类识别的初始影像，该遥感影像数据是由Landsat 5号卫星的TM主题成像仪于2007年1月2日所拍摄。实施例区域的图像大小为500×500像元，该区域位于海南岛西北部的新英湾红树林自然保护区。本实例具体实施采用标准C++编程语言在VC 6.0平台下实现，遥感影像数据的读写操作通过开源地理数据格式转换类库GDAL 1.60实现。具体实施步骤如下：

步骤1：以GDAL为影像数据读写工具，利用GDAL.Open方法读取遥感影像。根据海南岛的海洋海岸特征，实施中将地物划分为清洁水体、浑浊水体、植被、滩涂、建设用地、裸地和其他共7类，对应的这7类地物在遥感图像中各自选取相应的一个样本点(卫星遥感图像是对地物的实际探测与显示，每个地物在图像中将会有对应的像素进行描述，为此，可通过实地勘察并在其对应的遥感图像中选取相应的样本点)，本实例中所选取的相应的样本点在遥感图像中的(X，Y)空间坐标值分别为(151，248)、(457，208)、(330，274)、(185，55)、(295，266)、(113，189)、(130，437)。构建一个长度为1500000的一维数组(500行×500列×6个波段)，依次读取各波段遥感影像数据并存入该数组。

步骤2：①初始化所有像素点的区域对象属性，将整幅图像所有像素点属性初始化为-1(表示像素点还没有进行区域标号)；②以遥感影像的所有像素作为种子点标记对象，依次遍历遥感影像各像素，根据该像素点与其邻接的8个像素点的相似度来进行像素点的归类，其相似度的衡量准则以该像素点与其相邻的各像素点是否存在差异为基准，如果没有差别，则标记为同一区域；③对所标记的区域进行赋值，最终将像元组合为4524个区域对象，实现空间相邻且光谱信息相似的像素点合并成为一个对象；④求取各个对象中所有像素的光谱均值，并将该均值作为此对象的光谱数值赋值给所对应的区域对象，实现遥感影像的分割。

步骤3：初始化各类别密度分布待估计的参数Θ，包括各类别的比例、均值向量μ和协方差矩阵∑。根据步骤1中所选取的7个已标记样本作为聚类中心进行最近邻分类，即对数据集里的每个数据点(除了这7个样本点以外)，计算其与这7个样本点的欧氏距离，得到特征空间内的高维距离矩阵，进而计算出所有样本像元隶属于各特征地类的初估概率值。其中，初估概率的核心代码为：

int m_cntClasses；       //聚类数实施例区域

long m_cntData；         //数据数实施例区域

int m_cntAttributes；    //特征数实施例区域

double* m_segmented；    //分割结果实施例区域

double* *m_knowledge；   //知识矩阵实施例区域

double *dis；            //实施例区域内分割结果与知识矩阵间的距离

double**m_probs；        //_probs[i][j]表示第i个样本属于第j类的概率

long i；                 //中间过程参数

intj，k，sumDis，jWatch；//中间过程参数

double sumProbs，max；   //中间过程参数

for(i＝0；i＜m_cntData；i++)

{

     sumDis＝0；

     for(j＝0；j＜m_cntClasses；j++)

     {

        dis[j]＝0；

        for(k＝0；k＜m_cntAttributes；k++)

        {

             dis[j]+＝fabs(m_segmented[i+k*m_cntData]-m_knowledge[j][k])；

         }

         sumDis+＝dis[j]；

     }

     max＝0；

     sumProbs＝0；

     for(j＝0；j＜m_cntClasses；j++)

     {

          m_probs[i][j]＝1-(double)dis[j]/sumDis；

          if(m_probs[i][j]＞max)

          {

             max＝m_probs[i][j]；

             jWatch＝j；

          }

          sumProbs+＝m_probs[i][j]；

     }

     for(j＝0；j＜m_cntClasses；j++)

     {

          m_probs[i][j]/＝sumProbs；

     }

步骤4：遍历数据集，获取当前样本的特征向量及初估概率，对已具有初估概率的对象通过期望步(Expectation Step，简称E步)和极大似然步(Maximum LikelihoodStep，简称M步)进行迭代计算(本实例中设置迭代步长为20)，在依次E步和M步迭代操作完成后，计算出每个样本隶属于各类别的概率。具体描述如下：输入：X＝{x₁，x₂，…，x_N}(存放遥感图像像元样本)，K(遥感影像划分的特征地类个数).

过程：

初始化Z＝{z₁，z₂，…，z_N}，当Z_n∈{1，2，…，K}

估计通过

N_k＝n(z_n＝k)

π_k＝N_k/N

${\mu }_k={\Sigma }_{n=1}^{N_k}{x}_n/N_k$

${\Sigma }_k={\Sigma }_{n=1}^{N_k}(x_n-{\mu }_k){(x_n-{\mu }_k)}^T/N_k$

循环迭代

E步：对于每一个x_n∈X，计算p(z_n＝k|x_n)通过以下公式

$p(z=k|x)=\frac{{\pi }_{k}N\left(x\right|{\mu }_k,{\Sigma }_k)}{\underset{k=1}{\overset{K}{\Sigma }}{\pi }_{k}N\left(x\right|{\mu }_k,{\Sigma }_k)}$

M步：对于k＝1：K，计算通过以下公式

${\pi }_k=\frac{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)}{N}$

${\mu }_k=\frac{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)x_n}{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)}$

${\Sigma }_k=\frac{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)(x_n-{\mu }_k){(x_n-{\mu }_k)}^T}{\underset{n=1}{\overset{N}{\Sigma }}p(z_n=k|x_n)}$

Θ^old←Θ^new

直到ΔlnL(Θ)→0结束循环

1)E步是利用当前初始参数估计结果和遥感像元数据值计算完整样本数据的后验概率，其实质见公式①。首先，遍历数据集，获取当前样本的特征向量及初估概率；计算样本的对数高斯密度函数，并分别对各类特征向量累加对数高斯密度；求得计算样本的分量概率密度，对样本的分量概率密度与先验概率的对数加权相加；对所获得的对数概率转换为正规化的概率。具体的核心代码如下：

//遍历数据集

for(i＝0；i＜m_cntData；i++)

{

    //获取当前样本的特征向量

    for(k＝0；k＜m_cntAttributes；k++)

    {

        x[k]＝m_segmented[i+k*m_cntData]；

    }

//计算样本的Log联合概率密度

int j；

double*logComponentDensities；

double*priors；

logComponentDensities＝new double[m_cntClasses]；

priors＝new double[m_cntClasses]；

for(j＝0；j＜m_cntClasses；j++)

{

     logComponentDensities[j]＝0；

     //遍历特征

     for(k＝0；k＜m_cntAttributes；k++)

     {

        //对数高斯密度累加

        logComponentDensities[j]+＝-log(sqrt(2*3.14159265))

                -log(stdev)-(x-mean)*(x-mean)/(2*stdev*stdev)；

      }

}

//遍历类别

for(j＝0；j＜m_cntClasses；j++)

{

    //先验概率

    priors[j]＝m_priors[j]；//_cnts[j]/(double)m_cntData；

    //加权

    logJointDensities[j]＝log(priors[j])+logComponentDensities[j]；

}

    //将Log概率转换为正规化的概率

    double max＝0；

    double sum＝0；

    for(j＝0；j＜m_cntClasses；j++)

    {

         if(logs[j]＞max)

         {

            max＝logs[j]；

         }

    }

    for(j＝0；j＜m_cntClasses；j++)

    {

         probs[j]＝exp(logs[j]-max)；

         sum+＝probs[j]；

    }

    for(j＝0；j＜m_cntClasses；j++)

    {

         probs[j]/＝sum；

    }

    //保存后验概率

    sure＝-1；

    for(j＝0；j＜m_cntClasses；j++)

    {

        m_probs[i][j]＝probs[j]；

        if(m_probs[i][j]＜0)

        {

            m_probs[i][j]＝0；

         }

    }

}

2)M步是利用当前所有样本的后验概率重新进行参数估计，其实质见公式②。首先，遍历训练集，用各地类的标记样本同模型估计的各聚类中心进行加权叠加(各取二分之一相加)实现特征值的均值化操作，使聚类中心向标记样本漂移，通过公式⑤、⑥、⑦计算出高斯混合模型的参数θ_k＝{π_k，μ_k，∑_k}的值。具体的核心代码如下：

//计算先验概率

for(i＝0；i＜m_cntData；i++)

{

     //遍历类别

     for(j＝0；j＜m_cntClasses；j++)

    {

         //累加概率

         sum[j]+＝m_probs[i][j]；

    }

}

for(j＝0；j＜m_cntClasses；j++)

{

    m_priors[j]＝sum[j]/(double)m_cntData；

}

//计算均值

for(i＝0；i＜m_cntData；i++)

{

    //遍历类别

    for(j＝0；j＜m_cntClasses；j++)

    {

        //遍历特征

        for(k＝0；k＜m_cntAttributes；k++)

        {

            //累加特征值

            m_means[j][k]+＝m_probs[i][j]*m_segmented[i+k*m_cntData]；

         }

    }

}

for(j＝0；j＜m_cntClasses；j++)

{

    for(k＝0；k＜m_cntAttributes；k++)

    {

        m_means[j][k]/＝sum[j]；

        //与先验知识的加权叠加

        m_means[j][k]＝m_means[j][k]*0.5+m_knowledge[j][k]*0.5；

    }

}

//计算方差

for(i＝0；i＜m_cntData；i++)

{

    //遍历类别

    for(j＝0；j＜m_cntClasses；j++)

    {

         //遍历特征

         for(k＝0；k＜m_cntAttributes；k++)

         {

             m_devs[j][k]+＝m_probs[i][j]

          *(m_segmented[i+k*m_cntData]-m_means[j][k])

          *(m_segmented[i+k*m_cntData]-m_means[j][k])；

      }

   }

}

for(j＝0；j＜m_cntClasses；j++)

{

    for(k＝0；k＜m_cntAttributes；k++)

    {

        m_devs[j][k]/＝sum[j]；

    }

}

步骤4：与一般的彩色图像不同，遥感影像的光谱信息中隐含了大量可以利用的光谱特性，如水体的反射主要在蓝绿光波段，其他波段吸收都很强，特别到了近红外波段，吸收就更强，当水中含泥沙时，由于泥沙散射，可见光波段反射率会增加，峰值出现在黄红区，这些知识都可成为遥感影像分类的重要依据。为此，通过光谱特征空间中不同地类在不同指数的知识对其进行修正，通过指数计算实现图像增强，达到识别某些信息或去掉某些不必要信息的目的。本实例中通过计算可见光指数(公式⑧)、归一化植被指数(公式⑨)、归一化水体指数(公式⑩)来的构建新的光谱空间规律。这三个计算指数均反映了地物的亮度信息，VLI通常指向反射率较高的建设用地，NDVI指向植被信息，NDWI指向水体信息。对这三个指数进行拉伸和直方图均衡化操作，分别得出分别反映建设用地、植被和水体的地类概率空间结果图，分别如图3、图4、图5所示。

步骤5：根据比对所计算出的每个样本所隶属于各自类别的概率，判定所属概率最大的类别即为其所属特征地类，进行标记和赋值，利用GDAL的方法对分类结果进行图像输出，结果如图6所示。

为了检测基于半监督学习的遥感影像特征地类识别方法具体的识别效果，将该方法得到的结果同常用的最大似然识别方法及标准EM识别方法得到的结果进行比较，具体采用混淆矩阵的方法对识别结果进行精度评价。混淆矩阵是一个m行、m列的二维矩阵，其中，m代表类别的数量，列通常表示为地面的实值，行表示分类的结果。基于混淆矩阵，可以计算用以评价分类精度的生产者精度、使用者精度和总精度以及Kappa系数。随机选取256个检测点，建立混淆矩阵，计算相关指标进行评价，三种方法相应的精度评价分别见表1、表2和表3。

表1半监督学习分类的精度评价结果

表2最大似然分类的精度评价结果

表3标准EM分类的精度评价结果

精度评价结果表明，半监督学习分类算法有效提高了分类的精度。值得一提的是，最大似然分类算法对训练样本数量的要求较高，而半监督学习分类算法对于每个地类只需要选取一个样本。可见，基于半监督学习的遥感影像分类算法在效果和效率上均优于最大似然分类算法，对标准EM算法的改进效果也较显著，并具备一定的通用性。