CN105785322A - 一种波束形成方法 - Google Patents

Abstract

本发明涉及一种波束形成方法，包括：对接收阵所接收的数据做傅里叶变换分析，进行相位补偿、累加；对累加数据做共轭相乘，得到方位θ处的能量；计算出频率f处的各阵元接收数据协方差矩阵，记为第一协方差矩阵；对第一协方差矩阵做M阶次传感器互相关处理，得到虚拟后线阵；对虚拟后线阵做常规波束形成，做目标方位估计，得到该虚拟后线阵中各个阵元在方位θ处的能量；计算出频率f处的虚拟后线阵中各阵元接收数据协方差矩阵；判断当前的M阶次传感器互相关CBF在-3dB处的主瓣宽度是否达到用户所要求的分辨力，若尚未达到，令频率f处的虚拟后线阵中各阵元接收数据协方差矩阵作为第一协方差矩阵并将M的值增1后重新做M阶次传感器互相关处理，否则，结束操作。

Description

一种波束形成方法

技术领域

本发明涉及目标检测领域，特别涉及一种波束形成方法。

背景技术

在目标检测和方位估计中，波束形成(Beamforming，BF)对拖线阵具有重要作用。对于常规波束形成，角度分辨率无法突破“瑞利限”。随着实际情况的复杂多变化，对同时估计多个目标方位的需求越来越迫切。为了能够得到高分辨波束形成技术，相关学者提出了不同方法，大致可分为两类：

1、子空间类方法，该类方法是基于协方差矩阵特征分解理论，利用信号子空间和噪声子空间。正交特性使目标分辨能力显著提高，能够得到方位参数的渐近无偏估计，并使方位估计方差接近克拉美罗界(CRB)，突破了“瑞利限”制约，具有良好的参数估计性能。其中主要包括Pisarenko法、MUSIC法、Johnson法、Mini-Norm法、JESPRIT法、MODE法、加权子空间拟合(WSF)法等。这类方法己经成为阵列信号处理领域中研究的热点，工程应用前景良好。但该类方法对先验知识依赖性较强，对目标检测的最低门限信噪比要求较高。

2、参数模型方法，该类方法根据信号特点，利用参数模型(如AR模型、MA模型、ARMA模型等)来拟合信号过程，并进一步估计方位，较好地解决了加窗函数带来的不利影响。在信噪比不太低、数据长度不太短的情况下，这类方法具有比常规波束形成方法更好的估计性能，其中包括自相关法、协方差法、线性预测法等。但是这类方法估计性能受模型准确性及相应阶数的影响较大，对目标检测的最低门限信噪比要求较高。

随着降噪技术不断提高，拖线阵接收数据的信噪比会随之在不断降低，现有的高分辨波束形成技术已不能很好地满足拖线阵需要。为了能够得到同时具有高分辨和宽容性波束形成方法，可以通过深入研究以下3种方法得以实现：1.继续研究子空间类方法和参数模型方法，得到同时具备高分辨和宽容性的波束形成方法。2.通过合成孔径方法扩大线阵有效孔径，得到同时具备高分辨和宽容性的波束形成方法。3.通过虚拟阵元的方法，得到同时具备高分辨和宽容性的波束形成方法。

发明内容

本发明的目的在于克服现有技术中的波束形成技术无法同时满足高分辨率与宽容性的缺陷，从而提供一种同时具备高分辨和宽容性的波束形成方法。

为了实现上述目的，本发明提供了一种波束形成方法，包括：

步骤1)、对接收阵中的各个阵元所接收的数据做傅里叶变换分析，然后对傅里叶变换分析所得到的结果进行相位补偿、累加，各阵元在频率单元f处累加数据；最后对所述累加数据做共轭相乘，从而得到方位θ处的能量；由所述方位θ处的能量计算出频率f处的各阵元接收数据协方差矩阵，将该协方差矩阵记为第一协方差矩阵；

步骤2)、对第一协方差矩阵做M阶次传感器互相关处理，得到包括虚拟阵元与原始阵元的虚拟后线阵；其中，M的初始值为1；

步骤3)、对步骤2)所得到的虚拟后线阵做常规波束形成，得到M阶次传感器互相关CBF，对所述M阶次传感器互相关CBF做目标方位估计，得到该虚拟后线阵中各个阵元在方位θ处的能量；由虚拟后线阵中各个阵元在方位θ处的能量计算出频率f处的虚拟后线阵中各阵元接收数据协方差矩阵；

步骤4)、判断当前的M阶次传感器互相关CBF在-3dB处的主瓣宽度是否达到用户所要求的分辨力，若尚未达到，令步骤3)所得到的频率f处的虚拟后线阵中各阵元接收数据协方差矩阵作为第一协方差矩阵并将M的值增1后重新执行步骤2)，否则，结束操作。

上述技术方案中，所述步骤1)包括：

步骤1-1)、接收阵中的各个阵元所接收的数据x_i(t)表示为：

x_i(t)＝s(t+△τ_i)+n_i(t),0≤i≤N-1(1)；

式中，s(t)为目标辐射信号，t为时刻，△τ_i＝d·i·cos(θ₀)/c为目标辐射信号到第i阵元相比到第0号参考阵元的时延差，i为阵元号，d为阵间距，c为参考声速，n_i(t)为第i阵元在t时刻接收的背景噪声；

步骤1-2)、对所接收的数据x_i(t)做傅里叶变换分析得：

$X_i=\left(f\right)=\mathrm{FFT}\left(x_i\left(t\right)\right)=\mathrm{FFT}(s(t+{\mathrm{\Δ\τ}}_i)+n_i\left(t\right))=S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_i}+N_i\left(f\right),(0\≤i\≤N-1,f_l\≤f\≤f_h)---\left(2\right)$

式中，f_l为滤波器频带下限，f_h为滤波器频带上限；

步骤1-3)、然后对步骤1-2)所得到的X_i(f)进行相位补偿、累加，得各阵元在频率单元f处累加数据：

$Y\left(f\right)=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{X}_i\left(f\right)\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)},---\left(3\right)$

式中，τ_i(θ)＝d·i·cos(θ)/c为进行相位补偿所设第i阵元相比第0号参考阵元的时延差；

