基于稳健遗传算法的矿山开采沉陷预计参数反演

杨晓玉 朱晓峻

(1.合肥财经职业学院建筑工程学院,安徽 合肥 230601;2.安徽大学资源与环境工程学院,安徽 合肥 230601;3.中国矿业大学环境与测绘学院,江苏 徐州 221116)

随着我国经济的飞速发展,煤炭资源的产量也逐年增加。大规模的煤炭开采也带来了一系列严重的社会和环境问题[1],如地面沉降、建筑物破坏、地裂缝、重金属污染等[2-4],已引发广泛关注。地下煤炭资源开采引起的地表沉陷是矿区最常见的灾害之一,直接影响地表建筑物和基础设施,可能危及生命和财产损失。统计数据显示,中国开采沉陷总面积约为6×103km2,并且每年以约240 km2的速度递增[5]。矿区沉陷问题严重影响着矿区的可持续发展。为此,为实现资源开发与生态保护的和谐共存,提出了条带开采、充填开采、覆岩离层注浆等地表减沉技术,地表沉陷的准确预计是决定这些方法成功与否的一个因素。因此,准确预计地下开采引起的地表沉陷具有重要的现实意义。

概率积分法是目前国内应用最广泛的地表沉陷预计方法。概率积分法的应用在减少矿区地表沉陷灾害和经济损失方面发挥了积极作用。准确的开采沉陷预计应以准确的参数为基础,但是由于上覆岩层的复杂性,从上覆岩层的物理性质中很难获得可靠的概率积分法参数,目前主要通过实测地表沉降反演获得。目前,我国已经建立了大量的地表移动特征的观测站,对推进开采沉陷研究发挥了重要作用。

概率积分法模型是一个复杂的非线性模型。此外,上覆岩层的复杂性、概率积分法的模型误差、现场监测误差等因素导致概率积分法参数反演结果的可靠性和稳定性较差[6]。因此,有学者对参数反演方法进行了研究。目前常用的概率积分法参数反演方法有特征点法、最小二乘曲线拟合法[7]、模矢搜索法、遗传算法[8]、粒子群算法[9-10]等。特征点法是利用地表下沉盆地特征点计算地表下沉参数的一种近似方法。当测量点不是这些特征点时,该方法误差较大。最小二乘曲线拟合方法对初始参数和粗差非常敏感,这往往会导致反演结果失真,偏离参数的实际物理范围。郭广礼、汪云甲[11-12]提出了一种基于稳健估计和最小二乘曲线拟合的测量数据粗差参数反演方法,在一定程度上提高了参数反演的精度。吴侃等[13]使用模式搜索方法来获得准确可靠的参数,但该方法对初始参数的要求也很高,不恰当的初始参数往往会导致结果陷入局部最小陷阱,无法获得全局最优参数。查剑锋等[14]利用具有良好全局搜索能力的遗传算法解决初始参数问题,获得准确可靠的参数。但遗传算法参数反演过程中的适应度计算仍采用最小二乘法判断,难免受到实测数据误差较大的影响,降低参数反演精度。

为此,本文提出了一种稳健估计方法与遗传算法的融合方法,对概率积分法预计参数进行反演,解决了概率积分法模型误差和现场实测数据粗差造成的参数反演困难,使得到的参数更加稳健可靠。

1 概率积分法模型与稳健遗传反演算法

1.1 概率积分法模型

概率积分法是一种基于随机介质理论的开采沉陷预计方法[15]。根据概率积分法开采沉陷预计的原理(如图1),微小单元引起的地表沉陷计算公式如下:

图1 概率积分法模型计算的示意Fig.1 A diagram of surfer movement calculated by PIM

式中,We(x,y) 为小单元开采引起的地表沉陷;(x,y)为地表点坐标;r为主要影响半径,r=H/tanβ;H为开采深度;tanβ为主要影响角的正切。

对整个工作面进行积分,工作面开采引起的任意一点下沉值可计算为

式中,W(x,y) 为工作面开采沉陷量;W0为最大地面下沉值,W0=mqcosα;m为开采厚度;q为下沉系数;α为煤层倾角。

同样,任意位置的倾斜、曲率、水平移动和水平变形均可由二重积分方程计算获得:

式中,φ为地表点的移动变形方向,从x轴的正向逆时针计算到指定方向的角值。

1.2 基于遗传算法参数反演的基本步骤

遗传算法(简称GA)是一种借鉴生物界自然选择和自然遗传机制的高度并行、随机、自适应搜索算法,它主要用于处理最优化问题和机器学习。遗传算法对优化函数的限制很少,可以处理非常特殊的函数。同时对搜索空间的若干个点进行搜索,可以有效地利用全局信息,防止收敛于局部最优解。

