CN101694671A - 一种基于地学栅格图像的空间加权主成分分析的方法 - Google Patents

Abstract

本发明一种基于地学栅格图像的空间加权主成分分析的方法，包括选定要进行主成分分析的N个栅格图像或称N个变量、定义待处理图像的加权值W、求变量间加权相关系数矩阵、根据相关矩阵求解特征根、确定主成分进行综合评价以及相对应的一套统计分析图表。本发明是利用样本位置属性给样本中增加空间权重因素来来改进相关矩阵，可以提高那些大权重(权重值接近于1)像元处的重要性，减小那些权重较小点处(权重值接近于0)像元的影响。本发明适用于地学数据如地质数据、矿产数据、地球化学异常、地球物理数据和遥感数据等进行矿产勘查和资源评价、环境污染评价、自然灾害分析等地学信息提取和地学数据挖掘工作。

Description

一种基于地学栅格图像的空间加权主成分分析的方法

技术领域

本发明涉及一种地学信息处理方法，特别是一种基于地学栅格图像的空间加权主成分分析的方法，属于地球空间信息技术领域。

技术背景

主成分分析(principal components analysis，PCA)也称主分量分析，是一种简化数据集的技术，旨在利用降维的思想，把多指标转化为少数几个综合指标。它借助于一个正交变换，将其分量相关的原随机向量转化成其分量不相关的新随机向量，这在代数上表现为将原随机向量的协方差阵变换成对角形阵，在几何上表现为将原坐标系变换成新的正交坐标系，使之指向样本点散布最开的p个正交方向，然后对多维变量系统进行降维处理，使之能以一个较高的精度转换成低维变量系统，再通过构造适当的价值函数，进一步把低维系统转化成一维系统。PCA把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。

PCA广泛应用在地学数据处理中。尤其，它已经成为一个标准的图像处理方法。它主要用来实现这两个目标：①减少相关图像的个数，由少数相互独立的主要波段组合代表多波段图像数据最大变异性。②增加多波段图像组合各个部分的可解释性。

PCA的实现需要有着多个属性的随机样本。因此，挑选一些合适的子样本以代表特定地区往往是必要的。在地学应用中，如果研究的主要对象是某种类型岩石的特征，则可能需要对岩石有关的像素做一个限制，以确保它们中的大部分反映了该岩石的类型特征。目前，商业地理信息系统软件和图像处理软件中提供了两种标准的选取像元的方法，其一是利用设置一个一定尺度的矩形窗口来使得处理分析只应用于这个窗口内的像元上；其二是设置一个多边形掩模，这样所有的图像处理过程都只在这个掩模中的像元中进行。在地理信息系统中建立一个掩模是一件十分简单的事情，比如，某种岩层或者某种岩石类型都可以通过在地理信息系统中借助查询功能选中，然后进行将其定义为掩模即可。不过这类方法都可以被认为是利用了一个有着二值(0和1)的二值掩模。前一种方法可以看作是后一种方法当掩模形状变为规范形状后的一个特例。毫无疑问的是，在某些应用之中对这两种方法进行一些约束以达到去掉一些分析中不相关的样本是必要的。

发明内容

本发明的技术解决问题：克服现有的PCA方法的不足，提供一种基于地学栅格图像的空间加权主成分分析的方法，该方法利用样本位置属性给样本中增加空间权重因素来来改进相关矩阵，以提高样本中特殊位置的重要性。

本发明的技术方案：一种基于地学栅格图像的空间加权主成分分析的方法包括如下步骤：

第一步，选定要进行主成分分析的P个栅格图像；

第二步，定义待处理图像的加权值W；

第三步，构建变量间的加权相关系数矩阵；

由第一步和第二步所选定的P个栅格图像和加权值W，计算加权相关系数矩阵R，作为第四步计算特征根和特征向量的数据基础；

第四步，根据相关系数矩阵，求解特征根和特征向量；

根据第三步所获得的加权相关系数矩阵R，运用线性代数求解线性方程组的方法求解加权相关系数矩阵的特征根和特征向量，并获得所有的主成分表达式；

第五步，根据特征根和特征向量的统计特征，生成相关的统计分析图表；

由第四步计算所获得的特征根和特征向量，计算主成分的累积贡献率、因子载荷，并生成对应的特征根分布和累计贡献图、因子载荷图、主成分对栅格图像贡献图和主成分分析的统计报表；

第六步，确定并生成主成分综合评价图；

通过对第五步生成的图表进行分析，并与所研究的问题相结合，确定所选择的主成分，再根据第四步选择对应的主成分表达式计算因子得分，并生成主成分综合评价因子得分图。

