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

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

Info

Publication number
CN112036460A
CN112036460A CN202010855119.2A CN202010855119A CN112036460A CN 112036460 A CN112036460 A CN 112036460A CN 202010855119 A CN202010855119 A CN 202010855119A CN 112036460 A CN112036460 A CN 112036460A
Authority
CN
China
Prior art keywords
factors
wavelet
factor
time series
spring
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010855119.2A
Other languages
English (en)
Other versions
CN112036460B (zh
Inventor
谷秀粉
孙洪广
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202010855119.2A priority Critical patent/CN112036460B/zh
Publication of CN112036460A publication Critical patent/CN112036460A/zh
Application granted granted Critical
Publication of CN112036460B publication Critical patent/CN112036460B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/23Clustering techniques
    • G06F18/232Non-hierarchical techniques
    • G06F18/2321Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24Querying
    • G06F16/245Query processing
    • G06F16/2458Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2465Query processing support for facilitating data mining operations in structured databases
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24Querying
    • G06F16/245Query processing
    • G06F16/2458Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2474Sequence data queries, e.g. querying versioned data

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Linguistics (AREA)
  • Software Systems (AREA)
  • Mathematical Physics (AREA)
  • Fuzzy Systems (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
  • Complex Calculations (AREA)

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中动态因素分析具体为:
Figure BDA0002646159540000021
Fj(t)=Fj(t-1)+ηj (2);
其中,t表示观测时间,m表示考虑的驱动力的数量,Yi(t)、μi(t)、εi分别表示第i个观测时间序列的泉流量时间序列、水平参数和具体误差,Fj(t)是第j个共同的潜在因素,αi,j是第j个潜在因素对第i个泉流量时间序列的负载,是量化潜在因素对泉流量影响的重要指标,ηj是状态方程公式(2)的剩余量;并且假设εi和ηj为白噪声,分别用协方差矩阵R和Q来表示;动态因素被R的两种结构抽取:1)方差相等的对角阵;2)方差不相等的对角阵。公式(1)和(2)中的所有参数和动态因素模型通过极大似然法和结合卡尔曼滤波器和期望最大化算法进行测试;最优动态因素分析模型(潜在因素数量m和协方差矩阵R)是根据Akaike信息准则的小样本校正值AICc进行选择;根据不同动态因素模型中的潜在因素对泉流量的负载选择较为重要的潜在因素。
进一步的,步骤S4具体为:
小波分析使用下面的公式探索信号和小波函数之间的关系:
Figure BDA0002646159540000031
其中,ψ是母小波,k是时间序列的长度,a和b分别是分解尺度和平移系数;ψa,b(t)表示在分解尺度a,平移系数为b时的母小波。如果a和b以2进制间隔不连续取样,即a=2-p,b=2-pq,p和q分别表示小波分解尺度和平移系数,则公式(3)是离散小波变换对应的母小波,表示为:
ψp,q(t)=2p/2ψ(2pt-q) (4);
其中,ψp,q(t)表示在分解尺度p,平移系数为q时的母小波,用来表征高频信号的特征;
用下面公式表征低频信号域特征的小波函数即父小波表示为:
φp,q(t)=2p/2φ(2pt-q) (5);
为了表征信号整个区域的特征,则离散小波变换表示成:
Figure BDA0002646159540000032
其中,f(t)表示目标时间序列;k表示时间序列的长度;J是最大的分解尺度;θJ,q是在分解尺度J上近似小波系数,用于表示低频区域的信息,λp,q表示高频区域的小波系数;所有q对应的θJ,q和λp,q分别表示成近似分量AJ和细节分量Dp;并且J根据下面的公式选择:
Figure BDA0002646159540000041
进一步的,步骤S6中互相关性系数表示为:
Figure BDA0002646159540000042
Figure BDA0002646159540000043
其中,rxy(τ)表示互相关系数,τ是滞后时间,σx和σy分别是时间序列x和y的标准差,
Figure BDA0002646159540000044
Figure BDA0002646159540000045
表示x和y的均值,Cxy(τ)表示x和y的互协方差函数。
利用潜在因素的最大分解尺度上的近似小波系数进行聚类,不仅保留了潜在因素的波动特征,而且通过降低维数提高聚类的效率。
进一步的,步骤S7中首先对长度为k的潜在因素时间序列x和驱动因素时间序列y进行连续小波变换;为了下面方便表示,将x和y分别记为xi,i=1,2,...,k和yi,i=1,2,...,k;xi和yi为等间隔(△)取样,则在时间tj=j·△和尺度a上的连续小波变换Wi x(a)和Wi y(a)分别表示为:
Figure BDA0002646159540000046
Figure BDA0002646159540000047
其中“*”表示共轭算子。
潜在因素时间序列xi和驱动因素序列yi的小波相干性系数
Figure BDA0002646159540000048
表示为:
Figure BDA0002646159540000049
其中,a表示尺度,S表示光滑算子,Wi x(a)和Wi y(a)表示时间序列x和y的连续的小波变换,
Figure BDA0002646159540000051
“*”表示共轭算子。
并且时间序列的相位特征表示为:
Φi(a)=tan-1(Im(a-1Wi xy(a))/Re(a-1Wi xy(a))) (13);
其中,Im表示虚部,Re表示实部。
进一步的,步骤S8中不同的因素对泉流量的相对贡献率Def为:
Def=EDF*DFloading (14);
其中,EDF和DFloading分别表示在对应尺度上能量分布和动态因素对时间序列的负载。
发明原理:动态因素分析可以有效揭示时间序列的动态结构,识别潜在因素。小波分析是分解和量化时间序列的理想工具,本发明是结合动态因素分析和小波分析的优点发展了一种识别和量化潜在因素控制泉流量的技术。
有益效果:与现有技术相比,本发明利用动态因素分析揭示时间序列的动态结构特征,抽取泉流量时间序列的公共潜在因素,利用多分辨率互相关性和小波相干性分析有效识别潜在因素及其对泉流量影响的作用尺度。最后利用有效动态效率量化主要潜在因素对泉流量的相对贡献。结合动态因素分析和小波分析有效的解决现有方法无法定量评价潜在因素对泉流量影响的缺陷。第二,为研究复杂的水文过程提供了方便可靠的方法,具有重要的理论和工程意义。
附图说明
图1为本发明所述方法的流程图;
图2a为在组1中潜在因素对不同泉流量时间序列的负载图,2b为在组2中潜在因素对不同泉流量时间序列的负载图,2c为在组3中潜在因素对不同泉流量时间序列的负载图;
图3为潜在因素与潜在因素的最高负载对应的泉流量波动特征;
图4为潜在因素的小波聚类图;
图5a为组1中潜在因素与潜在因素最高负载对应的泉流量序列的小波相干性图;5b为组2中潜在因素与潜在因素最高负载对应的泉流量序列的小波相干性图;5c为组3中潜在因素与潜在因素最高负载对应的泉流量序列的小波相干性图。
具体实施方式
下面结合说明书附图,对本发明作进一步的说明。
如图1所示,本发明提出的识别量化控制泉流量潜在因素的方法,包括如下步骤:
S1、选定所研究的喀斯特地貌区域,获取区域内多组泉流量及驱动因素数据,并对获取的数据进行标准化。
本实施例选定的要研究的喀斯特区域为娘子关泉流域,在娘子关泉流域水文过程的主要驱动因素包括降水,气候变化和人类活动;数据标准化是指减去均值除以标准差。
S2、异常值检测,使用方框图扫描泉流量和驱动因素时间序列数据,判断是否存在异常值,如果存在异常值,则用异常值前后对应值的均值代替;如果不存在异常值,则直接执行步骤S3。
所述异常值是指时间序列中的值与平均值的偏差超过两倍标准差的值。
S3、利用动态因素分析提取多组泉流量数据公共的潜在因素,依据在不同动态因素模型中的潜在因素对泉流量的负载大小选择主要的潜在因素;
动态因素分析可以揭示时间序列的动态结构,该方法对具有记忆特征的非平稳时间序列仍然有效。
动态因素分析具体为:
Figure BDA0002646159540000061
Fj(t)=Fj(t-1)+ηj (2);
其中,t表示观测时间,m表示考虑的驱动力的数量,Yi(t)、μi(t)、εi分别表示第i个观测的泉流量时间序列、水平参数和具体误差,Fj(t)是第j个共同的潜在因素,αi,j是第j个潜在因素对第i个泉流量序列的负载,是量化潜在因素对泉流量影响的重要指标。潜在ηj是状态方程公式(2)的剩余量;并且假设εi和ηj为白噪声,分别用协方差矩阵R和Q来表示。动态因素被R的两种结构抽取:1)方差相等的对角阵;2)方差不相等的对角阵。公式(1)和(2)中的所有参数和动态因素模型通过极大似然法和卡尔曼滤波器并且结合期望最大化算法进行测试。最优动态因素模型(潜在因素的数量(m)和协方差矩阵(R))是根据Akaike信息准则的小样本校正值(AICc)进行选择。根据不同动态因素模型的潜在因素对泉流量的负载选择较为重要的潜在因素。
S4、对潜在因素,泉流量和驱动因素时间序列数据进行小波分解。小波分解可以将时间序列分解到不同的尺度水平上,揭示不同尺度下的波动特征;并且潜在因素的小波能量分布表征了在不同尺度上泉流量受驱动因素的影响程度,是表征潜在因素对泉流量相对贡献的重要参考指标。
小波分析分解泉流量,潜在因素和驱动因素时间序列,使用下面的方程探索信号和小波函数之间的关系:
Figure BDA0002646159540000071
其中,ψ是母小波,k是时间序列的长度,a表示分解尺度,b是平移系数,反映时间上的平移,ψa,b(t)表示在分解尺度a,平移系数为b时的母小波;如果这两个系数是连续的,公式(3)是连续小波变换对应的母小波;如果a和b以2进制间隔不连续取样(a=2-p,b=2-pq,p和q分别表示小波分解尺度和时间上的平移系数),则公式(3)是离散小波变换,其对应的母小波表示为:
ψp,q(t)=2p/2ψ(2pt-q) (4);
其中,ψp,q(t)为在分解尺度p和平移因子q下的母小波,主要表征高频信号的特征。用下面公式表征低频信号域特征的小波函数即父小波表示为:
φp,q(t)=2p/2φ(2pt-q) (5);
φ表示父小波,φp,q(t)表示在分解尺度为p,平移因子为q时的父小波。
为了表征泉流量,潜在因素和驱动因素时间序列在整个区域的特征,则离散小波变换可以表示为:
Figure BDA0002646159540000072
其中,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分别表示成近似分量AJ和细节分量Dp。并且J根据下面的公式选择:
Figure BDA0002646159540000081
其中f1表示母小波滤波器长度。
S5、基于最大尺度的近似小波系数对潜在因素进行K-均值聚类,潜在因素的最大尺度的小波系数不仅可以表征潜在因素的波动特征,而且通过减少潜在因素时间序列的维数来提高聚类效率。
S6、基于小波分解计算潜在因素和驱动因素的多分辨率互相关系数,并获得不同分辨率下的互相关系数绝对值的最大值,用来表征潜在因素和驱动因素时间序列之间的相似程度。互相关性系数表示为:
Figure BDA0002646159540000082
Figure BDA0002646159540000083
其中,rxy(τ)表示时间序列x和y互相关系数,τ是滞后时间,σx和σy分别是潜在因素时间序列x和驱动因素时间序列y的标准差,
Figure BDA0002646159540000084
Figure BDA0002646159540000085
分别表示x和y的均值,xt和yt分别表示x和y在t时刻的值,Cxy(τ)表示x和y的互协方差函数。
用多分辨率互相关系数绝对值的最大值表征两个时间序列的相似程度和联系,这里我们主要考虑的是两个时间序列的相关性的大小,忽略了关系的正负性。
S7、分析潜在因素和驱动因素的小波相干性,揭示他们时频域上的关系特征,目的是为了对多分辨率互相关方法识别潜在因素进行补充。并且小波相干性分析可以将时间序列之间的延迟反应和关系特征进行区分表征。在引入小波相干性之前,首先介绍连续的小波变换。对长度为k的潜在因素时间序列x和驱动因素时间序列y进行连续小波变换。为了下面公式的方便表达,把x和y进一步分别记为xi,i=1,2,...,k和yi,i=1,2,...,k。xi和yi为等间隔(记为△)取样,则在时间Tj=j·△和尺度a上的连续小波变换Wi x(a),Wi y(a)分别表示为:
Figure BDA0002646159540000091
Figure BDA0002646159540000092
其中“*”表示共轭算子。
潜在因素时间序列xi和驱动因素时间序列yi的小波相干性系数
Figure BDA0002646159540000093
表示为:
Figure BDA0002646159540000094
其中,a表示尺度,S表示光滑算子,Wi x(a)和Wi y(a)分别表示时间序列x和y的连续小波变换,
Figure BDA0002646159540000095
“*”表示共轭算子。
并且潜在因素和驱动因素时间序列的相对相位表示为:
Φi(a)=tan-1(Im(a-1Wi xy(a))/Re(a-1Wi xy(a))) (13);
其中,Im表示虚部,Re表示实部。
S8、根据小波能量分布和负载特征表示不同的潜在因素对泉流量的相对贡献。
Def=EDF*DFloading (14);
其中,EDF和DFloading分别表示在对应尺度上能量分布和动态因素对时间序列的负载,Def为有效动态效率。
实施例:
该实施例的数据来自中国北方最大的喀斯特泉—娘子关泉。娘子关泉群位于阳泉市平定县娘子关镇,由几十个单泉组成。娘子关泉接收大约7394平方千米集水区的补充。由于娘子关泉流域是北方地区的重工业区,所以在这个实施例中,除了考虑对泉流量有重要影响的气候因素外,还考虑了人类活动的影响。因此考虑的主要驱动因素具体包括季风(ISM:印度夏季风和WNPM:西太平洋北部季风),降水,地下水补充和人类活动。利用滑动窗对获得的1959-2007年泉流量数据进行分组,一共选取了3组泉流量数据进行具体的分析,如表1所示。
表1.泉流量时间序列的详细划分方式
Figure BDA0002646159540000096
Figure BDA0002646159540000101
表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中不同协方差矩阵下的动态因素分析结果
Figure BDA0002646159540000102
表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.潜在因素和驱动力因素的小波能量分布
Figure BDA0002646159540000111
Figure BDA0002646159540000121
表4展示了潜在因素和泉流量时间序列在不同分辨率水平下的互相关系数绝对值的最大值,组1中的DF41(4-16个月),组3中的DF41(2-16个月)和DF44(2-16个月)与泉流量显示出最大的互相干系数,考虑到地质介质对局部地下水的补给具有相对较小的阻尼作用,并且根据图2的结果,说明这三个因素主要表示局部地下水补给对泉流量的影响。可以根据类似的方法,对其他因素进行识别。
表4.潜在因素和泉流量时间序列在多分辨率水平上互相关系数绝对值的最大值
Figure BDA0002646159540000122
表4中主要潜在因素加粗倾斜表示,最大的相关性系数加粗强调。
图5a-5c表示组1中主要潜在因素(DF41和DF44)与最大潜在因素负载所对应的泉流量时间序列之间的小波相干性结果。对于DF41(代表局部地下水的补充),与泉流量之间显示出很强的相关性。并且表征了与泉流量的同相特征。DF44在2-4个月和≥16个月的周期上表示诸如煤矿开采,发电和化学工程等人类活动对泉流量时间序列的影响。同样的,对组2和组3进行分析。通过多分辨率互相关系数和小波相干性分析,可以对主要的潜在因素进行有效识别。并且通过有效动态效率表示主要潜在因素对泉流量的相对贡献。则主要的结论如下所示:
第一组:DF41表示地下水开发前的局部地下水的补给(4-16个月),|Rmax|=1.00,Def≤42.61%;DF44表示人类活动的影响(2-4个月),|Rmax|=0.67-0.89,Def≤0.36%和(≥16个月),|Rmax|=0.91-0.92,Def≤27.16%。
第二组:DF41表示降水和径流的影响(4-16个月),|Rmax|=0.88-0.99,Def≤29.87%;DF42表示区域地下水补充的影响(4-16个月),|Rmax|=0.93-1.00,Def≤29.90%。
第三组:DF43代表季风的影响(8-16个月),|Rmax|=0.74-0.78,Def≤3.02%;DF44表示地下水开采后地下水的补充的影响(2-16个月),|Rmax|=1.00,Def≤41.00%。
因此,本发明提出的识别和定量评价控制泉流量潜在因素的分析方法,利用动态因素分析提取泉流量公共的潜在因素。多分辨率互相关系数用来揭示主要潜在因素与驱动力之间的相似程度,从而对潜在因素进行识别。小波相干性分析用来揭示时频域上潜在因素和泉流量时间序列之间的相关关系,对潜在因素的识别进行补充。最后使用有效动态效率定量的评价潜在因素对泉流量的相对贡献。本发明针对复杂的水文过程,开发了一种具有在时频域表征非平稳信号动态结构特征的技术,适用于量化人类活动和气候变化等主要影响因素对泉流量在不同尺度上的影响和贡献。

