汪 凌，張富明，楊海燕，張 華，李 哲，楊夫杰

(1. 東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，330013，南昌；2. 山西省煤炭地質(zhì)物探測繪院，030699，山西，晉中；3. 中國礦業(yè)大學(xué)資源與地球科學(xué)學(xué)院，221116，江蘇，徐州)

0 引言

我國是世界第一產(chǎn)煤及消費大國，自20世紀以來，我國國民經(jīng)濟增長迅速，對煤炭資源的需求量不斷在增加[1]。近年來，時間域電磁探測方法被廣泛地應(yīng)用在工程地質(zhì)勘察和礦井超前地質(zhì)預(yù)報探測領(lǐng)域中[2]。由于有限環(huán)境的限制和目標層(體)深度的不確定性，小回線裝置在礦井巷道、工程隧道等地下有限空間及城市和工程淺層探測中通常被采用。該裝置通常以共軸形式(中心回線)和共面形式(偶極裝置)為主[3-6]。為了增大對地下目標體的探測深度以及探測目標體的響應(yīng)程度，通常采用增大發(fā)射電流或者發(fā)射匝數(shù)來增強發(fā)射磁矩。然而這使得線圈間的互感增大，關(guān)斷時間變長，淺部“盲區(qū)”增大?；诖?，在多匝小回線瞬變電磁裝置基礎(chǔ)上，課題組提出圓錐型場源裝置，如圖1所示。

圖1 圓錐型場源發(fā)射裝置理論模型圖

圓錐型發(fā)射場源模型可以近似地表示為由n個半徑介于r1和r2之間的單匝線圈組成。如圖1所示，單匝線圈內(nèi)電流強度為I，頂、底部線圈之間的垂直距離為D。各匝線圈的半徑及其中心點之間的距離可以分別表示為[7-9]：

(1)

從式(1)中可以看出，ri只受到圓錐型裝置頂?shù)装霃絩1和r2和裝置匝數(shù)n的約束。

經(jīng)研究結(jié)果表明，圓錐型場源裝置的關(guān)斷時間和電感系數(shù)約為多匝小回線的1/8[10]。通過傳統(tǒng)瞬變電磁理論疊加可以得到該裝置下的一維層狀模型瞬變電磁響應(yīng)[11-12]。

目前，關(guān)于圓錐型場源瞬變電磁反演方法研究較少，主要以“煙圈”快速成像方法為主。該方法是一種半定量定性的成像方法，不需要初始模型代入計算[13-14]。近年來，許多研究學(xué)者將其他領(lǐng)域的反演方法運用于圓錐型場源瞬變電磁反演方法上。其中包括基于最小二乘改進的“煙圈”反演法、遺傳算法、粒子群優(yōu)化算法等[15-16]。OCCAM反演方法以其平滑模型思想的優(yōu)勢性被廣泛應(yīng)用于各種勘探方法，包括大地電磁法、CSAMT等[17-18]?；诖?，本文將OCCAM反演方法應(yīng)用于圓錐型場源瞬變電磁法反演中。通過本文研究工作，期望獲得適用于圓錐型場源的反演算法，豐富圓錐型場源的理論基礎(chǔ)。

1 OCCAM反演基本原理

1987年，Constable等人首先提出了OCCAM反演理論。若反演想取得最優(yōu)解，則反演最終模型應(yīng)該盡可能簡單、光滑，模型粗糙度應(yīng)盡可能小[19]。OCCAM反演方法原理簡述如下。

定義一個粗糙度，粗糙度可以表示為：

(2)

式中m為模型電性參數(shù)，一般為電阻率或是電阻率對數(shù)。

本文采用矩陣形式進行表示為：

R1=‖?m‖2

(3)

式中?是N×N的矩陣，定義為：

(4)

正演模型的響應(yīng)表示如下：

dj=Fj[m],j=1,2,...,M

(5)

則目標函數(shù)的擬合差可寫為：

X2=‖Wd-WF[m]‖2

(6)

式中：W=diag{1/σ1,1/σ2,...,1/σm}。

在此基礎(chǔ)上，引入拉格朗日乘子μ-1，構(gòu)成目標函數(shù)：

(7)

在反演迭代過程中，為使目標函數(shù)取最小，通常令?mU等于零，則等式(7)可以寫成下列形式：

μ-1(WJ)TWJm-μ-1(WJ)TWd+?T?m=0

(8)

式中J是M×N階的雅可比矩陣：

J=?mF

(9)

假設(shè)第k次迭代已完成，定義向量：

mK+1(μ)=[μ?T?(WJK)TWJK]-1(WJK)TWdK

(10)

用一系列μ值計算模型mK+1(μ)的真正擬合差：

Xk+1(μ)=‖Wd-WF[mk+1(μ)]‖

(11)

最終找到一個μ值相對最小，擬合差相對最小，光滑度相對高的模型，即是反演最終結(jié)果。

Occam反演方法反演流程如圖2。

