传统GM(1, 1)模型固有的系统误差会给预测工作造成一定的负面影响，针对传统GM(1, 1)模型的背景值取紧邻均值的不足，周世健等^[1]、李大军等^[2]提出了加权的PGM(1, 1)模型；牛东晓等^[3]针对不同变化趋势的背景值取相同的参数值不能充分降低模型的预测误差，引入不同的参数p_i(i=1, 2, …, n-1)来改进GM(1, 1)模型背景值序列的计算公式，在不同的时刻对应不同的p_i(i=1, 2, …, n-1)；王忠桃^[4]、王宝强等^[5]将这种背景值构造方法和灰元N引入GM(1, 1)的预测模型建立了灰色非线性模型，弱化了累加序列须成某一指数变化规律，提高了拟合精度。但以上研究都只顾及了背景值而忽略初始值对预测模型精度的影响。本文通过分析GM(1, 1)和PGM(1, 1)的不足，在对初始值进行改正的同时引入王忠桃、王宝强等人提出的背景值构造方法和灰元N建立综合优化的灰色非线性模型，顾及引入的未知参数较多，用传统方法难以求解，而粒子群算法适用多目标优化且具有收敛速度快、结构简单的特点，故采用粒子群算法求解最优参数，从而建立基于粒子群算法和加权灰色组合的PSO-GM模型。通过两个实例对比分析，验证了优化模型的可靠性和实用性。

1 GM(1, 1)缺陷分析和PGM(1, 1)的建模机理 1.1 GM(1, 1)缺陷分析

由灰色模型的建模过程可知，传统GM(1, 1)模型背景值构造公式如下：

(1)

其中x⁽¹⁾(k)=，x⁽⁰⁾(i)为原始序列。文献[1-3]认为，当原始序列数据变化平缓时，构造出的背景值是合适的，模型的偏差也较小。然而，事实上即便在很短的时间Δt里仍不可避免地出现突变量或序列的变化不是以一种平均速率发生的，即当原始序列波动较大时，应用上式构造出来的背景值会出现较大的滞后误差。为此从预测的精度出发引入了一种新的背景值取值方法：

(2)

由公式推导可得p和a存在如下关系：

(3)

可证明，当|a|趋向0时，p趋向于0.5；当|a|较小时，p非常接近0.5；当|a|较大时，p与0.5偏离较大。因此，简单地取p等于0.5是导致传统GM(1, 1)模型|a|较大时预测失效的原因。

1.2 PGM(1,1)模型的建立

针对上述问题，利用式(2)的背景值取值方法建立PGM(1, 1)模型。具体算法为：p的确定是先取p=0.5作为近似值，以原始序列和模拟序列的平均相对误差最小为约束条件，按传统GM(1,1)模型建立方法，得到灰参数=[a⁰, u⁰]^T及模拟近似值序列L⁰=，将白化方程写成离散形式并进行泰勒级数展开，建立误差方程：

(4)

组成矩阵形式为：

(5)

其中，V=[v(2), v(3), …, v(n)]^T; δ_X=[δ_p, δ_a, δ_u]; l=[l₍₂₎, l₍₃₎, …, l_(n)]; l_(k+1)=L_(k+1)-, k=1, 2, …, n－1。利用迭代法求解灰参数和最佳背景值p^[7]。

2 基于粒子群算法的PSO-GM(1, 1, N, p, ξ)

文献[3]针对不同变化趋势的背景值取相同的参数值不能充分降低模型的预测误差，引入不同的参数p_i(i=1, 2, …, n－1)来改进GM(1, 1)模型背景值序列的计算公式，在不同时刻对应不同p_i(i=1, 2, …, n－1)。文献[4]将这种背景值构造方法和灰元N同时引入GM(1, 1)的预测模型，建立了灰色非线性模型。本文在这些基础上通过引入初始值改正和上面提到的背景值构造方法及灰元N，建立优化的灰色非线性模型，命名为GM(1, 1, N, p, ξ)模型。其中ξ(任意值)为x⁽⁰⁾(1)=x⁽⁰⁾(1)+ξ的修正值，并以此作为四阶龙格-库塔法求解累加序列的初始值。

2.1 PSO-GM(1, 1, N, p, ξ)建模机理

