CSAMT数据静态校正方法与流程

本发明涉及一种针对csamt数据的静态校正方法。

背景技术：

可控源音频大地电磁法(controlledsourceaudio-frequencymagnetotellurics,简称csamt卡尼亚电阻率测深曲线，因此又称可控源音频大地电磁测深法，具有工作效率高、勘探深度范围大、空间分辨率高、地形影响小、高阻电屏蔽作用小等优点。csamt在信号采集过程中，地表或近地表往往存在局部导电不均匀体，当电流流过不均匀体的表面时，就会形成电荷累积效应，导致电场发生畸变，产生静态效应。在“视电阻率-频率”双对数坐标系中相邻各测点的视电阻率曲线沿着视电阻率轴上下移动，反映在视电阻率拟断面图上，会出现横向范围不大的陡立状密集带，掩盖地下由矿体或目标物引起的异常，对实际地质情况的判断造成误导。

对静态效应的校正方法可以分为:(1)曲线平移法做静校正；(2)数值分析法做静校正；(3)联合反演法做静校正；(4)利用相位换算资料做静校正；(5)利用磁场数据做静校正；(6)直接的二维以及多维反演方法；(7)小波滤波方法等。曲线平移法效果的好坏很大程度上取决于处理人员对数据的认知水平,对静态效应的判断和处理经验起着很关键的作用，容易误将那些并非静态效应引起的异常当作静态效应给“校正”了；对于某些地形背景简单的静态效应,尤其是地形影响产生的静态效应,可以采用数值模拟的方法来进行识别并作出校正，而在实际野外条件下这些往往是未知的，所以只能是作为参考研究用；联合反演法相当于重复做了两种地球物理方法，成本高，实用性不强；小波变换是一种变窗口变换，将地表不均匀性引起的静态效应与大构造异常在数学上进行直观的区分，有较高的分辨率和校正精度，然而小波方法存在的问题就是小波基函数的选择，小波基选择的不同以及分解层数的不同都会导致分解结果出现较大差异，进而会影响后面的静态校正效果。

hilbert-huang变换(hilbert-huangtransform，简称hht)方法是美国国家宇航局的n.e.huang等人提出的一种新的信号处理方法,对于非平稳、非线性信号有清晰的物理意义，能够得到信号的振幅一时间一频率分布特征，并且具有自适应性。针对传统方法在压制静态效应时所存在的问题，本发明提出了基于hht的csamt数据的静态校正。

技术实现要素：

本发明要解决的技术问题是提供一种csamt数据静态校正方法，该方法不存在基函数选择的问题，在实现静态校正的同时能够在压制静态效应的同时保留大构造异常信息。

为了解决上述技术问题，本发明的csamt数据静态校正方法包括以下步骤：

步骤一、采集可控源音频大地电磁离散点视电阻率数据，拟合出可控源音频大地电磁原始信号s(t)，对其进行经验模态分解(empiricalmodedecomposition，简称emd)得到n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t)，并将各本征模态函数作为imf分量；

步骤二、根据观测点与参考点的全部频点的视电阻率值之间的互相关系数rxy识别静态效应,如果存在静态效应，则按照步骤三对原始信号s(t)进行重构；

步骤三、根据步骤一得到的各imf分量对原始信号s(t)进行重构，得到静态校正后的可控源音频大地电磁重构信号s'(t)；

对原始信号s(t)进行重构的方法如下：

(1)采用基于高斯白噪声的自适应阈值作为重构阈值，设第i个imf分量ci(t)的重构阈值为ti：

其中，c表示阈值系数，取值范围为0.5-0.7；m表示选取的信号长度，即各imf分量的离散点个数，ei表示第i个imf分量ci(t)的能量强度，其计算公式为：

其中，是第一个imf分量c1(t)的方差；

(2)用h(ik)(t)表示第i个imf分量ci(t)在时间tk处的离散值，k＝1,2,...m,则第i个imf分量ci(t)对应的重构分信号为

则可控源音频大地电磁重构信号为s'(t)：

所述的阈值系数优选c＝0.7。

所述的步骤一中，对csamt的视电阻率数据进行经验模态分解得到多个本征模态函数分量的方法如下：

(一)、找出可控源音频大地电磁原始信号s(t)所有的极大值点和极小值点；

