一种基于中心共轭零化的快速高精度阵列互耦自校正方法

1.本发明涉及一种基于中心共轭零化的快速高精度阵列互耦自校正方法，特别涉及存在阵元互耦情形下的阵列信号处理，属于电子信息技术领域，尤其是阵列信号处理技术领域。


背景技术：

2.阵列信号处理技术的一个分支是利用传感器阵列估计信号的波达方向，广泛应用于雷达、通信、声呐及医学诊断等军用和民用领域。基于参数建模的波达方向估计方法由于具有很高的分辨率而更受青睐。但这类方法的性能易受阵列的自身误差影响，例如阵元幅相误差，阵元位置误差以及互耦等。阵列互耦来源于波的传播特性，一般难以避免。因而，阵列的互耦校正问题受到广泛关注。
3.阵列互耦校正技术分为离线校正与在线自校正。离线校正一般在暗室内进行。通过在单径环境下离线测量阵列在多个方向的响应，计算出互耦校正系数，从而可在阵列实际应用时补偿互耦系数以实现互耦校正。这种方法可靠性高，但实现的设备成本和时间成本也高。在线自校正简称自校正，是在阵列进行波达方向估计时，同步估计出阵列互耦系数。这种方法成本较低，并且可跟踪因外部环境导致的阵列误差变化。因此这类方法近年来被广泛研究，并涌现出众多新的自校正算法。
4.文献1(潘玉剑,等,“一种基于零化约束的doa与互耦的联合估计方法”[p],202111172456.2,2021.12.03)提出了一种基于零化约束的互耦校正方法。该方法无需互耦自由度的先验信息，且计算量较低。但由于并未限定滤波器的零点位于单位圆上，导致所建模型与实际模型存在偏差，使得该方法的互耦系数估计精度较低。为了提高互耦的估计精度，文献2(y.pan,et al.,“joint estimation of doa and mutual coupling via imposing an annihilation constraint”[j],digital signal processing,130,103726,2022)提出，在零化滤波后，采用后处理的方法对互耦系数进行二次迭代估计，这种方法一定程度上提高了互耦系数的估计精度。但由于该方法需要两个迭代过程，使其一方面会累积误差，另一方面会增加计算复杂度，延长算法处理时间。另外文献2同样未限定滤波器的零点位于单位圆上，因而其参数估计精度仍有提升的空间。


技术实现要素：

