一种基于阵列信号数值计算的水下目标辐射噪声相干信源测向方法与流程

1.本发明属于水下目标辐射噪声测向定位技术领域，尤其涉及一种基于阵列信号数值计算的水下目标辐射噪声相干信源测向方法，适用于多途信道中相干信源的测向。


背景技术：

2.对水下目标辐射噪声进行测向定位具有重要意义。但由于水声信道的不均匀波导与有界空间结构特性，目标辐射噪声在海水中容易发生多途传播现象，从而形成多个相干信源。在众多测向定位算法中，常规波束形成(conditional beamforming,cbf)算法具有最好的稳健性，即对阵列安装误差和噪声特性变化包容性最好。但cbf算法受限于“傅里叶限”限制，无法分辨位于同一个波束宽度内的多个相干信源。若要提高cbf算法的方位分辨能力，需要扩大阵列有效孔径或阵元数目，但在水下机动平台上较难实施。其它测向算法，例如最小方差无畸变响应波束形成算法(minimum variance distortionless response，mvdr)、多重信号分类方法(multiple signal characteristics，music)和旋转不变子空间技术(estimation of signal parameters via rotational invariance techniques，esprit)等虽然具有更高的方位分辨能力，但它们普遍要求海洋背景噪声与目标辐射噪声不相关、多个目标辐射噪声之间不相关，并且要求背景噪声严格服从某些确定的分布，以及依赖阵列高精度安装和标定。在多途水声信道中，阵列接收的目标辐射噪声可能由多个相干信源构成，而且海洋背景噪声统计特性也十分复杂，往往受到洋流、航运和气象等时变空变因素影响。因此高分辨测向算法在实际应用中面临较大的局限性，最直接表现就是稳健性不足，远不如cbf方法。为了对多个相干信源进行测向，在这些高分辨测向算法的基础上产生了一些改进算法，即在测向之前采用空间平滑等方法做解相干预处理，但预处理会损失阵列孔径和方位分辨能力。除此之外，还有一些折中算法在方位分辨率和稳健性方面取得一定的平衡，例如波束域mvdr算法、对角加载算法和基于最差性能优化的稳健波束形成算法等。这些折中方法的方位分辨率和稳健性均介于cbf和其它高分辨算法之间。
3.与上述算法相比，以极大似然(maximum likelihood，ml)算法和加权子空间拟合(weight subspace fitting，wsf)为典型代表的空间拟合类算法性能更加突出：通过构造拟合关系得到一个需要进行多维求解的目标函数，目标方位估计精度高，尤其在低信噪比和小快拍数据的情况下估计性能优良，而且不需要进行解相干处理就可以对相关信号进行doa估计(安春莲.独立信号与相干信号并存的测向算法研究[d].哈尔滨工程大学.2013)。从性能上，空间拟合类算法十分适合多途条件下的水下目标辐射噪声的测向定位。但这类算法的不足在于目标函数求解过程涉及多维搜索，计算量繁重，很难满足工程应用中实时计算要求。


技术实现要素：

