CN106443801A - 一种大地电磁阻抗估计的时频分析方法 - Google Patents

Abstract

本发明公开了一种大地电磁阻抗估计的时频分析方法，通过采用Hilbert‑Huang变换中的经验模态分解技术，将大地电磁信号分解得到的各阶模态函数进行Hilbert变换，得到大地电磁信号4个分量的时频谱将得到的瞬时幅值代替传统傅立叶方法中的傅氏谱，带入到最小二次方程，求解时频点上的阻抗张量，计算出各个时频点的阻抗张量元素，把大地电磁信号的时频功率谱图转变为时频阻抗谱图；时频阻抗谱图中，将具有相同瞬时频率的各时频点的响应函数值累加起来，除以时频点的个数，其均值则为频率值的响应函数估计值。本发明有利于实现稳健估计，由此方法计算得到的阻抗元素更符合地层结构的实际；有更强的噪声抑制能力，估算的结果更具稳健性。

Description

一种大地电磁阻抗估计的时频分析方法

技术领域

本发明涉及大地电磁测深（Magnetotelluric，MT）资料处理领域，具体说是一种能提高大地电磁阻抗估计精度和可解释性的时频分析方法。

背景技术

大地电磁测深（MT）是一种以天然交变电磁场为场源的电磁勘探方法，其电磁场信号弱、频带宽、极易受各种噪声的干扰，是典型的非线性、非平稳信号。处理MT信号最传统的方法是Fourier变换，而Fourier变换是以平稳信号为理论基础的，与MT信号的非平稳特性不相适应，因此用Fourier变换分析MT信号有着明显的缺陷。在大地电磁测深法中，数据分析、定量反演和地质解释都必须建立在正确可靠的大地电磁响应函数的基础上才有意义，在没有确定MT响应函数是否可靠之前，去追求定量反演的精度意义不大。怎样获得无偏的阻抗估算，是学者们研究MT测深的重要课题。所幸MT信号的处理技术与其它学科的发展息息相关，现代科学技术每发展一步，都可能会带动MT处理技术前进一大步。

Hilbert-Huang变换（Hilbert-Huang transform，HHT）是最新发展起来的处理非线性、非平稳信号的时频分析技术，主要是利用经验模式分解( Empirical ModeDecomposition，简称 EMD)对Hilbert变换中的大多数非平稳信号进行平稳化的时频处理。HHT在许多领域得到了不同程度的研究与应用，已被证明是分析非平稳信号的强有力工具，也为分析非平稳的大地电磁信号提供了一种选择。提高大地电磁阻抗估计精度和可解释性提供条件。

发明内容

针对上述现有技术存在的不足，本发明的目的是提供一种基Hilbert-Huang变换的能提高大地电磁阻抗估计精度的的时频分析方法。

为实现上述目的，本发明采用的技术方案如下：

一种大地电磁阻抗估计的时频分析方法，包括如下步骤：

（1）将大地电磁信号的4个分量H _x 、E _y 、H _y 、E _x 分别都进行如下操作：采用Hilbert-Huang变换中的经验模态分解技术，将大地电磁信号的分量进行分解，自适应地得到各阶模态函数（intrinsic mode function，IMF），然后对各阶IMF进行Hilbert变换，计算得到各阶模态函数的瞬时幅度a _i (t)和瞬时频率ω _i (t)，综合各阶IMF的瞬时谱，得到大地电磁的时频谱H(ω, t )，这样就可以分别得到大地电磁信号4个分量的时频谱。

（2）根据时频谱H(ω, t )，求得任意时频点上信号的瞬时幅值，将得到的瞬时幅值代替传统傅立叶方法中的傅氏谱，带入到最小二次方程，求解时频点上的阻抗张量，数学模型如下：

（1），

（2），

上述， ,,,是具有相同瞬时频率的4个电磁场分量在时频点上的瞬时幅值；和是上述4个电磁场分量中任意两道数据的共轭复数；

（3）瞬时阻抗函数的分量根据公式（1）-（2）可求得， （3），

类同公式（3），其他几个分量,,也同理得到；然后计算出各个时频点的阻抗张量元素，把大地电磁信号的时频功率谱图转变为时频阻抗谱图；时频阻抗谱图中，将具有相同瞬时频率(ω₀)的各时频点的响应函数值累加起来，除以时频点的个数，其均值则为频率值(ω₀)的响应函数估计值：

（4），

然后利用、分别求得电阻率和相位，式中ij表示下标xx, xy, yx, yy；f表示频率。

