圓柱坐標(biāo)系下三維瞬變電磁過套管響應(yīng)特征研究

蔣海蓉 趙 寧*② 趙鴻圖 秦 策 肖占山 米曉利

(①河南理工大學(xué)物理與電子信息學(xué)院，河南焦作 454000; ②河南省瓦斯地質(zhì)與瓦斯治理重點實驗室——省部共建國家重點實驗室培育基地，河南焦作 454000; ③中國石油集團(tuán)測井有限公司，陜西西安 710077； ④東方地球物理公司綜合物化探處，河北涿州 072751)

0 引言

隨著油氣勘探的不斷深入，勘探目標(biāo)從傳統(tǒng)的構(gòu)造油氣藏逐步轉(zhuǎn)向復(fù)雜油氣藏和隱蔽油氣藏，勘探難度增加，對地球物理勘探技術(shù)的要求也越來越高。距井眼較遠(yuǎn)區(qū)域的地質(zhì)構(gòu)造、儲層橫向展布、裂縫發(fā)育及斷層分布等信息對于油氣藏探測具有重要意義。近些年，基于瞬變電磁理論的井孔時域電磁測井方法克服了傳統(tǒng)頻率域電磁測井儀器通過降低頻率、增加源距實現(xiàn)遠(yuǎn)距離探測的技術(shù)局限性，通過測量電流脈沖關(guān)斷過程激勵的電磁衰變信號和晚期電磁場，采用較短的源距就可以達(dá)到頻率域長源距才能實現(xiàn)的探測效果[1]。

瞬變電磁法經(jīng)過近四十年的發(fā)展，取得了不錯的效果，但目前該方法在測井領(lǐng)域的應(yīng)用還處于探索階段。Kaufman等[2]首先提出了時域感應(yīng)測井理論； Anderson等[3]利用FFT變換，將頻率域測井響應(yīng)轉(zhuǎn)換到時間域，模擬了井眼、薄層等因素的瞬態(tài)響應(yīng)； 周仕新等[4]對磁偶極源激發(fā)的井中瞬變電磁進(jìn)行了模擬； Banning等[5]通過井中瞬變電磁三分量感應(yīng)測井?dāng)?shù)值模擬及成像技術(shù)對地層傾角及視電阻率進(jìn)行了估算； Onegova[6]討論了垂直井眼環(huán)境下泥漿電阻率和線圈陣列長度發(fā)生變化時的裝置相對偏心率與視電阻率之間的關(guān)系； Hagiwara[7]采用三軸感應(yīng)測井方式，基于三軸響應(yīng)比值估算了地層傾角和各向異性； 李建慧等[8]開展了地—井瞬變電磁三維正演； Commer等[9]運用時域有限差分法，在直角坐標(biāo)系下進(jìn)行三維建模，計算了鋼套管條件下的瞬變電磁響應(yīng)； Swidinsky等[10]采用重疊回線模式的發(fā)射、接收線圈，研究了井中瞬變電磁波響應(yīng)特征； 李建慧等[11]采用非結(jié)構(gòu)化四面體網(wǎng)格對局部區(qū)域進(jìn)行網(wǎng)格加密，并運用矢量有限元方法計算了復(fù)雜回線源的脈沖響應(yīng)； 段書新等[12]研究了井中瞬變電磁法對井旁低阻介質(zhì)的探測能力； Meng等[13]基于全空間瞬變電磁場表達(dá)式和Halley優(yōu)化算法提出了一種視電阻率校正算法，對徑向?qū)訝钅Ｐ瓦M(jìn)行模擬，結(jié)果證明該算法可以消除徑向低阻內(nèi)層對瞬變電磁測井響應(yīng)的影響； 陳衛(wèi)營等[14]分析了三維地—井瞬變電磁法的電場和磁場全分量探測效果； 郭建磊等[15]基于有限差分法實現(xiàn)了地—井瞬變電磁三維各向異性正演，并分析了半空間模型和塊狀異常體模型的三分量響應(yīng)特征。

