CN106842113A - 高采样1比特量化情况下的信号到达角高精度估计方法 - Google Patents

Abstract

本发明公开了一种高采样1比特量化情况下的信号到达角高精度估计方法，首先产生由M个阵元构成的均匀阵列接收的一个快拍信号并将该信号1比特量化采样；然后建立以p范数为稀疏约束项并结合不敏感支持向量机回归的ε‑SVR模型的稀疏信号表示模型，利用该非凸优化模型可以稀疏表示原信号信息；最后运用基于交替方向乘子法(ADMM)求解该非凸优化模型，求出稀疏表示系数，由稀疏表示系数确定信号的波达方向。本发明方法，能在采样量少(如单快拍)，无关信源是否相干，不需计算自相关矩阵的情况下，高精度估计信号波达方向，并且基于1比特量化采样情况下的信号波达方向估计更易于在工程中快速实现。

Description

高采样1比特量化情况下的信号到达角高精度估计方法

技术领域

本发明属于空间谱估计领域，具体涉及一种高采样1比特量化情况下的信号到达角高精度估计方法。

背景技术

空间谱估计是利用在空间中，以不同排列方式组合在一起的接收阵列，对信源信息进行接收并处理，从而获得信源参数的技术。空间谱估计对信源空域参数，如信源个数，信号波达方向(DOA)及俯仰角的优越估计性能，已经使其广泛应用在通信，雷达，声呐及生物医学等领域。常见的空间阵列排列方式有，均匀线阵，均匀圆阵，平面阵。其中，均匀线阵在实际分析中应用广泛。本发明后续采用均匀线阵进行实例的仿真实验。

对于信号波达方向(DOA)的估计问题，有很多学者从不同方面提出了空间谱估计的方法，有ARMA谱分析法，最大似然法，熵谱分析法和特征分解法等，其中基于特征分解的子空间算法主要有MUSIC，ESPRIT，WSF等，但是这些算法在有些信号模型下会失效，如信源是相干的，并且在实际环境中信噪比低，采样数不够多，获取的自相关矩阵不准确时，这些算法的估计性能也会大大下降。近年来，随着稀疏理论的不断完善，利用稀疏信号表达的方法对未知信号波达方向的估计已经取得了进步，避免了自相关矩阵的运算，也不受信源相干性的影响，同时也提高了估计性能。

在实际应用中，阵列接收机接收到的数据需要被量化。1比特采样量化就是利用1位二进制数据来表示观测信号中的符号信息，在工程实现中，也不需要复杂的电路结构，只需要一个比较器就可以完成，大大降低了工程实现成本。

因此，采用1比特采样量化并稀疏重构原始信号信息的方法比传统意义上的DOA估计算法，估计精确要好，采样要求低，工程实现简单，快速。

发明内容

本发明的目的是提供一种高采样1比特量化情况下的信号到达角高精度估计方法，解决了在低信噪比下，采样量少，自相关矩阵获取不准确的问题。

本发明所采用的技术方案是，高采样1比特量化情况下的信号到达角高精度估计方法，具体包括以下步骤：

步骤1，产生由M个阵元构成的均匀阵列，其仅接收的一个快拍信号y，该快拍信号经过1比特量化采样形成；

步骤2，建立以p范数为稀疏信号表达方式和外加抗噪声干扰的不敏感支持向量机回归的ε-SVR模型，利用该非凸优化模型构建原始信号稀疏表示方式；

步骤3，运用基于交替方向乘子法(ADMM)求解该非凸优化模型，求出稀疏表示系数，由稀疏表示系数确定信号的波达方向。

本发明的有益效果是，本发明高采样1比特量化情况下的信号到达角高精度估计方法，能在采样量少(如单快拍)，无关信源是否相干，不需计算自相关矩阵的情况下，高精度估计信号波达方向，并且基于1比特量化采样情况下的信号波达方向估计更易于在工程中实现。

附图说明

图1是本发明用于接收单快拍信号的均匀阵列结构图；

图2是本发明方法得到的实际信号波达方向的估计图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明高采样1比特量化情况下的信号到达角高精度估计方法，具体包括以下步骤：

步骤1，产生由M个阵元构成的均匀阵列，其仅接收的一个快拍信号y，该快拍信号经过1比特量化采样形成；

