做了这么多年阵列信号处理每次带新人入门我最推荐的章节就是《阵列信号处理及MATLAB实现》里关于阵列响应矩阵的这部分。不管是做DOA估计还是波束形成最后都要落到阵列流型矩阵的构造上而不少同学恰恰是在这里开始出问题——线阵还算顺利一换成圆阵、L型或任意阵代码就开始乱。这篇文章把我自己在这个项目里反复验证过的实现思路、MATLAB代码和踩坑记录整理出来覆盖均匀线阵、均匀圆阵、L型阵列、平面阵列和任意阵列五种情况适合正在做阵列信号处理仿真、或者需要手写DOA/波束形成算法的朋友参考。1. 先搞清楚阵列响应矩阵到底在描述什么1.1 为什么几乎所有阵列算法都绕不开这个矩阵阵列信号处理的核心任务是从多个传感器接收到的信号中提取空间信息比如信号来自哪个方向、有几个信号、信号强度如何。无论是MUSIC、ESPRIT这类超分辨DOA算法还是MVDR、LCMV这类自适应波束形成算法第一步几乎都是构造阵列的响应矩阵也就是常说的阵列流型矩阵Array Manifold Matrix。它描述的是当单位幅度的平面波从某个方向入射时阵列各个阵元接收到的复信号之间幅度相同、相位差异是多少。这个相位差异完全由阵元位置和入射方向决定。换句话说阵列响应矩阵就是阵列几何结构与空间方向之间的一张“查找表”。有了它你才能把“方向”这个物理量换算成“阵元间的相位差”进而换算成协方差矩阵里的特征结构。我用一个生活化的类比如果把阵列比作一排话筒那响应矩阵就是提前记录好的“每个歌手站在台上不同位置时每个话筒收到的声音相位差”。没有这张表后面所有声源定位的算法都是空中楼阁。1.2 远场窄带假设是这一切的前提讨论阵列响应矩阵之前必须明确模型假设。绝大多数的阵列信号处理教材和工程实现都建立在两个前提上远场假设信号源离阵列足够远到达阵列的波前可以近似为平面波。这样各阵元接收信号的差异就只体现在波程差导致的相位差上而不是幅度衰减差上。窄带假设信号带宽远小于载波频率包络在跨越阵列孔径的时间内基本不变。这样信号可以用复包络表示阵元间的时延可以近似为相移。在这两个假设下对于包含 (N) 个阵元的阵列一个来自方向 (\theta)一维或 ((\theta, \varphi))二维的信号其阵列响应向量导向矢量可以写成[ \mathbf{a}(\theta) \left[ e^{j \cdot \frac{2\pi}{\lambda} \cdot \tau_1(\theta)},\ e^{j \cdot \frac{2\pi}{\lambda} \cdot \tau_2(\theta)},\ \cdots,\ e^{j \cdot \frac{2\pi}{\lambda} \cdot \tau_N(\theta)} \right]^T ]其中 (\tau_i(\theta)) 是第 (i) 个阵元相对于参考点的波程差(\lambda) 是载波波长。当有 (K) 个来自不同方向的信号时把每个方向的导向矢量按列拼在一起就得到阵列响应矩阵[ \mathbf{A} [\mathbf{a}(\theta_1),\ \mathbf{a}(\theta_2),\ \cdots,\ \mathbf{a}(\theta_K)] ]后面所有的MATLAB实现本质上都是在算这个矩阵。不同阵列的区别只是 (\tau_i(\theta)) 的表达式不同。2. 均匀线阵最基础也是最容易验证的案例2.1 均匀线阵导向矢量的推导与MATLAB实现均匀线阵ULA是最简单的阵列结构(N) 个全向阵元沿一条直线等间距排列间距为 (d)。假设信号从与阵列法线夹角为 (\theta) 的方向入射取第一个阵元为相位参考点则第 (n) 个阵元相对参考点的波程差为[ \tau_n (n-1) \cdot d \cdot \sin\theta ]注意这里 (\theta) 的定义不同教材可能不一样。有的定义成与端射方向的夹角有的定义成与法线的夹角这直接决定了后面是 (\sin\theta) 还是 (\cos\theta)。我习惯采用“与法线夹角”的定义这也是MATLAB中大多数工具箱默认的方式。对应的导向矢量为[ \mathbf{a}(\theta) \left[ 1,\ e^{-j \frac{2\pi}{\lambda} d \sin\theta},\ e^{-j \frac{2\pi}{\lambda} 2d \sin\theta},\ \cdots,\ e^{-j \frac{2\pi}{\lambda} (N-1)d \sin\theta} \right]^T ]MATLAB实现非常简洁而且我强烈建议用向量化写法不要用for循环function a ula_steering(theta_deg, d, N, lambda) % theta_deg: 入射角度单位度 % d: 阵元间距单位米 % N: 阵元个数 % lambda: 载波波长单位米 theta deg2rad(theta_deg); n (0:N-1).; a exp(-1j * 2 * pi * d / lambda * n * sin(theta)); end调用示例lambda 0.3; % 假设载波频率1GHz d lambda / 2; % 经典半波长间距 N 8; % 8元线阵 a_30deg ula_steering(30, d, N, lambda);这里有一个细节n * sin(theta)是一个 (N \times 1) 的列向量整个表达式用了矩阵乘法避免了循环。当需要同时计算多个方向时可以把输入改成角度向量利用MATLAB的隐式扩展直接得到矩阵效率更高。2.2 从导向矢量到响应矩阵的工程习惯实际仿真时我们很少只计算单一方向的导向矢量更多是构造一整张流型表。比如待扫描的角度范围是 (-90^\circ) 到 (90^\circ)步进 (0.1^\circ)那么响应矩阵就是 (N \times 1801) 的复数矩阵。我通常这样写theta_scan -90:0.1:90; A_scan zeros(N, length(theta_scan)); for idx 1:length(theta_scan) A_scan(:, idx) ula_steering(theta_scan(idx), d, N, lambda); end这里的循环只扫描了1801个点总共8个阵元耗时在毫秒级完全不用优化。但如果扫描角度很密、阵元数很大还是建议直接用矩阵运算一次生成。把函数改成支持向量输入function A ula_steering_multi(theta_vec_deg, d, N, lambda) theta_vec deg2rad(theta_vec_deg(:)).; % 1 x M n (0:N-1).; % N x 1 A exp(-1j * 2 * pi * d / lambda * n * sin(theta_vec)); % N x M end这样一次调用就能生成整个响应矩阵。顺便说一句生成响应矩阵后我习惯先画一下某个方向的波束图验证正确性。对于一个8元半波长间距线阵30度方向入射用 ( \mathbf{a}^H(\theta_{scan}) \cdot \mathbf{a}(30^\circ) ) 画出的幅度响应应该在30度处有明显峰值。如果峰值位置不对先检查角度定义是否统一再检查符号是否反了。2.3 均匀线阵的间距选择与相位参考点问题做线阵仿真时阵元间距 (d) 是一个必须认真对待的参数。最常见的设置是 (d \lambda/2)原因是空间采样定理当间距超过半波长时导向矢量会出现栅瓣导致多个角度对应的响应完全相同DOA估计就会出现模糊。我在初学阶段踩过一次坑把间距设成了 (d \lambda)结果MUSIC谱在 (-30^\circ) 和 (30^\circ) 同时出现峰值当时还以为是算法写错了排查了半天才发现是间距问题。另外一个容易被忽略的问题是相位参考点的选择。上面的公式默认取第一个阵元为参考点但有些文献取阵列中心为参考点。参考点变了导向矢量的整体相位会旋转但各阵元之间的相对相位差不变。对于MUSIC这类利用特征子空间正交性的算法整体相位旋转会被协方差矩阵的求逆或者特征分解过程吸收掉不影响最终结果。但如果你在做相位校准、或者把多个子阵拼接起来做联合估计就必须统一参考点否则相位对齐会出问题。3. 均匀圆阵、L型阵列与平面阵列二维角度估计的场景3.1 均匀圆阵的导向矢量实现与角度定义陷阱均匀圆阵UCA由 (M) 个全向阵元均匀分布在半径为 (R) 的圆周上。它相比线阵最大的优势是能够同时估计方位角和俯仰角而且方位角覆盖范围是完整的 (360^\circ)不存在线阵的“前视/后视模糊”问题。设第 (m) 个阵元所在的方位角为[ \varphi_m \frac{2\pi (m-1)}{M} ]当信号从方位角 (\varphi)、俯仰角 (\theta) 入射时这里的俯仰角我定义为与z轴正方向的夹角(\theta0) 表示从阵列正上方入射以圆心为相位参考点第 (m) 个阵元的波程差为[ \tau_m -R \cos(\varphi - \varphi_m) \sin\theta ]导向矢量为[ a_m(\varphi, \theta) e^{j \frac{2\pi}{\lambda} R \cos(\varphi - \varphi_m) \sin\theta} ]MATLAB实现如下function a uca_steering(phi_deg, theta_deg, R, M, lambda) phi deg2rad(phi_deg); theta deg2rad(theta_deg); phi_m 2 * pi * (0:M-1). / M; a exp(1j * 2 * pi * R / lambda * cos(phi - phi_m) * sin(theta)); end注意这里的符号和线阵版本不一样线阵我用的是-1j圆阵这里用的是1j。本质上取决于波程差的正负号约定只要在整个系统中保持一致即可。我见过不少同学把线阵的代码直接抄过来改成圆阵符号不一致导致波束指向完全错误。我的建议是写代码时在注释里明确标注“相位参考点为圆心”和“入射方向与z轴正向的夹角”避免自己下次也被绕晕。圆阵还有一个关键参数是半径 (R)。工程上常见的选择是 (R \lambda / (4 \sin(\pi/M)))这个值对应的是相邻阵元之间的弧长约为半波长。如果半径太大圆阵的空域采样会出现混叠太小则阵列孔径不够分辨率下降。3.2 L型阵列两个子阵的拼接与旋转不变性L型阵列由两条相互垂直的均匀线阵组成典型的配置是x轴方向有 (N_x) 个阵元y轴方向有 (N_y) 个阵元原点处共用同一个阵元。它的最大价值在于既保留了均匀线阵结构简单、处理方便的优点又天然具备二维角度估计能力而且非常适合ESPRIT算法——因为两个子阵之间存在旋转不变关系。假设阵元间距都是 (d)信号从方位角 (\varphi)、俯仰角 (\theta) 入射方位角定义为与x轴的夹角俯仰角定义为与z轴的夹角则x轴子阵的导向矢量为[ \mathbf{a}_x \left[ 1,\ e^{-j \frac{2\pi}{\lambda} d \sin\theta \cos\varphi},\ \cdots,\ e^{-j \frac{2\pi}{\lambda} (N_x-1)d \sin\theta \cos\varphi} \right]^T ]y轴子阵的导向矢量为[ \mathbf{a}_y \left[ 1,\ e^{-j \frac{2\pi}{\lambda} d \sin\theta \sin\varphi},\ \cdots,\ e^{-j \frac{2\pi}{\lambda} (N_y-1)d \sin\theta \sin\varphi} \right]^T ]整体阵列的响应向量一般有两种拼法一种是直接[a_x; a_y(2:end)]去掉原点处重复的阵元另一种是[a_x; a_y]保留原点重复。两种都能用但后续算法在构造协方差矩阵时如果保留了重复阵元协方差矩阵的噪声子空间维度会多一维对算法没有实质影响只是要注意维度匹配。MATLAB实现function a larray_steering(phi_deg, theta_deg, d, Nx, Ny, lambda) phi deg2rad(phi_deg); theta deg2rad(theta_deg); nx (0:Nx-1).; ny (1:Ny-1).; % 从1开始跳过原点处重复阵元 a_x exp(-1j * 2 * pi * d / lambda * nx * sin(theta) * cos(phi)); a_y exp(-1j * 2 * pi * d / lambda * ny * sin(theta) * sin(phi)); a [a_x; a_y]; end这里有一个实际工程中常遇到的问题L型阵列的两个子阵性能不一致因为x轴子阵只对 (\cos\varphi) 敏感y轴子阵只对 (\sin\varphi) 敏感。当信号从接近 (0^\circ) 方向入射时y轴子阵的导向矢量几乎全为1对俯仰角的估计能力会下降。这不是代码问题而是阵列结构本身的特性。如果项目对全角度范围内的估计一致性要求很高L型阵列未必是最优选我会建议考虑圆阵或者面阵。3.3 平面阵列用Kronecker积快速生成响应矩阵平面阵列URA是另一种常用的二维阵列阵元在x-y平面上按矩形网格排列x方向 (N_x) 个阵元、间距 (d_x)y方向 (N_y) 个阵元、间距 (d_y)。它和L型阵列最大的区别是平面阵是一个完整的二维孔径而不是两个一维子阵的简单拼接因此在二维波束形成时主瓣更窄、旁瓣更低。假设信号从方位角 (\varphi)、俯仰角 (\theta) 入射第 ((m, n)) 个阵元m对应x方向索引n对应y方向索引的波程差为[ \tau_{m,n} (m-1)d_x \sin\theta \cos\varphi (n-1)d_y \sin\theta \sin\varphi ]导向矢量可以写成两个一维导向矢量的Kronecker积[ \mathbf{a}_{URA}(\varphi, \theta) \mathbf{a}_y(\varphi, \theta) \otimes \mathbf{a}_x(\varphi, \theta) ]其中 (\mathbf{a}_x) 和 (\mathbf{a}_y) 分别是沿x轴和y轴的线阵导向矢量。这个性质非常重要它意味着平面阵的二维扫描可以拆成两个一维扫描计算量大幅降低。MATLAB实现function a ura_steering(phi_deg, theta_deg, dx, dy, Nx, Ny, lambda) phi deg2rad(phi_deg); theta deg2rad(theta_deg); nx (0:Nx-1).; ny (0:Ny-1).; a_x exp(-1j * 2 * pi * dx / lambda * nx * sin(theta) * cos(phi)); a_y exp(-1j * 2 * pi * dy / lambda * ny * sin(theta) * sin(phi)); a kron(a_y, a_x); % 注意kron的顺序 end关于kron的顺序我踩过坑。kron(a_y, a_x)出来的向量前 (N_x) 个元素对应y方向索引为0时的x方向各阵元这与我们习惯的“先x后y”索引顺序相反。如果你后续要按reshape成二维矩阵来画阵列激励分布必须先搞清楚这个顺序否则画出来的波束图方向是歪的。建议用一个小例子验证设 (N_x2, N_y3)手动算一下kron(a_y, a_x)的索引对应关系再对照你的后续处理需求调整顺序。还有一个细节当 (d_x d_y \lambda/2) 时矩阵运算中的sin(theta)*cos(phi)和sin(theta)*sin(phi)要确保角度单位一致。我在实际项目里为了避免单位混淆会在函数入口统一转成弧度所有内部计算只用弧度只在外部接口保留角度制。4. 任意阵列的通用实现一套代码通吃所有几何结构4.1 从位置矩阵到导向矢量的通用公式前几节讲的都是规则阵列每种阵列都推导了专门的解析表达式。但在实际项目中阵列结构往往不是标准的——可能是圆环的一部分、可能是随机布阵、可能是带有阵元位置误差的标称阵列。这时候最稳妥的做法是直接用阵元位置坐标计算导向矢量一套代码通吃所有阵列几何。通用公式并不复杂。设第 (i) 个阵元的空间坐标为 ((x_i, y_i, z_i))信号从方向 ((\theta, \varphi)) 入射入射方向的单位向量为[ \mathbf{u} [\sin\theta \cos\varphi,\ \sin\theta \sin\varphi,\ \cos\theta]^T ]以原点为相位参考点第 (i) 个阵元与参考点之间的波程差为[ \tau_i x_i \sin\theta \cos\varphi y_i \sin\theta \sin\varphi z_i \cos\theta ]导向矢量为[ a_i(\theta, \varphi) e^{j \frac{2\pi}{\lambda} \tau_i} ]注意符号约定。如果我用exp(j...)对应的DOA算法在构造协方差矩阵时也要用一致的约定。这个公式的价值在于任何阵元位置矩阵都可以直接代入计算不需要针对每种阵列单独推导。4.2 通用函数的MATLAB实现与验证方法function a array_steering(pos, phi_deg, theta_deg, lambda) % pos: N x 3 矩阵每行是一个阵元的 [x, y, z] 坐标 % phi_deg: 方位角度 % theta_deg: 俯仰角度 phi deg2rad(phi_deg); theta deg2rad(theta_deg); u [sin(theta) * cos(phi); sin(theta) * sin(phi); cos(theta)]; tau pos * u; % N x 1每个阵元的波程差 a exp(1j * 2 * pi / lambda * tau); end对于均匀线阵也可以用这个通用函数验证。比如8个阵元沿x轴排列间距半波长位置矩阵就是N 8; d lambda / 2; pos [(0:N-1). * d, zeros(N, 2)]; a_generic array_steering(pos, 0, 90 - 30, lambda);注意这里的角度换算。由于通用函数里的 (\theta) 定义是与z轴正向的夹角要模拟信号从x-y平面内、与x轴成30度方向入射需要让 (\theta 60^\circ)从z轴转到x-y平面然后 (\varphi 0^\circ)。这种角度变换是使用通用函数时最容易出错的地方我建议在写代码之前先画一个坐标系示意图把所有角度定义标注清楚。4.3 网格阵列与稀疏阵列的位置生成除了规则阵列任意阵列代码在稀疏阵列sparse array仿真中特别常用。比如我最近在做的一个项目需要在给定孔径内随机布置阵元然后评估不同布阵方案的角度估计精度。位置矩阵的生成很简单rng(42); N 16; aperture 4 * lambda; pos [rand(N, 1) * aperture, rand(N, 1) * aperture, zeros(N, 1)];生成之后先用刚才的通用函数算一个方向上的导向矢量再画一下各阵元的相位分布。如果相位分布和理论波程差对不上大概率是位置坐标的单位问题——记得所有坐标都要用米不要混用毫米或波长。对于网格阵列生成位置矩阵更直接[x, y] meshgrid(0:Nx-1, 0:Ny-1); pos [x(:) * dx, y(:) * dy, zeros(Nx*Ny, 1)];然后代入通用函数得到的结果和前面ura_steering用Kronecker积算出来的应该完全一致。如果对不上就可以用这个通用版来排查是哪个环节出了问题。这种“通用实现打底、专用函数做快速验证”的组合是我在实际项目中比较推荐的做法。5. 实操中的避坑经验与性能验证5.1 最容易翻车的5个细节对照表我在复现这本书的代码、以及自己写项目的过程里整理了阵列响应矩阵实现中几个高频踩坑点列成一张表供参考常见问题典型表现排查方向角度单位混用波束指向偏差极大甚至指向反方向统一在函数入口转弧度内部只用弧度指数符号反了波束指向变成镜像角度检查exp(1j...)还是exp(-1j...)与协方差矩阵构造保持一致间距超过半波长MUSIC谱出现两个相邻峰值检查 (d) 是否大于 (\lambda/2)圆阵俯仰角定义混乱仿真结果与理论波束图不一致明确俯仰角是与z轴夹角还是与x-y平面夹角Kronecker积顺序不对平面阵二维波束图方向偏转用 (N_x2, N_y3) 的小阵列手推验证这些问题的共同根源是定义不统一。我个人的经验是在项目初期建立一个“约定文档”把坐标系、角度定义、相位参考点、符号约定全部写清楚即使是自己一个人做项目也要写否则过两周回来看代码会非常痛苦。5.2 用协方差矩阵和MUSIC谱验证响应矩阵的正确性构造响应矩阵本身不难难的是确认它是对的。我的验证方法是用MUSIC算法做端到端测试如果MUSIC谱在真实信号方向出现峰值说明响应矩阵的构造是正确的如果峰值偏了、出现伪峰、或者干脆没有峰问题大概率出在响应矩阵上。完整的验证代码如下% 参数设置 lambda 0.3; d lambda / 2; N 8; K 2; % 两个信号 theta_true [-20, 35]; % 真实方向 SNR 20; % 信噪比dB snap 500; % 构造响应矩阵 A zeros(N, K); for idx 1:K A(:, idx) ula_steering(theta_true(idx), d, N, lambda); end % 生成接收数据 S (randn(K, snap) 1j * randn(K, snap)) / sqrt(2); X A * S; noise (randn(N, snap) 1j * randn(N, snap)) / sqrt(2) * 10^(-SNR/20); X X noise; % 协方差矩阵与特征分解 Rxx X * X / snap; [V, D] eig(Rxx); eigval diag(D); [~, idx_sort] sort(eigval); Vn V(:, idx_sort(1:N-K)); % 噪声子空间 % MUSIC谱扫描 theta_scan -90:0.1:90; P_music zeros(size(theta_scan)); for idx 1:length(theta_scan) a_theta ula_steering(theta_scan(idx), d, N, lambda); P_music(idx) 1 / (a_theta * (Vn * Vn) * a_theta); end P_music 10 * log10(P_music / max(P_music)); % 找峰值 [pks, locs] findpeaks(P_music, MinPeakHeight, -10); est_theta theta_scan(locs); disp(估计角度:); disp(est_theta);当响应矩阵正确时est_theta应该非常接近[-20, 35]。实测SNR在20dB、500快拍下误差通常在0.1度以内。如果估计角度偏差很大先检查响应矩阵的相位分布图再检查噪声子空间的维度是否正确。5.3 响应矩阵在分布式阵列与多子阵场景中的扩展思考最近很多同行开始讨论分布式阵列信号处理也就是把多个子阵分散布置在较大范围内通过协同处理获得更大的虚拟孔径。这时候响应矩阵的构造逻辑就有了新变化每个子阵都有自己的位置和姿态整体响应矩阵要考虑子阵间的相位中心差。我在测试这类方案时通用任意阵列实现几乎成了刚需——只需要把所有子阵的阵元坐标换算到同一个全局坐标系下然后直接调用一次array_steering就好不需要为每个子阵单独写公式。但这里有个前提子阵间的时间同步和相位校准必须做扎实否则响应矩阵算得再准实际数据也对不上。我的建议是仿真阶段就在接收数据中加入随机的子阵间相位偏移检验估计算法对校准误差的鲁棒性。这样能更早暴露工程实现中的隐患而不是到了外场试验才措手不及。 阵列流型矩阵MATLAB实现：线阵、圆阵、L型及任意阵列导向矢量代码详解
做了这么多年阵列信号处理每次带新人入门我最推荐的章节就是《阵列信号处理及MATLAB实现》里关于阵列响应矩阵的这部分。不管是做DOA估计还是波束形成最后都要落到阵列流型矩阵的构造上而不少同学恰恰是在这里开始出问题——线阵还算顺利一换成圆阵、L型或任意阵代码就开始乱。这篇文章把我自己在这个项目里反复验证过的实现思路、MATLAB代码和踩坑记录整理出来覆盖均匀线阵、均匀圆阵、L型阵列、平面阵列和任意阵列五种情况适合正在做阵列信号处理仿真、或者需要手写DOA/波束形成算法的朋友参考。1. 先搞清楚阵列响应矩阵到底在描述什么1.1 为什么几乎所有阵列算法都绕不开这个矩阵阵列信号处理的核心任务是从多个传感器接收到的信号中提取空间信息比如信号来自哪个方向、有几个信号、信号强度如何。无论是MUSIC、ESPRIT这类超分辨DOA算法还是MVDR、LCMV这类自适应波束形成算法第一步几乎都是构造阵列的响应矩阵也就是常说的阵列流型矩阵Array Manifold Matrix。它描述的是当单位幅度的平面波从某个方向入射时阵列各个阵元接收到的复信号之间幅度相同、相位差异是多少。这个相位差异完全由阵元位置和入射方向决定。换句话说阵列响应矩阵就是阵列几何结构与空间方向之间的一张“查找表”。有了它你才能把“方向”这个物理量换算成“阵元间的相位差”进而换算成协方差矩阵里的特征结构。我用一个生活化的类比如果把阵列比作一排话筒那响应矩阵就是提前记录好的“每个歌手站在台上不同位置时每个话筒收到的声音相位差”。没有这张表后面所有声源定位的算法都是空中楼阁。1.2 远场窄带假设是这一切的前提讨论阵列响应矩阵之前必须明确模型假设。绝大多数的阵列信号处理教材和工程实现都建立在两个前提上远场假设信号源离阵列足够远到达阵列的波前可以近似为平面波。这样各阵元接收信号的差异就只体现在波程差导致的相位差上而不是幅度衰减差上。窄带假设信号带宽远小于载波频率包络在跨越阵列孔径的时间内基本不变。这样信号可以用复包络表示阵元间的时延可以近似为相移。在这两个假设下对于包含 (N) 个阵元的阵列一个来自方向 (\theta)一维或 ((\theta, \varphi))二维的信号其阵列响应向量导向矢量可以写成[ \mathbf{a}(\theta) \left[ e^{j \cdot \frac{2\pi}{\lambda} \cdot \tau_1(\theta)},\ e^{j \cdot \frac{2\pi}{\lambda} \cdot \tau_2(\theta)},\ \cdots,\ e^{j \cdot \frac{2\pi}{\lambda} \cdot \tau_N(\theta)} \right]^T ]其中 (\tau_i(\theta)) 是第 (i) 个阵元相对于参考点的波程差(\lambda) 是载波波长。当有 (K) 个来自不同方向的信号时把每个方向的导向矢量按列拼在一起就得到阵列响应矩阵[ \mathbf{A} [\mathbf{a}(\theta_1),\ \mathbf{a}(\theta_2),\ \cdots,\ \mathbf{a}(\theta_K)] ]后面所有的MATLAB实现本质上都是在算这个矩阵。不同阵列的区别只是 (\tau_i(\theta)) 的表达式不同。2. 均匀线阵最基础也是最容易验证的案例2.1 均匀线阵导向矢量的推导与MATLAB实现均匀线阵ULA是最简单的阵列结构(N) 个全向阵元沿一条直线等间距排列间距为 (d)。假设信号从与阵列法线夹角为 (\theta) 的方向入射取第一个阵元为相位参考点则第 (n) 个阵元相对参考点的波程差为[ \tau_n (n-1) \cdot d \cdot \sin\theta ]注意这里 (\theta) 的定义不同教材可能不一样。有的定义成与端射方向的夹角有的定义成与法线的夹角这直接决定了后面是 (\sin\theta) 还是 (\cos\theta)。我习惯采用“与法线夹角”的定义这也是MATLAB中大多数工具箱默认的方式。对应的导向矢量为[ \mathbf{a}(\theta) \left[ 1,\ e^{-j \frac{2\pi}{\lambda} d \sin\theta},\ e^{-j \frac{2\pi}{\lambda} 2d \sin\theta},\ \cdots,\ e^{-j \frac{2\pi}{\lambda} (N-1)d \sin\theta} \right]^T ]MATLAB实现非常简洁而且我强烈建议用向量化写法不要用for循环function a ula_steering(theta_deg, d, N, lambda) % theta_deg: 入射角度单位度 % d: 阵元间距单位米 % N: 阵元个数 % lambda: 载波波长单位米 theta deg2rad(theta_deg); n (0:N-1).; a exp(-1j * 2 * pi * d / lambda * n * sin(theta)); end调用示例lambda 0.3; % 假设载波频率1GHz d lambda / 2; % 经典半波长间距 N 8; % 8元线阵 a_30deg ula_steering(30, d, N, lambda);这里有一个细节n * sin(theta)是一个 (N \times 1) 的列向量整个表达式用了矩阵乘法避免了循环。当需要同时计算多个方向时可以把输入改成角度向量利用MATLAB的隐式扩展直接得到矩阵效率更高。2.2 从导向矢量到响应矩阵的工程习惯实际仿真时我们很少只计算单一方向的导向矢量更多是构造一整张流型表。比如待扫描的角度范围是 (-90^\circ) 到 (90^\circ)步进 (0.1^\circ)那么响应矩阵就是 (N \times 1801) 的复数矩阵。我通常这样写theta_scan -90:0.1:90; A_scan zeros(N, length(theta_scan)); for idx 1:length(theta_scan) A_scan(:, idx) ula_steering(theta_scan(idx), d, N, lambda); end这里的循环只扫描了1801个点总共8个阵元耗时在毫秒级完全不用优化。但如果扫描角度很密、阵元数很大还是建议直接用矩阵运算一次生成。把函数改成支持向量输入function A ula_steering_multi(theta_vec_deg, d, N, lambda) theta_vec deg2rad(theta_vec_deg(:)).; % 1 x M n (0:N-1).; % N x 1 A exp(-1j * 2 * pi * d / lambda * n * sin(theta_vec)); % N x M end这样一次调用就能生成整个响应矩阵。顺便说一句生成响应矩阵后我习惯先画一下某个方向的波束图验证正确性。对于一个8元半波长间距线阵30度方向入射用 ( \mathbf{a}^H(\theta_{scan}) \cdot \mathbf{a}(30^\circ) ) 画出的幅度响应应该在30度处有明显峰值。如果峰值位置不对先检查角度定义是否统一再检查符号是否反了。2.3 均匀线阵的间距选择与相位参考点问题做线阵仿真时阵元间距 (d) 是一个必须认真对待的参数。最常见的设置是 (d \lambda/2)原因是空间采样定理当间距超过半波长时导向矢量会出现栅瓣导致多个角度对应的响应完全相同DOA估计就会出现模糊。我在初学阶段踩过一次坑把间距设成了 (d \lambda)结果MUSIC谱在 (-30^\circ) 和 (30^\circ) 同时出现峰值当时还以为是算法写错了排查了半天才发现是间距问题。另外一个容易被忽略的问题是相位参考点的选择。上面的公式默认取第一个阵元为参考点但有些文献取阵列中心为参考点。参考点变了导向矢量的整体相位会旋转但各阵元之间的相对相位差不变。对于MUSIC这类利用特征子空间正交性的算法整体相位旋转会被协方差矩阵的求逆或者特征分解过程吸收掉不影响最终结果。但如果你在做相位校准、或者把多个子阵拼接起来做联合估计就必须统一参考点否则相位对齐会出问题。3. 均匀圆阵、L型阵列与平面阵列二维角度估计的场景3.1 均匀圆阵的导向矢量实现与角度定义陷阱均匀圆阵UCA由 (M) 个全向阵元均匀分布在半径为 (R) 的圆周上。它相比线阵最大的优势是能够同时估计方位角和俯仰角而且方位角覆盖范围是完整的 (360^\circ)不存在线阵的“前视/后视模糊”问题。设第 (m) 个阵元所在的方位角为[ \varphi_m \frac{2\pi (m-1)}{M} ]当信号从方位角 (\varphi)、俯仰角 (\theta) 入射时这里的俯仰角我定义为与z轴正方向的夹角(\theta0) 表示从阵列正上方入射以圆心为相位参考点第 (m) 个阵元的波程差为[ \tau_m -R \cos(\varphi - \varphi_m) \sin\theta ]导向矢量为[ a_m(\varphi, \theta) e^{j \frac{2\pi}{\lambda} R \cos(\varphi - \varphi_m) \sin\theta} ]MATLAB实现如下function a uca_steering(phi_deg, theta_deg, R, M, lambda) phi deg2rad(phi_deg); theta deg2rad(theta_deg); phi_m 2 * pi * (0:M-1). / M; a exp(1j * 2 * pi * R / lambda * cos(phi - phi_m) * sin(theta)); end注意这里的符号和线阵版本不一样线阵我用的是-1j圆阵这里用的是1j。本质上取决于波程差的正负号约定只要在整个系统中保持一致即可。我见过不少同学把线阵的代码直接抄过来改成圆阵符号不一致导致波束指向完全错误。我的建议是写代码时在注释里明确标注“相位参考点为圆心”和“入射方向与z轴正向的夹角”避免自己下次也被绕晕。圆阵还有一个关键参数是半径 (R)。工程上常见的选择是 (R \lambda / (4 \sin(\pi/M)))这个值对应的是相邻阵元之间的弧长约为半波长。如果半径太大圆阵的空域采样会出现混叠太小则阵列孔径不够分辨率下降。3.2 L型阵列两个子阵的拼接与旋转不变性L型阵列由两条相互垂直的均匀线阵组成典型的配置是x轴方向有 (N_x) 个阵元y轴方向有 (N_y) 个阵元原点处共用同一个阵元。它的最大价值在于既保留了均匀线阵结构简单、处理方便的优点又天然具备二维角度估计能力而且非常适合ESPRIT算法——因为两个子阵之间存在旋转不变关系。假设阵元间距都是 (d)信号从方位角 (\varphi)、俯仰角 (\theta) 入射方位角定义为与x轴的夹角俯仰角定义为与z轴的夹角则x轴子阵的导向矢量为[ \mathbf{a}_x \left[ 1,\ e^{-j \frac{2\pi}{\lambda} d \sin\theta \cos\varphi},\ \cdots,\ e^{-j \frac{2\pi}{\lambda} (N_x-1)d \sin\theta \cos\varphi} \right]^T ]y轴子阵的导向矢量为[ \mathbf{a}_y \left[ 1,\ e^{-j \frac{2\pi}{\lambda} d \sin\theta \sin\varphi},\ \cdots,\ e^{-j \frac{2\pi}{\lambda} (N_y-1)d \sin\theta \sin\varphi} \right]^T ]整体阵列的响应向量一般有两种拼法一种是直接[a_x; a_y(2:end)]去掉原点处重复的阵元另一种是[a_x; a_y]保留原点重复。两种都能用但后续算法在构造协方差矩阵时如果保留了重复阵元协方差矩阵的噪声子空间维度会多一维对算法没有实质影响只是要注意维度匹配。MATLAB实现function a larray_steering(phi_deg, theta_deg, d, Nx, Ny, lambda) phi deg2rad(phi_deg); theta deg2rad(theta_deg); nx (0:Nx-1).; ny (1:Ny-1).; % 从1开始跳过原点处重复阵元 a_x exp(-1j * 2 * pi * d / lambda * nx * sin(theta) * cos(phi)); a_y exp(-1j * 2 * pi * d / lambda * ny * sin(theta) * sin(phi)); a [a_x; a_y]; end这里有一个实际工程中常遇到的问题L型阵列的两个子阵性能不一致因为x轴子阵只对 (\cos\varphi) 敏感y轴子阵只对 (\sin\varphi) 敏感。当信号从接近 (0^\circ) 方向入射时y轴子阵的导向矢量几乎全为1对俯仰角的估计能力会下降。这不是代码问题而是阵列结构本身的特性。如果项目对全角度范围内的估计一致性要求很高L型阵列未必是最优选我会建议考虑圆阵或者面阵。3.3 平面阵列用Kronecker积快速生成响应矩阵平面阵列URA是另一种常用的二维阵列阵元在x-y平面上按矩形网格排列x方向 (N_x) 个阵元、间距 (d_x)y方向 (N_y) 个阵元、间距 (d_y)。它和L型阵列最大的区别是平面阵是一个完整的二维孔径而不是两个一维子阵的简单拼接因此在二维波束形成时主瓣更窄、旁瓣更低。假设信号从方位角 (\varphi)、俯仰角 (\theta) 入射第 ((m, n)) 个阵元m对应x方向索引n对应y方向索引的波程差为[ \tau_{m,n} (m-1)d_x \sin\theta \cos\varphi (n-1)d_y \sin\theta \sin\varphi ]导向矢量可以写成两个一维导向矢量的Kronecker积[ \mathbf{a}_{URA}(\varphi, \theta) \mathbf{a}_y(\varphi, \theta) \otimes \mathbf{a}_x(\varphi, \theta) ]其中 (\mathbf{a}_x) 和 (\mathbf{a}_y) 分别是沿x轴和y轴的线阵导向矢量。这个性质非常重要它意味着平面阵的二维扫描可以拆成两个一维扫描计算量大幅降低。MATLAB实现function a ura_steering(phi_deg, theta_deg, dx, dy, Nx, Ny, lambda) phi deg2rad(phi_deg); theta deg2rad(theta_deg); nx (0:Nx-1).; ny (0:Ny-1).; a_x exp(-1j * 2 * pi * dx / lambda * nx * sin(theta) * cos(phi)); a_y exp(-1j * 2 * pi * dy / lambda * ny * sin(theta) * sin(phi)); a kron(a_y, a_x); % 注意kron的顺序 end关于kron的顺序我踩过坑。kron(a_y, a_x)出来的向量前 (N_x) 个元素对应y方向索引为0时的x方向各阵元这与我们习惯的“先x后y”索引顺序相反。如果你后续要按reshape成二维矩阵来画阵列激励分布必须先搞清楚这个顺序否则画出来的波束图方向是歪的。建议用一个小例子验证设 (N_x2, N_y3)手动算一下kron(a_y, a_x)的索引对应关系再对照你的后续处理需求调整顺序。还有一个细节当 (d_x d_y \lambda/2) 时矩阵运算中的sin(theta)*cos(phi)和sin(theta)*sin(phi)要确保角度单位一致。我在实际项目里为了避免单位混淆会在函数入口统一转成弧度所有内部计算只用弧度只在外部接口保留角度制。4. 任意阵列的通用实现一套代码通吃所有几何结构4.1 从位置矩阵到导向矢量的通用公式前几节讲的都是规则阵列每种阵列都推导了专门的解析表达式。但在实际项目中阵列结构往往不是标准的——可能是圆环的一部分、可能是随机布阵、可能是带有阵元位置误差的标称阵列。这时候最稳妥的做法是直接用阵元位置坐标计算导向矢量一套代码通吃所有阵列几何。通用公式并不复杂。设第 (i) 个阵元的空间坐标为 ((x_i, y_i, z_i))信号从方向 ((\theta, \varphi)) 入射入射方向的单位向量为[ \mathbf{u} [\sin\theta \cos\varphi,\ \sin\theta \sin\varphi,\ \cos\theta]^T ]以原点为相位参考点第 (i) 个阵元与参考点之间的波程差为[ \tau_i x_i \sin\theta \cos\varphi y_i \sin\theta \sin\varphi z_i \cos\theta ]导向矢量为[ a_i(\theta, \varphi) e^{j \frac{2\pi}{\lambda} \tau_i} ]注意符号约定。如果我用exp(j...)对应的DOA算法在构造协方差矩阵时也要用一致的约定。这个公式的价值在于任何阵元位置矩阵都可以直接代入计算不需要针对每种阵列单独推导。4.2 通用函数的MATLAB实现与验证方法function a array_steering(pos, phi_deg, theta_deg, lambda) % pos: N x 3 矩阵每行是一个阵元的 [x, y, z] 坐标 % phi_deg: 方位角度 % theta_deg: 俯仰角度 phi deg2rad(phi_deg); theta deg2rad(theta_deg); u [sin(theta) * cos(phi); sin(theta) * sin(phi); cos(theta)]; tau pos * u; % N x 1每个阵元的波程差 a exp(1j * 2 * pi / lambda * tau); end对于均匀线阵也可以用这个通用函数验证。比如8个阵元沿x轴排列间距半波长位置矩阵就是N 8; d lambda / 2; pos [(0:N-1). * d, zeros(N, 2)]; a_generic array_steering(pos, 0, 90 - 30, lambda);注意这里的角度换算。由于通用函数里的 (\theta) 定义是与z轴正向的夹角要模拟信号从x-y平面内、与x轴成30度方向入射需要让 (\theta 60^\circ)从z轴转到x-y平面然后 (\varphi 0^\circ)。这种角度变换是使用通用函数时最容易出错的地方我建议在写代码之前先画一个坐标系示意图把所有角度定义标注清楚。4.3 网格阵列与稀疏阵列的位置生成除了规则阵列任意阵列代码在稀疏阵列sparse array仿真中特别常用。比如我最近在做的一个项目需要在给定孔径内随机布置阵元然后评估不同布阵方案的角度估计精度。位置矩阵的生成很简单rng(42); N 16; aperture 4 * lambda; pos [rand(N, 1) * aperture, rand(N, 1) * aperture, zeros(N, 1)];生成之后先用刚才的通用函数算一个方向上的导向矢量再画一下各阵元的相位分布。如果相位分布和理论波程差对不上大概率是位置坐标的单位问题——记得所有坐标都要用米不要混用毫米或波长。对于网格阵列生成位置矩阵更直接[x, y] meshgrid(0:Nx-1, 0:Ny-1); pos [x(:) * dx, y(:) * dy, zeros(Nx*Ny, 1)];然后代入通用函数得到的结果和前面ura_steering用Kronecker积算出来的应该完全一致。如果对不上就可以用这个通用版来排查是哪个环节出了问题。这种“通用实现打底、专用函数做快速验证”的组合是我在实际项目中比较推荐的做法。5. 实操中的避坑经验与性能验证5.1 最容易翻车的5个细节对照表我在复现这本书的代码、以及自己写项目的过程里整理了阵列响应矩阵实现中几个高频踩坑点列成一张表供参考常见问题典型表现排查方向角度单位混用波束指向偏差极大甚至指向反方向统一在函数入口转弧度内部只用弧度指数符号反了波束指向变成镜像角度检查exp(1j...)还是exp(-1j...)与协方差矩阵构造保持一致间距超过半波长MUSIC谱出现两个相邻峰值检查 (d) 是否大于 (\lambda/2)圆阵俯仰角定义混乱仿真结果与理论波束图不一致明确俯仰角是与z轴夹角还是与x-y平面夹角Kronecker积顺序不对平面阵二维波束图方向偏转用 (N_x2, N_y3) 的小阵列手推验证这些问题的共同根源是定义不统一。我个人的经验是在项目初期建立一个“约定文档”把坐标系、角度定义、相位参考点、符号约定全部写清楚即使是自己一个人做项目也要写否则过两周回来看代码会非常痛苦。5.2 用协方差矩阵和MUSIC谱验证响应矩阵的正确性构造响应矩阵本身不难难的是确认它是对的。我的验证方法是用MUSIC算法做端到端测试如果MUSIC谱在真实信号方向出现峰值说明响应矩阵的构造是正确的如果峰值偏了、出现伪峰、或者干脆没有峰问题大概率出在响应矩阵上。完整的验证代码如下% 参数设置 lambda 0.3; d lambda / 2; N 8; K 2; % 两个信号 theta_true [-20, 35]; % 真实方向 SNR 20; % 信噪比dB snap 500; % 构造响应矩阵 A zeros(N, K); for idx 1:K A(:, idx) ula_steering(theta_true(idx), d, N, lambda); end % 生成接收数据 S (randn(K, snap) 1j * randn(K, snap)) / sqrt(2); X A * S; noise (randn(N, snap) 1j * randn(N, snap)) / sqrt(2) * 10^(-SNR/20); X X noise; % 协方差矩阵与特征分解 Rxx X * X / snap; [V, D] eig(Rxx); eigval diag(D); [~, idx_sort] sort(eigval); Vn V(:, idx_sort(1:N-K)); % 噪声子空间 % MUSIC谱扫描 theta_scan -90:0.1:90; P_music zeros(size(theta_scan)); for idx 1:length(theta_scan) a_theta ula_steering(theta_scan(idx), d, N, lambda); P_music(idx) 1 / (a_theta * (Vn * Vn) * a_theta); end P_music 10 * log10(P_music / max(P_music)); % 找峰值 [pks, locs] findpeaks(P_music, MinPeakHeight, -10); est_theta theta_scan(locs); disp(估计角度:); disp(est_theta);当响应矩阵正确时est_theta应该非常接近[-20, 35]。实测SNR在20dB、500快拍下误差通常在0.1度以内。如果估计角度偏差很大先检查响应矩阵的相位分布图再检查噪声子空间的维度是否正确。5.3 响应矩阵在分布式阵列与多子阵场景中的扩展思考最近很多同行开始讨论分布式阵列信号处理也就是把多个子阵分散布置在较大范围内通过协同处理获得更大的虚拟孔径。这时候响应矩阵的构造逻辑就有了新变化每个子阵都有自己的位置和姿态整体响应矩阵要考虑子阵间的相位中心差。我在测试这类方案时通用任意阵列实现几乎成了刚需——只需要把所有子阵的阵元坐标换算到同一个全局坐标系下然后直接调用一次array_steering就好不需要为每个子阵单独写公式。但这里有个前提子阵间的时间同步和相位校准必须做扎实否则响应矩阵算得再准实际数据也对不上。我的建议是仿真阶段就在接收数据中加入随机的子阵间相位偏移检验估计算法对校准误差的鲁棒性。这样能更早暴露工程实现中的隐患而不是到了外场试验才措手不及。