幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法与流程

本发明属于阵列信号处理领域，具体涉及一种幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法。

背景技术：

波达方向估计是阵列信号处理的一个重要研究方向并在军用和民用领域发挥着广泛的作用。众所周知，对于一个具有l个阵元的均匀线性阵列，采用music、esprit等经典波达方向估计方法最多只能有效估计l-1个信源。然而在实际的信源定位当中，往往存在着信源数大于阵元数，即欠定的情况，因此如何以更少的阵元数估计更多的信源数已成为信源定位领域的研究热点并吸引了国内外研究学者的广泛关注。解决该问题的一种有效方案是采用非均匀阵列在统计量或累积量域形成虚拟阵列，从而以更高的自由度数完成欠定情况下的波达方向估计。目前，常用的非均匀阵列为嵌套阵列和互质阵列。

基于嵌套阵列，p.pal和p.p.vaidyanathan于2010年在《nestedarrays:anovelapproachtoarrayprocessingwithenhanceddegreesoffreedom》提出了空间平滑music方法，该方法虽可实现欠定情况下的波达方向估计，但在构建满秩矩阵时没能用到虚拟阵列的所有阵元，从而导致一定的信息损失。

基于互质阵列，s.qi，y.d.zhang和m.g.amin于2015年在《generalizedcoprimearrayconfigurationsfordirection-of-arrivalestimation》提出了稀疏重构类欠定波达方向估计方法，该方法存在的主要问题是计算量巨大。

同时需要进一步指出的是，上述所有方法均只适用于理想的阵列接收条件。然而实际中由于本地振荡器时钟偏移和传感器位置误差等的影响，阵列接收数据往往存在一定的幅度和相位误差，即幅相误差。在这种背景条件下，已有欠定波达方向估计方法性能均会出现明显下降甚至失效。因此，研究幅相误差下的有效欠定波达方向估计方法对提高实际探测系统性能具有重要意义。

事实上，为了降低阵列幅相误差的影响，近年来研究学者们也相继提出了一些有效的方法，包括s.g.cao和z.f.ye等人在《ahadamardproductbasedmethodfordoaestimationandgain-phaseerrorcalibration》提出的hadamardproduct方法，j.liu，w.d.zhou等人在《covariancematrixbasedfastsmoothedsparsedoaestimationwithpartlycalibratedarray》提出的快速sl0方法等，但这些方法均不适用于欠定情况下的波达方向估计，同时大多需要大量的角度栅格搜索，复杂度高，且对非高斯白噪声非常敏感，噪声鲁棒性较弱。

技术实现要素：

本发明克服了现有技术中的不足，提供一种幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法，其用于解决现有波达方向估计方法在阵列存在幅相误差及信源数大于阵元数(即欠定)时的失效问题。

本发明的技术方案如下：

一种幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法，其包括如下步骤：

步骤一：根据给定阵列时域接收数据和信源数码，以形成最大连续虚拟孔径为原则，选择阵列输出数据，构建四阶累积量向量r；

步骤二：构建阈值判断函数f和以及连乘子函数f1，f2，在参考信源交替启停工作下联合应用阈值判决和连乘子函数在采样时间段为2q内完成阵列幅相误差估计；在前q采样时间段内，停止参考信源工作，当阈值判断函数f>0且进行代数运算完成幅相误差估计；当阈值判断函数不满足f>0且启动参考信源工作，并在后q采样时间段内进行代数运算完成幅相误差估计，其中q为时间长度，q>0,其单位为秒；

步骤三：利用幅相误差估计结果，修正四阶累积量向量r，得到并基于子阵划分策略获得toeplitz满秩矩阵；

步骤四：通过esprit算法完成波达方向估计。

优选地，步骤一中构建所述四阶累积量向量r的步骤如下：

s1：预先设置k个具有非零四阶峰值的信源入射到由2l个阵元组成的嵌套阵列上，阵元位置为{0d,1d,…,ld,2ld,…,l²d}，其中标号位置为{0d,1d,2d,…,ld}的阵元对应于子阵1，标号位置为{0d,ld,2ld,…,l²d}的阵元对应于子阵2，且标号位置为{0d,1d,ld}的阵元已校准，d＝λ/2，λ为载波波长；

s2：设x(t)为阵列时域接收数据，则第m,n,p,q∈{1,…,2l}个阵元观测数据的四阶累积量为：

其中，函数cum{g}代表进行四阶累积量运算；pm和γm分别代表第m(∈{1,…,2l})个阵元的位置和幅相误差，pn和γn分别代表第n(∈{1,…,2l})个阵元的位置和幅相误差，为γn的共轭；pp和γp分别代表第p(∈{1,…,2l})个阵元的位置和幅相误差；pq和γq分别代表第q(∈{1,…,2l})个阵元的位置和幅相误差，为γq的共轭；

xm(t)代表x(t)的第m∈{1,…,2l}行数据；

xn(t)代表x(t)的第n∈{1,…,2l}行数据，为对xn(t)的各组成元素进行共轭运算后的结果；

