0 概述

跳频扩频(Frequency-Hopping Spread Spectrum, FHSS)是指使用伪随机码进行频移键控, 使得信号发射的载波频率不断跳变同时频谱得到扩展的方法^[1]。FHSS技术具有保密性和组网能力强的优点, 将其用于军用通信, 能够有效避开干扰, 保证通信效能。

欠定盲源分离是指在信源数多于观测数的情况下, 从接收到的混合信号中分离、恢复出源信号的过程^[2]。跳频信号的盲源分离是进行信号侦察处理及电子对抗的一个关键步骤, 是开展情报获取、敌我识别、干扰引导和威胁等级分析等军事行动的前提。因此, 研究跳频信号盲源分离方法具有重要的理论意义和应用价值。

通过判断跳频网台之间跳变的时间是否同步, 可以将跳频网台组网方式分为同步组网和异步组网两种。异步组网是指在各跳频网台之间没有时间同步的情况下, 对跳频信号用时频分析工具得到每跳的到达时刻, 然后根据不同电台到达时刻和持续时长进行分离^[3]。该方法简单有效, 但其适用范围仅局限于异步跳频网台。同步组网则是基于通信系统内各跳频电台之间具有相同的时间同步, 即电台遵循相同跳变时间规律的情况。现有欠定条件下同步组网跳频信号的盲源分离方法主要采用“两步法”, 即第1步估计混合矩阵, 第2步恢复源信号^[4]。

混合矩阵的估计结果将直接影响后续源信号的恢复效果。文献[5-7]方法采用基于稀疏聚类的混合矩阵估计方法, 计算相对简单, 估计精度较高, 但是需要满足一定的稀疏条件, 应用范围有所限制。文献[8-10]基于张量分解的估计算法, 利用信号的高阶统计量构造张量, 将混合矩阵的估计问题转化为张量分解问题, 适用于信源相互独立并且非高斯分布的情况, 精度较高, 但是计算量较大。其中, 文献[8]方法将压缩感知理论与张量分析相结合, 在实现盲源分离的同时完成了波达方向(Direction of Arrival, DOA)估计, 但只适用于具有特殊表示形式的源信号。文献[9-10]方法把张量分析应用到盲源分离问题求解中, 将混合矩阵的估计问题转化为张量分解问题, 提高了估计的准确性和鲁棒性, 但其中用于CP(Canonical-Polyadic)分解的最小二乘(Alternating Least Square, ALS)算法存在对初值敏感、易陷入局部极值和运算时间较长的不足。

文献[11-13]针对源信号恢复问题进行了研究。文献[11]采用稀疏重构的方法恢复源信号。该方法可以实现稀疏信号的盲分离, 但无法有效分离非稀疏的源信号。文献[12]采用二元掩蔽法进行源信号的恢复。该方法计算简单, 适用于线性混合方式, 是目前较为常用的欠定盲分离算法, 源信号越稀疏其恢复效果越好, 但在实际应用中解决跳频信号分离问题的性能并不出色。针对非充分稀疏的混合信号, 文献[13]采用子空间投影法进行源信号的重构, 只要任意时频点同时存在的源信号个数小于阵元数, 即可利用子空间投影原理确定任意时频点上源信号对应的混合矩阵。

本文针对欠定条件下跳频信号的盲源分离问题, 按照上文所述的“两步法”, 提出一种基于平行因子分析模型与子空间投影法的盲源分离方法。将混合矩阵估计问题转化为张量CP分解问题, 利用直接三线性分解改进交替ALS算法对其求解, 同时在迭代过程中采用标准线搜索加速收敛得到混合矩阵。在此基础上, 通过子空间投影法完成盲源分离, 并剔除离散噪点进一步优化分离效果。

1 跳频信号模型

跳频信号s_n(t)的数学表达式^[14]描述为:

