井間電磁成像測井?dāng)?shù)值模式匹配算法

臧德福, 姬勇力, 李智強(qiáng), 晁永勝, 楊愛鋒

(1.中石化勝利石油工程有限公司測井公司, 山東 東營 257096;2.中國電子科技集團(tuán)公司第22研究所, 河南 新鄉(xiāng) 453003)

0 引 言

井間電磁成像技術(shù)主要用來對(duì)井間的電導(dǎo)率分布進(jìn)行成像。與電測井相比,它具有更大的探測范圍;與大地電磁及地面-井電磁法相比,它具有更高的精度[1]。利用該方法可進(jìn)一步求取孔隙度、流體飽和度分布和滲透率分布信息,從而發(fā)現(xiàn)漏失的油氣層,監(jiān)測油氣水的運(yùn)移情況,確定最佳的布井位置,提高采收率。

20世紀(jì)80年代以來,美國加州大學(xué)伯克利分校的研究人員對(duì)井間電磁成像技術(shù)進(jìn)行了大量的理論和實(shí)驗(yàn)研究[2-6],并取得了重要進(jìn)展,包括井間電磁系統(tǒng)的靈敏度分析和適用范圍、井間電磁方法的分辨率、儀器裝備研究、金屬套管的影響、井間電磁成像方法。自1997年勝利油田、中國石油大學(xué)(華東)與美國的EMI(Electromagnetic Instruments Inc.)公司合作,開展井間電磁成像系統(tǒng)的應(yīng)用研究,并在勝利油田進(jìn)行了3對(duì)井的測量試驗(yàn)。3次試驗(yàn)都采集到可靠的井間電磁數(shù)據(jù),利用這些數(shù)據(jù)進(jìn)行軸對(duì)稱二維近似反演得到的井間電導(dǎo)率成像圖在分析井間油氣分布、砂體展布等方面見到較好的效果,充分證明井間電磁成像技術(shù)是油藏研究中具有開創(chuàng)性的重要技術(shù)手段[7-8]。近年來,勝利測井公司又開展了更遠(yuǎn)距離、更高性能的井間電磁成像技術(shù)的研究[9]。

計(jì)算井間電磁場的正演方法包括微分方程法和積分方程法2類。求解微分方程可采用有限元法和有限差分法,其求解范圍是整個(gè)井間區(qū)域,經(jīng)離散后形成含大型稀疏矩陣的線性方程組,適用于復(fù)雜井間地層模型。由于井間電磁場的求解區(qū)域要比單井電磁場大得多,導(dǎo)致計(jì)算量巨大,只有超級(jí)計(jì)算機(jī)才能完成。積分方程的求解范圍只是井間電導(dǎo)率異常區(qū)域,因而與微分方程相比求解區(qū)域要小得多,適用于簡單的井間地層模型。積分方程法更適合于反演過程中的正演計(jì)算,因?yàn)槠溆?jì)算量相對(duì)較小,且成像區(qū)域可選為有限區(qū)域。但是,積分方程經(jīng)離散后形成含致密矩陣的線性方程組,若在計(jì)算時(shí)存儲(chǔ)整個(gè)矩陣,則所需內(nèi)存量很大,計(jì)算速度將變慢,從而降低計(jì)算效率。上述2種方法計(jì)算效率都不高,因此,開發(fā)出高效的井間電磁成像的正演模擬方法是該領(lǐng)域的一項(xiàng)重要研究內(nèi)容[10]。

Padensi等[11-12]在研究電磁散射時(shí),把波膜的概念與有限元結(jié)合起來,提出一種稱為數(shù)值模式匹配(NMM)的半解析、半有限元解法。該方法把二維數(shù)值問題轉(zhuǎn)化為一維解析解和一維數(shù)值解的結(jié)合,大大減少了計(jì)算量。Chew等[13-14]把該方法應(yīng)用到普通電阻率測井中,其效率是有限元的數(shù)百倍。張庚驥等[15-16]用NMM法研究了軸對(duì)稱條件下縱向成層、徑向不均勻的地層模型中普通電阻率測井的響應(yīng),在縱向上采用解析方法推導(dǎo)了反射矩陣和透射矩陣的遞推公式,在徑向上采用數(shù)值解時(shí)改進(jìn)了基函數(shù)(幅度基函數(shù)和斜度基函數(shù)),其結(jié)果不僅達(dá)到了有限元方法的精度,而且大大提高了運(yùn)算速度。仵杰等[17]針對(duì)NMM算法在復(fù)雜測井情況下的誤差詳細(xì)分析了基函數(shù)類型、個(gè)數(shù)以及徑向網(wǎng)格的剖分邊界位置,提出了相應(yīng)的改進(jìn)措施。本文在上述基礎(chǔ)上,將NMM法用于井間電磁成像正演計(jì)算,結(jié)合井間電磁成像數(shù)據(jù)測量的特點(diǎn),利用并行計(jì)算快速計(jì)算出所有數(shù)據(jù)。

