一种基于平方和与半定规划的最大似然波达方向估计方法与流程

本发明属于雷达技术领域，涉及阵列信号处理领域，尤其涉及一种基于平方和与半定规划的最大似然波达方向估计方法。

技术背景

远场窄带信号源的波达方向估计(direction of arrival，DOA)是雷达、声纳、卫星通信和无线通信等领域的一个热点问题。DOA估计的目标是根据阵列上的观测样本，估计空间中多个信号源的方位。近几十年来，研究人员提出了多种类型的方向估计方法，包括基于最大似然(Maximum Likelihood，ML)准则的方法和基于子空间的方法。其中，基于最大似然准则的方法能实现最优的估计性能。然而该类方法通常需要求解一个非线性非凸的多维优化问题，计算复杂极高。交替投影方法(Alternating Projection，AP)将多维优化问题转化为一系列一维优化子问题来近似求解。然而这些子问题仍然是非凸优化问题，难以求解，因此通常通过穷举法搜索最优值，计算复杂仍然很高。子空间类的方法，如多维信号分类算法(Multiple SIgnal Classification，MUSIC)、基于旋转不变性的信号参数估计技术(Estimation of Signal Parameters via Rotational Invariance Technique，ESPRIT)和MUSIC求根(Root-MUSIC)算法，能以较低的复杂度进行精确的方向估计。但是这类算法在估计高度相关的信号源方向时偏差较大。另外，通过利用均匀线阵(uniform linear array，ULA)阵列流型的特殊结构，迭代二次最大似然算法(iterative quadratic maximum likelihood，IQML)及其改进的方向估计方法MODE能够以闭式解迭代的方式逼近最大似然的最优性能，且其计算复杂度较低。最近，研究人员提出了基于压缩感知的DOA估计方法，如稀疏参数估计方法(Sparse and Parameter Algorithm)。压缩感知方法有很多特性，如对于不受信号相关性的影响，在一次快拍的情况下实现方向估计，不需要预知信号源数目。但是这类方法空间分辨率较低且容易产生错误的方向估计。

技术实现要素：

针对现有DOA估计方法的不足，如计算量大，不能直接估计相关信号源方位，分辨率不够高等，本发明提出了一种的基于平方和与半定规划的最大似然DOA估计方法。该方法基于最大似然准则，对传统的交替投影方法进行了改进，不需要通过穷举法搜索最优值，稳健性和计算效率都得到了提高。

为解决上述技术问题，本发明采用以下技术方案，包括以下步骤：

步骤1，根据阵列天线的接收信号模型和最大似然准则建立方向估计优化问题，利用坐标轮换法和交替投影原理将其转化为一系列优化子问题；

步骤2，将子问题中的方向角变量代换为实数变量t，并将子问题构造成分式多项式优化问题；

步骤3，利用平方和特性将多项式优化问题转化为半定规划问题进行求解；

步骤4，用所得最优解构建关于t的一元高次方程组并求解；解值t所对应的原方向变量即为信号源波达方向估计；当迭代收敛时，得到信号源波达方向。

以上技术方案的具体步骤如下，其中步骤2到步骤4为方案的主要特点和改进：

(1)步骤1具体包括以下子步骤：

1a)设定天线阵为均匀线阵，阵元数目为N，工作中心波长为λ，阵元间距为d；空间中有M个随机分布的远场窄带平稳信号s_m(k)分别从方向θ_m照射该天线阵列，θ＝[θ₁,θ₂,...,θ_M]；设阵列噪声为加性高斯白噪声，接收信号可表示为：

x(k)＝A(θ)s(k)+n(k),k＝1,2,…,N_t (1)

其中x(k)为N×1维阵列接收信号，s(k)＝[s₁(k),s₂(k),…,s_M(k)]^T为M×1维远场窄带信号矢量，M为发射端信源个数，k为时刻，n(k)为加性噪声向量，A(θ)为N×M维的阵列流形矩阵，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_M)]，T表示转置；

1b)根据以上信号模型，得到方向估计的似然函数

以和s(k)为优化变量，最大化似然函数L，得到它们关于θ的解析解，并代回到似然函数中，将最大似然估计等价转变为以下优化问题

其中，为接收信号自相关矩阵，P_A(θ)＝A(θ)(A^H(θ)A(θ))^-1A^H(θ)为阵列流形A(θ)的投影矩阵。