Claims (8)

1.一种识别量化控制泉流量潜在因素的方法,其特征在于,包括以下步骤:
S1、选定研究的喀斯特地貌区域,获取区域内多组泉流量数据及驱动因素数据,并对获取的数据进行标准化;
S2、异常值检测,判断是否存在异常值,如果存在异常值,则对异常值进行处理;如果不存在异常值,则直接执行步骤S3;
S3、利用动态因素分析提取多组泉流量数据公共的潜在因素,依据在不同动态因素模型中的潜在因素对泉流量的负载大小选择主要的潜在因素;
S4、对时间序列数据进行小波分解;
S5、基于最大分解水平上的近似小波系数对潜在因素进行K-均值聚类分析;
S6、基于小波分解计算潜在因素和驱动因素的多分辨率互相关性,用多分辨率互相关系数绝对值的最大值表征两个时间序列的相似程度和联系;
S7、分析潜在因素和驱动因素之间的小波相干性;
S8、利用对应尺度上的小波能量分布和潜在因素对泉流量负载的乘积表示不同的因素对泉流量的相对贡献。
2.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S1中数据标准化是指减去均值除以标准差对数据进行标准化。
3.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S2中异常值是指时间序列中的值与平均值的偏差超过两倍标准差的值,使用方框图扫描泉流量和驱动因素时间序列是否存在异常值,对异常值进行处理的具体方法为:利用异常值前后对应值的均值代替异常值。
4.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S3中动态因素分析具体为:
Figure FDA0002646159530000011
Fj(t)=Fj(t-1)+ηj (2);
其中,t表示观测时间,m表示考虑的驱动力的数量,Yi(t)、μi(t)、εi分别表示第i个观测时间序列的泉流量时间序列、水平参数和具体误差,Fj(t)是第j个共同的潜在因素,αi,j是第j个潜在因素对第i个泉流量时间序列的负载,ηj是状态方程公式(2)的剩余量;并且假设εi和ηj为白噪声,分别用协方差矩阵R和Q来表示;动态因素被R的两种结构抽取:1)方差相等的对角阵;2)方差不相等的对角阵。公式(1)和(2)中的所有参数和动态因素模型通过极大似然法和卡尔曼滤波器并且结合期望最大化算法进行测试;最优动态因素模型(潜在因素的数量m和协方差矩阵R)是根据Akaike信息准则的小样本校正值AICc进行选择;根据不同动态因素模型中的潜在因素对泉流量的负载选择较为重要的潜在因素。
5.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S4具体为:
时间序列和小波函数之间的关系为:
Figure FDA0002646159530000021
其中,ψ表示母小波,k是时间序列的长度,a和b分别是分解尺度和时间上的平移系数,ψa,b(t)表示在分解尺度a,平移因子为b时的母小波;如果a和b以2进制间隔不连续取样,即a=2-p,b=2-pq,p和q分别表示小波分解尺度和时间上的平移,则公式(3)是离散小波变换对应的母小波,表示为:
ψp,q(t)=2p/2ψ(2pt-q) (4);
其中,ψp,q(t)表示在分解尺度p,平移系数为q时的母小波,用来表征高频信号的特征。
表征低频信号域特征的小波函数即父小波表示为:
φp,q(t)=2p/2φ(2pt-q) (5);
则表征信号整个区域的特征的离散小波变换表示为:
Figure FDA0002646159530000022
其中,f(t)表示目标时间序列;k表示时间序列的长度;J是最大的分解尺度;θJ,q是在分解尺度J上的近似系数,用于表示低频区域的信息,λp,q表示高频区域的小波系数;所有q对应的θJ,q和λp,q分别表示成近似分量AJ和细节分量Dp;并且J根据下面的公式选择:
Figure FDA0002646159530000031
6.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S6中互相关性系数表示为:
Figure FDA0002646159530000032
Figure FDA0002646159530000033
其中,rxy(τ)表示互相关系数,τ是滞后时间,σx和σy分别是时间序列x和y的标准差,
Figure FDA0002646159530000034
Figure FDA0002646159530000035
表示x和y的均值,Cxy(τ)表示x和y的互协方差函数。
7.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S7中分析潜在因素和驱动因素的小波相干性,首先对长度为k的潜在因素x和驱动因素时间序列y进行连续小波变换,为了下面方便表示,并把x和y分别记为xi,i=1,2,...,k和yi,i=1,2,...,k;xi和yi为等间隔(△)取样,则在时间tj=j·△和尺度a上的连续小波变换Wi x(a)和Wi y(a)分别表示为:
Figure FDA0002646159530000036
Figure FDA0002646159530000037
其中“*”表示共轭算子;
潜在因素时间序列xi和驱动因素yi的小波相干性系数
Figure FDA0002646159530000038
表示为:
Figure FDA0002646159530000039
其中,a表示尺度,S表示光滑算子,Wi x(a)和Wi y(a)表示时间序列xi和yi的连续的小波变换,
Figure FDA00026461595300000310
“*”表示共轭算子。
并且时间序列的相位特征表示为:
Φi(a)=tan-1(Im(a-1Wi xy(a))/Re(a-1Wi xy(a))) (13);
其中,Im表示虚部,Re表示实部。
8.根据权利要求1所述的识别量化控制泉流量潜在因素的方法,其特征在于,步骤S8中不同的因素对泉流量的相对贡献率Def为:
Def=EDF*DFloading (14);
其中,EDF和DFloading分别表示在对应尺度上能量分布和动态因素对时间序列的负载。
CN202010855119.2A 2020-08-24 2020-08-24 一种识别量化控制泉流量潜在因素的方法 Active CN112036460B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010855119.2A CN112036460B (zh) 2020-08-24 2020-08-24 一种识别量化控制泉流量潜在因素的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010855119.2A CN112036460B (zh) 2020-08-24 2020-08-24 一种识别量化控制泉流量潜在因素的方法

