1.   引言

随着遥感技术的发展，人们可以获得越来越多的高维大规模遥感图像数据。在过去的几十年中，许多学者使用遥感图像数据在光谱融合^[1]、目标分类^[2]、目标检测^[3]、图像去噪^[4]、图像压缩^[5]等研究领域做了大量研究。高维大规模数据所蕴含的巨大的信息为后续遥感图像处理带来巨大挑战，过高的数据维度不利于后续的存储、传输以及运算，而特征提取成为解决这一问题的关键技术。

近几年研究人员提出了多种用于特征提取的算法，典型方法包括主成分分析方法(Principal Component Analysis, PCA)^[6]、奇异值分解(Singular Value Decomposition, SVD) 和独立成分分析方法(Independent Component Analysis, ICA)^[7]。PCA是基于数据的2阶相关统计特征找到一组正交向量集，这个向量集从最小二乘意义上对原始信号进行表示。SVD则是在最大能量的意义上选择较大的变量对原矩阵进行分解，以此来提取数据矩阵的主要特征。上述方法主要针对数据的2阶统计特征，并假设数据满足高斯分布。然而，许多真实数据集的分布通常不满足高斯分布。在这种情况下，基于高阶统计量的ICA方法近年来受到了越来越多的关注。ICA是对数据进行高阶统计分析以寻求非高斯量最大化来分离各个独立的信号源。FastICA是ICA的一种改进算法，它可以选择偏度、负熵、峭度等指标作为非高斯度量^[8]。虽然FastICA算法能达到3次收敛速度，优于其他大多数常用的ICA算法^[9]，但是其每次迭代都需要涉及所有的像素来寻找最优的投影方向，这导致其在处理高维数据时耗时大。

为了解决这一问题，Geng等人^[10]提出了一种主偏度分析方法(Principal Skewness Analysis, PSA)，该方法的求解可以转化为计算数据3阶统计张量特征对，因此PSA可以看作PCA的3阶泛化。当选择偏度作为非高斯指标时，可等效为FastICA。在此基础上，文献[11,12]进一步提出了动量PSA (Momentum Principal Skewness Analysis, MPSA) 和主峰度分析方法(Principal Kurtosis Analysis, PKA)。超对称张量并不具有天然的正交性^[13,14]，然而PSA和MPSA算法在求解过程中均采用了正交补策略，因此它们得到的解(特征对)是相互正交的。这导致除了第1个特征对外，得到的其他特征对都偏离了协偏度张量的实际特征对。针对这个问题，Geng等人^[15]又提出了一种新的非正交的PSA (Nonorthogonal Principal Skewness Analysis, NPSA)算法。NPSA算法通过引入克罗内克积提出了一种新的搜索策略，该策略可以在更大的空间内搜索解，因此可以获得比PSA算法更精确的解，但是NPSA算法得到的仍然是近似解。

针对上述局限，本文提出了一种新的高精度并行主偏度分析(Parallel Principle of Skewness Analysis, PPSA)算法，PPSA算法把协偏度张量切片的特征向量作为初始值进行迭代求解。与现有PSA算法施加约束再依次求解各特征对不同，PPSA算法不再施加约束，而是将各个特征对的初始值并行地、一次性地全部算出，有效地避免了施加约束带来的误差。相比以前只能得到近似解，PPSA算法可以准确得到实际解。

2.   PSA算法

在PSA算法中，构造协偏度张量，类似于构造PCA中的协方差矩阵，来计算图像在${\boldsymbol{u}}$方向上的偏度。本文假设图像数据集为${\boldsymbol{X}} = [{\boldsymbol{x}}{_1},{{\boldsymbol{x}}_i}_2, \cdots ,{{\boldsymbol{x}}_i}_N] \in {{\boldsymbol{R}}^{L \times N}}$，其中${{\boldsymbol{x}}_i}$是$L \times 1$向量，$N$是像素个数。图像首先应集中并满足${\boldsymbol{R}} = {{\boldsymbol{F}}^{\rm{T}}}({\boldsymbol{X}} - {\boldsymbol{m}})$，其中${\boldsymbol{m}} = (1/N)\displaystyle\sum\nolimits_{i = 1}^N {{x_i}}$是均值，${\boldsymbol{R}} = [{r_1},{r_2}, \cdots ,{r_N}]$为白化图像。${\boldsymbol{F}} = {\boldsymbol{E}}{{\boldsymbol{D}}^{ - 0.5}}$称为白化算子，其中${\boldsymbol{E}}$为协方差矩阵的特征向量矩阵，${\boldsymbol{D}}$为对应的特征值对角矩阵。然后，协偏度张量由式(1)计算得到。

