近场源浅地表多频电磁法正反演研究

Research on Forward Modeling and Inversion of Multi-Frequency Electromagnetic Method with Near-Field Source in Shallow Surface

摘要: 工程地质及环境问题日渐突出，浅地表电磁探测也日渐发展起来。相较于传统直流电法，近场源频率域电磁探测方法及其装置均有更显著的优势。论文首先介绍了近场源多频电磁法正演公式及数值解法，并对反演方法进行了阐述；然后对文献中的经典模型进行正演，分析了正演结果的精度；最后基于不同的地电模型、不同噪声水平的电磁数据进行已知与未知模型反演，并对比了反演结果。对比反演结果可知，近场源多频电磁法对三层模型中的H型地电模型在较强噪声水平下仍与理论模型吻合较好，且按不固定的未知模型来反演也能较好拟合。研究结果表明，本文正演过程精度高，反方法也表现出了很强的抗干扰能力，为近场源频率域电磁探测方法推广提供理论依据。

Abstract: Engineering geology and environmental problems are becoming more and more prominent, and shallow surface EM detection is also developing. Compared with the traditional direct current method, the method and device of near-field source and frequency domain EM detection have more significant advantages. In this paper, the formula and numerical solution of near-field source and multi-frequency EM method are introduced, and the inversion method is described. Then, the classical models in the literature are performed and the accuracy of the forward results is analyzed. Finally, based on the different electromagnetic models and different noise levels, the known and unknown models are retrieved, and the results are compared. The comparison inversion results show that multi-frequency electromagnetic method of the near-field source is consistent with the theoretical model for the H-type geoelectric model under the strong noise level, and it can also be well fitted by the unfixed unknown model inversion. The research shows that forward modeling has a high precision and the inversion calculation has a high anti-interference ability in this paper, which provides a theoretical basis for the promotion of near-field source and frequency domain EM method.

2. 基本理论

2.1. 磁偶源电磁响应正演公式

本文中采用的GEM-2仪器采用了双线圈，除却接收线圈还有一个参考线圈用来接收一次场 $H_0$ 从而消除接收线圈中的一次场。垂直磁偶源在半空间层状介质中的电磁响应由Frischknecht (1967), Ward (1967), and Ward和Hohmann (1988) [4] [5] [6] 等人给出。

如图1所示，当测量系统位于地面以上h位置时，接收线圈中接收到的二次场 $H_s$ 与一次场 $H_0$ 的比值 [6] 为：

$\frac{H_s}{H_0}=r^3{\displaystyle {\int }_0^{\infty }R\left(\lambda \right){\lambda }^2\exp \left(-2u_{0}h\right)\cdot J_0\left(\lambda r\right)\text{d}\lambda }$ (1)

2.2. Hankel变换

对于解决含有零阶贝塞尔函数的积分，Hankel变换是一种高效的数值方法。其算法 [7] [8] [9] 的描述如下：

若连续函数 $f\left(r\right)$ 表示为 $f\left(r\right)={\displaystyle {\int }_0^{\infty }K\left(\lambda \right)J_0\left(r\lambda \right)\text{d}\lambda }$ ，则 $f\left(r\right)$ 的数值解为：

$f\left(r\right)=\frac{1}{r}{\displaystyle {\sum }_{i=1}^{n}K\left({\lambda }_i\right)C_i}$ (2)

Figure 1. Classical layered geoelectric model

图1. 经典层状地电模型

其中 ${\lambda }_i=\frac{1}{r}\times {10}^{a+\left(i-1\right)s}$ ， 为滤波系数，常用的滤波系数有47点、61点、120点和140点，当 $n=120$ 点时， $a=-8.385$ 为Hankel积分中波数的初值， $s=9.04226468670e^{-2}$ 为波数间隔。

经由此变换，可以得出二次场 $H_s$ 与一次场 $H_0$ 的比值的数值解：