1 井間電磁成像中的NMM法

1.1 地層模型和NMM法原理

井間電磁成像是在2口井中分別放入發(fā)射與接收裝置,通過測量不同深度的磁場響應(yīng)反演2口井之間剖面的電導(dǎo)率分布(見圖1)。

圖1 井間電磁成像地層模型

(1)

在無發(fā)射源時(shí),式(1)等式右邊為0,即

(2)

其邊界條件為

1.2 NMM法基本原理

采用分離變量法對(duì)式(2)求解[15-16],令

(3)

而fm(r)沒有解析解,只能用數(shù)值解求解。fm(r)近似表示為

(4)

式中,Cm是本征向量cm組成的矩陣。式(4)等價(jià)于式(5)矩陣形式的本征值方程

(5)

1.3 NMM法計(jì)算步驟

(1) 建立地層模型。

(2) 選取基函數(shù),形成矩陣Am、B。

(3) 求解廣義特征值問題,得到Cm。

(4) 確定發(fā)射源所在位置。

(7) 確定接收線圈的位置,代入式(3)即可求出接收線圈處的電場強(qiáng)度。

2 井間電磁成像中NMM法的應(yīng)用

2.1 重構(gòu)NMM法計(jì)算步驟

從上面的推導(dǎo)可知,NMM法是1個(gè)發(fā)射點(diǎn)計(jì)算1個(gè)接收?qǐng)鰪?qiáng),井間電磁成像系統(tǒng)為確保成像效果,需盡可能多地采集數(shù)據(jù)。通常做法是接收線圈固定在某深度,發(fā)射線圈在整個(gè)測量深度區(qū)間和移動(dòng)以固定間距采集數(shù)據(jù),然后變換接收線圈深度,重復(fù)測量。例如井間距500 m,深度區(qū)間200 m,每2 m采集1個(gè)點(diǎn),則共需采集10 000個(gè)點(diǎn),逐個(gè)計(jì)算場強(qiáng)會(huì)導(dǎo)致一次正演所用時(shí)間非常長。分析NMM法的計(jì)算原理與井間電磁成像的測量方式,可知當(dāng)?shù)刭|(zhì)模型建立以后,正演中含有大量重復(fù)計(jì)算,可以省略,這樣大大縮短正演時(shí)間。

(2) 地層劃分要足夠細(xì),反演成像分辨率才高,Cm與Λm基本上每層模型都要用到,可以將2個(gè)參數(shù)在第1次計(jì)算時(shí)存儲(chǔ)下來,在后面的每個(gè)發(fā)射-接收點(diǎn)直接調(diào)用。只要地層電導(dǎo)率模型不變,Cm與Λm也不變。

(3) 每個(gè)接收點(diǎn)都對(duì)應(yīng)相同數(shù)量的發(fā)射點(diǎn)位置,可以將僅與發(fā)射點(diǎn)相關(guān)的公式全部計(jì)算出來,這樣就避免了每個(gè)接收點(diǎn)都重復(fù)計(jì)算。

從上述計(jì)算步驟中可以看出,在前3步中可提前計(jì)算出固定值用于后續(xù)計(jì)算,這樣節(jié)省大量運(yùn)算時(shí)間,而且其所占存儲(chǔ)空間相對(duì)于有限元與有限差分也非常小,普通計(jì)算機(jī)完全可以滿足需求。

2.2 數(shù)值模式匹配的并行計(jì)算

