方思南， 潘和平*， 杜婷， 王智， 鄧呈祥

1 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢　400074 2 湖北煤炭地質(zhì)勘查院，武漢　430070

考慮卷積完全匹配層數(shù)值色散的井間電磁三維正演

方思南1， 潘和平1*， 杜婷2， 王智1， 鄧呈祥1

1 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢400074 2 湖北煤炭地質(zhì)勘查院，武漢430070

摘要本文將以卷積完全匹配層為吸收邊界條件的時域有限差分法應(yīng)用到井間電磁的三維正演模擬中. 證明了卷積完全匹配層中的數(shù)值色散會因?yàn)橛行а由煲蜃佣a(chǎn)生，列舉常規(guī)有效延伸因子和網(wǎng)格間距對電磁波相速度各向異性的影響，并通過波場快照驗(yàn)證卷積完全匹配層中數(shù)值色散的存在；進(jìn)而推廣三維卷積完全匹配層中最大有效延伸因子、最大網(wǎng)格間距與激勵源主頻之間的約束，以此完善卷積完全匹配層的最優(yōu)參數(shù)設(shè)置方案. 在此基礎(chǔ)上，通過正演結(jié)果中二次場的垂直磁場分量和走時來展示靈敏度分布，以此劃定井間電磁勘探的優(yōu)勢區(qū)域，提出井間電磁正反演所需要的約束測井系列和最佳取井方案.

關(guān)鍵詞井間電磁(CWEM)； 三維正演； 時域有限差分(FDTD)； 卷積完全匹配層(CPML)； 數(shù)值色散

1引言

隨著現(xiàn)今油氣新區(qū)的勘探難度加大，油氣勘探方向越來越多的聚焦在已經(jīng)進(jìn)入到開發(fā)中后階段的油田上，這就需要更加準(zhǔn)確地獲取井間的地層構(gòu)造、儲層分布和油氣運(yùn)移等信息.井間電磁(Cross-Well Electromagnetic，CWEM)法，因其能低成本、高效率地完成數(shù)據(jù)采集工作，能通過反演電阻率來定量評價儲層的含油性，能將高豐度的一維測井信息延展到三維儲層空間中(沈金松等，2014)，能通過與地震數(shù)據(jù)聯(lián)合反演來監(jiān)控儲層動態(tài)信息(Carcione et al.，2012)，所以可以勝任油田二次開發(fā)的勘探要求.

石油勘探中井間電磁成像系統(tǒng)的發(fā)展，由20世紀(jì)90年代的先導(dǎo)性理論研究(Newman，1995；Wilt et al.，1995)，到EMI公司生產(chǎn)的XBH-2000系統(tǒng)(曾文沖等，2001)，再到近年Schlumberger公司改進(jìn)的Deeplook-EM系統(tǒng)(Al-Ali et al.，2009)，其中，Deeplook-EM井間電磁成像系統(tǒng)采用主要發(fā)射頻率為5～1000 Hz的磁偶極子發(fā)射線圈，和靈敏度達(dá)到10-6nT的陣列接收線圈.在井間距確定的情況下，井間電磁的最大發(fā)射頻率主要受套管、地層電阻率、發(fā)射線圈功率和接收線圈靈敏度的限制.

現(xiàn)階段關(guān)于井間電磁的數(shù)值模擬研究主要集中在金屬套管的影響規(guī)律(魏寶君等，2014)和理論模型的正演響應(yīng)(Maclennan et al.，2014；Donadille and Al-Ofi，2012)方面，本文主要關(guān)注的是時域有限差分(Finite-Difference Time-Domain，F(xiàn)DTD)法在三維井間電磁波正演中的應(yīng)用.

由于目前井間電磁勘探的最大井間距約900 m(沈金松等，2014)，縱向分辨率受發(fā)射頻率、地層電導(dǎo)率的影響，通常為井間距的2%～5%(Alumbaugh and Morrison，1995)，而常見儲層的厚度通常在10 m以內(nèi)，加上井周3 m范圍內(nèi)的地層電阻率受泥漿侵入影響嚴(yán)重，故針對井間電磁正演的時域有限差分需要在小區(qū)域、小網(wǎng)格中進(jìn)行，而這與低頻電磁波中較大的趨膚深度相矛盾.因此，低頻井間電磁波正演的關(guān)鍵問題是使吸收邊界在靠近目標(biāo)體的同時還能通過小網(wǎng)格來吸收低頻電磁波.

在常規(guī)PML(Perfectly Matched Layers)的基礎(chǔ)上(Bérenger，1994，1996；Chew and Weedon，1994)，具有嚴(yán)格因果關(guān)系、對低頻波有較好吸收效果的復(fù)頻移完全匹配層(Complex Frequency-Shifted Perfectly Matched Layers，CFS-PML)被提出來(Kuzuoglu and Mittra，1996)，CFS-PML可以被設(shè)置到緊鄰目標(biāo)體的位置以減少網(wǎng)格、節(jié)省內(nèi)存，可以保證迭代后期的反射較小(Bérenger，2012)，可以對某一頻率范圍內(nèi)的電磁波做最強(qiáng)吸收，因而被成功地應(yīng)用到低頻電磁波的正演中(李展輝和黃清華，2014).本文采用的卷積完全匹配層(Convolutional Perfectly Matched Layers，CPML)算法(Roden and Gedney，2000)，是CFS-PML的一種成熟有效的執(zhí)行方法，具有非常廣闊的發(fā)展前景(Bérenger，2007).

