層狀大地表面中心回線瞬變電磁響應特征

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

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

?

層狀大地表面中心回線瞬變電磁響應特征

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

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

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

漢克爾積分變換； 正弦變換； 余弦變換； 線性數字濾波； 瞬變電磁響應

0 引 言

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

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

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

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

1 理論基礎

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

(1)

(2)

(3)

當忽略位移電流，并在中心點接收時：

Eφ=Hr=0

(4)

(5)

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

(6)

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

(7)

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

(8)

其相應的時間導數為：

(9)

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

(10)

2 模型計算

2.1 均勻半空間

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

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

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

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

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

2.2 兩層模型

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

圖3 兩層模型示意圖

圖4 兩層模型時間域響應

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

2.3 三層模型

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

圖5 三層模型示意圖

圖6 三層模型時間域響應

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

2.4 改變中間層厚度模型

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

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

2.5 改變底層電阻率模型

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

由圖8(a)可見，第二層電阻率越大，G型曲線向下偏離均勻半空間響應曲線越明顯，但當第二層電阻率增大到一定程度后，其取值的不同在時間域響應曲線中難以辨別。由圖8(b)可見，第二層電阻率越小，D型曲線先向下后向上偏離均勻半空間響應曲線幅度越大，在晚期階段，第二層電阻率不同的三種D型時間域響應呈直線形態，且平行于均勻半空間情況下的晚期響應。由圖8(c)可見，第三層電阻率不同的三種A型時間域響應曲線基本重合，無明顯差異。由圖8(d)可見，第三層電阻率越大，H型曲線向下偏離D型響應曲線越明顯，但當第三層電阻率增大到一定程度后，其取值的不同在時間域響應曲線上反映得不再明顯。由圖8(e)可見，第三層電阻率越小，K型曲線偏離G型響應曲線越明顯，在晚期階段，第三層電阻率不同的三種K型時間域響應曲線幾乎平行。由圖8(f)可見，第三層電阻率越小，Q型曲線先向下后向上偏離D型響應曲線幅度越大，第三層電阻率不同的三種Q型時間域響應曲線差異明顯，并在晚期階段，呈相互平行的直線。

圖7 改變中間層厚度模型時間域響應

圖8 改變底層電阻率模型時間域響應

3 結 論

通過對均勻半空間、兩層和三層介質共七種地電斷面的電磁響應曲線進行對比分析，得到以下幾點認識和結論：

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

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

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

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

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

6)G型地電斷面的第二層電阻率越大，其響應曲線向下偏離均勻半空間響應曲線越明顯，但當第二層電阻率增大到一定程度后，其取值的不同在時間域響應曲線中難以辨別，在晚期階段，第二層電阻率不同的G型時間域響應曲線具有重合趨勢；D型地電斷面的第二層電阻率越小，其響應曲線先向下后向上偏離均勻半空間響應曲線幅度越大，在晚期階段，第二層電阻率不同的D型時間域響應呈直線形態，且平行于均勻半空間情況下的晚期響應。在中晚期階段，A型和H型時間域響應曲線隨第三層電阻率的變化規律與G型時間域響應曲線隨第二層電阻率的變化規律類似，K型和Q型時間域響應曲線隨第三層電阻率的變化規律與D型時間域響應曲線隨第二層電阻率的變化規律類似。

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

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

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

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

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

[5] 陳向斌, 胡社榮, 張超. 瞬變電磁場響應計算的頻-時域轉換方法綜述[J]. 工程地球物理學報, 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] 王華軍. 正余弦變換的數值濾波算法[J]. 工程地球物理學報, 2004, 1(4): 329-335.WANGHJ.Digitalfilteralgorithmofthesineandcosinetransform[J].ChineseJournalofEngineeringGeophysics, 2004, 1(4): 329-335. (InChinese)

[10]吳瓊. 大回線源電磁場正演與波場變換理論研究[D]. 西安: 長安大學，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

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

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

1001-1749(2015)05-0560-06

P 631.3

A

10.3969/j.issn.1001-1749.2015.05.03