随着科学技术、数字通信等的迅速进步, 高精度和低延迟的空间谱估计也得到了飞速的发展. DOA^[1]估计是空间谱估计中的重要分支之一, 其广泛应用在雷达^[2]、通信、声呐、地震^[3]、勘测以及生物医学等众多工程领域^[4], 吸引了广大的学者研究.

常用的DOA估计算法有非高分辨和高分辨算法两种. 非高分辨算法早期以相关干涉仪^[5–7]算法及常规波束形成法(Conventional Beam Forming, CBF)^[8]为代表, 后来Capon法和最大熵法得到了发展. 20世纪70年代后期, 高分辨算法以多重信号分类(Multiple Signal Classification, MUSIC)^[9]、旋转不变子空间(Estimation of Signal Parameters via Rotational Invariance Techniques ,ESPRIT)^[10,11]最为突出. 两种算法在一定程度上提高了测向的精度. 随后, 以最大似然算法(Maximum Likelihood, ML)^[12,13]和权重子空间拟合(Weighted Subspace Fitting, WSF)^[14]等最具代表性. 常用的天线阵列有均匀线阵^[15]和均匀圆阵两种. 根据上面的文献, CBF也称Bartlett波束形成法, 它是传统时域上傅里叶谱估计方法的一种延伸, 存在着瑞利限制, 只能通过扩大阵列的孔径来提高信号的测量精度, 但实验成本也会随之增加, 不利于在实际中的应用. 相关干涉仪算法不需要增加阵列孔径的大小就可以比CBF算法的精度高, 其实质^[16]是利用电磁波信号到达固定间距的天线阵中不同天线阵元之间的时间差所产生的相位关系来确定无线电信号的方位. 同时, 因高分辨算法在实际应用中不能满足高实时性要求, 所以在实际应用中也面临着巨大的挑战^[17], 虽然之后也有些改进, 但是计算复杂度仍是很高^[18–20]. 相比较于直线阵列, 圆阵可以实现信号的全方位测量, 测向结果不随信号的增益强弱而改变^[21]. 相关干涉仪算法及圆阵的优点, 决定了其在实际测向中是比较好的方案.

相关干涉仪算法大致可以分为两种, 一类是全局搜索拟合, 另一类是比值法, 即根据数据的功率比值求解信号的角度. 在工程应用中已经根据比值法研究出具体的仪器, 并已投入到实际的应用中^[22], 如德国的施瓦兹测向仪等. 但是该仪器的测向精度仍不是很高. 针对实际中的应用受限和精度低的问题, 广大学者也进行过无数的研究.

虽然相关干涉仪算法具有测向高实时性, 但是该算法依旧存在着解模糊问题^[23,24]. 文献[25]曾提出当基线长度与信号波长的比值大于0.5时, 相位的求解存在整周期的多值问题, 即存在相位模糊问题; 当基线长度与信号波长的比值小于0.5时, 相位和测向角是正余弦函数关系, 存在基线镜像对称模糊. 文献[26]改进是根据正余弦在最大值处的相反性进一步提高信号测量的相关性(本文称之为SCCIA), 以求解出精度更高的信号角度, 但是信号的测量精度仍不是很高.

针对上述存在的基线镜像对称、相位模糊以及精度低的问题, 本文提出了一种改进相关干涉仪算法(Improved Correlation Interferometer Algorithm, ICIA). 该算法先将接收到的数据进行处理, 然后根据象限分类法确定信号的所属象限, 这样成功的避免了基线镜像对称和相位模糊问题. 且提高了信号测量的精度; 除此之外本文还利用幅度比值法降低了计算复杂度, 提高了测向的实时性. 仿真实验证明了本文算法的可靠性与有效性.

1 空间谱估计数学模型和相关干涉仪算法 1.1 空间谱估计数学模型

假设有一个远场窄带信号源, 其传播速度为c, 中心频率为w₀, 信号的波长为λ, 从水平角度为θ的方向入射到一均匀平面圆阵上. 该均匀圆阵阵元数为M、半径为R, 则相邻阵元间角度为360°/M. 理想情况下, 假设该圆阵中的各阵元是各向同性的, 不存在信号之间的干扰、互耦影响等, 且信号在特定的信道系统中传输, 传输过程中加入的噪声是零均值高斯白噪声.