目前關(guān)于CPML優(yōu)勢吸收頻率的研究比較深入(Bérenger，2012)，而關(guān)于CPML網(wǎng)格數(shù)值色散的研究剛剛起步(de la Kethulle de Ryhove and Mittet，2014)，且主要對一維網(wǎng)格間距要求做的理論推導(dǎo)，并未考慮CPML所在的三維坐標(biāo)延伸空間，本文重點(diǎn)分析三維CPML中數(shù)值色散的來源、現(xiàn)象和限制方案.

2CPML的最優(yōu)參數(shù)

CPML算法基于坐標(biāo)延伸的麥克斯韋方程組，通過平面波在任意兩層介質(zhì)分界面上的零反射條件來提供本構(gòu)參數(shù)的約束關(guān)系.

2.1含CPML的FDTD時域方程

頻率域含 CFS-PML的無源麥克斯韋旋度方程(Bérenger，2007)中，常規(guī)三維歐幾里得空間的拉普拉斯算子因?yàn)镃FS-PML中的坐標(biāo)延伸因子而被改寫為一般形式(Chew and Weedon，1994)，如e分量(電場分量)v方向(v=x,y,z)的坐標(biāo)延伸因子sev可寫為(Kuzuoglu and Mittra，1996)：

(1)

其中，αev是電場在v方向的坐標(biāo)延伸因子的自由度，σev是CFS-PML中的電導(dǎo)率，κev是有效延伸因子(Bérenger，2007)，ε是背景場的介電常數(shù).

另外，麥克斯韋方程變換到時域中會因?yàn)樽鴺?biāo)延伸因子而產(chǎn)生卷積項(xiàng)，如式(2)為時域變換后正演區(qū)域和吸收邊界通用的Ex離散形式(Elsherbeni and Demir，2009)：

(2)

(3)

其中，

(4)

(5)

CPML中Ex的更新方程(式(2)—式(5))可以進(jìn)一步推廣到Ey、Ez、Hx、Hy、Hz，在常規(guī)FDTD(Elsherbeni and Demir，2009)的基礎(chǔ)上，CPML對電磁波的吸收效果主要受αuv、σuv(Bérenger，2012)和κuv(u=e,m，v=x,y,z)的影響.結(jié)合前人對于αuv、σuv的深入研究(李展輝和黃清華，2014；Bérenger，2012)，本文重點(diǎn)關(guān)注CPML中κuv的數(shù)值色散問題，并討論最優(yōu)參數(shù)的選取方案.

2.2α影響CPML的吸收頻率

參考式(1)，可以設(shè)置CPML的特征頻率fav為(Bérenger，2012)：

(6)

通常在CPML算法中，設(shè)置平行于CPML法方向上的參數(shù)αev、αmv由CPML內(nèi)邊界的αmin線性變化到CPML外邊界的αmax，以保證某一頻率范圍內(nèi)的入射波能夠被均勻吸收(Elsherbeni and Demir, 2009)：

(7)

其中，l為從CPML內(nèi)邊界開始數(shù)的網(wǎng)格層數(shù)，L為CPML從內(nèi)邊界到外邊界的總層數(shù).那么，對于α的設(shè)置，應(yīng)先結(jié)合式(6)中入射波的頻率范圍，來確定αmin和αmax，再通過式(7)來計算CPML中各層的αev和αmv.

2.3σ影響CPML的吸收效率

電導(dǎo)率σev和導(dǎo)磁率σmv直接影響CPML對電磁波的吸收效率：σev、σmv過小會降低CPML的吸收效率，過大會在CPML的層界面上產(chǎn)生顯著的反射(Bérenger，2012).在設(shè)置平行CPML法方向的σev和σmv時，采用式(8)來賦值最能提升CPML的吸收效率：

(8)

其中，冪m可以取2、3、4；σeb為背景場電導(dǎo)率，若背景場為真空環(huán)境，則令σeb=αmin；最外層CPML電導(dǎo)率σmax=ξσeb，可以取比例系數(shù)ξ<10.

2.4κ影響CPML的數(shù)值色散

與關(guān)注度較高的α和σ不同的是，目前CPML中κ的影響還沒有學(xué)者重點(diǎn)關(guān)注過，直觀的理解是，κ造成CPML的網(wǎng)格坐標(biāo)延伸，那么與電導(dǎo)率的作用一致，κ越大吸收越強(qiáng)烈.但在正演實(shí)驗(yàn)中，我們發(fā)現(xiàn)當(dāng)κ大于某一范圍后，外側(cè)CPML會產(chǎn)生持續(xù)的、震蕩的干擾波，這是典型的數(shù)值色散特征.根據(jù)這一發(fā)現(xiàn)，本節(jié)推導(dǎo)由κ產(chǎn)生的數(shù)值色散的原因，并給出差分近似后相速度的各向異性關(guān)系，以期研究控制數(shù)值色散的方法.

平行于CPML法方向上的κuv由式(9)給出(Elsherbeni and Demir, 2009)：

(9)

其中，κmax為CPML最外層的值，通?？梢匀?～11(Roden and Gedney，2000).

