相干信源DOA 估计的最大似然算法应用及改进

2022-10-20陈佳

电子技术应用 2022年9期

陈佳

(中国船舶集团有限公司第七二三研究所，江苏扬州 225000)

0 引言

雷达信号的波达方向（Direction of Arrival，DOA）估计在阵列信号处理领域中一直是一个非常重要的研究内容，主要通过提取出空间中按照一定样式排布的天线阵列接收回波信号的特征参数来估计目标的方位信息，因此国内外学者提出了各种各样的DOA 估计新算法[1-3]。其中最大似然估计法[4]很早就被提出，但由于其运算量较大而没有得到广泛应用。本文讨论了一种有效的算法，通过交替投影算法[5-6]把最大似然估计中多维非线性问题转化为一维问题来解决，并在搜索时采用斐波那契数列法提高搜索效率，使其运算复杂性降低同时有效减少了计算量，满足工程应用的实时性要求。

1 阵列信号处理的数学模型[7-9]

假定天线阵列由m 个阵元组成，后面接有M 个接收通道，d(d ＜M)个窄带远场、独立的目标信号，其方向角分别为θ1,θ2，…，θd，各个阵元噪声为：n1(t)，n2(t)，…，nm(t)，假定噪声之间相互独立、功率相同，且是空间平稳的高斯分布白噪声，第i 个阵元的输出为：

式中，sk(t)为第k 个目标回波的信号；τi(θk)为第k 个回波信号到达第i 个阵元时相对于参考阵元的时延。有：

上述各个矢量表示分别为：

得出阵列输出矢量形式为：

写成矩阵形式为：

式中，x(t)为阵列的M×1 维快拍数据矢量，s(t) 为目标信号的d×1 维数据矢量，A(θ)为M×d维阵列流型（导向矢量）矩阵，为信号向量张成的子空间（信号子空间），n(t)为阵列M×1 维噪声数据矢量。此式为DOA 估计的基本数学模型表达式。DOA 估计任务就是利用输出信号x (t) 估计出信源个数及所在方位θi，i=1,2,…,d。

2 最大似然估计算法[10-15]

令N 次采样数据为{x1,x2,…,xN}，信号源数据为{s1,s2,…,sN}，噪声数据为{n1,n2,…,nN}。如果噪声信号{ni}是均值为零，方差为σ2的各态历经复高斯过程，得出观测矢量L 次快拍联合概率密度函数为：

式中det[]表示矩阵的行列式，两边同时取负对数为：

求解上式最大似然估计问题，固定θ 和si，（忽略常数项）得出σ2的最大似然估计：

将上式代入（11）可得到θ 和si的最大似然估计：

根据对数函数的单调性原则，可将上式等效为：

固定θ,求上式极小得到si的估计值为：

定义A 的投影阵为：

将（15）代入（16）中可得到确定性最大似然的准则，即：

可得最大似然准则以迹的形式表示为：

综上所述，就可以得到关于θ 的最大似然估计器，上式一般用修正的牛顿多维优化方法求解。由于在计算最大似然估计中需要多维非线性搜索，并且在每次的迭代搜索都要计算A 的投影矩阵PA，能够看出这种算法的计算量相当大，因此提出了交替投影迭代（Alternating Projection iterative algorithm，AP）算法应用于最大似然估计算法研究中用以减少计算量。本文以两个信源的情况描述AP 算法，很容易就能够扩展到多个信源的情况。

AP 算法具体流程如下：

（1）估计第一个信源的情形，为：

（2）估计第二个信源，假定第一个信源的方位为θ1

（3）进行迭代得：

（4）继续迭代得：

（5）重复（3）和（4）直到收敛。

具体流程如图1 所示。

图1 两维交替投影过程

采用交替投影法将多维搜索转化到一维搜索上，但是随着测角精度提升，搜索次数也随之增大，再此基础上采用一维斐波那契数列法搜索极值点。

