CN112036460A - 一种识别量化控制泉流量潜在因素的方法 - Google Patents

Abstract

本发明公开了一种识别量化控制泉流量潜在因素的方法，包括选定研究的喀斯特地貌区域，获取区域内多组泉流量数据及驱动因素数据；异常值检测，如果存在异常值，则对异常值进行处理；否则直接执行下一步；利用动态因素分析提取多组泉流量数据公共的潜在因素；对时间序列数据进行小波分解；进行K‑均值聚类分析；基于小波分解计算潜在因素和驱动因素的多分辨率互相关性；分析潜在因素和驱动因素之间的小波相干性；对潜在因素进行识别；利用对应尺度上的有效动态效率表示不同因素对泉流量的相对贡献率。本发明开发了具有在时频域表征非平稳信号动态结构特征的技术，适用于量化人类活动和气候改变等主要影响因素对泉流量在不同尺度上的影响和贡献。

Description

一种识别量化控制泉流量潜在因素的方法

技术领域

本发明涉及复杂的水文过程探究技术，具体涉及一种识别和量化控制泉流量的潜在因素的方法。

背景技术

喀斯特地貌含水介质结构的多重性和形态的高度复杂性，使喀斯特区域的水文效应和水资源的评价与预测一直是水文研究的重点和难点。识别和量化控制泉流量的潜在因素，并且揭示不同尺度中的水文过程特征，可以改善对水文系统中复杂过程的理解。为水文过程的预测模型的建立提供重要参考，这对于解决喀斯特地区的水资源可持续利用问题具有重要意义。

目前已经发展了许多方法像水文模型，数值方法和时间序列分析对复杂的水文过程进行探索，但很少有方法可以定量的评价动态因素对泉流量的贡献(见文献1：蒙海花，王腊春，2010.岩溶流域水文模型研究进展[J].地理科学进展,29,1311-1318；见文献2:Berendrecht,W.L.,van Geer,F.C.,2016.A dynamic factor modeling framework foranalyzing multiple groundwater head series simultaneously.Journal ofHydrology 536,50-60)。并且由于地质介质的非均质性，水文循环的复杂性以及它们在野外条件下的相互作用，识别泉流量变化的驱动力和定量评价驱动力的影响也是十分困难的。

动态因素分析可以通过一些潜在因素描述对泉流量变化的影响，以反映泉流量的动态结构特征。动态因素分析可以减少大型数据的维数，并且对于相互依赖的非平稳数据同样适用。在水文地质学中，动态因素分析已经用来识别地下水位波动中的补给和抽取的发展趋势。研究者对动态因素分析在水文过程中的运用做出了很大的努力和贡献，但是仍然无法量化不同的潜在因素对水文过程驱动的影响。

发明内容

发明目的：本发明的目的是提供一种可以有效识别驱动泉流量变化的潜在因素并且可以量化不同潜在因素对泉流量波动的影响的识别和量化控制泉流量的潜在因素的方法，此外该方法可以在时频域表征时间序列之间的关系特征。

技术方案：本发明的一种识别量化控制泉流量潜在因素的方法，包括以下步骤：

S1、选定研究的喀斯特地貌区域，获取区域内多组泉流量及驱动因素数据，并对获取的数据进行标准化；

S2、异常值检测，判断是否存在异常值，如果存在异常值，则对异常值进行处理；如果不存在异常值，则直接执行步骤S3；

S3、利用动态因素分析提取多组泉流量数据公共的潜在因素，依据在不同动态因素模型中的潜在因素对泉流量的负载大小选择主要的潜在因素；动态因素分析可以揭示时间序列的动态结构，该方法对具有记忆特征的非平稳时间序列仍然有效；

S4、对时间序列数据进行小波分解；小波分解可以将时间序列分解到不同的时间尺度上，揭示不同时间尺度上的波动特征。小波分解是多分辨率互相关性分析的基础，并且小波能量分布表征了时间序列在不同尺度上的受控制因素的影响程度，是表征影响因素对响应因素相对贡献的重要参考指标；

S5、基于最大尺度上的近似小波系数进行K-均值聚类分析；潜在因素的最大尺度近似尺度的小波系数不仅保留了潜在因素的波动特征，还可以降低时间序列的维数提高聚类效率；