(二)、采用三次样条插值法分别将所有的极大值点和极小值点连接起来形成上包络线s1、下包络线s2：用原始信号s(t)各离散点视电阻率值ρ减去上包络线s1与下包络线s2的对应点均值,得到多个离散差值点并对其进行拟合得到第一个差值拟合曲线h11(t)；将第一个差值拟合曲线h11(t)所有的极大值点和极小值点连接起来形成差值信号上包络线h11、下包络线h12；考察第一个差值拟合曲线h11(t)是否满足imf的两个条件,即：极值点数目(包括极大值和极小值)和过零点数目相等或最多相差1个；任意点对应的差值信号上包络线h11上的离散点值与差值信号下包络线h12上的离散点值的平均值为0；如果满足则将第一个差值拟合曲线h11(t)作为第1个本征模态函数c1(t)，转到步骤三；否则用第一个差值拟合曲线h11(t)替代原始信号s(t)重复上述筛选过程，直至从原始信号s(t)中筛选出满足imf的两个条件的差值拟合曲线并将其作为第1个本征模态函数c1(t)；

(三)将剩下的信号r1(t)＝s(t)-c1(t)作为新的信号，按照步骤(一)、(二)，得到第二个本征模态函数c2(t)；依次类推，最后得到n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t)和一个不可分解的残余分量rn(t)：

n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t)作为imf分量。

所述步骤二中，识别静态效应方法如下:将观测点的全部频点的视电阻率值看成是一组数据系列xi，同样参考点的全部频点所对应的视电阻率值也是一组数据系列yi，将这两组数据进行相关匹配，求取两者之间的互相关系数rxy，

其中n为频点数量，为观测点的全部频点视电阻率值的平均值，参考点的全部频点视电阻率值的平均值；

若rxy≥0.85则认为存在静态效应。

本发明还可以包括下述步骤：

步骤四、对步骤一得到的各imf分量c1(t)、c2(t)、...ci(t)、...cn(t)作hilbert变换得到其瞬时频率和hilbert谱以及边际谱；针对任一imf分量ci(t)，其瞬时频率为ωi(τ)，原始信号s(t)的hilbert谱h(ω，τ)以及原始信号s(t)的边际谱h(ω)如下：

式中，yi(τ)是imf分量ci(t)经过希尔伯特变换后的信号，τ是希尔伯特变换后信号的时间变量，θi(τ)为相位，

h(ωi,τ)＝h(ωi(τ),τ)≡ai(τ)

其中h(ωi，τ)为imf分量ci(t)对应的hilbert谱中的一条谱线，ai(τ)为时变幅值；

本发明的有益效果是：(1)传统静态校正方法(例如：小波方法)都存在基函数选择困难的问题，而本发明提出的方法不存在这一问题，避免了这一限制条件。(2)本发明提出的csamt静态校正能够在压制静态效应的同时，保留原始数据中异常体的信息，使其不失真。(3)根据步骤四得到的hilbert谱和边际谱h(ω)与原始数据的hht谱和hht边际谱进行比对，可以直观地反映出静态效应的校正情况。时频谱中能够清楚显现信号各时段的不同频率，从边际谱中可以看出信号中存在的优势频率范围。

附图说明：

图1：本发明的csamt数据静态校正方法流程图。

图2：emd分解可控源音频大地电磁数据筛选imf分量的过程示意图。

图3：可控源音频大地电磁数据的hht谱。

图4：可控源音频大地电磁数据的hht边际谱。

图5：均匀半空间正演模型。

图6(a)、图6(b)、图6(c)分别为均匀半空间正演模型未经过静态校正处理、基于小波方法静态校正处理和采用本发明进行静态校正的视电阻率等值线图。

图7：含盖层的均匀半空间正演模型。

图8(a)、图8(b)、图8(c)分别为含盖层的均匀半空间正演模型未经过静态校正处理、基于小波方法静态校正处理和采用本发明进行静态校正的视电阻率等值线图。

具体实施方式：

下面结合附图与实施例对本发明作进一步详细说明。

如图1所示，本发明的csamt数据静态校正方法包括以下步骤：

步骤一、采集可控源音频大地电磁离散点视电阻率数据，拟合出可控源音频大地电磁原始信号s(t)，对其进行经验模态分解(empiricalmodedecomposition,简称emd)得到n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t)，并将各本征模态函数作为imf分量；如图2所示，具体分解步骤如下：

(一)、找出可控源音频大地电磁原始信号s(t)所有的极大值点和极小值点；

