基于最小泛函的CSAMT视电阻率计算方法

杨轮凯， 胡 建， 叶景艳， 肖明辉， 杨华忠

(华东有色地球化学勘查与海洋地质调查研究院，南京 210007)

基于最小泛函的CSAMT视电阻率计算方法

杨轮凯， 胡 建， 叶景艳， 肖明辉， 杨华忠

(华东有色地球化学勘查与海洋地质调查研究院，南京 210007)

可控源电磁时间域数据的视电阻率计算,因其原理简单，计算处理速度快，可以较快反映地下介质电阻率分布情况，而在实际生产中得到广泛应用。传统的视电阻率算法采用求解均匀半空间电磁场响应函数的根，没有从整体上考虑野外数据中存在的噪声对视电阻率地影响。利用最优化理论，考虑野外数据整体的拟合情况，采用求解特定目标泛函的变分问题，提出了一种可控源电磁勘探时间域数据视电阻率求解方法，在实际应用中得到了优于传统方法的计算结果。

CSAMT; 视电阻率; 时间域; 最小泛函视电阻率算法; 最优化理论

0 引言

可控源音频大地电磁测深法(CSAMT)是一种利用人工发射电流源激发交变电磁场，观测电磁场变化，进而研究地下介质电阻率分布特征的勘探方法，属于频率域电磁测深方法的一种。 Strangway等[1]针对大地电磁法场源的随机性和信号微弱，提出了一种改进方案(采用可以控制的人工场源),CSAMT在地下水资源、地热、工程物探、油气检测、金属矿、煤田地质灾害预测等领域得到了广泛地应用，并取得了较好的效果[2]。CSAMT是在大地电磁法(MT),音频大地电磁法(AMT)的基础上发展而来的,MT、AMT采用测量相互正交的电场和磁场来计算视电阻率—卡尼亚视电阻率。由于MT、AMT利用高空电离层和赤道上空的雷电激发的电磁波作为发射信号，当电磁波到达地面时，电磁波已经可以近似为平面波，按照平面波入射理论进行资料处理和解释。

图1 可控源电磁勘探标量装置示意图Fig.1 Controlled source electromagnetic prospecting scalar device schematic diagram

图1 可控源电磁勘探标量装置示意图Fig.1 Controlled source electromagnetic prospecting scalar device schematic diagram

目前，CSAMT的资料处理一般采用三种方法：①不用近场及过渡区的资料，只取远场数据进行反演；②Routh等[3-7]提出的全区资料一维反演；③进行近场校正，将CSAMT资料校正成MT资料后，再用MT方法反演。第一种方法，需要加大发射源与接收点的距离，使得CSAMT遇到和MT同样的问题，即信噪比降低，容易受到工区噪声的干扰，同时也限制了CSAMT在实际生产中的应用和造成较大的资料浪费；第二种方法，对初始模型比较敏感，同时反演过程中需要计算场源产生的电磁场，计算量大，处理速度慢；第三钟方法，罗延钟等[8-11]提出的采用近场校正，抑制由于场源影响产生的视电阻率畸变。上述校正方法都是基于均匀半空间中电偶源的响应公式，当地下不均匀时，近场校正的应用效果不理想。汤井田等[10]提出了对视电阻率的重新定义来克服上述困难，其基本原理是建立观测信号与同一观测装置参数的均匀半空间模型响应的关系方程，再求方程的根，即为所需的视电阻率。这种视电阻率计算方法存在两方面的问题：①实际观测的数据中不可避免的存在噪声干扰，即便采用平滑滤波等抑制噪声的影响依然难以避免视电阻率曲线包含假信息；②由于电磁勘探物理问题的多解性，需要人机交互挑选，不可避免地引入人为干扰。

笔者提出了一种可控源电磁勘探时间域数据视电阻率求解新方法，最小泛函拟合视电阻率计算方法，该方法构造满足特定条件目标泛函，通过迭代搜索目标泛函的极小，得到视电阻率光滑曲线函数。该方法综合考虑了数据噪声对视电阻率造成的影响，同时也克服了电磁勘探中的多解性问题。

1 最小泛函拟合视电阻率计算方法

传统的视电阻率计算是依据均匀半空间中电磁场响应与观察数据等式关系，求解正演公式的根，见式(1)。

d=f(m)±e

(1)

其中：d为实际观测数据；f(m)为均匀半空间的正演公式；e为数据误差。对式(1)求根便得到该数据对应的视电阻率结果，见式(2)。

ma=f-1(d±e)=f-1(d)±ε

(2)

其中：d为所求的视电阻率；f(m)为计算视电阻率的误差。通过式(1)和式(2)可以发现，观测数据中的噪声对视电阻率计算的影响，在强噪声环境中，对视电阻率结果造成不可预测的结果。在收发距较小，大地电阻率较高或工作频率较低时，观测的电磁场属于近场区或过渡区场，在此情况下，视电阻率将发生畸变，不能直观反映地下介质的电性分布，不利于资料地处理解释。传统的视电阻率计算方法没有整体考虑数据拟合，造成视电阻率曲线存在“毛刺”，使得后续反演出现假异常的情况。