S6、基于小波分解计算潜在因素和驱动因素的多分辨率互相关性，用多分辨率互相关系数绝对值的最大值表征两个时间序列的相似程度和联系，从而对潜在因素进行识别；需要注意的是，这里只考虑时间序列的相关程度，忽略了正负性，即考虑相关系数的绝对值；

S7、分析潜在因素和驱动因素之间的小波相干性；揭示潜在因素和驱动因素时频域上的关系特征，进一步识别潜在因素；并且小波相干性分析可以将时间序列之间的延迟反应和关系特征分别进行表征；

S8、利用对应尺度上的小波能量分布和负载的乘积表示不同的因素对泉流量的相对贡献。

进一步的，步骤S1中数据标准化是指减去均值除以标准差对数据进行标准化。

进一步的，步骤S2中异常值是指时间序列中的值与平均值的偏差超过两倍标准差的值，使用方框图扫描泉流量和驱动因素时间序列是否存在异常值，对异常值进行处理的具体方法为：利用异常值前后对应值的均值代替。

进一步的，步骤S3中动态因素分析具体为：

F_j(t)＝F_j(t-1)+η_j (2)；

其中，t表示观测时间，m表示考虑的驱动力的数量，Y_i(t)、μ_i(t)、ε_i分别表示第i个观测时间序列的泉流量时间序列、水平参数和具体误差，F_j(t)是第j个共同的潜在因素，α_i,j是第j个潜在因素对第i个泉流量时间序列的负载，是量化潜在因素对泉流量影响的重要指标，η_j是状态方程公式(2)的剩余量；并且假设ε_i和η_j为白噪声，分别用协方差矩阵R和Q来表示；动态因素被R的两种结构抽取：1)方差相等的对角阵；2)方差不相等的对角阵。公式(1)和(2)中的所有参数和动态因素模型通过极大似然法和结合卡尔曼滤波器和期望最大化算法进行测试；最优动态因素分析模型(潜在因素数量m和协方差矩阵R)是根据Akaike信息准则的小样本校正值AICc进行选择；根据不同动态因素模型中的潜在因素对泉流量的负载选择较为重要的潜在因素。

进一步的，步骤S4具体为：

小波分析使用下面的公式探索信号和小波函数之间的关系：

其中，ψ是母小波，k是时间序列的长度，a和b分别是分解尺度和平移系数；ψ_a,b(t)表示在分解尺度a，平移系数为b时的母小波。如果a和b以2进制间隔不连续取样，即a＝2^-p，b＝2^-pq，p和q分别表示小波分解尺度和平移系数，则公式(3)是离散小波变换对应的母小波，表示为：

ψ_p,q(t)＝2^p/2ψ(2^pt-q) (4)；

其中，ψ_p,q(t)表示在分解尺度p，平移系数为q时的母小波，用来表征高频信号的特征；

用下面公式表征低频信号域特征的小波函数即父小波表示为：

φ_p,q(t)＝2^p/2φ(2^pt-q) (5)；

为了表征信号整个区域的特征，则离散小波变换表示成：

其中，f(t)表示目标时间序列；k表示时间序列的长度；J是最大的分解尺度；θ_J,q是在分解尺度J上近似小波系数，用于表示低频区域的信息，λ_p,q表示高频区域的小波系数；所有q对应的θ_J,q和λ_p,q分别表示成近似分量A_J和细节分量D_p；并且J根据下面的公式选择：

进一步的，步骤S6中互相关性系数表示为：

其中，r_xy(τ)表示互相关系数，τ是滞后时间，σ_x和σ_y分别是时间序列x和y的标准差，

和

表示x和y的均值，C_xy(τ)表示x和y的互协方差函数。

利用潜在因素的最大分解尺度上的近似小波系数进行聚类，不仅保留了潜在因素的波动特征，而且通过降低维数提高聚类的效率。

进一步的，步骤S7中首先对长度为k的潜在因素时间序列x和驱动因素时间序列y进行连续小波变换；为了下面方便表示，将x和y分别记为x_i,i＝1,2,...,k和y_i,i＝1,2,...,k；x_i和y_i为等间隔(△)取样，则在时间t_j＝j·△和尺度a上的连续小波变换W_i ^x(a)和W_i ^y(a)分别表示为：

其中“*”表示共轭算子。

潜在因素时间序列x_i和驱动因素序列y_i的小波相干性系数

表示为：

其中，a表示尺度，S表示光滑算子，W_i ^x(a)和W_i ^y(a)表示时间序列x和y的连续的小波变换，