其中，“°”表示两个向量的外积，协偏度张量${\boldsymbol{S}}$是一个大小为$L \times L \times L$的超对称张量。那么图像在任意方向上的偏度${\boldsymbol{u}}$可以通过式(2)计算。

${\boldsymbol{u}} \in {{\boldsymbol{R}}^{L \times 1}}$是一个单位向量, ${{\boldsymbol{u}}^{\rm{T}}}{\boldsymbol{u}} = 1$。因此优化模型为

用拉格朗日方法求解上述方程得到

采用定点法计算每个单位的${\boldsymbol{u}}$, ${\boldsymbol{u}}$可以表示为

如果它存在一个不动点，则称解${\boldsymbol{u}}$为第1主偏度方向，$\lambda$为其偏度。

同样，$(\lambda ,{\boldsymbol{u}})$也称为一个张量的特征值/特征向量对，由Lim^[16]和Qi^[17]引入。

防止第2特征向量收敛于相同的特征向量，在第1个向量的正交补空间生成新的张量计算${\boldsymbol{u}}$

其中，${\boldsymbol{P}}_u^ \bot = {\boldsymbol{I}} - {\boldsymbol{u}}{({{\boldsymbol{u}}^{\rm{T}}}{\boldsymbol{u}})^{ - 1}}{{\boldsymbol{u}}^{\rm{T}}}$是${\boldsymbol{u}}$的正交补投影算子，${\boldsymbol{I}}$为$L \times L$单位矩阵。然后，对新张量${\boldsymbol{S}}$采用同样的迭代方法，得到第2个特征向量，接下来的过程也采用同样的方法进行，直至求解得到$L$个特征对。

3.   PPSA算法

3.1   现有PSA算法局限

现有的 PSA 算法主要采用了串行方式依次求解各个特征对，为了避免收敛到同一特征对而施加约束，但是这也将导致除了第1个特征对外，其他特征对都偏离了实际特征对^[15]。并且现有PSA算法初始化均采用随机生成初始值，这意味着给定不同的初始化，可能收敛到不同的解(特征对)^[18]。在这里举一个简单的例子来更直观地说明这两种现象。考虑一个2阶超对称张量${\boldsymbol{S}} \in {\mathbb{R}^{2 \times 2 \times 2}}$，其两个切片分别为

其解为${{\boldsymbol{u}}_1} = {\left[ {0.8812, - 0.4727} \right]^{\rm{T}}},{{\boldsymbol{u}}_2} = [0.3757, 0.9267]^{\rm{T}}$。它们的内积为${\boldsymbol{u}}_1^{\rm{T}}{{\boldsymbol{u}}_2} = - 0.1070$，即它们是非正交的。

使用PSA以及MPSA方法求解时，当初始值距离${{\boldsymbol{u}}_1}$更近时，其会先收敛到${\boldsymbol{u}}_1^{{\rm{PSA}}}$，再根据正交约束原则，生成初始值${\boldsymbol{u}}_2^{{\rm{PSA}}}$。其最终收敛结果为${\boldsymbol{u}}_1^{{\rm{PSA}}} = {[0.8812, - 0.4727]^{\rm{T}}}$,${\boldsymbol{u}}_2^{{\rm{PSA}}} = {[0.4727,0.8812]^{\rm{T}}}$。

