有理阵列DOA估计实战：从阵列设计到MUSIC谱估计的完整复现

阵列信号处理这个圈子，DOA估计是绕不开的经典命题。我最早接触这块是在做一套声源定位的小系统，当时用的是均匀线阵，半波长间距，教科书式的配置，跑起来确实稳。但后来场景一变，要求阵列孔径尽量大、阵元数量又不能太多，均匀阵的局限性一下子就暴露出来了——稀疏布阵会带来栅瓣和模糊，均匀密布又成本吃不消。正是在这个卡点上，我开始认真读稀疏阵列、非均匀阵列这条线的文献，其中有一篇让我印象很深，就是关于有理阵列用于DOA估计的工作。这篇东西的价值不在于提出了多么复杂的算法，而在于它用一种非常"工程师思维"的方式，把阵列几何的设计和估计算法解耦开来，给出了一个可计算、可复现的框架。下面我就按自己读论文、复现代码、踩坑调试的完整路径，把这篇论文的核心内容掰开揉碎讲一遍，适合做阵列信号处理、雷达/声呐/无线定位方向的朋友参考，也适合刚入门DOA、想搞清楚"阵列几何到底怎么影响估计性能"的读者。

1. 先搞清楚这篇论文到底在解决什么问题

1.1 从均匀线阵的天花板说起

要理解有理阵列的动机，得先明白均匀线阵（ULA）的痛点在哪。ULA的阵元间距通常取半波长，这个选择不是随便定的——它是空间采样的奈奎斯特条件。你可以把阵列想象成对空间中的来波方向做"采样"，阵元间距就是采样间隔，来波方向的正弦值就是被采样的信号。间距超过半波长，空间采样率不够，就会出现栅瓣，也就是在真实方向之外冒出一堆虚假的峰值，DOA估计直接失效。

问题在于，半波长的约束把阵列孔径和阵元数量死死绑在一起。想要大孔径提升分辨率，就得增加阵元数；想省阵元，孔径就上不去。这在工程上很要命，尤其是天线成本高、安装空间受限的场景。稀疏阵列的思路就是打破这个绑定：阵元间距拉大，但通过非均匀的排布方式，让阵列的"差分协阵"（difference co-array）依然保持连续、无空洞，从而在不增加物理阵元的前提下扩展有效孔径。

有理阵列属于稀疏阵列家族里比较优雅的一支。它的核心思想是：阵元位置不用整数倍半波长这种离散网格来定义，而是用有理数倍的基准间距来刻画，通过精心选择这些有理数，让差分协阵覆盖一段连续的整数格点。这个"有理"二字不是数学上的花架子，它直接决定了阵列的稀疏模式能不能用有限个参数描述清楚，进而决定了后续估计算法能不能高效实现。

1.2 有理阵列的核心承诺与适用边界

论文给出的核心承诺有两层。第一层是几何层面：有理阵列能够以远少于ULA的阵元数，构造出孔径大得多的连续差分协阵，这意味着可辨识的信源数上限（由连续协阵的自由度决定）显著提升。第二层是算法层面：由于阵元位置是有理数结构，协阵的构造和后续的谱估计可以借助数论和多项式工具来组织，避免了纯暴力搜索。

但这里必须说清楚适用边界，这也是我读论文时反复提醒自己的地方。有理阵列的优势建立在"差分协阵连续"这个前提上，一旦设计不当出现空洞，自由度就打折扣，甚至比ULA还差。另外，它假设的是窄带远场信源模型，信源之间不相干或者相干性可控，噪声是空间白噪声。实际场景里如果存在严重的多径相干、阵列校准误差、互耦效应，有理阵列的稀疏结构反而可能放大这些误差的影响——稀疏阵对校准精度的敏感度普遍高于密布阵，这一点后面会专门讲。

1.3 为什么值得花时间读这篇

市面上讲稀疏阵列的文献不少，嵌套阵、互质阵的教程也很多。这篇论文的差异化在于它把"有理"这个约束用得很彻底，给出了阵列位置参数和协阵连续性之间的显式关系，而不是像很多工程文献那样只给一个搜索算法然后甩一堆仿真图。对我这种喜欢从原理往下推的人来说，它的推导链条是完整的：从有理数的选取规则，到差分协阵的整数格点覆盖证明，再到基于协阵的MUSIC类谱估计，每一环都能对上。复现的时候你会有一种"每一步都知道为什么"的踏实感，而不是调参调到怀疑人生。

2. 有理阵列的几何设计与数论基础

