基于狼群算法的概率积分法模型参数反演方法研究

李靖宇，王 磊，朱尚军，滕超群，江克贵

(安徽理工大学测绘学院，安徽 淮南 232001)

我国近70%的天然能源来自于煤炭，近年来煤炭资源在能源消费中的占比逐渐降低，但短期内仍难以被取代。煤炭资源虽使用广泛、成本较低，但煤炭开采会造成开采区岩层发生移动、变形和破坏，为提前实施有效的预防措施，必须准确地预计开采沉陷的影响范围及变形值大小[1]。我国矿山开采沉陷预计的主要方法是概率积分法，由于模型简单、结果精准等优点，得到广泛应用[2]。提高开采沉陷预测精度[3]主要有两个途径：一是提高预测模型精度；二是提高预测参数的准确性与可靠性。在预测模型正确的情况下，概率积分法模型参数的精确性会直接影响到预计的结果[4]。因此，研究基于地表移动实测资料的概率积分参数稳健反演方法对开采沉陷精准预计具有重要意义。

概率积分法模型参数反演历程[5]中，经历了直接反演、实验设计方法反演、优化算法反演、智能算法反演的过程。智能算法是现阶段概率积分法参数反演的主要方法，查剑锋等[6]构建了基于遗传算法的概率积分法预计参数反演，抗干扰能力强，在求取参数方面具有一定的准确性、可靠性和稳定性，可以避免计算发散、陷入局部最优解；徐梦强等[7]提出采用粒子群算法反演概率积分法参数，搜索速度较快、效率高且算法简单；苏军明等[8]利用模拟退火算法反演概率积分法参数，精度较高，并且计算过程简单、通用；此后又有专家学者将果蝇算法[9]、量子遗传算法[10]引入到概率积分法参数反演中，取得了良好的效果。上述智能算法具有求参精度高、计算速度快、可用于非线性问题等良好的寻优特性。为获得更好的寻优特性，吴虎胜等[11]提出了狼群算法。狼群算法作为一种新智能算法，可以运用到复杂的、非线性的、多参数的寻优过程中，求取参数时能够保证较高准确性、可靠性和稳定性，同时能够以较大的速率找到最优解，可以从多个点出发进行搜索，在避免点与点之间的影响的同时，防止陷入局部最优解。目前，狼群算法被实际应用在求解多维背包问题、水电站水库优化调度以及旅行商等问题中[12-15]，在概率积分法参数反演中还未得到应用。

因此，本文拟将狼群算法引入到概率积分法参数反演中，构建基于WPA的概率积分法模型参数反演方法，通过模拟实验系统的讨论方法的可靠性，最后开展工程应用实验，并且讨论该算法在概率积分法参数反演中的可靠性和适用性。

1 基于WPA的概率积分法模型参数反演方法

1.1 概率积分法模型及其参数体系

如图1所示，考虑单元开采引起变形的三维情况，整个工作面开采引起坐标为(x,y)的地表任意点A的下沉值的计算过程见式(1)～(4)。

(1)

W0=mqcosα

(2)

W0(x)=

(3)

W0(y)=

(4)

地表任意点A沿φ方向的水平移动值，其计算公式见式(5)和式(6)。

U(x,y,φ)=bri(x,y,φ)

(5)

(6)

式中：φ为x轴的正向逆时针到指定方向的角值；r为主要影响半径；b为水平移动系数。

图1 单元开采地表任意点下沉预计原理图Fig.1 Schematic diagram of surface subsidenceprediction at any point of unit mining

综上可知，概率积分法模型的参数体系为B=[q,tanβ,b,θ,S1,S2,S3,S4]，参数体系B的解算精度直接决定开采沉陷的预计精度。

1.2 基于WPA的概率积分法模型参数反演方法

1.2.1 WPA理论基础