1c)根据坐标轮换法，将问题改造为一系列对θ_m,m＝1,2,…,M进行迭代估计的子问题：定义阈值ε>0，其中下标m＝1,2,…,M表示DOA的编号，上标k＝1,2,…,K表示第k次迭代，K为最大迭代次数，表示第k次迭代中第m个方向的估计值；在第k轮迭代中对θ_m进行估计时，固定其他M-1个DOA等于根据交替投影原理，定义如下矩阵和向量并将它们代入式中，得到k轮迭代中对θ_m进行更新的优化子问题

(2)步骤2具体包括以下子步骤：

2a)定义将问题中目标函数的母子和分母分别表示为

2b)定义变量代换并代入a(θ)中，可将a(θ)的第(k+1)个元素表示为如下形式：

再定义变量代换t＝tan(v)，并将三角变换和代入式得到

其中，而h_k^r(t)与h_kⁱ(t)分别为多项式(1-t²+2jt)²的实部和虚部；根据以上变量代换，当设定时，对于t∈R，双射随t单调递增，且θ的值域为若设定则对于t∈R，θ的值域变为

2c)将式分别代入式和式中得到

其中，表示R_i的第k行l列的元素，和分别表示m_i,k的实部与虚部；再将式代入以下定义得到

2d)将f₂(t)和f₁(t)代入多项式优化问题中，将其转化为最大化以下分式多项式问题：

其中R代表实数域；

(3)步骤3具体包括以下子步骤：

3a)将最大化分式多项式问题其等价转化为求解以下多项式的最小上界p：

根据定义f₁(t)>0，因此将上式转化为

3b)问题中的约束条件等价于pf₁(t)-f₂(t)可以表示成平方和的形式，即存在N维半正定矩阵Z，使得下面的等式恒成立

其中，t＝[1,t,...,t^N-1]^T；因此，可将优化问题表示为：

其中，H^(N,k)是N维汉克尔矩阵，并满足

3c)利用凸规划包解得p和Z的最优解p^*和Z^*；

(4)步骤4具体包括以下子步骤：

4a)根据步骤4所求得的最优解Z^*，建立以下一元高次方程组：

Z^*t＝0 (17)

4b)求解该方程组：定义Z^*的零空间为N(Z^*)，且N(Z^*)的秩为r_n，其中r_n＝M；用高斯消元法对式进行消元得到N-r_n个r_n阶方程，分别求这N-r_n个方程的根；其中每个方程的前M-1个根对应中的方向，第M个根可通过多项式因式分解后系数恒等关系求解得到，将所求第M个根表示为t^*，通过以下关系得到所估计的DOA为

4c)当均更新一次以后，令检验以及k>K是否成立；若两者都不成立，令迭代次数k＝k+1，重复步骤2到步骤5；否则迭代终止迭代，信号源的方向即为

有益效果

与现有技术相比，本发明具有突出的实质性特点和显著的进步。本发明与现有方法相比，具有以下优点：

1.传统的方向高分辨方法，例如MUSIC法、ESPRIT法、Capon法等，只能处理非相关信号，对相关信号需要首先进行去相关，限制了这些算法的应用。而本发明方法是基于最大似然表示的DOA估计方法，对信号的相关性不敏感，如图2(a)和图2(b)，因而能直接对任意相关性信号源的DOA进行有效估计，因此应用较广泛。

2.原有的交替投影算法虽然基于最大似然准则，但是在迭代过程中估计每一维的方向时需要通过穷举法来搜索方向，因此当方向估计精度增加时，搜索精度需要增加，从而计算复杂度也会大大增加；本发明将该问题转换为凸优化问题进行求解，能够以较为恒定的复杂度求得全局最优解，大大提高了DOA估计的精度和分辨率，如图3(a)和图3(b)所示。

本发明更进一步涉及阵列信号处理技术领域中的一种基于平方和与半定规划的最大似然波达方向估计方法，可用于解决阵列接收到的信号具有相关性且噪声功率未知情况下的多目标波达方向估计问题。

附图说明

下面结合附图和具体实施方式对本发明做进一步说明。

图1为本发明方法流程图；

图2为仿真实验一中本发明方法和IQML、MODE方法在不同信噪比条件下的性能图：图2(a)是均方根误差图，图2(b)是检测概率图；