$\frac{H_s}{H_0}=r^3{\displaystyle {\int }_0^{\infty }R\left(\lambda \right){\lambda }^2\exp \left(-2u_{0}h\right)J_0\left(\lambda r\right)\text{d}\lambda }=\frac{1}{r}{r}^3{\displaystyle {\sum }_{i=1}^{n}R\left(\lambda \right){\lambda }^2\exp \left(-2u_{0}h\right)C_i}$ (3)

2.3. 反射系数推导

式中 $R\left(\lambda \right)$ 即前式中的 $r_{TE}$ ，为反射系数 [4] ，可由下式给出：

$r_{TE}=\frac{Y_0-{\hat{Y}}_1}{Y_0+{\hat{Y}}_1}$ (4)

其中， $Y_0=\frac{u_0}{{\hat{z}}_0}$ 为自由空间的本征导纳， ${\hat{Y}}_1$ 为地表导纳， ${\hat{z}}_n=i\omega {\mu }_n$ 。

对于N层层状大地，地表导纳为

${\hat{Y}}_1=Y_1\frac{{\hat{Y}}_2+Y_1\tanh \left(u_1{l}_1\right)}{Y_1+{\hat{Y}}_2\tanh \left(u_1{l}_1\right)}$ (5)

${\hat{Y}}_n=Y_n\frac{{\hat{Y}}_{n+1}+Y_n\tanh \left(u_n{l}_n\right)}{Y_n+{\hat{Y}}_{n+1}\tanh \left(u_n{l}_n\right)}$ (6)

${\hat{Y}}_N=Y_N$ (7)

式中，

$u_n={\left({\lambda }^2-i\omega {\mu }_n{\sigma }_n\right)}^{1/2}$ (8)

由此，从最底层开始向上逐步递推可以得到地表导纳 ${\hat{Y}}_1$ 和第一层地层导纳 $Y_1$ ，从而得出水平层状地层的反射系数 $r_{TE}$ ，进而求取层状介质磁偶源的频率域电磁响应。

2.4. 反演理论

阻尼最小二乘反演问题 [10] 的目的可视为使目标函数趋于极小，对于模型正演公式F(m)，在反演计算中，求解问题为使目标函数φ取值极小：

$\phi ={\left[F\left(m\right)-d\right]}^{\prime }\cdot C\cdot \left[F\left(m\right)-d\right]$ (9)

其中，m为模型向量( $N\times 1$ )，d为数据向量( $M\times 1$ )，C为数据方差倒数的对角方阵( $M\times M$ )。

取 $e=F\left(m\right)-d$ 为当前拟合差向量( $M\times 1$ )，则存在以下关系：

$\text{g}=J^{\prime }\cdot C\cdot e$ (10)

其中， $J\left(M\times N\right)$ 为雅克比矩阵(灵敏度矩阵)， $\text{g}\left(N\times 1\right)$ 为梯度向量

取 $d_m$ 为模型修正量，则高斯牛顿修订量公式为：

$\left(J^{\prime }\cdot C\cdot J\right)\cdot d_m=-J^{\prime }\cdot C\cdot e=-\text{g}$ (11)

阻尼最小二成方程为：

$\left(J^{\prime }\cdot C\cdot J+\lambda I\right)\cdot d_m=-J^{\prime }\cdot C\cdot e=-\text{g}$ (12)

计算机实现反演过程采用迭代格式为：

$m_{k+1}=m_k+inv\left(J^{\prime }\ast C\ast J+\lambda \ast I\right)\ast \left(-\text{g}\right)$ (13)

3. 正演结果及精度分析