近些年來，瞬變電磁法數(shù)值模擬技術(shù)日趨成熟，目前常用的方法包括有限元法[16-17]、有限差分法[18-19]、積分方程法[20-21]和有限體積法[22-23]等。針對這些瞬變電磁數(shù)值模擬方法，國內(nèi)外學(xué)者大多在笛卡爾坐標(biāo)系中進(jìn)行數(shù)值模擬，但在測井問題中，井孔及地層非均勻介質(zhì)的電導(dǎo)率在小尺度(套管、泥漿等)上存在劇烈變化。傳統(tǒng)笛卡爾坐標(biāo)系為了能夠反映出小尺度上的介質(zhì)參數(shù)變化，需要采用很細(xì)的剖分網(wǎng)格，這大大增加了計算區(qū)域的網(wǎng)格數(shù)量，通常情況下可達(dá)千萬級(107)甚至超過億級(108)。雖然當(dāng)前計算機(jī)計算能力和性能得到了迅速提升，但若不改進(jìn)建模方式，仍然難以滿足實際計算需求。而圓柱坐標(biāo)系能較好地模擬井眼及周圍小尺度上的介質(zhì)屬性變化，可以有效減少網(wǎng)格剖分?jǐn)?shù)量，提升計算效率。因此，本文在圓柱坐標(biāo)系下開展模擬。

瞬變電磁法中發(fā)射波形是影響探測效果的重要指標(biāo)之一，常用的發(fā)射波形有梯形波、半正弦波、三角波等，發(fā)射波形可分為開啟(on-time)、穩(wěn)定、關(guān)斷(off-time)三個階段。目前，許多學(xué)者將發(fā)射波形相關(guān)參數(shù)對瞬變電磁響應(yīng)結(jié)果的影響作為研究重點[24]。齊彥福等[25]通過改變時間道上的瞬時電流強(qiáng)度，并運用有限元方法完成了三維航空電磁全波形正演模擬與響應(yīng)特征分析。趙越等[26]利用梯形波、三角波和半正弦波對淺海瞬變電磁場進(jìn)行了數(shù)值模擬，并分析了不同波形在開啟和關(guān)斷下的響應(yīng)特征。Zeng等[27]采用有限元方法，對不同電導(dǎo)率的半空間進(jìn)行了全波形正演模擬。Ji等[28]運用有限差分法分析了線性變化的斜坡關(guān)斷時間對淺層異常體的全波瞬變電磁響應(yīng)。

綜上所述，在井孔地層非均勻介質(zhì)條件下，如何有效減少網(wǎng)格數(shù)量進(jìn)而降低方程組的求解規(guī)模，以及如何設(shè)置合理的發(fā)射波形及相關(guān)參數(shù)，對目前瞬變電磁方法的研究具有重要意義。本文在圓柱坐標(biāo)系下對求解區(qū)域進(jìn)行離散，基于有限體積法實現(xiàn)了套管井三維瞬變電磁正演模擬[29]，重點對套管電導(dǎo)率、相對磁導(dǎo)率、厚度和內(nèi)徑變化對響應(yīng)結(jié)果的影響進(jìn)行了研究，對比了梯形波、三角波、半正弦波的響應(yīng)特征，并分析了梯形波各階段持續(xù)時間對電磁響應(yīng)特征的影響。

1 正演理論

1.1 控制方程

在準(zhǔn)靜態(tài)情況下，時間域麥克斯韋方程組可表示為

(1)

?×H-J=S2

(2)

式中：E表示電場強(qiáng)度；B表示磁感應(yīng)強(qiáng)度；H表示磁場強(qiáng)度；J表示電流密度；S1表示外加回線源的磁流密度，下標(biāo)“1”代表磁性源；S2表示外加線源的電流密度，下標(biāo)“2”代表電性源；t表示時間。上式中場量之間的本構(gòu)關(guān)系可表示為

