孟令博, 耿修瑞, 杨炜暾
中国科学院电子学研究所 中国科学院空间信息处理与应用系统技术重点实验室, 北京 100190; 中国科学院大学, 北京 100049
2017年11月28日 收稿; 2018年3月27日 收修改稿
基金项目: 国家自然科学基金（41601402，41701539）资助
通信作者: 耿修瑞, E-mail:xrgeng@mail.ie.ac.cn

摘要: 基于峭度指标的FastICA算法具有较快的收敛速度和较高的计算效率，被广泛应用于多光谱图像的特征提取。经典的FastICA算法基于固定点迭代法得到图像的各个独立成分，在迭代过程中，每一个独立成分的求解都需要所有像元的参与。因此，当数据量较大或图像中像元较多时，FastICA的计算量很大，此时它的速度优势就会大打折扣。遥感数据一般都具有较大的尺寸，因此如何将FastICA直接应用于遥感数据，是一个具有实际意义的问题。通过引入多光谱图像协峭度张量的概念，将FastICA的固定点迭代问题转化为代数形式的张量计算，避免每次迭代过程中需所有像元参与的缺陷，因而大大降低计算复杂度。多光谱图像实验结果表明，该算法明显快于传统的基于峭度指标的FastICA算法。
关键词: FastICA高阶统计特性峭度指标协峭度张量
An improved FastICA algorithm based on kurtosis index
MENG Lingbo, GENG Xiurui, YANG Weitun
Key Laboratory of Technology in Geo-spatial Information Processing and Application System of CAS, Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China; University of Chinese Academy of Sciences, Beijing 100190, China


Abstract: Independent component analysis (ICA) is a popular signal processing method based on the high-order statistical characteristics of data. It has been widely used in imagery processing. However, classical FastICA algorithm achieves independent components (ICs) of the data by the fixed-point iteration method. It needs cyclic iterations to get every IC, and in the iterative process the solution of each IC requires the participation of all pixels. Therefore, when the amount of data or the number of pixels is large, FastICA is time-consuming. At this point, its speed advantage is greatly reduced. The remote sensing data generally has a large size. So how to directly apply FastICA to remote sensing data is a practical problem. In this study, by introducing the cokurtosis tensor, the fixed point iteration problem of FastICA is transformed into tensor computation in algebraic form, and the participation of all pixels is avoided, thus greatly reducing the computational complexity. Experimental results for multispectral images show that the proposed algorithm is faster than the classical FastICA algorithm based on kurtosis.
Keywords: FastICAhigh order statisticskurtosis indexcokurtosis tensor
多/高光谱遥感图像通常包含多个光谱波段^[1]，可提供关于地物的丰富的光谱信息，但同时也会造成一定的数据冗余及数据处理过程的计算负担^[2]。因此，在很多情况下需要引入数据降维或者特征提取的手段对数据进行预处理^[3]。主成分分析^[4-5](principal component analysis，PCA)是一种最常用的特征提取方法，并常被应用于遥感图像的数据降维。PCA以图像方差为指标进行特征提取，而没有考虑噪声在光谱特征空间的分布，因此降维结果易受具有较大方差的噪声影响。Green等^[6]提出基于图像信噪比的特征提取方法，即最大噪声成分(maximum noise fraction，MNF)算法，以解决这一问题。Lee等^[7]提出噪声调整的主成分分析方法(NAPC)，也能够减少噪声对特征提取结果的影响，它基本与MNF等价。此外，还有考虑了图像空间特征的最大/最小自相关因子分析算法，它基于图像数据的自相关性而不是方差，因此不受向量幅值的影响。
以上这些特征提取算法都是基于数据的二阶统计特性，它一般假设图像数据服从高斯分布。但实际的遥感图像数据往往并不满足这一条件，所以基于二阶统计特性的特征提取方法很多情况下并不能有效反映数据中的感兴趣结构。比如，当感兴趣目标为小目标时，基于二阶统计特性的特征提取方法往往将其作为不重要的信息舍弃。而基于高阶统计特性的特征提取方法就能很好地解决这类问题^[8-9]。高阶统计特性在遥感图像处理中有很多应用，例如目标自动识别、异常检测、光谱解混、特征提取等^[10-12]。基于高阶统计特性的特征提取方法绝大多数是基于独立成分分析(independent component analysis，ICA)的方法及其改进算法^[13-15]，其中FastICA算法以其快速的收敛速度和较高的计算效率得到广泛应用。FastICA算法有4种优化指标，其中基于峭度指标的FastICA算法具有明确的物理意义，且相比基于偏度指标的FastICA算法有更好的收敛性能和更快的收敛速度，因此应用最为广泛^[16-20]。但是，在求解各个独立成分时每一步迭代过程均需所有像元的参与，所以当图像尺寸较大时往往计算量很大且耗时较长。然而，遥感图像尺寸一般较大，此时FastICA算法的计算优势大打折扣。
因此，本文提出一种基于峭度指标的FastICA算法的等价算法。通过引入协峭度张量，将FastICA的固定点迭代问题转化为代数形式的张量计算，在每次迭代中避免所有像元的参与，因而大大降低计算复杂度。
1 基于峭度指标的FastICA算法ICA算法是一种从盲源分离技术发展而来的算法^[15]，它能够将多维观测信号分解为尽量相互统计独立的独立分量，分离出的分量是源信号的一种近似，因而在图像处理、语音识别、远程通信等领域得到广泛应用。
为了更为方便和简洁地说明ICA，假设观测信号是源信号的线性变换。假设观测信号为