“*”表示共轭算子。

并且时间序列的相位特征表示为：

Φ_i(a)＝tan^-1(Im(a^-1W_i ^xy(a))/Re(a^-1W_i ^xy(a))) (13)；

其中，Im表示虚部，Re表示实部。

进一步的，步骤S8中不同的因素对泉流量的相对贡献率D_ef为：

D_ef＝E_DF*DF_loading (14)；

其中，E_DF和DF_loading分别表示在对应尺度上能量分布和动态因素对时间序列的负载。

发明原理：动态因素分析可以有效揭示时间序列的动态结构，识别潜在因素。小波分析是分解和量化时间序列的理想工具，本发明是结合动态因素分析和小波分析的优点发展了一种识别和量化潜在因素控制泉流量的技术。

有益效果：与现有技术相比，本发明利用动态因素分析揭示时间序列的动态结构特征，抽取泉流量时间序列的公共潜在因素，利用多分辨率互相关性和小波相干性分析有效识别潜在因素及其对泉流量影响的作用尺度。最后利用有效动态效率量化主要潜在因素对泉流量的相对贡献。结合动态因素分析和小波分析有效的解决现有方法无法定量评价潜在因素对泉流量影响的缺陷。第二，为研究复杂的水文过程提供了方便可靠的方法，具有重要的理论和工程意义。

附图说明

图1为本发明所述方法的流程图；

图2a为在组1中潜在因素对不同泉流量时间序列的负载图，2b为在组2中潜在因素对不同泉流量时间序列的负载图，2c为在组3中潜在因素对不同泉流量时间序列的负载图；

图3为潜在因素与潜在因素的最高负载对应的泉流量波动特征；

图4为潜在因素的小波聚类图；

图5a为组1中潜在因素与潜在因素最高负载对应的泉流量序列的小波相干性图；5b为组2中潜在因素与潜在因素最高负载对应的泉流量序列的小波相干性图；5c为组3中潜在因素与潜在因素最高负载对应的泉流量序列的小波相干性图。

具体实施方式

下面结合说明书附图，对本发明作进一步的说明。

如图1所示，本发明提出的识别量化控制泉流量潜在因素的方法，包括如下步骤：

S1、选定所研究的喀斯特地貌区域，获取区域内多组泉流量及驱动因素数据，并对获取的数据进行标准化。

本实施例选定的要研究的喀斯特区域为娘子关泉流域，在娘子关泉流域水文过程的主要驱动因素包括降水，气候变化和人类活动；数据标准化是指减去均值除以标准差。

S2、异常值检测，使用方框图扫描泉流量和驱动因素时间序列数据，判断是否存在异常值，如果存在异常值，则用异常值前后对应值的均值代替；如果不存在异常值，则直接执行步骤S3。

所述异常值是指时间序列中的值与平均值的偏差超过两倍标准差的值。

S3、利用动态因素分析提取多组泉流量数据公共的潜在因素，依据在不同动态因素模型中的潜在因素对泉流量的负载大小选择主要的潜在因素；

动态因素分析可以揭示时间序列的动态结构，该方法对具有记忆特征的非平稳时间序列仍然有效。

动态因素分析具体为：

F_j(t)＝F_j(t-1)+η_j (2)；

其中，t表示观测时间，m表示考虑的驱动力的数量，Y_i(t)、μ_i(t)、ε_i分别表示第i个观测的泉流量时间序列、水平参数和具体误差，F_j(t)是第j个共同的潜在因素，α_i,j是第j个潜在因素对第i个泉流量序列的负载，是量化潜在因素对泉流量影响的重要指标。潜在η_j是状态方程公式(2)的剩余量；并且假设ε_i和η_j为白噪声，分别用协方差矩阵R和Q来表示。动态因素被R的两种结构抽取：1)方差相等的对角阵；2)方差不相等的对角阵。公式(1)和(2)中的所有参数和动态因素模型通过极大似然法和卡尔曼滤波器并且结合期望最大化算法进行测试。最优动态因素模型(潜在因素的数量(m)和协方差矩阵(R))是根据Akaike信息准则的小样本校正值(AICc)进行选择。根据不同动态因素模型的潜在因素对泉流量的负载选择较为重要的潜在因素。

