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

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

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

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

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

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

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

1引言

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

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

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

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

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

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

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

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

2.1含CPML的FDTD時域方程

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

(1)

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

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

(2)

(3)

其中，

(4)

(5)

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

2.2α影響CPML的吸收頻率

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

(6)

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

(7)

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

2.3σ影響CPML的吸收效率

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

(8)

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

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

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

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

(9)

其中，κmax為CPML最外層的值，通常可以取5～11(Roden and Gedney，2000).

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

-ω(n+0.5)Δt).

(10)

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

(11)

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

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

(12)

(13)

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

(14)

(15)

(16)

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

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

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

圖1　CPML差分近似后vφ/c的各向異性關系(a) λ/δ=100時，vφ/c與θ、κ的關系；(b) κ=2和κ=5時，vφ/c與θ、λ/δ的關系.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介質中因為κ的存在，數(shù)值色散也必然存在，依照式(9)的κ分布原則，預估數(shù)值色散主要發(fā)生在CPML的外層，而最內層的vφ/c接近于1.又由于吸收邊界的目的是保證最終反射回內部正演區(qū)域的電磁波在一定百分比(γ)以內，所以對CPML的數(shù)值色散要求反而不用像內部區(qū)域中那樣嚴格，只需要保證外層CPML中的因數(shù)值色散而產(chǎn)生的高頻震蕩波被內層CPML吸收即可，這也解釋了為什么通常式(7)中的αmin取較大值的原因.

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

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

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

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

圖2　在與發(fā)射源同一深度平面上，CPML吸收效果的波場快照其他參數(shù)不變的情況下，(a1)—(a4)的參數(shù)設置為κmax=5、αmin=8.4×10-4、αmax=2.8×10-4； (b1)—(b4)的參數(shù)設置為κmax=5、αmin=2.8×10-4、αmax=8.4×10-4； (c1)—(c4)的參數(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ù)設置中隨時間分布圖對應到圖2中的坐標，這里(a)為(100, 40)點處的正演結果，(b)為(100, 4)點處的正演結果.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)所在的時刻內，由內部區(qū)域輻射到CPML中的電磁波已經(jīng)被大部分吸收，但在外層CPML內Hz的幅值依舊很大，這種波陣面已經(jīng)離開該FDTD網(wǎng)格而電磁場分量依舊不斷震蕩的現(xiàn)象是典型的數(shù)值色散.若此處的網(wǎng)格處于真空環(huán)境中，則電磁場分量會呈指數(shù)上漲；這里三種參數(shù)設置中，外層CPML的電導率相同，電磁波的吸收過程會與數(shù)值色散形成某種動態(tài)平衡，則此時Hz的幅值可以作為衡量數(shù)值色散大小的標準.較(a3)而言，(b3)中Hz的幅值較大，說明內層CPML的αmin取較大值對于數(shù)值色散的抑制作用更強；(c3)中Hz的幅值較大且色散范圍更廣，因數(shù)值色散產(chǎn)生了明顯的反射波，且CPML中κ越大數(shù)值色散產(chǎn)生的區(qū)域越靠內層，這都與圖1a中的結論相一致.

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

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

雖然根據(jù)圖1b中的結論，網(wǎng)格間距在δ<λ/10以后繼續(xù)減小不會顯著改善外層CPML的數(shù)值色散，但是CPML中的數(shù)值色散出現(xiàn)的同時意味著電磁波相速度vφ的減小，這會從高頻色散波的吸收效率的角度影響到CPML的吸收性能.結合各向同性介質中針對一維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)，我們可以進一步將垂直網(wǎng)格間距要求推廣到三維CPML中：

(18)

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

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

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

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

3.1異常體電導率的影響

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

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

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

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

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

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

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

4結論與建議

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

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

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

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

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

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

(3) 井間電磁勘探中最好的取井方案是采用三口(及以上)不在同一直線上的井做兩組(及以上)觀測，實際生產(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.

附中文參考文獻

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

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

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

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

張庚驥. 1982. 感應測井的高次幾何因子. 地球物理學報, 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

基金項目國家自然科學基金項目(41074086)資助．

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

doi:10.6038/cjg20160531 中圖分類號P631

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

方思南， 潘和平， 杜婷等. 2016. 考慮卷積完全匹配層數(shù)值色散的井間電磁三維正演.地球物理學報，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.