简介面向无线通信、雷达探测、声学成像等领域中的二维到达角DOA估计学习者这份MATLAB实现包聚焦基于增广矩阵束的L型阵列DOA估计方法。代码将相互垂直的两个子阵列接收数据融合为增广矩阵通过矩阵束处理联合求解信号的方位角与仰角适合需要理解2D DOA算法本质、并希望借助可运行程序验证理论的研究人员。压缩包内共3个文件包含2个MATLAB脚本.m和1个自动保存文件.asv总大小仅2KB体积轻量、结构清晰便于直接阅读与修改。已有262人浏览/学习。资源覆盖L型阵列建模、增广矩阵构造、矩阵束形成、DOA解算与结果可视化等关键步骤通过修改阵列间距、噪声水平、信号频率等参数即可观察估计性能变化也可作为进一步学习ESPRIT、MUSIC等DOA估计算法的实验起点具有较好的参考与扩展价值。1. 二维DOA估计里的L阵列与Matrix Pencil为什么这组合值得上手做阵列信号处理的人,碰到2D DOA估计的第一反应往往是上均匀矩形阵配MUSIC或者ESPRIT。但真到了项目里,你会发现阵元数量、安装空间、实时性三项一卡,平面阵方案经常是第一个被砍掉的。反而是L阵列加上Matrix Pencil(矩阵束)这条路线,用两条正交的均匀线阵,就能同时估出方位角和俯仰角,阵元数少、计算量小、单快拍也能跑,只是代价藏在配对这一步。这篇笔记把L阵列的信号建模、Hankel矩阵构造、SVD截断、极点提取完整拆开,给出可直接复现的MATLAB框架,再把镜像模糊、pencil参数选择、相干源这些容易翻车的点逐一说明。适合正在搭DOA仿真验证、或想用最小硬件成本验证二维测向算法的工程师。2. L型阵列的信号建模两条线阵如何拼出二维角度2.1 为什么选L阵列而不是平面阵L阵列的结构很简单:沿x轴放M个阵元,沿z轴再放M个阵元,原点共用,实际物理阵元数是2M-1。它和均匀矩形阵的本质区别是:矩形阵直接采样整个二维波前,L阵列只采样波前在两个正交方向上的投影。投影信息量少了,但换来两个实用优势——阵元少一半以上、边缘安装友好(比如机身蒙皮、车厢边缘),代价是角度解算多了配对这步。从自由度看,每个子阵的均匀线阵最多分辨M-1个信源,两条子阵合计也是这个量级,比同阵元数的平面阵少了约一个维度。所以L阵列适合信源数不大、安装受限、对实时性敏感的场合。另外,L阵列天然避开了均匀矩形阵在仰角接近天顶时的模糊问题——z轴子阵测的是sin(俯仰),俯仰在正负90度内是单调的。2.2 方向向量约定与信号模型约定信号来波方向的单位向量为u [cos(el)·cos(az), cos(el)·sin(az), sin(el)],其中az是相对x轴的方位角,el是相对水平面的俯仰角。x轴子阵第i个阵元相对原点的相位差是2π(d/λ)·i·cos(el)·cos(az),z轴子阵是2π(d/λ)·i·sin(el)。% 方向向量构造函数, d_lambda 阵元间距/波长 % az, el 输入弧度, 返回 M x 1 复数方向向量 function a steer_vec_x(az, el, M, d_lambda) a exp(1j * 2 * pi * d_lambda * (0:M-1). * cos(el) * cos(az)); end function a steer_vec_z(el, M, d_lambda) a exp(1j * 2 * pi * d_lambda * (0:M-1). * sin(el)); end这里用exp(j·k·r)的约定,后面所有相位提取都按这个符号走,别中途换约定。d_lambda取0.5是经典配置,一是满足空间采样定理,二是相位范围刚好落在[-π,π]内,angle()取相位时不会卷绕。2.3 两个子阵各自估计,配对才是主角x轴子阵的极点和z轴子阵的极点分别对应两个投影量:x轴极点相位 ∝ cos(el)·cos(az),方位和俯仰耦合在一起;z轴极点相位 ∝ sin(el),只含俯仰。问题是两组极点各自独立解出来,谁和谁属于同一个信源并不知道。这就是L阵列做2D DOA的核心难点,叫做配对。常见配对思路有三类:用每个分量的复幅度做最近邻、用特征向量联合对角化、直接枚举配对后做数据重构残差检验。第三种最稳,后面第5章会给完整实现。3. 用Matrix Pencil做2D DOA估计可复现的MATLAB实现3.1 矩阵束的一句话直觉Matrix Pencil的核心思想是:把一段观测序列看成若干个复指数分量的叠加,每个分量对应一个极点z_k。如果能从数据里分离出这些极点,角度就藏在极点相位里。做法是把观测序列重排成一个Hankel矩阵,利用Hankel矩阵相邻两行的平移关系,构造一个矩阵束,它的广义特征值就是极点。对阵列数据来说,每个阵元的观测就相当于一个采样点,阵元间距就是采样间隔。所以同一个pencil框架可以直接搬到阵列域:一个子阵的快拍数据,构造Hankel矩阵,求极点,极点的相位就是相邻阵元间的波程相位增量。这也是为什么pencil在单快拍下依然能工作——它靠的是Hankel矩阵的秩结构,而不是多快拍的统计平均。3.2 完整可跑的最小示例:仿真数据生成先定场景:8元子阵、2个不相关窄带信号、100快拍、信噪比15dB。clear; clc; rng(2024); M 8; % 每个子阵的阵元数(含原点,原点共用) d_lambda 0.5; % 阵元间距/波长 N 100; % 快拍数 K 2; % 信源数 snr_db 15; az_true [30, -20] * pi/180; % 方位角, 注意线阵对称模糊,后面单独讲 el_true [25, 40] * pi/180; % 俯仰角 S (randn(K, N) 1j*randn(K, N)) / sqrt(2); % 复高斯信源, 不相关 Xx steer_vec_x(az_true(1), el_true(1), M, d_lambda) * S(1,:) ... steer_vec_x(az_true(2), el_true(2), M, d_lambda) * S(2,:); Xz steer_vec_z(el_true(1), M, d_lambda) * S(1,:) ... steer_vec_z(el_true(2), M, d_lambda) * S(2,:); % 加复高斯白噪声 noise_power 10^(-snr_db/10); Xx Xx sqrt(noise_power/2)*(randn(M,N) 1j*randn(M,N)); Xz Xz sqrt(noise_power/2)*(randn(M,N) 1j*randn(M,N));逻辑说明:信源矩阵S是K×N,每一行是一个信源的复包络。两个方向向量分别作用在S的两行上,叠加后得到两个子阵的观测。注意这里方向向量用的是第一节的约定,信源之间不相关,pencil不需要协方差矩阵,直接吃原始快拍数据。参数说明:rng(2024)固定随机种子,保证每次跑出来的图一致;如果换信噪比,只需要改snr_db;信源数K要事先给准,这是所有子空间类方法的前提。3.3 核心函数:matrix_pencil_poles这个函数输入某个子阵的M×N快拍数据,输出K个极点。关键点在于多快拍的处理方式——对每个快拍构造Hankel块再水平拼接,拼接后取左奇异向量U而不是右奇异向量V。function poles matrix_pencil_poles(X, K, L) % X: M x N 子阵接收数据 % K: 信源个数(用于SVD截断) % L: pencil参数, 建议取 M/3 ~ M/2 之间 % 输出: 1 x K 极点, 模接近1, 相位含角度信息 [M, N] size(X); m M - L 1; % Hankel矩阵行数 Xe zeros(m, L * N); % 增强矩阵: 每个快拍的Hankel块水平拼接 for n 1:N x X(:, n); Hn zeros(m, L); for jj 1:L % 显式填充Hankel矩阵 Hn(:, jj) x(jj : jj m - 1); end Xe(:, (n-1)*L 1 : n*L) Hn; end [U, ~, ~] svd(Xe, econ); Us U(:, 1:K); % 信号子空间: 取左奇异向量 U1 Us(1:end-1, :); % 行方向平移不变性 U2 Us(2:end, :); Psi pinv(U1) * U2; % 旋转算子 poles eig(Psi).; end逻辑说明:Hankel矩阵的每一行是观测序列的一个窗口滑动,行索引对应阵元位移方向。拼接了N个快拍后,增强矩阵是m×(L·N),SVD得到的左奇异向量U的行数仍是m,行之间保留了平移结构,所以U1、U2的旋转关系成立。如果这时候去取V做平移,会得到L·N行,平移结构被多快拍打乱,解出来全是错的——这是实际实现里最常见的翻车点之一。参数说明:L就是pencil参数,也叫束参数。L太小时Hankel矩阵的列数不够,奇异值截断后子空间不稳;L太大时m变小,信号子空间的行数不足。经验区间是M/3到M/2,同时要满足L-1≥K且m-1≥K,否则旋转算子没有足够的行做最小二乘。我一般先用Lfloor(M/2)跑通,再扫一遍L验证稳定性。3.4 配对与角度解算对两个子阵分别调用matrix_pencil_poles,得到两组极点。配对用的是幅度一致性:同一个信源在x轴和z轴上的复包络幅度应该一致,先估计每个极点的幅度,枚举全排列选幅度差最小的组合。function [az_est, el_est] pair_by_amplitude(px, pz, Xx, Xz, d_lambda) % px, pz: 两个子阵估计出的极点(1 x K) % Xx, Xz: 两个子阵的观测数据, 用于幅度拟合 % 返回配对后的方位角和俯仰角(弧度) K numel(px); M size(Xx, 1); % 构造范德蒙德矩阵, 用最小二乘估计每个极点的复幅度 Zx zeros(M, K); Zz zeros(M, K); for k 1:K Zx(:, k) px(k).^(0:M-1).; Zz(:, k) pz(k).^(0:M-1).; end Rx (Zx * Zx) \ (Zx * Xx); % K x N, 每个快拍每个极点的复幅度 Rz (Zz * Zz) \ (Zz * Xz); ampx mean(abs(Rx), 2); % K x 1 ampz mean(abs(Rz), 2); % 枚举全排列, 让幅度差平方和最小 Perms perms(1:K); best_cost inf; best_perm []; for ii 1:size(Perms, 1) p Perms(ii, :); cost sum((ampx - ampz(p)).^2); if cost best_cost best_cost cost; best_perm p; end end % 按配对结果解算角度 el_est zeros(1, K); az_est zeros(1, K); for k 1:K el_est(k) asin(angle(pz(best_perm(k))) / (2*pi*d_lambda)); cos_el cos(el_est(k)); cos_az angle(px(k)) / (2*pi*d_lambda) / cos_el; az_est(k) acos(cos_az); end end逻辑说明:先从极点和数据里把幅度估计出来,幅度相近的极点就是同一个信源。这一步要求两个子阵在原点共用且幅相响应一致,实际工程中如果两个子阵的通道增益有差异,要先校准,否则幅度配对会系统性错配。角度解算时,俯仰角直接从z轴极点反解asin,方位角再用x轴极点的相位除以cos(俯仰)后反解acos。参数说明:K小时全排列枚举没问题,K超过6就不现实了,那时改用贪心或匈牙利算法做二分匹配。另外,幅度配对的精度受极点的条件数影响,如果两个信源的俯仰角很接近,z轴极点几乎重合,Zz矩阵接近奇异,幅度估计会炸,这种场景下配对应改用第5章的重构残差法。3.5 跑一遍并验证结果把上面两节拼起来,调用主流程:% 主流程 L floor(M/2); % pencil参数取4 px matrix_pencil_poles(Xx, K, L); pz matrix_pencil_poles(Xz, K, L); [az_est, el_est] pair_by_amplitude(px, pz, Xx, Xz, d_lambda); disp(真实值(角度):); disp([az_true; el_true] * 180/pi); disp(估计值(角度):); disp([az_est; el_est] * 180/pi);正常跑通时,打印出来的估计值应该和真实值接近,误差在0.5度以内(15dB、100快拍下)。如果对不上,先检查三件事:一是px和pz的极点模值是否都接近1,偏离太多说明L参数或信源数K不对;二是配对函数返回的配对结果,幅度差是否明显小于错配情况;三是snr_db是否太低,10dB以下幅度配对开始出现概率性错配。4. 常见问题排查与避坑:为什么你的结果总在抖动4.1 镜像模糊:方位角正负分不清现象:真实方位角是负30度,估计出来却是正30度,而且两个信源的角度估计相互交叉,时好时坏。原因:x轴均匀线阵对±az方向完全对称,cos(az)相同,投影信息丢失了符号。这跟算法无关,是阵列几何的固有问题,不是pencil能解决的。解决:三个办法,按成本排序。一,利用场景先验,比如雷达只扫前向,直接把acos主值映射到需要的象限;二,把x轴换成非均匀阵,打破对称性;三,在原点附近加一个错位的辅助阵元,用辅助阵元的相位差定符号。仿真阶段最省事的是先设az在0到90度之间,验证算法链路后再处理实际象限。4.2 pencil参数L乱选导致极点飞了现象:L取M-1时,估计角度剧烈抖动;L取2时,极点倒是稳定,但角度分辨率明显变差,两个靠近的信源分不开。原因:L决定Hankel矩阵的行数mM-L1和列数L。L太大,m太小,信号子空间行数不足,旋转矩阵的最小二乘解病态;L太小,Hankel矩阵的列数不够,奇异值截断后的子空间无法容纳足够的信息。解决:在L∈[3, M-2]范围内扫描,画估计角度的稳定性曲线,取一个平坦区间内的值。我常用的做法是Lfloor(M/2)起步,然后分别试Lfloor(M/3)和Lfloor(2M/3),三个结果一致就放心了。扫描时可以顺便监控eig(Psi)的模值,所有极点模值落在0.95到1.05之间才算正常。4.3 相干信源导致秩亏,极点少一个现象:两个信源是同一信号的多径,幅度相关,结果只估出一个角度,或者估计的角度是两个真实角的中间值。原因:pencil方法同样基于子空间,相干信源让Hankel增强矩阵的有效秩低于K,SVD截断后信号子空间塌缩,丢失信息。解决:做前后向平滑,把数据翻折拼接后再进pencil。在matrix_pencil_poles外面包一层:% 前后向平滑预处理: 对M x N数据 X_fb [X, conj(flipud(X))]; % M x 2N % 然后把X_fb传给matrix_pencil_poles逻辑说明:flipud把阵元顺序反转并取共轭,等同于把阵列的镜像观测补进来。相干信源在正反两个方向上的合成响应不同,拼接后协方差恢复满秩。代价是黑匣子式的自由度减半——等效快拍数翻倍,但对相干源的分辨能力有时反而不如加空间平滑彻底。工程上如果做的是多径场景,建议直接上空间平滑,每3到4个阵元一组做重叠子阵。4.4 多快拍拼接后取错奇异向量现象:单快拍数据跑MP结果挺正常,改成多快拍后输出全乱。很多人翻论文里单快拍的公式,把取V那一步原样搬到多快拍代码里,结果彻底翻车。原因:单快拍时,Hankel矩阵的右奇异向量V的行数等于L,行方向对应序列的滑动窗口,平移结构在V上成立。但多快拍水平拼接后,V的行数变成L×N,不同快拍的列混在一起,平移结构不复存在。此时只有左奇异向量U的行数仍是m,保住了阵元位移结构,取U才正确。解决:严格按第3章的写法——先对每个快拍构造Hankel块,水平拼接成增强矩阵,然后取左奇异向量U做U1/U2旋转。排查时打印size(Us),确认行数等于M-L1;如果是L或L×N,说明取错了。4.5 俯仰角接近90度或d/lambda太大,角度爆炸现象:其他角度的结果都正常,唯独某一组俯仰角接近90度时,方位角误差突然变成十几度,甚至出现相位跳变。原因:俯仰角接近90度时,来波方向接近天顶,x轴子阵的投影cos(el)·cos(az)趋近于0,几乎不含方位信息,噪声被cos(el)除法放大。另外,如果d/lambda超过0.5,相邻阵元相位差可能超过±π,angle()取相位时发生卷绕,角度直接跳变。解决:一是把d/lambda控制在0.4到0.5之间,别为了增大孔径硬上0.6;二是对解算结果做合理性检查,当cos(el)小于0.1时,直接标记该信源方位不可估,不输出数值;三是如果应用场景确实有高俯仰角目标,考虑改用双L阵列或加一条y轴,用y轴子阵补方位信息。仿真时想绕开这个坑,把el_true设在10到70度之间。5. 让配对更稳的重构残差法:一个可以直接替换的进阶技巧幅度配对在15dB以上、极点分离度好的时候够用,但信噪比掉到5dB,或者两个信源的俯仰角只差5度时,幅度估计本身就被噪声污染,幅度差最小准则开始出错。重构残差法的思路完全不同:对每一种候选配对,直接用估计出的角度构造联合方向矩阵,用最小二乘反演出信源波形,再用这个波形重构整个双轴观测数据,计算重构误差。正确配对必须能同时解释两条子阵的数据,残差一定最小。这个准则不依赖幅度估计,天然免疫幅度噪声。function [az_best, el_best] pair_by_reconstruction(px, pz, Xx, Xz, d_lambda) % 枚举所有配对, 用联合重构残差选最优 % 适合低信噪比或俯仰角接近的场景 K numel(px); [M, ~] size(Xx); D [Xx; Xz]; % 2M x N 联合观测 Perms perms(1:K); best_res inf; az_best zeros(1, K); el_best zeros(1, K); for ii 1:size(Perms, 1) p Perms(ii, :); az_tmp zeros(1, K); el_tmp zeros(1, K); for k 1:K el_tmp(k) asin(angle(pz(p(k))) / (2*pi*d_lambda)); cos_az angle(px(k)) / (2*pi*d_lambda) / cos(el_tmp(k)); az_tmp(k) acos(cos_az); end % 构造联合方向矩阵: 上半块x轴, 下半块z轴 A zeros(2*M, K); for k 1:K A(1:M, k) exp(1j*2*pi*d_lambda*(0:M-1).*cos(el_tmp(k))*cos(az_tmp(k))); A(M1:2*M, k) exp(1j*2*pi*d_lambda*(0:M-1).*sin(el_tmp(k))); end S_hat A \ D; % K x N 最小二乘信源波形 res norm(D - A * S_hat, fro); if res best_res best_res res; az_best az_tmp; el_best el_tmp; end end end逻辑说明:关键在于A和D的拼接方式,同一个信源在x轴和z轴的复振幅必须一致,这是联合重构能够成立的物理前提。S_hat是一次性从联合观测里反演出来的,它同时解释了x轴和z轴的数据,残差res越小说明这个配对与两个子阵的观测越自洽。枚举全排列时,K4以内运算量都可以接受;K更大时改成随机采样候选配对,或者先用幅度配对缩小候选集,再用重构残差决赛。一个小技巧:拿到重构残差最小的配对后,可以在az和el附近各做±2度的小网格搜索,把每个角度组合再代进重构残差,取最小值对应的角度作为最终输出。这就相当于做了一次简化版的最大似然精估计,通常能把误差再压下去0.2到0.3度,代价只是多跑几十次残差计算。我早期做L阵列的时候,一直用幅度配对,觉得简单直接,结果项目里遇到两个目标俯仰角只差4度的场景,输出错配率达到三成,当时以为是阵列标定问题,折腾了快两周。后来换成重构残差法,同一批数据错配率直接归零。现在我的习惯是仿真阶段两个方法都跑,幅度配对出结果,重构残差做仲裁,两个一致才认为这组角度可信。这个习惯帮我挡掉了不少低信噪比下的隐性错配,希望也能帮到你。本文还有配套的精品资源点击获取 L阵列结合Matrix Pencil的二维DOA估计：原理与MATLAB实现
简介面向无线通信、雷达探测、声学成像等领域中的二维到达角DOA估计学习者这份MATLAB实现包聚焦基于增广矩阵束的L型阵列DOA估计方法。代码将相互垂直的两个子阵列接收数据融合为增广矩阵通过矩阵束处理联合求解信号的方位角与仰角适合需要理解2D DOA算法本质、并希望借助可运行程序验证理论的研究人员。压缩包内共3个文件包含2个MATLAB脚本.m和1个自动保存文件.asv总大小仅2KB体积轻量、结构清晰便于直接阅读与修改。已有262人浏览/学习。资源覆盖L型阵列建模、增广矩阵构造、矩阵束形成、DOA解算与结果可视化等关键步骤通过修改阵列间距、噪声水平、信号频率等参数即可观察估计性能变化也可作为进一步学习ESPRIT、MUSIC等DOA估计算法的实验起点具有较好的参考与扩展价值。1. 二维DOA估计里的L阵列与Matrix Pencil为什么这组合值得上手做阵列信号处理的人,碰到2D DOA估计的第一反应往往是上均匀矩形阵配MUSIC或者ESPRIT。但真到了项目里,你会发现阵元数量、安装空间、实时性三项一卡,平面阵方案经常是第一个被砍掉的。反而是L阵列加上Matrix Pencil(矩阵束)这条路线,用两条正交的均匀线阵,就能同时估出方位角和俯仰角,阵元数少、计算量小、单快拍也能跑,只是代价藏在配对这一步。这篇笔记把L阵列的信号建模、Hankel矩阵构造、SVD截断、极点提取完整拆开,给出可直接复现的MATLAB框架,再把镜像模糊、pencil参数选择、相干源这些容易翻车的点逐一说明。适合正在搭DOA仿真验证、或想用最小硬件成本验证二维测向算法的工程师。2. L型阵列的信号建模两条线阵如何拼出二维角度2.1 为什么选L阵列而不是平面阵L阵列的结构很简单:沿x轴放M个阵元,沿z轴再放M个阵元,原点共用,实际物理阵元数是2M-1。它和均匀矩形阵的本质区别是:矩形阵直接采样整个二维波前,L阵列只采样波前在两个正交方向上的投影。投影信息量少了,但换来两个实用优势——阵元少一半以上、边缘安装友好(比如机身蒙皮、车厢边缘),代价是角度解算多了配对这步。从自由度看,每个子阵的均匀线阵最多分辨M-1个信源,两条子阵合计也是这个量级,比同阵元数的平面阵少了约一个维度。所以L阵列适合信源数不大、安装受限、对实时性敏感的场合。另外,L阵列天然避开了均匀矩形阵在仰角接近天顶时的模糊问题——z轴子阵测的是sin(俯仰),俯仰在正负90度内是单调的。2.2 方向向量约定与信号模型约定信号来波方向的单位向量为u [cos(el)·cos(az), cos(el)·sin(az), sin(el)],其中az是相对x轴的方位角,el是相对水平面的俯仰角。x轴子阵第i个阵元相对原点的相位差是2π(d/λ)·i·cos(el)·cos(az),z轴子阵是2π(d/λ)·i·sin(el)。% 方向向量构造函数, d_lambda 阵元间距/波长 % az, el 输入弧度, 返回 M x 1 复数方向向量 function a steer_vec_x(az, el, M, d_lambda) a exp(1j * 2 * pi * d_lambda * (0:M-1). * cos(el) * cos(az)); end function a steer_vec_z(el, M, d_lambda) a exp(1j * 2 * pi * d_lambda * (0:M-1). * sin(el)); end这里用exp(j·k·r)的约定,后面所有相位提取都按这个符号走,别中途换约定。d_lambda取0.5是经典配置,一是满足空间采样定理,二是相位范围刚好落在[-π,π]内,angle()取相位时不会卷绕。2.3 两个子阵各自估计,配对才是主角x轴子阵的极点和z轴子阵的极点分别对应两个投影量:x轴极点相位 ∝ cos(el)·cos(az),方位和俯仰耦合在一起;z轴极点相位 ∝ sin(el),只含俯仰。问题是两组极点各自独立解出来,谁和谁属于同一个信源并不知道。这就是L阵列做2D DOA的核心难点,叫做配对。常见配对思路有三类:用每个分量的复幅度做最近邻、用特征向量联合对角化、直接枚举配对后做数据重构残差检验。第三种最稳,后面第5章会给完整实现。3. 用Matrix Pencil做2D DOA估计可复现的MATLAB实现3.1 矩阵束的一句话直觉Matrix Pencil的核心思想是:把一段观测序列看成若干个复指数分量的叠加,每个分量对应一个极点z_k。如果能从数据里分离出这些极点,角度就藏在极点相位里。做法是把观测序列重排成一个Hankel矩阵,利用Hankel矩阵相邻两行的平移关系,构造一个矩阵束,它的广义特征值就是极点。对阵列数据来说,每个阵元的观测就相当于一个采样点,阵元间距就是采样间隔。所以同一个pencil框架可以直接搬到阵列域:一个子阵的快拍数据,构造Hankel矩阵,求极点,极点的相位就是相邻阵元间的波程相位增量。这也是为什么pencil在单快拍下依然能工作——它靠的是Hankel矩阵的秩结构,而不是多快拍的统计平均。3.2 完整可跑的最小示例:仿真数据生成先定场景:8元子阵、2个不相关窄带信号、100快拍、信噪比15dB。clear; clc; rng(2024); M 8; % 每个子阵的阵元数(含原点,原点共用) d_lambda 0.5; % 阵元间距/波长 N 100; % 快拍数 K 2; % 信源数 snr_db 15; az_true [30, -20] * pi/180; % 方位角, 注意线阵对称模糊,后面单独讲 el_true [25, 40] * pi/180; % 俯仰角 S (randn(K, N) 1j*randn(K, N)) / sqrt(2); % 复高斯信源, 不相关 Xx steer_vec_x(az_true(1), el_true(1), M, d_lambda) * S(1,:) ... steer_vec_x(az_true(2), el_true(2), M, d_lambda) * S(2,:); Xz steer_vec_z(el_true(1), M, d_lambda) * S(1,:) ... steer_vec_z(el_true(2), M, d_lambda) * S(2,:); % 加复高斯白噪声 noise_power 10^(-snr_db/10); Xx Xx sqrt(noise_power/2)*(randn(M,N) 1j*randn(M,N)); Xz Xz sqrt(noise_power/2)*(randn(M,N) 1j*randn(M,N));逻辑说明:信源矩阵S是K×N,每一行是一个信源的复包络。两个方向向量分别作用在S的两行上,叠加后得到两个子阵的观测。注意这里方向向量用的是第一节的约定,信源之间不相关,pencil不需要协方差矩阵,直接吃原始快拍数据。参数说明:rng(2024)固定随机种子,保证每次跑出来的图一致;如果换信噪比,只需要改snr_db;信源数K要事先给准,这是所有子空间类方法的前提。3.3 核心函数:matrix_pencil_poles这个函数输入某个子阵的M×N快拍数据,输出K个极点。关键点在于多快拍的处理方式——对每个快拍构造Hankel块再水平拼接,拼接后取左奇异向量U而不是右奇异向量V。function poles matrix_pencil_poles(X, K, L) % X: M x N 子阵接收数据 % K: 信源个数(用于SVD截断) % L: pencil参数, 建议取 M/3 ~ M/2 之间 % 输出: 1 x K 极点, 模接近1, 相位含角度信息 [M, N] size(X); m M - L 1; % Hankel矩阵行数 Xe zeros(m, L * N); % 增强矩阵: 每个快拍的Hankel块水平拼接 for n 1:N x X(:, n); Hn zeros(m, L); for jj 1:L % 显式填充Hankel矩阵 Hn(:, jj) x(jj : jj m - 1); end Xe(:, (n-1)*L 1 : n*L) Hn; end [U, ~, ~] svd(Xe, econ); Us U(:, 1:K); % 信号子空间: 取左奇异向量 U1 Us(1:end-1, :); % 行方向平移不变性 U2 Us(2:end, :); Psi pinv(U1) * U2; % 旋转算子 poles eig(Psi).; end逻辑说明:Hankel矩阵的每一行是观测序列的一个窗口滑动,行索引对应阵元位移方向。拼接了N个快拍后,增强矩阵是m×(L·N),SVD得到的左奇异向量U的行数仍是m,行之间保留了平移结构,所以U1、U2的旋转关系成立。如果这时候去取V做平移,会得到L·N行,平移结构被多快拍打乱,解出来全是错的——这是实际实现里最常见的翻车点之一。参数说明:L就是pencil参数,也叫束参数。L太小时Hankel矩阵的列数不够,奇异值截断后子空间不稳;L太大时m变小,信号子空间的行数不足。经验区间是M/3到M/2,同时要满足L-1≥K且m-1≥K,否则旋转算子没有足够的行做最小二乘。我一般先用Lfloor(M/2)跑通,再扫一遍L验证稳定性。3.4 配对与角度解算对两个子阵分别调用matrix_pencil_poles,得到两组极点。配对用的是幅度一致性:同一个信源在x轴和z轴上的复包络幅度应该一致,先估计每个极点的幅度,枚举全排列选幅度差最小的组合。function [az_est, el_est] pair_by_amplitude(px, pz, Xx, Xz, d_lambda) % px, pz: 两个子阵估计出的极点(1 x K) % Xx, Xz: 两个子阵的观测数据, 用于幅度拟合 % 返回配对后的方位角和俯仰角(弧度) K numel(px); M size(Xx, 1); % 构造范德蒙德矩阵, 用最小二乘估计每个极点的复幅度 Zx zeros(M, K); Zz zeros(M, K); for k 1:K Zx(:, k) px(k).^(0:M-1).; Zz(:, k) pz(k).^(0:M-1).; end Rx (Zx * Zx) \ (Zx * Xx); % K x N, 每个快拍每个极点的复幅度 Rz (Zz * Zz) \ (Zz * Xz); ampx mean(abs(Rx), 2); % K x 1 ampz mean(abs(Rz), 2); % 枚举全排列, 让幅度差平方和最小 Perms perms(1:K); best_cost inf; best_perm []; for ii 1:size(Perms, 1) p Perms(ii, :); cost sum((ampx - ampz(p)).^2); if cost best_cost best_cost cost; best_perm p; end end % 按配对结果解算角度 el_est zeros(1, K); az_est zeros(1, K); for k 1:K el_est(k) asin(angle(pz(best_perm(k))) / (2*pi*d_lambda)); cos_el cos(el_est(k)); cos_az angle(px(k)) / (2*pi*d_lambda) / cos_el; az_est(k) acos(cos_az); end end逻辑说明:先从极点和数据里把幅度估计出来,幅度相近的极点就是同一个信源。这一步要求两个子阵在原点共用且幅相响应一致,实际工程中如果两个子阵的通道增益有差异,要先校准,否则幅度配对会系统性错配。角度解算时,俯仰角直接从z轴极点反解asin,方位角再用x轴极点的相位除以cos(俯仰)后反解acos。参数说明:K小时全排列枚举没问题,K超过6就不现实了,那时改用贪心或匈牙利算法做二分匹配。另外,幅度配对的精度受极点的条件数影响,如果两个信源的俯仰角很接近,z轴极点几乎重合,Zz矩阵接近奇异,幅度估计会炸,这种场景下配对应改用第5章的重构残差法。3.5 跑一遍并验证结果把上面两节拼起来,调用主流程:% 主流程 L floor(M/2); % pencil参数取4 px matrix_pencil_poles(Xx, K, L); pz matrix_pencil_poles(Xz, K, L); [az_est, el_est] pair_by_amplitude(px, pz, Xx, Xz, d_lambda); disp(真实值(角度):); disp([az_true; el_true] * 180/pi); disp(估计值(角度):); disp([az_est; el_est] * 180/pi);正常跑通时,打印出来的估计值应该和真实值接近,误差在0.5度以内(15dB、100快拍下)。如果对不上,先检查三件事:一是px和pz的极点模值是否都接近1,偏离太多说明L参数或信源数K不对;二是配对函数返回的配对结果,幅度差是否明显小于错配情况;三是snr_db是否太低,10dB以下幅度配对开始出现概率性错配。4. 常见问题排查与避坑:为什么你的结果总在抖动4.1 镜像模糊:方位角正负分不清现象:真实方位角是负30度,估计出来却是正30度,而且两个信源的角度估计相互交叉,时好时坏。原因:x轴均匀线阵对±az方向完全对称,cos(az)相同,投影信息丢失了符号。这跟算法无关,是阵列几何的固有问题,不是pencil能解决的。解决:三个办法,按成本排序。一,利用场景先验,比如雷达只扫前向,直接把acos主值映射到需要的象限;二,把x轴换成非均匀阵,打破对称性;三,在原点附近加一个错位的辅助阵元,用辅助阵元的相位差定符号。仿真阶段最省事的是先设az在0到90度之间,验证算法链路后再处理实际象限。4.2 pencil参数L乱选导致极点飞了现象:L取M-1时,估计角度剧烈抖动;L取2时,极点倒是稳定,但角度分辨率明显变差,两个靠近的信源分不开。原因:L决定Hankel矩阵的行数mM-L1和列数L。L太大,m太小,信号子空间行数不足,旋转矩阵的最小二乘解病态;L太小,Hankel矩阵的列数不够,奇异值截断后的子空间无法容纳足够的信息。解决:在L∈[3, M-2]范围内扫描,画估计角度的稳定性曲线,取一个平坦区间内的值。我常用的做法是Lfloor(M/2)起步,然后分别试Lfloor(M/3)和Lfloor(2M/3),三个结果一致就放心了。扫描时可以顺便监控eig(Psi)的模值,所有极点模值落在0.95到1.05之间才算正常。4.3 相干信源导致秩亏,极点少一个现象:两个信源是同一信号的多径,幅度相关,结果只估出一个角度,或者估计的角度是两个真实角的中间值。原因:pencil方法同样基于子空间,相干信源让Hankel增强矩阵的有效秩低于K,SVD截断后信号子空间塌缩,丢失信息。解决:做前后向平滑,把数据翻折拼接后再进pencil。在matrix_pencil_poles外面包一层:% 前后向平滑预处理: 对M x N数据 X_fb [X, conj(flipud(X))]; % M x 2N % 然后把X_fb传给matrix_pencil_poles逻辑说明:flipud把阵元顺序反转并取共轭,等同于把阵列的镜像观测补进来。相干信源在正反两个方向上的合成响应不同,拼接后协方差恢复满秩。代价是黑匣子式的自由度减半——等效快拍数翻倍,但对相干源的分辨能力有时反而不如加空间平滑彻底。工程上如果做的是多径场景,建议直接上空间平滑,每3到4个阵元一组做重叠子阵。4.4 多快拍拼接后取错奇异向量现象:单快拍数据跑MP结果挺正常,改成多快拍后输出全乱。很多人翻论文里单快拍的公式,把取V那一步原样搬到多快拍代码里,结果彻底翻车。原因:单快拍时,Hankel矩阵的右奇异向量V的行数等于L,行方向对应序列的滑动窗口,平移结构在V上成立。但多快拍水平拼接后,V的行数变成L×N,不同快拍的列混在一起,平移结构不复存在。此时只有左奇异向量U的行数仍是m,保住了阵元位移结构,取U才正确。解决:严格按第3章的写法——先对每个快拍构造Hankel块,水平拼接成增强矩阵,然后取左奇异向量U做U1/U2旋转。排查时打印size(Us),确认行数等于M-L1;如果是L或L×N,说明取错了。4.5 俯仰角接近90度或d/lambda太大,角度爆炸现象:其他角度的结果都正常,唯独某一组俯仰角接近90度时,方位角误差突然变成十几度,甚至出现相位跳变。原因:俯仰角接近90度时,来波方向接近天顶,x轴子阵的投影cos(el)·cos(az)趋近于0,几乎不含方位信息,噪声被cos(el)除法放大。另外,如果d/lambda超过0.5,相邻阵元相位差可能超过±π,angle()取相位时发生卷绕,角度直接跳变。解决:一是把d/lambda控制在0.4到0.5之间,别为了增大孔径硬上0.6;二是对解算结果做合理性检查,当cos(el)小于0.1时,直接标记该信源方位不可估,不输出数值;三是如果应用场景确实有高俯仰角目标,考虑改用双L阵列或加一条y轴,用y轴子阵补方位信息。仿真时想绕开这个坑,把el_true设在10到70度之间。5. 让配对更稳的重构残差法:一个可以直接替换的进阶技巧幅度配对在15dB以上、极点分离度好的时候够用,但信噪比掉到5dB,或者两个信源的俯仰角只差5度时,幅度估计本身就被噪声污染,幅度差最小准则开始出错。重构残差法的思路完全不同:对每一种候选配对,直接用估计出的角度构造联合方向矩阵,用最小二乘反演出信源波形,再用这个波形重构整个双轴观测数据,计算重构误差。正确配对必须能同时解释两条子阵的数据,残差一定最小。这个准则不依赖幅度估计,天然免疫幅度噪声。function [az_best, el_best] pair_by_reconstruction(px, pz, Xx, Xz, d_lambda) % 枚举所有配对, 用联合重构残差选最优 % 适合低信噪比或俯仰角接近的场景 K numel(px); [M, ~] size(Xx); D [Xx; Xz]; % 2M x N 联合观测 Perms perms(1:K); best_res inf; az_best zeros(1, K); el_best zeros(1, K); for ii 1:size(Perms, 1) p Perms(ii, :); az_tmp zeros(1, K); el_tmp zeros(1, K); for k 1:K el_tmp(k) asin(angle(pz(p(k))) / (2*pi*d_lambda)); cos_az angle(px(k)) / (2*pi*d_lambda) / cos(el_tmp(k)); az_tmp(k) acos(cos_az); end % 构造联合方向矩阵: 上半块x轴, 下半块z轴 A zeros(2*M, K); for k 1:K A(1:M, k) exp(1j*2*pi*d_lambda*(0:M-1).*cos(el_tmp(k))*cos(az_tmp(k))); A(M1:2*M, k) exp(1j*2*pi*d_lambda*(0:M-1).*sin(el_tmp(k))); end S_hat A \ D; % K x N 最小二乘信源波形 res norm(D - A * S_hat, fro); if res best_res best_res res; az_best az_tmp; el_best el_tmp; end end end逻辑说明:关键在于A和D的拼接方式,同一个信源在x轴和z轴的复振幅必须一致,这是联合重构能够成立的物理前提。S_hat是一次性从联合观测里反演出来的,它同时解释了x轴和z轴的数据,残差res越小说明这个配对与两个子阵的观测越自洽。枚举全排列时,K4以内运算量都可以接受;K更大时改成随机采样候选配对,或者先用幅度配对缩小候选集,再用重构残差决赛。一个小技巧:拿到重构残差最小的配对后,可以在az和el附近各做±2度的小网格搜索,把每个角度组合再代进重构残差,取最小值对应的角度作为最终输出。这就相当于做了一次简化版的最大似然精估计,通常能把误差再压下去0.2到0.3度,代价只是多跑几十次残差计算。我早期做L阵列的时候,一直用幅度配对,觉得简单直接,结果项目里遇到两个目标俯仰角只差4度的场景,输出错配率达到三成,当时以为是阵列标定问题,折腾了快两周。后来换成重构残差法,同一批数据错配率直接归零。现在我的习惯是仿真阶段两个方法都跑,幅度配对出结果,重构残差做仲裁,两个一致才认为这组角度可信。这个习惯帮我挡掉了不少低信噪比下的隐性错配,希望也能帮到你。本文还有配套的精品资源点击获取