步骤1-4)、对步骤1-3)所得到的Y(f)进行共轭相乘，得方位θ处的能量为P(θ)，然后将P(θ)改写为A(θ)R(f)A^H(θ)的形式。

P(θ)的计算公式如下：

$\begin{array}{c}P\left(\mathrm{\θ}\right)=Y\left(f\right)\·Y^{*}\left(f\right)\\ =\{\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\π}{\mathrm{f\Δ\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}\}\·\{\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\π\Δf}{\mathrm{\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}{\}}^{*}\\ =A\left(\mathrm{\θ}\right)R_{\mathrm{SS}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{SN}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{NS}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{NN}}\left(f\right)A^H\left(\mathrm{\θ}\right)\\ =A\left(\mathrm{\θ}\right)R\left(f\right)A^H\left(\mathrm{\θ}\right),\end{array}---\left(4\right);$

式中，(·)^*为共轭，(·)^H为共轭转置，A(θ)、R_SS(f)、R_SN(f)、R_NS(f)、R_NN(f)如下式所示：

$A\left(\mathrm{\θ}\right)=[e^{-j2\mathrm{\πf}{\mathrm{\τ}}_0\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_1\left(\mathrm{\θ}\right)},...,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}\left(\mathrm{\θ}\right)}],---\left(5A\right)$

$R_{\mathrm{SS}}\left(f\right)={[S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}}]}^T\·{[S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^{*},---\left(5B\right)$

$R_{\mathrm{SN}}\left(f\right)={[S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N1}}]}^T\·{[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^{*},---\left(5C\right)$

$R_{\mathrm{NS}}\left(f\right)={[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^T\·{[S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^{*},---\left(5D\right)$

$R_{\mathrm{NN}}\left(f\right)={[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^T\·{[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^{*},---\left(5E\right)$

式中，(·)^T为转置，一般情况下，信号和背景噪声是不相关的，数据经分段处理后累加得：N_i(f)S^*(f)≈0,(0≤i≤N-1)、，则R_SN(f)≈R_NS(f)≈0，R(f)如式(6)所示：

$\begin{array}{c}R\left(f\right)=R_{\mathrm{SS}}\left(f\right)+R_{\mathrm{SN}}\left(f\right)+R_{\mathrm{NS}}\left(f\right)+R_{\mathrm{NN}}\left(f\right)\≈R_{\mathrm{SS}}\left(f\right)+R_{\mathrm{NN}}\left(f\right)\\ =\left[\begin{array}{c}S\left(f\right)S^{*}\left(f\right)+N_0\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1}+N_0\left(f\right)N_1^{*}\left(f\right),...,S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}}+N_0\left(f\right)N_{N-1}^{*}\left(f\right)\\ S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0}+N_1\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)+N_1\left(f\right)N_1^{*}\left(f\right),...,S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}}+N_1\left(f\right)N_{N-1}^{*}\left(f\right)\\ .\\ .\\ .\\ S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0}+N_{N-1}\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1}+N_{N-1}\left(f\right)N_1^{*}\left(f\right),...,S\left(f\right)S^{*}\left(f\right)+N_{N-1}\left(f\right)N_{N-1}^{*}\left(f\right)\end{array}\right]\end{array},---\left(6\right).$

上述技术方案中，在步骤2)中，对第一协方差矩阵R_f进行传感器互相关处理所得结果的表达式如下：

$Z_P\left(f\right)=\left\{\begin{array}{c}\frac{1}{N-\left|P\right|}\underset{n=-P}{\overset{N-1}{\mathrm{\Σ}}}{R}_f(n,n+P),1-N\≤P\≤0\\ \frac{1}{N-P}\underset{n=1}{\overset{N-P}{\mathrm{\Σ}}}{R}_f(n,n+P),1\≤P\≤N-1\end{array}\right.---\left(7\right)$

式中，R_f为频率f处的各阵元接收数据协方差矩阵，n、P为所取R_f中元素位置；Z_-N+1(f),Z_-N+2(f),…,Z₀(f),…,Z_N-2(f),Z_N-1(f)表示为：

$Z_{-N+1}\left(f\right)=S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0}+N_{N-1}\left(f\right)N_0\left(f\right)=S^{\′}\left(f\right)e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{-N+1}}+{N^{\′}}_{-N+1}\left(f\right),---\left(8A\right)$

$\begin{array}{c}{Z}_{-N+2}\left(f\right)=\frac{S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-2}-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0}+N_{N-2}\left(f\right)N_0\left(f\right)+S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1}+N_{N-1}\left(f\right)N_1\left(f\right)}{2}\\ =\frac{(S\left(f\right)S\left(f\right)+S\left(f\right)S\left(f\right))e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{-N+2}}+(N_{N-2}\left(f\right)N_0\left(f\right)+N_{N-1}\left(f\right)N_1\left(f\right))}{2}=S^{\′}\left(f\right)e^{j2\mathrm{\πf}\mathrm{\Δ}{{\mathrm{\τ}}^{\′}}_{-N+2}}+{N^{\′}}_{-N+2}\left(f\right),\end{array}---\left(8B\right)$

$Z_0\left(f\right)=\frac{S\left(f\right)S\left(f\right)+N_0\left(f\right)N_0\left(f\right)+...+S\left(f\right)S\left(f\right)+N_{N-1}\left(f\right)N_{N-1}\left(f\right)}{N}=S^{\′}\left(f\right)+{N^{\′}}_0\left(f\right),---\left(8C\right)$

$\begin{array}{c}{Z}_{N-2}\left(f\right)=\frac{S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-2}}+N_{N-2}\left(f\right)N_0\left(f\right)+S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_{N-1}\left(f\right)N_1\left(f\right)}{2}\\ =\frac{(S\left(f\right)S\left(f\right)+S\left(f\right)S\left(f\right))e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{N-2}}+(N_{N-2}\left(f\right)N_0\left(f\right)+N_{N-1}\left(f\right)N_1\left(f\right))}{2}=S^{\′}\left(f\right)e^{j2\mathrm{\πf}\mathrm{\Δ}{{\mathrm{\τ}}^{\′}}_{N-2}}+{N^{\′}}_{N-2}\left(f\right),\end{array}---\left(8D\right)$

