基于杜宾-沃森统计量的大地电磁一维反演方法❋

2023-02-21李予国段双敏

中国海洋大学学报（自然科学版） 2023年2期

杨雯，李予国,2❋❋，段双敏，韩波，罗鸣,2

(1. 中国海洋大学海洋地球科学学院，山东青岛 266100； 2. 青岛海洋科学与技术试点国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266237； 3. 中国地质大学(武汉) 地质调查研究院，湖北武汉 430074)

大地电磁法(Magnetotelluric,MT)是利用天然交变电磁信号研究地球内部电性结构的一种地球物理方法，现已广泛应用于油气资源勘探和俯冲带、地壳与地幔电性结构研究中[1-3]。大地电磁反演是获得地下介质电阻率信息的重要方法。然而，MT反演存在多解性，如何减少MT反演的多解性一直是重要的研究方向。自从Karl Pearson提出卡方检验以来，在地球物理反演中便一直使用均方根拟合差(RMS misfit)衡量反演模型响应与实测数据的拟合程度[4]。

1980年代末以来，为了改善大地电磁反演问题的不适定性，吉洪诺夫(Tikhonov)正则化方法开始被应用于MT反演中，并取得了迅速发展[5-7]。目前，常用的大地电磁反演方法包括OCCAM反演方法、非线性共轭梯度法(NLCG)、拟牛顿法和高斯-牛顿法等[7-10]。这些反演方法关于反演模型响应与观测数据拟合程度的描述大都依赖于目标函数中的均方根拟合差，却没有考虑反演模型响应与观测数据之差(常称为残差)的统计特征，这可能导致反演模型响应形态与观测数据形态出现偏离较大的情况，尤其是当观测数据存在较大误差时，甚至可能得到错误的反演模型[11]。于是用于衡量实测数据与反演模型响应之间关系的其它统计量开始被应用于大地电磁反演中。Smith 和 Booker将斯皮尔曼(Spearman)相关系数应用于大地电磁(MT)反演结果评价中[12]。Jones利用DW统计量衡量观测数据与模型响应之间残差的自相关性[11]，并将其应用于加拿大Great Slave Lake地区大地电磁测深数据反演中[13]。将统计学度量应用于MT反演中，一方面可以评价反演结果的合理性，另一方面可以为反演算法的改进提供依据，以便进一步优化反演方法。但是目前关于这方面的研究还不多见。本文将DW统计量应用于高斯-牛顿大地电磁一维反演中，通过修改目标函数，实现了基于DW统计量的大地电磁一维反演方法(以下简称为DW反演方法)。通过模型算例分析了该方法的有效性，并与传统的高斯-牛顿反演结果进行了对比。

1 反演方法

1.1 DW统计量

杜宾-沃森(Durbin-Watson，简称DW)统计方法是用于检验回归分析中残差一阶自相关性的一种方法，

该方法常用于检验计量经济学中的自相关问题[14-15]。一阶自相关DW统计量的定义式为[11,16-17]：

(1)

图1 序列自相关性与DW统计量的关系

表1 部分DW统计量临界值分布表(k=2)

通常地，DW统计量与自相关系数R具有如下关系[16,18]：

DW≈2(1-R)。

(2)

式中，R为一阶自相关系数，其表达式为：

(3)

一阶自相关系数R的取值范围为(-1，1)。通常认为，当|R|<0.2时，残差序列的自相关性极弱。

1.2 DW反演方法

大地电磁反演问题实质上是目标函数的最优化问题。基于吉洪诺夫正则化反演思想，反演目标函数中通常包含数据拟合差项和正则化项[19]：

φ(m)=φd(m)+λ1φm(m)。

(4)

式中：m=[m1,m2,···,mM]T为模型参数向量，由反演模型中各地层的电阻率和厚度构成；φd(m)为观测数据与理论模型响应的拟合差，可表示为：

(5)

(6)

式中：Wm为拉普拉斯算子的差分近似；mref为参考模型。

在传统的反演方法中，反演目标函数仅考虑了反演模型响应与观测数据的拟合差，但没有考虑其残差的统计学特征。本文将DW统计量加入反演目标函数中，通过反演倾向于寻找在拟合差和残差自相关性两个方面都可以接受的模型。经过修改后的反演目标函数具有如下形式[11]：

