1. 均匀直线阵单频信号方位估计基线实现
基线任务与信号模型
本案例建立一条最小可用的方位估计链路:均匀直线阵(ULA)接收单频窄带远场信号,用常规波束形成(CBF)在角度域扫描,取方位谱峰值作为来向估计。后续案例的栅格误差、导向矢量校验都复用这里的阵列参数与函数结构。
设定阵列沿 x 轴排布,M 个阵元间距 d,信号波长 λ = c/f,来向 θ 以阵列法线(y 轴)为 0°,逆时针为正。相邻阵元相对参考阵元的时延对应相位差:
φ = 2π·d·sin(θ)/λ
导向矢量(steering vector)为 a(θ) = [1, e^{-jφ}, e^{-j2φ}, …, e^{-j(M-1)φ}]^T。CBF 方位谱为各方向导向矢量与接收数据的匹配输出功率:
P(θ) = |a(θ)^H · x|² / M²
本案例取 f = 1 kHz、c = 1500 m/s(λ = 1.5 m)、M = 8、d = λ/2 = 0.75 m、真实来向 θ₀ = 20°,单快拍加复高斯噪声,信噪比约 20 dB。以下为完整可运行实现(Python + NumPy)。
可复用基线实现
import numpy as np
C = 1500.0 # 声速 m/s
F = 1000.0 # 信号频率 Hz
LAM = C / F # 波长 1.5 m
M = 8 # 阵元数
D = LAM / 2 # 阵元间距 0.75 m
THETA0 = 20.0 # 真实来向(度)
def steering(theta_deg, m=M, d=D, lam=LAM):
"""ULA 导向矢量,theta 以阵列法线为 0°,逆时针为正"""
theta = np.deg2rad(np.asarray(theta_deg, dtype=float))
idx = np.arange(m)
phi = 2 * np.pi * d * np.sin(theta) / lam
return np.exp(-1j * idx[..., None] * phi) if theta.ndim else np.exp(-1j * idx * phi)
def cbf_spectrum(x, thetas_deg, m=M, d=D, lam=LAM):
"""常规波束形成方位谱,x 为 M×1 接收数据"""
A = np.stack([steering(t, m, d, lam) for t in np.atleast_1d(thetas_deg)], axis=1)
return np.abs(A.conj().T @ x).ravel() ** 2 / m ** 2
def estimate_doa(x, thetas_deg):
P = cbf_spectrum(x, thetas_deg)
k = int(np.argmax(P))
return thetas_deg[k], P
# 生成单快拍接收数据
rng = np.random.default_rng(2024)
a0 = steering(THETA0)
s = np.exp(1j * rng.uniform(0, 2 * np.pi)) # 单位功率单频复信号
snr_db = 20.0
noise_power = 10 ** (-snr_db / 10)
n = np.sqrt(noise_power / 2) * (rng.standard_normal(M) + 1j * rng.standard_normal(M))
x = a0 * s + n
thetas = np.arange(-90.0, 90.0 + 0.1, 0.1) # 0.1° 粗栅格
theta_hat, P = estimate_doa(x, thetas)
print(f"估计来向: {theta_hat:.2f}°, 真实来向: {THETA0:.2f}°")
print(f"峰值功率: {P.max():.4f}")
运行后应看到估计值落在 20° 附近(示例随机种子下偏差通常在零点几度量级,具体数值取决于噪声实现)。若偏差明显超过 1°,先检查导向矢量的相位符号与 sin 定义是否一致,而不是急着换算法。
阵元数与间距对主瓣宽度的影响
主瓣宽度决定基线流程的角度分辨能力。对 ULA,3 dB 主瓣宽度近似为 Δθ ≈ 0.886·λ/(M·d·cosθ₀)(弧度)。保持 d = λ/2 时,阵元数翻倍主瓣宽度约减半;保持 M·d 不变时主瓣宽度基本不变,但 d > λ/2 会出现栅瓣,导致峰值搜索可能锁到错误方向。下面用同一函数扫描三组参数验证这一关系。
def mainlobe_width_deg(m, d, lam=LAM, theta0=THETA0):
return np.rad2deg(0.886 * lam / (m * d * np.cos(np.deg2rad(theta0))))
for m, d in [(8, LAM/2), (16, LAM/2), (8, LAM)]:
w = mainlobe_width_deg(m, d)
print(f"M={m:2d}, d={d/LAM:.2f}λ -> 理论3dB主瓣宽度 ≈ {w:.2f}°")
预期输出:M=8, d=0.50λ 约 8.5°,M=16, d=0.50λ 约 4.2°,M=8, d=1.00λ 约 4.2° 但伴随栅瓣。第三组说明:仅看主瓣宽度会误判性能,必须同时检查方位谱在 ±90° 内是否出现与主瓣等高的旁瓣峰。
d > λ/2 时,栅瓣出现条件为 sinθ_g = sinθ₀ ± k·λ/d。本例 d = λ 时在 θ ≈ -20° 附近会出现镜像栅瓣,峰值搜索可能返回错误来向。工程上通常取 d ≤ λ/2 作为无模糊基线配置。基线流程与判读关系
结果判读与基线验收
基线流程的验收不看单次估计是否“恰好等于 20°”,而看三点:其一,方位谱主峰位置与真实来向偏差在栅格步长与噪声引起的合理范围内(本例 0.1° 栅格、20 dB 信噪比下通常小于 1°);其二,主峰两侧旁瓣电平明显低于主峰,无等高峰;其三,改变 M 或 d 后主瓣宽度变化方向与理论公式一致。满足这三点,即可把 steering、cbf_spectrum、estimate_doa 三个函数作为后续案例的公共基线。
| 检查项 | 判据 | 不通过时的常见原因 |
|---|---|---|
| 峰值位置 | |θ̂ − θ₀| < 1° | 导向矢量相位符号反、sin/cos 混用 |
| 旁瓣结构 | 无与主峰等高的栅瓣 | d > λ/2 或频率设置与波长不一致 |
| 主瓣宽度趋势 | M 增大时 Δθ 减小 | 扫描范围或步长设置错误 |
snr_db 调到 0 dB 再运行一次,可观察峰值位置抖动与旁瓣抬升;这为后续案例讨论估计方差与栅格误差提供了对照基线。