小波消噪时间序列模型预测八五三农场年降水量

从红兴隆分局气象台收集到八五三农场1959～2005年的实测年降水序列资料。现采用小波消噪方法及前述的时间序列分析方法，以八五三农场为例，采用1959～1998年数据建立基于小波消噪的年降水时序模型，对年降水量进行模拟、预测，揭示该区降水变化规律，1999～2005年数据留作预留检验。

5.1.2.1　实测年降水序列消噪

由于年降水序列时间尺度为年，噪声不是很强，因此，在消噪时应避免掩盖真实信号。

将八五三农场实测年降水序列Pt进行1次Mallat算法分解，小波函数采用Symlets8正交小波，采用heursure启发式方法估算小波分解系数的阈值，利用软阈值方法进行消噪，再利用Mallat算法进行重构，即可得到消噪后年降水序列，如图5.1所示。

图5.1　八五三农场年降水序列启发式软阈值方式消噪（1959～1998）

（a）原始序列；（b）消噪序列；（c）噪声序列

5.1.2.2　年降水消噪序列的平稳化

从图5.1（b）可以看出，年降水消噪序列为一非平稳序列。要使转换为平稳时间序列，可以对其进行差分处理。由于在均值水平上不平稳，因此认为只需对其进行一次差分即可达到平稳，即：

式中　B——后移算子。

通过差分，使得年降水消噪序列由非平稳序列转化为平稳序列Xt，如图5.2所示。

图5.2　八三五农场差分后年降水消噪序列变化曲线

5.1.2.3　差分后年降水消噪序列的标准化处理

采用公式=（Xt-）/σX（=-5.3150，σX=143.9687）对差分后实测年降水消噪序列Xt进行标准化处理，如图5.3所示。

图5.3　八五三农场差分后年降水消噪标准化序列变化曲线

5.1.2.4　建立差分后年降水消噪标准化序列时间序列模型

1.建立趋势项模型

采用Kendall秩次相关法（丁晶等，1988；张新民等，2003）对八五三农场差分后年降水消噪标准化序列是否存在趋势成分进行检验。

构造统计量U：

式中　n——序列长度；

M——所有对偶值（，i＜j）（i=1，2，…，n-1；j=i+1，i+2，…，n）中＜出现的次数。

当n增加时，U很快收敛于标准正态分布，因此可用其进行趋势成分的显著性检验。

假定原序列无趋势。给定显著水平α后，查标准正态分布表得到Uα/2。当|U|＜Uα/2时，接受原假设，即原序列趋势不显著；反之，原序列趋势显著。

取显著水平α=0.05，查标准正态分布表得Uα/2=1.96。经过计算，对于八五三农场差分后年降水消噪标准化序列，M=394，U=0.569＜Uα/2=1.96，所以八五三农场差分后年降水消噪标准化序列无明显趋势，即趋势项Tt=0。

2.建立周期项模型

样本数n=39，则谐波数L=19。经过计算，得到序列的总方差=1.0000，均值=8.5402×10-18≈0。在显著性水平α=0.05时，Fα=3.27（n1=2，n2=36）。通过编程计算，达到显著性水平的谐波序号为12，计算结果见表5.1。所以，周期项模型为：

式中　——周期项Ct的估计值，mm。

表5.1　显著性检验计算结果

根据周期项模型绘制周期项曲线，如图5.4所示。

图5.4　周期项变化曲线

3.建立随机项模型

从差分后年降水消噪标准化序列中扣除趋势项和周期项后，就可以得到随机序列，即：

（1）正态性检验：(www.xing528.com)

经过计算，随机序列Rt的均值为=-2.2204×10-16≈0，方差为=0.7659，偏态系数CsR=-0.3827，所以随机序列Rt近似于正态分布，不必进行正态性转化。

（2）自相关分析：

