Optimal selection of time functions for describing coal mining-induced dynamic subsidence at single surface point using AIC criterion
-
摘要:
时间函数法是最常用的煤矿区地表动态形变预测方法之一。其中，描述地表单点“S”型沉降过程的数学模型选取直接影响时间函数法的形变预测精度和可靠性。现有研究大都以拟合残差最小化为目标，通过修正或引入“S”型增长数学模型提高时间函数法的预测精度。然而，该方式容易导致“过拟合”现象，增加模型复杂度和参数反演难度。为克服该问题，引入了拟合残差和模型复杂度2个关键的优化选取指标，以来自7个不同地质采矿条件矿区的103个观测点时序沉降值为样本，采用理论分析与赤池信息量准则探讨了12种常见“S”型增长模型在煤矿区地表单点沉降过程描述中的优化选取。结果表明：①12种模型中，5个四参数模型的平均拟合残差为3.51 cm，明显优于二参数Knothe模型(14.10 cm)，但仅略优于6个三参数模型(4.78 cm)；②在兼顾拟合残差和模型复杂度的准则下，三参数模型普遍优于四参数和二参数，说明三参数模型比较适合描述矿区地表单点动态沉降过程，而四参数和二参数模型则分别存在“过拟合”(过度参数化)和“欠拟合”现象；③在6个三参数模型中，模型的优化选取与覆岩岩性有关。其中，软和中硬覆岩条件下，目前尚未被引入时间函数法的三参数Hossfeld模型能够较好地兼顾拟合残差小和模型复杂度低2个关键指标，但在坚硬覆岩条件下，Weibull模型表现则优于Hossfeld模型。
Abstract:The time function method is one of the most commonly used methods for predicting surface dynamic displacements in coal mine areas. In which, the accuracy and reliability of the predicted displacements, to a large extent, depends on the selected mathematical functions for describing the “S”-typed dynamic subsidence at a single surface point (referred to as time functions). Nearly all of the existing studies primarily improve or introduce “S”-shaped growth functions with a single object to minimizing the fitting residuals between thein-situmonitored and the model-fitted subsidence. Such a strategy, however, would result in “overfitting” (or over-parameterization), thereby increasing the complexity of the constructed time function model and the difficulty of model parameter inversion. To this end, the optimal selection of time functions was analyzed in this paper using two indicators of fitting residual and model complexity, rather than the former one in existing studies. More specifically, time-series subsidence observations at 103 field points in seven coal mining areas with different geological mining conditions were selected to be observation samples for ensuring the applicability of the optimal time function. Then, 12 common “S-shaped” growth models were chosen to candidates, and the theoretical analysis and Akaike information criterion (AIC) were further used to analyze the optimal selection of time function from the chosen 12 “S”-shaped models. The results show that: ① Among the 12 selected models, the mean mis-fitting error of the five four-parameter models is about 3.51 cm, which is obviously smaller than that of the two-parameter Knothe model (14.10 cm), but just slightly smaller than the six three-parameter models (4.78 cm); ② In the view of making a trade-off between fitting residuals and model complexity (assessing by the AIC), the AICs of the six three-parameter models are smaller than those of the four-parameter and two-parameter models.This indicates that the three-parameters models are preferrable to describe the temporal evolution of subsidence at a single point, and the four-parameter and two-parameters models may be over-fitted and under-fitted, respectively; ③ Among the six selected models, the optimal selection of time function is related to the lithology of the overburden rock strata; that is, Hossfeld model, which has not been introduced into the time function method, is preferrable under soft and medium-hard overburden strata, whereas Weibull model is preferrable under hard overburden strata.
-
0. 引 言
地下煤矿开采极易导致岩体移动变形，诱发地面沉降、山体滑坡、基础设施破坏、生态恶化等系列地质环境问题[1-2]。因此，准确认知上覆岩层及地表复杂的移动变形过程，并精准预测其演化趋势是矿山地质灾害精准防控的关键[3-4]。目前，煤矿区地表动态移动变形预测方法主要有物理力学法和时间函数法[2-3]。在两类方法中，物理力学法能够顾及矿山岩体物理力学变化过程，物理意义明确，但其高度依赖矿区地质条件刻画精细度和上覆岩体力学参数的可靠性，很大程度限制了该方法的实际应用。时间函数法是一种将描述地表单点动态沉降过程的时间函数（单点“动态”）与开采导致的最终沉降盆地（面域“静态”）融合形成的数学方法。该方法可通过历史形变观测值估计数学模型参数，并据此预测地表潜在移动变形趋势。由于时间函数法无需已知上覆岩层力学参数，因此被广泛应用于矿区地表动态移动变形预测[5-9]。
在时间函数法中，用于描述地表单点动态沉降过程的时间函数选取对于预测结果精度和适用范围起至关重要的作用。已有研究表明：长壁采煤（目前应用较多的地下采煤方法）引起的地表单点动态沉降基本服从“S”型增长。因此，Knothe函数[10]、Weibull函数[11]、Logistic函数[12]、Richards函数[13]、Boltzmann函数[14]和Bertalanffy函数[15]等“S”型增长数学模型被相继用于描述长壁开采导致的地表单点动态沉降过程，并据此修正时间函数法，取得了一定的成效。
然而，“S”型增长数学模型众多，且形式各异[16]。现有研究大都以拟合残差最小化（即强拟合性）为指标，以单个矿区少量沉降观测为样本，修正或引入“S”型数学模型以期提高时间函数法的预测精度[15-17]。该方式因过度追求拟合残差最小而易出现“过拟合”或“过度参数化”现象，进而增加构建模型的复杂度，极大提高参数反演难度。理论上，最优化时间函数模型应同时具备以下特征：①较好地描述地表单点沉降过程（强拟合性）；②适用于不同的地质采矿条件（高适应性）；③模型参数尽可能少（低复杂性）。因此，从时间函数优化选取角度来讲，现有研究不仅容易出现“过拟合”，而且以单个矿区少量观测样本为基础验证模型可靠性也容易出现所选模型在特定矿区应用效果较好，而适应性不高（即未考虑“高适应性”准则）。
鉴于此，引入了拟合残差和模型复杂度2个关键指标，通过选取不同地质采矿条件的大样本分析时间函数的优化选取。具体地，首先选取12种常见的“S”型函数模型为研究对象，推导和分析各模型的特征值与煤矿区地表单点动态沉降理论特征的差异。12种模型包含了已被引入时间函数法的7种模型以及目前尚未被引入但在生物种群增长等领域广泛使用的5种模型。在此基础上，以来自7个不同地质采矿条件煤矿区的103个地表时序沉降观测值为样本，利用赤池信息量准则（Akaike Information Criterion, AIC）系统地对比分析各模型的适应性、拟合优度和复杂度差异，并据此给出煤矿区地表单点动态沉降时间函数优化选取建议。研究有助于拓展时间函数法的应用范围，提高动态沉降预测精度，推动矿山动态地质灾害精准防控，丰富矿区开采沉陷理论。
1. 长壁采煤区地表单点沉降演化过程
地下长壁采煤导致的地表单点动态沉降服从“S”型函数增长，即该点大致经历初始沉降、加速沉降与残余沉降3个阶段[18]。如图1所示，在初始阶段，当地下开采即将影响到指定地表点时（
$t = 0$ ），该点的沉降量$W(t)$ 、沉降速度$v(t)$ 和加速度$a(t)$ 均为0；当地下开采开始影响该地表点时（$t > 0$ ），沉降逐渐显现，沉降速度逐渐增大，且沉降加速度逐渐增加到最大值；随着地下工作面继续推进，该点进入加速阶段，$W(t)$ 继续增大，$v(t)$ 逐渐加快至峰值并逐渐减缓，且沉降曲线拐点处下沉速度达到最快，加速度最小（此时$v(t) = {v_{\max }}$ 、$a(t) = 0$ ）；之后，随着时间的增长（$t \to \infty $ ），该点逐渐进入衰减阶段，沉降量$W(t)$ 继续增大且趋于常数，$v(t)$ 和$a(t)$ 渐趋于0，直至该点最终稳定。![]()
图 1 长壁采煤导致的地表单点沉降量、沉降速度和加速度演化过程示意Figure 1. Schematic of temporal evolution of subsidence, velocity and acceleration at single surface point caused by longwall coal mining2. 时间函数优化选取方法
为分析各模型的可靠性和复杂性，首先选取12种模型作为研究对象，其中包括目前已被引入时间函数法的7个模型（即Knothe、Weibull、Logistic、Bertalanffy、Richards、Gompertz和Boltzmann）以及目前尚未被引入时间函数法但在生物增长等领域常用的5个模型（即Hossfeld、Leavokic-I、Leavokic-III、Yoshida和Sloboda）。在12种模型中按参数个数分为二参数、三参数和四参数模型，其中，Knothe为二参数模型，Weibull、Logistic、Bertalanffy、Gompertz、Hossfeld和Leavokic-III（简称“Lea-III”）为三参数模型，而Leavokic-I（简称“Lea-I”）、Sloboda、Yoshida、Richards和Boltzmann为四参数模型。之后，将从理论分析和AIC准则辅助指导2个角度分析时间函数的优化选取。
2.1 理论特征分析
长壁采煤导致的地表不同点位动态沉降虽然在量级和过程等方面存在差异（与地质采矿条件相关），但正如第一节所述，各点在初始时刻（
$t = 0$ ）和稳定后（$t = + \infty $ ）的关键沉降量、沉降速度和加速度特性是一致的，即$W(t = 0) \equiv 0$ 、$v(t = 0) \equiv 0$ 、$a(t = 0) \equiv 0$ 、$W(t = + \infty ) \equiv $ C（C为常数）、$v(t = + \infty ) \equiv 0$ 和$a(t = + \infty ) \equiv 0$ ，本节将据此分析所选12种模型在$t = 0$ 和$t = + \infty $ 时的主要特征。为此首先根据12种所选模型的动态沉降表达式推导其对应的速度和加速度的表达式，并计算了其在$t = 0$ 和$t = + \infty $ 时的特征值，结果见表1。表 1 12种“S”型增长模型在$t = 0$ 和$t \to \infty $ 时的沉降量、速度和加速度特征Table 1. Comparison of subsidence, velocity, and acceleration of the selected 12 “S-shaped” time functions at the dates$t = 0$ and ,$t \to \infty $ respectively时间函数 下沉量$ y $ 下沉速度$ y' $ 下沉加速度$y''$ $ t=0 $ $ t \rightarrow \infty $ y y' y'' y y' y'' Knothe $y=a\left(1-{\rm{e}}^{-c t}\right)$ $y'=a c {\rm{e}}^{-c t}$ $y''=-a c^{2} {\rm{e}}^{-c t}$ 0 ac $ -a c^{2} $ a 0 0 Weibull $y=a\left(1-{\rm{e}}^{-b t^c}\right)$ $y'=a b c t^{c-1} {\rm{e}}^{-b t^{c} }$ $y''=a b c {\rm{e}}^{-b t^c} t^{c-2}\left[(c-1)-b c t^c\right]$ 0 0 0 a 0 0 Logistic $y=\dfrac{a}{1+b {\rm{e}}^{-c t} }$ $y'=\dfrac{a b c}{{\rm{e}}^{c t}\left(b {\rm{e} }^{-c t}+1\right)^2}$ $y''=\dfrac{a b c^2}{{\rm{e}}^{c t}\left(b {\rm{e}}^{-c t}-1\right)^2}\left[\dfrac{2 b}{{\rm{e}}^{c t}\left(b {\rm{e}}^{-c t}+1\right)}-1\right]$ $ \dfrac{a}{1+b} $ $ \dfrac{a b c}{(1+b)^{2}} $ $ \dfrac{a b c^{2}(c-1)}{(1+b)^{3}} $ a 0 0 Gompertz $y=a {\rm{e}}^{1-b {\rm{e}}^{-c t} }$ $y'=a b c {\rm{e} }^{-c t} {\rm{e}}^{\left(1-b {\rm{e} }^{-c t}\right)}$ $ y''=a b c^2 e^{-c t} e^{\left(1-b e^{-c t}\right)}\left(b e^{-c t}-1\right)$ $a {\rm{e}}^{1-b}$ $a b c {\rm{e}}^{1-b}$ $a b c^2(b-1)^* {\rm{e}}^{1-b}$ ae 0 0 Bertalanffy $y=a\left(1-b {\rm{e}}^{-c t}\right)^{-3}$ $y'=-3 a b c {\rm{e}}^{-c t}\left(1-b {\rm{e}}^{-c t}\right)^{-4}$ $\begin{aligned}y''= & 3 a b c^2 {\rm{e} }^{-c t}\left(1-b {\rm{e} }^{-c t}\right)^{-4}+ \\& 12 a b^2 c^2 {\rm{e} }^{-2 c t}\left(1-{\rm{e} }^{-c t}\right)^{-5}\end{aligned}$ $ \dfrac{a}{(1-b)^{3}} $ $ \dfrac{-3 a b c}{(1-b)^{4}} $ $ \dfrac{3 a b c^{2}}{(1-b)^{4}} $ a 0 0 Hossfeld $ y=\dfrac{t^c}{b+t^c / a}$ $y'=\dfrac{b c t^{c-1}}{\left(b+t^c / a\right)^2}$ $y''=y'\left[\dfrac{c-1}{t}-\dfrac{2 c t^{c-1} }{a\left(b+t^c / a\right)}\right]$ 0 0 0 a 0 0 Leavokic-Ⅲ $y=a\left(\dfrac{t^2}{b+t^2}\right)^c $ $ y'=\dfrac{2 b c y}{t\left(b+t^2\right)} $ $y''=2 b c \dfrac{t\left(b+t^2\right) y'-\left(3 t^2+b\right) y}{t^2\left(b+t^2\right)^2}$ 0 0 0 a 0 0 Leavokic-Ⅰ $ y=a\left(\dfrac{t^d}{b+t^d}\right)^c $ $y'=b c d \dfrac{y}{t\left(b+t^d\right)}$ $y''=\dfrac{t\left(b+t^d\right) y'-\left[(d+1) t^d+b\right] y}{t^2\left(b+t^d\right)^2}$ 0 0 0 a 0 0 Yoshida $y=\dfrac{a t^d}{b+t^d}+c $ $ y'=\dfrac{t^{d-1}}{\left(b+t^d\right)^2} $ $y''=\dfrac{t^{d-2} }{\left(b+t^d\right)}\left[(d-1)-\dfrac{2 d t^d}{b+t^d}\right]$ c 0 0 a+c 0 0 Sloboda $y=a {\rm{e}}^{-b {\rm{e}}^{-c t^d} }$ $y'=b c d t^{d-1} {\rm{e}}^{-c t^d} y$ $\begin{array}{l}y'' = {t^{d - 1} }{{\rm{e}}^{ - c{t^d} } }y' + \\\left[ {(d - 1){t^{d - 2} } - dc{t^{2d - 2} } } \right]{{\rm{e}}^{ - c{t^d} } }y\end{array}$ $a {\rm{e}}^{-b}$ 0 0 a 0 0 Richards $y=a\left(1-b {\rm{e}}^{-c t}\right)^{\dfrac{1}{1-d} }$ $\begin{array}{l}y' = {\left( {1 - b{ {\rm{e} }^{ - ct} } } \right)^{\frac{d}{ {1 - d} } } }\\{ {\rm{e} }^{ - ct} } y/(1 - d)\end{array}$ $\begin{array}{l}y'' = abc{\left[ {b{ { {{\rm{e}}} }^{ - ct} } - (1 - d)} \right]\times}\\{ {\rm{e} }^{ - ct} }{\left( {1 - b{e^{ - ct} } } \right)^{\left( {\frac{1}{ {1 - d} } - 2} \right)} }/{(1 - d)^2}\end{array}$ $a(1-b)^{\dfrac{1}{1-d}}$ $\dfrac{a b c(1-b)^{\dfrac{d}{1-d}}}{(1-d)}$ $\dfrac{a b c(b-1+d)(1-b)^{\left(\dfrac{1}{1-d}-2\right)} }{(1-d)^2}$ a 0 0 Boltzmann $y=b+\dfrac{a-b}{1+{\rm{e} }^{(t-c) / d} }$ $y'=\dfrac{(a-b) {\rm{e}}^{[(t-c) / d)]} }{d\left[1+{\rm{e}}^{(t-c) / d}\right]^2}$ $y''=y' \dfrac{{\rm{e}}^{((c-t) / d)}-1}{d {\rm{e}}^{[(c-t) / d+1]} }$ $\dfrac{a+b {\rm{e}}^{-\dfrac{c}{d} } }{{\rm{e}}^{-\dfrac{c}{d} }+1}$ $\dfrac{(a-b) {\rm{e}}^{-\dfrac{c}{d} } }{d \left({\rm{e}}^{-\dfrac{c}{d} }+1\right)^2}$ $\dfrac{{\rm{e}}^{(c / d-1)} }{{\rm{e}}^{(c / d+1)^3} } \dfrac{(a-b) {\rm{e}}^{(c / d} )}{d^2}$ a 0 0 注：a、b、c、d为时间函数模型参数。 从表1中可以看出，当
$t = + \infty $ 时，所选的12种模型均能满足沉降量、沉降速度和加速度理论特征（即$W(t = + \infty ) \equiv $ C（C为常数），$v(t = + \infty ) \equiv 0$ 和$a(t = + \infty ) \equiv 0$ ）。然而，当$t = 0$ 时，Weibull、Hossfeld、Lea-I和Lea-III模型的沉降值、沉降速度和加速度同时为零，Yoshida和Sloboda模型的初始沉降以及Knothe模型的初始速度和加速度为0，而Logistic、Gompertz、Bertalanffy、Richards和Boltzmann的初始沉降量、沉降速度和加速度均不为0。该结果表明，仅Weibull、Hossfeld、Leavokic-I和Leavokic-III能严格满足煤矿区地表单点在$t = 0$ 和$t = + \infty $ 时的动态沉降、沉降速度和加速度特征。2.2 基于AIC准则的时间函数优化选取
如前所述，最优化时间函数应同时兼顾强拟合性、高适应性和低复杂性3个关键指标。然而，现有研究大都仅通过拟合误差（Fitting Error）和
${R^2}$ 等拟合指标指导时间函数优化选取（即仅追求强拟合性），并未同时考虑强拟合性、高适应性和低复杂性之间的平衡。事实上，非线性模型强拟合性与低复杂性（即模型参数个数）之间通常是矛盾的，即追求强拟合性可能以增加模型参数个数（即复杂性）为代价，导致过拟合或“过度参数化”现象，进而产生模型参数难以准确反演或参数反演不稳定等关键问题。因此，如何能够平衡强拟合性和低复杂性2个重要指标是模型优化选取的关键。AIC准则是一种权衡所估计模型的复杂度和模型拟合优度有效方法[19-21]。该方法首先假设模型拟合误差服从正态分布，并在“熵”概念的基础上建立了AIC指标估计表达式[22]：
$$ A_{\rm{IC}} = {\rm{ln}}\frac{{R_{\rm{SS}}}}{n} + 2k $$ (1) 式中，
$R_{\rm{SS}}$ （Residual Sum of Squares）为残差平方和；$n$ 为样本数量；$k$ 为模型参数个数。从式（1）可以看出，AIC指标在考虑模型拟合可靠性（用$R_{\rm{SS}}$ 表示）的同时附加参数惩罚项，以此避免对数据的过拟合（Over-fitting）。在AIC准则中，最小AIC值对应的模型即被认为是当前观测样本下拟合可靠性和复杂性的最佳平衡模型。然而，受观测样本数量、观测误差和拟合误差等因素影响，AIC估计值通常含有误差。为减小误差对模型优化选取的影响，依据统计分级思想，通过以不同模型AIC值与最小值之差
$\Delta $ AIC（即$\Delta $ AIC(i)=AIC(i)−min(AIC), i=1, 2, 3, ···, 12）为分级标准[23]（具体见表2），依次将各个模型划分至3个不同的推荐使用等级（Level-I、Level-II和Level-III，其中Level-I为最优化模型），并据此分析12种模型在煤矿区地表单点动态沉降过程描述中的优化选取。表 2 基于$\Delta $ AIC分类参照标准Table 2. Reference based on$\Delta $ AIC classification水平 $\Delta $AIC 差异强弱 推荐使用等级 Level-I 0～2 弱差异 最高 Level-II 2～6 中等差异 适中 Level-III >6 强差异 较低 3. 真实数据分析与结果
3.1 试验数据概况
为顾及时间函数优化选取中的高适应性指标，在7个不同地质采矿条件长壁采煤矿区（详细信息见表3）收集了103个地表动态沉降观测值作为“S”型时间函数优化选取的样本。其中，7个矿区涵盖了我国比较常见的软弱、中硬和坚硬覆岩力学特征，且松散层厚度也涵盖了薄松散层到巨厚松散层的特征变化。来自不同地质采矿条件的大观测样本有助于验证时间函数优化选取的强适应性。同时，AIC准则的应用有助于选取12种模型中具有强拟合性和低复杂性的最优平衡模型，从而实现煤矿区单点动态沉降时间函数的优化选取（即同时兼顾高可靠性、强适应性和低复杂性）。
表 3 103个沉降观测值样本所处矿区的松散层和覆岩力学特征对比Table 3. Lithology of the seven interesting mining areas where subsidence observations at 103 surface points were collected矿区信息 松散层厚度 覆岩岩性 观测点数量 东营市某矿 巨厚松散层 中硬 26 宿州市某矿 厚松散层 中硬 15 榆林市某矿 厚黄土层 中硬 7 淮北市某矿 厚松散层 软弱 7 唐山市某矿 厚松散层 中硬 16 阳泉市某矿 薄松散层 坚硬 17 丰城市某矿 薄松散层 中硬 15 3.2 试验结果
在利用AIC准则优化选取时间函数之前，准确估计12种时间函数模型的参数至关重要。考虑到12种模型均具有高度非线性特征，因此使用遗传算法[24]（Genetic Algorithm, GA）和Levenberg-Marqurdt算法（LM）串行的非线性全局优化算法估计模型参数。具体地，为增大全局选优能力，首先设置遗传算法种群数量为20 000，最大迭代次数为5 000，交叉概率为0.7。根据该参数设定并利用遗传算法分别反演103个样本点处12个“S”型增长模型参数。之后，以遗传算法估计值作为初始参数，利用LM算法精化参数，提升全局反演能力。
图2展示了7个矿区随机选取的一个观测样本点处12个模型的拟合对比。从图2可以看出，除二参数Knothe模型外，其余11个模型基本能够描述随机选取的7个样本点的动态沉降过程。然而，从图2还可看出，除Knothe外的11个模型对于所展示的观测样本拟合优度仍存在差异，因此需着重研究描述煤矿区地表单点动态沉降过程的时间函数优化选取。
![]()
图 2 12种模型在7个矿区地表单个观测样本点动态沉降的拟合效果对比Figure 2. Comparison of 12 models for fitting time-series subsidence observations at single point in seven different mining areas为实现时间函数优化选取，利用式（1）计算了12种模型在各观测样本点的AIC值以及模型间的
$ \Delta $ AIC值。之后，参照表2所列$\Delta $ AIC分级规则将同一观测样本点的12种模型分级（即Level-I、II和III）。图3展示了12种“S”型模型在103个观测样本点所对应的不同分级所占比例。由图3可知，12种模型的$\Delta $ AIC均存在Level-I、II和III的分级。换句话说，面对不同地质采矿条件，12种模型均难以完美描述所有样本点的动态沉降过程（即Level-I占比为100%）。![]()
图 3 12种模型的$\Delta $ AIC分级占比统计Figure 3. Percent of selected 12 models at Level-I, II, and III based on$\Delta $ AIC estimates然而，从图3中可以看出，Hossfeld模型被划入Level-I等级的比例最高，达到了57.3%，且Level-I和Level-II高达87.3%。该结果表明，在12种时间函数模型中，Hossfeld模型能够较好地兼顾高适应性、强拟合性和低复杂性3个重要的优化选取指标。此外，从2.1节的理论分析可知，Hossfeld模型在
$t = 0$ 和$t = + \infty $ 时的曲线特征与长壁采煤导致的地表沉降过程理论特征值完全一致。因此，从兼顾强适应性、高可靠性和低复杂性的角度而言，Hossfeld模型（目前尚未被引入）可作为时间函数法中描述地表单点动态沉降过程的优化模型。4. 讨 论
4.1 覆岩结构与岩性对时间函数选取的影响
由3.2节的分析可知，在
$\Delta $ AIC分级规则下，Hossfeld模型在大部分（占比87.3%）观测样本中被划分至Level-I或II等级，是所有模型中占比最高的。然而，需要注意的是，在103个样本点中，Hossfeld模型在13个观测样本点（约占12.6%）被划定Level-III等级（即其他模型的$\Delta $ AIC指标比Hossfeld模型更优）。为进一步分析Hossfeld模型被划分为Level-III的原因，按照上覆岩层的岩性（即软、中硬和坚硬）和松散层厚度（薄、厚和巨厚）将13个Level-III等级的观测样本点分类。分析发现在13个Level-III等级的观测样本中，绝大部分点（8个观测样本点）位于坚硬覆岩和薄松散层覆盖的阳泉某矿区。图4显示了阳泉某矿区两个典型动态沉降观测序列以及Hossfeld（处于Level-III水平）和Weibull（处于Level-I水平）2个模型的拟合结果，图中动态沉降为阳泉某矿区随机抽取2个样本点的动态沉降观测值。由图4可知，Hossfeld模型对于样本点的沉降初始阶段（50~100 d）和残余沉降阶段描述效果均不及Weibull模型。根据式（1）可知，在相同的模型参数个数（Hossfeld和Weibull模型均为三参数）和时序沉降观测值的前提下，Weibull模型的拟合残差越小，其对应的AIC值也越小。因此，在
$\Delta $ AIC分级规则下即被划入了更优等级。![]()
图 4 Level-III等级下Hossfeld和Weibull模型的拟合效果对比Figure 4. Comparison diagram of fitting effects of Hossfeld and Weibull models at Level-III Level造成Hossfled模型被划入Level-III的主要原因如下。相较于较软或中硬上覆岩层，坚硬覆岩和薄松散层的地质条件对于上覆岩层的支撑力更强，因此其地表点位达到最大下沉速度和最大沉降量的时间较为滞后。从表1中可以发现，Hossfeld模型动态沉降“拐点”（即最大沉降速度发生时刻）处的沉降值约为最大沉降的37%，而Weibull模型则可达50%。这就意味着，Weibull模型对于拐点出现较为滞后（比如坚硬覆岩和/或薄松散层）的描述效果比Hossfeld模型更好。因此，坚硬覆岩和薄松散层地质条件下可优选Weibull模型作为优化时间函数。
4.2 仅考虑拟合残差的时间函数选取
如前所述，AIC准则实际上是通过附加参数惩罚项以避免“过拟合”问题，基于AIC准则选取的时间函数能够在保证拟合精度的同时最简化函数模型，提高计算效率和参数反演稳健性。为对比不同准则对于时间函数模型选取的影响，本节仅考虑强拟合性（即绝大部分现有研究选取的模型优选准则）分析12种“S”型增长模型的拟合效果和模型选取。具体地，利用3.2节中所述的算法反演103个样本点处12种“S”型增长模型参数，计算观测值与拟合值之间的均方根误差（Root Mean Square Error，RMSE），该值越小，拟合优度越好，反之亦然。图5显示了12种函数模型在各样本点的拟合均方根误差统计直方图。
![]()
图 5 12种选取模型在选取的样本点的拟合均方根误差统计直方图Figure 5. Frequency of root mean square errors of the selected 12 “S-shaped” models for fitting the dynamic subsidence at the selected 103 observation points由图5可知，除Knothe函数外，其余函数的拟合均方根误差绝大部分集中在0~0.05 m，其中四参数模型（Lea-I、Yoshida、Sloboda、Richards、Boltzmann）拟合均方根误差在此区间的平均样本点数为85个，占比82.5%；三参数模型（Gompertz、Bertalanffy、Hossfeld、Lea-III、Weibull和Logisitic）拟合均方根误差在0~0.05 m区间平均样本点数约77个，占比74.8%，略低于四参数模型。为进一步分析12种模型的拟合优度差异，采用“分位数断点法”将每个点的拟合优度差异分为3个等级。具体地，对于任意样本点，分位数断点法首先将所有模型在该点的拟合均方根误差按照从小到大排序。然后以33.3%、66.6%和100%三个分位数值作为分段“断点”，并据此将各模型划定为Level-I（[0~33.3%)）、Level-II（[33.3%~66.6%)）和Level-III（[66.6%~100%]）三个等级。其中，Level-I表示该模型在12种“S”型增长模型中拟合优度最好，Level-II次之，Level-III较差。
图6展示了12种模型在103个观测样本点被划定为Level-I到III的相对比例。由图6可知，12种“S”型模型中四参数模型（即Lea-I、Yoshida、Slobada、Richards和Boltzmann模型）的总体拟合优度（Level-I的平均比例为42.7%）要高于三参数（Level-I平均比例为27.7%）和二参数模型（Level-I平均比例为5.08%）。事实上，该结果也是预期的，因为相较于二参数和三参数模型，四参数“S”型函数增加了未知数个数，并且可通过调整参数值减少样本的拟合残差。
![]()
图 6 12种模型拟合均方根误差划定的Level-I、Level-II和Level-III的占比对比Figure 6. Comparison of the RMSE levels (I, II, and III) of the selected 12 models此外，从图6中还可以看出，相较于其它四参数模型，Slobada整体拟合精度较高（Level-I和Level-II占比达到了87.5%）；相较于其他三参数模型，Hossfeld的整体拟合精度最高（Level-I和Level-II占比达到了87.2%）。然而，通过对比Slobada和Hossfeld模型拟合优度发现，2个模型被划定为Level-I和Level-II的比例相当（仅相差0.3%）。此外，从图5中可以看出，Slobada模型平均拟合均方根误差为0.0315 m，仅比Hossfeld模型（0.0368 m）低5.3 mm。换句话说，四参数Slobada模型以增加模型复杂性（新增一个未知参数）为代价，但其平均拟合精度仅提高了5.3 mm。考虑到采煤沉陷区地表最大沉降量级通常可达几十甚至数百厘米，因此，5.3 mm的精度改善对于拟合结果影响甚微。
5. 结 论
1）Weibull、Hossfeld、Lea-I和Lea-III模型能完全满足煤矿区地表单点在
$t = 0$ 和$t = + \infty $ 时的动态沉降、沉降速度和加速度理论极值特征。2）在AIC准则下，12种模型中，三参数的Hossfeld函数能够较好地兼顾强拟合性、高适应性和低复杂性3个优化选取关键指标。
3）当地质采矿未知时，可优选Hossfeld模型作为描述煤矿区地表单点动态沉降的时间函数；当地质采矿条件已知且上覆岩层坚硬时，优选Weibull模型，其他条件下可优选Hossfeld模型。
4）相较于三参数Hossfeld模型，四参数Slobada模型（四参数表现最优）在新增一个未知数的代价下拟合精度仅提高了5.3 mm，存在“过度参数化”现象。
-
![]()
图 1 长壁采煤导致的地表单点沉降量、沉降速度和加速度演化过程示意
Figure 1. Schematic of temporal evolution of subsidence, velocity and acceleration at single surface point caused by longwall coal mining
![]()
图 2 12种模型在7个矿区地表单个观测样本点动态沉降的拟合效果对比
Figure 2. Comparison of 12 models for fitting time-series subsidence observations at single point in seven different mining areas
![]()
图 3 12种模型的
$\Delta $ AIC分级占比统计Figure 3. Percent of selected 12 models at Level-I, II, and III based on
$\Delta $ AIC estimates![]()
图 4 Level-III等级下Hossfeld和Weibull模型的拟合效果对比
Figure 4. Comparison diagram of fitting effects of Hossfeld and Weibull models at Level-III Level
![]()
图 5 12种选取模型在选取的样本点的拟合均方根误差统计直方图
Figure 5. Frequency of root mean square errors of the selected 12 “S-shaped” models for fitting the dynamic subsidence at the selected 103 observation points
![]()
图 6 12种模型拟合均方根误差划定的Level-I、Level-II和Level-III的占比对比
Figure 6. Comparison of the RMSE levels (I, II, and III) of the selected 12 models
表 1 12种“S”型增长模型在
$t = 0$ 和$t \to \infty $ 时的沉降量、速度和加速度特征Table 1 Comparison of subsidence, velocity, and acceleration of the selected 12 “S-shaped” time functions at the dates
$t = 0$ and ,$t \to \infty $ respectively时间函数 下沉量$ y $ 下沉速度$ y' $ 下沉加速度$y''$ $ t=0 $ $ t \rightarrow \infty $ y y' y'' y y' y'' Knothe $y=a\left(1-{\rm{e}}^{-c t}\right)$ $y'=a c {\rm{e}}^{-c t}$ $y''=-a c^{2} {\rm{e}}^{-c t}$ 0 ac $ -a c^{2} $ a 0 0 Weibull $y=a\left(1-{\rm{e}}^{-b t^c}\right)$ $y'=a b c t^{c-1} {\rm{e}}^{-b t^{c} }$ $y''=a b c {\rm{e}}^{-b t^c} t^{c-2}\left[(c-1)-b c t^c\right]$ 0 0 0 a 0 0 Logistic $y=\dfrac{a}{1+b {\rm{e}}^{-c t} }$ $y'=\dfrac{a b c}{{\rm{e}}^{c t}\left(b {\rm{e} }^{-c t}+1\right)^2}$ $y''=\dfrac{a b c^2}{{\rm{e}}^{c t}\left(b {\rm{e}}^{-c t}-1\right)^2}\left[\dfrac{2 b}{{\rm{e}}^{c t}\left(b {\rm{e}}^{-c t}+1\right)}-1\right]$ $ \dfrac{a}{1+b} $ $ \dfrac{a b c}{(1+b)^{2}} $ $ \dfrac{a b c^{2}(c-1)}{(1+b)^{3}} $ a 0 0 Gompertz $y=a {\rm{e}}^{1-b {\rm{e}}^{-c t} }$ $y'=a b c {\rm{e} }^{-c t} {\rm{e}}^{\left(1-b {\rm{e} }^{-c t}\right)}$ $ y''=a b c^2 e^{-c t} e^{\left(1-b e^{-c t}\right)}\left(b e^{-c t}-1\right)$ $a {\rm{e}}^{1-b}$ $a b c {\rm{e}}^{1-b}$ $a b c^2(b-1)^* {\rm{e}}^{1-b}$ ae 0 0 Bertalanffy $y=a\left(1-b {\rm{e}}^{-c t}\right)^{-3}$ $y'=-3 a b c {\rm{e}}^{-c t}\left(1-b {\rm{e}}^{-c t}\right)^{-4}$ $\begin{aligned}y''= & 3 a b c^2 {\rm{e} }^{-c t}\left(1-b {\rm{e} }^{-c t}\right)^{-4}+ \\& 12 a b^2 c^2 {\rm{e} }^{-2 c t}\left(1-{\rm{e} }^{-c t}\right)^{-5}\end{aligned}$ $ \dfrac{a}{(1-b)^{3}} $ $ \dfrac{-3 a b c}{(1-b)^{4}} $ $ \dfrac{3 a b c^{2}}{(1-b)^{4}} $ a 0 0 Hossfeld $ y=\dfrac{t^c}{b+t^c / a}$ $y'=\dfrac{b c t^{c-1}}{\left(b+t^c / a\right)^2}$ $y''=y'\left[\dfrac{c-1}{t}-\dfrac{2 c t^{c-1} }{a\left(b+t^c / a\right)}\right]$ 0 0 0 a 0 0 Leavokic-Ⅲ $y=a\left(\dfrac{t^2}{b+t^2}\right)^c $ $ y'=\dfrac{2 b c y}{t\left(b+t^2\right)} $ $y''=2 b c \dfrac{t\left(b+t^2\right) y'-\left(3 t^2+b\right) y}{t^2\left(b+t^2\right)^2}$ 0 0 0 a 0 0 Leavokic-Ⅰ $ y=a\left(\dfrac{t^d}{b+t^d}\right)^c $ $y'=b c d \dfrac{y}{t\left(b+t^d\right)}$ $y''=\dfrac{t\left(b+t^d\right) y'-\left[(d+1) t^d+b\right] y}{t^2\left(b+t^d\right)^2}$ 0 0 0 a 0 0 Yoshida $y=\dfrac{a t^d}{b+t^d}+c $ $ y'=\dfrac{t^{d-1}}{\left(b+t^d\right)^2} $ $y''=\dfrac{t^{d-2} }{\left(b+t^d\right)}\left[(d-1)-\dfrac{2 d t^d}{b+t^d}\right]$ c 0 0 a+c 0 0 Sloboda $y=a {\rm{e}}^{-b {\rm{e}}^{-c t^d} }$ $y'=b c d t^{d-1} {\rm{e}}^{-c t^d} y$ $\begin{array}{l}y'' = {t^{d - 1} }{{\rm{e}}^{ - c{t^d} } }y' + \\\left[ {(d - 1){t^{d - 2} } - dc{t^{2d - 2} } } \right]{{\rm{e}}^{ - c{t^d} } }y\end{array}$ $a {\rm{e}}^{-b}$ 0 0 a 0 0 Richards $y=a\left(1-b {\rm{e}}^{-c t}\right)^{\dfrac{1}{1-d} }$ $\begin{array}{l}y' = {\left( {1 - b{ {\rm{e} }^{ - ct} } } \right)^{\frac{d}{ {1 - d} } } }\\{ {\rm{e} }^{ - ct} } y/(1 - d)\end{array}$ $\begin{array}{l}y'' = abc{\left[ {b{ { {{\rm{e}}} }^{ - ct} } - (1 - d)} \right]\times}\\{ {\rm{e} }^{ - ct} }{\left( {1 - b{e^{ - ct} } } \right)^{\left( {\frac{1}{ {1 - d} } - 2} \right)} }/{(1 - d)^2}\end{array}$ $a(1-b)^{\dfrac{1}{1-d}}$ $\dfrac{a b c(1-b)^{\dfrac{d}{1-d}}}{(1-d)}$ $\dfrac{a b c(b-1+d)(1-b)^{\left(\dfrac{1}{1-d}-2\right)} }{(1-d)^2}$ a 0 0 Boltzmann $y=b+\dfrac{a-b}{1+{\rm{e} }^{(t-c) / d} }$ $y'=\dfrac{(a-b) {\rm{e}}^{[(t-c) / d)]} }{d\left[1+{\rm{e}}^{(t-c) / d}\right]^2}$ $y''=y' \dfrac{{\rm{e}}^{((c-t) / d)}-1}{d {\rm{e}}^{[(c-t) / d+1]} }$ $\dfrac{a+b {\rm{e}}^{-\dfrac{c}{d} } }{{\rm{e}}^{-\dfrac{c}{d} }+1}$ $\dfrac{(a-b) {\rm{e}}^{-\dfrac{c}{d} } }{d \left({\rm{e}}^{-\dfrac{c}{d} }+1\right)^2}$ $\dfrac{{\rm{e}}^{(c / d-1)} }{{\rm{e}}^{(c / d+1)^3} } \dfrac{(a-b) {\rm{e}}^{(c / d} )}{d^2}$ a 0 0 注：a、b、c、d为时间函数模型参数。 
下载: 导出CSV
表 2 基于
$\Delta $ AIC分类参照标准Table 2 Reference based on
$\Delta $ AIC classification水平 $\Delta $AIC 差异强弱 推荐使用等级 Level-I 0～2 弱差异 最高 Level-II 2～6 中等差异 适中 Level-III >6 强差异 较低
下载: 导出CSV
表 3 103个沉降观测值样本所处矿区的松散层和覆岩力学特征对比
Table 3 Lithology of the seven interesting mining areas where subsidence observations at 103 surface points were collected
矿区信息 松散层厚度 覆岩岩性 观测点数量 东营市某矿 巨厚松散层 中硬 26 宿州市某矿 厚松散层 中硬 15 榆林市某矿 厚黄土层 中硬 7 淮北市某矿 厚松散层 软弱 7 唐山市某矿 厚松散层 中硬 16 阳泉市某矿 薄松散层 坚硬 17 丰城市某矿 薄松散层 中硬 15
下载: 导出CSV
-
[1] 汪云甲. 矿区生态扰动监测研究进展与展望[J]. 测绘学报,2017,46(10):1705−1716. doi: 10.11947/j.AGCS.2017.20170358 WANG Yunjia. Research progress and prospect on ecological disturbance monitoring in mining area[J]. Acta Geodaetica et Cartographica Sinica,2017,46(10):1705−1716. doi: 10.11947/j.AGCS.2017.20170358
[2] SENGUPTA M. Environmental impacts of mining: Monitoring, Restoration and Control[M]. Baco Roton: CRC Press, 2021: 1-30.
[3] SHI Y,ZHAO M,HAO J. Study on numerical models in predicting surface deformation caused by underground coal mining[J]. Geotechnical and Geological Engineering,2021,39(6):4457−4473. doi: 10.1007/s10706-021-01775-2
[4] 朱建军,杨泽发,李志伟. InSAR矿区地表三维形变监测与预计研究进展[J]. 测绘学报,2019,48(2):135−144. doi: 10.11947/j.AGCS.2019.20180188 ZHU Jianjun,YANG Zefa,LI Zhiwei. Recent progress in retrieving and predicting mining-induced 3D displacements using InSAR[J]. Acta Geodaetica et Cartographica Sinica,2019,48(2):135−144. doi: 10.11947/j.AGCS.2019.20180188
[5] 刘玉成,曹树刚,刘延保. 可描述地表沉陷动态过程的时间函数模型探讨[J]. 岩土力学,2010,31(3):925−931. doi: 10.3969/j.issn.1000-7598.2010.03.044 LIU Yucheng,CAO Shugang,LIU Yanbao. Discussion on some time functions for describing dynamic course of surface subsidence due to mining[J]. Rock and Soil Mechanics,2010,31(3):925−931. doi: 10.3969/j.issn.1000-7598.2010.03.044
[6] 崔希民,缪协兴,赵英利,等. 论地表移动过程的时间函数[J]. 煤炭学报,1999,24(5):453−456. doi: 10.3321/j.issn:0253-9993.1999.05.002 CUI Ximin,MIAO Xiexing,ZHAO Yingli,et al. Discussion on the time function of time dependent surface movement[J]. Journal of China Coal Society,1999,24(5):453−456. doi: 10.3321/j.issn:0253-9993.1999.05.002
[7] YANG Z,LI Z,ZHU J,et al. An InSAR-Based temporal probability integral method and its application for predicting mining-induced dynamic deformations and assessing progressive damage to surface buildings[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2018,11(2):472−484. doi: 10.1109/JSTARS.2018.2789341
[8] JIANG C,WANG L,YU X,et al. Prediction of 3D deformation due to large gradient mining subsidence based on InSAR and constraints of IDPIM model[J]. International Journal of Remote Sensing,2021,42(1):208−239. doi: 10.1080/01431161.2020.1804088
[9] CUI X,WANG J,LIU Y. Prediction of progressive surface subsidence above longwall coal mining using a time function[J]. International Journal of Rock Mechanics and Mining Sciences,2001,38(7):1057−1063. doi: 10.1016/S1365-1609(01)00061-2
[10] 郭旭炜,杨晓琴,柴双武. 分段Knothe函数优化及其动态求参[J]. 岩土力学,2020,41(6):2091−2097. GUO Xuwei,YANG Xiaoqin,CHAI Shuangwu. Optimization of the segmented Knothe function and its dynamic parameter calculation[J]. Rock and Soil Mechanics,2020,41(6):2091−2097.
[11] 刘东海,邓念东,姚 婷,等. 基于开采沉陷实测数据的Weibull时间函数模型参数研究[J]. 煤炭科学技术,2021,49(9):152−158. LIU Donghai,DENG Niandong,YAO Ting,et al. Study on parameters of Weibull time function model based on sited measured mining subsidence data[J]. Coal Science and Technology,2021,49(9):152−158.
[12] 李春意,赵 亮,李 铭,等. 基于Logistic时间函数地表动态沉陷预测及优化求参研究[J]. 安全与环境学报,2020,20(6):2202−2210. doi: 10.13637/j.issn.1009-6094.2019.1105 LI Chunyi,ZHAO Liang,LI Ming,et al. Prediction of surface progressive subsidence and optimization of predicting model parameters based on the Logistic time function[J]. Journal of Safety and Environment,2020,20(6):2202−2210. doi: 10.13637/j.issn.1009-6094.2019.1105
[13] 卢克东,徐良骥,牛亚超. 基于GA-PSO融合算法的开采沉陷Richards预计模型参数优化[J]. 金属矿山,2021(2):155−160. LU Kedong,XU Liangji,NIU Yachao. Parameter optimization on Richards model of mining subsidence based on GA-PSO Hybrid algorithm[J]. Metal Mine,2021(2):155−160.
[14] 王 宁,吴 侃,刘 锦,等. 基于Boltzmann函数的开采沉陷预测模型[J]. 煤炭学报,2013,38(8):1352−1356. doi: 10.13225/j.cnki.jccs.2013.08.018 WANG Ning,WU Kan,LIU Jin,et al. Model for mining subsidence prediction based on Boltzmann function[J]. Journal of China Coal Society,2013,38(8):1352−1356. doi: 10.13225/j.cnki.jccs.2013.08.018
[15] 高 超,徐乃忠,孙万明,等. 基于Bertalanffy时间函数的地表动态沉陷预测模型[J]. 煤炭学报,2020,45(8):2740−2748. GAO Chao,XU Naizhong,SUN Wanming,et al. Dynamic surface subsidence prediction model based on Bertalanffy time function[J]. Journal of China Coal Society,2020,45(8):2740−2748.
[16] BORIS Z. Analysis of Growth Equations[J]. Forest Science,1993(03):594−616.
[17] 杨泽发,易辉伟,朱建军,等. 基于InSAR时序形变的矿区全盆地沉降时空演化规律分析[J]. 中国有色金属学报,2016,26(7):1515−1522. doi: 10.19476/j.ysxb.1004.0609.2016.07.020 YANG Zefa,YI Huiwei,ZHU Jianjun,et al. Spatio-temporal evolution law analysis of whole mining subsidence basin based on InSAR-derived time-series deformation[J]. The Chinese Journal of NonferrousMetals,2016,26(7):1515−1522. doi: 10.19476/j.ysxb.1004.0609.2016.07.020
[18] 黄乐亭. 地表动态沉陷变形的三个阶段与规律[J]. 矿山测量,2003(3):18−20. doi: 10.3969/j.issn.1001-358X.2003.03.006 HUANG Leting. Research on three stages and law of dynamic surface subsidence deformation[J]. Mine Surveying,2003(3):18−20. doi: 10.3969/j.issn.1001-358X.2003.03.006
[19] 刘璋温. 赤池信息量准则 AIC 及其意义[J]. 数学的实践与认识,1980(3):64−72. LIU Zhangwen. Akaike’s Information Criterion and its significance[J]. Mathematics in Practice and Theory,1980(3):64−72.
[20] TjØRVE E,TjØRVE K M C. A unified approach to the Richards-model family for use in growth analyses: Why we need only two model forms[J]. Journal of Theoretical Biology,2010,267(3):417−425. doi: 10.1016/j.jtbi.2010.09.008
[21] KARLSSON P S,BEHRENZ L,SHUKUR G. Performances of model selection criteria when variables are Ill conditioned[J]. Computational Economics,2019,54(1):77−98. doi: 10.1007/s10614-017-9682-8
[22] KAMO K-I,YOSHIMOTO A. Comparative analysis on selecting growth function based on three different information criteria for the purpose of carbon estimation[J]. Forest Science and Technology,2013,9(2):65−71. doi: 10.1080/21580103.2013.801165
[23] VRANA J,REMSE V,MATYSIOKOVA B,et al. Choosing the right sigmoid growth function using the unified-models approach[J]. International Journal of Avian Science,2019,161(1):13−26.
[24] KATOCH S,CHAUHAN S S,KUMAR V. A review on genetic algorithm: past, present, and future[J]. Multimedia Tools and Applications,2021,80(5):8091−8126. doi: 10.1007/s11042-020-10139-6
-
期刊类型引用(1)
1. 赵月，王志伟，张国建，王翔，丁文壮. Hossfeld模型在矿区地表动态沉降预测应用的可行性分析. 中国矿业. 2024(03): 124-132 . 
百度学术
其他类型引用(3)
计量
- 文章访问数: 82
- HTML全文浏览量: 14
- PDF下载量: 27
- 被引次数: 4