Publications (2)

Publication Number Publication Date
CN112036460A true CN112036460A (zh) 2020-12-04
CN112036460B CN112036460B (zh) 2022-08-30

Family

ID=73580621

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010855119.2A Active CN112036460B (zh) 2020-08-24 2020-08-24 一种识别量化控制泉流量潜在因素的方法

Country Status (1)

Country Link
CN (1) CN112036460B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112988815A (zh) * 2021-03-16 2021-06-18 重庆工商大学 一种大规模高维高速流数据在线异常检测的方法及系统
CN114440758A (zh) * 2022-01-09 2022-05-06 西北大学 一种区域尺度上滑坡对降雨响应的分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5924048A (en) * 1997-03-14 1999-07-13 Mccormack; Michael D. Automated material balance system for hydrocarbon reservoirs using a genetic procedure
CN110120046A (zh) * 2019-03-27 2019-08-13 长安大学 一种融合dem、光学遥感和形变信息的潜在滑坡识别方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5924048A (en) * 1997-03-14 1999-07-13 Mccormack; Michael D. Automated material balance system for hydrocarbon reservoirs using a genetic procedure
CN110120046A (zh) * 2019-03-27 2019-08-13 长安大学 一种融合dem、光学遥感和形变信息的潜在滑坡识别方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112988815A (zh) * 2021-03-16 2021-06-18 重庆工商大学 一种大规模高维高速流数据在线异常检测的方法及系统
CN112988815B (zh) * 2021-03-16 2023-09-05 重庆工商大学 一种大规模高维高速流数据在线异常检测的方法及系统
CN114440758A (zh) * 2022-01-09 2022-05-06 西北大学 一种区域尺度上滑坡对降雨响应的分析方法
CN114440758B (zh) * 2022-01-09 2023-04-14 西北大学 一种区域尺度上滑坡对降雨响应的分析方法