綜合式(7)、式(8)和式(9)可以發(fā)現(xiàn)，CPML中的參數(shù)都是以垂直于邊界的坐標(biāo)為軸對稱分布，因此對于三維CPML而言，我們可以簡化考慮二維情況下的數(shù)值色散問題.這里以垂直于y坐標(biāo)的CPML為研究對象，采用TE波從垂直于z坐標(biāo)的各個角度入射，此時所有電場對于z坐標(biāo)而言都是橫向電場，即構(gòu)成TEz模式，不難推導(dǎo)TEz模式在自由空間中的平面波解：

-ω(n+0.5)Δt).

(10)

再由式(2)推導(dǎo)出CPML中TEz模式的二維離散麥克斯韋方程組，如式(11)所示(由于在垂直于y坐標(biāo)的CPML中，κex和κmx的值為1，Ψeyx、Ψmzx的值為0，所以這里直接給出簡化后的結(jié)果)：

(11)

將TEz模式的平面波解式(10)代入式(11)中，結(jié)合三角恒等式：

cos(p-q)-cos(p+q)=2sinpsinq,

(12)

(13)

正是因?yàn)棣实挠绊?，造成CPML中的二維數(shù)值色散關(guān)系與常規(guī)FDTD的不同.為了更直觀地表現(xiàn)這種數(shù)值色散現(xiàn)象，進(jìn)一步令常規(guī)網(wǎng)格均勻剖分Δx=Δy=δ，波矢量滿足：

(14)

(15)

(16)

CPML數(shù)值色散關(guān)系的最終表達(dá)式(式(16))說明TEz模式平面波的相速度vφ與θ、λ/δ、κ有關(guān)，將此關(guān)系繪制成二維情況下的極坐標(biāo)圖(圖1)，若圖中θ=0°方向代表的是y方向，θ=90°方向代表的是x方向，那么圖1可以理解為某一層的CPML介質(zhì)垂直于y方向平鋪，平面波可以從圓弧上的任意方向入射.圖1a限定λ/δ=100，以κ為參變量，給出相速度與光速之比vφ/c與入射角θ之間的關(guān)系；圖1b限定κ=2和κ=5，以λ/δ為參變量，給出vφ/c與θ之間的關(guān)系.

圖1a中設(shè)置的λ/δ=100條件，保證了來自均勻網(wǎng)格的FDTD的數(shù)值色散非常小，在θ=90°方向上，vφ/c=1，而隨著κ的增大和θ的減小，數(shù)值色散不斷加劇，且在θ=0°方向近似呈vφ/c=1/κ的關(guān)系.圖1b主要反映λ/δ對數(shù)值色散的影響：引起常規(guī)網(wǎng)格數(shù)值色散的因素在CPML中也適用，若λ/δ<10則會因?yàn)榫W(wǎng)格間距過大在各個方向產(chǎn)生不可忽略的數(shù)值色散；在λ/δ增大到10了之后，繼續(xù)減小網(wǎng)格間距不會對κ造成的數(shù)值色散產(chǎn)生明顯的影響，劇烈的色散現(xiàn)象在θ=0°方向依然存在.

因此，CPML中κmax的選擇需要兼顧吸收效率與高頻色散：當(dāng)激勵源主頻較高時趨膚效應(yīng)強(qiáng)烈，為了減小數(shù)值色散可以令κmax取較小值；當(dāng)激勵源主頻較低時CPML的吸收效率也會降低，可以在增大網(wǎng)格間距的同時令κmax取較大值.

圖1　CPML差分近似后vφ/c的各向異性關(guān)系(a) λ/δ=100時，vφ/c與θ、κ的關(guān)系；(b) κ=2和κ=5時，vφ/c與θ、λ/δ的關(guān)系.Fig.1　The anisotropy of vφ/c after difference approximate in CPML(a) While λ/δ=100, the relationship between vφ/c and θ, κ; (b) While κ=2 or κ=5, the relationship between vφ/c and θ, λ/δ.

2.5CPML的最優(yōu)參數(shù)討論

在CPML介質(zhì)中因?yàn)棣实拇嬖?，?shù)值色散也必然存在，依照式(9)的κ分布原則，預(yù)估數(shù)值色散主要發(fā)生在CPML的外層，而最內(nèi)層的vφ/c接近于1.又由于吸收邊界的目的是保證最終反射回內(nèi)部正演區(qū)域的電磁波在一定百分比(γ)以內(nèi)，所以對CPML的數(shù)值色散要求反而不用像內(nèi)部區(qū)域中那樣嚴(yán)格，只需要保證外層CPML中的因數(shù)值色散而產(chǎn)生的高頻震蕩波被內(nèi)層CPML吸收即可，這也解釋了為什么通常式(7)中的αmin取較大值的原因.

為了驗(yàn)證這一參數(shù)設(shè)置原則的可行性，并佐證CPML中數(shù)值色散的存在及影響因素，首先在網(wǎng)格步長為2 m的均勻三維真空空間中做高頻正演實(shí)驗(yàn)：在內(nèi)部區(qū)域中設(shè)置磁偶極子源，并發(fā)射截斷三余弦脈沖，脈沖的最大發(fā)射頻率接近均勻網(wǎng)格數(shù)值色散條件的最大頻率，這樣在限制常規(guī)網(wǎng)格的數(shù)值色散同時，突出坐標(biāo)延伸網(wǎng)格的數(shù)值色散影響.采用的截斷三余弦脈沖源較時諧場源和高斯脈沖源而言，在脈沖的兩個端點(diǎn)處對時間的前五階導(dǎo)數(shù)都為零，具更好的起始和終止時刻的波形平滑性，從而使其脈沖頻譜更加接近理論值，也就能更好地控制因?yàn)樵吹慕尤攵綆У母哳l成分，這一點(diǎn)對于會產(chǎn)生數(shù)值色散的CPML而言非常重要.外部區(qū)域以32層CPML和1層理想電導(dǎo)體(Perfect Electric Conductor，PEC)為吸收邊界，用以展示CPML中電磁波相速度的變化規(guī)律.并按照式(6)提供的關(guān)系，用截斷三余弦脈沖源的脈沖底座寬度(如2×10-7Hz-1)來計算CPML中的αmin和αmax，這樣可以對比不同的內(nèi)層CPML的αmin對外層CPML的數(shù)值色散的抑制作用.