S4、对潜在因素，泉流量和驱动因素时间序列数据进行小波分解。小波分解可以将时间序列分解到不同的尺度水平上，揭示不同尺度下的波动特征；并且潜在因素的小波能量分布表征了在不同尺度上泉流量受驱动因素的影响程度，是表征潜在因素对泉流量相对贡献的重要参考指标。

小波分析分解泉流量，潜在因素和驱动因素时间序列，使用下面的方程探索信号和小波函数之间的关系：

其中，ψ是母小波，k是时间序列的长度，a表示分解尺度，b是平移系数，反映时间上的平移，ψ_a,b(t)表示在分解尺度a，平移系数为b时的母小波；如果这两个系数是连续的，公式(3)是连续小波变换对应的母小波；如果a和b以2进制间隔不连续取样(a＝2^-p，b＝2^-pq，p和q分别表示小波分解尺度和时间上的平移系数)，则公式(3)是离散小波变换，其对应的母小波表示为：

ψ_p,q(t)＝2^p/2ψ(2^pt-q) (4)；

其中，ψ_p,q(t)为在分解尺度p和平移因子q下的母小波，主要表征高频信号的特征。用下面公式表征低频信号域特征的小波函数即父小波表示为：

φ_p,q(t)＝2^p/2φ(2^pt-q) (5)；

φ表示父小波，φ_p,q(t)表示在分解尺度为p，平移因子为q时的父小波。

为了表征泉流量，潜在因素和驱动因素时间序列在整个区域的特征，则离散小波变换可以表示为：

其中，f(t)表示泉流量，潜在因素和驱动因素时间序列；k表示时间序列的长度；J是最大的分解尺度；在对应的分解尺度p和时间平移因子q上，θ_J,q表示在分解尺度J上近似小波系数，用于表示低频区域的信息，λ_p,q表示高频区域的小波系数，ψ_p,q(t)为在分解尺度p和平移因子q下的母小波，φ_p,q(t)表示在分解尺度J，平移因子q上的父小波。所有q对应的θ_J,q和λ_p,q分别表示成近似分量A_J和细节分量D_p。并且J根据下面的公式选择：

其中f₁表示母小波滤波器长度。

S5、基于最大尺度的近似小波系数对潜在因素进行K-均值聚类，潜在因素的最大尺度的小波系数不仅可以表征潜在因素的波动特征，而且通过减少潜在因素时间序列的维数来提高聚类效率。

S6、基于小波分解计算潜在因素和驱动因素的多分辨率互相关系数，并获得不同分辨率下的互相关系数绝对值的最大值，用来表征潜在因素和驱动因素时间序列之间的相似程度。互相关性系数表示为：

其中，r_xy(τ)表示时间序列x和y互相关系数，τ是滞后时间，σ_x和σ_y分别是潜在因素时间序列x和驱动因素时间序列y的标准差，

和

分别表示x和y的均值，x_t和y_t分别表示x和y在t时刻的值，C_xy(τ)表示x和y的互协方差函数。

用多分辨率互相关系数绝对值的最大值表征两个时间序列的相似程度和联系，这里我们主要考虑的是两个时间序列的相关性的大小，忽略了关系的正负性。

S7、分析潜在因素和驱动因素的小波相干性，揭示他们时频域上的关系特征，目的是为了对多分辨率互相关方法识别潜在因素进行补充。并且小波相干性分析可以将时间序列之间的延迟反应和关系特征进行区分表征。在引入小波相干性之前，首先介绍连续的小波变换。对长度为k的潜在因素时间序列x和驱动因素时间序列y进行连续小波变换。为了下面公式的方便表达，把x和y进一步分别记为x_i,i＝1,2,...,k和y_i,i＝1,2,...,k。x_i和y_i为等间隔(记为△)取样，则在时间T_j＝j·△和尺度a上的连续小波变换W_i ^x(a)，W_i ^y(a)分别表示为：

其中“*”表示共轭算子。

潜在因素时间序列x_i和驱动因素时间序列y_i的小波相干性系数

表示为：

其中，a表示尺度，S表示光滑算子，W_i ^x(a)和W_i ^y(a)分别表示时间序列x和y的连续小波变换，

“*”表示共轭算子。

并且潜在因素和驱动因素时间序列的相对相位表示为：

Φ_i(a)＝tan^-1(Im(a^-1W_i ^xy(a))/Re(a^-1W_i ^xy(a))) (13)；

