CN101308544B - 一种基于栅格的空间异质模式识别方法 - Google Patents

Abstract

一种基于栅格的空间异质模式识别方法及分层方法，其识别方法步骤如下：(1)特征选择：依据目标变量确定影响因子变量集；(2)对影响因子变量集经过预处理，形成统一分辨率及坐标的多维栅格数据集；(3)对形成的统一分辨率及坐标的多维栅格数据集，采用空间共享最近邻方法，计算栅格单元间相似性，从而进行非监督分类，得到聚类结果；(4)对步骤(3)得到的聚类结果建立规则库，采用规则库中的规则运算于新的数据集，产生空间异质模式图，再将此方法所产生的空间异质模式，可以用于抽样调查中的样本分层方法。本发明对我国大范围的资源调查及区划、遥感精度评估及环境监测具有重要理论意义及实用价值。

Description

一种基于栅格的空间异质模式识别方法

技术领域

本发明一种基于栅格的空间异质模式识别方法，及利用该方法所产生的空间异质模式，用于抽样调查中的样本分层方法，属于地球信息科学技术领域。 

背景技术

空间异质(spatial heterogeneity)模式是指空间或者区域差异的特性，这种差异性主要是地区间各种不同的诸如自然地理条件、社会经济及专题属性等因子影响而产生的(Anselin，1992)。区域(空间)差异性在地域广阔、气候及经济发展水平跨度大的我国尤为明显。正确识别区域(空间)异质模式具有重要的意义，可以帮助人们把握区域之间的实际差异制定不同的调查方案或区域政策，直接的用途就是可以根据这些异质模式在资源环境调查中进行样本分层，提高估计的精度及调查的效率。 

最简单的空间异质模式是行政区划，从政治意义上划分行政边界，没有涉及到具体的区域条件；更复杂一些的则是根据一定的知识或者机理划分区域，如生态分区、气候分区，这些划分虽然相对客观，但划分过于笼统，缺乏区域的实际情况既在一定证据(数据)的支持，因此划分有可能存在较大的偏差。 

比较客观的空间异质模式识别及划分是采用机理同数据相结合的方法，目前主要是基于数据分析的单要素或者多要素方法。单要素是根据目标变量进行分级划分，方法简单可行，但主要的问题是一般很难直接获取即时覆盖区域的目标要素信息，而且有时候区域之间的差异是综合多个因素的变量(如东部与西部之间的发展差异)，因此单要素方法有很大的局限性。相对而言，多要素方法可以获取很多相关信息(如区域的自然地理属性与社会经济变量)，合成这些变量对区域在某一方面的差异性进行可行与客观的估计(尤其是涉及到综合性指标的区域差异性，如东西部发展差距)。 

区域差异的一个直接应用领域是区域抽样调查中的分层，可以直接将空间(区域)的异质模式作为分层框架，有助于提高总体变量估计值的精度。当前国内外的许多抽样分层方案(如采样)分为两类：一种直接采用已有分区单元对样本划分层次，如行政单元、生态分区及土地利用类型等；另外一种采用一定的分类方法根据单要素或其他变量进行分层，如McRoberts等报告了许多采用遥感图像进行分层的实验研究。 

无论是直接的行政区划、单要素方法还是多要素方法，现存的分区(分层)方法有一定的缺点： 

(1)单要素分层方法(Li et al.，2005；Wang et al.2002；刘兆德等，2001；Cochran，1977；冯士雍等，1998)过于简化，如采用行政单元、高程或土壤类型，对于估计目标未必适用(如估计目标为耕地面积比例)，有时候很难获取覆盖全区域的单要素因子进行分层，而采用插值算法得到覆盖全区的结果存在很大的误差。 

(2)当前多要素分析方法主要是基于非监督分类方法(McRoberts et al.，2006；Bergen et al.，2005；Gallego，2005；Tan et al.，2006)，包括分级判别树、ISODATE、K-Means等，这些方法简单易行，但主要的问题是对于涉及到包括遥感数据在内的高维空间数据集学习结果并不理想，尤其是这些算法都以传统的距离函数(如欧氏距离、夹角余弦等)衡量数据点之间的相似性，不能适应数据项中因属性值变化而引起的大小、形状、密度的变化，针对标准数据集的分类效果不理想；算法也没有融入空间相关性的信息，不能反映地球科学数据中“越近越相关”的Tobler(1979)第一定律。 

(3)当数据量巨大时，如果算法很复杂，将耗费大量的时间进行计算、空间搜索及比较等，算法运行耗时，不适合处理高维的海量数据。 

(4)现有的非监督分类算法(McRoberts et al.，2006；Bergen et al.，2005；田光进等，2001；Gallego，2005)用于识别空间异质模式或者分层框架时，使用的影响因子数量有限，主要局限于连续型变量，其中很多直接采用遥感传感器的光谱数据进行分类达到分层的目标，没有综合性地考虑影响区域的其他因子(如类别变量，包括地质、地貌、植被类型)，而这些因子可能对区域差异影响比较重要。 

发明内容

本发明的技术解决问题：克服现有技术的不足，提出一种基于栅格的识别空间异质模式的方法，该方法提高了识别空间异质模式的分类效果，此外还提高空间异质模式识别的计算效率。 

本发明还提供一种采用上述基于栅格的识别空间异质模式的方法所产生的空间异质模式图，用于抽样调查中的样本分层方法，该分层方法可为区域的资源环境调查提供信息量丰富的分层框架，以提高估计精度或分层的效率。 

本发明基于栅格的空间异质模式识别方法步骤如下： 

(1)特征选择，既依据目标变量y确定影响因子集X，方法为： 

a.主要采用四种方法来衡量采集的备选影响因子变量x和目标变量y之间的相关性； 