$Z_{N-1}\left(f\right)=S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_{N-1}\left(f\right)N_0\left(f\right)=S^{\′}\left(f\right)e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{N-1}}+{N^{\′}}_{N-1}\left(f\right),---\left(8E\right).$

上述技术方案中，在步骤3)中，虚拟后线阵中各个阵元在方位θ处的能量为：

P'(θ)＝A'(θ)R'_fA'^H(θ)(9)

其中， $A^{\′}\left(\mathrm{\θ}\right)=[e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{-N+1}\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{-N+2}\left(\mathrm{\θ}\right)},...,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_0\left(\mathrm{\θ}\right)},...,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-2}\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}\left(\mathrm{\θ}\right)}],$

上述技术方案中，在步骤3)中，M阶次传感器互相关CBF在-3dB处的主瓣宽度为：

其中，λ为波长。

本发明的优点在于：

1、本发明的方法通过对协方差矩阵进行M(M≥2)阶次CrossSensor处理，可虚拟出近(2^M-1)(N-1)个虚拟阵元(N为原始阵元数)。虚拟阵元的增加可扩大线阵有效孔径，降低波束形成主瓣宽度，提高方位分辨率。

2、在CrossSensor处理过程中，本发明的方法对协方差矩阵同一斜对角线上不同元素进行了叠加运算，可进一步削弱噪声对波束形成的影响，提高了波束形成鲁棒性。

3、数值仿真和海上试验结果表明，本发明的方法既能有效降低波束形成主瓣宽度，提高方位分辨率，又可削弱方位历程图干扰背景。相比已有的逆波束形成，本发明的方法具有较好的方位分辨率；相比已有的基于AR模型的高分辨逆波束形成，该方法对最低门限信噪比的要求较低，方位估计均方误差较小。

附图说明

图1是本发明方法所适用的接收阵的示意图；

图2是采用本发明的方法对接收阵做多次CrossSensor处理后，接收阵中包括原始阵元与虚拟阵元在内的阵元变化情况的示意图；

图3为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得正确分辨双目标概率曲线；

图4为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得空间谱(SNR＝-10dB)的示意图；

图5为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得正确分辨双目标概率曲线的示意图；

图6为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得空间谱(SNR＝-10dB)的示意图；

图7为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得方位估计均方误差曲线的示意图；

图8为在一个海试实验中所用水平拖曳线列阵以及目标运动的示意图；

图9为在一个数据长度为700s，采样率为f_s＝20kHz，滤波器频带为f＝75～80Hz的实验中，CBF所得方位历程图；

图10为在一个数据长度为700s，采样率为f_s＝20kHz，滤波器频带为f＝90～95Hz的实验中，CBF所得方位历程图；

图11(a)是6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得1s处的空间谱的示意图；

图11(b)是CBF方法所得方位历程图；

图11(c)是FIM-IBF方法所得方位历程图；

图11(d)是2阶次CrossSensorCBF法所得方位历程图；

图11(e)是3阶次CrossSensorCBF法所得方位历程图；

图11(f)是4阶次CrossSensorCBF法所得方位历程图；

图11(g)是AR-IBF方法所得方位历程图。

具体实施方式

现结合附图对本发明作进一步的描述。

在对本发明的方法做详细说明前，首先对本发明的方法所适用的接收阵加以描述。图1为一接收阵示意图，该接收阵为阵元数为N的等间距水平拖线阵，将这N个阵元依次记为第0、1、2、……、N-1号阵元；其中，第0号阵元为参考阵元。θ₀为目标相对水平拖线阵辐射信号方向。

基于图1所示接收阵，本发明的波束形成方法包括：

步骤1)、对接收阵中的各个阵元所接收的数据做傅里叶变换分析，然后对傅里叶变换分析所得到的结果进行相位补偿、累加，各阵元在频率单元f处累加数据；最后对所述累加数据做共轭相乘，从而得到方位θ处的能量。

下面对步骤1)中所完成的操作做详细描述。

步骤1-1)、接收阵中的各个阵元所接收的数据x_i(t)可表示为：

x_i(t)＝s(t+△τ_i)+n_i(t),0≤i≤N-1(1)

式中，s(t)为目标辐射信号，t为时刻，△τ_i＝d·i·cos(θ₀)/c为目标辐射信号到第i阵元相比到第0号参考阵元的时延差，i为阵元号，d为阵间距，c为参考声速，n_i(t)为第i阵元在t时刻接收的背景噪声。

步骤1-2)、对所接收的数据x_i(t)做傅里叶变换分析可得：

$X_i=\left(f\right)=\mathrm{FFT}\left(x_i\left(t\right)\right)=\mathrm{FFT}(s(t+{\mathrm{\Δ\τ}}_i)+n_i\left(t\right))=S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_i}+N_i\left(f\right),(0\≤i\≤N-1,f_l\≤f\≤f_h)---\left(2\right)$

式中，f_l为滤波器频带下限，f_h为滤波器频带上限。

步骤1-3)、然后对X_i(f)进行相位补偿、累加，可得各阵元在频率单元f处累加数据：

$Y\left(f\right)=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{X}_i\left(f\right)\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)},---\left(3\right)$

式中，τ_i(θ)＝d·i·cos(θ)/c为进行相位补偿所设第i阵元相比第0号参考阵元的时延差。

步骤1-4)、对步骤1-3)所得到的Y(f)进行共轭相乘，可得方位θ处的能量为P(θ)，然后将P(θ)改写为A(θ)R(f)A^H(θ)的形式。

P(θ)的计算公式如下：

$\begin{array}{c}P\left(\mathrm{\θ}\right)=Y\left(f\right)\·Y^{*}\left(f\right)\\ =\{\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\π}{\mathrm{f\Δ\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}\}\·\{\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}{\}}^{*}\\ =A\left(\mathrm{\θ}\right)R_{\mathrm{SS}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{SN}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{NS}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{NN}}\left(f\right)A^H\left(\mathrm{\θ}\right)\\ =A\left(\mathrm{\θ}\right)R\left(f\right)A^H\left(\mathrm{\θ}\right),\end{array}---\left(4\right);$