当初始值距离${{\boldsymbol{u}}_2}$更近时，其会先收敛到${\boldsymbol{u}}_2^{{\rm{PSA}}}$，再根据正交约束原则，生成初始值${\boldsymbol{u}}_1^{{\rm{PSA}}}$。其最终收敛结果为${\boldsymbol{u}}_1^{{\rm{PSA}}} = {[0.3756,0.9268]^{\rm{T}}}$, ${\boldsymbol{u}}_2^{{\rm{PSA}}} = [0.9268, - 0.3756]^{\rm{T}}$。上述结果如图1所示。

图 1  不同初始值，PSA, MPSA算法求解结果

下载: 全尺寸图片 幻灯片

使用NPSA求解上述2阶张量时，当初始值距离${{\boldsymbol{u}}_1}$更近时，其收敛到${\boldsymbol{u}}_1^{{\rm{NPSA}}}$，再根据其约束原则，收敛到${\boldsymbol{u}}_2^{{\rm{NPSA}}}$。其收敛结果为${\boldsymbol{u}}_1^{{\rm{NPSA}}} = [0.8812, - 0.4727]^{\rm{T}},{\boldsymbol{u}}_2^{{\rm{NPSA}}} = {[0.3351,0.9422]^{\rm{T}}}$。

当初始值距离${{\boldsymbol{u}}_2}$更近时，其收敛到${\boldsymbol{u}}_2^{{\rm{NPSA}}}$，再根据其约束原则，收敛到${\boldsymbol{u}}_1^{{\rm{NPSA}}}$。其收敛结果为${\boldsymbol{u}}_1^{{\rm{NPSA}}} = [0.3756, 0.9268]^{\rm{T}},{\boldsymbol{u}}_2^{{\rm{NPSA}}} = {[0.8799, - 0.4752]^{\rm{T}}}$。上述结果如图2所示。

图 2  不同初始值，NPSA算法求解结果

下载: 全尺寸图片 幻灯片

这个简单的例子直观地说明了现有的 PSA 算法并不能准确地得到协偏度张量的实际特征对，还证明了现有的 PSA算法随机生成初始值，导致每次运行都会收敛到不同的特征对。并且上述问题会随着张量维数的增加表现得更加明显，在3阶2维超对称张量求解中，可能会收敛到两个不同的特征对，在3阶$n$维超对称张量求解中，因为初始值的不同可能收敛到$n!$个不同的特征对。因此如何设计出计算求解精确、稳定的方法至关重要。

3.2   PPSA算法

针对上述问题，本文在PSA算法基础上进行改进，提出一种高精度并行PSA算法。PPSA算法不再施加约束，而是将各个特征对的初始值并行地全部算出，有效地避免了施加约束和串行计算带来的误差。

求解协偏度张量的实际特征对可以理解为求解张量的局部极大值解^[18]，只要保证初始值落在协偏度张量各个收敛区域，然后再使用不动点迭代法进行求解，就能保证精确地求解到所有极大值解。因此求解协偏度张量的特征对问题可转化为初始值选取的问题。常见的初始化方法有随机初始化，基于正交约束的随机初始化。随机初始化的方法虽然可以精确求解，但是它们受初始值影响很大，并且是不可重复的。因此设计出一种初始值可重复，可以稳定的精确的求解出协偏度张量的特征对的方法是很重要的。

PPSA算法充分考虑数据结构，选用协偏度张量的全部切片的特征向量作为迭代的初始值。首先，此方法是可重复的。其次，此方法是并行的。最后，对比随机初始化方法，此方法可以获得稳定精确的解的原因是我们相信它将提供在最优解附近的初始猜测。

假设${\boldsymbol{S}} \in {\mathbb{R}^{[3,n]}}$，如果$\lambda$是协偏度张量${\boldsymbol{S}}$的一个特征值，设${\boldsymbol{u}}$是${\boldsymbol{S}}$关于$\lambda$的单位特征向量，即${\boldsymbol{u}} = {[{a_1},{a_2}, \cdots ,{a_n}]^{\rm{T}}}$且${{\boldsymbol{u}}^{\rm{T}}}{\boldsymbol{u}} = 1$。求解${\boldsymbol{S}}{ \times _1}{\boldsymbol{u}}{ \times _3}{\boldsymbol{u}} = \lambda {\boldsymbol{u}}$，其等价求解

