基于张量分析的欠定混合矩阵估计算法
时间：2023-06-08 11:25:15 来源：雅意学习网 本文已影响 人 
马宝泽，李国军，向翠玲，徐阳
（1.重庆邮电大学光电工程学院，重庆 400065；
2.重庆邮电大学超视距可信信息传输研究所，重庆 400065；
3.重庆邮电大学光电信息感测与传输技术重庆市重点实验室博士后科研工作站，重庆 400065）
盲源分离（BSS,blind source separation）是指在仅以观测混合信号为先验信息的情况下，实现信道混合矩阵辨识和期望信源恢复的信号处理技术[1-4]。其中，信道混合矩阵估计是盲源分离问题研究的关键，直接影响着估计信源的分离精度[5-6]。传统的独立分量分析（ICA,independent component analysis）方法、独立矢量分析（IVA,independent vector analysis）方法能够处理超定和正定情况下的混合信道矩阵估计问题[2]，但不能用于欠定混合信道矩阵盲辨识。由于欠定线性混合系统中方程的个数小于未知数的个数，即观测信号数少于源信号数，这将导致估计混合矩阵不存在伪逆，给解决欠定混合矩阵估计问题带来挑战[7-9]。因此，欠定混合信道矩阵的盲辨识问题具有重要的研究价值，对实现欠定情况的源信号恢复至关重要。
目前，解决欠定混合矩阵估计的方法主要集中在观测信号的特定属性，如稀疏变换[6,10]、四阶累计量[11]、二阶统计量[12]、特征函数[13]、张量分析[14-16]等。其中，利用信号在变换域的稀疏性解决欠定盲辨识问题是常用的方案之一，代表算法包括TIFROM（time-frequency ratio of mixtures）[10]以及基于稀疏分量分析的欠定盲源分离（UBSS-SCA,underdetermined blind source separation-sparse component analysis）算法[6]。稀疏变换类方法利用信号在稀疏域的散点数据线性聚类特性依次估计混合信道矩阵的方向向量，然而观测信号个数以及数据的稀疏程度均会限制该类方法的应用范围。当源信号在变换域非稀疏但是满足统计独立性且具有时间结构特征时，可以根据观测信号的高阶统计特性完成欠定混合矩阵估计任务。文献[11]利用信号在不同时延下计算的四阶累积量提出了FOBIUM（fourth-order cumulant-based blind identification of underdetermined mixtures）算法，但是估计四阶累积量是比较复杂的工作且需要较大的数据量才能获取有效的累积量统计信息。为了提升欠定混合矩阵盲辨识的性能并降低算法计算复杂度，文献[12]通过观测信号时延数据得到的空间协方差矩阵提出了SOBIUM（second-order covariance based blind identification of underdetermined mixtures）方法。FOBIUM 和SOBIUM 适用于时域信息丰富的信号，可利用信号的时延统计信息实现欠定盲辨识。随后，文献[13]针对复值数据提出了基于观测信号特征函数构造张量模型的欠定盲辨识方法，为构造高阶张量数据提供了研究思路。此外，文献[14]利用信号分割策略提出了基于张量分析的大规模系统盲辨识方法，文献[15]将张量分析理论应用到欠定模态识别领域，文献[16]则提出在线稀疏张量重构方法可以作为多维信号处理的潜在技术。综上所述，欠定混合矩阵估计问题大致可分为基于稀疏变换[6,10]和高阶统计量[11,12,14-16]2 个解决思路。
本文从理论上阐述了信道盲辨识和张量分析的关系，将欠定混合矩阵估计问题转化为张量分解问题，并分析了由自协方差矩阵构成的三阶张量的对称性。为了提取观测信号的有效统计信息构造张量模型，同时降低计算复杂度并提高收敛速度，本文提出了一种新的欠定混合矩阵估计算法。首先，根据信号分割[1,14]策略计算观测信号每个子段的自协方差矩阵构造对称的三阶张量；
其次，通过截断多线性奇异值分解（MLSVD,multi-linear singular value decomposition）方法将原张量降维压缩成一个低维的核张量；
再次，引入增强的线性搜索（ELS,enhanced line search）技术[17]来加速交替最小二乘（ALS,alternating least square）算法[18]的收敛；
最后，利用对称核张量分解得到的因子矩阵作为估计混合矩阵测度，从而实现信道盲辨识。实验表明，本文算法在处理欠定混合矩阵估计问题时的估计误差均小于稀疏变换方法和传统的高阶统计量方法。
信道矩阵的盲辨识研究是源信号恢复的关键，特别是在欠定混合的条件下。假设n个源信号表示为，m个观测信号表示为，其中，k表示离散时间，N表示信号长度，(·)T表示转置运算。由于瞬时混合过程中传输信道仅考虑了源信号在传输过程中幅度衰减的影响，因此混合系统可以用随机生成的矩阵A∈ℝm×n表示。则第i个传感器接收到的观测信号可以表示为
其中，aij表示第j个源信号到达第i个传感器的衰减系数，即混合矩阵A的第(i,j)个元素，s j(k)表示第j个源信号，ni(k)表示第i个观测信号对应的加性噪声。
线性瞬时混合过程的矩阵形式可以表示为
瞬时混合信道盲辨识中的矩阵估计问题可以转化为对称张量的分解问题，为了阐释混合矩阵估计和张量分解的关系，需要先介绍有关张量分解的基本概念。应用广泛的张量分解技术主要包括CP（CANDECOMP/PARAFAC）分解和Tucker 分解，两者都是张量奇异值分解（SVD,singular value decomposition）的扩展版本。CP 分解的计算复杂度为仅与张量阶数有关的线性形式，而Tucker 分解则存在“维数灾难”的问题，即计算复杂度随着张量阶数以指数的形式增加[19]。因此，本文选择CP 分解估计欠定混合矩阵，建立欠定信道盲辨识与张量分解的关系。
CP 分解能够将一个高阶张量分解为一系列秩一分量之和的形式，不失一般性，给定一个三阶张量 X ∈ℝI×J×K，经CP 分解后可以表示为
利用分割策略[1,14]构造三阶对称张量实现欠定混合矩阵估计的流程如图1 所示，首先，将观测信号分割为不重叠的若干个子块；
然后，计算每个子块的自协方差矩阵并堆叠成张量的形式；
最后，根据对称张量的分解唯一性理论实现混合矩阵的估计。通过图1 可以看出，张量分解后由vi和bi构成的因子矩阵中包含信道矩阵信息，因此可以将因子矩阵V或者B作为估计混合矩阵，图1 中以V为例，即。此外，ci构成的因子矩阵包含信号的分块信息。分割策略计算自协方差矩阵构造张量的方式可以实现将混合矩阵估计问题转化为张量分解。同时，张量分解的唯一性能够保障混合矩阵估计的可行性和有效性。
图1 分割策略应用于欠定混合矩阵估计的流程
由上述分析可知，Cx中存在2 个相等的模态，故称Cx为对称张量。如果存在m=P，则认为三阶张量Cx具有超对称性。不失一般性，对称张量可以通过CP 分解获得2 个维度相同的因子矩阵，即
其中，2 个矩阵A可以视为不同的因子矩阵，分别记为AL和AR。AL和AR独立更新且不存在强制相等的显式约束，自协方差矩阵的固有对称性将使这2 个因子矩阵最终收敛，并且按照某个对角矩阵缩放。研究表明，只要CP 分解具有唯一性，就能够得到对称张量分解的解决方案。
通过MLSVD 对构造的三阶对称张量Cx降维，将原始张量压缩为维度较小的核张量。由三阶张量CP 分解唯一性的充要条件可知n≤P，不妨设定核张量的最大维度不超过n+1，即观测信号分割构造的自协方差张量第三个维度需要由P降到n+1，得到m×m× (n+1)的核张量。其中，P的选取是个开放性问题，需要根据具体的应用场景、信号属性以及张量分解的唯一性约束等确定。首先，三阶张量按照Mode-k形式展开的矩阵分别记为。然后，分别计算矩阵C1、C2和C3的SVD，即
传统的CP 分解问题通常采用ALS 算法[18]解决，每次迭代时每个步骤都包含将未知参数子集固定到当前估计中，然后对剩余的未知参数集进行修正，最后对互补参数集进行优化。因此，张量分解的代价函数可以表示为
令矩阵X∈ℝKJ×I表示张量X 的矩阵展开形式，记为XKJ×I，在给定因子矩阵V、B、C在第t-1和t-2次迭代的估计值情况下，第t次迭代中的最优松弛因子μ可以通过最小化代价函数得到，则式(11)可以改写为
其中，⊙表示Khatri-Rao 乘积。则式(13)可简化为
为了方便符号标记，将式(13)中的上标t和t-2省略。在式(14)中，T3、T2、T1、T0分别定义为
算法1本文算法流程
步骤1计算观测信号分割后每个子段的自协方差矩阵，利用式(8)构造三阶对称张量 Cx∈ℝm×m×P。
步骤2通过截断MLSVD 方法对张量Cx进行降维，利用式(10)得到左奇异矩阵计算三阶对称的核张量 G∈Km×m×（n+1）。
步骤3将核张量G 作为ELS+ALS 优化张量分解方法的输入，对因子矩阵进行初始化并更新。
①采取强制随机初始化的方式实现因子矩阵的初始化，即初始化
② 根据式(13)～式(15)将代价函数转化为多项式并求偏导计算最优的松弛因子μ，使式(14)取得极小值。
步骤4将步骤3 得到的初始化因子矩阵和原始对称张量Cx作为ALS 算法的输入，重复步骤3直到达到迭代终止条件，最终实现欠定瞬时混合信道矩阵的盲辨识，即不妨将因子矩阵V作为估计混合矩阵。
本文算法的计算复杂度分析从3 个方面展开，其一，所提算法和SOBIUM 算法构造的张量切片矩阵都属于二维代数结构，具有比FOBIUM 算法中的四阶累积量更低的复杂度。其二，对于三阶对称张量 Cx∈ℝm×m×P直接执行CP 标准分解时，令张量的秩为信源数n，则ALS 单次迭代的计算复杂度为[18]；
而压缩后的核张量G 经过ELS+ALS 分解的总复杂度可以表示为其三，与传统ALS 迭代算法相比，ELS+ALS 组合方法的复杂度相对较低，有效提高张量分解的收敛速度，降低算法的运算时间。
综上所述，本文算法在构造对称张量、张量压缩和增强线性搜索等方面能够降低复杂度，证明了所提算法的有效性和可行性。此外，本文的张量分析概念包含张量分解与估计混合矩阵的关系、构造对称三阶张量、张量压缩以及核张量分解4 个环节，构建了以因子矩阵作为欠定混合矩阵估计测度的信道矩阵盲辨识方法。
为了验证本文算法在实数域欠定混合矩阵估计中的可行性和有效性，分别采用LI-TIFROM[10]、UBSS-SCA[6]、SOBIUM[12]作为对比算法。源信号来自文献[2]的语音信号和音乐信号，其中，采样频率为16 kHz，语音时长为10 s。由于时域语音和音乐信号都是由实数组成的，因此该实验中的瞬时混合矩阵为随机生成的m×n维实数矩阵。从语音信号集中任意选取n个作为源信号，通过混合信道矩阵后得到m个观测信号，直接对时域观测信号进行张量化处理并借助ELS+ALS 实现混合信道矩阵的盲辨识。此外，为了验证本文算法对通信信号的有效性，本文采用多维正交调幅（QAM,quadrature amplitude modulation）信号作为信源估计混合矩阵。实验采取的性能指标为信道矩阵的平均相对误差，可以表示为error=[12,15]。其中，表示混合矩阵A最佳的排序和尺度估计。此外，平均相对误差越小表明估计混合矩阵越接近真实的混合矩阵，说明本文算法在信道矩阵盲辨识方面性能越优异。
5.1 实测语音信号的欠定混合矩阵估计
针对瞬时混合的实数域非平稳语音信号，分别在混合矩阵为3×4、4×5和5×6这3 种应用场景下验证各个算法在处理不同混合矩阵时的估计能力。由于分割子段长度大于或等于40 ms 时才能保证语音信号的非平稳性[3]，则P应小于或等于250。不妨令P=250，在每个混合信道场景下随机生成50 个实数域混合矩阵，每个混合矩阵中包含m×n个未知元素，不同信道矩阵下的估计误差（语音）如图2 所示。从图2 中可以看出，相较于对比算法，本文算法在3 种场景下均取得了最小的信道矩阵估计误差，说明本文算法在不同混合信道下取得的欠定混合矩阵估计性能优于其余对比算法。此外，随着信号个数的增加，LI-TIFROM 和本文算法的估计误差随之下降，而UBSS-SCA 和SOBIUM 算法的估计误差基本不变。
图2 不同信道矩阵下的估计误差（语音）
不失一般性，在混合矩阵为5×6的情况下分别验证步长为1、步长为1.25、步长为的3 种ALS算法以及ELS 优化的ALS 算法在最小化代价函数方面的迭代收敛快慢。图3 表示4 种算法计算代价函数随迭代次数的变化情况。从图3 可以看出，变步长的ALS 比传统的ALS 具有更快的收敛速度，且ELS+ALS 率先达到收敛状态，说明ELS 能够加快ALS 的收敛速度。
图3 代价函数随迭代次数的变化情况（语音）
为了验证本文算法在不同信噪比（SNR,signal to noise ratio）下的抗噪性能，在混合信道为3×4的情况下利用随机生成的高斯白噪声使SNR 取值范围为-5～50 dB 且间隔为5 dB。在不同SNR 条件下，4 种算法得到的欠定混合矩阵估计误差（语音）如图4 所示。从图4 可以看出，随着SNR 的增加，所有算法的估计误差都随之降低，当SNR>20 dB时，算法估计误差基本趋于稳定。本文算法在选取的SNR 区间内得到的矩阵估计误差最小，其次是SOBIUM 算法和LI-TIFROM 算法。
图4 不同SNR 下的估计误差（语音）
5.2 实测音乐信号的欠定混合矩阵估计
鉴于音乐信号通常包含重复的音色信息，故从由语音和音乐信号组成的数据库中选取对应的信源验证算法的欠定混合矩阵估计能力。为保证语音子段序列的非平稳性，选取P=250，分别在混合矩阵为3×4、4×5和5×6这3 种混合信道情况下随机生成50 个对应的混合矩阵，4 种算法在不同混合矩阵条件下的估计误差（音乐）如图5 所示。从图5 可以看出，本文算法取得了比其余算法更低的欠定矩阵估计误差，且可以取得与图2 基本一致的结论，即随着信号个数的增加，LI-TIFROM 和本文算法的估计误差逐渐降低。
图5 不同混合矩阵下的估计误差（音乐）
图6 代价函数随迭代次数的变化情况（音乐）
在不同SNR 情况下，4 种算法得到的欠定矩阵估计误差（音乐）如图7 所示。与图4 对比可知，本文算法取得的估计误差基本不变，说明音乐信号的引入没有改变本文算法的估计性能。此外，LI-TIFROM 和SOBIUM 的辨识性能受音乐信号的影响相对较大。
图7 不同SNR 下的估计误差（音乐）
5.3 QAM 信号的欠定混合矩阵估计
为了验证本文算法对以QAM 信号为代表的通信信号混合矩阵估计的可行性和有效性，分别在混合矩阵为3×4、4×5和5×6这3 种场景下随机生成相应的实数矩阵，模拟信号在传输过程中的幅度衰减情况。其中，QAM 信号的长度为1 000，采样频率为8 000 Hz，信道带宽为20 MHz，不妨令分割子段数P=200。由于源信号和观测信号均为复数序列，故得到的估计混合矩阵为复数矩阵。因此，提取估计混合矩阵的实部计算信道矩阵的平均相对误差，经过50 次仿真后的平均结果如表1所示。
表1 QAM 信号的混合矩阵平均相对误差
由表1 可知，随着混合矩阵形状增大，各算法的估计误差也随之增大。LI-TIFROM 和UBSS-SCA依赖信号的稀疏变换属性，在处理QAM 信号的复值观测数据时估计误差较大。本文算法估计性能优于对比算法，说明分割策略下的自协方差矩阵相较于SOBIUM 的时延协方差矩阵具有更明显的信号特征表达能力。
针对欠定信道矩阵盲辨识问题，提出了一种基于张量分析的欠定矩阵盲估计算法。首先，利用分割策略将观测信号分割为不重叠的子段，计算每个子段的自协方差矩阵构造对称的三阶张量；
然后，通过MLSVD 将原张量压缩为低维的核张量并进行张量分解；
最后，引入ELS 技术加速ALS 算法的收敛，得到的因子矩阵作为估计的欠定混合矩阵，但分割子段数的选取则受到信号属性和应用场景的约束。实验表明，本文算法在处理欠定混合矩阵估计问题时，相较于稀疏变换方法和传统的高阶统计量方法在估计误差方面具有优势。此外，从理论上研究复值信号的欠定混合矩阵估计问题将是未来工作的重点。