b.如果两者相关，则保留此备选影响因子变量x，多个备选影响因子变量x形成影响因子变量集X。 

上述判断备选影响因子变量x和目标变量y之间的相关性的方法为： 

a.同时采用条件直方图、二维散点图、基于信息熵的分类贡献函数和粗糙集的属性重要性函数四种相关性分析方法来判断备选影响因子变量x和目标变量y之间的相关性；当备选影响因子变量x为连续变量时，采用二维散点图来判断影响因子变量x和目标变量y之间的相关性；当备选影响因子变量x为类别变量或者离散变量，或者将连续变量转化成离散变量时，采用条件直方图来判断影响因子变量x和目标变量y之间的相关性，如果条件直方图的差异显著，则说明两者之间具有相关性，保留此备选影响因子变量x；对于连续变量，或类别变量，或离散变量均可以采用基于信息熵的分类贡献函数，即IG及GR，或粗糙集的属性重要性函数来判断影响因子变量x和目标变量y之间的相关性；所述的基于信息熵的分类贡献函数，即IG及GR的计算公式为： 

IG(y，x)＝H(y)-H(y|x) 

GR(y，x)＝(H(x)-H(y|x))/H(x) 

其中，H(y)是熵的定义，H(y|x)是对类别或离散化后以影响因子变量x为条件的y的条件熵定义； 

所述的粗糙集的属性重要性函数的计算公式如下： 

${\mathrm{\σ}}_{(x,y)}\left(x_i\right)=1-\frac{\mathrm{\γ}(X-x_i,y)}{\mathrm{\γ}(X,y)}$

其中，X为影响因子变量集，x_i为计算属性重要性的影响因子变量，γ(X，y)是粗糙集中目标变量y对影响因子变量集X的依赖性； 

b.IG和GR或σ_(X，y)(x_i)值越大，表明两者越相关，保留此备选影响因子变量，如果IG和GR或σ_(X，y)(x_i)值接近0，则表明两者不相关，移除此备选影响因子变量； 

c.对同时采用上述四种方法的所选择的影响因子变量x采用变差膨胀因子法进行多重共线性分析，以消除冗余信息，提高影响因子变量集X的信息含量。 

令变差膨胀因子为VIF，则有：当VIF＞10时，计算每个影响因子变量x的VIF值，去除具有最大VIF的影响因子变量x，再进行判断计算，直到所有的影响因子变量x的VIF＜10，所得的影响因子变量集即作为的影响变量因子集X。 

该步骤(2)中融合多种类型及多来源的空间数据进行空间异质模式的识别，每个因子变量对应了空间数据库的一个图层，格式可以是栅格或者矢量，数据类型可以是类别、离散或者连续，数据可根据目标变量从来自于气候气象(如降雨、温度)、自然地理条件(如高程、坡度、地貌、植被)、专题属性(如土地利用)、社会经济(如人口 分布、城镇分布及GDP)及遥感影像(如TM影像的特定波段、NDVI数据)等空间数据库中进行选择；而影响因子变量集的选择采用了四种相关性分析方法及多重共线形方法获取冗余信息含量少、信息量高的针对目标变量y的影响因子集X。 

(2)对每个选中的影响因子变量x经过预处理，形成统一分辨率及坐标的多维栅格数据集； 

不同的存储方式的影响因子变量x的预处理方法不同： 

a.对于矢量数据类型的预处理 

当影响因子变量x以矢量格式存储时，将矢量转变成栅格，形成统一的栅格数据；对于点矢量数据类型，采用k-最近邻方法或者Kriging方法插值算法，将点数据转换成栅格数据；对于线矢量数据类型采用数字微分分析法DDA法或Bresenham方法转换成栅格数据；对于面矢量数据类型采用扫描法转换成统一分辨率的栅格数据； 

b.对于栅格数据的预处理方法 

当影响因子变量x以与目标变量y的空间分辨率或坐标不一致的栅格格式存储时，采用最近邻法重采样栅格，形成统一分辨率的栅格数据； 

c.对步骤(1)和步骤(2)形成统一的栅格数据，如果为离散或类别数据，则将其数值化，即将这些数据转换成数值型，便于同其他的连续型数据共同聚类分析，转换方法采用n-1变量法。所述的采用n-1变量法是将有n个类的类别或离散数据划分成具有0-1二值的n-1个变量。 

对于离散或类别数据的数值化，即将这些数据转换成数值(连续)型，便于同其他的连续型数据共同处理与分析，采用n-1变量法将有n个类的类别或离散数据划分成具有0-1二值的n-1个变量的方法。 

(3)对形成的统一分辨率及坐标的多维栅格数据集，采用空间共享最近邻方法，计算栅格单元间相似性，从而进行非监督分类，得到聚类结果； 

其中空间共享最近邻算法Spatial SNN如下： 

步骤1计算栅格数据中每两个栅格单元之间影响因子变量集X的内积、或余玄函数或扩展Jaccard，计算的具体方式为：令任意两个栅格单元的影响因子矢量为x₁与x₂，则选择以下三种方法中一种计算两个栅格单元之间的基本的相似性系数(简称为基相似性系数)： 

a.采用内积的计算式为： 

s^(P)(x₁，x₂)＝(x₁)^Tx₂(T表示因子矢量x₂的转值) 

b.采用余弦函数的计算式为： 

s^(C)(x₁，x₂)＝((x₁)^Tx₂)/(|X₁|·|x₂|) 

c.采用扩展的Jaccard的计算式为： 

s^(J)(x₁，x₂)＝((x₁)^Tx₂)/(|x₁|²+|x₂|²-(x₁)^Tx₂) 

在实际计算中推荐使用后两种，即余玄函数或扩展Jaccard会取得更好的效果。 

步骤2根据步骤(1)中的计算结果构造相似性矩阵(相似性图)，具体的构造方法是：令栅格的单元数为n，则有以下步骤： 

