姜國慶,程久龍,孫曉云,陳 丁,王海彥

(1.江蘇省地質(zhì)調(diào)查研究院基礎(chǔ)地質(zhì)研究所,江蘇南京 210018;2.中國礦業(yè)大學(xué)(北京)煤炭資源與安全開采國家重點實驗室,北京 100083)

全空間瞬變電磁全區(qū)視電阻率優(yōu)化二分搜索算法

姜國慶1,2,程久龍2,孫曉云2,陳 丁2,王海彥2

(1.江蘇省地質(zhì)調(diào)查研究院基礎(chǔ)地質(zhì)研究所,江蘇南京 210018;2.中國礦業(yè)大學(xué)(北京)煤炭資源與安全開采國家重點實驗室,北京 100083)

將優(yōu)化二分搜索算法應(yīng)用于全區(qū)視電阻率計算,以期提高計算精度和運算效率。以均勻全空間回線源中心的瞬變響應(yīng)為基礎(chǔ),探討全空間磁感應(yīng)強度垂向分量Bz及其時間變化率?Bz/?T核函數(shù)曲線特征,推導(dǎo)晚期視電阻率表達式。給出全空間條件下由?Bz/?T數(shù)值積分計算Bz的方法,并通過模擬計算對轉(zhuǎn)換精度進行分析。對瞬變延時T一定時Bz和?Bz/?T相對于電阻率ρ的單調(diào)性進行研究,結(jié)果表明Bz為ρ的單值函數(shù),滿足對ρ二分搜索的條件,而?Bz/?T為ρ的雙值函數(shù),其二分搜索的前提是輔助參數(shù)Z在變量ρ的取值區(qū)間內(nèi)恒滿足Z≤1.225,并在此基礎(chǔ)上對二分搜索的初始區(qū)間和具體流程進行了優(yōu)化。模型算例和應(yīng)用實例表明,該優(yōu)化算法計算速度快、精度高,并可同時獲得由Bz和?Bz/?T定義的全區(qū)視電阻率ρall和ρ′all。

全空間瞬變電磁法;中心回線裝置;全區(qū)視電阻率;二分搜索算法;優(yōu)化

全空間瞬變電磁法屬于近距離觀測方法,具有體積效應(yīng)相對較小、方向性強、分辨率高、對低阻體敏感、施工效率高等優(yōu)勢,因此,在礦井水文地質(zhì)勘查中得到了廣泛應(yīng)用[1-5]。但是相對而言全空間瞬變電磁法理論的研究相對滯后,以中心回線和重疊回線裝置為例,目前,這兩種常用裝置的視電阻率計算還多是采用半空間計算公式或經(jīng)校正后的晚期經(jīng)驗公式[1,6],已有全空間全區(qū)視電阻率的研究[7]也僅對感應(yīng)電動勢定義全區(qū)視電阻率進行了重點探討,但是理論和實踐[8-9]都證明磁場與感應(yīng)電動勢定義全區(qū)視電阻率在數(shù)值和解釋深度上都存在差異,并且在反映地下結(jié)構(gòu)和解釋深度方面磁場要明顯優(yōu)于感應(yīng)電動勢。

二分搜索算法是求非線性方程根的常用方法,因其算法穩(wěn)定、易于實現(xiàn)而得到廣泛的應(yīng)用。Peter George[10]在磁測數(shù)據(jù)反演中,為了求得著名的吉洪諾夫曲線中合適的折衷參數(shù)β,分別使用了插值法和二分法,結(jié)果,二分搜索算法因不要求曲線的曲率而更為穩(wěn)健;Leonard[11]在研究時間域電磁法探測未爆炸的爆炸物時,使用了二分法來求解超越方程的根qs,該參數(shù)控制著均勻場中導(dǎo)電球體的時間域二次場響應(yīng)的表達式。張成范等[12]考慮到二分法對初始隔根區(qū)間的限制低而使用二分搜索算法來計算均勻半空間表面矩形大定源回線瞬變電磁測深數(shù)據(jù)的全區(qū)視電阻率;李文堯等[13]通過核函數(shù)變量Z的二分法計算中心回線瞬變電磁全期視電阻率。陳清禮等[14]、劉俊等[15]采用對電阻率ρ的二分搜索計算感應(yīng)電動勢定義全區(qū)視電阻率,但未對函數(shù)單調(diào)區(qū)間及二分法初始區(qū)間進行深入探討。筆者將二分搜索算法應(yīng)用到全空間瞬變電磁法全區(qū)視電阻率計算中,同時計算Bz和?Bz/?T定義的全區(qū)視電阻率,并在深入探討核函數(shù)單調(diào)性的基礎(chǔ)上對算法流程進行優(yōu)化,使其運算速度更快,適應(yīng)性更強。