其中，Im表示虚部，Re表示实部。

S8、根据小波能量分布和负载特征表示不同的潜在因素对泉流量的相对贡献。

D_ef＝E_DF*DF_loading (14)；

其中，E_DF和DF_loading分别表示在对应尺度上能量分布和动态因素对时间序列的负载，D_ef为有效动态效率。

实施例：

该实施例的数据来自中国北方最大的喀斯特泉—娘子关泉。娘子关泉群位于阳泉市平定县娘子关镇，由几十个单泉组成。娘子关泉接收大约7394平方千米集水区的补充。由于娘子关泉流域是北方地区的重工业区，所以在这个实施例中，除了考虑对泉流量有重要影响的气候因素外，还考虑了人类活动的影响。因此考虑的主要驱动因素具体包括季风(ISM：印度夏季风和WNPM：西太平洋北部季风)，降水，地下水补充和人类活动。利用滑动窗对获得的1959-2007年泉流量数据进行分组，一共选取了3组泉流量数据进行具体的分析，如表1所示。

表1.泉流量时间序列的详细划分方式

表2展示了组1在两种协方差矩阵(R)结构下的动态因素分析结果。结果显示，在方差相等的对角矩阵下获得的4个潜在因素具有最小的AICc值。所以组1中，在协方差矩阵为方差不相等的对角阵条件下获得的4个潜在因素最好的表征泉流量波动的特征。并且在协方差矩阵为方差不相等的对角阵的条件下，对于不同潜在因素的个数(2，3和4)，潜在因素对泉流量T1和T7的最高负载出现次数最多，并且DF41和DF44分别对T1和T7具有最高的负载，所以DF41和DF44被选择为组1的主要潜在因素。利用相似的方法，组2和组3泉流量时间序列分别获得了4个潜在因素DF41-DF44，并且组2的主要潜在因素为DF41和DF42，组3的主要潜在因素为DF43和DF44。

表2.组1中不同协方差矩阵下的动态因素分析结果

表2中DFmn表示m个潜在因素中的第n个潜在因素。T1-T8表示时间段1-8的泉流量时间序列。主要动态因素DF及其对应的协方差矩阵以粗体显示。

图2a-2c展示了三组泉流量时间序列的潜在因素对泉流量时间序列的负载。在组1-3中，动态因素对泉流量的负载都小于等于0.47。由于喀斯特地区地质的复杂性，这种载荷特征可能是由于地质介质对水文过程的阻尼作用导致了潜在因素对泉流量时间序列的相对较小的负载，但是图2a-2c中的结果足以说明泉流量时间序列存在公共的潜在因素。每组中的潜在因素DF44对T1，T2和T3的负载为0，考虑到娘子关泉流域是在1979年开始开发地下水，并且在1980s后加速了对地下水的开发，所以每组中的DF44可能表示人类活动(地下水抽取)或者受人类活动影响下的其他潜在因素对泉流量的影响。此外，组1和组3中的DF41对T1具有最大的负载而对于其他时间段的负载较小，所以这两个潜在因素主要表示地下水开采前潜在因素对泉流量时间序列的影响。

图3描述了每组中潜在因素与潜在因素最大负载所对应的泉流量时间序列的波动特征。每一组中的DF42均呈线性下降趋势，这可能反映了由于重工业区(煤炭开采，化学工程，发电和冶金)的建设而导致的泉流量的下降(区域地下水的补给的影响)。在娘子关流域，7月至9月的降水占全年降雨量的60-70％。在组1中，DF43的年最大值集中在8至10月份，考虑到泉流量高峰期至少比降雨高峰期延迟约一个月，所以DF43可能代表降水和泉流量之间复杂的相互作用。每组中的DF41的年最小值发生在5月至9月，因为考虑到这段时期相对充足的降雨，灌溉抽水量减少。因此，DF41可能反映了人类活动对泉流量的季节性影响。基于图2的分析，并且根据潜在因素DF44与DF44最大负载泉流量时间序列的波动特征，初步判断组1和组2中的DF44可能表示人类活动的影响，组3中的DF44表示地下水开采后的局部地下水的补充的影响。

