塑料具有比强度高、弹性高、耐磨与减摩性好、原料来源丰富、易加工、品种繁多、成本相对较低等优点，不但可以用作结构材料，而且作为功能材料亦有着广泛的发展前景^[1]。但塑料易老化，特别是在苛刻的工况下，塑料的老化常导致设备过早失效，材料大量流失，不但造成经济损失重大，资源浪费严重，甚至会因材料的失效造成环境污染。老化已成为限制塑料进一步发展和应用的关键问题之一^[2]。

由于老化原因的多样性与老化机理的复杂性，通过老化原因定量表述其引起的力学性能变化比较困难，目前主要通过老化试验来确定塑料的力学性能变化，进而得出其使用寿命。但自然老化试验耗时很长^[3]，可能数年，甚至数十年，不仅很难跟上材料研究的高速发展，而且费用相当昂贵，获得的试验数据也很有限，而人工加速老化试验与实际工况又存在不可避免的差异。针对这一状况，近年来，在已知相关试验数据的基础上，建立数学模型，对材料的老化行为进行预测已成为研究人员关注的热点之一^[4-7]。

灰色GM(1，1)模型是基于灰色系统理论的一种预测方法^[8]，具有要求历史数据少、运算方便、易于检验等优点^[9]，已在诸多领域得到了广泛应用^[10-14]。但实践中发现，数据的离散程度较大时，灰色GM(1，1)模型的预测精度较差^[15]。为解决这一问题，提出了灰色残差预测^[16-17]，但该方法在残差序列有正有负的情况下，无法确定预测值中残差修正值的正负。针对这一情况，笔者尝试引入马尔克夫过程，以其确定未来时刻(即预测值)残差修正值的正负，建立基于马尔克夫的灰色残差GM(1，1)模型，并以大气自然老化对LDPE棚模拉伸强度的影响为例，研究所建模型在塑料老化领域中的适用性。

1 模型的建立 1.1 灰色GM(1，1)模型

设原始时间序列X⁽⁰⁾有n个时刻的非负实际值，则X⁽⁰⁾ ={x⁽⁰⁾(1)，x⁽⁰⁾(2)，…，x⁽⁰⁾(n)}，由下式生成一次累加序列X⁽¹⁾，X⁽¹⁾ ={x⁽¹⁾(1)，x⁽¹⁾(2)，…，x⁽¹⁾(n)}。

$ {x^{\left( 1 \right)}}\left( k \right) = \sum\limits_{i = 1}^k {{x^{\left( 0 \right)}}\left( i \right)} ,k = 1,2 \cdots ,n; $

(1)

由下式生成X⁽¹⁾的紧邻均值序列Z⁽¹⁾，Z⁽¹⁾ ={z⁽¹⁾(2)，z⁽¹⁾(3)，…，z⁽¹⁾(n)}。

$ {z^{\left( 1 \right)}}\left( k \right) = 1/2\left( {{x^{\left( 1 \right)}}\left( k \right) + {x^{\left( 1 \right)}}\left( {k - 1} \right)} \right),k = 1,2, \cdots ,n; $

(2)

构建微分方程dX⁽¹⁾/dt +aX⁽¹⁾ =b，其中，参数a、b可由下式得到。

$ {\left[ {a,b} \right]^{\rm{T}}} = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Y}}, $

(3)

式中：$\boldsymbol{Y}=\left[\begin{matrix} {{x}^{(0)}}(2) \\ {{x}^{(0)}}(3) \\ \vdots \\ {{x}^{(0)}}(n) \\ \end{matrix} \right], \boldsymbol{B}=\left[\begin{matrix}-{{z}^{(1)}}(2) & 1 \\-{{z}^{(1)}}(3) & 1 \\ \vdots & \vdots \\-{{z}^{(1)}}(n) & 1 \\ \end{matrix} \right]$。

求解上述微分方程，得到如下GM(1，1)模型：

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

(4)

其还原值应满足如下公式：

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

(5)

1.2 灰色残差GM(1，1)模型

由下式建立残差(即原始值与模拟值之差)序列ε₁(0)，ε₁(0)={ε₁(0)(2)，…，ε₁(0)(n)}。

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

(6)

采用下式对ε₁(0)作非负处理，得到ε⁽⁰⁾，ε⁽⁰⁾= {ε⁽⁰⁾(2)，…，ε⁽⁰⁾(n)}。

$ {\varepsilon ^{\left( 0 \right)}}\left( k \right) = \left| {\varepsilon _1^{\left( 0 \right)}\left( k \right)} \right|,k = 2,3, \cdots ,n。$

(7)

对ε⁽⁰⁾进行GM(1，1)建模，具体方法同上1.1，得到如下预测模型。

$ {{\hat \varepsilon }^{\left( 1 \right)}}\left( {k + 1} \right) = \left( {{\varepsilon {\left( 0 \right)\left( 2 \right)}} - \frac{{{b_\varepsilon }}}{{{a_\varepsilon }}}} \right){{\rm{e}}^{ - {a_\varepsilon }k}} + \frac{{{b_\varepsilon }}}{{{a_\varepsilon }}},k = 2,3, \cdots ,n; $