相比现有技术而言，本发明的有益效果是：1、利用时频分析从瞬时谱上统计估算大地电磁的响应参量，比Fourier分析更有利于实现稳健估计；2、本发明没有要求信号频率成分在整个时间段都平稳，而其时频阻抗值是由真实存在的各时频点谱值计算得到的，相比傅立叶谱，减小了信号非平稳性对谱估计带来的误差，更接近大地电磁信号的真实谱，由此方法计算得到的阻抗元素更符合地层结构的实际；3、不用对大地电磁信号进行窄带滤波，瞬时频率的频率分辨率更高；4、在时频阻抗谱图中，将具有相同瞬时频率(ω₀)的各时频点的响应函数值累加起来，求其均值作为频率值(ω₀)的响应函数估计值，其过程其实质是个整体平均的过程，算法上就有抑制噪声的能力，所以基于Hilbert-Huang的时频阻抗估计方法比基于傅立叶谱的传统方法有更强的噪声抑制能力，估算的结果更具稳健性。

附图说明

图1为本发明工艺流程图

图2为本发明所述经验模态分解法流程图

图3为本发明所述大地电磁4个分量Hy、Ex、Hx、Ey的EMD分解图，其中（a）为Hy分量上的EMD分解图，(b)为Ex分量上的EMD分解图，(c)为Hx分量上的EMD分解图，(d)为Ey分量上的EMD分解图。

图4为大地电磁4个分量的Hilbert时频谱：其中（a）为H _y 分量上的时频谱，(b)为E _x 分量上的时频谱，(c)为H _x 分量上的时频谱，(d)为E _Y 分量上的时频谱。

图5为估算的电阻率ρ _xy：(a)为电阻率的时频分布图，表征了其时变特征；(b)为(a)的局部放大图(time: 0~0.0041s; frequency: 0~337.5Hz ); (c)是由(a)获得的电阻率曲线（由具有相同瞬时频率的各点视电阻率叠加再求平均得到）；(d)为从图(b)中选出来的47.5Hz瞬时频率的时变电阻率。

图6为估算的相位Φ _xy： (a) 为相位的时频分布，表征现了其时变特征； (b) 为(a)的局部放大图 (time: 0~0.0041s; frequency: 0~337.5Hz ); (c) 由图(a)获得的相位曲线(由具有相同瞬时频率的各点相位叠加再求平均得到)； (d) 从图 (b)中选出来的47.5Hz瞬时频率的时变相位。

图7为传统Fourier 变换方法和本发明的处理结果对比：(a) 视电阻率曲线 ρ _xy的对比； (b) 视电阻率曲线 ρ _yx 的对比；(c) 相位曲线 Φ _xy 的对比； (d) 相位曲线Φ _yx 的对比。

具体实施方式

现结合具体实施例，来对本发明作进一步描述。以下实施例仅为本发明的优选方式，并非只用于限制本发明的保护范围，任意与本发明相同或相似的技术方案均视为本发明的保护范围。

实施例一

如图1所示，大地电磁阻抗估计的时频分析方法，包括如下步骤：

（1）采用经验模态分解法，将大地电磁信号的4个分量H _x 、E _y 、H _y 、E _x 分别进行分解。分解过程如图2所示，分别设定大地电磁信号每个分量上的时间序列信号为，首先确定出上的所有极值点，然后将所有极值点用曲线拟合，得到信号的上、下包络曲线，分别设为和，则上、下包络的平均曲线为：； 用减去后，得到剩余部分设为，即：。

理想情况下，应该是一个固有模态函数，因为的构造过程就是使它满足固有模态函数的条件。但由于包络曲线逼近的过冲和俯冲作用，会产生新的极值，影响原来极值的位置与大小，因此，分解得到的并不完全满足IMF条件。虽然这种影响是间接的，因为参与筛选的是上、下包络线的平均曲线。即使拟合过程做得非常好，信号斜坡上的一个微小凸点或凹点，也会把局部零均值从直线变为曲线时被放大为新的局部极值点，这些新的极值点正是前一次筛选过程中遗漏的。

因此，需要反复筛选以恢复所有低幅值的叠加波。筛选过程的目的有两个：一方面消除信号上的骑波，另一方面使波形轮廓更加对称。用代替，用相应的上、下包络线和重复上述过程，即：

。

上述筛分过程重复k次，当hk(t)满足以下两个条件：(1)整个信号中，零点数与极点数相等或至多相差1；(2)信号上任意一点，由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为零，即信号关于时间轴局部对称，则hk(t)就可定义为原始信号的第一个IMF，设为C1（t），信号的剩余部分设为，即：

;

是从原始数据中处理得到的第1个固有模态函数，为原始信号中最短的周期分量，即频率最高的分量。为把从原始信号中分离出来后得到的剩余分量。 把剩余部分作为新的原信号重复以上过程，继续进行EMD分解，依次提取第2个固有模态函数，如此重复n次，直到所得的剩余部分为一单调信号或前后两个剩余项和的标准偏差(standard deviation，SD) 满足0.05<SD<0.3时分解完毕，SD定义如下：

此时，得到第n个固有模态函数和剩余分量：

…

；

这样，原始的信号可表示为所有的IMF及剩余量之和：。大地电磁信号的4个分量H _y 、E _x 、H _x 、E _y 的EMD分解结果如图3所示。