$ {\mathit{\boldsymbol{s}}_n}(t) = {v_n}(t)\sum\limits_{k = 0}^{K - 1} {{\rm{exp}}} [{\rm{j}}({f_{nk}}{t^\prime } + {\varphi _{nk}})] {\rm{rect}} \left( {\frac{{{t^\prime }}}{{{T_n}}}} \right) $

(1)

其中, T_n为跳频信号s_n(t)的跳周期, K为观测时间Δt内总跳数, f_nk为第k跳(k=1, 2, …, K)的载频, φ_nk为初始相位, 起始的非完整跳在观测时间内的持续时长为Δt_0n, v_n(t)是跳频信号s_n(t)的基带复包络, t′=t-(k-1)T_n-Δt_0n表示瞬时刻, rect表示单位矩形窗函数。

盲源分离是指从接收到的混合信号中分离、恢复出混合前源信号的过程。盲源分离系统模型如图 1所示。

设M个观测信号是来自于N个源信号的线性混合, 则盲源分离信号模型描述为:

$ \mathit{\boldsymbol{X}}(t) = \mathit{\boldsymbol{AS}}(t) + \mathit{\boldsymbol{N}}(t) $

(2)

其中, X(t)=[x₁(t), x₂(t), …, x_M(t)]^T表示接收的M个观测信号, S(t)=[s₁(t), s₂(t), …, s_N(t)]^T表示输入的N个未知源信号, N(t)表示与源信号独立的加性噪声, A为M×N维的混合矩阵。

在实际应用中, 根据信源数目和观测信号数目的关系, 可以将盲源分离模型分为超定、正定和欠定3类, 超定即观测数目大于源信号数目的情况, 正定即观测数目等于源信号数目的情况, 欠定即观测数目小于源信号数目的情况^[15]。

2 混合矩阵估计

本文基于平行因子模型PARAFAC进行混合矩阵估计。通过计算跳频信号时延相关矩阵构造三阶张量, 将混合矩阵估计问题转化为张量CP分解问题。同时改进ALS算法, 利用直接三线性分解方法粗估加载矩阵, 将其作为ALS算法的初始迭代矩阵, 并在迭代过程中采用标准线搜索加速收敛得到混合矩阵。

2.1 张量分析相关概念

定义1(向量的外积)^[16]向量a∈ℝ^I₁×1与b∈ℝ^I₂×1的外积是一个I₁×I₂维的矩阵C, 定义为:

$ \mathit{\boldsymbol{C}} = \mathit{\boldsymbol{a}} \circ \mathit{\boldsymbol{b}} = \mathit{\boldsymbol{a}}{\mathit{\boldsymbol{b}}^{\rm{T}}} $

(3)

定义2(秩-1张量)^[16]如果三阶张量χ等于3个向量的外积, 则其秩为1。

定义3(CP分解)^[16]任何张量χ都允许分解为秩-1张量的总和。在三阶张量的情况下, CP分解定义为:

$ \mathit{\boldsymbol{\chi }} = \sum\limits_{f = 1}^F {{\mathit{\boldsymbol{a}}_f}} \circ {\mathit{\boldsymbol{b}}_f} \circ {\mathit{\boldsymbol{c}}_f} $

(4)

其中, a∈ℝ^I₁×1, b∈ℝ^I₂×1, c∈ℝ^I₃×1, F为张量χ的秩, f=1, 2, …, F。

定义4(Kronecker积)^[16] I₁×I₂维矩阵A=[a₁, a₂, …, a_I2]与I₃×I₄维矩阵B的Kronecker积是一个I₁I₃×I₂I₄维的矩阵, 并有如下表达式:

