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

姜國(guó)慶,程久龍,孫曉云,陳 丁,王海彥

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

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

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

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

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

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

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

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

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

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

此時(shí),電阻率可以表示為

式(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)上升曲線,漸進(jìn)值Y(Z→∞)≈1。由圖1可知當(dāng)實(shí)測(cè)磁場(chǎng)值Bz小于一次場(chǎng)值(μI/2a)時(shí),由Bz總能求得單值的全區(qū)視電阻率值。?Bz/?T的核函數(shù)Y′(Z)為雙值函數(shù),其函數(shù)曲線存在一個(gè)轉(zhuǎn)折點(diǎn)Z=1.225,該轉(zhuǎn)折點(diǎn)處Y′(Z)取得極大值Y′max=0.41。由圖1可知當(dāng)實(shí)測(cè)歸一化Y′(Z)＞0.41時(shí),全區(qū)視電阻率ρ無(wú)定義;當(dāng)Y′(Z)= 0.41時(shí),全區(qū)視電阻率ρ存在惟一解;Y′(Z)＜0.41時(shí),全區(qū)視電阻率ρ存在分別對(duì)應(yīng)早期和晚期部分的兩個(gè)解。

晚期條件下Z→0,通過(guò)對(duì)誤差函數(shù)和e-Z2進(jìn)行泰勒級(jí)數(shù)展開(kāi),并省略Z的高次項(xiàng)可以得到

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

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

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

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

2 磁感應(yīng)強(qiáng)度Bz的計(jì)算

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

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

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

圖2 磁場(chǎng)Bz計(jì)算值與理論值對(duì)比曲線Fig.2 Contrastive curves of Bzcalculatedvalue and theoretical value

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5.1 模型算例

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

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

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

5.2 應(yīng)用實(shí)例

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

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

圖5 全空間瞬變電磁法0號(hào)測(cè)點(diǎn)視電阻率曲線Fig.5 Apparent resistivity curves of survey station number 0

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

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

6 結(jié) 論

(1)對(duì)于任意瞬變延時(shí)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ū)間跨度大的問(wèn)題,又可以減少運(yùn)算量。

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

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

[2] 郭 純,劉白宙,白登海.地下全空間瞬變電磁技術(shù)在煤礦巷道掘進(jìn)頭的連續(xù)跟蹤超前探測(cè)[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] 程久龍,邱 浩,葉云濤,等.礦井瞬變電磁法波場(chǎng)變換與數(shù)據(jù)處理方法研究[J].煤炭學(xué)報(bào),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] 于景邨,劉志新,劉樹(shù)才,等.深部采場(chǎng)突水構(gòu)造礦井瞬變電磁法探查理論及應(yīng)用[J].煤炭學(xué)報(bào),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] 姜國(guó)慶.礦井瞬變電磁法數(shù)據(jù)處理系統(tǒng)研究[D].青島:山東科技大學(xué),2011.

[6] 于景邨,劉樹(shù)才,王揚(yáng)州.巷道內(nèi)金屬體瞬變電磁響應(yīng)特征及處理技術(shù)[J].煤炭學(xué)報(bào),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é)報(bào),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] 許建榮,李愛(ài)勇,楊 生.TEM中心回線法瞬變磁場(chǎng)求取和全區(qū)視電阻率計(jì)算[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í)間域瞬變電磁法中心方式全程視電阻率的數(shù)值計(jì)算[J].地球物理學(xué)報(bào),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] 張成范,翁愛(ài)華,孫世棟,等.計(jì)算矩形大定源回線瞬變電磁測(cè)深全區(qū)視電阻率[J].吉林大學(xué)學(xué)報(bào)(地球科學(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é)報(bào)(自然科學(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] 陳清禮,嚴(yán)良俊,付志紅,等.長(zhǎng)偏移距瞬變電磁法全區(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é)報(bào),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] 姜志海,岳建華,劉樹(shù)才.多匝重疊小回線裝置的礦井瞬變電磁觀測(cè)系統(tǒng)[J].煤炭學(xué)報(bào),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] 李金銘.地電場(chǎng)與電法勘探[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] 牛之璉.時(shí)間域電磁法原理[M].長(zhǎng)沙:中南大學(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 ρa(bǔ)lland ρ′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é)任編輯:韓晉平

國(guó)家自然科學(xué)基金資助項(xiàng)目(51034003,51174210);高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目(20120023110014)

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

姜國(guó)慶,程久龍,孫曉云,等.全空間瞬變電磁全區(qū)視電阻率優(yōu)化二分搜索算法[J].煤炭學(xué)報(bào),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