$\mathit{\boldsymbol{X}} = \left( {\begin{array}{*{20}{c}}{{x_{11}}}& \cdots &{{x_{1n}}}\\ \vdots &{}& \vdots \\{{x_{p1}}}& \cdots &{{x_{pn}}}\end{array}} \right) = \left( {\begin{array}{*{20}{c}}{{{\mathit{\boldsymbol{\underline x}} }_1}}\\{{{\mathit{\boldsymbol{\underline x}} }_2}}\\ \vdots \\{{{\mathit{\boldsymbol{\underline x}} }_\mathit{\boldsymbol{p}}}}\end{array}} \right) = \left[ {{\mathit{\boldsymbol{x}}_1},{\mathit{\boldsymbol{x}}_2}, \cdots ,{\mathit{\boldsymbol{x}}_n}} \right],$

(1)

它由相互独立且均值为零的源信号

$\mathit{\boldsymbol{S}} = \left( {\begin{array}{*{20}{c}}{{s_{11}}}& \cdots &{{s_{1n}}}\\ \vdots &{}& \vdots \\{{s_{p1}}}& \cdots &{{s_{pn}}}\end{array}} \right) = \left( {\begin{array}{*{20}{c}}{{{\mathit{\boldsymbol{\underline s}} }_1}}\\{{{\mathit{\boldsymbol{\underline s}} }_2}}\\ \vdots \\{{{\mathit{\boldsymbol{\underline s}} }_\mathit{\boldsymbol{p}}}}\end{array}} \right)$

(2)

线性混合而成，其中$\underline{\boldsymbol{x}}_{1}, \underline{\boldsymbol{x}}_{2}, \cdots, \underline{\boldsymbol{x}}_{p}$代表X中的行向量，$\underline{\boldsymbol{s}}_{1}, \underline{\boldsymbol{s}}_{2}, \cdots, \underline{\boldsymbol{s}}_{p}$代表S中的行向量。ICA要解决的问题是根据观测信号得到源信号，即找到一个方阵W = [w₁, w₂, …, w_p]满足如下等式

$\mathit{\boldsymbol{S}} = {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{X}} = \left( {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{w}}_1^{\rm{T}}\mathit{\boldsymbol{X}}}\\{\mathit{\boldsymbol{w}}_2^{\rm{T}}\mathit{\boldsymbol{X}}}\\ \cdots \\{\mathit{\boldsymbol{w}}_{\rm{p}}^{\rm{T}}\mathit{\boldsymbol{X}}}\end{array}} \right).$

(3)

ICA可以描述为在源信号和混合矩阵都未知的情况下，根据观测矩阵X找到一个分离矩阵W能够尽可能地分离出源信号S。FastICA是ICA算法中最为流行，应用最广泛的一种算法。相比于其他算法，FastICA算法计算效率高，运算量小，相对容易实现。
对于零均值的随机变量x，其峭度可以表示为

${\mathop{\rm kurt}\nolimits} (x) = E\left\{ {{x^4}} \right\} - 3{\left( {E\left\{ {{x^2}} \right\}} \right)^2},$

(4)

显然，对于高斯信号，它的峭度值为0。
假设观测信号X是源信号S线性混合而成，经过向量w可以分离出一个源信号分量w^TX。一般情况下，在处理观测信号之前会进行数据白化，因此，不妨假设，观测信号X的均值向量为0，协方差矩阵为单位阵，令y = w^TX，则它的峭度计算可以进行如下简化^[20]

