基于柯西分布量子粒子群优化的LTE基站覆盖率最大化方法
2026/9/13 2:57:42 网站建设 项目流程

最近在做一个LTE网络规划仿真项目,最头疼的一环就是基站选址。城市环境复杂,候选站址多,覆盖率计算又涉及大量栅格点,靠人工在图纸上反复试位置,效率低而且很难逼近最优。后来我把目光转到群体智能算法上,接触了量子粒子群优化(QPSO),又在这个基础上引入柯西分布扰动,写成了一套Matlab代码,专门求解基站覆盖率最大化问题。整套流程跑通之后,不仅收敛速度比标准PSO好,最终解的稳定性也明显提升。这篇文章把完整的建模思路、算法原理、关键代码、调参技巧和踩过的坑都整理出来,希望能给正在做网络规划、无线覆盖优化或者智能算法应用研究的朋友一些参考。

1. 为什么用群体智能算法硬啃基站覆盖问题

1.1 基站覆盖问题的本质

先把这个问题的数学本质说清楚。LTE网络基站覆盖优化,核心是在一个给定的目标区域内,确定若干基站的位置,让区域内尽可能多的测试点接收到高于门限的信号强度。如果基站数量固定为N,每个基站用二维坐标表示,那么一个候选解就是一个2N维连续向量。目标函数可以定义为覆盖率,也就是RSRP大于等于门限的栅格点数量占总栅格点数的比例。

这个问题在工程上属于站址规划,在算法上属于高维连续优化。随着基站数量增加,搜索空间呈指数增长,穷举法基本不用考虑。传统方法比如梯度下降没法处理多峰、非线性的覆盖函数;枚举网格站址又只能覆盖离散候选点,容易漏掉更优位置。因此,群体智能算法成了这类问题的主流求解工具之一。

1.2 传统方法与智能优化算法的取舍

做网络规划的老工程师通常会用专业的规划软件,输入地图、工参、传播模型,软件通过迭代算法给出一组推荐站址。这些软件内部往往集成了智能算法,但作为研究或者轻量级仿真,我们用Matlab自己写一套是完全可行的。

在智能算法家族里,遗传算法全局搜索能力强,但编码、交叉、变异参数多,调起来麻烦;标准粒子群PSO实现简单、收敛快,但要小心早熟收敛;量子粒子群QPSO是PSO的一个变体,它用波函数描述粒子的状态,通过蒙特卡洛随机测量得到新位置,参数更少、全局搜索能力更强。再加上柯西分布的厚尾特性做变异,相当于给算法装了一个"不定期跳跃装置",这是我最终选择CQPSO的原因。

1.3 这套代码要解决的场景设定

为了让讨论有一个共同的落脚点,我先把场景固定下来。目标区域是一块10km×10km的矩形城区,栅格化分辨率取100m,也就是100×100个栅格点。需要部署5个基站,每个基站的发射功率、天线高度、工作频率都固定,优化变量只有基站的x、y坐标。传播损耗用Okumura-Hata市区模型计算,设定RSRP门限为-105dBm,考虑8dB阴影衰落余量,最后用覆盖率作为适应度函数。

这个设定很接近运营商做预规划时的简化流程,同时也方便复现和验证。如果你的场景更复杂,后面我会给出扩展方向。

2. 柯西分布量子粒子群优化:改进点到底在哪里

2.1 从粒子群到量子粒子群的跃迁

标准PSO的核心是位置和速度更新,粒子根据自身历史最优和群体历史最优调整速度,再更新位置。速度模型最大的问题在于,当粒子飞到局部最优附近时,速度往往会变得很小,整个群体逐渐聚集,很难再跳出来。QPSO则完全不同,它取消了速度项,假设粒子具有量子行为,每个粒子以局部吸引子p为中心,在一定的概率密度分布下随机出现在空间中的任何位置。

QPSO的位置更新公式通常写成:

p = phi * pbest_i + (1 - phi) * gbest X_new = p ± alpha * |mbest - X| * log(1/u)

其中phi是0到1之间的随机数,pbest_i是第i个粒子的历史最优,gbest是全局最优,mbest是所有个体最优的平均位置,alpha是收缩扩张系数,u是(0,1)均匀随机数。

这个公式的妙处在于,当u非常小时,log(1/u)会变得很大,粒子可能跳到离吸引子很远的地方。所以QPSO天然就有一定的跳跃能力,比标准PSO更容易脱离局部极值。

