CN104268624A - 一种基于单纯形萤火虫法反演大地电阻率测量数据的方法 - Google Patents

Abstract

一种基于单纯形萤火虫法反演大地电阻率测量数据的方法，包括：第一步，对每个萤火虫个体的位置和荧光素值进行初始化；第二步，更新萤火虫个体的荧光素值；第三步，萤火虫个体的邻域划分；第四步，确定第i只萤火虫向其邻居集合中第j只萤火虫移动的概率；第五步，利用轮盘赌法选择萤火虫个体j，然后对第i只萤火虫进行移动以更新萤火虫的位置；第六步：在N个邻域内的最优萤火虫按照当前位置为中心的搜索半径的变动范围进行单纯形法的局部搜索；第七步：若土壤反演的目标函数值已小于设定误差，则直接输出结果；若不满足，则最优萤火虫使用单纯形作局部搜索后，再返回第二步更新荧光素值。本发明可以增强算法的局部搜索能力和加快算法收敛速度。

Description

一种基于单纯形萤火虫法反演大地电阻率测量数据的方法

技术邻域

本发明涉及电力系统大地土壤电阻率反演的技术邻域，具体是一种基于单纯形萤火虫法反演大地电阻率测量数据的方法。

背景技术

目前电力系统广泛使用四极法作为大地电阻率的测量方法。由于现实中不存在均匀电阻率土壤的缘故，四极法测量取得视在电阻率是随极距变化的。由测量的视在电阻率推算土壤结构参数的过程称为反演。目前电力系统的土壤反演方法主要有最速下降法、最小二乘法、BFGS拟牛顿法和遗传算法，目前新型优化方法的应用尚不完善。

单纯形方法属于传统优化方法，其基本思想是比较一般单纯形的若干个顶点的目标函数的数值，并在迭代过程中逐步地把单纯形向最优点移动。人工萤火虫算法是最新的群智能优化方法之一，它利用萤火虫个体利用萤光素诱导其他萤火虫个体发光来吸引伴侣，各个萤火虫个体向荧光素值高的位置移动，再通过在动态决策域内寻找最高荧光素值的位置从而确定目标函数的最优解。

单纯形法虽然收敛较快，但需要指定初值，有时反演的结果与初值设置有关；人工萤火虫法虽然无需指定初值，但计算结果较慢，容易收敛于局部解。

发明内容

本发明提供一种基于单纯形萤火虫法反演大地电阻率测量数据的方法，其克服了经典优化方法存在反演结果只是初值范围的局部解问题，又弥补了人工萤火虫算法局部搜索能力不足的缺点，进而增强算法的局部搜索能力和加快算法收敛速度。

一种基于单纯形萤火虫法反演大地电阻率测量数据的方法，包括如下步骤：

第一步：设在n层水平多层土壤反演的2n-1维搜索空间内有p(p＝20n-20)个萤火虫，其中n为大于等于2的自然数，对每个萤火虫个体的位置和荧光素值进行初始化；

第二步：使用公式(1)把第i个萤火虫在第t次迭代的位置Θ_i(t)对应的目标函数值f_RMS-error(Θ_i(t))转化为荧光素值；

l_i(t)＝(1-η)l_i(t-1)+γf_RMS-error(Θ_i(t))   (1)

公式(1)中，t代表迭代次数，li(t)为第t代第i个萤火虫的萤光素值，η是荧光素的挥发因子，η∈(0，l]，γ是荧光素更新率；

第三步：从第二代开始，萤火虫算法根据各个萤火虫荧光素值的相近程度将该群体分成N个邻域： $N_i\left(t\right)=\{j:||x_j\left(t\right)-x_i\left(t\right)||<r_d^i(t+1);l_i\left(t\right)$ $<l_j\left(t\right)\},$ l_j(t)代表当前萤火虫j的荧光素值，l_i(t)代表当前萤火虫i的荧光素值，N_i(t)表示t时刻第i个邻域的密度；x_j(t)和x_i(t)分别表示萤火虫j和萤火虫i第t次迭代求解过程中对应于求解空间中的位置，||||代表向量范数，为局部决策域半径，表达式如下：

$r_d^i(t+1)=\frac{r_s}{1+{\mathrm{\&mu;N}}_i\left(t\right)}---\left(3\right)$

其中，rs是萤火虫感知半径；μ是动态域决策更新率；

第四步：根据公式(4)确定第t次迭代第i只萤火虫向其邻居集合中第j只萤火虫移动的概率P_ij(t)：

$P_{\mathrm{ij}}\left(t\right)=\frac{l_j\left(t\right)-l_i\left(t\right)}{\underset{k=1}{\overset{N_i\left(t\right)}{\mathrm{\&Sigma;}}}{l}_k\left(t\right)-l_i\left(t\right)}---\left(4\right)$