Also Published As

Publication number Publication date
CN112036460B (zh) 2022-08-30

Similar Documents

Publication Publication Date Title
Liu et al. Modeling the daily suspended sediment concentration in a hyperconcentrated river on the Loess Plateau, China, using the Wavelet–ANN approach
Saad et al. Scalodeep: A highly generalized deep learning framework for real‐time earthquake detection
Nalley et al. Using discrete wavelet transforms to analyze trends in streamflow and precipitation in Quebec and Ontario (1954–2008)
Li et al. Trend modeling for traffic time series analysis: An integrated study
Sang A practical guide to discrete wavelet decomposition of hydrologic time series
CN110361778B (zh) 一种基于生成对抗网络的地震数据重建方法
Zamani et al. Learning from data for wind–wave forecasting
CN103728551B (zh) 一种基于级联集成分类器的模拟电路故障诊断方法
Xu et al. Wavelet-denoising multiple echo state networks for multivariate time series prediction
CN112036460B (zh) 一种识别量化控制泉流量潜在因素的方法
Sugiartawan et al. Prediction by a hybrid of wavelet transform and long-short-term-memory neural network
CN103310789B (zh) 一种基于改进的并行模型组合的声音事件识别方法
Cai et al. Impacts of regional characteristics on improving the accuracy of groundwater level prediction using machine learning: The case of central eastern continental United States
CN101587186B (zh) 一种雷达脉内调制信号的特征提取方法
Avci et al. A new automatic target recognition system based on wavelet extreme learning machine
CN104200291A (zh) 一种基于小波变换和arma-svm的涌水量预测方法
CN108549908A (zh) 基于多采样概率核主成分模型的化工过程故障检测方法
Oh et al. The combined use of dynamic factor analysis and wavelet analysis to evaluate latent factors controlling complex groundwater level fluctuations in a riverside alluvial aquifer
CN113780242A (zh) 一种基于模型迁移学习的跨场景水声目标分类方法
Zaifoglu et al. Regional frequency analysis of precipitation using time series clustering approaches
CN115758082A (zh) 一种轨道交通变压器故障诊断方法
Tran et al. Data reformation–A novel data processing technique enhancing machine learning applicability for predicting streamflow extremes
CN117612025A (zh) 基于扩散模型的遥感图像屋顶识别方法及系统
Nicolet et al. A multi-criteria leave-two-out cross-validation procedure for max-stable process selection
CN104200472A (zh) 基于非局部小波信息的遥感图像变化检测方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant