简介面向无线通信、雷达探测与声纳系统的实际定位与测向需求这份MATLAB脚本面向有一定信号处理基础的读者演示了最大似然估计在波达角估计中的完整应用。脚本的核心是通过多天线接收信号间的相位差建立波达角与阵列传播模型之间的数学关系再利用最大似然准则寻找最可能产生观测数据的信号到达方向。压缩包内仅含一个.m文件大小约2KB代码量小但步骤完整覆盖了数据预处理与去噪、模型设定、似然函数构建、基于MATLAB优化函数的迭代求解以及均方误差评估等环节。已有184人学习下载适合正在学习阵列信号处理、最大似然理论或波达方向估计的初学者与工程师参考。脚本将参数估计转化为最优化问题类似梯度下降或牛顿法的迭代思路清晰便于读者结合不同信噪比条件下的误差分析快速理解最大似然估计波达角估计的原理与实现细节对算法验证和课程设计都有直接帮助。1. 最大似然估计在 MATLAB 里解决什么问题最大似然估计MLE是阵列测向里最能扛极端条件的方法。拿到均匀线阵的接收数据第一反应通常是 Capon 波束形成或 MUSIC 谱估计二者在信噪比高、快拍多、信源不相关时表现不错可快拍数掉到几十、信源相干、SNR 只有几 dB 时伪峰和子空间泄漏就会把结果带偏。MLE 的出发点完全不同把估计问题写成概率模型找一组角度让观测数据出现概率最大。它的渐近方差逼近 Cramér-Rao 下界小快拍和相干信源下也扛得住代价是计算量高。下面这套 Maximum_likelyhood_estimation.m 风格的实现从似然推导、MATLAB 仿真、网格搜索到结果校验跑通一条最大似然角估计闭环。2. 最大似然角估计的原理把观测写成概率模型2.1 阵列观测模型与复高斯噪声假设先固定模型。M 元均匀线阵ULA阵元间距 d来波方向 θ 以阵列法线为 0°K 个远场窄带信号同时入射接收机在 N 个快拍内采样观测矩阵 X 是 M×N 复矩阵X A(θ) S NA(θ) 是 M×K 方向矩阵第 k 列是导向矢量 a(θ_k)第 m 个分量为 e^(j 2π (m-1) d sinθ_k / λ)。噪声 N 假设为独立同分布的复圆高斯白噪声每个元素实部虚部独立方差各为 σ²/2。这个假设是整个最大似然估计的起点噪声分布写错后面的似然函数全错。MATLAB 里导向矢量通常写成归一化间距 d_lambda d/λ 的形式这样不用显式写波长和 ULA 波数图里的横轴习惯一致function a steervec(theta_deg, M, d_lambda) % 输入角度(度)、阵元数、归一化间距输出 Mx1 导向矢量 idx (0:M-1).; a exp(1j * 2 * pi * d_lambda * idx * sind(theta_deg)); end这里用 sind 而不是 sin因为角度以度为单位传入sind 内部做了一次度到弧度换算避免每次调用都写 deg2rad。d_lambda0.5 对应半波长布阵是 DOA 仿真默认值超过 0.5 会在可见区内出现栅瓣。2.2 对数似然函数与“浓缩”技巧把 S 当作确定性未知参数时观测 X 的条件概率密度是复高斯形式对数似然可以写成log p(X|θ,S,σ²) -MN log(πσ²) - (1/σ²) ‖X - AS‖²_F直接对这个函数做极大化要同时优化 K 个角度、K×N 个信号和噪声方差变量太多。标准做法是“浓缩”先固定 θ对 S 求最小二乘解 Ŝ A†XA† 是伪逆代回原式消掉 S。整理之后最大化对数似然等价于最大化J(θ) Tr{ P_A(θ) R̂ }其中 P_A A (A^H A)^(-1) A^H 是 A 列空间的投影矩阵R̂ XX^H/N 是样本协方差矩阵。这个迹的物理含义是把接收数据的能量投影到导向矢量张成的信号子空间上投影能量越大这组 θ 越可能解释观测。这个形式是最大似然角估计的核心实现目标。它和 MUSIC 的本质区别在于MUSIC 需要先估计噪声子空间信源相干时噪声子空间估计失真而 MLE 只依赖信号子空间的投影不要求信源不相关所以相干信源场景下 MLE 依然可用。2.3 确定性 MLE 和随机 MLE 怎么选信号模型还有一个分支如果假设 S 本身也是复高斯随机过程得到的是随机 MLE。两者代价函数和实现成本差别很大类型信号模型代价函数实现成本确定性 MLES 为未知常量Tr{P_A R̂}直接搜索角度网格搜索 局部优化随机 MLES 服从复高斯分布含 log det(A R_s A^H σ²I)需要迭代常用 EM 或牛顿法工程上绝大多数情况先做确定性 MLE。它不需要估计信号的统计特性N 足够大时两种 MLE 渐近等价只有很在意小样本性能时才值得上随机 MLE而且那时多半要用 matlab 优化工具箱里的 fmincon 做带约束迭代收敛过程需要额外调试。3. 用 MATLAB 写第一个最大似然角估计脚本3.1 先造一份能复现的仿真数据不用现成工具箱直接生成两份已知真值的数据方便后面核对结果rng(2024); M 8; N 200; d_lambda 0.5; theta_true [10, -20]; K numel(theta_true); A steervec(theta_true, M, d_lambda); S (randn(K, N) 1j*randn(K, N)) / sqrt(2); % 每个信源信号功率为 1 X A * S sqrt(0.1) * (randn(M, N) 1j*randn(M, N)) / sqrt(2); % 噪声功率 0.1信号项系数除以 sqrt(2) 是为了让 randn 生成的实部和虚部合成后功率为 1对应单信源每阵元 SNR 10 dB信号功率 1噪声功率 0.1。这里不用 awgn 函数是因为 awgn 按整个矩阵功率定标两个信源叠加时功率控制不如直接乘系数直观。3.2 网格搜索最大似然谱并画图有了 X先算样本协方差再在 -90° 到 90° 的网格上扫单目标代价函数Rxx (X * X) / N; theta_grid -90:0.1:90; spec zeros(size(theta_grid)); for i 1:numel(theta_grid) a steervec(theta_grid(i), M, d_lambda); P_a (a * a) / (a * a); % 导向矢量的投影矩阵 spec(i) real(trace(P_a * Rxx)); end [~, idx] max(spec); theta_coarse theta_grid(idx); plot(theta_grid, spec); % matlab 画最大似然角谱 xlabel(角度 (deg)); ylabel(Tr(P_A Rxx)); title(最大似然角谱);P_a 是单列导向矢量对应的投影矩阵归一化分母 a*a 正好是 M。每次循环构造一次 M×M 投影矩阵再做迹1801 个网格点大约几百毫秒作为验证够用要提速可以把 a 的共轭转置提前算好或者改成矩阵化批量计算。画出来的谱在两个真实角度附近会出现尖峰但因为是单目标模型去拟合双目标数据峰位有偏移这是预期行为不是 bug。3.3 用 fminbnd 做角度精化网格步长 0.1° 决定了估计分辨率想把精度推到 0.01° 级别就要把步长降到 0.01°网格点翻十倍。更聪明的做法是先粗扫再局部精化fun (th) -real(trace(((steervec(th, M, d_lambda) * steervec(th, M, d_lambda)) ... / (steervec(th, M, d_lambda) * steervec(th, M, d_lambda))) * Rxx)); theta_est fminbnd(fun, theta_coarse - 0.5, theta_coarse 0.5);匿名函数里重复调用了三次 steervec代码可读性好但速度慢精化阶段总共只跑几十次函数求值这点开销无所谓。注意 fminbnd 求的是最小值所以代价函数加了负号。搜索区间取网格峰值左右各 0.5°前提是粗网格步长不能大于 0.5°否则可能漏掉真正的全局峰。网格步长和精化方式的对应关系如下网格步长分辨率上限1801 点耗时量级适用场景1°1°几十 ms粗定位0.1°0.1°几百 ms配合 fminbnd 精化0.01°0.01°秒级不推荐单独使用对 K1 的多目标情形这个单目标精化不适用。常见做法是交替投影固定其他角度只扫一个迭代两三轮后收敛到局部极值或者先用 MUSIC 给出初值再用 fminsearch 在高维空间细化。多目标时初始值选不好会收敛到伪峰这是最大似然估计多维搜索绕不开的问题。4. 决定最大似然估计精度的参数与 Cramér-Rao 下界4.1 四个关键参数怎么定同样一套最大似然角估计代码换参数后误差可以差两个数量级。设计仿真时优先盯这四个量参数常用范围对结果的影响M 阵元数832方差按 M³ 下降收益最大但互耦和孔径成本上升N 快拍数1001000方差按 1/N 下降N 太小时协方差矩阵条件数恶化d/λ 阵元间距0.5超过 0.5 出现栅瓣网格搜索可能锁到镜像角度SNR020 dB低于门限时估计误差不是平滑增大而是跳变阵元数对精度的影响最明显。CRB 里 M 以三次方出现8 元改 16 元理论误差直接缩到八分之一左右所以阵列设计阶段优先加阵元而不是加快拍。快拍数提升是线性的100 快拍换 1000 快拍只换回 3.2 倍的方差改善。N 太小时样本协方差 R̂ 接近奇异投影矩阵计算在数值上很容易出问题常见做法是对 R̂ 做对角加载即 R̂ εIε 取迹的 0.01 倍左右。4.2 把 CRB 写进 MATLAB 校验理论极限Cramér-Rao 下界是判断实现有没有写对的重要参照。单个远场信源、复高斯噪声、高 SNR 近似下角度估计的 CRB 有闭式解function crb crb_ula(theta_deg, M, d_lambda, N, snr_db) % 单目标、高 SNR 近似下的角度估计 CRB返回弧度平方 snr 10^(snr_db / 10); th deg2rad(theta_deg); crb (6 / (N * M * (M^2 - 1) * snr)) ... * (1 / (2 * pi * d_lambda * cos(th)))^2; end std_theta sqrt(crb_ula(10, 8, 0.5, 200, 10)) * 180 / pi;这个组合跑出来大约 0.05°。注意公式里 cosθ 在分母上角度接近 ±90° 时 CRB 发散所以端射方向的估计误差天然就大。另外这是高 SNR 近似实际 CRB 在低 SNR 时要比公式高一点蒙特卡洛结果落在公式值 11.5 倍之间都算正常。做阵列设计时可以直接用这个函数画“误差 vs 角度”曲线看哪个扇区满足指标。5. 验证最大似然估计的蒙特卡洛技巧5.1 用 RMSE 和 CRB 对比判断实现有没有写错单次仿真跑出 0.04° 的误差没有意义可能是运气好。把数据生成挪进循环跑 500 次蒙特卡洛统计均方根误差nMc 500; errs zeros(1, nMc); M 8; N 200; d_lambda 0.5; theta_true 10; K 1; % 先只验单目标避免多峰混淆 A steervec(theta_true, M, d_lambda); for mc 1:nMc S (randn(K, N) 1j*randn(K, N)) / sqrt(2); X A * S sqrt(0.1) * (randn(M, N) 1j*randn(M, N)) / sqrt(2); % 这里插入 3.2 和 3.3 的网格搜索与 fminbnd 代码得到 theta_est errs(mc) theta_est - theta_true; end rmse sqrt(mean(errs.^2));校验规则RMSE 与 CRB 平方根转成度的比值在 11.5 之间说明似然函数、导向矢量和优化过程都写对了比值大于 3先怀疑粗网格步长太大导致 fminbnd 初始区间没包住峰值把 0.1° 步长改成 0.2° 再试确认问题方向。还有一个常被忽略的坑角度估计在 ±90° 边界附近fminbnd 可能把峰值推出定义域精化之前先对网格峰值做一次越界钳制clamp 到 [-89.9, 89.9]。最后一个小技巧是调试时把粗网格步长放宽到 0.5°先 plot 出来看谱峰形状确认单峰再加密比一上来就跑 0.01° 步长省两个数量级的时间。本文还有配套的精品资源点击获取