$ \mathit{\boldsymbol{A}} \otimes \mathit{\boldsymbol{B}} = \left( {\begin{array}{*{20}{c}} {{a_{11}}\mathit{\boldsymbol{B}}}&{{a_{12}}\mathit{\boldsymbol{B}}}& \cdots &{{a_{1n}}\mathit{\boldsymbol{B}}}\\ {{a_{21}}\mathit{\boldsymbol{B}}}&{{a_{22}}\mathit{\boldsymbol{B}}}&{}&{{a_{2n}}\mathit{\boldsymbol{B}}}\\ \vdots & \vdots &{}& \vdots \\ {{a_{m1}}\mathit{\boldsymbol{B}}}&{{a_{m2}}\mathit{\boldsymbol{B}}}& \cdots &{{a_{mn}}\mathit{\boldsymbol{B}}} \end{array}} \right) $

(5)

定义5(Khatri-Ra积)^[16] I₁×I₀维矩阵A=[a₁, a₂, …, a_I0]和I₂×I₀维矩阵B=[b₁, b₂, …, b_I0]的Khatri-Rao积是一个I₁I₂×I₀维的矩阵, 并有如下表达式:

$ \mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{B}} = [{\mathit{\boldsymbol{a}}_1} \otimes {\mathit{\boldsymbol{b}}_1}, {\mathit{\boldsymbol{a}}_2} \otimes {\mathit{\boldsymbol{b}}_2}, \cdots , {\mathit{\boldsymbol{a}}_{{I_0}}} \otimes {\mathit{\boldsymbol{b}}_{{I_0}}}] $

(6)

PARAFAC模型最重要的一个特征就是CP分解具有唯一性, 此为采用该模型进行数据分析的首要条件。下面给出PARAFAC模型分解唯一性的重要定理。

定理1(PARAFAC模型的分解唯一性)^[17]设χ∈ℝ^{I₁×I₂×I₃}是一个如式(4)所示的CP分解张量, 其中, A∈ℝ^I₁×F, B∈^I₂×F, C∈^I₃×F。若满足:

$ k(\mathit{\boldsymbol{A}}) + k(\mathit{\boldsymbol{B}}) + k(\mathit{\boldsymbol{C}}) \ge 2F + 2 $

(7)

则矩阵A、B和C在存在尺度模糊和列模糊的条件下是唯一的(也称为本质唯一)。

2.2 张量构造

根据文献[9], 对于互不相关的零均值非平稳实信号, 源信号在t时刻的协方差矩阵可以表示为:

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{C}}_1} = \mathit{\boldsymbol{E}}\{ {\mathit{\boldsymbol{x}}_t}{\mathit{\boldsymbol{x}}_{t + {\tau _1}}}\} = \mathit{\boldsymbol{A}} \cdot {\mathit{\boldsymbol{D}}_1} \cdot \mathit{\boldsymbol{A}}}\\ {{\mathit{\boldsymbol{C}}_2} = \mathit{\boldsymbol{E}}\{ {\mathit{\boldsymbol{x}}_t}{\mathit{\boldsymbol{x}}_{t + {\tau _1}}}\} = \mathit{\boldsymbol{A}} \cdot {\mathit{\boldsymbol{D}}_1} \cdot \mathit{\boldsymbol{A}}}\\ \vdots \\ {{\mathit{\boldsymbol{C}}_k} = \mathit{\boldsymbol{E}}\{ {\mathit{\boldsymbol{x}}_i}{\mathit{\boldsymbol{x}}_{t + {\tau _k}}}\} = \mathit{\boldsymbol{A}} \cdot {\mathit{\boldsymbol{D}}_k} \cdot \mathit{\boldsymbol{A}}} \end{array} $

(8)

其中, D_k=E{s_ts_t+tk}是对角矩阵, k=1, 2, …, K。为简化运算, 首先忽略噪声的影响, 而需要解决的问题是在M < N条件下, 从C_k中求解出A, 即实现欠定情况下基于二阶统计特征的盲源分离。

将矩阵C₁, C₂, …, C_k按如下方式压缩为张量C∈${\mathbb{C}}$^M×M×K, 其中, C_ijk=(C_k)_ij, i=1, 2, …, M, j=1, 2, …, M, k=1, 2, …, K。通过D_kf=(D_k)_ff定义矩阵D_kf, 得到:

