遗传算法在天线方向图综合与圆环阵列优化中的MATLAB实现
2026/9/16 13:23:26 网站建设 项目流程

简介:面向电磁仿真与天线设计方向的MATLAB开发者,这份资料聚焦用遗传算法优化36单元圆形天线阵列,内容涉及旁瓣电平与增益性能权衡,适合具备一定天线基础并希望掌握智能优化算法的进阶读者。压缩包共10个文件,包含2个m源码、1个txt优化系数、6个xls辐射性能表格及1个asv自动保存备份,m文件实现主优化算法与距离约束函数,xls记录多频点阵列因子结果,整体体积仅1.22MB。已有291人学习下载,说明该主题受到一定关注。利用该资源可直接运行遗传算法主程序,对照优化后系数文本和不同频率下的辐射性能表格,理解从适应度函数设计、选择交叉变异到迭代收敛的完整流程。对于想深入掌握MATLAB ga函数在天线布局优化中应用的研究者,这份轻量工程是较完整的示例参考。

1. 方向图综合的新解法:为什么把优化问题交给遗传算法

36 单元圆形天线阵列的优化,本质上是一个典型的多峰、非凸、高维参数寻优问题。传统方向图综合方法,比如泰勒综合、切比雪夫综合,设计思路大多是解析地控制旁瓣电平,方法是优美的,但一旦阵列几何变成圆环、需要考虑互耦修正和扫描角约束时,解析解便迅速失去掌控力。另一个极端是暴力扫参,但 36 个阵元的激励幅度和相位构成的设计空间,已经完全超出网格搜索的工程可行性。

遗传算法(GA)之所以成为这个问题的默认起点,不是因为它在所有场景下都最优,而是它具备两个关键特性,一是无需求导,方向图综合的目标函数经常带有不连续约束,比如要求主瓣内的波纹小于 0.5 dB、旁瓣区域必须低于 -20 dB,梯度类的 LMI 或凸优化很难把这些工程约束直接写成标准形式,而 GA 只关心适应度函数能不能算出来;二是在大搜索空间下能找到可接受解的性价比高,36 单元阵列的幅度和相位变量加起来至少有 36 维,GA 通过种群并行搜索,不需要任何初始猜测。

这篇文章就是顺着这条线往下走的。第一步先建立圆环阵列的方向图模型,第二步把优化目标翻译成 GA 能理解的适应度函数,第三步落在 MATLAB 的代码实现和参数校准上,最后讨论从仿真走向设计时必须做的几个修正。文中给的代码不加第三方工具箱,核心 GA 用 MATLAB 基础函数就能跑通,有 Global Optimization Toolbox 的可以对照着用ga原生接口提速。

2. 圆环阵列的阵列因子建模:从几何排布到方向图计算

2.1 36 单元均匀圆环阵列的几何与阵列因子公式

圆形天线阵列的定义很直接,36 个阵元等间距分布在半径为 r 的圆周上。单元 n 的角度位置为 φn = 2πn/36(n = 0,1,...,35)。假设阵列位于 x-y 平面,观察方向由仰角 θ 和方位角 φ 决定。总场方向图可以写作:

AF(θ,φ) = Σ_{n=1}^{36} I_n · exp( j·k·r·sinθ·cos(φ-φn) + j·βn )
  • I_n:第 n 个阵元的激励幅度
  • βn:第 n 个阵元的激励相位
  • k = 2π/λ:自由空间波数
  • r:阵列半径,通常用波长归一化表示,比如 r = 0.6λ 到 1.5λ

相比线阵,圆环阵列的阵列因子没有标准的 FFT 加速形式,方向图是方位角 φ 的周期函数。工程上最关心的两个设计自由度为:

  • 幅度分布 I_n,控制旁瓣和波束宽度
  • 相位分布 βn,控制波束指向

这两个参数全自由时,设计变量为 72 维。实际优化中常把幅度归一化到 [0,1](或等价的 dB 衰减),相位归一化到 [-π, π]。如果进一步考虑结构对称性,比如让激励幅度按环形周期重复,变量可以降到 18 维左右,收敛速度会有明显提升,但代价是可实现的方向图集会缩小。对于 36 单元的规模,一般不推荐强行降维,GA 的收敛速度完全能处理 72 维的搜索,细节会在第 4 章展开。