a.先构造一个空相似性图，其中节点表示每个栅格单元，节点之间的连接权重表示节点之间的基相似性； 

b.根据步骤1计算的基相似性给相应节点之间的连接权重赋值。 

步骤3根据步骤(2)构造的邻近矩阵，相似矩阵的稀疏化及构造共享最近邻图SNN，有以下具体步骤： 

a.用户根据目标及数据特点设定强相似性域值k，仅当x₁与x₂之间的基相似性s(x₁，x₂)≥k时，连接得以保留，否则无效予以删除，从而稀疏化相似矩阵； 

b.根据稀疏化后的矩阵构造最近邻图，即当两个栅格单元x₁与x₂之间彼此共享k个共同的最近邻节点时，二者之间构成共享最近邻居边界，否则二者之间的边界在共享最近邻图中移除。边界的共享最近邻的连接权重值即为共享的最近邻数。 

步骤4：根据前述得到的共享最近邻图SNN计算每点SpSNN_D，SpSNN_D为该点的融合了空间信息的共享最近邻密度函数，依据不同的情况，在一般的共享最近邻密度的基础上本发明融入了空间相关性信息，设计了以下的空间最近邻计算方法，分两种情况采用不同算式。 

令x₁，x₂为特定的两个栅格单元的影响因子矢量，NN(x₁)为由影响因子矢量经相似性计算的x₁的最近邻单元集合，size(A)为集合A中的元素个数，则该特征包括了以下所改进的相似性计算方法，即将空间相关性作为特殊的属性信息融入到相似性计算中，分两种不同的情况采用不同处理方法： 

a.若没有区域的目标变量的空间变异模式信息，或者无法通过相关数据采用地统计学方法求得目标变量的空间变异规律，可直接通过空间连接性将空间相关性融入到相似性计算中： 

$\mathrm{SpSim}(x_1,x_2)=\mathrm{size}(\mathrm{NN}\left(x_1\right)\∩\mathrm{NN}\left(x_2\right))+w_{x_1{x}_2}$

式中， 

为x₁与x₂的空间连接性：x₁与x₂共享一条边为1，只共享一个点为0.5，否则为不相邻(值为0)。对非空间属性，距离的计算可采用内积及余弦函数，可消除欧 氏距离带来的不利影响。 

b.如果掌握了研究区域的目标变量的空间变异模式信息，即得到空间变异模型γ及其相关的参数，则可将空间相关性按照如下方法融入到相似性计算中： 

$\mathrm{SpSim}(x_1,x_2)=\mathrm{size}(\mathrm{NN}\left(x_1\right)\∩\mathrm{NN}\left(x_2\right))+\tilde{C}(x_1,x_2)$

式中， 

体现了x₁与x₂之间的空间关联性，由Kriging方法的空间变异函数γ所确定， $\tilde{C}(x_1,x_2)={\widetilde{\mathrm{\σ}}}^2-\mathrm{\γ}(x_1,x_2).$

由以上(1)或(2)的相似性的计算可得到融合空间信息的空间共享最近邻密度函数： 

SpSNN_D＝|Q|，Q＝{x₂|SpSim(x₁，x₂)≥Eps} 

式中，Eps为用户设定的相似性域值，大于Eps为表明x₁与x₂相似，否则为不相似。 

步骤5根据栅格点的空间共享最近邻密度SpSNN_D识别核心点和边界点，具体包括以下两步： 

a.当一栅格点x的SpSNN_D(x)≥MinPts，则划分为核心点集，其中MinPts根据应用目标及数据特点确定的核心点划分域值； 

b.当一栅格点x的SpSNN_D(x)＜MinPts，但x属于另外一个核心点的共享最近邻，则将其划分为边界点； 

步骤6聚集核心点，把核心点分成类，按照以下步骤进行： 

a.扫描由步骤5所得到的核心点集，如果两个核心点属于彼此的共享最近邻点，则将其归为相同的类别； 

b.通过a步的彼此共享的相互连通的多个核心点构成一个类别，而不连通的不同的核心点集归为不同的类别； 

c.对形成的不同的类别进行标记，如标为1，2，…，n类，n为不连通的共享核心点集的个数。 

步骤7对非核心点，如果为边界点，则对边界点进行聚集，对于糙点，则移除，得到聚类结果。扫描非核心点集，分以下两种情况处理： 

a.如果该点不为核心点，但为边界点，则根据步骤5中边界点的定义，将其划分为其所归属的核心点所在的类别。 

b.如果该点不为边界点，则说明其为糙点，应该移除，不参与聚类。 

(4)对聚类结果建立规则库，采用规则库中的规则运算于新的数据集，产生空间异质模式图。 

具体分为两步： 

步骤1采用决策树与粗糙集建立规则库，主要有以下步骤： 

a.对步骤(3)得出的聚类类别，将每个栅格单元对应的影响因子，聚类类别作为一条训练样本，其中，影响因子作为归纳的各条件参数，而聚类类别作为目标的分类值。这些样本单元构成了用于训练用的样本数据集； 

b.采用10次10倍的交叉验证方法训练样本，采用决策回归树方法建立规则。决策回归树也叫C4.5算法，它通过两个步骤来建立决策树规则集：树的生成阶段和树的剪枝阶段，这种方法的优点是规则归纳、非参数分类、可容纳多种类型数据(类别或连续型)一起进行归纳。决策树归纳法的缺点是过于拟合、对噪音及边界点敏感及随机性选择训练样本影响，可通过元学习算法bootstrap、boosting及二者的混合加强这种算法。因为这是通用的算法，在此不在详细叙述该算法。 