圖2為o((a1)—(a4))、p((b1)—(b4))、q((c1)—(c4))三種CPML的參數(shù)設(shè)置方案在1～4時刻的Hz波場快照，xy平面與發(fā)射線圈在同一平面上，z軸和色標(biāo)均表示Hz的幅值大小.實(shí)驗(yàn)中，o和p參數(shù)設(shè)置主要是為了比較不同αmin對高頻反射波的吸收效果，o和q參數(shù)設(shè)置主要是為了比較不同的κmax產(chǎn)生的數(shù)值色散效果.

圖2(a1)—(c1)所在的時刻，截斷三余弦脈沖已經(jīng)完成發(fā)射，電磁波剛剛進(jìn)入內(nèi)層CPML，此時未發(fā)生明顯的吸收邊界反射，且三種參數(shù)設(shè)置對應(yīng)的波場完全一致.

圖2(a2)—(c2)所在的時刻，電磁波的波陣面進(jìn)入到外層CPML，此時(c2)在吸收邊界內(nèi)的波場與(a2)、(b2)有較大的差異：一方面在電磁波入射角較小的外層CPML處出現(xiàn)了較明顯的數(shù)值色散，而(c2)中數(shù)值色散范圍更廣、幅值更大，說明數(shù)值色散主要發(fā)生在入射角較小的CPML區(qū)域中，且κ越大數(shù)值色散越嚴(yán)重，這與圖1a中的vφ/c隨θ和κ變化的規(guī)律相呼應(yīng)；另一方面(a2)、(b2)的波陣面要比(c2)更靠近PEC，這在圖3b中表現(xiàn)的很明顯，圖3b是離CPML外邊界距離為4 m的網(wǎng)格處的Hz正演結(jié)果，不難發(fā)現(xiàn)電磁波的傳播速度在o、p兩種參數(shù)條件下一致，且比q參數(shù)條件下要快，這說明了電磁波的傳播速度主要受κ的影響，且κ越大vφ越小，這也證明了圖1a中κ與vφ/c的關(guān)系的正確性.

圖2　在與發(fā)射源同一深度平面上，CPML吸收效果的波場快照其他參數(shù)不變的情況下，(a1)—(a4)的參數(shù)設(shè)置為κmax=5、αmin=8.4×10-4、αmax=2.8×10-4； (b1)—(b4)的參數(shù)設(shè)置為κmax=5、αmin=2.8×10-4、αmax=8.4×10-4； (c1)—(c4)的參數(shù)設(shè)置為κmax=11、 αmin=8.4×10-4、αmax=2.8×10-4.其中，(a1)—(c1)的時刻為0.48×10-6s； (a2)—(c2)的時刻為0.74×10-6s； (a3)—(c3)的時刻為1.11×10-6s； (a4)—(c4)的時刻為1.39×10-6s.Fig.2　The snapshot of absorption effect in CPML in the same depth of emission sourceIn the condition that other parameter remain stationary, the parameter settings of (a1)—(a4) are κmax=5, αmin=8.4×10-4, αmax=2.8×10-4; the parameter settings of (b1)—(b4) are κmax=5, αmin=2.8×10-4, αmax=8.4×10-4; the parameter settings of (c1)—(c4) are κmax=11, αmin=8.4×10-4, αmax=2.8×10-4. And the time of (a1)—(c1) is 0.48×10-6s; the time of (a2)—(c2) is 0.74×10-6s; the time of (a3)—(c3) is 1.11×10-6s; the time of (a4)—(c4) is 1.39×10-6s.

圖3　Hz在o、p、q三種參數(shù)設(shè)置中隨時間分布圖對應(yīng)到圖2中的坐標(biāo)，這里(a)為(100, 40)點(diǎn)處的正演結(jié)果，(b)為(100, 4)點(diǎn)處的正演結(jié)果.Fig.3　In three parameter settings as o, p, q, Hz varying over timeCorresponding to the coordinate in Fig.2, (a) is the forward solution in (100,40), and (b) is the forward solution in (100,4).