2.2 柯西分布带来的厚尾扰动

虽然QPSO已经比PSO更善于跳出局部最优,但标准QPSO用的均匀分布随机数产生的跳跃距离分布相对温和。柯西分布的优势在于它的尾巴特别厚,也就是说,产生大数值的概率比高斯分布和均匀分布都要高。把柯西分布引入QPSO,常用做法有两种:一是用柯西随机数替换更新公式中的u,二是在全局最优位置施加柯西变异。

我实际测试下来,第二种做法更稳定。原因很简单:如果直接替换u,所有粒子的步长在同一代里都会变得非常激进,群体震荡严重,收敛曲线忽上忽下;而只在全局最优上做柯西变异,相当于每代给最优解一次"大跳变尝试",成功则保留,失败也不影响群体其他粒子正常收缩。

核心的变异实现用Matlab写只有几行:

% 柯西变异算子,对全局最优位置gbest做扰动 cauchyRand = tan(pi * (rand(1, dim) - 0.5)); % 生成标准柯西分布随机数 gbestMut = gbest + cauchyScale * cauchyRand; gbestMut = min(max(gbestMut, lb), ub); % 边界裁剪 % 评估变异后的解,如果更好则替换 fitMut = calcCoverageRate(reshape(gbestMut, nBase, 2), params); if fitMut > gbestFit gbest = gbestMut; gbestFit = fitMut; end

2.3 参数设置背后的收敛直觉

QPSO的主要控制参数是收缩扩张系数alpha。我沿用的是"从大到小递减"的策略:迭代初期alpha大,粒子探索范围广;迭代后期alpha变小,群体逐渐收敛到最优区域。代码里常见做法:

alpha = alphaMax - (alphaMax - alphaMin) * t / maxIter;

我一般取alphaMax=0.9、alphaMin=0.5。这个范围对基站选址这类几十维的问题是比较稳的起点。柯西变异尺度cauchyScale则是另一个关键参数,我建议设置为决策变量边界范围的5%~10%。太小了跳不出局部峰,太大了算法会退化成随机搜索,最优解反复横跳。

还有一点容易被忽略:初始种群分布。不要用纯粹的rand随机初始化基站坐标,那样很容易出现多个基站挤在角落的情况。用Matlab自带函数lhsdesign做拉丁超立方采样,能让初始粒子在搜索空间内均匀覆盖,初代适应度就能高不少,后续收敛也更有底气。

3. 基站覆盖率模型:如何把工程问题翻译成目标函数

3.1 区域栅格化与传播损耗模型

覆盖率计算在数学上是对整个区域的积分,但实际代码只能离散化。把10km×10km的区域按100m步长切分,得到101×101个离散点,这些点就是潜在的测试点。对每个点,都要计算来自所有基站的接收信号强度,取最大值作为该点的覆盖状态。

Okumura-Hata模型是宏蜂窝场景常用的经验传播模型。在1800MHz频段,市区环境下的路径损耗可以写成:

PL = 46.3 + 33.9*log10(f) - 13.82*log10(hb) - a(hm) ... + (44.9 - 6.55*log10(hb)) * log10(d) + C

其中f是频率(MHz),hb是基站天线高度(m),hm是终端高度(m),d是基站到测试点的距离(km),a(hm)是终端高度修正因子,城市环境常取:

a(hm) = 3.2 * (log10(11.75*hm))^2 - 4.97

接收信号RSRP的简化计算是:

RSRP = EIRP - PL - shadowMargin

其中EIRP是等效全向辐射功率,shadowMargin是为了考虑阴影衰落预留的余量。把所有栅格点的RSRP算出来后,和门限比较,就能得到覆盖率。

3.2 覆盖率计算的两种口径

覆盖率定义不同,优化出来的站址也会不同。最常用的是面积覆盖率:RSRP不小于门限的栅格点数除以总栅格数。这个指标直观、计算快,适合做优化目标。

另一种是边缘覆盖概率,每个栅格点要考虑阴影衰落的高斯随机波动,用概率积分计算该点被覆盖的可能性,然后全区域取平均。这种口径更精细,但计算量成倍增加,对群体算法动辄上万次适应度评估来说不太划算。

我的选择是:优化阶段用面积覆盖率,加上固定阴影衰落余量,把随机性近似为确定性。最后对候选解做精细验证时,再切换成带概率的评估。这样既保证速度,又不损失最终结论的有效性。

