杜子龙 仪梦婷 李志远

(1. 山东科技大学 测绘科学与工程学院, 山东 青岛 266590; 2. 中国海洋大学 文学与新闻传播学院, 山东 青岛 266100)

0 引言

矿区地下开采引起的地表沉陷是一个复杂的四维空间问题,沉陷预测为采空区上方的建筑物,铁路、水库等采后支护与治理提供了理论依据,很多研究中注重于开采沉陷的最终结果,即静态地表移动与变形预测。静态预测侧重于观测的整体性,对开采沉陷的动态过程研究尚不充分,这对地表建筑物的保护是不利的。实际上,由于矿山开采导致的地表移动与变形是一个动态的过程,变形与移动越剧烈,地表建筑物的损坏就越严重[1]。地表的下沉具有三个过程:缓慢发展期、加速发展期、发展衰减期[2]。对地面建筑物等的保护需要考虑变形随时间发展的全过程,目前常用动态预测地表下沉的时间函数有波兰学者Knothe[3]于1952年提出的Knothe时间函数;Sroka与 Schober于1982年提出的Sroka-Schober[4]时间函数,Kowalski于1999年提出的广义时间函数;Gonzalfz等在2007年研究使用的正态分布函数作为时间函数等;我国在采动过程中采用的预计函数多数为基于Knothe时间函数的改进预测模型,大量专家与学者针对Knothe时间函数无法反应下沉速度与加速度的缺陷提出了各种改进方法,常占强等提出的分段Knothe时间函数,提高了预测精度;张兵等改进了分段Knothe函数,使之居于更大的适用范围,并提高了精度[5-8]。但Knothe时间函数观测周期长,模型复杂,所需参数多。基于经验函数的动态预计方法公式计算复杂,且其精度受限于开采实际经验及岩移参数。基于Logistics模型的下沉预测时间函数能够利用较少的观测量提早预测开采沉陷的沉陷量,能够反映地表下沉的全过程,对于采动过程中引发的地表移动与变形的预测和控制的研究有着积极意义。 Logistic模型能够满足单个点动态下沉预测的要求,模型简单,参数少,应用方便。本文通过控制参与拟合的数据个数,一组使用72期的数据进行拟合,另一组使用81期的数据进行拟合,得到模型预测方程,与实测值的对比表明81期数据组的预计方程要优于72期数据组,通过预测方程对82期91期的沉降值进行预测,得到了较高的预测精度。

1 Logistics模型概述

在当前的社会学,生物学,计量经济学等方面,逻辑回归模型都有广泛的应用。Logistic模型在最初研究人口增长规律时由马尔萨斯提出, 之后荷兰数学家威赫尔斯特将其归纳为一般的数学公式,Logistic曲线在有限空间内的数值的增长会趋于一个稳定值。其发展具有三个过程:缓慢发展期、加速发展期、发展衰减期,它反映了事物的发生、发展、成熟并趋于稳定的过程[9]。

Logistic增长模型的微分方程为:

(1)

对该微分形式进行积分求解可得到该方程的解如下:

(2)

式(2)所表示的曲线即为Logistic模型,总体上呈现S型或反S型曲线(图1),其特点是单调递增(S型)或单调递减(反S型)。由于该曲线为S型,具有地表下沉的相关特征,曲线经历缓慢发展,加速发展,最终趋于稳定的三个阶段,可以将该函数看为Logistic地表下沉时间函数,其中y代表采动过程中该点任意时间的下沉值,参数K为该点最终沉降量W0,参数A跟b为下沉影响因子。其表达式如下:

(3)

图1 Logistic模型示意图

2 基于 Logistic模型的预计流程

Logistic模型曲线形态规律与采动引起的地表下沉规律相符合[10],运用Logistics模型可以预计地表下沉,通过预测的值和实测的值之间的对比来验证模型。Logistic模型预计流程如图2所示。

图2 Logistic模型预计流程

3 Logistic模型参数估计

Logistic曲线是以参数K为渐近线,连续的、单调递减(增)的曲线,曲线随时间T的变化速度由参数a和b确定。对于式(2)中的三个未知参数,本文将采用在倒数求和法的基础上,使用非线性回归的方法,估计模型参数。

(1)倒数求和法确定迭代初始值

为了提高精度,对W0,A,b初值的估计采用倒数求和法[11]。倒数求和法是将样本均分为三等分,要求时间间隔一样,设样本为(ti,Wi),i=1,2,…,n;令r=int(n/3),当样本数量N不能等分为三部分时,可以适当舍弃开始的1～2个数。

