灰色模型自提出以来，在理论研究和实践应用方面都得到了发展。灰色系统理论中目前应用最广泛的为GM(1, 1)模型，其具有所需样本少、原理简单、运算简便等优点。但由于其固有的缺陷，模型的建模精度存在较大偏差，稳定性也不足^[1]。很多学者对GM(1, 1)模型进行了改进，主要方式包括构建无偏模型^[2]、改进背景值的构造方式^[3-4]、调整或修改模型的边界条件^[5]、优化灰作用量^[6]、对原始数据进行函数变换^[7]、优化模型参数估计方法^[8]以及优化残差序列^[9]等。以上优化模型中，GM(1, 1)模型建模的主要方式为对原始数据序列进行一阶累加，求解模型，再一阶累减还原得到模型的拟合和预测数据。由于在实际问题中原始数据间存在不等价随机性的问题，吴利丰^[10]提出分数阶算子灰色预测模型，通过分数阶精确调节累加数之间的数量级，为不同的原始数据序列选择特定的分数阶算子构建灰色预测模型，能明显提高模型对数据的拟合和预测效果。

目前，灰色预测模型在变形监测领域中得到非常广泛的应用^[11]，但是基于分数阶建模思想构建分数阶算子的灰色预测模型在变形监测中的研究和应用还较少。鉴于此，本文构建了粒子群算法优化的分数阶算子EGM(1, 1)模型，通过典型的变形监测算例对比分析得出，优化模型对变形监测数据具有较高的拟合和预测效果，具有实际应用价值。

1 分数阶算子EGM(1, 1)模型

设X⁽⁰⁾={x⁽⁰⁾(1), x⁽⁰⁾(2), …, x⁽⁰⁾(n)}为原始数据序列，r∈R⁺，X^(r)={x^(r)(1), x^(r)(2), …, x^(r)(n)}是X⁽⁰⁾={x⁽⁰⁾(1), x⁽⁰⁾(2), …, x⁽⁰⁾(n)}的r阶累加生成算子，其中，${X^{(r)}}(k) = \sum\limits_{i = 1}^k {\frac{{\varGamma (r + k - i)}}{{\varGamma (k - i + 1)\Gamma (r)}}{x^{(0)}}(i), {\rm{ }}k = 1, 2, \cdots , n} $。z^(r)(k)=(z^(r)(2), z^(r)(3), …, z^(r)(n)), 其中, ${z^{(r)}}(k) = \frac{{{x^{(r)}}(k) + {x^{(r)}}(k - 1)}}{2}, {\rm{ }}k = 2, 3, \cdots , n$。X^(－r)=(x^(－r)(1), x^(－r)(2), …, x^(－r)(n))为X⁽⁰⁾的r阶累减生成算子，其中，$\begin{array}{l} {X^{({\rm{ - }}r)}}(k) = \sum\limits_{i = 1}^{k{\rm{ - }}1} {{{( - 1)}^i}\frac{{\varGamma (r + 1)}}{{\varGamma (i + 1)\varGamma (r - i + 1)}}{x^{(0)}}(k - i)} , {\rm{ }}k = 1, 2, \cdots , n \end{array}$。

分数阶算子EGM(1, 1)模型x^(r－1)(k)+az^(r)(k)=b中的参数向量$\boldsymbol{\hat a} = {[a, b]^{\rm T}}$可以运用最小二乘法估算：

$ \boldsymbol{\hat a} = {({\boldsymbol{B}^{\rm T}}\boldsymbol{B})^{ - 1}}{\boldsymbol{B}^{\rm T}}\boldsymbol{Y} $

(1)

式中，

$ \boldsymbol{Y} = \left[ {\begin{array}{*{20}{c}} {{x^{(r - 1)}}(2)}\\ {{x^{(r - 1)}}(3)}\\ \vdots \\ {{x^{(r - 1)}}(n)} \end{array}} \right]{\rm{ , }} \boldsymbol{B} = \left[ {\begin{array}{*{20}{c}} { - {z^{(r)}}(2)}&1\\ { - {z^{(r)}}(3)}&1\\ \vdots & \vdots \\ { - {z^{(r)}}(n)}&1 \end{array}} \right] $