假设接收信号满足窄带信号，即信号经过阵列长度的时间远小于信号的相干时间，如图1所示，来波方向为θ_i(i＝1,2,....I)，入射N根天线，阵元间隔为d的均匀线阵的阵列响应矢量为：

定义方向矩阵为：

相对应的信号幅度为则经采样1比特量化后的信号y：

其中V是高斯白噪声,Tr是给定阈值，若测量值实部(虚部)大于该阈值实部(虚部)，则y的实部(虚部)为1，否则为-1。

步骤2，建立以p范数为稀疏信号表达方式和外加抗噪声干扰的不敏感支持向量机回归的ε-SVR模型，利用该非凸优化模型构建原始信号稀疏表示方式；

考虑这N个信号的入射角取值范围为：[-90 90]，为此，将[-90 90]均匀细分为181个均匀空间并建立导向矢量过完备字典：

A＝[a(θ₁) a(θ₂) … a(θ₁₈₁)]∈C^M×181 (4)

相对应的信号幅度S＝(s₁,s₂,…s₁₈₁)；为了从单快拍y中估计θ_i，建立以p范数为稀疏约束项和不敏感ε-SVR关于角度估计的稀疏表示模型：

(·)_n,r和(·)_n,i分别表示·的第n分量的实部和虚部，ζ_n,ζ^* _n,ε>0。

步骤3，运用基于交替方向乘子法(ADMM)求解该非凸优化模型，求出稀疏表示系数，由稀疏表示系数确定信号的波达方向。

(5)式是一个非凸优化模型，使用交替方向乘子法(ADMM)迭代算法求解在过完备字典A(基)下的稀疏表示系数S，从而得出信号的波达方向。具体步骤如下：

首先，根据ADMM，引入辅助变量t＝S，构造如下的拉格朗日函数：

关于该模型p,ρ,c,ε的选取很重要，惩罚因子c>0，控制对错分样本的惩罚程度，c值过大过小都会影响系统的泛化能力；不敏感损失系数ε过小，则回归精度高，但支持向量增多，反之，回归精度降低；p＝1，则演变为S的1范数，p＝0，则演变为S的0范数，0范数表示向量非零元素的个数，实际上0范数不可解，所以，为了更稀疏约束信号，选取0<p<1，p越接近0，S越稀疏；ρ控制||S||^p和的均衡，以及控制ADMM算法迭代步长；

然后，分布优化如下：

①固定t,λ，求解S，则(6)式等价为：

为进一步简化(7)，即：

其中，I＝181，显然(8)式取最小值时，可以得出S与同相位，所以直接考虑复数模的运算，并把(8)式转化为I个子函数求解问题：

其中是t,S,λ的模，

因为在上分别单调递增和递减，所以的全局最小值出现在上，即(9)式等价于：

现在讨论在上的凹凸性，并分别对求1阶导，2阶导，3阶导：

对p分三种情况讨论(10)式的最小值点：

(1)1<p≤2

因为时，故是凸的，并且 所以在区间有唯一解，使用二分法对(11)求零点即为的最小值点；

(2)p＝1

当p＝1时(10)式等价于：

则的最小值点为：

(3)0<p<1

0<p<1时故是单增的，的凹凸性取决于的符号，记的解为：

当 是凹的，故的最小值点在0或之间；

当 可能是正的或负的，分和两个区间讨论：

i.在上， 是凹的，故的最小值点在0或之间；

ii.在上， 是凸的，且如果 是单增的，故上的局部最小值点如果使用二分法求取上的局部最小值点；

显然，的全局最小值点通过比较和上的局部最小值而获得；

②固定S,λ，求取t，则(6)式等价为：

目标函数进一步化简为：

其中令w＝t-r，则(17)式等价于：

其中进一步令real(·)和imag(·)分别表示·的实部与虚部，则(19)化简为：

其中，H＝[real(A)-imag(A)]∈C^M×362，H_n∈C^1×362为H的第n行，R＝[imag(A)real(A)]∈C^M×362，R_n∈C^1×362为R的第n行，D＝real(Ar+V-Tr)∈C^M×1，D1＝imag(Ar+V-Tr)∈C^M×1；

