CN109541572B - 一种基于线性环境噪声模型的子空间方位估计方法 - Google Patents

Abstract

本发明涉及一种基于线性环境噪声模型的子空间方位估计方法，对于实际的噪声环境，由于风成噪声和航船噪声的影响，其空间噪声强度分布是非均匀的，使用线性噪声模型进行建模，将空间噪声功率分布函数进行Fourier级数展开，利用有限项Fourier级数近似拟合环境噪声模型。提取出线性噪声模型条件下的信号子空间，利用信号子空间与噪声子空间的正交性得到目标的方位估计。

Description

一种基于线性环境噪声模型的子空间方位估计方法

技术领域

本发明属于阵列信号处理、声呐信号处理等领域，涉及一种基于线性环境噪声模型的子空间方位估计方法。

背景技术

基于传感器阵列的方位估计(Direction of Arrival,DOA)在声纳、雷达等领域得到广泛的应用。利用水平放置的水听器阵列对水下和水面多目标进行高分辨的方位估计是水下阵列信号处理的重要内容，也是进行水中目标探测和定位的关键技术。

声呐阵列放置在海洋环境中，阵列接收信号除了目标信号外还有海洋环境噪声。海洋环境噪声是海洋声场的背景场，对海洋环境噪声的研究可以降低或抑制环境噪声对水下设备的干扰，提高水下设备的性能和作用距离。海洋环境噪声的产生因素很多，通常包括潮夕、波浪所引起的压力波和湍流引起的压力脉动，以及地震活动、风动海面、降雨、分子热运动、海洋中生物群体的活动等。目前大多数方位估计方法均假设环境噪声为高斯白噪声，即噪声协方差矩阵可以写成噪声方差与单位矩阵的乘积形式。然而，实际各阵元接收的海洋环境噪声可能是相关的，并且由于风成噪声和远处航船噪声等因素的影响，环境噪声强度的空间分布也可能具有方向性，其在水平方向上是随着角度变化的慢变函数。因此，在复杂海洋环境条件下，白噪声假设的噪声模型可能造成严重的方位估计偏差。

目前主要的方位估计方法有常规波束形成算法和最小方差无失真波束形成算法等。常规波束形成算法是基于常规波束形成的扫描测向方法，具有结构简单计算量小的特点，并且环境噪声模型误差具有一定的稳健性，然而受到“瑞利准则”的限制，空间方位分辨能力较差；最小方差无失真波束形成算法利用CAPON自适应波束形成器实现空间扫描波束输出功率谱估计，空间方位分辨能力优于常规波束形成算法，然而该算法假设噪声模型为白噪声模型。S.D.Somasundaram研究了宽带稳健的最小方差无失真波束形成器设计方法(Somasundaram S D.Wideband Robust Capon Beamforming for Passive Sonar.IEEEJournal of Oceanic Engineering,2013,38(2):308-322.)，并将其应用于被动声纳的多目标探测，提高了目标的空间分辨能力，但其噪声模型仍假设为白噪声模型，对于复杂环境噪声条件的方位估计性能仍需检验。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于线性环境噪声模型的子空间方位估计方法，针对水平固定布放的水听器基阵，能够在复杂海洋环境噪声条件下实现高分辨水下目标方位估计。

技术方案

一种基于线性环境噪声模型的子空间方位估计方法，其特征在于：采用水平放置于水下的直线传感器基阵，称为水平阵列；阵元数为M，阵元间距为d，K个目标数分布在该水平阵列的远场范围内，估计步骤如下：

步骤1：x(n),n＝1,…,N为直线传感器基阵列的输出向量，计算阵列输出向量的采样协方差矩阵

其中：N是数据点数，上标“H”为共轭转置算子；

步骤2、估计模型阶数：首先在无水下探测目标K-0时，计算常规波束形成方位谱图：

其中：a(θ)为阵列流形向量a(θ)＝[1,exp(-j2πfτ),…,exp{-j2πf(M-1)τ}]^T，f为信号频率，τ＝dsinθ/c，d为阵元间距，c为声速，θ为空间方位角度；

对常规波束形成输出的方位谱CBF(θ)进行Fourier级数拟合，其拟合阶数即为环境噪声的建模阶数，记为J；

步骤3：计算模型参数矩阵