如图1, 将均匀圆阵中的M个阵元分成两个单元, 单元1包括阵元1一个阵元, 作为参考阵元; 单元2包括其他M–1个阵元, 作为普通阵元. 每次测量信号时选择普通阵元中的一个阵元和参考阵元进行数据测量, 则信号到达普通阵元与参考阵元之间的时间差为:

${\tau _i} = \frac{R}{c}\cos \left( {\frac{{360 * (i - 1)}}{M} - \theta } \right)$

(1)

式(1)中, ${\tau _i}$ 为普通阵元与参考阵元接收到的信号的时间差, $i = {\text{2}},\cdots,M$ .

则普通阵元与参考阵元接收到信号的相位差 ${\varphi _i}$ 为:

${\varphi _i} = 360 * \frac{R}{c}\cos \left( {\frac{{360 * (i - 1)}}{M} - \theta } \right)$

(2)

其中, 假设接收到的噪声是零均值高斯白噪声 ${N_i}(t),i = 2, \cdots ,M$ . 则同一时刻, 参考阵元接收到的数据模型为 ${X_i}(t)$ , 普通阵元接收到的数据模型为 ${Y_i}(t)$ , 则:

$\left\{\begin{gathered} {X_i}(t) = {A_i}\cos ({w_0}t - {\varphi _i}{\text{/2}}) + {N_i}(t) \\ {Y_i}(t) = {B_i}\cos ({w_0}t + {\varphi _i}{\text{/2}}) + {N_i}(t) \\ \end{gathered} \right.$

(3)

式(3)中, ${A_i}$ 为参考阵元接收到信号的幅值, ${A_i}\cos ({w_0}t + {\varphi _i}{\text{/2}})$ 为参考阵元接收到的信号部分, ${B_i}$ 为第 $i$ 个阵元接收到的信号的幅值, ${B_i}\cos ({w_0}t + {\varphi _i}{\text{/2}})$ 为第 $i$ 个阵元接收到的信号部分, ${N_i}(t)$ 为阵元接收到的噪声部分, 且该噪声均值为零, 方差为 ${\sigma ^2}$ . 则两个阵元接收到的数据为 ${{X}}{Y_i} = \left[ {{X_i},{Y_i}} \right]$ , 再根据接收到的数据 ${{X}}{Y_i}$ 进行处理求解信号的到达方位, 最后通过参考阵元与每一个非参考阵元接收数据得到 ${{X}}Y{\text{ = }}[{{X}}{Y_2}, \cdots ,{{X}}{Y_M}]$ . DOA估计问题即可表述为: 给定观测数据 ${{X}}Y$ , 反推出信号到达方向. 在本文中, $\hat \theta $ 表示 $\theta $ 的推导值.

1.2 相关干涉仪算法

相关干涉仪算法的原理是通过无线电信号相对于天线阵参考方向的方位与天线阵元间信号的相位分布关系进一步推知无线电信号的入射方位. 圆形阵列中所选阵元组成的基线长度与信号波长的比例关系会使测向中出现基线镜像对称和相位模糊问题. 以下是相关干涉仪算法的原理模型.

根据图1的模型, 现以均匀5阵元阵列为例(如图2). 首先将天线阵列分为单元1和单元2, 其中单元1包含一个天线阵元并将其作为参考阵元; 单元2包含阵列中的其他四个阵元, 作为普通阵元, 以顺时针编号, 参考阵元记作阵元1, 其他阵元分别为阵元2, 3, 4和5.

假设单元二中的阵元切换到第 $i$ 个阵元 $\left( {i{\text{ = }}2,3,4,5} \right)$ , 此时阵元1接收到的数据为:

${x_{ij}} = {A_{ij}}\cos ({w_0}t - {\varphi _{ij}}/2)$

(4)

其中, $i = 2,3,4,5;j = 1, \cdots ,NUM$ , $NUM$ 为抽样次数, ${A_{ij}}$ 为第 $j$ 次参考阵元接收到的幅值, ${\varphi _{ij}}$ 为第 $j$ 次阵元1与第 $i$ 个阵元的接收数据的相位差.

单元二中第 $i$ 个阵元接收到的数据为:

${y_{ij}} = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{ij}}/2)$

(5)

其中, ${B_{ij}}$ 为第 $j$ 次非参考阵元 $i$ 接收到的幅值.

如图3, 将单元二中阵元接收到的数据通过四步移相(0°/90°/180°/270°), 得到如下数据:

