煤炭资源开采引起的地表沉降在附加荷载的影响下会导致地表、建(构)筑物荷载、地基等发生不均匀沉降，为确保建筑物安全、反馈设计参数，需对其进行变形监测^[1]。依据监测数据，用模型表达变形时序与变化量之间的不确定性关系，对于施工质量监测和运行安全保障具有重要意义。邓聚龙^[2]提出的灰色系统理论被广泛应用于变形监测领域。

建筑物形变量一般会随时间发生非线性变化，将GM(1, 1)模型中的灰作用量替换为幂指数形式，能够在样本数量较少的情况下灵活反映数据序列之间的非线性关系^[3]。为提高幂模型的精度，相关学者从不同角度对其进行优化改进，常用的改进方法有：优化初始条件^[4]、改进模型背景值构造方式^[5]、优化参数估计方法^[6]、构建无偏模型^[7]等。但随着时间的推移，系统结构也在不断变化，若使用恒定参数会使模型拟合结果与预测结果之间的误差变大。王正新^[8]在灰作用量中引入系统延迟时间和时间作用参数，构建适合分析小样本振荡序列的非线性优化模型，削弱了原始序列随机波动对建模的影响。但上述模型是基于等间隔的原始观测序列，而在实际变形监测中很难实现连续等间距观测。近年来，有学者建立非等间距GM(1, 1)幂模型^[9]，但未考虑到背景值等权构造的缺陷，导致优化模型对变形监测原始数据的拟合和预测精度不够理想。此外，现有的非等间距GM(1, 1)幂模型的研究忽略了参数随时间的动态变化特性，使得系统内部的相关变量无法随观测序列的变动而更新。

本文提出一种优化背景值的非等间距线性时变参数GM(1, 1)幂模型。根据GM(1, 1)幂模型的白化微分方程推导模型的时间响应序列，基于积分原理建立灰色微分方程，求解模型的线性参数，并使用粒子群算法优化选取非线性参数，避免出现参数选取时目标函数适应度局部最优的问题，从而有效提高模型的拟合与预测精度。

1 线性时变参数非等间距GM(1, 1)幂模型

对同一序列而言，发展系数a保持不变，灰作用量b会随时间发生变化^[10]。在非等间距GM(1, 1)幂模型中，若灰作用量保持不变，则估计误差会持续增大。由文献[11]可知，将b^*=e^μtb通过级数展开为(1+μt)b，用关于时间t的线性函数(μbt+b)代替原灰作用量b，能够改进非线性系统内部小样本数据序列建模精度不足的问题。由于μ为常数，为方便表示，下文用(bt+c)代替b建立线性时变参数非等间距GM(1, 1)幂模型。

定义1^[12]   设原始非负光滑序列为X⁽⁰⁾={x⁽⁰⁾(t₁), x⁽⁰⁾(t₂), …, x⁽⁰⁾(t_n)}，时序间距Δt_i=t_i－t_i－1≠const, i=2, 3, …, n，Δt₁=1, 称

$ \frac{{{\rm{d}}{x^{(0)}}(t)}}{{{\rm{d}}t}} + a{x^{(0)}}(t) = (bt + c){\left[ {{x^{(0)}}(t)} \right]^\gamma } $

(1)

为线性时变参数非等间距GM(1, 1)幂模型的白化微分方程。

定理1  线性时变参数非等间距GM(1, 1)幂模型的时间响应序列为：

