郭嵩巍,劉小畔,鄭凱,張磊

(1.中國地質(zhì)大學(xué)(北京) 地球科學(xué)與資源學(xué)院，北京 100083； 2.浙江省地礦勘探院，浙江 杭州 310000； 3.長江大學(xué) 地球物理與石油資源學(xué)院，湖北 武漢 430100； 4.內(nèi)蒙古國土資源勘查開發(fā)有限責(zé)任公司,內(nèi)蒙古 呼和浩特 010020)

0 引言

瞬變電磁法作為一種常見的電磁測深物探方法，廣泛應(yīng)用于水文勘查、工程勘查以及礦產(chǎn)勘查[1-2]。在目前的工程應(yīng)用中，除“煙圈”反演快速成像法外，一維反演結(jié)果仍然是地質(zhì)解釋的主要依據(jù)。對于中心回線裝置，可以解析推導(dǎo)出一維正演響應(yīng)的計算公式[3]，其數(shù)值計算的核心是，內(nèi)層積分的漢克爾變換[4]、外層積分的余弦變換[5]均可采用數(shù)值濾波算法，隨著濾波系數(shù)的不斷推陳出新，計算量明顯減少，大大提高了正演的計算效率和精度。在反演方面，目前常見的反演方法，如阻尼最小二乘法(馬奎特法)[6]、Occam法[7-10]、正則化法[11-13]等，都需要計算正演響應(yīng)對模型的導(dǎo)數(shù)，即雅克比矩陣，因此雅克比矩陣的計算效率和效果就成為反演計算的關(guān)鍵因素之一。計算雅克比矩陣最簡便的方法是差分法，該算法程序?qū)崿F(xiàn)較為簡單，缺點是不僅調(diào)用正演次數(shù)多，計算效率低，而且有可能增加累積誤差。因此，本文詳細(xì)推導(dǎo)了中心回線裝置瞬變電磁一維正演全區(qū)視電阻率的雅克比矩陣的解析計算公式，并通過Occam反演，對比分析了理論模型和實測數(shù)據(jù)的反演效果，證明推導(dǎo)正確，值得采用。

1 瞬變電磁正演響應(yīng)之全區(qū)視電阻率

目前，瞬變電磁儀器大多采集的數(shù)據(jù)為感應(yīng)電動勢，對其歸一化后可轉(zhuǎn)換為晚期視電阻率，然而晚期視電阻率是假設(shè)在時間趨于無窮大時近似推導(dǎo)所得，因此在響應(yīng)早期，視電阻率曲線變形嚴(yán)重，呈明顯高阻假象，成為了地質(zhì)解釋重要的干擾因素之一，因此很多學(xué)者提出了全區(qū)視電阻率(也稱全期視電阻率)的概念[14]。由于通過磁場強度計算的全區(qū)視電阻率更為接近真正的視電阻率定義[15]，可通過感應(yīng)電動勢轉(zhuǎn)換為磁場強度計算得出。因此，本文將該全區(qū)視電阻率作為Occam反演的擬合參數(shù)。

對于中心回線瞬變電磁的一維正演，可直接求解磁場強度。假設(shè)n層地電模型，第j層電阻率為ρj、層厚度為hj，頻率域的磁場強度響應(yīng)為：

(1)

其中：a為發(fā)射線框等效半徑；u的遞推公式[3]為

(2)

対于發(fā)射波形的選擇，大多瞬變電磁儀普遍采用占空比1∶1的方波，如V8電法工作站的TD50方波，在正演模擬中只需考慮該波形的半個周期，可視為階躍波：

(3)

式中I0為發(fā)射電流。對上式進行傅里葉變換，得到頻率域中發(fā)射電流表達(dá)式：

(4)

于是，瞬變電磁響應(yīng)可化簡為：

(5)

(6)

對于均勻半空間地電模型，可以推導(dǎo)出Hz的解析表達(dá)式[3]：

(7)

其中erf(x)為誤差函數(shù)。歸一化Hz得到函數(shù)：

(8)

最后求取Z(x)的反函數(shù)x，得出基于磁場強度的全區(qū)視電阻率：

(9)

由于函數(shù)Z(x)為隱函數(shù)，無法推導(dǎo)出其反函數(shù)，因此無法通過解析算法求解，但可利用函數(shù)在0～1之間單調(diào)遞增的特點，通過二分查找的辦法得到數(shù)值解。筆者曾采用該方法實現(xiàn)了全區(qū)視電阻率的計算以及對應(yīng)的“煙圈”反演，取得了不錯的效果[15]。

可以看到，瞬變電磁正演響應(yīng)計算可簡化為：首先計算核函數(shù)，然后再對內(nèi)層積分(計算頻率域的電磁響應(yīng))做漢克爾變換，對外層積分(頻率域到時間域的轉(zhuǎn)換)做余弦變換(式(6))，得到Hz和Z(x)，最后采用二分查找法通過Z(x)求出x，從而得到正演響應(yīng)的全區(qū)視電阻率。

