999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

氣槍震源激發條件對走時變化觀測結果的影響

2018-11-05 10:56:26周青云陳俊磊
地震研究 2018年2期

周青云 陳俊磊

摘要:基于4支單槍容量2000 in3的氣槍組成的陣列,2017年2月在賓川地震信號發射臺大銀甸水庫內開展了不同激發氣槍壓力、沉放深度和水平位移情況下的對比試驗。利用發射臺周邊9個流動寬頻帶數字地震儀與震源附近3個參考臺的數據,使用互相關時延檢測的方法計算了不同激發條件下的氣槍信號的走時變化,并結合激發條件進行了分析。研究表明:(1)不同氣槍壓力或沉放深度下,激發的氣槍信號的波形和互相關系數存在明顯差異:(2)氣槍壓力從9MPa升為15MPa,走時變化可達0.06s:(3)走時變化的可能原因是主頻段(3~5Hz)中的相對高頻成分發生變化,相對低頻成分未變化,二者耦合在一起形成了偽走時變化:(4)由于氣槍信號不同時窗走時變化不盡相同,技術手段難以消除壓力的影響,所以建議未來氣槍激發時固定激發條件;利用已有資料做走時變化相關研究時選用震幅相近的數據。

關鍵詞:氣槍震源;走時變化;激發條件

中圖分類號:P315.3 文獻標識碼:A 文章編號:1000-0666(2018)02-0264-09

0 引言

板塊運動會導致地下介質應力狀態發生改變,地下介質應力狀態的改變與地震的孕育和發生密切相關,研究地殼介質應力隨時間的變化是研究地震規律的有效途徑。地震波是獲取地球內部結構和介質狀態變化等重要而常見的手段,人們對地球內部的知識絕大部分來源于對地震波的研究(陳颙,朱日祥,2005)。隨著數字地震臺數量的不斷增加、資料的逐漸積累和新方法的引入,利用重復地震和背景噪聲等方法研究地球介質變化取得了重要進展。

重復地震是指發生在斷層同一位置上具有近乎相似的震級和震源機制解,并且在同一臺站上波形記錄高度相似的一組天然地震(Nadeau et al,1995; Schaff,Beroza,2004)。利用重復地震,王鵬等(2016)研究了山東乳山震群前后的地殼介質波速變化。隨著地震背景噪聲成像技術的發展(Yao et al,2006:Bensen et al,2010:Zheng et al,2017;鄭定昌,王俊,2017),背景噪聲也被用于地下介質變化的監測研究中。王俊等(2016)利用背景噪聲研究了蘆山M7.0地震前后震源區波速變化時空特征。重復地震和背景噪聲的方法是被動的等待和利用地震信號,前者的空間重復性和空間分布局限,后者時間分辨率較低并且受噪聲源影響較大。

近年來,全國建立了多個氣槍震源基地,利用氣槍震源主動向地下發射重復性很高的地震波(王偉濤等,2017)。氣槍震源具有激發間隔短(最小激發間隔低于10min)、波形之間相關系數高(近場氣槍信號相關系數大于0.99)等特點(陳蒙,2014),有效彌補了被動源時間分辨率不足、空間和波形重復性不夠等不足之處。魏蕓蕓等(2016)利用新疆氣槍信號分析了新疆2次5級左右地震前后波速變化。

利用主動源信號開展走時變化的研究中,檢測到傳播路徑長為5~50km時的走時變化大多是10-3~10-2s數量級(楊微等,2010;劉自鳳等,2015;魏蕓蕓等,2016),大氣壓力變化、潮汐應力等因素的影響比走時變化小1個數量級左右(Wang et al,2008; Niu et al,2008;楊微等,2010)。這些研究考慮的都是傳播路徑上影響因素變化導致的走時變化,沒有分析震源變化對走時變化的影響。本文利用2017年2月賓川地震信號發射臺集中激發期間發射臺周邊流動臺的數據,計算了不同氣槍激發條件下的走時變化并分析了可能的成因,為以后的氣槍激發工作提出了建議。

