一种STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法

摘要

本发明属于一种地球物理数据处理方法，具体涉及一种STRATAGEM大地电磁测量系统(EH4测量系统)阻抗张量的稳健(Robust)估算方法。该方法把统计学中的最大似然估计(M估计)引入到STRATAGEM系统的独立观测中，使多次独立观测的数据均参与最终阻抗运算，但是每次观测参与计算的权值(重)不同。从而避免了在最小二乘估算中某次独立观测影响整体计算结果的可能性，提高了估算的可靠性与稳健性。

说明书

技术领域

本发明属于一种地球物理数据处理方法，具体涉及一种STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法。

背景技术

STRATAGEM(EH4)大地电磁测量系统是上世纪末期引入我国的先进大地电磁(音频)测量设备。该设备目前在国内有一百余套，在寻找地下水、资源勘探中有广泛的用途。但是该系统所配的数据处理软件(IMAGEM)仍然采用较原始的大地电磁阻抗估算方法--最小二乘方法。由于个别时段的干扰信号影响最小二乘估算结果，不同时段组合运算所得到的阻抗曲线可能有较大差异，从而导致两个方面结果：

(1)阻抗运算结果飞点数量多；

(2)阻抗运算结果不稳定。

而调整系统阻抗估算流程，并融合稳健统计方法所计算出的阻抗曲线则可以有效地抑制这种情况的发生，为大地电磁资料的解释提供较为可靠的依据，是长时间以来广大用户迫切希望解决的实际应用问题。

发明内容

本发明的目的在于针对现有大地电磁阻抗估算方法的缺陷，提供一种STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法，以提高估算的可靠性。

为实现上述目的，本发明的技术方案如下：一种STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法，包括如下步骤：

(1)读取原始时间序列数据，获取单个测点的多次独立观测时间序列原始结果；

(2)用FFT变换方法，实现由时间域到频率域的转化，得到每个频点多次独立观测的频谱值；

(3)采用最大似然估计的方法，对同一频点的多个独立观测结果给出不同的权重；

(4)进行原始频点谱的叠加，实现由FFT变换得到的原始频点到阻抗频点的叠加；

(5)由每个阻抗频点叠加频谱值计算不同测量场的互功率谱；

(6)通过融合标定使测量场转化为地球物理单位意义的电磁场数值，从而计算出阻抗值。

进一步，如上所述的STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法，步骤(3)中所述的最大似然估计的方法包括如下步骤：

①用最小二乘法计算同一频点的多个独立观测结果的初始阻抗；

②计算初始尺度估计参数；

③计算权函数w(r)，并计算迭代更新电场值数据；

④采用最小二乘法计算同一频点的多个独立观测结果的新的阻抗张量与误差尺度；

⑤判断残差是否达到设定值以及是否达到设定的迭代次数；如果否，则返回步骤③，如果是则进行后续的阻抗计算。

更进一步，如上所述的STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法，步骤⑤中残差r的表达式采用如下形式：

$r=\frac{E_i-H_i\hat{T}}{\sigma }$

其中，E_i为实测电场，H_i为实测磁场，为阻抗张量，σ为尺度估计参数。

更进一步，如上所述的STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法，步骤③中以权函数w(r)的权，迭代后的电场值表示为：

$E_{i1}={\hat{E}}_{i0}+w(r_{i0}/{\hat{\sigma }}_0)r_{i0}$

其中，预测电场值${\hat{E}}_{i0}={H^T}_i{\hat{T}}_0,$r_iO为预测电场值及其与实测电场值之间的残差，分别为阻抗张量与尺度估计参数σ初始值由最小二乘方法对N次独立观测所得的结果，H^T_i为估算阻抗张量。

本发明的有益效果如下：本发明把统计学中的最大似然估计(M估计)引入到STRATAGEM系统的独立观测中，使多次独立观测的数据均参与最终阻抗运算，但是每次观测参与计算的权值(重)不同。从而避免了在最小二乘估算中某次独立观测影响整体计算结果的可能性，提高了估算的可靠性。

