瞬變電磁全期視電阻率的弦截無限逼近算法

曾慶寧 羅 瀛 劉 帥 龍 超

(①桂林電子科技大學認知無線電與信息處理教育部重點實驗室，廣西桂林 541004；②桂林電子科技大學信息與通信學院，廣西桂林 541004)

0 引言

瞬變電磁法(TEM)在石油勘探 、礦產(chǎn)調(diào)查、地下水勘察、工程地質(zhì)勘探等諸多領(lǐng)域得到了廣泛應(yīng)用[1-6]。傳統(tǒng)的瞬變電磁反演是將測量的感應(yīng)電動勢轉(zhuǎn)化為地質(zhì)體的視電阻率[7-10]，從而得到探區(qū)真實的地質(zhì)圖像。

可通過近似公式法[1,3](包括后期近似公式與前期近似公式)或者改進的近似公式法[11-13]，由感應(yīng)電動勢求解視電阻率，這些方法計算相對簡單，但通常存在較大誤差，而且在某些時段，尤其是在反映中淺地層地質(zhì)的早期時段，這種誤差可能導致錯誤的地質(zhì)反演結(jié)果。特別是當目標體積較小或目標體的電阻率與圍巖的電阻率相近時，這種誤差更可能導致錯誤反演結(jié)果。為此，可使用后全期法和前全期法進行無限精確求解[14-17]，以克服這種誤差的影響，使反演成像更為精細和準確，但這種無限精確求解法需要使用復(fù)雜非線性函數(shù)的求根技術(shù)。

二分法和牛頓法是兩種行之有效的非線性函數(shù)無限精確求根方法[18-21]，前者已成功應(yīng)用于地面瞬變電磁數(shù)據(jù)的視電阻率任意精度求解[14-17]，后者成功應(yīng)用于井下瞬變電磁數(shù)據(jù)的視電阻率任意精度求解[18]。然而，二分法通常需要較多的迭代次數(shù)，收斂速度較慢；牛頓法需要首先求出瞬變電磁感應(yīng)電動勢函數(shù)的導數(shù)，對于井下瞬變電磁信號，該導數(shù)的表達式并不復(fù)雜，但對于地面瞬變電磁信號，其表達式則非常復(fù)雜，尤其是其計算誤差很大，往往導致算法失效。

本文引進非線性函數(shù)求根的弦截方法[21]，并將其用于地面重疊回線或中心回線瞬變電磁視電阻率的任意精度求解，給出了弦截無限逼近算法。文中還進一步說明了在同等精度要求下，弦截法比二分法所需的迭代次數(shù)更少，具有更快的收斂速度。同時，弦截法可避免類似牛頓法那樣需要對瞬變電磁信號求導的問題，為瞬變電磁數(shù)據(jù)的視電阻率反演提供了另一種滿足任意精度要求且收斂速度較快的計算方法。最后針對地下水倉的實測瞬變電磁數(shù)據(jù)，應(yīng)用弦截法求解其全期視電阻率的無限精確解，并據(jù)此進行了地質(zhì)反演成像，效果較好。

1 全期視電阻率

瞬變電磁法中，經(jīng)常采用重疊回線或中心回線法采集數(shù)據(jù)，接收線圈于時刻t測得的感應(yīng)電動勢為[1]

(1)

式中核函數(shù)

(2)

(3)

通過式(1)求解視電阻率ρ的方法稱為全期法?？梢宰C明，F(xiàn)(z)是z的上凸函數(shù)[20-21]，因此可以找到唯一的點z0，使F(z)在z≤z0時是z的單調(diào)遞增函數(shù)，在z≥z0時是單調(diào)遞減函數(shù)[17-18]。通過計算，得到z0=1.613630000000744。因此，對給定的時間t，由式(1)求解視電阻率ρ通常可獲得兩個解: 在區(qū)間(0,z0]獲得的解稱為“后全期解”，并稱相應(yīng)的方法為“后全期法”；在[z0,+∞)獲得的解稱為“前全期解”，并稱相應(yīng)的方法為“前全期法”。

從式(3)可以看出，z與ρt有關(guān)。如果t較大，則ρt較大，因而當z≤z0時用“后全期解”或“后全期法”較合理；反之，如果t較小，則應(yīng)該用“前全期解”或“前全期法”。

