用凸优化求解阵列信号处理：CVX方法在DOA估计中的应用

一、引言

在现代雷达、声呐、无线通信以及卫星导航等领域，波达方向（Direction of Arrival，DOA）估计 是阵列信号处理中的核心问题。它能够帮助我们确定多个信号源的入射方向，从而为目标定位、跟踪以及波束赋形提供关键信息。

传统的 DOA 估计算法主要包括 经典波束形成（Conventional Beamforming，CBF）、MUSIC（Multiple Signal Classification） 和 ESPRIT（Estimation of Signal Parameters via Rotational Invariance Techniques） 等。这些方法在实际应用中取得了广泛成功，但也存在一些局限性：

• 分辨率受限：当入射信号的角度间隔较小时，传统方法往往难以分辨。
• 对噪声与相干信号敏感：特征分解类方法在低信噪比或相干信号环境下性能显著下降。
• 需要信源数先验：如 MUSIC 需要准确估计信号源个数，否则结果会失真。

随着 稀疏表示理论 与 压缩感知（Compressed Sensing, CS） 的发展，研究者逐渐发现 DOA 估计问题可以建模为一个 稀疏重构问题：在一组离散化的角度栅格上，真实信号源数量远少于候选方向总数，因此可以通过稀疏优化方法恢复信号功率分布。

在这一思路下，凸优化（Convex Optimization） 成为解决 DOA 问题的重要工具。通过对稀疏信号的 L1 范数最小化约束，可以在噪声环境下获得稳定且分辨率更高的 DOA 估计结果。CVX 作为 MATLAB 平台上成熟的凸优化工具箱，为这种方法的实现提供了便捷途径。

本文将系统介绍如何将 DOA 估计转化为稀疏恢复问题，并利用 CVX 工具箱进行建模与求解，从而展示凸优化方法在阵列信号处理中的优势。

二、阵列信号模型

在讨论凸优化方法之前，我们首先需要建立阵列接收信号的数学模型。以最常见的 均匀线阵（Uniform Linear Array, ULA） 为例，假设阵列共有 $M$ 个阵元，阵元间距为 $d$，有 $P$ 个窄带平面波从方向 {θ1,θ2,…,θP}\{\theta_1, \theta_2, \dots, \theta_P\}{θ1​,θ2​,…,θP​} 入射，则第 $m$ 个阵元接收到的信号可以表示为：

$$x_m(t)=\sum _{p=1}^P{s}_p(t)e^{-j2\pi \frac{d}{\lambda }(m-1)\sin {\theta }_p}+n_m(t)$$

其中：

• $s_p(t)$ 表示第 $p$ 个入射信号；
• $\lambda$ 为载波波长；
• $n_m(t)$ 为噪声；
• 指数项代表 阵列导向矢量（Steering Vector） 对应的相位延迟。

将所有阵元的接收信号写成向量形式，可得阵列信号模型：

$$\mathbf{x}(t)=\mathbf{A}(\theta )\mathbf{s}(t)+\mathbf{n}(t)$$

• $\mathbf{x}(t)\in {\mathbb{C}}^{M\times 1}$：阵列接收向量；
• $\mathbf{A}(\theta )=\left[\mathbf{a}({\theta }_1),\mathbf{a}({\theta }_2),\dots ,\mathbf{a}({\theta }_P)\right]\in {\mathbb{C}}^{M\times P}$ 为 阵列流形矩阵；
• $\mathbf{a}({\theta }_p)={\left[1,e^{-j2\pi \frac{d}{\lambda }\sin {\theta }_p},\dots ,e^{-j2\pi \frac{d}{\lambda }(M-1)\sin {\theta }_p}\right]}^T$ 为方向 ${\theta }_p$ 的导向矢量；
• $\mathbf{s}(t)\in {\mathbb{C}}^{P\times 1}$ 为信号源向量；
• $\mathbf{n}(t)\in {\mathbb{C}}^{M\times 1}$ 为噪声向量，通常假设为复高斯白噪声。

协方差矩阵

对信号取快拍平均，可以得到样本协方差矩阵：

$${\mathbf{R}}_x=\mathbb{E}\left[\mathbf{x}(t){\mathbf{x}}^H(t)\right]\approx \frac{1}{L}\sum _{l=1}^L\mathbf{x}(l){\mathbf{x}}^H(l)$$

展开可得：

$${\mathbf{R}}_x=\mathbf{A}(\theta ){\mathbf{R}}_s{\mathbf{A}}^H(\theta )+{\sigma }_n^2{\mathbf{I}}_M$$

其中：