c.采用10次10倍的交叉验证方法训练样本，采用粗糙集方法生成规则。粗糙集是一种分类规则归纳工具，可以处理数据中的不确定性，在许多领域都取得了广泛的应用。粗糙集归纳算法步骤：离散化、相似性构造、属性约简及规则合成。离散化既根据条件变量同决策变量之间的关系离散化连续变量或者对类别变量重新分组；相似性构造既评估属性之间的依赖性及属性重要性便于特征变量的选择；约简既除去不重要的属性及一些属性的冗余值便于好的规则的生成；而规则的合成既从约简表中合成最小长度的决策规则。因为粗糙集也是通用的算法，在此不在详细叙述该算法。 

无论是决策回归树，还是粗糙集，都是非参数监督分类方。从该类方法得到的分类规则以谓词逻辑及决策树作为基本的推理机制，每条规则有前件与后件，彼此通过连接符(即：析取∨，合取∧，取否 )进行连接。前件描述了得出后件分类值的需要满足的条件。对一条规则的不确定性测量为支持度(Support)与置信度(Confidence)：前者指匹配规则前件或后件的对象在总对象中的百分数，而置信度指匹配整个规则的对象在满足规则前件的对象中的百分数。根据经验设置适当的置信度与支持度域值以便于产生高效的规则。 

d.对决策树及粗糙集生成针对每一个聚类类别的多条规则进行合并产生规则集，对于规则的置信度小于80％予以删除，保留可信度高的精简的规则集；此外，产生的规则可结合实地调查及领域专家经验进行检验及提炼(既比较、分析并确定最佳规则)，最后输入到库中作为对新数据分类的规则库。 

步骤2当从步骤1中得到规则库后，可应用这些规则对类似目标的新数据集进行聚类运算，完成识别空间异质模式。具体有以下两步： 

a.首先将目标数据集经过前面步骤(2)的预处理，形成同训练样本相同分辨率的新数据集； 

b.从新数据集中提取相关的因子变量，将这些因子变量作为规则集中的条件变量(规则集主要是IF-THEN的逻辑运算)，采用谓词逻辑运算器、或决策树运算器对这些条件因子变量进行运算，对其所构成的栅格单元进行分类，如对同一个栅格单元分类规则产生不同的类别，则采用“投票”原则，选取最大投票数目的类别作为结果。 

c.得到的各个栅格单元的类别即为所识别的空间异质类，不同类别的栅格单元在空间上的分布构成了空间异质模式图。该栅格单元经过ArcGIS等GIS软件的矢量化后形成了空间异质模式的多边形图。 

(5)对产生的空间异质模式可使用到抽样调查的分层中提高总体比例估计精度及分层效率。 

其步骤如下： 

步骤1得到样本的分层框架：采用ArcGIS等GIS软件中的矢量化功能，矢量化所得到的空间异质模式分类的栅格图像，产生不同空间异质特性的矢量多边形图，作为样本分层的框架。因为矢量化操作是许多GIS软件的功能，此处不在详细叙述该栅格矢量化过程，此处注意的是在矢量化过程中相同的类别码属于同一层，不同的代码形成了不同的层。 

步骤2根据样本分层的框架对样本进行分层：如果样本为多边形或线，当样本被不同层的多边形所覆盖时，需要将子样本沿着层边界划分成几个子样本，每个子样本归于其所在的层，具体包括以下两种不同的处理方法： 

a.如果样本为线，则沿着该线段与不同层交叉的边界处划分成两个或多个样本，如果线段被两个不同类别的层划分，则沿着切点划分到两个不同的层中；如果为3个或以上的层划分，则沿着这些层与该线段交叉的切点划分到不同层中，划分后的线段位于那层就属于该层。 

b.如果样本为多边型，则沿着该多边型与不同层交叉的边界线划分成两个或多个样本，如果该多边型被两个不同类别的层划分，则沿着切点划分到两个不同的层中；如果为3个或以上的层划分，则沿着这些层与该多边型交叉的切线划分到不同层中，划分后的多边型位于那层就属于该层。 

步骤3根据分层结果，对研究区域的总体比例参数进行估计：首先在每一层中采用Ripley的模型公式估算此层的均值和方差，然后再采用Cochran模型公式估计研究区域总体的比例的均值及方差；在估计时需要考虑样本规模对结果的影响，每个子样本的规 模不一样，若样本为多边形，则采用其面积作为其规模因子；若样本为线状，则采用其长度作为规模因子。主要包括两步： 

a.根据每层内的样本估计待估计比例的均值及方差： 

均值： ${\widehat{\mathrm{\β}}}_h\left(a\right)={\∫}_{N_h}{\mathrm{\β}}_h\left(a\right)\frac{S_h\left(a\right)}{{\∫}_{N_h}{S}_h\left(a\right)}\mathrm{da}={\∫}_{N_h}{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)\mathrm{da}=\frac{1}{N_h}{\∫}_{N_a}{n}_h{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)\mathrm{da}$

其中，β_h(a)＝β_ah，既第h层第a个样本单元大小(如面积或线段长度)比例；同样w_h(a)为既第h层第a个样本单元权重。 

均值的估计方差： 

${\widehat{\mathrm{\σ}}}_{{\widehat{\mathrm{\β}}}_h\left(a\right)}{\left(n_h\right)}^2=E_h{\left[{\widehat{\mathrm{\β}}}_h\left(a\right){\mathrm{\β}}_h\left(a\right)\right]}^2=E{[\frac{1}{n_h}\underset{a=1}{\overset{n_h}{\mathrm{\Σ}}}{n}_h{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)-\frac{1}{N_h}{\∫}_{N_h}\left(n_h{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)\mathrm{da}\right)]}^2$

$=\frac{1}{n_h}\{1-E_h[r(a-a^{\′})\left]\right\}{\widehat{\mathrm{\σ}}}_{{\mathrm{\β}}_h\left(N_h\right)}^2=F\left(n_h\right){\widehat{\mathrm{\σ}}}_{{\widehat{\mathrm{\β}}}_h\left(a\right)}\left(N_h\right)$