2 雅克比矩陣的解析推導(dǎo)

式(9)對模型參數(shù)求導(dǎo)，根據(jù)復(fù)合函數(shù)的鏈?zhǔn)角髮?dǎo)法則，層狀介質(zhì)模型任意層j的模型電阻率ρj及厚度hj的導(dǎo)數(shù)為：

(10)

(11)

利用反函數(shù)導(dǎo)數(shù)與原函數(shù)導(dǎo)數(shù)的倒數(shù)關(guān)系，得到：

(12)

根據(jù)Z(x)與Hz的關(guān)系式(8)，得到

(13)

將式(11)～式(13)帶入式(10)，化簡得到：

(14)

其中：

(15)

對核函數(shù)導(dǎo)數(shù)的計算，根據(jù)鏈?zhǔn)角髮?dǎo)法則得：

(16)

(17)

令：

(18)

(19)

其中，

(20)

于是模型導(dǎo)數(shù)式(14)可表示為：

(21)

(22)

(23)

(24)

根據(jù)正演公式(式(1)、式(2))，進而推出第n層的模型導(dǎo)數(shù)關(guān)系：

(25)

(26)

3 Occam反演

Constable等提出Occam反演法[7]，其基本原理是尋找在目標(biāo)函數(shù)相對極小情況下的光滑模型，它具有迭代過程穩(wěn)定、收斂較快且對初始模型依賴度較小等優(yōu)點。該方法在反演過程中每迭代一次需要對模型求一次導(dǎo)數(shù)，因此由模型導(dǎo)數(shù)構(gòu)成的雅克比矩陣的計算精度和計算效率將決定Occam反演的計算效果和計算效率。Occam反演原理簡述如下。

m為正演模型向量，d為正演響應(yīng)向量，表示為：

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

(27)

數(shù)據(jù)擬合差：

χ2=‖Wd-WF[m]‖2。

(28)

根據(jù)約束極小化理論，利用拉格朗日乘子形成一個無約束的函數(shù)U:

(29)

式中?為粗糙度矩陣，可表示為：

(30)

式中右端第一項是粗糙度，第二項是由拉格朗日乘子加權(quán)的擬合差；取梯度，引起U不變的向量m遵循

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

(31)

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

(32)

可表示為：

(33)

上式可通過式(21)、(22)、(26)計算得出。

最終，假設(shè)第k次反演迭代已完成，k+1次迭代的模型向量為：

mk+1(μ)=[μ?T?+(WJk)TWJk]-1(WJk)TWdk，

(34)

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

χk+1(μ)=‖Wd-WF[mk+1(μ)]‖，

(35)

其中：μ為拉格朗日乘子，χ為擬合差。

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

4 理論模型算例

為驗證解析算法的正確性，本文采用與中間差分算法

(36)

做對比。從理論上說，差分步長越短，結(jié)果越接近真實值，但步長過短，會由于正演的精度不夠而偏差增大，一般采用的經(jīng)驗差分步長(Δm/m)取0.001。

對典型三層地電模型做一維Occam反演：不失一般性，Occam反演的初始模型選擇均勻半空間，電阻率取響應(yīng)的平均值，模型層采用對數(shù)等間距遞增，起始厚度5 m，對數(shù)遞增間距為0.1，共分10層，最大反演深度687.2 m。

H型地電模型：層電阻率分別為100、10、100 Ω·m，層厚度分別為60、30 m，計算結(jié)果見圖1、表1。K型地電模型：層電阻率分別為10、100、10 Ω·m，層厚度分別為30、60 m，計算結(jié)果見圖2、表2。

從上述算例可以看出，差分算法不僅計算效率低，而且曲線擬合不佳，反演模型與真實模型差距大，而解析算法則表現(xiàn)優(yōu)秀。差分算法效果不佳的原因很可能與全區(qū)視電阻率的計算有關(guān)，因為全區(qū)視電阻率在計算隱函數(shù)Z(x)的反函數(shù)時，是通過二分查找后插值所得，本身誤差較大，再對其進行差分求導(dǎo)，放大了誤差。因為若差分步長較大，差分本身的誤差較大；而步長較小，正演響應(yīng)Hz變化不大，二分查找的結(jié)果非常接近而趨于相等，差分結(jié)果趨近于零，總之會放大誤差。隨著反演不斷迭代，誤差會進一步累積放大，最終導(dǎo)致曲線擬合困難、反演模型震蕩的結(jié)果(圖1、圖2)。因此，對于全區(qū)視電阻率的反演不建議使用差分算法。

表1 H型模型反演結(jié)果Table 1 Inversion result parameters of H-type layered model

