一种雷电流幅值累积概率分布曲线拟合计算方法与流程

本发明涉及雷电防护技术领域，更具体地说，涉及一种雷电流幅值累积概率分布曲线拟合计算方法。

背景技术：

雷电流幅值累积概率分布在宏观上体现雷电活动规律，也是开展雷击计算分析的重要参数之一。获取准确的雷电流幅值累积概率分布，对于认知雷电活动规律和特征、准确评价输电线路和杆塔的防雷性能、指导防雷工程设计至关重要。

早期由于雷电数据匮乏，我国利用新杭线磁钢棒记录的雷电流幅值数据，反演得出对数形式的雷电流幅值概率方程参数，并列入我国电力规程。目前在雷击计算分析中，常采用ieee推荐公式，如公式(8)。

公式(8)中，i为雷电流幅值，p(>i)为雷电流幅值大于i的概率，a、b为常数参数，ieee推荐值为a＝31、b＝2.6。a的物理意义为中值电流，即雷电流幅值大于a的概率为50％，b为表征雷电流幅值累积概率分布曲线陡度的参数，b值越大曲线越陡。

ieee推荐参数是综合上世纪全世界数百次自然雷电观测后拟合的总体结果，雷电活动随时间、空间变化的差异性非常大，不宜直接采用，实际应用中往往根据某个时间段内某个特定地区或线路走廊的雷电地闪监测数据统计拟合得到。

随着广域雷电地闪监测系统覆盖全国，且绝大部分雷电探测站运行时间超过10年，已积累了海量雷电地闪监测数据，具备了开展雷电参数时空差异性统计研究的条件。广域雷电地闪监测系统针对特定地区或线路走廊记录每一次地闪雷电流幅值及地闪总次数，由此可按定义计算出大于任意幅值的雷电流占地闪总次数百分比，实际应用中通常每隔数ka选择一个计算点，结果作为雷电流幅值累积概率分布的统计散点。由统计散点得到分布曲线的拟合过程，即曲线拟合，要保证拟合误差尽可能小，依赖于拟合方法的准确有效性。i的取值，理论上是(0,∞)，但自然雷电地闪中200ka以上的雷电流已属罕见，因此雷电流幅值累积概率公式的拟合，本质是利用一个区间如[0,200]或[0,300]等之上的观测统计散点近似代表(0,∞)上的分布曲线。实际观测统计中，a的取值在20～40之间较常见；b的值大于1，较常见是在2～3之间。

本发明旨在提出一种将雷电流幅值累积概率分布统计散点拟合为近似分布曲线的计算方法。

技术实现要素：

本发明要解决的技术问题在于，在n个统计散点前提下，以ieee推荐表达式为原型为原型函数、以最小二乘为准则，结合中值电流法与levenberg-marquardt(列文博格-马夸尔特)算法，获得一种雷电流幅值累积概率分布曲线拟合计算方法，使得拟合曲线与统计散点最为接近、残差平方和最小。

本发明解决其技术问题所采用的技术方案是：构造一种雷电流幅值累积概率分布曲线拟合计算方法，包括以下步骤：

步骤1)：将雷电地闪监测数据统计计算i取不同值时对应的p(>i)值所得n个散点按i从小至大排序，记此n个散点分别为(x1,y1)、(x2,y2)、....、(xn,yn)，(xj,yj)表示第j个散点，j为散点的序号，确定如公式(9)所示的原型函数和如公式(10)所示的最小二乘准则下拟合误差计算式，y(xj)为函数y(x)在xj处的值，r表示n个散点的拟合残差平方和，

；

步骤2)：采用中值电流法得到中值电流拟合值a1*，基于a1*反算出陡度参数的拟合值b1*；

步骤3)：构造如公式(11)所示的levenberg-marquardt残差向量r(p)与如公式(12)所示的目标函数f(p)，

；

针对公式(9)的二元参数，对应拟合向量p＝[a,b]^t、元数m＝2；残差函数向量r(p)表征在n个统计散点处的残差组成的向量，每个元素rj(p)为m元函数；在最小二乘准则下，使f(p)值最小的向量p*即为最优解，残差平方和亦最小；

步骤4)：levenberg-marquardt迭代初始化，设置初始点p＝p0、向量计算精度ε、阻尼因子μ、阻尼因子变化倍率ν，其中，p0由中值电流法的结果代入，即p0＝[a1*,b1*]^t；

步骤5)：残差向量r(p)及jacobi矩阵j(p)迭代，代入p值计算残差向量r(p)、jacobi矩阵j(p)，jacobi矩阵j(p)如公式(13)所示，其中p1表示向量p第1个元素、p2表示向量p第2个元素，

；

步骤6)：求解增量向量，代入p值计算矩阵s(p)＝j^t(p)j(p)+μdiag(j^t(p)j(p))，构造增量正规方程s(p)δp＝-j^t(p)r(p)，求解出δp，其中diag(j^t(p)j(p))表示由j^t(p)j(p)矩阵主对角元素组成的对角矩阵；