(二)、采用三次样条插值法分别将所有的极大值点和极小值点连接起来形成上包络线s1、下包络线s2；用各离散点视电阻率值ρ减去包络线s1、下包络线s2对应点的均值(图中m(t)为均值拟合曲线)得到多个离散差值点并对其进行拟合得到第一个差值拟合曲线h11(t)；采用三次样条插值法分别将第一个差值拟合曲线h11(t)所有的极大值点和极小值点连接起来形成差值信号上包络线h11、下包络线h12；考察第一个差值拟合曲线h11(t)是否满足imf的两个条件,即：1.极值点数目(包括极大值和极小值)和过零点数目相等或最多相差1个；2.任意点对应的差值信号上包络线h11上的点与差值信号下包络线h12上的点平均值为0；如果满足则将第一个差值拟合曲线h11(t)作为第1个本征模态函数c1(t)，转到步骤二；否则用第一个差值拟合曲线h11(t)替代原始信号s(t)重复上述筛选过程，直至从原始信号s(t)中筛选出满足imf的两个条件的差值拟合曲线并将其作为第1个本征模态函数c1(t)；c1(t)代表了原始信号s(t)中的高频成分；

(三)将剩下的信号r1(t)＝s(t)-c1(t)作为新的信号，按照步骤(一)、(二)继续筛选，得到第二个本征模态函数c2(t)；依次类推，最后得到n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t)和一个不可分解的残余分量rn(t)：残余分量rn(t)代表信号的均值或趋势项；这样原始信号s(t)就被分解成了n个imf分量(即n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t))和一个趋势项(即残余分量rn(t))之和：

经验模态分解得到的n个本征模态函数c1(t)、c2(t)、...ci(t)、...cn(t)按照频率从高到低的顺序排列。

步骤二、根据观测点与参考点的全部频点的视电阻率值之间的互相关系数rxy识别静态效应,具体方法如下：

首先，如果可控源音频大地电磁数据存在静态效应，则其等值线图中会在存在静态效应的观测点处出现明显的陡立带(如图6(a)和图8(a)中所示)。为了能够证实出现陡立带的观测点处确实存在静态效应，我们需要计算该观测点与其周围观测点的互相关系数。根据静态效应的特点，在双对数坐标中，受静态影响的观测点曲线与不受静态影响的参考点曲线形态不变，结合地下电性连续变化的特点，将观测点的全部频点的视电阻率值看成是一组数据系列xi，同样参考点的全部频点所对应的视电阻率值也是一组数据系列yi，将这两组数据进行相关匹配，求取两者之间的互相关系数rxy，认为如果曲线形态相同或者相近，且它们的互相关系数大，说明这是由静态效应引起的数据偏移，应当加以校正，反之，如果相关系数小，则判定这是由异常引起的反映地下电性的真实数据。一般选取临近观测点且有明显数值差异的视电阻率数据，或者是通过其他手段获得的区域背景视电阻率值数据作为参考点视电阻率数据。互相关系数rxy计算公式如下：

其中n为频点数量，为观测点的全部频点视电阻率值的平均值，参考点的全部频点视电阻率值的平均值；

若rxy≥0.85则认为存在静态效应，按照步骤三对原始信号s(t)进行重构；

步骤三、根据步骤一得到的各imf分量对原始信号s(t)进行重构，得到静态校正后的可控源音频大地电磁重构信号s'(t)；

对原始信号s(t)进行重构的方法如下：

(1)在重构阈值选择时，若重构阈值太小，则可能静态效应不能完全压制；若重构阈值太大，会将原始信号s(t)的大构造异常信息也压制掉，考虑到emd分解信号的特点，重构阈值采用基于高斯白噪声的自适应阈值。设第i个imf分量ci(t)的重构阈值为ti：

其中，c表示阈值系数，取值范围为0.5-0.7(经过大量实验分析，优选c＝0.7)；m表示选取的信号长度，即各imf分量选取的离散点个数，ei表示第i个imf分量ci(t)的能量强度，其计算公式为：

其中，是第一个imf分量c1(t)的方差；

(2)用h(ik)(t)表示第i个imf分量ci(t)在时间tk处的离散值，k＝1,2,...m,则第i个imf分量ci(t)对应的重构分信号为

则可控源音频大地电磁重构信号为s'(t)：

