層狀大地表面中心回線瞬變電磁響應(yīng)特征

吳 瓊， 李永博， 李 貅， 晉 達

(1.中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，廊坊 065000；2.長安大學(xué) 地質(zhì)工程與測繪學(xué)院，西安 710054；3.中石化中原油田分公司 物探研究院，濮陽 457001)

?

層狀大地表面中心回線瞬變電磁響應(yīng)特征

吳 瓊1， 李永博1， 李 貅2， 晉 達3

(1.中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所，廊坊 065000；2.長安大學(xué) 地質(zhì)工程與測繪學(xué)院，西安 710054；3.中石化中原油田分公司 物探研究院，濮陽 457001)

中心回線裝置下水平地層的時間域電磁場表達式為雙重積分，內(nèi)層為漢克爾型積分，外層為正弦或余弦積分，線性數(shù)字濾波法是求解此類特殊積分的有效方法。這里采用線性數(shù)字濾波法分別計算了均勻半空間、兩層和三層介質(zhì)共七種地電斷面的電磁響應(yīng)，分析了不同類型響應(yīng)曲線的形態(tài)特征，研究了中間層厚度變化和底層電阻率變化對時間域響應(yīng)曲線形態(tài)的影響。通過對比與分析，得出了一些規(guī)律性的認(rèn)識，對于瞬變電磁數(shù)據(jù)的處理解釋具有一定的指導(dǎo)意義。

漢克爾積分變換； 正弦變換； 余弦變換； 線性數(shù)字濾波； 瞬變電磁響應(yīng)

0 引 言

瞬變電磁法具有工作效率高、電性分辨力強、探測深度大、抗干擾能力強等顯著優(yōu)勢,因而被廣泛應(yīng)用于礦產(chǎn)勘探、工程勘察、環(huán)境災(zāi)害地質(zhì)調(diào)查等多個領(lǐng)域[1-4]。時間域瞬變電磁響應(yīng)的求取通常有兩種方案：①直接法，根據(jù)一定的初始條件和邊界條件在時間域中直接求解；②間接法，即頻-時轉(zhuǎn)換，根據(jù)頻譜分析理論，首先在頻率域求解給定場源的電磁響應(yīng)，再通過傅里葉反變換求得相應(yīng)的時間域電磁響應(yīng)[5]。因頻率域電磁場表達式較為簡單，故其應(yīng)用更為普遍。

實現(xiàn)頻-時轉(zhuǎn)換的主要方法有離散傅里葉變換法、延遲譜法、逆拉普拉斯變換法、數(shù)字濾波法等[5]，離散傅里葉變換法因其所需的頻率和核函數(shù)抽樣次數(shù)多、計算量大、速度慢而極少被采用，延遲譜法的晚期響應(yīng)不太穩(wěn)定，逆拉普拉斯變換法的計算量大、計算速度較慢，數(shù)字濾波法相比前三種方法，在層狀大地情況下具有計算精度高、速度快且晚期響應(yīng)穩(wěn)定的優(yōu)點，是當(dāng)前較常用的一種計算方法。

在線性數(shù)字濾波算法方面，諸多學(xué)者開展了系列的研究工作，D.Guptasarma[6]提出了一種用于計算時間域激電響應(yīng)的線性數(shù)字濾波算法，并給出了21點濾波系數(shù)；Aderson[7]提出了滯后褶積快速濾波算法；王華軍[8]通過快速漢克爾變換理論導(dǎo)出了正弦變換和余弦變換數(shù)字濾波算法，給出了兩組250點濾波系數(shù)，并驗證了這組濾波算法具有計算速度快、精度高的優(yōu)點[9]。

這里采用250點正弦變換濾波算法，分別計算了一層(均勻半空間)、兩層和三層水平介質(zhì)的時間域電磁響應(yīng)，總結(jié)了典型地電斷面的電磁響應(yīng)曲線特征，并研究了中間層厚度變化和底層電阻率變化對時間域響應(yīng)曲線形態(tài)的影響，得到了一些規(guī)律性的認(rèn)識。

1 理論基礎(chǔ)