Γ＝[vec(∑₀),vec(∑₁),…,vec(∑_j),…,vec(∑_J)]

其中：

J＝2L+1，

步骤4：噪声协方差矩阵Σ的初始估计为

其中上标括号中的数字表示迭代次数；

步骤5：利用第i-1次迭代的噪声协方差矩阵Σ^(i-1)计算第i次迭代的噪声预白化处理后的协方差矩阵

步骤6：对

进行特征分解，其中

是维度为K×K的主特征值矩阵，

为主特征向量即主特征值对应的特征向量；所述主特征值矩阵为对角线元素为数值较大的前K个特征值；

步骤7：重构信号子空间矩阵

其中I为单位对角矩阵；

步骤8：对于信号协方差矩阵进行向量化运算

同时对采样协方差矩阵

进行向量化运算，即

矩阵向量化运算vec{}是指把矩阵的列向量按列组合，形成一个长的列向量

步骤9：计算第i次迭代时线性噪声模型的系数向量η⁽ⁱ⁾，

pinv{}表示伪逆运算算子；

步骤10：计算第i次迭代时的噪声协方差矩阵Σ⁽ⁱ⁾＝vec^-1{Γη⁽ⁱ⁾}，其中vec^-1{}表示矩阵向量化的逆运算；

步骤11：判断迭代终止条件为abs(L⁽ⁱ⁾-L^(i-1))≤10^-3是否成立，成立则进行下一步骤，不成立则返回步骤4，并且迭代序号i＝i+1；

其中代价函数

||表示矩阵的行列式的值，tr{}表示矩阵求迹运算，abs()表示求绝对值运算；

步骤12：得到最终的重构信号子空间矩阵

上标It表示最终的迭代次数；

步骤13：利用信号子空间与噪声子空间的正交性，得到目标的方位谱为P(θ)＝1/{a^H(θ)Πa(θ)}，其中

方位谱的最大值为目标的方位。

有益效果

本发明提出的一种基于线性环境噪声模型的子空间方位估计方法，对于实际的噪声环境，由于风成噪声和航船噪声的影响，其空间噪声强度分布是非均匀的，使用线性噪声模型进行建模，将空间噪声功率分布函数进行Fourier级数展开，利用有限项Fourier级数近似拟合环境噪声模型。提取出线性噪声模型条件下的信号子空间，利用信号子空间与噪声子空间的正交性得到目标的方位估计。

附图说明

图1是一种基于线性环境噪声模型的子空间方位估计方法的总体流程框图；

图2是方位谱图。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明解决其技术问题所采用的技术方案包括以下几个方面：

1：一个水平放置于水下的直线传感器基阵的阵元数为M，阵元间距为d，此时阵列的总长度为的d(M-1)；

利用该阵列进行水下或水面多目标的方位估计，设目标数为K，且均分布在该水平阵列的远场范围内，入射方位为θ＝(θ₁,…,θ_K)，其中θ_k表示第k个信号的入射方位，阵列输出向量表示为x(n)＝As(n)+e(n),n＝1,…,N，其中，N是数据点数，s(n)＝[s₁(n),…,s_K(n)]^T为入射信号波形，其中上标“T”表示向量转置运算，e(n)是阵元接收噪声，A是阵列流形矩阵(steering matrix)。对于水平放置的传感器阵列，阵列流形向量a(θ)＝[1,exp(-j2πfτ),…,exp{-j2πf(M-1)τ}]^T，其中f为信号频率，τ＝d sinθ/c，d为阵元间距，c为声速。数据协方差矩阵为R＝E{x(n)x^H(n)}＝APA^H+Σ，其中E{}表示数学期望算子，P为对角矩阵，对角线元素为信号功率，上标“H”为共轭转置算子，Σ为噪声协方差矩阵，即Σ＝E{e(n)e^H(n)}，实际应用中使用采样协方差矩阵

表示近似数据协方差矩阵R，即

2：线性噪声模型的建立

在给定采样时间n时，噪声强度为一个随机变量v(θ,n)。因此，阵元接收噪声波形为

假设噪声强度v(θ,n)空间上不相关，时间上为零均值高斯白噪声，则噪声协方差矩阵表示为

其中ε(θ)为空间噪声功率密度，空间噪声功率密度ε(θ)为入射方位角θ的周期函数，2π为一个周期，将空间噪声功率密度ε(θ)按照Fourier级数展开得到

