高阻體模型瞬變電磁響應(yīng)特征及其反演淺析

陳 輝，杜興忠，孫永清

(中國電建集團(tuán) 貴陽勘測(cè)設(shè)計(jì)研究院有限公司，貴陽 550081)

0 引言

瞬變電磁法(Time domain electromagnetic methods，TEM)是利用線圈或接地電極觀測(cè)二次場(chǎng)(響應(yīng)場(chǎng))，通過對(duì)這些響應(yīng)信息提取和分析，從而達(dá)到探測(cè)地下地質(zhì)體的方法[1]。國內(nèi)、外許多學(xué)者對(duì)于瞬變電磁法從1維正反演，到2.5維正反演，再到3維正反演做了大量工作，均取得了很多有價(jià)值的研究成果，但這些研究大部分是基于低阻異常體的響應(yīng)特征出發(fā)，很少研究純高阻體地質(zhì)狀況的響應(yīng)特性。筆者從建立1維高阻層狀模型和2.5維高阻矩形柱模型出發(fā)，分別研究其正演響應(yīng)特征，并依據(jù)正演響應(yīng)結(jié)果，對(duì)反演做了初步分析。

1 TEM 2.5維正反演概述

1.1 正演模擬

2.5維正演模擬先利用拉氏變換，將2.5維電磁問題從時(shí)間域轉(zhuǎn)換到拉氏域，再利用傅氏變換從拉氏域轉(zhuǎn)換到傅氏域中，并構(gòu)建拉氏傅氏域中電磁場(chǎng)分量沿走向方向的對(duì)偶二階微分方程，由對(duì)偶二階微分方程和邊界條件形成邊值問題；然后，采用有限元法求解該邊值問題，得到電磁場(chǎng)的解。將求解的方程結(jié)果經(jīng)過逆拉氏傅氏變換，可以得到各節(jié)點(diǎn)在時(shí)間域三維空間的電磁場(chǎng)分量。

對(duì)于瞬變電磁2.5維的對(duì)偶二階方程的推導(dǎo)，前人已做了大量研究，這里直接引用，中心回線TEM 2.5維沿走向方向的電磁場(chǎng)二階對(duì)偶偏微分方程形式為[2-3]：

a4u+b4up

(1)

a5v

(2)

其中，

(3)

(4)

其中，系數(shù)ηu和ηv計(jì)算式為：

(5)

(6)

電磁場(chǎng)二階對(duì)偶微分式(1)和式(2)，邊界條件式(3)和式(4)四個(gè)方程一起構(gòu)成了在(s,m)域中沿走向方向(這里為y方向)電磁分量的邊值問題，采用有限元法進(jìn)行求解。

(7)

式中：求和下限t=i1是(m+1)/2的整數(shù)部分。

采用余弦變換數(shù)值濾波算法做逆傅氏變換，根據(jù)正余弦變換數(shù)值濾波算法[5]，得到式(8)。

(8)

1.2 反演分析

1.2.1 阻尼最小二乘法

瞬變電磁常采用的是阻尼最小二乘法進(jìn)行反演解釋，它是1963年馬奎特(Marquardt)[6-7]通過改進(jìn)高斯-牛頓法的法方程組而得到的一種新的反演方法。對(duì)于目標(biāo)函數(shù)式(9)的極小值問題，阻尼最小二乘法的處理方式為式(10)。

(9)

(JTJ+λI)δ=ε

(10)

其中：J表示雅可比矩陣；λ表示阻尼因子(它是一個(gè)正數(shù)，用以控制校正的方向和步長)；δ是修正量即δ=p(k+1)-p(k)，從式(10)中可以明顯看出它是λ的函數(shù)即δ=δ(λ)，ε是殘差向量。所以阻尼最小二乘法的迭代公式為：

p(k+1)=p(k)+(JTJ+λI)-1ε

(11)

在反演算法中，雅克比矩陣J的計(jì)算精度與反演結(jié)果的準(zhǔn)確性有直接關(guān)系，設(shè)從模型空間映射到數(shù)據(jù)空間的函數(shù)為φ，模型參數(shù)用p表示。對(duì)雅克比矩陣行列中的每個(gè)元素的一階偏導(dǎo)數(shù)的計(jì)算，采用有限差分方法中的一階向前差分為式(12)。

(12)

阻尼最小二乘反演法的反演是一個(gè)循環(huán)迭代的過程，迭代一次修改一次模型參數(shù)，直到由模型計(jì)算出的電磁響應(yīng)值與實(shí)測(cè)或者理論數(shù)據(jù)間的差值滿足設(shè)定的要求。

1.2.2 阻尼特征參數(shù)法

