雷电流幅值累积概率分布宏观上可体现雷电活动规律,是工程上雷击计算分析的重要参数之一[1-6]。获取准确的雷电流幅值累积概率分布,对于认知雷电活动规律和特征、准确评价输电线路和杆塔的防雷性能、指导防雷工程设计有重要意义[7-14]。

我国电力规程参考了前苏联对数形式的雷电流幅值概率方程,早期由于雷电数据匮乏,利用新杭线磁钢棒记录的雷电流幅值数据反演得出了方程的参数,且被长期使用[15-16]。国外学者提出指数分式形式的雷电流幅值累积概率分布公式,IEEE采用了这一形式并在对全球雷电流幅值研究基础上给出了推荐参数,这一结果被广泛接受[14,17]。随着广域雷电地闪监测系统覆盖全国,且绝大部分雷电探测站运行时间超过10年,已积累了海量雷电地闪监测数据,从而具备了开展雷电参数时空差异性统计研究的条件[17-19]。

广域雷电地闪监测系统针对特定地区或线路走廊记录每一次地闪雷电流幅值及地闪总次数,由此可按定义算出大于任意幅值的雷电流占地闪总次数的百分比,实际应用中通常每隔数kA选作一个计算点,计算结果作为雷电流幅值累积概率分布的统计散点。而由散点得到分布曲线的拟合过程,相关研究成果较少。文献[17]提出了将雷电流幅值分成3个区间、使用多项式分段拟合的方法,最后在分段点处进行连续化处理,得到与IEEE推荐形式较一致的曲线。

本文以IEEE推荐公式为原型,提出3种雷电流幅值累积概率分布曲线的拟合方法,结合实例对比分析拟合结果,Levenberg-Marquardt法拟合效果最佳,可切实提高分布公式拟合参数的准确度,同时所需运算量也最大,提出的混合拟合方法可在保持拟合节果不变的条件下进一步减少计算开销,运算时间较Levenberg -Marquardt法减少约26%。

IEEE推荐的雷电流幅值累积概率分布公式为

式中：I为雷电流幅值,kA;P(>I)为雷电流幅值超过I的概率;a、b为待定参数,IEEE推荐为a=31、b=2.6,对应的曲线为图1。

图1 IEEE推荐雷电流幅值累积概率分布曲线 Fig.1 Curve of cumulative probability distribution (CPD) of lightning current magnitude (LCM) recommended by IEEE

IEEE推荐参数是综合全球雷电流幅值的平均结果,雷电活动随时间空间变化的差异性很大,不宜直接采用该平均结果[10],实际应用中往往根据某个时间段内某个特定地区或线路走廊的雷电地闪监测数据拟合得到。I的取值,理论上是(0, ∞),但自然雷电地闪中200 kA以上的雷电流已属罕见,因此雷电流幅值累积概率公式的拟合,本质是利用一个区间如[0, 200]或[0, 300]等之上的观测统计散点近似代表(0, ∞)上的分布曲线。实际观测统计中,a的取值为10~40,较常见的是20~40;b的值大于1,较常见的是2~3。

定义函数,即本文要进行拟合的原型函数为

式中：自变量x表示雷电流幅值,x≥0;因变量y表示雷电流幅值超过x的概率;a、b为待定参数。在已获得一定数量的雷电地闪监测数据情况下,统计计算得出I从小至大顺序取不同值时对应的P(>I)值,即取得N个统计散点(X1, Y1)、(X2, Y2)、···、(XN, YN)。统计散点应尽可能接近最终拟合曲线,判断这种“接近程度”的准则,常用为最小二乘准则,其误差定义为

式中：R为N个散点的拟合残差平方和(RSS);(Xj, Yj)为第j个散点;j为散点的序号;y(Xj)为y(x)在Xj处的函数值。

记a的拟合值为a*,b的拟合值为b*,显然若a*、b*使R值越小,说明该拟合方法越好。

对于式(2),a的物理意义为中值电流,即P(>a)=0.5,以纵轴取值0.5画水平线与函数曲线相交,交点横坐标则为a。b的物理意义为曲线陡度,b值越大曲线越陡、衰减越快。“中值电流”法是基于a的物理意义,在N个散点中确定纵坐标最接近0.5的点,求解出a*值,再计算出b*值。

