3. 最小方差无失真响应（MVDR）波束形成

任务说明

本案例实现最小方差无失真响应（MVDR）波束形成器，在存在方向性干扰的声场环境中增强期望方向信号，并与延迟-求和（DSB）波束形成器进行输出信干比对比。输入为均匀线阵的仿真接收数据，输出为两种波束形成器的加权向量及波束输出信号，并给出可量化的性能指标。

信号模型与算法原理

考虑\(M\)元均匀线阵，阵元间距为半波长。期望信号从\(\theta_0\)方向入射，干扰从\(\theta_1\)方向入射。第\(m\)个阵元在时刻\(n\)的接收信号为：

\[ x_m(n) = s(n - \tau_m(\theta_0)) + j(n - \tau_m(\theta_1)) + v_m(n) \]

其中\(\tau_m(\theta) = (m-1)d\sin\theta / c\)为相对参考阵元的时延，\(v_m(n)\)为独立高斯白噪声。MVDR波束形成器在保证期望方向无失真响应的约束下最小化输出功率：

\[ \min_{\mathbf{w}} \mathbf{w}^H \mathbf{R}_{xx} \mathbf{w} \quad \text{s.t.} \quad \mathbf{w}^H \mathbf{a}(\theta_0) = 1 \]

其闭式解为\(\mathbf{w}_{\text{MVDR}} = \frac{\mathbf{R}_{xx}^{-1}\mathbf{a}(\theta_0)}{\mathbf{a}^H(\theta_0)\mathbf{R}_{xx}^{-1}\mathbf{a}(\theta_0)}\)。延迟-求和波束形成器则直接取\(\mathbf{w}_{\text{DSB}} = \mathbf{a}(\theta_0)/M\)。

代码实现

import numpy as np

def generate_array_data(theta_s, theta_i, M=8, N=2000, snr=10, inr=20, seed=42):
    """生成均匀线阵接收数据（半波长间距）"""
    rng = np.random.default_rng(seed)
    d_lambda = 0.5
    n = np.arange(N)
    # 期望信号与干扰的复包络
    s = np.exp(1j * 0.1 * n)                     # 期望信号
    j = np.exp(1j * 0.3 * n) * (10 ** (inr / 20)) # 干扰（功率高于信号）
    # 导向向量
    def steer(theta_deg):
        theta = np.deg2rad(theta_deg)
        return np.exp(-1j * 2 * np.pi * d_lambda * np.arange(M) * np.sin(theta))
    a_s = steer(theta_s).reshape(-1, 1)
    a_i = steer(theta_i).reshape(-1, 1)
    # 接收数据矩阵 X: M x N
    X = a_s @ s.reshape(1, -1) + a_i @ j.reshape(1, -1)
    # 加性噪声（功率归一化到1，SNR由信号幅度控制）
    noise_power = 10 ** (-snr / 10)
    noise = np.sqrt(noise_power / 2) * (rng.standard_normal((M, N)) + 1j * rng.standard_normal((M, N)))
    return X + noise, a_s.flatten()

def mvdr_weights(Rxx, a_s):
    """MVDR波束形成器加权向量"""
    Rinv = np.linalg.inv(Rxx)
    w = Rinv @ a_s / (a_s.conj() @ Rinv @ a_s)
    return w

def dsb_weights(a_s, M):
    """延迟-求和波束形成器加权向量"""
    return a_s / M

def beamformer_output(w, X):
    """计算波束输出信号"""
    return w.conj() @ X

def sinr_estimate(y, s_ref, theta_s, theta_i, M, d_lambda=0.5):
    """基于输出信号与参考信号相关性估计SINR（简化方法）"""
    # 期望信号分量估计（与参考信号相关部分）
    corr = np.sum(y * np.conj(s_ref)) / len(s_ref)
    s_power = np.abs(corr) ** 2
    total_power = np.mean(np.abs(y) ** 2)
    return 10 * np.log10(s_power / (total_power - s_power + 1e-12))

# ========== 参数设置 ==========
M = 8
theta_s, theta_i = 10.0, -30.0
N = 2000
snr, inr = 10, 20

# 生成数据
X, a_s = generate_array_data(theta_s, theta_i, M, N, snr, inr)
# 采样协方差矩阵
Rxx = (X @ X.conj().T) / N

# 计算两种波束形成器权重
w_mvdr = mvdr_weights(Rxx, a_s)
w_dsb = dsb_weights(a_s, M)

# 输出信号
y_mvdr = beamformer_output(w_mvdr, X)
y_dsb = beamformer_output(w_dsb, X)

# 参考信号（用于SINR估计）
n_ref = np.arange(N)
s_ref = np.exp(1j * 0.1 * n_ref)

# 性能对比
sinr_mvdr = sinr_estimate(y_mvdr, s_ref, theta_s, theta_i, M)
sinr_dsb = sinr_estimate(y_dsb, s_ref, theta_s, theta_i, M)

print(f"MVDR 输出SINR: {sinr_mvdr:.2f} dB")
print(f"DSB  输出SINR: {sinr_dsb:.2f} dB")
print(f"MVDR 相对DSB提升: {sinr_mvdr - sinr_dsb:.2f} dB")

运行与验证

依赖：Python 3.8+，NumPy库。

运行方式：将代码保存为mvdr_beamformer.py，在终端执行python mvdr_beamformer.py。

验证方法：程序输出两种波束形成器的输出SINR估计值。判定标准：

• MVDR的输出SINR应显著高于DSB（通常提升10dB以上，因MVDR在干扰方向形成零陷）；
• DSB输出SINR接近输入SNR（10dB）减去干扰泄漏损失；
• MVDR输出SINR应接近理论上限（由阵元数M和输入SNR决定，此处约17-20dB）。

若需验证波束方向图，可额外计算20*log10(abs(w.conj() @ steer(theta)))在-90°~90°扫描，观察MVDR在-30°处形成深零陷而DSB无此特性。