（2）然后对大地电磁信号的4个分量H _y 、E _x 、H _x 、E _y 分别分解得到的模态函数进行Hilbert变换：对于某一个模态函数，设为，首先求出其Hilbert变换及相应的解析信号：

；

和共同组合为一解析信号 ，采用极坐标形式表示：

，

这里瞬时幅度； 。

解析信号的极坐标形式反映了Hilbert变换的物理含义：它是通过一正弦曲线的频率和幅值调制获得信号局部最佳逼近。其相应的瞬时频率为：。对每一分量上的模态函数作上述Hilbert变换，并求解出相应解析信号的瞬时幅度a _i (t)和瞬时频率ω _i (t)，从而原始信号可以表示为：，其中Re表示取实部。这里不考虑剩余分量，因为它是一个平均趋势或者是一个常量。相同的数据若用Fourier形式展开是 。

（3）综合各分量IMF的瞬时谱，得到大地电磁的时频谱H(ω, t )，，如图4所示。

（4）根据时频谱H(ω, t )，求得任意时频点上信号的瞬时幅值，将得到的瞬时幅值代替传统傅立叶方法中的傅氏谱，带入到最小二次方程，求解时频点上的阻抗张量，公式如下：（1） ；

（2）；

上述， ,,,是具有相同瞬时频率的4个电磁场分量在时频点上的瞬时幅值；和是上述4个电磁场分量中任意两道数据的共轭复数；是和他的共轭复数在时频点的乘积, 定义为分量在时频点的瞬时自功率谱。相似的， 是和在时频点的瞬时互功率谱，对于其他分量，自功率谱和互功率谱的计算方法也是一致的。

（5）瞬时阻抗函数的分量根据公式（1）-（2）所得，（3）；

其他几个分量,,也同理得到。根据公式（3），[微软用户5]计算出各个时频点的阻抗张量元素，把大地电磁信号的时频功率谱图转变为时频阻抗谱图。

（6）时频阻抗谱图中，将具有相同瞬时频率(ω₀)的各时频点的响应函数值累加起来，除以时频点的个数，其均值则为频率值(ω₀)的响应函数估计值：

（4）；

然后利用、分别求得电阻率和相位，式中ij表示下标xx, xy, yx, yy；f表示频率。

从图5-图6上所知，用Hilbert 瞬时谱代替传统的Fourier谱，进行MT 响应参数的估算，将时频功率谱转化为时频阻抗谱、时频相位谱，可在时频域上来统计估算MT响应函数。这种方法并没有要求信号频率成分在整个时间段都平稳，时频响应参数是由真实存在的各时频点谱值计算得到的。相比傅立叶谱，减小了信号非平稳性对谱估计带来的误差，更接近大地电磁信号的真实谱，所以由此方法计算得到的响应函数更符合实际。从HHT瞬时谱上统计估算参量，比Fourier分析，更有利于实现稳健估计, 最小化了大地电磁信号非稳性带来的估算偏差。

从图7可知，本发明中的阻抗时频分析与传统Fourier方法相比：两种算法得到的曲线离散值的基本形态是一致的，而本发明处理结果比传统算法的数据更集中，曲线更平稳一些。就整个曲线形态来看，规律是明显的，估算的参量鲁棒性更好，曲线异常的情况得到缓解，缓解了个别频点鲁棒性差，异常突出的现象，估算的结果更具稳健性，资料的可解释性得到明显提高。

Claims (1)

1.一种大地电磁阻抗估计的时频分析方法，其特征在于，包括如下步骤：

（1）采用Hilbert-Huang变换中的经验模态分解技术，将大地电磁信号的4个分量H _x 、E _y 、 H _y 、E _x 分别进行分解，自适应地得到各阶模态函数，然后分别对相应分量的各阶模态函数进行Hilbert变换，计算得到相应分量上的各阶模态函数的瞬时幅度a _i (t)和瞬时频率ω _i (t)，综合相应分量上的各阶IMF的瞬时谱，得到大地电磁信号4个分量的时频谱H(ω, t )；

（2）根据时频谱H(ω, t )，求得任意时频点上信号的瞬时幅值，将得到的瞬时幅值代替传统傅立叶方法中的傅氏谱，带入到最小二次方程，求解时频点上的阻抗张量，数学模型如下：

（1），

（2）；

（3）瞬时阻抗函数的分量根据公式（1）-（2）可求得， （3），

类同公式（3），其他几个分量,,也同理得到；

计算出各个时频点的阻抗张量元素，把大地电磁信号的时频功率谱图转变为时频阻抗谱图；时频阻抗谱图中，将具有相同瞬时频率(ω₀)的各时频点的响应函数值累加起来，除以时频点的个数，其均值则为频率值(ω₀)的响应函数估计值：（4），

然后利用、分别求得电阻率和相位，式中ij表示下标xx, xy, yx, yy；f表示频率。