利用拉格朗日乘子法来解决这个约束优化问题，构造如下的拉格朗日函数：

其中，是拉格朗日乘子；

根据最优化理论，L分别对W,ξ,ξ^*求偏微分，并令所得等式为0，得到下式：

将(22)式带入(21)式得到对偶优化问题：

通过矩阵简化，可以转化(23)式为：

其中Q＝[H^T*diag(real(y)R^T*diag(imag(y)]∈C^362×2M，

P＝[D^T*diag(real(y)D1^T*diag(imag(y)]∈C^1×2M，T＝[real(y)^T imag(y)^T]

公式(24)用CVX优化工具包，求解出β，进而解出：

W＝Qβ (25)

由可以得出w，利用w＝t-r，求出：

t(k+1)＝w+r (26)

k为迭代次数；

③固定S,t，求λ，即：

λ(k+1)＝λ(k)+ρ(t-S) (27)

最后，①，②，③步循环迭代一定次数，直到S,λ,t收敛，即求得信号的波达方向的稀疏表示系数S，然后作出信号角度频谱图，即S的幅值曲线，尖峰处角度，即为估计的信号到达角。

实施例1

步骤1：有M＝20个阵元，其仅接受一个快拍信号y，噪声方差为0.01，信噪比为16.99db，设定阈值Tr为接受信号的最大值与最小值之间的随机值。假设信号元来波方向为-20°,-40°，经过(3)式1比特采样量化得到y。

步骤2：建立以p范数为稀疏约束项和不敏感ε-SVR关于角度估计的稀疏表示模型：

(·)_n,r和(·)_n,i分别表示·的第n分量的实部和虚部，

步骤3：运用基于交替方向乘子法(ADMM)求解该非凸优化模型，求出稀疏表示系数，由稀疏表示系数确定信号的波达方向。

根据ADMM，引入辅助变量t＝S,构造如下的拉格朗日函数：

因为ADMM算法是迭代算法，所以首先初始化所求参数t,λ，然后分布优化如下：

①固定t,λ，求解S，则(29)式等价为：

为进一步简化(30)，即：

其中，I＝181，显然(31)式取最小值时，可以得出S与同相位，所以直接考虑复数模的运算，并把(31)式转化为I个子函数求解问题：

其中是t,S,λ的模，

因为在上分别单调递增和递减，所以的全局最小值出现在上，即(32)式等价于：

现在讨论在上的凹凸性，并分别对求1阶导，2阶导，3阶导：

对p分三种情况讨论(33)式的最小值点：

(1)1<p≤2

因为时，故是凸的，并且 所以在区间有唯一解，使用二分法对(34)求零点即为的最小值点。

(2)p＝1

当p＝1时(33)式等价于：

则的最小值点为：

(3)0<p<1

0<p<1时故是单增的，的凹凸性取决于的符号，记的解为：

当 是凹的，故的最小值点在0或之间。

当 可能是正的或负的，分和两个区间讨论：

i.在上，是凹的，故的最小值点在0或之间。

ii.在上,是凸的，且如果 是单增的，故上的局部最小值点如果使用二分法求取上的局部最小值点。

显然，的全局最小值点通过比较和上的局部最小值而获得。

②固定S,λ，求取t,则(29)式等价为：

目标函数进一步化简为：

其中令w＝t-r，则(40)式等价于：

其中进一步令real(·)和imag(·)分别表示·的实部与虚部，则(42)化简为：

其中，H＝[real(A)-imag(A)]∈C^20×362，H_n∈C^1×362为H的第n行，R＝[imag(A)real(A)]∈C^20×362，R_n∈C^1×362为R的第n行，D＝real(Ar+V-Tr)∈C^20×1，D1＝imag(Ar+V-Tr)∈C^{20 ×1}。

利用拉格朗日乘子法来解决这个约束优化问题，构造如下的拉格朗日函数：

其中，是拉格朗日乘子。

根据最优化理论，L分别对W,ξ,ξ^*求偏微分，并令所得等式为0，得到下式：

将(45)式带入(44)式得到对偶优化问题：

通过矩阵简化，可以转化(46)式为：

其中Q＝[H^T*diag(real(y)R^T*diag(imag(y)]∈C^362×40，

P＝[D^T*diag(real(y)D1^T*diag(imag(y)]∈C^1×40，T＝[real(y)^T imag(y)^T]

公式(47)用CVX优化工具包，可以求解出β，进而解出：

W＝Qβ (48)

由可以得出w，利用w＝t-r，求出：

t(k+1)＝w+r (49)

③固定S,t，求λ，即：

λ(k+1)＝λ(k)+ρ(t-S) (50)

①，②，③步循环迭代一定次数k，直到S,λ,t收敛，即求得S，

作出S曲线，尖峰点处的角度即所对应的信号到达角。如图2所示，信号到达角为-20°,-40°，高精度估计了这两个信源方向。

Claims (4)

1.高采样1比特量化情况下的信号到达角高精度估计方法，其特征在于，具体包括以下步骤：

步骤1，产生由M个阵元构成的均匀阵列，其仅接收的一个快拍信号y，该快拍信号经过1比特量化采样形成；

步骤2，建立以p范数为稀疏信号表达方式和外加抗噪声干扰的不敏感支持向量机回归的ε-SVR模型，利用该非凸优化模型构建原始信号稀疏表示方式；

步骤3，运用基于交替方向乘子法(ADMM)求解该非凸优化模型，求出稀疏表示系数，由稀疏表示系数确定信号的波达方向。

2.根据权利要求1所述的高采样1比特量化情况下的信号到达角高精度估计方法，其特征在于，步骤1具体步骤为：

假设接收信号满足窄带信号，即信号经过阵列长度的时间远小于信号的相干时间，来波方向为θ_i(i＝1,2,....I)，入射N根天线，阵元间隔为d的均匀线阵的阵列响应矢量为：

$a\left({\mathrm{\&theta;}}_i\right)={\left[\begin{array}{cccc}1& \exp \left(-j2\mathrm{\&pi;}\frac{d}{\mathrm{\&lambda;}}{\mathrm{sin\&theta;}}_i\right)& \mathrm{...}& \exp \left(-j2\mathrm{\&pi;}\left(M-1\right)\frac{d}{\mathrm{\&lambda;}}{\mathrm{sin\&theta;}}_i\right)\end{array}\right]}^T---\left(1\right)$

定义方向矩阵为:

相对应的信号幅度为则经采样1比特量化后的信号y：

$y=sign(\tilde{A}\tilde{S}+V-Tr)---\left(3\right)$

其中V是高斯白噪声，Tr是给定阈值，若测量值实部(虚部)大于该阈值实部(虚部)，则y的实部(虚部)为1，否则为-1。

3.根据权利要求1或2所述的高采样1比特量化情况下的信号到达角高精度估计方法，其特征在于，步骤2具体步骤为：

考虑N个信号的入射角取值范围为：[-90 90]，为此，将[-90 90]均匀细分为181个均匀空间并建立导向矢量过完备字典：

A＝[a(θ₁) a(θ₂) … a(θ₁₈₁)]∈C^M×181 (4)

相对应的信号幅度S＝(s₁,s₂,…s₁₈₁)；为了从单快拍y中估计θ_i，建立以p范数为稀疏约束项和不敏感ε-SVR关于角度估计的稀疏表示模型：

$\begin{array}{cc}\min & \left|\right|S|{|}^p+c\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&xi;}}_n+{{\mathrm{\&xi;}}_n}^{*}\right)\\ s.t& y_{n,r}\&lsqb;{\left(AS+V-Tr\right)}_{n,r}+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0\\ & \begin{array}{cc}{y}_{n,i}\&lsqb;{\left(AS+V-Tr\right)}_{n,i}+{{\mathrm{\&zeta;}}^{*}}_n-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0& n=1,2,\mathrm{...}M\end{array}\end{array}---\left(5\right)$

(·)_n,r和(·)_n,i分别表示·的第n分量的实部和虚部，ζ_n,ζ^* _n,ε>0。

4.根据权利要求3所述的高采样1比特量化情况下的信号到达角高精度估计方法，其特征在于，步骤3具体步骤如下：

首先，根据ADMM，引入辅助变量t＝S，构造如下的拉格朗日函数：

$\begin{array}{cc}\min & \left|\right|S|{|}^p+c\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&zeta;}}_n+{{\mathrm{\&xi;}}_n}^{*}\right)+{\mathrm{\&lambda;}}^T\left(t-S\right)+\frac{\mathrm{\&rho;}}{2}||t-S|{|}^2\\ s.t& y_{n,r}\&lsqb;{\left(At+V-Tr\right)}_{n,r}+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0\\ & \begin{array}{cc}{y}_{n,i}\&lsqb;{\left(At+V-Tr\right)}_{n,i}+{{\mathrm{\&zeta;}}_n}^{*}-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0& n=1,2,\mathrm{...}M\end{array}\end{array}---\left(6\right)$

关于该模型p,ρ,c,ε的选取很重要，惩罚因子c>0，控制对错分样本的惩罚程度，c值过大过小都会影响系统的泛化能力；不敏感损失系数ε过小，则回归精度高，但支持向量增多，反之，回归精度降低；p＝1，则演变为S的1范数，p＝0，则演变为S的0范数，0范数表示向量非零元素的个数，实际上0范数不可解，所以，为了更稀疏约束信号，选取0<p<1，p越接近0，S越稀疏；ρ控制||S||^p和的均衡，以及控制ADMM算法迭代步长；

然后，分布优化如下：

①固定t,λ，求解S，则(6)式等价为：

$\begin{array}{cc}\min & \left|\right|S|{|}^p+{\mathrm{\&lambda;}}^T(t-S)+\frac{\mathrm{\&rho;}}{2}||t-S|{|}^2\end{array}---\left(7\right)$

为进一步简化(7)，即：

$\begin{array}{cc}min& \underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}\left(\right|S_i{|}^p+\frac{\mathrm{\&rho;}}{2}{\left(S_i-{\hat{S}}_i\right)}^2)\end{array}---\left(8\right)$

其中，I＝181，显然(8)式取最小值时，得出S与同相位，所以直接考虑复数模的运算，并把(8)式转化为I个子函数求解问题：

$\begin{array}{cc}\min & h\left({\stackrel{\&OverBar;}{S}}_i\right)=|{\stackrel{\&OverBar;}{S}}_i{|}^p+\frac{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}{2}{\left({\stackrel{\&OverBar;}{S}}_i-{\hat{S}}_i\right)}^2\end{array}---\left(9\right)$

其中是t,S,λ的模，

因为在上分别单调递增和递减，所以的全局最小值出现在上，即(9)式等价于：

$\begin{array}{cccc}\min & h\left({\stackrel{\&OverBar;}{S}}_i\right)={{\stackrel{\&OverBar;}{S}}_i}^p+\frac{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}{2}{\left({\stackrel{\&OverBar;}{S}}_i-{\hat{S}}_i\right)}^2,& s.t& {\stackrel{\&OverBar;}{S}}_i\&Element;\&lsqb;0,{\hat{S}}_i\&rsqb;\end{array}---\left(10\right)$

现在讨论在上的凹凸性，并分别对求1阶导，2阶导，3阶导：

$\frac{\&part;h\left({\stackrel{\&OverBar;}{S}}_i\right)}{\&part;{\stackrel{\&OverBar;}{S}}_i}=p{{\stackrel{\&OverBar;}{S}}_i}^{p-1}+\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left({\stackrel{\&OverBar;}{S}}_i-{\hat{S}}_i\right)---\left(11\right)$

$\frac{{\&part;}^{2}h\left({\stackrel{\&OverBar;}{S}}_i\right)}{\&part;{{\stackrel{\&OverBar;}{S}}_i}^2}=p(p-1){{\stackrel{\&OverBar;}{S}}_i}^{p-2}+\stackrel{\&OverBar;}{\mathrm{\&rho;}}---\left(12\right)$

$\frac{{\&part;}^{3}h\left({\stackrel{\&OverBar;}{S}}_i\right)}{\&part;{{\stackrel{\&OverBar;}{S}}_i}^3}=p(p-1)(p-2){{\stackrel{\&OverBar;}{S}}_i}^{p-3}---\left(13\right)$

对p分三种情况讨论(10)式的最小值点：

(1)1<p≤2

因为时，故是凸的，并且 所以在区间有唯一解，使用二分法对(11)求零点即为的最小值点；

(2)p＝1

当p＝1时(10)式等价于：

$\begin{array}{cccc}\min & \frac{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}{2}{({\stackrel{\&OverBar;}{S}}_i-({\hat{S}}_i-\frac{1}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}\left)\right)}^2& s.t& {\stackrel{\&OverBar;}{S}}_i\&Element;\&lsqb;0,{\hat{S}}_i\&rsqb;\end{array}---\left(14\right)$

