简介这是一份面向雷达通信与阵列信号处理学习者的Matlab仿真资源专注于相控阵雷达方向图的绘制与分析。资源包含可运行的.m源码、对应的.fig交互图以及3张运行结果截图方便读者直观对比仿真输出与理论结果。压缩包内共5个文件整体仅287KB适合本科、硕士阶段课程设计或科研入门时快速上手也适用于智能优化、信号处理等方向的Matlab仿真爱好者参考。目前已有223人学习使用代码结构简洁运行环境支持Matlab 2014/2019a有助于读者理解相控阵天线方向图形成原理、阵元布阵与波束指向等关键概念并可直接扩展用于更复杂的阵列仿真任务。1. 相控阵方向图仿真从公式到Matlab代码的最后一公里雷达系统工程师在评估相控阵天线时最关心的三个指标通常是主瓣宽度、旁瓣电平和波束扫描能力。教材上给出的方向图公式很简洁但真正动手用Matlab计算时阵元间距的选取、角度网格的划分、加权窗的引入方式每一处细节都会让最终画出来的方向图形状产生肉眼可见的差异。这套源码包里的xiangkongzhen1.m实现的是一套完整的相控阵方向图计算流程从阵列几何建模到方向图绘制一步到位运行结果中已经包含了三张不同参数下的方向图输出。适合本科和硕士阶段做雷达通信课程设计、毕业设计时快速验证阵列参数也适合刚接触阵列信号处理的工程师把教材公式落成可运行的代码。2. 从阵列流型到方向图函数相控阵波束成形的数学基础2.1 均匀线阵的远场方向图推导相控阵方向图的本质是各阵元辐射场在远场某一点的矢量叠加。考虑一个由N个各向同性阵元组成的均匀线阵阵元间距为d当平面波以入射角θ到达阵列时相邻阵元之间的波程差为d·sinθ对应的相位差为2π·d·sinθ/λ。如果阵列采用移相器对各阵元信号进行相位补偿使波束指向θ₀方向那么阵列的阵列因子可以写成AF(θ) Σ wₙ · exp(j·2π·n·d·(sinθ - sinθ₀)/λ)n 0, 1, ..., N-1其中wₙ是第n个阵元的幅度加权系数。当wₙ全为1时上式是一个等比数列求和可以闭式化简为sin(Nx)/sin(x)的形式这就是均匀加权的方向图函数。峰值出现在θ θ₀处第一零点位置由N·π·d·(sinθ - sinθ₀)/λ ±π决定因此主瓣宽度与N·d成反比。这个推导听起来简单但落到代码里有一个容易被忽略的细节Matlab的exp函数接受的是弧度制相位而方向图横轴通常习惯用角度显示。源码里如果直接把角度带入指数函数画出来的方向图会出现严重的畸变。正确做法是先建立角度网格再用deg2rad转换后再进入复数运算。2.2 阵元间距与栅瓣的边界条件阵元间距d的选择直接决定方向图是否会出现栅瓣。栅瓣的出现条件是当sinθ在[-1, 1]范围内时存在除主瓣方向外的其他角度使得各阵元相位差为2π的整数倍。数学上这个条件等价于d/λ ≥ 1/(1sinθ₀)。对于扫描到60°的极端情况d/λ必须小于约0.536才能避免栅瓣如果只做正侧向扫描d/λ 1即可保证不出现栅瓣。工程上通用的做法是取d λ/2这样在[-90°, 90°]的整个可视范围内都不会出现栅瓣同时阵元间距足够大阵元间的互耦影响也相对可控。xiangkongzhen1.m中如果默认使用的是半波长间距那么在修改参数时需要注意一旦把阵元间距调大方向图上会在偏离主瓣的位置冒出等高的虚假波峰这不是代码错误而是阵列设计本身违反了空间采样定律此时应该缩小d或者限制最大扫描角。2.3 方向图函数在Matlab中的向量化表达方向图计算的效率瓶颈在于N个阵元和M个角度网格的双重循环。对于常见的32阵元、1801个角度采样点0.1°分辨率双重循环需要约5.7万次复数指数运算Matlab跑起来并不慢但代码不够优雅而且循环版本在参数扫描场景下比如连续计算50组不同阵元数的方向图会明显变慢。向量化的思路是用一个N×M的二维矩阵一次性计算所有阵元和所有角度的相位延迟% 角度网格从-90到90度共1801个点 theta linspace(-90, 90, 1801); theta_rad deg2rad(theta); % 阵元位置N个阵元间距d单位波长 N 32; d 0.5; % 半波长间距 positions (0:N-1) * d; % 相位差矩阵size N x M phase_diff 2 * pi * positions * sin(theta_rad); % 导向矢量每列对应一个角度 steering exp(1j * phase_diff); % 波束指向30度时的加权向量 theta0_rad deg2rad(30); w exp(1j * 2 * pi * positions * sin(theta0_rad)); % 阵列因子对每个角度做加权求和 AF w * steering; AF_db 20 * log10(abs(AF) / max(abs(AF)));这段代码的关键在于positions * sin(theta_rad)用的是矩阵乘法而非逐元素乘法positions是N×1列向量sin(theta_rad)是1×M行向量两者相乘得到一个N×M矩阵矩阵的每一个元素(n, m)恰好是第n个阵元在角度theta(m)方向上的空间相位延迟。w是加权向量的共轭转置维度为1×N与N×M的steering矩阵相乘得到1×M的阵列因子向量。最后一步用max(abs(AF))归一化后转成dB单位保证方向图峰值在0 dB处。这段向量化代码在Matlab 2014和2019a上都能直接运行老版本不需要任何额外工具箱只需要基础MATLAB环境即可。3. xiangkongzhen1.m关键代码拆解阵列构建到方向图绘制全流程3.1 主程序结构与参数定义段xiangkongzhen1.m的主程序结构通常分为参数定义、阵列构建、方向图计算和结果可视化四段。参数段是调试时最常改动的部分需要明确每个变量的物理含义和单位参数变量典型值物理含义修改时的注意事项N16/32/64阵元数量增大则主瓣变窄计算量线性增加d0.5阵元间距波长倍数超过0.5且大角度扫描时会出现栅瓣theta00或30波束指向角度超过60°时主瓣明显展宽winones(1,N)幅度加权窗加窗后主瓣变宽旁瓣降低参数定义段在代码里通常集中排列在文件头部方便批量修改。一个常见问题是源码里可能直接写死了theta -90:0.1:90的角度网格如果改成linspace(-90,90,181)虽然采样点数一样但网格的分布方式不同——后者是闭区间均匀采样端点恰好包含±90°而前者是步进采样最后一个点可能是89.9°。这个差别在方向图横轴标注时几乎看不出影响但在计算第一零点位置或做高精度扫描角分析时闭区间采样的端点效应会更明显。3.2 阵列响应矩阵的构建方式主程序里最核心的一段是生成方向图数据。常见的实现方式有两种第一种是逐角度遍历对每个角度单独计算所有阵元的合成场强第二种是向量化的矩阵乘法方式。逐角度遍历的代码更贴近教材公式theta -90:0.1:90; AF zeros(size(theta)); for idx 1:length(theta) % 当前角度下的相位延迟向量 phase_delay 2 * pi * d * sin(deg2rad(theta(idx))) .* (0:N-1); % 阵列因子加权向量与导向向量的内积 AF(idx) sum(win .* exp(1j * phase_delay)); end % 归一化并转dB AF_db 20 * log10(abs(AF) / max(abs(AF)));phase_delay的计算中(0:N-1)是阵元索引向量sin(deg2rad(theta(idx)))是当前扫描角度下的正弦值两者逐元素相乘得到每个阵元的空间相位延迟。win如果全为1则等效于均匀加权如果换成切比雪夫窗则得到低旁瓣方向图。循环遍历方式的好处是调试方便可以在AF(idx)这一行设置断点检查任意角度下阵列因子的实部和虚部是否合理。向量化版本在第二节已经给出两者的计算结果完全一致。需要特别注意的是exp(1j * phase_delay)中phase_delay是长度为N的行向量exp函数返回同样长度的复数向量sum对这个向量求和得到该角度下的阵列因子复数值。如果误用了sum(exp(1j * phase_delay), 2)在行向量上不会有区别但如果代码里把(0:N-1)写成(0:N-1)变成列向量sum的方向就会出错导致结果变成每个阵元单独求和后的拼接方向图形状完全错误。这是移植代码时最常见的隐性bug。3.3 方向图绘制与坐标轴设置绘制方向图时线性和极坐标两种展示方式各有适用场景。极坐标图能直观展示方向图的覆盖范围适合展示波束扫描效果直角坐标图方便读取旁瓣电平的具体数值适合做参数对比% 直角坐标方向图 figure; plot(theta, AF_db, b-, LineWidth, 1.5); xlabel(角度 (deg)); ylabel(归一化幅度 (dB)); grid on; ylim([-60, 5]); xlim([-90, 90]); % 极坐标方向图角度转弧度后绘制 figure; polarplot(deg2rad(theta), max(AF_db, -60), r-, LineWidth, 1.5); title(极坐标方向图);极坐标绘制时特别注意polarplot函数要求数据非负所以需要把低于-60 dB的部分截断到-60 dB否则极小值会在极坐标图上产生大量毛刺。ylim设为[-60, 5]是因为方向图归一化后峰值在0 dB向下留60 dB的动态范围足够观察主瓣和旁瓣的相对关系再低的电平在实际雷达系统中会被接收机噪声底限淹没显示了也没有工程意义。4. 方向图参数调优阵元数、窗函数与扫描角的取舍4.1 阵元数对波束宽度与增益的影响阵元数N直接决定阵列孔径长度L (N-1)·d。对于均匀线阵半功率波束宽度近似为0.886·λ/(N·d·cosθ₀)弧度换算成角度后在侧射方向θ₀0°有BW ≈ 50.8/(N·d/λ)度。以d0.5λ为例16阵元时半功率波束宽度约为6.4°32阵元时约3.2°64阵元时约1.6°。这个关系在雷达设计中的含义是要获得1°的波束宽度在X波段波长约3cm需要约48个阵元对应的阵列物理长度约72cm。修改阵元数时除了主瓣变窄还需要关注增益变化。均匀线阵的定向增益近似为G ≈ 2·N·d/λ线性值dB形式为10·log10(2·N·d/λ)。阵元数从16增加到32增益理论上应增加约3 dB。如果仿真结果中增益变化不符合这个规律优先检查归一化方式方向图归一化用的是max(abs(AF))此时峰值始终为0 dB增益变化不会显式体现出来除非改用abs(AF)的绝对幅度来对比。4.2 窗函数加权切比雪夫窗与泰勒窗的选择均匀加权方向图的第一旁瓣电平约为-13.2 dB这在很多雷达场景下不满足低旁瓣要求。加窗是工程上最直接的旁瓣抑制手段。Matlab的signal工具箱提供了多种窗函数但对于阵列方向图加权最常用的是切比雪夫窗Dolph-Chebyshev和泰勒窗Taylor% 切比雪夫窗旁瓣电平可控但阵元数较大时边缘阵元激励剧烈跳变 N 32; SLL -30; % 目标旁瓣电平单位dB win_cheb chebwin(N, abs(SLL)); % 泰勒窗旁瓣电平从近区到远区逐渐衰减激励分布更平滑 nbar 4; % 等旁瓣区域内的零点数 win_taylor taylorwin(N, nbar, abs(SLL)); % 对比加窗前后的方向图 d 0.5; theta linspace(-90, 90, 1801); theta_rad deg2rad(theta); positions (0:N-1) * d; steering exp(1j * 2 * pi * positions * sin(theta_rad)); AF_uniform sum(steering, 1); AF_cheb win_cheb * steering; AF_taylor win_taylor * steering; norm_dB (x) 20 * log10(abs(x) / max(abs(x))); figure; plot(theta, norm_dB(AF_uniform), k-, LineWidth, 1); hold on; plot(theta, norm_dB(AF_cheb), r--, LineWidth, 1.2); plot(theta, norm_dB(AF_taylor), b-., LineWidth, 1.2); xlabel(角度 (deg)); ylabel(归一化幅度 (dB)); ylim([-60, 5]); legend(均匀加权, 切比雪夫窗 (-30dB), 泰勒窗 (-30dB)); grid on;chebwin(N, abs(SLL))中第二个参数必须是正数表示旁瓣电平的绝对值。taylorwin(N, nbar, abs(SLL))的nbar控制等旁瓣区的零点数量减小nbar会让方向图在远旁瓣区更快衰减但近区旁瓣会略高于设定值。运行这段代码后会看到切比雪夫窗方向图的所有旁瓣都严格压在-30 dB以下但主瓣宽度比均匀加宽了约1.4倍泰勒窗方向图的近区旁瓣接近-30 dB远区旁瓣逐渐降低主瓣展宽比切比雪夫窗略小工程上更常使用。4.3 波束扫描时的主瓣展宽与增益损失当波束从法线方向扫描到大角度时投影孔径变小主瓣按1/cosθ₀的关系展宽同时增益按cosθ₀的规律下降。仿真中设置theta0 60时方向图主瓣宽度会比侧射时展宽约2倍峰值增益下降约3 dB。这个现象在xiangkongzhen1.m中可以通过改变波束指向角参数直接观察到。figure; for theta0 [0, 30, 60] w_scan exp(1j * 2 * pi * positions * sin(deg2rad(theta0))); AF_scan w_scan * steering; AF_scan_dB 20 * log10(abs(AF_scan) / max(abs(AF_scan))); plot(theta, AF_scan_dB, LineWidth, 1.2, DisplayName, sprintf(\\theta_0 %d°, theta0)); hold on; end xlabel(角度 (deg)); ylabel(归一化幅度 (dB)); ylim([-60, 5]); xlim([-90, 90]); legend show; grid on;扫描到60°时主瓣峰值位置虽然正确落在60°方向但主瓣形状明显不对称——靠近法线一侧的旁瓣被压缩远离法线一侧的旁瓣抬升。这是因为均匀线阵在扫描状态下的有效孔径是N·d·cosθ₀主瓣两侧的零点位置不再对称。如果代码里对每个扫描角都独立做了max归一化那么扫描方向图之间无法直接比较增益差异要比较扫描损耗必须用一个固定的参考值比如侧射时的峰值幅度做归一化。5. 快速验证方向图代码正确性的几个实用技巧5.1 角度网格分辨率与栅瓣伪影的识别角度网格的步长如果太粗方向图的主瓣峰值可能落在网格点之间导致画出来的主瓣顶部是平的看起来像出现了分裂波束。检查方法很简单把网格步长从1°改成0.1°如果主瓣形状出现明显变化说明原来的分辨率不足。对于N32、d0.5λ的阵列主瓣宽度约3.2°至少需要0.1°的网格步长才能分辨主瓣轮廓对于N128的阵列主瓣宽度不到0.8°建议用0.01°的网格步长。在极坐标绘制时只要数据是线性间隔的polarplot会自动在相邻采样点间做插值但插值会产生虚假的方位角——如果主瓣峰值恰好落在两个采样点之间极坐标图上的主瓣看起来会比实际更宽。5.2 远旁瓣区域的数值噪声处理方向图的动态范围可以到-100 dB以下但在这个电平上abs(AF)的值接近epsMatlab浮点精度约2.2e-16转换到dB后会出现非物理的剧烈跳动。如果AF矩阵中存在微小的数值误差20*log10之后在-80 dB以下会出现大量锯齿状的伪旁瓣。处理方式有两个第一是在绘制时截断动态范围用max(AF_db, -60)限制显示下限第二是在计算过程中使用双精度复数运算避免中间结果被截断成单精度。5.3 用解析解验证数值计算均匀加权线阵在侧射状态下的方向图有解析表达式AF(θ) sin(N·π·d·sinθ/λ) / (N·sin(π·d·sinθ/λ))。在完成代码后可以用这个公式做一次快速验证检查数值计算的第一零点位置和旁瓣电平是否与理论值吻合。第一零点出现在sinθ λ/(N·d)处对于N16、d0.5λθ ≈ 7.18°第一旁瓣电平的理论值为-13.26 dB。如果代码计算结果与这两个值不匹配优先检查阵元间距是否用了波长归一化以及角度转换过程中deg2rad和rad2deg是否配对使用。% 解析解对比验证 N 16; d 0.5; theta_test linspace(-90, 90, 3601); theta_rad deg2rad(theta_test); % 解析公式侧射均匀加权 x pi * d * sin(theta_rad); AF_analytic sin(N * x) ./ (N * sin(x) eps); % 与数值计算结果对比 positions (0:N-1) * d; steering_test exp(1j * 2 * pi * positions * sin(theta_rad)); AF_numeric sum(steering_test, 1); AF_numeric_norm AF_numeric / max(abs(AF_numeric)); max_diff max(abs(abs(AF_analytic) - abs(AF_numeric_norm))); fprintf(最大偏差%e\n, max_diff);sin(N*x)在x接近0时的分子和分母都趋向0直接相除会出现0/0的NaN值。处理方式是在分母上加一个极小值eps或者对x0附近的点单独赋值1。这段验证代码跑通后可以确认主程序的方向图计算核心逻辑没有相位或坐标上的错误之后再调试加窗、扫描等扩展功能时每次改动都能用这个验证脚本做回归测试。本文还有配套的精品资源点击获取