对于给定的地电模型，其电磁响应是在每个频率上二次场与一次场比值的同相分量I与正交分量Q。如图2所示，为电导率为 $\sigma =0.01\text{S}/\text{m}$ 的均匀大地电磁响应图。为便于比较，直接采用频率值为横坐标(一般采用响应数，与频率值相对应)。当频率过低时，电磁响应值过小，且与频率相关性小，使得频率测深没有意义。当频率值过高时，同相分量I达到了最大值，而正交分量Q也会急速趋近于零。因此实际工作中所用测量频率也应该在该范围之内。对于GEM-2型多频电磁探测仪而言，其最高工作频率超过48 kHz，低频低于1 kHz，所测的分量值强度高，且与频率相关性高 [4] [11] [12] 。文中正演过程采用的频点值也均在此范围之内。

3.1. 频率域电磁响应计算

图1为典型的层状介质地电模型，第一层电导率取 ${\sigma }_1=0.01\text{S}/\text{m}$ ，接近地表第四系覆盖物的电导率值，厚度取 $l_1=10\text{m}$ 。第二层厚度取 $l_2=5\text{m}$ ，第三层为向下无限延伸。通过改变第二层和底层的电导率值，依次获得五组层状地电模型：均匀层状模型M1，三层Q型地电模型M2，H型地电模型M3，K型地电模型M4以及A型地电模型M5，如表1所示。

通过对这五组模型进行正演计算，得到如图3所示的结果。为便于与已有结果进行比较，正演结果图采用以频点为横坐标，电磁响应值为纵坐标，双对数形式显示。在所取频点范围内，均匀地电模型的同相分量I和正交分量Q均极近似直线变化，I值处于0.7 ppm到60 ppm之间，Q值处于30 ppm到800 ppm之间。相同频率处，正交分量Q值均大于同相分量I值，与已知电磁响应分量关系一致。比较图3中(a)、(b)、(c)与(d)，发现地层电阻率的改变对同相分量及正交分量均有影响。由数值可知，地层电阻率值对正

Figure 2. EM response components I & Q

图2. 电磁响应I & Q分量图

Table 1. Layered geoelectric models

表1. 层状介质模型分类

交分量的影响较小，而对同相分量的影响较大，在三层Q型模型中，I值变化区间为30 ppm~300 ppm，三层H型模型中，I值变化区间为1 ppm~200 ppm，三层K型模型中，I值变化区间为30 ppm~200 ppm，三层A型模型中，I值变化区间为0.09 ppm~30 ppm。分析地电模型，第一层电阻率相同，Q型模型以下为低阻地层，H型模型以下存在低阻薄层和高阻底层，K型模型以下存在高阻薄层和低阻底层，而A型模型以下地层为高阻，且有更高阻的底层。由此可知，地层电阻率越大，电磁响应同相分量I值越小。

3.2. 精度验证

Haoping Huang在之前的研究中 [3] 做过类似研究，我们将所得的结果与其结果进行比较，已验证其可靠性，如图4所示。

所采用模型为三层地电模型，第一层电导率为 ${\sigma }_1=0.01\text{S}/\text{m}$ ，地层厚度为 $l_1=15\text{m}$ ，第二层电导率为 ${\sigma }_2=0.2\text{S}/\text{m}$ ，地层厚度为 $l_2=5\text{m}$ ，底层电导率为 ${\sigma }_3=0.005\text{S}/\text{m}$ 。

由图4中(a)与(b)对比可知，在 $f=2000\text{Hz}$ 时，三层地电模型中对比图中频率响应值的两个分量参考值分别为： $I\approx 7.3\text{ppm}$ ， $Q\approx 88\text{ppm}$ ，正演结果分量值分别为： $I\approx 7.5\text{ppm}$ ， $Q\approx 90\text{ppm}$ ，误差值分别为2.73%和2.27%。误差值均在5%以内，从而验证了正演计算的准确性与可靠性。

4. 反演结果

为了探讨不同地层模型反演时抗干扰能力的强弱，采用对正演结果人为添加噪声，再反演的方式进行。反演过程基于模型参数又可分为两类：已知模型反演与未知模型反演。先进行已知模型反演，反演过程中选取的模型层数与理论一致；再进行未知模型反演，选取模型层数大于理论层数。