由式(7)可知，求解协偏度张量的局部极值解等价多个主成分分析联合求极大极小值。把协偏度张量的切片的特征向量作为迭代初始值可以理解为从单个主成分的极值点去逼近联合主成分的极值点。因此对比随机选取初始值，此初始化方法生成的初始值距离实际特征对更近。为了验证这一结论，本文在4.1节做了蒙特卡罗仿真实验。

为方便起见，总结PPSA算法如算法1所示。

算法1 PPSA算法
输入：输入数据${\boldsymbol{R}} \in {{\boldsymbol{R}}^{L \times N} }$
输出：输出变换矩阵${\boldsymbol{U}}$, ${\boldsymbol{Y}} = {{\boldsymbol{U}}^{\rm{T}}}\tilde {\boldsymbol{R}}$
(1) 白化数据，得到$\tilde {\boldsymbol{R}}$
(2) 根据式(1)计算m阶$L$维张量${\boldsymbol{S}}$
(3) 计算张量${\boldsymbol{S}}$切片的所有特征向量，记为V% main loop: %
(4) ${\text{for }}i = 1:{L^2}{\text{ do}}$
(5) 　　k = 0
(6) 　　${\boldsymbol{u}}_i^{(k)} = {\boldsymbol{V}}(:,i)$
(7) 　　while stop conditions are not meet do
(8) 　　　　${\boldsymbol{u}}_i^{(k + 1)} = {\boldsymbol{S}}{ \times _1}{\boldsymbol{u}}_i^{(k + 1)}{ \times _3}{\boldsymbol{u}}_i^{(k + 1)}$
(9) 　　　　${\boldsymbol{u}}_i^{(k + 1)} = {\boldsymbol{u}}_i^{(k + 1)}/{\left\| {{\boldsymbol{u}}_i^{(k + 1)} } \right\|_2}$
(10) 　　end while
(11) 　　${{\boldsymbol{U}}_{:i} } = {\boldsymbol{u}}_i^{(k + 1)}$
(12) end for

下载: 导出CSV 

| 显示表格

需要注意的是，PPSA算法中步骤7中的循环终止条件包含最小允许误差$\varepsilon$和最大循环次数${{K}}$。本文$\varepsilon$设置为0.0001，${{K}}$设置为1000。${\boldsymbol{U}}$是最终的非正交主偏度变换矩阵。

同样对上述3.1节的例子进行计算，用PPSA算法重复上述操作，求解得到的结果为

如图3所示，PPSA算法获得结果与${{\boldsymbol{u}}_1}$, ${{\boldsymbol{u}}_2}$完美重合，并且PPSA算法求解结果唯一。

图 3  PPSA算法求解结果

下载: 全尺寸图片 幻灯片

4.   实验与复杂度分析

本文将PPSA算法应用于盲图像分离、SAR图像去噪以及高光谱特征提取，并与几种经典方法开展了对比试验。所有的算法都是在1台CPU AMD Ryzen 7 5800H, 16 GB RAM, @3.20 GHz的笔记本电脑上完成的，所有程序在MATLAB R2021a上编程和实现。

4.1   蒙特卡罗仿真实验

为了检验本文所提出的初始化算法的效果，本文采用蒙特卡罗方法进行验证。实验采用有$n$个极大值特征对的3阶$n$维超对称张量进行求解，$n$分别取2～9。采用3种不同的初始化方法进行迭代求解，进行10万次蒙特卡罗仿真后，对比求解正确的结果，其结果如表1所示。

表 1  不同初始化方法的比较，选取n个初始值计算100 000次的结果

	随机选取$n$个初始值	随机选取$n$个正交初始值	单个切片的特征向量作为初始值
2	49027	89073	98270
3	12864	53080	98029
4	7450	35373	97222
5	2361	15565	85831
6	836	5171	50430
7	89	281	35851
8	27	129	19981
9	8	29	17004

下载: 导出CSV 

| 显示表格