圖2(a3)—(c3)所在的時刻內(nèi)，由內(nèi)部區(qū)域輻射到CPML中的電磁波已經(jīng)被大部分吸收，但在外層CPML內(nèi)Hz的幅值依舊很大，這種波陣面已經(jīng)離開該FDTD網(wǎng)格而電磁場分量依舊不斷震蕩的現(xiàn)象是典型的數(shù)值色散.若此處的網(wǎng)格處于真空環(huán)境中，則電磁場分量會呈指數(shù)上漲；這里三種參數(shù)設(shè)置中，外層CPML的電導(dǎo)率相同，電磁波的吸收過程會與數(shù)值色散形成某種動態(tài)平衡，則此時Hz的幅值可以作為衡量數(shù)值色散大小的標(biāo)準(zhǔn).較(a3)而言，(b3)中Hz的幅值較大，說明內(nèi)層CPML的αmin取較大值對于數(shù)值色散的抑制作用更強(qiáng)；(c3)中Hz的幅值較大且色散范圍更廣，因數(shù)值色散產(chǎn)生了明顯的反射波，且CPML中κ越大數(shù)值色散產(chǎn)生的區(qū)域越靠內(nèi)層，這都與圖1a中的結(jié)論相一致.

圖2(a4)—(c4)所在的時刻內(nèi)，數(shù)值色散在被CPML壓制的同時也遺留下了明顯的反射波，其中(a4)外層CPML的數(shù)值色散被壓制的最小.結(jié)合圖3a中內(nèi)層CPML處Hz的正演結(jié)果，發(fā)現(xiàn)p參數(shù)條件下高頻色散波的吸收情況較差，而q參數(shù)條件下的高頻色散波出現(xiàn)的較早，這也再次佐證了前面的結(jié)論.

總之，圖2和圖3都說明了CPML網(wǎng)格比常規(guī)均勻網(wǎng)格的數(shù)值色散要求要高，αmin設(shè)置為較大值能在一定程度上提升CPML對數(shù)值色散的抑制效果，但無法從根本上解決因?yàn)棣识傻臄?shù)值色散問題.

雖然根據(jù)圖1b中的結(jié)論，網(wǎng)格間距在δ<λ/10以后繼續(xù)減小不會顯著改善外層CPML的數(shù)值色散，但是CPML中的數(shù)值色散出現(xiàn)的同時意味著電磁波相速度vφ的減小，這會從高頻色散波的吸收效率的角度影響到CPML的吸收性能.結(jié)合各向同性介質(zhì)中針對一維PML的垂直網(wǎng)格間距δ的限制要求(de la Kethulle de Ryhove and Mittet，2014)：

(17)

其中，γ>0為數(shù)值色散帶來的電場或磁場誤差，這里可以取γ=1%；m為電磁波在CPML中傳播的總距離與趨膚深度S(ω)的比值，通常取小于2的值；cσ為電磁波在CPML中的傳播速度；Δt為滿足Courant穩(wěn)定性條件的時間步長.

將式(16)在垂直入射情況下的vφ代入式(17)，我們可以進(jìn)一步將垂直網(wǎng)格間距要求推廣到三維CPML中：

(18)

其中，c0為真空中的光速，ω為電磁波的角頻率.通過式(18)一方面可以限定網(wǎng)格間距δ的大小，另一方面可以在δ確定了的同時，通過γ來約束κmax.

綜合以上論述，選擇CPML的最優(yōu)參數(shù)時，首先應(yīng)該選擇起止時刻波形平滑的激勵源以減少源在接入時刻產(chǎn)生高頻電磁波，然后考慮激勵源的頻譜區(qū)間fmax和fmin，根據(jù)式(6)分別計算αmax和αmin，再考慮激勵源主頻的大小，通過式(18)來確定κmax和δmax，根據(jù)背景場電導(dǎo)率來給出σmax，最后參考正演結(jié)果的數(shù)值色散情況反饋調(diào)整κmax.在本節(jié)實(shí)驗(yàn)γ相同的情況下，按照CPML最優(yōu)參數(shù)設(shè)置方案，能增大43倍的網(wǎng)格間距，減少1/4時間節(jié)點(diǎn)，最終提升正演計算速度兩百多倍.

3井間電磁的數(shù)值模擬及討論

本文的正演模擬主要關(guān)注井間電磁波法的優(yōu)勢探測區(qū)域.根據(jù)法拉第定律，接收線圈處的垂直磁場分量Hz與感應(yīng)電壓成正比，在發(fā)射電流相同的情況下也與電阻抗成正比(Donadille and Al-Ofi，2012).因此，為了突出三維全空間中各網(wǎng)格的靈敏度分布，將單網(wǎng)格大小的高阻體或低阻體逐次嵌入到正演區(qū)域的網(wǎng)格中，并采用2.5節(jié)討論的CPML最優(yōu)參數(shù)分別進(jìn)行正演；令正演結(jié)果中異常體引起的垂直磁場分量為Hzd=Hza-Hzh(Hza為有高阻體或低阻體情況下接收線圈處的垂直磁場分量，Hzh為均勻全空間情況下接收線圈處的垂直磁場分量)，并將Hzd作為評價井間電磁波優(yōu)勢探測區(qū)域的標(biāo)準(zhǔn).

3.1異常體電導(dǎo)率的影響

設(shè)置含CPML的FDTD總計算區(qū)域?yàn)?00 m×200 m×200 m，均勻網(wǎng)格間距為8 m，CPML占據(jù)外部的8層網(wǎng)格.圖5所示為FDTD內(nèi)部區(qū)域，其中發(fā)射探管所在直井的xy坐標(biāo)為(128, 96)，接收探管所在直井的xy坐標(biāo)為(272, 96).

首先討論異常體電導(dǎo)率的變化Δσ對Hzd的影響：設(shè)置發(fā)射線圈坐標(biāo)(128,96,96)，接收線圈坐標(biāo)(272,96,104)，取背景電導(dǎo)率為0.02 S·m-1，并將單個網(wǎng)格大小的異常體設(shè)置到(200,104,104)處，異常體電導(dǎo)率在0.01～0.024 S·m-1之間變化.