4.1. 已知模型反演结果

对于文中采用的四组三层地电模型，均分别添加10%、20%、30%、40%及50%的噪声。反演结果如下列表格所示，为便于比较，选取两次迭代结果列于表格内，首个为最接近实际模型，且迭代次数最小的迭代结果。另外一个为迭代趋于稳定时的迭代结果。本文所用反演过程，一般在20次迭代之内均可趋于稳定，选取第20次迭代结果用于比较。当噪声水平达到50%时，迭代误差极大，迭代不收敛，但也选取第20次迭代结果作为参考。

由表2可知，噪声水平在30%以内时，阻尼最小二乘法对Q型模型的表层和底层反演效果较好，而对中间层效果不佳，即便中间过程可以接近其真实值，但由于误差值水平较高，不能作为最终反演结果，但反演结果能够准确确定地电模型分类。另外其对层厚的反演效果也不佳。

由表3可知，噪声水平在40%以内时，该反演对H型模型电阻率值反演效果较好，除底层电阻率数值外，均接近真实数值，但反演结果能够准确确定地电模型分类。对地层厚度的反演效果较好，同时迭代次数也较小。

Figure 3. EM response forward modeling results of different geoelectric models

图3. 不同地电模型电磁响应正演结果

Figure 4. Compare forward modeling results with the known results

图4. 正演结果与已知文献结果对比

Table 2. Inversion of Q-type geoelectric model

表2. Q型地层模型反演结果

Table 3. Inversion of H-type geoelectric model

表3. H型地层模型反演结果

Table 4. Inversion of K-type geoelectric model

表4. K型地层模型反演结果

由表4可知，噪声水平在10%以内时，该反演对K型模型表层和底层电阻率值反演结果接近真实值，而对中间高阻层反演结果差，完全偏离真实值，且对地层厚度反演效果也较差。其他噪声水平时，反演结果误差更大，且不能确定地电模型。

由表5可知，在各个不同噪声水平下，该反演均能准确确定地电模型。但除了表层电阻率值，其他各层电阻率值均与真实值偏离极大，且对地层厚度反演效果差。

对比上述结果可知，近场源多频电磁法对不同地电模型抗噪效果不尽相同，就三层地电模型而言，效果由好到差有H型>Q型>K型>A型。因此采用近场源多频电磁法完成不同的勘探任务时，大致了解勘探对象的地电特征及噪声干扰情况，为该方法的合理使用及成果解释提供科学指导。

4.2. 未知模型反演结果

我们采用添加10%噪声的三层H型模型的数据来进行未知模型反演。反演过程中采用的初始模型选取了五层地电模型，反演结果如表6所示。

由表6可知，选取初始模型层数大于理论层数进行反演，得出的结果依旧可以接近理论数值。上面三层因为等值现象的关系可以合并为H型模型的表层，而第四层的低阻薄层目标层十分接近理论数值，

Table 5. Inversion of A-type geoelectric model

表5. A型地层模型反演结果

Table 6. Inversion of H-type geoelectric model

表6. H型模型未知模型反演结果

而底层的高阻层也与理论模型一致。可以得出结论：即使选取的初始模型不固定，即在反演初始模型不清楚的情况下，反演结果的拟合效果依然相对较好。

5. 结论

通过与现有文献的结果对比，证明了本文的正演结果正确，反演方法也表现出了很强的抗干扰能力，特别是对于H型地电模型而言，反演结果好，精确度高，为近场源频率域电磁测法提供了可靠的理论依据。即便是初始模型不确定时，该反演依然能够有较好的拟合效果。近场源多频电磁法野外施工方便，资料处理可靠，该方法在浅地表电磁勘探中具有广阔的应用前景。

基金项目

国家科技基础性工作专项项目(2013FY110800)。

NOTES

^*通讯作者。