人工狼群狩猎区域是一个N×D的欧氏空间(N为人工狼总数，D为待求参数维度)，任意一匹人工狼的位置为Wn=(wn1,wn2,…,wnd,…,wnD)。在人工狼群中，每一匹人工狼都是随机生成的，随机生成人工狼群的公式见(7)和式(8)。

Wn=(wn1,wn2,…,wnd,…,wnD),

(1≤n≤N,1≤d≤D)

(7)

wnd=wmin+rand×(wmax-wmin)

(8)

式中：wnd为第n匹人工狼在第d维空间中的位置(即第n组随机参数中的第d个参数)；wmax、wmin分别为wnd中的最大值、最小值；rand为-1至1之间均匀分布的随机数。

头狼产生规则：在人工狼群中，头狼感知气味浓度最大，头狼不断变化且不参与游行、奔袭、围攻行为。

游走行为：探狼是除头狼外、感知气味浓度较大的人工狼。探狼每次游走向不同方向前进，每次前进后返回原来位置，若探狼前进后感知气味浓度超过头狼感知气味浓度，则更新头狼并进行奔袭行为，若不是，则继续游走直至更新头狼或达到最大游走次数；游走结束后，若探狼感知最大气味浓度大于当前位置感知气味浓度，则更新探狼。第k次前进后，探狼的位置计算见式(9)。

(9)

式中：wnd为开始搜寻猎物前，探狼在第d维空间中的位置；stepad为探狼在第d维空间中的游走步长。

奔袭行为：猛狼是距头狼较近的人工狼。猛狼在奔袭过程中，若感知气味浓度超过头狼感知气味浓度，则更新头狼且继续奔袭行为，直至猛狼与头狼之间的距离满足要求。猛狼奔袭后的位置见式(10)。

(10)

围攻行为：猛狼、探狼接收头狼的信号后，开始围攻猎物，若每次攻击后感知气味浓度超过原位置感知气味浓度，则更新该人工狼。第k+1次攻击后的位置见式(11)。

(11)

狼群更新规则：围攻结束后，淘汰感知气味浓度较小的人工狼，同时随机生成相同数量的人工狼。

1.2.2 参数反演方法构建

设开采区任意一点(x,y)的实测下沉值为Wxy实、水平移动值为Uxy实，将狼群算法反演参数带入概率积分法模型，得到的该点下沉值为Wxy预、水平移动值为Uxy预，以实测值与预计值之差平方和最小为准则，构建适应度函数为式(12)。

F=∑[(Wxy实-Wxy预)2+(Uxy实-Uxy预)2]

(12)

令概率积分法参数B=[q,tanβ,b,θ,S1,S2,S3,S4]为人工狼狩猎区域，以适应度函数值作为判断人工狼感知气味浓度大小的依据。在每一次迭代过程中，人工狼群在狩猎区域中感知气味浓度并对猎物(概率积分法参数)进行捕获，达到终止条件后，输出最优概率积分法参数。 求参主要流程如下所述。

1) 设置狼群数量N、最大迭代次数Dmax、最大游走次数Tmax等参数；输入观测点实测下沉值Wxy实、水平移动值Uxy实及概率积分法参数的波动范围ΔB=[Δq,Δtanβ,Δb,Δθ,ΔS1,ΔS2,ΔS3,ΔS4]，根据式(7)和式(8)生成人工狼群Bi=[qi,tanβi,bi,θi,S1i,S2i,S3i,S4i](i=1…N)。

2) 以适应度函数值为判断人工狼感知气味浓度Yn大小的依据，适应度函数值越小则人工狼感知气味浓度越大，选择感知气味浓度最大的人工狼为头狼，记其位置、感知气味浓度分别为wlead、Ylead。

3) 选择除头狼外、感知气味浓度较大的人工狼为探狼，探狼根据式(9)进行游走行为，直至达到最大游走次数Tmax或更新头狼。