1 全空間中心回線裝置全區(qū)視電阻率

均勻全空間介質(zhì)中水平放置的圓形發(fā)射框中心的瞬變響應(yīng)(磁感應(yīng)強度垂向分量Bz及其時間變化率?Bz/?T)可以表示為

此時,電阻率可以表示為

式(1),(2)可表示為

其中,Y(Z)和Y′(Z)分別為全空間中心回線裝置Bz和?Bz/?T的核函數(shù)。圖1為根據(jù)式(6),(7)所繪制的核函數(shù)曲線,從圖1可看出,Bz的核函數(shù)Y(Z)為單值函數(shù),其相應(yīng)函數(shù)曲線為單調(diào)上升曲線,漸進值Y(Z→∞)≈1。由圖1可知當(dāng)實測磁場值Bz小于一次場值(μI/2a)時,由Bz總能求得單值的全區(qū)視電阻率值。?Bz/?T的核函數(shù)Y′(Z)為雙值函數(shù),其函數(shù)曲線存在一個轉(zhuǎn)折點Z=1.225,該轉(zhuǎn)折點處Y′(Z)取得極大值Y′max=0.41。由圖1可知當(dāng)實測歸一化Y′(Z)＞0.41時,全區(qū)視電阻率ρ無定義;當(dāng)Y′(Z)= 0.41時,全區(qū)視電阻率ρ存在惟一解;Y′(Z)＜0.41時,全區(qū)視電阻率ρ存在分別對應(yīng)早期和晚期部分的兩個解。

晚期條件下Z→0,通過對誤差函數(shù)和e-Z2進行泰勒級數(shù)展開,并省略Z的高次項可以得到

由式(3),(4),(8)可得到由Bz計算的晚期視電阻率表達式

圖1 全空間中心回線裝置核函數(shù)曲線

Fig.1 Kernel function curves of full space central loop unit

由式(3),(5),(9)可得到由?Bz/?T計算的晚期視電阻率表達式

2 磁感應(yīng)強度Bz的計算

瞬變電磁系統(tǒng)一般以固定的采樣頻率進行信號采集(如TerraTEM系統(tǒng)f=500 kHz),記錄數(shù)據(jù)為時窗范圍內(nèi)采樣值的算數(shù)平均值。假設(shè)某瞬變電磁系統(tǒng)具有N個時窗,對應(yīng)N+1個時窗邊界,定義每個時窗中心所對應(yīng)延遲時間分別為T1,T2,…,TN,時窗邊界所對應(yīng)延遲時間分別為TB1,TB2,…,TBN+1。根據(jù)電磁學(xué)原理Bz和?Bz/?T的轉(zhuǎn)換關(guān)系可以表示為

其中,Bz(TBN+1)為TBN+1時刻的磁場值,磁場計算的關(guān)鍵就是根據(jù)觀測?Bz/?T變化趨勢準(zhǔn)確估計TBN+1時刻的磁場值。由式(5),(9)可以推導(dǎo)近區(qū)時均勻全空間條件下?Bz/?T具有與均勻半空間相同的衰減規(guī)律,即:?Bz/?T∝T-5/2。通過計算可知全空間瞬變電磁法采用對數(shù)域線性回歸方法[8]計算Bz(TBN+1)的精度較高,這里結(jié)合具體算例對轉(zhuǎn)換精度進行分析。

