一种平稳干扰环境下的贝叶斯稳健波束形成方法
技术领域
本发明属于阵列信号处理领域，涉及一种平稳干扰环境下的贝叶斯稳健波束形成方法。
背景技术
自适应波束形成是阵列信号处理的一项重要研究内容，其通过空域滤波增强期望入射方向的信号，抑制其它入射方向的干扰信号。最小方差无失真响应(Minimum variancedistortionless response,MVDR)算法是具有代表性的一类自适应波束形成方法，该算法在约束期望信号的响应恒为1的前提下，最小化波束输出方差，以达到抗干扰目的。MVDR波束形成器具有远高于传统数据独立波束形成器的输出信干噪比(Signal-to-interference-plus-noise ratio,SINR)。然而实际阵列观测数据中常混杂有期望信号，导致理想干扰+噪声协方差(Interference-plus-noise covariance,INC)矩阵不易获得，只能以采样协方差矩阵近似，此时，若预设的期望信号导向矢量存在偏差，MVDR波束形成器会将期望信号误认为干扰而进行抑制，引起信号自消效应。
稳健自适应波束形成(Robust adaptive beamforming，RAB)算法作为MVDR波束形成算法的一种改进形式，提供了平稳干扰环境下对抗导向矢量与协方差矩阵失配的一种有效技术手段。RAB算法的典型代表为最差(worst case)性能优化算法(S.A.Vorobyov,A.B.Gershman，Z.-Q.Luo，Robust adaptive beamforming using worst-caseperformance optimization:a solution to the signal mismatch problem,IEEETrans.Signal Process.51(2)(2003)313-324.)，其通过约束期望信号导向矢量不确定集内的阵列响应恒为1，展宽波束主瓣，消除信号自消效应。此算法的缺点为最优用户参数选取困难，计算复杂度高。A.Khabbazibasmenj等人提出一种消除worst case算法先验信息依赖性的RAB算法(A.Khabbazibasmenj,S.A.Vorobyov,A.Hassanien,Robust adaptivebeamforming based on steering vector estimation with as little as possibleprior information,IEEE Trans.Signal Process.60(6)(2012)2974-2987.)，其通过自适应修正期望信号导向矢量，使波束主瓣指向真实目标，减小了阵列自由度损耗。该算法同样基于内点法求解所构造的凸优化问题，因而计算复杂度较高，限制了其实际应用范围。Z.Zhang等人提出一种旨在解决协方差矩阵失配问题的基于空间功率谱采样(Spatialpower spectrum sampling,SPSS)的INC矩阵重构方法(Z.Zhang,W.Liu,W.Leng,A.Wang,H.Shi,Interference-plus-Noise covariance matrix reconstruction via spatialpower spectrum sampling for robust adaptive beamforming,IEEE SignalProcess.Lett.23(1)(2016)121-125.)，其利用空域功率谱采样技术估计INC矩阵，以剔除阵列快拍数据中所含的期望信号分量，然而此方法在阵元位置存在误差时性能损失严重。
此外，采用贝叶斯观点描述期望信号导向矢量误差的统计特性也是提高波束形成器稳健性的一条有效技术途径。K.L.Bell等人首先提出了贝叶斯波束形成器的概念(K.L.Bell,Y.Ephraim,H.L.Van Trees,A Bayesian approach to robust adaptivebeamforming,IEEE Trans.Signal Process.48(2)(2000)386-398.)，其可视作一组主瓣方位指向期望信号波达方向(Direction of arrival,DOA)不确定集内各离散点的MVDR波束形成器的加权和，权值由各DOA采样点的后验分布决定。然而，此算法未对协方差矩阵拟合误差进行考虑，仍以采样协方差矩阵替代INC矩阵计算波束形成权向量，因而在高输入信噪比(Signal-to-noise ratio,SNR)下性能损失严重。不同于K.L.Bell等人对期望信号DOA误差施加的均匀先验分布，C.J.Lam等人提出一种基于高斯DOA误差先验分布的贝叶斯波束形成算法(C.J.Lam,A.C.Singer,Bayesian beamforming for DOA uncertainty:Theory andImplementation,IEEE Trans.Signal Process.54(11)(2006)4435-4445.)，并采用傅里叶变换降低K.Bell所提算法实现过程中的计算复杂度。然而此算法仍是基于采样协方差矩阵计算贝叶斯权向量的，因而具有与K.L.Bell所提算法相同的性能局限。
发明内容
要解决的技术问题
为了避免现有技术的不足之处，本发明提出一种平稳干扰环境下的贝叶斯稳健波束形成方法。本发明面向模型失配、小样本、先验信息缺失等信号环境中的阵列波束形成需求，通过引入贝叶斯机器学习技术，建立基于观测数据多级概率模型的阵列处理理论与方法体系，并将其应用于解决该领域中的平稳环境干扰抑制这一关键问题，以突破阵列信号处理理论进一步发展完善所遇到的瓶颈、应对现有阵列信号处理方法应用于解决日趋复杂信号环境中的阵列波束形成问题时所受到的挑战。
技术方案
由于贝叶斯概率模型可灵活描述阵列接收信号中各分量的统计特性，因而有利于准确提取期望信号的波形信息。本发明利用变分推断(Variational inference)准则估计分层概率模型中各隐变量的后验分布，并根据估计结果修正预设的期望信号导向矢量，重构INC矩阵，以计算贝叶斯权向量。
假设期望信号和干扰信号均为远场窄带信号，其入射到一个由N个阵元组成的均匀线阵上，且阵元间距为入射信号波长的一半。当利用该线阵的接收信号进行自适应波束形成时，本发明所采用的技术方案包括以下步骤：
步骤1：建立阵列接收信号模型；
步骤2：建立贝叶斯概率模型，具体内容如下：
假设阵列接收信号模型中的各采样时刻的干扰+噪声向量nk相互独立，且服从均值为0，精度矩阵为Λ的复高斯分布，则nk的概率密度函数可表示为
式中，Λ为N×N维矩阵；
根据以上概率模型，可知阵列观测数据的似然函数服从如下多变量高斯分布
p(Y|s，v，Λ)＝π-NK|Λ|Ketr{-(Y-vsH)HΛ(Y-vsH)}
式中，Y为由K个采样时刻的阵列接收数据yk，k＝1，…，K组成的N×K维矩阵，s＝[s1，…，sK]T为由K个快拍的期望信号波形幅度信息组成的向量，etr{·}表示对矩阵迹求指数操作；
假设精度矩阵Λ服从如下Wishart先验分布
p(Λ)＝W(Λ|W，v)∝|Λ|v-Netr{-W-1Λ}
式中，∝表示正比于，W表示均值矩阵，v表示自由度；
期望信号波形幅度向量s内各元素假设相互独立，均服从如下均值为0，精度为的复高斯分布
式中，服从如下伽马先验分布
式中，Γ(·)表示伽马函数，a＞0表示形状参数，b＞0表示伸缩参数；
假设期望信号导向矢量v服从如下复Watson分布
式中，μ∈CSN为复单位球中的标称导向矢量，为聚焦参数；此分布关于μ空间对称；正则因子cp(λ)可表示为
式中，表示kummer汇合超几何函数，表示升因子；
令其先验概率服从如下Wishart分布
式中，W2表示均值矩阵，v2表示自由度；
用集合表示概率模型中全部未知变量，并将其统计特性与观测数据Y的分布结合起来得到如下联合概率密度
步骤3：使用变分推断方法，得到Θ中各变量的后验分布参数的更新公式，具体内容如下：
(1)初始化概率模型参数，即令v＝vp，其中vp为预设的期望信号导向矢量，v＝v2＝N，a＝b＝10-6，s＝1K×1，W＝W2＝10+6·IN×N，其中1K×1为K×1维的全1向量，IN×N为N×N维的单位矩阵；
(2)定义ln f(v)＝ln[p(v)/cp(λ)]并将此式在模v0处进行泰勒展开，得到
式中，为N×N维的海塞矩阵；将上式两端取指数并进行正则化，可知v近似服从如下高斯先验分布
式中，模v0可依牛顿下降法迭代计算，第t+1步的迭代公式为
式中，梯度和海塞矩阵的计算公式分别为
(3)期望信号导向矢量v的后验概率q(v)服从复高斯分布，其均值向量μv和精度矩阵Λv的迭代更新公式为
式中，<·>表示求数学期望操作；
(4)定义其后验概率服从Wishart分布，即此分布中自由度v3和尺度矩阵W3的更新公式为
v3＝v2+1
(5)期望信号波形幅度向量s的后验概率q(s)服从复高斯分布，其均值向量μs和精度矩阵Λs的迭代更新公式为
(6)干扰+噪声精度矩阵Λ的后验概率q(Λ)服从Wishart分布，其尺度矩阵W1和自由度v1的迭代更新公式分别为
v1＝K+v
(7)期望信号精度的后验概率服从伽马分布，即其中各分布参数的迭代更新公式为
a2＝K+a
(8)迭代(2)-(7)过程，直至满足收敛条件；
步骤4：将干扰+噪声精度矩阵及期望信号导向矢量的收敛解Λopt和vopt代入公式计算得到最优权矢量wopt。
本发明进一步的技术方案为：步骤1具体如下：
假设阵列采样快拍数为K，则各采样时刻的阵列接收数据可表示为：
式中，yk、v和nk均为N×1维向量，分别表示采样时刻k的阵列输出信号，目标分量及干扰+噪声分量；sk为采样时刻k下期望信号波形的幅度；目标、干扰及噪声分量假设相互独立；v为导向矢量，其表达式为其中λ0为入射信号波长，{d1，d2，…，dN}为阵元位置坐标，θ为期望信号的DOA。
本发明进一步的技术方案为：步骤3中的收敛条件为其中为当前步计算出的期望信号精度值，为上一步计算出的期望信号精度值。
一种计算机系统，其特征在于包括：一个或多个处理器，计算机可读存储介质，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述的方法。
一种计算机可读存储介质，其特征在于存储有计算机可执行指令，所述指令在被执行时用于实现上述的方法。
一种计算机程序，其特征在于包括计算机可执行指令，所述指令在被执行时用于实现上述的方法。
有益效果
本发明提出的一种平稳干扰环境下的贝叶斯稳健波束形成方法，依据贝叶斯机器学习技术在对混合信号中各个分量特征重构性能方面所具有的优势，建立以高精度空域滤波器参数估计为目的的稳健波束形成理论框架。该框架中引入了对接收信号各分量空域特征进行概率建模的步骤，其重要性体现在以下两个方面：对观测数据的分层概率建模过程实质上是利用一组合理先验分布达到对观测数据的最佳拟合，这一过程较好地利用了观测数据中各信号分量的结构特征，具有与极大似然类方法类似的原理，但高效的参数估计算法的使用能够显著改善波束形成过程的计算效率；各信号分量的结构信息包含在用于空域滤波器参数估计的分层概率模型中，该模型中先验分布参数能够依据实际信号环境进行自适应调整，从而使阵列波束形成方法更好地适应不同信号环境中的处理需求。有益效果如下：
(1)由于现有贝叶斯波束形成方法对期望信号导向矢量的高斯分布假设与实际概率模型并不吻合，因而难以获得最优的目标阵列流形估计精度。为此，本发明引入依赖于入射信号方位的Waston分布实现对期望信号导向矢量的精确估计，随后对干扰和噪声成分施加合理的先验分布假设，并利用观测数据对概率参数进行优化，以得到与真实模型相吻合的分布函数，从而有效地提高对干扰+噪声协方差矩阵的重构精度。相关方法具有显著优于对角加载类和已有贝叶斯类方法的阵列结构适应能力和输出SINR。
(2)针对现有贝叶斯波束形成方法在非共轭先验分布模型中无法推断参数的问题，本发明提出基于局域变分推断的滤波器参数估计方法，从而获得比传统概率采样方法更高的计算效率。具体地，将期望信号导向矢量的Watson先验在其模点处进行二阶泰勒展开，得到一高斯形式的概率密度下界，由于此近似先验分布与阵列接收数据似然函数满足共轭关系，故可应用变分准则精确推断出后验分布形式。上述概率近似与变分推断过程合称为局域变分推断方法。通过不断迭代最小化精确先验与近似先验之间的间隙，可获得与传统概率采样方法相当的计算精度。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图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}为阵元位置坐标，θ为期望信号的DOA。
步骤三：根据步骤二中建立的阵列接收信号模型，建立分层概率框架，包括以下子步骤：
子步骤一：假设各采样时刻的干扰+噪声向量nk相互独立，且服从均值为0，精度矩阵为Λ的复高斯分布，则nk的概率密度函数可表示为
式中，Λ为N×N维矩阵。
根据以上概率模型，可知阵列观测数据的似然函数服从如下多变量高斯分布：
p(Y|s，v，Λ)＝π-NK|Λ|Kctr{-(Y-vsH)HΛ(Y-vsH)}
式中，Y为由K个采样时刻的阵列接收数据yk，k＝1，…，K组成的N×K维矩阵，s＝[s1，…，sK]T为由K个快拍的期望信号波形幅度信息组成的向量，etr{·}表示对矩阵迹求指数操作；
子步骤二：假设精度矩阵Λ服从如下Wishart先验分布
p(Λ)＝W(Λ|W，v)∝|Λ|v-Netr{-W-1Λ}
式中，∝表示正比于，W表示均值矩阵，v表示自由度；
子步骤三：期望信号波形幅度向量s内各元素假设相互独立，均服从如下均值为0，精度为的复高斯分布
式中，服从如下伽马先验分布
式中，Γ(·)表示伽马函数，a＞0表示形状参数，b＞0表示伸缩参数；
子步骤四：假设期望信号导向矢量v服从如下复Watson分布
式中，μ∈CSN为复单位球中的标称导向矢量，为聚焦参数。此分布关于μ空间对称。正则因子cp(λ)可表示为
式中，表示kummer汇合超几何函数，表示升因子。
令其先验概率服从如下Wishart分布
式中，W2表示均值矩阵，v2表示自由度；
子步骤五：用集合表示概率模型中全部未知变量，并将其统计特性与观测数据Y的分布结合起来得到如下联合概率密度
步骤四：利用变分推断方法计算Θ中各变量的后验分布参数，包括以下子步骤：
子步骤一：初始化概率模型参数，即令v＝vp，其中vp为预设的期望信号导向矢量，v＝v2＝N，a＝b＝10-6，s＝1K×1，W＝W2＝10+6·IN×N，其中1K×1为K×1维的全1向量，IN×N为N×N维的单位矩阵；
子步骤二：定义ln f(v)＝ln[p(v)/cp(λ)]并将此式在模v0处进行泰勒展开，得到
式中，为N×N维的海塞矩阵。将上式两端取指数并进行正则化，可知v近似服从如下高斯先验分布
式中，模v0可依牛顿下降法迭代计算，第t+1步的迭代公式为
式中，梯度和海塞矩阵的计算公式分别为
子步骤三：期望信号导向矢量v的后验概率q(v)服从复高斯分布，其均值向量μv和精度矩阵Λv的迭代更新公式为
式中，<·>表示求数学期望操作；
子步骤四：定义其后验概率服从Wishart分布，即此分布中自由度v3和尺度矩阵W3的更新公式为
v3＝v2+1
子步骤五：期望信号波形幅度向量s的后验概率q(s)服从复高斯分布，其均值向量μs和精度矩阵Λs的迭代更新公式为
子步骤六：干扰+噪声精度矩阵Λ的后验概率q(Λ)服从Wishart分布，其尺度矩阵W1和自由度v1的迭代更新公式分别为
v1＝K+v
子步骤七：期望信号精度的后验概率服从伽马分布，即其中各分布参数的迭代更新公式为
a2＝K+a
子步骤八：迭代子步骤二至七过程，直至满足收敛条件，即其中为当前步计算出的期望信号精度值，为上一步计算出的期望信号精度值。
步骤五：将干扰+噪声精度矩阵及期望信号导向矢量的收敛解Λopt和vopt代入公式计算得到最优权矢量wopt。
利用计算机进行数值仿真，检验本发明所提方法的估计性能。
仿真采用阵元间距为半波长的10元均匀线阵，各阵元接收噪声服从均值为0，方差为1的复高斯分布，且相互独立。阵列入射信号总数为3，其中期望信号数为1，等功率干扰信号数为2，且各入射信号波形均服从复高斯分布。干噪比(Interference-TO-noise ratio，INR)设为30dB。干扰信号入射DOA分别为30°和50°，期望信号入射DOA的预设值为5°。蒙特卡洛实验次数设为100，每次实验中，各阵元位置误差服从[-0.05λ0，0.05λ0]范围内的均匀分布，其中λ0为入射信号波长。
1)存在DOA估计误差时，比较本发明所提方法与现有五种稳健波束形成方法的输出SINR结果
考虑由DOA估计误差引致的期望信号导向矢量失配情形，即每次蒙特卡洛实验中，期望信号DOA的预设值与真实值间的误差服从[-3°，3°]范围内的均匀分布。阵列采样快拍数为50。利用所提方法与现有五种方法，即WORST CASE 1方法(S.A.Vorobyov，A.B.Gershman，Z.-Q.Luo，Robust adaptive beamforming using worst-caseperformance optimization：a solution to the signal mismatch problem，IEEETrans.Signal Process.51(2)(2003)313-324.)，WORST CASE 2方法(A.Khabbazibasmenj，S.A.Vorobyov，A.Hassanien，Robust adaptive beamforming based on steering vectorestimation with as little as possible prior information，IEEE Trans.SignalProcess.60(6)(2012)2974-2987.)，SPSS-INC方法(Z.Zhang，W.Liu，W.Leng，A.Wang，H.Shi，Interference-plus-Noise covariance matrix reconstruction via spatialpower spectrum sampling for robust adaptive beamforming，IEEE SignalProcess.Lett.23(1)(2016)121-125.)，BAYESIAN 1方法(K.L.Bell，Y.Ephraim，H.L.VanTrees，A Bayesian approach to robust adaptive beamforming，IEEE Trans.SignalProcess.48(2)(2000)386-398.)，BAYESIAN 2方法(C.J.Lam，A.C.Singer，Bayesianbeamforming for DOA uncertainty：Theory and Implementation，IEEE Trans.SignalProcess.54(11)(2006)4435-4445.)进行波束形成，绘制输出SINR曲线。图1为六种方法的波束形成结果与最优SINR对比图。从图1所示结果可知，所提方法的阵列输出SINR最接近最优输出SINR，因而其干扰抑制性能优于现有方法。
2)存在相干局域散射时，比较本发明所提方法与现有五种稳健波束形成方法的输出SINR结果
真实期望信号导向矢量可表示为其中为预设的期望信号导向矢量，表示相干散射路径，ηi，i＝1，2，3，4为各路径上对应的相位。各次蒙特卡洛实验中，参数{θi}，{ηi}独立同分布，其中{θi}服从均值为3°，方差为1°的高斯分布，{ηi}服从[0，2π]内的均匀分布。设置快拍数为50，选择worst case 1、worst case 2、SPSS-INC、Bayesian 1和Bayesian 2五种方法作为本发明所提方法的比较对象，分别计算每种方法在不同输入SNR条件下获得的波束输出SINR。图2所示为上述六种方法的输出SINR及最优SINR随输入SNR的变化曲线。从图2所示结果可知，所提方法具有优于其他现有方法的干扰抑制性能，反映了所提方法可有效利用分层概率模型避免不同信号分量之间的串扰，在实现准确的期望信号波形重构的同时较好地抑制了干扰分量。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明公开的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。
