張艷龍，吳燕清，鄧劉洋，崔少北

(1.重慶大學 煤礦災害動力學與控制國家重點實驗室，重慶 400044； 2.中煤科工集團重慶研究院有限公司，重慶 400037)

瞬變電磁法通過不接地回線以脈沖電流為場源，激勵探測體感應二次電流，在脈沖間歇期間測量二次場隨時間變化的響應[1-2]。采集并處理該響應信號，根據衰減異常信息，分析目標體電性特征[3]。

采用大定源回線裝置測量響應信號是常用的工作方式，受限于數據源與處理方式，在實際應用中多以一維反演及相對視電阻率成像為主要成果體現。要實現三維正反演并成熟地在地面勘查中應用仍需進行大量的研究工作[4-6]。頻域電磁法三維正反演技術日趨成熟[7-8]，如何通過數據轉換進行三維正反演已成為熱點問題。以傅里葉變換為基礎進行數據轉換[9]，形成了時頻轉換理論，為瞬變電磁數據的三維正反演提供了一種新方法[10-11]。

基于時頻轉換理論，筆者通過正則化技術，結合有限差分法與NLCG反演算法，在采用大定源回線裝置進行瞬變電磁探測時，進行了理想模型及實測瞬變電磁數據的頻域正反演研究，論證了三維反演的可行性，同時驗證了其高分辨能力。

1 時頻轉換

從時域信號中獲取頻域響應時，時域信號需滿足采樣定律。只有具備時域和頻域電磁響應等價條件，才能實現相互轉化。瞬變電磁儀器是依據時間進行數據采集的，其采樣范圍有限，數據不滿足傅里葉變換要求，“混頻效應”會使頻域響應的可信度降低。為解決以上問題，直接對傅里葉變換分段積分進行離散，得到矩陣形式。采用最小二乘法，并對頻域響應加入光滑約束，完成時域到頻域的轉換。

1.1 關斷時間校正

時頻轉換理論以階躍發(fā)射波形為基礎，而野外實測的瞬變電磁數據的發(fā)射波通常是具有一定關斷時間的斜階躍波，因此在對時域信號進行頻域轉換之前要對關斷時間進行校正，得到理想階躍關斷條件下的瞬變電磁數據。關斷時間通過最小二乘擬合后沿估計矯正。

1.2 時頻轉換矩陣計算方程

瞬變電磁的時、頻域響應可以相互轉換，其具體步驟為：傅里葉變換→歐拉變換→折線逼近法離散→矩陣變換→正則化求解。

傅里葉變換：

(1)

式中：h(t)為時域電磁響應；t為時間；I(ω)為電流波形的頻譜；ω為角頻率，ω=2πf；H(ω)為頻域電磁響應。

在瞬變電磁法中，最常用的是下階躍波形，故式(1)可變換為式(2)：

(2)

通過歐拉公式，推導頻域—時域電磁場的正弦關系：

(3)

式中：Re為取實符號，表示頻率域諧變電磁場的實分量；h(0)為零時刻磁場的初始值[12]。

對式(3)用折線逼近法進行離散[13]得到式(4)：

[sin(ωi+1t)-sin(ωit)]

(4)

式中：Fs(ω)=Re[Hz(ω)]/ω；Hz(ω)為頻域電磁響應。

將式(4)用矩陣形式表達，令N離散的頻率響應Fs構成矢量D，由M個測道組成的時域響應矢量為d，對于第ti道響應di有：

(5)

令系數：

進一步令：

(6)

則式(5)可寫為：

用矩陣表示為：

d=LD

(7)

式中：Dj為第j個離散頻率響應構成的矢量；L為系數矩陣。

從式(5)和式(6)可以推出矩陣系數L：

(8)

其中矩陣M可以通過式(6)求出，只與余弦和角頻率有關。

1.3 正則化技術

為了提高分辨率，采用有利于解決不適定問題且收斂性更好的正則化反演方法[13]。由于矩陣L的奇異特性導致不能直接求解式(5)，解決該問題的經典方案是采用正則化技術，通過對解加以適當約束以得到穩(wěn)定的解[14]，將矩陣方程計算轉化為最優(yōu)正則化問題，見式(9):

(9)

式中:WD為模型光滑矩陣;λ為正則化參數。

最優(yōu)正則化參數的選擇方案很多，這里采用的是L-曲線法求取最優(yōu)正則化參數。L-曲線法是一種基于數據誤差水平未知的啟發(fā)式選取正則化因子λ的方法[15]。

1.4 數據計算驗證

為驗證時頻轉換的精度與可行性，進行三維理論模型正演數據及工程勘察數據的轉換，對比響應特征。

1)三維模型試驗

在一個電阻率為150 Ω·m均勻半空間，發(fā)射回線邊長分別為400、400 m,在框內中心位置下方 100 m 安設一個三維高阻異常體，異常體尺寸為 80 m ×80 m×80 m，發(fā)射電流3 A,三維模型布設如圖1 所示。

圖1 三維模型布設平面圖

