基于监督下降法的直流电测深曲线反演

Inversion of Direct Current Sounding Curve Based on Supervised Descent Method

摘要: 电阻率反演是直流电测深资料最重要的定量解释方法之一，但常规基于梯度信息的反演算法中灵活实现层厚和电阻率的区间约束比较困难。为此，本文将监督下降法引入电测深曲线反演以实现灵活的先验信息引入，进而重构地下模型参数。利用正演合成的数据，论证了训练集中不同初始模型的选择对SDM反演可行性的影响。此外，开展了两组数值实验进一步探讨SDM的泛化能力。结果表明，采用SDM融入先验信息具有可行性。该方法反演过程中不涉及偏导数的计算，不仅可以克服对初始模型的依赖达到快速收敛的目的，并且具有一定的泛化能力。

Abstract: Resistivity inversion is one of the most important quantitative interpretation methods for DC sounding data, but it is difficult to flexibly implement interval constraints of layer thickness and resistivity in conventional gradient information-based inversion algorithms. To this end, this pa-per introduces the supervised descent method into the inversion of electrical sounding curves to achieve flexible introduction of prior information, and then reconstruct the parameters of the underground model. Using the data synthesized by forward modeling, the influence of the selection of different initial models in the training set on the feasibility of SDM inversion is demonstrated. In addition, two sets of numerical experiments were carried out to further explore the generalization ability of SDM. The results show that it is feasible to use SDM to incorporate prior information. The inversion process of this method does not involve the calculation of partial derivatives, which can not only overcome the dependence on the initial model to achieve fast convergence, but also has certain generalization ability.

1. 引言

直流电阻率法被广泛应用于水文、工程、环境地质调查等多个领域 [1]。电阻率反演是直流电测深资料最重要的定量解释方法之一，但其本身是一个非线性问题，许多学者为此做了大量研究。早期电阻率探测中，通常采用量板模型人工解释的方法对探测区域进行地质推断 [2] [3]，随着Backus-Gilbert反演理论的提出，Inman [4] 和Parker等 [5] 将其推广至离散模型，之后反演理论及方法得到发展。目前已有的直流电阻率反演方法如：平滑约束的最小二乘法、共轭梯度法、拟牛顿法、高斯牛顿法等 [6] [7] [8] [9] 都发挥了重要的作用。以上方法可以归属于基于梯度下降的传统线性化迭代反演，其通过最小化目标函数迭代计算反演模型参数。但是此类方法涉及偏导数的计算，引入先验信息的方式不够灵活，反演中对初始模型的依赖性较强，导致迭代更新时易陷入局部极小。机器学习技术可以从有明确的输入与输出的任务中学习其中的映射关系，然后使用该映射关系实现给定的输入和输出。近年来，其作为一种“通用技术”成功应用于多个领域，如：自然语言处理、语音识别、无人驾驶、医学图像、游戏战略规划等 [10] [11] [12] [13] [14]。与此同时基于机器学习技术的算法在地球物理领域中也是研究的前沿，如：McCormack等 [15] 基于反向传播神经网络实现了地震数据的自动道编辑和初至提取，处理结果与人工方法具有高一致性；徐海浪等 [16] 采用BP神经网络实现了直流电阻率二维反演，得到了比传统方法更优的反演效果；罗飞等 [17] 提出带约束的Markov决策过程，为在信噪比较低情况下准确进行初至走时的自动拾取提供了方法。Shreedhar等 [18] 将OCSVM应用于探地雷达数据处理中，实现了对路面薄层的检测。

监督下降法(supervised descent method, SDM)属于机器学习方法中的一种，最早Xiong和De la Torre等 [19] 于2013年提出用于解决人脸识别问题。它首先在训练集中学习并记录目标函数的下降方向，然后利用此方向直接更新未知模型来解决优化问题。SDM被成功应用于胸部EIT成像、微波图像重建、睡意实时监测等 [20] [21] [22] 方面。在地球物理领域，Guo等 [23] [24] SDM应用于处理大地电磁数据，证实了该方法在大地电磁领域中的可行性。Hu等 [25] 在定向电磁随钻测井问题中成功应用SDM，实现了对地下信息的随钻实时更新反演。Peng Hao等 [26] 应用SDM实现了各向异性地层测井问题的反演。已有研究表明，SDM在地球物理反演领域具有较大应用潜力，为了论证SDM在直流电测深数据反演问题中的可行性，本文针对层状模型应用SDM。