$ {\hat x^{(0)}}\left( {{t_k}} \right) = \left\{ {\begin{array}{*{20}{l}} {{x^{(0)}}\left( {{t_1}} \right), k = 1}\\ {{{\left\{ {\left[ {{{\left[ {{x^{(0)}}} \right]}^{1 - \gamma }} - \frac{b}{a}{t_1} + \frac{b}{{{a^2}(1 - \gamma )}} - \frac{c}{a}} \right]{{\rm{e}}^{a(1 - \gamma )\left( {{t_1} - {t_k}} \right)}} + \frac{b}{a}{t_k} - \frac{b}{{{a^2}(1 - \gamma )}} + \frac{c}{a}} \right\}}^{\frac{1}{{1 - \gamma }}}}, k = 2, \cdots , n} \end{array}} \right. $

(2)

证明  根据伯努利微分方程求解法则，式(1)两端同时除以[x⁽⁰⁾(t)]^γ得：

$ {\left[ {{x^{(0)}}(t)} \right]^{ - \gamma }}\frac{{{\rm{d}}{x^{(0)}}(t)}}{{{\rm{d}}t}} + {\left[ {{x^{(0)}}(t)} \right]^{1 - \gamma }} = bt + c $

(3)

令y⁽⁰⁾(t)=[x⁽⁰⁾(t)]^1－γ，求导有$\frac{{{\rm{d}}{y^{(0)}}(t)}}{{{\rm{d}}t}} = (1 - \gamma ){\left[ {{x^{(0)}}(t)} \right]^{ - \gamma }}\frac{{{\rm{d}}{x^{(0)}}(t)}}{{{\rm{d}}t}}$代入式(3)并在方程两边同时乘以(1－γ)得：

$ {y^{(0)}}(t) = {C_1}{{\rm{e}}^{ - a(1 - \gamma )t}} + \frac{b}{a}t - \frac{b}{{{a^2}(1 - \gamma )}} + \frac{c}{a} $

(4)

式中，C₁为任意常数，由建模初始条件决定。

线性时变参数非等间距GM(1, 1)幂模型的时间响应函数可表示为：

$ {x^{(0)}}(t) = {\left\{ {{C_1}{{\rm{e}}^{ - a(1 - \gamma )t}} + \frac{b}{a}t - \frac{b}{{{a^2}(1 - \gamma )}} + \frac{c}{a}} \right\}^{\frac{1}{{1 - \gamma }}}} $

(5)

时间响应序列为：

$ \begin{array}{*{20}{c}} {{{\hat x}^{(0)}}\left( {{t_k}} \right) = }\\ {\left\{ {\left[ {{{\left[ {{x^{(0)}}\left( {{t_1}} \right)} \right]}^{1 - \gamma }} - \frac{b}{a}{t_1} + \frac{b}{{{a^2}(1 - \gamma )}} - \frac{c}{a}} \right] \cdot } \right.}\\ {{{\left. {{{\rm{e}}^{ - a(1 - \gamma )\left( {{t_k} - {t_1}} \right)}} + \frac{b}{a}{t_k} - \frac{b}{{{a^2}(1 - \gamma )}} + \frac{c}{a}} \right\}}^{\frac{1}{{1 - \gamma }}}}} \end{array} $

(6)

式中，t=t_k, k=2, 3, …, n，定理1得证。

此外，由${\hat x^{(0)}}\left( {{t_1}} \right) = {x^{(0)}}\left( {{t_1}} \right)$可推导参数C₁的解析式为：

$ \begin{gathered} C_{1}= \\ \left\{\left[x^{(0)}\left(t_{1}\right)\right]^{1-\gamma}-\frac{b}{a} t_{1}+\frac{b}{a^{2}(1-\gamma)}-\frac{c}{a}\right\} \mathrm{e}^{-a(\gamma-1) t_{1}} \end{gathered} $

(7)

对于原始数据本身呈近似幂次增长的光滑序列，利用上述模型的时间响应序列即可求得模型拟合值。相比于传统GM(1, 1)幂模型，上述模型无需进行差分还原处理，可简化建模过程。

2 线性时变参数非等间距GM(1, 1)幂模型背景值优化 2.1 模型最优背景值构建

传统GM(1, 1)幂模型背景值计算一般采用公式z⁽¹⁾(t_k)=0.5x⁽¹⁾(t_k－1)+0.5x⁽¹⁾(t_k)，k=2, 3, …, n，即认为前后紧邻值等权，这与白化方程中对应的$\int_{{t_{k - 1}}}^{{t_k}} {{x^{(1)}}} (t){\rm{d}}t$存在差异。在非等间距序列中，Δt_i=t_i－t_i－1≠const，即原始序列不以平均速率变形，传统GM(1, 1)幂模型背景值的构造方式在非等间距序列中可能会导致拟合与预测效果不佳，因此需要根据目标函数对背景值的前后紧邻序列赋以不同的权重。

改进后的背景值表达式为：

$ z^{(1)}\left(t_{k}\right)=\lambda x^{(1)}\left(t_{k-1}\right)+(1-\lambda) x^{(1)}\left(t_{k}\right) $

(8)

式中，λ∈[0, 1], k=2, 3, …, n。

2.2 优化背景值的时变参数非等间距GM(1, 1)幂模型参数估计

定义2  设非负原始光滑序列为：X⁽⁰⁾={x⁽⁰⁾(t₁), x⁽⁰⁾(t₂), …, x⁽⁰⁾(t_n)}，间距Δt_i=t_i－t_i－1≠const，Δt₁=1；X⁽⁰⁾的(1-γ)次幂生成序列为y⁽⁰⁾(t_k)=[x⁽⁰⁾(t_k)]^1－γ，y⁽⁰⁾(t₁)=0，Y⁽⁰⁾的一阶累减生成序列为Y^(－1)={y^(－1)(t₂), y^(－1)(t₃), …, y^(－1)(t_n)}，其中y^(－1)(t_k)=y⁽⁰⁾(t_k)－y⁽⁰⁾(t_k－1)/Δt_k；Z⁽⁰⁾={z⁽⁰⁾(t₂), z⁽⁰⁾(t₃), …, z⁽⁰⁾(t_n)}为Y⁽⁰⁾的紧邻均值序列，z⁽⁰⁾(t_i)=λx⁽⁰⁾(t_k－1)+(1－λ)x⁽⁰⁾(t_k)，其中，λ∈[0, 1], k=2, 3, …, n。

定理2  线性时变参数非等间距GM(1, 1)幂模型的灰微分方程为：

$ \begin{aligned} &y^{(-1)}\left(t_{k}\right) \Delta t_{k}+a(1-\gamma) z^{(0)}\left(t_{i}\right) \Delta t_{k}= \\ &\frac{1}{2} b(1-\gamma)\left(t_{k}^{2}-t_{k-1}^{2}\right)+c(1-\gamma) \Delta t_{k} \end{aligned} $

(9)

证明  式(1)两边同时乘以(1－γ)[x⁽⁰⁾]^－γ得：

$ \begin{gathered} (1-\gamma)\left[x^{(0)}(t)\right]^{-\gamma} \frac{\mathrm{d} x^{(0)}(t)}{\mathrm{d} t}+a(1-\gamma) \cdot \\ {\left[x^{(0)}(t)\right]^{1-\gamma}=(b t+c)(1-\gamma)} \end{gathered} $

(10)

令y⁽⁰⁾(t)=[x⁽⁰⁾]^1－γ，则式(10)可简化为：

$ \frac{\mathrm{d} y^{(0)}(t)}{\mathrm{d} t}+a(1-\gamma) y^{(0)}(t)=(b t+c)(1-\gamma) $

(11)

式(11)两边在区间[t_k－1, t_k]上积分得：

$ \begin{aligned} &y^{(-1)}\left(t_{k}\right) \Delta t_{k}+a(1-\gamma) z^{(0)}\left(t_{i}\right) \Delta t_{k}= \\ &\frac{1}{2} b(1-\gamma)\left(t_{k}^{2}-t_{k-1}^{2}\right)+c(1-\gamma) \Delta t_{k} \end{aligned} $

(12)

即定理2得证。

确定幂指数γ与背景值权重λ后，同样采用最小二乘法估计灰微分方程参数，即

$ {\left[ {\begin{array}{*{20}{l}} a&b&c \end{array}} \right]^{\rm{T}}} = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Y}} $

