MATLAB输电线路建模:从单位长度参数到频率扫描与瞬态模拟
2026/9/14 14:40:37 网站建设 项目流程

简介:这套MATLAB工具箱专注于架空与地下输电线路建模,完整覆盖单位长度参数计算、传播特性分析、频率扫描以及瞬态模拟等核心功能。用户输入线路尺寸与材料参数后,即可快速获得电阻、电感、电容、电导等关键参数,并通过频率扫描评估不同频段下的线路响应,为输变电工程设计与故障分析提供数据支撑。资源面向电力系统工程师、研究人员以及相关专业本硕学生,既可用于日常仿真分析,也可直接支撑课程设计、期末大作业或毕业设计。包内共265个文件,以M脚本/函数为主(232个m),配有13个mat数据文件、8个PDF说明文档及若干txt配置说明,压缩包整体约73.19MB。程序采用参数化编程,代码结构清晰、注释详细,附赠可直接运行的案例数据;同时,压缩包内还包含PDF等辅助学习材料,便于用户对照理解建模流程与仿真结果。目前已有120人学习/下载,适合希望快速掌握MATLAB输电线建模方法并自行扩展仿真场景的读者。

1. 输电线路建模工具箱:从单位长度参数到频率扫描与瞬态模拟的一条链路

MATLAB 的输电线路建模工具箱,最终要交付的是一套能把架空线和地下电缆从单位长度参数一路算到频率扫描和瞬态模拟的数据流。架空线路的导线多、几何位置开阔,地下电缆则要考虑绝缘层、护套和土壤回流,两者的单位长度阻抗和导纳形式差别很大,但到传播常数 γ 和特征阻抗 Zc 之后,频域求解的套路又统一起来了。这个工具箱适合两类人:一类是电磁暂态仿真前想自己校核线路参数的工程师,另一类是研究宽频等效模型、需要批量扫频和做相位修正的学生。我一般会先把参数计算按“几何输入—一次参数—二次参数—时域”分层,而不是把所有逻辑塞进一个脚本;这样改土壤电阻率、换电缆绝缘材料,只需要改输入参数结构体,后面三四层代码都不动。

2. 单位长度参数:架空线的复深度法与地下电缆的同轴分层公式

2.1 为什么单位长度参数不能当常数处理

输电线路的 R、L、C、G 在教科书里常给成 50Hz 的定值,但工具箱一旦做频率扫描,就绕不开集肤效应和大地回流。架空线的电阻不是 Rdc 一直不变,频率升高后电流向导体表面集中,交流电阻会往上走;大地回流路径也会随频率变浅,等效电感跟着降。地下电缆更明显:金属护套既是结构件又是电磁回路,绝缘层里的位移电流在几十 kHz 以上不能忽略,介质损耗对应的 G 随频率近似线性增长。

因此单位长度参数的计算函数,第一个入参必须是频率 f,而不是“默认工频”。串联阻抗写成 Z_series = R + jωL,并联导纳写成 Y_shunt = G + jωC。一次参数算准了,后面的 γ、Zc、频率扫描和瞬态模拟才谈得上可靠。

2.2 架空线:Carson 复深度公式的 MATLAB 实现

架空线路单位长度参数最常用的做法,是 Carson 大地回流公式。完整公式含无穷积分,写进工具箱有点重;工程上更常见的是复深度近似:把大地回流等效到地下一个复数深度 p 的地方,然后按镜像法计算。下面这个函数实现了单根导线-大地回路的最小版本。

function [R,L,C,G] = overhead_params(h, r, rho_c, rho_g, f) % 架空线单位长度参数(单根导线-大地回路) % h : 导线到地面平均高度, m % r : 导线半径, m % rho_c : 导线电阻率, Ω·m % rho_g : 土壤电阻率, Ω·m % f : 频率标量或向量, Hz mu0 = 4e-7*pi; eps0 = 8.8541878128e-12; w = 2*pi*f; idx0 = (w == 0); w(idx0) = 2*pi*1e-9; % 先用小频率占位,最后补 0Hz 极限 Rdc = rho_c / (pi*r^2); % 1) 圆截面导体的内部阻抗:Bessel 解 k = sqrt(1j*mu0*w / rho_c); Zi = zeros(size(w)); for m = 1:numel(w) kr = k(m)*r; Zi(m) = Rdc * kr * besseli(0, kr) / (2 * besseli(1, kr)); end % 2) 大地回路的复深度近似 p = sqrt(rho_g ./ (1j*mu0*w)); % 复数等效深度 Zg = 1j*mu0*w/(2*pi) .* log((2*(h+p))./r); Z = Zi + Zg; R = real(Z); L = imag(Z) ./ w; C = 2*pi*eps0 ./ log(2*h/r) * ones(size(w)); G = zeros(size(w)); % 3) 0Hz 直流极限 if any(idx0) R(idx0) = Rdc; L(idx0) = mu0/(2*pi)*(log(2*h/r) + 0.25); end end