如圖4所示為七次FDTD的正演結(jié)果，當(dāng)發(fā)射線圈采用三余弦脈沖時，接收線圈處的Hzd的波形接近正弦波波形.由于僅僅是異常體電導(dǎo)率的變化不會引起Hzd的相位變化，故每一條曲線均會在T1和T2時刻達(dá)到極值，且T1時刻對應(yīng)三余弦脈沖源的極大值所在的時刻，T1時刻Hzd的幅值大于T2時刻.統(tǒng)計正演計算的結(jié)果，每條曲線均有Hzd(T1)∶Hzd(T2)=-1.23，以及參考感應(yīng)測井幾何因子理論(張庚驥，1982)的擬合關(guān)系式Hzd(T1)=-18.099Δσ3+0.931Δσ2-0.052Δσ.這表明了井間電磁的測量靈敏度與異常體電導(dǎo)率的相對變化量正相關(guān)；那么在井間電磁波的三維正演中，可以用Hzd(T1)來代表異常體二次場的振幅，用T1來代表異常體二次場的相位.這樣一方面可以精簡反演中的觀測數(shù)據(jù)，又一方面可以將正演截斷時間設(shè)置到三維空間中的最大T1之后，從而提升2倍以上的正反演速度.

圖4　電導(dǎo)率漸變的異常體在發(fā)射線圈、接收線圈中間時，Hzd隨時間的變化曲線Fig.4　Hz varying over time while anomalous grid′s conductivity changes nearby the middle of transmitter coil and receiver coil

圖5　高阻體的二次場的Hzd(T1)在三維空間中的分布圖井軸上黑色圓環(huán)代表發(fā)射和接收線圈，(a1)、(a2)、(a3)分別為接收線圈的z=88、104、120時的Hzd(T1)，(b1)、(b2)、(b3)分別為接收線圈的z=88、104、120時的T1.Fig.5　Hzd(T1) of high resistance′s secondary field in three-dimensional spaceThe black circle in well axis means transmitter coil and receiver coil, (a1), (a2), (a3) is Hzd(T1) correspond with z=88,104,120, and (b1), (b2), (b3) is T1 correspond with z=88, 104,120.

3.2異常體位置的影響

基于圖4的認(rèn)識，在接下來的異常體位置的變化對Hzd的影響的討論中，我們只需要分析全空間中Hzd(T1)和T1的值即可.令電導(dǎo)率為0.01 S·m-1的高阻網(wǎng)格在背景電導(dǎo)率為0.02 S·m-1的三維內(nèi)部區(qū)域中移動，每移動一次做一次三維正演，將高阻網(wǎng)格所在的位置作為x、y、z，將每個高阻網(wǎng)格對應(yīng)的正演結(jié)果中接收線圈處的Hzd(T1)和T1分別作為圖5a和圖5b中的色標(biāo)，可以得到如圖5所示的三維分布圖.

從圖5(a1)—(a3)中可以看出，Hzd的靈敏度分布整體上為橢球體(幅值相對較大的紅、橙、藍(lán)、紫區(qū)域)，且長軸在收發(fā)線圈的連線上.井間內(nèi)側(cè)區(qū)域的Hzd為正值，外側(cè)區(qū)域的Hzd為負(fù)值；Hzd的靈敏度與該位置到接收、發(fā)射線圈距離的乘積近似呈負(fù)相關(guān).需要注意的是，在x方向上，高阻體在臨近線圈的內(nèi)側(cè)區(qū)域和外側(cè)區(qū)域分別產(chǎn)生了Hzd的極大值和極小值，在y方向和z方向上卻都接近零值，而靠近井軸的區(qū)域恰巧是泥漿侵入嚴(yán)重的地層，地層電阻率的空間變化相對復(fù)雜.因此，在井間電磁的正反演研究中，可以在收發(fā)線圈附近加密網(wǎng)格以實(shí)現(xiàn)井周復(fù)雜地層的精細(xì)描述，也可以結(jié)合陣列感應(yīng)測井或電阻率測井的泥漿侵入信息，和聲電成像測井或地層傾角測井的地層剖面信息，來提供線性化反演的初值及約束.

相比Hzd而言，三維空間中T1的靈敏度分布要簡單很多，圖5(b1)—(b3)中T1最小的位置都在收發(fā)線圈的連線上，且T1的空間等勢面呈橢球體向外逐漸增大；從T1的顏色標(biāo)尺上看，井間內(nèi)側(cè)區(qū)域的T1梯度較小，外側(cè)區(qū)域的T1梯度較大，說明T1能夠更好的提供外側(cè)區(qū)域的二次場走時信息.

由于T1只受異常體位置的影響，靈敏度在井間外側(cè)區(qū)域較大，而Hzd同時受異常體電導(dǎo)率和位置的影響，靈敏度在靠近接收線圈或發(fā)射線圈的區(qū)域較大，故井間電磁反演中T1貢獻(xiàn)了更多異常體位置的信息，Hzd貢獻(xiàn)了更多異常體電導(dǎo)率的信息.為了通過T1來反演異常體的位置，最好的方案是在接收、發(fā)射井的斜交方向上再做一組井間電磁觀測，結(jié)合兩口接收井中的T1信息，理論上就能通過三個橢球體等勢面的交點(diǎn)確定異常體的位置.而在油田生產(chǎn)中，對于在目標(biāo)儲層內(nèi)相距較近三口井不容易找到的情況，可以退而求其次，選用兩口并不平行的井段(如斜井段和直井段)，同樣可以實(shí)現(xiàn)異常體位置的標(biāo)定.