(13)

式中，

$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {(\gamma - 1){z^{(0)}}\left( {{t_2}} \right)\Delta {t_2}}&{0.5(1 - \gamma )\left( {t_2^2 - t_1^2} \right)}&{(1 - \gamma )\Delta {t_2}}\\ {(\gamma - 1){z^{(0)}}\left( {{t_3}} \right)\Delta {t_3}}&{0.5(1 - \gamma )\left( {t_3^2 - t_2^2} \right)}&{(1 - \gamma )\Delta {t_3}}\\ \vdots & \vdots & \vdots \\ {(\gamma - 1){z^{(0)}}\left( {{t_n}} \right)\Delta {t_n}}&{0.5(1 - \gamma )\left( {t_n^2 - t_{n - 1}^2} \right)}&{(1 - \gamma )\Delta {t_n}} \end{array}} \right], \mathit{\boldsymbol{Y}} = \left[ {\begin{array}{*{20}{c}} {{y^{( - 1)}}\left( {{t_2}} \right)\Delta {t_2}}\\ {{y^{( - 1)}}\left( {{t_3}} \right)\Delta {t_3}}\\ \vdots \\ {{y^{( - 1)}}\left( {{t_n}} \right)\Delta {t_n}} \end{array}} \right] $

(14)

2.3 粒子群算法优化幂指数和背景值权重