4) 选择距离头狼较近的人工狼为猛狼，猛狼根据式(10)向头狼奔袭，若奔袭过程中猛狼感知气味浓度Yn大于Ylead，则更新头狼，继续召唤行为；若Yn小于Ylead，则猛狼继续奔袭，直至猛狼与头狼间距离Dn小于要求距离D。

5) 根据式(11)，探狼联合猛狼对猎物进行围攻。

6) 对头狼进行更新，淘汰部分感知气味浓度较小的人工狼，并根据式(7)和式(8)随机生成等数量人工狼。

7) 判断是否达到精度要求或最大迭代次数Dmax，若达到，则输出头狼位置，即概率积分法预计参数，否则转至第二步骤。

基于WPA的概率积分模型参数反演方法程序流程图如图2所示。

图2 WPA参数反演程序流程图Fig.2 Flow chart of WPA parameter inversion program

2 模拟实验

2.1 模拟实验工作面概况

以淮南矿区煤系地层为背景，设计模拟实验工作面参数为：煤层采厚m=3.0 m，煤层倾角α=5°，走向长度D3=800 m，倾向长度D1=500 m，平均采深H=400 m，工作面采用全部垮落法管理顶板。假设模拟工作面的概率积分法模型参数为：下沉系数q=0.8，主要影响正切角tanβ=2.5，水平移动系数b=0.25，开采影响传播角θ=85°，上拐点、下拐点、左拐点、右拐点偏移距分别为S1=60 m、S2=60 m、S3=60 m、S4=60 m。本次模拟实验在移动盆地内设计了沿走向、倾向的两条主断面的移动和变形的观测线，相邻两个监测点之间相距20 m，走向观测线1 180 m，共60个观测点，倾向观测线880 m，共45个观测点。模拟工作面设计图如图3所示。

2.2 狼群算法参数反演准确性

使用WPA进行概率积分参数反演，并进行模拟实验参数设计值与模拟实验参数反演值的对比，通过参数平均值、反演参数相对误差以及反演参数中误差来体现狼群算法反演参数的准确性。智能算法在一定范围内具有随机性，为避免偶然误差的影响，在相同的实验条件下连续进行了10次实验。设置人工狼数量为50，最大迭代次数为50，实验结果见表1。

由表1可知,在参数相对误差方面，q、tanβ、b、θ的参数相对误差均不超过0.3%，S1、S2、S3、S4对求参敏感性较低,但参数相对误差均小于2.1%；WPA反演概率比积分法模型参数具有较高的准确度。在参数中误差方面，q、tanβ、b、θ的参数中误差均不超过1.03，S1、S2、S3、S4的参数中误差均小于8.9 m；WPA反演概率比积分法模型参数具有较高的精度。综上所述，WPA反演概率积分法模型参数具有较好的稳健性。

图3 模拟工作面设计图Fig.3 Simulated working face design drawing

表1 WPA反演参数值与设计值对比Table 1 Comparison of WPA inversion parameter valueswith design values

为更加直观地看出下沉与水平移动的拟合情况，得到10次模拟实验反演参数均值计算出的下沉值曲线和水平移动值曲线见图4和图5。

由图4和图5可以看出，下沉值和水平移动值的实测曲线与拟合曲线基本一致，点位实测值与拟合值间的绝对值误差均不超过5 mm，达到预期精度。

2.3 基于WPA的概率积分法模型参数反演方法的可靠性研究

为进一步验证基于WPA算法的概率积分法模型参数反演方法的可靠性，对模拟实验10次反演结果中的概率积分法模型参数分别进行波动性分析，对比WPA反演参数与设计参数，结果见图6。

图4 模拟实验WPA拟合下沉值与实测下沉值对比图Fig.4 The comparison diagram of the simulated WPAfitted subsidence value and the measuredsubsidence value

图5 模拟实验WPA拟合水平移动值与实测水平移动值对比图Fig.5 The comparison diagram of the simulated WPAfitted horizontal movement value and the measuredhorizontal movement value