• ${\mathbf{R}}_s=\mathbb{E}[\mathbf{s}(t){\mathbf{s}}^H(t)]$ 表示信号源协方差矩阵；
• ${\sigma }_n^2$ 为噪声功率；
• ${\mathbf{I}}_M$ 为 $M\times M$ 阵列单位矩阵。

通过上述模型，DOA 估计问题可以转化为：在候选角度集合中，找到能够最好解释观测数据 $\mathbf{x}(t)$ 或 ${\mathbf{R}}_x$ 的方向集合。在后续章节中，我们将进一步展示如何将这一问题建模为稀疏优化问题，并借助 CVX 工具箱求解。

三、稀疏建模与DOA估计

在前文的阵列信号模型中，我们已经得到接收信号的表达式：

$$\mathbf{x}(t)=\mathbf{A}(\theta )\mathbf{s}(t)+\mathbf{n}(t)$$

其中 $\mathbf{A}(\theta )$ 是阵列流形矩阵，包含了所有真实信源方向 {θ1,…,θP}\{\theta_1,\dots,\theta_P\}{θ1​,…,θP​} 对应的导向矢量。

然而，在实际应用中，我们并不知道真实的 ${\theta }_p$，只能在一个较密的角度网格 {θ~1,θ~2,…,θ~N}\{\tilde{\theta}_1,\tilde{\theta}_2,\dots,\tilde{\theta}_N\}{θ~1​,θ~2​,…,θ~N​} 上进行扫描。通常 $N\gg P$，此时得到的字典矩阵为：

$$\mathbf{A}=[\mathbf{a}({\tilde{\theta }}_1),\mathbf{a}({\tilde{\theta }}_2),\dots ,\mathbf{a}({\tilde{\theta }}_N)]\in {\mathbb{C}}^{M\times N}$$

这样，接收信号可以表示为：

$$\mathbf{x}(t)\approx \mathbf{A}\mathbf{s}(t)+\mathbf{n}(t)$$

其中 $\mathbf{s}(t)\in {\mathbb{C}}^{N\times 1}$ 在绝大多数角度对应的分量为零，只有在真实信源方向附近才有显著能量。换句话说，$\mathbf{s}(t)$ 是一个稀疏向量。

稀疏建模的核心思想

• 角度离散化：将连续的 DOA 参数空间离散为有限个候选角度点；
• 稀疏先验：假设实际信源数目远小于候选角度数目，因此 $\mathbf{s}(t)$ 在字典空间中是稀疏的；
• 观测与重构：通过观测 $\mathbf{x}(t)$ 或其统计量 ${\mathbf{R}}_x$，在约束条件下恢复稀疏的 $\mathbf{s}(t)$，其非零位置即为信源入射方向。

特征向量投影

在具体实现中，一个常用做法是取协方差矩阵 ${\mathbf{R}}_x$ 的最大特征向量 $\mathbf{u}$ 作为观测量。这样，可以将 DOA 估计问题转化为：

$$\mathbf{u}\approx \mathbf{A}\mathbf{s}$$

其中 $\mathbf{s}$ 是一个稀疏向量，代表不同扫描角度上的信号强度。

因此，DOA 问题被改写为一个稀疏信号恢复问题：在稀疏约束下寻找 $\mathbf{s}$，使得 $\mathbf{As}$ 尽可能接近 $\mathbf{u}$。

四、凸优化建模

在上一节中，我们将 DOA 问题建模为一个稀疏信号恢复问题：

$$\mathbf{u}\approx \mathbf{A}\mathbf{s}$$

其中：

• $\mathbf{u}\in {\mathbb{C}}^{M\times 1}$ 为观测向量（通常取协方差矩阵的最大特征向量）；
• $\mathbf{A}\in {\mathbb{C}}^{M\times N}$ 为扫描角度对应的阵列字典矩阵；
• $\mathbf{s}\in {\mathbb{C}}^{N\times 1}$ 为稀疏解，非零位置对应实际 DOA。

1. 理想的稀疏优化模型

在稀疏建模框架下，我们希望找到最稀疏的解：

$$\underset{\mathbf{s}}{\min }\Vert \mathbf{s}{\Vert }_0\text{s.t.}\Vert \mathbf{u}-\mathbf{As}{\Vert }_2^2\le ϵ$$

其中 $\Vert \mathbf{s}{\Vert }_0$ 表示向量中非零元素的个数，$ϵ$ 控制拟合误差范围。

然而，L0 范数优化是 NP-困难问题，在实际中无法直接求解。

2. L1 范数替代与凸优化

为了获得可解的形式，通常采用 L1 范数替代 L0 范数，从而将问题转化为一个凸优化问题：