基于遗传算法参数反演概率积分法预计参数的基本步骤如下:

步骤1:初始化遗传算法运行参数,设置种群大小、决策变量个数、编码方式、长度及进化终止准则等信息。

步骤2:随机生成初始种群。

步骤3:计算种群中各个个体适应度,个体适应度大小与个体被遗传到下一代的机会大小成正比。适应度必须大于或等于零。对于不同的问题,需要确定好由目标函数值到适应度之间的转化规则。

步骤4:执行选择运算,选择适应度较高的个体构成下一代种群;选择运算建立在对个体适应度的评价之上,针对具体的问题可以采用不同的选择算子。

步骤5:对选出的下一代个体进行交叉运算,交叉运算模拟自然界生物的基因重组现象,从配对后的2 个父代个体产生出新的子代个体,它是产生新个体的主要方法,决定了遗传算法的全局搜索能力。

步骤6:对选出的下一代个体进行变异运算,变异运算改变个体编码串中的某些基因值,从而生成新的个体,它是产生新个体的辅助方法,决定了遗传算法的局部搜索能力。

步骤7:终止条件判断,若当前步满足步骤1 中设置的终止准则,则算法终止,将进化过程中得到的具有最大适应度的个体作为最优解输出;否则,转到步骤3。

1.3 稳健估计方法

稳健估计(Robust Estimation)又称抗差估计,是一种极大似然数估计中的一种特殊方法,具备抗粗差的能力,保证最后的结果不受或较小受粗差的影响。稳健估计是在极大似然估计的基础上建立起来的,都是通过一定的原则影响对数据的选择,从而达到消除异值,保证数据稳健化。

稳健估计的方法有很多种,本文采用选权迭代法结合遗传算法进行反演参数,其模型如下:开始假设每个观测值都是等精度等权的即Pi=1,经过一次的遗传算法求参后,求出预计值与实测值的差值,通过对残差大小的判断,从而改变每个观测值的权重,则在原有的残差平方和的基础上加入权重的影响:

VTPV= min ,

本研究选用Huber 法和残差一次范数最小法进行试验。

(1)Huber 法的函数:

式中,为第k+1 次迭代时V所用的权;U=Vk/S,S为V的最大似然值即S=Median(V),Vk为第k次迭代后的第k个观测值的残差;C为很小的数。

(2)一次范数法也叫残差绝对最小法,其权函数:

1.4 稳健遗传算法

为了解决利用遗传算法求解开采沉陷概率积分法参数时抗粗差能力差的问题,本研究尝试将稳健估计方法与遗传算法相结合,具体步骤(如图2)如下:

图2 稳健遗传算法的基本步骤Fig.2 Flow chart of PIM parameters inversion method based on GA and robust estimation

首先根据已知地质采矿条件,准备开采工作面平均采深H、煤层采厚m、煤层倾角α、工作面走向长ZC、倾向长QC等参数,然后给定的概率积分法预计参数(下沉系数q、主要影响正切值tanβ、开采影响传播角θ、拐点偏移距S、水平移动系数b)的范围,生成初始种群,将这些生成的种群解码成概率积分法参数,预计出每个观测点的预计下沉值,再将预计值与实测值相比较,求出其残差平方和,根据残差平方和采用稳健估计的方法确定每个测点的权重,残差平方和越大,其适应度越低,即相互占的权重越低,最后根据稳健估计后的权重重新将种群进行选择操作、交叉操作和变异操作,直到预计值与实测值十分接近才停止循环。这样便通过稳健估计与遗传算法相结合获得一组最优的概率积分法参数,实测值中异常值或粗差通过稳健估计降权之后,降低对遗传算法最终结果的影响,保证获得的参数稳定准确。图中f(x) =Cmax-g(x) 为适应度函数,用以判断每次迭代个体的适应度,Cmax为迭代过程中取得的最大阈值。

2 基于稳健遗传算法的开采沉陷预计参数反演方法精度验证

2.1 工程实例

五沟煤矿位于中国安徽省淮北市,选取五沟煤矿1013 工作面作为研究对象进行分析。1013 工作面走向长度575 m,倾斜宽度150 m,采区面积约86 300 m2,平均采深388 m,平均倾斜角10°,平均采厚3.1 m。工作面上覆岩层以砂岩、泥岩为主,厚度约94 m,松散层约270 m。工作面表面平整,地表标高26.50～27.56 m。采用综合机械化长壁垮落采煤法开采。该工作面于2008 年5 月开始回采,至2009 年4 月回采结束。为了获得该工作面地表沉陷观测数据,在五沟煤矿1013 工作面上方设置了一条地表移动变形观测站,地表移动变形观测站示意如图3 所示,地面沉降监测采用四阶水准测量,每公里高差均方误差不大于10 mm[16]。