2. 水平层状地电模型正演问题

电测深水平地层正演模型如图1。假定水平地面下有n层水平层状地面，地层电阻率从上至下分别为 ${\rho }_1,{\rho }_2,\cdots ,{\rho }_n$，厚度分别为 $h_1,h_2,\cdots ,h_n$，每层底面到地面的距离分别为 $H_1,H_2,\cdots ,H_{n-1},H_n=\infty$。在地面有一个点电流源，电流强度为I。

电阻率转换函数的递推公式为：

$\{\begin{array}{l}{T}_i\left(\lambda \right)={\rho }_i\frac{1+\frac{T_{i+1}\left(\lambda \right)-{\rho }_i}{T_{i+1}\left(\lambda \right)+{\rho }_i}{\text{e}}^{-2\lambda h_i}}{1-\frac{T_{i+1}\left(\lambda \right)-{\rho }_i}{T_{i+1}\left(\lambda \right)+{\rho }_i}{\text{e}}^{-2\lambda h_i}}\\ T_n\left(\lambda \right)={\rho }_n\end{array}$ (1)

电测深视电阻率表达式为：

${\rho }_s\left(r\right)=r^2{\displaystyle {\int }_0^{\infty }{T}_1}\left(\lambda \right)J_1\left(\lambda r\right)\lambda \text{d}\lambda$ (2)

其中，r为供电极距AB/2， $T_1\left(\lambda \right)$ 为电阻率转换函数， $\lambda$ 为积分系数， $J_1$ 为第一类贝塞尔函数。正演时采用数字滤波法求解贝塞尔函数在 $\left(0,\infty \right)$ 上的数值积分。

Figure 1. Horizontal layered medium model

图1. 水平层状介质模型

3. SDM理论及反演算法

3.1. SDM理论

水平层状地层中直流电测深数据的反演问题可以描述为式(3)并使其最小化：

$S\left(m\right)={\Vert F\left(m\right)-d_{obs}\Vert }_2^2$ (3)

式中F为正演模型算子， $d_{obs}$ 为观测数据矢量， $m$ 为模型参数矢量， $S\left(m\right)$ 为目标函数。

对目标函数 $S\left(m\right)$ 在 $m=m_0+\Delta m$ 处进行二阶泰勒展开并忽略二阶以上高阶项，即

$S\left(m_0+\text{Δ}m\right)\approx S\left(m_0\right)+J_S{\left(m_0\right)}^{\text{T}}\text{Δ}m+\frac{1}{2}\text{Δ}{m}^{\text{T}}{H}_S\left(m_0\right)\text{Δ}m$ (4)

式中 $m_0$ 是模型参数的初始值， $J_S$ 与 $H_S$ 是 $S\left(m\right)$ 在 $m_0$ 处的雅可比与Hessian矩阵。对上式求极小可得到式(5)：

$\text{Δ}m=-H_S{\left(m\right)}^{-1}{J}_S\left(m\right)=-2H_S{\left(m\right)}^{-1}{J}_F{\left(m\right)}^{\text{T}}\left(F\left(m\right)-d_{obs}\right)=R\text{Δ}d$ (5)

其中 $J_F\left(m\right)$ 是正演响应函数的雅可比矩阵， $\text{Δ}d=F\left(m\right)-d_{obs}$ 是正演数据与观测数据之差， $R$ 是由 $\text{Δ}d$ 和 $\text{Δ}m$ 学习的下降方向。基于梯度下降的方法(如牛顿法)采用的是迭代的方法逼近局部最小值。相较于此类方法，SDM根据先验信息生成模型进行正演得到数据，由 $\text{Δ}m$ 与 $\text{Δ}d$ 迭代直接求得一系列下降方向，然后将其直接应用于反演中更新模型参数。训练过程中N个模型的下降方向可通过迭代求解式(6)得到：

$\underset{R_k}{\arg \min }\left(\underset{n=1}{\overset{N}{{\displaystyle \sum }}}{\Vert \text{Δ}{m}_k^n-R_k\text{Δ}{d}_k^n\Vert }_2^2\right)$ (6)

其中

$\begin{array}{l}\text{Δ}{m}_k^n=m_{train}^n-m_k^n\\ \text{Δ}{d}_k^n=d_{train}^n-d_k^n\end{array}$ (7)

式中N是训练样本的数量， $m_{train}^n$ 是训练集中第n个模型的参数， $m_k^n$ 是第k次迭代时第n个训练模型更新后的模型参数， $d_{train}^n$ 和 $d_k^n$ 是 $m_{train}^n$ 和 $m_k^n$ 的相应正演模拟数据。

通过求解式(7)，迭代时模型参数的更新如下式：

$m_k^n=m_{k-1}^n+R_{k-1}\text{Δ}{d}_{k-1}^n$ (8)

3.2. SDM反演

SDM反演可分为两个阶段：离线训练阶段与在线反演阶段。

3.2.1. 离线训练阶段

对式(6)以矩阵形式写成

$\underset{R_k}{\arg \min }\left({\Vert \text{Δ}{M}_k-\text{Δ}{D}_k{R}_k^{\text{T}}\Vert }_F^2\right)$ (9)

其中

$\text{Δ}{M}_k=\left[\begin{array}{c}\text{Δ}{m}_k^{1\text{T}}\\ \text{Δ}{m}_k^{2\text{T}}\\ ⋮\\ \text{Δ}{m}_k^{N\text{T}}\end{array}\right],\text{Δ}{D}_k=\left[\begin{array}{c}\text{Δ}{d}_k^{1\text{T}}\\ \text{Δ}{d}_k^{2\text{T}}\\ ⋮\\ \text{Δ}{d}_k^{N\text{T}}\end{array}\right]$ (10)

离线训练的第k步中，对式(9)求偏导得：

$R_k^{\text{T}}={\left(\text{Δ}{D}_k^{\text{T}}\text{Δ}{D}_k+\mu I\right)}^{-1}\text{Δ}{D}_k\text{Δ}{M}_k$ (11)

其中 $\mu$ 是用来稳定解的正则化因子，在此将其取为与 $\text{Δ}{D}_k$ 的最大特征值成比例 [24]。初始模型 $m_0$ 可选取为层状均匀模型。

迭代时模型更新为：

$M_{k+1}=M_k+\text{Δ}{D}_k\cdot R_k^{\text{T}}$ (12)

定义模型误差为：

${\text{rms}}_M=\frac{1}{N}\underset{n=1}{\overset{N}{{\displaystyle \sum }}}\frac{{\Vert \text{Δ}{m}_k^n\Vert }_2}{{\Vert m_T^n\Vert }_2}$ (13)

定义数据误差为：

${\text{rms}}_D=\frac{1}{N}\underset{n=1}{\overset{N}{{\displaystyle \sum }}}\frac{{\Vert \text{Δ}{d}_k^n\Vert }_2}{{\Vert d_T^n\Vert }_2}$ (14)

3.2.2. 在线反演阶段

在线反演阶段见式(15)：

$m_{k+1}=m_k+R_k\left(d_{obs}-F\left(m_k\right)\right)$ (15)

在线反演阶段初始模型 $m_0$ 与离线训练阶段需保持一致。

定义数据误差为：

${\text{rms}}_D=\frac{{\Vert d_{obs}-F\left(m_k\right)\Vert }_2}{{\Vert d_{obs}\Vert }_2}$ (16)

3.3. SDM流程

离线训练阶段：

1) 生成训练集，给定初始模型 $m_0$，设置循环终止条件等参数；