假設(shè)均勻全空間電阻率為100 Ω·m,發(fā)射線圈為邊長2 m的方框(等效半徑a=2/π),發(fā)射電流為1 A,采樣頻率f=500 kHz,記錄時窗為Terra TEM系統(tǒng)的中間時間序列(intermediate time series)的前100個窗口(Windows),舍棄延遲時間小于10 μs的時窗。圖2中?Bz/?T曲線為式(5)理論值按f采樣所得,采用對數(shù)域線性回歸方法計算Bz(TBN+1)值為3.431 95×10-18T,由式(4)計算理論值為3.438 94× 10-18T,二者相對誤差僅為0.203%。由圖2中Bz積分轉(zhuǎn)換值與理論值對比可以看出,二者基本一致。同時,由相對誤差曲線可知所有時窗積分轉(zhuǎn)換的相對誤差均在0.6%以下,并且表現(xiàn)為兩端高、中間低的特點,分析晚期相對誤差高值是由Bz(TBN+1)的估算誤差引起,早期相對誤差高值是由時窗窄、采樣相對不足引起。

圖2 磁場Bz計算值與理論值對比曲線Fig.2 Contrastive curves of Bzcalculatedvalue and theoretical value

3 函數(shù)單調(diào)性探討

由式(1),(2)可知磁感應(yīng)強度垂向分量Bz(T, ρ)及其對時間變化率?Bz(T,ρ)/?T均為電阻率ρ的非線性復(fù)雜函數(shù),要通過對ρ的二分搜索來計算全區(qū)視電阻率就必須對兩函數(shù)相對于變量ρ的單調(diào)性進行探討。

3.1 Bz(T,ρ)的單調(diào)性探討

根據(jù)Bz的核函數(shù)曲線特征,Y(Z)為Z的單調(diào)增函數(shù)。針對特定測點的瞬變電磁觀測過程,發(fā)射線圈半徑a和發(fā)射電流I一定時,由式(4)可知Bz亦為Z的單調(diào)增函數(shù)。同時,對于某一瞬變延遲時間T,Z與ρ成反比,綜合以上可以判斷磁感應(yīng)強度Bz為電阻率ρ的單調(diào)減函數(shù),滿足對ρ的二分搜索來計算全區(qū)視電阻率的條件。如圖3(a)為磁場Bz隨電阻率ρ的變化曲線,由圖中可以看出對應(yīng)任意瞬變延遲時間T,Bz均為ρ的單調(diào)減函數(shù)。

3.2 ?Bz(T,ρ)/?T的單調(diào)性探討

根據(jù)?Bz/?T的核函數(shù)曲線特征,Y′(Z)為Z的雙值函數(shù),因此,針對特定測點觀測的某一瞬變延遲時間T(即參數(shù)a,I和T均為常數(shù)時),?Bz/?T與Y′(Z)相對變量Z具有相同的函數(shù)單調(diào)性,即均為Z的雙值函數(shù),分區(qū)間單調(diào)。同時,由Z與ρ的反比關(guān)系,可以判定在上述條件下?Bz/?T亦為電阻率ρ的雙值函數(shù)。為了滿足二分搜索算法對函數(shù)單調(diào)性的要求,下面對函數(shù)?Bz/?T的單調(diào)區(qū)間進行分析。

圖3 磁場Bz和磁場對時間導(dǎo)數(shù)?Bz/?T隨電阻率ρ的變化曲線Fig.3 Curves of Bzand?Bz/?T changing with resistivity ρ