式中，(·)^*为共轭，(·)^H为共轭转置，A(θ)、R_SS(f)、R_SN(f)、R_NS(f)、R_NN(f)如下式所示：

$A\left(\mathrm{\θ}\right)=[e^{-j2\mathrm{\πf}{\mathrm{\τ}}_0\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_1\left(\mathrm{\θ}\right)},...,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}\left(\mathrm{\θ}\right)}],---\left(5A\right)$

$R_{\mathrm{SS}}\left(f\right)={[S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}}]}^T\·{[S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^{*},---\left(5B\right)$

$R_{\mathrm{SN}}\left(f\right)={[S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N1}}]}^T\·{[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^{*},---\left(5C\right)$

$R_{\mathrm{NS}}\left(f\right)={[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^T\·{[S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},...,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^{*},---\left(5D\right)$

$R_{\mathrm{NN}}\left(f\right)={[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^T\·{[N_0\left(f\right),N_1\left(f\right),...,N_{N-1}\left(f\right)]}^{*},---\left(5E\right)$

式中，(·)^T为转置，一般情况下，信号和背景噪声是不相关的，数据经分段处理后累加可得：N_i(f)S^*(f)≈0,(0≤i≤N-1)、，则R_SN(f)≈R_NS(f)≈0，R(f)如式(6)所示：

$\begin{array}{c}R\left(f\right)=R_{\mathrm{SS}}\left(f\right)+R_{\mathrm{SN}}\left(f\right)+R_{\mathrm{NS}}\left(f\right)+R_{\mathrm{NN}}\left(f\right)\≈R_{\mathrm{SS}}\left(f\right)+R_{\mathrm{NN}}\left(f\right)\\ =\left[\begin{array}{c}S\left(f\right)S^{*}\left(f\right)+N_0\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1}+N_0\left(f\right)N_1^{*}\left(f\right),...,S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}}+N_0\left(f\right)N_{N-1}^{*}\left(f\right)\\ S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0}+N_1\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)+N_1\left(f\right)N_1^{*}\left(f\right),...,S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}}+N_1\left(f\right)N_{N-1}^{*}\left(f\right)\\ .\\ .\\ .\\ S\left(f\right)S^{*}{e}^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0}+N_{N-1}\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1}{N}_{N-1}\left(f\right)N_1^{*}\left(f\right),...,S\left(f\right)S^{*}\left(f\right)+N_{N-1}\left(f\right)N_{N-1}^{*}\left(f\right)\end{array}\right]\end{array},---\left(6\right)$

步骤2)、令R_f＝R(f)，对R_f进行传感器互相关(CrossSensor)处理，生成N-1个虚拟阵元，这些虚拟阵元与原始线阵中的阵元一起组成一个新的虚拟后线阵。

对R_f进行CrossSensor处理所得结果的表达式如下：

$Z_P\left(f\right)=\left\{\begin{array}{c}\frac{1}{N-\left|P\right|}\underset{n=-P}{\overset{N-1}{\mathrm{\Σ}}}{R}_f(n,n+P),1-N\≤P\≤0\\ \frac{1}{N-P}\underset{n=1}{\overset{N-P}{\mathrm{\Σ}}}{R}_f(n,n+P),1\≤P\≤N-1\end{array}\right.---\left(7\right)$

式中，R_f为频率f处的各阵元接收数据协方差矩阵，n、P为所取R_f中元素位置。Z_-N+1(f),Z_-N+2(f),…,Z₀(f),…,Z_N-2(f),Z_N-1(f)可表示为：

$Z_{-N+1}\left(f\right)=S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0}+N_{N-1}\left(f\right)N_0\left(f\right)=S^{\′}\left(f\right)e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{-N+1}}+{N^{\′}}_{-N+1}\left(f\right),---\left(8A\right)$

$\begin{array}{c}{Z}_{-N+2}\left(f\right)=\frac{S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-2}-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_0}+N_{N-2}\left(f\right)N_0\left(f\right)+S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1}+N_{N-1}\left(f\right)N_1\left(f\right)}{2}\\ =\frac{(S\left(f\right)S\left(f\right)+S\left(f\right)S\left(f\right))e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{-N+2}}+(N_{N-2}\left(f\right)N_0\left(f\right)+N_{N-1}\left(f\right)N_1\left(f\right))}{2}=S^{\′}\left(f\right)e^{j2\mathrm{\πf}\mathrm{\Δ}{{\mathrm{\τ}}^{\′}}_{-N+2}}+{N^{\′}}_{-N+2}\left(f\right),\end{array}---\left(8B\right)$

$Z_0\left(f\right)=\frac{S\left(f\right)S\left(f\right)+N_0\left(f\right)N_0\left(f\right)+...+S\left(f\right)S\left(f\right)+N_{N-1}\left(f\right)N_{N-1}\left(f\right)}{N}=S^{\′}\left(f\right)+{N^{\′}}_0\left(f\right),---\left(8C\right)$

$\begin{array}{c}{Z}_{N-2}\left(f\right)=\frac{S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf\Δ}{\mathrm{\τ}}_{N-2}}+N_{N-2}\left(f\right)N_0\left(f\right)+S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_{N-1}\left(f\right)N_1\left(f\right)}{2}\\ =\frac{(S\left(f\right)S\left(f\right)+S\left(f\right)S\left(f\right))e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{N-2}}+(N_{N-2}\left(f\right)N_0\left(f\right)+N_{N-1}\left(f\right)N_1\left(f\right))}{2}=S^{\′}\left(f\right)e^{j2\mathrm{\πf}\mathrm{\Δ}{{\mathrm{\τ}}^{\′}}_{N-2}}+{N^{\′}}_{-N+2}\left(f\right),\end{array}---\left(8D\right)$

$Z_{N-1}\left(f\right)=S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_{N-1}\left(f\right)N_0\left(f\right)=S^{\′}\left(f\right)e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{N-1}}+{N^{\′}}_{N-1}\left(f\right),---\left(8E\right)$