2.2 MATLAB 中高效计算阵列因子的实现

阵列因子计算是适应度函数的最内层循环,优化过程中会被调用数百次甚至上千次,它的效率直接决定一次优化要跑几分钟还是几小时。常见做法是先用矩阵运算一次性生成所有阵元在全部扫描角上的相位差项,而不是用 for 循环逐点累加。

function AF = circular_array_factor(I, beta, theta, phi, R, lambda) % I : 36x1 归一化激励幅度,范围 [0,1] % beta : 36x1 激励相位,单位是弧度 % theta : 1xNt 俯仰角采样向量 % phi : 1xNp 方位角采样向量 % R : 阵列半径,单位是米 % lambda : 工作波长,单位是米 N = length(I); % 阵元的方位角位置 phi_n = (0:N-1) * 2 * pi / N; % 1xN % 构造 [Nt*Np, N] 维的观察方向矩阵 [TH, PH] = meshgrid(theta, phi); TH = TH(:)'; % 1 x (Nt*Np) PH = PH(:)'; % 1 x (Nt*Np) k = 2 * pi / lambda; % phase_matrix 的每一行对应一个观察方向,每一列对应一个阵元 % 这里使用了广播机制,实际内存占用约为 (Nt*Np)*N 的复数矩阵 phase = zeros(length(TH), N); for n = 1:N phase(:, n) = I(n) .* exp(1j * (k * R * sin(TH(:)) .* cos(PH(:) - phi_n(n)) + beta(n))); end AF = sum(phase, 2); % 每个观察方向的总场 AF = reshape(AF, length(theta), length(phi)); end

逻辑说明:这个函数把空间角网格摊平成列向量,对每个阵元乘上幅度、加上相位旋转,再按列累加得到总场。MATLAB 对矩阵运算的优化远好于循环,这里虽然保留了一个for n = 1:N的循环,但循环次数只有 36,内层已经是向量化操作,性能足够。如果要进一步提速,可以把sincos的计算也矩阵化,不过 36×Nt×Np 的三角运算在 MATLAB 里并不会成为瓶颈。

参数说明:Rlambda的单位一致性很容易被忽略。若用波长归一化,直接令R=0.5, lambda=1即可,此时k*R等于 π,方向图的主瓣宽度和旁瓣分布都是相对于波长的,物理上更直观。实际工程中如果需要改变工作频率,把Rlambda都乘以同样的缩放因子不改变结果,所以建议在脚本里统一用归一化值。

2.3 优化目标怎么定:旁瓣电平、方向性和波束指向约束

方向图综合的目标函数设置直接影响优化结果的可落地程度。最常见的目标组合是:

  • 最小化最大旁瓣电平(MSLL)
  • 保证最大辐射方向对准预定角度
  • 可选地,约束主瓣波纹在 0.5 dB 以内

把这三个目标合成一个适应度函数时,通常用加权罚函数法,因为 GA 本身不擅长处理硬约束。常见做法是:MSLL 作为主目标,波束指向偏差作为罚项,波纹作为约束。