$ {\mathit{\boldsymbol{C}}_{ijk}} = \sum\limits_{f = 1}^F {{a_{if}}} {a_{jf}}{d_{kf}} $

(9)

将式(9)改写为:

$ \mathit{\boldsymbol{C}} = \sum\limits_{f = 1}^f {{\mathit{\boldsymbol{a}}_f}} \circ {\mathit{\boldsymbol{a}}_f} \circ {\mathit{\boldsymbol{d}}_f} $

(10)

其中, a_f和d_f分别代表矩阵A和D的列向量。

张量沿3个方向的切片展开为矩阵:C₁∈${\mathbb{C}}$^M2×K, C₂∈${\mathbb{C}}$^KM×M和C₃∈${\mathbb{C}}$^MK×M, 满足如下关系:

$ {({\mathit{\boldsymbol{C}}_1})_{(i - 1)M + j, k}} = {({\mathit{\boldsymbol{C}}_2})_{(k - 1)K + i, j}} = {({\mathit{\boldsymbol{C}}_3})_{(j - 1)M + k, i}} = {\mathit{\boldsymbol{C}}_{ijk}} $

(11)

根据式(11), 张量的3个切片展开矩阵可表示为:

$ {{\mathit{\boldsymbol{C}}_1} = (\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{A}}) \cdot {\mathit{\boldsymbol{D}}^{\rm{T}}}} $

(12)

$ {{\mathit{\boldsymbol{C}}_2} = (\mathit{\boldsymbol{D}} \odot \mathit{\boldsymbol{A}}) \cdot {\mathit{\boldsymbol{A}}^{\rm{T}}}} $

(13)

$ {{\mathit{\boldsymbol{C}}_3} = (\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{D}}) \cdot {\mathit{\boldsymbol{A}}^{\rm{T}}}} $

(14)

由定理1可知, 满足一定条件时张量的CP分解唯一存在。文献[9]给出了传感器数目和源信号数目使CP分解唯一的关系, 即观测数M为2~8时, 所允许的源信号最大数目N_max分别为2、4、6、10、15、20、26。至此, 可将欠定混合矩阵的估计问题转化为三阶张量的CP分解问题。

2.3 基于CP分解的混合矩阵估计算法

作为解决PARAFAC模型的经典算法, ALS算法通过交替最小二乘误差γ来获得加载矩阵, 如式(15)所示:

$ \gamma = \left\| {{\mathit{\boldsymbol{C}}_1} - (\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{A}}){\mathit{\boldsymbol{D}}^{\rm{T}}}} \right\|_{\rm{F}}^2 $

(15)

其中, ‖·‖_F代表F-范数。矩阵A固定时, D的估计值可由最小二乘结果得到, 如式(16)所示:

$ \mathit{\boldsymbol{\hat D}} = {[{(\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{A}})^ + }{\mathit{\boldsymbol{C}}_1}]^{\rm{T}}} $

(16)

其中, (·)⁺表示广义逆矩阵。

对于式(13)与式(14), 通过同样的方式可得:

$ {\mathit{\boldsymbol{\hat A}} = {{[{{(\mathit{\boldsymbol{D}} \odot \mathit{\boldsymbol{A}})}^ + }{\mathit{\boldsymbol{C}}_2}]}^{\rm{T}}}} $

(17)

$ {\mathit{\boldsymbol{\hat A}} = {{[{{(\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{D}})}^ + }{\mathit{\boldsymbol{C}}_3}]}^{\rm{T}}}} $

(18)

由于初始矩阵的选择对ALS算法的准确性及收敛路径有较大影响, 因此本文通过直接三线性分解粗估加载矩阵, 将其作为ALS算法初始矩阵。

2.3.1 直接三线性分解

直接三线性分解^[18]是一种非迭代的求解PARAFAC模型的三维分解方法, 但其精度和可靠性比ALS算法要差, 在基于PARAFAC模型的ALS算法中可以用于粗估初始矩阵。算法具体步骤如下:

步骤1  根据式(13)~式(14), 分别对C₁、C₂和C₃进行奇异值分解, 其中, U、V取前F个奇异值矢量, W取前两个左奇异值矢量, 下标I、J、K分别表示张量的行、列、通道, R为信源数目。

$ {{\mathit{\boldsymbol{C}}_1} = {\mathit{\boldsymbol{U}}_I}{\mathit{\boldsymbol{S}}_I}{\mathit{\boldsymbol{V}}_I}, \mathit{\boldsymbol{U}} = {\mathit{\boldsymbol{U}}_I}(I \times R)} $

(19)

$ {{\mathit{\boldsymbol{C}}_2} = {\mathit{\boldsymbol{U}}_J}{\mathit{\boldsymbol{S}}_J}\mathit{\boldsymbol{V}}, \mathit{\boldsymbol{V}} = {\mathit{\boldsymbol{U}}_J}(J \times R)} $

(20)

$ {{\mathit{\boldsymbol{C}}_3} = {\mathit{\boldsymbol{U}}_K}{\mathit{\boldsymbol{S}}_K}{\mathit{\boldsymbol{V}}_K}, \mathit{\boldsymbol{W}} = {\mathit{\boldsymbol{U}}_K}(K \times 2)} $

(21)

步骤2  利用式(22)、式(23)构造样本矩阵G₁、G₂:

$ {{\mathit{\boldsymbol{G}}_1} = \sum\limits_{k = 1}^K {{w_{k1}}} {\mathit{\boldsymbol{U}}^{\rm{T}}}{\mathit{\boldsymbol{C}}_k}\mathit{\boldsymbol{V}}} $

(22)

$ {{\mathit{\boldsymbol{G}}_1} = \sum\limits_{k = 1}^K {{w_{k2}}} {\mathit{\boldsymbol{U}}^{\rm{T}}}{\mathit{\boldsymbol{C}}_k}\mathit{\boldsymbol{V}}} $

(23)

其中, w_ki为矩阵W中的元素, C_k∈ℝ^M×M为张量C的第k个切片矩阵。

步骤3  使用QZ分解G₁、G₂组成的特征方程, 令L、M分别为方程的特征向量, 矩阵A和D分别可由式(24)、式(25)得到^[13]:

$ {{\mathit{\boldsymbol{G}}_1}\mathit{\boldsymbol{L}}{\lambda _{2r}} = {\mathit{\boldsymbol{G}}_2}\mathit{\boldsymbol{L}}{\lambda _{1r}}, {\mathit{\boldsymbol{A}}_0} = \mathit{\boldsymbol{U}}{\mathit{\boldsymbol{L}}^{ - 1}}} $

(24)

$ {{\mathit{\boldsymbol{G}}_1}\mathit{\boldsymbol{M}}{\lambda _{2r}} = {\mathit{\boldsymbol{G}}_2}\mathit{\boldsymbol{M}}{\lambda _{1r}}, {\mathit{\boldsymbol{A}}_0} = \mathit{\boldsymbol{V}}{\mathit{\boldsymbol{M}}^{ - 1}}} $

(25)

至此, 通过直接三线性分解得到粗估的加载矩阵, 将其作为ALS算法的初始矩阵, 开始迭代。

2.3.2 ALS算法

ALS算法通过交替最小二乘误差来获得加载矩阵, 算法具体步骤如下:

步骤1  随机确定初始矩阵A^it-1、D^it-1, 利用直接三线性分解算法粗估加载矩阵, 其中, it表示迭代次数。

步骤2  利用式(16)和A^it-1估计D^it, 利用式(17)和A^it-1、D^it估计A₀^it, 利用式(18)和D^it、A₀^it估计A^it。