其中F(n_h)＝(1/n_h){1-E_p[r(a-a′)]}，E_p[r(a-a′)]为空间相关性结构，求解方法由Ripley(1981)及Rodriguez-Iturbe(1974)等提出。 

b.根据第一步所计算出的各层的估计均值及方差，估计总体比例均值及方差，采用Cochran提出的经典方法： 

总体比例参数β_STR的估计为： 

${\widehat{\stackrel{\‾}{\mathrm{\β}}}}_{\mathrm{STR}}=\frac{\underset{h}{\overset{L}{\mathrm{\Σ}}}{n}_h{\widehat{\stackrel{\‾}{\mathrm{\β}}}}_h}{n}=\underset{h=1}{\overset{L}{\mathrm{\Σ}}}{w}_h{\widehat{\stackrel{\‾}{\mathrm{\β}}}}_h$ 其中n＝n₁+n₂+…n_L； 

β_STR均值的估计方差为： $\hat{V}\left(\widehat{\stackrel{\‾}{\mathrm{\β}}}\right)=\underset{h=1}{\overset{L}{\mathrm{\Σ}}}{w}_h^2{\hat{V}}_h\left(\widehat{\stackrel{\‾}{\mathrm{\β}}}\right)$

其中， 

为第h层的方差。 

本发明与现有技术相比有如下优点： 

(1)本发明提出了空间共享最近邻(Spatial SNN)的概念来衡量数据点的相似性。共享最近邻以数据的共享最近邻项作为衡量数据点是否相似的指标，有效地克服常用的直接以最近邻作为衡量指标的不适合处理密度、形状及大小易变的高维数据的局限；而本发明所提出的空间共享最近邻的概念在共享最近邻的基础上融入了空间相关性因子，使得这种方法充分融合了地学数据中客观存在的空间相关性信息，使得不仅能处理密度、形状及大小易变的高维数据，也适合处理空间数据，高效地识别高维数据中的相似 性，能适应高维数据的形状及密度易变的特点，对提高识别空间异质模式的分类效果有重要的贡献。这个优点是常用的非监督聚类方法所没有的。此外，本发明的Spatial SNN对于糙点及边界点不敏感，使运算具有较强的鲁棒性。 

(2)对得到的聚类结果建立规则库，采用规则库中的规则运算于新的数据集，产生空间异质模式图，有效地提高了计算效率。 

(3)此外，本发明将离散或者类别因子变量通过转换同其他的连续型因子变量参与聚类，克服了现有方法因子变量有限与不能包括类别或离散变量进行分类的局限。 

(4)同时本发明通过四种相关性分析方法与多重共线形方法提取与目标变量相关的影响因子变量集，确保影响因子变量信息含量高而冗余信息少，进一步提高分类效果起重要作用。 

(5)针对非监督算法计算耗时的特点，本发明采用了决策树及粗糙集等非参数方法归纳所得分类的规则，这些规则通过简单的IF-THEN、谓词逻辑或者决策树推理算法用高效地用来分析其他的类似目标及背景的海量高维数据，极大地提高了非监督分类既空间异质模式识别的计算效率，也克服了直接使用复杂的聚类识别算法的效率低的缺点。 

(5)当前主要采用已有的区划图(如行政或生态)或者遥感分类图进行样本的分层，很少用多要素聚类结果进行分层，极大地限制了调查的效率或估计精度的提高。本发明将所产生的空间异质模式图用于分层极大地提高了总体估计的精度或调查的效率。 

(6)本发明提出的空间异质模式识别方法对我国大范围的资源调查及区划、遥感精度评估及环境监测具有重要理论意义及实用价值 

附图说明

图1为本发明主要流程图，包括了空间异质模式的识别及分层估计步骤； 

图2为本发明的空间共享最近邻(Spatial SNN)实现过程； 

图3为一个典型案例的三种不同非监督聚类方法结果对比，其中图3a为采用常用的k-means得到的分区，图3b为自组织映射聚类(SOM)方法得到的分区，而图3c为采用本发明得到的分区结果，图中可见，由本发明得到的分区细节丰富，有更多的信息含量； 

图4a和图4b为采用图3得到的空间异质模式图(多边形)对样本进行分层估计，对得到的估计误差及置信区间同其他方法的比较(即本发明提出的方法同k-means聚类、SOM聚类、土壤类型及地貌图比较)。 

具体实施方式

如图1所示，以耕地调查为案例，本发明的实施方式如下： 

1、特征选择，既依据目标变量y确定影响因子集X。 

此处主要以耕地面积比例作为目标变量，研究区域为我国的山东省中部地区，备选的地理空间因子覆盖山东全省，有四类因子：遥感光谱及派生数据，包括Landsat的ETM+波段数据(2-5波段)、NDVI、高程及坡度；自然地理数据，包括土壤类型、植被类型、地形；气象数据，包括年平均表面温度及降雨；社会经济数据，包括GDP与人口密度。总共包括了15个备选因子集。这15个因子通过四种相关性分析及多重共线形分析后，剔除其中的坡度及年平均表面温度两个因子变量，最后得到13个影响因子集作为非监督学习的数据集，即ETM+影像波段数据(2-5共四个波段)、NDVI、高程、坡度、土壤类型、植被类型、地表类型、年平均降雨、GDP及人口密度。 

2、预处理得到统一分辨率的多维栅格数据集。 

预处理涉及到几方面，首先是将矢量转换成栅格数据，对土壤类型及植被类型等多边形数据采用扫描法栅格化，而对于ETM+及其他栅格数据采用最近邻重采样法转换成统一分辨率的栅格数据(分辨率为0.456km×0.456km)集中；然后是对类别或离散数据转化成连续型变量，以便于同连续型数据进入下步的计算(如土壤类型共14类，衍生出13个变量；而地貌类型共11类，衍生出10个变量)。 

