一种基于因子分析模型的高光谱数据降维方法

摘要

一种基于因子分析模型的高光谱数据降维方法，其步骤如下：(1)高光谱数据读入；(2)建立高光谱数据降维的因子分析模型；(3)计算数据的均值、协方差矩阵以及相关矩阵；(4)计算数据相关矩阵的特征值和标准化特征向量；(5)由参数估计方法进行因子载荷矩阵的求解；(6)计算因子分析模型中特殊因子的协方差和数据变量的共同度；(7)计算基于方差最大的因子载荷旋转矩阵；(8)利用基于加权最小二乘方法计算因子得分；(9)得到表征高光谱数据的本征维数，实现高光谱数据降维。该方法是一种自动的高光谱数据降维方法，能够有效的去除高光谱数据波段间的相关性、增加地物不同类别间的可分性。

说明书

技术领域

本发明涉及一种基于因子分析模型的高光谱数据降维方法，属于高光谱数据处理方法与应用技术领域，适用于高光谱数据降维的理论方法和应用技术研究。

背景技术

高光谱成像仪是一种新型的遥感载荷，其光谱具有紧密、连续的特点，可以同时记录被测同一地物的光谱和空间信息特征，因此，高光谱数据维数高且数据量大；由于光谱分辨率较高，各个波段之间具有较高的相关性，存在了大量的冗余信息。并且随着波段数据的增加，数据处理量呈指数上升。因此，如何将高光谱数据降维，消除波段间的相关性，成为数据处理的关键环节。

目前针对高光谱数据降维已经提出了很多方法，主要分为波段选择和特征提取两大类。波段选择的目的是将n维高光谱数据用具有代表意义的m维(m<n)数据子集代替，该类方法的问题是虽然较好的保留了原图像数据的特性，但是损失了高光谱数据的信息量以及光谱细节信息；特征提取的目的是将n维高光谱数据压缩到m维(m<n)模式空间中，使类别间的可分性最大。因此，该类方法在数据降维的同时突出了图像数据间的差异性，但改变了数据的原有特性。目前数据降维方法主要存在的问题是：不能很好的保留数据提供的所有有用信息，不能将相同类型地物特征表示在相同的波段中。

发明内容

本发明的目的是针对现有高光谱数据降维方法损失数据信息等不足提出一种基于因子分析模型的高光谱数据降维方法。

本发明的技术解决方案是：本发明一种基于因子分析模型的高光谱数据降维方法，具体是指一种因子分析模型将高光谱高维空间数据降到表征高光谱数据本征特性的低维空间的方法。该方法主要是利用因子分析模型，并通过计算因子载荷矩阵、基于方差最大的因子载荷旋转、计算因子得分得到表征高光谱数据的本征维数，从而实现高光谱数据降维。

本发明一种基于因子分析模型的高光谱数据降维方法，其步骤如下：

(1)高光谱遥感数据读入；

(2)建立高光谱数据降维的因子分析模型；

(3)计算高光谱数据的均值、协方差矩阵以及相关矩阵；

(4)计算数据相关矩阵的特征值和标准化特征向量；

(5)由主成分解方法进行因子载荷矩阵的求解；

(6)计算因子分析模型中特殊因子的方差矩阵和数据变量的共同度；

(7)计算基于方差最大的因子载荷旋转矩阵；

(8)利用基于加权最小二乘方法计算因子得分；

(9)得到表征高光谱数据的本征维数，实现高光谱数据降维。

其中，步骤(2)中所述的“建立高光谱数据降维的因子分析模型”，其因子分析模型为：

X＝μ+AF+ε

式中，X＝(x₁，x₂，...，x_p)′是可观测的随机变量，均值μ＝E(X)，X的公共因子F＝(F₁，F₂，...，F_m)′(m<p)是不可观测的随机变量，且满足均值E(F)＝0、方差D(F)＝I_m(即F各个分量方差为1且互不相关)，特殊因子ε＝(ε₁，ε₂，...，ε_p)′与F互不相关，且满足均值E(ε)＝0、方差$D\left(ϵ\right)=\mathrm{diag}({{\sigma }_1}^2,{{\sigma }_2}^2,...,{{\sigma }_p}^2)\stackrel{\mathrm{def}}{=}D;$公共因子一般对X每一个分量都起作用，而特殊因子只对X的某一分量起作用，而且特殊因子之间以及特殊因子与公共因子之间互不相关；因子载荷矩阵A＝(a_ij)_p×m是待估计的系数矩阵，在因子分析中，特殊因子起残差或噪声的作用。

