基于多元周期平稳时序分析及灰色理论的水华预测及因素分析方法

摘要

本发明公开了一种基于多元周期平稳时序分析及灰色理论的水华预测及因素分析方法，属于环境工程技术领域。所述方法包括进行监测数据采集及预处理；特征因素筛选；特征因素时序建模；特征因素时序预测和特征因素分析的步骤。本发明解决了多元时序的多元周期平稳性检验问题；给出了水华特征因素的筛选方法，为水华特征因素多元时序的合理建模提供了必要的条件，也为水华特征因素的因素分析提供了基础。本发明对多特征因素的水华形成过程描述更加全面，从而提高了水华预测结果的准确度、可信度。根据本发明方法给出的水华预测结果与实测数据的误差，可以得到各影响因素与水华发生现象联系的密切程度。

说明书

技术领域

本发明涉及一种水华预测及因素分析方法，属于环境工程技术领域。

背景技术

水华是指出现在富营养化水体中，当具备适宜的光照、水温、气候及水文等有利于藻类 生长和聚集的环境条件时，藻类爆发性繁殖聚集并达到一定浓度的一种现象。大规模的水华 危害很大，它不仅会破坏水体生态系统的结构和功能，危害人类健康，而且还会降低水资源 利用效能，威胁水资源的可持续开发和利用，造成巨大的经济损失。从总体上看，目前还缺 少能在短期内有效治理水华的技术和手段。因此在水华未得到有效治理之前，对水华的发生 进行准确预测便于有关部门采取应对措施，降低危害，因此，水华防治具有重要的科学意义 和应用价值。

在水华的防治工作中，水华预测一直都是一个难点，这是因为水华形成的随机过程是一 个复杂的物理、化学和生物综合反应过程，与水华形成有关的特征因素较多，相互作用关系 密切，构成一个复杂的水生生态系统，由于生态过程内在机制的复杂性、人类活动的影响及 水质信息自动测报、接收等条件的限制，目前对水华形成机理和爆发的成因还没有完全清楚， 因而直接建立水华机理生态模型比较困难。随着水华预测研究的深入，部分学者尝试建立精 确可靠的数学模型来对水华形成的随机过程进行预测，取得了一定的进展。然而，现有的水 华预测方法仍存在预测精度不高，预测步长较短的问题。

此外，不同特征因素与水华形成的相关程度也不尽相同，因此水华形成的因素分析也是 水华防治工作中的一项重要任务。

时间序列分析是一种研究随机过程的重要数学预测工具。传统一元时间序列预测的基本 假设是可利用的信息都包含在历史数据中，历史数据决定了未来时间序列的变化，但是这种 假设过多的依赖于时间序列历史数据，而忽略了事物之间的关联性，也忽略了外部因素对时 间序列走势的影响。而多元时间序列的建模和预测则不仅考虑了时间序列历史数据对未来的 影响，更多的考虑到多元协方差和相关系数在建模过程中的重要作用，比一元时间序列建模 和预测更具备应用潜力，是一种适于描述和预测在多种特征因素作用下水华形成随机过程的 方法。因此采用多元时序分析方法，对水华形成的特征因素多元时序建模，从而进行水华预 测为一种有效途径。

然而，时序分析方法虽然可以充分挖掘数据中的随机性信息，该方法存在预测精度的保 证依赖于大样本量的局限。若采用常规的直接预测方法进行时序模型预测，则预测精度难以 保证，从而降低了水华预测结果的可信度。灰色理论是研究少数据不确定性的理论，即研究 在少数据不确定性的背景下，数据的处理、现象的分析、模型的建立、发展趋势的预测、事 物的决策以及系统的控制与状态的评估。采用灰色理论预测多元时序模型，可以在保留时序 建模优点的基础上，很好地克服时序模型预测的局限性，提高模型的预测精度。

另外，多元时序分析不仅适于建模预测，也是一种因素分析的有效工具。因此可利用多 元时序分析预测结果对水华形成过程中的各种特征因素进行因素分析。

应用多元时序分析及灰色理论进行水华预测及因素分析还需要解决以下问题：

1.由于水华的发生具有随机性，描述水华形成过程中各特征因素的随机性波动，不但需要考 虑了同一时刻多个特征因素之间的互相关性，还要考虑各特征因素在不同时刻的自相关性。

2.水华的发生与季节变化相关，因而水华的发生还具有周期性，描述水华形成过程中各特征 因素的周期性波动，不仅需要考虑周期性的环境季节变化对水华形成过程的影响，还要考虑 不同特征因素间多重周期性变化的交互影响。

3.水华形成过程中特征因素的随机性和周期性波动在不同特征因素间的交互影响也是需要考 虑的问题。

4.如何将适用于少数据预测灰色理论方法与所建的时序模型相结合，从而提高所建时序模型 的预测精度，需要具体研究。

5.如何基于所建的时序模型预测结果来进行水华特征因素分析。