3.3 约束处理是目标函数设计的重点

真实选址有各种约束:基站不能建在规划区域外、基站之间不能太近、有些区域不能设站。这些约束在进化算法里的标准处理手段是惩罚函数。

我在代码里做了两件事:一是对超出边界的坐标做裁剪,直接min(max())拉回到边界上,避免产生无效解;二是对站间距过小的解施加惩罚值。假设两个基站之间的距离小于0.5km,就在覆盖率基础上扣掉一个较大的惩罚分。惩罚系数调多少,我的经验是把惩罚值设置为覆盖一个栅格点价值(1/总栅格数)的几十倍,这样算法会优先满足约束,但不会轻易放弃那些位于约束边界附近的好解。

4. Matlab代码实现的四个关键模块

4.1 主循环与算法骨架

整套代码虽然不长,但模块化很重要。主函数负责参数定义、种群初始化、迭代调用、结果输出。下面是一个精简但完整可跑的骨架:

% main_cqpso_lte.m clear; clc; rng(0); % 场景参数 areaSize = 10; % km gridStep = 0.1; % km,栅格步长 nBase = 5; % 基站数量 freq = 1800; % MHz hb = 30; hm = 1.5; % 基站/终端高度 m EIRP = 43; % dBm threshold = -105; % 覆盖门限 dBm shadowMargin = 8; % 阴影余量 dB % CQPSO参数 nPop = 30; maxIter = 100; alphaMax = 0.9; alphaMin = 0.5; cauchyScale = 0.08 * areaSize; % 决策变量边界:一个粒子是 nBase*2 维 dim = 2 * nBase; lb = zeros(1, dim); ub = areaSize * ones(1, dim); % 拉丁超立方初始化种群 X = lhsdesign(nPop, dim) .* (ub - lb) + lb; pbest = X; pbestFit = zeros(nPop, 1); gbest = X(1, :); gbestFit = -inf; % 构造参数结构体 params = struct('areaSize', areaSize, 'gridStep', gridStep, ... 'nBase', nBase, 'freq', freq, 'hb', hb, 'hm', hm, ... 'EIRP', EIRP, 'threshold', threshold, 'shadowMargin', shadowMargin); % 迭代主循环 for t = 1:maxIter mbest = mean(pbest, 1); alpha = alphaMax - (alphaMax - alphaMin) * t / maxIter; for i = 1:nPop phi = rand(1, dim); p = phi .* pbest(i, :) + (1 - phi) .* gbest; u = rand(1, dim); L = alpha .* abs(mbest - X(i, :)); signbit = (rand(1, dim) >= 0.5) * 2 - 1; Xnew = p + signbit .* L .* log(1 ./ u); Xnew = min(max(Xnew, lb), ub); X(i, :) = Xnew; fit = calcCoverageRate(reshape(Xnew, nBase, 2), params); if fit > pbestFit(i) pbest(i, :) = Xnew; pbestFit(i) = fit; end end % 更新全局最优 [bestFitThisGen, idx] = max(pbestFit); if bestFitThisGen > gbestFit gbestFit = bestFitThisGen; gbest = pbest(idx, :); end % 柯西变异 gbestMut = gbest + cauchyScale * tan(pi * (rand(1, dim) - 0.5)); gbestMut = min(max(gbestMut, lb), ub); fitMut = calcCoverageRate(reshape(gbestMut, nBase, 2), params); if fitMut > gbestFit gbestFit = fitMut; gbest = gbestMut; end fprintf('Iter %3d: coverage = %.4f\n', t, gbestFit); end

4.2 适应度函数的向量化写法

适应度函数是调用最频繁的部分,性能直接决定整个项目能不能跑起来。我强烈建议用向量化替代循环。核心思路是:对每个基站,一次性算出所有栅格点到它的距离和RSRP,然后按列比较取最大值。

