CN110908001A - 一种大地电磁测深信号的重构方法及系统 - Google Patents

一种大地电磁测深信号的重构方法及系统 Download PDF

Info

Publication number
CN110908001A
CN110908001A CN201911292329.9A CN201911292329A CN110908001A CN 110908001 A CN110908001 A CN 110908001A CN 201911292329 A CN201911292329 A CN 201911292329A CN 110908001 A CN110908001 A CN 110908001A
Authority
CN
China
Prior art keywords
modal component
order
signal
modal
component sequence
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
CN201911292329.9A
Other languages
English (en)
Other versions
CN110908001B (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.)
Jilin University
Original Assignee
Jilin University
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 Jilin University filed Critical Jilin University
Priority to CN201911292329.9A priority Critical patent/CN110908001B/zh
Publication of CN110908001A publication Critical patent/CN110908001A/zh
Application granted granted Critical
Publication of CN110908001B publication Critical patent/CN110908001B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38Processing data, e.g. for analysis, for interpretation, for correction

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Investigating Or Analyzing Materials By The Use Of Magnetic Means (AREA)

Abstract

本发明涉及一种大地电磁测深信号的重构方法及系统。该方法包括获取初始的大地电磁测深信号;对初始的大地电磁测深信号进行互补集合经验模态分解,得到多个模态分量序列;对多个模态分量序列进行希尔伯特变换,得到每个模态分量序列的瞬时频谱;根据瞬时频谱的频点的分布情况,筛选瞬时频谱对应的模态分量序列;对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。本发明所提供得一种大地电磁测深信号的重构方法及系统,解决现有技术中对于“死频带”的大地电磁测深信号不能实现良好的去噪效果的问题。

Description