xp(t)代表x(t)的第p∈{1,…,2l}行数据；

xq(t)代表x(t)的第q∈{1,…,2l}行数据，为对xq(t)的各组成元素进行共轭运算后的结果；是第k(1≤k≤k)个信源信号的四阶累积量，sk代表第k个信源信号的复幅度；

φk＝2πsin(θk)/λ；其中φk表示基础相位差，θk为第k(1≤k≤k)个信源信号的入射角度；

s3：基于嵌套阵列，以形成最大连续虚拟孔径为原则，选择m,n,p,q值形成四阶累积量向量r，表示为

r＝fac

其中，f为对应于第m,n,p,q个阵元幅相误差积的对角矩阵；a＝[a(θ1),…,a(θk)]代表虚拟导向矩阵，代表对应于第k(1≤k≤k)个信源信号的虚拟导向矢量，l为实际阵元数的一半，4l²+1为形成的虚拟阵元数，d为阵元间距，满足d＝λ/2，λ为载波波长，φk＝2πsin(θk)/λ，其中φk表示基础相位差，θk为第k(1≤k≤k)个信源信号的入射角度，是第k(1≤k≤k)个信源信号的四阶累积量，sk代表第k个信源信号的复幅度，上标t表示转置操作。

优选地，在步骤二中所述阈值判决函数定义为

f＝q2-τ＝|c2,1,1,1+c2,1,2,2|/2-τ

其中，f为对应于子阵1的阈值判决函数，为对应于子阵2的阈值判决函数，

且

为用户设定的判决阈值，|g|表示取绝对值操作，表示对k个元素进行求和运算，是第k(1≤k≤k)个信源信号的四阶累积量，sk代表第k个信源信号的复幅度；

连乘子函数定义为

其中，f1和f2为对应于子阵1的连乘子函数，和为对应于子阵2的连乘子函数，且

代表对括号内k1个元素进行乘积运算。

优选地，在步骤二中完成阵列幅相误差估计的具体步骤为：

a)若f＞0且则子阵1和子阵2的幅相误差分别估计为

其中，和分别为第k1+1和l+k2+1个阵元的幅相误差，f1和f2为对应于子阵1的连乘子函数，和为对应于子阵2的连乘子函数，q2＝|c2,1,1,1+c2,1,2,2|/2，q3＝(c1,1,1,1+c2,2,2,2+cl+1,l+1,l+1,l+1)/3；且

b)若f≤0或则启动参考信源发射参考信号，并在下一个采样周期基于上述相同的处理方式完成阵列幅相误差估计。

优选地，在步骤三中获得toeplitz满秩矩阵的具体步骤为：

设代表f的估计值，取逆处理后得到修正的四阶累积量向量表示为

其中，为的逆矩阵，a代表虚拟导向矩阵，

将分成2l²+1份得到满秩矩阵

其中，0m×n代表m×n的零矩阵，代表(2l²+1)×(2l²+1)的单位矩阵。

规定反对角矩阵为

形成toeplitz满秩矩阵

优选地，选择m,n,p,q值形成四阶累积量向量，其具体步骤如下：

其中，r(l)代表四阶累积量向量的第l个元素。

与现有技术相比，本发明的有益效果是：

本发明方法通过联合利用嵌套阵列和四阶累积量来扩展阵列虚拟孔径，不仅具有更好的噪声鲁棒性，而且仅用2l个阵元有效估计2l²个信源。

本发明提供的方法在一个参考信源适时工作下，联合应用阈值判决和连乘子函数完成阵列幅相误差估计，无需噪声方差和信源数的先验信息，且求解过程仅是简单的代数运算，复杂度低、运行时间短。

本发明提供的方法通过子阵划分获得toeplitz满秩矩阵，进一步基于esprit获得波达方向估计的闭式解，有效避免了大运算量的角度栅格搜索。

附图说明

图1是本发明采用的部分校准嵌套传感器阵列的结构图；

图2是本发明提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法的流程图；

图3是本发明提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法可估计最大信源数的结果图；

图4a是本发明提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法在高斯白噪声下的均方根误差随信噪比变化比较图；

图4b是本发明提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法在高斯色噪声下的均方根误差随信噪比变化比较图；

图4c是本发明提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法在未知非均匀噪声下的均方根误差随信噪比变化比较图；

图5是本发明提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法的单次平均运行时间随采样快拍数变化比较图。

具体实施方式

以下结合附图对本发明进行详细描述。

参照图1和图2，本发明的实现步骤如下：

步骤一：合理选择阵列输出数据构建四阶累积量向量；

假设k个具有非零四阶峰值的信源入射到由2l个阵元组成的嵌套阵列上，如图1所示。阵元位置为{0d,1d,…,ld,2ld,…,l²d}，其中标号位置为{0d,1d,…,ld}的阵元对应于子阵1，标号位置为{0d,ld,…,l²d}的阵元对应于子阵2，且标号位置为{0d,1d,ld}的阵元已校准，d＝λ/2，λ为载波波长。设x(t)为阵列时域接收数据，则第m,n,p,q∈{1,…,2l}个阵元观测数据的四阶累积量可计算为：