一個半徑為a的圓形線圈水平放置在水平地層表面，通入電流的大小為I0，角頻率為ω，磁導(dǎo)率為μ。水平地層的電導(dǎo)率由上到下依次為σ1、σ2、…、σn-1、σn，各地層厚度依次為H1、H2、…、Hn-1、∞。采用柱坐標(biāo)系，設(shè)原點在圓回線中心點，取z軸向下為正。當(dāng)z=0時，接收線圈位于地表，此時各電磁場分量在頻率域的積分表達式為[1]：

(1)

(2)

(3)

當(dāng)忽略位移電流，并在中心點接收時：

Eφ=Hr=0

(4)

(5)

公式(5)為各向同性水平地層表面圓回線源中心點的頻率域電磁場表達式，在回線中心點只能觀測到磁場的垂直分量。在均勻半空間這一特殊情況下，回線源中心點頻率域垂直磁場的積分表達式可化簡為式(6)。

(6)

諧變電流條件下的電磁場量Hz(ω)與階躍電流條件下的電磁場量Hz(t)滿足如下關(guān)系：

(7)

利用歐拉公式、Hz(ω)的復(fù)數(shù)性質(zhì)、δ(ω)函數(shù)的積分性質(zhì)、積分區(qū)間的對稱性、ReHz(ω)和lmHz(ω)的奇偶性及其頻率特性等[1]，可得：

(8)

其相應(yīng)的時間導(dǎo)數(shù)為：

(9)

文中水平層狀介質(zhì)的時間域響應(yīng)，采用250點正弦變換濾波算法,計算公式為[9-10]式(10)。

(10)

2 模型計算

2.1 均勻半空間

為研究時間域響應(yīng)曲線隨回線半徑的變化規(guī)律，令I(lǐng)0=1A、ρ1=1 Ω·m，回線半徑a分別取五個不同的值，求得的時間域響應(yīng)曲線如圖1所示。

圖1 均勻半空間時間域響應(yīng)(不同回線半徑)

為研究時間域響應(yīng)曲線隨電阻率的變化規(guī)律，令I(lǐng)0=1 A、a=100 m，均勻半空間電阻率ρ1分別取五個不同的值，求得的時間域響應(yīng)曲線如圖2所示。

圖2 均勻半空間時間域響應(yīng)(不同電阻率)

均勻半空間時間域響應(yīng)曲線具有如下規(guī)律性特征：在雙對數(shù)坐標(biāo)下，早期響應(yīng)曲線近乎水平，晚期響應(yīng)曲線為向下傾斜的直線，與時間軸的夾角約為68.2°。回線源半徑越小，在中心點的早期響應(yīng)幅值越大，越早進入晚期，均勻半空間電阻率越小，在中心點的早期響應(yīng)幅值越小，越晚進入晚期。

2.2 兩層模型

兩層模型分為G型和D型兩種地電斷面，其參數(shù)取值如圖3所示，線性數(shù)字濾波法求得的時間域響應(yīng)見圖4。

圖3 兩層模型示意圖

圖4 兩層模型時間域響應(yīng)

從時間域響應(yīng)曲線形態(tài)來看，三者在早期階段基本重合，在雙對數(shù)坐標(biāo)下，為近水平的直線，隨著時間的推移，G型逐漸向下偏離均勻半空間響應(yīng)曲線，D型先向下后向上偏離均勻半空間響應(yīng)曲線。由于瞬變電磁場服從熱傳導(dǎo)方程的規(guī)律，隨著時間的增加，場向深處傳播過程中以“煙圈”形式逐漸向外擴散，時間域響應(yīng)逐漸衰減，而衰減曲線的斜率對應(yīng)著地層的信息，G型的第二層電阻率大于第一層電阻率，故時間域響應(yīng)衰減速率增大，表現(xiàn)為時間域響應(yīng)曲線在均勻半空間響應(yīng)曲線的下方，D型的第二層電阻率小于第一層電阻率，故時間域響應(yīng)衰減速率減小，表現(xiàn)為時間域響應(yīng)曲線在均勻半空間響應(yīng)曲線的上方。響應(yīng)曲線的兩種不同形態(tài)，是由于層狀模型高電阻率目標(biāo)或低電阻率目標(biāo)的電磁場衰減特性引起的。由此可見，瞬變電磁響應(yīng)曲線反映了地下介質(zhì)的電性變化特征，即當(dāng)電磁場由低阻層進入高阻層時衰減加快，當(dāng)電磁場由高阻層進入低阻層時衰減放慢。