3、采用空间共享最近邻方法对栅格单元进行聚类。 

采用空间共享最近邻公式分析识别栅格单元中的核心点，进行栅格单元的分类。聚类的过程首先是根据距离函数构造相似矩阵，然后是相似矩阵稀疏化、构造共享最近邻图、计算Spatial SNN密度、识别核心及边界数据点、聚类及糙点的移除。以下结合图2描述该算法主要步骤(x及x_i为像素单元)： 

(1)采用适当的距离函数构造最近邻矩阵(图2b，其中，距离函数的衡量标准以余玄函数及扩展Jaccard来衡量高维的时空数据的相似性效果较好)，我们此处采用了扩展Jaccard距离测量函数。 

(2)相似矩阵的稀疏化及共享最近邻图的构造(图2c，其中邻近点大小参数k被设置便于发现邻近的节点便于聚类)。 

(3)每点的SpSNN密度的计算(图5d，其中SpSNN_D(x)为x的融合了空间信息的共享最近邻密度函数，Eps为相似性域值用于判断两点是否相似。而SpSim(x₁，x₂)则为融合了空间信息的点x₁与x₂之间的相似性函数，既空间连接性相似性公式或空间变异性公式，NN(x₁)为x₁的邻近点数目)。 

(4)核心与边界点的识别(图2e，其中CS与BS分别为核心点集与边界点集，Nb(x₁)为根据共享最近邻密度标准确定的x₁点的邻近点集，而MinPts为核心点的域值，由用 户设定)。 

(5)核心点的聚集(图2f，其中Cl(x₁)为x₁对应的栅格点的聚类标记)。 

(6)边界点的聚集及糙点的识别与移除(图2g其中NS为糙点集)。 

在以上的学习算法中，各参数值为：k＝50，Eps＝30，MinPts＝20。 

4、根据聚类分析的结果采用决策树及粗糙集归纳分类规则，建立空间异构模式识别的规则库。聚类划分成6个不同的类别，采用监督归纳算法得到总共800条规则，其中，类别1共有219条规则，类别2共有176条规则，类别3有174条规则，类别4有210条规则，类别5有49条规则，类别5有5条规则。一条典型的规则如： 

$\→\mathrm{Class}(x_{i,j},3)(\mathrm{Support}=80\%,\mathrm{Confidence}=90\%)$

本条规则可解译为如果一单元(栅格单元)，其属性(因子变量)满足一定的条件，则其聚类类别为3类，属于中等水平。本条规则的支持度为80％，而信任度为90％。这些产生式知识规则通过框架编码方式保存到知识库中，使用的时候从知识库中调用，通过谓词推理算法运用这些规则从而实现对海量多维空间数据的聚类。 

5、采用k-means、自组织映射(SOM)聚类方法对相同的变量因子进行相似性计算识别空间异质模式，计算这些方法及空间共享最近邻规则方法所得结果的β指数。 

计算得到三种聚类方法的β指数，结果表明本发明所提出的方法的β指数最高，为2.57；而k-means的为1.32，SOM方法为2.17。β指数说明我们的方法效果是最佳的(相同条件下最大)。实验也考虑了变差因素的影响，在山东省内选择不同地块进行了10次实验，所提出的改进方法β平均值为2.478，变差0.2；而k-Means均值为1.201，变差为0.51；SOM均值为2.101，变差为0.32。消除了随机因素的影响，说明本发明可取得较好计算效果。 

图3显示了本案例的三种不同聚类方法结果对比，其中a图为采用k-means得到的分区，b图为SOM得到的分区，而c图为采用本发明提出的空间共享最近邻算法得到的规则经过推理运算得的分区结果。图中可见，由本发明方法得到的分区细节丰富，信息量丰富。 

6、以上几种不同方法得到的分层图对样本数据进行分层抽样及估计，根据结果及误差评估方法的效果及精度。用于分层及估计的样本为长方形的航片样本，每个航片样本的大小为28公里×23公里，而用于分层的栅格数据分辨率为0.456公里×0.456公里，满足了误差要求。航片的属性为每航片样本中的耕地面积比例，航片样本(包括被 划分成子样本)的面积作为样本的规模。估计目标为研究区域的耕地面积比率的总体均值。 

分层与估计的步骤： 

(1)将不同方法(既k-means、SOM及空间共享最近邻图)得到的空间异质模式栅格图转化成矢量图，采用了ARCGIS中的栅格矢量化功能，转化成Shape文件格式； 

(2)将不同方法所得到的分层框架与航片样本进行叠加分析(既Overlay)，并将被不同层覆盖的样本划分成多个子样本，每个子样本属于其所在层。 

(3)采用估计了样本规模既航片样本的面积的模型公式求解总体参数与变差，该步分两小步： 

第一步，根据每层内的样本估计耕地面积比例的均值及方差： 

均值： ${\widehat{\mathrm{\β}}}_h\left(a\right)={\∫}_{N_h}{\mathrm{\β}}_h\left(a\right)\frac{S_h\left(a\right)}{{\∫}_{N_h}{S}_h\left(a\right)}\mathrm{da}={\∫}_{N_h}{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)\mathrm{da}=\frac{1}{N_h}{\∫}_{N_a}{n}_h{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)\mathrm{da}$

其中，β_h(a)＝β_ah，既第h层第a个样本单元面积比例；同样w_h(a)为既第h层第a个样本单元权重。 

均值的估计方差： 

