简介互质面阵在相同阵元数下能提供更大的阵列孔径和更低的互耦影响为高精度二维DOA估计创造了条件在雷达探测、目标定位等场景中具有实用价值。对应的MATLAB实现方案面向信号处理与雷达专业学生及需要解决阵列测向模糊问题的研究人员代码涵盖互质面阵模型建立与基于二维MUSIC算法的解模糊处理采用参数化编程、步骤清晰、注释完整便于快速理解与二次开发整套代码按主流程、核心算子与结果分析分层组织既可直接运行也可修改参数进行算法对比。压缩包共4个文件包含3个m脚本和1张分析图m脚本分别对应主流程、Kronecker积与Khatri-Rao运算等关键模块分析图直观展示仿真结果及解模糊效果整体仅54KB轻量易用目前已有1739人学习下载。读者可通过源码与示意图系统掌握互质面阵信号建模、角度估计以及模糊消除的完整链路尤其适合课程设计、毕业设计等阵列信号处理场景。1. 互质面阵下的二维DOA估计为什么MUSIC会给出假方位做阵列信号处理的人第一次把MUSIC算法搬到互质面阵上十有八九会看到谱图上多出好几个长得一模一样的尖峰有的峰甚至比真实来波方向还亮。这不是MUSIC失效而是互质面阵的稀疏阵元间距触发了空间混叠也就是常说的模糊峰。互质面阵用两套间距互质的均匀子阵嵌套在一起用较小的物理孔径换来了远超阵元数目的自由度代价就是单一子阵的MUSIC谱没法直接读。本文要解决的就是这件事先把模糊产生的几何条件讲透再给出一套从双子阵估计、差集共阵重构到峰配对的去模糊流程最后用蒙特卡洛脚本验证结果可复现。2. 互质面阵的信号模型与二维MUSIC谱从子阵布阵到谱峰搜索2.1 互质面阵的两种常见布局与自由度收益互质面阵不是一种固定结构的阵列而是「两个均匀子阵嵌套/平移」的设计思路。最常见的布局是L型互质阵列x轴和y轴各放一组互质线阵共享原点处的参考阵元两组线阵分别负责方位角和俯仰角的信息采集最后通过配对得到二维角度。另一类常见布局是矩形互质面阵两个子阵本身都是稀疏均匀矩形栅格一个间距为Md另一个间距为NdM和N互质d取半波长。两种布局都能做二维DOA估计区别在于L型实现简单、配对麻烦矩形面阵的差集共阵更完整自由度更高但数据量和计算量也大。选择哪种布局取决于信号环境和硬件约束。如果你的阵元位置已经固化成矩形板那矩形互质面阵是顺理成章的选择如果只是想在现有均匀面阵基础上改造把原来的阵元抽稀成两套互质栅格L型方案改动更小。我一般建议先用矩形互质面阵做原理验证因为差集共阵的虚拟阵元是天然的矩形栅格后处理代码更直观。互质面阵的自由度收益其实是很多人的第一直觉陷阱。一个8阵元的均匀面阵最多分辨8个信号同样的8个物理阵元拆成44的互质子阵差集共阵能提供的虚拟阵元数远大于8理论上能同时估计十几个信号。当然这只是自由度上限实际还要看信噪比、快拍数和搜索网格密度。把方向向量当做「空间频率」来看互质结构等于在空间频率域做了非均匀采样再用两套采样网格的交集把模糊剔除。2.2 二维MUSIC谱的数学形式与网格搜索实现二维MUSIC的核心还是子空间分解。设阵列有P个物理阵元第k个快拍的接收数据为x(k) A s(k) n(k)其中A是P×D的阵列流形矩阵D是信号源数。协方差矩阵R E[x x^H]做特征分解后大特征值对应的特征向量张成信号子空间Es其余张成噪声子空间En。MUSIC谱函数写作P_music(θ, φ) 1 / (a^H(θ, φ) En En^H a(θ, φ))其中a(θ, φ)是来波方向为(θ, φ)时的导向矢量。注意二维MUSIC需要在θ和φ两个维度上都做网格扫描每个网格点计算一次谱值所以复杂度是O(Nθ × Nφ × P²)Nθ和Nφ分别是两个角度的网格点数。实际工程里常用先粗搜后精搜的两步策略先用5°步长找出谱峰大致位置再在峰附近用0.1°步长细化这样能把计算量降一到两个数量级。导向矢量的构造是容易被忽视的地方。矩形互质面阵的阵元坐标是两组整数对的集合子阵1的阵元坐标是{(mM, nN) d}子阵2是{(mN, nM) d}m和n取遍各自子阵的阵元序号。有了坐标集合导向矢量就是a(θ, φ) exp(j 2π/λ (x sinθ cosφ y sinθ sinφ))这里的θ是俯仰角φ是方位角。写成Python实现时务必把坐标单位统一到波长否则sinθ项会差一个倍数谱峰位置全部跑偏。2.3 用Python构造互质面阵并生成单快拍测试数据下面这段代码演示了如何构造一个M3、N5的矩形互质面阵并生成两个远场窄带信号的接收数据。代码里的坐标生成方式可以直接当成模板复用。import numpy as np def coprime_rect_array(M3, N5, d0.5): 构造矩形互质面阵的物理阵元坐标 M, N 为互质因子d 为半波长 返回坐标数组 shape(P, 2)单位是波长 sub1 [] sub2 [] # 子阵1x方向间距 M*dy方向间距 N*d for ix in range(N): for iy in range(M): sub1.append((ix * M * d, iy * N * d)) # 子阵2x方向间距 N*dy方向间距 M*d for ix in range(M): for iy in range(N): sub2.append((ix * N * d, iy * M * d)) # 去掉原点重复阵元 coords list(set(sub1) | set(sub2)) return np.array(coords) def steering_vector(coords, theta, phi, wavelength1.0): 计算某个来波方向的导向矢量 theta: 俯仰角 (0~90度)phi: 方位角 (0~360度) coords 单位是波长wavelength 归一化为1 k 2 * np.pi / wavelength x coords[:, 0] y coords[:, 1] phase k * (x * np.sin(theta) * np.cos(phi) y * np.sin(theta) * np.sin(phi)) return np.exp(1j * phase) # 生成阵列坐标P 是物理阵元数 coords coprime_rect_array(M3, N5) P coords.shape[0] print(物理阵元数:, P) # 两个真实来波方向 thetas [30, 50] # 俯仰角单位度 phis [45, 120] # 方位角单位度 D len(thetas) # 构造流形矩阵 A A np.zeros((P, D), dtypecomplex) for idx in range(D): A[:, idx] steering_vector(coords, np.deg2rad(thetas[idx]), np.deg2rad(phis[idx])) # 生成500个快拍信号是随机复包络加高斯白噪声 K 500 S np.random.randn(D, K) 1j * np.random.randn(D, K) X A S 0.1 * (np.random.randn(P, K) 1j * np.random.randn(P, K))这段代码里有三个关键点。第一子阵1的坐标是(ixMd, iyNd)子阵2是(ixNd, iyMd)这样x轴和y轴的间距正好互换了保证两个方向都能获得互质间距带来的自由度。第二用set去重是因为原点处的阵元(0,0)会被两个子阵同时包含真实硬件不可能放两个重叠的阵元去掉后物理阵元数会少于M×NN×M。第三信号用随机复包络生成快拍数K取500噪声方差0.1对应约20dB的输入信噪比这个设置在MUSIC里属于比较舒服的工作区间。有了X和坐标下一步就是把X的协方差矩阵做特征分解然后计算二维谱。这部分我会在第4章结合去模糊一起给出完整代码因为单子阵的MUSIC谱本来就是要被「怀疑」的中间产物。3. 模糊峰是怎么来的栅瓣周期与互质因子之间的关系3.1 空间混叠的数学条件均匀线阵的阵元间距如果大于半波长导向矢量在空间频率上会出现多个等价的峰值。原因不复杂导向矢量是exp(j 2π d sinθ / λ)当d超过λ/2时sinθ每变化一个周期相位就能绕回相同的值MUSIC谱就在多个θ位置出现同样的峰这些峰就是模糊峰。互质面阵的两个子阵间距分别是Md和Nd都远大于半波长所以每个子阵单独拿来做二维MUSIC谱面上必然有成片的模糊峰。模糊峰的位置不是随机的它服从栅瓣公式。以x轴方向的子阵为例设间距为L d则第k个栅瓣满足sinθ_k sinθ_true k λ / (L d)。L3时真实方向为30°的信号在sinθ域会额外出现对应k±1的两个栅瓣位置换算回角度就是大约70.5°和-30°如果允许负角度。二维情况下x和y两个方向独立产生栅瓣所以模糊峰在(θ, φ)平面上形成网格状分布网格间隔由M和N的数值决定。这就是为什么互质面阵不能直接套均匀阵列的MUSIC流程你没法判断哪个峰是真实的。有人会想「谱峰最高的那个总是真的吧」但在噪声影响下模糊峰的幅度和真实峰几乎一样因为导向矢量对这两个方向的响应是完全等价的。只有当信噪比非常高、快拍数非常大时真实峰才会因有限样本效应略微高一点拿这个当判据等于在赌运气完全不可靠。3.2 真实峰与模糊峰的判别特征从谱图上观察真实峰和模糊峰有一个几何上的区别真实峰在两个子阵的MUSIC谱里都出现且出现在同一个角度位置模糊峰虽然也出现在两个子阵的谱里但位置不重合。原因在于两个子阵的栅瓣周期不同子阵1的栅瓣间隔由M决定子阵2的间隔由N决定M和N互质所以两条栅瓣网格只在真实方向上有公共交点其余交点全部错开。这个性质就是「互质」名字的由来也是所有去模糊算法的出发点。你可以把它想象成两把刻度不同的尺子去量同一个长度两把尺子的刻度线只有在真实长度处才对齐其余地方对不齐。于是最朴素的去模糊方法就是分别对两个子阵做二维MUSIC谱各自找出前Q个峰值然后做峰配对找两组峰集中距离小于某个容差的组合剩下的峰值删除。实际操作中还要注意一个细节两个子阵的孔径不同角度分辨率也不同。子阵1的物理孔径大约是N×Md子阵2是M×Nd两者差异明显。孔径大的那个谱峰更锐利孔径小的谱峰更胖所以在做峰值匹配时容差不能设成一个固定值最好根据两个子阵各自的角度分辨率分别设定不然会把真实峰也滤掉。3.3 为什么不能用「子阵去模糊」替代「联合处理」有一种省事的思路既然互质面阵可以拆成两个均匀子阵那先让两个子阵各自做DOA估计然后取交集不就行了吗这个思路方向对但直接实现会踩到两个坑。第一个坑是子阵孔径太小分辨率不够。互质子阵的阵元间距大但阵元数量少阵列孔径通常比完整互质面阵的差集共阵小很多。两个信号角度接近时子阵MUSIC谱可能只有一个峰两个峰挨得太近没法分离。子阵估计给出的候选峰里压根没有真实方向后面取交集自然一场空。第二个坑是二维峰配对的组合爆炸。每个子阵谱提取20个峰两组峰做最近邻匹配最坏情况下要比较20×20个组合还要处理一个峰被多个峰匹配的歧义。工程上真正靠谱的做法是把两个子阵的数据放到一起做联合估计或者用差集共阵构造一个虚拟的均匀阵列。前者是第4章要讲的联合谱处理后者是利用差集共阵的虚拟阵元做传统MUSIC。两条路我都跑过差集共阵的路更稳因为它把互质结构彻底转化成一个虚拟均匀面阵后续处理完全不用关心模糊问题。4. 基于互质属性的去模糊实现MUSIC谱矫正的三步走4.1 第一步分别对两个子阵做二维MUSIC提取候选峰去模糊流程的第一步是把互质面阵拆回两个子阵各自独立地计算MUSIC谱并提取峰值。这里的峰值提取不能只取谱最大点而是要取局部极大值。二维谱峰提取可以用简单的滑窗比较某个网格点的谱值大于它周围8个邻居就记为一个候选峰然后按谱值从大到小排序取前Q个。from scipy.ndimage import maximum_filter def extract_peaks(spectrum, theta_grid, phi_grid, Q20): 从二维MUSIC谱中提取前Q个局部峰值 spectrum: (Ntheta, Nphi) 的谱值 theta_grid, phi_grid: 搜索网格单位是度 返回 (方位角列表, 俯仰角列表, 峰值列表) local_max maximum_filter(spectrum, size3) peak_mask (spectrum local_max) peak_indices np.argwhere(peak_mask) # 按谱值排序 peak_values spectrum[peak_mask] sorted_idx np.argsort(peak_values)[::-1][:Q] theta_peaks [] phi_peaks [] for idx in sorted_idx: r, c peak_indices[idx] theta_peaks.append(theta_grid[r]) phi_peaks.append(phi_grid[c]) return theta_peaks, phi_peaksmaximum_filter的size参数决定了峰与峰之间的最小间距。size3意味着两个相邻峰至少要隔一个网格点如果搜索步长是1°那两个角度差小于2°的峰会被合并成一个这可能把真实接近的两个信号吞掉。更稳妥的做法是先做一轮峰值提取再对每个峰周围做抛物线插值精化把峰值位置修正到亚网格精度。插值公式对二维谱可以分解成两个一维插值在θ方向上取峰附近的三个点做二次插值φ方向同样处理。Q的取值要大于信号源数D一般取D×5左右。取太少会把真实峰漏掉取太多则会引入噪声峰增加后续配对的误匹配率。我一般设D的上限为物理阵元数减1Q设成这个上限的两倍然后观察配对结果中是否有稳定的重复峰再调小Q。这个调参过程是手工的但做几次就有手感了。4.2 第二步差集共阵重构把互质结构变成虚拟均匀面阵差集共阵的思路是把互质面阵的非均匀物理阵元做自相关处理等效出一个阵元位置更密集的虚拟阵列。具体做法是对接收数据的协方差矩阵做向量化vec操作得到一个新的数据向量它的等价导向矢量由物理阵元的差集坐标给出。差集坐标集合定义为{x_i - x_j}其中x_i和x_j遍历所有物理阵元位置这个集合经过整理后会包含一个完整的均匀矩形栅格区域这个区域就是差集共阵的连续部分。def difference_coarray(coords, M, N, d0.5): 计算差集共阵的连续均匀栅格部分 返回虚拟阵元坐标列表以及对应的选择矩阵 # 所有两两差集坐标单位是d的整数倍 diff_set set() P coords.shape[0] for i in range(P): for j in range(P): dx int(round((coords[i, 0] - coords[j, 0]) / d)) dy int(round((coords[i, 1] - coords[j, 1]) / d)) diff_set.add((dx, dy)) # 找出能构成连续矩形栅格的最大矩形区域 dx_values [v[0] for v in diff_set] dy_values [v[1] for v in diff_set] dx_min, dx_max min(dx_values), max(dx_values) dy_min, dy_max min(dy_values), max(dy_values) virtual_coords [] for ix in range(dx_min, dx_max 1): for iy in range(dy_min, dy_max 1): if (ix, iy) in diff_set: virtual_coords.append((ix * d, iy * d)) return np.array(virtual_coords)虚拟阵列的连续栅格区域大小直接决定了可分辨信号数的上限。M3、N5时物理阵元大约有20个差集共阵的连续区域能达到一个较大的矩形栅格虚拟阵元数量通常超过80。这带来的直接好处是虚拟均匀面阵可以直接用传统MUSIC不需要再处理模糊峰因为虚拟阵元间距就是半波长的整数倍没有空间混叠。但差集共阵重构有一个代价向量化后的数据向量是单快拍的等效协方差矩阵的秩为1直接做特征分解只能得到一个非零特征值。解决办法是空间平滑把虚拟阵列划分成若干重叠子阵用子阵间的平移关系重建秩通常叫SS-MUSIC。空间平滑会把虚拟阵列的有效孔径缩小所以在平滑次数和孔径之间有个折中这个折中参数我会在第5章展开说。4.3 第三步谱峰配对与模糊峰剔除的完整代码最后一步是把两个子阵的候选峰做几何匹配保留两边的公共峰也就是真实方向。这里的关键是距离度量要放在方向余弦空间而不是直接比较角度。因为MUSIC谱栅瓣在sinθ域是等间隔的在角度域不是等间隔直接用角度差做匹配模糊峰也可能被误配进来。def disambiguate_by_pairing(theta1, phi1, theta2, phi2, tol0.03): 用方向余弦距离匹配两组峰返回匹配对 theta1/phi1: 子阵1的峰方向度 theta2/phi2: 子阵2的峰方向度 tol: 方向余弦距离容差 matched [] # 转方向余弦 cos1 np.cos(np.deg2rad(theta1)) * np.cos(np.deg2rad(phi1)) sin1 np.cos(np.deg2rad(theta1)) * np.sin(np.deg2rad(phi1)) cos2 np.cos(np.deg2rad(theta2)) * np.cos(np.deg2rad(phi2)) sin2 np.cos(np.deg2rad(theta2)) * np.sin(np.deg2rad(phi2)) for i in range(len(theta1)): dist np.sqrt((cos1[i] - cos2)**2 (sin1[i] - sin2)**2) min_idx np.argmin(dist) if dist[min_idx] tol: matched.append((theta1[i], phi1[i])) return matched配对完成后你得到的匹配峰数量应该等于真实信号数D因为两个子阵的栅瓣网格只在真实方向重合。但实际中噪声可能让某些假峰恰好落在容差范围内所以匹配结果里会出现多余峰。处理多余峰的办法不是直接删而是看匹配峰之间的谱值一致性真实方向在两个子阵里都有较高的谱值假峰至少有一边谱值偏低。可以设定一个谱值比门限比如真实峰在两个谱里的值都大于各自最大谱值的0.3倍否则剔除。这个门限需要根据实测噪声调0.2到0.4之间比较常见。最后一步是把匹配结果和差集共阵SS-MUSIC的结果做交叉验证。两种方法的物理原理不同误差来源也不同如果它们对同一信号输出角度差在1°以内基本可以确认这个方向是真的。我把这种交叉验证当成发布结果前的例行检查跑一次不到半分钟但能避免很多后续解释不清的问题。5. 互质面阵DOA估计避坑阵列几何、搜索步长与相干源的5个坑5.1 互质因子选得太小自由度优势直接归零现象M2、N3的互质面阵差集共阵的连续栅格区域很小虚拟阵元数量寥寥无几估计性能甚至不如同阵元数的均匀面阵。原因互质因子的乘积决定了差集共阵的连续区域范围。M和N太小差集集合里能形成连续矩形栅格的点不够多空间平滑后虚拟阵列孔径严重缩水。互质的好处是需要因子足够大才显现的M和N至少取3和5实际工程里我见过用5和7的。解决在设计阶段先算差集共阵的连续栅格尺寸不要只看物理阵元数。用上一章的difference_coarray函数把M和N遍历一遍找连续区域最大的组合。M3、N5和M4、N7是两组性价比不错的起步参数前者阵元少、适合原理验证后者孔径大、适合实际测向。5.2 二维搜索网格太粗MUSIC谱峰分裂成双峰现象搜索步长设为2°时谱峰顶部出现两个相邻的局部极大值峰值提取算法把它们当成两个峰配对时产生大量虚假匹配。原因MUSIC谱峰的真实宽度和阵列孔径成反比步长大于谱峰半宽时离散采样会落在峰两侧而不是峰顶形成双峰形态。这是离散化带来的伪峰不是真实的两个信号。解决先粗搜后精搜是成熟的规避方案。粗搜步长用5°定位峰的大致区域再对每个区域用0.1°步长重算谱。如果不想二次搜索可以给峰值提取算法加一个最小峰间距约束小于这个间距的相邻峰直接合并峰间距取预期角度分辨率的1/2。5.3 信号相干导致协方差矩阵秩亏现象两个来波信号完全相干比如同一个发射源的多径MUSIC谱只出一个峰或者两个峰幅度严重不对称去模糊后只恢复出一个方向。原因相干信号让协方差矩阵的秩降到1噪声子空间的维度判断出错MUSIC把两个信号当成一个来波。这不是互质结构的问题是所有子空间类算法的共同弱点。解决在信号模型中加入去相干处理常见做法有前向-后向空间平滑或直接在数据域用Toeplitz重构。差集共阵方法里空间平滑本身就有去相干作用所以相干源场景优先走SS-MUSIC路径不要用双峰配对路径。5.4 快拍数不足子空间泄漏导致谱峰偏移现象快拍数只有50时谱峰位置随机偏移约1°~2°多次实验的估计方差大甚至出现模糊峰比真实峰更稳定的假象。原因有限快拍让协方差矩阵的估计有误差特征向量不再精确张成信号子空间MUSIC谱峰会向噪声方向偏移。快拍越少偏移越随机而这种随机性对峰值配对的影响比绝对误差更大。解决快拍数低于100时先用差集共阵重构再对虚拟阵列做多次采样平均还能用前后向平均来加倍有效快拍数。如果实时性不允许攒快拍就一定要在配对容差里留出偏移余量方向余弦容差建议从0.03放宽到0.05。5.5 峰值配对时容差设置不当真实峰被误删现象配对容差设为0.01结果两个子阵的真实峰因为角度估计误差超出容差匹配失败最终输出结果里直接少了目标。原因容差写得太严苛没有考虑两个子阵孔径不同导致的分辨率差异。孔径小的子阵谱峰胖估计误差大真实方向可能偏离理论值好几度。解决把容差设成与两个子阵各自的角度分辨率线性相关。以M3、N5为例孔径小的子阵在30°俯仰角处的3dB谱宽约为6°容差取谱宽的一半比较安全。经验法则方向余弦空间容差0.05对应大约3°的角度误差足够覆盖中等信噪比下的估计偏差。6. 验证去模糊效果的一套指标与蒙特卡洛仿真脚本去模糊做得到底行不行不能靠肉眼数峰。三个指标最有说服力RMSE均方根误差、检测成功率和模糊峰残留率。RMSE衡量估计精度检测成功率衡量算法稳定性模糊峰残留率专门看去模糊是否彻底。三者合在一起能同时看出你的方案「偏不偏」「丢不丢目标」「干不干净」。def monte_carlo_doa(coords, M, N, trials200, snr_db20): 蒙特卡洛验证互质面阵去模糊性能 返回: 两个方向的RMSE、检测成功率、模糊峰残留率 true_theta, true_phi 30, 45 rmse_list [] success 0 residual_fog 0 for _ in range(trials): # 生成单信号数据信噪比由噪声方差控制 noise_var 10 ** (-snr_db / 10) A steering_vector(coords, np.deg2rad(true_theta), np.deg2rad(true_phi)) X A (np.random.randn(1, 500) 1j * np.random.randn(1, 500)) X np.sqrt(noise_var) * (np.random.randn(coords.shape[0], 500) 1j * np.random.randn(coords.shape[0], 500)) # 两个子阵分别估计配对去模糊 theta_hat, phi_hat doa_estimate(X, coords, M, N) if len(theta_hat) 1: err np.sqrt((theta_hat[0] - true_theta)**2 (phi_hat[0] - true_phi)**2) rmse_list.append(err) if err 1.0: success 1 else: residual_fog 1 rmse np.sqrt(np.mean(np.array(rmse_list)**2)) return rmse, success / trials, residual_fog / trials每次蒙特卡洛试验里doa_estimate返回的去模糊结果可能不止一个峰。只要输出峰数量大于1就计一次模糊残留说明该次去模糊没清理干净。正常参数下残留率应该低于5%如果残留率高于10%优先检查配对容差和峰值提取数量大概率是假峰被当成真实目标输出了。我自己的经验是把蒙特卡洛试验当成一种「回归测试」每改一次配对算法就跑一遍200次试验对比RMSE和残留率的变化。有一次我把容差从0.05改成0.03检测成功率从92%掉到71%但残留率也降了这才意识到容差是在「漏检」和「误检」之间做权衡没有绝对正确的值。现在我的做法是同时输出这两个指标用Pareto前沿的方式选择参数而不是单独盯着某一个数。还有一个便宜好用的验证技巧把谱峰图和配对结果画在同一张图上用圆点标出真实方向用叉号标出算法输出。目视检查能捕捉到很多指标发现不了的问题比如谱峰合并、栅瓣与真实峰在同一条等值线上等。这套验证流程跑通之后再去面对不同信噪比、不同来波数目的场景心里就有底了。希望帮到你。本文还有配套的精品资源点击获取