步骤3  利用式(15)和A^it、D^it计算误差函数γ^it。如果|γ^it-γ^it-1|>ξ, 则执行it=it+1并跳转到步骤2;否则, Â=A^it, $\mathit{\boldsymbol{\hat D}} = {\mathit{\boldsymbol{D}}^{{\rm{it}}}}$。

ALS算法对初值敏感且容易陷入局部极值, 同时运算时间较长, 因此, 本文使用标准线搜索加速收敛。

2.3.3 算法流程

在t次迭代, 利用算法预测一定的迭代步长ρ, 对矩阵A进行更新, 可以表示为:A=A+ρΔA。

本文采用的线加速方案只对矩阵D进行线搜索:

$ {\mathit{\boldsymbol{D}}^{{\rm{new}}}} = {\mathit{\boldsymbol{D}}^{it - 2}} + {\mathit{\boldsymbol{R}}_{{\rm{LS}}}}({\mathit{\boldsymbol{D}}^{it - 1}} - {\mathit{\boldsymbol{D}}^{it - 2}}) $

(26)

针对矩阵A, 利用式(18)采用最小二乘法做进一步估计。

算法步长的计算应较为简单, 如果它比相应的迭代需要更多的时间则没有意义。最简单的情况是R_LS被赋予固定值(在1.2和1.3之间)或者被设置为it^1/3[19]。本文中R_LS设置为it^1/3。改进ALS算法的具体步骤如下:

步骤1  利用直接三线性分解算法粗估加载矩阵A₀、D₀。

步骤2  利用式(16)和A₀估计D₁, 利用式(16)和A₀、D₁估计A₁, 令it=2, 则A^it-2=A₀, D^it-2=D₀, A^it-1=A₁, D^it-1=D₁。

步骤3  根据线加速公式(式(26))计算D^new, 利用式(18)和A^it-1、D^it-1计算A^new。

步骤4  利用式(15)和A^new、D^new、A^it-1、D^it-1计算误差函数γ^new和γ^it-1。

步骤5  若γ^new < γ^it-1, 则A^it-1=A^new, D^it-1=D^new, γ^it-1=γ^new; 否则, 直接执行步骤6。

步骤6  利用式(16)、式(17)和A^it-1、D^it-1计算A^it和D^it:

$ {\mathit{\boldsymbol{\hat D}} = {{[{{({\mathit{\boldsymbol{A}}^{it - 1}} \odot {\mathit{\boldsymbol{A}}^{it - 1}})}^ + }{\mathit{\boldsymbol{C}}_1}]}^{\rm{T}}}} $

(27)

$ {\mathit{\boldsymbol{\hat A}} = {{[{{({\mathit{\boldsymbol{D}}^{it}} \odot {\mathit{\boldsymbol{A}}^{it - 1}})}^ + }{\mathit{\boldsymbol{C}}_2}]}^{\rm{T}}}} $

(28)

步骤7  利用式(15)和A^it、D^it计算误差函数γ^it。如果|γ^it-γ^it-1|>ξ, 则执行it=it+1, 并跳转到步骤3;否则, Â=A^it, D=D^it。

3 源信号恢复

文献[14]提出了基于子空间投影法的源信号恢复算法, 在实际源信号恢复过程中设定一个与跳频信号信源功率有关的门限, 将其与混合矩阵的正交投影作比较来判断任意时频点同时存在的源信号个数是否小于阵元数, 并通过对离散噪点的进一步剔除, 达到更好的分离效果。跳频源信号的相对功率定义为:

$ P = \sum\limits_{i = 1}^{{L_k}} {{{\left\| {\mathit{\boldsymbol{X}}({t_{k, i}}, {f_{k, i}})} \right\|}_2}} /({L_k}\left\| {{\mathit{\boldsymbol{a}}_n}} \right\|) $

(29)

其中, ‖·‖₂代表二范数, L_k表示源信号的时频单源点集合所含有的向量个数。