发明内容

本发明的目的是为了解决现有的水华预测结果不够准确，水华形成过程中多种特征因素 间的交互影响建模困难，以及不同特征因素与水华发生的相关性程度的判定问题，采取基于 多元周期平稳时序分析以及灰色理论的技术手段，达到通过对水华形成过程中的特征因素时 序建模分析得到与实际情况更为相符的水华预测效果以及因素分析结果。

为便于说明，本说明书中所有未经解释的名词及字母含义均由下述假设解释：与水华现 象有关的特征因素分为两种：一种是影响水华发生的特征因素，例如氮、磷、pH值、溶解氧、 水温、光照度等，以下叫做影响因素；另一种是表征水华发生的特征因素，例如叶绿素浓度、 藻密度等，以下叫做表征因素。以Y_t表示t时刻的特征因素向量；以y_it表示第i个特征因素 在t时刻的量值，总采样时间为N，t＝1,2,…，N，共有n个特征因素，i＝1,2,…,n。

本发明提供的基于多元周期平稳时序分析及灰色理论的水华预测及因素分析方法主要包 括以下五个步骤：

步骤一、监测数据采集及预处理；

对可能影响或表征水华发生的多个特征因素进行监测。由监测设备采集到的多个特征因 素的监测数据随时间变化而变化，因此为多元时间序列。通常情况下，原始特征因素监测数 据包含有异常采样点、均值漂移以及部分采样点数据缺失等等现象，难以直接对其进行时序 分析。为消除这些现象对时序分析造成的影响，提高特征因素时序模型的拟合精度，并且统 一原始特征因素时序的初值，应对每个特征因素原始时序分别作预处理。

步骤二、特征因素筛选；

水华的发生具有随机性，并且与季节变化相关，因而水华的发生还具有周期性，因此其 特征因素的多元时序在自然界中呈现随机性、周期性变化，但其多个特征因素之间的变化机 理是相对稳定的，不会随时间变化，即特征因素的多元时序的方差不随时间而变化，因此特 征因素时序应为多元周期平稳时序。

为此需要对预处理后的特征因素多元时序进行多元周期平稳性检验，并根据多元周期平 稳性检验结果进行特征因素筛选，只有多元周期平稳性检验通过的特征因素才适合采用本发 明方法。本发明基于延迟交叉相关矩阵法给出多元周期平稳性检验方法，可见具体实施方式 中的步骤二。本发明的特征因素筛选方法如下：

1、每个特征因素的一元周期平稳性分析。由于一元周期平稳是多元周期平稳的必要条件， 因此首先对每个特征因素的一元周期时序进行一元周期平稳性检验。采用自相关函数检验方 法。将一元周期平稳性检验不通过的特征因素剔除；

2、对一元周期平稳性检验通过的所有表征因素，例如叶绿素浓度、藻密度等，采用本发 明给出的多元周期平稳性检验方法检验，找出能够通过多元周期平稳性检验的包含最多表征 因素的组合；

3、对通过多元周期平稳性检验的表征因素组合，加入尽可能多的通过一元周期稳定性检 验的影响因素，以组成新的特征因素组合，直至新组成的特征因素组合不能通过多元周期平 稳性检验或包含全部影响因素为止。

这样，根据特征因素筛选结果，对符合多元周期平稳条件的特征因素组合进行以下步骤 三和步骤四。

步骤三、特征因素时序建模；

1、确定特征因素时序结构；

水华特征因素的多元时序在自然界中呈现随机性、周期性变化，但不存在趋势性变化， 因此将t时刻的特征因素向量Y_t分解为周期项C_t和随机项R_t的叠加，

Y_t＝C_t+R_t    (1)

$Y_t=\left(\begin{array}{c}{y}_{1t}\\ y_{2t}\\ .\\ .\\ .\\ y_{\mathrm{nt}}\end{array}\right),$$C_t=\left(\begin{array}{c}{c}_{1t}\\ c_{2t}\\ .\\ .\\ .\\ c_{\mathrm{nt}}\end{array}\right),$$R_t=\left(\begin{array}{c}{r}_{1t}\\ r_{2t}\\ .\\ .\\ .\\ r_{\mathrm{nt}}\end{array}\right)$

其中c_it为第i个特征因素的周期项，r_it为第i个特征因素的随机项，i＝1,2,…,n，其中n为特 征因素的个数。

2、建立特征因素时序周期项模型；

由于多个特征因素之间周期性变化存在交互性影响，周期项C_t若采用一元周期模型描述 则无法反映这种交互性影响，本发明采用适用于挖掘数据潜在周期性规律以及反映多重周期 间交互性影响的多重潜周期模型描述，即

其中C(t)为多重潜周期模型的多重潜周期函数，q为潜周期角频率个数，ω_j为第j个角频率， a_ij为多重潜周期模型的第i个特征因素的第j个角频率对应的幅值，为多重潜周期模型的 第i个特征因素的第j个角频率对应的相位，i＝1,2,…,n。

