简介面向无线传感器网络学习者的 LEACH 路由协议 MATLAB 实现代码适合做课程实验、毕业设计或节能路由算法入门复现。LEACH 通过随机簇头选举、簇内数据聚合与轮换机制均衡节点负载常用于观察网络寿命、能量消耗与数据传输效率。压缩包仅 2KB内含 1 个 .m 文件属于精简的 MATLAB 源码包便于直接查看、运行和二次修改。已有 689 人学习。代码围绕节点随机分布、簇头概率选举、簇内与基站通信、能量消耗更新及多轮仿真展开可帮助读者理解协议主流程并在此基础上调整簇头比例、初始能量和网络规模记录并绘制存活节点数、剩余能量等性能曲线。单文件结构也便于逐段调试、参数对照和结果复现为算法对比、课程报告与论文实验提供可运行的基础脚本。1. 先看结论一份 LEACH 路由协议的 MATLAB 代码能省掉你两周的踩坑时间接手分簇路由这类课题时最耗时间的通常不是公式推导而是“把公式变成一条能跑出曲线的能量曲线”。这份 LEACH 路由协议的 MATLAB 代码解决的就是这件事它把簇头选举、成簇、TDMA 时隙分配、稳态传输和一阶无线电能量模型串成一条可执行的链路跑完一轮就能拿到节点剩余能量、簇头分布和死亡节点数。适合三类人要快速复现 LEACH 基线做对比实验的研究生、要在自家协议前先搭一个可信参照系的工程师、以及被“代码跑得通但结果明显不对”折腾过的人。它不解决协议创新但把基线做扎实后面的对比才有意义。2. 一阶无线电能量模型把“耗电”换算成可对比的焦耳数LEACH 之所以能画出那条经典的“轮次—存活节点”曲线核心不在选举公式而在能量模型足够诚实。换一套能量参数同一份协议代码能给出完全相反的结论所以动手写选举之前先把这个模型固定下来。2.1 为什么用一阶无线电模型而不是简单的 d² 衰减一阶无线电模型First Order Radio Model把发送端能耗拆成两部分射频电路消耗和功放消耗。功放部分按距离分成自由空间模型和多径衰落模型分界点是阈值 d0。近距离下信号按 d² 衰减远距离下按 d⁴ 衰减。这个分界不是装饰。它直接决定了 LEACH 的一个关键行为簇头到基站距离远时簇头必须用 d⁴ 发送汇聚数据这是簇头比普通节点死得快的物理原因也是后续所有“能量均衡”改进LEACH-C、SEP、DEEC的出发点。如果你图省事全程用 d² 算发送能耗会发现簇头寿命异常长、曲线漂亮得不像真的然后你的改进算法在仿真里“提升有限”。这不是算法问题是模型被你简化掉了关键非线性项。常见参数取值如下这份代码默认用的就是这一套参数符号典型取值说明电路能耗E_elec50 nJ/bit收发各计一次自由空间功放系数E_fs10 pJ/bit/m²d d0 时生效多径功放系数E_mp0.0013 pJ/bit/m⁴d ≥ d0 时生效分界距离d0√(E_fs/E_mp) ≈ 87.7 m由前两项算出来不要手填数据聚合能耗E_DA5 nJ/bit/signal簇头每聚合一个包加一次初始能量E00.5 J全网统一数据包长k_data4000 bit稳态阶段每个成员上行控制包长k_ctrl200 bit广播、入簇请求d0 一定要用公式算不要抄一个 87.7 进去。一旦你改 E_fs 或 E_mp 做敏感性分析手填的 d0 会和公式脱节出现“远距离反而比近距离省电”的诡异结果——这类翻车我在复现别人的对比实验时见过不止一次。2.2 发送、接收、聚合三个函数的落地写法把能量计算拆成三个独立函数是整个工程里最值得的一个决定。后面做多跳、做中继、做异构节点都只是换个地方调用它们不用动主循环。function E tx_energy(k, d, prm) % TX_ENERGY 发送 k bit 到距离 d 处的能耗J % k : 数据包长度单位 bit % d : 发送端到接收端的欧氏距离单位 m % prm : 参数结构体含 E_elec / E_fs / E_mp / d0 if d prm.d0 % 自由空间模型近距离功放按 d^2 增长 E k * prm.E_elec k * prm.E_fs * d^2; else % 多径衰落模型远距离功放按 d^4 增长 E k * prm.E_elec k * prm.E_mp * d^4; end end function E rx_energy(k, prm) % RX_ENERGY 接收 k bit 的能耗J与距离无关 E k * prm.E_elec; end function E da_energy(k, prm) % DA_ENERGY 簇头聚合 k bit 数据的额外能耗J E k * prm.E_DA; end逻辑说明发送能耗的第一项是电路开销与距离无关只要发就得付第二项是功放开销随距离呈平方或四次方增长。接收端只付电路开销因为接收不需要功放。聚合开销只发生在簇头且和“聚合了几个包”成正比不是和比特数成正比——严格写应该是E_DA * 成员数这里按比特折算是一种近似误差小但要在论文里说清楚。参数说明prm结构体统一管理所有物理量主脚本里一次性初始化后面所有函数只传引用式的结构体避免十几个全局变量满天飞。prm.E_elec 50e-9; % 50 nJ/bit prm.E_fs 10e-12; % 10 pJ/bit/m^2 prm.E_mp 0.0013e-12; % 0.0013 pJ/bit/m^4 prm.E_DA 5e-9; % 5 nJ/bit/signal prm.d0 sqrt(prm.E_fs / prm.E_mp); % 约 87.7 m由公式导出 prm.E0 0.5; % 初始能量 J prm.k_data 4000; % 数据包 bit prm.k_ctrl 200; % 控制包 bit2.3 用单点校验确认模型没写错模型写完别急着跑全流程先做一组手算对照。取 k4000 bit、d50 m自由空间项4000 × 50e-9 2.0e-4 J功放项4000 × 10e-12 × 2500 1.0e-4 J合计 3.0e-4 J。若 d100 m超过 d0功放项变成 4000 × 0.0013e-12 × 1e8 5.2e-4 J合计 7.2e-4 J。距离翻倍能耗从 3.0e-4 涨到 7.2e-4翻了 2.4 倍而不止 4 倍——这就是模型跨界切换的直接证据。写个五行的脚本对一遍比跑完 1000 轮发现曲线不像 LEACH 再去回溯要省事得多。% 手算对照验证 d0 两侧的能耗跳变 for d [50, 87, 88, 100] fprintf(d%5.1f m, E_tx%.4e J\n, d, tx_energy(4000, d, prm)); end输出的四个值应当呈现“88 m 处斜率明显变陡”的形态。看不到这个拐点说明prm.d0被覆盖过或者tx_energy里的判断写成了之外的别的比较。这个自检花两分钟能拦住后面一连串说不清的现象。3. 簇头选举把阈值公式落地成不会选出死节点的代码能量模型定死了接下来是 LEACH 的灵魂——每轮按概率 p 选簇头。公式在纸上只有一行落到代码里有三处容易出错G 集合的维护方式、轮次的取模边界、以及“能量耗尽节点必须排除”这条硬约束。3.1 阈值公式里每个符号的真实含义标准写法是节点 n 在第 r 轮被选为簇头的概率为T(n) p / (1 − p × (r mod 1/p))当 n ∈ G否则 T(n) 0。p 是期望簇头比例比如 0.05 表示每轮大约 5% 的节点当簇头。G 是“最近 1/p 轮内还没当过簇头”的节点集合。分母那一项的意义是补偿随着轮次推进G 集合里剩下的节点越来越少分母变小单节点的入选概率被抬高从而保证每 1/p 轮每个节点恰好当过一次簇头。很多复现代码把 G 集合写成mod(r, round(1/p))的派生判断能跑但有两个隐患一是当节点因能量耗尽退出网络时它仍被算进 G 集合的分母导致剩余节点被系统性低估概率二是浮点轮次或非整数 1/p 时取模结果不稳。稳妥做法是给每个节点挂一个lastCH记录“上次当簇头的轮次”用轮次差判断是否在冷却期内。3.2 带能量约束的选举实现function chFlag select_ch(E, lastCH, p, r) % SELECT_CH 按 LEACH 阈值公式选出本轮簇头 % E : 1×n 节点剩余能量J % lastCH : 1×n 每个节点上次当选簇头的轮次初始为 -inf % p : 期望簇头比例 % r : 当前轮次 % chFlag : 1×n 逻辑数组true 表示本轮为簇头 n numel(E); chFlag false(1, n); cool round(1 / p); % 冷却窗口长度p0.05 时为 20 轮 for i 1:n if E(i) 0 continue; % 能量耗尽永久退出选举 end % 冷却期内的节点不进 G 集合 if (r - lastCH(i)) cool continue; end T p / (1 - p * mod(r, cool)); if rand() T chFlag(i) true; end end end逻辑说明三层过滤依次生效——死节点排除、冷却期排除、剩余节点按阈值随机判定。冷却窗口cool用round(1/p)而不是1/p本身是为了兼容 p0.07、0.1 这类不能被 1 整除的取值避免mod作用在浮点上产生1e-16量级的残差把判断带偏。参数说明p不要超过 0.1。p 太大簇头过多聚合开销和簇头到基站的远距离传输会同时上升整网能耗反而变差p 太小则单簇成员过多簇头聚合压力大、TDMA 帧变长。0.05 是文献里最常用的起点做 p 的敏感性扫描时建议取 0.03 / 0.05 / 0.08 / 0.1 四档。还有一点值得单独强调选举完必须做一次兜底检查。if ~any(chFlag) any(E 0) % 极端情况本轮一个簇头都没选出来概率低但会发生 alive find(E 0); chFlag(alive(randi(numel(alive)))) true; end没有这段兜底小规模网络比如 20 个节点、p0.05偶尔会出现“全轮无簇头”此时所有节点保持睡眠状态、本轮零能耗曲线上表现为一个不该存在的平段。这种偶发 bug 最难查因为重跑一次就没了。3.3 轮次推进与“每轮重置什么”主循环里最容易搞混的是“哪些状态跨轮保留、哪些每轮清零”。给个明确的分工变量是否跨轮保留原因节点坐标保留静态网络部署后不变剩余能量 E保留逐轮递减lastCH保留冷却期判定的依据簇头标记 chFlag每轮重置只对本轮有效簇成员关系每轮重置成簇结果每轮重算TDMA 时隙表每轮重置随簇成员数变化死亡节点计数累计用于画存活曲线这张表看似琐碎但它对应的是无数人栽过的坑把chFlag定义在循环外第二轮开始所有节点都还是簇头曲线直接躺平。我一般会把网络状态封装成一个 struct主循环里只做“复制 → 修改 → 存回”三步减少这种变量作用域错误。net struct(pos, pos, E, E0 * ones(1, n), ... lastCH, -inf(1, n), dead, zeros(1, n)); for r 1:rounds state net; % 本轮工作副本 state.chFlag select_ch(state.E, state.lastCH, p, r); state.lastCH(state.chFlag) r; % 当选的更新冷却计时 [state, stat] run_round(state, prm, bs); net state; % 写回E 已扣减 record(r, :) stat; endlastCH初始化为-inf而不是 0是为了让第 1 轮所有节点都满足r - lastCH cool直接进入候选池。初始化为 0 时第 1 轮会算出1 - 0 1 20全网没人有资格当选然后靠兜底逻辑强选一个——能跑但不是 LEACH 的行为。4. 成簇、TDMA 时隙与稳态传输一次完整轮次的执行链路选举只决定“谁当簇头”真正的能量消耗发生在成簇和稳态传输里。这一章把一轮拆成四个阶段簇头广播、成员入簇、时隙分配、数据上传每个阶段都对应不同的包长和能量扣减。4.1 簇头广播与成员入簇距离矩阵怎么算才不吃亏节点数在 100 左右时直接建距离矩阵比在双重循环里反复调norm快得多而且代码可读性好function D dist_matrix(pos, bs) % DIST_MATRIX 返回 (n1)×(n1) 距离矩阵最后一行为基站 P [pos, bs]; % 把基站拼在末尾 D squareform(pdist(P)); % pdist 要求按列排注意转置 endpdist要求输入是“每列一个观测”所以P是必须的。少了这个转置矩阵会变成把每个坐标分量当成一个节点节点数直接翻倍得到的距离矩阵维度对不上后面索引全乱。这类错误不会报错只会让结果莫名其妙。簇头广播阶段每个簇头发一个控制包扣tx_energy(k_ctrl, d)非簇头收到所有簇头广播后选择距离最近的那个入簇扣若干次rx_energy(k_ctrl)。严格说非簇头应该先收到全部广播再决定所以接收次数等于簇头数chList find(state.chFlag); memberOf zeros(1, n); % 0 表示自己是簇头 for i find(~state.chFlag state.E 0) [~, idx] min(D(i, chList)); memberOf(i) chList(idx); state.E(i) state.E(i) - numel(chList) * rx_energy(prm.k_ctrl, prm); end成员入簇时向簇头发一个控制包扣一次发送能耗。这一步常被省略因为占比小但如果你要做“控制开销占比”这类分析漏掉它就是错的。我的习惯是把所有控制包都算进去宁可多算不要事后被问“你这广播开销去哪了”。4.2 TDMA 时隙分配与帧长计算LEACH 的稳态阶段用 TDMA簇头给自己簇内每个成员分配一个时隙成员只能在自己的时隙发数据其余时间休眠。时隙表本身不耗能计算开销忽略但它决定了稳态阶段的时长进而影响“每轮时间”这个指标。function slot build_tdma(memberOf, chList) % BUILD_TDMA 为每个簇头生成成员时隙表 % 返回 cell 数组第 k 个元素是第 k 个簇头的成员索引向量 slot cell(1, numel(chList)); for k 1:numel(chList) slot{k} find(memberOf chList(k)); end end稳态阶段每个成员在各自时隙里发k_data比特给簇头簇头接收并聚合。簇头自己也要把聚合后的数据发给基站——注意这一步的距离是“簇头到基站”不是簇内距离通常是整个网络里最长的一跳d⁴ 项在这里发力。for k 1:numel(chList) c chList(k); mem slot{k}; if isempty(mem) state.E(c) 0 continue; % 空簇且已死亡跳过 end % 成员上行 for m mem d D(m, c); state.E(m) state.E(m) - tx_energy(prm.k_data, d, prm); state.E(c) state.E(c) - rx_energy(prm.k_data, prm); end % 簇头聚合 state.E(c) state.E(c) - prm.E_DA * prm.k_data * max(numel(mem), 1); % 簇头向基站发送 dBS D(c, end); state.E(c) state.E(c) - tx_energy(prm.k_data, dBS, prm); end这里有个细节值得注意max(numel(mem), 1)保证空簇成员全死时仍有 1 个包的聚合开销而不是 0。空簇意味着簇头自己也要上传自己的数据聚合开销不该消失。4.3 死亡判定与每轮统计量能量扣完之后必须立刻判死且判死要放在所有节点的能量都更新完之后否则先判死的节点会影响后面对它的能量扣减function state check_death(state) newlyDead state.E 0 state.dead 0; state.dead(newlyDead) 1; state.E(state.dead 1) 0; % 防止负能量累积 state.nAlive sum(state.dead 0); endstate.E(state.dead 1) 0这一行看着多余实则是必须的。节点死亡后如果后续某轮它仍然被当成簇成员参与能量扣减能量会变成 -0.3 J 这样的负值累加总剩余能量时这些负值会把总量抵消掉画出来的“全网剩余能量”曲线出现反常抬升。把负能量截断为 0 是标准做法代价是可忽略的数值误差。每轮统计我一般记这几个量后面画图和做对比全靠它们stat [r, state.nAlive, sum(state.E), sum(state.chFlag), numel(chList)];依次是轮次、存活节点数、全网剩余能量、簇头数、实际有成员的簇数。第三个和第五个常被忽略——簇头数和“有效簇数”不是一回事当出现空簇时两者会分叉这个差值本身就是诊断信息。4.4 让整轮函数保持无副作用把上面四段拼成run_round时坚持一个原则只通过入参和返回值通信不读全局变量。好处是你可以随时拿某轮的中间状态单独重跑、对比、打日志而不用重跑前面几百轮。function [state, stat] run_round(state, prm, bs) D dist_matrix(state.pos, bs); chList find(state.chFlag); % ……成簇、TDMA、稳态、判死 state check_death(state); stat [sum(state.nAlive), sum(state.E), numel(chList)]; end这种写法在调试阶段的回报极高。有一次我发现第 340 轮之后曲线斜率突变直接把该轮的状态存成 mat 文件单独喂给run_round重跑十几遍十秒定位到是空簇的聚合开销写成了 0。要是每次都得从第 1 轮跑起这个 bug 能吃掉一整天。5. 避坑与排查五条踩过的血泪记录仿真代码的特点是“不报错但结果不对”下面五条都是我在复现和给别人看代码时反复遇到的按出现频率排序。5.1 曲线前几轮就掉节点现象第 3 到第 5 轮开始出现死亡节点和文献里“前 200 轮全部存活”对不上。原因初始能量 E0 单位搞混。文献里常写 0.5 J 或 2 J也有人写 500 mJ、写成 0.5 却按焦耳以外的量级在算更隐蔽的是能量模型里的 nJ 没转成 JE_elec 50而不是50e-9导致单次发送能耗是实际的一亿倍。解决跑之前先打印一次单轮总能耗按 100 节点、5 个簇头估算应在 0.10.5 J 量级如果打出来是 1e7 数量级一定是单位没转。把prm初始化单独抽成一个函数所有量都带e-9、e-12后缀杜绝手写漏乘。5.2 簇头数每轮都远大于预期现象设了 p0.05、100 个节点每轮却选出十几个甚至几十个簇头。原因阈值公式的分母写错或者mod(r, cool)里 cool 用了1/p的浮点值导致分母接近 0单节点概率被放大。另一种情况是冷却期判断写反r - lastCH cool当成每年节点每轮都重新参与竞争。解决在选举函数里加一行断言跑完检查sum(chFlag)是否落在[0.5×n×p, 2×n×p]区间超出就打日志。同时把cool强制取整mod的两个入参都保证是整数。5.3 存活曲线是一条直线下坠现象存活节点数从第一轮开始匀速下降画出来是标准斜线而不是文献里“先平后陡”的形态。原因通常是没有做能量耗尽节点退出选举的处理或者判死后没有把E截断为 0。前者让死节点继续被选为簇头后者让负能量在总和里抵消两种情况都会把曲线压成直线。解决确认select_ch第一行有if E(i) 0, continue, end确认check_death里有能量截断。再检查簇成员分配时是否还包含死节点——发现死节点仍在接收广播扣能量说明find(~chFlag E 0)漏了E 0这个条件。5.4 同样的参数两次运行曲线不一样现象不改任何参数连跑两次首节点死亡轮次差了二三十轮。原因rand()没有固定种子。这不是 bug但对比实验必须有可复现性尤其当你要证明“改进算法比 LEACH 好 8%”时随机波动可能就有 8%。解决每次实验开头rng(seed)且不同算法用同一组种子。做多组实验取平均时建议 20 组以上种子固定成1:20这样别人拿你的代码能跑出完全一致的图。seeds 1:20; res zeros(numel(seeds), 3); % FND / HND / LND for s 1:numel(seeds) rng(seeds(s)); res(s, :) run_experiment(prm, p, rounds); end fprintf(FND %.1f ± %.1f\n, mean(res(:,1)), std(res(:,1)));5.5 簇头到基站的能耗看起来不合理地低现象簇头明明离基站 150 m按 d⁴ 算应该有明显开销但日志里簇头每轮只掉一点点能量。原因d0被手填成了 200 或者忘了用公式算导致 150 m 仍走自由空间分支或者距离矩阵里基站坐标拼错了位置D(c, end)取到的其实是某个普通节点。解决打印prm.d0确认约等于 87.7单测一次tx_energy(4000, 150, prm)手算应为 4000×50e-9 4000×0.0013e-12×150⁴ ≈ 2.0e-4 2.63e-3 ≈ 2.83e-3 J如果打出来是 1e-4 量级就是走了自由空间分支。基站坐标建议在dist_matrix里最后一项拼入并用size(D)断言等于n1。提示这五条里有四条都能靠“单轮手算 断言”提前拦住别等跑完 1000 轮再回头查。6. 从单轮到多轮批量指标落盘与三条值得试的改法单轮跑通只是热身真正拿去做对比要的是批量实验和可复现的指标。我通常固定三件事种子、轮数、指标定义。指标就取三个——首节点死亡轮次FND、半数节点死亡轮次HND、最后节点死亡轮次LND再补一条全网剩余能量曲线。function [fnd, hnd, lnd, eCurve] run_experiment(prm, p, rounds, n) % 初始化网络并跑满 rounds 轮返回三个死亡轮次与能量曲线 pos 100 * rand(2, n); % 100m×100m 区域 bs [50; 150]; % 基站位于区域外 E prm.E0 * ones(1, n); lastCH -inf(1, n); eCurve zeros(1, rounds); fnd rounds; hnd rounds; lnd rounds; for r 1:rounds chFlag select_ch(E, lastCH, p, r); if ~any(chFlag) any(E 0) alive find(E 0); chFlag(alive(randi(numel(alive)))) true; end lastCH(chFlag) r; state struct(pos, pos, E, E, chFlag, chFlag, ... lastCH, lastCH, dead, zeros(1, n)); [state, ~] run_round(state, prm, bs); E state.E; nAlive sum(state.dead 0); eCurve(r) sum(E); if fnd rounds nAlive n, fnd r; end if hnd rounds nAlive n/2, hnd r; end if lnd rounds nAlive 0, lnd r; break; end end end逻辑上要注意fnd的初值设为rounds而不是 0这样“没死过节点”和“第 1 轮就死”能区分开避免把未死亡误判成第 0 轮。hnd用 n/2节点数为奇数时向下取整比 n/2稳。参数上rounds至少给到 LND 的两倍否则最后一段曲线被截断看不出全网死亡后的收敛形态。基站位置放在区域外比如 y150是常见设定它制造的远距离上行正是簇头负载集中的来源也是后面所有改进型算法的着力点。三条值得动手改的方向按性价比排序一是把簇头到基站改成多跳中继选离基站更近的簇头转发对 d⁴ 项的削减立竿见影二是引入剩余能量因子把阈值里的 p 换成 p×(E/E0)让高能量节点多当簇头这条改动只有一行但效果明显三是给每轮加上真实的时隙时长把“每轮时间”也统计出来因为很多场景比的是“单位时间内的存活轮数”而不是纯粹的轮数。最后说个习惯上的教训。早些年我做对比实验图一出来发现改进算法只在某些种子下领先换个种子就落后当时第一反应是“改进不够显著”差点把方向推翻。后来把种子固定成 20 组、把方差打进表格才发现原方案在 20 组里只有 3 组领先均值差被两组异常值撑起来的。从那以后我每次跑对比实验都强制走一遍“多种子 标准差 原始曲线截图留档”不再看单次结果下结论。希望帮到你。本文还有配套的精品资源点击获取