其中，步骤(3)所述的“计算高光谱数据的均值、协方差矩阵以及相关矩阵”，其计算方法如下：

$\bar{X}=\frac{1}{n}\underset{t=1}{\overset{n}{\Sigma }}{X}_{\left(t\right)}=({\bar{x}}_1,{\bar{x}}_2,...,{\bar{x}}_p)\prime$

$V=\frac{1}{n}\underset{t=1}{\overset{n}{\Sigma }}(X_{\left(t\right)}-\bar{X})(X_{\left(t\right)}-\bar{X})\prime \stackrel{\mathrm{def}}{=}{V}_{\mathrm{ij}}$

$R=\left(r_{\mathrm{ij}}\right)=\left(\frac{V_{\mathrm{ij}}}{\sqrt{V_{\mathrm{ii}}}\sqrt{V_{\mathrm{jj}}}}\right);(i,j=\mathrm{1,2},...,p)$

式中，X为数据的均值；V为协方差矩阵；R为相关矩阵；n表示像元数目；X_(t)表示第t个像元。

其中，步骤(4)中所述的“计算数据相关矩阵的特征值和标准化特征向量”，其含义说明如下：计算相关矩阵R的特征值为λ₁≥λ₂≥…≥λ_p≥0，以及各个特征值对应的施密特正交化的特征向量l₁，l₂，…，l_p。

其中，步骤(5)中所述的“由主成分解方法进行因子载荷矩阵的求解”，其计算方法如下：

首先根据$\frac{{\lambda }_1+{\lambda }_2+···+{\lambda }_m}{{\lambda }_1+···+{\lambda }_m+···{\lambda }_p}\ge 0.9$的最小整数，确定公共因子个数m；

然后利用主成分解的方法求解因子载荷矩阵：利用步骤(4)中计算得到的相关矩阵R的特征值为λ₁≥λ₂≥…≥λ_p≥0，对应的施密特正交化后的特征向量为l₁，l₂，…，l_p，则利用相关矩阵R的谱分解：$R=\underset{i=1}{\overset{p}{\Sigma }}{\lambda }_i{l}_i{l}_i^{\prime },$得到因子分析模型的一个解：$A=(\sqrt{{\lambda }_1}{l}_1,···,\sqrt{{\lambda }_m}{l}_m)\stackrel{\mathrm{def}}{=}{\left(a_{\mathrm{ij}}\right)}_{p\times m}.$

其中，步骤(6)中所述的“计算因子分析模型中特殊因子的方差矩阵和数据变量的共同度”，其含义说明如下：因子分析模型中特殊因子的方差为：${\sigma }_i^2=1-\underset{t=1}{\overset{m}{\Sigma }}{a}_{\mathrm{it}}^2,(i=\mathrm{1,2},...,p),$变量x_i的共同度h_i²估计为$h_i^2=\underset{t=1}{\overset{m}{\Sigma }}{a}_{\mathrm{it}}^2,(i=\mathrm{1,2},...,p).$

其中，步骤(7)中所述的“计算基于方差最大的因子载荷旋转矩阵”，其含义说明如下：对公共因子做正交旋转就是对载荷矩阵A作正交变换，右乘正交矩阵Г，使得AΓ具有更鲜明的意义，旋转后的公共因子向量为F^*＝Г′F，它的各个分量也是互不相关的公共因子；根据正交矩阵Г的不同选取方式，将构造出不同的正交旋转方法，本发明中采用的方法是最大方差的旋转方法；具体实现过程如下：

令$A^{*}=\mathrm{A\Gamma }={\left(a_{\mathrm{ij}}^{*}\right)}_{p\times m},$$d_{\mathrm{ij}}=a_{\mathrm{ij}}^{*}/h_i,$${\bar{d}}_j=\frac{1}{p}\underset{i=1}{\overset{p}{\Sigma }}{d}_{\mathrm{ij}}^2,$则A^*的第j列元素平方的相对方差可以定义为：$V_j=\frac{1}{p}\underset{i=1}{\overset{p}{\Sigma }}{(d_{\mathrm{ij}}^2-{\bar{d}}_j)}^2,$式中，取是为了消除d_ij符号不同的影响；用除以h_i是为了消除各个原始变量X_i对公共因子依赖程度不同的影响；最大方差旋转方法就是选择正交矩阵Г，使A^*所有m列元素平方的相对方差V达到最大：V＝V₁+V₂+…+V_m；逐次对每两个公共因子进行上述旋转，对公共因子F_l和F_k旋转，就是对A矩阵的第1列和第k列进行正交变换，使这两列的元素平方的相对方差之和最大，而其余列不变，旋转矩阵Г_lk可以表示为：

