楊輪凱， 胡 建， 葉景艷， 肖明輝， 楊華忠

(華東有色地球化學(xué)勘查與海洋地質(zhì)調(diào)查研究院，南京 210007)

基于最小泛函的CSAMT視電阻率計算方法

楊輪凱， 胡 建， 葉景艷， 肖明輝， 楊華忠

(華東有色地球化學(xué)勘查與海洋地質(zhì)調(diào)查研究院，南京 210007)

可控源電磁時間域數(shù)據(jù)的視電阻率計算,因其原理簡單，計算處理速度快，可以較快反映地下介質(zhì)電阻率分布情況，而在實際生產(chǎn)中得到廣泛應(yīng)用。傳統(tǒng)的視電阻率算法采用求解均勻半空間電磁場響應(yīng)函數(shù)的根，沒有從整體上考慮野外數(shù)據(jù)中存在的噪聲對視電阻率地影響。利用最優(yōu)化理論，考慮野外數(shù)據(jù)整體的擬合情況，采用求解特定目標泛函的變分問題，提出了一種可控源電磁勘探時間域數(shù)據(jù)視電阻率求解方法，在實際應(yīng)用中得到了優(yōu)于傳統(tǒng)方法的計算結(jié)果。

CSAMT; 視電阻率; 時間域; 最小泛函視電阻率算法; 最優(yōu)化理論

0 引言

可控源音頻大地電磁測深法(CSAMT)是一種利用人工發(fā)射電流源激發(fā)交變電磁場，觀測電磁場變化，進而研究地下介質(zhì)電阻率分布特征的勘探方法，屬于頻率域電磁測深方法的一種。 Strangway等[1]針對大地電磁法場源的隨機性和信號微弱，提出了一種改進方案(采用可以控制的人工場源),CSAMT在地下水資源、地?zé)?、工程物探、油氣檢測、金屬礦、煤田地質(zhì)災(zāi)害預(yù)測等領(lǐng)域得到了廣泛地應(yīng)用，并取得了較好的效果[2]。CSAMT是在大地電磁法(MT),音頻大地電磁法(AMT)的基礎(chǔ)上發(fā)展而來的,MT、AMT采用測量相互正交的電場和磁場來計算視電阻率—卡尼亞視電阻率。由于MT、AMT利用高空電離層和赤道上空的雷電激發(fā)的電磁波作為發(fā)射信號，當(dāng)電磁波到達地面時，電磁波已經(jīng)可以近似為平面波，按照平面波入射理論進行資料處理和解釋。

圖1 可控源電磁勘探標量裝置示意圖Fig.1 Controlled source electromagnetic prospecting scalar device schematic diagram

圖1 可控源電磁勘探標量裝置示意圖Fig.1 Controlled source electromagnetic prospecting scalar device schematic diagram

目前，CSAMT的資料處理一般采用三種方法：①不用近場及過渡區(qū)的資料，只取遠場數(shù)據(jù)進行反演；②Routh等[3-7]提出的全區(qū)資料一維反演；③進行近場校正，將CSAMT資料校正成MT資料后，再用MT方法反演。第一種方法，需要加大發(fā)射源與接收點的距離，使得CSAMT遇到和MT同樣的問題，即信噪比降低，容易受到工區(qū)噪聲的干擾，同時也限制了CSAMT在實際生產(chǎn)中的應(yīng)用和造成較大的資料浪費；第二種方法，對初始模型比較敏感，同時反演過程中需要計算場源產(chǎn)生的電磁場，計算量大，處理速度慢；第三鐘方法，羅延鐘等[8-11]提出的采用近場校正，抑制由于場源影響產(chǎn)生的視電阻率畸變。上述校正方法都是基于均勻半空間中電偶源的響應(yīng)公式，當(dāng)?shù)叵虏痪鶆驎r，近場校正的應(yīng)用效果不理想。湯井田等[10]提出了對視電阻率的重新定義來克服上述困難，其基本原理是建立觀測信號與同一觀測裝置參數(shù)的均勻半空間模型響應(yīng)的關(guān)系方程，再求方程的根，即為所需的視電阻率。這種視電阻率計算方法存在兩方面的問題：①實際觀測的數(shù)據(jù)中不可避免的存在噪聲干擾，即便采用平滑濾波等抑制噪聲的影響依然難以避免視電阻率曲線包含假信息；②由于電磁勘探物理問題的多解性，需要人機交互挑選，不可避免地引入人為干擾。