设有非负序列X⁰=(x⁰(1), x⁰(2), …, x⁰(n))，对该序列进行一次累加X⁽¹⁾=(x⁽¹⁾(1), x⁽¹⁾(2), …, x⁽¹⁾(n))，其中x⁽¹⁾(k)=，k=1, 2, …, n，则X⁽¹⁾可以用曲线f(t, a, u, N)拟合，N∈(0, 2)的待定常数^[5]，并且f(t, a, u, N)满足以下微分方程：

(6)

上式也可以表示为：

(7)

灰参数=[a, u]^T=(B^TB)^－1B^TY，其中

求出灰参数后，采用四阶龙格-库塔法以x⁽⁰⁾(1)=x⁽⁰⁾(1)+ξ(ξ为任意值)为初始值得到累加序列的模拟值和预测值，通过累减得到原始序列的模拟值和预测值^[11]。

由于引入的参数较多，一般的方法难以求解，而粒子群算法对求解非线性、不可微的复杂优化问题收敛速度快、解质量高，并且粒子群算法程序实现简单，适合多目标优化问题求解。具体算法可描述为：设搜索区域为D维空间，种群中共有M个粒子，粒子i在第t次迭代的位置为=，粒子的速度为v_i=。在每次迭代中，粒子i跟踪个体极值点P_i=[p_i1, p_i2, …, p_iD]^T和全局极值点P_g=[p_g1, p_g2, …, p_gD]^T来更新自己的位置和速度，准则为：

(8)

为使得粒子的速度在更新时不至过大，设置v∈(－2, 2)。公式的第一部分为粒子先前的速度；第二部分表示粒子自身的经验和记忆；第三部分表示粒子间信息共享与相互合作。式中，t为当前迭代次数；c₁为粒子更新自身最优值的学习因子，c₂为粒子更新全部粒子最优值的学习因子，c₁、c₂一般取值为2；r₁和r₂为[0, 1]区间的随机数；ω为惯性权重，它是保持原来速度大小的因子。ω较大时有利于跳出局部最优解，进行全局搜索；ω较小时有利于局部寻优，加快算法收敛。文中采用随迭代次数线性递减的ω值：

(9)

式中，ω_max=0.95, ω_min=0.4，i为当前迭代次数，i_max为最大迭代次数。

适应度函数为：

(10)

具体算法流程:

1) 初始化粒子群，包括种群规模m=30，粒子维数为n，最大迭代次数i_max=1 000，粒子位置y_i=[λ_i1, λ_i2…λ_in－1, N_i, ξ_i]和速度V_i=[v_{i, 1}, v_{i, 2}, …, v_{i, n}, v_{i, n+1}]，其中0.1≤λ_i≤1, N∈(0, 2)。根据适应度函数计算每个粒子的适应度^[3]。

2) 根据全部粒子适应度大小初始化个体极值P_i和计算全局极值P_g。

3) 迭代寻优，由式(9)更新惯性权重，由式(8)更新每个粒子的速度和位置，并限制其速度满足v∈(－2, 2)，计算每个粒子的适应度f_i，并和个体极值比较。若f_i < f(P_i)，则将f_i替换为当前粒子的最佳值，f(P_i)=f_i，同时更新当前粒子的最优位置P_i=y_i。

4) 将每个粒子的当前最优值和全局最优值比较，若f(P_i) < f(P_g)，则将f(P_i)替换为当前全局最优值，即f(P_g)=f(P_i)，同时更新全局最优位置。

5) 判断迭代次数是否达到最大迭代值，若是则输出全局极值和全局最优位置，以及计算出对应的模拟和预测序列。否则转入第2步，直至满足迭代条件。

2.2 优化模型的精度检验

优化模型的精度评定主要从3个方面着手：1)计算均方差比(S₁为原始序列的均方差，S₂为残差序列的均方差)。对于评定预测模型的好坏，C值越小越好，因为C小说明S₂小、S₁大，即模拟数据方差小；原始数据方差大，说明模拟数据比较集中，摆动幅度小；原始数据比较分散，摆动幅度大，为使模拟效果好，要求C值尽可能小些，一般要求C小于0.35，最大不超过0.65^[8]。2)计算相对误差均值MAPE=(∑(e/x⁽⁰⁾))/n。e为残差值，n为数据个数，相对误差均值越小越好，它也是评价模型精度的常用指标。3)计算灰色关联度。ρ为分辨系数，0 < ρ < 1，一般取ρ=0.5^[9]，拟合序列和原始序列的灰色关联度越大，说明拟合序列越好地保持了原始序列的内在变化规律，模型的拟合和预测效果越好。通过以上3个精度评价指标综合评价优化模型的预测效果。