由式(2),(5)可知?Bz/?T是電阻率ρ的非線性復(fù)雜函數(shù),通過前面對核函數(shù)Y′(Z)的探討,可以從理論上對函數(shù)單調(diào)特性的實質(zhì)進行分析。在參數(shù)a, I和T一定的條件下,?Bz/?T與Y′(Z)相對變量Z具有相同的函數(shù)單調(diào)性,因此,二者曲線具有相同的轉(zhuǎn)折點,即Z=1.225,以轉(zhuǎn)折點為界,Z＜1.225對應(yīng)瞬變延時的晚期部分,隨著Z值的減小?Bz/?T單調(diào)下降;Z＞1.225對應(yīng)瞬變延時早期部分,隨著Z值的增大?Bz/?T單調(diào)下降。由此,可以得到?Bz/?T相對于電阻率ρ的函數(shù)單調(diào)性的實質(zhì)是:在特定裝置(a為常數(shù))、選定瞬變延時(T為常數(shù))條件下,Z在變量ρ的取值區(qū)間內(nèi)是否恒滿足Z≤1.225或者Z≥1.225的條件,在實際勘探中筆者更關(guān)注瞬變延時晚期部分,即Z≤1.225的情況。

如果假設(shè)方形發(fā)射線框邊長L=2 m,全空間視電阻率ρ取值區(qū)間定為[10-2,104]Ω·m,由Z= aμ/4ρT 可以計算要使得Z在ρ的取值區(qū)間內(nèi)恒滿足Z≤1.225的條件,就必須使瞬變延時T≥26.66 μs。物理模擬實驗[16]和現(xiàn)場實測結(jié)果均表明多匝重疊小回線裝置的關(guān)斷時間一般為n×102μs,且關(guān)斷時間終點附近觀測數(shù)據(jù)會出現(xiàn)明顯震蕩,真正有效觀測數(shù)據(jù)的起始瞬變延時T遠大于上述臨界值,因此,二分搜索算法亦適用于?Bz/?T定義全區(qū)視電阻率的計算。圖3(b)為?Bz/?T隨電阻率ρ的變化曲線,可以看出,當(dāng)瞬變延遲時間T=0.001和0.01 ms時,?Bz/?T不滿足區(qū)間內(nèi)單調(diào)的條件,存在極大值點;當(dāng)瞬變延遲時間T≥0.1 ms時,?Bz/?T在ρ的取值區(qū)間內(nèi)單調(diào)下降;瞬變延遲時間T= 0.026 66 ms為臨界狀態(tài)。

4 優(yōu)化二分搜索算法的實現(xiàn)

二分搜索算法為數(shù)值分析中二分法的一種應(yīng)用,它的主要思想就是“取中”和“比較”。為了實現(xiàn)對電阻率ρ的二分搜索,首先必須確定其取值區(qū)間[ρa, ρb],文獻[14-15]采用固定區(qū)間[0.1,10 000],由文獻[17]可知常見巖、礦石的電阻率值變化區(qū)間為[10-n,10n],跨度達幾個數(shù)量級,因此,采用固定區(qū)間一方面可能導(dǎo)致視電阻率超出區(qū)間范圍使計算無解,另一方面會增加二分搜索運算量,因此,有必要對二分區(qū)間進行優(yōu)化。

已知Bz為ρ的單值函數(shù),由Bz求全區(qū)視電阻率的區(qū)間最小值ρa可以取為0;?Bz/?T為ρ的雙值函數(shù),ρa的選取必須滿足Z≤1.225的單調(diào)條件,因此,對于瞬變延時T,由Zmax=1.225可以計算其區(qū)間最小值ρZ=1.225。通過晚期核函數(shù)YL(Z)和Y′L(Z)的推導(dǎo)過程發(fā)現(xiàn),對于相同的Z值晚期核函數(shù)YL(Z)和Y′L(Z)均大于精確核函數(shù)Y(Z)和Y′(Z),從而由式(10),(11)所計算的晚期視電阻率值(ρL,ρ′L)總大于精確計算的全區(qū)視電阻率值(ρ,ρ′),并且隨著延遲時間T的增大,晚期值逐漸接近于全區(qū)值。據(jù)此,考慮采用晚期視電阻率(ρL,ρ′L)作為區(qū)間上限,從而極大地縮小搜索區(qū)間,提高運算速度。筆者將全空間Bz和?Bz/?T二分搜索的取值區(qū)間分別優(yōu)化為[0, ρL]和[ρZ=1.225,ρ′L]。以?Bz/?T計算全區(qū)視電阻率為例,經(jīng)優(yōu)化后的計算步驟如下:

(1)對T(i)時刻觀測值?Bz(i)/?T進行關(guān)斷時間校正。并由式(5),(7)計算Z=1.225時磁感應(yīng)強度變化率的理論極大值?Bzmax(i)/?T。

(2)比較?Bz(i)/?T與?Bzmax(i)/?T的相對大小。如果?Bz(i)/?T=?Bzmax(i)/?T,則全區(qū)視電阻率ρ′all(i)=ρ′Z=1.225;如果?Bz(i)/?T＞?Bzmax(i)/?T,則ρ′all(i)無解;如果?Bz(i)/?T＜?Bzmax(i)/?T,則進行下一步。

(3)由式(3)計算ρ′Z=1.225,由式(11)計算晚期視電阻率ρ′L,從而得到二分搜索初始取值區(qū)間ρa= ρ′Z=1.225,ρb=ρ′L。

(4)計算區(qū)間[ρa,ρb]的中點ρm=(ρa+ρb)/2,并將ρm帶入式(5),(7)計算中值電阻率所對應(yīng)的理論值?Bzm(i)/?T。

(5)計算觀測值?Bz(i)/?T與理論值?Bzm(i)/?T之間的相對誤差ε,如果相對誤差ε小于給定誤差限ε0,則ρ′all(i)=ρm,計算結(jié)束。否則,進行下一步。

(6)如果?Bz(i)/?T＞?Bzm(i)/?T則ρb=ρm,反之,ρa=ρm。計算搜索區(qū)間長度δ=ρb-ρa,如果δ小于給定最小區(qū)間長度δ0,則ρ′all(i)=ρm,計算結(jié)束。否則,返回步驟(4)繼續(xù)執(zhí)行。

5 模型算例與應(yīng)用實例

5.1 模型算例

均勻全空間條件下磁感應(yīng)強度垂向分量Bz及其時間變化率?Bz/?T具有解析解,并且求取的全區(qū)視電阻率具有標(biāo)準(zhǔn)值,因此最能檢驗計算結(jié)果的正確性和精度。這里選用均勻全空間模型對算法的精度及其運算效率進行分析。

假設(shè)均勻全空間電阻率為ρ=1 Ω·m,瞬變延遲時間T在對數(shù)域近似等間隔分布,共100個時窗。由式(10),(11)計算晚期視電阻率ρL和ρ′L;由二分搜索算法計算全區(qū)視電阻率ρall和ρ′all,誤差限ε0取為10-10。圖4為計算結(jié)果,由圖4可以看出,二分搜索算法計算的全區(qū)視電阻率ρall和ρ′all與模型電阻率ρ吻合較好,不同時窗計算全區(qū)視電阻率的相對誤差呈鋸齒狀波動,各時窗相對誤差最大值εmax僅為6.7× 10-8%。而晚期近似公式計算的ρL和ρ′L均大于全空間電阻率ρ,并且隨著T的增加ρL和ρ′L逐漸接近于全空間電阻率ρ,計算值與模型電阻率相對誤差在對數(shù)域呈近直線下降,最大相對誤差εmax為1%。同時,由于晚期視電阻率值與模型電阻率差異不大,采用晚期視電阻率ρL和ρ′L作為二分搜索區(qū)間上限可以有效減少搜索次數(shù),該模型計算結(jié)果表明,采用[0,ρL]和[ρZ=1.225,ρ′L]作為搜索區(qū)間與采用[0.01,100 000]相比分別可以減少二分搜索次數(shù)達32.1%和31.0%。

圖4 均勻全空間模型計算視電阻率及相對誤差曲線Fig.4 Apparent resistivity and relative error curves of homogeneous full space model

5.2 應(yīng)用實例

為查明某礦工作面截水巷道底板巖層富水性及含水構(gòu)造分布情況,開展了全空間瞬變電磁法探測研究,測線長度360 m,點距10 m,共37個測點。圖5為0號測點視電阻率曲線,從圖中可以看出全區(qū)視電阻率ρall和ρ′all的曲線類型一致,但是電阻率值存在差異,具體表現(xiàn)為在相對早期ρ′all明顯大于ρall,而在晚期階段二者重合或較為接近。同時可以看出ρ′all曲線變化較為劇烈,而ρall曲線變化則較為平緩。

