基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法技术

本发明专利技术为一种基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法，选取N个钻孔应变台站作为N个网络节点；设置一个时间窗口长度n和移动步长step，构造时间窗口内的钻孔应变大小为n×N的多维节点矩阵；多维节点矩阵进行多通道奇异谱分析，分解后得到特征向量和对应空间时间主成分；选取高频成分重构各网络节点的有效成分；计算节点两两之间的皮尔森相关系数；连通相关系数大于0.8的两网络节点，构造应变网络的边，写出应变网络的邻接矩阵；计算应变网络的拓扑结构属性；并移动时间窗口，构造时间窗口内的应变网络并计算拓扑结构属性；输出钻孔应变网络结构的演化图，判断地震相关异常。本发明专利技术实现了钻孔应变数据的多台站联合分析。

【技术实现步骤摘要】
基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法
本专利技术属于数据为钻孔应变观测数据处理领域，具体地而言为一种基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法。
技术介绍
受到美国“板块边界观测(PBO)”的启发，中国改造和建立了约100个钻孔应变观测点，其中约40套YRY-4型四元件钻孔应变仪和约60套TJ型钻孔体应变仪。YRY-4型钻孔应变仪的安装深度为40m左右，而TJ型钻孔体应变仪的深度约为60m。建立的钻孔应变观测点采样间隔为1分钟，经过多年运行，大部分观测点观测仪器工作状态稳定，可以清晰的记录到光滑的固体潮。特别是，观测到的钻孔应变数据在满足自检关系方面基本达到理想要求，为获取高质量的观测数据提供了坚实的基础。钻孔应变观测在人类对地壳活动研究中扮演着重要的角色。随着钻孔应变观测技术的日益成熟，研究人员将钻孔应变数据应用到地震应变阶、慢地震、火山喷发预测、地球自由震荡和地震前兆异常提取等研究领域，并取得了重要进展。目前，通过对钻孔应变数据震前异常的分析，我国学者针对唐山地震、汶川地震和芦山地震的研究发现了地震前兆的证据。但是绝大多数的研究是对单独台站的数据进行独立分析，这就大大增加了判断地震前兆异常的风险。而多个观测台站的联合分析可以消除这种局限，增加前兆异常判断的可靠性。在地震学中，地震活动具有复杂的时空分布特征，其孕育环境则是地壳内复杂的力学系统。复杂网络就是大量真实复杂系统的抽象，它能够刻画复杂系统内部的各种相互作用或关系。从统计物理学来看，网络是一个包含了大量个体及个体之间相互作用的系统。典型的网络是由许多节点与连接节点间的一些边组成的，其中节点用来代表真实系统中不同的个体，而边则用来表示个体间的关系。目前，对于孤立台站的钻孔应变观测数据，研究人员已经应用合适的方法提取与地震活动有关的信息。CN106918836A公开了一种基于主成分分析的钻孔应变数据异常提取方法，有效的利用主成分分析的方法对钻孔应变数据进行分析，利用主成分分析中的特征值和特征向量角度分别将地壳的微弱变化表征出来；实现了在有较强背景干扰的情况下对钻孔应变数据异常的精确提取。CN109031403B公开了一种基于S-K特征的钻孔应变数据异常提取方法，计算每一天的面应变差分数据的偏度和峰度；定义了一个地壳平稳时的背景并计算每一天与背景的偏移程度，通过该专利技术能够有效的对钻孔应变数据进行分析，对微小的地震前兆异常进行提取。然而，对多观测站点的联合分析，不仅需要提取有效的地壳形变信号，还要考虑多观测站点间的共同特征，因此并不适用于对多观测站点的联合分析。
技术实现思路
本专利技术所要解决的技术问题在于提供一种基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法，对多观测站点对钻孔应变数据进行分析。本专利技术是这样实现的，一种基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法，包括以下步骤：步骤a，选取N个钻孔应变台站作为N个网络节点；步骤b，设置一个时间窗口长度n和移动步长step，构造时间窗口内的钻孔应变大小为n×N的多维节点矩阵；步骤c，多维节点矩阵进行多通道奇异谱分析，分解后得到特征向量和对应空间时间主成分；步骤d，选取高频成分重构各网络节点的有效成分；步骤e，对N个网络节点之间的有效成分，计算节点两两之间的皮尔森相关系数；步骤f，连通相关系数大于0.8的两网络节点，构造应变网络的边，写出应变网络的邻接矩阵；步骤g，计算应变网络的拓扑结构属性；步骤h，移动时间窗口，重复步骤b-g，构造时间窗口内的应变网络并计算拓扑结构属性；步骤i，输出钻孔应变网络结构的演化图，判断地震相关异常。进一步地，步骤b，多维节点矩阵为其中，xij是长度为n的窗口内，N个台站的钻孔应变观测序列，i＝1,2,...,N,j＝1,2,...,n。进一步地，步骤c，多维节点矩阵进行多通道奇异谱分析，包括对多维节点矩阵做如下时迟排列：其中，矩阵元素xlt,(t＝1,2,...,n；l＝1,2,...,N)其下标t,l分别为样本时序和地理空间站点变量序号，m为最大时迟数；对Y矩阵作时空分解，其分量形式为特征向量Ek是N×m维，分量是第k个特征向量在l通道滞后j的分量，是N维向量的m个不同时间滞后的序列，特征向量Ek称为空间-时间域的EOF，时间系数是第i个状态Xi在Ek上的投影，序列称为第k个空间-时间主成分，其中1≤i≤N-M+1。进一步地，步骤d，选取选取高频成分重构各网络节点的有效成分包括：当仅使用第k个特征成分来重建振荡分量序列时，公式为上式中，表示场中第l个变量在第i个时序的第k个重建分量，根据钻孔应变数据特征，选择高频的特征因子的子集S来重构各节点的有效成分，进一步地，步骤e计算皮尔森相关系数包括：设步骤d得到的l个节点重构成分为计算节点两两之间的皮尔森相关系数rij,i,j∈l,公式中，Xip是第i个节点的第p个元素，Xjp是第j个节点的第p个元素，i,j∈l,p＝1,..,n。进一步地，步骤f具体包括：判断相关系数，假设相关系数大于0.8，则两节点强相关，连通两节点构造拓扑网络的边；如果相关系数小于0.8，则认为两节点弱相关，不能连通，节点间没有边；写出应变网络的邻接矩阵A，设A＝{aij},i,j∈{1,...,l}，若节点i和节点j之间，有边连通时，则aij＝1；反之，两节点间没有边连通时，aij＝0，其中aij是应变网络的邻接矩阵中的对应元素。进一步地，步骤g具体包括：利用三个静态的应变网络拓扑结构属性来评估钻孔应变网络的特征，分别是平均度聚类系数C和平均路径L；N个节点构成的网络平均节点度为各节点i的度ki的和，度ki为与该节点连接的边的总和：认为平均度表示钻孔应变各台站间的联系强度，其中aij是应变网络的邻接矩阵中的对应元素；聚类系数C用于描述网络中节点的聚类情况，各节点i的ki个邻居节点之间的实际边数Ei与总的可能边数的比就是节点i的聚类系数ci，那么网络的平均聚类系数C为：平均路径长度L为网络中任意2个节点之间距离dij的平均值：进一步地，步骤h具体包括：将长度为n的时间窗口向终止时刻方向移动，步长为step，每移动一步，重复步骤b～g，即取多通道奇异谱分析后，重构各节点的合适成分，根据节点间相关性，构造应变网络，并计算平均度聚类系数C和平均路径L。这样，当窗口从钻孔应变序列起始时刻连续移动到终止时刻，就获得了一系列随时间演化的平均度聚类系数C和平均路径L值。本专利技术与现有技术相比，有益效果在于：本专利技术将多站点钻孔应变观测网络化，提取观测网络中的地震相关异常。首先是将每个钻孔应变台站定义为应变网络的一个节点，对多个节点台站的钻孔应变序列进行多通道奇异谱分析；选取分解后的部分有效的特征向量与对应主成分重构；计算各节点之间本文档来自技高网...

