CN111273367B - 一种大地电磁阻抗的估算方法 - Google Patents

Abstract

本发明提供一种大地电磁阻抗的估算方法，包括将从观测站所测的时间域大地电磁信号进行傅里叶变换且变换到频率域、利用最小二乘法估算初始阻抗并得出初始残差、区分出各类噪声、得出对应磁场归一化的权重值、获取估算出的大地电磁阻抗、对衡量噪声是否接近高斯噪声的参数的数值进行判断、输出阻抗等步骤。本发明基于噪声本身加权的方式进行阻抗估计，能很好地利用噪声与信号之间的关系，并且利用这种相关性进行加权归一化，使随信号强度增加而增加的高斯噪声在加权归一化后还原成了高斯噪声；本发明方法与传统的最小二乘法和稳健方法对比，这种基于噪声本身加权的方法真正意义上解决了随信号强度增加而增加的高斯噪声的问题。

Description

一种大地电磁阻抗的估算方法

技术领域

本发明涉及大地电磁信号处理领域，具体涉及一种大地电磁阻抗的估算方法。

背景技术

大地电磁法在矿产勘查、地质勘查、深部构造、地震预报、地热等方面中得到了广泛的应用。大地电磁法主要利用天然场源研究地下电阻率从浅到深的变化，勘探深度大，工作便利。但利用天然场源也使大地电磁观测到的电磁场信号非常微弱，很容易受到各种噪声的干扰，极大地影响了大地电磁数据的质量和应用效果。因此，如何压制干扰，获取无偏的大地电磁阻抗估计，一直是大地电磁领域讨论和研究的热点。

在大地电磁信号处理中，主要采用以下步骤：首先，对大地电磁信号进行傅里叶变换，将时间域的大地电磁信号转换到频率域；其次，在频率域中利用平面波电磁场线性关系实现阻抗估计。目前，常用的大地电磁阻抗估计的方法大都基于线性回归，在线性回归中有最小二乘法，对于独立均一高斯分布噪声，最小二乘法可以获得无偏估计。但是，在实际测量获得的大地电磁数据中，只有极少的大地电磁数据的噪声满足独立均一高斯分布噪声这个条件，大部分的大地电磁信号存在离群值，因此，在最小二乘法估计方法的基础上提出robust方法，用来剔除一些离群值，减少离群值对阻抗估计的影响。但是研究表明robust方法估计的残差不服从高斯分布，其分布具有厚尾，即反映剔除离群数据后数据噪声不完全是高斯分布，robust估计的效率会受到影响。如何减少这种影响，获得最优的阻抗估计是大地电磁数据处理中亟需解决的问题。

发明内容

本发明的目的在于提供一种大地电磁阻抗的估算方法，直接利用噪声本身进行加权，用于解决实测大地电磁中数据中存在的问题(即同一频带内由于存在不同噪声源或噪声随信号强度增加而增加所带来的非高斯噪声的问题)，最终很好地获得了无偏的大地电磁阻抗估计值。具体技术方案如下：

一种大地电磁阻抗的估算方法，包括以下步骤：

步骤一、将从观测站所测的时间域大地电磁信号进行傅里叶变换，变换到频率域；

步骤二、利用最小二乘法估算初始阻抗并得出初始残差；

步骤三、对残差能量和输入信号能量进行聚类，区分出n类噪声，n取大于等于1的自然数；

步骤四、对每一类噪声进行线性拟合，得到线性函数f_i(x)，i＝1,2,...,n；将磁场能量带入线性关系，得出下式所示对应磁场归一化的权重值W：

W＝diag[1/f₁(B₁),1/f₂(B₂),...,1/f_n(B_n)]；

其中：diag表示对角矩阵，它的对角项都为w(r_i)；B₁、B₂、…、B_n为每一类噪声的磁场值；f_i(B_i)每一类噪声关于磁场的线性函数，i＝1,2,...,n；

步骤五、结合步骤四所得权重值W，通过下式得到估算出的大地电磁阻抗

并衡量噪声是否接近高斯噪声的参数的数值；

其中：E为傅里叶变换后的电场信号，B为傅里叶变换后的磁场信号，B^H为B的共轭转置矩阵；

步骤六、进行判断，若迭代次数满足要求或参数的数值满足要求，则进入下一步，否则，取n＝n+1，返回步骤三；