${\widehat{\mathrm{\σ}}}_{{\widehat{\mathrm{\β}}}_h\left(a\right)}{\left(n_h\right)}^2=E_h{\left[{\widehat{\mathrm{\β}}}_h\left(a\right){\mathrm{\β}}_h\left(a\right)\right]}^2=E{[\frac{1}{n_h}\underset{a=1}{\overset{n_h}{\mathrm{\Σ}}}{n}_h{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)-\frac{1}{N_h}{\∫}_{N_h}\left(n_h{\mathrm{\β}}_h\left(a\right)w_h\left(a\right)\mathrm{da}\right)]}^2$

$=\frac{1}{n_h}\{1-E_h[r(a-a^{\′})\left]\right\}{\widehat{\mathrm{\σ}}}_{{\mathrm{\β}}_h\left(N_h\right)}^2=F\left(n_h\right){\widehat{\mathrm{\σ}}}_{{\widehat{\mathrm{\β}}}_h\left(a\right)}\left(N_h\right)$

其中F(n_h)＝(1/n_h){1-E_h[r(a-a′)]}，E_h[r(a-a′)]为空间相关性结构，求解方法由Ripley(1981)及Rodriguez-Iturbe(1974)等提出。 

第二步：根据第一步所计算出的各层的估计均值及方差，估计总体比例均值及方差，采用Cochran提出的经典方法： 

分层的总体均值β_STR的估计值为： 

${\widehat{\stackrel{\‾}{\mathrm{\β}}}}_{\mathrm{STR}}=\frac{\underset{h}{\overset{L}{\mathrm{\Σ}}}{n}_h{\widehat{\stackrel{\‾}{\mathrm{\β}}}}_h}{n}=\underset{h=1}{\overset{L}{\mathrm{\Σ}}}{w}_h{\widehat{\stackrel{\‾}{\mathrm{\β}}}}_h$ 其中  n＝n₁+n₂+…n_L； 

均值的估计方差为： $\hat{V}\left(\widehat{\stackrel{\‾}{\mathrm{\β}}}\right)=\underset{h=1}{\overset{L}{\mathrm{\Σ}}}{w}_h^2{\hat{V}}_h\left(\widehat{\stackrel{\‾}{\mathrm{\β}}}\right)$

其中， 

为第h的方差。 

图4a显示了六种方法用于样本分层所得的估计方差随抽样比而变化的比较，这六种方法包括简单随机法、地貌图、土壤类型图及图3所示的三种非监督聚类方法(既k-means、SOM及空间共享最近邻)；而图4b则比较了三种聚类方法(既k-means、SOM及空间共享最近邻)参与分层与估计所得到的估计值的置信区间。图中可见，随抽样比从小到大变化，空间共享最近邻法得到的误差最小，而置信区间也最窄，表明估计精度最高。 

该实施例说明，采用本发明方法对现有因子进行聚类，划分空间异质模式，作为抽样调查中分层框架，可提高估计精度；而使用在抽样调查中，根据空间异质模式分层样本，在满足预定精度前提下也可少采样本，节省调查经费，提高分层效率。本发明提出的空间异质模式识别方法对我国大范围的资源调查及区划、遥感精度评估及环境监测具有重要理论意义及实用价值。 

Claims (8)

1.一种基于栅格的空间异质模式识别方法，其特征在于步骤如下：

(1.1)依据目标变量y确定影响因子变量集X；所述目标变量y是代表研究目标或估计目标的变量；所述影响因子变量集X为会对研究区域的空间异质模式产生影响的地理空间因子变量构成的集合，这些因子来自于气候气象因子，包括降雨、温度；或来自于自然地理条件，包括高程、坡度、地貌、植被；或来自于遥感影像，包括TM影像的特定波段、NDVI数据；

(1.2)对影响因子变量集X经过预处理，形成统一分辨率及坐标的多维栅格数据集；

(1.3)对形成的统一分辨率及坐标的多维栅格数据集，采用空间共享最近邻Spatial SNN算法，计算栅格单元间相似性，从而进行非监督分类，得到聚类结果；

(1.4)对步骤(1.3)得到的聚类结果建立规则库，采用规则库中的规则运算于新的多维栅格数据集，产生空间异质模式图。

2.根据权利要求1所述的基于栅格的空间异质模式识别方法，其特征在于：所述步骤(1.1)中依据目标变量y确定影响因子变量集X具体如下：

(2.1)确定所采集的备选影响因子变量x和目标变量y之间的相关性；

(2.2)如果两者相关，则保留此备选影响因子变量x，多个备选影响因子变量x形成影响因子变量集X。

3.根据权利要求2所述的基于栅格的空间异质模式识别方法，其特征在于：所述的步骤(2.1)判断备选影响因子变量x和目标变量y之间的相关性具体为：

(3.1)同时采用条件直方图、二维散点图、基于信息熵的分类贡献函数和粗糙集的属性重要性函数四种相关性分析方式来判断备选影响因子变量x和目标变量y之间的相关性；当备选影响因子变量x为连续变量时，采用二维散点图来判断x和目标变量y之间的相关性；当备选影响因子变量x为类别变量或者离散变量，或者将连续变量转化成离散变量时，采用条件直方图来判断影响因子变量x和目标变量y之间的相关性，如果条件直方图的差异显著，则说明两者之间具有相关性，保留此备选影响因子变量x；当备选影响因子为连续变量，或类别变量，或离散变量时均采用基于信息熵的分类贡献函数，即IG及GR，或粗糙集的属性重要性函数σ_(X，y)(x_i)来判断影响因子变量x和目标变量y之间的相关性；所述的基于信息熵的分类贡献函数，即IG及GR的计算公式为：

IG(y，x)＝H(y)-H(y|x)

GR(y，x)＝(H(x)-H(y|x))/H(x)

其中H(y)是熵的定义，H(y|x)是对类别化或离散化后影响因子变量x的条件熵的定义；

所述的粗糙集的属性重要性函数的计算公式如下：