function fitness = fitness_function(x, theta_scan, phi_scan) % x 是待优化的染色体向量,前 36 个为幅度,后 36 个为相位 N = 36; I = x(1:N); beta = x(N+1:2*N); % 幅度归一化:防止 GA 把幅度整体缩小来压低旁瓣 I = I / max(I); % 目标方向:例如 theta=30度, phi=45度 theta0 = 30 * pi / 180; phi0 = 45 * pi / 180; theta_range = linspace(0, pi, 361); phi_range = linspace(0, 2*pi, 721); AF = circular_array_factor(I, beta, theta_range, phi_range, 0.8, 1); AF_dB = 20 * log10(abs(AF) / max(abs(AF(:))) + eps); % 主瓣区域外的最大值,即旁瓣电平 mainlobe_radius = 8 * pi / 180; % 主瓣半径,单位弧度 [TH_mesh, PH_mesh] = meshgrid(theta_range, phi_range); mainlobe_mask = zeros(size(AF_dB)); for idx = 1:length(theta_range) for jdx = 1:length(phi_range) % 计算与目标方向的大圆距离 cos_angle = sin(theta0)*sin(TH_mesh(idx,jdx))*cos(PHI_mesh(idx,jdx)-phi0) ... + cos(theta0)*cos(TH_mesh(idx,jdx)); angle_dist = acos(min(1, max(-1, cos_angle))); if angle_dist < mainlobe_radius mainlobe_mask(idx, jdx) = 1; end end end sll_region = AF_dB(~mainlobe_mask); MSLL = max(sll_region); % 方向性系数的粗估计(数值积分) theta_w = theta_range(2) - theta_range(1); phi_w = phi_range(2) - phi_range(1); D = 4*pi * max(abs(AF(:)))^2 / sum(sum(abs(AF).^2 .* sin(TH_mesh) .* theta_w .* phi_w)); % 适应度:最小化旁瓣和波束指向误差,方向性作为奖励项 lambda_sll = 1.0; lambda_steer = 0.5; lambda_d = 0.01; % 波束指向误差:计算最大辐射方向 [max_val, max_idx] = max(abs(AF(:))); [max_theta_idx, max_phi_idx] = ind2sub(size(AF), max_idx); steer_error = acos( min(1, max(-1, ... sin(theta0)*sin(theta_range(max_theta_idx))*cos(phi_range(max_phi_idx)-phi0) + ... cos(theta0)*cos(theta_range(max_theta_idx)) ))); fitness = lambda_sll * MSLL + lambda_steer * steer_error * 180/pi - lambda_d * D; end

逻辑说明:这个适应度函数返回的数值越小代表个体越优。MSLL 是主要惩罚项,波束指向误差则确保 GA 不会通过歪斜主瓣来压低旁瓣,方向性作为奖励项防止主瓣过度展宽换取低旁瓣。主瓣半径设在 8 度,用户可以按工作频段和实际波束宽度调整。

参数说明:三个 lambda 系数决定了优化偏好。旁瓣约束很紧的场合,把lambda_sll加大到 2 到 3;波束指向精度要求高的,lambda_steer上调;如果只关心旁瓣不关心增益,lambda_d可以直接设 0。注意幅度归一化那句I = I/max(I),它防止 GA 通过整体缩小激励幅度来作弊,因为方向图的旁瓣电平是相对值,整体缩小不改变形状,但归一化后可以保证至少有一个阵元达到满激励,使旁瓣电平的优化代价真实反映幅度动态范围。

3. 遗传算法与 MATLAB 中的完整实现:编码、算子与主循环

3.1 编码方案设计:幅度和相位的染色体排布

遗传算法的编码决定了搜索空间的结构。36 单元阵列的每个单元有幅度和相位两个变量,染色体总长度为 72。常见做法是实数编码,每一位是一个双精度浮点数。幅度变量限定在 [0, 1],相位变量限定在 [-π, π]。

lb = [zeros(1, 36), -pi * ones(1, 36)]; % 幅度下限0,相位下限 -pi ub = [ones(1, 36), pi * ones(1, 36)]; % 幅度上限1,相位上限 pi

实数编码的优势在于交叉和变异算子可以直接在连续域操作,GA 的搜索粒度不受二进制位数的限制。一些资料里的工作使用二进制格雷编码做离散幅度,比如 6 bit 对应 64 级衰减,那是为了匹配实际衰减器的量化步进。如果你的设计最终要用移相器和衰减器实现,建议在适应度函数里加一个量化步骤,把连续变量舍入到可用位数后再计算方向图,这种“先连续搜索,最后再量化”的路径往往在量化后性能严重退化,因为 GA 没有感知量化误差。

3.2 选择、交叉、变异算子在代码里的具体形式

MATLAB 的 Global Optimization Toolbox 提供了现成的ga函数,但为了讲清机制,也为了在没有工具箱的机器上能跑,这里给一套手写 GA 的核心代码。选择用锦标赛选择,交叉用模拟二进制交叉(SBX),变异用多项式变异。