3、建立特征因素时序随机项模型；

从Y_t减去C_t后，对随机项，即Y_t的平稳随机性部分R_t采用应用最广泛、建模简单且适 于预测的多元自回归模型描述，即

$R_t=\underset{j=1}{\overset{p}{\Sigma }}{H}_j{R}_{t-j}+E_t---\left(3\right)$

$E_t=\left(\begin{array}{c}{ϵ}_{1t}\\ {ϵ}_{2t}\\ .\\ .\\ .\\ {ϵ}_{\mathrm{nt}}\end{array}\right)$

其中p为多元自回归阶数，H_j为n×n多元自回归系数矩阵，R_t-j为在t-j时刻下的随机项，E_t为n维白噪声向量，η_ikj为第i个特征因素对第k个特征因素的多元自回归系数，ε_it为第i个 特征因素的白噪声，i,k=1,2,…,n。

4、建立特征因素多元周期平稳时序模型。

将周期项C_t的多重潜周期模型公式(2)与随机项R_t的多元自回归模型公式(3)结合。由于 水华形成过程机理复杂，其特征因素的随机性和周期性变化之间也可能存在交互影响，因此 随机项R_t和周期项C_t之间也可能具有相关性，若仅将二者的模型简单相加，则会忽略这种相 关性。为此，为了提高模型拟合精度，本发明将周期项多重潜周期模型代入随机项多元自回 归模型进行迭代运算，提出特征因素的多重潜周期多元自回归混合模型

$Y_t=\underset{j=1}{\overset{p}{\Sigma }}{H}_j{Y}_{t-j}+C^{*}\left(t\right)+E_t---\left(4\right)$

其中Y_t-j为在t-j时刻下的特征因素向量，$C^{*}\left(t\right)=\left(\begin{array}{c}\underset{j=1}{\overset{q}{\Sigma }}{b}_{1j}\cos ({\omega }_{j}t+{\phi }_{1j})\\ \underset{j=1}{\overset{q}{\Sigma }}{b}_{2j}\cos ({\omega }_{j}t+{\phi }_{2j})\\ .\\ .\\ .\\ \underset{j=1}{\overset{q}{\Sigma }}{b}_{\mathrm{nj}}\cos ({\omega }_{j}t+{\phi }_{\mathrm{nj}})\end{array}\right)$为多重潜周期多元自回 归混合模型的多重潜周期函数，b_ij为多重潜周期多元自回归混合模型的第i个特征因素的第j 个角频率对应的幅值，φ_ij为多重潜周期多元自回归混合模型的第i个特征因素的第j个角频 率对应的相位，i＝1,2，…,n，其计算方法见具体实施方式步骤三。

步骤四、特征因素时序预测；

1、特征因素多重潜周期多元自回归混合模型预测；

根据最佳预测原理(即最小均方误差预测原理)，本发明给出特征因素多重潜周期多元自 回归混合模型在t时刻向前预测l步的最佳预测值Y_t+l计算公式：

$Y_{t+l}=\underset{j=1}{\overset{p}{\Sigma }}{H}_j{Y}_{t+l-j}+C^{*}(t+l)---\left(5\right)$

2、多重潜周期多元自回归混合模型灰色预测。

对于特征因素多重潜周期多元自回归混合模型，若用基于监测数据处理检验后得到的待 预测时序直接对未来预测，由于监测数据样本大小的限制，其精度往往不理想。灰色理论中 的系统灰色预测嵌套法，是针对某种构造的系统灰色预测模型，将GM（1，1）模型嵌入GM （1，N）模型求解，以获得各行为变量的预测值。本发明基于系统灰色预测嵌套法，将特征 因素多重潜周期多元自回归混合模型每一步预测值作为待预测时序的最后一个数据，并删去 待预测时序中第一个数据，再重新建立特征因素多重潜周期多元自回归混合模型，使模型时 刻得到修正，以此类推，直至预测到所要求的步数为止，这种嵌套预测法让待预测时序得以 吐故纳新，通常要比直接预测精度高。

步骤五、特征因素分析；

多元时序分析不仅适于预测，也是一种因素分析的有效工具。本发明利用前面的多元时 序分析预测结果给出水华特征因素的因素分析方法。

1、预测误差计算；

根据步骤二的特征因素筛选结果，对符合多元周期平稳条件的特征因素组合的多元周期 平稳时序采用本发明中的步骤三和步骤四建立多重潜周期多元自回归混合模型并进行灰色预 测后，将表征因素的预测数据与实测数据比较可得到预测误差，预测误差的计算方法见具体 实施方式中的步骤五。

2、因素分析；

因素分析的原则为：

(1)表征因素预测误差越小的特征因素组合，其组合内影响因素与表征因素的相关性越高， 反之则越低；

