尾矿库是存储固液混合态矿产废渣、具有高势能的人造泥石流危险源，其安全程度严重影响社会经济与生态环境^[1].而浸润线作为尾矿库的生命线，其位置高低与尾矿库安全密切相关.准确掌握浸润线高度，对确保尾矿库正常运行至关重要.然而浸润线高度受多种因素影响，致使其高度呈非线性变化，为准确、有效预测浸润线高度及变化趋势带来困难^[2].鉴于此，李全明等^[3]利用模糊理论分析尾矿库失事概率，开发风险评估软件预测浸润线高度变化，有效决策坝体安全稳定性，但预测精度尚待提高.李刚^[4]用灰色理论的方法研究坝体变形数据，建立坝体位移预测模型，预测坝体形变趋势，效果良好，但只是基于仿真，并未经过工程验证，而且预测精度有待提高.利用灰度模型GM(grey mode)针对小样本、贫信息数据的优势，结合BP网络非线性处理能力强的特点，建立了尾矿库浸润线预测模型，该模型仿真相对传统方法预测精度较高.

1 GM(1, 1)模型算法局限及优化

灰色系统理论是由我国学者邓聚龙教授于1982年提出的^[5]，是将离散随机数变为较有规律的生成数，该生成数随机性被显著衰弱，建立起的微分方程形式的模型，这样便于对其变化过程进行研究和描述.

1.1 GM(1, 1)模型算法

GM(1, 1)模型由一个单变量的一阶微分方程构成，用于复杂系统内某一主导因素特征值拟合和预测. GM(1, 1)模型为最基本GM，推导过程见文献[6].

设X⁽⁰⁾为原始正数据序列：

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

(1)

设X⁽¹⁾为X⁽⁰⁾的1-AGO序列：

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

(2)

其中${x^{(1)}}(k) = \sum\limits_{i = 1}^k {{x^{(0)}}} (i)$, k=1, 2, 3, …, n.

设Z⁽¹⁾是X⁽¹⁾的紧邻均值序列：

${Z^{(1)}} = \left( {{z^{(1)}}(2), {z^{(1)}}(3), {z^{(1)}}(4), \cdots , {z^{(1)}}(n)} \right) $

(3)

其中${z^{(1)}}(k) = \frac{1}{2}\left({{x^{(1)}}(k) + {x^{(1)}}(k - 1)} \right)$，k=2, 3, 4, …, n.

定义GM(1, 1)的基本形式为

${x^{(0)}}(k) + a{z^{(1)}}(k) = b $

(4)

其中：a、b为辨识参数，a为发展系数，反映待预测量${\hat x^{(1)}}、{\hat x^{(0)}}$的变化趋势；b为灰色作用量，反映数据间的变化关系.

设$\hat a$为待估参数，$\hat a = {(a, b)^{\rm{T}}}$，则

$B\hat a = {Y_n} $

由最小二乘估计法求$\hat a$的解为

$\hat a = {\left( {{B^{\rm{T}}}B} \right)^{ - 1}}{B^{\rm{T}}}{Y_n} $

(5)

其中：$B = \left[{\begin{array}{*{20}{c}} { - {z^{(1)}}(2)}&1\\ { - {z^{(1)}}(3)}&1\\ \vdots & \vdots \\ { - {z^{(1)}}(n)}&1 \end{array}} \right]$，${Y_n} = \left[{\begin{array}{*{20}{c}} {{x^{(0)}}(2)}\\ {{x^{(0)}}(3)}\\ \vdots \\ {{x^{(0)}}(n)} \end{array}} \right]$.

白化GM(1, 1)模型，其形式为

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

(6)

式(4)的时间响应序列为

${\hat x^{(1)}}(k + 1) = \left( {{x^{(0)}}(1) - \frac{b}{a}} \right){{\rm{e}}^{ - ak}} + \frac{b}{a} $

(7)

对式(7)做一次累减(1-IAGO)，得到原始序列的预测值为

