趙永成
(長江大學(xué)油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室, 武漢 430100)
隨著社會(huì)經(jīng)濟(jì)的持續(xù)發(fā)展,國內(nèi)外對(duì)油氣礦產(chǎn)資源的需求也與日俱增。對(duì)深地資源的迫切需求,地球深部鉆探技術(shù)與裝備、深部非常規(guī)能源開發(fā)以及深地空間資源開發(fā)均進(jìn)入了快速發(fā)展時(shí)期。井地電磁法在鉆井中布置激勵(lì)線源,該方法可以大大減少地表的電磁信號(hào)等有源噪聲的干擾,提高對(duì)油氣探測的精確度。相對(duì)于地震波勘探而言,井地電磁法更具有方便、快捷和成本低的優(yōu)點(diǎn),并且不會(huì)對(duì)自然環(huán)境造成破壞,有利于對(duì)目標(biāo)和相對(duì)低阻覆蓋區(qū)的探測,且勘查速度快、成本低、通行性好、可大區(qū)域快速勘探,尤其是在植被發(fā)育的覆蓋地區(qū),更具有一般勘探手段難以達(dá)到的效果。井地電磁法能更好地適應(yīng)我國山區(qū)大面積尋找礦產(chǎn)、干旱地區(qū)找地下水資源和油田剩余油氣的挖掘及其他方面的需要。井地電磁法的響應(yīng)分為頻率域和時(shí)間域。井地頻率域測井最先由俄羅斯科學(xué)家在20世紀(jì)70年代提出,頻率域測井主要是對(duì)井中線激勵(lì)源發(fā)射低頻信號(hào),在地表檢測電磁場信息,反演電磁響應(yīng)特征,提取地層信息,過去的幾年里該方法在俄國和我國新疆等多地區(qū)進(jìn)行了測試,取得了良好的效果[1]。井地電磁法時(shí)間域電磁測深一維正演理論比較成熟、水平電偶極子源,垂直磁偶極子源等瞬變電磁場國內(nèi)外研究得比較透徹,但由于實(shí)際條件、應(yīng)用等方面的局限性,對(duì)于垂直電偶源在特殊地層中瞬變電磁場國內(nèi)外研究尚少,對(duì)垂直電性源的場強(qiáng)變化還需進(jìn)一步研究。
最早研究井地電磁法的國家是俄羅斯和日本,并且他們研究出了發(fā)射和接收設(shè)備以用于油氣藏邊界的圈定。徐建華等[2]在1992年首次在文章中給出了層狀導(dǎo)電媒質(zhì)中垂直電偶極子的場,又在1994年提出多層介質(zhì)中偶極子場的系數(shù)遞推關(guān)系。針對(duì)計(jì)算難度較大的垂直線電源,何繼善在2001年首次給出了幾種不同的電條件下的垂直線電源電磁場的理論表達(dá)式,可用于地面電磁法和海洋電磁法的理論計(jì)算,王志剛等[3-4]在2007年采用積分方程的方法開發(fā)了井地電磁法三維正演的建模程序,其基于Born近似法的算法可以很好地?cái)M合實(shí)際儲(chǔ)層。羅毅翔[5]2011年對(duì)井中垂直長導(dǎo)線源層狀大地電磁響應(yīng)進(jìn)行了數(shù)值計(jì)算,采用D.Guptasarma和Singh提出的一種數(shù)字濾波法解決了貝塞爾函數(shù)的積分問題,計(jì)算結(jié)果和精度達(dá)到了要求。然后,對(duì)程序的正確性進(jìn)行了效驗(yàn)。采用幾種典型的大地模型,對(duì)井地電磁的響應(yīng)進(jìn)行了正演計(jì)算。李術(shù)才等[6]2014年建立多種典型的城市地鐵溶洞水槽模型,利用井間電阻率成像,可以保證數(shù)據(jù)的有效精度。曹輝等[7]2015年進(jìn)行了頻率域井地電磁法物理模擬,用實(shí)驗(yàn)的方法證明了井地電磁法對(duì)油氣區(qū)域范圍的確定。2018年張嘉琪[8]分析了井地電磁一維模型中影響垂直長導(dǎo)線源電磁場響應(yīng)的幾個(gè)關(guān)鍵因素,發(fā)現(xiàn)井地電磁法對(duì)高阻體反映更為靈敏,且當(dāng)發(fā)射源與異常體處于同一水平線或者測點(diǎn)位于井中時(shí),井地電磁法對(duì)異常體的反映更為明顯。2020年6月楊一凡[9]通過井地電位測試技術(shù)檢測出含油邊界,明確油水分布特征,規(guī)避低產(chǎn)風(fēng)險(xiǎn),優(yōu)選高阻、高含油區(qū)域。2021年4月趙云生[10]通過井地瞬變電磁三維正演軟件模擬了不同頻率下異常體的電磁場響應(yīng),研究了不同儲(chǔ)層規(guī)模、不同極化常熟下的異常體變化規(guī)律。2022年3月吳瓊等[11]根據(jù)畢奧-薩法爾定律推導(dǎo)出了自由空間圓形電流環(huán)的場強(qiáng)三分量表達(dá)式,針對(duì)常見的水平板狀體、垂直板狀體、傾斜板狀體三種不同的模型,進(jìn)行正演模擬,分析電流環(huán)場逼近純異常體的合理性和近似程度,為基于電流環(huán)的異常場特征和解釋方法的推演提供了借鑒。綜上,研究電磁場的分布特點(diǎn)是該勘探方法的基礎(chǔ),因此本文的垂直電性源時(shí)空分布特點(diǎn)的研究是必要的。
在準(zhǔn)靜態(tài)情況下,在有源區(qū)域的時(shí)間域麥克斯韋方程可表示為
(1)
?×H(r,t)=J(r,t)+Js(r,t)
(2)
?×[J(r,t)+Js(r,t)]=0
(3)
式中:r為位置矢量;E(r,t)、B(r,t)、H(r,t)和J(r,t)分別為r處、t時(shí)刻的電場強(qiáng)度磁感應(yīng)強(qiáng)度、磁場強(qiáng)度和傳導(dǎo)電流密度;Js(r,t)為r處、t時(shí)刻的外加源電流密度。基于宏觀介質(zhì)的本構(gòu)關(guān)系:
B=μH
(4)
(5)
(6)
為了將式(6)轉(zhuǎn)化為時(shí)域有限元(FETD)的近似方程,定義如下殘差矢量:
(7)
計(jì)算區(qū)域被離散成若干個(gè)有限元。在加權(quán)積分的意義上,每個(gè)元素的殘差矢量被強(qiáng)制為零。在幾種不同形狀的有限元單元中,選擇可以有效地處理復(fù)雜地質(zhì)模型的四面體單元。
設(shè)計(jì)算區(qū)域內(nèi)的殘差加權(quán)積分為
?ΩW(r)R(r,t)dV=0
(8)
式中:Ω為加權(quán)系數(shù)。加權(quán)余量法通過最小化權(quán)函數(shù)W(r)和殘差矢量R(r,t)的內(nèi)積來尋找最優(yōu)解,即W(r)和R(r,t)正交。將式(7)代入式(8)可得
(9)
應(yīng)用第一矢量格林定理,將式(9)第一項(xiàng)積分展開:
(10)
使用非結(jié)構(gòu)四面體網(wǎng)格離散計(jì)算區(qū)域,采用自動(dòng)滿足電場切向分量連續(xù)且無散的Nedelec矢量插值基函數(shù)近似單元內(nèi)電場線性分布[12],對(duì)于四面體單元內(nèi)任意位置的電場可以表示為
(11)
式中:Ej(t)為四面體單元第j條棱邊的切向電場,可稱為自由度;N(r)為矢量插值基函數(shù)。基于Galerkin法,當(dāng)采用矢量插值基函數(shù)N(r)作為加權(quán)系數(shù)時(shí),內(nèi)邊界面兩側(cè)的面積分相互抵消,外邊界距離發(fā)射源足夠遠(yuǎn),滿足Sommerfeld邊界條件,在忽略面積分后,結(jié)合(11)電場方程滿足:
(12)
將所有單元加權(quán)殘差累加,可寫成矩陣:
(13)
式中:S和M為整體剛度矩陣和質(zhì)量矩陣;J為源項(xiàng),具體表達(dá)式如下:
(14)
(15)
(16)
對(duì)于接地長導(dǎo)線源,將其分割為若干個(gè)小段,每段近似為電偶極子[13],將電偶極子置于四面體網(wǎng)格棱邊上,可以有效近似實(shí)際接地長導(dǎo)線源。對(duì)于每個(gè)電偶極子,其電流密度可表示為
Js(r,t)=js(t)uδ(r-rs)
(17)
式中:js(t)為電流密度強(qiáng)度;u為電流方向向量;rs為發(fā)射電偶極源的位置;δ(r)是脈沖函數(shù),可表示為
(18)
而js(t)可表示為[14]
js(t)=I(t)dl
(19)
式中:I(t)為t時(shí)刻電流強(qiáng)度;dl為電偶極子長度。
聯(lián)立式(13)~式(16),采用二階后推歐拉近似式中i+2時(shí)刻電場與電流密度對(duì)時(shí)間的導(dǎo)數(shù)[15]
(20)
(21)
將式(20)和式(21)分別代入式(13)和式(16)中可得:
(3M+2ΔtS)E(i+2)(t)=
M[4E(i+1)(t)-E(i)(t)]-2ΔtJ(i+2)(t)
(22)
2ΔtJ(i+2)(t)=3J(i+2)(t)-4J(i+1)(t)+J(i)(t)
(23)
由式(23)即可得到任意發(fā)射源電流波形的電流源項(xiàng),式(22)可簡寫為
KE=b
(24)
式中:K為系數(shù)矩陣;E為棱邊上的未知電場;b為已知發(fā)射源項(xiàng)。
如圖1所示,參考設(shè)計(jì)了均勻半空間條件下的井地瞬變電磁模型,發(fā)射源A極位于井口,坐標(biāo)為(0,0,0),B極分別位于900 m,1 000 m,1 100 m,發(fā)射電流為1 A,空氣電阻率為1×108Ω·m,背景電阻率為100 Ω·m,異常體位于坐標(biāo)軸的y方向,大小為200 m×200 m×100 m,電阻率為10 Ω·m,異常體中心點(diǎn)距離井孔500 m。計(jì)算地面接收點(diǎn)為(0,200,0),(0,500,0),(0,1 000,0)的所有電磁場分量的響應(yīng)。當(dāng)接收位置在(0,500,0)分別計(jì)算發(fā)射源B極位于900 m,1 000 m,1 100 m時(shí)電磁場分量Ex的響應(yīng)結(jié)果,如圖2所示。從計(jì)算結(jié)果可以看出,當(dāng)接收位置不變,只改變發(fā)射源的長度時(shí),Ex的響應(yīng)結(jié)果曲線形態(tài)不會(huì)發(fā)生變化,只是數(shù)值大小會(huì)有所改變。
圖1 理論驗(yàn)證模型示意圖
圖2 發(fā)射源500 m時(shí)電磁場分量(Ex)響應(yīng)結(jié)果
為了驗(yàn)證本文所用程序的正確性和可行性,按照?qǐng)D1中設(shè)計(jì)的發(fā)射源長度為1 000 m,接收位置分別為(0,200,0)、(0,500,0)、(0,1 000,0)時(shí)計(jì)算電磁場分量Ex的響應(yīng)結(jié)果,其計(jì)算結(jié)果與數(shù)值解進(jìn)行對(duì)比,對(duì)比結(jié)果如圖3所示。圖3(a)、圖3(b)、圖3(c)分別為接收點(diǎn)(0,200,0)、(0,500,0)、(0,1 000,0)時(shí)的Ex響應(yīng)值對(duì)比結(jié)果,從中可以看到兩種計(jì)算結(jié)果的擬合情況比較好,可以證明本文垂直電性源算法程序的正確性和可靠性。
圖3 理論驗(yàn)證數(shù)據(jù)對(duì)比
為進(jìn)行井地瞬變電磁法正演模擬,設(shè)計(jì)背景電阻率為100 Ω·m,無異常體情況下,均勻半空間模型,如圖4所示。發(fā)射源兩端位于地下900~1 400 m處,垂直長導(dǎo)線源在地面投影位于(0,0,0),接收測網(wǎng)位于地面橫坐標(biāo)是500~2 500 m,測點(diǎn)間隔為50 m,縱坐標(biāo)是-1 000~1 000 m,測點(diǎn)間隔為50 m,觀測時(shí)間分別為0.1 ms、5.17 ms、25 ms、50 ms。其不同時(shí)刻Ex和Ey的響應(yīng)結(jié)果如圖5和圖6所示,圖5和圖6中(a)、(b)、(c)、(d)分別對(duì)應(yīng)0.1 ms、5.17 ms、25 ms、50 ms。按照結(jié)果分析,電場分量Ex響應(yīng)圖像展現(xiàn)了對(duì)稱性,隨著電場的擴(kuò)散,極大值越來越小。表明隨著時(shí)間的推移,場值在減小,電場分量Ey在平面呈現(xiàn)了隨時(shí)間的分布情況,Ey擁有對(duì)稱性,極大值始終在原點(diǎn)發(fā)射源四周,隨時(shí)間推移Ey在減小。
圖4 均勻半空間模型收發(fā)位置示意圖
圖5 不同時(shí)刻Ex響應(yīng)分布
圖6 不同時(shí)刻Ey響應(yīng)分布
圖7 收發(fā)位置示意圖
同樣模型參數(shù)不變,只改變觀測平面的位置,再次計(jì)算其結(jié)果如圖8所示。圖8(a)、圖8(b)、圖8(c)、圖8(d)分別為0.1 ms、1 ms、5.17 ms、25 ms時(shí)刻電場分量Ex的響應(yīng)結(jié)果示意圖??梢娫趥?cè)面觀察垂直電性源Ex場強(qiáng)分布的特點(diǎn)是關(guān)于橫軸上下對(duì)稱,且隨著時(shí)間場強(qiáng)也是逐漸擴(kuò)散。
圖8 不同時(shí)刻Ex響應(yīng)分布
本文設(shè)計(jì)并分析均勻半空間模型的垂直電性源瞬變電磁場時(shí)空分布特征,驗(yàn)證了垂直電性源地井瞬變電磁正演程序的可行性,通過均勻半空間下的正演模擬,得到垂直電性源瞬變電磁場的平面分布特征及其隨著時(shí)間場值的擴(kuò)散規(guī)律,探究在本文觀測系統(tǒng)下的有效觀測區(qū)域與優(yōu)勢觀測分量。主要結(jié)論如下。
(1)分析不同時(shí)刻、不同深度的均勻半空間瞬變電磁場時(shí)空響應(yīng)平面分布特征,可知電場Ex分量響應(yīng)效果最佳。
(2)模擬結(jié)果表明電磁場分量Ex呈對(duì)稱分布,水平接收時(shí)以x軸對(duì)稱,垂直接收時(shí)以z軸對(duì)稱,隨時(shí)間逐漸向周圍擴(kuò)散,并且場值逐漸減小。同時(shí)Ey的響應(yīng)圖像也擁有對(duì)稱性,極大值始終在原點(diǎn)發(fā)射源四周,隨時(shí)間推移Ey在減小。