步骤七、将参数取值最优的一次的阻抗进行输出(找出α最大的那次迭代计算出的阻抗记为最终阻抗)。

以上技术方案中优选的，所述步骤一具体是：将实地采集的大地电磁时间序列依次通过预白、重叠、加窗和快速傅里叶变换，将时间域序列转换为频率域；

所述预白采用一阶差分或自回归模型；

所述重叠具体是：将信号分成多段，但相邻段与段之间是有重叠部分；

所述加窗采用Hamming窗口或Blackman窗口；

加上噪声后电场和磁场的关系采用下式表示：

E＝ZB+ε；

其中：E为傅里叶变换后的电场信号，B为傅里叶变换后的磁场信号，ε为噪声。

以上技术方案中优选的，所述步骤二具体是：利用下式所示的最小二乘法公式计算得到初始的大地电磁阻抗：

Z＝EB^H(BB^H)^-1；

采用下式计算残差r：

r＝E-ZB。

以上技术方案中优选的，所述步骤三中：采用kmeans方法进行聚类，根据欧氏距离来寻找聚类的中心点；n取大于等于1且小于等于20的自然数。

以上技术方案中优选的，所述步骤四中：线性拟合采用最小二乘法回归或稳健回归。

以上技术方案中优选的，所述步骤五中：用于衡量噪声是否接近高斯噪声的参数，是稳定分布中的α参数，α参数的取值越接近2，说明噪声越接近高斯分布。

应用本发明的技术方案，效果是：本发明基于噪声本身加权的方式进行阻抗估计，能很好地利用噪声与信号之间的关系，并且利用这种相关性进行加权归一化，使随信号强度增加而增加的高斯噪声在加权归一化后还原成了高斯噪声；本发明方法与传统的最小二乘法和稳健(robust)方法对比，这种基于噪声本身加权的方法真正意义上解决了随信号强度增加而增加的高斯噪声的问题。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明优选实施例的流程示意图；

图2是噪声α值随着噪声种类数的增加变化图；

图3是西藏一测站大地电磁数据进行噪声能量和磁道能量聚类和拟合示意图，其中：不同形状代表不同种类的噪声，黑线代表拟合出的权重值；

图4是西藏一测站大地电磁数据一频点的qq和pp图，其中：(a)为未经归一化robust方法处理后的qq图，(b)为经过基于噪声本身加权归一化后的qq图，(c)为未经归一化robust方法处理后的pp图，(d)为经过基于噪声本身加权归一化后的pp图；其中：+为残差数据点，虚线代表的是pp中的95％的置信区间；

图5是计算得出的视电阻率与真实电阻率的误差及相位与真实相位的误差对比图，其中：(a)是计算得出的视电阻率与真实电阻率的误差；(b)是计算得出的相位与真实相位的误差。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。

实施例：

一种大地电磁阻抗的估算方法，包括步骤如附图1，详情包括：

第一步、将时间域大地电磁信号进行傅里叶变换，变换到频率域，具体是：

将实地采集的大地电磁时间序列通过预白、重叠、加窗和快速傅里叶变换，将时间域序列转换为频率域中，其中：预白可选用的方式有一阶差分和自回归模型，此处优选一阶差分；重叠是指将信号分成多段，但相邻段与段之间是有重叠部分(具体重叠多少根据实际操作确定，原则上数据长度越长，重叠也越长)；加窗是为了防止频谱泄露，一般加的窗口有Hamming窗口和Blackman窗口，此处优选Hamming窗口，如表达式1)：

其中：ω(n)为Hamming系数，N为数据总数，n＝1,2,…,N。

快速傅里叶变换是离散傅里叶变换的一种快速的编程形式，其本质还是离散傅里叶变换，假设时间域的大地电磁信号为x(t)，变换后的频率域大地电磁信号为x(f)，通过离散傅里叶变换公式表达式2)进行变换：

其中：j为复数标志，f为频率。

在对大地电磁时间域信号进行傅里叶变换后，对于一个特定的频段，将傅里叶变换后的电场信号E和磁场信号B采用表达式3)所示的向量进行表示：

其中：Exi为电场x方向第i个数据，Eyi为电场y方向第i个数据，Bxi为磁场第x方向第i个数据，Byi为磁场y方向第i个数据；

加上噪声后电场和磁场的关系采用表达式4)表示：

E＝ZB+ε 4)；

其中：E为傅里叶变换后的电场信号，B为傅里叶变换后的磁场信号，ε为噪声。

