[0001] 本发明涉及土地分类技术领域，特别涉及一种土地利用类型分类的方法和系统。

[0002] 随着全球变化研究的深入，土地利用/覆被变化研究已经成为全球环境变化研究的核心领域内容。应用遥感技术结合地理信息系统(Geographic Information System，简称GIS)、计算机技术以及传统调查方法进行土地利用类型的解译和分类已经成为当前获得大尺度、高时效性土地利用时空数据的重要手段。

[0003] 在遥感技术的研究中，通过遥感影像判别各种目标是遥感技术发展的重要一环，遥感数据库的建立关乎专业信息提取、动态变化预测和专题地图制作等重要方面。

[0004] 土地利用遥感分类实际上是遥感图像自动分类识别的过程，也就是用计算机模拟人类知觉，完成遥感图像分析和理解的过程。土地利用遥感分类的核心问题即是一个对遥感图像特征分析提取、图像分割和聚类，进行分类识别的过程。土地利用遥感分类的具体过程是把遥感图像中每一个像元或区域归为土地利用类型分类系统中的一种类别，也就是通过对各类地物的光谱特征分析来选择特征参数，划分出特征空间，将遥感图像的像元划分到特征空间中。在现有技术中，常用的土地利用遥感分类方法有：目视解译法、监督分类法和非监督分类法。

[0005] 目视解译法，主要是通过遥感图像处理软件对遥感图像进行任意放大、缩小及图像增强处理，以达到最佳目视效果，判读人员直接沿影像特征边缘勾绘出地类界线。这种方法的缺点是主要依靠人工解译的方式进行分类，不仅耗费时间长，而且不同人员的解译结果不同，导致分类结果存在差异，无法实现自动分类。

[0006] 监督分类法，又称训练区分类法，基本特点是在分类之前通过实地抽样调查，配合人工目视判读，实现对遥感图像上某些抽样区域中的地物类别属性预先已知，即选取训练区，再使计算机按照已知类型特征去拟出判决函数，从而实现对土地利用类型的分类。此方法的缺陷与目视解译法类似，在此方法中选取训练区的过程同样是靠人工判断，所以工作量大，耗时长，而且结果存在差异，无法实现自动化。

[0007] 非监督分类法，即是仅凭遥感图像的光谱特征的分布规律进行分类，对不同的类型实现区分，是一种自动进行分类的方式。此方法的缺陷在于，按照该方法分类，只能区分出不同类型的土地，无法确定其土地利用类型，而且精确度较低，无法满足实际应用的需要。

[0077] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0078] 本发明方法主要包括：

[0079] 收集待分类地区的历史土地利用矢量数据，根据所述历史土地利用矢量数据建立若干土地利用类型三维特征空间模型，其中，所述历史土地利用矢量数据包括斑块和斑块土地利用类型信息，每个土地利用类型三维特征空间模型与一种土地利用类型对应；

[0080] 获得待分类地区的卫星遥感影像栅格数据，所述卫星遥感影像栅格数据包括若干像元，每个像元是将待分类地区的卫星遥感影像利用栅格划分而得的子单元；

[0081] 将所述历史土地利用矢量数据与卫星遥感影像栅格数据进行叠加处理，获得卫星遥感影像栅格数据中每个像元的土地利用类型信息；

[0082] 依次判断每个像元是否在与其属于同一土地利用类型的土地利用类型三维特征空间模型之内，获得落在土地利用类型三维特征空间模型之外的像元，在卫星遥感影像栅格数据中将所述落在土地利用类型三维特征空间模型之外的像元标记为第一类变化像元；

[0083] 将第一类变化像元当前的土地利用类型转换为其实际所属的土地利用类型，实现土地利用类型分类。

[0084] 参照图1所示，本发明所述方法，具体步骤如下：

[0085] s1、收集待分类地区的历史土地利用矢量数据，所述历史土地利用矢量数据包括斑块和斑块土地利用类型信息，所述斑块土地利用类型信息记录每个斑块对应的土地利用类型；为每种土地利用类型的选取多个样本斑块，并得到样本斑块对应的样本斑块土地利用类型信息，包括样板斑块与样本斑块土地利用类型信息的矢量数据为样本矢量数据；

[0086] s2、对所述样本矢量数据进行训练区分析，得到训练区矢量数据；所述训练区即斑块中纯净的、单一类型的、已知土地利用类型的图像区域；

[0087] s3、通过将训练区矢量数据与卫星遥感影像栅格数据叠加处理，为每一种土地利用类型建立一个土地利用类型三维特征空间模型；所述卫星遥感影像栅格数据是指将待分类地区的卫星遥感影像利用栅格划分为多个子单元而得到的栅格数据，所述子单元被称为像元；