图3为仿真实验二中本发明方法和IQML、MODE和SPA方法在信号源方位角间距变化时的性能图：图3(a)为均方根误差图，图3(b)是检测概率。

具体实施方式

以下结合附图具体说明。

参照图1，说明本发明一种基于平方和与半定规划的最大似然波达方向估计方法，其具体实施步骤如下：

步骤1，根据阵列天线的接收信号模型和最大似然准则建立方向估计优化问题，利用坐标轮换法和交替投影原理将其转化为一系列优化子问题；

1a)设定天线阵为均匀线阵，阵元数目为N，工作中心波长为λ，阵元间距为d；空间中有M个随机分布的远场窄带平稳信号s_m(k)分别从方向θ_m照射该天线阵列，θ＝[θ₁,θ₂,...,θ_M]；设阵列噪声为加性高斯白噪声，接收信号可表示为：

x(k)＝A(θ)s(k)+n(k),k＝1,2,…,N_t (1)

其中x(k)为N×1维阵列接收信号，s(k)＝[s₁(k),s₂(k),…,s_M(k)]^T为M×1维远场窄带信号矢量，M为发射端信源个数，k为时刻，n(k)为加性噪声向量，A(θ)为N×M维的阵列流形矩阵，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_M)]，T表示转置；

1b)根据以上信号模型，得到方向估计的似然函数

以和s(k)为优化变量，最大化似然函数L，得到它们关于θ的解析解，并代回到似然函数中，将最大似然估计等价转变为以下优化问题

其中，为接收信号自相关矩阵，P_A(θ)＝A(θ)(A^H(θ)A(θ))^-1A^H(θ)为阵列流形A(θ)的投影矩阵。

1c)将方向估计优化问题通过坐标轮换法改造为一系列的一维方向估计子问题，并利用交替投影原理对子问题的代价函数进行化简；

根据坐标轮换法，将问题改造为一系列对θ_m,m＝1,2,…,M进行迭代估计的子问题：定义阈值ε>0，其中下标m＝1,2,…,M表示DOA的编号，上标k＝1,2,…,K表示第k次迭代，K为最大迭代次数，表示第k次迭代中第m个方向的估计值；在第k轮迭代中对θ_m进行估计时，固定其他M-1个DOA等于根据交替投影原理，定义如下矩阵和向量并将它们代入式中，得到k轮迭代中对θ_m进行更新的优化子问题

本发明使用最大化似然函数L的方法来进行方向估计，对信号源的相关性不敏感，因此适用范围广。似然函数L中的噪声功率和信号波形s(k)是未知的，求其关于θ的解析解，并代回到似然函数L中，消除方向角以外的未知变量，避免对冗余参数的估计，从而提高了本发明方法的稳健性，但同时也使目标函数形式变得更为复杂，这是最大似然方向估计的难点所在；因此利用坐标轮换法和交替投影原理，将多维方向估计问题转化为一系列的一维方向估计子问题，以避免求解多维问题，降低算法复杂度；但所得子问题仍然不能高效的精确求解，目前的方法是通过穷举法对定义域进行搜索，但是搜索精度越高复杂度越大。

步骤2，将子问题中的方向角变量代换为实数变量t，并将子问题构造成分式多项式优化问题；

本发明首先将变量θ代换为实数变量t，这是因为：一、通过变量代换将未知变量从指数形式变为一元高次多项式形式，便于进行分析；二、现有的数学理论能够处理一元高次多项式优化问题。

2a)定义将问题中目标函数的母子和分母分别表示为

2b)定义变量代换并代入a(θ)中，可将a(θ)的第(k+1)个元素表示为如下形式：

再定义变量代换t＝tan(v)，并将三角变换和代入式得到

其中，而h_k^r(t)与h_kⁱ(t)分别为多项式(1-t²+2jt)²的实部和虚部；根据以上变量代换，当设定时，对于t∈R，双射随t单调递增，且θ的值域为若设定则对于t∈R，θ的值域变为

2c)将式分别代入式和式中得到

其中，表示R_i的第k行l列的元素，和分别表示m_i,k的实部与虚部；再将式代入以下多项式定义得到

2d)将f₂(t)和f₁(t)代入多项式优化问题中，将其转化为最大化以下分式多项式问题：

其中R代表实数域；

步骤3，利用平方和特性将多项式优化问题转化为半定规划问题进行求解；

3a)将最大化分式多项式问题其等价转化为求解以下多项式的最小上界p：