第二步、利用最小二乘法估算初始阻抗并得出初始残差，具体是：

利用最小二乘公式表达式5)计算得到初始阻抗(初始的大地电磁阻抗)：

Z＝EB^H(BB^H)^-1 5)；

其中：Z为初始的大地电磁阻抗，E为傅里叶变换后的电场信号，B为傅里叶变换后的磁场信号，H上标表示共轭转置，即B^H为B的共轭转置矩阵；

利用表达式6)计算残差r：

r＝E-ZB 6)。

第三步、对残差能量和输入信号(磁道)能量进行聚类，区分出n类噪声，n取大于等于1的自然数(优选n取大于等于1且小于等于10的自然数)，具体是：

此处的聚类，使用kmeans方法进行聚类，考虑到同一噪声来源之间的几何距离一般都是最近的，一般使用一范数根据欧氏距离来寻找聚类的中心点。分类的目的也就是区分出不同噪声来源，对不同噪声进行归一化；同时剔除数据量很少的种类中的数据，这种数据量很少的数据一般包含信号信息较少，故直接进行剔除。

第四步、先对每一类噪声进行线性拟合，得到线性函数f₁(x)、f₂(x)、...、f_n(x)，此处拟合线性关系时可采用最小二乘法回归和robust回归，此处优选最小二乘法回归；再将磁场能量带入线性关系，得出对应磁场归一化的权重值W，如表达式7)：

W＝diag[1/f₁(B₁),1/f₂(B₂),...,1/f_n(B_n)] 7)；

其中：diag表示对角矩阵，它的对角项都为w(r_i)(可参照现有技术中robust方法获得)；B₁、B₂、…、B_n为每一类噪声的磁场值；f_i(B_i)每一类噪声关于磁场的线性函数，i＝1,2,...,n；

权重函数可根据实际情况选择，如选用huber权重函数、bisquare权重函数等，使噪声还原到最为接近高斯噪声。

第五步、结合步骤四所得权重值W(即权重矩阵)，通过表达式8)得到估算出的大地电磁阻抗

并衡量噪声是否接近高斯噪声的参数的数值；

此处，为了判断噪声是否接近高斯噪声，选用统计学中的stable分布的α参数，噪声的α值越接近2说明噪声越接近高斯分布。

第六步、进行判断，若迭代次数满足要求或参数(α参数)的数值满足要求，则进入下一步，否则，取n＝n+1，返回第三步；

第七步、将参数(α参数)取值最优的一次的阻抗进行输出。

应用本实施例的技术方案，效果是：

本发明基于噪声本身加权的方式进行阻抗估计，这种方法很好地利用了噪声与信号之间的关系，并且利用这种相关性进行加权归一化，使随信号强度增加而增加的高斯噪声在加权归一化后还原成了高斯噪声。

本发明适用性强，具体是：无论是在高噪的环境中，还是在低噪的环境中，采集的大地电磁数据都能进行处理。一般来说，低噪环境n值较小时(一般取5)就可以取得很好的结果；而高噪环境n可以取得大一些(但一般不大于10，利于节约计算时间)。图2展示西藏观测站中一个大地电磁数据处理，可以看出基本随着噪声种类n增加，噪声的α值增大，意味着更加接近高斯噪声，可以得到更优的大地电磁阻抗估计。图2中n＝2时就已经超过1.8，满足条件，可以输出数据。图3展现了十类噪声(n＝10)时的聚类结果。

与传统的最小二乘法和robust方法对比，这种基于噪声本身加权的方法真正意义上解决了随信号强度增加而增加的高斯噪声的问题。该估计的残差，通过q-q图和p-p图进行后验概率验证，如果qq图和pp图中，残差对应高斯分布是一条直线，那么残差是服从高斯分布的，也就满足最小二乘法的前提条件，即可以得到无偏的阻抗估计。其中qq图中验证的为残差幅值的分布，如果残差实部虚部都满足高斯分布，那么其幅值满足瑞利分布；而pp图直接验证的是残差的实部是否满足高斯分布。从图4中西藏观测站的大地电磁数据的残差的qq图和pp图可以看出，经过robust方法处理后的大地电磁数据残差是存在厚尾现象的(详见(a)和(c))，明显远离高斯分布，因此，获得的阻抗估计是有偏的。而经过本发明基于噪声本身加权归一化后的大地电磁数据的残差在qq图和pp图中明显是一条直线((b)和(d))，解决了残差厚尾现象，很好地满足了高斯分布，这意味着本发明方法处理后的数据满足了最小二乘的前提条件，获得的阻抗是无偏的。

