MATLAB 实现 FOCUSS 算法：稀疏恢复 DOA 估计全流程

MATLAB 实现 FOCUSS 算法：稀疏恢复 DOA 估计全流程

一、引言

在现代雷达、声呐以及无线通信系统中，到达方向估计（Direction of Arrival, DOA） 是阵列信号处理中的核心问题之一。通过对接收信号的空间处理，可以估计目标信号的入射角度，从而实现目标探测、定位与跟踪。

传统的 DOA 估计算法包括 经典波束形成（Conventional Beamforming, CBF）、Capon（最小方差无畸变响应，MVDR） 以及 MUSIC（Multiple Signal Classification） 等方法。这些方法在快拍数充足、信噪比较高的情况下能够得到较好的性能，但在 低快拍数、强噪声或需要高分辨率 的场景下，其效果往往受到限制。

近年来，随着 稀疏表示与压缩感知理论 的发展，研究者们将稀疏恢复思想引入 DOA 估计问题。其基本出发点是：尽管扫描角度网格可能非常密集，但真实存在的信号源数量往往远小于栅格数，因此信号在角度域具有明显的稀疏性。基于此，可以通过稀疏优化的方法在高分辨率下重建信号源方向。

在众多稀疏恢复方法中，FOCUSS（Focal Underdetermined System Solver） 算法因其迭代加权最小二乘的形式而受到广泛关注。它通过不断调整加权矩阵，使得解逐渐收缩到稀疏方向，从而有效抑制旁瓣并提升分辨率。

本文将以 FOCUSS 算法为核心，结合公式推导与 MATLAB 仿真，完整展示一个稀疏恢复 DOA 估计的实现流程，并与简单的 伪逆（PINV）方法 进行对比分析。

二、阵列信号模型

在 DOA 问题中，我们通常考虑 均匀线阵（Uniform Linear Array, ULA） 的情况。假设阵列共有 $M$ 个阵元，阵元间距为 $d$，有 $P$ 个窄带远场信号从角度 {θ1,θ2,…,θP}\{\theta_1, \theta_2, \dots, \theta_P\}{θ1​,θ2​,…,θP​} 入射到阵列，则阵列接收信号可以建模为：

$$x(t)=\sum _{p=1}^P{s}_p(t)a({\theta }_p)+n(t)$$

其中：

• $s_p(t)$ 为第 $p$ 个信号源的基带信号；
• $a({\theta }_p)$ 为对应的阵列导向矢量；
• $n(t)$ 为噪声向量。

2.1 阵列导向矢量

对于 ULA，当信号以入射角 $\theta$ 到达时，其阵列导向矢量为：

$$a(\theta )=\left[\begin{array}{c}1\\ e^{-j\frac{2\pi d}{\lambda }\sin \theta }\\ e^{-j\frac{2\pi d}{\lambda }\cdot 2\sin \theta }\\ ⋮\\ e^{-j\frac{2\pi d}{\lambda }\cdot (M-1)\sin \theta }\end{array}\right]$$

其中 $\lambda =\frac{c}{f_0}$ 表示载波波长，$c$ 为光速，$f_0$ 为载波频率。

2.2 阵列流形矩阵（超完备字典）

若对一组扫描角度 {θ1,θ2,…,θN}\{\theta_1, \theta_2, \dots, \theta_N\}{θ1​,θ2​,…,θN​} 建立导向矢量，则可得到阵列流形矩阵：

$$A=[a({\theta }_1),a({\theta }_2),\dots ,a({\theta }_N)]$$

其中 $A$ 的尺寸为 $M\times N$，通常 $N\gg M$，因此 $A$ 被称为 超完备字典。

2.3 自相关矩阵

阵列接收信号的空间协方差矩阵定义为：

$$R_{xx}=E[x(t)x^H(t)]$$

在有限快拍 $L$ 的情况下，可以用样本协方差近似：

$${\hat{R}}_{xx}=\frac{1}{L}\sum _{l=1}^{L}x(l)x^H(l)$$

该矩阵包含了入射信号的空间特性，是 DOA 估计算法的核心输入之一。