表3展示了在组1中4个潜在因素与潜在因素最高负载所对应时间段的驱动因素的小波能量分布。潜在因素的能量分布主要集中在A4(16个月)上，说明在大尺度上的水文过程可以解释泉流量波动的大部分特征。降水和季风在D3(8个月)和A4(16个月)上具有相对大的能量分布，并且降水在较短的时间内(小于4个月)具有相对较大的能量分布，这说明降水对季节周期的泉流量波动的影响更大。在组2和组3中可以获得类似的结论。为了更准确的识别潜在因素，对3组泉流量潜在因素的A4对应的小波系数进行K-均值聚类，获得2个聚类组。图4展示了聚类组1和聚类组2中的潜在因素具有不同的初始值和波动趋势。

表3.潜在因素和驱动力因素的小波能量分布

表4展示了潜在因素和泉流量时间序列在不同分辨率水平下的互相关系数绝对值的最大值，组1中的DF41(4-16个月)，组3中的DF41(2-16个月)和DF44(2-16个月)与泉流量显示出最大的互相干系数，考虑到地质介质对局部地下水的补给具有相对较小的阻尼作用，并且根据图2的结果，说明这三个因素主要表示局部地下水补给对泉流量的影响。可以根据类似的方法，对其他因素进行识别。

表4.潜在因素和泉流量时间序列在多分辨率水平上互相关系数绝对值的最大值

表4中主要潜在因素加粗倾斜表示，最大的相关性系数加粗强调。

图5a-5c表示组1中主要潜在因素(DF41和DF44)与最大潜在因素负载所对应的泉流量时间序列之间的小波相干性结果。对于DF41(代表局部地下水的补充)，与泉流量之间显示出很强的相关性。并且表征了与泉流量的同相特征。DF44在2-4个月和≥16个月的周期上表示诸如煤矿开采，发电和化学工程等人类活动对泉流量时间序列的影响。同样的，对组2和组3进行分析。通过多分辨率互相关系数和小波相干性分析，可以对主要的潜在因素进行有效识别。并且通过有效动态效率表示主要潜在因素对泉流量的相对贡献。则主要的结论如下所示：

第一组：DF41表示地下水开发前的局部地下水的补给(4-16个月)，|R_max|＝1.00，D_ef≤42.61％；DF44表示人类活动的影响(2-4个月)，|R_max|＝0.67-0.89，D_ef≤0.36％和(≥16个月)，|R_max|＝0.91-0.92，D_ef≤27.16％。

第二组：DF41表示降水和径流的影响(4-16个月)，|R_max|＝0.88-0.99，D_ef≤29.87％；DF42表示区域地下水补充的影响(4-16个月)，|R_max|＝0.93-1.00，D_ef≤29.90％。

第三组：DF43代表季风的影响(8-16个月)，|R_max|＝0.74-0.78，D_ef≤3.02％；DF44表示地下水开采后地下水的补充的影响(2-16个月)，|R_max|＝1.00，D_ef≤41.00％。

因此，本发明提出的识别和定量评价控制泉流量潜在因素的分析方法，利用动态因素分析提取泉流量公共的潜在因素。多分辨率互相关系数用来揭示主要潜在因素与驱动力之间的相似程度，从而对潜在因素进行识别。小波相干性分析用来揭示时频域上潜在因素和泉流量时间序列之间的相关关系，对潜在因素的识别进行补充。最后使用有效动态效率定量的评价潜在因素对泉流量的相对贡献。本发明针对复杂的水文过程，开发了一种具有在时频域表征非平稳信号动态结构特征的技术，适用于量化人类活动和气候变化等主要影响因素对泉流量在不同尺度上的影响和贡献。

Claims (8)

1.一种识别量化控制泉流量潜在因素的方法，其特征在于，包括以下步骤：

S1、选定研究的喀斯特地貌区域，获取区域内多组泉流量数据及驱动因素数据，并对获取的数据进行标准化；

S2、异常值检测，判断是否存在异常值，如果存在异常值，则对异常值进行处理；如果不存在异常值，则直接执行步骤S3；

S3、利用动态因素分析提取多组泉流量数据公共的潜在因素，依据在不同动态因素模型中的潜在因素对泉流量的负载大小选择主要的潜在因素；

S4、对时间序列数据进行小波分解；

S5、基于最大分解水平上的近似小波系数对潜在因素进行K-均值聚类分析；

S6、基于小波分解计算潜在因素和驱动因素的多分辨率互相关性，用多分辨率互相关系数绝对值的最大值表征两个时间序列的相似程度和联系；

S7、分析潜在因素和驱动因素之间的小波相干性；