$\begin{aligned} \operatorname{kurt}(\boldsymbol{y})=& E\left\{\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{X}\right)^{4}\right\}-3\left(E\left\{\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{X}\right)^{2}\right\}\right)^{2} \\ &=E\left\{\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{X}\right)^{4}\right\}-3\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{X} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{w}\right)^{2} \\ &=E\left\{\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{X}\right)^{4}\right\}-3\|\boldsymbol{w}\|^{4}, \end{aligned}$

(5)

式中:w满足约束条件‖w‖ = 1，基于峭度指标的FastICA算法的目标函数可以最终表示为

$J(\boldsymbol{w})=E\left\{\left(\boldsymbol{w}^{\mathrm{T}} \boldsymbol{X}\right)^{4}\right\}-3\|\boldsymbol{w}\|^{4}+F\left(\|\boldsymbol{w}\|^{2}\right),$

(6)

式中:F是关于约束条件的惩罚函数。固定点算法求解式(6)可以得到迭代公式

$\left\{\begin{aligned} \boldsymbol{w}^{(k+1)} &=\frac{\boldsymbol{X}\left[\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{w}^{(k)}\right) \odot\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{w}^{(k)}\right) \odot\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{w}^{(k)}\right)\right]}{n}-3 \boldsymbol{w}^{(k)} \\ \boldsymbol{w}^{(k+1)} &=\frac{\boldsymbol{w}^{(k+1)}}{\left\|\boldsymbol{w}^{(k+1)}\right\|} \end{aligned}\right.,$

(7)

式中：“⊙”符号代表矩阵点乘，表示2个矩阵对应位置的每个元素相乘。下同。
综合上文，基于峭度指标的FastICA算法步骤表示如下
1) 首先随机生成一个模为1的向量w⁽⁰⁾，并令k = 0；
2) 将其代入迭代公式

$\begin{array}{l}{\mathit{\boldsymbol{w}}^{(k + 1)}} = \\\frac{{\mathit{\boldsymbol{X}}\left[ {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{w}}^{(k)}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{w}}^{(k)}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}{\mathit{\boldsymbol{w}}^{(k)}}} \right)} \right]}}{n} - 3{\mathit{\boldsymbol{w}}^{(k)}};\end{array}$

3)将w^(k+1)单位化${\boldsymbol{w}}^{(k+1)}=\frac{{\boldsymbol{w}}^{(k+1)}}{\left\|{\boldsymbol{w}}^{(k+1)}\right\|}$；
4) 如果|(w^(k+1))^Tw^(k)|的值不足够接近1，则令k = k+1，跳到步骤2)；否则输出向量w^(k+1)。
上述过程得到的向量w^(k+1)为分离矩阵W的一个列向量，(w^(k+1))^TX可以分离出一个源信号的近似，即一个独立分量。为了分离出p个独立分量，需要将上述步骤运行p次。同时为保证每次分离出的独立分量是不同的，要对w^(k+1)进行正交投影，使其与已求得的其他向量均正交，可通过下式进行正交化

$\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{w}}^{(k + 1)}} = {\mathit{\boldsymbol{w}}^{(k + 1)}} - {{\overline {\mathit{\boldsymbol{WW}}} }^{\rm{T}}}{\mathit{\boldsymbol{w}}^{(k + 1)}}}\\{ = \left( {\mathit{\boldsymbol{I}} - {{\overline {\mathit{\boldsymbol{WW}}} }^{\rm{T}}}} \right){\mathit{\boldsymbol{w}}^{(k + 1)}},}\end{array}$

(8)

式中：W即为已经求出的w向量的集合。
上述算法具有较快的收敛速度，一般情况下，它在10~100次迭代甚至几次迭代内就会得到最优解。但它的运算量依然很大，而且循环迭代求解过程需要全部数据的参与，因此下面提出一种改进的FastICA算法以减少运算量，新算法引入协峭度张量的概念使得循环迭代过程中不需要所有像元的参与。
2 改进的FastICA算法改进的FastICA算法是对原始FastICA算法的迭代公式的改进，通过引入协峭度张量，在降低算法计算复杂度的同时避免所有像元参与算法的迭代过程。下面，首先介绍协峭度张量的概念。
2.1 协峭度张量的计算仍然假设白化后的数据为$\boldsymbol{X = }\left(\begin{matrix} {{x}_{11}} & \cdots & {{x}_{1n}} \\ \vdots & {} & \vdots \\ {{x}_{p1}} & \cdots & {{x}_{pn}} \\\end{matrix} \right) = \left[{{\boldsymbol{x}}_{\text{1}}}, {{\boldsymbol{x}}_{\text{2}}}, \cdots {{\boldsymbol{x}}_{n}} \right]$, 那么协峭度张量可通过下式计算得到