图6 模拟实验10次反演参数的波动情况Fig.6 The fluctuation of inversion parameters was simulated for 10 times

从图6可以看出，WPA反演参数波动范围较小，反演参数较为接近于真值，基于WPA算法的概率积分法模型参数反演方法具有良好的可靠性。

3 工程实例

3.1 矿区概况

淮南顾桥矿南二采区1414(1)工作面是顾桥煤矿南区的第一个采面，工作面从矿床边界开始向巷道往回开采，采用机械化掘进，垮落法管理顶板，一次采全高。工作面长约2 120 m，宽约251 m，走向为充分采动，倾向为非充分采动，总体为非充分采动。平均采高为3.0 m，煤层倾角在平均为5°，工作面平均深度为735 m，在距离切眼1 144 m处以及停采线976 m处布置倾向观测线，倾向线长度为1 500 m，相邻点之间点间距为30 m，共布设了3个控制点和50个监测点。 在下山方向偏离工作面中心线的距离39 m处设置走向观测线，长度为3 480 m，相邻点之间点间距为30 m/60 m，共布设了3个控制点和95个监测点。

3.2 实验结果及分析

以淮南顾桥矿南二采区1414(1)工作面的实测数据为基础，将WPA反演参数均值计算所得的下沉值、水平移动值与监测点的下沉值、水平移动值进行对比,为避免实验的偶然性，在相同条件下进行10次实验，然后分别计算反演参数平均值、中误差，结果见表2。

由表2可知，在参数中误差方面，q、tanβ、b、θ的参数中误差均不超过0.89，S1、S2、S3、S4的参数中误差均不超过4.5 m；WPA反演顾桥南矿1414(1)工作面的概率积分法参数为：q=0.935 7，tanβ=1.988 1，b=0.422 8，θ=84.138 3°，S1=-9.136 6 m，S2=-21.283 9 m，S3=55.954 4 m，S4=34.111 3 m。

将10次工程实例反演参数均值带入概率积分法模型，计算出的下沉值曲线和水平移动值曲线见图7和图8。

表2 WPA反演参数工程实例结果Table 2 Engineering example results of WPAinversion parameters

图7 工程实例WPA拟合下沉值与实测下沉值对比图Fig.7 The comparison diagram of WPA fitted subsidencevalue and the measured subsidence value

图8 工程实例WPA拟合水平移动值与实测水平移动值对比图Fig.8 The comparison diagram of WPA fitted horizontalmoving value and measured horizontal moving value

从图7和图8可以看出，下沉值和水平移动值的实测曲线与拟合曲线基本一致，点位实测值与拟合值间的绝对值误差均不超过300 mm。WPA反演概率积分法参数取均值后，其下沉与水平移动值的拟合中误差为114.81 mm，其中，下沉值拟合中误差约为最大下沉值的2.7%，水平移动值拟合中误差约为最大水平移动值的8.6%，满足工程应用精度[18]。

4 结 论

1) 模拟实验结果表明，WPA反演参数q、tanβ、b、θ的参数相对误差均不超过0.3%，参数中误差均不超过1.03；S1、S2、S3、S4对求参敏感性较低，但参数相对误差均小于2.1%，S1、S2、S3、S4的参数中误差均不超过8.9 m。WPA反演参数波动范围较小，反演参数较为接近于真值。以上证明了WPA在概率积分法模型参数反演中，具有良好的稳健性，且反演参数具有较高的准确性、精度。

2) 利用本文提出的基于WPA的概率积分法模型参数反演方法对淮南顾桥矿南二采区1414(1)工作面概率积分法参数反演，将由反演结果计算的下沉值、水平移动拟合值与真实值对比，获得的下沉值与水平移动拟合中误差为114.81 mm，下沉值拟合中误差约为最大下沉值的2.7%，水平移动值拟合中误差约为最大水平移动值的8.6%，满足工程应用精度。