對實測的感應(yīng)電動勢V0(t)，由于噪聲和誤差的影響，可能導致利用式(1)求解視電阻率ρ時出現(xiàn)無解，這時可使用αV0(t)(α<1)代替V0(t)，重新求解視電阻率，參數(shù)α的定義和計算參見文獻[19]。

求解視電阻率ρ的后全期解與前全期解，獲得其精確解是瞬變電磁數(shù)據(jù)反演和解釋的關(guān)鍵環(huán)節(jié)。令

(4)

則ρ的求解問題轉(zhuǎn)化為求取非線性函數(shù)g(ρ)的解。由式(3)及F(z)的上凸特性，函數(shù)g(ρ)在(0,ρ0]內(nèi)單調(diào)遞增，在[ρ0,+∞)內(nèi)單調(diào)遞減，其中

(5)

因此視電阻率ρ的后全期解與前全期解可分別在[ρ0,+∞)與(0,ρ0]區(qū)間通過求g(ρ)的解獲得。

二分法、牛頓法是對非線性單調(diào)函數(shù)進行任意精確求根的方法，文獻[14-18]詳細描述了利用這兩種方法求解全期視電阻率的具體做法。弦截法也是對單調(diào)函數(shù)進行任意精確求根的方法，但弦截法與二分法相比，具有迭代次數(shù)少、收斂快的特點，而且無須象牛頓法那樣對函數(shù)求導。下面先介紹弦截求根法，然后給出全期視電阻率的弦截求解方法和步驟。

2 弦截求根法

設(shè)函數(shù)f(x)在區(qū)間[b,c]內(nèi)單調(diào)并且有根。不失一般性，不妨設(shè)f(x)為單調(diào)遞減函數(shù)。如圖1所示，任取初始值x1,x2∈[b,c]，逐步迭代計算f(x)過點(xi-1,f(xi-1))和(xi-2,f(xi-2))的弦截線AB在x軸的截距xi(i=3,4,…)，則xi將收斂于f(x)在區(qū)間[b,c]的唯一根x*，這就是弦截求根法的基本思想。

弦截線AB的表達式為

圖1 弦截法示意圖

設(shè)直線AB與x軸的截距為xi，則弦截法無限逼近求根的迭代公式為

(6)

3 全期視電阻率的弦截逼近算法

根據(jù)上述原理，對每個給定的時間點，都有兩個視電阻率ρ的全期解，分別為后全期解與前全期解。后全期解可在[ρ0,ρb]內(nèi)用弦截法求解，而前全期解可在[ρa,ρ0]內(nèi)用弦截法求解，其中ρb為視電阻率上界的一個估值，通常取ρb=1.0×106Ω·m;ρa為視電阻率下界的一個估值，通常取ρa=1.0×10-6Ω·m。

視電阻率ρ的后全期解的弦截解法步驟如下：

(1)任意給定視電阻率ρ的精度ε>0，任取ρ(1)、ρ(2)∈[ρ0,ρb]；令i=2，運用式(2)～式(4)計算g[ρ(i)]和g[ρ(i-1)]。

(2)令i=∶i+1，計算g[ρ(i)]。

(3)計算

ρ(i)=ρ(i-2)-g[ρ(i-2)]×

(4)若|ρ(i)-ρ(i-1)|<ε，則迭代停止，ρ(i)即為所求之解； 否則，轉(zhuǎn)步驟(2)繼續(xù)進行迭代計算。

前全期解計算步驟據(jù)此類推，不再贅述。

4 理論模型計算效果分析

視電阻率的計算方法可分為近似公式法和無限逼近法。近似公式法計算量小且算法穩(wěn)定，但得到的視電阻率誤差較大，尤其是在反映中淺層地質(zhì)狀況的早期采樣時間段，相對誤差有時可達5%以上。而無限逼近法可獲得誤差任意小的視電阻率，缺點是計算量大，穩(wěn)定性也不如近似公式法。二分法[14-17]、牛頓法[18]及本文提出的弦截法均為無限逼近法，理論上而言，牛頓法收斂最快，弦截法次之，二分法最慢[20-24]。牛頓法需要利用瞬變電磁信號的導數(shù)函數(shù)，對于式(1)～式(3)表示的地面瞬變電磁信號而言，其導數(shù)表達式過于復(fù)雜，其計算誤差影響也很大，因此，牛頓法并不適用于地面瞬變電磁法的視電阻率計算，而二分法和弦截法無需導數(shù)計算。下面通過兩個理論模型，比較二分法與弦截法的計算效果。

