井下瞬變電磁法中全期視電阻率的牛頓求解法

曾慶寧 ，張法全， 張海如

(桂林電子科技大學(xué) 信息與通信學(xué)院，桂林 541004)

近年發(fā)展起來的井下瞬變電磁法，因其相關(guān)設(shè)備體積小、成本低、地層分辨率高等優(yōu)勢，成為巷道超前探測地質(zhì)異常體的重要方法[1，7]。

井下瞬變電磁法原理與地面瞬變電磁法原理基本相同，但由于礦井瞬變電磁探測是在井下巷道內(nèi)進(jìn)行的，電磁場為全空間而非半空間。井下與地面一樣，主要通過反演地質(zhì)構(gòu)造體的視電阻率來發(fā)現(xiàn)異常地質(zhì)，而視電阻率的求解，可通過后期近似公式或前期近似公式進(jìn)行計算[1，7]，但難免存在誤差，而且在某些采樣時刻，這種誤差可能是無法容忍的。為此，可考慮使用后全期法和前全期法進(jìn)行精確求解[1-8]，但這種求解法需要使用復(fù)雜非線性方程求根技術(shù)。二分法是一種行之有效的非線性方程求根方法，已成功地用于對地面瞬變電磁法視電阻率的任意精度求解[2-3]，該方法亦可用于井下瞬變電磁法視電阻率的任意精度求解[1]。然而，二分法通常需要迭代的次數(shù)很多，收斂速度較慢。

本次研究引進(jìn)非線性方程求根的牛頓方法，并將其用于井下瞬變電磁法視電阻率的任意精度求解。通過實例，說明了牛頓法與二分法一樣，用于求解井下全期視電阻率是可行的，并且 在同等精度要求下，牛頓法比二分法所需的迭代次數(shù)要少許多，明顯具有更快的收斂速度，為井下瞬變電磁法基于視電阻率的地質(zhì)反演，提供了一種滿足任意精度要求且速度更快的運算方法。

1 全期視電阻率

在井下瞬變電磁法中，通常采用的是重疊回線或中心回線，其接收線圈于時刻t對二次場產(chǎn)生的感應(yīng)電壓為[1]

(1)

其中

F(u)=u3e-u2

(2)

(3)

由感應(yīng)電壓V(t)通過式(1)求解視電阻率ρ的方法稱為全期法。式(2)所示的F(u)稱為核函數(shù)，對其求導(dǎo)數(shù)可得

F′(u)=(3-2u2)u2e-u2

(4)

在(0,+∞)上F′(u)有唯一的零點

(5)

容易證明：F(u)在 0

對實測的感應(yīng)電壓V0(t)，由于各方面的噪聲與誤差效應(yīng)，有時難免導(dǎo)致式(1)求解視電阻率ρ可能無解，這時可使用αV0(t)代替V0(t)的方法(α<1)，通過重新定義視電阻率而獲得解決[4]。

求解視電阻率ρ的后全期解與前全期解，是瞬變電磁法進(jìn)行正確反演和資料解釋的關(guān)鍵環(huán)節(jié)，獲得視電阻率ρ的后全期精確值和前全期精確值是我們所期望的。對此，如果令

(6)

則ρ的求解問題轉(zhuǎn)化為函數(shù)g(ρ)的求根問題。由式(3)和F(u)的單調(diào)性，g(ρ)將在(0,ρ0]內(nèi)單調(diào)遞增，在[ρ0,+∞)內(nèi)單調(diào)遞減，其中

(7)

因此視電阻率ρ的后全期解與前全期解可分別在[ρ0,+∞)與(0,ρ0]中求g(ρ)的根而獲得。

二分法和牛頓法都是對單調(diào)函數(shù)進(jìn)行任意精確求根的常用方法。文獻(xiàn)[3]分別詳細(xì)描述了二分法求解半空間全期視電阻率的具體步驟，而全空間全期視電阻率的二分法算法步驟與半空間時是完全類似的[1]。

牛頓法與二分法相比，通常具有所需迭代次數(shù)少、收斂速度更快的特點。下面先介紹牛頓求根法，然后給出全期視電阻率的牛頓求解法。

2 牛頓求根法

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

圖1 牛頓法示意圖Fig．1 Diagrammatic sketch of Newton method

切線AB的斜率為f(x)在xi-1的導(dǎo)數(shù)f′(xi-1)，因此，切線AB的方程為

y=f′(xi-1)(x-xi-1)+f(xi-1)

(8)

設(shè)該直線與x軸的截距為xi，則

(9)

式(9)即為牛頓求根法的迭代公式。

3 全期視電阻率的牛頓求解法

根據(jù)上述第一節(jié)，對給定的時間，視電阻率ρ的全期解有兩個，分別為后全期解與前全期解。后全期解可在[ρ0,ρb]內(nèi)用牛頓法求解，而前全期解可在[ρa(bǔ),ρ0]內(nèi)用牛頓法求解，其中ρb為視電阻率上界的一個估值，通?？扇ˇ裝=105Ω·m，ρa(bǔ)為式電阻率下界的一個估值，通?？扇ˇ補(bǔ)=10-5Ω·m。

由式(9)可知，牛頓求解法中需用到函數(shù)的導(dǎo)數(shù)。因此由式(6)求解g(ρ)的根時，需用到g(ρ)的導(dǎo)數(shù)，容易得出g′(ρ)為

(10)

其中

F′(u)=(3-2u2)u2e-u2

(11)

(12)

(13)

于是視電阻率ρ的后(前)全期解的牛頓解法步驟可描述如下：

步驟1：任意給定視電阻率ρ的精度誤差ε>0，任取ρ(0)∈(ρ0,ρb)，令i=0。