從20世紀(jì)60年代起,人們就開始探索數(shù)值計(jì)算方法的并行化技術(shù)和計(jì)算結(jié)構(gòu)的并行化設(shè)計(jì)方法,解決單核計(jì)算機(jī)難以完成的工程與科學(xué)技術(shù)領(lǐng)域內(nèi)規(guī)模巨大、實(shí)現(xiàn)要求嚴(yán)格的數(shù)值計(jì)算問題。并行計(jì)算就是把一個(gè)需要非常巨大的計(jì)算能力才能解決的問題分成許多小問題,把這些小問題分配給不同計(jì)算機(jī)或處理器進(jìn)行并行處理,將計(jì)算結(jié)果綜合得到最終的結(jié)果。為加快運(yùn)算速度和解決大容量存儲(chǔ)的求解問題,并行編程需要具有多個(gè)內(nèi)部處理器的單計(jì)算機(jī)或者多個(gè)互聯(lián)的計(jì)算機(jī)[18]。以1999年NVIDIA公司提出GPU概念為分界點(diǎn),之前的并行計(jì)算都是基于多CPU并行處理,之后的并行計(jì)算向著多核CPU與GPU異構(gòu)協(xié)同的方向發(fā)展。CPU并行與GPU并行的區(qū)別在于單個(gè)GPU的多核心,重復(fù)計(jì)算能力強(qiáng),通過低投入的GPU計(jì)算陣列就可以達(dá)到以往大型CPU陣列并行系統(tǒng)的效率。但是GPU的每個(gè)處理核心在執(zhí)行指令的時(shí)候控制優(yōu)化少,速度不夠高,近些年出現(xiàn)的多核CPU其核心最多也就16核,并行度不高,因此未來的發(fā)展方向是CPU/GPU異構(gòu)協(xié)同發(fā)展[19]。

GPU并行大多用于超級(jí)計(jì)算機(jī),考慮到井間電磁成像系統(tǒng)數(shù)據(jù)處理的實(shí)際情況以及便攜性與經(jīng)濟(jì)性,采用多核CPU并行基本可以滿足需要,而當(dāng)前無論是臺(tái)式機(jī)還是便攜機(jī),都具有多核CPU,若采用并行計(jì)算將多核CPU充分利用,則重構(gòu)NMM法的計(jì)算速度可以大幅提高。圖2是重構(gòu)NMM法的并行處理流程圖。

圖2 NMM法在井間電磁成像正演中流程圖

2.3 模擬仿真

2.3.1 均勻地層對(duì)比驗(yàn)證

通過與均勻地層的解析解對(duì)比,驗(yàn)證NMM法的精確度。設(shè)均勻介質(zhì)中有一半徑為aT、匝數(shù)為NT的發(fā)射線圈,其中通過交變電流IT=I0ej ω t。它在空間形成的交變電場與磁場也隨時(shí)間按照正弦規(guī)律變化。將坐標(biāo)原點(diǎn)選在線圈中心,z軸垂直于線圈平面,在軸對(duì)稱的柱坐標(biāo)系(r,φ,z)中,有

(6)

建立地質(zhì)模型,縱向100層,徑向300層,網(wǎng)格間距2 m×2 m,地層電導(dǎo)率為1 S/m,幅度和斜度基函數(shù)分別取52點(diǎn)。圖3為10 Hz頻率下,隨井間距變化的誤差示意圖。NMM算法在頻率選取原則范圍內(nèi)的精度完全滿足系統(tǒng)需求[7]。

圖3 NMM法與解析解誤差隨井間距變化關(guān)系

為了說明本文算法的優(yōu)勢,在四核便攜式計(jì)算機(jī)(主頻2.5 GHz)上采用原NMM法及結(jié)合井間電磁測量的NMM法、并行NMM法計(jì)算,對(duì)其運(yùn)算時(shí)間進(jìn)行對(duì)比,在XP雙核(主頻2.5 GHz)和win7 四核系統(tǒng)的臺(tái)式機(jī)(主頻3.3 GHz)進(jìn)行并行運(yùn)算對(duì)比。圖4是運(yùn)算耗費(fèi)時(shí)間示意圖,從圖4中可知,原NMM方法一個(gè)點(diǎn)需要36 s,而且多個(gè)點(diǎn)重復(fù)調(diào)用單點(diǎn)計(jì)算程序,其所用時(shí)間與計(jì)算點(diǎn)數(shù)成線性正比關(guān)系。結(jié)合井間電磁成像工作原理改進(jìn)后的NMM法單點(diǎn)計(jì)算只需10 s,隨著計(jì)算點(diǎn)數(shù)的增加,優(yōu)化后的NMM法計(jì)算速度優(yōu)勢越來越明顯,對(duì)于1萬點(diǎn)計(jì)算,原方法需100 h,而優(yōu)化后的方法只需12 min。當(dāng)采用雙核并行時(shí),僅需7 min,若在主頻更高計(jì)算機(jī)上進(jìn)行4核并行計(jì)算,只需2.5 min,可以大大提高NMM法反演成像的效率。