第五步：利用轮盘赌法选择萤火虫个体j，然后根据公式(5)对第i只萤火虫进行移动以更新萤火虫的位置：

${\mathrm{\&Theta;}}_i(t+1)={\mathrm{\&Theta;}}_i\left(t\right)+s\left(\frac{{\mathrm{\&Theta;}}_j\left(t\right)-{\mathrm{\&Theta;}}_i\left(t\right)}{\left|\right|{\mathrm{\&Theta;}}_j\left(t\right)-{\mathrm{\&Theta;}}_i\left(t\right)\left|\right|}\right)---\left(5\right)$

公式(5)中s为步长；

第六步：在N个邻域内的最优萤火虫按照当前位置为中心的搜索半径w的变动范围进行单纯形法的局部搜索，

$w\left(t\right)=0.1[1.1-\frac{K\left(t\right)}{K_{\max }}]---\left(6\right)$

公式(6)中，K(t)为当前迭代次数，K_max为设定的最大迭代次数；

第七步：若土壤反演的目标函数值已小于设定误差，则直接输出结果；若不满足，则最优萤火虫使用单纯形作局部搜索后，再按返回第二步更新荧光素值。

进一步的，第二步中利用四极法测量浅层大地电阻率时，土壤反演的目标函数f_RMS-error用视在电阻率测量值与视在电阻率估算值的均方根误差来表示，其计算式为公式(2)：

$f_{\mathrm{RMS}-\mathrm{error}}\left(\mathrm{\&Theta;}\right)=\min \sqrt{\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left(\frac{{\mathrm{\&rho;}}_{\mathrm{ai}}-{\mathrm{\&rho;}}_{\mathrm{Mi}}}{{\mathrm{\&rho;}}_{\mathrm{Mi}}}\right)}^2}{m}}$

Θ＝(ρ₁,…,ρ_n,h₁,…,h_n-1)   (2)

公式(2)中，m为四极法的极距数量，ρ_Mi为第i组极距对应的视在电阻率测量值，ρ_ai为第i组极距对应的视在电阻率计算值，ρ_ai的结果与Θ有关，Θ是n层土壤模型结构，h_i(i＝1,2,...,n-1)和ρ_i(i＝1,2,...,n)分别代表第i层的厚度和电阻率，其结果由优化方法调整取得。

本发明的能准确有效的对四极法测量数据进行反演，准确取得大地电阻率分布，避免造成变电站和接地极的接地电阻的设计值与实际值相差较大，有利于准确有效的掌握变电站及接地极附近的跨步电势差和接触电势差分布。

附图说明

图1是利用单纯形法反演测量结果示意图；

图2是利用利用萤火虫法反演测量结果示意图；

图3是利用本发明单纯形萤火虫法反演测量结果示意图。

具体实施方式

下面将结合本发明中的附图，对本发明中的技术方案进行清楚、完整地描述。

本发明提供一种基于单纯形萤火虫法反演大地电阻率测量数据的方法，包括如下步骤：

第一步：设在n层水平多层土壤反演的2n-1维搜索空间内有p(推荐取p＝20n-20)个萤火虫，其中n为大于等于2的自然数，对每个萤火虫个体先进行初始化，即随机指定萤火虫的位置，并把每个萤火虫的荧光素值均指定为某一常数l₀。

第二步：使用公式(1)把第i个萤火虫在第t次迭代的位置Θ_i(t)对应的目标函数值f_RMS-error(Θ_i(t))转化为荧光素值；

l_i(t)＝(1-η)l_i(t-1)+γf_RMS-error(Θ_i(t))   (1)

公式(1)中，t代表迭代次数，li(t)为第t代第i个萤火虫的萤光素值；η是荧光素的挥发因子，η∈(0，l]；γ是荧光素更新率，本实施例取η＝0.5，γ＝0.2。利用四极法测量浅层大地电阻率时，土壤反演的目标函数f_RMS-error可用视在电阻率测量值与视在电阻率估算值的均方根误差，其计算式为公式(2)：

$f_{\mathrm{RMS}-\mathrm{error}}\left(\mathrm{\&Theta;}\right)=\min \sqrt{\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left(\frac{{\mathrm{\&rho;}}_{\mathrm{ai}}-{\mathrm{\&rho;}}_{\mathrm{Mi}}}{{\mathrm{\&rho;}}_{\mathrm{Mi}}}\right)}^2}{m}}$

Θ＝(ρ₁,…,ρ_n,h₁,…,h_n-1)   (2)

公式(2)中，m为四极法的极距数量，ρ_Mi为第i组极距对应的视在电阻率测量值(Ω·m)，ρ_ai为第i组极距对应的视在电阻率计算值(Ω·m)，ρ_ai的结果与Θ有关，Θ是n层土壤模型结构，h_i(i＝1,2,...,n-1)和ρ_i(i＝1,2,...,n)分别代表第i层的厚度和电阻率，其结果由优化方法调整取得。