在N个散点中,若恰好存在Yk=0.5（通常很难出现）则直接拟合得出a*=Xk。若不存在Yk=0.5的散点,则在0.5附近找出最近的(Xm, Ym)、(Xm+1, Ym+1)两点,且满足Ym>0.5>Ym+1,过(Xm, Ym)、(Xm+1, Ym+1)两点作直线,与y=0.5交点横坐标即认为是a*,即在局部线性化找出中值电流,过程如图2所示,线性方程由式(4)求解,即

由式(2)可知,在任意散点处代入坐标值和已拟合得出的a*值,均可计算出一个b的拟合值b*,综合考虑求取平均值较好,即

式中：j为散点序号;(Xj, Yj)为第j个散点坐标;a*为上一步通过中值电流得到的参数a的拟合值;b*为本方法所得参数b的拟合值。

曲线拟合中较常见的是线性拟合,即线性回归。对于雷电流幅值累积概率分布公式这种较复杂的非线性函数,可线性变换之后再进行拟合。

将式(2)变换得到式(6),定义线性函数如式(7)：

式中：x≥0;w为自变量;u为因变量;λ、η为待定参数。

对比式(6)和式(7)的相似性,将前文所述N个散点按照式(8)、(9)变元换算,得出变换后N个新的散点(W1, U1)、(W2, U2)、···、(WN, UN)。式(8)、(9)为：

式中：(Xj, Yj)为前文所述第j个散点坐标;(Wj, Uj)为经换元后第j个散点坐标。

然后以式(7)为原型进行线性拟合,计算得出拟

图2 求解中值电流示意图 Fig.2 Sketch of finding median current

合值λ*、η*,线性拟合过程不再赘述。最后由式(10)、(11)反变换得出a*、b*：

式中：λ*、η*分别为参数λ、η的拟合值;a*、b*为本方法所得参数a、b的拟合值;e为自然常数。

Levenberg-Marquardt法是较常用的一种非线性最小二乘问题求解方法[20-22],其在Gauss-Newton法基础上引入阻尼因子,加强算法的鲁棒性。

记M个拟合参数组成向量p,N个观测散点的横、纵坐标分别组成向量为X、Y,在每个观测点处与拟合函数值存在相应的残差,N个残差函数组成向量如式(12)、(13)：

rj(p)=Yj-y(p, Xj) (12)

r(p)=[r1(p), r2(p),···, rN(p)]T (13)

式中：y(p, Xj)为代入参数向量p和Xj后的函数值;rj(p)为第j个散点处拟合残差函数;r(p)为残差函数向量。

定义目标函数F(p)为

式中rT(p)为r(p)的转置。可见,F(p)恰为拟合残差平方和R的1/2,在最小二乘原则下,使F(p)值最小的向量p*即为最优解,残差平方和R亦最小。具体到式(2)所示原型函数,向量p=[a, b]T,M=2。

残差向量r(p)对p中每个变量的偏导数组成N×M阶矩阵,即Jacobi矩阵,在本文式(2)原型函数、p=[a, b]T、M=2条件下,其计算式为

J(p)=[J1(p) J2(p)] (15)

式中：J(p)为Jacobi矩阵;J1(p)、J2(p)为J(p)的两个子矩阵。

Levenberg-Marquardt算法的核心是：在给定参数向量初始值p0基础上,计算寻找合适的增量Δp使目标函数F(p)沿值减小方向变化,并令p=p+Δp迭代至逼近p*,F(p)达到最小值,其中增量Δp由方程(18)、(19)求解得出,即：

S(p)=JT(p)J(p)+μdiag(JT(p)J(p)) (18)

S(p)Δp= -JT(p)r(p) (19)

式中：μ为阻尼因子;diag(JT(p)J(p))表示由JT(p)J(p)主对角元素组成的对角矩阵;Δp为向量p的迭代增量;S(p)为中间变量。

Levenberg-Marquardt算法在每轮迭代中可有多种指标衡量本轮迭代质量的优劣,并据此调节μ值。本文采用最简单的做法,仅比较F(p+Δp)和F(p)大小,决定Δp是否可接受。Levenberg-Marquardt法参数向量p的解算流程如图3,实施步骤如下：

