CN108181611B - 基于子空间的压缩感知高分辨阵列处理方法 - Google Patents

Abstract

本发明公开一种基于子空间的压缩感知高分辨阵列处理方法，利用被测量信号在稀疏基下的稀疏性，构造凸优化计算函数。本发明能够在信号和噪声相干条件下识别信号，不同的信噪比条件下，本发明的均方根误差小于其它对比算法，且计算时间与现有方法相差无几。

Description

基于子空间的压缩感知高分辨阵列处理方法

技术领域

本发明涉及一种水声信号处理技术，具体涉及一种基于子空间的压缩感知高分辨阵列处理方法。

背景技术

浅海海洋中的目标检测和定位问题是海洋科学研究和技术应用的基本问题。在浅海海域中，来自不同的信号源的信号在传播过程中，信号幅度的衰弱，海洋中的噪声和信号之间的相干性对不同的信号源的检测和定位造成了影响，因此在浅海中需要阵列信号处理技术抑制干扰信号来实现对声信号的分离和定位。

波达方向(Direction of Arrival,DOA)估计是阵列信号处理领域的关键问题，在雷达、通讯、地震等众多领域有着广泛的应用。在最早的阵列信号DOA估计算法中最具代表性的当属波束成形(CBF)，而后又出现了一些以Capon算法为代表的的高分辨谱估计算法，这类算法较之CBF大大提高了分辨能力。上个世纪80年代中后期又相继出现了以多重信号分类(MUSIC)和旋转不变子空间(ESPRIT)为代表的的子空间类算法，它们主要是根据阵列接收信号的统计特性来估计目标的到达角，因此需要大量的独立同分布的测量数据，而且由于信号子空间与噪声子空间相互渗透，不能对信号源的方位进行有效测定。

声波在海洋中传播时，由于海水分层介质的折射和海面、海底的反射，声源和传感器阵列之间存在多个声传播路径的多径效应，并且每个声信号线路径均为声源信号的复本，相互之间为相干信号。在多个可以发出不相关的声信号的源的条件下，每一个源的声信号均存在多径效应现象。对于每一个源的多径中，只有直达线路径能够代表该声源方向，而其余线路径会对声源定位产生干扰，因此多径效应降低了声源定位的精度。在此情况下，对于声源定位而言，直达线路经的声信号可以作为声源定位算法中的有效信号，经过折射和反射的声信号复本可以理解为和信号相干的噪声。如图2所示为一个信号源和多个声线路径，其中黑色直线代表直达线路径，其余短划线为在海底和海面反射的线路径。直达线路径与阵列的夹角即为信号源的方位，但其余的线路径会对信号源定位产生干扰。

虽然常规波束形成算法，MUSIC算法等算法在水声信号处理，特别是声信号源定位方面有了较深入全面的理论和应用研究，但是这些算法的分辨力、实时性以及对海洋复杂噪声的鲁棒性仍然面临巨大的挑战。

压缩感知(Compressive Sensing,CS)在信号处理领域有了广泛而深入的研究，压缩感知理论提供了一种新的信号重构和参数估计方式，压缩感知理论指出如果信号是稀疏的或者是可压缩的，我们就可以利用信号的稀疏性，通过远远少于Nyquist采样定理规定的测量数据就可以对该稀疏信号进行采样并且恢复。国内外对于基于压缩感知的DOA估计已经有了一定的研究。Donoho等人指出稀疏信号具有受限等距条件(RIP)时，可以将压缩感知中的范数问题转换成范数求解问题。Tropp等人利用字典矩阵的原子集合的累积相关性，通过迭代运算求解稀疏向量，也就是经典的正交匹配追踪算法(OMP)。近年来，在水声学背景下，基于压缩感知的DOA估计有了较快的发展。Gerstoft等人利用了最小绝对值收敛和选择算子(LASSO)和最大后验相结合的方法提出了浅海环境下的基于压缩感知的DOA估计算法。Gurbuz等人提出了在未知信号数目的情况下的利用压缩感知的DOA估计算法。Northardt等人将压缩感知与预期可能性相结合提出DOA估计算法。