图3 地表移动变形观测站示意Fig.3 Observation point layout above 1013 working face

2.2 基于实测数据的精度验证

为了证明本文参数反演方法具有抗粗差的能力,在实际观测数据中加入了若干粗差。例如,N7 和N32 点的下沉值增加200 mm,N24 点的下沉值减少400 mm。然后利用这些下沉数据反演概率积分法参数,利用反演参数预计观测站的下沉值。最后,将预计下沉值与实测下沉值进行对比,分析参数反演方法的精度。图4 为N 线观测站在上述工作面存在粗差时的最终下沉曲线。从图4 可以看出,整个下沉曲线符合一般下沉特征。但下沉盆地边缘点N7 和N32的下沉值明显大于周边点,下沉盆地中心点N24 的下沉值明显小于周边监测点。这些数据是最终数据中的较大误差,会影响概率积分法的最终反演结果。

图4 地表移动下沉曲线Fig.4 Curve of measured values

为了分析这些粗差对不同方法反演结果的影响,本文采用3 种方法分别反演概率积分法的预计参数。第一种方法是采用遗传算法对无粗差数据的参数进行反演(从下沉数据中去除N7、N24、N32 点的下沉值);第二种方法是采用遗传算法对有粗差数据的参数进行反演;第三种方法是采用本文的方法(稳健遗传算法)对有粗差数据的参数进行反演。在实验结果的基础上,分析了遗传算法与稳健估计方法的融合方法抗粗差的能力。

3 种方法的实测值与预计值的对比曲线如图5所示。由于第一种方法的数据没有明显误差,可以认为概率积分法反演参数是准确的参数值,并将后一种方法得到的反演参数与第一种方法得到的反演参数进行比较分析。从图5 可以看出,3 种方法得到的预计下沉曲线与实测下沉值的变化趋势一致,即工作面中心下沉值较大,盆地边界下沉值较小。本文方法(一种有粗差的遗传算法和稳健估计方法的融合方法)与第一种方法(无粗差的遗传算法)计算的下沉曲线拟合结果基本一致,说明本文方法避免了粗差对反演结果的影响,反演参数准确。而第二种方法(误差较大的遗传算法)计算的下沉曲线与第一种方法拟合较差,下沉曲线明显左移。在N24 点处的粗差的影响导致拟合曲线移位。下沉曲线难以准确表示下沉中心的位置,可能会对后续灾害评估造成错误的指导。因此,遗传算法与稳健估计方法的融合方法对粗差具有抗干扰性,且对粗差的抗干扰性明显优于遗传算法。

图5 3 种方法实测值与预计值对比分析Fig.5 Contrast analysis diagram of the measured values and prediction values by three methods

为了进一步分析本文方法抗粗差的能力,将3 种方法得到的概率积分法反演参数及均方根误差进行比较,如表1 所示。同样,由于第一种方法的数据没有明显误差,因此将第一种方法得到的概率积分法反演参数视为准确的参数值,并将后2 种方法得到的反演参数与第一种方法得到的反演参数进行比较分析。由表1 可以看出,第二种方法得到的概率积分法反演参数与第一种方法相差较大,q的相对误差为2.73%,tanβ的相对误差为2.83%,θ的相对误差为0.62%,S的相对误差为50.33%。本文方法得到的概率积分法反演参数与第一种方法相差不大,q的相对误差为0.84%,tanβ的相对误差为1.61%,θ的相对误差为0.36%,S的相对误差为8.6%。因此,本文方法得到的参数更接近于精确的参数值,且没有较大的误差。

表1 3 种方法反演的概率积分法参数对比Table 1 The differences between PIM parameters obtained by three methods

同时,对3 种方法的均方根误差进行了分析。均方根误差(RMSE)的公式为

式中,n为观测点数;WMeasuredi为点i的实测下沉值;WInversedi为反演方法预计的点i的预计下沉值。

在表1 中,不存在粗差的遗传算法的RMSE为55 mm,存在粗差的遗传算法的RMSE为87 mm,本文方法(存在粗差的稳健遗传算法)的RMSE为58 mm。本文方法的RMSE明显小于第二种方法,更接近第一种方法。

从上述分析可以看出,遗传算法与稳健估计方法的融合方法避免了粗差对反演结果的影响,反演参数更加准确可靠。该方法解决了现场实测数据误差较大导致的参数反演困难的问题。

2.3 模拟实验分析

为了验证稳健遗传算法在沉陷参数反演中效果,设计如下模拟试验:

概率积分法预计参数大多根据地表移动监测结果反演得出,由于实测地表下沉盆地与概率积分法模型描述的下沉盆地存在一定的出入,根据实测数据反演概率积分法参数的准确性受到反演算法或模型误差的双重影响,很难根据实测数据反演结果的好坏判断出是模型误差或算法准确性的影响程度。因此,为了在验证算法准确性时不考虑模型误差的影响,实验方案设置为:首先根据预先设定的概率积分法参数(以此为真值),预计出某一具体工作面上方一些离散点的下沉和水平移动值;然后以此预计值为起算数据,反演其对应的概率积分法参数,计算得出相应的反演值,比较反演值和真值之间的差异来判断算法的准确性。

模拟试验区的地质采矿条件为:煤层采厚3.0 m,倾角12°,工作面倾向长D1=400 m,走向长D3=600 m,平均采深H=300 m,采用全部垮落法管理顶板。预计点间隔75 m,预计的离散点与工作面的相对位置关系见图6。

图6 模拟工作面示意Fig.6 Schematic of simulated working face

为了验证算法对于数据中的粗差是否具有抗干扰能力,在工作面不同位置对实验数据人为加入一些异值点,其实验如下:分别在下沉盆地边界、工作面拐点处、最大下沉点处设置粗差(预计值增加200 mm),分别利用不采用稳健估计的遗传算法和稳健遗传算法反演含有粗差数据的概率积分法参数,并与设置粗差前的概率积分法参数进行比较,结果见表2。参数反演值预计的拟合效果见图7。

表2 不同位置的粗差对比Table 2 Comparison of gross errors at different positions

图7 不同位置的粗差对拟合效果的影响Fig.7 Influence of different position gross error on fitting effect

从表2 和图7 对比看出,当异值点在工作面边缘时,无论采用单纯的遗传算法还是结合稳健估计的方法,求参的效果都比较好。而且在拟合图中可以看出,拟合线在异值点处拟合明显不好,趋于理论的值。但当异值点出现在拐点和最大下沉点附近时,单纯采用遗传算法求参的效果都受一定影响,与理论值有一定偏差,采用稳健遗传算法求参的效果较好,说明测量值中存在粗差时,通过稳健遗传算法可以有效避免粗差对反演结果的影响。

从表2 中可以看出,在人为加入异值点后,当异值点出现在下沉盆地边缘时,单纯采用遗传算法求参的效果较好,而且在拟合图中可以看出,拟合线在异值点处拟合明显不好,趋于理论的值,说明下沉盆地附近出现的粗差对于参数求解结果的精度影响较小,因此在平时实际移动变形观测时,在盆地边缘的点测量的要求可适当放宽。但当粗差出现在工作面拐点处和最大下沉点处时,其反演结果偏离理论值比较大、精度较低,尤其是tanβ的值受到很大的影响。一方面说明当实测值中出现异值时,采用遗传算法进行概率积分法参数反演效果较差,拟合线受异值影响偏离理论值;另一方面说明当异值在下沉盆地边缘时,对反演结果影响较小,当异值出现在工作面拐点处和下沉盆地中心时,对反演结果影响较大。

采用Huber 法与遗传算法进行结合时,反演出来的结果与单纯采用遗传算法反演的结果类似,说明采用Huber 法稳健估计并没有达到稳健反演结果的效果。

当采用一次范数法与遗传算法相结合时,无论异值出现在下沉盆地的什么位置,反演出来的结果离理论值都较近,比采用遗传算法和Huber 法求参效果都好,说明一次范数法一定程度上达到了稳健化的效果,降低了异值点或粗差的干扰。因此,表明一次范数法结合遗传算法适用于概率积分法参数反演,可以在一定程度上降低实测数据粗差对反演结果的影响,达到提高求解精度的效果。

3 结 论

针对利用遗传算法反演概率积分法参数时抗粗差能力差的问题,提出将稳健估计与遗传算法相结合的方法,提高反演参数过程中对异值点的抗干扰能力,保证参数求取结果稳健、精确。同时,设计实验对该方法的抗差能力进行检验,得出如下的结论。

(1)单纯采用遗传算法进行参数的反演时,当粗差出现在下沉盆地边缘时,其抗差能力较强;而当粗差出现在拐点处和下沉盆地中心时,其抗差能力较差,反演结果严重偏离理论值。说明当异值点在某种特定的位置时,遗传算法具有一定的抗差能力,但整体的抗差能力较差。以及说明拐点和最大下沉点的精度对于最终参数求取的结果影响较大,在实际野外测量时应提高这些特殊位置的测量精度,以保证求取的概率积分法参数的精度。

(2)采用一次范数法与遗传算法相结合时,其求取的结果比采用Huber 法与遗传算法进行结合时和单纯采用遗传算法的结果好,说明在一定程度上降低实测数据粗差对反演结果的影响,达到提高求解精度的效果。