式中，θ是F_l和F_k的旋转角度，矩阵其余元素为0；求θ使V达到最大，由微积分中求极值的方法，将V对θ求导，并令其为0：

$V=V_1+V_2+···+V_m=\frac{1}{p}\underset{i=1}{\overset{p}{\Sigma }}{(d_{i1}^2-{\bar{d}}_1)}^2+···+\frac{1}{p}\underset{i=1}{\overset{p}{\Sigma }}{(d_{\mathrm{im}}^2-{\bar{d}}_m)}^2$

所以，由$\frac{\partial V}{\partial \theta }=0,$得：$\tan 4\theta =\frac{M-2\mathrm{UV}/p}{W-(U^2-V^2)/p},$式中：$U=\underset{i=1}{\overset{p}{\Sigma }}{u}_i,$$V=\underset{i=1}{\overset{p}{\Sigma }}{v}_i,$$W=\underset{i=1}{\overset{p}{\Sigma }}({u_i}^2-v_i^2),$$M=2\underset{i=1}{\overset{p}{\Sigma }}{u}_i{v}_i,$$u_i={\left(\frac{a_{i1}}{h_i}\right)}^2+{\left(\frac{a_{i2}}{h_i}\right)}^2,$$v_i=2\frac{a_{i1}{a}_{i2}}{{h_i}^2},$a_i1、a_t2是因子载荷矩阵中的元素，h_i²是变量x_i的共同度。

第一轮旋转后的因子载荷矩阵为A⁽¹⁾，然后开始下一轮的旋转，得到一系列的因子载荷矩阵为：A⁽¹⁾，A⁽²⁾，...，A^(s)，...，则必有：V⁽¹⁾≤V⁽²⁾≤...≤V^(s)≤...，V^(s)为A^(s)各列元素平方的相对方差之和，实际应用中，当V^(s-1)-V^(s)≤ε时，即可停止旋转。

其中，步骤(8)中所述的“利用基于加权最小二乘方法计算因子得分”，其含义说明如下：在得到A和D的基础上，采用加权最小二乘对F进行求解：$\hat{F}={(\hat{A}\prime {\hat{D}}^{-1}\hat{A})}^{-1}A\prime D^{-1}X,$即得到了F的加权最小二乘估计，式中，是估计得到的因子载荷矩阵，是估计得到的特殊因子方差矩阵，X是读入数据。

其中，步骤(9)中所述的“得到表征高光谱数据的本征维数，实现高光谱数据降维”，其含义说明如下：利用因子得分估计即可得到表征高光谱数据的本征维数，从而实现了高光谱数据降维。

本发明与现有技术相比的优点在于：克服了现有高光谱数据降维方法损失数据信息、无法保留数据原有特性等局限，本方法利用了因子分析模型，并采用了基于方差最大的因子矩阵的正交旋转，得到了表征高光谱数据的本征维数，实现了高光谱数据降维。它具有以下的优点：(1)利用因子分析模型对高光谱数据进行建模求得反映高光谱数据结构与特征的本征维数，消除波段间的相关性的同时简化了数据结构，较好的保持了数据原有特性；(2)利用了基于方差最大的因子载荷旋转，易于各个因子所代表含义的解释，且每种地物只在一个因子中起主导作用，最大程度的增大了不同类型地物之间的可分性。

附图说明

图1为本发明涉及的一种基于因子分析模型的高光谱数据降维方法的实现流程

具体实施方式

为了更好的说明本发明涉及的基于因子分析模型的高光谱数据降维方法，利用PHI航空高光谱数据进行江苏方麓茶场地区农作物精细分类。本发明一种基于因子分析模型的高光谱数据降维方法，实现流程如图1所示，具体实现步骤如下：

(1)高光谱数据的读入：读入方麓茶场PHI高光谱数据，去掉信噪比低、大气吸收等波段，原数据大小为210×150×64；

(2)建立高光谱数据降维的因子分析模型；

高光谱数据降维的因子分析模型为：

X＝μ+AF+ε