步骤7)：精度判断，若满足|δp|<ε则停止迭代、计算结束，此时p即为p*，跳转至步骤9)，否则跳转至步骤8)；

步骤8)：迭代判断，计算f(p+δp)和f(p)的值并比较大小，若f(p+δp)<f(p)，则令p＝δp+p更新p、令μ＝μ/ν减小阻尼因子，并跳转至步骤5)；若f(p+δp)≥f(p)，则令μ＝μν放大阻尼因子，并跳转至步骤6)；

步骤9)：levenberg-marquardt法迭代计算结束时，a*＝p*1、b*＝p*2，其中p*1表示向量p*的第1个元素、p*2表示向量p*的第2个元素。

优选地，在所述步骤2)中，在n个散点中，通过线性插值寻找中值电流拟合值a1*，基于中值电流拟合值a1*反算陡度参数b的拟合值b1*。

优选地，在所述步骤2)中，在n个散点中，若存在yk＝0.5则直接得出中值电流拟合值a1*＝xk；否则，在0.5附近找出最近的两点(xm,ym)和(xm+1,ym+1)，且满足ym>0.5>ym+1，过(xm,ym)、(xm+1,ym+1)两点的直线与y＝0.5交点横坐标即为a1*，由公式(6)计算得出，

；

然后基于中值电流拟合值a1*反算陡度参数b的拟合值b1*，对公式(2)作变换，按公式(15)计算b的拟合值b1*，

。

实施本发明一种雷电流幅值累积概率分布曲线拟合计算方法，具有以下有益效果：

1、在相同条件下，本发明的计算结果较线性变换拟合、中值电流法拟合等常见简单方法的拟合结果，残差平方和显著减小，约减少1～2个数量级，拟合曲线与统计散点最为接近，在常见拟合方法的结果中最优。

2、本发明使用中值电流法的结果作为初始值开始levenberg-marquardt迭代，初始值计算简单易实现，且避免了选择初始值的盲目性。

3、在相同条件下，本发明由于使用中值电流法的结果作为初始值开始levenberg-marquardt迭代，较选择ieee推荐值作初始值开始levenberg-marquardt迭代，减少了迭代步数，减少约26％的运算时间。

4、本发明使用levenberg-marquardt迭代过程中，迭代判断使用了简化的评价方法，仅根据迭代后目标函数值大小变化决定是否接受本轮迭代计算结果，这种评价方法简单易实现；

5、本发明以ieee推荐表达式为原型，保证计算结果的通用性。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明利用中值电流法求解其拟合值a1*的示意图；

图2为一种雷电流幅值累积概率分布曲线拟合计算方法的流程框图；

图3为采用本发明描述方法拟合2015年全国地闪雷电流幅值累积概率分布散点得到的曲线图；

图4为采用本发明描述方法和中值电流法、线性变换法对2005～2015年全国地闪雷电流幅值累积概率分布散点进行拟合的误差对比图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

如图2所示，本发明提供一种雷电流幅值累积概率分布曲线拟合计算方法，包括以下步骤：

步骤1)：确定原型函数如式(16)，确定最小二乘准则下拟合误差计算式如式(17)；其中，r为n个散点的拟合残差平方和，j为散点的序号，(xj,yj)为第j个散点，y(xj)为函数y(x)在xj处的值，

。

步骤2)：在n个散点中，通过线性插值寻找中值电流拟合值a1*，基于中值电流拟合值a1*反算陡度参数b的拟合值b1*。如图1所示，在n个散点中，若存在yk＝0.5则直接得出中值电流拟合值a1*＝xk；否则，在0.5附近找出最近的两点(xm,ym)和(xm+1,ym+1)，且满足ym>0.5>ym+1，过(xm,ym)、(xm+1,ym+1)两点的直线与y＝0.5交点横坐标即为a1*，由公式(18)计算得出，

；

然后基于中值电流拟合值a1*反算陡度参数b的拟合值b1*，对公式(16)作变换，按公式(19)计算b的拟合值b1*，

。

步骤3)：构造如公式(20)所示的levenberg-marquardt残差向量r(p)与如公式(21)所示的目标函数f(p)，

；

针对公式(16)的二元参数，对应拟合向量p＝[a,b]^t、元数m＝2；残差函数向量r(p)表征在n个统计散点处的残差组成的向量，每个元素rj(p)为m元函数；在最小二乘准则下，使f(p)值最小的向量p*即为最优解，残差平方和亦最小；

步骤4)：levenberg-marquardt迭代初始化，设置初始点p＝p0、向量计算精度ε、阻尼因子μ、阻尼因子变化倍率ν，其中，p0由中值电流法的结果代入，即p0＝[a1*,b1*]^t。向量计算精度ε、阻尼因子μ、阻尼因子变化倍率ν的取值分别为ε＝1.0×10^-5、μ＝0.001、ν＝10。