J=σE

(3)

B=μH

(4)

式中σ和μ分別表示電導(dǎo)率和磁導(dǎo)率。在時間域進(jìn)行麥克斯韋方程離散，需要調(diào)用式(3)和式(4)的本構(gòu)關(guān)系表示離散系統(tǒng)。時域問題的電磁場公式為

(5)

(6)

其中

(7)

將式(7)代入式(5)可得

(8)

對上式進(jìn)行變形得到

(9)

式(5)～式(9)中：C表示邊旋度算子； 矩陣M表示內(nèi)積算子； 上標(biāo)“e”和“f”分別代表網(wǎng)格單元的邊和面； 上標(biāo)“T”表示轉(zhuǎn)置。

1.2 有限體積離散

有限體積法將待求解區(qū)域離散為互不重疊的體積單元，并對需要求解的微分方程進(jìn)行體積積分。本文采用交錯網(wǎng)格實現(xiàn)有限體積空間離散，將電場分量定義在網(wǎng)格單元的棱邊中心，磁場分量定義在網(wǎng)格單元的面中心?；趫A柱坐標(biāo)系可方便進(jìn)行井眼建模，不會在井眼部位產(chǎn)生大量的離散網(wǎng)格，可有效提高計算速度。圓柱坐標(biāo)系下網(wǎng)格單元離散形式如圖1所示，其中圖1a所示為未進(jìn)行角度離散的完全對稱網(wǎng)格單元，圖1b為包含角度離散的三維網(wǎng)格單元。本文選用圖1a所示剖分形式對求解區(qū)域進(jìn)行離散。

圖1 有限體積空間離散單元(a)圓柱對稱網(wǎng)格單元(未做角度離散)； (b)三維圓柱網(wǎng)格單元(角度離散后)。r、θ分別表示半徑和角度。

本文采用后向歐拉法對時間進(jìn)行離散，其一階差分形式為

(10)

式中：g(t)表示當(dāng)前時刻數(shù)據(jù)；g(t-Δt)表示前一時刻數(shù)據(jù)，其中Δt是時間步長。本文在電流變化劇烈的階段采用較小的時間步長，在平穩(wěn)階段采用稍大的時間步長，這樣既能保證數(shù)據(jù)的穩(wěn)定性又能提高計算效率。結(jié)合式(9)，可通過下式求解第n+1個時間步長的場值

(11)

即

(12)

式中：時間步長Δt(n)=t(n+1)-t(n)，t(n+1)表示當(dāng)前時刻，t(n)表示前一時刻；I表示單位矩陣。

假設(shè)在0時刻向磁性源中通入非階躍電流，通過改變電流強(qiáng)度實現(xiàn)發(fā)射波形的瞬態(tài)響應(yīng)模擬。這里分別采用梯形波、半正弦波和三角波三種發(fā)射波形，波形示意圖見圖2，其激勵函數(shù)分別為

(13)

(14)

(15)

圖2 發(fā)射波形示意圖(a)梯形波； (b)半正弦波； (c)三角波

為確保時域麥克斯韋方程組解的唯一性及合理性，需給出初始條件和邊界條件。由于下階躍波在0時刻前保持穩(wěn)定的電流，會在整個空間產(chǎn)生穩(wěn)定的磁場，因此0時刻的初始條件為H(0)=H0，這里H0表示0時刻的磁場分布，0時刻的電場分布為E(0)=0。 但是，梯形波、半正弦波和三角波在0時刻之前并未接通電流，所以在整個空間任意位置的初始磁場和電場分布均為0，即H(0)=0，E(0)=0。本文使用Dirichlet邊界條件，假設(shè)求解區(qū)域外邊界Г的切向電場分量為0，即E×n=0|Г，這里n表示Г的外法向方向向量。最后采用直接求解器PARDISO對離散后的方程組進(jìn)行計算。