[0005]
为了克服现有研究的不足，本发明提供了一种基于中心共轭零化的快速高精度阵列互耦自校正方法。
[0006]
一种基于中心共轭零化的快速高精度阵列互耦自校正方法的具体步骤如下：
[0007]
步骤一：阵列布阵及互耦存在时的信号建模：
[0008]
布设均匀线阵，阵元数量为m，相邻阵元间距为d；k个波长为λ，来自θ1,θ2,
…
,θk方向的窄带信号入射到均匀线阵，共接收到n个快拍，多快拍接收数据建模为：y＝cas+e，y＝[y[1]，y[2]，
…
，y[n]]，信号矩阵s＝[s[1]，s[2]，
…
，s[n]]，噪声矩阵e＝[ε[1]，ε[2]，
…
，ε
[n]]，c为互耦矩阵，c＝toeplitz([1,c1,c2,
…
,c
p
,0,
…
,0])，toeplitz(
·
)表示生成对称toeplitz矩阵，c1,c2,
…
,c
p
为互耦系数，p表示互耦自由度；
[0009]
其中，在第n个快拍时阵列的接收数据建模为：其中，在第n个快拍时阵列的接收数据建模为：其中，表示复数集合，流型矩阵a＝[a(θ1),a(θ2),
…
,a(θk)]，a(θ1),a(θ2),
…
,a(θk)表示导向矢量，第k个窄带信号导向矢量a(θk)的第m个元素为am(θk)＝exp[j(m-1)uk]，k＝1,2,
…
,k，m＝1,2,
…
,m，uk＝2πdcos(θk)/λ；s[n]为信号向量，ε[n]为噪声向量。
[0010]
步骤二：基于奇异值分解的数据维度压缩：
[0011]
对y进行奇异值分解，则y经维度压缩后得到vs表示由y的k个最大奇异值对应的右奇异矢量构成的矩阵。
[0012]
步骤三：基于中心共轭零化的优化问题建立：
[0013]
优化问题表示为：其中，||
·
||2表示2范数，z的向量化z＝vec(z)，vec(
·
)表示将矩阵按列堆叠向量化，ik表示k阶单位矩阵，表示kronecker积，互耦矩阵c
′
＝toeplitz([1,c
t
,0,...,0])，互耦向量c＝[c1,c2,...,c
p
′
]
t
，p
′
表示互耦自由度的估计值，辅助参数ηk为γ的第m(k-1)至第mk个元素构成的向量，k＝1,2,...,k，待优化向量υ＝[c
t
,γ
t
,h
t
]
t
；多快拍矩阵t-(γ)＝[t t
(η1),t t
(η2),
…
,t t
(ηk)]
t
，为构建toeplitz矩阵的算子，t(ηk)的第i行第j列元素为ηn的第i-j+k+1个元素，i＝1,2,
…
,(m-k)，j＝1,2,
…
,(k+1)，零化滤波器ω为常数向量；s.t.表示约束条件，(
·
)h、(
·
)
t
分别表示取共轭转置和转置，im(
·
)取对应变量的虚部，g为变换矩阵。im(ghh)＝0表征零化滤波器h具有中心共轭的性质。
[0014]
步骤四：基于实数sqp优化方法的优化问题求解：
[0015]
首先初始化其中c0＝0
p
′×1，0
p
′×1表示p
′
维零向量，γ的初始值γ0和h的初始值h0通过假设阵列无互耦并对z进行类tls-esptrit方法获得；然后进行迭代计算，在第q+1次迭代过程中得到υ
q+1
＝υq+δυ，q≥0，其中为第q次迭代计算得到的结果，δυ为搜索向量。其中i
p
′
表示p
′
阶单位矩阵，其余依此类推；δυ
′
为实数向量，通过求解如下实数线性二次优化问题得到：
[0016]
其中，f向量re(
·
)取对应变量的实部，第q次迭代得到的互耦矩阵b向量
d矩阵
[0017]
其中0
2mn
×
(k+1)
表示2mn行k+1列的零矩阵,其余依此类推；b矩阵其中，为构建toeplitz矩阵的另一种算子，r(hq)的第i行为j
k+1
为k+1阶反对角单位矩阵；多快拍q矩阵即为γq的第m(k-1)至第mk个元素构成的向量，[
·
]
·
,1:p
′
表示取矩阵的前p
′
列构成新的矩阵，toeplitz矩阵hankel矩阵hankel矩阵为的第1个元素，其余依此类推。
[0018]
步骤五：波达方向与互耦估计：
[0019]
迭代收敛后，得到收敛解迭代收敛后，得到收敛解为互耦向量的估计值，为γ的估计值，为h的估计值，将带入得到多项式，求解该多项式的k个根得到x1,x2,...,xk，则doa估计为，则doa估计为其中，∠(
·
)取复数的幅角主值；根据得到估计的互耦矩阵
[0020]
优选地，步骤三中的ω取值为其中，h0为步骤(4)迭代计算中h的初始值，为h0的复共轭。
[0021]
优选地，步骤三中的p
′
可等于或小于步骤(1)中的p。
[0022]
优选地，步骤三中变换矩阵g的设定方法如下：当k为奇数时，令n＝(k+1)/2，有其中，in表示n阶单位矩阵，jn表示n阶反对角单位矩阵；当k为偶数时，令n＝(k-1)/2，有
[0023]
与现有技术相比，本发明的有益效果在于：
[0024]
本发明通过赋予零化滤波器中心共轭的性质，从而让滤波器的零点位于单位圆上，使得所建模型与阵列信号处理模型相符，可提高互耦校正精度；其次，本发明方法无需
进行二次迭代，因而可提高算法执行速度。
附图说明
[0025]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0026]
图1为本发明基于中心共轭零化的快速高精度阵列互耦自校正方法的流程图；
[0027]
图2为本发明方法与其他方法在不同信噪比下波达方向估计比较示意图；
[0028]
图3为本发明方法与其他方法在不同信噪比下互耦估计性能比较示意图；
[0029]
图4为本发明方法与其他方法算法执行时间比较示意图；
[0030]
图5为本发明方法与其他方法的滤波器零点位置比较示意图。
具体实施方式
[0031]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0032]
如图1所示，一种基于中心共轭零化的快速高精度阵列互耦自校正方法，具体步骤如下：
[0033]
步骤一：阵列布阵及互耦存在时的信号建模：
[0034]
布设均匀线阵，阵元数量为m，相邻阵元间距为d；k个波长为λ，来自θ1,θ2,
…
,θk方向的窄带信号入射到均匀线阵，共接收到n个快拍，多快拍接收数据建模为：y＝cas+e，y＝[y[1]，y[2]，
…
，y[n]]，信号矩阵s＝[s[1]，s[2]，
…
，s[n]]，噪声矩阵e＝[ε[1]，ε[2]，
…
，ε[n]]，c为互耦矩阵，c＝toeplitz([1,c1,c2,
…
,c
p
,0,
…
,0])，toeplitz(
·
)表示生成对称toeplitz矩阵，c1,c2,
…
,c
p
为互耦系数，p表示互耦自由度；
[0035]
其中，在第n个快拍时阵列的接收数据建模为：其中，在第n个快拍时阵列的接收数据建模为：其中，表示复数集合，流型矩阵a＝[a(θ1),a(θ2),
…
,a(θk)]，a(θ1),a(θ2),
…
,a(θk)表示导向矢量，第k个窄带信号导向矢量a(θk)的第m个元素为am(θk)＝exp[j(m-1)uk]，k＝1,2,
…
,k，m＝1,2,
…
,m，uk＝2πdcos(θk)/λ；s[n]为信号向量，ε[n]为噪声向量，设定噪声为零均值高斯白噪声。
[0036]
具体的，互耦矩阵c可写为
[0037][0038]
步骤二：基于奇异值分解的数据维度压缩：
[0039]
对y进行奇异值分解，则y经维度压缩后得到vs表示由y的k个最大奇异值对应的右奇异矢量构成的矩阵。
[0040]
步骤三；基于中心共轭零化的优化问题建立：
[0041]
优化问题表示为：其中，||
·
||2表示2范数，z的向量化z＝vec(z)，vec(
·
)表示将矩阵按列堆叠向量化，ik表示k阶单位矩阵，表示kronecker积，互耦矩阵c
′
＝toeplitz([1,c
t
,0,...,0])，互耦向量c＝[c1,c2,...,c
p
′
]
t
，p
′
表示互耦自由度的估计值，p
′
可等于或小于p，即不要求互耦自由度精确的先验信息；辅助参数ηk为γ的第m(k-1)至第mk个元素构成的向量，k＝1,2,...,k，待优化向量υ＝[c
t
,γ
t
,h
t
]
t
；多快拍矩阵t-(γ)＝[t t
(η1),t t
(η2),
…
,t t
(ηk)]
t
，为构建toeplitz矩阵的算子，t(ηk)的第i行第j列元素为ηn的第i-j+k+1个元素，i＝1,2,
…
,(m-k)，j＝1,2,
…
,(k+1)，零化滤波器ω为常数向量，取值为其中，h0为下一步迭代计算中h的初始值，为h0的复共轭；s.t.表示约束条件，(
·
)h、(
·
)
t
分别表示取共轭转置和转置，im(
·
)取对应变量的虚部，g为变换矩阵。im(ghh)＝0表征零化滤波器h具有中心共轭的性质。所谓向量具有中心共轭的性质是指：向量的元素相对中心点共轭对称。即，对于h，有h＝j
k+1h*
。滤波器系数中心共轭对称是保证滤波器零点位于单位圆上的必要条件，使得所建模型与真实阵列信号处理模型相符。变换矩阵g的设定方法如下：当k为奇数时，令n＝(k+1)/2，有其中，in表示n阶单位矩阵，jn表示n阶反对角单位矩阵；当k为偶数时，令n＝(k-1)/2，有
[0042]
步骤四；基于实数sqp优化方法的优化问题求解：
[0043]
首先初始化其中c0＝0
p
′×1，0
p
′×1表示p
′
维零向量，γ的初始值γ0和h的初始值h0通过假设阵列无互耦并对z进行类tls-esptrit方法获得；然后进行迭代计算，在第q+1次迭代过程中得到υ
q+1
＝υq+δυ，q≥0，其中
为第q次迭代计算得到的结果，δυ为搜索向量。其中i
p
′
表示p
′
阶单位矩阵，其余依此类推；δυ
′
为实数向量，通过求解如下实数线性二次优化问题得到：
[0044]
其中，f向量re(
·
)取对应变量的实部，第q次迭代得到的互耦矩阵b向量d矩阵
[0045]
其中0
2mn
×
(k+1)
表示2mn行k+1列的零矩阵,其余依此类推；b矩阵其中，为构建toeplitz矩阵的另一种算子，r(hq)的第i行为j
k+1
为k+1阶反对角单位矩阵；多快拍q矩阵即为γq的第m(k-1)至第mk个元素构成的向量，[
·
]
·
,1:p
′
表示取矩阵的前p
′
列构成新的矩阵，toeplitz矩阵hankel矩阵hankel矩阵为的第1个元素，其余依此类推。
[0046]
所谓实数sqp优化方法是指：将非线性的复数优化问题转化实数迭代求解，每次迭代只需求解一个实数线性二次优化问题。对应的实数线性二次优化问题可通过拉格朗日数乘方法求解。
[0047]
步骤五；波达方向与互耦估计：
[0048]
迭代收敛后，得到收敛解迭代收敛后，得到收敛解为互耦向量的估计值，为γ的估计值，为h的估计值，将带入得到多项式，求解该多项式的k个根得到x1,x2,...,xk，则doa估计为，则doa估计为其中，∠(
·
)取复数的幅角主值；根据得到估计的互耦矩阵
[0049]
下面结合仿真实例对本发明方法的效果进行验证。
[0050]
仿真实例1：设置均匀线阵阵元数为m＝8，d＝λ/2，阵元互耦度p＝2，互耦参数为c1＝0.5+0.4j，c2＝0.1-0.03j，p
′
＝p。假设k＝2个入射信号，doa分别为θ1＝89
°
+δθ1，θ2＝105
°
+δθ2，δθ1与δθ2在[-0.5
°
,0.5
°
]内均匀分布。设置快拍数为50个，信噪比从-10db扫描至10db。比较背景技术中文献2的方法与本发明方法对波达方向和互耦系数的估计性能差异。为了提供性能参照，还提供了克拉美罗下限(crb)的结果。crb可通过文献：z.m.liu andy.y.zhou,“aunified framework and sparse bayesian perspective for direction-of-arrival estimation inthe presence ofarray imperfections,”ieee transactions on signal processing,vol.61,no.15,pp.3786
–
3798,2013推得。仿真结果如图2和图3所示。从图中可以看出，对于波达方向和互耦的估计，本发明方法在在低信噪比下的估计精度高于对比算法。
[0051]
仿真实例2：设置信噪比为15db，其余仿真条件与前述相同。比较各方法算法执行时间差异。仿真结果如图4所示。从图中可以看出，本发明方法的算法执行速度快于对比方法。单次执行时间小于1.2毫秒。
[0052]
仿真实例3：设置信噪比为-5db，其余仿真条件与前述相同。比较各方法的滤波器的零点位置。仿真结果如图5所示。从图中可以看出，只有本发明方法中的滤波器零点位于单位圆上，因此更符合真实阵列信号处理模型，可获得更高的参数估计精度。
[0053]
以上结合附图对本发明的实施方式作了详细说明，但本发明不限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明原理和精神的情况下，对这些实施方式进行多种变化、修改、替换和变型，仍落入本发明的保护范围内。