φ(m)=φd(m)+λ1φm(m)+λ2(DW-2)2。

(7)

上式中等号右端第三项是将DW统计量归一化为零。在反演中，正则化因子λ1从一个较大的初值开始，随迭代次数的增加逐渐减小；而与此相反，对于权重系数λ2，则先给定一个较小的初值，在迭代过程中缓慢递增。λ1递减和λ2递增的系数通常位于(1，2)之间。采用这种策略的基本思想是，在迭代反演的前期，注重数据拟合和模型约束，而在反演的后期通过增加DW项的权重，减少实测数据与反演模型响应之间残差的自相关性。

本文采用高斯-牛顿最优化方法求解目标函数式(7)的极小值。迭代反演过程中，模型更新量为：

(8)

其中gk和Hk分别为第k次迭代时目标函数φ(m)的梯度向量和Hessian矩阵，其表达式为：

(9)

(10)

在传统的高斯-牛顿反演方法中，gk和Hk分别仅由式(9)和式(10)右端的前两项构成，而除该两项以外的其它项反映了DW约束项对梯度向量和Hessian矩阵的贡献。其中：

(11)

Hdw=2QJTKTKJ-2PQ2LLT+8Q2JTKTKe(LE)T+
8PQ3LE(LE)T。

(12)

1.3 合成数据添加误差

大地电磁阻抗Z可以写成下列指数形式

Z=reiφ。

(13)

由式(13)可以得到阻抗微小变化dZ的表达式

dZ=dreiφ+reiφ·idφ。

(14)

式(14)除以式(13)，得

(15)

上式表明，大地电磁阻抗Z的相对变化相当于其幅值的相对变化和相位的绝对变化。大多数时间序列分析表明，大地电磁阻抗的实部和虚部通常具有相等的误差。于是，有

(16)

这意味着，对于大地电磁阻抗Z，其相位的绝对误差与幅值的相对误差相对应。也就是说，如果给大地电磁阻抗幅值添加相对误差er, 则对应的应该给阻抗相位添加er弧度或er·180/π角度的绝对误差。

(17)

大地电磁理论模型合成数据反演中，通常会给正演合成数据加入一定的随机噪声，并将其结果作为反演数据。然而，若直接给视电阻率添加相对误差，可能会使得所添加的误差具有自相关性，以2.1节将要讨论的三层高阻薄层模型为例，对此进行说明。

首先，在频率范围10-3～103Hz之间选取21个对数等间隔的频点，正演计算得到三层高阻薄层模型21个频点的视电阻率和阻抗相位。然后，在视电阻率数据中加入1%的随机相对误差噪声，相应地在阻抗相位数据中加入0.286 5°的随机绝对误差。用这样的方法得到的视电阻率误差和相位误差随周期的变化关系如图2所示。

图2 视电阻率误差(a)和阻抗相位误差(b)分布特征

从图2中可以看到，阻抗相位误差在周期范围(10-3～103s)内是随机分布的。阻抗相位误差的DW统计量DWpha=1.74和自相关系数Rpha=0.036，这表明阻抗相位误差不具有自相关性。长周期视电阻率的误差是在零线附近上下随机分布的，而短周期视电阻率的误差并非是随机分布的。视电阻率误差的DW统计量DWrho为 0.81，自相关系数Rrho为0.568，这表明视电阻率误差具有一定的自相关性。图3分别绘出了视电阻率误差序列和相位误差序列与其滞后一阶误差序列的散点图，从图3中也可以看到，视电阻率误差没有呈现出随机分布的状态，但相位误差呈现出了明显的随机分布特征。基于DW统计量的大地电磁反演方法需要满足噪声误差随机性的条件，因此需要对视电阻率进行转换。

图3 视电阻率误差序列(a)和阻抗相位误差序列(b)与其滞后一阶误差序列散点图