$\mathit{\boldsymbol{\underline K}} = \frac{1}{n}\sum\limits_{i = 1}^n {{\mathit{\boldsymbol{x}}_i}} \circ {\mathit{\boldsymbol{x}}_i} \circ {\mathit{\boldsymbol{x}}_i} \circ {\mathit{\boldsymbol{x}}_i}.$

(9)

式中:“°”代表外积，x_i°x_i°x_i°x_i可以生成一个p维的4阶张量。根据式(9)可知，

$\mathit{\boldsymbol{\underline K}} (i,j,s,t) = \frac{1}{n}\sum\limits_{k = 1}^n {{x_{ik}}} {x_{jk}}{x_{sk}}{x_{tk}},$

(10)

显然K是一个p×p×p×p的超对称张量。
2.2 改进的FastICA算法原理定理2.1 ?对于任意向量w = [w₁, …, w_p]^T，有下面的结论成立：

$\begin{array}{*{20}{c}}{\mathit{\boldsymbol{X}}\left[ {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right)} \right]/n}\\{ = \mathit{\boldsymbol{\underline K}} { \times _1}\mathit{\boldsymbol{w}}{ \times _2}\mathit{\boldsymbol{w}}{ \times _3}\mathit{\boldsymbol{w}}.}\end{array}$

(11)

则式(7)的迭代公式变为

$\left\{ \begin{array}{l}{\mathit{\boldsymbol{w}}^{(k + 1)}} = \mathit{\boldsymbol{\underline K}} { \times _1}{\mathit{\boldsymbol{w}}^{(k)}}{ \times _2}{\mathit{\boldsymbol{w}}^{(k)}}{ \times _3}{\mathit{\boldsymbol{w}}^{(k)}} - 3{\mathit{\boldsymbol{w}}^{(k)}}\\{\mathit{\boldsymbol{w}}^{(k + 1)}} = \frac{{{\mathit{\boldsymbol{w}}^{(k + 1)}}}}{{\left\| {{\mathit{\boldsymbol{w}}^{(k + 1)}}} \right\|}}\end{array} \right..$

(12)

式中：×₁×₂×₃代表模乘，下同。
对于大小为p×n的数据，X[(X^Tw)⊙(X^Tw)⊙(X^Tw)]/n的乘法运算量为2n(p+1)+p次，而K×₁w×₂w×₃w的乘法运算量为(p⁴+p³+p²)。一般情况下n$\gg $p，所以改进的FastICA算法将大大降低运算复杂度。将基于峭度指标的改进的FastICA算法应用于多/高光谱图像特征提取的过程可以描述如下：
Table

改进的FastICA算法伪代码 ??1)首先，计算图像的协峭度张量 $ \underline{\boldsymbol{K}}=\frac{1}{n} \sum\limits_{i=1}^{n} \boldsymbol{x}_{i} \circ \boldsymbol{x}_{i} \circ \boldsymbol{x}_{i} \circ \boldsymbol{x}_{i}$； ??2)令 P =I ，P 为正交投影矩阵; ??3)for i=1 :取的独立成分数； ??4)令k=0，并生成随机单位向量 wi(k)； ??5)将 wi(k)进行正交投影 wi(k)= Pwi(k)； ??6)将 wi(k)进行迭代更新 $\left\{\begin{aligned} \boldsymbol{w}_{i}^{(k+1)} &=\underline{\boldsymbol{K}} \times_{1} \boldsymbol{w}_{i}^{(k)} \times_{2} \boldsymbol{w}_{i}^{(k)} \times_{3} \boldsymbol{w}_{i}^{(k)}-3 \boldsymbol{w}_{i}^{(k)} \\ \boldsymbol{w}_{i}^{(k+1)} &=\frac{\boldsymbol{w}_{i}^{(k+1)}}{\left\|\boldsymbol{w}_{i}^{(k+1)}\right\|} \end{aligned}\right.;$ ??7)k=k+1； ??8)如果满足停止条件，则输出 wi(k+1)，否则跳到5)； ??9)令 P = I - WWT，W 为已求得的向量 wi 的集合 W =[ w1, …,wi ]； ??10)end； ?? Y =WTX即为特征提取后的结果。