根据定义f₁(t)>0，因此将上式转化为

3b)问题中的约束条件等价于pf₁(t)-f₂(t)可以表示成平方和的形式，即存在N维半正定矩阵Z，使得下面的等式恒成立

其中，t＝[1,t,...,t^N-1]^T；因此，可将优化问题表示为：

其中，H^(N,k)是N维汉克尔矩阵，并满足

3c)利用凸规划包解得p和Z的最优解p^*和Z^*；

在本发明中，凸规划包为本领域技术人员所公知的软件处理程序，例如SeDuMi和CVX等。

步骤4，用所得最优解构建关于t的一元高次方程组并求解；解值t所对应的原方向变量即为信号源波达方向估计；当迭代收敛时，得到信号源波达方向；

4a)根据步骤4所求得的最优解Z^*，建立以下一元高次方程组：

Z^*t＝0 (17)

4b)求解该方程组：定义Z^*的零空间为N(Z^*)，且N(Z^*)的秩为r_n，其中r_n＝M；用高斯消元法对式进行消元得到N-r_n个r_n阶方程，分别求这N-r_n个方程的根；其中每个方程的前M-1个根对应中的方向，第M个根可通过多项式因式分解后系数恒等关系求解得到，将所求第M个根表示为t^*，通过以下关系所估计的DOA为

4c)当均更新一次以后，令检验以及k>K是否成立；若两者都不成立，令迭代次数k＝k+1，重复步骤2到步骤5；否则迭代终止迭代，信号源的方向即为

本发明通过步骤3和步骤4两步求解最大化分式多项式问题，首先利用变量代换和平方和特性将一维方向估计问题转化为半定规划问题，求得问题目标函数的最优值p^*，然后利用该信息建立方程组求解与最优值p^*相对应的最优解t^*，从而在理论上保证了子问题求得全局最优解。相比于传统穷举法，本发明的复杂度恒定且求解精度和稳定性大大增加。

下面结合仿真实验对本发明的效果做进一步说明。

(1)实验条件：

为了进一步说明本发明的基于平方和与半定规划的最大似然波达方向估计方法较传统DOA估计方法的优越性，做如下两个仿真实验。

系统模型：采用阵元数为N的均匀线阵，阵元间距为半波长，平稳快拍观测数N_t＝100。本发明采用常用的针对DOA估计方法的性能评测指标，即方向估计的均方根误差(root mean square error，RMSE)和正确检测出信号的概率来评价不同方法的性能，仿真实验中的每一幅图中的每一个点均由1000次独立实验得到。

(2)实验结果分析

实验一：假设有2个远场窄带平稳相干信号源，其信号的相关系数为1，到达角为：θ₁＝△u和θ₂＝-△u，其中△u＝0.2165/2×BW_NN，BW_NN＝2arcsin(2/N)，BW_NN为阵列的主瓣宽度。设阵列阵元数N＝12，噪声为零均值的复高斯白噪声。

图2给出了本发明方法和IQML、MODE方法在不同信噪比条件下均方根误差和检测概率。如图2(a)所示，横坐标表示信噪比，纵坐标表示均方根误差，“克拉美罗界”为无偏估计方法的均方根误差的下界。图2(a)表明，三种方法的估计性能对信号的相关性都不敏感，本发明方法比其他方法更快的逼近克拉美罗界，因此信噪比较低时性能更加稳健。如图2(b)所示，横坐标表示信噪比，纵坐标表示检测概率。图2(b)表明，在信噪比较低时，本发明方法比其他方法更稳健。

实验二：假设有2个远场窄带平稳独立信号源，两个信号源的中心为0，间距从0.02BW_NN变化到0.2BW_NN。设阵列阵元数N＝10，噪声为零均值的复高斯白噪声，信噪比为10dB。

图3给出了本发明方法和IQML、MODE和SPA方法在信号源方位角间距变化时的均方根误差和检测概率。如图3(a)所示，横坐标表示归一化的信号源方位角间距，纵坐标表示均方根误差。图3(a)表明，随着信号源间距的增加，本发明方法比其他方法更快的逼近克拉美罗界，因此有更高的空间分辨率。如图3(b)所示，横坐标表示归一化的信号源方位角间距，纵坐标表示检测概率。图3(b)也表明，信号源方位角间距相同时，本发明方法比其他方法的检测概率更高，因此更加稳健。