企业官网建设流程全解析

简介：面向MATLAB信号处理学习者的最大似然波达角估计实用脚本，聚焦于利用最大似然估计（MLE）算法从多天线接收信号中解析波达角，适用于无线通信、雷达探测与声纳系统等领域的学习与研究。脚本涵盖数据预处理、模型设定、似然函数构建、优化求解及误差计算等完整流程，通过相位差模型将波达角与观测数据关联，并借助fminunc等优化工具进行迭代估计，帮助读者掌握MLE在阵列信号处理中的具体实现与性能评估方法。最大似然估计作为经典参数估计方法，其在波达角问题中的代码实现具有较强的工程参考价值。包内仅有1个m文件，文件大小仅2KB，无需额外数据即可运行，适合作为入门或课程设计的参考模板。已有184人学习，对于希望快速理解最大似然估计原理及波达角估计代码结构的开发者，是一份轻量但功能完整的示例。

1. 当MUSIC失效时，最大似然估计为什么值得自己写一遍

在很多阵列信号处理项目里，MUSIC、ESPRIT这类子空间算法是默认首选，因为它们算得快、不用迭代。但一旦遇到低信噪比、少快拍或者相干信源，子空间分解会直接丢掉信号子空间维度，谱峰分裂，估计偏差非常大。这时候最大似然估计MLE的价值就体现出来了：它在模型正确的前提下是渐进最优的，方差能逼近克拉美罗界。Maximum_likelyhood_estimation.m就是这样一个用MATLAB从零实现最大似然波达角估计的完整脚本，覆盖数据生成、似然函数构造、fminunc优化和MSE评估。适合做雷达测向、声纳定位和5G阵列通信的工程师，也适合写课程设计时想弄懂MLE而不是只调工具箱的人。这篇笔记我会从阵列模型讲到集中似然，再给出可直接跑的MATLAB代码和性能验证方法，最后把我调试时踩过的优化器坑一并说清。

2. 均匀线阵模型与集中似然函数：MLE为什么能逼近克拉美罗界

最大似然估计的第一步不是写优化代码，而是把接收信号模型写到能求导的程度。这里考虑最常见的均匀线阵ULA，假设有M个阵元，阵元间距d，信号波长为lambda，有一个远场窄带信号从角度theta入射。那么第t个快拍的接收向量可以写成：

x(t) = a(theta) * s(t) + n(t)

其中a(theta)是导向矢量，s(t)是信号复幅度，n(t)是零均值复高斯白噪声。这个模型看起来简单，但它决定了后面所有推导：导向矢量里只有角度这一个未知参数，而幅度、噪声功率都算扰动项。MLE的思路就是找到一个theta，让观测数据x的联合概率密度最大。

2.1 导向矢量与阵列流型矩阵的实现

对于ULA，第m个阵元相对参考点的相位延迟是2*pi*d*(m-1)*sin(theta)/lambda，所以导向矢量写出来是一个复指数向量。工程实现时要注意：角度要转成弧度，阵元间隔以波长为单位时，代码里可以直接用sin(theta)乘以系数，不需要真的去计算每条传播路径的延迟。下面这段MATLAB代码生成单信号源的导向矢量：

function a = steering_vector(theta_deg, M, d_over_lambda) % theta_deg: 入射角，单位度 % M: 阵元数量 % d_over_lambda: 阵元间距除以波长，通常取0.5 theta = deg2rad(theta_deg); m = (0:M-1).'; a = exp(1j * 2 * pi * d_over_lambda * m * sin(theta)); end

这段函数里，m是列向量，sin(theta)是标量，两者相乘得到M维相位差向量。exp(1j * ...)生成复指数导向矢量，它是后续所有似然函数计算的基础。如果信源不止一个，就把多个导向矢量横向拼成矩阵A = [a(theta1), a(theta2), ...]，这个矩阵叫阵列流型矩阵。在最大似然估计中，我们最终要优化的是这个矩阵里包含的角度参数。

2.2 集中似然：把多维优化降成一维

直接对theta、信号幅度S和噪声功率sigma^2求联合最大似然，参数太多，优化不现实。标准做法是先用最小二乘解把线性参数消掉。对固定theta，信号幅度的最大似然解是：

s_hat = (a^H * a)^(-1) * a^H * x