阻尼特征參數(shù)反演法(Damped Eigenparameter Inversion approach)，是由Jupp[8]和Vozoff[9]提出并用于地球物理數(shù)據(jù)聯(lián)合反演解釋中，它結(jié)合了馬奎特方法和奇異值截?cái)嗟膬?yōu)點(diǎn)，將經(jīng)奇異值分解后得到的特征值分為重要參數(shù)和不重要參數(shù)，截?cái)嘈〉钠娈愔?，通過改變奇異值來修正阻尼因子。相對(duì)于阻尼最小二乘方法從兩次計(jì)算結(jié)果來改變阻尼因子，該方法綜合所有的數(shù)據(jù)特性，從整體上指導(dǎo)反演過程的進(jìn)行，具有快速穩(wěn)定的特點(diǎn)。

對(duì)于時(shí)間域瞬變電磁勘探，假設(shè)M個(gè)測(cè)點(diǎn)感應(yīng)電動(dòng)勢(shì)的數(shù)據(jù)表示成向量為d=(d1,d2,…,dM)T，表示地下介質(zhì)電導(dǎo)率等電性參數(shù)的模型參數(shù)可以表示成P=(P1,P2,…,PN)T。

大多數(shù)地球物理問題并不是線性的，對(duì)于非線性反演問題的線性化用雅可比矩陣進(jìn)行描述，它是一種從N維參數(shù)空間到M維數(shù)據(jù)空間的線性映射。引入誤差向量ε，方程可寫成：

(13)

當(dāng)采用式(13)計(jì)算時(shí)，如果奇異值Sj較小，計(jì)算會(huì)出現(xiàn)不穩(wěn)定甚至不收斂的情況，為解決該問題，可以通過阻尼因子tj改善反演過程，于是引入一個(gè)阻尼修正步長：

(14)

每迭代一次就修改一次模型參數(shù)，由式(12)可得出，新模型參數(shù)的計(jì)算式為式(15)。

P(new)=P(old)+VδQ

(15)

新模型參數(shù)生成后，重新計(jì)算新模型的響應(yīng)值和誤差向量。反演過程就是一個(gè)循環(huán)迭代的過程，直到滿足預(yù)先設(shè)置的結(jié)束準(zhǔn)則迭代終止。

1.2.3 反演計(jì)算步驟

本文實(shí)現(xiàn)反演的計(jì)算步驟如下：

1)設(shè)置反演參數(shù)，讀取反演觀測(cè)數(shù)據(jù)。

2)設(shè)置初始模型，計(jì)算其正演理論響應(yīng)值。

3)計(jì)算模型理論響應(yīng)值和反演觀測(cè)數(shù)據(jù)的誤差。

4)判斷誤差是否滿足要求，不滿足要求跳到步驟5)，滿足要求直接跳到步驟9)。

5)計(jì)算雅克比矩陣。

6)計(jì)算新的模型參數(shù)及其理論響應(yīng)值。

7)計(jì)算新模型的理論響應(yīng)值和反演觀測(cè)數(shù)據(jù)的誤差。

8)再次判斷誤差是否滿足要求或者達(dá)到最大迭代次數(shù)，不滿足要求或未達(dá)到最大迭代次數(shù)則跳到步驟5)，滿足要求或已達(dá)到最大迭代次數(shù)則跳到步驟9)。

9)輸出得到迭代結(jié)束后的模型參數(shù)。

2 正演模擬

2.1 高阻層狀模型1維TEM響應(yīng)特征

根據(jù)TEM一維正演理論[10-12]，構(gòu)建一維層狀K型地電模型如圖1，具體模型參數(shù)見表1。

圖1 層狀模型

表1 層狀模型參數(shù)