$\begin{array}{*{20}{c}} {{{\hat x}^{(0)}}(k + 1) = {{\hat x}^{(1)}}(k + 1) - {{\hat x}^{(1)}}(k) = }\\ {\left( {1 - {{\rm{e}}^a}} \right)\left( {{x^{(0)}}(1) - \frac{b}{a}} \right){{\rm{e}}^{ - ak}}} \end{array} $

(8)

其中：${\hat X^{(0)}}$=(${\hat x^{(0)}}$(1), ${\hat x^{(0)}}$(2), ${\hat x^{(0)}}$(3), …, ${\hat x^{(0)}}$(k))，以此得到原始序列的预测值.

1.2 GM(1, 1)模型算法局限

实践发现，GM(1, 1)模型的拟合或预测效果参差不齐，通过分析GM(1, 1)模型发现: 1)灰色预测模型的预测精度与被预测对象的递变规律以及数据序列的光滑度有关；2)灰色微分拟合法建立的离散拟合方程是一个近似差分方程，很难保证拟合方程与待拟合系统的微分方程严格近似.

由式(5)知，辨识参数的求解与背景值Z⁽¹⁾密切相关.而原始定义中的背景值由紧邻均值得到，也可以采用权系数的方式求背景值.所以Z⁽¹⁾求取方式不同时，辨识参数a、b值的大小就会有所变化.准确合理地a、b值能够提高模型序列的预测精度^[6]，因此，提出一种回避背景值的辨识参数求解方法，避免求解a、b时对背景值试算选取的步骤以及选取不当造成模型预测精度低的问题.

除此之外，由式(8)可知，原始序列预测值是在初始条件为${\hat x^{(1)}}$(1)=x⁽¹⁾(1)=x⁽⁰⁾(1)下得到的，并且拟合曲线必过(1, x⁽⁰⁾(1))点.而原始定义中该拟合曲线不一定通过此点，这使得将${\hat x^{(1)}}$(1)=x⁽¹⁾(1)=x⁽⁰⁾(1)作为初始条件是不合理的^[7]，因此，提出用λx⁽⁰⁾(1)来代替x⁽⁰⁾(1)作为新的初始值，避免求解时初始值为固定值的不足，所以只要求出λ值就可以确定新初始值，避免预测精度受固定初始值的影响.

1.3 GM(1, 1)模型算法优化

GM(1, 1)模型的优化过程主要分为以下2步：首先采用改进的欧拉公式回避背景值的求解，以得到辨识参数；再采用自适应构造法构造初始值，克服初始条件为原序列初始值时为模型预测带来预测偏差大的问题.

1) 结合改进的欧拉公式回避背景值求解辨识参数，过程如下.

先用欧拉公式求初步预测值

${\bar y_{n + 1}} = {y_n} + hf\left( {{x_n}, {y_n}} \right) $

(9)

其中: h为步长，h=x_n+1-x_n；y_n+1的精度可能不够，采用梯形公式加以校正，梯形公式为

${y_{n + 1}} = {y_n} + \frac{h}{2}\left[ {f\left( {{x_n}, {y_n}} \right) + f\left( {{x_{n + 1}}, {y_{n + 1}}} \right)} \right] $

(10)

将式(10)代入式(11)，

$\left. \begin{array}{l} y_{n + 1}^{(0)} = {y_n} + hf\left( {{x_n}, {y_n}} \right)\\ y_{n + 1}^{(k + 1)} = {y_n} + \frac{h}{2}\left[ {f\left( {{x_n}, {y_n}} \right) + f\left( {{x_{n + 1}}, y_{n + 1}^{(k)}} \right)} \right], \\ \;\;\;\;\;\;\;k = 0, 1, 2, \cdots \end{array} \right\} $

(11)

得到校正值，即

${y_{n + 1}} = {y_n} + \frac{h}{2}\left[ {f\left( {{x_n}, {y_n}} \right) + f\left( {{x_{n + 1}}, {{\bar y}_{n + 1}}} \right)} \right] $

(12)

此时建立具有预测-校正功能的改进欧拉公式：