(2)

因

$ \begin{array}{c} {x^{(r - 1)}}(k) = {x^{(r)}}(k) - {x^{(r)}}(k - 1) =\\ \sum\limits_{i = 1}^k {\frac{{\varGamma (r + k - 1)}}{{\varGamma (k - i + 1)\varGamma (r)}}{x^{(0)}}(i)} - \\ \sum\limits_{i = 1}^{k - i} {\frac{{\varGamma (r + k - i - 1)}}{{\varGamma (k - i)\varGamma (r)}}{x^{(0)}}(i)} , k = 2, 3, \cdots , n \end{array} $

则:

$ {z^{(r)}}(k) = \frac{{\sum\limits_{i = 1}^k {\frac{{\varGamma (r + k - i)}}{{\varGamma (k - i + 1)\varGamma (r)}}{x^{(0)}}(i) + \sum\limits_{i = 1}^{k - 1} {\frac{{\varGamma (r + k - i)}}{{\varGamma (k - i + 1)\varGamma (r)}}{x^{(0)}}(i)} } }}{2} $

(3)

$ \boldsymbol{Y} = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^2 {\frac{{\varGamma (r + 2 - i)}}{{\varGamma (2 - i + 1)\varGamma (r)}}{x^{(0)}}(i) - \sum\limits_{i = 1}^1 {\frac{{\varGamma (r + 2 - 1 - i)}}{{\varGamma (2 - i)\varGamma (r)}}{x^{(0)}}(i)} } }\\ {\sum\limits_{i = 1}^3 {\frac{{\varGamma (r + 3 - i)}}{{\varGamma (3 - i + 1)\varGamma (r)}}{x^{(0)}}(i) - \sum\limits_{i = 1}^2 {\frac{{\varGamma (r + 3 - 1 - i)}}{{\varGamma (3 - i)\varGamma (r)}}{x^{(0)}}(i)} } }\\ \vdots \\ {\sum\limits_{i = 1}^n {\frac{{\varGamma (r + n - i)}}{{\varGamma (n - i + 1)\varGamma (r)}}{x^{(0)}}(i) - \sum\limits_{i = 1}^{n - 1} {\frac{{\varGamma (r + n - 1 - i)}}{{\varGamma (n - i)\varGamma (r)}}{x^{(0)}}(i)} } } \end{array}} \right] $

$ \boldsymbol{B} = \left[ {\begin{array}{*{20}{c}} { - \frac{1}{2}[\sum\limits_{i = 1}^2 {\frac{{\varGamma (r + 2 - i)}}{{\varGamma (2 - i + 1)\varGamma (r)}}{x^{(0)}}(i) + \sum\limits_{i = 1}^1 {\frac{{\varGamma (r + 1 - i)}}{{\varGamma (2 - i)\varGamma (r)}}{x^{(0)}}(i)} } ]}&1\\ { - \frac{1}{2}[\sum\limits_{i = 1}^3 {\frac{{\varGamma (r + 3 - i)}}{{\varGamma (3 - i + 1)\varGamma (r)}}{x^{(0)}}(i) + \sum\limits_{i = 1}^2 {\frac{{\varGamma (r + 2 - i)}}{{\varGamma (3 - i)\varGamma (r)}}{x^{(0)}}(i)} } ]}&1\\ \vdots & \vdots \\ { - \frac{1}{2}[\sum\limits_{i = 1}^n {\frac{{\varGamma (r + n - i)}}{{\varGamma (n - i + 1)\varGamma (r)}}{x^{(0)}}(i) + \sum\limits_{i = 1}^{n - 1} {\frac{{\varGamma (r + n - i - 1)}}{{\varGamma (n - i)\varGamma (r)}}{x^{(0)}}(i)} } ]}&1 \end{array}} \right] $

式(4)为分数阶算子EGM(1, 1)模型x^(r－1)(k)+az^(r)(k)=b的白化微分方程：

$ \frac{{{\rm d}{x^{(r)}}}}{{{\rm d}t}} + a{x^{(r)}} = b $

(4)