S8、利用对应尺度上的小波能量分布和潜在因素对泉流量负载的乘积表示不同的因素对泉流量的相对贡献。

2.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S1中数据标准化是指减去均值除以标准差对数据进行标准化。

3.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S2中异常值是指时间序列中的值与平均值的偏差超过两倍标准差的值，使用方框图扫描泉流量和驱动因素时间序列是否存在异常值，对异常值进行处理的具体方法为：利用异常值前后对应值的均值代替异常值。

4.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S3中动态因素分析具体为：

F_j(t)＝F_j(t-1)+η_j (2)；

其中，t表示观测时间，m表示考虑的驱动力的数量，Y_i(t)、μ_i(t)、ε_i分别表示第i个观测时间序列的泉流量时间序列、水平参数和具体误差，F_j(t)是第j个共同的潜在因素，α_i,j是第j个潜在因素对第i个泉流量时间序列的负载，η_j是状态方程公式(2)的剩余量；并且假设ε_i和η_j为白噪声，分别用协方差矩阵R和Q来表示；动态因素被R的两种结构抽取：1)方差相等的对角阵；2)方差不相等的对角阵。公式(1)和(2)中的所有参数和动态因素模型通过极大似然法和卡尔曼滤波器并且结合期望最大化算法进行测试；最优动态因素模型(潜在因素的数量m和协方差矩阵R)是根据Akaike信息准则的小样本校正值AICc进行选择；根据不同动态因素模型中的潜在因素对泉流量的负载选择较为重要的潜在因素。

5.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S4具体为：

时间序列和小波函数之间的关系为：

其中，ψ表示母小波，k是时间序列的长度，a和b分别是分解尺度和时间上的平移系数，ψ_a,b(t)表示在分解尺度a，平移因子为b时的母小波；如果a和b以2进制间隔不连续取样，即a＝2^-p，b＝2^-pq，p和q分别表示小波分解尺度和时间上的平移，则公式(3)是离散小波变换对应的母小波，表示为：

ψ_p,q(t)＝2^p/2ψ(2^pt-q) (4)；

其中，ψ_p,q(t)表示在分解尺度p，平移系数为q时的母小波，用来表征高频信号的特征。

表征低频信号域特征的小波函数即父小波表示为：

φ_p,q(t)＝2^p/2φ(2^pt-q) (5)；

则表征信号整个区域的特征的离散小波变换表示为：

其中，f(t)表示目标时间序列；k表示时间序列的长度；J是最大的分解尺度；θ_J,q是在分解尺度J上的近似系数，用于表示低频区域的信息，λ_p,q表示高频区域的小波系数；所有q对应的θ_J,q和λ_p,q分别表示成近似分量A_J和细节分量D_p；并且J根据下面的公式选择：

6.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S6中互相关性系数表示为：

其中，r_xy(τ)表示互相关系数，τ是滞后时间，σ_x和σ_y分别是时间序列x和y的标准差，

和

表示x和y的均值，C_xy(τ)表示x和y的互协方差函数。

7.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S7中分析潜在因素和驱动因素的小波相干性，首先对长度为k的潜在因素x和驱动因素时间序列y进行连续小波变换，为了下面方便表示，并把x和y分别记为x_i,i＝1,2,...,k和y_i,i＝1,2,...,k；x_i和y_i为等间隔(△)取样，则在时间t_j＝j·△和尺度a上的连续小波变换W_i ^x(a)和W_i ^y(a)分别表示为：

其中“*”表示共轭算子；

潜在因素时间序列x_i和驱动因素y_i的小波相干性系数

表示为：

其中，a表示尺度，S表示光滑算子，W_i ^x(a)和W_i ^y(a)表示时间序列x_i和y_i的连续的小波变换，

“*”表示共轭算子。

并且时间序列的相位特征表示为：

Φ_i(a)＝tan^-1(Im(a^-1W_i ^xy(a))/Re(a^-1W_i ^xy(a))) (13)；

其中，Im表示虚部，Re表示实部。

8.根据权利要求1所述的识别量化控制泉流量潜在因素的方法，其特征在于，步骤S8中不同的因素对泉流量的相对贡献率D_ef为：

D_ef＝E_DF*DF_loading (14)；

其中，E_DF和DF_loading分别表示在对应尺度上能量分布和动态因素对时间序列的负载。

Images