此外，为了进一步对比本发明方法的有效性，通过理论合成的大地电磁数据，加上与信号相关的噪声进行测试，通过估算出的大地电磁阻抗，通过简单换算成视电阻率和相位比较。详见图5，(a)为计算得出的视电阻率与真实电阻率的误差，(b)为计算得出的相位与真实相位的误差，从图5中(a)和(b)可知：理论数据的真实阻抗是已知的，因此，可以对比各种方法估算出的视电阻率与真实视电阻率的误差以及相位与真实相位的误差来量化每个方法的有效性。从测试结果可以看出：在不同频段的视电阻率估计上，绝大部分情况中基于噪声本身加权归一化的方法的误差是最小的，在有些频段视电阻率和相位估计的误差甚至远远小于传统的robust和最小二乘法，其他频段的视电阻率和相位也与传统robust方法基本相似，总体上比robust方法效果要更好。

由上可知，本发明方法不仅仅在理论上解决了存在不同噪声源或噪声随信号强度增加而增加带来的问题，同时，也在实际运用中能估算出更精确的阻抗，在大地电磁数据处理中具有十分重要的意义。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种大地电磁阻抗的估算方法，其特征在于，包括以下步骤：

步骤一、将从观测站所测的时间域大地电磁信号进行傅里叶变换，变换到频率域；

步骤二、利用最小二乘法估算初始阻抗并得出初始残差；

步骤三、对残差能量和输入信号能量进行聚类，区分出n类噪声，n取大于等于1的自然数；

步骤四、对每一类噪声进行线性拟合，得到线性函数f_i(x)，i＝1,2,...,n；将磁场能量带入线性关系，通过下式得出对应磁场归一化的权重值W：

W＝diag[1/f₁(B₁),1/f₂(B₂),...,1/f_n(B_n)]；

其中：diag表示对角矩阵，它的对角项都为w(r_i)；B₁、B₂、…、B_n为每一类噪声的磁场值；f_i(B_i)每一类噪声关于磁场的线性函数，i＝1,2,...,n；

步骤五、结合步骤四所得权重值W，通过下式得到估算出的大地电磁阻抗

并衡量噪声是否接近高斯噪声的参数的数值；

其中：E为傅里叶变换后的电场信号，B为傅里叶变换后的磁场信号，B^H为B的共轭转置矩阵；

步骤六、进行判断，若迭代次数满足要求或参数的数值满足要求，则进入下一步，否则，取n＝n+1，返回步骤三；

步骤七、将参数取值最优的一次的阻抗进行输出。

2.根据权利要求1所述的大地电磁阻抗的估算方法，其特征在于，所述步骤一具体是：将实地采集的大地电磁时间序列依次通过预白、重叠、加窗和快速傅里叶变换，将时间域序列转换为频率域；

所述预白采用一阶差分或自回归模型；

所述重叠具体是：将信号分成多段，但相邻段与段之间是有重叠部分；

所述加窗采用Hamming窗口或Blackman窗口；

加上噪声后电场和磁场的关系采用下式表示：

E＝ZB+ε；

其中：E为傅里叶变换后的电场信号，Β为傅里叶变换后的磁场信号，ε为噪声。

3.根据权利要求2所述的大地电磁阻抗的估算方法，其特征在于，所述步骤二具体是：利用下式所示的最小二乘法公式计算得到初始的大地电磁阻抗Z：

Z＝EB^H(BB^H)^-1；

采用下式计算残差r：

r＝E-ZB。

4.根据权利要求1所述的大地电磁阻抗的估算方法，其特征在于，所述步骤三中：采用kmeans方法进行聚类，根据欧氏距离来寻找聚类的中心点；n取大于等于1且小于等于20的自然数。

5.根据权利要求1所述的大地电磁阻抗的估算方法，其特征在于，所述步骤四中：线性拟合采用最小二乘法回归或稳健回归。

6.根据权利要求1所述的大地电磁阻抗的估算方法，其特征在于，所述步骤五中：用于衡量噪声是否接近高斯噪声的参数，是稳定分布中的α参数，α参数的取值越接近2，说明噪声越接近高斯分布。

Images