步骤5)：残差向量r(p)及jacobi矩阵j(p)迭代，代入p值计算残差向量r(p)、jacobi矩阵j(p)，jacobi矩阵j(p)如公式(22)所示，其中p1表示向量p第1个元素、p2表示向量p第2个元素，

。

步骤6)：求解增量向量，代入p值计算矩阵s(p)＝j^t(p)j(p)+μdiag(j^t(p)j(p))，构造增量正规方程s(p)δp＝-j^t(p)r(p)，求解出δp，其中diag(j^t(p)j(p))表示由j^t(p)j(p)矩阵主对角元素组成的对角矩阵。

步骤7)：精度判断，若满足|δp|<ε则停止迭代、计算结束，此时p即为p*，跳转至步骤9)，否则跳转至步骤8)。

步骤8)：迭代判断，计算f(p+δp)和f(p)的值并比较大小，若f(p+δp)<f(p)，则令p＝δp+p更新p、令μ＝μ/ν减小阻尼因子，并跳转至步骤5)；若f(p+δp)≥f(p)，则令μ＝μν放大阻尼因子，并跳转至步骤6)。

步骤9)：levenberg-marquardt法迭代计算结束时，a*＝p*1、b*＝p*2，其中p*1表示向量p*的第1个元素、p*2表示向量p*的第2个元素。

以2005～2015年全国范围内广域雷电监测系统探测的雷电地闪为例，使用本发明方法进行逐年雷电流幅值累积概率分布曲线拟合。在拟合之前，统计得到以5ka为间隔、雷电流幅值i取0、5、...、300ka处p(>i)的值，每年的雷电流幅值累积概率分布各61个散点(即n＝61)，并按雷电流幅值由小至大排序。本实施例参照本发明方法如图2所示的流程进行2015年全国地闪雷电流幅值累积概率分布散点的拟合计算，包括以下步骤：

步骤1)：散点个数n＝61，确定如公式(16)所示的原型函数，确定如公式(17)所示的拟合残差平方和作为最小二乘准则下拟合误差，用于评判拟合结果优劣。

步骤2)：在61个散点中，通过线性插值得出中值电流a1*＝23.52；基于中值电流a1*和公式(19)，计算陡度参数b的拟合值b1*＝2.81。

步骤3)：构造levenberg-marquardt残差向量与目标函数，针对公式(16)所示的二元参数，构造拟合向量p＝[a,b]^t、元数m＝2，如公式(20)所示的残差函数向量r(p)，如公式(21)所示的目标函数f(p)；在最小二乘准则下，使f(p)值最小的向量p*即为最优解，拟合误差r亦最小。

步骤4)：levenberg-marquardt迭代初始化，设置初始点p＝p0＝[23.52,2.81]^t、向量计算精度ε＝1.0×10^-5，阻尼因子μ＝0.001、阻尼因子变化倍率ν＝10。

步骤5)：残差函数向量r(p)及jacobi矩阵迭代，将p值代入公式(20)计算残差向量r(p)、代入公式(22)计算jacobi矩阵j(p)，其中p1表示向量p第1个元素、p2表示向量p第2个元素。

步骤6)：求解增量向量，代入p值计算矩阵s(p)＝j^t(p)j(p)+μdiag(j^t(p)j(p))，构造增量正规方程s(p)δp＝-j^t(p)r(p)，求解出δp，其中diag(j^t(p)j(p))表示由j^t(p)j(p)矩阵主对角元素组成的对角矩阵。

步骤7)：精度判断，若满足|δp|<ε则停止迭代、计算结束，此时p即为p*，跳转至步骤9)，否则跳转至步骤8)；

步骤8)：迭代判断，计算f(p+δp)和f(p)的值并比较大小，若f(p+δp)<f(p)，则令p＝δp+p更新p、令μ＝μ/ν减小阻尼因子，并跳转至步骤6；若f(p+δp)≥f(p)，则令μ＝μν放大阻尼因子，并跳转至步骤6)。

步骤9)：迭代计算结束，p*＝[23.11,2.42]^t，即a*＝23.11、b*＝2.42。

根据计算结果，可以绘制出2015年全国地闪雷电流幅值累积概率分布曲线，并与散点放在同一个坐标系中对比，如图3所示，同时可以按公式(17)算出拟合误差r＝0.0018。同理，可全部得出2005～2015年逐年全国地闪雷电流幅值累积概率分布结果，并计算出拟合误差。图3为采用本发明描述方法拟合2015年全国地闪雷电流幅值累积概率分布散点得到的曲线图；图4为采用本发明描述方法和中值电流法、线性变换法对2005～2015年全国地闪雷电流幅值累积概率分布散点进行拟合的误差对比图。从图3和图4可以看出，本发明方法的拟合曲线与统计散点极为贴近，拟合误差r最小，且较另外两种拟合计算方法的误差低1～3个数量级。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。