2.3 三層模型

三層模型分為A型、K型、H型和Q型四種地電斷面，其參數(shù)取值如圖5所示。線性數(shù)字濾波法求得的時間域響應(yīng)見圖6。

圖5 三層模型示意圖

圖6 三層模型時間域響應(yīng)

從時間域響應(yīng)曲線形態(tài)來看，五者形態(tài)各異。在雙對數(shù)坐標(biāo)下，隨著時間的推移，A型逐漸向下偏離均勻半空間響應(yīng)曲線，在中晚期衰減變快；K型和均勻半空間響應(yīng)曲線首尾重合，在中期衰減稍快；H型圍繞均勻半空間響應(yīng)曲線先“下凹”后“上凸”，在晚期逐漸向均勻半空間響應(yīng)曲線靠近；Q型呈“下凹”、“上凸”、“下凹”、“上凸”規(guī)律性變化，尤其在晚期，曲線衰減極為緩慢。總體來看，A型是GG型組合，H型是DG型組合，K型是GD型組合，Q型是DD型組合，因此在中晚期階段，曲線的形態(tài)變化仍是與電磁場是由低阻進入高阻還是由高阻進入低阻兩種不同電性界面性質(zhì)有關(guān)。

2.4 改變中間層厚度模型

令G型、D型地電斷面的第一層厚度與A型、H型、K型、Q型地電斷面的第二層厚度分別取25m、50m、100m，同類型模型的其他參數(shù)取值保持固定，計算出的時間域響應(yīng)曲線如圖7所示。

由圖7(a)、7(b)可見，G型、D型地電斷面的時間域響應(yīng)曲線在早期與均勻半空間響應(yīng)曲線基本重合，隨著時間的推移，G型曲線逐漸向下偏離均勻半空間響應(yīng)曲線，D型曲線先向下后向上偏離均勻半空間響應(yīng)曲線，G型和D型的第一層厚度越小，曲線偏離越明顯，在晚期階段，三種不同第一層厚度的G型和D型時間域響應(yīng)曲線均具有重合趨勢。由圖7(c)、7(e)可見，A型、K型地電斷面的時間域響應(yīng)曲線在早期與具有相同地電參數(shù)的G型時間域響應(yīng)曲線基本重合，隨著時間的推移，A型曲線逐漸向下偏離G型響應(yīng)曲線，K型曲線逐漸向上偏離G型響應(yīng)曲線，不同第二層厚度的A型和K型時間域響應(yīng)曲線均比較接近，很難分辨。由圖7(d)、7(f)可見，H型、Q型地電斷面的時間域響應(yīng)曲線在早期與具有相同地電參數(shù)的D型時間域響應(yīng)曲線基本重合，隨著時間的推移，H型曲線逐漸向下偏離D型響應(yīng)曲線，Q型曲線先向下后向上偏離D型響應(yīng)曲線，H型和Q型的第二層厚度越小，曲線偏離越明顯，在晚期階段，三種不同第二層厚度的H型和Q型時間域響應(yīng)曲線均具有重合趨勢。

2.5 改變底層電阻率模型

令G型、D型、A型、K型、H型、Q型六種地電斷面的底層電阻率分別取3個不同的值，同類型模型的其他參數(shù)取值保持固定，計算出的時間域響應(yīng)曲線如圖8所示。

由圖8(a)可見，第二層電阻率越大，G型曲線向下偏離均勻半空間響應(yīng)曲線越明顯，但當(dāng)?shù)诙与娮杪试龃蟮揭欢ǔ潭群螅淙≈档牟煌跁r間域響應(yīng)曲線中難以辨別。由圖8(b)可見，第二層電阻率越小，D型曲線先向下后向上偏離均勻半空間響應(yīng)曲線幅度越大，在晚期階段，第二層電阻率不同的三種D型時間域響應(yīng)呈直線形態(tài)，且平行于均勻半空間情況下的晚期響應(yīng)。由圖8(c)可見，第三層電阻率不同的三種A型時間域響應(yīng)曲線基本重合，無明顯差異。由圖8(d)可見，第三層電阻率越大，H型曲線向下偏離D型響應(yīng)曲線越明顯，但當(dāng)?shù)谌龑与娮杪试龃蟮揭欢ǔ潭群螅淙≈档牟煌跁r間域響應(yīng)曲線上反映得不再明顯。由圖8(e)可見，第三層電阻率越小，K型曲線偏離G型響應(yīng)曲線越明顯，在晚期階段，第三層電阻率不同的三種K型時間域響應(yīng)曲線幾乎平行。由圖8(f)可見，第三層電阻率越小，Q型曲線先向下后向上偏離D型響應(yīng)曲線幅度越大，第三層電阻率不同的三種Q型時間域響應(yīng)曲線差異明顯，并在晚期階段，呈相互平行的直線。

