4. 基于协方差矩阵的MVDR方位估计实现

协方差矩阵估计与MVDR谱计算链路

MVDR（Capon）方位谱的核心是：对每个扫描角度 θ，用导向矢量 a(θ) 对样本协方差矩阵 R̂ 做一次“白化—投影”，得到该方向的输出功率。有限快拍下 R̂ 秩亏或接近奇异，直接求逆会放大噪声，因此必须引入对角加载。本节用一段可运行的 Python 实现完整链路，并比较不同快拍数下的分辨表现。

实现要点与关键代码

下面给出一个自洽的最小实现：构造两源 ULA 接收数据，估计 R̂，做对角加载，用 Cholesky 分解求逆，再在角度栅格上计算 MVDR 谱。所有参数在本例中定义，可直接运行。

import numpy as np

def ula_steering(N, d_lambda, theta_deg):
    """均匀直线阵导向矢量，theta 为与阵列法线夹角（度）"""
    m = np.arange(N)
    return np.exp(-1j * 2 * np.pi * d_lambda * m * np.sin(np.deg2rad(theta_deg)))

def estimate_cov(X):
    """X: N×K，返回样本协方差 R̂ = X Xᴴ / K"""
    K = X.shape[1]
    return X @ X.conj().T / K

def mvdr_spectrum(R, N, d_lambda, grid_deg, eps=0.05):
    """对角加载 + Cholesky 求逆 + MVDR 谱"""
    delta = eps * np.trace(R).real / N          # 加载量
    Rl = R + delta * np.eye(N, dtype=complex)
    L = np.linalg.cholesky(Rl)                  # Rl = L Lᴴ
    Rl_inv = np.linalg.solve(L.conj().T, np.linalg.solve(L, np.eye(N)))
    P = np.empty(len(grid_deg))
    for i, th in enumerate(grid_deg):
        a = ula_steering(N, d_lambda, th)
        P[i] = 1.0 / np.real(a.conj() @ Rl_inv @ a)
    return P, delta

# ---- 教学示例参数 ----
np.random.seed(0)
N, d_lambda = 8, 0.5
src_deg = [-10.0, 8.0]          # 两源来向
snr_db = 10.0
grid = np.arange(-60, 60.01, 0.5)

def make_data(K):
    A = np.column_stack([ula_steering(N, d_lambda, s) for s in src_deg])
    S = (np.random.randn(2, K) + 1j*np.random.randn(2, K)) / np.sqrt(2)
    S *= 10**(snr_db/20)
    noise = (np.random.randn(N, K) + 1j*np.random.randn(N, K)) / np.sqrt(2)
    return A @ S + noise

for K in (32, 128, 512):
    X = make_data(K)
    R = estimate_cov(X)
    P, delta = mvdr_spectrum(R, N, d_lambda, grid, eps=0.05)
    peaks = grid[np.argsort(P)[-2:]]
    print(f"K={K:4d}  δ={delta:.4f}  谱峰≈{np.sort(peaks)}")

结果判读与快拍数对比

运行后应观察到：K=32 时两源谱峰可能合并或位置偏移，K=128 起两峰清晰分离，K=512 时峰位稳定在 −10° 与 8° 附近。判断依据是谱峰数量、峰位偏差（相对真实来向）与两峰之间的谷深。若 K 很小且 δ 过小，会出现多个虚假尖峰，这是 R̂ 病态导致求逆放大的典型现象。