一种在时频域获取大地电磁阻抗的方法技术

本发明专利技术提供了一种在时频域获取大地电磁阻抗的方法，包括以下步骤：步骤一：输入信号各道数据x

【技术实现步骤摘要】
一种在时频域获取大地电磁阻抗的方法


[0001]本专利技术涉及地球物理勘探
，特别是涉及一种在时频域获取大地电磁阻抗的方法。

技术介绍

[0002]大地电磁测深法是地球物理勘探中的一种重要方法，用于地下深部电性结构的勘探，在油气藏普查、火山活动监测、地热资源勘探以及大区域构造研究中具有广泛的应用。大地电磁测深采用天然场源，因而具有经济成本低、仪器轻便的优点，然而天然电磁场源的使用也产生了一个问题——信号与噪声水平是非平稳的，甚至高度是可变的。这使得在频率域内计算大地电磁响应时，非平稳的含噪MT信号可能产生畸变明显的或者剧烈震荡的视电阻率
‑
相位曲线。当测得的MT信号信噪比较低或时长较短时，响应曲线的畸变或震荡尤为严重，在实际测量中常常需要对电磁场信号进行长时间、低噪声的观测，这使得野外工作难度增大且人力物力成本增加。
[0003]自大地电磁测深法提出以来，如何获取可靠的阻抗估计一直困扰着广大研究者。前人对此进行了大量的研究工作，并提出了诸多的阻抗估计方法，根据这些方法处理对象的不同可分为频率域以及时间域阻抗估计方法。
[0004]频率域阻抗估计方法是在频率域内基于傅里叶频谱估计大地电磁阻抗的一类方法，是大地电磁阻抗估计研究中发展最早的，也是诸多方法中应用最广泛的一类方法。Swift(1967)提出了基于最小二乘的MT阻抗估计方法，为频率域阻抗估计的研究奠定了基础，但该方法抗噪声干扰能力差，难以适应实际数据采集中复杂的噪声环境。远参考(Remote Reference，RR)与robust阻抗估计方法的发展有效地改善了这一情况，其中Gamble et al.(1979))提出的RR方法通过设置远参考站能够有效压制相干噪声的影响，而Sutarno(1989，1990)、Larsen et al.(1996)、Egbert(1997)、Chave and Thomson(2004)以及Sutarno(2008)发展的多种自适应加权的robust方法，可以有效压制明显不合理的飞点对大地电磁阻抗估计结果的影响。此外Chave(2014，2017)对阻抗估计中的高斯核进行了研究，Eisel and Egbert(2001)、Wawrzyniak et al.(2013)研究了阻抗估计中的误差估计方法，Oettinger et al.(2001)、Weckmann et al.(2005)、Sato et al.(2018)以及Li et al.(2021)提出了多种大地电磁去噪方法，在频率域与时间域内的对异常谱值与异常波形进行压制。这些阻抗估计方法与去噪方法的发展使得对低信噪比MT数据的稳健阻抗估计成为可能。然而由于这类阻抗估计方法是在频率域中进行的，傅里叶变换的平稳性假设使得仍需足够的观测时间来获得可靠的电磁场频谱。这在中高频的阻抗估计中可能是微不足道的，但当涉及到低频段时，往往需要进行长周期的观测，从而产生较高的观测成本。
[0005]Kunetz(1972)提出了在时间域中进行阻抗估计的策略，通过直接对时间序列进行阻抗估计，有效避免傅里叶变换对电磁场序列的平稳性假设。在此基础上，Mecmechan and Barrodale(1985)、Yee et al.(1988)、Spagnolini(1994)以及(2004)相继展开了对时间域阻抗估计的研究，并逐步将时间域阻抗估计发展成了一套自成体系的处理方
法。然而由于时间域与频率域运算符的对应关系，频率域中进行的矩阵相乘在时间域中变成了序列卷积，这使得时间域阻抗估计方法计算量大幅大，计算复杂度高，同时阻抗估计结果受噪声类型影响大，在实际应用中难以取得理想的应用效果。
[0006]传统频率域阻抗估计方法对信噪比要求高、所需观测时长较长，但阻抗估计过程计算量小，收敛速度快。时间域阻抗估计方法降低了频率域阻抗估计对信噪比与观测时长的严格要求，但计算复杂度高难以适用于长序列数据低频阻抗的求解。

技术实现思路