1.3 正確性驗證

本文程序運行平臺為Inter i5四核處理器，主頻為1.6GHz，內(nèi)存為12GB。為了驗證本文程序的正確性，建立圖3所示均勻半空間模型，對模型正演模擬得到的數(shù)值解與解析解進(jìn)行對比。觀測裝置收—發(fā)縱向偏移距為5m，線圈面積為1m2，發(fā)射電流為1A。模擬區(qū)域徑向均勻網(wǎng)格的最小尺度為0.3m，最大為3.0m，縱向均勻網(wǎng)格最小尺度為1.0m，最大為3.0m，在擴(kuò)展區(qū)域以1.3為等比因子進(jìn)行擴(kuò)充。在電流斷開1×10-6～1×10-2s內(nèi)選取41個時間點進(jìn)行觀測，早期時間步長為1×10-8s，隨時間的推移，逐漸增大時間步長。一次正演計算耗時175s。數(shù)值解和解析解的計算結(jié)果如圖4所示，縱軸dBz/dt表示垂直磁感應(yīng)分量Bz關(guān)于時間的

圖3 瞬變電磁均勻半空間模型

圖4 均勻半空間模型瞬變電磁響應(yīng)

導(dǎo)數(shù)。分析圖4可知，解析解和數(shù)值解吻合較好，說明本文瞬變電磁計算方法是正確的。

2 響應(yīng)特征分析

2.1 套管響應(yīng)特征分析

為了研究套管對瞬變電磁響應(yīng)的影響，分別討論套管電導(dǎo)率、相對磁導(dǎo)率、厚度和內(nèi)徑的響應(yīng)特征。建立圖5a所示套管測井模型，采用圖5b所示套管井瞬變電磁裝置，發(fā)射電流為5A。

圖5 套管井模型(a)及瞬變電磁裝置剖面示意圖(b)

r1、r2分別表示井眼內(nèi)徑和外徑； Δr=r2-r1表示套管厚度；rT、rR分別表示發(fā)射線圈和接收線圈的半徑，這里rT=rR=9.6cm；σh、σw、σt分別表示圍巖、井液、目標(biāo)層的電導(dǎo)率。

首先模擬不同套管電導(dǎo)率(σc)和相對磁導(dǎo)率μr(μr=μ/μ0，μ表示介質(zhì)磁導(dǎo)率，μ0=4π×10-7H/m表示真空中的磁導(dǎo)率)條件下的瞬變電磁響應(yīng)。μr=100，σc=1×105、5×105、1×106S/m的計算結(jié)果見圖6。由圖可知，當(dāng)其他參數(shù)不變時，σc越高，早期時間道的dBz/dt值越低，衰減越慢，衰減的時間也越長。σc=1×105S/m，μr=100、150、200的計算結(jié)果見圖7?？梢钥闯觯坠艿摩蘲越高，早期時間道上的響應(yīng)信號幅值越低，信號衰減越慢。

討論了套管電導(dǎo)率和相對磁導(dǎo)率對瞬變電磁探測的影響后，假設(shè)其他參數(shù)不變，分析套管厚度Δr及內(nèi)徑r1對瞬變電磁響應(yīng)信號的影響。

設(shè)套管電導(dǎo)率σc=1×105S/m，相對磁導(dǎo)率μr=100，套管厚度Δr分別為0.8、1.0、1.5cm的電磁響應(yīng)計算結(jié)果見圖8。分析圖8可知，套管厚度越大，早期瞬變電磁響應(yīng)信號越弱，且電磁響應(yīng)信號衰減越慢。

在其他參數(shù)不變的條件下，假設(shè)套管厚度Δr為1.0cm，內(nèi)徑r1分別為20.32、22.86、25.40cm時的響應(yīng)曲線見圖9?？梢钥闯?，所有響應(yīng)曲線的變化趨勢一致，但套管內(nèi)徑r1越大，整個觀測過程的信號幅值越低。