步骤8)中的停止条件有2个：一是容错率(τ)，即要求w_i^(k+1)和w_i^(k)足够接近(如1－|w_i^(k+1)w_i^(k)|≤τ)；二是迭代次数，即设置一个最大迭代次数k_max，当到达最大迭代次数时跳出循环，以防止程序进入死循环。
2.3 算法复杂度分析比较基于峭度指标的FastICA算法及其改进算法的迭代公式可知，改进的FastICA算法的迭代公式更加简洁，且不需要全部像元点的参与，因此会节省更多的计算时间。为了更加准确地比较这两种算法，需要对它们进行详细的计算复杂度分析。
假设需要从p个波段的多/高光谱图像数据中提取出p个独立分量，且该图像包含n个像元点。由于改进后的FastICA算法与原始的FastICA算法相比仅有两处不同，一是迭代公式不同，二是是否需要计算协峭度张量。因此，只针对这两处不同进行运算量的比较。FastICA算法中的X[(X^Tw)⊙(X^Tw)⊙(X^Tw)]/n的乘法运算量为(2n(p+1)+p)，改进后的FastICA算法中K×₁w×₂w×₃w的乘法运算量为(p⁴+p³+p²)，如果提取一个独立分量平均需要k次迭代，则迭代过程需要运行kp次，则原始的FastICA算法相应部分需要kp[2np+2n+p]≈2knp²次乘法运算；而改进FastICA算法需要kp(p⁴+p³+p²)≈kp⁵次乘法运算，但它需要提前计算协峭度张量(np⁴/8次乘法运算)，所以它总共需要约(kp⁵+np⁴/8)次乘法运算。由于n$\gg $kp，所以改进的FastICA算法的计算复杂度可以简化为np⁴/8，那么两种算法的计算复杂度比值约为16k/p²(FastICA/改进的FastICA)。这说明当平均需要的迭代次数较多而波段数较少时，改进的FastICA算法会节省计算时间；反之，当波段数较大而需要的平均迭代次数较少时，FastICA算法较节省时间。
3 基于模拟数据的实验分析原始的FastICA算法和改进后的FastICA算法虽然都使用固定点迭代算法，但是二者有本质上的不同。原始的FastICA算法在迭代过程中需要全部像元的参与，改进后的FastICA算法则依赖预先计算好的协峭度张量。为了更加客观地比较两种算法，FastICA算法和改进的FastICA算法使用相同的初始向量。所有的仿真实验均在同一电脑上使用相同的操作系统完成。本文中的算法全部使用MATLAB7.1软件实现。
3.1 计算时间和数据维数之间的关系在本节中，设计一组实验考察数据的维数对两种算法计算时间的影响。首先用MATLAB软件中的rand函数产生10组模拟数据，每组数据尺寸为800×800，维数从2增加到15。考虑到初始值的影响，对两种算法设置相同的初始值，各运行100次取平均计算时间进行比较。图 1给出两种算法的计算时间。从图中可以看出随着数据维数的增加，改进的FastICA需要的计算时间急剧增加，当数据维数超过11时，改进的FastICA算法反而需要更多的计算时间。这也和上文分析的结果一致。图 2给出计算协峭度张量的时间占总的计算时间的比重，可以看出随着数据维数的增加它所占的比重越来越大，甚至占用绝大部分的计算时间。因此，如果能够准确快速地估计数据的协峭度张量将极大地提高改进的FastICA算法的计算效率。
Fig. 1

	Download: JPG larger image
图 1 原始FastICA算法和改进的FastICA算法的计算时间随波段数变化的趋势 Fig. 1 Variations in the computing time of the original and improved FastICA algorithms with the band number 图 1 原始FastICA算法和改进的FastICA算法的计算时间随波段数变化的趋势 Fig. 1 Variations in the computing time of the original and improved FastICA algorithms with the band number -->

Fig. 2

	Download: JPG larger image
图 2 不同数据维数条件下计算协峭度张量的时间占总的计算时间的比重 Fig. 2 Ratio of the calculation time for cokurtosis tensor to the total calculation time 图 2 不同数据维数条件下计算协峭度张量的时间占总的计算时间的比重 Fig. 2 Ratio of the calculation time for cokurtosis tensor to the total calculation time -->