圖1中外部虛線框為發(fā)射回線，中部實線框為接收范圍，內部有剖面線的線框為高阻異常體邊界，中心圓點為坐標原點，黑色圓點為瞬變電磁測點，線距為40 m,點距為20 m。

選取框內A、B兩個測點的數據進行時頻轉換對比，其中A點位于高阻異常體附近，B點位于高阻異常體較遠處。把在A、B點采集的時域信號通過時頻轉換矩陣計算方程進行時頻轉換、對比。三維模型數據時頻轉換對比如圖2所示。

(a)A點時域信號圖

(b)A點垂直磁場響應和時頻轉換對比圖

(c)B點時域信號圖

(d)B點垂直磁場響應和時頻轉換對比圖

圖2中A點、B點轉換計算的正則化因子λ分別為1.15×10-4、2.97×10-3。

數值試驗討論了時域到頻率的計算效果。對時域垂直磁場響應進行了計算，利用時頻轉換方法，將A、B兩個測點的瞬變電磁數據轉化為頻域響應，由圖2可知，不同形態(tài)、反映數據的時頻轉換曲線與理論值基本一致，表明利用正則化技術的時頻轉換理論是可行的，變換結果的正確性證明了時頻轉換方法的可靠性。

2)實測數據試驗

本次瞬變電磁法試驗采用40 m×20 m的網格，線距為40 m，點距為20 m。使用V8多功能電法儀中心回線裝置，發(fā)射線框尺寸為280 m×280 m，發(fā)射電流為10 A，采樣頻率為25 Hz，單點采樣時間 2 min，增益×4。接收區(qū)域位于發(fā)射線框中心120 m×120 m范圍。選取試驗點為L4測線160 m里程點，標記為實測點C。

實測數據時頻轉換對比如圖3所示。

(a)感應電動勢與時域磁場圖

(b)感應電動勢及其擬合曲線圖

(c)感應電動勢OCCAM反演模型圖

(d)頻域磁場對比圖

圖3中轉換計算的正則化因子λ為4.87×103。

圖3給出了28個時間窗口的感應電動勢值，從圖中可以看出實測曲線比較光滑，無需進行濾波等預處理，首個和最后一個時間道分別為2.320×10-5s和1.002×10-2s，在轉換之前，得到了理想階躍關斷條件下瞬變電磁響應。對關斷時間校正數據進行積分，得到時間域垂直磁場，原始數據和積分后的垂直磁場如圖3(a)所示，擬合差小于4%。反演后得到層狀模型如圖3(c)所示，應用時頻轉換方法得到的頻率域垂直磁場響應與圖3(c)模型正演得到的頻率域垂直磁場數據對比如圖3(d)所示。

通過對圖3 分析可知，時頻轉換曲線與理論值基本一致。理論模型數據及實測數據試驗證明了本文采用的時頻轉換方法的有效性、穩(wěn)定性及精確度，為瞬變電磁法由時域向頻域轉換，以及在頻域的三維正反演提供了依據。

2 三維模型正演

三維瞬變電磁正演研究及算法較多，如積分方程法[16-17]、時域有限差分等，頻域響應計算后轉換成時域[18-20]。為驗證時頻轉換技術的可行性，建立正演模型，開展正演研究，以期得到理論模型下的正演成果，為下一步反演提供數據。

2.1 基本原理

基于麥克斯韋方程組，地下空間內的電場、磁場信號可被視為離散參數，將三維模型解構為大型幾何空間集合體，每一組響應構成一個單元。求解幾何單元外棱交點上的電場信號，通過非線性插值方法及電磁感應原理，可求得模型中任何單點的頻域電場、磁場信號參數，最后求得三維模型的電磁場響應。

基于有限差分法實現三維瞬變電磁正演時還需注意幾個問題：

1)磁場的散度：約束計算與求解，約束條件—磁場散度等于0；

2)網格剖分：模型剖分、單元劃分時，網格設置應滿足狄利克雷邊界條件。

2.2 三維模型驗證

建立電阻率為150 Ω·m的均勻半空間，在深度100 m安設1個5 Ω·m三維低阻異常體，200 m的正方形發(fā)射線框敷設于水平地面，發(fā)射電流為3 A，設三維低阻異常體的中心點在水平地表的投影為三維坐標原點，z軸方向為沿地表垂直向下，主剖面為沿x軸的測線剖面。三維模型布設與觀測方式如圖4 所示。

圖4 三維模型布設與觀測示意圖

設置三維網格剖分參數如下：沿坐標軸方向 50 m×50 m×30 m，空間網格大小30 m×30 m×30 m。主剖面上0.01、20 ms的正演數據時域響應(方形點)如圖5所示。使用積分方程方法求得的時域信號用曲線繪制[21]。

圖5 三維模型正演數據時域響應曲線

由圖5可知，該模型的正演試驗表明，基于有限差分法的三維模型正演數據時域響應與Newman等[12]的結論基本一致，證明筆者采用的三維正演算法可靠并具有高精度，能為三維模型反演的研究提供依據。

3 三維模型反演