从实验结果上看，基于协偏度张量的单个切片的特征向量作为初始值的方法好于基于正交约束的随机初始化方法以及随机初始值的初始化方法。但是随着协偏度阶数的增加，此初始化得到全部解的概率会有所下降。为了改善这种情况，可以适当增加初始值的数量，选用全部切片的特征向量作为迭代初始值重复上述实验，其结果如表2所示。

表 2  不同初始化方法的比较，选取n × n个初始值计算100 000次的结果

	随机选取$n \times n$个初始值	随机选取$n \times n$正交初始值	全部切片的特征向量作为初始值
2	90833	98184	99997
3	88401	92830	99999
4	88225	90029	99989
5	85646	88451	99994
6	79673	70096	99988
7	79382	69156	99974
8	77576	68546	99991
9	64338	66570	99938

下载: 导出CSV 

| 显示表格

从实验结果上看，随着初值点数的增加，各种初始化方法得到所有特征对的概率都有明显提升。尤其是基于全部切片的特征向量作为初始值的方法，其在所有的实验中准确率均达到了99.9%。因此，对比随机初始化方法，选用协偏度张量的全部切片的特征向量作为迭代的初始值距离实际特征对更近。

4.2   盲图像分离实验

本文首先将PPSA算法应用于盲图像分离(Blind Image Separation, BIS) 问题上，BIS的目的是估计混合矩阵，记为${\boldsymbol{B}}$(或其逆矩阵，即分解矩阵，记为${\boldsymbol{U}}$)。本文将PPSA算法与FastICA^[8], PSA^[10], MPSA^[11], NPSA^[15]和MSDP^[18]等算法进行比较，评价各算法的分离性能。在本实验中，本文选取$n$幅大小为256像素 × 256像素的灰度图像作为源图像，其中$n$分别取值2～6。由于已知真实的源图像，因此可以将其与各算法获得的结果进行比较。由于篇幅的限制，本文没有显示所有实验结果，只选择$n = 3$的情况进行展示，结果如图4所示。为了确保结论的可靠性，本文还进行了其他两种组合实验，在每种组合中本文随机选择3种不同的源图像进行实验。

图 4  源图像以及各种算法处理结果

下载: 全尺寸图片 幻灯片

从实验结果来看，上述6种算法均可以从混合图像中分离出源图像，为了定量评估上述6种算法的性能，本文使用5个评估指标对上述6种算法获得的分离结果进行评价。5个评价指标分别为符号间干扰(InterSymbol-Interference, ISI)^[15]、总均方误差(Total Mean Square Error, TMSE)^[15]、相关系数指数(Correlation Coefficient, 记为$\rho$)^[15]、峰值信噪比(Peak Signal to Noise Ratio, PSNR)^[15]和运行时间(Time, T)。值得注意的是除PPSA算法与MSDP算法外，其余4种算法获得的结果均具有随机性^[18]，因此本文取10次运行的平均值作为最终结果。

表3详细地列出了3种不同组合下6种算法的符号间干扰、总均方误差、相关系数、峰值信噪比和运行时间。对于ISI和TMSE指标，指数值越小，代表该方法的性能越好；$\rho$和PSNR越大，代表该方法的性能更好。通过对表3进行分析可以得到如下结论：

表 3  FastICA, PSA, MPSA, NPSA, MSDP和PPSA算法评估结果

	指标	FastICA	PSA	MPSA	NPSA	MSDP	PPSA
1	ISI	0.6757	0.0681	0.0551	0.0253	0.0203	0.0203
	TMSE	3.2856×10–11	2.5547×10–11	2.5548×10–11	3.8168×10–11	1.8645×10–11	1.8645×10–11
	$\rho$	0.9544	0.9954	0.9954	0.9992	0.9998	0.9998
		0.9889	0.9992	0.9988	0.99997	0.9992	0.9992
		0.9968	1.0000	1.0000	1.0000	1.0000	1.0000
	PSNR	66.1696	76.6105	77.7449	83.6350	85.9129	85.9129
		70.0325	72.6448	73.7954	74.4358	74.4870	74.4870
		79.2579	90.9901	90.4324	91.1667	91.3285	91.3285
	T(s)	0.0042	0.0011	0.0010	0.0043	2.8594	0.0032