高阻層狀模型1～高阻層狀模型4和兩種均勻半空間模型(電阻率為10 Ω·m和103Ω·m)的瞬變電磁響應(yīng)曲線共同繪入圖2，從圖2中可以看到，對(duì)于相同電阻率的高阻層：同一層厚時(shí)(模型1和模型2、模型3和模型4)，圍巖電阻率相對(duì)越高，高阻層響應(yīng)曲線進(jìn)入衰減的時(shí)間相對(duì)越早。對(duì)于高阻層狀模型整體來看，在關(guān)斷電流后的某個(gè)時(shí)刻(設(shè)這個(gè)時(shí)刻為t1)后響應(yīng)曲線衰減速度相對(duì)較快，然后在某一時(shí)刻(設(shè)這個(gè)時(shí)刻為t2)達(dá)到一個(gè)最大值，其后衰減速度逐漸減慢。t1、t2值可通過模型1～模型4與對(duì)應(yīng)均勻半空間模型的相對(duì)比值來確定(圖3)，比值越遠(yuǎn)離1線表示模型響應(yīng)值和均勻半空間響應(yīng)值差異相對(duì)越大，反之亦然。圖3中比值遠(yuǎn)離1線的起始位置(可選擇與1線的相對(duì)誤差等于0.01%作為起始位置判定點(diǎn))對(duì)應(yīng)的時(shí)刻即為t1的值，離1線垂直距離最遠(yuǎn)的點(diǎn)對(duì)應(yīng)的時(shí)刻即為t2值。從圖2、圖3中還可以看到，同一圍巖電阻率時(shí)(模型1和模型3、模型2和模型4)，高阻層相對(duì)越厚，進(jìn)入衰減過程t1時(shí)刻后正演響應(yīng)隨時(shí)間衰減相對(duì)越快，達(dá)到衰減最大t2時(shí)刻相對(duì)較晚，反之亦然。

圖2 高阻層狀模型TEM響應(yīng)曲線

圖3 模型和均勻半空間TEM響應(yīng)的相對(duì)比值

層狀模型5～層狀模型8的高阻層電阻率差異不同倍數(shù)時(shí)TEM響應(yīng)見圖4、圖6、圖8、圖10。采用相鄰倍數(shù)響應(yīng)的比值曲線對(duì)比高倍數(shù)與前一倍數(shù)響應(yīng)的差異情況，見圖5、圖7、圖9、圖11，并計(jì)算各模型相鄰倍數(shù)響應(yīng)相對(duì)誤差均值見表2。從圖3～圖10和表2中可以看到，地層層厚相同時(shí)，同一模型在高阻層電阻率超過圍巖電阻率103倍后，相鄰倍數(shù)TEM響應(yīng)的相對(duì)誤差均值均小于0.03%，說明超過103倍后，即使高阻層電阻率較圍巖電阻率變化很大，但TEM響應(yīng)值變化卻很小(幾乎無變化)。

表2 相鄰倍數(shù)TEM響應(yīng)相對(duì)誤差平均值表

圖4 電阻率差異不同倍數(shù)時(shí)TEM響應(yīng)曲線(模型5)

圖5 相鄰倍數(shù)響應(yīng)的比值曲線(模型5)

圖6 電阻率差異不同倍數(shù)時(shí)TEM響應(yīng)曲線(模型6)

圖7 相鄰倍數(shù)響應(yīng)的比值曲線(模型6)

圖8 電阻率差異不同倍數(shù)時(shí)TEM響應(yīng)曲線(模型7)

圖9 相鄰倍數(shù)響應(yīng)的比值曲線(模型7)

圖11 相鄰倍數(shù)響應(yīng)的比值曲線(模型8)

2.2 高阻矩形柱模型2.5維TEM響應(yīng)特征

本次研究采用有限元法，將所要求解的連續(xù)場(chǎng)離散為很多相互連接的較小單元，把連續(xù)的求解域離散化，在每個(gè)單元內(nèi)用插值函數(shù)進(jìn)行插值，建立起每個(gè)單元網(wǎng)格各自的公式，最后再將這些單元集合起來求解，以各單元節(jié)點(diǎn)函數(shù)值為未知量的線性方程組，得出原來的場(chǎng)值問題的解[13]。建立離散網(wǎng)格，向下前5層每層層厚設(shè)置為10 m，其后的每層厚度比上層按6%遞增，深度范圍0 m～400 m；橫向采用等間距為20 m的均勻網(wǎng)格，寬度為-650 m～650 m；測(cè)點(diǎn)范圍為-590 m～590 m，測(cè)點(diǎn)間距為20 m。發(fā)射電流最大值為10 A，中心回線發(fā)射線圈半徑為13 m，匝數(shù)為4匝；發(fā)射電流波形為方波，電流平穩(wěn)時(shí)間為7.6 ms，7.6 ms關(guān)斷電流；二次場(chǎng)采樣時(shí)間為9.60 ms，26個(gè)時(shí)間采樣道。均勻半空間中建立矩形柱體模型，模型頂部離地面深度為50 m，寬度為120 m，模型深度為115 m(圖12)。設(shè)圍巖的電阻率為103Ω·m，高阻體的電阻率為105Ω·m，低阻體的電阻率為10 Ω·m。

圖12 地下單一矩形異常體模型(模型9)

矩形柱模型的2.5維TEM響應(yīng)曲線見圖13，從圖13中可以看到：均勻半空間中高阻、低阻矩形柱體的2.5維TEM響應(yīng)曲線，在異常體的正上方位置(測(cè)點(diǎn)0值位置)均存在明顯的異常形態(tài)，兩者的異常響應(yīng)形態(tài)相反[14]。