$$\underset{\mathbf{s}}{\min }\Vert \mathbf{s}{\Vert }_1\text{s.t.}\Vert \mathbf{u}-\mathbf{As}{\Vert }_2^2\le ϵ$$

• L1 范数 $\Vert \mathbf{s}{\Vert }_1={\sum }_i|s_i|$ 能够有效促进稀疏解；
• 约束 $\Vert \mathbf{u}-\mathbf{As}{\Vert }_2^2\le ϵ$ 保证解与观测一致性。

该问题是 凸优化问题，可通过标准凸优化工具（如 CVX）高效求解。

3. 一般化的 P-范数模型

更一般地，可以考虑 P 范数 ($p\ge 1$) 作为正则化目标：

$$\underset{\mathbf{s}}{\min }\Vert \mathbf{s}{\Vert }_p^p\text{s.t.}\Vert \mathbf{u}-\mathbf{As}{\Vert }_2^2\le ϵ$$

• 当 $p=1$ 时，得到经典的稀疏恢复模型（Basis Pursuit Denoising, BPDN）；
• 当 $p=2$ 时，优化更平滑，但稀疏性下降；
• 当 $1<p<2$ 时，提供稀疏性与稳定性的折中。

4. CVX 求解框架

在 MATLAB 的 CVX 工具箱中，可以将上述问题写作：

cvx_begin
    variable s(N) complex
    minimize( norm(s,1) )                 % L1 范数目标
    subject to
        norm(u - A*s,2) <= sqrt(epsilon)  % 残差约束
cvx_end

其中：

• s(N) 为待求的稀疏系数向量；
• norm(s,1) 对应 L1 范数目标；
• norm(u - A*s,2) <= sqrt(epsilon) 控制误差范围。

五、MATLAB 实现与仿真

%% ==============================================================
% 说明：
% - 需预先安装并配置 CVX（http://cvxr.com/cvx/）
% - 本脚本包含：数据仿真、字典构造、CVX求解、绘图
% - 不考虑相干信源（独立复高斯）
% ==============================================================

close all; clear all;clc;
rng(2); % 固定随机种子

%% ---------------- 参数定义区 ----------------
c = 299792458; % 光速 (m/s)，避免依赖 physconst
f0 = 10e9; % 载波频率 (Hz)
fg = 50e9; % 全局采样率 (Hz) —— 本窄带模型不直接用到
BW = 1e8; % 信号带宽 (Hz) —— 本窄带模型不直接用到
src_angle_deg = [-10, -20]; % 目标来向角 (度)
src_angle = deg2rad(src_angle_deg);
M = 32; % 阵元数
L = 128; % 快拍数
snr = 20; % 信噪比 (dB)
scan_angle_deg = -60:0.1:60; % 扫描角度 (度)
scan_angle = deg2rad(scan_angle_deg);
% 非相干信源
is_coherent = false;

%% ---------------- 参数计算区 ----------------
lambda = c / f0; % 波长 (m)
d = lambda / 2; % 阵元间距 (m)
P = numel(src_angle); % 信源数

%% ---------------- 回波生成 ----------------
% x_sig: M×L 快拍矩阵；sigma: 噪声方差；R_sig: 样本协方差矩阵
[x_sig, sigma, R_sig] = echo_generate(M, d, lambda, src_angle, L, snr, is_coherent);

%% ---------------- 字典与观测构造 ----------------
% 扫描角度的阵列流形字典（ULA，半波间距，窄带近似）
i = (0:M-1).';
A = exp(-1i*2*pi*(d/lambda) * i * sin(scan_angle)); % M×Ns