(8)

其中，参数a_ε、b_ε的求法与公式(3)类似。

通过下式进行累减还原，得到残差修正值${{\hat{\varepsilon }}^{(0)}}$。

$ \begin{array}{*{20}{c}} {{{\hat \varepsilon }^{\left( 0 \right)}}\left( {k + 1} \right) = {{\hat \varepsilon }^{\left( 1 \right)}}\left( {k + 1} \right) - {{\hat \varepsilon }^{\left( 1 \right)}}\left( k \right) = \left( {1 - {{\rm{e}}^{{a_\varepsilon }}}} \right),}\\ {\left( {{\varepsilon ^{\left( 0 \right)}} \left( 2 \right)- \frac{{{b_\varepsilon }}}{{{a_\varepsilon }}}} \right){{\rm{e}}^{ - {a_\varepsilon }k}},k = 2,3, \cdots ,n。} \end{array} $

(9)

用${{\hat{\varepsilon }}^{(0)}}$修正${{\hat{X }}^{(0)}}$，具体公式如下。

$ \begin{array}{*{20}{c}} {\hat x_\varepsilon ^{\left( 0 \right)}\left( {k + 1} \right) = {{\hat x}^{\left( 0 \right)}}\left( {k + 1} \right) \pm {{\hat \varepsilon }^{\left( 0 \right)}}\left( {k + 1} \right) = \left( {1 - {{\rm{e}}^a}} \right)\left( {{x^{\left( 0 \right)}}\left( 1 \right) - \frac{b}{a}} \right){{\rm{e}}^{ - ak}} \pm }\\ {\left( {1 - {{\rm{e}}^{{a_\varepsilon }}}} \right)\left( {{\varepsilon {\left( 0 \right)\left( 2 \right)}} - \frac{{{b_\varepsilon }}}{{{a_\varepsilon }}}} \right){{\rm{e}}^{ - {a_\varepsilon }k}},k = 2,3, \cdots ,n。} \end{array} $

(10)

上式称为灰色残差修正GM(1，1)模型，简称灰色残差GM(1，1)模型，其中每个残差修正值的符号应与ε₁(0)中相应数值的符号保持一致，并用下面的马尔克夫过程来确定预测值的残差修正值符号。

1.3 马尔克夫过程预测模型

马尔可夫过程基于“系统每一时刻的状态仅仅取决于前一时刻的状态，而与其过去的历史无关”这一概念，研究系统的状态及状态的转移，进而根据系统状态之间的转移概率(即从某一状态转换到另一状态的可能性)来预测系统状态未来的发展趋势^[18]。

采用马尔可夫过程确定未来某一时刻残差修正值正、负号的具体步骤如下：

1) 确定状态。定义两个状态，状态1表示残差为正值，状态2表示残差为负值。

2) 根据残差序列的正负状态，求出状态转移概率矩阵P。

$ \mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} {{p_{11}}}&{{p_{12}}}\\ {{p_{21}}}&{{p_{22}}} \end{array}} \right], $

其中：元素p_ij表示从状态i转移到状态j的概率，p_ij≥0，且$\sum\limits_{j=1}^{2}{{{p}_{ij}}}=1 $。p_ij可由下式确定：

$ {p_{ij}} = {s_{ij}}/{s_i},i = 1,2;j = 1,2。$

(11)

其中：s_ij为状态i转移到状态j的次数；s_i为状态i出现的次数。

3) 确定初始(即当前)状态向量。设s⁽⁰⁾=(s₁⁽⁰⁾，s₂⁽⁰⁾)为初始状态向量，s₁⁽⁰⁾表示处于状态1的概率，s₂⁽⁰⁾表示处于状态2的概率。

4) 根据如下所示的状态转移公式，求出下一时刻即第t时刻状态转移的结果s^(t)。

$ {s^{\left( t \right)}} = {\mathit{\boldsymbol{s}}^{\left( 0 \right)}} \cdot \mathit{\boldsymbol{P}}。$

(12)

以s^(t)中概率大的状态作为t时刻残差修正值的状态(即正负号)，若两状态的概率相等，则取上一时刻残差修正值的状态作为t时刻残差修正值的状态。

1.4 误差检验

由下式计算相对误差。

$ {\Delta _k}\varepsilon = \left| {{x^{\left( 0 \right)}}\left( k \right) - {{\hat x}^{\left( 0 \right)}}\left( k \right)} \right|/{x^{\left( 0 \right)}}\left( k \right),k = 1,2, \cdots ,n。$

(13)

2 应用实例

LDPE棚模在光、热、氧等大气自然因素的综合作用下，拉伸强度随老化时间的延长而下降，具体数据如表 1^[19]所示。下面以表 1中老化0到15个月对应的实际拉伸强度为原始数据，根据上述所建模型，预测老化18个月和21个月时LDPE棚模的拉伸强度，并与表 1中老化18个月和21个月对应的实际拉伸强度进行比较，研究所建模型在塑料老化行为预测中的适用性。