對瞬變電磁過套管響應(yīng)特征的分析結(jié)果表明，對套管井進(jìn)行瞬變電磁探測時，套管參數(shù)σc、μr、Δr、r1會對瞬態(tài)響應(yīng)結(jié)果產(chǎn)生較強(qiáng)的影響，研究其如何影響電磁場是校正金屬套管影響的基礎(chǔ)。

2.2 不同發(fā)射波形響應(yīng)特征分析

采用上節(jié)基本模型研究不同發(fā)射波形(梯形波、

圖6 套管不同電導(dǎo)率σc下瞬變電磁響應(yīng)曲線

圖7 套管不同相對磁導(dǎo)率μr下瞬變電磁響應(yīng)曲線

圖8 不同套管厚度Δr下瞬變電磁響應(yīng)曲線

半正弦波和三角波的波形)下的瞬變電磁套管井電磁響應(yīng)特征。假設(shè)套管電導(dǎo)率σc=1×105S/m，相對磁導(dǎo)率μr=100，發(fā)射波形通電時間為6ms，其中梯形波的上升沿(0～t1)和下降沿(t2～t3)時長均為0.01ms，分別計算目標(biāo)層電導(dǎo)率σt為0.20、0.01S/m時的電磁響應(yīng)，結(jié)果見圖10。

圖9 不同套管內(nèi)徑r1下瞬變電磁響應(yīng)曲線

圖10 不同σt下瞬變電磁響應(yīng)曲線(a)σt=0.20S/m； (b)σt=0.01S/m

對比圖10所示三種波形的響應(yīng)曲線可知，電磁響應(yīng)信號從高到低所對應(yīng)的發(fā)射波形為梯形波、半正弦波和三角波。不同發(fā)射波形的響應(yīng)信號特征受地下介質(zhì)電導(dǎo)率的影響顯著，且同種發(fā)射波形對高電導(dǎo)率目標(biāo)層的響應(yīng)信號強(qiáng)于低電導(dǎo)率目標(biāo)層的響應(yīng)信號。

2.3 梯形波各階段響應(yīng)特征分析

在實際工作中，由于很難實現(xiàn)對發(fā)射電流的即開即斷，所以分析梯形波各階段電磁響應(yīng)的差別對探測系統(tǒng)參數(shù)的設(shè)計具有重要意義。如圖2a所示，梯形波的供電過程包括三個階段：開啟階段(0～t1)、穩(wěn)定階段(t1～t2)和關(guān)斷階段(t2～t3)。開啟和關(guān)斷過程中電流均會在地層中激發(fā)渦流，產(chǎn)生瞬變電磁信號，這些信號可能會相互疊加或相互抵消。電流斷開后，利用觀測裝置測量地層中的感應(yīng)信號并對其進(jìn)行分析，可計算地層介質(zhì)的地電信息。

假設(shè)目標(biāo)層電導(dǎo)率σt為0.01S/m，開啟時間t0和關(guān)斷時長(t3-t2)均為0.1ms，其他參數(shù)不變，分別計算梯形波的穩(wěn)定時長(t2-t1)為1、2、6、20、50ms時的瞬變電磁響應(yīng)，結(jié)果見圖11。由圖可知，在波形開啟和關(guān)斷時長固定的條件下，當(dāng)波形穩(wěn)定時長小于6ms時，穩(wěn)定時長對響應(yīng)結(jié)果的影響明顯； 當(dāng)穩(wěn)定時長大于6ms時，各時間道上測得的電磁響應(yīng)結(jié)果基本不受影響。因此，在套管環(huán)境下開展瞬變電磁探測需設(shè)置合理的梯形波穩(wěn)定時長，一方面要考慮探測效果，另一方面，穩(wěn)定時長過長則需要更多的時間步長，會增加計算時間，在野外工作中也會增加能耗。