理論上講,測點附近全空間介質(zhì)的電性結(jié)構(gòu)是確定的,無論是磁感應(yīng)強度垂向分量Bz還是其時間變化率?Bz/?T所定義的全區(qū)視電阻率應(yīng)該是一致的,但是實測結(jié)果表明全區(qū)視電阻率ρall和ρ′all存在差異,特別是早期階段,差異較為明顯。分析可知差異的產(chǎn)生主要有兩方面的原因:一是磁感應(yīng)強度垂向分量Bz及其時間變化率?Bz/?T在地下不同介質(zhì)界面處的表現(xiàn)不同。ρ′all對地下界面的反映比較敏感,但容易引起假薄層,ρall對地下界面反映則相對平緩。二是磁感應(yīng)強度垂向分量Bz是?Bz/?T經(jīng)過數(shù)值積分得到,存在轉(zhuǎn)換誤差,因此,由Bz和?Bz/?T計算的全區(qū)視電阻率ρall和ρ′all也存在差異[18-20]。

圖5 全空間瞬變電磁法0號測點視電阻率曲線Fig.5 Apparent resistivity curves of survey station number 0

完成該測線37個測點ρall和ρ′all計算總耗時不到2 min,可見該算法運算速度較快。從圖6可以看出ρall和ρ′all的斷面特征基本一致,斷面上以150,220, 330 m樁號為中心分別存在①,②和③3個低阻異常區(qū),其中以①號異常區(qū)的范圍和低阻特征最明顯,推斷該異常對應(yīng)范圍巖層富水性相對較強或存在含水構(gòu)造;②和③低阻異常區(qū)的范圍較小、視電阻率略低,推斷②和③異常對應(yīng)范圍巖層為弱富水,對工作面回采不構(gòu)成威脅。同時,可以看出ρall斷面上低阻異常區(qū)的頂界埋深均在20 m左右,而ρ′all斷面上埋深約為30 m,說明ρall的解釋深度要相對小于ρ′all。為了驗證物探成果的準(zhǔn)確性,礦方在巷道145 m處布置了一個底板探測孔,該鉆孔在23 m深度時開始出水,但水量不大,鉆孔最大出水量為3 m3/h,物探低阻異常區(qū)與底板巖層富水區(qū)位置吻合較好,并且ρall解釋深度與實際情況更接近。

圖6 全空間瞬變電磁法視電阻率斷面Fig.6 Apparent resistivity sections of full space TEM

6 結(jié) 論

(1)對于任意瞬變延時T,Bz均為ρ的單調(diào)減函數(shù),均滿足二分搜索算法的條件,而?Bz/?T為電阻率ρ的雙值函數(shù),二分搜索算法成立的條件為Z在變量ρ的取值區(qū)間內(nèi)恒滿足Z≤1.225;

(2)將Bz和?Bz/?T二分搜索的初始區(qū)間分別優(yōu)化為[0,ρL]和[ρZ=1.225,ρ′L],既可以解決巖、礦石電阻率值區(qū)間跨度大的問題,又可以減少運算量。

(3)全空間條件下磁感應(yīng)強度Bz可以通過?Bz/?T的數(shù)值積分求得,文中采用的對數(shù)域線性回歸方法計算精度相對較高。由Bz和?Bz/?T計算的全區(qū)視電阻率ρall和ρ′all的曲線類型一致,但在電阻率值和解釋深度存在差異。實際探測中,同時計算ρall和ρ′all,分別反演獲得擬斷面圖,并結(jié)合實際地質(zhì)資料對比分析,可以獲得更好的探測結(jié)果。

[1] 于景邨.礦井瞬變電磁法勘探[M].徐州:中國礦業(yè)大學(xué)出版社, 2007.