${\mathrm{\σ}}_{(X,y)}\left(x_i\right)=1-\frac{\mathrm{\γ}(X-x_i,y)}{\mathrm{\γ}(X,y)}$

其中，X为影响因子变量集，x_i为计算属性重要性的影响因子变量，γ(X，y)是粗糙集中目标变量y对影响因子变量集X的依赖性；

IG、GR或σ_(X，y)(x_i)值越大，表明两者越相关，保留此备选影响因子变量，如果IG、GR或σ_(X，y)(x_i)值接近0，则表明两者不相关，移除此备选影响因子变量；

(3.2)对同时采用上述四种方式所选择的影响因子变量x采用变差膨胀因子法进行多重共线性分析，以消除冗余信息，提高影响因子变量x的信息含量。

4.根据权利要求3所述的基于栅格的空间异质模式识别方法，其特征在于：所述的采用变差膨胀因子法进行多重共线性分析，用VIF表示变差膨胀因子，则有：当VIF＞10时，计算每个影响因子变量x的VIF值，去除具有最大VIF的影响因子变量x，再进行判断计算，直到所有的影响因子变量x的VIF＜10，所得的影响因子变量集即作为影响变量因子集X。

5.根据权利要求1所述的基于栅格的空间异质模式识别方法，其特征在于：所述步骤(1.2)中预处理步骤为：不同的存储方式的影响因子变量x的预处理步骤不同：

(5.1)对于矢量数据的预处理

当影响因子变量x以矢量格式存储时，将矢量转变成栅格，形成统一的栅格数据；对于点矢量数据，采用k-最近邻或者Kriging插值算法，将点数据转换成栅格数据；对于线矢量数据采用数字微分分析法DDA或Bresenham算法转换成栅格数据；对于面矢量数据采用扫描法转换成统一分辨率的栅格数据；

(5.2)对于栅格数据的预处理

当影响因子变量x以与目标变量y的空间分辨率或坐标不一致的栅格格式存储时，采用最近邻法重采样栅格，形成统一分辨率的栅格数据；

(5.3)对步骤(5.1)和步骤(5.2)形成统一的栅格数据，如果为离散或类别数据，则将其数值化，即将这些数据转换成数值型，便于同其他的连续型数据共同聚类分析。

6.根据权利要求1所述的基于栅格的空间异质模式识别方法，其特征在于：所述步骤(1.3)所述的空间共享最近邻Spatial SNN算法如下：

(6.1)计算两栅格单元的影响因子变量集X的内积、或余弦函数或扩展Jaccard；

(6.2)根据步骤(6.1)中的计算结果构造邻近矩阵；

(6.3)根据步骤(6.2)构造的邻近矩阵，进行相似矩阵的稀疏化及构造共享最近邻图；

(6.4)根据前述得到的共享最近邻图计算每点的空间共享最近邻密度函数SpSNN_D，SpSNN_D为该点的融合了空间信息的共享最近邻密度函数；

(6.5)根据空间共享最近邻密度函数SpSNN_D识别核心点和边界点；

(6.6)聚集核心点，把核心点分成类；

(6.7)对非核心点，如果为边界点，则对边界点进行聚集，对于糙点，则移除，得到聚类结果。

7.根据权利要求6所述的基于栅格的空间异质模式识别方法，其特征在于：所述步骤(6.4)中的空间共享最近邻密度函数SpSNN_D将空间相关性作为特殊的属性信息融入到相似性计算中，具体计算如下：

x₁，x₂为两个栅格单元的影响因子矢量，NN(x₁)为由影响因子矢量经相似性计算的x₁的最近邻单元集合，size(A)为集合A中的元素个数，则分两种不同的情况采用不同的计算方式：

(7.1)若没有研究区域的目标变量的空间变异模式信息，或者无法通过相关数据采用统计学方式求得目标变量的空间变异规律，可直接通过空间连接性将空间相关性融入到相似性计算中：

$\mathrm{SpSim}(x_1,x_2)=\mathrm{size}(\mathrm{NN}\left(x_1\right)\∩\mathrm{NN}\left(x_2\right))+w_{x_1{x}_2}$

式中，

为x₁与x₂的空间连接性：x₁与x₂共享一条边为1，只共享一个点为0.5，否则为x₁与x₂不相邻，对非空间属性，距离的计算采用内积及余弦函数，可消除欧氏距离带来的不利影响；

(7.2)如果掌握了研究区域的目标变量的空间变异模式信息，即得到空间变异模型γ及其相关的参数，则可将空间相关性按照如下方式融入到相似性计算中：

$\mathrm{SpSim}(x_1,x_2)=\mathrm{size}(\mathrm{NN}\left(x_1\right)\∩\mathrm{NN}\left(x_2\right))+\tilde{C}(x_1,x_2)$

式中，

体现了x₁与x₂之间的空间关联性，由Kriging算法的空间变异模型γ所确定，

(7.3)由以上(7.1)或(7.2)的相似性的计算可得到融合空间信息的空间共享最近邻密度函数：

SpSNN_D＝|Q|，Q＝{x₂|SpSim(x₁，x₂)≥Eps}

式中，Eps为用户设定的相似性域值，大于Eps为表明x1与x2相似，否则为不相似。

8.根据权利要求1所述的基于栅格的空间异质模式识别方法，其特征在于：所述步骤(1.4)的对聚类结果建立规则库，采用规则库中的规则运算于新的数据集，产生空间异质模式图的步骤为：

(8.1)建立规则库

采用决策树、或粗糙集对步骤(1.3)中的聚类结果归纳每一类的识别规则集，并将其保存至规则库中；

(8.2)使用规则库中的规则，采用谓词逻辑运算器，或决策树运算器对新的影响因子变量集X进行运算，对其所构成的数据集进行分类，完成识别空间异质模式。

Images