假设时频点(t′, f′)上对应跳频源信号的混合矩阵列向量为A_R=[a_n1, a_n2, …, a_nR], R为每个时频点上同时存在的信源个数, 则在该时频点上的信号模型为:

$ \mathit{\boldsymbol{X}}({t^\prime }, {f^\prime }) = {\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{S}}_R}({t^\prime }, {f^\prime }) + \mathit{\boldsymbol{V}}({t^\prime }, {f^\prime }) $

(30)

其中, S_R(t′, f′)=[S_n1(t′, f′), S_n2(t′, f′), …, S_nR(t′, f′)]。

当时频点上的源信号个数小于观测阵元个数时, 混合矩阵列向量A_R的正交投影矩阵H为:

$ \mathit{\boldsymbol{H}} = \mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{A}}_R}{(\mathit{\boldsymbol{A}}_R^{\rm{H}}{\mathit{\boldsymbol{A}}_R})^{ - 1}}\mathit{\boldsymbol{A}}_R^{\rm{H}} $

(31)

其中, I为单位矩阵。由此可以得到:

$ \{ {a_{n1}}, {a_{n2}}, \cdots , {a_{nR}}\} = \mathop {{\rm{ argmin }}}\limits_{n1, n2, \cdots , nR} \{ \left\| {\mathit{\boldsymbol{HX}}({t^\prime }, {f^\prime })} \right\|{\mathit{\boldsymbol{A}}_R}\} $

(32)

将式(32)代入式(30), 在A_R已知的情况下可以计算出源信号S(t′, f′), 从而实现时频点上源信号个数小于接收阵元个数时的盲源分离。但当时频点(t′, f′)上同时存在的源信号个数等于接收阵元的个数M′时, 则不能按此来计算, 而是需要利用与源信号相对功率偏差相关的门限ε来判断何种情况下源信号数目等于接收阵元的个数。本文将门限设定为相对功率的0.2倍。

当R=M′时, 利用式(33)估计源信号S(t′, f′), 实现盲源分离:

$ \mathit{\boldsymbol{\hat S}}({t^\prime }, {f^\prime }) = {\mathit{\boldsymbol{A}}^{ - 1}}\mathit{\boldsymbol{X}}({t^\prime }, {f^\prime }) $

(33)

源信号恢复算法的具体步骤如下:

步骤1  设R=M-1, 基于式(30)计算式(31)。

步骤2  设定门限ε, 计算$\mathop {\min }\limits_{{\mathit{\boldsymbol{A}}_R}} \left\{ {\left\| {\mathit{\boldsymbol{HX}}\left( {{t^\prime }, {f^\prime }} \right)} \right\|{\mathit{\boldsymbol{A}}_R}} \right\}$。

步骤3  若$\mathop {\min }\limits_{{\mathit{\boldsymbol{A}}_R}} \left\{ {\left\| {\mathit{\boldsymbol{HX}}\left( {{t^\prime }, {f^\prime }} \right)} \right\|{\mathit{\boldsymbol{A}}_R}} \right\}$≤ε, 则利用式(30)计算估计跳频源信号; 否则, 利用式(33)计算估计跳频源信号。

4 仿真与结果分析

设定接收阵元数目为M=3, 阵元间距为2.5 m, 接收机处理带宽为[1 MHz, 50 MHz], 采样频率为100 MHz。跳频信号的调制方式均为BPSK, 频率集与跳频周期如表 1所示。

跳频源信号的时域波形图和STFT时频图如图 2和图 3所示, 混合跳频信号的时域波形图如图 4所示, 混合加噪跳频信号的STFT时频图如图 5所示。

4.1 混合矩阵估计

混合矩阵的估计结果将直接影响后续源信号的恢复效果。本文采用均方误差评价混合矩阵的估计性能, 如式(34)所示:

$ {E_1} = 10 \times {\rm{lg}}\left( {\frac{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {{{({{\hat a}_{ij}} - {a_{ij}})}^2}} } }}{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {a_{ij}^2} } }}} \right) $