2	ISI	0.4844	0.1442	0.1524	0.0787	0.0745	0.0745
	TMSE	3.2241×10–11	2.6798×10–11	2.6337×10–11	2.5154×10–11	1.9555×10–14	1.9555×10–14
	$\rho$	0.9749	0.9967	0.9967	0.9995	0.9997	0.9997
		0.9900	0.9905	0.9907	0.9920	0.9921	0.9921
		0.9970	0.9999	0.9998	1.0000	1.0000	1.0000
	PSNR	69.1886	74.9389	75.6708	83.4588	83.7674	83.7674
		73.7143	73.1186	74.3531	70.1218	70.6017	70.6017
		76.9175	89.0466	90.4683	94.0776	94.7823	94.7823
	T(s)	0.0052	0.0011	0.0010	0.0042	3.1046	0.0021
3	ISI	0.2965	0.0405	0.0452	0.0041	0.0006	0.0006
	TMSE	2.6603×10–10	2.6628×10–10	1.9986×10–10	3.4050×10–10	1.7480×10–10	1.7480×10–10
	$\rho$	0.9756	0.9966	0.9962	0.9997	1.0000	1.0000
		0.9950	0.9994	0.9995	1.0000	1.0000	1.0000
		0.9969	0.9998	0.9999	1.0000	1.0000	1.0000
	PSNR	69.9228	83.2559	78.8206	91.4669	96.6388	96.6388
		61.8123	61.1492	58.8872	61.2297	64.0212	64.0212
		82.3327	87.5416	90.4749	94.0424	95.1910	95.1910
	T(s)	0.0048	0.0012	0.0009	0.0033	2.9786	0.0022

下载: 导出CSV 

| 显示表格

(1)MSDP算法可以准确得到张量的特征对^[18]，通过观察MSDP算法与PPSA算法的评价指标，可以发现PPSA算法与MSDP算法获得的结果一致，这说明PPSA算法也可以准确得到张量的实际特征对，验证了本文所提方法的有效性。同时，PPSA算法的运行时间远小于MSDP算法的运行时间，这说明PPSA算法在精确求解的同时还可以做到快速求解。

(2)PPSA算法获得的ISI与TMSE指标始终小于其余4种方法的，ISI与TMSE指标越小，图像分离效果越好，这说明PPSA算法在BIS应用中可以获得更好的整体效果。

(3)PPSA算法获得的相关系数与峰值信噪比也在这6种算法中处于较高水平，相关系数值和峰值信噪比越大，分离图像质量越好，这说明PPSA算法性能稳定。

最后，结合上述5个指标的比较结果，说明PPSA算法在BIS应用中具有更准确和更鲁棒的性能。

4.3   SAR图像去噪实验

本文将PPSA算法用于SAR图像去噪问题上，实验使用的是高分三号拍摄的一幅美国旧金山图像，图像大小为448像素×364像素，极化方式为HV极化。对源图像施加噪声来模拟低信噪比图像，通常认为SAR图像的斑点噪声服从伽马分布或者高斯分布，因此分别在源图上施加均值为0，方差为 100的高斯噪声(记为${E_{{\text{gs}}}}$=0，$\delta _{{\text{gs}}}^2$=100)和均值为2，方差为80的伽马噪声(记为${E_{{\text{gm}}}}$=2，$\delta _{{\text{gm}}}^{\text{2}}$=80)来仿真低信噪比的混合图像。将原始图像按逐行的方式转化为相应的3个1维信号，产生随机混合矩阵 rand(3,3)。1维信号与随机混合矩阵混合后得到混合图像，再使用FastICA, PSA, NPSA, PPSA、均值滤波和中值滤波等算法对混合图像进行去噪处理，结果如图5所示。

图 5  原始数据以及各种算法处理结果

下载: 全尺寸图片 幻灯片

因为FastICA, PSA和NPSA的处理结果受随机初值影响，所以本文选择它们10次实验的最佳的处理结果进行对比以及展示。同样，本文也选择出均值滤波以及中值滤波最佳的结果进行对比展示，均值滤波与中值滤波均选择5×5大小的滤波器。