function [pop, fitness] = init_population(popsize, dim, lb, ub) % 标准种群初始化:均匀随机采样 pop = repmat(lb, popsize, 1) + rand(popsize, dim) .* repmat(ub - lb, popsize, 1); fitness = zeros(popsize, 1); for i = 1:popsize fitness(i) = fitness_function(pop(i, :)); end end
function [child1, child2] = sbx_crossover(parent1, parent2, lb, ub, eta_c) % SBX 模拟二进制交叉,eta_c 是分布指数,一般取 20 dim = length(parent1); u = rand(1, dim); beta = zeros(1, dim); beta(u <= 0.5) = (2 * u(u <= 0.5)).^(1/(eta_c+1)); beta(u > 0.5) = 1 ./ (2 * (1 - u(u > 0.5))).^(1/(eta_c+1)); child1 = 0.5 * ((1 + beta) .* parent1 + (1 - beta) .* parent2); child2 = 0.5 * ((1 - beta) .* parent1 + (1 - beta) .* parent2); % 边界修补 child1 = min(max(child1, lb), ub); child2 = min(max(child2, lb), ub); end
function child = polynomial_mutation(parent, lb, ub, eta_m, pm) % 多项式变异,eta_m 通常取 20,pm 是变异概率 child = parent; dim = length(parent); for i = 1:dim if rand < pm r = rand; delta = zeros(1, dim); if r < 0.5 delta(i) = (2*r)^(1/(eta_m+1)) - 1; else delta(i) = 1 - (2*(1-r))^(1/(eta_m+1)); end child(i) = parent(i) + delta(i) * (ub(i) - lb(i)); end end child = min(max(child, lb), ub); end

逻辑说明:SBX 交叉不同于简单的算术交叉,它会以较大概率生成靠近父代的子代,同时保留小概率的远距离探索,这在连续空间中比均匀交叉收敛快。多项式变异类似地偏向小步扰动,保证种群在后期有精细搜索能力。选择、交叉、变异的组合顺序是:先锦标赛选择出popsize/2对父代,交叉生成popsize个子代,然后对每个子代按概率变异,最后父代和子代合并,按适应度排序取前popsize个进入下一代,也就是精英保留策略。

参数说明:eta_ceta_m控制子代与父代的相似度。数值越大,子代越接近父代,搜索越保守,适合后期精调;数值小则探索性强,适合前期全局搜索。实践中可以动态调整,比如前 30% 迭代用eta_c=15,后面增大到30

3.3 主循环设计与收敛条件

主循环要控制的变量包括种群大小、最大迭代数、精英个数和终止条件。36 单元的问题设计变量 72 维,种群太小容易早熟,太大会让单次迭代极慢。

popsize = 120; % 种群大小 maxgen = 300; % 最大迭代代数 elite_count = 4; % 每代保留的精英个数 dim = 72; lb = [zeros(1,36), -pi*ones(1,36)]; ub = [ones(1,36), pi*ones(1,36)]; [pop, fitness] = init_population(popsize, dim, lb, ub); best_fitness_history = zeros(maxgen, 1); for gen = 1:maxgen % 锦标赛选择 selected = zeros(popsize, dim); for i = 1:popsize idx1 = randi([1, popsize]); idx2 = randi([1, popsize]); if fitness(idx1) <= fitness(idx2) selected(i, :) = pop(idx1, :); else selected(i, :) = pop(idx2, :); end end % 交叉生成子代 offspring = zeros(popsize, dim); for i = 1:2:popsize [offspring(i, :), offspring(i+1, :)] = sbx_crossover(... selected(i, :), selected(i+1, :), lb, ub, 20); end % 变异 for i = 1:popsize offspring(i, :) = polynomial_mutation(offspring(i, :), lb, ub, 20, 0.1); end % 子代适应度 offspring_fitness = zeros(popsize, 1); for i = 1:popsize offspring_fitness(i) = fitness_function(offspring(i, :)); end % 精英保留:合并父代和子代 combined_pop = [pop; offspring]; combined_fitness = [fitness; offspring_fitness]; [sorted_fitness, sort_idx] = sort(combined_fitness); pop = combined_pop(sort_idx(1:popsize), :); fitness = sorted_fitness(1:popsize); best_fitness_history(gen) = fitness(1); % 早停条件:连续 50 代适应度变化小于阈值 if gen > 50 && abs(best_fitness_history(gen) - best_fitness_history(gen-20)) < 0.05 fprintf('收敛于第 %d 代\n', gen); break; end end