[2] 郭 純,劉白宙,白登海.地下全空間瞬變電磁技術(shù)在煤礦巷道掘進頭的連續(xù)跟蹤超前探測[J].地震地質(zhì),2006,28(3):456-462.

Guo Chun,Liu Baizhou,Bai Denghai.Predication of water disasters ahead of tunneling in coal mine using continuous detection by UTEM [J].Seismology and Geology,2006,28(3):456-462.

[3] 程久龍,邱 浩,葉云濤,等.礦井瞬變電磁法波場變換與數(shù)據(jù)處理方法研究[J].煤炭學(xué)報,2013,38(9):1646-1650.

Cheng Jiulong,Qiu Hao,Ye Yuntao.Research on wave-field transformation and data processing of the mine transient electromagenetic method[J].Journal of China Coal Society,2013,38(9):1646-1650.

[4] 于景邨,劉志新,劉樹才,等.深部采場突水構(gòu)造礦井瞬變電磁法探查理論及應(yīng)用[J].煤炭學(xué)報,2007,32(8):818-821.

Yu Jingcun,Liu Zhixin,Liu Shucai,et al.Theoretical analysis of mine transient electromagnetic method and its application in detec-ting water burst structures in deep coal stope[J].Journal of China Coal Society,2007,32(8):818-821.

[5] 姜國慶.礦井瞬變電磁法數(shù)據(jù)處理系統(tǒng)研究[D].青島:山東科技大學(xué),2011.

[6] 于景邨,劉樹才,王揚州.巷道內(nèi)金屬體瞬變電磁響應(yīng)特征及處理技術(shù)[J].煤炭學(xué)報,2008,33(12):1403-1407.

Yu Jingcun,Liu Shucai,Wang Yangzhou.Respounse characteristic of transient electromagnetic to metallic facilities in coal mines and the disposal technology[J].Journal of China Coal Society,2008,33 (12):1403-1407.

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

Yang Haiyan,Deng Juzhi,Zhang Hua,et al.Research on full-space apparent resistivity interpretation technique in mine transient electromagnetic method[J].Chinese Journal of Geophysics,2010,53(3): 651-656.

[8] 許建榮,李愛勇,楊 生.TEM中心回線法瞬變磁場求取和全區(qū)視電阻率計算[J].地質(zhì)與勘探,2008,44(6):62-68.

Xu Jianrong,Li Aiyong,Yang Sheng.Calculation of transient magnetic field and all time apparent resistivity based on central TEM loops method[J].Geology and Prospecting,2008,44(6):62-68.

[9] 白登海,Maxwell A Meju,盧 健,等.時間域瞬變電磁法中心方式全程視電阻率的數(shù)值計算[J].地球物理學(xué)報,2003,46(5): 697-704.

Bai Denghai,Maxwell A Meju,Lu Jian,et al.Numerical calculation of all-time apparent resistivity for the central loop transient electromagnetic method[J].Chinese Journal of Geophysics,2003,46(5): 697-704.

[10] Peter G L.Forward modeling and inversion of geophysical magnetic data[D].Columbia:The University of British Columbia,2003: 121-126.

[11] Leonard R P.Detecting unexploded ordnance with time domain electromagnetic induction[D].Vancouver:Simon Fraser University,1996:15.

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

Zhang Chengfan,Weng Aihua,Sun Shidong,et al.Computation of whole-time apparent resistivity of large rectangular loop[J].Journal of Jilin University(Earth Science Edition),2009,39(4):755-758.

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

Li Wenyao,Yan Chongwei.Bisection algorithm of central loop time domain transient electromagnetism all time apparent resistivity[J].Journal of Kunming University of Science and Technology(Natural Science Edition),2013,38(2):26-33.

[14] 陳清禮,嚴良俊,付志紅,等.長偏移距瞬變電磁法全區(qū)視電阻率的二分搜索數(shù)值算法[J].石油地球物理勘探,2009,44(6): 779-782.

Chen Qingli,Yan Liangjun,Fu Zhihong,et al.Searching algorithm for full time apparent in long-offset transient electromagnetic sounding[J].Oil Geophysical Prospecting,2009,44(6):779-782.