3.2 计算时间和像元数目之间的关系为考察像元总的数目对两种算法的影响，生成一组维数为6，像元数目从100×100到900×900的随机模拟数据。同样地，对两种算法设置相同的初始值，统计其所需要的计算时间(100次平均)，如图 3所示。从图中可以看出随着像元数目的增多，两种算法所需的计算时间都在逐渐增加，但改进后的FastICA算法增加得更慢，且它需要的计算时间始终少于原始FastICA算法需要的计算时间。
Fig. 3

	Download: JPG larger image
图 3 原始FastICA算法和改进的FastICA算法在不同像元数目条件下的计算时间比较 Fig. 3 Comparison of the computing time between the original FastICA and improved FastICA at different pixel numbers 图 3 原始FastICA算法和改进的FastICA算法在不同像元数目条件下的计算时间比较 Fig. 3 Comparison of the computing time between the original FastICA and improved FastICA at different pixel numbers -->

4 基于真实数据的实验分析本节使用真实的多光谱图像数据比较FastICA算法及改进FastICA算法的计算效率。该多光谱数据为Landsat5 TM反射率数据，获取时间为2001年6月30日，该数据为江西北部的鄱阳湖部分区域数据，如图 4所示。该图像分辨率为30 m，尺寸为800×800，共包含6个波段，分别为485，569，660，840，1 676和2 223 nm。
Fig. 4

	Download: JPG larger image
图 4 鄱阳湖第一波段图像 Fig. 4 The first-band image of Poyang Lake 图 4 鄱阳湖第一波段图像 Fig. 4 The first-band image of Poyang Lake -->

将原始的基于峭度指标的FastICA算法和改进后的FastICA算法应用于该多光谱图像，将特征提取结果进行详细对比分析发现，原始的基于峭度指标的FastICA算法和改进后的FastICA算法提取结果一致，特征提取结果如图 5所示。同时基于峭度指标的FastICA算法需要的时间为8.529 2 s，而改进后的FastICA算法仅仅需要2.275 5 s。因此，改进后的FastICA算法在保持原有结果不变的同时计算速度远快于原始的FastICA算法。
Fig. 5

	Download: JPG larger image
图 5 两种FastICA算法提取结果 Fig. 5 The feature extraction results of the original FastICA and improved FastICA 图 5 两种FastICA算法提取结果 Fig. 5 The feature extraction results of the original FastICA and improved FastICA -->

5 总结本文提出一种改进的FastICA算法以降低FastICA的计算复杂度，它是原始FasICA算法的代数等价算法，能够产生与之完全相同的特征提取结果。改进的FastICA算法引入协峭度张量的概念，在迭代过程中不需要全部像元的参与，而仅仅需要协峭度张量的参与。当数据维数不是特别高时，改进的FastICA算法相比于原始的FastICA算法计算复杂度更低，能节省更多的计算时间。但当数据维数很高时，改进的FastICA算法实际上并不能节省计算时间。因此该方法更加适用于多光谱图像的特征提取。改进的FastICA算法中计算协峭度张量花费绝大部分时间，因此寻求一种快速计算协峭度张量的方法可以极大地提高算法的计算效率，缩短计算时间。
附录：定理2.1的证明
假设白化后的数据为$\boldsymbol{X = }\left(\begin{matrix} {{x}_{11}} & \cdots & {{x}_{1n}} \\ \vdots & {} & \vdots \\ {{x}_{p1}} & \cdots & {{x}_{pn}} \\\end{matrix} \right) = \left[{{\boldsymbol{x}}_{\text{1}}}, {{\boldsymbol{x}}_{\text{2}}}, \cdots {{\boldsymbol{x}}_{n}} \right]$，向量w = [w₁, …, w_p]^T.
首先计算X[(X^Tw)⊙(X^Tw)⊙(X^Tw)]/n.
令

$\mathit{\boldsymbol{b}} = \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) = {\left[ {{{\left( {\sum\limits_{i = 1}^p {{x_{i1}}} {w_i}} \right)}^3}, \cdots ,{{\left( {\sum\limits_{i = 1}^p {{x_{in}}} {w_i}} \right)}^3}} \right]^{\rm{T}}},$

(A1)

式(A1)中第k个元素表示为

$b(k) = {\left( {\sum\limits_{i = 1}^p {{x_{ik}}} {w_i}} \right)^3} = \sum\limits_{i = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{s = 1}^p {{x_{ik}}} } } {w_i}{x_{jk}}{w_j}{x_{sk}}{w_s},$