1）根据原型函数和已知散点,构造残差函数向量r(p)与目标函数F(p)。

2）初始化,设置初始点p=p0、向量计算精度ε,阻尼因子μ、阻尼因子变化倍率ν。

3）向量及矩阵计算,代入p值计算残差向量r(p)、Jacobi矩阵J(p),如式(10)和(12)。

4）求解增量向量,代入p值计算矩阵S(p)=JT(p)J(p)+μdiag(JT(p)J(p)),构造增量正规方程S(p)Δp= -JT(p)r(p),求解出Δp。

5）精度判断,若满足|Δp| 图3 Levenberg-Marquardt算法简要流程 Fig.3 Brief flow chart of Levenberg-Marquardt algorithm

≥F(p),则令μ=μν放大阻尼因子,并跳转至步骤4）。

计算结束时,p*=[a*, b*]T,即a*=\(p_{1}^{*}\)、b*=\(p_{2}^{*}\),其中\(p_{1}^{*}\)、\(p_{2}^{*}\)分别为向量p*的第1和第2个元素。

本文使用2005—2015年全国雷电地闪雷电流幅值和每年地闪总数,以5 kA为间隔,统计计算得出雷电流幅值I=0、5、···、300 kA处P(>I)的值,每年的雷电流幅值累积概率分布共61个散点（即N=61）,按雷电流幅值由小至大排序。再分别使用前文所述方法对这些统计散点进行曲线拟合。

“中值电流”法求解a*较容易,求解b*时引用ln(Yj-1–1)和ln(Xj/a*)的值,为使计算结果有意义以及数值计算时存在数据截断,需剔除掉Yj为1（对应的Xj为0或接近0）和Yj接近0的点,本文以1-Yj 图4 3种拟合方法的残差平方和 Fig.4 Comparison of residual square sum by 3 methods

图5 2010年雷电流数据散点及3种方法拟合曲线 Fig.5 Scatter data of CPD of LCM in 2010 and fitting curves by 3 methods

表1 3种方法拟合结果对比 Table 1 Comparison of results from 3 methods

2N次对数、2N次乘方、N次y(x)函数运算,求解增量向量Δp需1次矩阵加法、3次矩阵乘法和1次矩阵求逆运算,而矩阵乘法和求逆涉及的乘法或加法运算则为N2级次。“中值电流”法和线性变换法计算量较小,而Levenberg-Marquardt法中每一轮迭代的计算开销就已显著大于前两者,在具体实践中可将“中值电流”法的拟合结果作为初始值代入Levenberg-Marquardt法中,加速迭代过程,形成混合的拟合方法,减少计算开销。几种方法的运算量对比见表2。

对于2005—2015年每年全国雷电地闪数据统计算得到的61个散点,即N=61,在笔者Intel Core i3-2348 2.30GHz主频CPU、4GB内存、Windows 10操作系统配置及MATLAB运行环境下,上述几种拟合方法11次运算的平均时间均在ms级：“中值电流”法平均花费9.17 ms、线性变换法平均花费9.47 ms、Levenberg-Marquardt法平均花费17.15 ms,“中值电流”与Levenberg-Marquardt混合法平均花费12.68 ms,较Levenberg-Marquardt法减少约26%。

1）雷电活动随时间空间变化的差异性非常大,实际应用中需要利用雷电地闪监测数据的统计结果拟合得出雷电流幅值累积概率分布的参数,提出以拟合残差平方和大小作为衡量观测散点接近拟合曲线的标准,拟合结果残差平方和越小,该拟合方法越好。

2）提出3种可用于雷电流幅值累积概率分布曲线的拟合方法,即“中值电流”法、线性变换法和Levenberg-Marquardt法,前两者计算过程简单易于实现,Levenberg-Marquardt法鲁棒性强,计算开销也最大。

3）计算结果表明,Levenberg-Marquardt法的拟合残差平方和最小,方法最优,“中值电流”法次之,线性变换法的拟合残差平方和显著高于前两者,高出1~2个数量级,得出的a*、b*值也显著大于前两者,3种方法的拟合曲线及散点在雷电流幅值I>150 kA后几乎重合,差异主要存在于I 表2 不同方法计算开销对比 Table 2 Comparison of computation cost by different methods

