一种稳健的波达方向估计方法
【技术领域】
[0001] 本发明属于阵列信号处理技术领域，设及一种稳健的波达方向估计方法。
【背景技术】
[0002] 波达方向值irection-of-Arrival，D0A)是所接收来波信号的重要参数。利用一 组传感器阵列接收数据实现D0A估计，在雷达、无线通信、声响、语音信号处理、地震波检测 等阵列信号处理相关应用领域起着重要的作用，是信号处理领域的热点研究课题。
[0003]纵观大量现存的D0A估计方法，MUSIC(multiplesignalclassification)方法和 ESP民IT(estimationofsignalparametersviarotationalinvariancetechnique)方 法是最常用的DOA估计方法。MUSIC系列方法利用信号子空间与噪声子空间的正交性获得 MUSIC谱的峰值实现波达角度估计。但是，MUSIC方法无法获得解析解。ESPRIT系列方法将 传感器阵列划分为两个子阵，第一个子阵的信号子空间特征向量是第二个子阵的信号子空 间特征向量与包含D0A信息的旋转矩阵的乘积，通过计算旋转矩阵获得D0A估计的解析解。 但是，仅仅利用一个相关矩阵的策略限制了ESPRIT方法的稳健性及估计精度。为了提高 算法的稳健性，文献[3]提出了一种利用解盲信号分离问题的二阶盲辨识（second-order blindidentification,S0BI)实现DOA估计的新方法，本文将其称为S0BI-D0A方法，多个 相关矩阵的利用提高了该算法的稳健性和分辨能力。但是，S0BI-D0A算法要求目标矩阵中 至少有一个正定矩阵构建预白化矩阵进行预白化操作，但是，噪声的影响和有限样本数 的限制导致此目标矩阵的正定性难W保证，因此白化操作难免引入误差，而且这些在预白 化处理阶段引入的误差无法在后续处理中得W校正，这就限制了S0BI-D0A算法的估计精 度。
[0004]W下是本发明引用的相关参考文献：
[000引 [l]VinodVR，BoonP.N，AndyWΗK.Insi組tsIntoMUSIC-LikeAlgoi'ithmU]· IEEETransactionsonSignalProcessing，2013,61 (10) :2551-2556.
[0006] [2]TripathyP，SrivastavaSC，Sin組SN.AModifiedTLS-ESPRIT-BasedMethod forLow-FrequencyModeIdentificationinPowerSystemsUtilizingSynchrophasor Measurements[J].IEEETransactionsonPowerSystems，2011，26(2) :719-727.
[0007] [3]PourrostarnJ，ZekavatSA，PourkhaatounM.Super-resolution direction-of-arrivalestimationviablindsignalseparationmethods[J].IEEE RadarConference，2007,614-617.
[0008] [4]YeredorA.Non-orthogonaljointdiagonalizationintheleast-squares sensewithapplicationinblindsourceseparation[J].IEEETransactionson SignalProcessing,2002,50(7) :1545-1553.
[0009] [5]FengDZ，ZhengWX，CichockiA.Matrix-groupalgorithmviaimproved whiteningprocessforextractingstatisticallyindependentsourcesfromarray signals[J].IEEETransactionsonSignalProcessing，2007,55 (3) :962-977.
[0010] [6]Lee皿，WengrovitzM. Statistical characterization of the MUSIC null spectrum[J].IEEE Transactions on Signal Processing，1991，39 (6) :1333-1347.
【发明内容】
[0011] 针对上述现有技术中存在的问题或缺陷，本发明的目的在于，提供一种实现对波 达方向更加稳健、精度更高且解析的估计，有效解决波达方向问题的波达方向估计方法。
[0012] 为了实现上述目的，本发明采用如下技术方案：
[0013] 一种稳健的波达方向估计方法，具体包括W下步骤：
[0014]步骤一：根据阵元接收信号，给出传感器阵列的接收数据；
[001引步骤二：构造目标矩阵组C ;
[0016] 步骤Ξ:对目标矩阵组C进行降维处理，得到降维后的目标矩阵组C；
[0017] 步骤四：对目标矩阵组日进行实对称化，将复数域的目标矩阵组^>转化为实数域对 称矩阵护；
[0018] 步骤五：根据实数域对称矩阵护，采用非正交联合对角化-列交替更新与对角化 算法估计波达方向角。
[0019] 具体地，所述步骤一中的传感器阵列的接收数据为：
[0020] x(t) = ，…，XM(t)]T= A( Θ )s(t)+n(t)
[0021] 其中，Xi(t)，…，XmW分另懐示第1~Μ个阵元接收的信号；A(目）=[日（目1)，···， a( Θ J]表示列满秩的阵列流型，其列向量
v=l，·，·，Ν 为导向矢量，是第ν(ν=1，…，脚个信源到第u(u=l，···，]?)个阵元的波达角；s(t)表示均值为零，两两相互独立的源信号；n(t)表示所接收的均值为零、方差为σ"2且相互独 立的阵元噪声；t，t=l，···，Τ表示快拍。
[0022] 具体地，所述步骤二中的目标矩阵组C为：
[0023]
[0024]其中，
V= 1，…， Ν，阵元间隔为d，d《λ/2, λ表示信号波长；〇Ρ，ρ=1，-''^Ρ表示对角阵。
[00巧]具体地，所述步骤Ξ中的降维后的目标矩阵组C为：
[0026]
[0027] 其中，
I对Φ进行奇异值分解，取Φ的Ν个主特征矢量 组成降维矩阵Γ= [Ui，U2,…，叫].ΓΗ为的共辆转置，ii"为的共辆转置。
[0028] 具体地，所述步骤四中的实数域对称矩阵为：
[0029]
[0030]其中，
&p，f=i，一，2f为对角阵。
[0031] 具体地，所述步骤五的实现方式包括：
[0032] 建立表征联合对角化程度的代价函数，通过两个交替子步AC和DC实现代价函数 的最小化，得到对角阵此的估计值；
[0033] 根据公式（17)和（18)，利用对角阵1H勺估计值得到DP，P= 1，…，P:
[0036] 根据公式（1扣和
利用P= 1, . . . ,P得到Am,。似：
[0037] 0=[σ =Α,0·?^,Ρ = 1,...,p](13)
[0038]其中，Ρ= 4Κ:
[0039] 根据公式（12)，利用Am,。似得到Φ的估计值〇m(k)，m= 1，-..,4,k= 1，…， Κ:
[0040] Αι,ι化）=Rs化），Αι,2化）=氏化）巫Η，Αζ,ι化）=巫Rs化），Α2,2化）=巫Rs化）巫Η (12)
[0041] 波达方向角的统计平均值为：
[0042]
[0043] 其中，angle( ·)代表幅角主值；是〇m(k)的第V个对角元素，则D0A最终估 计为谷=(如...Λ)。
[0044] 与现有技术相比，本发明具有W下技术效果：
[0045] 1、本发明是一种解析的波达方向估计方法，无须后续对谱峰值的捜索求解步骤。
[0046] 2、本发明利用两个子阵接收数据的相关矩阵，构建了一组包含旋转矩阵信息的目 标矩阵，并据此实现对旋转矩阵的多次估计，进而实现对波达方向角统计平均值的求解，提 高了所提方法的稳健性和估计精度。
[0047] 3、本发明在求取旋转矩阵的过程中，无须引入预白化操作，避免了白化误差，提高 了所提方法的估计精度。
[004引 4、本发明所采用的非正交联合对角化算法-列交替更新与对角化（ACDC)，是性能 优异的解联合对角化问题的有效算法，本发明创新性地关注了其冗余分量，有效地解决了 旋转矩阵的估计问题。
【附图说明】
[0049] 图1是不同算法在不同信噪比下的PR曲线；
[0050] 图2是快拍数为Τ= 400时，不同算法在不同信噪比下的RISE曲线；
[005。图3是信噪比为SNR= 10地时，不同算法在不同快拍数下的RISE曲线。
[0052] 下面结合附图和实施例对本发明的方法做进一步详细地解释和说明。
【具体实施方式】
[0053] 遵从上述技术方案，本发明的稳健的波达方向估计方法，具体包括W下步骤：
[0054] 步骤一：根据阵元接收信号，给出传感器阵列的接收数据
[00巧]假设有Μ个全向阵元的等距线阵，阵元间隔为山（d《λ/2,λ表示信号波长）， 阵元接收Ν个入射方向为Θ= (01，…，ej的远场窄带信号，第U个阵元接收的信号可 表示为：
[0056]
(1)
[0057] 其中，Θ、是第v(v= 1，…，脚个信源到第u(u= 1，…，M)个阵元的波达角， Sy(t)是第V个信源发射的源信号，η。(t)是第U个阵元的阵元噪声。
[0058] 传感器阵列在快拍t，t= 1，…，T时所接收的数据为：
[0059] x(t) = [Xi(t)，…，XM(t)]T=A(Θ)s(t)+n(t) 似
[0060] 此时阵列流型A(θ) = [a(θ1)，…，a(θw)]是列满秩的，其列向量
v=l，·，·，Ν为导向矢量。s(t) = [Si(t)，···iSnW]τ表示均值为零，两两相互独立的源信号，即
[006dE[s(t)sH(t+kAτ)] =diag[p1似，…，PN(k)]=氏化）. 做
[00的]其中，diag[ ·]表示W向量?为对角线元素的对角阵。A脚=Ε[·α命成V= 1，…，Ν，k= 1，2，…，Κ，Δτ为非零延时步长。
[006引式似中n(t) = [ni(t)，…，前似^表示所接收的均值为零、方差为且相互独 立的阵元噪声，即
[0064]
p,q = \,...,M(4)
[006引其中，δP,为单位脉冲函数，当P=q时，δp,q= 1 ;当P声q时，δp,q= 0。为了 推导简便，在后续的公式推导中，不再体现出噪声的影响。
[0066] 步骤二：根据步骤一得到的传感器阵列的接收数据构造目标矩阵组C
[0067] 将步骤一中的传感器阵列分为两个子阵，分别包含的阵元编号为1，…，M-1和 2, -",Μ。则第一个子阵接收数据Xi(t)及第二个子阵接收数据X2(t)可W分别表示为： [006引 Xi(t) =[Xi(t)，…，Xm1(t)]T=A1(Θ)s(t) (5)
[0069] X2(t) = [X2(t)，…，XM(t)]T=Az(Θ)s(t) =Ai(Θ)巫s(t) 化）
[0070]其中，Ai (目）=[曰1(目1)，…，曰1(目N)]，
V = 1，·，·，Ν。Α2(Θ)的含