其中，pm和γm分别代表第m(∈{1,…,2l})个阵元的位置和幅相误差，是第k(1≤k≤k)个信源信号的四阶累积量，φk＝2πsin(θk)/λ；

上标*为共轭操作，函数cum{g}代表进行四阶累积量运算。

基于嵌套阵列的性质，选择合适的m,n,p,q值形成具有如下特性的四阶累积量向量r，表示为

r＝fac

其中，a＝[a(θ1),…,a(θk)]，f为对应于第m,n,p,q个阵元幅相误差积的对角矩阵。

上标t表示转置操作。

四阶累积量向量r的详细形成过程为

步骤二：在参考信源适时工作下联合应用阈值判决和连乘子函数完成阵列幅相误差估计；

定义阈值判决函数为：

f＝q2-τ＝|c2,1,1,1+c2,1,2,2|/2-τ

其中，为用户设定的判决阈值，|g|表示取绝对值操作。

定义连乘子函数为：

若f＞0且则子阵1和子阵2的幅相误差分别估计为

其中，q3＝(c1,1,1,1+c2,2,2,2+cl+1,l+1,l+1,l+1)/3。

若f≤0或则启动参考信源发射参考信号，并在下一个采样周期基于上述相同的处理方式完成阵列幅相误差估计。

步骤三：利用幅相误差估计结果修正四阶累积量向量并基于子阵划分策略获得toeplitz满秩矩阵；

设代表f的估计值，取逆处理后得到修正的四阶累积量向量表示为

将分成2l²+1份得到满秩矩阵

其中，0m×n代表m×n的零矩阵，代表(2l²+1)×(2l²+1)的单位矩阵。

规定反对角矩阵为

形成toeplitz满秩矩阵

在本发明的具体实施例中，toeplitz满秩矩阵可表示为

其中，

上标h为共轭转置操作，diag(·)代表向量对角化操作。

对进行特征值分解获得(2l²+1)×k的信号子空间es，根据esprit算法的核心思想，一定存在一个非奇异值矩阵t，满足

将es的前2l²行和后2l²行分别用es1和es2表示，则有

其中，上标代表伪逆运算。

对进行特征值分解并获得特征值qk,k＝1,k,k，则波达方向可估计为

其中，arg(·)代表取角操作。

下面通过仿真实验分析本发明所提出的幅相误差下基于部分校准嵌套阵列的欠定波达方向估计方法的估计性能及计算有效性。仿真过程均通过matlab软件进行，计算机配置为英特尔酷睿i3cpu，3.4g主频，4g内存。

仿真实验1：该实验用以论证本发明所提方法用2l个阵元有效估计2l²个信源的能力。总阵元个数为2l＝4，8个信源由{-60°,-43°,-26°,-9°,8°,25°,42°,59°}入射到阵列，信噪比为20分贝，采样快拍数为2000，阵列幅度误差和相位误差由和产生，且ξl和ζl为均匀分布在[-0.5,0.5]范围内的随机变量。仿真结果如图3所示，可以看出本发明方法估计的结果与真实值仅存在很小的偏差，证明了本发明方法确实可以利用2l个阵元有效估计2l²个信源。

仿真实验2：该实验用以分析本发明所提方法在高斯白噪声/色噪声和未知非均匀噪声下的波达方向估计性能。总阵元个数为2l＝10，两个信源由{-15.2°,9.8°}入射到阵列，采样快拍数为2000，信噪比由0分贝变化到20分贝，阵列幅相误差与仿真实验1相同。图4a～图4c依次为高斯白噪声、高斯色噪声和未知非均匀噪声下的波达方向估计均方根误差随信噪比的变化关系图，其中高斯色噪声方差满足

未知非均匀噪声方差为

n＝diag([0.53.81.28.02.30.75151.95.22.2]^t)

由仿真结果可以看出，本发明所提方法性能明显优于对比方法，且对阵列幅相误差和各种噪声均具有很好的鲁棒性。

仿真实验3：该实验用以评价本发明所提方法的计算复杂度。总阵元个数为2l＝10，两个信源由{-15.2°,9.8°}入射到阵列，信噪比固定为20分贝，采样快拍数由200变换到2000，阵列幅相误差与仿真实验1相同。用以对比的hadamardproduct方法的角度栅格搜索步长为0.1，搜索范围为[-30°,30°]。仿真所用时间如图5所示，可以发现本发明所提供的方法的运行时间由于无需角度栅格搜索而远远低于hadamardproduct方法。

上述仿真实验结果充分证明了本发明所提供的相位误差下的波达方向估计方法的有效性、鲁棒性和高效性。

可以理解的是，以上实施方式仅仅是为了说明本发明的原理而采用的示例性实施方式，然而本发明并不局限于此。本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。