2. 常规波束形成方位谱的栅格与插值误差分析

问题定位：峰值搜索的“台阶效应”

常规波束形成（CBF）的方位谱是连续角度 θ 的函数，但工程实现中只能在离散栅格上计算：P(θ_k)，θ_k = -90° + k·Δθ。峰值搜索返回的是栅格点，而真实来向几乎不会恰好落在栅格上，于是产生栅格失配误差，其上限约为 ±Δθ/2。本节用同一组阵列数据，对比粗栅格、细栅格与三点抛物线插值三种峰值定位方式，量化这一误差。

数据与基线设置

沿用单频窄带远场模型：均匀直线阵 M=16，半波长间距 d=λ/2，载频 f0=1 kHz，采样率 fs=8 kHz，快拍数 N=512，单目标真实来向设为 θ0 = 12.3°（刻意取非整栅格值，以暴露失配）。噪声为复高斯白噪声，信噪比约 10 dB。CBF 方位谱按导向矢量扫描计算：

P(θ) = a(θ)^H R a(θ),  a(θ) = [1, e^{-jπ sinθ}, ..., e^{-jπ(M-1) sinθ}]^T

其中 R = (1/N)·X·X^H 为样本协方差矩阵。注意导向矢量按 d=λ/2 化简后相位项为 e^{-jπ m sinθ}，与频率无关，这是半波长阵的便利之处。

三种峰值定位方式的实现

核心差异只在“如何从离散谱值得到角度估计”，谱计算部分完全共用：

import numpy as np

M, N, f0, fs = 16, 512, 1000.0, 8000.0
theta0 = 12.3                      # 真实来向（度）
rng = np.random.default_rng(7)

def steering(theta_deg):
    th = np.deg2rad(np.asarray(theta_deg, dtype=float))
    m = np.arange(M)
    return np.exp(-1j * np.pi * m[:, None] * np.sin(th)[None, :])  # M x K

# 生成单目标阵列数据
t = np.arange(N) / fs
s = np.exp(1j * 2 * np.pi * f0 * t)
A = steering([theta0])                       # M x 1
X = A @ s[None, :] + 0.3 * (rng.standard_normal((M, N))
                            + 1j * rng.standard_normal((M, N))) / np.sqrt(2)
R = X @ X.conj().T / N

def spectrum(grid_deg):
    A = steering(grid_deg)
    return np.real(np.einsum('mk,mn,nk->k', A.conj(), R, A))

def peak_coarse(grid_deg):
    P = spectrum(grid_deg)
    k = int(np.argmax(P))
    return grid_deg[k]

def peak_interp(grid_deg):
    P = spectrum(grid_deg)
    k = int(np.argmax(P))
    if k == 0 or k == len(grid_deg) - 1:
        return grid_deg[k]
    y0, y1, y2 = P[k-1], P[k], P[k+1]
    denom = (y0 - 2*y1 + y2)
    delta = 0.5 * (y0 - y2) / denom if denom != 0 else 0.0
    return grid_deg[k] + delta * (grid_deg[1] - grid_deg[0])

grids = {
    '粗栅格 2.0°': np.arange(-90, 90 + 1e-9, 2.0),
    '中栅格 0.5°': np.arange(-90, 90 + 1e-9, 0.5),
    '细栅格 0.1°': np.arange(-90, 90 + 1e-9, 0.1),
}
for name, g in grids.items():
    est = peak_coarse(g)
    print(f"{name:12s} 峰值={est:7.2f}°  偏差={est - theta0:+6.2f}°")

g = grids['中栅格 0.5°']
print(f"插值(0.5°栅格) 峰值={peak_interp(g):7.2f}°  偏差={peak_interp(g) - theta0:+6.2f}°")

三点抛物线插值假设主瓣顶部近似抛物线，用峰值点及其左右邻点拟合顶点，偏移量 δ = 0.5·(y₋₁ − y₊₁)/(y₋₁ − 2y₀ + y₊₁)，单位为栅格间隔。它把有效分辨率从 Δθ 提升到远小于 Δθ 的量级，代价仅是三次谱值运算。

栅格失配误差的量级与判读

下表为按上述设置运行得到的示例结果（随机种子固定，属可复现的教学示例，非通用实测结论）：

判读要点：粗栅格的偏差被限制在半个栅格内，但已接近 0.3°，对后续波束域处理或跟踪滤波而言不可忽略；细栅格把偏差压到 0.05° 以内，但计算量随栅格数线性增长（0.1° 栅格需 1801 个扫描点，是 2° 栅格的 20 倍）；插值在 0.5° 栅格上就达到与 0.1° 细栅格相当的精度，是精度与计算量的折中点。

栅格失配误差的抑制手段

• 细化栅格：最直接，误差上限随 Δθ 线性下降，但计算量同比上升，适合离线或小规模阵。
• 峰值插值：抛物线（三点）或高斯拟合，几乎零额外开销，是工程默认首选。
• 两级搜索：先粗栅格定位主瓣，再在峰值邻域 ±Δθ 内以 0.01° 细扫，兼顾速度与精度。
• 解析/迭代优化：对 P(θ) 在峰值附近做梯度上升或牛顿迭代，可逼近连续谱的真实极值。

验证与通过条件

复现时按以下条件判断实现是否正确：

1. 三种栅格的 argmax 估计值必须落在各自栅格点上，且偏差绝对值不超过 Δθ/2；
2. 细栅格（0.1°）的偏差应显著小于粗栅格（2°），体现误差随步长线性下降；
3. 插值结果应位于峰值点与其较近邻点之间，且偏差小于所用栅格的 Δθ/2；
4. 若把 theta0 改为恰好落在栅格上的值（如 12.0°），粗栅格偏差应接近 0，说明误差来源确为栅格失配而非算法错误。

最后一条是关键的对照实验：它把“栅格失配”与“谱估计本身偏差”区分开，避免把插值或细栅格带来的改善误判为算法精度提升。