在大地电磁反演中，通常会将视电阻率的对数值作为反演数据。若采用给视电阻率对数值添加绝对误差的方式，则所加误差便不具有自相关性。依据前面的讨论，如果给大地电磁阻抗添加相对误差er, 则应给视电阻率加入2er的相对误差，这相当于给视电阻率对数添加绝对误差log10(1+2er)。其原因解释如下。

假如给视电阻率ρa添加相对误差2er，则绝对误差Δρa=2ρa·er。若视电阻率取对数，y=log10ρa，则其绝对误差为：

(18)

还是以高阻薄层模型为例，在视电阻率对数数据中加入log10(1+1%)的随机绝对误差。图4给出了视电阻率对数的绝对误差随周期的变化关系以及散点图。由图4可见，视电阻率对数的误差呈现出了随机分布状态。视电阻率对数误差的DW统计量DWrho= 1.74和自相关系数Rrho=0.036，这表明各个频点视电阻率对数的绝对误差之间不存在自相关性。

图4 视电阻率对数值的绝对误差分布情况(a)和视电阻率对数值的误差与其滞后一阶误差序列散点图(b)

考虑到理想情况下，各个频点大地电磁数据的误差之间应各自独立，不存在自相关。因此，在大地电磁反演中应将视电阻率对数值和阻抗相位作为反演数据，在添加误差时应加入相应的绝对误差。

2 反演算例

2.1 高阻薄层模型

首先考虑一个含有高阻薄层的三层模型(如图5中黑色虚线)。该模型是根据文献[11,20]建立的页岩油气高阻薄层模型。第一层的电阻率和厚度分别为30 Ω·m和120 m；中间层为高阻薄层，其电阻率和厚度分别为120 Ω·m和10 m，高阻层下方是电阻率为2.5 Ω·m的均匀半空间。假设21个频点对数等间隔地均匀分布在频率范围10-3～103Hz之间。在各个频点的视电阻率对数和阻抗相位数据中，分别加入0.004 3(log10(1+0.01)=0.004 3，相当于给视电阻率添加1%的相对误差)和0.005弧度的随机噪声，从而构成反演数据。

图5 反演结果对比

反演初始模型为100 Ω·m的均匀半空间。正则化因子和DW权重因子的初值分别为λ1=105和λ2=10-4。在迭代过程中，λ1逐渐减小，而λ2逐渐增大。在本例中，λ1的递减系数为1.23，λ2的递增系数为1.60。图5展示了经过第50次迭代后获得的反演结果。从反演结果中可以看到，高阻薄层厚度及其电阻率均得到了很好的重构，其厚度非常接近真实值，其电阻率为151.30 Ω·m，比其真实值偏高。图6展示了第47次至第50次迭代的反演结果。从图6可以看出DW反演方法重构高阻薄层的过程。为了比较起见，图5也给出了传统高斯-牛顿法反演结果。高斯-牛顿法反演方法很好地恢复了第一层和第三层的电阻率,但是对高阻薄层几乎完全没有反映。

由于大地电磁法存在等值性现象，含有高阻薄层的正演响应曲线与不含高阻薄层的正演响应曲线十分接近，使得大地电磁法对高阻薄层的分辨率极低。海洋大地电磁方法对高阻薄层分辨率低的特性制约了其在高阻薄层探测方面的发展，而本例则表明含DW项的反演方法可以通过减小反演拟合残差的自相关性，在一定程度上提高海洋大地电磁方法对高阻薄层的探测能力。

((a)第47次迭代，(b)第48次迭代，(c)第49次迭代，(d)第50次迭代。(a) The 47th iteration， (b) The 48th iteration， (c) The 49th iteration， (d) The 50th iteration.)

图7给出了DW反演迭代过程中均方根RMS、视电阻率对数的DW统计量(DWrho)及相位DW统计量(DWpha)的变化情况。从图7中可以看出，均方根RMS随着迭代次数的增加不断下降，而DWrho和DWpha随迭代次数逐渐增大，并逐渐趋近于它们的期望值。

图7 DW反演迭代过程中RMS与DW值的变化