圖13 矩形柱體2.5維TEM響應(yīng)曲線

參考模型9各個(gè)參數(shù)，構(gòu)建地下高阻、低阻同時(shí)存在的2.5維柱狀體模型(圖14)，測(cè)點(diǎn)范圍為-590 m～590 m，測(cè)點(diǎn)間距為20 m，圍巖的電阻率為103Ω·m，高阻體的電阻率為105Ω·m，低阻體的電阻率為10 Ω·m。得到其瞬變電磁響應(yīng)曲線如圖15所示。從圖14中可以看到，當(dāng)同時(shí)存在低阻和高阻模型時(shí)，低阻體TEM響應(yīng)曲線非常明顯，而高阻體響應(yīng)曲線盡管也有一定異常形態(tài)，但相對(duì)低阻體響應(yīng)曲線而言并不明顯。

圖14 地下兩個(gè)矩形異常體模型(模型10)

圖15 地下兩個(gè)矩形異常體TEM響應(yīng)曲線

3 反演淺析

根據(jù)高阻體模型的2.5維正演響應(yīng)特征，采用阻尼最小二乘反演方法和阻尼特征參數(shù)反演法，分別對(duì)模型10進(jìn)行反演分析，兩種反演方法結(jié)果數(shù)據(jù)對(duì)比見表3，反演結(jié)果圖見圖16。從表3中可以看到，阻尼特征參數(shù)反演法均方根相對(duì)誤差(RMS)和噪信比(RSR)較低，迭代次數(shù)較少，收斂較快，反演耗時(shí)較少；從圖16中可以看出，兩種反演方法的效果均能較準(zhǔn)確地反演出低阻體的特征，也能在一定程度上反演出高阻體的形態(tài)，但阻尼最小二乘法對(duì)高阻體的反演結(jié)果信息并沒有阻尼特征參數(shù)法豐富，說明阻尼特征參數(shù)反演法對(duì)高阻體的反演信息量及效果相對(duì)較好。

圖16 高阻-低阻矩形柱體模型反演電導(dǎo)率結(jié)果

表3 兩種反演方法結(jié)果數(shù)據(jù)對(duì)比

4 結(jié)果與認(rèn)識(shí)

通過建立1維高阻層狀模型和2.5維高阻矩形柱模型，模擬了地下高阻體地電模型的瞬變電磁正演響應(yīng)，討論了圍巖電阻率差異、高阻體電阻率差異等情況下的1維響應(yīng)特征，分析了高阻矩形柱體的2.5維響應(yīng)特征，并對(duì)其響應(yīng)結(jié)果進(jìn)行反演分析，得到了一些認(rèn)識(shí)如下：

1)對(duì)1維高阻層狀模型正演模擬可知：同一層厚時(shí)圍巖電阻率相對(duì)越高，高阻層響應(yīng)曲線進(jìn)入衰減的時(shí)間相對(duì)越早；高阻層越厚，進(jìn)入衰減過程t1時(shí)刻后TEM響應(yīng)隨時(shí)間衰減相對(duì)越快，達(dá)到衰減最大t2時(shí)刻較晚，反之亦然。地層層厚相同時(shí)，同一模型在高阻層電阻率超過圍巖電阻率103倍后，其TEM響應(yīng)值變化很小(幾乎無變化)。

2)對(duì)2.5維高阻矩形柱狀體模型正演模擬可知：當(dāng)?shù)叵聠为?dú)存在高阻或低阻異常體時(shí)，在異常體位置TEM響應(yīng)明顯，兩者的異常響應(yīng)形態(tài)相反。當(dāng)?shù)叵峦瑫r(shí)存在低阻和高阻異常體時(shí)，低阻體TEM響應(yīng)明顯，而高阻體響應(yīng)并不明顯。

3)分別采用阻尼最小二乘反演方法和阻尼特征參數(shù)反演法,對(duì)2.5維高阻矩形柱狀體模型進(jìn)行反演分析可知：阻尼特征參數(shù)反演法均方根相對(duì)誤差和噪信比較低，迭代次數(shù)較少，收斂較快，反演耗時(shí)較少，對(duì)高阻體的反演結(jié)果信息相對(duì)豐富、效果較明顯。

4)由于影響瞬變電磁響應(yīng)的因素復(fù)雜，筆者在各向同性介質(zhì)中所建模型存在局限性，在實(shí)際情況下，即使相同條件，也很難達(dá)到相同效果，這也是今后需要進(jìn)一步研究的方向。