一种大地电磁测深信号的重构方法及系统
技术领域
本发明涉及信号处理领域,特别是涉及一种大地电磁测深信号的重构方法及系统。
背景技术
大地电磁测深法(Magnetotelluric sounding,简称MT)使用天然交变电磁场作为场源。其中,0.1Hz-10Hz频段的天然大地电磁场信号较弱,该频段测得的大地电磁测深数据易受到各类电磁噪声的干扰,因此称该频段为大地电磁测深法的“死频段”。随着深地资源勘查开采项目的推进,要求大地电磁测深法具有更高的数据质量,但矿集区的大地电磁测深测量越来越多的受到矿山采矿、电网输电、信号发射塔等的影响,往往使得原本就较弱的“死频段”信号严重被干扰,强噪声干扰甚至会淹没有效信号,影响了大地电磁测深(即MT)的应用效果。因此,对“死频段”电磁噪声进行有效压制,提高干扰区的MT数据质量,成为MT方法应用的重要挑战。
目前用于大地电磁去噪的方法主要包括以下几种:①高通滤波、低通滤波、带通滤波、简单功率谱、远参考、robust阻抗估计等传统方法;②Rhoplus校正、数学形态滤波、小波变换、Hilbert-Huang变换(简称HHT)、S变换、压缩感知重构等新兴数字滤波方法。
传统方法往往受限于待处理数据的信噪比,以及噪声的相关性等,当处理信噪比较低或含有相关噪声的大地电磁数据时难以取得令人满意的滤波效果。在新兴的方法中,Rhoplus多应用于音频大地电磁测深(Audiomagnetotelluric sounding,简称AMT)。而在“死频带”噪声压制中,AMT与MT“死频带”范围不一致,且Rhoplus受到观测时间与地区的影响。即Rhoplus需要大量的观测数据作为支撑;但是,数据观测需要大量人力、物力,造成经济成本较高,不适用于大地电磁测深“死频带”噪声压制。小波去噪、S变换、数学形态滤波、压缩感知重构等其他数字滤波方法受到各自适用噪声类型,分解正交基的选择,信噪比等限制,在噪声类型复杂,信噪比低的大地电磁测深“死频段”往往不能取得通用性(适用于各类噪声压制)的滤波效果;HHT法具有很大的优越性,它通过EMD(Empirical ModeDecomposition)自适应地提取非平稳的大地电磁信号中的局部均值信息,然后通过Hilbert变换得到MT信号瞬时频谱信息,不存在正交基选择的困难,并且该方法在理论上几乎不会因频谱泄漏而产生偏差。但矿集区噪声往往种类多、形态复杂,当信号中混有突变噪声(如方波、三角波等)时,对信号进行EMD会产生明显的模态混叠效应(不同频率的信号存在于同一分解分量中或者同一频率信号存在于不同分解分量中),尤其是在“死频段”很难保证测得原始信号的信噪比。当测得噪声形态复杂且信噪比低时,造成视电阻率、相位曲线会严重离散,通常处理手段难以取得较好的去噪效果,不能保证“死频带”的大地电磁测深信号的准确性。
发明内容
本发明的目的是提供一种大地电磁测深信号的重构方法及系统,解决现有技术中对于“死频带”的大地电磁测深信号不能实现良好的去噪效果的问题。
为实现上述目的,本发明提供了如下方案:
一种大地电磁测深信号的重构方法,包括:
获取初始的大地电磁测深信号;
对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个模态分量序列;
对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱;
根据所述瞬时频谱的频点的分布情况,筛选所述瞬时频谱对应的模态分量序列;
若所述瞬时频谱频点的分布情况是连续分布且均方差小于均方差阈值,则保留所述瞬时频谱对应的模态分量序列;
若所述瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值,则剔除所述瞬时频谱对应的模态分量序列;
对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。
可选的,所述对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个同一阶模态分量序列,具体包括:
获取m组正负成对的高斯白噪声;
将所述m组正负成对的高斯白噪声与所述初始的大地电磁测深信号进行组合,得到2m个合成信号;
对于第i个合成信号,获取第i个合成信号的第j阶的模态分量;
将所述第i个合成信号减去所述第j阶的模态分量,得到更新后的第i个合成信号;
判断所述更新后的第i个合成信号的极值数目是否大于2,得到第一判断结果;
若所述第一判断结果表示所述更新后的第i个合成信号的极值数目不大于2,则停止分解,得到所述第i个合成信号的不同阶的模态分量;
若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t);
利用公式
Figure BDA0002319453510000031
计算所述更新后的第i个合成信号的第j+1阶对应的包络平均线;
将所述更新后的第i个合成信号减去所述第j+1阶对应的包络平均线,得到第j+1阶的中间线;
判断所述第j+1阶的中间线是否符合确定标准,得到第二判断结果;所述确定标准为所述第j+1阶的中间线的上包络线和下包络线的均值为0且所述中间线的极值点数目与过零点的数目的差值不大于1;
若所述第二判断结果表示所述第j+1阶的中间线不符合确定标准,则将所述更新后的第i个合成信号替换为所述第j+1阶的中间线,返回若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)的步骤;
若所述第二判断结果表示所述第j+1阶的中间线符合确定标准,则将所述第j+1阶的中间线确定为第j+1阶的模态分量;
确定每个合成信号的所有阶的模态分量;
将相同阶的所有模态分量加权平均,得到同一阶的模态分量序列。
可选的,所述对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱,具体包括:
利用公式
Figure BDA0002319453510000041
得到希尔伯特变换后的第j阶的模态分量序列;IMF'j(t)为希尔伯特变换后的第j阶的模态分量序列,IMFj(τ)为第j阶的模态分量序列,π为常数,τ为时间,t为时间;
利用
Figure BDA0002319453510000042
确定第j阶的模态分量序列的瞬时振幅;
利用公式
Figure BDA0002319453510000043
确定第j阶的模态分量序列的瞬时频率;φj(t)瞬时相位,
Figure BDA0002319453510000044
根据所述第j阶的模态分量序列的瞬时频率和所述第j阶的模态分量序列的瞬时振幅,确定第j阶的模态分量序列的瞬时频谱。
可选的,所述对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号,具体包括:
利用公式
Figure BDA0002319453510000045
确定重构后的大地电磁测深信号;y1(t)为重构后的大地电磁测深信号,n为保留的模态分量序列的数目,IMFi为第i个保留的模态分量序列。
一种大地电磁测深信号的重构系统,包括:
初始信号获取模块,用于获取初始的大地电磁测深信号;
模态分量序列确定模块,用于对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个模态分量序列;
瞬时频谱确定模块,用于对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱;
模态分量序列筛选模块,用于根据所述瞬时频谱的频点的分布情况,筛选所述瞬时频谱对应的模态分量序列;
模态分量序列保留模块,用于若所述瞬时频谱频点的分布情况是连续分布且均方差小于均方差阈值,则保留所述瞬时频谱对应的模态分量序列;
模态分量序列剔除模块,用于若所述瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值,则剔除所述瞬时频谱对应的模态分量序列;
信号重构模块,用于对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。
可选的,所述模态分量序列确定模块具体包括:
高斯白噪声获取单元,用于获取m组正负成对的高斯白噪声;
信号合成单元,用于将所述m组正负成对的高斯白噪声与所述初始的大地电磁测深信号进行组合,得到2m个合成信号;
第j阶的模态分量获取单元,用于对于第i个合成信号,获取第i个合成信号的第j阶的模态分量;
合成信号更新单元,用于将所述第i个合成信号减去所述第j阶的模态分量,得到更新后的第i个合成信号;
第一判断模块,用于判断所述更新后的第i个合成信号的极值数目是否大于2,得到第一判断结果;
不同阶的模态分量确定单元,用于若所述第一判断结果表示所述更新后的第i个合成信号的极值数目不大于2,则停止分解,得到所述第i个合成信号的不同阶的模态分量;
包线确定单元,用于若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t);
包络平均线确定单元,用于利用公式
Figure BDA0002319453510000051
计算所述更新后的第i个合成信号的第j+1阶对应的包络平均线;
第一中间线确定单元,用于将所述更新后的第i个合成信号减去所述第j+1阶对应的包络平均线,得到第j+1阶的中间线;
第二判断单元,用于判断所述第j+1阶的中间线是否符合确定标准,得到第二判断结果;所述确定标准为所述第j+1阶的中间线的上包络线和下包络线的均值为0且所述中间线的极值点数目与过零点的数目的差值不大于1;
第二中间线确定单元,用于若所述第二判断结果表示所述第j+1阶的中间线不符合确定标准,则将所述更新后的第i个合成信号替换为所述第j+1阶的中间线,返回若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)的步骤;
第j+1阶的模态分量确定单元,用于若所述第二判断结果表示所述第j+1阶的中间线符合确定标准,则将所述第j+1阶的中间线确定为第j+1阶的模态分量;
所有阶的模态分量确定单元,用于确定每个合成信号的所有阶的模态分量;
模态分量序列确定单元,用于将相同阶的所有模态分量加权平均,得到同一阶的模态分量序列。
可选的,所述瞬时频谱确定模块具体包括:
模态分量序列变换单元,用于利用公式
Figure BDA0002319453510000061
得到希尔伯特变换后的第j阶的模态分量序列;IMF'j(t)为希尔伯特变换后的第j阶的模态分量序列,IMFj(τ)为第j阶的模态分量序列,π为常数,τ为时间常数,t为时间;
瞬时振幅确定单元,用于利用
Figure BDA0002319453510000062
确定第j阶的模态分量序列的瞬时振幅;
瞬时频率确定单元,用于利用公式
Figure BDA0002319453510000063
确定第j阶的模态分量序列的瞬时频率;φj(t)瞬时相位,
Figure BDA0002319453510000064
瞬时频谱确定单元,用于根据所述第j阶的模态分量序列的瞬时频率和所述第j阶的模态分量序列的瞬时振幅,确定第j阶的模态分量序列的瞬时频谱。
可选的,所述信号重构模块具体包括:
信号重构单元,用于利用公式
Figure BDA0002319453510000071
确定重构后的大地电磁测深信号;y1(t)为重构后的大地电磁测深信号,n为保留的模态分量序列的数目,IMFi为第i个保留的模态分量序列。
根据本发明提供的具体实施例,本发明公开了以下技术效果:
本发明所提供的一种大地电磁测深信号的重构方法及系统,通过互补集合经验模态分解,改变初始的大地电磁测深信号的极值分布,从而压制模态混叠效应。互补集合经验模态分解初始的大地电磁测深信号中的噪声部分的突变部分以及剩余的大尺度相对平滑部分可被分解到不同模态分量序列中。再通过希尔伯特变换得到瞬时频谱,得到的瞬时频谱能够精细刻画出不同模态分量序列的频谱随时间变换的特征。当瞬时频谱的频点的分布情况是连续分布且均方差小于均方差阈值时,保留所述瞬时频谱对应的模态分量序列,当瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值时,剔除所述瞬时频谱对应的模态分量序列,通过瞬时频谱的频点的分布情况对模态分量序列进行筛选,可以实现噪声突变部分与大尺度平缓部分的分离以及压制造成电阻率和相位离散的噪声成分。即通过瞬时频谱的频点的分布情况去除初始的大地电磁测深信号的噪声。进而再对所有的保留的模态分量序列进行重构,重构后的大地电磁测深信号取得较好的去噪效果,保证了“死频带”的大地电磁测深信号的准确性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明所提供的一种大地电磁测深信号的重构方法的流程示意图;
图2为本发明所提供的一种大地电磁测深信号的重构系统的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本发明的目的是提供一种大地电磁测深信号的重构方法及系统,解决现有技术中对于“死频带”的大地电磁测深信号不能实现良好的去噪效果的问题。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明所提供的一种大地电磁测深信号的重构方法的流程示意图,如图1所示,本发明所提供的一种大地电磁测深信号的重构方法,包括:
S101,获取初始的大地电磁测深信号。
S102,对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个模态分量序列。
获取m组正负成对的高斯白噪声。
将所述m组正负成对的高斯白噪声与所述初始的大地电磁测深信号进行组合,得到2m个合成信号。通过引入正负成对的高斯白噪声,改变初始的大地电磁测深信号的极值分布,从而压制模态混叠效应。
当高斯白噪声为n(t),初始的大地电磁测深信号为y0时,利用公式
Figure BDA0002319453510000081
得到合成信号。其中,M1和M2分别为添加了正、负高斯白噪声的合成信号。
对于第i个合成信号,获取第i个合成信号的第j阶的模态分量。
将所述第i个合成信号减去所述第j阶的模态分量,得到更新后的第i个合成信号。
判断所述更新后的第i个合成信号的极值数目是否大于2,得到第一判断结果。
若所述第一判断结果表示所述更新后的第i个合成信号的极值数目不大于2,则停止分解,得到所述第i个合成信号的不同阶的模态分量。
其中,当第i个合成信号的分解次数大于设定次数阈值时,第i个合成信号也停止分解。
若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)。
利用公式
Figure BDA0002319453510000091
计算所述更新后的第i个合成信号的第j+1阶对应的包络平均线。
将所述更新后的第i个合成信号减去所述第j+1阶对应的包络平均线,得到第j+1阶的中间线。
判断所述第j+1阶的中间线是否符合确定标准,得到第二判断结果;所述确定标准为所述第j+1阶的中间线的上包络线和下包络线的均值为0且所述中间线的极值点数目与过零点的数目的差值不大于1。
若所述第二判断结果表示所述第j+1阶的中间线不符合确定标准,则将所述更新后的第i个合成信号替换为所述第j+1阶的中间线,返回若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)的步骤。
若所述第二判断结果表示所述第j+1阶的中间线符合确定标准,则将所述第j+1阶的中间线确定为第j+1阶的模态分量。
确定每个合成信号的所有阶的模态分量。
将相同阶的所有模态分量加权平均,得到同一阶的模态分量序列。
其中,利用公式
Figure BDA0002319453510000092
确定第j阶的模态分量序列。
S103,对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱。使用希尔伯特变换得到的瞬时频谱对截断噪声和窄带正、余弦噪声进行拾取。
利用公式
Figure BDA0002319453510000101
得到希尔伯特变换后的第j阶的模态分量序列;IMF'j(t)为希尔伯特变换后的第j阶的模态分量序列,IMFj(τ)为第j阶的模态分量序列,π为常数,τ为时间常数,t为时间。
利用
Figure BDA0002319453510000102
确定第j阶的模态分量序列的瞬时振幅。
利用公式
Figure BDA0002319453510000103
确定第j阶的模态分量序列的瞬时频率;φj(t)瞬时相位,
Figure BDA0002319453510000104
根据所述第j阶的模态分量序列的瞬时频率和所述第j阶的模态分量序列的瞬时振幅,确定第j阶的模态分量序列的瞬时频谱。
模态分量序列的瞬时频谱是以时间t为横轴,以瞬时频率为纵轴,建立坐标系;并将瞬时振幅以散点的形式投在坐标上,并以颜色深浅为表示振幅大小,得到各阶的模态分量序列的瞬时频谱。
S104,根据所述瞬时频谱的频点的分布情况,筛选所述瞬时频谱对应的模态分量序列。
S105,若所述瞬时频谱频点的分布情况是连续分布且均方差小于均方差阈值,则保留所述瞬时频谱对应的模态分量序列。
S106,若所述瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值,则剔除所述瞬时频谱对应的模态分量序列。
S107,对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。
利用公式
Figure BDA0002319453510000105
确定重构后的大地电磁测深信号;y1(t)为重构后的大地电磁测深信号,n为保留的模态分量序列的数目,IMFi为第i个保留的模态分量序列。
本发明获取初始的大地电磁测深信号的噪声经过互补集合经验模态分解多次分解,使其突变部分以及剩余的大尺度相对平滑部分可被分解到不同模态分量中。其中大尺度相对平滑部分能通过窄带正、余弦波进行拟合,由于窄带正、余弦波只对视电阻率、相位曲线的极少数频点产生影响,且受影响的频点不在“死频段”内,只需滤去比重较少的噪声突变部分,这使得即便在信噪比较低的情况下仍能基本恢复视电阻率、相位曲线的变化趋势,且取得圆滑、连续变化的曲线。
本发明基于希尔伯特变换得到的瞬时频谱能够精细刻画出各阶模态分量序列频谱随时间变换的特征,根据频点在时、频轴上的分布状态(大尺度噪声在时间轴上分布宽,在频率轴上分布窄,截断噪声)以及频点的离散情况(有效信号频点连续分布,噪声频点离散程度大、散乱分布)对模态分量序列进行拾取,并使用拾取的模态分量序列进行重构,可以实现噪声突变部分与大尺度平缓部分的分离、压制能造成视电阻率、相位离散的噪声成分。
因此,本发明使用互补集合经验模态分解和希尔伯特变换相结合的方法对大地电磁“死频段”含噪信号进行处理,相比于传统大地电磁去噪方该方法对信噪比要求低,不受噪声相关性的制约,相比于S变换、小波变换、数学形态滤波、压缩感知重构以及HHT没有正交基选择的困难,本发明对信噪比要求低、且适用噪声类型广(不受矿集区噪声复杂形态的影响),相比于Rhoplus校正,本发明所提供的一种大地电磁测深信号的重构方法无需大量观测数据作为支撑以及人机交互操作,因而具有操作简单、结果客观、经济成本低等优点。
对应本发明所提供的一种大地电磁测深信号的重构方法,本发明还提供一种大地电磁测深信号的重构系统,如图2所示,本发明所提供的一种大地电磁测深信号的重构系统包括:初始信号获取模块201、模态分量序列确定模块202、瞬时频谱确定模块203、模态分量序列筛选模块204、模态分量序列保留模块205、模态分量序列剔除模块206和信号重构模块207。
初始信号获取模块201用于获取初始的大地电磁测深信号。
模态分量序列确定模块202用于对所述初始的大地电磁测深信号进行互补集合经验模态分解得到多个模态分量序列。
瞬时频谱确定模块203用于对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱。
模态分量序列筛选模块204用于根据所述瞬时频谱的频点的分布情况,筛选所述瞬时频谱对应的模态分量序列。
模态分量序列保留模块205用于若所述瞬时频谱频点的分布情况是连续分布且均方差小于均方差阈值,则保留所述瞬时频谱对应的模态分量序列。
模态分量序列剔除模块206用于若所述瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值,则剔除所述瞬时频谱对应的模态分量序列。
信号重构模块207用于对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。
所述模态分量序列确定模块202具体包括:高斯白噪声获取单元、信号合成单元、第j阶的模态分量获取单元、合成信号更新单元、第一判断模块、不同阶的模态分量确定单元、包线确定单元、包络平均线确定单元、第一中间线确定单元、第二判断单元、第二中间线确定单元、第j+1阶的模态分量确定单元、所有阶的模态分量确定单元和模态分量序列确定单元。
高斯白噪声获取单元用于获取m组正负成对的高斯白噪声。
信号合成单元用于将所述m组正负成对的高斯白噪声与所述初始的大地电磁测深信号进行组合,得到2m个合成信号。
第j阶的模态分量获取单元用于对于第i个合成信号,获取第i个合成信号的第j阶的模态分量。
合成信号更新单元用于将所述第i个合成信号减去所述第j阶的模态分量,得到更新后的第i个合成信号。
第一判断模块用于判断所述更新后的第i个合成信号的极值数目是否大于2,得到第一判断结果。
不同阶的模态分量确定单元用于若所述第一判断结果表示所述更新后的第i个合成信号的极值数目不大于2,则停止分解,得到所述第i个合成信号的不同阶的模态分量。
包线确定单元用于若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)。
包络平均线确定单元用于利用公式
Figure BDA0002319453510000131
计算所述更新后的第i个合成信号的第j+1阶对应的包络平均线。
第一中间线确定单元用于将所述更新后的第i个合成信号减去所述第j+1阶对应的包络平均线,得到第j+1阶的中间线。
第二判断单元用于判断所述第j+1阶的中间线是否符合确定标准,得到第二判断结果;所述确定标准为所述第j+1阶的中间线的上包络线和下包络线的均值为0且所述中间线的极值点数目与过零点的数目的差值不大于1。
第二中间线确定单元用于若所述第二判断结果表示所述第j+1阶的中间线不符合确定标准,则将所述更新后的第i个合成信号替换为所述第j+1阶的中间线,返回若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)的步骤。
第j+1阶的模态分量确定单元用于若所述第二判断结果表示所述第j+1阶的中间线符合确定标准,则将所述第j+1阶的中间线确定为第j+1阶的模态分量。
所有阶的模态分量确定单元用于确定每个合成信号的所有阶的模态分量。
模态分量序列确定单元用于将相同阶的所有模态分量加权平均,得到同一阶的模态分量序列。
所述瞬时频谱确定模块203具体包括:模态分量序列变换单元、瞬时振幅确定单元、瞬时频率确定单元和瞬时频谱确定单元。
模态分量序列变换单元用于利用公式
Figure BDA0002319453510000132
得到希尔伯特变换后的第j阶的模态分量序列;IMF'j(t)为希尔伯特变换后的第j阶的模态分量序列,IMFj(τ)为第j阶的模态分量序列,π为常数,τ为时间常数,t为时间。
瞬时振幅确定单元用于利用
Figure BDA0002319453510000141
确定第j阶的模态分量序列的瞬时振幅。
瞬时频率确定单元用于利用公式
Figure BDA0002319453510000142
确定第j阶的模态分量序列的瞬时频率;φj(t)瞬时相位,
Figure BDA0002319453510000143
瞬时频谱确定单元用于根据所述第j阶的模态分量序列的瞬时频率和所述第j阶的模态分量序列的瞬时振幅,确定第j阶的模态分量序列的瞬时频谱。
所述信号重构模块207具体包括:信号重构单元。
信号重构单元用于利用公式
Figure BDA0002319453510000144
确定重构后的大地电磁测深信号;y1(t)为重构后的大地电磁测深信号,n为保留的模态分量序列的数目,IMFi为第i个保留的模态分量序列。
本说明书中各个实施例采用递进的方式描述,每个实施例重点说明的都是与其他实施例的不同之处,各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本发明中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处。综上所述,本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种大地电磁测深信号的重构方法,其特征在于,包括:
获取初始的大地电磁测深信号;
对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个模态分量序列;
对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱;
根据所述瞬时频谱的频点的分布情况,筛选所述瞬时频谱对应的模态分量序列;
若所述瞬时频谱频点的分布情况是连续分布且均方差小于均方差阈值,则保留所述瞬时频谱对应的模态分量序列;
若所述瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值,则剔除所述瞬时频谱对应的模态分量序列;
对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。
2.根据权利要求1所述的一种大地电磁测深信号的重构方法,其特征在于,所述对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个同一阶模态分量序列,具体包括:
获取m组正负成对的高斯白噪声;
将所述m组正负成对的高斯白噪声与所述初始的大地电磁测深信号进行组合,得到2m个合成信号;
对于第i个合成信号,获取第i个合成信号的第j阶的模态分量;
将所述第i个合成信号减去所述第j阶的模态分量,得到更新后的第i个合成信号;
判断所述更新后的第i个合成信号的极值数目是否大于2,得到第一判断结果;
若所述第一判断结果表示所述更新后的第i个合成信号的极值数目不大于2,则停止分解,得到所述第i个合成信号的不同阶的模态分量;
若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t);
利用公式
Figure FDA0002319453500000021
计算所述更新后的第i个合成信号的第j+1阶对应的包络平均线;
将所述更新后的第i个合成信号减去所述第j+1阶对应的包络平均线,得到第j+1阶的中间线;
判断所述第j+1阶的中间线是否符合确定标准,得到第二判断结果;所述确定标准为所述第j+1阶的中间线的上包络线和下包络线的均值为0且所述中间线的极值点数目与过零点的数目的差值不大于1;
若所述第二判断结果表示所述第j+1阶的中间线不符合确定标准,则将所述更新后的第i个合成信号替换为所述第j+1阶的中间线,返回若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)的步骤;
若所述第二判断结果表示所述第j+1阶的中间线符合确定标准,则将所述第j+1阶的中间线确定为第j+1阶的模态分量;
确定每个合成信号的所有阶的模态分量;
将相同阶的所有模态分量加权平均,得到同一阶的模态分量序列。
3.根据权利要求1所述的一种大地电磁测深信号的重构方法,其特征在于,所述对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱,具体包括:
利用公式
Figure FDA0002319453500000022
得到希尔伯特变换后的第j阶的模态分量序列;IMF'j(t)为希尔伯特变换后的第j阶的模态分量序列,IMFj(τ)为第j阶的模态分量序列,π为常数,τ为时间常数,t为时间;
利用
Figure FDA0002319453500000023
确定第j阶的模态分量序列的瞬时振幅;
利用公式
Figure FDA0002319453500000024
确定第j阶的模态分量序列的瞬时频率;φj(t)瞬时相位,
Figure FDA0002319453500000031
根据所述第j阶的模态分量序列的瞬时频率和所述第j阶的模态分量序列的瞬时振幅,确定第j阶的模态分量序列的瞬时频谱。
4.根据权利要求1所述的一种大地电磁测深信号的重构方法,其特征在于,所述对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号,具体包括:
利用公式
Figure FDA0002319453500000032
确定重构后的大地电磁测深信号;y1(t)为重构后的大地电磁测深信号,n为保留的模态分量序列的数目,IMFi为第i个保留的模态分量序列。
5.一种大地电磁测深信号的重构系统,其特征在于,包括:
初始信号获取模块,用于获取初始的大地电磁测深信号;
模态分量序列确定模块,用于对所述初始的大地电磁测深信号进行互补集合经验模态分解,得到多个模态分量序列;
瞬时频谱确定模块,用于对多个所述模态分量序列进行希尔伯特变换,得到每个所述模态分量序列的瞬时频谱;
模态分量序列筛选模块,用于根据所述瞬时频谱的频点的分布情况,筛选所述瞬时频谱对应的模态分量序列;
模态分量序列保留模块,用于若所述瞬时频谱频点的分布情况是连续分布且均方差小于均方差阈值,则保留所述瞬时频谱对应的模态分量序列;
模态分量序列剔除模块,用于若所述瞬时频谱频点的分布情况不是连续分布或均方差不小于均方差阈值,则剔除所述瞬时频谱对应的模态分量序列;
信号重构模块,用于对所有的保留的模态分量序列进行重构,得到重构后的大地电磁测深信号。
6.根据权利要求5所述的一种大地电磁测深信号的重构系统,其特征在于,所述模态分量序列确定模块具体包括:
高斯白噪声获取单元,用于获取m组正负成对的高斯白噪声;
信号合成单元,用于将所述m组正负成对的高斯白噪声与所述初始的大地电磁测深信号进行组合,得到2m个合成信号;
第j阶的模态分量获取单元,用于对于第i个合成信号,获取第i个合成信号的第j阶的模态分量;
合成信号更新单元,用于将所述第i个合成信号减去所述第j阶的模态分量,得到更新后的第i个合成信号;
第一判断模块,用于判断所述更新后的第i个合成信号的极值数目是否大于2,得到第一判断结果;
不同阶的模态分量确定单元,用于若所述第一判断结果表示所述更新后的第i个合成信号的极值数目不大于2,则停止分解,得到所述第i个合成信号的不同阶的模态分量;
包线确定单元,用于若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t);
包络平均线确定单元,用于利用公式
Figure FDA0002319453500000041
计算所述更新后的第i个合成信号的第j+1阶对应的包络平均线;
第一中间线确定单元,用于将所述更新后的第i个合成信号减去所述第j+1阶对应的包络平均线,得到第j+1阶的中间线;
第二判断单元,用于判断所述第j+1阶的中间线是否符合确定标准,得到第二判断结果;所述确定标准为所述第j+1阶的中间线的上包络线和下包络线的均值为0且所述中间线的极值点数目与过零点的数目的差值不大于1;
第二中间线确定单元,用于若所述第二判断结果表示所述第j+1阶的中间线不符合确定标准,则将所述更新后的第i个合成信号替换为所述第j+1阶的中间线,返回若所述第一判断结果表示所述更新后的第i个合成信号的极值数目大于2,则确定更新后的第i个合成信号的上包络线xmax(t)和下包络线xmin(t)的步骤;
第j+1阶的模态分量确定单元,用于若所述第二判断结果表示所述第j+1阶的中间线符合确定标准,则将所述第j+1阶的中间线确定为第j+1阶的模态分量;
所有阶的模态分量确定单元,用于确定每个合成信号的所有阶的模态分量;
模态分量序列确定单元,用于将相同阶的所有模态分量加权平均,得到同一阶的模态分量序列。
7.根据权利要求5所述的一种大地电磁测深信号的重构系统,其特征在于,所述瞬时频谱确定模块具体包括:
模态分量序列变换单元,用于利用公式
Figure FDA0002319453500000051
得到希尔伯特变换后的第j阶的模态分量序列;IMF'j(t)为希尔伯特变换后的第j阶的模态分量序列,IMFj(τ)为第j阶的模态分量序列,π为常数,τ为时间常数,t为时间;
瞬时振幅确定单元,用于利用
Figure FDA0002319453500000052
确定第j阶的模态分量序列的瞬时振幅;
瞬时频率确定单元,用于利用公式
Figure FDA0002319453500000053
确定第j阶的模态分量序列的瞬时频率;φj(t)瞬时相位,
Figure FDA0002319453500000054
瞬时频谱确定单元,用于根据所述第j阶的模态分量序列的瞬时频率和所述第j阶的模态分量序列的瞬时振幅,确定第j阶的模态分量序列的瞬时频谱。
8.根据权利要求5所述的一种大地电磁测深信号的重构系统,其特征在于,所述信号重构模块具体包括:
信号重构单元,用于利用公式
Figure FDA0002319453500000055
确定重构后的大地电磁测深信号;y1(t)为重构后的大地电磁测深信号,n为保留的模态分量序列的数目,IMFi为第i个保留的模态分量序列。
CN201911292329.9A 2019-12-16 2019-12-16 一种大地电磁测深信号的重构方法及系统 Expired - Fee Related CN110908001B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911292329.9A CN110908001B (zh) 2019-12-16 2019-12-16 一种大地电磁测深信号的重构方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911292329.9A CN110908001B (zh) 2019-12-16 2019-12-16 一种大地电磁测深信号的重构方法及系统