(34)

其中, a_ij表示混合矩阵A中第i行第j列元素, â_ij表示估计得到的混合矩阵Â中第i行第j列元素。E₁的值越小, 表示混合矩阵的估计精度越高。

文献[1]采用K均值聚类算法迭代估计出混合矩阵, 文献[2]采用ALS算法进行估计, 文献[3]采用的AP聚类方法将每个样本都作为可能的聚类中心, 最后通过迭代确定聚类中心。本文算法先用非迭代的方法粗估混合矩阵, 再进行迭代, 确定混合矩阵。4种算法均方误差随信噪比的变化曲线如图 6所示。可以看出:文献[1]直接采用K均值聚类算法, 当接收阵元信噪比达到5 dB时, E₁极限值约为-15.5 dB, 与其他算法相比估计准确性较差; 在低信噪比条件下, 本文算法性能明显优于文献[2]算法, 与文献[3]算法接近; 在高信噪比条件下, 本文算法性能明显优于文献[3]算法, 略优于文献[2]算法。

4.2 同步组网跳频信号盲源分离

图 7和图 8分别给出了源信号数目为4、阵元数为3、接收阵元信噪比为3 dB时, 文献[20]基于聚类与子空间投影法的盲源分离结果与本文算法盲源分离结果的STFT时频图。可以看出, 本文算法分选出的跳频源信号时频图较文献[20]算法噪点更少, 跳频信号的时频图中能量更集中, 说明本文算法能更有效地实现盲源分离。

本文采用均方误差衡量多跳频信号网台分选的性能, 如式(35)所示:

$ {E_2} = 10 \times {\rm{lg}}\left( {\frac{{\sum\limits_{i = 1}^N E \{ {{({\mathit{\boldsymbol{s}}_i}(t) - {{\mathit{\boldsymbol{\hat s}}}_i}(t))}^2}\} }}{{\sum\limits_{i = 1}^N E \{ \mathit{\boldsymbol{s}}_i^2(t)\} }}} \right) $

(35)

其中, ŝ_i(t)为源信号s_i(t)的估计值。E₂的值越小, 表明分离效果越好, 分离后的信号与源信号越接近。在源信号保持不变的情况下, N=4, 阵元数目M=3、4、5分别对应欠定、正定、超定3种情况, 信噪比以5 dB步长从-5 dB递增至35 dB, 每点进行100次蒙特卡洛实验, 得到本文算法与文献[20]算法的信干比变化曲线, 如图 9所示, 可以看出:在低信噪比条件下, 文献[20]算法略优于欠定条件下本文算法, 但是在正定和超定条件下本文算法要明显优于文献[20]算法; 在高信噪比条件下, 本文算法在欠定、正定及超定条件下均优于文献[20]算法, 均方误差提升约1 dB, 本文算法在M=3、4、5时分选性能依次增加, 低信噪比条件下分离性能增加明显, 高信噪比条件下性能接近。

5 结束语

在现代战争中, 争夺电磁频谱使用和控制权的军事斗争日益激烈, 通信组网趋向密集和频繁。作为抗干扰通信的重要方式, 欠定条件下同步组网跳频信号的盲源分离具有重要意义。本文将欠定盲源分离中混合矩阵的估计问题转化为观测信号统计量所组成张量的CP分解问题, 并通过改进的ALS算法进行求解, 以满足稀疏性要求。在源信号恢复时采用子空间投影法, 并通过对离散噪点的进一步剔除, 达到更好的分离效果。仿真结果表明本文方法是一种有效的欠定盲源分离算法, 但其在建立跳频信号混合模型中采用了最简单的线性瞬时混合模型, 适用范围有限。后续将考虑时延、多径干扰等更多影响因素, 建立适用范围更广的混合模型。此外, 本文方法虽然降低了运算的迭代次数, 但是算法复杂度仍较高, 运算时间较长, 下一步对此进行改进。