$\left. \begin{array}{l} 预测:{y_{n + 1}} = {y_n} + hf\left( {{x_n}, {y_n}} \right)\\ 校正:{y_{n + 1}} = {y_n} + \frac{h}{2}\left[ {f\left( {{x_n}, {y_n}} \right) + f\left( {{x_{n + 1}}, {{\bar y}_{n + 1}}} \right)} \right] \end{array} \right\} $

(13)

将其改写为平均化形式：

$\left. \begin{array}{l} {y_p} = {y_n} + hf\left( {{x_n}, {y_n}} \right)\\ {y_c} = {y_n} + hf\left( {{x_{n + 1}}, {y_p}} \right)\\ {y_{n + 1}} = \frac{1}{2}\left( {{y_p} + {y_c}} \right) \end{array} \right\} $

(14)

变换式(6)，即

$\frac{{{\rm{d}}{x^{(1)}}(t)}}{{{\rm{d}}t}} = b - a{x^{(1)}}(t) $

(15)

令

$f\left( {t, {x^{(1)}}(t)} \right) = \frac{{{\rm{d}}{x^{(1)}}(t)}}{{{\rm{d}}t}} = b - a{x^{(1)}}(t) $

(16)

将(16)代入(14)，得

$\begin{array}{*{20}{c}} {{x^{(1)}}(t + 1) + \left( {2 - {a^2}{h^2} + ah} \right) = }\\ {{x^{(1)}}(t)(2 - ah) + 2bh - ab{h^2}} \end{array} $

(17)

因为x⁽¹⁾(t+1)=x⁽¹⁾(t)+x⁽⁰⁾(t+1)，则

${x^{(0)}}(t + 1) = - {x^{(1)}}(t)\frac{{ah}}{{1 + ah}} + \frac{{bh}}{{1 + ah}} $

(18)

令

$M = - \frac{{ah}}{{1 + ah}}, N = \frac{{bh}}{{1 + ah}} $

(19)

取h=1，所以只要求出M、N的值，就能求出a、b.将式(18)表示为

$C\hat b = Y $

(20)

令

$C = \left[ {\begin{array}{*{20}{c}} {{x^{(1)}}(1)}&1\\ {{x^{(1)}}(2)}&1\\ \vdots & \vdots \\ {{x^{(1)}}(n - 1)}&1 \end{array}} \right], \quad Y = \left[ {\begin{array}{*{20}{c}} {{x^{(0)}}(2)}\\ {{x^{(0)}}(3)}\\ \vdots \\ {{x^{(0)}}(n)} \end{array}} \right] $

通过最小二乘法得到M，N，即

$\hat b = {\left( {{C^{\rm{T}}}C} \right)^{ - 1}}{C^{\rm{T}}}Y $

(21)

由于C中只涉及X⁽¹⁾，Y只涉及原始序列，只要求出M、N，再将M、N代入式(19)，就可以得到a、b.此过程避开了背景值的求解得到了辨识参数.至此，回避背景值求解辨识参数的过程介绍完毕.

2) 提出构造自适应初始值

提出用λx⁽⁰⁾(1)代替x⁽⁰⁾(1)，将初始条件设为

$ \hat{x}^{(1)}(1)=x^{(1)}(1)=\lambda x^{(0)}(1) $

(22)

将式(22)代入式(7)，得

$ \hat{x}^{(1)}(1)=x^{(1)}(1)=\lambda x^{(0)}(1) $

(23)

将指标函数定义为

$\begin{array}{*{20}{c}} {W = \sum\limits_{k = 1}^n {{{\left[ {{{\hat x}^{(1)}}(k) - {x^{(1)}}(k)} \right]}^2}} = }\\ {\sum\limits_{k = 1}^n {{{\left\{ {\left[ {\lambda {x^{(0)}}(1) - \frac{b}{a}} \right]{{\rm{e}}^{ - a(k - 1)}} + \frac{b}{a} - {x^{(1)}}(k)} \right\}}^2}} } \end{array} $

(24)

求取λ值时，将式(24)看作是W关于λ的函数.由于W为预测序列误差平方和，只要求得使W取最小值时的λ值，就能满足误差最小的目标，得到较小误差的预测值，过程如下：