根据地球物理电磁勘探的原理，视电阻率随时间变化的函数是一个光滑的函数，即视电阻率随时间的函数的二阶导数积分最小，如式(3)所示。

(3)

其中：t0、tn为野外观测的起止时刻。

同时，要求正演数据和野外观测数据的拟合度小于等于给定的门限，如式(4)所示。

(4)

其中：f(ma(ti))为均匀半空间中电磁时间域正演计算公式；S2为给定的拟合度门限值。

引入拉格朗日乘子P，构造一个目标泛函U，如式(5)所示。

(5)

引入二阶导数算法，式(5)可以改写为式(6)。

(6)

∂t∂ma+p(-j)t(d0-jma)=0

(7)

具体算法的流程如下：

1)设置当前迭代次数为k，如果k=1，初始化视电阻率；否则，将上次迭代的视电阻率结果作为初始值。计算初始数据拟合度，记为rms0; 计算初始视电阻率粗糙度，记为ruf0；计算Jacobian矩阵，Jacobian矩阵(具体步骤引用文献[12])。设置拉格朗日乘子的初始值，一般令p=105。令搜索步长β=1。

3)令m等于当前新的视电阻率计算m的数据拟合，记为rms1；计算当前新的视电阻率m的粗糙度，记为ruf1。

4)如果rms1小于等于给定的拟合度门限，转步骤5)；否则，如果rms1小于rms0，转步骤6)，否则，令k=k+1，转步骤2)。

5)搜索一个最优的拉格朗日乘子P，使得rms1等于S2，且P尽可能的大。

6)如果ruf1大于ruf0，令k=k+1，转步骤2)。

7)本次迭代结束，m即为最优的视电阻率。令k=k+1，转步骤1)。

图2 基于最小泛涵的视电阻率计算流程图Fig.2 Based on the minimum the apparent resistivity calculation flow chart of the culvert

2 均匀半空间时间域电磁场正演

可控源电磁方法一般常采用接地电偶极源，笔者着重讨论水平电偶极源的均与半空间响应的计算。纳比吉安[13]采用TE分量、TM分量分解的方法，推导了电偶极源的电磁场表达式，形如式(8)、式(9)。

(8)

(9)

由于实际施工中采用的接地电源是长电流源，可以看做多个小电偶极源的线积分。为了兼顾计算精度和速度，可以采用高斯积分公式，模拟野外长电流源的电磁响应。

3 野外数据测试

我们选取了某工区一条测线的可控源电磁实际施工采集数据进行测试。反演方法采用Constable[12]提出的Occam反演方法，Occam方法是一种的正则化反演方法，对初始模型要求不高。图3为传统视电阻率算法结果的Occam反演剖面，从图3中可以看出，传统方法求得的视电阻率的Occam反演剖面存在很多“挂面条”现象，这是因为实际采集数据中存在噪声，而传统的视电阻率求解方法没有考虑到噪声的影响，从而对视电阻率数据造成干扰。图4是最小泛函视电阻率算法结果的反演剖面与工区内地震、测井推断的地质信息叠合图。由图4可以发现，反演剖面的横向连续性得到改善，没有明显的“挂面条”现象，说明笔者提出的视电阻率求解方法较好抑制原始数据中的噪声对视电阻率的影响，且19km处的断裂与已知吻合，而在传统的视电阻率算法结果的反演剖面中，该断裂没有得到较好地反映，主要是因为受到邻近测点的从浅至深的高阻条带的影响。

4 结论

通过测试，我们发现最小泛函视电阻率算法，通过整体拟合，考虑各个测点数据的噪声情况，有选择地保留有效信息、压制干扰，从而最大限度获取地下介质的电阻率分布情况。并取得了以下几点认识：

1)笔者提出的方法对数据中的噪声有一定地抑制作用，在求解方法中将数据中的噪声作为数据的权重(噪声越大，权重越小)，一定程度上克服噪声对视电阻率影响。

图3 某工区一条测线的传统视电阻率数据的反演结果Fig.3 The traditional apparent resistivity data inversion results of a line in a work area

图4 某工区一条测线的最小泛函视电阻率数据的视模反演结果与地震、测井数据推断的地质信息对比情况Fig.4 The minimum functional visual apparent resistivity data of a line in a work area, overlay with the geological information interpolated by seismic, well logging data

2)因为该方法采用迭代搜索，无论从什么初始条件出发，总能搜索到是观测数据拟合误差达到最小的视电阻率，即该方法具有一定的数值稳定性。

3)方法暂未考虑地形起伏的影响，在实际工作中，如果工区地形起伏较大时，需要用二维的方法进行处理。

4)可以利用现有的MT反演程序，避免全区资料反演方法计算量大的缺点，可以实现野外现场反演。

[1] GOLDSTEIN M.A.，STRANGWAY D .W. Audio-frequency magnetotelluric with a grounded electric dipole source[J]. Geophysics, 1975, 40 (1):669-683.