4.1 均勻半空間

假設(shè)均勻半空間的電阻率為80Ω·m。采用中心回線探測方式，發(fā)射天線有效面積為200m2，接收天線的有效面積為60m2，發(fā)射電流為2A，關(guān)斷時間為130μs，采樣間隔為20μs，共采集128個數(shù)據(jù)。通過式(1)～式(3)可獲得正演數(shù)據(jù)V(ti),i=1,2,…,128。

表1列出了在均勻半空間全期視電阻率計算過程中不同精度條件下，基于V(ti)分別用二分法和弦截法對所有點求解全期視電阻率所需的平均迭代次數(shù)。這里假設(shè)視電阻率ρ的變化范圍為1.0×10-8～1.0×106Ω·m， 弦截法的初始值ρ(1)=40Ω·m,ρ(2)=50Ω·m。

4.2 層狀模型

假設(shè)地下包含三個電性層，各層電阻率分別為50、8、70Ω·m，層厚分別為150、50m、+∞，其他各項參數(shù)同均勻半空間模型。表2為不同方法在不同精度要求下的迭代次數(shù)統(tǒng)計。

表1 均勻半空間模型二分法與弦截法計算視電阻率的平均迭代次數(shù)統(tǒng)計

表2 層狀模型二分法與弦截法計算視電阻率的平均迭代次數(shù)統(tǒng)計

對比表1與表2可以看出，在同等精度要求時，弦截法所需的迭代次數(shù)比二分法小得多。由于每次迭代兩種算法的運算量基本相當，因此弦截法所需的總運算量也比二分法少得多。

需特別指出的是，弦截法與牛頓法對初值都較敏感，實際計算后/前全期視電阻率時，應(yīng)在前述的步驟4中增加對ρ(i)是否仍落在區(qū)間[ρ0,ρb]或[ρa,ρ0]進行判斷； 如果超出范圍，則應(yīng)停止迭代，改用其他方法(比如近似公式法)。此外，在極后期或極前期，由于感應(yīng)電動勢太小且隨機噪聲較大，與二分法和牛頓法一樣，如果利用弦截法計算視電阻率，解的誤差會較大，這種情況下，采用近似公式法計算視電阻率，誤差會小一些且計算過程更穩(wěn)健。

5 應(yīng)用實例

實測數(shù)據(jù)采集自廣西州景礦區(qū)，地下300m處有一水倉。使用地面大定源回線TEM探測法，大回線設(shè)計在水倉的正上方，呈矩形狀，矩形邊長分別為500m和400m；在回線內(nèi)部設(shè)計4條測線，間距為30m，每條測線有8個測點，間距為25m；每個測點的二次場測道為62道，測道先密后疏；接收天線有效面積為100m2，發(fā)射電流為15A，采樣頻率為6.25Hz，關(guān)斷時間為250μs，數(shù)據(jù)疊加次數(shù)為16。任一測點與大回線任意一邊的最近距離均大于120m，因此可使用瞬變電磁二次感應(yīng)電動勢的中心回線公式(式(1))進行反演。

反演時對淺層盲區(qū)進行填充，視電阻率采用后全期解，且均使用弦截法進行求解，初值均取ρ(1)=30Ω·m,ρ(2)=40Ω·m,ρa=1.0×10-5Ω·m,ρb=1.0×105Ω·m，逼近精度ε=1.0×10-3。圖2為其中一條過水倉測線的視電阻率反演剖面，可見反演結(jié)果很好地展現(xiàn)了地下300m深度水倉的準確位置，甚至水倉的大致形狀及大小(水倉為低阻，即圖中深藍色區(qū)域)也清晰可見。

圖2 過水倉測線的視電阻率反演剖面

6 結(jié)束語

對瞬變電磁法中全期視電阻率的求解，給出了弦截無限逼近求解法及其具體的算法步驟。弦截求解法與二分求解法、牛頓求解法一樣，可用于求解后全期和前全期視電阻率任意精度解，但弦截法無需像牛頓法那樣計算瞬變電磁信號的導數(shù)，避免了復(fù)雜的導數(shù)計算，且所需的迭代次數(shù)比二分法少得多，具有更快的收斂速度。