则的最小值点为：

${\stackrel{\&OverBar;}{S}}_i=\left\{\begin{array}{ccc}{\hat{S}}_i-\frac{1}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}& if& {\hat{S}}_i-\frac{1}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}>0\\ 0& if& {\hat{S}}_i-\frac{1}{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}\&le;0\end{array}\right.---\left(15\right)$

(3)0<p<1

0<p<1时故是单增的，的凹凸性取决于的符号，记的解为：

${\widetilde{\stackrel{~}{S}}}_i={\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&rho;}}}{p(1-p)}\right)}^{\frac{1}{p-2}}---\left(16\right)$

当 是凹的，故的最小值点在0或之间；

当 是正的或负的，分和两个区间讨论：

i.在上， 是凹的，故的最小值点在0或之间；

ii.在上， 是凸的，且如果 是单增的，故上的局部最小值点如果使用二分法求取上的局部最小值点；

显然，的全局最小值点通过比较和上的局部最小值而获得；

②固定S,λ，求取t，则(6)式等价为：

$\begin{array}{cc}\min & c\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&zeta;}}_n+{{\mathrm{\&xi;}}_n}^{*}\right)+{\mathrm{\&lambda;}}^T\left(t-S\right)+\frac{\mathrm{\&rho;}}{2}\left|\right|t-S|{|}^2\\ s.t& y_{n,r}\&lsqb;{\left(At+V-Tr\right)}_{n,r}+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0\\ & \begin{array}{cc}{y}_{n,i}\&lsqb;{\left(At+V-Tr\right)}_{n,i}+{\mathrm{\&zeta;}}_n^{*}-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0& n=1,2,\mathrm{...}M\end{array}\end{array}---\left(17\right)$