这个函数里,p是复数,所以log的结果也是复数;取real(Z)得到考虑土壤损耗后的电阻,imag(Z)/w得到包含大地回流的电感。besseli(0,kr)/besseli(1,kr)在高频极限下趋近 1,不会因为集肤效应发散。要注意两个参数:h是导线对地等效高度,不是相间距离;rho_g是土壤电阻率,不是接地电阻,常见干燥土壤在 10~1000 Ω·m 之间,这个值对低频扫描的实部影响很大。

2.3 地下电缆:同轴绝缘层和护套回路

地下电缆和架空线的几何模型不同。常用建模方式是把单芯电缆拆成芯线、绝缘层、金属护套、外护套和大地几层,其中“芯线-绝缘层-护套”构成一个同轴回路,护套又是下一次大地回流的导体。先算主回路,参数公式如下:

function [R,L,C,G] = cable_params(rc, rs, epr, tand, f, rho_c) % 单芯同轴电缆:芯线-绝缘层-金属护套 % rc : 芯线半径, m % rs : 绝缘层外半径/护套内半径, m % epr : 绝缘材料相对介电常数 % tand : 介质损耗角正切 % f : 频率, Hz % rho_c: 芯线电阻率, Ω·m mu0 = 4e-7*pi; eps0 = 8.8541878128e-12; w = 2*pi*f; Rdc = rho_c/(pi*rc^2); % 芯线内部阻抗,Bessel 解与架空线相同 k = sqrt(1j*w*mu0/rho_c); Zi = zeros(size(w)); for m = 1:numel(w) kr = k(m)*rc; Zi(m) = Rdc * kr*besseli(0,kr)/(2*besseli(1,kr)); end % 芯-护套回路的电感与电容 L = mu0/(2*pi)*log(rs/rc); C = 2*pi*eps0*epr/log(rs/rc); % 介质损耗 G = w .* C * tand; R = real(Zi); end

这段代码把护套看成理想导体,实际工程里金属护套也有电阻和集肤损耗,严格做法是把护套电阻串进 R,再按 Pollaczek 公式叠加土壤回流阻抗;但对于 XLPE 电缆几十 kHz 以内的频率扫描,上述“芯-护套主回路”已经能反映主要相位特性。绝缘材料参数可以查表得到:XLPE 的 epr 约 2.3,tand 约 0.001;油纸绝缘的 epr 约 3.5~4.0,tand 稍大。

对象主要公式注意点
架空线电容C = 2πε0 / ln(2h/r)h 用等效对地高度
架空线大地回流p = sqrt(ρg/(jωμ0))f 接近 0 时取直流极限
电缆电感L = μ0/(2π) ln(rs/rc)rs 是护套内半径
电缆介质损耗G = ωC tanδ几十 kHz 后不能忽略

2.4 多相线路:从标量到 Z、Y 矩阵

