最近在啃《阵列信号处理及MATLAB实现》啃到阵列响应矩阵这一章的时候有种豁然开朗的感觉。以前做波束成形和DOA估计公式来回套总觉得差一口气后来发现差的就是对“阵列响应矩阵”本身的理解。你可以把它理解成阵列的“指纹”——来波方向一变这个矩阵就跟着变而MVDR、MUSIC、ESPRIT这些经典算法干的事情本质上都是从接收数据里反推这个“指纹”对应的方向。这篇文章主要聊均匀线阵、均匀圆阵、L型阵列、平面阵列和任意阵列这五种构型的阵列响应矩阵怎么建模以及用MATLAB实现的完整代码和调试心得。产品经理会告诉你“这个功能很强大”工程师会告诉你“跑通了你才知道有多爽”。适合刚上手阵列信号处理的研究生、转算法岗的工程师以及所有想把公式变成可运行代码的同学。1. 阵列响应矩阵到底是什么先搞懂这个“指纹”模型先说两个前提假设因为后面所有数学推导都建立在它们之上。第一个是窄带假设。信号带宽远小于载频这样不同阵元之间的时延可以近似等效成相移。如果信号是宽带信号比如语音信号单纯用相移描述就不够了得考虑频域处理或聚焦变换。第二个是远场假设。目标距离足够远入射波近似为平面波方向只由角度决定和距离无关。工程上只要阵元孔径远小于目标距离这个假设基本成立。比如麦克风阵列收音距离一两米以上就可以当作远场处理。在这两个假设下设阵列有M个阵元来波方向为θ以第一个阵元或阵列中心为参考点第m个阵元相对参考点的波程差是Δ_m那么对应的相移就是a_m(θ) exp(-j * 2π * Δ_m / λ)把M个阵元的响应排成列向量a(θ) [a_1(θ), a_2(θ), ..., a_M(θ)]ᵀ这就是导向矢量有的书叫方向向量或steering vector。当空间中有K个来波方向θ_1到θ_K时把这些列向量并排拼接就得到一个M×K的矩阵A [a(θ_1), a(θ_2), ..., a(θ_K)]这个A就是阵列响应矩阵也叫阵列流型矩阵。接收数据的模型可以写成X A * S N其中X是M×L的接收数据矩阵L为快拍数S是K×L的信号矩阵N是M×L的噪声矩阵。这个式子基本贯穿了所有阵列信号处理算法。所以求阵列响应矩阵并不是最终目的它是你搭建仿真环境、验证算法的基础设施。实际写代码的时候要注意一点不同教材对角度定义不一样有的用与阵列法线的夹角有的用与阵列轴向的夹角有的用方位角加俯仰角。我建议在一开始就固定一套角度定义并且注释写清楚否则后面换阵列构型时很容易乱。2. 五种阵列构型怎么建模从ULA到任意阵列2.1 均匀线阵ULA一维角度的经典入门均匀线阵是最简单也最常用的阵列构型M个阵元等间距排在一条直线上间距记为d来波方向θ通常定义为与阵列法线的夹角。以第1个阵元为参考第m个阵元的位置是(m - 1) * d波程差为(m - 1) * d * sin(θ)所以导向矢量是a_ULA(θ) exp(-j * 2π * d * sin(θ) / λ * [0, 1, ..., M - 1]ᵀ)MATLAB实现非常简洁function a ula_steering(M, d, theta_deg, lambda) % 均匀线阵导向矢量 % 输入: % M : 阵元数 % d : 阵元间距与lambda同单位 % theta_deg: 来波方向度可以为向量 % lambda : 波长 % 输出: % a : M×K 导向矢量矩阵 theta deg2rad(theta_deg(:).); a exp(-1j * 2 * pi * d / lambda * (0:M-1). * sin(theta)); end这里最关键的一个参数就是阵元间距d。常规工程中取d λ/2也就是半波长。为什么因为当d超过λ/2时sin(θ)的取值范围[-1,1]内可能出现多个角度对应相同的相位差造成方向估计模糊也就是栅瓣。取半波长时空间采样满足“空间Nyquist条件”相位差与方向是一一对应的。实际项目里我也会遇到必须用d λ/2的情况比如两个阵元之间要放结构件物理空间不够。这时单纯用均匀线阵就会出问题一般会考虑稀疏阵列设计或者用解模糊算法但这都是后话了。2.2 均匀圆阵UCA360度全方位覆盖的几何直觉均匀圆阵是把M个阵元均匀分布在半径为R的圆周上第m个阵元的角度位置为φ_m 2π * (m - 1) / M。以圆心为参考点来波方向θ通常取方位角从x轴正方向逆时针在圆上的投影关系会相对绕一些。第m个阵元与参考点之间的波程差是R * cos(θ - φ_m)所以导向矢量是a_UCA(θ) exp(j * 2π * R / λ * cos(θ - φ_m))注意这里我用了正号有的资料会写成负号。符号问题不影响理论推导但一定要和接收数据模型保持一致否则画出来的方向图会左右颠倒。这也是做圆阵时最容易踩的坑之一。function a uca_steering(M, R, theta_deg, lambda) % 均匀圆阵导向矢量方位角 % 输入: % M : 阵元数 % R : 圆阵半径 % theta_deg: 方位角度可以为向量 % lambda : 波长 % 输出: % a : M×K 导向矢量矩阵 theta deg2rad(theta_deg(:).); phi_m 2 * pi * (0:M-1). / M; a exp(1j * 2 * pi * R / lambda * cos(theta - phi_m)); end圆阵和线阵最大的区别在于线阵的导向矢量是一个标准的复指数序列具有范德蒙德结构很多算法可以直接利用这种结构圆阵则没有这种结构因为各阵元的投影不是均匀的。圆阵的优势是水平面内360度全方位覆盖各方向性能比较均匀不存在线阵那样的端射模糊问题。但它在DOA估计里不能直接套用MUSIC因为MUSIC假设阵列流型是满秩的、各导向矢量线性独立圆阵在某些角度和半径组合下流型矩阵的秩会退化。常用的做法是进行相位模式激励Phase Mode Excitation把圆阵变换到模式空间再套用常规算法。这个部分书里写得比较简略我实战中调试了很久才把模式空间的导向矢量拼对。2.3 L型阵列一个简单的二维扩展思路L型阵列实际上就是把两个均匀线阵垂直放置一个沿x轴一个沿y轴原点处共用同一个阵元。假设每个臂有M个阵元含原点那么总阵元数为2M - 1。来波方向用方位角θ和俯仰角ϕ共同描述。入射方向单位向量可以写成u [cos(θ) * cos(ϕ), sin(θ) * cos(ϕ), sin(ϕ)]ᵀx轴臂上第m个阵元的位置是((m - 1) * d, 0, 0)它和入射方向的内积就是(m - 1) * d * cos(θ) * cos(ϕ)所以x臂的导向矢量为a_x exp(-j * 2π * d / λ * cos(θ) * cos(ϕ) * [0, 1, ..., M - 1]ᵀ)y臂类似只是投影从cos(θ) * cos(ϕ)变成了sin(θ) * cos(ϕ)a_y exp(-j * 2π * d / λ * sin(θ) * cos(ϕ) * [0, 1, ..., M - 1]ᵀ)function A Larray_steering(M, d, theta_deg, phi_deg, lambda) % L型阵列导向矢量矩阵 % 输入: % M : 每个臂的阵元数含原点共用阵元 % d : 阵元间距 % theta_deg: 方位角度可以为向量 % phi_deg : 俯仰角度可以为向量 % lambda : 波长 % 输出: % A : (2*M-1)×K 导向矢量矩阵 theta deg2rad(theta_deg(:).); phi deg2rad(phi_deg(:).); temp_x exp(-1j * 2 * pi * d / lambda * (0:M-1). * (cos(theta) .* cos(phi))); temp_y exp(-1j * 2 * pi * d / lambda * (0:M-1). * (sin(theta) .* cos(phi))); % 注意原点阵元重复拼接时需要去掉一个 A [temp_x; temp_y(2:end, :)]; endL型阵列的核心优势是结构简单、孔径利用率高两个臂就能额外获得俯仰角信息硬件成本只增加了一个维度的线性阵列。但我不建议在L型阵列上直接做二维MUSIC然后分别估计方位角和俯仰角因为两个臂独立估计出来的角度需要配对配对错误在目标个数较多时非常头痛。实际工程里常用“联合估计算法”或者“基于旋转不变性的方法ESPRIT类”这样角度对是自动配好的。2.4 平面阵列UPA二维角度估计的标准答案平面阵列Uniform Planar Array, UPA是阵列信号处理里最常用的二维阵列构型。阵元在x-y平面上按矩形网格排列x方向有N_x个阵元间距dxy方向有N_y个阵元间距dy总共M N_x * N_y个阵元。来波方向同样用方位角θ和俯仰角ϕ表示。入射方向单位向量u [cos(θ) * sin(ϕ), sin(θ) * sin(ϕ), cos(ϕ)]ᵀ注意我这里用的俯仰角定义是“从z轴正方向往下偏转的角度”有的书定义是从x-y平面往上仰角两种定义会差一个余角。我在代码里用sin(ϕ)还是cos(ϕ)完全取决于这一定义。所以每次仿真前先把角度定义和坐标轴关系写清楚能省很多debug时间。阵元(mx, my)的位置是((mx - 1) * dx, (my - 1) * dy)导向矢量写成a(mx, my) exp(-j * 2π / λ * ((mx - 1) * dx * cos(θ) * sin(ϕ) (my - 1) * dy * sin(θ) * sin(ϕ)))MATLAB里可以先把x方向和y方向的一维导向矢量算出来再用Kronecker积合成二维function A upa_steering(Nx, Ny, dx, dy, theta_deg, phi_deg, lambda) % 均匀平面阵导向矢量矩阵 % 输入: % Nx, Ny : x和y方向的阵元数 % dx, dy : x和y方向的阵元间距 % theta_deg: 方位角度可以为向量 % phi_deg : 俯仰角度可以为向量 % lambda : 波长 % 输出: % A : (Nx*Ny)×K 导向矢量矩阵 theta deg2rad(theta_deg(:).); phi deg2rad(phi_deg(:).); ax exp(-1j * 2 * pi * dx / lambda * (0:Nx-1). * (cos(theta) .* sin(phi))); ay exp(-1j * 2 * pi * dy / lambda * (0:Ny-1). * (sin(theta) .* sin(phi))); % 张量积生成二维导向矢量矩阵 A kron(ax, ay); end用kron生成的时候务必搞清楚向量化顺序。我习惯把x方向作为内层、y方向作为外层这样生成的A每一行对应一个物理阵元位置可视化的时候方便reshape成Nx×Ny的二维数组检查方向图。平面阵列的优势在于可以获得完整的方向余弦信息因而可以用2D-MUSIC等算法直接在二维角度域上搜索。计算量会大不少因为扫描网格从一维变成二维了。如果[θ, ϕ]网格各取180个点和90个点那么扫描就是16200个方向每个方向都要计算一次谱值。工程上通常先用ESPRIT或波束形成做粗估再在局部区域做细网格搜索这样能把计算时间降低两个数量级。2.5 任意阵列一条代码通吃所有构型任意阵列是把问题抽象到最一般的层面已知每个阵元的空间坐标阵元可以放在三维空间中任何位置。来波方向用方向余弦向量u表示则第m个阵元的导向矢量为a_m exp(-j * 2π / λ * (x_m * u_x y_m * u_y z_m * u_z))对应MATLAB实现如下function A arbitrary_array_steering(pos, theta_deg, phi_deg, lambda) % 任意阵列导向矢量矩阵 % 输入: % pos : M×3 阵元坐标矩阵每一行为[x,y,z] % theta_deg : 方位角度可以为向量 % phi_deg : 俯仰角度可以为向量 % lambda : 波长 % 输出: % A : M×K 导向矢量矩阵 theta deg2rad(theta_deg(:).); phi deg2rad(phi_deg(:).); u [cos(theta) .* sin(phi); sin(theta) .* sin(phi); cos(phi)]; A exp(-1j * 2 * pi / lambda * (pos * u)); end这行代码里pos * u是矩阵乘法结果维度是M×K每个元素表示对应阵元和对应方向之间的空间相位差。任意阵列模型可以覆盖前面四种情况把pos设置成线阵坐标就是ULA设置成圆环坐标就是UCA设置成矩形网格就是UPA。实际项目里我一般会优先写这个通用函数再单独写各个构型的初始化函数用来生成pos坐标。这样做的好处是后续换阵列构型时只需要改坐标生成逻辑导向矢量计算、协方差矩阵、谱搜索这些代码完全不用动。3. 阵列响应矩阵的关键性质方向图、栅瓣和秩3.1 方向图看懂主瓣、栅瓣和旁瓣方向图描述的是“当来波方向是θ_0时如果我把波束指向θ阵列输出的幅度响应是多少”。数学上就是导向矢量在θ_0方向和θ方向的内积模值F(θ) |a(θ)ᵀ * a(θ_0)| / M当然这里没有加窗加窗后的方向图会改变主瓣宽度和旁瓣电平。用均匀线阵举例假设阵元间距d λ/2M 8来波方向为0度方向图会在0度出现一个主瓣。如果阵元间距增大到d λ你会发现在±30度和±90度附近会出现与主瓣等高的“栅瓣”。这在DOA估计里非常致命因为真实方向对应的谱峰和栅瓣对应的谱峰高度相同算法会分不清到底哪个才是真实来波。解释栅瓣的物理来源其实很简单。对阵元间距为d的均匀线阵两个方向θ_1和θ_2产生的导向矢量完全相同当且仅当满足sin(θ_1) - sin(θ_2) n * λ / d, n为非零整数如果d λ/2右边最小非零值是2而左边取值范围只有[-2, 2]所以等号只有在n 0时成立也就是θ_1 θ_2不会出现栅瓣。一旦d λ/2n ±1甚至更大的整数都可能落在取值范围内栅瓣就出现了。这里有个很实际的经验写仿真代码时先画一下方向图确认主瓣位置正确、没有栅瓣干扰再往下做算法验证。否则你后面调MUSIC谱峰找了好久都找不到错在哪里大概率就是阵列参数设置出了问题。3.2 秩与分辨力为什么阵列流型必须满列秩MUSIC算法能分辨K个信号的前提之一是阵列响应矩阵A满列秩即rank(A) K。这样接收数据协方差矩阵的信号子空间才是K维噪声子空间才是M-K维。对于均匀线阵不同来波方向对应的导向矢量是线性独立的理论上只要K ≤ MA就是满列秩。但数值上要注意当两个角度非常接近时对应的导向矢量也会非常接近协方差矩阵的条件数会变大MUSIC在这种情况下的分辨力受限于快拍数和信噪比。实际处理中我会用奇异值分解SVD看一下A的奇异值分布如果前K个奇异值和后面的差距不够大说明角度间隔太近或者阵列构型本身病态。均匀圆阵有一个比较隐蔽的问题当半径较小、阵元数较多时某些角度组合会导致阵列流型矩阵接近奇异。我自己试过R λ/8的圆阵在信号个数超过3个时MUSIC谱就乱七八糟后来换成R λ/2才好转。所以在用圆阵做DOA估计前最好先画一下“所有来波方向组合下A的最小奇异值”分布确认阵列构型不病态再放心跑算法。4. 实操把阵列响应矩阵用起来跑一个完整的MUSIC定位仿真这一节我们走一遍完整流程从生成阵列坐标和导向矢量到构造接收数据最后用MUSIC算法反解来波方向。这段代码你可以直接拷下去改参数用。4.1 数据生成从阵列响应矩阵到接收数据我以均匀线阵为例设M 8d λ/2两个来波方向分别是-20度和10度快拍数L 500信噪比约10 dB。clear; clc; close all; M 8; % 阵元数 d 0.5; % 阵元间距以波长为单位0.5表示半波长 lambda 1; % 波长归一化为1 theta_true [-20, 10]; % 真实来波方向度 K length(theta_true); % 信号源个数 L 500; % 快拍数 SNR 10; % 信噪比(dB) % 生成阵列响应矩阵 A exp(-1j * 2 * pi * d * (0:M-1). * sin(deg2rad(theta_true)) / lambda); % 生成信号复高斯随机信号相互独立 S (randn(K, L) 1j * randn(K, L)) / sqrt(2); % 生成噪声并按照SNR叠加 sigma_n sqrt(10^(-SNR/10)); N sigma_n * (randn(M, L) 1j * randn(M, L)) / sqrt(2); % 接收数据 X A * S N;这段代码里有两个细节值得注意。第一信号协方差矩阵不需要是满秩的两个信号如果完全相干比如多径环境MUSIC会失效我后面会再提。第二复噪声功率归一化后的系数是1/sqrt(2)是为了让噪声实部和虚部总功率为1这样SNR的定义比较直观。4.2 协方差矩阵估计一个容易被低估的细节用接收数据估计协方差矩阵时常用的是样本协方差Rxx X * X / L;严格来说当快拍数L有限时样本协方差和真实协方差有偏差偏差会导致MUSIC谱峰位置偏移。工程上L建议至少大于2M最好比M大5到10倍。如果快拍数实在有限可以做对角线加载在Rxx上加一个很小的正则项Rxx Rxx 1e-6 * eye(M);我见过不少同学直接用原始数据算MUSIC结果峰值乱跳加上对角线加载后瞬间稳定。这个技巧在低快拍、低信噪比场景下非常有效。4.3 谱搜索与DOA估计% 特征分解 [E, D] eig(Rxx); [~, idx] sort(diag(D), descend); E E(:, idx); En E(:, K1:end); % 噪声子空间 % 扫描角度范围 theta_scan -90:0.1:90; M_scan length(theta_scan); P_music zeros(1, M_scan); % 生成扫描导向矢量矩阵 A_scan exp(-1j * 2 * pi * d * (0:M-1). * sin(deg2rad(theta_scan)) / lambda); % MUSIC谱 for ii 1:M_scan a_theta A_scan(:, ii); P_music(ii) 1 / (a_theta * (En * En) * a_theta); end % 归一化并画图 P_music 10 * log10(P_music / max(P_music)); figure; plot(theta_scan, P_music, LineWidth, 1.2); xlabel(角度 (deg)); ylabel(MUSIC谱 (dB)); title(均匀线阵 MUSIC 测向仿真); grid on;跑完这段代码你能在-20度和10度附近看到两个清晰的谱峰。如果谱峰不明显优先检查A和A_scan的符号约定是否一致、En是否取对信号子空间和噪声子空间有没有搞反以及K的取值是否正确。4.4 用阵列响应矩阵画方向图对比画方向图其实也是调用阵列响应矩阵。比如我想看波束指向0度时的阵列方向图theta_scan -90:0.1:90; A_scan exp(-1j * 2 * pi * d * (0:M-1). * sin(deg2rad(theta_scan)) / lambda); w exp(-1j * 2 * pi * d * (0:M-1). * sin(deg2rad(0)) / lambda); % 指向0度的加权向量 pattern w * A_scan; pattern_dB 20 * log10(abs(pattern) / max(abs(pattern))); figure; plot(theta_scan, pattern_dB, LineWidth, 1.2); xlabel(角度 (deg)); ylabel(归一化方向图 (dB)); title(ULA 8元阵方向图波束指向0度); ylim([-40, 0]); grid on;这套代码稍加改动就能给任意阵列用只需要把每个角度下算导向矢量的方式和坐标匹配好。我强烈建议你把这个方向图打印出来和书上的图对照一下数值对得上再去跑后面的DOA算法会省掉很多排查时间。5. 常见问题与排查技巧实录我在调阵列响应矩阵相关代码时踩过不少坑这里整理几条最典型的希望能帮大家少走弯路。5.1 方向图左右颠倒或谱峰位置不对这是最常见的符号问题。有的公式写exp(-j * 2π * Δ / λ)有的写exp(j * 2π * Δ / λ)。两种写法都是自洽的关键在于发射信号模型、阵列响应矩阵、加权向量三者必须统一。如果方向图整体左右翻转大概率是符号约定不统一。排查方法是把来波方向设到30度看看波束主瓣是出现在30度还是-30度。5.2 栅瓣满天飞栅瓣的根源几乎都是阵元间距d λ/2。快速验证办法把d改成0.5看栅瓣是否消失。如果项目要求必须大间距那就得仔细分析哪些角度会模糊并在DOA估计阶段做解模糊处理。5.3 角度定义和坐标系混乱均匀线阵相对还好到了平面阵列和L型阵列方位角、俯仰角定义稍微不同导向矢量公式就完全变了。我建议每个工程目录里放一个“坐标系定义.md”写清楚x轴指向哪、y轴指向哪、方位角从哪个轴开始转、俯仰角是仰角还是俯视角。每次写新函数前先对照这个定义。5.4 圆阵MUSIC失效或不稳定圆阵直接套MUSIC需要满足阵列流型矩阵满秩的条件。解决思路是先用相位模式变换把圆阵变成虚拟的均匀线阵结构再跑MUSIC。如果不做变换圆阵半径要尽量接近λ/2并且目标个数不要太多否则谱峰会出现偏移。5.5 多径相干信号下MUSIC失效当两个信号完全相干时协方差矩阵的信号子空间会退化MUSIC没法分辨。这个问题的标准解法是空间平滑Spatial Smoothing也就是把均匀线阵划分成若干重叠子阵用子阵协方差矩阵的平均来“去相关”。代码上就是循环截取子阵数据取平均后再做特征分解。问题现象可能原因排查/解决建议方向图左右翻转导向矢量符号约定不一致统一使用exp(-j2πΔ/λ)并检查坐标定义出现等高的额外峰值阵元间距过大导致栅瓣将d改为λ/2验证如需大间距则设计稀疏阵列MUSIC谱峰偏移有限快拍导致协方差估计不准增大快拍数或进行对角线加载圆阵MUSIC失效阵列流型矩阵奇异使用相位模式变换或调整半径/阵元数两个相干信号分不开信号协方差矩阵秩亏使用前后向空间平滑预处理我把这些坑总结成一张速查表调试的时候先对照看一遍很多时候能直接定位问题。结尾最后分享一点个人体会阵列响应矩阵这东西看着是纯数学定义实际上和阵列几何、坐标系、波长、阵元位置这些硬件参数紧紧绑在一起。我刚学的时候以为记住a exp(-j * 2π * p * u / λ)就完事了后来发现真正困难的是把所有参数的单位、角度定义、坐标方向统一好让代码生成的结果和物理实际对得上。建议你先拿均匀线阵把整套仿真链路跑通再逐步扩展到圆阵、L型阵列和平面阵列最后把通用的任意阵列函数写出来。这样每一步都有对照验证出问题能快速定位。后面如果你开始做子阵划分、互耦校正或者稀疏阵列设计回头再看这个基础模型会有完全不一样的理解。