粒子群优化算法(particle swarm optimization，PSO)是一种全局优化算法。首先在搜索空间域内初始化一群随机粒子，然后经过不断迭代、调整得到最优解。PSO具有设置参数少、精度高、收敛速度快等优点，被广泛应用于函数优化、计算机大数据和神经网络等众多领域^[13]。最优解的选择方式有目标函数小于给定限差和设定迭代次数2种，给定限差可能会导致粒子陷入局部最优，设定迭代次数可保证每个初始种群粒子调整位置与速度，再通过迭代次数与适应度的收敛关系调整迭代次数，确保粒子达到全局最优。因此，采用设定迭代次数的方式来选取优化背景值的线性时变参数非等间距GM(1, 1)幂模型的最优幂指数γ与背景值权重λ。

目标函数为最小平均绝对百分比误差(MAPE)，即

$ \begin{array}{*{20}{c}} {\min f(\gamma , \lambda ) = \frac{1}{{n - 1}}\sum\limits_{k = 2}^n {\left| {\frac{{{{\hat x}^{(0)}}\left( {{t_k}} \right) - {x^{(0)}}\left( {{t_k}} \right)}}{{{x^{(0)}}\left( {{t_k}} \right)}}} \right|} , }\\ {\gamma \in [ - 3, 3], \lambda \in [0, 1]} \end{array} $

(15)

首先根据PSO的常规参数设定方法设置实验运行参数；然后通过沉降观测序列适应度收敛趋势调整群体规模与迭代次数，加快计算速度；最后确定本文PSO的初始化参数。具体运行参数初始化设置如下：学习因子c₁=2, c₂=2；初始速度V_min=－1, V_max=1；动态惯性权重因子ω=0.8；粒子搜索域γ∈[－3, 3], λ∈[0, 1]；群体规模size=100；进化迭代次数k=1 500；加权速度系数α₁, α₂∈[0, 1]。

模型非线性参数求解步骤如下：

1) 设置运行参数的初始化值，计算初始种群在搜索域内每个粒子的位置X_id¹和速度V_id¹。

2) 根据目标函数计算初始化粒子群各粒子的适应度fitness，寻找初始搜索域内参数γ与λ的个体最佳位置P_best1、P_best2和全局最佳位置G_best1、G_best2。

3) 在种群规模和进化次数的搜索条件下，通过迭代不断更新粒子速度V_id¹和位置X_id¹，即

$ \begin{array}{l} V_{id}^{k + 1} = \omega V_{id}^k + {c_1}{\alpha _1}\left( {{P_{{\rm{best }}}} - X_{id}^k} \right) + \\ {c_2}{\alpha _2}\left( {{G_{{\rm{best }}}} - X_{id}^k} \right)\\ X_{id}^{k + 1} = X_{id}^k + V_{id}^{k + 1} \end{array} $

(16)

4) 若某粒子个体位置适应度优于P_best1、P_best2，则设置新位置为P_best1、P_best2；若粒子全局位置适应度优于G_best1、G_best2，则设置新位置为G_best1、G_best2。

5) 输出运行结果G_best1为幂指数γ的最优解；参数G_best2为背景值权重λ的最优解。此时优化背景值的线性时变非等间距GM(1, 1)幂模型的MAPE最小。

2.4 模型精度验证

以MAPE为检验指标，比较不同优化角度下GM(1, 1)幂模型的拟合与预测精度。设原始序列观测值为X⁽⁰⁾={x⁽⁰⁾(t₁), x⁽⁰⁾(t₂), …, x⁽⁰⁾(t_n)}；模型求得原始序列的拟合值为${\hat X^{(0)}} = \left\{ {{{\hat x}^{(0)}}\left( {{t_1}} \right), } \right.\left. {{{\hat x}^{(0)}}\left( {{t_2}} \right), \cdots , {{\hat x}^{(0)}}\left( {{t_n}} \right)} \right\}$，则残差序列为：

$ \varepsilon\left(k_{k}\right)=x^{(0)}\left(k_{k}\right)-\hat{x}^{(0)}\left(k_{k}\right), k=2, 3, \cdots, n $

(17)

相对误差序列为：

$ P\left(k_{k}\right)=\frac{\varepsilon\left(k_{k}\right)}{x^{(0)}\left(k_{k}\right)} \times 100 \%, k=2, 3, \cdots, n $

(18)

MAPE为：

$ {\rm{MAPE}} = \frac{1}{{n - 1}}\sum\limits_{k = 2}^n {\left| {P\left( {{k_k}} \right)} \right|} $

(19)

3 工程应用及结果验证

为证明优化背景值的非等间距线性时变参数GM(1, 1)幂模型的适用性与优越性，利用2组矿区变形监测数据，分别建立非等间距幂指数优化模型^[3](模型1)、初始值加权平均的非等间距GM(1, 1)幂模型^[14](模型2)、线性时变参数GM (1, 1)幂模型^[15] (模型3)与本文模型进行对比分析。