[0088] s4、将所述历史土地利用矢量数据与卫星遥感影像栅格数据叠加处理，叠加后像元落在历史土地利用矢量数据中的一个斑块上，该斑块的土地利用类型信息即为落在该斑块中的所有像元的土地利用类型信息，再判断所有像元是否在其对应类型的土地利用类型三维特征空间模型之内，找到落在土地利用类型三维特征空间模型之外的像元，在卫星遥感影像栅格数据中标记落在土地利用类型三维特征空间模型之外的像元为第一类变化像元；将所述标记第一类变化像元的栅格数据作为第一类栅格数据；

[0089] s5、将第一类变化像元当前的土地利用类型转换为其实际的土地利用类型，转换后在卫星遥感影像栅格数据中标记转换土地利用类型的第一类变化像元为第二类变化像元，将标记过第二类变化像元的卫星遥感影像栅格数据作为第二类栅格数据，由此通过转换第一类变化像元土地利用类型而实现土地利用类型分类。

[0090] 在步骤s1当中，所述选取每种土地利用类型的样本斑块，在具体实施例中具体为，针对一种土地利用类型，将属于该土地利用类型所有斑块按照面积从大到小进行排序，所述斑块中包含矢量数据和土地利用类型信息，矢量数据表示该斑块的形状和面积，根据公式：

[0091]

[0092] 其中Pa表示第一阈值，或称之为选取样本阈值；

[0093] 代表按照面积由大到小排序后的斑块按顺序累加到第i个斑块之和；

[0094] As代表该土地利用类型所有斑块的总面积。

[0095] 实施例中设置第一阈值Pa取值在55％-65％，也就是说，当较大的i个斑块面积之和占所有斑块总面积55％-65％时，则选取这较大的i个斑块作为样本，样本斑块中包含的矢量数据即样本矢量数据；针对每一种土地利用类型重复进行上述操作，得到每一种土地利用类型的样本斑块和样本矢量数据。

[0096] 在步骤s2中，对所述样本矢量数据进行训练区分析，得到训练区矢量数据在具体实施例中为，取样本矢量数据中的样本斑块，所述样本斑块包括斑块形状信息和斑块面积信息；保持样本斑块形状信息不变，将面积缩小，使缩小后的样本斑块与原样本斑块的面积之比达到预先设置的第二阈值范围；取面积缩小后与原矢量数据的面积之比达到预先设置的第二阈值范围的斑块为训练区，训练区的土地利用类型信息与样本斑块土地利用类型信息一致，训练区与训练区土地利用类型信息构成训练区矢量数据。

[0097] 参考公式：

[0098]

[0099] 其中Ad为将样本矢量数据按比例缩小后的面积，A为样本矢量数据的原始面积，当二者之比，即斑块训练区分析值Pbuffer达到预先设定第二阈值的范围，则认为Ad为经过训练区分析后得到的训练区；所述Pbuffer在实施例中可以根据需要认为设定；按照上述方法对每个样本矢量数据进行训练区分析，得到相应的训练区。

[0100] 步骤s3中所述，为每一种土地利用类型建立土地利用类型三维特征空间模型在实施例中的具体方法是，针对一种土地利用类型，将该土地利用类型的训练区与经过标准化处理的多光谱或高光谱卫星遥感影像栅格数据进行叠加处理，得到含有多个波段的样本光谱数据；在本实施例中采用环境1号小卫星遥感影像进行叠加，所述多个波段的样本光谱数据具体包括可见光1、可见光2、可见光3、近红外线共4个波段，表示方式如下矩阵所示：

[0101]

[0102] 此矩阵中n代表本土地利用类型中，环境1号小卫星遥感影像包含像元的数量，矩阵中每一行代表一个像元1、2、3、4四个波段的光谱值；

[0103] 再对上述矩阵进行主成分分析，主成分分析是设法将原来众多具有一定相关性的多个波段的光谱值，重新构造成一组新的互相无关的综合指标值作为样本光谱数据的主成分，各主成分之间不包括重复特征，并且各主成分包含的特征依次减少。

[0104] 本实施例中所述主成分分析的数学方法为正交分解转换，具体为利用四个波段的光谱值来构建四个互相无关联的新变量，所述新变量被称为主成分。所述正交分解转换的方法为目前本领域广泛使用的数学手段，再此不做赘述。

[0105] 计算得到各主成分的方差，所述计算得到的主成分方差数值的大小即为该主成分代表的特征，方差数据越大表示该主成分代表的特征越多，也就是该主成分解释所述样本光谱数据的能力越强。因此选取方差最大的主成分作为第一主成分，如果第一主成分不足以代表原来多个光谱值所包括的特征，则需要再进一步选取多个主成分反映原有特征。