1 實驗概況

1.1 賓川地震信號發射臺

賓川地震信號發射臺位于云南省賓川縣城以西的大銀甸水庫內,建成于2011年,2012年正式投入使用。現擁有4支容量為2000in3的BoltLonglife型氣槍及相關配套設施,4支氣槍沿水平方向排列成正方形陣列沉放于水面以下10m。氣槍激發時,4支氣槍同時向水庫噴射高壓氣體產生地震波。配合氣槍發射臺建立的流動地震臺觀測臺網主要使用Reftec130B和GuralpCMG-40T地震儀,采樣率為100Hz(王彬等,2015)。在氣槍震源南西架設了3個參考臺,將參考臺記錄的波形當作是氣槍子波,與接收臺站信號反卷積可以去除震源影響。3個參考臺中CKT2更遠一些,CKT4與CKTO位置相同但使用了低靈敏度的擺(圖1)。已有研究表明,水庫氣槍陣列單槍信號傳播距離超過35km,疊加后波形可見距離超過200km;并中激發的信號單槍傳播距離達7km,疊加后波形可見距離超過27km(王彬等,2016)。

1.2 實驗概況

2016年2月14-20日賓川地震信號發射臺進行的氣槍主動源實驗中,氣槍陣列總共激發161次。筆者就氣槍壓力、沉放深度和水平位置進行了不同情況下的測試,其中:氣槍壓力分別為9MPa、12MPa和15MPa(圖2a):沉放深度(氣槍陣列到水面的距離)分別為8m、10m、12m(圖2b):水平位置分別為0m、-7m、7m(圖2c)。實驗共耗時約7d,為保證氣槍信號信噪比,盡量選在夜間開展實驗。

2 資料的選取及處理

2.1 資料的選取

當研究時間范圍內氣槍激發次數足夠多時,可以利用疊加的方法來提高信號的信噪比,增大氣槍信號的接收范圍,提高走時變化計算的精度等(劉自鳳等,2014;魏蕓蕓等,2016),但疊加會降低時間分辨率。考慮到本文研究對象的特征及氣槍激發次數的限制,結合各流動臺記錄到單槍波形的信噪比,最終選取氣槍周邊10km范圍內的8個流動臺站及較遠的1個流動臺(圖1)。

2.2 資料處理

2.2.1 利用參考臺數據提取激發時刻

氣槍控制臺在氣槍激發的過程中會生成激發日志記錄了激發時刻,但是由于槍控設備時鐘校準、機械延遲、時間記錄精度等因素的影響,槍控日志中的激發時刻與真實激發時刻存在較大差異。因此,筆者利用氣槍旁邊參考臺(CKTO)記錄到的氣槍信號初至時刻當作氣槍激發時刻。首先,去除參考臺連續記錄中的非零均值和長周期線性趨勢;再對信號進行主頻段3~5Hz(欒奕等,2016)的Butterworth帶通濾波;最后,選取參考波形窗口,利用互相關檢測方法提取激發時刻及截取實驗波形記錄。

震中距較遠的臺站信噪比較低,個別波形可能無法識別。為提高后續計算結果的可信度,我們利用互相關檢測的方法剔除了一些信噪比較低的波形:將3~5Hz帶通濾波后各臺站的161槍信號線性疊加后與每一槍計算相關系數,剔除掉相關系數小于0.6的波形。

2.2.2 利用互相關時延檢測技術提取走時變化

利用互相關時延檢測技術可以得到亞采樣級別(10-3s)的走時變化(Wang et al,2008;劉自鳳等,2015)。首先對參考臺和流動臺的氣槍信號進行3~5Hz的帶通濾波,降低噪聲干擾;隨后將同一槍的流動臺數據與參考臺數據反卷積,對齊到時并去除震源影響,得到臺站到參考臺的格林函數;將各流動臺的格林函數進行疊加,疊加結果作為模板信號;將單槍格林函數與模板干涉,對干涉結果余弦插值后提取得到走時變化。