3 算例实验与分析

以某大坝边坡变形监测线上的一观测点设为A点的位移值(数据取自文献[10])，取前7期分别用GM(1, 1)、PGM(1, 1)、PSO-GM(1, 1, N, p, ξ)建立模型，PSD-GM得到的最优参数y=[0.363, 0.100, 1.000, 0.168, 0.100, 0.823, 0.997, 0.210]，具体数据对比分析见表 1、表 2。

表 1、表 2中各模型拟合和预测结果表明，PGM(1, 1)模型和GM(1, 1)模型模拟与预测精度一致，说明当序列呈一定单调性变化时，使得PGM(1, 1)最优的p刚好为0.5，PGM(1, 1)便还原为GM(1, 1)，而PSD-GM通过粒子群算法优化参数，使得模型的模拟和预测精度在3种模型中最高。PSD-GM得到的灰色关联度R=0.8，而GM(1, 1)和PGM(1, 1)的灰色关联度R=0.629，R越大说明拟合值序列越好地保持了原始序列的内在变化规律，模型的模拟和预测效果越好。对比3种模型的精度评定结果表明，本文模型的均方差比仅为0.003，相对误差均值为0.000 67，而GM(1, 1)、PGM(1, 1)的均方差比和相对误差均值分别为0.015和0.008，相比减少了一个数量级。究其原因，是本文模型对初始值和背景值都进行了改进，相对传统模型和PGM(1, 1)模型精度都得到较大提高。

为了进一步说明模型的可靠性和适用性，引用文献[1]中实例，即湖北秭归县链子崖危岩体某监测点1978~1993年的监测资料。首先以1978~1987年实测值为原始数据序列，分别建立10维的GM(1, 1)、PGM(1, 1)和PSO-GM(1, 1, N, p, ξ)，本文模型寻得的最优参数y=[0.100, 0.100, 0.100, 0.100, 0.100, 0.100, 1.000, 0.100, 0.100, 0.580, 2.520]。同时用这3种模型来预测1988~1993年变形值并与实测值比较，具体数据见表 3、表 4、图 1。

表 3、表 4、图 1的结果显示，GM(1, 1)模型得到的均方差比为0.278，相对误差均值为0.136，灰色关联度为0.532；PGM(1, 1)模型得到的均方差比为0.289，相对误差均值为0.125，灰色关联度为0.6；PSD-GM的均方差比为0.244，相对误差均值为0.086，灰色关联度为0.673。说明PGM(1, 1)模型的预测精度总体略优于传统GM(1, 1)模型，验证了优化背景值能提高模型的模拟和预测精度。而PSD-GM模型得到的3种精度指标均优于PGM(1, 1)和GM(1, 1)模型，说明针对初始值和背景值同时改进的PSD-GM模型，既能更好地保存原始序列的内在变化规律，又进一步提高了模型的模拟和预测精度。从3种模型的预测结果对比分析看，随着预测时间的增长，3种模型的预测误差都大致呈递增趋势。通过3种模型的残差对比分析，随着预测区间的增大，除了前两期预测误差略大于PGM(1, 1)和GM(1, 1)，其余残差均优于PGM(1, 1)和GM(1, 1)模型，个别甚至提高了一个数量级。两次实验结果表明，基于初始值和背景值改进的PSO-GM(1, 1, N, p, ξ)模型兼具有GM(1, 1)和PGM(1, 1)的特性，同时也克服了两种模型存在的缺陷，明显提高了前两种模型的预测精度。

4 结语

本文基于初始值和背景值同时改进的灰色模型，通过引入不同背景权重p和灰指数N，并用改进的初始值由龙格-库塔法求解累加序列的模拟值和预测值，通过粒子群算法求解最优参数。优化模型克服了传统GM(1, 1)和PGM(1, 1)的不足，同时又兼具两种模型的特性，即当N=1, p₁=p₂=…=p_n, ξ=0时即为PGM(1, 1)模型；当N=1，p_i=0.5，ξ=0时即为GM(1, 1)模型。通过两个工程实例证明，优化模型的拟合精度高，预测效果好，预测周期长，能够很好地反映形变发生的内在规律，相比其他两种模型有更高的实用价值。