[0106] 得到主成分矩阵记做：

[0107]

[0108] 本实施例中选取包含特征最多的三个主成分，即前三个主成分；

[0109] 在样本光谱数据经过主成分分析之后，可以认为前三个主成分代表了样本光谱数据中大部分的特征，所以取前三个主成分进行建模即降低了维度，即实现了三位特征空间的建立，又在建立的模型中包含了样本光谱数据中大部分的特征；

[0110] 所述土地利用类型三维特征空间模型表达公式如下：

[0111]

[0112] 其中Pj代表代入模型像元的第j主成分的值，MPj代表该土地利用类型中，第j主2
成分的均值， 代表该土地利用类型中第j主成分的标准差，c 为常量，代表特征空间指数，在实施例中按照具体需要自行设定；

[0113] 根据上述方法，对每种土地利用类型建立对应的土地利用类型三维特征空间模型。

[0114] 按照本实施例中步骤s3所述方法，在步骤s4中，将所述历史土地利用矢量数据与GMS遥感影像叠加处理，得到一个普通光谱数据，叠加后像元落在土地利用矢量数据的一个斑块中，该斑块土地利用类型信息即所有落在该斑块中像元土地利用类型信息，按照上述主成分分析的方法再对得到的普通光谱数据进行主成分分析，得到普通光谱数据每个像元的前三个主成分，并逐个像元的将前三个主成分代入该像元对应的土地利用类型三维特征空间模型公式；判断像元是否落在土地利用类型三维特征空间模型之内，找到落在土地利用类型三维特征空间模型之外的像元。

[0115] 对于落在土地利用类型三维特征空间模型之内的像元，记做是不变像元，认为该像元所述土地利用类型为原有类型，不发生改变；对于落在土地利用类型三维特征空间模型之外的像元，在卫星遥感影像栅格数据中标记为第一类变化像元。

[0116] 另外，在此处需要说明的是，按照步骤s3中所述公式所制定的土地利用类型三维特征空间模型形状为一个椭球形空间，所述椭球形空间的大小是通过特征空间指数确定。具体应用中，第一类变化像元本质上即是落在椭球形空间之外的像元。所以为控制第一类变化像元的数量及第一类变化像元在所有像元所占比例，在制定某一种土地利用类型三维特征空间模型的时候，应当按照实际情况设置特征空间指数，致使计算结果在合理范围之内。

[0117] 本实施例中，根据上述过程，则步骤s5具体为：收集步骤s4中得到的所有第一类变化像元，将每个第一类变化像元的三个主成分再次代入到所有类型的土地利用类型三维特征空间模型，计算出该第一类变化像元与每一个土地利用类型三维特征空间模型的距离，并找到与该第一类变化像元距离最近的一个土地利用类型三维特征空间模型；

[0118] 所述计算距离公式如下：

[0119]

[0120] min(di)

[0121] 其中di代表代入的像元与第i土地利用类型三维特征空间模型的距离；Pj代表代入像元第j主成分的值，MPji代表第i土地利用类型中第j主成分的均值， 代表第i土地利用类型中第j主成分的标准差，代表第i土地利用类型三维特征空间模型中的特征空间指数；min(di)代表该第一类变化像元与所有土地利用类型三维特征空间模型的距离中最小的距离。

[0122] 按照上述计算，将所有第一类变化像元的土地利用类型转换为与每个第一类变化像元距离最近的土地利用类型三维特征空间模型对应的土地利用类型，在卫星遥感影像栅格数据中标记转换土地利用类型的像元为第二类变化像元，将标记过第二类变化像元的卫星遥感影像栅格数据作为第二类栅格数据。

[0123] 为了使分类更加精确，本实施例中，在转换第一类变化像元土地利用类型之前，进一步设置了一个转移阻力矩阵，所述转移阻力矩阵中包括了任意两种土地利用类型之间转换需要计算在内的阻力系数，所述阻力系数根据两种土地利用类型发生转换的概率大小而规定；按照以下公式，将步骤s5中得到的第一类变化像元与每一种土地利用类型三维特征空间模型的距离进一步乘以相应的阻力系数：

[0124] dLi＝ρi*di

[0125] min(dLi)

[0126] 其中ρi代表代入像元原始的土地利用类型向第i类土地利用类型转换的阻力系数；

[0127] dLi则代表基于转换规则之下代入像元与第i土地利用类型三维特征空间模型的距离；

[0128] min(dLi)则代表基于转换规则之下，该第一类变化像元与所有土地利用类型三维特征空间模型的距离中最小的距离。