Publications (2)

Publication Number Publication Date
CN110908001A true CN110908001A (zh) 2020-03-24
CN110908001B CN110908001B (zh) 2020-09-29

Family

ID=69825603

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911292329.9A Expired - Fee Related CN110908001B (zh) 2019-12-16 2019-12-16 一种大地电磁测深信号的重构方法及系统

Country Status (1)

Country Link
CN (1) CN110908001B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2022061596A1 (zh) * 2020-09-23 2022-03-31 深圳市速腾聚创科技有限公司 一种信号噪声滤除方法、装置、存储介质及激光雷达
CN115563791A (zh) * 2022-10-14 2023-01-03 吉林大学 基于压缩感知重构的大地电磁数据反演方法
CN117538944A (zh) * 2023-11-24 2024-02-09 长江大学 一种大地电磁功率谱的挑谱方法
CN117706646A (zh) * 2024-02-06 2024-03-15 山东万洋石油科技有限公司 一种测井仪电阻率测量优化校准方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5983162A (en) * 1996-08-12 1999-11-09 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Computer implemented empirical mode decomposition method, apparatus and article of manufacture
US7941298B2 (en) * 2006-09-07 2011-05-10 DynaDx Corporation Noise-assisted data analysis method, system and program product therefor
US20120089372A1 (en) * 2010-10-06 2012-04-12 Industrial Technology Research Institute Apparatus and method for adaptive time-frequency analysis
CN104568024A (zh) * 2015-01-21 2015-04-29 山东理工大学 振动式流量计特征信号提取方法
CN106650218A (zh) * 2016-10-19 2017-05-10 上海电机学院 基于ceemd算法和希尔伯特变换的谐波分析方法
CN107102356A (zh) * 2017-06-02 2017-08-29 成都理工大学 基于ceemd的地震信号高分辨率处理方法
CN107133606A (zh) * 2017-05-26 2017-09-05 广东工业大学 一种时间序列趋势提取的方法及装置
CN108345033A (zh) * 2018-01-26 2018-07-31 中国石油大学(华东) 一种微地震信号时频域初至检测方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5983162A (en) * 1996-08-12 1999-11-09 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Computer implemented empirical mode decomposition method, apparatus and article of manufacture
US7941298B2 (en) * 2006-09-07 2011-05-10 DynaDx Corporation Noise-assisted data analysis method, system and program product therefor
US20120089372A1 (en) * 2010-10-06 2012-04-12 Industrial Technology Research Institute Apparatus and method for adaptive time-frequency analysis
CN104568024A (zh) * 2015-01-21 2015-04-29 山东理工大学 振动式流量计特征信号提取方法
CN106650218A (zh) * 2016-10-19 2017-05-10 上海电机学院 基于ceemd算法和希尔伯特变换的谐波分析方法
CN107133606A (zh) * 2017-05-26 2017-09-05 广东工业大学 一种时间序列趋势提取的方法及装置
CN107102356A (zh) * 2017-06-02 2017-08-29 成都理工大学 基于ceemd的地震信号高分辨率处理方法
CN108345033A (zh) * 2018-01-26 2018-07-31 中国石油大学(华东) 一种微地震信号时频域初至检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
鄢高韩 等: ""CEEMD高分辨率时频分析方法研究与应用"", 《地球物理学进展》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2022061596A1 (zh) * 2020-09-23 2022-03-31 深圳市速腾聚创科技有限公司 一种信号噪声滤除方法、装置、存储介质及激光雷达
CN115563791A (zh) * 2022-10-14 2023-01-03 吉林大学 基于压缩感知重构的大地电磁数据反演方法
CN117538944A (zh) * 2023-11-24 2024-02-09 长江大学 一种大地电磁功率谱的挑谱方法
CN117538944B (zh) * 2023-11-24 2024-04-05 长江大学 一种大地电磁功率谱的挑谱方法
CN117706646A (zh) * 2024-02-06 2024-03-15 山东万洋石油科技有限公司 一种测井仪电阻率测量优化校准方法
CN117706646B (zh) * 2024-02-06 2024-04-19 山东万洋石油科技有限公司 一种测井仪电阻率测量优化校准方法