(2)在一个特征因素组合中，若加入一个或多个影响因素后，其表征因素预测误差与加入 前相比显著变小，则说明加入的影响因素与表征因素显著相关，其显著性判断标准见具体实 施方式中的步骤五；

(3)在一个特征因素组合中，若加入一个或多个影响因素后，其表征因素预测误差与影响 因素加入前相比未显著变小，则说明加入的影响因素与表征因素基本无关。

通过这种因素分析，可以定性地得到水华特征因素中各影响因素与表征因素的相关程度， 即各影响因素与水华发生现象联系的密切程度。

本发明的优点在于：

1.本发明提出了基于延迟交叉相关矩阵法的多元周期平稳性检验方法，解决了多元时序 的多元周期平稳性检验问题。

2.本发明根据所提出的多元周期平稳性检验方法，给出了水华特征因素的筛选方法，为 水华特征因素多元时序的合理建模提供了必要的条件，也为水华特征因素的因素分析提供了 基础。

3.本发明采用多重潜周期模型描述水华形成过程中各特征因素的周期性波动，考虑了周 期性的环境季节变化对水华形成过程的影响，其与一元潜周期模型相比，还考虑了不同特征 因素间多重周期性变化的交互影响。

4.本发明采用多元自回归模型描述水华形成过程中各特征因素的随机性波动，其与一元 自回归模型相比，不但考虑了同一时刻多个特征因素之间的互相关性，同时考虑了各特征因 素在不同时刻的自相关性。

5.本发明提出了水华特征因素的多重潜周期多元自回归混合模型，该模型将多重潜周期 模型代入多元自回归模型的迭代运算，并对两种模型参数联合求解，从而考虑了水华形成过 程中特征因素的随机性和周期性波动在不同特征因素间的交互影响，与仅将两种模型简单相 加的方法相比，本发明对多特征因素的水华形成过程描述更加全面，从而提高了水华预测结 果的准确度。

6.本发明提出了多重潜周期自回归模型的灰色预测方法，克服了多重潜周期多元自回归 混合模型预测精度依赖于大样本量的局限，提高了模型的预测精度，从而提高了水华预测结 果的可信度。

7.本发明提出了基于多元时序分析的水华特征因素分析方法，根据本发明方法给出的水 华预测结果与实测数据的误差，可以得到各影响因素与水华发生现象联系的密切程度。

附图说明

图1是本发明基于多元周期平稳时序及灰色理论的水华预测及因素分析方法的流程图；

图2是本发明多重潜周期多元自回归混合模型灰色预测的流程图；

图3是本发明特征因素组合的因素分析的流程图；

图4是实施例1中pH值的原始时序及预处理后时序；

图5是实施例1中耗氧量的原始时序及预处理后时序；

图6是实施例1中水温的原始时序及预处理后时序；

图7是实施例1中浊度的原始时序及预处理后时序；

图8是实施例1中氨氮的原始时序及预处理后时序；

图9是实施例1中总氮的原始时序及预处理后时序；

图10是实施例1中总磷的原始时序及预处理后时序；

图11是实施例1中溶解氧的原始时序及预处理后时序；

图12是实施例1中叶绿素的原始时序及预处理后时序；

图13是实施例1中藻密度的原始时序及预处理后时序；

图14是实施例1中采用四种方法对叶绿素的预测数据及实测数据；

图15是实施例1中采用四种方法对藻密度的预测数据及实测数据；

图中曲线编号分别为：1-实施例1中各相应特征因素的原始时序，2-实施例1中各相 应特征因素的预处理后时序，3-实施例1中各相应表征因素的实测数据，4-实施例1中各 相应表征因素采用一元周期平稳时序分析的预测数据，5-实施例1中各相应表征因素采用一 元周期平稳时序分析及灰色理论的预测数据，6-实施例1中各相应表征因素采用多元周期平 稳时序分析的预测数据，7-实施例1中各相应表征因素采用多元周期平稳时序分析及灰色理 论的预测数据。

具体实施方式

下面将结合附图和实施例1对本发明作进一步的详细说明。

本发明是一种基于多元周期平稳时序分析及灰色理论的水华预测及因素分析方法。

为便于说明，本说明书中所有未经解释的名词及字母含义均由下述假设解释：与水华现 象有关的特征因素分为两种：一种是影响水华发生的特征因素，例如氮、磷、pH值、溶解氧、 水温、光照度等，以下叫做影响因素；另一种是表征水华发生的特征因素，例如叶绿素浓度、 藻密度等，以下叫做表征因素。以Y_t表示t时刻的特征因素向量；以y_it表示第i个特征因素 在t时刻的量值，共有n个特征因素，i＝1,2,…,n。

具体方法实施流程如图1所示，通过如下步骤实现:

步骤一、监测数据采集及预处理；

