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

    復(fù)頻率參數(shù)完全匹配層吸收邊界在瞬變電磁法正演中的應(yīng)用

    2014-10-09 10:00:14李展輝黃清華
    地球物理學(xué)報(bào) 2014年4期
    關(guān)鍵詞:全空間邊界條件時(shí)域

    李展輝,黃清華

    北京大學(xué)地球與空間科學(xué)學(xué)院地球物理學(xué)系,北京 100871

    1 引言

    時(shí)域有限差分(Finite difference time domain,F(xiàn)DTD)是瞬變電磁法(Transient electromagnetic,TEM)正演常用方法之一.Goldman和Stoyer(1983)首次使用隱式FDTD模擬三維軸對(duì)稱模型,隨后Oristaglio和Hohmann(1984)使用修改的Du Fort-Frankel格式模擬了TEM二維模型,Wang和Hohmann(1993)在Du Fort-Frankel格式的基礎(chǔ)上發(fā)展了一種TEM三維FDTD正演方法,為后續(xù)部分學(xué)者借用(孫懷鳳等,2013;許洋鋮等,2012;Commer and Newman,2004;Endo and Noguchi,2002).FDTD對(duì)TEM的三維正演是在有限模型空間內(nèi)模擬能擴(kuò)散到更大區(qū)域的瞬變場(chǎng),因而需要引入截?cái)噙吔鐥l件.正演中通常采用狄利克雷邊界條件(Dirichlet boundary condition,DBC)(Wang and Hohmann,1993),并設(shè)置邊界足夠遠(yuǎn).DBC在擴(kuò)散場(chǎng)未抵達(dá)或剛抵達(dá)邊界時(shí)效果較好,但隨著延遲時(shí)間的推移,反射場(chǎng)將對(duì)原擴(kuò)散場(chǎng)產(chǎn)生較大的干擾.要消除這種干擾,則需要更大的模型,這將增加內(nèi)存和計(jì)算量.因而一些學(xué)者開始研究其他的吸收邊界來代替DBC,比如Mur吸收邊界(秦臻和胡文寶,2004;姬金祖,2008;岳建華,2007)以及在其基礎(chǔ)上修改而來的新吸收邊界條件(楊海燕和岳建華,2009).但是他們并沒有對(duì)吸收邊界的效果作詳細(xì)的討論,也沒有給出與理論解的對(duì)比.李墩柱(2010)將電磁波正演中應(yīng)用最普遍的完全匹配層(Perfectly Matched Layer,PML)吸收邊界(Berenger,1994,1996;李展輝等,2009;Chew and Weedon,1994)應(yīng)用到TEM半空間模型正演中,然而文中模型的水平范圍為[-3000m,3000m],半空間電導(dǎo)率為0.01S·m-1,但是最長延遲時(shí)間僅為8ms,渦流電場(chǎng)形成的“煙圈”還遠(yuǎn)未到達(dá)邊界,晚期吸收效果并不明確,在參數(shù)設(shè)置上也僅僅說明需要多次嘗試.本文深入研究了PML在TEM正演中的吸收效果,并最終采用了其中的一種——復(fù)頻率參數(shù)完全匹配層(Complex Frequency Shifted Perfectly Matched Layer,CFS-PML)(Kuzuoglu and Mittra,1996;Roden and Gedney,2000;姜永金等,2004;Berenger,2012).首先通過推導(dǎo)擴(kuò)散場(chǎng)在常規(guī)PML內(nèi)部的平面波解的形式,指出了常規(guī)PML在TEM正演中失效的原因,然后給出了擴(kuò)散場(chǎng)在CFS-PML內(nèi)部的平面波解,研究了CFS-PML各項(xiàng)參數(shù)對(duì)吸收效果的影響,并給出了CFS-PML的參數(shù)選擇方案.最后將CFSPML應(yīng)用到全空間和半空間模型中以驗(yàn)證CFS-PML的有效性.同時(shí)為了提高半空間模型的正演效率,本文采用了李展輝和黃清華(2011)提出的新方法,而全空間模型正演則采用普遍使用的Du Fort-Frankel格式.

    2 理論推導(dǎo)

    2.1 含PML的Maxwell頻域方程組

    頻率域帶PML的旋度方程可以表示為(Chew and Weedon,1994)

    其中

    常規(guī)PML中(Chew and Weedon,1994),

    其中κu為網(wǎng)格延拓因子,由于本文采用非等距網(wǎng)格,κu可直接包含于網(wǎng)格中.因此在下文中所有的κu值都設(shè)置為1.令ε′=εΔ-jσ/ω,可得

    2.2 常規(guī)PML和CFS-PML的吸收性質(zhì)

    設(shè)式(5)的平面波解為

    定義復(fù)波速

    可得到頻散方程

    對(duì)于常規(guī)PML,式(6)可寫為

    其中

    而TEM的低頻近似滿足

    那么復(fù)波速可以近似為

    將式(13)代入到式(10)中,并且考慮在x方向吸收,設(shè)置 (σx,σy,σz)= (σx,0,0),可得

    設(shè)

    為PML的吸收系數(shù),APML越小,代表吸收越強(qiáng)烈.從式(15)中可以看出,隨著角頻率ω的減小,APML也會(huì)隨著減小,而當(dāng)ω→0時(shí),亦有APML→0.這意味著PML對(duì)低頻場(chǎng)的吸收隨著頻率的降低而更為強(qiáng)烈,過強(qiáng)的吸收在數(shù)值計(jì)算中會(huì)引起大的異常(Berenger,2000,1996).

    對(duì)于CFS-PML,經(jīng)過類似的推導(dǎo)之后可以得到

    從式(16)可以看出,αx的出現(xiàn)使得CFS-PML避免了在ω接近0時(shí)出現(xiàn)吸收非常強(qiáng)烈的現(xiàn)象.

    設(shè)置CFS-PML的一個(gè)參考頻率:

    把式(17)代入式(16),得

    當(dāng)f?fαx時(shí),式(18)近似為式(14),CFS-PML退化成常規(guī)PML;當(dāng)f?fαx時(shí),式(18)近似為

    這僅僅是x方向的網(wǎng)格延拓.

    瞬變場(chǎng)的優(yōu)勢(shì)頻率會(huì)隨著時(shí)間逐漸變低.當(dāng)優(yōu)勢(shì)頻率足夠低并進(jìn)入PML內(nèi)部時(shí),PML會(huì)對(duì)這些優(yōu)勢(shì)頻率產(chǎn)生過強(qiáng)的吸收進(jìn)而使整個(gè)瞬變場(chǎng)產(chǎn)生異常.而CFS-PML能避免這種情況,但是σx/αx這個(gè)比例非常重要,從式(19)可以看出,對(duì)于f?fαx頻率部分只是一個(gè)網(wǎng)格延拓因子,如果過大同樣會(huì)引起數(shù)值異常.為此我們可以設(shè)置另一個(gè)參考頻率:

    2.3 含CFS-PML的Maxwell時(shí)域方程組

    時(shí)域方程組可以用卷積的形式表達(dá)出來(Roden and Gedney,2000),以Ex為例:

    其中ˉsy(t)和ˉsz(t)分別是1/sy(ω)和1/sz(ω)的拉普拉斯逆變換:

    其中δ(t)是狄拉克函數(shù),ζu(t)可以表述如下:

    其中u(t)為單位階躍函數(shù).那么最終式可以表示為

    其中

    其他分量可以類似地表達(dá)出來.

    2.4 CFS-PML參數(shù)選取研究

    其中L為PML的厚度.假設(shè)場(chǎng)沿x方向傳播,X=x.定義CFS-PML吸收系數(shù)為

    式中

    式(28)中存在兩個(gè)可控參數(shù):fαx和L.吸收系數(shù)A越小代表額外吸收越強(qiáng)烈.很明顯,吸收系數(shù)A是CFS-PML厚度L的遞減函數(shù),L越大,吸收越強(qiáng)烈.下面將對(duì)參數(shù)fαx選取作具體的分析.

    假設(shè)介質(zhì)是均勻的,且電導(dǎo)率為σ=1×10-3S·m-1,場(chǎng)源到PML內(nèi)邊界的距離為L1=1400m,PML厚度為L2=600m,關(guān)注的最大時(shí)間為tmax=1ms,那么抵達(dá)邊界并穿過PML返回場(chǎng)內(nèi)的最低頻率為

    圖1 不同參考頻率下CFS-PML對(duì)擴(kuò)散場(chǎng)的額外吸收系數(shù)比較Fig.1 The additional absorbing coefficient of CFS-PML under different reference frequencies

    3 數(shù)值模擬結(jié)果及討論

    本文除巷道帶異常體模型外,所有正演均采用相同的網(wǎng)格:網(wǎng)格數(shù)量為128×128×128,中央網(wǎng)格大小為Δxmin=Δymin=Δzmin=10m,從中央到四周逐漸增加網(wǎng)格步長,最大處為Δxmax=10Δxmin,Δymax=10Δymin,Δzmax=7Δzmin.圖2描述了網(wǎng)格的主要特征.CFS-PML的層數(shù)為8層,已經(jīng)包含在上述網(wǎng)格中.

    圖2 瞬變電磁三維正演網(wǎng)格示意圖Fig.2 A sketch map of the mesh for TEM 3Dforward modeling

    3.1 均勻全空間

    為了方便與理論解對(duì)比,我們采用了均勻全空間模型.同時(shí)為了探討CFS-PML對(duì)不同介質(zhì)的適應(yīng)性,我們選取了兩種不同的介質(zhì):σ=1×10-2S·m-1和σ=1×10-3S·m-1.以水平正方形線框作為發(fā)射線框,線框邊長L=70m,垂直方向位于z=0處,水平方向位于x-y平面的正中央.采用脈沖形式的發(fā)射電流,最終比對(duì)垂直磁場(chǎng)Hz.

    選擇的測(cè)線為z=0和y=0平面的交線.我們定義在CFS-PML與正常區(qū)域交接處Hz易號(hào)的時(shí)間點(diǎn)為早期、晚期的分界點(diǎn).首先對(duì)σ=1×10-2S·m-1的三個(gè)早期時(shí)間點(diǎn)進(jìn)行對(duì)比,根據(jù)牛之鏈(2007)中式(2-33)計(jì)算得到早晚期分界延遲時(shí)間約為4.5ms,那么首先選取t=0.35ms,1.23ms和2.64ms三個(gè)時(shí)間點(diǎn).圖3為在上述三個(gè)時(shí)間點(diǎn)測(cè)線上分別使用CFS-PML、DBC得到的Hz值與理論值的對(duì)比.可以看出,在早期三者幾乎一致.這是因?yàn)閿U(kuò)散場(chǎng)尚未或者剛剛抵達(dá)邊界,邊界條件的區(qū)別未能體現(xiàn)出來.因此在下文中,本文將著重?cái)U(kuò)散晚期的對(duì)比.

    圖4為t=6.80ms,12.50ms,和20.00ms時(shí)分別使用CFS-PML,DBC得到的Hz值和理論解的對(duì)比.可以看出t=6.8ms時(shí)使用DBC的結(jié)果在邊界處已經(jīng)出現(xiàn)了較大的異常,隨后這種異常開始擴(kuò)散到場(chǎng)的中央,結(jié)果越來越偏離理論解.而使用CFS-PML吸收邊界時(shí),不管在中央還是在邊界,Hz值依然和理論解符合得很好.CFS-PML的吸收作用在晚期非常顯著.

    圖5為σ=1×10-3S·m-1情況下不同時(shí)刻使用CFS-PML、DBC得到的Hz值和理論解之間的對(duì)比,同時(shí)加入了CFS-PML中設(shè)置fαx=∞所得的結(jié)果.σ=1×10-3S·m-1模型中場(chǎng)的擴(kuò)散速度要明顯快于σ=1×10-2S·m-1的情況,但是兩種情況下的場(chǎng)在一定的條件下具有相似性(王華軍,2008).根據(jù)這種相似性,t=2ms時(shí)場(chǎng)的形態(tài)類似于σ=1×10-2S·m-1時(shí)t=20ms情況下場(chǎng)的形態(tài),因此本文將σ=1×10-3S·m-1情況下最晚的一個(gè)時(shí)間點(diǎn)設(shè)定為t=2ms.從圖5可以看出,在晚期CFS-PML的優(yōu)勢(shì)依然非常明顯,而且經(jīng)過適當(dāng)設(shè)置后的CFS-PML也明顯優(yōu)于設(shè)置fαx=∞的CFS-PML.

    圖3 全空間σ=1×10-2S·m-1模型中,t=0.35ms,1.23ms和2.64ms時(shí)使用CFS-PML吸收邊界、DBC的Hz分量和理論解的對(duì)比Fig.3 The comparison among the Hzcomponents obtained from CFS-PML,DBC,and the analytical solutions at the times of 0.35ms,1.23ms,and 2.64ms of a whole space model withσ=1×10-2S·m-1

    圖4 全空間σ=1×10-2S·m-1模型中,t=6.80ms,12.50ms和20.00ms時(shí)使用CFS-PML吸收邊界條件、DBC的Hz分量和理論解的對(duì)比Fig.4 The comparison among the Hzcomponents obtained from CFS-PML,DBC,and the analytical solutions at the times of 6.80ms,12.50ms,and 20.00msof a whole space model withσ=1×1 0-2 S·m-1

    圖5 全空間σ=1×10-3S·m-1模型中,t=0.680ms和2.00ms時(shí)使用CFS-PML吸收邊界條件、DBC的Hz分量和理論解的對(duì)比,對(duì)比中還加入了在CFS-PML中設(shè)置fαx= ∞ 所得的Hz分量Fig.5 The comparison among the Hzcomponents obtained from CFS-PML,DBC,and the analytical solutions at the times of 0.68ms and 2.00ms of a whole space model withσ=1×10-3S·m-1.The Hz component obtained from CFS-PML with fαx= ∞is also included for reference

    為顯示CFS-PML內(nèi)部場(chǎng)的情況,我們將σ=1×10-3S·m-1模型中2ms時(shí)CFS-PML內(nèi)外的場(chǎng)都描繪出來并與理論解進(jìn)行對(duì)比,如圖6所示.可以看出,CFS-PML內(nèi)部的場(chǎng)與理論解截然不同.CFS-PML內(nèi)部出現(xiàn)了Hz的正負(fù)變化,意味著渦流的存在,而理論解并沒有Hz符號(hào)的變化.CFS-PML雖然不能將擴(kuò)散場(chǎng)完全吸收,但是其在吸收的同時(shí)以某種方式將渦流限制在了CFS-PML內(nèi)部,使得正常區(qū)域的擴(kuò)散場(chǎng)不受影響.

    3.2 均勻半空間

    本文采用了李展輝和黃清華(2011)介紹的方法在半空間模型中引入空氣層.由于空氣層中波速遠(yuǎn)大于地下場(chǎng)的擴(kuò)散速度,由式(10)可知,針對(duì)地下介質(zhì)設(shè)計(jì)的CFS-PML在空氣層中對(duì)波的吸收基本為零,因此在空氣層中需要在CFS-PML介質(zhì)內(nèi)設(shè)置一定的電導(dǎo)率進(jìn)行吸收.半空間模型中,地下介質(zhì)的電導(dǎo)率設(shè)置為σ=1×10-2S·m-1,空氣層中電導(dǎo)率設(shè)置為1×10-6S·m-1,并在靠近邊界處逐漸增加.地空邊界設(shè)置在z=0處,激勵(lì)源的位置與形式均與全空間模型一致.

    圖7顯示了使用CFS-PML吸收邊界、DBC所得到的Hz分量和理論解之間的對(duì)比.從圖中可以看出,t=6.80ms時(shí)使用DBC的Hz分量在兩端已經(jīng)偏離了理論解.隨著時(shí)間的推移,誤差越來越大,在20ms時(shí),誤差已經(jīng)在一個(gè)量級(jí)以上.而使用CFS-PML的Hz與理論解的吻合度大幅優(yōu)于使用DBC的情況,但是也有一定的誤差,而且這種誤差也隨時(shí)間的推移緩慢地增加.這可能是CFS-PML不能對(duì)空氣層中的場(chǎng)產(chǎn)生有效的吸收造成的,越到晚期越明顯.圖8詳細(xì)展示了線圈中心Hz相對(duì)誤差隨時(shí)間變化.從圖中可以看出,線圈中心相對(duì)誤差隨著時(shí)間逐漸增加,在20ms時(shí)誤差達(dá)到了13%.實(shí)際應(yīng)用中需要根據(jù)精度要求調(diào)整模型的大小以控制相對(duì)誤差.

    圖6 全空間σ=1×10-3S·m-1模型CFS-PML內(nèi)外Hz分量展示,并與理論解比較以體現(xiàn)CFS-PML內(nèi)部場(chǎng)特征Fig.6 The Hzcomponent in and out the CFS-PML of a whole space model withσ=1×10-3S·m-1 at 2ms

    3.3 巷道模型

    全空間瞬變電磁法典型的應(yīng)用之一為煤礦巷道內(nèi)地下水的探測(cè)(李宇等,2012;于景 等,2011;郭純等,2006).模型如圖9所示.地下水所在的位置電導(dǎo)率較高,因此模型中設(shè)定圍巖電導(dǎo)率為0.05S·m-1,異常體電導(dǎo)率為0.5S·m-1.巷道內(nèi)通常很狹小,因此發(fā)射線圈設(shè)定為3m×3m,并采用重疊回線的方式記錄感應(yīng)電壓.由于該模型不存在理論解,為了獲取類似理論解的參考解,我們將原網(wǎng)格x,y方向各增加60個(gè)網(wǎng)格,在z方向增加40個(gè)網(wǎng)格.總網(wǎng)格數(shù)量達(dá)到188×188×168,以保證在我們關(guān)心的時(shí)間窗內(nèi),邊界反射場(chǎng)沒有進(jìn)入觀測(cè)區(qū)域.考慮到模型尺寸問題,最大的延遲時(shí)間設(shè)置為0.6ms,因?yàn)楦鶕?jù)擴(kuò)散深度公式(牛之璉,2007)

    圖7 半空間模型中,t=6.80ms,12.50ms和20.00ms時(shí)使用CFS-PML吸收邊界條件、DBC的Hz分量和理論解的對(duì)比Fig.7 The comparisons among the Hzcomponents obtained from CFS-PML,DBC,and the analytical solutions at the times of t=6.80ms,12.50ms,and20.00ms of a half space model

    圖8 回線中心Hz的相對(duì)誤差隨時(shí)間分布圖Fig.8 The relative errors of Hzvarying over time at the center of the rectangular loop

    計(jì)算可得此時(shí)的擴(kuò)散深度為124m,已經(jīng)超出了本模型范圍.式(31)中ρ為電阻率,t為延遲時(shí)間.結(jié)果如圖10所示.從圖中可以看出CFS-PML的結(jié)果與參考解一直符合得很好,最終的相對(duì)誤差不超過6%.

    圖9 巷道內(nèi)含異常體模型,模型關(guān)于y=0對(duì)稱Fig.9 Underground roadway model with an anomalous body,which is symmetric about the y=0plane

    圖10 巷道模型的重疊回線測(cè)量垂直感應(yīng)電動(dòng)勢(shì)結(jié)果圖,左下角的小圖為相對(duì)于參考解的相對(duì)誤差Fig.10 Vertical emf of the underground roadway model.The relative errors to the reference solution are shown in the embedded figure

    4 結(jié)論

    本文將CFS-PML應(yīng)用到TEM三維FDTD正演中,通過適當(dāng)?shù)卦O(shè)置CFS-PML參數(shù)保證了數(shù)值計(jì)算的穩(wěn)定性和良好的吸收效果.全空間和半空間數(shù)值結(jié)果表明CFS-PML在瞬變場(chǎng)早期與傳統(tǒng)的狄利克雷邊界條件的結(jié)果一致,晚期則明顯優(yōu)于狄利克雷邊界條件的結(jié)果,并且能有效地節(jié)省計(jì)算時(shí)間和存儲(chǔ)空間.但是針對(duì)地下擴(kuò)散場(chǎng)設(shè)計(jì)的CFSPML并不能有效地吸收空氣中的電磁場(chǎng),導(dǎo)致半空間模型結(jié)果到后期也產(chǎn)生一定的誤差,并且隨著時(shí)間緩慢的增長.實(shí)際應(yīng)用中需要根據(jù)TEM正演的最大延遲時(shí)間以及對(duì)誤差的容忍度適當(dāng)調(diào)整模型的大小,以確保誤差在要求之內(nèi).本文最后將CFSPML應(yīng)用到與實(shí)際情況貼近的巷道模型中,采用重疊回線測(cè)量感應(yīng)電動(dòng)勢(shì),并與參考解進(jìn)行對(duì)比和誤差分析,再一次證明了CFS-PML的有效性.

    Bérenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics,114(2):185-200.

    Bérenger J P.1996.Three-dimensional perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics,127(2):363-379.

    Bérenger J P.1996.Perfectly matched layer for the FDTD solution of wave-structure interaction problems.IEEE Transactions on Antennas and Propagation,44(1):110-117.

    Bérenger J P.2000.Numerical reflection of evanescent waves by PMLs:origin and interpretation in the FDTD case.Expected consequences to other finite methods.Int.J.Numer.Model,13(2-3):103-114.

    Bérenger J P.2012.An optimized CFS-PML for wave-structure interaction problems.IEEE Transactions on Electromagnetic Compatibility,54(2):351-358.

    Chew W C,Weedon W H.1994.A 3Dperfectly matched medium from modified Maxwell′s equations with stretched coordinates.Microwave and Optical Technology Letters,7(13):599-604.

    Commer M,Newman G.2004.A parallel finite-difference approach for 3Dtransient electromagnetic modeling with galvanic sources.Geophysics,69(5):1192-1202.

    Endo M,Noguchi K.2002.Three-dimensional modeling considering the topography for the case of the time-domain electromagnetic method.Amsterdam:Elsevier Science BV,35:85-107.

    Goldman M M,Stoyer C H.1983.Finite-difference calculations of the transient field of an axially-symmetric earth for vertical magnetic dipole excitation.Geophysics,48(7):953-963.

    Guo C,Liu B Z,Bai D H.2006.Prediction of water disasters ahead of tunneling in coal mine using continuous detection by UWTEM.Seismology and Geology(in Chinese),28(3):456-462.

    Ji J Z,Gao X,Liu Z H.2008.Simplified format application of 2D two-rank Mur absorbing boundary condition.Journal of Shenyang Institute of Aeronautical Engineering(in Chinese),25(2):10-13.

    Jiang Y J,Mao J J,Cai S L.2004.Implementation and analysis of the perfectly matched layer media with CFS for the PSTD method.Journal of Microwaves(in Chinese),20(4):36-39.

    Kuzuoglu M,Mittra R.1996.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave and Guided Wave Letters,6(12):447-449.

    Li D Z.2010.Numerical Simulation of Transient Electromagnetic Fields(in Chinese)[Master′s thesis].Beijing:Peking University.

    Li Y,Tong B,Guo C F,etal.2012.Study on application and denoising of mine full space transient electromagnetic method.Coal Engineering(in Chinese),(11):29-31.

    Li Z H,Huang Q H,Wang Y B.2009.A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media.Chinese J.Geophys.(in Chinese),52(7):1915-1922.

    Li Z H,Huang Q H.2011.A 3Dcurvilinear staggered-grid finite difference time-domain method with an impulse source for transient electromagnetic method.Chinese Geophysics(in Chinese),240.

    Niu Z L.2007.Principle of Time Domain Electromagnetic Method(in Chinese).Changsha:Central South University Press.

    Oristaglio M L,Hohmann G W.1984.Diffusion of electromagnetic fields into a two-dimensional earth-a finite-difference approach.Geophysics,49(7):870-894.

    Qin Z,Hu W B.2004.Modified Mur absorbing boundary conditions for lossy medium FDTD calculation.Journal of Jianghan Petroleum Institute(in Chinese),26(2):76-77.

    Roden J A,Gedney S D.2000.Convolution PML(CPML):An efficient FDTD implementation of the CFS-PML for arbitrary media.Microwave and Optical Technology Letters,27(5):334-339.

    Sun H F,Li X,Li S C,etal.2013.Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time.Chinese J.Geophys.(in Chinese),56(3):1049-1064.

    Wang H J.2008.Time domain transient electromagnetism all time apparent resistivity translation algorithm.Chinese J.Geophys.(in Chinese),51(6):1936-1942.

    Wang T,Hohmann G W.1993.A finite-difference,time-domain solution for 3-dimensional electromagnetic modeling.Geophysics,58(6):797-809.

    Xu Y C,Lin J,Li S Y,etal.2012.Calculation of full-waveform airborne electromagnetic response with three-dimension finitedifference solution in time-domain.Chinese J.Geophys.(in Chinese),55(6):2105-2114.

    Yang H Y,Yue J H.2009.Application of absorbing boundary condition in whole-space computation of transient electromagnetic response.Journal of China University of Mining &Technology(in Chinese),38(2):263-268.

    Yu J C,Liu Z Q,Liao J J,etal.2011.Application of full space transient electromagnetic method to mine water prevention and control.Coal Science and Technology(in Chinese),39(9):110-113.

    Yue J H,Yang H Y,Hu B.2007.3Dfinite difference time domain numerical simulation for TEM in-mine.Progress in Geophysics(in Chinese),22(6):1904-1909.

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

    郭純,劉白宙,白登海.2006.地下全空間瞬變電磁技術(shù)在煤礦巷道掘進(jìn)頭的連續(xù)跟蹤超前探測(cè).地震地質(zhì),28(3):456-462.

    姬金祖,高旭,劉戰(zhàn)合.2008.二維Mur二階吸收邊界條件簡化格式的應(yīng)用.沈陽航空工業(yè)學(xué)院學(xué)報(bào),25(2):10-13.

    姜永金,毛鈞杰,柴舜連.2004.CFS-PML邊界條件在PSTD算法中的實(shí)現(xiàn)與性能分析.微波學(xué)報(bào),20(4):36-39.

    李墩柱.2010.地形起伏下的瞬變電磁場(chǎng)數(shù)值模擬[碩士論文].北京:北京大學(xué).

    李宇,童兵,郭昌放等.2012.礦井全空間瞬變電磁法應(yīng)用及其消噪研究.煤炭工程,(11):29-31.

    李展輝,黃清華,王彥賓.2009.三維錯(cuò)格時(shí)域偽譜法在頻散介質(zhì)井中雷達(dá)模擬中的應(yīng)用.地球物理學(xué)報(bào),52(7):1915-1922.

    李展輝,黃清華.2011.瞬變電磁正演之帶激勵(lì)源的三維曲線交錯(cuò)網(wǎng)格時(shí)域有限差分方法.中國地球物理,240.

    牛之璉.2007.時(shí)間域電磁法原理.長沙:中南大學(xué)出版社.

    秦臻,胡文寶.2004.用于有耗介質(zhì)FDTD計(jì)算的改進(jìn)的Mur吸收邊界條件.江漢石油學(xué)院學(xué)報(bào),26(2):76-77.

    孫懷鳳,李貅,李術(shù)才等.2013.考慮關(guān)斷時(shí)間的回線源激發(fā)TEM三維時(shí)域有限差分正演.地球物理學(xué)報(bào),56(3):1049-1064.

    王華軍.2008.時(shí)間域瞬變電磁法全區(qū)視電阻率的平移算法.地球物理學(xué)報(bào),51(6):1936-1942.

    許洋鋮,林君,李肅義等.2012.全波形時(shí)間域航空電磁響應(yīng)三維有限差分?jǐn)?shù)值計(jì)算.地球物理學(xué)報(bào),56(6):2105-2114.

    楊海燕,岳建華.2009.吸收邊界條件在全空間瞬變電磁計(jì)算中的應(yīng)用.中國礦業(yè)大學(xué)學(xué)報(bào),38(2):263-268.

    于景邨,劉振慶,廖俊杰等.2011.全空間瞬變電磁法在煤礦防治水中的應(yīng)用.煤炭科學(xué)技術(shù),39(9):110-113.

    岳建華,楊海燕,胡搏.2007.礦井瞬變電磁法三維時(shí)域有限差分?jǐn)?shù)值模擬.地球物理學(xué)進(jìn)展,22(6):1904-1909.

    猜你喜歡
    全空間邊界條件時(shí)域
    一類帶有Stieltjes積分邊界條件的分?jǐn)?shù)階微分方程邊值問題正解
    帶有積分邊界條件的奇異攝動(dòng)邊值問題的漸近解
    基于時(shí)域信號(hào)的三電平逆變器復(fù)合故障診斷
    基于極大似然準(zhǔn)則與滾動(dòng)時(shí)域估計(jì)的自適應(yīng)UKF算法
    全空間瞬變電磁場(chǎng)低阻層屏蔽效應(yīng)數(shù)值模擬研究
    中國煤炭(2016年1期)2016-05-17 06:11:33
    多層圓線圈在電磁計(jì)算中全空間磁感應(yīng)強(qiáng)度B的分布
    基于時(shí)域逆濾波的寬帶脈沖聲生成技術(shù)
    基于時(shí)域波形特征的輸電線雷擊識(shí)別
    全空間上一類半線性雙調(diào)和方程正解的衰減
    帶Robin邊界條件的2維隨機(jī)Ginzburg-Landau方程的吸引子
    蜜桃亚洲精品一区二区三区| 国产三级在线视频| 99riav亚洲国产免费| 偷拍熟女少妇极品色| 黄片wwwwww| 免费在线观看影片大全网站| 日日摸夜夜添夜夜添小说| 久久久国产成人精品二区| 精品久久久久久久末码| 香蕉av资源在线| 又紧又爽又黄一区二区| 欧美丝袜亚洲另类 | 日本-黄色视频高清免费观看| 91麻豆av在线| 91在线精品国自产拍蜜月| av中文乱码字幕在线| 国产精品一区www在线观看 | 久久久久免费精品人妻一区二区| 日韩大尺度精品在线看网址| 最新中文字幕久久久久| 成人无遮挡网站| 成人国产麻豆网| 黄色欧美视频在线观看| 国产久久久一区二区三区| 免费人成视频x8x8入口观看| 色在线成人网| 色视频www国产| 国产麻豆成人av免费视频| a级毛片a级免费在线| 自拍偷自拍亚洲精品老妇| 国产精品亚洲一级av第二区| 人人妻人人看人人澡| 桃红色精品国产亚洲av| 亚洲欧美日韩卡通动漫| 婷婷精品国产亚洲av在线| .国产精品久久| 免费看光身美女| 看片在线看免费视频| 久久精品国产鲁丝片午夜精品 | 精品国内亚洲2022精品成人| 国产精品嫩草影院av在线观看 | 十八禁网站免费在线| 婷婷精品国产亚洲av| 欧美激情国产日韩精品一区| 中文字幕久久专区| 韩国av在线不卡| 亚洲精华国产精华精| 亚洲国产欧美人成| 国产免费av片在线观看野外av| 欧美最新免费一区二区三区| 国产大屁股一区二区在线视频| а√天堂www在线а√下载| 国产高清不卡午夜福利| 联通29元200g的流量卡| 日韩欧美精品免费久久| 国产一区二区亚洲精品在线观看| 亚洲精品亚洲一区二区| 亚洲国产高清在线一区二区三| 成人毛片a级毛片在线播放| 精品一区二区三区视频在线观看免费| 亚洲av二区三区四区| 欧美+日韩+精品| 一区二区三区免费毛片| 哪里可以看免费的av片| 亚洲乱码一区二区免费版| 亚洲成人免费电影在线观看| 国产av一区在线观看免费| 国产亚洲欧美98| 狂野欧美白嫩少妇大欣赏| 丝袜美腿在线中文| 偷拍熟女少妇极品色| a级毛片免费高清观看在线播放| 亚洲精品日韩av片在线观看| 成年女人看的毛片在线观看| 三级男女做爰猛烈吃奶摸视频| av在线亚洲专区| 在线观看舔阴道视频| 91久久精品国产一区二区成人| 精品人妻1区二区| 少妇熟女aⅴ在线视频| 91狼人影院| 精品人妻熟女av久视频| 久久婷婷人人爽人人干人人爱| 成年女人看的毛片在线观看| 国产午夜福利久久久久久| 欧美最新免费一区二区三区| 精品福利观看| 18+在线观看网站| 国产精品久久视频播放| 日韩中文字幕欧美一区二区| 亚洲欧美精品综合久久99| 女人十人毛片免费观看3o分钟| 少妇丰满av| av专区在线播放| 日韩欧美国产在线观看| 热99re8久久精品国产| 亚洲最大成人手机在线| 国产成人影院久久av| 99久久无色码亚洲精品果冻| 搡老妇女老女人老熟妇| 国产精品三级大全| 日韩欧美一区二区三区在线观看| 国产不卡一卡二| 欧美中文日本在线观看视频| 国内精品久久久久久久电影| 国产三级中文精品| 亚洲色图av天堂| 亚洲午夜理论影院| 国内精品久久久久精免费| 成人鲁丝片一二三区免费| 欧美成人性av电影在线观看| 国产精品久久久久久精品电影| 男人舔奶头视频| 免费观看在线日韩| 久久久久久久久久成人| 国产综合懂色| 久久久久久久久久久丰满 | 国产精品人妻久久久影院| 国产欧美日韩精品亚洲av| 国产91精品成人一区二区三区| 精品国产三级普通话版| av国产免费在线观看| 色在线成人网| 国产乱人视频| 欧美不卡视频在线免费观看| 婷婷丁香在线五月| 91在线精品国自产拍蜜月| 他把我摸到了高潮在线观看| 国产一级毛片七仙女欲春2| 国产精品一区二区三区四区久久| 免费无遮挡裸体视频| 精品免费久久久久久久清纯| 三级国产精品欧美在线观看| 亚洲avbb在线观看| 校园人妻丝袜中文字幕| 色哟哟哟哟哟哟| 精品久久久噜噜| 久久午夜亚洲精品久久| 国产精品久久久久久精品电影| 亚洲精品在线观看二区| 老司机深夜福利视频在线观看| 国产伦在线观看视频一区| 国产精品电影一区二区三区| 哪里可以看免费的av片| 男女边吃奶边做爰视频| 在线免费观看的www视频| 日本-黄色视频高清免费观看| av在线观看视频网站免费| 男女做爰动态图高潮gif福利片| 精品午夜福利在线看| 国产蜜桃级精品一区二区三区| 欧美一区二区国产精品久久精品| 久久久久久久久久黄片| 久久亚洲精品不卡| 国内少妇人妻偷人精品xxx网站| 亚洲欧美日韩东京热| 亚洲精华国产精华液的使用体验 | 亚洲三级黄色毛片| 99久久成人亚洲精品观看| 欧美高清性xxxxhd video| 在现免费观看毛片| 嫁个100分男人电影在线观看| 国模一区二区三区四区视频| 亚洲欧美清纯卡通| 毛片女人毛片| 桃红色精品国产亚洲av| 露出奶头的视频| 极品教师在线免费播放| 床上黄色一级片| 黄色女人牲交| 亚洲精品一区av在线观看| 日本一本二区三区精品| 午夜免费男女啪啪视频观看 | 国产精品国产三级国产av玫瑰| 久久国产乱子免费精品| 99热这里只有精品一区| 俄罗斯特黄特色一大片| 日本色播在线视频| 在线观看av片永久免费下载| 人人妻人人看人人澡| aaaaa片日本免费| netflix在线观看网站| 尤物成人国产欧美一区二区三区| 蜜桃久久精品国产亚洲av| 99久久成人亚洲精品观看| 国内精品宾馆在线| 一边摸一边抽搐一进一小说| 欧美bdsm另类| 免费av毛片视频| 国产午夜精品久久久久久一区二区三区 | 久99久视频精品免费| 99热6这里只有精品| 美女高潮的动态| 成人亚洲精品av一区二区| 欧美高清性xxxxhd video| 欧美最黄视频在线播放免费| 天堂av国产一区二区熟女人妻| 日日摸夜夜添夜夜添小说| 中文字幕av在线有码专区| 天堂影院成人在线观看| 欧美成人一区二区免费高清观看| 免费观看人在逋| 久久精品夜夜夜夜夜久久蜜豆| 精品午夜福利在线看| 国产真实伦视频高清在线观看 | 国产亚洲91精品色在线| 偷拍熟女少妇极品色| 久久久久国内视频| 免费搜索国产男女视频| 欧美zozozo另类| 国产精品日韩av在线免费观看| 久久精品国产亚洲av香蕉五月| 午夜影院日韩av| 日本-黄色视频高清免费观看| av在线亚洲专区| 日日摸夜夜添夜夜添小说| 非洲黑人性xxxx精品又粗又长| 亚洲狠狠婷婷综合久久图片| 日本黄色视频三级网站网址| 亚洲人成伊人成综合网2020| 欧美中文日本在线观看视频| 免费看a级黄色片| 国产三级中文精品| 18禁黄网站禁片免费观看直播| 男女做爰动态图高潮gif福利片| 国产一区二区三区av在线 | 极品教师在线视频| 免费在线观看成人毛片| 又爽又黄无遮挡网站| 尤物成人国产欧美一区二区三区| 免费在线观看日本一区| 日韩在线高清观看一区二区三区 | 国产精品98久久久久久宅男小说| 欧美不卡视频在线免费观看| 色哟哟哟哟哟哟| 午夜福利18| 最新中文字幕久久久久| 国产爱豆传媒在线观看| 国产一区二区三区av在线 | 久久国产乱子免费精品| 22中文网久久字幕| 欧美激情在线99| 一边摸一边抽搐一进一小说| 高清日韩中文字幕在线| 亚洲一级一片aⅴ在线观看| 国产黄色小视频在线观看| 亚洲国产精品sss在线观看| 日韩欧美三级三区| 久久久久久久久大av| 变态另类丝袜制服| 亚洲精品色激情综合| 久久人人爽人人爽人人片va| 日韩欧美免费精品| 免费观看的影片在线观看| 给我免费播放毛片高清在线观看| 欧美+亚洲+日韩+国产| 一进一出好大好爽视频| 俺也久久电影网| 十八禁国产超污无遮挡网站| 国产视频一区二区在线看| 欧美黑人欧美精品刺激| 美女cb高潮喷水在线观看| 亚洲人成网站高清观看| 999久久久精品免费观看国产| 亚洲一区二区三区色噜噜| АⅤ资源中文在线天堂| 男人的好看免费观看在线视频| 免费观看的影片在线观看| 国内精品久久久久精免费| 熟妇人妻久久中文字幕3abv| 中文字幕久久专区| 亚洲在线观看片| 真人做人爱边吃奶动态| 看黄色毛片网站| 18+在线观看网站| 丝袜美腿在线中文| 成年人黄色毛片网站| 51国产日韩欧美| 国产精品女同一区二区软件 | 国产视频一区二区在线看| 久久久久久九九精品二区国产| 久久欧美精品欧美久久欧美| 国产爱豆传媒在线观看| 日本与韩国留学比较| 午夜福利在线在线| 国产男人的电影天堂91| 午夜a级毛片| 噜噜噜噜噜久久久久久91| 免费观看人在逋| 久久久久久大精品| 亚洲va日本ⅴa欧美va伊人久久| 日本在线视频免费播放| 成人鲁丝片一二三区免费| av黄色大香蕉| 亚洲一级一片aⅴ在线观看| 又黄又爽又免费观看的视频| 中文字幕av成人在线电影| 看黄色毛片网站| 韩国av在线不卡| 精品99又大又爽又粗少妇毛片 | 一级av片app| 黄色日韩在线| 免费观看在线日韩| 一本精品99久久精品77| 18禁在线播放成人免费| 色精品久久人妻99蜜桃| eeuss影院久久| 日本五十路高清| 亚洲美女黄片视频| 毛片一级片免费看久久久久 | 欧美国产日韩亚洲一区| 欧美人与善性xxx| 亚洲av.av天堂| 给我免费播放毛片高清在线观看| 亚洲成av人片在线播放无| 欧美bdsm另类| 97人妻精品一区二区三区麻豆| 国产高清有码在线观看视频| av在线老鸭窝| 久久国内精品自在自线图片| 丰满的人妻完整版| 精品免费久久久久久久清纯| 久99久视频精品免费| 直男gayav资源| 国产欧美日韩精品一区二区| 变态另类成人亚洲欧美熟女| 少妇熟女aⅴ在线视频| 欧美区成人在线视频| 国产v大片淫在线免费观看| 免费av不卡在线播放| 日日夜夜操网爽| 国产伦在线观看视频一区| 两个人的视频大全免费| 草草在线视频免费看| 少妇猛男粗大的猛烈进出视频 | www.www免费av| 日日夜夜操网爽| 超碰av人人做人人爽久久| 岛国在线免费视频观看| 久久人人爽人人爽人人片va| 国产爱豆传媒在线观看| 日韩,欧美,国产一区二区三区 | 亚洲五月天丁香| 一个人观看的视频www高清免费观看| 国产亚洲av嫩草精品影院| 在线播放国产精品三级| av在线亚洲专区| 婷婷精品国产亚洲av| 黄片wwwwww| 午夜福利视频1000在线观看| 亚洲人成网站在线播| 免费在线观看日本一区| 国产精品亚洲一级av第二区| 悠悠久久av| 草草在线视频免费看| 亚洲七黄色美女视频| 一进一出抽搐动态| 国产 一区 欧美 日韩| 男插女下体视频免费在线播放| 精品一区二区三区av网在线观看| 国产精品永久免费网站| 淫秽高清视频在线观看| 国产精品一区二区三区四区久久| 日本 欧美在线| 很黄的视频免费| 1000部很黄的大片| 久久久久久久久久黄片| 国产乱人伦免费视频| 国产高清视频在线播放一区| 免费av毛片视频| 99国产精品一区二区蜜桃av| bbb黄色大片| 日韩欧美在线二视频| 久久久国产成人精品二区| 亚洲欧美日韩高清专用| 日本a在线网址| 国产精品美女特级片免费视频播放器| 亚洲第一电影网av| 亚洲色图av天堂| 午夜福利高清视频| 亚洲av美国av| 亚洲精品粉嫩美女一区| 免费搜索国产男女视频| 国产精品久久视频播放| 在现免费观看毛片| 国产精品1区2区在线观看.| 久久久久久九九精品二区国产| 免费高清视频大片| 美女高潮喷水抽搐中文字幕| 香蕉av资源在线| 麻豆成人午夜福利视频| 中文字幕久久专区| 女同久久另类99精品国产91| av国产免费在线观看| 69av精品久久久久久| 97超级碰碰碰精品色视频在线观看| 可以在线观看的亚洲视频| 亚洲色图av天堂| 观看免费一级毛片| 亚洲av电影不卡..在线观看| 久久久久久国产a免费观看| 成熟少妇高潮喷水视频| netflix在线观看网站| 国产精品福利在线免费观看| 亚洲熟妇熟女久久| 精品欧美国产一区二区三| 老女人水多毛片| 又爽又黄a免费视频| 内射极品少妇av片p| 国内精品美女久久久久久| 亚洲成人中文字幕在线播放| 欧美性猛交黑人性爽| 亚洲黑人精品在线| 午夜福利成人在线免费观看| 亚洲人成网站高清观看| 99热这里只有是精品50| 身体一侧抽搐| 国产精品久久久久久精品电影| 精品久久国产蜜桃| 日本 欧美在线| 国产大屁股一区二区在线视频| 成年版毛片免费区| 老司机深夜福利视频在线观看| 日日啪夜夜撸| 99久久精品国产国产毛片| 国产又黄又爽又无遮挡在线| 日本-黄色视频高清免费观看| 国产av不卡久久| 毛片女人毛片| 午夜激情福利司机影院| 少妇的逼水好多| av在线观看视频网站免费| 亚洲熟妇熟女久久| 99在线视频只有这里精品首页| 亚洲色图av天堂| 香蕉av资源在线| 嫩草影院入口| 成人亚洲精品av一区二区| 亚洲在线自拍视频| 日本 av在线| 亚洲黑人精品在线| 久久久久久久久久成人| 人妻丰满熟妇av一区二区三区| 美女大奶头视频| 精品无人区乱码1区二区| 亚洲中文字幕日韩| 色播亚洲综合网| 内地一区二区视频在线| 国产精品无大码| 亚洲成人久久爱视频| 精品久久久久久成人av| 日本五十路高清| 色综合站精品国产| 天天一区二区日本电影三级| 91在线精品国自产拍蜜月| 午夜精品在线福利| 特大巨黑吊av在线直播| 51国产日韩欧美| 亚洲中文日韩欧美视频| 亚洲最大成人中文| 午夜福利18| 深夜精品福利| 亚洲18禁久久av| or卡值多少钱| 97超级碰碰碰精品色视频在线观看| 男人狂女人下面高潮的视频| 国产亚洲欧美98| 蜜桃亚洲精品一区二区三区| 欧美色欧美亚洲另类二区| 如何舔出高潮| 最后的刺客免费高清国语| 久久国内精品自在自线图片| videossex国产| 国产精品女同一区二区软件 | 亚洲成人免费电影在线观看| 中出人妻视频一区二区| 精品午夜福利在线看| 无人区码免费观看不卡| 日韩欧美国产一区二区入口| 观看免费一级毛片| 日本 av在线| 中国美白少妇内射xxxbb| 22中文网久久字幕| 一本精品99久久精品77| 国产精品亚洲美女久久久| 欧美潮喷喷水| 日本免费一区二区三区高清不卡| 国产精品福利在线免费观看| 最近中文字幕高清免费大全6 | 日日干狠狠操夜夜爽| 可以在线观看的亚洲视频| 男女那种视频在线观看| 久久久久久久亚洲中文字幕| 少妇熟女aⅴ在线视频| 久久亚洲精品不卡| 日韩欧美在线乱码| 亚洲真实伦在线观看| 色哟哟哟哟哟哟| 精品午夜福利在线看| 欧美xxxx黑人xx丫x性爽| a在线观看视频网站| 亚洲三级黄色毛片| 露出奶头的视频| 国产午夜精品久久久久久一区二区三区 | 91麻豆av在线| 一进一出抽搐gif免费好疼| 免费av毛片视频| 少妇猛男粗大的猛烈进出视频 | 校园人妻丝袜中文字幕| 日本a在线网址| 国产久久久一区二区三区| 又紧又爽又黄一区二区| 亚洲 国产 在线| 亚洲av熟女| 97碰自拍视频| 毛片一级片免费看久久久久 | 99riav亚洲国产免费| 日韩大尺度精品在线看网址| 免费黄网站久久成人精品| 国产精品亚洲一级av第二区| 午夜a级毛片| 久久6这里有精品| 99热6这里只有精品| 校园人妻丝袜中文字幕| 毛片女人毛片| 亚洲精品乱码久久久v下载方式| 熟妇人妻久久中文字幕3abv| avwww免费| 露出奶头的视频| 精品欧美国产一区二区三| 国产精品一区www在线观看 | 别揉我奶头 嗯啊视频| 伦理电影大哥的女人| 精品乱码久久久久久99久播| 国产人妻一区二区三区在| 日韩欧美在线乱码| 成人国产综合亚洲| 日本与韩国留学比较| 午夜福利视频1000在线观看| 久久久久久久久久久丰满 | 亚洲国产精品久久男人天堂| 嫩草影院新地址| 精品人妻一区二区三区麻豆 | 欧美+亚洲+日韩+国产| 非洲黑人性xxxx精品又粗又长| 日本一本二区三区精品| 乱码一卡2卡4卡精品| 国产又黄又爽又无遮挡在线| 精品乱码久久久久久99久播| 国产欧美日韩精品亚洲av| eeuss影院久久| 亚洲18禁久久av| 国产精品久久久久久av不卡| 春色校园在线视频观看| 99热6这里只有精品| 一级黄色大片毛片| 熟女电影av网| 国产精品亚洲美女久久久| 91av网一区二区| 亚洲18禁久久av| 日韩精品有码人妻一区| 嫩草影院新地址| 内地一区二区视频在线| 国产黄片美女视频| 日韩精品青青久久久久久| av国产免费在线观看| 亚洲精华国产精华精| 欧美日本视频| 国产熟女欧美一区二区| 男女边吃奶边做爰视频| 成人综合一区亚洲| 九色国产91popny在线| 久久久久精品国产欧美久久久| 国产又黄又爽又无遮挡在线| 最近最新中文字幕大全电影3| 欧美日韩国产亚洲二区| 亚洲一区二区三区色噜噜| 淫秽高清视频在线观看| 很黄的视频免费| 真人做人爱边吃奶动态| 日韩精品有码人妻一区| 色在线成人网| 久久精品影院6| 女生性感内裤真人,穿戴方法视频| 有码 亚洲区| 色精品久久人妻99蜜桃| 精品久久久久久久久久久久久| 日本五十路高清| 久久久久国内视频| 麻豆精品久久久久久蜜桃| 日本黄色片子视频| 观看免费一级毛片| 久久久色成人| 婷婷丁香在线五月| 国产精品一区二区三区四区免费观看 | 国产精品久久久久久久电影| 亚洲精品日韩av片在线观看| 99国产精品一区二区蜜桃av| 成人一区二区视频在线观看| 午夜免费成人在线视频| 精华霜和精华液先用哪个| 国产一级毛片七仙女欲春2| 国产精品久久久久久久久免| 黄色视频,在线免费观看| 国产黄片美女视频| 国产 一区精品| 亚洲精品久久国产高清桃花| 91狼人影院| 久久久久久久久大av| 精品久久久久久久久久久久久| 全区人妻精品视频| 亚州av有码| 亚洲中文字幕一区二区三区有码在线看| 国产精品自产拍在线观看55亚洲| 日日摸夜夜添夜夜添av毛片 |