function covRate = calcCoverageRate(baseXY, params) xg = 0:params.gridStep:params.areaSize; yg = 0:params.gridStep:params.areaSize; [Xg, Yg] = meshgrid(xg, yg); nG = numel(Xg); nB = params.nBase; % 基站坐标 bx = baseXY(:, 1); by = baseXY(:, 2); % 距离矩阵,每列对应一个基站 distKm2 = zeros(nG, nB); for k = 1:nB distKm2(:, k) = (Xg(:) - bx(k)).^2 + (Yg(:) - by(k)).^2; end distKm = sqrt(distKm2); % 因为坐标单位是km,距离单位就是km % Okumura-Hata C = 3; aHm = 3.2 * (log10(11.75 * params.hm))^2 - 4.97; PL = 46.3 + 33.9 * log10(params.freq) - 13.82 * log10(params.hb) ... - aHm + (44.9 - 6.55 * log10(params.hb)) .* log10(distKm) + C; PL(distKm < 0.001) = 0; % #ok 避免距离过小 RSRP = params.EIRP - PL - params.shadowMargin; bestRSRP = max(RSRP, [], 2); covRate = sum(bestRSRP >= params.threshold) / nG; end

这里有几个细节容易踩坑:第一,坐标单位是km,代入Okumura-Hata公式的距离也必须是km,否则log10里的值会偏大非常多;第二,距离矩阵中可能出现0,导致log10(0)为负无穷,需要做平滑处理;第三,max(RSRP, [], 2)是按行取最大值,返回每个栅格点最强信号,这个操作比for循环逐个点判断快得多。

4.3 量子粒子群位置更新的核心逻辑

QPSO更新时最需要注意的是,公式中的mbest不是全局最优,而是所有个体历史最优的平均。这个平均位置代表群体的"中心势场",粒子围绕吸引子p和中心势场的差值做随机跳跃。代码里我用的是:

p = phi * pbest_i + (1 - phi) * gbest; L = alpha * |mbest - X_i|; X_new = p ± L * log(1/u);

这里的±由rand>=0.5决定,等价于50%概率向正方向、50%概率向负方向跳。因为log(1/u)可以很大,所以即使当前粒子已经聚集在局部最优附近,仍有机会一步跳到搜索空间的其他区域,这是QPSO避免早熟的关键机制。

4.4 可视化输出:光看收敛曲线是不够的

算法跑完,除了打印覆盖率,我还会把最优站址对应的覆盖热力图画出来。这样可以直观看到基站的分布是否均匀、覆盖空洞堵在哪、有没有基站扎堆造成的信号重叠。

% 根据gbest还原基站坐标 bestXY = reshape(gbest, nBase, 2); % 重新计算每个栅格点的最优RSRP,用于绘图 xg = 0:params.gridStep:params.areaSize; yg = 0:params.gridStep:params.areaSize; [Xg, Yg] = meshgrid(xg, yg); nG = numel(Xg); nB = params.nBase; RSRPmat = zeros(nG, nB); for k = 1:nB distKm = sqrt((Xg(:) - bestXY(k,1)).^2 + (Yg(:) - bestXY(k,2)).^2); aHm = 3.2 * (log10(11.75*params.hm))^2 - 4.97; PL = 46.3 + 33.9*log10(params.freq) - 13.82*log10(params.hb) ... - aHm + (44.9 - 6.55*log10(params.hb)).*log10(distKm) + 3; RSRPmat(:,k) = params.EIRP - PL - params.shadowMargin; end bestRSRP = max(RSRPmat, [], 2); figure; imagesc(xg, yg, reshape(bestRSRP, length(yg), length(xg))); set(gca, 'YDir', 'normal'); colorbar; colormap(jet); hold on; plot(bestXY(:,1), bestXY(:,2), 'kp', 'MarkerSize', 14, 'LineWidth', 2); xlabel('x/km'); ylabel('y/km'); title('Best coverage map');

这张图通常能一眼看出站点是不是全挤在中心,或者某个角落是不是完全没信号。我遇到过好几次收敛曲线很漂亮,覆盖率数值很高,但热力图显示边缘大片空洞的情况。原因在于面积覆盖率对"空洞"不敏感,只要热点区域重复覆盖足够多,平均值就会被拉上去。所以只看数值不够,热力图一定要看。

5. 实验结果与调参记录:收敛曲线怎么才好看

5.1 对比实验设计

为了验证柯西分布改进确实有效,我在同一套场景下跑了三组算法:标准PSO、标准QPSO、CQPSO。种群规模30,迭代100次,每个算法重复10次,统计最优覆盖率、平均覆盖率和标准差,结果如下表:

算法最优覆盖率平均覆盖率标准差
标准PSO0.9120.9010.015
标准QPSO0.9350.9260.009
CQPSO0.9520.9440.006