通过观察图5，发现 PSA算法、NPSA算法、PPSA算法的去噪效果优于FastICA算法、中值滤波和均值滤波。为了客观评价分离滤波结果，本文分别计算图像去噪后的峰值信噪比(PSNR)^[15]和平均绝对误差(MAE)。对像素为M×N的图像，其MAE 定义为

式中，$h(x,y)$和$g(x,y)$分别表示原始图像与待评价图像的第$(m,n)$个像素的灰度值。为了保证结论的可靠性，本文还选择不同参数的伽马噪声与高斯噪声进行了2次实验，计算结果如表4所示。

表 4  不同算法的平均峰值信噪比和平均绝对误差(计算10次运行的平均结果)

组合	评估指标	FastICA	PSA	NPSA	PPSA	均值滤波	中值滤波
1	PSNR	64.3293	83.6292	83.9879	84.1115	62.0892	60.9464
	MAE	0.0719	0.0127	0.0120	0.0119	0.1749	0.1957
2	PSNR	61.3562	81.4705	81.8815	81.9711	62.7697	62.3189
	MAE	0.1631	0.0171	0.0163	0.0161	0.1574	0.1659
3	PSNR	68.6684	78.8200	78.8824	78.9555	64.0314	64.6669
	MAE	0.0666	0.0245	0.0243	0.0241	0.1306	0.1233
(注释：组合1中${E_{{\text{gs}}}}$=0, $\delta _{{\text{gs}}}^{\text{2}}$=100, ${E_{{\text{gm}}}}$=2, $\delta _{{\text{gm}}}^{\text{2}}$=80；组合2中${E_{{\text{gs}}}}$=0, $\delta _{{\text{gs}}}^{\text{2}}$=80, ${E_{{\text{gm}}}}$=2, $\delta _{{\text{gm}}}^{\text{2}}$=60；组合3中${E_{{\text{gs}}}}$=0, $\delta _{{\text{gs}}}^{\text{2}}$=60, ${E_{{\text{gm}}}}$=2, $\delta _{{\text{gm}}}^{\text{2}}$=40)

下载: 导出CSV 

| 显示表格

MAE越小，算法性能越好。由表4可看出，对于相同强度的高斯噪声和伽马噪声, 均值滤波法和中值滤波法的去噪效果相对较差，而 FastICA 方法相对较好，PSA方法在所有组合中均获得了不错的效果，尤其是PPSA方法，能显著地去除噪声，提高复原图像的质量。结合上述两个指标的比较结果，说明PPSA算法在SAR图像去噪处理中具有更强的去噪能力。

4.4   高光谱特征提取实验

本节使用高光谱数据集Cuprite来评价PPSA算法的性能。该数据集是由机载可见红外成像光谱仪于1997年6月在内华达州铜矿矿区获得的。将数据嵌入到可视化图像环境软件中，图像包含50个波段，每个波段有350像素×400像素，伪彩色图像如图6所示。

图 6  高光谱Cuprite数据彩色图像

下载: 全尺寸图片 幻灯片

本文直接选择数据集Cuprite全部的50个波段进行特征提取。当选用全部50个波段进行特征提取时，MSDP和FastICA算法的时间消耗是巨大的^[18]，因此本文使用PSA, MPSA, NPSA和PPSA算法进行特征提取。与4.2节的情况不同，在实际数据处理中，通常不知道独立分量的数量，所以PSA, MPSA和NPSA算法的独立分量的数量总是手动设置为频带的数量。因此用上述3种算法进行求解，它们最终全都获得50个独立分量。相比之下，PPSA算法是参数自适应的，它不需要事先设置独立分量的数量，它可以自动确定独立分量的数目，用PPSA算法进行求解，最终获得了6个独立分量。由于不知道真实数据的参考独立分量，因此4.2节仿真实验部分使用的评价指标不可用。因为篇幅的限制，本文不一一显示PSA, MPSA和NPSA这3种算法获得的结果。本文仅显示PPSA算法获得的结果，如图8所示。

图 8  PPSA算法对高光谱图像Cuprite的偏度比较

下载: 全尺寸图片 幻灯片