样本数n=39，取最大时移m=20。通过编程计算随机序列Rt的自相关系数，并绘制自相关图，在自相关图上加绘置信水平为95%（显著性水平α=0.05）的容许限，如图5.5所示。由图5.5可以看出，随机序列Rt的自相关图具有截尾性，即在k=1时rk出现最大值，k＞1时，rk落在95%容许限范围内，因此，可以初步判定随机序列Rt不属于AR（p）模型，具体采用什么模型还应进一步进行偏相关分析。

图5.5　随机序列自相关图

图5.6　随机序列偏相关图

（3）偏相关分析：

采用尤尔——沃尔克估计法，编程计算随机序列Rt的偏相关系数φk，并绘制偏相关图，在偏相关图上加绘置信水平为95%的容许限，如图5.6所示。由图5.6可以看出，偏相关图具有明显的拖尾性，即在k=1时φk出现最大值，k＞1时，φk随阶数的增加而逐渐减小，呈拖尾状，逐渐衰减地趋向于0。再根据自相关图具有截尾性特征，可以初步判定模型形式为滑动平均模型MA（q）。

（4）模型阶数判定：

取M=≈6，采用计算机编程，经过计算可知，q=1时，P，即rk在q=1处截尾，模型阶数为1，属于MA（1）模型（丁晶等，1997；付强，2005）。

（5）模型参数估计：

对于MA（1）模型，。经过计算，求得θ1=0.4098，=0.6558。因此，所需构建的滑动平均模型为：

式中　——随机项Rt的估计值，mm。

（6）模型的进一步识别：

采用AIC准则来进一步识别MA（q）模型的阶数是否合适。通过计算，当q=1时，AIC（0，1）=39ln0.6558+2×1=-14.4541，AIC达到最小值。这说明，初步确定的模型阶数为1阶是合适的。

（7）模型检验：

由式（5.5）经过整理后，得到残差序列表达式为：

对于MA（q）模型，残差序列εt应采用递推法来确定，即假定ε1959为0，然后利用式（5.6）递推求得ε1960、ε1961、…、ε1998。利用Matlab7.0编程计算，求得残差序列εt的均值=0.0044≈0，方差=0.6179，并绘制残差序列变化曲线，如图5.7所示。

图5.7　残差序列变化曲线

计算残差序列的自相关系数rk（ε），采用自相关系数综合检验法，检验残差项εt是否为独立序列。n=39，取m=10，经过计算，统计量3.7479，根据自由度（m-q=10-1=9）查χ2表得=16.919，Q＜，所以εt为独立随机序列。

5.1.2.5　模型组合

将趋势项、周期项和随机项模型叠加再还原，就可得到八五三农场基于小波消噪的降水时序预测模型，即：

5.1.2.6　模型拟合

采用已经建立的基于小波消噪的降水时序预测模型对八五三农场1960～1998年的年降水量进行拟合，如图5.8所示。

图5.8　基于小波消噪的八五三农场年降水时序模型拟合曲线

经过计算，拟合效果评价指标后验差比值C=0.0080，小误差频率p=1.0000，最大相对误差emax=0.53%，平均相对误差=0.29%，拟合准确率为100%（相对误差e＜20%为合格），其中，C、p均达到1级标准（陈南祥，1999）。从图5.8中可以看出基于小波消噪的八五三农场年降水时序模型拟合效果良好。

5.1.2.7　模型预测

利用式（5.7）得到1999年的年降水标准化预测值，再将还原，即可得到八五三农场1999年的年降水量预测值。将加入原始序列（t=1959～1999），按照前述方法可得到八五三农场2000年的年降水量预测值，其余以此类推。然后采用未参加建模的1999～2005年的实测年降水资料进行后验预测检验，见表5.2和图5.9。

经过计算，预测检验合格率为85.71%（相对误差e＜20%为合格），达到1级标准。因此，所建的基于小波消噪的八五三农场年降水时序模型可靠性和预测精度较高，可用于预测八五三农场未来年降水量。

表5.2　基于小波消噪的八五三农场年降水时序模型后验预测误差（1999～2005）

图5.9　基于小波消噪的八五三农场年降水时序模型预测曲线（1999～2005）