本发明的优点和积极效果是：本发明利用样本位置属性给样本中增加空间权重因素来来改进相关矩阵，可以提高那些大权重像元处的重要性，减小那些权重较小点处像元的影响，同时，一整套统计图表为主成分的分析和解释提供了极大的方便。本发明适用于地学数据如地质数据、矿产数据、地球化学异常、地球物理数据和遥感数据等进行矿产勘查和资源评价、环境污染评价、自然灾害分析等地学信息提取和地学数据挖掘工作。

附图说明

图1为本发明的主流程图；

图2加权主成分分析主界面和栅格图像选定图；

图3选定加权图像界面；

图4输出参数设置界面图；

图5特征值分布和累计贡献；

图6因子载荷图；

图7主成分对栅格图像贡献图；

图8第二主成分因子得分图。

具体实施方式

如图1所示，本发明的具体实施步骤如下：

第一步，选定要进行主成分分析的N个栅格图像或称N个变量；

首先进行数据准备，选择要分析的地学数据，将采集地学数据，如点数据，图像数据(TIF，JPG，BMP)等文件格式的数据通过内插和格式转换到栅格图像文件格式，这些方法在一些常用的商业GIS软件中都可以获得，如ArcGIS等。例如，Au_IDW是金属元素金的湖泊沉积物采样内插得到的栅格图像。图2是栅格图像加权主成分分析主界面及栅格图像选择界面，在图2一开始左边的列表框显示所有的栅格图像，将要进行主成分分析的栅格图像选定移到右边的列表框中，右边的列表框所显示的栅格图像，如图中的，w_idw，Cu_idw，au_idw和as_idw分别是地球化学元素W、Cu、Au和As的插值后元素含量的分布，以此四个元素为例进行主成分分析。

第二步，定义待处理图像的加权值W，在本发明中对于W加权证的赋值不给出具体的规定，因为加权值的取值与具体的问题有密切的关系，但对于加权值给出了必要满足的要求：即加权值是一幅栅格图像，它的取值在0到1的范围内(0≤W_ij≤1)。当定义了加权栅格图像就可以通过图3选择该栅格作为加权图像，图3是选择加权图像的界面，如weightl是用作加权的图像，该图层将被用来在第四步进行相关系数计算。

第三步，设置加权主成分分析的输出结果参数，在加权主成分分析中有一整套统计分析图表，设置输出参数可选择性地输出部分图表，图4是输出参数选择见面。在该界面中可以根据进行主成分分析的目的有选择地输出相关的统计图表，如在地球化学数据处理时，希望通过主成分分析方法寻找有利于成矿的地球化学元素组合，那么因子载荷图是必须输出的图表，通过该图表可以清楚的选择有利于成矿的地球化学元素组合。

第四步，根据选定的栅格图像和加权栅格图像，按照加权相关系数的定义构建变量间的加权相关系数矩阵，并根据加权相关系数矩阵求解特征值和对应的特征向量，其步骤如下：

(1)设X₁，X₂，…，X_p是p个选定的栅格图像，x_ijk(i＝1，…，p，l＝1，…，N，k＝1，..，M)表示第i个栅格图像在位置(l，k)处的取值，N为栅格图像的宽度，M为栅格图像的高度。W为加权栅格图像，w_lk为加权图像在位置(l，k)处的加权值。

(2)计算加权相关系数R

加权相关系数R＝(r_ij)_N×M，其中。

$r_{\mathrm{ij}}=\frac{\mathrm{\Σ}{w}_{\mathrm{ij}}(x_{\mathrm{ilk}}-\stackrel{\‾}{x_i})(x_{\mathrm{jlk}}-\stackrel{\‾}{x_j})}{\sqrt{\mathrm{\Σ}{w}_{\mathrm{ij}}{(x_{\mathrm{ilk}}-\stackrel{\‾}{x_i})}^2}\sqrt{\mathrm{\Σ}{w}_{\mathrm{ij}}{(x_{\mathrm{ilk}}-\stackrel{\‾}{x_j})}^2}}$ i，j＝1，...，p

(3)求解特征根和特征向量

设加权相关系数矩阵R的特征值为λ₁≥λ₂≥…≥λ_p≥0，相应的单位特征向量分别为(a₁，a₂，…a_p)。这时第i个主成份为

i＝1，…p，且Var(Y_i)＝λ_i，i＝1，…，p，λ_i反映了第i个主成分的方差，其顺序反映了该主成分在所有主成分中重要性的顺序，比值说明了第i个主成分的贡献率，也可以说是第i个主成分包含了原始变量信息量的百分比。单位特征向量a_i＝(a_i1，a_i2，…，a_ip)^τ为第i个主成分中各个变量的加权系数，a_ik(k＝1，…，p)的值的大小反映了第i个主成分对第k个变量X_k的依赖程度。第i个主成分可表示为：