依据阵元间协方差矩阵同一斜对角线上不同元素具有相同相位差的特点，即,0≤p,q≤N-1，可令△τ'_-N+1＝△τ_N-1-△τ₀、△τ'_-N+2＝△τ_N-2-△τ₀＝△τ_N-1-△τ₁、……、△τ'_N-2＝△τ₀-△τ_N-2＝△τ₁-△τ_N-1、△τ'_N-1＝△τ₀-△τ_N-1作为虚拟后线阵中各阵元(包括原始阵元与虚拟阵元)接收信号时延差，并令S'(f)＝S(f)S(f)、S'(f)＝(S(f)S(f)+S(f)S(f))/²、……、S'(f)＝(S(f)S(f)+S(f)S(f))/²、S'(f)＝S(f)S(f)作为虚拟后线阵中各阵元接收的信号，令N'_-N+1(f)＝N_N-1(f)N₀(f)、N'_-N+2(f)＝(N_N-2(f)N₀(f)+N_N-1(f)N₁(f))/²、……、N'_N-2(f)_＝(N_N-2(f)N₀(f)+N_N-1(f)N₁(f))/²、N'_N-1(f)＝N_N-1(f)N₀(f)作为虚拟后线阵中各阵元接收的背景噪声。

由式(8A)—式(8E)所示可知，Z_P(f),1-N≤P≤N-1已具有线列阵接收信号特性。

步骤3)、对步骤2)所得到的Z_P(f),1-N≤P≤N-1做常规波束形成(ConventionalBeam-Forming，CBF)，得到M阶次传感器互相关CBF，然后对所述M阶次传感器互相关CBF做目标方位估计，可得方位θ处的能量为：

P'(θ)＝A'(θ)R'_fA'^H(θ)(9)

其中， $A^{\′}\left(\mathrm{\θ}\right)=[e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{-N+1}\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{-N+2}\left(\mathrm{\θ}\right)},...,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_0\left(\mathrm{\θ}\right)},...,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-2}\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}\left(\mathrm{\θ}\right)}],$

步骤4)、令R'_f为虚拟线阵后各阵元接收数据的协方差矩阵R_f，重复M(M≥2)次步骤2)、步骤3)的处理过程，可依次完成M(M≥2)阶次CrossSensorCBF。

在步骤4)中，可完成M(M≥2)阶次CrossSensorCBF。下面对这一处理的必要性进行陈述。

由式(7)可知，Z_P(f)所含阵元数相比原始线阵阵元数多了近一倍。对Z_P(f),1-N≤P≤N-1进行波束形成，-3dB处的主瓣宽度可表示为：

式中，λ为波长，L'＝2L＝2(N-1)d为有效孔径，L＝(N-1)d为原始线阵有效孔径。

经过一次CrossSensor处理可增加近一倍的阵元数，阵元数的增加可扩大线阵有效孔径，由式(10)可知，在λ不变情况下，线阵有效孔径的增加可有效降低主瓣宽度，提高波束形成方位分辨率。相比式(4)，式(9)所得波束形成主瓣宽度缩减了近1/3。

经过M(M≥2)阶次重复处理式(7)-式(9)的过程，可虚拟出近(2^M-1)(N-1)个阵元，(2^M-1)(N-1)个阵元数的增加可使有效孔径扩大近(2^M-1)(N-1)倍(如图2所示)，所得波束形成空间谱主瓣宽度会缩减近波束形成方位分辨率得到有效提高，同时可进一步削弱噪声对波束形成结果的影响，进一步提高了波束形成鲁棒性。

令信号和背景噪声不相关，各阵元接收背景噪声是独立高斯噪声。在理论上，M(M≥2)阶次CrossSensorCBF在-3dB处的主瓣宽度，M(M≥2)阶次CrossSensorCBF的空间增益分别为：

G_SM≤10log(2^M·N-2^M+1)dB(12)

一旦由式(11)所得到的主瓣宽度达到用户所要求的分辨力，则步骤4)中所提到的M阶次CrossSensor处理即可终止。由式(11)和式(12)可知，相比CBF，M(M≥2)阶次CrossSensorCBF可使波束形成主瓣宽度减小近(2^M-1)(N-1)/(2^MN-2^M+1)(θ＝90°)，有效削弱噪声对波束形成的影响。

由式(7)及式(10)—(12)可知，M(M≥2)阶次CrossSensorCBF的优越性为：

1.相比CBF，会缩减近噪声抑制能力得到提高，具有更高的方位分辨力，在目标谱估计失真方面优于CBF。

2.M(M≥2)阶次CrossSensorCBF本质上仍然是线性空间谱估计，具有线性空间谱估计算法的优点。

本发明的方法(简称M(M≥2)阶次CrossSensorCBF)与现有技术中的方法(简称CBF、FIM-IBF、AR-IBF)相比具有明显的优点。下面结合实例，对本发明方法与现有技术中相关方法的效果进行比较。

在一个比较例中，采用水平拖线阵阵元数为N＝16，阵间距为d＝c/2f_c，f_c＝120Hz为中心频率，c＝1500m/s为参考声速。采样率为f_s＝5kHz，每次采样时间为T＝1s，并假设背景噪声为各向同性的宽带高斯噪声，频带为f＝40～200Hz，得到的仿真结果如下：

1)、非相干情况下分辨多目标

目标1辐射信号频率为f₁＝69Hz、幅度为A₁＝1v，相对水平拖线阵方位角为θ₁＝87°(t＝1～500s)；目标2辐射信号，目标2辐射信号频率为f₂＝107Hz、幅度为A₂＝0.8v，相对水平拖线阵方位角为θ₂＝90°(t＝1～500s)。

图3为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得正确分辨双目标概率曲线，该图中CBF、FIM-IBF、2阶次CrossSensorCBF法的分辨概率为0，三条曲线重合在一起。

图4为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得空间谱(SNR＝-10dB)。

2)、相干情况下分辨多目标

目标1辐射信号频率为f₁＝100Hz、幅度为A₁＝1v，相对水平拖线阵方位角为θ₁＝86.5°(t＝1～500s)；目标2辐射信号频率为f₂＝100Hz、幅度为A₂＝0.8v，相对水平拖线阵方位角为θ₂＝90°(t＝1～500s)。