2) 第k次迭代时，计算 $\text{Δ}{M}_k$ 与 $\text{Δ}{D}_k$，利用正则化因子求得稳定解 $R_k$ ；

3) 更新模型： $m_k^n=m_{k-1}^n+R_{k-1}\text{Δ}{d}_{k-1}^n$ ；

4) 重复步骤2)~3) 再次循环，直至满足循环终止条件，离线训练阶段结束。

在线反演阶段：

1) 设置与离线训练阶段一致的初始模型 $m_0$，设置循环终止条件等参数；

2) 第k次迭代时，计算 $d_{obs}-F\left(m_k\right)$ 与 $m_{k+1}$ ；

3) 重复步骤2) 再次循环，直至满足循环终止条件，在线反演阶段结束。

4. 模型算例

构建三层地层模型，本节设置了三组实验以测试反演方法的可行性，收敛性与泛化能力。所有数据均为正演模拟数据，反演的模型参数可表示为 $m=\left[\rho ,h\right]$。反演的参数为每一层的边界位置和电阻率。

4.1. 可行性分析

在可行性分析中，研究不同初始模型对SDM反演的影响并与最小二乘反演结果作对比。 ${\rho }_1$， ${\rho }_2$ 和 ${\rho }_3$ 分别为第一层至第三层的电阻率。首先假定大致圈定地层参数(层厚与视电阻率)范围即先验信息，然后在地层参数范围内随机改变电阻率值与层厚值，共生成100个模型得到一个训练集，在离线训练阶段得到下降矩阵 $R$。地层参数范围及初始模型的选取见表1，其中采用真实模型 $m_1$ 对训练过程中初始模型的选取进行测试，表中 $h_1$ 和 $h_2$ 分别为第一层和第二层的厚度。