3 不同激發條件對波形及重復性的影響

3.1 對波形的影響

筆者在實驗中改變了氣槍壓力、沉放深度和水平位移這3個影響因素。理論上氣槍壓力越大,釋放的能量越大,波形的振幅越大;沉放深度和波形的振幅存在一定的關聯;平行于氣槍水平移動方向的的臺站氣槍到時會有更明顯變化。圖3為53265臺在不同氣槍壓力下互相關對齊后的波形,可以看出氣槍壓力越大,激發地震波的振幅越大(圖3a)。圖3b~d是波形局部放大圖,圖中黑色倒三角為波形最低點,到時差異肉眼可見,壓力越小,到時越早。激發時間上,3槍僅相差約4h,可以認為外部環境因素及地殼介質狀態基本不變,“波速變化”的原因應僅來自于震源處。圖4展示的是不同沉放深度的波形圖,可以看出與不同壓力的波形圖具有相似的結果,沉放深度越大,到時越早。

3.2 對互相關系數的影響

相較于炸藥、夯車等主動震源,氣槍震源優勢之一是氣槍震源的重復性很好,10km以內臺站記錄到的信號可以達到0.9以上,震源附近可達到0.99以上(Wang et al.,2012;陳蒙,2014;王彬等,2016)。賓川地震信號發射臺氣槍陣列激發地震波的主頻段大約為3-5Hz(劉自鳳等,2012;陳佳等,2016h秕奕等,2016),因此我們將53265臺的氣槍進行3~5Hz的Butterworth帶通濾波后,兩兩求互相關系數,結果如圖5所示。從圖5可以看出:不同激發條件相關系數之間形成不同的條帶,壓力和沉放深度對相關系數均有影響;相關系數條帶受壓力影響大,受沉放深度影響較小;對角線上形成一系列相關系數大于0.99的正方形(圖5中A1、A2、A3等),表明同一種激發條件下激發的波形相關系數很高;高壓和低壓下激發的波形相關系數可能低于0.80,最低可達到0.73(圖5中B1、 B2等):相關系數超過0.9的占88.7%。

4 不同激發條件對走時變化計算的影響

4.1 走時變化檢測原理

如果x(t),y(t)是能量有限信號,則它們的互相關函數可以定義為:式中:“*”表示復共扼。顯然,互相關函數是2個信號之間時延:的函數。對于地震儀記錄到的2個波形x(t)和y(t),如果二者波形高度相似,只存在時間上的延遲,則可通過計算2個波形的互相關系數函數來得到時延τ。x(t)和y(t)的時延相關系數定義為:式中:T為地震波持續時間。當相關系數取得最大值時,對應的τ就是2個信號的時延。

本文使用的流動臺站的采樣頻率是100Hz,求得的相關函數的采樣間隔為0.01s,相關函數的峰值并不會恰好位于采樣點上,而計算相關函數時又是一次移動一個采樣點,因而真實的走時變化可能和測得的走時變化存在一定差異。因此,對相關函數進行余弦插值,使新得到的相關函數滿足精度要求。

4.2 數據處理流程

根據氣槍源附近流動臺數據的完整性和單槍信號的信噪比,筆者選擇實驗期間9個流動臺的數據和氣槍震源附近的3個參考臺進行分析(圖1)。數據處理流程如下(陳佳等,2017):(1)去除相關系數較小的氣槍信號記錄,閾值取0.6;(2)數據預處理:3~5Hz帶通濾波,去趨勢,去均值;(3)將流動臺氣槍數據與參考臺數據反卷積得到格林函數,以達到對齊到時和壓縮子波的目的;(4)格林函數線性疊加作為模板,將每一條格林函數與模板進行互相關時延檢測,得到相關函數;(5)對相關函數作余弦插值;(6)提取目標窗口的走時變化。

4.3 走時變化計算結果