斐波那契数列法，其基本思想是通过试探点函数值比较，使包含极大点的搜索区间不断缩小。该方法仅需要计算函数值，适用范围广，使用方便。

斐波那契数列法是建立在区间消去法原理基础上的试探方法，即在搜索区间[θmin,θmax]内适当插入两点θ1，θ2，并计算其估计值。

θ1，θ2将整个搜索区间分为三段，应用函数的单峰性质，通过函数值大小的比较，删去其中一段，使搜索区间得以缩小。然后再在保留下来的区间上作同样的处理，如此迭代下去，搜索区间无限缩小，从而得到极大点的数值近似解。

具体步骤如下：

（1）选定初始区间[θmin,θmax]，及ε＞0，利用式（b1-a1）求出计算函数值的次数n。并设eps=1×106。接着由下式：

计算试探点p1和q1。令k=1。

（2）如果f(pk)＜f(qk)，转步骤（3）；否则转步骤（4）。

3 仿真结果

本文以8 阵元为例，进行最大似然算法测向仿真。均匀线阵，阵元数为M=8，阵元之间间隔为λ/2，信号频率为f=24.3 GHz (λ=12.3 mm)，快拍点数为1 个，信源数为d=2，信噪比为5 dB，设噪声服从高斯分布的白噪声。信源入射角为-5°和5°，搜索角范围为（-18°，18°），以0.01°为一个步进。采用AP 方法估计来波信号方向角的仿真结果如图2 所示。

图2 列出4 次迭代效果图，由图2 可知第一次估计目标1 的初始角度为-0.93°，然后根据-0.93°估计出目标2 的初始角度为1.75°，第一次迭代根据目标2 的1.75°估计出目标1 角度为-4.56°，再根据-4.56°估计出目标2 为5.08°，直到估计目标1 为-4.98°，目标2 为5.01°满足收敛条件结束迭代，最终测出角度为-4.98°和5.01°。

图2 部分AP_DML 算法仿真结果1

其他条件与上述仿真一致，当信源入射角为-1.5°和8°时，仿真结果如图3 所示。

从图3 可以得出估计第二种仿真模型的计算的测角结果目标1 为-1.52°，目标2 为7.91°，满足收敛条件结束迭代。如果需要更精确的角度值，可将收敛条件提高，那么迭代次数也会增加，计算角度会更加精确。

图3 部分AP_DML 算法仿真结果2

从上文中以第一种仿真结果为例，能轻易估算出搜索一次角度为3 600 次，本次迭代次数为6 次，那么总共需要14 次角度搜索，总搜索次数为50 400，如果迭代次数增加，搜索总次数也随之增加，耗时量十分巨大，因此在AP 方法的基础上采用斐波那契数列法极值点。斐波那契数列法求极值点的前提是整个搜索区间内只有一个极值点，即一个波峰或者波谷。以估计第一个目标的角度为例（搜索其他估计角度的极值点类似），根据上述仿真发现，当搜索角度区间足够大时，会出现多个波峰情况，因此针对这种情况本文限定了搜索区间，除此之外可以看出除初始估计第一个信源角度外，其他估计时会出现一个凹口，因此也不能完全满足斐波那契数列法求极值点要求，但是根据算法得出凹口位置为本次迭代角度值，即初始估计第二个信源角度时凹口为初始估计第一个信源角度的值2.25°，因此在用斐波那契数列法求极值时可跳过本次迭代的角度值。具体搜索过程如图4 所示（以第一个仿真模型初始估计第一个信源角度为例）。

图4 中黑色点为整个斐波那契数列法的参考点，统计搜索出峰值点位置的总计算周期次数为16 次，与原先3 600 次相比大大减少周期计算次数，且最终结果与之前一致，显著提高了整个角度估计的效率。根据仿真验证得出当满足斐波那契数列法求极值的前提条件时，该算法用于其他角度搜索的最大次数一般不超过20次，所以能够显著提高峰值搜索效率。