第三步：从第二代开始，萤火虫算法根据各个萤火虫荧光素值的相近程度将该群体分成N个邻域： $N_i\left(t\right)=\{j:||x_j\left(t\right)-x_i\left(t\right)||<r_d^i(t+1);l_i\left(t\right)$ $<l_j\left(t\right)\},$ l_j(t)代表当前萤火虫j的荧光素值，l_i(t)代表当前萤火虫i的荧光素值，N_i(t)表示t时刻第i个邻域的密度；x_j(t)和x_i(t)分别表示萤火虫j和萤火虫i第t次迭代求解过程中对应于求解空间中的位置，||||代表向量范数，本实施例推荐使用欧几里德范数；为局部决策域半径，表达式如下：

$r_d^i(t+1)=\frac{r_s}{1+{\mathrm{\&mu;N}}_i\left(t\right)}---\left(3\right)$

其中，r_s是萤火虫感知半径，推荐取0.1；μ是动态域决策更新率，推荐取1.0。公式(3)表明：当萤火虫j的荧光素值大于萤火虫i的荧光素值，且萤火虫j与萤火虫i之间的距离小于萤火虫i所在邻域的决策域半径时，将萤火虫j划分到萤火虫i所在的邻域。

第四步：萤火虫在运动过程中，根据其邻居集合中各萤火虫的荧光素浓度来决定其移动方向，根据公式(4)确定t时刻第i只萤火虫向其邻居集合中第j只萤火虫移动的概率P_ij(t)：

$P_{\mathrm{ij}}\left(t\right)=\frac{l_j\left(t\right)-l_i\left(t\right)}{\underset{k=1}{\overset{N_i\left(t\right)}{\mathrm{\&Sigma;}}}{l}_k\left(t\right)-l_i\left(t\right)}---\left(4\right)$

公式(4)中，第i只萤火虫移动概率P_ij(t)根据轮盘赌法选择其移动的方向。

第五步：利用轮盘赌法选择萤火虫个体j，然后根据公式(5)对第i只萤火虫进行移动以更新萤火虫的位置：

${\mathrm{\&Theta;}}_i(t+1)={\mathrm{\&Theta;}}_i\left(t\right)+s\left(\frac{{\mathrm{\&Theta;}}_j\left(t\right)-{\mathrm{\&Theta;}}_i\left(t\right)}{\left|\right|{\mathrm{\&Theta;}}_j\left(t\right)-{\mathrm{\&Theta;}}_i\left(t\right)\left|\right|}\right)---\left(5\right)$

公式(5)中s为步长，本实施例推荐取0.1。

第六步：在N个邻域内的最优萤火虫按照当前位置为中心的搜索半径w的变动范围进行单纯形法的局部搜索。

$w\left(t\right)=0.1[1.1-\frac{K\left(t\right)}{K_{\max }}]---\left(6\right)$

公式(6)中，K(t)为当前迭代次数，K_max为设定的最大迭代次数。

单纯形方法的基本思想是比较一般单纯形的k+1个顶点的目标函数值，并在迭代过程中逐步地把单纯形向最优点移动。单纯形的移动是通过反射、收缩和扩大的三种运算来实现，这样就通过单纯形加强了萤火虫算法的局部搜索能力。公式(6)表明，在人工萤火虫迭代计算的初始阶段，单纯形局部搜索的范围较大以提高全局搜索能力，在迭代计算的最后阶段，搜索范围收窄以提高局部搜索的能力。

在传统萤火虫法中，动态决策域半径更新按公式(7)进行。

$r_d^t(t+1)=\min \{r_s,\max \{0,r_d^t\left(t\right)+\mathrm{\&beta;}(n_i-N_i\left(t\right))\left\}\right\}---\left(7\right)$

每只萤火虫在动态决策域半径内，选择荧光素值比自己高的个体组成邻域集N_i(t)，其中0<r_id(t)≤r_s，r_s为萤火虫个体的感知半径。

相较现有技术，本发明可以加快公式(7)中动态决策域半径动态决策域半径的更新过程。

第七步：若土壤反演的目标函数值已小于设定误差，例如公式(2)中视在电阻率测量值与视在电阻率估算值的均方根误差f_RMS-error已小于设定误差，则直接输出结果；若不满足，则最优萤火虫使用单纯形作局部搜索后，再按返回第二步更新荧光素值。

利用四极法对大地电阻率进行测量的测量结果如表1所示(a为电极间隔距离,ρ_Mi为极距对应的视在电阻率测量值电阻率)：

表1 四极法测量结果

