线阵、面阵、圆阵方向图仿真：从导向矢量到MATLAB实现

简介面向信号处理与无线通信领域学习者的一份小型MATLAB代码包内含线阵、面阵与圆阵三种典型天线阵列的方向图计算脚本适合初学者理解阵列因子、相位控制与辐射方向图生成原理。压缩包共4个文件其中3个为.m脚本分别对应均匀线阵、均匀面阵和均匀圆阵方向图另附1个txt说明文件整体体积仅2KB代码精简、便于快速编译和修改。已有553人学习过该资源从实践角度看用户可通过运行脚本直观对比不同阵列布局的方向性差异并在此基础上调整天线间距、排列方式或相位分布进一步掌握MATLAB在阵列天线建模与分析中的实际用法。配合描述中关于线阵、面阵、圆阵的知识点可辅助雷达、无线通信与声学应用场景的课程教学或入门研究。1. 你拿到的线阵、面阵、圆阵方向图代码问题出在哪写阵列方向图仿真很多人拿到 Patern.rar 这类代码包之前觉得无非是把公式代进 MATLAB最后 plot 一张图。真打开代码会发现线阵脚本用 sin(theta) 循环面阵脚本用 meshgrid 加二维求和圆阵脚本又是一套 cos(phi - phi_n) 的逻辑三个文件各自为政改了阵元数就出栅瓣改了波束指向就找不到峰值。方向是对的但哪一行代码对应哪一项物理量经常没人讲透。这篇按线阵、面阵、圆阵三条线走先把共用的导向矢量模型立住再给最小可运行代码和必调参数最后把三类脚本收进同一个阵列流形函数。适合雷达、通信、声呐、导航方向手上有代码但改不明白或者准备从头做阵列仿真的工程师。2. 方向图仿真的共同底座导向矢量、阵因子与 MATLAB 验证2.1 方向图先拆成“阵因子 × 单元因子”再谈代码在写任何一行 MATLAB 之前先确认一件物理上的事总方向图是单元因子与阵因子相乘。单元因子描述单个阵元的辐射模式阵因子描述阵元之间的相位叠加两者分开处理能省很多事。很多代码包默认单元是理想各向同性点源单元因子等于 1方向图就只取决于阵因子。真实工程里若用微带贴片、偶极子或喇叭单元因子是随角度变化的仿真脚本里要额外乘上对应方向图但坐标系和角度语义往往从这里开始跑偏。阵因子的计算本质是对每个扫描方向累加各阵元贡献AF(θ) Σ w_n exp(j·2π/λ·p_n·u(θ))。其中 p_n 是第 n 个阵元的位置u(θ) 是扫描方向的单位矢量w_n 是复权值决定波束指向与副瓣电平。这句话是整篇的地基后续线阵、面阵、圆阵都只是把 p_n 放进不同类型坐标里。搞不清这一步后面改参数时就会反复踩同一类坑以为在调副瓣实际在调相位参考点。2.2 线阵、面阵、圆阵的导向矢量差异对照三类阵列的差异集中在阵元坐标和导向矢量表达式上。把这三行表达式摆在一起看比分别翻三个脚本更能建立整体感。阵列类型阵元位置 p_n导向矢量对应扫描方向 θ/φ自由度特点常用 MATLAB 对象均匀线阵 ULAd·(n-1)·x̂exp(j·2πd·(n-1)·sinθ)只在一条线上控制波束phased.ULA矩形面阵 URA(m·dx, n·dy, 0)exp(j·2π(dx·m·u dy·n·v))方位、仰角独立可控phased.URA均匀圆阵 UCA(R·cosφ_n, R·sinφ_n, 0)exp(j·2πR·sinθ·cos(φ-φ_n))方位无边界360° 各向扫描phased.UCA表里 u sinθ·cosφv sinθ·sinφθ 从阵面法向z 轴算起φ 从 x 轴算起。这个角度约定在 2.3 和第四章代码里保持一致否则面阵和圆阵的方向图会出现镜像或旋转偏移。线阵那一行里 sinθ 直接描述了波程差随角度变化的全部信息面阵用两个方向余弦 u、v 把二维波束投影到 x、y 两轴圆阵则多了 cos(φ-φ_n)这个余弦项是圆阵区别于另外两类阵列的关键第五章会专门展开。2.3 用 ArrayResponse 验证公式三行出图如果 MATLAB 装了 Phased Array System Toolbox最快验证公式的方法是直接拿系统对象算不用手写循环。以 8 元半波长间距线阵为例下面代码得到等幅加权的阵列响应并与理论公式交叉验证。% 用 Phased Array System Toolbox 验证 ULA 方向图公式 fc 3e9; % 载波频率单位 Hz lambda physconst(LightSpeed) / fc; % 波长 ula phased.ULA(NumElements, 8, ... ElementSpacing, lambda/2); % 半波长间距 arr phased.ArrayResponse(SensorArray, ula, ... WeightsInputPort, true); % 打开权重输入口 ang -90:0.1:90; % 扫描方位角单位度 w ones(8, 1); % 等幅权重 resp arr(fc, ang, w); % 计算阵列响应 plot(ang, mag2db(abs(resp) / max(abs(resp)))); grid on;这里 ang 可以传一行向量也可以传 2×M 矩阵第二行放俯仰角线阵只传方位角即可。WeightsInputPort 设为 true 后第三个入参是复权重向量后面做切比雪夫锥削或圆阵相模法时都要靠这个口。输出 resp 是复数画图前要取幅值再归一化不然直接 mag2db 会得到一堆波纹。这个对象的价值是当参照物自己手写代码时如果两端结果对不上优先检查波束指向处的符号和相位参考点而不是怀疑公式本身。3. 线阵方向图从零手写最小代码盯住栅瓣和波数图3.1 最小 ULA 方向图代码不依赖工具箱不用工具箱12 行以内就能算完 8 元线阵的方向图。这一步的价值在于把相位参考点、共轭转置、sinθ 这些细节落回代码本身。% ULA 方向图等幅加权不依赖任何工具箱 N 8; % 阵元数 d 0.5; % 阵元间距单位波长 theta0 30; % 波束指向相对阵面法向单位度 theta -90:0.1:90; % 扫描范围单位度 s exp(1j * 2*pi * d * (0:N-1). * sin(theta*pi/180)); % s: N×L每列是某个扫描角的导向矢量 b exp(1j * 2*pi * d * (0:N-1). * sin(theta0*pi/180)); % 目标方向的导向矢量 AF b * s; % 阵列响应共轭转置让峰值对准 theta0 AF abs(AF) / max(abs(AF)); plot(theta, mag2db(AF)); grid on;核心在第 6 行的外积展开向量 (0:N-1). 与 sin(theta*pi/180) 做矩阵乘得到的每一列就是该扫描角下 N 个阵元的相位因子。第 8 行的 b 是波束指向处的导向矢量b * s 相当于把每个扫描方向与目标方向做相关峰值自然落在 theta0。若把 b 改成 b.不带共轭峰值会变成 -30 度这就是相位参考点选错的典型症状。代码本身简单但把它改成非等幅加权、改成俯仰角、改成复数权值时很多人会在这个符号上卡住。3.2 线阵的三个必调参数阵元数、间距、波束指向线阵方向图调参绝大多数场景跑不出这三个参数阵元数 N、阵元间距 d/λ、波束指向 θ0。它们各自的作用边界值得列出来参数对方向图的影响改坏的典型症状常用取值阵元数 N决定主瓣宽度与阵列增益N 太小主瓣宽、副瓣粗糙8~64间距 d/λ决定栅瓣位置与有效孔径d/λ 超过 0.5 后可见区内出现栅瓣0.3~0.5波束指向 θ0决定主瓣位置与展宽程度偏离法向过远时主瓣畸变±60° 以内栅瓣出现条件常被简写成 d/λ ≤ 0.5这只在 θ00 时成立。实际判断要用完整形式sinθ_g sinθ0 - λ/d。当 θ_g 落入 [-1,1] 区间可见区里就出现第二个峰值。所以做 30° 波束指向时d/λ0.5 已经离边界不远指向到 60° 时即使 d/λ0.5 也有栅瓣风险。这也是为什么很多雷达阵列的阵元间距取 0.4λ 而不是贴着 0.5λ 用。副瓣控制通常靠窗函数。线阵里最常用的是切比雪夫窗和 Taylor 窗前者把副瓣压到指定电平但主瓣会按比例拓宽。切换方法是在 3.1 代码里把 b 替换成 w.*b其中 w chebwin(N, 30)30 表示副瓣目标 -30dB。注意窗函数必须作用在阵元域而不是加在方向图上否则只会改变曲线的显示形状物理上完全无效。3.3 波数图盯栅瓣横轴换成 usinθ混叠藏不住角度域方向图有个缺点sin 函数把角度区域的间距压缩到两端栅瓣在 ±90° 附近时从方向图上很难判断它是不是刚进入可见区。把横轴换成 usinθ 画波数图这个问题就解决了。% 波数图画法横轴为 usin(theta) u -1:0.001:1; AF_u b * exp(1j * 2*pi * d * (0:N-1). * u); plot(u, mag2db(abs(AF_u) / max(abs(AF_u)))); grid on; xlabel(u sin\theta);波数图里峰值严格周期分布周期是 Δu λ/d。d/λ0.5 时 Δu2可见区 [-1,1] 内恰好放一个完整周期d/λ0.8 时 Δu1.25第二个峰必然挤进可见区。网上搜 ula 波数图 matlab 能拿到不少画法片段但大多数只给 plot 一行没解释这个周期关系。自己画时留意一点周期只与 d/λ 有关与 N 无关N 只决定每个峰的宽度和副瓣起伏。判断栅瓣先看波数图再看角度图顺序反过来容易误判。4. 面阵方向图矩形阵拆成两个线阵再合成二维波束4.1 矩形面阵方向图等于两个线阵方向图的乘积矩形栅格面阵有一个重要性质导向矢量可分离。x 方向 Nx 个阵元、y 方向 Ny 个阵元合成导向矢量是 a_x ⊗ a_y阵列响应也就能拆成两个线阵响应相乘AF(θ,φ) AF_x(u) · AF_y(v)其中 u sinθ·cosφv sinθ·sinφ。这个性质直接决定了面阵代码的写法不要写双层循环做二维求和先分别算两个一维响应再按元素相乘。下面代码以 8×8 面阵为例主波束指向仰角 30°、方位角 45°。% 矩形面阵方向图可分离写法AF AFx .* AFy Nx 8; Ny 8; dx 0.5; dy 0.5; % x、y 方向间距单位波长 theta0 30; phi0 45; % 波束指向度theta 从法向算起 theta 0:0.5:90; phi 0:1:360; [TH, PH] meshgrid(theta, phi); u sin(deg2rad(TH)) .* cos(deg2rad(PH)); v sin(deg2rad(TH)) .* sin(deg2rad(PH)); u0 sind(theta0)*cosd(phi0); v0 sind(theta0)*sind(phi0); sx exp(1j*2*pi*dx*(0:Nx-1). * u(:).); % Nx × L sy exp(1j*2*pi*dy*(0:Ny-1). * v(:).); % Ny × L wx exp(-1j*2*pi*dx*(0:Nx-1). * u0); % x 向波束指向权 wy exp(-1j*2*pi*dy*(0:Ny-1). * v0); % y 向波束指向权 AF (wx * sx) .* (wy * sy); % 可分离按元素乘 AF reshape(abs(AF), size(TH)); AF_db mag2db(AF / max(AF(:)) eps); surf(TH, PH, AF_db, EdgeColor, none); colormap(jet); colorbar; view(2); xlabel(仰角 theta度); ylabel(方位角 phi度);这段代码的巧妙之处在于把二维问题降成了两次一维矩阵乘法。sx 是 Nx×Lsy 是 Ny×LL 是所有角度网格点数8×8 阵元算 0.5°×1° 的全空间网格在 MATLAB 里只要几十毫秒。如果按老写法对每个角度 kron 一次再循环同样的网格可能要等好几秒而且代码一长错的地方不好找。4.2 面阵代码中角度定义与网格化的三个坑第一个坑是角度定义。MATLAB 相控阵工具箱的方位角/仰角惯例是天线上常用的 az/el 定义方向余弦是 cos(el)·sin(az) 与 sin(el)而许多教材和代码包从球坐标抄来的是俯仰角 θ 从法向算起的定义。两套定义混用时同一个 (30°, 45°) 会画出完全不同的方向图甚至出现两个主瓣。拿到别人的 Patern.rar 这类代码包第一件事是看文件头注释用的哪套约定而不是先跑图。第二个坑是网格疏密。画花瓣图 0.5° 步进够用但扫描范围到全空间时网格点数量会超过十万。此时角度域方向图容易在副瓣区域出现孤立的深零点画 surf 时表现为一条条细缝看起来像错误实际上是网格分辨率问题。可以先画波数图或只切几个固定剖面确认结构后再上全 3D。第三个坑是 dB 显示的 -Inf。mag2db(0) 是负无穷AF 在很多方向等于零surf 会出现黑色空洞。上面代码里加了 eps 作保护或者用 max(AF, 1e-6) 截断。这个小问题在 3.1 的一维 plot 里不明显到了 surf 上非常扎眼第一次画面阵的人几乎都会遇到。4.3 从线阵迁移到面阵的改动清单如果手里只有线阵代码扩成面阵时按下面这个表逐项核对比对着报错信息猜要快迁移点线阵写法面阵写法阵元位置d*(0:N-1) 一维数组(mdx, ndy) 二维坐标导向矢量exp(j*2πdn·sinθ)exp(j*2π(dx·m·u dy·n·v))波束指向θ0 一个参数θ0 与 φ0 两个参数副瓣控制一维窗函数二维窗 两个一维窗的外积绘图plot 一维曲线surf view(2)二维窗的写法是 w2d (w_x * w_y)两个一维切比雪夫窗外积。若 dx 和 dy 不相等两个方向的主瓣宽度和栅瓣条件要独立判断x 方向看 Nxdx/λy 方向看 Nydy/λ。迁移时常犯的错是把 Nx 和 Ny 写反导致方向图两个轴的主瓣宽度互换这种错从图形上非常难一眼看出来。5. 圆阵方向图为什么不能照搬线阵以及 360° 波束5.1 圆阵的导向矢量先写出这一行相位项圆阵的每个阵元分布在半径为 R 的圆上第 n 个阵元的方位角是 φ_n 2πn/N。波束方向由仰角 θ 和方位角 φ 描述阵元位置与波束方向的波程差投影为 R·sinθ·cos(φ-φ_n)。这个余弦项是圆阵与线阵、面阵的本质区别——线阵和面阵的导向矢量都能写成两个独立方向的乘积圆阵的 sinθ 和 cos(φ-φ_n) 耦合在一起没法拆开。% 均匀圆阵方向图阵元位置在 x-y 平面 N 16; R 2; % 阵元数、圆半径单位波长 theta0 30; phi0 60; % 波束指向度 phi_n (0:N-1) * 2*pi / N; % 阵元方位角弧度 theta 0:1:90; phi 0:1:360; [TH, PH] meshgrid(theta, phi); w exp(-1j * 2*pi*R * cos(phi0 - phi_n) * sind(theta0)); % N×1 权重 AF zeros(size(TH)); for k 1:numel(TH) s exp(1j * 2*pi*R * cos(PH(k) - phi_n) * ... sin(TH(k)*pi/180)).; AF(k) w * s; end AF abs(AF) / max(abs(AF(:))); surf(TH, PH, mag2db(AF eps), EdgeColor, none); view(2); colorbar; xlabel(theta度); ylabel(phi度);循环写法每轮只处理一个扫描角N16、网格 91×361 时大概一秒多出结果可读性优于向量化。注意加权向量 w 的相位项与导向矢量 s 的相位项互为共轭符号写反会导致主瓣出现在 (-30°, 240°) 而不是 (30°, 60°)。当 N 比较大、网格更密时代价是循环次数上升这时再考虑把所有扫描角组装成矩阵一次求积但方向图仿真的瓶颈通常在 surf 渲染而非计算。5.2 均匀加权圆阵副瓣高相模法的改进与边界同一数量级阵元数下均匀加权的圆阵第一副瓣通常在 -8dB 上下明显高于线阵的 -13.3dB。原因是圆阵的投影孔径随方位角变化等幅激励时不同方向累积的能量分布不均匀。相模法phase mode是常用改进思路先把阵元域权重做 DFT 变换到模式域在模式域加窗后再变回阵元域。% 相模法阵元域 - 模式域 - 加窗 - 回阵元域 modes ifft(w, N); % 到模式域 win chebwin(N, 30); % 模式域锥削窗 modes_w modes .* win; % 模式域加权 w_w fft(modes_w); % 回到阵元域相模法能压低副瓣但它的有效模式数有限。阵列流型在模式域的秩取决于半径与波长之比可用模式数大致在 2πR/λ 附近超过这个带宽的模式增益很低直接参与波束形成会把主瓣拉宽。这意味着圆阵的可用扫描范围也存在边界接近水平面θ 接近 90°时方向图退化最明显。工程里做圆阵波束形成先用这段代码看模式域分布比直接在阵元域调权重更能解释问题。5.3 圆阵仿真中常被忽略的互耦因素相邻阵元间距 2πR/N。R2λ、N16 时间距约 0.785λ互耦尚可R1λ、N16 时只有约 0.39λ低于半波长互耦效应会明显抬高低仰角方向的副瓣。很多教学代码完全没有互耦模型方向图仿真结果在理想条件下成立拿到真实阵列上会对不上。更接近工程的做法是用互耦矩阵 C 修正权值或者直接把全波仿真得到的 S 参数折算成互耦校正系数但这超出普通阵列方向图代码包的范畴。从 Patern.rar 这类包里拿到圆阵脚本时先确认它有没有互耦相关的测试用例没有的话就清楚它给的是一个理想上界。6. 三个自检用例与一个统一工具函数6.1 三个自检用例写完代码后先用最小用例验证再跑全参数扫描。以下三个用例能覆盖绝大多数实现错误检查项最小用例期望结果主瓣位置N16、d/λ0.5、θ030角度图峰值在 30°±1°栅瓣边界d/λ0.8、θ00波数图可见区出现第二个峰圆阵旋转对称N16、R2、θ00、φ0 任意方向图不随 φ0 旋转其中第三项最容易暴露圆阵代码的符号问题。θ00 时波束指向法向圆阵的加权向量各阵元相位相同此时改变 φ0 方向图不变。如果结果显示方向图随 φ0 变化说明加权或导向矢量里混入了绝对相位参考。6.2 一个统一工具函数线阵、面阵、圆阵的差异只在阵元坐标和权值方向图计算本身有一个统一的数学形式AF(u) w · exp(j·2π/λ·pos·u)。可以写成通用函数function AF array_pattern(pos, lambda, w, u) % pos: 3×N 阵元坐标单位与 lambda 保持一致 % lambda: 载波波长 % w: N×1 复权重 % u: 3×M 扫描方向单位向量每列一个方向 AF w * exp(1j*2*pi/lambda * (pos * u)); end调用时按阵列类型生成 pos% 线阵阵元在 x 轴 pos_ula [d*(0:N-1); zeros(2,N)]; % 面阵阵元在 x-y 平面 [m, n] meshgrid(0:Ny-1, 0:Nx-1); pos_ura [dx*m(:).; dy*n(:).; zeros(1,Nx*Ny)]; % 圆阵阵元在 x-y 平面圆周上 pos_uca [R*cos(phi_n); R*sin(phi_n); zeros(1,N)]; % 扫描方向单位向量 u u [sin(theta(:)).*cos(phi(:)).; ... sin(theta(:)).*sin(phi(:)).; ... cos(theta(:)).*ones(1,numel(phi))]; AF array_pattern(pos_uca, lambda, w, u);用这个函数验证代码包时有一个好处三个脚本的坐标生成部分独立检查任何一个数组维度对不上报错信息会精确指向 pos 和 w 不匹配。最后提醒一句所有代码里都要统一单位pos 用波长归一化后lambda 直接传 1 即可这会省掉一批量级错误。下次拿到任何阵列方向图代码包把 pos 坐标列出来这个函数就是最快的验收入口。本文还有配套的精品资源点击获取