对于单信源，a^H * a就是个复数值，等于M，所以a^H * x就是匹配滤波输出。把S的解析解代回原似然函数，就得到只关于theta的集中对数似然函数。忽略常数项后，最大化似然等价于最大化：

J(theta) = x^H * a * (a^H * a)^(-1) * a^H * x

这个表达式在DOA文献里常写成tr(P_A * R_xx)，其中P_A = A*(A^H*A)^(-1)*A^H是导向矢量张成空间的投影矩阵，R_xx是样本协方差矩阵。为什么要用多个快拍呢？单快拍噪声影响太大，用N个快拍平均得到协方差矩阵后，投影能量能被平滑出来。下面的代码计算负对数似然代价函数，方便后面用fminunc求最小值：

function val = neg_log_likelihood(theta_deg, X, M, d_over_lambda) % X: M x N 复数接收矩阵，N为快拍数 % 返回标量代价，优化时需要最小化该值 a = steering_vector(theta_deg, M, d_over_lambda); P_A = a * (a' * a)^(-1) * a'; % 投影矩阵 R_xx = X * X' / size(X, 2); % 样本协方差矩阵 val = -real(trace(P_A * R_xx)); % 取负号是因为fminunc找最小值 end

这里P_A是M维复方阵，trace(P_A * R_xx)计算信号在导向矢量方向上的能量。由于我们要求最小值，所以把要最大化的目标函数取负。real()是为了把浮点运算中微小的虚部误差去掉，不影响梯度方向。注意：a'是共轭转置，而不是普通转置，这是复信号处理最容易出错的地方。

下表列出了集中似然推导中几个关键参数的角色，弄混了会在优化时得到奇怪的结果：

在低快拍场景下，R_xx的估计误差会直接传导到代价函数里，导致MLE出现多个局部极值。这也是为什么后面必须做多起点搜索，而不是只调用一次fminunc。

3. 从仿真数据到fminunc求解：最大似然角估计的完整MATLAB实现

这一章直接给出能跑的脚本。设计思路是：先用已知角度生成仿真数据，再用最大似然估计把它还原出来，最后和真实值对比。这样你能清楚看到每一步数据在哪、参数怎么改。

3.1 生成带噪接收数据

仿真数据质量直接影响优化收敛。这里设置M=8个阵元，目标信号来自theta_true=10度，快拍数N=200，信噪比设为10dB。信号幅度用复高斯随机变量，噪声是独立复高斯白噪声。代码如下：

rng(0); % 固定随机种子，方便复现 M = 8; % 阵元数量 N = 200; % 快拍数 theta_true = 10; % 真实波达角（度） d_over_lambda = 0.5; % 半波长间距 snr = 10; % 信噪比dB a = steering_vector(theta_true, M, d_over_lambda); s = sqrt(10^(snr/10)) * (randn(1, N) + 1j*randn(1, N)) / sqrt(2); n = (randn(M, N) + 1j*randn(M, N)) / sqrt(2); X = a * s + n;

这里信号幅度sqrt(10^(snr/10))是把dB信噪比转成线性幅度比。噪声功率归一化为1，所以信号功率直接由这个幅度系数控制。randn(1,N)+1j*randn(1,N)得到复高斯信号，除以sqrt(2)保证实部和虚部总功率为1。X的每一列是一个快拍，每一行是一个阵元的采样。

3.2 一维网格扫描加fminunc精估计

一维角度搜索的初值比二维、多维场景容易处理。我的做法是先做粗网格扫描，用扫描结果作为fminunc的初始点。这样比随机初值稳定得多，也不会陷到远离真实角度的局部极值。下面这段代码完成网格扫描和精估计：

theta_grid = -90:0.5:90; % 粗网格 cost_grid = zeros(size(theta_grid)); for k = 1:length(theta_grid) cost_grid(k) = neg_log_likelihood(theta_grid(k), X, M, d_over_lambda); end [~, idx] = min(cost_grid); theta_init = theta_grid(idx); % 网格最优值作为初值 options = optimoptions('fminunc', ... 'Algorithm', 'quasi-newton', ... 'Display', 'off', ... 'MaxIterations', 200, ... 'OptimalityTolerance', 1e-8, ... 'StepTolerance', 1e-8); theta_hat = fminunc(@(t) neg_log_likelihood(t, X, M, d_over_lambda), ... theta_init, options); fprintf('估计角度: %.4f deg, 真实角度: %.4f deg\n', theta_hat, theta_true);

这段代码里，theta_grid步长0.5度，网格扫描的协方差矩阵重复计算了181次，对于单信源MLE完全够快。fminunc的Algorithm选quasi-newton，因为代价函数光滑，不需要提供解析梯度。MaxIterations设200次，对一维问题足够了。如果优化结果停在网格初值附近，说明网格分辨率不够，需要加密。

3.3 优化器参数怎么选

下表是我在调试这个脚本时常用的参数组合，不同值会影响收敛速度和精度，但不是越大越好：

如果你的MATLAB没有优化工具箱，也可以换成一维黄金分割搜索或fminbnd。fminbnd不需要初值，只需要角度区间，但收敛速度比fminunc慢一点。实际测试中，fminunc在半波长间距、单信源情况下基本无偏，偏差主要来自有限快拍和噪声实现。

提示：不要用fminunc默认的trust-region算法，因为它要求目标函数返回梯度值，这里我们只给了函数值，会直接报错。

4. 用MSE和克拉美罗界评估最大似然估计的性能

写完了优化流程，下一步是量化估计算法到底准不准。只看一次运行结果没有说服力，必须做蒙特卡洛统计。我会在这个章节给出MSE计算代码和克拉美罗界CRB的数值实现，方便你画性能曲线。

4.1 蒙特卡洛MSE统计

固定theta_true=10度，信噪比从0dB到20dB，每个信噪比下重复L=200次独立实验，记录每次估计结果，最后计算均方根误差RMSE。脚本主体是一个大循环，注意每次实验都要重新生成信号和噪声，不能复用同一组数据。代码如下：

L = 200; snr_list = 0:5:20; rmse_list = zeros(size(snr_list)); for snr_idx = 1:length(snr_list) snr = snr_list(snr_idx); err_sum = 0; for trial = 1:L % 重新生成数据和噪声 s = sqrt(10^(snr/10)) * (randn(1, N) + 1j*randn(1, N)) / sqrt(2); n = (randn(M, N) + 1j*randn(M, N)) / sqrt(2); X = a * s + n; % 网格初值 cost_grid = arrayfun(@(th) neg_log_likelihood(th, X, M, d_over_lambda), theta_grid); [~, idx] = min(cost_grid); theta_init = theta_grid(idx); % 精细优化 theta_hat = fminunc(@(t) neg_log_likelihood(t, X, M, d_over_lambda), ... theta_init, options); err_sum = err_sum + (theta_hat - theta_true)^2; end rmse_list(snr_idx) = sqrt(err_sum / L); end disp(table(snr_list.', rmse_list.', 'VariableNames', {'SNR_dB', 'RMSE_deg'}));

这里arrayfun遍历整个网格，比for循环简洁，但二者等价。err_sum累加的是角度误差平方，最后除以L再开方得到RMSE。当你运行这段代码时，会看到RMSE随信噪比升高而下降，但在高信噪比段会进入一个平台，这个平台通常是由网格扫描分辨率或优化容差导致的，而不是MLE本身的问题。

4.2 克拉美罗界的数值实现

最大似然估计的理论下界是CRB。单信源ULA的CRB有闭式表达式，可以直接计算。这里我给出一个数值实现，方便和蒙特卡洛结果画在同一张图上：

function crb = compute_crb(theta_deg, M, N, snr, d_over_lambda) theta = deg2rad(theta_deg); m = (0:M-1).'; % 导向矢量对角度的一阶导数 a = exp(1j * 2 * pi * d_over_lambda * m * sin(theta)); a_der = 1j * 2 * pi * d_over_lambda * m .* cos(theta) .* a; snr_linear = 10^(snr/10); % 单信源CRB公式 crb = 1 / (2 * N * snr_linear * (norm(a_der)^2 - abs(a'*a_der)^2 / M)); end