利用前述方法處理了9個臺站的數據,得到2017年2月14-17日最大振幅附近走時變化圖。圖6為53265臺的走時變化圖,其它8個臺站的走時變化與53265臺一致。從圖6可以看出:(1)走時變化和氣槍壓力高度呈正相關:(2)CKT2對氣槍壓力的變化更加敏感:(3)難以分辨出沉放深度對走時變化的影響:(4)水平位移對走時變化有影響:(5)氣槍氣壓改變可能導致的走時變化最大可達0.06s(53278臺)。

波在水中的傳播速度約為1.5km/s,當氣槍信號人射角為60°時,水平視速度約為1.3km/s,氣槍水平移動14m,走時變化應約為0.01s,與圖6顯示結果相近。與水平移動相比,沉放深度的改變除影響傳播路徑的幾何響應外,還會影響到震蕩氣泡與水體的耦合(陳惠芳等,2016)。本次實驗設計較不合理,沉放深度改變的同時,其它因素也發生了改變,未能剝離出沉放深度的影響。

5 分析與討論

5.1 全時窗的走時變化

從以上分析可以看出,選取振幅峰值附近的時間窗,所有臺站的走時變化在特定時窗內均呈現與氣槍壓力呈正相關的特征,表明空間域內這種相關性均存在。那么時間域內,氣槍信號的各震相的走時變化是否均存在與氣槍壓力正相關的特征?為此筆者將53265臺0~10.2s的氣槍信號以0.3s為窗長,分割為34個窗口,按照前述方法計算了走時變化。為減少影響因素和增加可信度,我們著重關注沉放深度和水平位移均不變、僅氣槍壓力改變的第70~154槍之間的走時變化。氣槍氣壓劇烈改變前后(第114槍)走時變化的趨勢可以簡單的分為上升、不變和下降3類,因此將時間上相鄰、趨勢上一致的走時變化曲線劃分為一類,結果見圖7。34個窗口共分為10類,其中:走時增加的2類(圖7e,i,窗長共3.3 s):走時減小的3類(圖7b,d,g,窗長共2.1s):走時基本不變的5類(窗長共4.8s)。分析表明,峰值附近震相的走時變化與氣槍壓力正相關,其它時間窗走時變化與氣槍壓力變化并非呈現絕對的正相關,也有可能不相關和負相關。

5.2 干涉模板的選擇

前面的研究是將格林函數疊加后作為干涉模板計算得到的走時變化結果。為研究不同干涉模板對走時變化計算的影響,筆者使用了前一槍的格林函數及不同槍壓下某一槍的格林函數作為干涉模板,計算了53265臺時間窗為3.3~3.6s的走時變化(圖8)。從圖8a可以看出,將前一槍格林函數作為干涉模板,其結果等于將格林函數疊加作為干涉模板的結果的微分,即:

Δt(PreOne)=d[Δt(Stack)](3)

從圖8b可以看出,使用不同激發壓力下的格林函數作為模板,其走時變化的結果除平臺高度外基本相同。以上分析表明,選擇不同的干涉模板計算出來的走時變化基本一致,即干涉模板的選擇不是走時變化的成因。

5.3 反卷積的影響

將流動臺記錄到的氣槍信號與參考臺的信號反卷積得到的格林函數具有對齊到時、壓制噪音的作用,但是會弱化震源信號不重復帶來的影響。筆者采用2.2節中的數據處理流程但不進行反卷積,利用波形記錄提取走時變化,結果見圖9。從圖9中可以看出,反卷積后格林函數提取的走時變化幅度約0.02s,遠小于波形記錄提取的走時變化幅度0.15s,表明反卷積會弱化震源不重復的影響。波形記錄提取的走時變化與激發壓力的相關系數較低,3個周期分別為0.30、0.29和0.23,而格林函數的相關系數為0.92。綜合分析表明,使用格林函數計算走時變化會弱化震源不重復性的影響,但能提高計算結果和激發條件的相關程度。

5.4 成因推測