假设时域中分别有信号和噪声s(t)和n(t)，其频域表达式为S(υ)和N(υ)，那么信号和噪声的相关系数为：

当信号与噪声不相干时，上述表达式的分子为0，即S(υ)N(υ)^H＝0。

当有P个信号时，即有SN^H为零矩阵，同时有NS^H为零矩阵。

其中S＝[S₁,S₂,…,S_P]^T,N＝[n₁,n₂,…,n_P]^T。

反之，当SN^H或者NS^H不为零矩阵时，表示信号与噪声相干。

考虑信号重构问题，已知某一个测量矩阵

以及未知信号s在该矩阵下的线性测量值

满足：

y＝Φs

由于y的维数低于s的维数，该问题为非定方程，有无穷多个解。如果信号s是K稀疏的，即||s||₀＝K，并且Φ满足有限等距性质(RIP)条件，信号s的重构可以通过求解

范数问题来实现：

上式中||·||₀为求向量的

范数。

如果上述稀疏信号s可以在稀疏基Ψ下进行稀疏表示，即s＝Ψx，x为信号s在Ψ变换域的稀疏表示，则有

y＝Φs＝ΦΨx＝Ax

其中，A＝ΦΨ，称之为传感矩阵。

但

范数求解问题是NP难问题，无法在多项式时间内求解。而x在足够稀疏的条件下，

范数求解可以用

范数求解来代替，那么就有如下的计算表达式：

当考虑到噪声误差之后，上式可以改写成

其中，δ为与噪声相关的参数。

稀疏信号s的求解问题就转化为凸优化问题，但该算法对于浅海环境下，信号与噪声存在相干性时的DOA估计存在分辨力不够和对噪声的鲁棒性差等问题。

发明内容

发明目的：本发明的目的在于解决现有技术中存在的不足，提供一种基于子空间的压缩感知高分辨阵列处理方法。

技术方案：

、一种基于子空间的压缩感知高分辨阵列处理方法，其特征在于：包括以下步骤：

步骤1)、假设有P个信号源，并且基阵为M个传感器构成的线性等距阵列，基阵中每一个阵元的输出表示为P个入射信号的线性组合，如果用s_p(t)表示第p个信号源发射的信号，则阵列中第m个阵元上的测量输出y_m(t)表示为

其中，p＝1,2,…,P，y_m(t)是第m个传感器的接收信号，g_p(θ_p)是第p个信号在第m个传感器上的增益，n_m(t)为干扰信号和背景噪声在第m个传感器上的构成的附加噪声，τ_mp为第p个信号在第m个传感器和参考传感器间的传播时间延迟；

将增益归一化后，在频域中，表达式(1)改写成：

其中，υ表示频率，S_p(υ)为s_p(t)的傅里叶变换，Y_m(υ)为y_m(t)的傅里叶变换，N_m(υ)为n_m(t)的傅里叶变换，

d为相邻传感器间的距离，c为波速，θ_p为第p个信号源方向和参考传感器之间的夹角；

表达式(2)描述的为第m个阵元，即单个阵元；对基阵接收端M个阵元的输出信号的采样是同时进行的，因此所有M个阵元上的输出均写成一个M维的向量的表达式：

Y＝GS+N (3)

其中，Y(υ)＝[Y₁(υ),Y₂(υ),…,Y_M(υ)]^T，Y(υ)表示了M个传感器上接收信号的频域表达式；G＝[G₁,G₂,…,G_M]^T，

N(υ)＝[n₁(υ),n₂(υ),…,n_M(υ)]^T；S＝[S₁(υ),S₂(υ),…,S_P(υ)]^T，p＝1,…,P。

步骤2)、接收信号的频谱矩阵为：

其中，Y为传感器接收数据，

E{·}表示期望，·^H表示共轭转置；

然后将频谱矩阵分解为信号子空间和噪声子空间之和，即：

步骤3)、考虑信号重构问题，已知某个测量矩阵

以及所求信号s在这个测量矩阵下的测量值