图5为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得正确分辨双目标概率曲线，该图中CBF、FIM-IBF、AR-IBF法分辨概率为0，三条曲线重合在一起。

图6为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得空间谱(SNR＝-10dB)。

3)、方位估计的均方误差

目标辐射信号频率为f₁＝100Hz、幅度为A₁＝1v，相对水平拖线阵方位角为θ₁＝90°(t＝1～500s)。

图7为不同信噪比情况下，6种方法(即CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法)所得方位估计均方误差曲线。

由图3和图5可知，

1.因原始线阵阵元数有限，多目标所在方位相对接收线阵的角度间隔小于“瑞利限”时，受“瑞利限”限制，无论目标辐射信号相关与否，阵元接收数据信噪比如何，CBF和FIM-IBF已不能正确分辨双目标；

2.通过高阶次CrossSensor处理，可虚拟出近(2^M-1)(N-1)个阵元，线阵有效孔径可增加近(2^M-1)(N-1)倍，因虚拟出阵元数和线阵有效孔径得到了有效增加，相比原始线阵已突破原始线阵的“瑞利限”限制。在仿真条件下，a)非相干情况，在同一正确分辨双目标概率下，3阶次CrossSensorCBF法相比AR-IBF方法对最低门限信噪比的要求下降了近5dB，4阶次CrossSensorCBF法相比AR-IBF方法对最低门限信噪比的要求下降了近6dB；b)相干情况，无论阵元接收数据信噪比如何，AR-IBF正确分辨双目标概率小于10％；而2阶次CrossSensorCBF在SNR≥-10dB时，正确分辨双目标概率大于50％，3阶次CrossSensorCBF在SNR≥-16dB时，正确分辨双目标概率大于50％，4阶次CrossSensorCBF在SNR≥-17dB时，正确分辨双目标概率大于50％。图3和图5仿真结果表明，经高阶次CrossSensor处理的波束形成具有较窄主瓣宽度和较弱噪声影响，经3阶次以上CrossSensor处理的波束形成在分辨双目标时，对最低门限信噪比的要求低于AR-IBF法，具有较好的鲁棒性。

由图7可知，因AR-IBF法只用到线阵接收数据协方差矩阵逆的部分信息，而较多信息被舍弃，致使其具有较大的方位估计均方误差；而CBF、FIM-IBF、高阶次CrossSensorCBF利用了线阵接收数据协方差矩阵的全部信息，方位估计均方误差较小，相比AR-IBF高分辨波束形成方法，高阶次CrossSensorCBF具有较好的鲁棒性。

在一个比较例中，采用海试数据进行分析验证，海试数据为2012年4月在南海进行目标检测试验所得。实验所用水平拖曳线列阵以及目标运动示意图如图8所示，实验采用32元拖曳线列阵接收信号，相邻阵元间隔为8m，拖曳线列阵的端向方位设为0°。

本次处理实验数据长度为700s，所用采样率为f_s＝20kHz。

1)、令滤波器频带为f＝75～80Hz

图9为CBF所得方位历程图；

2)、令滤波器频带为f＝90～95Hz

图10为CBF所得方位历程图。

由图9可知，在θ＝71°附近存在1个目标，在θ＝75°附近存在1个微弱目标；由图10可知，进一步验证了在θ＝75°附近存在1个微弱目标。

3)、令滤波器频带为f＝60～120Hz，分60个子带单元按不同方法进行窄带波束形成，然后再将窄带波束形成结果进行累加合成得到如下结果。图11(a)是不同方法所得1s处的空间谱。图11(b)—图11(g)分别为CBF方法、FIM-IBF方法、2阶次CrossSensorCBF法、3阶次CrossSensorCBF法、4阶次CrossSensorCBF法、AR-IBF方法所得方位历程图。

由图11可知，在由6种方法所得空间谱中，由目标运动轨迹的连续性可知，在t＝1s～100s内，在θ＝50°～60°、θ＝70°～75°内分别存在多目标。经2阶次以上CrossSensor处理的波束形成所得空间谱主瓣较窄，可以实现对θ＝50°～60°、θ＝70°～75°内双目标的有效检测，而CBF、FIM-IBF及AR-IBF法所得空间谱已不能实现对θ＝50°～60°、θ＝70°～75°内双目标的有效检测检测和分辨，且噪声对AR-IBF法波束形成影响较大，在其所得空间谱中背景能量与目标方位处能量差别较小，不利于目标检测和方位估计。同样对比图11(b)～图11(g)可知，经2阶次以上CrossSensor处理的波束形成所得波束主瓣较窄、方位历程图中目标方位显示效果较好，可以实现对θ＝50°～60°、θ＝70°～75°内双目标的辨别；而CBF、FIM-IBF及AR-IBF法所得方位历程图中目标方位显示效果较差，已不能实现对θ＝50°～60°、θ＝70°～75°内双目标的有效辨别，且受噪声影响，在AR-IBF法所得波束形成方位历程图中背景最亮，不利于实现对目标的有效辨别。海上试验数据处理结果表明，经2阶次以上CrossSensor处理的波束形成在分辨多目标时，对最低门限信噪比的要求低于AR-IBF法，具有较好的鲁棒性。

数值仿真结果和海试数据处理结果表明，在不同的信噪比下情况下，不同阶次CrossSensorCBF法所得空间谱具有不同主瓣宽度、方位分辨率，并受不同程度的噪声影响；并验证了应用本文所提出的基于高阶次CrossSensor处理的鲁棒性高分辨波束形成方法可以扩大线阵有效孔径，降低波束形成主瓣宽度，提高方位分辨率，有效降低方位历程图中的背景干扰对波束形成的影响，扩大目标方位与其他方位处的能量差别，改善了波束形成方位历程图显示效果。相比已有的逆波束形成(FIM-IBF)，该方法具有较好的方位分辨率；相比已有的基于AR模型的高分辨逆波束形成(AR-IBF)，该方法对最低门限信噪比的要求较低，方位估计均方误差较小，且可有效分辨相干信源，具有较好鲁棒性。另外，在本发明方法基础上，可联合子空间分类法、参数模型法、合成孔径法等技术手段进一步优化波束形成所得空间谱。