$\left\{\begin{gathered} y1 = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{ij}}/2 + 0^\circ ) \\ y2 = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{ij}}/2 + 90^\circ ) \\ y3 = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{ij}}/2 + 180^\circ ) \\ y4 = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{ij}}/2 + 270^\circ ) \\ \end{gathered} \right.$

(6)

其中, $i = 2,3,4,{\text{5}};j = {\text{1}}, \cdots ,NUM$ , $y1,y2,y3,y4$ 分别表示单元二中接收到的数据经过四步移相之后的结果.

然后将移相之后的数据结果与参考阵元接收到的数据进行矢量相加, 得到的数据为:

$\left\{\begin{gathered} A1 = {x_{ij}} + y1 \\ A2 = {x_{ij}} + y2 \\ A3 = {x_{ij}} + y3 \\ A4 = {x_{ij}} + y4 \\ \end{gathered} \right.$

(7)

再进行A/D转化得到四个对应的绝对幅值矢量(A1/A2/A3/A4), 最后将4个幅度值(A1/A2/A3/A4)的结果代入下式(8)中, 即可得到参考阵元与第 $i$ 个阵元间的相位差.

${\varphi _i} = - \arctan \left( {\frac{{|A2{|^2} - |A4{|^2}}}{{|A1{|^2} - |A3{|^2}}}} \right)$

(8)

根据公式(8)和下式(9)可以求解出信号的入射角 ${\theta _i}'$ :

${\theta _i}'{\text{ = }}\frac{{360 * (i - 1)}}{M} \pm {\rm{arc}}\cos \left( {\frac{{{\varphi _i}}}{{360}} \times \frac{\lambda }{R}} \right)$

(9)

单元二中的每个阵元都要进行上述的同样求解过程得到相应的测向角, 从而求出 $\hat {\hat \theta} = [{\theta _2}',\theta _3',\theta _4',\theta _5']$ , 求解的 $\hat {\hat \theta} $ 即为得到的信号角度, 再经过 $NUM$ 次求解以减小偶然性误差, 得到 $\hat {\hat \theta} $ 的平均值 $\hat \theta $ 为最后的信号求解角度.

2 两类误差模型和改进相关干涉仪算法 2.1 两类误差模型

在波达方位估计中, 阵列长短基线导致了相关干涉仪算法在测向中存在基线镜像对称和相位模糊问题. 其一, 当基线长度与信号波长的比值大于0.5时, 由公式(8)、(9)可知, 当求解出 ${\varphi _i}$ 时即可得到信号的方位角 $\hat \theta $ , 但是, 由于 ${\varphi _i}$ 的周期性, 将会引起信号方位角 $\hat \theta $ 的误差, 即存在相位模糊问题. 其二, 当基线长度与信号波长的比值小于0.5时, 由公式(9)知, 相位和测向角是三角函数关系, 信号的方位角无论是正负, 求解出的信号角度都将在(0~180°)范围内, 故当信号为负时, 信号的求解便存在误差, 即存在基线镜像对称模糊问题.

两类问题可以简单的描述如表1.

传统的算法是通过定义归一化的相关系数 ${\rho _i}$ , 并通过求解其大小确定信号的来波方向. 其中归一化的相关系数 ${\rho _i}$ 为:

${\rho _i} = \sum {\cos ({\varphi _i} - {\varphi _0})} $

(10)

式(10)中, ${\varphi _0}$ 为理论上求解出的相位差, ${\varphi _i}$ 为实际求解出的相位差.

传统做法通过一维搜索过程虽然可以解决基线镜像对称和相位模糊两类问题. 但是信号的角度测量误差比较大, 不利于在实际系统中的应用.

之后文献[26]对其进行了改进, 即根据正余弦值的相反性通过一维搜索过程进行改进, 具体代价函数为:

${\rho _i} = \sum {\cos ({\varphi _i} - {\varphi _0}) - } {\text{|sin(}}{\varphi _i} - {\varphi _0}){\text{|}}$

(11)

但是该做法仍不能满足测向的精度要求. 为了解决两类问题, 并提高测向的精度, 本文提出了一种新的改进相关干涉仪算法.

2.2 改进相关干涉仪算法

象限分类法顾名思义就是根据角度的不同将其分至四个象限中, 如150°则为第二象限. 本文处理的方式是将接收的数据先化简至(0°~360°), 然后再分至四个象限中, 这样可以成功的避免基线镜像对称问题和相位模糊问题.

改进相关干涉仪算法具体方法是: 将普通阵元与参考阵元接收到的信号相位差与360°作比取余, 将余数采用象限分类法分至不同的象限, 再将数据进行求解得到信号的绝对幅值, 利用比值法求解信号的入射角, 再加以还原信号的绝对入射角, 就可以得到信号的真实估计角度, 从而避免在测向中存在的基线镜像对称和相位模糊问题, 实现信号的准确估计. 具体步骤如下(以上述均匀五阵元为例):

改进相关干涉仪算法如下:

算法1. 改进相关干涉仪算法.

1) 将天线阵列分成两个单元, 单元1包含一个阵元作为参考阵元, 单元2包括其他四个阵元作为普通阵元, 如图2.