其中Fourier系数

式中l为Fourier级数展开阶数，c_l和s_l对应不同阶数下Fourier系数；

使用有限阶数的Fourier级数对空间噪声功率密度ε(θ)近似表示，设该阶数为L，则数据协方差矩阵

其中噪声协方差矩阵为

其中J＝2L+1，

和

当L＝0时

并且当阵元间距为半波长时，Σ＝c₀I，其中I为单位对角矩阵，此时环境噪声等同于零均值高斯白噪声。对噪声协方差矩阵

进行矩阵向量化运算(矩阵向量化运算vec{}是指把矩阵的列向量按列组合，形成一个长的列向量)可得到vec{Σ}＝Γη，式中矩阵Γ＝[vec(Σ₀),vec(Σ₁),…,vec(Σ_J)]，向量η＝[η₁,…,η_J]。

3：信号子空间提取

信号子空间通过迭代方式提取。首先给出噪声协方差矩阵Σ的初始估计为

其中上标括号中的数字表示迭代次数。利用第i-1次迭代的噪声协方差矩阵Σ^(i-1)计算第i次迭代的相关参数：①计算噪声预白化处理后的协方差矩阵

②对

进行特征分解，其中

是维度为K×K的主特征值矩阵(对角线元素为数值较大的前K个特征值)，

为主特征向量(即主特征值对应的特征向量)。③重构信号子空间矩阵

④对于信号协方差矩阵进行向量化运算

同时对采样协方差矩阵

进行向量化运算，即

⑤计算第i次迭代时线性噪声模型的系数向量

其中

为线性噪声模型系数，J为模型阶数，

pinv{}表示伪逆运算算子。⑥计算第i次迭代时的噪声协方差矩阵Σ⁽ⁱ⁾＝vec^-1{Γη⁽ⁱ⁾}，其中vec^-1{}表示矩阵向量化的逆运算。⑦判断迭代终止条件为abs(L⁽ⁱ⁾-L^(i-1))≤10^-3是否成立，成立则进行下一步骤，不成立则返回第①步操作，并且迭代序号i＝i+1。其中代价函数：

||表示矩阵的行列式的值，tr{}表示矩阵求迹运算，abs()表示求绝对值运算。⑧最终的重构信号子空间矩阵为

上标I表示最终的迭代次数。

4：目标方位估计

利用信号子空间与噪声子空间的正交性，目标的方位谱可表示为P(θ)＝1/{a^H(θ)Πa(θ)}，其中

搜索方位谱的最大值，获得目标的方位。

具体步骤为：

(1)布置一条水平放置于水下的直线传感器基阵，阵元数为M，阵元间距为d。利用该阵列接收和记录水声信号，x(n)为阵列输出向量；

(2)计算阵列采样协方差矩阵

(3)估计模型阶数。首先在无水下探测目标时，计算常规波束形成方位谱图，即

其中a(θ)为阵列流形向量a(θ)＝[1,exp(-j2πfτ),…,exp{-j2πf(M-1)τ}]^T，f为信号频率，τ＝dsinθ/c，d为阵元间距，c为声速，θ为空间方位角度。对常规波束形成输出的方位谱进行Fourier级数拟合，其拟合阶数即为环境噪声的建模阶数，记为J；

(4)计算模型参数矩阵Γ＝[vec(Σ₀),vec(Σ₁),…,vec(Σ_J)]，其中J＝2L+1，

(5)噪声协方差矩阵Σ的初始估计为

其中上标括号中的数字表示迭代次数；

(6)利用第i-1次迭代的噪声协方差矩阵Σ^(i-1)计算第i次迭代的噪声预白化处理后的协方差矩阵

(7)对

进行特征分解，其中

是维度为K×K的主特征值矩阵(对角线元素为数值较大的前K个特征值)，

为主特征向量(即主特征值对应的特征向量)；

(8)重构信号子空间矩阵

(9)对于信号协方差矩阵进行向量化运算

同时对采样协方差矩阵

进行向量化运算，即

(10)计算第i次迭代时线性噪声模型的系数向量

其中

为线性噪声模型系数，J为模型阶数，

pinv{}表示伪逆运算算子；