目标函数进一步化简为：

$\begin{array}{cc}\min & \frac{1}{2}\left|\right|t-r|{|}^2+\frac{c}{\mathrm{\&rho;}}\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&zeta;}}_n+{{\mathrm{\&xi;}}_n}^{*}\right)\end{array}---\left(18\right)$

其中令w＝t-r，则(17)式等价于：

$\begin{array}{cc}min& \frac{1}{2}\left|\right|w|{|}^2+c^{\&prime;}\underset{n=1}{\overset{M}{\&Sigma;}}({\mathrm{\&zeta;}}_n+{{\mathrm{\&xi;}}_n}^{*})\\ s.t& \left\{\begin{array}{c}{y}_{n,r}\&lsqb;{\left(Aw\right)}_{n,r}+{(Ar+V-Tr)}_{n,r}+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0\\ y_{n,i}\&lsqb;{\left(Aw\right)}_{n,i}+{(Ar+V-Tr)}_{n,i}+{{\mathrm{\&zeta;}}^{*}}_n-\mathrm{\&epsiv;}\&rsqb;\&GreaterEqual;0\end{array}\right.\end{array}---\left(19\right)$

其中进一步令real(·)和imag(·)分别表示·的实部与虚部，则(19)化简为：

$\begin{array}{cc}min& \frac{1}{2}\left|\right|W|{|}^2+c^{\&prime;}\underset{n=1}{\overset{M}{\&Sigma;}}({\mathrm{\&zeta;}}_n+{{\mathrm{\&xi;}}_n}^{*})\\ s.t& \left\{\begin{array}{c}{y}_{n,r}(H_{n}W+D_n+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;})\&GreaterEqual;0\\ y_{n,i}(R_{n}W+D{1}_n+{\mathrm{\&zeta;}}_n^{*}-\mathrm{\&epsiv;})\&GreaterEqual;0\end{array}\right.\end{array}---\left(20\right)$

其中，H＝[real(A)-imag(A)]∈C^M×362，H_n∈C^1×362为H的第n行，R＝[imag(A)real(A)]∈C^M×362，R_n∈C^1×362为R的第n行，D＝real(Ar+V-Tr)∈C^M×1，D1＝imag(Ar+V-Tr)∈C^M×1；

利用拉格朗日乘子法来解决这个约束优化问题，构造如下的拉格朗日函数：

$\begin{array}{c}L=\frac{1}{2}\left|\right|W|{|}^2+c^{\&prime;}\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&zeta;}}_n+{{\mathrm{\&xi;}}_n}^{*}\right)-\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&eta;}}_n{\mathrm{\&xi;}}_n+{\mathrm{\&eta;}}_n^{*}{\mathrm{\&xi;}}_n^{*}\right)-\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_n{y}_{n,r}\left(H_{n}W+D_n+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;}\right)\\ -\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_n^{*}{y}_{n,i}\left(R_{n}W+D{1}_n+{\mathrm{\&zeta;}}_n-\mathrm{\&epsiv;}\right)\end{array}---\left(21\right)$