氣槍產生的地震近場子波可分為壓力脈沖和氣泡脈沖,前者源于高壓氣體噴出瞬間的高頻(>20Hz)脈沖,后者源于氣泡在水中的低頻(<10Hz)震蕩(夏季等,2016),本文使用的是氣泡震蕩產生的低頻段信號。有學者認為氣槍氣壓越大,氣泡震蕩周期越大,主頻越低(何漢漪,2001;陳蒙,2014),也有學者發現氣槍壓力對主脈沖與第1個氣泡脈沖的時間間隔影響顯著,對主頻影響不顯著(夏季等,2016)。本文研究結果與后者一致。

筆者計算了9MPa和15MPa氣槍信號的歸一化傅氏譜及時頻圖,將二者的時頻譜相減,發現:(1)傅氏譜的峰值均為3.45Hz:(2)3~3.7Hz頻段、4.1~4.3頻段二者基本一致,3.7~4.1Hz頻段、4.3~5Hz頻段差異較大(圖10a);(3)二者的歸一化時頻圖4.3Hz以下形態基本一致,4.3Hz以上形態存在差異(圖10b~c);(4)圖10d中存在豎向條紋。這說明氣槍壓力的改變影響的并不是信號的主頻,而是更高的一些頻段,因此我們認為走時變化的成因是:氣槍壓力改變導致氣槍子波主頻段中相對高頻成分發生變化,變化后的高頻成分與未變化的低頻成分耦合在一起,形成了偽走時變化。

5.5 建議

利用50km范圍內臺站接收到的主動源信號監測到的走時變化一般在10-2~10-3s數量級,例如:魏蕓蕓等(2016)觀測到沙灣M5.0地震前后存在0.05~0.06s的走時變化;楊微等(2010)觀測到綿竹M5.6地震前后存在5~9ms的走時變化;劉自鳳等(2015)觀測到賓川主動源無地震期間走時變化波動幅度±0.02s;王寶善等(2016)利用賓川主動源計算的走時日變化數量級為10-3s。這些走時變化與氣槍氣壓導致的走時變化大小相當或更小,因此走時變化的研究中需去除氣槍壓力變化的影響。

通過前面的分析可以看出,氣槍震源的主頻段(3~5Hz)受氣槍壓力影響明顯且無明顯規律,接收臺站的位置、用于反卷積的參考臺的位置、選擇的時間窗口的位置和長度對走時變化計算均有影響,因此難以給出一個可靠的去除氣槍壓力變化影響的方法。建議后續激發實驗中,選用恒定的激發壓力;對于已有的數據,建議使用震幅基本一致的氣槍信號。

6 結論

本文利用2017年2月賓川地震信號發射臺集中激發的數據,研究了不同激發壓力、沉放深度和水平位移條件下氣槍信號的走時變化特征,取得了以下研究結果:

(1)不同壓力和沉放深度下的波形存在肉眼可見的差別;相同激發條件下激發的地震波的主頻段相關系數可達到0.99,不同激發條件下的相關系數最低為0.73。

(2)當激發條件發生變化時,即使波形相關系數大于0.9,也有可能產生一些偽走時變化,這種走時變化反應了震源的變化而非介質或外界因素的變化,反卷積并不能完全去掉震源的影響。氣槍氣壓變化對走時變化影響較大,氣槍信號震幅峰值附近震相的走時變化與壓力變化成正比,對于1~15km范圍內的臺站,最大可達0.06s;無法分辨出沉放深度對走時變化的影響,但根據其對波形的影響可以推測出對走時的影響可能達到10-2s級別;水平位移對走時變化的影響主要體現在傳播幾何路徑改變上。

(3)氣槍信號不同震相的走時變化與壓力變化的相關性存在差異:峰值附近的走時變化與氣壓變化正相關,其它震相有可能正相關,也可能不相關或負相關。據傅氏譜及時頻譜推測,其成因可能是氣槍壓力改變后,信號主頻段中的相對高頻成分發生了變化,變化后的高頻成分與未變化的低頻成分耦合在一起,形成了偽走時變化。