假設(shè)梯形波的關(guān)斷時長(t3-t2)為0.10ms，開啟和穩(wěn)定時長之和t2為6ms，其他參數(shù)不變，分別計

圖11 梯形波不同穩(wěn)定時長的瞬變電磁響應(yīng)曲線

算開啟時間t1為0.01、0.10、1.00、2.00ms的電磁響應(yīng)，結(jié)果見圖12。當(dāng)然，波形穩(wěn)定時長會隨開啟時間的增加而減少。由圖12可見，在波形開啟和穩(wěn)定時長總和不變的條件下，波形開啟時間對電磁響應(yīng)結(jié)果的影響差異很小，因此，在電流關(guān)斷后進(jìn)行觀測時，可以忽略開啟時間對電磁響應(yīng)的影響。

發(fā)射波形的關(guān)斷時間(t3-t2)是判斷電磁探測設(shè)備性能的重要指標(biāo)之一，因此，分析關(guān)斷時間對電磁信號的影響非常重要。假設(shè)波形開啟和穩(wěn)定時間總和t2為6ms，開啟時間t1為0.1ms，分別計算關(guān)斷時間(t3-t2)為0.001、0.010、0.100、0.200ms的電磁響應(yīng)，結(jié)果見圖13。分析圖中衰減曲線可知：電流關(guān)斷時間為0.010、0.100、0.200ms時，關(guān)斷時間越短，早期時段觀測信號越強(qiáng)，晚期時段觀測信號曲線基本重合； 關(guān)斷時間為0.001ms與0.010ms的響應(yīng)曲線幾乎完全重合?？梢?，關(guān)斷時間對響應(yīng)信號的影響主要集中在早期，對晚期的影響很小，而早期電磁信號反映的是淺層信息，因此開展瞬變電磁勘探時，設(shè)計合理的梯形波關(guān)斷時間尤為重要。

圖12 發(fā)射波形不同開啟時間的瞬變電磁響應(yīng)結(jié)果

圖13 梯形波不同關(guān)斷時間的瞬變電磁響應(yīng)曲線

3 結(jié)論

本文基于圓柱坐標(biāo)系開展套管井瞬變電磁建模，運用有限體積法和后向歐拉法分別進(jìn)行空間和時間離散，最終通過求解時域麥克斯韋方程組得到三維瞬變電磁過套管電磁響應(yīng)。分別針對套管的電導(dǎo)率、相對磁導(dǎo)率、厚度和內(nèi)徑等參數(shù)對電磁場響應(yīng)的影響特征進(jìn)行分析，為瞬變電磁過套管數(shù)據(jù)的處理及探測裝置參數(shù)的設(shè)計提供了技術(shù)支撐，得到如下結(jié)論。

(1)套管的電導(dǎo)率越高或相對磁導(dǎo)率越大，三維瞬變電磁過套管響應(yīng)信號的衰減時間越長，早期的響應(yīng)幅值越低； 套管的內(nèi)徑或厚度越大，早期的響應(yīng)信號幅值越低。

(2)發(fā)射波采用梯形波形、半正弦波形和三角波形，所對應(yīng)的套管井瞬變電磁響應(yīng)信號逐次減弱。

(3)發(fā)射波采用梯形波形且波形開啟和關(guān)斷時間固定時，隨穩(wěn)定時間的增加，早期時間道上的電磁響應(yīng)信號減弱，當(dāng)穩(wěn)定時間足夠長時，穩(wěn)定時間的變化對電磁響應(yīng)的影響逐漸減弱； 波形開啟和穩(wěn)定時間之和固定時，波形開啟時間的長短對電磁響應(yīng)的影響不大，可以忽略； 發(fā)射波形的關(guān)斷時間對電磁響應(yīng)的影響主要集中在早期，對晚期電磁響應(yīng)的影響較小。