4結(jié)論與建議

針對井間電磁勘探中發(fā)射頻率低、網(wǎng)格間距小、計算區(qū)域小的情況，采用CPML作為FDTD的吸收邊界條件，為了發(fā)揮CPML的低頻吸收優(yōu)勢，本文首先就CPML的最優(yōu)參數(shù)選擇做了討論，然后將其應(yīng)用到井間電磁正演中，得到結(jié)論為：

(1) 在影響CPML吸收效果的三個關(guān)鍵參數(shù)中，α決定CPML的優(yōu)勢吸收頻率，σ和κ影響CPML的吸收效率.但是，κ>1會使CPML中出現(xiàn)不可避免的數(shù)值色散，κ越大、θ越小則數(shù)值色散越嚴(yán)重，CPML的數(shù)值色散主要發(fā)生在外層，如不能很好的壓制，外層電磁場會出現(xiàn)持久的高頻震蕩并影響到內(nèi)部正演區(qū)域.

(2) 異常體的垂直磁場分量Hzd和Hzd的第一極值時刻T1可以作為評價優(yōu)勢探測區(qū)域的標(biāo)準(zhǔn).正演實(shí)驗(yàn)結(jié)果顯示，僅僅是異常體電導(dǎo)率的變化不會引起T1的變化，同一位置處Hzd(T1)比上Hzd(T2)為定值；Hzd的優(yōu)勢探測區(qū)域是以收發(fā)線圈的連線為長軸的橢球體內(nèi)側(cè)區(qū)域，越靠近收發(fā)線圈Hzd的幅值越大，高阻異常體對應(yīng)的內(nèi)側(cè)區(qū)域?yàn)闃O大值，外側(cè)區(qū)域?yàn)闃O小值；T1的優(yōu)勢探測區(qū)域則是以收發(fā)線圈的連線為長軸的橢球體外側(cè)區(qū)域，越靠近收發(fā)線圈連線的中點(diǎn)T1的值越小.

針對以上結(jié)論，提出CPML與井間電磁勘探中的建議如下：

(1) 為了兼顧低頻電磁波的吸收和高頻數(shù)值色散的壓制，CPML的最優(yōu)參數(shù)設(shè)置模式為：首先選擇起止時刻波形平滑的激勵源；然后根據(jù)激勵源的頻譜范圍來計算α，αmin>αmax能有效吸收數(shù)值色散帶來的高頻電磁波；根據(jù)激勵源主頻的大小來確定κmax和δmax，當(dāng)激勵源主頻較高時令κmax取較小值，當(dāng)激勵源主頻較低時令κmax取較大值；再根據(jù)背景場電導(dǎo)率選擇合適的σmax；最后參考正演結(jié)果反饋調(diào)整CPML的參數(shù).

(2) 能有效提高井間電磁正反演信噪比的方法有：加密收發(fā)線圈及異常體附近的網(wǎng)格；增加陣列感應(yīng)測井或電阻率測井以獲得泥漿侵入剖面；增加聲電成像測井或地層傾角測井以獲取地層剖面信息.

(3) 井間電磁勘探中最好的取井方案是采用三口(及以上)不在同一直線上的井做兩組(及以上)觀測，實(shí)際生產(chǎn)中也可以選擇兩口不平行的井段，以減小異常體位置的多解性.

References

Al-Ali Z A, Al-Buali M H, AlRuwaili S, et al. 2009. Looking deep into the reservoir.OilfieldReview, 21(2): 38-47.

Alumbaugh D L, Morrison H F. 1995. Theoretical and practical considerations for crosswell electromagnetic tomography assuming a cylindrical geometry.Geophysics, 60(3): 846-870. Bérenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves.JournalofComputationalPhysics, 114(2): 185-200. Bérenger J P. 1996. Perfectly matched layer for the FDTD solution of wave-structure interaction problems.IEEETransactionsonAntennasandPropagation, 44(1): 110-117. Bérenger J P. 2007. Perfectly Matched Layer (PML) for Computational Electromagnetics. Arizona: Morgan & Claypool. Bérenger J P. 2012. An optimized CFS-PML for wave-structure interaction problems.IEEETransactionsonElectromagneticCompatibility, 54(2): 351-358.Carcione J M, Gei D, Picotti S, et al. 2012. Cross-hole electromagnetic and seismic modeling for CO2detection and monitoring in a saline aquifer.JournalofPetroleumScienceandEngineering, 100: 162-172.

Chew W C, Weedon W H. 1994. A 3D perfectly matched medium from modified Maxwell′s equations with stretched coordinates.MicrowaveandOpticalTechnologyLetters, 7(13): 559-604. de la Kethulle de Ryhove S, Mittet R. 2014. 3D marine magnetotelluric modeling and inversion with the finite-difference time-domain method.Geophysics, 79(6): E269-E286.Donadille J M, Al-Ofi S M. 2012. Crosswell electromagnetic response in a fractured medium.Geophysics, 77(3): D53-D61.