(A2)

则

$\begin{array}{*{20}{l}}{\mathit{\boldsymbol{X}}\left[ {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right)} \right] = \mathit{\boldsymbol{Xb}} = {{\left[ {\sum\limits_{k = 1}^n {{x_{1k}}} b(k), \cdots ,\sum\limits_{k = 1}^n {{x_{pk}}} b(k)} \right]}^{\rm{T}}} = }\\{{{\left[ {\sum\limits_{k = 1}^n {{x_{1k}}} \sum\limits_{i = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{s = 1}^p {{x_{ik}}} } } {w_i}{x_{jk}}{w_j}{x_{sk}}{w_s}, \cdots ,\sum\limits_{k = 1}^n {{x_{pk}}} \sum\limits_{i = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{s = 1}^p {{x_{ik}}} } } {w_i}{x_{jk}}{w_j}{x_{sk}}{w_s}} \right]}^{\rm{T}}},}\end{array}$

(A3)

即

$\begin{array}{*{20}{c}}{\mathit{\boldsymbol{X}}\left[ {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right) \odot \left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{w}}} \right)} \right]/n = }\\{\frac{1}{n}{{\left[ {\sum\limits_{k = 1}^n {\sum\limits_{i = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{j = 1}^p {{x_{ik}}} } } } {w_i}{x_{jk}}{w_j}{x_{sk}}{w_s}{x_{1k}}, \cdots ,\sum\limits_{k = 1}^n {\sum\limits_{i = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{s = 1}^p {{x_{ik}}} } } } {w_i}{x_{jk}}{w_j}{x_{sk}}{w_s}{x_{pk}}} \right]}^{\rm{T}}}.}\end{array}$

(A4)

要计算K×₁w×₂w×₃w，首先要计算协峭度张量K，已知

$\underline K (i,j,s,t) = \frac{1}{n}\sum\limits_{k = 1}^n {{x_{ik}}} {x_{jk}}{x_{sk}}{x_{tk}},\mathit{\boldsymbol{\underline K}} \;为超对称张量.$

(A5)

先计算K?w，令S=K×₁w，
则S的任一元素为

$S(j,s,t) = \sum\limits_{i = 1}^p {\underline K } (i,j,s,t){w_i},$

(A6)

令D = K×₁w×₂w = S×₂w为p×p的矩阵，则D的任一元素为

$D(s,t) = \sum\limits_{j = 1}^p S (j,s,t){w_j} = \sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {\underline K } } (i,j,s,t){w_i}{w_j}.$

(A7)

令y = K×₁w×₂w×₃w，则y = D×₃w为p×1的向量，其任一元素为

$y(t) = \sum\limits_{s = 1}^p D (s,t){w_s} = \sum\limits_{s = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {\underline K } } } (i,j,s,t){w_i}{w_j}{w_s}.$

(A8)

又$\underline{K}(i, j, s, t) = \frac{1}{n} \sum\limits_{k = 1}^{n} x_{i k} x_{j k} x_{s k} x_{t k}$，则

$y(t) = \sum\limits_{s = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {\frac{1}{N}} } } \sum\limits_{k = 1}^n {{x_{ik}}} {x_{jk}}{x_{sk}}{x_{tk}}{w_i}{w_j}{w_s} = \frac{1}{n}\sum\limits_{s = 1}^p {\sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {\sum\limits_{k = 1}^n {{x_{ik}}} } } } {x_{jk}}{x_{sk}}{x_{tk}}{w_i}{w_j}{w_s}.$

(A9)

那么，

$\mathit{\boldsymbol{\underline K}} { \times _1}\mathit{\boldsymbol{w}}{ \times _2}\mathit{\boldsymbol{w}}{ \times _3}\mathit{\boldsymbol{w}} = \frac{1}{n}{\left[ {\sum\limits_{i,j,s = 1}^p {\sum\limits_{k = 1}^n {{w_i}} } {w_j}{w_s}{x_{ik}}{x_{jk}}{x_{sk}}{x_{1k}}, \cdots ,\sum\limits_{i,j,s = 1}^p {\sum\limits_{k = 1}^n {{w_i}} } {w_j}{w_s}{x_{ik}}{x_{jk}}{x_{sk}}{x_{pk}}} \right]^{\rm{T}}}.$

(A10)

比较式(A4)和式(A10)可知K×₁w×₂w×₃w与X[(X^Tw)⊙(X^Tw)⊙(X^Tw)]/n的结果是相等的。
参考文献