$\begin{array}{*{20}{c}} {\frac{{{\rm{d}}W}}{{{\rm{d}}\lambda }} = 2\sum\limits_{k = 1}^n {\left\{ {\left[ {\lambda {x^{(0)}}(1) - \frac{b}{a}} \right]{{\rm{e}}^{ - a(k - 1)}} + } \right.} }\\ {\frac{b}{a} - {x^{(1)}}(k)\} {x^{(0)}}(1){{\rm{e}}^{ - a(k - 1)}} = }\\ {2\left\{ {\lambda {{\left[ {{x^{(0)}}(1)} \right]}^2}\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2a(k - 1)}}} + } \right.}\\ {\begin{array}{*{20}{c}} {\frac{b}{a}{x^{(0)}}(1)\sum\limits_{k = 1}^n {{{\rm{e}}^{ - a(k - 1)}}} - }\\ {{x^{(0)}}(1)\sum\limits_{k = 1}^n {{x^{(1)}}} (k){{\rm{e}}^{ - a(k - 1)}} - }\\ {\frac{b}{a}{x^{(0)}}(1)\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2a(k - 1)}}} \} }\\ {} \end{array}} \end{array} $

(25)

设$\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2a(k - 1)}}} = A$，$\sum\limits_{k = 1}^n {{e^{ - a(k - 1)}}} = B$，$\sum\limits_{k = 1}^n {{x^{(1)}}} (k){{\rm{e}}^{ - a(k - 1)}} = C$

令$\frac{{{\rm{d}}W}}{{{\rm{d}}\lambda }} = 0$，则

$\lambda A{x^{(0)}}(1) + \frac{b}{a}(B - A) - C = 0 $

(26)

解之，得

$\lambda = \frac{{\frac{b}{a}(A - B) + C}}{{A{x^{(0)}}(1)}} $

(27)

W对λ的二阶导数为

$\frac{{{{\rm{d}}^2}W}}{{{\rm{d}}{\lambda ^2}}} = 2{\left[ {{x^{(0)}}(1)} \right]^2}\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2a(k - 1)}}} > 0 $

(28)

由于式(28)恒大于0，式(26)成立时，W在式(27)下取得极小值，此时的λ值即为所求值，λx⁽⁰⁾(1)就可作式(8)中的初始值.

此时，将a、b及λ值代入式(8)中，得

$\begin{array}{*{20}{c}} {{{\hat x}^{(0)}}(k + 1) = {{\hat x}^{(1)}}(k + 1) - {{\hat x}^{(1)}}(k) = }\\ {\left( {1 - {e^a}} \right)\left( {\lambda {x^{(0)}}(1) - \frac{b}{a}} \right){{\rm{e}}^{ - ak}}} \end{array} $

(29)

求解式(29)，当k≤n，${\hat x^{(0)}}(k)$为原始模型序列拟合值；当k>n，${\hat x^{(0)}}(k)$为原始模型序列预测值.

1.4 GM(1, 1)模型算法优化前后比较

将山东某尾矿库2号剖面3号监测点2015-06-15—2015-06-24每天10:00的浸润线高度作为样本数据(共10个)(见表 1)，可以看出，采用振弦式渗压计自动监测的浸润线高度数据与人工监测数据较接近，可作为安全评价的样本数据.

1) 将表 1中前8组数据作为训练数据，后两组数据作为测试数据进行计算

$\begin{array}{*{20}{c}} {{X^{(0)}} = (8.15, 8.35, 8.52, 8.85, 9.13, }\\ {9.22, 9.14, 8.95)}\\ {{X^{(1)}} = (8.15, 16.5, 25.02, 33.87, 43, 52.22, }\\ {61.36, 70.31)}\\ {{Z^{(1) = }}(12.325, 20.76, 29.445, 38.435, 47.61, }\\ {56.79, 65.835)} \end{array} $

由式(5)求得传统GM(1, 1)的a≈－0.013 5，b≈8.357 7.