2.1 阵元位置的有理数表示

有理阵列的阵元位置定义为基准间距d的有理数倍。设基准间距取d = λ/2（半波长），阵元位置集合可以写成：

S = { q_i * d | q_i ∈ Q, i = 0, 1, ..., N-1 }

其中Q是有理数集合，q_i通常表示为p_i / r_i的形式，p_i、r_i为整数。论文里为了便于分析，往往把所有q_i统一到同一个分母下，比如取公共分母R，则q_i = m_i / R，m_i为整数。这样一来，物理位置就映射到了一个整数格点m_i上，只不过这个格点的"单位"是d/R，比半波长细了R倍。

这个统一分母的操作是理解全文的关键。它把"有理数位置"这个看似抽象的概念，转化成了"在一个更细的整数网格上挑点"的问题。你可以类比成：原本ULA是在厘米刻度上等间距摆点，有理阵列是在毫米刻度上非均匀挑点，但最终物理间距还是厘米量级。细网格给了设计自由度，粗间距保证了稀疏性。

2.2 差分协阵与连续自由度

差分协阵是稀疏阵列理论的基石。给定阵元位置集合S，差分协阵定义为所有阵元位置两两之差构成的集合（含自身差，即零）：

D = { q_i - q_j | q_i, q_j ∈ S }

DOA估计的自由度取决于D中有多少个连续的唯一整数格点。为什么是差分？因为阵列接收信号的协方差矩阵，其元素本质上对应的是阵元对的相位差，而相位差正比于位置差。协方差矩阵经过向量化（vectorization）之后，等效于一个"虚拟阵列"的接收数据，这个虚拟阵列的阵元位置就是差分协阵D。所以差分协阵连续段越长，等效虚拟孔径越大，能分辨的信源数越多、分辨率越高。

ULA的差分协阵是连续的，从-(N-1)到(N-1)，自由度2N-1。有理阵列的目标就是：用N个物理阵元，构造出比2N-1长得多的连续差分协阵。论文给出的设计规则，本质上是在整数格点m_i的选取上做文章，让所有可能的差值m_i - m_j尽可能密集地覆盖一段连续整数区间。

2.3 有理数选取的约束条件

这里是最容易翻车的地方。不是随便挑一组有理数就能得到连续协阵。论文给出的约束大致可以归纳为几条：

• 公共分母R的选择要与阵元数匹配：R太小，细网格不够细，挑不出足够的稀疏模式；R太大，格点太密，物理间距趋近均匀，稀疏性丧失，而且数值条件变差。
• 整数格点集合M = {m_i}必须满足差集覆盖：即{m_i - m_j}要包含一段连续的整数区间[-(L), L]，L越大越好。这是一个组合设计问题，论文用了一些构造性的方法（类似差集、Singer差集的思想）来保证覆盖。
• 端点阵元要撑开孔径：通常把最小和最大的m_i放在区间两端，中间的点用来填差集的空洞。

我复现时踩的第一个坑就在这里。一开始我随手取了一组看起来"挺稀疏"的有理数，结果画出差分协阵一看，中间断了好几个格点，自由度还不如同阵元数的ULA。后来老老实实按论文的构造规则来，先定R，再按差集覆盖的要求反推m_i，才把连续段撑起来。这个教训很直接：有理阵列的设计是"先算后摆"，不能凭感觉。

2.4 一个可手算的小例子

为了让你有直观感受，我构造一个简化例子（数值做了调整，便于手算，不代表论文原始参数）。取R = 4，阵元数N = 4，整数格点选M = {0, 1, 4, 6}。物理位置就是{0, 0.25d, d, 1.5d}。

算一下差分集合：所有m_i - m_j，i, j遍历。列出非负差值：0, 1, 2, 3, 4, 5, 6。你会发现从0到6全部覆盖，连续段长度7，对应自由度（双边）约13。而同样4阵元的ULA，连续差分只有0到3，自由度7。也就是说，这个小小的有理阵列用4个阵元拿到了接近翻倍的连续自由度。当然这个例子很粗糙，实际设计中R和M的选取要满足更严格的覆盖证明，但它足以说明"有理+稀疏"的威力在哪。

注意：上面这个例子只用于说明差分协阵的覆盖逻辑，实际工程中R的取值、M的构造需要按论文的差集条件严格验证，切勿直接套用。

3. 从协阵到DOA谱估计的完整链路

3.1 接收信号模型与协方差矩阵

先把信号模型摆清楚。N个阵元，K个远场窄带信源，来波方向θ_k，接收信号写成：