2) 任意选择一组基线: 选择单元1中的参考阵元和单元2中的第 $i$ 个阵元.

针对单元1, 在特定的时间记录接收到的数据, 如公式(4); 在同一时间, 针对单元2中的阵元, 利用高速射频开关(一般保证信号的奈奎斯特抽样准则以便于信号的还原)在单元二中的阵元之间进行高速切换, 并记录第 $i$ 个阵元接收到的数据, 如公式(5).

3) 将单元二中接收到的数据通过四步移相(0°/90°/180°/270°), 得到信号数据如公式(6).

4) 将得到的 $\scriptstyle y1,y2,y3,y4$ 中的相位与360°作比取余数 $\scriptstyle {\varphi _{0ij}}$ , 利用象限分类法将 $\scriptstyle {\varphi _{0ij}}$ 分到四个象限中, 这样避免了基线镜像对称和相位模糊问题, 并将 $\scriptstyle {\varphi _{0ij}}$ 带回公式(6)得到 $\scriptstyle y{\text{1}}1,y{\text{2}}2,y3{\text{3}},y4{\text{4}}$ .

$\scriptstyle \left\{\begin{gathered}\scriptstyle y{\text{1}}1 = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{{\text{0}}ij}}/2 + 0^\circ ) \\\scriptstyle y2{\text{2}} = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{{\text{0}}ij}}/2 + 90^\circ ) \\\scriptstyle y3{\text{3}} = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{{\text{0}}ij}}/2 + 180^\circ ) \\\scriptstyle y4{\text{4}} = {B_{ij}}\cos ({w_0}t{\text{ + }}{\varphi _{{\text{0}}ij}}/2 + 270^\circ ) \\ \end{gathered} \right.$

(12)

其中, $\scriptstyle {\varphi _{0ij}}{\text{ = }}{\varphi _{ij}} \pm 360 * k$ , $\scriptstyle k$ 为整数;

5) 将得到的 $\scriptstyle y11, y22, y33, y44$ 与参考阵元接收到的数据 $\scriptstyle {x_{ij}}$ 分别进行图3中的矢量相加得到 $\scriptstyle A1, A2, A3, A4$ , 即:

$\left\{\begin{gathered}\scriptstyle A1 = {x_{ij}} + y11 \\\scriptstyle A2 = {x_{ij}} + y22 \\\scriptstyle A3 = {x_{ij}} + y33 \\\scriptstyle A4 = {x_{ij}} + y44 \\ \end{gathered}\right. $

(13)

6)将上式得到的 $\scriptstyle A1,A2,A3,A4$ 分别取模平方代入公式(8)、(9)求解 $\scriptstyle\hat \theta $ ;

7)如达到迭代次数 $\scriptstyle NUM$ ( $\scriptstyle NUM$ 的取值只要能排除实验的偶然性即可), 则退出停止, 否则返回步骤2).

通过上述分析, 与传统相关干涉仪算法相比, 本文提出的改进算法避开了一维全局搜索的过程可以降低在测向过程中的计算复杂度, 同时采用象限分类法避免了在求解过程中由于正余弦值的对称性而引起的基线镜像对称和相位模糊问题; 同时本文采用的是相关干涉仪算法中的比值法求解信号的角度, 提高了信号测量的精确度和实时性.

3 实验分析

在试验中, 假设阵列为均匀五阵元平面阵(即相邻阵元之间的角度差为72°), 则阵元间的相位差可以简化为:

${\varphi _i} = 360 * \frac{R}{c}\cos \left( {\frac{{360(i - 1)}}{{\text{5}}} - \theta } \right)$

(14)

在实验中信噪比(Signal-to-Noise Ratio, SNR)的定义为:

${\rm{SNR}} = 10{\log _{10}}\frac{{E\left[ {|x(t){|^2}} \right]}}{{{\sigma ^2}}}$