步驟2：i=i+1，按式(6)、式(2)、式(3)計算g[ρ(i-1)]，并按式(10)、式(11)、式(12)、式(13)計算g′[ρ(i-1)]。

步驟3：按下式計算ρ(i):

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

4 牛頓法與二分法的比較

二分法與牛頓法均為對井下全期視電阻率進(jìn)行任意精確求解的方法。二分法原理簡單，實現(xiàn)容易，但收斂速度較慢；而牛頓法通常具有比二分法快得多的收斂速度。

通過試?yán)容^二分法與牛頓法兩者的收斂速度。

假設(shè)井下全空間為均勻介質(zhì)空間，介質(zhì)的視電阻率80 Ωm。采用同點探測方式，發(fā)射天線有效面積200 m2，接收天線有效面積60 m2，發(fā)射電流2A，關(guān)斷時間為130 μs，時間均勻采樣的采樣間隔為20 μs，測道數(shù)為128。于是，通過式(1)-式(3)可獲得正演數(shù)據(jù)V(ti),i=1、2、…、128。

表1列出了對全期視電阻率不同精度要求下，通過V(ti)，用二分法和牛頓法在所有測道求解后全期視電阻率所需的平均迭代次數(shù)。

表1 二分法與牛頓法平均迭代次數(shù)

在表1數(shù)據(jù)的計算過程中，假設(shè)視電阻率ρ在10-8Ω·m～106Ω·m之間，牛頓法的初始值取為ρ(0)=50。計算結(jié)果表明：二分法和牛頓法均能很好地收斂于給定的視電阻率。

從表1中容易看出，牛頓法僅需幾次迭代即可達(dá)到很高的精度，在同等精度要求下，牛頓法所需的迭代次數(shù)比二分法要少許多。

應(yīng)當(dāng)注意：牛頓法對函數(shù)導(dǎo)數(shù)及初值較為敏感，而且前述的g′(ρ)的絕對值在很多地方很小，因此，實際計算后(前)全期視電阻率時，應(yīng)在前述的步驟4中增加對ρ(i)是否仍落在區(qū)間[ρ0,ρb]([ρa(bǔ),ρ0])的判斷，如果超出區(qū)間范圍，則應(yīng)停止迭代。此外，牛頓法與二分法一樣，對特別后期(或特別前期)的采樣時刻，由于感應(yīng)電壓的值太小且實際測量中又難免隨機(jī)噪聲的影響，容易導(dǎo)致解的較大誤差，這時最好用后期近似公式法(或前期近似公式法)代替。

5 應(yīng)用實例

實測數(shù)據(jù)來源于對陜西陳家山煤礦某井下巷道掘進(jìn)頭處的TEM探測，目的是探測采煤掘進(jìn)前方是否有積水，以避免透水事故發(fā)生。探測時，分別對前方上傾15°扇面、水平方向扇面和下傾15°扇面進(jìn)行探測。每個扇面的探測范圍從左邊-45°至右邊45°，每15°一個測點，共7個測點。發(fā)射天線有效面積8 m2，接收天線有效面積40 m2，采用共軸方式，發(fā)射與接收天線相距7 m，發(fā)射電流2 A，采樣頻率25 Hz，每測點30測道，關(guān)斷時間為140 μs。

圖2為上述水平方向扇面反演剖面圖。反演中，視電阻率采用后全期解，且均使用牛頓法求解，初值均取ρ(0)=50 Ω·m，而ρa(bǔ)=10-5Ω·m，ρb=105Ω·m，ε=10-3。由圖2可見，在掘進(jìn)頭水平扇形探測面的右偏25°～31°、距離80 m～150 m存在低阻異常區(qū)域。事后證明，此區(qū)域正是該煤礦已經(jīng)充水的一個采空區(qū)，與反演結(jié)果完全吻合。

圖2 水平扇形反演圖Fig．2 Horizontal fanlight inversion picture

6 結(jié)束語

對井下全期視電阻率的求解，給出了牛頓求解法及其具體的算法步驟。牛頓求解法與二分法一樣，可用于求后全期和前全期任意精度要求的解，但牛頓法所需的迭代次數(shù)比二分法少得多，具有更快的收斂速度。

參考文獻(xiàn)：

[1] 楊海燕, 鄧居智, 張華,等. 礦井瞬變電磁法全空間視電阻率解釋方法研究[J]. 地球物理學(xué)報, 2010, 53(3): 651-656.

[2] 張成范,翁愛華,孫世棟,等.計算矩形大定源回線瞬變電磁測深全區(qū)視電阻率[J].吉林大學(xué)學(xué)報：地球科學(xué)版, 2009,39(4):755-758.

[3] 李文堯，晏沖為. 中心回線瞬變電磁法全期視電阻率的二分法求解[J]. 昆明理工大學(xué)學(xué)報：自然科學(xué)版, 2013, 38(2):26-33.

[4] 白登海，MAXWELL A MEJU, 盧健，等. 時間域瞬變電磁法中心方式全區(qū)視電阻率的數(shù)值計算[J]. 地球物理學(xué)報，2003,46(9)： 697-704.

[5] 翁愛華，陸冬華，劉國興. 利用連分式定義瞬變電磁法全區(qū)視電阻率研究[J]. 煤田地質(zhì)與勘探，2003，31(3): 56-59.

[6] 蘇朱劉，胡文寶. 中心回線方式瞬變電磁測深虛擬全區(qū)式電阻率和一維反演方法[J]. 石油物探，2002，41(2): 216-221.

[7] 李好，胡運兵，吳燕清. 應(yīng)用礦井瞬變電磁法超前探測煤礦井下含水體[J]. 礦業(yè)安全與環(huán)保，2012，39(5): 49-52.

[8] 牛之璉.時域電磁法原理[M].長沙：中南大學(xué)出版社，2007.