[15] 劉 俊,楊海燕,林 云.地下瞬變電磁法全區(qū)視電阻率二分搜索算法[J].工程地球物理學(xué)報,2012,9(5):526-530.

Liu Jun,Yang Haiyan,Lin Yun.Binary search algorithm of all-time apparent resistivity in underground transient electromagnetic method[J].Chinese Journal of Engineering Geophysics,2012,9(5): 526-530.

[16] 姜志海,岳建華,劉樹才.多匝重疊小回線裝置的礦井瞬變電磁觀測系統(tǒng)[J].煤炭學(xué)報,2007,32(11):1152-1156.

Jiang Zhihai,Yue Jianhua,Liu Shucai.Mine transient electromagnetic observation system of small multi-turn coincident configuration [J].Journal of China Coal Society,2007,32(11):1152-1156.

[17] 李金銘.地電場與電法勘探[M].北京:地質(zhì)出版社,2005.

[18] Fullagar P K.Generation of conductivity-depth pseudo-section from coincident and in-loop TEM data[J].Exploration Geophysics, 1995,26:43-55.

[19] Raiche A P.The use and misuse of apparent resistivity in electromagnetic methods[J].Geophysics,1984,15:37-41.

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

Optimized binary search algorithm of full space transient electromagnetic method all-time apparent resistivity

JIANG Guo-qing1,2,CHENG Jiu-long2,SUN Xiao-yun2,CHEN Ding2,WANG Hai-yan2

(1.Basic Geological Research Institute,Geological Survey of Jiangsu Province,Nanjing 210018,China;2.State Key Laboratory of Coal Resources and Safe Mining,China University of Mining and Technology(Beijing),Beijing 100083,China)

For the purpose of higher precision and efficiency,the optimized binary search algorithm was used for the calculation of all-time apparent resistivity.Based on the transient response of loop center in homogeneous full space, the kernel function curve feature of full space was discussed and the expression of late period apparent resistivity was derived.Method for calculating Bzfrom?Bz/?T by numerical integration in full space was given,and the transition precision was analyzed through simulation.The function monotonicity of Bzand?Bz/?T according to resistivity ρ at certain transient delay T was studied,and the results indicate that Bzis a single value function of ρ and it constantly satisfies the condition of binary search for ρ,but?Bz/?T is a double value function of ρ and the precondition for binary search is instrumental variable Z constantly meeting the formula Z≤1.225 during the value interval of variable ρ,on the basis of systematic study both the initial value interval and detailed flow chart of binary search algorithm were optimized.Model calculation and applied case indicate that the optimized binary search algorithm proposed in this paper has the characteristics of fast calculation speed,high calculation precision and can obtain all-time apparent resistivity ρalland ρ′allat the same time respectively from Bzand?Bz/?T,it has a better practicability and better detection effect in fullspace transient electromagnetic method.

full space transient electromagnetic method;central loop unit;all-time apparent resistivity;binary search algorithm;optimized

P631.3

A

0253-9993(2014)12-2482-07

2014-04-11 責(zé)任編輯:韓晉平

國家自然科學(xué)基金資助項目(51034003,51174210);高等學(xué)校博士學(xué)科點專項科研基金資助項目(20120023110014)

姜國慶(1986—),男,山東濟寧人,工程師,碩士研究生。E-mail:jiangguoqing220@126.com。通訊作者:程久龍(1965—),男,安徽安慶人,教授,博士生導(dǎo)師。E-mail:jlcheng@126.com

姜國慶,程久龍,孫曉云,等.全空間瞬變電磁全區(qū)視電阻率優(yōu)化二分搜索算法[J].煤炭學(xué)報,2014,39(12):2482-2488.

10.13225/j.cnki.jccs.2014.0477

Jiang Guoqing,Cheng Jiulong,Sun Xiaoyun,et al.Optimized binary search algorithm of full space transient electromagnetic method all-time apparent resistivity[J].Journal of China Coal Society,2014,39(12):2482-2488.doi:10.13225/j.cnki.jccs.2014.0477