[2] 何继善.可控源音频大地电磁法[M].长沙:中南工业大学出版社，1990. HE J S.Control source audio frequency magnetotelluric[M].Changsha:Central South University of Technology Press,1990.(In Chinese)

[3] ROUTH P S, OLDENBURG D W. Inversion of controlled-source audio-frequency magnetotellurics data for a horizontally layered earth[J]. Geophysics, 1999, 64(6): 1689-1697.

[4] 王若，王妙月. 一维全资料CSAMT反演[J]. 石油地球物理勘探，2007,42(1): 107-114. WANG R, WANG M Y. Inversion of 1-D full CSAMT data[J] . OGP, 2007, 42(1):107-114.(In Chinese)

[5] 朱威，李桐林，尚通晓，等. CSAMT一维全区反演对比研究[J]. 吉林大学学报，2008,38(增刊):12-14 ZHU W, LI T L, SHANG T X, et al. Campared study of 1-D full region inversion of CSAMT[J]. Journal of Jilin University, 2008,38(s):12-14.(In Chinese)

[6] 何兴梅，胡祥云，陈玉萍，等. 奥克姆一维反演的应用[J]. 工程地球物理学报，2008,5(4): 439-443. HE M M,HU X Y,CHENG Y P, et al. Application of 1D Occam' s inversion to CSAMT[J]. Chinese Journal of Engineering Geophysics, 2008,5(4): 439-443.(In Chinese)

[7] 汤井田，张林成，肖晓，等. 基于频点CSAMT一维最小构造反演[J]. 物探化探计算技术，2011,33(6): 577-581. TANG J T，ZHANG L C，XIAO X, et al . One dimension CSAMT minimum structure inversion based on the frequency[J]. Computing Techniques for Geophysical and Geochemical Exploration，2011,33(6): 577-581.(In Chinese)

[8] YAMASHITA,HOLLOF.CSAMT case histories with a multi-channel CSAMT system and discussion of near-field data correction[J].Phoenix Geophys Ltd,1985(1)：276-278.

[9] 罗延钟，周玉水，万乐. 一种新的CSAMT资料近场校正方法，勘查地球物理勘查地球化学文集第20集[M]. 北京:地质出版社，1996. LUO Y Z, ZHOU Y S, WAN L. A new method for CSAMT near-field correction, 20thcorpus of geophysical exploration&geochemical Exploration[M]. Beijing:Geological Publishing House,1996.(In Chinese)

[10]汤井田，何继善. 水平电偶源频率域测深中全区视电阻率定义的新方法[J]. 地球物理学报，1994，37(04): 543-552. TANG J T, HE J S . A new method to define the full-zone resistivity in horizontal electric dipole frequency soundings on a layerd earth[J]. Chinese Journal Geophysics,1994,37(04): 543-552.(In Chinese)

[11]汤井田，何继善，水平多层介质上水平电偶源频率域电磁测深的阻抗实部等效电阻率[J].物探与化探，1994(3): 92-96. TANG J T,HE J S. Impedance reaL part equivalent resistivity in frequency electromagetic sunding of horizontal electric double source on horizontal multilayer media[J]. Geophysical &Geochemical Exploration，1994(3): 92-96.(In Chinese)

[12]CONSTABLE, S., R PARKER, C. G. CONSTABLE. Occam`s inversion : A practical algorithm for generating smooth models from electromagnetic sounding data[J].Geophysics,1989,52(3): 289-300.

[13]米萨克 N,纳比吉安.勘查地球物理电磁法 第一卷 理论[M].北京：地质出版社，1992. MISAC N,NANIGHIAN Electromagnetic methods in applied geophysics, vol1, theory[M].Beijing:Geological publishing house,1992.

Based on the minimum functional calculation method of the CSAMT apparent resistivity

YANG Lunkai, HU Jian, YE Jingyan, XIAO Minghui, YANG Huazhong

(Institute of Geochemical Exploration and Marine Geological Survey,ECE ,Nanjing 210007， China)

The calculation of the apparent resistivity in time domain of controlled source electromagnetic method is widely used in geophysical exploration for its low time costing, which could quickly acquire the distribution of the ground resistivity. The conventional apparent resistivity calculating algorithm uses the root of the response function of the uniform half-space electromagnetic field without considering the effect of the noise in the raw data on the apparent resistivity as a whole. In this paper, we proposed a method to calculate the apparent resistivity in time domain for controlled-source electromagnetic method by using the optimization theory to solve the variational problem of specific target functional. The fitting accuracy of field data is considered as a whole in one section. The method has been used in field operation, and acquired more acceptable result than those of the conventional calculation.

CSAMT ; apparent resistivity; time domain; minimal functional apparent resistivity method; optimization theory

2016-12-25 改回日期：2017-02-27

江苏省地质勘查专项资金项目(苏财建[2016]140)

杨轮凯(1971-)，男，高级工程师，主要从事电磁法及应用研究，E-mail：584945583@qq.com。

胡建(1980-)，男，高级工程师，主要从事矿产地质调查工作，E-mail：65153743@qq.com。

1001-1749(2017)03-0301-05

P 631.4

A

10.3969/j.issn.1001-1749.2017.03.01