在表1中初始模型1参数远离真实模型 $m_1$ 参数，相比之下初始模型2参数接近于真实模型 $m_1$ 参数值。图2为不同初始模型离线阶段模型误差与数据误差。从图2中不难看出，在不同的初始模型的选取下，SDM离线阶段均具有很快的收敛速度，十步之内模型误差与数据误差均可达到一个很小的值。图3为不同初始模型下SDM在线反演阶段数据误差，收敛速度同样很快，其反演结果在十步内收敛。此外，从表2中不难看出，当初始模型较准确与较差时，SDM均可对真实模型较为准确的估计。以上表明，SDM反演几乎不受初始模型的影响，其能够较好的解决该非线性反演问题。

Table 1. The K-type geoelectric model: parameters during the offline stage and the real model parameters during the online inversion stage

表1. K型地电模型：离线训练阶段模型参数及在线反演阶段真实模型参数

Figure 2. The training based on different initial models during the SDM offline stage. (a) Model misfit; (b) Data misfit

图2. SDM离线阶段基于不同初始模型训练。(a) 模型误差；(b) 数据误差

Figure 3. The data misfit during the SDM online stage based on different initial models

图3. SDM在线阶段基于不同初始模型的数据误差

Table 2. Inversion results of model m1 based on different initial models

表2. 基于不同初始模型的模型m₁反演结果

4.2. 泛化能力分析

文中“泛化能力”指的是对新鲜样本的反演能力。为了研究SDM的泛化能力，设置了两组模拟实验。两组实验中的训练集参数仍采用表1，训练时采用初始模型1，反演最大迭代次数仍为十次。第一组实验中对表3中给定的前四个真实模型进行反演，以进一步分析反演算例与训练集相同层数时，SDM对训练集内参数与训练集外参数泛化能力。模型 $m_2$ 的电阻率值和层厚值都在训练集范围内，模型 $m_3$ 的电阻率值在训练集范围内但层厚值不在，模型 $m_4$ 的层厚值在训练集范围内但电阻率值不在，模型 $m_5$ 的电阻率和层厚都不在训练集范围内。四个真实模型最终反演结果见图4，反演结果与实际有偏差但均在可接受范围内。这是因为SDM的训练集中包含了模型电阻率与层厚值分布的弱先验信息，SDM优化了训练数据的下降方向，以实现全局优化。

Table 3. Generalization ability test

表3. 泛化能力测试

Figure 4. The inversion model has the same number of layers as the training set: SDM inversion results for four models (a) m₂ (b) m₃ (c) m₄ (d) m₅

图4. 反演模型与训练集层数相同：四种模型SDM反演结果 (a) m₂ (b) m₃ (c) m₄ (d) m₅

第二组实验中对给定层数为两层与四层的真实模型进行反演，以进一步分析由三层模型的训练集训练得到的下降方向是否对多层地层具有泛化能力，表3中模型 $m_6$ 与 $m_7$ 反演结果见图5，反演参数在真实模型的合理范围内。结果表明，由三层模型训练的下降方向能够向多层地层推广。

Figure 5. Two-layer model and four-layer model SDM inversion results (a) m₆ (b) m₇

图5. 二层模型与四层模型SDM反演结果 (a) m₆ (b) m₇

5. 结论

本文采用监督下降法对电测深曲线进行反演，通过设置几组数值实验得出以下几点结论：

1) 将SDM应用于直流电测深数据反演中具有可行性。一维层状K型地电模型测试结果显示，SDM可以将约定层厚与电阻率范围的先验信息灵活地在训练阶段引入并应用于反演中。其能够避免雅可比矩阵的求解且不依赖于初始模型的选择。

2) SDM具有泛化能力，能够在弱先验信息的情况下得到相对合理的反演结果，但其泛化能力边界仍有待研究。