其余3种算法虽然均获得了50个独立分量，但是所提取到50个特征对中大多是没有信息的冗余结果，并且提取到的特征对除了第1个外其余都为近似解。把PSA, MPSA和NPSA这3种算法得到的50个独立分量分别作为初始值进行固定点迭代求解，最终全都收敛为6个独立分量，这6个独立分量与PPSA算法得到的结果完全相同，说明PPSA算法在对高维数据进行特征提取时具有更好的自适应性和准确性。

此外，本文还对上述4种算法的运行速度进行比较。PPSA算法可以并行地求解各特征对，因此PPSA算法可以使用Parfor并行加速计算。本文将使用Parfor并行加速计算的PPSA算法命名为FastPPSA。

图7绘制了不同数量波段下4种算法的运行时间曲线，根据Cuprite数据各波段的偏度值，从小到大排列依次选取$n$个波段进行特征提取，$n$的范围为2～50。

图 7  4种算法运行时间对比图

下载: 全尺寸图片 幻灯片

通过观察图7可以看出在处理相同数量的数据样本时(样本数大于20时)，MPSA算法运行时间小于FastPPSA算法运行时间小于PSA算法运行时间小于NPSA算法运行时间小于PPSA算法运行时间。同时，还可以看出在处理高阶数据时FastPPSA算法运行时间远小于PPSA算法运行时间。然而在处理低阶数据时FastPPSA运行时间要略大于PPSA算法的运行时间。这是因为当计算机处理的数据量较小时，Parfor 并行加速在初始化过程中需要消耗大量的时间，这部分时间在总的运行时间里占了很大的一部分。因此，当所要处理图像规模较小时，利用Parfor并行运算的加速效果并不明显。然而，随着张量维数的逐渐增大，耗费在并行计算额外开销上的时间占算法执行总时间的比例越来越小，所以并行加速效果明显。FastPPSA采用Parfor并行计算导致其在处理高维数据时运算速度很快，处理50个样本数据，其运行速度为3.1028 s, NPSA运行速度为17.2797 s, PPSA算法运行时间为28.1447 s。

结合上述实验结果与分析可知，本文所提出的PPSA算法不仅能够准确地对高光谱数据进行特征提取，而且还可以做到快速提取。这说明PPSA算法在高光谱数据特征提取上具有更好的准确性、自适应性和鲁棒性。

4.5   复杂度分析

本节从理论上比较PSA算法和PPSA算法的计算复杂度。由于这两种方法基本流程都一样，只是在求解张量特征对时主循环的次数不同，因此只需考虑主循环求解的计算复杂性。对于$L \times M$大小的数据，PSA算法只需要循环$L$次，假设每次循环需要迭代$K$次，其中$K$是平均迭代次数，则 PSA算法需要进行$KL$次操作。同样对于$L \times M$大小的数据，PPSA算法则需要循环${L^2}$次，同样也假设每次循环需要迭代$K$次，因此PPSA算法最终需要$K{L^2}$次操作。虽然本文PPSA算法的计算复杂度略高于PSA算法，但PPSA算法相对于PSA算法的主要优势在于PPSA算法是并行的，所以它可以采用并行计算来提升算法运行速度。由前面的高光谱运行速度对比实验可知，在处理高维数据时，采用并行加速计算后的PPSA算法运行速度更快。

5.   结束语

本文针对现有PSA算法不能准确得到实际解这一问题，在原有的PSA算法的基础上进行改进，提出了一种高精度并行PSA算法。与现有PSA算法串行求解特征对不同，本文所提算法不再施加任何约束，而是选择张量全部切片的特征向量作为初始值进行自由迭代求解，有效地避免了施加约束带来的误差。通过将PPSA算法与现有PSA算法进行实验对比，验证了PPSA算法在盲图像分离，SAR图像去噪以及高光谱数据特征提取上具有更好的准确性和鲁棒性。PSA, NPSA以及PPSA都关注数据集的3阶偏度。对于其他数据集，偏度可能并不总是最佳选择。在这种情况下，可以选用其他指标，如4阶峰度或其他高阶统计量替代。