筆者提出了一種可控源電磁勘探時間域數(shù)據(jù)視電阻率求解新方法，最小泛函擬合視電阻率計算方法，該方法構(gòu)造滿足特定條件目標泛函，通過迭代搜索目標泛函的極小，得到視電阻率光滑曲線函數(shù)。該方法綜合考慮了數(shù)據(jù)噪聲對視電阻率造成的影響，同時也克服了電磁勘探中的多解性問題。

1 最小泛函擬合視電阻率計算方法

傳統(tǒng)的視電阻率計算是依據(jù)均勻半空間中電磁場響應(yīng)與觀察數(shù)據(jù)等式關(guān)系，求解正演公式的根，見式(1)。

d=f(m)±e

(1)

其中：d為實際觀測數(shù)據(jù)；f(m)為均勻半空間的正演公式；e為數(shù)據(jù)誤差。對式(1)求根便得到該數(shù)據(jù)對應(yīng)的視電阻率結(jié)果，見式(2)。

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

(2)

其中：d為所求的視電阻率；f(m)為計算視電阻率的誤差。通過式(1)和式(2)可以發(fā)現(xiàn)，觀測數(shù)據(jù)中的噪聲對視電阻率計算的影響，在強噪聲環(huán)境中，對視電阻率結(jié)果造成不可預(yù)測的結(jié)果。在收發(fā)距較小，大地電阻率較高或工作頻率較低時，觀測的電磁場屬于近場區(qū)或過渡區(qū)場，在此情況下，視電阻率將發(fā)生畸變，不能直觀反映地下介質(zhì)的電性分布，不利于資料地處理解釋。傳統(tǒng)的視電阻率計算方法沒有整體考慮數(shù)據(jù)擬合，造成視電阻率曲線存在“毛刺”，使得后續(xù)反演出現(xiàn)假異常的情況。

根據(jù)地球物理電磁勘探的原理，視電阻率隨時間變化的函數(shù)是一個光滑的函數(shù)，即視電阻率隨時間的函數(shù)的二階導(dǎo)數(shù)積分最小，如式(3)所示。

(3)

其中：t0、tn為野外觀測的起止時刻。

同時，要求正演數(shù)據(jù)和野外觀測數(shù)據(jù)的擬合度小于等于給定的門限，如式(4)所示。

(4)

其中：f(ma(ti))為均勻半空間中電磁時間域正演計算公式；S2為給定的擬合度門限值。

引入拉格朗日乘子P，構(gòu)造一個目標泛函U，如式(5)所示。

(5)

引入二階導(dǎo)數(shù)算法，式(5)可以改寫為式(6)。

(6)

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

(7)

具體算法的流程如下：

1)設(shè)置當(dāng)前迭代次數(shù)為k，如果k=1，初始化視電阻率；否則，將上次迭代的視電阻率結(jié)果作為初始值。計算初始數(shù)據(jù)擬合度，記為rms0; 計算初始視電阻率粗糙度，記為ruf0；計算Jacobian矩陣，Jacobian矩陣(具體步驟引用文獻[12])。設(shè)置拉格朗日乘子的初始值，一般令p=105。令搜索步長β=1。

3)令m等于當(dāng)前新的視電阻率計算m的數(shù)據(jù)擬合，記為rms1；計算當(dāng)前新的視電阻率m的粗糙度，記為ruf1。

4)如果rms1小于等于給定的擬合度門限，轉(zhuǎn)步驟5)；否則，如果rms1小于rms0，轉(zhuǎn)步驟6)，否則，令k=k+1，轉(zhuǎn)步驟2)。

5)搜索一個最優(yōu)的拉格朗日乘子P，使得rms1等于S2，且P盡可能的大。

6)如果ruf1大于ruf0，令k=k+1，轉(zhuǎn)步驟2)。

7)本次迭代結(jié)束，m即為最優(yōu)的視電阻率。令k=k+1，轉(zhuǎn)步驟1)。

圖2 基于最小泛涵的視電阻率計算流程圖Fig.2 Based on the minimum the apparent resistivity calculation flow chart of the culvert

2 均勻半空間時間域電磁場正演

可控源電磁方法一般常采用接地電偶極源，筆者著重討論水平電偶極源的均與半空間響應(yīng)的計算。納比吉安[13]采用TE分量、TM分量分解的方法，推導(dǎo)了電偶極源的電磁場表達式，形如式(8)、式(9)。

(8)

(9)