其中，η_n,α_n,是拉格朗日乘子；

根据最优化理论，L分别对W,ξ,ξ^*求偏微分，并令所得等式为0，得到下式：

$\left\{\begin{array}{c}W=\underset{n=1}{\overset{M}{\&Sigma;}}({\mathrm{\&alpha;}}_n{y}_{n,r}{H}_n^T+{\mathrm{\&alpha;}}_n^{*}{y}_{n,i}{R}_n^T)\\ c^{\&prime;}-{\mathrm{\&eta;}}_n^{*}-{\mathrm{\&alpha;}}_n^{*}{y}_{n,i}=0\\ c^{\&prime;}-{\mathrm{\&eta;}}_n-{\mathrm{\&alpha;}}_n{y}_{n,r}=0\end{array}\right.---\left(22\right)$

将(22)式带入(21)式得到对偶优化问题：

$\begin{array}{cc}\min & \frac{1}{2}\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}{\left({\mathrm{\&alpha;}}_n{y}_{n,r}{H}_n^T+{\mathrm{\&alpha;}}_n^{*}{y}_{n,i}{R}_n^T\right)}^T\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&alpha;}}_m{y}_{m,r}{H}_m^T+{\mathrm{\&alpha;}}_m^{*}{y}_{m,i}{R}_m^T\right)+\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_n{y}_{n,r}{D}_n+{\mathrm{\&alpha;}}_n^{*}{y}_{n,i}D{1}_n\\ & -\mathrm{\&epsiv;}\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}\left({\mathrm{\&alpha;}}_n{y}_{n,r}+{\mathrm{\&alpha;}}_n^{*}{y}_{n,i}\right)\\ s.t& \begin{array}{cc}0\&le;{\mathrm{\&alpha;}}_n\&le;c^{\&prime;},0\&le;{\mathrm{\&alpha;}}_n^{*}\&le;c^{\&prime;}& n=1,\mathrm{2...}M\end{array}\end{array}---\left(23\right)$

通过矩阵简化，转化(23)式为：

$\begin{array}{cc}min& \frac{1}{2}{\mathrm{\&beta;}}^T{Q}^{T}Q\mathrm{\&beta;}+P\mathrm{\&beta;}-\mathrm{\&epsiv;}T\mathrm{\&beta;}\\ s.t& 0\&le;\mathrm{\&beta;}\&le;c^{\&prime;}\end{array}---\left(24\right)$

其中Q＝[H^T*diag(real(y)R^T*diag(imag(y)]∈C^362×2M，

P＝[D^T*diag(real(y)D1^T*diag(imag(y)]∈C^1×2M，T＝[real(y)^T imag(y)^T]

公式(24)用CVX优化工具包，求解出β，进而解出：

W＝Qβ (25)

由得出w，利用w＝t-r，求出：

t(k+1)＝w+r (26)

k为迭代次数；

③固定S,t，求λ，即：

λ(k+1)＝λ(k)+ρ(t-S) (27)

最后，①，②，③步循环迭代一定次数，直到S,λ,t收敛，即求得信号的波达方向的稀疏表示系数S，然后作出信号角度频谱图，即S的幅值曲线，尖峰处角度，即为估计的信号到达角。