圖7 改變中間層厚度模型時間域響應(yīng)

圖8 改變底層電阻率模型時間域響應(yīng)

3 結(jié) 論

通過對均勻半空間、兩層和三層介質(zhì)共七種地電斷面的電磁響應(yīng)曲線進行對比分析，得到以下幾點認(rèn)識和結(jié)論：

1)在雙對數(shù)坐標(biāo)下，早期的時間域響應(yīng)曲線為平行于時間軸的直線；晚期的時間域響應(yīng)曲線為向下傾斜的直線，它與時間軸的夾角約為68.2°。

2)回線源半徑越小，在中心點的早期響應(yīng)幅值越大，越早進入晚期；地層電阻率越小，在中心點的早期響應(yīng)幅值越小，越晚進入晚期。

3)當(dāng)由低阻層進入高阻層(G型)時，電磁響應(yīng)衰減加快；當(dāng)由高阻層進入低阻層(D型)時，電磁響應(yīng)衰減放慢。

4)A型與G型地電斷面的電磁響應(yīng)非常接近，不易分辨；K型和均勻半空間響應(yīng)曲線首尾重合，在中期衰減稍快，反映了中間層的高阻特征；H型圍繞均勻半空間響應(yīng)曲線先“下凹”后“上凸”，在中期衰減較慢，反映了中間層的低阻特征；Q型呈“下凹”、“上凸”、“下凹”、“上凸”規(guī)律性變化，尤其在晚期，曲線衰減極為緩慢。時間域電磁法對高阻層的分辨能力遠沒有對低阻層的分辨能力強。

5)G型和D型地電斷面的第一層厚度越小，電磁響應(yīng)曲線偏離均勻半空間響應(yīng)曲線越明顯，在晚期階段，不同第一層厚度的G型和D型時間域響應(yīng)曲線均具有重合趨勢；不同第二層厚度的A型和K型時間域響應(yīng)曲線均比較接近，很難分辨；H型和Q型的第二層厚度越小，曲線偏離具有相同地電參數(shù)的D型時間域響應(yīng)曲線越明顯，在晚期階段，不同第二層厚度的H型和Q型時間域響應(yīng)曲線均具有重合趨勢。

6)G型地電斷面的第二層電阻率越大，其響應(yīng)曲線向下偏離均勻半空間響應(yīng)曲線越明顯，但當(dāng)?shù)诙与娮杪试龃蟮揭欢ǔ潭群螅淙≈档牟煌跁r間域響應(yīng)曲線中難以辨別，在晚期階段，第二層電阻率不同的G型時間域響應(yīng)曲線具有重合趨勢；D型地電斷面的第二層電阻率越小，其響應(yīng)曲線先向下后向上偏離均勻半空間響應(yīng)曲線幅度越大，在晚期階段，第二層電阻率不同的D型時間域響應(yīng)呈直線形態(tài)，且平行于均勻半空間情況下的晚期響應(yīng)。在中晚期階段，A型和H型時間域響應(yīng)曲線隨第三層電阻率的變化規(guī)律與G型時間域響應(yīng)曲線隨第二層電阻率的變化規(guī)律類似，K型和Q型時間域響應(yīng)曲線隨第三層電阻率的變化規(guī)律與D型時間域響應(yīng)曲線隨第二層電阻率的變化規(guī)律類似。

上述規(guī)律性認(rèn)識對于瞬變電磁資料的處理和解釋，具有一定的參考價值。