CRB公式里的分母由两部分组成：导向矢量导数的能量norm(a_der)^2，以及投影到导向矢量方向后剩下的部分。注意这里的a'*a_der是复数内积，abs()取模后平方。利用恒等式a'*a = M，可以简化成a_der在导向矢量正交方向的能量。这个值越大，CRB越小，说明阵列对该方向的角度敏感度越高。

4.3 扫信噪比时的性能图表

在matlab里通常用semilogy画RMSE和CRB曲线，纵轴用对数刻度。下面是一个快速绘图代码：

crb_list = arrayfun(@(s) compute_crb(theta_true, M, N, s, d_over_lambda), snr_list); figure; semilogy(snr_list, rmse_list, 'o-', 'LineWidth', 1.5); hold on; semilogy(snr_list, sqrt(crb_list), 'r--', 'LineWidth', 1.5); xlabel('SNR (dB)'); ylabel('RMSE / sqrt(CRB) (deg)'); legend('MLE RMSE', 'sqrt(CRB)'); grid on;

用semilogy而不直接用plot是因为RMSE在低信噪比可能是好几度，在高信噪比会小于0.1度，线性坐标下会压扁。sqrt(crb_list)是把CRB转成角度标准差，与RMSE同量纲才能对比。如果MLE结果在高信噪比区域没有贴着CRB曲线，通常问题出在初值网格没有覆盖全局最优，或者优化器在平坦区域提前停止。