[0007]为了克服现有技术的不足，本专利技术的目的是提供一种在时频域获取大地电磁阻抗的方法，具有计算速度快、对电磁场信号观测时长要求低以及抗噪声能力强的优点。
[0008]为实现上述目的，本专利技术提供了一种在时频域获取大地电磁阻抗的方法，包括以下步骤：
[0009]步骤一：输入信号各道数据x
i
，其中i＝1，2，3，4，5；
[0010]步骤二：使用变分模态分解处理x
i
得到u
ik
；
[0011]步骤三：对u
ik
进行Hilbert变换得到
[0012]步骤四：基于u
ik
与求得瞬时振幅A
ik
、瞬时相位以及瞬时频率ω
ik
；
[0013]步骤五：根据A
ik
与ω
ik
求得信号边际谱H
i
(ω)；
[0014]步骤六：使用H
i
(ω)计算功率谱<X，Y>；
[0015]步骤七：将功率谱代入不同阻抗表达式中求得阻抗分量Z
xy
和Z
yx
；
[0016]步骤八：对不同阻抗表达式所得结果进行平均求得阻抗估计值和
[0017]步骤九：基于所得阻抗估计值，计算视电阻率与相位。
[0018]优选的，在所述步骤一中，从大地电磁测深仪器中导出大地电磁信号到计算机中并以n行5列矩阵的形式储存数据，其中n＝1，2，3，4
……
。
[0019]优选的，在所述步骤二中，取步骤一所得的矩阵的第i列记为x
i
(i＝1，2，3，4)并对其进行变分模态分解，得到各阶IMF分量记为u
ik
，其中u
ik
为第i列信号的第k阶分量，k＝1，
……
，K；K为预设的分解分量数，为了均衡计算量和分解效果，K取值为4。
[0020]优选的，在所述步骤三中，取步骤二所得分量u
ik
进行Hilbert变换得到
[0021][0022]式中，m为积分变量，t为时间。
[0023]优选的，在所述步骤四中，利用以下公式求得：
[0024][0025][0026][0027]式中，d为微分算子，arctan为通过反正切获取角度，u
ik
'与分别为u
ik
与对t
的导数，即du
ik
/dt与其中，求得的瞬时振幅A
ik
(t)与瞬时频率ω
ik
(t)都随t同步变化，将A
ik
看作是t与ω
ik
的共同函数，即A
ik
(t,ω)。
[0028]优选的，在所述步骤五中，取步骤四所得结果中i与ω相同而t不同的A
ik
(t,ω)记为A
i本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种在时频域获取大地电磁阻抗的方法，其特征在于，包括以下步骤：步骤一：输入信号各道数据x
i
，其中i＝1，2，3，4，5；步骤二：使用变分模态分解处理x
i
得到u
ik
；步骤三：对u
ik
进行Hilbert变换得到步骤四：基于u
ik
与求得瞬时振幅A
ik
、瞬时相位以及瞬时频率ω
ik
；步骤五：根据A
ik
与ω
ik
求得信号边际谱H
i
(ω)；步骤六：使用H
i
(ω)计算功率谱<X，Y>；步骤七：将功率谱代入不同阻抗表达式中求得阻抗分量Z
xy
和Z
yx
；步骤八：对不同阻抗表达式所得结果进行平均求得阻抗估计值和步骤九：基于所得阻抗估计值，计算视电阻率与相位。2.根据权利要求1所述的一种在时频域获取大地电磁阻抗的方法，其特征在于，在所述步骤一中，从大地电磁测深仪器中导出大地电磁信号到计算机中并以n行5列矩阵的形式储存数据，其中n＝1，2，3，4
……
。3.根据权利要求1所述的一种在时频域获取大地电磁阻抗的方法，其特征在于，在所述步骤二中，取步骤一所得的矩阵的第i列记为x
i
(i＝1，2，3，4)并对其进行变分模态分解，得到各阶IMF分量记为u
ik
，其中u
ik
为第i列信号的第k阶分量，k＝1，
……
，K；K为预设的分解分量数，为了均衡计算量和分解效果，K取值为4。4.根据权利要求1所述的一种在时频域获取大地电磁阻抗的方法，其特征在于，在所述步骤三中，取步骤二所得分量u
ik
进行Hilbert变换得到进行Hilbert变换得到式中，m为积分变量，t为时间。5.根据权利要求1所述的一种在时频域获取大地电磁阻抗的方法，其特征在于，在所述步骤四中，利用以下公式求得：步骤四中，利用以下公式求得：步骤四中，利用以下公式求得：式中，d为微分算子，arctan为通过反正切获取角度，u
ik
'与分别为u
ik
与对t的导数，即du
ik
/dt与其中，求得的瞬时振幅A
ik
(t)与瞬时频率ω
ik
(t)都随t同步变化，将A
ik
看作是t与ω
ik
的共同函数，即A
ik
(t,ω)。6.根据权利要求1所述的...

【专利技术属性】
技术研发人员：辛中华，韩江涛，刘立家，郭振宇，
申请(专利权)人：吉林大学，
类型：发明
国别省市：