圖4 NMM法計(jì)算速度對(duì)比

2.3.2 縱向成層一維井間地層模型

通過與均勻地層的對(duì)比證明了該方法的正確性以及高效率,然而均勻地質(zhì)模型很難見到,因此,為進(jìn)一步證明該方法的有效,按表1建立縱向成層一維井間模型,該模型有解析解[20]。地層共3層,邊界分別為0、50 m和100 m,每層電導(dǎo)率見表1。模型邊界以外電導(dǎo)率設(shè)為1 S/m,設(shè)頻率為10 Hz,計(jì)算接收?qǐng)鰪?qiáng)隨井間距變化的誤差(見圖5)。計(jì)算結(jié)果表明,對(duì)于縱向成層一維井間地層,NMM算法的計(jì)算誤差也非常小。

表1 水平成層一維井間地層電阻率

圖5 NMM法計(jì)算速度對(duì)比

2.3.3 二維軸對(duì)稱地層模型

目前,井間電磁三維成像難度太大,不易實(shí)現(xiàn)。對(duì)井間電磁進(jìn)行二維或2.5維成像是研究最為廣泛且成果最多的,魏寶君等[1]利用積分方程法針對(duì)軸對(duì)稱二維非均勻地層的井間電磁成像問題做了大量的正演、反演工作。本文利用NMM與積分方程法在相同軸對(duì)稱二維非均勻模型下開展正演計(jì)算,比較兩者的計(jì)算精度以驗(yàn)證NMM方法用于井間電磁成像的可行性。

圖6 軸對(duì)稱二維非均勻地層模型

首先建立軸對(duì)稱二維非均勻地質(zhì)模型(見圖6),數(shù)值代表電導(dǎo)率(單位S/m)。地層分10層,每層厚度為10 m,徑向500 m,縱向?yàn)?00～400 m,即測量深度范圍為100 m。第1層起始邊界300 m,第2層起始邊界310 m,以此類推,第10層起始邊界390 m。發(fā)射井以及模型邊界外的電導(dǎo)率設(shè)為1 S/m。不考慮接收井,只需計(jì)算接收線圈在距發(fā)射線圈徑向50、100、300、500 m的磁場強(qiáng)度。以發(fā)射線圈中心為徑向零點(diǎn),發(fā)射井半徑為0.101 6 m,其他各區(qū)塊對(duì)應(yīng)邊界如表2所示。發(fā)射點(diǎn)、接收點(diǎn)均為深度350 m處,頻率為10 Hz時(shí),計(jì)算磁場強(qiáng)度隨井間距變化的結(jié)果,利用積分方程法(IE)、有限元法(FEM)與NMM這3種方法分別計(jì)算,其結(jié)果如圖7所示。從圖7可知,3種方法的計(jì)算結(jié)果基本一致,且NMM與FEM方法隨著井間距的增加,其精度要高于IE方法。

表2 水平成層一維井間地層電阻率

圖7 NMM與積分方程法、有限元法井間電磁成像正演計(jì)算結(jié)果對(duì)比

3 結(jié) 論

(1) NMM法可在常規(guī)配置的計(jì)算機(jī)上進(jìn)行較大數(shù)據(jù)量的處理,適用于井間電磁成像這些復(fù)雜的正演數(shù)據(jù)處理。

(2) 結(jié)合井間電磁成像系統(tǒng)工作原理,通過重構(gòu)獲得改進(jìn)NMM法。在改進(jìn)NMM法的基礎(chǔ)上采用并行優(yōu)化技術(shù),最終獲得并行NMM法。