分别利用三种算法进行反演，利用单纯形法反演如图1所示，相对误差百分数为7.2％，利用萤火虫法反演如图2所示，相对误差百分数为9.1％，利用本发明单纯形萤火虫法反演如图3所示，相对误差百分数为1.38％。由此可见，利用单纯形萤火虫法能够有效提高反演精度。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何属于本技术邻域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (2)

1.一种基于单纯形萤火虫法反演大地电阻率测量数据的方法，其特征在于包括如下步骤：

第一步：设在n层水平多层土壤反演的2n-1维搜索空间内有p(p＝20n-20)个萤火虫，其中n为大于等于2的自然数，对每个萤火虫个体的位置和荧光素值进行初始化；

第二步：使用公式(1)把第i个萤火虫在第t次迭代的位置Θ_i(t)对应的目标函数值f_RMS-error(Θ_i(t))转化为荧光素值；

l_i(t)＝(1-η)l_i(t-1)+γf_RMS-error(Θ_i(t))   (1)

公式(1)中，t代表迭代次数，li(t)为第t代第i个萤火虫的萤光素值，η是荧光素的挥发因子，η∈(0，l]，γ是荧光素更新率；

第三步：从第二代开始，萤火虫算法根据各个萤火虫荧光素值的相近程度将该群体分成N个邻域： $N_i\left(t\right)=\{j:||x_j\left(t\right)-x_i\left(t\right)||<r_d^i(t+1);l_i\left(t\right)$ $<l_j\left(t\right)\},$ l_j(t)代表当前萤火虫j的荧光素值，l_i(t)代表当前萤火虫i的荧光素值，N_i(t)表示t时刻第i个邻域的密度；x_j(t)和x_i(t)分别表示萤火虫j和萤火虫i第t次迭代求解过程中对应于求解空间中的位置，||||代表向量范数，为局部决策域半径，表达式如下：

$r_d^i(t+1)=\frac{r_s}{1+{\mathrm{\&mu;N}}_i\left(t\right)}---\left(3\right)$

其中，rs是萤火虫感知半径；μ是动态域决策更新率；

第四步：根据公式(4)确定第t次迭代第i只萤火虫向其邻居集合中第j只萤火虫移动的概率P_ij(t)：

$P_{\mathrm{ij}}\left(t\right)=\frac{l_j\left(t\right)-l_i\left(t\right)}{\underset{k=1}{\overset{N_i\left(t\right)}{\mathrm{\&Sigma;}}}{l}_k\left(t\right)-l_i\left(t\right)}---\left(4\right)$

第五步：利用轮盘赌法选择萤火虫个体j，然后根据公式(5)对第i只萤火虫进行移动以更新萤火虫的位置：

${\mathrm{\&Theta;}}_i(t+1)={\mathrm{\&Theta;}}_i\left(t\right)+s\left(\frac{{\mathrm{\&Theta;}}_j\left(t\right)-{\mathrm{\&Theta;}}_i\left(t\right)}{\left|\right|{\mathrm{\&Theta;}}_j\left(t\right)-{\mathrm{\&Theta;}}_i\left(t\right)\left|\right|}\right)---\left(5\right)$

公式(5)中s为步长；

第六步：在N个邻域内的最优萤火虫按照当前位置为中心的搜索半径w的变动范围进行单纯形法的局部搜索，

$w\left(t\right)=0.1[1.1-\frac{K\left(t\right)}{K_{\max }}]---\left(6\right)$

公式(6)中，K(t)为当前迭代次数，K_max为设定的最大迭代次数；

第七步：若土壤反演的目标函数值已小于设定误差，则直接输出结果；若不满足，则最优萤火虫使用单纯形作局部搜索后，再按返回第二步更新荧光素值。

2.如权利要求1所述的基于单纯形萤火虫法反演大地电阻率测量数据的方法，其特征在于：第二步中利用四极法测量浅层大地电阻率时，土壤反演的目标函数f_RMS-error用视在电阻率测量值与视在电阻率估算值的均方根误差来表示，其计算式为公式(2)：

$f_{\mathrm{RMS}-\mathrm{error}}\left(\mathrm{\&Theta;}\right)=\min \sqrt{\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left(\frac{{\mathrm{\&rho;}}_{\mathrm{ai}}-{\mathrm{\&rho;}}_{\mathrm{Mi}}}{{\mathrm{\&rho;}}_{\mathrm{Mi}}}\right)}^2}{m}}$

Θ＝(ρ₁,…,ρ_n,h₁,…,h_n-1)   (2)

公式(2)中，m为四极法的极距数量，ρ_Mi为第i组极距对应的视在电阻率测量值，ρ_ai为第i组极距对应的视在电阻率计算值，ρ_ai的结果与Θ有关，Θ是n层土壤模型结构，h_i(i＝1,2,...,n-1)和ρ_i(i＝1,2,...,n)分别代表第i层的厚度和电阻率，其结果由优化方法调整取得。