满足y＝Φs，若信号s在稀疏基Ψ下是稀疏的，即有s＝Ψx，x为稀疏系数向量，Ψ为稀疏基矩阵，x在足够稀疏条件下，即||x||₀＜＜N，那么稀疏系数向量x的求解通过如下的

范数求解得到：

其中A为传感矩阵；当考虑噪声之后，表达式(6)改写成：

利用被测量信号在稀疏基下的稀疏性，构造凸优化计算函数：

所求得的向量

的峰值即为不同线路径的波达方向，矩阵G′为相关的系数矩阵，向量R_V为接收数据的信号子空间的所有元素的一维有序排列，即：

其中，

是

的一维有序排列，

为

的第i行第j列元素，

为接收数据的信号子空间，

所述步骤2)中对接收数据的频谱矩阵R进行特征值分解后，由较大的P个特征值和对应的特征向量构成矩阵

计算如下：

则信号子空间为：

为

的特征值构成的斜对角矩阵，U为(M×M)×(M×M)大小的单位方阵，该方阵的每一列μ₁,...,μ_M均为

的特征值λ₁,...,λ_M分别对应的特征向量，并且有λ₁≥λ₂≥…λ_M≥0；.^*表示共轭转置。

所述步骤2)中，根据公式(3)Y＝GS+N，可得：

E{YY^H}＝E{GSS^HG^H}+E{GSN^H}+E{NS^HG^H}+E{NN^H}，E{·}表示期望；

其中，E{GSS^HG^H}为纯信号谱矩阵，E{GSN^H}+E{NS^HG^H}为信号和噪声互谱矩阵，E{NN^H}为噪声谱矩阵。

由于频谱矩阵由信号子空间和噪声子空间构成，因此：

假设信号与信号之间不相关，即

那么有：

假设W为GSS^HG^H的一维有序排列，即有：

其中，

另外，

假设U为GSN^H+NS^HG^H的一维有序排列：

即：

将表达式(9)～(14)相结合可以得到：

R_V＝G′P+U

与表达式(7)相结合，

可得：

有益效果：与现有技术相比，本发明的分辨率得到极大的提高。传统方法中，5个分别相差5°的信号到达阵列时，在信号与噪声的相干性较强，即较低的信噪比(-17dB)下无法识别出所有的信号，而本发明方法可以识别出所有的信号。在不同的信噪比条件下，本发明方法的均方根误差小于其它对比算法，并且计算时间与现有方法相差无几。

附图说明

图1是本发明和现有技术在信噪比-17dB下的对比示意图；

图2为现有技术中浅海中声信号多径传播示意图。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

假设有P个信号源，并且基阵为M个传感器构成的线性等距阵列，基阵中每一个阵元的输出可以表示为P个入射信号的线性组合。如果用s_p(t)表示第p个(p＝1,2,…,P)信号源发射的信号，则阵列中第m个阵元上的测量输出y_m(t)可以表示为

y_m(t)是第m个传感器的接收信号，g_p(θ_p)是第p个信号在第m个传感器上的增益，n_m(t)为干扰信号和背景噪声在第m个传感器上的构成的附加噪声，τ_mp为第p个信号在第m个传感器和参考传感器间的传播时间延迟。

将增益归一化后，在频域中，表达式(1)可以改写成：

其中，

d为相邻传感器间的距离，c为波速，θ_p为第p个信号源方向和参考传感器之间的夹角。

表达式(2)描述的为第m个阵元，即单个阵元。对基阵接收端M个阵元的输出信号的采样是同时进行的，因此所有M个阵元上的输出可以写成一个M维的向量的表达式：

Y＝GS+N (3)

其中，Y(υ)＝[Y₁(υ),Y₂(υ),…,Y_M(υ)]^T，Y(υ)表示了M个传感器上接收信号的频域表达式；G＝[G₁,G₂,…,G_M]^T，

N(υ)＝[n₁(υ),n₂(υ),…,n_M(υ)]^T；S＝[S₁(υ),S₂(υ),…,S_P(υ)]^T，p＝1,…,P。