Also Published As

Publication number Publication date
CN110908001B (zh) 2020-09-29

Similar Documents

Publication Publication Date Title
CN110908001B (zh) 一种大地电磁测深信号的重构方法及系统
Chen et al. Distributed acoustic sensing coupling noise removal based on sparse optimization
CN108845352B (zh) 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法
CN102998706B (zh) 一种衰减地震数据随机噪声的方法及系统
CN101482617A (zh) 基于非下采样轮廓波的合成孔径雷达图像去噪方法
CN106199532B (zh) 基于混合傅立叶-小波分析的探地雷达信号降噪方法
CN102681014A (zh) 基于多项式拟合的规则线性干扰压制方法
Li et al. Wavelet-based higher order correlative stacking for seismic data denoising in the curvelet domain
CN113297987B (zh) 一种基于双目标函数优化的变分模态分解信号降噪方法
CN115700544A (zh) 一种联合经验模态分解及小波软阈值的色谱信号去噪方法
CN109359633B (zh) 基于希尔伯特-黄变换和小波脊线的信号联合分类方法
CN113723171A (zh) 基于残差生成对抗网络的脑电信号去噪方法
Ahmadi et al. Types of EMD algorithms
CN110989020A (zh) 一种音频大地电磁数据噪声干扰的滤波方法及系统
Li et al. Magnetotelluric signal-noise separation method based on SVM–CEEMDWT
Li et al. A hybrid filtering method based on a novel empirical mode decomposition for friction signals
Luo et al. Denoising of medical images using a reconstruction-average mechanism
CN109460614A (zh) 基于瞬时带宽的信号时间-频率分解方法
CN106125132B (zh) 含单频干扰地震道的迭代识别和压制方法
CN114449410A (zh) 一种多通道声纹信号同步采集系统及方法
Tayade et al. Medical image denoising and enhancement using DTCWT and Wiener filter
CN116029131B (zh) 一种海上风电场水下金属结构磁信号模拟和处理方法
Meng et al. Design of seismic data acquisition system based on improved wavelet threshold de-noising
Feng et al. Binarization method for on-line ferrograph image based on uniform curvelet transformation
Liu et al. The application of Hilbert-Huang Transform in speech enhancement

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200929

CF01 Termination of patent right due to non-payment of annual fee