王廣君,段雪影,黃鯨琿,胡祥云
(1.中國地質(zhì)大學(xué) 自動化學(xué)院,湖北 武漢 430074;2.中國地質(zhì)大學(xué) 復(fù)雜系統(tǒng)先進(jìn)控制與智能自動化湖北省重點實驗室,湖北 武漢 430074;3.中國地質(zhì)大學(xué) 地球物理與空間信息學(xué)院,湖北 武漢 430074;4.中國地質(zhì)大學(xué) 地質(zhì)過程與礦產(chǎn)資源國家重點實驗室,湖北 武漢 430074)
瞬變電磁法是一種高效率、無損、低成本的物探方法。通常通過發(fā)射機發(fā)射一個脈沖信號激勵地下地質(zhì)體產(chǎn)生二次場,然后對接收機測得的感應(yīng)電壓進(jìn)行反演從而得到地下介質(zhì)的信息。瞬變電磁法自1956年Wait[1]提出以來已經(jīng)成為了地表地質(zhì)調(diào)查中[2]最重要的方法。例如對銅、鉬、鉛鋅、鋁土礦、鈾等金屬資源的勘察[3];對巖土工程勘察中的采空區(qū)(空洞)、含水構(gòu)造、土石分界及巖層富水情況進(jìn)行地球物理勘查[4]。近年來,由于城市的快速發(fā)展,淺層地質(zhì)結(jié)構(gòu)的探測與成像在城市地下空間利用、防災(zāi)減災(zāi)等方面的作用越來越重要,城市地下空間的地質(zhì)構(gòu)造的精確勘察可以有效地避免地質(zhì)災(zāi)害和巨大的經(jīng)濟損失[5-7]。
傳統(tǒng)瞬變電磁法的淺層探測能力受限于發(fā)射機的關(guān)斷時間和接收機的過渡過程[8,9]。由于發(fā)射機的關(guān)斷時間不為零,在早期會激發(fā)一次場,導(dǎo)致探測到的早期信號是一次場和二次場疊加的總場。此外,非零關(guān)斷的磁場響應(yīng)與理想關(guān)斷的磁場響應(yīng)之間存在差異,導(dǎo)致后期計算視電阻率存在誤差,從而影響瞬變電磁法淺層探測中對二次場的反演解釋精度[10]。為此,孫天財?shù)韧茖?dǎo)了斜階躍場源激勵下的瞬變響應(yīng)解析式,并提出了一種不需要進(jìn)行關(guān)斷效應(yīng)校正直接計算全程視電阻率的方法,可減小視電阻率的誤差[11]。中南大學(xué)席振銖教授于2016年提出了等值反磁通瞬變電磁法[12],可以有效地剔除一次場的影響,該方法采用上下平行共軸兩組相同線圈為發(fā)射源,從而保證在該雙線圈合成的一次場零磁通平面上一次場的大小為零。在此基礎(chǔ)上,周超等提出了一種將等值反磁通瞬變電磁法和微動法聯(lián)合起來的聯(lián)合方法,為城市地下空間調(diào)查工作提供了一種行之有效的物探組合手段[13]。史存煥等提出了一種共軸偶極的瞬變電磁裝置方法,可減小多匝線圈的互感效應(yīng)[14]。但是目前對于瞬變電磁系統(tǒng)中的過渡過程還沒有充分的研究,在傳統(tǒng)的瞬變電磁法的淺層探測中,通常通過延時接收的方案來減小過渡過程對信號造成的畸變,但是這樣會造成淺層探測的盲區(qū)。此外,接收線圈的大小和匝數(shù)也會對過渡過程造成影響[15],所以有必要對瞬變電磁系統(tǒng)中的過渡過程進(jìn)行進(jìn)一步的研究。
本研究基于Simulink對瞬變電磁接收機的過渡過程進(jìn)行建模分析,通過與自制的瞬變電磁系統(tǒng)“WLZTEM-2022”的均勻半空間實測數(shù)據(jù)結(jié)合,在已知均勻半空間條件下反演得到接收線圈的等效電路的rLC參數(shù),即內(nèi)阻r、電感L、分布電容C,得到的rLC參數(shù)可以用于未知地下介質(zhì)淺層探測中過渡過程的矯正,且該方案易在接收機的FPGA(Field Programmable Gate Array)上實現(xiàn),實驗結(jié)果驗證了基于Simulink對瞬變電磁接收機的過渡過程建模的正確性,以及對淺層探測中過渡過程的矯正的必要性,并且提供了一種高精度的瞬變電磁法淺層探測的方案。
由于瞬變電磁每個接收線圈的分布參數(shù)都不一樣,且分布電容C的計算非常困難且復(fù)雜,在實際中,也沒有辦法直接測量線圈的分布電容,因此本文提出了一種全新計算rLC參數(shù)的方案,即在正演的過程中考慮過渡過程的影響,通過基于Simulink對接收機過渡過程的二階電路的建模分析,把接收線圈的rLC參數(shù)也加入到傳統(tǒng)的正演模型中。當(dāng)瞬變電磁系統(tǒng)的發(fā)射線圈和接收線圈大小和匝數(shù)設(shè)計好之后,通過傳統(tǒng)的正演模型和Simulink對接收機過渡過程的二階電路的建模,結(jié)合已知均勻半空間條件下的瞬變電磁系統(tǒng)“WLZTEM-2022”的實測數(shù)據(jù),通過蒙特卡洛反演得到二階電路的rLC參數(shù)(圖1)。所以本文的正演過程分為兩個部分:一個是基于傳統(tǒng)瞬變電磁法的非零關(guān)斷的電磁響應(yīng);另一個是基于Simulink對過渡過程二階電路的建模。
圖1 基于Simulink對瞬變電磁系統(tǒng)的過渡過程的rLC參數(shù)反演流程Fig.1 rLC parameter inversion flow of transient electromagnetic system transition process based on Simulink
針對淺層的瞬變電磁系統(tǒng)中發(fā)射機的關(guān)斷時間不可能嚴(yán)格等于零的問題,本文的解決方案是在瞬變電磁系統(tǒng)“WLZTEM-2022”的發(fā)射機中記錄發(fā)射電流I(t)的數(shù)據(jù),推導(dǎo)一次場的表達(dá)式,同時在傳統(tǒng)的瞬變電磁響應(yīng)的結(jié)果上做修正得到非零關(guān)斷電磁響應(yīng)的二次場信號。本課題組研制的瞬變電磁系統(tǒng)“WLZTEM-2022”的發(fā)射機發(fā)射電流的函數(shù)形式是e指數(shù)下降的I(t)=I0eα t。其中,I0為發(fā)射電流關(guān)斷時的大小(單位:A);α為衰減因子,決定關(guān)斷的時間,可以由設(shè)計的電路來調(diào)節(jié)。
非零關(guān)斷下一次場的電磁響應(yīng)Hp(單位:A/m)可以由Biotsafar定理得到:
(1)
式(1)中,d代表當(dāng)前積分元與單位向量dl之間的距離(m);θ是當(dāng)前積分元與單位向量dl之間的夾角(°);β為當(dāng)前積分元與接收線圈平面的夾角(°)。當(dāng)計算的感應(yīng)電壓在中心回線的時候,參數(shù)滿足θ=π/2,β=0且d=a,a為接收線圈的半徑(m),此時的一次場的大小為:Hp=I(t)/2a。在實際的應(yīng)用中需要記錄發(fā)射機的發(fā)射電流,以便理論計算一次場。
由于在實際應(yīng)用中,瞬變電磁系統(tǒng)中并不是直接測量磁場強度的,而是根據(jù)法拉第電磁定理,通過構(gòu)建磁環(huán)測量磁通量的變化,進(jìn)一步測量得到感應(yīng)電壓,因此計算得到一次場的感應(yīng)電壓Up(單位:V)為:
(2)
其中S為接收線圈的面積(m2);μ0為真空中的磁導(dǎo)率(H/m)。
(3)
其中,Hz(t-s)對應(yīng)理想關(guān)斷電磁響應(yīng)。對地下一維N層模型的理想關(guān)斷電磁響應(yīng)的理論計算已經(jīng)十分成熟。對于理想關(guān)斷響應(yīng)的計算,通常采用先計算頻域的響應(yīng),再通過傅里葉反變換轉(zhuǎn)換到時域的方法。對于半徑為a的發(fā)射線圈所激發(fā)的二次場,其對應(yīng)地表的中心回線處垂直磁場可以表達(dá)為:
(4)
公式(9)中的Kj是Gaver-Stehfest變換的系數(shù),wn對應(yīng)Hankel變換的系數(shù),h0和s分別對應(yīng)Hankel變換中對λ采樣的起始點和采樣間隔。最后結(jié)合杜哈梅爾積分(3)得到非零關(guān)斷下二次場的電磁響應(yīng):
(11)
(12)
其中,S對應(yīng)接收線圈的面積(m2)。實際中,接收機接收到的信號是一次場和二次場感應(yīng)電壓疊加的總場的信號:
U0(t)=Up(t)+Us(t)
(13)
下面一節(jié)會繼續(xù)討論接收機接收到的總場信號U0(t)經(jīng)過過渡過程之后的變化。
由于瞬變電磁儀的接收線圈中存在電感和分布電容,因此接收線圈在采集瞬變電磁響應(yīng)信號時會存在充電和放電的過程,從而導(dǎo)致采集到的瞬變電磁響應(yīng)信號與實際的瞬變電磁響應(yīng)信號存在一定的偏差,如圖2所示,該現(xiàn)象被稱之為過渡過程。
圖2 過渡過程示意圖Fig.2 Diagram of transition process
瞬變電磁系統(tǒng)的接收機需要使用接收線圈觀測瞬變電磁感應(yīng)的電壓,線圈和前置電路中的電容、電阻、電感構(gòu)成二階電路,如圖3(a)所示是簡化的接收線圈的等效電路: 其中R是阻尼電阻(Ω),r是接收線圈的內(nèi)阻(Ω),L為接收線圈的電感(H),C為接收線圈的電容(F,包含了分布電容),U0(t)是接收線圈的感應(yīng)電動勢(V),U(t)是信號經(jīng)過過渡過程之后的輸出電壓(V)。
圖3 簡化的接收線圈等效二階電路和Simulink仿真Fig.3 Simplified receiver coil equivalent second order circuit and Simulink simulation
根據(jù)基爾霍夫定理可得:
將式(15)帶入式(14)可得到微分方程:
(16)
式(16)可簡化為:
(17)
其中,ωb=ω0(1+r/R)1/2是電路的諧振頻率;ω0=(LC)-1/2是線圈的固有頻率;ζ=(r/(2L)+1/(2RC))/ωb是阻尼系數(shù);R為阻尼電阻(Ω)。當(dāng)ζ=1時,接收線圈工作與臨界阻尼狀態(tài);當(dāng)ζ>1時,為過阻尼狀態(tài);當(dāng)ζ<1時,為欠阻尼狀態(tài)。由阻尼系數(shù)可推導(dǎo)得到阻尼電阻:
(18)
本文基于Simulink仿真軟件求解二階電路的頻域方程,對式(14)做拉普拉斯變化后得到頻域的傳遞函數(shù)(圖3b):
(19)
式(19)表明過渡過程依賴于線圈的電阻r,電感L和電容C。由于電容的測量和計算非常困難,其中分布電容的值與線圈的大小、匝數(shù)、形狀、材質(zhì)等有關(guān),傳統(tǒng)的方法可以用諧振法[19]測量得到C參數(shù)的值,但是這種方法通常需要改變電路結(jié)構(gòu)用于觀察信號的波形,不利于電磁系統(tǒng)后期的調(diào)試和維護。本文提出一種新的測量方法,即用實測數(shù)據(jù)反演得到過渡過程對應(yīng)的特性參數(shù)。
傳統(tǒng)的瞬變電磁的數(shù)據(jù)反演中,反演模型對應(yīng)的是地下地質(zhì)體結(jié)構(gòu)的信息,由于未知參數(shù)較多,需要用到復(fù)雜的非線性反演的方法[20]。對于接收線圈分布參數(shù)的反演雖然也屬于非線性反演,但由于其反演的參數(shù)較少,所以可以用簡單的蒙特卡洛方法反演。對線圈的電阻r,電感L和電容C參數(shù)反演的步驟可以分為:
1)在特定的實驗條件下,設(shè)發(fā)射線圈的半徑a,接收線圈的面積S,地下均勻半空間的電導(dǎo)率ρ(Ω·m),理論正演計算得到非零關(guān)斷電磁響應(yīng)下的總場的感應(yīng)電壓U0T(ti), 同時記錄瞬變電磁系統(tǒng)“WLZTEM-2022”實測的發(fā)射電流I(t)和接收的總場的感應(yīng)電壓MT(ti),其中i=1,2,3 …,N代表測量的第i個數(shù)據(jù)點。
2)設(shè)定接收線圈分布參數(shù)的初始模型M=(r,L,C)。
3)基于Simulink過渡過程的二階電路建模得到矯正后的感應(yīng)電壓UT(ti)。
5)如果誤差函數(shù)滿足收斂條件ψ<ε,則反演過程結(jié)束;否則產(chǎn)生新的模型M=(r,L,C)+(Rand×r,Rand×L,Rand×C)返回步驟3)的計算。其中Rand表示隨機數(shù)。
實驗室自制的瞬變電磁系統(tǒng)“WLZTEM-2022”(圖4)包含發(fā)射機和接收機部分,瞬變電磁儀的探測深度與精度不僅取決于發(fā)射機的性能,也取決于接收機的數(shù)據(jù)采集能力。發(fā)射機的發(fā)射電流越大,最大探測深度越深,信噪比也越大,發(fā)射機的關(guān)斷時間越短,最小探測深度越淺;接收機的同步疊加次數(shù)越多、ADC(Analog-to-Digital Concerter)分辨率越高、放大器設(shè)計的越合理,也會使得信噪比越大,探測精度越高。
圖4 瞬變電磁系統(tǒng)“WLZTEM-2022”實物Fig.4 The real “WLZTEM-2022” transient electromagnetic system
發(fā)射機和接收機之間通過“同步”接口實現(xiàn)時間上的同步,通過同步單元將同步信號發(fā)送給接收機用于確定開始采樣的時刻。發(fā)射機中的電阻為可調(diào)電阻,用于改善發(fā)射電流波形,電流為可旋調(diào)節(jié)電流,逆時針旋轉(zhuǎn)電流變小,發(fā)射接口與發(fā)射線圈相連,接收機的通道與接收線圈相連,USB接口用于插入U盤存儲數(shù)據(jù)。該系統(tǒng)的操作流程如下:
1)打開開關(guān);
2)設(shè)置發(fā)射機和接收機的參數(shù),波形、頻率、占空比等;
3)設(shè)置發(fā)射電流,根據(jù)電流波形來調(diào)節(jié)匹配電阻旋鈕;
4)接收機采集數(shù)據(jù)并存儲;
5)后期的數(shù)據(jù)處理和反演。
利用上述實驗室自制的瞬變電磁系統(tǒng)“WLZTEM-2022”,于2021年8月31日在校內(nèi)進(jìn)行了淺層探測實驗,野外實驗的地點為中國地質(zhì)大學(xué)(武漢)西區(qū)數(shù)學(xué)物理學(xué)院后小樹林,該實驗采用中心回線裝置,發(fā)射線圈為5 m×5 m,接收線圈的面積為9 m2, 發(fā)射電流I0=6 A,采樣頻率為2.5 MHz,疊加次數(shù)為1 024次。小樹林的土壤在實驗當(dāng)天處于較濕狀態(tài),根據(jù)查閱的資料,較濕狀態(tài)下的土壤電阻率大約為100 Ω·m。
在上述野外實驗的條件下,實時記錄了瞬變電磁系統(tǒng)“WLZTEM-2022”的發(fā)射電流I(t)如圖5(a)所示。圖中顯示的是發(fā)射電流的對數(shù)logI(t)關(guān)于時間的t的實測數(shù)據(jù)。根據(jù)最小二乘法線性擬合得到發(fā)射電流的函數(shù)為:I(t)=6e-182 240t,瞬變電磁系統(tǒng)“WLZTEM-2022”的接收的實測數(shù)據(jù)MT(ti)如圖5(b)所示,數(shù)據(jù)的時間間隔為0.4 μs,由接收機的采樣率(2.5 MHz)確定。
圖5 瞬變電磁系統(tǒng)“WLZTEM-2022”野外實測數(shù)據(jù)Fig.5 Field measurement data of transient electromagnetic system “WLZTEM-2022”
在對過渡過程的rLC參數(shù)反演之前,需要理論計算的非零關(guān)斷電磁響應(yīng)U0(t),如圖6所示。圖6(a)是根據(jù)發(fā)射機發(fā)射的實際電流理論計算的一次場感應(yīng)電動勢Up(t),圖6(b)是根據(jù)野外實際的均勻半空間的電阻率ρ所對應(yīng)的非零關(guān)斷條件下的二次場感應(yīng)電動勢Us(t)。
圖6 理論計算的非零關(guān)斷電磁響應(yīng)(未經(jīng)過過渡過程)Fig.6 Theoretically calculated non-zero turn off electromagnetic reponse (without transition process)
本文考慮到了均勻半空間的電阻率ρ測量的不準(zhǔn)確性,在模型中計算了與理想值偏差±10 Ω·m 的結(jié)果。從圖6(b)中可以看出,均勻半空間的電阻率ρ的變化對二次場的影響很小。最后根據(jù)理論正演計算得到非零關(guān)斷電磁響應(yīng)下的總場的感應(yīng)電壓U0(t)=Up(t)+Us(t),挑選不同時段的理論數(shù)據(jù)U0T(ti)與接收機實測的數(shù)據(jù)MT(ti)反演過渡過程的rLC參數(shù)。
根據(jù)2.3節(jié)接收線圈分布參數(shù)的蒙特卡洛反演方法,得到不同的均勻半空間的電阻率ρ和數(shù)據(jù)點數(shù)N條件下的反演結(jié)果,如表1和圖7所示。表1中的誤差ψ0和誤差ψ分別表示未經(jīng)過過渡過程矯正前和經(jīng)過過渡過程矯正后的誤差函數(shù)。 從表中可以得知,矯正后的誤差ψ比誤差ψ0更小,矯正后的誤差定量上比矯正前的誤差約減小了3倍。 圖7(a)是未矯正的總場與實測數(shù)據(jù)的對比,從結(jié)果中可以看出,未經(jīng)過過渡過程處理的總場與實測數(shù)據(jù)在早期存在較大的差異,而且總場關(guān)于時間變化的函數(shù)在定性上就有很大的區(qū)別:未矯正的理論正演的曲線U0T(ti)關(guān)于時間是單調(diào)降低的;但是瞬變電磁系統(tǒng)“WLZTEM-2022”實測接收到的總場的感應(yīng)電壓MT(ti)在早期有一個上升的過程。圖7(b)是利用實測數(shù)據(jù)的前20個數(shù)據(jù)點反演的結(jié)果,反演得到的曲線與實測的數(shù)據(jù)在早期(t<5 μs)和晚期(27 μs
表1 不同的均勻半空間的電阻率ρ和數(shù)據(jù)點數(shù)N條件下的反演結(jié)果
圖7 過渡過程矯正前后的總場與實測數(shù)據(jù)的對比Fig.7 Comparison between the corrected total field by the transition process and measured data
對比不同數(shù)據(jù)長度N的反演結(jié)果可以發(fā)現(xiàn)數(shù)據(jù)長度N對電阻r的影響較小,數(shù)據(jù)長度N對電感和電容的值影響較大。對于淺層的探測應(yīng)該用早期的數(shù)據(jù)反演,從而得保證早期的信號與實測數(shù)據(jù)更吻合;對于傳統(tǒng)的深部探測應(yīng)該用晚期的數(shù)據(jù)反演,從而得保證晚期的信號與實測數(shù)據(jù)更吻合。此外數(shù)據(jù)長度N越大,反演得到的電感和電容就越小,根據(jù)電路理論分析,電容和電感越小,對應(yīng)二階電路的時間常數(shù)越小,信號變化越劇烈,這個結(jié)論也可以從反演的曲線得到,從而從理論自洽的方面驗證反演和過渡過程建模的正確性。
3.2節(jié)通過對接收線圈二階電路的建模得到了確定測量環(huán)境下rLC的實驗值,接下來用反演得到的rLC的值計算不同的地下層狀地質(zhì)模型的總場的響應(yīng),從而研究過渡過程對TEM淺層探測精度的影響。具體地,地下層狀地質(zhì)模型采用三層模型。第一層參數(shù)分別為h1=1 m,2 m,3 m,4 m,ρ1=1 000 Ω·m;第二層參數(shù)為h2=1 m,ρ2=1 Ω·m;第三層參數(shù)為ρ3=1 000 Ω·m。其模型的參數(shù)設(shè)置對應(yīng)于探測淺層埋藏的非電導(dǎo)體和電導(dǎo)體的差異。其中過渡過程的rLC參數(shù)取值取自上一節(jié)的蒙特卡洛反演得到的一組參數(shù)值:r=1.434 Ω,L=3.966e-6 H,C=2.891e-8 F。
圖8是過渡過程矯正前后的總場隨地下地質(zhì)模型(第一層探測深度h1)的變化。圖8(a)感應(yīng)電壓U0T(t)和UT(t)的對數(shù)圖,可以看出過渡過程矯正前后早期信號有較大差異,并且差異隨著時間的增加而減小。以理論曲線U0T(h1=4 m)和UT(h1=4 m)為例,兩條曲線在時間t>70 μs之后重合;同時從數(shù)值模擬的結(jié)果來看,過渡過程可以矯正早期70 μs內(nèi)的數(shù)據(jù)。這里值得注意的是,過渡過程對信號的矯正對整個測量時間內(nèi)的信號都會矯正,只是晚期的信號受過渡過程的影響較小。
圖8 過渡過程矯正前后的總場對淺層探測的影響Fig.8 The influence of the total field before and after the transition correction on shallow detection
另一方面從圖8(b)中得到過渡過程矯正前后淺層探測的靈敏度的變化。通常對于淺層探測的瞬變電磁法系統(tǒng)的接收機的測量時間的要求為小于1 μs。當(dāng)接收機能夠真實記錄早期信號(t<1 μs)的時候,測試的儀器才能夠滿足探測埋藏在淺層(0~5 m)的非電導(dǎo)體和電導(dǎo)體之間的差異。值得注意的是,在實際的測量中儀器測量得到的曲線是經(jīng)過過渡過程(rLC電路)影響之后的UT(t)曲線,其測量的結(jié)果(t<1 μs)對地質(zhì)模型(h1)的變化不敏感,此時的靈敏度幾乎為零,存在淺層探測(0~5 m)的測量盲區(qū)。然而如果對其數(shù)據(jù)通過反向的過渡過程的矯正,可以得到理論的總場的數(shù)據(jù)U0T(t)的曲線,可以看到曲線U0T(h1=1 m), 曲線U0T(h1=2 m),曲線U0T(h1=3 m)和曲線U0T(h1=4 m)在早期μs有明顯的差異,此時可以得到有效的可測量的靈敏度。綜上所述,理論上通過對測量數(shù)據(jù)的反向過渡過程的矯正可以提高瞬變電磁系統(tǒng)“WLZTEM-2022”的淺層(0~5 m)的探測能力。
1)通過理論數(shù)據(jù)與瞬變電磁測量的實驗數(shù)據(jù)對比,發(fā)現(xiàn)未經(jīng)過過渡過程處理的總場與實驗數(shù)據(jù)在早期有較大差異,驗證了對淺層探測中過渡過程的矯正的必要性。
2)通過對瞬變電磁儀器接收機的過渡過程建模及對接收線圈分布參數(shù)的蒙特卡洛反演,可以得到線圈的rLC參數(shù),經(jīng)過過渡過程后的理論曲線與實測的曲線在早期吻合,通過計算本文定義的誤差函數(shù)ψ,經(jīng)過過渡過程校正后的誤差定量上比矯正前的誤差約減小了3倍。
3)反演的結(jié)果表明反演的數(shù)據(jù)長度會影響反演的rLC參數(shù),反演的數(shù)據(jù)長度越大,反演得到的rLC參數(shù)越小,得到的感應(yīng)電壓曲線變化越劇烈;本文提出的反演算法可以直接移植到接收機的FPGA上。
4)把反演得到的過渡過程中rLC參數(shù)加入到傳統(tǒng)的TEM正演的算法中,可以得到對過渡過程準(zhǔn)確建模后的早期(t<1 μs)的理論總場曲線,從而提高瞬變電磁系統(tǒng)淺層(0~5 m)的探測能力。