逻辑说明:精英保留策略保证最优个体不会因为交叉变异而丢失,这是 GA 收敛性的基本保障。早停条件是 50 代的窗口对比,阈值 0.05 dB 根据适应度函数的数值范围设定,如果你的 MSLL 优化的典型跨度是 -10 到 -30 dB,0.05 是合理的精度。

参数说明:这里pop=120不是拍脑袋定的。72 维新问题常常需要变量数的 1.5 到 2 倍种群规模才能保证初期探索的多样性,即 108 到 144。600 代以内的典型收敛代数也不固定,当你看到 50 次迭代内最优适应度不再下降时,通常说明算法已经陷入了局部极值,此时要调的是变异概率而不是代数。

4. 结果分析与参数敏感性:GA 调优的几个正确姿势

4.1 一看收敛曲线,二看方向图,三看动态范围

优化跑完后最重要的验证步骤是画出收敛曲线和最终方向图。收敛曲线能判断早停条件是否触发、是否还有下降空间;方向图能直观看到旁瓣分布是否符合预期;激励幅度分布则暴露了一个非常实际的问题,阵元的动态范围。

figure; plot(best_fitness_history, 'LineWidth', 1.5); xlabel('迭代代数'); ylabel('适应度值'); title('遗传算法收敛曲线'); grid on;
% 绘制优化后的方向图 best_x = pop(1, :); I_opt = best_x(1:36) / max(best_x(1:36)); beta_opt = best_x(37:72); theta_plot = linspace(0, 180, 361) * pi/180; phi_plot = linspace(0, 360, 721) * pi/180; AF_opt = circular_array_factor(I_opt, beta_opt, theta_plot, phi_plot, 0.8, 1); AF_dB = 20*log10(abs(AF_opt)/max(abs(AF_opt(:))) + eps); figure; [H, P] = meshgrid(theta_plot*180/pi, phi_plot*180/pi); surf(H, P, AF_dB, 'EdgeColor', 'none'); xlabel('theta (deg)'); ylabel('phi (deg)'); zlabel('增益 (dB)'); colorbar; view(3);

注意一个常见陷阱:只看最大旁瓣电平的改善是不够的。因为单一直方图可能掩盖方向图中其他角度的旁瓣峰值略低于设定值但面积很大的情况。实际工程中还要看平均旁瓣电平和方位的对称性,特别是圆环阵列存在栅瓣簇,某些方位角的旁瓣可能比主瓣区域的峰值高出不少。

4.2 交叉率和变异率的敏感区间

遗传算法的超参数调优通常比问题本身的数学性质更让初学者头疼。对于一个 72 维的实数编码问题,参数可调的范围大致如下:

参数推荐范围数值偏小的影响数值偏大的影响
种群大小100~150早熟,种群多样性不够单代计算成本线性上升,过 200 对 72 维问题意义不大
交叉分布指数 eta_c15~25子代偏离父代过远,破坏收敛倾向子代过似父代,搜索变慢
变异概率 pm0.05~0.2陷入局部最优难以跳出变成随机搜索,最优解无法稳定保留
变异分布指数 eta_m10~30扰动太大,精英附近难以精调扰动太小,对局部最优的逃离能力不足

一个容易忽略的点是变异概率应当随迭代进程变化。收敛后期,种群高度同质化,此时维持一个高变异概率(比如 0.2)有助于跳出局部极值;相反,如果前期就设 0.2,算法会变成随机搜索,最优个体的适应度波动很大。常见做法是把变异概率从 0.15 到 0.05 线性衰减,或用自适应策略按种群多样度动态调整。

4.3 幅度动态范围的现实约束与惩罚函数修正