x(t) = A(θ) s(t) + n(t)

A(θ)是N×K的阵列流形矩阵，第(i,k)个元素是exp(-j2π(q_i*d/λ)*sin(θ_k))。注意这里的q_i就是有理数位置系数，这是有理阵列和ULA在模型上唯一的区别——流形矩阵的相位项带上了q_i这个非整数因子。

协方差矩阵R_x = E[x x^H]，实际中用T个快拍的样本协方差估计：R_x ≈ (1/T) Σ x(t) x(t)^H。这一步和ULA完全一样，没有任何特殊处理。论文在这里没有玩花样，这点我很欣赏——它把复杂度都放在了阵列几何设计上，信号处理前端保持标准。

3.2 向量化与虚拟阵列构造

关键一步是向量化。把R_x按列堆叠成vec(R_x)，利用Kronecker积的性质，可以得到：

vec(R_x) = (A* ⊙ A) p + σ² vec(I)

其中⊙是Khatri-Rao积，p是信源功率向量，σ²是噪声功率。A* ⊙ A的每一列对应一个"虚拟阵元"的导向矢量，其位置正是差分协阵的元素。这就是为什么差分协阵决定了自由度——虚拟阵列的孔径就是差分协阵的跨度。

实际操作中，vec(R_x)里对应重复差分位置的元素会有多个，需要做去重和平均（类似空间平滑的预处理）。论文的处理方式是提取连续差分段对应的那部分元素，重排成一个等效的ULA接收向量，然后对这个等效向量跑标准MUSIC。这个"重排"步骤是代码实现里最容易出错的地方，索引映射必须精确，错一位整个谱就废了。

3.3 等效ULA上的MUSIC谱估计

拿到连续差分段对应的等效接收向量y之后，构造等效协方差R_y = y y^H（单快拍形式，实际中配合空间平滑恢复秩），做特征分解，取噪声子空间E_n，谱函数：

P(θ) = 1 / (a_v(θ)^H E_n E_n^H a_v(θ))

a_v(θ)是等效虚拟ULA的导向矢量，间距为d/R（细网格单位）。扫描θ，找峰值即为DOA估计。

这里有个细节值得说：由于等效虚拟阵列的"阵元间距"是d/R，它的无模糊视场对应的正弦范围是±(R/2)左右（按细网格奈奎斯特），实际物理来波方向仍限制在±90度，所以扫描范围不用变，但谱的分辨率因为虚拟孔径变大而显著提升。我实测下来，同样4个物理阵元，有理阵列的谱峰宽度明显比ULA窄，两个相近角度的信源能分开，ULA则糊成一团。

3.4 快拍数与相干性的处理

单快拍向量化带来的问题是等效协方差秩亏，K个信源秩最多为1（如果不做处理）。标准解法是空间平滑（spatial smoothing）：把等效ULA再切成若干重叠子阵，对各子阵协方差求平均，恢复秩到子阵数。子阵长度的选择要在"秩恢复能力"和"孔径损失"之间权衡——子阵越长孔径越大但子阵数少、秩恢复弱；子阵越短秩恢复强但孔径缩水。

论文对这块给了指导但不算特别细，我复现时是自己按经验调的：子阵长度取等效ULA长度的一半左右，子阵数取K+1以上，保证秩够。快拍数T方面，样本协方差估计误差随T增大而减小，但有理阵列因为稀疏，协方差矩阵的条件数比ULA差，T太小的时候估计误差会被放大。我一般取T不低于阵元数的10倍，信噪比低的时候还要再加。

4. 实操复现：从参数设计到谱峰输出

4.1 环境与依赖准备

我用的是Python，核心依赖就三个：numpy做矩阵运算，scipy做特征分解和峰值检测，matplotlib画谱。没有用任何专门的阵列处理库，因为这篇论文的流程足够透明，自己写反而可控。

pip install numpy scipy matplotlib

代码结构我分成四个模块：阵列几何生成、信号仿真、协阵重排、MUSIC谱计算。分开写的好处是每步都能单独验证，出问题好定位。

4.2 阵列几何参数的设计流程

按论文思路，参数设计顺序是：先定阵元数N和目标连续自由度L，再选公共分母R，最后构造整数格点M。

我实际用的流程是这样的：先根据应用需求确定要分辨的最大信源数K_max，连续自由度至少要2K_max，反推需要的连续差分长度L ≥ K_max。然后选R，经验上R取N到2N之间比较平衡，太小稀疏性不够，太大数值条件差。最后用贪心+回溯的方式构造M：从{0}开始，每次加入一个格点，检查当前差集覆盖是否扩展，直到连续段达到L或阵元数用完。