3.1 宁东矿区GNSS沉降监测工程

以宁东矿区GNSS持续沉降监测数据为例，首先选取金凤20025监测站点前7期累积沉降量作为建模数据，然后利用第8期、第9期数据对模型进行预测结果验证，最后比较本文模型与其他3种幂模型的拟合精度和预测精度(表 1)。借助MATLAB软件得到模型参数：a=0.157 981、b=0.028 582、c=4.300 362、γ=-0.067 794、λ=0.777 521，非线性参数与MAPE的关系如图 1所示。对应的优化背景值的线性时变非等间距GM(1, 1)幂模型时间响应关系式为：

$ \begin{array}{*{20}{c}} {{{\hat x}^{(0)}}\left( {{t_k}} \right) = \left\{ { - 13.642\;18{{\rm{e}}^{0.168691{t_k}}} + } \right.}\\ {{{\left. {0.180\;92{t_k} + 26.148\;27} \right\}}^{\frac{1}{{1.067794}}}}, k = 2, 3, \cdots , n} \end{array} $

(20)

宁东矿区GNSS沉降监测工程计算结果表明，模型1、模型2和模型3拟合数值的MAPE分别为3.31%、8.62%、4.18%；预测数值的MAPE分别为9.97%、4.70%、3.75%。优化背景值的非等间距线性时变参数GM(1, 1)幂模型拟合与预测的MAPE分别为2.33%、2.10%，相比于上述3种优化模型，本文模型的精度更高，拟合值和预测值与原始观测序列更加吻合。

图 1可反映PSO的寻优过程，由图可见，在该沉降监测工程中，随着幂指数γ在搜索区间内逐渐靠近0，背景值权重λ在搜索区间内不断接近0.7，MAPE持续减小。当幂指数为-0.067 794、背景值权重为0.777 521时，MAPE取得最小值2.33%。

3.2 梁家矿铁路岩移沉降监测工程

以梁家矿4112回采工作面铁路岩移沉降数据为例，利用前8期原始沉降数据构建模型，后2期数据进行模型预测。4种模型的拟合值和预测值精度比较见表 2。本文模型参数分别为：a=0.327 970、b=0.006 893、c=1.681 693、γ=0.465 813、λ=0.792 941，非线性参数优化选取与MAPE的关系如图 2所示。对应的优化背景值的线性时变非等间距GM(1, 1)幂模型时间响应关系式为：

$ \begin{array}{*{20}{c}} {{{\hat x}^{(0)}}\left( {{t_k}} \right) = \left\{ { - 2.438\;10{{\rm{e}}^{0.175198{t_k}}} + } \right.}\\ {{{\left. {0.021\;02{t_k} + 5.007\;57} \right\}}^{\frac{1}{{0.534187}}}}, k = 2, 3, \cdots , n} \end{array} $

(21)

由表 2可见，在梁家矿铁路岩移沉降监测工程中，优化背景值的非等间距线性时变参数GM(1, 1)幂模型拟合值与预测值的MAPE分别为4.70%和6.38%，拟合值精度较模型3提高6.37%，预测值精度提高2.63%，进一步体现出本文模型在沉降监测工程建模与分析中的优越性。

由图 2可知，在搜索阈γ∈[－3, 3], λ∈[0, 1]内，背景值权重λ的收敛速度较幂指数γ更快，当γ=0.465 813、λ=0.792 941时，MAPE取到最小值4.70%。

4 结语

1) 针对沉降监测小样本非等间距原始数据序列建模问题，构建优化背景值的非等间距线性时变参数GM(1, 1)幂模型。基于粒子群算法，以最小MAPE为目标函数对非线性参数进行优化选取，可有效避免局部极小值陷阱，增强非线性参数选取的可靠性。

2) 以2组煤矿开采变形监测数据为例，对4种模型分别进行拟合与预测。计算结果表明，优化背景值的非等间距线性时变参数GM(1, 1) 幂模型在宁东矿区开采沉陷GNSS监测工程和梁家矿铁路岩移沉降监测工程的MAPE分别为2.33%、4.70%，预测误差分别为2.10%、6.38%，均优于其他3种常用的非等间距GM(1, 1)优化模型。

3) 相比于传统方法，优化背景值的非等间距线性时变参数GM(1, 1)幂模型可提升变形监测数据的拟合与预测精度，适用于短期沉降监测数据的预测分析，可拓宽GM(1, 1)幂模型在沉降监测数据处理中的应用范围。