仿真里的最优激励分布如果出现极端值,例如某些阵元幅度不足 0.1,方向图的旁瓣水平确实可以压得很低,但工程上这意味着那部分阵元几乎被关断,对硬件误差极其敏感,馈电网络的设计难度也陡增。解决思路是在适应度函数中增加一个动态范围惩罚项:

function fitness = fitness_with_dr(I, beta) % 在原有 fitness 基础上增加动态范围惩罚 base_fitness = fitness_function([I, beta]); I_normalized = I / max(I); dynamic_range = 20 * log10(max(I_normalized) / min(I_normalized) + eps); % 动态范围超过 15 dB 时开始惩罚 dr_threshold = 15; % dB lambda_dr = 0.2; dr_penalty = lambda_dr * max(0, dynamic_range - dr_threshold); fitness = base_fitness + dr_penalty; end

4.4 多轮优化与随机种子管理

遗传算法是随机算法,单次运行的最优结果不一定是全局最优。成熟做法是固定随机种子,跑 3 到 5 次独立优化,每次都从不同初始种群出发。如果多次结果的方向图结构相似,可认为结果可信;如果差异很大,说明搜索空间不够充分,需要增大种群或迭代代数。MATLAB 里设置rng(2024)即可复现某次优化结果。调参时用一个固定种子,验证鲁棒性时换多个种子,这是识别“偶然最优”和“结构最优”最直接的手段。

注意:同一种子下微调参数产生的方向图差异,可能只是初始种群的偶然性变化,不代表参数优劣。正确做法是每个参数组合用至少 3 个种子各跑一遍,取最优值的均值做对比。

5. 从 36 单元走向工程:耦合修正、扫描角鲁棒性与混合优化

上述流程解决的是“理想阵列因子”层面的优化,即所有阵元都是理想点源、互耦为零、通道幅度相位完全精确。但实际相控阵天线设计中,这三点都不成立。单元间的互耦会显著改变阵列的有源方向图,尤其对于圆环阵列,单元间的方位关系决定了互耦矩阵不是简单的平移不变矩阵,用理想点源算出的最优激励分布上机测试后旁瓣水平通常会恶化 3 dB 以上。常见做法是分两步走:先用 GA 在理想模型上找到合理的设计空间,再用全波仿真软件(如 HFSS、CST)提取有源单元方向图数据,代入 MATLAB 重建方向图,最后用 GA 做一次本地的精细校正。

扫描角鲁棒性也是一个常被忽略的工程需求。如果阵列需要在 ±45° 范围内扫描,每个扫描角都有不同的最优激励分布。把扫描角直接作为适应度函数中的附加维度,可以在一次优化中得到兼顾多个扫描角的方向图。实现时不需要修改算法的核心循环,只需要在计算适应度时计算多个角度的方向图并加权求和。

另一个实用的进阶技巧是 GA 与局部搜索算法的混合。GA 擅长在大空间中找到好的区域,但后期收敛速度远慢于梯度类算法。一个直接的做法是用fmincon把 GA 得到的最优解做局部精修,约束条件和代码结构几乎不用改动,只需把 GA 最优解作为初始点传入:

options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp'); x_refined = fmincon(@(x) fitness_function(x), best_x, ... [], [], [], [], lb, ub, @constraint_function, options);

约束函数constraint_function可以写成统一的数组输出形式,适合用来强约束主瓣波纹或指定频率点。混合策略的实际收益一般在 1 到 2 dB 的旁瓣改善,但要注意混合优化会把变量推向设计空间的边界,幅度参数可能出现大量接近 0 的值,需要要配合第 4 章的动态范围惩罚一起使用。

最后说一个验证 GA 优化结果有效性的小技巧:把幅度分布强制设为 1(满激励),只优化相位,再和幅度相位同时优化的结果对比。如果两种策略的方向图差异不大,说明你的阵列模型和约束条件对幅度的依赖度很低;如果幅度优化带来了显著旁瓣改善,那就要检查是不是旁瓣惩罚项太容易被幅度调节所利用了。这一步能有效识别过度优化,也是评估阵列能否用更简单的移相器实现的关键依据,方向图综合的工程验证,到这一步才算真正闭环。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询