式中，X＝(x₁，x₂，...，x₆₄)′是可观测的随机变量，均值μ＝E(X)，X的公共因子F＝(F₁，F₂，...，F_m)′(m<64)是不可观测的随机变量，且满足均值E(F)＝0、方差D(F)＝I_m(即F各个分量方差为1且互不相关)，特殊因子ε＝(ε₁，ε₂，...，ε₆₄)′与F互不相关，且满足均值E(ε)＝0、方差$D\left(ϵ\right)=\mathrm{diag}({{\sigma }_1}^2,{{\sigma }_2}^2,...,{{\sigma }_{64}}^2)\stackrel{\mathrm{def}}{=}D;$公共因子一般对X每一个分量都起作用，而特殊因子只对X的某一分量起作用，而且特殊因子之间以及特殊因子与公共因子之间互不相关；矩阵A＝(a_ij)_p×m是待估计的系数矩阵，称为因子载荷矩阵，在因子分析中，特殊因子起残差或噪声的作用；

(3)计算高光谱数据的均值、协方差矩阵以及相关矩阵；

210×150×64的高光谱数据均值X、协方差矩阵V、相关矩阵R的计算方法如下：

$\bar{X}=\frac{1}{n}\underset{t=1}{\overset{n}{\Sigma }}{X}_{\left(t\right)}=({\bar{x}}_1,{\bar{x}}_2,...,{\bar{x}}_{64})\prime$

$V=\frac{1}{n}\underset{t=1}{\overset{n}{\Sigma }}(X_{\left(t\right)}-\bar{X})(X_{\left(t\right)}-\bar{X})\prime \stackrel{\mathrm{def}}{=}{V}_{\mathrm{ij}}$

$R=\left(r_{\mathrm{ij}}\right)=\left(\frac{V_{\mathrm{ij}}}{\sqrt{V_{\mathrm{ii}}}\sqrt{V_{\mathrm{jj}}}}\right);(i,j=\mathrm{1,2},...,64)$

式中，X为数据的均值；V为协方差矩阵；R为相关矩阵；n＝210×150表示像元数目；X_(t)表示第t个像元。

(4)计算数据相关矩阵的特征值和标准化特征向量；

计算相关矩阵R的特征值为λ₁≥λ₂≥…≥λ₆₄≥0，对应的经施密特正交化后的特征向量为l₁，l₂，…，l₆₄；

(5)由主成分解方法进行因子载荷矩阵的求解；

首先根据$\frac{{\lambda }_1+{\lambda }_2+···+{\lambda }_m}{{\lambda }_1+···+{\lambda }_m+···{\lambda }_{64}}\ge 0.9$的最小整数确定公共因子个数m＝6；

然后利用主成分解的方法求解因子载荷矩阵：利用步骤(4)中计算得到的相关矩阵R的特征值为λ₁≥λ₂≥…≥λ₆₄≥0，对应的经施密特正交化后的特征向量为l₁，l₂，…，l₆₄，则利用相关矩阵R的谱分解：$R=\underset{i=1}{\overset{64}{\Sigma }}{\lambda }_i{l}_i{l}_i^{\prime },$得到因子分析模型的一个解，$A=(\sqrt{{\lambda }_1}{l}_1,···,\sqrt{{\lambda }_6}{l}_6)\stackrel{\mathrm{def}}{=}{\left(a_{\mathrm{ij}}\right)}_{64\times 6};$

(6)计算因子分析模型中特殊因子的方差矩阵和数据变量的共同度；

因子分析模型中特殊因子的方差为：${\sigma }_i^2=1-\underset{t=1}{\overset{m}{\Sigma }}{a}_{\mathrm{it}}^2,(i=\mathrm{1,2},...,64),$变量x_i的共同度h_i²估计为$h_i^2=\underset{t=1}{\overset{m}{\Sigma }}{a}_{\mathrm{it}}^2,(i=\mathrm{1,2},...,64);$

(7)计算基于方差最大的因子载荷旋转矩阵；

对公共因子做正交旋转就是对载荷矩阵A作正交变换，右乘正交矩阵Γ，使得AΓ具有更鲜明的意义，旋转后的公共因子向量为F^*＝Г′F，它的各个分量也是互不相关的公共因子；根据正交矩阵Г的不同选取方式，将构造出不同的正交旋转方法，本发明中采用的方法是最大方差的旋转方法；具体实现过程如下：