考虑信号重构问题，已知某个测量矩阵

以及所求信号s在这个测量矩阵下的测量值

满足y＝Φs，如果信号s在稀疏基Ψ下是稀疏的，即有s＝Ψx，x为稀疏系数向量，Ψ为稀疏基矩阵，x在足够稀疏条件下，即||x||₀＜＜N，那么稀疏系数向量x的求解可以通过如下的

范数求解得到：

当考虑噪声之后，表达式(4)改写成：

本发明中，信号s为空域内稀疏信号。则接收信号的频谱矩阵为：

其中，

E{·}表示期望，·^H表示共轭转置。

频谱矩阵可以分解为信号子空间和噪声子空间的和。

对R进行特征值分解：

其中

为

的特征值构成的斜对角矩阵，U为(M×M)×(M×M)大小的单位方阵，其中该方阵的每一列μ₁,...,μ_M均为

的特征值λ₁,...,λ_M分别对应的特征向量。并且有

λ₁≥λ₂≥…λ_M≥0

此时，信号子空间为：

即取最大的P个特征值和这P个特征值对应的特征向量构成信号子空间

.^*表示共轭转置。

将

中的每一个元素按照顺序进行前后排列得到:

其中，

是

的一维有序排列，

为

的第i行第j列元素，

由表达式(3)得到：

E{YY^H}＝E{GSS^HG^H}+E{GSN^H}+E{NS^HG^H}+E{NN^H} (11)

其中，E{·}表示期望。

其中，E{GSS^HG^H}为纯信号谱矩阵，E{GSN^H}+E{NS^HG^H}为信号和噪声互谱矩阵，E{NN^H}为噪声谱矩阵。

因为频谱矩阵由信号子空间和噪声子空间构成，所以有：

假设信号与信号之间不相关，即

那么有：

假设W为GSS^HG^H的一维有序排列，即有：

其中，

另外，

假设U为GSN^H+NS^HG^H的一维有序排列：

即：

将表达式(12)～(17)相结合可以得到：

R_V＝G′P+U

上式与表达式(5)相结合后，本发明算法可以通过如下的计算式实现：

上述整体过程中，通过将频谱矩阵R分成信号子空间和纯噪声子空间，并且计算只使用信号子空间来实现(也就是表达式(7)～(9))使得分辨率和鲁棒性得到提高和优化。

实施例1：

本实施例的具体实验条件描述如下：

来自5个不同信号源的远场信号被由20个传感器构成的线性等距阵列接收，该5个信号的理论达到角度分别为5°，10°，15°，20°和25°，每个传感器的采样次数为200次。

传感器阵列上的附加噪声由这5个信号的相干信号和高斯噪声(背景噪声)叠加构成。每个信号源除了发射出直达阵列的信号之外，还发射出另外4个相干信号(属于附加噪声)。假设直达信号的到达角为θ_i，那么这4个对应的相干信号的到达角分别为θ_i+1°,θ_i-1°,θ_i+2°,θ_i-2°。

如图1所示，本实施例中将本发明与几个现有技术进行对比，可见本发明能够识别出所有的信号源；并且，如表1和表2所示，在不同的信噪比条件下，本发明的均方根误差小于其它对比算法，并且计算时间与现有方法相差无几。

表1本发明和现有技术的均方根误差的对比

其中，表1的对比方案中所有的方法均在同一个信噪比条件下均运行50次，然后求在同一个信噪比条件下，每个算法的50次实验的均方根误差值。

表2本发明和现有技术的计算时间对比

Claims (3)

1.一种基于子空间的压缩感知高分辨阵列处理方法，其特征在于：包括以下步骤：

步骤1)、假设有P个信号源，并且基阵为M个传感器构成的线性等距阵列，基阵中每一个传感器的输出表示为P个入射信号的线性组合，如果用S_p(t)表示第p个信号源发射的信号，则阵列中第m个传感器上的测量输出y_m(t)表示为