(11)计算第i次迭代时的噪声协方差矩阵Σ⁽ⁱ⁾＝vec^-1{Γη⁽ⁱ⁾}，其中vec^-1{}表示矩阵向量化的逆运算；

(12)判断迭代终止条件为abs(L⁽ⁱ⁾-L^(i-1))≤10^-3是否成立，成立则进行下一步骤，不成立则返回步骤(5)，并且迭代序号i＝i+1。其中代价函数

||表示矩阵的行列式的值，tr{}表示矩阵求迹运算，abs()表示求绝对值运算。

(13)输出最终的重构信号子空间矩阵

(14)利用信号子空间与噪声子空间的正交性，目标的方位谱可表示为：

P(θ)＝1/{a^H(θ)Πa(θ)}，其中

搜索方位谱的最大值，获得目标的方位。

具体实施例：

利用一条10元水听器阵列接收远场三个点源信号，信号方位为55°、80°和90°，阵元间距为半波长，环境噪声为空间非均匀噪声环境，分别使用常规波束形成方法，MUSIC方法以及本发明给出的方法计算目标的方位谱图，如附图2所示，其中MUSIC方法对噪声环境较为敏感，在此环境下已经失效，常规波束形成方法分辨率较低，而本发明给出的方法可以准确给出三个目标的方位，并且对于空间非均匀噪声环境的影响最稳健。

Claims (1)

1.一种基于线性环境噪声模型的子空间方位估计方法，其特征在于：采用水平放置于水下的直线传感器基阵，称为水平阵列；阵元数为M，阵元间距为d，K个目标数分布在水平阵列的远场范围内，估计步骤如下：

步骤1：x(n),n＝1,…,N为直线传感器基阵列的输出向量，计算阵列输出向量的采样协方差矩阵

其中：N是数据点数，上标“H”为共轭转置算子；

步骤2、估计模型阶数：首先在无水下探测目标K-0时，计算常规波束形成方位谱图：

其中：a(θ)为阵列流形向量a(θ)＝[1,exp(-j2πfτ),…,exp{-j2πf(M-1)τ}]^T，f为信号频率，τ＝dsinθ/c，d为阵元间距，c为声速，θ为空间方位角度；

对常规波束形成输出的方位谱CBF(θ)进行Fourier级数拟合，其拟合阶数即为环境噪声的建模阶数，记为J；

步骤3：计算模型参数矩阵

Γ＝[vec(∑₀),vec(∑₁),…,vec(∑_j),…,vec(∑_J)]

其中：

J＝2L+1，

步骤4：噪声协方差矩阵Σ的初始估计为

其中上标括号中的数字表示迭代次数；

步骤5：利用第i-1次迭代的噪声协方差矩阵Σ^(i-1)计算第i次迭代的噪声预白化处理后的协方差矩阵

步骤6：对

进行特征分解，其中

是维度为K×K的主特征值矩阵，

为主特征向量即主特征值对应的特征向量；所述主特征值矩阵为对角线元素为数值较大的前K个特征值；

步骤7：重构信号子空间矩阵

其中I为单位对角矩阵；

步骤8：对于信号协方差矩阵进行向量化运算

同时对采样协方差矩阵

进行向量化运算，即

矩阵向量化运算vec{}是指把矩阵的列向量按列组合，形成一个长的列向量

步骤9：计算第i次迭代时线性噪声模型的系数向量η⁽ⁱ⁾，

pinv{}表示伪逆运算算子；

步骤10：计算第i次迭代时的噪声协方差矩阵Σ⁽ⁱ⁾＝vec^-1{Γη⁽ⁱ⁾}，其中vec^-1{}表示矩阵向量化的逆运算；

步骤11：判断迭代终止条件为abs(L⁽ⁱ⁾-L^(i-1))≤10^-3是否成立，成立则进行下一步骤，不成立则返回步骤4，并且迭代序号i＝i+1；

其中代价函数

||表示矩阵的行列式的值，tr{}表示矩阵求迹运算，abs()表示求绝对值运算；

步骤12：得到最终的重构信号子空间矩阵

上标It表示最终的迭代次数；

步骤13：利用信号子空间与噪声子空间的正交性，得到目标的方位谱为P(θ)＝1/{a^H(θ)Πa(θ)}，其中

方位谱的最大值为目标的方位。

Images