对可能影响或表征水华发生的多个特征因素进行监测。通常情况下，由监测设备采集到 的原始特征因素监测数据包含有异常采样点、均值漂移以及部分采样点数据缺失等等现象， 应对每个特征因素原始时序分别作预处理。本发明特征因素监测数据预处理方法为：

1、去除异常采样点，本发明去除绝对值大于3倍原始时序标准差的异常采样点。

2、去除均值漂移，本发明将原始时序减去其均值，使去除后的时序均值为零。

3、填补缺失数据，本发明以缺失数据的采样点前一个采样数据和缺失数据的采样点后一 个采样数据的均值来代替缺失数据的采样点。

步骤二、特征因素筛选；

对预处理后的特征因素多元时序进行多元周期平稳性检验，并根据多元周期平稳性检验 结果进行特征因素筛选。其多元周期平稳性检验方法如下：

记一组待检验的特征因素向量时序Y_t,t＝1,2,…,N的协方差矩阵为

Γ₀＝E[(Y_t-E(Y_t))(Y_t-E(Y_t))^T]

则Γ₀为n×n方阵。Y_t的交叉相关矩阵为

y₀＝D^-1Γ₀D^-1

其中由Γ₀的对角元素(Γ₁₁(0),…,Γ_nn(0))组成。则y₀中元素可 表示为

$y_{\mathrm{ij}}\left(0\right)=\frac{{\Gamma }_{\mathrm{ij}}\left(0\right)}{\sqrt{{\Gamma }_{\mathrm{ii}}\left(0\right){\Gamma }_{\mathrm{jj}}\left(0\right)}},i,j=\mathrm{1,2},···,n$

其中Γ_ij(0)表示Γ₀中第i行第j列元素。

将Y_t的延迟k阶的协方差矩阵记为

Γ_k＝E[(Y_t-E(Y_t))(Y_t-k-E(Y_t))^T],t＝1,2,…,N

其中Y_t-k表示t-k时刻的特征因素向量。Γ_k同样为n×n方阵。Y_t的延迟k阶交叉相关矩阵为

y_k＝D^-1Γ_kD^-1

则y_k中元素为

$y_{\mathrm{ij}}\left(k\right)=\frac{{\Gamma }_{\mathrm{ij}}\left(k\right)}{\sqrt{{\Gamma }_{\mathrm{ii}}\left(0\right){\Gamma }_{\mathrm{jj}}\left(0\right)}},i,j=\mathrm{1,2},···,n$

其中，Γ_ij(k)表示Γ_k中第i行第j列元素。

Γ_k可由

${\Gamma }_k=\frac{1}{N-k}\underset{t=k+1}{\overset{N}{\Sigma }}(Y_t-\bar{Y}){(Y_{t-k}-\bar{Y})}^T$

近似估计。其中$\bar{Y}=\frac{1}{N}\underset{t=k+1}{\overset{N}{\Sigma }}{Y}_t.$

当Y_t为多元周期性时序时，设其周期为T_i,i＝1,2,…,q，则有

$Y_t=Y_{t+T_i}$

则Y_t的延迟k+T_i阶的协方差矩阵为

${\Gamma }_{k+T_i}=E\left[(Y_t-E\left(Y_t\right)){(Y_{t-(k+T_i)}-E\left(Y_t\right))}^T\right]=E\left[(Y_t-E\left(Y_t\right)){(Y_{t-k}-E\left(Y_t\right))}^T\right]={\Gamma }_k$

可知，当Y_t为多元周期性时序时，其延迟k阶的协方差矩阵Γ_k也具有周期性。

因此若前k个y＝(y₀,y₁,…,y_k),k≤N-1具有周期性波动的性质，则Y_t的多元周期平稳性 检验通过。

步骤三、特征因素时序建模；

1、确定特征因素时序结构；

将t时刻的特征因素向量Y_t按公式(1)分解为周期项C_t和随机项R_t的叠加。

2、建立特征因素时序周期项模型；

本发明采用多重潜周期模型公式(2)描述特征因素时序周期项C_t。

3、建立特征因素时序随机项模型；

从Y_t减去C_t后，对随机项即Y_t的平稳随机性部分R_t采用多元自回归模型公式(3)描述。

4、建立特征因素多元周期平稳时序模型；

将周期项C_t的多重潜周期模型公式(2)与随机项R_t的多元自回归模型公式(3)结合为多重 潜周期多元自回归混合模型公式(4)。其结合方法如下：

设矩阵系数多项式

$H\left(z\right)=I-\underset{j=1}{\overset{p}{\Sigma }}{H}_j{z}^j$

其中z为矩阵系数多项式的根，I为n×n单位阵。

则根据公式(3)，有

$H\left(B\right)R_t=(I-\underset{j=1}{\overset{p}{\Sigma }}{H}_j{B}^j)R_t=R_t-\underset{j=1}{\overset{p}{\Sigma }}{H}_j{R}_{t-j}=E_t$

其中B为时间t的向后推移算子。

上式两边同乘H^-1(B)，则随机项R_t可表示为

R_t＝H^-1(B)Ε_t

且有

$\det (I-\underset{j=1}{\overset{p}{\Sigma }}{H}_j{z}^j)\ne 0,\left|z\right|\le 1$

将周期项C_t的多重潜周期模型的实值形式公式(2)写成复值的形式

$C_t=C\left(t\right)=\left(\begin{array}{c}\underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_{1j}\exp \left(i{\lambda }_{j}t\right)\\ \underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_{2j}\exp \left(i{\lambda }_{j}t\right)\\ .\\ .\\ .\\ \underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_{\mathrm{nj}}\exp \left(i{\lambda }_{j}t\right)\end{array}\right)=\underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_j\exp \left(i{\lambda }_{j}t\right)---\left(6\right)$

其中α_kj为多重潜周期模型复值形式的第k个特征因素的第j个角频率对应的幅值，λ_j为多重 潜周期模型复值形式的第j个角频率，${\alpha }_j=\left(\begin{array}{c}{\alpha }_{1j}\\ {\alpha }_{2j}\\ .\\ .\\ .\\ {\alpha }_{\mathrm{nj}}\end{array}\right).$

利用公式(6)可将公式(1)写成

$Y_t=C_t+R_t=\left(\begin{array}{c}\underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_{1j}\exp \left(i{\lambda }_{j}t\right)\\ \underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_{2j}\exp \left(i{\lambda }_{j}t\right)\\ .\\ .\\ .\\ \underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_{\mathrm{nj}}\exp \left(i{\lambda }_{j}t\right)\end{array}\right)+H^{-1}\left(B\right)E_t=\underset{j=1}{\overset{2q}{\Sigma }}{\alpha }_j\exp \left(i{\lambda }_{j}t\right)+H^{-1}\left(B\right)E_t---\left(7\right)$

在上式(7)两边同时乘上H(B)，取Η₀＝-I，得到

$H\left(B\right)Y_t=\underset{j=1}{\overset{2q}{\Sigma }}H\left(B\right){\alpha }_j\exp \left(i{\lambda }_{j}t\right)+E_t$

$=-\underset{j=1}{\overset{2q}{\Sigma }}\left(\underset{l=0}{\overset{p}{\Sigma }}{H}_l{\alpha }_j\exp \left(i{\lambda }_j(t-l)\right)\right)+E_t$

$=\underset{j=1}{\overset{2q}{\Sigma }}(-\underset{l=0}{\overset{p}{\Sigma }}{H}_l{\alpha }_j\exp \left(i{\lambda }_{j}l\right))\exp \left(i{\lambda }_{j}t\right)+E_t$

$=\underset{j=1}{\overset{2q}{\Sigma }}{\beta }_j\exp \left(i{\lambda }_{j}t\right)+E_t$

其中

β_kj≠0,k＝1,2,…,n表示多重潜周期多元自回归混合模型复值形式的第k个特征因素的第j个 角频率对应的幅值。

于是，可以将公式(7)写成

$Y_t=\underset{j=1}{\overset{p}{\Sigma }}{H}_j{Y}_{t-j}+\underset{j=1}{\overset{2q}{\Sigma }}{\beta }_j\exp \left(i{\lambda }_{j}t\right)+E_t---\left(8\right)$

公式(8)是多重潜周期多元自回归混合模型的复值形式。其对应的实值形式即为多重潜周期多 元自回归混合模型公式(4)。

实值形式的多重潜周期多元自回归混合模型的幅值b_ij和相位φ_ij,i＝1,2,…,n,j=1,2,…,q的 估计为：

如果λ_q＝π，取b_iq＝β_iq,φ_iq＝0,ω_q＝π，

如果λ_j∈(0,π)，取b_ij＝2|β_ij|,ω_j＝λ_j,φ_ij＝arg(β_ij)。

步骤四、特征因素时序预测；

1、特征因素多元周期平稳时序模型预测；

根据最小均方误差预测原理，得到特征因素多重潜周期多元自回归混合模型在t时刻向 前预测l步的最小均方误差预测公式(5)。

2、特征因素多元周期平稳时序灰色预测；

采用系统灰色预测嵌套法对特征因素多元周期平稳时序预测的流程如图2所示，具体步 骤如下：

(1)设定预测总步数H，令u＝1,2,…,H表示预测的步数。

(2)设{Y_ut},t＝1,2,…,N表示第u步预测时一组待预测的特征因素多元周期平稳时序，建 立多重潜周期模型

其中C_ut,C_u(t),q_u,a_uij,ω_uj,i=1,2,…,n为进行第u步预测时的周期项、多重潜周期模型的多 重潜周期函数及相应的模型参数。

(3)从Y_ut中减去周期项C_ut，对剩余的数据建立多元自回归模型

$R_{\mathrm{ut}}=\underset{j=1}{\overset{p_u}{\Sigma }}{H}_{\mathrm{uj}}{R}_{u(t-j)}+E_{\mathrm{ut}}---\left(10\right)$

其中R_ut,R_u(t-j),p_u,Η_uj为第u步预测时的随机项及相应的模型参数。Ε_ut为第u步预测时的白 噪声。

(4)建立多重潜周期多元自回归混合模型

$Y_{\mathrm{ut}}=\underset{j=1}{\overset{p_u}{\Sigma }}{H}_{\mathrm{uj}}{Y}_{u(t-j)}+C_u^{*}\left(t\right)+E_{\mathrm{ut}}---\left(11\right)$

其中p_u,Η_uj,为第u步预测时的模型参数及多重潜周期多元自回归混合模型的多重潜周 期函数。

(5)采用多重潜周期多元自回归混合模型的最小均方误差预测公式(5)对Y_uN进行向前一步 预测

$Y_{u(N+1)}=\underset{j=1}{\overset{p_u}{\Sigma }}{H}_{\mathrm{uj}}{Y}_{u(N+1-j)}+C_u^{*}(N+1)=Y_{N+u}---\left(12\right)$

其中Y_N+u为N+u时刻特征因素向量。至此完成特征因素多元周期平稳时序的一步预测。

(6)去掉该组时序的第一个数据Y_u1，加入时序的一步预测值Y_u(N+1)，得到一组新的时序 {Y_ut},t＝2,3,…,N+1，将该时序作为第u+1步预测时的待预测时序{Y_(u+1)t},t＝1,2,…,N。

(7)令u=u+1，则待预测时序表示为{Y_ut},t＝1,2,…,N，重复上述步骤(2)～(6)，可逐次预测 出特征因素多元时序{Y_N+u},u＝1,2,…,H。

步骤五、特征因素分析；

1、预测误差计算；

根据步骤二的特征因素筛选结果，对符合多元周期平稳条件的特征因素组合的多元周期 平稳时序采用本发明中的步骤三和步骤四建立多重潜周期多元自回归混合模型并进行灰色预 测后，将表征因素的预测数据与其同时刻的实测数据比较可得到预测误差，预测误差的表示 方式为平均误差绝对值，即

其中，H为预测数据总个数。

2、因素分析；

特征因素组合的因素分析的流程如图3所示，具体步骤为：

(1)计算表征因素预测误差，将表征因素组合作为已分析组合。

(2)从特征因素组合中任意选取一个未分析过的影响因素加入到已分析组合中，再对已分 析组合中的表征因素计算预测误差，若与加入前相比误差显著变小，则说明该影响因素与表 征因素相关性高，将其保留在已分析组合。否则说明该影响因素与表征因素相关性低，将其 从已分析组合中剔除。其显著性判断标准为

加入该影响因素后预测误差减小量≥加入该影响因素前预测误差的1%

若特征因素组合中含有多个表征因素，则以多个表征因素预测误差之和来进行显著性判断。

(3)重复上述步骤(2)，直至特征因素组合中所有影响因素均被分析过一遍，则已分析组合 中所保留的影响因素即为与表征因素或与水华发生关联密切的影响因素，将已分析组合中的 影响因素按其加入对表征因素预测误差改变的大小排序，定性地得到各影响因素与表征因素 的相关程度。

实施例1：

步骤一、特征因素监测数据采集及预处理；

对江苏省太湖2009年6月至2012年6月的10个水华特征因素进行监测，具体见表1。

表1水华特征因素监测名单

  名称   pH值   耗氧量   水温   浊度   氨氮   总氮   总磷   溶解氧   叶绿素   藻密度   单位   无   mg/L   ℃   NTU   mg/L   mg/L   mg/L   mg/L   mg/L   个/L

其中叶绿素和藻密度这两个特征因素为表征因素，其余的8个特征因素为影响因素。监 测设备一共记录了1104天的水华特征因素数据，其10个特征因素原始时序分别见图4至图 13中的灰色曲线1。对每个特征因素的原始时序分别作去除异常采样点、去除均值漂移以及 填补缺失数据的采样点的预处理，预处理后的特征因素时序见图4至图13中的黑色曲线2。

步骤二、特征因素筛选；

对预处理后的特征因素多元时序进行多元周期平稳性检验，并根据多元周期平稳性检验 结果进行特征因素筛选。特征因素筛选结果如下：

1、对所有特征因素采用自相关函数检验方法进行一元周期平稳性检验，可知浊度、氨氮、 总磷这三个特征因素的一元周期平稳性检验不通过，故应将其去掉不予考虑。

2、对一元周期平稳性检验通过的表征因素，即叶绿素浓度和藻密度，采用本发明给出的 多元周期平稳性检验方法，这两种表征因素的组合通过检验。

3、对表征因素的组合，分别加入通过一元周期稳定性检验的影响因素，即pH值、耗氧 量、水温、总氮、溶解氧这5个影响因素以组成新的特征因素组合，经检验可知由这5个影 响因素和2个表征因素组成的特征因素组合可以通过多元周期平稳性检验。

步骤三、特征因素时序建模；

取前1095天的特征因素时序用于建模，从而可根据前1095天数据对第1096天至第1104 天这9天的数据进行预测。

1、确定特征因素时序结构；

将特征因素时序按公式(1)分解为周期项和随机项的叠加。

2、建立特征因素时序周期项模型；

采用多重潜周期模型公式(2)描述特征因素时序周期项。

3、建立特征因素时序随机项模型；

从特征因素时序减去周期项后，对随机项采用多元自回归模型公式(3)描述。

4、建立特征因素多元周期平稳时序模型；

将周期项的多重潜周期模型与随机项的多元自回归模型结合为多重潜周期多元自回归混 合模型公式(4)。

步骤四、水华特征因素多元时序预测；

1、特征因素多元周期平稳时序模型预测；

根据最小均方误差预测原理，得到特征因素多重潜周期多元自回归混合模型的最小均方 误差预测公式(5)。

2、特征因素多元周期平稳时序灰色预测；

采用系统灰色预测嵌套法对前1095天的特征因素时序预测到第1104天的9步预测具体 步骤如下：

(1)设定预测总步数9，令u＝1,2,…,9表示预测的步数。

(2)设{Y_ut},t＝1,2,…,1095表示第u步预测时一组待预测的特征因素多元周期平稳时序， 建立多重潜周期模型公式(9)。

(3)从Y_ut中去掉周期项C_ut，对剩余的数据建立多元自回归模型公式(10)。

(4)建立多重潜周期多元自回归混合模型公式(11)。

(5)采用多重潜周期多元自回归混合模型的最小均方误差预测公式(5)对Y_u1095进行向前一 步预测公式(12)，得到第1095+u天的特征因素时序预测值Y_1095+u。至此完成特征因素时序的 一步预测。

(6)去掉该组时序的第一个数据Y_u1，加入时序的一步预测值Y_u1096，得到一组新的时序 {Y_ut},t＝2,3,…,1096，将该时序作为第u+1步预测时的待预测时序{Y_(u+1)t},t＝1,2,…,1095。

(7)令u=u+1，则待预测时序表示为{Y_ut},t＝1,2,…,1095，重复上述步骤(2)～(6)，可逐次预 测出特征因素多元时序{Y_1095+u},u＝1,2,…,9，即第1096天至第1104天的特征因素数据。

步骤五、特征因素分析；

1、预测误差计算；

根据步骤二的特征因素筛选结果，对符合多元周期平稳条件的特征因素组合采用本发明 方法建模预测后，其表征因素预测数据与实测数据比较可得到相应的预测误差。

2、因素分析如表2所示。

表2表征因素预测平均误差绝对值

从表2中可知，叶绿素和藻密度这两个表征因素的预测误差分别为0.3857和1.3910，预 测误差之和为1.7767；当加入影响因素溶解氧后，两个表征因素的预测误差分别为0.3055和 1.3153，其预测误差之和相比加入前减小了0.1559，减小量大于加入前预测误差之和的1%， 故可保留；当加入影响因素pH值后，表征因素预测误差之和相比加入前反而增加了0.0038， 故应剔除；当加入影响因素水温后，表征因素预测误差之和相比加入前总共减小了0.4398， 减小量大于加入前预测误差之和的1%，故可保留；当加入影响因素总氮后，表征因素预测误 差之和相比加入前总共减小了0.3191，减小量大于加入前预测误差之和的1%，故可保留；当 加入影响因素耗氧量后，表征因素预测误差之和相比加入前总共减小了0.04，减小量不到加 入前预测误差之和的1％，故应剔除。

因此，最终保留的影响因素为溶解氧、水温、总氮，说明这三个影响因素与叶绿素和藻 密度这两个表征因素的相关程度较高，即这三个影响因素与水华发生的关系最密切，而这三 个影响因素按其预测误差减小量从大到小，也即与水华发生的密切程度从高到低排序依次为 水温、总氮、溶解氧。

此外，为更直观地说明本发明采用多元周期平稳时序分析及灰色理论对水华预测的优势， 本发明另给出采用一元周期平稳时序分析、一元周期平稳时序分析及灰色理论、多元周期平 稳时序分析对水华特征因素多元时序的预测结果，以与本发明方法结果进行比较。

采用四种方法对第1096天至第1104天的水华特征因素多元时序进行预测，其2个表征 因素预测数据与实测数据分别如图14和图15所示，其7个特征因素的预测平均误差绝对值 如表3所示。

表3四种方法预测平均误差绝对值

可见，基于本发明方法的水华特征因素多元时序预测结果相比其它三种方法与实测结果更相 符、预测平均误差绝对值更小。