三、稀疏恢复思想与 FOCUSS 算法

3.1 稀疏恢复思想

在 DOA 问题中，阵列流形矩阵 $A\in {\mathbb{C}}^{M\times N}$ 通常包含大量扫描角度的导向矢量，其中 $N\gg M$。然而，实际存在的信号源数量 $P$ 远小于 $N$，这意味着在角度域中，信号具有 稀疏性。

因此，可以将 DOA 估计转化为如下稀疏表示问题：

$$u\approx As$$

其中：

• $u\in {\mathbb{C}}^{M\times 1}$ 为观测向量（例如由协方差矩阵特征分解得到的特征向量）；
• $A$ 为超完备字典；
• $s\in {\mathbb{C}}^{N\times 1}$ 为稀疏系数向量，其非零位置对应实际来波方向。

优化目标为：

$$\underset{s}{\min }\Vert s{\Vert }_{p}s.t.u\approx As$$

其中 $0<p\le 1$，通常选用 $p$ 接近 0 来增强稀疏性。

3.2 FOCUSS 算法推导

FOCUSS（Focal Underdetermined System Solver）算法的核心思想是：通过 迭代加权最小二乘，逐步增强解的稀疏性。其更新公式为：

$$s^{(k+1)}=W^{(k)}{A}^H(A{W}^{(k)}{A}^H+\lambda I{)}^{-1}u$$

其中：

• $W^{(k)}$ 为权值矩阵；
• $\lambda$ 为正则化因子，用于提高数值稳定性；
• $k$ 为迭代次数。

权值矩阵的定义为：

$$W^{(k)}=\text{diag}\left({\left|s^{(k)}\right|}^{1-\frac{p}{2}}\right)$$

初始解 $s^{(0)}$ 通常由伪逆解给出：

$$s^{(0)}=A^{†}u$$

3.3 收敛条件与终止准则

在迭代过程中，当满足如下条件时停止迭代：

$$\frac{\Vert s^{(k+1)}-s^{(k)}{\Vert }_2}{\Vert s^{(k)}{\Vert }_2}<ϵ$$

其中 $ϵ$ 为迭代误差阈值（如 $10^{-4}$）。

最终得到的 $s$ 向量中，幅度较大的分量对应实际来波方向，其模平方即可作为 DOA 功率谱估计。

3.4 参数影响

• 稀疏因子 $p$：控制结果的稀疏性，$p$ 越小，解越稀疏；
• 正则化参数 $\lambda$：过大时解趋近于零，过小时可能数值不稳定；
• 误差阈值 $ϵ$：决定迭代精度与计算开销。

四、PINV 方法对比

在介绍 FOCUSS 算法之前，我们先来看一个更为直接的基线方法——伪逆（Pseudo-Inverse, PINV）。

4.1 PINV 解的推导

对于观测向量 $u$ 和阵列流形矩阵 $A$，我们希望找到系数向量 $s$ 使得：

$$u\approx As$$

其最小二乘解可由伪逆给出：

$$s_{\text{pinv}}=A^{†}u$$

其中 $A^{†}$ 表示 $A$ 的摩尔–彭若斯伪逆（Moore–Penrose Pseudo-Inverse）。

在 MATLAB 中，可以直接利用 pinv 函数实现：

s_pinv = pinv(A) * u;

4.2 PINV 方法的特点

• 优点：实现简单、计算速度快，常作为基准方法使用。

• 缺点：没有稀疏性约束，容易受到噪声和字典冗余的影响，导致分辨率不足。

4.3 与 FOCUSS 的差异

由此可见，PINV 方法更适合作为对比基准，而 FOCUSS 算法在高分辨率 DOA 估计中展现出明显优势。

五、MATLAB 实现与仿真

close all; clear; clc;
rng(2);

%% 参数定义区
c = physconst('LightSpeed'); % 光速 m/s
f0 = 10e9; % 载波频率 Hz
fg = 50e9; % 系统全局采样率
BW = 1e8; % 信号带宽
src_angle = deg2rad([-10,-20]); % 目标来向角 (°->rad)
M = 32; % 阵元数量
L = 128; % 快拍数
snr = 20; % 信噪比 dB
scan_angle = deg2rad(-60:0.1:60); % 扫描角度范围