其中，p＝1,2,…,P，y_m(t)是第m个传感器的接收信号，g_p(θ_p)是第p个信号在第m个传感器上的增益，n_m(t)为干扰信号和背景噪声在第m个传感器上的构成的附加噪声，τ_mp为第p个信号在第m个传感器和参考传感器间的传播时间延迟；

将增益归一化后，在频域中，表达式(1)改写成：

其中，υ表示频率，S_p(υ)为S_p(t)的傅里叶变换，Y_m(υ)为y_m(t)的傅里叶变换，N_m(υ)为n_m(t)的傅里叶变换，

d为相邻传感器间的距离，c为波速，θ_p为第p个信号源方向和参考传感器之间的夹角，τ_p为第p个信号在相邻传感器间的传播时间延迟；

表达式(2)描述的为第m个传感器，即单个传感器；对基阵接收端M个传感器的输出信号的采样是同时进行的，因此所有M个传感器上的输出均写成一个M维的向量的表达式：

Y＝GS+N (3)

其中，Y＝[Y₁(υ)，Y₂(υ)，…，Y_M(υ)]^T，Y表示了M个传感器上接收信号的频域表达式；G＝[G₁,G₂,…,G_M]^T，

N＝[N₁(υ)，N₂(υ)，…，N_M(υ)]^T；S＝[S₁(υ),S₂(υ),…,S_P(υ)]^T，p＝1,…,P

步骤2)、接收信号的频谱矩阵为：

其中，

E{·}表示期望，·^H表示共轭转置；

然后将频谱矩阵分解为信号子空间和噪声子空间之和，即：

其中，

为频谱矩阵的信号子空间部分，

为频谱矩阵的噪声子空间部分；

步骤3)、考虑信号重构问题，已知某个测量矩阵

以及所求信号s在这个测量矩阵下的测量值

满足y＝Φs，若信号s在稀疏基Ψ下是稀疏的，即有s＝Ψx，x为稀疏系数向量，Ψ为稀疏基矩阵，x在足够稀疏条件下，即||x||₀＜＜P，其中P为矩阵列的数量，那么稀疏系数向量x的求解通过如下的l₁范数求解得到：

其中A为传感矩阵；当考虑噪声之后，表达式(6)改写成：

δ为噪声二范数的上界；

其中，利用被测量信号在稀疏基下的稀疏性，构造凸优化计算函数：

其中

所求得的向量

的峰值即为不同线路径的波达方向，矩阵G′为相关的系数矩阵，向量R_V为接收数据的信号子空间的所有元素的一维有序排列，即：

其中，

是

的一维有序排列，

为

的第i行第j列元素，

2.根据权利要求1所述的基于子空间的压缩感知高分辨阵列处理方法，其特征在于：步骤2)中对接收数据的频谱矩阵R进行特征值分解后，由较大的P个特征值和对应的特征向量构成矩阵

计算如下：

则信号子空间为：

为

的特征值构成的斜对角矩阵，U为(M×M)×(M×M)大小的单位方阵，该方阵的每一列μ₁,…,μ_M均为

的特征值λ₁,...,λ_M分别对应的特征向量，并且有λ₁≥λ₂≥…λ_M≥0；.*表示共轭转置。

3.根据权利要求1所述的基于子空间的压缩感知高分辨阵列处理方法，其特征在于：步骤2)中，根据公式(3)Y＝GS+N，可得：

E{YY^H}＝E{GSS^HG^H}+Ε{GSN^H}+E{NS^HG^H}+E{NN^H}，E{·}表示期望；

其中，E{GSS^HG^H}为纯信号谱矩阵，Ε{GSN^H}+E{NS^HG^H}为信号和噪声互谱矩阵，E{NN^H}为噪声谱矩阵；

由于频谱矩阵由信号子空间和噪声子空间构成，因此：

假设信号与信号之间不相关，即

那么有：

假设W为GSS^HG^H的一维有序排列，即有：

其中，

另外，

假设U为GSN^H+NS^HG^H的一维有序排列：

即：

将表达式(9)～(14)相结合得到：R_V＝G′P+U

再加上

可得：

Images