Y_i＝a_i1X₁+a_i2X₂+…+a_ipX_p，i＝1，…，p。

记X＝(X₁，X₂，…，X_p)，A＝(a₁，a₂，…a_p)，

Y＝(Y₁，…Y_p)，则A为正交阵，写成矩阵形式为：

Y＝A′X，Var(Y)＝Λ，且

主成份Y_i与原始变量X_k的相关系数ρ(Y_i，X_k)称为因子载荷(factor loading)，反映了第k变量在第i主成分上的重要性，可以证明

k，i＝1，…，p， $\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\ρ}}^2(Y_i,X_k)={\mathrm{\λ}}_i,$ $\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\ρ}}^2(Y_i,X_k)=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{a}_{\mathrm{ki}}^2.$

前i个主成分的累计贡献率为：

i＝1，…p，反映了前i个主成分对原始数据的贡献率。

第五步，根据设置的输出参数显示并分析各种统计图表。这些图表从不同侧面对主成分分析的结果给予分析和解释。

(1)特征值分布和累计贡献图，该图绘制了所有特征值从大到小的分布和累计贡献，以利于主成分个数的确定，如图5，该图由一个柱形图和一个线型图组成，柱形图由第四步中所得到的特征根λ₁，λ₂，…λ_p构成，反应了各个主成分的方差重要性，线型图由前i个主成分的累计贡献率v_i，…v_p组成，在该实施例中，第一主成分的方差为1.8，其贡献为44.8％，第二主成分的方差为1.02，第一、二主成分的累计贡献为70.3％。该图通常用于确定主成分的个数q，其遵循的原则为：

要求λ_q＞1，或者v_q≥70％，或者根据分析该图人为指定q的值。

(2)因子载荷图，该图定义了各主成分与各个地学栅格图像柱形图，主成分为柱形图的变量，因子载荷是其数据，清晰地反映了各主成分与各个地学栅格图像柱形图的关系，如图6，该图是一个层叠式多系列柱形图，其值由第四步得到的因子载荷构成，每一个主成分与所有变量(栅格图层)的相关系数组成一个系列，系列数由主成分的个数确定，在该实施例中，图中红色为第一主成分与所有变量的因子载荷组成，说明第一主成分与所有栅格图像都是正相关，其中as_idw和W_iwd的因子载荷较大；淡粉色为第二主成分与所有变量的因子载荷组成，说明第二主成分与as_idw和W_iwd为负相关，而与cu_idw和au_idw为正相关，由第二主成分所获得的因子得分其高值反映了cu_idw和au_idw组合，低值反映了as_idw和W_iwd组合。

(3)主成分对栅格图像贡献图，该图定义了各个地学栅格图像与主成分的柱形图，地学栅格图像为柱形图的变量，各主成分对地学栅格图像的贡献是数据，该图反映了各个主成分对地学栅格图像的贡献，如图7，该图是一个层叠式多系列柱形图，每个变量(栅格图像)与各个主成分的相关系数组成一个系列，系列数是变量数p，在该实施例中，该图中两个柱子表示第一主成分和第二主成分对四个栅格图像的贡献，例如，从图中可以看出，红色表示栅格图像W_idw与主成分的关系，W_idw与第一主成分是正相关，与第二主成分是负相关；绿色表示au_idw与主成分的关系，au_idw与两个主成分都是正相关，但第二主成分更好地反映了au_idw。

(4)进行加权主成分分析的统计报表，该报表包含在进行主成分分析的过程中所有的统计参数，包括，所有的特征值及其对应的特征向量、因子载荷等，见表1，该表把进行主成分分析的统计参数列在表中，以便于分析。

表1主成分分析统计报表

变量	w_idw	cu_idw	au_idw	as_idw
					均值	1.45384	10.24989	3.80004	9.74379
方差	8.37268	3.90419	2.57756	11.8657
相关系数
						w_idw	cu_idw	au_idw	as_idw
w_idw	1	0.11409	0.06240	0.53990
					cu_idw	0.11409	1.00000	0.17418	0.42897
au_idw	0.06240	0.17418	1.00000	0.10753
					as_idw	0.53990	0.42897	0.10753	1
					主成分的重要性
特征值(主成分方差)	1.7941	1.0193	0.8244	0.3622
					标准差	1.3395	1.0096	0.9080	0.6018
主成分的相对贡献	0.4485	0.2548	0.2061	0.0906
					主成分的累计贡献	0.4485	0.7033	0.9095	1.0000
					因子载荷
	第一主成分	第二主成分
					w_idw	0.5327	-0.4536