% 取 Rxx 的主特征向量作为观测向量 u
[U, D] = eig((R_sig+R_sig')/2); % 数值稳健：强制Hermitian
[~, idx] = max(real(diag(D)));
u = U(:, idx);

%% ---------------- CVX 稀疏恢复 ----------------
cvx_check_or_error(); % 没装CVX会抛错并提示安装地址
tor_lim = 1e-1; % 残差容许上限（发散时可适当增大）
P_cvx_1 = DOA_CVX(u, A, 1.0, tor_lim); % L1 范数
P_cvx_1_5 = DOA_CVX(u, A, 1.5, tor_lim); % L1.5 范数（p>=1 保持凸）
P_cvx_2 = DOA_CVX(u, A, 2.0, tor_lim); % L2 范数（更平滑，稀疏性弱）

%% 转换为 dB 值
P_cvx_1_db = 10*log10(P_cvx_1 + eps);
P_cvx_1_5_db = 10*log10(P_cvx_1_5 + eps);
P_cvx_2_db = 10*log10(P_cvx_2 + eps);

%% ---------------- 结果绘图 ----------------
figure('Color','w'); hold on; grid on;
plot(scan_angle_deg, P_cvx_1_db, 'LineWidth', 1.6);
plot(scan_angle_deg, P_cvx_1_5_db, 'LineWidth', 1.6);
plot(scan_angle_deg, P_cvx_2_db, 'LineWidth', 1.6);
yline(0,'k:');
xlabel('Angle (deg)'); ylabel('Normalized spectrum(dB)');
title('DOA via Sparse Recovery (CVX)');
legend({'p=1','p=1.5','p=2'}, 'Location','best');
xlim([min(scan_angle_deg) max(scan_angle_deg)]);
set(gca,'FontName','Arial','FontSize',11);
% 标注真实来向
for ang = src_angle_deg
xline(ang,'r--','LineWidth',1);
end
text(src_angle_deg+0.5, 0.95*ones(size(src_angle_deg)), ...
compose('%d° (true)', src_angle_deg), 'Color',[0.85 0 0], 'FontSize',10);
fprintf('完成：CVX-DOA 估计绘制完毕。\n');

%%函数区
function [x_sig, sigma, R_sig] = echo_generate(M, d, lambda, src_angle, L, snr_dB, is_coherent)
% 窄带 ULA 数据生成（非相干或相干，默认非相干）
% x = A s + n，s 为 P×L 独立复高斯（非相干），n ~ CN(0, sigma I)
P = numel(src_angle);
i = (0:M-1).';
A = exp(-1i*2*pi*(d/lambda) * i * sin(src_angle(:).')); % M×P
if is_coherent
% 相干：所有信源共享同一基带序列（仅相位/幅度不同）
base = (randn(1, L) + 1i*randn(1, L)) / sqrt(2);
s = repmat(base, P, 1);
else
% 非相干：各源独立复高斯，单位功率
s = (randn(P, L) + 1i*randn(P, L)) / sqrt(2);
end
% 设每个源单位功率，总接收信号（不含噪声）每通道功率 ~ P
% 依据 SNR 定噪声方差 sigma
SNR = 10.^(snr_dB/10);
sigma = P / SNR; % 每通道噪声方差
n = sqrt(sigma/2) * (randn(M, L) + 1i*randn(M, L));
x_sig = A*s + n;
R_sig = (x_sig * x_sig') / L; % 样本协方差
end

function s = DOA_CVX(u, A, p_norm, tor_lim)
% 基于 CVX 的稀疏恢复：
% min ||s||_p s.t. ||u - A s||_2^2 <= tor_lim
% 返回值做幅度归一化，便于绘图
Ns = size(A, 2);
% CVX 求解（允许复变量）
cvx_begin quiet
variable s_cvx(Ns) complex
minimize( sum( pow_abs(s_cvx, p_norm) ) )
subject to
sum( pow_abs(u - A * s_cvx, 2) ) <= tor_lim
cvx_end
s = abs(s_cvx);
if max(s) > 0
s = s / max(s);
end
end

function cvx_check_or_error()
% 检查 CVX 是否可用，否则抛出友好错误
if ~(exist('cvx_begin','file') == 2 || exist('cvx_setup','file') == 2)
error(['未检测到 CVX。请先安装并 addpath：\n' ...
' 下载： http://cvxr.com/cvx/\n' ...
' 安装： 解压 -> addpath(genpath(''cvx'')) -> cvx_setup']);
end
end

六、总结

本文围绕 凸优化方法在 DOA 估计中的应用 展开，从阵列信号模型出发，将传统 DOA 问题转化为 稀疏恢复问题，并进一步通过 L1 范数最小化 构建凸优化模型。在该框架下，利用 CVX 工具箱可以方便地实现求解，获得比传统方法更高分辨率和更强鲁棒性的方向估计结果。

与 MUSIC、ESPRIT 等特征分解类方法相比，基于凸优化的 DOA 估计具有以下特点：

• 高分辨率：能够区分角度间隔较小的信源；
• 鲁棒性强：在低信噪比或有限快拍数下仍能稳定工作；
• 无需信源数先验：稀疏建模自然地反映了信号稀疏性；
• 计算复杂度高：优化问题求解耗时较长，不适合实时系统。

总体而言，凸优化方法为 DOA 估计提供了一种新的思路和工具，尤其适合需要高精度离线处理的场景。未来，可以进一步结合 快速稀疏重构算法（如 OMP、FOCUSS、ADMM 等）来降低计算开销，从而推动该方法在工程中的实际应用。