图8绘出了DW方法第50次迭代反演模型与真实模型的视电阻率对数和相位之间的残差。从图8可以看到，视电阻率对数残差和相位残差在各频段上基本呈现出随机波动的状态。由表 1可知，在数据个数为21的情形下，DW统计量位于(1.54, 2.46)内时，反演数据和反演模型响应之间的残差将不存在自相关性。图5所示DW反演模型与真实模型视电阻率对数残差和相位残差的DW统计量，分别为DWrho= 1.638和DWpha= 1.915，它们均位于1.54～2.46之间，这表明DW反演模型的视电阻率对数残差和相位残差之间都不存在自相关性。

图8 DW反演模型与真实模型的视电阻率残差(a)和相位残差(b)

2.2 八层模型反演

图9中灰色实线所示一个八层地电模型。该模型常被用于验证大地电磁反演方法的效果[21-22]。先在0.000 1～1 000 Hz之间呈对数等间距选取40个频点，经正演计算得到40个频点的大地电磁响应，对视电阻率对数和相位分别添加0.004 3和0.005弧度的随机误差，将其作为反演数据。然后分别使用传统高斯-牛顿反演方法(GN)、OCCAM反演方法和DW反演方法进行反演计算。反演时，初始模型均为100 Ω·m均匀半空间。基于三种反演方法的反演结果如图9所示。从反演结果中可以看出，相较于传统高斯-牛顿反演方法和OCCAM反演方法，DW方法反演结果更接近真实模型，对高阻层和低阻层的识别更加准确。DW反演模型视电阻率对数和相位的DW统计量分别为DWrho=1.208和DWpha= 1.978。由表 1可知，对于40个样本点，当DW统计量位于(1.60, 2.40)区间内时，反演数据和反演模型响应之间的残差将不存在自相关性。而当DW统计量小于1.39时，则认为残差具有正一阶自相关。据此认为，DW反演模型的视电阻率残差之间具有一定的正自相关性，而相位残差不具有自相关性。

图9 DW反演结果与传统GN和OCCAM反演结果对比

图10 DW反演模型与真实模型的视电阻率残差(a)与相位残差(b)

2.3 实测数据反演

在本节中，研究将DW反演方法应用于南黄海实测海洋大地电磁数据反演中。南黄海中部隆起区由于高速海相碳酸盐岩的屏蔽作用，导致地震勘探方法应用不理想，成为制约海相油气勘探的关键[23]。海洋大地电磁方法由于不受高阻层屏蔽、对低阻层敏感度高的特点，已应用于该海域，初步探明海相碳酸盐岩层的厚度[24]。由于研究区水深仅21 m，海水运动(如波浪，潮汐等)使得大地电磁视电阻率和相位在1～10 s周期受到较大干扰，一定程度降低了数据质量和分辨率，因此，在数据处理时需对受海水运动影响较大频点的大地电磁数据进行剔除。下面，给出利用S4站位27个频点的视电阻率数据和19个频点的相位数据的反演结果。

对S4站位实测MT数据进行了DW 反演，并将反演结果与测点附近的测井电阻率数据进行对比(该站位位于CSDP-2井东南方向4 km处[25])。反演时，初始模型为10 Ω·m均匀半空间。图11为该测点实测MT数据的反演结果。由图11可以看出，DW反演结果与传统高斯-牛顿法反演结果[24]相比，对六百多米深处高阻界面的识别更为清晰。图12为DW反演迭代过程中均方根RMS与DW统计量的变化情况。随着迭代次数的增加，均方根误差RMS逐渐下降，而视电阻率与相位的DW值逐渐增大，最后趋于稳定。图13为反演模型的视电阻率和相位与实测数据的对比。图14为它们之间的残差。由表 1可知，在样本数量分别为27与19时，如果DW统计量分别位于区间(1.56, 2.44)与(1.53，2.47)内，则反演数据和反演模型响应之间的残差不存在自相关性。通过对比分析DW值以及拟合曲线可以得到，DWrho= 0.81,相位的DWpha= 0.61, 反演后残差具有正自相关性，这是因为实测数据经过处理后的误差无法满足随机性的特点，由于受海水运动、仪器噪声等影响，实测数据的误差具有较强的正自相关性。但是，较于传统的高斯-牛顿法反演，DW反演结果对浅层高阻层与620 m处的高阻层界面识别更清晰。