下面是典型输出格式，实际数值取决于你的随机种子和参数设置，但趋势是一致的：

可以看出，当信噪比大于5dB时，MLE已经接近克拉美罗界，继续增加信噪比，误差下降速度变缓。这个现象在波达角估计中非常典型。

5. 工程落地：初值、多起点搜索与参数校验

最后这部分分享几个我从调试这个脚本里总结出来的实用技巧，特别是当你把单信源扩展成多信源或者低信噪比场景时，这些小改动会直接决定算法是否收敛。

5.1 多起点搜索解决局部极值问题

集中似然函数在低信噪比时会出现很多毛刺，单次fminunc容易停在旁瓣位置。我的办法是取网格扫描中最小的三个局部极小值点作为初值，分别运行fminunc，最终取代价最小的结果。实现时不需要同时开线程，顺序跑就行，因为一维优化很快。大致伪代码如下：

[~, idxs] = islocalmin(cost_grid); local_min_idx = find(idxs); [~, sort_idx] = sort(cost_grid(local_min_idx)); best_starts = theta_grid(local_min_idx(sort_idx(1:min(3, length(sort_idx))))); best_theta = NaN; best_cost = inf; for start = best_starts [th, val] = fminunc(@(t) neg_log_likelihood(t, X, M, d_over_lambda), start, options); if val < best_cost best_cost = val; best_theta = th; end end

islocalmin是MATLAB R2017b之后提供的函数，它会返回逻辑索引。取出局部极小值后按代价排序，只取前三个。这个逻辑与全局最优点可能被局部极小值包围时特别有效。注意fminunc返回的val要与网格扫描的代价一致，因为代价函数已经取了负号，所以这里用最小值比较是合理的。

5.2 高精度网格和优化容差的配合

网格扫描步长和优化容差不是无关的。如果你把网格步长设置为0.1度，但StepTolerance仍然是1e-8，后续优化会做很多次微小迭代，其实没有意义。我一般让网格步长在0.5度到1度之间，然后让fminunc在最优网格附近做局部精化。反过来，如果网格步长太粗，比如2度，遇到靠近0度的曲线平坦区，初值误差会让fminunc多跑几十步，甚至漂移到另一个旁瓣。经验值：半波长阵元间距、8阵元、单信源时，网格步长1度足够，阵元数增加到16时，可以放大到1.5度。

5.3 用零度附近的角度测试初始化

一个很实用的调试手段是把theta_true设成0度，查看网格扫描代价曲线。因为0度附近sin(theta)对角度变化最敏感，但代价函数关于0度对称，如果初值落在正负两侧，优化结果可能都是对的。利用这个特性可以验证你的梯度方向和坐标定义是否正确：如果0度时估计值经常跳变到±90度，说明导向矢量的相位因子符号写反了，或者d_over_lambda用了负数。这种错误在MLE脚本里非常隐蔽，因为MUSIC可能对这种符号不敏感，但fminunc对梯度方向是敏感的。

另外，当你要把这段脚本用于实测数据时，别忘了先校准阵列的幅度和相位一致性。仿真里我们假设每个阵元通道响应完全一致，实际系统里会有通道失配，这时最大似然估计会失真，但不会完全失效。一个临时对策是在估计前对X做幅度归一化，让每个阵元的平均功率相同，然后用归一化后的数据走同一套似然函数。这个技巧能救回很多硬件误差，但也只能作为权宜之计，严谨的做法还是用校准源估计出各通道增益相位。

本文还有配套的精品资源，点击获取