简介：面向阵列信号处理研究者与学习者，围绕面阵二维酉矩阵参数估计问题提供一套仿真实现，重点解决多径传播、多个信号源共存时的方位与频率参数估计，适合具备矩阵论和信号处理基础、希望通过代码理解算法原理的读者。压缩包共2个文件，均为m脚本，整体仅1KB，代码精简可直接运行；两个脚本分别承担算法主流程与辅助数据生成演示功能，便于快速复现经典旋转不变技术与二维酉矩阵结合的处理流程并观察输出。已有548人学习下载，广泛用于无线通信、雷达与声学成像等领域的算法研究。学习者可从中获得完整的仿真骨架，包括信号模型设定、面阵配置、二维酉矩阵构建、旋转不变参数估计与结果验证等关键环节，覆盖接收数据生成到参数对比评估的完整链条，有助于理解从一维方法向二维场景扩展的数学原理，为后续改进或工程落地提供简明代码参考。
1. 面阵二维测向为什么绕不开 Unitary-ESPRIT：从复域到实域的降维打击
在均匀矩形面阵上同时估计方位角和俯仰角，是雷达、5G 波束管理、声学阵列测向里最常见的需求。传统 ESPRIT 要在复数域做特征分解，每个矩阵元素都带着实部和虚部，乘法量是实数的四倍，快拍一多整个仿真就慢得让人想摔键盘。Unitary-ESPRIT 利用面阵中心对称的结构，先用酉矩阵把复数接收数据变换到实数域，特征分解从复数域变成实数域，运算量直接省掉约四分之三；同时变换过程自动完成前后向平滑，等效快拍数翻倍，低信噪比下比普通 ESPRIT 更稳。这就是为什么做面阵二维 DOA 的团队，最后几乎都会落到 Unitary-ESPRIT 这个阵列信号处理算法上。适合正在被复数运算量、参数配对和信噪比门限折磨的工程师。
2. 先把面阵的二维数学模型立住：阵列流型、旋转不变性与酉矩阵的作用
2.1 面阵不是线阵的简单堆叠：从 URA 到 L 形阵的流型构造
面阵二维估计的第一步不是写代码，而是把阵列几何和数据排列约定清楚。最常见的布局是均匀矩形面阵（URA），沿 x 方向摆 M 个阵元、沿 y 方向摆 N 个阵元，阵元间距分别记为 dx、dy，总阵元数 MN 个。URA 对 Unitary-ESPRIT 特别友好，因为它的二维导向矢量可以写成 x、y 两个方向导向矢量的 Kronecker 积，流型矩阵结构清晰，酉变换矩阵也能照着同样的 Kronecker 积形式构造。
假设有 K 个远场窄带信号，第 k 个信号的俯仰角为 θk、方位角为 φk，定义方向余弦：
- u_k = sin(θk) cos(φk)
- v_k = sin(θk) sin(φk)
沿 x 方向的导向矢量为：
a_x(u_k) = exp(-j 2π dx x_index^T u_k / λ)
其中 x_index 是阵元的 x 坐标向量。沿 y 方向同理得到 a_y(v_k)。整个面阵对第 k 个信号的导向矢量是：
a_k = kron(a_y(v_k), a_x(u_k))
这里要特别注意 kron 的顺序。我习惯固定用 kron(ay, ax)，即先按 y 方向分块、块内按 x 方向排阵元。这样接收数据矩阵 X 的第 ((n-1)M + m) 行，对应坐标 (x_m, y_n)。这个约定后面所有选择矩阵、酉矩阵都要跟着它走，错一处就全错。
L 形阵和十字阵也能做二维估计，但 Unitary-ESPRIT 的前后向平滑依赖阵列关于中心点对称。L 形阵在转角处天然不对称，强行用酉变换会破坏实域化的前提，所以用 L 形阵的团队通常退化成两个一维估计再配对，或者只在两条臂上分别做酉 ESPRIT，配对成本比 URA 高不少。如果你的项目还没定阵型，选 URA 能省掉后续大量麻烦。
2.2 旋转不变方程怎么在二维上成立：两个方向的子阵选择
ESPRIT 的核心是"两个结构相同的子阵，输出之间差一个固定相位旋转"。一维线阵只有一组旋转不变关系，面阵有两组，分别对应 x 方向和 y 方向，这也是二维估计的信息来源。
沿 x 方向看，把每个 y 行块内的前 M-1 个阵元挑出来组成子阵一，后 M-1 个阵元挑出来组成子阵二。这两个子阵在 x 方向相差一个阵元间距 dx，所以同一信号在子阵二中的导向矢量，是子阵一中的导向矢量乘以 exp(j 2π dx u_k / λ)。写成矩阵形式就是：
A_x2 = A_x1 · Φ_x
其中 Φ_x 是 K×K 的对角矩阵，第 k 个对角元就是 exp(j 2π dx u_k / λ)。沿 y 方向做同样操作，得到：
A_y2 = A_y1 · Φ_y
Φ_y 的第 k 个对角元是 exp(j 2π dy v_k / λ)。只要能从数据里估计出 Φ_x 和 Φ_y，就能解出 u_k 和 v_k，再换算成俯仰角和方位角。
子阵选择用什么实现？在代码里用选择矩阵 J。Jx1 挑每个 y 块内的前 M-1 个阵元，Jx2 挑后 M-1 个：
Jx1 = kron(eye(N), [eye(M-1), zeros(M-1,1)]) Jx2 = kron(eye(N), [zeros(M-1,1), eye(M-1)])
注意这里 eye(N) 在左边，表示对 N 个 y 块每个都做同样的 x 方向选择。y 方向的选择矩阵则把 eye(N-1) 放在左边、eye(M) 放在右边：
Jy1 = kron([eye(N-1), zeros(N-1,1)], eye(M)) Jy2 = kron([zeros(N-1,1), eye(N-1)], eye(M))
这四个选择矩阵的维度分别是 (M-1)N × MN 和 (N-1)M × MN。构造它们时必须和 2.1 节约定的数据排列一致，否则提取到的子阵根本不是平移关系。
2.3 酉变换矩阵：中心对称约束下实域化的原理
酉变换是 Unitary-ESPRIT 区别于标准 ESPRIT 的关键。目标是把复数接收数据矩阵 X（MN×L 维）变换成实数矩阵，之后的特征分解全部在实数域进行。
前提条件：阵列必须是中心对称的，即阵列几何在旋转 180° 后与自身重合。URA 显然满足，这也解释了为什么 2.1 节强调阵元坐标要用中心对称编号。对 M 为偶数的情况，x 方向坐标取：
x_index = -M/2+0.5 : M/2-0.5
例如 M=6 时坐标是 [-2.5, -1.5, -0.5, 0.5, 1.5, 2.5]。这样导向矢量满足 a(u) = Π_M · conj(a(u)) 的共轭对称关系，其中 Π_M 是 M 阶反单位矩阵。
构造酉矩阵 Q 的标准方式是按奇偶分情况。偶数阶：
Q(2n) = (1/√2) · [ I_n, jI_n ; Π_n, -jΠ_n ]
奇数阶：
Q(2n+1) = (1/√2) · [ I_n, 0, jI_n ; 0, √2, 0 ; Π_n, 0, -jΠ_n ]
其中 Π_n 是 n 阶反单位矩阵。这个 Q 满足 Q^H Q = I，且 Π_n · conj(Q) = Q，即列共轭对称。整个面阵的酉矩阵用 Kronecker 积组合：
Q = kron(Q_N, Q_M)
维度是 MN×MN。有了 Q 之后，把原始数据做前后向增强：
Xaug = [X, Π_MN · conj(X) · Π_L]
其中 Π_MN 是 MN 阶反单位矩阵，Π_L 是 L 阶反单位矩阵。这一步把快拍从 L 变成了 2L，相当于把共轭反转后的数据拼接在原数据后面，实现了前后向平滑，噪声得到平均，低信噪比下自相关矩阵更稳。然后做酉变换：
Y = Q^H · Xaug
利用 Q 的共轭对称性质，可以证明 Y 在理论上是实矩阵，数值上用 real() 取实部截断即可。之后对 Y 做特征分解，得到的就是实值信号子空间 Es。相比标准 ESPRIT 在复数域直接对自相关矩阵做特征分解，这里每一步都是实矩阵运算，存储减半、乘法量降到约四分之一。这也是实域化带来的最直接收益。
3. 用 MATLAB 复现面阵 Unitary-ESPRIT：从流型生成到方位/俯仰角输出的完整代码
3.1 生成面阵的数据模型：阵元坐标、快拍与二维角度信号
先写数据生成部分。这一步把接收数据整理成二维数组的过程，和你把 MATLAB 图像转为二维矩阵的预处理很像，关键在于行列索引的顺序不能乱。下面这段代码生成一个 6×8 的 URA 接收数据，包含两个信号源，每个源有独立的俯仰角和方位角。
% 面阵参数 M = 6; % x 方向阵元数 N = 8; % y 方向阵元数 dx = 0.5; % x 方向阵元间距，以波长为单位 dy = 0.5; % y 方向阵元间距 K = 2; % 信号源个数 L = 200; % 快拍数 SNR_dB = 10; % 信噪比，单位 dB % 真实角度：俯仰角 theta、方位角 phi（度） theta = [30, 50]; phi = [20, 60]; % 方向余弦 u = sind(theta) .* cosd(phi); v = sind(theta) .* sind(phi); % 中心对称坐标，保证阵列满足共轭对称前提 x_idx = -M/2+0.5 : M/2-0.5; % M=6 -> [-2.5 -1.5 -0.5 0.5 1.5 2.5] y_idx = -N/2+0.5 : N/2-0.5; % 两个方向的导向矢量，并组合成二维流型矩阵 ax = exp(-1j*2*pi*dx * (x_idx.' * u)); % M x K ay = exp(-1j*2*pi*dy * (y_idx.' * v)); % N x K A = kron(ay, ax); % (M*N) x K，列对应各信号源 % 信号与高斯白噪声 S = (randn(K, L) + 1j*randn(K, L)) / sqrt(2); sigma = 10^(-SNR_dB/20); Nn = (randn(M*N, L) + 1j*randn(M*N, L)) / sqrt(2); X = A*S + sigma * Nn;这段代码里最容易被忽略的是 x_idx 和 y_idx 的取值。如果用 0:M-1 这种从零开始的编号，阵列的对称中心就不在原点，后面构造 Q 矩阵时共轭对称性质不成立，Y = real(Y) 这一步会丢掉有效信息。另一个要注意的是 A 用 kron(ay, ax) 而不是 kron(ax, ay)，这个顺序和后续选择矩阵、酉矩阵的构造顺序必须完全统一。信号 S 用复高斯分布生成，功率归一化到单位功率，噪声功率由 SNR_dB 计算得到。
3.2 核心实现：酉变换、实域特征分解与二维参数估计
核心估计函数放在独立文件里，输入接收数据 X 和面阵参数，输出估计的俯仰角与方位角。完整实现如下：
function [theta_est, phi_est] = unitary_esprit_2d(X, M, N, K, dx, dy) [MN, L] = size(X); % 构造酉矩阵 Q，kron 顺序与流型 kron(ay, ax) 保持一致 Q = kron(get_Q(N), get_Q(M)); % 前后向增强：把数据扩展成 2L 列 Pi_MN = fliplr(eye(MN)); Pi_L = fliplr(eye(L)); Xaug = [X, Pi_MN * conj(X) * Pi_L]; % 酉变换并取实部，理论上变换后为实矩阵 Y = Q' * Xaug; Y = real(Y); % 实域特征分解，取最大的 K 个特征向量作为信号子空间 R = Y * Y.'; [Evec, Eval] = eig(R); [~, idx] = sort(diag(Eval), 'descend'); Es = Evec(:, idx(1:K)); % 选择矩阵：x 方向取前/后 M-1 个阵元，y 方向取前/后 N-1 个阵元 Jx1 = kron(eye(N), [eye(M-1), zeros(M-1, 1)]); Jx2 = kron(eye(N), [zeros(M-1, 1), eye(M-1)]); Jy1 = kron([eye(N-1), zeros(N-1, 1)], eye(M)); Jy2 = kron([zeros(N-1, 1), eye(N-1)], eye(M)); % 把选择矩阵变换到酉域，得到实值的广义选择矩阵 Q_Mm = get_Q(M-1); Q_Nn = get_Q(N-1); Kx1 = 2 * real(kron(Q_N, Q_Mm)' * Jx1 * Q); Kx2 = 2 * real(kron(Q_N, Q_Mm)' * Jx2 * Q); Ky1 = 2 * real(kron(Q_Nn, Q_M)' * Jy1 * Q); Ky2 = 2 * real(kron(Q_Nn, Q_M)' * Jy2 * Q); % 两个方向的旋转矩阵，用最小二乘解 Psi_x = (Kx1*Es) \ (Kx2*Es); Psi_y = (Ky1*Es) \ (Ky2*Es); % 配对关键：用 Psi_x 的特征向量去对角化 Psi_y [T, Dx] = eig(Psi_x); Dy = T \ Psi_y * T; % 提取旋转相位 gamma_x = angle(diag(Dx)); gamma_y = angle(diag(Dy)); % 换算方向余弦与角度 u_est = gamma_x / (2*pi*dx); v_est = gamma_y / (2*pi*dy); phi_est = atan2d(v_est, u_est); temp = u_est.^2 + v_est.^2; temp(temp > 1) = NaN; % 超出单位圆视为无效 theta_est = asind(sqrt(temp)); end配套的酉矩阵生成函数：
function Q = get_Q(n) if mod(n, 2) == 0 n2 = n / 2; Pi = fliplr(eye(n2)); Q = 1/sqrt(2) * [eye(n2), 1j*eye(n2); Pi, -1j*Pi]; else n2 = (n-1) / 2; Pi = fliplr(eye(n2)); Q = 1/sqrt(2) * [eye(n2), zeros(n2,1), 1j*eye(n2); ... zeros(1,n2), sqrt(2), zeros(1,n2); ... Pi, zeros(n2,1), -1j*Pi]; end end逻辑上分为五段。第一段构造酉矩阵 Q，Xaug 是前后向增强后的数据，这一段直接决定实数化是否成立。第二段对实矩阵 R 做特征分解，信号子空间 Es 的列数必须等于信号源数 K，K 给多了会把噪声子空间带进来，给少了则旋转矩阵不满秩。第三段把四个选择矩阵变换到酉域，Kx1 的维度是 (M-1)N × MN，左乘 Es 后得到 (M-1)N × K 的矩阵，这一步不必显式计算阵列流型，而是直接在信号子空间上提取旋转不变关系。第四段的 Psi_x、Psi_y 是 K×K 矩阵，它们的特征值就携带了方向余弦信息。第五段的配对技巧会在 3.4 节单独解释。
3.3 参数设置对照表：阵元数、快拍数、信噪比对结果的影响
参数不是越大越好，不同参数对结果的约束方向完全不同。下面这张表是我在调试时总结的对应关系。
| 参数 | 常见取值 | 对结果的影响 |
|---|---|---|
| M、N 阵元数 | 4～16 | 阵元越多分辨率越高、抗噪越强，但互耦越明显，特征分解复杂度按 MN 的平方增长 |
| dx、dy 间距 | 0.5λ | 超过 0.5λ 会出现栅瓣，小于 0.5λ 会压缩方向余弦的可估计范围 |
| L 快拍数 | 100～1000 | L 越大自相关矩阵越稳，低信噪比下至少要 200 以上 |
| SNR | 0～20 dB | 低于 0 dB 时建议增大快拍或增加阵元，否则阈值效应会突然拉高 RMSE |
| K 信号源数 | 事先已知或 MDL 估计 | K 偏大会导致旋转矩阵出现伪特征值，K 偏小则漏掉真实角度 |
一个常见的误用是把阵元数加得很大而快拍数不变。阵元数增大后，同一信噪比下自相关矩阵的特征值散布更宽，需要的快拍数反而要更多。我一般会按"每增加一倍阵元，快拍数至少增加 50%"来粗调。
3.4 输出配对：为什么特征值配对错了角度就乱了
面阵二维估计的最后一个关键步骤是配对，即保证解出的第 k 个 u_k 和 v_k 属于同一个信号源。很多初版实现是分别对 Psi_x 和 Psi_y 做 eig，再按特征值大小排序对应起来。这个做法在低信噪比下几乎必错，因为两个矩阵的特征值排序互不相关，一旦噪声让某个特征值的模发生波动，序号就错位，最后画出来的散点图会出现多个错误的 (θ, φ) 组合。
正确做法利用了 Psi_x 和 Psi_y 共享同一组特征向量这个性质。在无噪声条件下，Psi_x 和 Psi_y 可以被同一个非奇异矩阵 T 对角化，所以代码里先对 Psi_x 做一次 eig 得到 T，再用 T 去相似变换 Psi_y：
[T, Dx] = eig(Psi_x); Dy = T \ Psi_y * T;Dx 和 Dy 的第 k 个对角元天然对应同一个信号源，配对自动完成，不需要排序。只有一行代码的差别，却是整个算法里最容易被忽略的翻车点。如果实测中仍然发现配对偶发错乱，可以对 Psi_x + 1j*Psi_y 做联合对角化近似，但大多数场景下上面的相似变换已经足够。
4. 面阵 Unitary-ESPRIT 的 5 个避坑点：配对、栅瓣、维数与边界
4.1 方位角和俯仰角配对错乱：现象、原因与一次相似变换修复
现象：估计出的角度个数是对的，但散点图上出现若干个位置完全离谱的点，真实源在 (30°, 20°)，结果却跑出 (30°, 60°) 这种组合，且每次蒙特卡洛跑出来的错配位置都不一样。
原因：对 Psi_x 和 Psi_y 分别做特征分解后按特征值排序强行对应。噪声会让两个矩阵特征值的模和相位发生不同程度扰动，排序会错位，导致 x 方向第 k 个分量和 y 方向第 j 个分量被拼在一起。
解决：只对 Psi_x 做一次特征分解得到特征向量矩阵 T，然后用 Dy = T \ Psi_y * T 去对角化 Psi_y。取 diag(Dy) 时天然和 diag(Dx) 按同一顺序排列。这一步是配对问题的后悔药，改一行代码就能稳定输出。
4.2 阵元间距超过半波长：相位卷绕让二维估计直接翻车
现象：把 dx、dy 改成 0.8λ 后，中高信噪比下原本正确的角度估计突然出现大偏差，有些源的角度直接跑到视场外，而且改变随机种子后偏差位置剧烈变化。
原因：旋转相位 gamma = 2π d u / λ 的真实值超过 π 后，MATLAB 的 angle() 只能返回 [-π, π] 区间内的主值，相位发生卷绕，解出来的方向余弦丢失了 2π 整数倍的信息，这就是栅瓣效应。
解决：阵元间距严格取 0.5λ 以内。0.5λ 是临界值，实际工程中我会取 0.45λ，留一点余量避免频率波动导致有效间距越界。如果因为孔径约束必须用大间距，就需要额外的解模糊步骤，比如用波束扫描先粗估角度范围，再对卷绕数做搜索，复杂度会上升一个量级。
4.3 kron 顺序与数据排列不一致：流型对的但角度全错
现象：代码逻辑检查了几遍都找不出问题，酉矩阵、选择矩阵、特征分解每一步都很规范，但估计出的角度里 x 方向的信息和 y 方向的信息互换了，或者干脆完全错乱。
原因：A = kron(ay, ax)，但构造 Q 时写成了 kron(Q_M, Q_N)，或者 Jx1 写成了 kron([eye(M-1), zeros(M-1,1)], eye(N))。Kronecker 积的顺序一旦和流型约定不一致，选择矩阵选出来的就不是物理上平移的子阵，旋转矩阵自然毫无意义。
解决：固定一套约定并贯彻到底。我习惯在代码文件头部写注释：数据排列 kron(ay, ax)，阵元序号从 1 到 MN，先变 x 索引。然后所有矩阵构造都对照这个注释检查。这个坑我踩过一次之后就再也不敢凭感觉写 kron 参数了。
4.4 快拍太少时特征值散裂：信号子空间维数怎么判断
现象：L 只有 20，SNR 也不高，特征分解后特征值谱没有明显的台阶，取最大的 K 个特征向量后，Psi_x 的奇异值发散，估计角度波动非常大。
原因：自相关矩阵是用样本估计的，快拍太少时小特征值被噪声抬高、大特征值被压低，信号子空间和噪声子空间的边界变模糊。虽然 Unitary-ESPRIT 的前后向增强把等效快拍变成了 2L，但 L 本身太小时增强也救不回来。
解决：先做信源数估计，用 MDL 或 AIC 准则判断 K，而不是拍脑袋填。当 L < 2K 时，直接放弃特征分解类方法，改用稀疏恢复或波束扫描。另外可以在 Y 变换前对数据做一次空间平滑，代价是有效孔径变小，但低信噪比下更耐用。
4.5 角度接近 0° 或 90°：asin 越界与误差放大
现象：真实俯仰角是 85° 时，估计结果经常落在 80°～90° 之间，多次蒙特卡洛的方差比其他角度大好几倍；有时 u_est^2 + v_est^2 超过 1，程序输出 NaN。
原因：θ 接近 90° 时，sinθ 对角度变化的导数趋于零，同样的方向余弦误差换算成角度误差会被放大。同时方向余弦接近 1，噪声稍大平方和就超界，asin 定义域失效。
解决：如果目标场景有 60° 以上的大俯仰角，改用天顶角定义或者直接在 u-v 域输出结果，不要强行转换。工程上我会在代码里保留 temp > 1 的判据，把所有越界点标记为异常而不是静默丢弃，这样后续数据分析时能区分"算法失效"和"角度本身靠近边界"。
5. 从仿真到实测：验证二维角度估计性能的三个习惯
5.1 用 RMSE 对比曲线确认算法真的收敛到 CRB
跑通一次仿真只说明代码没报错，不说明算法是对的。我的第一个验证动作是固定阵型和信噪比，跑 100 次蒙特卡洛，统计每个源的 RMSE，画成随 SNR 变化的曲线，再叠加上 CRB（克拉美-罗界）曲线。中高信噪比下 RMSE 应该贴着 CRB 走，低信噪比下陡然抬高属于正常的阈值效应。如果曲线整体偏离 CRB 好几倍，先回头检查数据排列和配对逻辑。注意计算角度 RMSE 时要做差值卷绕处理，否则 179° 和 -179° 的平均差会算成 358°。
5.2 多信源场景先看散点图，再谈统计指标
K 大于等于 3 之后，单一 RMSE 指标会掩盖配对错误，因为配错的点既可能出现在某一个源的统计里，也可能被平均掉。我现在的习惯是每次蒙特卡洛都画一幅 u-v 域散点图，三个真实源的位置标成叉号，估计点标成圆点。如果每个真实源周围都聚着一簇点，配对基本可靠；只要看到有圆点跑到两个簇中间，直接判定这次实现有配对问题。散点图也方便发现系统偏差，比如所有点整体向某个方向偏移 2°，这种偏差 RMSE 曲线未必看得清楚。
5.3 与波束扫描结果交叉验证，排除系统偏差
仿真环境里流型矩阵是精确已知的，但实测定标误差、阵元位置误差、互耦都会造成系统偏差。我的做法是用同一批数据跑一遍常规复数域 ESPRIT，再做一次二维波束扫描（或 MUSIC 谱峰搜索），对比三种方法给出的峰值位置。三者偏差在 1° 以内，基本可以认定实现没有大问题；如果 Unitary-ESPRIT 和波束扫描一致、但和复数域 ESPRIT 不一致，多半是酉变换里的共轭对称处理有误。如果三者互不一致，先怀疑数据采集时的阵元编号映射，而不是算法本身。
这套验证流程的最后一个好处是让实测阶段有据可依。我现在每接手一个新阵列配置，都会先用 3.1 节的生成器造一组已知角度的数据，跑一遍完整流程确认 RMSE 贴 CRB，再上实测数据。没有这一步兜底，实测数据里的任何异常都分不清是算法问题、通道问题还是标定问题，那就只能靠玄学调参了。希望这个从仿真到验证的最小闭环能帮到你，至少能让你少走几趟弯路。
本文还有配套的精品资源，点击获取