由表 1中LDPE棚模0到15个月的拉伸强度，可得原始数据数列X⁽⁰⁾={18.8，18.2，17.1，14.8，13.9，13.3}；由式(1)、式(2)分别得到：X⁽¹⁾={18.800 0，37.000 0，54.100 0，68.900 0，82.800 0，96.100 0}，Z⁽¹⁾={27.900 0，45.550 0，61.500 0，75.850 0，89.450 0}；由式(3)可得a=0.085 2，b=20.574 0，并由式(4)可得预测模型为${{\hat{x}}^{(1)}}$(k+1)=-222.678 9e^{-0.085 2k}+241.478 8，即${{\hat{x}}^{(1)}}$= {36.986 5，53.687 7，69.024 8，83.109 4，96.043 7}；根据式(5)做累减还原，可得灰色GM(1，1)模拟值${{\hat{x}}^{(1)}}$(见表 2)；由式(6)可得残差序列ε₁(0)(见表 2)；通过式(7)将ε₁(0)取绝对值，带入式(8)得到残差模型${{\hat{\varepsilon }}^{(1)}}$(k+1)=-3.825 2e^{-0.122 8k}+3.838 8，即${{\hat{\varepsilon }}^{(1)}}$={0.846 5，1.192 3，1.498 1，1.768 6}；根据式(9)做累减还原，可得ε⁽⁰⁾的灰色残差GM(1，1)模拟值${{\hat{\varepsilon }}^{(1)}}$(见表 2)；由式(10)得到修正后的模拟值(见表 2)；由式(13)得到相对误差(见表 2)。

由表 2可知，基于马尔克夫的灰色残差GM(1，1)模型的模拟结果中，拉伸强度的最大相对误差为1.29%，平均误差仅为0.73%，模拟精度大大高于文献[19]中采用灰色GM(1，1)模型得到的结果(最大相对误差为3.62%，平均误差为2.02%)。大气自然老化条件下，LDPE棚模拉伸强度随时间变化的离散程度较大，因此灰色GM(1，1)模型得到的模拟结果相对误差较大，而灰色残差GM(1，1)模型采用残差序列对原点附近的误差进行了修正，因此模拟结果更符合实际情况。

应用上述模型对老化18个月时的拉伸强度进行预测，其中残差修正值的符号确定过程如下。

由表 2中的残差序列可知，残差由正向正转移的次数是2，正值出现的次数为3，因此，由式(11)得到正向正转移的概率为p₁₁=2/3，正向负转移的概率为p₁₂=1/3；同理，负向正的转移次数为1，负值出现的次数为2。因此，负向负转移的概率为p₂₁=1/2，负向负转移的概率为p₂₂=1/2。由此得到马尔克夫状态转移概率矩阵P为

$ \mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} {{p_{11}}}&{{p_{12}}}\\ {{p_{21}}}&{{p_{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {2/3}&{1/3}\\ {1/2}&{1/2} \end{array}} \right]。$

将老化15个月的拉伸强度作为初始状态，由表 2可知其残差值为正，因而初始向量s⁽⁰⁾=(1，0)。

由式(12)得到第18个月状态转移的结果为

$ {s^{\left( {18} \right)}} = {\mathit{\boldsymbol{s}}^{\left( 0 \right)}} \cdot \mathit{\boldsymbol{P}} = \left( {1,0} \right) \cdot \left[ {\begin{array}{*{20}{c}} {2/3}&{1/3}\\ {1/2}&{1/2} \end{array}} \right] = \left( {2/3,1/3} \right)。$

即第18个月残差修正值为正的概率为2/3，为负的概率为1/3，因此取正值。由基于马尔克夫的灰色残差GM(1，1)模型得到第18个月的拉伸强度预测结果(见表 3)。表 3还列出了第21个月的预测结果，具体过程略。

由表 3可知，由基于马尔克夫的灰色残差GM(1，1)模型得到的第18个月和21个月的拉伸强度预测值与实际值的相对误差分别为1.49%和4.96%，预测精度比文献[19]中采用灰色GM(1，1)模型得到的结果(相对误差分别为3.40%和6.75%)有明显提高。上述结果表明：在根据灰色残差GM(1，1)模型得到未来时刻预测值的残差修正值的基础上，采用马尔克夫过程确定残差修正值的正负取值，得到的预测结果更符合实际情况。适用于塑料老化行为的预测。

3 结论

1) 在采用马尔克夫过程确定预测值残差修正值的正负号的基础上，以残差修正值对灰色GM(1，1)模型的模拟值与预测值进行修正，构建了基于马尔克夫的灰色残差GM(1，1)模型。

2) 以所建模型对大气自然老化环境下LDPE棚模的拉伸强度的变化进行预测。结果表明：基于马尔克夫的灰色残差GM(1，1)模型的模拟精度与预测精度明显优于灰色GM(1，1)模型，可用于塑料老化行为的预测。

3) 基于马尔克夫的灰色残差GM(1，1)模型在保留灰色GM(1，1)模型所需实验数据少、运算方便、易于检验等优点的基础上，大大提高了预测精度，拓宽了灰色理论的应用范围，为塑料老化行为的预测提供了一种简易而可靠的新途径。