这个构造过程有点组合搜索的味道，N小的时候（比如8以内）可以暴力枚举验证，N大了就得用论文给的构造性规则。我建议复现时先用小N把流程跑通，确认差分协阵画出来是连续的，再放大N。

4.3 信号仿真与协方差估计

仿真部分要注意信源角度别取在栅瓣敏感区。我一般取[-60°, 60°]范围内的随机角度，间隔至少大于理论分辨率的一半，避免两个源太近导致谱峰合并、看不出算法好坏。

import numpy as np def simulate_signal(M_positions, angles_deg, snr_db, T, d_over_lambda=0.5): N = len(M_positions) K = len(angles_deg) theta = np.deg2rad(angles_deg) # 阵列流形，位置为有理数倍d A = np.exp(-1j * 2 * np.pi * d_over_lambda * np.outer(M_positions, np.sin(theta))) S = (np.random.randn(K, T) + 1j * np.random.randn(K, T)) / np.sqrt(2) noise = (np.random.randn(N, T) + 1j * np.random.randn(N, T)) / np.sqrt(2) sig_power = 10 ** (snr_db / 10) X = A @ S * np.sqrt(sig_power) + noise return X, A

这段里M_positions就是q_i（有理数系数），注意流形用的是q_i而不是整数索引，这是有理阵列和ULA代码唯一的实质差别。SNR按单阵元信噪比定义，噪声功率归一化为1。

4.4 协阵重排的索引映射实现

这是整个复现里最烧脑的一段。思路是：先算出所有差分值q_i - q_j，统一到细网格整数（乘以R取整），建立"差分值 → 协方差矩阵元素索引"的映射字典，然后提取连续段，按顺序排成等效接收向量。

def coarray_rearrange(Rx, q, R): N = len(q) m = np.round(np.array(q) * R).astype(int) diff_map = {} for i in range(N): for j in range(N): diff = m[i] - m[j] diff_map.setdefault(diff, []).append((i, j)) # 找最长连续段 keys = sorted(diff_map.keys()) best_start, best_len, cur_start, cur_len = keys[0], 1, keys[0], 1 for k in keys[1:]: if k == cur_start + cur_len: cur_len += 1 else: if cur_len > best_len: best_start, best_len = cur_start, cur_len cur_start, cur_len = k, 1 if cur_len > best_len: best_start, best_len = cur_start, cur_len # 提取并平均重复项 y = [] for d in range(best_start, best_start + best_len): idxs = diff_map[d] val = np.mean([Rx[i, j] for (i, j) in idxs]) y.append(val) return np.array(y), best_start, best_len

这段代码里"平均重复项"很关键。同一个差分值可能对应多对(i,j)，它们的协方差元素理论上相等（理想无噪），实际有噪时取平均能降方差。我一开始偷懒只取第一个，结果谱底噪明显偏高，改成平均之后干净不少。

4.5 MUSIC谱计算与峰值提取

拿到等效向量y后，做空间平滑恢复秩，再特征分解。平滑子阵长度我取len(y)的2/3左右，子阵数保证大于信源数。

def music_spectrum(y, num_sources, sub_len, scan_deg, R): L = len(y) num_sub = L - sub_len + 1 Rsm = np.zeros((sub_len, sub_len), dtype=complex) for s in range(num_sub): seg = y[s:s+sub_len].reshape(-1, 1) Rsm += seg @ seg.conj().T Rsm /= num_sub eigvals, eigvecs = np.linalg.eigh(Rsm) En = eigvecs[:, :sub_len - num_sources] # 等效虚拟ULA间距为d/R a = np.exp(-1j * 2 * np.pi * (1.0 / R) * np.outer(np.arange(sub_len), np.sin(np.deg2rad(scan_deg)))) P = 1.0 / np.sum(np.abs(En.conj().T @ a) ** 2, axis=0) return P

注意导向矢量里的1.0/R，这是等效虚拟阵元的间距系数，写错这个谱峰位置会整体偏移，我第一次跑就栽在这，谱峰全歪了，排查半天才发现是间距单位没换算。

4.6 一次完整运行的参数记录

我记录一组实测参数供参考：N=6，R=8，M={0,1,3,6,10,13}（这组是我按差集覆盖调出来的，连续段覆盖0到13），信源2个，角度-20°和15°，SNR=10dB，快拍T=200。跑出来的谱在两个角度处峰值清晰，谱峰宽度比同阵元ULA窄约40%，两源间隔35°轻松分开。把SNR降到0dB，谱峰还在但底噪抬高，峰值位置基本不偏，说明有理阵列在中低信噪比下鲁棒性尚可，但前提是阵列校准准确。