[1]	Castaldi F, Palombo A, Santini F, et al. Evaluation of the potential of the current and forthcoming multispectral and hyperspectral imagers to estimate soil texture and organic carbon[J]. Remote Sensing of Environment, 2016, 179: 54-65.
[2]	Tian Y, Guo P, Lyu M R. Comparative studies on feature extraction methods for multispectral remote sensing image classification[C]//Systems Man and Cybernetics. Waikoloa HI: IEEE, 2006: 1275-1279.
[3]	Kaur R, Sehgal S. Dimension reduction of multispectral data using canonical analysis[J]. International Journal of Computer Applications, 2014, 70(21): 18-21.
[4]	Jaime Z, Ren J, Ren J. Structured covariance principal component analysis for real-time onsite feature extraction and dimensionality reduction in hyperspectral imaging[J]. Applied Optics, 2014, 53(20): 4440-4449. DOI:10.1364/AO.53.004440
[5]	Wold S, Esbensen K, Geladi P. Principal component analysis[J]. Chemometrics and Intelligent Laboratory Systems, 1987, 2(1): 37-52.
[6]	Green A A, Berman M, Switzer P, et al. A transformation for ordering multispectral data in terms of image quality with implications for noise removal[J]. IEEE Transactions on Geoscience and Remote Sensing, 1988, 26(1): 65-74. DOI:10.1109/36.3001
[7]	Lee J B, Woodyatt A S, Berman M. Enhancement of high spectral resolution remote-sensing data by a noise-adjusted principal components transform[J]. IEEE Transactions on Geoscience and Remote Sensing, 1990, 28(3): 295-304. DOI:10.1109/36.54356
[8]	Wang C, Zhang J, Gu Y. Target detection for hyperspectral images using ICA-based feature extraction[C]//Geoscience and Remote Sensing Symposium. Denver CO: IEEE, 2006: 850-853.
[9]	Ren H, Du Q, Wang J, et al. Automatic target recognition for hyperspectral imagery using high-order statistics[J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42(4): 1372-1385.
[10]	Gu Y, Liu Y, Zhang Y. A selective KPCA algorithm based on high-order statistics for anomaly detection in hyperspectral imagery[J]. IEEE Geoscience and Remote Sensing Letters, 2008, 5(1): 43-47. DOI:10.1109/LGRS.2007.907304
[11]	Wang N, Du B, Zhang L. An abundance caracteristic-based independent component analysis for hyperspectralunmixing[J]. IEEE Geoscience and Remote Sensing, 2015, 3(1): 416-428.
[12]	Wang J, Chang C. Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery[J]. International Society for Optical Engineering, 2006, 44(9): 2601-2616.
[13]	Lennon M, Mercier G, Mouchot MC, et al. Independent component analysis as a tool for the dimensionality reduction and the representation of hyperspectral images[C]//Geoscience and Remote Sensing Symposium. Sydney NSW: IEEE, 2001, 6(3): 2893-2895.
[14]	Hyvarinen A, Oja E, Hoyer P, et al. Image feature extraction by sparse coding and Independent component analysis[C]//Pattern Recognition. Brisbane: IEEE, 1998, 2: 1268-1273.
[15]	Comon P. Independent component analysis, a new concept?[J]. Signal Processing, 1994, 36(3): 287-314.
[16]	Chen P, Hung H, Komori O, et al. Robust independent component analysis via minimum-divergence estimation[J]. IEEE Journal of Selected Topics in Signal Processing, 2013, 7(4): 614-624. DOI:10.1109/JSTSP.2013.2247024
[17]	Shen H, Gretton A. Fast kernel-based independent component analysis[J]. IEEE Transactions on Signal Processing, 2009, 57(9): 3498-3511. DOI:10.1109/TSP.2009.2022857
[18]	Karvanen J, Koivunen V. Independent component analysis via optimum combining of kurtosis and skewness-based criteria[J]. Journal of the Franklin Institute, 2004, 341(5): 401-418. DOI:10.1016/j.jfranklin.2004.04.001
[19]	Hyvarinen A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE Transactions on Neural Networks, 1999, 10(3): 626-634. DOI:10.1109/72.761722
[20]	Hyv?rinen A, K?ster U. A fast fixed-point algorithm for independent component analysis[J]. Neural Computation, 1997, 9(7): 1483-1492. DOI:10.1162/neco.1997.9.7.1483