(3) 通過與均勻地層、水平成層一維地層解析解對(duì)比,以及與軸對(duì)稱非均勻地層積分方程法、有限元法對(duì)比,驗(yàn)證了改進(jìn)并行NMM法應(yīng)用于大范圍的井間電磁成像正演計(jì)算是可行的,而且效率更高。

參考文獻(xiàn):

[1] 魏寶君, 張庚驥, 曾文沖. 井間電磁成像的迭代反演算法 [J]. 地球物理學(xué)報(bào), 1999, 42(5): 711-719.

[2] Wm R Petrick. Three Dimensional Resistivity Inversion Using Alpha Centers [J]. Geophysics, 1981, 46(8): 1148-1162.

[3] Tripp A C. Two-dimensional Resistivity Inversion [J]. Geophysics, 1984, 49(10): 1708-1717.

[4] Zhang Z, Zhou Z. Real-time Quasi-2-D Inversion of Array Resistivity Logging Data Neural Networks [J]. Geophysics, 2002, 67(2): 517-524.

[5] Xiong Zong-hou. Electromagnetic Modeling of 3-D Structures the Method of System Iteration Using Integal Equations [J]. Geophysics, 1992, 57(12): 1556-1561.

[6] Wannamaker P E, Hohmann G W, Sanlilipo WA. Electromagnetic Modeling of Three-dimensional Bodies in Layered Earths Using Integral Equations [J]. Geophysics, 1984, 49(1): 60-74.

[7] 曾文沖, 趙文杰, 臧德福. 井間電磁成像系統(tǒng)應(yīng)用研究 [J]. 地球物理學(xué)報(bào), 2001, 44(3): 411-420.

[8] 郭紅旗, 臧德福, 王群力, 等. XBH2000井間電磁成像測井系統(tǒng) [J]. 石油儀器, 2002, 16(3): 25-27.

[9] 臧德福, 郭紅旗, 晁永勝, 等. 井間電磁成像測井系統(tǒng)分析與研究 [J]. 測井技術(shù), 2013, 37(2): 177-182.

[10] 栗建軍. 井間電磁測井原理、方法及套管對(duì)井間電磁測井影響規(guī)律的研究 [D]. 青島: 中國海洋大學(xué), 2004.

[11] Padensi M A, Ferreira L G. Method to Calculate the Reflection and Transmission of Guided Waves [J]. J. Opt. Soc. An. , 1982, 72(1): 126-130.

[12] Chew W C, et al. Diffraction of Waves by Discontinuities in Open Cylindrical Structures [M]. Houston: Proc. Internet. Symp. Institute of Elect. and Electron. Eng. Antennas and Propagation, 1983: 503-506.

[13] Chew W C, Barone S, Anderson B, et al. Diffrection of Axisymmetric Waves in a Borehole by Bed Boundary Discontinuities [J]. Geophysics, 1984, 49(24): 1586-1595.

[14] 聶在平, Chew W C, Liu Q H. 電磁波對(duì)軸對(duì)稱二維層狀介質(zhì)的散射 [J]. 地球物理學(xué)報(bào), 1992, 35(4): 479-489.

[15] 張庚驥, 汪涵明, 汪功禮. 成層介質(zhì)中交流電測井響應(yīng) [J]. 地球物理學(xué)報(bào), 1995, 38(6): 840-849.

[16] 張庚驥, 汪涵明. 普通電阻率測井的數(shù)值模式匹配算法 [J]. 石油大學(xué)學(xué)報(bào): 自然科學(xué)版, 1996, 20(2): 23-29.

[17] 仵杰. 感應(yīng)測井中的NMM法及其在咸水泥漿井中應(yīng)用 [J]. 國外測井技術(shù), 2011(1): 15-19.

[18] 李鵬, 邵明剛. 并行計(jì)算技術(shù) [J]. 中國科技信息, 2006(7): 24-25.

[19] 盧風(fēng)順, 宋君強(qiáng), 銀?？? 等. CPU/GPU協(xié)同并行計(jì)算研究綜述 [J]. 計(jì)算機(jī)科學(xué), 2011, 38(3): 5-9.

[20] 田子立, 孫以睿, 劉桂蘭. 感應(yīng)測井理論及其應(yīng)用 [M]. 北京: 石油工業(yè)出版社, 1984.