圖1 H型反演結(jié)果對比Fig.1 Comparison of H-type layered Model inversion results

圖2 K型反演結(jié)果對比Fig.2 Comparison of K-type layered Model inversion results

表2 K型模型反演結(jié)果Table 2 Inversion result parameters of K-type layered model

5 實測數(shù)據(jù)實例

解析算法對理論模型的反演表現(xiàn)良好。為進一步驗證解析算法Occam反演對實測數(shù)據(jù)的有效性，采用對實測剖面數(shù)據(jù)進行逐點一維反演的方法，即擬二維反演，繪制成電阻率反演斷面圖。與全區(qū)視電阻率擬斷面、“煙圈”反演的結(jié)果進行效果對比，并通過地質(zhì)鉆探來驗證方法的有效性。

5.1 地質(zhì)概況

原始數(shù)據(jù)采集自內(nèi)蒙古自治區(qū)涼城縣西水塘鎮(zhèn)附近，如圖3所示，為內(nèi)蒙古涼城縣某地?zé)峥辈轫椖恐械囊粭l剖面(400線)，測線位于山前傾斜平原，布置該測線的目的是為劃分工作區(qū)地層，推斷基底埋深，并進一步驗證地質(zhì)推斷的構(gòu)造，為地?zé)嵴宜峁┪锾揭罁?jù)。工作區(qū)地形平坦，空曠且開闊，無工業(yè)電流干擾，非常有利于開展地面瞬變電磁工作。

根據(jù)區(qū)域地質(zhì)資料顯示[16]，工作區(qū)地層主要有太古宇界桑干群，中生界侏羅系及新生界古近系、新近系、第四系(圖3)。區(qū)域地質(zhì)構(gòu)造屬陰山東西向構(gòu)造帶的南緣，主要有NE和NW兩個方向，其中NE向構(gòu)造屬高角度壓性斷層，斷層連續(xù)性強，是地下熱水和溫泉形成的主要通道。侵入巖主要分布在蠻漢山區(qū)及平原區(qū)深部，主要為太古宙代早期侵入巖，巖性以蘇長巖和似斑狀花崗巖為主。

1—中低山; 2—山前傾斜平原; 3—第四系全新統(tǒng)沖洪積物; 4—侏羅系火山巖; 5—太古宇桑干群片麻巖; 6—太古宙似斑狀花崗巖; 7—斷裂破碎帶; 8—測試瞬變電磁測深測線; 9—地質(zhì)界線; 10—地貌界線; 11—已有鉆孔；12—驗證地?zé)徙@孔(DR2); 13—瞬變電磁發(fā)射線框位置; 14—瞬變電磁接收線圈位置1—middle low mountain; 2—piedmont inclined plain; 3—Quaternary Holocene alluvial proluvial; 4—Jurassic volcanic rock; 5—Archaean Sanggan group gneiss; 6—Archaean porphyry like granite; 7—fracture fracture Fracture Zone; 8—test TEM sounding line; 9—geological boundary; 10—geomorphic boundary; 11—existing drilling；12—verification of geothermal drilling (DR2); 13—TEM emission line frame position; 14—TEM receiving coil position圖3 工作區(qū)地質(zhì)圖Fig.3 Geological map of work area

測線所在工作區(qū)全部為第四系覆蓋。據(jù)區(qū)內(nèi)已有鉆孔顯示(圖3，圖4a)，巖性從下至上分別為：太古宙似斑狀花崗巖(γ)，多為灰色中粗粒似斑狀花崗巖；古近系(E)為深棕紅色泥巖、泥質(zhì)砂巖、砂質(zhì)泥巖夾砂礫巖；新近系(N)上部為淺棕紅色泥巖、泥質(zhì)砂巖、砂質(zhì)泥巖夾砂礫巖，下部為土黃色、灰白色砂礫巖，與古近系地層呈假整合接觸；第四系(Q)巖性為淺黃棕色、黃褐色粉土、粉質(zhì)黏土與砂礫卵石、含卵礫中粗砂、中粗砂等。

圖4 工作區(qū)鉆孔柱狀圖Fig.4 Borehole histogram of work area

5.2 Occam反演及地質(zhì)解釋

工作裝置的設(shè)置：根據(jù)工作區(qū)內(nèi)已知鉆孔顯示(圖4a)，工作區(qū)地層分層明顯，侵入巖、沉積巖具有明顯的縱向電性差異，具備了瞬變電磁測深施工的物性條件。測線所在位置基底埋深推斷為在200～300 m之間，且地層中包括較厚的一套泥巖地層，方法試驗也證實了這套低阻泥巖的存在。 為穿透這套低阻地層探測到基底花崗巖，在充分利用場地寬度的情況下，設(shè)置發(fā)射裝置中心回線為400 m×600 m。

