• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2016-06-30 07:38:47方思南潘和平杜婷王智鄧呈祥
    地球物理學(xué)報 2016年5期
    關(guān)鍵詞:井間參數(shù)設(shè)置色散

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

    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.

    猜你喜歡
    井間參數(shù)設(shè)置色散
    “光的折射”“光的色散”知識鞏固
    “光的折射”“光的色散”知識鞏固
    “光的折射”“光的色散”知識鞏固
    『光的折射』『光的色散』隨堂練
    煤層氣井間抽機(jī)理及故障處理方法研究及應(yīng)用
    中國煤層氣(2019年4期)2019-11-23 08:42:50
    井間示蹤劑監(jiān)測在復(fù)雜斷塊油藏描述中的應(yīng)用
    錄井工程(2017年1期)2017-07-31 17:44:42
    蟻群算法求解TSP中的參數(shù)設(shè)置
    動車環(huán)境下U900異頻切換參數(shù)設(shè)置探討
    斜井井間地震三維射線追蹤方法
    基于MATLAB仿真的井下變壓器參數(shù)設(shè)置研究
    国产精品一国产av| 久久99精品国语久久久| 精品国产乱码久久久久久男人| 人体艺术视频欧美日本| 久久人妻熟女aⅴ| 大型av网站在线播放| 亚洲一区二区三区欧美精品| 欧美精品一区二区免费开放| 18禁黄网站禁片午夜丰满| 成人国产av品久久久| 午夜免费鲁丝| 丝袜美足系列| 国产片特级美女逼逼视频| 欧美成人精品欧美一级黄| 精品人妻在线不人妻| 人体艺术视频欧美日本| 大码成人一级视频| 高清欧美精品videossex| 黑人巨大精品欧美一区二区蜜桃| 黑丝袜美女国产一区| 在线观看免费午夜福利视频| 日日爽夜夜爽网站| 亚洲精品一区蜜桃| 国产三级黄色录像| 日韩伦理黄色片| 天堂8中文在线网| 国产亚洲av高清不卡| 成年动漫av网址| 精品国产一区二区三区久久久樱花| 女警被强在线播放| 国产99久久九九免费精品| 亚洲中文字幕日韩| 久久鲁丝午夜福利片| 香蕉国产在线看| 成人18禁高潮啪啪吃奶动态图| 久热爱精品视频在线9| 男女无遮挡免费网站观看| 国产伦人伦偷精品视频| 各种免费的搞黄视频| 国产高清国产精品国产三级| 丰满迷人的少妇在线观看| 国产免费一区二区三区四区乱码| 超碰97精品在线观看| a级毛片黄视频| 美女视频免费永久观看网站| 亚洲精品国产av成人精品| 国产1区2区3区精品| 在线观看免费高清a一片| 又黄又粗又硬又大视频| 午夜福利影视在线免费观看| 99国产精品一区二区三区| 亚洲精品第二区| 日韩制服丝袜自拍偷拍| 精品视频人人做人人爽| 国产在线视频一区二区| 精品国产乱码久久久久久男人| 欧美亚洲 丝袜 人妻 在线| 国产极品粉嫩免费观看在线| 久久久久国产一级毛片高清牌| 两个人看的免费小视频| 人妻人人澡人人爽人人| 激情视频va一区二区三区| 日韩电影二区| 国产精品久久久久久人妻精品电影 | 一区二区三区激情视频| 十八禁网站网址无遮挡| 日韩一卡2卡3卡4卡2021年| 大片免费播放器 马上看| 老熟女久久久| 亚洲一码二码三码区别大吗| 久久午夜综合久久蜜桃| 久久久久久人人人人人| 国产精品.久久久| 日本一区二区免费在线视频| 免费日韩欧美在线观看| 宅男免费午夜| 午夜老司机福利片| 国产精品熟女久久久久浪| 搡老岳熟女国产| 亚洲欧美一区二区三区黑人| 黑人猛操日本美女一级片| 久久国产精品影院| 男女床上黄色一级片免费看| bbb黄色大片| 亚洲国产欧美网| av天堂在线播放| 国产av精品麻豆| 在线av久久热| 亚洲欧美日韩另类电影网站| 午夜福利,免费看| 波多野结衣一区麻豆| 国产av国产精品国产| 男人舔女人的私密视频| 久久精品国产亚洲av高清一级| 国产激情久久老熟女| 一级a爱视频在线免费观看| 午夜免费男女啪啪视频观看| 操出白浆在线播放| 2021少妇久久久久久久久久久| 欧美在线黄色| a级毛片在线看网站| 中文字幕最新亚洲高清| 天天躁夜夜躁狠狠躁躁| 亚洲欧美精品自产自拍| 丰满人妻熟妇乱又伦精品不卡| 一级黄片播放器| 欧美老熟妇乱子伦牲交| 在线观看免费午夜福利视频| 老司机亚洲免费影院| 国产主播在线观看一区二区 | 欧美成人精品欧美一级黄| 亚洲av日韩精品久久久久久密 | 亚洲精品在线美女| 天天添夜夜摸| 久久久精品94久久精品| 亚洲专区中文字幕在线| 国产熟女欧美一区二区| av天堂久久9| 丝袜在线中文字幕| 精品免费久久久久久久清纯 | 久久毛片免费看一区二区三区| 免费人妻精品一区二区三区视频| 午夜日韩欧美国产| 亚洲第一av免费看| 亚洲五月婷婷丁香| 一级毛片我不卡| 在线观看一区二区三区激情| 欧美精品一区二区免费开放| 一本久久精品| 成年动漫av网址| 欧美亚洲 丝袜 人妻 在线| xxxhd国产人妻xxx| 在线观看一区二区三区激情| 大香蕉久久成人网| 午夜91福利影院| 欧美在线黄色| 欧美人与性动交α欧美精品济南到| 久久 成人 亚洲| videos熟女内射| 天天躁夜夜躁狠狠久久av| 91老司机精品| 欧美av亚洲av综合av国产av| 欧美成狂野欧美在线观看| 丁香六月欧美| 两个人免费观看高清视频| 日本av免费视频播放| 中文字幕av电影在线播放| av在线app专区| 国产成人精品无人区| av一本久久久久| 亚洲精品国产区一区二| 成人亚洲精品一区在线观看| 一级片'在线观看视频| 99热全是精品| 99精国产麻豆久久婷婷| 中文字幕最新亚洲高清| 国产深夜福利视频在线观看| 99国产精品一区二区蜜桃av | 久久精品国产亚洲av涩爱| 午夜老司机福利片| 777久久人妻少妇嫩草av网站| 亚洲,欧美精品.| 中文字幕人妻熟女乱码| 九色亚洲精品在线播放| 久久精品人人爽人人爽视色| 免费久久久久久久精品成人欧美视频| 成人国产一区最新在线观看 | 黄色怎么调成土黄色| 免费在线观看视频国产中文字幕亚洲 | 色94色欧美一区二区| 色精品久久人妻99蜜桃| 七月丁香在线播放| 亚洲成人手机| 亚洲欧美激情在线| 波多野结衣av一区二区av| 99re6热这里在线精品视频| 日韩一本色道免费dvd| 国产男人的电影天堂91| 中文字幕精品免费在线观看视频| av天堂在线播放| 五月开心婷婷网| 日本a在线网址| 99热国产这里只有精品6| 亚洲精品国产色婷婷电影| 男人爽女人下面视频在线观看| 亚洲 国产 在线| 国产一卡二卡三卡精品| 国产女主播在线喷水免费视频网站| 亚洲av男天堂| 国产1区2区3区精品| 丝袜在线中文字幕| 亚洲第一青青草原| 天天躁夜夜躁狠狠久久av| 国产精品欧美亚洲77777| 亚洲欧美中文字幕日韩二区| 搡老岳熟女国产| 欧美xxⅹ黑人| 啦啦啦在线观看免费高清www| 热99国产精品久久久久久7| 丁香六月天网| 精品国产一区二区三区久久久樱花| 日韩大片免费观看网站| 每晚都被弄得嗷嗷叫到高潮| 黄频高清免费视频| 国产99久久九九免费精品| 美女中出高潮动态图| 视频在线观看一区二区三区| 久久国产精品人妻蜜桃| 久久久久久免费高清国产稀缺| 亚洲国产看品久久| 香蕉丝袜av| 国产97色在线日韩免费| 超碰成人久久| 国产精品麻豆人妻色哟哟久久| 国产精品一区二区在线观看99| 免费观看a级毛片全部| 国产精品 欧美亚洲| 在现免费观看毛片| 国产真人三级小视频在线观看| 精品国产乱码久久久久久男人| 一本—道久久a久久精品蜜桃钙片| 亚洲自偷自拍图片 自拍| 少妇精品久久久久久久| 男女下面插进去视频免费观看| 亚洲精品国产一区二区精华液| 一边摸一边抽搐一进一出视频| av欧美777| 视频在线观看一区二区三区| 国产日韩一区二区三区精品不卡| 激情视频va一区二区三区| 老汉色av国产亚洲站长工具| 多毛熟女@视频| 亚洲伊人色综图| 久久精品成人免费网站| 又黄又粗又硬又大视频| 国产精品一区二区在线观看99| 国产视频首页在线观看| 国产精品麻豆人妻色哟哟久久| 波多野结衣一区麻豆| 亚洲国产精品成人久久小说| av有码第一页| 精品少妇久久久久久888优播| 亚洲国产看品久久| 亚洲欧美成人综合另类久久久| 精品少妇一区二区三区视频日本电影| 免费观看av网站的网址| 精品国产国语对白av| 日韩电影二区| 大香蕉久久成人网| 欧美日韩视频高清一区二区三区二| 亚洲精品第二区| 人妻 亚洲 视频| 国产成人精品无人区| 两人在一起打扑克的视频| 中文字幕av电影在线播放| 纯流量卡能插随身wifi吗| 91麻豆精品激情在线观看国产 | 又紧又爽又黄一区二区| 尾随美女入室| 麻豆国产av国片精品| 一级毛片女人18水好多 | 日本色播在线视频| 国产99久久九九免费精品| 日韩制服丝袜自拍偷拍| 男女之事视频高清在线观看 | 国产老妇伦熟女老妇高清| 另类亚洲欧美激情| 亚洲一卡2卡3卡4卡5卡精品中文| 好男人视频免费观看在线| 亚洲九九香蕉| 精品一品国产午夜福利视频| 欧美黑人欧美精品刺激| 国产在线一区二区三区精| 9色porny在线观看| 亚洲五月婷婷丁香| 国产成人av教育| 免费高清在线观看视频在线观看| 日韩人妻精品一区2区三区| 欧美日韩成人在线一区二区| 国产91精品成人一区二区三区 | 午夜激情av网站| 亚洲精品自拍成人| 国产老妇伦熟女老妇高清| 免费一级毛片在线播放高清视频 | 久久毛片免费看一区二区三区| 在现免费观看毛片| 在线观看国产h片| 一区在线观看完整版| 午夜福利一区二区在线看| 在现免费观看毛片| 中文精品一卡2卡3卡4更新| 精品国产一区二区三区久久久樱花| 国产精品久久久久久精品古装| 久久综合国产亚洲精品| 男人操女人黄网站| 大话2 男鬼变身卡| avwww免费| 国产亚洲av高清不卡| 亚洲国产精品一区三区| 久9热在线精品视频| 欧美日本中文国产一区发布| 中文字幕亚洲精品专区| 丝袜美足系列| 国产精品人妻久久久影院| 国产精品久久久久久人妻精品电影 | av在线老鸭窝| 麻豆乱淫一区二区| 50天的宝宝边吃奶边哭怎么回事| 日本色播在线视频| 视频区欧美日本亚洲| 97在线人人人人妻| 亚洲精品中文字幕在线视频| 日日爽夜夜爽网站| e午夜精品久久久久久久| 天堂中文最新版在线下载| 国产又爽黄色视频| 亚洲精品成人av观看孕妇| 久久热在线av| 国产成人av教育| 国产成人精品久久二区二区91| 尾随美女入室| 亚洲精品国产av蜜桃| 午夜91福利影院| 婷婷丁香在线五月| 只有这里有精品99| 久久女婷五月综合色啪小说| 丁香六月天网| 日韩 亚洲 欧美在线| 涩涩av久久男人的天堂| 超碰成人久久| 男女边吃奶边做爰视频| 大码成人一级视频| 精品国产乱码久久久久久男人| 日本91视频免费播放| 97人妻天天添夜夜摸| 国产av国产精品国产| 国产男人的电影天堂91| 欧美精品一区二区免费开放| 午夜福利免费观看在线| 18禁国产床啪视频网站| 欧美黑人精品巨大| 国产欧美日韩精品亚洲av| 岛国毛片在线播放| 免费在线观看完整版高清| 日韩电影二区| 久久久国产精品麻豆| 老司机影院毛片| 欧美xxⅹ黑人| 中文字幕人妻熟女乱码| 18禁裸乳无遮挡动漫免费视频| 国产一区亚洲一区在线观看| 国产精品免费大片| 亚洲av男天堂| 一本—道久久a久久精品蜜桃钙片| 桃花免费在线播放| 在线观看国产h片| 国产精品久久久久久精品电影小说| 老汉色∧v一级毛片| 欧美97在线视频| 老鸭窝网址在线观看| 日日夜夜操网爽| 精品国产一区二区三区四区第35| 嫁个100分男人电影在线观看 | 国产免费又黄又爽又色| 久久久精品国产亚洲av高清涩受| 极品少妇高潮喷水抽搐| 国产成人免费观看mmmm| av在线app专区| 国产成人一区二区三区免费视频网站 | 亚洲第一青青草原| 亚洲情色 制服丝袜| 欧美老熟妇乱子伦牲交| 久久人人爽av亚洲精品天堂| 高清视频免费观看一区二区| 久久久久久久精品精品| videosex国产| 久久精品国产综合久久久| 成人亚洲精品一区在线观看| 国产免费一区二区三区四区乱码| 香蕉丝袜av| 女人爽到高潮嗷嗷叫在线视频| 叶爱在线成人免费视频播放| www日本在线高清视频| 欧美xxⅹ黑人| 国产黄色视频一区二区在线观看| 亚洲五月婷婷丁香| 久久久久久人人人人人| 亚洲伊人久久精品综合| av在线老鸭窝| 国产精品三级大全| 亚洲国产精品999| av天堂久久9| 亚洲av美国av| 亚洲 欧美一区二区三区| 99国产精品免费福利视频| 我的亚洲天堂| 久久久国产精品麻豆| 精品国产国语对白av| 久久国产精品大桥未久av| 精品国产国语对白av| 精品国产乱码久久久久久男人| 欧美 亚洲 国产 日韩一| 国产伦人伦偷精品视频| 曰老女人黄片| 1024视频免费在线观看| 国产一级毛片在线| 亚洲国产av影院在线观看| 午夜免费成人在线视频| 中文字幕色久视频| 久久久精品免费免费高清| 丰满人妻熟妇乱又伦精品不卡| 1024香蕉在线观看| 一区二区三区四区激情视频| 人体艺术视频欧美日本| 在线 av 中文字幕| 亚洲精品自拍成人| 美女视频免费永久观看网站| 久久精品久久久久久噜噜老黄| 嫩草影视91久久| 久久久精品国产亚洲av高清涩受| 天天影视国产精品| svipshipincom国产片| 久久久欧美国产精品| 成人18禁高潮啪啪吃奶动态图| www.999成人在线观看| 真人做人爱边吃奶动态| 国产精品av久久久久免费| 又黄又粗又硬又大视频| 交换朋友夫妻互换小说| 99热全是精品| 欧美日韩一级在线毛片| av在线app专区| 亚洲综合色网址| 国产午夜精品一二区理论片| 精品一区二区三区av网在线观看 | 又紧又爽又黄一区二区| 免费观看a级毛片全部| 18禁观看日本| 亚洲免费av在线视频| 欧美精品高潮呻吟av久久| √禁漫天堂资源中文www| 国产又色又爽无遮挡免| xxxhd国产人妻xxx| 美女脱内裤让男人舔精品视频| 国产精品久久久久成人av| 母亲3免费完整高清在线观看| 天堂俺去俺来也www色官网| 日本色播在线视频| 国产野战对白在线观看| www.精华液| 人人妻人人爽人人添夜夜欢视频| 日日夜夜操网爽| 国产精品av久久久久免费| 一边摸一边抽搐一进一出视频| 国产成人影院久久av| 欧美 日韩 精品 国产| 99热网站在线观看| 日韩电影二区| 亚洲欧美日韩高清在线视频 | 好男人视频免费观看在线| 又黄又粗又硬又大视频| 黄色视频不卡| 久久久久精品国产欧美久久久 | 久久这里只有精品19| 国产免费又黄又爽又色| 久久久精品国产亚洲av高清涩受| 美女视频免费永久观看网站| 亚洲精品久久久久久婷婷小说| 黄色视频在线播放观看不卡| 亚洲欧美精品自产自拍| 精品国产国语对白av| 色94色欧美一区二区| 精品少妇内射三级| 在线精品无人区一区二区三| 国产片内射在线| 久久毛片免费看一区二区三区| 成人三级做爰电影| 国产在线观看jvid| 久久久精品国产亚洲av高清涩受| 亚洲av成人不卡在线观看播放网 | 青草久久国产| 日韩一区二区三区影片| 老鸭窝网址在线观看| 欧美精品啪啪一区二区三区 | 一级黄片播放器| 一级毛片黄色毛片免费观看视频| 国产1区2区3区精品| 久久精品国产亚洲av涩爱| 亚洲av电影在线观看一区二区三区| 国产真人三级小视频在线观看| 亚洲国产av新网站| 亚洲国产欧美日韩在线播放| 只有这里有精品99| 国产精品99久久99久久久不卡| 亚洲九九香蕉| 少妇猛男粗大的猛烈进出视频| 一本一本久久a久久精品综合妖精| 夫妻午夜视频| 男人舔女人的私密视频| 亚洲激情五月婷婷啪啪| 久久精品亚洲熟妇少妇任你| 久久精品国产亚洲av涩爱| 黑人猛操日本美女一级片| 成人国产一区最新在线观看 | 国产1区2区3区精品| 午夜福利,免费看| 国产日韩一区二区三区精品不卡| 日韩 欧美 亚洲 中文字幕| 狂野欧美激情性xxxx| 精品人妻一区二区三区麻豆| 99精品久久久久人妻精品| 日韩一区二区三区影片| 天天操日日干夜夜撸| 真人做人爱边吃奶动态| 久久天堂一区二区三区四区| 少妇人妻 视频| 国产黄色视频一区二区在线观看| 满18在线观看网站| 亚洲精品久久午夜乱码| 在线观看一区二区三区激情| 日本av免费视频播放| 欧美人与性动交α欧美软件| 欧美日韩综合久久久久久| 天天躁日日躁夜夜躁夜夜| 十八禁网站网址无遮挡| 伦理电影免费视频| 国产欧美日韩一区二区三 | av视频免费观看在线观看| 啦啦啦在线观看免费高清www| www.自偷自拍.com| 日本午夜av视频| 久久久久久人人人人人| 久久99精品国语久久久| 亚洲精品一卡2卡三卡4卡5卡 | 大片免费播放器 马上看| av在线老鸭窝| 日本一区二区免费在线视频| 久久久久久人人人人人| 亚洲av成人精品一二三区| 十八禁高潮呻吟视频| av网站在线播放免费| 国产免费又黄又爽又色| 男女边吃奶边做爰视频| 国产又色又爽无遮挡免| 国产视频首页在线观看| 一区二区日韩欧美中文字幕| 国产激情久久老熟女| 日本91视频免费播放| 桃花免费在线播放| 国产成人系列免费观看| 一边亲一边摸免费视频| 赤兔流量卡办理| av天堂久久9| 深夜精品福利| 亚洲av电影在线观看一区二区三区| 亚洲精品一二三| 看免费av毛片| 十八禁网站网址无遮挡| 国产1区2区3区精品| 99国产综合亚洲精品| 亚洲欧美清纯卡通| 黄频高清免费视频| 99国产精品一区二区三区| 欧美精品一区二区大全| 老司机在亚洲福利影院| 色精品久久人妻99蜜桃| 侵犯人妻中文字幕一二三四区| 超色免费av| 欧美成人精品欧美一级黄| 最黄视频免费看| 亚洲 国产 在线| e午夜精品久久久久久久| 巨乳人妻的诱惑在线观看| 男女边吃奶边做爰视频| 亚洲第一青青草原| av国产精品久久久久影院| 亚洲成色77777| 少妇裸体淫交视频免费看高清 | 久久精品亚洲av国产电影网| 久久精品亚洲熟妇少妇任你| av网站免费在线观看视频| 亚洲欧美一区二区三区国产| 男女边摸边吃奶| 亚洲av综合色区一区| 黄色片一级片一级黄色片| 新久久久久国产一级毛片| 色综合欧美亚洲国产小说| 免费在线观看完整版高清| 啦啦啦视频在线资源免费观看| 久久女婷五月综合色啪小说| 最近手机中文字幕大全| 最新的欧美精品一区二区| 狠狠精品人妻久久久久久综合| 久久久久久亚洲精品国产蜜桃av| 一区二区日韩欧美中文字幕| 两个人看的免费小视频| 美女高潮到喷水免费观看| 国语对白做爰xxxⅹ性视频网站| 老鸭窝网址在线观看| 国产精品久久久久久精品古装| 99香蕉大伊视频| 97在线人人人人妻| xxx大片免费视频| 青春草视频在线免费观看| 免费久久久久久久精品成人欧美视频| 成年动漫av网址| 国产日韩一区二区三区精品不卡| 亚洲中文av在线| 欧美性长视频在线观看| 亚洲av欧美aⅴ国产| 激情视频va一区二区三区| 中文字幕人妻熟女乱码| 新久久久久国产一级毛片| 2018国产大陆天天弄谢| 午夜免费男女啪啪视频观看|