(4)

令D1=S1-S2,D2=S2-S3,则

(5)

(6)

(7)

(2)非线性回归

为了提高拟合精度,迭代的方法采用Levenberg-Marquardt算法(简称L-M算法)。L-M算法是一种以最小二乘原理对多个参数进行拟合的方法,它是介于牛顿法与梯度下降法之间的一种非线性优化方法。L-M算法结合了最速下降法与线性方法,能有效处理冗余参数问题,保证在值域范围内快速收敛[12]。用L-M算法求得各参数的最佳估计值可以作为Logistic模型参数的最小二乘无偏性估计值。

4 实例验证

对某矿区观测线上的某监测点经过91个周期的地表沉降监测,该监测点的累计沉降量达到38.7 cm。借助Logistic模型进行拟合,将A组使用72期的数据进行拟合并预测沉降,B组使用81期的数据进行模型验证,并分别预测82期至91期的沉降值,本文采用计算机模拟分析参数。分别得到A、B组的预测函数与实测值的拟合方程及成果图(图3)及根据两组拟合方程得到的82期至91期的沉降值预测曲线图(图4)。

图3 某个地表沉降监测点实测值与A组与B组Logistic拟合结果图

图4 某个地表沉降监测点实测值与Logistic模型预测值图

WA=-0.368 555/(1+14.790 1·
exp(-0.0787 091·t))

WB=-0.381 006/(1+13.702 2·
exp(-0.074 0703·t))

如图3所示,拟合曲线WA、WB与该点地表下沉实测值在三个阶段内的吻合度较高,曲线WB较WA拟合程度更高,原因是WB较WA选取了较多的实测数据,并且具有较多处于下沉第三个阶段的实测数据,从预测结果来看,如图4,WB也比WA更接近真实值,拟合程度高相应的预测精度也高,这说明在本例中,单点下沉符合Logistic下沉时间曲线,两者具有一致的规律性。

根据统计结果显示出B组采用的81期数据拟合方程要优于A组采用72期数据拟合方程,其误差平方和及误差均方都较A组更接近于0,表明其拟合程度要优于采用较少观测点的拟合方程。B组参数估计成果表如表1,结果显示运用L-M算法对三点法估计值进行修正,得到了较好的效果。B组拟合方程统计结果显示采用该拟合方法的误差平方和SSE为0.006 621 8,误差均方MSE为0.000 087 1。通过A、B组Logistic预测方程分别得出82期到91期沉降预测值,并与实测沉降值进行比较,统计误差,结果见表2。

表1 B组参数估计成果表

(1)表1、表2结果表明,采用81期数据的预测结果要优于采用72天数据的预测结果,两组的预测差在1 cm左右,结果表明采用该方法时,选用的数据越多,拟合的程度就越高,预测的结果也会更准确。

表2 实测值与预测值对比 单位：m

(2)在采用72期的实测数据进行拟合预测的结果表明,采用较少的数据进行预测时预测值与实测值的误差为2 cm左右,同样具有比较高的精度,表明在应用该模型进行预测时,可以采用较少的数据来进行地表下沉预测,选取的地表沉降值数据组将对预测结果产生直接影响。

(3)B组实验组采用该拟合方法的误差平方和SSE为0.006 621 8,误差均方MSE为0.000 087 1,接近于0,拟合精度较高,说明该方法对于动态预测具有一定的精度可靠性。

(4)由表2结果知B组预测函数与实测值的最大残差出现在第91期,为0.012 067 3 m,最小残差出现第85期,为0.007 882 1 m,平均误差仅为0.01 m,具有较高的预测精度。

5 结束语

(1)通过Logistics模型进行参数拟合,运用L-M算法对参数进行修正,统计结果显示效果较好,并得到了动态预测方程。

(2)从实测数据与拟合曲线图可以看出,单点下沉关系曲线为反S型。A、B组误差平方和SSE与均方MSE均接近于0,拟合精度较高。

(3)在大量的实测数据的基础上,Logistic模型对于单点的地表沉降具有较好的拟合结果,预测精度较高,B组模型预测差最大为0.012 067 3 m,平均误差为0.01 m。对矿区地表建筑物损坏的实时预防和维护具有实际的指导意义。

(4)对比信息表明采用该方法进行预测时,选用的数据越多,拟合程度越优,沉陷预测结果也更接近于真实值,选取不同的沉降值参与拟合会直接影响预测结果。

(5)Logistics模型简单,参数少,可以适当减少观测次数,对地表下沉的动态过程具有一定的使用价值。