设$\boldsymbol{Y}、\boldsymbol{B} 、\boldsymbol{\hat a}$如上文所述，$ \boldsymbol{\hat a} = {[a, b]^{\rm T}} = {({ \boldsymbol{B} ^{\rm T}} \boldsymbol{B})^{ - 1}}{ \boldsymbol{B}^{\rm T}} \boldsymbol{Y}$，则分数阶算子EGM(1, 1)模型的白化微分方程$\frac{{{\rm d}{x^{(r)}}}}{{{\rm d}t}} + a{x^{(r)}} = b$的时间响应函数为：

$ {\hat x^{(r)}}(t) = [{x^{(r)}}(1) - \frac{b}{a}]{{\rm e}^{ - at}} + \frac{b}{a} $

(5)

分数阶算子EGM(1, 1)模型x^(r－1)(k)+az^(r)(k)=b的时间响应序列为：

$ \begin{array}{c} {\hat x^{(r)}}(k) = [{x^{(0)}}(1) - \frac{b}{a}]{{\rm e}^{ - a(k - 1)}} + \frac{b}{a} \\ k = 2, 3, \cdots , n \end{array} $

(6)

模型还原值：

$ \begin{array}{l} {{\hat x}^{(0)}}(k) = {({{\hat x}^{(r)}}){( - r)}}(k)=\\ \sum\limits_{i = 0}^{k - 1} {{{( - 1)}^i}} \frac{{\varGamma (r + 1)}}{{\varGamma (i + 1)\varGamma (r - i + 1)}}\\{{\hat x}^{(r)}}(k - i) , k = 2, 3, \cdots , n\\ {\hat x^{(0)}}(1) = {x^{(0)}}(1) \end{array} $

(7)

2 粒子群算法选择EGM(1, 1)模型最优分数阶次

粒子群(particle swarm optimization, PSO)算法具有概念易理解、调整参数少、编程易实现等优点，在神经网络训练和函数优化等领域得到广泛的应用^[12]。

PSO算法的基本思路为：在搜索空间中初始化为一群随机粒子(随机解)，然后通过迭代找到最优解。每一次迭代中，所有粒子都有一个由被优化的函数决定的适应度。基于适应度，粒子本身找到的最优解称为个体极值(P_best), 整个种群目前找到的最优解为全局极值(G_best)。在找到这2个最优值时，每个粒子基于个体和全局极值更新其速度和位置。每个粒子的速度和位置的更新表示为：

$ V_{id}^{K + 1} = wV_{id}^k + {c_1}{\alpha _1}(P_{\rm best} - X_{id}^k) + \\ {c_2}{\alpha _2}(G_{\rm best} - X_{id}^k),X_{id}^{k + 1} = X_{id}^k + V_{id}^{k + 1} $

式中，w为惯性权重因子，决定粒子先前速度对当前速度的影响程度，起到平衡算法全局搜索和局部搜索的作用，c₁和c₂为控制粒子速度的加速因子，随机取2左右的值，α₁和α₂为随机产生的在[0, 1]之间变化的加速度权重系数。每个粒子跟踪其先前的最佳位置和全局最佳位置，不断更新其速度和位置，直到迭代次数超过最大值。

最小平均相对误差下对分数阶算子EGM(1, 1)模型的最优分数阶次最优化问题为：

$ \begin{array}{c} \min f(r) = \\ \frac{1}{{n - 1}}\sum\limits_{k = 2}^n {\left| {\frac{{{{\hat x}^{(0)}}(k) - {x^{(0)}}(k)}}{{{x^{(0)}}(k)}}} \right|}, {\rm{ }}r \in {R^ + } \end{array} $

(8)

本文粒子群算法的运行参数设定如下：学习因子c₁、c₂=2，动态惯性权重因子w=0.8，算法最大迭代次数为3 000，微粒范围r∈[0, 1]，群体个体数目为50，精度设定为0.000 001。

粒子群算法优化的分数阶EGM(1, 1)模型的建模求解步骤如下：

1) 随机初始化粒子群中每个粒子的位置X_id¹和速度V_id¹, 可取P_best=1，即EGM(1, 1)模型；

2) 将粒子中的P_best设置为当前位置，G_best设置为初始群体中最佳粒子的位置；