[0004]
本发明针对空间拟合类算法存在的目标函数求解过程涉及多维搜索，计算量繁
重，很难满足工程应用中实时计算要求的问题，提出一种基于阵列信号数值计算的水下目标辐射噪声相干信源测向方法。
[0005]
为了实现上述目的，本发明采用以下技术方案：
[0006]
一种基于阵列信号数值计算的水下目标辐射噪声相干信源测向方法，包括：
[0007]
步骤1：生成矩阵和转换矩阵的非实时计算存储：预先将水下目标辐射噪声不同频率线谱的到达方向做数值化离散，然后依据不同信源数目分别做组合排列，最后计算每个组合下的生成矩阵和转换矩阵；
[0008]
步骤2：目标辐射噪声相干信源数目的实时估计：实时计算各个阵元通道采集的目标辐射噪声中线谱分量的样本协方差矩阵，然后通过盖氏圆信源估计方法估计目标辐射噪声相干信源数目；
[0009]
步骤3：目标辐射噪声相干信源到达方向分布范围的实时估计：计算各阵元通道采集的目标辐射噪声的线谱分量的dft值，然后通过常规波束形成方法估计目标辐射噪声相干信源到达方向分布范围；
[0010]
步骤4：目标辐射噪声相干信源到达方向的实时精确计算：依据目标辐射噪声相干信源数目估计结果和到达方向分布范围，确定目标辐射噪声相干信源到达方向的所有可能组合以及阵列信号数学模型；然后选取预先存储的生成矩阵、转换矩阵与频域阵列信号向量相乘，计算模型系数与模型偏差；最后将模型偏差最小时对应的到达方向组合确定为目标辐射噪声的相干信源真实到达方向，并将此时的模型系数作为相干信源强度。
[0011]
进一步地，所述步骤1包括：
[0012]
将水下目标辐射噪声的不同频率线谱到达方向做数值化离散，离散间隔为δθ分辨率，离散波达方向的数目为m；
[0013]
相干辐射噪声到达方向依据信源组合从1至l
‑
3做组合排列，组合个数分别为其中l为阵元数目；
[0014]
对于每个频率的每个组合，计算转换矩阵：
[0015]
d(f
i
，θ1，
…
，θ
m
)＝[g
h
(f
i
，θ1，
…
，θ
m
)
·
g(f
i
，θ1，
…
，θ
m
)]
‑1·
g
h
(f
i
，θ1，
…
，θ
m
)
[0016]
其中，1＜m≤m，g(f
i
，θ1，
…
，θ
m
)表示在线谱频率f
i
以及波达方向组合(θ1，
…
，θ
m
)条件下，阵列信号模型中的生成矩阵；阵列信号模型为：
[0017][0018]
记
[0019][0020]
s(f
i
)＝(s(f
i
，θ1) s2(f
i
，θ2)...s
m
(f
i
，θ
m
))
t
[0021]
x(f
i
)＝(x1(f
i
) x2(f
i
)...x
l
(f
i
))
t
[0022]
其中，s(f
i
)为待求的相干信源dft值；x(f
i
)为阵列实际接收的目标辐射噪声的线
谱dft值；g(f
i
，θ1，
…
，θ
m
)的第一行元素取值均为1，表示以第一个阵元为参考阵元，其它行元素表示其它阵元接收信号线谱dft值相对参考阵元接收信号线谱dft值的增量；g(f
i
，θ1，θ2，
…
，θ
m
)的第l行第k列元素取值为：
[0023][0024]
其中，d为阵元间距，满足d≤λ
i
/2，λ
i
为线谱频率f
i
对应的载波波长；
[0025]
将生成矩阵g(f
i
，θ1，θ2，
…
，θ
m
)和转换矩阵d(f
i
，θ1，
…
，θ
m
)进行非实时计算，并存储；
[0026]
每个频率的线谱需要预先存储的生成矩阵、转换矩阵的元素个数总数均为：
[0027][0028]
其中，l为阵元数目。
[0029]
进一步地，所述步骤2包括：
[0030]
首先，对水下目标辐射噪声进行观测接收和实时采集，做dft变换，实时计算l个阵元通道的基频线谱分量dft值的样本协方差矩阵：
[0031]
c＝x(f
i
)
h
·
x(f
i
)
[0032]
其中x(f
i
)＝(x1(f
i
)，x2(f
i
)，
…
，x
l
(f
i
))，为阵列实际接收的目标辐射噪声的线谱dft值；l为阵元数目；f
i
为线谱频率；
[0033]
然后，根据盖氏圆信源估计方法，由c的前(l
‑
1)行和前(l
‑
1)列构成新矩阵c1，对c1做特征值分解，计算c1的特征值和特征向量；将所述特征值从大到小排序，构造对角矩阵d1＝diag(λ
′1，λ
′2，
…
，λ
′
k
，
…
，λ
′
l
‑1)，其中k是待估计的相干辐射噪声信源个数，λ
′1、λ
′2…
λ
′
k
…
λ
′
l
‑11分别为从大到小排序的c1的特征值；根据所述特征向量构造酉矩阵u1＝[q
′1，q
′2，
…
，q
′
m
，
…
，q
′
l
‑1]，q
′1、q
′2…
q
′
m
…
q
′
l
‑1分别为λ
′1、λ
′2…
λ
′
k
…
λ
′
l
‑1对应的特征向量，满足再重构酉矩阵用生成的酉矩阵u对样本协方差矩阵c做酉变换：
[0034][0035]
其中，ρ1、ρ2、ρ3…
ρ
l
‑1、c
ll
为u
h
cu计算后矩阵中的元素；
[0036]
则按照下式估计相干信源数目：
[0037][0038]
其中，0≤d(n)≤1，是关于观测时间内共采集的噪声点数n的递减函数；gde(k)是用来帮助判断信源数目的一个过程量，当k＝k0时，gde(k0)为第一个比零小的值，则信源数的估计值取k＝k0‑
1。
[0039]
进一步地，所述步骤3包括：
[0040]
首先，依据x(f
i
)＝(x1(f
i
)，x2(f
i
)，
…
，x
l
(f
i
))做常规波束形成：
[0041]
d(θ
m
)＝a(θ
m
)
·
x
h
(k)，m＝1，2，
…
，m.
[0042]
其中，θ
m
为波达方向离散后第m个波达方向，a(θ
m
)为该波达方向的阵列导向矢量；
[0043]
基于d(θ
m
)计算功率，总共计算得到m个e(θ
m
)：
[0044]
e(θ
m
)＝d
*
(θ
m
)
·
d(θ
m
)，m＝1，2，
…
，m.
[0045]
然后，再根据步骤2估计得到的相干信源数目k，从m个e(θ
m
)中搜索峰值和次峰值，将峰值和次峰值所在方位记作当κ＝k时，以的每个方位值为中心，取峰值或次峰值
‑
3db宽度区间以内的离散方位值集合作为相干辐射噪声到达方向分布范围，记作当κ＜k时，以的每个方位值为中心，取峰值或次峰值
‑
6db宽度区间以内的离散方位值集合作为相干辐射噪声到达方向分布范围，记作
[0046]
进一步地，所述步骤4包括：
[0047]
首先，依据步骤2估计的相干信源数目k以及步骤3估计的目标辐射噪声相干信源到达方向确定相干信源到达方向的所有可能组合；记的元素个数分别为m1，m2，
…
，m
κ
；当κ＝k时，从κ个集合中各分配一个信源波达方向进行组合，组合数目为每个集合长度的乘积，即当κ＜k时，重新合并集合新集合元素个数为在中做k个相干信源的组合，组合数目为
[0048]
然后，在个波达方向组合中，对每一个相干信源波达方向的组合，记任意一个波达方向组合为查找预先存储的生成矩阵和转换矩阵计算k个相干信源的dft值：
[0049][0050]
以及模型偏差：
[0051][0052]
最后，对于个波达方向组合，总共有个模型偏差，从中搜索模型偏差最小值，将最小值对应的波达方向组合作为k个相干信源的真实波达方向，并将对应的dft值中每个元素的幅度平方值作为相干信源强度。
[0053]
与现有技术相比，本发明具有的有益效果：
[0054]
本发明通过非实时计算，预先存储阵列信号的生成矩阵和转换矩阵，大幅度降低了线谱分量信号空间的多维搜索求解复杂度；并通过盖氏圆方法和常规波束形成方法分别估计相干信源数目和波达方向分布范围，进一步降低线谱分量信号空间的多维搜索求解范
围。本发明充分利用了空间拟合类算法对相干信源波达方位估计精度高的优点，又大幅降低了实时计算复杂度，适用于多途条件下对水下目标辐射噪声的相干信源测向定位。
附图说明
[0055]
图1为本发明实施例一种基于阵列信号数值计算的水下目标辐射噪声相干信源测向方法的流程图。
具体实施方式
[0056]
下面结合附图和具体的实施例对本发明做进一步的解释说明：
[0057]
如图1所示，一种基于阵列信号数值计算的水下目标辐射噪声相干信源测向方法，具体步骤如下：
[0058]
s1：生成矩阵和转换矩阵的非实时计算存储：预先将水下目标辐射噪声不同频率线谱的到达方向做数值化离散，然后依据不同信源(多途)数目分别做组合排列，最后计算每个组合下的生成矩阵和转换矩阵；
[0059]
s2：目标辐射噪声相干信源数目的实时估计：实时计算各个阵元通道采集的目标辐射噪声中线谱分量的样本协方差矩阵，然后通过盖氏圆信源估计方法估计目标辐射噪声相干信源数目；
[0060]
s3：目标辐射噪声相干信源到达方向分布范围的实时估计：计算各阵元通道采集的目标辐射噪声的线谱分量的dft值，然后通过常规波束形成方法估计目标辐射噪声相干信源到达方向分布范围；
[0061]
s4：目标辐射噪声相干信源到达方向的实时精确计算：依据目标辐射噪声相干信源数目估计结果和到达方向分布范围，确定目标辐射噪声相干信源到达方向的所有可能组合以及阵列信号数学模型；然后选取预先存储的生成矩阵、转换矩阵与频域阵列信号向量相乘，计算模型系数与模型偏差；最后将模型偏差最小时对应的到达方向组合确定为目标辐射噪声的相干信源真实到达方向，并将此时的模型系数作为相干信源强度。
[0062]
具体地，所述s1包括：
[0063]
首先，将水下目标辐射噪声的不同频率线谱到达方向做数值化离散。当前，水下目标普遍采用螺旋桨推动的方式进行机动，在定速巡航状态下，目标的辐射噪声中一般包含了与螺旋桨转速成整倍数关系的多条线谱分量，线谱的强度往往与水下目标的排水量和航速成正比，因此可对目标辐射噪声的线谱分量进行观测。例如，对水下目标辐射噪声中的基频线谱分量进行观测。基频大小与水下目标的航速成正比，即由螺旋桨转速决定。当目标处于定速潜航状态时，该基频取值可以稳定分布在某个区间，记作f
basee
∈(f
a
，f
b
)，按照一定分辨率δf将其离散化，记作f∈{f1，f2，
…
，f
γ
}。进一步，对水下目标辐射噪声的波达方向进行离散。由于阵列接收具有指向性，能够观测的目标辐射噪声波达方向的范围有限，记作θ∈(θ
a
，θ
b
)，按照一定分辨率δθ将其离散化，记作θ∈{θ1，θ2，
…
，θ
m
}。
[0064]
然后，对于每一个线谱频率f，按照信源数目从1增加至(l
‑
3)，对其波达方向进行组合排列，组合个数分别为其中l为阵元数目，m为离散波达方向的数目。以频率f1的线谱为例，在1个信源的情况下，波达方向组合为(θ1)，(θ2)，
…
，(θ
m
)，共
个。在2个信源的情况下，波达方向组合为：
[0065]
(θ1，θ2)，(θ1，θ3)，
…
，(θ1，θ
m
)；(θ2，θ3)，(θ2，θ4)，
…
，(θ2，θ
m
)；
…
；(θ
m
‑1，θ
m
).共
[0066]
最后，对于每个频率的每个组合，计算转换矩阵：
[0067]
d(f
i
，θ1，
…
，θ
m
)＝[g
h
(f
i
，θ1，
…
，θ
m
)
·
g(f
i
，θ1，
…
，θ
m
)]
‑1·
g
h
(f
i
，θ1，
…
，θ
m
)
[0068]
其中，1＜m≤m，g(f
i
，θ1，
…
，θ
m
)表示在线谱频率f
i
以及波达方向组合(θ1，
…
，θ
m
)条件下，阵列信号模型中的生成矩阵。阵列信号模型为：
[0069][0070]
记
[0071][0072]
s(f
i
)＝(s(f
i
，θ1) s2(f
i
，θ2)
…
s
m
(f
i
，θ
m
))
t
[0073]
x(f
i
)＝(x1(f
i
) x2(f
i
)...x
l
(f
i
))
t
[0074]
其中，s(f
i
)为待求的相干信源dft值，x(f
i
)为阵列实际接收的目标辐射噪声的线谱dft值。g(f
i
，θ1，
…
，θ
m
)的第一行元素取值均为1，表示以第一个阵元为参考阵元，其它行元素表示其它阵元接收信号线谱dft值相对参考阵元接收信号线谱dft值的增量。对于构型固定和通道一致性经过校准的小规模接收阵列，g(f
i
，θ1，
…
，θ
m
)中的其它元素只与阵列构型、线谱频率f
i
以及波达方向组合(θ1，
…
，θ
m
)有关，可以预先计算。以均匀线阵为例，阵元数目为l，阵元间距为d，满足d≤λ
i
/2，λ
i
为线谱频率f
i
对应的载波波长，g(f
i
，θ1，θ2，
…
，θ
m
)的第l行第k列元素取值为：
[0075][0076]
因此，φ
l，k
、生成矩阵g(f
i
，θ1，θ2，
…
，θ
m
)和转换矩阵d(f
i
，θ1，
…
，θ
m
)都可非实时计算，并存储，计算量可忽略不计。每个频率的线谱需要预先存储的生成矩阵、转换矩阵的元素个数总数均分别为：
[0077][0078]
具体地，所述s2包括：
[0079]
首先，采用构型固定和通道一致性经过校准的小规模接收阵列，对水下目标辐射噪声进行观测接收和实时采集。假设阵元数目为l，采样频率f
s
，在观测时间内共采集噪声点数为n。第l个阵元通道采集的噪声样本记作x
l
＝(x
l
[1]，x
l
[2]，
…
，x
l
[n])，均为复数向量或iq数据。做dft变换，计算l个阵元通道的基频线谱分量dft值的样本协方差矩阵：
[0080]
c＝x(f
i
)
h
·
x(f
i
)
[0081]
其中x(f
i
)＝(x1(f
i
)，x2(f
i
)，
…
，x
l
(f
i
))，x
l
(f
i
)为第l个阵元实际接收的目标辐射噪声的线谱dft值。
[0082]
然后，参考盖氏圆信源估计方法[hsien
‑
tsai wu，jar
‑
fen yang and fwu
‑
kuen chen.source number estimation using transformed geschgorin radii[j].ieeetransactions on signal processing，1995.vol.43(6)：1325
‑
1333x]，由c的前(l
‑
1)行和前(l
‑
1)列构成新矩阵c1；对c1做特征值分解，计算c1的特征值和特征向量；将特征值从大到小排序，构造对角矩阵，记作d1＝diag(λ
′1，λ
′2，
…
，λ
′
k
，
…
，λ
′
l
‑1)，其中k是待估计的相干辐射噪声信源个数，λ
′1、λ
′2…
λ
′
k
…
λ
′
l
‑1分别为从大到小排序的c1的特征值；用所述特征向量构造酉矩阵，记作u1＝[q
′1，q
′2，
…
，q
′
m
，
…
，q
′
l
‑1]，q
′1、q
′2…
q
′
m
…
q
′
l
‑1分别为λ
′1、λ
′2…
λ
′
k
…
λ
′
l
‑1对应的特征向量，则满足再重构酉矩阵用生成的酉矩阵u对样本协方差矩阵c做酉变换：
[0083][0084]
其中，ρ1、ρ2、ρ3…
ρ
l
‑1、c
ll
为u
h
cu计算后矩阵中的元素，分别为ρ1、ρ2、ρ3…
ρ
l
‑1的共轭；
[0085]
进一步，计算
[0086][0087]
其中，0≤d(n)≤1，是关于n的递减函数，例如取值d(n)＝1/n。gde(k)是用来帮助判断信源数目的一个过程量，如果当k＝k0时，gde(k0)为第一个比零小的值，则信源数的估计值取k＝k0‑
1。显然k≤l
‑
3，即该盖氏圆方法适用于相干辐射噪声信源数目少于阵元数目3个以上的场合。例如阵元数为9，那么采用该盖氏圆方法可以估计信源数目不超过6个的相干辐射噪声。考虑到海洋中水声信道的多途现象具有随机性，所有多途噪声同时被阵列接收的概率较低，可近似为稀疏多径信道条件，因此该信源估计方法可以满足大部分水下目标测向定位需求。并且该盖氏圆信源估计方法在低信噪比和相干条件下具有良好的估计性能。
[0088]
具体地，所述s3包括：
[0089]
首先，依据所述s2计算的基频线谱分量的dft值向量x(f
i
)＝(x1(f
i
)，x2(f
i
)，
…
，x
l
(f
i
))做常规波束形成：
[0090]
d(θ
m
)＝a(θ
m
)
·
x
h
(k)，m＝1，2，
…
，m.
[0091]
其中，θ
m
为波达方向离散后第m个波达方向，a(θ
m
)为该波达方向的阵列导向矢量。进一步计算功率e(θ
m
)，总共计算得到m个e(θ
m
)。
[0092]
e(θ
m
)＝d
*
(θ
m
)
·
d(θ
m
)，m＝1，2，
…
，m.
[0093]
其中，d
*
(θ
m
)表示d(θ
m
)的共轭。
[0094]
然后，再根据s2估计得到的相干信源数目k，从m个e(θ
i
)中搜索峰值和其它次峰值，将峰值和次峰值所在方位记作当κ＝k时，以的每个方位值为中心，取峰值或次峰值
‑
3db宽度区间以内的离散方位值集合作为相干辐射噪声到达方向分布范围，记作当κ＜k时，以的每个方位值为中心，取峰值或次峰值
‑
6db宽度区间以内的离散方位值集合作为相干辐射噪声到达方向分布范围，记作
[0095]
具体地，所述s4包括：
[0096]
首先，依据s2估计的相干信源数目k以及s3估计的目标辐射噪声相干信源到达方向确定相干信源到达方向的所有可能组合。记的元素个数分别为m1，m2，
…
，m
κ
。当κ＝k时，从κ个集合中各分配一个信源波达方向进行组合，组合数目为每个集合长度的乘积，即当κ＜k时，重新合并集合新集合元素个数为在中做k个相干信源的组合，组合数目为
[0097]
然后，在个波达方向组合中，对每一个相干信源波达方向的组合，记任意一个波达方向组合为查找预先存储的生成矩阵和转换矩阵计算k个相干信源的dft值：
[0098][0099]
以及模型偏差：
[0100][0101]
最后，对于个波达方向组合，总共有个模型偏差，从中搜索模型偏差最小值，将最小值对应的波达方向组合作为k个相干信源的真实波达方向，并将对应的dft值中每个元素的幅度平方值作为相干信源强度。作为一种可实施方式，信源强度最大的波达方向可参考视为水下目标辐射噪声的最短路径到达方向。
[0102]
综上，本发明通过非实时计算，预先存储阵列信号的生成矩阵和转换矩阵，大幅度降低了线谱分量信号空间的多维搜索求解复杂度；并通过盖氏圆方法和常规波束形成方法分别估计相干信源数目和波达方向分布范围，进一步降低线谱分量信号空间的多维搜索求解范围。本发明充分利用了空间拟合类算法对相干信源波达方位估计精度高的优点，又大幅降低了实时计算复杂度，适用于多途条件下对水下目标辐射噪声的相干信源测向定位。
[0103]
以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应
视为本发明的保护范围。