实际使用时可根据式(11)来选择所需阶数M(M≥2)，以便满足方位估计分辨力需求；同时也要折中考虑高阶次CrossSensor处理所增加的运算量。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种波束形成方法，包括：

步骤1)、对接收阵中的各个阵元所接收的数据做傅里叶变换分析，然后对傅里叶变换分析所得到的结果进行相位补偿、累加，各阵元在频率单元f处累加数据；最后对所述累加数据做共轭相乘，从而得到方位θ处的能量；由所述方位θ处的能量计算出频率f处的各阵元接收数据协方差矩阵，将该协方差矩阵记为第一协方差矩阵；

步骤2)、对第一协方差矩阵做M阶次传感器互相关处理，得到包括虚拟阵元与原始阵元的虚拟后线阵；其中，M的初始值为1；

步骤3)、对步骤2)所得到的虚拟后线阵做常规波束形成，得到M阶次传感器互相关CBF，对所述M阶次传感器互相关CBF做目标方位估计，得到该虚拟后线阵中各个阵元在方位θ处的能量；由虚拟后线阵中各个阵元在方位θ处的能量计算出频率f处的虚拟后线阵中各阵元接收数据协方差矩阵；

步骤4)、判断当前的M阶次传感器互相关CBF在-3dB处的主瓣宽度是否达到用户所要求的分辨力，若尚未达到，令步骤3)所得到的频率f处的虚拟后线阵中各阵元接收数据协方差矩阵作为第一协方差矩阵并将M的值增1后重新执行步骤2)，否则，结束操作。

2.根据权利要求1所述的波束形成方法，其特征在于，所述步骤1)包括：

步骤1-1)、接收阵中的各个阵元所接收的数据x_i(t)表示为：

x_i(t)＝s(t+Δτ_i)+n_i(t),0≤i≤N-1(1)；

式中，s(t)为目标辐射信号，t为时刻，Δτ_i＝d·i·cos(θ₀)/c为目标辐射信号到第i阵元相比到第0号参考阵元的时延差，i为阵元号，d为阵间距，c为参考声速，n_i(t)为第i阵元在t时刻接收的背景噪声；

步骤1-2)、对所接收的数据x_i(t)做傅里叶变换分析得：

$X_i\left(f\right)=\mathrm{FFT}\left(x_i\left(t\right)\right)=\mathrm{FFT}(s(t+{\mathrm{\Δ\τ}}_i)+n_i\left(t\right))=S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_i}+N_i\left(f\right),(0\≤i\≤N-1,f_l\≤f\≤f_h)---\left(2\right)$

式中，f_l为滤波器频带下限，f_h为滤波器频带上限；

步骤1-3)、然后对步骤1-2)所得到的X_i(f)进行相位补偿、累加，得各阵元在频率单元f处累加数据：

$Y\left(f\right)=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{X}_i\left(f\right)\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)},---\left(3\right)$

式中，τ_i(θ)＝d·i·cos(θ)/c为进行相位补偿所设第i阵元相比第0号参考阵元的时延差；

步骤1-4)、对步骤1-3)所得到的Y(f)进行共轭相乘，得方位θ处的能量为P(θ)，然后将P(θ)改写为A(θ)R(f)A^H(θ)的形式；

P(θ)的计算公式如下：

$\begin{array}{c}P\left(\mathrm{\θ}\right)=Y\left(f\right)\·Y^{*}\left(f\right)\\ =\{\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}\}\·{\{\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_i}\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}+\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{N}_i\left(f\right)e^{-j2\mathrm{\πf}{\mathrm{\τ}}_i\left(\mathrm{\θ}\right)}\}}^{*}\\ =A\left(\mathrm{\θ}\right)R_{\mathrm{SS}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{SN}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{NS}}\left(f\right)A^H\left(\mathrm{\θ}\right)+A\left(\mathrm{\θ}\right)R_{\mathrm{NN}}\left(f\right)A^H\left(\mathrm{\θ}\right)\\ =A\left(\mathrm{\θ}\right)R\left(f\right)A^H\left(\mathrm{\θ}\right),\end{array}---\left(4\right);$

式中，(·)^*为共轭，(·)^H为共轭转置，A(θ)、R_SS(f)、R_SN(f)、R_NS(f)、R_NN(f)如下式所示：

$A\left(\mathrm{\θ}\right)=[e^{-j2\mathrm{\πf}{\mathrm{\τ}}_0\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_1\left(\mathrm{\θ}\right)},\·\·\·,e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}\left(\mathrm{\θ}\right)}],---\left(5A\right)$

$R_{\mathrm{SS}}\left(f\right)={[S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1},\·\·\·,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^T\·{[S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0},S\left(f\right)e^{j2f{\mathrm{\Δ\τ}}_1},\·\·\·,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^{*},---\left(5B\right)$

$R_{\mathrm{SN}}\left(f\right)={[S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1},\·\·\·,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N1}}]}^T\·{[N_0\left(f\right),N_1\left(f\right),\·\·\·,N_{N-1}\left(f\right)]}^{*},---\left(5C\right)$

$R_{\mathrm{NS}}\left(f\right)={[N_0\left(f\right),N_1\left(f\right),\·\·\·,N_{N-1}\left(f\right)]}^T\·{[S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0},S\left(f\right)e^{j2\mathrm{\πf\Δ}{\mathrm{\τ}}_1},\·\·\·,S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}]}^{*},---\left(5D\right)$

$R_{\mathrm{NN}}\left(f\right)={[N_0\left(f\right),N_1\left(f\right),\·\·\·,N_{N-1}\left(f\right)]}^T\·{[N_0\left(f\right),N_1\left(f\right),\·\·\·,N_{N-1}\left(f\right)]}^{*},---\left(5E\right)$

式中，(·)^T为转置，一般情况下，信号和背景噪声是不相关的，数据经分段处理后累加得：N_i(f)S^*(f)≈0,(0≤i≤N-1)、则R_SN(f)≈R_NS(f)≈0，R(f)如式(6)所示：