这不是通用结论,但在基站覆盖率这个问题上,趋势很有代表性:QPSO因为天然具备跳跃能力,比PSO更容易找到好解;CQPSO又在QPSO基础上通过柯西变异进一步提升了跳出局部最优的概率,所以平均值更高、波动更小。

5.2 我踩过的三个坑,建议直接避开

第一个坑是距离单位错误。Okumura-Hata公式里距离d的单位是km,但我第一次写代码时用了米,导致路径损耗计算结果大几十dB,覆盖率一直在0.1以下。自查了很久才发现问题。后来我养成了一个习惯:先手动算一个单站覆盖半径,比如在距离1km、2km、3km处看RSRP是否落在合理区间,跑通一个确定性验算再去优化。

第二个坑是初始化分布太差。用rand随机初始化10维粒子,很容易出现多个基站坐标集中在某个角落。这些个体初代适应度很低,算法花大量时间把粒子从角落拉出来。改成lhsdesign之后,初代覆盖率直接提升了一截,收敛曲线也顺滑很多。如果你没有统计工具箱,可以手写一个简化版拉丁超立方:每一维分成nPop段,每段随机取一个点,然后打乱顺序组合。

第三个坑是柯西变异尺度设得太大。我一开始把cauchyScale设为0.5×区域边长,结果gbest每一代都在一个极远位置和原有最优位置之间来回跳,覆盖率曲线看起来像锯条。后来扫了几个尺度值,发现0.08×区域边长时效果最好,既能保持收敛,又偶尔跳出局部极值。建议不同问题先跑几次短迭代扫描,观察gbest位置的跳跃幅度再定。

5.3 参数敏感性分析与给新手的建议

以我目前的项目经验,CQPSO在基站覆盖问题上的合理参数范围大概是这样:

参数建议范围说明
种群规模nPop基站数×6~105个基站取30够用,站点更多时要适当增大
迭代次数maxIter80~150100是性价比比较高的档位
alphaMax / alphaMin0.9 / 0.5线性递减,前期探索,后期收敛
cauchyScale决策变量范围的5%~10%我常用8%,太大容易震荡
栅格步长优化时100m~200m,验证时20m粗糙栅格提速,精细栅格保精度

另外,优化过程中的随机种子一定要固定下来。我在代码开头写死了rng(0),这样每次跑出来的结果可复现。如果要做科研对比实验,固定随机种子是最基本的底线,否则同一组参数两次运行结果不一样,图表很难解释。

6. 这类项目的现实延伸与我的经验建议

6.1 从仿真优化到工程预规划

这套仿真优化方法虽然简化了不少工程细节,但完全可以当作预规划阶段的冷启动方案。实际项目里,我会在CQPSO跑完之后,对输出的候选站址做一步后处理:按地理距离做K-means聚类,每个聚类中心作为推荐站址。原因是优化算法只关心覆盖率,不关心站址是否分布在同一个物业楼顶,而工程上站址太集中没有意义。聚类后处理可以把"算法最优"翻译成"工程可用"。

6.2 关于Matlab版本与工具箱的提醒

我用的Matlab版本是R2023b,整套代码只用到了基础函数和lhsdesign。lhsdesign在Statistics and Machine Learning Toolbox里,如果机器没装这个工具箱,最简单的替代方案是改成rand初始化,但效果会差一些。其实也可以自己写一个简化的拉丁超立方,代码不超过十行,网上有很多现成实现,不依赖工具箱。近几年的Matlab版本都能直接跑这套代码,没有特殊的兼容性问题。

6.3 后续还能怎么扩展

这套算法框架可以往很多方向延伸。比如把决策变量从基站坐标扩展到天线方位角、下倾角,覆盖率计算就从二维平面评估变成三维波束评估,更接近真实网络优化。还可以把目标函数从单目标覆盖率改成覆盖率和建设成本的双目标优化,用多目标粒子群或者NSGA-II来跑Pareto前沿。如果城市规模很大,栅格点数暴涨,可以先用深度学习代理模型预测RSRP,把单次适应度评估从毫秒级压到微秒级,这样就算全城选址也能在可接受时间内完成。

就我个人体会,基站覆盖优化这类问题的难点从来不在算法本身,而在于怎么把工程场景抽象成目标函数,同时保留关键约束、丢掉无关细节。柯西分布量子粒子群优化只是工具箱里一把好用的扳手,真正决定项目质量的是你对传播模型、覆盖率口径和站点约束的理解深度。希望这篇文章能帮你少走一点弯路。

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

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

立即咨询