在保持拟合效果不变的条件下较Levenberg -Marquardt法减少约26%的运算时间。

[1] 赵淳,雷梦飞,陈家宏,等. 输电线路走廊雷电流幅值分布统计方法[J]. 高电压技术,2017,43(5):1609-1614. ZHAO Chun, LEI Mengfei, CHEN Jiahong, et al.Statistical method of lightning current amplitude distribution along transmission line corridor[J]. High Voltage Engineering, 2017, 43(5): 1609-1614.

[2] 陈家宏,赵淳,谷山强,等. 我国电网雷电监测与防护技术现状及发展趋势[J]. 高电压技术,2016,42(11):3361-3375. CHEN Jiahong, ZHAO Chun, GU Shanqiang, et al.Present status and development trend of lightning detection and protection technology of power grid in China[J]. High Voltage Engineering, 2016, 42(11): 3361-3375.

[3] 李瑞芳,吴广宁,曹晓斌,等. 输电线路雷电绕击率的三维计算方法[J]. 电工技术学报,2009,24(10):134-138. LI Ruifang, WU Guangning, CAO Xiaobin, et al.Three-dimensional calculation method on shielding failure rate of transmission lines[J]. Transactions of China Electrotechnical Society, 2009, 24(10): 134-138.

[4] 赵淳,任华,黄修乾,等. 500 kV紧凑型耐张杆塔绕击耐雷性能精细化分析方法[J]. 高电压技术,2017,43(8):2746-2753. ZHAO Chun, REN Hua, HUANG Xiuqian, et al.Elaborate method for analyzing shielding lightning performance of 500 kV tension tower in compact transmission line[J]. High Voltage Engineering, 2017, 43(8): 2746-2753.

[5] 岳平,吴敏,谭进,等. 基于“杆塔圆法”的雷电参数统计与应用[J]. 电瓷避雷器,2013,35(1):107-111. YUE Ping, WU Min, TAN Jin, et al.Statistics and application of lightning parameters based on tower circle method[J]. Insulators and Surge Arresters, 2013, 35(1): 107-111.

[6] 江安烽,李锐海,牛萍,等. 后续雷击对输电线路绕击耐雷性能的影响研究[J]. 上海交通大学学报,2015,49(4):411-417. JIANG Anfeng, LI Ruihai, NIU Ping, et al.Impact of subsequent strokes on lightning shielding failure outage rate of transmission lines[J]. Journal of Shanghai Jiao Tong University, 2015, 49(4): 411-417.

[7] 李如虎. 输电线路雷电绕击及其防雷研究[J]. 南方电网技术,2009,3(1):53-57. LI Ruhu.Study on shielding failure of the power transmission line and its lightning strike proof[J]. Southern Power System Technology, 2009, 3(1): 53-57.

[8] 张勇军,黄慧,刘海洋,等. 含分布式水电的电网雷击风险评估[J]. 华南理工大学学报:自然科学版,2011,39(10):127-131. ZHANG Yongjun, HUANG Hui, LIU Haiyang, et al.Lightning strike risk evaluation of power grid with distributed hydropower[J]. Journal of South China University of Technology (Natural Science Edition), 2011, 39(10): 127-131.

[9] 赵芝,石季英,袁启海,等. 输电线路的雷击跳闸概率预测计算新方法[J]. 电力系统自动化,2015,39(3):51-58,139. ZHAO Zhi, SHI Jiying, YUAN Qihai, et al.A new prediction calculation method of transmission lines lightning trip probability[J]. Automation of Electric Power Systems, 2015, 39(3): 51-58, 139.

[10] 何金良,杨滚,余占清,等. 用于雷电防护的雷电流波形参数研究[J]. 建筑电气,2017,36(3):3-8. HE Jinliang, YANG Gun, YU Zhanqing, et al.A study on lightning current waveform parameters for lightning protection[J]. Building Electricity, 2017, 36(3): 3-8.

[11] 苏杰,谷山强,赵淳,等. 高速铁路牵引网暴露弧面法直击雷跳闸率计算[J]. 中国电机工程学报,2015,35(19):5103-5113. SU Jie, GU Shanqiang, ZHAO Chun, et al.Calculation of direct lightning stroke trip-out rate with exposure globoid method for high-speed railway traction network[J]. Proceedings of the CSEE, 2015, 35(19): 5103-5113.