頻域電磁法數據三維反演主要有擬牛頓法、高斯牛頓法、共軛梯度法和非線性共軛梯度法(NLCG)[22-24]。NLCG算法是非線性反演中確定性、統(tǒng)計性綜合的方法，可將復雜計算線性化。通過精密的網格剖分與地質建模，基于函數的梯度，不斷迭代計算，求得三維模型最優(yōu)解。

在視電阻率150 Ω·m的均勻半空間中設計2個10 Ω·m的低阻異常體，埋深為100 m，其三維形態(tài)為100 m正方體，地面敷設2個400 m方形發(fā)射線框，發(fā)射電流為3 A，測區(qū)邊長為280 m，測線共8條，每條測線布置10個測點，測框布置如圖6所示。

圖6 測框布置示意圖

三維模型反演網格剖分密度為30 m×25 m×20 m，共計反演迭代計算120次，三維模型反演成果切片如圖7所示。

圖7 三維模型反演成果切片示意圖

由圖7可知，三維模型反演成果的低阻異常體范圍與理論模型高度一致，表明NLCG反演算法可實現1條測線上不同發(fā)射源、不同測區(qū)的數據反演，同時證明NLCG反演算法在三維模型反演時的可行性并驗證了其精度。

4 工程數據三維反演

4.1 地質概況

勘探井田位于太行山脈西部沁水盆地南部邊緣，為低山丘陵地帶。井田內地表大部分被第四系黃土覆蓋，局部基巖出露。石炭系上統(tǒng)太原組與二疊系下統(tǒng)山西組為本勘探區(qū)主要含煤地層。

根據勘探區(qū)內鉆孔地質信息，得到該勘探區(qū)煤、巖層電性特征，如表1所示。

表1 地層電性特征

上覆風化砂泥巖，二疊系分布有砂巖、泥巖、煤層；石炭系分布有石灰?guī)r、煤層、砂巖、鐵(鋁)巖；奧陶系以泥灰?guī)r、石灰?guī)r為主。

由表1可以看出，該區(qū)地層電性特征有較大區(qū)別，上覆風化層及泥巖視電阻率值較低；二疊系底部、石炭系頂部煤層地層表現為較高視電阻率；奧陶系又呈高阻反映。表明勘探區(qū)是典型的三層H 模型。

4.2 數據處理

瞬變電磁法試驗采用40 m×20 m的網格，線距為40 m，點距為20 m。使用V8多功能電法儀中心回線裝置，發(fā)射線框尺寸為280 m×280 m，發(fā)射電流為10 A，采樣頻率25 Hz，單點采樣時間為2 min，增益×4。接收測框尺寸為120 m×120 m，接收面積小于1/3發(fā)射面積，測區(qū)內共布設南北向12條測線，工程布置如圖8所示。

圖8 工程布置圖

對在第5線300 m里程點采集的數據進行單點反演，該試驗點衰減曲線如圖9所示(B為磁感應強度)。結合附近鉆孔資料，此處3#煤層埋深為98 m。

圖9 試驗點衰減曲線

試驗點一維反演成果如圖10所示。

圖10 試驗點一維反演成果

由圖10可知，隨著埋藏深度增加二疊系(P)地層視電阻率逐步升高，進入石炭系(C)地層后視電阻率進一步升高，到達奧陶系(O)后視電阻率依然表現為高電阻率值特征。與鉆孔資料對比可知，反演結果基本反映出了鉆孔附近地層的電性特征。

試驗點數據表明,該次瞬變電磁探測參數合理、數據質量較高。將原始數據濾波去噪及關斷時間校正后，利用正則化技術進行了本區(qū)數據的時頻轉換。

4.3 頻域三維反演

經數據處理得到頻域數據，在反演前，結合地質信息建立約束模型。本次反演依據鉆孔地層信息進行物性層狀建模，深度方向設置不同物性層狀巖層共計15層，網格剖分密度為20 m×20 m×20 m，邊界擴充系數為1.5，滿足狄利克雷邊界條件。全部433個測點數據均參與三維反演，迭代擬合差RMS=1.25，反演深度為400 m。三維反演垂直切片成果如圖11 所示。

圖11 三維反演垂直切片圖

三維反演將多測區(qū)、多測線數據進行聯(lián)合反演，豐富了地下空間信息，沿煤層走向上進行密集垂直切片，切片間距為15 m(傳統(tǒng)一維反演只能沿測線生成單一縱向斷面圖)，編號102至107。

由圖11可以看出，該區(qū)域淺部存在風化砂巖相對高阻層，在煤層附近存在泥巖及二疊系砂巖(裂隙含水)相對低阻層，符合H模型與區(qū)域地質特征；低阻異常在空間上的變化趨勢為沿編號從小至大方向上由淺變深、由兩側向中部收斂，低阻異常體空間分布與形態(tài)更直觀。

5 結論

1)利用正則化技術解決時頻轉換過程中矩陣奇異性求解的問題，實現了瞬變電磁有限時域響應向頻率響應的轉換。

2)基于有限差分法實現了瞬變電磁頻域的三維正演，為三維正反演方法提供了參考。

3)通過理論模型三維反演，驗證了NLCG算法在瞬變頻域三維反演的可行性與精度。

4)在實際工程中三維反演多切片成果更直觀，地質解釋更全面。