2) 优化辨识参数a、b时，由式(22)可得M≈0.013 5，N≈8.416 4，将M，N代入式(20)，得a≈－0.013 3，b≈8.304 2，将a、b和x⁽⁰⁾(1)代入式(9)，求得优化辨识参数GM(1, 1)算法浸润线预测值，如表 2所示.此时的预测值仍以x⁽⁰⁾(1)作为初始值时得到的，并未涉及对初始值的优化，所以对比的只是优化辨识参数前后的效果.

3) 初始值优化时，由第②步可知a≈－0.013 3，b≈8.304 2，所以A=8.80，B=8.38，C=330.65，将A、B、C代入式(28)得λ≈0.94，这时将a≈－0.013 3，b≈8.304 2，λ≈0.94代入式(30)，计算得到新的初始值，求得优化辨识参数和初始值的GM(1, 1)模型预测值，如表 2所示. 图 1所示为GM(1, 1)优化前后浸润线预测值对比. 图 2所示为误差对比.

由表 2可知，优化辨识参数的GM(1, 1)算法比传统GM(1, 1)算法有更好的预测效果，预测值的相对误差有所减小，平均相对误差有所降低，能更加逼近浸润线实测值，说明优化后求解的a、b值更准确；优化辨识参数和初始值GM(1, 1)模型的平均相对误差值最小，并且23、24日的预测值平均相对误差较前两者明显减小，说明优化初始值和辨识参数的GM(1, 1)模型有较好的预测能力.

由图 1可知，优化辨识参数和初始值GM(1, 1)模型预测值最接近原始值，说明优化方法得当，对传统算法性能有所改善.由图 2可以看出，经优化辨识参数和初始值的模型预测误差相比其他两者更小，更接近0值，在0值上下波动幅度小，说明辨识参数的求解以及初始值的设置更加合理，从而得到了较好的预测效果.

2 构建改进灰色神经网络尾矿库浸润线预测模型

传统灰色神经网络模型，即GM(1, 1)模型浸润线预测时只使用了浸润线历史样本数据，采用了单变量主因素的主导思想预测因变量变化，导致样本数据种类较局限，不能准确反映相关因素所带来的变化^[8-11].为了考虑不同因素对浸润线高度变化造成的影响，克服GM(1, 1)模型只依靠单一信息样本，即只采用历史浸润线数据变化来预测未来浸润线高度变化的时间序列模型所带来的预测精度低的问题，提出加入GM(0, N)模型. 图 3所示为改进的灰色神经网络组合模型结构图.

3 基于改进灰色神经网络的尾矿库浸润线预测模型实验验证及应用

为了分析构建的浸润线预测模型是否合理，采用2015年6月山东某金矿尾矿库主坝2号剖面3号监测点实际浸润线历史数据及对应时间的其他监测量，如表 3所示.采用传统GM(1, 1)模型、传统BP神经网络模型、传统灰色神经网络模型以及尾矿库浸润线预测模型4种模型预测浸润线高度.

表 4所示为4种模型对6月26~30日浸润线预测值.从表中可知，尾矿库浸润线预测模型的浸润线预测值与实际值较接近，其平均相对误差为1.55%，与各传统模型的相比浸润线预测值平均相对误差降低明显，满足了预期平均相对误差保持在3%以内的目标；浸润线预测值相对误差和传统模型的相比也有所降低，保证了浸润线预测值相对误差在5%以内的预期目标.

选择后验差对尾矿库浸润线模型精度等级进行检验，经计算，其均方差比值为0.324小于一级精度检验标准的0.35，在精度等级上达到了一级，满足精度要求.综上所述，所建立的尾矿库浸润线预测模型短期预测误差较小，精度等级较高，能较好预测短期内浸润线高度趋势.

4 结束语

建立了改进灰色神经网络的尾矿库浸润线预测模型.该模型在尾矿库浸润线高度短期预测时精度较高，在短期预测中，尾矿库浸润线预测模型的浸润线预测值与实际值较接近，其平均相对误差为1.55%.与各传统模型相比，采用笔者提出的预测方法进行浸润线预测时，浸润线预测值平均相对误差降低明显，并且浸润线预测值平均相对误差能保持在3%以内，浸润线预测模型预测值相对误差能保持在5%以内，预测精度达到了一级标准，取得了较好的效果.