(4)氣槍壓力變化導致的走時變化最大可達0.06s,與已有研究中地震導致的走時變化相當或更大,因此相關研究中應考慮激發壓力的影響。建議以后氣槍的激發壓力恒定;已有數據研究走時變化時,建議使用震幅基本一致的氣槍信號。

兩位審稿老師給出了寶貴意見,云南省地震局主動源創新團隊各成員為開展實驗所付出的艱辛努力,本文使用了中國地震局地球物理研究所王寶善老師提供的走時變化計算程序,在此一并表示感謝。

參考文獻:

陳惠芳,林彬華,金星,等.2016,水庫大容量氣槍震源激發條件優化實驗研究[J].中國地震,32(2):241-248.

陳佳,李孝賓,楊軍,等.2016.云南賓川大容量氣槍震源波形頻譜特征分析[J].中國地震,32(2):216-221.

陳佳,葉泵,高瓊,等.2017.利用氣槍震源信號研究2016年云龍MS5.0地震前后波速變化特征[J].地震研究,40(4):550-556.

陳蒙.2014.利用水庫大容量非調制氣槍陣列進行區域尺度地下結構探測和監測[D].北京:中國地震局地球物理研究所.

陳颙,朱日祥.2005.設立“地下明燈研究計劃”的建議[J].地球科學進展,20(5):485-489.

何漢漪.2001.海上高分辨率地震技術及其應用[M].北京:地質出版社,52-77

劉自鳳,蘇有錦,王寶善,等.2015.賓川主動源地震波走時變化分析方法研究[J].地震研究,38(4):591-597

欒奕,楊宏峰,王寶善,等.2016.大容量氣槍主動源波形資料處理(一):云南賓川[J].中國地震,32(2):305-318.

王寶善,葛洪魁,王彬,等.2016.利用人工重復震源進行地下介質結構及其變化研究的探索和進展[J].中國地震,32(2):168-179.

王彬,李孝賓,劉自鳳,等.2016.賓川地震信號發射臺的震源系統、觀測系統和觀測結果[J].中國地震,32(2):193-201.

王彬,吳國華,蘇有錦,等.2015.賓川地震信號發射臺的選址、建設及初步觀測結果[J].地震研究,38(1):1-6.

王俊,鄭定昌,鄭江蓉,等.2016:利用背景噪聲自相關研究蘆山M7.0地震震源區地殼相對波速的時空變化特征[J].地震地質,38(1):15-168

王鵬,鄭建常,譚毅培.2016.利用重復地震研究山東乳山地區地殼介質波速變化[J].地震學報,38(5):728-738.

王偉濤,王寶善,蔣生森,等.2017.利用氣槍震源探測大陸淺部的地震學研究回顧與展望[J].地震研究,40(4):514-524.

魏蕓蕓,王海濤,蘇金波,等.2016.新疆2次中強地震前氣槍震源反射波震相走時異常變化初步研究[J].中國地震,32(2):270-281.

夏季,金星,蔡輝騰,等.2016.大容量氣槍震源子波時頻特性及其影響因素[J].中國地震,32(2):249-260.

楊微,葛洪魁,王寶善,等.2010.由精密控制人工震源觀測到的綿竹5.6級地震前后波速變化[J].地球物理學報,53(5):1149-1157.

鄭定昌,王俊.2017.基于背景噪聲的川滇地區勒夫波層析成像[J].地震學報,39(5):633-647.

Bensen G D,Ritzwoller M H,Barmin M P,et at.2010.Processing seismicambient noise data to obtain rdiphle broad-band surface wave dis-persion measurements[J].Geophysical Journal of the Royal Astrv-nomibal Society,169(3):1239-1260.

Nadeau R M,Foxall W,McEVIlly T V.1995.clustering and periodic。-currence of microearthquakes on the San Andreas faults at Parkfield,California[J].Science,267(5179):503-507.