令$A^{*}=\mathrm{A\Gamma }={\left(a_{\mathrm{ij}}^{*}\right)}_{64\times 6},$$d_{\mathrm{ij}}=a_{\mathrm{ij}}^{*}/h_i,$${\bar{d}}_j=\frac{1}{64}\underset{i=1}{\overset{64}{\Sigma }}{d}_{\mathrm{ij}}^2,$则A^*的第j列元素平方的相对方差可以定义为：$V_j=\frac{1}{64}\underset{i=1}{\overset{64}{\Sigma }}{(d_{\mathrm{ij}}^2-{\bar{d}}_j)}^2,$式中，取是为了消除d_ij符号不同的影响；用除以h_i是为了消除各个原始变量X_i对公共因子依赖程度不同的影响；选择除以h_i是因为A^*的第i行平方和：

$h_i^{*2}=\underset{j=1}{\overset{64}{\Sigma }}{a}_{\mathrm{ij}}^{*2}=(a_{i1}^{*},a_{i2}^{*},...,a_{i64}^{*})\left(\begin{array}{c}{a}_{i1}^{*}\\ a_{i2}^{*}\\ ·\\ ·\\ ·\\ a_{i64}^{*}\end{array}\right)=(a_{i1},a_{i2},...,a_{i64})\mathrm{\Gamma \Gamma }\prime \left(\begin{array}{c}{a}_{i1}\\ a_{i2}\\ ·\\ ·\\ ·\\ a_{i64}\end{array}\right)=\underset{j=1}{\overset{64}{\Sigma }}{a}_{\mathrm{ij}}^2=h_i^2$

最大方差旋转方法就是选择正交矩阵Г，使A^*所有6列元素平方的相对方差V达到最大：V＝V₁+V₂+…+V₆；逐次对每两个公共因子进行上述旋转，对公共因子F_l和F_k旋转，就是对A矩阵的第1列和第k列进行正交变换，使这两列的元素平方的相对方差之和最大，而其余列不变，旋转矩阵Г_lk可以表示为：

式中，θ是F_l和F_k的旋转角度，矩阵其余元素为0；求θ使V达到最大，由微积分中求极值的方法，将V对θ求导，并令其为0：

$V=V_1+V_2+···+V_6=\frac{1}{64}\underset{i=1}{\overset{64}{\Sigma }}{(d_{i1}^2-{\bar{d}}_1)}^2+···+\frac{1}{64}\underset{i=1}{\overset{64}{\Sigma }}{(d_{i6}^2-{\bar{d}}_6)}^2$

所以，由$\frac{\partial V}{\partial \theta }=0,$得：$\tan 4\theta =\frac{M-2\mathrm{UV}/64}{W-(U^2-V^2)/64}$，式中：$U=\underset{i=1}{\overset{64}{\Sigma }}{u}_i,$$V=\underset{i=1}{\overset{64}{\Sigma }}{v}_i,$$W=\underset{i=1}{\overset{64}{\Sigma }}({u_i}^2-v_i^2),$$M=2\underset{i=1}{\overset{64}{\Sigma }}{u}_i{v}_i,$$u_i={\left(\frac{a_{i1}}{h_i}\right)}^2+{\left(\frac{a_{i2}}{h_i}\right)}^2,$$v_i=2\frac{a_{i1}{a}_{i2}}{{h_i}^2},$a_t1、a_i2是因子载荷矩阵中的元素，h_i²是变量x_i的共同度。

第一轮旋转后的因子载荷矩阵为A⁽¹⁾，然后开始下一轮的旋转，得到一系列的因子载荷矩阵为：A⁽¹⁾，A⁽²⁾，...，A^(s)，...，则必有：V⁽¹⁾≤V⁽²⁾≤...≤V^(s)≤...，V^(s)为A^(s)各列元素平方的相对方差之和，实际应用中，当V^(s-1)-V^(s)≤ε(本发明中ε≤0.001)时，即可停止旋转；

(8)利用基于加权最小二乘方法计算因子得分；

在得到因子载荷矩阵A和特殊因子方差矩阵D的基础上，采用加权最小二乘对F进行求解：

$\hat{F}={(\hat{A}\prime {\hat{D}}^{-1}\hat{A})}^{-1}A\prime D^{-1}X,$

即得到了F的加权最小二乘估计；

(9)得到表征高光谱数据的本征维数，实现高光谱数据降维。

利用因子得分估计即可得到表征高光谱数据的本征维数，从而实现高光谱数据降维；得到了6个表征所使用高光谱数据的所有信息的本征波段。