[0129] 按照上述方法，步骤s5中所述将第一类变化像元当前的土地利用类型转换为其实际的土地利用类型在具体实施例中可以具体为，将所有第一类变化像元的土地利用类型转换为基于规则之下与每个第一类变化像元距离最近的土地利用类型三维特征空间模型对应的土地利用类型，在卫星遥感影像栅格数据中标记基于规则转换土地利用类型的像元为第三类变化像元，将所述标记过第三类变化像元的卫星遥感影像栅格数据作为第三类栅格数据。

[0130] 为更进一步的提高本发明所述方法的分类精确程度，所述方法还包括如下步骤：

[0131] 将第一类栅格数据与第二类栅格数据进行叠加分析对比，依据溶蚀算法对不同土地利用类型中面积较小，分布离散的斑块进行修正。

[0132] 或将第一类栅格数据与第三类栅格数据进行叠加分析对比，依据溶蚀算法对不同土地利用类型中面积较小，分布离散的斑块进行修正。

[0133] 对应上述方法，本发明中还公开了一种实现上述方法的土地利用分类系统，参照图2所示。

[0134] 所述系统包括：

[0135] 收集单元、模型建立单元，计算单元和转换单元；

[0136] 收集单元，用于收集待分类地区的历史土地利用矢量数据和卫星遥感影像栅格数据；

[0137] 模型建立单元，连接收集单元，用于根据所述历史土地利用矢量数据建立若干土地利用类型三维特征空间模型；

[0138] 计算单元，连接模型建立单元，用于将所述历史土地利用矢量数据与卫星遥感影像栅格数据叠加处理，获得卫星遥感影像栅格数据中每个像元的土地利用类型信息；依次判断每个像元是否在与其属于同一土地利用类型的土地利用类型三维特征空间模型之内，获得落在土地利用类型三维特征空间模型之外的像元，在卫星遥感影像栅格数据中将所述落在土地利用类型三维特征空间模型之外的像元标记为第一类变化像元；

[0139] 转换单元，连接计算单元，用于将第一类变化像元当前的土地利用类型转换为其实际所属的土地利用类型。

[0140] 所述模型建立单元包括：

[0141] 样本采集单元，训练区分析单元和构建单元；

[0142] 样本采集单元，用于为每种土地利用类型选取多个样本斑块，并得到样本斑块对应的土地利用类型信息，其中，包括样本斑块与样本斑块土地利用类型信息的矢量数据为样本矢量数据；

[0143] 训练区分析单元，用于对所述样本矢量数据进行训练区分析，得到训练区矢量数据；

[0144] 构建单元，用于将训练区矢量数据与卫星遥感影像栅格数据叠加处理，为每一种土地利用类型建立一个土地利用类型三维特征空间模型。

[0145] 所述转换单元包括：

[0146] 距离计算单元和类型信息变更单元；

[0147] 距离计算单元，用于计算第一类变化像元与所有土地利用类型三维特征空间模型的距离，并找到与该第一类变化像元距离最近的土地利用类型三维特征空间模型；

[0148] 类型信息变更单元，用于将第一类变化像元的土地利用类型转换为与该变化像元距离最近的土地利用类型三维特征空间模型对应的土地利用类型。

[0149] 所述转换单元进一步包括：

[0150] 转移阻力矩阵单元，用于设置转移阻力矩阵，所述转移阻力矩阵中包括了任意两种土地利用类型之间转换需要计算在内的阻力系数，所述阻力系数根据两种土地利用类型发生转换的概率大小而规定；转移阻力矩阵单元将第一类变化像元与每一种土地利用类型三维特征空间模型的距离进一步乘以相应的阻力系数，得到基于转换规则之下第一类变化像元与每一种土地利用类型三维特征空间模型的距离，获得基于转换规则之下与第一类变化像元距离最近的土地利用类型三维特征空间模型；

[0151] 所述类型信息变更单元将第一类变化像元的土地利用类型转换为基于规则之下与该变化像元距离最近的土地利用类型三维特征空间模型对应的土地利用类型。

[0152] 所述系统进一步包括：

[0153] 所述转换单元还用于将第一类变化像元当前的土地利用类型转换为其实际的土地利用类型后，在卫星遥感影像栅格数据中将经过土地利用类型转换的第一类变化像元标记为第二类变化像元；

[0154] 所述系统还包括溶蚀算法单元，连接计算单元与转换单元，用于提取计算单元标记过第一类变化像元的卫星遥感影像栅格数据作为第一类栅格数据；提取标记过第二类变化像元的卫星遥感影像栅格数据作为第二类栅格数据；对比第一类栅格数据和第二类栅格数据，依据溶蚀算法对转换土地利用类型之后分布离散的像元进行分类信息的校正。

[0155] 以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。