圖2 反演流程圖

2 理論模型算例

2.1 一維層狀介質(zhì)模型

層狀地電模型更能檢驗反演算法的有效性，因此分別建立3層(H型和K型)以及4層(KH型)地電模型(見表1)。發(fā)射電流為10 A，圓錐頂?shù)装霃椒謩e為0.5 m和1 m，接收面積為1 m2，匝數(shù)為 11。

表1 模型參數(shù)表

OCCAM反演的初始模型選擇均勻半空間，電阻率取響應(yīng)的平均值，模型的初始厚度為4 m，增加比例為1.1，一共30層。初始拉格朗日乘子設(shè)為0.5，擬合差閾值設(shè)為5%，迭代次數(shù)設(shè)為20次。計算結(jié)果見圖3、圖4、圖5。

(a)反演曲線；(b)反演結(jié)果的正演響應(yīng)擬合曲線

(a)反演曲線；(b)反演結(jié)果的正演響應(yīng)擬合曲線

(a)反演曲線；(b)反演結(jié)果的正演響應(yīng)擬合曲線

為了驗證本方法不依賴初始模型，分別采用不同初始電阻率(1 Ω·m和50 Ω·m)對同一數(shù)據(jù)進行反演，反演結(jié)果如圖6所示。

從圖3～圖6反演結(jié)果可以看出：1)OCCAM反演電阻率更接近各層介質(zhì)的真實電阻率，反演效果較好；2)OCCAM反演深度更接近層界面，對目標層的底界面分辨能力更強，分層效果更加明顯；3)該方法不依賴初始模型，結(jié)果準確。由于初始設(shè)置的反演迭代次數(shù)為20，在實際情況中，當(dāng)?shù)螖?shù)為5次時，就會到達誤差閾值，反演所需時間較少，故反演效率高。

圖6 不同初始電阻率反演結(jié)果對比圖

2.2 擬二維反演

建立如圖7所示的二維理論模型，設(shè)置31個測點，點距為10 m，發(fā)射電流為10 A，圓錐型場源頂?shù)装霃椒謩e為0.5 m和1 m，接收面積為1 m2，匝數(shù)為11。將正演計算數(shù)據(jù)進行OCCAM反演，反演結(jié)果如圖8所示。

圖7 擬二維理論模型圖

圖8 擬二維反演結(jié)果圖

從圖8中可以清晰看出，OCCAM反演方法在擬二維模型中，反演結(jié)果與實際模型對應(yīng)較好，結(jié)果較為準確，進一步證明該方法的有效性和可行性。

3 實測數(shù)據(jù)反演

采用圓錐型場源裝置在東北新肇某地區(qū)開展了瞬變電磁剖面實驗。實驗設(shè)備采用澳大利亞生產(chǎn)的Terra TEM瞬變電磁系統(tǒng)，發(fā)射電流為8 A，發(fā)基頻為50 Hz，采樣時間為1.687—433.487 ms，疊加次數(shù)為32次。圓錐型發(fā)射場源頂、底半徑分別為0.15 m 和0.325 m，高為0.5 m，繞制127 匝，接收線圈有效面積為13.448 m2。由于處于平原地區(qū)，地形起伏的影響可以忽略不計，又前文驗證了不依賴初始模型，故初始模型的建立不影響結(jié)果，反演深度控制在250 m內(nèi)。

為了驗證OCCAM反演方法的實際效果，將瞬變電磁法反演結(jié)果與CSAMT方法進行對比。如圖9所示，其中圖9(a)表示CSAMT反演結(jié)果，圖9(b)表示圓錐型場源瞬變電磁反演結(jié)果。圖9(b)為一段長650 m的剖面反演斷面圖，對應(yīng)圖9(a)中藍色虛線框的位置。

(b)OCCAM反演

從圖9中可以看出，經(jīng)過OCCAM反演方法處理得到的瞬變電磁法反演結(jié)果與CSAMT反演結(jié)果大概一致，地層電阻率差異特征明顯，說明該方法在實測數(shù)據(jù)處理是可行的。研究內(nèi)容為圓錐型場源瞬變電磁法數(shù)據(jù)處理提供新的方法和技術(shù)。

4 結(jié)論

1)圓錐形場源裝置可以等效為不同半徑不同高度線圈疊加而成，其正演響應(yīng)理論可以等效為不同場源半徑的中心回線瞬變電磁響應(yīng)疊加得到。最終通過140點漢克爾變換和250點余弦變換推導(dǎo)出正演響應(yīng)。

2)OCCAM反演方法不依賴初始模型選擇，反演結(jié)果更接近地層模型的真實電阻率。

3)OCCAM反演方法分層效果明顯。經(jīng)過實測數(shù)據(jù)處理表明，該方法可以作為圓錐型場源瞬變電磁反演方法的新手段。

致謝：重慶大學(xué)電氣工程學(xué)院徐正玉博士對文中算法設(shè)計提供了幫助！