[1] 方文藻，李予國，李貅. 瞬變電磁測深法原理[M]. 西安: 西北工業(yè)大學(xué)出版社，1993.FANGWZ,LIYG,LIX.Theprincipleoftransientelectromagneticsoundingmethod[M].Xi’an:NorthwesternPolytechnicUniversityPress, 1993. (InChinese)

[2] 樸化榮. 電磁測深原理[M]. 北京: 地質(zhì)出版社，1990.PIAOHR.Principleofelectromagneticsoundings[M].Beijing:GeologicalPublishingHouse, 1990. (InChinese)

[3] 李貅. 瞬變電磁測深的理論與應(yīng)用[M]. 西安: 陜西科學(xué)技術(shù)出版社，2002.LIX.Thetheoryandapplicationoftransientelectromagneticsounding[M].Xi’an:ShaanxiScienceandTechnologyPress, 2002. (InChinese)

[4] 牛之璉. 時間域電磁法原理[M]. 長沙: 中南大學(xué)出版社，2007.NIUZHL.Timedomainelectromagneticmethods[M].Changsha:CentralSouthUniversityPress, 2007. (InChinese)

[5] 陳向斌, 胡社榮, 張超. 瞬變電磁場響應(yīng)計算的頻-時域轉(zhuǎn)換方法綜述[J]. 工程地球物理學(xué)報, 2008, 5(2)：242-246.CHENXB,HUSR,ZHANGCH.Reviewoffrequency-domainandtime-domaintransformationmethodintransientelectromagneticfieldresponsecomputation[J].ChineseJournalofEngineeringGeophysics, 2008, 5(2): 242-246. (InChinese)

[6]GUPTASARMAD.Computationofthetime-domainresponseofapolarizableground[J].Geophysics, 1982, 47: 1574-1576.

[7]ANDERSONWL.FastHankeltransformsusingrelatedandlaggedconvolutions[J].Assn.Comp.Math.,Trans.,Math.Software, 1982, 8: 344-368.

[8]JOHANSENHK,SORENSENK.FastHankeltransforms[J].GeophysicalProspecting, 1979, 27: 876-901.

[9] 王華軍. 正余弦變換的數(shù)值濾波算法[J]. 工程地球物理學(xué)報, 2004, 1(4): 329-335.WANGHJ.Digitalfilteralgorithmofthesineandcosinetransform[J].ChineseJournalofEngineeringGeophysics, 2004, 1(4): 329-335. (InChinese)

[10]吳瓊. 大回線源電磁場正演與波場變換理論研究[D]. 西安: 長安大學(xué)，2012.WUQ.Forwardmodelingoflargelooptransientelectromagneticfieldandwave-fieldtransformtheory’sstudy[D].Xi’an:Chang’anUniversity, 2012. (InChinese)

Response features of central-loop transient electromagnetic field on the layered earth

WU Qiong1, LI Yong-bo1, LI Xiu2， JIN Da3

(1.Institute of Geophysical and Geochemical Exploration, Langfang 065000,China;2. School of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China；3. Geophysical Research Institute, Zhongyuan Oilfield Company, SINOPEC, Puyang 457001, China)

The expressions of time-domain electromagnetic fields with a central-loop on the layered earth are double integrals. The inner integral is Hankel transform while the outer integral is sine transform or cosine transform. Linear digital filtering is an effective method to solve such special integrals. The electromagnetic responses for half space, two-layer and three-layer models have been calculated with linear digital filtering method, and the response curves for the seven types (half space, G, D, A, K, H and Q) have been analyzed, and time-domain responses of layered models with changing interlayer-thickness and underlayer-resistivity have been studied. By comparing and analyzing, some regular cognitions are drawn, which have certain guiding significance for transient electromagnetic data’s processing and interpretation.

Hankel integral transform; sine transform; cosine transform; linear digital filtering; transient electromagnetic response

2015-05-11改回日期：2015-06-08

中央級公益性科研院所基本科研業(yè)務(wù)費項目(AS2013P01)

吳瓊(1986-)，女，碩士，主要從事電磁勘探理論與方法技術(shù)研究，E-mail：wuqiong8668@163.com。

1001-1749(2015)05-0560-06

P 631.3

A

10.3969/j.issn.1001-1749.2015.05.03