5. 踩坑实录与常见问题排查

5.1 差分协阵不连续导致自由度崩塌

这是最高频的问题。表现是谱估计能跑，但分辨率远低于预期，或者信源数一多就失效。根因是M的选取没满足差集覆盖，连续段中间有洞。排查方法很直接：把diff_map的键排序打印出来，看连续段到底多长，和理论值对比。我建议每次改M之后都画一张差分协阵覆盖图，横轴差分值、纵轴是否被覆盖，一眼就能看出洞在哪。

修复思路是调整M里的格点，优先补洞。比如连续段在某个值断了，就往M里加一个能产生该差值的格点，或者挪动现有格点。这个过程手工调比较费劲，N大了建议写个小搜索脚本自动补洞。

5.2 谱峰偏移与间距单位错误

前面提过，等效虚拟ULA的间距是d/R，导向矢量里必须体现这个1/R因子。如果误用1.0（当成半波长间距），谱峰会整体压缩或拉伸，角度估计系统性偏差。这个错误很隐蔽，因为谱形看起来还是"正常"的，只是峰位置不对。排查技巧：用一个已知角度的单信源做标定，看估计角度和真实角度是否吻合，不吻合先查间距系数。

5.3 相干信源下的秩亏与平滑参数

两个信源高度相干时，等效协方差的秩恢复不充分，MUSIC谱会出现伪峰或者漏峰。这时候要加大空间平滑的子阵数（即减小子阵长度），牺牲一点孔径换秩。但子阵太短孔径损失大，分辨率下降，需要权衡。我的经验是子阵长度不低于等效ULA长度的1/3，子阵数至少K+1，相干越强子阵数越要多。

5.4 稀疏阵对校准误差的放大效应

这是有理阵列的固有短板，必须正视。稀疏布阵意味着阵元间距大，相位误差在阵元间的累积更明显，互耦、位置误差、通道不一致都会被放大。我做过对比：给ULA和有理阵列加同样的位置扰动（比如0.01λ的随机偏移），ULA谱峰基本不动，有理阵列的谱峰明显展宽、旁瓣抬高。所以工程上如果有理阵列要用，校准精度要求比ULA高一个档次，或者引入稳健算法（比如对导向矢量做不确定性约束）。

5.5 常见问题速查表

5.6 数值条件与分母R的取舍

R选得太大，细网格太密，等效虚拟阵列的导向矢量列之间相关性增强，协方差矩阵条件数变差，特征分解的数值稳定性下降。我实测R超过2N之后，谱底噪开始明显抬升，低信噪比下更糟。所以R不是越大越好，N到2N之间是比较稳的区间。这个取舍论文里提了但没给硬指标，属于要靠实测调的经验参数。

6. 我对这套方法的一点实际体会

复现完这篇论文的流程，我最大的感受是：有理阵列的精髓在"设计"而不在"算法"。它的信号处理链路和标准MUSIC几乎一样，复杂度全压在阵列几何的构造上——怎么选R、怎么排M、怎么保证差集连续，这些才是决定成败的地方。这也意味着，如果你只是拿一组现成的有理数位置去跑代码，很可能跑不出论文的效果，因为位置参数和你的阵元数、目标自由度不匹配。

另外一点体会是关于工程落地的现实感。有理阵列在仿真里很漂亮，自由度翻倍、分辨率提升，但一旦引入校准误差、互耦、多径相干，优势会被侵蚀。所以我在实际项目里用它的策略是：优先用在信源数多、孔径受限、且校准条件可控的场景，比如实验室里的声源定位、短距离高精度测向；对于野外、低成本、免校准要求高的场景，还是老老实实用ULA或者成熟的互质阵更稳妥。这不是否定有理阵列，而是任何稀疏阵都得认清自己的适用边界。

最后分享一个调试小技巧：每次改阵列参数，先别急着跑完整MUSIC，先画差分协阵覆盖图和等效虚拟阵列的导向矢量相关性热图。覆盖图确认连续性，相关性热图确认数值条件（列间相关性太高说明R太大或M太密）。这两张图没问题，再跑谱估计，基本一次就能出结果，省去大量盲调时间。这个习惯是我踩了无数次"谱不对但不知道哪错"的坑之后养成的，比事后排查高效得多。