实际输电线是三相加地线,单根导线公式只能算自参数。常见扩展方式是把每根导线的自阻抗放主对角线,互阻抗用“导体到另一根导体镜像的几何均距”计算:Zij = jωμ0/(2π) ln(Dij'/dij)。把导体按层排序,最后得到一个 N×N 的串联阻抗矩阵和并联导纳矩阵。工具箱里不要对每个元素写死公式,用一个双层for循环配合导线坐标数组生成矩阵,后续频率扫描时矩阵的每个频率点都要重算一次。

3. 传播特性与频率扫描:用 γ 和 Zc 把线路变成频域网络

3.1 二次参数:传播常数、特征阻抗、波速的最小实现

有了单位长度 R、L、C、G,线路在频域里就可以压缩成两个二次参数。传播常数 γ 决定幅值衰减和相位滞后,特征阻抗 Zc 决定端口阻抗匹配。MATLAB 里只需要三行,但输出变量要分开,方便后续画图和计算。

function [gamma, Zc, alpha_db, v] = secondary_params(R,L,C,G,f) % 根据单位长度参数计算传播特性 % 返回: gamma 传播常数(1/m), Zc 特征阻抗(Ω), % alpha_db 衰减(dB/m), v 相速度(m/s) w = 2*pi*f; Zser = R + 1j*w.*L; Yshu = G + 1j*w.*C; gamma = sqrt(Zser .* Yshu); Zc = sqrt(Zser ./ Yshu); alpha_db = real(gamma) * 8.685889638; % Np/m 转 dB/m v = w ./ imag(gamma); end

real(gamma)的单位是 Np/m,乘 8.6859 转成 dB/m,便于和厂商数据、实测曲线对比;imag(gamma)是相位常数,单位 rad/m,用它算相速度时不能只看频率,还要确认 β 没有跨过 π 分支。

提示:MATLAB 的sqrt对复数取实部非负的分支,正好满足被动线路 γ 实部为正的习惯。不要手动把相位翻到负半平面,否则扫出来的衰减会变负。

3.2 频率扫描主程序:logspace 选点与开路输入阻抗

频率扫描不需要均匀取点,线路参数在低频段随频率变化快,高频段相对平缓,用logspace按数量级分布更合理。

% 0.1 Hz 到 1 MHz,共 801 个点 f = logspace(-1, 6, 801); [R,L,C,G] = overhead_params(15, 0.01, 1.7e-8, 100, f); [gamma, Zc, alpha_db, v] = secondary_params(R,L,C,G,f); figure; subplot(3,1,1); loglog(f, alpha_db); grid on; ylabel('衰减 dB/m'); subplot(3,1,2); semilogx(f, abs(Zc)); grid on; ylabel('|Zc| Ω'); subplot(3,1,3); semilogx(f, v); grid on; ylabel('相速度 m/s'); xlabel('频率 Hz');

这个扫描脚本能回答三件事:衰减是否随频率单调上升、特征阻抗是否趋近某个高频常数、相速度在高频段是否接近光速。如果只需要看谐振点,再加一条开路输入阻抗曲线:

l = 50e3; % 线路长度 50 km Zin = Zc .* coth(gamma .* l); % 末端开路输入阻抗 semilogx(f, abs(Zin));

末端开路时,Zin 的谐振峰出现在imag(gamma)*l = n*pi的位置。这个判断对后续做保护测距、故障滤波都非常有用。

3.3 扫频场景里三个最常翻车的地方

第一个是低频段零点。频率取 0.1Hz 时,ω 很小但不为零,复深度 p 会很大,计算仍然成立;但很多函数在 f=0 处直接除零,所以工具箱入口处要有直流极限分支。第二个是alpha_db的单位。返回的是衰减常数,不是整条线路总衰减;要看 100km 线路的总衰减,得乘长度,再换算成 dB。第三个是 Zc 在高频可能不是纯阻性,它带很小的虚部,画图时最好同时画abs(Zc)angle(Zc),否则看不到频变特征。

4. 瞬态模拟:用频域传递函数反变换回到时域

4.1 为什么不在时域直接解偏微分方程

时域里架空线和电缆的偏微分方程带频变参数,直接差分会出现稳定性问题,还要处理电容、电感矩阵的非对角耦合。常见工具箱做法是留在频域:线路就是一个两端口网络,单位长度参数在每个频率点算好以后,传播项就是exp(-γl),加上源阻抗和负载阻抗的反射系数,可以写出端口代数方程。最后用 IFFT 回到时域。这个方法对任意频率相关的 R、L、C、G 都适用,代价是每个频率点都要重算一次 γ。

4.2 支持任意端接的频域-时域函数

下面这个函数实现了带反射的通用频域阶跃响应。开路由大电阻代替,短路由小电阻代替,匹配负载直接设为 Zc。

function [t, y] = line_transient_fd(h, r, rho_c, rho_g, l, dt, N, Zs, ZL) % 频域法求线路阶跃响应 % h,r,rho_c,rho_g : 架空线几何与材料参数 % l : 线路长度, m % dt: 采样间隔, s; N: FFT 点数 % Zs: 源内阻, Ω; ZL: 负载阻抗, Ω t = (0:N-1)*dt; f = (0:N-1)/(N*dt); f(1) = 1e-9; % DC 点占位 [R,L,C,G] = overhead_params(h, r, rho_c, rho_g, f); Zser = R + 1j*2*pi*f.*L; Yshu = G + 1j*2*pi*f.*C; gamma = sqrt(Zser .* Yshu); Zc = sqrt(Zser ./ Yshu); % 反射系数:从线路看向源端和负载端 Gamma_s = (Zs - Zc)./(Zs + Zc); Gamma_r = (ZL - Zc)./(ZL + Zc); T = Zc./(Zs + Zc); % 源端分压 H = T .* (1 + Gamma_r) .* exp(-gamma*l) ... ./ (1 - Gamma_s .* Gamma_r .* exp(-2*gamma*l)); u = ones(N,1); % 单位阶跃 U = fft(u); y = real(ifft(H(:) .* U)); end

调用时注意端接值:理想开路不要写Inf,用ZL = 1e12代替,否则(Inf-Zc)/(Inf+Zc)会出现 NaN;源内阻小到可以忽略时,用Zs = 1e-3而不是 0,既贴近实际断路器回路,也避免分压公式分母出现零。返回的波形可以看到首行波延时、反射叠加和稳态值三个特征。

端接方式ZsZL期望结果
源匹配-负载匹配ZcZc单程波,幅值为 E/2
零源-开路1e-31e12电压加倍,末端趋近 E
零源-短路1e-31e-3末端电压趋近 0,电流有冲击

4.3 瞬态结果自检:先看时延,再看稳态

写完瞬态函数,不要直接拿去算绝缘配合。先用 50Hz 的单位长度参数估一下波速和传播时延,再从仿真波形里找首波半幅值点,两者误差在 10% 以内才说明频率轴和端接方向没接反。

[R0,L0,C0,G0] = overhead_params(h, r, rho_c, rho_g, 50); v0 = 1/sqrt(L0*C0); expected_delay = l / v0; [~, idx] = min(abs(y - 0.5*max(y))); fprintf('预计时延 %.3f ms,波头半幅点 %.3f ms\n', ... expected_delay*1e3, t(idx)*1e3);

如果首个半幅点提前或延后很多,优先检查 FFT 的频率向量是否与输入信号长度对齐,其次检查负载端接方向。反射系数方向反了的表现是波头极性反,波形前半段出现下凹而不是上升。

5. 把散脚本封装成可维护的 MATLAB 工具箱

5.1 一个最小自检脚本:把理论极限变成断言

工具箱和散脚本的区别,在于每个函数都能单独验证。我用一个selftest函数把高频极限、直流极限和传播时延三条理论判据写成断言,每次改参数后先跑一遍。

function ok = selftest() f = logspace(0, 6, 101); [R,L,C,G] = overhead_params(15, 0.01, 1.7e-8, 100, f); % 直流极限 assert(abs(R(1) - 1.7e-8/(pi*0.01^2)) < 1e-12); % 高频极限下 Zc 应趋近 sqrt(L/C),架空线接近 300 欧级别 [~, Zc] = secondary_params(R,L,C,G,f); Lhf = 4e-7*pi/(2*pi) * log(2*15/0.01); Chf = 2*pi*8.854e-12 / log(2*15/0.01); zinf = sqrt(Lhf/Chf); assert(abs(abs(Zc(end)) - zinf)/zinf < 0.1); disp('selftest passed'); end

高频极限误差留 10% 的余量,是因为 Bessel 内部电感和大地损耗仍然让 Zc 带微小虚部。直流极限那一条卡死 Rdc,能及时发现输入单位写错。

5.2 用 Package 目录组织函数,避免命名冲突

把所有函数放进一个+lineToolbox目录,外部统一用lineToolbox.overhead_params(...)调用。这样不会和 Simscape、MATLAB 优化工具箱里的同名函数冲突,也方便整个目录分发给同事。发布时在 MATLAB 的 App 选项卡里选择 Package Toolbox,指定主函数和说明文件,生成的 .mltbx 可以直接安装,比打包成 .rar 更干净。

5.3 三个高频踩坑点最后确认一遍

第一,FFT 频率轴必须从 0 开始,写成(1:N)/(N*dt)会让所有频率偏一个 bin,瞬态波形会整体偏移。第二,零频点不能直接删,f(1)=1e-9只是占位,最终输出里 R、L 要在 f=0 处补直流极限。第三,开路端接用有限大电阻代替无限大,既避免 Inf/Inf,也让反射系数在数值上连续。

这套结构把单位长度参数、频率扫描和瞬态模拟拆成三层,但共用同一份输入参数结构体,测试入口只剩一个selftest。每次改完导体半径或土壤电阻率,先跑一遍自检,再去看扫描曲线和暂态波形,问题定位会快很多。

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

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

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

立即咨询