一种运动干扰环境下的贝叶斯稳健波束形成方法
技术领域
本发明属于阵列信号处理领域，涉及一种运动干扰环境下的贝叶斯稳健波束形成方法。
背景技术
提高波束形成器干扰抑制能力的一种有效措施是采用自适应波束形成技术，在干扰方位自适应产生凹陷，从而提高信干噪比(Signal-to-interference-plus-noiseratio，SINR)。然而，由于自适应方向图通常具有陡峭的凹口，故运动干扰源易移出方向图“零陷”，导致传统稳健自适应波束形成算法性能于非平稳干扰环境中急剧下降。在针对运动干扰源的波束形成系统中，由于干扰信号方位随时间快速变化，当干扰方位先验未知时，无法提供融合各离散时间点上的波束输出结果以最终得到运动干扰高效抑制的有效策略。在这种情况下，如何有效地综合多个离散时间点处的数据，从而实现针对非平稳干扰环境的波束形成，以获得比已有方法更强的非理想信号环境适应能力，正是本发明所要解决的一个关键问题。
运动干扰抑制问题广泛存在于阵列波束形成领域，该问题可视为多模型多观测问题在阵列信号处理领域的延伸。受模型差异的影响，用于解决多模型多观测系统中运动干扰抑制问题的方法与已有单模型系统中的静态干扰抑制方法不尽相同。J.R.Guerci(J.R.Guerci,Theory and application of covariance matrix tapers for robustadaptive beamforming,IEEE Trans.Signal Process.47(4)(1999)977-985.)提出一种基于采样协方差矩阵锥化思路的运动干扰抑制算法。该协方差矩阵锥化(Covariance MatrixTaper,CMT)算法通过展宽干扰方位处的凹口实现增强波束形成器稳健性的目的。然而此算法仅关注运动干扰源的抑制问题，对期望信号导向矢量误差不具有稳健性。采用接收数据的各阶导作为波束形成算法的限制条件，A.B.Gershman等人(A.B.Gershman,U.Nickel,J.F.Adaptive beamforming algorithms with robustness against jammermotion,IEEE Trans.Signal Process.45(7)(1997)1878-1885.)(A.B.Gershman,E.Németh,J.F.Experimental performance of adaptive beamforming in a sonarenvironment with a towed array and moving interfering sources,IEEETrans.Signal Process.48(1)(2000)246-250.)将传统加载采样协方差矩阵求逆算法和特征向量投影算法加以改进，推广到运动干扰环境中。然而，上述两种算法面临对角加载因子难以精确计算、低SNR下噪声及信号子空间混叠等技术实现难点。S.A.Vorobyov等人(S.A.Vorobyov,A.B.Gershman,Z.Luo,N.Ma,Adaptive beamforming with jointrobustness against mismatched signal steering vector and interferencenonstationarity,IEEE Signal Process.Lett.11(2)(2004)108-111.)将Worst Case算法推广至非平稳干扰环境中，然而此算法实现过程中面临阵列接收数据协方差矩阵及期望信号导向矢量不确定集先验未知等困难。L.Zhang等人(L.Zhang,B.Li,L.Huang,T.Kirubarajan,H.C.So,Robust minimum dispersion distortionless responsebeamforming against fast-moving interferences,Signal Process.140(2017)190-197.)提出一种基于最小离差准则的稳健波束形成算法，其通过约束干扰动态方位集内的平均功率为0，达到从波束形成器输出中消除干扰能量的目的。该算法的缺陷为全局收敛性难以保证，易发散到局部最优解。
鉴于贝叶斯方法的多级概率结构能够消除不同信号分量之间的耦合，因此可以较好地保留各信号分量的局部特征，从而有助于获得较高的输出SINR，本发明提出一种结合狄利克雷过程综合利用信号整个观测时间内方位信息的非参数化贝叶斯波束形成方法。所提方法能充分利用干扰方位的时变特征消除传统方法波束图中可能出现的凹口偏移，同时较好地避免干扰信号阵列流形实时变化带来的负面影响。
发明内容
要解决的技术问题
为了避免现有技术的不足之处，本发明提出一种运动干扰环境下的贝叶斯稳健波束形成方法。
技术方案
由于非参数化贝叶斯方法是一种综合利用不同模型所得观测数据提取感兴趣目标信息的贝叶斯联合重构方法，因而可用于解决运动干扰抑制这一包含多个观测模型的波束形成问题。本发明采用狄利克雷过程依据实际信号环境自适应调整权矢量和跟踪干扰方位，并基于局域变分推断准则，着力解决分层概率模型下由非共轭分布对导致的隐变量后验分布难以估算的问题。
假设期望信号和干扰信号均为远场窄带信号，其入射到一个由N个阵元组成的均匀线阵上，且阵元间距为入射信号波长的一半。当利用该线阵的接收信号进行自适应波束形成时，本发明所采用的技术方案包括以下步骤：
步骤1：建立阵列接收信号模型；
步骤2：建立贝叶斯分层概率模型，具体内容如下：
以狄利克雷过程描述干扰+噪声分量的方位时变特性，即
nk～G
式中，～表示服从于，为离散随机变量，且在其支撑集上的概率和为1，δ(·)为狄拉克函数，πf～Beta(1,γ)∝γ(1-πf)γ-1，∝表示正比于，γ表示伸缩因子，且服从如下伽马分布：
第f组干扰+噪声分量共享一精度矩阵，即：
式中G0表示基准分布，Λf表示第f组高斯分布对应的精度矩阵，服从如下Wishart分布
p(Λf)＝W(Λf|W，v)∝|Λf|v-Netr{-W-1Λf}
式中，表示均值矩阵，v表示自由度；
引入分配向量zk，其为从一多值分布中抽取的随机变量(多值分布的参数集为{ωf}f＝1,…,K)，即：
式中zk是仅有一个元素为1的K×1维向量；
定义1[zk＝f]表示zk中的第f个元素为1，即将nk分配到第f类，则阵列观测数据的似然函数服从如下复高斯分布：
式中是对应的精度矩阵；
以复高斯分布表示期望信号幅度的统计特性：
式中s＝[s1,…,sK]T为由K个快拍的期望信号波形幅度组成的向量，为期望信号精度；
期望信号精度服从如下伽马分布
式中，Γ(·)表示伽马函数，a＞0表示形状参数，b＞0表示伸缩参数；
期望信号导向矢量服从如下Watson分布：
式中，为正则化参数，λ为聚焦参数；cp(λ)中的超几何函数定义为如下形式：
超参数μ,λ服从如下Watson-Gamma联合分布
用集合表示概率模型中全部未知变量，并将其统计特性与观测数据Y的分布结合起来得到如下联合概率密度：
步骤3：利用变分推断方法，得到各隐变量的后验分布参数的更新公式，具体内容如下：
(1)初始化概率模型参数，即令v＝μ＝m0＝vp，其中vp为预设的期望信号导向矢量，Λf＝IN×N(f＝1，…,K)，πf(f＝1,…,K)＝0.5，v＝N，a＝a0＝b＝b0＝c＝d＝10-6，s＝1K×1，W＝10+6·IN×N，其中1K×1为K×1维的全1向量，IN×N为N×N维的单位矩阵；
(2)第f组干扰+噪声分量Λf的后验概率q(Λf)服从Wishart分布，其自由度和均值矩阵的更新公式分别为：
式中，定义〈q(zk×f)〉＝φk,f，〈·〉表示求数学期望操作；φk,f的计算公式为：
式中，
(3)πf的后验概率服从贝塔分布，即q(πf)＝Beta(πf|cf,df)，其中超参数cf和df的更新公式分别为：
式中，
(4)伸缩因子γ的后验分布为伽马分布，即其中形状参数和伸缩参数的更新公式分别为：
式中〈ln(1-πf)〉＝Ψ(df)-Ψ(cf+df)，Ψ(·)为digamma函数；
(5)第k次快拍中期望信号波形幅度的后验概率q(sk)服从复高斯分布，即其均值和精度的迭代更新公式为
(6)期望信号导向矢量v的后验概率服从Watson分布，即其中超参数和分别是矩阵的主特征向量和主特征值，且
(7)μ的后验分布为Watson分布，即其中是正定矩阵(βmmH+vvH)的最大特征值，是其对应的特征向量；
(8)λ的后验分布为Gamma分布，即q(λ)＝Gam(λ|a1，b1)，其中形状参数a1和伸缩参数b1的更新公式分别为：
上式中，从上一迭代步骤中计算得到；
(9)期望信号精度的后验概率服从伽马分布，即其中形状参数a2和伸缩参数b2的更新公式分别为：
a2＝K+a
(10)迭代(2)-(9)过程，直至满足收敛条件；
步骤4：将第f组干扰+噪声精度矩阵及期望信号导向矢量的收敛解和vopt代入公式计算得到第f组阵列接收数据对应的最优权矢量其余各组阵列接收数据对应的最优权矢量采用同样计算方法得到。
本发明进一步的技术方案为：所述的步骤1阵列接收信号模型具体如下：
以一阵元数为N的线阵为例，其第k个采样时刻的接收数据可表示为：
式中分别表示k时刻的阵列输出，期望信号波形幅度及干扰+噪声分量，v表示期望信号导向矢量，Y＝[y1,…,yK]表示观测时间内的阵列接收数据集，K是采样快拍数；在线阵信号模型下，导向矢量一般可写为如下形式：
式中λ0为载波波长，[d1，…,dN]为各阵元物理坐标，θ0为期望信号的空间方位。
本发明进一步的技术方案为：所述的步骤3中的收敛条件为其中为当前步计算出的期望信号精度值，为上一步计算出的期望信号精度值。
一种计算机系统，其特征在于包括：一个或多个处理器，计算机可读存储介质，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述的方法。
一种计算机可读存储介质，其特征在于存储有计算机可执行指令，所述指令在被执行时用于实现上述的方法。
一种计算机程序，其特征在于包括计算机可执行指令，所述指令在被执行时用于实现上述的方法。
有益效果
贝叶斯波束形成领域的现有理论成果还无法很好地支持对运动干扰的有效抑制。现有波束形成算法中关于观测数据中所包含的信号空域分布信息的静态先验假设与阵列实际模型之间还存在较大的距离，它们并不能直接用于解决方位时变信号这一包含多个观测模型的波束形成问题。然而，在一些典型的雷达和声呐等应用环境中，尽管干扰入射方向未知且随时间快速变化，但其在一段时间间隔内的方位特征却是稳定的，对这一特征的利用有助于简化波束形成方法的实现过程。
因此，针对上述具体技术背景，本发明首先提出综合利用不同模型所得观测数据提取感兴趣目标信息的贝叶斯联合重构方法，即非参数化贝叶斯方法-狄利克雷过程，以增强对阵列流形失配、小样本等信号环境的适应能力，有效利用不同时刻观测数据之间共同的空域结构这一内在联系，进而改善基于单一观测模型的数据重构方法的精度。随后，本发明采用交替迭代优化的方法，实现滤波器参数的联合估计。变分推断算法可用于实现对多组参数的同步优化，其有效性在处理包括非共轭先验分布在内的各种复杂分层概率模型中得到了验证。变分推断方法的这一优势在参数估计方面的表现是由凸近似所带来的局部极值解减少，以及在适宜样本数和信噪比条件下，其在保证全局收敛的基础上估计精度更好地逼近相应的理论下界。与传统概率采样方法相比，本发明所采用的变分推断方法具有更高的计算效率。
本发明提出综合利用不同模型所得观测数据提取感兴趣目标信息的贝叶斯联合重构方法，即非参数化贝叶斯方法-狄利克雷过程，以增强对运动干扰环境的适应能力，有效改善基于单一观测模型的数据重构方法的精度。在参数估计过程中，本发明采用交替迭代优化的方法，实现滤波器参数的联合估计，其优势在参数估计方面的表现是由凸近似所带来的局部极值解减少，以及在适宜样本数和信噪比条件下，其在保证全局收敛的基础上估计精度更好地逼近相应的理论下界。本发明所采用的技术途径可根据观测数据的结构特性对其进行自动聚类(类别属性由狄利克雷过程中的混合因子标识)，因而能够有效利用不同时刻观测数据之间共同的空域结构这一内在联系，改善贝叶斯权系数的重构精度。
由于本发明所采用的技术途径可根据观测数据的结构特性对其进行自动聚类(类别属性由狄利克雷过程中的混合因子标识)，因而所得技术成果对多任务学习、非均匀采样条件下的频率估计、基于多观测平台的目标检测与跟踪等领域中的相关研究也具有较大的参考价值。
本发明的基本原理和实施方案经过了计算机数值仿真的验证。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图1为由DOA估计误差引致的期望信号导向矢量失配情形下，实施实例中分别利用本发明所提方法与其它五种稳健波束形成方法获得的阵列输出SINR随输入SNR变化的曲线，以及它们与最优输出SINR曲线的对比图。
图2为由相干局域散射引致的期望信号导向矢量失配情形下，实施实例中分别利用本发明所提方法与其它五种稳健波束形成方法获得的阵列输出SINR随输入SNR变化的曲线，以及它们与最优输出SINR曲线的对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明解决其技术问题所采用的技术方案可以分为以下5个步骤：
步骤一：使用阵元位置坐标为d1,…,dN的N元均匀线阵作为接收阵列，接收远场窄带期望与干扰信号。线阵上各个传感器阵元将接收到的物理信号转换为电信号，并通过放大电路和数据采集器得到离散时域信号。
步骤二：假设阵列采样快拍数为K，则各采样时刻的阵列接收数据可表示为：
式中，yk，v和nk均为N×1维向量，分别表示采样时刻k的阵列输出信号，目标分量及干扰+噪声分量。sk为采样时刻k下期望信号波形的幅度。目标、干扰及噪声分量假设相互独立。v为期望信号导向矢量，其表达式为其中λ0为入射信号波长，{d1,d2,…,dN}为阵元位置坐标，θ0为期望信号的DOA，Y＝[y1,…,yK]表示观测时间内的阵列接收数据集
步骤三：根据步骤二中建立的阵列接收信号模型，建立分层概率框架，包括以下子步骤：
子步骤一：以狄利克雷过程描述干扰+噪声分量的方位时变特性，即
nk～G
式中，～表示服从于，为离散随机变量，且在其支撑集上的概率和为1，δ(·)为狄拉克函数，πf～Beta(1,γ)∝γ(1-πf)γ-1，∝表示正比于，γ表示伸缩因子，且服从如下伽马分布：
子步骤二：第f组干扰+噪声分量共享一精度矩阵，即：
式中G0表示基准分布，Λf表示第f组高斯分布对应的精度矩阵，服从如下Wishart分布
p(Λf)＝W(Λf|W，v)∝|Λf|v-Netr{-W-1Λf}
式中，表示均值矩阵，v表示自由度；
子步骤三：引入分配向量zk，其为从一多值分布中抽取的随机变量(多值分布的参数集为{ωf}f＝1，…，K)，即：
式中zk是仅有一个元素为1的K×1维向量；
子步骤四：定义1[zk＝f]表示zk中的第f个元素为1，即将nk分配到第f类，则阵列观测数据的似然函数服从如下复高斯分布：
式中是对应的精度矩阵；
子步骤五：以复高斯分布表示期望信号幅度的统计特性：
式中s＝[s1,…,sK]T为由K个快拍的期望信号波形幅度组成的向量，为期望信号精度；
子步骤六：期望信号精度服从如下伽马分布
式中，Γ(·)表示伽马函数，a＞0表示形状参数，b＞0表示伸缩参数；
子步骤七：期望信号导向矢量服从如下Watson分布：
式中，为正则化参数，λ为聚焦参数。cp(λ)中的超几何函数定义为如下形式：
超参数μ,λ服从如下Watson-Gamma联合分布
子步骤八：以集合表示概率模型中全部未知变量，并将其统计特性与观测数据Y的分布结合起来得到如下联合概率密度：
步骤四：利用变分推断方法计算Θ中各变量的后验分布参数，包括以下子步骤：
子步骤一：初始化概率模型参数，即令v＝μ＝m0＝vp，其中vp为预设的期望信号导向矢量，Λf＝IN×N(f＝1,…,K)，πf(f＝1,…,K)＝0.5，v＝N，a＝a0＝b＝b0＝c＝d＝10-6，s＝1K×1，W＝10+6·IN×N，其中1K×1为K×1维的全1向量，IN×N为N×N维的单位矩阵；
子步骤二：第f组干扰+噪声分量Λf的后验概率q(Λf)服从Wishart分布，其自由度和均值矩阵的更新公式分别为：
式中，定义〈q(zk＝f)〉＝φk,f，〈·〉表示求数学期望操作。φk,f的计算公式为：
式中，
子步骤三：πf的后验概率服从贝塔分布，即q(πf)＝Beta(πf|cf,df)，其中超参数cf和df的更新公式分别为：
式中，
子步骤四：伸缩因子γ的后验分布为伽马分布，即其中形状参数和伸缩参数的更新公式分别为：
式中〈ln(1-πf)〉＝Ψ(df)-Ψ(cf+df)，Ψ(·)为digamma函数；
子步骤五：第k次快拍中期望信号波形幅度的后验概率q(sk)服从复高斯分布，即其均值和精度的迭代更新公式为
子步骤六：期望信号导向矢量v的后验概率服从Watson分布，即其中超参数和分别是矩阵的主特征向量和主特征值，且
子步骤七：μ的后验分布为Watson分布，即其中是正定矩阵(βmmH+vvH)的最大特征值，是其对应的特征向量；
子步骤八：λ的后验分布为Gamma分布，即q(λ)＝Gam(λ|a1,b1)，其中形状参数a1和伸缩参数b1的更新公式分别为：
上式中，从上一迭代步骤中计算得到；
子步骤九：期望信号精度的后验概率服从伽马分布，即其中形状参数a2和伸缩参数b2的更新公式分别为：
a2＝K+a
子步骤十：迭代子步骤二至九过程，直至满足收敛条件，即其中为当前步计算出的期望信号精度值，为上一步计算出的期望信号精度值；
步骤五：将第f组干扰+噪声精度矩阵及期望信号导向矢量的收敛解和vopt代入公式计算得到第f组阵列接收数据对应的最优权矢量其余各组阵列接收数据对应的最优权矢量采用同样计算方法得到。
利用计算机进行数值仿真，检验本发明所提方法的估计性能。
仿真采用阵元间距为半波长的10元均匀线阵，各阵元接收噪声服从均值为0，方差为1的复高斯分布，且相互独立。阵列入射信号总数为3，其中期望信号数为1，等功率运动干扰信号数为2，且各入射信号波形均服从复高斯分布。干噪比(Interference-TO-noiseratio,INR)设为30dB，干扰信号入射DOA随采样时刻的变化规律分别为θ1(k)＝30°+0.1°k和θ2(k)＝60°-0.1°k。期望信号入射DOA的预设值为5°。蒙特卡洛实验次数设为100，每次实验中，各阵元位置误差服从[-0.05λ0,0.05λ0]范围内的均匀分布，其中λ0为入射信号波长。第k个采样时刻阵列输出SINR的计算公式为其中wk和Ri+n(k)分别表示采样时刻k的权矢量和干扰+噪声协方差(Interference-plus-noisecovariance,INC)矩阵，采样时刻k的最优阵列输出SINR可将最优权矢量代入上式求解得到。
1)存在DOA估计误差时，比较本发明所提方法与现有五种稳健波束形成方法的输出SINR结果
考虑由DOA估计误差引致的期望信号导向矢量失配情形，即每次蒙特卡洛实验中，期望信号DOA的预设值与真实值间的误差服从[-3°,3°]范围内的均匀分布。阵列采样快拍数为50。利用所提方法与现有五种方法，即CMT方法(J.R.Guerci,Theory and applicationof covariance matrix tapers for robust adaptive beamforming,IEEE Trans.SignalProcess.47(4)(1999)977-985.),ROBUST LSMI方法(A.B.Gershman,U.Nickel,J.F.Adaptive beamforming algorithms with robustness against jammer motion,IEEE Trans.Signal Process.45(7)(1997)1878-1885.),ROBUST EP方法(A.B.Gershman,E.Németh,J.F.Experimental performance of adaptive beamforming in asonar environment with a towed array and moving interfering sources,IEEETrans.Signal Process.48(1)(2000)246-250.),WORST CASE方法(S.A.Vorobyov,A.B.Gershman,Z.Luo,N.Ma,Adaptive beamforming with joint robustness againstmismatched signal steering vector and interference nonstationarity,IEEESignal Process.Lett.11(2)(2004)108-111.),MDDR方法(L.Zhang,B.Li,L.Huang,T.Kirubarajan,H.C.So,Robust minimum dispersion distortionless responsebeamforming against fast-moving interferences,Signal Process.140(2017)190-197.)进行波束形成，绘制输出SINR曲线。图1为六种方法的波束形成结果与最优SINR对比图。从图1所示结果可知，所提方法由于具有更高的INC矩阵和期望信号导向矢量估计精度，因而其干扰抑制性能优于现有方法。
2)存在相干局域散射时，比较本发明所提方法与现有五种稳健波束形成方法的输出SINR结果
真实期望信号导向矢量可表示为其中为预设的期望信号导向矢量，表示相干散射路径，ηi,i＝1,2,3,4为各路径上对应的相位。各次蒙特卡洛实验中，参数{θ0i}，{ηi}独立同分布，其中{θ0i}服从均值为3°，方差为1°的高斯分布，(ηi}服从[0,2π]内的均匀分布。设置快拍数为50，选择CMT、ROBUST LSMI、ROBUST EP、WORST CASE和MDDR五种方法作为本发明所提方法的比较对象，分别计算每种方法在不同输入SNR条件下获得的波束输出SINR。图2所示为上述六种方法的输出SINR及最优SINR随输入SNR的变化曲线。从图2所示结果可知，所提方法具有优于其他现有方法的干扰抑制性能，反映了在建模过程中就对信号特征加以利用的思路是具有优势的，该思路可有效利用运动干扰的空域非平稳性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明公开的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。