Elsherbeni A, Demir V. 2009. The Finite-Difference Time-Domain Method for Electromagnetics with MATLAB Simulations. Raleigh: SciTech Publishing.

Kuzuoglu M, Mirrta R. 1996. Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEEMicrowaveandGuidedWaveLetters, 6(12): 447-449.

Li Z H, Huang Q H. 2014. Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in

transient electromagnetic method modeling.ChineseJ.Geophys. (in Chinese), 57(4): 1292-1299, doi: 10.6038/cjg20140426.Maclennan K, Karaoulis M, Revil A. 2014. Complex conductivity tomography using low-frequency crosswell electromagnetic data.Geophysics, 79(1): E23-E38. Newman G. 1995. Crosswell electromagnetic inversion using integral and differential equations.Geophysics, 60(3): 899-911.Roden J A, Gedney S D. 2000. Convolution PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media.MicrowaveandOpticalTechnologyLetters, 27(5): 334-339.

Shen J S, Wang Z G, Ma C, et al. 2014. Application of the cross-hole electromagnetic method (CHEM) in hydrocarbon reservoir monitoring.OilGeophysicalProspecting(in Chinese), 49(1): 213-224.

Wei B J, Chen T, Hou X L, et al. 2014. Simulating cross-hole electromagnetic field′s response with metal casing using Green′s functions of radial-layered medium and integral equations.JournalofChinaUniversityofPetroleum(in Chinese), 38(1): 57-63.

Wilt M J, Alumbaugh D L, Morrison H F, et al. 1995. Crosswell electromagnetic tomography: System design considerations and field results.Geophysics, 60(3): 871-885.

Zeng W C, Zhao W J, Zang D F. 2001. Application research of crosshole electromagnetic tomography.ChineseJ.Geophys. (in Chinese), 44(3): 411-420.

Zhang G J. 1982. High order geometrical factors of induction logging.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 25(4): 370-378.

附中文參考文獻(xiàn)

李展輝, 黃清華. 2014. 復(fù)頻率參數(shù)完全匹配層吸收邊界在瞬變電磁法正演中的應(yīng)用. 地球物理學(xué)報, 57(4): 1292-1299, doi: 10.6038/cjg20140426.

沈金松, 王志剛, 馬超等. 2014. 井間電磁油氣儲層監(jiān)測技術(shù)的發(fā)展和應(yīng)用. 石油地球物理勘探, 49(1): 213-224.

魏寶君, 陳濤, 侯學(xué)理等. 2014. 利用徑向成層介質(zhì)的Green函數(shù)和積分方程模擬含金屬套管井間電磁場的響應(yīng). 中國石油大學(xué)學(xué)報(自然科學(xué)版), 38(1): 57-63.

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

張庚驥. 1982. 感應(yīng)測井的高次幾何因子. 地球物理學(xué)報, 25(4): 370-378.

(本文編輯何燕)

Three-dimensional cross-well electromagnetic modeling considering numerical dispersion in convolutional perfectly matched layers

FANG Si-Nan1, PAN He-Ping1*, DU Ting2, WANG Zhi1, DENG Cheng-Xiang1

1InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China2HubeiCoalGeologicalSurveyInstitution,ChinaNationalAdministrationofCoalGeology,Wuhan430070,China

AbstractTaking Convolutional Perfectly Matched Layers (CPML) as the absorbing boundary conditions, this paper applied Finite-Difference Time-Domain (FDTD) into 3D modeling of cross-well electromagnetic (CWEM). The study proved that it was the real stretch which induced numerical dispersion in CPML, given the anisotropy effect of phase velocity of electromagnetic wave from the conventional real stretch and grid interval, and verified numerical dispersion in CPML through wave field snapshot; further on, generalized the restriction between the dominant frequency of source and maximum of real stretch, maximum of grid interval in CPML, so as to improve optimal parameter settings in CPML. Based on optimal CPML, we utilized the vertical magnetic field and the traveling time of second field to describe the distribution of sensitivity, so as to designate the dominant area of cross-well electromagnetic prospecting, and proposed the logging suite for restraint and preferred plan for selecting well which were required in the process of the modeling and inversion of cross-well electromagnetic.

KeywordsCross-Well Electromagnetic (CWEM); 3D modeling; Finite-Difference Time-Domain (FDTD); Convolutional Perfectly Matched Layers (CPML); Numerical dispersion

基金項(xiàng)目國家自然科學(xué)基金項(xiàng)目(41074086)資助．

作者簡介方思南，男，1987年生，博士研究生，主要從事測井處理解釋和井間電磁正反演方面的研究. E-mail: fangsinan@163.com *通訊作者潘和平，男，1953年生，教授，博士生導(dǎo)師，主要從事測井與井中物探的教學(xué)和科研工作.E-mail: panpinge@163.com

doi:10.6038/cjg20160531 中圖分類號P631

收稿日期2015-02-02，2016-04-06收修定稿

方思南， 潘和平， 杜婷等. 2016. 考慮卷積完全匹配層數(shù)值色散的井間電磁三維正演.地球物理學(xué)報，59(5):1888-1897,doi:10.6038/cjg20160531.

Fang S N, Pan H P, Du T, et al. 2016. Three-dimensional cross-well electromagnetic modeling considering numerical dispersion in convolutional perfectly matched layers.ChineseJ.Geophys. (in Chinese),59(5):1888-1897,doi:10.6038/cjg20160531.