步骤四：hht时频谱又称hht能量谱(简称hht谱，如图3所示)，能够显示出可控源音频大地电磁信号能量随时频的具体分布规律。它相当于傅里叶谱，但比傅里叶谱具有更高的频率分辨率。hilbert边际谱(如图4所示)是通过对hilbert谱积分得到的，是一种信号谱估计方法，提供了对频率总的幅值(能量)分布，代表了整个信号在时间跨度上的幅值累积效应，突破了傅里叶变换只对平稳信号有效的不足之处。

对步骤一得到的各imf分量c1(t)、c2(t)、...ci(t)、...cn(t)作hilbert变换得到其瞬时频率和hilbert谱以及边际谱；对imf分量ci(t)作hilbert变换的数学表达式如下：

式中，yi(τ)是imf分量ci(t)经过希尔伯特变换后的信号，τ是希尔伯特变换后信号的时间变量，yi(τ)的解析信号如下：

zi(τ)＝ci(t)+jyi(τ)

j是解析信号的虚部单位。

由欧拉方程：e^jθ＝cosθ+jsinθ和yi(τ)的解析信号表达式可知：

令则可以得出：

zi(τ)＝ci(t)+jyi(τ)＝a(τ)·e^jθ(τ)，

其中ai(τ)为时变的幅值，且

同理，可推导出相位θ(τ)的表达式如下：

同样，瞬时频率ωi(τ)可定义为：

时变幅值ai(τ)的时频分布被定义为：

h(ωi，τ)＝h(ωi(τ)，τ)≡ai(τ)

h(ω，τ)是hilbert谱的其中一条谱线。

最后，把所有分量的hilbert谱累积到一起，即得原始信号的hilbert谱：

hilbert边际谱h(ω)同样可以被表示如下：

其中，h(ω)为单位频率内幅度分布的叠加。

具体实施例：

图5为建立的均匀半空间模型，模型厚度为900米，共计14个频点，分别为1至8192hz,测线距离ab为2000米，共计41个测点，测点间距为50米。在该均匀半空间中，设置了一个导电体也就是会引起静态效应的静态体。该静态体位于测点20附近，顶部埋深为0.5米，尺寸大小为50m*50m，电阻率大小为500ω·m。在该静态体的正下方还设置了一个异常体，该异常体的顶部埋深为600米，尺寸大小为200m*300m，电阻率值为1000ω·m，该模型的背景电阻率为200ω·m。

图6(a)为图5中所示模型未经过静态校正处理的视电阻率等值线图，如图中所示，在测点20至测点22之间存在陡立状的等值线，通过计算这几个测点处的互相关系数可以确定静态效应的存在。图6(b)为基于小波方法进行静态校正的视电阻率等值线图，从图中可以看出，静态效应得到了很好的抑制，然而，异常体的位置却发生了失真现象，与模型中所设置的异常体位置有所偏差。图6(c)为基于本发明进行静态校正的视电阻率等值线图，从图中可以明显看出，在静态效应得到了有效的抑制的同时，异常体的信息也准确地体现了出来。

图7为一表面覆有盖层的均匀半空间模型，模型盖层厚度为50米，盖层下面的均匀半空间厚度为900米，共计14个频点，分别为1至8192hz,测线距离为2000米，共计41个测点，测点间距为50米。在该模型中设置了两个静态体，分别位于测点10和测点16附近。静态体的顶部埋深均为50米，尺寸大小均为50m*50m，电阻率大小分别10ω·m和500ω·m。在测点17至测点25之间还设置了一个异常体，该异常体的顶部埋深为650米，尺寸大小为200m*300m，电阻率值为100ω·m，该模型的背景电阻率为50ω·m。

图8(a)为图7中所示模型未经过静态校正处理的视电阻率等值线图，如图中所示，在测点10附近存在陡立状的等值线，通过计算该测点处的互相关系数可以确定静态效应的存在。而在测点16附近，盖层与下方的高阻体连在了一起，无法进行信息的有效判别，严重阻碍了下一步的数据解释。图8(b)为基于小波方法进行静态校正的视电阻率等值线图，从图中可以看出，静态效应得到了很好的抑制，然而，地下异常体的信息也同样被抑制掉了。图8(c)为基于本发明进行静态校正的视电阻率等值线图，从图中可以明显看出，在静态效应得到了有效的抑制的同时，异常体的信息也准确地体现了出来。