Niu F,Silver P G,Daley T M,et al.2008.Preseismic velocity changes ob-served from active source monitoring at the Parkfield SAFOD drillsite[J].Nature,454:204-208,doi:10.1038/nature07111.

Schaff D P,Beroza G C.2004.Coseismic and postseismic velocity changesmeasured by repeating earthquakes[J] .journal of Geophysical Re-search,109,B10302,doi:10.1029/2004JB003011.

Wang Rs,Ge H K,Yang W,et al.2012.Transmitting seismic station mo-nitors fault zone at depth[J].Eos Transactions American Geophysi-cal Union,93(5):49-50.

Wang Bs,Zhu P,Chen Y,et al.2008.Continuous subsurface velocitymeasurement with coda wave interferometry[J].Geophys Res,113,B12313,doi:10.1029/2007JB005023.

Yao H,Vander Hilst RD,De Hoop M V.2006.Surface -wave arraytomography in SE Tibet from ambient seismic noise and two -stationanalysis -I,Phase velocity maps[J].Geophysical Journal Interna-

Zheng D C,Saygin E,Cummins P,et al.2017.Tranadimensional Bayesianseismic ambient noise tomography across SE Tibet[J].Journal of A-sian Earth Sciences,134:86 -93.

主站蜘蛛池模板: 国产成人精品2021欧美日韩 | 成人午夜亚洲影视在线观看| 亚洲国产第一区二区香蕉| 午夜啪啪福利| 亚洲一区二区无码视频| 91香蕉国产亚洲一二三区| 欧美专区日韩专区| 国产av无码日韩av无码网站| 国产综合色在线视频播放线视| 久久99精品久久久大学生| 久久a级片| 国产精品成人一区二区| 91免费在线看| 免费人成又黄又爽的视频网站| 久久香蕉国产线看精品| 国产亚洲欧美在线专区| 五月婷婷导航| 欧美午夜理伦三级在线观看| 亚洲无码电影| 国产第八页| 97在线国产视频| 日韩精品专区免费无码aⅴ| 国产日韩欧美在线视频免费观看| 黄色福利在线| 成人av手机在线观看| 幺女国产一级毛片| 日韩二区三区| AV不卡国产在线观看| 亚洲欧洲日本在线| 蜜臀AV在线播放| 免费高清a毛片| 强乱中文字幕在线播放不卡| 亚洲天堂区| 丰满少妇αⅴ无码区| 亚洲国语自产一区第二页| 综合网久久| 久久精品人人做人人爽电影蜜月| 无码内射中文字幕岛国片 | 欧美视频在线观看第一页| 天堂岛国av无码免费无禁网站| 青青操视频在线| 亚洲精品日产精品乱码不卡| 久久国产乱子| 免费无码在线观看| 黄色在线不卡| 国产日韩欧美一区二区三区在线| 97青青青国产在线播放| 久久一本日韩精品中文字幕屁孩| 国产精品无码翘臀在线看纯欲| 99视频在线免费| 国产精品私拍在线爆乳| 国产日韩欧美精品区性色| 波多野结衣无码中文字幕在线观看一区二区| 国产在线拍偷自揄拍精品| 久热中文字幕在线| 精品丝袜美腿国产一区| 亚洲国产精品一区二区高清无码久久| 看你懂的巨臀中文字幕一区二区| 欧美19综合中文字幕| 日韩中文字幕免费在线观看| 亚洲码一区二区三区| 欧美精品啪啪一区二区三区| 日本免费a视频| 最新国产精品第1页| 国产成人一二三| 青草精品视频| 国产成人精品免费av| 亚洲人成亚洲精品| 亚洲无码视频图片| 国产欧美自拍视频| 四虎国产永久在线观看| 国产精品手机在线播放| 欧美啪啪精品| 一区二区无码在线视频| 情侣午夜国产在线一区无码| 亚洲人成网站色7777| 免费啪啪网址| 精品国产自| 欧美精品另类| 91在线视频福利| 欧美日韩中文国产| 99免费视频观看|