[12] 何金良. ±800 kV云广特高压直流线路雷电防护特性[J]. 南方电网技术,2013,7(1):21-27. HE Jinliang.The lightning protection characteristics of ±800 kV Yunnan-Guangdong UHVDC transmission line[J]. Southern Power System Technology, 2013, 7(1): 21-27.

[13] 陆俊杰,徐秀峰,梁聪. 基于正态分布雷电参数的220 kV线路雷击跳闸模拟[J]. 电瓷避雷器,2015,37(5):106-110. LU Junjie, XU Xiufeng, LIANG Cong.Simulation of lightning trip-out of 220 kV transmission lines based on normal distribution of lightning parameters[J]. Insulators and Surge Arresters, 2015, 37(5): 106-110.

[14] 李瑞芳,吴广宁,曹晓斌,等. 雷电流幅值概率计算公式[J]. 电工技术学报,2011,26(4):161-167. LI Ruifang, WU Guangning, CAO Xiaobin, et al.Formula for probability of lightning current amplitude[J]. Transactions of China Electrotechnical Society, 2011, 26(4): 161-167.

[15] 交流电气装置的过电压保护和绝缘配合:DL/T 620—1997[S]流电气装置的过电压保护和绝缘配合:DL/T 620—1997[S]. 北京:中国电力出版社,1997. Overvoltage protection and insulation coordination for AC electrical installations: DL/T 620—1997[S] DL/T 620—1997[S]. Beijing, China: China Power Press, 1997.

[16] 孙萍,郑庆军,吴璞三,等. 220 kV新杭线I回路27年雷电流幅值实测结果的技术分析[J]. 中国电力,2006,39(7):74-76. SUN Ping, ZHENG Qingjun, WU Pusan, et al.Measurement and calculation of lightning flashover rate of 220 kV Xinhang no.1 transmission line[J]. Electric Power, 2006, 39(7): 74-76.

[17] 刘刚,张弦,陈锡阳,等. 雷电流幅值概率分布函数的分段拟合方法[J]. 华南理工大学学报:自然科学版,2014,42(4):40-45. LIU Gang, ZHANG Xian, CHEN Xiyang, et al.Sectioned fitting method of probability distribution function of lightning current amplitude[J]. Journal of South China University of Technology (Natural Science Edition) , 2014, 42(4): 40-45.

[18] 余宏桥,周强明,邵立政,等. 湖北电网雷电流幅值概率分布特性分析[J]. 陕西电力,2014,42(7):49-51. YU Hongqiao, ZHOU Qiangming, SHAO Lizheng, et al.Analysis of distribution characteristics of lightning current magnitude in Hubei power network[J]. Shaanxi Electric Power, 2014, 42(7): 49-51.

[19] 李家启,王劲松,申双和,等.基于ADTD系统监测的雷电流幅值累积概率特征分析[J]. 气象,2011,37(2):226-231. LI Jiaqi, WANG Jingsong, SHEN Shuanghe, et al.Analysis of the cumulative probability of lightning current amplitudes based on the statistical data obtained from the ADTD system[J]. Meteorological Monthly, 2011, 37(2): 226-231.

[20] MADSEN K, NIELSEN H B, TINGLEFF O. Methods for non-linear least squares problems, 2nd edition[E/OL]. [2004-08-31]. .

[21] 严正,范翔,赵文恺,等. 自适应Levenberg-Marquardt方法提高潮流计算收敛性[J]. 中国电机工程学报,2015,35(8):1909-1918. YAN Zheng, FAN Xiang, ZHAO Wenkai, et al.Improving the convergence of power flow calculation by a self-adaptive levenberg-marquardt method[J]. Proceedings of the CSEE, 2015, 35(8): 1909-1918.

[22] 齐丽岩,肖现涛,张立卫. 求解半光滑方程组的LM方法收敛性分析[J]. 大连理工大学学报,2015,55(5):548-552. QI Liyan, XIAO Xiantao, ZHANG Liwei.Convergence analysis of LM method for semismooth equations[J]. Journal of Dalian University of Technology, 2015, 55(5): 548-552.