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

曾慶寧 ，張法全， 張海如

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

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

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

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

1 全期視電阻率

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

(1)

其中

F(u)=u3e-u2

(2)

(3)

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

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

(4)

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

(5)

容易證明：F(u)在 0

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

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

(6)

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

(7)

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

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

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

2 牛頓求根法

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

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

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

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

(8)

設該直線與x軸的截距為xi，則

(9)

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

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

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

由式(9)可知，牛頓求解法中需用到函數的導數。因此由式(6)求解g(ρ)的根時，需用到g(ρ)的導數，容易得出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)即為所求之解，否則轉步驟2繼續進行迭代。

4 牛頓法與二分法的比較

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

通過試例，比較二分法與牛頓法兩者的收斂速度。

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

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

表1 二分法與牛頓法平均迭代次數

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

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

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

5 應用實例

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

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

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

6 結束語

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

參考文獻：

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

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

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

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

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

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

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

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