$\begin{array}{c}R\left(f\right)=R_{\mathrm{SS}}\left(f\right)+R_{\mathrm{SN}}\left(f\right)+R_{\mathrm{NS}}\left(f\right)+R_{\mathrm{NN}}\left(f\right)\≈R_{\mathrm{SS}}\left(f\right)+R_{\mathrm{NN}}\left(f\right)\\ \left[\begin{array}{c}S\left(f\right)S^{*}\left(f\right)+N_0\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1}+N_0\left(f\right)N_1^{*}\left(f\right),\·\·\·,S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_0\left(f\right)N_{N-1}^{*}\left(f\right)\\ S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0}+N_1\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)+N_1\left(f\right)N_1^{*}\left(f\right),\·\·\·,S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_1\left(f\right)N_{N-1}^{*}\left(f\right)\\ \·\\ \·\\ \·\\ \mathrm{Sf}{S}^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2f{\mathrm{\Δ\τ}}_0}+N_{N-1}\left(f\right)N_0^{*}\left(f\right),S\left(f\right)S^{*}\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1}+N_{N-1}\left(f\right)N_1^{*}\left(f\right),\·\·\·,S\left(f\right)S^{*}\left(f\right)+N_{N-1}\left(f\right)N_{N-1}^{*}\left(f\right)\end{array}\right]\end{array},---\left(6\right).$

3.根据权利要求2所述的波束形成方法，其特征在于，在步骤2)中，对第一协方差矩阵R_f进行传感器互相关处理所得结果的表达式如下：

$Z_P\left(f\right)=\left\{\begin{array}{c}\frac{1}{N-\left|P\right|}\underset{n=-P}{\overset{N-1}{\mathrm{\Σ}}}{R}_f(n,n+P),1-N\≤P\≤0\\ \frac{1}{N-P}\underset{n=1}{\overset{N-P}{\mathrm{\Σ}}}{R}_f(n,n+P),1\≤P\≤N-1\end{array}\right.---\left(7\right)$

式中，R_f为频率f处的各阵元接收数据协方差矩阵，n、P为所取R_f中元素位置；Z_-N+1(f),Z_-N+2(f),…,Z₀(f),…,Z_N-2(f),Z_N-1(f)表示为：

$Z_{-N+1}\left(f\right)=S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0}+N_{N-1}\left(f\right)N_0\left(f\right)=S^{\′}\left(f\right)e^{j2\mathrm{\πf}\mathrm{\Δ}{{\mathrm{\τ}}_{-N+1}^{\′}}_{\text{\′}}}+N_{-N+1}^{\′}\left(f\right),---\left(8A\right)$

$\begin{array}{c}{Z}_{-N+2}\left(f\right)=\frac{S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-2}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0+N_{N-2}\left(f\right)N_0\left(f\right)+S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1}+N_{N-1}\left(f\right)N_1\left(f\right)}}{2}\\ =\frac{(S\left(f\right)S\left(f\right)+S\left(f\right)S\left(f\right))e^{j2\mathrm{\πf}\mathrm{\Δ}{\mathrm{\τ}}_{-N+2}^{\′}}+(N_{N-2}\left(f\right)N_0\left(f\right)+N_{N-1}\left(f\right)N_1\left(f\right))}{2}=S^{\′}\left(f\right)e^{j2\mathrm{\πf}{{\mathrm{\Δ\τ}}^{\′}}_{-N+2}}+{N^{\′}}_{-N+2}\left(f\right),\end{array}---\left(8B\right)$

$Z_0\left(f\right)=\frac{S\left(f\right)S\left(f\right)+N_0\left(f\right)N_0\left(f\right)+\·\·\·+S\left(f\right)S\left(f\right)+N_{N-1}\left(f\right)N_{N-1}\left(f\right)}{N}=S^{\′}\left(f\right)+{N^{\′}}_0\left(f\right),---\left(8C\right)$

$\begin{array}{c}{Z}_{N-2}\left(f\right)=\frac{S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-2}+N_{N-2}\left(f\right)N_0\left(f\right)+S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_1-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_{N-1}\left(f\right)N_1\left(f\right)}}{2}\\ =\frac{(S\left(f\right)S\left(f\right)+S\left(f\right)S\left(f\right))e^{j2\mathrm{\πf}\mathrm{\Δ}{\mathrm{\τ}}_{N-2}^{\′}}+(N_{N-2}\left(f\right)N_0\left(f\right)+N_{N-1}\left(f\right)N_1\left(f\right))}{2}=S^{\′}\left(f\right)e^{j2\mathrm{\πf}{{\mathrm{\Δ\τ}}^{\′}}_{N-2}}+{N^{\′}}_{N-2}\left(f\right)\end{array},---\left(8D\right)$

$Z_{N-1}\left(f\right)=S\left(f\right)S\left(f\right)e^{j2\mathrm{\πf}{\mathrm{\Δ\τ}}_0-j2\mathrm{\πf}{\mathrm{\Δ\τ}}_{N-1}}+N_{N-1}\left(f\right)N_0\left(f\right)=S^{\′}\left(f\right)e^{j2\mathrm{\πf\Δ}{{\mathrm{\τ}}^{\′}}_{N-1}}+{N^{\′}}_{N-1}\left(f\right),---\left(8E\right).$

4.根据权利要求3所述的波束形成方法，其特征在于，在步骤3)中，虚拟后线阵中各个阵元在方位θ处的能量为：

P′(θ)＝A′(θ)R′_fA^′H(θ)(9)

其中， $A^{\′}\left(\mathrm{\θ}\right)=[e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{-N+1}\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{\_N+2}\left(\mathrm{\θ}\right)},\·\·\·e^{-j2\mathrm{\πf}{\mathrm{\τ}}_0\left(\mathrm{\θ}\right)},e^{-j2\mathrm{\πf}{\mathrm{\τ}}_{N-2}\left(\mathrm{\θ}\right)},e^{j2\mathrm{\πf}{\mathrm{\τ}}_{N-1}\left(\mathrm{\θ}\right)}],$

5.根据权利要求4所述的波束形成方法，其特征在于，在步骤3)中，M阶次传感器互相关CBF在-3dB处的主瓣宽度为：

其中，λ为波长。