变量	w_idw	cu_idw	au_idw	as_idw
					cu_idw	0.4818	0.3641
au_idw	0.4818	0.4818
					as_idw	0.6544	-0.1847

第六步，生成主成分综合评价图，通过对第六步所生成的统计图表进行分析，根据所研究的问题确定要选择的主成分，由选定的主成分计算因子得分如第二主成分是选定的主成分i0＝2，

并主成分综合评价生成栅格图像。在实际应用中，可能会根据不同的分析目的生成多个主成分因子得分图。由于在本实施例中主要对金元素进行评价，而第二主成分的主要反映了金和铜组合，所以选择第二主成分对金元素进行综合评价。图8是第二主成分的因子得分图。

总之，本发明是利用样本位置属性给样本中增加空间权重因素来来改进相关矩阵，可以提高那些大权重(权重值接近于1)像元处的重要性，减小那些权重较小点处(权重值接近于0)像元的影响。本发明适用于地学数据如地质数据、矿产数据、地球化学异常、地球物理数据和遥感数据等进行矿产勘查和资源评价、环境污染评价、自然灾害分析等地学信息提取和地学数据挖掘工作。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (3)

1.一种基于地学栅格图像的空间加权主成分分析的方法，其特征在于包括如下步骤：

第一步，选定要进行主成分分析的P个栅格图像；

第二步，定义待处理图像的加权值W；

第三步，构建变量间的加权相关系数矩阵

由第一步和第二步所选定的P个栅格图像和加权值W，计算加权相关系数矩阵R，作为第四步计算特征根和特征向量的数据基础；

第四步，根据相关系数矩阵，求解特征根和特征向量；

根据第三步所获得的加权相关系数矩阵R，运用线性代数求解线性方程组的方法求解加权相关系数矩阵的特征根和特征向量，并获得所有的主成分表达式；

第五步，根据特征根和特征向量的统计特征，生成相关的统计分析图表；

由第四步计算所获得的特征根和特征向量，计算主成分的累积贡献率、因子载荷，并生成对应的特征根分布和累计贡献图、因子载荷图、主成分对栅格图像贡献图和主成分分析的统计报表；

第六步，确定并生成主成分综合评价图；

通过对第五步生成的图表进行分析，并与所研究的问题相结合，确定所选择的主成分，再根据第四步选择对应的主成分表达式计算因子得分，并生成主成分综合评价因子得分图。

2.根据权利要求1所述的基于地学栅格图像的空间加权主成分分析的方法，其特征在于：所述第二步中定义待处理图像的加权值W是距离某个矿床的距离大小，或者矿床的密度，距离某个断裂带的距离，或者某微量元素的浓度值。

3.根据权利要求1所述的基于地学栅格图像的空间加权主成分分析的方法，其特征在于：所述第三步中构建变量间的加权相关系数矩阵，加权值W被定义成了一幅加权值图像W，它的数值在0到1的范围内，假定A，B和C是三幅要进行处理的图像，它们在坐标点(i，j)点的值分别为A_ij，B_ij，C_ij，加权值W被定义成了一幅加权值图像W_ij，它的数值在0到1的范围内，0≤W_ij≤1，两幅图像，即图像A和图像B的加权后的相关系数矩阵被定义为：

$R(A,B)=\frac{\mathrm{\Σ}{W}_{\mathrm{ij}}(A_{\mathrm{ij}}-\stackrel{\‾}{A})(B_{\mathrm{ij}}-\stackrel{\‾}{B})}{\sqrt{\mathrm{\Σ}{W}_{\mathrm{ij}}{(A_{\mathrm{ij}}-\stackrel{\‾}{A})}^2}\sqrt{\mathrm{\Σ}{W}_{\mathrm{ij}}{(B_{\mathrm{ij}}-\stackrel{\‾}{B})}^2}}$

其中字母A，B分别代表的是A_ij，，B_ij的平均值，可以看出增加了权重的相关系数具有下面的属性：

(1)R(A，B)是对称的；

(2)-1≤R(A，B)≤1；

(3)R(A，B)＝1或者-1是当A与B成正线性相关或者成负线性相关时的情况；

(4)如果W_ij是常数的话，R(A，B)则就变成了普通的相关系数；

(5)如果W_ij是一幅有着0和1两个数值的二值图，那么W就相当于普通的一个掩模；

前三个属性保证该定义(1)满足相关系数的必要要求，后两个属性则表明定义(1)是一个比普通的相关系数公式和应用了一个二值掩模的相关系数公式的更为一般的形式，普通的掩模处理变成了这个空间加权情况的一个特例。

Images