本著從已知到未知的推斷原則，設(shè)定Occam反演初始模型。初始模型采用深度對數(shù)等間隔模型，讓反演模型與瞬變電磁響應(yīng)的采樣間隔方式一致，這樣一方面有利于曲線擬合，另一方面可以盡可能多分層，提高反演的縱向分辨能力。經(jīng)過多次實驗，以兼顧反演計算效率和突出局部異常為原則，確定初始模型起始深度為20 m，避開瞬變電磁響應(yīng)早期的盲區(qū)，對數(shù)間隔1.15，分25層，最小層厚3 m，最大層厚65 m，最大反演深度498 m；初始模型電阻率參考全區(qū)視電阻率的均值，定為60 Ω·m。

瞬變電磁測深視電阻率曲線大多表現(xiàn)為KH、KQH型，反演結(jié)果總體可分為5個大電性層(圖5)。第一層為中低阻，電阻率在 50 Ω·m左右，反演深度在0～30 m左右，與鉆孔揭露的第四系粉土、黏土對應(yīng)。第二層為中高阻層，電阻率在50～80 Ω·m之間，反演深度在30～80 m之間，對應(yīng)以砂巖、砂礫巖、粉砂巖等為主的第四系地層，透水性好，為含水層。第三層為低阻層，電阻率緩慢下降，在10～50 Ω·m之間，反演深度80～150 m左右，與第四系地層粉質(zhì)黏土對應(yīng)，隔水性較好，為隔水層。第四層為低層，電阻率在10 Ω·m左右，反演深度150～300 m之間，與新近系，古近系的紅色泥巖對應(yīng)，為隔水層。第五層為高阻，電阻率大于100 Ω·m，對應(yīng)基底太古代似斑狀花崗巖。

圖5 鉆孔投影點820點反演結(jié)果Fig.5 Inversion result of site 820 where is projectionposition of borehole

上述地質(zhì)推斷，經(jīng)鉆孔DR2驗證(圖4b)，投影到剖面位置為820點，與實際情況基本吻合。該井定位于構(gòu)造中，鉆探深度 308.6 m，成井后自流，自流量為4 483.64 m3/d，水溫 36 ℃，礦化度 1.01 g/L，是內(nèi)蒙古涼城縣地?zé)峥辈橹腥〉玫耐怀龀晒鸞17]。

5.3 反演結(jié)果對比分析

利用本實測數(shù)據(jù)地層分層明確、縱向電性差異明顯的特點，可測試反演算法在縱向上的分辨能力，因此繪制了全區(qū)視電阻率擬斷面和“煙圈”反演斷面圖，與Occam反演結(jié)果對比分析(圖6)：擬斷面和“煙圈”反演總體表現(xiàn)為“高—低—高”三層模型，中間層表現(xiàn)為巨厚的低阻層，對地層的分辨能力有限，而Occam反演則反演出第四系次高阻層，對應(yīng)含砂巖、粉砂巖較多的地層，該類地層滲透性好，一般為含水層，對找水工作意義重大；對于高阻基底，上述兩種的方法的反演結(jié)果基本一致，與鉆探結(jié)果基本吻合。

a—擬斷面;b—“煙圈”反演結(jié)果；c—Occam反演結(jié)果a—pseudo section; b—smoke ring inversion section; c—Occam inversion section圖6 400線瞬變電磁反演電阻率斷面對比Fig.6 Comparison of TEM inversion sections of line 400

6 結(jié)論

采用基于瞬變電磁全區(qū)視電阻率的Occam反演，有效避免了晚期視電阻率早期變形嚴(yán)重而造成的假高阻。在理論正演模型算例中顯示：雅克比矩陣的解析算法在計算效果和計算效率兩方面都優(yōu)于差分算法；在實測數(shù)據(jù)算例中可以看到：Occam反演結(jié)果與鉆孔更為吻合，對地層的縱向分辨能力高于煙圈反演，特別是在第四系低阻地層中次高阻砂巖層的反映，對瞬變電磁找水工作具有較高的實用價值。上述算例證明，解析算法公式推導(dǎo)正確，基于全區(qū)視電阻率的Occam反演不僅反演效率高，而且效果好。

此外，本文不僅推導(dǎo)了Occam反演中雅克比矩陣所需的對模型層電阻率的求導(dǎo)公式，還推導(dǎo)出了對模型層厚度的求導(dǎo)公式，可應(yīng)用于對模型層厚度反演的方法，如馬奎特反演等。

致謝：感謝內(nèi)蒙古國土資源勘查開發(fā)院以及內(nèi)蒙古自治區(qū)第一水文地質(zhì)工程地質(zhì)勘查院提供的寶貴實測數(shù)據(jù)及相關(guān)地質(zhì)資料;感謝審稿老師提出的寶貴意見。