【技术保护点】
1.一种基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法，其特征在于，包括以下步骤：/n步骤a，选取N个钻孔应变台站作为N个网络节点；/n步骤b，设置一个时间窗口长度n和移动步长step，构造时间窗口内的钻孔应变大小为n×N的多维节点矩阵；/n步骤c，多维节点矩阵进行多通道奇异谱分析，分解后得到特征向量和对应空间时间主成分；/n步骤d，选取高频成分重构各网络节点的有效成分；/n步骤e，对N个网络节点之间的有效成分，计算节点两两之间的皮尔森相关系数；/n步骤f，连通相关系数大于0.8的两网络节点，构造应变网络的边，写出应变网络的邻接矩阵；/n步骤g，计算应变网络的拓扑结构属性；/n步骤h，移动时间窗口，重复步骤b-g，构造时间窗口内的应变网络并计算拓扑结构属性；/n步骤i，输出钻孔应变网络结构的演化图，判断地震相关异常。/n

【技术特征摘要】
1.一种基于多通道奇异谱的钻孔应变网络拓扑结构异常检测方法，其特征在于，包括以下步骤：
步骤a，选取N个钻孔应变台站作为N个网络节点；
步骤b，设置一个时间窗口长度n和移动步长step，构造时间窗口内的钻孔应变大小为n×N的多维节点矩阵；
步骤c，多维节点矩阵进行多通道奇异谱分析，分解后得到特征向量和对应空间时间主成分；
步骤d，选取高频成分重构各网络节点的有效成分；
步骤e，对N个网络节点之间的有效成分，计算节点两两之间的皮尔森相关系数；
步骤f，连通相关系数大于0.8的两网络节点，构造应变网络的边，写出应变网络的邻接矩阵；
步骤g，计算应变网络的拓扑结构属性；
步骤h，移动时间窗口，重复步骤b-g，构造时间窗口内的应变网络并计算拓扑结构属性；
步骤i，输出钻孔应变网络结构的演化图，判断地震相关异常。


2.按照权利要求1所述的方法，其特征在于，
步骤b，多维节点矩阵为



其中，xij是长度为n的窗口内，N个台站的钻孔应变观测序列，i＝1,2,...,N,j＝1,2,...,n。


3.按照权利要求1所述的方法，其特征在于，
步骤c，多维节点矩阵进行多通道奇异谱分析，包括对多维节点矩阵做如下时迟排列：



其中，矩阵元素xlt,(t＝1,2,...,n；l＝1,2,...,N)其下标t,l分别为样本时序和地理空间站点变量序号，m为最大时迟数；
对Y矩阵作时空分解，其分量形式为






特征向量Ek是N×m维，分量是第k个特征向量在l通道滞后j的分量，是N维向量的m个不同时间滞后的序列，特征向量Ek称为空间-时间域的EOF，时间系数是第i个状态Xi在Ek上的投影，序列称为第k个空间-时间主成分，其中1≤i≤N-M+1。


4.按照权利要求3所述的方法，其特征在于，
步骤d，选取选取高频成分重构各网络节点的有效成分包括：当仅使用第k个特征成分来重建振荡分量序列时，公式为



上式中，表示场中第l个变量在第i个时序的第k个重建分量，根...

【专利技术属性】
技术研发人员：朱凯光，于紫凝，池成全，樊蒙璇，贺小丹，孙慧慧，
申请(专利权)人：吉林大学，
类型：发明
国别省市：吉林;22