用于气象等复杂系统的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法

autoregressive model[j].ieee transactions on signal processing，2016，64(7):1759-1773.”提出一个方案将最近发展起来的反向时滞变量选择方法将滞后变量的选取与条件格兰杰因果关系指数(conditional granger causality index，cgci)结合，进而实现精确的预测建模。这些方法将经典的格兰杰因果关系分析扩展到多维时间序列，适用于高维数据的因果关系探究。
[0006]
本发明针对复杂高维系统的因果关系分析问题，提出一种基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，用于对气象污染等研究领域变量间的因果分析建模。相较于传统预测分析模型，本发明方法模型预测精度更高。


技术实现要素：

[0007]
本发明要解决的技术问题是，针对传统的格兰杰因果分析方法因为结构限制不能适用于高维时间序列，且不能对过程中变量间的变量关系进行准确地识别，对原有的格兰杰因果分析模型进行了拓展，提出了一种基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，实现对高维的数据进行准确的因果关系探究和对过程中的信息进行实时的展示。进而实现对复杂系统进行精准建模。本方法的目的在于拓展原始的方法适用于高维且能够展示更多的动态信息。
[0008]
为了实现上述目的，本发明采用的技术方案如下：
[0009]
一种基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，本发明面对气象等复杂系统，提出了探究系统间各变量与其主要污染物aqi之间的因果分析方法。首先，将采集到的数据进行预处理，对缺失数据采用平均值插补方法进行补全。然后，对补全的数据进行平稳性检验及处理，以满足建立模型的假设。之后，将数据进行归一化，以消除不同变量量纲带来的影响。最后，建立基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，实现准确探究变量间因果关系的目的，同时展现不同变量间的动态因果关系指数，以达到定量、明确地分析系统间各变量与aqi的联系，实现对aqi影响因素的分析。具体步骤如下：
[0010]
步骤1：获取空气质量指标观测数据，对多维时间序列数据进行预处理，采用单位根检验方法对时间序列数据进行平稳性检验，再对时间序列数据进行归一化处理。
[0011]
步骤2：采取孤立森林算法进行数据异常点的检测，并按照分数高低进行剔除处理。对于时间序列代表n
×
n维矩阵，共n维变量，时间序列长度为n。假设时间序列共有z个异常值，那么经过孤立森林算法处理过的时间序列变为
[0012][0013]
式中，h(k)＝ln(k)+ζ,ζ＝0.5772156649(欧拉常数)。其中，s(x,n)就是记录x个样本的训练数据构成树的异常指数，s(x,n)取值范围为[0,1]，越接近0，该数值为正常数值，越接近于1该数值为异常值。
[0014]
步骤3：选用最大相关最小冗余算法为建立的条件格兰杰因果分析模型选择合适的条件变量，将变量之间的最大相关性和最小冗余性结合起来，计算一个综合评分，选择具有最高综合评分的特征作为下一个被选择的特征，迭代计算直到从经过孤立森林算法处理过的时间序列x
*
中选取m个变量作为条件变量，定义为目标变量xi的特征子集综合评分计算方式为：
[0015][0016]
式中，m为特征子集数，sm为特征子集，xj,xk为特征子集变量。
[0017]
步骤4：选取上述步骤处理后数据集(包括目标变量、来源变量、条件变量)的最优解释向量，具体为：
[0018]
步骤4.1：首先，设置最大延迟数p
max
，代表模型选取终止。贝叶斯信息准则值(bic)为目标变量xi的方差。开始选择前，解释向量为空，定义maxlag为一个辅助数组，用于跟踪已搜索的每个变量的延迟数，因此对于每个变量，它最初被设置为零。
[0019]
步骤4.2：对于每个处理后数据集的变量，其延迟数增加1，将这个延迟变量添加到当前解释向量这样就形成候选解释向量。对于第一个周期，当前解释向量为空，候选解释向量为[x
i，t-1
，x
j，t-1
，x
1，t-1
，...，x
m，t-1
]，若来源变量xj包含在特征子集s中，则候选解释向量为[x
i,t-1
,x
1,t-1
,...,x
j,t-1
...,x
m,t-1
]。对于第l+1循环来说，解释向量有l个滞后变量，并让每个变量的延迟数为τ(l
m+2
)，如果还没有出现，则其值设置为零。则候选解释向量为或计算由候选m+2个候选解释向量组成的m+2阶动态回归模型的bic值。
[0020]
步骤4.3：比较当前bic值与下个循环增加解释向量的bic值。选择最小bic值的解释向量更新当前的解释向量组，然后进入下一步。如果增加候选解释向量后的bic值没有比当前bic值更小，保留当前解释向量，然后将延迟数增加1，然后进行下一步。
[0021]
步骤4.4：如果数据集中所有变量的延迟数达到最大延迟数p
max
(τ(l
m+2
)＝p
max
)，则终止，否则转到步骤4.1。循环结束后，算法给出模型的最优解释向量。
[0022]
步骤5：建立条件格兰杰因果分析模型，将上述解释向量带入模型，得到变量间的因果关系值，绘制条件格兰杰因果分析的结果得到变量间的因果关系热力图；
[0023]
步骤6：根据因果关系热力图选择与目标变量相关性最大的变量。将上述相关性变量带入回声状态网络模型得到目标变量预测值和实际值的对比情况以及预测误差。通过预测误差值比较对比方法与本发明方法的优劣。
[0024]
进一步的，所述的条件格兰杰因果分析模型在得到变量间因果分析结果基础上，进一步探究变量之间的动态特性；所述的条件格兰杰因果分析模型包括两个var模型，也称为动态回归模型：
[0025]
第一个var模型是不受限制的模型，称为unrestricted模型，该模型包括所有滞后
变量的最优解释向量，表示为:
[0026][0027]
式中，a
i,k
是unrestricted模型的参数；m+2是定义时滞的模型阶数；u
i,k
是unrestricted模型的预测误差。
[0028]
第二个var模型是restricted模型，该模型通过排除解释变量中的来源变量得到，表示为：
[0029][0030]
式中，b
i,k
是restricted模型的参数；e
i,t
是restricted模型的预测误差。
[0031]
变量i,j之间条件格兰杰因果关系指数(cgci)的定义：
[0032][0033]
式中，是根据u
i,t
和e
i,t
得到的残差方差。
[0034]
与现有技术相比，本发明具有以下有益效果：
[0035]
本发明着眼于对高维数据变量间的因果分析定量和准确分析。首先，本发明将其原始数据进行预处理，然后基于变量选择和特征选择方法建立了因果分析算法，这种算法可以区别于一般的基于向量自回归框架的因果分析方法，实现对高维数据的因果分析。然后，将其最优解释向量送入条件格兰杰因果分析模型中，得到定量的变量间因果关系值。除此之外，条件格兰杰因果分析方法还得到变量间的因果关系热力图，为探究变量间的因果关系提供更多的信息。即本发明不仅针对高维数据得到了定量的因果关系度量，还可以得到变量间的因果关系。
附图说明
[0036]
图1为格兰杰因果图展示，节点表示向量，箭头表示向量间存在因果关系；
[0037]
图2为本发明的流程图；
[0038]
图3为展示的各变量之间的动态关系图；图3(a)为变量aqi随时间变化图；图3(b)为变量pm2.5随时间变化图；图3(c)为变量pm10随时间变化图；图3(d)为变量so2随时间变化图；图3(e)为变量co随时间变化图；图3(f)为变量no2随时间变化图；图3(g)为变量o3随时间变化图；
[0039]
图4为各方法变量之间的因果关系图；图4(a)为fblg方法得到的变量间因果关系图；图4(b)为lassogc方法得到的变量间因果关系图；图4(c)为copulagc方法得到的变量间因果关系图；图4(d)为本发明方法得到的变量间因果关系图；
[0040]
图5(a1)为fblg对aqi的目标值与预测值的对比图；5(a2)为fblg对aqi的预测误差图；
[0041]
图5(b1)为lassogc对aqi的目标值与预测值的对比图；图5(b2)为lassogc对aqi的预测误差图；
[0042]
图5(c1)为copulagc对aqi的目标值与预测值的对比图；图5(c2)为copulagc对aqi的预测误差图；
[0043]
图5(d1)为本发明对aqi的目标值与预测值的对比图；图5(d2)为本发明对aqi的预测误差图。
具体实施方式
[0044]
以下结合具体实例和仿真图，对本发明进行进一步详细说明。
[0045]
本发明所使用硬件设备包括pc机器一台。
[0046]
图2为本发明提供的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法流程图，具体包括以下步骤：
[0047]
步骤1：本小节选用上海气象数据集进行仿真实验。该数据集记录了从2013年12月2日至2021年4月13日的每天数据，一共2690个样本,每个样本包含7维时间序列，其各变量含义编号详细描述见表1。然后对多维aqi及气象数据集的缺失值进行插补和异常值进行分析处理；采用单位根检验方法对数据进行平稳性检验，并根据检验结果对数据进行一次差分平稳化处理；对时间序列数据进行归一化处理；
[0048]
表1上海气象数据编号对应表
[0049][0050]
步骤2：采取孤立森林算法进行数据异常点的检测，并按照分数高低进行剔除处理。对于时间序列代表n
×
n维矩阵，共n维变量，时间序列长度为n。假设时间序列共有z个异常值，那么经过孤立森林算法处理过的时间序列变为这里剔除异常点比例为10％，处理后样本数为2421。
[0051][0052]
式中，h(k)＝ln(k)+ζ,ζ＝0.5772156649(欧拉常数)。s(x,n)就是记录x个样本的训练数据构成树的异常指数，s(x,n)取值范围为[0,1]，越接近0，该数值为正常数值，越接近于1该数值为异常值。
[0053]
步骤3：选用最大相关最小冗余算法为建立的条件格兰杰因果分析模型选择合适的条件变量，本方法条件变量选择为所有变量的30％，即本实验数据集有7个变量，对应条件变量为2个。将变量之间的最大相关性和最小冗余性结合起来，计算一个综合评分，选择具有最高综合评分的特征作为下一个被选择的特征，迭代计算直到从经过孤立森林算法处理过的时间序列x
*
中选取m个变量作为条件变量，定义为目标变量xi的特征子集综合评分计算方式为：
[0054][0055]
式中，m为特征子集数，sm为特征子集，xj,xk为特征子集变量。
[0056]
步骤4：选取上述步骤处理后数据集(包括目标变量、来源变量、条件变量)的最优解释向量，具体为：
[0057]
步骤4.1：首先，设置最大延迟数p
max
，代表模型选取终止。贝叶斯信息准则值(bic)为目标变量xi的方差。开始选择前，解释向量为空，定义maxlag为一个辅助数组，用于跟踪已搜索的每个变量的延迟数，因此对于每个变量，它最初被设置为零。
[0058]
步骤4.2：对于每个处理后数据集的变量，其延迟数增加1，将这个延迟变量添加到当前解释向量这样就形成候选解释向量。对于第一个周期，当前解释向量为空，候选解释向量为[x
i，t-1
，x
j，t-1
，x
1，t-1
，...，x
m，t-1
]，若来源变量xj包含在特征子集s中，则候选解释向量为[x
i,t-1
,x
1,t-1
,...,x
j,t-1
...,x
m,t-1
]。对于第l+1循环来说，解释向量有l个滞后变量，并让每个变量的延迟数为τ(l
m+2
)，如果还没有出现，则其值设置为零。则候选解释向量为或计算由候选m+2个候选解释向量组成的m+2阶动态回归模型的bic值。
[0059]
步骤4.3：比较当前bic值与下个循环增加解释向量的bic值。选择最小bic值的解释向量更新当前的解释向量组，然后进入下一步。如果增加候选解释向量后的bic值没有比当前bic值更小，保留当前解释向量，然后将延迟数增加1，然后进行下一步。
[0060]
步骤4.4：如果数据集中所有变量的延迟数达到最大延迟数p
max
(τ(l
m+2
)＝p
max
)，则终止，否则转到步骤4.1。循环结束后，算法给出模型的最优解释向量。
[0061]
步骤5：建立条件格兰杰因果分析模型，将上述解释向量带入模型，得到变量间定量的因果关系值，绘制条件格兰杰因果分析的结果得到变量间的因果关系热力图；
[0062]
步骤6：根据因果关系热力图选择与目标变量相关性最大的变量。本实验选择aqi作为目标变量，其中fblg选择的相关变量是3(pm10)，lassogc选择的相关变量是3、5(pm10、co)，copulagc选择的相关变量是7(o3)，fblg选择的相关变量是2、3、5、6(pm2.5、pm10、co、no2)。将上述相关性变量带入回声状态网络模型得到目标变量预测值和实际值的对比情况以及预测误差。通过预测误差值比较对比方法与本发明方法的优劣。
[0063]
通过对数据集时间序列进行绘制，得到如图3所示的其他变量与aqi之间的动态因果关系，从中得到变量间在时间分布上的动态信息。另外，对时间序列结果取均值，得到变量间定量的因果值，为后续的确定变量间因果关系提供量化的指标。从图3给出的动态曲线可以看出变量之间在时间分布上的动态关系。如果需要对某个特定时期进行研究，需要详细的从曲线中得到更多的信息；根据因果分析的结果得到目标变量最相关的变量，各方法因果分析热力图如图4所示。然后，利用回声状态网络建立预测模型，对因果关系的结果进行检验。其中，回声状态网络的储备池参数设置如下：储备池维数40、稀疏度0.05、谱半径0.99和连接输入权值0.05。一共进行10次独立实验，将10次的均值作为最终的结果，以消除
误差的影响。预测指标有均方根误差(rmse)、标准均方根误差(nrmse)、平均绝对百分误差(mape)和对称平均绝对百分比误差(smape)。四者定义如下：
[0064][0065][0066][0067][0068]
式中，xi和分别表示真实值和预测值，n表示样本数量。
[0069]
预测结果如图5所示。从图中可以看出，本发明的方法(fs-btscg)选择的变量子集能够使预测的aqi结果即使的反应真实的变化趋势。10次独立重复实验的误差如表2所示。与前后向套索格兰杰因果分析方法(fblg)和套索格兰杰方法(lassogc)以及科普拉格兰杰方法(copulagc)继续对比，本发明得到的误差指标都取得了最好的结果，这也进一步说明了本发明的有效性。
[0070]
表2各种方法预测结果
[0071][0072]
以上所述实例仅为表达该发明的实施方式，不能理解为对本发明使用范围的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出相应改进，这些均属于本发明的保护范围。

技术特征：
1.一种用于气象等复杂系统的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，所述条件格兰杰因果分析方法面对气象或其他复杂系统，提出探究系统间各变量与其主要污染物aqi之间的因果分析方法；其特征在于，包括以下步骤：首先，将采集到的数据进行预处理，对缺失数据采用平均值插补方法进行补全；然后，对补全的数据进行平稳性检验及处理，以满足建立模型的假设；之后，将数据进行归一化，以消除不同变量量纲带来的影响；最后，建立基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，实现准确探究变量间因果关系的目的，同时展现不同变量间的动态因果关系指数，以达到定量、明确地分析系统间各变量与aqi的联系，实现对aqi影响因素的分析。2.根据权利要求1所述的一种用于气象等复杂系统的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，其特征在于，具体步骤如下：步骤1：获取空气质量指标观测数据，对多维时间序列数据进行预处理，并对时间序列数据进行平稳性检验，再对时间序列数据进行归一化处理；步骤2：采取孤立森林算法进行数据异常点的检测，按照分数高低进行剔除处理；对于时间序列代表n
×
n维矩阵，共n维变量，时间序列长度为n；假设时间序列共有z个异常值，那么经过孤立森林算法处理过的时间序列变为序列共有z个异常值，那么经过孤立森林算法处理过的时间序列变为式中，h(k)＝ln(k)+ζ,ζ＝0.5772156649(欧拉常数)；其中，s(x,n)就是记录x个样本的训练数据构成树的异常指数，s(x,n)取值范围为[0,1]，越接近0，该数值为正常数值，越接近于1该数值为异常值；步骤3：选用最大相关最小冗余算法为建立的条件格兰杰因果分析模型选择合适的条件变量，将变量之间的最大相关性和最小冗余性结合起来，计算一个综合评分，选择具有最高综合评分的特征作为下一个被选择的特征，迭代计算直到从经过孤立森林算法处理过的时间序列x
*
中选取m个变量作为条件变量，定义为目标变量x
i
的特征子集综合评分计算方式为：式中，m为特征子集数，s
m
为特征子集，x
j
,x
k
为特征子集变量；步骤4：选取上述步骤处理后数据集的最优解释向量，其中数据集中包括目标变量、来源变量、条件变量；步骤5：建立条件格兰杰因果分析模型，将上述解释向量带入模型，得到变量间的因果关系值，绘制条件格兰杰因果分析的结果得到变量间的因果关系热力图；步骤6：根据因果关系热力图选择与目标变量相关性最大的变量；将上述相关性变量带
入回声状态网络模型得到目标变量预测值和实际值的对比情况以及预测误差；通过预测误差值比较对比方法与本发明方法的优劣。3.根据权利要求2所述的一种用于气象等复杂系统的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，其特征在于，所述的步骤4具体如下：步骤4.1：首先，设置最大延迟数p
max
，代表模型选取终止；贝叶斯信息准则值bic为目标变量x
i
的方差；开始选择前，解释向量为空，定义maxlag为一个辅助数组，用于跟踪已搜索的每个变量的延迟数，因此对于每个变量，它最初被设置为零；步骤4.2：对于每个处理后数据集的变量，其延迟数增加1，将这个延迟变量添加到当前解释向量这样就形成候选解释向量；对于第一个周期，当前解释向量为空，候选解释向量为[x
i，t-1
，x
j，t-1
，x
1，t-1
，...，x
m，t-1
]，若来源变量x
j
包含在特征子集s中，则候选解释向量为[x
i,t-1
,x
1,t-1
,...,x
j,t-1
...,x
m,t-1
]；对于第l+1循环来说，解释向量有l个滞后变量，并让每个变量的延迟数为τ(l
m+2
)，如果还没有出现，则其值设置为零；则候选解释向量为或计算由候选m+2个候选解释向量组成的m+2阶动态回归模型的bic值；步骤4.3：比较当前bic值与下个循环增加解释向量的bic值；选择最小bic值的解释向量更新当前的解释向量组，然后进入下一步；如果增加候选解释向量后的bic值没有比当前bic值更小，保留当前解释向量，然后将延迟数增加1，然后进行下一步；步骤4.4：如果数据集中所有变量的延迟数达到最大延迟数p
max
(τ(l
m+2
)＝p
max
)，则终止，否则转到步骤4.1；循环结束后，算法给出模型的最优解释向量。4.根据权利要求2所述的一种用于气象等复杂系统的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，其特征在于，所述的条件格兰杰因果分析模型在得到变量间因果分析结果基础上，进一步探究变量之间的动态特性；所述的条件格兰杰因果分析模型包括两个var模型，也称为动态回归模型：第一个var模型是不受限制的模型，称为unrestricted模型，该模型包括所有滞后变量的最优解释向量，表示为:式中，a
i,k
是unrestricted模型的参数；m+2是定义时滞的模型阶数；u
i,k
是unrestricted模型的预测误差；第二个var模型是restricted模型，该模型通过排除解释变量中的来源变量得到，表示为：式中，b
i,k
是restricted模型的参数；e
i,t
是restricted模型的预测误差；变量i,j之间条件格兰杰因果关系指数(cgci)的定义：
式中，是根据u
i,t
和e
i,t
得到的残差方差。

技术总结
本发明提出一种用于气象等复杂系统的基于变量选择和反向时滞特征选择的条件格兰杰因果分析方法，属于数据挖掘技术领域。本发明先将采集到的数据进行预处理，对缺失数据采用平均值插补方法进行补全，并对数据进行平稳性检验及处理，以满足建立模型的假设。之后将数据进行归一化，以消除不同变量量纲带来的影响。最后建立基于变量选择和反向时滞特征选择的条件格兰杰因果分析模型，实现准确探究变量间因果关系的目的，同时展现不同变量间的因果关系指数，以达到定量、准确地分析系统间各变量间的因果关系。本发明能够实现对复杂系统进行精准建模，目的在于拓展原始的方法适用于高维且能够展示更多的动态信息。维且能够展示更多的动态信息。维且能够展示更多的动态信息。


技术研发人员：李蒙 那晓栋 韩敏
受保护的技术使用者：大连理工大学
技术研发日：2023.08.01
技术公布日：2023/10/15