附图说明

图1为本发明的总体方法流程图；

图2为最大似然估计(M估计)的方法流程图；

图3为对存在干扰的STRATAGEM系统的时序数据采用本发明的处理方法得到的阻抗曲线；

图4为STRATAGEM系统配置的IMAGEM处理方法得到的阻抗曲线。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明的主要创新点是把统计学中的最大似然估计(M估计)引入到STRATAGEM系统的独立观测中，使多次独立观测的数据均参与最终阻抗运算，其它有关FFT变换(傅里叶变换)，原始频点谱的叠加，互功率谱的计算以及阻抗值的计算均为本领域的公知技术。

如图1所示，一种STRATAGEM大地电磁测量系统阻抗张量的稳健估算方法，包括如下步骤：

(1)读取原始时间序列数据，获取单个测点的多次独立观测时间序列原始结果；

(2)用FFT变换方法，实现由时间域到频率域的转化，得到每个频点多次独立观测的频谱值；

(3)采用最大似然估计(M估计)的方法，对同一频点的多个独立观测结果给出不同的权重；

(4)进行原始频点谱的叠加，实现由FFT变换得到的原始频点到阻抗频点的叠加；

(5)由每个阻抗频点叠加频谱值计算不同测量场(电、磁场)的互功率谱；

(6)通过融合标定使测量场转化为地球物理单位意义的电磁场数值，从而计算出阻抗值。

如图2所示，最大似然估计(M估计)的方法包括如下步骤：

①用LS方法(最小二乘法)计算同一频点的多个独立观测结果的初始阻抗；

②计算初始尺度估计参数；

③计算权函数w(r)，并计算迭代更新电场值数据；

④采用LS方法(最小二乘法)计算同一频点的多个独立观测结果的新的阻抗张量与误差尺度；

⑤判断残差是否达到设定值以及是否达到设定的迭代次数；如果否，则返回步骤③，如果是则进行后续的阻抗计算。

在最大似然估计(M估计)过程中主要函数选取如下：

●损失函数

损失函数采用Huber(1981)的分段函数，残差r的表达式如下式所示：

$r=\frac{E_i-H_i\hat{T}}{\sigma }$

其中，E_i为实测电场，H_i为实测磁场，为阻抗张量，σ为尺度估计参数。

阻抗张量与σ尺度估计参数初始值由最小二乘方法对N次独立观测给出，分别为则，对于第i次观测值而言，其预测电场值及其与实测电场值之间的残差r_i0为：

$r_{i0}=\mathrm{Ei}-{\hat{E}}_{i0}$

预测电场值${\hat{E}}_{i0}={H^T}_i{\hat{T}}_0,$其中，H^T_i为估算阻抗张量。

权函数w(r)形式为损失函数的导数与残差r的比值。

则以权函数w(r)的权，迭代后的电场值可以表示为：

$E_{i1}={\hat{E}}_{i0}+w(r_{i0}/{\hat{\sigma }}_0)r_{i0}$

●尺度参数选取

尺度参数的选取与残差r的计算紧密联系，选取残差的均方根值计算方法，具体公式如下：

$\sigma ={\left(\frac{1}{n-p}\underset{i=1}{\overset{n}{\Sigma }}{r}_i^2\right)}^{\frac{1}{2}}·{\left(\sqrt{1-m_{\mathrm{jj}}}\right)}^{-1}$

其中，m_jj磁场帽子矩阵对角元素，n表示独立观测次数，p表示自由度。

图3为对存在干扰的STRATAGEM系统的时序数据采用本发明的处理方法得到的阻抗曲线，图4为STRATAGEM系统配置的IMAGEM处理方法得到的阻抗曲线。从两幅图的比较可以看出本发明提供的估算方法更稳健，估算的可靠性更高。