由于實際施工中采用的接地電源是長電流源，可以看做多個小電偶極源的線積分。為了兼顧計算精度和速度，可以采用高斯積分公式，模擬野外長電流源的電磁響應(yīng)。

3 野外數(shù)據(jù)測試

我們選取了某工區(qū)一條測線的可控源電磁實際施工采集數(shù)據(jù)進行測試。反演方法采用Constable[12]提出的Occam反演方法，Occam方法是一種的正則化反演方法，對初始模型要求不高。圖3為傳統(tǒng)視電阻率算法結(jié)果的Occam反演剖面，從圖3中可以看出，傳統(tǒng)方法求得的視電阻率的Occam反演剖面存在很多“掛面條”現(xiàn)象，這是因為實際采集數(shù)據(jù)中存在噪聲，而傳統(tǒng)的視電阻率求解方法沒有考慮到噪聲的影響，從而對視電阻率數(shù)據(jù)造成干擾。圖4是最小泛函視電阻率算法結(jié)果的反演剖面與工區(qū)內(nèi)地震、測井推斷的地質(zhì)信息疊合圖。由圖4可以發(fā)現(xiàn)，反演剖面的橫向連續(xù)性得到改善，沒有明顯的“掛面條”現(xiàn)象，說明筆者提出的視電阻率求解方法較好抑制原始數(shù)據(jù)中的噪聲對視電阻率的影響，且19km處的斷裂與已知吻合，而在傳統(tǒng)的視電阻率算法結(jié)果的反演剖面中，該斷裂沒有得到較好地反映，主要是因為受到鄰近測點的從淺至深的高阻條帶的影響。

4 結(jié)論

通過測試，我們發(fā)現(xiàn)最小泛函視電阻率算法，通過整體擬合，考慮各個測點數(shù)據(jù)的噪聲情況，有選擇地保留有效信息、壓制干擾，從而最大限度獲取地下介質(zhì)的電阻率分布情況。并取得了以下幾點認識：

1)筆者提出的方法對數(shù)據(jù)中的噪聲有一定地抑制作用，在求解方法中將數(shù)據(jù)中的噪聲作為數(shù)據(jù)的權(quán)重(噪聲越大，權(quán)重越小)，一定程度上克服噪聲對視電阻率影響。

圖3 某工區(qū)一條測線的傳統(tǒng)視電阻率數(shù)據(jù)的反演結(jié)果Fig.3 The traditional apparent resistivity data inversion results of a line in a work area

圖4 某工區(qū)一條測線的最小泛函視電阻率數(shù)據(jù)的視模反演結(jié)果與地震、測井?dāng)?shù)據(jù)推斷的地質(zhì)信息對比情況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)因為該方法采用迭代搜索，無論從什么初始條件出發(fā)，總能搜索到是觀測數(shù)據(jù)擬合誤差達到最小的視電阻率，即該方法具有一定的數(shù)值穩(wěn)定性。

3)方法暫未考慮地形起伏的影響，在實際工作中，如果工區(qū)地形起伏較大時，需要用二維的方法進行處理。

4)可以利用現(xiàn)有的MT反演程序，避免全區(qū)資料反演方法計算量大的缺點，可以實現(xiàn)野外現(xiàn)場反演。

[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].長沙:中南工業(yè)大學(xué)出版社，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一維全區(qū)反演對比研究[J]. 吉林大學(xué)學(xué)報，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] 何興梅，胡祥云，陳玉萍，等. 奧克姆一維反演的應(yīng)用[J]. 工程地球物理學(xué)報，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一維最小構(gòu)造反演[J]. 物探化探計算技術(shù)，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資料近場校正方法，勘查地球物理勘查地球化學(xué)文集第20集[M]. 北京:地質(zhì)出版社，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]湯井田，何繼善. 水平電偶源頻率域測深中全區(qū)視電阻率定義的新方法[J]. 地球物理學(xué)報，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]湯井田，何繼善，水平多層介質(zhì)上水平電偶源頻率域電磁測深的阻抗實部等效電阻率[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].北京：地質(zhì)出版社，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

江蘇省地質(zhì)勘查專項資金項目(蘇財建[2016]140)

楊輪凱(1971-)，男，高級工程師，主要從事電磁法及應(yīng)用研究，E-mail：584945583@qq.com。

胡建(1980-)，男，高級工程師，主要從事礦產(chǎn)地質(zhì)調(diào)查工作，E-mail：65153743@qq.com。

1001-1749(2017)03-0301-05

P 631.4

A

10.3969/j.issn.1001-1749.2017.03.01