%% 参数计算区
lambda = c / f0; % 载波波长 m
d = lambda / 2; % 阵元间距 m
P = length(src_angle); % 目标数

%% 回波生成
[x_sig, ~, R_sig] = echo_generate(M, d, lambda, src_angle, L, f0, BW, fg, snr);

%% 计算DOA栅格矩阵（完备字典）
i = 0:M-1;
a = exp(-1i*2*pi*d/lambda.*i'.*sin(scan_angle));

%% 计算Rxx的主特征向量
[U,~] = eig(R_sig);
u = U(:,end);

%% FOCUSS法
lamda_reg = 1e-4; % 正则化因子
lamda_err = 1e-4; % 迭代结束误差
lamda_spe = 0; % 稀疏因子
P_focuss = DOA_FOCUSS(a, u, lamda_spe, lamda_reg, lamda_err);

%% 伪逆法
P_pinv = DOA_PINV(u, a);

%% 绘图
figure;
plot(rad2deg(scan_angle), 10*log10(P_focuss), 'r-', 'LineWidth',1.5); hold on;
plot(rad2deg(scan_angle), 10*log10(P_pinv), 'b-', 'LineWidth',1.5);
xlabel('扫描角度 (°)'); ylabel('归一化功率 (dB)');
title('DOA 估计：FOCUSS vs PINV');
legend('FOCUSS','PINV');
grid on;

%% ================= 函数定义 =================
function [x, sigma, Rxx] = echo_generate(M, d, lambda, src_angle, L, f0, BW, fg, snr)
% 生成ULA接收的快拍数据和自相关矩阵
c = physconst('LightSpeed');
t = (0:L-1)/fg;
P = length(src_angle);
% 信号生成
x = zeros(M,L);
for p = 1:P
f_sig = f0 + (rand-0.5)*BW; % 随机选频带内信号
s = exp(1i*2*pi*f_sig*t); % 基带信号
a = exp(-1i*2*pi*d/lambda*(0:M-1)'*sin(src_angle(p))); % 导向矢量
x = x + a*s;
end
% 加噪声
x = awgn(x, snr, 'measured');
% 自相关矩阵
Rxx = (x*x')/L;
sigma = diag(Rxx);
end

function P_focuss = DOA_FOCUSS(scan_a, u, lamda_spe, lamda_reg, lamda_err)
% 基于FOCUSS的稀疏恢复DOA估计
rec_len = size(u,1);
Dg = scan_a;
s0 = Dg' / (Dg * Dg') * u;
for i = 1:1000
W = diag(s0.^(1 - lamda_spe/2));
s = W*W'*Dg' / (Dg*(W*W')*Dg' + lamda_reg*eye(rec_len)) * u;
if norm(s-s0,2)/norm(s0,2) < lamda_err
break;
end
s0 = s;
end
P_focuss = abs(s);
P_focuss = P_focuss ./ max(P_focuss);
P_focuss(P_focuss < 1e-4) = 1e-4;
end

function s_pinv = DOA_PINV(u, scan_a)
% 基于伪逆的稀疏恢复DOA估计
s_pinv = u' * pinv(scan_a)';
s_pinv = abs(s_pinv);
s_pinv = s_pinv ./ max(s_pinv);
end

六、总结

本文围绕 FOCUSS 算法 在 DOA 估计中的应用展开，完整介绍了从阵列信号模型、稀疏恢复思想，到算法公式推导与 MATLAB 仿真实现的全过程。通过与 PINV 方法 的对比，可以得到以下结论：

1. 稀疏恢复优势明显：FOCUSS 算法利用稀疏约束，使得角度域功率谱更为尖锐，具备更高的分辨率。
2. 抗噪声性能更好：相比于简单的伪逆方法，FOCUSS 能有效抑制噪声和旁瓣，适合在低信噪比条件下应用。
3. 计算复杂度适中：虽然需要迭代计算，但在常规阵列规模下仍具备可行性。