3) 计算模型分数阶r=P_best时的分数阶算子EGM(1, 1)模型的拟合平均相对误差；

4) 更新粒子群中每个粒子的V_id¹和X_id¹；

5) 对每个粒子，如果其粒子适应度优于P_best的适应度，将P_best设置为新位置；如果粒子适应度优于G_best的适应度，将G_best设置为新位置；

6) 输出G_best，即r的最优值，分数阶算子EGM(1, 1)模型的分数阶算子为r的最优值时，构建的模型拟合原始数据序列的平均相对误差值最小。

3 模型精度检验

由以上模型可以求出原始序列的拟合值为：

$ {\hat X^{(0)}}(k) = \left\{ {{{\hat x}^{(0)}}(1), {{\hat x}^{(0)}}(2), \cdots , {{\hat x}^{(0)}}(n)} \right\} $

(9)

令$\varepsilon (k) = {x^{(0)}}(k) - {\hat x^{(0)}}(k), {\rm{ }}k = 2, \cdots , n$。残差可表示为：

$ E(k) = \left\{ {\varepsilon (2), \varepsilon (3), \cdots , \varepsilon (n)} \right\} $

(10)

相对误差序列为：

$ {\rm PE}({T_k}) = \frac{{\varepsilon (k)}}{{{x^{(0)}}(k)}} \times 100\% , {\rm{ }}k = 2, 3, \cdots , n $

(11)

平均相对误差为：

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

(12)

4 优化模型实例分析应用 4.1 算例1

以新建的呼和浩特至准格尔铁路工程的沉降观测数据为例。以DK41+550号沉降监测点2015-07-02~13的12期累计沉降量数据作为建模数据，分别构建EGM(1, 1)模型和优化分数阶的EGM(1, 1)模型，并用构建的模型预测2015-07-14~17的4期沉降情况。2种模型的拟合和预测精度如表 1所示。

EGM(1, 1)模型的时间响应序列为：

$ \begin{array}{c} {{\hat x}^{(1)}}(k) = - 106.207{{\rm e}^{0.063(k - 1)}} - 112.157\\ k = 2, 3, \cdots , n \end{array} $

(13)

优化分数阶EGM(1, 1)模型(r=0.480 492)的时间响应序列为：

$ \begin{array}{c} {{\hat x}^{(0.480 \;492)}}(k) = - 268.466{{\rm e}^{ - 0.013(k - 1)}} - 274.416\\ k = 2, 3, \cdots , n \end{array} $

(14)

从结果可以看出, EGM(1, 1)模型的最优阶次为0.480 492。优化模型的拟合和预测结果的精度相比EGM(1, 1)模型有了很大程度的提高，其构建的模型更符合原始数据数列的发展规律。

4.2 算例2

向家坡滑坡位于渝黔高速公路K13+500~960段。以文献[13]中滑坡监测点JB5的监测值为例，对其空间位移进行建模和预测。用前7期数据构建模型，预测后2期滑坡的位移量。精度评价结果如表 2所示。

EGM(1, 1)模型的时间响应序列为：

$ \begin{array}{c} {{\hat x}^{(1)}}(k) = 168.694{{\rm e}^{0.102(k - 1)}} - 136.394\\ k = 2, 3, \cdots , n \end{array} $

(15)

优化分数阶EGM(1, 1)模型(r=0.244 999)的时间响应序列为：

$ \begin{array}{c} {{\hat x}^{(0.244\;999)}}(k) = 123.777{{\rm e}^{0.046(k - 1)}} - 156.077\\ k = 2,3, \cdots ,n \end{array} $

(16)

可以看出, 优化分数阶的EGM(1, 1)模型的拟合和预测精度明显优于EGM(1, 1)模型。

5 结语

最优阶次的选取是分数阶算子灰色预测模型的关键。本文通过粒子群算法选择拟合原始变形监测数据平均相对误差最小的分数阶次，构建粒子群算法优化的分数阶算子EGM(1, 1)模型。典型的变形监测算例分析表明，相比传统的灰色预测模型，优化模型对变形监测数据的拟合效果和预测精度更优，在变形监测领域中具有可行性和应用价值。