(15)

其中, $x(t)$ 天线阵列接收到的信号, ${\sigma ^2}$ 是噪声的功率. 同时定义均方根误差(Root-Mean-Square-Error, RMSE)为:

${\rm{RMSE}}= \sqrt {\frac{1}{{NUM}}\sum\limits_{m = 1}^{NUM} {|{{\hat {\hat \theta} }_m} - {\theta _m}{|^2}} } $

(16)

其中, ${\hat {\hat \theta} _m}$ 是 ${\theta _m}$ 在第m次实验所得出的推导值, 采用蒙特卡洛的实验方法, 每一种情况进行 $NUM$ 次实验.

实验一: 验证本文提出的改进相关干涉仪算法

在解决“基线镜像对称和相位模糊问题”上的可行性, 本文从以下几个方面进行仿真分析. 仿真一、仿真二和仿真三的实验条件如下: 信号速度为 $c = 3 \times {10^8}$ m/s, 信号频率为 $f = {\text{3}}0$ MHz, 信号的波长为 $\lambda {\text{ = }}c/f$ , 采用均匀5阵元阵列, Monte-Carlo实验为 $NUM = {\text{100}}$ 次, 信噪比 $SNR = 5\sim 20$ dB.

仿真1. 短基线测向性能比较(基线长度与信号波长比小于0.5即为短基线)

根据图4可以看出, 在该条件下, ICIA在性能上远优于CIA和SCCIA的性能, 在信噪比为5~20 dB条件下, ICIA均方根误差均小于0.5°, 而CIA和SCCIA均方根误差较大; 在信噪比大于15 dB时, 三种算法均方根误差均小于1°, 本文的改进算法性能更好一些, 有利于在实际工程中的应用.

仿真2. 长基线测向性能比较(基线长度与信号波长比大于0.5即为长基线)

根据图5可以看出, 在该条件下, ICIA在性能上远优于CIA和SCCIA的性能, 在信噪比为5~20 dB条件下, ICIA的均方根误差均小于0.5°, 而CIA和SCCIA均方根误差较大; 在信噪比大于15 dB时, 三种算法均方根误差均不大于0.5°, 本文的改进算法性能更好一些, 更利于在实际工程中的应用, 本图与图4相比较可知, 该算法的天线基线较大, 产生的误差较小, 在工程实践中对天线阵列的设计具有参考意义.

仿真3. 本文提出的算法在不同基线长度与信号波长比的条件下性能的比较

根据图6可以看出, 改进相关干涉仪算法在不同基线长度与信号波长比条件下的均方根误差均小于0.4°. 在基线长度与信号波长比相同的条件下, 均方根误差随着信噪比的增加而逐渐降低. 在信噪比相同时, 基线长度与信号波长比为2时, 均方根误差最小, 基线长度与信号波长比为0.3的均方根误差也限定在0.4°以内, 可认为是较为成功的消除了基线镜像对称和相位模糊问题. 因此, 本文改进的算法具有很高的参考价值, 有利于实际工程中的应用.

实验二: 验证本文提出算法的计算复杂度如表2. 实验条件: 信号速度为 $c = 3 \times {10^8}$ m/s, 信号频率为 $f = {\text{3}}0$ MHz, 信号的波长为 $\lambda {\text{ = }}c/f$ , 采用均匀5阵元阵列, Monte-Carlo实验为 $NUM = {\text{100}}$ 次, 信噪比 $SNR = {\text{1}}$ dB.

从表2可以看出, CIA和SCCIA在精度为0.1°时迭代的次数为3600, 运行时间比较长, 本文改进的算法使用接收信号的幅值进行计算, 降低了计算复杂度, 减小了运行的时间, 提高了测向的实时性. 在实际的应用中具有很好的指导意义.

4 结论与展望

本文针对相关干涉仪算法求解时遇到的基线镜像对称和相位模糊问题提出了一种改进相关干涉仪算法. 该算法根据信号到达阵列信号的相位差与信号的角度具有一一对应的关系, 先对相位求整取余, 再将余数通过象限分类, 确定信号的所处象限, 接下来根据传统相关干涉仪算法求解信号的初始角度, 最后通过逆运算求解出信号的最终角度. 在成功的解决基线镜像对称和相位模糊问题的基础上提高了信号的估计精度, 降低了计算复杂度, 提高了测向的实时性, 在波达方位估计中具有重要意义.